/[escript]/branches/arrayview_from_1695_trunk/finley/src/Mesh_rec8.c
ViewVC logotype

Contents of /branches/arrayview_from_1695_trunk/finley/src/Mesh_rec8.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1781 - (show annotations)
Thu Sep 11 05:03:14 2008 UTC (10 years, 10 months ago) by jfenwick
File MIME type: text/plain
File size: 15083 byte(s)
Branch commit

Merged changes from trunk version 1695 up to and including version 1779.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26