/[escript]/trunk/finley/src/Mesh_readGmsh.c
ViewVC logotype

Diff of /trunk/finley/src/Mesh_readGmsh.c

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

revision 3258 by gross, Fri Feb 12 05:09:06 2010 UTC revision 3259 by jfenwick, Mon Oct 11 01:48:14 2010 UTC
# Line 45  Finley_Mesh* Finley_Mesh_readGmsh(char* Line 45  Finley_Mesh* Finley_Mesh_readGmsh(char*
45    double time0=Finley_timer();    double time0=Finley_timer();
46  #endif  #endif
47    FILE * fileHandle_p = NULL;    FILE * fileHandle_p = NULL;
48    ElementTypeId* element_type=NULL;    Finley_ElementTypeId* element_type=NULL;
49    
50    
51    Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );    Esys_MPIInfo *mpi_info = Esys_MPIInfo_alloc( MPI_COMM_WORLD );
52    Finley_resetError();    Finley_resetError();
53    if (mpi_info->size>1) {    if (mpi_info->size>1) {
54      Finley_setError(IO_ERROR,"reading GMSH with MPI is not supported yet.");      Finley_setError(IO_ERROR,"reading GMSH with MPI is not supported yet.");
55      Paso_MPIInfo_free( mpi_info );      Esys_MPIInfo_free( mpi_info );
56      return NULL;      return NULL;
57    } else {    } else {
58    
# Line 66  Finley_Mesh* Finley_Mesh_readGmsh(char* Line 66  Finley_Mesh* Finley_Mesh_readGmsh(char*
66       if (fileHandle_p==NULL) {       if (fileHandle_p==NULL) {
67         sprintf(error_msg,"Opening Gmsh file %s for reading failed.",fname);         sprintf(error_msg,"Opening Gmsh file %s for reading failed.",fname);
68         Finley_setError(IO_ERROR,error_msg);         Finley_setError(IO_ERROR,error_msg);
69         Paso_MPIInfo_free( mpi_info );         Esys_MPIInfo_free( mpi_info );
70         return NULL;         return NULL;
71       }       }
72        
# Line 125  Finley_Mesh* Finley_Mesh_readGmsh(char* Line 125  Finley_Mesh* Finley_Mesh_readGmsh(char*
125         else if(!strncmp(&line[1], "ELM", 3) ||         else if(!strncmp(&line[1], "ELM", 3) ||
126          !strncmp(&line[1], "Elements", 8)) {          !strncmp(&line[1], "Elements", 8)) {
127        
128           ElementTypeId final_element_type = NoRef;           Finley_ElementTypeId final_element_type = Finley_NoRef;
129           ElementTypeId final_face_element_type = NoRef;           Finley_ElementTypeId final_face_element_type = Finley_NoRef;
130           ElementTypeId contact_element_type = NoRef;           Finley_ElementTypeId contact_element_type = Finley_NoRef;
131           numElements=0;           numElements=0;
132           numFaceElements=0;           numFaceElements=0;
133           scan_ret = fscanf(fileHandle_p, "%d", &totalNumElements);           scan_ret = fscanf(fileHandle_p, "%d", &totalNumElements);
# Line 137  Finley_Mesh* Finley_Mesh_readGmsh(char* Line 137  Finley_Mesh* Finley_Mesh_readGmsh(char*
137           tag=TMPMEMALLOC(totalNumElements,index_t);           tag=TMPMEMALLOC(totalNumElements,index_t);
138        
139        
140           element_type=TMPMEMALLOC(totalNumElements,ElementTypeId);           element_type=TMPMEMALLOC(totalNumElements,Finley_ElementTypeId);
141           vertices=TMPMEMALLOC(totalNumElements*MAX_numNodes_gmsh,index_t);           vertices=TMPMEMALLOC(totalNumElements*MAX_numNodes_gmsh,index_t);
142           if (! (Finley_checkPtr(id) || Finley_checkPtr(tag) || Finley_checkPtr(element_type) || Finley_checkPtr(vertices) ) ) {           if (! (Finley_checkPtr(id) || Finley_checkPtr(tag) || Finley_checkPtr(element_type) || Finley_checkPtr(vertices) ) ) {
143              /* read all in */              /* read all in */
# Line 146  Finley_Mesh* Finley_Mesh_readGmsh(char* Line 146  Finley_Mesh* Finley_Mesh_readGmsh(char*
146            FSCANF_CHECK(scan_ret, "fscanf: Finley_Mesh_readGmsh");            FSCANF_CHECK(scan_ret, "fscanf: Finley_Mesh_readGmsh");
147                switch (gmsh_type) {                switch (gmsh_type) {
148                    case 1:  /* line order 1 */                    case 1:  /* line order 1 */
149                        element_type[e]=Line2;                        element_type[e]=Finley_Line2;
150                        element_dim=1;                        element_dim=1;
151                        numNodesPerElement=2;                          numNodesPerElement=2;  
152                        break;                        break;
153                    case 2:  /* traingle order 1 */                    case 2:  /* traingle order 1 */
154                        element_type[e]=Tri3;                        element_type[e]=Finley_Tri3;
155                        numNodesPerElement= 3;                        numNodesPerElement= 3;
156                        element_dim=2;                        element_dim=2;
157                        break;                        break;
158                    case 3:  /* quadrilateral order 1 */                    case 3:  /* quadrilateral order 1 */
159                        element_type[e]=Rec4;                        element_type[e]=Finley_Rec4;
160                        numNodesPerElement= 4;                        numNodesPerElement= 4;
161                        element_dim=2;                        element_dim=2;
162                        break;                        break;
163                    case 4:  /* tetrahedron order 1 */                    case 4:  /* tetrahedron order 1 */
164                        element_type[e]=Tet4;                        element_type[e]=Finley_Tet4;
165                        numNodesPerElement= 4;                        numNodesPerElement= 4;
166                        element_dim=3;                        element_dim=3;
167                        break;                        break;
168                    case 5:  /* hexahedron order 1 */                    case 5:  /* hexahedron order 1 */
169                        element_type[e]=Hex8;                        element_type[e]=Finley_Hex8;
170                        numNodesPerElement= 8;                        numNodesPerElement= 8;
171                        element_dim=3;                        element_dim=3;
172                        break;                        break;
173                    case 8:  /* line order 2 */                    case 8:  /* line order 2 */
174                        if (useMacroElements) {                        if (useMacroElements) {
175                            element_type[e]=Line3Macro;                            element_type[e]=Finley_Line3Macro;
176                        } else {                        } else {
177                            element_type[e]=Line3;                            element_type[e]=Finley_Line3;
178                        }                        }
179                        numNodesPerElement= 3;                        numNodesPerElement= 3;
180                        element_dim=1;                        element_dim=1;
181                        break;                        break;
182                    case 9:  /* traingle order 2 */                    case 9:  /* traingle order 2 */
183                        if (useMacroElements) {                        if (useMacroElements) {
184                             element_type[e]=Tri6Macro;                             element_type[e]=Finley_Tri6Macro;
185                        } else {                        } else {
186                             element_type[e]=Tri6;                             element_type[e]=Finley_Tri6;
187                        }                        }
188                        numNodesPerElement= 6;                        numNodesPerElement= 6;
189                        element_dim=2;                        element_dim=2;
190                        break;                        break;
191                    case 10:  /* quadrilateral order 2 */                    case 10:  /* quadrilateral order 2 */
192                        if (useMacroElements) {                        if (useMacroElements) {
193                            element_type[e]=Rec9Macro;                            element_type[e]=Finley_Rec9Macro;
194                        } else {                        } else {
195                            element_type[e]=Rec9;                            element_type[e]=Finley_Rec9;
196                        }                        }
197                        numNodesPerElement= 9;                        numNodesPerElement= 9;
198                        element_dim=2;                        element_dim=2;
199                        break;                        break;
200                    case 11:  /* tetrahedron order 2 */                    case 11:  /* tetrahedron order 2 */
201                        if (useMacroElements) {                        if (useMacroElements) {
202                            element_type[e]=Tet10Macro;                            element_type[e]=Finley_Tet10Macro;
203                        } else {                        } else {
204                            element_type[e]=Tet10;                            element_type[e]=Finley_Tet10;
205                        }                        }
206                        numNodesPerElement= 10;                        numNodesPerElement= 10;
207                        element_dim=3;                        element_dim=3;
208                        break;                        break;
209                    case 16:  /* rectangular order 2 */                    case 16:  /* rectangular order 2 */
210                        element_type[e]=Rec8;                        element_type[e]=Finley_Rec8;
211                        numNodesPerElement= 8;                        numNodesPerElement= 8;
212                        element_dim=2;                        element_dim=2;
213                        break;                        break;
214                    case 17:  /* hexahedron order 2 */                    case 17:  /* hexahedron order 2 */
215                        element_type[e]=Hex20;                        element_type[e]=Finley_Hex20;
216                        numNodesPerElement= 20;                        numNodesPerElement= 20;
217                        element_dim=3;                        element_dim=3;
218                        break;                        break;
219                    case 15 :  /* point */                    case 15 :  /* point */
220                        element_type[e]=Point1;                        element_type[e]=Finley_Point1;
221                        numNodesPerElement= 1;                        numNodesPerElement= 1;
222                        element_dim=0;                        element_dim=0;
223                        break;                        break;
224                    default:                    default:
225                       element_type[e]=NoRef;                       element_type[e]=Finley_NoRef;
226                       sprintf(error_msg,"Unexected gmsh element type %d in mesh file %s.",gmsh_type,fname);                       sprintf(error_msg,"Unexected gmsh element type %d in mesh file %s.",gmsh_type,fname);
227                       Finley_setError(IO_ERROR,error_msg);                       Finley_setError(IO_ERROR,error_msg);
228                }                }
229                if (element_dim == numDim) {                if (element_dim == numDim) {
230                   if (final_element_type == NoRef) {                   if (final_element_type == Finley_NoRef) {
231                      final_element_type = element_type[e];                      final_element_type = element_type[e];
232                   } else if (final_element_type != element_type[e]) {                   } else if (final_element_type != element_type[e]) {
233                       sprintf(error_msg,"Finley can handle a single type of internal elements only.");                       sprintf(error_msg,"Finley can handle a single type of internal elements only.");
# Line 236  Finley_Mesh* Finley_Mesh_readGmsh(char* Line 236  Finley_Mesh* Finley_Mesh_readGmsh(char*
236                   }                   }
237                   numElements++;                   numElements++;
238                } else if (element_dim == numDim-1) {                } else if (element_dim == numDim-1) {
239                   if (final_face_element_type == NoRef) {                   if (final_face_element_type == Finley_NoRef) {
240                      final_face_element_type = element_type[e];                      final_face_element_type = element_type[e];
241                   } else if (final_face_element_type != element_type[e]) {                   } else if (final_face_element_type != element_type[e]) {
242                       sprintf(error_msg,"Finley can handle a single type of face elements only.");                       sprintf(error_msg,"Finley can handle a single type of face elements only.");
# Line 278  Finley_Mesh* Finley_Mesh_readGmsh(char* Line 278  Finley_Mesh* Finley_Mesh_readGmsh(char*
278              FSCANF_CHECK(scan_ret, "fscanf: Finley_Mesh_readGmsh");              FSCANF_CHECK(scan_ret, "fscanf: Finley_Mesh_readGmsh");
279            }            }
280                /* for tet10 the last two nodes need to be swapped */                /* for tet10 the last two nodes need to be swapped */
281                if ((element_type[e]==Tet10) || (element_type[e]==Tet10Macro)) {                if ((element_type[e]==Finley_Tet10) || (element_type[e]==Finley_Tet10Macro)) {
282                     itmp=vertices[INDEX2(9,e,MAX_numNodes_gmsh)];                     itmp=vertices[INDEX2(9,e,MAX_numNodes_gmsh)];
283                     vertices[INDEX2(9,e,MAX_numNodes_gmsh)]=vertices[INDEX2(8,e,MAX_numNodes_gmsh)];                     vertices[INDEX2(9,e,MAX_numNodes_gmsh)]=vertices[INDEX2(8,e,MAX_numNodes_gmsh)];
284                     vertices[INDEX2(8,e,MAX_numNodes_gmsh)]=itmp;                     vertices[INDEX2(8,e,MAX_numNodes_gmsh)]=itmp;
# Line 288  Finley_Mesh* Finley_Mesh_readGmsh(char* Line 288  Finley_Mesh* Finley_Mesh_readGmsh(char*
288                    
289              if (Finley_noError()) {              if (Finley_noError()) {
290                /* first we have to identify the elements to define Elementis and FaceElements */                /* first we have to identify the elements to define Elementis and FaceElements */
291                if (final_element_type == NoRef) {                if (final_element_type == Finley_NoRef) {
292                   if (numDim==1) {                   if (numDim==1) {
293                      final_element_type=Line2;                      final_element_type=Finley_Line2;
294                   } else if (numDim==2) {                   } else if (numDim==2) {
295                      final_element_type=Tri3;                      final_element_type=Finley_Tri3;
296                   } else if (numDim==3) {                   } else if (numDim==3) {
297                      final_element_type=Tet4;                      final_element_type=Finley_Tet4;
298                   }                   }
299                }                }
300                if (final_face_element_type == NoRef) {                if (final_face_element_type == Finley_NoRef) {
301                   if (numDim==1) {                   if (numDim==1) {
302                      final_face_element_type=Point1;                      final_face_element_type=Finley_Point1;
303                   } else if (numDim==2) {                   } else if (numDim==2) {
304                      final_face_element_type=Line2;                      final_face_element_type=Finley_Line2;
305                   } else if (numDim==3) {                   } else if (numDim==3) {
306                      final_face_element_type=Tri3;                      final_face_element_type=Finley_Tri3;
307                   }                   }
308                }                }
309                if (final_face_element_type == Line2) {                if (final_face_element_type == Finley_Line2) {
310                    contact_element_type=Line2_Contact;                    contact_element_type=Finley_Line2_Contact;
311                } else  if ( (final_face_element_type == Line3) || (final_face_element_type == Line3Macro) ) {                } else  if ( (final_face_element_type == Finley_Line3) || (final_face_element_type == Finley_Line3Macro) ) {
312                    contact_element_type=Line3_Contact;                    contact_element_type=Finley_Line3_Contact;
313                } else  if (final_face_element_type == Tri3) {                } else  if (final_face_element_type == Finley_Tri3) {
314                    contact_element_type=Tri3_Contact;                    contact_element_type=Finley_Tri3_Contact;
315                } else  if ( (final_face_element_type == Tri6) || (final_face_element_type == Tri6Macro)) {                } else  if ( (final_face_element_type == Finley_Tri6) || (final_face_element_type == Finley_Tri6Macro)) {
316                    contact_element_type=Tri6_Contact;                    contact_element_type=Finley_Tri6_Contact;
317                } else {                } else {
318                    contact_element_type=Point1_Contact;                    contact_element_type=Finley_Point1_Contact;
319                }                }
320                refElements= Finley_ReferenceElementSet_alloc(final_element_type,order, reduced_order);                refElements= Finley_ReferenceElementSet_alloc(final_element_type,order, reduced_order);
321                refFaceElements=Finley_ReferenceElementSet_alloc(final_face_element_type,order, reduced_order);                refFaceElements=Finley_ReferenceElementSet_alloc(final_face_element_type,order, reduced_order);
322                refContactElements= Finley_ReferenceElementSet_alloc(contact_element_type,order, reduced_order);                refContactElements= Finley_ReferenceElementSet_alloc(contact_element_type,order, reduced_order);
323                refPoints= Finley_ReferenceElementSet_alloc(Point1,order, reduced_order);                refPoints= Finley_ReferenceElementSet_alloc(Finley_Point1,order, reduced_order);
324                mesh_p->Elements=Finley_ElementFile_alloc(refElements, mpi_info);                mesh_p->Elements=Finley_ElementFile_alloc(refElements, mpi_info);
325                mesh_p->FaceElements=Finley_ElementFile_alloc(refFaceElements, mpi_info);                mesh_p->FaceElements=Finley_ElementFile_alloc(refFaceElements, mpi_info);
326                mesh_p->ContactElements=Finley_ElementFile_alloc(refContactElements, mpi_info);                mesh_p->ContactElements=Finley_ElementFile_alloc(refContactElements, mpi_info);
# Line 402  Finley_Mesh* Finley_Mesh_readGmsh(char* Line 402  Finley_Mesh* Finley_Mesh_readGmsh(char*
402       Finley_ReferenceElementSet_dealloc(refContactElements);       Finley_ReferenceElementSet_dealloc(refContactElements);
403       Finley_ReferenceElementSet_dealloc(refFaceElements);       Finley_ReferenceElementSet_dealloc(refFaceElements);
404       Finley_ReferenceElementSet_dealloc(refElements);       Finley_ReferenceElementSet_dealloc(refElements);
405       Paso_MPIInfo_free( mpi_info );       Esys_MPIInfo_free( mpi_info );
406       if (! Finley_noError()) {       if (! Finley_noError()) {
407          Finley_Mesh_free(mesh_p);          Finley_Mesh_free(mesh_p);
408          return NULL;          return NULL;

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

  ViewVC Help
Powered by ViewVC 1.1.26