/[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 1156 by gross, Mon May 21 06:45:14 2007 UTC revision 1887 by ksteube, Wed Oct 15 03:26:25 2008 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  /**************************************************************/  /*******************************************************
3    *
4    * Copyright (c) 2003-2008 by University of Queensland
5    * Earth Systems Science Computational Center (ESSCC)
6    * http://www.uq.edu.au/esscc
7    *
8    * Primary Business: Queensland, Australia
9    * Licensed under the Open Software License version 3.0
10    * http://www.opensource.org/licenses/osl-3.0.php
11    *
12    *******************************************************/
13    
 /*   Finley: read mesh */  
14    
15  /**************************************************************/  /**************************************************************/
16    
17  /*   Author: gross@access.edu.au */  /*   Finley: read mesh */
 /*   Version: $Id$ */  
18    
19  /**************************************************************/  /**************************************************************/
20    
21    #include <ctype.h>
22  #include "Mesh.h"  #include "Mesh.h"
23    
24    #define FSCANF_CHECK(scan_ret, reason) { if (scan_ret == EOF) perror(reason); return NULL; }
25    
26  /**************************************************************/  /**************************************************************/
27    
28  /*  reads a mesh from a Finley file of name fname */  /*  reads a mesh from a Finley file of name fname */
29    
30  Finley_Mesh* Finley_Mesh_read(char* fname,index_t order, index_t reduced_order,  bool_t optimize_labeling) {  Finley_Mesh* Finley_Mesh_read(char* fname,index_t order, index_t reduced_order,  bool_t optimize)
31    
32    {
33    
34      Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
35    dim_t numNodes, numDim, numEle, i0, i1;    dim_t numNodes, numDim, numEle, i0, i1;
36    index_t tag_key;    index_t tag_key;
37    Finley_Mesh *mesh_p=NULL;    Finley_Mesh *mesh_p=NULL;
# Line 37  Finley_Mesh* Finley_Mesh_read(char* fnam Line 40  Finley_Mesh* Finley_Mesh_read(char* fnam
40    double time0=Finley_timer();    double time0=Finley_timer();
41    FILE *fileHandle_p = NULL;    FILE *fileHandle_p = NULL;
42    ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;    ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;
43      int scan_ret;
44      
45    Finley_resetError();    Finley_resetError();
 #ifdef PASO_MPI  
   /* TODO */  
   Finley_setError(SYSTEM_ERROR,"Finley_Mesh_read: MPI is not suporrted yet.");`  
   return NULL;  
 #endif  
46    
47    /* get file handle */    if (mpi_info->size > 1) {
48    fileHandle_p = fopen(fname, "r");       Finley_setError(SYSTEM_ERROR,"Finley_Mesh_read: MPI is not suporrted yet.");
49    if (fileHandle_p==NULL) {    } else {
50      sprintf(error_msg,"Finley_Mesh_read: Opening file %s for reading failed.",fname);       /* get file handle */
51      Finley_setError(IO_ERROR,error_msg);       fileHandle_p = fopen(fname, "r");
52      return NULL;       if (fileHandle_p==NULL) {
53           sprintf(error_msg,"Finley_Mesh_read: Opening file %s for reading failed.",fname);
54           Finley_setError(IO_ERROR,error_msg);
55           Paso_MPIInfo_free( mpi_info );
56           return NULL;
57         }
58      
59         /* read header */
60         sprintf(frm,"%%%d[^\n]",LenString_MAX-1);
61         scan_ret = fscanf(fileHandle_p, frm, name);
62         FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
63      
64         /* get the nodes */
65      
66         scan_ret = fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);
67         FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
68         /* allocate mesh */
69         mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
70         if (Finley_noError()) {
71      
72            /* read nodes */
73            Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
74            if (Finley_noError()) {
75               if (1 == numDim) {
76                   for (i0 = 0; i0 < numNodes; i0++) {
77                    scan_ret = fscanf(fileHandle_p, "%d %d %d %le\n", &mesh_p->Nodes->Id[i0],
78                           &mesh_p->Nodes->globalDegreesOfFreedom[i0], &mesh_p->Nodes->Tag[i0],
79                           &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)]);
80                FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
81               }
82               } else if (2 == numDim) {
83                        for (i0 = 0; i0 < numNodes; i0++) {
84                              scan_ret = fscanf(fileHandle_p, "%d %d %d %le %le\n", &mesh_p->Nodes->Id[i0],
85                                     &mesh_p->Nodes->globalDegreesOfFreedom[i0], &mesh_p->Nodes->Tag[i0],
86                                     &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],
87                                     &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)]);
88                      FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
89                    }
90               } else if (3 == numDim) {
91                        for (i0 = 0; i0 < numNodes; i0++) {
92                              scan_ret = fscanf(fileHandle_p, "%d %d %d %le %le %le\n", &mesh_p->Nodes->Id[i0],
93                                     &mesh_p->Nodes->globalDegreesOfFreedom[i0], &mesh_p->Nodes->Tag[i0],
94                                     &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],
95                                     &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)],
96                                     &mesh_p->Nodes->Coordinates[INDEX2(2,i0,numDim)]);
97                      FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
98                    }
99               } /* if else else */
100            }
101            /* read elements */
102            if (Finley_noError()) {
103      
104               scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
105           FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
106               typeID=Finley_RefElement_getTypeId(element_type);
107               if (typeID==NoType) {
108                 sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s",element_type);
109                 Finley_setError(VALUE_ERROR,error_msg);
110               } else {
111                 /* read the elements */
112                 mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
113                 if (Finley_noError()) {
114                     Finley_ElementFile_allocTable(mesh_p->Elements, numEle);
115                     mesh_p->Elements->minColor=0;
116                     mesh_p->Elements->maxColor=numEle-1;
117                     if (Finley_noError()) {
118                        for (i0 = 0; i0 < numEle; i0++) {
119                          scan_ret = fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);
120                  FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
121                          mesh_p->Elements->Color[i0]=i0;
122                          mesh_p->Elements->Owner[i0]=0;
123                          for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {
124                               scan_ret = fscanf(fileHandle_p, " %d",
125                                  &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);
126                   FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
127                          } /* for i1 */
128                          scan_ret = fscanf(fileHandle_p, "\n");
129                  FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
130                        } /* for i0 */
131                     }
132                 }
133              }
134            }
135            /* get the face elements */
136            if (Finley_noError()) {
137                 scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
138             FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
139                 faceTypeID=Finley_RefElement_getTypeId(element_type);
140                 if (faceTypeID==NoType) {
141                   sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s for face elements",element_type);
142                   Finley_setError(VALUE_ERROR,error_msg);
143                 } else {
144                    mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
145                    if (Finley_noError()) {
146                       Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);
147                       if (Finley_noError()) {
148                          mesh_p->FaceElements->minColor=0;
149                          mesh_p->FaceElements->maxColor=numEle-1;
150                          for (i0 = 0; i0 < numEle; i0++) {
151                            scan_ret = fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);
152                FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
153                            mesh_p->FaceElements->Color[i0]=i0;
154                            mesh_p->FaceElements->Owner[i0]=0;
155                            for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {
156                                 scan_ret = fscanf(fileHandle_p, " %d",
157                                    &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);
158                     FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
159                            }   /* for i1 */
160                            scan_ret = fscanf(fileHandle_p, "\n");
161                FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
162                          } /* for i0 */
163                       }
164                    }
165                 }
166            }
167            /* get the Contact face element */
168            if (Finley_noError()) {
169                 scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
170             FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
171                 contactTypeID=Finley_RefElement_getTypeId(element_type);
172                 if (contactTypeID==NoType) {
173                   sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);
174                   Finley_setError(VALUE_ERROR,error_msg);
175                 } else {
176                   mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
177                   if (Finley_noError()) {
178                       Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);
179                       if (Finley_noError()) {
180                          mesh_p->ContactElements->minColor=0;
181                          mesh_p->ContactElements->maxColor=numEle-1;
182                          for (i0 = 0; i0 < numEle; i0++) {
183                            scan_ret = fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);
184                FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
185                            mesh_p->ContactElements->Color[i0]=i0;
186                            mesh_p->ContactElements->Owner[i0]=0;
187                            for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {
188                                scan_ret = fscanf(fileHandle_p, " %d",
189                                   &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);
190                    FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
191                            }   /* for i1 */
192                            scan_ret = fscanf(fileHandle_p, "\n");
193                FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
194                          } /* for i0 */
195                      }
196                   }
197                 }
198            }  
199            /* get the nodal element */
200            if (Finley_noError()) {
201                 scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
202             FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
203                 pointTypeID=Finley_RefElement_getTypeId(element_type);
204                 if (pointTypeID==NoType) {
205                   sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type);
206                   Finley_setError(VALUE_ERROR,error_msg);
207                 }
208                 mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
209                 if (Finley_noError()) {
210                    Finley_ElementFile_allocTable(mesh_p->Points, numEle);
211                    if (Finley_noError()) {
212                       mesh_p->Points->minColor=0;
213                       mesh_p->Points->maxColor=numEle-1;
214                       for (i0 = 0; i0 < numEle; i0++) {
215                         scan_ret = fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);
216                 FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
217                         mesh_p->Points->Color[i0]=i0;
218                         mesh_p->Points->Owner[i0]=0;
219                         for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {
220                             scan_ret = fscanf(fileHandle_p, " %d",
221                                &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);
222                 FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
223                         }  /* for i1 */
224                         scan_ret = fscanf(fileHandle_p, "\n");
225                 FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
226                       } /* for i0 */
227                    }
228                 }
229            }
230            /* get the name tags */
231            if (Finley_noError()) {
232               if (feof(fileHandle_p) == 0) {
233                  scan_ret = fscanf(fileHandle_p, "%s\n", name);
234              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
235                  while (feof(fileHandle_p) == 0) {
236                       scan_ret = fscanf(fileHandle_p, "%s %d\n", name, &tag_key);
237               FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
238                       Finley_Mesh_addTagMap(mesh_p,name,tag_key);
239                  }
240               }
241            }
242         }
243         /* close file */
244         fclose(fileHandle_p);
245      
246         /*   resolve id's : */
247         /* rearrange elements: */
248      
249         if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
250         if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
251      
252         /* that's it */
253         #ifdef Finley_TRACE
254         printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
255         #endif
256    }    }
257    
258    /* read header */    /* that's it */
259    sprintf(frm,"%%%d[^\n]",LenString_MAX-1);    if (! Finley_noError()) {
260    fscanf(fileHandle_p, frm, name);         Finley_Mesh_free(mesh_p);
   
   /* get the nodes */  
   
   fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);  
   /* allocate mesh */  
   mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order);  
   if (! Finley_noError()) return NULL;  
   
   Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);  
   if (! Finley_noError()) return NULL;  
   
   if (1 == numDim) {  
       for (i0 = 0; i0 < numNodes; i0++)  
     fscanf(fileHandle_p, "%d %d %d %le\n", &mesh_p->Nodes->Id[i0],  
            &mesh_p->Nodes->degreeOfFreedom[i0], &mesh_p->Nodes->Tag[i0],  
            &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)]);  
   } else if (2 == numDim) {  
       for (i0 = 0; i0 < numNodes; i0++)  
     fscanf(fileHandle_p, "%d %d %d %le %le\n", &mesh_p->Nodes->Id[i0],  
            &mesh_p->Nodes->degreeOfFreedom[i0], &mesh_p->Nodes->Tag[i0],  
            &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],  
            &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)]);  
   } else if (3 == numDim) {  
       for (i0 = 0; i0 < numNodes; i0++)  
     fscanf(fileHandle_p, "%d %d %d %le %le %le\n", &mesh_p->Nodes->Id[i0],  
            &mesh_p->Nodes->degreeOfFreedom[i0], &mesh_p->Nodes->Tag[i0],  
            &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],  
            &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)],  
            &mesh_p->Nodes->Coordinates[INDEX2(2,i0,numDim)]);  
   } /* if else else */  
   
   /* get the element type */  
   
   fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);  
   typeID=Finley_RefElement_getTypeId(element_type);  
   if (typeID==NoType) {  
     sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s",element_type);  
     Finley_setError(VALUE_ERROR,error_msg);  
     return NULL;  
   }  
   /* read the elements */  
   mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order);  
   Finley_ElementFile_allocTable(mesh_p->Elements, numEle);  
   mesh_p->Elements->minColor=0;  
   mesh_p->Elements->maxColor=numEle-1;  
   for (i0 = 0; i0 < numEle; i0++) {  
     fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);  
     mesh_p->Elements->Color[i0]=i0;  
     for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {  
          fscanf(fileHandle_p, " %d",  
             &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);  
     }   /* for i1 */  
     fscanf(fileHandle_p, "\n");  
   } /* for i0 */  
   
   /* get the face elements */  
   fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);  
   faceTypeID=Finley_RefElement_getTypeId(element_type);  
   faceTypeID=Finley_RefElement_getTypeId(element_type);  
   if (faceTypeID==NoType) {  
     sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s for face elements",element_type);  
     Finley_setError(VALUE_ERROR,error_msg);  
     return NULL;  
   }  
   mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order, mesh_p->reduced_order);  
   Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);  
   mesh_p->FaceElements->minColor=0;  
   mesh_p->FaceElements->maxColor=numEle-1;  
   for (i0 = 0; i0 < numEle; i0++) {  
     fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);  
     mesh_p->FaceElements->Color[i0]=i0;  
     for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {  
          fscanf(fileHandle_p, " %d",  
             &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);  
     }   /* for i1 */  
     fscanf(fileHandle_p, "\n");  
   } /* for i0 */  
   
   /* get the Contact face element */  
   fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);  
   contactTypeID=Finley_RefElement_getTypeId(element_type);  
   if (contactTypeID==NoType) {  
     sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);  
     Finley_setError(VALUE_ERROR,error_msg);  
     return NULL;  
   }  
   mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order, mesh_p->reduced_order);  
   Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);  
   mesh_p->ContactElements->minColor=0;  
   mesh_p->ContactElements->maxColor=numEle-1;  
   for (i0 = 0; i0 < numEle; i0++) {  
     fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);  
     mesh_p->ContactElements->Color[i0]=i0;  
     for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {  
         fscanf(fileHandle_p, " %d",  
            &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);  
     }   /* for i1 */  
     fscanf(fileHandle_p, "\n");  
   } /* for i0 */  
   
   /* get the nodal element */  
   fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);  
   pointTypeID=Finley_RefElement_getTypeId(element_type);  
   if (pointTypeID==NoType) {  
     sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type);  
     Finley_setError(VALUE_ERROR,error_msg);  
     return NULL;  
261    }    }
262    mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order, mesh_p->reduced_order);    Paso_MPIInfo_free( mpi_info );
263    Finley_ElementFile_allocTable(mesh_p->Points, numEle);    return mesh_p;
264    mesh_p->Points->minColor=0;  }
265    mesh_p->Points->maxColor=numEle-1;  
266    for (i0 = 0; i0 < numEle; i0++) {  Finley_Mesh* Finley_Mesh_read_MPI(char* fname,index_t order, index_t reduced_order,  bool_t optimize)
267      fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);  
268      mesh_p->Points->Color[i0]=i0;  {
269      for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {  
270          fscanf(fileHandle_p, " %d",    Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
271             &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);    dim_t numNodes, numDim, numEle, i0, i1;
272      }   /* for i1 */    Finley_Mesh *mesh_p=NULL;
273      fscanf(fileHandle_p, "\n");    char name[LenString_MAX],element_type[LenString_MAX],frm[20];
274    } /* for i0 */    char error_msg[LenErrorMsg_MAX];
275    /* get the name tags */    double time0=Finley_timer();
276    if (feof(fileHandle_p) == 0) {    FILE *fileHandle_p = NULL;
277       fscanf(fileHandle_p, "%s\n", name);    ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;
278       while (feof(fileHandle_p) == 0) {    Finley_TagMap* tag_map;
279         fscanf(fileHandle_p, "%s %d\n", name, &tag_key);    index_t tag_key;
280         Finley_Mesh_addTagMap(mesh_p,name,tag_key);    int scan_ret;
281    
282      Finley_resetError();
283    
284      if (mpi_info->rank == 0) {
285         /* get file handle */
286         fileHandle_p = fopen(fname, "r");
287         if (fileHandle_p==NULL) {
288           sprintf(error_msg,"Finley_Mesh_read: Opening file %s for reading failed.",fname);
289           Finley_setError(IO_ERROR,error_msg);
290           Paso_MPIInfo_free( mpi_info );
291           return NULL;
292       }       }
293    
294         /* read header */
295         sprintf(frm,"%%%d[^\n]",LenString_MAX-1);
296         scan_ret = fscanf(fileHandle_p, frm, name);
297         FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
298    
299         /* get the number of nodes */
300         scan_ret = fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);
301         FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
302    }    }
   /* close file */  
303    
304    fclose(fileHandle_p);  #ifdef PASO_MPI
305      /* MPI Broadcast numDim, numNodes, name if there are multiple MPI procs*/
306      if (mpi_info->size > 1) {
307        int temp1[3], error_code;
308        temp1[0] = numDim;
309        temp1[1] = numNodes;
310        temp1[2] = strlen(name) + 1;
311        error_code = MPI_Bcast (temp1, 3, MPI_INT,  0, mpi_info->comm);
312        if (error_code != MPI_SUCCESS) {
313          Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of temp1 failed");
314          return NULL;
315        }
316        numDim = temp1[0];
317        numNodes = temp1[1];
318        error_code = MPI_Bcast (name, temp1[2], MPI_CHAR, 0, mpi_info->comm);
319        if (error_code != MPI_SUCCESS) {
320          Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of name failed");
321          return NULL;
322        }
323      }
324    #endif
325    
326         /* allocate mesh */
327         mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
328         if (Finley_noError()) {
329        /* Each CPU will get at most chunkSize nodes so the message has to be sufficiently large */
330        int chunkSize = numNodes / mpi_info->size + 1, totalNodes=0, chunkNodes=0, chunkEle=0, nextCPU=1;
331        int *tempInts = TMPMEMALLOC(chunkSize*3+1, index_t);        /* Stores the integer message data */
332        double *tempCoords = TMPMEMALLOC(chunkSize*numDim, double); /* Stores the double message data */
333    
334        /*
335          Read chunkSize nodes, send it in a chunk to worker CPU which copies chunk into its local mesh_p
336          It doesn't matter that a CPU has the wrong nodes for its elements, this is sorted out later
337          First chunk sent to CPU 1, second to CPU 2, ...
338          Last chunk stays on CPU 0 (the master)
339          The three columns of integers (Id, gDOF, Tag) are gathered into a single array tempInts and sent together in a single MPI message
340        */
341    
342        if (mpi_info->rank == 0) {  /* Master */
343          for (;;) {            /* Infinite loop */
344            for (i0=0; i0<chunkSize*3+1; i0++) tempInts[i0] = -1;
345            for (i0=0; i0<chunkSize*numDim; i0++) tempCoords[i0] = -1.0;
346            chunkNodes = 0;
347            for (i1=0; i1<chunkSize; i1++) {
348              if (totalNodes >= numNodes) break;    /* End of inner loop */
349                  if (1 == numDim) {
350            scan_ret = fscanf(fileHandle_p, "%d %d %d %le\n",
351              &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
352              &tempCoords[i1*numDim+0]);
353            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
354              }
355                  if (2 == numDim) {
356            scan_ret = fscanf(fileHandle_p, "%d %d %d %le %le\n",
357              &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
358              &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1]);
359            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
360              }
361                  if (3 == numDim) {
362            scan_ret = fscanf(fileHandle_p, "%d %d %d %le %le %le\n",
363              &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
364              &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1], &tempCoords[i1*numDim+2]);
365            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
366              }
367              totalNodes++; /* When do we quit the infinite loop? */
368              chunkNodes++; /* How many nodes do we actually have in this chunk? It may be smaller than chunkSize. */
369            }
370            if (chunkNodes > chunkSize) {
371                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: error reading chunks of mesh, data too large for message size");
372                  return NULL;
373            }
374    #ifdef PASO_MPI
375            /* Eventually we'll send chunkSize nodes to each CPU numbered 1 ... mpi_info->size-1, here goes one of them */
376            if (nextCPU < mpi_info->size) {
377                  int mpi_error;
378              tempInts[chunkSize*3] = chunkNodes;   /* The message has one more int to send chunkNodes */
379              mpi_error = MPI_Send(tempInts, chunkSize*3+1, MPI_INT, nextCPU, 81720, mpi_info->comm);
380              if ( mpi_error != MPI_SUCCESS ) {
381                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts failed");
382                    return NULL;
383              }
384              mpi_error = MPI_Send(tempCoords, chunkSize*numDim, MPI_DOUBLE, nextCPU, 81721, mpi_info->comm);
385              if ( mpi_error != MPI_SUCCESS ) {
386                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempCoords failed");
387                    return NULL;
388              }
389            }
390    #endif
391            nextCPU++;
392            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
393            if (nextCPU > mpi_info->size) break; /* End infinite loop */
394          } /* Infinite loop */
395        }   /* End master */
396        else {  /* Worker */
397    #ifdef PASO_MPI
398          /* Each worker receives two messages */
399          MPI_Status status;
400              int mpi_error;
401          mpi_error = MPI_Recv(tempInts, chunkSize*3+1, MPI_INT, 0, 81720, mpi_info->comm, &status);
402          if ( mpi_error != MPI_SUCCESS ) {
403                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts failed");
404                return NULL;
405          }
406          mpi_error = MPI_Recv(tempCoords, chunkSize*numDim, MPI_DOUBLE, 0, 81721, mpi_info->comm, &status);
407          if ( mpi_error != MPI_SUCCESS ) {
408                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempCoords failed");
409                return NULL;
410          }
411          chunkNodes = tempInts[chunkSize*3];   /* How many nodes are in this workers chunk? */
412    #endif
413        }   /* Worker */
414    
415        /* Copy node data from tempMem to mesh_p */
416            Finley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes);
417            if (Finley_noError()) {
418          for (i0=0; i0<chunkNodes; i0++) {
419            mesh_p->Nodes->Id[i0]               = tempInts[0+i0];
420            mesh_p->Nodes->globalDegreesOfFreedom[i0]       = tempInts[chunkSize+i0];
421            mesh_p->Nodes->Tag[i0]              = tempInts[chunkSize*2+i0];
422            for (i1=0; i1<numDim; i1++) {
423              mesh_p->Nodes->Coordinates[INDEX2(i1,i0,numDim)]  = tempCoords[i0*numDim+i1];
424            }
425              }
426            }
427    
428        TMPMEMFREE(tempInts);
429        TMPMEMFREE(tempCoords);
430    
431            /* read elements */
432    
433        /* Read the element typeID */
434            if (Finley_noError()) {
435          if (mpi_info->rank == 0) {
436                scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
437            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
438                typeID=Finley_RefElement_getTypeId(element_type);
439          }
440    #ifdef PASO_MPI
441          if (mpi_info->size > 1) {
442            int temp1[2], mpi_error;
443            temp1[0] = (int) typeID;
444            temp1[1] = numEle;
445            mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
446            if (mpi_error != MPI_SUCCESS) {
447              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
448              return NULL;
449            }
450            typeID = (ElementTypeId) temp1[0];
451            numEle = temp1[1];
452          }
453    #endif
454              if (typeID==NoType) {
455                sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
456                Finley_setError(VALUE_ERROR, error_msg);
457              }
458        }
459    
460          /* Allocate the ElementFile */
461          mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
462          numNodes = mesh_p->Elements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
463    
464          /* Read the element data */
465          if (Finley_noError()) {
466        int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
467        int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
468        /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
469        if (mpi_info->rank == 0) {  /* Master */
470          for (;;) {            /* Infinite loop */
471            for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
472            chunkEle = 0;
473            for (i0=0; i0<chunkSize; i0++) {
474              if (totalEle >= numEle) break; /* End inner loop */
475              scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
476              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
477              for (i1 = 0; i1 < numNodes; i1++) {
478            scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
479                FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
480              }
481              scan_ret = fscanf(fileHandle_p, "\n");
482              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
483              totalEle++;
484              chunkEle++;
485            }
486    #ifdef PASO_MPI
487            /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
488            if (nextCPU < mpi_info->size) {
489                  int mpi_error;
490              tempInts[chunkSize*(2+numNodes)] = chunkEle;
491              mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81722, mpi_info->comm);
492              if ( mpi_error != MPI_SUCCESS ) {
493                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Elements failed");
494                    return NULL;
495              }
496            }
497    #endif
498            nextCPU++;
499            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
500            if (nextCPU > mpi_info->size) break; /* End infinite loop */
501          } /* Infinite loop */
502        }   /* End master */
503        else {  /* Worker */
504    #ifdef PASO_MPI
505          /* Each worker receives one message */
506          MPI_Status status;
507          int mpi_error;
508          mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81722, mpi_info->comm, &status);
509          if ( mpi_error != MPI_SUCCESS ) {
510                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for Elements failed");
511                return NULL;
512          }
513          chunkEle = tempInts[chunkSize*(2+numNodes)];
514    #endif
515        }   /* Worker */
516    
517    /*   resolve id's : */      /* Copy Element data from tempInts to mesh_p */
518        Finley_ElementFile_allocTable(mesh_p->Elements, chunkEle);
519            mesh_p->Elements->minColor=0;
520            mesh_p->Elements->maxColor=chunkEle-1;
521            if (Finley_noError()) {
522              #pragma omp parallel for private (i0, i1)
523          for (i0=0; i0<chunkEle; i0++) {
524            mesh_p->Elements->Id[i0]    = tempInts[i0*(2+numNodes)+0];
525            mesh_p->Elements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];
526                mesh_p->Elements->Owner[i0]  =mpi_info->rank;
527                mesh_p->Elements->Color[i0] = i0;
528            for (i1 = 0; i1 < numNodes; i1++) {
529              mesh_p->Elements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
530            }
531              }
532            }
533    
534        TMPMEMFREE(tempInts);
535          } /* end of Read the element data */
536    
537            /* read face elements */
538    
539        /* Read the element typeID */
540            if (Finley_noError()) {
541          if (mpi_info->rank == 0) {
542                scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
543            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
544                typeID=Finley_RefElement_getTypeId(element_type);
545          }
546    #ifdef PASO_MPI
547          if (mpi_info->size > 1) {
548            int temp1[2], mpi_error;
549            temp1[0] = (int) typeID;
550            temp1[1] = numEle;
551            mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
552            if (mpi_error != MPI_SUCCESS) {
553              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
554              return NULL;
555            }
556            typeID = (ElementTypeId) temp1[0];
557            numEle = temp1[1];
558          }
559    #endif
560              if (typeID==NoType) {
561                sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
562                Finley_setError(VALUE_ERROR, error_msg);
563              }
564        }
565    
566          /* Allocate the ElementFile */
567          mesh_p->FaceElements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
568          numNodes = mesh_p->FaceElements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
569    
570          /* Read the face element data */
571          if (Finley_noError()) {
572        int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
573        int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
574        /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
575        if (mpi_info->rank == 0) {  /* Master */
576          for (;;) {            /* Infinite loop */
577            for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
578            chunkEle = 0;
579            for (i0=0; i0<chunkSize; i0++) {
580              if (totalEle >= numEle) break; /* End inner loop */
581              scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
582              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
583              for (i1 = 0; i1 < numNodes; i1++) {
584            scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
585                FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
586              }
587              scan_ret = fscanf(fileHandle_p, "\n");
588              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
589              totalEle++;
590              chunkEle++;
591            }
592    #ifdef PASO_MPI
593            /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
594            if (nextCPU < mpi_info->size) {
595                  int mpi_error;
596              tempInts[chunkSize*(2+numNodes)] = chunkEle;
597              mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81723, mpi_info->comm);
598              if ( mpi_error != MPI_SUCCESS ) {
599                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for FaceElements failed");
600                    return NULL;
601              }
602            }
603    #endif
604            nextCPU++;
605            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
606            if (nextCPU > mpi_info->size) break; /* End infinite loop */
607          } /* Infinite loop */
608        }   /* End master */
609        else {  /* Worker */
610    #ifdef PASO_MPI
611          /* Each worker receives one message */
612          MPI_Status status;
613          int mpi_error;
614          mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81723, mpi_info->comm, &status);
615          if ( mpi_error != MPI_SUCCESS ) {
616                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for FaceElements failed");
617                return NULL;
618          }
619          chunkEle = tempInts[chunkSize*(2+numNodes)];
620    #endif
621        }   /* Worker */
622    
623    if (Finley_noError()) {      /* Copy Element data from tempInts to mesh_p */
624       Finley_Mesh_resolveNodeIds(mesh_p);      Finley_ElementFile_allocTable(mesh_p->FaceElements, chunkEle);
625    }          mesh_p->FaceElements->minColor=0;
626            mesh_p->FaceElements->maxColor=chunkEle-1;
627            if (Finley_noError()) {
628              #pragma omp parallel for private (i0, i1)
629          for (i0=0; i0<chunkEle; i0++) {
630            mesh_p->FaceElements->Id[i0]    = tempInts[i0*(2+numNodes)+0];
631            mesh_p->FaceElements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];
632                mesh_p->FaceElements->Owner[i0]  =mpi_info->rank;
633                mesh_p->FaceElements->Color[i0] = i0;
634            for (i1 = 0; i1 < numNodes; i1++) {
635              mesh_p->FaceElements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
636            }
637              }
638            }
639    
640        TMPMEMFREE(tempInts);
641          } /* end of Read the face element data */
642    
643            /* read contact elements */
644    
645        /* Read the element typeID */
646            if (Finley_noError()) {
647          if (mpi_info->rank == 0) {
648                scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
649            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
650                typeID=Finley_RefElement_getTypeId(element_type);
651          }
652    #ifdef PASO_MPI
653          if (mpi_info->size > 1) {
654            int temp1[2], mpi_error;
655            temp1[0] = (int) typeID;
656            temp1[1] = numEle;
657            mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
658            if (mpi_error != MPI_SUCCESS) {
659              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
660              return NULL;
661            }
662            typeID = (ElementTypeId) temp1[0];
663            numEle = temp1[1];
664          }
665    #endif
666              if (typeID==NoType) {
667                sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
668                Finley_setError(VALUE_ERROR, error_msg);
669              }
670        }
671    
672          /* Allocate the ElementFile */
673          mesh_p->ContactElements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
674          numNodes = mesh_p->ContactElements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
675    
676          /* Read the contact element data */
677          if (Finley_noError()) {
678        int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
679        int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
680        /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
681        if (mpi_info->rank == 0) {  /* Master */
682          for (;;) {            /* Infinite loop */
683            for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
684            chunkEle = 0;
685            for (i0=0; i0<chunkSize; i0++) {
686              if (totalEle >= numEle) break; /* End inner loop */
687              scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
688              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
689              for (i1 = 0; i1 < numNodes; i1++) {
690            scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
691                FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
692              }
693              scan_ret = fscanf(fileHandle_p, "\n");
694              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
695              totalEle++;
696              chunkEle++;
697            }
698    #ifdef PASO_MPI
699            /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
700            if (nextCPU < mpi_info->size) {
701                  int mpi_error;
702              tempInts[chunkSize*(2+numNodes)] = chunkEle;
703              mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81724, mpi_info->comm);
704              if ( mpi_error != MPI_SUCCESS ) {
705                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for ContactElements failed");
706                    return NULL;
707              }
708            }
709    #endif
710            nextCPU++;
711            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
712            if (nextCPU > mpi_info->size) break; /* End infinite loop */
713          } /* Infinite loop */
714        }   /* End master */
715        else {  /* Worker */
716    #ifdef PASO_MPI
717          /* Each worker receives one message */
718          MPI_Status status;
719          int mpi_error;
720          mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81724, mpi_info->comm, &status);
721          if ( mpi_error != MPI_SUCCESS ) {
722                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for ContactElements failed");
723                return NULL;
724          }
725          chunkEle = tempInts[chunkSize*(2+numNodes)];
726    #endif
727        }   /* Worker */
728    
729    /* rearrange elements: */      /* Copy Element data from tempInts to mesh_p */
730        Finley_ElementFile_allocTable(mesh_p->ContactElements, chunkEle);
731            mesh_p->ContactElements->minColor=0;
732            mesh_p->ContactElements->maxColor=chunkEle-1;
733            if (Finley_noError()) {
734              #pragma omp parallel for private (i0, i1)
735          for (i0=0; i0<chunkEle; i0++) {
736            mesh_p->ContactElements->Id[i0] = tempInts[i0*(2+numNodes)+0];
737            mesh_p->ContactElements->Tag[i0]    = tempInts[i0*(2+numNodes)+1];
738                mesh_p->ContactElements->Owner[i0]  =mpi_info->rank;
739                mesh_p->ContactElements->Color[i0] = i0;
740            for (i1 = 0; i1 < numNodes; i1++) {
741              mesh_p->ContactElements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
742            }
743              }
744            }
745    
746        TMPMEMFREE(tempInts);
747          } /* end of Read the contact element data */
748    
749            /* read nodal elements */
750    
751        /* Read the element typeID */
752            if (Finley_noError()) {
753          if (mpi_info->rank == 0) {
754                scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
755            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
756                typeID=Finley_RefElement_getTypeId(element_type);
757          }
758    #ifdef PASO_MPI
759          if (mpi_info->size > 1) {
760            int temp1[2], mpi_error;
761            temp1[0] = (int) typeID;
762            temp1[1] = numEle;
763            mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
764            if (mpi_error != MPI_SUCCESS) {
765              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
766              return NULL;
767            }
768            typeID = (ElementTypeId) temp1[0];
769            numEle = temp1[1];
770          }
771    #endif
772              if (typeID==NoType) {
773                sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
774                Finley_setError(VALUE_ERROR, error_msg);
775              }
776        }
777    
778          /* Allocate the ElementFile */
779          mesh_p->Points=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
780          numNodes = mesh_p->Points->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
781    
782          /* Read the nodal element data */
783          if (Finley_noError()) {
784        int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
785        int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
786        /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
787        if (mpi_info->rank == 0) {  /* Master */
788          for (;;) {            /* Infinite loop */
789            for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
790            chunkEle = 0;
791            for (i0=0; i0<chunkSize; i0++) {
792              if (totalEle >= numEle) break; /* End inner loop */
793              scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
794              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
795              for (i1 = 0; i1 < numNodes; i1++) {
796            scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
797                FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
798              }
799              scan_ret = fscanf(fileHandle_p, "\n");
800              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
801              totalEle++;
802              chunkEle++;
803            }
804    #ifdef PASO_MPI
805            /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
806            if (nextCPU < mpi_info->size) {
807                  int mpi_error;
808              tempInts[chunkSize*(2+numNodes)] = chunkEle;
809              mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81725, mpi_info->comm);
810              if ( mpi_error != MPI_SUCCESS ) {
811                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Points failed");
812                    return NULL;
813              }
814            }
815    #endif
816            nextCPU++;
817            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
818            if (nextCPU > mpi_info->size) break; /* End infinite loop */
819          } /* Infinite loop */
820        }   /* End master */
821        else {  /* Worker */
822    #ifdef PASO_MPI
823          /* Each worker receives one message */
824          MPI_Status status;
825          int mpi_error;
826          mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81725, mpi_info->comm, &status);
827          if ( mpi_error != MPI_SUCCESS ) {
828                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for Points failed");
829                return NULL;
830          }
831          chunkEle = tempInts[chunkSize*(2+numNodes)];
832    #endif
833        }   /* Worker */
834    
835    if (Finley_noError()) {      /* Copy Element data from tempInts to mesh_p */
836       Finley_Mesh_prepare(mesh_p);      Finley_ElementFile_allocTable(mesh_p->Points, chunkEle);
837    }          mesh_p->Points->minColor=0;
838            mesh_p->Points->maxColor=chunkEle-1;
839            if (Finley_noError()) {
840              #pragma omp parallel for private (i0, i1)
841          for (i0=0; i0<chunkEle; i0++) {
842            mesh_p->Points->Id[i0]  = tempInts[i0*(2+numNodes)+0];
843            mesh_p->Points->Tag[i0] = tempInts[i0*(2+numNodes)+1];
844                mesh_p->Points->Owner[i0]  =mpi_info->rank;
845                mesh_p->Points->Color[i0] = i0;
846            for (i1 = 0; i1 < numNodes; i1++) {
847              mesh_p->Points->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
848            }
849              }
850            }
851    
852        TMPMEMFREE(tempInts);
853          } /* end of Read the nodal element data */
854    
855          /* get the name tags */
856          if (Finley_noError()) {
857            char *remainder, *ptr;
858            int tag_key, error_code;
859            size_t len;
860            long cur_pos, end_pos;
861            if (mpi_info->rank == 0) {  /* Master */
862          /* Read the word 'Tag' */
863          if (! feof(fileHandle_p)) {
864            scan_ret = fscanf(fileHandle_p, "%s\n", name);
865            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
866          }
867          /* Read rest of file in one chunk, after using seek to find length */
868    
869    #if defined(_WIN32)  /* windows ftell lies on unix formatted text files */
870    
871          remainder = NULL;
872              len=0;
873          while (1)
874              {
875                 size_t MALLOC_CHUNK = 1024;
876                 size_t buff_size = 0;
877                 int ch;
878    
879                 ch = fgetc(fileHandle_p);
880                 if( ch == '\r' )
881                 {
882                    continue;
883                 }
884                 if( len+1 > buff_size )
885                 {
886                    TMPMEMREALLOC(remainder,remainder,buff_size+MALLOC_CHUNK,char);
887                 }
888                 if( ch == EOF )
889                 {
890                    /* hit EOF */
891                    remainder[len] = (char)0;
892                    break;
893                 }
894                 remainder[len] = (char)ch;
895                 len++;
896              }
897    #else
898              cur_pos = ftell(fileHandle_p);
899              fseek(fileHandle_p, 0L, SEEK_END);
900              end_pos = ftell(fileHandle_p);
901              fseek(fileHandle_p, (long)cur_pos, SEEK_SET);
902          remainder = TMPMEMALLOC(end_pos-cur_pos+1, char);
903          if (! feof(fileHandle_p)) {
904            scan_ret = fread(remainder, (size_t) end_pos-cur_pos, sizeof(char), fileHandle_p);
905                FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
906          }
907          remainder[end_pos-cur_pos] = 0;
908    #endif
909          len = strlen(remainder);
910          while ((--len)>0 && isspace(remainder[len])) remainder[len]=0;
911          len = strlen(remainder);
912              TMPMEMREALLOC(remainder,remainder,len+1,char);
913            }
914    #ifdef PASO_MPI
915            error_code = MPI_Bcast (&len, 1, MPI_INT,  0, mpi_info->comm);
916            if (error_code != MPI_SUCCESS) {
917              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of tag len failed");
918              return NULL;
919            }
920        if (mpi_info->rank != 0) {
921          remainder = TMPMEMALLOC(len+1, char);
922          remainder[0] = 0;
923        }
924            error_code = MPI_Bcast (remainder, len+1, MPI_CHAR,  0, mpi_info->comm);
925            if (error_code != MPI_SUCCESS) {
926              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of tags failed");
927              return NULL;
928            }
929    #endif
930        if (remainder[0]) {
931              ptr = remainder;
932              do {
933                sscanf(ptr, "%s %d\n", name, &tag_key);
934                if (*name) Finley_Mesh_addTagMap(mesh_p,name,tag_key);
935                ptr++;
936              } while(NULL != (ptr = strchr(ptr, '\n')) && *ptr);
937              TMPMEMFREE(remainder);
938        }
939          }
940    
941    /* optimize node labeling*/       }
942    
943    if (Finley_noError()) {       /* close file */
944        if (optimize_labeling) Finley_Mesh_optimizeNodeLabeling(mesh_p);       if (mpi_info->rank == 0) fclose(fileHandle_p);
945    }  
946         /*   resolve id's : */
947         /* rearrange elements: */
948    
949         if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
950         if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
951    
952         /* that's it */
953         #ifdef Finley_TRACE
954         printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
955         #endif
956    
957    /* that's it */    /* that's it */
958    #ifdef Finley_TRACE    if (! Finley_noError()) {
959    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);  
960    }    }
961      Paso_MPIInfo_free( mpi_info );
962    return mesh_p;    return mesh_p;
963  }  }

Legend:
Removed from v.1156  
changed lines
  Added in v.1887

  ViewVC Help
Powered by ViewVC 1.1.26