/[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 1732 - (show annotations)
Wed Aug 27 05:47:03 2008 UTC (11 years, 1 month 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
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, faceTypeID, contactTypeID, pointTypeID;
245 index_t tag_key;
246
247 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
259 /* read header */
260 sprintf(frm,"%%%d[^\n]",LenString_MAX-1);
261 fscanf(fileHandle_p, frm, name);
262
263 /* get the number of nodes */
264 fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);
265 }
266
267 #ifdef PASO_MPI
268 /* MPI Broadcast numDim, numNodes, name if there are multiple MPI procs*/
269 if (mpi_info->size > 1) {
270 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 /* Each CPU will get at most chunkSize nodes so the message has to be sufficiently large */
293 int chunkSize = numNodes / mpi_info->size + 1, totalNodes=0, chunkNodes=0, chunkEle=0, nextCPU=1;
294 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
297 /*
298 Read chunkSize nodes, send it in a chunk to worker CPU which copies chunk into its local mesh_p
299 It doesn't matter that a CPU has the wrong nodes for its elements, this is sorted out later
300 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 for (i0=0; i0<chunkSize*3+1; i0++) tempInts[i0] = -1;
308 for (i0=0; i0<chunkSize*numDim; i0++) tempCoords[i0] = -1.0;
309 chunkNodes = 0;
310 for (i1=0; i1<chunkSize; i1++) {
311 if (totalNodes >= numNodes) break; /* End of inner loop */
312 if (1 == numDim)
313 fscanf(fileHandle_p, "%d %d %d %le\n",
314 &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
315 &tempCoords[i1*numDim+0]);
316 if (2 == numDim)
317 fscanf(fileHandle_p, "%d %d %d %le %le\n",
318 &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
319 &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1]);
320 if (3 == numDim)
321 fscanf(fileHandle_p, "%d %d %d %le %le %le\n",
322 &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
323 &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1], &tempCoords[i1*numDim+2]);
324 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 }
327 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 if (nextCPU < mpi_info->size) {
334 int mpi_error;
335 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 if ( mpi_error != MPI_SUCCESS ) {
338 Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts failed");
339 return NULL;
340 }
341 mpi_error = MPI_Send(tempCoords, chunkSize*numDim, MPI_DOUBLE, nextCPU, 81721, mpi_info->comm);
342 if ( mpi_error != MPI_SUCCESS ) {
343 Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempCoords failed");
344 return NULL;
345 }
346 }
347 #endif
348 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 } /* 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 /* 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 mesh_p->Nodes->globalDegreesOfFreedom[i0] = tempInts[chunkSize+i0];
378 mesh_p->Nodes->Tag[i0] = tempInts[chunkSize*2+i0];
379 for (i1=0; i1<numDim; i1++) {
380 mesh_p->Nodes->Coordinates[INDEX2(i1,i0,numDim)] = tempCoords[i0*numDim+i1];
381 }
382 }
383 }
384
385 TMPMEMFREE(tempInts);
386 TMPMEMFREE(tempCoords);
387
388 /* read elements */
389
390 /* 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 if (mpi_info->size > 1) {
398 int temp1[2], mpi_error;
399 temp1[0] = (int) typeID;
400 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 typeID = (ElementTypeId) temp1[0];
407 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 /* Allocate the ElementFile */
417 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 /* Read the element data */
421 if (Finley_noError()) {
422 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 if (mpi_info->rank == 0) { /* Master */
426 for (;;) { /* Infinite loop */
427 for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
428 chunkEle = 0;
429 for (i0=0; i0<chunkSize; i0++) {
430 if (totalEle >= numEle) break; /* End inner loop */
431 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 #ifdef PASO_MPI
438 /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
439 if (nextCPU < mpi_info->size) {
440 int mpi_error;
441 tempInts[chunkSize*(2+numNodes)] = chunkEle;
442 mpi_error = MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81722, mpi_info->comm);
443 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 #endif
449 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 } /* Infinite loop */
453 } /* End master */
454 else { /* Worker */
455 #ifdef PASO_MPI
456 /* Each worker receives one message */
457 MPI_Status status;
458 int mpi_error;
459 mpi_error = MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81722, mpi_info->comm, &status);
460 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 chunkEle = tempInts[chunkSize*(2+numNodes)];
465 #endif
466 } /* Worker */
467
468 /* 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 mesh_p->Elements->Color[i0] = i0;
478 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 } /* 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 }
787
788 /* close file */
789 if (mpi_info->rank == 0) fclose(fileHandle_p);
790
791 /* resolve id's : */
792 /* rearrange elements: */
793
794 /* return mesh_p; */ /* ksteube temp return for debugging */
795
796 if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
797 if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
798
799 /* 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