1 /*
2 XLiFE++ is an extended library of finite elements written in C++
3     Copyright (C) 2014  Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin
4 
5     This program is free software: you can redistribute it and/or modify
6     it under the terms of the GNU General Public License as published by
7     the Free Software Foundation, either version 3 of the License, or
8     (at your option) any later version.
9     This program is distributed in the hope that it will be useful,
10     but WITHOUT ANY WARRANTY; without even the implied warranty of
11     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12     GNU General Public License for more details.
13     You should have received a copy of the GNU General Public License
14     along with this program.  If not, see <http://www.gnu.org/licenses/>.
15  */
16 
17 /*!
18   \file loadVtk.cpp
19   \author N. Kielbasiewicz
20   \since 5 may 2015
21   \date 6 may 2015
22 
23   \brief Read a file containing a mesh in VTK/VTU format
24 
25   The useful function is loadPly, which is a member function of the class Mesh.
26 */
27 
28 #include "../Mesh.hpp"
29 
30 namespace xlifepp {
31 
32 //! loads mesh from input stream according to POLYDATA dataset in vtk mesh format
loadPolyDataVtk(std::istream & data,number_t nodesDim)33 dimen_t Mesh::loadPolyDataVtk(std::istream& data, number_t nodesDim)
34 {
35   string_t strVal;
36   number_t numVal;
37 
38   data >> strVal;
39   if (strVal != "POINTS") { error("vtk_bad_block", words("First"), "POLYDATA", "POINTS"); }
40   number_t nbPoints;
41   data >> nbPoints;
42   data >> strVal; // data type of points (ignored)
43   // Definition of nodes
44   std::vector<std::vector<real_t> > coords(nbPoints,std::vector<real_t>(3));
45   nodes.resize(nbPoints);
46   vertices_.resize(nbPoints);
47   for (number_t i=0; i < nbPoints; ++i)
48   {
49     data >> coords[i][0] >> coords[i][1] >> coords[i][2];
50     vertices_[i]=i+1; // every node is a vertex
51   }
52   data >> strVal;
53   if (strVal != "POLYGONS") { error("vtk_bad_block", words("Second"), "POLYDATA", "POLYGONS"); }
54   number_t nbElems;
55   data >> nbElems;
56   data >> numVal; // total number of integers to read
57   std::vector<std::vector<number_t> > cells(nbElems);
58 
59   for (number_t i =0; i < nbElems; ++i)
60   {
61     number_t nbVertices;
62     data >> nbVertices;
63     cells[i].resize(nbVertices);
64     for (number_t j=0; j < nbVertices; ++j)
65     {
66       data >> numVal;
67       cells[i][j]=numVal+1;
68     }
69   }
70   dimen_t spaceDim=2;
71   if (nbPoints > 4)
72   {
73     // we take 4 vertices randomly and check their coplanarity
74     number_t v1=0, v2=(nbPoints-1)/3, v3=(2*(nbPoints-1))/3, v4=nbPoints-1;
75     if ( std::abs(dot(coords[v2]-coords[v1],crossProduct(coords[v3]-coords[v1],coords[v4]-coords[v1]))) > theEpsilon ) { spaceDim=3; }
76   }
77   if (nodesDim != 0 && nodesDim > spaceDim) { spaceDim=nodesDim; }
78 
79   nodes.resize(nbPoints,Point(std::vector<real_t>(spaceDim, 0.)));
80   for (number_t i=0; i < nbPoints; ++i)
81   {
82     if (spaceDim == 2) { nodes[i]=Point(coords[i][0], coords[i][1]); }
83     else { nodes[i]=Point(coords[i][0], coords[i][1], coords[i][2]); }
84   }
85 
86   elements_.resize(nbElems);
87   for (number_t i =0; i < nbElems; ++i)
88   {
89     number_t nbVertices = cells[i].size();
90     Interpolation* interp_p;
91     RefElement* re_p=0;
92     if (nbVertices == 3) // element is a triangle
93     {
94       interp_p=findInterpolation(_Lagrange, _standard, _P1, _H1);
95       re_p= findRefElement(_triangle, interp_p);
96     }
97     else if (nbVertices == 4) // element is a quadrangle
98     {
99       interp_p=findInterpolation(_Lagrange, _standard, _Q1, _H1);
100       re_p= findRefElement(_quadrangle, interp_p);
101     }
102     else { error("index_out_of_range", "nbVertices", 3, 4); }
103     // create GeomElement
104     elements_[i]=new GeomElement(this, re_p, spaceDim, i+1);
105     MeshElement * melt=elements_[i]->meshElement();
106     melt->vertexNumbers.resize(nbVertices);
107     melt->nodeNumbers.resize(nbVertices);
108     for (number_t j=0; j < nbVertices; ++j)
109     {
110       melt->vertexNumbers[j] = melt->nodeNumbers[j] = cells[i][j];
111     }
112     melt->setNodes(nodes);
113   }
114 
115   return 2;
116 }
117 
118 //! loads mesh from input stream according to RECTILINEAR_GRID dataset in vtk mesh format
loadRectilinearVtk(std::istream & data,number_t nodesDim)119 dimen_t Mesh::loadRectilinearVtk(std::istream& data, number_t nodesDim)
120 {
121   string_t strVal;
122   number_t numVal;
123 
124   data >> strVal;
125   if (strVal != "DIMENSIONS") { error("vtk_bad_block", words("First"), "RECTILINEAR_GRID", "DIMENSIONS"); }
126   number_t nx, ny, nz;
127   data >> nx >> ny >> nz;
128   std::vector<std::vector<real_t> > coords(nx*ny*nz,std::vector<real_t>(3));
129 
130   data >> strVal;
131   if (strVal != "X_COORDINATES") { error("vtk_bad_block", words("Second"), "RECTILINEAR_GRID", "X_COORDINATES"); }
132   data >> numVal;
133   if (numVal != nx) { error("bad_size", "", nx, numVal); }
134   data >> strVal; // data type of coordinates (ignored)
135   for (number_t i=0; i < nx; ++i) { data >> coords[i][0]; }
136 
137   data >> strVal;
138   if (strVal != "Y_COORDINATES") { error("vtk_bad_block", words("Third"), "RECTILINEAR_GRID", "Y_COORDINATES"); }
139   data >> numVal;
140   if (numVal != ny) { error("bad_size", "", ny, numVal); }
141   data >> strVal; // data type of coordinates (ignored)
142   for (number_t i=0; i < ny; ++i) { data >> coords[i][1]; }
143 
144   data >> strVal;
145   if (strVal != "Z_COORDINATES") { error("vtk_bad_block", words("Fourth"), "RECTILINEAR_GRID", "Z_COORDINATES"); }
146   data >> numVal;
147   if (numVal != nz) { error("bad_size", "", nz, numVal); }
148   data >> strVal; // data type of coordinates (ignored)
149   for (number_t i=0; i < nz; ++i) { data >> coords[i][2]; }
150 
151   // building vertices
152   vertices_.resize(nx*ny*nz);
153   for (number_t k=0; k < nz; ++k)
154   {
155     for (number_t j=0; j < ny; ++j)
156     {
157       for (number_t i=0; i < nx; ++i)
158       {
159         number_t l=k*nx*ny+j*nx+i;
160         vertices_[l]=l+1;
161       }
162     }
163   }
164   // if nx=1 and ny=1 -> segments along the z axis
165   // if nx=1 and nz=1 -> segments along the y axis
166   // if ny=1 and nz=1 -> segments along the x axis
167   // if nx=1 -> quadrangles along the yz plane
168   // if ny=1 -> quadrangles along the xz plane
169   // if nz=1 -> quadrangles along the xy plane
170   // else hexahedra
171 
172   Interpolation* interp_p;
173   if (ny*nz==1 || nx*nz==1 || nx*ny==1)
174   {
175     // segments are P1
176     interp_p=findInterpolation(_Lagrange, _standard, _P1, _H1);
177   }
178   else
179   {
180     // quadrangles and hexahedra are Q1
181     interp_p=findInterpolation(_Lagrange, _standard, _Q1, _H1);
182   }
183 
184   dimen_t spaceDim=3;
185   if (nx == 1 && ny == 1) { spaceDim=1; } // segments along the z axis
186   else if (nx == 1 && nz == 1) { spaceDim=1; } // segments along the y axis
187   else if (ny == 1 && nz == 1) { spaceDim=1; } // segments along the x axis
188   else if (nx == 1) { spaceDim=2; } // quadrangles along the yz plane
189   else if (ny == 1) { spaceDim=2; } // quadrangles along the xz plane
190   else if (nz == 1) { spaceDim=2; } // quadrangles along the xy plane
191   if (nodesDim != 0 && nodesDim > spaceDim) { spaceDim=nodesDim; }
192 
193   nodes.resize(nx*ny*nz,Point(std::vector<real_t>(spaceDim, 0.)));
194   for (number_t i=0; i < nx*ny*nz; ++i)
195   {
196     if (spaceDim == 1) { nodes[i]=Point(coords[i][0]); }
197     else if (spaceDim == 2) { nodes[i]=Point(coords[i][0], coords[i][1]); }
198     else { nodes[i]=Point(coords[i][0], coords[i][1], coords[i][2]); }
199   }
200 
201   if (nx == 1 && ny == 1) // segments along the z axis
202   {
203     RefElement* re_p = findRefElement(_segment, interp_p);
204     elements_.resize(nz-1);
205     for (number_t k=0; k < nz-1; ++k)
206     {
207       elements_[k]=new GeomElement(this, re_p, 1, k+1);
208       MeshElement* melt=elements_[k]->meshElement();
209       melt->nodeNumbers.resize(2);
210       melt->vertexNumbers.resize(2);
211       melt->vertexNumbers[0] = melt->nodeNumbers[0] = k+1;
212       melt->vertexNumbers[1] = melt->nodeNumbers[1] = k+2;
213       melt->setNodes(nodes);
214     }
215   }
216   else if (nx == 1 && nz == 1) // segments along the y axis
217   {
218     RefElement* re_p = findRefElement(_segment, interp_p);
219     elements_.resize(ny-1);
220     for (number_t j=0; j < ny-1; ++j)
221     {
222       elements_[j]=new GeomElement(this, re_p, 1, j+1);
223       MeshElement* melt=elements_[j]->meshElement();
224       melt->nodeNumbers.resize(2);
225       melt->vertexNumbers.resize(2);
226       melt->vertexNumbers[0] = melt->nodeNumbers[0] = j+1;
227       melt->vertexNumbers[1] = melt->nodeNumbers[1] = j+2;
228       melt->setNodes(nodes);
229     }
230   }
231   else if (ny == 1 && nz == 1) // segments along the x axis
232   {
233     RefElement* re_p = findRefElement(_segment, interp_p);
234     elements_.resize(nx-1);
235     for (number_t i=0; i < nx-1; ++i)
236     {
237       elements_[i]=new GeomElement(this, re_p, 1, i+1);
238       MeshElement* melt=elements_[i]->meshElement();
239       melt->nodeNumbers.resize(2);
240       melt->vertexNumbers.resize(2);
241       melt->vertexNumbers[0] = melt->nodeNumbers[0] = i+1;
242       melt->vertexNumbers[1] = melt->nodeNumbers[1] = i+2;
243       melt->setNodes(nodes);
244     }
245   }
246   else if (nx == 1) // quadrangles along the yz plane
247   {
248     RefElement* re_p = findRefElement(_quadrangle, interp_p);
249     elements_.resize((ny-1)*(nz-1));
250     for (number_t k=0; k < nz-1; ++k)
251     {
252       for (number_t j=0; j < ny-1; ++j)
253       {
254         number_t l=k*ny+j;
255         elements_[l]=new GeomElement(this, re_p, 2, l+1);
256         MeshElement* melt=elements_[l]->meshElement();
257         melt->nodeNumbers.resize(4);
258         melt->vertexNumbers.resize(4);
259         melt->vertexNumbers[0] = melt->nodeNumbers[0] = l+1;
260         melt->vertexNumbers[1] = melt->nodeNumbers[1] = l+2;
261         melt->vertexNumbers[2] = melt->nodeNumbers[2] = l+ny+1;
262         melt->vertexNumbers[3] = melt->nodeNumbers[3] = l+ny+2;
263         melt->setNodes(nodes);
264       }
265     }
266   }
267   else if (ny == 1) // quadrangles along the xz plane
268   {
269     RefElement* re_p = findRefElement(_quadrangle, interp_p);
270     elements_.resize((nx-1)*(nz-1));
271     for (number_t k=0; k < nz-1; ++k)
272     {
273       for (number_t i=0; i < nx-1; ++i)
274       {
275         number_t l=k*nx+i;
276         elements_[l]=new GeomElement(this, re_p, 2, l+1);
277         MeshElement* melt=elements_[l]->meshElement();
278         melt->nodeNumbers.resize(4);
279         melt->vertexNumbers.resize(4);
280         melt->vertexNumbers[0] = melt->nodeNumbers[0] = l+1;
281         melt->vertexNumbers[1] = melt->nodeNumbers[1] = l+2;
282         melt->vertexNumbers[2] = melt->nodeNumbers[2] = l+nx+1;
283         melt->vertexNumbers[3] = melt->nodeNumbers[3] = l+nx+2;
284         melt->setNodes(nodes);
285       }
286     }
287   }
288   else if (nz == 1) // quadrangles along the xy plane
289   {
290     RefElement* re_p = findRefElement(_quadrangle, interp_p);
291     elements_.resize((nx-1)*(ny-1));
292     for (number_t j=0; j < ny-1; ++j)
293     {
294       for (number_t i=0; i < nx-1; ++i)
295       {
296         number_t l=j*nx+i;
297         elements_[l]=new GeomElement(this, re_p, 2, l+1);
298         MeshElement* melt=elements_[l]->meshElement();
299         melt->nodeNumbers.resize(4);
300         melt->vertexNumbers.resize(4);
301         melt->vertexNumbers[0] = melt->nodeNumbers[0] = l+1;
302         melt->vertexNumbers[1] = melt->nodeNumbers[1] = l+2;
303         melt->vertexNumbers[2] = melt->nodeNumbers[2] = l+nx+1;
304         melt->vertexNumbers[3] = melt->nodeNumbers[3] = l+nx+2;
305         melt->setNodes(nodes);
306       }
307     }
308   }
309   else // hexahedra
310   {
311     RefElement* re_p = findRefElement(_hexahedron, interp_p);
312     elements_.resize((nx-1)*(ny-1)*(nz-1));
313     for (number_t k=0; k < nz-1; ++k)
314     {
315       for (number_t j=0; j < ny-1; ++j)
316       {
317         for (number_t i=0; i < nx-1; ++i)
318         {
319           number_t l=k*nx+ny+j*ny+i;
320           elements_[l]=new GeomElement(this, re_p, 3, l+1);
321           MeshElement* melt=elements_[l]->meshElement();
322           melt->nodeNumbers.resize(8);
323           melt->vertexNumbers.resize(8);
324           melt->vertexNumbers[0] = melt->nodeNumbers[0] = l+1;
325           melt->vertexNumbers[1] = melt->nodeNumbers[1] = l+2;
326           melt->vertexNumbers[2] = melt->nodeNumbers[2] = l+nx+1;
327           melt->vertexNumbers[3] = melt->nodeNumbers[3] = l+nx+2;
328           melt->vertexNumbers[4] = melt->nodeNumbers[4] = l+nx*ny+1;
329           melt->vertexNumbers[5] = melt->nodeNumbers[5] = l+nx*ny+2;
330           melt->vertexNumbers[6] = melt->nodeNumbers[6] = l+nx*ny+nx+1;
331           melt->vertexNumbers[7] = melt->nodeNumbers[7] = l+nx*ny+nx+2;
332           melt->setNodes(nodes);
333         }
334       }
335     }
336   }
337   return spaceDim;
338 }
339 
340 //! loads mesh from input stream according to STRUCTURED_GRID dataset in vtk mesh format
loadStructuredGridVtk(std::istream & data,number_t nodesDim)341 dimen_t Mesh::loadStructuredGridVtk(std::istream& data, number_t nodesDim)
342 {
343   string_t strVal;
344 
345   data >> strVal;
346   if (strVal != "DIMENSIONS") { error("vtk_bad_block", words("First"), "STRUCTURED_GRID", "DIMENSIONS"); }
347   number_t nx, ny, nz;
348   data >> nx >> ny >> nz;
349 
350   data >> strVal;
351   if (strVal != "POINTS") { error("vtk_bad_block", words("Second"), "STRUCTURED_GRID", "POINTS"); }
352   number_t nbPoints;
353   data >> nbPoints;
354   if (nbPoints != nx*ny*nz) { error("bad_size", "", nx*ny*nz, nbPoints); }
355   data >> strVal; // data type of coordinates (ignored)
356 
357   std::vector<std::vector<real_t> > coords(nbPoints,std::vector<real_t>(3));
358   vertices_.resize(nbPoints);
359   for (number_t i=0; i < nbPoints; ++i)
360   {
361     data >> coords[i][0] >> coords[i][1] >> coords[i][2];
362     vertices_[i]=i+1;
363   }
364 
365   // if nx=1 and ny=1 -> segments along the z axis
366   // if nx=1 and nz=1 -> segments along the y axis
367   // if ny=1 and nz=1 -> segments along the x axis
368   // if nx=1 -> quadrangles along the yz plane
369   // if ny=1 -> quadrangles along the xz plane
370   // if nz=1 -> quadrangles along the xy plane
371   // else hexahedra
372 
373   Interpolation* interp_p;
374   if (ny*nz==1 || nx*nz==1 || nx*ny==1)
375   {
376     // segments are P1
377     interp_p=findInterpolation(_Lagrange, _standard, _P1, _H1);
378   }
379   else
380   {
381     // quadrangles and hexahedra are Q1
382     interp_p=findInterpolation(_Lagrange, _standard, _Q1, _H1);
383   }
384 
385   dimen_t spaceDim=3;
386   if (nx == 1 && ny == 1) { spaceDim=1; } // segments along the z axis
387   else if (nx == 1 && nz == 1) { spaceDim=1; } // segments along the y axis
388   else if (ny == 1 && nz == 1) { spaceDim=1; } // segments along the x axis
389   else if (nx == 1) { spaceDim=2; } // quadrangles along the yz plane
390   else if (ny == 1) { spaceDim=2; } // quadrangles along the xz plane
391   else if (nz == 1) { spaceDim=2; } // quadrangles along the xy plane
392   if (nodesDim != 0 && nodesDim > spaceDim) { spaceDim=nodesDim; }
393 
394   nodes.resize(nbPoints,Point(std::vector<real_t>(spaceDim, 0.)));
395   for (number_t i=0; i < nbPoints; ++i)
396   {
397     if (spaceDim == 1) { nodes[i]=Point(coords[i][0]); }
398     else if (spaceDim == 2) { nodes[i]=Point(coords[i][0], coords[i][1]); }
399     else { nodes[i]=Point(coords[i][0], coords[i][1], coords[i][2]); }
400   }
401 
402   if (nx == 1 && ny == 1) // segments along the z axis
403   {
404     RefElement* re_p = findRefElement(_segment, interp_p);
405     elements_.resize(nz-1);
406     for (number_t k=0; k < nz-1; ++k)
407     {
408       elements_[k]=new GeomElement(this, re_p, 1, k+1);
409       MeshElement* melt=elements_[k]->meshElement();
410       melt->nodeNumbers.resize(2);
411       melt->vertexNumbers.resize(2);
412       melt->vertexNumbers[0] = melt->nodeNumbers[0] = k+1;
413       melt->vertexNumbers[1] = melt->nodeNumbers[1] = k+2;
414       melt->setNodes(nodes);
415     }
416   }
417   else if (nx == 1 && nz == 1) // segments along the y axis
418   {
419     RefElement* re_p = findRefElement(_segment, interp_p);
420     elements_.resize(ny-1);
421     for (number_t j=0; j < ny-1; ++j)
422     {
423       elements_[j]=new GeomElement(this, re_p, 1, j+1);
424       MeshElement* melt=elements_[j]->meshElement();
425       melt->nodeNumbers.resize(2);
426       melt->vertexNumbers.resize(2);
427       melt->vertexNumbers[0] = melt->nodeNumbers[0] = j+1;
428       melt->vertexNumbers[1] = melt->nodeNumbers[1] = j+2;
429       melt->setNodes(nodes);
430     }
431   }
432   else if (ny == 1 && nz == 1) // segments along the x axis
433   {
434     RefElement* re_p = findRefElement(_segment, interp_p);
435     elements_.resize(nx-1);
436     for (number_t i=0; i < nx-1; ++i)
437     {
438       elements_[i]=new GeomElement(this, re_p, 1, i+1);
439       MeshElement* melt=elements_[i]->meshElement();
440       melt->nodeNumbers.resize(2);
441       melt->vertexNumbers.resize(2);
442       melt->vertexNumbers[0] = melt->nodeNumbers[0] = i+1;
443       melt->vertexNumbers[1] = melt->nodeNumbers[1] = i+2;
444       melt->setNodes(nodes);
445     }
446   }
447   else if (nx == 1) // quadrangles along the yz plane
448   {
449     RefElement* re_p = findRefElement(_quadrangle, interp_p);
450     elements_.resize((ny-1)*(nz-1));
451     for (number_t k=0; k < nz-1; ++k)
452     {
453       for (number_t j=0; j < ny-1; ++j)
454       {
455         number_t l=k*ny+j;
456         elements_[l]=new GeomElement(this, re_p, 2, l+1);
457         MeshElement* melt=elements_[l]->meshElement();
458         melt->nodeNumbers.resize(4);
459         melt->vertexNumbers.resize(4);
460         melt->vertexNumbers[0] = melt->nodeNumbers[0] = l+1;
461         melt->vertexNumbers[1] = melt->nodeNumbers[1] = l+2;
462         melt->vertexNumbers[2] = melt->nodeNumbers[2] = l+ny+1;
463         melt->vertexNumbers[3] = melt->nodeNumbers[3] = l+ny+2;
464         melt->setNodes(nodes);
465       }
466     }
467   }
468   else if (ny == 1) // quadrangles along the xz plane
469   {
470     RefElement* re_p = findRefElement(_quadrangle, interp_p);
471     elements_.resize((nx-1)*(nz-1));
472     for (number_t k=0; k < nz-1; ++k)
473     {
474       for (number_t i=0; i < nx-1; ++i)
475       {
476         number_t l=k*nx+i;
477         elements_[l]=new GeomElement(this, re_p, 2, l+1);
478         MeshElement* melt=elements_[l]->meshElement();
479         melt->nodeNumbers.resize(4);
480         melt->vertexNumbers.resize(4);
481         melt->vertexNumbers[0] = melt->nodeNumbers[0] = l+1;
482         melt->vertexNumbers[1] = melt->nodeNumbers[1] = l+2;
483         melt->vertexNumbers[2] = melt->nodeNumbers[2] = l+nx+1;
484         melt->vertexNumbers[3] = melt->nodeNumbers[3] = l+nx+2;
485         melt->setNodes(nodes);
486       }
487     }
488   }
489   else if (nz == 1) // quadrangles along the xy plane
490   {
491     RefElement* re_p = findRefElement(_quadrangle, interp_p);
492     elements_.resize((nx-1)*(ny-1));
493     for (number_t j=0; j < ny-1; ++j)
494     {
495       for (number_t i=0; i < nx-1; ++i)
496       {
497         number_t l=j*nx+i;
498         elements_[l]=new GeomElement(this, re_p, 2, l+1);
499         MeshElement* melt=elements_[l]->meshElement();
500         melt->nodeNumbers.resize(4);
501         melt->vertexNumbers.resize(4);
502         melt->vertexNumbers[0] = melt->nodeNumbers[0] = l+1;
503         melt->vertexNumbers[1] = melt->nodeNumbers[1] = l+2;
504         melt->vertexNumbers[2] = melt->nodeNumbers[2] = l+nx+1;
505         melt->vertexNumbers[3] = melt->nodeNumbers[3] = l+nx+2;
506         melt->setNodes(nodes);
507       }
508     }
509   }
510   else // hexahedra
511   {
512     RefElement* re_p = findRefElement(_hexahedron, interp_p);
513     elements_.resize((nx-1)*(ny-1)*(nz-1));
514     for (number_t k=0; k < nz-1; ++k)
515     {
516       for (number_t j=0; j < ny-1; ++j)
517       {
518         for (number_t i=0; i < nx-1; ++i)
519         {
520           number_t l=k*nx+ny+j*ny+i;
521           elements_[l]=new GeomElement(this, re_p, 3, l+1);
522           MeshElement* melt=elements_[l]->meshElement();
523           melt->nodeNumbers.resize(8);
524           melt->vertexNumbers.resize(8);
525           melt->vertexNumbers[0] = melt->nodeNumbers[0] = l+1;
526           melt->vertexNumbers[1] = melt->nodeNumbers[1] = l+2;
527           melt->vertexNumbers[2] = melt->nodeNumbers[2] = l+nx+1;
528           melt->vertexNumbers[3] = melt->nodeNumbers[3] = l+nx+2;
529           melt->vertexNumbers[4] = melt->nodeNumbers[4] = l+nx*ny+1;
530           melt->vertexNumbers[5] = melt->nodeNumbers[5] = l+nx*ny+2;
531           melt->vertexNumbers[6] = melt->nodeNumbers[6] = l+nx*ny+nx+1;
532           melt->vertexNumbers[7] = melt->nodeNumbers[7] = l+nx*ny+nx+2;
533           melt->setNodes(nodes);
534         }
535       }
536     }
537   }
538   return spaceDim;
539 }
540 
541 //! loads mesh from input stream according to UNSTRUCTURED_GRID dataset in vtk mesh format
loadUnstructuredGridVtk(std::istream & data,number_t nodesDim)542 dimen_t Mesh::loadUnstructuredGridVtk(std::istream& data, number_t nodesDim)
543 {
544   string_t strVal;
545   number_t numVal;
546 
547   data >> strVal;
548   if (strVal != "POINTS") { error("vtk_bad_block", words("First"), "UNSTRUCTURED_GRID", "POINTS"); }
549   number_t nbPoints;
550   data >> nbPoints;
551   data >> strVal; // data type of points (ignored)
552 
553   // Definition of nodes
554   std::vector<std::vector<real_t> > coords(nbPoints,std::vector<real_t>(3));
555   vertices_.resize(nbPoints);
556   for (number_t i=0; i < nbPoints; ++i)
557   {
558     data >> coords[i][0] >> coords[i][1] >> coords[i][2];
559     vertices_[i]=i+1;
560   }
561 
562   data >> strVal;
563   if (strVal != "CELLS") { error("vtk_bad_block", words("Second"), "UNSTRUCTURED_GRID", "CELLS"); }
564   number_t nbElems;
565   data >> nbElems; // number of cells = number of elements
566   data >> numVal; // number of data to be read (useless)
567   std::vector<std::vector<number_t> > cells(nbElems,std::vector<number_t>(8));
568   std::vector<ShapeType> celltypes(nbElems);
569 
570   for (number_t c=0; c < nbElems; ++c)
571   {
572     number_t cellSize;
573     data >> cellSize;
574     cells[c].resize(cellSize);
575     for (number_t i=0; i < cellSize; ++i)
576     {
577       data >> numVal;
578       cells[c][i]=numVal+1;
579     }
580   }
581 
582   data >> strVal;
583   if (strVal != "CELL_TYPES") { error("vtk_bad_block", words("Third"), "UNSTRUCTURED_GRID", "CELL_TYPES"); }
584   data >> numVal; // number of cell types = number of cells
585   if (numVal != nbElems) { error("bad_size", "", nbElems, numVal); }
586   data >> numVal; // number of data to be read (useless)
587 
588   std::map<number_t, ShapeType> vtkShape;
589   std::map<number_t, number_t> vtkDim;
590   vtkShape[1]=_point;         vtkDim[1]=0;
591   vtkShape[3]=_segment;       vtkDim[3]=1;
592   vtkShape[5]=_triangle;      vtkDim[5]=2;
593   vtkShape[9]=_quadrangle;    vtkDim[9]=2;
594   vtkShape[10]=_tetrahedron;  vtkDim[10]=3;
595   vtkShape[12]=_hexahedron;   vtkDim[12]=3;
596   vtkShape[13]=_prism;        vtkDim[13]=3;
597   vtkShape[14]=_pyramid;      vtkDim[14]=3;
598 
599   dimen_t spaceDim=1;
600   for (number_t c=0; c < nbElems; ++c)
601   {
602     data >> numVal;
603     celltypes[c]=vtkShape[numVal];
604     if (vtkDim[numVal] > spaceDim) { spaceDim = vtkDim[numVal]; }
605   }
606   if (nodesDim != 0 && nodesDim > spaceDim) { spaceDim=nodesDim; }
607 
608   if (spaceDim == 2)
609   {
610     // we have to check if the geometry is 2D or 3D
611     if (coords.size() > 4)
612     {
613       // we take 4 vertices randomly and check their coplanarity
614       number_t v1=0, v2=(nbPoints-1)/3, v3=(2*(nbPoints-1))/3, v4=nbPoints-1;
615       if ( std::abs(dot(coords[v2]-coords[v1],crossProduct(coords[v3]-coords[v1],coords[v4]-coords[v1]))) > theEpsilon ) { spaceDim=3; }
616     }
617   }
618 
619   nodes.resize(nbPoints,Point(std::vector<real_t>(spaceDim, 0.)));
620   for (number_t i=0; i < nbPoints; ++i)
621   {
622     if (spaceDim == 1) { nodes[i]=Point(coords[i][0]); }
623     else if (spaceDim == 2) { nodes[i]=Point(coords[i][0], coords[i][1]); }
624     else { nodes[i]=Point(coords[i][0], coords[i][1], coords[i][2]); }
625   }
626 
627   elements_.resize(nbElems);
628   for (number_t c =0; c < nbElems; ++c)
629   {
630     Interpolation* interp_p;
631     if (celltypes[c] == _quadrangle || celltypes[c] == _hexahedron)
632     {
633       interp_p=findInterpolation(_Lagrange, _standard, _Q1, _H1);
634     }
635     else
636     {
637       interp_p=findInterpolation(_Lagrange, _standard, _P1, _H1);
638     }
639     RefElement* re_p=findRefElement(celltypes[c], interp_p);
640 
641     // create GeomElement
642     elements_[c]=new GeomElement(this, re_p, spaceDim, c+1);
643     MeshElement * melt=elements_[c]->meshElement();
644     number_t nbVertices = cells[c].size();
645     melt->vertexNumbers.resize(nbVertices);
646     melt->nodeNumbers.resize(nbVertices);
647     for (number_t i=0; i < nbVertices; ++i)
648     {
649       melt->vertexNumbers[i] = melt->nodeNumbers[i] = cells[c][i];
650     }
651     melt->setNodes(nodes);
652   }
653 
654   return spaceDim;
655 }
656 
657 
658 //! loads mesh from input file filename according to vtk mesh format
loadVtk(const string_t & filename,number_t nodesDim)659 void Mesh::loadVtk(const string_t& filename, number_t nodesDim)
660 {
661   trace_p->push("Mesh::loadVtk");
662   // error("not_yet_implemented", "Mesh::loadVtk(String filename)");
663 
664   std::ifstream data(filename.c_str());
665   if ( !data ) { error("file_failopen","Mesh::loadVtk",filename); }
666 
667   comment_ = string_t("VTK mesh format read from file ") + basenameWithExtension(filename);
668   if ( theVerboseLevel > 1) { info("loadFile_info", "VTK", filename); }
669 
670   string_t strVal;
671 
672   data >> strVal;
673   bool firstLineError = false;
674 
675   if (strVal != "#") { firstLineError=true; }
676   if (!firstLineError) { data >> strVal; }
677   if (!firstLineError && strVal != "vtk") { firstLineError=true; }
678   if (!firstLineError) { data >> strVal; }
679   if (!firstLineError && strVal != "DataFile") { firstLineError=true; }
680   if (!firstLineError) { data >> strVal; }
681   if (!firstLineError && strVal != "Version") { firstLineError=true; }
682   if (firstLineError) { error("vtk_format_block_first_line"); }
683   real_t version;
684   data >> version;
685 
686   // getting until end of line (normally it is empty)
687   getline(data,strVal);
688   // next line is a comment line
689   getline(data,strVal);
690 
691   data >> strVal;
692   if (strVal != "ASCII") { error("non_ascii"); }
693   data >> strVal;
694   if (strVal != "DATASET") { error("vtk_bad_format"); }
695   string_t dataset;
696   data >> dataset;
697   dimen_t spaceDim=0;
698   if (dataset == "POLYDATA") { spaceDim=loadPolyDataVtk(data, nodesDim); }
699   else if (dataset == "RECTILINEAR_GRID") { spaceDim=loadRectilinearVtk(data, nodesDim); }
700   else if (dataset == "STRUCTURED_GRID") { spaceDim=loadStructuredGridVtk(data, nodesDim); }
701   else if (dataset == "UNSTRUCTURED_GRID") { spaceDim=loadUnstructuredGridVtk(data, nodesDim); }
702   else { error("vtk_bad_format"); }
703 
704   // management of domain
705   MeshDomain* meshdom_p = (new GeomDomain(*this, "Omega", spaceDim))->meshDomain();
706   meshdom_p->geomElements.assign(elements_.begin(), elements_.end());
707   domains_.push_back(meshdom_p);
708   order_=1;
709   lastIndex_ = elements_.size();
710 
711   //compute measures and orientation of mesh elements
712   buildGeomData();
713   setShapeTypes();
714 
715   //sides and sidesOfSides lists are not built (see buildSides and buildSideOfSides)
716 
717   firstOrderMesh_p = this;
718 
719   trace_p->pop();
720 }
721 
722 //! loads mesh from input file filename according to vtp mesh format
loadVtp(const string_t & filename,number_t nodesDim)723 void Mesh::loadVtp(const string_t& filename, number_t nodesDim)
724 {
725   trace_p->push("Mesh::loadVtp");
726   error("not_yet_implemented", "Mesh::loadVtp(String filename)");
727 
728   comment_ = string_t("VTP mesh format read from file ") + basenameWithExtension(filename);
729   if ( theVerboseLevel > 1) { info("loadFile_info", "VTP", filename); }
730 
731   trace_p->pop();
732 }
733 
734 //! loads mesh from input file filename according to vtr mesh format
loadVtr(const string_t & filename,number_t nodesDim)735 void Mesh::loadVtr(const string_t& filename, number_t nodesDim)
736 {
737   trace_p->push("Mesh::loadVtr");
738   error("not_yet_implemented", "Mesh::loadVtr(String filename)");
739 
740   comment_ = string_t("VTR mesh format read from file ") + basenameWithExtension(filename);
741   if ( theVerboseLevel > 1) { info("loadFile_info", "VTR", filename); }
742 
743   trace_p->pop();
744 }
745 
746 //! loads mesh from input file filename according to vts mesh format
loadVts(const string_t & filename,number_t nodesDim)747 void Mesh::loadVts(const string_t& filename, number_t nodesDim)
748 {
749   trace_p->push("Mesh::loadVts");
750   error("not_yet_implemented", "Mesh::loadVts(String filename)");
751 
752   comment_ = string_t("VTS mesh format read from file ") + basenameWithExtension(filename);
753   if ( theVerboseLevel > 1) { info("loadFile_info", "VTS", filename); }
754 
755   trace_p->pop();
756 }
757 
758 //! loads mesh from input file filename according to vtu mesh format
loadVtu(const string_t & filename,number_t nodesDim)759 void Mesh::loadVtu(const string_t& filename, number_t nodesDim)
760 {
761   trace_p->push("Mesh::loadVtu");
762   error("not_yet_implemented", "Mesh::loadVtu(String filename)");
763 
764   comment_ = string_t("VTU mesh format read from file ") + basenameWithExtension(filename);
765   if ( theVerboseLevel > 1) { info("loadFile_info", "VTU", filename); }
766 
767   trace_p->pop();
768 }
769 
770 } // end of namespace xlifepp
771