finley/src/Mesh_read.c

/*
 ************************************************************
 *          Copyright 2006 by ACcESS MNRF                   *
 *                                                          *
 *              http://www.access.edu.au                    *
 *       Primary Business: Queensland, Australia            *
 *  Licensed under the Open Software License version 3.0    *
 *     http://www.opensource.org/licenses/osl-3.0.php       *
 *                                                          *
 ************************************************************
*/

/**************************************************************/

/*   Finley: read mesh */

/**************************************************************/

/*   Author: gross@access.edu.au */
/*   Version: $Id$ */

/**************************************************************/

#include "Mesh.h"

/**************************************************************/

/*  reads a mesh from a Finley file of name fname */

Finley_Mesh* Finley_Mesh_read(char* fname,index_t order, index_t reduced_order,  bool_t optimize_labeling) {

  dim_t numNodes, numDim, numEle, i0, i1;
  index_t tag_key;
  Finley_Mesh *mesh_p=NULL;
  char name[LenString_MAX],element_type[LenString_MAX],frm[20];
  char error_msg[LenErrorMsg_MAX];
  double time0=Finley_timer();
  FILE *fileHandle_p = NULL;
  ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;

  Finley_resetError();
#ifdef PASO_MPI
  /* TODO */
  Finley_setError(SYSTEM_ERROR,"Finley_Mesh_read: MPI is not suporrted yet.");`
  return NULL;
#endif

  /* get file handle */
  fileHandle_p = fopen(fname, "r");
  if (fileHandle_p==NULL) {
    sprintf(error_msg,"Finley_Mesh_read: Opening file %s for reading failed.",fname);
    Finley_setError(IO_ERROR,error_msg);
    return NULL;
  }

  /* read header */
  sprintf(frm,"%%%d[^\n]",LenString_MAX-1);
  fscanf(fileHandle_p, frm, name);

  /* get the nodes */

  fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);
  /* allocate mesh */
  mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order);
  if (! Finley_noError()) return NULL;

  Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
  if (! Finley_noError()) return NULL;

  if (1 == numDim) {
      for (i0 = 0; i0 < numNodes; i0++)
    fscanf(fileHandle_p, "%d %d %d %le\n", &mesh_p->Nodes->Id[i0],
           &mesh_p->Nodes->degreeOfFreedom[i0], &mesh_p->Nodes->Tag[i0],
           &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)]);
  } else if (2 == numDim) {
      for (i0 = 0; i0 < numNodes; i0++)
    fscanf(fileHandle_p, "%d %d %d %le %le\n", &mesh_p->Nodes->Id[i0],
           &mesh_p->Nodes->degreeOfFreedom[i0], &mesh_p->Nodes->Tag[i0],
           &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],
           &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)]);
  } else if (3 == numDim) {
      for (i0 = 0; i0 < numNodes; i0++)
    fscanf(fileHandle_p, "%d %d %d %le %le %le\n", &mesh_p->Nodes->Id[i0],
           &mesh_p->Nodes->degreeOfFreedom[i0], &mesh_p->Nodes->Tag[i0],
           &mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)],
           &mesh_p->Nodes->Coordinates[INDEX2(1,i0,numDim)],
           &mesh_p->Nodes->Coordinates[INDEX2(2,i0,numDim)]);
  } /* if else else */

  /* get the element type */

  fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
  typeID=Finley_RefElement_getTypeId(element_type);
  if (typeID==NoType) {
    sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s",element_type);
    Finley_setError(VALUE_ERROR,error_msg);
    return NULL;
  }
  /* read the elements */
  mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order);
  Finley_ElementFile_allocTable(mesh_p->Elements, numEle);
  mesh_p->Elements->minColor=0;
  mesh_p->Elements->maxColor=numEle-1;
  for (i0 = 0; i0 < numEle; i0++) {
    fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);
    mesh_p->Elements->Color[i0]=i0;
    for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {
         fscanf(fileHandle_p, " %d",
            &mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);
    }   /* for i1 */
    fscanf(fileHandle_p, "\n");
  } /* for i0 */

  /* get the face elements */
  fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
  faceTypeID=Finley_RefElement_getTypeId(element_type);
  faceTypeID=Finley_RefElement_getTypeId(element_type);
  if (faceTypeID==NoType) {
    sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s for face elements",element_type);
    Finley_setError(VALUE_ERROR,error_msg);
    return NULL;
  }
  mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order, mesh_p->reduced_order);
  Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);
  mesh_p->FaceElements->minColor=0;
  mesh_p->FaceElements->maxColor=numEle-1;
  for (i0 = 0; i0 < numEle; i0++) {
    fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);
    mesh_p->FaceElements->Color[i0]=i0;
    for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {
         fscanf(fileHandle_p, " %d",
            &mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);
    }   /* for i1 */
    fscanf(fileHandle_p, "\n");
  } /* for i0 */

  /* get the Contact face element */
  fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
  contactTypeID=Finley_RefElement_getTypeId(element_type);
  if (contactTypeID==NoType) {
    sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);
    Finley_setError(VALUE_ERROR,error_msg);
    return NULL;
  }
  mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order, mesh_p->reduced_order);
  Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);
  mesh_p->ContactElements->minColor=0;
  mesh_p->ContactElements->maxColor=numEle-1;
  for (i0 = 0; i0 < numEle; i0++) {
    fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);
    mesh_p->ContactElements->Color[i0]=i0;
    for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {
        fscanf(fileHandle_p, " %d",
           &mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);
    }   /* for i1 */
    fscanf(fileHandle_p, "\n");
  } /* for i0 */

  /* get the nodal element */
  fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
  pointTypeID=Finley_RefElement_getTypeId(element_type);
  if (pointTypeID==NoType) {
    sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type);
    Finley_setError(VALUE_ERROR,error_msg);
    return NULL;
  }
  mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order, mesh_p->reduced_order);
  Finley_ElementFile_allocTable(mesh_p->Points, numEle);
  mesh_p->Points->minColor=0;
  mesh_p->Points->maxColor=numEle-1;
  for (i0 = 0; i0 < numEle; i0++) {
    fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);
    mesh_p->Points->Color[i0]=i0;
    for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {
        fscanf(fileHandle_p, " %d",
           &mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);
    }   /* for i1 */
    fscanf(fileHandle_p, "\n");
  } /* for i0 */
  /* get the name tags */
  if (feof(fileHandle_p) == 0) {
     fscanf(fileHandle_p, "%s\n", name);
     while (feof(fileHandle_p) == 0) {
       fscanf(fileHandle_p, "%s %d\n", name, &tag_key);
       Finley_Mesh_addTagMap(mesh_p,name,tag_key);
     }
  }
  /* close file */

  fclose(fileHandle_p);

  /*   resolve id's : */

  if (Finley_noError()) {
     Finley_Mesh_resolveNodeIds(mesh_p);
  }

  /* rearrange elements: */

  if (Finley_noError()) {
     Finley_Mesh_prepare(mesh_p);
  }

  /* optimize node labeling*/

  if (Finley_noError()) {
      if (optimize_labeling) Finley_Mesh_optimizeNodeLabeling(mesh_p);
  }

  /* that's it */
  #ifdef Finley_TRACE
  printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
  #endif
  if (Finley_noError()) {
       if ( ! Finley_Mesh_isPrepared(mesh_p)) {
          Finley_setError(SYSTEM_ERROR,"Mesh is not prepared for calculation. Contact the programmers.");
       }
  } else {
       Finley_Mesh_dealloc(mesh_p);
  }
  return mesh_p;
}
1	jgs	150	/*
2	elspeth	616	************************************************************
3			* Copyright 2006 by ACcESS MNRF *
4			* *
5			* http://www.access.edu.au *
6			* Primary Business: Queensland, Australia *
7			* Licensed under the Open Software License version 3.0 *
8			* http://www.opensource.org/licenses/osl-3.0.php *
9			* *
10			************************************************************
11	jgs	150	*/
12
13	jgs	82	/**************************************************************/
14
15			/* Finley: read mesh */
16
17			/**************************************************************/
18
19			/* Author: gross@access.edu.au */
20			/* Version: $Id$ */
21
22			/**************************************************************/
23
24			#include "Mesh.h"
25
26			/**************************************************************/
27
28			/* reads a mesh from a Finley file of name fname */
29
30	gross	1062	Finley_Mesh* Finley_Mesh_read(char* fname,index_t order, index_t reduced_order, bool_t optimize_labeling) {
31	jgs	82
32	jgs	123	dim_t numNodes, numDim, numEle, i0, i1;
33	gross	1044	index_t tag_key;
34	bcumming	730	Finley_Mesh *mesh_p=NULL;
35	jgs	150	char name[LenString_MAX],element_type[LenString_MAX],frm[20];
36			char error_msg[LenErrorMsg_MAX];
37	jgs	82	double time0=Finley_timer();
38	gross	1028	FILE *fileHandle_p = NULL;
39			ElementTypeId typeID, faceTypeID, contactTypeID, pointTypeID;
40	bcumming	730
41	jgs	150	Finley_resetError();
42	gross	1044	#ifdef PASO_MPI
43			/* TODO */
44			Finley_setError(SYSTEM_ERROR,"Finley_Mesh_read: MPI is not suporrted yet.");`
45	gross	1156	return NULL;
46	gross	1044	#endif
47	jgs	82
48			/* get file handle */
49	gross	1028	fileHandle_p = fopen(fname, "r");
50	jgs	82	if (fileHandle_p==NULL) {
51	gross	1028	sprintf(error_msg,"Finley_Mesh_read: Opening file %s for reading failed.",fname);
52	jgs	150	Finley_setError(IO_ERROR,error_msg);
53	jgs	82	return NULL;
54			}
55
56			/* read header */
57	jgs	150	sprintf(frm,"%%%d[^\n]",LenString_MAX-1);
58			fscanf(fileHandle_p, frm, name);
59	gross	934
60	gross	939	/* get the nodes */
61
62			fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);
63	gross	934	/* allocate mesh */
64	gross	1062	mesh_p = Finley_Mesh_alloc(name,numDim,order,reduced_order);
65	gross	939	if (! Finley_noError()) return NULL;
66	jgs	82
67	gross	934	Finley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
68	jgs	150	if (! Finley_noError()) return NULL;
69	jgs	82
70			if (1 == numDim) {
71			for (i0 = 0; i0 < numNodes; i0++)
72	gross	934	fscanf(fileHandle_p, "%d %d %d %le\n", &mesh_p->Nodes->Id[i0],
73			&mesh_p->Nodes->degreeOfFreedom[i0], &mesh_p->Nodes->Tag[i0],
74			&mesh_p->Nodes->Coordinates[INDEX2(0,i0,numDim)]);
75	jgs	82	} else if (2 == numDim) {
76			for (i0 = 0; i0 < numNodes; i0++)
77	gross	934	fscanf(fileHandle_p, "%d %d %d %le %le\n", &mesh_p->Nodes->Id[i0],
78			&mesh_p->Nodes->degreeOfFreedom[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	jgs	82	} else if (3 == numDim) {
82			for (i0 = 0; i0 < numNodes; i0++)
83	gross	934	fscanf(fileHandle_p, "%d %d %d %le %le %le\n", &mesh_p->Nodes->Id[i0],
84			&mesh_p->Nodes->degreeOfFreedom[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	jgs	82	} /* if else else */
89
90			/* get the element type */
91
92			fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
93	gross	1028	typeID=Finley_RefElement_getTypeId(element_type);
94	jgs	82	if (typeID==NoType) {
95	gross	1028	sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s",element_type);
96	jgs	150	Finley_setError(VALUE_ERROR,error_msg);
97	jgs	82	return NULL;
98			}
99			/* read the elements */
100	gross	1062	mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order, mesh_p->reduced_order);
101	jgs	82	Finley_ElementFile_allocTable(mesh_p->Elements, numEle);
102	jgs	123	mesh_p->Elements->minColor=0;
103			mesh_p->Elements->maxColor=numEle-1;
104	jgs	82	for (i0 = 0; i0 < numEle; i0++) {
105			fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);
106			mesh_p->Elements->Color[i0]=i0;
107			for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {
108	jgs	126	fscanf(fileHandle_p, " %d",
109	jgs	82	&mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);
110			} /* for i1 */
111			fscanf(fileHandle_p, "\n");
112			} /* for i0 */
113
114			/* get the face elements */
115			fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
116	gross	1028	faceTypeID=Finley_RefElement_getTypeId(element_type);
117			faceTypeID=Finley_RefElement_getTypeId(element_type);
118	jgs	82	if (faceTypeID==NoType) {
119	gross	1028	sprintf(error_msg,"Finley_Mesh_read :Unidentified element type %s for face elements",element_type);
120	jgs	150	Finley_setError(VALUE_ERROR,error_msg);
121	jgs	82	return NULL;
122			}
123	gross	1062	mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order, mesh_p->reduced_order);
124	jgs	82	Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);
125	jgs	123	mesh_p->FaceElements->minColor=0;
126			mesh_p->FaceElements->maxColor=numEle-1;
127	jgs	82	for (i0 = 0; i0 < numEle; i0++) {
128			fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);
129			mesh_p->FaceElements->Color[i0]=i0;
130			for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {
131	jgs	126	fscanf(fileHandle_p, " %d",
132	jgs	82	&mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);
133			} /* for i1 */
134			fscanf(fileHandle_p, "\n");
135			} /* for i0 */
136
137			/* get the Contact face element */
138			fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
139	gross	1028	contactTypeID=Finley_RefElement_getTypeId(element_type);
140	jgs	82	if (contactTypeID==NoType) {
141	gross	1028	sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for contact elements",element_type);
142	jgs	150	Finley_setError(VALUE_ERROR,error_msg);
143	jgs	82	return NULL;
144			}
145	gross	1062	mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order, mesh_p->reduced_order);
146	jgs	82	Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);
147	jgs	123	mesh_p->ContactElements->minColor=0;
148			mesh_p->ContactElements->maxColor=numEle-1;
149	jgs	82	for (i0 = 0; i0 < numEle; i0++) {
150			fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);
151			mesh_p->ContactElements->Color[i0]=i0;
152			for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {
153	jgs	126	fscanf(fileHandle_p, " %d",
154	jgs	82	&mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);
155			} /* for i1 */
156			fscanf(fileHandle_p, "\n");
157			} /* for i0 */
158
159			/* get the nodal element */
160			fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
161	gross	1028	pointTypeID=Finley_RefElement_getTypeId(element_type);
162	jgs	82	if (pointTypeID==NoType) {
163	gross	1028	sprintf(error_msg,"Finley_Mesh_read: Unidentified element type %s for points",element_type);
164	jgs	150	Finley_setError(VALUE_ERROR,error_msg);
165	jgs	82	return NULL;
166			}
167	gross	1062	mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order, mesh_p->reduced_order);
168	jgs	82	Finley_ElementFile_allocTable(mesh_p->Points, numEle);
169	jgs	123	mesh_p->Points->minColor=0;
170			mesh_p->Points->maxColor=numEle-1;
171	jgs	82	for (i0 = 0; i0 < numEle; i0++) {
172			fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);
173			mesh_p->Points->Color[i0]=i0;
174			for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {
175	jgs	126	fscanf(fileHandle_p, " %d",
176	jgs	82	&mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);
177			} /* for i1 */
178			fscanf(fileHandle_p, "\n");
179			} /* for i0 */
180	gross	1044	/* get the name tags */
181			if (feof(fileHandle_p) == 0) {
182			fscanf(fileHandle_p, "%s\n", name);
183			while (feof(fileHandle_p) == 0) {
184			fscanf(fileHandle_p, "%s %d\n", name, &tag_key);
185			Finley_Mesh_addTagMap(mesh_p,name,tag_key);
186			}
187			}
188	jgs	82	/* close file */
189
190			fclose(fileHandle_p);
191
192			/* resolve id's : */
193	jgs	126
194	gross	1156	if (Finley_noError()) {
195			Finley_Mesh_resolveNodeIds(mesh_p);
196			}
197	jgs	126
198	jgs	82	/* rearrange elements: */
199	jgs	126
200	gross	1156	if (Finley_noError()) {
201			Finley_Mesh_prepare(mesh_p);
202			}
203	jgs	126
204	gross	1156	/* optimize node labeling*/
205
206			if (Finley_noError()) {
207			if (optimize_labeling) Finley_Mesh_optimizeNodeLabeling(mesh_p);
208			}
209
210	jgs	82	/* that's it */
211	jgs	149	#ifdef Finley_TRACE
212	jgs	82	printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
213	jgs	149	#endif
214	gross	964	if (Finley_noError()) {
215			if ( ! Finley_Mesh_isPrepared(mesh_p)) {
216			Finley_setError(SYSTEM_ERROR,"Mesh is not prepared for calculation. Contact the programmers.");
217			}
218			} else {
219			Finley_Mesh_dealloc(mesh_p);
220			}
221	jgs	82	return mesh_p;
222			}
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision