/[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

branches/domexper/dudley/src/Mesh_hex8.c revision 3086 by jfenwick, Thu Aug 5 05:07:58 2010 UTC branches/domexper/dudley/src/Mesh_tet4.c revision 3114 by jfenwick, Fri Aug 27 05:26:25 2010 UTC
# Line 22  Line 22 
22    
23  /**************************************************************/  /**************************************************************/
24    
25  #include "RectangularMesh.h"  #include "TriangularMesh.h"
26    
27  Dudley_Mesh* Dudley_RectangularMesh_Hex8(dim_t* numElements,  
28    Dudley_Mesh* Dudley_TriangularMesh_Tet4(dim_t* numElements,
29                                            double* Length,                                            double* Length,
30                                            bool_t* periodic,                        index_t order,
                                           index_t order,  
31                                            index_t reduced_order,                                            index_t reduced_order,
                                           bool_t useElementsOnFace,  
                                           bool_t useFullElementOrder,  
32                                            bool_t optimize)                                            bool_t optimize)
33  {  {
34    #define N_PER_E 1    #define N_PER_E 1
# Line 38  Dudley_Mesh* Dudley_RectangularMesh_Hex8 Line 36  Dudley_Mesh* Dudley_RectangularMesh_Hex8
36    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;
37    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;
38    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;
39    Dudley_ReferenceElementSet *refPoints=NULL, *refContactElements=NULL, *refFaceElements=NULL, *refElements=NULL;    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 46  Dudley_Mesh* Dudley_RectangularMesh_Hex8 Line 44  Dudley_Mesh* Dudley_RectangularMesh_Hex8
44    double time0=Dudley_timer();    double time0=Dudley_timer();
45    #endif    #endif
46    
47      const int LEFTTAG=1;      /* boundary x1=0 */
48      const int RIGHTTAG=2;     /* boundary x1=1 */
49      const int BOTTOMTAG=100;  /* boundary x3=1 */
50      const int TOPTAG=200;     /* boundary x3=0 */
51      const int FRONTTAG=10;    /* boundary x2=0 */
52      const int BACKTAG=20;     /* boundary x2=1 */
53    
54    
55    
56    
57    /* get MPI information */    /* get MPI information */
58    mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );    mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
59    if (! Dudley_noError()) {    if (! Dudley_noError()) {
# Line 63  Dudley_Mesh* Dudley_RectangularMesh_Hex8 Line 71  Dudley_Mesh* Dudley_RectangularMesh_Hex8
71    N2=N_PER_E*NE2+1;    N2=N_PER_E*NE2+1;
72    
73    /*  allocate mesh: */      /*  allocate mesh: */  
74    sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);    sprintf(name,"Triangular %d x %d x %d (x 5) mesh",N0,N1,N2);
75    out=Dudley_Mesh_alloc(name,DIM, mpi_info);    out=Dudley_Mesh_alloc(name,DIM, mpi_info);
76    if (! Dudley_noError()) {    if (! Dudley_noError()) {
77        Paso_MPIInfo_free( mpi_info );        Paso_MPIInfo_free( mpi_info );
78        return NULL;        return NULL;
79    }    }
80    refElements= Dudley_ReferenceElementSet_alloc(Hex8,order,reduced_order);    refElements= Dudley_ReferenceElementSet_alloc(Tet4,order,reduced_order);
81    if (useElementsOnFace) {    refFaceElements=Dudley_ReferenceElementSet_alloc(Tri3, order, reduced_order);
         refFaceElements=Dudley_ReferenceElementSet_alloc(Hex8Face, order, reduced_order);  
         refContactElements=Dudley_ReferenceElementSet_alloc(Hex8Face_Contact, order, reduced_order);  
   } else {  
         refFaceElements=Dudley_ReferenceElementSet_alloc(Rec4, order, reduced_order);  
         refContactElements=Dudley_ReferenceElementSet_alloc(Rec4_Contact, order, reduced_order);  
   }  
82    refPoints=Dudley_ReferenceElementSet_alloc(Point1, order, reduced_order);    refPoints=Dudley_ReferenceElementSet_alloc(Point1, order, reduced_order);
83        
84    
85    if ( Dudley_noError()) {    if ( Dudley_noError()) {
86        
87        Dudley_Mesh_setPoints(out,Dudley_ElementFile_alloc(refPoints, mpi_info));        Dudley_Mesh_setPoints(out,Dudley_ElementFile_alloc(refPoints, mpi_info));
       Dudley_Mesh_setContactElements(out,Dudley_ElementFile_alloc(refContactElements, mpi_info));  
88        Dudley_Mesh_setFaceElements(out,Dudley_ElementFile_alloc(refFaceElements, mpi_info));        Dudley_Mesh_setFaceElements(out,Dudley_ElementFile_alloc(refFaceElements, mpi_info));
89        Dudley_Mesh_setElements(out,Dudley_ElementFile_alloc(refElements, mpi_info));        Dudley_Mesh_setElements(out,Dudley_ElementFile_alloc(refElements, mpi_info));
90    
# Line 126  Dudley_Mesh* Dudley_RectangularMesh_Hex8 Line 127  Dudley_Mesh* Dudley_RectangularMesh_Hex8
127        /* get the number of surface elements */        /* get the number of surface elements */
128    
129        NFaceElements=0;        NFaceElements=0;
130        if (!periodic[2] && (local_NE2>0)) {        if (local_NE2>0) {
131            NDOF2=N2;            NDOF2=N2;
132            if (offset2==0) NFaceElements+=local_NE1*local_NE0;            if (offset2==0) NFaceElements+=2*local_NE1*local_NE0;     // each face is split
133            if (local_NE2+e_offset2 == NE2) NFaceElements+=local_NE1*local_NE0;            if (local_NE2+e_offset2 == NE2) NFaceElements+=2*local_NE1*local_NE0;
134        } else {        } else {
135            NDOF2=N2-1;            NDOF2=N2-1;
136        }        }
137    
138        if (!periodic[0] && (local_NE0>0)) {        if (local_NE0>0) {
139            NDOF0=N0;            NDOF0=N0;
140            if (e_offset0 == 0) NFaceElements+=local_NE1*local_NE2;            if (e_offset0 == 0) NFaceElements+=2*local_NE1*local_NE2;
141            if (local_NE0+e_offset0 == NE0) NFaceElements+=local_NE1*local_NE2;            if (local_NE0+e_offset0 == NE0) NFaceElements+=2*local_NE1*local_NE2;
142        } else {        } else {
143            NDOF0=N0-1;            NDOF0=N0-1;
144        }        }
145           if (!periodic[1] && (local_NE1>0)) {          
146               NDOF1=N1;        if (local_NE1>0) {
147               if (e_offset1 == 0) NFaceElements+=local_NE0*local_NE2;          NDOF1=N1;
148               if (local_NE1+e_offset1 == NE1) NFaceElements+=local_NE0*local_NE2;          if (e_offset1 == 0) NFaceElements+=2*local_NE0*local_NE2;
149           } else {          if (local_NE1+e_offset1 == NE1) NFaceElements+=2*local_NE0*local_NE2;
150               NDOF1=N1-1;        } else {
151           }          NDOF1=N1-1;
152          }
153    }    }
154        
155        
# Line 155  Dudley_Mesh* Dudley_RectangularMesh_Hex8 Line 157  Dudley_Mesh* Dudley_RectangularMesh_Hex8
157    if (Dudley_noError()) {    if (Dudley_noError()) {
158    
159        Dudley_NodeFile_allocTable(out->Nodes,local_N0*local_N1*local_N2);        Dudley_NodeFile_allocTable(out->Nodes,local_N0*local_N1*local_N2);
160        Dudley_ElementFile_allocTable(out->Elements,local_NE0*local_NE1*local_NE2);      /* we split the rectangular prism this code used to produce into 5 tetrahedrons */
161          Dudley_ElementFile_allocTable(out->Elements,local_NE0*local_NE1*local_NE2*5);
162        /* each border face will be split in half */
163        Dudley_ElementFile_allocTable(out->FaceElements,NFaceElements);        Dudley_ElementFile_allocTable(out->FaceElements,NFaceElements);
164    }    }
165        
# Line 178  Dudley_Mesh* Dudley_RectangularMesh_Hex8 Line 182  Dudley_Mesh* Dudley_RectangularMesh_Hex8
182             out->Nodes->globalDegreesOfFreedom[k]=Nstride0*(global_i0%NDOF0)             out->Nodes->globalDegreesOfFreedom[k]=Nstride0*(global_i0%NDOF0)
183                                                 +Nstride1*(global_i1%NDOF1)                                                 +Nstride1*(global_i1%NDOF1)
184                                                 +Nstride2*(global_i2%NDOF2);                                                 +Nstride2*(global_i2%NDOF2);
185           printf("N=%d: %f,%f,%f\n", k, out->Nodes->Coordinates[INDEX2(0,k,DIM)],
186            out->Nodes->Coordinates[INDEX2(1,k,DIM)],
187            out->Nodes->Coordinates[INDEX2(2,k,DIM)]);
188           }           }
189         }         }
190       }       }
191    // printf("Now listing elements\n");
192       /*   set the elements: */       /*   set the elements: */
193    
194         int global_adjustment=0;   // If we are not the only rank we may need to shift our pattern to match neighbours
195    
196       NN=out->Elements->numNodes;       NN=out->Elements->numNodes;
197       #pragma omp parallel for private(i0,i1,i2,k,node0)       #pragma omp parallel for private(i0,i1,i2,k,node0)
198       for (i2=0;i2<local_NE2;i2++) {       for (i2=0;i2<local_NE2;i2++) {
199         for (i1=0;i1<local_NE1;i1++) {         for (i1=0;i1<local_NE1;i1++) {
200           for (i0=0;i0<local_NE0;i0++) {           for (i0=0;i0<local_NE0;i0++) {
201                        
202             k=i0+local_NE0*i1+local_NE0*local_NE1*i2;                     k=5*(i0+local_NE0*i1+local_NE0*local_NE1*i2);
203             node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);             node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
204      
205             out->Elements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1)+NE0*NE1*(i2+e_offset2);         index_t res=5*((i0+e_offset0)+NE0*(i1+e_offset1)+NE0*NE1*(i2+e_offset2));
206             out->Elements->Tag[k]=0;         for (int j=0;j<5;++j)
207             out->Elements->Owner[k]=myRank;         {
208                   out->Elements->Id[k+j]=res+j;
209             out->Elements->Nodes[INDEX2(0,k,NN)] =node0                             ;                 out->Elements->Tag[k+j]=0;
210             out->Elements->Nodes[INDEX2(1,k,NN)] =node0+                    Nstride0;                 out->Elements->Owner[k+j]=myRank;
211             out->Elements->Nodes[INDEX2(2,k,NN)] =node0+           Nstride1+Nstride0;         }
212             out->Elements->Nodes[INDEX2(3,k,NN)] =node0+           Nstride1;  
213             out->Elements->Nodes[INDEX2(4,k,NN)] =node0+Nstride2                   ;         index_t v[8];
214             out->Elements->Nodes[INDEX2(5,k,NN)] =node0+Nstride2         +Nstride0;  /*     in non-rotated orientation the points are numbered as follows:
215             out->Elements->Nodes[INDEX2(6,k,NN)] =node0+Nstride2+Nstride1+Nstride0;         The top face (clockwise= 0,1,3,2), bottom face (4,5,7,6)*/
216             out->Elements->Nodes[INDEX2(7,k,NN)] =node0+Nstride2+Nstride1           ;         if ((global_adjustment+i0+i1+i2)%2==0)
217           {
218    printf("Type0 %d, %d, %d\n",i0,i1,i2);
219           v[0]=node0;
220           v[1]=node0+Nstride0;
221           v[2]=node0+Nstride1;
222           v[3]=node0+Nstride1+Nstride0;
223           v[4]=node0+Nstride2;
224           v[5]=node0+Nstride0+Nstride2;
225           v[6]=node0+Nstride1+Nstride2;
226           v[7]=node0+Nstride2+Nstride1+Nstride0;
227           }
228           else
229           {
230    printf("Type1 %d, %d, %d\n",i0,i1,i2);
231           // this form is rotated around the 0,2,4,6 face clockwise 90 degrees
232    
233           v[0]=node0+Nstride2;
234           v[2]=node0;
235           v[6]=node0+Nstride1;
236           v[4]=node0+Nstride1+Nstride2;
237    
238           v[1]=node0+Nstride0+Nstride2;
239           v[3]=node0+Nstride0;
240           v[7]=node0+Nstride1+Nstride0;
241           v[5]=node0+Nstride2+Nstride1+Nstride0;
242           }
243    //     index_t a=node0, b=node0+Nstride0, c=node0+Nstride1+Nstride0, d=node0+Nstride1;
244    //     index_t e=node0+Nstride2, f=node0+Nstride2+Nstride0, g=node0+Nstride2+Nstride1+Nstride0,
245    //       h=node0+Nstride2+Nstride1;
246    //    
247    //
248    //
249    //     a=0, b=1, c=3, d=2
250    //     e=4, f=5, g=7, h=6
251    
252           // elements nodes are numbered: centre, x, y, z
253    
254           out->Elements->Nodes[INDEX2(0,k,NN)] =v[4];
255           out->Elements->Nodes[INDEX2(1,k,NN)] =v[5];
256           out->Elements->Nodes[INDEX2(2,k,NN)] =v[6];
257           out->Elements->Nodes[INDEX2(3,k,NN)] =v[0];
258    
259           out->Elements->Nodes[INDEX2(0,k+1,NN)] =v[7];
260           out->Elements->Nodes[INDEX2(1,k+1,NN)] =v[6];
261           out->Elements->Nodes[INDEX2(2,k+1,NN)] =v[5];
262           out->Elements->Nodes[INDEX2(3,k+1,NN)] =v[3];
263    
264           out->Elements->Nodes[INDEX2(0,k+2,NN)] =v[2];
265           out->Elements->Nodes[INDEX2(1,k+2,NN)] =v[3];
266           out->Elements->Nodes[INDEX2(2,k+2,NN)] =v[0];
267           out->Elements->Nodes[INDEX2(3,k+2,NN)] =v[6];
268    
269           out->Elements->Nodes[INDEX2(0,k+3,NN)] =v[1];
270           out->Elements->Nodes[INDEX2(1,k+3,NN)] =v[0];
271           out->Elements->Nodes[INDEX2(2,k+3,NN)] =v[3];
272           out->Elements->Nodes[INDEX2(3,k+3,NN)] =v[5];
273    
274           // I can't work out where the center is for this one
275           out->Elements->Nodes[INDEX2(0,k+4,NN)] =v[5];
276           out->Elements->Nodes[INDEX2(1,k+4,NN)] =v[0];
277           out->Elements->Nodes[INDEX2(2,k+4,NN)] =v[6];
278           out->Elements->Nodes[INDEX2(3,k+4,NN)] =v[3];
279    
280    
281    
282    for (int j=0;j<4;++j)
283    {
284    
285      printf("Elt %d",j);
286      for (int m=0;m<4;++m)
287      {
288         printf(" %d",out->Elements->Nodes[INDEX2(m,k+j,NN)]);
289      }
290      printf("\n");
291    }
292    
293    
294    
295           }           }
296    
297         }         }
298       }       } /* end for */
299       /* face elements */       /* face elements */
300    // printf("Starting face elements\n");
301       NN=out->FaceElements->numNodes;       NN=out->FaceElements->numNodes;
302       totalNECount=NE0*NE1*NE2;       totalNECount=5*NE0*NE1*NE2;
303       faceNECount=0;       faceNECount=0;
304    printf("Top\n");
305       /*   these are the quadrilateral elements on boundary 1 (x3=0): */       /*   these are the quadrilateral elements on boundary 1 (x3=0): */
306       if (!periodic[2]  && (local_NE2>0)) {       if (local_NE2>0)
307         {
308         /* **  elements on boundary 100 (x3=0): */         /* **  elements on boundary 100 (x3=0): */
309         if (e_offset2==0) {         if (e_offset2==0)
310           {
311            #pragma omp parallel for private(i0,i1,k,node0)            #pragma omp parallel for private(i0,i1,k,node0)
312            for (i1=0;i1<local_NE1;i1++) {            for (i1=0;i1<local_NE1;i1++)
313              for (i0=0;i0<local_NE0;i0++) {        {
314                          for (i0=0;i0<local_NE0;i0++)
315                k=i0+local_NE0*i1+faceNECount;          {
316                  k=2*(i0+local_NE0*i1)+faceNECount;
317                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);
318            
319                out->FaceElements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1)+totalNECount;            index_t res=2*((i0+e_offset0)+NE0*(i1+e_offset1))+totalNECount;
320                out->FaceElements->Tag[k]=100;                out->FaceElements->Id[k]=res;
321                  out->FaceElements->Tag[k]=TOPTAG;
322                out->FaceElements->Owner[k]=myRank;                out->FaceElements->Owner[k]=myRank;
323                          out->FaceElements->Id[k+1]=res+1;
324                if  (useElementsOnFace) {                out->FaceElements->Tag[k+1]=TOPTAG;
325                   out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0                                 ;                out->FaceElements->Owner[k+1]=myRank;
326                   out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0         +Nstride1           ;            // in the element generation this face is a,b,c,d which is split by a-c
327                   out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0         +Nstride1+Nstride0;  
328                   out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+          Nstride0           ;            index_t v[4];
329                   out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+Nstride2                      ;            if ((global_adjustment+i0+i1)%2==0)
330                   out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+Nstride2+Nstride1           ;            {
331                   out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+Nstride2+Nstride1+Nstride0;          v[0]=node0+Nstride1+Nstride2;   // node 6
332                   out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+Nstride2          +Nstride0;          v[1]=node0+Nstride2+Nstride1+Nstride0; // node 7
333                } else {          v[2]=node0+Nstride2;    // node 4
334                   out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0                                 ;          v[3]=node0+Nstride0+Nstride2;   // node 5
335                   out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+           Nstride1           ;  
336                   out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+           Nstride1+Nstride0;              out->FaceElements->Nodes[INDEX2(0,k,NN)]=v[0];
337                   out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+                     Nstride0;              out->FaceElements->Nodes[INDEX2(1,k,NN)]=v[1];
338                }              out->FaceElements->Nodes[INDEX2(2,k,NN)]=v[3];
339    
340                out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=v[0];
341                out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=v[3];
342                out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=v[2];
343    
344              }
345              else
346              {
347            v[0]=node0+Nstride2;        // node4
348            v[1]=node0+Nstride0+Nstride2;   // node5
349            v[2]=node0;         // node0
350            v[3]=node0+Nstride0;        // node1
351    
352                out->FaceElements->Nodes[INDEX2(0,k,NN)]=v[0];
353                out->FaceElements->Nodes[INDEX2(1,k,NN)]=v[1];
354                out->FaceElements->Nodes[INDEX2(2,k,NN)]=v[2];
355    
356                out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=v[2];
357                out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=v[3];
358                out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=v[1];
359              }
360    printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)],  out->FaceElements->Nodes[INDEX2(1,k,NN)],
361    out->FaceElements->Nodes[INDEX2(2,k,NN)]);
362    printf("%d: %d, %d, %d\n", k+1,out->FaceElements->Nodes[INDEX2(0,k+1,NN)],  out->FaceElements->Nodes[INDEX2(1,k+1,NN)],
363    out->FaceElements->Nodes[INDEX2(2,k+1,NN)]);
364    
365              }              }
366            }            }
367            faceNECount+=local_NE1*local_NE0;  // printf("\n");
368              faceNECount+=2*local_NE1*local_NE0;
369         }         }
370         totalNECount+=NE1*NE0;         totalNECount+=2*NE1*NE0;
371    
372    printf("Bottom\n");
373         /* **  elements on boundary 200 (x3=1) */         /* **  elements on boundary 200 (x3=1) */
374         if (local_NE2+e_offset2 == NE2) {         if (local_NE2+e_offset2 == NE2) {
375            #pragma omp parallel for private(i0,i1,k,node0)            #pragma omp parallel for private(i0,i1,k,node0)
376            for (i1=0;i1<local_NE1;i1++) {            for (i1=0;i1<local_NE1;i1++) {
377              for (i0=0;i0<local_NE0;i0++) {              for (i0=0;i0<local_NE0;i0++) {
378                
379                k=i0+local_NE0*i1+faceNECount;                k=2*(i0+local_NE0*i1)+faceNECount;
380                node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(NE2-1);                node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(NE2-1);
381            
382                out->FaceElements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1)+totalNECount;            index_t res=2*((i0+e_offset0)+NE0*(i1+e_offset1))+totalNECount;
383                out->FaceElements->Tag[k]=200;                out->FaceElements->Id[k]=res;
384                  out->FaceElements->Tag[k]=BOTTOMTAG;
385                out->FaceElements->Owner[k]=myRank;                out->FaceElements->Owner[k]=myRank;
386                if  (useElementsOnFace) {                out->FaceElements->Id[k+1]=res+1;
387                   out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0+Nstride2                      ;                out->FaceElements->Tag[k+1]=BOTTOMTAG;
388                   out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+Nstride2+           Nstride0;                out->FaceElements->Owner[k+1]=myRank;
389                   out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+Nstride2+Nstride1+Nstride0;  
390                   out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+Nstride2+Nstride1           ;            index_t v[4];
391                    if ((global_adjustment+i0+i1+local_NE2-1)%2==0)
392                   out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0                                 ;            {
393                   out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+Nstride0                      ;          v[0]=node0;     // node 0
394                   out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+           Nstride1+Nstride0;          v[1]=node0+Nstride0;    // node 1
395                   out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+           Nstride1;          v[2]=node0+Nstride1;    // node 2
396                } else {          v[3]=node0+Nstride1+Nstride0;   // node 3
397                   out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0+Nstride2                      ;  
398                   out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+Nstride2           +Nstride0;              out->FaceElements->Nodes[INDEX2(0,k,NN)]=v[0];
399                   out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+Nstride2+Nstride1+Nstride0;              out->FaceElements->Nodes[INDEX2(1,k,NN)]=v[3];
400                   out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+Nstride2+Nstride1           ;              out->FaceElements->Nodes[INDEX2(2,k,NN)]=v[2];
401                }  
402                out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=v[0];
403                out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=v[1];
404                out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=v[3];
405    
406    
407    
408              }
409              else
410              {
411            v[0]=node0+Nstride2;    // 4
412            v[1]=node0+Nstride0+Nstride2;   //5
413            v[2]=node0;     // 0
414            v[3]=node0+Nstride0;    // 1
415    
416    
417                out->FaceElements->Nodes[INDEX2(0,k,NN)]=v[0];
418                out->FaceElements->Nodes[INDEX2(1,k,NN)]=v[2];
419                out->FaceElements->Nodes[INDEX2(2,k,NN)]=v[1];
420    
421                out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=v[1];
422                out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=v[2];
423                out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=v[3];
424              }
425    
426    printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)],  out->FaceElements->Nodes[INDEX2(1,k,NN)],
427    out->FaceElements->Nodes[INDEX2(2,k,NN)]);
428    printf("%d: %d, %d, %d\n", k+1,out->FaceElements->Nodes[INDEX2(0,k+1,NN)],  out->FaceElements->Nodes[INDEX2(1,k+1,NN)],
429    out->FaceElements->Nodes[INDEX2(2,k+1,NN)]);
430    
431              }              }
432            }            }
433            faceNECount+=local_NE1*local_NE0;  /*printf("\n");*/
434              faceNECount+=2*local_NE1*local_NE0;
435         }         }
436         totalNECount+=NE1*NE0;         totalNECount+=2*NE1*NE0;
437       }       }
438       if (!periodic[0]  && (local_NE0>0)) {  printf("Left\n");
439         if (local_NE0>0) {
440          /* **  elements on boundary 001 (x1=0): */          /* **  elements on boundary 001 (x1=0): */
441            
442          if (e_offset0 == 0) {          if (e_offset0 == 0) {
# Line 287  Dudley_Mesh* Dudley_RectangularMesh_Hex8 Line 444  Dudley_Mesh* Dudley_RectangularMesh_Hex8
444             for (i2=0;i2<local_NE2;i2++) {             for (i2=0;i2<local_NE2;i2++) {
445               for (i1=0;i1<local_NE1;i1++) {               for (i1=0;i1<local_NE1;i1++) {
446                
447                 k=i1+local_NE1*i2+faceNECount;                 k=2*(i1+local_NE1*i2)+faceNECount;
448    
449    //printf("%d, %d ",k,k+1);
450                 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);
451                 out->FaceElements->Id[k]=(i1+e_offset1)+NE1*(i2+e_offset2)+totalNECount;             index_t res=2*((i1+e_offset1)+NE1*(i2+e_offset2))+totalNECount;
452                 out->FaceElements->Tag[k]=1;                 out->FaceElements->Id[k]=res;
453                   out->FaceElements->Tag[k]=LEFTTAG;
454                 out->FaceElements->Owner[k]=myRank;                 out->FaceElements->Owner[k]=myRank;
455                   out->FaceElements->Id[k+1]=res+1;
456                   out->FaceElements->Tag[k+1]=LEFTTAG;
457                   out->FaceElements->Owner[k+1]=myRank;
458                
459                 if  (useElementsOnFace) {            index_t v[4];
460                    out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0                                 ;             if ((global_adjustment+0+i1+i2)%2==0)
461                    out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+Nstride2                      ;         {
462                    out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+Nstride2+Nstride1           ;         v[0]=node0+Nstride1+Nstride2;    // 6
463                    out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+Nstride1                      ;         v[1]=node0+Nstride2;         // 4
464                    out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+Nstride0                      ;         v[2]=node0+Nstride1;         // 2
465                    out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+Nstride2+Nstride0           ;         v[3]=node0;
466                    out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+Nstride2+Nstride1+Nstride0;              out->FaceElements->Nodes[INDEX2(0,k,NN)]=v[0];
467                    out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+Nstride1+Nstride0           ;              out->FaceElements->Nodes[INDEX2(1,k,NN)]=v[3];
468                 } else {              out->FaceElements->Nodes[INDEX2(2,k,NN)]=v[1];
469                    out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0                                 ;  
470                    out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+Nstride2                      ;              out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=v[0];
471                    out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+Nstride2+Nstride1           ;              out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=v[2];
472                    out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+           Nstride1           ;              out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=v[3];
473                 }         }
474           else
475           {
476           // this form is rotated around the 0,2,4,6 face clockwise 90 degrees
477           v[0]=node0+Nstride1;     // 2
478           v[1]=node0+Nstride1+Nstride2;    // 6
479           v[2]=node0;      // 0
480           v[3]=node0+Nstride2; // 4
481    
482                out->FaceElements->Nodes[INDEX2(0,k,NN)]=v[0];
483                out->FaceElements->Nodes[INDEX2(1,k,NN)]=v[2];
484                out->FaceElements->Nodes[INDEX2(2,k,NN)]=v[1];
485    
486                out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=v[1];
487                out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=v[2];
488                out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=v[3];
489           }
490    printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)],  out->FaceElements->Nodes[INDEX2(1,k,NN)],
491    out->FaceElements->Nodes[INDEX2(2,k,NN)]);
492    printf("%d: %d, %d, %d\n", k+1,out->FaceElements->Nodes[INDEX2(0,k+1,NN)],  out->FaceElements->Nodes[INDEX2(1,k+1,NN)],
493    out->FaceElements->Nodes[INDEX2(2,k+1,NN)]);
494               }               }
495             }             }
496             faceNECount+=local_NE1*local_NE2;             faceNECount+=2*local_NE1*local_NE2;
497          }          }
498          totalNECount+=NE1*NE2;  /*printf("\n");*/
499    printf("Right\n");
500            totalNECount+=2*NE1*NE2;
501          /* **  elements on boundary 002 (x1=1): */          /* **  elements on boundary 002 (x1=1): */
502          if (local_NE0+e_offset0 == NE0) {          if (local_NE0+e_offset0 == NE0) {
503             #pragma omp parallel for private(i1,i2,k,node0)             #pragma omp parallel for private(i1,i2,k,node0)
504             for (i2=0;i2<local_NE2;i2++) {             for (i2=0;i2<local_NE2;i2++) {
505               for (i1=0;i1<local_NE1;i1++) {               for (i1=0;i1<local_NE1;i1++) {
506                 k=i1+local_NE1*i2+faceNECount;                 k=2*(i1+local_NE1*i2)+faceNECount;
507    
508    //printf("%d, %d ",k,k+1);
509                 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);
510                 out->FaceElements->Id[k]=(i1+e_offset1)+NE1*(i2+e_offset2)+totalNECount;             index_t res=2*((i1+e_offset1)+NE1*(i2+e_offset2))+totalNECount;
511                 out->FaceElements->Tag[k]=2;                 out->FaceElements->Id[k]=res;
512                   out->FaceElements->Tag[k]=RIGHTTAG;
513                 out->FaceElements->Owner[k]=myRank;                 out->FaceElements->Owner[k]=myRank;
514                     out->FaceElements->Id[k+1]=res+1;
515                 if  (useElementsOnFace) {                 out->FaceElements->Tag[k+1]=RIGHTTAG;
516                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+                      Nstride0;                 out->FaceElements->Owner[k+1]=myRank;
517                    out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+           Nstride1+Nstride0;  
518                    out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2+Nstride1+Nstride0;            index_t v[4];
519                    out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride2+           Nstride0;             if ((global_adjustment+local_NE0-1+i1+i2)%2==0)
520                 {
521                    out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0                                 ;         v[0]=node0+Nstride0+Nstride2;    // 5
522                    out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+           Nstride1           ;         v[1]=node0+Nstride2+Nstride1+Nstride0; // 7
523                    out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+Nstride2+Nstride1           ;         v[2]=node0+Nstride0;     // 1
524                    out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+Nstride2                      ;         v[3]=node0+Nstride1+Nstride0;    // 3
525          
526                 } else {  
527                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0                      +Nstride0;              out->FaceElements->Nodes[INDEX2(0,k,NN)]=v[0];
528                    out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+           Nstride1+Nstride0;              out->FaceElements->Nodes[INDEX2(1,k,NN)]=v[3];
529                    out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2+Nstride1+Nstride0;              out->FaceElements->Nodes[INDEX2(2,k,NN)]=v[1];
530                    out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride2+           Nstride0;  
531                 }              out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=v[0];
532                        out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=v[2];
533                out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=v[3];
534           }
535           else
536           {
537           // this form is rotated around the 0,2,4,6 face clockwise 90 degrees
538           v[0]=node0+Nstride2;         // 4
539           v[1]=node0+Nstride1+Nstride0;    // 3
540           v[2]=node0+Nstride0+Nstride2;    // 5
541           v[3]=node0+Nstride0;         // 1
542    
543                out->FaceElements->Nodes[INDEX2(0,k,NN)]=v[0];
544                out->FaceElements->Nodes[INDEX2(1,k,NN)]=v[2];
545                out->FaceElements->Nodes[INDEX2(2,k,NN)]=v[1];
546    
547                out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=v[1];
548                out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=v[2];
549                out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=v[3];
550           }
551    printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)],  out->FaceElements->Nodes[INDEX2(1,k,NN)],
552    out->FaceElements->Nodes[INDEX2(2,k,NN)]);
553    printf("%d: %d, %d, %d\n", k+1,out->FaceElements->Nodes[INDEX2(0,k+1,NN)],  out->FaceElements->Nodes[INDEX2(1,k+1,NN)],
554    out->FaceElements->Nodes[INDEX2(2,k+1,NN)]);
555    
556               }               }
557             }             }
558             faceNECount+=local_NE1*local_NE2;             faceNECount+=2*local_NE1*local_NE2;
559           }           }
560           totalNECount+=NE1*NE2;           totalNECount+=2*NE1*NE2;
561       }       }
562       if (!periodic[1] && (local_NE1>0)) {  //printf("\n");
563    printf("Front\n");
564         if (local_NE1>0) {
565          /* **  elements on boundary 010 (x2=0): */          /* **  elements on boundary 010 (x2=0): */
566          if (e_offset1 == 0) {          if (e_offset1 == 0) {
567             #pragma omp parallel for private(i0,i2,k,node0)             #pragma omp parallel for private(i0,i2,k,node0)
568             for (i2=0;i2<local_NE2;i2++) {             for (i2=0;i2<local_NE2;i2++) {
569               for (i0=0;i0<local_NE0;i0++) {               for (i0=0;i0<local_NE0;i0++) {
570                 k=i0+local_NE0*i2+faceNECount;                 k=2*(i0+local_NE0*i2)+faceNECount;
571                 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);
572            //printf("%d, %d ",k,k+1);        
573                 out->FaceElements->Id[k]=(i2+e_offset2)+NE2*(e_offset0+i0)+totalNECount;             index_t res=2*((i2+e_offset2)+NE2*(e_offset0+i0))+totalNECount;
574                 out->FaceElements->Tag[k]=10;                 out->FaceElements->Id[k]=res;
575                   out->FaceElements->Tag[k]=FRONTTAG;
576                 out->FaceElements->Owner[k]=myRank;                 out->FaceElements->Owner[k]=myRank;
577                 if  (useElementsOnFace) {                 out->FaceElements->Id[k+1]=res+1;
578                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0                                 ;                 out->FaceElements->Tag[k+1]=FRONTTAG;
579                    out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+                      Nstride0;                 out->FaceElements->Owner[k+1]=myRank;
580                    out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2           +Nstride0;  
581                    out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride2                      ;            index_t v[4];
582                     if ((global_adjustment+i0+0+i2)%2==0)
583                    out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+           Nstride1           ;         {
584                    out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+Nstride1+           Nstride0;         v[0]=node0+Nstride2;         // 4
585                    out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+Nstride2+Nstride1+Nstride0;         v[1]=node0+Nstride0+Nstride2;    // 5
586                    out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+Nstride2+Nstride1           ;         v[2]=node0;              // 0
587                 } else {         v[3]=node0+Nstride0;         // 1
588                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0                                 ;  
589                    out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+                      Nstride0;  
590                    out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2+           Nstride0;              out->FaceElements->Nodes[INDEX2(0,k,NN)]=v[0];
591                    out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride2                      ;              out->FaceElements->Nodes[INDEX2(1,k,NN)]=v[2];
592                 }              out->FaceElements->Nodes[INDEX2(2,k,NN)]=v[1];
593    
594                out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=v[2];
595                out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=v[3];
596                out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=v[1];
597    
598    
599           }
600           else
601           {
602           // this form is rotated around the 0,2,4,6 face clockwise 90 degrees
603           v[0]=node0+Nstride1+Nstride2;    // 6
604           v[1]=node0+Nstride2+Nstride1+Nstride0;   // 7
605           v[2]=node0+Nstride2;         // 4
606           v[3]=node0+Nstride0+Nstride2;    // 5
607    
608                out->FaceElements->Nodes[INDEX2(0,k,NN)]=v[0];
609                out->FaceElements->Nodes[INDEX2(1,k,NN)]=v[3];
610                out->FaceElements->Nodes[INDEX2(2,k,NN)]=v[2];
611    
612                out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=v[0];
613                out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=v[3];
614                out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=v[1];
615    
616           }
617    printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)],  out->FaceElements->Nodes[INDEX2(1,k,NN)],
618    out->FaceElements->Nodes[INDEX2(2,k,NN)]);
619    printf("%d: %d, %d, %d\n", k+1,out->FaceElements->Nodes[INDEX2(0,k+1,NN)],  out->FaceElements->Nodes[INDEX2(1,k+1,NN)],
620    out->FaceElements->Nodes[INDEX2(2,k+1,NN)]);
621               }               }
622             }             }
623             faceNECount+=local_NE0*local_NE2;             faceNECount+=2*local_NE0*local_NE2;
624          }          }
625          totalNECount+=NE0*NE2;  // printf("\n");
626            totalNECount+=2*NE0*NE2;
627    printf("Back\n");
628          /* **  elements on boundary 020 (x2=1): */          /* **  elements on boundary 020 (x2=1): */
629          if (local_NE1+e_offset1 == NE1) {          if (local_NE1+e_offset1 == NE1) {
630             #pragma omp parallel for private(i0,i2,k,node0)             #pragma omp parallel for private(i0,i2,k,node0)
631             for (i2=0;i2<local_NE2;i2++) {             for (i2=0;i2<local_NE2;i2++) {
632               for (i0=0;i0<local_NE0;i0++) {               for (i0=0;i0<local_NE0;i0++) {
633                 k=i0+local_NE0*i2+faceNECount;                 k=2*(i0+local_NE0*i2)+faceNECount;
634                 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);
635      // printf("%d, %d ",k,k+1);  
636                 out->FaceElements->Id[k]=(i2+e_offset2)+NE2*(i0+e_offset0)+totalNECount;             index_t res=2*((i2+e_offset2)+NE2*(i0+e_offset0))+totalNECount;
637                 out->FaceElements->Tag[k]=20;                 out->FaceElements->Id[k]=res;
638                   out->FaceElements->Tag[k]=BACKTAG;
639                 out->FaceElements->Owner[k]=myRank;                 out->FaceElements->Owner[k]=myRank;
640                         out->FaceElements->Id[k+1]=res+1;
641                 if  (useElementsOnFace) {                 out->FaceElements->Tag[k+1]=BACKTAG;
642                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+           Nstride1           ;                 out->FaceElements->Owner[k+1]=myRank;
643                    out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride2+Nstride1           ;  
644                    out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2+Nstride1+Nstride0;            index_t v[4];
645                    out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride1+Nstride0           ;             if ((global_adjustment+i0+local_NE1-1+i2)%2==0)
646                 {
647                    out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0                                 ;         v[0]=node0+Nstride2+Nstride1+Nstride0;   // 7
648                    out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+Nstride2                      ;         v[1]=node0+Nstride1+Nstride2;        // 6
649                    out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+Nstride2+           Nstride0;         v[2]=node0+Nstride1+Nstride0;        // 3
650                    out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+                      Nstride0;         v[3]=node0+Nstride1;             // 2
651                 } else {  
652                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+           Nstride1           ;  
653                    out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride2+Nstride1           ;              out->FaceElements->Nodes[INDEX2(0,k,NN)]=v[0];
654                    out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2+Nstride1+Nstride0;              out->FaceElements->Nodes[INDEX2(1,k,NN)]=v[2];
655                    out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+           Nstride1+Nstride0;              out->FaceElements->Nodes[INDEX2(2,k,NN)]=v[1];
656                 }  
657                out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=v[2];
658                out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=v[3];
659                out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=v[1];
660    
661    
662           }
663           else
664           {
665           // this form is rotated around the 0,2,4,6 face clockwise 90 degrees
666           v[0]=node0+Nstride1+Nstride0;    // 3
667           v[1]=node0+Nstride1;         // 2
668           v[2]=node0+Nstride0;         // 1
669           v[3]=node0;              // 0
670    
671                out->FaceElements->Nodes[INDEX2(0,k,NN)]=v[0];
672                out->FaceElements->Nodes[INDEX2(1,k,NN)]=v[2];
673                out->FaceElements->Nodes[INDEX2(2,k,NN)]=v[3];
674    
675                out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=v[0];
676                out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=v[3];
677                out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=v[1];
678    
679           }
680    printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)],  out->FaceElements->Nodes[INDEX2(1,k,NN)],
681    out->FaceElements->Nodes[INDEX2(2,k,NN)]);
682    printf("%d: %d, %d, %d\n", k+1,out->FaceElements->Nodes[INDEX2(0,k+1,NN)],  out->FaceElements->Nodes[INDEX2(1,k+1,NN)],
683    out->FaceElements->Nodes[INDEX2(2,k+1,NN)]);
684               }               }
685             }             }
686             faceNECount+=local_NE0*local_NE2;             faceNECount+=2*local_NE0*local_NE2;
687          }          }
688          totalNECount+=NE0*NE2;          totalNECount+=2*NE0*NE2;
689       }       }
690    }    }
691    // printf("\n");
692    if (Dudley_noError()) {    if (Dudley_noError()) {
693       /* add tag names */       /* add tag names */
694       Dudley_Mesh_addTagMap(out,"top", 200);       Dudley_Mesh_addTagMap(out,"top", TOPTAG);
695       Dudley_Mesh_addTagMap(out,"bottom", 100);       Dudley_Mesh_addTagMap(out,"bottom", BOTTOMTAG);
696       Dudley_Mesh_addTagMap(out,"left", 1);       Dudley_Mesh_addTagMap(out,"left", LEFTTAG);
697       Dudley_Mesh_addTagMap(out,"right", 2);       Dudley_Mesh_addTagMap(out,"right", RIGHTTAG);
698       Dudley_Mesh_addTagMap(out,"front", 10);       Dudley_Mesh_addTagMap(out,"front", FRONTTAG);
699       Dudley_Mesh_addTagMap(out,"back", 20);       Dudley_Mesh_addTagMap(out,"back", BACKTAG);
700    }    }
701    /* prepare mesh for further calculatuions:*/    /* prepare mesh for further calculatuions:*/
702    if (Dudley_noError()) {    if (Dudley_noError()) {
# Line 438  Dudley_Mesh* Dudley_RectangularMesh_Hex8 Line 711  Dudley_Mesh* Dudley_RectangularMesh_Hex8
711    }    }
712      /* free up memory */      /* free up memory */
713    Dudley_ReferenceElementSet_dealloc(refPoints);    Dudley_ReferenceElementSet_dealloc(refPoints);
   Dudley_ReferenceElementSet_dealloc(refContactElements);  
714    Dudley_ReferenceElementSet_dealloc(refFaceElements);    Dudley_ReferenceElementSet_dealloc(refFaceElements);
715    Dudley_ReferenceElementSet_dealloc(refElements);    Dudley_ReferenceElementSet_dealloc(refElements);
716    Paso_MPIInfo_free( mpi_info );      Paso_MPIInfo_free( mpi_info );  
717    
718    return out;    return out;
719  }  }
720    
721    

Legend:
Removed from v.3086  
changed lines
  Added in v.3114

  ViewVC Help
Powered by ViewVC 1.1.26