1 // Copyright (C) 2007-2024 CEA, EDF
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, or (at your option) any later version.
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 namespace INTERP_KERNEL
41 template<class MyMeshType>
42 const double SplitterTetra<MyMeshType>::SPARSE_TRUNCATION_LIMIT=1.0e-14;
45 * output is expected to be allocated with 24*sizeof(void*) in order to store the 24 tetras.
46 * These tetras have to be deallocated.
48 template<class MyMeshType>
49 void SplitterTetra<MyMeshType>::splitIntoDualCells(SplitterTetra<MyMeshType> **output)
52 const double *tmp2[4]={tmp,tmp+3,tmp+6,tmp+9};
53 typename MyMeshType::MyConnType conn[4]={-1,-1,-1,-1};
56 splitMySelfForDual(tmp,i,conn[0]);
57 output[i]=new SplitterTetra<MyMeshType>(_src_mesh,tmp2,conn);
62 * SplitterTetra class computes for a list of cell ids of a given mesh \a srcMesh (badly named) the intersection with a
63 * single TETRA4 cell given by \a tetraCorners (of length 4) and \a nodesId (of length 4 too). \a nodedIds is given only to establish
64 * if a partial computation of a triangle has already been performed (to increase performance).
66 * The \a srcMesh can contain polyhedron cells.
69 * Constructor creating object from the four corners of the tetrahedron.
71 * @param srcMesh mesh containing the source elements
72 * @param tetraCorners array of four pointers to double[3] arrays containing the coordinates of the
73 * corners of the tetrahedron
75 template<class MyMeshType>
76 SplitterTetra<MyMeshType>::SplitterTetra(const MyMeshType& srcMesh, const double** tetraCorners, const typename MyMeshType::MyConnType *nodesId)
77 : _t(0), _src_mesh(srcMesh)
79 std::copy(nodesId,nodesId+4,_conn);
80 _coords[0]=tetraCorners[0][0]; _coords[1]=tetraCorners[0][1]; _coords[2]=tetraCorners[0][2];
81 _coords[3]=tetraCorners[1][0]; _coords[4]=tetraCorners[1][1]; _coords[5]=tetraCorners[1][2];
82 _coords[6]=tetraCorners[2][0]; _coords[7]=tetraCorners[2][1]; _coords[8]=tetraCorners[2][2];
83 _coords[9]=tetraCorners[3][0]; _coords[10]=tetraCorners[3][1]; _coords[11]=tetraCorners[3][2];
84 // create the affine transform
85 _t=new TetraAffineTransform(_coords);
89 * SplitterTetra class computes for a list of cell ids of a given mesh \a srcMesh (badly named) the intersection with a
90 * single TETRA4 cell given by \a tetraCorners (of length 4) and \a nodesId (of length 4 too). \a nodedIds is given only to establish
91 * if a partial computation of a triangle has already been performed (to increase performance).
93 * The \a srcMesh can contain polyhedron cells.
96 * Constructor creating object from the four corners of the tetrahedron.
98 * \param [in] srcMesh mesh containing the source elements
99 * \param [in] tetraCorners array 4*3 doubles containing corners of input tetrahedron (P0X,P0Y,P0Y,P1X,P1Y,P1Z,P2X,P2Y,P2Z,P3X,P3Y,P3Z).
101 template<class MyMeshType>
102 SplitterTetra<MyMeshType>::SplitterTetra(const MyMeshType& srcMesh, const double tetraCorners[12], const ConnType *conn): _t(0),_src_mesh(srcMesh)
105 { _conn[0]=0; _conn[1]=1; _conn[2]=2; _conn[3]=3; }
107 { _conn[0]=conn[0]; _conn[1]=conn[1]; _conn[2]=conn[2]; _conn[3]=conn[3]; }
108 _coords[0]=tetraCorners[0]; _coords[1]=tetraCorners[1]; _coords[2]=tetraCorners[2]; _coords[3]=tetraCorners[3]; _coords[4]=tetraCorners[4]; _coords[5]=tetraCorners[5];
109 _coords[6]=tetraCorners[6]; _coords[7]=tetraCorners[7]; _coords[8]=tetraCorners[8]; _coords[9]=tetraCorners[9]; _coords[10]=tetraCorners[10]; _coords[11]=tetraCorners[11];
110 // create the affine transform
111 _t=new TetraAffineTransform(_coords);
117 * Deletes _t and the coordinates (double[3] vectors) in _nodes
120 template<class MyMeshType>
121 SplitterTetra<MyMeshType>::~SplitterTetra()
124 for(typename HashMap< ConnType, double* >::iterator iter = _nodes.begin(); iter != _nodes.end() ; ++iter)
125 delete[] iter->second;
129 * \Forget already calculated triangles, which is crucial for calculation of barycenter of intersection
131 template<class MyMeshType>
132 void SplitterTetra<MyMeshType>::clearVolumesCache()
138 * This method destroys the 4 pointers pointed by tetraCorners[0],tetraCorners[1],tetraCorners[2] and tetraCorners[3]
139 * @param i is in 0..23 included.
140 * @param output is expected to be sized of 12 in order to.
142 template<class MyMeshType>
143 void SplitterTetra<MyMeshType>::splitMySelfForDual(double* output, int i, typename MyMeshType::MyConnType& nodeId)
147 nodeId=_conn[offset];
148 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;
150 int case1=caseToTreat/2;
151 int case2=caseToTreat%2;
152 const int tab[3][2]={{1,2},{3,2},{1,3}};
153 const int *curTab=tab[case1];
154 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.;
155 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.;
156 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.;
157 std::copy(pt1,pt1+3,output+case2*3);
158 std::copy(pt0,pt0+3,output+(abs(case2-1))*3);
159 std::copy(pt2,pt2+3,output+2*3);
160 std::copy(tmp[0],tmp[0]+3,output+3*3);
164 * Calculates the volume of intersection of an element in the source mesh and the target element.
165 * It first calculates the transformation that takes the target tetrahedron into the unit tetrahedron. After that, the
166 * faces of the source element are triangulated and the calculated transformation is applied
167 * to each triangle. The algorithm of Grandy, implemented in INTERP_KERNEL::TransformedTriangle is used
168 * to calculate the contribution to the volume from each triangle. The volume returned is the sum of these contributions
169 * divided by the determinant of the transformation.
171 * The class will cache the intermediary calculations of transformed nodes of source cells and volumes associated
172 * with triangulated faces to avoid having to recalculate these.
174 * @param element global number of the source element in C mode.
176 template<class MyMeshType>
177 double SplitterTetra<MyMeshType>::intersectSourceCell(typename MyMeshType::MyConnType element,
180 typedef typename MyMeshType::MyConnType ConnType;
181 const NumberingPolicy numPol=MyMeshType::My_numPol;
182 //{ could be done on outside?
183 // check if we have planar tetra element
184 if(_t->determinant() == 0.0)
187 LOG(2, "Planar tetra -- volume 0");
192 NormalizedCellType normCellType=_src_mesh.getTypeOfElement(OTT<ConnType,numPol>::indFC(element));
193 const CellModel& cellModelCell=CellModel::GetCellModel(normCellType);
194 ConnType nbOfNodes4Type=cellModelCell.isDynamic() ? _src_mesh.getNumberOfNodesOfElement(OTT<ConnType,numPol>::indFC(element)) : cellModelCell.getNumberOfNodes();
195 // halfspace filtering
196 bool isOutside[8] = {true, true, true, true, true, true, true, true};
197 bool isTargetOutside = false;
199 // calculate the coordinates of the nodes
200 ConnType *cellNodes=new ConnType[nbOfNodes4Type];
201 for(ConnType i = 0;i<nbOfNodes4Type;++i)
203 // we could store mapping local -> global numbers too, but not sure it is worth it
204 const ConnType globalNodeNum = getGlobalNumberOfNode(i, OTT<ConnType,numPol>::indFC(element), _src_mesh);
205 cellNodes[i]=globalNodeNum;
206 if(_nodes.find(globalNodeNum) == _nodes.end())
208 //for(HashMap< int , double* >::iterator iter3=_nodes.begin();iter3!=_nodes.end();iter3++)
209 // std::cout << (*iter3).first << " ";
210 //std::cout << std::endl << "*** " << globalNodeNum << std::endl;
211 calculateNode(globalNodeNum);
213 CheckIsOutside(_nodes[globalNodeNum], isOutside);
216 // halfspace filtering check
217 // NB : might not be beneficial for caching of triangles
218 for(int i = 0; i < 8; ++i)
220 isTargetOutside = true;
222 double totalVolume = 0.0;
226 /// calculator of intersection barycentre
227 UnitTetraIntersectionBary baryCalculator( _t->determinant() < 0.);
229 // get nb of sons of a cell
230 const ConnType* rawCellConn = _src_mesh.getConnectivityPtr() + OTT<ConnType,numPol>::conn2C( _src_mesh.getConnectivityIndexPtr()[ element ]);
231 const ConnType rawNbCellNodes = _src_mesh.getConnectivityIndexPtr()[ element+1 ] - _src_mesh.getConnectivityIndexPtr()[ element ];
232 unsigned nbOfSons = cellModelCell.getNumberOfSons2(rawCellConn, rawNbCellNodes);
234 for(unsigned ii = 0 ; ii < nbOfSons; ++ii)
236 // get sons connectivity
237 NormalizedCellType faceType;
238 ConnType *faceNodes, nbFaceNodes=-1;
239 if ( cellModelCell.isDynamic() )
241 faceNodes=new ConnType[nbOfNodes4Type];
242 nbFaceNodes = cellModelCell.fillSonCellNodalConnectivity2(ii,rawCellConn,rawNbCellNodes,faceNodes,faceType);
243 for ( ConnType i = 0; i < nbFaceNodes; ++i )
244 faceNodes[i] = OTT<ConnType,numPol>::coo2C(faceNodes[i]);
248 faceType = cellModelCell.getSonType(ii);
249 assert(CellModel::GetCellModel(faceType).getDimension() == 2);
250 nbFaceNodes = cellModelCell.getNumberOfNodesConstituentTheSon(ii);
251 faceNodes = new ConnType[nbFaceNodes];
252 cellModelCell.fillSonCellNodalConnectivity(ii,cellNodes,faceNodes);
254 // intersect a son with the unit tetra
259 // create the face key
260 TriangleFaceKey key = TriangleFaceKey(faceNodes[0], faceNodes[1], faceNodes[2]);
262 // calculate the triangle if needed
263 if(_volumes.find(key) == _volumes.end())
265 TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[1]], _nodes[faceNodes[2]]);
266 calculateVolume(tri, key);
267 totalVolume += _volumes[key];
269 baryCalculator.addSide( tri );
271 // count negative as face has reversed orientation
272 totalVolume -= _volumes[key];
279 // simple triangulation of faces along a diagonal :
290 //? not sure if this always works
292 // calculate the triangles if needed
294 // local nodes 1, 2, 3
295 TriangleFaceKey key1 = TriangleFaceKey(faceNodes[0], faceNodes[1], faceNodes[2]);
296 if(_volumes.find(key1) == _volumes.end())
298 TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[1]], _nodes[faceNodes[2]]);
299 calculateVolume(tri, key1);
300 totalVolume += _volumes[key1];
302 // count negative as face has reversed orientation
303 totalVolume -= _volumes[key1];
306 // local nodes 1, 3, 4
307 TriangleFaceKey key2 = TriangleFaceKey(faceNodes[0], faceNodes[2], faceNodes[3]);
308 if(_volumes.find(key2) == _volumes.end())
310 TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[2]], _nodes[faceNodes[3]]);
311 calculateVolume(tri, key2);
312 totalVolume += _volumes[key2];
316 // count negative as face has reversed orientation
317 totalVolume -= _volumes[key2];
324 ConnType nbTria = nbFaceNodes - 2; // split polygon into nbTria triangles
325 for ( ConnType iTri = 0; iTri < nbTria; ++iTri )
327 TriangleFaceKey key = TriangleFaceKey(faceNodes[0], faceNodes[1+iTri], faceNodes[2+iTri]);
328 if(_volumes.find(key) == _volumes.end())
330 TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[1+iTri]], _nodes[faceNodes[2+iTri]]);
331 calculateVolume(tri, key);
332 totalVolume += _volumes[key];
336 totalVolume -= _volumes[key];
343 std::cout << "+++ Error : Only elements with triangular and quadratilateral faces are supported at the moment." << std::endl;
350 baryCalculator.getBary( baryCentre );
351 _t->reverseApply( baryCentre, baryCentre );
355 // reset if it is very small to keep the matrix sparse
356 // is this a good idea?
357 if(epsilonEqual(totalVolume, 0.0, SPARSE_TRUNCATION_LIMIT))
362 LOG(2, "Volume = " << totalVolume << ", det= " << _t->determinant());
364 // NB : fault in article, Grandy, [8] : it is the determinant of the inverse transformation
365 // that should be used (which is equivalent to dividing by the determinant)
366 return std::fabs(1.0 / _t->determinant() * totalVolume) ;
370 * Calculates the intersection surface of two coplanar triangles.
372 * @param planeNormal normal of the plane for the first triangle
373 * @param planeConstant constant of the equation of the plane for the first triangle
374 * @param p1 coordinates of the first node of the first triangle
375 * @param p2 coordinates of the second node of the first triangle
376 * @param p3 coordinates of the third node of the first triangle
377 * @param p4 coordinates of the first node of the second triangle
378 * @param p5 coordinates of the second node of the second triangle
379 * @param p6 coordinates of the third node of the second triangle
380 * @param dimCaracteristic characteristic size of the meshes containing the triangles
381 * @param precision precision for double float data used for comparison
383 template<class MyMeshType>
384 double SplitterTetra<MyMeshType>::CalculateIntersectionSurfaceOfCoplanarTriangles(const double *const planeNormal,
385 const double planeConstant,
386 const double *const p1, const double *const p2, const double *const p3,
387 const double *const p4, const double *const p5, const double *const p6,
388 const double dimCaracteristic, const double precision)
390 typedef typename MyMeshType::MyConnType ConnType;
391 typedef double Vect2[2];
392 typedef double Triangle2[3][2];
394 const double *const tri0[3] = {p1, p2, p3};
395 const double *const tri1[3] = {p4, p5, p6};
397 // Plane of the first triangle defined by the normal of the triangle and the constant
398 // Project triangles onto coordinate plane most aligned with plane normal
400 double fmax = std::abs(planeNormal[0]);
401 double absMax = std::abs(planeNormal[1]);
407 absMax = std::abs(planeNormal[2]);
413 Triangle2 projTri0, projTri1;
418 // Project onto yz-plane.
419 for (i = 0; i < 3; ++i)
421 projTri0[i][0] = tri0[i][1];
422 projTri0[i][1] = tri0[i][2];
423 projTri1[i][0] = tri1[i][1];
424 projTri1[i][1] = tri1[i][2];
427 else if (maxNormal == 1)
429 // Project onto xz-plane.
430 for (i = 0; i < 3; ++i)
432 projTri0[i][0] = tri0[i][0];
433 projTri0[i][1] = tri0[i][2];
434 projTri1[i][0] = tri1[i][0];
435 projTri1[i][1] = tri1[i][2];
440 // Project onto xy-plane.
441 for (i = 0; i < 3; ++i)
443 projTri0[i][0] = tri0[i][0];
444 projTri0[i][1] = tri0[i][1];
445 projTri1[i][0] = tri1[i][0];
446 projTri1[i][1] = tri1[i][1];
450 // 2D triangle intersection routines require counterclockwise ordering.
454 for (int ii = 0; ii < 2; ++ii)
456 edge0[ii] = projTri0[1][ii] - projTri0[0][ii];
457 edge1[ii] = projTri0[2][ii] - projTri0[0][ii];
459 if ((edge0[0] * edge1[1] - edge0[1] * edge1[0]) < (double) 0.)
461 // Triangle is clockwise, reorder it.
462 for (int ii = 0; ii < 2; ++ii)
464 save[ii] = projTri0[1][ii];
465 projTri0[1][ii] = projTri0[2][ii];
466 projTri0[2][ii] = save[ii];
470 for (int ii = 0; ii < 2; ++ii)
472 edge0[ii] = projTri1[1][ii] - projTri1[0][ii];
473 edge1[ii] = projTri1[2][ii] - projTri1[0][ii];
475 if ((edge0[0] * edge1[1] - edge0[1] * edge1[0]) < (double) 0.)
477 // Triangle is clockwise, reorder it.
478 for (int ii = 0; ii < 2; ++ii)
480 save[ii] = projTri1[1][ii];
481 projTri1[1][ii] = projTri1[2][ii];
482 projTri1[2][ii] = save[ii];
486 std::vector<double> inter2;
487 intersec_de_triangle(projTri0[0], projTri0[1], projTri0[2],
488 projTri1[0], projTri1[1], projTri1[2],
490 dimCaracteristic, precision);
491 ConnType nb_inter=((ConnType)inter2.size())/2;
493 if(nb_inter >3) inter2=reconstruct_polygon(inter2);
496 std::vector<double> inter3;
497 inter3.resize(3 * nb_inter);
498 // Map 2D intersections back to the 3D triangle space.
501 double invNX = ((double) 1.) / planeNormal[0];
502 for (i = 0; i < nb_inter; i++)
504 inter3[3 * i + 1] = inter2[2 * i];
505 inter3[3 * i + 2] = inter2[2 * i + 1];
506 inter3[3 * i] = invNX * (planeConstant - planeNormal[1] * inter3[3 * i + 1] - planeNormal[2] * inter3[3 * i + 2]);
509 else if (maxNormal == 1)
511 double invNY = ((double) 1.) / planeNormal[1];
512 for (i = 0; i < nb_inter; i++)
514 inter3[3 * i] = inter2[2 * i];
515 inter3[3 * i + 2] = inter2[2 * i + 1];
516 inter3[3 * i + 1] = invNY * (planeConstant - planeNormal[0] * inter3[3 * i] - planeNormal[2] * inter3[3 * i + 2]);
521 double invNZ = ((double) 1.) / planeNormal[2];
522 for (i = 0; i < nb_inter; i++)
524 inter3[3 * i] = inter2[2 * i];
525 inter3[3 * i + 1] = inter2[2 * i + 1];
526 inter3[3 * i + 2] = invNZ * (planeConstant - planeNormal[0] * inter3[3 * i] - planeNormal[1] * inter3[3 * i + 1]);
529 surface = polygon_area<3>(inter3);
535 * Determine if a face is coplanar with a triangle.
536 * The first face is characterized by the equation of her plane
538 * @param planeNormal normal of the plane for the first triangle
539 * @param planeConstant constant of the equation of the plane for the first triangle
540 * @param coordsFace coordinates of the triangle face
541 * @param precision precision for double float data used for comparison
543 template<class MyMeshType>
544 bool SplitterTetra<MyMeshType>::IsFacesCoplanar(const double *const planeNormal,
545 const double planeConstant,
546 const double *const *const coordsFace,
547 const double precision)
549 // Compute the signed distances of triangle vertices to the plane. Use an epsilon-thick plane test.
550 // For faces not left
552 for (int i = 0; i < 3; ++i)
554 const double distance = dot(planeNormal, coordsFace[i]) - planeConstant;
555 if (epsilonEqual(distance, precision))
564 * Calculates the surface of intersection of a polygon face in the source mesh and a cell of the target mesh.
565 * It first calculates the transformation that takes the target tetrahedron into the unit tetrahedron. After that, the
566 * faces of the source element are triangulated and the calculated transformation is applied
568 * The algorithm is based on the algorithm of Grandy used in intersectSourceCell to compute
569 * the volume of intersection of two cell elements.
570 * The case with a source face colinear to one of the face of tetrahedrons is taking into account:
571 * the contribution of the face must not be counted two times.
573 * The class will cache the intermediary calculations of transformed nodes of source faces and surfaces associated
574 * with triangulated faces to avoid having to recalculate these.
576 * @param polyType type of the polygon source face
577 * @param polyNodesNbr number of the nodes of the polygon source face
578 * @param polyNodes numbers of the nodes of the polygon source face
579 * @param polyCoords coordinates of the nodes of the polygon source face
580 * @param dimCaracteristic characteristic size of the meshes containing the triangles
581 * @param precision precision for double float data used for comparison
582 * @param listOfTetraFacesTreated list of tetra faces treated
583 * @param listOfTetraFacesColinear list of tetra faces colinear with the polygon source faces
585 template<class MyMeshType>
586 double SplitterTetra<MyMeshType>::intersectSourceFace(const NormalizedCellType polyType,
587 const ConnType polyNodesNbr,
588 const ConnType *const polyNodes,
589 const double *const *const polyCoords,
590 const double dimCaracteristic,
591 const double precision,
592 std::multiset<TriangleFaceKey>& listOfTetraFacesTreated,
593 std::set<TriangleFaceKey>& listOfTetraFacesColinear)
595 double totalSurface = 0.0;
597 // check if we have planar tetra element
598 if(_t->determinant() == 0.0)
601 LOG(2, "Planar tetra -- volume 0");
605 // halfspace filtering
606 bool isOutside[8] = {true, true, true, true, true, true, true, true};
607 bool isStrictlyOutside[8] = {true, true, true, true, true, true, true, true};
608 bool isTargetStrictlyOutside = false;
609 bool isTargetOutside = false;
611 // calculate the coordinates of the nodes
612 for(ConnType i = 0;i<polyNodesNbr;++i)
614 const ConnType globalNodeNum = polyNodes[i];
615 if(_nodes.find(globalNodeNum) == _nodes.end())
617 calculateNode2(globalNodeNum, polyCoords[i]);
620 CheckIsStrictlyOutside(_nodes[globalNodeNum], isStrictlyOutside, precision);
621 CheckIsOutside(_nodes[globalNodeNum], isOutside, precision);
624 // halfspace filtering check
625 // NB : might not be beneficial for caching of triangles
626 for(int i = 0; i < 8; ++i)
628 if(isStrictlyOutside[i])
630 isTargetStrictlyOutside = true;
633 else if (isOutside[i])
635 isTargetOutside = true;
639 if (!isTargetStrictlyOutside)
644 // Faces are parallel
645 const int tetraFacesNodesConn[4][3] = {
650 double planeNormal[3];
651 for (int iTetraFace = 0; iTetraFace < 4; ++iTetraFace)
653 const int * const tetraFaceNodesConn = tetraFacesNodesConn[iTetraFace];
654 TriangleFaceKey key = TriangleFaceKey(_conn[tetraFaceNodesConn[0]],
655 _conn[tetraFaceNodesConn[1]],
656 _conn[tetraFaceNodesConn[2]]);
657 if (listOfTetraFacesTreated.find(key) == listOfTetraFacesTreated.end())
659 const double * const coordsTetraTriNode1 = _coords + tetraFaceNodesConn[0] * MyMeshType::MY_SPACEDIM;
660 const double * const coordsTetraTriNode2 = _coords + tetraFaceNodesConn[1] * MyMeshType::MY_SPACEDIM;
661 const double * const coordsTetraTriNode3 = _coords + tetraFaceNodesConn[2] * MyMeshType::MY_SPACEDIM;
662 calculateNormalForTria(coordsTetraTriNode1, coordsTetraTriNode2, coordsTetraTriNode3, planeNormal);
663 const double normOfTetraTriNormal = norm(planeNormal);
664 if (epsilonEqual(normOfTetraTriNormal, 0.))
666 for (int i = 0; i < 3; ++i)
673 const double invNormOfTetraTriNormal = 1. / normOfTetraTriNormal;
674 for (int i = 0; i < 3; ++i)
676 planeNormal[i] *= invNormOfTetraTriNormal;
679 double planeConstant = dot(planeNormal, coordsTetraTriNode1);
680 if (IsFacesCoplanar(planeNormal, planeConstant, polyCoords, precision))
682 ConnType nbrPolyTri = polyNodesNbr - 2; // split polygon into nbrPolyTri triangles
683 for (ConnType iTri = 0; iTri < nbrPolyTri; ++iTri)
685 double volume = CalculateIntersectionSurfaceOfCoplanarTriangles(planeNormal,
688 polyCoords[1 + iTri],
689 polyCoords[2 + iTri],
695 if (!epsilonEqual(volume, 0.))
697 totalSurface += volume;
698 listOfTetraFacesColinear.insert(key);
703 listOfTetraFacesTreated.insert(key);
708 // intersect a son with the unit tetra
713 // create the face key
714 TriangleFaceKey key = TriangleFaceKey(polyNodes[0], polyNodes[1], polyNodes[2]);
716 // calculate the triangle if needed
717 if (_volumes.find(key) == _volumes.end())
719 TransformedTriangle tri(_nodes[polyNodes[0]], _nodes[polyNodes[1]], _nodes[polyNodes[2]]);
720 calculateSurface(tri, key);
721 totalSurface += _volumes[key];
725 // count negative as face has reversed orientation
726 totalSurface -= _volumes[key];
733 // simple triangulation of faces along a diagonal :
744 //? not sure if this always works
746 // calculate the triangles if needed
748 // local nodes 1, 2, 3
749 TriangleFaceKey key1 = TriangleFaceKey(polyNodes[0], polyNodes[1], polyNodes[2]);
750 if (_volumes.find(key1) == _volumes.end())
752 TransformedTriangle tri(_nodes[polyNodes[0]], _nodes[polyNodes[1]], _nodes[polyNodes[2]]);
753 calculateSurface(tri, key1);
754 totalSurface += _volumes[key1];
758 // count negative as face has reversed orientation
759 totalSurface -= _volumes[key1];
762 // local nodes 1, 3, 4
763 TriangleFaceKey key2 = TriangleFaceKey(polyNodes[0], polyNodes[2], polyNodes[3]);
764 if (_volumes.find(key2) == _volumes.end())
766 TransformedTriangle tri(_nodes[polyNodes[0]], _nodes[polyNodes[2]], _nodes[polyNodes[3]]);
767 calculateSurface(tri, key2);
768 totalSurface += _volumes[key2];
772 // count negative as face has reversed orientation
773 totalSurface -= _volumes[key2];
780 ConnType nbrPolyTri = polyNodesNbr - 2; // split polygon into nbrPolyTri triangles
781 for (ConnType iTri = 0; iTri < nbrPolyTri; ++iTri)
783 TriangleFaceKey key = TriangleFaceKey(polyNodes[0], polyNodes[1 + iTri], polyNodes[2 + iTri]);
784 if (_volumes.find(key) == _volumes.end())
786 TransformedTriangle tri(_nodes[polyNodes[0]], _nodes[polyNodes[1 + iTri]],
787 _nodes[polyNodes[2 + iTri]]);
788 calculateSurface(tri, key);
789 totalSurface += _volumes[key];
793 totalSurface -= _volumes[key];
801 << "+++ Error : Only elements with triangular and quadratilateral faces are supported at the moment."
809 // reset if it is very small to keep the matrix sparse
810 // is this a good idea?
811 if(epsilonEqual(totalSurface, 0.0, SPARSE_TRUNCATION_LIMIT))
816 LOG(2, "Volume = " << totalSurface << ", det= " << _t->determinant());
822 * Calculates the volume of intersection of this tetrahedron with another one.
824 template<class MyMeshType>
825 double SplitterTetra<MyMeshType>::intersectTetra(const double** tetraCorners)
827 //{ could be done on outside?
828 // check if we have planar tetra element
829 if(_t->determinant() == 0.0)
832 LOG(2, "Planar tetra -- volume 0");
836 const unsigned nbOfNodes4Type=4;
837 // halfspace filtering
838 bool isOutside[8] = {true, true, true, true, true, true, true, true};
839 bool isTargetOutside = false;
841 // calculate the transformed coordinates of the nodes
842 double nodes[nbOfNodes4Type][3];
843 for(int i = 0;i<(int)nbOfNodes4Type;++i)
845 _t->apply(nodes[i], tetraCorners[i]);
846 CheckIsOutside(nodes[i], isOutside);
849 // halfspace filtering check
850 // NB : might not be beneficial for caching of triangles
851 for(int i = 0; i < 8; ++i)
855 isTargetOutside = true;
859 double totalVolume = 0.0;
863 const CellModel& cellModelCell=CellModel::GetCellModel(NORM_TETRA4);
864 ConnType cellNodes[4] = { 0, 1, 2, 3 }, faceNodes[3];
866 for(unsigned ii = 0 ; ii < 4 ; ++ii)
868 cellModelCell.fillSonCellNodalConnectivity(ii,cellNodes,faceNodes);
870 TransformedTriangle tri(nodes[faceNodes[0]], nodes[faceNodes[1]], nodes[faceNodes[2]]);
871 double vol = tri.calculateIntersectionVolume();
872 LOG(1, "ii = " << ii << " Volume=" << vol)
876 // reset if it is very small to keep the matrix sparse
877 // is this a good idea?
878 if(epsilonEqual(totalVolume, 0.0, SPARSE_TRUNCATION_LIMIT))
883 LOG(2, "Volume = " << totalVolume << ", det= " << _t->determinant());
885 // NB : fault in article, Grandy, [8] : it is the determinant of the inverse transformation
886 // that should be used (which is equivalent to dividing by the determinant)
887 return std::fabs(1.0 / _t->determinant() * totalVolume) ;
890 ////////////////////////////////////////////////////////
892 template<class MyMeshTypeT, class MyMeshTypeS>
893 SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::SplitterTetra2(const MyMeshTypeT& targetMesh, const MyMeshTypeS& srcMesh, SplittingPolicy policy)
894 :_target_mesh(targetMesh),_src_mesh(srcMesh),_splitting_pol(policy)
898 template<class MyMeshTypeT, class MyMeshTypeS>
899 SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::~SplitterTetra2()
904 template<class MyMeshTypeT, class MyMeshTypeS>
905 void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::releaseArrays()
907 // free potential sub-mesh nodes that have been allocated
908 if(_nodes.size()>=/*8*/_node_ids.size())
910 typename MyMeshTypeT::MyConnType nbOfNodesT = static_cast<typename MyMeshTypeT::MyConnType>(_node_ids.size());
911 std::vector<const double*>::iterator iter = _nodes.begin() + /*8*/nbOfNodesT;
912 while(iter != _nodes.end())
922 * \param [in] targetCell in C mode.
923 * \param [out] tetra is the output result tetra containers.
925 template<class MyMeshTypeT, class MyMeshTypeS>
926 void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::splitTargetCell2(typename MyMeshTypeT::MyConnType targetCell, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
928 typedef typename MyMeshTypeT::MyConnType TConnType;
929 const TConnType *refConn(_target_mesh.getConnectivityPtr());
930 const TConnType *cellConn(refConn+_target_mesh.getConnectivityIndexPtr()[targetCell]);
931 INTERP_KERNEL::NormalizedCellType gt(_target_mesh.getTypeOfElement(targetCell));
932 std::vector<TConnType> tetrasNodalConn;
933 std::vector<double> addCoords;
934 const double *coords(_target_mesh.getCoordinatesPtr());
935 SplitIntoTetras(_splitting_pol,gt,cellConn,refConn+_target_mesh.getConnectivityIndexPtr()[targetCell+1],coords,tetrasNodalConn,addCoords);
936 std::size_t nbTetras(tetrasNodalConn.size()/4); tetra.resize(nbTetras);
938 typename MyMeshTypeS::MyConnType tmp2[4];
939 for(std::size_t i=0;i<nbTetras;i++)
943 TConnType cellId(tetrasNodalConn[4*i+j]);
947 tmp[j*3+0]=coords[3*cellId+0];
948 tmp[j*3+1]=coords[3*cellId+1];
949 tmp[j*3+2]=coords[3*cellId+2];
953 tmp[j*3+0]=addCoords[3*(-cellId-1)+0];
954 tmp[j*3+1]=addCoords[3*(-cellId-1)+1];
955 tmp[j*3+2]=addCoords[3*(-cellId-1)+2];
958 tetra[i]=new SplitterTetra<MyMeshTypeS>(_src_mesh,tmp,tmp2);
963 * @param targetCell in C mode.
964 * @param tetra is the output result tetra containers.
966 template<class MyMeshTypeT, class MyMeshTypeS>
967 void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::splitTargetCell(typename MyMeshTypeT::MyConnType targetCell,
968 typename MyMeshTypeT::MyConnType nbOfNodesT,
969 typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
971 typedef typename MyMeshTypeT::MyConnType ConnType;
972 const NumberingPolicy numPol=MyMeshTypeT::My_numPol;
973 const int numTetra = static_cast<int>(_splitting_pol);
979 const double *nodes[4];
981 for(int node = 0; node < 4 ; ++node)
983 nodes[node]=getCoordsOfNode2(node, OTT<ConnType,numPol>::indFC(targetCell),_target_mesh,conn[node]);
985 std::copy(conn,conn+4,_node_ids.begin());
986 SplitterTetra<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_src_mesh, nodes,conn);
990 // Issue 0020634. To pass nbOfNodesT to calculateSubNodes (don't want to add an arg)
991 _node_ids.resize(nbOfNodesT);
993 // pre-calculate nodes
994 calculateSubNodes(_target_mesh, OTT<ConnType,numPol>::indFC(targetCell));
996 tetra.reserve(numTetra);
997 _nodes.reserve(30); // we never have more than this
999 switch ( nbOfNodesT )
1003 switch(_splitting_pol)
1007 const int subZone[8] = { 0, 1, 2, 3, 4, 5, 6, 7 };
1008 fiveSplit(subZone,tetra);
1014 const int subZone[8] = { 0, 1, 2, 3, 4, 5, 6, 7 };
1015 sixSplit(subZone,tetra);
1021 calculateGeneral24TetraOld(tetra);
1027 calculateGeneral48Tetra(tetra);
1042 splitConvex(targetCell, tetra);
1048 * Splits the hexahedron into five tetrahedra.
1049 * This method adds five SplitterTetra objects to the vector tetra.
1051 * @param subZone the local node numbers corresponding to the hexahedron corners - these are mapped onto {0,..,7}. Providing this allows the
1052 * splitting to be reused on the subzones of the GENERAL_* types of splitting
1054 template<class MyMeshTypeT, class MyMeshTypeS>
1055 void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::fiveSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
1057 // create tetrahedra
1058 for(int i = 0; i < 5; ++i)
1060 const double* nodes[4];
1061 typename MyMeshTypeS::MyConnType conn[4];
1062 for(int j = 0; j < 4; ++j)
1064 conn[j] = subZone[ SPLIT_NODES_5[4*i+j] ];
1065 typename MyMeshTypeS::MyConnType realConn;
1066 nodes[j] = getCoordsOfSubNode2(conn[j],realConn);
1069 SplitterTetra<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_src_mesh, nodes,conn);
1074 template<class MyMeshTypeT, class MyMeshTypeS>
1075 void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::sixSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
1077 sixSplitGen(subZone,tetra,[](SplitterTetra2& obj, typename MyMeshTypeS::MyConnType& conn, const double *&coords)
1079 typename MyMeshTypeS::MyConnType realConn;
1080 coords = obj.getCoordsOfSubNode2(conn,realConn);
1086 * Splits the hexahedron into six tetrahedra.
1087 * This method adds six SplitterTetra objects to the vector tetra.
1089 * @param subZone the local node numbers corresponding to the hexahedron corners - these are mapped onto {0,..,7}. Providing this allows the
1090 * splitting to be reused on the subzones of the GENERAL_* types of splitting
1092 template<class MyMeshTypeT, class MyMeshTypeS>
1093 void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::sixSplitGen(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra, std::function<void(SplitterTetra2& , typename MyMeshTypeS::MyConnType&, const double*&)> func)
1095 for(int i = 0; i < 6; ++i)
1097 const double* nodes[4];
1098 typename MyMeshTypeS::MyConnType conn[4];
1099 for(int j = 0; j < 4; ++j)
1101 conn[j] = subZone[SPLIT_NODES_6[4*i+j]];
1102 func(*this,conn[j],nodes[j]);
1104 SplitterTetra<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_src_mesh, nodes,conn);
1110 * Version of calculateGeneral24Tetra connectivity aware for P1P0 and P0P1
1112 template<class MyMeshTypeT, class MyMeshTypeS>
1113 void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calculateGeneral24Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
1115 calculateGeneral24TetraGen(tetra,[](SplitterTetra2& obj, typename MyMeshTypeS::MyConnType conn[4], const double* nodes[4]) {
1116 typename MyMeshTypeS::MyConnType realConn;
1117 nodes[2] = obj.getCoordsOfSubNode2(conn[2],realConn); conn[2] = realConn;
1118 nodes[3] = obj.getCoordsOfSubNode2(conn[3],realConn); conn[3] = realConn;
1123 * Version for 3D2DP0P0
1125 template<class MyMeshTypeT, class MyMeshTypeS>
1126 void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calculateGeneral24TetraOld(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
1128 calculateGeneral24TetraGen(tetra,[](SplitterTetra2& obj, typename MyMeshTypeS::MyConnType conn[4], const double* nodes[4]) {
1129 nodes[2] = obj.getCoordsOfSubNode(conn[2]);
1130 nodes[3] = obj.getCoordsOfSubNode(conn[3]);
1135 * Splits the hexahedron into 24 tetrahedra.
1136 * The splitting is done by combining the barycenter of the tetrahedron, the barycenter of each face
1137 * and the nodes of each edge of the face. This creates 6 faces * 4 edges / face = 24 tetrahedra.
1138 * The submesh nodes introduced are the barycenters of the faces and the barycenter of the cell.
1140 template<class MyMeshTypeT, class MyMeshTypeS>
1141 void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calculateGeneral24TetraGen(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra, std::function<void(SplitterTetra2& , typename MyMeshTypeS::MyConnType[4], const double*[4])> func)
1143 // The two nodes of the original mesh cell used in each tetrahedron.
1144 // The tetrahedra all have nodes (cellCenter, faceCenter, edgeNode1, edgeNode2)
1145 // For the correspondence of the nodes, see the GENERAL_48_SUB_NODES table in calculateSubNodes
1147 // nodes to use for tetrahedron
1148 const double* nodes[4];
1149 typename MyMeshTypeS::MyConnType conn[4];
1150 // get the cell center
1152 nodes[0] = getCoordsOfSubNode(conn[0]);
1154 for(int faceCenterNode = 8; faceCenterNode < 14; ++faceCenterNode)
1156 // get the face center
1157 conn[1] = faceCenterNode;
1158 nodes[1] = getCoordsOfSubNode(conn[1]);
1159 for(int j = 0; j < 4; ++j)
1161 const int row = 4*(faceCenterNode - 8) + j;
1162 conn[2] = TETRA_EDGES_GENERAL_24[2*row];
1163 conn[3] = TETRA_EDGES_GENERAL_24[2*row + 1];
1164 func(*this,conn,nodes);
1165 SplitterTetra<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_src_mesh, nodes, conn);
1172 * Splits the hexahedron into 48 tetrahedra.
1173 * The splitting is done by introducing the midpoints of all the edges
1174 * and the barycenter of the element as submesh nodes. The 8 hexahedral subzones thus defined
1175 * are then split into 6 tetrahedra each, as in Grandy, p. 449. The division of the subzones
1176 * is done by calling sixSplit().
1179 template<class MyMeshTypeT, class MyMeshTypeS>
1180 void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calculateGeneral48Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
1182 for(int i = 0; i < 8; ++i)
1184 sixSplitGen(GENERAL_48_SUBZONES+8*i,tetra,[](SplitterTetra2& obj, typename MyMeshTypeS::MyConnType& conn, const double *&coords){
1185 coords = obj.getCoordsOfSubNode(conn);
1191 * Splits the NORM_PYRA5 into 2 tetrahedra.
1193 template<class MyMeshTypeT, class MyMeshTypeS>
1194 void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::splitPyram5(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
1196 static const int SPLIT_PYPA5[2][4] =
1206 // create tetrahedra
1207 const double* nodes[4];
1208 typename MyMeshTypeS::MyConnType conn[4];
1209 for(int i = 0; i < 2; ++i)
1211 for(int j = 0; j < 4; ++j)
1212 nodes[j] = getCoordsOfSubNode2(SPLIT_PYPA5[i][j],conn[j]);
1213 SplitterTetra<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_src_mesh, nodes,conn);
1219 * Splits a convex cell into tetrahedra.
1221 template<class MyMeshTypeT, class MyMeshTypeS>
1222 void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::splitConvex(typename MyMeshTypeT::MyConnType targetCell,
1223 typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
1225 // Each face of a cell is split into triangles and
1226 // each of triangles and a cell barycenter form a tetrahedron.
1228 typedef typename MyMeshTypeT::MyConnType ConnType;
1229 const NumberingPolicy numPol=MyMeshTypeT::My_numPol;
1231 // get type of cell and nb of cell nodes
1232 NormalizedCellType normCellType=_target_mesh.getTypeOfElement(OTT<ConnType,numPol>::indFC(targetCell));
1233 const CellModel& cellModelCell=CellModel::GetCellModel(normCellType);
1234 ConnType nbOfCellNodes=cellModelCell.isDynamic() ? _target_mesh.getNumberOfNodesOfElement(OTT<ConnType,numPol>::indFC(targetCell)) : cellModelCell.getNumberOfNodes();
1236 // get nb of cell sons (faces)
1237 const ConnType* rawCellConn = _target_mesh.getConnectivityPtr() + OTT<ConnType,numPol>::conn2C( _target_mesh.getConnectivityIndexPtr()[ targetCell ]);
1238 const ConnType rawNbCellNodes = _target_mesh.getConnectivityIndexPtr()[ targetCell+1 ] - _target_mesh.getConnectivityIndexPtr()[ targetCell ];
1239 unsigned nbOfSons = cellModelCell.getNumberOfSons2(rawCellConn, rawNbCellNodes);
1241 // indices of nodes of a son
1242 static std::vector<ConnType> allNodeIndices; // == 0,1,2,...,nbOfCellNodes-1
1243 while ( allNodeIndices.size() < (std::size_t)nbOfCellNodes )
1244 allNodeIndices.push_back( static_cast<ConnType>(allNodeIndices.size()) );
1245 std::vector<ConnType> classicFaceNodes(4);
1246 if(cellModelCell.isQuadratic())
1247 throw INTERP_KERNEL::Exception("SplitterTetra2::splitConvex : quadratic 3D cells are not implemented yet !");
1248 ConnType* faceNodes = cellModelCell.isDynamic() ? &allNodeIndices[0] : &classicFaceNodes[0];
1250 // nodes of tetrahedron
1251 typename MyMeshTypeS::MyConnType conn[4];
1252 const double* nodes[4];
1253 nodes[3] = getCoordsOfSubNode2( nbOfCellNodes,conn[3]); // barycenter
1255 for(unsigned ii = 0 ; ii < nbOfSons; ++ii)
1257 // get indices of son's nodes: it's just next portion of allNodeIndices for polyhedron
1258 // and some of allNodeIndices accodring to cell model for a classsic cell
1259 unsigned nbFaceNodes = cellModelCell.getNumberOfNodesConstituentTheSon2(ii, rawCellConn, rawNbCellNodes);
1260 if ( normCellType != NORM_POLYHED )
1261 cellModelCell.fillSonCellNodalConnectivity(ii,&allNodeIndices[0],faceNodes);
1263 int nbTetra = nbFaceNodes - 2; // split polygon into nbTetra triangles
1265 // create tetrahedra
1266 for(int i = 0; i < nbTetra; ++i)
1268 nodes[0] = getCoordsOfSubNode2( faceNodes[0], conn[0]);
1269 nodes[1] = getCoordsOfSubNode2( faceNodes[1+i],conn[1]);
1270 nodes[2] = getCoordsOfSubNode2( faceNodes[2+i],conn[2]);
1271 SplitterTetra<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_src_mesh, nodes,conn);
1275 if ( normCellType == NORM_POLYHED )
1276 faceNodes += nbFaceNodes; // go to the next face
1281 * Precalculates all the nodes.
1282 * Retrieves the mesh nodes and allocates the necessary sub-mesh
1283 * nodes according to the splitting policy used.
1284 * This method is meant to be called once by the constructor.
1286 * @param targetMesh the target mesh
1287 * @param targetCell the global number of the cell that the object represents, in targetMesh mode.
1290 template<class MyMeshTypeT, class MyMeshTypeS>
1291 void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calculateSubNodes(const MyMeshTypeT& targetMesh, typename MyMeshTypeT::MyConnType targetCell)
1293 // retrieve real mesh nodes
1295 typename MyMeshTypeT::MyConnType nbOfNodesT = static_cast<typename MyMeshTypeT::MyConnType>(_node_ids.size());// Issue 0020634. _node_ids.resize(8);
1296 for(int node = 0; node < nbOfNodesT ; ++node)
1298 // calculate only normal nodes
1299 _nodes.push_back(getCoordsOfNode2(node, targetCell, targetMesh,_node_ids[node]));
1302 switch ( nbOfNodesT )
1306 // create sub-mesh nodes if needed
1307 switch(_splitting_pol)
1311 for(int i = 0; i < 7; ++i)
1313 double* barycenter = new double[3];
1314 calcBarycenter(4, barycenter, &GENERAL_24_SUB_NODES[4*i]);
1315 _nodes.push_back(barycenter);
1322 for(int i = 0; i < 19; ++i)
1324 double* barycenter = new double[3];
1325 calcBarycenter(2, barycenter, &GENERAL_48_SUB_NODES[2*i]);
1326 _nodes.push_back(barycenter);
1336 case 5: // NORM_PYRA5
1339 default: // convex 3d cell
1341 // add barycenter of a cell
1342 std::vector<int> allIndices(nbOfNodesT);
1343 for ( int i = 0; i < nbOfNodesT; ++i ) allIndices[i] = i;
1344 double* barycenter = new double[3];
1345 calcBarycenter(nbOfNodesT, barycenter, &allIndices[0]);
1346 _nodes.push_back(barycenter);