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

Diff of /branches/domexper/dudley/src/Mesh_tet4.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3220 by jfenwick, Wed Sep 29 00:33:16 2010 UTC revision 3221 by jfenwick, Wed Sep 29 01:00:21 2010 UTC
# Line 37  Dudley_Mesh* Dudley_TriangularMesh_Tet4( Line 37  Dudley_Mesh* Dudley_TriangularMesh_Tet4(
37    dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,Nstride0=0, Nstride1=0,Nstride2=0, local_NE0, local_NE1, local_NE2, local_N0=0, local_N1=0, local_N2=0;    dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,Nstride0=0, Nstride1=0,Nstride2=0, local_NE0, local_NE1, local_NE2, local_N0=0, local_N1=0, local_N2=0;
38    dim_t totalNECount,faceNECount,NDOF0=0,NDOF1=0,NDOF2=0,NFaceElements=0, NN;    dim_t totalNECount,faceNECount,NDOF0=0,NDOF1=0,NDOF2=0,NFaceElements=0, NN;
39    index_t node0, myRank, e_offset2, e_offset1, e_offset0=0, offset1=0, offset2=0, offset0=0, global_i0, global_i1, global_i2;    index_t node0, myRank, e_offset2, e_offset1, e_offset0=0, offset1=0, offset2=0, offset0=0, global_i0, global_i1, global_i2;
 //  Dudley_ReferenceElementSet *refPoints=NULL, *refFaceElements=NULL, *refElements=NULL;  
40    Dudley_Mesh* out;    Dudley_Mesh* out;
41    Paso_MPIInfo *mpi_info = NULL;    Paso_MPIInfo *mpi_info = NULL;
42    char name[50];    char name[50];
# Line 78  Dudley_Mesh* Dudley_TriangularMesh_Tet4( Line 77  Dudley_Mesh* Dudley_TriangularMesh_Tet4(
77        Paso_MPIInfo_free( mpi_info );        Paso_MPIInfo_free( mpi_info );
78        return NULL;        return NULL;
79    }    }
 //  refElements= Dudley_ReferenceElementSet_alloc(Tet4,order,reduced_order);  
 //  refFaceElements=Dudley_ReferenceElementSet_alloc(Tri3, order, reduced_order);  
 //  refPoints=Dudley_ReferenceElementSet_alloc(Point1, order, reduced_order);  
     
   
80    if ( Dudley_noError()) {    if ( Dudley_noError()) {
81        
82        Dudley_Mesh_setPoints(out,Dudley_ElementFile_alloc(Point1, mpi_info));        Dudley_Mesh_setPoints(out,Dudley_ElementFile_alloc(Point1, mpi_info));
# Line 183  Dudley_Mesh* Dudley_TriangularMesh_Tet4( Line 177  Dudley_Mesh* Dudley_TriangularMesh_Tet4(
177             out->Nodes->globalDegreesOfFreedom[k]=Nstride0*(global_i0%NDOF0)             out->Nodes->globalDegreesOfFreedom[k]=Nstride0*(global_i0%NDOF0)
178                                                 +Nstride1*(global_i1%NDOF1)                                                 +Nstride1*(global_i1%NDOF1)
179                                                 +Nstride2*(global_i2%NDOF2);                                                 +Nstride2*(global_i2%NDOF2);
 //     printf("N=%d: %f,%f,%f\n", k, out->Nodes->Coordinates[INDEX2(0,k,DIM)],  
 //      out->Nodes->Coordinates[INDEX2(1,k,DIM)],  
 //      out->Nodes->Coordinates[INDEX2(2,k,DIM)]);  
180           }           }
181         }         }
182       }       }
 // printf("Now listing elements\n");  
183       /*   set the elements: */       /*   set the elements: */
184    
185       int global_adjustment=(offset0+offset1+offset2)%2; // If we are not the only rank we may need to shift our pattern to match neighbours       int global_adjustment=(offset0+offset1+offset2)%2; // If we are not the only rank we may need to shift our pattern to match neighbours
# Line 216  Dudley_Mesh* Dudley_TriangularMesh_Tet4( Line 206  Dudley_Mesh* Dudley_TriangularMesh_Tet4(
206         The bottom face (anticlockwise= 0,1,3,2), top face (anticlockwise 4,5,7,6)*/         The bottom face (anticlockwise= 0,1,3,2), top face (anticlockwise 4,5,7,6)*/
207         if ((global_adjustment+i0+i1+i2)%2==0)         if ((global_adjustment+i0+i1+i2)%2==0)
208         {         {
 // printf("Type0 %d, %d, %d\n",i0,i1,i2);  
209         v[0]=node0;         v[0]=node0;
210         v[1]=node0+Nstride0;         v[1]=node0+Nstride0;
211         v[2]=node0+Nstride1;         v[2]=node0+Nstride1;
# Line 228  Dudley_Mesh* Dudley_TriangularMesh_Tet4( Line 217  Dudley_Mesh* Dudley_TriangularMesh_Tet4(
217         }         }
218         else         else
219         {         {
 // printf("Type1 %d, %d, %d\n",i0,i1,i2);  
220         // this form is rotated around the 0,2,4,6 face clockwise 90 degrees         // this form is rotated around the 0,2,4,6 face clockwise 90 degrees
221    
222         v[0]=node0+Nstride1;         // node 0 ends up in position 2         v[0]=node0+Nstride1;         // node 0 ends up in position 2
# Line 242  Dudley_Mesh* Dudley_TriangularMesh_Tet4( Line 230  Dudley_Mesh* Dudley_TriangularMesh_Tet4(
230         v[5]=node0+Nstride0;         // node 5 -> pos 1         v[5]=node0+Nstride0;         // node 5 -> pos 1
231         }         }
232    
 // for (int z=0;z<8;++z)  
 // {  
 //     printf("z[%d]=%d\n", z, v[z]);  
 //  
 // }  
   
   
 //     index_t a=node0, b=node0+Nstride0, c=node0+Nstride1+Nstride0, d=node0+Nstride1;  
 //     index_t e=node0+Nstride2, f=node0+Nstride2+Nstride0, g=node0+Nstride2+Nstride1+Nstride0,  
 //       h=node0+Nstride2+Nstride1;  
 //      
 //  
 //  
 //     a=0, b=1, c=3, d=2  
 //     e=4, f=5, g=7, h=6  
   
233         // elements nodes are numbered: centre, x, y, z         // elements nodes are numbered: centre, x, y, z
234    
235         out->Elements->Nodes[INDEX2(0,k,NN)] =v[4];         out->Elements->Nodes[INDEX2(0,k,NN)] =v[4];
# Line 286  Dudley_Mesh* Dudley_TriangularMesh_Tet4( Line 258  Dudley_Mesh* Dudley_TriangularMesh_Tet4(
258         out->Elements->Nodes[INDEX2(1,k+4,NN)] =v[0];         out->Elements->Nodes[INDEX2(1,k+4,NN)] =v[0];
259         out->Elements->Nodes[INDEX2(2,k+4,NN)] =v[6];         out->Elements->Nodes[INDEX2(2,k+4,NN)] =v[6];
260         out->Elements->Nodes[INDEX2(3,k+4,NN)] =v[3];         out->Elements->Nodes[INDEX2(3,k+4,NN)] =v[3];
   
   
 // for (int z=0;z<5;++z)  
 // {  
 // printf("E %d:",z);  
 //    for (int q=0;q<4;++q)  
 //    {  
 //  index_t id=out->Elements->Nodes[INDEX2(q,z,NN)];  
 //  printf("   %d = %f, %f, %f\n", id, out->Nodes->Coordinates[INDEX2(0,id,DIM)],  
 //      out->Nodes->Coordinates[INDEX2(1,id,DIM)], out->Nodes->Coordinates[INDEX2(2,id,DIM)]);  
 //    }  
 // printf("\n");}  
   
   
 /*  
 for (int j=0;j<4;++j)  
 {  
   
   printf("Elt %d",j);  
   for (int m=0;m<4;++m)  
   {  
      printf(" %d",out->Elements->Nodes[INDEX2(m,k+j,NN)]);  
   }  
   printf("\n");  
 }*/  
   
   
   
261           }           }
262    
263         }         }
264       } /* end for */       } /* end for */
265       /* face elements */       /* face elements */
 // printf("Starting face elements\n");  
266       NN=out->FaceElements->numNodes;       NN=out->FaceElements->numNodes;
267       totalNECount=5*NE0*NE1*NE2;       totalNECount=5*NE0*NE1*NE2;
268       faceNECount=0;       faceNECount=0;
 //printf("Bottom\n");  
269       /*   these are the quadrilateral elements on boundary 1 (x3=0): */       /*   these are the quadrilateral elements on boundary 1 (x3=0): */
270       if (local_NE2>0)       if (local_NE2>0)
271       {       {
# Line 374  for (int j=0;j<4;++j) Line 316  for (int j=0;j<4;++j)
316    
317    
318            }            }
 // printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)],  out->FaceElements->Nodes[INDEX2(1,k,NN)],  
 // out->FaceElements->Nodes[INDEX2(2,k,NN)]);  
 // printf("%d: %d, %d, %d\n", k+1,out->FaceElements->Nodes[INDEX2(0,k+1,NN)],  out->FaceElements->Nodes[INDEX2(1,k+1,NN)],  
 // out->FaceElements->Nodes[INDEX2(2,k+1,NN)]);  
   
319              }              }
320            }            }
 // printf("\n");  
321            faceNECount+=2*local_NE1*local_NE0;            faceNECount+=2*local_NE1*local_NE0;
322         }         }
323         totalNECount+=2*NE1*NE0;         totalNECount+=2*NE1*NE0;
324           /* **  elements on boundary 200 (x3=1) - Top*/
 //printf("Top\n");  
        /* **  elements on boundary 200 (x3=1) */  
325         if (local_NE2+e_offset2 == NE2) {         if (local_NE2+e_offset2 == NE2) {
326            #pragma omp parallel for private(i0,i1,k,node0)            #pragma omp parallel for private(i0,i1,k,node0)
327            for (i1=0;i1<local_NE1;i1++) {            for (i1=0;i1<local_NE1;i1++) {
# Line 430  for (int j=0;j<4;++j) Line 364  for (int j=0;j<4;++j)
364              out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n7;              out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n7;
365              out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n6;              out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n6;
366            }            }
   
 // printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)],  out->FaceElements->Nodes[INDEX2(1,k,NN)],  
 // out->FaceElements->Nodes[INDEX2(2,k,NN)]);  
 // printf("%d: %d, %d, %d\n", k+1,out->FaceElements->Nodes[INDEX2(0,k+1,NN)],  out->FaceElements->Nodes[INDEX2(1,k+1,NN)],  
 // out->FaceElements->Nodes[INDEX2(2,k+1,NN)]);  
   
367              }              }
368            }            }
 /*printf("\n");*/  
369            faceNECount+=2*local_NE1*local_NE0;            faceNECount+=2*local_NE1*local_NE0;
370         }         }
371         totalNECount+=2*NE1*NE0;         totalNECount+=2*NE1*NE0;
372       }       }
 //printf("Left\n");  
373       if (local_NE0>0) {       if (local_NE0>0) {
374          /* **  elements on boundary 001 (x1=0): */          /* **  elements on boundary 001 (x1=0): - Left*/
375            
376          if (e_offset0 == 0) {          if (e_offset0 == 0) {
377             #pragma omp parallel for private(i1,i2,k,node0)             #pragma omp parallel for private(i1,i2,k,node0)
# Line 453  for (int j=0;j<4;++j) Line 379  for (int j=0;j<4;++j)
379               for (i1=0;i1<local_NE1;i1++) {               for (i1=0;i1<local_NE1;i1++) {
380                
381                 k=2*(i1+local_NE1*i2)+faceNECount;                 k=2*(i1+local_NE1*i2)+faceNECount;
   
 //printf("%d, %d ",k,k+1);  
382                 node0=Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);                 node0=Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
383             index_t res=2*((i1+e_offset1)+NE1*(i2+e_offset2))+totalNECount;             index_t res=2*((i1+e_offset1)+NE1*(i2+e_offset2))+totalNECount;
384                 out->FaceElements->Id[k]=res;                 out->FaceElements->Id[k]=res;
# Line 491  for (int j=0;j<4;++j) Line 415  for (int j=0;j<4;++j)
415              out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n6;              out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n6;
416              out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n2;              out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n2;
417         }         }
 // printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)],  out->FaceElements->Nodes[INDEX2(1,k,NN)],  
 // out->FaceElements->Nodes[INDEX2(2,k,NN)]);  
 // printf("%d: %d, %d, %d\n", k+1,out->FaceElements->Nodes[INDEX2(0,k+1,NN)],  out->FaceElements->Nodes[INDEX2(1,k+1,NN)],  
 // out->FaceElements->Nodes[INDEX2(2,k+1,NN)]);  
418               }               }
419             }             }
420             faceNECount+=2*local_NE1*local_NE2;             faceNECount+=2*local_NE1*local_NE2;
421          }          }
 /*printf("\n");*/  
 //printf("Right\n");  
422          totalNECount+=2*NE1*NE2;          totalNECount+=2*NE1*NE2;
423          /* **  elements on boundary 002 (x1=1): */          /* **  elements on boundary 002 (x1=1): - Right*/
424          if (local_NE0+e_offset0 == NE0) {          if (local_NE0+e_offset0 == NE0) {
425             #pragma omp parallel for private(i1,i2,k,node0)             #pragma omp parallel for private(i1,i2,k,node0)
426             for (i2=0;i2<local_NE2;i2++) {             for (i2=0;i2<local_NE2;i2++) {
427               for (i1=0;i1<local_NE1;i1++) {               for (i1=0;i1<local_NE1;i1++) {
428                 k=2*(i1+local_NE1*i2)+faceNECount;                 k=2*(i1+local_NE1*i2)+faceNECount;
429    
 //printf("%d, %d ",k,k+1);  
430                 node0=Nstride0*N_PER_E*(NE0-1)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);                 node0=Nstride0*N_PER_E*(NE0-1)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
431             index_t res=2*((i1+e_offset1)+NE1*(i2+e_offset2))+totalNECount;             index_t res=2*((i1+e_offset1)+NE1*(i2+e_offset2))+totalNECount;
432                 out->FaceElements->Id[k]=res;                 out->FaceElements->Id[k]=res;
# Line 545  for (int j=0;j<4;++j) Line 462  for (int j=0;j<4;++j)
462              out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n3;              out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n3;
463              out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n7;              out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n7;
464         }         }
 // printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)],  out->FaceElements->Nodes[INDEX2(1,k,NN)],  
 // out->FaceElements->Nodes[INDEX2(2,k,NN)]);  
 // printf("%d: %d, %d, %d\n", k+1,out->FaceElements->Nodes[INDEX2(0,k+1,NN)],  out->FaceElements->Nodes[INDEX2(1,k+1,NN)],  
 // out->FaceElements->Nodes[INDEX2(2,k+1,NN)]);  
   
465               }               }
466             }             }
467             faceNECount+=2*local_NE1*local_NE2;             faceNECount+=2*local_NE1*local_NE2;
468           }           }
469           totalNECount+=2*NE1*NE2;           totalNECount+=2*NE1*NE2;
470       }       }
 //printf("\n");  
 //printf("Front\n");  
471       if (local_NE1>0) {       if (local_NE1>0) {
472          /* **  elements on boundary 010 (x2=0): */          /* **  elements on boundary 010 (x2=0): -Front*/
473          if (e_offset1 == 0) {          if (e_offset1 == 0) {
474             #pragma omp parallel for private(i0,i2,k,node0)             #pragma omp parallel for private(i0,i2,k,node0)
475             for (i2=0;i2<local_NE2;i2++) {             for (i2=0;i2<local_NE2;i2++) {
476               for (i0=0;i0<local_NE0;i0++) {               for (i0=0;i0<local_NE0;i0++) {
477                 k=2*(i0+local_NE0*i2)+faceNECount;                 k=2*(i0+local_NE0*i2)+faceNECount;
478                 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride2*N_PER_E*(i2+e_offset2);                 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride2*N_PER_E*(i2+e_offset2);
 //printf("%d, %d ",k,k+1);          
479             index_t res=2*((i2+e_offset2)+NE2*(e_offset0+i0))+totalNECount;             index_t res=2*((i2+e_offset2)+NE2*(e_offset0+i0))+totalNECount;
480                 out->FaceElements->Id[k]=res;                 out->FaceElements->Id[k]=res;
481                 out->FaceElements->Tag[k]=FRONTTAG;                 out->FaceElements->Tag[k]=FRONTTAG;
# Line 605  for (int j=0;j<4;++j) Line 514  for (int j=0;j<4;++j)
514              out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n4;              out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n4;
515    
516         }         }
 /*printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)],  out->FaceElements->Nodes[INDEX2(1,k,NN)],  
 out->FaceElements->Nodes[INDEX2(2,k,NN)]);  
 printf("%d: %d, %d, %d\n", k+1,out->FaceElements->Nodes[INDEX2(0,k+1,NN)],  out->FaceElements->Nodes[INDEX2(1,k+1,NN)],  
 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]);*/  
517               }               }
518             }             }
519             faceNECount+=2*local_NE0*local_NE2;             faceNECount+=2*local_NE0*local_NE2;
520          }          }
 // printf("\n");  
521          totalNECount+=2*NE0*NE2;          totalNECount+=2*NE0*NE2;
522  //printf("Back\n");          /* **  elements on boundary 020 (x2=1): - Back*/
         /* **  elements on boundary 020 (x2=1): */  
523          if (local_NE1+e_offset1 == NE1) {          if (local_NE1+e_offset1 == NE1) {
524             #pragma omp parallel for private(i0,i2,k,node0)             #pragma omp parallel for private(i0,i2,k,node0)
525             for (i2=0;i2<local_NE2;i2++) {             for (i2=0;i2<local_NE2;i2++) {
526               for (i0=0;i0<local_NE0;i0++) {               for (i0=0;i0<local_NE0;i0++) {
527                 k=2*(i0+local_NE0*i2)+faceNECount;                 k=2*(i0+local_NE0*i2)+faceNECount;
528                 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(NE1-1)+Nstride2*N_PER_E*(i2+e_offset2);                 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(NE1-1)+Nstride2*N_PER_E*(i2+e_offset2);
 // printf("%d, %d ",k,k+1);    
529             index_t res=2*((i2+e_offset2)+NE2*(i0+e_offset0))+totalNECount;             index_t res=2*((i2+e_offset2)+NE2*(i0+e_offset0))+totalNECount;
530                 out->FaceElements->Id[k]=res;                 out->FaceElements->Id[k]=res;
531                 out->FaceElements->Tag[k]=BACKTAG;                 out->FaceElements->Tag[k]=BACKTAG;
# Line 661  out->FaceElements->Nodes[INDEX2(2,k+1,NN Line 563  out->FaceElements->Nodes[INDEX2(2,k+1,NN
563              out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n3;              out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n3;
564    
565         }         }
 /*printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)],  out->FaceElements->Nodes[INDEX2(1,k,NN)],  
 out->FaceElements->Nodes[INDEX2(2,k,NN)]);  
 printf("%d: %d, %d, %d\n", k+1,out->FaceElements->Nodes[INDEX2(0,k+1,NN)],  out->FaceElements->Nodes[INDEX2(1,k+1,NN)],  
 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]);*/  
566               }               }
567             }             }
568             faceNECount+=2*local_NE0*local_NE2;             faceNECount+=2*local_NE0*local_NE2;
# Line 672  out->FaceElements->Nodes[INDEX2(2,k+1,NN Line 570  out->FaceElements->Nodes[INDEX2(2,k+1,NN
570          totalNECount+=2*NE0*NE2;          totalNECount+=2*NE0*NE2;
571       }       }
572    }    }
 // printf("\n");  
573    if (Dudley_noError()) {    if (Dudley_noError()) {
574       /* add tag names */       /* add tag names */
575       Dudley_Mesh_addTagMap(out,"top", TOPTAG);       Dudley_Mesh_addTagMap(out,"top", TOPTAG);
# Line 694  out->FaceElements->Nodes[INDEX2(2,k+1,NN Line 591  out->FaceElements->Nodes[INDEX2(2,k+1,NN
591        Dudley_Mesh_free(out);        Dudley_Mesh_free(out);
592    }    }
593      /* free up memory */      /* free up memory */
 //  Dudley_ReferenceElementSet_dealloc(refPoints);  
 //  Dudley_ReferenceElementSet_dealloc(refFaceElements);  
 //  Dudley_ReferenceElementSet_dealloc(refElements);  
594    Paso_MPIInfo_free( mpi_info );      Paso_MPIInfo_free( mpi_info );  
595    
596    return out;    return out;

Legend:
Removed from v.3220  
changed lines
  Added in v.3221

  ViewVC Help
Powered by ViewVC 1.1.26