-// Copyright (C) 2007-2021 CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2024 CEA, EDF, OPEN CASCADE
//
// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
//
#include "StdMeshers_Cartesian_3D.hxx"
#include "StdMeshers_CartesianParameters3D.hxx"
-
-#include "ObjectPool.hxx"
-#include "SMDS_MeshNode.hxx"
-#include "SMDS_VolumeTool.hxx"
-#include "SMESHDS_Mesh.hxx"
-#include "SMESH_Block.hxx"
-#include "SMESH_Comment.hxx"
-#include "SMESH_ControlsDef.hxx"
-#include "SMESH_Mesh.hxx"
-#include "SMESH_MeshAlgos.hxx"
-#include "SMESH_MeshEditor.hxx"
-#include "SMESH_MesherHelper.hxx"
-#include "SMESH_subMesh.hxx"
-#include "SMESH_subMeshEventListener.hxx"
+#include "StdMeshers_Cartesian_VL.hxx"
#include "StdMeshers_FaceSide.hxx"
+#include "StdMeshers_ViscousLayers.hxx"
+
+#include <ObjectPool.hxx>
+#include <SMDS_LinearEdge.hxx>
+#include <SMDS_MeshNode.hxx>
+#include <SMDS_VolumeOfNodes.hxx>
+#include <SMDS_VolumeTool.hxx>
+#include <SMESHDS_Mesh.hxx>
+#include <SMESH_Block.hxx>
+#include <SMESH_Comment.hxx>
+#include <SMESH_ControlsDef.hxx>
+#include <SMESH_Mesh.hxx>
+#include <SMESH_MeshAlgos.hxx>
+#include <SMESH_MeshEditor.hxx>
+#include <SMESH_MesherHelper.hxx>
+#include <SMESH_subMesh.hxx>
+#include <SMESH_subMeshEventListener.hxx>
#include <utilities.h>
#include <Utils_ExceptHandlers.hxx>
#include <gp_Sphere.hxx>
#include <gp_Torus.hxx>
+//STD
#include <limits>
+#include <mutex>
+#include <thread>
#include <boost/container/flat_map.hpp>
#define WINVER 0x0A00
#define _WIN32_WINNT 0x0A00
#endif
-
+#include <algorithm>
#include <tbb/parallel_for.h>
-//#include <tbb/enumerable_thread_specific.h>
#endif
using namespace std;
using namespace SMESH;
+std::mutex _eMutex;
+std::mutex _bMutex;
//=============================================================================
/*!
{
_name = "Cartesian_3D";
_shapeType = (1 << TopAbs_SOLID); // 1 bit /shape type
- _compatibleHypothesis.push_back("CartesianParameters3D");
+ _compatibleHypothesis.push_back( "CartesianParameters3D" );
+ _compatibleHypothesis.push_back( StdMeshers_ViscousLayers::GetHypType() );
_onlyUnaryInput = false; // to mesh all SOLIDs at once
_requireDiscreteBoundary = false; // 2D mesh not needed
{
aStatus = SMESH_Hypothesis::HYP_MISSING;
- const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape);
+ const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape, /*skipAux=*/false);
list <const SMESHDS_Hypothesis* >::const_iterator h = hyps.begin();
if ( h == hyps.end())
{
return false;
}
+ _hyp = nullptr;
+ _hypViscousLayers = nullptr;
+ _isComputeOffset = false;
+
for ( ; h != hyps.end(); ++h )
{
- if (( _hyp = dynamic_cast<const StdMeshers_CartesianParameters3D*>( *h )))
+ if ( !_hyp && ( _hyp = dynamic_cast<const StdMeshers_CartesianParameters3D*>( *h )))
{
aStatus = _hyp->IsDefined() ? HYP_OK : HYP_BAD_PARAMETER;
- break;
+ }
+ else
+ {
+ _hypViscousLayers = dynamic_cast<const StdMeshers_ViscousLayers*>( *h );
}
}
namespace
{
- typedef int TGeomID; // IDs of sub-shapes
+ /*!
+ * \brief Temporary mesh to hold
+ */
+ struct TmpMesh: public SMESH_Mesh
+ {
+ TmpMesh() {
+ _isShapeToMesh = (_id = 0);
+ _meshDS = new SMESHDS_Mesh( _id, true );
+ }
+ };
+
+ typedef int TGeomID; // IDs of sub-shapes
+ typedef TopTools_ShapeMapHasher TShapeHasher; // non-oriented shape hasher
+ typedef std::array< int, 3 > TIJK;
const TGeomID theUndefID = 1e+9;
mutable vector< TGeomID > _faceIDs;
B_IntersectPoint(): _node(NULL) {}
- bool Add( const vector< TGeomID >& fIDs, const SMDS_MeshNode* n=0 ) const;
+ bool Add( const vector< TGeomID >& fIDs, const SMDS_MeshNode* n=NULL ) const;
TGeomID HasCommonFace( const B_IntersectPoint * other, TGeomID avoidFace=-1 ) const;
size_t GetCommonFaces( const B_IntersectPoint * other, TGeomID * commonFaces ) const;
bool IsOnFace( TGeomID faceID ) const;
{
_curInd[0] = i; _curInd[1] = j; _curInd[2] = k;
}
+ void SetLineIndex(size_t i)
+ {
+ _curInd[_iVar2] = i / _size[_iVar1];
+ _curInd[_iVar1] = i % _size[_iVar1];
+ }
void operator++()
{
if ( ++_curInd[_iVar1] == _size[_iVar1] )
size_t LineIndex01 () const { return _curInd[_iVar1] + (_curInd[_iVar2] + 1 )* _size[_iVar1]; }
size_t LineIndex11 () const { return (_curInd[_iVar1] + 1 ) + (_curInd[_iVar2] + 1 )* _size[_iVar1]; }
void SetIndexOnLine (size_t i) { _curInd[ _iConst ] = i; }
+ bool IsValidIndexOnLine (size_t i) const { return i < _size[ _iConst ]; }
size_t NbLines() const { return _size[_iVar1] * _size[_iVar2]; }
};
+ struct FaceGridIntersector;
// --------------------------------------------------------------------------
/*!
* \brief Container of GridLine's
// index shift within _nodes of nodes of a cell from the 1st node
int _nodeShift[8];
- vector< const SMDS_MeshNode* > _nodes; // mesh nodes at grid nodes
+ vector< const SMDS_MeshNode* > _nodes; // mesh nodes at grid nodes
+ vector< const SMDS_MeshNode* > _allBorderNodes; // mesh nodes between the bounding box and the geometry boundary
+
vector< const F_IntersectPoint* > _gridIntP; // grid node intersection with geometry
ObjectPool< E_IntersectPoint > _edgeIntPool; // intersections with EDGEs
ObjectPool< F_IntersectPoint > _extIntPool; // intersections with extended INTERNAL FACEs
bool _toConsiderInternalFaces;
bool _toUseThresholdForInternalFaces;
double _sizeThreshold;
+ bool _toUseQuanta;
+ double _quanta;
SMESH_MesherHelper* _helper;
{
return i + j*_coords[0].size() + k*_coords[0].size()*_coords[1].size();
}
+ size_t NodeIndex( const TIJK& ijk ) const
+ {
+ return NodeIndex( ijk[0], ijk[1], ijk[2] );
+ }
size_t NodeIndexDX() const { return 1; }
size_t NodeIndexDY() const { return _coords[0].size(); }
size_t NodeIndexDZ() const { return _coords[0].size() * _coords[1].size(); }
LineIndexer GetLineIndexer(size_t iDir) const;
+ size_t GetLineDir( const GridLine* line, size_t & index ) const;
E_IntersectPoint* Add( const E_IntersectPoint& ip )
{
// --------------------------------------------------------------------------------
struct _Node //!< node either at a hexahedron corner or at intersection
{
- const SMDS_MeshNode* _node; // mesh node at hexahedron corner
+ const SMDS_MeshNode* _node; // mesh node at hexahedron corner
+ const SMDS_MeshNode* _boundaryCornerNode; // missing mesh node due to hex truncation on the boundary
const B_IntersectPoint* _intPoint;
const _Face* _usedInFace;
char _isInternalFlags;
:_node(n), _intPoint(ip), _usedInFace(0), _isInternalFlags(0) {}
const SMDS_MeshNode* Node() const
{ return ( _intPoint && _intPoint->_node ) ? _intPoint->_node : _node; }
+ const SMDS_MeshNode* BoundaryNode() const
+ { return _node ? _node : _boundaryCornerNode; }
const E_IntersectPoint* EdgeIntPnt() const
{ return static_cast< const E_IntersectPoint* >( _intPoint ); }
const F_IntersectPoint* FaceIntPnt() const
}
void Add( const E_IntersectPoint* ip )
{
+ const std::lock_guard<std::mutex> lock(_eMutex);
// Possible cases before Add(ip):
/// 1) _node != 0 --> _Node at hex corner ( _intPoint == 0 || _intPoint._node == 0 )
/// 2) _node == 0 && _intPoint._node != 0 --> link intersected by FACE
TGeomID _solidID;
double _size;
const SMDS_MeshElement* _volume; // new volume
+ std::vector<const SMDS_MeshElement*> _brotherVolume; // produced due to poly split
vector< SMESH_Block::TShapeID > _names; // name of side a polygon originates from
int _origNodeInd; // index of _hexNodes[0] node within the _grid
size_t _i,_j,_k;
bool _hasTooSmall;
-
-#ifdef _DEBUG_
int _cellID;
-#endif
public:
Hexahedron(Grid* grid);
bool isInHole() const;
bool hasStrangeEdge() const;
bool checkPolyhedronSize( bool isCutByInternalFace, double & volSize ) const;
+ int checkPolyhedronValidity( _volumeDef* volDef, std::vector<std::vector<int>>& splitQuantities,
+ std::vector<std::vector<const SMDS_MeshNode*>>& splitNodes );
+ const SMDS_MeshElement* addPolyhedronToMesh( _volumeDef* volDef, SMESH_MesherHelper& helper, const std::vector<const SMDS_MeshNode*>& nodes,
+ const std::vector<int>& quantities );
bool addHexa ();
bool addTetra();
bool addPenta();
if ( solids.size() == 2 )
{
if ( solids == solidsBef )
- return theUndefID; //solids.contain( prevID ) ? solids.otherThan( prevID ) : theUndefID;
+ return solids.contain( prevID ) ? solids.otherThan( prevID ) : theUndefID; // bos #29212
}
return solids.oneCommon( solidsBef );
}
bool B_IntersectPoint::Add( const vector< TGeomID >& fIDs,
const SMDS_MeshNode* n) const
{
+ const std::lock_guard<std::mutex> lock(_bMutex);
size_t prevNbF = _faceIDs.size();
if ( _faceIDs.empty() )
if ( it == _faceIDs.end() )
_faceIDs.push_back( fIDs[i] );
}
- if ( !_node )
+ if ( !_node && n != NULL )
_node = n;
return prevNbF < _faceIDs.size();
s[indices[iDir*3]], s[indices[iDir*3+1]], s[indices[iDir*3+2]]);
return li;
}
+ //================================================================================
+ /*
+ * Return direction [0,1,2] of a GridLine
+ */
+ size_t Grid::GetLineDir( const GridLine* line, size_t & index ) const
+ {
+ for ( size_t iDir = 0; iDir < 3; ++iDir )
+ if ( &_lines[ iDir ][0] <= line && line <= &_lines[ iDir ].back() )
+ {
+ index = line - &_lines[ iDir ][0];
+ return iDir;
+ }
+ return -1;
+ }
//=============================================================================
/*
* Creates GridLine's of the grid
const size_t nbGridNodes = _coords[0].size() * _coords[1].size() * _coords[2].size();
vector< TGeomID > shapeIDVec( nbGridNodes, theUndefID );
_nodes.resize( nbGridNodes, 0 );
+ _allBorderNodes.resize( nbGridNodes, 0 );
_gridIntP.resize( nbGridNodes, NULL );
SMESHDS_Mesh* mesh = helper.GetMeshDS();
if ( ++nodeCoord < coordEnd )
nodeParam = *nodeCoord - *coord0;
else
- break;
+ break;
}
if ( nodeCoord == coordEnd ) break;
}
+
// create a mesh node on a GridLine at ip if it does not coincide with a grid node
if ( nodeParam > ip->_paramOnLine + _tol )
{
SetOnShape( _nodes[ nodeIndex ], *_gridIntP[ nodeIndex ], & v );
UpdateFacesOfVertex( *_gridIntP[ nodeIndex ], v );
}
+ else if ( _toUseQuanta && !_allBorderNodes[ nodeIndex ] /*add all nodes outside the body. Used to reconstruct the hexahedrals when polys are not desired!*/)
+ {
+ gp_XYZ xyz = ( _coords[0][x] * _axes[0] +
+ _coords[1][y] * _axes[1] +
+ _coords[2][z] * _axes[2] );
+ _allBorderNodes[ nodeIndex ] = mesh->AddNode( xyz.X(), xyz.Y(), xyz.Z() );
+ mesh->SetNodeInVolume( _allBorderNodes[ nodeIndex ], shapeIDVec[ nodeIndex ]);
+ }
}
#ifdef _MY_DEBUG_
{
if ( intPnts.begin()->_transition != Trans_TANGENT &&
intPnts.begin()->_transition != Trans_APEX )
- throw SMESH_ComputeError (COMPERR_ALGO_FAILED,
- SMESH_Comment("Wrong SOLE transition of GridLine (")
- << li._curInd[li._iVar1] << ", " << li._curInd[li._iVar2]
- << ") along " << li._nameConst
- << ": " << trName[ intPnts.begin()->_transition] );
+ throw SMESH_ComputeError (COMPERR_ALGO_FAILED,
+ SMESH_Comment("Wrong SOLE transition of GridLine (")
+ << li._curInd[li._iVar1] << ", " << li._curInd[li._iVar2]
+ << ") along " << li._nameConst
+ << ": " << trName[ intPnts.begin()->_transition] );
}
else
{
SMESH_Comment("Wrong END transition of GridLine (")
<< li._curInd[li._iVar1] << ", " << li._curInd[li._iVar2]
<< ") along " << li._nameConst
- << ": " << trName[ intPnts.rbegin()->_transition ]);
+ << ": " << trName[ intPnts.rbegin()->_transition ]);
}
}
}
#endif
+
+ return;
}
//=============================================================================
tgtLink._link = _hexLinks + ( srcLink._link - other._hexLinks );
}
}
-#ifdef _DEBUG_
- _cellID = cellID;
-#else
- (void)cellID; // unused in release mode
-#endif
+
+ if (SALOME::VerbosityActivated())
+ _cellID = cellID;
}
//================================================================================
{
_hexNodes[iN]._isInternalFlags = 0;
+ // Grid node
_hexNodes[iN]._node = _grid->_nodes [ _origNodeInd + _grid->_nodeShift[iN] ];
_hexNodes[iN]._intPoint = _grid->_gridIntP[ _origNodeInd + _grid->_nodeShift[iN] ];
+ if ( _grid->_allBorderNodes[ _origNodeInd + _grid->_nodeShift[iN] ] )
+ _hexNodes[iN]._boundaryCornerNode = _grid->_allBorderNodes [ _origNodeInd + _grid->_nodeShift[iN] ];
+
if ( _hexNodes[iN]._node && !solid->Contains( _hexNodes[iN]._node->GetShapeID() ))
_hexNodes[iN]._node = 0;
+
if ( _hexNodes[iN]._intPoint && !solid->ContainsAny( _hexNodes[iN]._intPoint->_faceIDs ))
_hexNodes[iN]._intPoint = 0;
if ( _nbFaceIntNodes + _eIntPoints.size() > 0 &&
_nbFaceIntNodes + _eIntPoints.size() + _nbCornerNodes > 3)
{
- _intNodes.reserve( 3 * _nbBndNodes + _nbFaceIntNodes + _eIntPoints.size() );
+ _intNodes.reserve( 3 * ( _nbBndNodes + _nbFaceIntNodes + _eIntPoints.size() ));
// this method can be called in parallel, so use own helper
SMESH_MesherHelper helper( *_grid->_helper->GetMesh() );
// 1) add this->_eIntPoints to _Face::_eIntNodes
// 2) fill _intNodes and _vIntNodes
//
- const double tol2 = _grid->_tol * _grid->_tol;
+ const double tol2 = _grid->_tol * _grid->_tol * 4;
int facets[3], nbFacets, subEntity;
for ( int iF = 0; iF < 6; ++iF )
solid = _grid->GetSolid();
if ( !_grid->_geometry.IsOneSolid() )
{
- TGeomID solidIDs[20];
+ TGeomID solidIDs[20] = { 0 };
size_t nbSolids = getSolids( solidIDs );
if ( nbSolids > 1 )
{
_volumeDefs._nodes.clear();
_volumeDefs._quantities.clear();
_volumeDefs._names.clear();
-
// create a classic cell if possible
int nbPolygons = 0;
return !_volumeDefs._nodes.empty();
}
+
+ template<typename Type>
+ void computeHexa(Type& hex)
+ {
+ if ( hex )
+ hex->computeElements();
+ }
+
+ // Implement parallel computation of Hexa with c++ thread implementation
+ template<typename Iterator, class Function>
+ void parallel_for(const Iterator& first, const Iterator& last, Function&& f, const int nthreads = 1)
+ {
+ const unsigned int group = ((last-first))/std::abs(nthreads);
+
+ std::vector<std::thread> threads;
+ threads.reserve(nthreads);
+ Iterator it = first;
+ for (; it < last-group; it += group) {
+ // to create a thread
+ // Pass iterators by value and the function by reference!
+ auto lambda = [=,&f](){ std::for_each(it, std::min(it+group, last), f);};
+
+ // stack the threads
+ threads.push_back( std::thread( lambda ) );
+ }
+
+ std::for_each(it, last, f); // last steps while we wait for other threads
+ std::for_each(threads.begin(), threads.end(), [](std::thread& x){x.join();});
+ }
//================================================================================
/*!
* \brief Create elements in the mesh
// compute definitions of volumes resulted from hexadron intersection
#ifdef WITH_TBB
- tbb::parallel_for ( tbb::blocked_range<size_t>( 0, intHexa.size() ),
- ParallelHexahedron( intHexa ),
- tbb::simple_partitioner()); // computeElements() is called here
+ auto numOfThreads = std::thread::hardware_concurrency();
+ numOfThreads = (numOfThreads != 0) ? numOfThreads : 1;
+ parallel_for(intHexa.begin(), intHexa.end(), computeHexa<Hexahedron*>, numOfThreads );
#else
for ( size_t i = 0; i < intHexa.size(); ++i )
if ( Hexahedron * hex = intHexa[ i ] )
h->_eIntPoints.reserve(2);
h->_eIntPoints.push_back( ip );
added = true;
-#ifdef _DEBUG_
+
// check if ip is really inside the hex
- if ( h->isOutParam( ip->_uvw ))
+ if (SALOME::VerbosityActivated() && h->isOutParam( ip->_uvw ))
throw SALOME_Exception("ip outside a hex");
-#endif
}
}
return added;
{
F_IntersectPoint noIntPnt;
const bool toCheckNodePos = _grid->IsToCheckNodePos();
+ const bool useQuanta = _grid->_toUseQuanta;
int nbAdded = 0;
// add elements resulted from hexahedron intersection
}
} // loop to get nodes
- const SMDS_MeshElement* v = 0;
+ const SMDS_MeshElement* v = 0;
if ( !volDef->_quantities.empty() )
- {
- v = helper.AddPolyhedralVolume( nodes, volDef->_quantities );
- volDef->_size = SMDS_VolumeTool( v ).GetSize();
- if ( volDef->_size < 0 ) // invalid polyhedron
+ {
+ if ( !useQuanta )
{
- if ( ! SMESH_MeshEditor( helper.GetMesh() ).Reorient( v ) || // try to fix
- SMDS_VolumeTool( v ).GetSize() < 0 )
+ // split polyhedrons of with disjoint volumes
+ std::vector<std::vector<int>> splitQuantities;
+ std::vector<std::vector< const SMDS_MeshNode* > > splitNodes;
+ if ( checkPolyhedronValidity( volDef, splitQuantities, splitNodes ) == 1 )
+ v = addPolyhedronToMesh( volDef, helper, nodes, volDef->_quantities );
+ else
{
- helper.GetMeshDS()->RemoveFreeElement( v, /*sm=*/nullptr, /*fromGroups=*/false );
- v = nullptr;
- //_hasTooSmall = true;
-#ifdef _DEBUG_
- std::cout << "Remove INVALID polyhedron, _cellID = " << _cellID
- << " ijk = ( " << _i << " " << _j << " " << _k << " ) "
- << " solid " << volDef->_solidID << std::endl;
-#endif
+ int counter = -1;
+ for (size_t id = 0; id < splitQuantities.size(); id++)
+ {
+ v = addPolyhedronToMesh( volDef, helper, splitNodes[ id ], splitQuantities[ id ] );
+ if ( id < splitQuantities.size()-1 )
+ volDef->_brotherVolume.push_back( v );
+ counter++;
+ }
+ nbAdded += counter;
}
}
+ else
+ {
+ const double quanta = _grid->_quanta;
+ double polyVol = volDef->_size;
+ double hexaVolume = _sideLength[0] * _sideLength[1] * _sideLength[2];
+ if ( hexaVolume > 0.0 && polyVol/hexaVolume >= quanta /*set the volume if the relation is satisfied*/)
+ v = helper.AddVolume( _hexNodes[0].BoundaryNode(), _hexNodes[2].BoundaryNode(),
+ _hexNodes[3].BoundaryNode(), _hexNodes[1].BoundaryNode(),
+ _hexNodes[4].BoundaryNode(), _hexNodes[6].BoundaryNode(),
+ _hexNodes[7].BoundaryNode(), _hexNodes[5].BoundaryNode() );
+
+ }
}
else
{
if ( volDef->_volume )
{
helper.GetMeshDS()->SetMeshElementOnShape( volDef->_volume, volDef->_solidID );
+ for (auto broVol : volDef->_brotherVolume )
+ {
+ helper.GetMeshDS()->SetMeshElementOnShape( broVol, volDef->_solidID );
+ }
}
}
//================================================================================
/*!
* \brief Return true if the element is in a hole
+ * \remark consider a cell to be in a hole if all links in any direction
+ * comes OUT of geometry
*/
bool Hexahedron::isInHole() const
{
return volume > initVolume / _grid->_sizeThreshold;
}
+
+ //================================================================================
+ /*!
+ * \brief Check that all faces in polyhedron are connected so a unique volume is defined.
+ * We test that it is possible to go from any node to all nodes in the polyhedron.
+ * The set of nodes that can be visit within then defines a unique element.
+ * In case more than one polyhedron is detected. The function return the set of quantities and nodes defining separates elements.
+ * Reference to issue #bos[38521][EDF] Generate polyhedron with separate volume.
+ */
+ int Hexahedron::checkPolyhedronValidity( _volumeDef* volDef, std::vector<std::vector<int>>& splitQuantities,
+ std::vector<std::vector<const SMDS_MeshNode*>>& splitNodes )
+ {
+ int mySet = 1;
+ std::map<int,int> numberOfSets; // define set id with the number of faces associated!
+ if ( !volDef->_quantities.empty() )
+ {
+ auto connectivity = volDef->_quantities;
+ int accum = 0;
+ std::vector<bool> allFaces( connectivity.size(), false );
+ std::set<int> elementSet;
+ allFaces[ 0 ] = true; // the first node below to the first face
+ size_t connectedFaces = 1;
+ // Start filling the set with the nodes of the first face
+ splitQuantities.push_back( { connectivity[ 0 ] } );
+ splitNodes.push_back( { volDef->_nodes[ 0 ].Node() } );
+ elementSet.insert( volDef->_nodes[ 0 ].Node()->GetID() );
+ for (int n = 1; n < connectivity[ 0 ]; n++)
+ {
+ elementSet.insert( volDef->_nodes[ n ].Node()->GetID() );
+ splitNodes.back().push_back( volDef->_nodes[ n ].Node() );
+ }
+
+ numberOfSets.insert( std::pair<int,int>(mySet,1) );
+ while ( connectedFaces != allFaces.size() )
+ {
+ for (size_t innerId = 1; innerId < connectivity.size(); innerId++)
+ {
+ if ( innerId == 1 )
+ accum = connectivity[ 0 ];
+
+ if ( !allFaces[ innerId ] )
+ {
+ int faceCounter = 0;
+ for (int n = 0; n < connectivity[ innerId ]; n++)
+ {
+ int nodeId = volDef->_nodes[ accum + n ].Node()->GetID();
+ if ( elementSet.count( nodeId ) != 0 )
+ faceCounter++;
+ }
+ if ( faceCounter >= 2 ) // found coincidences nodes
+ {
+ for (int n = 0; n < connectivity[ innerId ]; n++)
+ {
+ int nodeId = volDef->_nodes[ accum + n ].Node()->GetID();
+ // insert new nodes so other faces can be identified as belowing to the element
+ splitNodes.back().push_back( volDef->_nodes[ accum + n ].Node() );
+ elementSet.insert( nodeId );
+ }
+ allFaces[ innerId ] = true;
+ splitQuantities.back().push_back( connectivity[ innerId ] );
+ numberOfSets[ mySet ]++;
+ connectedFaces++;
+ innerId = 0; // to restart searching!
+ }
+ }
+ accum += connectivity[ innerId ];
+ }
+
+ if ( connectedFaces != allFaces.size() )
+ {
+ // empty the set, and fill it with nodes of a unvisited face!
+ elementSet.clear();
+ accum = connectivity[ 0 ];
+ for (size_t faceId = 1; faceId < connectivity.size(); faceId++)
+ {
+ if ( !allFaces[ faceId ] )
+ {
+ splitNodes.push_back( { volDef->_nodes[ accum ].Node() } );
+ elementSet.insert( volDef->_nodes[ accum ].Node()->GetID() );
+ for (int n = 1; n < connectivity[ faceId ]; n++)
+ {
+ elementSet.insert( volDef->_nodes[ accum + n ].Node()->GetID() );
+ splitNodes.back().push_back( volDef->_nodes[ accum + n ].Node() );
+ }
+
+ splitQuantities.push_back( { connectivity[ faceId ] } );
+ allFaces[ faceId ] = true;
+ connectedFaces++;
+ break;
+ }
+ accum += connectivity[ faceId ];
+ }
+ mySet++;
+ numberOfSets.insert( std::pair<int,int>(mySet,1) );
+ }
+ }
+
+ if ( numberOfSets.size() > 1 )
+ {
+ bool allMoreThan2Faces = true;
+ for( auto k : numberOfSets )
+ {
+ if ( k.second <= 2 )
+ allMoreThan2Faces &= false;
+ }
+
+ if ( allMoreThan2Faces )
+ {
+ // The separate objects are suspect to be closed
+ return numberOfSets.size();
+ }
+ else
+ {
+ // Have to index the last face nodes to the final set
+ // contrary case return as it were a valid polyhedron for backward compatibility
+ return 1;
+ }
+ }
+ }
+ return numberOfSets.size();
+ }
+
+
+ //================================================================================
+ /*!
+ * \brief add original or separated polyhedrons to the mesh
+ */
+ const SMDS_MeshElement* Hexahedron::addPolyhedronToMesh( _volumeDef* volDef, SMESH_MesherHelper& helper, const std::vector<const SMDS_MeshNode*>& nodes,
+ const std::vector<int>& quantities )
+ {
+ const SMDS_MeshElement* v = helper.AddPolyhedralVolume( nodes, quantities );
+
+ volDef->_size = SMDS_VolumeTool( v ).GetSize();
+ if ( volDef->_size < 0 ) // invalid polyhedron
+ {
+ if ( ! SMESH_MeshEditor( helper.GetMesh() ).Reorient( v ) || // try to fix
+ SMDS_VolumeTool( v ).GetSize() < 0 )
+ {
+ helper.GetMeshDS()->RemoveFreeElement( v, /*sm=*/nullptr, /*fromGroups=*/false );
+ v = nullptr;
+ //_hasTooSmall = true;
+
+ if (SALOME::VerbosityActivated())
+ {
+ std::cout << "Remove INVALID polyhedron, _cellID = " << _cellID
+ << " ijk = ( " << _i << " " << _j << " " << _k << " ) "
+ << " solid " << volDef->_solidID << std::endl;
+ }
+ }
+ }
+ return v;
+ }
+
//================================================================================
/*!
* \brief Tries to create a hexahedron
*/
bool Hexahedron::debugDumpLink( Hexahedron::_Link* link )
{
-#ifdef _DEBUG_
- gp_Pnt p1 = link->_nodes[0]->Point(), p2 = link->_nodes[1]->Point();
- cout << "BUG: not shared link. IKJ = ( "<< _i << " " << _j << " " << _k << " )" << endl
- << "n1 (" << p1.X() << ", "<< p1.Y() << ", "<< p1.Z() << " )" << endl
- << "n2 (" << p2.X() << ", "<< p2.Y() << ", "<< p2.Z() << " )" << endl;
-#else
- (void)link; // unused in release mode
-#endif
+ if (SALOME::VerbosityActivated())
+ {
+ gp_Pnt p1 = link->_nodes[0]->Point(), p2 = link->_nodes[1]->Point();
+ cout << "BUG: not shared link. IKJ = ( "<< _i << " " << _j << " " << _k << " )" << endl
+ << "n1 (" << p1.X() << ", "<< p1.Y() << ", "<< p1.Z() << " )" << endl
+ << "n2 (" << p2.X() << ", "<< p2.Y() << ", "<< p2.Z() << " )" << endl;
+ }
+
return false;
}
//================================================================================
SMESH_MeshEditor::ElemFeatures face( SMDSAbs_Face );
SMESHDS_Mesh* meshDS = helper.GetMeshDS();
+ bool isQuantaSet = _grid->_toUseQuanta;
// check if there are internal or shared FACEs
bool hasInternal = ( !_grid->_geometry.IsOneSolid() ||
- _grid->_geometry._soleSolid.HasInternalFaces() );
+ _grid->_geometry._soleSolid.HasInternalFaces() );
for ( size_t iV = 0; iV < boundaryVolumes.size(); ++iV )
{
if ( !vTool.Set( boundaryVolumes[ iV ]))
continue;
-
TGeomID solidID = vTool.Element()->GetShapeID();
Solid * solid = _grid->GetOneOfSolids( solidID );
// find boundary facets
-
bndFacets.clear();
for ( int iF = 0, n = vTool.NbFaces(); iF < n; iF++ )
{
const SMDS_MeshElement* otherVol;
- bool isBoundary = vTool.IsFreeFace( iF, &otherVol );
+ bool isBoundary = isQuantaSet ? vTool.IsFreeFaceCheckAllNodes( iF, &otherVol ) : vTool.IsFreeFace( iF, &otherVol );
if ( isBoundary )
{
bndFacets.push_back( iF );
continue;
// create faces
-
if ( !vTool.IsPoly() )
vTool.SetExternalNormal();
for ( size_t i = 0; i < bndFacets.size(); ++i ) // loop on boundary facets
if ( nn[ iN ]->GetPosition()->GetDim() == 2 )
faceID = nn[ iN ]->GetShapeID();
}
- if ( faceID == 0 )
+ if ( faceID == 0 && !isQuantaSet /*if quanta is set boundary nodes at boundary does not coincide with any geometrical face */ )
faceID = findCommonFace( face.myNodes, helper.GetMesh() );
bool toCheckFace = faceID && (( !isBoundary ) ||
// if ( !faceID && !isBoundary )
// continue;
}
- if ( !faceID && !isBoundary )
+ if ( !faceID && !isBoundary && !isQuantaSet )
continue;
}
continue;
gp_Dir direction(1,0,0);
- const SMDS_MeshElement* anyFace = *facesToOrient.begin();
- editor.Reorient2D( facesToOrient, direction, anyFace );
+ TIDSortedElemSet refFaces;
+ editor.Reorient2D( facesToOrient, direction, refFaces, /*allowNonManifold=*/true );
}
}
return;
for ( size_t iN = 0; iN < volDef->_nodes.size(); ++iN )
volDef->_nodes[iN].Node()->setIsMarked( false );
}
+ if ( volDef->_brotherVolume.size() > 0 )
+ {
+ for (auto _bro : volDef->_brotherVolume )
+ {
+ _bro->setIsMarked( true );
+ boundaryElems.push_back( _bro );
+ }
+ }
}
}
bool StdMeshers_Cartesian_3D::Compute(SMESH_Mesh & theMesh,
const TopoDS_Shape & theShape)
{
+ if ( _hypViscousLayers )
+ {
+ const StdMeshers_ViscousLayers* hypViscousLayers = _hypViscousLayers;
+ _hypViscousLayers = nullptr;
+
+ StdMeshers_Cartesian_VL::ViscousBuilder builder( hypViscousLayers, theMesh, theShape );
+
+ std::string error;
+ TopoDS_Shape offsetShape = builder.MakeOffsetShape( theShape, theMesh, error );
+ if ( offsetShape.IsNull() )
+ throw SALOME_Exception( error );
+
+ SMESH_Mesh* offsetMesh = new TmpMesh();
+ offsetMesh->ShapeToMesh( offsetShape );
+ offsetMesh->GetSubMesh( offsetShape )->DependsOn();
+
+ this->_isComputeOffset = true;
+ if ( ! this->Compute( *offsetMesh, offsetShape ))
+ return false;
+
+ return builder.MakeViscousLayers( *offsetMesh, theMesh, theShape );
+ }
+
// The algorithm generates the mesh in following steps:
// 1) Intersection of grid lines with the geometry boundary.
grid._toConsiderInternalFaces = _hyp->GetToConsiderInternalFaces();
grid._toUseThresholdForInternalFaces = _hyp->GetToUseThresholdForInternalFaces();
grid._sizeThreshold = _hyp->GetSizeThreshold();
+ grid._toUseQuanta = _hyp->GetToUseQuanta();
+ grid._quanta = _hyp->GetQuanta();
+ if ( _isComputeOffset )
+ {
+ grid._toAddEdges = true;
+ grid._toCreateFaces = true;
+ }
grid.InitGeometry( theShape );
vector< TopoDS_Shape > faceVec;
grid._nodes[i]->setIsMarked( true );
}
+ for ( size_t i = 0; i < grid._allBorderNodes.size(); ++i )
+ if ( grid._allBorderNodes[i] &&
+ !grid._allBorderNodes[i]->IsNull() &&
+ grid._allBorderNodes[i]->NbInverseElements() == 0 )
+ {
+ nodesToRemove.push_back( grid._allBorderNodes[i] );
+ grid._allBorderNodes[i]->setIsMarked( true );
+ }
+
// do remove
for ( size_t i = 0; i < nodesToRemove.size(); ++i )
meshDS->RemoveFreeNode( nodesToRemove[i], /*smD=*/0, /*fromGroups=*/false );
