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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1811 - (show annotations)
Thu Sep 25 23:11:13 2008 UTC (11 years, 2 months ago) by ksteube
File MIME type: text/plain
File size: 34858 byte(s)
Copyright updated in all files

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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26