1 |
|
2 |
/* $Id$ */ |
3 |
|
4 |
/******************************************************* |
5 |
* |
6 |
* Copyright 2003-2007 by ACceSS MNRF |
7 |
* Copyright 2007 by University of Queensland |
8 |
* |
9 |
* http://esscc.uq.edu.au |
10 |
* Primary Business: Queensland, Australia |
11 |
* Licensed under the Open Software License version 3.0 |
12 |
* http://www.opensource.org/licenses/osl-3.0.php |
13 |
* |
14 |
*******************************************************/ |
15 |
|
16 |
/**************************************************************/ |
17 |
|
18 |
/* Finley: read mesh */ |
19 |
|
20 |
/**************************************************************/ |
21 |
|
22 |
#include "Mesh.h" |
23 |
|
24 |
/**************************************************************/ |
25 |
|
26 |
/* reads a mesh from a Finley file of name fname */ |
27 |
|
28 |
Finley_Mesh* Finley_Mesh_read(char* fname,index_t order, index_t reduced_order, bool_t optimize) |
29 |
|
30 |
{ |
31 |
|
32 |
Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD ); |
33 |
dim_t numNodes, numDim, numEle, i0, i1; |
34 |
index_t tag_key; |
35 |
Finley_Mesh *mesh_p=NULL; |
36 |
char name[LenString_MAX],element_type[LenString_MAX],frm[20]; |
37 |
char error_msg[LenErrorMsg_MAX]; |
38 |
double time0=Finley_timer(); |
39 |
FILE *fileHandle_p = NULL; |
40 |
ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID; |
41 |
|
42 |
Finley_resetError(); |
43 |
|
44 |
if (mpi_info->size > 1) { |
45 |
Finley_setError(SYSTEM_ERROR,"Finley_Mesh_read: MPI is not suporrted yet."); |
46 |
} else { |
47 |
/* get file handle */ |
48 |
fileHandle_p = fopen(fname, "r"); |
49 |
if (fileHandle_p==NULL) { |
50 |
sprintf(error_msg,"Finley_Mesh_read: Opening file %s for reading failed.",fname); |
51 |
Finley_setError(IO_ERROR,error_msg); |
52 |
Paso_MPIInfo_free( mpi_info ); |
53 |
return NULL; |
54 |
} |
55 |
|
56 |
/* read header */ |
57 |
sprintf(frm,"%%%d[^\n]",LenString_MAX-1); |
58 |
fscanf(fileHandle_p, frm, name); |
59 |
|
60 |
/* get the nodes */ |
61 |
|
62 |
fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes); |
63 |
/* allocate mesh */ |
64 |
mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order,mpi_info); |
65 |
if (Finley_noError()) { |
66 |
|
67 |
/* read nodes */ |
68 |
Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes); |
69 |
if (Finley_noError()) { |
70 |
if (1 == numDim) { |
71 |
for (i0 = 0; i0 < numNodes; i0++) |
72 |
fscanf(fileHandle_p, "%d %d %d %le\n", &mesh_p->Nodes->Id[i0], |
73 |
&mesh_p->Nodes->globalDegreesOfFreedom[i0], &mesh_p->Nodes->Tag[i0], |
74 |
&mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)]); |
75 |
} else if (2 == numDim) { |
76 |
for (i0 = 0; i0 < numNodes; i0++) |
77 |
fscanf(fileHandle_p, "%d %d %d %le %le\n", &mesh_p->Nodes->Id[i0], |
78 |
&mesh_p->Nodes->globalDegreesOfFreedom[i0], &mesh_p->Nodes->Tag[i0], |
79 |
&mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)], |
80 |
&mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)]); |
81 |
} else if (3 == numDim) { |
82 |
for (i0 = 0; i0 < numNodes; i0++) |
83 |
fscanf(fileHandle_p, "%d %d %d %le %le %le\n", &mesh_p->Nodes->Id[i0], |
84 |
&mesh_p->Nodes->globalDegreesOfFreedom[i0], &mesh_p->Nodes->Tag[i0], |
85 |
&mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)], |
86 |
&mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)], |
87 |
&mesh_p->Nodes->Coordinates[INDEX2(2,i0,numDim)]); |
88 |
} /* if else else */ |
89 |
} |
90 |
/* read elements */ |
91 |
if (Finley_noError()) { |
92 |
|
93 |
fscanf(fileHandle_p, "%s %d\n", element_type, &numEle); |
94 |
typeID=Finley_RefElement_getTypeId(element_type); |
95 |
if (typeID==NoType) { |
96 |
sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s",element_type); |
97 |
Finley_setError(VALUE_ERROR,error_msg); |
98 |
} else { |
99 |
/* read the elements */ |
100 |
mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info); |
101 |
if (Finley_noError()) { |
102 |
Finley_ElementFile_allocTable(mesh_p->Elements, numEle); |
103 |
mesh_p->Elements->minColor=0; |
104 |
mesh_p->Elements->maxColor=numEle-1; |
105 |
if (Finley_noError()) { |
106 |
for (i0 = 0; i0 < numEle; i0++) { |
107 |
fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]); |
108 |
mesh_p->Elements->Color[i0]=i0; |
109 |
for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) { |
110 |
fscanf(fileHandle_p, " %d", |
111 |
&mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]); |
112 |
} /* for i1 */ |
113 |
fscanf(fileHandle_p, "\n"); |
114 |
} /* for i0 */ |
115 |
} |
116 |
} |
117 |
} |
118 |
} |
119 |
/* get the face elements */ |
120 |
if (Finley_noError()) { |
121 |
fscanf(fileHandle_p, "%s %d\n", element_type, &numEle); |
122 |
faceTypeID=Finley_RefElement_getTypeId(element_type); |
123 |
if (faceTypeID==NoType) { |
124 |
sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s for face elements",element_type); |
125 |
Finley_setError(VALUE_ERROR,error_msg); |
126 |
} else { |
127 |
mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info); |
128 |
if (Finley_noError()) { |
129 |
Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle); |
130 |
if (Finley_noError()) { |
131 |
mesh_p->FaceElements->minColor=0; |
132 |
mesh_p->FaceElements->maxColor=numEle-1; |
133 |
for (i0 = 0; i0 < numEle; i0++) { |
134 |
fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]); |
135 |
mesh_p->FaceElements->Color[i0]=i0; |
136 |
for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) { |
137 |
fscanf(fileHandle_p, " %d", |
138 |
&mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]); |
139 |
} /* for i1 */ |
140 |
fscanf(fileHandle_p, "\n"); |
141 |
} /* for i0 */ |
142 |
} |
143 |
} |
144 |
} |
145 |
} |
146 |
/* get the Contact face element */ |
147 |
if (Finley_noError()) { |
148 |
fscanf(fileHandle_p, "%s %d\n", element_type, &numEle); |
149 |
contactTypeID=Finley_RefElement_getTypeId(element_type); |
150 |
if (contactTypeID==NoType) { |
151 |
sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type); |
152 |
Finley_setError(VALUE_ERROR,error_msg); |
153 |
} else { |
154 |
mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info); |
155 |
if (Finley_noError()) { |
156 |
Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle); |
157 |
if (Finley_noError()) { |
158 |
mesh_p->ContactElements->minColor=0; |
159 |
mesh_p->ContactElements->maxColor=numEle-1; |
160 |
for (i0 = 0; i0 < numEle; i0++) { |
161 |
fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]); |
162 |
mesh_p->ContactElements->Color[i0]=i0; |
163 |
for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) { |
164 |
fscanf(fileHandle_p, " %d", |
165 |
&mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]); |
166 |
} /* for i1 */ |
167 |
fscanf(fileHandle_p, "\n"); |
168 |
} /* for i0 */ |
169 |
} |
170 |
} |
171 |
} |
172 |
} |
173 |
/* get the nodal element */ |
174 |
if (Finley_noError()) { |
175 |
fscanf(fileHandle_p, "%s %d\n", element_type, &numEle); |
176 |
pointTypeID=Finley_RefElement_getTypeId(element_type); |
177 |
if (pointTypeID==NoType) { |
178 |
sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type); |
179 |
Finley_setError(VALUE_ERROR,error_msg); |
180 |
} |
181 |
mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info); |
182 |
if (Finley_noError()) { |
183 |
Finley_ElementFile_allocTable(mesh_p->Points, numEle); |
184 |
if (Finley_noError()) { |
185 |
mesh_p->Points->minColor=0; |
186 |
mesh_p->Points->maxColor=numEle-1; |
187 |
for (i0 = 0; i0 < numEle; i0++) { |
188 |
fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]); |
189 |
mesh_p->Points->Color[i0]=i0; |
190 |
for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) { |
191 |
fscanf(fileHandle_p, " %d", |
192 |
&mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]); |
193 |
} /* for i1 */ |
194 |
fscanf(fileHandle_p, "\n"); |
195 |
} /* for i0 */ |
196 |
} |
197 |
} |
198 |
} |
199 |
/* get the name tags */ |
200 |
if (Finley_noError()) { |
201 |
if (feof(fileHandle_p) == 0) { |
202 |
fscanf(fileHandle_p, "%s\n", name); |
203 |
while (feof(fileHandle_p) == 0) { |
204 |
fscanf(fileHandle_p, "%s %d\n", name, &tag_key); |
205 |
Finley_Mesh_addTagMap(mesh_p,name,tag_key); |
206 |
} |
207 |
} |
208 |
} |
209 |
} |
210 |
/* close file */ |
211 |
fclose(fileHandle_p); |
212 |
|
213 |
/* resolve id's : */ |
214 |
/* rearrange elements: */ |
215 |
|
216 |
if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p); |
217 |
if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize); |
218 |
|
219 |
/* that's it */ |
220 |
#ifdef Finley_TRACE |
221 |
printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0); |
222 |
#endif |
223 |
} |
224 |
|
225 |
/* that's it */ |
226 |
if (! Finley_noError()) { |
227 |
Finley_Mesh_free(mesh_p); |
228 |
} |
229 |
Paso_MPIInfo_free( mpi_info ); |
230 |
return mesh_p; |
231 |
} |
232 |
|
233 |
Finley_Mesh* Finley_Mesh_read_MPI(char* fname,index_t order, index_t reduced_order, bool_t optimize) |
234 |
|
235 |
{ |
236 |
|
237 |
Paso_MPIInfo *mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD ); |
238 |
dim_t numNodes, numDim, numEle, i0, i1; |
239 |
Finley_Mesh *mesh_p=NULL; |
240 |
char name[LenString_MAX],element_type[LenString_MAX],frm[20]; |
241 |
char error_msg[LenErrorMsg_MAX]; |
242 |
double time0=Finley_timer(); |
243 |
FILE *fileHandle_p = NULL; |
244 |
ElementTypeId typeID; |
245 |
|
246 |
#if 0 /* comment out the rest of the un-implemented crap for now */ |
247 |
/* See below */ |
248 |
ElementTypeId faceTypeID, contactTypeID, pointTypeID; |
249 |
index_t tag_key; |
250 |
#endif |
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 */ |
274 |
if (mpi_info->size > 0) { |
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 |
int chunkSize = numNodes / mpi_info->size + 1, totalNodes=0, chunkNodes=0, chunkEle=0, nextCPU=1; |
298 |
int *tempInts = TMPMEMALLOC(numNodes*3+1, index_t); |
299 |
double *tempCoords = TMPMEMALLOC(numNodes*numDim, double); |
300 |
|
301 |
/* |
302 |
Read a chunk of nodes, send to worker CPU if available, copy chunk into local mesh_p |
303 |
It doesn't matter that a CPU has the wrong nodes for its elements, this is sorted out later |
304 |
First chunk sent to CPU 1, second to CPU 2, ... |
305 |
Last chunk stays on CPU 0 (the master) |
306 |
The three columns of integers (Id, gDOF, Tag) are gathered into a single array tempInts and sent together in a single MPI message |
307 |
*/ |
308 |
|
309 |
if (mpi_info->rank == 0) { /* Master */ |
310 |
for (;;) { /* Infinite loop */ |
311 |
chunkNodes = 0; |
312 |
for (i0=0; i0<numNodes*3+1; i0++) tempInts[i0] = -1; |
313 |
for (i0=0; i0<numNodes*numDim; i0++) tempCoords[i0] = -1.0; |
314 |
for (i1=0; i1<chunkSize; i1++) { |
315 |
if (totalNodes >= numNodes) break; |
316 |
if (1 == numDim) |
317 |
fscanf(fileHandle_p, "%d %d %d %le\n", |
318 |
&tempInts[0+i1], &tempInts[numNodes+i1], &tempInts[numNodes*2+i1], |
319 |
&tempCoords[i1*numDim+0]); |
320 |
if (2 == numDim) |
321 |
fscanf(fileHandle_p, "%d %d %d %le %le\n", |
322 |
&tempInts[0+i1], &tempInts[numNodes+i1], &tempInts[numNodes*2+i1], |
323 |
&tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1]); |
324 |
if (3 == numDim) |
325 |
fscanf(fileHandle_p, "%d %d %d %le %le %le\n", |
326 |
&tempInts[0+i1], &tempInts[numNodes+i1], &tempInts[numNodes*2+i1], |
327 |
&tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1], &tempCoords[i1*numDim+2]); |
328 |
totalNodes++; |
329 |
chunkNodes++; |
330 |
} |
331 |
/* Eventually we'll send chunk of nodes to each CPU numbered 1 ... mpi_info->size-1, here goes one of them */ |
332 |
if (nextCPU < mpi_info->size) { |
333 |
#ifdef PASO_MPI |
334 |
int mpi_error; |
335 |
|
336 |
tempInts[numNodes*3] = chunkNodes; |
337 |
/* ksteube The size of this message can and should be brought down to chunkNodes*3+1, must re-org tempInts */ |
338 |
mpi_error = MPI_Send(tempInts, numNodes*3+1, MPI_INT, nextCPU, 81720, mpi_info->comm); |
339 |
if ( mpi_error != MPI_SUCCESS ) { |
340 |
Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts failed"); |
341 |
return NULL; |
342 |
} |
343 |
mpi_error = MPI_Send(tempCoords, numNodes*numDim, MPI_DOUBLE, nextCPU, 81721, mpi_info->comm); |
344 |
if ( mpi_error != MPI_SUCCESS ) { |
345 |
Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempCoords failed"); |
346 |
return NULL; |
347 |
} |
348 |
#endif |
349 |
nextCPU++; |
350 |
} |
351 |
if (totalNodes >= numNodes) break; |
352 |
} /* Infinite loop */ |
353 |
} /* End master */ |
354 |
else { /* Worker */ |
355 |
#ifdef PASO_MPI |
356 |
/* Each worker receives two messages */ |
357 |
MPI_Status status; |
358 |
mpi_error = MPI_Recv(tempInts, numNodes*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, numNodes*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[numNodes*3]; |
369 |
#endif |
370 |
} /* Worker */ |
371 |
|
372 |
#if 0 |
373 |
/* Display the temp mem for debugging */ |
374 |
printf("ksteube tempInts totalNodes=%d:\n", totalNodes); |
375 |
for (i0=0; i0<numNodes*3; i0++) { |
376 |
printf(" %2d", tempInts[i0]); |
377 |
if (i0%numNodes==numNodes-1) printf("\n"); |
378 |
} |
379 |
printf("ksteube tempCoords:\n"); |
380 |
for (i0=0; i0<chunkNodes*numDim; i0++) { |
381 |
printf(" %20.15e", tempCoords[i0]); |
382 |
if (i0%numDim==numDim-1) printf("\n"); |
383 |
} |
384 |
#endif |
385 |
|
386 |
printf("ksteube chunkNodes=%d numNodes=%d\n", chunkNodes, numNodes); |
387 |
/* Copy node data from tempMem to mesh_p */ |
388 |
Finley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes); |
389 |
if (Finley_noError()) { |
390 |
for (i0=0; i0<chunkNodes; i0++) { |
391 |
mesh_p->Nodes->Id[i0] = tempInts[0+i0]; |
392 |
mesh_p->Nodes->globalDegreesOfFreedom[i0] = tempInts[numNodes+i0]; |
393 |
mesh_p->Nodes->Tag[i0] = tempInts[numNodes*2+i0]; |
394 |
for (i1=0; i1<numDim; i1++) { |
395 |
mesh_p->Nodes->Coordinates[INDEX2(i1,i0,numDim)] = tempCoords[i0*numDim+i1]; |
396 |
} |
397 |
} |
398 |
} |
399 |
|
400 |
TMPMEMFREE(tempInts); |
401 |
TMPMEMFREE(tempCoords); |
402 |
|
403 |
/* read elements */ |
404 |
|
405 |
/* Read the element typeID */ |
406 |
if (Finley_noError()) { |
407 |
if (mpi_info->rank == 0) { |
408 |
fscanf(fileHandle_p, "%s %d\n", element_type, &numEle); |
409 |
typeID=Finley_RefElement_getTypeId(element_type); |
410 |
} |
411 |
#ifdef PASO_MPI |
412 |
if (mpi_info->size > 0) { |
413 |
int temp1[3]; |
414 |
temp1[0] = typeID; |
415 |
temp1[1] = numEle; |
416 |
mpi_error = MPI_Bcast (temp1, 2, MPI_INT, 0, mpi_info->comm); |
417 |
if (mpi_error != MPI_SUCCESS) { |
418 |
Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed"); |
419 |
return NULL; |
420 |
} |
421 |
typeID = temp1[0]; |
422 |
numEle = temp1[1]; |
423 |
} |
424 |
#endif |
425 |
if (typeID==NoType) { |
426 |
sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type); |
427 |
Finley_setError(VALUE_ERROR, error_msg); |
428 |
} |
429 |
} |
430 |
|
431 |
/* Read the element data */ |
432 |
mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info); |
433 |
numNodes = mesh_p->Elements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */ |
434 |
|
435 |
if (Finley_noError()) { |
436 |
int *tempInts = TMPMEMALLOC(numEle*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */ |
437 |
int chunkSize = numEle / mpi_info->size, totalEle=0, nextCPU=1; |
438 |
if (numEle % mpi_info->size != 0) chunkSize++; /* Remainder from numEle / mpi_info->size will be spread out one-per-CPU */ |
439 |
if (mpi_info->rank == 0) { /* Master */ |
440 |
for (;;) { /* Infinite loop */ |
441 |
chunkEle = 0; |
442 |
for (i0=0; i0<numEle*(2+numNodes)+1; i0++) tempInts[i0] = -1; |
443 |
for (i0=0; i0<chunkSize; i0++) { |
444 |
if (totalEle >= numEle) break; /* End infinite loop */ |
445 |
fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]); |
446 |
for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]); |
447 |
fscanf(fileHandle_p, "\n"); |
448 |
totalEle++; |
449 |
chunkEle++; |
450 |
} |
451 |
/* Eventually we'll send chunk of nodes to each CPU except 0 itself, here goes one of them */ |
452 |
if (nextCPU < mpi_info->size) { |
453 |
#ifdef PASO_MPI |
454 |
int mpi_error; |
455 |
|
456 |
tempInts[numEle*(2+numNodes)] = chunkEle; |
457 |
printf("ksteube CPU=%d/%d send to %d\n", mpi_info->rank, mpi_info->size, nextCPU); |
458 |
mpi_error = MPI_Send(tempInts, numEle*(2+numNodes)+1, MPI_INT, nextCPU, 81722, mpi_info->comm); |
459 |
if ( mpi_error != MPI_SUCCESS ) { |
460 |
Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Elements failed"); |
461 |
return NULL; |
462 |
} |
463 |
#endif |
464 |
nextCPU++; |
465 |
} |
466 |
if (totalEle >= numEle) break; /* End infinite loop */ |
467 |
} /* Infinite loop */ |
468 |
} /* End master */ |
469 |
else { /* Worker */ |
470 |
#ifdef PASO_MPI |
471 |
/* Each worker receives two messages */ |
472 |
MPI_Status status; |
473 |
printf("ksteube CPU=%d/%d recv on %d\n", mpi_info->rank, mpi_info->size, mpi_info->rank); |
474 |
mpi_error = MPI_Recv(tempInts, numEle*(2+numNodes)+1, MPI_INT, 0, 81722, mpi_info->comm, &status); |
475 |
if ( mpi_error != MPI_SUCCESS ) { |
476 |
Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for Elements failed"); |
477 |
return NULL; |
478 |
} |
479 |
chunkEle = tempInts[numEle*(2+numNodes)]; |
480 |
#endif |
481 |
} /* Worker */ |
482 |
#if 1 |
483 |
/* Display the temp mem for debugging */ |
484 |
printf("ksteube tempInts numEle=%d chunkEle=%d AAA:\n", numEle, chunkEle); |
485 |
for (i0=0; i0<numEle*(numNodes+2); i0++) { |
486 |
printf(" %2d", tempInts[i0]); |
487 |
if (i0%(numNodes+2)==numNodes+2-1) printf("\n"); |
488 |
} |
489 |
#endif |
490 |
|
491 |
/* Copy Element data from tempInts to mesh_p */ |
492 |
Finley_ElementFile_allocTable(mesh_p->Elements, chunkEle); |
493 |
mesh_p->Elements->minColor=0; |
494 |
mesh_p->Elements->maxColor=chunkEle-1; |
495 |
if (Finley_noError()) { |
496 |
#pragma omp parallel for private (i0, i1) |
497 |
for (i0=0; i0<chunkEle; i0++) { |
498 |
mesh_p->Elements->Id[i0] = tempInts[i0*(2+numNodes)+0]; |
499 |
mesh_p->Elements->Tag[i0] = tempInts[i0*(2+numNodes)+1]; |
500 |
for (i1 = 0; i1 < numNodes; i1++) { |
501 |
mesh_p->Elements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1]; |
502 |
} |
503 |
} |
504 |
} |
505 |
|
506 |
TMPMEMFREE(tempInts); |
507 |
} |
508 |
|
509 |
|
510 |
#if 0 /* this is the original code for reading elements */ |
511 |
/* read elements */ |
512 |
if (Finley_noError()) { |
513 |
|
514 |
fscanf(fileHandle_p, "%s %d\n", element_type, &numEle); |
515 |
typeID=Finley_RefElement_getTypeId(element_type); |
516 |
if (typeID==NoType) { |
517 |
sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s",element_type); |
518 |
Finley_setError(VALUE_ERROR,error_msg); |
519 |
} else { |
520 |
/* read the elements */ |
521 |
mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info); |
522 |
if (Finley_noError()) { |
523 |
Finley_ElementFile_allocTable(mesh_p->Elements, numEle); |
524 |
mesh_p->Elements->minColor=0; |
525 |
mesh_p->Elements->maxColor=numEle-1; |
526 |
if (Finley_noError()) { |
527 |
for (i0 = 0; i0 < numEle; i0++) { |
528 |
fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]); |
529 |
mesh_p->Elements->Color[i0]=i0; |
530 |
for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) { |
531 |
fscanf(fileHandle_p, " %d", |
532 |
&mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]); |
533 |
} /* for i1 */ |
534 |
fscanf(fileHandle_p, "\n"); |
535 |
} /* for i0 */ |
536 |
} |
537 |
} |
538 |
} |
539 |
} |
540 |
#endif |
541 |
|
542 |
printf("ksteube CPU=%d/%d Element typeID=%d\n", mpi_info->rank, mpi_info->size, typeID); |
543 |
|
544 |
#if 1 |
545 |
|
546 |
/* Define other structures to keep mesh_write from crashing */ |
547 |
/* Change the typeid from NoType later */ |
548 |
|
549 |
mesh_p->FaceElements=Finley_ElementFile_alloc(NoType, mesh_p->order, mesh_p->reduced_order, mpi_info); |
550 |
Finley_ElementFile_allocTable(mesh_p->FaceElements, 0); |
551 |
|
552 |
mesh_p->ContactElements=Finley_ElementFile_alloc(NoType, mesh_p->order, mesh_p->reduced_order, mpi_info); |
553 |
Finley_ElementFile_allocTable(mesh_p->ContactElements, 0); |
554 |
|
555 |
mesh_p->Points=Finley_ElementFile_alloc(NoType, mesh_p->order, mesh_p->reduced_order, mpi_info); |
556 |
Finley_ElementFile_allocTable(mesh_p->Points, 0); |
557 |
|
558 |
#endif |
559 |
|
560 |
#if 0 /* comment out the rest of the un-implemented crap for now */ |
561 |
/* get the face elements */ |
562 |
if (Finley_noError()) { |
563 |
fscanf(fileHandle_p, "%s %d\n", element_type, &numEle); |
564 |
faceTypeID=Finley_RefElement_getTypeId(element_type); |
565 |
if (faceTypeID==NoType) { |
566 |
sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s for face elements",element_type); |
567 |
Finley_setError(VALUE_ERROR,error_msg); |
568 |
} else { |
569 |
mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info); |
570 |
if (Finley_noError()) { |
571 |
Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle); |
572 |
if (Finley_noError()) { |
573 |
mesh_p->FaceElements->minColor=0; |
574 |
mesh_p->FaceElements->maxColor=numEle-1; |
575 |
for (i0 = 0; i0 < numEle; i0++) { |
576 |
fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]); |
577 |
mesh_p->FaceElements->Color[i0]=i0; |
578 |
for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) { |
579 |
fscanf(fileHandle_p, " %d", |
580 |
&mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]); |
581 |
} /* for i1 */ |
582 |
fscanf(fileHandle_p, "\n"); |
583 |
} /* for i0 */ |
584 |
} |
585 |
} |
586 |
} |
587 |
} |
588 |
/* get the Contact face element */ |
589 |
if (Finley_noError()) { |
590 |
fscanf(fileHandle_p, "%s %d\n", element_type, &numEle); |
591 |
contactTypeID=Finley_RefElement_getTypeId(element_type); |
592 |
if (contactTypeID==NoType) { |
593 |
sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type); |
594 |
Finley_setError(VALUE_ERROR,error_msg); |
595 |
} else { |
596 |
mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info); |
597 |
if (Finley_noError()) { |
598 |
Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle); |
599 |
if (Finley_noError()) { |
600 |
mesh_p->ContactElements->minColor=0; |
601 |
mesh_p->ContactElements->maxColor=numEle-1; |
602 |
for (i0 = 0; i0 < numEle; i0++) { |
603 |
fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]); |
604 |
mesh_p->ContactElements->Color[i0]=i0; |
605 |
for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) { |
606 |
fscanf(fileHandle_p, " %d", |
607 |
&mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]); |
608 |
} /* for i1 */ |
609 |
fscanf(fileHandle_p, "\n"); |
610 |
} /* for i0 */ |
611 |
} |
612 |
} |
613 |
} |
614 |
} |
615 |
/* get the nodal element */ |
616 |
if (Finley_noError()) { |
617 |
fscanf(fileHandle_p, "%s %d\n", element_type, &numEle); |
618 |
pointTypeID=Finley_RefElement_getTypeId(element_type); |
619 |
if (pointTypeID==NoType) { |
620 |
sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type); |
621 |
Finley_setError(VALUE_ERROR,error_msg); |
622 |
} |
623 |
mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info); |
624 |
if (Finley_noError()) { |
625 |
Finley_ElementFile_allocTable(mesh_p->Points, numEle); |
626 |
if (Finley_noError()) { |
627 |
mesh_p->Points->minColor=0; |
628 |
mesh_p->Points->maxColor=numEle-1; |
629 |
for (i0 = 0; i0 < numEle; i0++) { |
630 |
fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]); |
631 |
mesh_p->Points->Color[i0]=i0; |
632 |
for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) { |
633 |
fscanf(fileHandle_p, " %d", |
634 |
&mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]); |
635 |
} /* for i1 */ |
636 |
fscanf(fileHandle_p, "\n"); |
637 |
} /* for i0 */ |
638 |
} |
639 |
} |
640 |
} |
641 |
/* get the name tags */ |
642 |
if (Finley_noError()) { |
643 |
if (feof(fileHandle_p) == 0) { |
644 |
fscanf(fileHandle_p, "%s\n", name); |
645 |
while (feof(fileHandle_p) == 0) { |
646 |
fscanf(fileHandle_p, "%s %d\n", name, &tag_key); |
647 |
Finley_Mesh_addTagMap(mesh_p,name,tag_key); |
648 |
} |
649 |
} |
650 |
} |
651 |
#endif /* comment out the rest of the un-implemented crap for now */ |
652 |
} |
653 |
/* close file */ |
654 |
if (mpi_info->rank == 0) fclose(fileHandle_p); |
655 |
|
656 |
/* resolve id's : */ |
657 |
/* rearrange elements: */ |
658 |
|
659 |
if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p); |
660 |
if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize); |
661 |
return mesh_p; /* ksteube temp return for debugging */ |
662 |
|
663 |
/* that's it */ |
664 |
#ifdef Finley_TRACE |
665 |
printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0); |
666 |
#endif |
667 |
|
668 |
/* that's it */ |
669 |
if (! Finley_noError()) { |
670 |
Finley_Mesh_free(mesh_p); |
671 |
} |
672 |
Paso_MPIInfo_free( mpi_info ); |
673 |
return mesh_p; |
674 |
} |