_tetra.clear();
}
+ template<class RowType, class ConnType>
+ 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.
template<class MyMeshType, class MyMatrix>
void PolyhedronIntersectorP1P0<MyMeshType,MyMatrix>::intersectCells(ConnType targetCell, const std::vector<ConnType>& srcCells, MyMatrix& res)
{
- SplitterTetra<MyMeshType>* 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<MyMeshType>* subTetras[24];
for(typename std::vector<ConnType>::const_iterator iterCellS=srcCells.begin();iterCellS!=srcCells.end();iterCellS++)
+ {
+ releaseArrays();
+ ConnType nbOfNodesS=this->_src_mesh.getNumberOfNodesOfElement(OTT<ConnType,numPol>::indFC(*iterCellS));
+ _split.splitTargetCell(*iterCellS,nbOfNodesS,_tetra);
+ INTERP_KERNEL::NormalizedCellType srcType = this->_src_mesh.getTypeOfElement( OTT<ConnType,numPol>::indFC(*iterCellS) );
+ if( srcType == NORM_TETRA4 || (srcType == NORM_HEXA8 && sp != GENERAL_24 ))
{
- releaseArrays();
- ConnType nbOfNodesS=Intersector3D<MyMeshType,MyMatrix>::_src_mesh.getNumberOfNodesOfElement(OTT<ConnType,numPol>::indFC(*iterCellS));
- _split.splitTargetCell(*iterCellS,nbOfNodesS,_tetra);
- for(typename std::vector<SplitterTetra<MyMeshType>*>::iterator iter = _tetra.begin(); iter != _tetra.end(); ++iter)
+ for(typename std::vector<SplitterTetra<MyMeshType>*>::const_iterator iter = _tetra.cbegin(); iter != _tetra.cend(); ++iter)
{
(*iter)->splitIntoDualCells(subTetras);
+ double vol2 = 0.;
for(int i=0;i<24;i++)
{
SplitterTetra<MyMeshType> *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<ConnType,numPol>::indFC(sourceNode));
- if(iterRes==resRow.end())
- resRow.insert(std::make_pair(OTT<ConnType,numPol>::indFC(sourceNode),volume));
- else
- {
- double val=(*iterRes).second+volume;
- resRow.erase(OTT<ConnType,numPol>::indFC(sourceNode));
- resRow.insert(std::make_pair(OTT<ConnType,numPol>::indFC(sourceNode),val));
- }
- }
+ AddContributionInRow(resRow,OTT<ConnType,numPol>::indFC(sourceNode),volume);
delete tmp;
}
}
}
+ else
+ {// for HEXA and GENERAL_24 no need to use subsplitting into dual mesh
+ for(typename std::vector<ConnType>::const_iterator iterCellS=srcCells.begin();iterCellS!=srcCells.end();iterCellS++)
+ {
+ releaseArrays();
+ ConnType nbOfNodesS=Intersector3D<MyMeshType,MyMatrix>::_src_mesh.getNumberOfNodesOfElement(OTT<ConnType,numPol>::indFC(*iterCellS));
+ _split.splitTargetCell2(*iterCellS,_tetra);
+ for(typename std::vector<SplitterTetra<MyMeshType>*>::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<ConnType,numPol>::indFC(sourceNode0),volume/2.);
+ AddContributionInRow(resRow,OTT<ConnType,numPol>::indFC(sourceNode1),volume/2.);
+ }
+ }
+ }
+ }
}
}
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<MyMeshTypeS>* >& tetra);
void splitTargetCell(typename MyMeshTypeT::MyConnType targetCell, typename MyMeshTypeT::MyConnType nbOfNodesT,
typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);//to suppress
void fiveSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);//to suppress
void sixSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);//to suppress
void calculateGeneral24Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);//to suppress
+ void calculateGeneral24TetraOld(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
void calculateGeneral48Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);//to suppress
void splitPyram5(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);//to suppress
void splitConvex(typename MyMeshTypeT::MyConnType targetCell,//to suppress
inline const double* getCoordsOfSubNode2(typename MyMeshTypeT::MyConnType node, typename MyMeshTypeT::MyConnType& nodeId);//to suppress
//template<int n>
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<MyMeshTypeS>* >& tetra, std::function<void(SplitterTetra2& , typename MyMeshTypeS::MyConnType&, const double*&)> func);
+ void calculateGeneral24TetraGen(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra, std::function<void(SplitterTetra2& , typename MyMeshTypeS::MyConnType[4], const double*[4])> func);
private:
const MyMeshTypeT& _target_mesh;
const MyMeshTypeS& _src_mesh;
case GENERAL_24:
{
- calculateGeneral24Tetra(tetra);
+ calculateGeneral24TetraOld(tetra);
}
break;
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<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_src_mesh, nodes,conn);
tetra.push_back(t);
}
}
+ template<class MyMeshTypeT, class MyMeshTypeS>
+ void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::sixSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& 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.
* splitting to be reused on the subzones of the GENERAL_* types of splitting
*/
template<class MyMeshTypeT, class MyMeshTypeS>
- void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::sixSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
+ 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)
{
for(int i = 0; i < 6; ++i)
{
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<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_src_mesh, nodes,conn);
tetra.push_back(t);
}
}
+ /**
+ * Version of calculateGeneral24Tetra connectivity aware for P1P0 and P0P1
+ */
+ template<class MyMeshTypeT, class MyMeshTypeS>
+ void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calculateGeneral24Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& 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<class MyMeshTypeT, class MyMeshTypeS>
+ void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calculateGeneral24TetraOld(typename std::vector< SplitterTetra<MyMeshTypeS>* >& 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<class MyMeshTypeT, class MyMeshTypeS>
- void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calculateGeneral24Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
+ void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calculateGeneral24TetraGen(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra, std::function<void(SplitterTetra2& , typename MyMeshTypeS::MyConnType[4], const double*[4])> func)
{
// The two nodes of the original mesh cell used in each tetrahedron.
// The tetrahedra all have nodes (cellCenter, faceCenter, edgeNode1, edgeNode2)
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)
{
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<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_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
{
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);
+ });
}
}
INTERP_KERNEL::Interpolation3D myInterpolator;
std::vector<std::map<mcIdType,double> > 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();
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))
