-// Copyright (C) 2007-2014 CEA/DEN, EDF R&D
+// Copyright (C) 2007-2024 CEA, EDF
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
* \param [in] tetraCorners array 4*3 doubles containing corners of input tetrahedron (P0X,P0Y,P0Y,P1X,P1Y,P1Z,P2X,P2Y,P2Z,P3X,P3Y,P3Z).
*/
template<class MyMeshType>
- SplitterTetra<MyMeshType>::SplitterTetra(const MyMeshType& srcMesh, const double tetraCorners[12], const int *conn): _t(0),_src_mesh(srcMesh)
+ SplitterTetra<MyMeshType>::SplitterTetra(const MyMeshType& srcMesh, const double tetraCorners[12], const ConnType *conn): _t(0),_src_mesh(srcMesh)
{
if(!conn)
{ _conn[0]=0; _conn[1]=1; _conn[2]=2; _conn[3]=3; }
SplitterTetra<MyMeshType>::~SplitterTetra()
{
delete _t;
- for(HashMap< int, double* >::iterator iter = _nodes.begin(); iter != _nodes.end() ; ++iter)
+ for(typename HashMap< ConnType, double* >::iterator iter = _nodes.begin(); iter != _nodes.end() ; ++iter)
delete[] iter->second;
}
// get type of cell
NormalizedCellType normCellType=_src_mesh.getTypeOfElement(OTT<ConnType,numPol>::indFC(element));
const CellModel& cellModelCell=CellModel::GetCellModel(normCellType);
- unsigned nbOfNodes4Type=cellModelCell.isDynamic() ? _src_mesh.getNumberOfNodesOfElement(OTT<ConnType,numPol>::indFC(element)) : cellModelCell.getNumberOfNodes();
+ ConnType nbOfNodes4Type=cellModelCell.isDynamic() ? _src_mesh.getNumberOfNodesOfElement(OTT<ConnType,numPol>::indFC(element)) : cellModelCell.getNumberOfNodes();
// halfspace filtering
bool isOutside[8] = {true, true, true, true, true, true, true, true};
bool isTargetOutside = false;
// calculate the coordinates of the nodes
- int *cellNodes=new int[nbOfNodes4Type];
- for(int i = 0;i<(int)nbOfNodes4Type;++i)
+ ConnType *cellNodes=new ConnType[nbOfNodes4Type];
+ for(ConnType i = 0;i<nbOfNodes4Type;++i)
{
// we could store mapping local -> global numbers too, but not sure it is worth it
- const int globalNodeNum = getGlobalNumberOfNode(i, OTT<ConnType,numPol>::indFC(element), _src_mesh);
+ const ConnType globalNodeNum = getGlobalNumberOfNode(i, OTT<ConnType,numPol>::indFC(element), _src_mesh);
cellNodes[i]=globalNodeNum;
if(_nodes.find(globalNodeNum) == _nodes.end())
{
// halfspace filtering check
// NB : might not be beneficial for caching of triangles
for(int i = 0; i < 8; ++i)
- {
- if(isOutside[i])
- {
- isTargetOutside = true;
- }
- }
+ if(isOutside[i])
+ isTargetOutside = true;
double totalVolume = 0.0;
// get nb of sons of a cell
const ConnType* rawCellConn = _src_mesh.getConnectivityPtr() + OTT<ConnType,numPol>::conn2C( _src_mesh.getConnectivityIndexPtr()[ element ]);
- const int rawNbCellNodes = _src_mesh.getConnectivityIndexPtr()[ element+1 ] - _src_mesh.getConnectivityIndexPtr()[ element ];
+ const ConnType rawNbCellNodes = _src_mesh.getConnectivityIndexPtr()[ element+1 ] - _src_mesh.getConnectivityIndexPtr()[ element ];
unsigned nbOfSons = cellModelCell.getNumberOfSons2(rawCellConn, rawNbCellNodes);
for(unsigned ii = 0 ; ii < nbOfSons; ++ii)
{
// get sons connectivity
NormalizedCellType faceType;
- int *faceNodes, nbFaceNodes=-1;
+ ConnType *faceNodes, nbFaceNodes=-1;
if ( cellModelCell.isDynamic() )
{
- faceNodes=new int[nbOfNodes4Type];
+ faceNodes=new ConnType[nbOfNodes4Type];
nbFaceNodes = cellModelCell.fillSonCellNodalConnectivity2(ii,rawCellConn,rawNbCellNodes,faceNodes,faceType);
- for ( int i = 0; i < nbFaceNodes; ++i )
+ for ( ConnType i = 0; i < nbFaceNodes; ++i )
faceNodes[i] = OTT<ConnType,numPol>::coo2C(faceNodes[i]);
}
else
{
faceType = cellModelCell.getSonType(ii);
- const CellModel& faceModel=CellModel::GetCellModel(faceType);
- assert(faceModel.getDimension() == 2);
- faceNodes=new int[faceModel.getNumberOfNodes()];
+ assert(CellModel::GetCellModel(faceType).getDimension() == 2);
+ nbFaceNodes = cellModelCell.getNumberOfNodesConstituentTheSon(ii);
+ faceNodes = new ConnType[nbFaceNodes];
cellModelCell.fillSonCellNodalConnectivity(ii,cellNodes,faceNodes);
}
// intersect a son with the unit tetra
case NORM_POLYGON:
{
- int nbTria = nbFaceNodes - 2; // split polygon into nbTria triangles
- for ( int iTri = 0; iTri < nbTria; ++iTri )
+ ConnType nbTria = nbFaceNodes - 2; // split polygon into nbTria triangles
+ for ( ConnType iTri = 0; iTri < nbTria; ++iTri )
{
TriangleFaceKey key = TriangleFaceKey(faceNodes[0], faceNodes[1+iTri], faceNodes[2+iTri]);
if(_volumes.find(key) == _volumes.end())
/**
* Calculates the intersection surface of two coplanar triangles.
*
- * @param palneNormal normal of the plane for the first triangle
+ * @param planeNormal normal of the plane for the first triangle
* @param planeConstant constant of the equation of the plane for the first triangle
* @param p1 coordinates of the first node of the first triangle
* @param p2 coordinates of the second node of the first triangle
{
typedef typename MyMeshType::MyConnType ConnType;
typedef double Vect2[2];
- typedef double Vect3[3];
typedef double Triangle2[3][2];
const double *const tri0[3] = {p1, p2, p3};
* Determine if a face is coplanar with a triangle.
* The first face is characterized by the equation of her plane
*
- * @param palneNormal normal of the plane for the first triangle
+ * @param planeNormal normal of the plane for the first triangle
* @param planeConstant constant of the equation of the plane for the first triangle
* @param coordsFace coordinates of the triangle face
* @param precision precision for double float data used for comparison
*/
template<class MyMeshType>
double SplitterTetra<MyMeshType>::intersectSourceFace(const NormalizedCellType polyType,
- const int polyNodesNbr,
- const int *const polyNodes,
+ const ConnType polyNodesNbr,
+ const ConnType *const polyNodes,
const double *const *const polyCoords,
const double dimCaracteristic,
const double precision,
std::multiset<TriangleFaceKey>& listOfTetraFacesTreated,
std::set<TriangleFaceKey>& listOfTetraFacesColinear)
{
- typedef typename MyMeshType::MyConnType ConnType;
-
double totalSurface = 0.0;
// check if we have planar tetra element
bool isTargetOutside = false;
// calculate the coordinates of the nodes
- for(int i = 0;i<(int)polyNodesNbr;++i)
+ for(ConnType i = 0;i<polyNodesNbr;++i)
{
- const int globalNodeNum = polyNodes[i];
+ const ConnType globalNodeNum = polyNodes[i];
if(_nodes.find(globalNodeNum) == _nodes.end())
{
calculateNode2(globalNodeNum, polyCoords[i]);
double planeConstant = dot(planeNormal, coordsTetraTriNode1);
if (IsFacesCoplanar(planeNormal, planeConstant, polyCoords, precision))
{
- int nbrPolyTri = polyNodesNbr - 2; // split polygon into nbrPolyTri triangles
- for (int iTri = 0; iTri < nbrPolyTri; ++iTri)
+ ConnType nbrPolyTri = polyNodesNbr - 2; // split polygon into nbrPolyTri triangles
+ for (ConnType iTri = 0; iTri < nbrPolyTri; ++iTri)
{
double volume = CalculateIntersectionSurfaceOfCoplanarTriangles(planeNormal,
planeConstant,
case NORM_POLYGON:
{
- int nbrPolyTri = polyNodesNbr - 2; // split polygon into nbrPolyTri triangles
- for (int iTri = 0; iTri < nbrPolyTri; ++iTri)
+ ConnType nbrPolyTri = polyNodesNbr - 2; // split polygon into nbrPolyTri triangles
+ for (ConnType iTri = 0; iTri < nbrPolyTri; ++iTri)
{
TriangleFaceKey key = TriangleFaceKey(polyNodes[0], polyNodes[1 + iTri], polyNodes[2 + iTri]);
if (_volumes.find(key) == _volumes.end())
if(!isTargetOutside)
{
const CellModel& cellModelCell=CellModel::GetCellModel(NORM_TETRA4);
- int cellNodes[4] = { 0, 1, 2, 3 }, faceNodes[3];
+ ConnType cellNodes[4] = { 0, 1, 2, 3 }, faceNodes[3];
for(unsigned ii = 0 ; ii < 4 ; ++ii)
{
TransformedTriangle tri(nodes[faceNodes[0]], nodes[faceNodes[1]], nodes[faceNodes[2]]);
double vol = tri.calculateIntersectionVolume();
+ LOG(1, "ii = " << ii << " Volume=" << vol)
totalVolume += vol;
}
void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::releaseArrays()
{
// free potential sub-mesh nodes that have been allocated
- typename MyMeshTypeT::MyConnType nbOfNodesT = _node_ids.size();// Issue 0020634.
- if((int)_nodes.size()>=/*8*/nbOfNodesT)
+ if(_nodes.size()>=/*8*/_node_ids.size())
{
+ typename MyMeshTypeT::MyConnType nbOfNodesT = static_cast<typename MyMeshTypeT::MyConnType>(_node_ids.size());
std::vector<const double*>::iterator iter = _nodes.begin() + /*8*/nbOfNodesT;
while(iter != _nodes.end())
{
template<class MyMeshTypeT, class MyMeshTypeS>
void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::splitTargetCell2(typename MyMeshTypeT::MyConnType targetCell, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra)
{
- const int *refConn(_target_mesh.getConnectivityPtr());
- const int *cellConn(refConn+_target_mesh.getConnectivityIndexPtr()[targetCell]);
+ typedef typename MyMeshTypeT::MyConnType TConnType;
+ const TConnType *refConn(_target_mesh.getConnectivityPtr());
+ const TConnType *cellConn(refConn+_target_mesh.getConnectivityIndexPtr()[targetCell]);
INTERP_KERNEL::NormalizedCellType gt(_target_mesh.getTypeOfElement(targetCell));
- std::vector<int> tetrasNodalConn;
+ std::vector<TConnType> tetrasNodalConn;
std::vector<double> addCoords;
const double *coords(_target_mesh.getCoordinatesPtr());
SplitIntoTetras(_splitting_pol,gt,cellConn,refConn+_target_mesh.getConnectivityIndexPtr()[targetCell+1],coords,tetrasNodalConn,addCoords);
std::size_t nbTetras(tetrasNodalConn.size()/4); tetra.resize(nbTetras);
double tmp[12];
- int tmp2[4];
+ typename MyMeshTypeS::MyConnType tmp2[4];
for(std::size_t i=0;i<nbTetras;i++)
{
for(int j=0;j<4;j++)
{
- int cellId(tetrasNodalConn[4*i+j]);
+ TConnType cellId(tetrasNodalConn[4*i+j]);
tmp2[j]=cellId;
if(cellId>=0)
{
_node_ids.resize(8);
tetra.reserve(1);
const double *nodes[4];
- int conn[4];
+ ConnType conn[4];
for(int node = 0; node < 4 ; ++node)
{
nodes[node]=getCoordsOfNode2(node, OTT<ConnType,numPol>::indFC(targetCell),_target_mesh,conn[node]);
case GENERAL_24:
{
- calculateGeneral24Tetra(tetra);
+ calculateGeneral24TetraOld(tetra);
}
break;
for(int i = 0; i < 5; ++i)
{
const double* nodes[4];
- int conn[4];
+ typename MyMeshTypeS::MyConnType conn[4];
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)
{
const double* nodes[4];
- int conn[4];
+ typename MyMeshTypeS::MyConnType conn[4];
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)
- // For the correspondance of the nodes, see the GENERAL_48_SUB_NODES table in calculateSubNodes
+ // For the correspondence of the nodes, see the GENERAL_48_SUB_NODES table in calculateSubNodes
// nodes to use for tetrahedron
const double* nodes[4];
- int conn[4];
+ 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);
+ });
}
}
// create tetrahedra
const double* nodes[4];
- int conn[4];
+ typename MyMeshTypeS::MyConnType conn[4];
for(int i = 0; i < 2; ++i)
{
for(int j = 0; j < 4; ++j)
// get type of cell and nb of cell nodes
NormalizedCellType normCellType=_target_mesh.getTypeOfElement(OTT<ConnType,numPol>::indFC(targetCell));
const CellModel& cellModelCell=CellModel::GetCellModel(normCellType);
- unsigned nbOfCellNodes=cellModelCell.isDynamic() ? _target_mesh.getNumberOfNodesOfElement(OTT<ConnType,numPol>::indFC(targetCell)) : cellModelCell.getNumberOfNodes();
+ ConnType nbOfCellNodes=cellModelCell.isDynamic() ? _target_mesh.getNumberOfNodesOfElement(OTT<ConnType,numPol>::indFC(targetCell)) : cellModelCell.getNumberOfNodes();
// get nb of cell sons (faces)
const ConnType* rawCellConn = _target_mesh.getConnectivityPtr() + OTT<ConnType,numPol>::conn2C( _target_mesh.getConnectivityIndexPtr()[ targetCell ]);
- const int rawNbCellNodes = _target_mesh.getConnectivityIndexPtr()[ targetCell+1 ] - _target_mesh.getConnectivityIndexPtr()[ targetCell ];
+ const ConnType rawNbCellNodes = _target_mesh.getConnectivityIndexPtr()[ targetCell+1 ] - _target_mesh.getConnectivityIndexPtr()[ targetCell ];
unsigned nbOfSons = cellModelCell.getNumberOfSons2(rawCellConn, rawNbCellNodes);
// indices of nodes of a son
- static std::vector<int> allNodeIndices; // == 0,1,2,...,nbOfCellNodes-1
- while ( allNodeIndices.size() < nbOfCellNodes )
- allNodeIndices.push_back( allNodeIndices.size() );
- std::vector<int> classicFaceNodes(4);
- int* faceNodes = cellModelCell.isDynamic() ? &allNodeIndices[0] : &classicFaceNodes[0];
+ static std::vector<ConnType> allNodeIndices; // == 0,1,2,...,nbOfCellNodes-1
+ while ( allNodeIndices.size() < (std::size_t)nbOfCellNodes )
+ allNodeIndices.push_back( static_cast<ConnType>(allNodeIndices.size()) );
+ std::vector<ConnType> classicFaceNodes(4);
+ if(cellModelCell.isQuadratic())
+ throw INTERP_KERNEL::Exception("SplitterTetra2::splitConvex : quadratic 3D cells are not implemented yet !");
+ ConnType* faceNodes = cellModelCell.isDynamic() ? &allNodeIndices[0] : &classicFaceNodes[0];
// nodes of tetrahedron
- int conn[4];
+ typename MyMeshTypeS::MyConnType conn[4];
const double* nodes[4];
nodes[3] = getCoordsOfSubNode2( nbOfCellNodes,conn[3]); // barycenter
*
* @param targetMesh the target mesh
* @param targetCell the global number of the cell that the object represents, in targetMesh mode.
- * @param policy the splitting policy of the object
*
*/
template<class MyMeshTypeT, class MyMeshTypeS>
{
// retrieve real mesh nodes
- typename MyMeshTypeT::MyConnType nbOfNodesT = _node_ids.size();// Issue 0020634. _node_ids.resize(8);
+ typename MyMeshTypeT::MyConnType nbOfNodesT = static_cast<typename MyMeshTypeT::MyConnType>(_node_ids.size());// Issue 0020634. _node_ids.resize(8);
for(int node = 0; node < nbOfNodesT ; ++node)
{
// calculate only normal nodes
default:
break;
}
+ break;
case 5: // NORM_PYRA5
break;