/[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 4496 - (show annotations)
Mon Jul 15 06:53:44 2013 UTC (6 years ago) by caltinay
File size: 33645 byte(s)
finley (WIP):
-moved all of finley into its namespace
-introduced some shared pointers
-Mesh is now a class
-other bits and pieces...

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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /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 /release/3.0/finley/src/Mesh_hex20.cpp:2591-2601 /trunk/finley/src/Mesh_hex20.cpp:4257-4344

  ViewVC Help
Powered by ViewVC 1.1.26