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

Legend:
Removed from v.1170  
changed lines
  Added in v.1907

  ViewVC Help
Powered by ViewVC 1.1.26