| OptFEM2D Toolbox for Matlab V1.2b1 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 %
45 % OptFEM2DP1 [V1.2b1] - Copyright (C) 2013 CJS (LAGA)
46 %
47 % This file is part of OptFEM2DP1.
48 % OptFEM2DP1 is free software: you can redistribute it and/or modify
49 % it under the terms of the GNU General Public License as published by
50 % the Free Software Foundation, either version 3 of the License, or
51 % (at your option) any later version.
52 %
53 % OptFEM2DP1 is distributed in the hope that it will be useful,
54 % but WITHOUT ANY WARRANTY; without even the implied warranty of
55 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
56 % GNU General Public License for more details.
57 %
58 % You should have received a copy of the GNU General Public License
59 % along with this program. If not, see <http://www.gnu.org/licenses/>.
60 h=1/N;t=0:h:1;
61 if isOctave()
62 mesh=msh2m_structured_mesh(t,t,1,1:4); % package msh
63 me=mesh.t(1:3,:);
64 q=mesh.p;
65 else
66 [x,y] = meshgrid(t,t);
67 q=[x(:) y(:)]';
68 tri = delaunay(x,y);
69 me=tri';
70 end
71
72 nq=length(q);
73 nme=length(me);
74
75 ql=zeros(1,nq);
76 be=zeros(2,4*N);
77 bel=zeros(1,4*N);
78 % Points du bord
79 ql=zeros(1,nq);
80 I=find(q(2,:)==0);
81 ql(I)=1;
82 be(:,1:N)=[I(1:end-1);I(2:end)];
83 bel(:,1:N)=1;
84 I=find(q(1,:)==1);
85 ql(I)=2;
86 be(:,N+1:2*N)=[I(1:end-1);I(2:end)];
87 bel(:,N+1:2*N)=2;
88 I=find(q(2,:)==1);
89 ql(I)=3;
90 I=fliplr(I);
91 be(:,2*N+1:3*N)=[I(1:end-1);I(2:end)];
92 bel(:,2*N+1:3*N)=3;
93 I=find(q(1,:)==0);
94 ql(I)=4;
95 I=fliplr(I);
96 be(:,3*N+1:4*N)=[I(1:end-1);I(2:end)];
97 bel(:,3*N+1:4*N)=4;
98
99 Th=struct('q',q,'me',me,'ql',ql,'mel',zeros(1,nme),'be',be,'bel',bel, ...
100 'nq',size(q,2), ...
101 'nme',size(me,2), ...
102 'nbe',size(be,2), ...
103 'areas',ComputeAreaOpt(q,me),...
104 'lbe',EdgeLengthOpt(be,q));
105