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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4492 - (show annotations)
Tue Jul 2 01:44:11 2013 UTC (5 years, 9 months ago) by caltinay
File size: 14370 byte(s)
Finley changes that were held back while in release mode - moved more stuff
into finley namespace.

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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /branches/lapack2681/finley/src/Mesh_rec8.cpp:2682-2741 /branches/pasowrap/finley/src/Mesh_rec8.cpp:3661-3674 /branches/py3_attempt2/finley/src/Mesh_rec8.cpp:3871-3891 /branches/restext/finley/src/Mesh_rec8.cpp:2610-2624 /branches/ripleygmg_from_3668/finley/src/Mesh_rec8.cpp:3669-3791 /branches/stage3.0/finley/src/Mesh_rec8.cpp:2569-2590 /branches/symbolic_from_3470/finley/src/Mesh_rec8.cpp:3471-3974 /release/3.0/finley/src/Mesh_rec8.cpp:2591-2601 /trunk/finley/src/Mesh_rec8.cpp:4257-4344

  ViewVC Help
Powered by ViewVC 1.1.26