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

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

  ViewVC Help
Powered by ViewVC 1.1.26