/[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 1402 - (show annotations)
Thu Jan 31 00:17:49 2008 UTC (11 years, 3 months ago) by ksteube
Original Path: trunk/finley/src/Mesh_read.c
File MIME type: text/plain
File size: 28671 byte(s)
now implemented to read/fan out elements

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26