/[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

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

Legend:
Removed from v.149  
changed lines
  Added in v.2059

  ViewVC Help
Powered by ViewVC 1.1.26