ForkMaster/smesh
1
0
Fork 0
mirror of https://git.salome-platform.org/gitpub/modules/smesh.git synced 2025-02-05 03:54:17 +05:00
Code Issues Packages Projects Releases Wiki Activity

52499: Prismatic mesh is not computed on a prismatic shape

Browse Source 
22701: EDF SMESH: Crash when creating mesh

 1) Search of transformation using a least-square approximation
 2) Compute internal prism nodes using the found transformation
This commit is contained in:
eap 2014-09-24 17:38:27 +04:00
parent 214d7e4bdc
commit b385a0679d
10 changed files with 1238 additions and 364 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
src/SMESHUtils/SMESH_TypeDefs.hxx
View File

@ -43,7 +43,8 @@ typedef std::map<const SMDS_MeshElement*,
std::list<const SMDS_MeshElement*>, TIDCompare > TElemOfElemListMap; std::list<const SMDS_MeshElement*>, TIDCompare > TElemOfElemListMap;
typedef std::map<const SMDS_MeshElement*, typedef std::map<const SMDS_MeshElement*,
std::list<const SMDS_MeshNode*>, TIDCompare > TElemOfNodeListMap; std::list<const SMDS_MeshNode*>, TIDCompare > TElemOfNodeListMap;
typedef std::map<const SMDS_MeshNode*, const SMDS_MeshNode*> TNodeNodeMap; typedef std::map<const SMDS_MeshNode*,
const SMDS_MeshNode*, TIDCompare> TNodeNodeMap;
//!< Set of elements sorted by ID, to be used to assure predictability of edition //!< Set of elements sorted by ID, to be used to assure predictability of edition
typedef std::set< const SMDS_MeshElement*, TIDCompare > TIDSortedElemSet; typedef std::set< const SMDS_MeshElement*, TIDCompare > TIDSortedElemSet;

585
src/StdMeshers/StdMeshers_Prism_3D.cxx
View File

@ -77,7 +77,7 @@ using namespace std;
#define SHOWYXZ(msg, xyz) #define SHOWYXZ(msg, xyz)
#endif #endif
namespace TAssocTool = StdMeshers_ProjectionUtils; namespace NSProjUtils = StdMeshers_ProjectionUtils;
typedef SMESH_Comment TCom; typedef SMESH_Comment TCom;
@ -157,6 +157,10 @@ namespace {
fatherAlgo->GetGen() ); fatherAlgo->GetGen() );
return algo; return algo;
} }
const NSProjUtils::TNodeNodeMap& GetNodesMap()
{
return _src2tgtNodes;
}
}; };
//======================================================================= //=======================================================================
/*! /*!
@ -539,7 +543,7 @@ StdMeshers_Prism_3D::StdMeshers_Prism_3D(int hypId, int studyId, SMESH_Gen* gen)
{ {
_name = "Prism_3D"; _name = "Prism_3D";
_shapeType = (1 << TopAbs_SOLID); // 1 bit per shape type _shapeType = (1 << TopAbs_SOLID); // 1 bit per shape type
_onlyUnaryInput = false; // accept all SOLIDs at once _onlyUnaryInput = false; // mesh all SOLIDs at once
_requireDiscreteBoundary = false; // mesh FACEs and EDGEs by myself _requireDiscreteBoundary = false; // mesh FACEs and EDGEs by myself
_supportSubmeshes = true; // "source" FACE must be meshed by other algo _supportSubmeshes = true; // "source" FACE must be meshed by other algo
_neededLowerHyps[ 1 ] = true; // suppress warning on hiding a global 1D algo _neededLowerHyps[ 1 ] = true; // suppress warning on hiding a global 1D algo
@ -578,12 +582,12 @@ bool StdMeshers_Prism_3D::CheckHypothesis(SMESH_Mesh& a
for ( ; exp.More(); exp.Next() ) { for ( ; exp.More(); exp.Next() ) {
++nbFace; ++nbFace;
const TopoDS_Shape& face = exp.Current(); const TopoDS_Shape& face = exp.Current();
nbEdge = TAssocTool::Count( face, TopAbs_EDGE, 0 ); nbEdge = NSProjUtils::Count( face, TopAbs_EDGE, 0 );
nbWire = TAssocTool::Count( face, TopAbs_WIRE, 0 ); nbWire = NSProjUtils::Count( face, TopAbs_WIRE, 0 );
if ( nbEdge!= 4 || nbWire!= 1 ) { if ( nbEdge!= 4 || nbWire!= 1 ) {
if ( !notQuadFaces.empty() ) { if ( !notQuadFaces.empty() ) {
if ( TAssocTool::Count( notQuadFaces.back(), TopAbs_EDGE, 0 ) != nbEdge || if ( NSProjUtils::Count( notQuadFaces.back(), TopAbs_EDGE, 0 ) != nbEdge ||
TAssocTool::Count( notQuadFaces.back(), TopAbs_WIRE, 0 ) != nbWire ) NSProjUtils::Count( notQuadFaces.back(), TopAbs_WIRE, 0 ) != nbWire )
RETURN_BAD_RESULT("Different not quad faces"); RETURN_BAD_RESULT("Different not quad faces");
} }
notQuadFaces.push_back( face ); notQuadFaces.push_back( face );
@ -595,7 +599,7 @@ bool StdMeshers_Prism_3D::CheckHypothesis(SMESH_Mesh& a
RETURN_BAD_RESULT("Bad nb not quad faces: " << notQuadFaces.size()); RETURN_BAD_RESULT("Bad nb not quad faces: " << notQuadFaces.size());
// check total nb faces // check total nb faces
nbEdge = TAssocTool::Count( notQuadFaces.back(), TopAbs_EDGE, 0 ); nbEdge = NSProjUtils::Count( notQuadFaces.back(), TopAbs_EDGE, 0 );
if ( nbFace != nbEdge + 2 ) if ( nbFace != nbEdge + 2 )
RETURN_BAD_RESULT("Bad nb of faces: " << nbFace << " but must be " << nbEdge + 2); RETURN_BAD_RESULT("Bad nb of faces: " << nbFace << " but must be " << nbEdge + 2);
} }
@ -1057,6 +1061,8 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism)
return false; return false;
// Analyse mesh and geometry to find all block sub-shapes and submeshes // Analyse mesh and geometry to find all block sub-shapes and submeshes
// (after fixing IPAL52499 myBlock is used only as a holder of boundary nodes
// and location of internal nodes is computed by StdMeshers_Sweeper)
if ( !myBlock.Init( myHelper, thePrism )) if ( !myBlock.Init( myHelper, thePrism ))
return toSM( error( myBlock.GetError())); return toSM( error( myBlock.GetError()));
@ -1067,10 +1073,10 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism)
// Try to get gp_Trsf to get all nodes from bottom ones // Try to get gp_Trsf to get all nodes from bottom ones
vector<gp_Trsf> trsf; vector<gp_Trsf> trsf;
gp_Trsf bottomToTopTrsf; gp_Trsf bottomToTopTrsf;
if ( !myBlock.GetLayersTransformation( trsf, thePrism )) // if ( !myBlock.GetLayersTransformation( trsf, thePrism ))
trsf.clear(); // trsf.clear();
else if ( !trsf.empty() ) // else if ( !trsf.empty() )
bottomToTopTrsf = trsf.back(); // bottomToTopTrsf = trsf.back();
// To compute coordinates of a node inside a block, it is necessary to know // To compute coordinates of a node inside a block, it is necessary to know
// 1. normalized parameters of the node by which // 1. normalized parameters of the node by which
@ -1092,31 +1098,30 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism)
// Create nodes inside the block // Create nodes inside the block
// try to use transformation (issue 0020680) // use transformation (issue 0020680, IPAL0052499)
if ( !trsf.empty() ) StdMeshers_Sweeper sweeper;
// load boundary nodes
bool dummy;
list< TopoDS_Edge >::const_iterator edge = thePrism.myBottomEdges.begin();
for ( ; edge != thePrism.myBottomEdges.end(); ++edge )
{ {
// loop on nodes inside the bottom face int edgeID = meshDS->ShapeToIndex( *edge );
TParam2ColumnMap* u2col = const_cast<TParam2ColumnMap*>
( myBlock.GetParam2ColumnMap( edgeID, dummy ));
TParam2ColumnMap::iterator u2colIt = u2col->begin();
for ( ; u2colIt != u2col->end(); ++u2colIt )
sweeper.myBndColumns.push_back( & u2colIt->second );
}
// load node columns inside the bottom face
TNode2ColumnMap::iterator bot_column = myBotToColumnMap.begin(); TNode2ColumnMap::iterator bot_column = myBotToColumnMap.begin();
for ( ; bot_column != myBotToColumnMap.end(); ++bot_column ) for ( ; bot_column != myBotToColumnMap.end(); ++bot_column )
{ sweeper.myIntColumns.push_back( & bot_column->second );
const Prism_3D::TNode& tBotNode = bot_column->first; // bottom TNode
if ( tBotNode.GetPositionType() != SMDS_TOP_FACE )
continue; // node is not inside face
// column nodes; middle part of the column are zero pointers const double tol = getSweepTolerance( thePrism );
TNodeColumn& column = bot_column->second;
TNodeColumn::iterator columnNodes = column.begin();
for ( int z = 0; columnNodes != column.end(); ++columnNodes, ++z)
{
const SMDS_MeshNode* & node = *columnNodes;
if ( node ) continue; // skip bottom or top node
gp_XYZ coords = tBotNode.GetCoords(); if ( sweeper.ComputeNodes( *myHelper, tol ))
trsf[z-1].Transforms( coords ); {
node = meshDS->AddNode( coords.X(), coords.Y(), coords.Z() );
meshDS->SetNodeInVolume( node, volumeID );
}
} // loop on bottom nodes
} }
else // use block approach else // use block approach
{ {
@ -1853,8 +1858,8 @@ void StdMeshers_Prism_3D::AddPrisms( vector<const TNodeColumn*> & columns,
bool StdMeshers_Prism_3D::assocOrProjBottom2Top( const gp_Trsf & bottomToTopTrsf, bool StdMeshers_Prism_3D::assocOrProjBottom2Top( const gp_Trsf & bottomToTopTrsf,
const Prism_3D::TPrismTopo& thePrism) const Prism_3D::TPrismTopo& thePrism)
{ {
SMESH_subMesh * botSM = myBlock.SubMesh( ID_BOT_FACE ); SMESH_subMesh * botSM = myHelper->GetMesh()->GetSubMesh( thePrism.myBottom );
SMESH_subMesh * topSM = myBlock.SubMesh( ID_TOP_FACE ); SMESH_subMesh * topSM = myHelper->GetMesh()->GetSubMesh( thePrism.myTop );
SMESHDS_SubMesh * botSMDS = botSM->GetSubMeshDS(); SMESHDS_SubMesh * botSMDS = botSM->GetSubMeshDS();
SMESHDS_SubMesh * topSMDS = topSM->GetSubMeshDS(); SMESHDS_SubMesh * topSMDS = topSM->GetSubMeshDS();
@ -1885,18 +1890,22 @@ bool StdMeshers_Prism_3D::assocOrProjBottom2Top( const gp_Trsf & bottomToTopTrsf
<<" and #"<< topSM->GetId() << " seems different" )); <<" and #"<< topSM->GetId() << " seems different" ));
///RETURN_BAD_RESULT("Need to project but not allowed"); ///RETURN_BAD_RESULT("Need to project but not allowed");
NSProjUtils::TNodeNodeMap n2nMap;
const NSProjUtils::TNodeNodeMap* n2nMapPtr = & n2nMap;
if ( needProject ) if ( needProject )
{ {
return projectBottomToTop( bottomToTopTrsf ); if ( !projectBottomToTop( bottomToTopTrsf, thePrism ))
return false;
n2nMapPtr = & TProjction2dAlgo::instance( this )->GetNodesMap();
} }
TopoDS_Face botFace = TopoDS::Face( myBlock.Shape( ID_BOT_FACE )); if ( !n2nMapPtr || n2nMapPtr->size() < botSMDS->NbNodes() )
TopoDS_Face topFace = TopoDS::Face( myBlock.Shape( ID_TOP_FACE )); {
// associate top and bottom faces // associate top and bottom faces
TAssocTool::TShapeShapeMap shape2ShapeMap; NSProjUtils::TShapeShapeMap shape2ShapeMap;
const bool sameTopo = const bool sameTopo =
TAssocTool::FindSubShapeAssociation( botFace, myBlock.Mesh(), NSProjUtils::FindSubShapeAssociation( thePrism.myBottom, myHelper->GetMesh(),
topFace, myBlock.Mesh(), thePrism.myTop, myHelper->GetMesh(),
shape2ShapeMap); shape2ShapeMap);
if ( !sameTopo ) if ( !sameTopo )
for ( size_t iQ = 0; iQ < thePrism.myWallQuads.size(); ++iQ ) for ( size_t iQ = 0; iQ < thePrism.myWallQuads.size(); ++iQ )
@ -1908,9 +1917,9 @@ bool StdMeshers_Prism_3D::assocOrProjBottom2Top( const gp_Trsf & bottomToTopTrsf
{ {
for ( int iE = 0; iE < botSide->NbEdges(); ++iE ) for ( int iE = 0; iE < botSide->NbEdges(); ++iE )
{ {
TAssocTool::InsertAssociation( botSide->Edge( iE ), NSProjUtils::InsertAssociation( botSide->Edge( iE ),
topSide->Edge( iE ), shape2ShapeMap ); topSide->Edge( iE ), shape2ShapeMap );
TAssocTool::InsertAssociation( myHelper->IthVertex( 0, botSide->Edge( iE )), NSProjUtils::InsertAssociation( myHelper->IthVertex( 0, botSide->Edge( iE )),
myHelper->IthVertex( 0, topSide->Edge( iE )), myHelper->IthVertex( 0, topSide->Edge( iE )),
shape2ShapeMap ); shape2ShapeMap );
} }
@ -1926,9 +1935,9 @@ bool StdMeshers_Prism_3D::assocOrProjBottom2Top( const gp_Trsf & bottomToTopTrsf
if ( vb.IsSame( sideB->FirstVertex() ) && if ( vb.IsSame( sideB->FirstVertex() ) &&
vt.IsSame( sideT->LastVertex() )) vt.IsSame( sideT->LastVertex() ))
{ {
TAssocTool::InsertAssociation( botSide->Edge( 0 ), NSProjUtils::InsertAssociation( botSide->Edge( 0 ),
topSide->Edge( 0 ), shape2ShapeMap ); topSide->Edge( 0 ), shape2ShapeMap );
TAssocTool::InsertAssociation( vb, vt, shape2ShapeMap ); NSProjUtils::InsertAssociation( vb, vt, shape2ShapeMap );
} }
vb = myHelper->IthVertex( 1, botSide->Edge( botSide->NbEdges()-1 )); vb = myHelper->IthVertex( 1, botSide->Edge( botSide->NbEdges()-1 ));
vt = myHelper->IthVertex( 1, topSide->Edge( topSide->NbEdges()-1 )); vt = myHelper->IthVertex( 1, topSide->Edge( topSide->NbEdges()-1 ));
@ -1937,18 +1946,18 @@ bool StdMeshers_Prism_3D::assocOrProjBottom2Top( const gp_Trsf & bottomToTopTrsf
if ( vb.IsSame( sideB->FirstVertex() ) && if ( vb.IsSame( sideB->FirstVertex() ) &&
vt.IsSame( sideT->LastVertex() )) vt.IsSame( sideT->LastVertex() ))
{ {
TAssocTool::InsertAssociation( botSide->Edge( botSide->NbEdges()-1 ), NSProjUtils::InsertAssociation( botSide->Edge( botSide->NbEdges()-1 ),
topSide->Edge( topSide->NbEdges()-1 ), topSide->Edge( topSide->NbEdges()-1 ),
shape2ShapeMap ); shape2ShapeMap );
TAssocTool::InsertAssociation( vb, vt, shape2ShapeMap ); NSProjUtils::InsertAssociation( vb, vt, shape2ShapeMap );
} }
} }
} }
// Find matching nodes of top and bottom faces // Find matching nodes of top and bottom faces
TNodeNodeMap n2nMap; n2nMapPtr = & n2nMap;
if ( ! TAssocTool::FindMatchingNodesOnFaces( botFace, myBlock.Mesh(), if ( ! NSProjUtils::FindMatchingNodesOnFaces( thePrism.myBottom, myHelper->GetMesh(),
topFace, myBlock.Mesh(), thePrism.myTop, myHelper->GetMesh(),
shape2ShapeMap, n2nMap )) shape2ShapeMap, n2nMap ))
{ {
if ( sameTopo ) if ( sameTopo )
@ -1958,13 +1967,13 @@ bool StdMeshers_Prism_3D::assocOrProjBottom2Top( const gp_Trsf & bottomToTopTrsf
return toSM( error(TCom("Topology of faces #") << botSM->GetId() return toSM( error(TCom("Topology of faces #") << botSM->GetId()
<<" and #"<< topSM->GetId() << " seems different" )); <<" and #"<< topSM->GetId() << " seems different" ));
} }
}
// Fill myBotToColumnMap // Fill myBotToColumnMap
int zSize = myBlock.VerticalSize(); int zSize = myBlock.VerticalSize();
//TNode prevTNode; TNodeNodeMap::const_iterator bN_tN = n2nMapPtr->begin();
TNodeNodeMap::iterator bN_tN = n2nMap.begin(); for ( ; bN_tN != n2nMapPtr->end(); ++bN_tN )
for ( ; bN_tN != n2nMap.end(); ++bN_tN )
{ {
const SMDS_MeshNode* botNode = bN_tN->first; const SMDS_MeshNode* botNode = bN_tN->first;
const SMDS_MeshNode* topNode = bN_tN->second; const SMDS_MeshNode* topNode = bN_tN->second;
@ -1984,17 +1993,22 @@ bool StdMeshers_Prism_3D::assocOrProjBottom2Top( const gp_Trsf & bottomToTopTrsf
//================================================================================ //================================================================================
/*! /*!
* \brief Remove quadrangles from the top face and * \brief Remove faces from the top face and re-create them by projection from the bottom
* create triangles there by projection from the bottom
* \retval bool - a success or not * \retval bool - a success or not
*/ */
//================================================================================ //================================================================================
bool StdMeshers_Prism_3D::projectBottomToTop( const gp_Trsf & bottomToTopTrsf ) bool StdMeshers_Prism_3D::projectBottomToTop( const gp_Trsf & bottomToTopTrsf,
const Prism_3D::TPrismTopo& thePrism )
{ {
SMESHDS_Mesh* meshDS = myBlock.MeshDS(); if ( project2dMesh( thePrism.myBottom, thePrism.myTop ))
SMESH_subMesh * botSM = myBlock.SubMesh( ID_BOT_FACE ); {
SMESH_subMesh * topSM = myBlock.SubMesh( ID_TOP_FACE ); return true;
}
SMESHDS_Mesh* meshDS = myHelper->GetMeshDS();
SMESH_subMesh * botSM = myHelper->GetMesh()->GetSubMesh( thePrism.myBottom );
SMESH_subMesh * topSM = myHelper->GetMesh()->GetSubMesh( thePrism.myTop );
SMESHDS_SubMesh * botSMDS = botSM->GetSubMeshDS(); SMESHDS_SubMesh * botSMDS = botSM->GetSubMeshDS();
SMESHDS_SubMesh * topSMDS = topSM->GetSubMeshDS(); SMESHDS_SubMesh * topSMDS = topSM->GetSubMeshDS();
@ -2002,9 +2016,9 @@ bool StdMeshers_Prism_3D::projectBottomToTop( const gp_Trsf & bottomToTopTrsf )
if ( topSMDS && topSMDS->NbElements() > 0 ) if ( topSMDS && topSMDS->NbElements() > 0 )
topSM->ComputeStateEngine( SMESH_subMesh::CLEAN ); topSM->ComputeStateEngine( SMESH_subMesh::CLEAN );
const TopoDS_Face& botFace = TopoDS::Face( myBlock.Shape( ID_BOT_FACE )); // oriented within const TopoDS_Face& botFace = thePrism.myBottom; // oriented within
const TopoDS_Face& topFace = TopoDS::Face( myBlock.Shape( ID_TOP_FACE )); // the 3D SHAPE const TopoDS_Face& topFace = thePrism.myTop; // the 3D SHAPE
int topFaceID = meshDS->ShapeToIndex( topFace ); int topFaceID = meshDS->ShapeToIndex( thePrism.myTop );
SMESH_MesherHelper botHelper( *myHelper->GetMesh() ); SMESH_MesherHelper botHelper( *myHelper->GetMesh() );
botHelper.SetSubShape( botFace ); botHelper.SetSubShape( botFace );
@ -2071,6 +2085,9 @@ bool StdMeshers_Prism_3D::projectBottomToTop( const gp_Trsf & bottomToTopTrsf )
column.resize( zSize ); column.resize( zSize );
column.front() = botNode; column.front() = botNode;
column.back() = topNode; column.back() = topNode;
if ( _computeCanceled )
return toSM( error( SMESH_ComputeError::New(COMPERR_CANCELED)));
} }
// Create top faces // Create top faces
@ -2135,6 +2152,77 @@ bool StdMeshers_Prism_3D::projectBottomToTop( const gp_Trsf & bottomToTopTrsf )
return true; return true;
} }
//=======================================================================
//function : getSweepTolerance
//purpose : Compute tolerance to pass to StdMeshers_Sweeper
//=======================================================================
double StdMeshers_Prism_3D::getSweepTolerance( const Prism_3D::TPrismTopo& thePrism )
{
SMESHDS_Mesh* meshDS = myHelper->GetMeshDS();
SMESHDS_SubMesh * sm[2] = { meshDS->MeshElements( thePrism.myBottom ),
meshDS->MeshElements( thePrism.myTop ) };
double minDist = 1e100;
vector< SMESH_TNodeXYZ > nodes;
for ( int iSM = 0; iSM < 2; ++iSM )
{
if ( !sm[ iSM ]) continue;
SMDS_ElemIteratorPtr fIt = sm[ iSM ]->GetElements();
while ( fIt->more() )
{
const SMDS_MeshElement* face = fIt->next();
const int nbNodes = face->NbCornerNodes();
SMDS_ElemIteratorPtr nIt = face->nodesIterator();
nodes.resize( nbNodes + 1 );
for ( int iN = 0; iN < nbNodes; ++iN )
nodes[ iN ] = nIt->next();
nodes.back() = nodes[0];
// loop on links
double dist2;
for ( int iN = 0; iN < nbNodes; ++iN )
{
if ( nodes[ iN ]._node->GetPosition()->GetDim() < 2 &&
nodes[ iN+1 ]._node->GetPosition()->GetDim() < 2 )
{
// it's a boundary link; measure distance of other
// nodes to this link
gp_XYZ linkDir = nodes[ iN ] - nodes[ iN+1 ];
double linkLen = linkDir.Modulus();
bool isDegen = ( linkLen < numeric_limits<double>::min() );
if ( !isDegen ) linkDir /= linkLen;
for ( int iN2 = 0; iN2 < nbNodes; ++iN2 ) // loop on other nodes
{
if ( nodes[ iN2 ] == nodes[ iN ] ||
nodes[ iN2 ] == nodes[ iN+1 ]) continue;
if ( isDegen )
{
dist2 = ( nodes[ iN ] - nodes[ iN2 ]).SquareModulus();
}
else
{
dist2 = linkDir.CrossSquareMagnitude( nodes[ iN ] - nodes[ iN2 ]);
}
if ( dist2 > numeric_limits<double>::min() )
minDist = Min ( minDist, dist2 );
}
}
// measure length link
else if ( nodes[ iN ]._node < nodes[ iN+1 ]._node ) // not to measure same link twice
{
dist2 = ( nodes[ iN ] - nodes[ iN+1 ]).SquareModulus();
if ( dist2 > numeric_limits<double>::min() )
minDist = Min ( minDist, dist2 );
}
}
}
}
return 0.1 * Sqrt ( minDist );
}
//======================================================================= //=======================================================================
//function : project2dMesh //function : project2dMesh
//purpose : Project mesh faces from a source FACE of one prism (theSrcFace) //purpose : Project mesh faces from a source FACE of one prism (theSrcFace)
@ -2261,6 +2349,12 @@ namespace // utils used by StdMeshers_Prism_3D::IsApplicable()
if ( E.IsSame( Edge( i ))) return i; if ( E.IsSame( Edge( i ))) return i;
return -1; return -1;
} }
bool IsSideFace( const TopoDS_Shape& face ) const
{
if ( _faces->Contains( face )) // avoid returning true for a prism top FACE
return ( !_face.IsNull() || !( face.IsSame( _faces->FindKey( _faces->Extent() ))));
return false;
}
}; };
//-------------------------------------------------------------------------------- //--------------------------------------------------------------------------------
/*! /*!
@ -2392,12 +2486,11 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA
TEdgeWithNeighborsVec& botEdges = faceEdgesVec[ iF ]; TEdgeWithNeighborsVec& botEdges = faceEdgesVec[ iF ];
if ( botEdges.empty() ) if ( botEdges.empty() )
{
if ( !getEdges( botF, botEdges, /*noHoles=*/false )) if ( !getEdges( botF, botEdges, /*noHoles=*/false ))
break; break;
if ( allFaces.Extent()-1 <= (int) botEdges.size() ) if ( allFaces.Extent()-1 <= (int) botEdges.size() )
continue; // all faces are adjacent to botF - no top FACE continue; // all faces are adjacent to botF - no top FACE
}
// init data of side FACEs // init data of side FACEs
vector< PrismSide > sides( botEdges.size() ); vector< PrismSide > sides( botEdges.size() );
for ( int iS = 0; iS < botEdges.size(); ++iS ) for ( int iS = 0; iS < botEdges.size(); ++iS )
@ -2411,7 +2504,7 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA
} }
bool isOK = true; // ok for a current botF bool isOK = true; // ok for a current botF
bool isAdvanced = true; bool isAdvanced = true; // is new data found in a current loop
int nbFoundSideFaces = 0; int nbFoundSideFaces = 0;
for ( int iLoop = 0; isOK && isAdvanced; ++iLoop ) for ( int iLoop = 0; isOK && isAdvanced; ++iLoop )
{ {
@ -2420,7 +2513,8 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA
{ {
PrismSide& side = sides[ iS ]; PrismSide& side = sides[ iS ];
if ( side._face.IsNull() ) if ( side._face.IsNull() )
continue; continue; // probably the prism top face is the last of side._faces
if ( side._topEdge.IsNull() ) if ( side._topEdge.IsNull() )
{ {
// find vertical EDGEs --- EGDEs shared with neighbor side FACEs // find vertical EDGEs --- EGDEs shared with neighbor side FACEs
@ -2434,7 +2528,7 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA
if ( side._isCheckedEdge[ iE ] ) continue; if ( side._isCheckedEdge[ iE ] ) continue;
const TopoDS_Edge& vertE = side.Edge( iE ); const TopoDS_Edge& vertE = side.Edge( iE );
const TopoDS_Shape& neighborF = getAnotherFace( side._face, vertE, facesOfEdge ); const TopoDS_Shape& neighborF = getAnotherFace( side._face, vertE, facesOfEdge );
bool isEdgeShared = adjSide->_faces->Contains( neighborF ); bool isEdgeShared = adjSide->IsSideFace( neighborF );
if ( isEdgeShared ) if ( isEdgeShared )
{ {
isAdvanced = true; isAdvanced = true;
@ -2480,13 +2574,13 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA
{ {
if ( side._leftSide->_faces->Contains( f )) if ( side._leftSide->_faces->Contains( f ))
{ {
stop = true; stop = true; // probably f is the prism top face
side._leftSide->_face.Nullify(); side._leftSide->_face.Nullify();
side._leftSide->_topEdge.Nullify(); side._leftSide->_topEdge.Nullify();
} }
if ( side._rightSide->_faces->Contains( f )) if ( side._rightSide->_faces->Contains( f ))
{ {
stop = true; stop = true; // probably f is the prism top face
side._rightSide->_face.Nullify(); side._rightSide->_face.Nullify();
side._rightSide->_topEdge.Nullify(); side._rightSide->_topEdge.Nullify();
} }
@ -2805,8 +2899,9 @@ bool StdMeshers_Prism_3D::initPrism(Prism_3D::TPrismTopo& thePrism,
thePrism.myShape3D = shape3D; thePrism.myShape3D = shape3D;
if ( thePrism.myBottom.IsNull() ) if ( thePrism.myBottom.IsNull() )
thePrism.myBottom = TopoDS::Face( botSM->GetSubShape() ); thePrism.myBottom = TopoDS::Face( botSM->GetSubShape() );
thePrism.myBottom.Orientation( myHelper->GetSubShapeOri( shape3D, thePrism.myBottom.Orientation( myHelper->GetSubShapeOri( shape3D, thePrism.myBottom ));
thePrism.myBottom )); thePrism.myTop. Orientation( myHelper->GetSubShapeOri( shape3D, thePrism.myTop ));
// Get ordered bottom edges // Get ordered bottom edges
TopoDS_Face reverseBottom = // to have order of top EDGEs as in the top FACE TopoDS_Face reverseBottom = // to have order of top EDGEs as in the top FACE
TopoDS::Face( thePrism.myBottom.Reversed() ); TopoDS::Face( thePrism.myBottom.Reversed() );
@ -4221,3 +4316,355 @@ gp_Pnt2d StdMeshers_PrismAsBlock::TPCurveOnHorFaceAdaptor::Value(const Standard_
double r = ( U - i1->first ) / ( i2->first - i1->first ); double r = ( U - i1->first ) / ( i2->first - i1->first );
return i1->second * ( 1 - r ) + i2->second * r; return i1->second * ( 1 - r ) + i2->second * r;
} }
//================================================================================
/*!
* \brief Projects internal nodes using transformation found by boundary nodes
*/
//================================================================================
bool StdMeshers_Sweeper::projectIntPoints(const vector< gp_XYZ >& fromBndPoints,
const vector< gp_XYZ >& toBndPoints,
const vector< gp_XYZ >& fromIntPoints,
vector< gp_XYZ >& toIntPoints,
NSProjUtils::TrsfFinder3D& trsf,
vector< gp_XYZ > * bndError)
{
// find transformation
if ( trsf.IsIdentity() && !trsf.Solve( fromBndPoints, toBndPoints ))
return false;
// compute internal points using the found trsf
for ( size_t iP = 0; iP < fromIntPoints.size(); ++iP )
{
toIntPoints[ iP ] = trsf.Transform( fromIntPoints[ iP ]);
}
// compute boundary error
if ( bndError )
{
bndError->resize( fromBndPoints.size() );
gp_XYZ fromTrsf;
for ( size_t iP = 0; iP < fromBndPoints.size(); ++iP )
{
fromTrsf = trsf.Transform( fromBndPoints[ iP ] );
(*bndError)[ iP ] = toBndPoints[ iP ] - fromTrsf;
}
}
return true;
}
//================================================================================
/*!
* \brief Add boundary error to ineternal points
*/
//================================================================================
void StdMeshers_Sweeper::applyBoundaryError(const vector< gp_XYZ >& bndPoints,
const vector< gp_XYZ >& bndError1,
const vector< gp_XYZ >& bndError2,
const double r,
vector< gp_XYZ >& intPoints,
vector< double >& int2BndDist)
{
// fix each internal point
const double eps = 1e-100;
for ( size_t iP = 0; iP < intPoints.size(); ++iP )
{
gp_XYZ & intPnt = intPoints[ iP ];
// compute distance from intPnt to each boundary node
double int2BndDistSum = 0;
for ( size_t iBnd = 0; iBnd < bndPoints.size(); ++iBnd )
{
int2BndDist[ iBnd ] = 1 / (( intPnt - bndPoints[ iBnd ]).SquareModulus() + eps );
int2BndDistSum += int2BndDist[ iBnd ];
}
// apply bndError
for ( size_t iBnd = 0; iBnd < bndPoints.size(); ++iBnd )
{
intPnt += bndError1[ iBnd ] * ( 1 - r ) * int2BndDist[ iBnd ] / int2BndDistSum;
intPnt += bndError2[ iBnd ] * r * int2BndDist[ iBnd ] / int2BndDistSum;
}
}
}
//================================================================================
/*!
* \brief Creates internal nodes of the prism
*/
//================================================================================
bool StdMeshers_Sweeper::ComputeNodes( SMESH_MesherHelper& helper,
const double tol)
{
const size_t zSize = myBndColumns[0]->size();
const size_t zSrc = 0, zTgt = zSize-1;
if ( zSize < 3 ) return true;
vector< vector< gp_XYZ > > intPntsOfLayer( zSize ); // node coodinates to compute
// set coordinates of src and tgt nodes
for ( size_t z = 0; z < intPntsOfLayer.size(); ++z )
intPntsOfLayer[ z ].resize( myIntColumns.size() );
for ( size_t iP = 0; iP < myIntColumns.size(); ++iP )
{
intPntsOfLayer[ zSrc ][ iP ] = intPoint( iP, zSrc );
intPntsOfLayer[ zTgt ][ iP ] = intPoint( iP, zTgt );
}
// compute coordinates of internal nodes by projecting (transfroming) src and tgt
// nodes towards the central layer
vector< NSProjUtils::TrsfFinder3D > trsfOfLayer( zSize );
vector< vector< gp_XYZ > > bndError( zSize );
// boundary points used to compute an affine transformation from a layer to a next one
vector< gp_XYZ > fromSrcBndPnts( myBndColumns.size() ), fromTgtBndPnts( myBndColumns.size() );
vector< gp_XYZ > toSrcBndPnts ( myBndColumns.size() ), toTgtBndPnts ( myBndColumns.size() );
for ( size_t iP = 0; iP < myBndColumns.size(); ++iP )
{
fromSrcBndPnts[ iP ] = bndPoint( iP, zSrc );
fromTgtBndPnts[ iP ] = bndPoint( iP, zTgt );
}
size_t zS = zSrc + 1;
size_t zT = zTgt - 1;
for ( ; zS < zT; ++zS, --zT ) // vertical loop on layers
{
for ( size_t iP = 0; iP < myBndColumns.size(); ++iP )
{
toSrcBndPnts[ iP ] = bndPoint( iP, zS );
toTgtBndPnts[ iP ] = bndPoint( iP, zT );
}
if (! projectIntPoints( fromSrcBndPnts, toSrcBndPnts,
intPntsOfLayer[ zS-1 ], intPntsOfLayer[ zS ],
trsfOfLayer [ zS-1 ], & bndError[ zS-1 ]))
return false;
if (! projectIntPoints( fromTgtBndPnts, toTgtBndPnts,
intPntsOfLayer[ zT+1 ], intPntsOfLayer[ zT ],
trsfOfLayer [ zT+1 ], & bndError[ zT+1 ]))
return false;
// if ( zT == zTgt - 1 )
// {
// for ( size_t iP = 0; iP < myBndColumns.size(); ++iP )
// {
// gp_XYZ fromTrsf = trsfOfLayer [ zT+1].Transform( fromTgtBndPnts[ iP ] );
// cout << "mesh.AddNode( "
// << fromTrsf.X() << ", "
// << fromTrsf.Y() << ", "
// << fromTrsf.Z() << ") " << endl;
// }
// for ( size_t iP = 0; iP < myIntColumns.size(); ++iP )
// cout << "mesh.AddNode( "
// << intPntsOfLayer[ zT ][ iP ].X() << ", "
// << intPntsOfLayer[ zT ][ iP ].Y() << ", "
// << intPntsOfLayer[ zT ][ iP ].Z() << ") " << endl;
// }
fromTgtBndPnts.swap( toTgtBndPnts );
fromSrcBndPnts.swap( toSrcBndPnts );
}
// Compute two projections of internal points to the central layer
// in order to evaluate an error of internal points
bool centerIntErrorIsSmall;
vector< gp_XYZ > centerSrcIntPnts( myIntColumns.size() );
vector< gp_XYZ > centerTgtIntPnts( myIntColumns.size() );
for ( size_t iP = 0; iP < myBndColumns.size(); ++iP )
{
toSrcBndPnts[ iP ] = bndPoint( iP, zS );
toTgtBndPnts[ iP ] = bndPoint( iP, zT );
}
if (! projectIntPoints( fromSrcBndPnts, toSrcBndPnts,
intPntsOfLayer[ zS-1 ], centerSrcIntPnts,
trsfOfLayer [ zS-1 ], & bndError[ zS-1 ]))
return false;
if (! projectIntPoints( fromTgtBndPnts, toTgtBndPnts,
intPntsOfLayer[ zT+1 ], centerTgtIntPnts,
trsfOfLayer [ zT+1 ], & bndError[ zT+1 ]))
return false;
// evaluate an error of internal points on the central layer
centerIntErrorIsSmall = true;
if ( zS == zT ) // odd zSize
{
for ( size_t iP = 0; ( iP < myIntColumns.size() && centerIntErrorIsSmall ); ++iP )
centerIntErrorIsSmall =
(centerSrcIntPnts[ iP ] - centerTgtIntPnts[ iP ]).SquareModulus() < tol*tol;
}
else // even zSize
{
for ( size_t iP = 0; ( iP < myIntColumns.size() && centerIntErrorIsSmall ); ++iP )
centerIntErrorIsSmall =
(intPntsOfLayer[ zS-1 ][ iP ] - centerTgtIntPnts[ iP ]).SquareModulus() < tol*tol;
}
// Evaluate an error of boundary points
bool bndErrorIsSmall = true;
for ( size_t iP = 0; ( iP < myBndColumns.size() && bndErrorIsSmall ); ++iP )
{
double sumError = 0;
for ( size_t z = 1; z < zS; ++z ) // loop on layers
sumError += ( bndError[ z-1 ][ iP ].Modulus() +
bndError[ zSize-z ][ iP ].Modulus() );
bndErrorIsSmall = ( sumError < tol );
}
// compute final points on the central layer
std::vector< double > int2BndDist( myBndColumns.size() ); // work array of applyBoundaryError()
double r = zS / ( zSize - 1.);
if ( zS == zT )
{
for ( size_t iP = 0; iP < myIntColumns.size(); ++iP )
{
intPntsOfLayer[ zS ][ iP ] =
( 1 - r ) * centerSrcIntPnts[ iP ] + r * centerTgtIntPnts[ iP ];
}
if ( !bndErrorIsSmall )
{
applyBoundaryError( toSrcBndPnts, bndError[ zS-1 ], bndError[ zS+1 ], r,
intPntsOfLayer[ zS ], int2BndDist );
}
}
else
{
for ( size_t iP = 0; iP < myIntColumns.size(); ++iP )
{
intPntsOfLayer[ zS ][ iP ] =
r * intPntsOfLayer[ zS ][ iP ] + ( 1 - r ) * centerSrcIntPnts[ iP ];
intPntsOfLayer[ zT ][ iP ] =
r * intPntsOfLayer[ zT ][ iP ] + ( 1 - r ) * centerTgtIntPnts[ iP ];
}
if ( !bndErrorIsSmall )
{
applyBoundaryError( toSrcBndPnts, bndError[ zS-1 ], bndError[ zS+1 ], r,
intPntsOfLayer[ zS ], int2BndDist );
applyBoundaryError( toTgtBndPnts, bndError[ zT+1 ], bndError[ zT-1 ], r,
intPntsOfLayer[ zT ], int2BndDist );
}
}
//centerIntErrorIsSmall = true;
//bndErrorIsSmall = true;
if ( !centerIntErrorIsSmall )
{
// Compensate the central error; continue adding projection
// by going from central layer to the source and target ones
vector< gp_XYZ >& fromSrcIntPnts = centerSrcIntPnts;
vector< gp_XYZ >& fromTgtIntPnts = centerTgtIntPnts;
vector< gp_XYZ > toSrcIntPnts( myIntColumns.size() );
vector< gp_XYZ > toTgtIntPnts( myIntColumns.size() );
vector< gp_XYZ > srcBndError( myBndColumns.size() );
vector< gp_XYZ > tgtBndError( myBndColumns.size() );
fromTgtBndPnts.swap( toTgtBndPnts );
fromSrcBndPnts.swap( toSrcBndPnts );
for ( ++zS, --zT; zS < zTgt; ++zS, --zT ) // vertical loop on layers
{
// invert transformation
if ( !trsfOfLayer[ zS+1 ].Invert() )
trsfOfLayer[ zS+1 ] = NSProjUtils::TrsfFinder3D(); // to recompute
if ( !trsfOfLayer[ zT-1 ].Invert() )
trsfOfLayer[ zT-1 ] = NSProjUtils::TrsfFinder3D();
// project internal nodes and compute bnd error
for ( size_t iP = 0; iP < myBndColumns.size(); ++iP )
{
toSrcBndPnts[ iP ] = bndPoint( iP, zS );
toTgtBndPnts[ iP ] = bndPoint( iP, zT );
}
projectIntPoints( fromSrcBndPnts, toSrcBndPnts,
fromSrcIntPnts, toSrcIntPnts,
trsfOfLayer[ zS+1 ], & srcBndError );
projectIntPoints( fromTgtBndPnts, toTgtBndPnts,
fromTgtIntPnts, toTgtIntPnts,
trsfOfLayer[ zT-1 ], & tgtBndError );
// if ( zS == zTgt - 1 )
// {
// cout << "mesh2 = smesh.Mesh()" << endl;
// for ( size_t iP = 0; iP < myBndColumns.size(); ++iP )
// {
// gp_XYZ fromTrsf = trsfOfLayer [ zS+1].Transform( fromSrcBndPnts[ iP ] );
// cout << "mesh2.AddNode( "
// << fromTrsf.X() << ", "
// << fromTrsf.Y() << ", "
// << fromTrsf.Z() << ") " << endl;
// }
// for ( size_t iP = 0; iP < myIntColumns.size(); ++iP )
// cout << "mesh2.AddNode( "
// << toSrcIntPnts[ iP ].X() << ", "
// << toSrcIntPnts[ iP ].Y() << ", "
// << toSrcIntPnts[ iP ].Z() << ") " << endl;
// }
// sum up 2 projections
r = zS / ( zSize - 1.);
vector< gp_XYZ >& zSIntPnts = intPntsOfLayer[ zS ];
vector< gp_XYZ >& zTIntPnts = intPntsOfLayer[ zT ];
for ( size_t iP = 0; iP < myIntColumns.size(); ++iP )
{
zSIntPnts[ iP ] = r * zSIntPnts[ iP ] + ( 1 - r ) * toSrcIntPnts[ iP ];
zTIntPnts[ iP ] = r * zTIntPnts[ iP ] + ( 1 - r ) * toTgtIntPnts[ iP ];
}
// compensate bnd error
if ( !bndErrorIsSmall )
{
applyBoundaryError( toSrcBndPnts, srcBndError, bndError[ zS+1 ], r,
intPntsOfLayer[ zS ], int2BndDist );
applyBoundaryError( toTgtBndPnts, tgtBndError, bndError[ zT-1 ], r,
intPntsOfLayer[ zT ], int2BndDist );
}
fromSrcBndPnts.swap( toSrcBndPnts );
fromSrcIntPnts.swap( toSrcIntPnts );
fromTgtBndPnts.swap( toTgtBndPnts );
fromTgtIntPnts.swap( toTgtIntPnts );
}
} // if ( !centerIntErrorIsSmall )
else if ( !bndErrorIsSmall )
{
zS = zSrc + 1;
zT = zTgt - 1;
for ( ; zS < zT; ++zS, --zT ) // vertical loop on layers
{
for ( size_t iP = 0; iP < myBndColumns.size(); ++iP )
{
toSrcBndPnts[ iP ] = bndPoint( iP, zS );
toTgtBndPnts[ iP ] = bndPoint( iP, zT );
}
// compensate bnd error
applyBoundaryError( toSrcBndPnts, bndError[ zS-1 ], bndError[ zS-1 ], 0.5,
intPntsOfLayer[ zS ], int2BndDist );
applyBoundaryError( toTgtBndPnts, bndError[ zT+1 ], bndError[ zT+1 ], 0.5,
intPntsOfLayer[ zT ], int2BndDist );
}
}
// cout << "centerIntErrorIsSmall = " << centerIntErrorIsSmall<< endl;
// cout << "bndErrorIsSmall = " << bndErrorIsSmall<< endl;
// Create nodes
for ( size_t iP = 0; iP < myIntColumns.size(); ++iP )
{
vector< const SMDS_MeshNode* > & nodeCol = *myIntColumns[ iP ];
for ( size_t z = zSrc + 1; z < zTgt; ++z ) // vertical loop on layers
{
const gp_XYZ & xyz = intPntsOfLayer[ z ][ iP ];
if ( !( nodeCol[ z ] = helper.AddNode( xyz.X(), xyz.Y(), xyz.Z() )))
return false;
}
}
return true;
}

49
src/StdMeshers/StdMeshers_Prism_3D.hxx
View File

@ -56,6 +56,10 @@ namespace Prism_3D
struct TNode; struct TNode;
struct TPrismTopo; struct TPrismTopo;
} }
namespace StdMeshers_ProjectionUtils
{
class TrsfFinder3D;
}
class SMESHDS_SubMesh; class SMESHDS_SubMesh;
class TopoDS_Edge; class TopoDS_Edge;
@ -398,7 +402,42 @@ private:
}; // class StdMeshers_PrismAsBlock }; // class StdMeshers_PrismAsBlock
// ============================================= // ===============================================
/*!
* \brief Tool building internal nodes in a prism
*/
struct StdMeshers_Sweeper
{
std::vector< TNodeColumn* > myBndColumns; // boundary nodes
std::vector< TNodeColumn* > myIntColumns; // internal nodes
bool ComputeNodes( SMESH_MesherHelper& helper,
const double tol );
private:
gp_XYZ bndPoint( int iP, int z ) const
{ return SMESH_TNodeXYZ( (*myBndColumns[ iP ])[ z ]); }
gp_XYZ intPoint( int iP, int z ) const
{ return SMESH_TNodeXYZ( (*myIntColumns[ iP ])[ z ]); }
static bool projectIntPoints(const std::vector< gp_XYZ >& fromBndPoints,
const std::vector< gp_XYZ >& toBndPoints,
const std::vector< gp_XYZ >& fromIntPoints,
std::vector< gp_XYZ >& toIntPoints,
StdMeshers_ProjectionUtils::TrsfFinder3D& trsf,
std::vector< gp_XYZ > * bndError);
static void applyBoundaryError(const std::vector< gp_XYZ >& bndPoints,
const std::vector< gp_XYZ >& bndError1,
const std::vector< gp_XYZ >& bndError2,
const double r,
std::vector< gp_XYZ >& toIntPoints,
std::vector< double >& int2BndDist);
};
// ===============================================
/*! /*!
* \brief Algo building prisms on a prism shape * \brief Algo building prisms on a prism shape
*/ */
@ -480,7 +519,13 @@ public:
* create triangles there by projection from the bottom * create triangles there by projection from the bottom
* \retval bool - a success or not * \retval bool - a success or not
*/ */
bool projectBottomToTop( const gp_Trsf & bottomToTopTrsf ); bool projectBottomToTop( const gp_Trsf & bottomToTopTrsf,
const Prism_3D::TPrismTopo& thePrism );
/*!
* \brief Compute tolerance to pass to StdMeshers_Sweeper
*/
double getSweepTolerance( const Prism_3D::TPrismTopo& thePrism );
/*! /*!
* \brief Project mesh faces from a source FACE of one prism to * \brief Project mesh faces from a source FACE of one prism to

229
src/StdMeshers/StdMeshers_ProjectionUtils.cxx
View File

@ -66,6 +66,7 @@
#include <TopoDS_Shape.hxx> #include <TopoDS_Shape.hxx>
#include <gp_Pnt.hxx> #include <gp_Pnt.hxx>
#include <gp_Vec.hxx> #include <gp_Vec.hxx>
#include <math_Gauss.hxx>
#include <numeric> #include <numeric>
#include <limits> #include <limits>
@ -98,7 +99,7 @@ namespace HERE = StdMeshers_ProjectionUtils;
namespace { namespace {
static SMESHDS_Mesh* theMeshDS[2] = { 0, 0 }; // used to debug only static SMESHDS_Mesh* theMeshDS[2] = { 0, 0 }; // used for debug only
long shapeIndex(const TopoDS_Shape& S) long shapeIndex(const TopoDS_Shape& S)
{ {
if ( theMeshDS[0] && theMeshDS[1] ) if ( theMeshDS[0] && theMeshDS[1] )
@ -2402,3 +2403,229 @@ void StdMeshers_ProjectionUtils::SetEventListener(SMESH_subMesh* subMesh,
} }
} }
} }
namespace StdMeshers_ProjectionUtils
{
//================================================================================
/*!
* \brief Computes transformation beween two sets of 2D points using
* a least square approximation
*
* See "Surface Mesh Projection For Hexahedral Mesh Generation By Sweeping"
* by X.Roca, J.Sarrate, A.Huerta. (2.2)
*/
//================================================================================
bool TrsfFinder2D::Solve( const vector< gp_XY >& srcPnts,
const vector< gp_XY >& tgtPnts )
{
// find gravity centers
gp_XY srcGC( 0,0 ), tgtGC( 0,0 );
for ( size_t i = 0; i < srcPnts.size(); ++i )
{
srcGC += srcPnts[i];
tgtGC += tgtPnts[i];
}
srcGC /= srcPnts.size();
tgtGC /= tgtPnts.size();
// find trsf
math_Matrix mat (1,4,1,4, 0.);
math_Vector vec (1,4, 0.);
// cout << "m1 = smesh.Mesh('src')" << endl
// << "m2 = smesh.Mesh('tgt')" << endl;
double xx = 0, xy = 0, yy = 0;
for ( size_t i = 0; i < srcPnts.size(); ++i )
{
gp_XY srcUV = srcPnts[i] - srcGC;
gp_XY tgtUV = tgtPnts[i] - tgtGC;
xx += srcUV.X() * srcUV.X();
yy += srcUV.Y() * srcUV.Y();
xy += srcUV.X() * srcUV.Y();
vec( 1 ) += srcUV.X() * tgtUV.X();
vec( 2 ) += srcUV.Y() * tgtUV.X();
vec( 3 ) += srcUV.X() * tgtUV.Y();
vec( 4 ) += srcUV.Y() * tgtUV.Y();
// cout << "m1.AddNode( " << srcUV.X() << ", " << srcUV.Y() << ", 0 )" << endl
// << "m2.AddNode( " << tgtUV.X() << ", " << tgtUV.Y() << ", 0 )" << endl;
}
mat( 1,1 ) = mat( 3,3 ) = xx;
mat( 2,2 ) = mat( 4,4 ) = yy;
mat( 1,2 ) = mat( 2,1 ) = mat( 3,4 ) = mat( 4,3 ) = xy;
math_Gauss solver( mat );
if ( !solver.IsDone() )
return false;
solver.Solve( vec );
if ( vec.Norm2() < gp::Resolution() )
return false;
// cout << vec( 1 ) << "\t " << vec( 2 ) << endl
// << vec( 3 ) << "\t " << vec( 4 ) << endl;
_trsf.SetTranslation( tgtGC );
_srcOrig = srcGC;
gp_Mat2d& M = const_cast< gp_Mat2d& >( _trsf.HVectorialPart());
M( 1,1 ) = vec( 1 );
M( 2,1 ) = vec( 2 );
M( 1,2 ) = vec( 3 );
M( 2,2 ) = vec( 4 );
return true;
}
//================================================================================
/*!
* \brief Transforms a 2D points using a found transformation
*/
//================================================================================
gp_XY TrsfFinder2D::Transform( const gp_Pnt2d& srcUV ) const
{
gp_XY uv = srcUV.XY() - _srcOrig ;
_trsf.Transforms( uv );
return uv;
}
//================================================================================
/*!
* \brief Computes transformation beween two sets of 3D points using
* a least square approximation
*
* See "Surface Mesh Projection For Hexahedral Mesh Generation By Sweeping"
* by X.Roca, J.Sarrate, A.Huerta. (2.4)
*/
//================================================================================
bool TrsfFinder3D::Solve( const vector< gp_XYZ > & srcPnts,
const vector< gp_XYZ > & tgtPnts )
{
// find gravity center
gp_XYZ srcGC( 0,0,0 ), tgtGC( 0,0,0 );
for ( size_t i = 0; i < srcPnts.size(); ++i )
{
srcGC += srcPnts[i];
tgtGC += tgtPnts[i];
}
srcGC /= srcPnts.size();
tgtGC /= tgtPnts.size();
gp_XYZ srcOrig = 2 * srcGC - tgtGC;
gp_XYZ tgtOrig = srcGC;
// find trsf
math_Matrix mat (1,9,1,9, 0.);
math_Vector vec (1,9, 0.);
double xx = 0, yy = 0, zz = 0;
double xy = 0, xz = 0, yz = 0;
for ( size_t i = 0; i < srcPnts.size(); ++i )
{
gp_XYZ src = srcPnts[i] - srcOrig;
gp_XYZ tgt = tgtPnts[i] - tgtOrig;
xx += src.X() * src.X();
yy += src.Y() * src.Y();
zz += src.Z() * src.Z();
xy += src.X() * src.Y();
xz += src.X() * src.Z();
yz += src.Y() * src.Z();
vec( 1 ) += src.X() * tgt.X();
vec( 2 ) += src.Y() * tgt.X();
vec( 3 ) += src.Z() * tgt.X();
vec( 4 ) += src.X() * tgt.Y();
vec( 5 ) += src.Y() * tgt.Y();
vec( 6 ) += src.Z() * tgt.Y();
vec( 7 ) += src.X() * tgt.Z();
vec( 8 ) += src.Y() * tgt.Z();
vec( 9 ) += src.Z() * tgt.Z();
}
mat( 1,1 ) = mat( 4,4 ) = mat( 7,7 ) = xx;
mat( 2,2 ) = mat( 5,5 ) = mat( 8,8 ) = yy;
mat( 3,3 ) = mat( 6,6 ) = mat( 9,9 ) = zz;
mat( 1,2 ) = mat( 2,1 ) = mat( 4,5 ) = mat( 5,4 ) = mat( 7,8 ) = mat( 8,7 ) = xy;
mat( 1,3 ) = mat( 3,1 ) = mat( 4,6 ) = mat( 6,4 ) = mat( 7,9 ) = mat( 9,7 ) = xz;
mat( 2,3 ) = mat( 3,2 ) = mat( 5,6 ) = mat( 6,5 ) = mat( 8,9 ) = mat( 9,8 ) = yz;
math_Gauss solver( mat );
if ( !solver.IsDone() )
return false;
solver.Solve( vec );
if ( vec.Norm2() < gp::Resolution() )
return false;
// cout << endl
// << vec( 1 ) << "\t " << vec( 2 ) << "\t " << vec( 3 ) << endl
// << vec( 4 ) << "\t " << vec( 5 ) << "\t " << vec( 6 ) << endl
// << vec( 7 ) << "\t " << vec( 8 ) << "\t " << vec( 9 ) << endl;
_srcOrig = srcOrig;
_trsf.SetTranslation( tgtOrig );
gp_Mat& M = const_cast< gp_Mat& >( _trsf.HVectorialPart() );
M.SetRows( gp_XYZ( vec( 1 ), vec( 2 ), vec( 3 )),
gp_XYZ( vec( 4 ), vec( 5 ), vec( 6 )),
gp_XYZ( vec( 7 ), vec( 8 ), vec( 9 )));
return true;
}
//================================================================================
/*!
* \brief Transforms a 3D point using a found transformation
*/
//================================================================================
gp_XYZ TrsfFinder3D::Transform( const gp_Pnt& srcP ) const
{
gp_XYZ p = srcP.XYZ() - _srcOrig;
_trsf.Transforms( p );
return p;
}
//================================================================================
/*!
* \brief Transforms a 3D vector using a found transformation
*/
//================================================================================
gp_XYZ TrsfFinder3D::TransformVec( const gp_Vec& v ) const
{
return v.XYZ().Multiplied( _trsf.HVectorialPart() );
}
//================================================================================
/*!
* \brief Inversion
*/
//================================================================================
bool TrsfFinder3D::Invert()
{
if (( _trsf.Form() == gp_Translation ) &&
( _srcOrig.X() != 0 || _srcOrig.Y() != 0 || _srcOrig.Z() != 0 ))
{
// seems to be defined via Solve()
gp_XYZ newSrcOrig = _trsf.TranslationPart();
gp_Mat& M = const_cast< gp_Mat& >( _trsf.HVectorialPart() );
const double D = M.Determinant();
if ( D < 1e-3 * ( newSrcOrig - _srcOrig ).Modulus() )
{
#ifdef _DEBUG_
cerr << "TrsfFinder3D::Invert()"
<< "D " << M.Determinant() << " IsSingular " << M.IsSingular() << endl;
#endif
return false;
}
gp_Mat Minv = M.Inverted();
_trsf.SetTranslation( _srcOrig );
_srcOrig = newSrcOrig;
M = Minv;
}
else
{
_trsf.Invert();
}
return true;
}
}

55
src/StdMeshers/StdMeshers_ProjectionUtils.hxx
View File

@ -30,10 +30,14 @@
#include "SMESH_StdMeshers.hxx" #include "SMESH_StdMeshers.hxx"
#include "SMDS_MeshElement.hxx"
#include <TopTools_DataMapOfShapeShape.hxx> #include <TopTools_DataMapOfShapeShape.hxx>
#include <TopoDS_Edge.hxx> #include <TopoDS_Edge.hxx>
#include <TopoDS_Vertex.hxx>
#include <TopoDS_Face.hxx> #include <TopoDS_Face.hxx>
#include <TopoDS_Vertex.hxx>
#include <gp_Trsf.hxx>
#include <gp_Trsf2d.hxx>
#include <list> #include <list>
#include <map> #include <map>
@ -77,7 +81,54 @@ namespace StdMeshers_ProjectionUtils
{ {
typedef StdMeshers_ShapeShapeBiDirectionMap TShapeShapeMap; typedef StdMeshers_ShapeShapeBiDirectionMap TShapeShapeMap;
typedef TopTools_IndexedDataMapOfShapeListOfShape TAncestorMap; typedef TopTools_IndexedDataMapOfShapeListOfShape TAncestorMap;
typedef std::map<const SMDS_MeshNode*, const SMDS_MeshNode*> TNodeNodeMap; typedef std::map<const SMDS_MeshNode*, const SMDS_MeshNode*,
TIDCompare> TNodeNodeMap;
/*!
* \brief Finds transformation beween two sets of 2D points using
* a least square approximation
*/
class TrsfFinder2D
{
gp_Trsf2d _trsf;
gp_XY _srcOrig;
public:
TrsfFinder2D(): _srcOrig(0,0) {}
void Set( const gp_Trsf2d& t ) { _trsf = t; } // it's an alternative to Solve()
bool Solve( const std::vector< gp_XY >& srcPnts,
const std::vector< gp_XY >& tgtPnts );
gp_XY Transform( const gp_Pnt2d& srcUV ) const;
bool IsIdentity() const { return ( _trsf.Form() == gp_Identity ); }
};
/*!
* \brief Finds transformation beween two sets of 3D points using
* a least square approximation
*/
class TrsfFinder3D
{
gp_Trsf _trsf;
gp_XYZ _srcOrig;
public:
TrsfFinder3D(): _srcOrig(0,0,0) {}
void Set( const gp_Trsf& t ) { _trsf = t; } // it's an alternative to Solve()
bool Solve( const std::vector< gp_XYZ > & srcPnts,
const std::vector< gp_XYZ > & tgtPnts );
gp_XYZ Transform( const gp_Pnt& srcP ) const;
gp_XYZ TransformVec( const gp_Vec& v ) const;
bool IsIdentity() const { return ( _trsf.Form() == gp_Identity ); }
bool Invert();
};
/*! /*!
* \brief Looks for association of all sub-shapes of two shapes * \brief Looks for association of all sub-shapes of two shapes

525
src/StdMeshers/StdMeshers_Projection_2D.cxx
View File

@ -47,8 +47,10 @@
#include "utilities.h" #include "utilities.h"
#include <BRepAdaptor_Surface.hxx>
#include <BRep_Tool.hxx> #include <BRep_Tool.hxx>
#include <Bnd_B2d.hxx> #include <Bnd_B2d.hxx>
#include <GeomAPI_ProjectPointOnSurf.hxx>
#include <TopExp.hxx> #include <TopExp.hxx>
#include <TopExp_Explorer.hxx> #include <TopExp_Explorer.hxx>
#include <TopTools_DataMapIteratorOfDataMapOfShapeShape.hxx> #include <TopTools_DataMapIteratorOfDataMapOfShapeShape.hxx>
@ -372,125 +374,197 @@ namespace {
} // bool getBoundaryNodes() } // bool getBoundaryNodes()
//================================================================================
/*!
* \brief Check if two consecutive EDGEs are connected in 2D
* \param [in] E1 - a well oriented non-seam EDGE
* \param [in] E2 - a possibly well oriented seam EDGE
* \param [in] F - a FACE
* \return bool - result
*/
//================================================================================
bool are2dConnected( const TopoDS_Edge & E1,
const TopoDS_Edge & E2,
const TopoDS_Face & F )
{
double f,l;
Handle(Geom2d_Curve) c1 = BRep_Tool::CurveOnSurface( E1, F, f, l );
gp_Pnt2d uvLast1 = c1->Value( E1.Orientation() == TopAbs_REVERSED ? f : l );
Handle(Geom2d_Curve) c2 = BRep_Tool::CurveOnSurface( E2, F, f, l );
gp_Pnt2d uvFirst2 = c2->Value( f );
gp_Pnt2d uvLast2 = c2->Value( l );
double tol2 = 1e-5 * uvLast2.SquareDistance( uvFirst2 );
return (( uvLast1.SquareDistance( uvFirst2 ) < tol2 ) ||
( uvLast1.SquareDistance( uvLast2 ) < tol2 ));
}
//================================================================================
/*!
* \brief Compose TSideVector for both FACEs keeping matching order of EDGEs
*/
//================================================================================
bool getWires(const TopoDS_Face& tgtFace,
const TopoDS_Face& srcFace,
SMESH_Mesh * tgtMesh,
SMESH_Mesh * srcMesh,
const TAssocTool::TShapeShapeMap& shape2ShapeMap,
TSideVector& srcWires,
TSideVector& tgtWires,
const bool is1DComputed)
{
// get ordered src EDGEs
TError err;
srcWires = StdMeshers_FaceSide::GetFaceWires( srcFace, *srcMesh,/*skipMediumNodes=*/0, err);
if ( err && !err->IsOK() )
return false;
// make corresponding sequence of tgt EDGEs
tgtWires.resize( srcWires.size() );
for ( size_t iW = 0; iW < srcWires.size(); ++iW )
{
list< TopoDS_Edge > tgtEdges;
StdMeshers_FaceSidePtr srcWire = srcWires[iW];
TopTools_IndexedMapOfShape edgeMap; // to detect seam edges
for ( int iE = 0; iE < srcWire->NbEdges(); ++iE )
{
TopoDS_Edge E = TopoDS::Edge( shape2ShapeMap( srcWire->Edge( iE ), /*isSrc=*/true));
// reverse a seam edge encountered for the second time
const int index = edgeMap.Add( E );
if ( index < edgeMap.Extent() ) // E is a seam
{
// check which of edges to reverse, E or one already being in tgtEdges
if ( are2dConnected( tgtEdges.back(), E, tgtFace ))
{
list< TopoDS_Edge >::iterator eIt = tgtEdges.begin();
std::advance( eIt, index-1 );
eIt->Reverse();
}
else
{
E.Reverse();
}
}
if ( srcWire->NbEdges() == 1 && tgtMesh == srcMesh ) // circle
{
// try to verify ori by propagation
pair<int,TopoDS_Edge> nE = StdMeshers_ProjectionUtils::GetPropagationEdge
( srcMesh, E, srcWire->Edge( iE ));
if ( !nE.second.IsNull() )
E = nE.second;
}
tgtEdges.push_back( E );
}
tgtWires[ iW ].reset( new StdMeshers_FaceSide( tgtFace, tgtEdges, tgtMesh,
/*theIsForward = */ true,
/*theIgnoreMediumNodes = */false));
if ( is1DComputed &&
srcWires[iW]->GetUVPtStruct().size() !=
tgtWires[iW]->GetUVPtStruct().size())
return false;
}
return true;
}
//================================================================================ //================================================================================
/*! /*!
* \brief Preform projection in case if tgtFace.IsPartner( srcFace ) and in case * \brief Preform projection in case if tgtFace.IsPartner( srcFace ) and in case
* if projection by transformation is possible * if projection by 3D transformation is possible
*/ */
//================================================================================ //================================================================================
bool projectPartner(const TopoDS_Face& tgtFace, bool projectPartner(const TopoDS_Face& tgtFace,
const TopoDS_Face& srcFace, const TopoDS_Face& srcFace,
SMESH_Mesh * tgtMesh, const TSideVector& tgtWires,
SMESH_Mesh * srcMesh, const TSideVector& srcWires,
const TAssocTool::TShapeShapeMap& shape2ShapeMap) const TAssocTool::TShapeShapeMap& shape2ShapeMap,
TAssocTool::TNodeNodeMap& src2tgtNodes)
{ {
SMESH_Mesh * tgtMesh = tgtWires[0]->GetMesh();
SMESH_Mesh * srcMesh = srcWires[0]->GetMesh();
SMESHDS_Mesh* tgtMeshDS = tgtMesh->GetMeshDS(); SMESHDS_Mesh* tgtMeshDS = tgtMesh->GetMeshDS();
SMESHDS_Mesh* srcMeshDS = srcMesh->GetMeshDS(); SMESHDS_Mesh* srcMeshDS = srcMesh->GetMeshDS();
SMESH_MesherHelper helper( *tgtMesh );
const double tol = 1.e-7 * srcMeshDS->getMaxDim(); const double tol = 1.e-7 * srcMeshDS->getMaxDim();
gp_Trsf trsf; // transformation to get location of target nodes from source ones // transformation to get location of target nodes from source ones
StdMeshers_ProjectionUtils::TrsfFinder3D trsf;
if ( tgtFace.IsPartner( srcFace )) if ( tgtFace.IsPartner( srcFace ))
{ {
gp_Trsf srcTrsf = srcFace.Location(); gp_Trsf srcTrsf = srcFace.Location();
gp_Trsf tgtTrsf = tgtFace.Location(); gp_Trsf tgtTrsf = tgtFace.Location();
trsf = srcTrsf.Inverted() * tgtTrsf; trsf.Set( srcTrsf.Inverted() * tgtTrsf );
} }
else else
{ {
// Try to find the transformation // Try to find the 3D transformation
// make any local coord systems of src and tgt faces const int totNbSeg = 50;
vector<gp_Pnt> srcPP, tgtPP; // 3 points on face boundaries to make axes of CS vector< gp_XYZ > srcPnts, tgtPnts;
int tgtNbVert = SMESH_MesherHelper::Count( tgtFace, TopAbs_VERTEX, /*ignoreSame=*/true ); srcPnts.reserve( totNbSeg );
int srcNbVert = SMESH_MesherHelper::Count( srcFace, TopAbs_VERTEX, /*ignoreSame=*/true ); tgtPnts.reserve( totNbSeg );
SMESH_subMesh * srcSM = srcMesh->GetSubMesh( srcFace ); for ( size_t iW = 0; iW < srcWires.size(); ++iW )
SMESH_subMeshIteratorPtr smIt = srcSM->getDependsOnIterator(/*includeSelf=*/false,false);
srcSM = smIt->next(); // sm of a vertex
while ( smIt->more() && srcPP.size() < 3 )
{ {
srcSM = smIt->next(); const double minSegLen = srcWires[iW]->Length() / totNbSeg;
SMESHDS_SubMesh* srcSmds = srcSM->GetSubMeshDS(); for ( int iE = 0; iE < srcWires[iW]->NbEdges(); ++iE )
if ( !srcSmds ) continue;
SMDS_NodeIteratorPtr nIt = srcSmds->GetNodes();
while ( nIt->more() )
{ {
SMESH_TNodeXYZ p ( nIt->next()); int nbSeg = Max( 1, int( srcWires[iW]->EdgeLength( iE ) / minSegLen ));
bool pOK = false; double srcU = srcWires[iW]->FirstParameter( iE );
switch ( srcPP.size() ) double tgtU = tgtWires[iW]->FirstParameter( iE );
double srcDu = ( srcWires[iW]->LastParameter( iE )- srcU ) / nbSeg;
double tgtDu = ( tgtWires[iW]->LastParameter( iE )- tgtU ) / nbSeg;
for ( size_t i = 0; i < nbSeg; ++i )
{ {
case 0: pOK = true; break; srcPnts.push_back( srcWires[iW]->Value3d( srcU ).XYZ() );
tgtPnts.push_back( tgtWires[iW]->Value3d( tgtU ).XYZ() );
case 1: pOK = ( srcPP[0].SquareDistance( p ) > 10*tol ); break; srcU += srcDu;
tgtU += tgtDu;
case 2:
{
gp_Vec p0p1( srcPP[0], srcPP[1] ), p0p( srcPP[0], p );
// pOK = !p0p1.IsParallel( p0p, tol );
pOK = !p0p1.IsParallel( p0p, 3.14/20 ); // angle min 18 degrees
break;
}
}
if ( !pOK )
continue;
// find corresponding point on target shape
pOK = false;
gp_Pnt tgtP;
const TopoDS_Shape& tgtShape = shape2ShapeMap( srcSM->GetSubShape(), /*isSrc=*/true );
if ( tgtShape.ShapeType() == TopAbs_VERTEX )
{
tgtP = BRep_Tool::Pnt( TopoDS::Vertex( tgtShape ));
if ( srcNbVert == tgtNbVert || tgtPP.empty() )
pOK = true;
else
pOK = (( tgtP.Distance( tgtPP[0] ) > tol*tol ) &&
( tgtPP.size() == 1 || tgtP.Distance( tgtPP[1] ) > tol*tol ));
//cout << "V - nS " << p._node->GetID() << " - nT " << SMESH_Algo::VertexNode(TopoDS::Vertex( tgtShape),tgtMeshDS)->GetID() << endl;
}
else if ( tgtPP.size() > 0 )
{
if ( SMESHDS_SubMesh* tgtSmds = tgtMeshDS->MeshElements( tgtShape ))
{
double srcDist = srcPP[0].Distance( p );
double eTol = BRep_Tool::Tolerance( TopoDS::Edge( tgtShape ));
if (eTol < tol) eTol = tol;
SMDS_NodeIteratorPtr nItT = tgtSmds->GetNodes();
while ( nItT->more() && !pOK )
{
const SMDS_MeshNode* n = nItT->next();
tgtP = SMESH_TNodeXYZ( n );
pOK = ( fabs( srcDist - tgtPP[0].Distance( tgtP )) < 2*eTol );
//cout << "E - nS " << p._node->GetID() << " - nT " << n->GetID()<< " OK - " << pOK<< " " << fabs( srcDist - tgtPP[0].Distance( tgtP ))<< " tol " << eTol<< endl;
} }
} }
} }
if ( !pOK ) if ( !trsf.Solve( srcPnts, tgtPnts ))
continue;
srcPP.push_back( p );
tgtPP.push_back( tgtP );
}
}
if ( srcPP.size() != 3 )
return false; return false;
// make transformation // check trsf
gp_Trsf fromTgtCS, toSrcCS; // from/to global CS
gp_Ax2 srcCS( srcPP[0], gp_Vec( srcPP[0], srcPP[1] ), gp_Vec( srcPP[0], srcPP[2]));
gp_Ax2 tgtCS( tgtPP[0], gp_Vec( tgtPP[0], tgtPP[1] ), gp_Vec( tgtPP[0], tgtPP[2]));
toSrcCS .SetTransformation( gp_Ax3( srcCS ));
fromTgtCS.SetTransformation( gp_Ax3( tgtCS ));
fromTgtCS.Invert();
trsf = fromTgtCS * toSrcCS; bool trsfIsOK = true;
const int nbTestPnt = 20;
const size_t iStep = Max( 1, int( srcPnts.size() / nbTestPnt ));
// check boundary
for ( size_t i = 0; ( i < srcPnts.size() && trsfIsOK ); i += iStep )
{
gp_Pnt trsfTgt = trsf.Transform( srcPnts[i] );
trsfIsOK = ( trsfTgt.SquareDistance( tgtPnts[i] ) < tol*tol );
}
// check an in-FACE point
if ( trsfIsOK )
{
BRepAdaptor_Surface srcSurf( srcFace );
gp_Pnt srcP =
srcSurf.Value( 0.5 * ( srcSurf.FirstUParameter() + srcSurf.LastUParameter() ),
0.5 * ( srcSurf.FirstVParameter() + srcSurf.LastVParameter() ));
gp_Pnt tgtTrsfP = trsf.Transform( srcP );
TopLoc_Location loc;
GeomAPI_ProjectPointOnSurf& proj = helper.GetProjector( tgtFace, loc, 0.1*tol );
if ( !loc.IsIdentity() )
tgtTrsfP.Transform( loc.Transformation().Inverted() );
proj.Perform( tgtTrsfP );
trsfIsOK = ( proj.IsDone() &&
proj.NbPoints() > 0 &&
proj.LowerDistance() < tol );
}
if ( !trsfIsOK )
return false;
} }
// Fill map of src to tgt nodes with nodes on edges // Fill map of src to tgt nodes with nodes on edges
map<const SMDS_MeshNode* , const SMDS_MeshNode*> src2tgtNodes; src2tgtNodes.clear();
map<const SMDS_MeshNode* , const SMDS_MeshNode*>::iterator srcN_tgtN; TAssocTool::TNodeNodeMap::iterator srcN_tgtN;
bool tgtEdgesMeshed = false; bool tgtEdgesMeshed = false;
for ( TopExp_Explorer srcExp( srcFace, TopAbs_EDGE); srcExp.More(); srcExp.Next() ) for ( TopExp_Explorer srcExp( srcFace, TopAbs_EDGE); srcExp.More(); srcExp.Next() )
@ -518,31 +592,6 @@ namespace {
) )
return false; return false;
if ( !tgtEdge.IsPartner( srcEdge ))
{
if ( tgtNodes.empty() )
return false;
// check that transformation is OK by three nodes
gp_Pnt p0S = SMESH_TNodeXYZ( (srcNodes.begin()) ->second);
gp_Pnt p1S = SMESH_TNodeXYZ( (srcNodes.rbegin()) ->second);
gp_Pnt p2S = SMESH_TNodeXYZ( (++srcNodes.begin())->second);
gp_Pnt p0T = SMESH_TNodeXYZ( (tgtNodes.begin()) ->second);
gp_Pnt p1T = SMESH_TNodeXYZ( (tgtNodes.rbegin()) ->second);
gp_Pnt p2T = SMESH_TNodeXYZ( (++tgtNodes.begin())->second);
// transform source points, they must coincide with target ones
if ( p0T.SquareDistance( p0S.Transformed( trsf )) > tol ||
p1T.SquareDistance( p1S.Transformed( trsf )) > tol ||
p2T.SquareDistance( p2S.Transformed( trsf )) > tol )
{
//cout << "KO trsf, 3 dist: "
//<< p0T.SquareDistance( p0S.Transformed( trsf ))<< ", "
//<< p1T.SquareDistance( p1S.Transformed( trsf ))<< ", "
//<< p2T.SquareDistance( p2S.Transformed( trsf ))<< ", "<<endl;
return false;
}
}
if ( !tgtNodes.empty() ) if ( !tgtNodes.empty() )
{ {
map< double, const SMDS_MeshNode* >::iterator u_tn = tgtNodes.begin(); map< double, const SMDS_MeshNode* >::iterator u_tn = tgtNodes.begin();
@ -563,16 +612,6 @@ namespace {
if ( !tgtN || tgtV.ShapeType() != TopAbs_VERTEX ) if ( !tgtN || tgtV.ShapeType() != TopAbs_VERTEX )
return false; return false;
if ( !tgtV.IsPartner( srcV ))
{
// check that transformation is OK by three nodes
gp_Pnt p0S = SMESH_TNodeXYZ( srcN );
gp_Pnt p0T = SMESH_TNodeXYZ( tgtN );
if ( p0T.SquareDistance( p0S.Transformed( trsf )) > tol )
{
return false;
}
}
src2tgtNodes.insert( make_pair( srcN, tgtN )); src2tgtNodes.insert( make_pair( srcN, tgtN ));
} }
@ -580,7 +619,6 @@ namespace {
// Make new faces // Make new faces
// prepare the helper to adding quadratic elements if necessary // prepare the helper to adding quadratic elements if necessary
SMESH_MesherHelper helper( *tgtMesh );
helper.SetSubShape( tgtFace ); helper.SetSubShape( tgtFace );
helper.IsQuadraticSubMesh( tgtFace ); helper.IsQuadraticSubMesh( tgtFace );
@ -594,7 +632,7 @@ namespace {
const SMDS_MeshNode* nullNode = 0; const SMDS_MeshNode* nullNode = 0;
// indices of nodes to create properly oriented faces // indices of nodes to create properly oriented faces
bool isReverse = ( trsf.Form() != gp_Identity ); bool isReverse = ( !trsf.IsIdentity() );
int tri1 = 1, tri2 = 2, quad1 = 1, quad3 = 3; int tri1 = 1, tri2 = 2, quad1 = 1, quad3 = 3;
if ( isReverse ) if ( isReverse )
std::swap( tri1, tri2 ), std::swap( quad1, quad3 ); std::swap( tri1, tri2 ), std::swap( quad1, quad3 );
@ -614,7 +652,7 @@ namespace {
if ( srcN_tgtN->second == nullNode ) if ( srcN_tgtN->second == nullNode )
{ {
// create a new node // create a new node
gp_Pnt tgtP = gp_Pnt( SMESH_TNodeXYZ( srcNode )).Transformed( trsf ); gp_Pnt tgtP = trsf.Transform( SMESH_TNodeXYZ( srcNode ));
SMDS_MeshNode* n = helper.AddNode( tgtP.X(), tgtP.Y(), tgtP.Z() ); SMDS_MeshNode* n = helper.AddNode( tgtP.X(), tgtP.Y(), tgtP.Z() );
srcN_tgtN->second = n; srcN_tgtN->second = n;
switch ( srcNode->GetPosition()->GetTypeOfPosition() ) switch ( srcNode->GetPosition()->GetTypeOfPosition() )
@ -698,33 +736,6 @@ namespace {
} // bool projectPartner() } // bool projectPartner()
//================================================================================
/*!
* \brief Check if two consecutive EDGEs are connected in 2D
* \param [in] E1 - a well oriented non-seam EDGE
* \param [in] E2 - a possibly well oriented seam EDGE
* \param [in] F - a FACE
* \return bool - result
*/
//================================================================================
bool are2dConnected( const TopoDS_Edge & E1,
const TopoDS_Edge & E2,
const TopoDS_Face & F )
{
double f,l;
Handle(Geom2d_Curve) c1 = BRep_Tool::CurveOnSurface( E1, F, f, l );
gp_Pnt2d uvLast1 = c1->Value( E1.Orientation() == TopAbs_REVERSED ? f : l );
Handle(Geom2d_Curve) c2 = BRep_Tool::CurveOnSurface( E2, F, f, l );
gp_Pnt2d uvFirst2 = c2->Value( f );
gp_Pnt2d uvLast2 = c2->Value( l );
double tol2 = 1e-5 * uvLast2.SquareDistance( uvFirst2 );
return (( uvLast1.SquareDistance( uvFirst2 ) < tol2 ) ||
( uvLast1.SquareDistance( uvLast2 ) < tol2 ));
}
//================================================================================ //================================================================================
/*! /*!
* \brief Preform projection in case if the faces are similar in 2D space * \brief Preform projection in case if the faces are similar in 2D space
@ -733,60 +744,21 @@ namespace {
bool projectBy2DSimilarity(const TopoDS_Face& tgtFace, bool projectBy2DSimilarity(const TopoDS_Face& tgtFace,
const TopoDS_Face& srcFace, const TopoDS_Face& srcFace,
SMESH_Mesh * tgtMesh, const TSideVector& tgtWires,
SMESH_Mesh * srcMesh, const TSideVector& srcWires,
const TAssocTool::TShapeShapeMap& shape2ShapeMap, const TAssocTool::TShapeShapeMap& shape2ShapeMap,
const bool is1DComputed) const bool is1DComputed,
TAssocTool::TNodeNodeMap& src2tgtNodes)
{ {
// 1) Preparation SMESH_Mesh * tgtMesh = tgtWires[0]->GetMesh();
SMESH_Mesh * srcMesh = srcWires[0]->GetMesh();
// get ordered src EDGEs // WARNING: we can have problems if the FACE is symmetrical in 2D,
TError err; // then the projection can be mirrored relating to what is expected
TSideVector srcWires =
StdMeshers_FaceSide::GetFaceWires( srcFace, *srcMesh,/*ignoreMediumNodes = */false, err);
if ( err && !err->IsOK() )
return false;
// make corresponding sequence of tgt EDGEs // 1) Find 2D transformation
TSideVector tgtWires( srcWires.size() );
for ( size_t iW = 0; iW < srcWires.size(); ++iW )
{
list< TopoDS_Edge > tgtEdges;
StdMeshers_FaceSidePtr srcWire = srcWires[iW];
TopTools_IndexedMapOfShape edgeMap; // to detect seam edges
for ( int iE = 0; iE < srcWire->NbEdges(); ++iE )
{
TopoDS_Edge E = TopoDS::Edge( shape2ShapeMap( srcWire->Edge( iE ), /*isSrc=*/true));
// reverse a seam edge encountered for the second time
const int index = edgeMap.Add( E );
if ( index < edgeMap.Extent() ) // E is a seam
{
// check which of edges to reverse, E or one already being in tgtEdges
if ( are2dConnected( tgtEdges.back(), E, tgtFace ))
{
list< TopoDS_Edge >::iterator eIt = tgtEdges.begin();
std::advance( eIt, index-1 );
eIt->Reverse();
}
else
{
E.Reverse();
}
}
tgtEdges.push_back( E );
}
tgtWires[ iW ].reset( new StdMeshers_FaceSide( tgtFace, tgtEdges, tgtMesh,
/*theIsForward = */ true,
/*theIgnoreMediumNodes = */false));
if ( is1DComputed &&
srcWires[iW]->GetUVPtStruct().size() !=
tgtWires[iW]->GetUVPtStruct().size())
return false;
}
// 2) Find transformation StdMeshers_ProjectionUtils::TrsfFinder2D trsf;
gp_Trsf2d trsf;
{ {
// get 2 pairs of corresponding UVs // get 2 pairs of corresponding UVs
gp_Pnt2d srcP0 = srcWires[0]->Value2d(0.0); gp_Pnt2d srcP0 = srcWires[0]->Value2d(0.0);
@ -801,34 +773,74 @@ namespace {
toSrcCS .SetTransformation( srcCS ); toSrcCS .SetTransformation( srcCS );
fromTgtCS.SetTransformation( tgtCS ); fromTgtCS.SetTransformation( tgtCS );
fromTgtCS.Invert(); fromTgtCS.Invert();
trsf.Set( fromTgtCS * toSrcCS );
trsf = fromTgtCS * toSrcCS;
// check transformation // check transformation
bool trsfIsOK = true;
const double tol = 1e-5 * gp_Vec2d( srcP0, srcP1 ).Magnitude(); const double tol = 1e-5 * gp_Vec2d( srcP0, srcP1 ).Magnitude();
for ( double u = 0.12; u < 1.; u += 0.1 ) for ( double u = 0.12; ( u < 1. && trsfIsOK ); u += 0.1 )
{ {
gp_Pnt2d srcUV = srcWires[0]->Value2d( u ); gp_Pnt2d srcUV = srcWires[0]->Value2d( u );
gp_Pnt2d tgtUV = tgtWires[0]->Value2d( u ); gp_Pnt2d tgtUV = tgtWires[0]->Value2d( u );
gp_Pnt2d tgtUV2 = srcUV.Transformed( trsf ); gp_Pnt2d tgtUV2 = trsf.Transform( srcUV );
if ( tgtUV.Distance( tgtUV2 ) > tol ) trsfIsOK = ( tgtUV.Distance( tgtUV2 ) < tol );
}
// Find trsf using a least-square approximation
if ( !trsfIsOK )
{
// find trsf
const int totNbSeg = 50;
vector< gp_XY > srcPnts, tgtPnts;
srcPnts.resize( totNbSeg );
tgtPnts.resize( totNbSeg );
for ( size_t iW = 0; iW < srcWires.size(); ++iW )
{
const double minSegLen = srcWires[iW]->Length() / totNbSeg;
for ( int iE = 0; iE < srcWires[iW]->NbEdges(); ++iE )
{
int nbSeg = Max( 1, int( srcWires[iW]->EdgeLength( iE ) / minSegLen ));
double srcU = srcWires[iW]->FirstParameter( iE );
double tgtU = tgtWires[iW]->FirstParameter( iE );
double srcDu = ( srcWires[iW]->LastParameter( iE )- srcU ) / nbSeg;
double tgtDu = ( tgtWires[iW]->LastParameter( iE )- tgtU ) / nbSeg;
for ( size_t i = 0; i < nbSeg; ++i, srcU += srcDu, tgtU += tgtDu )
{
srcPnts.push_back( srcWires[iW]->Value2d( srcU ).XY() );
tgtPnts.push_back( tgtWires[iW]->Value2d( tgtU ).XY() );
}
}
}
if ( !trsf.Solve( srcPnts, tgtPnts ))
return false;
// check trsf
trsfIsOK = true;
const int nbTestPnt = 10;
const size_t iStep = Max( 1, int( srcPnts.size() / nbTestPnt ));
for ( size_t i = 0; ( i < srcPnts.size() && trsfIsOK ); i += iStep )
{
gp_Pnt2d trsfTgt = trsf.Transform( srcPnts[i] );
trsfIsOK = ( trsfTgt.Distance( tgtPnts[i] ) < tol );
}
if ( !trsfIsOK )
return false; return false;
} }
} } // "Find transformation" block
// 3) Projection // 2) Projection
typedef map<const SMDS_MeshNode* , const SMDS_MeshNode*, TIDCompare> TN2NMap; src2tgtNodes.clear();
TN2NMap src2tgtNodes; TAssocTool::TNodeNodeMap::iterator srcN_tgtN;
TN2NMap::iterator srcN_tgtN;
// fill src2tgtNodes in with nodes on EDGEs // fill src2tgtNodes in with nodes on EDGEs
for ( unsigned iW = 0; iW < srcWires.size(); ++iW ) for ( size_t iW = 0; iW < srcWires.size(); ++iW )
if ( is1DComputed ) if ( is1DComputed )
{ {
const vector<UVPtStruct>& srcUVs = srcWires[iW]->GetUVPtStruct(); const vector<UVPtStruct>& srcUVs = srcWires[iW]->GetUVPtStruct();
const vector<UVPtStruct>& tgtUVs = tgtWires[iW]->GetUVPtStruct(); const vector<UVPtStruct>& tgtUVs = tgtWires[iW]->GetUVPtStruct();
for ( unsigned i = 0; i < srcUVs.size(); ++i ) for ( size_t i = 0; i < srcUVs.size(); ++i )
src2tgtNodes.insert( make_pair( srcUVs[i].node, tgtUVs[i].node )); src2tgtNodes.insert( make_pair( srcUVs[i].node, tgtUVs[i].node ));
} }
else else
@ -880,7 +892,7 @@ namespace {
// create a new node // create a new node
gp_Pnt2d srcUV = srcHelper.GetNodeUV( srcFace, srcNode, gp_Pnt2d srcUV = srcHelper.GetNodeUV( srcFace, srcNode,
elem->GetNode( helper.WrapIndex(i+1,nbN)), &uvOK); elem->GetNode( helper.WrapIndex(i+1,nbN)), &uvOK);
gp_Pnt2d tgtUV = srcUV.Transformed( trsf ); gp_Pnt2d tgtUV = trsf.Transform( srcUV );
gp_Pnt tgtP = tgtSurface->Value( tgtUV.X(), tgtUV.Y() ); gp_Pnt tgtP = tgtSurface->Value( tgtUV.X(), tgtUV.Y() );
SMDS_MeshNode* n = tgtMeshDS->AddNode( tgtP.X(), tgtP.Y(), tgtP.Z() ); SMDS_MeshNode* n = tgtMeshDS->AddNode( tgtP.X(), tgtP.Y(), tgtP.Z() );
switch ( srcNode->GetPosition()->GetTypeOfPosition() ) switch ( srcNode->GetPosition()->GetTypeOfPosition() )
@ -920,6 +932,74 @@ namespace {
} // bool projectBy2DSimilarity(...) } // bool projectBy2DSimilarity(...)
//================================================================================
/*!
* \brief Fix bad faces by smoothing
*/
//================================================================================
void fixDistortedFaces( SMESH_MesherHelper& helper )
{
// Detect bad faces
bool haveBadFaces = false;
const TopoDS_Face& F = TopoDS::Face( helper.GetSubShape() );
SMESHDS_SubMesh* smDS = helper.GetMeshDS()->MeshElements( F );
if ( !smDS || smDS->NbElements() == 0 ) return;
SMDS_ElemIteratorPtr faceIt = smDS->GetElements();
double prevArea2D = 0;
vector< const SMDS_MeshNode* > nodes;
vector< gp_XY > uv;
while ( faceIt->more() && !haveBadFaces )
{
const SMDS_MeshElement* face = faceIt->next();
// get nodes
nodes.resize( face->NbCornerNodes() );
SMDS_MeshElement::iterator n = face->begin_nodes();
for ( size_t i = 0; i < nodes.size(); ++n, ++i )
nodes[ i ] = *n;
// get UVs
const SMDS_MeshNode* inFaceNode = 0;
if ( helper.HasSeam() )
for ( size_t i = 0; ( i < nodes.size() && !inFaceNode ); ++i )
if ( !helper.IsSeamShape( nodes[ i ]->getshapeId() ))
inFaceNode = nodes[ i ];
uv.resize( nodes.size() );
for ( size_t i = 0; i < nodes.size(); ++i )
uv[ i ] = helper.GetNodeUV( F, nodes[ i ], inFaceNode );
// compare orientation of triangles
for ( int iT = 0, nbT = nodes.size()-2; iT < nbT; ++iT )
{
gp_XY v1 = uv[ iT+1 ] - uv[ 0 ];
gp_XY v2 = uv[ iT+2 ] - uv[ 0 ];
double area2D = v2 ^ v1;
if (( haveBadFaces = ( area2D * prevArea2D < 0 )))
break;
prevArea2D = area2D;
}
}
// Fix faces
if ( haveBadFaces )
{
SMESH_MeshEditor editor( helper.GetMesh() );
TIDSortedElemSet faces;
for ( faceIt = smDS->GetElements(); faceIt->more(); )
faces.insert( faces.end(), faceIt->next() );
set<const SMDS_MeshNode*> fixedNodes;
editor.Smooth( faces, fixedNodes, SMESH_MeshEditor::CENTROIDAL, 5 );
}
}
} // namespace } // namespace
@ -930,6 +1010,8 @@ namespace {
bool StdMeshers_Projection_2D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theShape) bool StdMeshers_Projection_2D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theShape)
{ {
_src2tgtNodes.clear();
MESSAGE("Projection_2D Compute"); MESSAGE("Projection_2D Compute");
if ( !_sourceHypo ) if ( !_sourceHypo )
return false; return false;
@ -1001,17 +1083,25 @@ bool StdMeshers_Projection_2D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape&
is1DComputed = sm->IsMeshComputed(); is1DComputed = sm->IsMeshComputed();
} }
// get ordered src and tgt EDGEs
TSideVector srcWires, tgtWires;
if ( !getWires( tgtFace, srcFace, tgtMesh, srcMesh,
shape2ShapeMap, srcWires, tgtWires, is1DComputed ))
return false;
bool done = false; bool done = false;
if ( !done ) if ( !done )
{ {
// try to project from the same face with different location // try to project from the same face with different location
done = projectPartner( tgtFace, srcFace, tgtMesh, srcMesh, shape2ShapeMap ); done = projectPartner( tgtFace, srcFace, tgtWires, srcWires,
shape2ShapeMap, _src2tgtNodes );
} }
if ( !done ) if ( !done )
{ {
// projection in case if the faces are similar in 2D space // projection in case if the faces are similar in 2D space
done = projectBy2DSimilarity( tgtFace, srcFace, tgtMesh, srcMesh, shape2ShapeMap, is1DComputed); done = projectBy2DSimilarity( tgtFace, srcFace, tgtWires, srcWires,
shape2ShapeMap, is1DComputed, _src2tgtNodes);
} }
SMESH_MesherHelper helper( theMesh ); SMESH_MesherHelper helper( theMesh );
@ -1019,6 +1109,7 @@ bool StdMeshers_Projection_2D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape&
if ( !done ) if ( !done )
{ {
_src2tgtNodes.clear();
// -------------------- // --------------------
// Prepare to mapping // Prepare to mapping
// -------------------- // --------------------
@ -1278,6 +1369,14 @@ bool StdMeshers_Projection_2D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape&
if ( nbFaceBeforeMerge != nbFaceAtferMerge && !helper.HasDegeneratedEdges() ) if ( nbFaceBeforeMerge != nbFaceAtferMerge && !helper.HasDegeneratedEdges() )
return error(COMPERR_BAD_INPUT_MESH, "Probably invalid node parameters on geom faces"); return error(COMPERR_BAD_INPUT_MESH, "Probably invalid node parameters on geom faces");
// ----------------------------------------------------------------
// The mapper can create distorted faces by placing nodes out of the FACE
// boundary -- fix bad faces by smoothing
// ----------------------------------------------------------------
fixDistortedFaces( helper );
// ---------------------------------------------------------------- // ----------------------------------------------------------------
// The mapper can't create quadratic elements, so convert if needed // The mapper can't create quadratic elements, so convert if needed
// ---------------------------------------------------------------- // ----------------------------------------------------------------

6
src/StdMeshers/StdMeshers_Projection_2D.hxx
View File

@ -30,6 +30,7 @@
#include "SMESH_StdMeshers.hxx" #include "SMESH_StdMeshers.hxx"
#include "SMESH_Algo.hxx" #include "SMESH_Algo.hxx"
#include "StdMeshers_ProjectionUtils.hxx"
class StdMeshers_ProjectionSource2D; class StdMeshers_ProjectionSource2D;
@ -59,10 +60,13 @@ public:
*/ */
virtual void SetEventListener(SMESH_subMesh* whenSetToSubMesh); virtual void SetEventListener(SMESH_subMesh* whenSetToSubMesh);
protected:
protected:
const StdMeshers_ProjectionSource2D* _sourceHypo; const StdMeshers_ProjectionSource2D* _sourceHypo;
StdMeshers_ProjectionUtils::TNodeNodeMap _src2tgtNodes;
}; };
#endif #endif

2
src/StdMeshers/StdMeshers_Projection_3D.cxx
View File

@ -307,7 +307,7 @@ bool StdMeshers_Projection_3D::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aS
} }
} }
// Find matching nodes of tgt and src faces // Find matching nodes of tgt and src faces
TNodeNodeMap faceMatchingNodes; TAssocTool::TNodeNodeMap faceMatchingNodes;
if ( ! TAssocTool::FindMatchingNodesOnFaces( srcFace, srcMesh, tgtFace, tgtMesh, if ( ! TAssocTool::FindMatchingNodesOnFaces( srcFace, srcMesh, tgtFace, tgtMesh,
shape2ShapeMap, faceMatchingNodes )) shape2ShapeMap, faceMatchingNodes ))
return error(COMPERR_BAD_INPUT_MESH,SMESH_Comment("Mesh on faces #") return error(COMPERR_BAD_INPUT_MESH,SMESH_Comment("Mesh on faces #")

2
src/StdMeshers/StdMeshers_RadialPrism_3D.cxx
View File

@ -239,7 +239,7 @@ bool StdMeshers_RadialPrism_3D::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& a
} }
// Find matching nodes of in and out faces // Find matching nodes of in and out faces
TNodeNodeMap nodeIn2OutMap; ProjectionUtils::TNodeNodeMap nodeIn2OutMap;
if ( ! ProjectionUtils::FindMatchingNodesOnFaces( inFace, &aMesh, outFace, &aMesh, if ( ! ProjectionUtils::FindMatchingNodesOnFaces( inFace, &aMesh, outFace, &aMesh,
shape2ShapeMap, nodeIn2OutMap )) shape2ShapeMap, nodeIn2OutMap ))
return error(COMPERR_BAD_INPUT_MESH,SMESH_Comment("Mesh on faces #") return error(COMPERR_BAD_INPUT_MESH,SMESH_Comment("Mesh on faces #")

2
src/StdMeshers/StdMeshers_ViscousLayers.cxx
View File

@ -101,7 +101,7 @@ namespace VISCOUS_3D
const double theSmoothThickToElemSizeRatio = 0.3; const double theSmoothThickToElemSizeRatio = 0.3;
// what part of thickness is allowed till intersection // what part of thickness is allowed till intersection
// defined by SALOME_TESTS/Grids/smesh/viscous_layers_00/A5 // (defined by SALOME_TESTS/Grids/smesh/viscous_layers_00/A5)
const double theThickToIntersection = 1.5; const double theThickToIntersection = 1.5;
bool needSmoothing( double cosin, double tgtThick, double elemSize ) bool needSmoothing( double cosin, double tgtThick, double elemSize )