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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1713 - (show annotations)
Wed Aug 20 07:03:45 2008 UTC (11 years, 2 months ago) by ksteube
File MIME type: text/plain
File size: 29540 byte(s)
MPI read of mesh.fly: reduced memory requirements

1
2 /* $Id$ */
3
4 /*******************************************************
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
16 /**************************************************************/
17
18 /* Finley: read mesh */
19
20 /**************************************************************/
21
22 #include "Mesh.h"
23
24 /**************************************************************/
25
26 /* reads a mesh from a Finley file of name fname */
27
28 Finley_Mesh* Finley_Mesh_read(char* fname,index_t order, index_t reduced_order, bool_t optimize)
29
30 {
31
32 Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
33 dim_t numNodes, numDim, numEle, i0, i1;
34 index_t tag_key;
35 Finley_Mesh *mesh_p=NULL;
36 char name[LenString_MAX],element_type[LenString_MAX],frm[20];
37 char error_msg[LenErrorMsg_MAX];
38 double time0=Finley_timer();
39 FILE *fileHandle_p = NULL;
40 ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;
41
42 Finley_resetError();
43
44 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 }
55
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 }
224
225 /* that's it */
226 if (! Finley_noError()) {
227 Finley_Mesh_free(mesh_p);
228 }
229 Paso_MPIInfo_free( mpi_info );
230 return mesh_p;
231 }
232
233 Finley_Mesh* Finley_Mesh_read_MPI(char* fname,index_t order, index_t reduced_order, bool_t optimize)
234
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 ElementTypeId typeID;
245 ElementTypeId faceTypeID, contactTypeID, pointTypeID;
246 index_t tag_key;
247
248 Finley_resetError();
249
250 if (mpi_info->rank == 0) {
251 /* get file handle */
252 fileHandle_p = fopen(fname, "r");
253 if (fileHandle_p==NULL) {
254 sprintf(error_msg,"Finley_Mesh_read: Opening file %s for reading failed.",fname);
255 Finley_setError(IO_ERROR,error_msg);
256 Paso_MPIInfo_free( mpi_info );
257 return NULL;
258 }
259
260 /* read header */
261 sprintf(frm,"%%%d[^\n]",LenString_MAX-1);
262 fscanf(fileHandle_p, frm, name);
263
264 /* get the number of nodes */
265 fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);
266 }
267
268 #ifdef PASO_MPI
269 /* MPI Broadcast numDim, numNodes, name if there are multiple MPI procs*/
270 if (mpi_info->size > 1) {
271 int temp1[3], error_code;
272 temp1[0] = numDim;
273 temp1[1] = numNodes;
274 temp1[2] = strlen(name) + 1;
275 error_code = MPI_Bcast (temp1, 3, MPI_INT, 0, mpi_info->comm);
276 if (error_code != MPI_SUCCESS) {
277 Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of temp1 failed");
278 return NULL;
279 }
280 numDim = temp1[0];
281 numNodes = temp1[1];
282 error_code = MPI_Bcast (name, temp1[2], MPI_CHAR, 0, mpi_info->comm);
283 if (error_code != MPI_SUCCESS) {
284 Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of name failed");
285 return NULL;
286 }
287 }
288 #endif
289
290 /* allocate mesh */
291 mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info);
292 if (Finley_noError()) {
293 /* Each CPU will get at most chunkSize nodes so the message has to be sufficiently large for that */
294 int chunkSize = numNodes / mpi_info->size + 1, totalNodes=0, chunkNodes=0, chunkEle=0, nextCPU=1;
295 int *tempInts = TMPMEMALLOC(chunkSize*3+1, index_t); /* Stores the integer message data */
296 double *tempCoords = TMPMEMALLOC(chunkSize*numDim, double); /* Stores the double message data */
297
298 /*
299 Read chunkSize nodes, send it in a chunk to worker CPU which copies chunk into its local mesh_p
300 It doesn't matter that a CPU has the wrong nodes for its elements, this is sorted out later
301 First chunk sent to CPU 1, second to CPU 2, ...
302 Last chunk stays on CPU 0 (the master)
303 The three columns of integers (Id, gDOF, Tag) are gathered into a single array tempInts and sent together in a single MPI message
304 */
305
306 if (mpi_info->rank == 0) { /* Master */
307 for (;;) { /* Infinite loop */
308 for (i0=0; i0<chunkSize*3+1; i0++) tempInts[i0] = -1;
309 for (i0=0; i0<chunkSize*numDim; i0++) tempCoords[i0] = -1.0;
310 chunkNodes = 0;
311 for (i1=0; i1<chunkSize; i1++) {
312 if (totalNodes >= numNodes) break; /* Maybe end the infinite loop */
313 if (1 == numDim)
314 fscanf(fileHandle_p, "%d %d %d %le\n",
315 &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
316 &tempCoords[i1*numDim+0]);
317 if (2 == numDim)
318 fscanf(fileHandle_p, "%d %d %d %le %le\n",
319 &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
320 &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1]);
321 if (3 == numDim)
322 fscanf(fileHandle_p, "%d %d %d %le %le %le\n",
323 &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
324 &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1], &tempCoords[i1*numDim+2]);
325 totalNodes++; /* When do we quit the infinite loop? */
326 chunkNodes++; /* How many nodes do we actually have in this chunk? It may be smaller than chunkSize. */
327 }
328 if (chunkNodes > chunkSize) {
329 Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: error reading chunks of mesh, data too large for message size");
330 return NULL;
331 }
332 #ifdef PASO_MPI
333 /* Eventually we'll send chunkSize nodes to each CPU numbered 1 ... mpi_info->size-1, here goes one of them */
334 if (nextCPU < mpi_info->size) {
335 int mpi_error;
336 tempInts[chunkSize*3] = chunkNodes; /* The message has one more int to send chunkNodes */
337 mpi_error = MPI_Send(tempInts, chunkSize*3+1, MPI_INT, nextCPU, 81720, mpi_info->comm);
338 if ( mpi_error != MPI_SUCCESS ) {
339 Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts failed");
340 return NULL;
341 }
342 mpi_error = MPI_Send(tempCoords, chunkSize*numDim, MPI_DOUBLE, nextCPU, 81721, mpi_info->comm);
343 if ( mpi_error != MPI_SUCCESS ) {
344 Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempCoords failed");
345 return NULL;
346 }
347 nextCPU++;
348 }
349 #endif
350 if (totalNodes >= numNodes) break; /* Maybe end the infinite loop */
351 } /* Infinite loop */
352 } /* End master */
353 else { /* Worker */
354 #ifdef PASO_MPI
355 /* Each worker receives two messages */
356 MPI_Status status;
357 int mpi_error;
358 mpi_error = MPI_Recv(tempInts, chunkSize*3+1, MPI_INT, 0, 81720, mpi_info->comm, &status);
359 if ( mpi_error != MPI_SUCCESS ) {
360 Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts failed");
361 return NULL;
362 }
363 mpi_error = MPI_Recv(tempCoords, chunkSize*numDim, MPI_DOUBLE, 0, 81721, mpi_info->comm, &status);
364 if ( mpi_error != MPI_SUCCESS ) {
365 Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempCoords failed");
366 return NULL;
367 }
368 chunkNodes = tempInts[chunkSize*3]; /* How many nodes are in this workers chunk? */
369 #endif
370 } /* Worker */
371
372 #if 1
373 /* Display the temp mem for debugging */
374 printf("ksteube Nodes tempInts\n");
375 for (i0=0; i0<chunkSize*3; i0++) {
376 printf(" %2d", tempInts[i0]);
377 if (i0%chunkSize==chunkSize-1) printf("\n");
378 }
379 printf("ksteube tempCoords:\n");
380 for (i0=0; i0<chunkSize*numDim; i0++) {
381 printf(" %20.15e", tempCoords[i0]);
382 if (i0%numDim==numDim-1) printf("\n");
383 }
384 #endif
385
386 printf("ksteube numDim=%d numNodes=%d mesh name='%s' chunkNodes=%d numNodes=%d\n", numDim, numNodes, name, chunkNodes, numNodes);
387
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[chunkSize+i0];
394 mesh_p->Nodes->Tag[i0] = tempInts[chunkSize*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 }
400
401 TMPMEMFREE(tempInts);
402 TMPMEMFREE(tempCoords);
403
404 /* read elements */
405
406 /* 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 > 1) {
414 int temp1[3], mpi_error;
415 temp1[0] = (int) 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 = (ElementTypeId) 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 int chunkSize = numEle / mpi_info->size, totalEle=0, nextCPU=1;
439 if (numEle % mpi_info->size != 0) chunkSize++; /* Remainder from numEle / mpi_info->size will be spread out one-per-CPU */
440 if (mpi_info->rank == 0) { /* Master */
441 for (;;) { /* Infinite loop */
442 chunkEle = 0;
443 for (i0=0; i0<numEle*(2+numNodes)+1; i0++) tempInts[i0] = -1;
444 for (i0=0; i0<chunkSize; i0++) {
445 if (totalEle >= numEle) break; /* End infinite loop */
446 fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
447 for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
448 fscanf(fileHandle_p, "\n");
449 totalEle++;
450 chunkEle++;
451 }
452 /* Eventually we'll send chunk of nodes to each CPU except 0 itself, here goes one of them */
453 if (nextCPU < mpi_info->size) {
454 #ifdef PASO_MPI
455 int mpi_error;
456
457 tempInts[numEle*(2+numNodes)] = chunkEle;
458 printf("ksteube CPU=%d/%d send to %d\n", mpi_info->rank, mpi_info->size, nextCPU);
459 mpi_error = MPI_Send(tempInts, numEle*(2+numNodes)+1, MPI_INT, nextCPU, 81722, mpi_info->comm);
460 if ( mpi_error != MPI_SUCCESS ) {
461 Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Elements failed");
462 return NULL;
463 }
464 #endif
465 nextCPU++;
466 }
467 if (totalEle >= numEle) break; /* End infinite loop */
468 } /* Infinite loop */
469 } /* End master */
470 else { /* Worker */
471 #ifdef PASO_MPI
472 /* Each worker receives two messages */
473 MPI_Status status;
474 int mpi_error;
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 0
485 /* Display the temp mem for debugging */
486 printf("ksteube Elements 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 /* read elements */
514 if (Finley_noError()) {
515
516 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 #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 /* Change the typeid from NoType later */
550
551 mesh_p->FaceElements=Finley_ElementFile_alloc(NoType, mesh_p->order, mesh_p->reduced_order, mpi_info);
552 Finley_ElementFile_allocTable(mesh_p->FaceElements, 0);
553
554 mesh_p->ContactElements=Finley_ElementFile_alloc(NoType, mesh_p->order, mesh_p->reduced_order, mpi_info);
555 Finley_ElementFile_allocTable(mesh_p->ContactElements, 0);
556
557 mesh_p->Points=Finley_ElementFile_alloc(NoType, mesh_p->order, mesh_p->reduced_order, mpi_info);
558 Finley_ElementFile_allocTable(mesh_p->Points, 0);
559
560 #endif
561
562 #if 0 /* comment out the rest of the un-implemented crap for now */
563 /* get the face elements */
564 if (Finley_noError()) {
565 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
566 faceTypeID=Finley_RefElement_getTypeId(element_type);
567 if (faceTypeID==NoType) {
568 sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s for face elements",element_type);
569 Finley_setError(VALUE_ERROR,error_msg);
570 } else {
571 mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
572 if (Finley_noError()) {
573 Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);
574 if (Finley_noError()) {
575 mesh_p->FaceElements->minColor=0;
576 mesh_p->FaceElements->maxColor=numEle-1;
577 for (i0 = 0; i0 < numEle; i0++) {
578 fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);
579 mesh_p->FaceElements->Color[i0]=i0;
580 for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {
581 fscanf(fileHandle_p, " %d",
582 &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);
583 } /* for i1 */
584 fscanf(fileHandle_p, "\n");
585 } /* for i0 */
586 }
587 }
588 }
589 }
590 /* get the Contact face element */
591 if (Finley_noError()) {
592 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
593 contactTypeID=Finley_RefElement_getTypeId(element_type);
594 if (contactTypeID==NoType) {
595 sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);
596 Finley_setError(VALUE_ERROR,error_msg);
597 } else {
598 mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
599 if (Finley_noError()) {
600 Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);
601 if (Finley_noError()) {
602 mesh_p->ContactElements->minColor=0;
603 mesh_p->ContactElements->maxColor=numEle-1;
604 for (i0 = 0; i0 < numEle; i0++) {
605 fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);
606 mesh_p->ContactElements->Color[i0]=i0;
607 for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {
608 fscanf(fileHandle_p, " %d",
609 &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);
610 } /* for i1 */
611 fscanf(fileHandle_p, "\n");
612 } /* for i0 */
613 }
614 }
615 }
616 }
617 /* get the nodal element */
618 if (Finley_noError()) {
619 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
620 pointTypeID=Finley_RefElement_getTypeId(element_type);
621 if (pointTypeID==NoType) {
622 sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type);
623 Finley_setError(VALUE_ERROR,error_msg);
624 }
625 mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
626 if (Finley_noError()) {
627 Finley_ElementFile_allocTable(mesh_p->Points, numEle);
628 if (Finley_noError()) {
629 mesh_p->Points->minColor=0;
630 mesh_p->Points->maxColor=numEle-1;
631 for (i0 = 0; i0 < numEle; i0++) {
632 fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);
633 mesh_p->Points->Color[i0]=i0;
634 for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {
635 fscanf(fileHandle_p, " %d",
636 &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);
637 } /* for i1 */
638 fscanf(fileHandle_p, "\n");
639 } /* for i0 */
640 }
641 }
642 }
643 /* get the name tags */
644 if (Finley_noError()) {
645 if (feof(fileHandle_p) == 0) {
646 fscanf(fileHandle_p, "%s\n", name);
647 while (feof(fileHandle_p) == 0) {
648 fscanf(fileHandle_p, "%s %d\n", name, &tag_key);
649 Finley_Mesh_addTagMap(mesh_p,name,tag_key);
650 }
651 }
652 }
653 #endif /* comment out the rest of the un-implemented crap for now */
654 }
655 /* close file */
656 if (mpi_info->rank == 0) fclose(fileHandle_p);
657
658 /* resolve id's : */
659 /* rearrange elements: */
660
661 if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
662 if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
663 return mesh_p; /* ksteube temp return for debugging */
664
665 /* that's it */
666 #ifdef Finley_TRACE
667 printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
668 #endif
669
670 /* that's it */
671 if (! Finley_noError()) {
672 Finley_Mesh_free(mesh_p);
673 }
674 Paso_MPIInfo_free( mpi_info );
675 return mesh_p;
676 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26