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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1732 - (hide annotations)
Wed Aug 27 05:47:03 2008 UTC (11 years, 3 months ago) by ksteube
File MIME type: text/plain
File size: 33821 byte(s)
MPI read of finley mesh file working except for a bug related to mesh_prepare

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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26