finley/src/Mesh_rec4.c


/*******************************************************
*
* Copyright (c) 2003-2009 by University of Queensland
* Earth Systems Science Computational Center (ESSCC)
* http://www.uq.edu.au/esscc
*
* 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 first order elements (Rec4) 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_Rec4(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 2
  dim_t N0,N1,NE0,NE1,i0,i1,k,Nstride0=0,Nstride1=0, local_NE0, local_NE1, local_N0=0, local_N1=0, global_i0, global_i1;
  index_t offset0=0, offset1=0, e_offset0=0, e_offset1=0;
  dim_t totalNECount,faceNECount,NDOF0=0,NDOF1=0,NFaceElements,NN;
  Finley_ReferenceElementSet *refPoints=NULL, *refContactElements=NULL, *refFaceElements=NULL, *refElements=NULL;
  index_t node0, myRank;
  Finley_Mesh* out;
  Paso_MPIInfo *mpi_info = NULL;
  char name[50];
  #ifdef Finley_TRACE
  double time0=Finley_timer();
  #endif

  /* 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;
  }
  refElements= Finley_ReferenceElementSet_alloc(Rec4,out->order,out->reduced_order);
  if (useElementsOnFace) {
        refFaceElements=Finley_ReferenceElementSet_alloc(Rec4Face, out->order, out->reduced_order);
        refContactElements=Finley_ReferenceElementSet_alloc(Rec4Face_Contact, out->order, out->reduced_order);
  } else {
        refFaceElements=Finley_ReferenceElementSet_alloc(Line2, out->order, out->reduced_order);
        refContactElements=Finley_ReferenceElementSet_alloc(Line2_Contact, out->order, out->reduced_order);
  }
  refPoints=Finley_ReferenceElementSet_alloc(Point1, out->order, out->reduced_order);
  
  if ( Finley_noError()) {
  
      Finley_Mesh_setPoints(out,Finley_ElementFile_alloc(refPoints, mpi_info));
      Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(refContactElements, mpi_info));
      Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(refFaceElements, mpi_info));
      Finley_Mesh_setElements(out,Finley_ElementFile_alloc(refElements, mpi_info));
 
    /* 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+Nstride0;
           out->Elements->Nodes[INDEX2(2,k,NN)]=node0+Nstride1+Nstride0;
           out->Elements->Nodes[INDEX2(3,k,NN)]=node0+Nstride1;
         }
     }
     /* 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+Nstride1;
                  out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0;
                  out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride0;
                  out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride1+Nstride0;
               } else {
                  out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride1;
                  out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0;
               }
           }
           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+Nstride0;
                  out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride1+Nstride0;
                  out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride1;
                  out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0;
               } else {
                  out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride0;
                  out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride1+Nstride0;
               }
           }
           faceNECount+=local_NE1;
         }
         totalNECount+=NE1;
     }
     if (!periodic[1]  && (local_NE1>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+Nstride0;
                   out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride1+Nstride0;
                   out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride1;
               } else {
                   out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0;
                   out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+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+Nstride1+Nstride0;
                    out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride1;
                    out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0;
                    out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride0;
               } else {
                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride1+Nstride0;
                    out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride1;
               }
           }
           faceNECount+=local_NE0;
        }
        totalNECount+=NE0;
     }
  } 
  if (Finley_noError()) {    
     /* 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);
  }

  /* free up memory */
  Finley_ReferenceElementSet_dealloc(refPoints);
  Finley_ReferenceElementSet_dealloc(refContactElements);
  Finley_ReferenceElementSet_dealloc(refFaceElements);
  Finley_ReferenceElementSet_dealloc(refElements);
  Paso_MPIInfo_free( mpi_info );  

  return out;
}
1
2	/*******************************************************
3	*
4	* Copyright (c) 2003-2009 by University of Queensland
5	* Earth Systems Science Computational Center (ESSCC)
6	* http://www.uq.edu.au/esscc
7	*
8	* Primary Business: Queensland, Australia
9	* Licensed under the Open Software License version 3.0
10	* http://www.opensource.org/licenses/osl-3.0.php
11	*
12	*******************************************************/
13
14
15	/**************************************************************/
16
17	/* Finley: generates rectangular meshes */
18
19	/* Generates a numElements[0] x numElements[1] mesh with first order elements (Rec4) in the rectangle */
20	/* [0,Length[0]] x [0,Length[1]]. order is the desired accuracy of the integration scheme. */
21
22
23	/**************************************************************/
24
25	#include "RectangularMesh.h"
26
27
28	Finley_Mesh* Finley_RectangularMesh_Rec4(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 1
38	#define DIM 2
39	dim_t N0,N1,NE0,NE1,i0,i1,k,Nstride0=0,Nstride1=0, local_NE0, local_NE1, local_N0=0, local_N1=0, global_i0, global_i1;
40	index_t offset0=0, offset1=0, e_offset0=0, e_offset1=0;
41	dim_t totalNECount,faceNECount,NDOF0=0,NDOF1=0,NFaceElements,NN;
42	Finley_ReferenceElementSet refPoints=NULL, refContactElements=NULL, refFaceElements=NULL, refElements=NULL;
43	index_t node0, myRank;
44	Finley_Mesh* out;
45	Paso_MPIInfo *mpi_info = NULL;
46	char name[50];
47	#ifdef Finley_TRACE
48	double time0=Finley_timer();
49	#endif
50
51	/* get MPI information */
52	mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
53	if (! Finley_noError()) {
54	return NULL;
55	}
56	myRank=mpi_info->rank;
57
58	/* set up the global dimensions of the mesh */
59
60	NE0=MAX(1,numElements[0]);
61	NE1=MAX(1,numElements[1]);
62	N0=N_PER_E*NE0+1;
63	N1=N_PER_E*NE1+1;
64
65	/* allocate mesh: */
66	sprintf(name,"Rectangular %d x %d mesh",N0,N1);
67	out=Finley_Mesh_alloc(name,DIM,order, reduced_order, mpi_info);
68	if (! Finley_noError()) {
69	Paso_MPIInfo_free( mpi_info );
70	return NULL;
71	}
72	refElements= Finley_ReferenceElementSet_alloc(Rec4,out->order,out->reduced_order);
73	if (useElementsOnFace) {
74	refFaceElements=Finley_ReferenceElementSet_alloc(Rec4Face, out->order, out->reduced_order);
75	refContactElements=Finley_ReferenceElementSet_alloc(Rec4Face_Contact, out->order, out->reduced_order);
76	} else {
77	refFaceElements=Finley_ReferenceElementSet_alloc(Line2, out->order, out->reduced_order);
78	refContactElements=Finley_ReferenceElementSet_alloc(Line2_Contact, out->order, out->reduced_order);
79	}
80	refPoints=Finley_ReferenceElementSet_alloc(Point1, out->order, out->reduced_order);
81
82	if ( Finley_noError()) {
83
84	Finley_Mesh_setPoints(out,Finley_ElementFile_alloc(refPoints, mpi_info));
85	Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(refContactElements, mpi_info));
86	Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(refFaceElements, mpi_info));
87	Finley_Mesh_setElements(out,Finley_ElementFile_alloc(refElements, mpi_info));
88
89	/* work out the largest dimension */
90	if (N1==MAX(N0,N1)) {
91	Nstride0=1;
92	Nstride1=N0;
93	local_NE0=NE0;
94	e_offset0=0;
95	Paso_MPIInfo_Split(mpi_info,NE1,&local_NE1,&e_offset1);
96	} else {
97	Nstride0=N1;
98	Nstride1=1;
99	Paso_MPIInfo_Split(mpi_info,NE0,&local_NE0,&e_offset0);
100	local_NE1=NE1;
101	e_offset1=0;
102	}
103	offset0=e_offset0*N_PER_E;
104	offset1=e_offset1*N_PER_E;
105	local_N0=local_NE0>0 ? local_NE0*N_PER_E+1 : 0;
106	local_N1=local_NE1>0 ? local_NE1*N_PER_E+1 : 0;
107
108	/* get the number of surface elements */
109
110	NFaceElements=0;
111	if (!periodic[0] && (local_NE0>0)) {
112	NDOF0=N0;
113	if (e_offset0 == 0) NFaceElements+=local_NE1;
114	if (local_NE0+e_offset0 == NE0) NFaceElements+=local_NE1;
115	} else {
116	NDOF0=N0-1;
117	}
118	if (!periodic[1] && (local_NE1>0)) {
119	NDOF1=N1;
120	if (e_offset1 == 0) NFaceElements+=local_NE0;
121	if (local_NE1+e_offset1 == NE1) NFaceElements+=local_NE0;
122	} else {
123	NDOF1=N1-1;
124	}
125
126	/* allocate tables: */
127
128
129	Finley_NodeFile_allocTable(out->Nodes,local_N0*local_N1);
130	Finley_ElementFile_allocTable(out->Elements,local_NE0*local_NE1);
131	Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
132
133	}
134	if (Finley_noError()) {
135	/* create nodes */
136	#pragma omp parallel for private(i0,i1,k,global_i0,global_i1)
137	for (i1=0;i1<local_N1;i1++) {
138	for (i0=0;i0<local_N0;i0++) {
139	k=i0+local_N0*i1;
140	global_i0=i0+offset0;
141	global_i1=i1+offset1;
142	out->Nodes->Coordinates[INDEX2(0,k,DIM)]=DBLE(global_i0)/DBLE(N0-1)*Length[0];
143	out->Nodes->Coordinates[INDEX2(1,k,DIM)]=DBLE(global_i1)/DBLE(N1-1)*Length[1];
144	out->Nodes->Id[k]=Nstride0global_i0+Nstride1global_i1;
145	out->Nodes->Tag[k]=0;
146	out->Nodes->globalDegreesOfFreedom[k]=Nstride0*(global_i0%NDOF0)
147	+Nstride1*(global_i1%NDOF1);
148	}
149	}
150	/* set the elements: */
151	NN=out->Elements->numNodes;
152	#pragma omp parallel for private(i0,i1,k,node0)
153	for (i1=0;i1<local_NE1;i1++) {
154	for (i0=0;i0<local_NE0;i0++) {
155
156	k=i0+local_NE0*i1;
157	node0=Nstride0N_PER_E(i0+e_offset0)+Nstride1N_PER_E(i1+e_offset1);
158
159	out->Elements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1);
160	out->Elements->Tag[k]=0;
161	out->Elements->Owner[k]=myRank;
162
163	out->Elements->Nodes[INDEX2(0,k,NN)]=node0;
164	out->Elements->Nodes[INDEX2(1,k,NN)]=node0+Nstride0;
165	out->Elements->Nodes[INDEX2(2,k,NN)]=node0+Nstride1+Nstride0;
166	out->Elements->Nodes[INDEX2(3,k,NN)]=node0+Nstride1;
167	}
168	}
169	/* face elements */
170	NN=out->FaceElements->numNodes;
171	totalNECount=NE0*NE1;
172	faceNECount=0;
173	if (!periodic[0] && (local_NE0>0)) {
174	/* ** elements on boundary 001 (x1=0): */
175
176	if (e_offset0 == 0) {
177	#pragma omp parallel for private(i1,k,node0)
178	for (i1=0;i1<local_NE1;i1++) {
179
180	k=i1+faceNECount;
181	node0=Nstride1N_PER_E(i1+e_offset1);
182
183	out->FaceElements->Id[k]=i1+e_offset1+totalNECount;
184	out->FaceElements->Tag[k]=1;
185	out->FaceElements->Owner[k]=myRank;
186	if (useElementsOnFace) {
187	out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride1;
188	out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0;
189	out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride0;
190	out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride1+Nstride0;
191	} else {
192	out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride1;
193	out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0;
194	}
195	}
196	faceNECount+=local_NE1;
197	}
198	totalNECount+=NE1;
199	/* ** elements on boundary 002 (x1=1): */
200	if (local_NE0+e_offset0 == NE0) {
201	#pragma omp parallel for private(i1,k,node0)
202	for (i1=0;i1<local_NE1;i1++) {
203	k=i1+faceNECount;
204	node0=Nstride0N_PER_E(NE0-1)+Nstride1N_PER_E(i1+e_offset1);
205
206	out->FaceElements->Id[k]=(i1+e_offset1)+totalNECount;
207	out->FaceElements->Tag[k]=2;
208	out->FaceElements->Owner[k]=myRank;
209
210	if (useElementsOnFace) {
211	out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride0;
212	out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride1+Nstride0;
213	out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride1;
214	out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0;
215	} else {
216	out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride0;
217	out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride1+Nstride0;
218	}
219	}
220	faceNECount+=local_NE1;
221	}
222	totalNECount+=NE1;
223	}
224	if (!periodic[1] && (local_NE1>0)) {
225	/* ** elements on boundary 010 (x2=0): */
226	if (e_offset1 == 0) {
227	#pragma omp parallel for private(i0,k,node0)
228	for (i0=0;i0<local_NE0;i0++) {
229	k=i0+faceNECount;
230	node0=Nstride0N_PER_E(i0+e_offset0);
231
232	out->FaceElements->Id[k]=e_offset0+i0+totalNECount;
233	out->FaceElements->Tag[k]=10;
234	out->FaceElements->Owner[k]=myRank;
235
236	if (useElementsOnFace) {
237	out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0;
238	out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride0;
239	out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride1+Nstride0;
240	out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride1;
241	} else {
242	out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0;
243	out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride0;
244	}
245	}
246	faceNECount+=local_NE0;
247	}
248	totalNECount+=NE0;
249	/* ** elements on boundary 020 (x2=1): */
250	if (local_NE1+e_offset1 == NE1) {
251	#pragma omp parallel for private(i0,k,node0)
252	for (i0=0;i0<local_NE0;i0++) {
253	k=i0+faceNECount;
254	node0=Nstride0N_PER_E(i0+e_offset0)+Nstride1N_PER_E(NE1-1);
255
256	out->FaceElements->Id[k]=i0+e_offset0+totalNECount;
257	out->FaceElements->Tag[k]=20;
258	out->FaceElements->Owner[k]=myRank;
259	if (useElementsOnFace) {
260	out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride1+Nstride0;
261	out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride1;
262	out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0;
263	out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride0;
264	} else {
265	out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride1+Nstride0;
266	out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride1;
267	}
268	}
269	faceNECount+=local_NE0;
270	}
271	totalNECount+=NE0;
272	}
273	}
274	if (Finley_noError()) {
275	/* add tag names */
276	Finley_Mesh_addTagMap(out,"top", 20);
277	Finley_Mesh_addTagMap(out,"bottom", 10);
278	Finley_Mesh_addTagMap(out,"left", 1);
279	Finley_Mesh_addTagMap(out,"right", 2);
280	}
281	/* prepare mesh for further calculatuions:*/
282	if (Finley_noError()) {
283	Finley_Mesh_resolveNodeIds(out);
284	}
285	if (Finley_noError()) {
286	Finley_Mesh_prepare(out, optimize);
287	}
288
289	/* free up memory */
290	Finley_ReferenceElementSet_dealloc(refPoints);
291	Finley_ReferenceElementSet_dealloc(refContactElements);
292	Finley_ReferenceElementSet_dealloc(refFaceElements);
293	Finley_ReferenceElementSet_dealloc(refElements);
294	Paso_MPIInfo_free( mpi_info );
295
296	return out;
297	}
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision