/[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

revision 3221 by jfenwick, Wed Sep 29 01:00:21 2010 UTC revision 3224 by jfenwick, Wed Sep 29 05:19:37 2010 UTC
# Line 11  Line 11 
11  *  *
12  *******************************************************/  *******************************************************/
13    
   
14  /**************************************************************/  /**************************************************************/
15    
16  /*   Dudley: read mesh */  /*   Dudley: read mesh */
# Line 29  Line 28 
28    
29  #define MAX_numNodes_gmsh 20  #define MAX_numNodes_gmsh 20
30    
31  Dudley_Mesh* Dudley_Mesh_readGmsh(char* fname ,index_t numDim, index_t order, index_t reduced_order, bool_t optimize, bool_t useMacroElements) {  Dudley_Mesh *Dudley_Mesh_readGmsh(char *fname, index_t numDim, index_t order, index_t reduced_order, bool_t optimize,
32                      bool_t useMacroElements)
33    double version = 1.0;  {
34    int format = 0, size = sizeof(double), scan_ret;  
35    dim_t numNodes, totalNumElements=0, numTags=0, numNodesPerElement=0, numNodesPerElement2, element_dim=0;      double version = 1.0;
36    index_t e, i0, j, gmsh_type, partition_id, itmp, elementary_id;      int format = 0, size = sizeof(double), scan_ret;
37    index_t numElements=0, numFaceElements=0, *id=NULL, *tag=NULL, *vertices=NULL;      dim_t numNodes, totalNumElements = 0, numTags = 0, numNodesPerElement = 0, numNodesPerElement2, element_dim = 0;
38    Dudley_Mesh *mesh_p=NULL;      index_t e, i0, j, gmsh_type, partition_id, itmp, elementary_id;
39    char line[LenString_MAX+1];      index_t numElements = 0, numFaceElements = 0, *id = NULL, *tag = NULL, *vertices = NULL;
40    char error_msg[LenErrorMsg_MAX];      Dudley_Mesh *mesh_p = NULL;
41    double rtmp0, rtmp1;      char line[LenString_MAX + 1];
42        char error_msg[LenErrorMsg_MAX];
43        double rtmp0, rtmp1;
44  #ifdef Dudley_TRACE  #ifdef Dudley_TRACE
45    double time0=Dudley_timer();      double time0 = Dudley_timer();
46  #endif  #endif
47    FILE * fileHandle_p = NULL;      FILE *fileHandle_p = NULL;
48    ElementTypeId* element_type=NULL;      ElementTypeId *element_type = NULL;
   
49    
50    Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );      Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc(MPI_COMM_WORLD);
51    Dudley_resetError();      Dudley_resetError();
52    if (mpi_info->size>1) {      if (mpi_info->size > 1)
53      Dudley_setError(IO_ERROR,"reading GMSH with MPI is not supported yet.");      {
54      Paso_MPIInfo_free( mpi_info );      Dudley_setError(IO_ERROR, "reading GMSH with MPI is not supported yet.");
55      return NULL;      Paso_MPIInfo_free(mpi_info);
56    } else {      return NULL;
57        }
58       /* allocate mesh */      else
59          {
60       mesh_p = Dudley_Mesh_alloc(fname,numDim, mpi_info);  
61       if (! Dudley_noError()) return NULL;      /* allocate mesh */
62      
63       /* get file handle */      mesh_p = Dudley_Mesh_alloc(fname, numDim, mpi_info);
64       fileHandle_p = fopen(fname, "r");      if (!Dudley_noError())
65       if (fileHandle_p==NULL) {          return NULL;
66         sprintf(error_msg,"Opening Gmsh file %s for reading failed.",fname);  
67         Dudley_setError(IO_ERROR,error_msg);      /* get file handle */
68         Paso_MPIInfo_free( mpi_info );      fileHandle_p = fopen(fname, "r");
69         return NULL;      if (fileHandle_p == NULL)
70       }      {
71              sprintf(error_msg, "Opening Gmsh file %s for reading failed.", fname);
72       /* start reading */          Dudley_setError(IO_ERROR, error_msg);
73       while(1) {          Paso_MPIInfo_free(mpi_info);
74         if (! Dudley_noError()) break;          return NULL;
75         /* find line staring with $ */      }
76         do {  
77           if( ! fgets(line, sizeof(line), fileHandle_p) ) break;      /* start reading */
78           if(feof(fileHandle_p)) break;      while (1)
79         } while(line[0] != '$');      {
80              if (!Dudley_noError())
81         if (feof(fileHandle_p)) break;          break;
82              /* find line staring with $ */
83              do
84         /* format */          {
85         if (!strncmp(&line[1], "MeshFormat", 10)) {          if (!fgets(line, sizeof(line), fileHandle_p))
86           scan_ret = fscanf(fileHandle_p, "%lf %d %d\n", &version, &format, &size);              break;
87       FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");          if (feof(fileHandle_p))
88         }              break;
89         /* nodes are read */          }
90         if ( !strncmp(&line[1], "NOD", 3)   ||          while (line[0] != '$');
91              !strncmp(&line[1], "NOE", 3)   ||  
92              !strncmp(&line[1], "Nodes", 5)    ) {          if (feof(fileHandle_p))
93                    break;
94           scan_ret = fscanf(fileHandle_p, "%d", &numNodes);  
95       FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");          /* format */
96           if (! Dudley_noError()) break;          if (!strncmp(&line[1], "MeshFormat", 10))
97           Dudley_NodeFile_allocTable(mesh_p->Nodes, numNodes);          {
98           if (! Dudley_noError()) break;          scan_ret = fscanf(fileHandle_p, "%lf %d %d\n", &version, &format, &size);
99           for (i0 = 0; i0 < numNodes; i0++) {          FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
100              if (1 == numDim) {          }
101             scan_ret = fscanf(fileHandle_p, "%d %le %le %le\n", &mesh_p->Nodes->Id[i0],          /* nodes are read */
102                    &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],          if (!strncmp(&line[1], "NOD", 3) || !strncmp(&line[1], "NOE", 3) || !strncmp(&line[1], "Nodes", 5))
103                           &rtmp0,          {
104                           &rtmp1);  
105             FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");          scan_ret = fscanf(fileHandle_p, "%d", &numNodes);
106              } else if (2 == numDim) {          FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
107             scan_ret = fscanf(fileHandle_p, "%d %le %le %le\n", &mesh_p->Nodes->Id[i0],          if (!Dudley_noError())
108                    &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],              break;
109                    &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)],          Dudley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
110                           &rtmp0);          if (!Dudley_noError())
111             FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");              break;
112              } else if (3 == numDim) {          for (i0 = 0; i0 < numNodes; i0++)
113             scan_ret = fscanf(fileHandle_p, "%d %le %le %le\n", &mesh_p->Nodes->Id[i0],          {
114                    &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],              if (1 == numDim)
115                    &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)],              {
116                    &mesh_p->Nodes->Coordinates[INDEX2(2,i0,numDim)]);              scan_ret = fscanf(fileHandle_p, "%d %le %le %le\n", &mesh_p->Nodes->Id[i0],
117             FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");                        &mesh_p->Nodes->Coordinates[INDEX2(0, i0, numDim)], &rtmp0, &rtmp1);
118              }              FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
119              mesh_p->Nodes->globalDegreesOfFreedom[i0]=mesh_p->Nodes->Id[i0];              }
120              mesh_p->Nodes->Tag[i0]=0;              else if (2 == numDim)
121           }              {
122         }              scan_ret = fscanf(fileHandle_p, "%d %le %le %le\n", &mesh_p->Nodes->Id[i0],
123         /* elements */                        &mesh_p->Nodes->Coordinates[INDEX2(0, i0, numDim)],
124         else if(!strncmp(&line[1], "ELM", 3) ||                        &mesh_p->Nodes->Coordinates[INDEX2(1, i0, numDim)], &rtmp0);
125          !strncmp(&line[1], "Elements", 8)) {              FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
126                  }
127           ElementTypeId final_element_type = NoRef;              else if (3 == numDim)
128           ElementTypeId final_face_element_type = NoRef;              {
129           numElements=0;              scan_ret = fscanf(fileHandle_p, "%d %le %le %le\n", &mesh_p->Nodes->Id[i0],
130           numFaceElements=0;                        &mesh_p->Nodes->Coordinates[INDEX2(0, i0, numDim)],
131           scan_ret = fscanf(fileHandle_p, "%d", &totalNumElements);                        &mesh_p->Nodes->Coordinates[INDEX2(1, i0, numDim)],
132       FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");                        &mesh_p->Nodes->Coordinates[INDEX2(2, i0, numDim)]);
133                  FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
134           id=TMPMEMALLOC(totalNumElements,index_t);              }
135           tag=TMPMEMALLOC(totalNumElements,index_t);              mesh_p->Nodes->globalDegreesOfFreedom[i0] = mesh_p->Nodes->Id[i0];
136                  mesh_p->Nodes->Tag[i0] = 0;
137              }
138           element_type=TMPMEMALLOC(totalNumElements,ElementTypeId);          }
139           vertices=TMPMEMALLOC(totalNumElements*MAX_numNodes_gmsh,index_t);          /* elements */
140           if (! (Dudley_checkPtr(id) || Dudley_checkPtr(tag) || Dudley_checkPtr(element_type) || Dudley_checkPtr(vertices) ) ) {          else if (!strncmp(&line[1], "ELM", 3) || !strncmp(&line[1], "Elements", 8))
141              /* read all in */          {
142              for(e = 0; e < totalNumElements; e++) {  
143                scan_ret = fscanf(fileHandle_p, "%d %d", &id[e], &gmsh_type);          ElementTypeId final_element_type = NoRef;
144            FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");          ElementTypeId final_face_element_type = NoRef;
145                switch (gmsh_type) {          numElements = 0;
146                    case 1:  /* line order 1 */          numFaceElements = 0;
147                        element_type[e]=Line2;          scan_ret = fscanf(fileHandle_p, "%d", &totalNumElements);
148                        element_dim=1;          FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
149                        numNodesPerElement=2;    
150                        break;          id = TMPMEMALLOC(totalNumElements, index_t);
151                    case 2:  /* traingle order 1 */          tag = TMPMEMALLOC(totalNumElements, index_t);
152                        element_type[e]=Tri3;  
153                        numNodesPerElement= 3;          element_type = TMPMEMALLOC(totalNumElements, ElementTypeId);
154                        element_dim=2;          vertices = TMPMEMALLOC(totalNumElements * MAX_numNodes_gmsh, index_t);
155                        break;          if (!
156                    case 4:  /* tetrahedron order 1 */              (Dudley_checkPtr(id) || Dudley_checkPtr(tag) || Dudley_checkPtr(element_type)
157                        element_type[e]=Tet4;               || Dudley_checkPtr(vertices)))
158                        numNodesPerElement= 4;          {
159                        element_dim=3;              /* read all in */
160                        break;              for (e = 0; e < totalNumElements; e++)
161                    case 15 :  /* point */              {
162                        element_type[e]=Point1;              scan_ret = fscanf(fileHandle_p, "%d %d", &id[e], &gmsh_type);
163                        numNodesPerElement= 1;              FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
164                        element_dim=0;              switch (gmsh_type)
165                        break;              {
166                    default:              case 1: /* line order 1 */
167                       element_type[e]=NoRef;                  element_type[e] = Line2;
168                       sprintf(error_msg,"Unexected gmsh element type %d in mesh file %s.",gmsh_type,fname);                  element_dim = 1;
169                       Dudley_setError(IO_ERROR,error_msg);                  numNodesPerElement = 2;
170                }                  break;
171                if (element_dim == numDim) {              case 2: /* traingle order 1 */
172                   if (final_element_type == NoRef) {                  element_type[e] = Tri3;
173                      final_element_type = element_type[e];                  numNodesPerElement = 3;
174                   } else if (final_element_type != element_type[e]) {                  element_dim = 2;
175                       sprintf(error_msg,"Dudley can handle a single type of internal elements only.");                  break;
176                       Dudley_setError(IO_ERROR,error_msg);              case 4: /* tetrahedron order 1 */
177                       break;                  element_type[e] = Tet4;
178                   }                  numNodesPerElement = 4;
179                   numElements++;                  element_dim = 3;
180                } else if (element_dim == numDim-1) {                  break;
181                   if (final_face_element_type == NoRef) {              case 15:    /* point */
182                      final_face_element_type = element_type[e];                  element_type[e] = Point1;
183                   } else if (final_face_element_type != element_type[e]) {                  numNodesPerElement = 1;
184                       sprintf(error_msg,"Dudley can handle a single type of face elements only.");                  element_dim = 0;
185                       Dudley_setError(IO_ERROR,error_msg);                  break;
186                       break;              default:
187                   }                  element_type[e] = NoRef;
188                   numFaceElements++;                  sprintf(error_msg, "Unexected gmsh element type %d in mesh file %s.", gmsh_type, fname);
189                }                  Dudley_setError(IO_ERROR, error_msg);
190                              }
191         if(version <= 1.0){              if (element_dim == numDim)
192           scan_ret = fscanf(fileHandle_p, "%d %d %d", &tag[e], &elementary_id, &numNodesPerElement2);              {
193           FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");                  if (final_element_type == NoRef)
194           partition_id = 1;                  {
195                  if (numNodesPerElement2 != numNodesPerElement) {                  final_element_type = element_type[e];
196                       sprintf(error_msg,"Illegal number of nodes for element %d in mesh file %s.",id[e],fname);                  }
197                       Dudley_setError(IO_ERROR,error_msg);                  else if (final_element_type != element_type[e])
198                  }                  {
199         } else {                  sprintf(error_msg, "Dudley can handle a single type of internal elements only.");
200           scan_ret = fscanf(fileHandle_p, "%d", &numTags);                  Dudley_setError(IO_ERROR, error_msg);
201           FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");                  break;
202           elementary_id = tag[e] = partition_id = 1;                  }
203                  numNodesPerElement2=-1;                  numElements++;
204           for(j = 0; j < numTags; j++){              }
205             scan_ret = fscanf(fileHandle_p, "%d", &itmp);                      else if (element_dim == numDim - 1)
206             FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");              {
207             if (j == 0) {                  if (final_face_element_type == NoRef)
208               tag[e] = itmp;                  {
209             } else if (j == 1) {                  final_face_element_type = element_type[e];
210               elementary_id = itmp;                  }
211             } else if (j == 2) {                  else if (final_face_element_type != element_type[e])
212               partition_id = itmp;                  {
213                    }                  sprintf(error_msg, "Dudley can handle a single type of face elements only.");
214             /* ignore any other tags */                  Dudley_setError(IO_ERROR, error_msg);
215           }                  break;
216         }                  }
217                if (! Dudley_noError()) break;                  numFaceElements++;
218                for(j = 0; j < numNodesPerElement; j++) {              }
219          scan_ret = fscanf(fileHandle_p, "%d", &vertices[INDEX2(j,e,MAX_numNodes_gmsh)]);  
220              FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");              if (version <= 1.0)
221            }              {
222              }                  scan_ret = fscanf(fileHandle_p, "%d %d %d", &tag[e], &elementary_id, &numNodesPerElement2);
223              /* all elements have been read, now we have to identify the elements for dudley */                  FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
224                            partition_id = 1;
225              if (Dudley_noError()) {                  if (numNodesPerElement2 != numNodesPerElement)
226                /* first we have to identify the elements to define Elementis and FaceElements */                  {
227                if (final_element_type == NoRef) {                  sprintf(error_msg, "Illegal number of nodes for element %d in mesh file %s.", id[e],
228                   if (numDim==1) {                      fname);
229                      final_element_type=Line2;                  Dudley_setError(IO_ERROR, error_msg);
230                   } else if (numDim==2) {                  }
231                      final_element_type=Tri3;              }
232                   } else if (numDim==3) {              else
233                      final_element_type=Tet4;              {
234                   }                  scan_ret = fscanf(fileHandle_p, "%d", &numTags);
235                }                  FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
236                if (final_face_element_type == NoRef) {                  elementary_id = tag[e] = partition_id = 1;
237                   if (numDim==1) {                  numNodesPerElement2 = -1;
238                      final_face_element_type=Point1;                  for (j = 0; j < numTags; j++)
239                   } else if (numDim==2) {                  {
240                      final_face_element_type=Line2;                  scan_ret = fscanf(fileHandle_p, "%d", &itmp);
241                   } else if (numDim==3) {                  FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
242                      final_face_element_type=Tri3;                  if (j == 0)
243                   }                  {
244                }                      tag[e] = itmp;
245                mesh_p->Elements=Dudley_ElementFile_alloc(final_element_type, mpi_info);                  }
246                mesh_p->FaceElements=Dudley_ElementFile_alloc(final_face_element_type, mpi_info);                  else if (j == 1)
247                mesh_p->Points=Dudley_ElementFile_alloc(Point1, mpi_info);                  {
248                if (Dudley_noError()) {                      elementary_id = itmp;
249                    Dudley_ElementFile_allocTable(mesh_p->Elements, numElements);                  }
250                    Dudley_ElementFile_allocTable(mesh_p->FaceElements, numFaceElements);                  else if (j == 2)
251                    Dudley_ElementFile_allocTable(mesh_p->Points, 0);                  {
252                    if (Dudley_noError()) {                      partition_id = itmp;
253                        mesh_p->Elements->minColor=0;                  }
254                        mesh_p->Elements->maxColor=numElements-1;                  /* ignore any other tags */
255                        mesh_p->FaceElements->minColor=0;                  }
256                        mesh_p->FaceElements->maxColor=numFaceElements-1;              }
257                        mesh_p->Points->minColor=0;              if (!Dudley_noError())
258                        mesh_p->Points->maxColor=0;                  break;
259                        numElements=0;              for (j = 0; j < numNodesPerElement; j++)
260                        numFaceElements=0;              {
261                        for(e = 0; e < totalNumElements; e++) {                  scan_ret = fscanf(fileHandle_p, "%d", &vertices[INDEX2(j, e, MAX_numNodes_gmsh)]);
262                           if (element_type[e] == final_element_type) {                  FSCANF_CHECK(scan_ret, "fscanf: Dudley_Mesh_readGmsh");
263                              mesh_p->Elements->Id[numElements]=id[e];              }
264                              mesh_p->Elements->Tag[numElements]=tag[e];              }
265                              mesh_p->Elements->Color[numElements]=numElements;              /* all elements have been read, now we have to identify the elements for dudley */
266                              mesh_p->Elements->Owner[numElements]=0;  
267                              for (j = 0; j<  mesh_p->Elements->/*referenceElementSet->*/numNodes; ++j)  {              if (Dudley_noError())
268                                    mesh_p->Elements->Nodes[INDEX2(j, numElements, mesh_p->Elements->/*referenceElementSet->*/numNodes)]=vertices[INDEX2(j,e,MAX_numNodes_gmsh)];              {
269                              }              /* first we have to identify the elements to define Elementis and FaceElements */
270                              numElements++;              if (final_element_type == NoRef)
271                           } else if (element_type[e] == final_face_element_type) {              {
272                              mesh_p->FaceElements->Id[numFaceElements]=id[e];                  if (numDim == 1)
273                              mesh_p->FaceElements->Tag[numFaceElements]=tag[e];                  {
274                              mesh_p->FaceElements->Color[numFaceElements]=numFaceElements;                  final_element_type = Line2;
275                              mesh_p->FaceElements->Owner[numFaceElements]=0;                  }
276                              for (j = 0; j<  mesh_p->FaceElements->/*referenceElementSet->*/numNodes; ++j) {                  else if (numDim == 2)
277                                       mesh_p->FaceElements->Nodes[INDEX2(j, numFaceElements, mesh_p->FaceElements->/*referenceElementSet->*/numNodes)]=vertices[INDEX2(j,e,MAX_numNodes_gmsh)];                  {
278                              }                  final_element_type = Tri3;
279                              numFaceElements++;                  }
280                           }                  else if (numDim == 3)
281                        }                  {
282                   }                  final_element_type = Tet4;
283                }                  }
284              }              }
285           }              if (final_face_element_type == NoRef)
286           /* and clean up */              {
287           TMPMEMFREE(id);                  if (numDim == 1)
288           TMPMEMFREE(tag);                  {
289           TMPMEMFREE(element_type);                  final_face_element_type = Point1;
290           TMPMEMFREE(vertices);                  }
291        }                  else if (numDim == 2)
292        /* serach for end of data block */                  {
293        do {                  final_face_element_type = Line2;
294           if (!fgets(line, sizeof(line), fileHandle_p)) {                  }
295              sprintf(error_msg,"Unexected end of file in %s",fname);                  else if (numDim == 3)
296              Dudley_setError(IO_ERROR,error_msg);                  {
297           }                  final_face_element_type = Tri3;
298           if (feof(fileHandle_p)) {                  }
299              sprintf(error_msg,"Unexected end of file in %s",fname);              }
300              Dudley_setError(IO_ERROR,error_msg);              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           if (! Dudley_noError()) break;              mesh_p->Points = Dudley_ElementFile_alloc(Point1, mpi_info);
303         } while(line[0] != '$');              if (Dudley_noError())
304       }              {
305                      Dudley_ElementFile_allocTable(mesh_p->Elements, numElements);
306       /* close file */                  Dudley_ElementFile_allocTable(mesh_p->FaceElements, numFaceElements);
307       fclose(fileHandle_p);                  Dudley_ElementFile_allocTable(mesh_p->Points, 0);
308       /* clean up */                  if (Dudley_noError())
309       if (! Dudley_noError()) {                  {
310          Dudley_Mesh_free(mesh_p);                  mesh_p->Elements->minColor = 0;
311          return NULL;                  mesh_p->Elements->maxColor = numElements - 1;
312       }                  mesh_p->FaceElements->minColor = 0;
313       /*   resolve id's : */                  mesh_p->FaceElements->maxColor = numFaceElements - 1;
314       if (Dudley_noError()) Dudley_Mesh_resolveNodeIds(mesh_p);                  mesh_p->Points->minColor = 0;
315       /* rearrange elements: */                  mesh_p->Points->maxColor = 0;
316       if (Dudley_noError()) Dudley_Mesh_prepare(mesh_p, optimize);                  numElements = 0;
317       /* free up memory */                  numFaceElements = 0;
318       Paso_MPIInfo_free( mpi_info );                  for (e = 0; e < totalNumElements; e++)
319       if (! Dudley_noError()) {                  {
320          Dudley_Mesh_free(mesh_p);                      if (element_type[e] == final_element_type)
321          return NULL;                      {
322       } else {                      mesh_p->Elements->Id[numElements] = id[e];
323           return mesh_p;                      mesh_p->Elements->Tag[numElements] = tag[e];
324       }                      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                        {
328                            mesh_p->Elements->Nodes[INDEX2
329                                        (j, numElements,
330                                         mesh_p->
331                                         Elements-> /*referenceElementSet-> */ numNodes)] =
332                            vertices[INDEX2(j, e, MAX_numNodes_gmsh)];
333                        }
334                        numElements++;
335                        }
336                        else if (element_type[e] == final_face_element_type)
337                        {
338                        mesh_p->FaceElements->Id[numFaceElements] = id[e];
339                        mesh_p->FaceElements->Tag[numFaceElements] = tag[e];
340                        mesh_p->FaceElements->Color[numFaceElements] = numFaceElements;
341                        mesh_p->FaceElements->Owner[numFaceElements] = 0;
342                        for (j = 0; j < mesh_p->FaceElements-> /*referenceElementSet-> */ numNodes; ++j)
343                        {
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        Paso_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.3221  
changed lines
  Added in v.3224

  ViewVC Help
Powered by ViewVC 1.1.26