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

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

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

revision 1028 by gross, Wed Mar 14 00:15:24 2007 UTC revision 1312 by ksteube, Mon Sep 24 06:18:44 2007 UTC
# Line 1  Line 1 
 /*  
  ************************************************************  
  *          Copyright 2006 by ACcESS MNRF                   *  
  *                                                          *  
  *              http://www.access.edu.au                    *  
  *       Primary Business: Queensland, Australia            *  
  *  Licensed under the Open Software License version 3.0    *  
  *     http://www.opensource.org/licenses/osl-3.0.php       *  
  *                                                          *  
  ************************************************************  
 */  
1    
2  /**************************************************************/  /* $Id$ */
3    
4  /*   Finley: read mesh */  /*******************************************************
5     *
6     *           Copyright 2003-2007 by ACceSS MNRF
7     *       Copyright 2007 by University of Queensland
8     *
9     *                http://esscc.uq.edu.au
10     *        Primary Business: Queensland, Australia
11     *  Licensed under the Open Software License version 3.0
12     *     http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15    
16  /**************************************************************/  /**************************************************************/
17    
18  /*   Author: gross@access.edu.au */  /*   Finley: read mesh */
 /*   Version: $Id$ */  
19    
20  /**************************************************************/  /**************************************************************/
21    
# Line 27  Line 25 
25    
26  /*  reads a mesh from a Finley file of name fname */  /*  reads a mesh from a Finley file of name fname */
27    
28  Finley_Mesh* Finley_Mesh_read(char* fname,index_t order) {  Finley_Mesh* Finley_Mesh_read(char* fname,index_t order, index_t reduced_order,  bool_t optimize)
29    
30    {
31    
32      Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
33    dim_t numNodes, numDim, numEle, i0, i1;    dim_t numNodes, numDim, numEle, i0, i1;
34      index_t tag_key;
35    Finley_Mesh *mesh_p=NULL;    Finley_Mesh *mesh_p=NULL;
36    char name[LenString_MAX],element_type[LenString_MAX],frm[20];    char name[LenString_MAX],element_type[LenString_MAX],frm[20];
37    char error_msg[LenErrorMsg_MAX];    char error_msg[LenErrorMsg_MAX];
38    double time0=Finley_timer();    double time0=Finley_timer();
39    FILE *fileHandle_p = NULL;    FILE *fileHandle_p = NULL;
40    ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;    ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;
41      
42    Finley_resetError();    Finley_resetError();
43    
44    /* get file handle */    if (mpi_info->size > 1) {
45    fileHandle_p = fopen(fname, "r");       Finley_setError(SYSTEM_ERROR,"Finley_Mesh_read: MPI is not suporrted yet.");
46    if (fileHandle_p==NULL) {    } else {
47      sprintf(error_msg,"Finley_Mesh_read: Opening file %s for reading failed.",fname);       /* get file handle */
48      Finley_setError(IO_ERROR,error_msg);       fileHandle_p = fopen(fname, "r");
49      return NULL;       if (fileHandle_p==NULL) {
50    }         sprintf(error_msg,"Finley_Mesh_read: Opening file %s for reading failed.",fname);
51           Finley_setError(IO_ERROR,error_msg);
52    /* read header */         Paso_MPIInfo_free( mpi_info );
53    sprintf(frm,"%%%d[^\n]",LenString_MAX-1);         return NULL;
54    fscanf(fileHandle_p, frm, name);       }
55      
56    /* get the nodes */       /* read header */
57         sprintf(frm,"%%%d[^\n]",LenString_MAX-1);
58    fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);       fscanf(fileHandle_p, frm, name);
59    /* allocate mesh */    
60  #ifndef PASO_MPI       /* get the nodes */
61    mesh_p = Finley_Mesh_alloc(name,numDim,order);    
62    if (! Finley_noError()) return NULL;       fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);
63  #else       /* allocate mesh */
64    /* TODO */       mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
65  #endif       if (Finley_noError()) {
66      
67  #ifndef PASO_MPI          /* read nodes */
68    Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);          Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
69    if (! Finley_noError()) return NULL;          if (Finley_noError()) {
70  #else             if (1 == numDim) {
71    /* TODO */                 for (i0 = 0; i0 < numNodes; i0++)
72  #endif                  fscanf(fileHandle_p, "%d %d %d %le\n", &mesh_p->Nodes->Id[i0],
73                           &mesh_p->Nodes->globalDegreesOfFreedom[i0], &mesh_p->Nodes->Tag[i0],
74    if (1 == numDim) {                         &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)]);
75        for (i0 = 0; i0 < numNodes; i0++)             } else if (2 == numDim) {
76      fscanf(fileHandle_p, "%d %d %d %le\n", &mesh_p->Nodes->Id[i0],                      for (i0 = 0; i0 < numNodes; i0++)
77             &mesh_p->Nodes->degreeOfFreedom[i0], &mesh_p->Nodes->Tag[i0],                            fscanf(fileHandle_p, "%d %d %d %le %le\n", &mesh_p->Nodes->Id[i0],
78             &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)]);                                   &mesh_p->Nodes->globalDegreesOfFreedom[i0], &mesh_p->Nodes->Tag[i0],
79    } else if (2 == numDim) {                                   &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],
80        for (i0 = 0; i0 < numNodes; i0++)                                   &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)]);
81      fscanf(fileHandle_p, "%d %d %d %le %le\n", &mesh_p->Nodes->Id[i0],             } else if (3 == numDim) {
82             &mesh_p->Nodes->degreeOfFreedom[i0], &mesh_p->Nodes->Tag[i0],                      for (i0 = 0; i0 < numNodes; i0++)
83             &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],                            fscanf(fileHandle_p, "%d %d %d %le %le %le\n", &mesh_p->Nodes->Id[i0],
84             &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)]);                                   &mesh_p->Nodes->globalDegreesOfFreedom[i0], &mesh_p->Nodes->Tag[i0],
85    } else if (3 == numDim) {                                   &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],
86        for (i0 = 0; i0 < numNodes; i0++)                                   &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)],
87      fscanf(fileHandle_p, "%d %d %d %le %le %le\n", &mesh_p->Nodes->Id[i0],                                   &mesh_p->Nodes->Coordinates[INDEX2(2,i0,numDim)]);
88             &mesh_p->Nodes->degreeOfFreedom[i0], &mesh_p->Nodes->Tag[i0],             } /* if else else */
89             &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],          }
90             &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)],          /* read elements */
91             &mesh_p->Nodes->Coordinates[INDEX2(2,i0,numDim)]);          if (Finley_noError()) {
92    } /* if else else */    
93               fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
94    /* get the element type */             typeID=Finley_RefElement_getTypeId(element_type);
95               if (typeID==NoType) {
96    fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);               sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s",element_type);
97    typeID=Finley_RefElement_getTypeId(element_type);               Finley_setError(VALUE_ERROR,error_msg);
98    if (typeID==NoType) {             } else {
99      sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s",element_type);               /* read the elements */
100      Finley_setError(VALUE_ERROR,error_msg);               mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
101      return NULL;               if (Finley_noError()) {
102    }                   Finley_ElementFile_allocTable(mesh_p->Elements, numEle);
103    /* read the elements */                   mesh_p->Elements->minColor=0;
104  #ifndef PASO_MPI                   mesh_p->Elements->maxColor=numEle-1;
105    mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order);                   if (Finley_noError()) {
106  #else                      for (i0 = 0; i0 < numEle; i0++) {
107    /* TODO */                        fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);
108  #endif                        mesh_p->Elements->Color[i0]=i0;
109    Finley_ElementFile_allocTable(mesh_p->Elements, numEle);                        for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {
110    mesh_p->Elements->minColor=0;                             fscanf(fileHandle_p, " %d",
111    mesh_p->Elements->maxColor=numEle-1;                                &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);
112    for (i0 = 0; i0 < numEle; i0++) {                        } /* for i1 */
113      fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);                        fscanf(fileHandle_p, "\n");
114      mesh_p->Elements->Color[i0]=i0;                      } /* for i0 */
115      for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {                   }
116           fscanf(fileHandle_p, " %d",               }
117              &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);            }
118      }   /* for i1 */          }
119      fscanf(fileHandle_p, "\n");          /* get the face elements */
120    } /* for i0 */          if (Finley_noError()) {
121                 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
122    /* get the face elements */               faceTypeID=Finley_RefElement_getTypeId(element_type);
123    fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);               if (faceTypeID==NoType) {
124    faceTypeID=Finley_RefElement_getTypeId(element_type);                 sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s for face elements",element_type);
125    faceTypeID=Finley_RefElement_getTypeId(element_type);                 Finley_setError(VALUE_ERROR,error_msg);
126    if (faceTypeID==NoType) {               } else {
127      sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s for face elements",element_type);                  mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
128      Finley_setError(VALUE_ERROR,error_msg);                  if (Finley_noError()) {
129      return NULL;                     Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);
130    }                     if (Finley_noError()) {
131  #ifndef PASO_MPI                        mesh_p->FaceElements->minColor=0;
132    mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order);                        mesh_p->FaceElements->maxColor=numEle-1;
133  #else                        for (i0 = 0; i0 < numEle; i0++) {
134    /* TODO */                          fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);
135  #endif                          mesh_p->FaceElements->Color[i0]=i0;
136    Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);                          for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {
137    mesh_p->FaceElements->minColor=0;                               fscanf(fileHandle_p, " %d",
138    mesh_p->FaceElements->maxColor=numEle-1;                                  &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);
139    for (i0 = 0; i0 < numEle; i0++) {                          }   /* for i1 */
140      fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);                          fscanf(fileHandle_p, "\n");
141      mesh_p->FaceElements->Color[i0]=i0;                        } /* for i0 */
142      for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {                     }
143           fscanf(fileHandle_p, " %d",                  }
144              &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);               }
145      }   /* for i1 */          }
146      fscanf(fileHandle_p, "\n");          /* get the Contact face element */
147    } /* for i0 */          if (Finley_noError()) {
148                 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
149    /* get the Contact face element */               contactTypeID=Finley_RefElement_getTypeId(element_type);
150    fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);               if (contactTypeID==NoType) {
151    contactTypeID=Finley_RefElement_getTypeId(element_type);                 sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);
152    if (contactTypeID==NoType) {                 Finley_setError(VALUE_ERROR,error_msg);
153      sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);               } else {
154      Finley_setError(VALUE_ERROR,error_msg);                 mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
155      return NULL;                 if (Finley_noError()) {
156    }                     Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);
157  #ifndef PASO_MPI                     if (Finley_noError()) {
158    mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order);                        mesh_p->ContactElements->minColor=0;
159  #else                        mesh_p->ContactElements->maxColor=numEle-1;
160    /* TODO */                        for (i0 = 0; i0 < numEle; i0++) {
161  #endif                          fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);
162    Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);                          mesh_p->ContactElements->Color[i0]=i0;
163    mesh_p->ContactElements->minColor=0;                          for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {
164    mesh_p->ContactElements->maxColor=numEle-1;                              fscanf(fileHandle_p, " %d",
165    for (i0 = 0; i0 < numEle; i0++) {                                 &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);
166      fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);                          }   /* for i1 */
167      mesh_p->ContactElements->Color[i0]=i0;                          fscanf(fileHandle_p, "\n");
168      for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {                        } /* for i0 */
169          fscanf(fileHandle_p, " %d",                    }
170             &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);                 }
171      }   /* for i1 */               }
172      fscanf(fileHandle_p, "\n");          }  
173    } /* for i0 */          /* get the nodal element */
174            if (Finley_noError()) {
175    /* get the nodal element */               fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
176    fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);               pointTypeID=Finley_RefElement_getTypeId(element_type);
177    pointTypeID=Finley_RefElement_getTypeId(element_type);               if (pointTypeID==NoType) {
178    if (pointTypeID==NoType) {                 sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type);
179      sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type);                 Finley_setError(VALUE_ERROR,error_msg);
180      Finley_setError(VALUE_ERROR,error_msg);               }
181      return NULL;               mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
182                 if (Finley_noError()) {
183                    Finley_ElementFile_allocTable(mesh_p->Points, numEle);
184                    if (Finley_noError()) {
185                       mesh_p->Points->minColor=0;
186                       mesh_p->Points->maxColor=numEle-1;
187                       for (i0 = 0; i0 < numEle; i0++) {
188                         fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);
189                         mesh_p->Points->Color[i0]=i0;
190                         for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {
191                             fscanf(fileHandle_p, " %d",
192                                &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);
193                         }  /* for i1 */
194                         fscanf(fileHandle_p, "\n");
195                       } /* for i0 */
196                    }
197                 }
198            }
199            /* get the name tags */
200            if (Finley_noError()) {
201               if (feof(fileHandle_p) == 0) {
202                  fscanf(fileHandle_p, "%s\n", name);
203                  while (feof(fileHandle_p) == 0) {
204                       fscanf(fileHandle_p, "%s %d\n", name, &tag_key);
205                       Finley_Mesh_addTagMap(mesh_p,name,tag_key);
206                  }
207               }
208            }
209         }
210         /* close file */
211         fclose(fileHandle_p);
212      
213         /*   resolve id's : */
214         /* rearrange elements: */
215      
216         if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
217         if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
218      
219         /* that's it */
220         #ifdef Finley_TRACE
221         printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
222         #endif
223    }    }
 #ifndef PASO_MPI  
   mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order);  
 #else  
   /* TODO */  
 #endif  
   Finley_ElementFile_allocTable(mesh_p->Points, numEle);  
   mesh_p->Points->minColor=0;  
   mesh_p->Points->maxColor=numEle-1;  
   for (i0 = 0; i0 < numEle; i0++) {  
     fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);  
     mesh_p->Points->Color[i0]=i0;  
     for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {  
         fscanf(fileHandle_p, " %d",  
            &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);  
     }   /* for i1 */  
     fscanf(fileHandle_p, "\n");  
   } /* for i0 */  
   
   
   /* close file */  
   
   fclose(fileHandle_p);  
   
   /*   resolve id's : */  
   
   Finley_Mesh_resolveNodeIds(mesh_p);  
   
   /* rearrange elements: */  
   
   Finley_Mesh_prepare(mesh_p);  
224    
225    /* that's it */    /* that's it */
226    #ifdef Finley_TRACE    if (! Finley_noError()) {
227    printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);         Finley_Mesh_free(mesh_p);
   #endif  
   if (Finley_noError()) {  
        if ( ! Finley_Mesh_isPrepared(mesh_p)) {  
           Finley_setError(SYSTEM_ERROR,"Mesh is not prepared for calculation. Contact the programmers.");  
        }  
   } else {  
        Finley_Mesh_dealloc(mesh_p);  
228    }    }
229      Paso_MPIInfo_free( mpi_info );
230    return mesh_p;    return mesh_p;
231  }  }

Legend:
Removed from v.1028  
changed lines
  Added in v.1312

  ViewVC Help
Powered by ViewVC 1.1.26