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
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"
30 #include "UnitTetraIntersectionBary.hxx"
31 #include "VolSurfFormulae.hxx"
39 /// Smallest volume of the intersecting elements in the transformed space that will be returned as non-zero.
40 /// 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.
41 #define SPARSE_TRUNCATION_LIMIT 1.0e-14
43 namespace INTERP_KERNEL
47 * Constructor creating object from target cell global number
49 * @param srcMesh mesh containing the source elements
50 * @param targetMesh mesh containing the target elements
51 * @param targetCell global number of the target cell
54 /*template<class MyMeshType>
55 SplitterTetra<MyMeshType>::SplitterTetra(const MyMeshType& srcMesh, const MyMeshType& targetMesh, typename MyMeshType::MyConnType targetCell)
56 : _src_mesh(srcMesh), _t(0)
58 // get array of corners of target tetraedron
59 const double* tetraCorners[4];
60 for(int i = 0 ; i < 4 ; ++i)
61 tetraCorners[i] = getCoordsOfNode(i, targetCell, targetMesh);
62 // create the affine transform
63 createAffineTransform(tetraCorners);
67 * output is expected to be allocated with 24*sizeof(void*) in order to store the 24 tetras.
68 * These tetras have to be deallocated.
70 template<class MyMeshType>
71 void SplitterTetra<MyMeshType>::splitIntoDualCells(SplitterTetra<MyMeshType> **output)
74 const double *tmp2[4]={tmp,tmp+3,tmp+6,tmp+9};
75 typename MyMeshType::MyConnType conn[4]={-1,-1,-1,-1};
78 splitMySelfForDual(tmp,i,conn[0]);
79 output[i]=new SplitterTetra<MyMeshType>(_src_mesh,tmp2,conn);
84 * Constructor creating object from the four corners of the tetrahedron.
86 * @param srcMesh mesh containing the source elements
87 * @param tetraCorners array of four pointers to double[3] arrays containing the coordinates of the
88 * corners of the tetrahedron
90 template<class MyMeshType>
91 SplitterTetra<MyMeshType>::SplitterTetra(const MyMeshType& srcMesh, const double** tetraCorners, const typename MyMeshType::MyConnType *nodesId)
92 : _t(0), _src_mesh(srcMesh)
94 std::copy(nodesId,nodesId+4,_conn);
95 _coords[0]=tetraCorners[0][0]; _coords[1]=tetraCorners[0][1]; _coords[2]=tetraCorners[0][2];
96 _coords[3]=tetraCorners[1][0]; _coords[4]=tetraCorners[1][1]; _coords[5]=tetraCorners[1][2];
97 _coords[6]=tetraCorners[2][0]; _coords[7]=tetraCorners[2][1]; _coords[8]=tetraCorners[2][2];
98 _coords[9]=tetraCorners[3][0]; _coords[10]=tetraCorners[3][1]; _coords[11]=tetraCorners[3][2];
99 // create the affine transform
100 createAffineTransform(tetraCorners);
104 * This contructor is used to build part of 1/24th dual cell of tetraCorners.
105 * @param i is in 0..23 included.
106 * @param nodeId is the id of first node of this in target mesh in C mode.
108 /*template<class MyMeshType>
109 SplitterTetra<MyMeshType>::SplitterTetra(const MyMeshType& srcMesh, const double** tetraCorners, int i, typename MyMeshType::MyConnType nodeId)
110 : _t(0), _src_mesh(srcMesh), _conn(nodeId)
112 double *newCoords[4];
113 splitMySelfForDual(tetraCorners,newCoords,i);
114 createAffineTransform(newCoords);
120 * Deletes _t and the coordinates (double[3] vectors) in _nodes
123 template<class MyMeshType>
124 SplitterTetra<MyMeshType>::~SplitterTetra()
127 for(HashMap< int, double* >::iterator iter = _nodes.begin(); iter != _nodes.end() ; ++iter)
128 delete[] iter->second;
132 * \Forget already calculated triangles, which is crucial for calculation of barycenter of intersection
134 template<class MyMeshType>
135 void SplitterTetra<MyMeshType>::clearVolumesCache()
141 * This method destroys the 4 pointers pointed by tetraCorners[0],tetraCorners[1],tetraCorners[2] and tetraCorners[3]
142 * @param i is in 0..23 included.
143 * @param output is expected to be sized of 12 in order to.
145 template<class MyMeshType>
146 void SplitterTetra<MyMeshType>::splitMySelfForDual(double* output, int i, typename MyMeshType::MyConnType& nodeId)
150 nodeId=_conn[offset];
151 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;
153 int case1=caseToTreat/2;
154 int case2=caseToTreat%2;
155 const int tab[3][2]={{1,2},{3,2},{1,3}};
156 const int *curTab=tab[case1];
157 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.;
158 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[2]=(tmp[0][2]+tmp[curTab[0]][2]+tmp[curTab[1]][2])/3.;
159 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.;
160 std::copy(pt1,pt1+3,output+case2*3);
161 std::copy(pt0,pt0+3,output+(abs(case2-1))*3);
162 std::copy(pt2,pt2+3,output+2*3);
163 std::copy(tmp[0],tmp[0]+3,output+3*3);
167 * Calculates the volume of intersection of an element in the source mesh and the target element.
168 * It first calculates the transformation that takes the target tetrahedron into the unit tetrahedron. After that, the
169 * faces of the source element are triangulated and the calculated transformation is applied
170 * to each triangle. The algorithm of Grandy, implemented in INTERP_KERNEL::TransformedTriangle is used
171 * to calculate the contribution to the volume from each triangle. The volume returned is the sum of these contributions
172 * divided by the determinant of the transformation.
174 * The class will cache the intermediary calculations of transformed nodes of source cells and volumes associated
175 * with triangulated faces to avoid having to recalculate these.
177 * @param element global number of the source element in C mode.
179 template<class MyMeshType>
180 double SplitterTetra<MyMeshType>::intersectSourceCell(typename MyMeshType::MyConnType element,
183 typedef typename MyMeshType::MyConnType ConnType;
184 const NumberingPolicy numPol=MyMeshType::My_numPol;
185 //{ could be done on outside?
186 // check if we have planar tetra element
187 if(_t->determinant() == 0.0)
190 LOG(2, "Planar tetra -- volume 0");
195 NormalizedCellType normCellType=_src_mesh.getTypeOfElement(OTT<ConnType,numPol>::indFC(element));
196 const CellModel& cellModelCell=CellModel::GetCellModel(normCellType);
197 unsigned nbOfNodes4Type=cellModelCell.isDynamic() ? _src_mesh.getNumberOfNodesOfElement(OTT<ConnType,numPol>::indFC(element)) : cellModelCell.getNumberOfNodes();
198 // halfspace filtering
199 bool isOutside[8] = {true, true, true, true, true, true, true, true};
200 bool isTargetOutside = false;
202 // calculate the coordinates of the nodes
203 int *cellNodes=new int[nbOfNodes4Type];
204 for(int i = 0;i<(int)nbOfNodes4Type;++i)
206 // we could store mapping local -> global numbers too, but not sure it is worth it
207 const int globalNodeNum = getGlobalNumberOfNode(i, OTT<ConnType,numPol>::indFC(element), _src_mesh);
208 cellNodes[i]=globalNodeNum;
209 if(_nodes.find(globalNodeNum) == _nodes.end())
211 //for(HashMap< int , double* >::iterator iter3=_nodes.begin();iter3!=_nodes.end();iter3++)
212 // std::cout << (*iter3).first << " ";
213 //std::cout << std::endl << "*** " << globalNodeNum << std::endl;
214 calculateNode(globalNodeNum);
217 checkIsOutside(_nodes[globalNodeNum], isOutside);
220 // halfspace filtering check
221 // NB : might not be beneficial for caching of triangles
222 for(int i = 0; i < 8; ++i)
226 isTargetOutside = true;
230 double totalVolume = 0.0;
234 /// calculator of intersection barycentre
235 UnitTetraIntersectionBary baryCalculator( _t->determinant() < 0.);
237 // get nb of sons of a cell
238 const ConnType* rawCellConn = _src_mesh.getConnectivityPtr() + OTT<ConnType,numPol>::conn2C( _src_mesh.getConnectivityIndexPtr()[ element ]);
239 const int rawNbCellNodes = _src_mesh.getConnectivityIndexPtr()[ element+1 ] - _src_mesh.getConnectivityIndexPtr()[ element ];
240 unsigned nbOfSons = cellModelCell.getNumberOfSons2(rawCellConn, rawNbCellNodes);
242 for(unsigned ii = 0 ; ii < nbOfSons; ++ii)
244 // get sons connectivity
245 NormalizedCellType faceType;
246 int *faceNodes, nbFaceNodes=-1;
247 if ( cellModelCell.isDynamic() )
249 faceNodes=new int[nbOfNodes4Type];
250 nbFaceNodes = cellModelCell.fillSonCellNodalConnectivity2(ii,rawCellConn,rawNbCellNodes,faceNodes,faceType);
251 for ( int i = 0; i < nbFaceNodes; ++i )
252 faceNodes[i] = OTT<ConnType,numPol>::coo2C(faceNodes[i]);
256 faceType = cellModelCell.getSonType(ii);
257 const CellModel& faceModel=CellModel::GetCellModel(faceType);
258 assert(faceModel.getDimension() == 2);
259 faceNodes=new int[faceModel.getNumberOfNodes()];
260 cellModelCell.fillSonCellNodalConnectivity(ii,cellNodes,faceNodes);
262 // intersect a son with the unit tetra
267 // create the face key
268 TriangleFaceKey key = TriangleFaceKey(faceNodes[0], faceNodes[1], faceNodes[2]);
270 // calculate the triangle if needed
271 if(_volumes.find(key) == _volumes.end())
273 TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[1]], _nodes[faceNodes[2]]);
274 calculateVolume(tri, key);
275 totalVolume += _volumes[key];
277 baryCalculator.addSide( tri );
279 // count negative as face has reversed orientation
280 totalVolume -= _volumes[key];
287 // simple triangulation of faces along a diagonal :
298 //? not sure if this always works
300 // calculate the triangles if needed
302 // local nodes 1, 2, 3
303 TriangleFaceKey key1 = TriangleFaceKey(faceNodes[0], faceNodes[1], faceNodes[2]);
304 if(_volumes.find(key1) == _volumes.end())
306 TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[1]], _nodes[faceNodes[2]]);
307 calculateVolume(tri, key1);
308 totalVolume += _volumes[key1];
310 // count negative as face has reversed orientation
311 totalVolume -= _volumes[key1];
314 // local nodes 1, 3, 4
315 TriangleFaceKey key2 = TriangleFaceKey(faceNodes[0], faceNodes[2], faceNodes[3]);
316 if(_volumes.find(key2) == _volumes.end())
318 TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[2]], _nodes[faceNodes[3]]);
319 calculateVolume(tri, key2);
320 totalVolume += _volumes[key2];
324 // count negative as face has reversed orientation
325 totalVolume -= _volumes[key2];
332 int nbTria = nbFaceNodes - 2; // split polygon into nbTria triangles
333 for ( int iTri = 0; iTri < nbTria; ++iTri )
335 TriangleFaceKey key = TriangleFaceKey(faceNodes[0], faceNodes[1+iTri], faceNodes[2+iTri]);
336 if(_volumes.find(key) == _volumes.end())
338 TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[1+iTri]], _nodes[faceNodes[2+iTri]]);
339 calculateVolume(tri, key);
340 totalVolume += _volumes[key];
344 totalVolume -= _volumes[key];
351 std::cout << "+++ Error : Only elements with triangular and quadratilateral faces are supported at the moment." << std::endl;
358 baryCalculator.getBary( baryCentre );
359 _t->reverseApply( baryCentre, baryCentre );
363 // reset if it is very small to keep the matrix sparse
364 // is this a good idea?
365 if(epsilonEqual(totalVolume, 0.0, SPARSE_TRUNCATION_LIMIT))
370 LOG(2, "Volume = " << totalVolume << ", det= " << _t->determinant());
372 // NB : fault in article, Grandy, [8] : it is the determinant of the inverse transformation
373 // that should be used (which is equivalent to dividing by the determinant)
374 return std::fabs(1.0 / _t->determinant() * totalVolume) ;
378 * Calculates the intersection surface of two coplanar triangles.
380 * @param palneNormal normal of the plane for the first triangle
381 * @param planeConstant constant of the equation of the plane for the first triangle
382 * @param p1 coordinates of the first node of the first triangle
383 * @param p2 coordinates of the second node of the first triangle
384 * @param p3 coordinates of the third node of the first triangle
385 * @param p4 coordinates of the first node of the second triangle
386 * @param p5 coordinates of the second node of the second triangle
387 * @param p6 coordinates of the third node of the second triangle
388 * @param dimCaracteristic characteristic size of the meshes containing the triangles
389 * @param precision precision for double float data used for comparison
391 template<class MyMeshType>
392 double SplitterTetra<MyMeshType>::CalculateIntersectionSurfaceOfCoplanarTriangles(const double *const planeNormal,
393 const double planeConstant,
394 const double *const p1, const double *const p2, const double *const p3,
395 const double *const p4, const double *const p5, const double *const p6,
396 const double dimCaracteristic, const double precision)
398 typedef typename MyMeshType::MyConnType ConnType;
399 typedef double Vect2[2];
400 typedef double Vect3[3];
401 typedef double Triangle2[3][2];
403 const double *const tri0[3] = {p1, p2, p3};
404 const double *const tri1[3] = {p4, p5, p6};
406 // Plane of the first triangle defined by the normal of the triangle and the constant
407 // Project triangles onto coordinate plane most aligned with plane normal
409 double fmax = std::abs(planeNormal[0]);
410 double absMax = std::abs(planeNormal[1]);
416 absMax = std::abs(planeNormal[2]);
422 Triangle2 projTri0, projTri1;
427 // Project onto yz-plane.
428 for (i = 0; i < 3; ++i)
430 projTri0[i][0] = tri0[i][1];
431 projTri0[i][1] = tri0[i][2];
432 projTri1[i][0] = tri1[i][1];
433 projTri1[i][1] = tri1[i][2];
436 else if (maxNormal == 1)
438 // Project onto xz-plane.
439 for (i = 0; i < 3; ++i)
441 projTri0[i][0] = tri0[i][0];
442 projTri0[i][1] = tri0[i][2];
443 projTri1[i][0] = tri1[i][0];
444 projTri1[i][1] = tri1[i][2];
449 // Project onto xy-plane.
450 for (i = 0; i < 3; ++i)
452 projTri0[i][0] = tri0[i][0];
453 projTri0[i][1] = tri0[i][1];
454 projTri1[i][0] = tri1[i][0];
455 projTri1[i][1] = tri1[i][1];
459 // 2D triangle intersection routines require counterclockwise ordering.
463 for (int ii = 0; ii < 2; ++ii)
465 edge0[ii] = projTri0[1][ii] - projTri0[0][ii];
466 edge1[ii] = projTri0[2][ii] - projTri0[0][ii];
468 if ((edge0[0] * edge1[1] - edge0[1] * edge1[0]) < (double) 0.)
470 // Triangle is clockwise, reorder it.
471 for (int ii = 0; ii < 2; ++ii)
473 save[ii] = projTri0[1][ii];
474 projTri0[1][ii] = projTri0[2][ii];
475 projTri0[2][ii] = save[ii];
479 for (int ii = 0; ii < 2; ++ii)
481 edge0[ii] = projTri1[1][ii] - projTri1[0][ii];
482 edge1[ii] = projTri1[2][ii] - projTri1[0][ii];
484 if ((edge0[0] * edge1[1] - edge0[1] * edge1[0]) < (double) 0.)
486 // Triangle is clockwise, reorder it.
487 for (int ii = 0; ii < 2; ++ii)
489 save[ii] = projTri1[1][ii];
490 projTri1[1][ii] = projTri1[2][ii];
491 projTri1[2][ii] = save[ii];
495 std::vector<double> inter2;
496 intersec_de_triangle(projTri0[0], projTri0[1], projTri0[2],
497 projTri1[0], projTri1[1], projTri1[2],
499 dimCaracteristic, precision);
500 ConnType nb_inter=((ConnType)inter2.size())/2;
502 if(nb_inter >3) inter2=reconstruct_polygon(inter2);
505 std::vector<double> inter3;
506 inter3.resize(3 * nb_inter);
507 // Map 2D intersections back to the 3D triangle space.
510 double invNX = ((double) 1.) / planeNormal[0];
511 for (i = 0; i < nb_inter; i++)
513 inter3[3 * i + 1] = inter2[2 * i];
514 inter3[3 * i + 2] = inter2[2 * i + 1];
515 inter3[3 * i] = invNX * (planeConstant - planeNormal[1] * inter3[3 * i + 1] - planeNormal[2] * inter3[3 * i + 2]);
518 else if (maxNormal == 1)
520 double invNY = ((double) 1.) / planeNormal[1];
521 for (i = 0; i < nb_inter; i++)
523 inter3[3 * i] = inter2[2 * i];
524 inter3[3 * i + 2] = inter2[2 * i + 1];
525 inter3[3 * i + 1] = invNY * (planeConstant - planeNormal[0] * inter3[3 * i] - planeNormal[2] * inter3[3 * i + 2]);
530 double invNZ = ((double) 1.) / planeNormal[2];
531 for (i = 0; i < nb_inter; i++)
533 inter3[3 * i] = inter2[2 * i];
534 inter3[3 * i + 1] = inter2[2 * i + 1];
535 inter3[3 * i + 2] = invNZ * (planeConstant - planeNormal[0] * inter3[3 * i] - planeNormal[1] * inter3[3 * i + 1]);
538 surface = polygon_area<3>(inter3);
544 * Determine if a face is coplanar with a triangle.
545 * The first face is characterized by the equation of her plane
547 * @param palneNormal normal of the plane for the first triangle
548 * @param planeConstant constant of the equation of the plane for the first triangle
549 * @param coordsFace coordinates of the triangle face
550 * @param precision precision for double float data used for comparison
552 template<class MyMeshType>
553 bool SplitterTetra<MyMeshType>::IsFacesCoplanar(const double *const planeNormal,
554 const double planeConstant,
555 const double *const *const coordsFace,
556 const double precision)
558 // Compute the signed distances of triangle vertices to the plane. Use an epsilon-thick plane test.
559 // For faces not left
561 for (int i = 0; i < 3; ++i)
563 const double distance = dot(planeNormal, coordsFace[i]) - planeConstant;
564 if (epsilonEqual(distance, precision))
576 * Calculates the surface of intersection of a polygon face in the source mesh and a cell of the target mesh.
577 * It first calculates the transformation that takes the target tetrahedron into the unit tetrahedron. After that, the
578 * faces of the source element are triangulated and the calculated transformation is applied
580 * The algorithm is based on the algorithm of Grandy used in intersectSourceCell to compute
581 * the volume of intersection of two cell elements.
582 * The case with a source face colinear to one of the face of tetrahedrons is taking into account:
583 * the contribution of the face must not be counted two times.
585 * The class will cache the intermediary calculations of transformed nodes of source faces and surfaces associated
586 * with triangulated faces to avoid having to recalculate these.
588 * @param polyType type of the polygon source face
589 * @param polyNodesNbr number of the nodes of the polygon source face
590 * @param polyNodes numbers of the nodes of the polygon source face
591 * @param polyCoords coordinates of the nodes of the polygon source face
592 * @param polyCoords coordinates of the nodes of the polygon source face
593 * @param dimCaracteristic characteristic size of the meshes containing the triangles
594 * @param precision precision for double float data used for comparison
595 * @param listOfTetraFacesTreated list of tetra faces treated
596 * @param listOfTetraFacesColinear list of tetra faces colinear with the polygon source faces
598 template<class MyMeshType>
599 double SplitterTetra<MyMeshType>::intersectSourceFace(const NormalizedCellType polyType,
600 const int polyNodesNbr,
601 const int *const polyNodes,
602 const double *const *const polyCoords,
603 const double dimCaracteristic,
604 const double precision,
605 std::multiset<TriangleFaceKey>& listOfTetraFacesTreated,
606 std::set<TriangleFaceKey>& listOfTetraFacesColinear)
608 typedef typename MyMeshType::MyConnType ConnType;
610 double totalSurface = 0.0;
612 // check if we have planar tetra element
613 if(_t->determinant() == 0.0)
616 LOG(2, "Planar tetra -- volume 0");
620 // halfspace filtering
621 bool isOutside[8] = {true, true, true, true, true, true, true, true};
622 bool isStrictlyOutside[8] = {true, true, true, true, true, true, true, true};
623 bool isTargetStrictlyOutside = false;
624 bool isTargetOutside = false;
626 // calculate the coordinates of the nodes
627 for(int i = 0;i<(int)polyNodesNbr;++i)
629 const int globalNodeNum = polyNodes[i];
630 if(_nodes.find(globalNodeNum) == _nodes.end())
632 //for(HashMap< int , double* >::iterator iter3=_nodes.begin();iter3!=_nodes.end();iter3++)
633 // std::cout << (*iter3).first << " ";
634 //std::cout << std::endl << "*** " << globalNodeNum << std::endl;
635 calculateNode2(globalNodeNum, polyCoords[i]);
638 checkIsStrictlyOutside(_nodes[globalNodeNum], isStrictlyOutside, precision);
639 checkIsOutside(_nodes[globalNodeNum], isOutside, precision);
642 // halfspace filtering check
643 // NB : might not be beneficial for caching of triangles
644 for(int i = 0; i < 8; ++i)
646 if(isStrictlyOutside[i])
648 isTargetStrictlyOutside = true;
651 else if (isOutside[i])
653 isTargetOutside = true;
657 if (!isTargetStrictlyOutside)
662 // Faces are parallel
663 const int tetraFacesNodesConn[4][3] = {
668 double planeNormal[3];
669 for (int iTetraFace = 0; iTetraFace < 4; ++iTetraFace)
671 const int * const tetraFaceNodesConn = tetraFacesNodesConn[iTetraFace];
672 TriangleFaceKey key = TriangleFaceKey(_conn[tetraFaceNodesConn[0]],
673 _conn[tetraFaceNodesConn[1]],
674 _conn[tetraFaceNodesConn[2]]);
675 if (listOfTetraFacesTreated.find(key) == listOfTetraFacesTreated.end())
677 const double * const coordsTetraTriNode1 = _coords + tetraFaceNodesConn[0] * MyMeshType::MY_SPACEDIM;
678 const double * const coordsTetraTriNode2 = _coords + tetraFaceNodesConn[1] * MyMeshType::MY_SPACEDIM;
679 const double * const coordsTetraTriNode3 = _coords + tetraFaceNodesConn[2] * MyMeshType::MY_SPACEDIM;
680 calculateNormalForTria(coordsTetraTriNode1, coordsTetraTriNode2, coordsTetraTriNode3, planeNormal);
681 const double normOfTetraTriNormal = norm(planeNormal);
682 if (epsilonEqual(normOfTetraTriNormal, 0.))
684 for (int i = 0; i < 3; ++i)
691 const double invNormOfTetraTriNormal = 1. / normOfTetraTriNormal;
692 for (int i = 0; i < 3; ++i)
694 planeNormal[i] *= invNormOfTetraTriNormal;
697 double planeConstant = dot(planeNormal, coordsTetraTriNode1);
698 if (IsFacesCoplanar(planeNormal, planeConstant, polyCoords, precision))
700 int nbrPolyTri = polyNodesNbr - 2; // split polygon into nbrPolyTri triangles
701 for (int iTri = 0; iTri < nbrPolyTri; ++iTri)
703 double volume = CalculateIntersectionSurfaceOfCoplanarTriangles(planeNormal,
706 polyCoords[1 + iTri],
707 polyCoords[2 + iTri],
713 if (!epsilonEqual(volume, 0.))
715 totalSurface += volume;
716 listOfTetraFacesColinear.insert(key);
721 listOfTetraFacesTreated.insert(key);
726 // intersect a son with the unit tetra
731 // create the face key
732 TriangleFaceKey key = TriangleFaceKey(polyNodes[0], polyNodes[1], polyNodes[2]);
734 // calculate the triangle if needed
735 if (_volumes.find(key) == _volumes.end())
737 TransformedTriangle tri(_nodes[polyNodes[0]], _nodes[polyNodes[1]], _nodes[polyNodes[2]]);
738 calculateSurface(tri, key);
739 totalSurface += _volumes[key];
743 // count negative as face has reversed orientation
744 totalSurface -= _volumes[key];
751 // simple triangulation of faces along a diagonal :
762 //? not sure if this always works
764 // calculate the triangles if needed
766 // local nodes 1, 2, 3
767 TriangleFaceKey key1 = TriangleFaceKey(polyNodes[0], polyNodes[1], polyNodes[2]);
768 if (_volumes.find(key1) == _volumes.end())
770 TransformedTriangle tri(_nodes[polyNodes[0]], _nodes[polyNodes[1]], _nodes[polyNodes[2]]);
771 calculateSurface(tri, key1);
772 totalSurface += _volumes[key1];
776 // count negative as face has reversed orientation
777 totalSurface -= _volumes[key1];
780 // local nodes 1, 3, 4
781 TriangleFaceKey key2 = TriangleFaceKey(polyNodes[0], polyNodes[2], polyNodes[3]);
782 if (_volumes.find(key2) == _volumes.end())
784 TransformedTriangle tri(_nodes[polyNodes[0]], _nodes[polyNodes[2]], _nodes[polyNodes[3]]);
785 calculateSurface(tri, key2);
786 totalSurface += _volumes[key2];
790 // count negative as face has reversed orientation
791 totalSurface -= _volumes[key2];
798 int nbrPolyTri = polyNodesNbr - 2; // split polygon into nbrPolyTri triangles
799 for (int iTri = 0; iTri < nbrPolyTri; ++iTri)
801 TriangleFaceKey key = TriangleFaceKey(polyNodes[0], polyNodes[1 + iTri], polyNodes[2 + iTri]);
802 if (_volumes.find(key) == _volumes.end())
804 TransformedTriangle tri(_nodes[polyNodes[0]], _nodes[polyNodes[1 + iTri]],
805 _nodes[polyNodes[2 + iTri]]);
806 calculateSurface(tri, key);
807 totalSurface += _volumes[key];
811 totalSurface -= _volumes[key];
819 << "+++ Error : Only elements with triangular and quadratilateral faces are supported at the moment."
827 // reset if it is very small to keep the matrix sparse
828 // is this a good idea?
829 if(epsilonEqual(totalSurface, 0.0, SPARSE_TRUNCATION_LIMIT))
834 LOG(2, "Volume = " << totalSurface << ", det= " << _t->determinant());
840 * Calculates the volume of intersection of this tetrahedron with another one.
842 template<class MyMeshType>
843 double SplitterTetra<MyMeshType>::intersectTetra(const double** tetraCorners)
845 //{ could be done on outside?
846 // check if we have planar tetra element
847 if(_t->determinant() == 0.0)
850 LOG(2, "Planar tetra -- volume 0");
854 const unsigned nbOfNodes4Type=4;
855 // halfspace filtering
856 bool isOutside[8] = {true, true, true, true, true, true, true, true};
857 bool isTargetOutside = false;
859 // calculate the transformed coordinates of the nodes
860 double nodes[nbOfNodes4Type][3];
861 for(int i = 0;i<(int)nbOfNodes4Type;++i)
863 _t->apply(nodes[i], tetraCorners[i]);
864 checkIsOutside(nodes[i], isOutside);
867 // halfspace filtering check
868 // NB : might not be beneficial for caching of triangles
869 for(int i = 0; i < 8; ++i)
873 isTargetOutside = true;
877 double totalVolume = 0.0;
881 const CellModel& cellModelCell=CellModel::GetCellModel(NORM_TETRA4);
882 int cellNodes[4] = { 0, 1, 2, 3 }, faceNodes[3];
884 for(unsigned ii = 0 ; ii < 4 ; ++ii)
886 cellModelCell.fillSonCellNodalConnectivity(ii,cellNodes,faceNodes);
888 TransformedTriangle tri(nodes[faceNodes[0]], nodes[faceNodes[1]], nodes[faceNodes[2]]);
889 double vol = tri.calculateIntersectionVolume();
893 // reset if it is very small to keep the matrix sparse
894 // is this a good idea?
895 if(epsilonEqual(totalVolume, 0.0, SPARSE_TRUNCATION_LIMIT))
900 LOG(2, "Volume = " << totalVolume << ", det= " << _t->determinant());
902 // NB : fault in article, Grandy, [8] : it is the determinant of the inverse transformation
903 // that should be used (which is equivalent to dividing by the determinant)
904 return std::fabs(1.0 / _t->determinant() * totalVolume) ;
907 ////////////////////////////////////////////////////////
909 template<class MyMeshTypeT, class MyMeshTypeS>
910 SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::SplitterTetra2(const MyMeshTypeT& targetMesh, const MyMeshTypeS& srcMesh, SplittingPolicy policy)
911 :_target_mesh(targetMesh),_src_mesh(srcMesh),_splitting_pol(policy)
915 template<class MyMeshTypeT, class MyMeshTypeS>
916 SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::~SplitterTetra2()
921 template<class MyMeshTypeT, class MyMeshTypeS>
922 void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::releaseArrays()
924 // free potential sub-mesh nodes that have been allocated
925 typename MyMeshTypeT::MyConnType nbOfNodesT = _node_ids.size();// Issue 0020634.
926 if((int)_nodes.size()>=/*8*/nbOfNodesT)
928 std::vector<const double*>::iterator iter = _nodes.begin() + /*8*/nbOfNodesT;
929 while(iter != _nodes.end())
939 * @param targetCell in C mode.
940 * @param tetra is the output result tetra containers.
942 template<class MyMeshTypeT, class MyMeshTypeS>
943 void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::splitTargetCell(typename MyMeshTypeT::MyConnType targetCell,
944 typename MyMeshTypeT::MyConnType nbOfNodesT,
945 typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
947 typedef typename MyMeshTypeT::MyConnType ConnType;
948 const NumberingPolicy numPol=MyMeshTypeT::My_numPol;
949 const int numTetra = static_cast<int>(_splitting_pol);
955 const double *nodes[4];
957 for(int node = 0; node < 4 ; ++node)
959 nodes[node]=getCoordsOfNode2(node, OTT<ConnType,numPol>::indFC(targetCell),_target_mesh,conn[node]);
961 std::copy(conn,conn+4,_node_ids.begin());
962 SplitterTetra<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_src_mesh, nodes,conn);
966 // Issue 0020634. To pass nbOfNodesT to calculateSubNodes (don't want to add an arg)
967 _node_ids.resize(nbOfNodesT);
969 // pre-calculate nodes
970 calculateSubNodes(_target_mesh, OTT<ConnType,numPol>::indFC(targetCell));
972 tetra.reserve(numTetra);
973 _nodes.reserve(30); // we never have more than this
975 switch ( nbOfNodesT )
979 switch(_splitting_pol)
983 const int subZone[8] = { 0, 1, 2, 3, 4, 5, 6, 7 };
984 fiveSplit(subZone,tetra);
990 const int subZone[8] = { 0, 1, 2, 3, 4, 5, 6, 7 };
991 sixSplit(subZone,tetra);
997 calculateGeneral24Tetra(tetra);
1003 calculateGeneral48Tetra(tetra);
1018 splitConvex(targetCell, tetra);
1024 * Splits the hexahedron into five tetrahedra.
1025 * This method adds five SplitterTetra objects to the vector tetra.
1027 * @param subZone the local node numbers corresponding to the hexahedron corners - these are mapped onto {0,..,7}. Providing this allows the
1028 * splitting to be reused on the subzones of the GENERAL_* types of splitting
1030 template<class MyMeshTypeT, class MyMeshTypeS>
1031 void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::fiveSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
1033 // create tetrahedra
1034 for(int i = 0; i < 5; ++i)
1036 const double* nodes[4];
1038 for(int j = 0; j < 4; ++j)
1040 conn[j] = subZone[ SPLIT_NODES_5[4*i+j] ];
1041 nodes[j] = getCoordsOfSubNode(conn[j]);
1043 SplitterTetra<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_src_mesh, nodes,conn);
1049 * Splits the hexahedron into six tetrahedra.
1050 * This method adds six SplitterTetra objects to the vector tetra.
1052 * @param subZone the local node numbers corresponding to the hexahedron corners - these are mapped onto {0,..,7}. Providing this allows the
1053 * splitting to be reused on the subzones of the GENERAL_* types of splitting
1055 template<class MyMeshTypeT, class MyMeshTypeS>
1056 void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::sixSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
1058 for(int i = 0; i < 6; ++i)
1060 const double* nodes[4];
1062 for(int j = 0; j < 4; ++j)
1064 conn[j] = subZone[SPLIT_NODES_6[4*i+j]];
1065 nodes[j] = getCoordsOfSubNode(conn[j]);
1067 SplitterTetra<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_src_mesh, nodes,conn);
1073 * Splits the hexahedron into 24 tetrahedra.
1074 * The splitting is done by combining the barycenter of the tetrahedron, the barycenter of each face
1075 * and the nodes of each edge of the face. This creates 6 faces * 4 edges / face = 24 tetrahedra.
1076 * The submesh nodes introduced are the barycenters of the faces and the barycenter of the cell.
1079 template<class MyMeshTypeT, class MyMeshTypeS>
1080 void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calculateGeneral24Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
1082 // The two nodes of the original mesh cell used in each tetrahedron.
1083 // The tetrahedra all have nodes (cellCenter, faceCenter, edgeNode1, edgeNode2)
1084 // For the correspondance of the nodes, see the GENERAL_48_SUB_NODES table in calculateSubNodes
1086 // nodes to use for tetrahedron
1087 const double* nodes[4];
1089 // get the cell center
1091 nodes[0] = getCoordsOfSubNode(conn[0]);
1093 for(int faceCenterNode = 8; faceCenterNode < 14; ++faceCenterNode)
1095 // get the face center
1096 conn[1] = faceCenterNode;
1097 nodes[1] = getCoordsOfSubNode(conn[1]);
1098 for(int j = 0; j < 4; ++j)
1100 const int row = 4*(faceCenterNode - 8) + j;
1101 conn[2] = TETRA_EDGES_GENERAL_24[2*row];
1102 conn[3] = TETRA_EDGES_GENERAL_24[2*row + 1];
1103 nodes[2] = getCoordsOfSubNode(conn[2]);
1104 nodes[3] = getCoordsOfSubNode(conn[3]);
1106 SplitterTetra<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_src_mesh, nodes, conn);
1114 * Splits the hexahedron into 48 tetrahedra.
1115 * The splitting is done by introducing the midpoints of all the edges
1116 * and the barycenter of the element as submesh nodes. The 8 hexahedral subzones thus defined
1117 * are then split into 6 tetrahedra each, as in Grandy, p. 449. The division of the subzones
1118 * is done by calling sixSplit().
1121 template<class MyMeshTypeT, class MyMeshTypeS>
1122 void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calculateGeneral48Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
1124 // Define 8 hexahedral subzones as in Grandy, p449
1125 // the values correspond to the nodes that correspond to nodes 1,2,3,4,5,6,7,8 in the subcell
1126 // For the correspondance of the nodes, see the GENERAL_48_SUB_NODES table in calculateSubNodes
1127 static const int subZones[64] =
1129 0,8,21,12,9,20,26,22,
1130 8,1,13,21,20,10,23,26,
1131 12,21,16,3,22,26,25,17,
1132 21,13,2,16,26,23,18,25,
1133 9,20,26,22,4,11,24,14,
1134 20,10,23,26,11,5,15,24,
1135 22,26,25,17,14,24,19,7,
1136 26,23,18,25,24,15,6,19
1139 for(int i = 0; i < 8; ++i)
1141 sixSplit(&subZones[8*i],tetra);
1146 * Splits the NORM_PYRA5 into 2 tetrahedra.
1148 template<class MyMeshTypeT, class MyMeshTypeS>
1149 void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::splitPyram5(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
1151 static const int SPLIT_PYPA5[2][4] =
1161 // create tetrahedra
1162 const double* nodes[4];
1164 for(int i = 0; i < 2; ++i)
1166 for(int j = 0; j < 4; ++j)
1167 nodes[j] = getCoordsOfSubNode2(SPLIT_PYPA5[i][j],conn[j]);
1168 SplitterTetra<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_src_mesh, nodes,conn);
1174 * Splits a convex cell into tetrahedra.
1176 template<class MyMeshTypeT, class MyMeshTypeS>
1177 void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::splitConvex(typename MyMeshTypeT::MyConnType targetCell,
1178 typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
1180 // Each face of a cell is split into triangles and
1181 // each of triangles and a cell barycenter form a tetrahedron.
1183 typedef typename MyMeshTypeT::MyConnType ConnType;
1184 const NumberingPolicy numPol=MyMeshTypeT::My_numPol;
1186 // get type of cell and nb of cell nodes
1187 NormalizedCellType normCellType=_target_mesh.getTypeOfElement(OTT<ConnType,numPol>::indFC(targetCell));
1188 const CellModel& cellModelCell=CellModel::GetCellModel(normCellType);
1189 unsigned nbOfCellNodes=cellModelCell.isDynamic() ? _target_mesh.getNumberOfNodesOfElement(OTT<ConnType,numPol>::indFC(targetCell)) : cellModelCell.getNumberOfNodes();
1191 // get nb of cell sons (faces)
1192 const ConnType* rawCellConn = _target_mesh.getConnectivityPtr() + OTT<ConnType,numPol>::conn2C( _target_mesh.getConnectivityIndexPtr()[ targetCell ]);
1193 const int rawNbCellNodes = _target_mesh.getConnectivityIndexPtr()[ targetCell+1 ] - _target_mesh.getConnectivityIndexPtr()[ targetCell ];
1194 unsigned nbOfSons = cellModelCell.getNumberOfSons2(rawCellConn, rawNbCellNodes);
1196 // indices of nodes of a son
1197 static std::vector<int> allNodeIndices; // == 0,1,2,...,nbOfCellNodes-1
1198 while ( allNodeIndices.size() < nbOfCellNodes )
1199 allNodeIndices.push_back( allNodeIndices.size() );
1200 std::vector<int> classicFaceNodes(4);
1201 int* faceNodes = cellModelCell.isDynamic() ? &allNodeIndices[0] : &classicFaceNodes[0];
1203 // nodes of tetrahedron
1205 const double* nodes[4];
1206 nodes[3] = getCoordsOfSubNode2( nbOfCellNodes,conn[3]); // barycenter
1208 for(unsigned ii = 0 ; ii < nbOfSons; ++ii)
1210 // get indices of son's nodes: it's just next portion of allNodeIndices for polyhedron
1211 // and some of allNodeIndices accodring to cell model for a classsic cell
1212 unsigned nbFaceNodes = cellModelCell.getNumberOfNodesConstituentTheSon2(ii, rawCellConn, rawNbCellNodes);
1213 if ( normCellType != NORM_POLYHED )
1214 cellModelCell.fillSonCellNodalConnectivity(ii,&allNodeIndices[0],faceNodes);
1216 int nbTetra = nbFaceNodes - 2; // split polygon into nbTetra triangles
1218 // create tetrahedra
1219 for(int i = 0; i < nbTetra; ++i)
1221 nodes[0] = getCoordsOfSubNode2( faceNodes[0], conn[0]);
1222 nodes[1] = getCoordsOfSubNode2( faceNodes[1+i],conn[1]);
1223 nodes[2] = getCoordsOfSubNode2( faceNodes[2+i],conn[2]);
1224 SplitterTetra<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_src_mesh, nodes,conn);
1228 if ( normCellType == NORM_POLYHED )
1229 faceNodes += nbFaceNodes; // go to the next face
1234 * Precalculates all the nodes.
1235 * Retrieves the mesh nodes and allocates the necessary sub-mesh
1236 * nodes according to the splitting policy used.
1237 * This method is meant to be called once by the constructor.
1239 * @param targetMesh the target mesh
1240 * @param targetCell the global number of the cell that the object represents, in targetMesh mode.
1241 * @param policy the splitting policy of the object
1244 template<class MyMeshTypeT, class MyMeshTypeS>
1245 void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calculateSubNodes(const MyMeshTypeT& targetMesh, typename MyMeshTypeT::MyConnType targetCell)
1247 // retrieve real mesh nodes
1249 typename MyMeshTypeT::MyConnType nbOfNodesT = _node_ids.size();// Issue 0020634. _node_ids.resize(8);
1250 for(int node = 0; node < nbOfNodesT ; ++node)
1252 // calculate only normal nodes
1253 _nodes.push_back(getCoordsOfNode2(node, targetCell, targetMesh,_node_ids[node]));
1256 switch ( nbOfNodesT )
1260 // create sub-mesh nodes if needed
1261 switch(_splitting_pol)
1265 for(int i = 0; i < 7; ++i)
1267 double* barycenter = new double[3];
1268 calcBarycenter(4, barycenter, &GENERAL_24_SUB_NODES[4*i]);
1269 _nodes.push_back(barycenter);
1276 // Each sub-node is the barycenter of two other nodes.
1277 // For the edges, these lie on the original mesh.
1278 // For the faces, these are the edge sub-nodes.
1279 // For the cell these are two face sub-nodes.
1280 static const int GENERAL_48_SUB_NODES[38] =
1282 0,1, // sub-node 9 (edge)
1283 0,4, // sub-node 10 (edge)
1284 1,5, // sub-node 11 (edge)
1285 4,5, // sub-node 12 (edge)
1286 0,3, // sub-node 13 (edge)
1287 1,2, // sub-node 14 (edge)
1288 4,7, // sub-node 15 (edge)
1289 5,6, // sub-node 16 (edge)
1290 2,3, // sub-node 17 (edge)
1291 3,7, // sub-node 18 (edge)
1292 2,6, // sub-node 19 (edge)
1293 6,7, // sub-node 20 (edge)
1294 8,11, // sub-node 21 (face)
1295 12,13, // sub-node 22 (face)
1296 9,17, // sub-node 23 (face)
1297 10,18, // sub-node 24 (face)
1298 14,15, // sub-node 25 (face)
1299 16,19, // sub-node 26 (face)
1300 20,25 // sub-node 27 (cell)
1303 for(int i = 0; i < 19; ++i)
1305 double* barycenter = new double[3];
1306 calcBarycenter(2, barycenter, &GENERAL_48_SUB_NODES[2*i]);
1307 _nodes.push_back(barycenter);
1316 case 5: // NORM_PYRA5
1319 default: // convex 3d cell
1321 // add barycenter of a cell
1322 std::vector<int> allIndices(nbOfNodesT);
1323 for ( int i = 0; i < nbOfNodesT; ++i ) allIndices[i] = i;
1324 double* barycenter = new double[3];
1325 calcBarycenter(nbOfNodesT, barycenter, &allIndices[0]);
1326 _nodes.push_back(barycenter);