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

  ViewVC Help
Powered by ViewVC 1.1.26