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

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

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

revision 4491 by caltinay, Fri May 31 07:09:03 2013 UTC revision 4492 by caltinay, Tue Jul 2 01:44:11 2013 UTC
# Line 42  Finley_Mesh* Finley_Mesh_readGmsh(char* Line 42  Finley_Mesh* Finley_Mesh_readGmsh(char*
42    char line[LenString_MAX+1], name[LenString_MAX+1];    char line[LenString_MAX+1], name[LenString_MAX+1];
43    char error_msg[LenErrorMsg_MAX];    char error_msg[LenErrorMsg_MAX];
44    double rtmp0, rtmp1;    double rtmp0, rtmp1;
45    Finley_ReferenceElementSet *refPoints=NULL, *refContactElements=NULL, *refFaceElements=NULL, *refElements=NULL;    ReferenceElementSet *refPoints=NULL, *refContactElements=NULL, *refFaceElements=NULL, *refElements=NULL;
46  #ifdef Finley_TRACE  #ifdef Finley_TRACE
47    double time0=Finley_timer();    double time0=Finley_timer();
48  #endif  #endif
49    FILE * fileHandle_p = NULL;    FILE * fileHandle_p = NULL;
50    Finley_ElementTypeId* element_type=NULL;    ElementTypeId* element_type=NULL;
51    
52    
53    Esys_MPIInfo *mpi_info = Esys_MPIInfo_alloc( MPI_COMM_WORLD );    Esys_MPIInfo *mpi_info = Esys_MPIInfo_alloc( MPI_COMM_WORLD );
# Line 126  Finley_Mesh* Finley_Mesh_readGmsh(char* Line 126  Finley_Mesh* Finley_Mesh_readGmsh(char*
126         else if(!strncmp(&line[1], "ELM", 3) ||         else if(!strncmp(&line[1], "ELM", 3) ||
127          !strncmp(&line[1], "Elements", 8)) {          !strncmp(&line[1], "Elements", 8)) {
128        
129           Finley_ElementTypeId final_element_type = Finley_NoRef;           ElementTypeId final_element_type = NoRef;
130           Finley_ElementTypeId final_face_element_type = Finley_NoRef;           ElementTypeId final_face_element_type = NoRef;
131           Finley_ElementTypeId contact_element_type = Finley_NoRef;           ElementTypeId contact_element_type = NoRef;
132           numElements=0;           numElements=0;
133           numFaceElements=0;           numFaceElements=0;
134           scan_ret = fscanf(fileHandle_p, "%d", &totalNumElements);           scan_ret = fscanf(fileHandle_p, "%d", &totalNumElements);
# Line 138  Finley_Mesh* Finley_Mesh_readGmsh(char* Line 138  Finley_Mesh* Finley_Mesh_readGmsh(char*
138           tag=new index_t[totalNumElements];           tag=new index_t[totalNumElements];
139        
140        
141           element_type=new Finley_ElementTypeId[totalNumElements];           element_type=new ElementTypeId[totalNumElements];
142           vertices=new index_t[totalNumElements*MAX_numNodes_gmsh];           vertices=new index_t[totalNumElements*MAX_numNodes_gmsh];
143           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) ) ) {
144              /* read all in */              /* read all in */
# Line 147  Finley_Mesh* Finley_Mesh_readGmsh(char* Line 147  Finley_Mesh* Finley_Mesh_readGmsh(char*
147            FSCANF_CHECK(scan_ret, "fscanf: Finley_Mesh_readGmsh");            FSCANF_CHECK(scan_ret, "fscanf: Finley_Mesh_readGmsh");
148                switch (gmsh_type) {                switch (gmsh_type) {
149                    case 1:  /* line order 1 */                    case 1:  /* line order 1 */
150                        element_type[e]=Finley_Line2;                        element_type[e]=Line2;
151                        element_dim=1;                        element_dim=1;
152                        numNodesPerElement=2;                          numNodesPerElement=2;  
153                        break;                        break;
154                    case 2:  /* triangle order 1 */                    case 2:  /* triangle order 1 */
155                        element_type[e]=Finley_Tri3;                        element_type[e]=Tri3;
156                        numNodesPerElement= 3;                        numNodesPerElement= 3;
157                        element_dim=2;                        element_dim=2;
158                        break;                        break;
159                    case 3:  /* quadrilateral order 1 */                    case 3:  /* quadrilateral order 1 */
160                        element_type[e]=Finley_Rec4;                        element_type[e]=Rec4;
161                        numNodesPerElement= 4;                        numNodesPerElement= 4;
162                        element_dim=2;                        element_dim=2;
163                        break;                        break;
164                    case 4:  /* tetrahedron order 1 */                    case 4:  /* tetrahedron order 1 */
165                        element_type[e]=Finley_Tet4;                        element_type[e]=Tet4;
166                        numNodesPerElement= 4;                        numNodesPerElement= 4;
167                        element_dim=3;                        element_dim=3;
168                        break;                        break;
169                    case 5:  /* hexahedron order 1 */                    case 5:  /* hexahedron order 1 */
170                        element_type[e]=Finley_Hex8;                        element_type[e]=Hex8;
171                        numNodesPerElement= 8;                        numNodesPerElement= 8;
172                        element_dim=3;                        element_dim=3;
173                        break;                        break;
174                    case 8:  /* line order 2 */                    case 8:  /* line order 2 */
175                        if (useMacroElements) {                        if (useMacroElements) {
176                            element_type[e]=Finley_Line3Macro;                            element_type[e]=Line3Macro;
177                        } else {                        } else {
178                            element_type[e]=Finley_Line3;                            element_type[e]=Line3;
179                        }                        }
180                        numNodesPerElement= 3;                        numNodesPerElement= 3;
181                        element_dim=1;                        element_dim=1;
182                        break;                        break;
183                    case 9:  /* triangle order 2 */                    case 9:  /* triangle order 2 */
184                        if (useMacroElements) {                        if (useMacroElements) {
185                             element_type[e]=Finley_Tri6Macro;                             element_type[e]=Tri6Macro;
186                        } else {                        } else {
187                             element_type[e]=Finley_Tri6;                             element_type[e]=Tri6;
188                        }                        }
189                        numNodesPerElement= 6;                        numNodesPerElement= 6;
190                        element_dim=2;                        element_dim=2;
191                        break;                        break;
192                    case 10:  /* quadrilateral order 2 */                    case 10:  /* quadrilateral order 2 */
193                        if (useMacroElements) {                        if (useMacroElements) {
194                            element_type[e]=Finley_Rec9Macro;                            element_type[e]=Rec9Macro;
195                        } else {                        } else {
196                            element_type[e]=Finley_Rec9;                            element_type[e]=Rec9;
197                        }                        }
198                        numNodesPerElement= 9;                        numNodesPerElement= 9;
199                        element_dim=2;                        element_dim=2;
200                        break;                        break;
201                    case 11:  /* tetrahedron order 2 */                    case 11:  /* tetrahedron order 2 */
202                        if (useMacroElements) {                        if (useMacroElements) {
203                            element_type[e]=Finley_Tet10Macro;                            element_type[e]=Tet10Macro;
204                        } else {                        } else {
205                            element_type[e]=Finley_Tet10;                            element_type[e]=Tet10;
206                        }                        }
207                        numNodesPerElement= 10;                        numNodesPerElement= 10;
208                        element_dim=3;                        element_dim=3;
209                        break;                        break;
210                    case 16:  /* rectangular order 2 */                    case 16:  /* rectangular order 2 */
211                        element_type[e]=Finley_Rec8;                        element_type[e]=Rec8;
212                        numNodesPerElement= 8;                        numNodesPerElement= 8;
213                        element_dim=2;                        element_dim=2;
214                        break;                        break;
215                    case 17:  /* hexahedron order 2 */                    case 17:  /* hexahedron order 2 */
216                        element_type[e]=Finley_Hex20;                        element_type[e]=Hex20;
217                        numNodesPerElement= 20;                        numNodesPerElement= 20;
218                        element_dim=3;                        element_dim=3;
219                        break;                        break;
220                    case 15 :  /* point */                    case 15 :  /* point */
221                        element_type[e]=Finley_Point1;                        element_type[e]=Point1;
222                        numNodesPerElement= 1;                        numNodesPerElement= 1;
223                        element_dim=0;                        element_dim=0;
224                        break;                        break;
225                    default:                    default:
226                       element_type[e]=Finley_NoRef;                       element_type[e]=NoRef;
227                       sprintf(error_msg,"Unexpected gmsh element type %d in mesh file %s.",gmsh_type,fname);                       sprintf(error_msg,"Unexpected gmsh element type %d in mesh file %s.",gmsh_type,fname);
228                       Finley_setError(IO_ERROR,error_msg);                       Finley_setError(IO_ERROR,error_msg);
229                }                }
230                if (element_dim == numDim) {                if (element_dim == numDim) {
231                   if (final_element_type == Finley_NoRef) {                   if (final_element_type == NoRef) {
232                      final_element_type = element_type[e];                      final_element_type = element_type[e];
233                   } else if (final_element_type != element_type[e]) {                   } else if (final_element_type != element_type[e]) {
234                       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 237  Finley_Mesh* Finley_Mesh_readGmsh(char* Line 237  Finley_Mesh* Finley_Mesh_readGmsh(char*
237                   }                   }
238                   numElements++;                   numElements++;
239                } else if (element_dim == numDim-1) {                } else if (element_dim == numDim-1) {
240                   if (final_face_element_type == Finley_NoRef) {                   if (final_face_element_type == NoRef) {
241                      final_face_element_type = element_type[e];                      final_face_element_type = element_type[e];
242                   } else if (final_face_element_type != element_type[e]) {                   } else if (final_face_element_type != element_type[e]) {
243                       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 279  Finley_Mesh* Finley_Mesh_readGmsh(char* Line 279  Finley_Mesh* Finley_Mesh_readGmsh(char*
279              FSCANF_CHECK(scan_ret, "fscanf: Finley_Mesh_readGmsh");              FSCANF_CHECK(scan_ret, "fscanf: Finley_Mesh_readGmsh");
280            }            }
281                /* for tet10 the last two nodes need to be swapped */                /* for tet10 the last two nodes need to be swapped */
282                if ((element_type[e]==Finley_Tet10) || (element_type[e]==Finley_Tet10Macro)) {                if ((element_type[e]==Tet10) || (element_type[e]==Tet10Macro)) {
283                     itmp=vertices[INDEX2(9,e,MAX_numNodes_gmsh)];                     itmp=vertices[INDEX2(9,e,MAX_numNodes_gmsh)];
284                     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)];
285                     vertices[INDEX2(8,e,MAX_numNodes_gmsh)]=itmp;                     vertices[INDEX2(8,e,MAX_numNodes_gmsh)]=itmp;
# Line 289  Finley_Mesh* Finley_Mesh_readGmsh(char* Line 289  Finley_Mesh* Finley_Mesh_readGmsh(char*
289                    
290              if (Finley_noError()) {              if (Finley_noError()) {
291                /* first we have to identify the elements to define Elements and FaceElements */                /* first we have to identify the elements to define Elements and FaceElements */
292                if (final_element_type == Finley_NoRef) {                if (final_element_type == NoRef) {
293                   if (numDim==1) {                   if (numDim==1) {
294                      final_element_type=Finley_Line2;                      final_element_type=Line2;
295                   } else if (numDim==2) {                   } else if (numDim==2) {
296                      final_element_type=Finley_Tri3;                      final_element_type=Tri3;
297                   } else if (numDim==3) {                   } else if (numDim==3) {
298                      final_element_type=Finley_Tet4;                      final_element_type=Tet4;
299                   }                   }
300                }                }
301                if (final_face_element_type == Finley_NoRef) {                if (final_face_element_type == NoRef) {
302                   if (numDim==1) {                   if (numDim==1) {
303                      final_face_element_type=Finley_Point1;                      final_face_element_type=Point1;
304                   } else if (numDim==2) {                   } else if (numDim==2) {
305                      final_face_element_type=Finley_Line2;                      final_face_element_type=Line2;
306                   } else if (numDim==3) {                   } else if (numDim==3) {
307                      final_face_element_type=Finley_Tri3;                      final_face_element_type=Tri3;
308                   }                   }
309                }                }
310                if (final_face_element_type == Finley_Line2) {                if (final_face_element_type == Line2) {
311                    contact_element_type=Finley_Line2_Contact;                    contact_element_type=Line2_Contact;
312                } else  if ( (final_face_element_type == Finley_Line3) || (final_face_element_type == Finley_Line3Macro) ) {                } else  if ( (final_face_element_type == Line3) || (final_face_element_type == Line3Macro) ) {
313                    contact_element_type=Finley_Line3_Contact;                    contact_element_type=Line3_Contact;
314                } else  if (final_face_element_type == Finley_Tri3) {                } else  if (final_face_element_type == Tri3) {
315                    contact_element_type=Finley_Tri3_Contact;                    contact_element_type=Tri3_Contact;
316                } else  if ( (final_face_element_type == Finley_Tri6) || (final_face_element_type == Finley_Tri6Macro)) {                } else  if ( (final_face_element_type == Tri6) || (final_face_element_type == Tri6Macro)) {
317                    contact_element_type=Finley_Tri6_Contact;                    contact_element_type=Tri6_Contact;
318                } else {                } else {
319                    contact_element_type=Finley_Point1_Contact;                    contact_element_type=Point1_Contact;
320                }                }
321                refElements= Finley_ReferenceElementSet_alloc(final_element_type,order, reduced_order);                refElements= ReferenceElementSet_alloc(final_element_type,order, reduced_order);
322                refFaceElements=Finley_ReferenceElementSet_alloc(final_face_element_type,order, reduced_order);                refFaceElements=ReferenceElementSet_alloc(final_face_element_type,order, reduced_order);
323                refContactElements= Finley_ReferenceElementSet_alloc(contact_element_type,order, reduced_order);                refContactElements= ReferenceElementSet_alloc(contact_element_type,order, reduced_order);
324                refPoints= Finley_ReferenceElementSet_alloc(Finley_Point1,order, reduced_order);                refPoints= ReferenceElementSet_alloc(Point1,order, reduced_order);
325                mesh_p->Elements=new ElementFile(refElements, mpi_info);                mesh_p->Elements=new ElementFile(refElements, mpi_info);
326                mesh_p->FaceElements=new ElementFile(refFaceElements, mpi_info);                mesh_p->FaceElements=new ElementFile(refFaceElements, mpi_info);
327                mesh_p->ContactElements=new ElementFile(refContactElements, mpi_info);                mesh_p->ContactElements=new ElementFile(refContactElements, mpi_info);
# Line 416  Finley_Mesh* Finley_Mesh_readGmsh(char* Line 416  Finley_Mesh* Finley_Mesh_readGmsh(char*
416       /* rearrange elements: */       /* rearrange elements: */
417       if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);       if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
418       /* free up memory */       /* free up memory */
419       Finley_ReferenceElementSet_dealloc(refPoints);       ReferenceElementSet_dealloc(refPoints);
420       Finley_ReferenceElementSet_dealloc(refContactElements);       ReferenceElementSet_dealloc(refContactElements);
421       Finley_ReferenceElementSet_dealloc(refFaceElements);       ReferenceElementSet_dealloc(refFaceElements);
422       Finley_ReferenceElementSet_dealloc(refElements);       ReferenceElementSet_dealloc(refElements);
423       Esys_MPIInfo_free( mpi_info );       Esys_MPIInfo_free( mpi_info );
424       if (! Finley_noError()) {       if (! Finley_noError()) {
425          Finley_Mesh_free(mesh_p);          Finley_Mesh_free(mesh_p);

Legend:
Removed from v.4491  
changed lines
  Added in v.4492

  ViewVC Help
Powered by ViewVC 1.1.26