/[escript]/branches/doubleplusgood/finley/src/Mesh_read.cpp
ViewVC logotype

Contents of /branches/doubleplusgood/finley/src/Mesh_read.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1811 - (show annotations)
Thu Sep 25 23:11:13 2008 UTC (11 years, 1 month ago) by ksteube
Original Path: trunk/finley/src/Mesh_read.c
File MIME type: text/plain
File size: 35980 byte(s)
Copyright updated in all files

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26