/[escript]/branches/domexper/dudley/src/Mesh_read.c
ViewVC logotype

Contents of /branches/domexper/dudley/src/Mesh_read.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3220 - (show annotations)
Wed Sep 29 00:33:16 2010 UTC (8 years, 7 months ago) by jfenwick
File MIME type: text/plain
File size: 21862 byte(s)
Removing references to Quadrature.? and ShapeFns

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14
15 /**************************************************************/
16
17 /* Dudley: read mesh */
18
19 /**************************************************************/
20
21 #include <ctype.h>
22 #include "Mesh.h"
23
24 #define FSCANF_CHECK(scan_ret, reason) { if (scan_ret == EOF) { perror(reason); Dudley_setError(IO_ERROR,"scan error while reading dudley file"); return NULL;} }
25
26
27 Dudley_Mesh* Dudley_Mesh_read(char* fname,index_t order, index_t reduced_order, bool_t optimize)
28 {
29 Paso_MPIInfo *mpi_info = NULL;
30 dim_t numNodes, numDim, numEle, i0, i1;
31 Dudley_Mesh *mesh_p=NULL;
32 // Dudley_ReferenceElementSet *refPoints=NULL, *refFaceElements=NULL, *refElements=NULL;
33 char name[LenString_MAX],element_type[LenString_MAX],frm[20];
34 char error_msg[LenErrorMsg_MAX];
35 FILE *fileHandle_p = NULL;
36 ElementTypeId typeID=NoRef;
37 int scan_ret;
38
39 Dudley_resetError();
40 mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
41
42 if (mpi_info->rank == 0) {
43 /* get file handle */
44 fileHandle_p = fopen(fname, "r");
45 if (fileHandle_p==NULL) {
46 sprintf(error_msg,"Dudley_Mesh_read: Opening file %s for reading failed.",fname);
47 Dudley_setError(IO_ERROR,error_msg);
48 Paso_MPIInfo_free( mpi_info );
49 return NULL;
50 }
51
52 /* read header */
53 sprintf(frm,"%%%d[^\n]",LenString_MAX-1);
54 scan_ret = fscanf(fileHandle_p, frm, name);
55 FSCANF_CHECK(scan_ret, "Dudley_Mesh_read")
56
57 /* get the number of nodes */
58 scan_ret = fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);
59 FSCANF_CHECK(scan_ret, "Dudley_Mesh_read")
60 }
61
62 #ifdef PASO_MPI
63 /* MPI Broadcast numDim, numNodes, name if there are multiple MPI procs*/
64 if (mpi_info->size > 1)
65 {
66 int temp1[3];
67 if (mpi_info->rank == 0)
68 {
69 temp1[0] = numDim;
70 temp1[1] = numNodes;
71 temp1[2] = strlen(name) + 1;
72 } else {
73 temp1[0] = 0;
74 temp1[1] = 0;
75 temp1[2] = 1;
76 }
77 MPI_Bcast (temp1, 3, MPI_INT, 0, mpi_info->comm);
78 numDim = temp1[0];
79 numNodes = temp1[1];
80 MPI_Bcast (name, temp1[2], MPI_CHAR, 0, mpi_info->comm);
81 }
82 #endif
83
84 /* allocate mesh */
85 mesh_p = Dudley_Mesh_alloc(name,numDim,mpi_info);
86
87 if (Dudley_noError())
88 {
89 /* Each CPU will get at most chunkSize nodes so the message has to be sufficiently large */
90 int chunkSize = numNodes / mpi_info->size + 1, totalNodes=0, chunkNodes=0, nextCPU=1;
91 int *tempInts = TMPMEMALLOC(chunkSize*3+1, index_t); /* Stores the integer message data */
92 double *tempCoords = TMPMEMALLOC(chunkSize*numDim, double); /* Stores the double message data */
93
94 /*
95 Read chunkSize nodes, send it in a chunk to worker CPU which copies chunk into its local mesh_p
96 It doesn't matter that a CPU has the wrong nodes for its elements, this is sorted out later
97 First chunk sent to CPU 1, second to CPU 2, ...
98 Last chunk stays on CPU 0 (the master)
99 The three columns of integers (Id, gDOF, Tag) are gathered into a single array tempInts and sent
100 together in a single MPI message
101 */
102
103 if (mpi_info->rank == 0) /* Master */
104 {
105 for (;;) /* Infinite loop */
106 {
107 #pragma omp parallel for private (i0) schedule(static)
108 for (i0=0; i0<chunkSize*3+1; i0++) tempInts[i0] = -1;
109
110 #pragma omp parallel for private (i0) schedule(static)
111 for (i0=0; i0<chunkSize*numDim; i0++) tempCoords[i0] = -1.0;
112
113 chunkNodes = 0;
114 for (i1=0; i1<chunkSize; i1++)
115 {
116 if (totalNodes >= numNodes) break; /* End of inner loop */
117 if (1 == numDim) {
118 scan_ret = fscanf(fileHandle_p, "%d %d %d %le\n",
119 &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
120 &tempCoords[i1*numDim+0]);
121 FSCANF_CHECK(scan_ret, "Dudley_Mesh_read")
122 }
123 if (2 == numDim)
124 {
125 scan_ret = fscanf(fileHandle_p, "%d %d %d %le %le\n",
126 &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
127 &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1]);
128 FSCANF_CHECK(scan_ret, "Dudley_Mesh_read")
129 }
130 if (3 == numDim)
131 {
132 scan_ret = fscanf(fileHandle_p, "%d %d %d %le %le %le\n",
133 &tempInts[0+i1], &tempInts[chunkSize+i1], &tempInts[chunkSize*2+i1],
134 &tempCoords[i1*numDim+0], &tempCoords[i1*numDim+1], &tempCoords[i1*numDim+2]);
135 FSCANF_CHECK(scan_ret, "Dudley_Mesh_read")
136 }
137 totalNodes++; /* When do we quit the infinite loop? */
138 chunkNodes++; /* How many nodes do we actually have in this chunk? It may be smaller than chunkSize. */
139 }
140 if (chunkNodes > chunkSize)
141 {
142 Dudley_setError(PASO_MPI_ERROR, "Dudley_Mesh_read: error reading chunks of mesh, data too large for message size");
143 return NULL;
144 }
145 #ifdef PASO_MPI
146 /* Eventually we'll send chunkSize nodes to each CPU numbered 1 ... mpi_info->size-1, here goes one of them */
147 if (nextCPU < mpi_info->size)
148 {
149 tempInts[chunkSize*3] = chunkNodes; /* The message has one more int to send chunkNodes */
150 MPI_Send(tempInts, chunkSize*3+1, MPI_INT, nextCPU, 81720, mpi_info->comm);
151 MPI_Send(tempCoords, chunkSize*numDim, MPI_DOUBLE, nextCPU, 81721, mpi_info->comm);
152 }
153 #endif
154 nextCPU++;
155 /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
156 if (nextCPU > mpi_info->size) break; /* End infinite loop */
157 } /* Infinite loop */
158 } /* End master */
159 else /* Worker */
160 {
161 #ifdef PASO_MPI
162 /* Each worker receives two messages */
163 MPI_Status status;
164 MPI_Recv(tempInts, chunkSize*3+1, MPI_INT, 0, 81720, mpi_info->comm, &status);
165 MPI_Recv(tempCoords, chunkSize*numDim, MPI_DOUBLE, 0, 81721, mpi_info->comm, &status);
166 chunkNodes = tempInts[chunkSize*3]; /* How many nodes are in this workers chunk? */
167 #endif
168 } /* Worker */
169
170 /* Copy node data from tempMem to mesh_p */
171 Dudley_NodeFile_allocTable(mesh_p->Nodes, chunkNodes);
172 if (Dudley_noError())
173 {
174 #pragma omp parallel for private (i0, i1) schedule(static)
175 for (i0=0; i0<chunkNodes; i0++)
176 {
177 mesh_p->Nodes->Id[i0] = tempInts[0+i0];
178 mesh_p->Nodes->globalDegreesOfFreedom[i0] = tempInts[chunkSize+i0];
179 mesh_p->Nodes->Tag[i0] = tempInts[chunkSize*2+i0];
180 for (i1=0; i1<numDim; i1++)
181 {
182 mesh_p->Nodes->Coordinates[INDEX2(i1,i0,numDim)] = tempCoords[i0*numDim+i1];
183 }
184 }
185 }
186 TMPMEMFREE(tempInts);
187 TMPMEMFREE(tempCoords);
188 }
189
190 /* *********************************** read elements ******************************************************************/
191 if (Dudley_noError())
192 {
193 /* Read the element typeID */
194 if (mpi_info->rank == 0)
195 {
196 scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
197 FSCANF_CHECK(scan_ret, "Dudley_Mesh_read")
198 typeID=eltTypeFromString(element_type);
199 }
200 #ifdef PASO_MPI
201 if (mpi_info->size > 1)
202 {
203 int temp1[2], mpi_error;
204 temp1[0] = (int) typeID;
205 temp1[1] = numEle;
206 mpi_error = MPI_Bcast (temp1, 2, MPI_INT, 0, mpi_info->comm);
207 if (mpi_error != MPI_SUCCESS)
208 {
209 Dudley_setError(PASO_MPI_ERROR, "Dudley_Mesh_read: broadcast of Element typeID failed");
210 return NULL;
211 }
212 typeID = (ElementTypeId) temp1[0];
213 numEle = temp1[1];
214 }
215 #endif
216 if (typeID==NoRef) {
217 sprintf(error_msg, "Dudley_Mesh_read: Unidentified element type %s", element_type);
218 Dudley_setError(VALUE_ERROR, error_msg);
219 }
220 }
221
222 /* Allocate the ElementFile */
223 if (Dudley_noError()) {
224 // refElements= Dudley_ReferenceElementSet_alloc(typeID,order, reduced_order);
225 mesh_p->Elements=Dudley_ElementFile_alloc(typeID, mpi_info);
226 numNodes = mesh_p->Elements->numNodes; /* New meaning for numNodes: num nodes per element */
227 }
228
229 /* *************************** Read the element data *******************************************************************/
230 if (Dudley_noError())
231 {
232 int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1, chunkEle=0;
233 int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
234 /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
235 if (mpi_info->rank == 0) /* Master */
236 {
237 for (;;) /* Infinite loop */
238 {
239 #pragma omp parallel for private (i0) schedule(static)
240 for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
241 chunkEle = 0;
242 for (i0=0; i0<chunkSize; i0++)
243 {
244 if (totalEle >= numEle) break; /* End inner loop */
245 scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
246 FSCANF_CHECK(scan_ret, "Dudley_Mesh_read")
247 for (i1 = 0; i1 < numNodes; i1++)
248 {
249 scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
250 FSCANF_CHECK(scan_ret, "Dudley_Mesh_read")
251 }
252 scan_ret = fscanf(fileHandle_p, "\n");
253 FSCANF_CHECK(scan_ret, "Dudley_Mesh_read")
254 totalEle++;
255 chunkEle++;
256 }
257 #ifdef PASO_MPI
258 /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
259 if (nextCPU < mpi_info->size) {
260 tempInts[chunkSize*(2+numNodes)] = chunkEle;
261 MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81722, mpi_info->comm);
262 }
263 #endif
264 nextCPU++;
265 /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
266 if (nextCPU > mpi_info->size) break; /* End infinite loop */
267 } /* Infinite loop */
268 } /* End master */
269 else
270 { /* Worker */
271 #ifdef PASO_MPI
272 /* Each worker receives one message */
273 MPI_Status status;
274 MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81722, mpi_info->comm, &status);
275 chunkEle = tempInts[chunkSize*(2+numNodes)];
276 #endif
277 } /* Worker */
278 Dudley_ElementFile_allocTable(mesh_p->Elements, chunkEle);
279
280 /* Copy Element data from tempInts to mesh_p */
281 if (Dudley_noError())
282 {
283 mesh_p->Elements->minColor=0;
284 mesh_p->Elements->maxColor=chunkEle-1;
285 #pragma omp parallel for private (i0, i1) schedule(static)
286 for (i0=0; i0<chunkEle; i0++)
287 {
288 mesh_p->Elements->Id[i0] = tempInts[i0*(2+numNodes)+0];
289 mesh_p->Elements->Tag[i0] = tempInts[i0*(2+numNodes)+1];
290 mesh_p->Elements->Owner[i0] =mpi_info->rank;
291 mesh_p->Elements->Color[i0] = i0;
292 for (i1 = 0; i1 < numNodes; i1++)
293 {
294 mesh_p->Elements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
295 }
296 }
297 }
298 TMPMEMFREE(tempInts);
299 }
300 /* ******************** end of Read the element data ******************************************************/
301
302 /* ********************* read face elements ***************************************************************/
303 if (Dudley_noError())
304 {
305 /* Read the element typeID */
306 if (mpi_info->rank == 0)
307 {
308 scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
309 FSCANF_CHECK(scan_ret, "Dudley_Mesh_read")
310 typeID=eltTypeFromString(element_type);
311 }
312 #ifdef PASO_MPI
313 if (mpi_info->size > 1)
314 {
315 int temp1[2];
316 temp1[0] = (int) typeID;
317 temp1[1] = numEle;
318 MPI_Bcast (temp1, 2, MPI_INT, 0, mpi_info->comm);
319 typeID = (ElementTypeId) temp1[0];
320 numEle = temp1[1];
321 }
322 #endif
323 if (typeID==NoRef)
324 {
325 sprintf(error_msg, "Dudley_Mesh_read: Unidentified element type %s", element_type);
326 Dudley_setError(VALUE_ERROR, error_msg);
327 }
328 if (Dudley_noError())
329 {
330 /* Allocate the ElementFile */
331 // refFaceElements= Dudley_ReferenceElementSet_alloc(typeID,order, reduced_order);
332 mesh_p->FaceElements=Dudley_ElementFile_alloc(typeID, mpi_info);
333 numNodes = mesh_p->FaceElements->numNodes; /* New meaning for numNodes: num nodes per element */
334 }
335 }
336 /* ********************** Read the face element data ********************************************************* */
337
338 if (Dudley_noError())
339 {
340 int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1, chunkEle=0;
341 int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
342 /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
343 if (mpi_info->rank == 0) /* Master */
344 {
345 for (;;) /* Infinite loop */
346 {
347 #pragma omp parallel for private (i0) schedule(static)
348 for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
349 chunkEle = 0;
350 for (i0=0; i0<chunkSize; i0++)
351 {
352 if (totalEle >= numEle) break; /* End inner loop */
353 scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
354 FSCANF_CHECK(scan_ret, "Dudley_Mesh_read")
355 for (i1 = 0; i1 < numNodes; i1++)
356 {
357 scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
358 FSCANF_CHECK(scan_ret, "Dudley_Mesh_read")
359 }
360 scan_ret = fscanf(fileHandle_p, "\n");
361 FSCANF_CHECK(scan_ret, "Dudley_Mesh_read")
362 totalEle++;
363 chunkEle++;
364 }
365 #ifdef PASO_MPI
366 /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
367 if (nextCPU < mpi_info->size)
368 {
369 tempInts[chunkSize*(2+numNodes)] = chunkEle;
370 MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81723, mpi_info->comm);
371 }
372 #endif
373 nextCPU++;
374 /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
375 if (nextCPU > mpi_info->size) break; /* End infinite loop */
376 } /* Infinite loop */
377 } /* End master */
378 else /* Worker */
379 {
380 #ifdef PASO_MPI
381 /* Each worker receives one message */
382 MPI_Status status;
383 MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81723, mpi_info->comm, &status);
384 chunkEle = tempInts[chunkSize*(2+numNodes)];
385 #endif
386 } /* Worker */
387 Dudley_ElementFile_allocTable(mesh_p->FaceElements, chunkEle);
388 if (Dudley_noError())
389 {
390 /* Copy Element data from tempInts to mesh_p */
391 mesh_p->FaceElements->minColor=0;
392 mesh_p->FaceElements->maxColor=chunkEle-1;
393 #pragma omp parallel for private (i0, i1)
394 for (i0=0; i0<chunkEle; i0++)
395 {
396 mesh_p->FaceElements->Id[i0] = tempInts[i0*(2+numNodes)+0];
397 mesh_p->FaceElements->Tag[i0] = tempInts[i0*(2+numNodes)+1];
398 mesh_p->FaceElements->Owner[i0] =mpi_info->rank;
399 mesh_p->FaceElements->Color[i0] = i0;
400 for (i1 = 0; i1 < numNodes; i1++)
401 {
402 mesh_p->FaceElements->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
403 }
404 }
405 }
406
407 TMPMEMFREE(tempInts);
408 }
409 /* ************************************* end of Read the face element data *************************************** */
410 /* ************************************* read contact elements ************************************************** */
411
412 // Dudley does not have contact elements
413
414 /* ********************************* read nodal elements ****************************************************** */
415 /* ******************************* Read the element typeID */
416 if (Dudley_noError())
417 {
418 if (mpi_info->rank == 0)
419 {
420 scan_ret = fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
421 FSCANF_CHECK(scan_ret, "Dudley_Mesh_read")
422 typeID=eltTypeFromString(element_type);
423 }
424 #ifdef PASO_MPI
425 if (mpi_info->size > 1)
426 {
427 int temp1[2];
428 temp1[0] = (int) typeID;
429 temp1[1] = numEle;
430 MPI_Bcast (temp1, 2, MPI_INT, 0, mpi_info->comm);
431 typeID = (ElementTypeId) temp1[0];
432 numEle = temp1[1];
433 }
434 #endif
435 if (typeID==NoRef)
436 {
437 sprintf(error_msg, "Dudley_Mesh_read: Unidentified element type %s", element_type);
438 Dudley_setError(VALUE_ERROR, error_msg);
439 }
440 }
441 if (Dudley_noError())
442 {
443 /* Allocate the ElementFile */
444 // refPoints= Dudley_ReferenceElementSet_alloc(typeID,order, reduced_order);
445 mesh_p->Points=Dudley_ElementFile_alloc(typeID, mpi_info);
446 numNodes = mesh_p->Points->numNodes; /* New meaning for numNodes: num nodes per element */
447 }
448 /********************************** Read the nodal element data **************************************************/
449 if (Dudley_noError())
450 {
451 int chunkSize = numEle / mpi_info->size + 1, totalEle=0, nextCPU=1, chunkEle=0;
452 int *tempInts = TMPMEMALLOC(chunkSize*(2+numNodes)+1, index_t); /* Store Id + Tag + node list (+ one int at end for chunkEle) */
453 /* Elements are specified as a list of integers...only need one message instead of two as with the nodes */
454 if (mpi_info->rank == 0) /* Master */
455 {
456 for (;;) /* Infinite loop */
457 {
458 #pragma omp parallel for private (i0) schedule(static)
459 for (i0=0; i0<chunkSize*(2+numNodes)+1; i0++) tempInts[i0] = -1;
460 chunkEle = 0;
461 for (i0=0; i0<chunkSize; i0++)
462 {
463 if (totalEle >= numEle) break; /* End inner loop */
464 scan_ret = fscanf(fileHandle_p, "%d %d", &tempInts[i0*(2+numNodes)+0], &tempInts[i0*(2+numNodes)+1]);
465 FSCANF_CHECK(scan_ret, "Dudley_Mesh_read")
466 for (i1 = 0; i1 < numNodes; i1++)
467 {
468 scan_ret = fscanf(fileHandle_p, " %d", &tempInts[i0*(2+numNodes)+2+i1]);
469 FSCANF_CHECK(scan_ret, "Dudley_Mesh_read")
470 }
471 scan_ret = fscanf(fileHandle_p, "\n");
472 FSCANF_CHECK(scan_ret, "Dudley_Mesh_read")
473 totalEle++;
474 chunkEle++;
475 }
476 #ifdef PASO_MPI
477 /* Eventually we'll send chunk of elements to each CPU except 0 itself, here goes one of them */
478 if (nextCPU < mpi_info->size)
479 {
480 tempInts[chunkSize*(2+numNodes)] = chunkEle;
481 MPI_Send(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, nextCPU, 81725, mpi_info->comm);
482 }
483 #endif
484 nextCPU++;
485 /* Infinite loop ends when I've read a chunk for each of the worker nodes plus one more chunk for the master */
486 if (nextCPU > mpi_info->size) break; /* End infinite loop */
487 } /* Infinite loop */
488 } /* End master */
489 else /* Worker */
490 {
491 #ifdef PASO_MPI
492 /* Each worker receives one message */
493 MPI_Status status;
494 MPI_Recv(tempInts, chunkSize*(2+numNodes)+1, MPI_INT, 0, 81725, mpi_info->comm, &status);
495 chunkEle = tempInts[chunkSize*(2+numNodes)];
496 #endif
497 } /* Worker */
498
499 /* Copy Element data from tempInts to mesh_p */
500 Dudley_ElementFile_allocTable(mesh_p->Points, chunkEle);
501 if (Dudley_noError())
502 {
503 mesh_p->Points->minColor=0;
504 mesh_p->Points->maxColor=chunkEle-1;
505 #pragma omp parallel for private (i0, i1) schedule(static)
506 for (i0=0; i0<chunkEle; i0++)
507 {
508 mesh_p->Points->Id[i0] = tempInts[i0*(2+numNodes)+0];
509 mesh_p->Points->Tag[i0] = tempInts[i0*(2+numNodes)+1];
510 mesh_p->Points->Owner[i0] =mpi_info->rank;
511 mesh_p->Points->Color[i0] = i0;
512 for (i1 = 0; i1 < numNodes; i1++)
513 {
514 mesh_p->Points->Nodes[INDEX2(i1, i0, numNodes)] = tempInts[i0*(2+numNodes)+2+i1];
515 }
516 }
517 }
518
519 TMPMEMFREE(tempInts);
520 } /* ******************************** end of Read the nodal element data ************************************ */
521 /****************** get the name tags *****************************************/
522 if (Dudley_noError())
523 {
524 char *remainder=0, *ptr;
525 size_t len=0;
526 #ifdef PASO_MPI
527 int len_i;
528 #endif
529 int tag_key;
530 if (mpi_info->rank == 0) /* Master */
531 {
532 /* Read the word 'Tag' */
533 if (! feof(fileHandle_p))
534 {
535 scan_ret = fscanf(fileHandle_p, "%s\n", name);
536 FSCANF_CHECK(scan_ret, "Dudley_Mesh_read")
537 }
538 #if defined(_WIN32) /* windows ftell lies on unix formatted text files */
539 remainder = NULL;
540 len=0;
541 while (1)
542 {
543 size_t malloc_chunk = 1024;
544 size_t buff_size = 0;
545 int ch;
546 ch = fgetc(fileHandle_p);
547 if( ch == '\r' )
548 {
549 continue;
550 }
551 if( len+1 > buff_size )
552 {
553 TMPMEMREALLOC(remainder,remainder,buff_size+malloc_chunk,char);
554 }
555 if( ch == EOF )
556 {
557 /* hit EOF */
558 remainder[len] = (char)0;
559 break;
560 }
561 remainder[len] = (char)ch;
562 len++;
563 }
564 #else
565 /* Read rest of file in one chunk, after using seek to find length */
566 {
567 long cur_pos, end_pos;
568 cur_pos = ftell(fileHandle_p);
569 fseek(fileHandle_p, 0L, SEEK_END);
570 end_pos = ftell(fileHandle_p);
571 fseek(fileHandle_p, (long)cur_pos, SEEK_SET);
572 remainder = TMPMEMALLOC(end_pos-cur_pos+1, char);
573 if (! feof(fileHandle_p))
574 {
575 scan_ret = fread(remainder, (size_t) end_pos-cur_pos,
576 sizeof(char), fileHandle_p);
577 FSCANF_CHECK(scan_ret, "Dudley_Mesh_read")
578 remainder[end_pos-cur_pos] = 0;
579 }
580 }
581 #endif
582 len = strlen(remainder);
583 while ((len>1) && isspace(remainder[--len])) {remainder[len]=0;}
584 len = strlen(remainder);
585 TMPMEMREALLOC(remainder,remainder,len+1,char);
586 } /* Master */
587 #ifdef PASO_MPI
588
589 len_i=(int) len;
590 MPI_Bcast (&len_i, 1, MPI_INT, 0, mpi_info->comm);
591 len=(size_t) len_i;
592 if (mpi_info->rank != 0)
593 {
594 remainder = TMPMEMALLOC(len+1, char);
595 remainder[0] = 0;
596 }
597 if (MPI_Bcast (remainder, len+1, MPI_CHAR, 0, mpi_info->comm)!=MPI_SUCCESS)
598 Dudley_setError(PASO_MPI_ERROR, "Dudley_Mesh_read: broadcast of remainder failed");
599 #endif
600
601 if (remainder[0])
602 {
603 ptr = remainder;
604 do
605 {
606 sscanf(ptr, "%s %d\n", name, &tag_key);
607 if (*name) Dudley_Mesh_addTagMap(mesh_p,name,tag_key);
608 ptr++;
609 } while(NULL != (ptr = strchr(ptr, '\n')) && *ptr);
610 }
611 if (remainder)
612 TMPMEMFREE(remainder);
613 }
614
615 /* close file */
616 if (mpi_info->rank == 0) fclose(fileHandle_p);
617
618 /* resolve id's : */
619 /* rearrange elements: */
620 if (Dudley_noError()) Dudley_Mesh_resolveNodeIds(mesh_p);
621 if (Dudley_noError()) Dudley_Mesh_prepare(mesh_p, optimize);
622
623 /* that's it */
624 if (! Dudley_noError())
625 {
626 Dudley_Mesh_free(mesh_p);
627 }
628 /* free up memory */
629 // Dudley_ReferenceElementSet_dealloc(refPoints);
630 // Dudley_ReferenceElementSet_dealloc(refFaceElements);
631 // Dudley_ReferenceElementSet_dealloc(refElements);
632 Paso_MPIInfo_free( mpi_info );
633 return mesh_p;
634 }
635

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26