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 "InterpKernelException.hxx"
27 #include "InterpKernelHashMap.hxx"
28 #include "VectorUtils.hxx"
36 namespace INTERP_KERNEL
38 // Schema according to which the splitting is performed.
39 // Each line represents one tetrahedron. The numbering is as follows :
51 static const int SPLIT_NODES_5[20] = /* WHY not all well oriented ???? */
60 static const int SPLIT_NODES_5_WO[20] = /* WO for well oriented !!! normals of 3 first points are OUTSIDE the TETRA4 */
69 static const int SPLIT_NODES_6[24] = /* WHY all badly oriented ???? */
79 static const int SPLIT_NODES_6_WO[24] = /* WO for well oriented !!! normals of 3 first points are OUTSIDE the TETRA4 */
89 // Each sub-node is the barycenter of 4 other nodes.
90 // For the faces, these are on the orignal mesh.
91 // For the barycenter, the four face sub-nodes are used.
92 static const int GENERAL_24_SUB_NODES[28] =
94 0,1,4,5,// sub-node 8 (face)
95 0,1,2,3,// sub-node 9 (face)
96 0,3,4,7,// sub-node 10 (face)
97 1,2,5,6,// sub-node 11 (face)
98 4,5,6,7,// sub-node 12 (face)
99 2,3,6,7,// sub-node 13 (face)
100 8,9,10,11// sub-node 14 (cell)
103 static const int GENERAL_24_SUB_NODES_WO[28] =
105 0,4,5,1,// sub-node 8 (face)
106 0,1,2,3,// sub-node 9 (face)
107 0,3,7,4,// sub-node 10 (face)
108 1,5,6,2,// sub-node 11 (face)
109 4,7,6,5,// sub-node 12 (face)
110 2,6,7,3,// sub-node 13 (face)
111 8,9,10,11// sub-node 14 (cell)
114 static const int TETRA_EDGES_GENERAL_24[48] =
116 // face with center 8
121 // face with center 9
126 // face with center 10
131 // face with center 11
136 // face with center 12
141 // face with center 13
148 // Each sub-node is the barycenter of two other nodes.
149 // For the edges, these lie on the original mesh.
150 // For the faces, these are the edge sub-nodes.
151 // For the cell these are two face sub-nodes.
152 static const int GENERAL_48_SUB_NODES[38] =
154 0,1, // sub-node 8 (edge)
155 0,4, // sub-node 9 (edge)
156 1,5, // sub-node 10 (edge)
157 4,5, // sub-node 11 (edge)
158 0,3, // sub-node 12 (edge)
159 1,2, // sub-node 13 (edge)
160 4,7, // sub-node 14 (edge)
161 5,6, // sub-node 15 (edge)
162 2,3, // sub-node 16 (edge)
163 3,7, // sub-node 17 (edge)
164 2,6, // sub-node 18 (edge)
165 6,7, // sub-node 19 (edge)
166 8,11, // sub-node 20 (face)
167 12,13, // sub-node 21 (face)
168 9,17, // sub-node 22 (face)
169 10,18, // sub-node 23 (face)
170 14,15, // sub-node 24 (face)
171 16,19, // sub-node 25 (face)
172 20,25 // sub-node 26 (cell)
175 // Define 8 hexahedral subzones as in Grandy, p449
176 // the values correspond to the nodes that correspond to nodes 1,2,3,4,5,6,7,8 in the subcell
177 // For the correspondance of the nodes, see the GENERAL_48_SUB_NODES table in calculateSubNodes
178 static const int GENERAL_48_SUBZONES[64] =
180 0,8,21,12,9,20,26,22,
181 8,1,13,21,20,10,23,26,
182 12,21,16,3,22,26,25,17,
183 21,13,2,16,26,23,18,25,
184 9,20,26,22,4,11,24,14,
185 20,10,23,26,11,5,15,24,
186 22,26,25,17,14,24,19,7,
187 26,23,18,25,24,15,6,19
190 static const int GENERAL_48_SUBZONES_2[64] =
192 0,-1,-14,-5,-2,-13,-19,-15,
193 -1,1,-6,-14,-13,-3,-16,-19,
194 -5,-14,-9,3,-15,-19,-18,-10,
195 -14,-6,2,-9,-19,-16,-11,-18,
196 -2,-13,-19,-15,4,-4,-17,-7,
197 -13,-3,-16,-19,-4,5,-8,-17,
198 -15,-19,-18,-10,-7,-17,-12,7,
199 -19,-16,-11,-18,-17,-8,6,-12};
201 void SplitHexa8IntoTetras(SplittingPolicy policy, const int *nodalConnBg, const int *nodalConnEnd, const double *coords,
202 std::vector<int>& tetrasNodalConn, std::vector<double>& addCoords) throw(INTERP_KERNEL::Exception);
204 void SplitIntoTetras(SplittingPolicy policy, NormalizedCellType gt, const int *nodalConnBg, const int *nodalConnEnd, const double *coords,
205 std::vector<int>& tetrasNodalConn, std::vector<double>& addCoords) throw(INTERP_KERNEL::Exception);
208 * \brief Class representing a triangular face, used as key in caching hash map in SplitterTetra.
211 class TriangleFaceKey
217 * Sorts the given nodes (so that the order in which they are passed does not matter) and
218 * calculates a hash value for the key.
220 * @param node1 global number of the first node of the face
221 * @param node2 global number of the second node of the face
222 * @param node3 global number of the third node of the face
224 TriangleFaceKey(int node1, int node2, int node3)
226 Sort3Ints(_nodes, node1, node2, node3);
227 _hashVal = ( _nodes[0] + _nodes[1] + _nodes[2] ) % 29;
231 * Equality comparison operator.
232 * Compares this TriangleFaceKey object to another and determines if they represent the same face.
234 * @param key TriangleFaceKey with which to compare
235 * @return true if key has the same three nodes as this object, false if not
237 bool operator==(const TriangleFaceKey& key) const
239 return _nodes[0] == key._nodes[0] && _nodes[1] == key._nodes[1] && _nodes[2] == key._nodes[2];
243 * Less than operator.
245 * @param key TriangleFaceKey with which to compare
246 * @return true if this object has the three nodes less than the nodes of the key object, false if not
248 bool operator<(const TriangleFaceKey& key) const
250 for (int i = 0; i < 3; ++i)
252 if (_nodes[i] < key._nodes[i])
256 else if (_nodes[i] > key._nodes[i])
265 * Returns a hash value for the object, based on its three nodes.
266 * This value is not unique for each face.
268 * @return a hash value for the object
275 inline static void Sort3Ints(int* sorted, int node1, int node2, int node3);
278 /// global numbers of the three nodes, sorted in ascending order
281 /// hash value for the object, calculated in the constructor
286 * Method to sort three integers in ascending order
288 * @param sorted int[3] array in which to store the result
289 * @param x1 first integer
290 * @param x2 second integer
291 * @param x3 third integer
293 inline void TriangleFaceKey::Sort3Ints(int* sorted, int x1, int x2, int x3)
301 sorted[1] = x2 < x3 ? x2 : x3;
302 sorted[2] = x2 < x3 ? x3 : x2;
318 sorted[1] = x1 < x3 ? x1 : x3;
319 sorted[2] = x1 < x3 ? x3 : x1;
332 * \brief Template specialization of INTERP_KERNEL::hash<T> function object for use with a
333 * with TriangleFaceKey as key class.
336 template<> class hash<INTERP_KERNEL::TriangleFaceKey>
340 * Operator() that returns the precalculated hashvalue of a TriangleFaceKey object.
342 * @param key a TriangleFaceKey object
343 * @return an integer hash value for key
345 int operator()(const INTERP_KERNEL::TriangleFaceKey& key) const
347 return key.hashVal();
352 namespace INTERP_KERNEL
355 * \brief Class calculating the volume of intersection between a tetrahedral target element and
356 * source elements with triangular or quadratilateral faces.
359 template<class MyMeshType>
364 SplitterTetra(const MyMeshType& srcMesh, const double** tetraCorners, const typename MyMeshType::MyConnType *nodesId);
366 SplitterTetra(const MyMeshType& srcMesh, const double tetraCorners[12]);
370 double intersectSourceCell(typename MyMeshType::MyConnType srcCell, double* baryCentre=0);
371 double intersectSourceFace(const NormalizedCellType polyType,
372 const int polyNodesNbr,
373 const int *const polyNodes,
374 const double *const *const polyCoords,
375 const double dimCaracteristic,
376 const double precision,
377 std::multiset<TriangleFaceKey>& listOfTetraFacesTreated,
378 std::set<TriangleFaceKey>& listOfTetraFacesColinear);
380 double intersectTetra(const double** tetraCorners);
382 typename MyMeshType::MyConnType getId(int id) { return _conn[id]; }
384 void splitIntoDualCells(SplitterTetra<MyMeshType> **output);
386 void splitMySelfForDual(double* output, int i, typename MyMeshType::MyConnType& nodeId);
388 void clearVolumesCache();
391 inline void checkIsOutside(const double* pt, bool* isOutside, const double errTol = DEFAULT_ABS_TOL) const;
392 inline void checkIsStrictlyOutside(const double* pt, bool* isStrictlyOutside, const double errTol = DEFAULT_ABS_TOL) const;
393 inline void calculateNode(typename MyMeshType::MyConnType globalNodeNum);
394 inline void calculateNode2(typename MyMeshType::MyConnType globalNodeNum, const double* node);
395 inline void calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key);
396 inline void calculateSurface(TransformedTriangle& tri, const TriangleFaceKey& key);
398 static inline bool IsFacesCoplanar(const double *const planeNormal, const double planeConstant,
399 const double *const *const coordsFace, const double precision);
400 static inline double CalculateIntersectionSurfaceOfCoplanarTriangles(const double *const planeNormal,
401 const double planeConstant,
402 const double *const p1, const double *const p2, const double *const p3,
403 const double *const p4, const double *const p5, const double *const p6,
404 const double dimCaracteristic, const double precision);
407 SplitterTetra(const SplitterTetra& t);
409 /// disallow assignment
410 SplitterTetra& operator=(const SplitterTetra& t);
413 /// affine transform associated with this target element
414 TetraAffineTransform* _t;
416 /// HashMap relating node numbers to transformed nodes, used for caching
417 HashMap< int , double* > _nodes;
419 /// HashMap relating triangular faces to calculated volume contributions, used for caching
420 HashMap< TriangleFaceKey, double > _volumes;
422 /// reference to the source mesh
423 const MyMeshType& _src_mesh;
425 // node id of the first node in target mesh in C mode.
426 typename MyMeshType::MyConnType _conn[4];
430 /// Smallest volume of the intersecting elements in the transformed space that will be returned as non-zero.
431 /// 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.
432 static const double SPARSE_TRUNCATION_LIMIT;
436 * Function used to filter out elements by checking if they belong to one of the halfspaces
437 * x <= 0, x >= 1, y <= 0, y >= 1, z <= 0, z >= 1, (indexed 0 - 7). The function updates an array of boolean variables
438 * which indicates whether the points that have already been checked are all in a halfspace. For each halfspace,
439 * 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.
441 * @param pt double[3] containing the coordiantes of a transformed point
442 * @param isOutside bool[8] which indicate the results of earlier checks.
444 template<class MyMeshType>
445 inline void SplitterTetra<MyMeshType>::checkIsOutside(const double* pt, bool* isOutside, const double errTol) const
447 isOutside[0] = isOutside[0] && (pt[0] < errTol);
448 isOutside[1] = isOutside[1] && (pt[0] > (1.0-errTol) );
449 isOutside[2] = isOutside[2] && (pt[1] < errTol);
450 isOutside[3] = isOutside[3] && (pt[1] > (1.0-errTol));
451 isOutside[4] = isOutside[4] && (pt[2] < errTol);
452 isOutside[5] = isOutside[5] && (pt[2] > (1.0-errTol));
453 isOutside[6] = isOutside[6] && (1.0 - pt[0] - pt[1] - pt[2] < errTol);
454 isOutside[7] = isOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] > (1.0-errTol) );
457 template<class MyMeshType>
458 inline void SplitterTetra<MyMeshType>::checkIsStrictlyOutside(const double* pt, bool* isStrictlyOutside, const double errTol) const
460 isStrictlyOutside[0] = isStrictlyOutside[0] && (pt[0] < -errTol);
461 isStrictlyOutside[1] = isStrictlyOutside[1] && (pt[0] > (1.0 + errTol));
462 isStrictlyOutside[2] = isStrictlyOutside[2] && (pt[1] < -errTol);
463 isStrictlyOutside[3] = isStrictlyOutside[3] && (pt[1] > (1.0 + errTol));
464 isStrictlyOutside[4] = isStrictlyOutside[4] && (pt[2] < -errTol);
465 isStrictlyOutside[5] = isStrictlyOutside[5] && (pt[2] > (1.0 + errTol));
466 isStrictlyOutside[6] = isStrictlyOutside[6] && (1.0 - pt[0] - pt[1] - pt[2] < -errTol);
467 isStrictlyOutside[7] = isStrictlyOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] > (1.0 + errTol));
471 * Calculates the transformed node with a given global node number.
472 * Gets the coordinates for the node in _src_mesh with the given global number and applies TetraAffineTransform
473 * _t to it. Stores the result in the cache _nodes. The non-existance of the node in _nodes should be verified before
476 * @param globalNodeNum global node number of the node in the mesh _src_mesh
479 template<class MyMeshType>
480 inline void SplitterTetra<MyMeshType>::calculateNode(typename MyMeshType::MyConnType globalNodeNum)
482 const double* node = _src_mesh.getCoordinatesPtr()+MyMeshType::MY_SPACEDIM*globalNodeNum;
483 double* transformedNode = new double[MyMeshType::MY_SPACEDIM];
484 assert(transformedNode != 0);
485 _t->apply(transformedNode, node);
486 _nodes[globalNodeNum] = transformedNode;
491 * Calculates the transformed node with a given global node number.
492 * Applies TetraAffineTransform * _t to it.
493 * Stores the result in the cache _nodes. The non-existence of the node in _nodes should be verified before * calling.
494 * The only difference with the previous method calculateNode is that the coordinates of the node are passed in arguments
495 * and are not recalculated in order to optimize the method.
497 * @param globalNodeNum global node number of the node in the mesh _src_mesh
500 template<class MyMeshType>
501 inline void SplitterTetra<MyMeshType>::calculateNode2(typename MyMeshType::MyConnType globalNodeNum, const double* node)
503 double* transformedNode = new double[MyMeshType::MY_SPACEDIM];
504 assert(transformedNode != 0);
505 _t->apply(transformedNode, node);
506 _nodes[globalNodeNum] = transformedNode;
510 * Calculates the volume contribution from the given TransformedTriangle and stores it with the given key in .
511 * Calls TransformedTriangle::calculateIntersectionVolume to perform the calculation.
513 * @param tri triangle for which to calculate the volume contribution
514 * @param key key associated with the face
516 template<class MyMeshType>
517 inline void SplitterTetra<MyMeshType>::calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key)
519 const double vol = tri.calculateIntersectionVolume();
520 _volumes.insert(std::make_pair(key, vol));
524 * Calculates the surface contribution from the given TransformedTriangle and stores it with the given key in.
525 * Calls TransformedTriangle::calculateIntersectionSurface to perform the calculation.
527 * @param tri triangle for which to calculate the surface contribution
528 * @param key key associated with the face
530 template<class MyMeshType>
531 inline void SplitterTetra<MyMeshType>::calculateSurface(TransformedTriangle& tri, const TriangleFaceKey& key)
533 const double surf = tri.calculateIntersectionSurface(_t);
534 _volumes.insert(std::make_pair(key, surf));
537 template<class MyMeshTypeT, class MyMeshTypeS=MyMeshTypeT>
541 SplitterTetra2(const MyMeshTypeT& targetMesh, const MyMeshTypeS& srcMesh, SplittingPolicy policy);
543 void releaseArrays();
544 void splitTargetCell2(typename MyMeshTypeT::MyConnType targetCell, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
545 void splitTargetCell(typename MyMeshTypeT::MyConnType targetCell, typename MyMeshTypeT::MyConnType nbOfNodesT,
546 typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);//to suppress
547 void fiveSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);//to suppress
548 void sixSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);//to suppress
549 void calculateGeneral24Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);//to suppress
550 void calculateGeneral48Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);//to suppress
551 void splitPyram5(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);//to suppress
552 void splitConvex(typename MyMeshTypeT::MyConnType targetCell,//to suppress
553 typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);//to suppress
554 void calculateSubNodes(const MyMeshTypeT& targetMesh, typename MyMeshTypeT::MyConnType targetCell);//to suppress
555 inline const double* getCoordsOfSubNode(typename MyMeshTypeT::MyConnType node);//to suppress
556 inline const double* getCoordsOfSubNode2(typename MyMeshTypeT::MyConnType node, typename MyMeshTypeT::MyConnType& nodeId);//to suppress
558 inline void calcBarycenter(int n, double* barycenter, const typename MyMeshTypeT::MyConnType* pts);//to suppress
560 const MyMeshTypeT& _target_mesh;
561 const MyMeshTypeS& _src_mesh;
562 SplittingPolicy _splitting_pol;
563 /// vector of pointers to double[3] containing the coordinates of the
564 /// (sub) - nodes of split target cell
565 std::vector<const double*> _nodes;
566 std::vector<typename MyMeshTypeT::MyConnType> _node_ids;
570 * Calculates the barycenter of n (sub) - nodes
572 * @param n number of nodes for which to calculate barycenter
573 * @param barycenter pointer to double[3] array in which to store the result
574 * @param pts pointer to int[n] array containing the (sub)-nodes for which to calculate the barycenter
576 template<class MyMeshTypeT, class MyMeshTypeS>
578 inline void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calcBarycenter(int n, double* barycenter, const typename MyMeshTypeT::MyConnType* pts)
580 barycenter[0] = barycenter[1] = barycenter[2] = 0.0;
581 for(int i = 0; i < n ; ++i)
583 const double* pt = getCoordsOfSubNode(pts[i]);
584 barycenter[0] += pt[0];
585 barycenter[1] += pt[1];
586 barycenter[2] += pt[2];
595 * Accessor to the coordinates of a given (sub)-node
597 * @param node local number of the (sub)-node 0,..,7 are the elements nodes, sub-nodes are numbered from 8,..
598 * @return pointer to double[3] containing the coordinates of the nodes
600 template<class MyMeshTypeT, class MyMeshTypeS>
601 inline const double* SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::getCoordsOfSubNode(typename MyMeshTypeT::MyConnType node)
603 // replace "at()" with [] for unsafe but faster access
604 return _nodes.at(node);
608 * Accessor to the coordinates of a given (sub)-node
610 * @param node local number of the (sub)-node 0,..,7 are the elements nodes, sub-nodes are numbered from 8,..
611 * @param nodeId is an output that is node id in target whole mesh in C mode.
612 * @return pointer to double[3] containing the coordinates of the nodes
614 template<class MyMeshTypeT, class MyMeshTypeS>
615 const double* SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::getCoordsOfSubNode2(typename MyMeshTypeT::MyConnType node, typename MyMeshTypeT::MyConnType& nodeId)
617 const double *ret=_nodes.at(node);
619 nodeId=_node_ids[node];