1 // Copyright (C) 2007-2024 CEA, EDF
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, or (at your option) any later version.
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 "INTERPKERNELDefines.hxx"
24 #include "TransformedTriangle.hxx"
25 #include "TetraAffineTransform.hxx"
26 #include "InterpolationOptions.hxx"
27 #include "InterpKernelException.hxx"
28 #include "InterpKernelHashMap.hxx"
29 #include "VectorUtils.hxx"
30 #include "MCIdType.hxx"
38 namespace INTERP_KERNEL
40 // Schema according to which the splitting is performed.
41 // Each line represents one tetrahedron. The numbering is as follows :
53 static const int SPLIT_NODES_5[20] = /* WHY not all well oriented ???? */
62 static const int SPLIT_NODES_5_WO[20] = /* WO for well oriented !!! normals of 3 first points are OUTSIDE the TETRA4 */
71 static const int SPLIT_NODES_6[24] = /* WHY all badly oriented ???? */
81 static const int SPLIT_NODES_6_WO[24] = /* WO for well oriented !!! normals of 3 first points are OUTSIDE the TETRA4 */
91 // Each sub-node is the barycenter of 4 other nodes.
92 // For the faces, these are on the original mesh.
93 // For the barycenter, the four face sub-nodes are used.
94 static const int GENERAL_24_SUB_NODES[28] =
96 0,1,4,5,// sub-node 8 (face)
97 0,1,2,3,// sub-node 9 (face)
98 0,3,4,7,// sub-node 10 (face)
99 1,2,5,6,// sub-node 11 (face)
100 4,5,6,7,// sub-node 12 (face)
101 2,3,6,7,// sub-node 13 (face)
102 8,9,10,11// sub-node 14 (cell)
105 static const int GENERAL_24_SUB_NODES_WO[28] =
107 0,4,5,1,// sub-node 8 (face)
108 0,1,2,3,// sub-node 9 (face)
109 0,3,7,4,// sub-node 10 (face)
110 1,5,6,2,// sub-node 11 (face)
111 4,7,6,5,// sub-node 12 (face)
112 2,6,7,3,// sub-node 13 (face)
113 8,9,10,11// sub-node 14 (cell)
116 static const int TETRA_EDGES_GENERAL_24[48] =
118 // face with center 8
123 // face with center 9
128 // face with center 10
133 // face with center 11
138 // face with center 12
143 // face with center 13
150 // Each sub-node is the barycenter of two other nodes.
151 // For the edges, these lie on the original mesh.
152 // For the faces, these are the edge sub-nodes.
153 // For the cell these are two face sub-nodes.
154 static const int GENERAL_48_SUB_NODES[38] =
156 0,1, // sub-node 8 (edge)
157 0,4, // sub-node 9 (edge)
158 1,5, // sub-node 10 (edge)
159 4,5, // sub-node 11 (edge)
160 0,3, // sub-node 12 (edge)
161 1,2, // sub-node 13 (edge)
162 4,7, // sub-node 14 (edge)
163 5,6, // sub-node 15 (edge)
164 2,3, // sub-node 16 (edge)
165 3,7, // sub-node 17 (edge)
166 2,6, // sub-node 18 (edge)
167 6,7, // sub-node 19 (edge)
168 8,11, // sub-node 20 (face)
169 12,13, // sub-node 21 (face)
170 9,17, // sub-node 22 (face)
171 10,18, // sub-node 23 (face)
172 14,15, // sub-node 24 (face)
173 16,19, // sub-node 25 (face)
174 20,25 // sub-node 26 (cell)
177 // Define 8 hexahedral subzones as in Grandy, p449
178 // the values correspond to the nodes that correspond to nodes 1,2,3,4,5,6,7,8 in the subcell
179 // For the correspondence of the nodes, see the GENERAL_48_SUB_NODES table in calculateSubNodes
180 static const int GENERAL_48_SUBZONES[64] =
182 0,8,21,12,9,20,26,22,
183 8,1,13,21,20,10,23,26,
184 12,21,16,3,22,26,25,17,
185 21,13,2,16,26,23,18,25,
186 9,20,26,22,4,11,24,14,
187 20,10,23,26,11,5,15,24,
188 22,26,25,17,14,24,19,7,
189 26,23,18,25,24,15,6,19
192 static const mcIdType GENERAL_48_SUBZONES_2[64] =
194 0,-1,-14,-5,-2,-13,-19,-15,
195 -1,1,-6,-14,-13,-3,-16,-19,
196 -5,-14,-9,3,-15,-19,-18,-10,
197 -14,-6,2,-9,-19,-16,-11,-18,
198 -2,-13,-19,-15,4,-4,-17,-7,
199 -13,-3,-16,-19,-4,5,-8,-17,
200 -15,-19,-18,-10,-7,-17,-12,7,
201 -19,-16,-11,-18,-17,-8,6,-12};
203 void SplitHexa8IntoTetras(SplittingPolicy policy, const mcIdType *nodalConnBg, const mcIdType *nodalConnEnd, const double *coords,
204 std::vector<mcIdType>& tetrasNodalConn, std::vector<double>& addCoords);
206 INTERPKERNEL_EXPORT void SplitIntoTetras(SplittingPolicy policy, NormalizedCellType gt, const mcIdType *nodalConnBg, const mcIdType *nodalConnEnd, const double *coords,
207 std::vector<mcIdType>& tetrasNodalConn, std::vector<double>& addCoords);
210 * \brief Class representing a triangular face, used as key in caching hash map in SplitterTetra.
213 class TriangleFaceKey
219 * Sorts the given nodes (so that the order in which they are passed does not matter) and
220 * calculates a hash value for the key.
222 * @param node1 global number of the first node of the face
223 * @param node2 global number of the second node of the face
224 * @param node3 global number of the third node of the face
226 TriangleFaceKey(mcIdType node1, mcIdType node2, mcIdType node3)
228 Sort3Ints(_nodes, node1, node2, node3);
229 _hashVal = ( _nodes[0] + _nodes[1] + _nodes[2] ) % 29;
233 * Equality comparison operator.
234 * Compares this TriangleFaceKey object to another and determines if they represent the same face.
236 * @param key TriangleFaceKey with which to compare
237 * @return true if key has the same three nodes as this object, false if not
239 bool operator==(const TriangleFaceKey& key) const
241 return _nodes[0] == key._nodes[0] && _nodes[1] == key._nodes[1] && _nodes[2] == key._nodes[2];
245 * Less than operator.
247 * @param key TriangleFaceKey with which to compare
248 * @return true if this object has the three nodes less than the nodes of the key object, false if not
250 bool operator<(const TriangleFaceKey& key) const
252 for (int i = 0; i < 3; ++i)
254 if (_nodes[i] < key._nodes[i])
258 else if (_nodes[i] > key._nodes[i])
267 * Returns a hash value for the object, based on its three nodes.
268 * This value is not unique for each face.
270 * @return a hash value for the object
272 mcIdType hashVal() const
277 inline static void Sort3Ints(mcIdType* sorted, mcIdType node1, mcIdType node2, mcIdType node3);
280 /// global numbers of the three nodes, sorted in ascending order
283 /// hash value for the object, calculated in the constructor
288 * Method to sort three integers in ascending order
290 * @param sorted mcIdType[3] array in which to store the result
291 * @param x1 first integer
292 * @param x2 second integer
293 * @param x3 third integer
295 inline void TriangleFaceKey::Sort3Ints(mcIdType* sorted, mcIdType x1, mcIdType x2, mcIdType x3)
303 sorted[1] = x2 < x3 ? x2 : x3;
304 sorted[2] = x2 < x3 ? x3 : x2;
320 sorted[1] = x1 < x3 ? x1 : x3;
321 sorted[2] = x1 < x3 ? x3 : x1;
334 * \brief Template specialization of INTERP_KERNEL::hash<T> function object for use with a
335 * with TriangleFaceKey as key class.
338 template<> class hash<INTERP_KERNEL::TriangleFaceKey>
342 * Operator() that returns the precalculated hashvalue of a TriangleFaceKey object.
344 * @param key a TriangleFaceKey object
345 * @return an integer hash value for key
347 mcIdType operator()(const INTERP_KERNEL::TriangleFaceKey& key) const
349 return key.hashVal();
354 namespace INTERP_KERNEL
357 * \brief Class calculating the volume of intersection between a tetrahedral target element and
358 * source elements with triangular or quadratilateral faces.
361 template<class MyMeshType>
365 typedef typename MyMeshType::MyConnType ConnType;
367 SplitterTetra(const MyMeshType& srcMesh, const double** tetraCorners, const typename MyMeshType::MyConnType *nodesId);
369 SplitterTetra(const MyMeshType& srcMesh, const double tetraCorners[12], const ConnType *conn = 0);
373 double intersectSourceCell(typename MyMeshType::MyConnType srcCell, double* baryCentre=0);
374 double intersectSourceFace(const NormalizedCellType polyType,
375 const ConnType polyNodesNbr,
376 const ConnType *const polyNodes,
377 const double *const *const polyCoords,
378 const double dimCaracteristic,
379 const double precision,
380 std::multiset<TriangleFaceKey>& listOfTetraFacesTreated,
381 std::set<TriangleFaceKey>& listOfTetraFacesColinear);
383 double intersectTetra(const double** tetraCorners);
385 ConnType getId(mcIdType id) { return _conn[id]; }
387 void splitIntoDualCells(SplitterTetra<MyMeshType> **output);
389 void splitMySelfForDual(double* output, int i, ConnType& nodeId);
391 void clearVolumesCache();
394 inline static void CheckIsOutside(const double* pt, bool* isOutside, const double errTol = DEFAULT_ABS_TOL);
395 inline static void CheckIsStrictlyOutside(const double* pt, bool* isStrictlyOutside, const double errTol = DEFAULT_ABS_TOL);
396 inline void calculateNode(ConnType globalNodeNum);
397 inline void calculateNode2(ConnType globalNodeNum, const double* node);
398 inline void calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key);
399 inline void calculateSurface(TransformedTriangle& tri, const TriangleFaceKey& key);
401 static inline bool IsFacesCoplanar(const double *const planeNormal, const double planeConstant,
402 const double *const *const coordsFace, const double precision);
403 static inline double CalculateIntersectionSurfaceOfCoplanarTriangles(const double *const planeNormal,
404 const double planeConstant,
405 const double *const p1, const double *const p2, const double *const p3,
406 const double *const p4, const double *const p5, const double *const p6,
407 const double dimCaracteristic, const double precision);
410 SplitterTetra(const SplitterTetra& t);
412 /// disallow assignment
413 SplitterTetra& operator=(const SplitterTetra& t);
416 /// affine transform associated with this target element
417 TetraAffineTransform* _t;
419 /// HashMap relating node numbers to transformed nodes, used for caching
420 HashMap< ConnType , double* > _nodes;
422 /// HashMap relating triangular faces to calculated volume contributions, used for caching
423 HashMap< TriangleFaceKey, double > _volumes;
425 /// reference to the source mesh
426 const MyMeshType& _src_mesh;
428 // node id of the first node in target mesh in C mode.
433 /// Smallest volume of the intersecting elements in the transformed space that will be returned as non-zero.
434 /// 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.
435 static const double SPARSE_TRUNCATION_LIMIT;
439 * Function used to filter out elements by checking if they belong to one of the halfspaces
440 * x <= 0, x >= 1, y <= 0, y >= 1, z <= 0, z >= 1, (indexed 0 - 7). The function updates an array of boolean variables
441 * which indicates whether the points that have already been checked are all in a halfspace. For each halfspace,
442 * 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.
444 * @param pt double[3] containing the coordinates of a transformed point
445 * @param isOutside bool[8] which indicate the results of earlier checks.
447 template<class MyMeshType>
448 inline void SplitterTetra<MyMeshType>::CheckIsOutside(const double* pt, bool* isOutside, const double errTol)
450 isOutside[0] = isOutside[0] && (pt[0] < errTol);
451 isOutside[1] = isOutside[1] && (pt[0] > (1.0-errTol) );
452 isOutside[2] = isOutside[2] && (pt[1] < errTol);
453 isOutside[3] = isOutside[3] && (pt[1] > (1.0-errTol));
454 isOutside[4] = isOutside[4] && (pt[2] < errTol);
455 isOutside[5] = isOutside[5] && (pt[2] > (1.0-errTol));
456 isOutside[6] = isOutside[6] && (1.0 - pt[0] - pt[1] - pt[2] < errTol);
457 isOutside[7] = isOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] > (1.0-errTol) );
460 template<class MyMeshType>
461 inline void SplitterTetra<MyMeshType>::CheckIsStrictlyOutside(const double* pt, bool* isStrictlyOutside, const double errTol)
463 isStrictlyOutside[0] = isStrictlyOutside[0] && (pt[0] < -errTol);
464 isStrictlyOutside[1] = isStrictlyOutside[1] && (pt[0] > (1.0 + errTol));
465 isStrictlyOutside[2] = isStrictlyOutside[2] && (pt[1] < -errTol);
466 isStrictlyOutside[3] = isStrictlyOutside[3] && (pt[1] > (1.0 + errTol));
467 isStrictlyOutside[4] = isStrictlyOutside[4] && (pt[2] < -errTol);
468 isStrictlyOutside[5] = isStrictlyOutside[5] && (pt[2] > (1.0 + errTol));
469 isStrictlyOutside[6] = isStrictlyOutside[6] && (1.0 - pt[0] - pt[1] - pt[2] < -errTol);
470 isStrictlyOutside[7] = isStrictlyOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] > (1.0 + errTol));
474 * Calculates the transformed node with a given global node number.
475 * Gets the coordinates for the node in _src_mesh with the given global number and applies TetraAffineTransform
476 * _t to it. Stores the result in the cache _nodes. The non-existance of the node in _nodes should be verified before
479 * @param globalNodeNum global node number of the node in the mesh _src_mesh
482 template<class MyMeshType>
483 inline void SplitterTetra<MyMeshType>::calculateNode(typename MyMeshType::MyConnType globalNodeNum)
485 const double* node = _src_mesh.getCoordinatesPtr()+MyMeshType::MY_SPACEDIM*globalNodeNum;
486 double* transformedNode = new double[MyMeshType::MY_SPACEDIM];
487 assert(transformedNode != 0);
488 _t->apply(transformedNode, node);
489 _nodes[globalNodeNum] = transformedNode;
494 * Calculates the transformed node with a given global node number.
495 * Applies TetraAffineTransform * _t to it.
496 * Stores the result in the cache _nodes. The non-existence of the node in _nodes should be verified before * calling.
497 * The only difference with the previous method calculateNode is that the coordinates of the node are passed in arguments
498 * and are not recalculated in order to optimize the method.
500 * @param globalNodeNum global node number of the node in the mesh _src_mesh
503 template<class MyMeshType>
504 inline void SplitterTetra<MyMeshType>::calculateNode2(typename MyMeshType::MyConnType globalNodeNum, const double* node)
506 double* transformedNode = new double[MyMeshType::MY_SPACEDIM];
507 assert(transformedNode != 0);
508 _t->apply(transformedNode, node);
509 _nodes[globalNodeNum] = transformedNode;
513 * Calculates the volume contribution from the given TransformedTriangle and stores it with the given key in .
514 * Calls TransformedTriangle::calculateIntersectionVolume to perform the calculation.
516 * @param tri triangle for which to calculate the volume contribution
517 * @param key key associated with the face
519 template<class MyMeshType>
520 inline void SplitterTetra<MyMeshType>::calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key)
522 const double vol = tri.calculateIntersectionVolume();
523 _volumes.insert(std::make_pair(key, vol));
527 * Calculates the surface contribution from the given TransformedTriangle and stores it with the given key in.
528 * Calls TransformedTriangle::calculateIntersectionSurface to perform the calculation.
530 * @param tri triangle for which to calculate the surface contribution
531 * @param key key associated with the face
533 template<class MyMeshType>
534 inline void SplitterTetra<MyMeshType>::calculateSurface(TransformedTriangle& tri, const TriangleFaceKey& key)
536 const double surf = tri.calculateIntersectionSurface(_t);
537 _volumes.insert(std::make_pair(key, surf));
540 template<class MyMeshTypeT, class MyMeshTypeS=MyMeshTypeT>
544 SplitterTetra2(const MyMeshTypeT& targetMesh, const MyMeshTypeS& srcMesh, SplittingPolicy policy);
546 void releaseArrays();
547 SplittingPolicy getSplittingPolicy() const { return _splitting_pol; }
548 void splitTargetCell2(typename MyMeshTypeT::MyConnType targetCell, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
549 void splitTargetCell(typename MyMeshTypeT::MyConnType targetCell, typename MyMeshTypeT::MyConnType nbOfNodesT,
550 typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);//to suppress
551 void fiveSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);//to suppress
552 void sixSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);//to suppress
553 void calculateGeneral24Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);//to suppress
554 void calculateGeneral24TetraOld(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
555 void calculateGeneral48Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);//to suppress
556 void splitPyram5(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);//to suppress
557 void splitConvex(typename MyMeshTypeT::MyConnType targetCell,//to suppress
558 typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);//to suppress
559 void calculateSubNodes(const MyMeshTypeT& targetMesh, typename MyMeshTypeT::MyConnType targetCell);//to suppress
560 inline const double* getCoordsOfSubNode(typename MyMeshTypeT::MyConnType node);//to suppress
561 inline const double* getCoordsOfSubNode2(typename MyMeshTypeT::MyConnType node, typename MyMeshTypeT::MyConnType& nodeId);//to suppress
563 inline void calcBarycenter(typename MyMeshTypeT::MyConnType n, double* barycenter, const int* pts);//to suppress
565 void sixSplitGen(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra, std::function<void(SplitterTetra2& , typename MyMeshTypeS::MyConnType&, const double*&)> func);
566 void calculateGeneral24TetraGen(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra, std::function<void(SplitterTetra2& , typename MyMeshTypeS::MyConnType[4], const double*[4])> func);
568 const MyMeshTypeT& _target_mesh;
569 const MyMeshTypeS& _src_mesh;
570 SplittingPolicy _splitting_pol;
571 /// vector of pointers to double[3] containing the coordinates of the
572 /// (sub) - nodes of split target cell
573 std::vector<const double*> _nodes;
574 std::vector<typename MyMeshTypeT::MyConnType> _node_ids;
578 * Calculates the barycenter of n (sub) - nodes
580 * @param n number of nodes for which to calculate barycenter
581 * @param barycenter pointer to double[3] array in which to store the result
582 * @param pts pointer to int[n] array containing the (sub)-nodes for which to calculate the barycenter
584 template<class MyMeshTypeT, class MyMeshTypeS>
586 inline void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calcBarycenter(typename MyMeshTypeT::MyConnType n, double* barycenter, const int* pts)
588 barycenter[0] = barycenter[1] = barycenter[2] = 0.0;
589 for(int i = 0; i < n ; ++i)
591 const double* pt = getCoordsOfSubNode(pts[i]);
592 barycenter[0] += pt[0];
593 barycenter[1] += pt[1];
594 barycenter[2] += pt[2];
597 barycenter[0] /= (double)n;
598 barycenter[1] /= (double)n;
599 barycenter[2] /= (double)n;
603 * Accessor to the coordinates of a given (sub)-node
605 * @param node local number of the (sub)-node 0,..,7 are the elements nodes, sub-nodes are numbered from 8,..
606 * @return pointer to double[3] containing the coordinates of the nodes
608 template<class MyMeshTypeT, class MyMeshTypeS>
609 inline const double* SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::getCoordsOfSubNode(typename MyMeshTypeT::MyConnType node)
611 // replace "at()" with [] for unsafe but faster access
612 return _nodes.at(node);
616 * Accessor to the coordinates of a given (sub)-node
618 * @param node local number of the (sub)-node 0,..,7 are the elements nodes, sub-nodes are numbered from 8,..
619 * @param nodeId is an output that is node id in target whole mesh in C mode.
620 * @return pointer to double[3] containing the coordinates of the nodes
622 template<class MyMeshTypeT, class MyMeshTypeS>
623 const double* SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::getCoordsOfSubNode2(typename MyMeshTypeT::MyConnType node, typename MyMeshTypeT::MyConnType& nodeId)
625 const double *ret=_nodes.at(node);
627 nodeId=_node_ids[node];