/[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 3259 - (show annotations)
Mon Oct 11 01:48:14 2010 UTC (8 years, 11 months ago) by jfenwick
File MIME type: text/plain
File size: 14480 byte(s)
Merging dudley and scons updates from branches

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26