/[escript]/branches/trilinos_from_5897/dudley/src/Mesh_readGmsh.cpp
ViewVC logotype

Contents of /branches/trilinos_from_5897/dudley/src/Mesh_readGmsh.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6009 - (show annotations)
Wed Mar 2 04:13:26 2016 UTC (3 years, 1 month ago) by caltinay
File size: 15388 byte(s)
Much needed sync with trunk...

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2016 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17 #include "Mesh.h"
18
19 #include <cstdio>
20
21 namespace dudley {
22
23 #define MAX_numNodes_gmsh 20
24
25 /* reads a mesh from a Dudley file of name fname */
26 Dudley_Mesh *Dudley_Mesh_readGmsh(char *fname, index_t numDim, index_t order, index_t reduced_order, bool optimize)
27 {
28 double version = 1.0;
29 int format = 0, size = sizeof(double), scan_ret;
30 dim_t numNodes, totalNumElements = 0, numTags = 0, numNodesPerElement = 0, numNodesPerElement2, element_dim = 0;
31 index_t e, i0, j, gmsh_type, partition_id, itmp, elementary_id;
32 index_t numElements = 0, numFaceElements = 0, *id = NULL, *tag = NULL, *vertices = NULL;
33 Dudley_Mesh *mesh_p = NULL;
34 char line[1024];
35 char error_msg[1024];
36 double rtmp0, rtmp1;
37 FILE *fileHandle_p = NULL;
38 Dudley_ElementTypeId *element_type = NULL;
39
40 /* No! Bad! take a parameter for this */
41 escript::JMPI mpi_info = escript::makeInfo(MPI_COMM_WORLD);
42 if (mpi_info->size > 1)
43 throw DudleyException("reading GMSH with MPI is not supported yet.");
44
45 /* allocate mesh */
46 mesh_p = Dudley_Mesh_alloc(fname, numDim, mpi_info);
47
48 /* get file handle */
49 fileHandle_p = fopen(fname, "r");
50 if (fileHandle_p == NULL) {
51 sprintf(error_msg, "Opening Gmsh file %s for reading failed.", fname);
52 throw DudleyException(error_msg);
53 }
54
55 /* start reading */
56 while (1) {
57 /* find line staring with $ */
58 do {
59 if (!fgets(line, sizeof(line), fileHandle_p))
60 break;
61 if (feof(fileHandle_p))
62 break;
63 }
64 while (line[0] != '$');
65
66 if (feof(fileHandle_p))
67 break;
68
69 /* format */
70 if (!strncmp(&line[1], "MeshFormat", 10)) {
71 scan_ret = fscanf(fileHandle_p, "%lf %d %d\n", &version, &format, &size);
72 if (scan_ret == EOF)
73 throw DudleyException("readGmsh: early EOF while reading file");
74 }
75 /* nodes are read */
76 if (!strncmp(&line[1], "NOD", 3) || !strncmp(&line[1], "NOE", 3) || !strncmp(&line[1], "Nodes", 5)) {
77
78 scan_ret = fscanf(fileHandle_p, "%d", &numNodes);
79 if (scan_ret == EOF)
80 throw DudleyException("readGmsh: early EOF while reading file");
81 Dudley_NodeFile_allocTable(mesh_p->Nodes, numNodes);
82 for (i0 = 0; i0 < numNodes; i0++)
83 {
84 if (1 == numDim)
85 {
86 scan_ret = fscanf(fileHandle_p, "%d %le %le %le\n", &mesh_p->Nodes->Id[i0],
87 &mesh_p->Nodes->Coordinates[INDEX2(0, i0, numDim)], &rtmp0, &rtmp1);
88 if (scan_ret == EOF)
89 throw DudleyException("readGmsh: early EOF while reading file");
90 }
91 else if (2 == numDim)
92 {
93 scan_ret = fscanf(fileHandle_p, "%d %le %le %le\n", &mesh_p->Nodes->Id[i0],
94 &mesh_p->Nodes->Coordinates[INDEX2(0, i0, numDim)],
95 &mesh_p->Nodes->Coordinates[INDEX2(1, i0, numDim)], &rtmp0);
96 if (scan_ret == EOF)
97 throw DudleyException("readGmsh: early EOF while reading file");
98 }
99 else if (3 == numDim)
100 {
101 scan_ret = fscanf(fileHandle_p, "%d %le %le %le\n", &mesh_p->Nodes->Id[i0],
102 &mesh_p->Nodes->Coordinates[INDEX2(0, i0, numDim)],
103 &mesh_p->Nodes->Coordinates[INDEX2(1, i0, numDim)],
104 &mesh_p->Nodes->Coordinates[INDEX2(2, i0, numDim)]);
105 if (scan_ret == EOF)
106 throw DudleyException("readGmsh: early EOF while reading file");
107 }
108 mesh_p->Nodes->globalDegreesOfFreedom[i0] = mesh_p->Nodes->Id[i0];
109 mesh_p->Nodes->Tag[i0] = 0;
110 }
111 }
112 /* elements */
113 else if (!strncmp(&line[1], "ELM", 3) || !strncmp(&line[1], "Elements", 8))
114 {
115 Dudley_ElementTypeId final_element_type = Dudley_NoRef;
116 Dudley_ElementTypeId final_face_element_type = Dudley_NoRef;
117 numElements = 0;
118 numFaceElements = 0;
119 scan_ret = fscanf(fileHandle_p, "%d", &totalNumElements);
120 if (scan_ret == EOF)
121 throw DudleyException("readGmsh: early EOF while reading file");
122
123 id = new index_t[totalNumElements];
124 tag = new index_t[totalNumElements];
125
126 element_type = new Dudley_ElementTypeId[totalNumElements];
127 vertices = new index_t[totalNumElements * MAX_numNodes_gmsh];
128 /* read all in */
129 for (e = 0; e < totalNumElements; e++)
130 {
131 scan_ret = fscanf(fileHandle_p, "%d %d", &id[e], &gmsh_type);
132 if (scan_ret == EOF)
133 throw DudleyException("readGmsh: early EOF while reading file");
134 switch (gmsh_type)
135 {
136 case 1: /* line order 1 */
137 element_type[e] = Dudley_Line2;
138 element_dim = 1;
139 numNodesPerElement = 2;
140 break;
141 case 2: /* triangle order 1 */
142 element_type[e] = Dudley_Tri3;
143 numNodesPerElement = 3;
144 element_dim = 2;
145 break;
146 case 4: /* tetrahedron order 1 */
147 element_type[e] = Dudley_Tet4;
148 numNodesPerElement = 4;
149 element_dim = 3;
150 break;
151 case 15: /* point */
152 element_type[e] = Dudley_Point1;
153 numNodesPerElement = 1;
154 element_dim = 0;
155 break;
156 default:
157 element_type[e] = Dudley_NoRef;
158 sprintf(error_msg, "Unexpected gmsh element type %d in mesh file %s.", gmsh_type, fname);
159 throw DudleyException(error_msg);
160 }
161 if (element_dim == numDim)
162 {
163 if (final_element_type == Dudley_NoRef)
164 {
165 final_element_type = element_type[e];
166 }
167 else if (final_element_type != element_type[e])
168 {
169 sprintf(error_msg, "Dudley can handle a single type of internal elements only.");
170 throw DudleyException(error_msg);
171 }
172 numElements++;
173 }
174 else if (element_dim == numDim - 1)
175 {
176 if (final_face_element_type == Dudley_NoRef) {
177 final_face_element_type = element_type[e];
178 } else if (final_face_element_type != element_type[e]) {
179 sprintf(error_msg, "Dudley can handle a single type of face elements only.");
180 throw DudleyException(error_msg);
181 }
182 numFaceElements++;
183 }
184
185 if (version <= 1.0) {
186 scan_ret = fscanf(fileHandle_p, "%d %d %d", &tag[e], &elementary_id, &numNodesPerElement2);
187 if (scan_ret == EOF)
188 throw DudleyException("readGmsh: early EOF while reading file");
189 partition_id = 1;
190 if (numNodesPerElement2 != numNodesPerElement) {
191 sprintf(error_msg, "Illegal number of nodes for element %d in mesh file %s.", id[e],
192 fname);
193 throw DudleyException(error_msg);
194 }
195 } else {
196 scan_ret = fscanf(fileHandle_p, "%d", &numTags);
197 if (scan_ret == EOF)
198 throw DudleyException("readGmsh: early EOF while reading file");
199 elementary_id = tag[e] = partition_id = 1;
200 numNodesPerElement2 = -1;
201 for (j = 0; j < numTags; j++) {
202 scan_ret = fscanf(fileHandle_p, "%d", &itmp);
203 if (scan_ret == EOF)
204 throw DudleyException("readGmsh: early EOF while reading file");
205 if (j == 0) {
206 tag[e] = itmp;
207 } else if (j == 1) {
208 elementary_id = itmp;
209 } else if (j == 2) {
210 partition_id = itmp;
211 }
212 /* ignore any other tags */
213 }
214 }
215 for (j = 0; j < numNodesPerElement; j++)
216 {
217 scan_ret = fscanf(fileHandle_p, "%d", &vertices[INDEX2(j, e, MAX_numNodes_gmsh)]);
218 if (scan_ret == EOF)
219 throw DudleyException("readGmsh: early EOF while reading file");
220 }
221 }
222 /* all elements have been read, now we have to identify the elements for dudley */
223
224 /* first we have to identify the elements to define Elements and FaceElements */
225 if (final_element_type == Dudley_NoRef) {
226 if (numDim == 1) {
227 final_element_type = Dudley_Line2;
228 } else if (numDim == 2) {
229 final_element_type = Dudley_Tri3;
230 } else if (numDim == 3) {
231 final_element_type = Dudley_Tet4;
232 }
233 }
234 if (final_face_element_type == Dudley_NoRef) {
235 if (numDim == 1) {
236 final_face_element_type = Dudley_Point1;
237 } else if (numDim == 2) {
238 final_face_element_type = Dudley_Line2;
239 } else if (numDim == 3) {
240 final_face_element_type = Dudley_Tri3;
241 }
242 }
243 mesh_p->Elements = Dudley_ElementFile_alloc(final_element_type, mpi_info);
244 mesh_p->FaceElements = Dudley_ElementFile_alloc(final_face_element_type, mpi_info);
245 mesh_p->Points = Dudley_ElementFile_alloc(Dudley_Point1, mpi_info);
246 Dudley_ElementFile_allocTable(mesh_p->Elements, numElements);
247 Dudley_ElementFile_allocTable(mesh_p->FaceElements, numFaceElements);
248 Dudley_ElementFile_allocTable(mesh_p->Points, 0);
249 mesh_p->Elements->minColor = 0;
250 mesh_p->Elements->maxColor = numElements - 1;
251 mesh_p->FaceElements->minColor = 0;
252 mesh_p->FaceElements->maxColor = numFaceElements - 1;
253 mesh_p->Points->minColor = 0;
254 mesh_p->Points->maxColor = 0;
255 numElements = 0;
256 numFaceElements = 0;
257 for (e = 0; e < totalNumElements; e++)
258 {
259 if (element_type[e] == final_element_type)
260 {
261 mesh_p->Elements->Id[numElements] = id[e];
262 mesh_p->Elements->Tag[numElements] = tag[e];
263 mesh_p->Elements->Color[numElements] = numElements;
264 mesh_p->Elements->Owner[numElements] = 0;
265 for (j = 0; j < mesh_p->Elements-> /*referenceElementSet-> */ numNodes; ++j)
266 {
267 mesh_p->Elements->Nodes[INDEX2
268 (j, numElements,
269 mesh_p->
270 Elements-> /*referenceElementSet-> */ numNodes)] =
271 vertices[INDEX2(j, e, MAX_numNodes_gmsh)];
272 }
273 numElements++;
274 }
275 else if (element_type[e] == final_face_element_type)
276 {
277 mesh_p->FaceElements->Id[numFaceElements] = id[e];
278 mesh_p->FaceElements->Tag[numFaceElements] = tag[e];
279 mesh_p->FaceElements->Color[numFaceElements] = numFaceElements;
280 mesh_p->FaceElements->Owner[numFaceElements] = 0;
281 for (j = 0; j < mesh_p->FaceElements-> /*referenceElementSet-> */ numNodes; ++j)
282 {
283 mesh_p->FaceElements->Nodes[INDEX2
284 (j, numFaceElements,
285 mesh_p->
286 FaceElements-> /*referenceElementSet-> */
287 numNodes)] =
288 vertices[INDEX2(j, e, MAX_numNodes_gmsh)];
289 }
290 numFaceElements++;
291 }
292 }
293 /* and clean up */
294 delete[] id;
295 delete[] tag;
296 delete[] element_type;
297 delete[] vertices;
298 }
299 /* name tags (thanks to Antoine Lefebvre, antoine.lefebvre2@mail.mcgill.ca ) */
300 else if (!strncmp(&line[1], "PhysicalNames", 13)) {
301 char name[1024];
302 index_t tag_key;
303 scan_ret = fscanf(fileHandle_p, "%d", &numTags);
304 if (scan_ret == EOF)
305 throw DudleyException("readGmsh: early EOF while reading file");
306 for (i0 = 0; i0 < numTags; i0++) {
307 scan_ret = fscanf(fileHandle_p, "%d %d %s\n", &itmp, &tag_key, name);
308 if (scan_ret == EOF)
309 throw DudleyException("readGmsh: early EOF while reading file");
310 if (itmp != 2)
311 throw DudleyException("Dudley_Mesh_readGmsh: expecting two entries per physical name.");
312 if (strlen(name) < 3)
313 throw DudleyException("Dudley_Mesh_readGmsh: illegal tagname (\" missing?)");
314 name[strlen(name)-1]='\0';
315 Dudley_Mesh_addTagMap(mesh_p,&name[1],tag_key);
316 }
317 }
318 /* search for end of data block */
319 do
320 {
321 if (!fgets(line, sizeof(line), fileHandle_p)) {
322 sprintf(error_msg, "Unexpected end of file in %s", fname);
323 throw DudleyException(error_msg);
324 }
325 if (feof(fileHandle_p)) {
326 sprintf(error_msg, "Unexpected end of file in %s", fname);
327 throw DudleyException(error_msg);
328 }
329 }
330 while (line[0] != '$');
331 }
332
333 fclose(fileHandle_p);
334 Dudley_Mesh_resolveNodeIds(mesh_p);
335 Dudley_Mesh_prepare(mesh_p, optimize);
336 return mesh_p;
337 }
338
339 } // namespace dudley
340

Properties

Name Value
svn:mergeinfo /branches/4.0fordebian/dudley/src/Mesh_readGmsh.cpp:5567-5588 /branches/lapack2681/finley/src/Mesh_readGmsh.cpp:2682-2741 /branches/pasowrap/dudley/src/Mesh_readGmsh.cpp:3661-3674 /branches/py3_attempt2/dudley/src/Mesh_readGmsh.cpp:3871-3891 /branches/restext/finley/src/Mesh_readGmsh.cpp:2610-2624 /branches/ripleygmg_from_3668/dudley/src/Mesh_readGmsh.cpp:3669-3791 /branches/stage3.0/finley/src/Mesh_readGmsh.cpp:2569-2590 /branches/symbolic_from_3470/dudley/src/Mesh_readGmsh.cpp:3471-3974 /branches/symbolic_from_3470/ripley/test/python/dudley/src/Mesh_readGmsh.cpp:3517-3974 /release/3.0/finley/src/Mesh_readGmsh.cpp:2591-2601 /release/4.0/dudley/src/Mesh_readGmsh.cpp:5380-5406 /trunk/dudley/src/Mesh_readGmsh.cpp:4257-4344,5898-6007 /trunk/ripley/test/python/dudley/src/Mesh_readGmsh.cpp:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26