_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];
- 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.splitTargetCell(*iterCellS,nbOfNodesS,_tetra);
- for(typename std::vector<SplitterTetra<MyMeshType>*>::iterator iter = _tetra.begin(); iter != _tetra.end(); ++iter)
- {
- (*iter)->splitIntoDualCells(subTetras);
- for(int i=0;i<24;i++)
- {
- SplitterTetra<MyMeshType> *tmp=subTetras[i];
- double volume = tmp->intersectSourceCell(targetCell);
- 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));
- }
- }
- delete tmp;
- }
- }
- }
+ if( _split.getSplittingPolicy() != GENERAL_24)
+ {
+ SplitterTetra<MyMeshType>* subTetras[24];
+ 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.splitTargetCell(*iterCellS,nbOfNodesS,_tetra);
+ 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);
+ AddContributionInRow(resRow,OTT<ConnType,numPol>::indFC(sourceNode),volume);
+ delete tmp;
+ }
+ }
+ }
+ }
+ else
+ {
+ 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.splitTargetCell(*iterCellS,nbOfNodesS,_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(2) ), sourceNode1( (*iter)->getId(3) );
+ AddContributionInRow(resRow,OTT<ConnType,numPol>::indFC(sourceNode0),volume/2.);
+ AddContributionInRow(resRow,OTT<ConnType,numPol>::indFC(sourceNode1),volume/2.);
+ }
+ }
+ }
}
}
typename MyMeshTypeS::MyConnType conn[4];
// get the cell center
conn[0] = 14;
- nodes[0] = getCoordsOfSubNode(conn[0]);
+ typename MyMeshTypeS::MyConnType realConn;
+ 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]);
+ nodes[2] = getCoordsOfSubNode2(conn[2],realConn); conn[2] = realConn;
+ nodes[3] = getCoordsOfSubNode2(conn[3],realConn); conn[3] = realConn;
SplitterTetra<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_src_mesh, nodes, conn);
tetra.push_back(t);