1 // Copyright (C) 2007-2013 CEA/DEN, EDF R&D
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License.
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 // Lesser General Public License for more details.
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
20 #ifndef __SPLITTERTETRA_HXX__
21 #define __SPLITTERTETRA_HXX__
23 #include "TransformedTriangle.hxx"
24 #include "TetraAffineTransform.hxx"
25 #include "InterpolationOptions.hxx"
26 #include "InterpKernelHashMap.hxx"
27 #include "VectorUtils.hxx"
35 namespace INTERP_KERNEL
37 // Schema according to which the splitting is performed.
38 // Each line represents one tetrahedron. The numbering is as follows :
50 static const int SPLIT_NODES_5[20] = /* WHY not all well oriented ???? */
59 static const int SPLIT_NODES_5_WO[20] = /* WO for well oriented !!! normals of 3 first points are OUTSIDE the TETRA4 */
68 static const int SPLIT_NODES_6[24] = /* WHY all badly oriented ???? */
78 static const int SPLIT_NODES_6_WO[24] = /* WO for well oriented !!! normals of 3 first points are OUTSIDE the TETRA4 */
88 // Each sub-node is the barycenter of 4 other nodes.
89 // For the faces, these are on the orignal mesh.
90 // For the barycenter, the four face sub-nodes are used.
91 static const int GENERAL_24_SUB_NODES[28] =
93 0,1,4,5,// sub-node 8 (face)
94 0,1,2,3,// sub-node 9 (face)
95 0,3,4,7,// sub-node 10 (face)
96 1,2,5,6,// sub-node 11 (face)
97 4,5,6,7,// sub-node 12 (face)
98 2,3,6,7,// sub-node 13 (face)
99 8,9,10,11// sub-node 14 (cell)
102 static const int TETRA_EDGES_GENERAL_24[48] =
104 // face with center 8
109 // face with center 9
114 // face with center 10
119 // face with center 11
124 // face with center 12
129 // face with center 13
137 * \brief Class representing a triangular face, used as key in caching hash map in SplitterTetra.
140 class TriangleFaceKey
146 * Sorts the given nodes (so that the order in which they are passed does not matter) and
147 * calculates a hash value for the key.
149 * @param node1 global number of the first node of the face
150 * @param node2 global number of the second node of the face
151 * @param node3 global number of the third node of the face
153 TriangleFaceKey(int node1, int node2, int node3)
155 sort3Ints(_nodes, node1, node2, node3);
156 _hashVal = ( _nodes[0] + _nodes[1] + _nodes[2] ) % 29;
160 * Equality comparison operator.
161 * Compares this TriangleFaceKey object to another and determines if they represent the same face.
163 * @param key TriangleFaceKey with which to compare
164 * @return true if key has the same three nodes as this object, false if not
166 bool operator==(const TriangleFaceKey& key) const
168 return _nodes[0] == key._nodes[0] && _nodes[1] == key._nodes[1] && _nodes[2] == key._nodes[2];
172 * Less than operator.
174 * @param key TriangleFaceKey with which to compare
175 * @return true if this object has the three nodes less than the nodes of the key object, false if not
177 bool operator<(const TriangleFaceKey& key) const
179 for (int i = 0; i < 3; ++i)
181 if (_nodes[i] < key._nodes[i])
185 else if (_nodes[i] > key._nodes[i])
194 * Returns a hash value for the object, based on its three nodes.
195 * This value is not unique for each face.
197 * @return a hash value for the object
204 inline void sort3Ints(int* sorted, int node1, int node2, int node3);
207 /// global numbers of the three nodes, sorted in ascending order
210 /// hash value for the object, calculated in the constructor
215 * Method to sort three integers in ascending order
217 * @param sorted int[3] array in which to store the result
218 * @param x1 first integer
219 * @param x2 second integer
220 * @param x3 third integer
222 inline void TriangleFaceKey::sort3Ints(int* sorted, int x1, int x2, int x3)
230 sorted[1] = x2 < x3 ? x2 : x3;
231 sorted[2] = x2 < x3 ? x3 : x2;
247 sorted[1] = x1 < x3 ? x1 : x3;
248 sorted[2] = x1 < x3 ? x3 : x1;
261 * \brief Template specialization of INTERP_KERNEL::hash<T> function object for use with a
262 * with TriangleFaceKey as key class.
265 template<> class hash<INTERP_KERNEL::TriangleFaceKey>
269 * Operator() that returns the precalculated hashvalue of a TriangleFaceKey object.
271 * @param key a TriangleFaceKey object
272 * @return an integer hash value for key
274 int operator()(const INTERP_KERNEL::TriangleFaceKey& key) const
276 return key.hashVal();
281 namespace INTERP_KERNEL
285 * \brief Class calculating the volume of intersection between a tetrahedral target element and
286 * source elements with triangular or quadratilateral faces.
289 template<class MyMeshType>
294 SplitterTetra(const MyMeshType& srcMesh, const double** tetraCorners, const typename MyMeshType::MyConnType *nodesId);
298 double intersectSourceCell(typename MyMeshType::MyConnType srcCell, double* baryCentre=0);
299 double intersectSourceFace(const NormalizedCellType polyType,
300 const int polyNodesNbr,
301 const int *const polyNodes,
302 const double *const *const polyCoords,
303 const double dimCaracteristic,
304 const double precision,
305 std::multiset<TriangleFaceKey>& listOfTetraFacesTreated,
306 std::set<TriangleFaceKey>& listOfTetraFacesColinear);
308 double intersectTetra(const double** tetraCorners);
310 typename MyMeshType::MyConnType getId(int id) { return _conn[id]; }
312 void splitIntoDualCells(SplitterTetra<MyMeshType> **output);
314 void splitMySelfForDual(double* output, int i, typename MyMeshType::MyConnType& nodeId);
316 void clearVolumesCache();
320 inline void createAffineTransform(const double** corners);
321 inline void checkIsOutside(const double* pt, bool* isOutside, const double errTol = DEFAULT_ABS_TOL) const;
322 inline void checkIsStrictlyOutside(const double* pt, bool* isStrictlyOutside, const double errTol = DEFAULT_ABS_TOL) const;
323 inline void calculateNode(typename MyMeshType::MyConnType globalNodeNum);
324 inline void calculateNode2(typename MyMeshType::MyConnType globalNodeNum, const double* node);
325 inline void calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key);
326 inline void calculateSurface(TransformedTriangle& tri, const TriangleFaceKey& key);
328 static inline bool IsFacesCoplanar(const double *const planeNormal, const double planeConstant,
329 const double *const *const coordsFace, const double precision);
330 static inline double CalculateIntersectionSurfaceOfCoplanarTriangles(const double *const planeNormal,
331 const double planeConstant,
332 const double *const p1, const double *const p2, const double *const p3,
333 const double *const p4, const double *const p5, const double *const p6,
334 const double dimCaracteristic, const double precision);
337 SplitterTetra(const SplitterTetra& t);
339 /// disallow assignment
340 SplitterTetra& operator=(const SplitterTetra& t);
343 /// affine transform associated with this target element
344 TetraAffineTransform* _t;
346 /// HashMap relating node numbers to transformed nodes, used for caching
347 HashMap< int , double* > _nodes;
349 /// HashMap relating triangular faces to calculated volume contributions, used for caching
350 HashMap< TriangleFaceKey, double
352 // , hash_compare<TriangleFaceKey,TriangleFaceKeyComparator>
356 /// reference to the source mesh
357 const MyMeshType& _src_mesh;
359 // node id of the first node in target mesh in C mode.
360 typename MyMeshType::MyConnType _conn[4];
364 /// Smallest volume of the intersecting elements in the transformed space that will be returned as non-zero.
365 /// Since the scale is always the same in the transformed space (the target tetrahedron is unitary), this number is independent of the scale of the meshes.
366 static const double SPARSE_TRUNCATION_LIMIT;
370 * Creates the affine transform _t from the corners of the tetrahedron. Used by the constructors
372 * @param corners double*[4] array containing pointers to four double[3] arrays with the
373 * coordinates of the corners of the tetrahedron
375 template<class MyMeshType>
376 inline void SplitterTetra<MyMeshType>::createAffineTransform(const double** corners)
378 // create AffineTransform from tetrahedron
379 _t = new TetraAffineTransform( corners );
383 * Function used to filter out elements by checking if they belong to one of the halfspaces
384 * x <= 0, x >= 1, y <= 0, y >= 1, z <= 0, z >= 1, (indexed 0 - 7). The function updates an array of boolean variables
385 * which indicates whether the points that have already been checked are all in a halfspace. For each halfspace,
386 * the corresponding array element will be true if and only if it was true when the method was called and pt lies in the halfspace.
388 * @param pt double[3] containing the coordiantes of a transformed point
389 * @param isOutside bool[8] which indicate the results of earlier checks.
391 template<class MyMeshType>
392 inline void SplitterTetra<MyMeshType>::checkIsOutside(const double* pt, bool* isOutside, const double errTol) const
394 isOutside[0] = isOutside[0] && (pt[0] < errTol);
395 isOutside[1] = isOutside[1] && (pt[0] > (1.0-errTol) );
396 isOutside[2] = isOutside[2] && (pt[1] < errTol);
397 isOutside[3] = isOutside[3] && (pt[1] > (1.0-errTol));
398 isOutside[4] = isOutside[4] && (pt[2] < errTol);
399 isOutside[5] = isOutside[5] && (pt[2] > (1.0-errTol));
400 isOutside[6] = isOutside[6] && (1.0 - pt[0] - pt[1] - pt[2] < errTol);
401 isOutside[7] = isOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] > (1.0-errTol) );
404 template<class MyMeshType>
405 inline void SplitterTetra<MyMeshType>::checkIsStrictlyOutside(const double* pt, bool* isStrictlyOutside, const double errTol) const
407 isStrictlyOutside[0] = isStrictlyOutside[0] && (pt[0] < -errTol);
408 isStrictlyOutside[1] = isStrictlyOutside[1] && (pt[0] > (1.0 + errTol));
409 isStrictlyOutside[2] = isStrictlyOutside[2] && (pt[1] < -errTol);
410 isStrictlyOutside[3] = isStrictlyOutside[3] && (pt[1] > (1.0 + errTol));
411 isStrictlyOutside[4] = isStrictlyOutside[4] && (pt[2] < -errTol);
412 isStrictlyOutside[5] = isStrictlyOutside[5] && (pt[2] > (1.0 + errTol));
413 isStrictlyOutside[6] = isStrictlyOutside[6] && (1.0 - pt[0] - pt[1] - pt[2] < -errTol);
414 isStrictlyOutside[7] = isStrictlyOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] > (1.0 + errTol));
418 * Calculates the transformed node with a given global node number.
419 * Gets the coordinates for the node in _src_mesh with the given global number and applies TetraAffineTransform
420 * _t to it. Stores the result in the cache _nodes. The non-existance of the node in _nodes should be verified before
423 * @param globalNodeNum global node number of the node in the mesh _src_mesh
426 template<class MyMeshType>
427 inline void SplitterTetra<MyMeshType>::calculateNode(typename MyMeshType::MyConnType globalNodeNum)
429 const double* node = _src_mesh.getCoordinatesPtr()+MyMeshType::MY_SPACEDIM*globalNodeNum;
430 double* transformedNode = new double[MyMeshType::MY_SPACEDIM];
431 assert(transformedNode != 0);
432 _t->apply(transformedNode, node);
433 _nodes[globalNodeNum] = transformedNode;
438 * Calculates the transformed node with a given global node number.
439 * Applies TetraAffineTransform * _t to it.
440 * Stores the result in the cache _nodes. The non-existence of the node in _nodes should be verified before * calling.
441 * The only difference with the previous method calculateNode is that the coordinates of the node are passed in arguments
442 * and are not recalculated in order to optimize the method.
444 * @param globalNodeNum global node number of the node in the mesh _src_mesh
447 template<class MyMeshType>
448 inline void SplitterTetra<MyMeshType>::calculateNode2(typename MyMeshType::MyConnType globalNodeNum, const double* node)
450 double* transformedNode = new double[MyMeshType::MY_SPACEDIM];
451 assert(transformedNode != 0);
452 _t->apply(transformedNode, node);
453 _nodes[globalNodeNum] = transformedNode;
457 * Calculates the volume contribution from the given TransformedTriangle and stores it with the given key in .
458 * Calls TransformedTriangle::calculateIntersectionVolume to perform the calculation.
460 * @param tri triangle for which to calculate the volume contribution
461 * @param key key associated with the face
463 template<class MyMeshType>
464 inline void SplitterTetra<MyMeshType>::calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key)
466 const double vol = tri.calculateIntersectionVolume();
467 _volumes.insert(std::make_pair(key, vol));
471 * Calculates the surface contribution from the given TransformedTriangle and stores it with the given key in.
472 * Calls TransformedTriangle::calculateIntersectionSurface to perform the calculation.
474 * @param tri triangle for which to calculate the surface contribution
475 * @param key key associated with the face
477 template<class MyMeshType>
478 inline void SplitterTetra<MyMeshType>::calculateSurface(TransformedTriangle& tri, const TriangleFaceKey& key)
480 const double surf = tri.calculateIntersectionSurface(_t);
481 _volumes.insert(std::make_pair(key, surf));
484 template<class MyMeshTypeT, class MyMeshTypeS=MyMeshTypeT>
488 SplitterTetra2(const MyMeshTypeT& targetMesh, const MyMeshTypeS& srcMesh, SplittingPolicy policy);
490 void releaseArrays();
491 void splitTargetCell(typename MyMeshTypeT::MyConnType targetCell, typename MyMeshTypeT::MyConnType nbOfNodesT,
492 typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
493 void fiveSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
494 void sixSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
495 void calculateGeneral24Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
496 void calculateGeneral48Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
497 void splitPyram5(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
498 void splitConvex(typename MyMeshTypeT::MyConnType targetCell,
499 typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
500 void calculateSubNodes(const MyMeshTypeT& targetMesh, typename MyMeshTypeT::MyConnType targetCell);
501 inline const double* getCoordsOfSubNode(typename MyMeshTypeT::MyConnType node);
502 inline const double* getCoordsOfSubNode2(typename MyMeshTypeT::MyConnType node, typename MyMeshTypeT::MyConnType& nodeId);
504 inline void calcBarycenter(int n, double* barycenter, const typename MyMeshTypeT::MyConnType* pts);
506 const MyMeshTypeT& _target_mesh;
507 const MyMeshTypeS& _src_mesh;
508 SplittingPolicy _splitting_pol;
509 /// vector of pointers to double[3] containing the coordinates of the
510 /// (sub) - nodes of split target cell
511 std::vector<const double*> _nodes;
512 std::vector<typename MyMeshTypeT::MyConnType> _node_ids;
516 * Calculates the barycenter of n (sub) - nodes
518 * @param n number of nodes for which to calculate barycenter
519 * @param barycenter pointer to double[3] array in which to store the result
520 * @param pts pointer to int[n] array containing the (sub)-nodes for which to calculate the barycenter
522 template<class MyMeshTypeT, class MyMeshTypeS>
524 inline void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calcBarycenter(int n, double* barycenter, const typename MyMeshTypeT::MyConnType* pts)
526 barycenter[0] = barycenter[1] = barycenter[2] = 0.0;
527 for(int i = 0; i < n ; ++i)
529 const double* pt = getCoordsOfSubNode(pts[i]);
530 barycenter[0] += pt[0];
531 barycenter[1] += pt[1];
532 barycenter[2] += pt[2];
541 * Accessor to the coordinates of a given (sub)-node
543 * @param node local number of the (sub)-node 0,..,7 are the elements nodes, sub-nodes are numbered from 8,..
544 * @return pointer to double[3] containing the coordinates of the nodes
546 template<class MyMeshTypeT, class MyMeshTypeS>
547 inline const double* SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::getCoordsOfSubNode(typename MyMeshTypeT::MyConnType node)
549 // replace "at()" with [] for unsafe but faster access
550 return _nodes.at(node);
554 * Accessor to the coordinates of a given (sub)-node
556 * @param node local number of the (sub)-node 0,..,7 are the elements nodes, sub-nodes are numbered from 8,..
557 * @param nodeId is an output that is node id in target whole mesh in C mode.
558 * @return pointer to double[3] containing the coordinates of the nodes
560 template<class MyMeshTypeT, class MyMeshTypeS>
561 const double* SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::getCoordsOfSubNode2(typename MyMeshTypeT::MyConnType node, typename MyMeshTypeT::MyConnType& nodeId)
563 const double *ret=_nodes.at(node);
565 nodeId=_node_ids[node];