src/finleyC/Mesh_read.c

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

/*   Finley: read mesh */

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

/*   Copyrights by ACcESS Australia 2003/04 */
/*   Author: gross@access.edu.au */
/*   Version: $Id$ */

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

#include "Finley.h"
#include "Mesh.h"

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

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

Finley_Mesh* Finley_Mesh_read(char* fname,index_t order) {

  dim_t numNodes, numDim, numEle, i0, i1;
  char name[LenString_MAX],element_type[LenString_MAX];
  Finley_NodeFile *nodes_p=NULL;
  double time0=Finley_timer();

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

  /* read header */
  fscanf(fileHandle_p, "%63[^\n]", name);

  /* get the nodes */

  fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);
  nodes_p=Finley_NodeFile_alloc(numDim);
  if (Finley_ErrorCode!=NO_ERROR) return NULL;
  Finley_NodeFile_allocTable(nodes_p, numNodes);
  if (Finley_ErrorCode!=NO_ERROR) return NULL;

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

  /* get the element type */

  fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
  ElementTypeId typeID=Finley_RefElement_getTypeId(element_type);
  if (typeID==NoType) {
    Finley_ErrorCode=VALUE_ERROR;
    sprintf(Finley_ErrorMsg,"Unidentified element type %s",element_type);
    return NULL;
  }

  /* allocate mesh */

  Finley_Mesh * mesh_p =Finley_Mesh_alloc(name,numDim,order);
  if (Finley_ErrorCode!=NO_ERROR) return NULL;
  mesh_p->Nodes=nodes_p;

  /* read the elements */
  mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->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);
  ElementTypeId faceTypeID=Finley_RefElement_getTypeId(element_type);
  if (faceTypeID==NoType) {
    Finley_ErrorCode=VALUE_ERROR;
    sprintf(Finley_ErrorMsg,"Unidentified element type %s for face elements",element_type);
    return NULL;
  }
  mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->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);
  ElementTypeId contactTypeID=Finley_RefElement_getTypeId(element_type);
  if (contactTypeID==NoType) {
    Finley_ErrorCode=VALUE_ERROR;
    sprintf(Finley_ErrorMsg,"Unidentified element type %s for contact elements",element_type);
    return NULL;
  }
  mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->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);
  ElementTypeId pointTypeID=Finley_RefElement_getTypeId(element_type);
  if (pointTypeID==NoType) {
    Finley_ErrorCode=VALUE_ERROR;
    sprintf(Finley_ErrorMsg,"Unidentified element type %s for points",element_type);
    return NULL;
  }
  mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->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 */


  /* close file */

  fclose(fileHandle_p);

  /*   resolve id's : */

  Finley_Mesh_resolveNodeIds(mesh_p);

  /* rearrange elements: */

  Finley_Mesh_prepare(mesh_p);
  printf ("nodes read!\n");

  /* that's it */
  printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
  if (Finley_ErrorCode!=NO_ERROR) Finley_Mesh_dealloc(mesh_p);
  return mesh_p;
}
/*
* $Log$
* Revision 1.3  2005/07/22 03:53:08  jgs
* Merge of development branch back to main trunk on 2005-07-22
*
* Revision 1.2  2005/07/08 04:07:54  jgs
* Merge of development branch back to main trunk on 2005-07-08
*
* Revision 1.1.1.1.2.2  2005/07/18 10:34:54  gross
* some informance improvements when reading meshes
*
* Revision 1.1.1.1.2.1  2005/06/29 02:34:53  gross
* some changes towards 64 integers in finley
*
* Revision 1.1.1.1  2004/10/26 06:53:57  jgs
* initial import of project esys2
*
* Revision 1.1.1.1  2004/06/24 04:00:40  johng
* Initial version of eys using boost-python.
*
*
*/

1	jgs	82	/**************************************************************/
2
3			/* Finley: read mesh */
4
5			/**************************************************************/
6
7			/* Copyrights by ACcESS Australia 2003/04 */
8			/* Author: gross@access.edu.au */
9			/* Version: $Id$ */
10
11			/**************************************************************/
12
13			#include "Finley.h"
14			#include "Mesh.h"
15
16			/**************************************************************/
17
18			/* reads a mesh from a Finley file of name fname */
19
20	jgs	123	Finley_Mesh* Finley_Mesh_read(char* fname,index_t order) {
21	jgs	82
22	jgs	123	dim_t numNodes, numDim, numEle, i0, i1;
23	jgs	82	char name[LenString_MAX],element_type[LenString_MAX];
24			Finley_NodeFile *nodes_p=NULL;
25			double time0=Finley_timer();
26
27			/* get file handle */
28			FILE * fileHandle_p = fopen(fname, "r");
29			if (fileHandle_p==NULL) {
30			Finley_ErrorCode=IO_ERROR;
31			sprintf(Finley_ErrorMsg,"Opening file %s for reading failed.",fname);
32			return NULL;
33			}
34
35			/* read header */
36			fscanf(fileHandle_p, "%63[^\n]", name);
37
38			/* get the nodes */
39
40			fscanf(fileHandle_p, "%1d%*s %d\n", &numDim,&numNodes);
41			nodes_p=Finley_NodeFile_alloc(numDim);
42			if (Finley_ErrorCode!=NO_ERROR) return NULL;
43			Finley_NodeFile_allocTable(nodes_p, numNodes);
44			if (Finley_ErrorCode!=NO_ERROR) return NULL;
45
46			if (1 == numDim) {
47			for (i0 = 0; i0 < numNodes; i0++)
48			fscanf(fileHandle_p, "%d %d %d %le\n", &nodes_p->Id[i0],
49			&nodes_p->degreeOfFreedom[i0], &nodes_p->Tag[i0],
50			&nodes_p->Coordinates[INDEX2(0,i0,numDim)]);
51			} else if (2 == numDim) {
52			for (i0 = 0; i0 < numNodes; i0++)
53			fscanf(fileHandle_p, "%d %d %d %le %le\n", &nodes_p->Id[i0],
54			&nodes_p->degreeOfFreedom[i0], &nodes_p->Tag[i0],
55			&nodes_p->Coordinates[INDEX2(0,i0,numDim)],
56			&nodes_p->Coordinates[INDEX2(1,i0,numDim)]);
57			} else if (3 == numDim) {
58			for (i0 = 0; i0 < numNodes; i0++)
59			fscanf(fileHandle_p, "%d %d %d %le %le %le\n", &nodes_p->Id[i0],
60			&nodes_p->degreeOfFreedom[i0], &nodes_p->Tag[i0],
61			&nodes_p->Coordinates[INDEX2(0,i0,numDim)],
62			&nodes_p->Coordinates[INDEX2(1,i0,numDim)],
63			&nodes_p->Coordinates[INDEX2(2,i0,numDim)]);
64			} /* if else else */
65
66			/* get the element type */
67
68			fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
69			ElementTypeId typeID=Finley_RefElement_getTypeId(element_type);
70			if (typeID==NoType) {
71			Finley_ErrorCode=VALUE_ERROR;
72			sprintf(Finley_ErrorMsg,"Unidentified element type %s",element_type);
73			return NULL;
74			}
75
76			/* allocate mesh */
77
78			Finley_Mesh * mesh_p =Finley_Mesh_alloc(name,numDim,order);
79			if (Finley_ErrorCode!=NO_ERROR) return NULL;
80			mesh_p->Nodes=nodes_p;
81
82			/* read the elements */
83			mesh_p->Elements=Finley_ElementFile_alloc(typeID,mesh_p->order);
84			Finley_ElementFile_allocTable(mesh_p->Elements, numEle);
85	jgs	123	mesh_p->Elements->minColor=0;
86			mesh_p->Elements->maxColor=numEle-1;
87	jgs	82	for (i0 = 0; i0 < numEle; i0++) {
88			fscanf(fileHandle_p, "%d %d", &mesh_p->Elements->Id[i0], &mesh_p->Elements->Tag[i0]);
89			mesh_p->Elements->Color[i0]=i0;
90			for (i1 = 0; i1 < mesh_p->Elements->ReferenceElement->Type->numNodes; i1++) {
91	jgs	126	fscanf(fileHandle_p, " %d",
92	jgs	82	&mesh_p->Elements->Nodes[INDEX2(i1, i0, mesh_p->Elements->ReferenceElement->Type->numNodes)]);
93			} /* for i1 */
94			fscanf(fileHandle_p, "\n");
95			} /* for i0 */
96
97			/* get the face elements */
98			fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
99			ElementTypeId faceTypeID=Finley_RefElement_getTypeId(element_type);
100			if (faceTypeID==NoType) {
101			Finley_ErrorCode=VALUE_ERROR;
102			sprintf(Finley_ErrorMsg,"Unidentified element type %s for face elements",element_type);
103			return NULL;
104			}
105			mesh_p->FaceElements=Finley_ElementFile_alloc(faceTypeID,mesh_p->order);
106			Finley_ElementFile_allocTable(mesh_p->FaceElements, numEle);
107	jgs	123	mesh_p->FaceElements->minColor=0;
108			mesh_p->FaceElements->maxColor=numEle-1;
109	jgs	82	for (i0 = 0; i0 < numEle; i0++) {
110			fscanf(fileHandle_p, "%d %d", &mesh_p->FaceElements->Id[i0], &mesh_p->FaceElements->Tag[i0]);
111			mesh_p->FaceElements->Color[i0]=i0;
112			for (i1 = 0; i1 < mesh_p->FaceElements->ReferenceElement->Type->numNodes; i1++) {
113	jgs	126	fscanf(fileHandle_p, " %d",
114	jgs	82	&mesh_p->FaceElements->Nodes[INDEX2(i1, i0, mesh_p->FaceElements->ReferenceElement->Type->numNodes)]);
115			} /* for i1 */
116			fscanf(fileHandle_p, "\n");
117			} /* for i0 */
118
119			/* get the Contact face element */
120			fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
121			ElementTypeId contactTypeID=Finley_RefElement_getTypeId(element_type);
122			if (contactTypeID==NoType) {
123			Finley_ErrorCode=VALUE_ERROR;
124			sprintf(Finley_ErrorMsg,"Unidentified element type %s for contact elements",element_type);
125			return NULL;
126			}
127			mesh_p->ContactElements=Finley_ElementFile_alloc(contactTypeID,mesh_p->order);
128			Finley_ElementFile_allocTable(mesh_p->ContactElements, numEle);
129	jgs	123	mesh_p->ContactElements->minColor=0;
130			mesh_p->ContactElements->maxColor=numEle-1;
131	jgs	82	for (i0 = 0; i0 < numEle; i0++) {
132			fscanf(fileHandle_p, "%d %d", &mesh_p->ContactElements->Id[i0], &mesh_p->ContactElements->Tag[i0]);
133			mesh_p->ContactElements->Color[i0]=i0;
134			for (i1 = 0; i1 < mesh_p->ContactElements->ReferenceElement->Type->numNodes; i1++) {
135	jgs	126	fscanf(fileHandle_p, " %d",
136	jgs	82	&mesh_p->ContactElements->Nodes[INDEX2(i1, i0, mesh_p->ContactElements->ReferenceElement->Type->numNodes)]);
137			} /* for i1 */
138			fscanf(fileHandle_p, "\n");
139			} /* for i0 */
140
141			/* get the nodal element */
142			fscanf(fileHandle_p, "%s %d\n", element_type, &numEle);
143			ElementTypeId pointTypeID=Finley_RefElement_getTypeId(element_type);
144			if (pointTypeID==NoType) {
145			Finley_ErrorCode=VALUE_ERROR;
146			sprintf(Finley_ErrorMsg,"Unidentified element type %s for points",element_type);
147			return NULL;
148			}
149			mesh_p->Points=Finley_ElementFile_alloc(pointTypeID,mesh_p->order);
150
151			Finley_ElementFile_allocTable(mesh_p->Points, numEle);
152	jgs	123	mesh_p->Points->minColor=0;
153			mesh_p->Points->maxColor=numEle-1;
154	jgs	82	for (i0 = 0; i0 < numEle; i0++) {
155			fscanf(fileHandle_p, "%d %d", &mesh_p->Points->Id[i0], &mesh_p->Points->Tag[i0]);
156			mesh_p->Points->Color[i0]=i0;
157			for (i1 = 0; i1 < mesh_p->Points->ReferenceElement->Type->numNodes; i1++) {
158	jgs	126	fscanf(fileHandle_p, " %d",
159	jgs	82	&mesh_p->Points->Nodes[INDEX2(i1, i0, mesh_p->Points->ReferenceElement->Type->numNodes)]);
160			} /* for i1 */
161			fscanf(fileHandle_p, "\n");
162			} /* for i0 */
163
164
165			/* close file */
166
167			fclose(fileHandle_p);
168
169			/* resolve id's : */
170	jgs	126
171	jgs	82	Finley_Mesh_resolveNodeIds(mesh_p);
172	jgs	126
173	jgs	82	/* rearrange elements: */
174	jgs	126
175	jgs	82	Finley_Mesh_prepare(mesh_p);
176	jgs	126	printf ("nodes read!\n");
177
178	jgs	82	/* that's it */
179			printf("timing: reading mesh: %.4e sec\n",Finley_timer()-time0);
180			if (Finley_ErrorCode!=NO_ERROR) Finley_Mesh_dealloc(mesh_p);
181			return mesh_p;
182			}
183			/*
184			* $Log$
185	jgs	126	* Revision 1.3 2005/07/22 03:53:08 jgs
186			* Merge of development branch back to main trunk on 2005-07-22
187			*
188	jgs	123	* Revision 1.2 2005/07/08 04:07:54 jgs
189			* Merge of development branch back to main trunk on 2005-07-08
190	jgs	82	*
191	jgs	126	* Revision 1.1.1.1.2.2 2005/07/18 10:34:54 gross
192			* some informance improvements when reading meshes
193			*
194	jgs	123	* Revision 1.1.1.1.2.1 2005/06/29 02:34:53 gross
195			* some changes towards 64 integers in finley
196			*
197			* Revision 1.1.1.1 2004/10/26 06:53:57 jgs
198			* initial import of project esys2
199			*
200	jgs	82	* Revision 1.1.1.1 2004/06/24 04:00:40 johng
201			* Initial version of eys using boost-python.
202			*
203			*
204			*/
205
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision