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