1 // Copyright (C) 2007-2008 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
19 #ifndef __SPLITTERTETRA_HXX__
20 #define __SPLITTERTETRA_HXX__
22 #include "TransformedTriangle.hxx"
23 #include "TetraAffineTransform.hxx"
32 # include <ext/hash_map>
36 using __gnu_cxx::hash_map;
38 using stdext::hash_map;
39 using stdext::hash_compare;
42 namespace INTERP_KERNEL
45 * \brief Class representing a triangular face, used as key in caching hash map in SplitterTetra.
54 * Sorts the given nodes (so that the order in which they are passed does not matter) and
55 * calculates a hash value for the key.
57 * @param node1 global number of the first node of the face
58 * @param node2 global number of the second node of the face
59 * @param node3 global number of the third node of the face
61 TriangleFaceKey(int node1, int node2, int node3)
63 sort3Ints(_nodes, node1, node2, node3);
64 _hashVal = ( _nodes[0] + _nodes[1] + _nodes[2] ) % 29;
68 * Equality comparison operator.
69 * Compares this TriangleFaceKey object to another and determines if they represent the same face.
71 * @param key TriangleFaceKey with which to compare
72 * @return true if key has the same three nodes as this object, false if not
74 bool operator==(const TriangleFaceKey& key) const
76 return _nodes[0] == key._nodes[0] && _nodes[1] == key._nodes[1] && _nodes[2] == key._nodes[2];
80 * Returns a hash value for the object, based on its three nodes.
81 * This value is not unique for each face.
83 * @return a hash value for the object
91 operator size_t () const
97 inline void sort3Ints(int* sorted, int node1, int node2, int node3);
100 /// global numbers of the three nodes, sorted in ascending order
103 /// hash value for the object, calculated in the constructor
108 * Method to sort three integers in ascending order
110 * @param sorted int[3] array in which to store the result
111 * @param x1 first integer
112 * @param x2 second integer
113 * @param x3 third integer
115 inline void TriangleFaceKey::sort3Ints(int* sorted, int x1, int x2, int x3)
123 sorted[1] = x2 < x3 ? x2 : x3;
124 sorted[2] = x2 < x3 ? x3 : x2;
140 sorted[1] = x1 < x3 ? x1 : x3;
141 sorted[2] = x1 < x3 ? x3 : x1;
159 * \brief Template specialization of __gnu_cxx::hash<T> function object for use with a __gnu_cxx::hash_map
160 * with TriangleFaceKey as key class.
164 class hash<INTERP_KERNEL::TriangleFaceKey>
169 * Operator() that returns the precalculated hashvalue of a TriangleFaceKey object.
171 * @param key a TriangleFaceKey object
172 * @return an integer hash value for key
174 int operator()(const INTERP_KERNEL::TriangleFaceKey& key) const
176 return key.hashVal();
181 struct TriangleFaceKeyComparator
183 bool operator()(const INTERP_KERNEL::TriangleFaceKey& key1,
184 const INTERP_KERNEL::TriangleFaceKey& key2 ) const
186 return key1.hashVal() < key2.hashVal();
191 namespace INTERP_KERNEL
195 * \brief Class calculating the volume of intersection between a tetrahedral target element and
196 * source elements with triangular or quadratilateral faces.
199 template<class MyMeshType>
204 SplitterTetra(const MyMeshType& srcMesh, const double** tetraCorners, const typename MyMeshType::MyConnType *nodesId);
208 double intersectSourceCell(typename MyMeshType::MyConnType srcCell);
210 typename MyMeshType::MyConnType getId(int id) { return _conn[id]; }
212 void splitIntoDualCells(SplitterTetra<MyMeshType> **output);
214 void splitMySelfForDual(double* output, int i, typename MyMeshType::MyConnType& nodeId);
218 inline void createAffineTransform(const double** corners);
219 inline void checkIsOutside(const double* pt, bool* isOutside) const;
220 inline void calculateNode(typename MyMeshType::MyConnType globalNodeNum);
221 inline void calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key);
225 SplitterTetra(const SplitterTetra& t);
227 /// disallow assignment
228 SplitterTetra& operator=(const SplitterTetra& t);
231 /// affine transform associated with this target element
232 TetraAffineTransform* _t;
234 /// hash_map relating node numbers to transformed nodes, used for caching
235 hash_map< int , double* > _nodes;
237 /// hash_map relating triangular faces to calculated volume contributions, used for caching
238 hash_map< TriangleFaceKey, double
240 , hash_compare<TriangleFaceKey,TriangleFaceKeyComparator>
244 /// reference to the source mesh
245 const MyMeshType& _src_mesh;
247 // node id of the first node in target mesh in C mode.
248 typename MyMeshType::MyConnType _conn[4];
254 * Creates the affine transform _t from the corners of the tetrahedron. Used by the constructors
256 * @param corners double*[4] array containing pointers to four double[3] arrays with the
257 * coordinates of the corners of the tetrahedron
259 template<class MyMeshType>
260 inline void SplitterTetra<MyMeshType>::createAffineTransform(const double** corners)
262 // create AffineTransform from tetrahedron
263 _t = new TetraAffineTransform( corners );
267 * Function used to filter out elements by checking if they belong to one of the halfspaces
268 * x <= 0, x >= 1, y <= 0, y >= 1, z <= 0, z >= 1, (indexed 0 - 7). The function updates an array of boolean variables
269 * which indicates whether the points that have already been checked are all in a halfspace. For each halfspace,
270 * 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.
272 * @param pt double[3] containing the coordiantes of a transformed point
273 * @param isOutside bool[8] which indicate the results of earlier checks.
275 template<class MyMeshType>
276 inline void SplitterTetra<MyMeshType>::checkIsOutside(const double* pt, bool* isOutside) const
278 isOutside[0] = isOutside[0] && (pt[0] <= 0.0);
279 isOutside[1] = isOutside[1] && (pt[0] >= 1.0);
280 isOutside[2] = isOutside[2] && (pt[1] <= 0.0);
281 isOutside[3] = isOutside[3] && (pt[1] >= 1.0);
282 isOutside[4] = isOutside[4] && (pt[2] <= 0.0);
283 isOutside[5] = isOutside[5] && (pt[2] >= 1.0);
284 isOutside[6] = isOutside[6] && (1.0 - pt[0] - pt[1] - pt[2] <= 0.0);
285 isOutside[7] = isOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] >= 1.0);
289 * Calculates the transformed node with a given global node number.
290 * Gets the coordinates for the node in _src_mesh with the given global number and applies TetraAffineTransform
291 * _t to it. Stores the result in the cache _nodes. The non-existance of the node in _nodes should be verified before
294 * @param globalNodeNum global node number of the node in the mesh _src_mesh
297 template<class MyMeshType>
298 inline void SplitterTetra<MyMeshType>::calculateNode(typename MyMeshType::MyConnType globalNodeNum)
300 const double* node = _src_mesh.getCoordinatesPtr()+MyMeshType::MY_SPACEDIM*globalNodeNum;
301 double* transformedNode = new double[MyMeshType::MY_SPACEDIM];
302 assert(transformedNode != 0);
303 _t->apply(transformedNode, node);
304 _nodes[globalNodeNum] = transformedNode;
308 * Calculates the volume contribution from the given TransformedTriangle and stores it with the given key in .
309 * Calls TransformedTriangle::calculateIntersectionVolume to perform the calculation.
311 * @param tri triangle for which to calculate the volume contribution
312 * @param key key associated with the face
314 template<class MyMeshType>
315 inline void SplitterTetra<MyMeshType>::calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key)
317 const double vol = tri.calculateIntersectionVolume();
318 _volumes.insert(std::make_pair(key, vol));
321 template<class MyMeshType>
325 SplitterTetra2(const MyMeshType& targetMesh, const MyMeshType& srcMesh, SplittingPolicy policy);
327 void releaseArrays();
328 void splitTargetCell(typename MyMeshType::MyConnType targetCell, typename MyMeshType::MyConnType nbOfNodesT,
329 typename std::vector< SplitterTetra<MyMeshType>* >& tetra);
330 void fiveSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshType>* >& tetra);
331 void sixSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshType>* >& tetra);
332 void calculateGeneral24Tetra(typename std::vector< SplitterTetra<MyMeshType>* >& tetra);
333 void calculateGeneral48Tetra(typename std::vector< SplitterTetra<MyMeshType>* >& tetra);
334 void calculateSubNodes(const MyMeshType& targetMesh, typename MyMeshType::MyConnType targetCell);
335 inline const double* getCoordsOfSubNode(typename MyMeshType::MyConnType node);
336 inline const double* getCoordsOfSubNode2(typename MyMeshType::MyConnType node, typename MyMeshType::MyConnType& nodeId);
338 inline void calcBarycenter(double* barycenter, const typename MyMeshType::MyConnType* pts);
340 const MyMeshType& _target_mesh;
341 const MyMeshType& _src_mesh;
342 SplittingPolicy _splitting_pol;
343 /// vector of pointers to double[3] containing the coordinates of the
344 /// (sub) - nodes of split target cell
345 std::vector<const double*> _nodes;
346 std::vector<typename MyMeshType::MyConnType> _node_ids;
350 * Calculates the barycenter of n (sub) - nodes
352 * @param n number of nodes for which to calculate barycenter
353 * @param barycenter pointer to double[3] array in which to store the result
354 * @param pts pointer to int[n] array containing the (sub)-nodes for which to calculate the barycenter
356 template<class MyMeshType>
358 inline void SplitterTetra2<MyMeshType>::calcBarycenter(double* barycenter, const typename MyMeshType::MyConnType* pts)
360 barycenter[0] = barycenter[1] = barycenter[2] = 0.0;
361 for(int i = 0; i < n ; ++i)
363 const double* pt = getCoordsOfSubNode(pts[i]);
364 barycenter[0] += pt[0];
365 barycenter[1] += pt[1];
366 barycenter[2] += pt[2];
375 * Accessor to the coordinates of a given (sub)-node
377 * @param node local number of the (sub)-node 0,..,7 are the elements nodes, sub-nodes are numbered from 8,..
378 * @return pointer to double[3] containing the coordinates of the nodes
380 template<class MyMeshType>
381 inline const double* SplitterTetra2<MyMeshType>::getCoordsOfSubNode(typename MyMeshType::MyConnType node)
383 // replace "at()" with [] for unsafe but faster access
384 return _nodes.at(node);
388 * Accessor to the coordinates of a given (sub)-node
390 * @param node local number of the (sub)-node 0,..,7 are the elements nodes, sub-nodes are numbered from 8,..
391 * @param nodeId is an output that is node id in target whole mesh in C mode.
392 * @return pointer to double[3] containing the coordinates of the nodes
394 template<class MyMeshType>
395 const double* SplitterTetra2<MyMeshType>::getCoordsOfSubNode2(typename MyMeshType::MyConnType node, typename MyMeshType::MyConnType& nodeId)
397 const double *ret=_nodes.at(node);
399 nodeId=_node_ids[node];