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;
  Finley_Mesh *mesh_p=NULL;
  char name[LenString_MAX],element_type[LenString_MAX],frm[20], line[LenString_MAX+1];
  char error_msg[LenErrorMsg_MAX];
  double rtmp0, rtmp1;
  double time0=Finley_timer();

  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);
  if (! Finley_noError()) return NULL;

  /* get file handle */
  FILE * 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;
      index_t numElements=0;
      index_t numFaceElements=0;
      fscanf(fileHandle_p, "%d", &totalNumElements);

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