/[escript]/trunk/finley/src/Mesh_hex20.cpp
ViewVC logotype

Contents of /trunk/finley/src/Mesh_hex20.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6939 - (show annotations)
Mon Jan 20 03:37:18 2020 UTC (4 months, 1 week ago) by uqaeller
File size: 31354 byte(s)
Updated the copyright header.


1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2020 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Apache License, version 2.0
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014-2017 by Centre for Geoscience Computing (GeoComp)
14 * Development from 2019 by School of Earth and Environmental Sciences
15 **
16 *****************************************************************************/
17
18
19 /****************************************************************************
20
21 Finley: generates rectangular meshes
22
23 Generates a NE0 x NE1 x NE2 mesh with second order elements (Hex20) in the
24 brick [0,l0] x [0,l1] x [0,l2].
25 order is the desired accuracy of the integration scheme.
26
27 *****************************************************************************/
28
29 #include "FinleyDomain.h"
30
31 #include <escript/index.h>
32
33 #define MAX3(_arg1_,_arg2_,_arg3_) std::max(_arg1_,std::max(_arg2_,_arg3_))
34
35 using escript::DataTypes::real_t;
36
37 namespace finley {
38
39 escript::Domain_ptr FinleyDomain::createHex20(dim_t NE0, dim_t NE1, dim_t NE2,
40 double l0, double l1, double l2,
41 bool periodic0, bool periodic1, bool periodic2,
42 int order, int reduced_order,
43 bool useElementsOnFace, bool useFullElementOrder,
44 bool useMacroElements, bool optimize,
45 escript::JMPI mpiInfo)
46 {
47 const int N_PER_E = 2;
48 const int DIM = 3;
49 dim_t Nstride0=0, Nstride1=0, Nstride2=0, local_NE0, local_NE1, local_NE2;
50 index_t e_offset0, e_offset1, e_offset2;
51
52 const int myRank = mpiInfo->rank;
53
54 // set up the global dimensions of the mesh
55 NE0 = std::max(dim_t(1), NE0);
56 NE1 = std::max(dim_t(1), NE1);
57 NE2 = std::max(dim_t(1), NE2);
58 const dim_t N0 = N_PER_E*NE0+1;
59 const dim_t N1 = N_PER_E*NE1+1;
60 const dim_t N2 = N_PER_E*NE2+1;
61
62 // allocate mesh
63 std::stringstream name;
64 name << "Brick " << N0 << " x " << N1 << " x " << N2;
65 FinleyDomain* out = new FinleyDomain(name.str(), DIM, mpiInfo);
66
67 const_ReferenceElementSet_ptr refPoints, refContactElements, refFaceElements, refElements;
68 bool generateAllNodes = useFullElementOrder || useMacroElements;
69
70 if (generateAllNodes) {
71 if (useMacroElements) {
72 refElements.reset(new ReferenceElementSet(Hex27Macro, order, reduced_order));
73 } else {
74 refElements.reset(new ReferenceElementSet(Hex27, order, reduced_order));
75 }
76 if (useElementsOnFace) {
77 throw escript::NotImplementedError("rich elements for Hex27 elements are not supported.");
78 } else {
79 if (useMacroElements) {
80 refFaceElements.reset(new ReferenceElementSet(Rec9Macro, order, reduced_order));
81 } else {
82 refFaceElements.reset(new ReferenceElementSet(Rec9, order, reduced_order));
83 }
84 refContactElements.reset(new ReferenceElementSet(Rec9_Contact, order, reduced_order));
85 }
86 } else { // !generateAllNodes
87 refElements.reset(new ReferenceElementSet(Hex20, order, reduced_order));
88 if (useElementsOnFace) {
89 refFaceElements.reset(new ReferenceElementSet(Hex20Face, order, reduced_order));
90 refContactElements.reset(new ReferenceElementSet(Hex20Face_Contact, order, reduced_order));
91 } else {
92 refFaceElements.reset(new ReferenceElementSet(Rec8, order, reduced_order));
93 refContactElements.reset(new ReferenceElementSet(Rec8_Contact, order, reduced_order));
94 }
95 }
96 refPoints.reset(new ReferenceElementSet(Point1, order, reduced_order));
97
98 ElementFile* points = new ElementFile(refPoints, mpiInfo);
99 out->setPoints(points);
100 ElementFile* contacts = new ElementFile(refContactElements, mpiInfo);
101 out->setContactElements(contacts);
102 ElementFile* faces = new ElementFile(refFaceElements, mpiInfo);
103 out->setFaceElements(faces);
104 ElementFile* elements = new ElementFile(refElements, mpiInfo);
105 out->setElements(elements);
106
107 // work out the largest dimension
108 if (N2==MAX3(N0,N1,N2)) {
109 Nstride0 = 1;
110 Nstride1 = N0;
111 Nstride2 = N0*N1;
112 local_NE0 = NE0;
113 e_offset0 = 0;
114 local_NE1 = NE1;
115 e_offset1 = 0;
116 mpiInfo->split(NE2, &local_NE2, &e_offset2);
117 } else if (N1==MAX3(N0,N1,N2)) {
118 Nstride0 = N2;
119 Nstride1 = N0*N2;
120 Nstride2 = 1;
121 local_NE0 = NE0;
122 e_offset0 = 0;
123 mpiInfo->split(NE1, &local_NE1, &e_offset1);
124 local_NE2 = NE2;
125 e_offset2 = 0;
126 } else {
127 Nstride0 = N1*N2;
128 Nstride1 = 1;
129 Nstride2 = N1;
130 mpiInfo->split(NE0, &local_NE0, &e_offset0);
131 local_NE1 = NE1;
132 e_offset1 = 0;
133 local_NE2 = NE2;
134 e_offset2 = 0;
135 }
136 const index_t offset0 = e_offset0*N_PER_E;
137 const index_t offset1 = e_offset1*N_PER_E;
138 const index_t offset2 = e_offset2*N_PER_E;
139 const dim_t local_N0 = local_NE0>0 ? local_NE0*N_PER_E+1 : 0;
140 const dim_t local_N1 = local_NE1>0 ? local_NE1*N_PER_E+1 : 0;
141 const dim_t local_N2 = local_NE2>0 ? local_NE2*N_PER_E+1 : 0;
142 dim_t NDOF0=0, NDOF1=0, NDOF2=0;
143
144 // get the number of surface elements
145 dim_t NFaceElements = 0;
146 if (!periodic2 && local_NE2 > 0) {
147 NDOF2=N2;
148 if (offset2==0) NFaceElements+=local_NE1*local_NE0;
149 if (local_NE2+e_offset2 == NE2) NFaceElements+=local_NE1*local_NE0;
150 } else {
151 NDOF2=N2-1;
152 }
153
154 if (!periodic0 && local_NE0 > 0) {
155 NDOF0=N0;
156 if (e_offset0 == 0) NFaceElements+=local_NE1*local_NE2;
157 if (local_NE0+e_offset0 == NE0) NFaceElements+=local_NE1*local_NE2;
158 } else {
159 NDOF0=N0-1;
160 }
161 if (!periodic1 && local_NE1 > 0) {
162 NDOF1=N1;
163 if (e_offset1 == 0) NFaceElements+=local_NE0*local_NE2;
164 if (local_NE1+e_offset1 == NE1) NFaceElements+=local_NE0*local_NE2;
165 } else {
166 NDOF1=N1-1;
167 }
168
169 // allocate tables
170 NodeFile* nodes = out->getNodes();
171 nodes->allocTable(local_N0*local_N1*local_N2);
172 elements->allocTable(local_NE0*local_NE1*local_NE2);
173 faces->allocTable(NFaceElements);
174
175 // create nodes
176 #pragma omp parallel for
177 for (index_t i2=0; i2<local_N2; i2++) {
178 for (index_t i1=0; i1<local_N1; i1++) {
179 for (index_t i0=0; i0<local_N0; i0++) {
180 const dim_t k = i0+local_N0*i1+local_N0*local_N1*i2;
181 const index_t global_i0 = i0+offset0;
182 const index_t global_i1 = i1+offset1;
183 const index_t global_i2 = i2+offset2;
184 nodes->Coordinates[INDEX2(0,k,DIM)] = (real_t)global_i0/(real_t)(N0-1)*l0;
185 nodes->Coordinates[INDEX2(1,k,DIM)] = (real_t)global_i1/(real_t)(N1-1)*l1;
186 nodes->Coordinates[INDEX2(2,k,DIM)] = (real_t)global_i2/(real_t)(N2-1)*l2;
187 nodes->Id[k] = Nstride0*global_i0+Nstride1*global_i1+Nstride2*global_i2;
188 nodes->Tag[k] = 0;
189 nodes->globalDegreesOfFreedom[k] = Nstride0*(global_i0%NDOF0)
190 + Nstride1*(global_i1%NDOF1)
191 + Nstride2*(global_i2%NDOF2);
192 }
193 }
194 }
195
196 // set the elements
197 dim_t NN = elements->numNodes;
198 index_t* eNodes = elements->Nodes;
199 #pragma omp parallel for
200 for (index_t i2=0; i2<local_NE2; i2++) {
201 for (index_t i1=0; i1<local_NE1; i1++) {
202 for (index_t i0=0; i0<local_NE0; i0++) {
203 const dim_t k = i0+local_NE0*i1+local_NE0*local_NE1*i2;
204 const index_t node0 = Nstride0*N_PER_E*(i0+e_offset0)
205 + Nstride1*N_PER_E*(i1+e_offset1)
206 + Nstride2*N_PER_E*(i2+e_offset2);
207
208 elements->Id[k] = (i0+e_offset0) + NE0*(i1+e_offset1)
209 + NE0*NE1*(i2+e_offset2);
210 elements->Tag[k] = 0;
211 elements->Owner[k] = myRank;
212
213 eNodes[INDEX2(0,k,NN)] =node0;
214 eNodes[INDEX2(1,k,NN)] =node0+ 2*Nstride0;
215 eNodes[INDEX2(2,k,NN)] =node0+ 2*Nstride1+2*Nstride0;
216 eNodes[INDEX2(3,k,NN)] =node0+ 2*Nstride1;
217 eNodes[INDEX2(4,k,NN)] =node0+2*Nstride2;
218 eNodes[INDEX2(5,k,NN)] =node0+2*Nstride2 +2*Nstride0;
219 eNodes[INDEX2(6,k,NN)] =node0+2*Nstride2+2*Nstride1+2*Nstride0;
220 eNodes[INDEX2(7,k,NN)] =node0+2*Nstride2+2*Nstride1;
221 eNodes[INDEX2(8,k,NN)] =node0+ 1*Nstride0;
222 eNodes[INDEX2(9,k,NN)] =node0+ 1*Nstride1+2*Nstride0;
223 eNodes[INDEX2(10,k,NN)]=node0+ 2*Nstride1+1*Nstride0;
224 eNodes[INDEX2(11,k,NN)]=node0+ 1*Nstride1;
225 eNodes[INDEX2(12,k,NN)]=node0+1*Nstride2;
226 eNodes[INDEX2(13,k,NN)]=node0+1*Nstride2 +2*Nstride0;
227 eNodes[INDEX2(14,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
228 eNodes[INDEX2(15,k,NN)]=node0+1*Nstride2+2*Nstride1;
229 eNodes[INDEX2(16,k,NN)]=node0+2*Nstride2 +1*Nstride0;
230 eNodes[INDEX2(17,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
231 eNodes[INDEX2(18,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
232 eNodes[INDEX2(19,k,NN)]=node0+2*Nstride2+1*Nstride1;
233 if (generateAllNodes) {
234 eNodes[INDEX2(20,k,NN)]=node0+ 1*Nstride1+1*Nstride0;
235 eNodes[INDEX2(21,k,NN)]=node0+1*Nstride2 +1*Nstride0;
236 eNodes[INDEX2(22,k,NN)]=node0+1*Nstride2+1*Nstride1+2*Nstride0;
237 eNodes[INDEX2(23,k,NN)]=node0+1*Nstride2+2*Nstride1+1*Nstride0;
238 eNodes[INDEX2(24,k,NN)]=node0+1*Nstride2+1*Nstride1;
239 eNodes[INDEX2(25,k,NN)]=node0+2*Nstride2+1*Nstride1+1*Nstride0;
240 eNodes[INDEX2(26,k,NN)]=node0+1*Nstride2+1*Nstride1+1*Nstride0;
241 }
242 }
243 }
244 }
245
246 // face elements
247 NN = faces->numNodes;
248 dim_t totalNECount=NE0*NE1*NE2;
249 dim_t faceNECount = 0;
250 eNodes = faces->Nodes;
251
252 // these are the quadrilateral elements on boundary 1 (x3=0):
253 if (!periodic2 && local_NE2 > 0) {
254 // ** elements on boundary 100 (x3=0):
255 if (offset2==0) {
256 #pragma omp parallel for
257 for (index_t i1=0; i1<local_NE1; i1++) {
258 for (index_t i0=0; i0<local_NE0; i0++) {
259 const dim_t k = i0+local_NE0*i1+faceNECount;
260 const index_t node0 = Nstride0*N_PER_E*(i0+e_offset0)
261 + Nstride1*N_PER_E*(i1+e_offset1);
262
263 faces->Id[k] = (i0+e_offset0) + NE0*(i1+e_offset1)
264 + totalNECount;
265 faces->Tag[k] = 100;
266 faces->Owner[k] = myRank;
267
268 if (useElementsOnFace) {
269 eNodes[INDEX2(0,k,NN)] =node0;
270 eNodes[INDEX2(1,k,NN)] =node0 +2*Nstride1;
271 eNodes[INDEX2(2,k,NN)] =node0 +2*Nstride1+2*Nstride0;
272 eNodes[INDEX2(3,k,NN)] =node0+ 2*Nstride0;
273 eNodes[INDEX2(4,k,NN)] =node0+2*Nstride2;
274 eNodes[INDEX2(5,k,NN)] =node0+2*Nstride2+2*Nstride1;
275 eNodes[INDEX2(6,k,NN)] =node0+2*Nstride2+2*Nstride1+2*Nstride0;
276 eNodes[INDEX2(7,k,NN)] =node0+2*Nstride2 +2*Nstride0;
277 eNodes[INDEX2(8,k,NN)] =node0+ 1*Nstride1;
278 eNodes[INDEX2(9,k,NN)] =node0+ 2*Nstride1+1*Nstride0;
279 eNodes[INDEX2(10,k,NN)]=node0+ 1*Nstride1+2*Nstride0;
280 eNodes[INDEX2(11,k,NN)]=node0+ 1*Nstride0;
281 eNodes[INDEX2(12,k,NN)]=node0+1*Nstride2;
282 eNodes[INDEX2(13,k,NN)]=node0+1*Nstride2+2*Nstride1;
283 eNodes[INDEX2(14,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
284 eNodes[INDEX2(15,k,NN)]=node0+1*Nstride2 +2*Nstride0;
285 eNodes[INDEX2(16,k,NN)]=node0+2*Nstride2+1*Nstride1;
286 eNodes[INDEX2(17,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
287 eNodes[INDEX2(18,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
288 eNodes[INDEX2(19,k,NN)]=node0+2*Nstride2 +1*Nstride0;
289 } else { // !useElementsOnFace
290 eNodes[INDEX2(0,k,NN)] =node0;
291 eNodes[INDEX2(1,k,NN)] =node0+ 2*Nstride1;
292 eNodes[INDEX2(2,k,NN)] =node0+ 2*Nstride1+2*Nstride0;
293 eNodes[INDEX2(3,k,NN)] =node0+ 2*Nstride0;
294 eNodes[INDEX2(4,k,NN)] =node0+ 1*Nstride1;
295 eNodes[INDEX2(5,k,NN)] =node0+ 2*Nstride1+1*Nstride0;
296 eNodes[INDEX2(6,k,NN)] =node0+ 1*Nstride1+2*Nstride0;
297 eNodes[INDEX2(7,k,NN)] =node0+ 1*Nstride0;
298 if (generateAllNodes) {
299 eNodes[INDEX2(8,k,NN)]=node0+ 1*Nstride1+1*Nstride0;
300 }
301 }
302 }
303 }
304 faceNECount+=local_NE1*local_NE0;
305 }
306 totalNECount+=NE1*NE0;
307
308 // ** elements on boundary 200 (x3=1):
309 if (local_NE2+e_offset2 == NE2) {
310 #pragma omp parallel for
311 for (index_t i1=0; i1<local_NE1; i1++) {
312 for (index_t i0=0; i0<local_NE0; i0++) {
313 const dim_t k = i0+local_NE0*i1+faceNECount;
314 const index_t node0 = Nstride0*N_PER_E*(i0+e_offset0)
315 + Nstride1*N_PER_E*(i1+e_offset1)
316 + Nstride2*N_PER_E*(NE2-1);
317
318 faces->Id[k] = (i0+e_offset0) + NE0*(i1+e_offset1)
319 + totalNECount;
320 faces->Tag[k] = 200;
321 faces->Owner[k] = myRank;
322 if (useElementsOnFace) {
323 eNodes[INDEX2(0,k,NN)] =node0+2*Nstride2;
324 eNodes[INDEX2(1,k,NN)] =node0+2*Nstride2+ 2*Nstride0;
325 eNodes[INDEX2(2,k,NN)] =node0+2*Nstride2+2*Nstride1+2*Nstride0;
326 eNodes[INDEX2(3,k,NN)] =node0+2*Nstride2+2*Nstride1;
327
328 eNodes[INDEX2(4,k,NN)] =node0;
329 eNodes[INDEX2(5,k,NN)] =node0+2*Nstride0;
330 eNodes[INDEX2(6,k,NN)] =node0+ 2*Nstride1+2*Nstride0;
331 eNodes[INDEX2(7,k,NN)] =node0+ 2*Nstride1;
332
333 eNodes[INDEX2(8,k,NN)] =node0+2*Nstride2+ 1*Nstride0;
334 eNodes[INDEX2(9,k,NN)] =node0+2*Nstride2+1*Nstride1+2*Nstride0;
335 eNodes[INDEX2(10,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
336 eNodes[INDEX2(11,k,NN)]=node0+2*Nstride2+1*Nstride1;
337
338 eNodes[INDEX2(12,k,NN)]=node0+1*Nstride2;
339 eNodes[INDEX2(13,k,NN)]=node0+1*Nstride2 +2*Nstride0;
340 eNodes[INDEX2(14,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
341 eNodes[INDEX2(15,k,NN)]=node0+1*Nstride2+2*Nstride1;
342
343 eNodes[INDEX2(16,k,NN)]=node0+ 1*Nstride0;
344 eNodes[INDEX2(17,k,NN)]=node0+ 1*Nstride1+2*Nstride0;
345 eNodes[INDEX2(18,k,NN)]=node0+ 2*Nstride1+1*Nstride0;
346 eNodes[INDEX2(19,k,NN)]=node0+ 1*Nstride1;
347 } else { // !useElementsOnFace
348 eNodes[INDEX2(0,k,NN)] =node0+2*Nstride2;
349 eNodes[INDEX2(1,k,NN)] =node0+2*Nstride2 +2*Nstride0;
350 eNodes[INDEX2(2,k,NN)] =node0+2*Nstride2+2*Nstride1+2*Nstride0;
351 eNodes[INDEX2(3,k,NN)] =node0+2*Nstride2+2*Nstride1;
352 eNodes[INDEX2(4,k,NN)] =node0+2*Nstride2 +1*Nstride0;
353 eNodes[INDEX2(5,k,NN)] =node0+2*Nstride2+1*Nstride1+2*Nstride0;
354 eNodes[INDEX2(6,k,NN)] =node0+2*Nstride2+2*Nstride1+1*Nstride0;
355 eNodes[INDEX2(7,k,NN)] =node0+2*Nstride2+1*Nstride1;
356 if (generateAllNodes) {
357 eNodes[INDEX2(8,k,NN)]=node0+2*Nstride2+1*Nstride1+1*Nstride0;
358 }
359 }
360 }
361 }
362 faceNECount+=local_NE1*local_NE0;
363 }
364 totalNECount+=NE1*NE0;
365 } // !periodic2 && local_NE2 > 0
366
367 if (!periodic0 && local_NE0 > 0) {
368 // ** elements on boundary 001 (x1=0):
369 if (e_offset0 == 0) {
370 #pragma omp parallel for
371 for (index_t i2=0; i2<local_NE2; i2++) {
372 for (index_t i1=0; i1<local_NE1; i1++) {
373 const dim_t k = i1+local_NE1*i2+faceNECount;
374 const index_t node0 = Nstride1*N_PER_E*(i1+e_offset1)
375 + Nstride2*N_PER_E*(i2+e_offset2);
376 faces->Id[k] = (i1+e_offset1) + NE1*(i2+e_offset2)
377 + totalNECount;
378 faces->Tag[k] = 1;
379 faces->Owner[k] = myRank;
380
381 if (useElementsOnFace) {
382 eNodes[INDEX2(0,k,NN)] =node0;
383 eNodes[INDEX2(1,k,NN)] =node0+2*Nstride2;
384 eNodes[INDEX2(2,k,NN)] =node0+2*Nstride2+2*Nstride1;
385 eNodes[INDEX2(3,k,NN)] =node0+2*Nstride1;
386
387 eNodes[INDEX2(4,k,NN)] =node0+2*Nstride0;
388 eNodes[INDEX2(5,k,NN)] =node0+2*Nstride2+2*Nstride0;
389 eNodes[INDEX2(6,k,NN)] =node0+2*Nstride2+2*Nstride1+2*Nstride0;
390 eNodes[INDEX2(7,k,NN)] =node0+2*Nstride1+2*Nstride0;
391
392 eNodes[INDEX2(8,k,NN)] =node0+1*Nstride2;
393 eNodes[INDEX2(9,k,NN)] =node0+2*Nstride2+1*Nstride1;
394 eNodes[INDEX2(10,k,NN)]=node0+1*Nstride2+2*Nstride1;
395 eNodes[INDEX2(11,k,NN)]=node0+ 1*Nstride1;
396
397 eNodes[INDEX2(12,k,NN)]=node0+ 1*Nstride0;
398 eNodes[INDEX2(13,k,NN)]=node0+2*Nstride2 +1*Nstride0;
399 eNodes[INDEX2(14,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
400 eNodes[INDEX2(15,k,NN)]=node0+2*Nstride1+ 1*Nstride0;
401
402 eNodes[INDEX2(16,k,NN)]=node0+1*Nstride2+ 2*Nstride0;
403 eNodes[INDEX2(17,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
404 eNodes[INDEX2(18,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
405 eNodes[INDEX2(19,k,NN)]=node0+1*Nstride1+ 2*Nstride0;
406 } else { // !useElementsOnFace
407 eNodes[INDEX2(0,k,NN)] =node0;
408 eNodes[INDEX2(1,k,NN)] =node0+2*Nstride2;
409 eNodes[INDEX2(2,k,NN)] =node0+2*Nstride2+2*Nstride1;
410 eNodes[INDEX2(3,k,NN)] =node0+ 2*Nstride1;
411 eNodes[INDEX2(4,k,NN)] =node0+1*Nstride2;
412 eNodes[INDEX2(5,k,NN)] =node0+2*Nstride2+1*Nstride1;
413 eNodes[INDEX2(6,k,NN)] =node0+1*Nstride2+2*Nstride1;
414 eNodes[INDEX2(7,k,NN)] =node0+ 1*Nstride1;
415 if (generateAllNodes) {
416 eNodes[INDEX2(8,k,NN)] =node0+1*Nstride2+1*Nstride1;
417 }
418 }
419 }
420 }
421 faceNECount+=local_NE1*local_NE2;
422 } // e_offset0 == 0
423 totalNECount+=NE1*NE2;
424
425 // ** elements on boundary 002 (x1=1):
426 if (local_NE0+e_offset0 == NE0) {
427 #pragma omp parallel for
428 for (index_t i2=0; i2<local_NE2; i2++) {
429 for (index_t i1=0; i1<local_NE1; i1++) {
430 const dim_t k = i1+local_NE1*i2+faceNECount;
431 const index_t node0 = Nstride0*N_PER_E*(NE0-1)
432 + Nstride1*N_PER_E*(i1+e_offset1)
433 + Nstride2*N_PER_E*(i2+e_offset2);
434 faces->Id[k] = (i1+e_offset1) + NE1*(i2+e_offset2)
435 + totalNECount;
436 faces->Tag[k] = 2;
437 faces->Owner[k] = myRank;
438
439 if (useElementsOnFace) {
440 eNodes[INDEX2(0,k,NN)]=node0+ 2*Nstride0;
441 eNodes[INDEX2(1,k,NN)]=node0+ 2*Nstride1+2*Nstride0;
442 eNodes[INDEX2(2,k,NN)]=node0+2*Nstride2+2*Nstride1+2*Nstride0;
443 eNodes[INDEX2(3,k,NN)]=node0+2*Nstride2+ 2*Nstride0;
444
445 eNodes[INDEX2(4,k,NN)]=node0;
446 eNodes[INDEX2(5,k,NN)]=node0+ 2*Nstride1;
447 eNodes[INDEX2(6,k,NN)]=node0+2*Nstride2+2*Nstride1;
448 eNodes[INDEX2(7,k,NN)]=node0+2*Nstride2;
449
450 eNodes[INDEX2(8,k,NN)]=node0+ 1*Nstride1+2*Nstride0;
451 eNodes[INDEX2(9,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
452 eNodes[INDEX2(10,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
453 eNodes[INDEX2(11,k,NN)]=node0+1*Nstride2+ 2*Nstride0;
454
455 eNodes[INDEX2(12,k,NN)]=node0+ 1*Nstride0;
456 eNodes[INDEX2(13,k,NN)]=node0+ 2*Nstride1+1*Nstride0;
457 eNodes[INDEX2(14,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
458 eNodes[INDEX2(15,k,NN)]=node0+2*Nstride2+ 1*Nstride0;
459
460 eNodes[INDEX2(16,k,NN)]=node0+ 1*Nstride1;
461 eNodes[INDEX2(17,k,NN)]=node0+1*Nstride2+2*Nstride1;
462 eNodes[INDEX2(18,k,NN)]=node0+2*Nstride2+1*Nstride1;
463 eNodes[INDEX2(19,k,NN)]=node0+1*Nstride2;
464 } else { // !useElementsOnFace
465 eNodes[INDEX2(0,k,NN)]=node0 +2*Nstride0;
466 eNodes[INDEX2(1,k,NN)]=node0+ 2*Nstride1+2*Nstride0;
467 eNodes[INDEX2(2,k,NN)]=node0+2*Nstride2+2*Nstride1+2*Nstride0;
468 eNodes[INDEX2(3,k,NN)]=node0+2*Nstride2+ 2*Nstride0;
469 eNodes[INDEX2(4,k,NN)]=node0+ 1*Nstride1+2*Nstride0;
470 eNodes[INDEX2(5,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
471 eNodes[INDEX2(6,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
472 eNodes[INDEX2(7,k,NN)]=node0+1*Nstride2 +2*Nstride0;
473 if (generateAllNodes) {
474 eNodes[INDEX2(8,k,NN)]=node0+1*Nstride2+1*Nstride1+2*Nstride0;
475 }
476 }
477 }
478 }
479 faceNECount+=local_NE1*local_NE2;
480 }
481 totalNECount+=NE1*NE2;
482 } // !periodic0 && local_NE0 > 0
483
484 if (!periodic1 && local_NE1 > 0) {
485 // ** elements on boundary 010 (x2=0):
486 if (e_offset1 == 0) {
487 #pragma omp parallel for
488 for (index_t i2=0; i2<local_NE2; i2++) {
489 for (index_t i0=0; i0<local_NE0; i0++) {
490 const dim_t k = i0+local_NE0*i2+faceNECount;
491 const index_t node0 = Nstride0*N_PER_E*(i0+e_offset0)
492 + Nstride2*N_PER_E*(i2+e_offset2);
493
494 faces->Id[k] = (i2+e_offset2) + NE2*(e_offset0+i0)
495 + totalNECount;
496 faces->Tag[k] = 10;
497 faces->Owner[k] = myRank;
498 if (useElementsOnFace) {
499 eNodes[INDEX2(0,k,NN)]=node0;
500 eNodes[INDEX2(1,k,NN)]=node0+ 2*Nstride0;
501 eNodes[INDEX2(2,k,NN)]=node0+2*Nstride2 +2*Nstride0;
502 eNodes[INDEX2(3,k,NN)]=node0+2*Nstride2;
503
504 eNodes[INDEX2(4,k,NN)]=node0+ 2*Nstride1;
505 eNodes[INDEX2(5,k,NN)]=node0+2*Nstride1+ 2*Nstride0;
506 eNodes[INDEX2(6,k,NN)]=node0+2*Nstride2+2*Nstride1+2*Nstride0;
507 eNodes[INDEX2(7,k,NN)]=node0+2*Nstride2+2*Nstride1;
508
509 eNodes[INDEX2(8,k,NN)]=node0+ 1*Nstride0;
510 eNodes[INDEX2(9,k,NN)]=node0+1*Nstride2+ 2*Nstride0;
511 eNodes[INDEX2(10,k,NN)]=node0+2*Nstride2+ 1*Nstride0;
512 eNodes[INDEX2(11,k,NN)]=node0+1*Nstride2;
513
514 eNodes[INDEX2(12,k,NN)]=node0+ 1*Nstride1;
515 eNodes[INDEX2(13,k,NN)]=node0+ 1*Nstride1+2*Nstride0;
516 eNodes[INDEX2(14,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
517 eNodes[INDEX2(15,k,NN)]=node0+2*Nstride2+1*Nstride1;
518
519 eNodes[INDEX2(16,k,NN)]=node0+ 2*Nstride1+1*Nstride0;
520 eNodes[INDEX2(17,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
521 eNodes[INDEX2(18,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
522 eNodes[INDEX2(19,k,NN)]=node0+1*Nstride2+2*Nstride1;
523 } else { // !useElementsOnFace
524 eNodes[INDEX2(0,k,NN)]=node0;
525 eNodes[INDEX2(1,k,NN)]=node0+ 2*Nstride0;
526 eNodes[INDEX2(2,k,NN)]=node0+2*Nstride2+ 2*Nstride0;
527 eNodes[INDEX2(3,k,NN)]=node0+2*Nstride2;
528 eNodes[INDEX2(4,k,NN)]=node0+ 1*Nstride0;
529 eNodes[INDEX2(5,k,NN)]=node0+1*Nstride2+ 2*Nstride0;
530 eNodes[INDEX2(6,k,NN)]=node0+2*Nstride2+ 1*Nstride0;
531 eNodes[INDEX2(7,k,NN)]=node0+1*Nstride2;
532 if (generateAllNodes) {
533 eNodes[INDEX2(8,k,NN)]=node0+1*Nstride2+ 1*Nstride0;
534 }
535 }
536 }
537 }
538 faceNECount+=local_NE0*local_NE2;
539 } // e_offset1==0
540 totalNECount+=NE0*NE2;
541
542 // ** elements on boundary 020 (x2=1):
543 if (local_NE1+e_offset1 == NE1) {
544 #pragma omp parallel for
545 for (index_t i2=0; i2<local_NE2; i2++) {
546 for (index_t i0=0; i0<local_NE0; i0++) {
547 const dim_t k = i0+local_NE0*i2+faceNECount;
548 const index_t node0 = Nstride0*N_PER_E*(i0+e_offset0)
549 + Nstride1*N_PER_E*(NE1-1)
550 + Nstride2*N_PER_E*(i2+e_offset2);
551
552 faces->Id[k] = (i2+e_offset2) + NE2*(i0+e_offset0)
553 + totalNECount;
554 faces->Tag[k] = 20;
555 faces->Owner[k] = myRank;
556
557 if (useElementsOnFace) {
558 eNodes[INDEX2(0,k,NN)]=node0+ 2*Nstride1;
559 eNodes[INDEX2(1,k,NN)]=node0+2*Nstride2+2*Nstride1;
560 eNodes[INDEX2(2,k,NN)]=node0+2*Nstride2+2*Nstride1+2*Nstride0;
561 eNodes[INDEX2(3,k,NN)]=node0+2*Nstride1+2*Nstride0;
562
563 eNodes[INDEX2(4,k,NN)]=node0;
564 eNodes[INDEX2(5,k,NN)]=node0+2*Nstride2;
565 eNodes[INDEX2(6,k,NN)]=node0+2*Nstride2+ 2*Nstride0;
566 eNodes[INDEX2(7,k,NN)]=node0+ 2*Nstride0;
567
568 eNodes[INDEX2(8,k,NN)]=node0+1*Nstride2+2*Nstride1;
569 eNodes[INDEX2(9,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
570 eNodes[INDEX2(10,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
571 eNodes[INDEX2(11,k,NN)]=node0+ 2*Nstride1+1*Nstride0;
572
573 eNodes[INDEX2(12,k,NN)]=node0+ 1*Nstride1;
574 eNodes[INDEX2(13,k,NN)]=node0+2*Nstride2+1*Nstride1;
575 eNodes[INDEX2(14,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
576 eNodes[INDEX2(15,k,NN)]=node0+ 1*Nstride1+2*Nstride0;
577
578 eNodes[INDEX2(16,k,NN)]=node0+1*Nstride2;
579 eNodes[INDEX2(17,k,NN)]=node0+2*Nstride2 +1*Nstride0;
580 eNodes[INDEX2(18,k,NN)]=node0+1*Nstride2 +2*Nstride0;
581 eNodes[INDEX2(19,k,NN)]=node0+ 1*Nstride0;
582 } else { // !useElementsOnFace
583 eNodes[INDEX2(0,k,NN)]=node0+ 2*Nstride1;
584 eNodes[INDEX2(1,k,NN)]=node0+2*Nstride2+2*Nstride1;
585 eNodes[INDEX2(2,k,NN)]=node0+2*Nstride2+2*Nstride1+2*Nstride0;
586 eNodes[INDEX2(3,k,NN)]=node0+ 2*Nstride1+2*Nstride0;
587 eNodes[INDEX2(4,k,NN)]=node0+1*Nstride2+2*Nstride1;
588 eNodes[INDEX2(5,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
589 eNodes[INDEX2(6,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
590 eNodes[INDEX2(7,k,NN)]=node0+ 2*Nstride1+1*Nstride0;
591 if (generateAllNodes) {
592 eNodes[INDEX2(8,k,NN)]=node0+1*Nstride2+2*Nstride1+1*Nstride0;
593 }
594 }
595 }
596 }
597 faceNECount+=local_NE0*local_NE2;
598 }
599 totalNECount+=NE0*NE2;
600 }
601
602 // add tag names
603 out->setTagMap("top", 200);
604 out->setTagMap("bottom", 100);
605 out->setTagMap("left", 1);
606 out->setTagMap("right", 2);
607 out->setTagMap("front", 10);
608 out->setTagMap("back", 20);
609
610 // prepare mesh for further calculations
611 out->resolveNodeIds();
612 out->prepare(optimize);
613 return out->getPtr();
614 }
615
616 } // namespace finley
617

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /branches/4.0fordebian/finley/src/Mesh_hex20.cpp:5567-5588 /branches/lapack2681/finley/src/Mesh_hex20.cpp:2682-2741 /branches/pasowrap/finley/src/Mesh_hex20.cpp:3661-3674 /branches/py3_attempt2/finley/src/Mesh_hex20.cpp:3871-3891 /branches/restext/finley/src/Mesh_hex20.cpp:2610-2624 /branches/ripleygmg_from_3668/finley/src/Mesh_hex20.cpp:3669-3791 /branches/stage3.0/finley/src/Mesh_hex20.cpp:2569-2590 /branches/symbolic_from_3470/finley/src/Mesh_hex20.cpp:3471-3974 /branches/trilinos_from_5897/finley/src/Mesh_hex20.cpp:5898-6118 /release/3.0/finley/src/Mesh_hex20.cpp:2591-2601 /release/4.0/finley/src/Mesh_hex20.cpp:5380-5406 /trunk/finley/src/Mesh_hex20.cpp:4257-4344

  ViewVC Help
Powered by ViewVC 1.1.26