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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26