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

Annotation of /trunk/finley/src/Mesh_rec8.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1733 - (hide annotations)
Thu Aug 28 01:46:07 2008 UTC (11 years, 1 month ago) by gross
File MIME type: text/plain
File size: 15083 byte(s)
bug fix: if the numer of elements to be generated suppose to be zero
still face elements and nodes are generated. this is fixed now.


1 jgs 150
2 ksteube 1312 /* $Id$ */
3    
4     /*******************************************************
5     *
6     * Copyright 2003-2007 by ACceSS MNRF
7     * Copyright 2007 by University of Queensland
8     *
9     * http://esscc.uq.edu.au
10     * Primary Business: Queensland, Australia
11     * Licensed under the Open Software License version 3.0
12     * http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15    
16 jgs 82 /**************************************************************/
17    
18     /* Finley: generates rectangular meshes */
19    
20     /* Generates a numElements[0] x numElements[1] mesh with second order elements (Rec8) in the rectangle */
21     /* [0,Length[0]] x [0,Length[1]]. order is the desired accuracy of the integration scheme. */
22    
23     /**************************************************************/
24    
25     #include "RectangularMesh.h"
26    
27 bcumming 782
28 ksteube 1312 Finley_Mesh* Finley_RectangularMesh_Rec8(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 bcumming 782 {
37 ksteube 1312 #define N_PER_E 2
38     #define DIM 2
39     dim_t N0,N1,NE0,NE1,i0,i1,k,Nstride0,Nstride1;
40     dim_t totalNECount,faceNECount,NDOF0,NDOF1,NFaceElements,NN, local_NE0, local_NE1, local_N0, local_N1;
41     index_t e_offset1, e_offset0, offset0, offset1, global_i0, global_i1;
42     index_t node0, myRank;
43     Finley_Mesh* out;
44     Paso_MPIInfo *mpi_info = NULL;
45     char name[50];
46     double time0=Finley_timer();
47 bcumming 782
48 ksteube 1312 /* get MPI information */
49     mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
50     if (! Finley_noError()) {
51     return NULL;
52 bcumming 782 }
53 ksteube 1312 myRank=mpi_info->rank;
54 bcumming 782
55 ksteube 1312 /* set up the global dimensions of the mesh */
56 bcumming 782
57 jgs 82 NE0=MAX(1,numElements[0]);
58     NE1=MAX(1,numElements[1]);
59 ksteube 1312 N0=N_PER_E*NE0+1;
60     N1=N_PER_E*NE1+1;
61 jgs 82
62 ksteube 1312 /* allocate mesh: */
63 jgs 82 sprintf(name,"Rectangular %d x %d mesh",N0,N1);
64 ksteube 1312 out=Finley_Mesh_alloc(name,DIM,order, reduced_order, mpi_info);
65     if (! Finley_noError()) {
66     Paso_MPIInfo_free( mpi_info );
67 jgs 82 return NULL;
68     }
69 bcumming 782
70 ksteube 1312 if (useFullElementOrder) {
71     Finley_Mesh_setElements(out,Finley_ElementFile_alloc(Rec9,
72     out->order,
73     out->reduced_order,
74     mpi_info));
75     if (useElementsOnFace) {
76     Finley_setError(SYSTEM_ERROR,"rich elements for Rec9 elements is not supported yet.");
77     } else {
78     Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Line3,
79     out->order,
80     out->reduced_order,
81     mpi_info));
82     Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Line3_Contact,
83     out->order,
84     out->reduced_order,
85     mpi_info));
86 jgs 82 }
87 ksteube 1312 } else {
88     Finley_Mesh_setElements(out,Finley_ElementFile_alloc(Rec8,out->order,out->reduced_order,mpi_info));
89     if (useElementsOnFace) {
90     Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Rec8Face,
91     out->order,
92     out->reduced_order,
93     mpi_info));
94     Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Rec8Face_Contact,
95     out->order,
96     out->reduced_order,
97     mpi_info));
98     } else {
99     Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Line3,
100     out->order,
101     out->reduced_order,
102     mpi_info));
103     Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Line3_Contact,
104     out->order,
105     out->reduced_order,
106     mpi_info));
107 jgs 82 }
108     }
109 ksteube 1312 Finley_Mesh_setPoints(out,Finley_ElementFile_alloc(Point1,
110     out->order,
111     out->reduced_order,
112     mpi_info));
113 bcumming 787 if (! Finley_noError()) {
114 ksteube 1312 Paso_MPIInfo_free( mpi_info );
115     Finley_Mesh_free(out);
116 bcumming 787 return NULL;
117     }
118 jgs 82
119 ksteube 1312 /* work out the largest dimension */
120     if (N1==MAX(N0,N1)) {
121     Nstride0=1;
122     Nstride1=N0;
123     local_NE0=NE0;
124     e_offset0=0;
125     Paso_MPIInfo_Split(mpi_info,NE1,&local_NE1,&e_offset1);
126 gross 964 } else {
127 ksteube 1312 Nstride0=N1;
128     Nstride1=1;
129     Paso_MPIInfo_Split(mpi_info,NE0,&local_NE0,&e_offset0);
130     local_NE1=NE1;
131     e_offset1=0;
132 jgs 82 }
133 ksteube 1312 offset0=e_offset0*N_PER_E;
134     offset1=e_offset1*N_PER_E;
135 gross 1733 local_N0=local_NE0>0 ? local_NE0*N_PER_E+1 : 0;
136     local_N1=local_NE1>0 ? local_NE1*N_PER_E+1 : 0;
137 bcumming 782
138 ksteube 1312 /* get the number of surface elements */
139 bcumming 782
140     NFaceElements=0;
141 gross 1733 if (!periodic[0] && (local_NE0>0)) {
142 ksteube 1312 NDOF0=N0;
143     if (e_offset0 == 0) NFaceElements+=local_NE1;
144     if (local_NE0+e_offset0 == NE0) NFaceElements+=local_NE1;
145 bcumming 782 } else {
146     NDOF0=N0-1;
147     }
148 gross 1733 if (!periodic[1] && (local_NE1>0)) {
149 ksteube 1312 NDOF1=N1;
150     if (e_offset1 == 0) NFaceElements+=local_NE0;
151     if (local_NE1+e_offset1 == NE1) NFaceElements+=local_NE0;
152 bcumming 782 } else {
153     NDOF1=N1-1;
154     }
155    
156 ksteube 1312 /* allocate tables: */
157 bcumming 782
158 ksteube 1312 Finley_NodeFile_allocTable(out->Nodes,local_N0*local_N1);
159     Finley_ElementFile_allocTable(out->Elements,local_NE0*local_NE1);
160 bcumming 782 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
161    
162 ksteube 1312 if (Finley_noError()) {
163     /* create nodes */
164    
165     #pragma omp parallel for private(i0,i1,k,global_i0,global_i1)
166     for (i1=0;i1<local_N1;i1++) {
167     for (i0=0;i0<local_N0;i0++) {
168     k=i0+local_N0*i1;
169     global_i0=i0+offset0;
170     global_i1=i1+offset1;
171     out->Nodes->Coordinates[INDEX2(0,k,DIM)]=DBLE(global_i0)/DBLE(N0-1)*Length[0];
172     out->Nodes->Coordinates[INDEX2(1,k,DIM)]=DBLE(global_i1)/DBLE(N1-1)*Length[1];
173     out->Nodes->Id[k]=Nstride0*global_i0+Nstride1*global_i1;
174     out->Nodes->Tag[k]=0;
175     out->Nodes->globalDegreesOfFreedom[k]=Nstride0*(global_i0%NDOF0)
176     +Nstride1*(global_i1%NDOF1);
177     }
178     }
179     /* set the elements: */
180     NN=out->Elements->numNodes;
181     #pragma omp parallel for private(i0,i1,k,node0)
182     for (i1=0;i1<local_NE1;i1++) {
183     for (i0=0;i0<local_NE0;i0++) {
184    
185     k=i0+local_NE0*i1;
186     node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1);
187 bcumming 782
188 ksteube 1312 out->Elements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1);
189     out->Elements->Tag[k]=0;
190     out->Elements->Owner[k]=myRank;
191 bcumming 782
192 ksteube 1312 out->Elements->Nodes[INDEX2(0,k,NN)]=node0;
193     out->Elements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride0;
194     out->Elements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride1+2*Nstride0;
195     out->Elements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride1;
196     out->Elements->Nodes[INDEX2(4,k,NN)]=node0+1*Nstride0;
197     out->Elements->Nodes[INDEX2(5,k,NN)]=node0+Nstride1+2*Nstride0;
198     out->Elements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride1+1*Nstride0;
199     out->Elements->Nodes[INDEX2(7,k,NN)]=node0+Nstride1;
200     if (useFullElementOrder) {
201     out->Elements->Nodes[INDEX2(8,k,NN)]=node0+1*Nstride1+1*Nstride0;
202     }
203     }
204 bcumming 782 }
205 ksteube 1312 /* face elements */
206     NN=out->FaceElements->numNodes;
207     totalNECount=NE0*NE1;
208     faceNECount=0;
209 gross 1733 if (!periodic[0] && (local_NE0>0)) {
210 ksteube 1312 /* ** elements on boundary 001 (x1=0): */
211    
212     if (e_offset0 == 0) {
213     #pragma omp parallel for private(i1,k,node0)
214     for (i1=0;i1<local_NE1;i1++) {
215    
216     k=i1+faceNECount;
217     node0=Nstride1*N_PER_E*(i1+e_offset1);
218 bcumming 782
219 ksteube 1312 out->FaceElements->Id[k]=i1+e_offset1+totalNECount;
220     out->FaceElements->Tag[k]=1;
221     out->FaceElements->Owner[k]=myRank;
222     if (useElementsOnFace) {
223     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride1;
224     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0;
225     out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride0;
226     out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride1+2*Nstride0;
227     out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+Nstride1;
228     out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+1*Nstride0;
229     out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+Nstride1+2*Nstride0;
230     out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+2*Nstride1+1*Nstride0;
231     } else {
232     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride1;
233     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0;
234     out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride1;
235     }
236     }
237     faceNECount+=local_NE1;
238     }
239     totalNECount+=NE1;
240     /* ** elements on boundary 002 (x1=1): */
241     if (local_NE0+e_offset0 == NE0) {
242     #pragma omp parallel for private(i1,k,node0)
243     for (i1=0;i1<local_NE1;i1++) {
244     k=i1+faceNECount;
245     node0=Nstride0*N_PER_E*(NE0-1)+Nstride1*N_PER_E*(i1+e_offset1);
246 bcumming 782
247 ksteube 1312 out->FaceElements->Id[k]=(i1+e_offset1)+totalNECount;
248     out->FaceElements->Tag[k]=2;
249     out->FaceElements->Owner[k]=myRank;
250 bcumming 782
251 ksteube 1312 if (useElementsOnFace) {
252     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride0;
253     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride1+2*Nstride0;
254     out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride1;
255     out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0;
256     out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+Nstride1+2*Nstride0;
257     out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+2*Nstride1+1*Nstride0;
258     out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+Nstride1;
259     out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+1*Nstride0;
260     } else {
261     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride0;
262     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride1+2*Nstride0;
263     out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride1+2*Nstride0;
264     }
265     }
266     faceNECount+=local_NE1;
267     }
268     totalNECount+=NE1;
269 bcumming 782 }
270 gross 1733 if (!periodic[1] && (local_NE0>0)) {
271 ksteube 1312 /* ** elements on boundary 010 (x2=0): */
272     if (e_offset1 == 0) {
273     #pragma omp parallel for private(i0,k,node0)
274     for (i0=0;i0<local_NE0;i0++) {
275     k=i0+faceNECount;
276     node0=Nstride0*N_PER_E*(i0+e_offset0);
277    
278     out->FaceElements->Id[k]=e_offset0+i0+totalNECount;
279     out->FaceElements->Tag[k]=10;
280     out->FaceElements->Owner[k]=myRank;
281    
282     if (useElementsOnFace) {
283     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0;
284     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride0;
285     out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride1+2*Nstride0;
286     out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride1;
287     out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+1*Nstride0;
288     out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+Nstride1+2*Nstride0;
289     out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride1+1*Nstride0;
290     out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+Nstride1;
291     } else {
292     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0;
293     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride0;
294     out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+1*Nstride0;
295     }
296     }
297     faceNECount+=local_NE0;
298     }
299     totalNECount+=NE0;
300     /* ** elements on boundary 020 (x2=1): */
301     if (local_NE1+e_offset1 == NE1) {
302     #pragma omp parallel for private(i0,k,node0)
303     for (i0=0;i0<local_NE0;i0++) {
304     k=i0+faceNECount;
305     node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(NE1-1);
306 bcumming 782
307 ksteube 1312 out->FaceElements->Id[k]=i0+e_offset0+totalNECount;
308     out->FaceElements->Tag[k]=20;
309     out->FaceElements->Owner[k]=myRank;
310     if (useElementsOnFace) {
311     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride1+2*Nstride0;
312     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride1;
313     out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0;
314     out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride0;
315     out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+2*Nstride1+1*Nstride0;
316     out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+Nstride1;
317     out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+1*Nstride0;
318     out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+Nstride1+2*Nstride0;
319     } else {
320     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride1+2*Nstride0;
321     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride1;
322     out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride1+1*Nstride0;
323     }
324     }
325     faceNECount+=local_NE0;
326     }
327     totalNECount+=NE0;
328     }
329     /* add tag names */
330     Finley_Mesh_addTagMap(out,"top", 20);
331     Finley_Mesh_addTagMap(out,"bottom", 10);
332     Finley_Mesh_addTagMap(out,"left", 1);
333     Finley_Mesh_addTagMap(out,"right", 2);
334 bcumming 782
335 ksteube 1312 /* prepare mesh for further calculatuions:*/
336     if (Finley_noError()) {
337     Finley_Mesh_resolveNodeIds(out);
338 bcumming 782 }
339 ksteube 1312 if (Finley_noError()) {
340     Finley_Mesh_prepare(out, optimize);
341     }
342 bcumming 782 }
343 ksteube 1312 if (!Finley_noError()) {
344     Finley_Mesh_free(out);
345 bcumming 782 }
346 ksteube 1312 /* free up memory */
347     Paso_MPIInfo_free( mpi_info );
348 bcumming 782 #ifdef Finley_TRACE
349     printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
350     #endif
351    
352 gross 964 return out;
353 bcumming 782 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26