/[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 1628 - (hide annotations)
Fri Jul 11 13:12:46 2008 UTC (11 years, 3 months ago) by phornby
File MIME type: text/plain
File size: 28803 byte(s)

Merge in /branches/windows_from_1456_trunk_1620_merged_in branch.

You will find a preserved pre-merge trunk in tags under tags/trunk_at_1625.
That will be useful for diffing & checking on my stupidity.

Here is a list of the conflicts and their resolution at this
point in time.


=================================================================================
(LLWS == looks like white space).

finley/src/Assemble_addToSystemMatrix.c - resolve to branch - unused var. may be wrong.....
finley/src/CPPAdapter/SystemMatrixAdapter.cpp - resolve to branch - LLWS
finley/src/CPPAdapter/MeshAdapter.cpp - resolve to branch - LLWS
paso/src/PCG.c - resolve to branch - unused var fixes.
paso/src/SolverFCT.c - resolve to branch - LLWS
paso/src/FGMRES.c - resolve to branch - LLWS
paso/src/Common.h - resolve to trunk version. It's omp.h's include... not sure it's needed,
but for the sake of saftey.....
paso/src/Functions.c - resolve to branch version, indentation/tab removal and return error
on bad unimplemented Paso_FunctionCall.
paso/src/SolverFCT_solve.c - resolve to branch version, unused vars
paso/src/SparseMatrix_MatrixVector.c - resolve to branch version, unused vars.
escript/src/Utils.cpp - resloved to branch, needs WinSock2.h
escript/src/DataExpanded.cpp - resolved to branch version - LLWS
escript/src/DataFactory.cpp - resolve to branch version
=================================================================================

This currently passes tests on linux (debian), but is not checked on windows or Altix yet.

This checkin is to make a trunk I can check out for windows to do tests on it.

Known outstanding problem is in the operator=() method of exceptions
causing warning messages on the intel compilers.

May the God of doughnuts have mercy on my soul.


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 phornby 1628 ElementTypeId typeID;
245 ksteube 1402
246 phornby 1628 /* these are in unimplemented code below */
247     #if 0
248     index_t tag_key;
249     ElementTypeId faceTypeID, contactTypeID, pointTypeID;
250     #endif
251    
252 ksteube 1360 Finley_resetError();
253    
254     if (mpi_info->rank == 0) {
255     /* get file handle */
256     fileHandle_p = fopen(fname, "r");
257     if (fileHandle_p==NULL) {
258     sprintf(error_msg,"Finley_Mesh_read: Opening file %s for reading failed.",fname);
259     Finley_setError(IO_ERROR,error_msg);
260     Paso_MPIInfo_free( mpi_info );
261     return NULL;
262     }
263 ksteube 1402
264 ksteube 1360 /* read header */
265     sprintf(frm,"%%%d[^\n]",LenString_MAX-1);
266     fscanf(fileHandle_p, frm, name);
267 ksteube 1402
268     /* get the number of nodes */
269 ksteube 1360 fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);
270     }
271    
272     #ifdef PASO_MPI
273     /* MPI Broadcast numDim, numNodes, name */
274     if (mpi_info->size > 0) {
275     int temp1[3], error_code;
276     temp1[0] = numDim;
277     temp1[1] = numNodes;
278     temp1[2] = strlen(name) + 1;
279     error_code = MPI_Bcast (temp1, 3, MPI_INT, 0, mpi_info->comm);
280     if (error_code != MPI_SUCCESS) {
281     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of temp1 failed");
282     return NULL;
283     }
284     numDim = temp1[0];
285     numNodes = temp1[1];
286     error_code = MPI_Bcast (name, temp1[2], MPI_CHAR, 0, mpi_info->comm);
287     if (error_code != MPI_SUCCESS) {
288     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of name failed");
289     return NULL;
290     }
291     }
292     #endif
293    
294     /* allocate mesh */
295     mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
296     if (Finley_noError()) {
297 phornby 1628 int chunkSize = numNodes / mpi_info->size + 1, totalNodes=0, chunkNodes=0, chunkEle=0, nextCPU=1;
298     #ifdef PASO_MPI
299     int mpi_error;
300     #endif
301 ksteube 1360 int *tempInts = TMPMEMALLOC(numNodes*3+1, index_t);
302     double *tempCoords = TMPMEMALLOC(numNodes*numDim, double);
303    
304     /*
305 ksteube 1402 Read a chunk of nodes, send to worker CPU if available, copy chunk into local mesh_p
306     It doesn't matter that a CPU has the wrong nodes for its elements, this is sorted out later
307 ksteube 1360 First chunk sent to CPU 1, second to CPU 2, ...
308     Last chunk stays on CPU 0 (the master)
309     The three columns of integers (Id, gDOF, Tag) are gathered into a single array tempInts and sent together in a single MPI message
310     */
311    
312     if (mpi_info->rank == 0) { /* Master */
313     for (;;) { /* Infinite loop */
314     chunkNodes = 0;
315 ksteube 1402 for (i0=0; i0<numNodes*3+1; i0++) tempInts[i0] = -1;
316 ksteube 1360 for (i0=0; i0<numNodes*numDim; i0++) tempCoords[i0] = -1.0;
317     for (i1=0; i1<chunkSize; i1++) {
318     if (totalNodes >= numNodes) break;
319     if (1 == numDim)
320     fscanf(fileHandle_p, "%d %d %d %le\n",
321     &tempInts[0+i1], &tempInts[numNodes+i1], &tempInts[numNodes*2+i1],
322     &tempCoords[i1*numDim+0]);
323     if (2 == numDim)
324     fscanf(fileHandle_p, "%d %d %d %le %le\n",
325     &tempInts[0+i1], &tempInts[numNodes+i1], &tempInts[numNodes*2+i1],
326     &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1]);
327     if (3 == numDim)
328     fscanf(fileHandle_p, "%d %d %d %le %le %le\n",
329     &tempInts[0+i1], &tempInts[numNodes+i1], &tempInts[numNodes*2+i1],
330     &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1], &tempCoords[i1*numDim+2]);
331     totalNodes++;
332     chunkNodes++;
333     }
334     /* Eventually we'll send chunk of nodes to each CPU numbered 1 ... mpi_info->size-1, here goes one of them */
335     if (nextCPU < mpi_info->size) {
336     #ifdef PASO_MPI
337     tempInts[numNodes*3] = chunkNodes;
338 ksteube 1402 /* ksteube The size of this message can and should be brought down to chunkNodes*3+1, must re-org tempInts */
339     mpi_error = MPI_Send(tempInts, numNodes*3+1, MPI_INT, nextCPU, 81720, mpi_info->comm);
340 ksteube 1360 if ( mpi_error != MPI_SUCCESS ) {
341     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts failed");
342     return NULL;
343     }
344 ksteube 1402 mpi_error = MPI_Send(tempCoords, numNodes*numDim, MPI_DOUBLE, nextCPU, 81721, mpi_info->comm);
345 ksteube 1360 if ( mpi_error != MPI_SUCCESS ) {
346     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempCoords failed");
347     return NULL;
348     }
349     #endif
350     nextCPU++;
351     }
352     if (totalNodes >= numNodes) break;
353     } /* Infinite loop */
354     } /* End master */
355     else { /* Worker */
356     #ifdef PASO_MPI
357     /* Each worker receives two messages */
358     MPI_Status status;
359 ksteube 1402 mpi_error = MPI_Recv(tempInts, numNodes*3+1, MPI_INT, 0, 81720, mpi_info->comm, &status);
360 ksteube 1360 if ( mpi_error != MPI_SUCCESS ) {
361     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts failed");
362     return NULL;
363     }
364 ksteube 1402 mpi_error = MPI_Recv(tempCoords, numNodes*numDim, MPI_DOUBLE, 0, 81721, mpi_info->comm, &status);
365 ksteube 1360 if ( mpi_error != MPI_SUCCESS ) {
366     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempCoords failed");
367     return NULL;
368     }
369     chunkNodes = tempInts[numNodes*3];
370     #endif
371     } /* Worker */
372 ksteube 1402
373 ksteube 1360 #if 0
374 ksteube 1402 /* Display the temp mem for debugging */
375 ksteube 1360 printf("ksteube tempInts totalNodes=%d:\n", totalNodes);
376     for (i0=0; i0<numNodes*3; i0++) {
377     printf(" %2d", tempInts[i0]);
378     if (i0%numNodes==numNodes-1) printf("\n");
379     }
380     printf("ksteube tempCoords:\n");
381     for (i0=0; i0<chunkNodes*numDim; i0++) {
382     printf(" %20.15e", tempCoords[i0]);
383     if (i0%numDim==numDim-1) printf("\n");
384     }
385     #endif
386    
387 ksteube 1402 printf("ksteube chunkNodes=%d numNodes=%d\n", chunkNodes, numNodes);
388     /* Copy node data from tempMem to mesh_p */
389     Finley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes);
390     if (Finley_noError()) {
391     for (i0=0; i0<chunkNodes; i0++) {
392     mesh_p->Nodes->Id[i0] = tempInts[0+i0];
393     mesh_p->Nodes->globalDegreesOfFreedom[i0] = tempInts[numNodes+i0];
394     mesh_p->Nodes->Tag[i0] = tempInts[numNodes*2+i0];
395     for (i1=0; i1<numDim; i1++) {
396     mesh_p->Nodes->Coordinates[INDEX2(i1,i0,numDim)] = tempCoords[i0*numDim+i1];
397     }
398     }
399 ksteube 1360 }
400    
401     TMPMEMFREE(tempInts);
402     TMPMEMFREE(tempCoords);
403    
404 ksteube 1402 /* read elements */
405 ksteube 1360
406 ksteube 1402 /* Read the element typeID */
407     if (Finley_noError()) {
408     if (mpi_info->rank == 0) {
409     fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
410     typeID=Finley_RefElement_getTypeId(element_type);
411     }
412     #ifdef PASO_MPI
413     if (mpi_info->size > 0) {
414     int temp1[3];
415     temp1[0] = typeID;
416     temp1[1] = numEle;
417     mpi_error = MPI_Bcast (temp1, 2, MPI_INT, 0, mpi_info->comm);
418     if (mpi_error != MPI_SUCCESS) {
419     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
420     return NULL;
421     }
422     typeID = temp1[0];
423     numEle = temp1[1];
424     }
425     #endif
426     if (typeID==NoType) {
427     sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
428     Finley_setError(VALUE_ERROR, error_msg);
429     }
430     }
431    
432     /* Read the element data */
433     mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
434     numNodes = mesh_p->Elements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
435    
436     if (Finley_noError()) {
437     int *tempInts = TMPMEMALLOC(numEle*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
438 phornby 1628 int chunkSize = numEle / mpi_info->size, totalEle=0, nextCPU=1;
439     #ifdef PASO_MPI
440     int mpi_error;
441     #endif
442 ksteube 1402 if (numEle % mpi_info->size != 0) chunkSize++; /* Remainder from numEle / mpi_info->size will be spread out one-per-CPU */
443     if (mpi_info->rank == 0) { /* Master */
444     for (;;) { /* Infinite loop */
445     chunkEle = 0;
446     for (i0=0; i0<numEle*(2+numNodes)+1; i0++) tempInts[i0] = -1;
447     for (i0=0; i0<chunkSize; i0++) {
448     if (totalEle >= numEle) break; /* End infinite loop */
449     fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
450     for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
451     fscanf(fileHandle_p, "\n");
452     totalEle++;
453     chunkEle++;
454     }
455     /* Eventually we'll send chunk of nodes to each CPU except 0 itself, here goes one of them */
456     if (nextCPU < mpi_info->size) {
457     #ifdef PASO_MPI
458     tempInts[numEle*(2+numNodes)] = chunkEle;
459     printf("ksteube CPU=%d/%d send to %d\n", mpi_info->rank, mpi_info->size, nextCPU);
460     mpi_error = MPI_Send(tempInts, numEle*(2+numNodes)+1, MPI_INT, nextCPU, 81722, mpi_info->comm);
461     if ( mpi_error != MPI_SUCCESS ) {
462     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Elements failed");
463     return NULL;
464     }
465     #endif
466     nextCPU++;
467     }
468     if (totalEle >= numEle) break; /* End infinite loop */
469     } /* Infinite loop */
470     } /* End master */
471     else { /* Worker */
472     #ifdef PASO_MPI
473     /* Each worker receives two messages */
474     MPI_Status status;
475     printf("ksteube CPU=%d/%d recv on %d\n", mpi_info->rank, mpi_info->size, mpi_info->rank);
476     mpi_error = MPI_Recv(tempInts, numEle*(2+numNodes)+1, MPI_INT, 0, 81722, mpi_info->comm, &status);
477     if ( mpi_error != MPI_SUCCESS ) {
478     Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for Elements failed");
479     return NULL;
480     }
481     chunkEle = tempInts[numEle*(2+numNodes)];
482     #endif
483     } /* Worker */
484     #if 1
485     /* Display the temp mem for debugging */
486     printf("ksteube tempInts numEle=%d chunkEle=%d AAA:\n", numEle, chunkEle);
487     for (i0=0; i0<numEle*(numNodes+2); i0++) {
488     printf(" %2d", tempInts[i0]);
489     if (i0%(numNodes+2)==numNodes+2-1) printf("\n");
490     }
491     #endif
492    
493     /* Copy Element data from tempInts to mesh_p */
494     Finley_ElementFile_allocTable(mesh_p->Elements, chunkEle);
495     mesh_p->Elements->minColor=0;
496     mesh_p->Elements->maxColor=chunkEle-1;
497     if (Finley_noError()) {
498     #pragma omp parallel for private (i0, i1)
499     for (i0=0; i0<chunkEle; i0++) {
500     mesh_p->Elements->Id[i0] = tempInts[i0*(2+numNodes)+0];
501     mesh_p->Elements->Tag[i0] = tempInts[i0*(2+numNodes)+1];
502     for (i1 = 0; i1 < numNodes; i1++) {
503     mesh_p->Elements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
504     }
505     }
506     }
507    
508     TMPMEMFREE(tempInts);
509     }
510    
511    
512     #if 0 /* this is the original code for reading elements */
513 ksteube 1360 /* read elements */
514     if (Finley_noError()) {
515 ksteube 1402
516 ksteube 1360 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
517     typeID=Finley_RefElement_getTypeId(element_type);
518     if (typeID==NoType) {
519     sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s",element_type);
520     Finley_setError(VALUE_ERROR,error_msg);
521     } else {
522     /* read the elements */
523     mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
524     if (Finley_noError()) {
525     Finley_ElementFile_allocTable(mesh_p->Elements, numEle);
526     mesh_p->Elements->minColor=0;
527     mesh_p->Elements->maxColor=numEle-1;
528     if (Finley_noError()) {
529     for (i0 = 0; i0 < numEle; i0++) {
530     fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);
531     mesh_p->Elements->Color[i0]=i0;
532     for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {
533     fscanf(fileHandle_p, " %d",
534     &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);
535     } /* for i1 */
536     fscanf(fileHandle_p, "\n");
537     } /* for i0 */
538     }
539     }
540     }
541     }
542 ksteube 1402 #endif
543    
544     printf("ksteube CPU=%d/%d Element typeID=%d\n", mpi_info->rank, mpi_info->size, typeID);
545    
546     #if 1
547    
548     /* Define other structures to keep mesh_write from crashing */
549    
550     mesh_p->FaceElements=Finley_ElementFile_alloc(0, mesh_p->order, mesh_p->reduced_order, mpi_info);
551     Finley_ElementFile_allocTable(mesh_p->FaceElements, 0);
552    
553     mesh_p->ContactElements=Finley_ElementFile_alloc(0, mesh_p->order, mesh_p->reduced_order, mpi_info);
554     Finley_ElementFile_allocTable(mesh_p->ContactElements, 0);
555    
556     mesh_p->Points=Finley_ElementFile_alloc(0, mesh_p->order, mesh_p->reduced_order, mpi_info);
557     Finley_ElementFile_allocTable(mesh_p->Points, 0);
558    
559     #endif
560    
561     #if 0 /* comment out the rest of the un-implemented crap for now */
562 ksteube 1360 /* get the face elements */
563     if (Finley_noError()) {
564     fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
565     faceTypeID=Finley_RefElement_getTypeId(element_type);
566     if (faceTypeID==NoType) {
567     sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s for face elements",element_type);
568     Finley_setError(VALUE_ERROR,error_msg);
569     } else {
570     mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
571     if (Finley_noError()) {
572     Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);
573     if (Finley_noError()) {
574     mesh_p->FaceElements->minColor=0;
575     mesh_p->FaceElements->maxColor=numEle-1;
576     for (i0 = 0; i0 < numEle; i0++) {
577     fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);
578     mesh_p->FaceElements->Color[i0]=i0;
579     for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {
580     fscanf(fileHandle_p, " %d",
581     &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);
582     } /* for i1 */
583     fscanf(fileHandle_p, "\n");
584     } /* for i0 */
585     }
586     }
587     }
588     }
589     /* get the Contact face element */
590     if (Finley_noError()) {
591     fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
592     contactTypeID=Finley_RefElement_getTypeId(element_type);
593     if (contactTypeID==NoType) {
594     sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);
595     Finley_setError(VALUE_ERROR,error_msg);
596     } else {
597     mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
598     if (Finley_noError()) {
599     Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);
600     if (Finley_noError()) {
601     mesh_p->ContactElements->minColor=0;
602     mesh_p->ContactElements->maxColor=numEle-1;
603     for (i0 = 0; i0 < numEle; i0++) {
604     fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);
605     mesh_p->ContactElements->Color[i0]=i0;
606     for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {
607     fscanf(fileHandle_p, " %d",
608     &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);
609     } /* for i1 */
610     fscanf(fileHandle_p, "\n");
611     } /* for i0 */
612     }
613     }
614     }
615 ksteube 1402 }
616 ksteube 1360 /* get the nodal element */
617     if (Finley_noError()) {
618     fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
619     pointTypeID=Finley_RefElement_getTypeId(element_type);
620     if (pointTypeID==NoType) {
621     sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type);
622     Finley_setError(VALUE_ERROR,error_msg);
623     }
624     mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
625     if (Finley_noError()) {
626     Finley_ElementFile_allocTable(mesh_p->Points, numEle);
627     if (Finley_noError()) {
628     mesh_p->Points->minColor=0;
629     mesh_p->Points->maxColor=numEle-1;
630     for (i0 = 0; i0 < numEle; i0++) {
631     fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);
632     mesh_p->Points->Color[i0]=i0;
633     for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {
634     fscanf(fileHandle_p, " %d",
635     &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);
636     } /* for i1 */
637     fscanf(fileHandle_p, "\n");
638     } /* for i0 */
639     }
640     }
641     }
642     /* get the name tags */
643     if (Finley_noError()) {
644     if (feof(fileHandle_p) == 0) {
645     fscanf(fileHandle_p, "%s\n", name);
646     while (feof(fileHandle_p) == 0) {
647     fscanf(fileHandle_p, "%s %d\n", name, &tag_key);
648     Finley_Mesh_addTagMap(mesh_p,name,tag_key);
649     }
650     }
651     }
652 ksteube 1402 #endif /* comment out the rest of the un-implemented crap for now */
653 ksteube 1360 }
654     /* close file */
655 ksteube 1402 if (mpi_info->rank == 0) fclose(fileHandle_p);
656    
657 ksteube 1360 /* resolve id's : */
658     /* rearrange elements: */
659 ksteube 1402
660 ksteube 1360 if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
661     if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
662 ksteube 1402 return mesh_p; /* ksteube temp return for debugging */
663    
664 ksteube 1360 /* that's it */
665     #ifdef Finley_TRACE
666     printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
667     #endif
668    
669     /* that's it */
670     if (! Finley_noError()) {
671     Finley_Mesh_free(mesh_p);
672     }
673     Paso_MPIInfo_free( mpi_info );
674     return mesh_p;
675     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26