/[escript]/branches/arrayview_from_1695_trunk/finley/src/Mesh_read.c
ViewVC logotype

Contents of /branches/arrayview_from_1695_trunk/finley/src/Mesh_read.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1781 - (show annotations)
Thu Sep 11 05:03:14 2008 UTC (10 years, 9 months ago) by jfenwick
File MIME type: text/plain
File size: 36015 byte(s)
Branch commit

Merged changes from trunk version 1695 up to and including version 1779.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26