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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26