1 /*
2 XLiFE++ is an extended library of finite elements written in C++
3 Copyright (C) 2014 Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin
4
5 This program is free software: you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13 You should have received a copy of the GNU General Public License
14 along with this program. If not, see <http://www.gnu.org/licenses/>.
15 */
16
17 /*!
18 \file SurfMeshQuaCube.cpp
19 \author Y. Lafranche
20 \since 28 Apr 2014
21 \date 28 Apr 2014
22
23 \brief Implementation of xlifepp::subdivision::SurfMeshQuaCube class members and related functions
24 */
25
26 #include "SurfMeshQuaCube.hpp"
27
28 #include <sstream>
29
30 namespace xlifepp {
31 namespace subdivision {
32 using std::map;
33 using std::pair;
34 using std::vector;
35
36 //-------------------------------------------------------------------------------
37 // Constructors, Destructor
38 //-------------------------------------------------------------------------------
39
40 /*!
41 Build a mesh of quadrangles by successive subdivisions from an initial mesh of quadrangles.
42
43 \param rots : rotations to be applied to the cube to get its final position.
44 Each rotation, if any, is defined by an angle in degrees and the number of the
45 absolute axis (1, 2 or 3) around which the rotation is made.
46 Each rotation is applied in turn to the cube starting from the canonical initial
47 position where the cube is centered at the origin and its edges are parallel to the axes.
48
49 \param nboctants : number of octants to be filled
50 nboctants may take a value in [1, 8].
51 For the big cube, the initial mesh may consist of 6 quadrangles instead of
52 24 when this big cube is the juxtatposition of 8 small cubes, which reduces
53 the number of quadrangles of the final mesh by a factor 4.
54 In order to activate this behaviour, nboctants must take the value -8.
55
56 \param nbsubdiv : subdivision level (0 by default)
57 nbsubdiv = 0 corresponds to the initial mesh.
58 For nbsubdiv > 0, each quadrangle is subdivided into 4 quadrangles with the
59 same orientation as the original one: if the first quadrangle is (1,2,3,4), the
60 vector 12x13 defines the exterior normal to the face.
61
62 \param order : order of the quadrangles in the final mesh (1 by default)
63 The default value is 1, which leads to a Q1 mesh, in which case each
64 quadrangle is defined by its 4 vertices.
65 For higher orders, the supplemental vertices correspond to the regular
66 Lagrange mesh.
67
68 \param edgeLength : edge length of the cube (1. by default)
69 \param Center : center of the cube ((0,0,0) by default)
70 \param minVertexNum : minimum number associated to the vertices of the mesh (1 by default)
71 \param minElementNum : minimum number associated to the elements of the mesh (1 by default)
72 */
SurfMeshQuaCube(const vector<pair<real_t,dimen_t>> & rots,const int nboctants,const number_t nbsubdiv,const number_t order,const real_t edgeLength,const Point Center,const number_t minVertexNum,const number_t minElementNum)73 SurfMeshQuaCube::SurfMeshQuaCube(const vector<pair<real_t, dimen_t> >& rots, const int nboctants,
74 const number_t nbsubdiv, const number_t order,
75 const real_t edgeLength, const Point Center,
76 const number_t minVertexNum, const number_t minElementNum)
77 : QuadrangleMesh(nbsubdiv, order, 0, minVertexNum, minElementNum)
78 {
79 // Initialization (nbsubdiv=0)
80 number_t VertexNum, ElementNum;
81 initMesh(rots,nboctants,edgeLength,Center,VertexNum,ElementNum);
82 buildNcheck(VertexNum);
83 }
84
85 //-------------------------------------------------------------------------------
86 // Private member functions
87 //-------------------------------------------------------------------------------
88
89 /*!
90 Create the initial mesh (or "seed" of the mesh) to be subdivided.
91 In each octant, the boundary of the cube is a set of 3 quadrangles sharing a vertex.
92 */
initMesh(const vector<pair<real_t,dimen_t>> & rots,const int nboctants,const real_t edLen,const Point & Center,number_t & VertexNum,number_t & ElementNum)93 void SurfMeshQuaCube::initMesh(const vector<pair<real_t, dimen_t> >& rots, const int nboctants,
94 const real_t edLen, const Point& Center,
95 number_t& VertexNum, number_t& ElementNum){
96 // Define the 27 characteristic points of the cube near the origin Pt[4] (corners,
97 // middle of the edges and faces) in order to define the selected octants.
98 vector<Point> Pt(cubePoints(edLen, rots, Center));
99
100 DefaultGeometry *DG = new DefaultGeometry();
101 // to help numbering: i = index i vector Pt, V[i] = rank in the global list of vertices
102 // 0,1,2,3, 4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27
103 const int V[] = {-1,0,1,2,-1,3,4,5,6,7, 8, 9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25};
104 // Here, Pt[4] is unused and thus, so is V[4]. As a consequence, ranks of the 26 points used
105 // are shifted to take this into account: they are consecutive in the interval [0,25].
106 int iPh = 0;
107 ElementNum = minElementNum_ - 1;
108 VertexNum = minVertexNum_ - 1;
109 number_t mVN = minVertexNum_;
110 refnum_t sBx, sBy, sBz, sDOM, sIx, sIy, sIz;
111 switch (int(std::abs(float(nboctants)))) {
112 default:
113 // Whole cube = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
114 case 0:
115 title_ = "Cube - Quadrangle mesh over the cube";
116 { // Definition of 1 subdomain (no boundary on the surface of the cube, no interface)
117 number_t B_patch[] = {0}, B_dimen[] = {0};
118 number_t I_patch[] = {0}, I_dimen[] = {0};
119 number_t D_patch[] = {1}, D_dimen[] = {2};
120 PatchGeometry *P_geom[] = {DG}; // patch # 1
121 TG_ = TopoGeom( 0,B_patch,B_dimen,
122 0,I_patch,I_dimen,
123 sizeof(D_patch)/sizeof(number_t),D_patch,D_dimen, P_geom);
124
125 // Description of the subdomain patches (# 1)
126 TG_.setDescription(++iPh) = "The surface of the cube"; sDOM = TG_.sigma(iPh);
127
128 // 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27
129 const int W[]={-1,-1,-1,-1,-1,-1, 0,-1,-1,-1,-1, 1,-1, 2,-1,-1,-1, 3,-1,-1, 4,-1,-1,-1, 5, 6,-1, 7};
130 // Definition of the initial vertices, taking into account the boundaries,
131 // interfaces and subdomains they belong to.
132 listV_.push_back(Vertex(W[ 6]+mVN,sDOM,Pt[ 6]));
133 listV_.push_back(Vertex(W[11]+mVN,sDOM,Pt[11]));
134 listV_.push_back(Vertex(W[13]+mVN,sDOM,Pt[13]));
135 listV_.push_back(Vertex(W[17]+mVN,sDOM,Pt[17]));
136 listV_.push_back(Vertex(W[20]+mVN,sDOM,Pt[20]));
137 listV_.push_back(Vertex(W[24]+mVN,sDOM,Pt[24]));
138 listV_.push_back(Vertex(W[25]+mVN,sDOM,Pt[25]));
139 listV_.push_back(Vertex(W[27]+mVN,sDOM,Pt[27]));
140 VertexNum += 8;
141 listT_.push_back(Quadrangle(++ElementNum,W[27],W[20],W[6],W[17]));
142 listT_.push_back(Quadrangle(++ElementNum,W[24],W[25],W[13],W[11]));
143 listT_.push_back(Quadrangle(++ElementNum,W[24],W[11],W[6],W[20]));
144 listT_.push_back(Quadrangle(++ElementNum,W[27],W[17],W[13],W[25]));
145 listT_.push_back(Quadrangle(++ElementNum,W[13],W[17],W[6],W[11]));
146 listT_.push_back(Quadrangle(++ElementNum,W[27],W[25],W[24],W[20]));
147 }
148 break;
149 // 1 octant (1/8 cube) = = = = = = = = = = = = = = = = = = = = = = = = =
150 case 1:
151 title_ = "Cube - Quadrangle mesh over 1 octant X>0, Y>0, Z>0";
152 { // Definition of 3 boundaries and 1 subdomain
153 number_t B_patch[] = {1,2,3}, B_dimen[] = {1,1,1};
154 number_t I_patch[] = {0}, I_dimen[] = {0};
155 number_t D_patch[] = {4}, D_dimen[] = {2};
156 PatchGeometry *P_geom[] = {DG,DG,DG, DG}; // patch # 1,2,3, 4
157 TG_ = TopoGeom(sizeof(B_patch)/sizeof(number_t),B_patch,B_dimen,
158 0,I_patch,I_dimen,
159 sizeof(D_patch)/sizeof(number_t),D_patch,D_dimen, P_geom);
160 }
161 // Description of the boundary patches (# 1,2,3)
162 TG_.setDescription(++iPh) = "Boundary: edges in YZ plane containing the center"; sBx = TG_.sigma(iPh);
163 TG_.setDescription(++iPh) = "Boundary: edges in XZ plane containing the center"; sBy = TG_.sigma(iPh);
164 TG_.setDescription(++iPh) = "Boundary: edges in XY plane containing the center"; sBz = TG_.sigma(iPh);
165 // Description of the subdomain patches (# 4)
166 TG_.setDescription(++iPh) = "Faces in octant #1"; sDOM = TG_.sigma(iPh);
167
168 // Definition of the initial vertices, taking into account the boundaries,
169 // interfaces and subdomains they belong to.
170 listV_.push_back(Vertex(V[1]+mVN, sBy|sBz|sDOM,Pt[1]));
171 listV_.push_back(Vertex(V[2]+mVN, sBz|sDOM,Pt[2]));
172 listV_.push_back(Vertex(V[3]+mVN,sBx |sBz|sDOM,Pt[3]));
173 listV_.push_back(Vertex(V[5]+mVN, sBy |sDOM,Pt[5]));
174 listV_.push_back(Vertex(V[6]+mVN, sDOM,Pt[6]));
175 listV_.push_back(Vertex(V[7]+mVN,sBx |sDOM,Pt[7]));
176 listV_.push_back(Vertex(V[8]+mVN,sBx|sBy |sDOM,Pt[8]));
177 VertexNum += 7;
178 listT_.push_back(Quadrangle(++ElementNum,V[1],V[2],V[6],V[5]));
179 listT_.push_back(Quadrangle(++ElementNum,V[3],V[7],V[6],V[2]));
180 listT_.push_back(Quadrangle(++ElementNum,V[8],V[5],V[6],V[7]));
181 break;
182 // 2 octants (1/4 cube) = = = = = = = = = = = = = = = = = = = = = = = = =
183 case 2:
184 title_ = "Cube - Quadrangle mesh over 2 octants Y>0, Z>0";
185 { // Definition of 2 boundaries, 1 interface and 1 subdomain
186 number_t B_patch[] = {1,2}, B_dimen[] = {1,1};
187 number_t I_patch[] = {4}, I_dimen[] = {1};
188 number_t D_patch[] = {3}, D_dimen[] = {2};
189 PatchGeometry *P_geom[] = {DG,DG, DG, DG}; // patch # 1,2, 3, 4
190 TG_ = TopoGeom(sizeof(B_patch)/sizeof(number_t),B_patch,B_dimen,
191 sizeof(I_patch)/sizeof(number_t),I_patch,I_dimen,
192 sizeof(D_patch)/sizeof(number_t),D_patch,D_dimen, P_geom);
193 }
194 // Description of the boundary patches (# 1,2)
195 TG_.setDescription(++iPh) = "Boundary: edges in XZ plane containing the center"; sBy = TG_.sigma(iPh);
196 TG_.setDescription(++iPh) = "Boundary: edges in XY plane containing the center"; sBz = TG_.sigma(iPh);
197 // Description of the subdomain patches (# 3)
198 TG_.setDescription(++iPh) = "Faces in octants #1 and 2"; sDOM = TG_.sigma(iPh);
199
200 // Description of the interface patches (# 4)
201 TG_.setDescription(++iPh) = "Interface: YZ plane containing the center"; sIx = TG_.sigma(iPh);
202
203 listV_.push_back(Vertex(V[ 1]+mVN,sBy|sBz|sDOM ,Pt[ 1]));
204 listV_.push_back(Vertex(V[ 2]+mVN, sBz|sDOM ,Pt[ 2]));
205 listV_.push_back(Vertex(V[ 3]+mVN, sBz|sDOM|sIx,Pt[ 3]));
206 listV_.push_back(Vertex(V[ 5]+mVN,sBy |sDOM ,Pt[ 5]));
207 listV_.push_back(Vertex(V[ 6]+mVN, sDOM ,Pt[ 6]));
208 listV_.push_back(Vertex(V[ 7]+mVN, sDOM|sIx,Pt[ 7]));
209 listV_.push_back(Vertex(V[ 8]+mVN,sBy |sDOM|sIx,Pt[ 8]));
210 listV_.push_back(Vertex(V[ 9]+mVN,sBy|sBz|sDOM ,Pt[ 9]));
211 listV_.push_back(Vertex(V[10]+mVN, sBz|sDOM ,Pt[10]));
212 listV_.push_back(Vertex(V[11]+mVN, sDOM ,Pt[11]));
213 listV_.push_back(Vertex(V[12]+mVN,sBy |sDOM ,Pt[12]));
214 VertexNum += 11;
215 listT_.push_back(Quadrangle(++ElementNum,V[1],V[2],V[6],V[5]));
216 listT_.push_back(Quadrangle(++ElementNum,V[3],V[7],V[6],V[2]));
217 listT_.push_back(Quadrangle(++ElementNum,V[8],V[5],V[6],V[7]));
218 listT_.push_back(Quadrangle(++ElementNum,V[10],V[9],V[12],V[11]));
219 listT_.push_back(Quadrangle(++ElementNum,V[10],V[11],V[7],V[3]));
220 listT_.push_back(Quadrangle(++ElementNum,V[12],V[8],V[7],V[11]));
221 break;
222 // 3 octants (3/8 cube) = = = = = = = = = = = = = = = = = = = = = = = = =
223 case 3:
224 title_ = "Cube - Quadrangle mesh over 3 octants Z>0, minus the octant X>0, Y<0";
225 { // Definition of 3 boundaries, 2 interfaces and 1 subdomain
226 number_t B_patch[] = {1,2,3}, B_dimen[] = {1,1,1};
227 number_t I_patch[] = {5,6}, I_dimen[] = {1,1};
228 number_t D_patch[] = {4}, D_dimen[] = {2};
229 PatchGeometry *P_geom[] = {DG,DG,DG, DG, DG,DG}; // patch # 1,2,3, 4, 5,6
230 TG_ = TopoGeom(sizeof(B_patch)/sizeof(number_t),B_patch,B_dimen,
231 sizeof(I_patch)/sizeof(number_t),I_patch,I_dimen,
232 sizeof(D_patch)/sizeof(number_t),D_patch,D_dimen, P_geom);
233 }
234 // Description of the boundary patches (# 1,2,3)
235 TG_.setDescription(++iPh) = "Boundary: edges in YZ plane containing the center"; sBx = TG_.sigma(iPh);
236 TG_.setDescription(++iPh) = "Boundary: edges in XZ plane containing the center"; sBy = TG_.sigma(iPh);
237 TG_.setDescription(++iPh) = "Boundary: edges in XY plane containing the center"; sBz = TG_.sigma(iPh);
238 // Description of the subdomain patches (# 4)
239 TG_.setDescription(++iPh) = "Faces in octants #1, 2 and 3"; sDOM = TG_.sigma(iPh);
240 // Description of the interface patches (# 5,6)
241 TG_.setDescription(++iPh) = "Interface: YZ plane containing the center"; sIx = TG_.sigma(iPh);
242 TG_.setDescription(++iPh) = "Interface: XZ plane containing the center"; sIy = TG_.sigma(iPh);
243
244 listV_.push_back(Vertex(V[ 1]+mVN, sBy|sBz|sDOM ,Pt[ 1]));
245 listV_.push_back(Vertex(V[ 2]+mVN, sBz|sDOM ,Pt[ 2]));
246 listV_.push_back(Vertex(V[ 3]+mVN, sBz|sDOM|sIx ,Pt[ 3]));
247 listV_.push_back(Vertex(V[ 5]+mVN, sBy |sDOM ,Pt[ 5]));
248 listV_.push_back(Vertex(V[ 6]+mVN, sDOM ,Pt[ 6]));
249 listV_.push_back(Vertex(V[ 7]+mVN, sDOM|sIx ,Pt[ 7]));
250 listV_.push_back(Vertex(V[ 8]+mVN,sBx|sBy |sDOM|sIx|sIy,Pt[ 8]));
251 listV_.push_back(Vertex(V[ 9]+mVN, sBz|sDOM |sIy,Pt[ 9]));
252 listV_.push_back(Vertex(V[10]+mVN, sBz|sDOM ,Pt[10]));
253 listV_.push_back(Vertex(V[11]+mVN, sDOM ,Pt[11]));
254 listV_.push_back(Vertex(V[12]+mVN, sDOM |sIy,Pt[12]));
255 listV_.push_back(Vertex(V[13]+mVN, sDOM ,Pt[13]));
256 listV_.push_back(Vertex(V[14]+mVN, sBz|sDOM ,Pt[14]));
257 listV_.push_back(Vertex(V[15]+mVN,sBx |sBz|sDOM ,Pt[15]));
258 listV_.push_back(Vertex(V[16]+mVN,sBx |sDOM ,Pt[16]));
259 VertexNum += 15;
260 listT_.push_back(Quadrangle(++ElementNum,V[1],V[2],V[6],V[5]));
261 listT_.push_back(Quadrangle(++ElementNum,V[3],V[7],V[6],V[2]));
262 listT_.push_back(Quadrangle(++ElementNum,V[8],V[5],V[6],V[7]));
263 listT_.push_back(Quadrangle(++ElementNum,V[10],V[9],V[12],V[11]));
264 listT_.push_back(Quadrangle(++ElementNum,V[10],V[11],V[7],V[3]));
265 listT_.push_back(Quadrangle(++ElementNum,V[12],V[8],V[7],V[11]));
266 listT_.push_back(Quadrangle(++ElementNum,V[9],V[14],V[13],V[12]));
267 listT_.push_back(Quadrangle(++ElementNum,V[15],V[16],V[13],V[14]));
268 listT_.push_back(Quadrangle(++ElementNum,V[13],V[16],V[8],V[12]));
269 break;
270 // 4 octants (1/2 cube) = = = = = = = = = = = = = = = = = = = = = = = = =
271 case 4:
272 title_ = "Cube - Quadrangle mesh over 4 octants Z>0";
273 { // Definition of 1 boundary, 2 interfaces and 1 subdomain
274 number_t B_patch[] = {1}, B_dimen[] = {1};
275 number_t I_patch[] = {3,4}, I_dimen[] = {1,1};
276 number_t D_patch[] = {2}, D_dimen[] = {2};
277 PatchGeometry *P_geom[] = {DG, DG, DG,DG}; // patch # 1, 2, 3,4
278 TG_ = TopoGeom(sizeof(B_patch)/sizeof(number_t),B_patch,B_dimen,
279 sizeof(I_patch)/sizeof(number_t),I_patch,I_dimen,
280 sizeof(D_patch)/sizeof(number_t),D_patch,D_dimen, P_geom);
281 }
282 // Description of the boundary patches (# 1)
283 TG_.setDescription(++iPh) = "Boundary: edges in XY plane containing the center"; sBz = TG_.sigma(iPh);
284 // Description of the subdomain patches (# 2)
285 TG_.setDescription(++iPh) = "Faces in octants #1 to 4, in the half-space Z>0"; sDOM = TG_.sigma(iPh);
286 // Description of the interface patches (# 3,4)
287 TG_.setDescription(++iPh) = "Interface: YZ plane containing the center"; sIx = TG_.sigma(iPh);
288 TG_.setDescription(++iPh) = "Interface: XZ plane containing the center"; sIy = TG_.sigma(iPh);
289
290 listV_.push_back(Vertex(V[ 1]+mVN,sBz|sDOM |sIy,Pt[ 1]));
291 listV_.push_back(Vertex(V[ 2]+mVN,sBz|sDOM ,Pt[ 2]));
292 listV_.push_back(Vertex(V[ 3]+mVN,sBz|sDOM|sIx ,Pt[ 3]));
293 listV_.push_back(Vertex(V[ 5]+mVN, sDOM |sIy,Pt[ 5]));
294 listV_.push_back(Vertex(V[ 6]+mVN, sDOM ,Pt[ 6]));
295 listV_.push_back(Vertex(V[ 7]+mVN, sDOM|sIx ,Pt[ 7]));
296 listV_.push_back(Vertex(V[ 8]+mVN, sDOM|sIx|sIy,Pt[ 8]));
297 listV_.push_back(Vertex(V[ 9]+mVN,sBz|sDOM |sIy,Pt[ 9]));
298 listV_.push_back(Vertex(V[10]+mVN,sBz|sDOM ,Pt[10]));
299 listV_.push_back(Vertex(V[11]+mVN, sDOM ,Pt[11]));
300 listV_.push_back(Vertex(V[12]+mVN, sDOM |sIy,Pt[12]));
301 listV_.push_back(Vertex(V[13]+mVN, sDOM ,Pt[13]));
302 listV_.push_back(Vertex(V[14]+mVN,sBz|sDOM ,Pt[14]));
303 listV_.push_back(Vertex(V[15]+mVN,sBz|sDOM|sIx ,Pt[15]));
304 listV_.push_back(Vertex(V[16]+mVN, sDOM|sIx ,Pt[16]));
305 listV_.push_back(Vertex(V[17]+mVN, sDOM ,Pt[17]));
306 listV_.push_back(Vertex(V[18]+mVN,sBz|sDOM ,Pt[18]));
307 VertexNum += 17;
308 listT_.push_back(Quadrangle(++ElementNum,V[1],V[2],V[6],V[5]));
309 listT_.push_back(Quadrangle(++ElementNum,V[3],V[7],V[6],V[2]));
310 listT_.push_back(Quadrangle(++ElementNum,V[8],V[5],V[6],V[7]));
311 listT_.push_back(Quadrangle(++ElementNum,V[10],V[9],V[12],V[11]));
312 listT_.push_back(Quadrangle(++ElementNum,V[10],V[11],V[7],V[3]));
313 listT_.push_back(Quadrangle(++ElementNum,V[12],V[8],V[7],V[11]));
314 listT_.push_back(Quadrangle(++ElementNum,V[9],V[14],V[13],V[12]));
315 listT_.push_back(Quadrangle(++ElementNum,V[15],V[16],V[13],V[14]));
316 listT_.push_back(Quadrangle(++ElementNum,V[13],V[16],V[8],V[12]));
317 listT_.push_back(Quadrangle(++ElementNum,V[18],V[1],V[5],V[17]));
318 listT_.push_back(Quadrangle(++ElementNum,V[18],V[17],V[16],V[15]));
319 listT_.push_back(Quadrangle(++ElementNum,V[16],V[17],V[5],V[8]));
320 break;
321 // 5 octants (5/8 cube) = = = = = = = = = = = = = = = = = = = = = = = = =
322 case 5:
323 title_ = "Cube - Quadrangle mesh over 5 octants, 4 with Z>0, plus the octant X>0, Y>0, Z<0";
324 { // Definition of 3 boundaries, 3 interfaces and 1 subdomain
325 number_t B_patch[] = {1,2,3}, B_dimen[] = {1,1,1};
326 number_t I_patch[] = {5,6,7}, I_dimen[] = {1,1,1};
327 number_t D_patch[] = {4}, D_dimen[] = {2};
328 number_t nbBound = sizeof(B_patch)/sizeof(number_t);
329 PatchGeometry *P_geom[] = {DG,DG,DG, DG, DG,DG,DG}; // patch # 1,2,3, 4, 5,6,7
330 TG_ = TopoGeom( nbBound,B_patch,B_dimen,
331 sizeof(I_patch)/sizeof(number_t),I_patch,I_dimen,
332 sizeof(D_patch)/sizeof(number_t),D_patch,D_dimen, P_geom);
333 }
334 // Description of the boundary patches (# 1,2,3)
335 TG_.setDescription(++iPh) = "Boundary: edges in YZ plane containing the center"; sBx = TG_.sigma(iPh);
336 TG_.setDescription(++iPh) = "Boundary: edges in XZ plane containing the center"; sBy = TG_.sigma(iPh);
337 TG_.setDescription(++iPh) = "Boundary: edges in XY plane containing the center"; sBz = TG_.sigma(iPh);
338 // Description of the subdomain patches (# 4)
339 TG_.setDescription(++iPh) = "Faces in octants #1 to 5"; sDOM = TG_.sigma(iPh);
340 // Description of the interface patches (# 5,6,7)
341 TG_.setDescription(++iPh) = "Interface: YZ plane containing the center"; sIx = TG_.sigma(iPh);
342 TG_.setDescription(++iPh) = "Interface: XZ plane containing the center"; sIy = TG_.sigma(iPh);
343 TG_.setDescription(++iPh) = "Interface: XY plane containing the center"; sIz = TG_.sigma(iPh);
344
345 listV_.push_back(Vertex(V[ 1]+mVN, sBy|sBz|sDOM |sIy|sIz,Pt[ 1]));
346 listV_.push_back(Vertex(V[ 2]+mVN, sDOM |sIz,Pt[ 2]));
347 listV_.push_back(Vertex(V[ 3]+mVN,sBx |sBz|sDOM|sIx |sIz,Pt[ 3]));
348 listV_.push_back(Vertex(V[ 5]+mVN, sDOM |sIy ,Pt[ 5]));
349 listV_.push_back(Vertex(V[ 6]+mVN, sDOM ,Pt[ 6]));
350 listV_.push_back(Vertex(V[ 7]+mVN, sDOM|sIx ,Pt[ 7]));
351 listV_.push_back(Vertex(V[ 8]+mVN, sDOM|sIx|sIy ,Pt[ 8]));
352 listV_.push_back(Vertex(V[ 9]+mVN, sBz|sDOM |sIy ,Pt[ 9]));
353 listV_.push_back(Vertex(V[10]+mVN, sBz|sDOM ,Pt[10]));
354 listV_.push_back(Vertex(V[11]+mVN, sDOM ,Pt[11]));
355 listV_.push_back(Vertex(V[12]+mVN, sDOM |sIy ,Pt[12]));
356 listV_.push_back(Vertex(V[13]+mVN, sDOM ,Pt[13]));
357 listV_.push_back(Vertex(V[14]+mVN, sBz|sDOM ,Pt[14]));
358 listV_.push_back(Vertex(V[15]+mVN, sBz|sDOM|sIx ,Pt[15]));
359 listV_.push_back(Vertex(V[16]+mVN, sDOM|sIx ,Pt[16]));
360 listV_.push_back(Vertex(V[17]+mVN, sDOM ,Pt[17]));
361 listV_.push_back(Vertex(V[18]+mVN, sBz|sDOM ,Pt[18]));
362 listV_.push_back(Vertex(V[19]+mVN, sBy |sDOM ,Pt[19]));
363 listV_.push_back(Vertex(V[20]+mVN, sDOM ,Pt[20]));
364 listV_.push_back(Vertex(V[21]+mVN,sBx |sDOM ,Pt[21]));
365 listV_.push_back(Vertex(V[22]+mVN,sBx|sBy |sDOM ,Pt[22]));
366 VertexNum += 21;
367 // Cubes with Z>0
368 listT_.push_back(Quadrangle(++ElementNum,V[1],V[2],V[6],V[5]));
369 listT_.push_back(Quadrangle(++ElementNum,V[3],V[7],V[6],V[2]));
370 listT_.push_back(Quadrangle(++ElementNum,V[8],V[5],V[6],V[7]));
371 listT_.push_back(Quadrangle(++ElementNum,V[10],V[9],V[12],V[11]));
372 listT_.push_back(Quadrangle(++ElementNum,V[10],V[11],V[7],V[3]));
373 listT_.push_back(Quadrangle(++ElementNum,V[12],V[8],V[7],V[11]));
374 listT_.push_back(Quadrangle(++ElementNum,V[9],V[14],V[13],V[12]));
375 listT_.push_back(Quadrangle(++ElementNum,V[15],V[16],V[13],V[14]));
376 listT_.push_back(Quadrangle(++ElementNum,V[13],V[16],V[8],V[12]));
377 listT_.push_back(Quadrangle(++ElementNum,V[18],V[1],V[5],V[17]));
378 listT_.push_back(Quadrangle(++ElementNum,V[18],V[17],V[16],V[15]));
379 listT_.push_back(Quadrangle(++ElementNum,V[16],V[17],V[5],V[8]));
380 // Cubes with Z<0
381 listT_.push_back(Quadrangle(++ElementNum,V[19],V[20],V[2],V[1]));
382 listT_.push_back(Quadrangle(++ElementNum,V[21],V[3],V[2],V[20]));
383 listT_.push_back(Quadrangle(++ElementNum,V[19],V[22],V[21],V[20]));
384 break;
385 // 6 octants (3/4 cube) = = = = = = = = = = = = = = = = = = = = = = = = =
386 case 6:
387 title_ = "Cube - Quadrangle mesh over 6 octants, 4 with Z>0, plus 2 octants Y>0, Z<0";
388 { // Definition of 2 boundaries, 3 interfaces and 1 subdomain
389 number_t B_patch[] = {1,2}, B_dimen[] = {1,1};
390 number_t I_patch[] = {4,5,6}, I_dimen[] = {1,1,1};
391 number_t D_patch[] = {3}, D_dimen[] = {2};
392 PatchGeometry *P_geom[] = {DG,DG, DG, DG,DG,DG}; // patch # 1,2, 3, 4,5,6
393 TG_ = TopoGeom(sizeof(B_patch)/sizeof(number_t),B_patch,B_dimen,
394 sizeof(I_patch)/sizeof(number_t),I_patch,I_dimen,
395 sizeof(D_patch)/sizeof(number_t),D_patch,D_dimen, P_geom);
396 }
397 // Description of the boundary patches
398 TG_.setDescription(++iPh) = "Boundary: edges in XZ plane containing the center"; sBy = TG_.sigma(iPh);
399 TG_.setDescription(++iPh) = "Boundary: edges in XY plane containing the center"; sBz = TG_.sigma(iPh);
400 // Description of the subdomain patches
401 TG_.setDescription(++iPh) = "Faces in octants #1 to 6"; sDOM = TG_.sigma(iPh);
402 // Description of the interface patches
403 TG_.setDescription(++iPh) = "Interface: YZ plane containing the center"; sIx = TG_.sigma(iPh);
404 TG_.setDescription(++iPh) = "Interface: XZ plane containing the center"; sIy = TG_.sigma(iPh);
405 TG_.setDescription(++iPh) = "Interface: XY plane containing the center"; sIz = TG_.sigma(iPh);
406
407 listV_.push_back(Vertex(V[ 1]+mVN,sBy|sBz|sDOM |sIy|sIz,Pt[ 1]));
408 listV_.push_back(Vertex(V[ 2]+mVN, sDOM |sIz,Pt[ 2]));
409 listV_.push_back(Vertex(V[ 3]+mVN, sDOM|sIx |sIz,Pt[ 3]));
410 listV_.push_back(Vertex(V[ 5]+mVN, sDOM |sIy ,Pt[ 5]));
411 listV_.push_back(Vertex(V[ 6]+mVN, sDOM ,Pt[ 6]));
412 listV_.push_back(Vertex(V[ 7]+mVN, sDOM|sIx ,Pt[ 7]));
413 listV_.push_back(Vertex(V[ 8]+mVN, sDOM|sIx|sIy ,Pt[ 8]));
414 listV_.push_back(Vertex(V[ 9]+mVN,sBy|sBz|sDOM |sIy|sIz,Pt[ 9]));
415 listV_.push_back(Vertex(V[10]+mVN, sDOM |sIz,Pt[10]));
416 listV_.push_back(Vertex(V[11]+mVN, sDOM ,Pt[11]));
417 listV_.push_back(Vertex(V[12]+mVN, sDOM |sIy ,Pt[12]));
418 listV_.push_back(Vertex(V[13]+mVN, sDOM ,Pt[13]));
419 listV_.push_back(Vertex(V[14]+mVN, sBz|sDOM ,Pt[14]));
420 listV_.push_back(Vertex(V[15]+mVN, sBz|sDOM|sIx ,Pt[15]));
421 listV_.push_back(Vertex(V[16]+mVN, sDOM|sIx ,Pt[16]));
422 listV_.push_back(Vertex(V[17]+mVN, sDOM ,Pt[17]));
423 listV_.push_back(Vertex(V[18]+mVN, sBz|sDOM ,Pt[18]));
424 listV_.push_back(Vertex(V[19]+mVN,sBy |sDOM ,Pt[19]));
425 listV_.push_back(Vertex(V[20]+mVN, sDOM ,Pt[20]));
426 listV_.push_back(Vertex(V[21]+mVN, sDOM|sIx ,Pt[21]));
427 listV_.push_back(Vertex(V[22]+mVN,sBy |sDOM|sIx ,Pt[22]));
428 listV_.push_back(Vertex(V[23]+mVN,sBy |sDOM ,Pt[23]));
429 listV_.push_back(Vertex(V[24]+mVN, sDOM ,Pt[24]));
430 VertexNum += 23;
431 // Cubes with Z>0
432 listT_.push_back(Quadrangle(++ElementNum,V[1],V[2],V[6],V[5]));
433 listT_.push_back(Quadrangle(++ElementNum,V[3],V[7],V[6],V[2]));
434 listT_.push_back(Quadrangle(++ElementNum,V[8],V[5],V[6],V[7]));
435 listT_.push_back(Quadrangle(++ElementNum,V[10],V[9],V[12],V[11]));
436 listT_.push_back(Quadrangle(++ElementNum,V[10],V[11],V[7],V[3]));
437 listT_.push_back(Quadrangle(++ElementNum,V[12],V[8],V[7],V[11]));
438 listT_.push_back(Quadrangle(++ElementNum,V[9],V[14],V[13],V[12]));
439 listT_.push_back(Quadrangle(++ElementNum,V[15],V[16],V[13],V[14]));
440 listT_.push_back(Quadrangle(++ElementNum,V[13],V[16],V[8],V[12]));
441 listT_.push_back(Quadrangle(++ElementNum,V[18],V[1],V[5],V[17]));
442 listT_.push_back(Quadrangle(++ElementNum,V[18],V[17],V[16],V[15]));
443 listT_.push_back(Quadrangle(++ElementNum,V[16],V[17],V[5],V[8]));
444 // Cubes with Z<0
445 listT_.push_back(Quadrangle(++ElementNum,V[19],V[20],V[2],V[1]));
446 listT_.push_back(Quadrangle(++ElementNum,V[21],V[3],V[2],V[20]));
447 listT_.push_back(Quadrangle(++ElementNum,V[19],V[22],V[21],V[20]));
448 listT_.push_back(Quadrangle(++ElementNum,V[24],V[23],V[9],V[10]));
449 listT_.push_back(Quadrangle(++ElementNum,V[24],V[10],V[3],V[21]));
450 listT_.push_back(Quadrangle(++ElementNum,V[22],V[23],V[24],V[21]));
451 break;
452 // 7 octants (7/8 cube) = = = = = = = = = = = = = = = = = = = = = = = = =
453 case 7:
454 title_ = "Cube - Quadrangle mesh over 7 octants, 8 octants minus the octant X>0, Y<0, Z<0";
455 { // Definition of 3 boundaries, 3 interfaces and 1 subdomain
456 number_t B_patch[] = {1,2,3}, B_dimen[] = {1,1,1};
457 number_t I_patch[] = {5,6,7}, I_dimen[] = {1,1,1};
458 number_t D_patch[] = {4}, D_dimen[] = {2};
459 number_t nbBound = sizeof(B_patch)/sizeof(number_t);
460 PatchGeometry *P_geom[] = {DG,DG,DG, DG, DG,DG,DG}; // patch # 1,2,3, 4, 5,6,7
461 TG_ = TopoGeom( nbBound,B_patch,B_dimen,
462 sizeof(I_patch)/sizeof(number_t),I_patch,I_dimen,
463 sizeof(D_patch)/sizeof(number_t),D_patch,D_dimen, P_geom);
464 }
465 // Description of the boundary patches (# 1,2,3)
466 TG_.setDescription(++iPh) = "Boundary: edges in YZ plane containing the center"; sBx = TG_.sigma(iPh);
467 TG_.setDescription(++iPh) = "Boundary: edges in XZ plane containing the center"; sBy = TG_.sigma(iPh);
468 TG_.setDescription(++iPh) = "Boundary: edges in XY plane containing the center"; sBz = TG_.sigma(iPh);
469 // Description of the subdomain patches (# 4)
470 TG_.setDescription(++iPh) = "Faces in octants #1 to 7"; sDOM = TG_.sigma(iPh);
471 // Description of the interface patches (# 5,6,7)
472 TG_.setDescription(++iPh) = "Interface: YZ plane containing the center"; sIx = TG_.sigma(iPh);
473 TG_.setDescription(++iPh) = "Interface: XZ plane containing the center"; sIy = TG_.sigma(iPh);
474 TG_.setDescription(++iPh) = "Interface: XY plane containing the center"; sIz = TG_.sigma(iPh);
475
476 listV_.push_back(Vertex(V[ 1]+mVN, sBy|sBz|sDOM |sIy|sIz,Pt[ 1]));
477 listV_.push_back(Vertex(V[ 2]+mVN, sDOM |sIz,Pt[ 2]));
478 listV_.push_back(Vertex(V[ 3]+mVN, sDOM|sIx |sIz,Pt[ 3]));
479 listV_.push_back(Vertex(V[ 5]+mVN, sDOM |sIy ,Pt[ 5]));
480 listV_.push_back(Vertex(V[ 6]+mVN, sDOM ,Pt[ 6]));
481 listV_.push_back(Vertex(V[ 7]+mVN, sDOM|sIx ,Pt[ 7]));
482 listV_.push_back(Vertex(V[ 8]+mVN, sDOM|sIx|sIy ,Pt[ 8]));
483 listV_.push_back(Vertex(V[ 9]+mVN, sDOM |sIy|sIz,Pt[ 9]));
484 listV_.push_back(Vertex(V[10]+mVN, sDOM |sIz,Pt[10]));
485 listV_.push_back(Vertex(V[11]+mVN, sDOM ,Pt[11]));
486 listV_.push_back(Vertex(V[12]+mVN, sDOM |sIy ,Pt[12]));
487 listV_.push_back(Vertex(V[13]+mVN, sDOM ,Pt[13]));
488 listV_.push_back(Vertex(V[14]+mVN, sDOM |sIz,Pt[14]));
489 listV_.push_back(Vertex(V[15]+mVN,sBx |sBz|sDOM|sIx |sIz,Pt[15]));
490 listV_.push_back(Vertex(V[16]+mVN, sDOM|sIx ,Pt[16]));
491 listV_.push_back(Vertex(V[17]+mVN, sDOM ,Pt[17]));
492 listV_.push_back(Vertex(V[18]+mVN, sBz|sDOM ,Pt[18]));
493 listV_.push_back(Vertex(V[19]+mVN, sBy |sDOM ,Pt[19]));
494 listV_.push_back(Vertex(V[20]+mVN, sDOM ,Pt[20]));
495 listV_.push_back(Vertex(V[21]+mVN, sDOM|sIx ,Pt[21]));
496 listV_.push_back(Vertex(V[22]+mVN,sBx|sBy |sDOM|sIx|sIy ,Pt[22]));
497 listV_.push_back(Vertex(V[23]+mVN, sDOM |sIy ,Pt[23]));
498 listV_.push_back(Vertex(V[24]+mVN, sDOM ,Pt[24]));
499 listV_.push_back(Vertex(V[25]+mVN, sDOM ,Pt[25]));
500 listV_.push_back(Vertex(V[26]+mVN,sBx |sDOM ,Pt[26]));
501 VertexNum += 25;
502 // Cubes with Z>0
503 listT_.push_back(Quadrangle(++ElementNum,V[1],V[2],V[6],V[5]));
504 listT_.push_back(Quadrangle(++ElementNum,V[3],V[7],V[6],V[2]));
505 listT_.push_back(Quadrangle(++ElementNum,V[8],V[5],V[6],V[7]));
506 listT_.push_back(Quadrangle(++ElementNum,V[10],V[9],V[12],V[11]));
507 listT_.push_back(Quadrangle(++ElementNum,V[10],V[11],V[7],V[3]));
508 listT_.push_back(Quadrangle(++ElementNum,V[12],V[8],V[7],V[11]));
509 listT_.push_back(Quadrangle(++ElementNum,V[9],V[14],V[13],V[12]));
510 listT_.push_back(Quadrangle(++ElementNum,V[15],V[16],V[13],V[14]));
511 listT_.push_back(Quadrangle(++ElementNum,V[13],V[16],V[8],V[12]));
512 listT_.push_back(Quadrangle(++ElementNum,V[18],V[1],V[5],V[17]));
513 listT_.push_back(Quadrangle(++ElementNum,V[18],V[17],V[16],V[15]));
514 listT_.push_back(Quadrangle(++ElementNum,V[16],V[17],V[5],V[8]));
515 // Cubes with Z<0
516 listT_.push_back(Quadrangle(++ElementNum,V[19],V[20],V[2],V[1]));
517 listT_.push_back(Quadrangle(++ElementNum,V[21],V[3],V[2],V[20]));
518 listT_.push_back(Quadrangle(++ElementNum,V[19],V[22],V[21],V[20]));
519 listT_.push_back(Quadrangle(++ElementNum,V[24],V[23],V[9],V[10]));
520 listT_.push_back(Quadrangle(++ElementNum,V[24],V[10],V[3],V[21]));
521 listT_.push_back(Quadrangle(++ElementNum,V[22],V[23],V[24],V[21]));
522 listT_.push_back(Quadrangle(++ElementNum,V[23],V[25],V[14],V[9]));
523 listT_.push_back(Quadrangle(++ElementNum,V[26],V[15],V[14],V[25]));
524 listT_.push_back(Quadrangle(++ElementNum,V[26],V[25],V[23],V[22]));
525 break;
526 // 8 octants (cube) = = = = = = = = = = = = = = = = = = = = = = = = = = =
527 case 8:
528 title_ = "Cube - Quadrangle mesh over 8 octants";
529 { // Definition of 3 interfaces and 1 subdomain (no boundary on the surface of the cube)
530 number_t B_patch[] = {0}, B_dimen[] = {0};
531 number_t I_patch[] = {2,3,4}, I_dimen[] = {1,1,1};
532 number_t D_patch[] = {1}, D_dimen[] = {2};
533 PatchGeometry *P_geom[] = {DG, DG,DG,DG}; // patch # 1, 2,3,4
534 TG_ = TopoGeom( 0,B_patch,B_dimen,
535 sizeof(I_patch)/sizeof(number_t),I_patch,I_dimen,
536 sizeof(D_patch)/sizeof(number_t),D_patch,D_dimen, P_geom);
537 }
538 // Description of the subdomain patches (# 1)
539 TG_.setDescription(++iPh) = "The surface of the cube"; sDOM = TG_.sigma(iPh);
540 // Description of the interface patches (# 2,3,4)
541 TG_.setDescription(++iPh) = "Interface: YZ plane containing the center"; sIx = TG_.sigma(iPh);
542 TG_.setDescription(++iPh) = "Interface: XZ plane containing the center"; sIy = TG_.sigma(iPh);
543 TG_.setDescription(++iPh) = "Interface: XY plane containing the center"; sIz = TG_.sigma(iPh);
544
545 listV_.push_back(Vertex(V[ 1]+mVN,sDOM |sIy|sIz,Pt[ 1]));
546 listV_.push_back(Vertex(V[ 2]+mVN,sDOM |sIz,Pt[ 2]));
547 listV_.push_back(Vertex(V[ 3]+mVN,sDOM|sIx |sIz,Pt[ 3]));
548 listV_.push_back(Vertex(V[ 5]+mVN,sDOM |sIy ,Pt[ 5]));
549 listV_.push_back(Vertex(V[ 6]+mVN,sDOM ,Pt[ 6]));
550 listV_.push_back(Vertex(V[ 7]+mVN,sDOM|sIx ,Pt[ 7]));
551 listV_.push_back(Vertex(V[ 8]+mVN,sDOM|sIx|sIy ,Pt[ 8]));
552 listV_.push_back(Vertex(V[ 9]+mVN,sDOM |sIy|sIz,Pt[ 9]));
553 listV_.push_back(Vertex(V[10]+mVN,sDOM |sIz,Pt[10]));
554 listV_.push_back(Vertex(V[11]+mVN,sDOM ,Pt[11]));
555 listV_.push_back(Vertex(V[12]+mVN,sDOM |sIy ,Pt[12]));
556 listV_.push_back(Vertex(V[13]+mVN,sDOM ,Pt[13]));
557 listV_.push_back(Vertex(V[14]+mVN,sDOM |sIz,Pt[14]));
558 listV_.push_back(Vertex(V[15]+mVN,sDOM|sIx |sIz,Pt[15]));
559 listV_.push_back(Vertex(V[16]+mVN,sDOM|sIx ,Pt[16]));
560 listV_.push_back(Vertex(V[17]+mVN,sDOM ,Pt[17]));
561 listV_.push_back(Vertex(V[18]+mVN,sDOM |sIz,Pt[18]));
562 listV_.push_back(Vertex(V[19]+mVN,sDOM |sIy ,Pt[19]));
563 listV_.push_back(Vertex(V[20]+mVN,sDOM ,Pt[20]));
564 listV_.push_back(Vertex(V[21]+mVN,sDOM|sIx ,Pt[21]));
565 listV_.push_back(Vertex(V[22]+mVN,sDOM|sIx|sIy ,Pt[22]));
566 listV_.push_back(Vertex(V[23]+mVN,sDOM |sIy ,Pt[23]));
567 listV_.push_back(Vertex(V[24]+mVN,sDOM ,Pt[24]));
568 listV_.push_back(Vertex(V[25]+mVN,sDOM ,Pt[25]));
569 listV_.push_back(Vertex(V[26]+mVN,sDOM|sIx ,Pt[26]));
570 listV_.push_back(Vertex(V[27]+mVN,sDOM ,Pt[27]));
571 VertexNum += 26;
572 // Cubes with Z>0
573 listT_.push_back(Quadrangle(++ElementNum,V[1],V[2],V[6],V[5]));
574 listT_.push_back(Quadrangle(++ElementNum,V[3],V[7],V[6],V[2]));
575 listT_.push_back(Quadrangle(++ElementNum,V[8],V[5],V[6],V[7]));
576 listT_.push_back(Quadrangle(++ElementNum,V[10],V[9],V[12],V[11]));
577 listT_.push_back(Quadrangle(++ElementNum,V[10],V[11],V[7],V[3]));
578 listT_.push_back(Quadrangle(++ElementNum,V[12],V[8],V[7],V[11]));
579 listT_.push_back(Quadrangle(++ElementNum,V[9],V[14],V[13],V[12]));
580 listT_.push_back(Quadrangle(++ElementNum,V[15],V[16],V[13],V[14]));
581 listT_.push_back(Quadrangle(++ElementNum,V[13],V[16],V[8],V[12]));
582 listT_.push_back(Quadrangle(++ElementNum,V[18],V[1],V[5],V[17]));
583 listT_.push_back(Quadrangle(++ElementNum,V[18],V[17],V[16],V[15]));
584 listT_.push_back(Quadrangle(++ElementNum,V[16],V[17],V[5],V[8]));
585 // Cubes with Z<0
586 listT_.push_back(Quadrangle(++ElementNum,V[19],V[20],V[2],V[1]));
587 listT_.push_back(Quadrangle(++ElementNum,V[21],V[3],V[2],V[20]));
588 listT_.push_back(Quadrangle(++ElementNum,V[19],V[22],V[21],V[20]));
589 listT_.push_back(Quadrangle(++ElementNum,V[24],V[23],V[9],V[10]));
590 listT_.push_back(Quadrangle(++ElementNum,V[24],V[10],V[3],V[21]));
591 listT_.push_back(Quadrangle(++ElementNum,V[22],V[23],V[24],V[21]));
592 listT_.push_back(Quadrangle(++ElementNum,V[23],V[25],V[14],V[9]));
593 listT_.push_back(Quadrangle(++ElementNum,V[26],V[15],V[14],V[25]));
594 listT_.push_back(Quadrangle(++ElementNum,V[26],V[25],V[23],V[22]));
595 listT_.push_back(Quadrangle(++ElementNum,V[27],V[19],V[1],V[18]));
596 listT_.push_back(Quadrangle(++ElementNum,V[27],V[18],V[15],V[26]));
597 listT_.push_back(Quadrangle(++ElementNum,V[27],V[26],V[22],V[19]));
598 break;
599 }
600 }
601
602 } // end of namespace subdivision
603 } // end of namespace xlifepp
604