finley/src/Mesh_hex8.c


/* $Id$ */

/*******************************************************
 *
 *           Copyright 2003-2007 by ACceSS MNRF
 *       Copyright 2007 by University of Queensland
 *
 *                http://esscc.uq.edu.au
 *        Primary Business: Queensland, Australia
 *  Licensed under the Open Software License version 3.0
 *     http://www.opensource.org/licenses/osl-3.0.php
 *
 *******************************************************/

/**************************************************************/

/*   Finley: generates rectangular meshes */

/*   Generates a numElements[0] x numElements[1] x numElements[2] mesh with first order elements (Hex8) in the brick */
/*   [0,Length[0]] x [0,Length[1]] x [0,Length[2]]. order is the desired accuracy of the */
/*   integration scheme. */

/**************************************************************/

#include "RectangularMesh.h"

Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,
                                          double* Length,
                                          bool_t* periodic,
                                          index_t order, 
                                          index_t reduced_order, 
                                          bool_t useElementsOnFace,
                                          bool_t useFullElementOrder,
                                          bool_t optimize)
{
  #define N_PER_E 1
  #define DIM 3
  dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,Nstride0,Nstride1,Nstride2, local_NE0, local_NE1, local_NE2, local_N0, local_N1, local_N2;
  dim_t totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements, NN;
  index_t node0, myRank, e_offset2, e_offset1, e_offset0, offset1, offset2, offset0, global_i0, global_i1, global_i2;
  Finley_Mesh* out;
  Paso_MPIInfo *mpi_info = NULL;
  char name[50];
  double time0=Finley_timer();

  /* get MPI information */
  mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
  if (! Finley_noError()) {
        return NULL;
  }
  myRank=mpi_info->rank;

  /* set up the global dimensions of the mesh */

  NE0=MAX(1,numElements[0]);
  NE1=MAX(1,numElements[1]);
  NE2=MAX(1,numElements[2]);
  N0=N_PER_E*NE0+1;
  N1=N_PER_E*NE1+1;
  N2=N_PER_E*NE2+1;

  /*  allocate mesh: */  
  sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
  out=Finley_Mesh_alloc(name,DIM,order, reduced_order, mpi_info);
  if (! Finley_noError()) { 
      Paso_MPIInfo_free( mpi_info );
      return NULL;
  }

 Finley_Mesh_setElements(out,Finley_ElementFile_alloc(Hex8,out->order,out->reduced_order,mpi_info));
 if (useElementsOnFace) {
         Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Hex8Face,
                                                                  out->order,
                                                                  out->reduced_order,
                                                                  mpi_info));
         Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Hex8Face_Contact,
                                                                    out->order,
                                                                    out->reduced_order,
                                                                    mpi_info));
  } else {
         Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Rec4,
                                                                  out->order,
                                                                  out->reduced_order,
                                                                  mpi_info));
         Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Rec4_Contact,
                                                                     out->order,
                                                                     out->reduced_order,
                                                                     mpi_info));
  }
  Finley_Mesh_setPoints(out,Finley_ElementFile_alloc(Point1,
                                                 out->order,
                                                 out->reduced_order,
                                                 mpi_info));
  if (! Finley_noError()) {
      Paso_MPIInfo_free( mpi_info );
      Finley_Mesh_free(out);
      return NULL;
  }

  /* work out the largest dimension */
  if (N2==MAX3(N0,N1,N2)) {
     Nstride0=1;
     Nstride1=N0;
     Nstride2=N0*N1;
     local_NE0=NE0;
     e_offset0=0;
     local_NE1=NE1;
     e_offset1=0;
     Paso_MPIInfo_Split(mpi_info,NE2,&local_NE2,&e_offset2);
  } else if (N1==MAX3(N0,N1,N2)) {
     Nstride0=N2;
     Nstride1=N0*N2;
     Nstride2=1;
     local_NE0=NE0;
     e_offset0=0;
     Paso_MPIInfo_Split(mpi_info,NE1,&local_NE1,&e_offset1);
     local_NE2=NE2;
     e_offset2=0;
  } else {
     Nstride0=N1*N2;
     Nstride1=1;
     Nstride2=N1;
     Paso_MPIInfo_Split(mpi_info,NE0,&local_NE0,&e_offset0);
     local_NE1=NE1;
     e_offset1=0;
     local_NE2=NE2;
     e_offset2=0;
  }
  offset0=e_offset0*N_PER_E;
  offset1=e_offset1*N_PER_E;
  offset2=e_offset2*N_PER_E;
  local_N0=local_NE0*N_PER_E+1;
  local_N1=local_NE1*N_PER_E+1;
  local_N2=local_NE2*N_PER_E+1;

  /* get the number of surface elements */

  NFaceElements=0;
  if (!periodic[2]) {
    NDOF2=N2;
    if (offset2==0) NFaceElements+=local_NE1*local_NE0;
    if (local_NE2+e_offset2 == NE2) NFaceElements+=local_NE1*local_NE0;
  } else {
      NDOF2=N2-1;
  }
 
  if (!periodic[0]) {
     NDOF0=N0;
     if (e_offset0 == 0) NFaceElements+=local_NE1*local_NE2;
     if (local_NE0+e_offset0 == NE0) NFaceElements+=local_NE1*local_NE2;
  } else {
      NDOF0=N0-1;
  }
  if (!periodic[1]) {
     NDOF1=N1;
     if (e_offset1 == 0) NFaceElements+=local_NE0*local_NE2;
     if (local_NE1+e_offset1 == NE1) NFaceElements+=local_NE0*local_NE2;
  } else {
      NDOF1=N1-1;
  }

  /*  allocate tables: */

  Finley_NodeFile_allocTable(out->Nodes,local_N0*local_N1*local_N2);
  Finley_ElementFile_allocTable(out->Elements,local_NE0*local_NE1*local_NE2);
  Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);

  if (Finley_noError()) {
     /* create nodes */
   
     #pragma omp parallel for private(i0,i1,i2,k,global_i0,global_i1,global_i2)
     for (i2=0;i2<local_N2;i2++) {
       for (i1=0;i1<local_N1;i1++) {
         for (i0=0;i0<local_N0;i0++) {
           k=i0+local_N0*i1+local_N0*local_N1*i2;
           global_i0=i0+offset0;
           global_i1=i1+offset1;
           global_i2=i2+offset2;
           out->Nodes->Coordinates[INDEX2(0,k,DIM)]=DBLE(global_i0)/DBLE(N0-1)*Length[0];
           out->Nodes->Coordinates[INDEX2(1,k,DIM)]=DBLE(global_i1)/DBLE(N1-1)*Length[1];
           out->Nodes->Coordinates[INDEX2(2,k,DIM)]=DBLE(global_i2)/DBLE(N2-1)*Length[2];
           out->Nodes->Id[k]=Nstride0*global_i0+Nstride1*global_i1+Nstride2*global_i2;
           out->Nodes->Tag[k]=0;
           out->Nodes->globalDegreesOfFreedom[k]=Nstride0*(global_i0%NDOF0) 
                                               +Nstride1*(global_i1%NDOF1) 
                                               +Nstride2*(global_i2%NDOF2);
         }
       }
     }
     /*   set the elements: */
     NN=out->Elements->numNodes;
     #pragma omp parallel for private(i0,i1,i2,k,node0) 
     for (i2=0;i2<local_NE2;i2++) {
       for (i1=0;i1<local_NE1;i1++) {
         for (i0=0;i0<local_NE0;i0++) {
           
           k=i0+local_NE0*i1+local_NE0*local_NE1*i2;        
           node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
   
           out->Elements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1)+NE0*NE1*(i2+e_offset2);
           out->Elements->Tag[k]=0;
           out->Elements->Owner[k]=myRank;

           out->Elements->Nodes[INDEX2(0,k,NN)] =node0                             ;
           out->Elements->Nodes[INDEX2(1,k,NN)] =node0+                    Nstride0;
           out->Elements->Nodes[INDEX2(2,k,NN)] =node0+           Nstride1+Nstride0;
           out->Elements->Nodes[INDEX2(3,k,NN)] =node0+           Nstride1;
           out->Elements->Nodes[INDEX2(4,k,NN)] =node0+Nstride2                   ;
           out->Elements->Nodes[INDEX2(5,k,NN)] =node0+Nstride2         +Nstride0;
           out->Elements->Nodes[INDEX2(6,k,NN)] =node0+Nstride2+Nstride1+Nstride0;
           out->Elements->Nodes[INDEX2(7,k,NN)] =node0+Nstride2+Nstride1           ;
         }
       }
     }
     /* face elements */
     NN=out->FaceElements->numNodes;
     totalNECount=NE0*NE1*NE2;
     faceNECount=0;
     /*   these are the quadrilateral elements on boundary 1 (x3=0): */
     if (!periodic[2]) {
       /* **  elements on boundary 100 (x3=0): */
       if (e_offset2==0) {
          #pragma omp parallel for private(i0,i1,k,node0) 
          for (i1=0;i1<local_NE1;i1++) {
            for (i0=0;i0<local_NE0;i0++) {
           
              k=i0+local_NE0*i1+faceNECount;
              node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1);
     
              out->FaceElements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1)+totalNECount;
              out->FaceElements->Tag[k]=100;
              out->FaceElements->Owner[k]=myRank;
        
              if  (useElementsOnFace) {
                 out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0                                 ;
                 out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0         +Nstride1           ;
                 out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0         +Nstride1+Nstride0;
                 out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+          Nstride0           ;
                 out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+Nstride2                      ;
                 out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+Nstride2+Nstride1           ;
                 out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+Nstride2+Nstride1+Nstride0;
                 out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+Nstride2          +Nstride0;
              } else {
                 out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0                                 ;
                 out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+           Nstride1           ;
                 out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+           Nstride1+Nstride0;
                 out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+                     Nstride0;
              }
            }
          }
          faceNECount+=local_NE1*local_NE0;
       }
       totalNECount+=NE1*NE0;
       /* **  elements on boundary 200 (x3=1) */
       if (local_NE2+e_offset2 == NE2) {
          #pragma omp parallel for private(i0,i1,k,node0) 
          for (i1=0;i1<local_NE1;i1++) {
            for (i0=0;i0<local_NE0;i0++) {
      
              k=i0+local_NE0*i1+faceNECount;
              node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(NE2-1);
        
              out->FaceElements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1)+totalNECount;
              out->FaceElements->Tag[k]=200;
              out->FaceElements->Owner[k]=myRank;
              if  (useElementsOnFace) {
                 out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0+Nstride2                      ;
                 out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+Nstride2+           Nstride0;
                 out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+Nstride2+Nstride1+Nstride0;
                 out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+Nstride2+Nstride1           ;
      
                 out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0                                 ;
                 out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+Nstride0                      ;
                 out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+           Nstride1+Nstride0;
                 out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+           Nstride1;
              } else {
                 out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0+Nstride2                      ;
                 out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+Nstride2           +Nstride0;
                 out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+Nstride2+Nstride1+Nstride0;
                 out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+Nstride2+Nstride1           ;
              }
            }
          }
          faceNECount+=local_NE1*local_NE0;
       }
       totalNECount+=NE1*NE0;
     }
     if (!periodic[0]) {
        /* **  elements on boundary 001 (x1=0): */
     
        if (e_offset0 == 0) {
           #pragma omp parallel for private(i1,i2,k,node0) 
           for (i2=0;i2<local_NE2;i2++) {
             for (i1=0;i1<local_NE1;i1++) {
      
               k=i1+local_NE1*i2+faceNECount;
               node0=Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
               out->FaceElements->Id[k]=(i1+e_offset1)+NE1*(i2+e_offset2)+totalNECount;
               out->FaceElements->Tag[k]=1;
               out->FaceElements->Owner[k]=myRank;
      
               if  (useElementsOnFace) {
                  out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0                                 ;
                  out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+Nstride2                      ;
                  out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+Nstride2+Nstride1           ;
                  out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+Nstride1                      ;
                  out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+Nstride0                      ;
                  out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+Nstride2+Nstride0           ;
                  out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+Nstride2+Nstride1+Nstride0;
                  out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+Nstride1+Nstride0           ;
               } else {
                  out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0                                 ;
                  out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+Nstride2                      ;
                  out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+Nstride2+Nstride1           ;
                  out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+           Nstride1           ;
               }
             }
           }
           faceNECount+=local_NE1*local_NE2;
        }
        totalNECount+=NE1*NE2;
        /* **  elements on boundary 002 (x1=1): */
        if (local_NE0+e_offset0 == NE0) {
           #pragma omp parallel for private(i1,i2,k,node0) 
           for (i2=0;i2<local_NE2;i2++) {
             for (i1=0;i1<local_NE1;i1++) {
               k=i1+local_NE1*i2+faceNECount;
               node0=Nstride0*N_PER_E*(NE0-1)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
               out->FaceElements->Id[k]=(i1+e_offset1)+NE1*(i2+e_offset2)+totalNECount;
               out->FaceElements->Tag[k]=2;
               out->FaceElements->Owner[k]=myRank;
   
               if  (useElementsOnFace) {
                  out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+                      Nstride0;
                  out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+           Nstride1+Nstride0;
                  out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2+Nstride1+Nstride0;
                  out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride2+           Nstride0;
      
                  out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0                                 ;
                  out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+           Nstride1           ;
                  out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+Nstride2+Nstride1           ;
                  out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+Nstride2                      ;
      
               } else {
                  out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0                      +Nstride0;
                  out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+           Nstride1+Nstride0;
                  out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2+Nstride1+Nstride0;
                  out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride2+           Nstride0;
               }
         
             }
           }
           faceNECount+=local_NE1*local_NE2;
         }
         totalNECount+=NE1*NE2;
     }
     if (!periodic[1]) {
        /* **  elements on boundary 010 (x2=0): */
        if (e_offset1 == 0) {
           #pragma omp parallel for private(i0,i2,k,node0) 
           for (i2=0;i2<local_NE2;i2++) {
             for (i0=0;i0<local_NE0;i0++) {
               k=i0+local_NE0*i2+faceNECount;
               node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride2*N_PER_E*(i2+e_offset2);
         
               out->FaceElements->Id[k]=(i2+e_offset2)+NE2*(e_offset0+i0)+totalNECount;
               out->FaceElements->Tag[k]=10;
               out->FaceElements->Owner[k]=myRank;
               if  (useElementsOnFace) {
                  out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0                                 ;
                  out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+                      Nstride0;
                  out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2           +Nstride0;
                  out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride2                      ;
      
                  out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+           Nstride1           ;
                  out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+Nstride1+           Nstride0;
                  out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+Nstride2+Nstride1+Nstride0;
                  out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+Nstride2+Nstride1           ;
               } else {
                  out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0                                 ;
                  out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+                      Nstride0;
                  out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2+           Nstride0;
                  out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride2                      ;
               }
             }
           }
           faceNECount+=local_NE0*local_NE2;
        }
        totalNECount+=NE0*NE2;
        /* **  elements on boundary 020 (x2=1): */
        if (local_NE1+e_offset1 == NE1) {
           #pragma omp parallel for private(i0,i2,k,node0) 
           for (i2=0;i2<local_NE2;i2++) {
             for (i0=0;i0<local_NE0;i0++) {
               k=i0+local_NE0*i2+faceNECount;
               node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(NE1-1)+Nstride2*N_PER_E*(i2+e_offset2);
   
               out->FaceElements->Id[k]=(i2+e_offset2)+NE2*(i0+e_offset0)+totalNECount;
               out->FaceElements->Tag[k]=20;
               out->FaceElements->Owner[k]=myRank;
      
               if  (useElementsOnFace) {
                  out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+           Nstride1           ;
                  out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride2+Nstride1           ;
                  out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2+Nstride1+Nstride0;
                  out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride1+Nstride0           ;
      
                  out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0                                 ;
                  out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+Nstride2                      ;
                  out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+Nstride2+           Nstride0;
                  out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+                      Nstride0;
               } else {
                  out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+           Nstride1           ;
                  out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride2+Nstride1           ;
                  out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2+Nstride1+Nstride0;
                  out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+           Nstride1+Nstride0;
               }
             }
           }
           faceNECount+=local_NE0*local_NE2;
        }
        totalNECount+=NE0*NE2;
     }
     /* add tag names */
     Finley_Mesh_addTagMap(out,"top", 200);
     Finley_Mesh_addTagMap(out,"bottom", 100);
     Finley_Mesh_addTagMap(out,"left", 1);
     Finley_Mesh_addTagMap(out,"right", 2);
     Finley_Mesh_addTagMap(out,"front", 10);
     Finley_Mesh_addTagMap(out,"back", 20);
   
     /* prepare mesh for further calculatuions:*/
     if (Finley_noError()) {
         Finley_Mesh_resolveNodeIds(out);
     }
     if (Finley_noError()) {
         Finley_Mesh_prepare(out, optimize);
     }
  }

  if (!Finley_noError()) {
      Finley_Mesh_free(out);
  }
  /* free up memory */
  Paso_MPIInfo_free( mpi_info );  
  #ifdef Finley_TRACE
  printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
  #endif

  return out;
}
1	jgs	150
2	ksteube	1312	/* $Id$ */
3
4			/*******************************************************
5			*
6			* Copyright 2003-2007 by ACceSS MNRF
7			* Copyright 2007 by University of Queensland
8			*
9			* http://esscc.uq.edu.au
10			* 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	jgs	82	/**************************************************************/
17
18			/* Finley: generates rectangular meshes */
19
20			/* Generates a numElements[0] x numElements[1] x numElements[2] mesh with first order elements (Hex8) in the brick */
21			/* [0,Length[0]] x [0,Length[1]] x [0,Length[2]]. order is the desired accuracy of the */
22			/* integration scheme. */
23
24			/**************************************************************/
25
26			#include "RectangularMesh.h"
27
28	ksteube	1312	Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,
29			double* Length,
30			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	bcumming	751	{
37	ksteube	1312	#define N_PER_E 1
38			#define DIM 3
39			dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,Nstride0,Nstride1,Nstride2, local_NE0, local_NE1, local_NE2, local_N0, local_N1, local_N2;
40			dim_t totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements, NN;
41			index_t node0, myRank, e_offset2, e_offset1, e_offset0, offset1, offset2, offset0, global_i0, global_i1, global_i2;
42			Finley_Mesh* out;
43			Paso_MPIInfo *mpi_info = NULL;
44			char name[50];
45			double time0=Finley_timer();
46	bcumming	751
47	ksteube	1312	/* get MPI information */
48			mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
49			if (! Finley_noError()) {
50			return NULL;
51	bcumming	751	}
52	ksteube	1312	myRank=mpi_info->rank;
53	bcumming	751
54	ksteube	1312	/* set up the global dimensions of the mesh */
55	bcumming	751
56	jgs	82	NE0=MAX(1,numElements[0]);
57			NE1=MAX(1,numElements[1]);
58			NE2=MAX(1,numElements[2]);
59	ksteube	1312	N0=N_PER_E*NE0+1;
60			N1=N_PER_E*NE1+1;
61			N2=N_PER_E*NE2+1;
62	jgs	82
63	ksteube	1312	/* allocate mesh: */
64			sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
65			out=Finley_Mesh_alloc(name,DIM,order, reduced_order, mpi_info);
66			if (! Finley_noError()) {
67			Paso_MPIInfo_free( mpi_info );
68			return NULL;
69	jgs	153	}
70
71	ksteube	1312	Finley_Mesh_setElements(out,Finley_ElementFile_alloc(Hex8,out->order,out->reduced_order,mpi_info));
72			if (useElementsOnFace) {
73			Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Hex8Face,
74			out->order,
75			out->reduced_order,
76			mpi_info));
77			Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Hex8Face_Contact,
78			out->order,
79			out->reduced_order,
80			mpi_info));
81	jgs	82	} else {
82	ksteube	1312	Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Rec4,
83			out->order,
84			out->reduced_order,
85			mpi_info));
86			Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Rec4_Contact,
87			out->order,
88			out->reduced_order,
89			mpi_info));
90	jgs	82	}
91	ksteube	1312	Finley_Mesh_setPoints(out,Finley_ElementFile_alloc(Point1,
92			out->order,
93			out->reduced_order,
94			mpi_info));
95			if (! Finley_noError()) {
96			Paso_MPIInfo_free( mpi_info );
97			Finley_Mesh_free(out);
98			return NULL;
99	jgs	82	}
100	ksteube	1312
101			/* work out the largest dimension */
102			if (N2==MAX3(N0,N1,N2)) {
103			Nstride0=1;
104			Nstride1=N0;
105			Nstride2=N0*N1;
106			local_NE0=NE0;
107			e_offset0=0;
108			local_NE1=NE1;
109			e_offset1=0;
110			Paso_MPIInfo_Split(mpi_info,NE2,&local_NE2,&e_offset2);
111			} else if (N1==MAX3(N0,N1,N2)) {
112			Nstride0=N2;
113			Nstride1=N0*N2;
114			Nstride2=1;
115			local_NE0=NE0;
116			e_offset0=0;
117			Paso_MPIInfo_Split(mpi_info,NE1,&local_NE1,&e_offset1);
118			local_NE2=NE2;
119			e_offset2=0;
120	jgs	82	} else {
121	ksteube	1312	Nstride0=N1*N2;
122			Nstride1=1;
123			Nstride2=N1;
124			Paso_MPIInfo_Split(mpi_info,NE0,&local_NE0,&e_offset0);
125			local_NE1=NE1;
126			e_offset1=0;
127			local_NE2=NE2;
128			e_offset2=0;
129	jgs	82	}
130	ksteube	1312	offset0=e_offset0*N_PER_E;
131			offset1=e_offset1*N_PER_E;
132			offset2=e_offset2*N_PER_E;
133			local_N0=local_NE0*N_PER_E+1;
134			local_N1=local_NE1*N_PER_E+1;
135			local_N2=local_NE2*N_PER_E+1;
136	bcumming	751
137	ksteube	1312	/* get the number of surface elements */
138	jgs	82
139	ksteube	1312	NFaceElements=0;
140			if (!periodic[2]) {
141			NDOF2=N2;
142			if (offset2==0) NFaceElements+=local_NE1*local_NE0;
143			if (local_NE2+e_offset2 == NE2) NFaceElements+=local_NE1*local_NE0;
144	bcumming	751	} else {
145	ksteube	1312	NDOF2=N2-1;
146	bcumming	751	}
147	ksteube	1312
148			if (!periodic[0]) {
149			NDOF0=N0;
150			if (e_offset0 == 0) NFaceElements+=local_NE1*local_NE2;
151			if (local_NE0+e_offset0 == NE0) NFaceElements+=local_NE1*local_NE2;
152	jgs	82	} else {
153	ksteube	1312	NDOF0=N0-1;
154	jgs	82	}
155			if (!periodic[1]) {
156	ksteube	1312	NDOF1=N1;
157			if (e_offset1 == 0) NFaceElements+=local_NE0*local_NE2;
158			if (local_NE1+e_offset1 == NE1) NFaceElements+=local_NE0*local_NE2;
159			} else {
160			NDOF1=N1-1;
161	jgs	82	}
162	jgs	153
163	ksteube	1312	/* allocate tables: */
164	jgs	82
165	ksteube	1312	Finley_NodeFile_allocTable(out->Nodes,local_N0local_N1local_N2);
166			Finley_ElementFile_allocTable(out->Elements,local_NE0local_NE1local_NE2);
167			Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
168	jgs	82
169	ksteube	1312	if (Finley_noError()) {
170			/* create nodes */
171	jgs	82
172	ksteube	1312	#pragma omp parallel for private(i0,i1,i2,k,global_i0,global_i1,global_i2)
173			for (i2=0;i2<local_N2;i2++) {
174			for (i1=0;i1<local_N1;i1++) {
175			for (i0=0;i0<local_N0;i0++) {
176			k=i0+local_N0i1+local_N0local_N1*i2;
177			global_i0=i0+offset0;
178			global_i1=i1+offset1;
179			global_i2=i2+offset2;
180			out->Nodes->Coordinates[INDEX2(0,k,DIM)]=DBLE(global_i0)/DBLE(N0-1)*Length[0];
181			out->Nodes->Coordinates[INDEX2(1,k,DIM)]=DBLE(global_i1)/DBLE(N1-1)*Length[1];
182			out->Nodes->Coordinates[INDEX2(2,k,DIM)]=DBLE(global_i2)/DBLE(N2-1)*Length[2];
183			out->Nodes->Id[k]=Nstride0global_i0+Nstride1global_i1+Nstride2*global_i2;
184			out->Nodes->Tag[k]=0;
185			out->Nodes->globalDegreesOfFreedom[k]=Nstride0*(global_i0%NDOF0)
186			+Nstride1*(global_i1%NDOF1)
187			+Nstride2*(global_i2%NDOF2);
188	jgs	82	}
189			}
190			}
191	ksteube	1312	/* set the elements: */
192			NN=out->Elements->numNodes;
193			#pragma omp parallel for private(i0,i1,i2,k,node0)
194			for (i2=0;i2<local_NE2;i2++) {
195			for (i1=0;i1<local_NE1;i1++) {
196			for (i0=0;i0<local_NE0;i0++) {
197
198			k=i0+local_NE0i1+local_NE0local_NE1*i2;
199			node0=Nstride0N_PER_E(i0+e_offset0)+Nstride1N_PER_E(i1+e_offset1)+Nstride2N_PER_E(i2+e_offset2);
200	jgs	82
201	ksteube	1312	out->Elements->Id[k]=(i0+e_offset0)+NE0(i1+e_offset1)+NE0NE1*(i2+e_offset2);
202			out->Elements->Tag[k]=0;
203			out->Elements->Owner[k]=myRank;
204	jgs	82
205	ksteube	1312	out->Elements->Nodes[INDEX2(0,k,NN)] =node0 ;
206			out->Elements->Nodes[INDEX2(1,k,NN)] =node0+ Nstride0;
207			out->Elements->Nodes[INDEX2(2,k,NN)] =node0+ Nstride1+Nstride0;
208			out->Elements->Nodes[INDEX2(3,k,NN)] =node0+ Nstride1;
209			out->Elements->Nodes[INDEX2(4,k,NN)] =node0+Nstride2 ;
210			out->Elements->Nodes[INDEX2(5,k,NN)] =node0+Nstride2 +Nstride0;
211			out->Elements->Nodes[INDEX2(6,k,NN)] =node0+Nstride2+Nstride1+Nstride0;
212			out->Elements->Nodes[INDEX2(7,k,NN)] =node0+Nstride2+Nstride1 ;
213	jgs	82	}
214			}
215			}
216	ksteube	1312	/* face elements */
217			NN=out->FaceElements->numNodes;
218			totalNECount=NE0NE1NE2;
219			faceNECount=0;
220			/* these are the quadrilateral elements on boundary 1 (x3=0): */
221			if (!periodic[2]) {
222			/* ** elements on boundary 100 (x3=0): */
223			if (e_offset2==0) {
224			#pragma omp parallel for private(i0,i1,k,node0)
225			for (i1=0;i1<local_NE1;i1++) {
226			for (i0=0;i0<local_NE0;i0++) {
227
228			k=i0+local_NE0*i1+faceNECount;
229			node0=Nstride0N_PER_E(i0+e_offset0)+Nstride1N_PER_E(i1+e_offset1);
230	jgs	82
231	ksteube	1312	out->FaceElements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1)+totalNECount;
232			out->FaceElements->Tag[k]=100;
233			out->FaceElements->Owner[k]=myRank;
234
235			if (useElementsOnFace) {
236			out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0 ;
237			out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0 +Nstride1 ;
238			out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0 +Nstride1+Nstride0;
239			out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+ Nstride0 ;
240			out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+Nstride2 ;
241			out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+Nstride2+Nstride1 ;
242			out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+Nstride2+Nstride1+Nstride0;
243			out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+Nstride2 +Nstride0;
244			} else {
245			out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0 ;
246			out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+ Nstride1 ;
247			out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+ Nstride1+Nstride0;
248			out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+ Nstride0;
249			}
250			}
251			}
252			faceNECount+=local_NE1*local_NE0;
253	jgs	82	}
254	ksteube	1312	totalNECount+=NE1*NE0;
255			/* ** elements on boundary 200 (x3=1) */
256			if (local_NE2+e_offset2 == NE2) {
257			#pragma omp parallel for private(i0,i1,k,node0)
258			for (i1=0;i1<local_NE1;i1++) {
259			for (i0=0;i0<local_NE0;i0++) {
260
261			k=i0+local_NE0*i1+faceNECount;
262			node0=Nstride0N_PER_E(i0+e_offset0)+Nstride1N_PER_E(i1+e_offset1)+Nstride2N_PER_E(NE2-1);
263
264			out->FaceElements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1)+totalNECount;
265			out->FaceElements->Tag[k]=200;
266			out->FaceElements->Owner[k]=myRank;
267			if (useElementsOnFace) {
268			out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0+Nstride2 ;
269			out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+Nstride2+ Nstride0;
270			out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+Nstride2+Nstride1+Nstride0;
271			out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+Nstride2+Nstride1 ;
272
273			out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0 ;
274			out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+Nstride0 ;
275			out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+ Nstride1+Nstride0;
276			out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+ Nstride1;
277			} else {
278			out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0+Nstride2 ;
279			out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+Nstride2 +Nstride0;
280			out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+Nstride2+Nstride1+Nstride0;
281			out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+Nstride2+Nstride1 ;
282			}
283			}
284			}
285			faceNECount+=local_NE1*local_NE0;
286	jgs	82	}
287	ksteube	1312	totalNECount+=NE1*NE0;
288	jgs	82	}
289	ksteube	1312	if (!periodic[0]) {
290			/* ** elements on boundary 001 (x1=0): */
291
292			if (e_offset0 == 0) {
293			#pragma omp parallel for private(i1,i2,k,node0)
294			for (i2=0;i2<local_NE2;i2++) {
295			for (i1=0;i1<local_NE1;i1++) {
296
297			k=i1+local_NE1*i2+faceNECount;
298			node0=Nstride1N_PER_E(i1+e_offset1)+Nstride2N_PER_E(i2+e_offset2);
299			out->FaceElements->Id[k]=(i1+e_offset1)+NE1*(i2+e_offset2)+totalNECount;
300			out->FaceElements->Tag[k]=1;
301			out->FaceElements->Owner[k]=myRank;
302
303			if (useElementsOnFace) {
304			out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0 ;
305			out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+Nstride2 ;
306			out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+Nstride2+Nstride1 ;
307			out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+Nstride1 ;
308			out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+Nstride0 ;
309			out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+Nstride2+Nstride0 ;
310			out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+Nstride2+Nstride1+Nstride0;
311			out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+Nstride1+Nstride0 ;
312			} else {
313			out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0 ;
314			out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+Nstride2 ;
315			out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+Nstride2+Nstride1 ;
316			out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+ Nstride1 ;
317			}
318			}
319			}
320			faceNECount+=local_NE1*local_NE2;
321			}
322			totalNECount+=NE1*NE2;
323			/* ** elements on boundary 002 (x1=1): */
324			if (local_NE0+e_offset0 == NE0) {
325			#pragma omp parallel for private(i1,i2,k,node0)
326			for (i2=0;i2<local_NE2;i2++) {
327			for (i1=0;i1<local_NE1;i1++) {
328			k=i1+local_NE1*i2+faceNECount;
329			node0=Nstride0N_PER_E(NE0-1)+Nstride1N_PER_E(i1+e_offset1)+Nstride2N_PER_E(i2+e_offset2);
330			out->FaceElements->Id[k]=(i1+e_offset1)+NE1*(i2+e_offset2)+totalNECount;
331			out->FaceElements->Tag[k]=2;
332			out->FaceElements->Owner[k]=myRank;
333	jgs	82
334	ksteube	1312	if (useElementsOnFace) {
335			out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+ Nstride0;
336			out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+ Nstride1+Nstride0;
337			out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2+Nstride1+Nstride0;
338			out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride2+ Nstride0;
339
340			out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0 ;
341			out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+ Nstride1 ;
342			out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+Nstride2+Nstride1 ;
343			out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+Nstride2 ;
344
345			} else {
346			out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0 +Nstride0;
347			out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+ Nstride1+Nstride0;
348			out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2+Nstride1+Nstride0;
349			out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride2+ Nstride0;
350			}
351
352			}
353			}
354			faceNECount+=local_NE1*local_NE2;
355	jgs	82	}
356	ksteube	1312	totalNECount+=NE1*NE2;
357	jgs	82	}
358	ksteube	1312	if (!periodic[1]) {
359			/* ** elements on boundary 010 (x2=0): */
360			if (e_offset1 == 0) {
361			#pragma omp parallel for private(i0,i2,k,node0)
362			for (i2=0;i2<local_NE2;i2++) {
363			for (i0=0;i0<local_NE0;i0++) {
364			k=i0+local_NE0*i2+faceNECount;
365			node0=Nstride0N_PER_E(i0+e_offset0)+Nstride2N_PER_E(i2+e_offset2);
366
367			out->FaceElements->Id[k]=(i2+e_offset2)+NE2*(e_offset0+i0)+totalNECount;
368			out->FaceElements->Tag[k]=10;
369			out->FaceElements->Owner[k]=myRank;
370			if (useElementsOnFace) {
371			out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0 ;
372			out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+ Nstride0;
373			out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2 +Nstride0;
374			out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride2 ;
375
376			out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+ Nstride1 ;
377			out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+Nstride1+ Nstride0;
378			out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+Nstride2+Nstride1+Nstride0;
379			out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+Nstride2+Nstride1 ;
380			} else {
381			out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0 ;
382			out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+ Nstride0;
383			out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2+ Nstride0;
384			out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride2 ;
385			}
386			}
387			}
388			faceNECount+=local_NE0*local_NE2;
389			}
390			totalNECount+=NE0*NE2;
391			/* ** elements on boundary 020 (x2=1): */
392			if (local_NE1+e_offset1 == NE1) {
393			#pragma omp parallel for private(i0,i2,k,node0)
394			for (i2=0;i2<local_NE2;i2++) {
395			for (i0=0;i0<local_NE0;i0++) {
396			k=i0+local_NE0*i2+faceNECount;
397			node0=Nstride0N_PER_E(i0+e_offset0)+Nstride1N_PER_E(NE1-1)+Nstride2N_PER_E(i2+e_offset2);
398	jgs	82
399	ksteube	1312	out->FaceElements->Id[k]=(i2+e_offset2)+NE2*(i0+e_offset0)+totalNECount;
400			out->FaceElements->Tag[k]=20;
401			out->FaceElements->Owner[k]=myRank;
402
403			if (useElementsOnFace) {
404			out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+ Nstride1 ;
405			out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride2+Nstride1 ;
406			out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2+Nstride1+Nstride0;
407			out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride1+Nstride0 ;
408
409			out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0 ;
410			out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+Nstride2 ;
411			out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+Nstride2+ Nstride0;
412			out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+ Nstride0;
413			} else {
414			out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+ Nstride1 ;
415			out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride2+Nstride1 ;
416			out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2+Nstride1+Nstride0;
417			out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+ Nstride1+Nstride0;
418			}
419			}
420			}
421			faceNECount+=local_NE0*local_NE2;
422			}
423			totalNECount+=NE0*NE2;
424	bcumming	751	}
425	ksteube	1312	/* add tag names */
426			Finley_Mesh_addTagMap(out,"top", 200);
427			Finley_Mesh_addTagMap(out,"bottom", 100);
428			Finley_Mesh_addTagMap(out,"left", 1);
429			Finley_Mesh_addTagMap(out,"right", 2);
430			Finley_Mesh_addTagMap(out,"front", 10);
431			Finley_Mesh_addTagMap(out,"back", 20);
432	bcumming	751
433	ksteube	1312	/* prepare mesh for further calculatuions:*/
434			if (Finley_noError()) {
435			Finley_Mesh_resolveNodeIds(out);
436	bcumming	751	}
437	ksteube	1312	if (Finley_noError()) {
438			Finley_Mesh_prepare(out, optimize);
439	bcumming	751	}
440			}
441
442	ksteube	1312	if (!Finley_noError()) {
443			Finley_Mesh_free(out);
444	bcumming	751	}
445	bcumming	782	/* free up memory */
446	ksteube	1312	Paso_MPIInfo_free( mpi_info );
447	bcumming	751	#ifdef Finley_TRACE
448			printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
449			#endif
450	ksteube	1312
451			return out;
452	bcumming	751	}
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision