finley/src/Mesh_readGmsh.c

/*
 ************************************************************
 *          Copyright 2006, 2007 by ACcESS MNRF             *
 *                                                          *
 *              http://www.access.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: read mesh */

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

/*   Author: gross@access.edu.au */
/*   Version: $Id: Mesh_read.c 730 2006-05-15 04:03:49Z bcumming $ */

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

#include "Mesh.h"
#include <stdio.h>


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

/*  reads a mesh from a Finley file of name fname */

#define MAX_numNodes_gmsh 10

Finley_Mesh* Finley_Mesh_readGmsh(char* fname ,index_t numDim, index_t order, index_t reduced_order, bool_t optimize_labeling) {

  double version = 1.0;
  int format = 0, size = sizeof(double);
  dim_t numNodes, totalNumElements=0, numTags=0, numNodesPerElement, numNodesPerElement2, element_dim;
  index_t e, i0, j, gmsh_type, partition_id, itmp, final_element_type,  elementary_id;
  index_t numElements=0, numFaceElements=0, *id=NULL, *tag=NULL, *vertices=NULL;
  Finley_Mesh *mesh_p=NULL;
  char line[LenString_MAX+1];
  char error_msg[LenErrorMsg_MAX];
  double rtmp0, rtmp1;
  double time0=Finley_timer();
  FILE * fileHandle_p = NULL;
  ElementTypeId* element_type=NULL;

  Finley_resetError();
 #if PASO_MPI
    sprintf(error_msg,"reading GMSH with MPI is not supported yet.");
    Finley_setError(IO_ERROR,error_msg);
    return NULL;
#else
  /* allocate mesh */

  mesh_p = Finley_Mesh_alloc(fname,numDim,order, reduced_order);
  if (! Finley_noError()) return NULL;

  /* get file handle */
  fileHandle_p = fopen(fname, "r");
  if (fileHandle_p==NULL) {
    sprintf(error_msg,"Opening Gmsh file %s for reading failed.",fname);
    Finley_setError(IO_ERROR,error_msg);
    return NULL;
  }

  /* start reading */
  while(1) {
    if (! Finley_noError()) break;
    /* find line staring with $ */
    do {
      if( ! fgets(line, sizeof(line), fileHandle_p) ) break;
      if(feof(fileHandle_p)) break;
    } while(line[0] != '$');

    if (feof(fileHandle_p)) break;


    /* format */
    if (!strncmp(&line[1], "MeshFormat", 10)) {
      fscanf(fileHandle_p, "%lf %d %d\n", &version, &format, &size);
    }
    /* nodes are read */
    if ( !strncmp(&line[1], "NOD", 3)   ||
         !strncmp(&line[1], "NOE", 3)   ||
         !strncmp(&line[1], "Nodes", 5)    ) {
     
      fscanf(fileHandle_p, "%d", &numNodes);
      if (! Finley_noError()) break;
      Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
      if (! Finley_noError()) break;
      for (i0 = 0; i0 < numNodes; i0++) {
         if (1 == numDim) {
           fscanf(fileHandle_p, "%d %le %le %le\n", &mesh_p->Nodes->Id[i0],
                  &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)], 
                      &rtmp0, 
                      &rtmp1);
         } else if (2 == numDim) {
           fscanf(fileHandle_p, "%d %le %le %le\n", &mesh_p->Nodes->Id[i0],
                  &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],
                  &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)],
                      &rtmp0);
         } else if (3 == numDim) {
           fscanf(fileHandle_p, "%d %le %le %le\n", &mesh_p->Nodes->Id[i0],
                  &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],
                  &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)],
                  &mesh_p->Nodes->Coordinates[INDEX2(2,i0,numDim)]);
   
       
         }
         mesh_p->Nodes->degreeOfFreedom[i0]=mesh_p->Nodes->Id[i0];
         mesh_p->Nodes->Tag[i0]=0;
      }
    }
    /* elements */
    else if(!strncmp(&line[1], "ELM", 3) ||
        !strncmp(&line[1], "Elements", 8)) {

      ElementTypeId final_element_type = NoType;
      ElementTypeId final_face_element_type = NoType;
      numElements=0;
      numFaceElements=0;
      fscanf(fileHandle_p, "%d", &totalNumElements);

      id=TMPMEMALLOC(totalNumElements,index_t);
      tag=TMPMEMALLOC(totalNumElements,index_t);


      element_type=TMPMEMALLOC(totalNumElements,ElementTypeId);
      vertices=TMPMEMALLOC(totalNumElements*MAX_numNodes_gmsh,index_t);
      if (! (Finley_checkPtr(id) || Finley_checkPtr(tag) || Finley_checkPtr(element_type) || Finley_checkPtr(vertices) ) ) {
         /* read all in */
         for(e = 0; e < totalNumElements; e++) {
           fscanf(fileHandle_p, "%d %d", &id[e], &gmsh_type);
           switch (gmsh_type) {
               case 1:  /* line order 1 */
                   element_type[e]=Line2;
                   element_dim=1;
                   numNodesPerElement=2;  
                   break;
               case 2:  /* traingle order 1 */
                   element_type[e]=Tri3;
                   numNodesPerElement= 3;
                   element_dim=2;
                   break;
               case 3:  /* quadrilateral order 1 */
                   element_type[e]=Rec4;
                   numNodesPerElement= 4;
                   element_dim=2;
                   break;
               case 4:  /* tetrahedron order 1 */
                   element_type[e]=Tet4;
                   numNodesPerElement= 4;
                   element_dim=3;
                   break;
               case 5:  /* hexahedron order 1 */
                   element_type[e]=Hex8;
                   numNodesPerElement= 8;
                   element_dim=3;
                   break;
               case 8:  /* line order 2 */
                   element_type[e]=Line3;
                   numNodesPerElement= 3;
                   element_dim=1;
                   break;
               case 9:  /* traingle order 2 */
                   element_type[e]=Tri6;
                   numNodesPerElement= 6;
                   element_dim=2;
                   break;
               case 10:  /* quadrilateral order 2 */
                   element_type[e]=Rec9;
                   numNodesPerElement= 9;
                   element_dim=2;
                   break;
               case 11:  /* tetrahedron order 2 */
                   element_type[e]=Tet10;
                   numNodesPerElement= 10;
                   element_dim=3;
                   break;
               case 15 :  /* point */
                   element_type[e]=Point1;
                   numNodesPerElement= 1;
                   element_dim=0;
                   break;
               default:
                  element_type[e]=NoType;
                  sprintf(error_msg,"Unexected gmsh element type %d in mesh file %s.",gmsh_type,fname);
                  Finley_setError(IO_ERROR,error_msg);
           }
           if (element_dim == numDim) {
              if (final_element_type == NoType) {
                 final_element_type = element_type[e];
              } else if (final_element_type != element_type[e]) {
                  sprintf(error_msg,"Finley can handle a single type of internal elements only.");
                  Finley_setError(IO_ERROR,error_msg);
                  break;
              }
              numElements++;
           } else if (element_dim == numDim-1) {
              if (final_face_element_type == NoType) {
                 final_face_element_type = element_type[e];
              } else if (final_face_element_type != element_type[e]) {
                  sprintf(error_msg,"Finley can handle a single type of face elements only.");
                  Finley_setError(IO_ERROR,error_msg);
                  break;
              }
              numFaceElements++;
           }
           
       if(version <= 1.0){
         fscanf(fileHandle_p, "%d %d %d", &tag[e], &elementary_id, &numNodesPerElement2);
         partition_id = 1;
             if (numNodesPerElement2 != numNodesPerElement) {
                  sprintf(error_msg,"Illegal number of nodes for element %d in mesh file %s.",id[e],fname);
                  Finley_setError(IO_ERROR,error_msg);
             }
       } else {
         fscanf(fileHandle_p, "%d", &numTags);
         elementary_id = tag[e] = partition_id = 1;
             numNodesPerElement2=-1;
         for(j = 0; j < numTags; j++){
           fscanf(fileHandle_p, "%d", &itmp);       
           if (j == 0) {
             tag[e] = itmp;
           } else if (j == 1) {
             elementary_id = itmp;
           } else if (j == 2) {
             partition_id = itmp;
               }
           /* ignore any other tags */
         }
       }
           if (! Finley_noError()) break;
           for(j = 0; j < numNodesPerElement; j++) fscanf(fileHandle_p, "%d", &vertices[INDEX2(j,e,MAX_numNodes_gmsh)]);
         }
         /* all elements have been read, now we have to identify the elements for finley */
     
         if (Finley_noError()) {
           /* first we have to identify the elements to define Elementis and FaceElements */
           if (final_element_type == NoType) {
              if (numDim==1) {
                 final_element_type=Line2;
              } else if (numDim==2) {
                 final_element_type=Tri3;
              } else if (numDim==3) {
                 final_element_type=Tet4;
              }
           }
           if (final_face_element_type == NoType) {
              if (numDim==1) {
                 final_face_element_type=Point1;
              } else if (numDim==2) {
                 final_face_element_type=Line2;
              } else if (numDim==3) {
                 final_face_element_type=Tri3;
              }
           }
printf("ELEMENTS READ\n");
           mesh_p->Elements=Finley_ElementFile_alloc(final_element_type,mesh_p->order, mesh_p->reduced_order);
           mesh_p->FaceElements=Finley_ElementFile_alloc(final_face_element_type,mesh_p->order, mesh_p->reduced_order);
           mesh_p->ContactElements=Finley_ElementFile_alloc(Point1_Contact,mesh_p->order, mesh_p->reduced_order);
           mesh_p->Points=Finley_ElementFile_alloc(Point1,mesh_p->order, mesh_p->reduced_order);
           if (Finley_noError()) {
               Finley_ElementFile_allocTable(mesh_p->Elements, numElements);
               Finley_ElementFile_allocTable(mesh_p->FaceElements, numFaceElements);
               Finley_ElementFile_allocTable(mesh_p->ContactElements, 0);
               Finley_ElementFile_allocTable(mesh_p->Points, 0);
               if (Finley_noError()) {
                   mesh_p->Elements->minColor=0;
                   mesh_p->Elements->maxColor=numElements-1;
                   mesh_p->FaceElements->minColor=0;
                   mesh_p->FaceElements->maxColor=numFaceElements-1;
                   mesh_p->ContactElements->minColor=0;
                   mesh_p->ContactElements->maxColor=0;
                   mesh_p->Points->minColor=0;
                   mesh_p->Points->maxColor=0;
                   numElements=0;
                   numFaceElements=0;
                   for(e = 0; e < totalNumElements; e++) {
                      if (element_type[e] == final_element_type) {
                         mesh_p->Elements->Id[numElements]=id[e];
                         mesh_p->Elements->Tag[numElements]=tag[e];
                         mesh_p->Elements->Color[numElements]=numElements;
                         for (j = 0; j<  mesh_p->Elements->ReferenceElement->Type->numNodes; ++j) 
                                  mesh_p->Elements->Nodes[INDEX2(j, numElements, mesh_p->Elements->ReferenceElement->Type->numNodes)]=
                                                                                                  vertices[INDEX2(j,e,MAX_numNodes_gmsh)];
                         numElements++;
                      } else if (element_type[e] == final_face_element_type) {
                         mesh_p->FaceElements->Id[numFaceElements]=id[e];
                         mesh_p->FaceElements->Tag[numFaceElements]=tag[e];
                         mesh_p->FaceElements->Color[numFaceElements]=numFaceElements;
                         for (j = 0; j<  mesh_p->FaceElements->ReferenceElement->Type->numNodes; ++j) 
                                  mesh_p->FaceElements->Nodes[INDEX2(j, numFaceElements, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]=
                                                                                                  vertices[INDEX2(j,e,MAX_numNodes_gmsh)];
                         numFaceElements++;
                      }
                   }
              }
           }
         }
      }
printf("ELEMENTS reordered\n");
      /* and clean up */
      TMPMEMFREE(id);
      TMPMEMFREE(tag);
      TMPMEMFREE(element_type);
      TMPMEMFREE(vertices);
   }
   /* serach for end of data block */
   do {
      if (!fgets(line, sizeof(line), fileHandle_p)) {
         sprintf(error_msg,"Unexected end of file in %s",fname);
         Finley_setError(IO_ERROR,error_msg);
      }
      if (feof(fileHandle_p)) {
         sprintf(error_msg,"Unexected end of file in %s",fname);
         Finley_setError(IO_ERROR,error_msg);
      }
      if (! Finley_noError()) break;
    } while(line[0] != '$');
  }

  /* close file */
  fclose(fileHandle_p);
  /*   resolve id's : */
  if (Finley_noError()) {
      Finley_Mesh_resolveNodeIds(mesh_p);
  }
  /* rearrange elements: */
  if (Finley_noError()) {
     Finley_Mesh_prepare(mesh_p);
  }

  /* optimize node labeling*/
  if (Finley_noError()) {
      if (optimize_labeling) Finley_Mesh_optimizeNodeLabeling(mesh_p);
  }
  /* that's it */
  #ifdef Finley_TRACE
  printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
  #endif
  #endif
  if (Finley_noError()) {
       if (  ! Finley_Mesh_isPrepared(mesh_p)) {
          Finley_setError(SYSTEM_ERROR,"Mesh is not prepared for calculation. Contact the programmers.");
       }
  } else {
     Finley_Mesh_dealloc(mesh_p);
  }
  return mesh_p;
}
1	gross	935	/*
2			************************************************************
3			* Copyright 2006, 2007 by ACcESS MNRF *
4			* *
5			* http://www.access.edu.au *
6			* Primary Business: Queensland, Australia *
7			* Licensed under the Open Software License version 3.0 *
8			* http://www.opensource.org/licenses/osl-3.0.php *
9			* *
10			************************************************************
11			*/
12
13			/**************************************************************/
14
15			/* Finley: read mesh */
16
17			/**************************************************************/
18
19			/* Author: gross@access.edu.au */
20			/* Version: $Id: Mesh_read.c 730 2006-05-15 04:03:49Z bcumming $ */
21
22			/**************************************************************/
23
24			#include "Mesh.h"
25			#include <stdio.h>
26
27
28			/**************************************************************/
29
30			/* reads a mesh from a Finley file of name fname */
31
32			#define MAX_numNodes_gmsh 10
33
34			Finley_Mesh* Finley_Mesh_readGmsh(char* fname ,index_t numDim, index_t order, index_t reduced_order, bool_t optimize_labeling) {
35
36			double version = 1.0;
37			int format = 0, size = sizeof(double);
38			dim_t numNodes, totalNumElements=0, numTags=0, numNodesPerElement, numNodesPerElement2, element_dim;
39			index_t e, i0, j, gmsh_type, partition_id, itmp, final_element_type, elementary_id;
40	gross	1028	index_t numElements=0, numFaceElements=0, id=NULL, tag=NULL, *vertices=NULL;
41	gross	935	Finley_Mesh *mesh_p=NULL;
42	gross	1028	char line[LenString_MAX+1];
43	gross	935	char error_msg[LenErrorMsg_MAX];
44			double rtmp0, rtmp1;
45			double time0=Finley_timer();
46	gross	1028	FILE * fileHandle_p = NULL;
47			ElementTypeId* element_type=NULL;
48	gross	935
49			Finley_resetError();
50			#if PASO_MPI
51			sprintf(error_msg,"reading GMSH with MPI is not supported yet.");
52			Finley_setError(IO_ERROR,error_msg);
53			return NULL;
54			#else
55			/* allocate mesh */
56	gross	997
57	gross	1062	mesh_p = Finley_Mesh_alloc(fname,numDim,order, reduced_order);
58	gross	935	if (! Finley_noError()) return NULL;
59
60			/* get file handle */
61	gross	1028	fileHandle_p = fopen(fname, "r");
62	gross	935	if (fileHandle_p==NULL) {
63			sprintf(error_msg,"Opening Gmsh file %s for reading failed.",fname);
64			Finley_setError(IO_ERROR,error_msg);
65			return NULL;
66			}
67
68			/* start reading */
69			while(1) {
70			if (! Finley_noError()) break;
71			/* find line staring with $ */
72			do {
73			if( ! fgets(line, sizeof(line), fileHandle_p) ) break;
74			if(feof(fileHandle_p)) break;
75			} while(line[0] != '$');
76
77			if (feof(fileHandle_p)) break;
78
79
80			/* format */
81			if (!strncmp(&line[1], "MeshFormat", 10)) {
82			fscanf(fileHandle_p, "%lf %d %d\n", &version, &format, &size);
83			}
84			/* nodes are read */
85			if ( !strncmp(&line[1], "NOD", 3) \|\|
86			!strncmp(&line[1], "NOE", 3) \|\|
87			!strncmp(&line[1], "Nodes", 5) ) {
88
89			fscanf(fileHandle_p, "%d", &numNodes);
90			if (! Finley_noError()) break;
91			Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
92			if (! Finley_noError()) break;
93			for (i0 = 0; i0 < numNodes; i0++) {
94			if (1 == numDim) {
95			fscanf(fileHandle_p, "%d %le %le %le\n", &mesh_p->Nodes->Id[i0],
96			&mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],
97			&rtmp0,
98			&rtmp1);
99			} else if (2 == numDim) {
100			fscanf(fileHandle_p, "%d %le %le %le\n", &mesh_p->Nodes->Id[i0],
101			&mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],
102			&mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)],
103			&rtmp0);
104			} else if (3 == numDim) {
105			fscanf(fileHandle_p, "%d %le %le %le\n", &mesh_p->Nodes->Id[i0],
106			&mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],
107			&mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)],
108			&mesh_p->Nodes->Coordinates[INDEX2(2,i0,numDim)]);
109
110
111			}
112			mesh_p->Nodes->degreeOfFreedom[i0]=mesh_p->Nodes->Id[i0];
113			mesh_p->Nodes->Tag[i0]=0;
114			}
115			}
116			/* elements */
117			else if(!strncmp(&line[1], "ELM", 3) \|\|
118			!strncmp(&line[1], "Elements", 8)) {
119
120			ElementTypeId final_element_type = NoType;
121			ElementTypeId final_face_element_type = NoType;
122	gross	1028	numElements=0;
123			numFaceElements=0;
124	gross	935	fscanf(fileHandle_p, "%d", &totalNumElements);
125
126	gross	1028	id=TMPMEMALLOC(totalNumElements,index_t);
127			tag=TMPMEMALLOC(totalNumElements,index_t);
128
129
130			element_type=TMPMEMALLOC(totalNumElements,ElementTypeId);
131			vertices=TMPMEMALLOC(totalNumElements*MAX_numNodes_gmsh,index_t);
132	gross	935	if (! (Finley_checkPtr(id) \|\| Finley_checkPtr(tag) \|\| Finley_checkPtr(element_type) \|\| Finley_checkPtr(vertices) ) ) {
133			/* read all in */
134			for(e = 0; e < totalNumElements; e++) {
135			fscanf(fileHandle_p, "%d %d", &id[e], &gmsh_type);
136			switch (gmsh_type) {
137			case 1: /* line order 1 */
138			element_type[e]=Line2;
139			element_dim=1;
140			numNodesPerElement=2;
141			break;
142			case 2: /* traingle order 1 */
143			element_type[e]=Tri3;
144			numNodesPerElement= 3;
145			element_dim=2;
146			break;
147			case 3: /* quadrilateral order 1 */
148			element_type[e]=Rec4;
149			numNodesPerElement= 4;
150			element_dim=2;
151			break;
152			case 4: /* tetrahedron order 1 */
153			element_type[e]=Tet4;
154			numNodesPerElement= 4;
155			element_dim=3;
156			break;
157			case 5: /* hexahedron order 1 */
158			element_type[e]=Hex8;
159			numNodesPerElement= 8;
160			element_dim=3;
161			break;
162			case 8: /* line order 2 */
163			element_type[e]=Line3;
164			numNodesPerElement= 3;
165			element_dim=1;
166			break;
167			case 9: /* traingle order 2 */
168			element_type[e]=Tri6;
169			numNodesPerElement= 6;
170			element_dim=2;
171			break;
172			case 10: /* quadrilateral order 2 */
173			element_type[e]=Rec9;
174			numNodesPerElement= 9;
175			element_dim=2;
176			break;
177			case 11: /* tetrahedron order 2 */
178			element_type[e]=Tet10;
179			numNodesPerElement= 10;
180			element_dim=3;
181			break;
182			case 15 : /* point */
183			element_type[e]=Point1;
184			numNodesPerElement= 1;
185			element_dim=0;
186			break;
187			default:
188			element_type[e]=NoType;
189			sprintf(error_msg,"Unexected gmsh element type %d in mesh file %s.",gmsh_type,fname);
190			Finley_setError(IO_ERROR,error_msg);
191			}
192			if (element_dim == numDim) {
193			if (final_element_type == NoType) {
194			final_element_type = element_type[e];
195			} else if (final_element_type != element_type[e]) {
196			sprintf(error_msg,"Finley can handle a single type of internal elements only.");
197			Finley_setError(IO_ERROR,error_msg);
198			break;
199			}
200			numElements++;
201			} else if (element_dim == numDim-1) {
202			if (final_face_element_type == NoType) {
203			final_face_element_type = element_type[e];
204			} else if (final_face_element_type != element_type[e]) {
205			sprintf(error_msg,"Finley can handle a single type of face elements only.");
206			Finley_setError(IO_ERROR,error_msg);
207			break;
208			}
209			numFaceElements++;
210			}
211
212			if(version <= 1.0){
213			fscanf(fileHandle_p, "%d %d %d", &tag[e], &elementary_id, &numNodesPerElement2);
214			partition_id = 1;
215			if (numNodesPerElement2 != numNodesPerElement) {
216			sprintf(error_msg,"Illegal number of nodes for element %d in mesh file %s.",id[e],fname);
217			Finley_setError(IO_ERROR,error_msg);
218			}
219			} else {
220			fscanf(fileHandle_p, "%d", &numTags);
221			elementary_id = tag[e] = partition_id = 1;
222			numNodesPerElement2=-1;
223			for(j = 0; j < numTags; j++){
224			fscanf(fileHandle_p, "%d", &itmp);
225			if (j == 0) {
226			tag[e] = itmp;
227			} else if (j == 1) {
228			elementary_id = itmp;
229			} else if (j == 2) {
230			partition_id = itmp;
231			}
232	gross	1028	/* ignore any other tags */
233	gross	935	}
234			}
235			if (! Finley_noError()) break;
236			for(j = 0; j < numNodesPerElement; j++) fscanf(fileHandle_p, "%d", &vertices[INDEX2(j,e,MAX_numNodes_gmsh)]);
237			}
238			/* all elements have been read, now we have to identify the elements for finley */
239	gross	997
240	gross	935	if (Finley_noError()) {
241	gross	1003	/* first we have to identify the elements to define Elementis and FaceElements */
242			if (final_element_type == NoType) {
243			if (numDim==1) {
244			final_element_type=Line2;
245			} else if (numDim==2) {
246			final_element_type=Tri3;
247	gross	1092	} else if (numDim==3) {
248	gross	1003	final_element_type=Tet4;
249			}
250			}
251			if (final_face_element_type == NoType) {
252			if (numDim==1) {
253			final_face_element_type=Point1;
254			} else if (numDim==2) {
255			final_face_element_type=Line2;
256	gross	1092	} else if (numDim==3) {
257	gross	1003	final_face_element_type=Tri3;
258			}
259			}
260	gross	1156	printf("ELEMENTS READ\n");
261	gross	1062	mesh_p->Elements=Finley_ElementFile_alloc(final_element_type,mesh_p->order, mesh_p->reduced_order);
262			mesh_p->FaceElements=Finley_ElementFile_alloc(final_face_element_type,mesh_p->order, mesh_p->reduced_order);
263			mesh_p->ContactElements=Finley_ElementFile_alloc(Point1_Contact,mesh_p->order, mesh_p->reduced_order);
264			mesh_p->Points=Finley_ElementFile_alloc(Point1,mesh_p->order, mesh_p->reduced_order);
265	gross	935	if (Finley_noError()) {
266			Finley_ElementFile_allocTable(mesh_p->Elements, numElements);
267			Finley_ElementFile_allocTable(mesh_p->FaceElements, numFaceElements);
268			Finley_ElementFile_allocTable(mesh_p->ContactElements, 0);
269			Finley_ElementFile_allocTable(mesh_p->Points, 0);
270			if (Finley_noError()) {
271			mesh_p->Elements->minColor=0;
272			mesh_p->Elements->maxColor=numElements-1;
273			mesh_p->FaceElements->minColor=0;
274			mesh_p->FaceElements->maxColor=numFaceElements-1;
275			mesh_p->ContactElements->minColor=0;
276			mesh_p->ContactElements->maxColor=0;
277			mesh_p->Points->minColor=0;
278			mesh_p->Points->maxColor=0;
279			numElements=0;
280			numFaceElements=0;
281			for(e = 0; e < totalNumElements; e++) {
282			if (element_type[e] == final_element_type) {
283			mesh_p->Elements->Id[numElements]=id[e];
284			mesh_p->Elements->Tag[numElements]=tag[e];
285			mesh_p->Elements->Color[numElements]=numElements;
286			for (j = 0; j< mesh_p->Elements->ReferenceElement->Type->numNodes; ++j)
287			mesh_p->Elements->Nodes[INDEX2(j, numElements, mesh_p->Elements->ReferenceElement->Type->numNodes)]=
288			vertices[INDEX2(j,e,MAX_numNodes_gmsh)];
289			numElements++;
290			} else if (element_type[e] == final_face_element_type) {
291			mesh_p->FaceElements->Id[numFaceElements]=id[e];
292			mesh_p->FaceElements->Tag[numFaceElements]=tag[e];
293			mesh_p->FaceElements->Color[numFaceElements]=numFaceElements;
294			for (j = 0; j< mesh_p->FaceElements->ReferenceElement->Type->numNodes; ++j)
295			mesh_p->FaceElements->Nodes[INDEX2(j, numFaceElements, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]=
296			vertices[INDEX2(j,e,MAX_numNodes_gmsh)];
297			numFaceElements++;
298			}
299			}
300			}
301			}
302			}
303			}
304	gross	1156	printf("ELEMENTS reordered\n");
305	gross	935	/* and clean up */
306			TMPMEMFREE(id);
307			TMPMEMFREE(tag);
308			TMPMEMFREE(element_type);
309			TMPMEMFREE(vertices);
310			}
311			/* serach for end of data block */
312			do {
313			if (!fgets(line, sizeof(line), fileHandle_p)) {
314			sprintf(error_msg,"Unexected end of file in %s",fname);
315			Finley_setError(IO_ERROR,error_msg);
316			}
317			if (feof(fileHandle_p)) {
318	gross	940	sprintf(error_msg,"Unexected end of file in %s",fname);
319	gross	935	Finley_setError(IO_ERROR,error_msg);
320			}
321			if (! Finley_noError()) break;
322			} while(line[0] != '$');
323			}
324
325			/* close file */
326			fclose(fileHandle_p);
327	gross	1156	/* resolve id's : */
328			if (Finley_noError()) {
329			Finley_Mesh_resolveNodeIds(mesh_p);
330	gross	935	}
331			/* rearrange elements: */
332	gross	1156	if (Finley_noError()) {
333			Finley_Mesh_prepare(mesh_p);
334			}
335
336			/* optimize node labeling*/
337			if (Finley_noError()) {
338			if (optimize_labeling) Finley_Mesh_optimizeNodeLabeling(mesh_p);
339			}
340	gross	935	/* that's it */
341			#ifdef Finley_TRACE
342			printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
343			#endif
344			#endif
345	gross	964	if (Finley_noError()) {
346			if ( ! Finley_Mesh_isPrepared(mesh_p)) {
347			Finley_setError(SYSTEM_ERROR,"Mesh is not prepared for calculation. Contact the programmers.");
348			}
349	gross	1156	} else {
350			Finley_Mesh_dealloc(mesh_p);
351	gross	964	}
352	gross	935	return mesh_p;
353			}