finley/src/Mesh_rec8.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] mesh with second order elements (Rec8) in the rectangle */
/*   [0,Length[0]] x [0,Length[1]]. order is the desired accuracy of the integration scheme. */

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

#include "RectangularMesh.h"


Finley_Mesh* Finley_RectangularMesh_Rec8(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 2
  #define DIM 2
  dim_t N0,N1,NE0,NE1,i0,i1,k,Nstride0,Nstride1;
  dim_t totalNECount,faceNECount,NDOF0,NDOF1,NFaceElements,NN, local_NE0, local_NE1, local_N0, local_N1;
  index_t e_offset1, e_offset0, offset0, offset1, global_i0, global_i1;
  index_t node0, myRank;
  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]);
  N0=N_PER_E*NE0+1;
  N1=N_PER_E*NE1+1;

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

  if (useFullElementOrder) {
     Finley_Mesh_setElements(out,Finley_ElementFile_alloc(Rec9,
                                            out->order,
                                            out->reduced_order,
                                            mpi_info));
     if (useElementsOnFace) {
         Finley_setError(SYSTEM_ERROR,"rich elements for Rec9 elements is not supported yet.");
     } else {
         Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Line3,
                                                    out->order,
                                                    out->reduced_order,
                                                    mpi_info));
         Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Line3_Contact,
                                                       out->order,
                                                       out->reduced_order,
                                                       mpi_info));
     }
  } else  {
     Finley_Mesh_setElements(out,Finley_ElementFile_alloc(Rec8,out->order,out->reduced_order,mpi_info));
     if (useElementsOnFace) {
         Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Rec8Face,
                                                                  out->order,
                                                                  out->reduced_order,
                                                                  mpi_info));
         Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Rec8Face_Contact,
                                                                    out->order,
                                                                    out->reduced_order,
                                                                    mpi_info));
     } else {
         Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Line3,
                                                                  out->order,
                                                                  out->reduced_order,
                                                                  mpi_info));
         Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Line3_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 (N1==MAX(N0,N1)) {
     Nstride0=1;
     Nstride1=N0;
     local_NE0=NE0;
     e_offset0=0;
     Paso_MPIInfo_Split(mpi_info,NE1,&local_NE1,&e_offset1);
  } else {
     Nstride0=N1;
     Nstride1=1;
     Paso_MPIInfo_Split(mpi_info,NE0,&local_NE0,&e_offset0);
     local_NE1=NE1;
     e_offset1=0;
  }
  offset0=e_offset0*N_PER_E;
  offset1=e_offset1*N_PER_E;
  local_N0=local_NE0>0 ? local_NE0*N_PER_E+1 : 0;
  local_N1=local_NE1>0 ? local_NE1*N_PER_E+1 : 0;

  /* get the number of surface elements */

  NFaceElements=0;
  if (!periodic[0] &&  (local_NE0>0)) {
     NDOF0=N0;
     if (e_offset0 == 0) NFaceElements+=local_NE1;
     if (local_NE0+e_offset0 == NE0) NFaceElements+=local_NE1;
  } else {
      NDOF0=N0-1;
  }
  if (!periodic[1] && (local_NE1>0)) {
     NDOF1=N1;
     if (e_offset1 == 0) NFaceElements+=local_NE0;
     if (local_NE1+e_offset1 == NE1) NFaceElements+=local_NE0;
  } else {
      NDOF1=N1-1;
  }

  /*  allocate tables: */

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

  if (Finley_noError()) {
     /* create nodes */
   
     #pragma omp parallel for private(i0,i1,k,global_i0,global_i1)
     for (i1=0;i1<local_N1;i1++) {
       for (i0=0;i0<local_N0;i0++) {
           k=i0+local_N0*i1;
           global_i0=i0+offset0;
           global_i1=i1+offset1;
           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->Id[k]=Nstride0*global_i0+Nstride1*global_i1;
           out->Nodes->Tag[k]=0;
           out->Nodes->globalDegreesOfFreedom[k]=Nstride0*(global_i0%NDOF0) 
                                               +Nstride1*(global_i1%NDOF1);
       }
     }
     /*   set the elements: */
     NN=out->Elements->numNodes;
     #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;        
           node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1);

           out->Elements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1);
           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+2*Nstride0;
           out->Elements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride1+2*Nstride0;
           out->Elements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride1;
           out->Elements->Nodes[INDEX2(4,k,NN)]=node0+1*Nstride0;
           out->Elements->Nodes[INDEX2(5,k,NN)]=node0+Nstride1+2*Nstride0;
           out->Elements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride1+1*Nstride0;
           out->Elements->Nodes[INDEX2(7,k,NN)]=node0+Nstride1;
           if (useFullElementOrder) {
              out->Elements->Nodes[INDEX2(8,k,NN)]=node0+1*Nstride1+1*Nstride0;
           }
         }
     }
     /* face elements */
     NN=out->FaceElements->numNodes;
     totalNECount=NE0*NE1;
     faceNECount=0;
     if (!periodic[0] && (local_NE0>0)) {
        /* **  elements on boundary 001 (x1=0): */
     
        if (e_offset0 == 0) {
           #pragma omp parallel for private(i1,k,node0) 
           for (i1=0;i1<local_NE1;i1++) {
      
               k=i1+faceNECount;
               node0=Nstride1*N_PER_E*(i1+e_offset1);

               out->FaceElements->Id[k]=i1+e_offset1+totalNECount;
               out->FaceElements->Tag[k]=1;
               out->FaceElements->Owner[k]=myRank;
               if (useElementsOnFace) {
                  out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride1;
                  out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0;
                  out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride0;
                  out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride1+2*Nstride0;
                  out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+Nstride1;
                  out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+1*Nstride0;
                  out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+Nstride1+2*Nstride0;
                  out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+2*Nstride1+1*Nstride0;
               } else {
                  out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride1;
                  out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0;
                  out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride1;
               }
           }
           faceNECount+=local_NE1;
        }
        totalNECount+=NE1;
        /* **  elements on boundary 002 (x1=1): */
        if (local_NE0+e_offset0 == NE0) {
           #pragma omp parallel for private(i1,k,node0) 
           for (i1=0;i1<local_NE1;i1++) {
               k=i1+faceNECount;
               node0=Nstride0*N_PER_E*(NE0-1)+Nstride1*N_PER_E*(i1+e_offset1);

               out->FaceElements->Id[k]=(i1+e_offset1)+totalNECount;
               out->FaceElements->Tag[k]=2;
               out->FaceElements->Owner[k]=myRank;

               if (useElementsOnFace) {
                  out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride0;
                  out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride1+2*Nstride0;
                  out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride1;
                  out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0;
                  out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+Nstride1+2*Nstride0;
                  out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+2*Nstride1+1*Nstride0;
                  out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+Nstride1;
                  out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+1*Nstride0;
               } else {
                  out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride0;
                  out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride1+2*Nstride0;
                  out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride1+2*Nstride0;
               }
           }
           faceNECount+=local_NE1;
         }
         totalNECount+=NE1;
     }
     if (!periodic[1] && (local_NE0>0)) {
        /* **  elements on boundary 010 (x2=0): */
        if (e_offset1 == 0) {
           #pragma omp parallel for private(i0,k,node0) 
           for (i0=0;i0<local_NE0;i0++) {
               k=i0+faceNECount;
               node0=Nstride0*N_PER_E*(i0+e_offset0);
         
               out->FaceElements->Id[k]=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+2*Nstride0;
                   out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride1+2*Nstride0;
                   out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride1;
                   out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+1*Nstride0;
                   out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+Nstride1+2*Nstride0;
                   out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride1+1*Nstride0;
                   out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+Nstride1;
               } else {
                   out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0;
                   out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride0;
                   out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+1*Nstride0;
               }
           }
           faceNECount+=local_NE0;
        }
        totalNECount+=NE0;
        /* **  elements on boundary 020 (x2=1): */
        if (local_NE1+e_offset1 == NE1) {
           #pragma omp parallel for private(i0,k,node0) 
           for (i0=0;i0<local_NE0;i0++) {
               k=i0+faceNECount;
               node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(NE1-1);
   
               out->FaceElements->Id[k]=i0+e_offset0+totalNECount;
               out->FaceElements->Tag[k]=20;
               out->FaceElements->Owner[k]=myRank;
               if (useElementsOnFace) {
                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride1+2*Nstride0;
                    out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride1;
                    out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0;
                    out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride0;
                    out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+2*Nstride1+1*Nstride0;
                    out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+Nstride1;
                    out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+1*Nstride0;
                    out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+Nstride1+2*Nstride0;
               } else {
                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride1+2*Nstride0;
                    out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride1;
                    out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride1+1*Nstride0;
               }
           }
           faceNECount+=local_NE0;
        }
        totalNECount+=NE0;
     }
     /* add tag names */
     Finley_Mesh_addTagMap(out,"top", 20);
     Finley_Mesh_addTagMap(out,"bottom", 10);
     Finley_Mesh_addTagMap(out,"left", 1);
     Finley_Mesh_addTagMap(out,"right", 2);
   
     /* 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
2	/* $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	/**************************************************************/
17
18	/* Finley: generates rectangular meshes */
19
20	/* Generates a numElements[0] x numElements[1] mesh with second order elements (Rec8) in the rectangle */
21	/* [0,Length[0]] x [0,Length[1]]. order is the desired accuracy of the integration scheme. */
22
23	/**************************************************************/
24
25	#include "RectangularMesh.h"
26
27
28	Finley_Mesh* Finley_RectangularMesh_Rec8(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	{
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;
44	Paso_MPIInfo *mpi_info = NULL;
45	char name[50];
46	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]);
58	NE1=MAX(1,numElements[1]);
59	N0=N_PER_E*NE0+1;
60	N1=N_PER_E*NE1+1;
61
62	/* allocate mesh: */
63	sprintf(name,"Rectangular %d x %d mesh",N0,N1);
64	out=Finley_Mesh_alloc(name,DIM,order, reduced_order, mpi_info);
65	if (! Finley_noError()) {
66	Paso_MPIInfo_free( mpi_info );
67	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>0 ? local_NE0*N_PER_E+1 : 0;
136	local_N1=local_NE1>0 ? local_NE1*N_PER_E+1 : 0;
137
138	/* get the number of surface elements */
139
140	NFaceElements=0;
141	if (!periodic[0] && (local_NE0>0)) {
142	NDOF0=N0;
143	if (e_offset0 == 0) NFaceElements+=local_NE1;
144	if (local_NE0+e_offset0 == NE0) NFaceElements+=local_NE1;
145	} else {
146	NDOF0=N0-1;
147	}
148	if (!periodic[1] && (local_NE1>0)) {
149	NDOF1=N1;
150	if (e_offset1 == 0) NFaceElements+=local_NE0;
151	if (local_NE1+e_offset1 == NE1) NFaceElements+=local_NE0;
152	} else {
153	NDOF1=N1-1;
154	}
155
156	/* allocate tables: */
157
158	Finley_NodeFile_allocTable(out->Nodes,local_N0*local_N1);
159	Finley_ElementFile_allocTable(out->Elements,local_NE0*local_NE1);
160	Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
161
162	if (Finley_noError()) {
163	/* create nodes */
164
165	#pragma omp parallel for private(i0,i1,k,global_i0,global_i1)
166	for (i1=0;i1<local_N1;i1++) {
167	for (i0=0;i0<local_N0;i0++) {
168	k=i0+local_N0*i1;
169	global_i0=i0+offset0;
170	global_i1=i1+offset1;
171	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]=Nstride0global_i0+Nstride1global_i1;
174	out->Nodes->Tag[k]=0;
175	out->Nodes->globalDegreesOfFreedom[k]=Nstride0*(global_i0%NDOF0)
176	+Nstride1*(global_i1%NDOF1);
177	}
178	}
179	/* set the elements: */
180	NN=out->Elements->numNodes;
181	#pragma omp parallel for private(i0,i1,k,node0)
182	for (i1=0;i1<local_NE1;i1++) {
183	for (i0=0;i0<local_NE0;i0++) {
184
185	k=i0+local_NE0*i1;
186	node0=Nstride0N_PER_E(i0+e_offset0)+Nstride1N_PER_E(i1+e_offset1);
187
188	out->Elements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1);
189	out->Elements->Tag[k]=0;
190	out->Elements->Owner[k]=myRank;
191
192	out->Elements->Nodes[INDEX2(0,k,NN)]=node0;
193	out->Elements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride0;
194	out->Elements->Nodes[INDEX2(2,k,NN)]=node0+2Nstride1+2Nstride0;
195	out->Elements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride1;
196	out->Elements->Nodes[INDEX2(4,k,NN)]=node0+1*Nstride0;
197	out->Elements->Nodes[INDEX2(5,k,NN)]=node0+Nstride1+2*Nstride0;
198	out->Elements->Nodes[INDEX2(6,k,NN)]=node0+2Nstride1+1Nstride0;
199	out->Elements->Nodes[INDEX2(7,k,NN)]=node0+Nstride1;
200	if (useFullElementOrder) {
201	out->Elements->Nodes[INDEX2(8,k,NN)]=node0+1Nstride1+1Nstride0;
202	}
203	}
204	}
205	/* face elements */
206	NN=out->FaceElements->numNodes;
207	totalNECount=NE0*NE1;
208	faceNECount=0;
209	if (!periodic[0] && (local_NE0>0)) {
210	/* ** elements on boundary 001 (x1=0): */
211
212	if (e_offset0 == 0) {
213	#pragma omp parallel for private(i1,k,node0)
214	for (i1=0;i1<local_NE1;i1++) {
215
216	k=i1+faceNECount;
217	node0=Nstride1N_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+2Nstride1+2Nstride0;
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+2Nstride1+1Nstride0;
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=Nstride0N_PER_E(NE0-1)+Nstride1N_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+2Nstride1+2Nstride0;
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+2Nstride1+1Nstride0;
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+2Nstride1+2Nstride0;
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] && (local_NE0>0)) {
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=Nstride0N_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+2Nstride1+2Nstride0;
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+2Nstride1+1Nstride0;
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=Nstride0N_PER_E(i0+e_offset0)+Nstride1N_PER_E(NE1-1);
306
307	out->FaceElements->Id[k]=i0+e_offset0+totalNECount;
308	out->FaceElements->Tag[k]=20;
309	out->FaceElements->Owner[k]=myRank;
310	if (useElementsOnFace) {
311	out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2Nstride1+2Nstride0;
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+2Nstride1+1Nstride0;
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+2Nstride1+2Nstride0;
321	out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride1;
322	out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2Nstride1+1Nstride0;
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	/* prepare mesh for further calculatuions:*/
336	if (Finley_noError()) {
337	Finley_Mesh_resolveNodeIds(out);
338	}
339	if (Finley_noError()) {
340	Finley_Mesh_prepare(out, optimize);
341	}
342	}
343	if (!Finley_noError()) {
344	Finley_Mesh_free(out);
345	}
346	/* free up memory */
347	Paso_MPIInfo_free( mpi_info );
348	#ifdef Finley_TRACE
349	printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
350	#endif
351
352	return out;
353	}
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision