/[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 3117 by jfenwick, Mon Aug 30 02:10:26 2010 UTC revision 3118 by jfenwick, Mon Aug 30 04:26:17 2010 UTC
# Line 25  Line 25 
25  #include "TriangularMesh.h"  #include "TriangularMesh.h"
26    
27    
28    /* Be careful reading this function. The X? and NStride? are 1,2,3 but the loop vars are 0,1,2 */
29  Dudley_Mesh* Dudley_TriangularMesh_Tet4(dim_t* numElements,  Dudley_Mesh* Dudley_TriangularMesh_Tet4(dim_t* numElements,
30                                            double* Length,                                            double* Length,
31                        index_t order,                        index_t order,
# Line 182  Dudley_Mesh* Dudley_TriangularMesh_Tet4( Line 183  Dudley_Mesh* Dudley_TriangularMesh_Tet4(
183             out->Nodes->globalDegreesOfFreedom[k]=Nstride0*(global_i0%NDOF0)             out->Nodes->globalDegreesOfFreedom[k]=Nstride0*(global_i0%NDOF0)
184                                                 +Nstride1*(global_i1%NDOF1)                                                 +Nstride1*(global_i1%NDOF1)
185                                                 +Nstride2*(global_i2%NDOF2);                                                 +Nstride2*(global_i2%NDOF2);
186  /*     printf("N=%d: %f,%f,%f\n", k, out->Nodes->Coordinates[INDEX2(0,k,DIM)],  //     printf("N=%d: %f,%f,%f\n", k, out->Nodes->Coordinates[INDEX2(0,k,DIM)],
187          out->Nodes->Coordinates[INDEX2(1,k,DIM)],  //      out->Nodes->Coordinates[INDEX2(1,k,DIM)],
188          out->Nodes->Coordinates[INDEX2(2,k,DIM)]);*/  //      out->Nodes->Coordinates[INDEX2(2,k,DIM)]);
189           }           }
190         }         }
191       }       }
# Line 322  for (int j=0;j<4;++j) Line 323  for (int j=0;j<4;++j)
323       NN=out->FaceElements->numNodes;       NN=out->FaceElements->numNodes;
324       totalNECount=5*NE0*NE1*NE2;       totalNECount=5*NE0*NE1*NE2;
325       faceNECount=0;       faceNECount=0;
326  //printf("Top\n");  //printf("Bottom\n");
327       /*   these are the quadrilateral elements on boundary 1 (x3=0): */       /*   these are the quadrilateral elements on boundary 1 (x3=0): */
328       if (local_NE2>0)       if (local_NE2>0)
329       {       {
# Line 336  for (int j=0;j<4;++j) Line 337  for (int j=0;j<4;++j)
337          {          {
338                k=2*(i0+local_NE0*i1)+faceNECount;                k=2*(i0+local_NE0*i1)+faceNECount;
339                node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1);                node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1);
       
340            index_t res=2*((i0+e_offset0)+NE0*(i1+e_offset1))+totalNECount;            index_t res=2*((i0+e_offset0)+NE0*(i1+e_offset1))+totalNECount;
341                out->FaceElements->Id[k]=res;                out->FaceElements->Id[k]=res;
342                out->FaceElements->Tag[k]=TOPTAG;                out->FaceElements->Tag[k]=BOTTOMTAG;
343                out->FaceElements->Owner[k]=myRank;                out->FaceElements->Owner[k]=myRank;
344                out->FaceElements->Id[k+1]=res+1;                out->FaceElements->Id[k+1]=res+1;
345                out->FaceElements->Tag[k+1]=TOPTAG;                out->FaceElements->Tag[k+1]=BOTTOMTAG;
346                out->FaceElements->Owner[k+1]=myRank;                out->FaceElements->Owner[k+1]=myRank;
           // in the element generation this face is a,b,c,d which is split by a-c  
347    
348            index_t n4=node0+Nstride2;            index_t n0=node0;
349            index_t n5=node0+Nstride0+Nstride2;            index_t n1=node0+Nstride0;
350            index_t n6=node0+Nstride1+Nstride2;            index_t n2=node0+Nstride1;
351            index_t n7=node0+Nstride0+Nstride1+Nstride2;            index_t n3=node0+Nstride0+Nstride1;
352    
353            if ((global_adjustment+i0+i1)%2==0)            if ((global_adjustment+i0+i1)%2==0)
354            {            {
355              out->FaceElements->Nodes[INDEX2(0,k,NN)]=n4;              out->FaceElements->Nodes[INDEX2(0,k,NN)]=n0;
356              out->FaceElements->Nodes[INDEX2(1,k,NN)]=n5;              out->FaceElements->Nodes[INDEX2(1,k,NN)]=n3;
357              out->FaceElements->Nodes[INDEX2(2,k,NN)]=n6;              out->FaceElements->Nodes[INDEX2(2,k,NN)]=n1;
   
             out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n5;  
             out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n7;  
             out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n6;  
358    
359                out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n0;
360                out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n2;
361                out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n3;
362    
363            }            }
364            else            else
365            {            {
366              out->FaceElements->Nodes[INDEX2(0,k,NN)]=n4;              out->FaceElements->Nodes[INDEX2(0,k,NN)]=n0;
367              out->FaceElements->Nodes[INDEX2(1,k,NN)]=n5;              out->FaceElements->Nodes[INDEX2(1,k,NN)]=n2;
368              out->FaceElements->Nodes[INDEX2(2,k,NN)]=n7;              out->FaceElements->Nodes[INDEX2(2,k,NN)]=n1;
369    
370              out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n4;              out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n1;
371              out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n7;              out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n2;
372              out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n6;              out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n3;
373    
374    
375    
# Line 388  for (int j=0;j<4;++j) Line 386  for (int j=0;j<4;++j)
386         }         }
387         totalNECount+=2*NE1*NE0;         totalNECount+=2*NE1*NE0;
388    
389  //printf("Bottom\n");  //printf("Top\n");
390         /* **  elements on boundary 200 (x3=1) */         /* **  elements on boundary 200 (x3=1) */
391         if (local_NE2+e_offset2 == NE2) {         if (local_NE2+e_offset2 == NE2) {
392            #pragma omp parallel for private(i0,i1,k,node0)            #pragma omp parallel for private(i0,i1,k,node0)
# Line 400  for (int j=0;j<4;++j) Line 398  for (int j=0;j<4;++j)
398    
399            index_t res=2*((i0+e_offset0)+NE0*(i1+e_offset1))+totalNECount;            index_t res=2*((i0+e_offset0)+NE0*(i1+e_offset1))+totalNECount;
400                out->FaceElements->Id[k]=res;                out->FaceElements->Id[k]=res;
401                out->FaceElements->Tag[k]=BOTTOMTAG;                out->FaceElements->Tag[k]=TOPTAG;
402                out->FaceElements->Owner[k]=myRank;                out->FaceElements->Owner[k]=myRank;
403                out->FaceElements->Id[k+1]=res+1;                out->FaceElements->Id[k+1]=res+1;
404                out->FaceElements->Tag[k+1]=BOTTOMTAG;                out->FaceElements->Tag[k+1]=TOPTAG;
405                out->FaceElements->Owner[k+1]=myRank;                out->FaceElements->Owner[k+1]=myRank;
406    
407            index_t n0=node0;            index_t n4=node0+Nstride2;
408            index_t n1=node0+Nstride0;            index_t n5=node0+Nstride0+Nstride2;
409            index_t n2=node0+Nstride1;            index_t n6=node0+Nstride1+Nstride2;
410            index_t n3=node0+Nstride1+Nstride0;            index_t n7=node0+Nstride1+Nstride0+Nstride2;
411    
412            if ((global_adjustment+i0+i1+local_NE2-1)%2==0)            if ((global_adjustment+i0+i1+local_NE2-1)%2==0)
413            {            {
414  //          out->FaceElements->Nodes[INDEX2(0,k,NN)]=v[0];              out->FaceElements->Nodes[INDEX2(0,k,NN)]=n4;
415  //          out->FaceElements->Nodes[INDEX2(1,k,NN)]=v[3];              out->FaceElements->Nodes[INDEX2(1,k,NN)]=n5;
416  //          out->FaceElements->Nodes[INDEX2(2,k,NN)]=v[2];              out->FaceElements->Nodes[INDEX2(2,k,NN)]=n6;
 //  
 //          out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=v[0];  
 //          out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=v[1];  
 //          out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=v[3];  
   
             out->FaceElements->Nodes[INDEX2(0,k,NN)]=n0;  
             out->FaceElements->Nodes[INDEX2(1,k,NN)]=n3;  
             out->FaceElements->Nodes[INDEX2(2,k,NN)]=n2;  
   
             out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n0;  
             out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n1;  
             out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n3;  
   
417    
418                out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n5;
419                out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n7;
420                out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n6;
421            }            }
422            else            else
423            {            {
424    
425              out->FaceElements->Nodes[INDEX2(0,k,NN)]=n0;              out->FaceElements->Nodes[INDEX2(0,k,NN)]=n4;
426              out->FaceElements->Nodes[INDEX2(1,k,NN)]=n1;              out->FaceElements->Nodes[INDEX2(1,k,NN)]=n5;
427              out->FaceElements->Nodes[INDEX2(2,k,NN)]=n2;              out->FaceElements->Nodes[INDEX2(2,k,NN)]=n7;
428    
429              out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n1;              out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n4;
430              out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n3;              out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n7;
431              out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n2;              out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n6;
432            }            }
433    
434  // printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)],  out->FaceElements->Nodes[INDEX2(1,k,NN)],  // printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)],  out->FaceElements->Nodes[INDEX2(1,k,NN)],

Legend:
Removed from v.3117  
changed lines
  Added in v.3118

  ViewVC Help
Powered by ViewVC 1.1.26