1 // -*- C++ -*-
2 /**
3 * @brief Mesh is defined as a collection of n-D vertices and simplices
4 *
5 * Copyright 2005-2021 Airbus-EDF-IMACS-ONERA-Phimeca
6 *
7 * This library is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU Lesser General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * This library is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU Lesser General Public License for more details.
16 *
17 * You should have received a copy of the GNU Lesser General Public License
18 * along with this library. If not, see <http://www.gnu.org/licenses/>.
19 *
20 */
21 #include <fstream>
22 #include <algorithm>
23
24 #include "openturns/Mesh.hxx"
25 #include "openturns/PersistentObjectFactory.hxx"
26 #include "openturns/Log.hxx"
27 #include "openturns/Os.hxx"
28 #include "openturns/Exception.hxx"
29 #include "openturns/Graph.hxx"
30 #include "openturns/Cloud.hxx"
31 #include "openturns/Curve.hxx"
32 #include "openturns/Polygon.hxx"
33 #include "openturns/PolygonArray.hxx"
34 #include "openturns/Collection.hxx"
35 #include "openturns/SpecFunc.hxx"
36 #include "openturns/PlatformInfo.hxx"
37
38 BEGIN_NAMESPACE_OPENTURNS
39
40 CLASSNAMEINIT(Mesh)
41
42 static const Factory<Mesh> Factory_Mesh;
43
44 /* Default constructor */
Mesh(const UnsignedInteger dimension)45 Mesh::Mesh(const UnsignedInteger dimension)
46 : PersistentObject()
47 , dimension_(dimension)
48 , hasBeenChecked_(false)
49 , vertices_(1, dimension) // At least one point
50 , simplices_()
51 {
52 // Nothing to do
53 if (vertices_.getDescription().isBlank()) vertices_.setDescription(Description::BuildDefault(dimension, "t"));
54 }
55
56 /* Parameters constructor, simplified interface for 1D case */
Mesh(const Sample & vertices)57 Mesh::Mesh(const Sample & vertices)
58 : PersistentObject()
59 , dimension_(vertices.getDimension())
60 , hasBeenChecked_(false)
61 , vertices_(0, vertices.getDimension())
62 , simplices_()
63 {
64 // Use the vertices accessor to initialize the kd-tree
65 setVertices(vertices);
66 }
67
68 /* Parameters constructor, simplified interface for 1D case */
Mesh(const Sample & vertices,const IndicesCollection & simplices,const Bool checkMeshValidity)69 Mesh::Mesh(const Sample & vertices,
70 const IndicesCollection & simplices,
71 const Bool checkMeshValidity)
72 : PersistentObject()
73 , dimension_(vertices.getDimension())
74 , hasBeenChecked_(false)
75 , vertices_(0, vertices.getDimension())
76 , simplices_(simplices)
77 {
78 // Use the vertices accessor to initialize the kd-tree
79 setVertices(vertices);
80 if (checkMeshValidity)
81 checkValidity();
82 }
83
84 /* Clone method */
clone() const85 Mesh * Mesh::clone() const
86 {
87 return new Mesh(*this);
88 }
89
90 /* Dimension accessor */
getDimension() const91 UnsignedInteger Mesh::getDimension() const
92 {
93 return dimension_;
94 }
95
96 /* Description of the vertices accessor */
setDescription(const Description & description)97 void Mesh::setDescription(const Description & description)
98 {
99 vertices_.setDescription(description);
100 }
101
getDescription() const102 Description Mesh::getDescription() const
103 {
104 return vertices_.getDescription();
105 }
106
107 /* Vertices accessor */
getVertices() const108 Sample Mesh::getVertices() const
109 {
110 return vertices_;
111 }
112
setVertices(const Sample & vertices)113 void Mesh::setVertices(const Sample & vertices)
114 {
115 vertices_ = vertices;
116 if (vertices_.getDescription().isBlank()) vertices_.setDescription(Description::BuildDefault(vertices_.getDimension(), "t"));
117 hasBeenChecked_ = false;
118 }
119
120 /* Vertex accessor */
getVertex(const UnsignedInteger index) const121 Point Mesh::getVertex(const UnsignedInteger index) const
122 {
123 if (!(index < getVerticesNumber())) throw InvalidArgumentException(HERE) << "Error: the vertex index=" << index << " must be less than the number of vertices=" << getVerticesNumber();
124 return vertices_[index];
125 }
126
setVertex(const UnsignedInteger index,const Point & vertex)127 void Mesh::setVertex(const UnsignedInteger index,
128 const Point & vertex)
129 {
130 vertices_[index] = vertex;
131 hasBeenChecked_ = false;
132 }
133
134 /* Simplices accessor */
getSimplices() const135 IndicesCollection Mesh::getSimplices() const
136 {
137 return simplices_;
138 }
139
setSimplices(const IndicesCollection & simplices)140 void Mesh::setSimplices(const IndicesCollection & simplices)
141 {
142 if (!(simplices == simplices_))
143 {
144 simplices_ = simplices;
145 }
146 hasBeenChecked_ = false;
147 }
148
149 /* Simplex accessor */
getSimplex(const UnsignedInteger index) const150 Indices Mesh::getSimplex(const UnsignedInteger index) const
151 {
152 if (!(index < getSimplicesNumber())) throw InvalidArgumentException(HERE) << "Error: the simplex index=" << index << " must be less than the number of simplices=" << getSimplicesNumber();
153 return Indices(simplices_.cbegin_at(index), simplices_.cend_at(index));
154 }
155
156 /* Check the mesh validity */
checkValidity() const157 void Mesh::checkValidity() const
158 {
159
160 if (hasBeenChecked_) return;
161 // Check the vertices: no duplicate, no unused vertex
162 // Check the simplices: no simplex with duplicate vertices, no simplex with unknown vertex, no simplex with a number of vertices different from dimension+1
163 for (UnsignedInteger i = 0; i < getSimplicesNumber(); ++ i)
164 {
165 Indices simplex(getSimplex(i));
166 if (simplex.getSize() != getDimension() + 1)
167 throw InvalidArgumentException(HERE) << "Error: mesh has dimension " << getDimension() << " but simplex #" << i << " has size" << simplex.getSize();
168
169 if (!simplex.check(getVerticesNumber()))
170 throw InvalidArgumentException(HERE) << "Error: mesh has " << getVerticesNumber() << " vertices but simplex #" << i << " refers to an unknown vertex";
171 }
172 // Check that no ball can be included into the intersection of two simplices
173 // One it has been checked everything is ok
174 hasBeenChecked_ = true;
175 }
176
isValid() const177 Bool Mesh::isValid() const
178 {
179 Bool result = true;
180 try
181 {
182 checkValidity();
183 }
184 catch (Exception &)
185 {
186 result = false;
187 }
188 return result;
189 }
190
191 /* Build the affine matrix associated to the simplex at the given index*/
buildSimplexMatrix(const UnsignedInteger index,SquareMatrix & matrix) const192 void Mesh::buildSimplexMatrix(const UnsignedInteger index,
193 SquareMatrix & matrix) const
194 {
195 if (!(index < getSimplicesNumber())) throw InvalidArgumentException(HERE) << "Error: the simplex index=" << index << " must be less than the number of simplices=" << getSimplicesNumber();
196 if (matrix.getDimension() != dimension_ + 1)
197 matrix = SquareMatrix(dimension_ + 1);
198 // Loop over the vertices of the simplex
199 for (UnsignedInteger j = 0; j <= dimension_; ++j)
200 {
201 const UnsignedInteger vertexJ = simplices_(index, j);
202 for (UnsignedInteger i = 0; i < dimension_; ++i) matrix(i, j) = vertices_(vertexJ, i);
203 matrix(dimension_, j) = 1.0;
204 }
205 }
206
207 /* Check if the given point is in the given simplex and returns its barycentric coordinates */
checkPointInSimplexWithCoordinates(const Point & point,const UnsignedInteger index,Point & coordinates) const208 Bool Mesh::checkPointInSimplexWithCoordinates(const Point & point,
209 const UnsignedInteger index,
210 Point & coordinates) const
211 {
212 if (!(index < getSimplicesNumber())) return false;
213 const Scalar epsilon = ResourceMap::GetAsScalar("Mesh-VertexEpsilon");
214 if (dimension_ == 1)
215 {
216 const Scalar x = point[0];
217 const Scalar x0 = vertices_(simplices_(index, 0), 0);
218 const Scalar x1 = vertices_(simplices_(index, 1), 0);
219 if ((x - x0) * (x - x1) > epsilon)
220 return false;
221 coordinates = Point(2);
222 if (x0 == x1)
223 {
224 // x, x0 and x1 are amost at the same position, any value would work.
225 coordinates[0] = 1.0;
226 coordinates[1] = 0.0;
227 }
228 else
229 {
230 const Scalar alpha = (x1 - x) / (x1 - x0);
231 coordinates[0] = alpha;
232 coordinates[1] = 1.0 - alpha;
233 }
234 return true;
235 }
236 else if (dimension_ == 2)
237 {
238 const Scalar x0 = vertices_(simplices_(index, 0), 0);
239 const Scalar y0 = vertices_(simplices_(index, 0), 1);
240 const Scalar x01 = vertices_(simplices_(index, 1), 0) - x0;
241 const Scalar y01 = vertices_(simplices_(index, 1), 1) - y0;
242 const Scalar x02 = vertices_(simplices_(index, 2), 0) - x0;
243 const Scalar y02 = vertices_(simplices_(index, 2), 1) - y0;
244 const Scalar det = (x02 * y01 - y02 * x01);
245 if (det == 0.0)
246 {
247 return false;
248 }
249 const Scalar x = point[0] - x0;
250 const Scalar y = point[1] - y0;
251 coordinates = Point(3);
252 coordinates[1] = (x02 * y - y02 * x) / det;
253 coordinates[2] = (x * y01 - y * x01) / det;
254 coordinates[0] = 0.5 + (0.5 - coordinates[1] - coordinates[2]);
255 return coordinates[0] >= -epsilon && coordinates[0] <= 1.0 + epsilon && coordinates[1] >= -epsilon && coordinates[1] <= 1.0 + epsilon && coordinates[2] >= -epsilon && coordinates[2] <= 1.0 + epsilon;
256 }
257 SquareMatrix matrix(dimension_ + 1);
258 buildSimplexMatrix(index, matrix);
259 Point v(point);
260 v.add(1.0);
261 coordinates = matrix.solveLinearSystem(v, false);
262 for (UnsignedInteger i = 0; i <= dimension_; ++i) if ((coordinates[i] < -epsilon) || (coordinates[i] > 1.0 + epsilon)) return false;
263 return true;
264 }
265
266 /* Get the number of vertices */
getVerticesNumber() const267 UnsignedInteger Mesh::getVerticesNumber() const
268 {
269 return vertices_.getSize();
270 }
271
272 /* Get the number of simplices */
getSimplicesNumber() const273 UnsignedInteger Mesh::getSimplicesNumber() const
274 {
275 return simplices_.getSize();
276 }
277
computeSimplicesVolume() const278 Point Mesh::computeSimplicesVolume() const
279 {
280 const UnsignedInteger nrSimplices = getSimplicesNumber();
281 Point result(nrSimplices);
282 if (nrSimplices == 0) return result;
283 if (getDimension() == 1)
284 {
285 for (UnsignedInteger index = 0; index < nrSimplices; ++index)
286 {
287 const Scalar x0 = vertices_(simplices_(index, 0), 0);
288 const Scalar x1 = vertices_(simplices_(index, 1), 0);
289 result[index] = std::abs(x1 - x0);
290 }
291 }
292 else if (getDimension() == 2)
293 {
294 for (UnsignedInteger index = 0; index < nrSimplices; ++index)
295 {
296 const Scalar x0 = vertices_(simplices_(index, 0), 0);
297 const Scalar y0 = vertices_(simplices_(index, 0), 1);
298 const Scalar x01 = vertices_(simplices_(index, 1), 0) - x0;
299 const Scalar y01 = vertices_(simplices_(index, 1), 1) - y0;
300 const Scalar x02 = vertices_(simplices_(index, 2), 0) - x0;
301 const Scalar y02 = vertices_(simplices_(index, 2), 1) - y0;
302 result[index] = 0.5 * std::abs(x02 * y01 - x01 * y02);
303 }
304 }
305 else if (getDimension() == 3)
306 {
307 for (UnsignedInteger index = 0; index < nrSimplices; ++index)
308 {
309 const Scalar x0 = vertices_(simplices_(index, 0), 0);
310 const Scalar y0 = vertices_(simplices_(index, 0), 1);
311 const Scalar z0 = vertices_(simplices_(index, 0), 2);
312 const Scalar x01 = vertices_(simplices_(index, 1), 0) - x0;
313 const Scalar y01 = vertices_(simplices_(index, 1), 1) - y0;
314 const Scalar z01 = vertices_(simplices_(index, 1), 2) - z0;
315 const Scalar x02 = vertices_(simplices_(index, 2), 0) - x0;
316 const Scalar y02 = vertices_(simplices_(index, 2), 1) - y0;
317 const Scalar z02 = vertices_(simplices_(index, 2), 2) - z0;
318 const Scalar x03 = vertices_(simplices_(index, 3), 0) - x0;
319 const Scalar y03 = vertices_(simplices_(index, 3), 1) - y0;
320 const Scalar z03 = vertices_(simplices_(index, 3), 2) - z0;
321 result[index] = std::abs(x01 * (y02 * z03 - z02 * y03) + y01 * (z02 * x03 - x02 * z03) + z01 * (x02 * y03 - y02 * x03)) / 6.0;
322 }
323 }
324 else
325 {
326 // General case
327 SquareMatrix matrix(getDimension() + 1);
328 Scalar sign = 0.0;
329 const Scalar logGamma = SpecFunc::LogGamma(getDimension() + 1);
330 for (UnsignedInteger index = 0; index < nrSimplices; ++index)
331 {
332 buildSimplexMatrix(index, matrix);
333 result[index] = exp(matrix.computeLogAbsoluteDeterminant(sign, false) - logGamma);
334 }
335 }
336 return result;
337 }
338
339 /* Compute P1 gram matrix */
computeP1Gram() const340 CovarianceMatrix Mesh::computeP1Gram() const
341 {
342 // If no simplex, the P1 gram matrix is null
343 if (simplices_.getSize() == 0) return CovarianceMatrix(0);
344 const UnsignedInteger simplexSize = getVertices().getDimension() + 1;
345 SquareMatrix elementaryGram(simplexSize, Point(simplexSize * simplexSize, 1.0 / (simplexSize * (simplexSize + 1.0))));
346 for (UnsignedInteger i = 0; i < simplexSize; ++i) elementaryGram(i, i) *= 2.0;
347 const UnsignedInteger verticesSize = vertices_.getSize();
348 const UnsignedInteger simplicesSize = simplices_.getSize();
349 const Point simplexVolume(computeSimplicesVolume());
350 SquareMatrix gram(verticesSize);
351 for (UnsignedInteger i = 0; i < simplicesSize; ++i)
352 {
353 const Indices simplex(getSimplex(i));
354 const Scalar delta = simplexVolume[i];
355 for (UnsignedInteger j = 0; j < simplexSize; ++j)
356 {
357 const UnsignedInteger newJ = simplex[j];
358 for (UnsignedInteger k = 0; k < simplexSize; ++k)
359 {
360 const UnsignedInteger newK = simplex[k];
361 gram(newJ, newK) += delta * elementaryGram(j, k);
362 } // Loop over second vertex
363 } // Loop over first vertex
364 } // Loop over simplices
365 return CovarianceMatrix(gram.getImplementation());
366 }
367
368 /* Get the numerical volume of the domain */
getVolume() const369 Scalar Mesh::getVolume() const
370 {
371 const Point volumes(computeSimplicesVolume());
372 return volumes.norm1();
373 }
374
375 /* Check if the domain is empty, i.e if its numerical volume is zero */
isEmpty() const376 Bool Mesh::isEmpty() const
377 {
378 return isNumericallyEmpty();
379 }
380
isNumericallyEmpty() const381 Bool Mesh::isNumericallyEmpty() const
382 {
383 return getVolume() <= ResourceMap::GetAsScalar("Domain-SmallVolume");
384 }
385
386 /* Tells if the mesh is regular */
isRegular() const387 Bool Mesh::isRegular() const
388 {
389 // For now, only 1D regular meshes are considered
390 if (getDimension() != 1) return false;
391 const UnsignedInteger size = getSimplicesNumber();
392 if (size <= 1) return true;
393 Bool regular = true;
394 const Scalar epsilon = ResourceMap::GetAsScalar("Mesh-VertexEpsilon");
395 const Scalar step = vertices_(simplices_(0, 1), 0) - vertices_(simplices_(0, 0), 0);
396 const Scalar absStep = std::abs(step);
397 for (UnsignedInteger i = 1; (i < size) && regular; ++ i)
398 {
399 const Scalar delta = vertices_(simplices_(i, 1), 0) - vertices_(simplices_(i, 0), 0);
400 regular = (std::abs(delta - step) <= absStep * epsilon);
401 }
402 return regular;
403 }
404
405 /* Lower bound of the bounding box */
getLowerBound() const406 Point Mesh::getLowerBound() const
407 {
408 return vertices_.getMin();
409 }
410
411 /* Upper bound of the bounding box */
getUpperBound() const412 Point Mesh::getUpperBound() const
413 {
414 return vertices_.getMax();
415 }
416
417 /* Orientation management */
fixOrientation()418 void Mesh::fixOrientation()
419 {
420 SquareMatrix matrix(dimension_ + 1);
421 const UnsignedInteger size = getSimplicesNumber();
422 for (UnsignedInteger i = 0; i < size; ++i)
423 fixOrientation(i, matrix);
424 }
425
fixOrientation(const UnsignedInteger & index,SquareMatrix & matrix)426 void Mesh::fixOrientation(const UnsignedInteger & index,
427 SquareMatrix & matrix)
428 {
429 if (getDimension() == 1)
430 {
431 IndicesCollection::iterator cit = simplices_.begin_at(index);
432 if (vertices_(*(cit + 1), 0) < vertices_(*cit, 0)) std::swap(*cit, *(cit + 1));
433 return;
434 }
435 if (getDimension() == 2)
436 {
437 IndicesCollection::iterator cit = simplices_.begin_at(index);
438 const Scalar v1x = vertices_(*cit, 0);
439 const Scalar v1y = vertices_(*cit, 1);
440 const Scalar v2x = vertices_(*(cit + 1), 0);
441 const Scalar v2y = vertices_(*(cit + 1), 1);
442 const Scalar v3x = vertices_(*(cit + 2), 0);
443 const Scalar v3y = vertices_(*(cit + 2), 1);
444 if ((v3x - v2x) * (v1y - v2y) < (v1x - v2x) * (v3y - v2y)) std::swap(*cit, *(cit + 1));
445 return;
446 }
447 if (getDimension() == 3)
448 {
449 IndicesCollection::iterator cit = simplices_.begin_at(index);
450 const Scalar v1x = vertices_(*cit, 0);
451 const Scalar v1y = vertices_(*cit, 1);
452 const Scalar v1z = vertices_(*cit, 2);
453 const Scalar v2x = vertices_(*(cit + 1), 0);
454 const Scalar v2y = vertices_(*(cit + 1), 1);
455 const Scalar v2z = vertices_(*(cit + 1), 2);
456 const Scalar v3x = vertices_(*(cit + 2), 0);
457 const Scalar v3y = vertices_(*(cit + 2), 1);
458 const Scalar v3z = vertices_(*(cit + 2), 2);
459 const Scalar v4x = vertices_(*(cit + 3), 0);
460 const Scalar v4y = vertices_(*(cit + 3), 1);
461 const Scalar v4z = vertices_(*(cit + 3), 2);
462 if ((v1x - v4x) * ((v2y - v4y) * (v3z - v4z) - (v3y - v4y) * (v2z - v4z)) + (v3x - v4x) * ((v1y - v4y) * (v2z - v4z) - (v2y - v4y) * (v1z - v4z)) < (v2x - v4x) * ((v1y - v4y) * (v3z - v4z) - (v3y - v4y) * (v1z - v4z)) ) std::swap(*cit, *(cit + 1));
463 return;
464 }
465 buildSimplexMatrix(index, matrix);
466 Scalar sign = 0.0;
467 (void) matrix.computeLogAbsoluteDeterminant(sign, false);
468 // In odd dimension the positive orientation is for a negative determinant of
469 // the simplex matrix
470 if ((sign > 0.0) == (dimension_ % 2 == 1)) return;
471 IndicesCollection::iterator cit = simplices_.begin_at(index);
472 std::swap(*cit, *(cit + 1));
473 }
474
475 /* Compute weights such that an integral of a function over the mesh
476 * is a weighted sum of its values at the vertices
477 */
computeWeights() const478 Point Mesh::computeWeights() const
479 {
480 // Compute the weights of the vertices by distributing the volume of each simplex among its vertices
481 checkValidity();
482 const UnsignedInteger numVertices = getVerticesNumber();
483 const UnsignedInteger numSimplices = getSimplicesNumber();
484 Point weights(numVertices, 0.0);
485 Point simplexVolume(computeSimplicesVolume());
486 for (UnsignedInteger simplex = 0; simplex < numSimplices; ++simplex)
487 {
488 const Scalar weight = simplexVolume[simplex];
489 for (IndicesCollection::const_iterator cit = simplices_.cbegin_at(simplex); cit != simplices_.cend_at(simplex); ++cit)
490 {
491 weights[*cit] += weight;
492 }
493 }
494 // Normalize the weights: each simplex has dim+1 vertices, so each vertex
495 // get 1/(dim+1) of the volume of the simplices it belongs to
496 weights /= (dimension_ + 1.0);
497 return weights;
498 }
499
500 /* Comparison operator */
operator ==(const Mesh & other) const501 Bool Mesh::operator == (const Mesh & other) const
502 {
503 if (this == &other) return true;
504 return (vertices_ == other.vertices_) && (simplices_ == other.simplices_);
505 }
506
507 /* String converter */
__repr__() const508 String Mesh::__repr__() const
509 {
510 return OSS(true) << "class=" << GetClassName()
511 << " name=" << getName()
512 << " dimension=" << getDimension()
513 << " vertices=" << vertices_.__repr__()
514 << " simplices=" << simplices_.__repr__();
515 }
516
__str__(const String &) const517 String Mesh::__str__(const String & ) const
518 {
519 return __repr__();
520 }
521
522 /* Drawing method */
draw() const523 Graph Mesh::draw() const
524 {
525 if (!(dimension_ <= 3)) throw InvalidArgumentException(HERE) << "Error: cannot draw a mesh of dimension > 3, here dimension=" << dimension_;
526 if (dimension_ == 1) return draw1D();
527 if (dimension_ == 2) return draw2D();
528 if (dimension_ == 3) return draw3D();
529 return Graph();
530 }
531
draw1D() const532 Graph Mesh::draw1D() const
533 {
534 checkValidity();
535 if (dimension_ != 1) throw InvalidArgumentException(HERE) << "Error: cannot draw a mesh of dimension different from 1 with the draw1D() method, here dimension=" << dimension_;
536 const UnsignedInteger verticesSize = getVerticesNumber();
537 const UnsignedInteger simplicesSize = getSimplicesNumber();
538 if (!(verticesSize > 0)) throw InvalidArgumentException(HERE) << "Error: cannot draw a mesh with no vertex.";
539 Graph graph(String(OSS() << "Mesh " << getName()), "", getDescription()[0], true, "topright");
540 // The vertices
541 Cloud vertices(vertices_, Sample(verticesSize, Point(1, 0.0)));
542 vertices.setColor("red");
543 vertices.setLegend(String(OSS() << verticesSize << " node" << (verticesSize > 1 ? "s" : "")));
544 // The simplices
545 for (UnsignedInteger i = 0; i < simplicesSize; ++i)
546 {
547 Sample data(2, 2);
548 data(0, 0) = vertices_(simplices_(i, 0), 0);
549 data(1, 0) = vertices_(simplices_(i, 1), 0);
550 Curve simplex(data);
551 simplex.setColor("blue");
552 if (i == 0) simplex.setLegend(String(OSS() << simplicesSize << " element" << (simplicesSize > 1 ? "s" : "")));
553 graph.add(simplex);
554 }
555 graph.add(vertices);
556 return graph;
557 }
558
draw2D() const559 Graph Mesh::draw2D() const
560 {
561 checkValidity();
562 const UnsignedInteger verticesSize = getVerticesNumber();
563 const UnsignedInteger simplicesSize = getSimplicesNumber();
564 if (!(verticesSize > 0)) throw InvalidArgumentException(HERE) << "Error: cannot draw a mesh with no vertex.";
565 Graph graph(String(OSS() << "Mesh " << getName()), getDescription()[0], getDescription()[1], true, "topright");
566 // The vertices
567 Cloud vertices(vertices_);
568 vertices.setColor("red");
569 if (vertices_.getSize() > ResourceMap::GetAsUnsignedInteger("Mesh-LargeSize")) vertices.setPointStyle("dot");
570 vertices.setLegend(String(OSS() << verticesSize << " node" << (verticesSize > 1 ? "s" : "")));
571 // The simplices
572 for (UnsignedInteger i = 0; i < simplicesSize; ++i)
573 {
574 Sample data(4, 2);
575 data[0] = vertices_[simplices_(i, 0)];
576 data[1] = vertices_[simplices_(i, 1)];
577 data[2] = vertices_[simplices_(i, 2)];
578 data[3] = vertices_[simplices_(i, 0)];
579 Curve simplex(data);
580 simplex.setColor("blue");
581 if (i == 0) simplex.setLegend(String(OSS() << simplicesSize << " element" << (simplicesSize > 1 ? "s" : "")));
582 graph.add(simplex);
583 }
584 graph.add(vertices);
585 return graph;
586 }
587
draw3D(const Bool drawEdge,const Scalar thetaX,const Scalar thetaY,const Scalar thetaZ,const Bool shading,const Scalar rho) const588 Graph Mesh::draw3D(const Bool drawEdge,
589 const Scalar thetaX,
590 const Scalar thetaY,
591 const Scalar thetaZ,
592 const Bool shading,
593 const Scalar rho) const
594 {
595 SquareMatrix R(3);
596 const Scalar sinThetaX = sin(thetaX);
597 const Scalar cosThetaX = cos(thetaX);
598 const Scalar sinThetaY = sin(thetaY);
599 const Scalar cosThetaY = cos(thetaY);
600 const Scalar sinThetaZ = sin(thetaZ);
601 const Scalar cosThetaZ = cos(thetaZ);
602 R(0, 0) = cosThetaY * cosThetaZ;
603 R(1, 0) = -cosThetaY * sinThetaZ;
604 R(2, 0) = sinThetaY;
605 R(0, 1) = cosThetaX * sinThetaZ + sinThetaX * sinThetaY * cosThetaZ;
606 R(1, 1) = cosThetaX * cosThetaZ - sinThetaX * sinThetaY * sinThetaZ;
607 R(2, 1) = -sinThetaX * cosThetaY;
608 R(0, 2) = sinThetaX * sinThetaZ - cosThetaX * sinThetaY * cosThetaZ;
609 R(1, 2) = sinThetaX * cosThetaZ + cosThetaX * sinThetaY * sinThetaZ;
610 R(2, 2) = cosThetaX * cosThetaY;
611 return draw3D(drawEdge, R, shading, rho);
612 }
613
614 namespace
615 {
616 // Check whether a face of a simplex is inner or on a boundary. Arguments must be sorted.
Mesh_isInnerFace(const Indices & simplices0,const Indices & simplices1,const Indices & simplices2)617 Bool Mesh_isInnerFace(const Indices & simplices0, const Indices & simplices1, const Indices & simplices2)
618 {
619 std::vector<UnsignedInteger> common01;
620 std::set_intersection(simplices0.begin(), simplices0.end(), simplices1.begin(), simplices1.end(), std::back_inserter(common01));
621 if (common01.size() < 2) return false;
622 std::vector<UnsignedInteger> common012;
623 std::set_intersection(simplices2.begin(), simplices2.end(), common01.begin(), common01.end(), std::back_inserter(common012));
624 return common012.size() > 1;
625 }
626
627 // Check whether a face of a simplex is oriented toward the front or back
Mesh_isVisible(const Point & visuVertex0,const Point & visuVertex1,const Point & visuVertex2)628 Bool Mesh_isVisible(const Point & visuVertex0, const Point & visuVertex1, const Point & visuVertex2)
629 {
630 return (visuVertex1[0] - visuVertex0[0]) * (visuVertex2[1] - visuVertex0[1]) <= (visuVertex1[1] - visuVertex0[1]) * (visuVertex2[0] - visuVertex0[0]);
631 }
632 } // anonymous namespace
633
draw3D(const Bool drawEdge,const SquareMatrix & rotation,const Bool shading,const Scalar rho) const634 Graph Mesh::draw3D(const Bool drawEdge,
635 const SquareMatrix & rotation,
636 const Bool shading,
637 const Scalar rho) const
638 {
639 checkValidity();
640 // First, check if the matrix is a rotation matrix of R^3
641 if (rotation.getDimension() != 3) throw InvalidArgumentException(HERE) << "Error: the matrix is not a 3d square matrix.";
642 if (!(Point((rotation * rotation.transpose() - IdentityMatrix(3)).getImplementation()).norm() <= 1e-5)) throw InvalidArgumentException(HERE) << "Error: the matrix is not a rotation matrix.";
643 const UnsignedInteger verticesSize = getVerticesNumber();
644 const UnsignedInteger simplicesSize = getSimplicesNumber();
645 if (!(verticesSize > 0)) throw InvalidArgumentException(HERE) << "Error: cannot draw a mesh with no vertex or no simplex.";
646 // We use a basic Painter algorithm for the visualization
647 // Second, transform the vertices if needed
648 const Bool noRotation = rotation.isDiagonal();
649 const Point verticesCenter(noRotation ? Point(3) : vertices_.computeMean());
650 const Sample visuVertices(noRotation ? vertices_ : rotation.getImplementation()->genSampleProd(vertices_ - verticesCenter, true, false, 'R') + verticesCenter);
651
652 // Third, split all the simplices into triangles and compute their mean depth
653 Collection< std::pair<Scalar, Indices> > trianglesAndDepth(4 * simplicesSize);
654 UnsignedInteger triangleIndex = 0;
655 Indices triangle(3);
656 const UnsignedInteger numSimplices = getSimplicesNumber();
657 const UnsignedInteger numVertices = getVerticesNumber();
658 Collection<Indices> mapVerticesToSimplices(numVertices, Indices(0));
659 for (UnsignedInteger i = 0; i < numSimplices; ++i)
660 {
661 for (IndicesCollection::const_iterator cit = simplices_.cbegin_at(i), guard = simplices_.cend_at(i); cit != guard; ++cit)
662 {
663 mapVerticesToSimplices[*cit].add(i);
664 }
665 } // Loop over simplices
666 IndicesCollection verticesToSimplices(mapVerticesToSimplices);
667 const Bool backfaceCulling = ResourceMap::GetAsBool("Mesh-BackfaceCulling");
668 for (UnsignedInteger i = 0; i < simplicesSize; ++i)
669 {
670 const UnsignedInteger i0 = simplices_(i, 0);
671 const UnsignedInteger i1 = simplices_(i, 1);
672 const UnsignedInteger i2 = simplices_(i, 2);
673 const UnsignedInteger i3 = simplices_(i, 3);
674 const Indices simplicesVertex0(verticesToSimplices.cbegin_at(i0), verticesToSimplices.cend_at(i0));
675 const Indices simplicesVertex1(verticesToSimplices.cbegin_at(i1), verticesToSimplices.cend_at(i1));
676 const Indices simplicesVertex2(verticesToSimplices.cbegin_at(i2), verticesToSimplices.cend_at(i2));
677 const Indices simplicesVertex3(verticesToSimplices.cbegin_at(i3), verticesToSimplices.cend_at(i3));
678 const Point visuVertex0(visuVertices[i0]);
679 const Point visuVertex1(visuVertices[i1]);
680 const Point visuVertex2(visuVertices[i2]);
681 const Point visuVertex3(visuVertices[i3]);
682 // First face: AB=p0p1, AC=p0p2.
683 if (((!backfaceCulling) || Mesh_isVisible(visuVertex0, visuVertex1, visuVertex2)) && (!Mesh_isInnerFace(simplicesVertex0, simplicesVertex1, simplicesVertex2)))
684 {
685 triangle = {i0, i1, i2};
686 trianglesAndDepth[triangleIndex].first = visuVertices(i0, 2) + visuVertices(i1, 2) + visuVertices(i2, 2);
687 trianglesAndDepth[triangleIndex].second = triangle;
688 ++triangleIndex;
689 }
690
691 // Second face: AB=p0p2, AC=p0p3.
692 if (((!backfaceCulling) || Mesh_isVisible(visuVertex0, visuVertex2, visuVertex3)) && (!Mesh_isInnerFace(simplicesVertex0, simplicesVertex2, simplicesVertex3)))
693 {
694 triangle[0] = i0;
695 triangle[1] = i2;
696 triangle[2] = i3;
697 trianglesAndDepth[triangleIndex].first = visuVertices(i0, 2) + visuVertices(i2, 2) + visuVertices(i3, 2);
698 trianglesAndDepth[triangleIndex].second = triangle;
699 ++triangleIndex;
700 }
701
702 // Third face: AB=p0p3, AC=p0p1.
703 if (((!backfaceCulling) || Mesh_isVisible(visuVertex0, visuVertex3, visuVertex1)) && (!Mesh_isInnerFace(simplicesVertex0, simplicesVertex3, simplicesVertex1)))
704 {
705 triangle[0] = i0;
706 triangle[1] = i3;
707 triangle[2] = i1;
708 trianglesAndDepth[triangleIndex].first = visuVertices(i0, 2) + visuVertices(i3, 2) + visuVertices(i1, 2);
709 trianglesAndDepth[triangleIndex].second = triangle;
710 ++triangleIndex;
711 }
712
713 // Fourth face: AB=p1p3, AC=p1p2.
714 if (((!backfaceCulling) || Mesh_isVisible(visuVertex1, visuVertex3, visuVertex2)) && (!Mesh_isInnerFace(simplicesVertex1, simplicesVertex3, simplicesVertex2)))
715 {
716 triangle[0] = i1;
717 triangle[1] = i3;
718 triangle[2] = i2;
719 trianglesAndDepth[triangleIndex].first = visuVertices(i1, 2) + visuVertices(i3, 2) + visuVertices(i2, 2);
720 trianglesAndDepth[triangleIndex].second = triangle;
721 ++triangleIndex;
722 }
723 }
724 trianglesAndDepth.resize(triangleIndex);
725
726 // Fourth, draw the triangles in decreasing depth
727 Graph graph(String(OSS() << "Mesh " << getName()), getDescription()[0], getDescription()[1], true, "topright");
728 std::sort(trianglesAndDepth.begin(), trianglesAndDepth.end());
729 Scalar clippedRho = std::min(1.0, std::max(0.0, rho));
730 if (rho != clippedRho) LOGWARN(OSS() << "The shrinking factor must be in (0,1), here rho=" << rho);
731 Sample face(3, 2);
732 Sample data;
733 Description palette;
734 if (!drawEdge)
735 {
736 data = Sample(3 * trianglesAndDepth.getSize(), 2);
737 palette = Description(trianglesAndDepth.getSize());
738 }
739 UnsignedInteger base = 0;
740 UnsignedInteger index = 0;
741
742 const Scalar kSpecular = ResourceMap::GetAsScalar("Mesh-SpecularFactor");
743 const Scalar kDiffuse = ResourceMap::GetAsScalar("Mesh-DiffuseFactor");
744 const Scalar kAmbient = ResourceMap::GetAsScalar("Mesh-AmbientFactor");
745 const Scalar shininess = ResourceMap::GetAsScalar("Mesh-Shininess");
746
747 const Scalar redAmbient = 1.0;
748 const Scalar greenAmbient = 1.0;
749 const Scalar blueAmbient = 0.0;
750 Point Iambient(3);
751 Iambient[0] = kAmbient * redAmbient;
752 Iambient[1] = kAmbient * greenAmbient;
753 Iambient[2] = kAmbient * blueAmbient;
754
755 const Scalar redFace = 0.0;
756 const Scalar greenFace = 0.0;
757 const Scalar blueFace = 1.0;
758
759 const Scalar redEdge = 1.0;
760 const Scalar greenEdge = 0.0;
761 const Scalar blueEdge = 0.0;
762
763 const Scalar redLight = 1.0;
764 const Scalar greenLight = 1.0;
765 const Scalar blueLight = 1.0;
766 Point Ilight(3);
767
768 // Will be modified if shading == true
769 String faceColor = Drawable::ConvertFromRGB(redFace, greenFace, blueFace);
770 String edgeColor = Drawable::ConvertFromRGB(redEdge, greenEdge, blueEdge);
771
772 for (UnsignedInteger i = trianglesAndDepth.getSize(); i > 0; --i)
773 {
774 const UnsignedInteger i0 = trianglesAndDepth[i - 1].second[0];
775 const UnsignedInteger i1 = trianglesAndDepth[i - 1].second[1];
776 const UnsignedInteger i2 = trianglesAndDepth[i - 1].second[2];
777 if (clippedRho < 1.0)
778 {
779 const Point faceCenter((visuVertices[i0] + visuVertices[i1] + visuVertices[i2]) / 3.0);
780 face(0, 0) = faceCenter[0];
781 face(0, 1) = faceCenter[1];
782 face(1, 0) = faceCenter[0];
783 face(1, 1) = faceCenter[1];
784 face(2, 0) = faceCenter[0];
785 face(2, 1) = faceCenter[1];
786 if (clippedRho > 0.0)
787 {
788 face(0, 0) += clippedRho * (visuVertices(i0, 0) - faceCenter[0]);
789 face(0, 1) += clippedRho * (visuVertices(i0, 1) - faceCenter[1]);
790 face(1, 0) += clippedRho * (visuVertices(i1, 0) - faceCenter[0]);
791 face(1, 1) += clippedRho * (visuVertices(i1, 1) - faceCenter[1]);
792 face(2, 0) += clippedRho * (visuVertices(i2, 0) - faceCenter[0]);
793 face(2, 1) += clippedRho * (visuVertices(i2, 1) - faceCenter[1]);
794 }
795 }
796 else
797 {
798 face(0, 0) = visuVertices(i0, 0);
799 face(0, 1) = visuVertices(i0, 1);
800 face(1, 0) = visuVertices(i1, 0);
801 face(1, 1) = visuVertices(i1, 1);
802 face(2, 0) = visuVertices(i2, 0);
803 face(2, 1) = visuVertices(i2, 1);
804 }
805
806 if (shading)
807 {
808 // The light source is behind the observer
809 const Point ab(visuVertices[i1] - visuVertices[i0]);
810 const Point ac(visuVertices[i2] - visuVertices[i0]);
811 Point N(3);
812 // The normal is vect(ab, ac)
813 N[0] = ab[1] * ac[2] - ab[2] * ac[1];
814 N[1] = ab[2] * ac[0] - ab[0] * ac[2];
815 N[2] = ab[0] * ac[1] - ab[1] * ac[0];
816 N /= N.norm();
817 // Flip the normal if it is pointing backward
818 if (N[2] < 0.0) N *= -1.0;
819 const Scalar cosTheta = N[2];
820 // R is a unit vector by construction
821 Point R(N * (2.0 * cosTheta));
822 R[2] -= 1.0;
823 const Scalar cosPhi = std::abs(R[2]);
824 const Scalar Idiffuse = kDiffuse * cosTheta;
825 const Scalar Ispecular = kSpecular * pow(cosPhi, shininess);
826 Ilight[0] = Ispecular * redLight;
827 Ilight[1] = Ispecular * greenLight;
828 Ilight[2] = Ispecular * blueLight;
829 // Face color using Phong model
830 faceColor = Drawable::ConvertFromRGB(Iambient[0] + Idiffuse * redFace + Ilight[0], Iambient[1] + Idiffuse * greenFace + Ilight[1], Iambient[2] + Idiffuse * blueFace + Ilight[2]);
831 edgeColor = Drawable::ConvertFromRGB(Iambient[0] + Idiffuse * redEdge + Ilight[0], Iambient[1] + Idiffuse * greenEdge + Ilight[1], Iambient[2] + Idiffuse * blueEdge + Ilight[2]);
832 } // shading
833 if (drawEdge)
834 {
835 Polygon faceAndEdge(face);
836 faceAndEdge.setColor(faceColor);
837 faceAndEdge.setEdgeColor(edgeColor);
838 graph.add(faceAndEdge);
839 } // drawEdge
840 else
841 {
842 data(base, 0) = face(0, 0);
843 data(base, 1) = face(0, 1);
844 data(base + 1, 0) = face(1, 0);
845 data(base + 1, 1) = face(1, 1);
846 data(base + 2, 0) = face(2, 0);
847 data(base + 2, 1) = face(2, 1);
848 base += 3;
849 palette[index] = faceColor;
850 ++index;
851 } // !drawEdge
852 }
853 if (!drawEdge) graph.add(PolygonArray(data, 3, palette));
854 return graph;
855 }
856
857 /* Import mesh from FreeFem 2D mesh files */
ImportFromMSHFile(const String & fileName)858 Mesh Mesh::ImportFromMSHFile(const String & fileName)
859 {
860 std::ifstream file(fileName.c_str(), std::ios::in);
861 if (!file) throw FileNotFoundException(HERE) << "Error: can't open file " << fileName;
862 // Bording case: empty file
863 if (file.eof())
864 {
865 Log::Info(OSS() << "File " << fileName << " is empty.");
866 return Mesh();
867 }
868 // First, the header: it is made of 3 integers, the number of vertices, the number of simplices and the number of elements on the boundary, currently not used by OT
869 UnsignedInteger verticesNumber = 0;
870 UnsignedInteger simplicesNumber = 0;
871 UnsignedInteger scratch = 0;
872 file >> verticesNumber;
873 file >> simplicesNumber;
874 file >> scratch;
875 LOGINFO(OSS() << "Number of vertices=" << verticesNumber << ", number of simplices=" << simplicesNumber);
876 // Parse the vertices
877 Sample vertices(verticesNumber, 2);
878 for (UnsignedInteger i = 0; i < verticesNumber; ++i)
879 {
880 file >> vertices(i, 0);
881 file >> vertices(i, 1);
882 file >> scratch;
883 LOGINFO(OSS() << "vertex " << i << "=" << vertices[i]);
884 }
885 // Parse the simplices
886 IndicesCollection simplices(simplicesNumber, 3);
887 for (UnsignedInteger i = 0; i < simplicesNumber; ++i)
888 {
889 file >> simplices(i, 0);
890 file >> simplices(i, 1);
891 file >> simplices(i, 2);
892 --simplices(i, 0);
893 --simplices(i, 1);
894 --simplices(i, 2);
895 file >> scratch;
896 LOGINFO(OSS() << "simplex " << i << "=" << simplices(i, 0) << " " << simplices(i, 1) << " " << simplices(i, 2));
897 }
898 file.close();
899 Mesh result(vertices, simplices);
900 result.fixOrientation();
901 return result;
902 }
903
904 /* VTK export */
streamToVTKFormat() const905 String Mesh::streamToVTKFormat() const
906 {
907 return streamToVTKFormat(simplices_);
908 }
909
streamToVTKFormat(const IndicesCollection & simplices) const910 String Mesh::streamToVTKFormat(const IndicesCollection & simplices) const
911 {
912 if (!(dimension_ <= 3)) throw InvalidDimensionException(HERE) << "Error: cannot export a mesh of dimension=" << dimension_ << " into the VTK format. Maximum dimension is 3.";
913 const UnsignedInteger oldPrecision = PlatformInfo::GetNumericalPrecision();
914 PlatformInfo::SetNumericalPrecision(16);
915 OSS oss;
916 // First, the file version and identifier
917 oss << "# vtk DataFile Version 3.0\n";
918 // Second, the header
919 oss << getName() << "\n";
920 // Third, the format
921 oss << "ASCII\n";
922 oss << "\n";
923 // Fourth, the data set
924 oss << "DATASET UNSTRUCTURED_GRID\n";
925 // Fifth, the geometrical and topological data
926 // The vertices
927 const UnsignedInteger numVertices = getVerticesNumber();
928 oss << "POINTS " << numVertices << " float\n";
929 for (UnsignedInteger i = 0; i < numVertices; ++i)
930 {
931 String separator("");
932 for (UnsignedInteger j = 0; j < dimension_; ++j, separator = " ")
933 oss << separator << vertices_(i, j);
934 for (UnsignedInteger j = dimension_; j < 3; ++j)
935 oss << separator << "0.0";
936 oss << "\n";
937 }
938 // The simplices
939 oss << "\n";
940 const UnsignedInteger numSimplices = simplices.getSize();
941 // If no simplex, assume that it is a cloud of points
942 if (numSimplices == 0)
943 {
944 oss << "CELLS " << numVertices << " " << 2 * numVertices << "\n";
945 for (UnsignedInteger i = 0; i < numVertices; ++i) oss << "1 " << i << "\n";
946 oss << "\n";
947 oss << "CELL_TYPES " << numVertices << "\n";
948 for (UnsignedInteger i = 0; i < numVertices; ++i) oss << "1\n";
949 PlatformInfo::SetNumericalPrecision(oldPrecision);
950 return oss;
951 }
952 // There is at least one simplex. Assume homogeneous simplices,
953 // ie all the simplices are of the same kind as the first one
954 UnsignedInteger verticesPerSimplex = 1;
955 UnsignedInteger lastIndex = simplices(0, 0);
956 while ((verticesPerSimplex <= dimension_) && (simplices(0, verticesPerSimplex) != lastIndex))
957 {
958 lastIndex = simplices(0, verticesPerSimplex);
959 ++verticesPerSimplex;
960 }
961 oss << "CELLS " << numSimplices << " " << (verticesPerSimplex + 1) * numSimplices << "\n";
962 for (UnsignedInteger i = 0; i < numSimplices; ++i)
963 {
964 oss << verticesPerSimplex;
965 for (UnsignedInteger j = 0; j < verticesPerSimplex; ++j) oss << " " << simplices(i, j);
966 oss << "\n";
967 }
968 oss << "\n";
969 // If no simplices, assume vertices type
970 oss << "CELL_TYPES " << numSimplices << "\n";
971 UnsignedInteger cellType = 0;
972 // Cell type is:
973 // 1 for vertex
974 // 3 for line
975 // 5 for triangle
976 // 10 for tetrahedron
977 if (verticesPerSimplex == 1) cellType = 1;
978 if (verticesPerSimplex == 2) cellType = 3;
979 if (verticesPerSimplex == 3) cellType = 5;
980 if (verticesPerSimplex == 4) cellType = 10;
981 for (UnsignedInteger i = 0; i < numSimplices; ++i) oss << cellType << "\n";
982 PlatformInfo::SetNumericalPrecision(oldPrecision);
983 return oss;
984 }
985
exportToVTKFile(const String & fileName) const986 void Mesh::exportToVTKFile(const String & fileName) const
987 {
988 exportToVTKFile(fileName, simplices_);
989 }
990
exportToVTKFile(const String & fileName,const IndicesCollection & simplices) const991 void Mesh::exportToVTKFile(const String & fileName,
992 const IndicesCollection & simplices) const
993 {
994 std::ofstream file(fileName.c_str(), std::ios::out);
995 if (!file) throw FileNotFoundException(HERE) << "Error: can't open file " << fileName;
996 const String content(streamToVTKFormat(simplices));
997 file << content;
998 file.close();
999 }
1000
1001 /* Method save() stores the object through the StorageManager */
save(Advocate & adv) const1002 void Mesh::save(Advocate & adv) const
1003 {
1004 PersistentObject::save(adv);
1005 adv.saveAttribute("dimension_", dimension_);
1006 adv.saveAttribute("hasBeenChecked_", hasBeenChecked_);
1007 adv.saveAttribute("vertices_", vertices_);
1008 adv.saveAttribute("simplices_", simplices_);
1009 }
1010
1011 /* Method load() reloads the object from the StorageManager */
load(Advocate & adv)1012 void Mesh::load(Advocate & adv)
1013 {
1014 PersistentObject::load(adv);
1015 adv.loadAttribute("dimension_", dimension_);
1016 adv.loadAttribute("hasBeenChecked_", hasBeenChecked_);
1017 adv.loadAttribute("vertices_", vertices_);
1018 adv.loadAttribute("simplices_", simplices_);
1019 }
1020
1021 END_NAMESPACE_OPENTURNS
1022