/[escript]/branches/doubleplusgood/dudley/src/Mesh_readGmsh.cpp
ViewVC logotype

Diff of /branches/doubleplusgood/dudley/src/Mesh_readGmsh.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

temp_trunk_copy/finley/src/Mesh_readGmsh.c revision 1384 by phornby, Fri Jan 11 02:29:38 2008 UTC trunk/dudley/src/Mesh_readGmsh.c revision 3259 by jfenwick, Mon Oct 11 01:48:14 2010 UTC
# Line 1  Line 1 
1    
 /* $Id$ */  
   
2  /*******************************************************  /*******************************************************
3   *  *
4   *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2010 by University of Queensland
5   *       Copyright 2007 by University of Queensland  * Earth Systems Science Computational Center (ESSCC)
6   *  * http://www.uq.edu.au/esscc
7   *                http://esscc.uq.edu.au  *
8   *        Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
9   *  Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
10   *     http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
11   *  *
12   *******************************************************/  *******************************************************/
13    
14  /**************************************************************/  /**************************************************************/
15    
16  /*   Finley: read mesh */  /*   Dudley: read mesh */
17    
18  /**************************************************************/  /**************************************************************/
19    
20  #include "Mesh.h"  #include "Mesh.h"
21  #include <stdio.h>  #include <stdio.h>
22    
23    #define FSCANF_CHECK(scan_ret, reason) { if (scan_ret == EOF) { perror(reason); Dudley_setError(IO_ERROR,"scan error while reading dudley file"); return NULL;} }
24    
25  /**************************************************************/  /**************************************************************/
26    
27  /*  reads a mesh from a Finley file of name fname */  /*  reads a mesh from a Dudley file of name fname */
   
 #define MAX_numNodes_gmsh 10  
28    
29  Finley_Mesh* Finley_Mesh_readGmsh(char* fname ,index_t numDim, index_t order, index_t reduced_order, bool_t optimize) {  #define MAX_numNodes_gmsh 20
30    
31    double version = 1.0;  Dudley_Mesh *Dudley_Mesh_readGmsh(char *fname, index_t numDim, index_t order, index_t reduced_order, bool_t optimize,
32    int format = 0, size = sizeof(double);                    bool_t useMacroElements)
33    dim_t numNodes, totalNumElements=0, numTags=0, numNodesPerElement, numNodesPerElement2, element_dim;  {
34    index_t e, i0, j, gmsh_type, partition_id, itmp, final_element_type,  elementary_id;  
35    index_t numElements=0, numFaceElements=0, *id=NULL, *tag=NULL, *vertices=NULL;      double version = 1.0;
36    Finley_Mesh *mesh_p=NULL;      int format = 0, size = sizeof(double), scan_ret;
37    char line[LenString_MAX+1];      dim_t numNodes, totalNumElements = 0, numTags = 0, numNodesPerElement = 0, numNodesPerElement2, element_dim = 0;
38    char error_msg[LenErrorMsg_MAX];      index_t e, i0, j, gmsh_type, partition_id, itmp, elementary_id;
39    double rtmp0, rtmp1;      index_t numElements = 0, numFaceElements = 0, *id = NULL, *tag = NULL, *vertices = NULL;
40    double time0=Finley_timer();      Dudley_Mesh *mesh_p = NULL;
41    FILE * fileHandle_p = NULL;      char line[LenString_MAX + 1];
42    ElementTypeId* element_type=NULL;      char error_msg[LenErrorMsg_MAX];
43        double rtmp0, rtmp1;
44    Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );  #ifdef Dudley_TRACE
45    Finley_resetError();      double time0 = Dudley_timer();
46    if (mpi_info->size>1) {  #endif
47      Finley_setError(IO_ERROR,"reading GMSH with MPI is not supported yet.");      FILE *fileHandle_p = NULL;
48      Paso_MPIInfo_free( mpi_info );      Dudley_ElementTypeId *element_type = NULL;
49      return NULL;  
50    } else {      Esys_MPIInfo *mpi_info = Esys_MPIInfo_alloc(MPI_COMM_WORLD);
51        Dudley_resetError();
52       /* allocate mesh */      if (mpi_info->size > 1)
53          {
54       mesh_p = Finley_Mesh_alloc(fname,numDim,order, reduced_order, mpi_info);      Dudley_setError(IO_ERROR, "reading GMSH with MPI is not supported yet.");
55       if (! Finley_noError()) return NULL;      Esys_MPIInfo_free(mpi_info);
56          return NULL;
57       /* get file handle */      }
58       fileHandle_p = fopen(fname, "r");      else
59       if (fileHandle_p==NULL) {      {
60         sprintf(error_msg,"Opening Gmsh file %s for reading failed.",fname);  
61         Finley_setError(IO_ERROR,error_msg);      /* allocate mesh */
62         Paso_MPIInfo_free( mpi_info );  
63         return NULL;      mesh_p = Dudley_Mesh_alloc(fname, numDim, mpi_info);
64       }      if (!Dudley_noError())
65              return NULL;
66       /* start reading */  
67       while(1) {      /* get file handle */
68         if (! Finley_noError()) break;      fileHandle_p = fopen(fname, "r");
69         /* find line staring with $ */      if (fileHandle_p == NULL)
70         do {      {
71           if( ! fgets(line, sizeof(line), fileHandle_p) ) break;          sprintf(error_msg, "Opening Gmsh file %s for reading failed.", fname);
72           if(feof(fileHandle_p)) break;          Dudley_setError(IO_ERROR, error_msg);
73         } while(line[0] != '$');          Esys_MPIInfo_free(mpi_info);
74              return NULL;
75         if (feof(fileHandle_p)) break;      }
76      
77          /* start reading */
78         /* format */      while (1)
79         if (!strncmp(&line[1], "MeshFormat", 10)) {      {
80           fscanf(fileHandle_p, "%lf %d %d\n", &version, &format, &size);          if (!Dudley_noError())
81         }          break;
82         /* nodes are read */          /* find line staring with $ */
83         if ( !strncmp(&line[1], "NOD", 3)   ||          do
84              !strncmp(&line[1], "NOE", 3)   ||          {
85              !strncmp(&line[1], "Nodes", 5)    ) {          if (!fgets(line, sizeof(line), fileHandle_p))
86                        break;
87           fscanf(fileHandle_p, "%d", &numNodes);          if (feof(fileHandle_p))
88           if (! Finley_noError()) break;              break;
89           Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);          }
90           if (! Finley_noError()) break;          while (line[0] != '$');
91           for (i0 = 0; i0 < numNodes; i0++) {  
92              if (1 == numDim) {          if (feof(fileHandle_p))
93             fscanf(fileHandle_p, "%d %le %le %le\n", &mesh_p->Nodes->Id[i0],          break;
94                    &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],  
95                           &rtmp0,          /* format */
96                           &rtmp1);          if (!strncmp(&line[1], "MeshFormat", 10))
97              } else if (2 == numDim) {          {
98             fscanf(fileHandle_p, "%d %le %le %le\n", &mesh_p->Nodes->Id[i0],          scan_ret = fscanf(fileHandle_p, "%lf %d %d\n", &version, &format, &size);
99                    &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],          FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
100                    &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)],          }
101                           &rtmp0);          /* nodes are read */
102              } else if (3 == numDim) {          if (!strncmp(&line[1], "NOD", 3) || !strncmp(&line[1], "NOE", 3) || !strncmp(&line[1], "Nodes", 5))
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)],          scan_ret = fscanf(fileHandle_p, "%d", &numNodes);
106                    &mesh_p->Nodes->Coordinates[INDEX2(2,i0,numDim)]);          FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
107                  if (!Dudley_noError())
108                          break;
109              }          Dudley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
110              mesh_p->Nodes->globalDegreesOfFreedom[i0]=mesh_p->Nodes->Id[i0];          if (!Dudley_noError())
111              mesh_p->Nodes->Tag[i0]=0;              break;
112           }          for (i0 = 0; i0 < numNodes; i0++)
113         }          {
114         /* elements */              if (1 == numDim)
115         else if(!strncmp(&line[1], "ELM", 3) ||              {
116          !strncmp(&line[1], "Elements", 8)) {              scan_ret = fscanf(fileHandle_p, "%d %le %le %le\n", &mesh_p->Nodes->Id[i0],
117                            &mesh_p->Nodes->Coordinates[INDEX2(0, i0, numDim)], &rtmp0, &rtmp1);
118           ElementTypeId final_element_type = NoType;              FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
119           ElementTypeId final_face_element_type = NoType;              }
120           numElements=0;              else if (2 == numDim)
121           numFaceElements=0;              {
122           fscanf(fileHandle_p, "%d", &totalNumElements);              scan_ret = fscanf(fileHandle_p, "%d %le %le %le\n", &mesh_p->Nodes->Id[i0],
123                            &mesh_p->Nodes->Coordinates[INDEX2(0, i0, numDim)],
124           id=TMPMEMALLOC(totalNumElements,index_t);                        &mesh_p->Nodes->Coordinates[INDEX2(1, i0, numDim)], &rtmp0);
125           tag=TMPMEMALLOC(totalNumElements,index_t);              FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
126                  }
127                  else if (3 == numDim)
128           element_type=TMPMEMALLOC(totalNumElements,ElementTypeId);              {
129           vertices=TMPMEMALLOC(totalNumElements*MAX_numNodes_gmsh,index_t);              scan_ret = fscanf(fileHandle_p, "%d %le %le %le\n", &mesh_p->Nodes->Id[i0],
130           if (! (Finley_checkPtr(id) || Finley_checkPtr(tag) || Finley_checkPtr(element_type) || Finley_checkPtr(vertices) ) ) {                        &mesh_p->Nodes->Coordinates[INDEX2(0, i0, numDim)],
131              /* read all in */                        &mesh_p->Nodes->Coordinates[INDEX2(1, i0, numDim)],
132              for(e = 0; e < totalNumElements; e++) {                        &mesh_p->Nodes->Coordinates[INDEX2(2, i0, numDim)]);
133                fscanf(fileHandle_p, "%d %d", &id[e], &gmsh_type);              FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
134                switch (gmsh_type) {              }
135                    case 1:  /* line order 1 */              mesh_p->Nodes->globalDegreesOfFreedom[i0] = mesh_p->Nodes->Id[i0];
136                        element_type[e]=Line2;              mesh_p->Nodes->Tag[i0] = 0;
137                        element_dim=1;          }
138                        numNodesPerElement=2;            }
139                        break;          /* elements */
140                    case 2:  /* traingle order 1 */          else if (!strncmp(&line[1], "ELM", 3) || !strncmp(&line[1], "Elements", 8))
141                        element_type[e]=Tri3;          {
142                        numNodesPerElement= 3;  
143                        element_dim=2;          Dudley_ElementTypeId final_element_type = Dudley_NoRef;
144                        break;          Dudley_ElementTypeId final_face_element_type = Dudley_NoRef;
145                    case 3:  /* quadrilateral order 1 */          numElements = 0;
146                        element_type[e]=Rec4;          numFaceElements = 0;
147                        numNodesPerElement= 4;          scan_ret = fscanf(fileHandle_p, "%d", &totalNumElements);
148                        element_dim=2;          FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
149                        break;  
150                    case 4:  /* tetrahedron order 1 */          id = TMPMEMALLOC(totalNumElements, index_t);
151                        element_type[e]=Tet4;          tag = TMPMEMALLOC(totalNumElements, index_t);
152                        numNodesPerElement= 4;  
153                        element_dim=3;          element_type = TMPMEMALLOC(totalNumElements, Dudley_ElementTypeId);
154                        break;          vertices = TMPMEMALLOC(totalNumElements * MAX_numNodes_gmsh, index_t);
155                    case 5:  /* hexahedron order 1 */          if (!
156                        element_type[e]=Hex8;              (Dudley_checkPtr(id) || Dudley_checkPtr(tag) || Dudley_checkPtr(element_type)
157                        numNodesPerElement= 8;               || Dudley_checkPtr(vertices)))
158                        element_dim=3;          {
159                        break;              /* read all in */
160                    case 8:  /* line order 2 */              for (e = 0; e < totalNumElements; e++)
161                        element_type[e]=Line3;              {
162                        numNodesPerElement= 3;              scan_ret = fscanf(fileHandle_p, "%d %d", &id[e], &gmsh_type);
163                        element_dim=1;              FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
164                        break;              switch (gmsh_type)
165                    case 9:  /* traingle order 2 */              {
166                        element_type[e]=Tri6;              case 1: /* line order 1 */
167                        numNodesPerElement= 6;                  element_type[e] = Dudley_Line2;
168                        element_dim=2;                  element_dim = 1;
169                        break;                  numNodesPerElement = 2;
170                    case 10:  /* quadrilateral order 2 */                  break;
171                        element_type[e]=Rec9;              case 2: /* triangle order 1 */
172                        numNodesPerElement= 9;                  element_type[e] = Dudley_Tri3;
173                        element_dim=2;                  numNodesPerElement = 3;
174                        break;                  element_dim = 2;
175                    case 11:  /* tetrahedron order 2 */                  break;
176                        element_type[e]=Tet10;              case 4: /* tetrahedron order 1 */
177                        numNodesPerElement= 10;                  element_type[e] = Dudley_Tet4;
178                        element_dim=3;                  numNodesPerElement = 4;
179                        break;                  element_dim = 3;
180                    case 15 :  /* point */                  break;
181                        element_type[e]=Point1;              case 15:    /* point */
182                        numNodesPerElement= 1;                  element_type[e] = Dudley_Point1;
183                        element_dim=0;                  numNodesPerElement = 1;
184                        break;                  element_dim = 0;
185                    default:                  break;
186                       element_type[e]=NoType;              default:
187                       sprintf(error_msg,"Unexected gmsh element type %d in mesh file %s.",gmsh_type,fname);                  element_type[e] = Dudley_NoRef;
188                       Finley_setError(IO_ERROR,error_msg);                  sprintf(error_msg, "Unexected gmsh element type %d in mesh file %s.", gmsh_type, fname);
189                }                  Dudley_setError(IO_ERROR, error_msg);
190                if (element_dim == numDim) {              }
191                   if (final_element_type == NoType) {              if (element_dim == numDim)
192                      final_element_type = element_type[e];              {
193                   } else if (final_element_type != element_type[e]) {                  if (final_element_type == Dudley_NoRef)
194                       sprintf(error_msg,"Finley can handle a single type of internal elements only.");                  {
195                       Finley_setError(IO_ERROR,error_msg);                  final_element_type = element_type[e];
196                       break;                  }
197                   }                  else if (final_element_type != element_type[e])
198                   numElements++;                  {
199                } else if (element_dim == numDim-1) {                  sprintf(error_msg, "Dudley can handle a single type of internal elements only.");
200                   if (final_face_element_type == NoType) {                  Dudley_setError(IO_ERROR, error_msg);
201                      final_face_element_type = element_type[e];                  break;
202                   } else if (final_face_element_type != element_type[e]) {                  }
203                       sprintf(error_msg,"Finley can handle a single type of face elements only.");                  numElements++;
204                       Finley_setError(IO_ERROR,error_msg);              }
205                       break;              else if (element_dim == numDim - 1)
206                   }              {
207                   numFaceElements++;                  if (final_face_element_type == Dudley_NoRef)
208                }                  {
209                                  final_face_element_type = element_type[e];
210         if(version <= 1.0){                  }
211           fscanf(fileHandle_p, "%d %d %d", &tag[e], &elementary_id, &numNodesPerElement2);                  else if (final_face_element_type != element_type[e])
212           partition_id = 1;                  {
213                  if (numNodesPerElement2 != numNodesPerElement) {                  sprintf(error_msg, "Dudley can handle a single type of face elements only.");
214                       sprintf(error_msg,"Illegal number of nodes for element %d in mesh file %s.",id[e],fname);                  Dudley_setError(IO_ERROR, error_msg);
215                       Finley_setError(IO_ERROR,error_msg);                  break;
216                  }                  }
217         } else {                  numFaceElements++;
218           fscanf(fileHandle_p, "%d", &numTags);              }
219           elementary_id = tag[e] = partition_id = 1;  
220                  numNodesPerElement2=-1;              if (version <= 1.0)
221           for(j = 0; j < numTags; j++){              {
222             fscanf(fileHandle_p, "%d", &itmp);                        scan_ret = fscanf(fileHandle_p, "%d %d %d", &tag[e], &elementary_id, &numNodesPerElement2);
223             if (j == 0) {                  FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
224               tag[e] = itmp;                  partition_id = 1;
225             } else if (j == 1) {                  if (numNodesPerElement2 != numNodesPerElement)
226               elementary_id = itmp;                  {
227             } else if (j == 2) {                  sprintf(error_msg, "Illegal number of nodes for element %d in mesh file %s.", id[e],
228               partition_id = itmp;                      fname);
229                    }                  Dudley_setError(IO_ERROR, error_msg);
230             /* ignore any other tags */                  }
231           }              }
232         }              else
233                if (! Finley_noError()) break;              {
234                for(j = 0; j < numNodesPerElement; j++) fscanf(fileHandle_p, "%d", &vertices[INDEX2(j,e,MAX_numNodes_gmsh)]);                  scan_ret = fscanf(fileHandle_p, "%d", &numTags);
235                /* for tet10 the last two nodes need to be swapped */                  FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
236                if (element_type[e]==Tet10) {                  elementary_id = tag[e] = partition_id = 1;
237                     itmp=vertices[INDEX2(9,e,MAX_numNodes_gmsh)];                  numNodesPerElement2 = -1;
238                     vertices[INDEX2(9,e,MAX_numNodes_gmsh)]=vertices[INDEX2(8,e,MAX_numNodes_gmsh)];                  for (j = 0; j < numTags; j++)
239                     vertices[INDEX2(8,e,MAX_numNodes_gmsh)]=itmp;                  {
240                }                  scan_ret = fscanf(fileHandle_p, "%d", &itmp);
241              }                  FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
242              /* all elements have been read, now we have to identify the elements for finley */                  if (j == 0)
243                            {
244              if (Finley_noError()) {                      tag[e] = itmp;
245                /* first we have to identify the elements to define Elementis and FaceElements */                  }
246                if (final_element_type == NoType) {                  else if (j == 1)
247                   if (numDim==1) {                  {
248                      final_element_type=Line2;                      elementary_id = itmp;
249                   } else if (numDim==2) {                  }
250                      final_element_type=Tri3;                  else if (j == 2)
251                   } else if (numDim==3) {                  {
252                      final_element_type=Tet4;                      partition_id = itmp;
253                   }                  }
254                }                  /* ignore any other tags */
255                if (final_face_element_type == NoType) {                  }
256                   if (numDim==1) {              }
257                      final_face_element_type=Point1;              if (!Dudley_noError())
258                   } else if (numDim==2) {                  break;
259                      final_face_element_type=Line2;              for (j = 0; j < numNodesPerElement; j++)
260                   } else if (numDim==3) {              {
261                      final_face_element_type=Tri3;                  scan_ret = fscanf(fileHandle_p, "%d", &vertices[INDEX2(j, e, MAX_numNodes_gmsh)]);
262                   }                  FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
263                }              }
264                mesh_p->Elements=Finley_ElementFile_alloc(final_element_type,mesh_p->order, mesh_p->reduced_order, mpi_info);              }
265                mesh_p->FaceElements=Finley_ElementFile_alloc(final_face_element_type,mesh_p->order, mesh_p->reduced_order, mpi_info);              /* all elements have been read, now we have to identify the elements for dudley */
266                mesh_p->ContactElements=Finley_ElementFile_alloc(Point1_Contact,mesh_p->order, mesh_p->reduced_order, mpi_info);  
267                mesh_p->Points=Finley_ElementFile_alloc(Point1,mesh_p->order, mesh_p->reduced_order, mpi_info);              if (Dudley_noError())
268                if (Finley_noError()) {              {
269                    Finley_ElementFile_allocTable(mesh_p->Elements, numElements);              /* first we have to identify the elements to define Elementis and FaceElements */
270                    Finley_ElementFile_allocTable(mesh_p->FaceElements, numFaceElements);              if (final_element_type == Dudley_NoRef)
271                    Finley_ElementFile_allocTable(mesh_p->ContactElements, 0);              {
272                    Finley_ElementFile_allocTable(mesh_p->Points, 0);                  if (numDim == 1)
273                    if (Finley_noError()) {                  {
274                        mesh_p->Elements->minColor=0;                  final_element_type = Dudley_Line2;
275                        mesh_p->Elements->maxColor=numElements-1;                  }
276                        mesh_p->FaceElements->minColor=0;                  else if (numDim == 2)
277                        mesh_p->FaceElements->maxColor=numFaceElements-1;                  {
278                        mesh_p->ContactElements->minColor=0;                  final_element_type = Dudley_Tri3;
279                        mesh_p->ContactElements->maxColor=0;                  }
280                        mesh_p->Points->minColor=0;                  else if (numDim == 3)
281                        mesh_p->Points->maxColor=0;                  {
282                        numElements=0;                  final_element_type = Dudley_Tet4;
283                        numFaceElements=0;                  }
284                        for(e = 0; e < totalNumElements; e++) {              }
285                           if (element_type[e] == final_element_type) {              if (final_face_element_type == Dudley_NoRef)
286                              mesh_p->Elements->Id[numElements]=id[e];              {
287                              mesh_p->Elements->Tag[numElements]=tag[e];                  if (numDim == 1)
288                              mesh_p->Elements->Color[numElements]=numElements;                  {
289                              for (j = 0; j<  mesh_p->Elements->ReferenceElement->Type->numNodes; ++j)  {                  final_face_element_type = Dudley_Point1;
290                                    mesh_p->Elements->Nodes[INDEX2(j, numElements, mesh_p->Elements->ReferenceElement->Type->numNodes)]=vertices[INDEX2(j,e,MAX_numNodes_gmsh)];                  }
291                              }                  else if (numDim == 2)
292                              numElements++;                  {
293                           } else if (element_type[e] == final_face_element_type) {                  final_face_element_type = Dudley_Line2;
294                              mesh_p->FaceElements->Id[numFaceElements]=id[e];                  }
295                              mesh_p->FaceElements->Tag[numFaceElements]=tag[e];                  else if (numDim == 3)
296                              mesh_p->FaceElements->Color[numFaceElements]=numFaceElements;                  {
297                              for (j = 0; j<  mesh_p->FaceElements->ReferenceElement->Type->numNodes; ++j) {                  final_face_element_type = Dudley_Tri3;
298                                       mesh_p->FaceElements->Nodes[INDEX2(j, numFaceElements, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]=vertices[INDEX2(j,e,MAX_numNodes_gmsh)];                  }
299                              }              }
300                              numFaceElements++;              mesh_p->Elements = Dudley_ElementFile_alloc(final_element_type, mpi_info);
301                           }              mesh_p->FaceElements = Dudley_ElementFile_alloc(final_face_element_type, mpi_info);
302                        }              mesh_p->Points = Dudley_ElementFile_alloc(Dudley_Point1, mpi_info);
303                   }              if (Dudley_noError())
304                }              {
305              }                  Dudley_ElementFile_allocTable(mesh_p->Elements, numElements);
306           }                  Dudley_ElementFile_allocTable(mesh_p->FaceElements, numFaceElements);
307           /* and clean up */                  Dudley_ElementFile_allocTable(mesh_p->Points, 0);
308           TMPMEMFREE(id);                  if (Dudley_noError())
309           TMPMEMFREE(tag);                  {
310           TMPMEMFREE(element_type);                  mesh_p->Elements->minColor = 0;
311           TMPMEMFREE(vertices);                  mesh_p->Elements->maxColor = numElements - 1;
312        }                  mesh_p->FaceElements->minColor = 0;
313        /* serach for end of data block */                  mesh_p->FaceElements->maxColor = numFaceElements - 1;
314        do {                  mesh_p->Points->minColor = 0;
315           if (!fgets(line, sizeof(line), fileHandle_p)) {                  mesh_p->Points->maxColor = 0;
316              sprintf(error_msg,"Unexected end of file in %s",fname);                  numElements = 0;
317              Finley_setError(IO_ERROR,error_msg);                  numFaceElements = 0;
318           }                  for (e = 0; e < totalNumElements; e++)
319           if (feof(fileHandle_p)) {                  {
320              sprintf(error_msg,"Unexected end of file in %s",fname);                      if (element_type[e] == final_element_type)
321              Finley_setError(IO_ERROR,error_msg);                      {
322           }                      mesh_p->Elements->Id[numElements] = id[e];
323           if (! Finley_noError()) break;                      mesh_p->Elements->Tag[numElements] = tag[e];
324         } while(line[0] != '$');                      mesh_p->Elements->Color[numElements] = numElements;
325       }                      mesh_p->Elements->Owner[numElements] = 0;
326                          for (j = 0; j < mesh_p->Elements-> /*referenceElementSet-> */ numNodes; ++j)
327       /* close file */                      {
328       fclose(fileHandle_p);                          mesh_p->Elements->Nodes[INDEX2
329       /* clean up */                                      (j, numElements,
330       if (! Finley_noError()) {                                       mesh_p->
331          Finley_Mesh_free(mesh_p);                                       Elements-> /*referenceElementSet-> */ numNodes)] =
332          return NULL;                          vertices[INDEX2(j, e, MAX_numNodes_gmsh)];
333       }                      }
334       /*   resolve id's : */                      numElements++;
335       Finley_Mesh_resolveNodeIds(mesh_p);                      }
336       /* rearrange elements: */                      else if (element_type[e] == final_face_element_type)
337       Finley_Mesh_prepare(mesh_p, optimize);                      {
338       /* that's it */                      mesh_p->FaceElements->Id[numFaceElements] = id[e];
339       #ifdef Finley_TRACE                      mesh_p->FaceElements->Tag[numFaceElements] = tag[e];
340       printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);                      mesh_p->FaceElements->Color[numFaceElements] = numFaceElements;
341       #endif                      mesh_p->FaceElements->Owner[numFaceElements] = 0;
342       Paso_MPIInfo_free( mpi_info );                      for (j = 0; j < mesh_p->FaceElements-> /*referenceElementSet-> */ numNodes; ++j)
343       return mesh_p;                      {
344    }                          mesh_p->FaceElements->Nodes[INDEX2
345                                        (j, numFaceElements,
346                                         mesh_p->
347                                         FaceElements-> /*referenceElementSet-> */
348                                         numNodes)] =
349                            vertices[INDEX2(j, e, MAX_numNodes_gmsh)];
350                        }
351                        numFaceElements++;
352                        }
353                    }
354                    }
355                }
356                }
357            }
358            /* and clean up */
359            TMPMEMFREE(id);
360            TMPMEMFREE(tag);
361            TMPMEMFREE(element_type);
362            TMPMEMFREE(vertices);
363            }
364            /* serach for end of data block */
365            do
366            {
367            if (!fgets(line, sizeof(line), fileHandle_p))
368            {
369                sprintf(error_msg, "Unexected end of file in %s", fname);
370                Dudley_setError(IO_ERROR, error_msg);
371            }
372            if (feof(fileHandle_p))
373            {
374                sprintf(error_msg, "Unexected end of file in %s", fname);
375                Dudley_setError(IO_ERROR, error_msg);
376            }
377            if (!Dudley_noError())
378                break;
379            }
380            while (line[0] != '$');
381        }
382    
383        /* close file */
384        fclose(fileHandle_p);
385        /* clean up */
386        if (!Dudley_noError())
387        {
388            Dudley_Mesh_free(mesh_p);
389            return NULL;
390        }
391        /*   resolve id's : */
392        if (Dudley_noError())
393            Dudley_Mesh_resolveNodeIds(mesh_p);
394        /* rearrange elements: */
395        if (Dudley_noError())
396            Dudley_Mesh_prepare(mesh_p, optimize);
397        /* free up memory */
398        Esys_MPIInfo_free(mpi_info);
399        if (!Dudley_noError())
400        {
401            Dudley_Mesh_free(mesh_p);
402            return NULL;
403        }
404        else
405        {
406            return mesh_p;
407        }
408        }
409  }  }

Legend:
Removed from v.1384  
changed lines
  Added in v.3259

  ViewVC Help
Powered by ViewVC 1.1.26