#include "utilities.h"
-#include <BRep_Tool.hxx>
+#include <BRepAdaptor_Surface.hxx>
#include <BRepClass3d_SolidClassifier.hxx>
+#include <BRep_Tool.hxx>
#include <ElCLib.hxx>
#include <Extrema_GenExtPS.hxx>
#include <Extrema_POnSurf.hxx>
return;
}
- issimple[iNode] = (listNewNodes.size()==nbSteps);
+ issimple[iNode] = (listNewNodes.size()==nbSteps); // is node medium
itNN[ iNode ] = listNewNodes.begin();
prevNod[ iNode ] = node;
nextNod[ iNode ] = listNewNodes.front();
- if( !issimple[iNode] ) {
+ if( !elem->IsQuadratic() || !issimple[iNode] ) {
if ( prevNod[ iNode ] != nextNod [ iNode ])
iNotSameNode = iNode;
else {
//cout<<" nbSame = "<<nbSame<<endl;
if ( nbSame == nbNodes || nbSame > 2) {
- //MESSAGE( " Too many same nodes of element " << elem->GetID() );
- INFOS( " Too many same nodes of element " << elem->GetID() );
+ MESSAGE( " Too many same nodes of element " << elem->GetID() );
+ //INFOS( " Too many same nodes of element " << elem->GetID() );
return;
}
while ( nodeIt->more() )
xyz.push_back( TNodeXYZ( cast2Node( nodeIt->next() )));
+ int i, nbNodes = element->NbNodes();
+
if ( element->GetType() == SMDSAbs_Face ) // --------------------------------------------------
{
- // gravity center
- gp_XYZ gc(0,0,0);
- gc = accumulate( xyz.begin(), xyz.end(), gc );
- gc /= element->NbNodes();
-
- // compute face normal using gc
- gp_Vec normal(0,0,0);
+ // compute face normal
+ gp_Vec faceNorm(0,0,0);
xyz.push_back( xyz.front() );
- for ( int i = 0; i < element->NbNodes(); ++i )
+ for ( i = 0; i < nbNodes; ++i )
{
- gp_Vec edge( xyz[i], xyz[i+1]);
- gp_Vec n2gc( xyz[i], gc );
- normal += edge ^ n2gc;
+ gp_Vec edge1( xyz[i+1], xyz[i]);
+ gp_Vec edge2( xyz[i+1], xyz[(i+2)%nbNodes] );
+ faceNorm += edge1 ^ edge2;
}
- double faceDoubleArea = normal.Magnitude();
- if ( faceDoubleArea <= numeric_limits<double>::min() )
- return true; // invalid face
- normal /= faceDoubleArea;
+ double normSize = faceNorm.Magnitude();
+ if ( normSize <= tol )
+ {
+ // degenerated face: point is out if it is out of all face edges
+ for ( i = 0; i < nbNodes; ++i )
+ {
+ SMDS_MeshNode n1( xyz[i].X(), xyz[i].Y(), xyz[i].Z() );
+ SMDS_MeshNode n2( xyz[i+1].X(), xyz[i+1].Y(), xyz[i+1].Z() );
+ SMDS_MeshEdge edge( &n1, &n2 );
+ if ( !isOut( &edge, point, tol ))
+ return false;
+ }
+ return true;
+ }
+ faceNorm /= normSize;
// check if the point lays on face plane
gp_Vec n2p( xyz[0], point );
- if ( fabs( n2p * normal ) > tol )
+ if ( fabs( n2p * faceNorm ) > tol )
return true; // not on face plane
- // check if point is out of face boundary
- int i, out = false;
- for ( i = 0; !out && i < element->NbNodes(); ++i )
+ // check if point is out of face boundary:
+ // define it by closest transition of a ray point->infinity through face boundary
+ // on the face plane.
+ // First, find normal of a plane perpendicular to face plane, to be used as a cutting tool
+ // to find intersections of the ray with the boundary.
+ gp_Vec ray = n2p;
+ gp_Vec plnNorm = ray ^ faceNorm;
+ normSize = plnNorm.Magnitude();
+ if ( normSize <= tol ) return false; // point coincides with the first node
+ plnNorm /= normSize;
+ // for each node of the face, compute its signed distance to the plane
+ vector<double> dist( nbNodes + 1);
+ for ( i = 0; i < nbNodes; ++i )
{
- gp_Vec edge( xyz[i], xyz[i+1]);
- gp_Vec n2p ( xyz[i], point );
- gp_Vec cross = edge ^ n2p;
- out = ( cross * normal < -tol );
- }
- if ( out && element->IsPoly() )
+ gp_Vec n2p( xyz[i], point );
+ dist[i] = n2p * plnNorm;
+ }
+ dist.back() = dist.front();
+ // find the closest intersection
+ int iClosest = -1;
+ double rClosest, distClosest = 1e100;;
+ gp_Pnt pClosest;
+ for ( i = 0; i < nbNodes; ++i )
{
- // define point position by the closest edge
- double minDist = numeric_limits<double>::max();
- int iMinDist;
- for ( i = 0; i < element->NbNodes(); ++i )
+ double r;
+ if ( fabs( dist[i]) < tol )
+ r = 0.;
+ else if ( fabs( dist[i+1]) < tol )
+ r = 1.;
+ else if ( dist[i] * dist[i+1] < 0 )
+ r = dist[i] / ( dist[i] - dist[i+1] );
+ else
+ continue; // no intersection
+ gp_Pnt pInt = xyz[i] * (1.-r) + xyz[i+1] * r;
+ gp_Vec p2int ( point, pInt);
+ if ( p2int * ray > -tol ) // right half-space
{
- gp_Vec edge( xyz[i], xyz[i+1]);
- gp_Vec n1p ( xyz[i], point);
- double dist = ( edge ^ n1p ).Magnitude() / edge.Magnitude();
- if ( dist < minDist )
- iMinDist = i;
+ double intDist = p2int.SquareMagnitude();
+ if ( intDist < distClosest )
+ {
+ iClosest = i;
+ rClosest = r;
+ pClosest = pInt;
+ distClosest = intDist;
+ }
}
- gp_Vec edge( xyz[iMinDist], xyz[iMinDist+1]);
- gp_Vec n2p ( xyz[iMinDist], point );
- gp_Vec cross = edge ^ n2p;
- out = ( cross * normal < -tol );
}
- return out;
+ if ( iClosest < 0 )
+ return true; // no intesections - out
+
+ // analyse transition
+ gp_Vec edge( xyz[iClosest], xyz[iClosest+1] );
+ gp_Vec edgeNorm = -( edge ^ faceNorm ); // normal to intersected edge pointing out of face
+ gp_Vec p2int ( point, pClosest );
+ bool out = (edgeNorm * p2int) < -tol;
+ if ( rClosest > 0. && rClosest < 1. ) // not node intersection
+ return out;
+
+ // ray pass through a face node; analyze transition through an adjacent edge
+ gp_Pnt p1 = xyz[ (rClosest == 0.) ? ((iClosest+nbNodes-1) % nbNodes) : (iClosest+1) ];
+ gp_Pnt p2 = xyz[ (rClosest == 0.) ? iClosest : ((iClosest+2) % nbNodes) ];
+ gp_Vec edgeAdjacent( p1, p2 );
+ gp_Vec edgeNorm2 = -( edgeAdjacent ^ faceNorm );
+ bool out2 = (edgeNorm2 * p2int) < -tol;
+
+ bool covexCorner = ( edgeNorm * edgeAdjacent * (rClosest==1. ? 1. : -1.)) < 0;
+ return covexCorner ? (out || out2) : (out && out2);
}
if ( element->GetType() == SMDSAbs_Edge ) // --------------------------------------------------
{
- for ( int i = 1; i < element->NbNodes(); ++i )
+ // point is out of edge if it is NOT ON any straight part of edge
+ // (we consider quadratic edge as being composed of two straight parts)
+ for ( i = 1; i < nbNodes; ++i )
{
gp_Vec edge( xyz[i-1], xyz[i]);
gp_Vec n1p ( xyz[i-1], point);
double dist = ( edge ^ n1p ).Magnitude() / edge.Magnitude();
if ( dist > tol )
- return true;
+ continue;
gp_Vec n2p( xyz[i], point );
if ( fabs( edge.Magnitude() - n1p.Magnitude() - n2p.Magnitude()) > tol )
- return true;
+ continue;
+ return false; // point is ON this part
}
- return false;
+ return true;
}
// Node or 0D element -------------------------------------------------------------------------
{
case 4:
NewElem = theHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3], id, theForce3d);
break;
+ case 5:
+ NewElem = theHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3], aNds[4], id, theForce3d);
+ break;
case 6:
NewElem = theHelper.AddVolume(aNds[0], aNds[1], aNds[2], aNds[3], aNds[4], aNds[5], id, theForce3d);
break;
SMESH_subMesh* sm = smIt->next();
if ( SMESHDS_SubMesh *smDS = sm->GetSubMeshDS() ) {
aHelper.SetSubShape( sm->GetSubShape() );
- if ( !theForce3d) aHelper.SetCheckNodePosition(true);
nbCheckedElems += convertElemToQuadratic(smDS, aHelper, theForce3d);
}
}
NewVolume = aHelper.AddVolume(aNds[0], aNds[1], aNds[2],
aNds[3], id, theForce3d );
break;
+ case 5:
+ NewVolume = aHelper.AddVolume(aNds[0], aNds[1], aNds[2],
+ aNds[3], aNds[4], id, theForce3d);
+ break;
case 6:
NewVolume = aHelper.AddVolume(aNds[0], aNds[1], aNds[2],
aNds[3], aNds[4], aNds[5], id, theForce3d);
return SEW_OK;
}
+//================================================================================
/*!
\brief Creates a hole in a mesh by doubling the nodes of some particular elements
\param theElems - the list of elements (edges or faces) to be replicated
replicated nodes should be associated to.
\return TRUE if operation has been completed successfully, FALSE otherwise
*/
+//================================================================================
+
bool SMESH_MeshEditor::DoubleNodes( const TIDSortedElemSet& theElems,
const TIDSortedElemSet& theNodesNot,
const TIDSortedElemSet& theAffectedElems )
return res;
}
+//================================================================================
/*!
\brief Creates a hole in a mesh by doubling the nodes of some particular elements
\param theMeshDS - mesh instance
\param theIsDoubleElem - flag os to replicate element or modify
\return TRUE if operation has been completed successfully, FALSE otherwise
*/
+//================================================================================
+
bool SMESH_MeshEditor::doubleNodes( SMESHDS_Mesh* theMeshDS,
const TIDSortedElemSet& theElems,
const TIDSortedElemSet& theNodesNot,
theNodeNodeMap[ aCurrNode ] = aNewNode;
myLastCreatedNodes.Append( aNewNode );
}
- isDuplicate |= (aCurrNode == aNewNode);
+ isDuplicate |= (aCurrNode != aNewNode);
newNodes[ ind++ ] = aNewNode;
}
if ( !isDuplicate )
return res;
}
+//================================================================================
/*!
- \brief Check if element located inside shape
- \return TRUE if IN or ON shape, FALSE otherwise
+ \brief Creates a hole in a mesh by doubling the nodes of some particular elements
+ \param theNodes - identifiers of nodes to be doubled
+ \param theModifiedElems - identifiers of elements to be updated by the new (doubled)
+ nodes. If list of element identifiers is empty then nodes are doubled but
+ they not assigned to elements
+ \return TRUE if operation has been completed successfully, FALSE otherwise
*/
+//================================================================================
-static bool isInside(const SMDS_MeshElement* theElem,
- BRepClass3d_SolidClassifier& theBsc3d,
- const double theTol)
+bool SMESH_MeshEditor::DoubleNodes( const std::list< int >& theListOfNodes,
+ const std::list< int >& theListOfModifiedElems )
{
- gp_XYZ centerXYZ (0, 0, 0);
- SMDS_ElemIteratorPtr aNodeItr = theElem->nodesIterator();
- while (aNodeItr->more())
+ myLastCreatedElems.Clear();
+ myLastCreatedNodes.Clear();
+
+ if ( theListOfNodes.size() == 0 )
+ return false;
+
+ SMESHDS_Mesh* aMeshDS = GetMeshDS();
+ if ( !aMeshDS )
+ return false;
+
+ // iterate through nodes and duplicate them
+
+ std::map< const SMDS_MeshNode*, const SMDS_MeshNode* > anOldNodeToNewNode;
+
+ std::list< int >::const_iterator aNodeIter;
+ for ( aNodeIter = theListOfNodes.begin(); aNodeIter != theListOfNodes.end(); ++aNodeIter )
{
- SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next();
- centerXYZ += gp_XYZ(aNode->X(), aNode->Y(), aNode->Z());
+ int aCurr = *aNodeIter;
+ SMDS_MeshNode* aNode = (SMDS_MeshNode*)aMeshDS->FindNode( aCurr );
+ if ( !aNode )
+ continue;
+
+ // duplicate node
+
+ const SMDS_MeshNode* aNewNode = aMeshDS->AddNode( aNode->X(), aNode->Y(), aNode->Z() );
+ if ( aNewNode )
+ {
+ anOldNodeToNewNode[ aNode ] = aNewNode;
+ myLastCreatedNodes.Append( aNewNode );
+ }
}
- gp_Pnt aPnt(centerXYZ);
- theBsc3d.Perform(aPnt, theTol);
- TopAbs_State aState = theBsc3d.State();
- return (aState == TopAbs_IN || aState == TopAbs_ON );
+
+ // Create map of new nodes for modified elements
+
+ std::map< SMDS_MeshElement*, vector<const SMDS_MeshNode*> > anElemToNodes;
+
+ std::list< int >::const_iterator anElemIter;
+ for ( anElemIter = theListOfModifiedElems.begin();
+ anElemIter != theListOfModifiedElems.end(); ++anElemIter )
+ {
+ int aCurr = *anElemIter;
+ SMDS_MeshElement* anElem = (SMDS_MeshElement*)aMeshDS->FindElement( aCurr );
+ if ( !anElem )
+ continue;
+
+ vector<const SMDS_MeshNode*> aNodeArr( anElem->NbNodes() );
+
+ SMDS_ElemIteratorPtr anIter = anElem->nodesIterator();
+ int ind = 0;
+ while ( anIter->more() )
+ {
+ SMDS_MeshNode* aCurrNode = (SMDS_MeshNode*)anIter->next();
+ if ( aCurr && anOldNodeToNewNode.find( aCurrNode ) != anOldNodeToNewNode.end() )
+ {
+ const SMDS_MeshNode* aNewNode = anOldNodeToNewNode[ aCurrNode ];
+ aNodeArr[ ind++ ] = aNewNode;
+ }
+ else
+ aNodeArr[ ind++ ] = aCurrNode;
+ }
+ anElemToNodes[ anElem ] = aNodeArr;
+ }
+
+ // Change nodes of elements
+
+ std::map< SMDS_MeshElement*, vector<const SMDS_MeshNode*> >::iterator
+ anElemToNodesIter = anElemToNodes.begin();
+ for ( ; anElemToNodesIter != anElemToNodes.end(); ++anElemToNodesIter )
+ {
+ const SMDS_MeshElement* anElem = anElemToNodesIter->first;
+ vector<const SMDS_MeshNode*> aNodeArr = anElemToNodesIter->second;
+ if ( anElem )
+ aMeshDS->ChangeElementNodes( anElem, &aNodeArr[ 0 ], anElem->NbNodes() );
+ }
+
+ return true;
+}
+
+namespace {
+
+ //================================================================================
+ /*!
+ \brief Check if element located inside shape
+ \return TRUE if IN or ON shape, FALSE otherwise
+ */
+ //================================================================================
+
+ template<class Classifier>
+ bool isInside(const SMDS_MeshElement* theElem,
+ Classifier& theClassifier,
+ const double theTol)
+ {
+ gp_XYZ centerXYZ (0, 0, 0);
+ SMDS_ElemIteratorPtr aNodeItr = theElem->nodesIterator();
+ while (aNodeItr->more())
+ centerXYZ += TNodeXYZ(cast2Node( aNodeItr->next()));
+
+ gp_Pnt aPnt = centerXYZ / theElem->NbNodes();
+ theClassifier.Perform(aPnt, theTol);
+ TopAbs_State aState = theClassifier.State();
+ return (aState == TopAbs_IN || aState == TopAbs_ON );
+ }
+
+ //================================================================================
+ /*!
+ * \brief Classifier of the 3D point on the TopoDS_Face
+ * with interaface suitable for isInside()
+ */
+ //================================================================================
+
+ struct _FaceClassifier
+ {
+ Extrema_ExtPS _extremum;
+ BRepAdaptor_Surface _surface;
+ TopAbs_State _state;
+
+ _FaceClassifier(const TopoDS_Face& face):_extremum(),_surface(face),_state(TopAbs_OUT)
+ {
+ _extremum.Initialize( _surface,
+ _surface.FirstUParameter(), _surface.LastUParameter(),
+ _surface.FirstVParameter(), _surface.LastVParameter(),
+ _surface.Tolerance(), _surface.Tolerance() );
+ }
+ void Perform(const gp_Pnt& aPnt, double theTol)
+ {
+ _state = TopAbs_OUT;
+ _extremum.Perform(aPnt);
+ if ( _extremum.IsDone() )
+ for ( int iSol = 1; iSol <= _extremum.NbExt() && _state == TopAbs_OUT; ++iSol)
+ _state = ( _extremum.Value(iSol) <= theTol ? TopAbs_IN : TopAbs_OUT );
+ }
+ TopAbs_State State() const
+ {
+ return _state;
+ }
+ };
}
+//================================================================================
/*!
\brief Creates a hole in a mesh by doubling the nodes of some particular elements
\param theElems - group of of elements (edges or faces) to be replicated
- \param theNodesNot - group of nodes not to replicated
+ \param theNodesNot - group of nodes not to replicate
\param theShape - shape to detect affected elements (element which geometric center
located on or inside shape).
The replicated nodes should be associated to affected elements.
\return TRUE if operation has been completed successfully, FALSE otherwise
*/
+//================================================================================
bool SMESH_MeshEditor::DoubleNodesInRegion( const TIDSortedElemSet& theElems,
const TIDSortedElemSet& theNodesNot,
return false;
const double aTol = Precision::Confusion();
- BRepClass3d_SolidClassifier bsc3d(theShape);
- bsc3d.PerformInfinitePoint(aTol);
+ auto_ptr< BRepClass3d_SolidClassifier> bsc3d;
+ auto_ptr<_FaceClassifier> aFaceClassifier;
+ if ( theShape.ShapeType() == TopAbs_SOLID )
+ {
+ bsc3d.reset( new BRepClass3d_SolidClassifier(theShape));;
+ bsc3d->PerformInfinitePoint(aTol);
+ }
+ else if (theShape.ShapeType() == TopAbs_FACE )
+ {
+ aFaceClassifier.reset( new _FaceClassifier(TopoDS::Face(theShape)));
+ }
// iterates on indicated elements and get elements by back references from their nodes
TIDSortedElemSet anAffected;
SMDS_ElemIteratorPtr nodeItr = anElem->nodesIterator();
while ( nodeItr->more() )
{
- const SMDS_MeshNode* aNode = static_cast<const SMDS_MeshNode*>(nodeItr->next());
+ const SMDS_MeshNode* aNode = cast2Node(nodeItr->next());
if ( !aNode || theNodesNot.find(aNode) != theNodesNot.end() )
continue;
SMDS_ElemIteratorPtr backElemItr = aNode->GetInverseElementIterator();
while ( backElemItr->more() )
{
- SMDS_MeshElement* curElem = (SMDS_MeshElement*)backElemItr->next();
+ const SMDS_MeshElement* curElem = backElemItr->next();
if ( curElem && theElems.find(curElem) == theElems.end() &&
- isInside( curElem, bsc3d, aTol ) )
+ ( bsc3d.get() ?
+ isInside( curElem, *bsc3d, aTol ) :
+ isInside( curElem, *aFaceClassifier, aTol )))
anAffected.insert( curElem );
}
}
return DoubleNodes( theElems, theNodesNot, anAffected );
}
+//================================================================================
/*!
* \brief Generated skin mesh (containing 2D cells) from 3D mesh
* The created 2D mesh elements based on nodes of free faces of boundary volumes
* \return TRUE if operation has been completed successfully, FALSE otherwise
*/
+//================================================================================
bool SMESH_MeshEditor::Make2DMeshFrom3D()
{
{
const SMDS_MeshVolume* volume = vIt->next();
SMDS_VolumeTool vTool( volume );
+ vTool.SetExternalNormal();
const bool isPoly = volume->IsPoly();
const bool isQuad = volume->IsQuadratic();
for ( int iface = 0, n = vTool.NbFaces(); iface < n; iface++ )
vector<const SMDS_MeshNode *> nodes;
int nbFaceNodes = vTool.NbFaceNodes(iface);
const SMDS_MeshNode** faceNodes = vTool.GetFaceNodes(iface);
- if (vTool.IsFaceExternal(iface))
- {
- int inode = 0;
- for ( ; inode < nbFaceNodes; inode += isQuad ? 2 : 1)
+ int inode = 0;
+ for ( ; inode < nbFaceNodes; inode += isQuad ? 2 : 1)
+ nodes.push_back(faceNodes[inode]);
+ if (isQuad)
+ for ( inode = 1; inode < nbFaceNodes; inode += 2)
nodes.push_back(faceNodes[inode]);
- if (isQuad)
- for ( inode = 1; inode < nbFaceNodes; inode += 2)
- nodes.push_back(faceNodes[inode]);
- }
- else
- {
- int inode = nbFaceNodes-1;
- for ( ; inode >=0; inode -= isQuad ? 2 : 1)
- nodes.push_back(faceNodes[inode]);
- if (isQuad)
- for ( inode = nbFaceNodes-2; inode >=0; inode -= 2)
- nodes.push_back(faceNodes[inode]);
- }
// add new face based on volume nodes
if (aMesh->FindFace( nodes ) )
