/[escript]/branches/arrayview_from_1695_trunk/finley/src/Mesh_read.c
ViewVC logotype

Contents of /branches/arrayview_from_1695_trunk/finley/src/Mesh_read.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1646 - (show annotations)
Tue Jul 15 05:27:39 2008 UTC (10 years, 10 months ago) by phornby
Original Path: trunk/finley/src/Mesh_read.c
File MIME type: text/plain
File size: 28882 byte(s)
Unused vars due to the #if 0's in the code.
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 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26