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

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

  ViewVC Help
Powered by ViewVC 1.1.26