/[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 1660 - (show annotations)
Mon Jul 21 03:23:46 2008 UTC (11 years, 2 months ago) by ksteube
File MIME type: text/plain
File size: 28958 byte(s)
Some variables were not declared as a result of recent cleanup

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 int mpi_error;
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], mpi_error;
415 temp1[0] = (int) 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 = (ElementTypeId) 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 if (numEle % mpi_info->size != 0) chunkSize++; /* Remainder from numEle / mpi_info->size will be spread out one-per-CPU */
440 if (mpi_info->rank == 0) { /* Master */
441 for (;;) { /* Infinite loop */
442 chunkEle = 0;
443 for (i0=0; i0<numEle*(2+numNodes)+1; i0++) tempInts[i0] = -1;
444 for (i0=0; i0<chunkSize; i0++) {
445 if (totalEle >= numEle) break; /* End infinite loop */
446 fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
447 for (i1 = 0; i1 < numNodes; i1++) fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
448 fscanf(fileHandle_p, "\n");
449 totalEle++;
450 chunkEle++;
451 }
452 /* Eventually we'll send chunk of nodes to each CPU except 0 itself, here goes one of them */
453 if (nextCPU < mpi_info->size) {
454 #ifdef PASO_MPI
455 int mpi_error;
456
457 tempInts[numEle*(2+numNodes)] = chunkEle;
458 printf("ksteube CPU=%d/%d send to %d\n", mpi_info->rank, mpi_info->size, nextCPU);
459 mpi_error = MPI_Send(tempInts, numEle*(2+numNodes)+1, MPI_INT, nextCPU, 81722, mpi_info->comm);
460 if ( mpi_error != MPI_SUCCESS ) {
461 Finley_setError(PASO_MPI_ERROR, "Finley_Mesh_read: send of tempInts for Elements failed");
462 return NULL;
463 }
464 #endif
465 nextCPU++;
466 }
467 if (totalEle >= numEle) break; /* End infinite loop */
468 } /* Infinite loop */
469 } /* End master */
470 else { /* Worker */
471 #ifdef PASO_MPI
472 /* Each worker receives two messages */
473 MPI_Status status;
474 int mpi_error;
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 /* Change the typeid from NoType later */
550
551 mesh_p->FaceElements=Finley_ElementFile_alloc(NoType, mesh_p->order, mesh_p->reduced_order, mpi_info);
552 Finley_ElementFile_allocTable(mesh_p->FaceElements, 0);
553
554 mesh_p->ContactElements=Finley_ElementFile_alloc(NoType, mesh_p->order, mesh_p->reduced_order, mpi_info);
555 Finley_ElementFile_allocTable(mesh_p->ContactElements, 0);
556
557 mesh_p->Points=Finley_ElementFile_alloc(NoType, mesh_p->order, mesh_p->reduced_order, mpi_info);
558 Finley_ElementFile_allocTable(mesh_p->Points, 0);
559
560 #endif
561
562 #if 0 /* comment out the rest of the un-implemented crap for now */
563 /* get the face elements */
564 if (Finley_noError()) {
565 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
566 faceTypeID=Finley_RefElement_getTypeId(element_type);
567 if (faceTypeID==NoType) {
568 sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s for face elements",element_type);
569 Finley_setError(VALUE_ERROR,error_msg);
570 } else {
571 mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
572 if (Finley_noError()) {
573 Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);
574 if (Finley_noError()) {
575 mesh_p->FaceElements->minColor=0;
576 mesh_p->FaceElements->maxColor=numEle-1;
577 for (i0 = 0; i0 < numEle; i0++) {
578 fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);
579 mesh_p->FaceElements->Color[i0]=i0;
580 for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {
581 fscanf(fileHandle_p, " %d",
582 &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);
583 } /* for i1 */
584 fscanf(fileHandle_p, "\n");
585 } /* for i0 */
586 }
587 }
588 }
589 }
590 /* get the Contact face element */
591 if (Finley_noError()) {
592 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
593 contactTypeID=Finley_RefElement_getTypeId(element_type);
594 if (contactTypeID==NoType) {
595 sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);
596 Finley_setError(VALUE_ERROR,error_msg);
597 } else {
598 mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
599 if (Finley_noError()) {
600 Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);
601 if (Finley_noError()) {
602 mesh_p->ContactElements->minColor=0;
603 mesh_p->ContactElements->maxColor=numEle-1;
604 for (i0 = 0; i0 < numEle; i0++) {
605 fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);
606 mesh_p->ContactElements->Color[i0]=i0;
607 for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {
608 fscanf(fileHandle_p, " %d",
609 &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);
610 } /* for i1 */
611 fscanf(fileHandle_p, "\n");
612 } /* for i0 */
613 }
614 }
615 }
616 }
617 /* get the nodal element */
618 if (Finley_noError()) {
619 fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
620 pointTypeID=Finley_RefElement_getTypeId(element_type);
621 if (pointTypeID==NoType) {
622 sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type);
623 Finley_setError(VALUE_ERROR,error_msg);
624 }
625 mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order, mesh_p->reduced_order, mpi_info);
626 if (Finley_noError()) {
627 Finley_ElementFile_allocTable(mesh_p->Points, numEle);
628 if (Finley_noError()) {
629 mesh_p->Points->minColor=0;
630 mesh_p->Points->maxColor=numEle-1;
631 for (i0 = 0; i0 < numEle; i0++) {
632 fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);
633 mesh_p->Points->Color[i0]=i0;
634 for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {
635 fscanf(fileHandle_p, " %d",
636 &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);
637 } /* for i1 */
638 fscanf(fileHandle_p, "\n");
639 } /* for i0 */
640 }
641 }
642 }
643 /* get the name tags */
644 if (Finley_noError()) {
645 if (feof(fileHandle_p) == 0) {
646 fscanf(fileHandle_p, "%s\n", name);
647 while (feof(fileHandle_p) == 0) {
648 fscanf(fileHandle_p, "%s %d\n", name, &tag_key);
649 Finley_Mesh_addTagMap(mesh_p,name,tag_key);
650 }
651 }
652 }
653 #endif /* comment out the rest of the un-implemented crap for now */
654 }
655 /* close file */
656 if (mpi_info->rank == 0) fclose(fileHandle_p);
657
658 /* resolve id's : */
659 /* rearrange elements: */
660
661 if (Finley_noError()) Finley_Mesh_resolveNodeIds(mesh_p);
662 if (Finley_noError()) Finley_Mesh_prepare(mesh_p, optimize);
663 return mesh_p; /* ksteube temp return for debugging */
664
665 /* that's it */
666 #ifdef Finley_TRACE
667 printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
668 #endif
669
670 /* that's it */
671 if (! Finley_noError()) {
672 Finley_Mesh_free(mesh_p);
673 }
674 Paso_MPIInfo_free( mpi_info );
675 return mesh_p;
676 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26