/[escript]/trunk/finley/src/Mesh_read.c
ViewVC logotype

Contents of /trunk/finley/src/Mesh_read.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1628 - (show annotations)
Fri Jul 11 13:12:46 2008 UTC (11 years, 2 months ago) by phornby
File MIME type: text/plain
File size: 28803 byte(s)

Merge in /branches/windows_from_1456_trunk_1620_merged_in branch.

You will find a preserved pre-merge trunk in tags under tags/trunk_at_1625.
That will be useful for diffing & checking on my stupidity.

Here is a list of the conflicts and their resolution at this
point in time.


=================================================================================
(LLWS == looks like white space).

finley/src/Assemble_addToSystemMatrix.c - resolve to branch - unused var. may be wrong.....
finley/src/CPPAdapter/SystemMatrixAdapter.cpp - resolve to branch - LLWS
finley/src/CPPAdapter/MeshAdapter.cpp - resolve to branch - LLWS
paso/src/PCG.c - resolve to branch - unused var fixes.
paso/src/SolverFCT.c - resolve to branch - LLWS
paso/src/FGMRES.c - resolve to branch - LLWS
paso/src/Common.h - resolve to trunk version. It's omp.h's include... not sure it's needed,
but for the sake of saftey.....
paso/src/Functions.c - resolve to branch version, indentation/tab removal and return error
on bad unimplemented Paso_FunctionCall.
paso/src/SolverFCT_solve.c - resolve to branch version, unused vars
paso/src/SparseMatrix_MatrixVector.c - resolve to branch version, unused vars.
escript/src/Utils.cpp - resloved to branch, needs WinSock2.h
escript/src/DataExpanded.cpp - resolved to branch version - LLWS
escript/src/DataFactory.cpp - resolve to branch version
=================================================================================

This currently passes tests on linux (debian), but is not checked on windows or Altix yet.

This checkin is to make a trunk I can check out for windows to do tests on it.

Known outstanding problem is in the operator=() method of exceptions
causing warning messages on the intel compilers.

May the God of doughnuts have mercy on my soul.


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 /* these are in unimplemented code below */
247 #if 0
248 index_t tag_key;
249 ElementTypeId faceTypeID, contactTypeID, pointTypeID;
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 #ifdef PASO_MPI
299 int mpi_error;
300 #endif
301 int *tempInts = TMPMEMALLOC(numNodes*3+1, index_t);
302 double *tempCoords = TMPMEMALLOC(numNodes*numDim, double);
303
304 /*
305 Read a chunk of nodes, send to worker CPU if available, copy chunk into local mesh_p
306 It doesn't matter that a CPU has the wrong nodes for its elements, this is sorted out later
307 First chunk sent to CPU 1, second to CPU 2, ...
308 Last chunk stays on CPU 0 (the master)
309 The three columns of integers (Id, gDOF, Tag) are gathered into a single array tempInts and sent together in a single MPI message
310 */
311
312 if (mpi_info->rank == 0) { /* Master */
313 for (;;) { /* Infinite loop */
314 chunkNodes = 0;
315 for (i0=0; i0<numNodes*3+1; i0++) tempInts[i0] = -1;
316 for (i0=0; i0<numNodes*numDim; i0++) tempCoords[i0] = -1.0;
317 for (i1=0; i1<chunkSize; i1++) {
318 if (totalNodes >= numNodes) break;
319 if (1 == numDim)
320 fscanf(fileHandle_p, "%d %d %d %le\n",
321 &tempInts[0+i1], &tempInts[numNodes+i1], &tempInts[numNodes*2+i1],
322 &tempCoords[i1*numDim+0]);
323 if (2 == numDim)
324 fscanf(fileHandle_p, "%d %d %d %le %le\n",
325 &tempInts[0+i1], &tempInts[numNodes+i1], &tempInts[numNodes*2+i1],
326 &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1]);
327 if (3 == numDim)
328 fscanf(fileHandle_p, "%d %d %d %le %le %le\n",
329 &tempInts[0+i1], &tempInts[numNodes+i1], &tempInts[numNodes*2+i1],
330 &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1], &tempCoords[i1*numDim+2]);
331 totalNodes++;
332 chunkNodes++;
333 }
334 /* Eventually we'll send chunk of nodes to each CPU numbered 1 ... mpi_info->size-1, here goes one of them */
335 if (nextCPU < mpi_info->size) {
336 #ifdef PASO_MPI
337 tempInts[numNodes*3] = chunkNodes;
338 /* ksteube The size of this message can and should be brought down to chunkNodes*3+1, must re-org tempInts */
339 mpi_error = MPI_Send(tempInts, numNodes*3+1, MPI_INT, nextCPU, 81720, mpi_info->comm);
340 if ( mpi_error != MPI_SUCCESS ) {
341 Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts failed");
342 return NULL;
343 }
344 mpi_error = MPI_Send(tempCoords, numNodes*numDim, MPI_DOUBLE, nextCPU, 81721, mpi_info->comm);
345 if ( mpi_error != MPI_SUCCESS ) {
346 Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempCoords failed");
347 return NULL;
348 }
349 #endif
350 nextCPU++;
351 }
352 if (totalNodes >= numNodes) break;
353 } /* Infinite loop */
354 } /* End master */
355 else { /* Worker */
356 #ifdef PASO_MPI
357 /* Each worker receives two messages */
358 MPI_Status status;
359 mpi_error = MPI_Recv(tempInts, numNodes*3+1, MPI_INT, 0, 81720, mpi_info->comm, &status);
360 if ( mpi_error != MPI_SUCCESS ) {
361 Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts failed");
362 return NULL;
363 }
364 mpi_error = MPI_Recv(tempCoords, numNodes*numDim, MPI_DOUBLE, 0, 81721, mpi_info->comm, &status);
365 if ( mpi_error != MPI_SUCCESS ) {
366 Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempCoords failed");
367 return NULL;
368 }
369 chunkNodes = tempInts[numNodes*3];
370 #endif
371 } /* Worker */
372
373 #if 0
374 /* Display the temp mem for debugging */
375 printf("ksteube tempInts totalNodes=%d:\n", totalNodes);
376 for (i0=0; i0<numNodes*3; i0++) {
377 printf(" %2d", tempInts[i0]);
378 if (i0%numNodes==numNodes-1) printf("\n");
379 }
380 printf("ksteube tempCoords:\n");
381 for (i0=0; i0<chunkNodes*numDim; i0++) {
382 printf(" %20.15e", tempCoords[i0]);
383 if (i0%numDim==numDim-1) printf("\n");
384 }
385 #endif
386
387 printf("ksteube chunkNodes=%d numNodes=%d\n", chunkNodes, numNodes);
388 /* Copy node data from tempMem to mesh_p */
389 Finley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes);
390 if (Finley_noError()) {
391 for (i0=0; i0<chunkNodes; i0++) {
392 mesh_p->Nodes->Id[i0] = tempInts[0+i0];
393 mesh_p->Nodes->globalDegreesOfFreedom[i0] = tempInts[numNodes+i0];
394 mesh_p->Nodes->Tag[i0] = tempInts[numNodes*2+i0];
395 for (i1=0; i1<numDim; i1++) {
396 mesh_p->Nodes->Coordinates[INDEX2(i1,i0,numDim)] = tempCoords[i0*numDim+i1];
397 }
398 }
399 }
400
401 TMPMEMFREE(tempInts);
402 TMPMEMFREE(tempCoords);
403
404 /* read elements */
405
406 /* Read the element typeID */
407 if (Finley_noError()) {
408 if (mpi_info->rank == 0) {
409 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
410 typeID=Finley_RefElement_getTypeId(element_type);
411 }
412 #ifdef PASO_MPI
413 if (mpi_info->size > 0) {
414 int temp1[3];
415 temp1[0] = typeID;
416 temp1[1] = numEle;
417 mpi_error = MPI_Bcast (temp1, 2, MPI_INT, 0, mpi_info->comm);
418 if (mpi_error != MPI_SUCCESS) {
419 Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: broadcast of Element typeID failed");
420 return NULL;
421 }
422 typeID = temp1[0];
423 numEle = temp1[1];
424 }
425 #endif
426 if (typeID==NoType) {
427 sprintf(error_msg, "Finley_Mesh_read: Unidentified element type %s", element_type);
428 Finley_setError(VALUE_ERROR, error_msg);
429 }
430 }
431
432 /* Read the element data */
433 mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
434 numNodes = mesh_p->Elements->ReferenceElement->Type->numNodes; /* New meaning for numNodes: num nodes per element */
435
436 if (Finley_noError()) {
437 int *tempInts = TMPMEMALLOC(numEle*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
438 int chunkSize = numEle / mpi_info->size, totalEle=0, nextCPU=1;
439 #ifdef PASO_MPI
440 int mpi_error;
441 #endif
442 if (numEle % mpi_info->size != 0) chunkSize++; /* Remainder from numEle / mpi_info->size will be spread out one-per-CPU */
443 if (mpi_info->rank == 0) { /* Master */
444 for (;;) { /* Infinite loop */
445 chunkEle = 0;
446 for (i0=0; i0<numEle*(2+numNodes)+1; i0++) tempInts[i0] = -1;
447 for (i0=0; i0<chunkSize; i0++) {
448 if (totalEle >= numEle) break; /* End infinite loop */
449 fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
450 for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
451 fscanf(fileHandle_p, "\n");
452 totalEle++;
453 chunkEle++;
454 }
455 /* Eventually we'll send chunk of nodes to each CPU except 0 itself, here goes one of them */
456 if (nextCPU < mpi_info->size) {
457 #ifdef PASO_MPI
458 tempInts[numEle*(2+numNodes)] = chunkEle;
459 printf("ksteube CPU=%d/%d send to %d\n", mpi_info->rank, mpi_info->size, nextCPU);
460 mpi_error = MPI_Send(tempInts, numEle*(2+numNodes)+1, MPI_INT, nextCPU, 81722, mpi_info->comm);
461 if ( mpi_error != MPI_SUCCESS ) {
462 Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Elements failed");
463 return NULL;
464 }
465 #endif
466 nextCPU++;
467 }
468 if (totalEle >= numEle) break; /* End infinite loop */
469 } /* Infinite loop */
470 } /* End master */
471 else { /* Worker */
472 #ifdef PASO_MPI
473 /* Each worker receives two messages */
474 MPI_Status status;
475 printf("ksteube CPU=%d/%d recv on %d\n", mpi_info->rank, mpi_info->size, mpi_info->rank);
476 mpi_error = MPI_Recv(tempInts, numEle*(2+numNodes)+1, MPI_INT, 0, 81722, mpi_info->comm, &status);
477 if ( mpi_error != MPI_SUCCESS ) {
478 Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: receive of tempInts for Elements failed");
479 return NULL;
480 }
481 chunkEle = tempInts[numEle*(2+numNodes)];
482 #endif
483 } /* Worker */
484 #if 1
485 /* Display the temp mem for debugging */
486 printf("ksteube tempInts numEle=%d chunkEle=%d AAA:\n", numEle, chunkEle);
487 for (i0=0; i0<numEle*(numNodes+2); i0++) {
488 printf(" %2d", tempInts[i0]);
489 if (i0%(numNodes+2)==numNodes+2-1) printf("\n");
490 }
491 #endif
492
493 /* Copy Element data from tempInts to mesh_p */
494 Finley_ElementFile_allocTable(mesh_p->Elements, chunkEle);
495 mesh_p->Elements->minColor=0;
496 mesh_p->Elements->maxColor=chunkEle-1;
497 if (Finley_noError()) {
498 #pragma omp parallel for private (i0, i1)
499 for (i0=0; i0<chunkEle; i0++) {
500 mesh_p->Elements->Id[i0] = tempInts[i0*(2+numNodes)+0];
501 mesh_p->Elements->Tag[i0] = tempInts[i0*(2+numNodes)+1];
502 for (i1 = 0; i1 < numNodes; i1++) {
503 mesh_p->Elements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
504 }
505 }
506 }
507
508 TMPMEMFREE(tempInts);
509 }
510
511
512 #if 0 /* this is the original code for reading elements */
513 /* read elements */
514 if (Finley_noError()) {
515
516 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
517 typeID=Finley_RefElement_getTypeId(element_type);
518 if (typeID==NoType) {
519 sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s",element_type);
520 Finley_setError(VALUE_ERROR,error_msg);
521 } else {
522 /* read the elements */
523 mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
524 if (Finley_noError()) {
525 Finley_ElementFile_allocTable(mesh_p->Elements, numEle);
526 mesh_p->Elements->minColor=0;
527 mesh_p->Elements->maxColor=numEle-1;
528 if (Finley_noError()) {
529 for (i0 = 0; i0 < numEle; i0++) {
530 fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);
531 mesh_p->Elements->Color[i0]=i0;
532 for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {
533 fscanf(fileHandle_p, " %d",
534 &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);
535 } /* for i1 */
536 fscanf(fileHandle_p, "\n");
537 } /* for i0 */
538 }
539 }
540 }
541 }
542 #endif
543
544 printf("ksteube CPU=%d/%d Element typeID=%d\n", mpi_info->rank, mpi_info->size, typeID);
545
546 #if 1
547
548 /* Define other structures to keep mesh_write from crashing */
549
550 mesh_p->FaceElements=Finley_ElementFile_alloc(0, mesh_p->order, mesh_p->reduced_order, mpi_info);
551 Finley_ElementFile_allocTable(mesh_p->FaceElements, 0);
552
553 mesh_p->ContactElements=Finley_ElementFile_alloc(0, mesh_p->order, mesh_p->reduced_order, mpi_info);
554 Finley_ElementFile_allocTable(mesh_p->ContactElements, 0);
555
556 mesh_p->Points=Finley_ElementFile_alloc(0, mesh_p->order, mesh_p->reduced_order, mpi_info);
557 Finley_ElementFile_allocTable(mesh_p->Points, 0);
558
559 #endif
560
561 #if 0 /* comment out the rest of the un-implemented crap for now */
562 /* get the face elements */
563 if (Finley_noError()) {
564 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
565 faceTypeID=Finley_RefElement_getTypeId(element_type);
566 if (faceTypeID==NoType) {
567 sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s for face elements",element_type);
568 Finley_setError(VALUE_ERROR,error_msg);
569 } else {
570 mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
571 if (Finley_noError()) {
572 Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);
573 if (Finley_noError()) {
574 mesh_p->FaceElements->minColor=0;
575 mesh_p->FaceElements->maxColor=numEle-1;
576 for (i0 = 0; i0 < numEle; i0++) {
577 fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);
578 mesh_p->FaceElements->Color[i0]=i0;
579 for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {
580 fscanf(fileHandle_p, " %d",
581 &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);
582 } /* for i1 */
583 fscanf(fileHandle_p, "\n");
584 } /* for i0 */
585 }
586 }
587 }
588 }
589 /* get the Contact face element */
590 if (Finley_noError()) {
591 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
592 contactTypeID=Finley_RefElement_getTypeId(element_type);
593 if (contactTypeID==NoType) {
594 sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);
595 Finley_setError(VALUE_ERROR,error_msg);
596 } else {
597 mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
598 if (Finley_noError()) {
599 Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);
600 if (Finley_noError()) {
601 mesh_p->ContactElements->minColor=0;
602 mesh_p->ContactElements->maxColor=numEle-1;
603 for (i0 = 0; i0 < numEle; i0++) {
604 fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);
605 mesh_p->ContactElements->Color[i0]=i0;
606 for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {
607 fscanf(fileHandle_p, " %d",
608 &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);
609 } /* for i1 */
610 fscanf(fileHandle_p, "\n");
611 } /* for i0 */
612 }
613 }
614 }
615 }
616 /* get the nodal element */
617 if (Finley_noError()) {
618 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
619 pointTypeID=Finley_RefElement_getTypeId(element_type);
620 if (pointTypeID==NoType) {
621 sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type);
622 Finley_setError(VALUE_ERROR,error_msg);
623 }
624 mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
625 if (Finley_noError()) {
626 Finley_ElementFile_allocTable(mesh_p->Points, numEle);
627 if (Finley_noError()) {
628 mesh_p->Points->minColor=0;
629 mesh_p->Points->maxColor=numEle-1;
630 for (i0 = 0; i0 < numEle; i0++) {
631 fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);
632 mesh_p->Points->Color[i0]=i0;
633 for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {
634 fscanf(fileHandle_p, " %d",
635 &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);
636 } /* for i1 */
637 fscanf(fileHandle_p, "\n");
638 } /* for i0 */
639 }
640 }
641 }
642 /* get the name tags */
643 if (Finley_noError()) {
644 if (feof(fileHandle_p) == 0) {
645 fscanf(fileHandle_p, "%s\n", name);
646 while (feof(fileHandle_p) == 0) {
647 fscanf(fileHandle_p, "%s %d\n", name, &tag_key);
648 Finley_Mesh_addTagMap(mesh_p,name,tag_key);
649 }
650 }
651 }
652 #endif /* comment out the rest of the un-implemented crap for now */
653 }
654 /* close file */
655 if (mpi_info->rank == 0) fclose(fileHandle_p);
656
657 /* resolve id's : */
658 /* rearrange elements: */
659
660 if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
661 if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
662 return mesh_p; /* ksteube temp return for debugging */
663
664 /* that's it */
665 #ifdef Finley_TRACE
666 printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
667 #endif
668
669 /* that's it */
670 if (! Finley_noError()) {
671 Finley_Mesh_free(mesh_p);
672 }
673 Paso_MPIInfo_free( mpi_info );
674 return mesh_p;
675 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26