/[escript]/trunk/finley/src/Mesh_read.c
ViewVC logotype

Diff of /trunk/finley/src/Mesh_read.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

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

  ViewVC Help
Powered by ViewVC 1.1.26