From d5d777a7f913b46b8068cb79a41b5ad0b97c3f7a Mon Sep 17 00:00:00 2001 From: eap Date: Fri, 17 Jun 2022 14:23:44 +0300 Subject: [PATCH] bos #17015 [CEA 17008] Body fitting with Viscous Layers for CFD meshing --- doc/gui/input/cartesian_algo.rst | 18 +- resources/StdMeshers.xml.in | 2 + src/SMDS/SMDS_LinearEdge.cxx | 7 + src/SMDS/SMDS_LinearEdge.hxx | 3 +- src/SMESH/SMESH_MeshEditor.cxx | 3 + src/SMESHDS/SMESHDS_Mesh.cxx | 12 + src/SMESHDS/SMESHDS_Mesh.hxx | 1 + src/SMESHUtils/SMESH_MeshAlgos.cxx | 119 +++ src/SMESHUtils/SMESH_MeshAlgos.hxx | 7 + src/SMESHUtils/SMESH_TypeDefs.hxx | 1 + src/StdMeshers/CMakeLists.txt | 2 + src/StdMeshers/StdMeshers_Cartesian_3D.cxx | 120 ++- src/StdMeshers/StdMeshers_Cartesian_3D.hxx | 4 + src/StdMeshers/StdMeshers_Cartesian_VL.cxx | 880 ++++++++++++++++++++ src/StdMeshers/StdMeshers_Cartesian_VL.hxx | 62 ++ src/StdMeshers/StdMeshers_ViscousLayers.cxx | 20 +- src/StdMeshers/StdMeshers_ViscousLayers.hxx | 5 + 17 files changed, 1227 insertions(+), 39 deletions(-) create mode 100644 src/StdMeshers/StdMeshers_Cartesian_VL.cxx create mode 100644 src/StdMeshers/StdMeshers_Cartesian_VL.hxx diff --git a/doc/gui/input/cartesian_algo.rst b/doc/gui/input/cartesian_algo.rst index 76b150e7e..7cf7695a0 100644 --- a/doc/gui/input/cartesian_algo.rst +++ b/doc/gui/input/cartesian_algo.rst @@ -7,7 +7,9 @@ Body Fitting 3D meshing algorithm Body Fitting algorithm generates hexahedrons of a Cartesian grid in the internal part of geometry and polyhedrons and other types of elements at the intersection of Cartesian cells with the geometrical -boundary. +boundary. The algorithm supports construction of +:ref:`viscous boundary layers` (use +:ref:`Viscous Layers hypothesis ` to define them). .. image:: ../images/cartesian3D_sphere.png :align: center @@ -15,8 +17,6 @@ boundary. .. centered:: A sphere meshed by Body Fitting algorithm -.. note:: The algorithm creates only 3D elements. To add 2D elements use :ref:`Generate boundary elements ` operation. - The meshing algorithm is as follows. #. Lines of a Cartesian structured grid defined by :ref:`Body Fitting Parameters ` hypothesis are intersected with the geometry boundary, thus nodes lying on the boundary are found. This step also allows finding out for each node of the Cartesian grid if it is inside or outside the geometry. @@ -27,7 +27,17 @@ The meshing algorithm is as follows. * add a hexahedron in the mesh, if all nodes are inside * add a polyhedron or another cell type in the mesh, if some nodes are inside and some outside. -To apply this algorithm when you define your mesh, select **Body Fitting** in the list of 3D algorithms and add **Body Fitting Parameters** hypothesis. The following dialog will appear: +.. _cartesian_VL_anchor: + +**Viscous boundary layers** are constructed as follows: + + * create an offset geometry with offset value equal to full layers thickness by using BRepOffset_MakeOffset OCCT class; + * mesh the offset geometry with the Body Fitting algorithm; + * create prisms of the layers by projecting boundary nodes of offset geometry onto the boundary of initial geometry. + +.. note:: Viscous layers can't be constructed on geometry with **shared/internal** faces. + +To apply the algorithm when you define your mesh, select **Body Fitting** in the list of 3D algorithms and add **Body Fitting Parameters** hypothesis. The following dialog will appear: .. _cartesian_hyp_anchor: diff --git a/resources/StdMeshers.xml.in b/resources/StdMeshers.xml.in index c8f8b5f3d..53899c8f3 100644 --- a/resources/StdMeshers.xml.in +++ b/resources/StdMeshers.xml.in @@ -555,12 +555,14 @@ group-id ="0" priority ="20" hypos ="CartesianParameters3D" + opt-hypos ="ViscousLayers" support-submeshes="false" output ="HEXA" need-hyp ="true" dim ="3"> Cartesian_3D=BodyFitted() + ViscousLayers=ViscousLayers(SetTotalThickness(),SetNumberLayers(),SetStretchFactor(),SetFaces(1),SetFaces(2),SetMethod(),SetGroupName()) diff --git a/src/SMDS/SMDS_LinearEdge.cxx b/src/SMDS/SMDS_LinearEdge.cxx index 38a5260de..42107daf7 100644 --- a/src/SMDS/SMDS_LinearEdge.cxx +++ b/src/SMDS/SMDS_LinearEdge.cxx @@ -45,6 +45,13 @@ SMDS_LinearEdge::SMDS_LinearEdge(const SMDS_MeshNode * node1, myNodes[1] = node2; } +void SMDS_LinearEdge::Init(const SMDS_MeshNode * node1, + const SMDS_MeshNode * node2) +{ + myNodes[0] = node1; + myNodes[1] = node2; +} + int SMDS_LinearEdge::NbNodes() const { return 2; diff --git a/src/SMDS/SMDS_LinearEdge.hxx b/src/SMDS/SMDS_LinearEdge.hxx index 61641a0f6..577f7e955 100644 --- a/src/SMDS/SMDS_LinearEdge.hxx +++ b/src/SMDS/SMDS_LinearEdge.hxx @@ -32,7 +32,8 @@ class SMDS_EXPORT SMDS_LinearEdge: public SMDS_CellOfNodes { public: - SMDS_LinearEdge(const SMDS_MeshNode * node1, const SMDS_MeshNode * node2); + SMDS_LinearEdge(const SMDS_MeshNode * node1=nullptr, const SMDS_MeshNode * node2=nullptr); + void Init( const SMDS_MeshNode * node1, const SMDS_MeshNode * node2); virtual SMDSAbs_ElementType GetType() const { return SMDSAbs_Edge; } virtual SMDSAbs_GeometryType GetGeomType() const { return SMDSGeom_EDGE; } diff --git a/src/SMESH/SMESH_MeshEditor.cxx b/src/SMESH/SMESH_MeshEditor.cxx index be874359b..52addcd4f 100644 --- a/src/SMESH/SMESH_MeshEditor.cxx +++ b/src/SMESH/SMESH_MeshEditor.cxx @@ -7281,6 +7281,7 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes, { const SMDS_MeshElement* elem = *eIt; SMESHDS_SubMesh* sm = mesh->MeshElements( elem->getshapeId() ); + bool marked = elem->isMarked(); bool keepElem = applyMerge( elem, newElemDefs, nodeNodeMap, /*noHoles=*/false ); if ( !keepElem ) @@ -7317,6 +7318,8 @@ void SMESH_MeshEditor::MergeNodes (TListOfListOfNodes & theGroupsOfNodes, sm->AddElement( newElem ); if ( elem != newElem ) ReplaceElemInGroups( elem, newElem, mesh ); + if ( marked && newElem ) + newElem->setIsMarked( true ); } } } diff --git a/src/SMESHDS/SMESHDS_Mesh.cxx b/src/SMESHDS/SMESHDS_Mesh.cxx index 563e2de98..18a7d570a 100644 --- a/src/SMESHDS/SMESHDS_Mesh.cxx +++ b/src/SMESHDS/SMESHDS_Mesh.cxx @@ -1161,6 +1161,18 @@ void SMESHDS_Mesh::UnSetNodeOnShape(const SMDS_MeshNode* aNode) sm->RemoveNode(aNode); } +//======================================================================= +//function : UnSetElementOnShape +//purpose : +//======================================================================= +void SMESHDS_Mesh::UnSetElementOnShape(const SMDS_MeshElement * anElement) +{ + int shapeId = anElement->getshapeId(); + if (shapeId > 0) + if ( SMESHDS_SubMesh* sm = MeshElements( shapeId )) + sm->RemoveElement(anElement); +} + //======================================================================= //function : SetMeshElementOnShape //purpose : diff --git a/src/SMESHDS/SMESHDS_Mesh.hxx b/src/SMESHDS/SMESHDS_Mesh.hxx index dcbd43604..f9eb25b69 100644 --- a/src/SMESHDS/SMESHDS_Mesh.hxx +++ b/src/SMESHDS/SMESHDS_Mesh.hxx @@ -609,6 +609,7 @@ class SMESHDS_EXPORT SMESHDS_Mesh : public SMDS_Mesh void SetNodeOnEdge (const SMDS_MeshNode * aNode, const TopoDS_Edge& S, double u=0.); void SetNodeOnVertex(const SMDS_MeshNode * aNode, const TopoDS_Vertex & S); void UnSetNodeOnShape(const SMDS_MeshNode * aNode); + void UnSetElementOnShape(const SMDS_MeshElement * anElt); void SetMeshElementOnShape (const SMDS_MeshElement * anElt, const TopoDS_Shape & S); void UnSetMeshElementOnShape(const SMDS_MeshElement * anElt, const TopoDS_Shape & S); void SetNodeInVolume(const SMDS_MeshNode * aNode, int Index); diff --git a/src/SMESHUtils/SMESH_MeshAlgos.cxx b/src/SMESHUtils/SMESH_MeshAlgos.cxx index b12d6065d..5dca55c66 100644 --- a/src/SMESHUtils/SMESH_MeshAlgos.cxx +++ b/src/SMESHUtils/SMESH_MeshAlgos.cxx @@ -253,6 +253,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() void init(const SMDS_MeshElement* elem, double tolerance); }; std::vector< ElementBox* > _elements; + //std::string _name; typedef ObjectPool< ElementBox > TElementBoxPool; @@ -290,6 +291,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() if ( theElemIt && !theElemIt->more() ) std::cout << "WARNING: ElementBndBoxTree constructed on empty iterator!" << std::endl; #endif + //_name = "0"; SMDS_ElemIteratorPtr elemIt = theElemIt ? theElemIt : mesh.elementsIterator( elemType ); while ( elemIt->more() ) @@ -342,6 +344,15 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() if ( child->isLeaf() && child->_elements.capacity() > child->_elements.size() ) SMESHUtils::CompactVector( child->_elements ); + + // child->_name = _name + char('0' + j); + // if ( child->isLeaf() && child->_elements.size() ) + // { + // cout << child->_name << " "; + // for ( size_t i = 0; i < child->_elements.size(); ++i ) + // cout << child->_elements[i]->_element->GetID() << " "; + // cout << endl; + // } } } @@ -2566,3 +2577,111 @@ SMESH_ElementSearcher* SMESH_MeshAlgos::GetElementSearcher(SMDS_Mesh& { return new SMESH_ElementSearcherImpl( mesh, tolerance, elemIt ); } + + +//================================================================================ +/*! + * \brief Intersect a ray with a convex volume + * \param [in] ray - the ray + * \param [in] rayLen - ray length + * \param [in] vol - the volume + * \param [out] tMin - return a ray parameter where the ray enters the volume + * \param [out] tMax - return a ray parameter where the ray exit the volume + * \param [out] iFacetMin - facet index where the ray enters the volume + * \param [out] iFacetMax - facet index where the ray exit the volume + * \return bool - true if the ray intersects the volume + */ +//================================================================================ + +bool SMESH_MeshAlgos::IntersectRayVolume( const gp_Ax1& ray, + const double rayLen, + const SMDS_MeshElement* vol, + double & tMin, + double & tMax, + int & iFacetMin, + int & iFacetMax) +{ + /* Ray-Convex Polyhedron Intersection Test by Eric Haines, erich@eye.com + * + * This test checks the ray against each face of a polyhedron, checking whether + * the set of intersection points found for each ray-plane intersection + * overlaps the previous intersection results. If there is no overlap (i.e. + * no line segment along the ray that is inside the polyhedron), then the + * ray misses and returns false; else true. + */ + SMDS_VolumeTool vTool; + if ( !vTool.Set( vol )) + return false; + + tMin = -Precision::Infinite() ; + tMax = rayLen ; + + /* Test each plane in polyhedron */ + for ( int iF = 0; iF < vTool.NbFaces(); ++iF ) + { + const SMDS_MeshNode** fNodes = vTool.GetFaceNodes( iF ); + gp_XYZ normal; + vTool.GetFaceNormal( iF, + normal.ChangeCoord(1), + normal.ChangeCoord(2), + normal.ChangeCoord(3) ); + double D = - ( normal * SMESH_NodeXYZ( fNodes[0] )); + + /* Compute intersection point T and sidedness */ + double vd = ray.Direction().XYZ() * normal; + double vn = ray.Location().XYZ() * normal + D; + if ( vd == 0.0 ) { + /* ray is parallel to plane - check if ray origin is inside plane's + half-space */ + if ( vn > 0.0 ) + /* ray origin is outside half-space */ + return false; + } + else + { + /* ray not parallel - get distance to plane */ + double t = -vn / vd ; + if ( vd < 0.0 ) + { + /* front face - T is a near point */ + if ( t > tMax ) return false; + if ( t > tMin ) { + /* hit near face */ + tMin = t ; + iFacetMin = iF; + } + } + else + { + /* back face - T is a far point */ + if ( t < tMin ) return false; + if ( t < tMax ) { + /* hit far face */ + tMax = t ; + iFacetMax = iF; + } + } + } + } + + /* survived all tests */ + /* Note: if ray originates on polyhedron, may want to change 0.0 to some + * epsilon to avoid intersecting the originating face. + */ + if ( tMin >= 0.0 ) { + /* outside, hitting front face */ + return true; + } + else + { + if ( tMax < rayLen ) { + /* inside, hitting back face */ + return true; + } + else + { + /* inside, but back face beyond tmax */ + return false; + } + } +} diff --git a/src/SMESHUtils/SMESH_MeshAlgos.hxx b/src/SMESHUtils/SMESH_MeshAlgos.hxx index ea5a1ee09..f0c7c2d80 100644 --- a/src/SMESHUtils/SMESH_MeshAlgos.hxx +++ b/src/SMESHUtils/SMESH_MeshAlgos.hxx @@ -168,6 +168,13 @@ namespace SMESH_MeshAlgos const gp_XY& t0, const gp_XY& t1, const gp_XY& t2, double & bc0, double & bc1); + /*! + * \brief Intersect volume by a ray + */ + bool IntersectRayVolume( const gp_Ax1& ray, const double rayLen, + const SMDS_MeshElement* vol, + double & tMin, double & tMax, + int & iFacetMin, int & iFacetMax); /*! * Return a face having linked nodes n1 and n2 and which is * - not in avoidSet, diff --git a/src/SMESHUtils/SMESH_TypeDefs.hxx b/src/SMESHUtils/SMESH_TypeDefs.hxx index 8da7c6a9b..964495bca 100644 --- a/src/SMESHUtils/SMESH_TypeDefs.hxx +++ b/src/SMESHUtils/SMESH_TypeDefs.hxx @@ -202,6 +202,7 @@ struct SMESH_TNodeXYZ : public gp_XYZ } return false; } + void SetXYZ( const gp_XYZ& p ) { SetCoord( p.X(), p.Y(), p.Z() ); } const SMDS_MeshNode* Node() const { return _node; } double Distance(const SMDS_MeshNode* n) const { return (SMESH_TNodeXYZ( n )-*this).Modulus(); } double SquareDistance(const SMDS_MeshNode* n) const { return (SMESH_TNodeXYZ( n )-*this).SquareModulus(); } diff --git a/src/StdMeshers/CMakeLists.txt b/src/StdMeshers/CMakeLists.txt index f3f4c7c27..907727b5b 100644 --- a/src/StdMeshers/CMakeLists.txt +++ b/src/StdMeshers/CMakeLists.txt @@ -120,6 +120,7 @@ SET(StdMeshers_HEADERS StdMeshers_Projection_1D2D.hxx StdMeshers_CartesianParameters3D.hxx StdMeshers_Cartesian_3D.hxx + StdMeshers_Cartesian_VL.hxx StdMeshers_QuadFromMedialAxis_1D2D.hxx StdMeshers_PolygonPerFace_2D.hxx StdMeshers_PolyhedronPerSolid_3D.hxx @@ -183,6 +184,7 @@ SET(StdMeshers_SOURCES StdMeshers_Projection_1D2D.cxx StdMeshers_CartesianParameters3D.cxx StdMeshers_Cartesian_3D.cxx + StdMeshers_Cartesian_VL.cxx StdMeshers_Adaptive1D.cxx StdMeshers_QuadFromMedialAxis_1D2D.cxx StdMeshers_PolygonPerFace_2D.cxx diff --git a/src/StdMeshers/StdMeshers_Cartesian_3D.cxx b/src/StdMeshers/StdMeshers_Cartesian_3D.cxx index 6c50428ff..0d9bfb3d8 100644 --- a/src/StdMeshers/StdMeshers_Cartesian_3D.cxx +++ b/src/StdMeshers/StdMeshers_Cartesian_3D.cxx @@ -24,21 +24,25 @@ // #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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include #include #include @@ -130,7 +134,8 @@ StdMeshers_Cartesian_3D::StdMeshers_Cartesian_3D(int hypId, SMESH_Gen * gen) { _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 @@ -149,19 +154,26 @@ bool StdMeshers_Cartesian_3D::CheckHypothesis (SMESH_Mesh& aMesh, { aStatus = SMESH_Hypothesis::HYP_MISSING; - const list& hyps = GetUsedHypothesis(aMesh, aShape); + const list& hyps = GetUsedHypothesis(aMesh, aShape, /*skipAux=*/false); list ::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( *h ))) + if ( !_hyp && ( _hyp = dynamic_cast( *h ))) { aStatus = _hyp->IsDefined() ? HYP_OK : HYP_BAD_PARAMETER; - break; + } + else + { + _hypViscousLayers = dynamic_cast( *h ); } } @@ -170,7 +182,9 @@ bool StdMeshers_Cartesian_3D::CheckHypothesis (SMESH_Mesh& aMesh, namespace { - typedef int TGeomID; // IDs of sub-shapes + 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; @@ -408,6 +422,11 @@ namespace { _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] ) @@ -419,8 +438,10 @@ namespace 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 @@ -460,11 +481,16 @@ namespace { 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 ) { @@ -1358,6 +1384,20 @@ namespace 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 @@ -1959,11 +1999,11 @@ namespace { 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 { @@ -1978,11 +2018,13 @@ namespace 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; } //============================================================================= @@ -2639,7 +2681,7 @@ namespace 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() ); @@ -6356,6 +6398,27 @@ namespace 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 ); + + std::string error; + TopoDS_Shape offsetShape = builder.MakeOffsetShape( theShape, error ); + if ( offsetShape.IsNull() ) + throw SALOME_Exception( error ); + + SMESH_Mesh* offsetMesh = builder.MakeOffsetMesh(); + + this->_isComputeOffset = true; + if ( ! this->Compute( *offsetMesh, offsetShape )) + return false; + + return builder.MakeViscousLayers( theMesh, theShape ); + } + // The algorithm generates the mesh in following steps: // 1) Intersection of grid lines with the geometry boundary. @@ -6383,6 +6446,11 @@ bool StdMeshers_Cartesian_3D::Compute(SMESH_Mesh & theMesh, grid._toConsiderInternalFaces = _hyp->GetToConsiderInternalFaces(); grid._toUseThresholdForInternalFaces = _hyp->GetToUseThresholdForInternalFaces(); grid._sizeThreshold = _hyp->GetSizeThreshold(); + if ( _isComputeOffset ) + { + grid._toAddEdges = true; + grid._toCreateFaces = true; + } grid.InitGeometry( theShape ); vector< TopoDS_Shape > faceVec; diff --git a/src/StdMeshers/StdMeshers_Cartesian_3D.hxx b/src/StdMeshers/StdMeshers_Cartesian_3D.hxx index 91605135d..db71da85c 100644 --- a/src/StdMeshers/StdMeshers_Cartesian_3D.hxx +++ b/src/StdMeshers/StdMeshers_Cartesian_3D.hxx @@ -37,8 +37,10 @@ * Issue 0021336 * Issue #16523: Treatment of internal faces * Issue #17237: Body fitting on sub-mesh + * Issue #17015: Body fitting with Viscous Layers */ class StdMeshers_CartesianParameters3D; +class StdMeshers_ViscousLayers; class STDMESHERS_EXPORT StdMeshers_Cartesian_3D : public SMESH_3D_Algo { @@ -61,6 +63,8 @@ public: void setSubmeshesComputed(SMESH_Mesh& aMesh, const TopoDS_Shape& theShape ); const StdMeshers_CartesianParameters3D* _hyp; + const StdMeshers_ViscousLayers* _hypViscousLayers; + bool _isComputeOffset; }; #endif diff --git a/src/StdMeshers/StdMeshers_Cartesian_VL.cxx b/src/StdMeshers/StdMeshers_Cartesian_VL.cxx new file mode 100644 index 000000000..b47880198 --- /dev/null +++ b/src/StdMeshers/StdMeshers_Cartesian_VL.cxx @@ -0,0 +1,880 @@ +// Copyright (C) 2007-2021 CEA/DEN, EDF R&D, OPEN CASCADE +// +// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, +// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// File : StdMeshers_Cartesian_VL.cxx +// Created : Tue May 24 13:03:09 2022 +// Author : Edward AGAPOV (eap) + +#include "StdMeshers_Cartesian_VL.hxx" +#include "StdMeshers_ViscousLayers.hxx" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace StdMeshers_Cartesian_VL; + +namespace +{ + typedef int TGeomID; // IDs of sub-shapes + + /*! + * \brief Temporary mesh + */ + struct TmpMesh: public SMESH_Mesh + { + TmpMesh() { + _isShapeToMesh = (_id = 0); + _meshDS = new SMESHDS_Mesh( _id, true ); + } + }; + // -------------------------------------------------------------------------- + /*! + * \brief Edge of viscous layers; goes from a grid node by normal to boundary + */ + struct VLEdge + { + std::vector< SMESH_NodeXYZ > _nodes; + double _uv[2]; // position of TgtNode() on geometry + double _length; + + const SMDS_MeshNode* TgtNode() const { return _nodes.back().Node(); } + }; + typedef NCollection_DataMap< const SMDS_MeshNode*, VLEdge* > TNode2VLEdge; + // -------------------------------------------------------------------------- + /*! + * \brief VLEdge's of one shape + */ + struct VLEdgesOnShape + { + std::vector< VLEdge > _edges; + TopoDS_Shape _initShape; + TopoDS_Shape _offsetShape; + int _initShapeID; + bool _toCheckCoinc; // to check if nodes on shape boundary coincide + }; + + //================================================================================ + /*! + * \brief Project a point on offset FACE to EDGEs of an initial FACE + * \param [in] offP - point to project + * \param [in] initFace - FACE to project on + * \return gp_Pnt - projection point + */ + //================================================================================ + + gp_Pnt projectToWire( const gp_Pnt& offP, + const TopoDS_Face& initFace, + gp_Pnt2d & uv ) + { + double minDist = Precision::Infinite(); + const double tol = Precision::Confusion(); + gp_Pnt proj, projction; + Standard_Real param; + ShapeAnalysis_Curve projector; + + for ( TopExp_Explorer eExp( initFace, TopAbs_EDGE ); eExp.More(); eExp.Next() ) + { + BRepAdaptor_Curve initCurve( TopoDS::Edge( eExp.Current() )); + //const double f = initCurve.FirstParameter(), l = initCurve.LastParameter(); + double dist = projector.Project( initCurve, offP, tol, proj, param, /*adjustToEnds=*/false ); + if ( dist < minDist ) + { + projction = proj; + minDist = dist; + } + } + uv.SetCoord(0.,0.); // !!!!!!! + return projction; + } + + //================================================================================ + /*! + * \brief Project nodes from offset FACE to initial FACE + * \param [inout] theEOS - VL edges on a geom FACE + * \param [inout] theOffsetMDS - offset mesh to fill in + */ + //================================================================================ + + void projectToFace( VLEdgesOnShape & theEOS, + SMESHDS_Mesh* theOffsetMDS, + TNode2VLEdge & theN2E ) + { + SMESHDS_SubMesh* sm = theOffsetMDS->MeshElements( theEOS._offsetShape ); + if ( !sm || sm->NbElements() == 0 || sm->NbNodes() == 0 ) + return; + theEOS._edges.resize( sm->NbNodes() ); + + const TopoDS_Face& initFace = TopoDS::Face( theEOS._initShape ); + ShapeAnalysis_Surface projector( BRep_Tool::Surface( initFace )); + const double clsfTol = 1e2 * BRep_Tool::MaxTolerance( initFace, TopAbs_VERTEX ); + BRepTopAdaptor_FClass2d classifier( initFace, clsfTol ); + + const double tol = Precision::Confusion(); + //const double f = initCurve.FirstParameter(), l = initCurve.LastParameter(); + gp_Pnt proj; + + int iN = 0; + for ( SMDS_NodeIteratorPtr nIt = sm->GetNodes(); nIt->more(); ++iN ) + { + SMESH_NodeXYZ offP = nIt->next(); + gp_Pnt2d uv = projector.ValueOfUV( offP, tol ); + TopAbs_State st = classifier.Perform( uv ); + if ( st == TopAbs_IN || st == TopAbs_ON ) + { + proj = projector.Value( uv ); + } + else + { + proj = projectToWire( offP, initFace, uv ); + } + if ( st == TopAbs_ON || st == TopAbs_OUT ) + theEOS._toCheckCoinc = true; + + VLEdge & vlEdge = theEOS._edges[ iN ]; + vlEdge._nodes.push_back( offP.Node() ); + vlEdge._nodes.push_back( theOffsetMDS->AddNode( proj.X(), proj.Y(), proj.Z() )); + vlEdge._uv[0] = uv.X(); + vlEdge._uv[1] = uv.Y(); + //vlEdge._length = proj.Distance( offP ); + + theN2E.Bind( offP.Node(), &vlEdge ); + } + return; + + } + + //================================================================================ + /*! + * \brief Project nodes from offset EDGE to initial EDGE + * \param [inout] theEOS - VL edges on a geom EDGE + * \param [inout] theOffsetMDS - offset mesh to add new nodes to + */ + //================================================================================ + + void projectToEdge( VLEdgesOnShape & theEOS, + SMESHDS_Mesh* theOffsetMDS, + TNode2VLEdge & theN2E ) + { + SMESHDS_SubMesh* sm = theOffsetMDS->MeshElements( theEOS._offsetShape ); + if ( !sm || sm->NbElements() == 0 ) + return; + theEOS._edges.resize( sm->NbNodes() ); + + ShapeAnalysis_Curve projector; + BRepAdaptor_Curve initCurve( TopoDS::Edge( theEOS._initShape )); + const double tol = Precision::Confusion(); + const double f = initCurve.FirstParameter(), l = initCurve.LastParameter(); + gp_Pnt proj; + Standard_Real param; + + int iN = 0; + for ( SMDS_NodeIteratorPtr nIt = sm->GetNodes(); nIt->more(); ++iN ) + { + SMESH_NodeXYZ offP = nIt->next(); + double dist = projector.Project( initCurve, offP, tol, proj, param, /*adjustToEnds=*/false ); + bool paramOK = ( f < param && param < l ); + if ( !paramOK ) + { + if ( param < f ) + param = f; + else if ( param > l ) + param = l; + theEOS._toCheckCoinc = true; + } + proj = initCurve.Value( param ); + + VLEdge & vlEdge = theEOS._edges[ iN ]; + vlEdge._nodes.push_back( offP.Node() ); + vlEdge._nodes.push_back( theOffsetMDS->AddNode( proj.X(), proj.Y(), proj.Z() )); + vlEdge._uv[0] = param; + vlEdge._length = paramOK ? dist : proj.Distance( offP ); + + theN2E.Bind( offP.Node(), &vlEdge ); + } + return; + } + + //================================================================================ + /*! + * \brief Compute heights of viscous layers and finds FACEs with VL + * \param [in] hyp - viscous layers hypothesis + * \param [in] mesh - the main mesh + * \param [in] shape - the main shape + * \param [out] vlH - heights of viscous layers to compute + * \param [out] shapesWVL - IDs of shapes with VL + * \return bool - isOK + */ + //================================================================================ + + bool computeVLHeight( const StdMeshers_ViscousLayers* hyp, + SMESHDS_Mesh * mesh, + const TopoDS_Shape & shape, + std::vector< double > & vlH, + std::set< TGeomID > & shapesWVL ) + { + const double T = hyp->GetTotalThickness(); + const double f = hyp->GetStretchFactor(); + const int N = hyp->GetNumberLayers(); + const double h0 = hyp->Get1stLayerThickness( T, f, N ); + + vlH.reserve( hyp->GetNumberLayers() ); + vlH.push_back( h0 ); + for ( double h = h0 * f; (int)vlH.size() < N-1; h *= f ) + vlH.push_back( h + vlH.back() ); + vlH.push_back( T ); + + + TopTools_IndexedMapOfShape faces; + TopExp::MapShapes( shape, TopAbs_FACE, faces ); + + if ( hyp->IsToIgnoreShapes() ) + { + for ( int i = 1; i <= faces.Size(); ++i ) + shapesWVL.insert( mesh->ShapeToIndex( faces( i ))); + + for ( TGeomID& sID : hyp->GetBndShapes() ) + shapesWVL.erase( sID ); + } + else + { + for ( TGeomID& sID : hyp->GetBndShapes() ) + shapesWVL.insert( sID ); + } + + for ( TGeomID fID : shapesWVL ) + { + const TopoDS_Shape & face = mesh->IndexToShape( fID ); + for ( TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next() ) + shapesWVL.insert( mesh->ShapeToIndex( exp.Current() )); + for ( TopExp_Explorer exp( face, TopAbs_VERTEX ); exp.More(); exp.Next() ) + shapesWVL.insert( mesh->ShapeToIndex( exp.Current() )); + } + + return !shapesWVL.empty(); + } + + //================================================================================ + /*! + * \brief Create intermediate nodes on VLEdge + * \param [inout] vlEdge - VLEdge to divide + * \param [in] vlH - thicknesses of layers + * \param [inout] mesh - the mesh no fill in + */ + //================================================================================ + + void divideVLEdge( VLEdge* vlEdge, + const std::vector< double >& vlH, + SMESHDS_Mesh* mesh ) + { + SMESH_NodeXYZ lastNode = vlEdge->_nodes.back(); + SMESH_NodeXYZ frstNode = vlEdge->_nodes[0]; + gp_XYZ dir = frstNode - lastNode; + + vlEdge->_nodes.resize( vlH.size() + 1 ); + + for ( size_t i = 1; i < vlH.size(); ++i ) + { + double r = vlH[ i-1 ] / vlH.back(); + gp_XYZ p = lastNode + dir * r; + vlEdge->_nodes[ vlH.size() - i ] = mesh->AddNode( p.X(), p.Y(), p.Z() ); + } + vlEdge->_nodes.back() = lastNode; + } + + //================================================================================ + /*! + * \brief Create a polyhedron from nodes of VLEdge's + * \param [in] edgesVec - the VLEdge's + * \param [in] vNodes - node buffer + * \param [in] editor - editor of offset mesh + */ + //================================================================================ + + bool makePolyhedron( const std::vector< VLEdge*> & edgesVec, + std::vector< const SMDS_MeshNode* > & vNodes, + SMESH_MeshEditor& editor, + SMESH_MeshEditor::ElemFeatures elemType) + { + elemType.SetPoly( true ); + elemType.myPolyhedQuantities.clear(); + elemType.myNodes.clear(); + + // add facets of walls + size_t nbBaseNodes = edgesVec.size(); + for ( size_t iN = 0; iN < nbBaseNodes; ++iN ) + { + VLEdge* e0 = edgesVec[ iN ]; + VLEdge* e1 = edgesVec[( iN + 1 ) % nbBaseNodes ]; + size_t nbN0 = e0->_nodes.size(); + size_t nbN1 = e1->_nodes.size(); + if ( nbN0 == nbN1 ) + { + for ( size_t i = 1; i < nbN0; ++i ) + { + elemType.myNodes.push_back( e1->_nodes[ i - 1 ].Node()); + elemType.myNodes.push_back( e1->_nodes[ i ].Node()); + elemType.myNodes.push_back( e0->_nodes[ i ].Node()); + elemType.myNodes.push_back( e0->_nodes[ i - 1 ].Node()); + elemType.myPolyhedQuantities.push_back( 4 ); + } + } + else + { + for ( size_t i = 0; i < nbN1; ++i ) + elemType.myNodes.push_back( e1->_nodes[ i ].Node() ); + for ( size_t i = 0; i < nbN0; ++i ) + elemType.myNodes.push_back( e0->_nodes[ nbN0 - 1 - i ].Node() ); + elemType.myPolyhedQuantities.push_back( nbN0 + nbN1 ); + } + } + + // add facets of top + vNodes.clear(); + for ( size_t iN = 0; iN < nbBaseNodes; ++iN ) + { + VLEdge* e0 = edgesVec[ iN ]; + elemType.myNodes.push_back( e0->_nodes.back().Node() ); + vNodes.push_back( e0->_nodes[ 0 ].Node()); + } + elemType.myPolyhedQuantities.push_back( nbBaseNodes ); + + // add facets of bottom + elemType.myNodes.insert( elemType.myNodes.end(), vNodes.rbegin(), vNodes.rend() ); + elemType.myPolyhedQuantities.push_back( nbBaseNodes ); + + const SMDS_MeshElement* vol = editor.AddElement( elemType.myNodes, elemType ); + vol->setIsMarked( true ); // to add to group + + return vol; + } + + //================================================================================ + /*! + * \brief Transform prism nodal connectivity to that of polyhedron + * \param [inout] nodes - nodes of prism, return nodes of polyhedron + * \param [inout] elemType - return quantities of polyhedron, + */ + //================================================================================ + + void prism2polyhedron( std::vector< const SMDS_MeshNode* > & nodes, + SMESH_MeshEditor::ElemFeatures & elemType ) + { + elemType.myPolyhedQuantities.clear(); + elemType.myNodes.clear(); + + // walls + size_t nbBaseNodes = nodes.size() / 2; + for ( size_t i = 0; i < nbBaseNodes; ++i ) + { + int iNext = ( i + 1 ) % nbBaseNodes; + elemType.myNodes.push_back( nodes[ iNext ]); + elemType.myNodes.push_back( nodes[ i ]); + elemType.myNodes.push_back( nodes[ i + nbBaseNodes ]); + elemType.myNodes.push_back( nodes[ iNext + nbBaseNodes ]); + elemType.myPolyhedQuantities.push_back( 4 ); + } + + // base + elemType.myNodes.insert( elemType.myNodes.end(), nodes.begin(), nodes.begin() + nbBaseNodes ); + elemType.myPolyhedQuantities.push_back( nbBaseNodes ); + + // top + elemType.myNodes.insert( elemType.myNodes.end(), nodes.rbegin(), nodes.rbegin() + nbBaseNodes ); + elemType.myPolyhedQuantities.push_back( nbBaseNodes ); + + nodes.swap( elemType.myNodes ); + } + + //================================================================================ + /*! + * \brief Create prisms from faces + * \param [in] theEOS - shape to treat + * \param [inout] theMesh - offset mesh to fill in + */ + //================================================================================ + + bool makePrisms( VLEdgesOnShape & theEOS, + SMESH_Mesh* theMesh, + TNode2VLEdge & theN2E ) + { + SMESHDS_SubMesh* sm = theMesh->GetMeshDS()->MeshElements( theEOS._offsetShape ); + if ( !sm || sm->NbElements() == 0 ) + return true; + + SMESH_MeshEditor editor( theMesh ); + SMESH_MeshEditor::ElemFeatures volumElem( SMDSAbs_Volume ); + + std::vector< const SMDS_MeshNode* > vNodes; + std::vector< VLEdge*> edgesVec; + for ( SMDS_ElemIteratorPtr eIt = sm->GetElements(); eIt->more(); ) + { + const SMDS_MeshElement* face = eIt->next(); + if ( face->GetType() != SMDSAbs_Face ) + continue; + + const int nbNodes = face->NbCornerNodes(); + edgesVec.resize( nbNodes ); + vNodes.resize( nbNodes * 2 ); + int maxNbLayer = 0, minNbLayer = IntegerLast(); + for ( int i = 0; i < nbNodes; ++i ) + { + const SMDS_MeshNode* n = face->GetNode( i ); + if ( !theN2E.IsBound( n )) + return false; + edgesVec[ i ] = theN2E( n ); + maxNbLayer = Max( maxNbLayer, edgesVec[i]->_nodes.size() - 1 ); + minNbLayer = Min( minNbLayer, edgesVec[i]->_nodes.size() - 1 ); + } + if ( maxNbLayer == minNbLayer ) + { + size_t nbPrism = edgesVec[0]->_nodes.size() - 1; + for ( size_t iP = 0; iP < nbPrism; ++iP ) + { + vNodes.resize( nbNodes * 2 ); + for ( int i = 0; i < nbNodes; ++i ) + { + vNodes[ i ] = edgesVec[ i ]->_nodes[ iP + 1 ].Node(); + vNodes[ i + nbNodes ] = edgesVec[ i ]->_nodes[ iP ].Node(); + } + volumElem.SetPoly( nbNodes > 4 ); + if ( volumElem.myIsPoly ) + prism2polyhedron( vNodes, volumElem ); + + if ( const SMDS_MeshElement* vol = editor.AddElement( vNodes, volumElem )) + vol->setIsMarked( true ); // to add to group + } + } + else // at inlet/outlet + { + makePolyhedron( edgesVec, vNodes, editor, volumElem ); + } + editor.ClearLastCreated(); + + // move the face to the top of prisms, on mesh boundary + //theMesh->GetMeshDS()->ChangeElementNodes( face, fNodes.data(), nbNodes ); + } + return true; + } + + //================================================================================ + /*! + * \brief Append the offset mesh to the initial one + */ + //================================================================================ + + void copyMesh( SMESH_Mesh* theFromMesh, + SMESH_Mesh* theToMesh, + int theSolidID) + { + SMESHDS_Mesh* offsetMDS = theFromMesh->GetMeshDS(); + SMESHDS_Mesh* initMDS = theToMesh->GetMeshDS(); + + const smIdType nShift = initMDS->NbNodes(); + + for ( SMDS_NodeIteratorPtr nIt = offsetMDS->nodesIterator(); nIt->more(); ) + { + SMESH_NodeXYZ pOff = nIt->next(); + const SMDS_MeshNode* n = initMDS->AddNodeWithID( pOff.X(), pOff.Y(), pOff.Z(), + nShift + pOff.Node()->GetID() ); + initMDS->SetNodeInVolume( n, theSolidID ); + } + + SMESH_MeshEditor editor( theToMesh ); + SMESH_MeshEditor::ElemFeatures elemType; + std::vector nIniIDs; + + for ( SMDS_ElemIteratorPtr eIt = offsetMDS->elementsIterator(); eIt->more(); ) + { + const SMDS_MeshElement* offElem = eIt->next(); + const int nbNodes = offElem->NbNodes(); + nIniIDs.resize( nbNodes ); + for ( int i = 0; i < nbNodes; ++i ) + { + const SMDS_MeshNode* offNode = offElem->GetNode( i ); + nIniIDs[ i ] = offNode->GetID() + nShift; + } + elemType.Init( offElem, /*basicOnly=*/false ); + const SMDS_MeshElement* iniElem = editor.AddElement( nIniIDs, elemType ); + initMDS->SetMeshElementOnShape( iniElem, theSolidID ); + iniElem->setIsMarked( offElem->isMarked() ); + editor.ClearLastCreated(); + } + return; + } + + //================================================================================ + /*! + * \brief set elements of layers boundary to sub-meshes + * \param [in] theEOS - vlEdges of a sub-shape + * \param [inout] theMesh - the mesh + * \param [in] theN2E - map of node to vlEdge + * \param [in] theNShift - nb of nodes in the mesh before adding the offset mesh + */ + //================================================================================ + + void setBnd2Sub( VLEdgesOnShape & theEOS, + SMESH_Mesh* theInitMesh, + SMESH_Mesh* theOffsMesh, + const TNode2VLEdge & theN2E, + const smIdType theNShift, + TIDSortedNodeSet& theNodesToCheckCoinc ) + { + SMESHDS_Mesh* offsetMDS = theOffsMesh->GetMeshDS(); + SMESHDS_Mesh* initMDS = theInitMesh->GetMeshDS(); + + SMESHDS_SubMesh* sm = offsetMDS->MeshElements( theEOS._offsetShape ); + if ( !sm || ( sm->NbElements() + sm->NbNodes() == 0 )) + return; + + for ( SMDS_NodeIteratorPtr nIt = sm->GetNodes(); nIt->more(); ) + { + const SMDS_MeshNode* offNode = nIt->next(); + VLEdge* vlEd = theN2E( offNode ); + const SMDS_MeshNode* nOffBnd = vlEd->_nodes.back().Node(); + const SMDS_MeshNode* nIniBnd = initMDS->FindNode( nOffBnd->GetID() + theNShift ); + initMDS->UnSetNodeOnShape( nIniBnd ); + + switch( theEOS._initShape.ShapeType() ) { + case TopAbs_FACE: initMDS->SetNodeOnFace( nIniBnd, TopoDS::Face( theEOS._initShape ), + vlEd->_uv[0], vlEd->_uv[1] ); + break; + case TopAbs_EDGE: initMDS->SetNodeOnEdge( nIniBnd, TopoDS::Edge( theEOS._initShape ), + vlEd->_uv[0]); + break; + case TopAbs_VERTEX: initMDS->SetNodeOnVertex( nIniBnd, TopoDS::Vertex( theEOS._initShape )); + break; + default:; + } + } + + std::vector< const SMDS_MeshNode* > iniNodes, iniNodesBnd; + std::vector< VLEdge*> edgesVec; + for ( SMDS_ElemIteratorPtr eIt = sm->GetElements(); eIt->more(); ) + { + const SMDS_MeshElement* offElem = eIt->next(); + if ( offElem->GetType() != SMDSAbs_Face && + offElem->GetType() != SMDSAbs_Edge ) + continue; + + const int nbNodes = offElem->NbCornerNodes(); + iniNodes.resize( nbNodes ); + iniNodesBnd.resize( nbNodes ); + //edgesVec.resize( nbNodes ); + for ( int i = 0; i < nbNodes; ++i ) + { + const SMDS_MeshNode* nOff = offElem->GetNode( i ); + const SMDS_MeshNode* nIni = initMDS->FindNode( nOff->GetID() + theNShift ); + iniNodes[ i ] = nIni; + VLEdge* vlEd = theN2E( nOff ); + const SMDS_MeshNode* nOffBnd = vlEd->_nodes.back().Node(); + const SMDS_MeshNode* nIniBnd = initMDS->FindNode( nOffBnd->GetID() + theNShift ); + iniNodesBnd[ i ] = nIniBnd; + } + if ( const SMDS_MeshElement* iniElem = initMDS->FindElement( iniNodes )) + { + initMDS->UnSetElementOnShape( iniElem ); + initMDS->SetMeshElementOnShape( iniElem, theEOS._initShape ); + + // move the face to the top of prisms, on init mesh boundary + initMDS->ChangeElementNodes( iniElem, iniNodesBnd.data(), nbNodes ); + + if ( theEOS._toCheckCoinc ) + theNodesToCheckCoinc.insert( iniNodesBnd.begin(), iniNodesBnd.end() ); + } + } + return; + } +} // namespace + + +//================================================================================ +/*! + * \brief Create a temporary mesh used to hold elements of the offset shape + */ +//================================================================================ + +SMESH_Mesh* StdMeshers_Cartesian_VL::ViscousBuilder::MakeOffsetMesh() +{ + return _offsetMesh; +} + +//================================================================================ +/*! + * \brief Constructor + */ +//================================================================================ + +StdMeshers_Cartesian_VL::ViscousBuilder::ViscousBuilder( const StdMeshers_ViscousLayers* hyp ) + : _hyp( hyp ), _offsetMesh( new TmpMesh ) +{ +} + +//================================================================================ +/*! + * \brief Destructor + */ +//================================================================================ + +StdMeshers_Cartesian_VL::ViscousBuilder::~ViscousBuilder() +{ + delete _offsetMesh; _offsetMesh = 0; +} + +//================================================================================ +/*! + * \brief Create an offset shape from a given one + * \param [in] theShape - input shape + * \param [out] theError - error description + * \return TopoDS_Shape - result offset shape + */ +//================================================================================ + +TopoDS_Shape +StdMeshers_Cartesian_VL::ViscousBuilder::MakeOffsetShape(const TopoDS_Shape & theShape, + std::string & theError ) +{ + // TopoDS_Shape S = theShape; + // BRep_Builder aBuilder; + // TopoDS_Compound comp; + // aBuilder.MakeCompound( comp ); + // for ( TopExp_Explorer e( theShape, TopAbs_FACE); e.More(); e.Next() ) + // { + // aBuilder.Add( comp, e.Current()); + // e.Next(); + // if ( !e.More() ) + // break; + // } + // S = comp; + + double offset = -_hyp->GetTotalThickness(); + double tol = Precision::Confusion(); + _makeOffset.Initialize( theShape, offset, tol, BRepOffset_Skin, /*Intersection=*/false, + /*selfInter=*/false, GeomAbs_Intersection ); + _makeOffset.MakeOffsetShape(); + if ( _makeOffset.IsDone() ) + { + _offsetShape = _makeOffset.Shape(); + _offsetMesh->ShapeToMesh( _offsetShape ); + _offsetMesh->GetSubMesh( _offsetShape )->DependsOn(); + //SMESH_MesherHelper::WriteShape( _offsetShape ); + return _offsetShape; + } + + switch ( _makeOffset.Error() ) + { + case BRepOffset_NoError: + theError = "OK. Offset performed successfully.";break; + case BRepOffset_BadNormalsOnGeometry: + theError = "Degenerated normal on input data.";break; + case BRepOffset_C0Geometry: + theError = "C0 continuity of input data.";break; + case BRepOffset_NullOffset: + theError = "Null offset of all faces.";break; + case BRepOffset_NotConnectedShell: + theError = "Incorrect set of faces to remove, the remaining shell is not connected.";break; + case BRepOffset_CannotTrimEdges: + theError = "Can not trim edges.";break; + case BRepOffset_CannotFuseVertices: + theError = "Can not fuse vertices.";break; + case BRepOffset_CannotExtentEdge: + theError = "Can not extent edge.";break; + default: + theError = "operation not done."; + } + theError = "BRepOffset_MakeOffset error: " + theError; + + return TopoDS_Shape(); +} + +//================================================================================ +/*! + * \brief Return a sub-shape of the offset shape generated from a given initial sub-shape + */ +//================================================================================ + +TopoDS_Shape StdMeshers_Cartesian_VL::ViscousBuilder::getOffsetSubShape( const TopoDS_Shape& S ) +{ + const TopTools_ListOfShape& newShapes = _makeOffset.Generated( S ); + for ( const TopoDS_Shape& ns : newShapes ) + if ( ns.ShapeType() == S.ShapeType() ) + return ns; + return TopoDS_Shape(); +} + +//================================================================================ +/*! + * \brief Create prismatic mesh between _offsetShape and theShape + * \param [out] theMesh - mesh to fill in + * \param [in] theShape - initial shape + * \return bool - is Ok + */ +//================================================================================ + +bool StdMeshers_Cartesian_VL::ViscousBuilder::MakeViscousLayers( SMESH_Mesh & theMesh, + const TopoDS_Shape & theShape ) +{ + SMESHDS_Mesh* offsetMDS = _offsetMesh->GetMeshDS(); + SMESHDS_Mesh* initMDS = theMesh.GetMeshDS(); + offsetMDS->SetAllCellsNotMarked(); + + // Compute heights of viscous layers and finds FACEs with VL + std::vector< double > vlH; + std::set< int > shapesWVL; + if ( !computeVLHeight( _hyp, initMDS, theShape, vlH, shapesWVL )) + return false; + + std::vector< VLEdgesOnShape > edgesOnShape; + edgesOnShape.reserve( offsetMDS->MaxShapeIndex() + 1 ); + TNode2VLEdge n2e; + + // loop on sub-shapes to project nodes from offset boundary to initial boundary + TopAbs_ShapeEnum types[3] = { TopAbs_VERTEX, TopAbs_EDGE, TopAbs_FACE }; + for ( TopAbs_ShapeEnum shType : types ) + { + TopTools_IndexedMapOfShape shapes; + TopExp::MapShapes( theShape, shType, shapes ); + for ( int i = 1; i <= shapes.Size(); ++i ) + { + edgesOnShape.resize( edgesOnShape.size() + 1 ); + VLEdgesOnShape& EOS = edgesOnShape.back(); + + EOS._initShape = shapes( i ); + EOS._offsetShape = getOffsetSubShape( EOS._initShape ); + EOS._initShapeID = initMDS->ShapeToIndex( EOS._initShape ); + EOS._toCheckCoinc = false; + + // project boundary nodes of offset mesh to boundary of init mesh + // (new nodes are created in the offset mesh) + switch( EOS._offsetShape.ShapeType() ) { + case TopAbs_VERTEX: + { + EOS._edges.resize( 1 ); + EOS._edges[0]._nodes.resize( 2 ); + EOS._edges[0]._nodes[0] = SMESH_Algo::VertexNode( TopoDS::Vertex( EOS._offsetShape ), + offsetMDS ); + gp_Pnt offP = BRep_Tool::Pnt( TopoDS::Vertex( EOS._initShape )); + EOS._edges[0]._nodes[1] = offsetMDS->AddNode( offP.X(), offP.Y(), offP.Z() ); + //EOS._edges[0]._length = offP.Distance( EOS._edges[0]._nodes[0] ); + n2e.Bind( EOS._edges[0]._nodes[0].Node(), & EOS._edges[0] ); + break; + } + case TopAbs_EDGE: + { + projectToEdge( EOS, offsetMDS, n2e ); + break; + } + case TopAbs_FACE: + { + projectToFace( EOS, offsetMDS, n2e ); + break; + } + default:; + } + + // create nodes of layers + if ( _hyp->GetNumberLayers() > 1 ) + { + if ( shapesWVL.count( EOS._initShapeID )) + for ( size_t i = 0; i < EOS._edges.size(); ++i ) + { + divideVLEdge( &EOS._edges[ i ], vlH, offsetMDS ); + } + } + } // loop on shapes + } // loop on shape types + + // create prisms + bool prismsOk = true; + for ( size_t i = 0; i < edgesOnShape.size(); ++i ) + { + VLEdgesOnShape& EOS = edgesOnShape[ i ]; + if ( EOS._initShape.ShapeType() == TopAbs_FACE ) + { + if ( !makePrisms( EOS, _offsetMesh, n2e )) + prismsOk = false; + } + } + + // copy offset mesh to the main one + initMDS->Modified(); + initMDS->CompactMesh(); + smIdType nShift = initMDS->NbNodes(); + TGeomID solidID = initMDS->ShapeToIndex( theShape ); + copyMesh( _offsetMesh, & theMesh, solidID ); + + + if ( !prismsOk ) + { + if ( SMESH_subMesh * sm = theMesh.GetSubMesh( theShape )) + { + sm->GetComputeError() = + SMESH_ComputeError::New( COMPERR_ALGO_FAILED, + "Viscous layers construction error: bad mesh on offset geometry" ); + } + return prismsOk; + } + + // set elements of layers boundary to sub-meshes + TIDSortedNodeSet nodesToCheckCoinc; + for ( size_t i = 0; i < edgesOnShape.size(); ++i ) + { + VLEdgesOnShape& EOS = edgesOnShape[ i ]; + setBnd2Sub( EOS, &theMesh, _offsetMesh, n2e, nShift, nodesToCheckCoinc ); + } + + // merge coincident nodes + SMESH_MeshEditor editor( &theMesh ); + SMESH_MeshEditor::TListOfListOfNodes nodesToMerge; + editor.FindCoincidentNodes( nodesToCheckCoinc, vlH[0]/20., nodesToMerge, false ); + editor.MergeNodes( nodesToMerge ); + + // create a group + if ( !_hyp->GetGroupName().empty() ) + { + SMDS_MeshGroup* group = _hyp->CreateGroup( _hyp->GetGroupName(), theMesh, SMDSAbs_Volume ); + + for ( SMDS_ElemIteratorPtr it = initMDS->elementsIterator(); it->more(); ) + { + const SMDS_MeshElement* el = it->next(); + if ( el->isMarked() ) + group->Add( el ); + } + } + + return prismsOk; +} diff --git a/src/StdMeshers/StdMeshers_Cartesian_VL.hxx b/src/StdMeshers/StdMeshers_Cartesian_VL.hxx new file mode 100644 index 000000000..ddcfe0bfa --- /dev/null +++ b/src/StdMeshers/StdMeshers_Cartesian_VL.hxx @@ -0,0 +1,62 @@ +// Copyright (C) 2007-2021 CEA/DEN, EDF R&D, OPEN CASCADE +// +// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, +// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// File : StdMeshers_Cartesian_VL.hxx +// Created : Tue May 24 12:32:01 2022 +// Author : Edward AGAPOV (eap) + +#ifndef __StdMeshers_Cartesian_VL_HXX__ +#define __StdMeshers_Cartesian_VL_HXX__ + +#include + +class StdMeshers_ViscousLayers; +class SMESH_Mesh; + +namespace StdMeshers_Cartesian_VL +{ + class ViscousBuilder + { + public: + + ViscousBuilder( const StdMeshers_ViscousLayers* hypViscousLayers ); + ~ViscousBuilder(); + + TopoDS_Shape MakeOffsetShape(const TopoDS_Shape & theShape, + std::string & theError ); + + SMESH_Mesh* MakeOffsetMesh(); + + bool MakeViscousLayers( SMESH_Mesh & theMesh, + const TopoDS_Shape & theShape ); + + private: + + TopoDS_Shape getOffsetSubShape( const TopoDS_Shape& S ); + + const StdMeshers_ViscousLayers* _hyp; + BRepOffset_MakeOffset _makeOffset; + SMESH_Mesh* _offsetMesh; + TopoDS_Shape _offsetShape; + }; +} + +#endif diff --git a/src/StdMeshers/StdMeshers_ViscousLayers.cxx b/src/StdMeshers/StdMeshers_ViscousLayers.cxx index 680708d7e..7b38e23f8 100644 --- a/src/StdMeshers/StdMeshers_ViscousLayers.cxx +++ b/src/StdMeshers/StdMeshers_ViscousLayers.cxx @@ -635,13 +635,7 @@ namespace VISCOUS_3D const double T = ( realThickness > 0 ) ? realThickness : GetTotalThickness(); const double f = GetStretchFactor(); const int N = GetNumberLayers(); - const double fPowN = pow( f, N ); - double h0; - if ( fPowN - 1 <= numeric_limits::min() ) - h0 = T / N; - else - h0 = T * ( f - 1 )/( fPowN - 1 ); - return h0; + return StdMeshers_ViscousLayers::Get1stLayerThickness( T, f, N ); } bool UseSurfaceNormal() const @@ -1445,7 +1439,17 @@ bool StdMeshers_ViscousLayers::IsShapeWithLayers(int shapeIndex) const ( std::find( _shapeIds.begin(), _shapeIds.end(), shapeIndex ) != _shapeIds.end() ); return IsToIgnoreShapes() ? !isIn : isIn; } - +// -------------------------------------------------------------------------------- +double StdMeshers_ViscousLayers::Get1stLayerThickness( double T, double f, int N ) +{ + const double fPowN = pow( f, N ); + double h0; + if ( fPowN - 1 <= numeric_limits::min() ) + h0 = T / N; + else + h0 = T * ( f - 1 )/( fPowN - 1 ); + return h0; +} // -------------------------------------------------------------------------------- SMDS_MeshGroup* StdMeshers_ViscousLayers::CreateGroup( const std::string& theName, SMESH_Mesh& theMesh, diff --git a/src/StdMeshers/StdMeshers_ViscousLayers.hxx b/src/StdMeshers/StdMeshers_ViscousLayers.hxx index 96e6b53c5..ec0b35c43 100644 --- a/src/StdMeshers/StdMeshers_ViscousLayers.hxx +++ b/src/StdMeshers/StdMeshers_ViscousLayers.hxx @@ -94,6 +94,11 @@ public: const TopoDS_Shape& aShape, SMESH_Hypothesis::Hypothesis_Status& aStatus); + // Compute thickness of the 1st layer + static double Get1stLayerThickness( double thickTotal, + double factor, + int nbLayers ); + // Checks if viscous layers should be constructed on a shape bool IsShapeWithLayers(int shapeIndex) const; -- 2.39.2