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

temp/finley/src/Mesh_read.c revision 1387 by trankine, Fri Jan 11 07:45:26 2008 UTC trunk/finley/src/Mesh_read.c revision 2754 by gross, Wed Nov 18 07:44:26 2009 UTC
# Line 1  Line 1 
1    
 /* $Id$ */  
   
2  /*******************************************************  /*******************************************************
3   *  *
4   *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2009 by University of Queensland
5   *       Copyright 2007 by University of Queensland  * Earth Systems Science Computational Center (ESSCC)
6   *  * http://www.uq.edu.au/esscc
7   *                http://esscc.uq.edu.au  *
8   *        Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
9   *  Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
10   *     http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
11   *  *
12   *******************************************************/  *******************************************************/
13    
14    
15  /**************************************************************/  /**************************************************************/
16    
# Line 19  Line 18 
18    
19  /**************************************************************/  /**************************************************************/
20    
21    #include <ctype.h>
22  #include "Mesh.h"  #include "Mesh.h"
23    
24  /**************************************************************/  #define FSCANF_CHECK(scan_ret, reason) { if (scan_ret == EOF) { perror(reason); Finley_setError(IO_ERROR,"scan error while reading finley file"); return NULL;} }
   
 /*  reads a mesh from a Finley file of name fname */  
25    
 Finley_Mesh* Finley_Mesh_read(char* fname,index_t order, index_t reduced_order,  bool_t optimize)  
26    
27    Finley_Mesh* Finley_Mesh_read(char* fname,index_t order, index_t reduced_order,  bool_t optimize)
28  {  {
29        Paso_MPIInfo *mpi_info = NULL;
30    Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );      dim_t numNodes, numDim, numEle, i0, i1;
31    dim_t numNodes, numDim, numEle, i0, i1;      Finley_Mesh *mesh_p=NULL;
32    index_t tag_key;      Finley_ReferenceElementSet *refPoints=NULL, *refContactElements=NULL, *refFaceElements=NULL, *refElements=NULL;
33    Finley_Mesh *mesh_p=NULL;      char name[LenString_MAX],element_type[LenString_MAX],frm[20];
34    char name[LenString_MAX],element_type[LenString_MAX],frm[20];      char error_msg[LenErrorMsg_MAX];
35    char error_msg[LenErrorMsg_MAX];      FILE *fileHandle_p = NULL;
36    double time0=Finley_timer();      ElementTypeId typeID=NoRef;
37    FILE *fileHandle_p = NULL;      int scan_ret;
38    ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;  
39          Finley_resetError();
40    Finley_resetError();      mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
41    
42    if (mpi_info->size > 1) {      if (mpi_info->rank == 0) {
43       Finley_setError(SYSTEM_ERROR,"Finley_Mesh_read: MPI is not suporrted yet.");          /* get file handle */
44    } else {          fileHandle_p = fopen(fname, "r");
45       /* get file handle */          if (fileHandle_p==NULL) {
46       fileHandle_p = fopen(fname, "r");              sprintf(error_msg,"Finley_Mesh_read: Opening file %s for reading failed.",fname);
47       if (fileHandle_p==NULL) {              Finley_setError(IO_ERROR,error_msg);
48         sprintf(error_msg,"Finley_Mesh_read: Opening file %s for reading failed.",fname);              Paso_MPIInfo_free( mpi_info );
49         Finley_setError(IO_ERROR,error_msg);              return NULL;
        Paso_MPIInfo_free( mpi_info );  
        return NULL;  
      }  
     
      /* read header */  
      sprintf(frm,"%%%d[^\n]",LenString_MAX-1);  
      fscanf(fileHandle_p, frm, name);  
     
      /* get the nodes */  
     
      fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);  
      /* allocate mesh */  
      mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);  
      if (Finley_noError()) {  
     
         /* read nodes */  
         Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);  
         if (Finley_noError()) {  
            if (1 == numDim) {  
                for (i0 = 0; i0 < numNodes; i0++)  
                 fscanf(fileHandle_p, "%d %d %d %le\n", &mesh_p->Nodes->Id[i0],  
                        &mesh_p->Nodes->globalDegreesOfFreedom[i0], &mesh_p->Nodes->Tag[i0],  
                        &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)]);  
            } else if (2 == numDim) {  
                     for (i0 = 0; i0 < numNodes; i0++)  
                           fscanf(fileHandle_p, "%d %d %d %le %le\n", &mesh_p->Nodes->Id[i0],  
                                  &mesh_p->Nodes->globalDegreesOfFreedom[i0], &mesh_p->Nodes->Tag[i0],  
                                  &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],  
                                  &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)]);  
            } else if (3 == numDim) {  
                     for (i0 = 0; i0 < numNodes; i0++)  
                           fscanf(fileHandle_p, "%d %d %d %le %le %le\n", &mesh_p->Nodes->Id[i0],  
                                  &mesh_p->Nodes->globalDegreesOfFreedom[i0], &mesh_p->Nodes->Tag[i0],  
                                  &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],  
                                  &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)],  
                                  &mesh_p->Nodes->Coordinates[INDEX2(2,i0,numDim)]);  
            } /* if else else */  
50          }          }
51          /* read elements */  
52          if (Finley_noError()) {          /* read header */
53              sprintf(frm,"%%%d[^\n]",LenString_MAX-1);
54             fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);          scan_ret = fscanf(fileHandle_p, frm, name);
55             typeID=Finley_RefElement_getTypeId(element_type);          FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
56             if (typeID==NoType) {  
57               sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s",element_type);          /* get the number of nodes */
58               Finley_setError(VALUE_ERROR,error_msg);          scan_ret = fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);
59             } else {          FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
60               /* read the elements */      }
61               mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);  
62               if (Finley_noError()) {      #ifdef PASO_MPI
63                   Finley_ElementFile_allocTable(mesh_p->Elements, numEle);          /* MPI Broadcast numDim, numNodes, name if there are multiple MPI procs*/
64                   mesh_p->Elements->minColor=0;          if (mpi_info->size > 1) {
65                   mesh_p->Elements->maxColor=numEle-1;              int temp1[3];
66                   if (Finley_noError()) {              if (mpi_info->rank == 0) {
67                      for (i0 = 0; i0 < numEle; i0++) {                  temp1[0] = numDim;
68                        fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);                  temp1[1] = numNodes;
69                        mesh_p->Elements->Color[i0]=i0;                  temp1[2] = strlen(name) + 1;
70                        for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {              } else {
71                             fscanf(fileHandle_p, " %d",                  temp1[0] = 0;
72                                &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);                  temp1[1] = 0;
73                        } /* for i1 */                  temp1[2] = 1;
74                        fscanf(fileHandle_p, "\n");              }
75                      } /* for i0 */              MPI_Bcast (temp1, 3, MPI_INT,  0, mpi_info->comm);
76                   }              numDim = temp1[0];
77               }              numNodes = temp1[1];
78            }              error_code = MPI_Bcast (name, temp1[2], MPI_CHAR, 0, mpi_info->comm);
79          }          }
80          /* get the face elements */      #endif
81          if (Finley_noError()) {  
82               fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);      /* allocate mesh */
83               faceTypeID=Finley_RefElement_getTypeId(element_type);      mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
84               if (faceTypeID==NoType) {      
85                 sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s for face elements",element_type);      if (Finley_noError()) {
86                 Finley_setError(VALUE_ERROR,error_msg);          /* Each CPU will get at most chunkSize nodes so the message has to be sufficiently large */
87               } else {          int chunkSize = numNodes / mpi_info->size + 1, totalNodes=0, chunkNodes=0,  nextCPU=1;
88                  mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);          int *tempInts = TMPMEMALLOC(chunkSize*3+1, index_t);        /* Stores the integer message data */
89                  if (Finley_noError()) {          double *tempCoords = TMPMEMALLOC(chunkSize*numDim, double); /* Stores the double message data */
90                     Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);  
91                     if (Finley_noError()) {          /*
92                        mesh_p->FaceElements->minColor=0;          Read chunkSize nodes, send it in a chunk to worker CPU which copies chunk into its local mesh_p
93                        mesh_p->FaceElements->maxColor=numEle-1;          It doesn't matter that a CPU has the wrong nodes for its elements, this is sorted out later
94                        for (i0 = 0; i0 < numEle; i0++) {          First chunk sent to CPU 1, second to CPU 2, ...
95                          fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);          Last chunk stays on CPU 0 (the master)
96                          mesh_p->FaceElements->Color[i0]=i0;          The three columns of integers (Id, gDOF, Tag) are gathered into a single array tempInts and sent together in a single MPI message
97                          for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {          */
98                               fscanf(fileHandle_p, " %d",  
99                                  &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);          if (mpi_info->rank == 0) {  /* Master */
100                          }   /* for i1 */              for (;;) {            /* Infinite loop */
101                          fscanf(fileHandle_p, "\n");                  #pragma omp parallel for private (i0) schedule(static)
102                        } /* for i0 */                  for (i0=0; i0<chunkSize*3+1; i0++) tempInts[i0] = -1;
103                     }                  
104                  }                  #pragma omp parallel for private (i0) schedule(static)
105               }                  for (i0=0; i0<chunkSize*numDim; i0++) tempCoords[i0] = -1.0;
106                    
107                    chunkNodes = 0;
108                    for (i1=0; i1<chunkSize; i1++) {
109                        if (totalNodes >= numNodes) break;    /* End of inner loop */
110                        if (1 == numDim) {
111                            scan_ret = fscanf(fileHandle_p, "%d %d %d %le\n",
112                                                &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
113                                                &tempCoords[i1*numDim+0]);
114                            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
115                        }
116                        if (2 == numDim) {
117                            scan_ret = fscanf(fileHandle_p, "%d %d %d %le %le\n",
118                                                &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
119                                                &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1]);
120                            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
121                        }
122                        if (3 == numDim) {
123                            scan_ret = fscanf(fileHandle_p, "%d %d %d %le %le %le\n",
124                                                &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
125                                                &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1], &tempCoords[i1*numDim+2]);
126                            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
127                        }
128                        totalNodes++; /* When do we quit the infinite loop? */
129                        chunkNodes++; /* How many nodes do we actually have in this chunk? It may be smaller than chunkSize. */
130                    }
131                    if (chunkNodes > chunkSize) {
132                        Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: error reading chunks of mesh, data too large for message size");
133                        return NULL;
134                    }
135                    #ifdef PASO_MPI
136                        /* Eventually we'll send chunkSize nodes to each CPU numbered 1 ... mpi_info->size-1, here goes one of them */
137                        if (nextCPU < mpi_info->size) {
138                            tempInts[chunkSize*3] = chunkNodes;   /* The message has one more int to send chunkNodes */
139                            MPI_Send(tempInts, chunkSize*3+1, MPI_INT, nextCPU, 81720, mpi_info->comm);
140                            MPI_Send(tempCoords, chunkSize*numDim, MPI_DOUBLE, nextCPU, 81721, mpi_info->comm);
141                        }
142                    #endif
143                    nextCPU++;
144                    /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
145                    if (nextCPU > mpi_info->size) break; /* End infinite loop */
146                } /* Infinite loop */
147            }   /* End master */
148            else {  /* Worker */
149                #ifdef PASO_MPI
150                    /* Each worker receives two messages */
151                    MPI_Status status;
152                    MPI_Recv(tempInts, chunkSize*3+1, MPI_INT, 0, 81720, mpi_info->comm, &status);
153                    MPI_Recv(tempCoords, chunkSize*numDim, MPI_DOUBLE, 0, 81721, mpi_info->comm, &status);
154                    chunkNodes = tempInts[chunkSize*3];   /* How many nodes are in this workers chunk? */
155                #endif
156            }   /* Worker */
157    
158            /* Copy node data from tempMem to mesh_p */
159            Finley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes);
160            
161            if (Finley_noError()) {
162                #pragma omp parallel for private (i0, i1) schedule(static)
163                for (i0=0; i0<chunkNodes; i0++) {
164                    mesh_p->Nodes->Id[i0]               = tempInts[0+i0];
165                    mesh_p->Nodes->globalDegreesOfFreedom[i0]       = tempInts[chunkSize+i0];
166                    mesh_p->Nodes->Tag[i0]              = tempInts[chunkSize*2+i0];
167                    for (i1=0; i1<numDim; i1++) {
168                        mesh_p->Nodes->Coordinates[INDEX2(i1,i0,numDim)]  = tempCoords[i0*numDim+i1];
169                    }
170                }
171          }          }
172          /* get the Contact face element */          TMPMEMFREE(tempInts);
173          if (Finley_noError()) {          TMPMEMFREE(tempCoords);
174               fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);      }
175               contactTypeID=Finley_RefElement_getTypeId(element_type);  
176               if (contactTypeID==NoType) {      /* ***********************************  read elements ******************************************************************/
177                 sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);      if (Finley_noError()) {
178                 Finley_setError(VALUE_ERROR,error_msg);  
179               } else {          /* Read the element typeID */
180                 mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);          if (mpi_info->rank == 0) {
181                 if (Finley_noError()) {              scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
182                     Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
183                     if (Finley_noError()) {              typeID=Finley_ReferenceElement_getTypeId(element_type);
184                        mesh_p->ContactElements->minColor=0;          }
185                        mesh_p->ContactElements->maxColor=numEle-1;          #ifdef PASO_MPI
186                        for (i0 = 0; i0 < numEle; i0++) {              if (mpi_info->size > 1) {
187                          fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);                  int temp1[2], mpi_error;
188                          mesh_p->ContactElements->Color[i0]=i0;                  temp1[0] = (int) typeID;
189                          for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {                  temp1[1] = numEle;
190                              fscanf(fileHandle_p, " %d",                  mpi_error = MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
191                                 &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);                  if (mpi_error != MPI_SUCCESS) {
192                          }   /* for i1 */                      Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
193                          fscanf(fileHandle_p, "\n");                      return NULL;
194                        } /* for i0 */                  }
195                    }                  typeID = (ElementTypeId) temp1[0];
196                 }                  numEle = temp1[1];
197               }              }
198          }            #endif
199          /* get the nodal element */          if (typeID==NoRef) {
200          if (Finley_noError()) {              sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
201               fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);              Finley_setError(VALUE_ERROR, error_msg);
202               pointTypeID=Finley_RefElement_getTypeId(element_type);            }
203               if (pointTypeID==NoType) {      }
204                 sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type);  
205                 Finley_setError(VALUE_ERROR,error_msg);      /* Allocate the ElementFile */
206               }      if (Finley_noError()) {
207               mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);          refElements= Finley_ReferenceElementSet_alloc(typeID,mesh_p->order, mesh_p->reduced_order);
208               if (Finley_noError()) {          mesh_p->Elements=Finley_ElementFile_alloc(refElements, mpi_info);
209                  Finley_ElementFile_allocTable(mesh_p->Points, numEle);          numNodes = mesh_p->Elements->referenceElementSet->numNodes; /* New meaning for numNodes: num nodes per element */
210                  if (Finley_noError()) {      }
211                     mesh_p->Points->minColor=0;  
212                     mesh_p->Points->maxColor=numEle-1;      /* *************************** Read the element data *******************************************************************/
213                     for (i0 = 0; i0 < numEle; i0++) {      if (Finley_noError()) {
214                       fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);  
215                       mesh_p->Points->Color[i0]=i0;          int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1, chunkEle=0;
216                       for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {          int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
217                           fscanf(fileHandle_p, " %d",          /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
218                              &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);          if (mpi_info->rank == 0) {  /* Master */
219                       }  /* for i1 */              for (;;) {            /* Infinite loop */
220                       fscanf(fileHandle_p, "\n");                  #pragma omp parallel for private (i0) schedule(static)
221                     } /* for i0 */                  for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
222                  }          
223               }                  chunkEle = 0;
224                    for (i0=0; i0<chunkSize; i0++) {
225                        if (totalEle >= numEle) break; /* End inner loop */
226                        scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
227                        FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
228                        for (i1 = 0; i1 < numNodes; i1++) {
229                            scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
230                            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
231                        }
232                        scan_ret = fscanf(fileHandle_p, "\n");
233                        FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
234                        totalEle++;
235                        chunkEle++;
236                    }
237                    #ifdef PASO_MPI
238                        /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
239                        if (nextCPU < mpi_info->size) {
240                            tempInts[chunkSize*(2+numNodes)] = chunkEle;
241                            MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81722, mpi_info->comm);
242                        }
243                    #endif
244                    nextCPU++;
245                    /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
246                    if (nextCPU > mpi_info->size) break; /* End infinite loop */
247                } /* Infinite loop */
248            }   /* End master */
249            else {  /* Worker */
250                #ifdef PASO_MPI
251                    /* Each worker receives one message */
252                    MPI_Status status;
253                    MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81722, mpi_info->comm, &status);
254                    chunkEle = tempInts[chunkSize*(2+numNodes)];
255                #endif
256            }   /* Worker */
257    
258            
259            Finley_ElementFile_allocTable(mesh_p->Elements, chunkEle);
260    
261        
262            /* Copy Element data from tempInts to mesh_p */
263            if (Finley_noError()) {
264    
265                mesh_p->Elements->minColor=0;
266                mesh_p->Elements->maxColor=chunkEle-1;
267                #pragma omp parallel for private (i0, i1) schedule(static)
268                for (i0=0; i0<chunkEle; i0++) {
269                    mesh_p->Elements->Id[i0]    = tempInts[i0*(2+numNodes)+0];
270                    mesh_p->Elements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];
271                    mesh_p->Elements->Owner[i0]  =mpi_info->rank;
272                    mesh_p->Elements->Color[i0] = i0;
273                    for (i1 = 0; i1 < numNodes; i1++) {
274                        mesh_p->Elements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
275                    }
276                }
277            }
278    
279            TMPMEMFREE(tempInts);
280        }
281        /* ******************** end of Read the element data ******************************************************/
282    
283        /* ********************* read face elements ***************************************************************/
284        if (Finley_noError()) {
285            /* Read the element typeID */
286    
287            if (mpi_info->rank == 0) {
288                 scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
289                 FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
290                 typeID=Finley_ReferenceElement_getTypeId(element_type);
291            }
292            #ifdef PASO_MPI
293                if (mpi_info->size > 1) {
294                    int temp1[2];
295                    temp1[0] = (int) typeID;
296                    temp1[1] = numEle;
297                    MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
298                    typeID = (ElementTypeId) temp1[0];
299                    numEle = temp1[1];
300                }
301            #endif
302            if (typeID==NoRef) {
303                sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
304                Finley_setError(VALUE_ERROR, error_msg);
305          }          }
306          /* get the name tags */          if (Finley_noError()) {
307          if (Finley_noError()) {              /* Allocate the ElementFile */
308             if (feof(fileHandle_p) == 0) {              refFaceElements= Finley_ReferenceElementSet_alloc(typeID,mesh_p->order, mesh_p->reduced_order);
309                fscanf(fileHandle_p, "%s\n", name);              mesh_p->FaceElements=Finley_ElementFile_alloc(refFaceElements, mpi_info);
310                while (feof(fileHandle_p) == 0) {              numNodes = mesh_p->FaceElements->referenceElementSet->numNodes; /* New meaning for numNodes: num nodes per element */
311                     fscanf(fileHandle_p, "%s %d\n", name, &tag_key);          }
312                     Finley_Mesh_addTagMap(mesh_p,name,tag_key);  
313                }      }
314             }  
315        /* ********************** Read the face element data ********************************************************* */
316    
317        if (Finley_noError()) {
318            int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1, chunkEle=0;
319            int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
320            /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
321            if (mpi_info->rank == 0) {  /* Master */
322                for (;;) {            /* Infinite loop */
323                    #pragma omp parallel for private (i0) schedule(static)
324                    for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
325            
326                    chunkEle = 0;
327                    for (i0=0; i0<chunkSize; i0++) {
328                        if (totalEle >= numEle) break; /* End inner loop */
329                        scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
330                        FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
331                        for (i1 = 0; i1 < numNodes; i1++) {
332                            scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
333                            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
334                        }
335                        scan_ret = fscanf(fileHandle_p, "\n");
336                        FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
337                        totalEle++;
338                        chunkEle++;
339                    }
340                    #ifdef PASO_MPI
341                        /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
342                        if (nextCPU < mpi_info->size) {
343                            tempInts[chunkSize*(2+numNodes)] = chunkEle;
344                            MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81723, mpi_info->comm);
345                        }
346                    #endif
347                    nextCPU++;
348                    /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
349                    if (nextCPU > mpi_info->size) break; /* End infinite loop */
350                } /* Infinite loop */
351            }   /* End master */
352            else {  /* Worker */
353                #ifdef PASO_MPI
354                    /* Each worker receives one message */
355                    MPI_Status status;
356                    MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81723, mpi_info->comm, &status);
357                    chunkEle = tempInts[chunkSize*(2+numNodes)];
358                    #endif
359            }   /* Worker */
360            
361            Finley_ElementFile_allocTable(mesh_p->FaceElements, chunkEle);
362            
363            if (Finley_noError()) {
364                /* Copy Element data from tempInts to mesh_p */
365                
366                mesh_p->FaceElements->minColor=0;
367                mesh_p->FaceElements->maxColor=chunkEle-1;
368                #pragma omp parallel for private (i0, i1)
369                for (i0=0; i0<chunkEle; i0++) {
370                    mesh_p->FaceElements->Id[i0]    = tempInts[i0*(2+numNodes)+0];
371                    mesh_p->FaceElements->Tag[i0]   = tempInts[i0*(2+numNodes)+1];
372                    mesh_p->FaceElements->Owner[i0]  =mpi_info->rank;
373                    mesh_p->FaceElements->Color[i0] = i0;
374                    for (i1 = 0; i1 < numNodes; i1++) {
375                        mesh_p->FaceElements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
376                    }
377                }
378          }          }
      }  
      /* close file */  
      fclose(fileHandle_p);  
     
      /*   resolve id's : */  
      /* rearrange elements: */  
     
      if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);  
      if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);  
     
      /* that's it */  
      #ifdef Finley_TRACE  
      printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);  
      #endif  
   }  
   
   /* that's it */  
   if (! Finley_noError()) {  
        Finley_Mesh_free(mesh_p);  
   }  
   Paso_MPIInfo_free( mpi_info );  
   return mesh_p;  
 }  
379    
380  Finley_Mesh* Finley_Mesh_read_MPI(char* fname,index_t order, index_t reduced_order,  bool_t optimize)          TMPMEMFREE(tempInts);
381        }
382  {      /* ************************************* end of Read the face element data *************************************** */
383    
384    Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );  
385    dim_t numNodes, numDim, numEle, i0, i1;      /* ************************************* read contact elements ************************************************** */
386    index_t tag_key;  
387    Finley_Mesh *mesh_p=NULL;      /* Read the element typeID */
388    char name[LenString_MAX],element_type[LenString_MAX],frm[20];      if (Finley_noError()) {
389    char error_msg[LenErrorMsg_MAX];          if (mpi_info->rank == 0) {
390    double time0=Finley_timer();              scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
391    FILE *fileHandle_p = NULL;              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
392    ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;              typeID=Finley_ReferenceElement_getTypeId(element_type);
393              }
394    Finley_resetError();          #ifdef PASO_MPI
395                if (mpi_info->size > 1) {
396    if (mpi_info->rank == 0) {                  int temp1[2];
397       /* get file handle */                  temp1[0] = (int) typeID;
398       fileHandle_p = fopen(fname, "r");                  temp1[1] = numEle;
399       if (fileHandle_p==NULL) {                  MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
400         sprintf(error_msg,"Finley_Mesh_read: Opening file %s for reading failed.",fname);                  typeID = (ElementTypeId) temp1[0];
401         Finley_setError(IO_ERROR,error_msg);                  numEle = temp1[1];
402         Paso_MPIInfo_free( mpi_info );              }
403         return NULL;          #endif
404       }          if (typeID==NoRef) {
405                  sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
406       /* read header */              Finley_setError(VALUE_ERROR, error_msg);
407       sprintf(frm,"%%%d[^\n]",LenString_MAX-1);           }
      fscanf(fileHandle_p, frm, name);  
     
      /* get the nodes */  
     
      fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);  
   }  
   
 #ifdef PASO_MPI  
   /* MPI Broadcast numDim, numNodes, name */  
   if (mpi_info->size > 0) {  
     int temp1[3], error_code;  
     temp1[0] = numDim;  
     temp1[1] = numNodes;  
     temp1[2] = strlen(name) + 1;  
     error_code = MPI_Bcast (temp1, 3, MPI_INT,  0, mpi_info->comm);  
     if (error_code != MPI_SUCCESS) {  
       Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of temp1 failed");  
       return NULL;  
     }  
     numDim = temp1[0];  
     numNodes = temp1[1];  
     error_code = MPI_Bcast (name, temp1[2], MPI_CHAR, 0, mpi_info->comm);  
     if (error_code != MPI_SUCCESS) {  
       Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of name failed");  
       return NULL;  
408      }      }
409    }  
410  #endif      if (Finley_noError()) {
411    printf("ksteube CPU=%d/%d name='%s' numDim=%d numNodes=%d\n", mpi_info->rank, mpi_info->size, name, numDim, numNodes);          /* Allocate the ElementFile */
412            refContactElements= Finley_ReferenceElementSet_alloc(typeID,mesh_p->order, mesh_p->reduced_order);
413       /* allocate mesh */          mesh_p->ContactElements=Finley_ElementFile_alloc(refContactElements, mpi_info);
414       mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);          numNodes = mesh_p->ContactElements->referenceElementSet->numNodes; /* New meaning for numNodes: num nodes per element */
415       if (Finley_noError()) {      }
416      int chunkSize = numNodes / mpi_info->size + 1, totalNodes=0, chunkNodes=0, nextCPU=1, mpi_error;      /* *************************** Read the contact element data ************************************************* */
417      int *tempInts = TMPMEMALLOC(numNodes*3+1, index_t);      if (Finley_noError()) {
418      double *tempCoords = TMPMEMALLOC(numNodes*numDim, double);          int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1, chunkEle=0;
419            int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
420      /*          /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
421        Read a chunk of nodes, send to worker CPU if necessary, copy chunk into local mesh_p          if (mpi_info->rank == 0) {  /* Master */
422        It doesn't matter that a CPU has the wrong nodes, this is sorted out later              for (;;) {            /* Infinite loop */
423        First chunk sent to CPU 1, second to CPU 2, ...                  #pragma omp parallel for private (i0) schedule(static)
424        Last chunk stays on CPU 0 (the master)                  for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
425        The three columns of integers (Id, gDOF, Tag) are gathered into a single array tempInts and sent together in a single MPI message        
426      */                  chunkEle = 0;
427                    for (i0=0; i0<chunkSize; i0++) {
428      if (mpi_info->rank == 0) {  /* Master */                      if (totalEle >= numEle) break; /* End inner loop */
429        for (;;) {            /* Infinite loop */                      scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
430          chunkNodes = 0;                      FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
431          for (i0=0; i0<numNodes*3; i0++) tempInts[i0] = -1;                      for (i1 = 0; i1 < numNodes; i1++) {
432          for (i0=0; i0<numNodes*numDim; i0++) tempCoords[i0] = -1.0;                          scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
433          for (i1=0; i1<chunkSize; i1++) {                          FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
434            if (totalNodes >= numNodes) break;                      }
435                if (1 == numDim)                      scan_ret = fscanf(fileHandle_p, "\n");
436          fscanf(fileHandle_p, "%d %d %d %le\n",                      FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
437            &tempInts[0+i1], &tempInts[numNodes+i1], &tempInts[numNodes*2+i1],                      totalEle++;
438            &tempCoords[i1*numDim+0]);                      chunkEle++;
439                if (2 == numDim)                  }
440          fscanf(fileHandle_p, "%d %d %d %le %le\n",                  #ifdef PASO_MPI
441            &tempInts[0+i1], &tempInts[numNodes+i1], &tempInts[numNodes*2+i1],                      /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
442            &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1]);                      if (nextCPU < mpi_info->size) {
443                if (3 == numDim)                          tempInts[chunkSize*(2+numNodes)] = chunkEle;
444          fscanf(fileHandle_p, "%d %d %d %le %le %le\n",                          MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81724, mpi_info->comm);
445            &tempInts[0+i1], &tempInts[numNodes+i1], &tempInts[numNodes*2+i1],                      }
446            &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1], &tempCoords[i1*numDim+2]);                  #endif
447            totalNodes++;                  nextCPU++;
448            chunkNodes++;                  /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
449          }                  if (nextCPU > mpi_info->size) break; /* End infinite loop */
450          /* Eventually we'll send chunk of nodes to each CPU numbered 1 ... mpi_info->size-1, here goes one of them */              } /* Infinite loop */
451          if (nextCPU < mpi_info->size) {          }   /* End master */
452  #ifdef PASO_MPI          else {  /* Worker */
453            tempInts[numNodes*3] = chunkNodes;              #ifdef PASO_MPI
454            mpi_error = MPI_Send(tempInts, numNodes*3+1, MPI_INT, nextCPU, 8727, mpi_info->comm);                  /* Each worker receives one message */
455            if ( mpi_error != MPI_SUCCESS ) {                  MPI_Status status;
456                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts failed");                  MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81724, mpi_info->comm, &status);
457                  return NULL;                  chunkEle = tempInts[chunkSize*(2+numNodes)] ;
458            }              #endif
459            mpi_error = MPI_Send(tempCoords, numNodes*numDim, MPI_DOUBLE, nextCPU, 8728, mpi_info->comm);          }   /* Worker */
460            if ( mpi_error != MPI_SUCCESS ) {  
461                  Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempCoords failed");          /* Copy Element data from tempInts to mesh_p */
462                  return NULL;           Finley_ElementFile_allocTable(mesh_p->ContactElements, chunkEle);
463            }          
464  #endif           if (Finley_noError()) {        
465            nextCPU++;               mesh_p->ContactElements->minColor=0;
466          }               mesh_p->ContactElements->maxColor=chunkEle-1;
467          if (totalNodes >= numNodes) break;               #pragma omp parallel for private (i0, i1)
468        } /* Infinite loop */               for (i0=0; i0<chunkEle; i0++) {
469      }   /* End master */                   mesh_p->ContactElements->Id[i0] = tempInts[i0*(2+numNodes)+0];
470      else {  /* Worker */                   mesh_p->ContactElements->Tag[i0]    = tempInts[i0*(2+numNodes)+1];
471  #ifdef PASO_MPI                   mesh_p->ContactElements->Owner[i0]  =mpi_info->rank;
472        /* Each worker receives two messages */                   mesh_p->ContactElements->Color[i0] = i0;
473        MPI_Status status;                   for (i1 = 0; i1 < numNodes; i1++) {
474        mpi_error = MPI_Recv(tempInts, numNodes*3+1, MPI_INT, 0, 8727, mpi_info->comm, &status);                       mesh_p->ContactElements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
475        if ( mpi_error != MPI_SUCCESS ) {                   }
476              Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts failed");               }
             return NULL;  
       }  
       mpi_error = MPI_Recv(tempCoords, numNodes*numDim, MPI_DOUBLE, 0, 8728, mpi_info->comm, &status);  
       if ( mpi_error != MPI_SUCCESS ) {  
             Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempCoords failed");  
             return NULL;  
       }  
       chunkNodes = tempInts[numNodes*3];  
 #endif  
     }   /* Worker */  
 printf("ksteube CPU=%d/%d name='%s' numDim=%d numNodes=%d chunkNodes=%d\n", mpi_info->rank, mpi_info->size, name, numDim, numNodes, chunkNodes);  
 #if 0  
         printf("ksteube tempInts totalNodes=%d:\n", totalNodes);  
         for (i0=0; i0<numNodes*3; i0++) {  
           printf(" %2d", tempInts[i0]);  
           if (i0%numNodes==numNodes-1) printf("\n");  
         }  
         printf("ksteube tempCoords:\n");  
         for (i0=0; i0<chunkNodes*numDim; i0++) {  
           printf(" %20.15e", tempCoords[i0]);  
           if (i0%numDim==numDim-1) printf("\n");  
         }  
 #endif  
   
     /* Copy tempMem into mesh_p */  
         Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);  
     for (i0=0; i0<numNodes; i0++) {  
       mesh_p->Nodes->Id[i0]                 = tempInts[0+i0];  
       mesh_p->Nodes->globalDegreesOfFreedom[i0]     = tempInts[numNodes+i0];  
       mesh_p->Nodes->Tag[i0]                = tempInts[numNodes*2+i0];  
       for (i1=0; i1<numDim; i1++) {  
         mesh_p->Nodes->Coordinates[INDEX2(i1,i0,numDim)]    = tempCoords[i0*numDim+i1];  
       }  
477          }          }
478            TMPMEMFREE(tempInts);
479        } /* end of Read the contact element data */
480    
481      TMPMEMFREE(tempInts);      /* ********************************* read nodal elements ****************************************************** */
     TMPMEMFREE(tempCoords);  
482    
483          return mesh_p; /* ksteube temp for debugging */      /* *******************************  Read the element typeID */
484    
485          /* read elements */      if (Finley_noError()) {
486          if (Finley_noError()) {          if (mpi_info->rank == 0) {
487                  scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
488             fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);              FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
489             typeID=Finley_RefElement_getTypeId(element_type);              typeID=Finley_ReferenceElement_getTypeId(element_type);
490             if (typeID==NoType) {          }
491               sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s",element_type);          #ifdef PASO_MPI
492               Finley_setError(VALUE_ERROR,error_msg);              if (mpi_info->size > 1) {
493             } else {                  int temp1[2];
494               /* read the elements */                  temp1[0] = (int) typeID;
495               mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);                  temp1[1] = numEle;
496               if (Finley_noError()) {                  MPI_Bcast (temp1, 2, MPI_INT,  0, mpi_info->comm);
497                   Finley_ElementFile_allocTable(mesh_p->Elements, numEle);                  typeID = (ElementTypeId) temp1[0];
498                   mesh_p->Elements->minColor=0;                  numEle = temp1[1];
499                   mesh_p->Elements->maxColor=numEle-1;              }
500                   if (Finley_noError()) {          #endif
501                      for (i0 = 0; i0 < numEle; i0++) {          if (typeID==NoRef) {
502                        fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);              sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
503                        mesh_p->Elements->Color[i0]=i0;              Finley_setError(VALUE_ERROR, error_msg);
504                        for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {           }
505                             fscanf(fileHandle_p, " %d",      }
506                                &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);  
507                        } /* for i1 */      if (Finley_noError()) {
508                        fscanf(fileHandle_p, "\n");          /* Allocate the ElementFile */
509                      } /* for i0 */          refPoints= Finley_ReferenceElementSet_alloc(typeID,mesh_p->order, mesh_p->reduced_order);
510                   }          mesh_p->Points=Finley_ElementFile_alloc(refPoints, mpi_info);
511               }          numNodes = mesh_p->Points->referenceElementSet->numNodes; /* New meaning for numNodes: num nodes per element */
512            }      }
513          }  
514          /* get the face elements */      /**********************************  Read the nodal element data **************************************************/
515          if (Finley_noError()) {      if (Finley_noError()) {
516               fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);          int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1, chunkEle=0;
517               faceTypeID=Finley_RefElement_getTypeId(element_type);          int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
518               if (faceTypeID==NoType) {          /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
519                 sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s for face elements",element_type);          if (mpi_info->rank == 0) {  /* Master */
520                 Finley_setError(VALUE_ERROR,error_msg);              for (;;) {            /* Infinite loop */
521               } else {                  #pragma omp parallel for private (i0) schedule(static)
522                  mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);                  for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
523                  if (Finley_noError()) {          
524                     Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);                  chunkEle = 0;
525                     if (Finley_noError()) {                  for (i0=0; i0<chunkSize; i0++) {
526                        mesh_p->FaceElements->minColor=0;                      if (totalEle >= numEle) break; /* End inner loop */
527                        mesh_p->FaceElements->maxColor=numEle-1;                      scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
528                        for (i0 = 0; i0 < numEle; i0++) {                      FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
529                          fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);                      for (i1 = 0; i1 < numNodes; i1++) {
530                          mesh_p->FaceElements->Color[i0]=i0;                          scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
531                          for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {                          FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
532                               fscanf(fileHandle_p, " %d",                      }
533                                  &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);                      scan_ret = fscanf(fileHandle_p, "\n");
534                          }   /* for i1 */                      FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
535                          fscanf(fileHandle_p, "\n");                      totalEle++;
536                        } /* for i0 */                      chunkEle++;
537                     }                  }
538                  }                  #ifdef PASO_MPI
539               }                      /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
540          }                      if (nextCPU < mpi_info->size) {
541          /* get the Contact face element */                          tempInts[chunkSize*(2+numNodes)] = chunkEle;
542          if (Finley_noError()) {                          MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81725, mpi_info->comm);
543               fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);                      }
544               contactTypeID=Finley_RefElement_getTypeId(element_type);                  #endif
545               if (contactTypeID==NoType) {                  nextCPU++;
546                 sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);                  /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
547                 Finley_setError(VALUE_ERROR,error_msg);                  if (nextCPU > mpi_info->size) break; /* End infinite loop */
548               } else {              } /* Infinite loop */
549                 mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);          }   /* End master */
550                 if (Finley_noError()) {          else {  /* Worker */
551                     Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);              #ifdef PASO_MPI
552                     if (Finley_noError()) {                  /* Each worker receives one message */
553                        mesh_p->ContactElements->minColor=0;                  MPI_Status status;
554                        mesh_p->ContactElements->maxColor=numEle-1;                  MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81725, mpi_info->comm, &status);
555                        for (i0 = 0; i0 < numEle; i0++) {                  chunkEle = tempInts[chunkSize*(2+numNodes)];
556                          fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);              #endif
557                          mesh_p->ContactElements->Color[i0]=i0;          }   /* Worker */
558                          for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {  
559                              fscanf(fileHandle_p, " %d",          /* Copy Element data from tempInts to mesh_p */
560                                 &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);          Finley_ElementFile_allocTable(mesh_p->Points, chunkEle);
561                          }   /* for i1 */          
562                          fscanf(fileHandle_p, "\n");          if (Finley_noError()) {
563                        } /* for i0 */              mesh_p->Points->minColor=0;
564                    }              mesh_p->Points->maxColor=chunkEle-1;
565                 }              #pragma omp parallel for private (i0, i1) schedule(static)
566               }              for (i0=0; i0<chunkEle; i0++) {
567          }                    mesh_p->Points->Id[i0]  = tempInts[i0*(2+numNodes)+0];
568          /* get the nodal element */                  mesh_p->Points->Tag[i0] = tempInts[i0*(2+numNodes)+1];
569          if (Finley_noError()) {                  mesh_p->Points->Owner[i0]  =mpi_info->rank;
570               fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);                  mesh_p->Points->Color[i0] = i0;
571               pointTypeID=Finley_RefElement_getTypeId(element_type);                  for (i1 = 0; i1 < numNodes; i1++) {
572               if (pointTypeID==NoType) {                      mesh_p->Points->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
573                 sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type);                  }
574                 Finley_setError(VALUE_ERROR,error_msg);              }
575               }          }
576               mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);  
577               if (Finley_noError()) {          TMPMEMFREE(tempInts);
578                  Finley_ElementFile_allocTable(mesh_p->Points, numEle);      } /* ******************************** end of Read the nodal element data ************************************************************* */
579                  if (Finley_noError()) {  
580                     mesh_p->Points->minColor=0;      
581                     mesh_p->Points->maxColor=numEle-1;      /******************  get the name tags *****************************************/
582                     for (i0 = 0; i0 < numEle; i0++) {      if (Finley_noError()) {
583                       fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);          char *remainder=0, *ptr;
584                       mesh_p->Points->Color[i0]=i0;          size_t len=0;
585                       for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {          #ifdef PASO_MPI
586                           fscanf(fileHandle_p, " %d",              int len_i;
587                              &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);          #endif
588                       }  /* for i1 */          int tag_key;
589                       fscanf(fileHandle_p, "\n");          if (mpi_info->rank == 0) {  /* Master */
590                     } /* for i0 */              /* Read the word 'Tag' */
591                  }              if (! feof(fileHandle_p)) {
592               }                  scan_ret = fscanf(fileHandle_p, "%s\n", name);
593          }                  FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
594          /* get the name tags */              }
595          if (Finley_noError()) {  
596             if (feof(fileHandle_p) == 0) {              #if defined(_WIN32)  /* windows ftell lies on unix formatted text files */
597                fscanf(fileHandle_p, "%s\n", name);                  remainder = NULL;
598                while (feof(fileHandle_p) == 0) {                  len=0;
599                     fscanf(fileHandle_p, "%s %d\n", name, &tag_key);                  while (1)
600                     Finley_Mesh_addTagMap(mesh_p,name,tag_key);                  {
601                }                      size_t malloc_chunk = 1024;
602             }                      size_t buff_size = 0;
603          }                      int ch;
604       }  
605       /* close file */                      ch = fgetc(fileHandle_p);
606       fclose(fileHandle_p);                      if( ch == '\r' )
607                          {
608       /*   resolve id's : */                          continue;
609       /* rearrange elements: */                      }
610                          if( len+1 > buff_size )
611       if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);                      {
612       if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);                          TMPMEMREALLOC(remainder,remainder,buff_size+malloc_chunk,char);
613                          }
614       /* that's it */                      if( ch == EOF )
615       #ifdef Finley_TRACE                      {
616       printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);                          /* hit EOF */
617       #endif                          remainder[len] = (char)0;
618                            break;
619    /* that's it */                      }
620    if (! Finley_noError()) {                      remainder[len] = (char)ch;
621         Finley_Mesh_free(mesh_p);                      len++;
622    }                  }
623    Paso_MPIInfo_free( mpi_info );              #else
624    return mesh_p;                  /* Read rest of file in one chunk, after using seek to find length */
625                    {
626                        long cur_pos, end_pos;
627    
628                        cur_pos = ftell(fileHandle_p);
629                        fseek(fileHandle_p, 0L, SEEK_END);
630                        end_pos = ftell(fileHandle_p);
631                        fseek(fileHandle_p, (long)cur_pos, SEEK_SET);
632                        remainder = TMPMEMALLOC(end_pos-cur_pos+1, char);
633                        if (! feof(fileHandle_p))
634                        {
635                            scan_ret = fread(remainder, (size_t) end_pos-cur_pos,
636                                                 sizeof(char), fileHandle_p);
637    
638                            FSCANF_CHECK(scan_ret, "Finley_Mesh_read")
639                            remainder[end_pos-cur_pos] = 0;
640                        }
641                    }
642                #endif
643                len = strlen(remainder);
644                while ((len>1) && isspace(remainder[--len])) {remainder[len]=0;}
645                len = strlen(remainder);
646                TMPMEMREALLOC(remainder,remainder,len+1,char);
647            } /* Master */
648            #ifdef PASO_MPI
649    
650                len_i=(int) len;
651                MPI_Bcast (&len_i, 1, MPI_INT,  0, mpi_info->comm);
652                len=(size_t) len_i;
653                if (mpi_info->rank != 0) {
654                    remainder = TMPMEMALLOC(len+1, char);
655                    remainder[0] = 0;
656                }
657                if (MPI_Bcast (remainder, len+1, MPI_CHAR,  0, mpi_info->comm) !=
658                        MPI_SUCCESS)
659                    Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of remainder failed");
660            #endif
661    
662            if (remainder[0]) {
663                ptr = remainder;
664                do {
665                    sscanf(ptr, "%s %d\n", name, &tag_key);
666                    if (*name) Finley_Mesh_addTagMap(mesh_p,name,tag_key);
667                    ptr++;
668                } while(NULL != (ptr = strchr(ptr, '\n')) && *ptr);
669            }
670            if (remainder)
671                TMPMEMFREE(remainder);
672        }
673    
674        /* close file */
675        if (mpi_info->rank == 0) fclose(fileHandle_p);
676    
677        /*   resolve id's : */
678        /* rearrange elements: */
679        if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
680        if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
681    
682        /* that's it */
683        if (! Finley_noError()) {
684            Finley_Mesh_free(mesh_p);
685        }
686        /* free up memory */
687        Finley_ReferenceElementSet_dealloc(refPoints);
688        Finley_ReferenceElementSet_dealloc(refContactElements);
689        Finley_ReferenceElementSet_dealloc(refFaceElements);
690        Finley_ReferenceElementSet_dealloc(refElements);
691        Paso_MPIInfo_free( mpi_info );
692        return mesh_p;
693  }  }
694    

Legend:
Removed from v.1387  
changed lines
  Added in v.2754

  ViewVC Help
Powered by ViewVC 1.1.26