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

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;i1Nodes->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;i1Elements->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;qElements->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;i1FaceElements->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;i1FaceElements->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;i0FaceElements->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;i0FaceElements->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