| OptFEM2DP1 Toolbox V1.2b2 Matlab/Octave Optimized P1-Lagrange Finite Element Method in 2D

- Main Page
- Related Pages
- Files

- File List
- File Members

All Files Functions Pages
SquareMesh.m
Go to the documentation of this file.
1 function Th=SquareMesh(N)
2 % function Th=SquareMesh(N)
3 % Initialization of the Mesh structure for a square domain
4 %
5 % Square domain is `[0,1]\times[0,1]`.<br/>
6 %
7 % This mesh has 4 boundary labels :
8 % - label 1 : boundary `y=0`
9 % - label 2 : boundary `x=1`
10 % - label 3 : boundary `y=1`
11 % - label 4 : boundary `x=0` <br/>
12 % There are N+1 points on each boundary.
13 %
14 % Parameters:
15 % N: integer, number of elements on a boundary
16 %
17 % Return values:
18 % Th: mesh structure
19 %
20 % Generated fields of mesh structure:
21 % q: Array of vertices coordinates, `2\times\nq` array. <br/>
22 % `{\q}(\il,j)` is the
23 % `\il`-th coordinate of the `j`-th vertex, `\il\in\{1,2\}` and
24 % `j\in\ENS{1}{\nq}`
25 % me: Connectivity array, `3\times\nme` array. <br/>
26 % `\me(\jl,k)` is the storage index of the
27 % `\jl`-th vertex of the `k`-th triangle in the array `\q` of vertices coordinates, `\jl\in\{1,2,3\}` and
28 % `k\in{\ENS{1}{\nme}}`.
29 % ql: Array of vertices labels, `1\times\nq` array.
30 % mel: Array of elements labels, `1\times\nme` array.
31 % be: Connectivity array for boundary edges, `2\times\nbe` array.<br/>
32 % `\be(\il,l)` is the storage index of the
33 % `\il`-th vertex of the `l`-th edge in the array `\q` of vertices coordinates, `\il\in\{1,2\}` and
34 % `l\in{\ENS{1}{\nbe}}`.
35 % bel: Array of boundary edges labels, `1\times\nbe` array.
36 % nq: total number of vertices, also denoted by `\nq`
37 % nme: total number of elements, also denoted by `\nme`
38 % nbe: total number of boundary edges, also denoted by `\nbe`
39 % areas: Array of areas, `1\times\nme` array. areas(k) is the area of the `k`-th triangle.
40 % lbe: Array of edges lengths, `1\times\nbe` array. `\lbe(j)` is the length of the `j`-th edge.
41 %
42 % See also
43 % #ComputeAreaOpt, #EdgeLengthOpt, #GetMeshOpt
44 % Copyright:
45 % See \ref license
46 h=1/N;t=0:h:1;
47 if isOctave()
48 mesh=msh2m_structured_mesh(t,t,1,1:4); % package msh
49 me=mesh.t(1:3,:);
50 q=mesh.p;
51 else
52 [x,y] = meshgrid(t,t);
53 q=[x(:) y(:)]';
54 tri = delaunay(x,y);
55 me=tri';
56 end
57
58 nq=length(q);
59 nme=length(me);
60
61 ql=zeros(1,nq);
62 be=zeros(2,4*N);
63 bel=zeros(1,4*N);
64 % Points du bord
65 ql=zeros(1,nq);
66 I=find(q(2,:)==0);
67 ql(I)=1;
68 be(:,1:N)=[I(1:end-1);I(2:end)];
69 bel(:,1:N)=1;
70 I=find(q(1,:)==1);
71 ql(I)=2;
72 be(:,N+1:2*N)=[I(1:end-1);I(2:end)];
73 bel(:,N+1:2*N)=2;
74 I=find(q(2,:)==1);
75 ql(I)=3;
76 I=fliplr(I);
77 be(:,2*N+1:3*N)=[I(1:end-1);I(2:end)];
78 bel(:,2*N+1:3*N)=3;
79 I=find(q(1,:)==0);
80 ql(I)=4;
81 I=fliplr(I);
82 be(:,3*N+1:4*N)=[I(1:end-1);I(2:end)];
83 bel(:,3*N+1:4*N)=4;
84
85 Th=struct('q',q,'me',me,'ql',ql,'mel',zeros(1,nme),'be',be,'bel',bel, ...
86 'nq',size(q,2), ...
87 'nme',size(me,2), ...
88 'nbe',size(be,2), ...
89 'areas',ComputeAreaOpt(q,me),...
90 'lbe',EdgeLengthOpt(be,q));
91