[bos #42217][EDF 28921] Horseshoe with bodyfitting.

Intermediate commit - parallel compute is disabled, some unrelated code commented out for debug!
Hexahedron::compute() was splitted to smaller methods.
Added debug output.
SMDS_MeshElement::Print() now is a virtual method, otherwise only parent method was called.
This commit is contained in:
Konstantin Leontev 2024-07-24 18:06:11 +01:00
parent 4b39d980ce
commit 9f7018e42a
6 changed files with 687 additions and 450 deletions

View File

@ -32,7 +32,7 @@
class SMDS_EXPORT SMDS_FaceOfNodes: public SMDS_CellOfNodes
{
public:
void Print(std::ostream & OS) const;
virtual void Print(std::ostream & OS) const override;
SMDS_FaceOfNodes(const SMDS_MeshNode* node1,
const SMDS_MeshNode* node2,
const SMDS_MeshNode* node3);

View File

@ -142,7 +142,7 @@ public:
SMDS_Mesh* GetMesh() const;
void Print(std::ostream & OS) const;
virtual void Print(std::ostream & OS) const;
friend SMDS_EXPORT std::ostream & operator <<(std::ostream & OS, const SMDS_MeshElement *);
friend class SMDS_ElementFactory;

View File

@ -65,7 +65,7 @@ class SMDS_EXPORT SMDS_MeshNode: public SMDS_MeshElement
virtual bool IsMediumNode(const SMDS_MeshNode* /*node*/) const { return false; }
virtual int NbCornerNodes() const { return 1; }
void Print(std::ostream & OS) const;
virtual void Print(std::ostream & OS) const override;
private:

View File

@ -47,7 +47,7 @@ class SMDS_EXPORT SMDS_PolygonalFaceOfNodes : public SMDS_CellOfNodes
virtual int NbEdges() const;
virtual int NbFaces() const;
virtual void Print (std::ostream & OS) const;
virtual void Print (std::ostream & OS) const override;
virtual const SMDS_MeshNode* GetNode(const int ind) const;

View File

@ -62,7 +62,7 @@ class SMDS_EXPORT SMDS_VolumeOfNodes: public SMDS_CellOfNodes
const int nbNodes);
~SMDS_VolumeOfNodes();
void Print(std::ostream & OS) const;
virtual void Print(std::ostream & OS) const override;
int NbFaces() const;
int NbNodes() const;
int NbEdges() const;

View File

@ -773,6 +773,16 @@ namespace
if ( node )
_node = _intPoint->_node = node;
}
friend std::ostream& operator<<(std::ostream& os, const _Node& node)
{
if (node._node)
node._node->Print(os);
else
os << "mesh node at hexahedron corner is null" << '\n';
return os;
}
};
// --------------------------------------------------------------------------------
struct _Link // link connecting two _Node's
@ -783,6 +793,21 @@ namespace
vector< _Node* > _fIntNodes; // _Node's at _fIntPoints
vector< _Link > _splits;
_Link(): _faces{ 0, 0 } {}
friend std::ostream& operator<<(std::ostream& os, const _Link& link)
{
os << "Link:\n";
for (int i = 0; i < 2; ++i)
{
if (link._nodes[i])
os << *link._nodes[i];
else
os << "link node with index " << i << " is null" << '\n';
}
return os;
}
};
// --------------------------------------------------------------------------------
struct _OrientedLink
@ -852,6 +877,16 @@ namespace
}
}
}
friend std::ostream& operator<<(std::ostream& os, const _OrientedLink& link)
{
if (link._link)
os << "Oriented " << *link._link;
else
os << "Oriented link is null" << '\n';
return os;
}
};
// --------------------------------------------------------------------------------
struct _SplitIterator //! set to _hexLinks splits on one side of INTERNAL FACEs
@ -1114,6 +1149,26 @@ namespace
TID2Nb::iterator id2nb = id2nbMap.insert( s0 ).first;
id2nb->second++;
}
// Called by compute()
//================================================================================
_Face* createPolygon(const SMESH_Block::TShapeID& name);
void closePolygon(
_Face* polygon, _Node* n2, _Node* nFirst, _Face& quad, vector<_Node*>& chainNodes, size_t& nbUsedEdgeNodes, _Face* prevPolyg);
void connectPolygonLinks(
const Solid* solid, _Face* polygon, _Face& quad, vector<_Node*>& chainNodes, std::vector<_OrientedLink>& splits, const bool toCheckSideDivision);
bool collectSplits(std::vector<_OrientedLink>& splits, const _Face& quad, _Face* polygon, int quadIndex);
std::set<TGeomID> getConcaveFaces(const Solid* solid);
std::vector<_Node*> getChainNodes(const Solid* solid, const IsInternalFlag intFlag);
void clearHexUsedInFace();
void clearIntUsedInFace();
void addPolygonsToLinks();
std::vector<_OrientedLink*> getFreeLinks();
int notUsedIntersectionNodesToVInt();
bool createPolygons(const bool hasEdgeIntersections, const IsInternalFlag intFlag);
void setNamesForNoNamePolygons();
bool createVolume(const Solid* solid);
//================================================================================
}; // class Hexahedron
#ifdef WITH_TBB
@ -2974,23 +3029,20 @@ namespace
//================================================================================
/*!
* \brief Compute mesh volumes resulted from intersection of the Hexahedron
* \brief Collected faces can be used to avoid connecting nodes laying on them
*/
bool Hexahedron::compute( const Solid* solid, const IsInternalFlag intFlag )
std::set<TGeomID> Hexahedron::getConcaveFaces(const Solid* solid)
{
_polygons.clear();
_polygons.reserve( 20 );
MESSAGE("Collect concave faces...");
for ( int iN = 0; iN < 8; ++iN )
_hexNodes[iN]._usedInFace = 0;
if ( intFlag & IS_CUT_BY_INTERNAL_FACE && !_grid->_toAddEdges ) // Issue #19913
preventVolumesOverlapping();
std::set< TGeomID > concaveFaces; // to avoid connecting nodes laying on them
if ( solid->HasConcaveVertex() )
if (!solid->HasConcaveVertex())
{
MESSAGE("There's no concave faces here. Return.");
return {};
}
std::set< TGeomID > concaveFaces;
for ( const E_IntersectPoint* ip : _eIntPoints )
{
if ( const ConcaveFace* cf = solid->GetConcave( ip->_shapeID ))
@ -3005,26 +3057,34 @@ namespace
if ( this->hasEdgesAround( cf ))
concaveFaces.insert( cf->_concaveFace );
}
MESSAGE("Collected concave faces: " << concaveFaces.size());
return concaveFaces;
}
// Create polygons from quadrangles
// --------------------------------
vector< _OrientedLink > splits;
vector<_Node*> chainNodes;
_Face* coplanarPolyg;
const bool hasEdgeIntersections = !_eIntPoints.empty();
const bool toCheckSideDivision = isImplementEdges() || intFlag & IS_CUT_BY_INTERNAL_FACE;
for ( int iF = 0; iF < 6; ++iF ) // loop on 6 sides of a hexahedron
//================================================================================
/*!
* \brief Creates a new polygone with a given name
*/
Hexahedron::_Face* Hexahedron::createPolygon(const SMESH_Block::TShapeID& name)
{
_Face& quad = _hexQuads[ iF ] ;
MESSAGE("Create a polygon with a name: " << name);
_polygons.resize( _polygons.size() + 1 );
_Face* polygon = &_polygons.back();
polygon->_polyLinks.reserve( 20 );
polygon->_name = quad._name;
polygon->_name = name;
return polygon;
}
//================================================================================
/*!
* \brief Collects split links on 4 sides of a quadrangle.
* Returns false if the quad on FACE is not split.
*/
bool Hexahedron::collectSplits(std::vector<_OrientedLink>& splits, const _Face& quad, _Face* polygon, int quadIndex)
{
MESSAGE("Collect splits...");
splits.clear();
for ( int iE = 0; iE < 4; ++iE ) // loop on 4 sides of a quadrangle
@ -3032,15 +3092,52 @@ namespace
splits.push_back( quad._links[ iE ].ResultLink( iS ));
if ( splits.size() == 4 &&
isQuadOnFace( iF )) // check if a quad on FACE is not split
isQuadOnFace( quadIndex )) // check if a quad on FACE is not split
{
MESSAGE("The quad on FACE is not split. Swap splits with polygon links.");
polygon->_links.swap( splits );
continue; // goto the next quad
return false;
}
MESSAGE("Collected splits: " << splits.size());
return true;
}
void Hexahedron::closePolygon(
_Face* polygon, _Node* n2, _Node* nFirst, _Face& quad, vector<_Node*>& chainNodes, size_t& nbUsedEdgeNodes, _Face* prevPolyg)
{
MESSAGE("Try to close polygon. Number of used edge nodes: " << nbUsedEdgeNodes);
if ( nFirst == n2 )
return;
if ( !findChain( n2, nFirst, quad, chainNodes ))
{
if ( !closePolygon( polygon, chainNodes ))
if ( !isImplementEdges() )
chainNodes.push_back( nFirst );
}
for ( size_t i = 1; i < chainNodes.size(); ++i )
{
polygon->AddPolyLink( chainNodes[i-1], chainNodes[i], prevPolyg );
nbUsedEdgeNodes += bool( chainNodes[i]->IsUsedInFace( polygon ));
MESSAGE("Added link for nodes:\n" << *chainNodes[i-1] << *chainNodes[i]);
}
MESSAGE("Polygon was closed. Number of used edge nodes: " << nbUsedEdgeNodes);
}
void Hexahedron::connectPolygonLinks(
const Solid* solid, _Face* polygon, _Face& quad, vector<_Node*>& chainNodes, std::vector<_OrientedLink>& splits, const bool toCheckSideDivision)
{
MESSAGE("Connect polygon links...");
const std::set<TGeomID> concaveFaces = getConcaveFaces(solid); // to avoid connecting nodes laying on them
// add splits of links to a polygon and add _polyLinks to make
// polygon's boundary closed
int nbSplits = splits.size();
if (( nbSplits == 1 ) &&
( quad._eIntNodes.empty() ||
@ -3055,6 +3152,7 @@ namespace
size_t nbUsedEdgeNodes = 0;
_Face* prevPolyg = 0; // polygon previously created from this quad
MESSAGE("Number of splits: " << nbSplits);
while ( nbSplits > 0 )
{
size_t iS = 0;
@ -3063,10 +3161,7 @@ namespace
if ( !polygon->_links.empty() )
{
_polygons.resize( _polygons.size() + 1 );
polygon = &_polygons.back();
polygon->_polyLinks.reserve( 20 );
polygon->_name = quad._name;
polygon = createPolygon(quad._name);
}
polygon->_links.push_back( splits[ iS ] );
splits[ iS++ ]._link = 0;
@ -3077,15 +3172,18 @@ namespace
for ( ; nFirst != n2 && iS < splits.size(); ++iS )
{
_OrientedLink& split = splits[ iS ];
MESSAGE("Current split: " << split);
if ( !split ) continue;
n1 = split.FirstNode();
MESSAGE("\n\tnFirst: " << *nFirst << "\tn1: " << *n1 << "\tn2: " << *n2);
if ( n1 == n2 &&
n1->_intPoint &&
(( n1->_intPoint->_faceIDs.size() > 1 && toCheckSideDivision ) ||
( n1->_isInternalFlags )))
{
// n1 is at intersection with EDGE
MESSAGE("n1 is at intersection with EDGE");
if ( findChainOnEdge( splits, polygon->_links.back(), split, concaveFaces,
iS, quad, chainNodes ))
{
@ -3102,6 +3200,7 @@ namespace
else if ( n1 != n2 )
{
// try to connect to intersections with EDGEs
MESSAGE("n1 != n2 -> try to connect to intersections with EDGEs");
if ( quad._eIntNodes.size() > nbUsedEdgeNodes &&
findChain( n2, n1, quad, chainNodes ))
{
@ -3120,6 +3219,7 @@ namespace
// try to connect to a split ending on the same FACE
else
{
MESSAGE("try to connect to a split ending on the same FACE");
_OrientedLink foundSplit;
for ( size_t i = iS; i < splits.size() && !foundSplit; ++i )
if (( foundSplit = splits[ i ]) &&
@ -3131,6 +3231,8 @@ namespace
{
foundSplit._link = 0;
}
MESSAGE("Found split: " << foundSplit);
if ( foundSplit )
{
if ( n2 != foundSplit.FirstNode() )
@ -3147,7 +3249,11 @@ namespace
polygon->AddPolyLink( n2, n1, prevPolyg );
}
}
} // if ( n1 != n2 )
}
else
{
MESSAGE("n1 == n2, split will be added as is");
}// if ( n1 != n2 )
polygon->_links.push_back( split );
split._link = 0;
@ -3156,20 +3262,7 @@ namespace
} // loop on splits
if ( nFirst != n2 ) // close a polygon
{
if ( !findChain( n2, nFirst, quad, chainNodes ))
{
if ( !closePolygon( polygon, chainNodes ))
if ( !isImplementEdges() )
chainNodes.push_back( nFirst );
}
for ( size_t i = 1; i < chainNodes.size(); ++i )
{
polygon->AddPolyLink( chainNodes[i-1], chainNodes[i], prevPolyg );
nbUsedEdgeNodes += bool( chainNodes[i]->IsUsedInFace( polygon ));
}
}
closePolygon(polygon, n2, nFirst, quad, chainNodes, nbUsedEdgeNodes, prevPolyg);
if ( polygon->_links.size() < 3 && nbSplits > 0 )
{
@ -3178,31 +3271,81 @@ namespace
}
} // while ( nbSplits > 0 )
MESSAGE("Polygon's links were connected");
}
//================================================================================
/*!
* \brief Collects chain nodes
*/
std::vector<Hexahedron::_Node*> Hexahedron::getChainNodes(const Solid* solid, const IsInternalFlag intFlag)
{
std::vector< _OrientedLink > splits;
std::vector<_Node*> chainNodes;
const bool toCheckSideDivision = isImplementEdges() || intFlag & IS_CUT_BY_INTERNAL_FACE;
for ( int iF = 0; iF < 6; ++iF ) // loop on 6 sides of a hexahedron
{
_Face& quad = _hexQuads[ iF ] ;
_Face* polygon = createPolygon(quad._name);
if (!collectSplits(splits, quad, polygon, iF))
continue; // goto the next quad
connectPolygonLinks(solid, polygon, quad, chainNodes, splits, toCheckSideDivision);
if ( polygon->_links.size() < 3 )
{
_polygons.pop_back();
}
} // loop on 6 hexahedron sides
// Create polygons closing holes in a polyhedron
// ----------------------------------------------
return chainNodes;
}
// clear _usedInFace
//================================================================================
/*!
* \brief Reset used in faces for hex nodes
*/
void Hexahedron::clearHexUsedInFace()
{
for ( int iN = 0; iN < 8; ++iN )
_hexNodes[iN]._usedInFace = 0;
}
//================================================================================
/*!
* \brief Reset used in faces for int nodes
*/
void Hexahedron::clearIntUsedInFace()
{
for ( size_t iN = 0; iN < _intNodes.size(); ++iN )
_intNodes[ iN ]._usedInFace = 0;
}
// add polygons to their links and mark used nodes
for ( size_t iP = 0; iP < _polygons.size(); ++iP )
//================================================================================
/*!
* \brief Add polygons to their links and mark used nodes
*/
void Hexahedron::addPolygonsToLinks()
{
_Face& polygon = _polygons[ iP ];
for ( size_t iL = 0; iL < polygon._links.size(); ++iL )
for (_Face& polygon : _polygons)
{
polygon._links[ iL ].AddFace( &polygon );
polygon._links[ iL ].FirstNode()->_usedInFace = &polygon;
for (_OrientedLink& link : polygon._links)
{
link.AddFace( &polygon );
link.FirstNode()->_usedInFace = &polygon;
}
}
// find free links
vector< _OrientedLink* > freeLinks;
}
//================================================================================
/*!
* \brief Finds free links in polygons
*/
std::vector<Hexahedron::_OrientedLink*> Hexahedron::getFreeLinks()
{
std::vector< _OrientedLink* > freeLinks;
freeLinks.reserve(20);
for ( size_t iP = 0; iP < _polygons.size(); ++iP )
{
@ -3211,29 +3354,55 @@ namespace
if ( polygon._links[ iL ].NbFaces() < 2 )
freeLinks.push_back( & polygon._links[ iL ]);
}
return freeLinks;
}
//================================================================================
/*!
* \brief Put not used intersection nodes to _vIntNodes
*/
int Hexahedron::notUsedIntersectionNodesToVInt()
{
int nbVertexNodes = 0;
// TEMP: the loops below can be commented out for simplified testing
// for ( size_t iN = 0; iN < _vIntNodes.size(); ++iN )
// nbVertexNodes += ( !_vIntNodes[ iN ]->IsUsedInFace() );
// const double tol = 1e-3 * Min( Min( _sideLength[0], _sideLength[1] ), _sideLength[0] );
// for ( size_t iN = _nbFaceIntNodes; iN < _intNodes.size(); ++iN )
// {
// if ( _intNodes[ iN ].IsUsedInFace() ) continue;
// if ( dynamic_cast< const F_IntersectPoint* >( _intNodes[ iN ]._intPoint )) continue;
// _Node* equalNode =
// findEqualNode( _vIntNodes, _intNodes[ iN ].EdgeIntPnt(), tol*tol );
// if ( !equalNode )
// {
// _vIntNodes.push_back( &_intNodes[ iN ]);
// ++nbVertexNodes;
// }
// }
return nbVertexNodes;
}
//================================================================================
/*!
* \brief Creates polygons from free links
*/
bool Hexahedron::createPolygons(const bool hasEdgeIntersections, const IsInternalFlag intFlag)
{
MESSAGE("Create polygons from free links...");
// find free links
vector< _OrientedLink* > freeLinks = getFreeLinks();
int nbFreeLinks = freeLinks.size();
SCRUTE(nbFreeLinks);
if ( nbFreeLinks == 1 ) return false;
// put not used intersection nodes to _vIntNodes
int nbVertexNodes = 0; // nb not used vertex nodes
{
for ( size_t iN = 0; iN < _vIntNodes.size(); ++iN )
nbVertexNodes += ( !_vIntNodes[ iN ]->IsUsedInFace() );
const double tol = 1e-3 * Min( Min( _sideLength[0], _sideLength[1] ), _sideLength[0] );
for ( size_t iN = _nbFaceIntNodes; iN < _intNodes.size(); ++iN )
{
if ( _intNodes[ iN ].IsUsedInFace() ) continue;
if ( dynamic_cast< const F_IntersectPoint* >( _intNodes[ iN ]._intPoint )) continue;
_Node* equalNode =
findEqualNode( _vIntNodes, _intNodes[ iN ].EdgeIntPnt(), tol*tol );
if ( !equalNode )
{
_vIntNodes.push_back( &_intNodes[ iN ]);
++nbVertexNodes;
}
}
}
int nbVertexNodes = notUsedIntersectionNodesToVInt();
std::set<TGeomID> usedFaceIDs;
std::vector< TGeomID > faces;
@ -3246,6 +3415,7 @@ namespace
size_t iPolygon = _polygons.size();
while ( nbFreeLinks > 0 )
{
SCRUTE(nbFreeLinks);
if ( iPolygon == _polygons.size() )
{
_polygons.resize( _polygons.size() + 1 );
@ -3259,6 +3429,7 @@ namespace
if (( !hasEdgeIntersections ) ||
( nbFreeLinks < 4 && nbVertexNodes == 0 ))
{
MESSAGE("!hasEdgeIntersections || (nbFreeLinks < 4 && nbVertexNodes == 0), get a remaining link to start from...");
// get a remaining link to start from
for ( size_t iL = 0; iL < freeLinks.size() && !curLink; ++iL )
if (( curLink = freeLinks[ iL ] ))
@ -3282,150 +3453,153 @@ namespace
}
else // there are intersections with EDGEs
{
MESSAGE("there are intersections with EDGEs...");
// get a remaining link to start from, one lying on minimal nb of FACEs
{
typedef pair< TGeomID, int > TFaceOfLink;
TFaceOfLink faceOfLink( -1, -1 );
TFaceOfLink facesOfLink[3] = { faceOfLink, faceOfLink, faceOfLink };
for ( size_t iL = 0; iL < freeLinks.size(); ++iL )
if ( freeLinks[ iL ] )
{
faces = freeLinks[ iL ]->GetNotUsedFace( usedFaceIDs );
if ( faces.size() == 1 )
{
faceOfLink = TFaceOfLink( faces[0], iL );
if ( !freeLinks[ iL ]->HasEdgeNodes() )
break;
facesOfLink[0] = faceOfLink;
}
else if ( facesOfLink[0].first < 0 )
{
faceOfLink = TFaceOfLink(( faces.empty() ? -1 : faces[0]), iL );
facesOfLink[ 1 + faces.empty() ] = faceOfLink;
}
}
for ( int i = 0; faceOfLink.first < 0 && i < 3; ++i )
faceOfLink = facesOfLink[i];
// {
// typedef pair< TGeomID, int > TFaceOfLink;
// TFaceOfLink faceOfLink( -1, -1 );
// TFaceOfLink facesOfLink[3] = { faceOfLink, faceOfLink, faceOfLink };
// for ( size_t iL = 0; iL < freeLinks.size(); ++iL )
// if ( freeLinks[ iL ] )
// {
// faces = freeLinks[ iL ]->GetNotUsedFace( usedFaceIDs );
// if ( faces.size() == 1 )
// {
// faceOfLink = TFaceOfLink( faces[0], iL );
// if ( !freeLinks[ iL ]->HasEdgeNodes() )
// break;
// facesOfLink[0] = faceOfLink;
// }
// else if ( facesOfLink[0].first < 0 )
// {
// faceOfLink = TFaceOfLink(( faces.empty() ? -1 : faces[0]), iL );
// facesOfLink[ 1 + faces.empty() ] = faceOfLink;
// }
// }
// for ( int i = 0; faceOfLink.first < 0 && i < 3; ++i )
// faceOfLink = facesOfLink[i];
if ( faceOfLink.first < 0 ) // all faces used
{
for ( size_t iL = 0; iL < freeLinks.size() && faceOfLink.first < 1; ++iL )
if (( curLink = freeLinks[ iL ]))
{
faceOfLink.first =
curLink->FirstNode()->IsLinked( curLink->LastNode()->_intPoint );
faceOfLink.second = iL;
}
usedFaceIDs.clear();
}
curFace = faceOfLink.first;
curLink = freeLinks[ faceOfLink.second ];
freeLinks[ faceOfLink.second ] = 0;
}
usedFaceIDs.insert( curFace );
polygon._links.push_back( *curLink );
--nbFreeLinks;
// if ( faceOfLink.first < 0 ) // all faces used
// {
// for ( size_t iL = 0; iL < freeLinks.size() && faceOfLink.first < 1; ++iL )
// if (( curLink = freeLinks[ iL ]))
// {
// faceOfLink.first =
// curLink->FirstNode()->IsLinked( curLink->LastNode()->_intPoint );
// faceOfLink.second = iL;
// }
// usedFaceIDs.clear();
// }
// curFace = faceOfLink.first;
// curLink = freeLinks[ faceOfLink.second ];
// freeLinks[ faceOfLink.second ] = 0;
// }
// usedFaceIDs.insert( curFace );
// polygon._links.push_back( *curLink );
// --nbFreeLinks;
// find all links lying on a curFace
do
{
// go forward from curLink
curNode = curLink->LastNode();
curLink = 0;
for ( size_t iL = 0; iL < freeLinks.size() && !curLink; ++iL )
if ( freeLinks[ iL ] &&
freeLinks[ iL ]->FirstNode() == curNode &&
freeLinks[ iL ]->LastNode()->IsOnFace( curFace ))
{
curLink = freeLinks[ iL ];
freeLinks[ iL ] = 0;
polygon._links.push_back( *curLink );
--nbFreeLinks;
}
} while ( curLink );
// do
// {
// // go forward from curLink
// curNode = curLink->LastNode();
// curLink = 0;
// for ( size_t iL = 0; iL < freeLinks.size() && !curLink; ++iL )
// if ( freeLinks[ iL ] &&
// freeLinks[ iL ]->FirstNode() == curNode &&
// freeLinks[ iL ]->LastNode()->IsOnFace( curFace ))
// {
// curLink = freeLinks[ iL ];
// freeLinks[ iL ] = 0;
// polygon._links.push_back( *curLink );
// --nbFreeLinks;
// }
// } while ( curLink );
std::reverse( polygon._links.begin(), polygon._links.end() );
// std::reverse( polygon._links.begin(), polygon._links.end() );
curLink = & polygon._links.back();
do
{
// go backward from curLink
curNode = curLink->FirstNode();
curLink = 0;
for ( size_t iL = 0; iL < freeLinks.size() && !curLink; ++iL )
if ( freeLinks[ iL ] &&
freeLinks[ iL ]->LastNode() == curNode &&
freeLinks[ iL ]->FirstNode()->IsOnFace( curFace ))
{
curLink = freeLinks[ iL ];
freeLinks[ iL ] = 0;
polygon._links.push_back( *curLink );
--nbFreeLinks;
}
} while ( curLink );
// curLink = & polygon._links.back();
// do
// {
// // go backward from curLink
// curNode = curLink->FirstNode();
// curLink = 0;
// for ( size_t iL = 0; iL < freeLinks.size() && !curLink; ++iL )
// if ( freeLinks[ iL ] &&
// freeLinks[ iL ]->LastNode() == curNode &&
// freeLinks[ iL ]->FirstNode()->IsOnFace( curFace ))
// {
// curLink = freeLinks[ iL ];
// freeLinks[ iL ] = 0;
// polygon._links.push_back( *curLink );
// --nbFreeLinks;
// }
// } while ( curLink );
curNode = polygon._links.back().FirstNode();
// curNode = polygon._links.back().FirstNode();
if ( polygon._links[0].LastNode() != curNode )
{
if ( nbVertexNodes > 0 )
{
// add links with _vIntNodes if not already used
chainNodes.clear();
for ( size_t iN = 0; iN < _vIntNodes.size(); ++iN )
if ( !_vIntNodes[ iN ]->IsUsedInFace() &&
_vIntNodes[ iN ]->IsOnFace( curFace ))
{
_vIntNodes[ iN ]->_usedInFace = &polygon;
chainNodes.push_back( _vIntNodes[ iN ] );
}
if ( chainNodes.size() > 1 &&
curFace != _grid->PseudoIntExtFaceID() ) /////// TODO
{
sortVertexNodes( chainNodes, curNode, curFace );
}
for ( size_t i = 0; i < chainNodes.size(); ++i )
{
polygon.AddPolyLink( chainNodes[ i ], curNode );
curNode = chainNodes[ i ];
freeLinks.push_back( &polygon._links.back() );
++nbFreeLinks;
}
nbVertexNodes -= chainNodes.size();
}
// if ( polygon._links.size() > 1 )
{
polygon.AddPolyLink( polygon._links[0].LastNode(), curNode );
freeLinks.push_back( &polygon._links.back() );
++nbFreeLinks;
}
}
// if ( polygon._links[0].LastNode() != curNode )
// {
// if ( nbVertexNodes > 0 )
// {
// // add links with _vIntNodes if not already used
// chainNodes.clear();
// for ( size_t iN = 0; iN < _vIntNodes.size(); ++iN )
// if ( !_vIntNodes[ iN ]->IsUsedInFace() &&
// _vIntNodes[ iN ]->IsOnFace( curFace ))
// {
// _vIntNodes[ iN ]->_usedInFace = &polygon;
// chainNodes.push_back( _vIntNodes[ iN ] );
// }
// if ( chainNodes.size() > 1 &&
// curFace != _grid->PseudoIntExtFaceID() ) /////// TODO
// {
// sortVertexNodes( chainNodes, curNode, curFace );
// }
// for ( size_t i = 0; i < chainNodes.size(); ++i )
// {
// polygon.AddPolyLink( chainNodes[ i ], curNode );
// curNode = chainNodes[ i ];
// freeLinks.push_back( &polygon._links.back() );
// ++nbFreeLinks;
// }
// nbVertexNodes -= chainNodes.size();
// }
// // if ( polygon._links.size() > 1 )
// {
// polygon.AddPolyLink( polygon._links[0].LastNode(), curNode );
// freeLinks.push_back( &polygon._links.back() );
// ++nbFreeLinks;
// }
// }
} // if there are intersections with EDGEs
if ( polygon._links.size() < 2 ||
polygon._links[0].LastNode() != polygon._links.back().FirstNode() )
{
_polygons.clear();
break; // closed polygon not found -> invalid polyhedron
}
// if ( polygon._links.size() < 2 ||
// polygon._links[0].LastNode() != polygon._links.back().FirstNode() )
// {
// _polygons.clear();
// break; // closed polygon not found -> invalid polyhedron
// }
if ( polygon._links.size() == 2 )
{
if ( freeLinks.back() == &polygon._links.back() )
{
freeLinks.pop_back();
--nbFreeLinks;
}
if ( polygon._links.front().NbFaces() > 0 )
polygon._links.back().AddFace( polygon._links.front()._link->_faces[0] );
if ( polygon._links.back().NbFaces() > 0 )
polygon._links.front().AddFace( polygon._links.back()._link->_faces[0] );
MESSAGE("polygon._links.size() == 2...");
// if ( freeLinks.back() == &polygon._links.back() )
// {
// freeLinks.pop_back();
// --nbFreeLinks;
// }
// if ( polygon._links.front().NbFaces() > 0 )
// polygon._links.back().AddFace( polygon._links.front()._link->_faces[0] );
// if ( polygon._links.back().NbFaces() > 0 )
// polygon._links.front().AddFace( polygon._links.back()._link->_faces[0] );
if ( iPolygon == _polygons.size()-1 )
_polygons.pop_back();
// if ( iPolygon == _polygons.size()-1 )
// _polygons.pop_back();
}
else // polygon._links.size() >= 2
{
MESSAGE("polygon._links.size() >= 2, add polygon to its links...");
// add polygon to its links
for ( size_t iL = 0; iL < polygon._links.size(); ++iL )
{
@ -3435,26 +3609,27 @@ namespace
if ( /*hasEdgeIntersections &&*/ iPolygon == _polygons.size() - 1 )
{
// check that a polygon does not lie on a hexa side
coplanarPolyg = 0;
for ( size_t iL = 0; iL < polygon._links.size() && !coplanarPolyg; ++iL )
{
if ( polygon._links[ iL ].NbFaces() < 2 )
continue; // it's a just added free link
// look for a polygon made on a hexa side and sharing
// two or more haxa links
size_t iL2;
coplanarPolyg = polygon._links[ iL ]._link->_faces[0];
for ( iL2 = iL + 1; iL2 < polygon._links.size(); ++iL2 )
if ( polygon._links[ iL2 ]._link->_faces[0] == coplanarPolyg &&
!coplanarPolyg->IsPolyLink( polygon._links[ iL ]) &&
!coplanarPolyg->IsPolyLink( polygon._links[ iL2 ]) &&
coplanarPolyg < & _polygons[ nbQuadPolygons ])
break;
if ( iL2 == polygon._links.size() )
coplanarPolyg = 0;
}
_Face* coplanarPolyg = nullptr;
// for ( size_t iL = 0; iL < polygon._links.size() && !coplanarPolyg; ++iL )
// {
// if ( polygon._links[ iL ].NbFaces() < 2 )
// continue; // it's a just added free link
// // look for a polygon made on a hexa side and sharing
// // two or more haxa links
// size_t iL2;
// coplanarPolyg = polygon._links[ iL ]._link->_faces[0];
// for ( iL2 = iL + 1; iL2 < polygon._links.size(); ++iL2 )
// if ( polygon._links[ iL2 ]._link->_faces[0] == coplanarPolyg &&
// !coplanarPolyg->IsPolyLink( polygon._links[ iL ]) &&
// !coplanarPolyg->IsPolyLink( polygon._links[ iL2 ]) &&
// coplanarPolyg < & _polygons[ nbQuadPolygons ])
// break;
// if ( iL2 == polygon._links.size() )
// coplanarPolyg = 0;
// }
if ( coplanarPolyg ) // coplanar polygon found
{
MESSAGE("coplanar polygon found");
freeLinks.resize( freeLinks.size() - polygon._polyLinks.size() );
nbFreeLinks -= polygon._polyLinks.size();
@ -3531,16 +3706,16 @@ namespace
} // end of case ( polygon._links.size() > 2 )
} // while ( nbFreeLinks > 0 )
for ( auto & face_ip : tmpAddedFace )
{
curFace = face_ip.first;
for ( const B_IntersectPoint* ip : face_ip.second )
{
auto it = std::find( ip->_faceIDs.begin(), ip->_faceIDs.end(), curFace );
if ( it != ip->_faceIDs.end() )
ip->_faceIDs.erase( it );
}
}
// for ( auto & face_ip : tmpAddedFace )
// {
// curFace = face_ip.first;
// for ( const B_IntersectPoint* ip : face_ip.second )
// {
// auto it = std::find( ip->_faceIDs.begin(), ip->_faceIDs.end(), curFace );
// if ( it != ip->_faceIDs.end() )
// ip->_faceIDs.erase( it );
// }
// }
if ( _polygons.size() < 3 )
return false;
@ -3548,18 +3723,29 @@ namespace
// check volume size
double volSize = 0;
_hasTooSmall = ! checkPolyhedronSize( intFlag & IS_CUT_BY_INTERNAL_FACE, volSize );
_volumeDefs._size = volSize;
for ( size_t i = 0; i < 8; ++i )
if ( _hexNodes[ i ]._intPoint == &ipTmp )
_hexNodes[ i ]._intPoint = 0;
if ( _hasTooSmall )
return false; // too small volume
return !_hasTooSmall; // too small volume
}
/*!
* \brief Sets names for no-name polygons
*/
void Hexahedron::setNamesForNoNamePolygons()
{
MESSAGE("Try to find out names of no-name polygons...");
// Try to find out names of no-name polygons (issue # 19887)
if ( _grid->IsToRemoveExcessEntities() && _polygons.back()._name == SMESH_Block::ID_NONE )
{
// TODO: why do we check only last face in vector?
if (!_grid->IsToRemoveExcessEntities() || _polygons.back()._name != SMESH_Block::ID_NONE )
return;
MESSAGE("No-name was found! It will be named after grid coordinates.");
gp_XYZ uvwCenter =
0.5 * ( _grid->_coords[0][_i] + _grid->_coords[0][_i+1] ) * _grid->_axes[0] +
0.5 * ( _grid->_coords[1][_j] + _grid->_coords[1][_j+1] ) * _grid->_axes[1] +
@ -3586,6 +3772,13 @@ namespace
}
}
/*!
* \brief Creates a volume as a classic element (Hexa, Tetra, Penta, Pyra) or non-specific polyhedron
*/
bool Hexahedron::createVolume(const Solid* solid)
{
MESSAGE("Create a volume...");
_volumeDefs._nodes.clear();
_volumeDefs._quantities.clear();
_volumeDefs._names.clear();
@ -3609,20 +3802,63 @@ namespace
else if ( nbNodes == 5 && nbPolygons == 5 ) isClassicElem = addPyra ();
if ( !isClassicElem )
{
MESSAGE("It's not a classic element. Create a volume by polygons and links...");
for ( size_t iF = 0; iF < _polygons.size(); ++iF )
{
const size_t nbLinks = _polygons[ iF ]._links.size();
SCRUTE(nbLinks);
if ( nbLinks < 3 ) continue;
_volumeDefs._quantities.push_back( nbLinks );
_volumeDefs._names.push_back( _polygons[ iF ]._name );
SCRUTE(_polygons[ iF ]._name);
for ( size_t iL = 0; iL < nbLinks; ++iL )
_volumeDefs._nodes.push_back( _polygons[ iF ]._links[ iL ].FirstNode() );
}
}
_volumeDefs._solidID = solid->ID();
_volumeDefs._size = volSize;
return !_volumeDefs._nodes.empty();
}
/*!
* \brief Compute mesh volumes resulted from intersection of the Hexahedron
*/
bool Hexahedron::compute( const Solid* solid, const IsInternalFlag intFlag )
{
MESSAGE("Compute solid " << solid->ID() << " with internal flag " << intFlag << " ...");
// Reset polygons
_polygons.clear();
_polygons.reserve( 20 );
clearHexUsedInFace();
// Check volumes
if ( intFlag & IS_CUT_BY_INTERNAL_FACE && !_grid->_toAddEdges ) // Issue #19913
preventVolumesOverlapping();
// Create polygons from quadrangles
// --------------------------------
vector<_Node*> chainNodes = getChainNodes(solid, intFlag);
// TODO: check if we can safely move it inside createPolygons()
const bool hasEdgeIntersections = !_eIntPoints.empty();
// Create polygons closing holes in a polyhedron
// ----------------------------------------------
clearIntUsedInFace();
addPolygonsToLinks();
// Create polygons
if (!createPolygons(hasEdgeIntersections, intFlag))
return false;
// Fix polygons' names
setNamesForNoNamePolygons();
return createVolume(solid);
}
template<typename Type>
@ -3651,11 +3887,12 @@ namespace
{
// A simple case - just one task executed in one thread.
// TODO: check if it's faster to do it sequentially
threads.reserve(numTasksTotal);
for (; it < last; ++it)
{
threads.emplace_back(f, std::ref(*it));
}
// threads.reserve(numTasksTotal);
// for (; it < last; ++it)
// {
// threads.emplace_back(f, std::ref(*it));
// }
std::for_each(it, last, f);
}
else
{
@ -3818,16 +4055,16 @@ namespace
#endif
// simplify polyhedrons
if ( _grid->IsToRemoveExcessEntities() )
{
for ( size_t i = 0; i < intHexa.size(); ++i )
if ( Hexahedron * hex = intHexa[ i ] )
hex->removeExcessSideDivision( allHexa );
// if ( _grid->IsToRemoveExcessEntities() )
// {
// for ( size_t i = 0; i < intHexa.size(); ++i )
// if ( Hexahedron * hex = intHexa[ i ] )
// hex->removeExcessSideDivision( allHexa );
for ( size_t i = 0; i < intHexa.size(); ++i )
if ( Hexahedron * hex = intHexa[ i ] )
hex->removeExcessNodes( allHexa );
}
// for ( size_t i = 0; i < intHexa.size(); ++i )
// if ( Hexahedron * hex = intHexa[ i ] )
// hex->removeExcessNodes( allHexa );
// }
// add volumes
for ( size_t i = 0; i < intHexa.size(); ++i )