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_TXX__
20 #define __SPLITTERTETRA_TXX__
22 #include "SplitterTetra.hxx"
24 #include "TetraAffineTransform.hxx"
25 #include "TransformedTriangle.hxx"
26 #include "MeshUtils.hxx"
27 #include "VectorUtils.hxx"
28 #include "CellModel.hxx"
36 /// Smallest volume of the intersecting elements in the transformed space that will be returned as non-zero.
37 /// 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.
38 #define SPARSE_TRUNCATION_LIMIT 1.0e-14
40 namespace INTERP_KERNEL
44 * Constructor creating object from target cell global number
46 * @param srcMesh mesh containing the source elements
47 * @param targetMesh mesh containing the target elements
48 * @param targetCell global number of the target cell
51 /*template<class MyMeshType>
52 SplitterTetra<MyMeshType>::SplitterTetra(const MyMeshType& srcMesh, const MyMeshType& targetMesh, typename MyMeshType::MyConnType targetCell)
53 : _src_mesh(srcMesh), _t(0)
55 // get array of corners of target tetraedron
56 const double* tetraCorners[4];
57 for(int i = 0 ; i < 4 ; ++i)
58 tetraCorners[i] = getCoordsOfNode(i, targetCell, targetMesh);
59 // create the affine transform
60 createAffineTransform(tetraCorners);
64 * output is expected to be allocated with 24*sizeof(void*) in order to store the 24 tetras.
65 * These tetras have to be deallocated.
67 template<class MyMeshType>
68 void SplitterTetra<MyMeshType>::splitIntoDualCells(SplitterTetra<MyMeshType> **output)
71 const double *tmp2[4]={tmp,tmp+3,tmp+6,tmp+9};
72 typename MyMeshType::MyConnType conn[4]={-1,-1,-1,-1};
75 splitMySelfForDual(tmp,i,conn[0]);
76 output[i]=new SplitterTetra<MyMeshType>(_src_mesh,tmp2,conn);
81 * Constructor creating object from the four corners of the tetrahedron.
83 * @param srcMesh mesh containing the source elements
84 * @param tetraCorners array of four pointers to double[3] arrays containing the coordinates of the
85 * corners of the tetrahedron
87 template<class MyMeshType>
88 SplitterTetra<MyMeshType>::SplitterTetra(const MyMeshType& srcMesh, const double** tetraCorners, const typename MyMeshType::MyConnType *nodesId)
89 : _t(0), _src_mesh(srcMesh)
91 std::copy(nodesId,nodesId+4,_conn);
92 _coords[0]=tetraCorners[0][0]; _coords[1]=tetraCorners[0][1]; _coords[2]=tetraCorners[0][2];
93 _coords[3]=tetraCorners[1][0]; _coords[4]=tetraCorners[1][1]; _coords[5]=tetraCorners[1][2];
94 _coords[6]=tetraCorners[2][0]; _coords[7]=tetraCorners[2][1]; _coords[8]=tetraCorners[2][2];
95 _coords[9]=tetraCorners[3][0]; _coords[10]=tetraCorners[3][1]; _coords[11]=tetraCorners[3][2];
96 // create the affine transform
97 createAffineTransform(tetraCorners);
101 * This contructor is used to build part of 1/24th dual cell of tetraCorners.
102 * @param i is in 0..23 included.
103 * @param nodeId is the id of first node of this in target mesh in C mode.
105 /*template<class MyMeshType>
106 SplitterTetra<MyMeshType>::SplitterTetra(const MyMeshType& srcMesh, const double** tetraCorners, int i, typename MyMeshType::MyConnType nodeId)
107 : _t(0), _src_mesh(srcMesh), _conn(nodeId)
109 double *newCoords[4];
110 splitMySelfForDual(tetraCorners,newCoords,i);
111 createAffineTransform(newCoords);
117 * Deletes _t and the coordinates (double[3] vectors) in _nodes
120 template<class MyMeshType>
121 SplitterTetra<MyMeshType>::~SplitterTetra()
124 for(hash_map< int, double* >::iterator iter = _nodes.begin(); iter != _nodes.end() ; ++iter)
125 delete[] iter->second;
129 * This method destroys the 4 pointers pointed by tetraCorners[0],tetraCorners[1],tetraCorners[2] and tetraCorners[3]
130 * @param i is in 0..23 included.
131 * @param output is expected to be sized of 12 in order to.
133 template<class MyMeshType>
134 void SplitterTetra<MyMeshType>::splitMySelfForDual(double* output, int i, typename MyMeshType::MyConnType& nodeId)
138 nodeId=_conn[offset];
139 tmp[0]=_coords+3*offset; tmp[1]=_coords+((offset+1)%4)*3; tmp[2]=_coords+((offset+2)%4)*3; tmp[3]=_coords+((offset+3)%4)*3;
141 int case1=caseToTreat/2;
142 int case2=caseToTreat%2;
143 const int tab[3][2]={{1,2},{3,2},{1,3}};
144 const int *curTab=tab[case1];
145 double pt0[3]; pt0[0]=(tmp[curTab[case2]][0]+tmp[0][0])/2.; pt0[1]=(tmp[curTab[case2]][1]+tmp[0][1])/2.; pt0[2]=(tmp[curTab[case2]][2]+tmp[0][2])/2.;
146 double pt1[3]; pt1[0]=(tmp[0][0]+tmp[curTab[0]][0]+tmp[curTab[1]][0])/3.; pt1[1]=(tmp[0][1]+tmp[curTab[0]][1]+tmp[curTab[1]][1])/3.; pt1[1]=(tmp[0][2]+tmp[curTab[0]][2]+tmp[curTab[1]][2])/3.;
147 double pt2[3]; pt2[0]=(tmp[0][0]+tmp[1][0]+tmp[2][0]+tmp[3][0])/4.; pt2[1]=(tmp[0][1]+tmp[1][1]+tmp[2][1]+tmp[3][1])/4.; pt2[2]=(tmp[0][2]+tmp[1][2]+tmp[2][2]+tmp[3][2])/4.;
148 std::copy(pt1,pt1+3,output+case2*3);
149 std::copy(pt0,pt0+3,output+(abs(case2-1))*3);
150 std::copy(pt2,pt2+3,output+2*3);
151 std::copy(tmp[0],tmp[0]+3,output+3*3);
155 * Calculates the volume of intersection of an element in the source mesh and the target element.
156 * It first calculates the transformation that takes the target tetrahedron into the unit tetrahedron. After that, the
157 * faces of the source element are triangulated and the calculated transformation is applied
158 * to each triangle. The algorithm of Grandy, implemented in INTERP_KERNEL::TransformedTriangle is used
159 * to calculate the contribution to the volume from each triangle. The volume returned is the sum of these contributions
160 * divided by the determinant of the transformation.
162 * The class will cache the intermediary calculations of transformed nodes of source cells and volumes associated
163 * with triangulated faces to avoid having to recalculate these.
165 * @param element global number of the source element in C mode.
167 template<class MyMeshType>
168 double SplitterTetra<MyMeshType>::intersectSourceCell(typename MyMeshType::MyConnType element)
170 typedef typename MyMeshType::MyConnType ConnType;
171 const NumberingPolicy numPol=MyMeshType::My_numPol;
172 //{ could be done on outside?
173 // check if we have planar tetra element
174 if(_t->determinant() == 0.0)
177 LOG(2, "Planar tetra -- volume 0");
182 NormalizedCellType normCellType=_src_mesh.getTypeOfElement(OTT<ConnType,numPol>::indFC(element));
183 const CellModel& cellModelCell=CellModel::getCellModel(normCellType);
184 unsigned nbOfNodes4Type=cellModelCell.getNumberOfNodes();
185 // halfspace filtering
186 bool isOutside[8] = {true, true, true, true, true, true, true, true};
187 bool isTargetOutside = false;
189 // calculate the coordinates of the nodes
190 int *cellNodes=new int[nbOfNodes4Type];
191 for(int i = 0;i<nbOfNodes4Type;++i)
193 // we could store mapping local -> global numbers too, but not sure it is worth it
194 const int globalNodeNum = getGlobalNumberOfNode(i, OTT<ConnType,numPol>::indFC(element), _src_mesh);
195 cellNodes[i]=globalNodeNum;
196 if(_nodes.find(globalNodeNum) == _nodes.end())
198 //for(hash_map< int , double* >::iterator iter3=_nodes.begin();iter3!=_nodes.end();iter3++)
199 // std::cout << (*iter3).first << " ";
200 //std::cout << std::endl << "*** " << globalNodeNum << std::endl;
201 calculateNode(globalNodeNum);
204 checkIsOutside(_nodes[globalNodeNum], isOutside);
207 // halfspace filtering check
208 // NB : might not be beneficial for caching of triangles
209 for(int i = 0; i < 8; ++i)
213 isTargetOutside = true;
217 double totalVolume = 0.0;
221 for(unsigned ii = 0 ; ii < cellModelCell.getNumberOfSons() ; ++ii)
223 NormalizedCellType faceType = cellModelCell.getSonType(ii);
224 const CellModel& faceModel=CellModel::getCellModel(faceType);
225 assert(faceModel.getDimension() == 2);
226 int *faceNodes=new int[faceModel.getNumberOfNodes()];
227 cellModelCell.fillSonCellNodalConnectivity(ii,cellNodes,faceNodes);
232 // create the face key
233 TriangleFaceKey key = TriangleFaceKey(faceNodes[0], faceNodes[1], faceNodes[2]);
235 // calculate the triangle if needed
236 if(_volumes.find(key) == _volumes.end())
238 TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[1]], _nodes[faceNodes[2]]);
239 calculateVolume(tri, key);
240 totalVolume += _volumes[key];
242 // count negative as face has reversed orientation
243 totalVolume -= _volumes[key];
250 // simple triangulation of faces along a diagonal :
261 //? not sure if this always works
263 // calculate the triangles if needed
265 // local nodes 1, 2, 3
266 TriangleFaceKey key1 = TriangleFaceKey(faceNodes[0], faceNodes[1], faceNodes[2]);
267 if(_volumes.find(key1) == _volumes.end())
269 TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[1]], _nodes[faceNodes[2]]);
270 calculateVolume(tri, key1);
271 totalVolume += _volumes[key1];
273 // count negative as face has reversed orientation
274 totalVolume -= _volumes[key1];
277 // local nodes 1, 3, 4
278 TriangleFaceKey key2 = TriangleFaceKey(faceNodes[0], faceNodes[2], faceNodes[3]);
279 if(_volumes.find(key2) == _volumes.end())
281 TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[2]], _nodes[faceNodes[3]]);
282 calculateVolume(tri, key2);
283 totalVolume += _volumes[key2];
287 // count negative as face has reversed orientation
288 totalVolume -= _volumes[key2];
294 std::cout << "+++ Error : Only elements with triangular and quadratilateral faces are supported at the moment." << std::endl;
301 // reset if it is very small to keep the matrix sparse
302 // is this a good idea?
303 if(epsilonEqual(totalVolume, 0.0, SPARSE_TRUNCATION_LIMIT))
308 LOG(2, "Volume = " << totalVolume << ", det= " << _t->determinant());
310 // NB : fault in article, Grandy, [8] : it is the determinant of the inverse transformation
311 // that should be used (which is equivalent to dividing by the determinant)
312 return std::fabs(1.0 / _t->determinant() * totalVolume) ;
315 ////////////////////////////////////////////////////////
317 template<class MyMeshType>
318 SplitterTetra2<MyMeshType>::SplitterTetra2(const MyMeshType& targetMesh, const MyMeshType& srcMesh, SplittingPolicy policy):_target_mesh(targetMesh),_src_mesh(srcMesh),
319 _splitting_pol(policy)
323 template<class MyMeshType>
324 SplitterTetra2<MyMeshType>::~SplitterTetra2()
329 template<class MyMeshType>
330 void SplitterTetra2<MyMeshType>::releaseArrays()
332 // free potential sub-mesh nodes that have been allocated
335 std::vector<const double*>::iterator iter = _nodes.begin() + 8;
336 while(iter != _nodes.end())
346 * @param targetCell in C mode.
347 * @param tetra is the output result tetra containers.
349 template<class MyMeshType>
350 void SplitterTetra2<MyMeshType>::splitTargetCell(typename MyMeshType::MyConnType targetCell, typename MyMeshType::MyConnType nbOfNodesT,
351 typename std::vector< SplitterTetra<MyMeshType>* >& tetra)
353 typedef typename MyMeshType::MyConnType ConnType;
354 const NumberingPolicy numPol=MyMeshType::My_numPol;
355 const int numTetra = static_cast<int>(_splitting_pol);
361 const double *nodes[4];
363 for(int node = 0; node < 4 ; ++node)
365 nodes[node]=getCoordsOfNode2(node, OTT<ConnType,numPol>::indFC(targetCell),_target_mesh,conn[node]);
367 std::copy(conn,conn+4,_node_ids.begin());
368 SplitterTetra<MyMeshType>* t = new SplitterTetra<MyMeshType>(_src_mesh, nodes,conn);
373 // pre-calculate nodes
374 calculateSubNodes(_target_mesh, OTT<ConnType,numPol>::indFC(targetCell));
376 tetra.reserve(numTetra);
377 _nodes.reserve(30); // we never have more than this
379 switch(_splitting_pol)
383 const int subZone[8] = { 0, 1, 2, 3, 4, 5, 6, 7 };
384 fiveSplit(subZone,tetra);
390 const int subZone[8] = { 0, 1, 2, 3, 4, 5, 6, 7 };
391 sixSplit(subZone,tetra);
397 calculateGeneral24Tetra(tetra);
403 calculateGeneral48Tetra(tetra);
412 * Splits the hexahedron into five tetrahedra.
413 * This method adds five SplitterTetra objects to the vector tetra.
415 * @param subZone the local node numbers corresponding to the hexahedron corners - these are mapped onto {0,..,7}. Providing this allows the
416 * splitting to be reused on the subzones of the GENERAL_* types of splitting
418 template<class MyMeshType>
419 void SplitterTetra2<MyMeshType>::fiveSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshType>* >& tetra)
421 // Schema according to which the splitting is performed.
422 // Each line represents one tetrahedron. The numbering is as follows :
435 static const int SPLIT_NODES_5[20] =
445 for(int i = 0; i < 5; ++i)
447 const double* nodes[4];
449 for(int j = 0; j < 4; ++j)
451 nodes[j] = getCoordsOfSubNode2(subZone[ SPLIT_NODES_5[4*i+j] ],conn[j]);
453 SplitterTetra<MyMeshType>* t = new SplitterTetra<MyMeshType>(_src_mesh, nodes,conn);
459 * Splits the hexahedron into six tetrahedra.
460 * This method adds six SplitterTetra objects to the vector tetra.
462 * @param subZone the local node numbers corresponding to the hexahedron corners - these are mapped onto {0,..,7}. Providing this allows the
463 * splitting to be reused on the subzones of the GENERAL_* types of splitting
465 template<class MyMeshType>
466 void SplitterTetra2<MyMeshType>::sixSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshType>* >& tetra)
468 // Schema according to which the splitting is performed.
469 // Each line represents one tetrahedron. The numbering is as follows :
481 static const int SPLIT_NODES_6[24] =
491 for(int i = 0; i < 6; ++i)
493 const double* nodes[4];
495 for(int j = 0; j < 4; ++j)
497 nodes[j] = getCoordsOfSubNode2(subZone[ SPLIT_NODES_6[4*i+j] ],conn[j]);
499 SplitterTetra<MyMeshType>* t = new SplitterTetra<MyMeshType>(_src_mesh, nodes,conn);
505 * Splits the hexahedron into 24 tetrahedra.
506 * The splitting is done by combining the barycenter of the tetrahedron, the barycenter of each face
507 * and the nodes of each edge of the face. This creates 6 faces * 4 edges / face = 24 tetrahedra.
508 * The submesh nodes introduced are the barycenters of the faces and the barycenter of the cell.
511 template<class MyMeshType>
512 void SplitterTetra2<MyMeshType>::calculateGeneral24Tetra(typename std::vector< SplitterTetra<MyMeshType>* >& tetra)
514 // The two nodes of the original mesh cell used in each tetrahedron.
515 // The tetrahedra all have nodes (cellCenter, faceCenter, edgeNode1, edgeNode2)
516 // For the correspondance of the nodes, see the GENERAL_48_SUB_NODES table in calculateSubNodes
517 static const int TETRA_EDGES[48] =
519 // face with center 9
524 // face with center 10
529 // face with center 11
534 // face with center 12
539 // face with center 13
544 // face with center 14
551 // nodes to use for tetrahedron
552 const double* nodes[4];
554 // get the cell center
555 nodes[0] = getCoordsOfSubNode2(14,conn[0]);
557 for(int faceCenterNode = 8; faceCenterNode < 14; ++faceCenterNode)
559 // get the face center
560 nodes[1] = getCoordsOfSubNode2(faceCenterNode,conn[1]);
561 for(int j = 0; j < 4; ++j)
563 const int row = 4*(faceCenterNode - 9) + j;
564 nodes[2] = getCoordsOfSubNode2(TETRA_EDGES[2*row],conn[2]);
565 nodes[3] = getCoordsOfSubNode2(TETRA_EDGES[2*row + 1],conn[3]);
567 SplitterTetra<MyMeshType>* t = new SplitterTetra<MyMeshType>(_src_mesh, nodes, conn);
575 * Splits the hexahedron into 48 tetrahedra.
576 * The splitting is done by introducing the midpoints of all the edges
577 * and the barycenter of the element as submesh nodes. The 8 hexahedral subzones thus defined
578 * are then split into 6 tetrahedra each, as in Grandy, p. 449. The division of the subzones
579 * is done by calling sixSplit().
582 template<class MyMeshType>
583 void SplitterTetra2<MyMeshType>::calculateGeneral48Tetra(typename std::vector< SplitterTetra<MyMeshType>* >& tetra)
585 // Define 8 hexahedral subzones as in Grandy, p449
586 // the values correspond to the nodes that correspond to nodes 1,2,3,4,5,6,7,8 in the subcell
587 // For the correspondance of the nodes, see the GENERAL_48_SUB_NODES table in calculateSubNodes
588 static const int subZones[64] =
590 0,8,21,12,9,20,26,22,
591 8,1,13,21,20,10,23,26,
592 12,21,16,3,22,26,25,17,
593 21,13,2,16,26,23,18,25,
594 9,20,26,22,4,11,24,14,
595 20,10,23,26,11,5,15,24,
596 22,26,25,17,14,24,19,7,
597 26,23,18,25,24,15,6,19
600 for(int i = 0; i < 8; ++i)
602 sixSplit(&subZones[8*i],tetra);
607 * Precalculates all the nodes.
608 * Retrieves the mesh nodes and allocates the necessary sub-mesh
609 * nodes according to the splitting policy used.
610 * This method is meant to be called once by the constructor.
612 * @param targetMesh the target mesh
613 * @param targetCell the global number of the cell that the object represents, in targetMesh mode.
614 * @param policy the splitting policy of the object
617 template<class MyMeshType>
618 void SplitterTetra2<MyMeshType>::calculateSubNodes(const MyMeshType& targetMesh, typename MyMeshType::MyConnType targetCell)
620 // retrieve real mesh nodes
622 for(int node = 0; node < 8 ; ++node)
624 // calculate only normal nodes
625 _nodes.push_back(getCoordsOfNode2(node, targetCell, targetMesh,_node_ids[node]));
628 // create sub-mesh nodes if needed
629 switch(_splitting_pol)
633 // Each sub-node is the barycenter of 4 other nodes.
634 // For the faces, these are on the orignal mesh.
635 // For the barycenter, the four face sub-nodes are used.
636 static const int GENERAL_24_SUB_NODES[28] =
638 0,1,4,5,// sub-node 9 (face)
639 0,1,2,3,// sub-node 10 (face)
640 0,3,4,7,// sub-node 11 (face)
641 1,2,5,6,// sub-node 12 (face)
642 4,5,6,7,// sub-node 13 (face)
643 2,3,6,7,// sub-node 14 (face)
644 9,10,11,12// sub-node 15 (cell)
647 for(int i = 0; i < 7; ++i)
649 double* barycenter = new double[3];
650 calcBarycenter<4>(barycenter, &GENERAL_24_SUB_NODES[4*i]);
651 _nodes.push_back(barycenter);
658 // Each sub-node is the barycenter of two other nodes.
659 // For the edges, these lie on the original mesh.
660 // For the faces, these are the edge sub-nodes.
661 // For the cell these are two face sub-nodes.
662 static const int GENERAL_48_SUB_NODES[38] =
664 0,1, // sub-node 9 (edge)
665 0,4, // sub-node 10 (edge)
666 1,5, // sub-node 11 (edge)
667 4,5, // sub-node 12 (edge)
668 0,3, // sub-node 13 (edge)
669 1,2, // sub-node 14 (edge)
670 4,7, // sub-node 15 (edge)
671 5,6, // sub-node 16 (edge)
672 2,3, // sub-node 17 (edge)
673 3,7, // sub-node 18 (edge)
674 2,6, // sub-node 19 (edge)
675 6,7, // sub-node 20 (edge)
676 8,11, // sub-node 21 (face)
677 12,13, // sub-node 22 (face)
678 9,17, // sub-node 23 (face)
679 10,18, // sub-node 24 (face)
680 14,15, // sub-node 25 (face)
681 16,19, // sub-node 26 (face)
682 20,25 // sub-node 27 (cell)
685 for(int i = 0; i < 19; ++i)
687 double* barycenter = new double[3];
688 calcBarycenter<2>(barycenter, &GENERAL_48_SUB_NODES[2*i]);
689 _nodes.push_back(barycenter);