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

Contents of /trunk/finley/src/Mesh_hex8.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: 19295 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 numElements[0] x numElements[1] x numElements[2] mesh with
24 first order elements (Hex8) in the brick
25 [0,Length[0]] x [0,Length[1]] x [0,Length[2]].
26 order is the desired accuracy of the integration scheme.
27
28 *****************************************************************************/
29
30 #include "FinleyDomain.h"
31
32 #include <escript/index.h>
33
34 #define MAX3(_arg1_,_arg2_,_arg3_) std::max(_arg1_,std::max(_arg2_,_arg3_))
35
36 using escript::DataTypes::real_t;
37
38 namespace finley {
39
40 escript::Domain_ptr FinleyDomain::createHex8(dim_t NE0, dim_t NE1, dim_t NE2,
41 double l0, double l1, double l2,
42 bool periodic0, bool periodic1, bool periodic2,
43 int order, int reduced_order,
44 bool useElementsOnFace,
45 bool optimize, escript::JMPI mpiInfo)
46 {
47 const int N_PER_E = 1;
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 << "Rectangular " << N0 << " x " << N1 << " x " << N2 << " mesh";
65 FinleyDomain* out = new FinleyDomain(name.str(), DIM, mpiInfo);
66
67 const_ReferenceElementSet_ptr refPoints, refContactElements, refFaceElements, refElements;
68 if (useElementsOnFace) {
69 refFaceElements.reset(new ReferenceElementSet(Hex8Face, order, reduced_order));
70 refContactElements.reset(new ReferenceElementSet(Hex8Face_Contact, order, reduced_order));
71 } else {
72 refFaceElements.reset(new ReferenceElementSet(Rec4, order, reduced_order));
73 refContactElements.reset(new ReferenceElementSet(Rec4_Contact, order, reduced_order));
74 }
75 refElements.reset(new ReferenceElementSet(Hex8, order, reduced_order));
76 refPoints.reset(new ReferenceElementSet(Point1, order, reduced_order));
77
78 out->setPoints(new ElementFile(refPoints, mpiInfo));
79 out->setContactElements(new ElementFile(refContactElements, mpiInfo));
80 ElementFile* faces = new ElementFile(refFaceElements, mpiInfo);
81 out->setFaceElements(faces);
82 ElementFile* elements = new ElementFile(refElements, mpiInfo);
83 out->setElements(elements);
84
85 // work out the largest dimension
86 if (N2==MAX3(N0,N1,N2)) {
87 Nstride0 = 1;
88 Nstride1 = N0;
89 Nstride2 = N0*N1;
90 local_NE0 = NE0;
91 e_offset0 = 0;
92 local_NE1 = NE1;
93 e_offset1 = 0;
94 mpiInfo->split(NE2, &local_NE2, &e_offset2);
95 } else if (N1==MAX3(N0,N1,N2)) {
96 Nstride0 = N2;
97 Nstride1 = N0*N2;
98 Nstride2 = 1;
99 local_NE0 = NE0;
100 e_offset0 = 0;
101 mpiInfo->split(NE1, &local_NE1, &e_offset1);
102 local_NE2 = NE2;
103 e_offset2 = 0;
104 } else {
105 Nstride0 = N1*N2;
106 Nstride1 = 1;
107 Nstride2 = N1;
108 mpiInfo->split(NE0, &local_NE0, &e_offset0);
109 local_NE1 = NE1;
110 e_offset1 = 0;
111 local_NE2 = NE2;
112 e_offset2 = 0;
113 }
114 const index_t offset0 = e_offset0*N_PER_E;
115 const index_t offset1 = e_offset1*N_PER_E;
116 const index_t offset2 = e_offset2*N_PER_E;
117 const dim_t local_N0 = local_NE0>0 ? local_NE0*N_PER_E+1 : 0;
118 const dim_t local_N1 = local_NE1>0 ? local_NE1*N_PER_E+1 : 0;
119 const dim_t local_N2 = local_NE2>0 ? local_NE2*N_PER_E+1 : 0;
120 dim_t NDOF0=0, NDOF1=0, NDOF2=0;
121
122 // get the number of surface elements
123 dim_t NFaceElements = 0;
124 if (!periodic2 && local_NE2 > 0) {
125 NDOF2=N2;
126 if (offset2==0)
127 NFaceElements+=local_NE1*local_NE0;
128 if (local_NE2+e_offset2 == NE2)
129 NFaceElements+=local_NE1*local_NE0;
130 } else {
131 NDOF2=N2-1;
132 }
133
134 if (!periodic0 && local_NE0 > 0) {
135 NDOF0=N0;
136 if (e_offset0 == 0)
137 NFaceElements+=local_NE1*local_NE2;
138 if (local_NE0+e_offset0 == NE0)
139 NFaceElements+=local_NE1*local_NE2;
140 } else {
141 NDOF0=N0-1;
142 }
143 if (!periodic1 && local_NE1 > 0) {
144 NDOF1=N1;
145 if (e_offset1 == 0)
146 NFaceElements+=local_NE0*local_NE2;
147 if (local_NE1+e_offset1 == NE1)
148 NFaceElements+=local_NE0*local_NE2;
149 } else {
150 NDOF1=N1-1;
151 }
152
153 // allocate tables
154 NodeFile* nodes = out->getNodes();
155 nodes->allocTable(local_N0*local_N1*local_N2);
156 elements->allocTable(local_NE0*local_NE1*local_NE2);
157 faces->allocTable(NFaceElements);
158
159 // create nodes
160 #pragma omp parallel for
161 for (index_t i2=0; i2<local_N2; i2++) {
162 for (index_t i1=0; i1<local_N1; i1++) {
163 for (index_t i0=0; i0<local_N0; i0++) {
164 const dim_t k = i0+local_N0*i1+local_N0*local_N1*i2;
165 const index_t global_i0 = i0+offset0;
166 const index_t global_i1 = i1+offset1;
167 const index_t global_i2 = i2+offset2;
168 nodes->Coordinates[INDEX2(0,k,DIM)]=(real_t)global_i0/(real_t)(N0-1)*l0;
169 nodes->Coordinates[INDEX2(1,k,DIM)]=(real_t)global_i1/(real_t)(N1-1)*l1;
170 nodes->Coordinates[INDEX2(2,k,DIM)]=(real_t)global_i2/(real_t)(N2-1)*l2;
171 nodes->Id[k] = Nstride0*global_i0+Nstride1*global_i1+Nstride2*global_i2;
172 nodes->Tag[k] = 0;
173 nodes->globalDegreesOfFreedom[k] = Nstride0*(global_i0%NDOF0)
174 + Nstride1*(global_i1%NDOF1)
175 + Nstride2*(global_i2%NDOF2);
176 }
177 }
178 }
179
180 // set the elements
181 dim_t NN = elements->numNodes;
182 index_t* eNodes = elements->Nodes;
183 #pragma omp parallel for
184 for (index_t i2 = 0; i2 < local_NE2; i2++) {
185 for (index_t i1 = 0; i1 < local_NE1; i1++) {
186 for (index_t i0 = 0; i0 < local_NE0; i0++) {
187 const dim_t k = i0+local_NE0*i1+local_NE0*local_NE1*i2;
188 const index_t node0 = Nstride0*N_PER_E*(i0+e_offset0)
189 + Nstride1*N_PER_E*(i1+e_offset1)
190 + Nstride2*N_PER_E*(i2+e_offset2);
191
192 elements->Id[k] = (i0+e_offset0) + NE0*(i1+e_offset1)
193 + NE0*NE1*(i2+e_offset2);
194 elements->Tag[k] = 0;
195 elements->Owner[k] = myRank;
196
197 eNodes[INDEX2(0,k,NN)] = node0;
198 eNodes[INDEX2(1,k,NN)] = node0+Nstride0;
199 eNodes[INDEX2(2,k,NN)] = node0+Nstride1+Nstride0;
200 eNodes[INDEX2(3,k,NN)] = node0+Nstride1;
201 eNodes[INDEX2(4,k,NN)] = node0+Nstride2;
202 eNodes[INDEX2(5,k,NN)] = node0+Nstride2+Nstride0;
203 eNodes[INDEX2(6,k,NN)] = node0+Nstride2+Nstride1+Nstride0;
204 eNodes[INDEX2(7,k,NN)] = node0+Nstride2+Nstride1;
205 }
206 }
207 }
208
209 // face elements
210 NN = faces->numNodes;
211 dim_t totalNECount = NE0*NE1*NE2;
212 dim_t faceNECount = 0;
213 eNodes = faces->Nodes;
214
215 // these are the quadrilateral elements on boundary 1 (x3=0):
216 if (!periodic2 && local_NE2 > 0) {
217 // ** elements on boundary 100 (x3=0):
218 if (e_offset2 == 0) {
219 #pragma omp parallel for
220 for (index_t i1 = 0; i1 < local_NE1; i1++) {
221 for (index_t i0=0; i0 < local_NE0; i0++) {
222 const dim_t k = i0 + local_NE0*i1 + faceNECount;
223 const index_t node0 = Nstride0*N_PER_E*(i0+e_offset0)
224 + Nstride1*N_PER_E*(i1+e_offset1);
225
226 faces->Id[k] = (i0+e_offset0) + NE0*(i1+e_offset1)
227 + totalNECount;
228 faces->Tag[k] = 100;
229 faces->Owner[k] = myRank;
230
231 if (useElementsOnFace) {
232 eNodes[INDEX2(0,k,NN)] = node0;
233 eNodes[INDEX2(1,k,NN)] = node0+Nstride1;
234 eNodes[INDEX2(2,k,NN)] = node0+Nstride1+Nstride0;
235 eNodes[INDEX2(3,k,NN)] = node0+Nstride0;
236 eNodes[INDEX2(4,k,NN)] = node0+Nstride2;
237 eNodes[INDEX2(5,k,NN)] = node0+Nstride2+Nstride1;
238 eNodes[INDEX2(6,k,NN)] = node0+Nstride2+Nstride1+Nstride0;
239 eNodes[INDEX2(7,k,NN)] = node0+Nstride2+Nstride0;
240 } else {
241 eNodes[INDEX2(0,k,NN)] = node0;
242 eNodes[INDEX2(1,k,NN)] = node0+Nstride1;
243 eNodes[INDEX2(2,k,NN)] = node0+Nstride1+Nstride0;
244 eNodes[INDEX2(3,k,NN)] = node0+Nstride0;
245 }
246 }
247 }
248 faceNECount += local_NE1*local_NE0;
249 }
250 totalNECount += NE1*NE0;
251
252 // ** elements on boundary 200 (x3=1):
253 if (local_NE2+e_offset2 == NE2) {
254 #pragma omp parallel for
255 for (index_t i1 = 0; i1 < local_NE1; i1++) {
256 for (index_t i0 = 0; i0 < local_NE0; i0++) {
257 const dim_t k = i0 + local_NE0*i1 + faceNECount;
258 const index_t node0 = Nstride0*N_PER_E*(i0+e_offset0)
259 + Nstride1*N_PER_E*(i1+e_offset1)
260 + Nstride2*N_PER_E*(NE2-1);
261
262 faces->Id[k] = (i0+e_offset0) + NE0*(i1+e_offset1)
263 + totalNECount;
264 faces->Tag[k] = 200;
265 faces->Owner[k] = myRank;
266 if (useElementsOnFace) {
267 eNodes[INDEX2(0,k,NN)] = node0+Nstride2;
268 eNodes[INDEX2(1,k,NN)] = node0+Nstride2+ Nstride0;
269 eNodes[INDEX2(2,k,NN)] = node0+Nstride2+Nstride1+Nstride0;
270 eNodes[INDEX2(3,k,NN)] = node0+Nstride2+Nstride1;
271
272 eNodes[INDEX2(4,k,NN)] = node0;
273 eNodes[INDEX2(5,k,NN)] = node0+Nstride0;
274 eNodes[INDEX2(6,k,NN)] = node0+ Nstride1+Nstride0;
275 eNodes[INDEX2(7,k,NN)] = node0+ Nstride1;
276 } else {
277 eNodes[INDEX2(0,k,NN)] = node0+Nstride2;
278 eNodes[INDEX2(1,k,NN)] = node0+Nstride2 +Nstride0;
279 eNodes[INDEX2(2,k,NN)] = node0+Nstride2+Nstride1+Nstride0;
280 eNodes[INDEX2(3,k,NN)] = node0+Nstride2+Nstride1;
281 }
282 }
283 }
284 faceNECount += local_NE1*local_NE0;
285 }
286 totalNECount += NE1*NE0;
287 } // !periodic2 && local_NE2 > 0
288
289 if (!periodic0 && local_NE0 > 0) {
290 // ** elements on boundary 001 (x1=0):
291 if (e_offset0 == 0) {
292 #pragma omp parallel for
293 for (index_t i2 = 0; i2 < local_NE2; i2++) {
294 for (index_t i1 = 0; i1 < local_NE1; i1++) {
295 const dim_t k = i1 + local_NE1*i2 + faceNECount;
296 const index_t node0 = Nstride1*N_PER_E*(i1+e_offset1)
297 + Nstride2*N_PER_E*(i2+e_offset2);
298 faces->Id[k] = (i1+e_offset1) + NE1*(i2+e_offset2)
299 + totalNECount;
300 faces->Tag[k] = 1;
301 faces->Owner[k] = myRank;
302
303 if (useElementsOnFace) {
304 eNodes[INDEX2(0,k,NN)] = node0;
305 eNodes[INDEX2(1,k,NN)] = node0+Nstride2;
306 eNodes[INDEX2(2,k,NN)] = node0+Nstride2+Nstride1;
307 eNodes[INDEX2(3,k,NN)] = node0+Nstride1;
308 eNodes[INDEX2(4,k,NN)] = node0+Nstride0;
309 eNodes[INDEX2(5,k,NN)] = node0+Nstride2+Nstride0;
310 eNodes[INDEX2(6,k,NN)] = node0+Nstride2+Nstride1+Nstride0;
311 eNodes[INDEX2(7,k,NN)] = node0+Nstride1+Nstride0;
312 } else {
313 eNodes[INDEX2(0,k,NN)] = node0;
314 eNodes[INDEX2(1,k,NN)] = node0+Nstride2;
315 eNodes[INDEX2(2,k,NN)] = node0+Nstride2+Nstride1;
316 eNodes[INDEX2(3,k,NN)] = node0+Nstride1;
317 }
318 }
319 }
320 faceNECount += local_NE1*local_NE2;
321 }
322 totalNECount += NE1*NE2;
323
324 // ** elements on boundary 002 (x1=1):
325 if (local_NE0+e_offset0 == NE0) {
326 #pragma omp parallel for
327 for (index_t i2 = 0; i2 < local_NE2; i2++) {
328 for (index_t i1 = 0; i1 < local_NE1; i1++) {
329 const dim_t k = i1 + local_NE1*i2 + faceNECount;
330 const index_t node0 = Nstride0*N_PER_E*(NE0-1)
331 + Nstride1*N_PER_E*(i1+e_offset1)
332 + Nstride2*N_PER_E*(i2+e_offset2);
333 faces->Id[k] = (i1+e_offset1) + NE1*(i2+e_offset2)
334 + totalNECount;
335 faces->Tag[k] = 2;
336 faces->Owner[k] = myRank;
337
338 if (useElementsOnFace) {
339 eNodes[INDEX2(0,k,NN)] = node0+Nstride0;
340 eNodes[INDEX2(1,k,NN)] = node0+Nstride1+Nstride0;
341 eNodes[INDEX2(2,k,NN)] = node0+Nstride2+Nstride1+Nstride0;
342 eNodes[INDEX2(3,k,NN)] = node0+Nstride2+Nstride0;
343
344 eNodes[INDEX2(4,k,NN)] = node0;
345 eNodes[INDEX2(5,k,NN)] = node0+Nstride1;
346 eNodes[INDEX2(6,k,NN)] = node0+Nstride2+Nstride1;
347 eNodes[INDEX2(7,k,NN)] = node0+Nstride2;
348 } else {
349 eNodes[INDEX2(0,k,NN)] = node0+Nstride0;
350 eNodes[INDEX2(1,k,NN)] = node0+Nstride1+Nstride0;
351 eNodes[INDEX2(2,k,NN)] = node0+Nstride2+Nstride1+Nstride0;
352 eNodes[INDEX2(3,k,NN)] = node0+Nstride2+Nstride0;
353 }
354 }
355 }
356 faceNECount += local_NE1*local_NE2;
357 }
358 totalNECount += NE1*NE2;
359 } // !periodic0 && local_NE0 > 0
360
361 if (!periodic1 && local_NE1 > 0) {
362 // ** elements on boundary 010 (x2=0):
363 if (e_offset1 == 0) {
364 #pragma omp parallel for
365 for (index_t i2 = 0; i2 < local_NE2; i2++) {
366 for (index_t i0 = 0; i0 < local_NE0; i0++) {
367 const dim_t k = i0 + local_NE0*i2 + faceNECount;
368 const index_t node0 = Nstride0*N_PER_E*(i0+e_offset0)
369 + Nstride2*N_PER_E*(i2+e_offset2);
370
371 faces->Id[k] = (i2+e_offset2) + NE2*(e_offset0+i0)
372 + totalNECount;
373 faces->Tag[k] = 10;
374 faces->Owner[k] = myRank;
375 if (useElementsOnFace) {
376 eNodes[INDEX2(0,k,NN)] = node0;
377 eNodes[INDEX2(1,k,NN)] = node0+Nstride0;
378 eNodes[INDEX2(2,k,NN)] = node0+Nstride2+Nstride0;
379 eNodes[INDEX2(3,k,NN)] = node0+Nstride2;
380
381 eNodes[INDEX2(4,k,NN)] = node0+Nstride1;
382 eNodes[INDEX2(5,k,NN)] = node0+Nstride1+Nstride0;
383 eNodes[INDEX2(6,k,NN)] = node0+Nstride2+Nstride1+Nstride0;
384 eNodes[INDEX2(7,k,NN)] = node0+Nstride2+Nstride1;
385 } else {
386 eNodes[INDEX2(0,k,NN)] = node0;
387 eNodes[INDEX2(1,k,NN)] = node0+Nstride0;
388 eNodes[INDEX2(2,k,NN)] = node0+Nstride2+Nstride0;
389 eNodes[INDEX2(3,k,NN)] = node0+Nstride2;
390 }
391 }
392 }
393 faceNECount += local_NE0*local_NE2;
394 }
395 totalNECount += NE0*NE2;
396
397 // ** elements on boundary 020 (x2=1):
398 if (local_NE1+e_offset1 == NE1) {
399 #pragma omp parallel for
400 for (index_t i2 = 0; i2 < local_NE2; i2++) {
401 for (index_t i0 = 0; i0 < local_NE0; i0++) {
402 const dim_t k = i0 + local_NE0*i2 + faceNECount;
403 const index_t node0 = Nstride0*N_PER_E*(i0+e_offset0)
404 + Nstride1*N_PER_E*(NE1-1)
405 + Nstride2*N_PER_E*(i2+e_offset2);
406
407 faces->Id[k] = (i2+e_offset2) + NE2*(i0+e_offset0)
408 + totalNECount;
409 faces->Tag[k] = 20;
410 faces->Owner[k] = myRank;
411
412 if (useElementsOnFace) {
413 eNodes[INDEX2(0,k,NN)] = node0+Nstride1;
414 eNodes[INDEX2(1,k,NN)] = node0+Nstride2+Nstride1;
415 eNodes[INDEX2(2,k,NN)] = node0+Nstride2+Nstride1+Nstride0;
416 eNodes[INDEX2(3,k,NN)] = node0+Nstride1+Nstride0;
417
418 eNodes[INDEX2(4,k,NN)] = node0;
419 eNodes[INDEX2(5,k,NN)] = node0+Nstride2;
420 eNodes[INDEX2(6,k,NN)] = node0+Nstride2+Nstride0;
421 eNodes[INDEX2(7,k,NN)] = node0+Nstride0;
422 } else {
423 eNodes[INDEX2(0,k,NN)] = node0+Nstride1;
424 eNodes[INDEX2(1,k,NN)] = node0+Nstride2+Nstride1;
425 eNodes[INDEX2(2,k,NN)] = node0+Nstride2+Nstride1+Nstride0;
426 eNodes[INDEX2(3,k,NN)] = node0+Nstride1+Nstride0;
427 }
428 }
429 }
430 faceNECount += local_NE0*local_NE2;
431 }
432 totalNECount += NE0*NE2;
433 }
434
435 // add tag names
436 out->setTagMap("top", 200);
437 out->setTagMap("bottom", 100);
438 out->setTagMap("left", 1);
439 out->setTagMap("right", 2);
440 out->setTagMap("front", 10);
441 out->setTagMap("back", 20);
442
443 // prepare mesh for further calculations
444 out->resolveNodeIds();
445 out->prepare(optimize);
446 return out->getPtr();
447 }
448
449 } // namespace finley
450

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26