1 // Copyright (C) 2007-2012 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];
366 * Creates the affine transform _t from the corners of the tetrahedron. Used by the constructors
368 * @param corners double*[4] array containing pointers to four double[3] arrays with the
369 * coordinates of the corners of the tetrahedron
371 template<class MyMeshType>
372 inline void SplitterTetra<MyMeshType>::createAffineTransform(const double** corners)
374 // create AffineTransform from tetrahedron
375 _t = new TetraAffineTransform( corners );
379 * Function used to filter out elements by checking if they belong to one of the halfspaces
380 * x <= 0, x >= 1, y <= 0, y >= 1, z <= 0, z >= 1, (indexed 0 - 7). The function updates an array of boolean variables
381 * which indicates whether the points that have already been checked are all in a halfspace. For each halfspace,
382 * 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.
384 * @param pt double[3] containing the coordiantes of a transformed point
385 * @param isOutside bool[8] which indicate the results of earlier checks.
387 template<class MyMeshType>
388 inline void SplitterTetra<MyMeshType>::checkIsOutside(const double* pt, bool* isOutside, const double errTol) const
390 isOutside[0] = isOutside[0] && (pt[0] < errTol);
391 isOutside[1] = isOutside[1] && (pt[0] > (1.0-errTol) );
392 isOutside[2] = isOutside[2] && (pt[1] < errTol);
393 isOutside[3] = isOutside[3] && (pt[1] > (1.0-errTol));
394 isOutside[4] = isOutside[4] && (pt[2] < errTol);
395 isOutside[5] = isOutside[5] && (pt[2] > (1.0-errTol));
396 isOutside[6] = isOutside[6] && (1.0 - pt[0] - pt[1] - pt[2] < errTol);
397 isOutside[7] = isOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] > (1.0-errTol) );
400 template<class MyMeshType>
401 inline void SplitterTetra<MyMeshType>::checkIsStrictlyOutside(const double* pt, bool* isStrictlyOutside, const double errTol) const
403 isStrictlyOutside[0] = isStrictlyOutside[0] && (pt[0] < -errTol);
404 isStrictlyOutside[1] = isStrictlyOutside[1] && (pt[0] > (1.0 + errTol));
405 isStrictlyOutside[2] = isStrictlyOutside[2] && (pt[1] < -errTol);
406 isStrictlyOutside[3] = isStrictlyOutside[3] && (pt[1] > (1.0 + errTol));
407 isStrictlyOutside[4] = isStrictlyOutside[4] && (pt[2] < -errTol);
408 isStrictlyOutside[5] = isStrictlyOutside[5] && (pt[2] > (1.0 + errTol));
409 isStrictlyOutside[6] = isStrictlyOutside[6] && (1.0 - pt[0] - pt[1] - pt[2] < -errTol);
410 isStrictlyOutside[7] = isStrictlyOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] > (1.0 + errTol));
414 * Calculates the transformed node with a given global node number.
415 * Gets the coordinates for the node in _src_mesh with the given global number and applies TetraAffineTransform
416 * _t to it. Stores the result in the cache _nodes. The non-existance of the node in _nodes should be verified before
419 * @param globalNodeNum global node number of the node in the mesh _src_mesh
422 template<class MyMeshType>
423 inline void SplitterTetra<MyMeshType>::calculateNode(typename MyMeshType::MyConnType globalNodeNum)
425 const double* node = _src_mesh.getCoordinatesPtr()+MyMeshType::MY_SPACEDIM*globalNodeNum;
426 double* transformedNode = new double[MyMeshType::MY_SPACEDIM];
427 assert(transformedNode != 0);
428 _t->apply(transformedNode, node);
429 _nodes[globalNodeNum] = transformedNode;
434 * Calculates the transformed node with a given global node number.
435 * Applies TetraAffineTransform * _t to it.
436 * Stores the result in the cache _nodes. The non-existence of the node in _nodes should be verified before * calling.
437 * The only difference with the previous method calculateNode is that the coordinates of the node are passed in arguments
438 * and are not recalculated in order to optimize the method.
440 * @param globalNodeNum global node number of the node in the mesh _src_mesh
443 template<class MyMeshType>
444 inline void SplitterTetra<MyMeshType>::calculateNode2(typename MyMeshType::MyConnType globalNodeNum, const double* node)
446 double* transformedNode = new double[MyMeshType::MY_SPACEDIM];
447 assert(transformedNode != 0);
448 _t->apply(transformedNode, node);
449 _nodes[globalNodeNum] = transformedNode;
453 * Calculates the volume contribution from the given TransformedTriangle and stores it with the given key in .
454 * Calls TransformedTriangle::calculateIntersectionVolume to perform the calculation.
456 * @param tri triangle for which to calculate the volume contribution
457 * @param key key associated with the face
459 template<class MyMeshType>
460 inline void SplitterTetra<MyMeshType>::calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key)
462 const double vol = tri.calculateIntersectionVolume();
463 _volumes.insert(std::make_pair(key, vol));
467 * Calculates the surface contribution from the given TransformedTriangle and stores it with the given key in.
468 * Calls TransformedTriangle::calculateIntersectionSurface to perform the calculation.
470 * @param tri triangle for which to calculate the surface contribution
471 * @param key key associated with the face
473 template<class MyMeshType>
474 inline void SplitterTetra<MyMeshType>::calculateSurface(TransformedTriangle& tri, const TriangleFaceKey& key)
476 const double surf = tri.calculateIntersectionSurface(_t);
477 _volumes.insert(std::make_pair(key, surf));
480 template<class MyMeshTypeT, class MyMeshTypeS=MyMeshTypeT>
484 SplitterTetra2(const MyMeshTypeT& targetMesh, const MyMeshTypeS& srcMesh, SplittingPolicy policy);
486 void releaseArrays();
487 void splitTargetCell(typename MyMeshTypeT::MyConnType targetCell, typename MyMeshTypeT::MyConnType nbOfNodesT,
488 typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
489 void fiveSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
490 void sixSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
491 void calculateGeneral24Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
492 void calculateGeneral48Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
493 void splitPyram5(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
494 void splitConvex(typename MyMeshTypeT::MyConnType targetCell,
495 typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
496 void calculateSubNodes(const MyMeshTypeT& targetMesh, typename MyMeshTypeT::MyConnType targetCell);
497 inline const double* getCoordsOfSubNode(typename MyMeshTypeT::MyConnType node);
498 inline const double* getCoordsOfSubNode2(typename MyMeshTypeT::MyConnType node, typename MyMeshTypeT::MyConnType& nodeId);
500 inline void calcBarycenter(int n, double* barycenter, const typename MyMeshTypeT::MyConnType* pts);
502 const MyMeshTypeT& _target_mesh;
503 const MyMeshTypeS& _src_mesh;
504 SplittingPolicy _splitting_pol;
505 /// vector of pointers to double[3] containing the coordinates of the
506 /// (sub) - nodes of split target cell
507 std::vector<const double*> _nodes;
508 std::vector<typename MyMeshTypeT::MyConnType> _node_ids;
512 * Calculates the barycenter of n (sub) - nodes
514 * @param n number of nodes for which to calculate barycenter
515 * @param barycenter pointer to double[3] array in which to store the result
516 * @param pts pointer to int[n] array containing the (sub)-nodes for which to calculate the barycenter
518 template<class MyMeshTypeT, class MyMeshTypeS>
520 inline void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calcBarycenter(int n, double* barycenter, const typename MyMeshTypeT::MyConnType* pts)
522 barycenter[0] = barycenter[1] = barycenter[2] = 0.0;
523 for(int i = 0; i < n ; ++i)
525 const double* pt = getCoordsOfSubNode(pts[i]);
526 barycenter[0] += pt[0];
527 barycenter[1] += pt[1];
528 barycenter[2] += pt[2];
537 * Accessor to the coordinates of a given (sub)-node
539 * @param node local number of the (sub)-node 0,..,7 are the elements nodes, sub-nodes are numbered from 8,..
540 * @return pointer to double[3] containing the coordinates of the nodes
542 template<class MyMeshTypeT, class MyMeshTypeS>
543 inline const double* SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::getCoordsOfSubNode(typename MyMeshTypeT::MyConnType node)
545 // replace "at()" with [] for unsafe but faster access
546 return _nodes.at(node);
550 * Accessor to the coordinates of a given (sub)-node
552 * @param node local number of the (sub)-node 0,..,7 are the elements nodes, sub-nodes are numbered from 8,..
553 * @param nodeId is an output that is node id in target whole mesh in C mode.
554 * @return pointer to double[3] containing the coordinates of the nodes
556 template<class MyMeshTypeT, class MyMeshTypeS>
557 const double* SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::getCoordsOfSubNode2(typename MyMeshTypeT::MyConnType node, typename MyMeshTypeT::MyConnType& nodeId)
559 const double *ret=_nodes.at(node);
561 nodeId=_node_ids[node];