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

Annotation of /trunk/finley/src/Mesh_rec4.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: 11958 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 first order elements (Rec4) in the rectangle */
21     /* [0,Length[0]] x [0,Length[1]]. order is the desired accuracy of the integration scheme. */
22    
23    
24     /**************************************************************/
25    
26     #include "RectangularMesh.h"
27    
28    
29 ksteube 1312 Finley_Mesh* Finley_RectangularMesh_Rec4(dim_t* numElements,
30     double* Length,
31     bool_t* periodic,
32     index_t order,
33     index_t reduced_order,
34     bool_t useElementsOnFace,
35     bool_t useFullElementOrder,
36     bool_t optimize)
37 bcumming 730 {
38 ksteube 1312 #define N_PER_E 1
39     #define DIM 2
40     dim_t N0,N1,NE0,NE1,i0,i1,k,Nstride0,Nstride1, local_NE0, local_NE1, local_N0, local_N1, global_i0, global_i1;
41     index_t offset0, offset1, e_offset0, e_offset1;
42     dim_t totalNECount,faceNECount,NDOF0,NDOF1,NFaceElements,NN;
43     index_t node0, myRank;
44     Finley_Mesh* out;
45     Paso_MPIInfo *mpi_info = NULL;
46     char name[50];
47     double time0=Finley_timer();
48 bcumming 730
49 ksteube 1312 /* get MPI information */
50     mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
51     if (! Finley_noError()) {
52     return NULL;
53 bcumming 730 }
54 ksteube 1312 myRank=mpi_info->rank;
55 bcumming 730
56 ksteube 1312 /* set up the global dimensions of the mesh */
57 bcumming 730
58 jgs 82 NE0=MAX(1,numElements[0]);
59     NE1=MAX(1,numElements[1]);
60 ksteube 1312 N0=N_PER_E*NE0+1;
61     N1=N_PER_E*NE1+1;
62 jgs 82
63 ksteube 1312 /* allocate mesh: */
64 jgs 82 sprintf(name,"Rectangular %d x %d mesh",N0,N1);
65 ksteube 1312 out=Finley_Mesh_alloc(name,DIM,order, reduced_order, mpi_info);
66     if (! Finley_noError()) {
67     Paso_MPIInfo_free( mpi_info );
68 bcumming 782 return NULL;
69     }
70 bcumming 787
71 ksteube 1312 Finley_Mesh_setElements(out,Finley_ElementFile_alloc(Rec4,out->order,out->reduced_order,mpi_info));
72 jgs 82 if (useElementsOnFace) {
73 ksteube 1312 Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Rec4Face,
74     out->order,
75     out->reduced_order,
76     mpi_info));
77     Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Rec4Face_Contact,
78     out->order,
79     out->reduced_order,
80     mpi_info));
81 jgs 82 } else {
82 ksteube 1312 Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Line2,
83     out->order,
84     out->reduced_order,
85     mpi_info));
86     Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Line2_Contact,
87     out->order,
88     out->reduced_order,
89     mpi_info));
90 jgs 82 }
91 ksteube 1312 Finley_Mesh_setPoints(out,Finley_ElementFile_alloc(Point1,
92     out->order,
93     out->reduced_order,
94     mpi_info));
95 jgs 150 if (! Finley_noError()) {
96 ksteube 1312 Paso_MPIInfo_free( mpi_info );
97     Finley_Mesh_free(out);
98 jgs 82 return NULL;
99     }
100 jgs 153
101 ksteube 1312 /* work out the largest dimension */
102     if (N1==MAX(N0,N1)) {
103     Nstride0=1;
104     Nstride1=N0;
105     local_NE0=NE0;
106     e_offset0=0;
107     Paso_MPIInfo_Split(mpi_info,NE1,&local_NE1,&e_offset1);
108 jgs 82 } else {
109 ksteube 1312 Nstride0=N1;
110     Nstride1=1;
111     Paso_MPIInfo_Split(mpi_info,NE0,&local_NE0,&e_offset0);
112     local_NE1=NE1;
113     e_offset1=0;
114 jgs 82 }
115 ksteube 1312 offset0=e_offset0*N_PER_E;
116     offset1=e_offset1*N_PER_E;
117 gross 1733 local_N0=local_NE0>0 ? local_NE0*N_PER_E+1 : 0;
118     local_N1=local_NE1>0 ? local_NE1*N_PER_E+1 : 0;
119 jgs 82
120 ksteube 1312 /* get the number of surface elements */
121 jgs 82
122 ksteube 1312 NFaceElements=0;
123 gross 1733 if (!periodic[0] && (local_NE0>0)) {
124 ksteube 1312 NDOF0=N0;
125     if (e_offset0 == 0) NFaceElements+=local_NE1;
126     if (local_NE0+e_offset0 == NE0) NFaceElements+=local_NE1;
127 gross 964 } else {
128 ksteube 1312 NDOF0=N0-1;
129 jgs 82 }
130 gross 1733 if (!periodic[1] && (local_NE1>0)) {
131 ksteube 1312 NDOF1=N1;
132     if (e_offset1 == 0) NFaceElements+=local_NE0;
133     if (local_NE1+e_offset1 == NE1) NFaceElements+=local_NE0;
134 bcumming 730 } else {
135     NDOF1=N1-1;
136     }
137    
138 ksteube 1312 /* allocate tables: */
139 bcumming 730
140 ksteube 1312 Finley_NodeFile_allocTable(out->Nodes,local_N0*local_N1);
141     Finley_ElementFile_allocTable(out->Elements,local_NE0*local_NE1);
142 bcumming 730 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
143    
144 ksteube 1312 if (Finley_noError()) {
145     /* create nodes */
146     #pragma omp parallel for private(i0,i1,k,global_i0,global_i1)
147     for (i1=0;i1<local_N1;i1++) {
148     for (i0=0;i0<local_N0;i0++) {
149     k=i0+local_N0*i1;
150     global_i0=i0+offset0;
151     global_i1=i1+offset1;
152     out->Nodes->Coordinates[INDEX2(0,k,DIM)]=DBLE(global_i0)/DBLE(N0-1)*Length[0];
153     out->Nodes->Coordinates[INDEX2(1,k,DIM)]=DBLE(global_i1)/DBLE(N1-1)*Length[1];
154     out->Nodes->Id[k]=Nstride0*global_i0+Nstride1*global_i1;
155     out->Nodes->Tag[k]=0;
156     out->Nodes->globalDegreesOfFreedom[k]=Nstride0*(global_i0%NDOF0)
157     +Nstride1*(global_i1%NDOF1);
158     }
159     }
160     /* set the elements: */
161     NN=out->Elements->numNodes;
162     #pragma omp parallel for private(i0,i1,k,node0)
163     for (i1=0;i1<local_NE1;i1++) {
164     for (i0=0;i0<local_NE0;i0++) {
165    
166     k=i0+local_NE0*i1;
167     node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1);
168 bcumming 787
169 ksteube 1312 out->Elements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1);
170     out->Elements->Tag[k]=0;
171     out->Elements->Owner[k]=myRank;
172 bcumming 782
173 ksteube 1312 out->Elements->Nodes[INDEX2(0,k,NN)]=node0;
174     out->Elements->Nodes[INDEX2(1,k,NN)]=node0+Nstride0;
175     out->Elements->Nodes[INDEX2(2,k,NN)]=node0+Nstride1+Nstride0;
176     out->Elements->Nodes[INDEX2(3,k,NN)]=node0+Nstride1;
177     }
178     }
179     /* face elements */
180     NN=out->FaceElements->numNodes;
181     totalNECount=NE0*NE1;
182     faceNECount=0;
183 gross 1733 if (!periodic[0] && (local_NE0>0)) {
184 ksteube 1312 /* ** elements on boundary 001 (x1=0): */
185    
186     if (e_offset0 == 0) {
187     #pragma omp parallel for private(i1,k,node0)
188     for (i1=0;i1<local_NE1;i1++) {
189    
190     k=i1+faceNECount;
191     node0=Nstride1*N_PER_E*(i1+e_offset1);
192 bcumming 730
193 ksteube 1312 out->FaceElements->Id[k]=i1+e_offset1+totalNECount;
194     out->FaceElements->Tag[k]=1;
195     out->FaceElements->Owner[k]=myRank;
196     if (useElementsOnFace) {
197     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride1;
198     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0;
199     out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride0;
200     out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride1+Nstride0;
201     } else {
202     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride1;
203     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0;
204     }
205     }
206     faceNECount+=local_NE1;
207     }
208     totalNECount+=NE1;
209     /* ** elements on boundary 002 (x1=1): */
210     if (local_NE0+e_offset0 == NE0) {
211     #pragma omp parallel for private(i1,k,node0)
212     for (i1=0;i1<local_NE1;i1++) {
213     k=i1+faceNECount;
214     node0=Nstride0*N_PER_E*(NE0-1)+Nstride1*N_PER_E*(i1+e_offset1);
215 bcumming 730
216 ksteube 1312 out->FaceElements->Id[k]=(i1+e_offset1)+totalNECount;
217     out->FaceElements->Tag[k]=2;
218     out->FaceElements->Owner[k]=myRank;
219 bcumming 730
220 ksteube 1312 if (useElementsOnFace) {
221     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride0;
222     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride1+Nstride0;
223     out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride1;
224     out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0;
225     } else {
226     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride0;
227     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride1+Nstride0;
228     }
229     }
230     faceNECount+=local_NE1;
231     }
232     totalNECount+=NE1;
233     }
234 gross 1733 if (!periodic[1] && (local_NE1>0)) {
235 ksteube 1312 /* ** elements on boundary 010 (x2=0): */
236     if (e_offset1 == 0) {
237     #pragma omp parallel for private(i0,k,node0)
238     for (i0=0;i0<local_NE0;i0++) {
239     k=i0+faceNECount;
240     node0=Nstride0*N_PER_E*(i0+e_offset0);
241    
242     out->FaceElements->Id[k]=e_offset0+i0+totalNECount;
243     out->FaceElements->Tag[k]=10;
244     out->FaceElements->Owner[k]=myRank;
245    
246     if (useElementsOnFace) {
247     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0;
248     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride0;
249     out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride1+Nstride0;
250     out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride1;
251     } else {
252     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0;
253     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride0;
254     }
255     }
256     faceNECount+=local_NE0;
257     }
258     totalNECount+=NE0;
259     /* ** elements on boundary 020 (x2=1): */
260     if (local_NE1+e_offset1 == NE1) {
261     #pragma omp parallel for private(i0,k,node0)
262     for (i0=0;i0<local_NE0;i0++) {
263     k=i0+faceNECount;
264     node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(NE1-1);
265    
266     out->FaceElements->Id[k]=i0+e_offset0+totalNECount;
267     out->FaceElements->Tag[k]=20;
268     out->FaceElements->Owner[k]=myRank;
269     if (useElementsOnFace) {
270     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride1+Nstride0;
271     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride1;
272     out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0;
273     out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride0;
274     } else {
275     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride1+Nstride0;
276     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride1;
277     }
278     }
279     faceNECount+=local_NE0;
280     }
281     totalNECount+=NE0;
282     }
283     /* add tag names */
284     Finley_Mesh_addTagMap(out,"top", 20);
285     Finley_Mesh_addTagMap(out,"bottom", 10);
286     Finley_Mesh_addTagMap(out,"left", 1);
287     Finley_Mesh_addTagMap(out,"right", 2);
288    
289     /* prepare mesh for further calculatuions:*/
290     if (Finley_noError()) {
291     Finley_Mesh_resolveNodeIds(out);
292     }
293     if (Finley_noError()) {
294     Finley_Mesh_prepare(out, optimize);
295     }
296 bcumming 730 }
297    
298 ksteube 1312 if (!Finley_noError()) {
299     Finley_Mesh_free(out);
300 bcumming 730 }
301 bcumming 782 /* free up memory */
302 ksteube 1312 Paso_MPIInfo_free( mpi_info );
303 bcumming 782 #ifdef Finley_TRACE
304     printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
305     #endif
306 ksteube 1312
307     return out;
308 bcumming 730 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26