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 |
} |