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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26