/[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 3114 - (show annotations)
Fri Aug 27 05:26:25 2010 UTC (9 years, 5 months ago) by jfenwick
File MIME type: text/plain
File size: 12663 byte(s)
It doesn't pass all tests but this is major progress

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(refPoints, mpi_info));
87 Dudley_Mesh_setFaceElements(out,Dudley_ElementFile_alloc(refFaceElements, mpi_info));
88 Dudley_Mesh_setElements(out,Dudley_ElementFile_alloc(refElements, 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 #pragma omp parallel for private(i0,i1)
198 for (i1=0;i1<local_NE1;i1++)
199 {
200 for (i0=0;i0<local_NE0;i0++)
201 {
202 /* we will split this "rectangle" into two triangles */
203 dim_t k=2*(i0+local_NE0*i1);
204 index_t node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1);
205
206 out->Elements->Id[k]=2*((i0+e_offset0)+NE0*(i1+e_offset1));
207 out->Elements->Tag[k]=0;
208 out->Elements->Owner[k]=myRank;
209 out->Elements->Id[k+1]=out->Elements->Id[k]+1;
210 out->Elements->Tag[k+1]=0;
211 out->Elements->Owner[k+1]=myRank;
212
213 /* a,b,c,d gives the nodes in the rectangle in clockwise order*/
214 index_t a=node0,b=node0+Nstride0,c=node0+Nstride1+Nstride0,d=node0+Nstride1;
215
216 //printf("node0=%d, Nstride0=%d, Nstride1=%d\n", node0, Nstride0, Nstride1);
217 //printf("a=%d, b=%d, c=%d, d=%d\n",a,b,c,d);
218
219
220 /* For a little bit of variety */
221 if (node0%2)
222 {
223 out->Elements->Nodes[INDEX2(0,k,NN)]=a;
224 out->Elements->Nodes[INDEX2(1,k,NN)]=b;
225 out->Elements->Nodes[INDEX2(2,k,NN)]=d;
226 out->Elements->Nodes[INDEX2(0,k+1,NN)]=b;
227 out->Elements->Nodes[INDEX2(1,k+1,NN)]=c;
228 out->Elements->Nodes[INDEX2(2,k+1,NN)]=d;
229 }
230 else
231 {
232 out->Elements->Nodes[INDEX2(0,k,NN)]=a;
233 out->Elements->Nodes[INDEX2(1,k,NN)]=b;
234 out->Elements->Nodes[INDEX2(2,k,NN)]=c;
235 out->Elements->Nodes[INDEX2(0,k+1,NN)]=a;
236 out->Elements->Nodes[INDEX2(1,k+1,NN)]=c;
237 out->Elements->Nodes[INDEX2(2,k+1,NN)]=d;
238 }
239 /*
240 for (int q=k;q<k+2;++q)
241 {
242 printf("E=%d: %d, %d, %d\n", q, out->Elements->Nodes[INDEX2(0,q,NN)],
243 out->Elements->Nodes[INDEX2(1,q,NN)],
244 out->Elements->Nodes[INDEX2(2,q,NN)]);
245 }
246 */
247 }
248 }
249 //printf("Starting face elements\n");
250 /* face elements */
251 NN=out->FaceElements->numNodes;
252 totalNECount=2*NE0*NE1; /* because we have split the rectangles */
253 faceNECount=0;
254 if (local_NE0>0) {
255 /* ** elements on boundary 001 (x1=0): */
256
257 if (e_offset0 == 0) {
258 #pragma omp parallel for private(i1)
259 for (i1=0;i1<local_NE1;i1++) {
260
261 dim_t k=i1+faceNECount;
262 index_t node0=Nstride1*N_PER_E*(i1+e_offset1);
263
264 out->FaceElements->Id[k]=i1+e_offset1+totalNECount;
265 out->FaceElements->Tag[k]=LEFTTAG;
266 out->FaceElements->Owner[k]=myRank;
267 out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride1;
268 out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0;
269
270
271
272 //printf("E=%d: %d=%d %d=%d\n",k,INDEX2(0,k,NN),out->FaceElements->Nodes[INDEX2(0,k,NN)],
273 //INDEX2(1,k,NN),out->FaceElements->Nodes[INDEX2(1,k,NN)]);
274 }
275 faceNECount+=local_NE1;
276 }
277 totalNECount+=NE1;
278 /* ** elements on boundary 002 (x1=1): */
279 if (local_NE0+e_offset0 == NE0) {
280 #pragma omp parallel for private(i1)
281 for (i1=0;i1<local_NE1;i1++) {
282 dim_t k=i1+faceNECount;
283 index_t node0=Nstride0*N_PER_E*(NE0-1)+Nstride1*N_PER_E*(i1+e_offset1);
284
285 out->FaceElements->Id[k]=(i1+e_offset1)+totalNECount;
286 out->FaceElements->Tag[k]=RIGHTTAG;
287 out->FaceElements->Owner[k]=myRank;
288 out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride0;
289 out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride1+Nstride0;
290
291 //printf("E=%d: %d=%d %d=%d\n",k,INDEX2(0,k,NN),out->FaceElements->Nodes[INDEX2(0,k,NN)],
292 //INDEX2(1,k,NN),out->FaceElements->Nodes[INDEX2(1,k,NN)]);
293 }
294 faceNECount+=local_NE1;
295 }
296 totalNECount+=NE1;
297 }
298 if (local_NE1>0) {
299 /* ** elements on boundary 010 (x2=0): */
300 if (e_offset1 == 0) {
301 #pragma omp parallel for private(i0)
302 for (i0=0;i0<local_NE0;i0++) {
303 dim_t k=i0+faceNECount;
304 index_t node0=Nstride0*N_PER_E*(i0+e_offset0);
305
306 out->FaceElements->Id[k]=e_offset0+i0+totalNECount;
307 out->FaceElements->Tag[k]=BOTTOMTAG;
308 out->FaceElements->Owner[k]=myRank;
309
310 out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0;
311 out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride0;
312
313 //printf("E=%d: %d=%d %d=%d\n",k,INDEX2(0,k,NN),out->FaceElements->Nodes[INDEX2(0,k,NN)],
314 //INDEX2(1,k,NN),out->FaceElements->Nodes[INDEX2(1,k,NN)]);
315 }
316 faceNECount+=local_NE0;
317 }
318 totalNECount+=NE0;
319 /* ** elements on boundary 020 (x2=1): */
320 if (local_NE1+e_offset1 == NE1) {
321 #pragma omp parallel for private(i0)
322 for (i0=0;i0<local_NE0;i0++) {
323 dim_t k=i0+faceNECount;
324 index_t node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(NE1-1);
325
326 out->FaceElements->Id[k]=i0+e_offset0+totalNECount;
327 out->FaceElements->Tag[k]=TOPTAG;
328 out->FaceElements->Owner[k]=myRank;
329 out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride1+Nstride0;
330 out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride1;
331 //printf("E=%d: %d=%d %d=%d\n",k,INDEX2(0,k,NN),out->FaceElements->Nodes[INDEX2(0,k,NN)],
332 //INDEX2(1,k,NN),out->FaceElements->Nodes[INDEX2(1,k,NN)]);
333 }
334 faceNECount+=local_NE0;
335 }
336 totalNECount+=NE0;
337 }
338 }
339 if (Dudley_noError()) {
340 /* add tag names */
341 Dudley_Mesh_addTagMap(out,"top", TOPTAG);
342 Dudley_Mesh_addTagMap(out,"bottom", BOTTOMTAG);
343 Dudley_Mesh_addTagMap(out,"left", LEFTTAG);
344 Dudley_Mesh_addTagMap(out,"right", RIGHTTAG);
345 }
346 /* prepare mesh for further calculatuions:*/
347 if (Dudley_noError()) {
348 Dudley_Mesh_resolveNodeIds(out);
349 }
350 if (Dudley_noError()) {
351 Dudley_Mesh_prepare(out, optimize);
352 }
353
354 /* free up memory */
355 Dudley_ReferenceElementSet_dealloc(refPoints);
356 Dudley_ReferenceElementSet_dealloc(refFaceElements);
357 Dudley_ReferenceElementSet_dealloc(refElements);
358 Paso_MPIInfo_free( mpi_info );
359
360 return out;
361 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26