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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26