/[escript]/trunk/finley/src/Mesh_rec8.c
ViewVC logotype

Diff of /trunk/finley/src/Mesh_rec8.c

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

revision 471 by jgs, Fri Jan 27 01:33:02 2006 UTC revision 1388 by trankine, Fri Jan 11 07:45:58 2008 UTC
# Line 1  Line 1 
1  /*  
2   ******************************************************************************  /* $Id$ */
3   *                                                                            *  
4   *       COPYRIGHT  ACcESS 2003,2004,2005 -  All Rights Reserved              *  /*******************************************************
5   *                                                                            *   *
6   * This software is the property of ACcESS. No part of this code              *   *           Copyright 2003-2007 by ACceSS MNRF
7   * may be copied in any form or by any means without the expressed written    *   *       Copyright 2007 by University of Queensland
8   * consent of ACcESS.  Copying, use or modification of this software          *   *
9   * by any unauthorised person is illegal unless that person has a software    *   *                http://esscc.uq.edu.au
10   * license agreement with ACcESS.                                             *   *        Primary Business: Queensland, Australia
11   *                                                                            *   *  Licensed under the Open Software License version 3.0
12   ******************************************************************************   *     http://www.opensource.org/licenses/osl-3.0.php
13  */   *
14     *******************************************************/
15    
16  /**************************************************************/  /**************************************************************/
17    
# Line 21  Line 22 
22    
23  /**************************************************************/  /**************************************************************/
24    
 /*   Author: gross@access.edu.au */  
 /*   Version: $Id$ */  
   
 /**************************************************************/  
   
25  #include "RectangularMesh.h"  #include "RectangularMesh.h"
26    
 /**************************************************************/  
27    
28  Finley_Mesh* Finley_RectangularMesh_Rec8(int* numElements,double* Length,int* periodic,int order,int useElementsOnFace) {  Finley_Mesh* Finley_RectangularMesh_Rec8(dim_t* numElements,
29    dim_t N0,N1,NE0,NE1,i0,i1,totalNECount,faceNECount,NDOF0,NDOF1,NFaceElements,NUMNODES,M0,M1;                                            double* Length,
30    index_t k,node0;                                            bool_t* periodic,
31                                              index_t order,
32                                              index_t reduced_order,
33                                              bool_t useElementsOnFace,
34                                              bool_t useFullElementOrder,
35                                              bool_t optimize)
36    {
37      #define N_PER_E 2
38      #define DIM 2
39      dim_t N0,N1,NE0,NE1,i0,i1,k,Nstride0,Nstride1;
40      dim_t totalNECount,faceNECount,NDOF0,NDOF1,NFaceElements,NN, local_NE0, local_NE1, local_N0, local_N1;
41      index_t e_offset1, e_offset0, offset0, offset1, global_i0, global_i1;
42      index_t node0, myRank;
43    Finley_Mesh* out;    Finley_Mesh* out;
44      Paso_MPIInfo *mpi_info = NULL;
45    char name[50];    char name[50];
46    double time0=Finley_timer();    double time0=Finley_timer();
47    
48      /* get MPI information */
49      mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
50      if (! Finley_noError()) {
51            return NULL;
52      }
53      myRank=mpi_info->rank;
54    
55      /* set up the global dimensions of the mesh */
56    
57    NE0=MAX(1,numElements[0]);    NE0=MAX(1,numElements[0]);
58    NE1=MAX(1,numElements[1]);    NE1=MAX(1,numElements[1]);
59    N0=2*NE0+1;    N0=N_PER_E*NE0+1;
60    N1=2*NE1+1;    N1=N_PER_E*NE1+1;
61    
62    if (N0<=N1) {    /*  allocate mesh: */  
63       M0=1;    sprintf(name,"Rectangular %d x %d mesh",N0,N1);
64       M1=N0;    out=Finley_Mesh_alloc(name,DIM,order, reduced_order, mpi_info);
65    } else {    if (! Finley_noError()) {
66       M0=N1;        Paso_MPIInfo_free( mpi_info );
67       M1=1;        return NULL;
68      }
69    
70      if (useFullElementOrder) {
71         Finley_Mesh_setElements(out,Finley_ElementFile_alloc(Rec9,
72                                                out->order,
73                                                out->reduced_order,
74                                                mpi_info));
75         if (useElementsOnFace) {
76             Finley_setError(SYSTEM_ERROR,"rich elements for Rec9 elements is not supported yet.");
77         } else {
78             Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Line3,
79                                                        out->order,
80                                                        out->reduced_order,
81                                                        mpi_info));
82             Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Line3_Contact,
83                                                           out->order,
84                                                           out->reduced_order,
85                                                           mpi_info));
86         }
87      } else  {
88         Finley_Mesh_setElements(out,Finley_ElementFile_alloc(Rec8,out->order,out->reduced_order,mpi_info));
89         if (useElementsOnFace) {
90             Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Rec8Face,
91                                                                      out->order,
92                                                                      out->reduced_order,
93                                                                      mpi_info));
94             Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Rec8Face_Contact,
95                                                                        out->order,
96                                                                        out->reduced_order,
97                                                                        mpi_info));
98         } else {
99             Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Line3,
100                                                                      out->order,
101                                                                      out->reduced_order,
102                                                                      mpi_info));
103             Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Line3_Contact,
104                                                                         out->order,
105                                                                         out->reduced_order,
106                                                                         mpi_info));
107         }
108      }
109      Finley_Mesh_setPoints(out,Finley_ElementFile_alloc(Point1,
110                                                     out->order,
111                                                     out->reduced_order,
112                                                     mpi_info));
113      if (! Finley_noError()) {
114          Paso_MPIInfo_free( mpi_info );
115          Finley_Mesh_free(out);
116          return NULL;
117    }    }
118    
119      /* work out the largest dimension */
120      if (N1==MAX(N0,N1)) {
121         Nstride0=1;
122         Nstride1=N0;
123         local_NE0=NE0;
124         e_offset0=0;
125         Paso_MPIInfo_Split(mpi_info,NE1,&local_NE1,&e_offset1);
126      } else {
127         Nstride0=N1;
128         Nstride1=1;
129         Paso_MPIInfo_Split(mpi_info,NE0,&local_NE0,&e_offset0);
130         local_NE1=NE1;
131         e_offset1=0;
132      }
133      offset0=e_offset0*N_PER_E;
134      offset1=e_offset1*N_PER_E;
135      local_N0=local_NE0*N_PER_E+1;
136      local_N1=local_NE1*N_PER_E+1;
137    
138      /* get the number of surface elements */
139    
140    NFaceElements=0;    NFaceElements=0;
141    if (!periodic[0]) {    if (!periodic[0]) {
142        NDOF0=N0;       NDOF0=N0;
143        NFaceElements+=2*NE1;       if (e_offset0 == 0) NFaceElements+=local_NE1;
144         if (local_NE0+e_offset0 == NE0) NFaceElements+=local_NE1;
145    } else {    } else {
146        NDOF0=N0-1;        NDOF0=N0-1;
147    }    }
148    if (!periodic[1]) {    if (!periodic[1]) {
149        NDOF1=N1;       NDOF1=N1;
150        NFaceElements+=2*NE0;       if (e_offset1 == 0) NFaceElements+=local_NE0;
151         if (local_NE1+e_offset1 == NE1) NFaceElements+=local_NE0;
152    } else {    } else {
153        NDOF1=N1-1;        NDOF1=N1-1;
154    }    }
155    
   /*  allocate mesh: */  
     
   sprintf(name,"Rectangular %d x %d mesh",N0,N1);  
   out=Finley_Mesh_alloc(name,2,order);  
   if (! Finley_noError()) return NULL;  
   
   out->Elements=Finley_ElementFile_alloc(Rec8,out->order);  
   if (useElementsOnFace) {  
      out->FaceElements=Finley_ElementFile_alloc(Rec8Face,out->order);  
      out->ContactElements=Finley_ElementFile_alloc(Rec8Face_Contact,out->order);  
   } else {  
      out->FaceElements=Finley_ElementFile_alloc(Line3,out->order);  
      out->ContactElements=Finley_ElementFile_alloc(Line3_Contact,out->order);  
   }  
   out->Points=Finley_ElementFile_alloc(Point1,out->order);  
   if (! Finley_noError()) {  
        Finley_Mesh_dealloc(out);  
        return NULL;  
   }  
     
156    /*  allocate tables: */    /*  allocate tables: */
     
   Finley_NodeFile_allocTable(out->Nodes,N0*N1);  
   Finley_ElementFile_allocTable(out->Elements,NE0*NE1);  
   Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);  
   if (! Finley_noError()) {  
       Finley_Mesh_dealloc(out);  
       return NULL;  
   }  
   /*  set nodes: */  
                                                                                                                                                                                                       
   #pragma omp parallel for private(i0,i1,k)  
   for (i1=0;i1<N1;i1++) {  
     for (i0=0;i0<N0;i0++) {  
       k=M0*i0+M1*i1;  
       out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(i0)/DBLE(N0-1)*Length[0];  
       out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];  
       out->Nodes->Id[k]=i0+N0*i1;  
       out->Nodes->Tag[k]=0;  
       out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1);  
     }  
   }  
   /* tags for the faces: */  
   if (!periodic[1]) {  
     for (i0=0;i0<N0;i0++) {  
       out->Nodes->Tag[M0*i0+M1*0]+=10;  
       out->Nodes->Tag[M0*i0+M1*(N1-1)]+=20;  
     }  
   }  
   if (!periodic[0]) {  
     for (i1=0;i1<N1;i1++) {  
       out->Nodes->Tag[M0*0+M1*i1]+=1;  
       out->Nodes->Tag[M0*(N0-1)+M1*i1]+=2;  
     }  
   }  
157    
158    /*   set the elements: */    Finley_NodeFile_allocTable(out->Nodes,local_N0*local_N1);
159        Finley_ElementFile_allocTable(out->Elements,local_NE0*local_NE1);
160    #pragma omp parallel for private(i0,i1,k,node0)    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
   for (i1=0;i1<NE1;i1++) {  
     for (i0=0;i0<NE0;i0++) {  
       k=i0+NE0*i1;  
       node0=2*i0+2*i1*N0;  
   
       out->Elements->Id[k]=k;  
       out->Elements->Tag[k]=0;  
       out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1);  
   
       out->Elements->Nodes[INDEX2(0,k,8)]=node0;  
       out->Elements->Nodes[INDEX2(1,k,8)]=node0+2;  
       out->Elements->Nodes[INDEX2(2,k,8)]=node0+2*N0+2;  
       out->Elements->Nodes[INDEX2(3,k,8)]=node0+2*N0;  
       out->Elements->Nodes[INDEX2(4,k,8)]=node0+1;  
       out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0+2;  
       out->Elements->Nodes[INDEX2(6,k,8)]=node0+2*N0+1;  
       out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0;  
161    
162      }    if (Finley_noError()) {
163    }       /* create nodes */
164    out->Elements->minColor=0;    
165    out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0);       #pragma omp parallel for private(i0,i1,k,global_i0,global_i1)
166           for (i1=0;i1<local_N1;i1++) {
167    /*   face elements: */         for (i0=0;i0<local_N0;i0++) {
168    if (useElementsOnFace) {             k=i0+local_N0*i1;
169       NUMNODES=8;             global_i0=i0+offset0;
170    } else {             global_i1=i1+offset1;
171       NUMNODES=3;             out->Nodes->Coordinates[INDEX2(0,k,DIM)]=DBLE(global_i0)/DBLE(N0-1)*Length[0];
172    }             out->Nodes->Coordinates[INDEX2(1,k,DIM)]=DBLE(global_i1)/DBLE(N1-1)*Length[1];
173                 out->Nodes->Id[k]=Nstride0*global_i0+Nstride1*global_i1;
174                 out->Nodes->Tag[k]=0;
175    totalNECount=NE0*NE1;             out->Nodes->globalDegreesOfFreedom[k]=Nstride0*(global_i0%NDOF0)
176    faceNECount=0;                                                 +Nstride1*(global_i1%NDOF1);
     
   if (!periodic[1]) {  
      /* **  elements on boundary 010 (x2=0): */  
     
      #pragma omp parallel for private(i0,k,node0)  
      for (i0=0;i0<NE0;i0++) {  
        k=i0+faceNECount;  
        node0=2*i0;  
   
        out->FaceElements->Id[k]=i0+totalNECount;  
        out->FaceElements->Tag[k]=10;  
        out->FaceElements->Color[k]=i0%2;  
   
        if (useElementsOnFace) {  
           out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
           out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2;  
           out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0+2;  
           out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0;  
           out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;  
           out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0+2;  
           out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0+1;  
           out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0;  
        } else {  
           out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
           out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2;  
           out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;  
        }  
      }  
      totalNECount+=NE0;  
      faceNECount+=NE0;  
     
      /* **  elements on boundary 020 (x2=1): */  
     
      #pragma omp parallel for private(i0,k,node0)  
      for (i0=0;i0<NE0;i0++) {  
        k=i0+faceNECount;  
        node0=2*i0+2*(NE1-1)*N0;  
   
        out->FaceElements->Id[k]=i0+totalNECount;  
        out->FaceElements->Tag[k]=20;  
        out->FaceElements->Color[k]=i0%2+2;  
        if (useElementsOnFace) {  
           out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0+2;  
           out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0;  
           out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0;  
           out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2;  
           out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+2*N0+1;  
           out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0;  
           out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+1;  
           out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0+2;  
        } else {  
           out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0+2;  
           out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0;  
           out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0+1;  
177         }         }
178       }       }
179       totalNECount+=NE0;       /*   set the elements: */
180       faceNECount+=NE0;       NN=out->Elements->numNodes;
181    }       #pragma omp parallel for private(i0,i1,k,node0)
182    if (!periodic[0]) {       for (i1=0;i1<local_NE1;i1++) {
183       /* **  elements on boundary 001 (x1=0): */           for (i0=0;i0<local_NE0;i0++) {
184                
185       #pragma omp parallel for private(i1,k,node0)             k=i0+local_NE0*i1;        
186       for (i1=0;i1<NE1;i1++) {             node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1);
187         k=i1+faceNECount;  
188         node0=2*i1*N0;             out->Elements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1);
189               out->Elements->Tag[k]=0;
190         out->FaceElements->Id[k]=i1+totalNECount;             out->Elements->Owner[k]=myRank;
191         out->FaceElements->Tag[k]=1;  
192         out->FaceElements->Color[k]=(i1%2)+4;             out->Elements->Nodes[INDEX2(0,k,NN)]=node0;
193               out->Elements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride0;
194         if (useElementsOnFace) {             out->Elements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride1+2*Nstride0;
195            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0;             out->Elements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride1;
196            out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;             out->Elements->Nodes[INDEX2(4,k,NN)]=node0+1*Nstride0;
197            out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2;             out->Elements->Nodes[INDEX2(5,k,NN)]=node0+Nstride1+2*Nstride0;
198            out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0+2;             out->Elements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride1+1*Nstride0;
199            out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0;             out->Elements->Nodes[INDEX2(7,k,NN)]=node0+Nstride1;
200            out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;             if (useFullElementOrder) {
201            out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0+2;                out->Elements->Nodes[INDEX2(8,k,NN)]=node0+1*Nstride1+1*Nstride0;
202            out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0+1;             }
203         } else {           }
           out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0;  
           out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;  
           out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0;  
        }  
     
204       }       }
205       totalNECount+=NE1;       /* face elements */
206       faceNECount+=NE1;       NN=out->FaceElements->numNodes;
207           totalNECount=NE0*NE1;
208       /* **  elements on boundary 002 (x1=1): */       faceNECount=0;
209         if (!periodic[0]) {
210            /* **  elements on boundary 001 (x1=0): */
211            
212       #pragma omp parallel for private(i1,k,node0)          if (e_offset0 == 0) {
213       for (i1=0;i1<NE1;i1++) {             #pragma omp parallel for private(i1,k,node0)
214         k=i1+faceNECount;             for (i1=0;i1<local_NE1;i1++) {
215         node0=2*(NE0-1)+2*i1*N0;        
216                   k=i1+faceNECount;
217                   node0=Nstride1*N_PER_E*(i1+e_offset1);
218    
219                   out->FaceElements->Id[k]=i1+e_offset1+totalNECount;
220                   out->FaceElements->Tag[k]=1;
221                   out->FaceElements->Owner[k]=myRank;
222                   if (useElementsOnFace) {
223                      out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride1;
224                      out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0;
225                      out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride0;
226                      out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride1+2*Nstride0;
227                      out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+Nstride1;
228                      out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+1*Nstride0;
229                      out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+Nstride1+2*Nstride0;
230                      out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+2*Nstride1+1*Nstride0;
231                   } else {
232                      out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride1;
233                      out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0;
234                      out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride1;
235                   }
236               }
237               faceNECount+=local_NE1;
238            }
239            totalNECount+=NE1;
240            /* **  elements on boundary 002 (x1=1): */
241            if (local_NE0+e_offset0 == NE0) {
242               #pragma omp parallel for private(i1,k,node0)
243               for (i1=0;i1<local_NE1;i1++) {
244                   k=i1+faceNECount;
245                   node0=Nstride0*N_PER_E*(NE0-1)+Nstride1*N_PER_E*(i1+e_offset1);
246    
247                   out->FaceElements->Id[k]=(i1+e_offset1)+totalNECount;
248                   out->FaceElements->Tag[k]=2;
249                   out->FaceElements->Owner[k]=myRank;
250    
251                   if (useElementsOnFace) {
252                      out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride0;
253                      out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride1+2*Nstride0;
254                      out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride1;
255                      out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0;
256                      out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+Nstride1+2*Nstride0;
257                      out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+2*Nstride1+1*Nstride0;
258                      out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+Nstride1;
259                      out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+1*Nstride0;
260                   } else {
261                      out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride0;
262                      out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride1+2*Nstride0;
263                      out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride1+2*Nstride0;
264                   }
265               }
266               faceNECount+=local_NE1;
267             }
268             totalNECount+=NE1;
269         }
270         if (!periodic[1]) {
271            /* **  elements on boundary 010 (x2=0): */
272            if (e_offset1 == 0) {
273               #pragma omp parallel for private(i0,k,node0)
274               for (i0=0;i0<local_NE0;i0++) {
275                   k=i0+faceNECount;
276                   node0=Nstride0*N_PER_E*(i0+e_offset0);
277            
278                   out->FaceElements->Id[k]=e_offset0+i0+totalNECount;
279                   out->FaceElements->Tag[k]=10;
280                   out->FaceElements->Owner[k]=myRank;
281          
282                   if (useElementsOnFace) {
283                       out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0;
284                       out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride0;
285                       out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride1+2*Nstride0;
286                       out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride1;
287                       out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+1*Nstride0;
288                       out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+Nstride1+2*Nstride0;
289                       out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride1+1*Nstride0;
290                       out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+Nstride1;
291                   } else {
292                       out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0;
293                       out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride0;
294                       out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+1*Nstride0;
295                   }
296               }
297               faceNECount+=local_NE0;
298            }
299            totalNECount+=NE0;
300            /* **  elements on boundary 020 (x2=1): */
301            if (local_NE1+e_offset1 == NE1) {
302               #pragma omp parallel for private(i0,k,node0)
303               for (i0=0;i0<local_NE0;i0++) {
304                   k=i0+faceNECount;
305                   node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(NE1-1);
306        
307         out->FaceElements->Id[k]=i1+totalNECount;                 out->FaceElements->Id[k]=i0+e_offset0+totalNECount;
308         out->FaceElements->Tag[k]=2;                 out->FaceElements->Tag[k]=20;
309         out->FaceElements->Color[k]=(i1%2)+6;                 out->FaceElements->Owner[k]=myRank;
310                   if (useElementsOnFace) {
311                        out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride1+2*Nstride0;
312                        out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride1;
313                        out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0;
314                        out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride0;
315                        out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+2*Nstride1+1*Nstride0;
316                        out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+Nstride1;
317                        out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+1*Nstride0;
318                        out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+Nstride1+2*Nstride0;
319                   } else {
320                        out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride1+2*Nstride0;
321                        out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride1;
322                        out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride1+1*Nstride0;
323                   }
324               }
325               faceNECount+=local_NE0;
326            }
327            totalNECount+=NE0;
328         }
329         /* add tag names */
330         Finley_Mesh_addTagMap(out,"top", 20);
331         Finley_Mesh_addTagMap(out,"bottom", 10);
332         Finley_Mesh_addTagMap(out,"left", 1);
333         Finley_Mesh_addTagMap(out,"right", 2);
334        
335         if (useElementsOnFace) {       /* prepare mesh for further calculatuions:*/
336            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2;       if (Finley_noError()) {
337            out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0+2;           Finley_Mesh_resolveNodeIds(out);
338            out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0;       }
339            out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;       if (Finley_noError()) {
340            out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0+2;           Finley_Mesh_prepare(out, optimize);
           out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0+1;  
           out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0;  
           out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;  
        } else {  
           out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2;  
           out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0+2;  
           out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+2;  
        }  
341       }       }
      totalNECount+=NE1;  
      faceNECount+=NE1;  
   
342    }    }
343    out->FaceElements->minColor=0;    if (!Finley_noError()) {
344    out->FaceElements->maxColor=7;        Finley_Mesh_free(out);
345      }
346    /*  face elements done: */    /* free up memory */
347        Paso_MPIInfo_free( mpi_info );  
   /*   condense the nodes: */  
     
   Finley_Mesh_resolveNodeIds(out);  
   
   /* prepare mesh for further calculatuions:*/  
   
   Finley_Mesh_prepare(out) ;  
   
348    #ifdef Finley_TRACE    #ifdef Finley_TRACE
349    printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);    printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
350    #endif    #endif
351    
   if (! Finley_noError()) {  
       Finley_Mesh_dealloc(out);  
       return NULL;  
   }  
352    return out;    return out;
353  }  }

Legend:
Removed from v.471  
changed lines
  Added in v.1388

  ViewVC Help
Powered by ViewVC 1.1.26