From 009cf624f2bb7a09c1031af917bb9086fe520853 Mon Sep 17 00:00:00 2001 From: Anthony Geay Date: Mon, 7 Aug 2023 16:29:10 +0200 Subject: [PATCH] [EDF27859] : Correct bug in case of HEXA/HEXA in P1P0 mode with PLANAR_FACE5 / PLANAR_FACE6 / GENERAL_24 --- .../PolyhedronIntersectorP1P0.txx | 66 ++++++++++++++----- src/INTERP_KERNEL/SplitterTetra.hxx | 5 ++ src/INTERP_KERNEL/SplitterTetra.txx | 60 +++++++++++++---- .../Test/MEDCouplingBasicsTestInterp.cxx | 4 +- .../MEDCouplingRemapperTest.py | 29 +++++++- 5 files changed, 132 insertions(+), 32 deletions(-) diff --git a/src/INTERP_KERNEL/PolyhedronIntersectorP1P0.txx b/src/INTERP_KERNEL/PolyhedronIntersectorP1P0.txx index 5741db4ac..844c3982d 100644 --- a/src/INTERP_KERNEL/PolyhedronIntersectorP1P0.txx +++ b/src/INTERP_KERNEL/PolyhedronIntersectorP1P0.txx @@ -66,6 +66,23 @@ namespace INTERP_KERNEL _tetra.clear(); } + template + void AddContributionInRow(RowType& row, ConnType colId, double value) + { + if(value != 0.) + { + typename RowType::const_iterator iterRes=row.find(colId); + if(iterRes==row.end()) + row.insert(std::make_pair(colId,value)); + else + { + double val=(*iterRes).second+value; + row.erase(colId); + row.insert(std::make_pair(colId,val)); + } + } + } + /** * Calculates the volume of intersection of an element in the source mesh and the target element * represented by the object. @@ -80,37 +97,52 @@ namespace INTERP_KERNEL template void PolyhedronIntersectorP1P0::intersectCells(ConnType targetCell, const std::vector& srcCells, MyMatrix& res) { - SplitterTetra* subTetras[24]; typename MyMatrix::value_type& resRow=res[targetCell]; + INTERP_KERNEL::SplittingPolicy sp( _split.getSplittingPolicy() ); + if( sp == GENERAL_48 ) + THROW_IK_EXCEPTION("GENERAL_28 spliting is not supported for P1P0 interpolation"); + SplitterTetra* subTetras[24]; for(typename std::vector::const_iterator iterCellS=srcCells.begin();iterCellS!=srcCells.end();iterCellS++) + { + releaseArrays(); + ConnType nbOfNodesS=this->_src_mesh.getNumberOfNodesOfElement(OTT::indFC(*iterCellS)); + _split.splitTargetCell(*iterCellS,nbOfNodesS,_tetra); + INTERP_KERNEL::NormalizedCellType srcType = this->_src_mesh.getTypeOfElement( OTT::indFC(*iterCellS) ); + if( srcType == NORM_TETRA4 || (srcType == NORM_HEXA8 && sp != GENERAL_24 )) { - releaseArrays(); - ConnType nbOfNodesS=Intersector3D::_src_mesh.getNumberOfNodesOfElement(OTT::indFC(*iterCellS)); - _split.splitTargetCell(*iterCellS,nbOfNodesS,_tetra); - for(typename std::vector*>::iterator iter = _tetra.begin(); iter != _tetra.end(); ++iter) + for(typename std::vector*>::const_iterator iter = _tetra.cbegin(); iter != _tetra.cend(); ++iter) { (*iter)->splitIntoDualCells(subTetras); + double vol2 = 0.; for(int i=0;i<24;i++) { SplitterTetra *tmp=subTetras[i]; double volume = tmp->intersectSourceCell(targetCell); + vol2 += volume; ConnType sourceNode=tmp->getId(0); - if(volume!=0.) - { - typename MyMatrix::value_type::const_iterator iterRes=resRow.find(OTT::indFC(sourceNode)); - if(iterRes==resRow.end()) - resRow.insert(std::make_pair(OTT::indFC(sourceNode),volume)); - else - { - double val=(*iterRes).second+volume; - resRow.erase(OTT::indFC(sourceNode)); - resRow.insert(std::make_pair(OTT::indFC(sourceNode),val)); - } - } + AddContributionInRow(resRow,OTT::indFC(sourceNode),volume); delete tmp; } } } + else + {// for HEXA and GENERAL_24 no need to use subsplitting into dual mesh + for(typename std::vector::const_iterator iterCellS=srcCells.begin();iterCellS!=srcCells.end();iterCellS++) + { + releaseArrays(); + ConnType nbOfNodesS=Intersector3D::_src_mesh.getNumberOfNodesOfElement(OTT::indFC(*iterCellS)); + _split.splitTargetCell2(*iterCellS,_tetra); + for(typename std::vector*>::const_iterator iter = _tetra.cbegin(); iter != _tetra.cend(); ++iter) + { + double volume = std::abs( (*iter)->intersectSourceCell(targetCell) ); + // node #0 is for internal node node #1 is for the node at the middle of the face + ConnType sourceNode0( (*iter)->getId(0) ), sourceNode1( (*iter)->getId(1) ); + AddContributionInRow(resRow,OTT::indFC(sourceNode0),volume/2.); + AddContributionInRow(resRow,OTT::indFC(sourceNode1),volume/2.); + } + } + } + } } } diff --git a/src/INTERP_KERNEL/SplitterTetra.hxx b/src/INTERP_KERNEL/SplitterTetra.hxx index 30c46ddcb..5a5e3a559 100644 --- a/src/INTERP_KERNEL/SplitterTetra.hxx +++ b/src/INTERP_KERNEL/SplitterTetra.hxx @@ -544,12 +544,14 @@ namespace INTERP_KERNEL SplitterTetra2(const MyMeshTypeT& targetMesh, const MyMeshTypeS& srcMesh, SplittingPolicy policy); ~SplitterTetra2(); void releaseArrays(); + SplittingPolicy getSplittingPolicy() const { return _splitting_pol; } void splitTargetCell2(typename MyMeshTypeT::MyConnType targetCell, typename std::vector< SplitterTetra* >& tetra); void splitTargetCell(typename MyMeshTypeT::MyConnType targetCell, typename MyMeshTypeT::MyConnType nbOfNodesT, typename std::vector< SplitterTetra* >& tetra);//to suppress void fiveSplit(const int* const subZone, typename std::vector< SplitterTetra* >& tetra);//to suppress void sixSplit(const int* const subZone, typename std::vector< SplitterTetra* >& tetra);//to suppress void calculateGeneral24Tetra(typename std::vector< SplitterTetra* >& tetra);//to suppress + void calculateGeneral24TetraOld(typename std::vector< SplitterTetra* >& tetra); void calculateGeneral48Tetra(typename std::vector< SplitterTetra* >& tetra);//to suppress void splitPyram5(typename std::vector< SplitterTetra* >& tetra);//to suppress void splitConvex(typename MyMeshTypeT::MyConnType targetCell,//to suppress @@ -559,6 +561,9 @@ namespace INTERP_KERNEL inline const double* getCoordsOfSubNode2(typename MyMeshTypeT::MyConnType node, typename MyMeshTypeT::MyConnType& nodeId);//to suppress //template inline void calcBarycenter(typename MyMeshTypeT::MyConnType n, double* barycenter, const int* pts);//to suppress + private: + void sixSplitGen(const int* const subZone, typename std::vector< SplitterTetra* >& tetra, std::function func); + void calculateGeneral24TetraGen(typename std::vector< SplitterTetra* >& tetra, std::function func); private: const MyMeshTypeT& _target_mesh; const MyMeshTypeS& _src_mesh; diff --git a/src/INTERP_KERNEL/SplitterTetra.txx b/src/INTERP_KERNEL/SplitterTetra.txx index 77b54d2cc..94b5b1624 100644 --- a/src/INTERP_KERNEL/SplitterTetra.txx +++ b/src/INTERP_KERNEL/SplitterTetra.txx @@ -1021,7 +1021,7 @@ namespace INTERP_KERNEL case GENERAL_24: { - calculateGeneral24Tetra(tetra); + calculateGeneral24TetraOld(tetra); } break; @@ -1065,13 +1065,26 @@ namespace INTERP_KERNEL for(int j = 0; j < 4; ++j) { conn[j] = subZone[ SPLIT_NODES_5[4*i+j] ]; - nodes[j] = getCoordsOfSubNode(conn[j]); + typename MyMeshTypeS::MyConnType realConn; + nodes[j] = getCoordsOfSubNode2(conn[j],realConn); + conn[j] = realConn; } SplitterTetra* t = new SplitterTetra(_src_mesh, nodes,conn); tetra.push_back(t); } } + template + void SplitterTetra2::sixSplit(const int* const subZone, typename std::vector< SplitterTetra* >& tetra) + { + sixSplitGen(subZone,tetra,[](SplitterTetra2& obj, typename MyMeshTypeS::MyConnType& conn, const double *&coords) + { + typename MyMeshTypeS::MyConnType realConn; + coords = obj.getCoordsOfSubNode2(conn,realConn); + conn = realConn; + }); + } + /** * Splits the hexahedron into six tetrahedra. * This method adds six SplitterTetra objects to the vector tetra. @@ -1080,7 +1093,7 @@ namespace INTERP_KERNEL * splitting to be reused on the subzones of the GENERAL_* types of splitting */ template - void SplitterTetra2::sixSplit(const int* const subZone, typename std::vector< SplitterTetra* >& tetra) + void SplitterTetra2::sixSplitGen(const int* const subZone, typename std::vector< SplitterTetra* >& tetra, std::function func) { for(int i = 0; i < 6; ++i) { @@ -1089,22 +1102,46 @@ namespace INTERP_KERNEL for(int j = 0; j < 4; ++j) { conn[j] = subZone[SPLIT_NODES_6[4*i+j]]; - nodes[j] = getCoordsOfSubNode(conn[j]); + func(*this,conn[j],nodes[j]); } SplitterTetra* t = new SplitterTetra(_src_mesh, nodes,conn); tetra.push_back(t); } } + /** + * Version of calculateGeneral24Tetra connectivity aware for P1P0 and P0P1 + */ + template + void SplitterTetra2::calculateGeneral24Tetra(typename std::vector< SplitterTetra* >& tetra) + { + calculateGeneral24TetraGen(tetra,[](SplitterTetra2& obj, typename MyMeshTypeS::MyConnType conn[4], const double* nodes[4]) { + typename MyMeshTypeS::MyConnType realConn; + nodes[2] = obj.getCoordsOfSubNode2(conn[2],realConn); conn[2] = realConn; + nodes[3] = obj.getCoordsOfSubNode2(conn[3],realConn); conn[3] = realConn; + }); + } + + /*! + * Version for 3D2DP0P0 + */ + template + void SplitterTetra2::calculateGeneral24TetraOld(typename std::vector< SplitterTetra* >& tetra) + { + calculateGeneral24TetraGen(tetra,[](SplitterTetra2& obj, typename MyMeshTypeS::MyConnType conn[4], const double* nodes[4]) { + nodes[2] = obj.getCoordsOfSubNode(conn[2]); + nodes[3] = obj.getCoordsOfSubNode(conn[3]); + }); + } + /** * Splits the hexahedron into 24 tetrahedra. * The splitting is done by combining the barycenter of the tetrahedron, the barycenter of each face * and the nodes of each edge of the face. This creates 6 faces * 4 edges / face = 24 tetrahedra. * The submesh nodes introduced are the barycenters of the faces and the barycenter of the cell. - * */ template - void SplitterTetra2::calculateGeneral24Tetra(typename std::vector< SplitterTetra* >& tetra) + void SplitterTetra2::calculateGeneral24TetraGen(typename std::vector< SplitterTetra* >& tetra, std::function func) { // The two nodes of the original mesh cell used in each tetrahedron. // The tetrahedra all have nodes (cellCenter, faceCenter, edgeNode1, edgeNode2) @@ -1115,7 +1152,7 @@ namespace INTERP_KERNEL typename MyMeshTypeS::MyConnType conn[4]; // get the cell center conn[0] = 14; - nodes[0] = getCoordsOfSubNode(conn[0]); + nodes[0] = getCoordsOfSubNode(conn[0]); for(int faceCenterNode = 8; faceCenterNode < 14; ++faceCenterNode) { @@ -1127,16 +1164,13 @@ namespace INTERP_KERNEL const int row = 4*(faceCenterNode - 8) + j; conn[2] = TETRA_EDGES_GENERAL_24[2*row]; conn[3] = TETRA_EDGES_GENERAL_24[2*row + 1]; - nodes[2] = getCoordsOfSubNode(conn[2]); - nodes[3] = getCoordsOfSubNode(conn[3]); - + func(*this,conn,nodes); SplitterTetra* t = new SplitterTetra(_src_mesh, nodes, conn); tetra.push_back(t); } } } - /** * Splits the hexahedron into 48 tetrahedra. * The splitting is done by introducing the midpoints of all the edges @@ -1150,7 +1184,9 @@ namespace INTERP_KERNEL { for(int i = 0; i < 8; ++i) { - sixSplit(GENERAL_48_SUBZONES+8*i,tetra); + sixSplitGen(GENERAL_48_SUBZONES+8*i,tetra,[](SplitterTetra2& obj, typename MyMeshTypeS::MyConnType& conn, const double *&coords){ + coords = obj.getCoordsOfSubNode(conn); + }); } } diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTestInterp.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTestInterp.cxx index 1977fabd8..a188997b6 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTestInterp.cxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTestInterp.cxx @@ -1287,8 +1287,8 @@ void MEDCouplingBasicsTestInterp::test3DInterpP1P0_1() INTERP_KERNEL::Interpolation3D myInterpolator; std::vector > res; myInterpolator.setPrecision(1e-12); - INTERP_KERNEL::SplittingPolicy sp[] = { INTERP_KERNEL::PLANAR_FACE_5, INTERP_KERNEL::PLANAR_FACE_6, INTERP_KERNEL::GENERAL_24, INTERP_KERNEL::GENERAL_48 }; - for ( int i = 0; i < 4; ++i ) + INTERP_KERNEL::SplittingPolicy sp[] = { INTERP_KERNEL::PLANAR_FACE_5, INTERP_KERNEL::PLANAR_FACE_6, INTERP_KERNEL::GENERAL_24}; + for ( int i = 0; i < 3; ++i ) { myInterpolator.setSplittingPolicy( sp[i] ); res.clear(); diff --git a/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py b/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py index 4f1551dba..c5221bbaa 100644 --- a/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py @@ -1651,7 +1651,34 @@ class MEDCouplingBasicsTest(unittest.TestCase): except InterpKernelException as e: self.assertTrue("fail to locate point" in e.what()) else: - self.assertTrue(false,"") + self.assertTrue(false,"") + + def testP1P0OnHexa_1(self): + """ + See EDF27859 : This test focuses on P1P0 interpolation with source containing HEXA. So P1P0 intersector is going to split into tetras + the source cell. + """ + trgMesh = MEDCouplingUMesh("mesh",3) + trgMesh.setCoords( DataArrayDouble([18500.0, 0.0, 0.0, 18544.0, 0.0, 0.0, 18544.0, 0.0, 200.0, 18500.0, 0.0, 200.0, 18497.96424104365, 274.44295777156043, 0.0, 18541.959399238567, 275.0956869684225, 0.0, 18541.959399238567, 275.0956869684225, 200.0, 18497.96424104365, 274.44295777156043, 200.0],8,3) ) + firstPts = DataArrayDouble(3*10) ; firstPts[:] = 0. ; firstPts.rearrange(3) + trgMesh.setCoords( DataArrayDouble.Aggregate(firstPts,trgMesh.getCoords()) ) # this line is important to check that correct ids are taken into account + trgMesh.allocateCells(1) + trgMesh.insertNextCell(NORM_HEXA8,[10,11,12,13,14,15,16,17]) + trgMesh.writeVTK("trgMeshNonReg.vtu") + + srcMesh = trgMesh.deepCopy() + cc = trgMesh.computeCellCenterOfMass()[0] + trgMesh.scale(cc,1.01) # This line is to workaround the EDF28414 bug inside 3D intersector + + expectedMatrix0 = [{10: 503624.09065889206, 11: 100868.41855508549, 12: 503863.42469772906, 13: 100629.0845162416, 14: 100629.08451623631, 15: 503863.4246977626, 16: 100868.418555101, 17: 503624.0906588909}] + expectedMatrix1 = [{10: 604492.509213978, 11: 201736.8371101737, 12: 201736.83711016813, 13: 201497.50307132734, 14: 201258.16903247262, 15: 201497.50307133005, 16: 604492.5092140044, 17: 201258.16903247265}] + expectedMatrix2 = [{10: 302066.754077835, 11: 302425.7551361466, 12: 302425.7551361466, 13: 302066.754077835, 14: 302066.7540778395, 15: 302425.7551361595, 16: 302425.7551361595, 17: 302066.75407783955}] + for sp,expectedMatrix in [ (PLANAR_FACE_5,expectedMatrix0),(PLANAR_FACE_6,expectedMatrix1),(GENERAL_24,expectedMatrix2)]: + remap = MEDCouplingRemapper() + remap.setSplittingPolicy( sp ) + remap.prepare(srcMesh,trgMesh,"P1P0") + mat = remap.getCrudeMatrix() + self.checkMatrix(expectedMatrix,mat,18,1.0) def checkMatrix(self,mat1,mat2,nbCols,eps): self.assertEqual(len(mat1),len(mat2)) -- 2.39.2