/[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 2310 by gross, Mon Mar 16 04:07:43 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-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); 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        temp1[0] = numDim;
313        temp1[1] = numNodes;
314        temp1[2] = strlen(name) + 1;
315        error_code = MPI_Bcast (temp1, 3, MPI_INT,  0, mpi_info->comm);
316        if (error_code != MPI_SUCCESS) {
317          Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of temp1 failed");
318          return NULL;
319        }
320        numDim = temp1[0];
321        numNodes = temp1[1];
322        error_code = MPI_Bcast (name, temp1[2], MPI_CHAR, 0, mpi_info->comm);
323        if (error_code != MPI_SUCCESS) {
324          Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of name failed");
325          return NULL;
326        }
327      }
328    #endif
329    
   /*   resolve id's : */  
330    
331    if (Finley_noError()) {       /* allocate mesh */
332       Finley_Mesh_resolveNodeIds(mesh_p);       mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
333    }       if (Finley_noError()) {
334        /* Each CPU will get at most chunkSize nodes so the message has to be sufficiently large */
335        int chunkSize = numNodes / mpi_info->size + 1, totalNodes=0, chunkNodes=0, chunkEle=0, nextCPU=1;
336        int *tempInts = TMPMEMALLOC(chunkSize*3+1, index_t);        /* Stores the integer message data */
337        double *tempCoords = TMPMEMALLOC(chunkSize*numDim, double); /* Stores the double message data */
338    
339        /*
340          Read chunkSize nodes, send it in a chunk to worker CPU which copies chunk into its local mesh_p
341          It doesn't matter that a CPU has the wrong nodes for its elements, this is sorted out later
342          First chunk sent to CPU 1, second to CPU 2, ...
343          Last chunk stays on CPU 0 (the master)
344          The three columns of integers (Id, gDOF, Tag) are gathered into a single array tempInts and sent together in a single MPI message
345        */
346    
347        if (mpi_info->rank == 0) {  /* Master */
348          for (;;) {            /* Infinite loop */
349            for (i0=0; i0<chunkSize*3+1; i0++) tempInts[i0] = -1;
350            for (i0=0; i0<chunkSize*numDim; i0++) tempCoords[i0] = -1.0;
351            chunkNodes = 0;
352            for (i1=0; i1<chunkSize; i1++) {
353              if (totalNodes >= numNodes) break;    /* End of inner loop */
354                  if (1 == numDim) {
355            scan_ret = fscanf(fileHandle_p, "%d %d %d %le\n",
356              &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
357              &tempCoords[i1*numDim+0]);
358            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
359              }
360                  if (2 == numDim) {
361            scan_ret = fscanf(fileHandle_p, "%d %d %d %le %le\n",
362              &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
363              &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1]);
364            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
365              }
366                  if (3 == numDim) {
367            scan_ret = fscanf(fileHandle_p, "%d %d %d %le %le %le\n",
368              &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
369              &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1], &tempCoords[i1*numDim+2]);
370            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
371              }
372              totalNodes++; /* When do we quit the infinite loop? */
373              chunkNodes++; /* How many nodes do we actually have in this chunk? It may be smaller than chunkSize. */
374            }
375            if (chunkNodes > chunkSize) {
376                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: error reading chunks of mesh, data too large for message size");
377                  return NULL;
378            }
379    #ifdef PASO_MPI
380            /* Eventually we'll send chunkSize nodes to each CPU numbered 1 ... mpi_info->size-1, here goes one of them */
381            if (nextCPU < mpi_info->size) {
382                  int mpi_error;
383              tempInts[chunkSize*3] = chunkNodes;   /* The message has one more int to send chunkNodes */
384              mpi_error = MPI_Send(tempInts, chunkSize*3+1, MPI_INT, nextCPU, 81720, mpi_info->comm);
385              if ( mpi_error != MPI_SUCCESS ) {
386                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts failed");
387                    return NULL;
388              }
389              mpi_error = MPI_Send(tempCoords, chunkSize*numDim, MPI_DOUBLE, nextCPU, 81721, mpi_info->comm);
390              if ( mpi_error != MPI_SUCCESS ) {
391                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempCoords failed");
392                    return NULL;
393              }
394            }
395    #endif
396            nextCPU++;
397            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
398            if (nextCPU > mpi_info->size) break; /* End infinite loop */
399          } /* Infinite loop */
400        }   /* End master */
401        else {  /* Worker */
402    #ifdef PASO_MPI
403          /* Each worker receives two messages */
404          MPI_Status status;
405              int mpi_error;
406          mpi_error = MPI_Recv(tempInts, chunkSize*3+1, MPI_INT, 0, 81720, mpi_info->comm, &status);
407          if ( mpi_error != MPI_SUCCESS ) {
408                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts failed");
409                return NULL;
410          }
411          mpi_error = MPI_Recv(tempCoords, chunkSize*numDim, MPI_DOUBLE, 0, 81721, mpi_info->comm, &status);
412          if ( mpi_error != MPI_SUCCESS ) {
413                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempCoords failed");
414                return NULL;
415          }
416          chunkNodes = tempInts[chunkSize*3];   /* How many nodes are in this workers chunk? */
417    #endif
418        }   /* Worker */
419    
420    /* rearrange elements: */      /* Copy node data from tempMem to mesh_p */
421            Finley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes);
422            if (Finley_noError()) {
423          for (i0=0; i0<chunkNodes; i0++) {
424            mesh_p->Nodes->Id[i0]               = tempInts[0+i0];
425            mesh_p->Nodes->globalDegreesOfFreedom[i0]       = tempInts[chunkSize+i0];
426            mesh_p->Nodes->Tag[i0]              = tempInts[chunkSize*2+i0];
427            for (i1=0; i1<numDim; i1++) {
428              mesh_p->Nodes->Coordinates[INDEX2(i1,i0,numDim)]  = tempCoords[i0*numDim+i1];
429            }
430              }
431            }
432    
433        TMPMEMFREE(tempInts);
434        TMPMEMFREE(tempCoords);
435    
436            /* read elements */
437    
438        /* Read the element typeID */
439            if (Finley_noError()) {
440          if (mpi_info->rank == 0) {
441                scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
442            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
443                typeID=Finley_RefElement_getTypeId(element_type);
444          }
445    #ifdef PASO_MPI
446          if (mpi_info->size > 1) {
447            int temp1[2], mpi_error;
448            temp1[0] = (int) typeID;
449            temp1[1] = numEle;
450            mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
451            if (mpi_error != MPI_SUCCESS) {
452              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
453              return NULL;
454            }
455            typeID = (ElementTypeId) temp1[0];
456            numEle = temp1[1];
457          }
458    #endif
459              if (typeID==NoType) {
460                sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
461                Finley_setError(VALUE_ERROR, error_msg);
462              }
463        }
464    
465          /* Allocate the ElementFile */
466          mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
467          numNodes = mesh_p->Elements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
468    
469          /* Read the element data */
470          if (Finley_noError()) {
471        int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
472        int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
473        /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
474        if (mpi_info->rank == 0) {  /* Master */
475          for (;;) {            /* Infinite loop */
476            for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
477            chunkEle = 0;
478            for (i0=0; i0<chunkSize; i0++) {
479              if (totalEle >= numEle) break; /* End inner loop */
480              scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
481              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
482              for (i1 = 0; i1 < numNodes; i1++) {
483            scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
484                FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
485              }
486              scan_ret = fscanf(fileHandle_p, "\n");
487              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
488              totalEle++;
489              chunkEle++;
490            }
491    #ifdef PASO_MPI
492            /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
493            if (nextCPU < mpi_info->size) {
494                  int mpi_error;
495              tempInts[chunkSize*(2+numNodes)] = chunkEle;
496              mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81722, mpi_info->comm);
497              if ( mpi_error != MPI_SUCCESS ) {
498                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Elements failed");
499                    return NULL;
500              }
501            }
502    #endif
503            nextCPU++;
504            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
505            if (nextCPU > mpi_info->size) break; /* End infinite loop */
506          } /* Infinite loop */
507        }   /* End master */
508        else {  /* Worker */
509    #ifdef PASO_MPI
510          /* Each worker receives one message */
511          MPI_Status status;
512          int mpi_error;
513          mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81722, mpi_info->comm, &status);
514          if ( mpi_error != MPI_SUCCESS ) {
515                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for Elements failed");
516                return NULL;
517          }
518          chunkEle = tempInts[chunkSize*(2+numNodes)];
519    #endif
520        }   /* Worker */
521    
522    if (Finley_noError()) {      /* Copy Element data from tempInts to mesh_p */
523       Finley_Mesh_prepare(mesh_p);      Finley_ElementFile_allocTable(mesh_p->Elements, chunkEle);
524    }          mesh_p->Elements->minColor=0;
525            mesh_p->Elements->maxColor=chunkEle-1;
526            if (Finley_noError()) {
527              #pragma omp parallel for private (i0, i1)
528          for (i0=0; i0<chunkEle; i0++) {
529            mesh_p->Elements->Id[i0]    = tempInts[i0*(2+numNodes)+0];
530            mesh_p->Elements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];
531                mesh_p->Elements->Owner[i0]  =mpi_info->rank;
532                mesh_p->Elements->Color[i0] = i0;
533            for (i1 = 0; i1 < numNodes; i1++) {
534              mesh_p->Elements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
535            }
536              }
537            }
538    
539        TMPMEMFREE(tempInts);
540          } /* end of Read the element data */
541    
542            /* read face elements */
543    
544        /* Read the element typeID */
545            if (Finley_noError()) {
546          if (mpi_info->rank == 0) {
547                scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
548            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
549                typeID=Finley_RefElement_getTypeId(element_type);
550          }
551    #ifdef PASO_MPI
552          if (mpi_info->size > 1) {
553            int temp1[2], mpi_error;
554            temp1[0] = (int) typeID;
555            temp1[1] = numEle;
556            mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
557            if (mpi_error != MPI_SUCCESS) {
558              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
559              return NULL;
560            }
561            typeID = (ElementTypeId) temp1[0];
562            numEle = temp1[1];
563          }
564    #endif
565              if (typeID==NoType) {
566                sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
567                Finley_setError(VALUE_ERROR, error_msg);
568              }
569        }
570    
571          /* Allocate the ElementFile */
572          mesh_p->FaceElements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
573          numNodes = mesh_p->FaceElements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
574    
575          /* Read the face element data */
576          if (Finley_noError()) {
577        int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
578        int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
579        /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
580        if (mpi_info->rank == 0) {  /* Master */
581          for (;;) {            /* Infinite loop */
582            for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
583            chunkEle = 0;
584            for (i0=0; i0<chunkSize; i0++) {
585              if (totalEle >= numEle) break; /* End inner loop */
586              scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
587              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
588              for (i1 = 0; i1 < numNodes; i1++) {
589            scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
590                FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
591              }
592              scan_ret = fscanf(fileHandle_p, "\n");
593              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
594              totalEle++;
595              chunkEle++;
596            }
597    #ifdef PASO_MPI
598            /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
599            if (nextCPU < mpi_info->size) {
600                  int mpi_error;
601              tempInts[chunkSize*(2+numNodes)] = chunkEle;
602              mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81723, mpi_info->comm);
603              if ( mpi_error != MPI_SUCCESS ) {
604                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for FaceElements failed");
605                    return NULL;
606              }
607            }
608    #endif
609            nextCPU++;
610            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
611            if (nextCPU > mpi_info->size) break; /* End infinite loop */
612          } /* Infinite loop */
613        }   /* End master */
614        else {  /* Worker */
615    #ifdef PASO_MPI
616          /* Each worker receives one message */
617          MPI_Status status;
618          int mpi_error;
619          mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81723, mpi_info->comm, &status);
620          if ( mpi_error != MPI_SUCCESS ) {
621                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for FaceElements failed");
622                return NULL;
623          }
624          chunkEle = tempInts[chunkSize*(2+numNodes)];
625    #endif
626        }   /* Worker */
627    
628    /* optimize node labeling*/      /* Copy Element data from tempInts to mesh_p */
629        Finley_ElementFile_allocTable(mesh_p->FaceElements, chunkEle);
630            mesh_p->FaceElements->minColor=0;
631            mesh_p->FaceElements->maxColor=chunkEle-1;
632            if (Finley_noError()) {
633              #pragma omp parallel for private (i0, i1)
634          for (i0=0; i0<chunkEle; i0++) {
635            mesh_p->FaceElements->Id[i0]    = tempInts[i0*(2+numNodes)+0];
636            mesh_p->FaceElements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];
637                mesh_p->FaceElements->Owner[i0]  =mpi_info->rank;
638                mesh_p->FaceElements->Color[i0] = i0;
639            for (i1 = 0; i1 < numNodes; i1++) {
640              mesh_p->FaceElements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
641            }
642              }
643            }
644    
645        TMPMEMFREE(tempInts);
646          } /* end of Read the face element data */
647    
648            /* read contact elements */
649    
650        /* Read the element typeID */
651            if (Finley_noError()) {
652          if (mpi_info->rank == 0) {
653                scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
654            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
655                typeID=Finley_RefElement_getTypeId(element_type);
656          }
657    #ifdef PASO_MPI
658          if (mpi_info->size > 1) {
659            int temp1[2], mpi_error;
660            temp1[0] = (int) typeID;
661            temp1[1] = numEle;
662            mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
663            if (mpi_error != MPI_SUCCESS) {
664              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
665              return NULL;
666            }
667            typeID = (ElementTypeId) temp1[0];
668            numEle = temp1[1];
669          }
670    #endif
671              if (typeID==NoType) {
672                sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
673                Finley_setError(VALUE_ERROR, error_msg);
674              }
675        }
676    
677          /* Allocate the ElementFile */
678          mesh_p->ContactElements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
679          numNodes = mesh_p->ContactElements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
680    
681          /* Read the contact element data */
682          if (Finley_noError()) {
683        int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
684        int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
685        /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
686        if (mpi_info->rank == 0) {  /* Master */
687          for (;;) {            /* Infinite loop */
688            for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
689            chunkEle = 0;
690            for (i0=0; i0<chunkSize; i0++) {
691              if (totalEle >= numEle) break; /* End inner loop */
692              scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
693              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
694              for (i1 = 0; i1 < numNodes; i1++) {
695            scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
696                FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
697              }
698              scan_ret = fscanf(fileHandle_p, "\n");
699              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
700              totalEle++;
701              chunkEle++;
702            }
703    #ifdef PASO_MPI
704            /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
705            if (nextCPU < mpi_info->size) {
706                  int mpi_error;
707              tempInts[chunkSize*(2+numNodes)] = chunkEle;
708              mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81724, mpi_info->comm);
709              if ( mpi_error != MPI_SUCCESS ) {
710                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for ContactElements failed");
711                    return NULL;
712              }
713            }
714    #endif
715            nextCPU++;
716            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
717            if (nextCPU > mpi_info->size) break; /* End infinite loop */
718          } /* Infinite loop */
719        }   /* End master */
720        else {  /* Worker */
721    #ifdef PASO_MPI
722          /* Each worker receives one message */
723          MPI_Status status;
724          int mpi_error;
725          mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81724, mpi_info->comm, &status);
726          if ( mpi_error != MPI_SUCCESS ) {
727                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for ContactElements failed");
728                return NULL;
729          }
730          chunkEle = tempInts[chunkSize*(2+numNodes)];
731    #endif
732        }   /* Worker */
733    
734    if (Finley_noError()) {      /* Copy Element data from tempInts to mesh_p */
735        if (optimize_labeling) Finley_Mesh_optimizeNodeLabeling(mesh_p);      Finley_ElementFile_allocTable(mesh_p->ContactElements, chunkEle);
736    }          mesh_p->ContactElements->minColor=0;
737            mesh_p->ContactElements->maxColor=chunkEle-1;
738            if (Finley_noError()) {
739              #pragma omp parallel for private (i0, i1)
740          for (i0=0; i0<chunkEle; i0++) {
741            mesh_p->ContactElements->Id[i0] = tempInts[i0*(2+numNodes)+0];
742            mesh_p->ContactElements->Tag[i0]    = tempInts[i0*(2+numNodes)+1];
743                mesh_p->ContactElements->Owner[i0]  =mpi_info->rank;
744                mesh_p->ContactElements->Color[i0] = i0;
745            for (i1 = 0; i1 < numNodes; i1++) {
746              mesh_p->ContactElements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
747            }
748              }
749            }
750    
751        TMPMEMFREE(tempInts);
752          } /* end of Read the contact element data */
753    
754            /* read nodal elements */
755    
756        /* Read the element typeID */
757            if (Finley_noError()) {
758          if (mpi_info->rank == 0) {
759                scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
760            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
761                typeID=Finley_RefElement_getTypeId(element_type);
762          }
763    #ifdef PASO_MPI
764          if (mpi_info->size > 1) {
765            int temp1[2], mpi_error;
766            temp1[0] = (int) typeID;
767            temp1[1] = numEle;
768            mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
769            if (mpi_error != MPI_SUCCESS) {
770              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
771              return NULL;
772            }
773            typeID = (ElementTypeId) temp1[0];
774            numEle = temp1[1];
775          }
776    #endif
777              if (typeID==NoType) {
778                sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
779                Finley_setError(VALUE_ERROR, error_msg);
780              }
781        }
782    
783          /* Allocate the ElementFile */
784          mesh_p->Points=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
785          numNodes = mesh_p->Points->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
786    
787          /* Read the nodal element data */
788          if (Finley_noError()) {
789        int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1;
790        int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
791        /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
792        if (mpi_info->rank == 0) {  /* Master */
793          for (;;) {            /* Infinite loop */
794            for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
795            chunkEle = 0;
796            for (i0=0; i0<chunkSize; i0++) {
797              if (totalEle >= numEle) break; /* End inner loop */
798              scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
799              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
800              for (i1 = 0; i1 < numNodes; i1++) {
801            scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
802                FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
803              }
804              scan_ret = fscanf(fileHandle_p, "\n");
805              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
806              totalEle++;
807              chunkEle++;
808            }
809    #ifdef PASO_MPI
810            /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
811            if (nextCPU < mpi_info->size) {
812                  int mpi_error;
813              tempInts[chunkSize*(2+numNodes)] = chunkEle;
814              mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81725, mpi_info->comm);
815              if ( mpi_error != MPI_SUCCESS ) {
816                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Points failed");
817                    return NULL;
818              }
819            }
820    #endif
821            nextCPU++;
822            /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
823            if (nextCPU > mpi_info->size) break; /* End infinite loop */
824          } /* Infinite loop */
825        }   /* End master */
826        else {  /* Worker */
827    #ifdef PASO_MPI
828          /* Each worker receives one message */
829          MPI_Status status;
830          int mpi_error;
831          mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81725, mpi_info->comm, &status);
832          if ( mpi_error != MPI_SUCCESS ) {
833                Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for Points failed");
834                return NULL;
835          }
836          chunkEle = tempInts[chunkSize*(2+numNodes)];
837    #endif
838        }   /* Worker */
839    
840        /* Copy Element data from tempInts to mesh_p */
841        Finley_ElementFile_allocTable(mesh_p->Points, chunkEle);
842            mesh_p->Points->minColor=0;
843            mesh_p->Points->maxColor=chunkEle-1;
844            if (Finley_noError()) {
845              #pragma omp parallel for private (i0, i1)
846          for (i0=0; i0<chunkEle; i0++) {
847            mesh_p->Points->Id[i0]  = tempInts[i0*(2+numNodes)+0];
848            mesh_p->Points->Tag[i0] = tempInts[i0*(2+numNodes)+1];
849                mesh_p->Points->Owner[i0]  =mpi_info->rank;
850                mesh_p->Points->Color[i0] = i0;
851            for (i1 = 0; i1 < numNodes; i1++) {
852              mesh_p->Points->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
853            }
854              }
855            }
856    
857        TMPMEMFREE(tempInts);
858          } /* end of Read the nodal element data */
859    
860          /* get the name tags */
861          if (Finley_noError()) {
862            char *remainder=0, *ptr;
863            size_t len=0;
864    #ifdef PASO_MPI
865            int len_i;
866    #endif
867            int tag_key;
868            if (mpi_info->rank == 0) {  /* Master */
869          /* Read the word 'Tag' */
870          if (! feof(fileHandle_p)) {
871            scan_ret = fscanf(fileHandle_p, "%s\n", name);
872            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
873          }
874    
875    #if defined(_WIN32)  /* windows ftell lies on unix formatted text files */
876    
877          remainder = NULL;
878              len=0;
879          while (1)
880              {
881                 size_t malloc_chunk = 1024;
882                 size_t buff_size = 0;
883                 int ch;
884    
885                 ch = fgetc(fileHandle_p);
886                 if( ch == '\r' )
887                 {
888                    continue;
889                 }
890                 if( len+1 > buff_size )
891                 {
892                    TMPMEMREALLOC(remainder,remainder,buff_size+malloc_chunk,char);
893                 }
894                 if( ch == EOF )
895                 {
896                    /* hit EOF */
897                    remainder[len] = (char)0;
898                    break;
899                 }
900                 remainder[len] = (char)ch;
901                 len++;
902              }
903    #else
904          /* Read rest of file in one chunk, after using seek to find length */
905              {
906                 long cur_pos, end_pos;
907    
908                 cur_pos = ftell(fileHandle_p);
909                 fseek(fileHandle_p, 0L, SEEK_END);
910                 end_pos = ftell(fileHandle_p);
911                 fseek(fileHandle_p, (long)cur_pos, SEEK_SET);
912                 remainder = TMPMEMALLOC(end_pos-cur_pos+1, char);
913                 if (! feof(fileHandle_p))
914                 {
915                    scan_ret = fread(remainder, (size_t) end_pos-cur_pos,
916                                     sizeof(char), fileHandle_p);
917    
918                    FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
919                    remainder[end_pos-cur_pos] = 0;
920                }
921              }
922    #endif
923          len = strlen(remainder);
924          while ((--len)>0 && isspace(remainder[len])) remainder[len]=0;
925          len = strlen(remainder);
926              TMPMEMREALLOC(remainder,remainder,len+1,char);
927            }
928    #ifdef PASO_MPI
929    
930            len_i=(int) len;
931            if (MPI_Bcast (&len_i, 1, MPI_INT,  0, mpi_info->comm) != MPI_SUCCESS)
932            {
933               Finley_setError(PASO_MPI_ERROR,
934                               "Finley_Mesh_read: broadcast of tag len failed");
935               return NULL;
936            }
937            len=(size_t) len_i;
938        if (mpi_info->rank != 0) {
939          remainder = TMPMEMALLOC(len+1, char);
940          remainder[0] = 0;
941        }
942            if (MPI_Bcast (remainder, len+1, MPI_CHAR,  0, mpi_info->comm) !=
943                MPI_SUCCESS)
944            {
945               Finley_setError(PASO_MPI_ERROR,
946                               "Finley_Mesh_read: broadcast of tags failed");
947               return NULL;
948            }
949    #endif
950    
951        if (remainder[0]) {
952              ptr = remainder;
953              do {
954                sscanf(ptr, "%s %d\n", name, &tag_key);
955                if (*name) Finley_Mesh_addTagMap(mesh_p,name,tag_key);
956                ptr++;
957              } while(NULL != (ptr = strchr(ptr, '\n')) && *ptr);
958              TMPMEMFREE(remainder);
959        }
960          }
961    
962         }
963    
964         /* close file */
965         if (mpi_info->rank == 0) fclose(fileHandle_p);
966    
967         /*   resolve id's : */
968         /* rearrange elements: */
969    
970         if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
971         if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
972    
973         /* that's it */
974         #ifdef Finley_TRACE
975         printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
976         #endif
977    
978    /* that's it */    /* that's it */
979    #ifdef Finley_TRACE    if (! Finley_noError()) {
980    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);  
981    }    }
982      Paso_MPIInfo_free( mpi_info );
983    return mesh_p;    return mesh_p;
984  }  }

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

  ViewVC Help
Powered by ViewVC 1.1.26