/[escript]/branches/domexper/dudley/src/Mesh_tri3.c
ViewVC logotype

Contents of /branches/domexper/dudley/src/Mesh_tri3.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3220 - (show annotations)
Wed Sep 29 00:33:16 2010 UTC (8 years, 6 months ago) by jfenwick
File MIME type: text/plain
File size: 12754 byte(s)
Removing references to Quadrature.? and ShapeFns

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 /* Dudley: generates triangular meshes by splitting rectangles */
18
19 /* Generates a numElements[0] x numElements[1] x 2 mesh with first order elements (Tri3) in the rectangle */
20 /* [0,Length[0]] x [0,Length[1]]. order is the desired accuracy of the integration scheme. */
21
22
23 /**************************************************************/
24
25 #include "TriangularMesh.h"
26
27
28 Dudley_Mesh* Dudley_TriangularMesh_Tri3(dim_t* numElements,
29 double* Length,
30 index_t order,
31 index_t reduced_order,
32 bool_t optimize)
33 {
34 #define N_PER_E 1
35 #define DIM 2
36 dim_t N0,N1,NE0,NE1,i0,i1, Nstride0=0,Nstride1=0, local_NE0, local_NE1, local_N0=0, local_N1=0;
37 index_t offset0=0, offset1=0, e_offset0=0, e_offset1=0;
38 dim_t totalNECount,faceNECount,NDOF0=0,NDOF1=0,NFaceElements;
39 // Dudley_ReferenceElementSet *refPoints=NULL, *refFaceElements=NULL, *refElements=NULL;
40 index_t myRank;
41 Dudley_Mesh* out;
42 Paso_MPIInfo *mpi_info = NULL;
43 char name[50];
44 const int LEFTTAG=1; /* boundary x1=0 */
45 const int RIGHTTAG=2; /* boundary x1=1 */
46 const int BOTTOMTAG=10; /* boundary x2=0 */
47 const int TOPTAG=20; /* boundary x2=1 */
48
49 #ifdef Dudley_TRACE
50 double time0=Dudley_timer();
51 #endif
52
53 /* get MPI information */
54 mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
55 if (! Dudley_noError()) {
56 return NULL;
57 }
58 myRank=mpi_info->rank;
59
60 /* set up the global dimensions of the mesh */
61
62 NE0=MAX(1,numElements[0]);
63 NE1=MAX(1,numElements[1]);
64 N0=N_PER_E*NE0+1;
65 N1=N_PER_E*NE1+1;
66
67 /* This code was originally copied from Finley's Rec4 constructor.
68 NE? refers to the number of rectangular elements in each direction.
69 The number of nodes produced is the same but the number of non-face elements
70 will double.
71 */
72
73 /* allocate mesh: */
74 sprintf(name,"Triangular %d x %d (x 2) mesh",N0,N1);
75 out=Dudley_Mesh_alloc(name,DIM, mpi_info);
76 if (! Dudley_noError()) {
77 Paso_MPIInfo_free( mpi_info );
78 return NULL;
79 }
80 // refElements= Dudley_ReferenceElementSet_alloc(Tri3,order,reduced_order);
81 // refFaceElements=Dudley_ReferenceElementSet_alloc(Line2, order, reduced_order);
82 // refPoints=Dudley_ReferenceElementSet_alloc(Point1, order, reduced_order);
83
84 if ( Dudley_noError()) {
85
86 Dudley_Mesh_setPoints(out,Dudley_ElementFile_alloc(Point1, mpi_info));
87 Dudley_Mesh_setFaceElements(out,Dudley_ElementFile_alloc(Line2, mpi_info));
88 Dudley_Mesh_setElements(out,Dudley_ElementFile_alloc(Tri3, mpi_info));
89
90 //printf("N0=%d, N1=%d\n",N0,N1);
91
92 /* work out the largest dimension */
93 // if (N1==MAX(N0,N1)) {
94 // Nstride0=1;
95 // Nstride1=N0;
96 // local_NE0=NE0;
97 // e_offset0=0;
98 // Paso_MPIInfo_Split(mpi_info,NE1,&local_NE1,&e_offset1);
99 // printf("local_NE0=%d, local_NE1=%d, e_offset0=%d, e_offset1=%d\n", local_NE0, local_NE1, e_offset0, e_offset1);
100 //
101 // } else {
102 // Nstride0=N1;
103 // Nstride1=1;
104 // Paso_MPIInfo_Split(mpi_info,NE0,&local_NE0,&e_offset0);
105 // local_NE1=NE1;
106 // e_offset1=0;
107 //
108 // printf("local_NE0=%d, local_NE1=%d, e_offset0=%d, e_offset1=%d\n", local_NE0, local_NE1, e_offset0, e_offset1);
109 //
110 // }
111
112
113 // I'm trying this with always walking the same way
114 Nstride0=1;
115 Nstride1=N0;
116 if (N1==MAX(N0,N1)) {
117 local_NE0=NE0;
118 e_offset0=0;
119 Paso_MPIInfo_Split(mpi_info,NE1,&local_NE1,&e_offset1);
120 //printf("local_NE0=%d, local_NE1=%d, e_offset0=%d, e_offset1=%d\n", local_NE0, local_NE1, e_offset0, e_offset1);
121
122 } else {
123 Paso_MPIInfo_Split(mpi_info,NE0,&local_NE0,&e_offset0);
124 local_NE1=NE1;
125 e_offset1=0;
126
127 //printf("local_NE0=%d, local_NE1=%d, e_offset0=%d, e_offset1=%d\n", local_NE0, local_NE1, e_offset0, e_offset1);
128
129 }
130
131 //printf("Nstride0=%d, Nstride1=%d\n", Nstride0, Nstride1);
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 (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 (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
158 Dudley_NodeFile_allocTable(out->Nodes,local_N0*local_N1);
159 /* Dudley_ElementFile_allocTable(out->Elements,local_NE0*local_NE1);
160 Dudley_ElementFile_allocTable(out->FaceElements,NFaceElements);
161 */
162
163 /* This code was oringinally copied from Finley's rec4 generator
164 We double these numbers because each "rectangle" will be split into
165 two triangles. So the number of nodes is the same but the
166 number of elements will double */
167 Dudley_ElementFile_allocTable(out->Elements,local_NE0*local_NE1*2);
168 Dudley_ElementFile_allocTable(out->FaceElements,NFaceElements);
169
170
171
172 }
173 if (Dudley_noError()) {
174 /* create nodes */
175 #pragma omp parallel for private(i0,i1)
176 for (i1=0;i1<local_N1;i1++) {
177 for (i0=0;i0<local_N0;i0++) {
178 dim_t k=i0+local_N0*i1;
179 dim_t global_i0=i0+offset0;
180 dim_t global_i1=i1+offset1;
181 out->Nodes->Coordinates[INDEX2(0,k,DIM)]=DBLE(global_i0)/DBLE(N0-1)*Length[0];
182 out->Nodes->Coordinates[INDEX2(1,k,DIM)]=DBLE(global_i1)/DBLE(N1-1)*Length[1];
183 out->Nodes->Id[k]=Nstride0*global_i0+Nstride1*global_i1;
184 out->Nodes->Tag[k]=0;
185 out->Nodes->globalDegreesOfFreedom[k]=Nstride0*(global_i0%NDOF0)
186 +Nstride1*(global_i1%NDOF1);
187
188
189 // printf("N=%d: %f,%f\n", k, out->Nodes->Coordinates[INDEX2(0,k,DIM)],
190 // out->Nodes->Coordinates[INDEX2(1,k,DIM)]);
191
192 }
193 }
194 //printf("Now listing elements\n");
195 /* set the elements: */
196 dim_t NN=out->Elements->numNodes;
197 index_t global_adjustment=(offset0+offset1)%2;
198 #pragma omp parallel for private(i0,i1)
199 for (i1=0;i1<local_NE1;i1++)
200 {
201 for (i0=0;i0<local_NE0;i0++)
202 {
203 /* we will split this "rectangle" into two triangles */
204 dim_t k=2*(i0+local_NE0*i1);
205 index_t node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1);
206
207 out->Elements->Id[k]=2*((i0+e_offset0)+NE0*(i1+e_offset1));
208 out->Elements->Tag[k]=0;
209 out->Elements->Owner[k]=myRank;
210 out->Elements->Id[k+1]=out->Elements->Id[k]+1;
211 out->Elements->Tag[k+1]=0;
212 out->Elements->Owner[k+1]=myRank;
213
214 /* a,b,c,d gives the nodes in the rectangle in clockwise order*/
215 index_t a=node0,b=node0+Nstride0,c=node0+Nstride1+Nstride0,d=node0+Nstride1;
216
217 // printf("node0=%d, Nstride0=%d, Nstride1=%d\n", node0, Nstride0, Nstride1);
218 // printf("a=%d, b=%d, c=%d, d=%d\n",a,b,c,d);
219
220
221 /* For a little bit of variety */
222 if ((global_adjustment+node0)%2)
223 {
224 out->Elements->Nodes[INDEX2(0,k,NN)]=a;
225 out->Elements->Nodes[INDEX2(1,k,NN)]=b;
226 out->Elements->Nodes[INDEX2(2,k,NN)]=d;
227 out->Elements->Nodes[INDEX2(0,k+1,NN)]=b;
228 out->Elements->Nodes[INDEX2(1,k+1,NN)]=c;
229 out->Elements->Nodes[INDEX2(2,k+1,NN)]=d;
230 }
231 else
232 {
233 out->Elements->Nodes[INDEX2(0,k,NN)]=a;
234 out->Elements->Nodes[INDEX2(1,k,NN)]=b;
235 out->Elements->Nodes[INDEX2(2,k,NN)]=c;
236 out->Elements->Nodes[INDEX2(0,k+1,NN)]=a;
237 out->Elements->Nodes[INDEX2(1,k+1,NN)]=c;
238 out->Elements->Nodes[INDEX2(2,k+1,NN)]=d;
239 }
240
241 // for (int q=k;q<k+2;++q)
242 // {
243 // printf("E=%d: %d, %d, %d\n", q, out->Elements->Nodes[INDEX2(0,q,NN)],
244 // out->Elements->Nodes[INDEX2(1,q,NN)],
245 // out->Elements->Nodes[INDEX2(2,q,NN)]);
246 // }
247
248 }
249 }
250 //printf("Starting face elements\n");
251 /* face elements */
252 NN=out->FaceElements->numNodes;
253 totalNECount=2*NE0*NE1; /* because we have split the rectangles */
254 faceNECount=0;
255 if (local_NE0>0) {
256 /* ** elements on boundary 001 (x1=0): */
257
258 if (e_offset0 == 0) {
259 #pragma omp parallel for private(i1)
260 for (i1=0;i1<local_NE1;i1++) {
261
262 dim_t k=i1+faceNECount;
263 index_t node0=Nstride1*N_PER_E*(i1+e_offset1);
264
265 out->FaceElements->Id[k]=i1+e_offset1+totalNECount;
266 out->FaceElements->Tag[k]=LEFTTAG;
267 out->FaceElements->Owner[k]=myRank;
268 out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride1;
269 out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0;
270
271
272
273 // printf("E=%d: %d=%d %d=%d\n",k,INDEX2(0,k,NN),out->FaceElements->Nodes[INDEX2(0,k,NN)],
274 // INDEX2(1,k,NN),out->FaceElements->Nodes[INDEX2(1,k,NN)]);
275 }
276 faceNECount+=local_NE1;
277 }
278 totalNECount+=NE1;
279 /* ** elements on boundary 002 (x1=1): */
280 if (local_NE0+e_offset0 == NE0) {
281 #pragma omp parallel for private(i1)
282 for (i1=0;i1<local_NE1;i1++) {
283 dim_t k=i1+faceNECount;
284 index_t node0=Nstride0*N_PER_E*(NE0-1)+Nstride1*N_PER_E*(i1+e_offset1);
285
286 out->FaceElements->Id[k]=(i1+e_offset1)+totalNECount;
287 out->FaceElements->Tag[k]=RIGHTTAG;
288 out->FaceElements->Owner[k]=myRank;
289 out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride0;
290 out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride1+Nstride0;
291
292 // printf("E=%d: %d=%d %d=%d\n",k,INDEX2(0,k,NN),out->FaceElements->Nodes[INDEX2(0,k,NN)],
293 // INDEX2(1,k,NN),out->FaceElements->Nodes[INDEX2(1,k,NN)]);
294 }
295 faceNECount+=local_NE1;
296 }
297 totalNECount+=NE1;
298 }
299 if (local_NE1>0) {
300 /* ** elements on boundary 010 (x2=0): */
301 if (e_offset1 == 0) {
302 #pragma omp parallel for private(i0)
303 for (i0=0;i0<local_NE0;i0++) {
304 dim_t k=i0+faceNECount;
305 index_t node0=Nstride0*N_PER_E*(i0+e_offset0);
306
307 out->FaceElements->Id[k]=e_offset0+i0+totalNECount;
308 out->FaceElements->Tag[k]=BOTTOMTAG;
309 out->FaceElements->Owner[k]=myRank;
310
311 out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0;
312 out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride0;
313
314 // printf("E=%d: %d=%d %d=%d\n",k,INDEX2(0,k,NN),out->FaceElements->Nodes[INDEX2(0,k,NN)],
315 // INDEX2(1,k,NN),out->FaceElements->Nodes[INDEX2(1,k,NN)]);
316 }
317 faceNECount+=local_NE0;
318 }
319 totalNECount+=NE0;
320 /* ** elements on boundary 020 (x2=1): */
321 if (local_NE1+e_offset1 == NE1) {
322 #pragma omp parallel for private(i0)
323 for (i0=0;i0<local_NE0;i0++) {
324 dim_t k=i0+faceNECount;
325 index_t node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(NE1-1);
326
327 out->FaceElements->Id[k]=i0+e_offset0+totalNECount;
328 out->FaceElements->Tag[k]=TOPTAG;
329 out->FaceElements->Owner[k]=myRank;
330
331 out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride1+Nstride0;
332 out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride1;
333 /*printf("E=%d: %d=%d %d=%d\n",k,INDEX2(0,k,NN),out->FaceElements->Nodes[INDEX2(0,k,NN)],
334 INDEX2(1,k,NN),out->FaceElements->Nodes[INDEX2(1,k,NN)]); */
335 }
336 faceNECount+=local_NE0;
337 }
338 totalNECount+=NE0;
339 }
340 }
341 if (Dudley_noError()) {
342 /* add tag names */
343 Dudley_Mesh_addTagMap(out,"top", TOPTAG);
344 Dudley_Mesh_addTagMap(out,"bottom", BOTTOMTAG);
345 Dudley_Mesh_addTagMap(out,"left", LEFTTAG);
346 Dudley_Mesh_addTagMap(out,"right", RIGHTTAG);
347 }
348 /* prepare mesh for further calculatuions:*/
349 if (Dudley_noError()) {
350 Dudley_Mesh_resolveNodeIds(out);
351 }
352 if (Dudley_noError()) {
353 Dudley_Mesh_prepare(out, optimize);
354 }
355
356 /* free up memory */
357 // Dudley_ReferenceElementSet_dealloc(refPoints);
358 // Dudley_ReferenceElementSet_dealloc(refFaceElements);
359 // Dudley_ReferenceElementSet_dealloc(refElements);
360 Paso_MPIInfo_free( mpi_info );
361
362 return out;
363 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26