1 // Copyright (C) 2007-2023 CEA, EDF, OPEN CASCADE
3 // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License, or (at your option) any later version.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
22 // File : StdMeshers_Cartesian_VL.cxx
23 // Created : Tue May 24 13:03:09 2022
24 // Author : Edward AGAPOV (eap)
26 #include "StdMeshers_Cartesian_VL.hxx"
27 #include "StdMeshers_ViscousLayers.hxx"
29 #include <SMDS_MeshGroup.hxx>
30 #include <SMESHDS_Mesh.hxx>
31 #include <SMESHDS_SubMesh.hxx>
32 #include <SMESH_Algo.hxx>
33 #include <SMESH_Mesh.hxx>
34 #include <SMESH_MeshEditor.hxx>
35 #include <SMESH_MesherHelper.hxx>
36 #include <SMESH_TypeDefs.hxx>
37 #include <SMESH_subMesh.hxx>
39 #include <BRepAdaptor_Curve.hxx>
40 #include <BRepTopAdaptor_FClass2d.hxx>
41 #include <BRep_Builder.hxx>
42 #include <BRep_Tool.hxx>
43 #include <ShapeAnalysis_Curve.hxx>
44 #include <ShapeAnalysis_Surface.hxx>
46 #include <TopExp_Explorer.hxx>
49 using namespace StdMeshers_Cartesian_VL;
53 typedef int TGeomID; // IDs of sub-shapes
56 * \brief Temporary mesh
58 struct TmpMesh: public SMESH_Mesh
61 _isShapeToMesh = (_id = 0);
62 _meshDS = new SMESHDS_Mesh( _id, true );
65 // --------------------------------------------------------------------------
67 * \brief Edge of viscous layers; goes from a grid node by normal to boundary
71 std::vector< SMESH_NodeXYZ > _nodes;
72 double _uv[2]; // position of TgtNode() on geometry
75 const SMDS_MeshNode* TgtNode() const { return _nodes.back().Node(); }
77 typedef NCollection_DataMap< const SMDS_MeshNode*, VLEdge* > TNode2VLEdge;
78 // --------------------------------------------------------------------------
80 * \brief VLEdge's of one shape
84 std::vector< VLEdge > _edges;
85 TopoDS_Shape _initShape;
86 TopoDS_Shape _offsetShape;
89 bool _toCheckCoinc; // to check if nodes on shape boundary coincide
92 //================================================================================
94 * \brief Project a point on offset FACE to EDGEs of an initial FACE
95 * \param [in] offP - point to project
96 * \param [in] initFace - FACE to project on
97 * \return gp_Pnt - projection point
99 //================================================================================
101 gp_Pnt projectToWire( const gp_Pnt& offP,
102 const TopoDS_Face& initFace,
105 double minDist = Precision::Infinite();
106 const double tol = Precision::Confusion();
107 gp_Pnt proj, projction;
109 ShapeAnalysis_Curve projector;
111 for ( TopExp_Explorer eExp( initFace, TopAbs_EDGE ); eExp.More(); eExp.Next() )
113 BRepAdaptor_Curve initCurve( TopoDS::Edge( eExp.Current() ));
114 //const double f = initCurve.FirstParameter(), l = initCurve.LastParameter();
115 double dist = projector.Project( initCurve, offP, tol, proj, param, /*adjustToEnds=*/false );
116 if ( dist < minDist )
122 uv.SetCoord(0.,0.); // !!!!!!!
126 //================================================================================
128 * \brief Project nodes from offset FACE to initial FACE
129 * \param [inout] theEOS - VL edges on a geom FACE
130 * \param [inout] theOffsetMDS - offset mesh to fill in
132 //================================================================================
134 void projectToFace( VLEdgesOnShape & theEOS,
135 SMESHDS_Mesh* theOffsetMDS,
136 TNode2VLEdge & theN2E )
138 SMESHDS_SubMesh* sm = theOffsetMDS->MeshElements( theEOS._offsetShape );
139 if ( !sm || sm->NbElements() == 0 || sm->NbNodes() == 0 )
141 theEOS._edges.resize( sm->NbNodes() );
143 const TopoDS_Face& initFace = TopoDS::Face( theEOS._initShape );
144 ShapeAnalysis_Surface projector( BRep_Tool::Surface( initFace ));
145 const double clsfTol = 1e2 * BRep_Tool::MaxTolerance( initFace, TopAbs_VERTEX );
146 BRepTopAdaptor_FClass2d classifier( initFace, clsfTol );
148 const double tol = Precision::Confusion();
149 //const double f = initCurve.FirstParameter(), l = initCurve.LastParameter();
153 for ( SMDS_NodeIteratorPtr nIt = sm->GetNodes(); nIt->more(); ++iN )
155 SMESH_NodeXYZ offP = nIt->next();
156 gp_Pnt2d uv = projector.ValueOfUV( offP, tol );
157 TopAbs_State st = classifier.Perform( uv );
158 if ( st == TopAbs_IN || st == TopAbs_ON )
160 proj = projector.Value( uv );
164 proj = projectToWire( offP, initFace, uv );
166 if ( st == TopAbs_ON || st == TopAbs_OUT )
167 theEOS._toCheckCoinc = true;
169 VLEdge & vlEdge = theEOS._edges[ iN ];
170 vlEdge._nodes.push_back( offP.Node() );
171 vlEdge._nodes.push_back( theOffsetMDS->AddNode( proj.X(), proj.Y(), proj.Z() ));
172 vlEdge._uv[0] = uv.X();
173 vlEdge._uv[1] = uv.Y();
174 //vlEdge._length = proj.Distance( offP );
176 theN2E.Bind( offP.Node(), &vlEdge );
182 //================================================================================
184 * \brief Project nodes from offset EDGE to initial EDGE
185 * \param [inout] theEOS - VL edges on a geom EDGE
186 * \param [inout] theOffsetMDS - offset mesh to add new nodes to
188 //================================================================================
190 void projectToEdge( VLEdgesOnShape & theEOS,
191 SMESHDS_Mesh* theOffsetMDS,
192 TNode2VLEdge & theN2E )
194 SMESHDS_SubMesh* sm = theOffsetMDS->MeshElements( theEOS._offsetShape );
195 if ( !sm || sm->NbElements() == 0 )
197 theEOS._edges.resize( sm->NbNodes() );
199 ShapeAnalysis_Curve projector;
200 BRepAdaptor_Curve initCurve( TopoDS::Edge( theEOS._initShape ));
201 const double tol = Precision::Confusion();
202 const double f = initCurve.FirstParameter(), l = initCurve.LastParameter();
207 for ( SMDS_NodeIteratorPtr nIt = sm->GetNodes(); nIt->more(); ++iN )
209 SMESH_NodeXYZ offP = nIt->next();
210 double dist = projector.Project( initCurve, offP, tol, proj, param, /*adjustToEnds=*/false );
211 bool paramOK = ( f < param && param < l );
216 else if ( param > l )
218 theEOS._toCheckCoinc = true;
220 proj = initCurve.Value( param );
222 VLEdge & vlEdge = theEOS._edges[ iN ];
223 vlEdge._nodes.push_back( offP.Node() );
224 vlEdge._nodes.push_back( theOffsetMDS->AddNode( proj.X(), proj.Y(), proj.Z() ));
225 vlEdge._uv[0] = param;
226 vlEdge._length = paramOK ? dist : proj.Distance( offP );
228 theN2E.Bind( offP.Node(), &vlEdge );
233 //================================================================================
235 * \brief Compute heights of viscous layers and finds FACEs with VL
236 * \param [in] hyp - viscous layers hypothesis
237 * \param [in] mesh - the main mesh
238 * \param [in] shape - the main shape
239 * \param [out] vlH - heights of viscous layers to compute
240 * \param [out] shapesWVL - IDs of shapes with VL
241 * \return bool - isOK
243 //================================================================================
245 void computeVLHeight( const StdMeshers_ViscousLayers* hyp,
246 std::vector< double > & vlH )
248 const double T = hyp->GetTotalThickness();
249 const double f = hyp->GetStretchFactor();
250 const int N = hyp->GetNumberLayers();
251 const double h0 = hyp->Get1stLayerThickness( T, f, N );
253 vlH.reserve( hyp->GetNumberLayers() );
255 for ( double h = h0 * f; (int)vlH.size() < N-1; h *= f )
256 vlH.push_back( h + vlH.back() );
260 //================================================================================
262 * \brief Create intermediate nodes on VLEdge
263 * \param [inout] vlEdge - VLEdge to divide
264 * \param [in] vlH - thicknesses of layers
265 * \param [inout] mesh - the mesh no fill in
267 //================================================================================
269 void divideVLEdge( VLEdge* vlEdge,
270 const std::vector< double >& vlH,
273 SMESH_NodeXYZ lastNode = vlEdge->_nodes.back();
274 SMESH_NodeXYZ frstNode = vlEdge->_nodes[0];
275 gp_XYZ dir = frstNode - lastNode;
277 vlEdge->_nodes.resize( vlH.size() + 1 );
279 for ( size_t i = 1; i < vlH.size(); ++i )
281 double r = vlH[ i-1 ] / vlH.back();
282 gp_XYZ p = lastNode + dir * r;
283 vlEdge->_nodes[ vlH.size() - i ] = mesh->AddNode( p.X(), p.Y(), p.Z() );
285 vlEdge->_nodes.back() = lastNode;
288 //================================================================================
290 * \brief Create a polyhedron from nodes of VLEdge's
291 * \param [in] edgesVec - the VLEdge's
292 * \param [in] vNodes - node buffer
293 * \param [in] editor - editor of offset mesh
295 //================================================================================
297 bool makePolyhedron( const std::vector< VLEdge*> & edgesVec,
298 std::vector< const SMDS_MeshNode* > & vNodes,
299 SMESH_MeshEditor& editor,
300 SMESH_MeshEditor::ElemFeatures elemType)
302 elemType.SetPoly( true );
303 elemType.myPolyhedQuantities.clear();
304 elemType.myNodes.clear();
306 // add facets of walls
307 size_t nbBaseNodes = edgesVec.size();
308 for ( size_t iN = 0; iN < nbBaseNodes; ++iN )
310 VLEdge* e0 = edgesVec[ iN ];
311 VLEdge* e1 = edgesVec[( iN + 1 ) % nbBaseNodes ];
312 size_t nbN0 = e0->_nodes.size();
313 size_t nbN1 = e1->_nodes.size();
316 for ( size_t i = 1; i < nbN0; ++i )
318 elemType.myNodes.push_back( e1->_nodes[ i - 1 ].Node());
319 elemType.myNodes.push_back( e1->_nodes[ i ].Node());
320 elemType.myNodes.push_back( e0->_nodes[ i ].Node());
321 elemType.myNodes.push_back( e0->_nodes[ i - 1 ].Node());
322 elemType.myPolyhedQuantities.push_back( 4 );
327 for ( size_t i = 0; i < nbN1; ++i )
328 elemType.myNodes.push_back( e1->_nodes[ i ].Node() );
329 for ( size_t i = 0; i < nbN0; ++i )
330 elemType.myNodes.push_back( e0->_nodes[ nbN0 - 1 - i ].Node() );
331 elemType.myPolyhedQuantities.push_back( nbN0 + nbN1 );
337 for ( size_t iN = 0; iN < nbBaseNodes; ++iN )
339 VLEdge* e0 = edgesVec[ iN ];
340 elemType.myNodes.push_back( e0->_nodes.back().Node() );
341 vNodes.push_back( e0->_nodes[ 0 ].Node());
343 elemType.myPolyhedQuantities.push_back( nbBaseNodes );
345 // add facets of bottom
346 elemType.myNodes.insert( elemType.myNodes.end(), vNodes.rbegin(), vNodes.rend() );
347 elemType.myPolyhedQuantities.push_back( nbBaseNodes );
349 const SMDS_MeshElement* vol = editor.AddElement( elemType.myNodes, elemType );
350 vol->setIsMarked( true ); // to add to group
355 //================================================================================
357 * \brief Transform prism nodal connectivity to that of polyhedron
358 * \param [inout] nodes - nodes of prism, return nodes of polyhedron
359 * \param [inout] elemType - return quantities of polyhedron,
361 //================================================================================
363 void prism2polyhedron( std::vector< const SMDS_MeshNode* > & nodes,
364 SMESH_MeshEditor::ElemFeatures & elemType )
366 elemType.myPolyhedQuantities.clear();
367 elemType.myNodes.clear();
370 size_t nbBaseNodes = nodes.size() / 2;
371 for ( size_t i = 0; i < nbBaseNodes; ++i )
373 int iNext = ( i + 1 ) % nbBaseNodes;
374 elemType.myNodes.push_back( nodes[ iNext ]);
375 elemType.myNodes.push_back( nodes[ i ]);
376 elemType.myNodes.push_back( nodes[ i + nbBaseNodes ]);
377 elemType.myNodes.push_back( nodes[ iNext + nbBaseNodes ]);
378 elemType.myPolyhedQuantities.push_back( 4 );
382 elemType.myNodes.insert( elemType.myNodes.end(), nodes.begin(), nodes.begin() + nbBaseNodes );
383 elemType.myPolyhedQuantities.push_back( nbBaseNodes );
386 elemType.myNodes.insert( elemType.myNodes.end(), nodes.rbegin(), nodes.rbegin() + nbBaseNodes );
387 elemType.myPolyhedQuantities.push_back( nbBaseNodes );
389 nodes.swap( elemType.myNodes );
392 //================================================================================
394 * \brief Create prisms from faces
395 * \param [in] theEOS - shape to treat
396 * \param [inout] theMesh - offset mesh to fill in
398 //================================================================================
400 bool makePrisms( VLEdgesOnShape & theEOS,
402 TNode2VLEdge & theN2E )
404 SMESHDS_SubMesh* sm = theMesh->GetMeshDS()->MeshElements( theEOS._offsetShape );
405 if ( !sm || sm->NbElements() == 0 )
408 SMESH_MeshEditor editor( theMesh );
409 SMESH_MeshEditor::ElemFeatures volumElem( SMDSAbs_Volume );
411 std::vector< const SMDS_MeshNode* > vNodes;
412 std::vector< VLEdge*> edgesVec;
413 for ( SMDS_ElemIteratorPtr eIt = sm->GetElements(); eIt->more(); )
415 const SMDS_MeshElement* face = eIt->next();
416 if ( face->GetType() != SMDSAbs_Face )
419 const int nbNodes = face->NbCornerNodes();
420 edgesVec.resize( nbNodes );
421 vNodes.resize( nbNodes * 2 );
422 int maxNbLayer = 0, minNbLayer = IntegerLast();
423 for ( int i = 0; i < nbNodes; ++i )
425 const SMDS_MeshNode* n = face->GetNode( i );
426 if ( !theN2E.IsBound( n ))
428 edgesVec[ i ] = theN2E( n );
429 maxNbLayer = Max( maxNbLayer, edgesVec[i]->_nodes.size() - 1 );
430 minNbLayer = Min( minNbLayer, edgesVec[i]->_nodes.size() - 1 );
432 if ( maxNbLayer == minNbLayer )
434 size_t nbPrism = edgesVec[0]->_nodes.size() - 1;
435 for ( size_t iP = 0; iP < nbPrism; ++iP )
437 vNodes.resize( nbNodes * 2 );
438 for ( int i = 0; i < nbNodes; ++i )
440 vNodes[ i ] = edgesVec[ i ]->_nodes[ iP + 1 ].Node();
441 vNodes[ i + nbNodes ] = edgesVec[ i ]->_nodes[ iP ].Node();
443 volumElem.SetPoly( nbNodes > 4 );
444 if ( volumElem.myIsPoly )
445 prism2polyhedron( vNodes, volumElem );
447 if ( const SMDS_MeshElement* vol = editor.AddElement( vNodes, volumElem ))
448 vol->setIsMarked( true ); // to add to group
451 else // at inlet/outlet
453 makePolyhedron( edgesVec, vNodes, editor, volumElem );
455 editor.ClearLastCreated();
457 // move the face to the top of prisms, on mesh boundary
458 //theMesh->GetMeshDS()->ChangeElementNodes( face, fNodes.data(), nbNodes );
463 //================================================================================
465 * \brief Create faces from edges
466 * \param [inout] theEOS - shape to treat
467 * \param [inout] theMesh - offset mesh to fill in
468 * \param [inout] theN2E - map of node to VLEdge
469 * \param [inout] theFaceID - ID of WOVL FACE for new faces to set on
472 //================================================================================
474 bool makeFaces( VLEdgesOnShape & theEOS,
476 TNode2VLEdge & theN2E,
477 const TGeomID theFaceID)
479 SMESHDS_SubMesh* sm = theMesh->GetMeshDS()->MeshElements( theEOS._offsetShape );
480 if ( !sm || sm->NbElements() == 0 )
483 SMESH_MeshEditor editor( theMesh );
484 SMESH_MeshEditor::ElemFeatures elem2D( SMDSAbs_Face );
486 std::vector< const SMDS_MeshNode* > fNodes( 4 );
487 std::vector<const SMDS_MeshElement *> foundVolum;
488 std::vector< VLEdge*> edgesVec;
489 for ( SMDS_ElemIteratorPtr eIt = sm->GetElements(); eIt->more(); )
491 const SMDS_MeshElement* edge = eIt->next();
492 if ( edge->GetType() != SMDSAbs_Edge )
495 const int nbNodes = 2; //edge->NbCornerNodes();
496 edgesVec.resize( nbNodes );
497 fNodes.resize( nbNodes * 2 );
498 for ( int i = 0; i < nbNodes; ++i )
500 const SMDS_MeshNode* n = edge->GetNode( i );
501 if ( !theN2E.IsBound( n ))
503 edgesVec[ i ] = theN2E( n );
505 size_t nbFaces = edgesVec[0]->_nodes.size() - 1;
506 for ( size_t iF = 0; iF < nbFaces; ++iF )
508 fNodes[ 0 ] = edgesVec[ 0 ]->_nodes[ iF ].Node();
509 fNodes[ 1 ] = edgesVec[ 1 ]->_nodes[ iF ].Node();
510 fNodes[ 2 ] = edgesVec[ 1 ]->_nodes[ iF + 1 ].Node();
511 fNodes[ 3 ] = edgesVec[ 0 ]->_nodes[ iF + 1 ].Node();
513 if ( const SMDS_MeshElement* face = editor.AddElement( fNodes, elem2D ))
515 theMesh->GetMeshDS()->SetMeshElementOnShape( face, theFaceID );
517 if ( SMDS_Mesh::GetElementsByNodes( fNodes, foundVolum, SMDSAbs_Volume ))
519 TIDSortedElemSet faces = { face }, volumes = { foundVolum[0] };
520 editor.Reorient2DBy3D( faces, volumes, /*outside=*/true );
524 editor.ClearLastCreated();
526 // move the face to the top of prisms, on mesh boundary
527 //theMesh->GetMeshDS()->ChangeElementNodes( face, fNodes.data(), nbNodes );
532 //================================================================================
534 * \brief Append the offset mesh to the initial one
536 //================================================================================
538 void copyMesh( SMESH_Mesh* theFromMesh,
539 SMESH_Mesh* theToMesh,
542 SMESHDS_Mesh* offsetMDS = theFromMesh->GetMeshDS();
543 SMESHDS_Mesh* initMDS = theToMesh->GetMeshDS();
545 const smIdType nShift = initMDS->NbNodes();
547 for ( SMDS_NodeIteratorPtr nIt = offsetMDS->nodesIterator(); nIt->more(); )
549 SMESH_NodeXYZ pOff = nIt->next();
550 const SMDS_MeshNode* n = initMDS->AddNodeWithID( pOff.X(), pOff.Y(), pOff.Z(),
551 nShift + pOff.Node()->GetID() );
552 initMDS->SetNodeInVolume( n, theSolidID );
555 SMESH_MeshEditor editor( theToMesh );
556 SMESH_MeshEditor::ElemFeatures elemType;
557 std::vector<smIdType> nIniIDs;
559 for ( SMDS_ElemIteratorPtr eIt = offsetMDS->elementsIterator(); eIt->more(); )
561 const SMDS_MeshElement* offElem = eIt->next();
562 const int nbNodes = offElem->NbNodes();
563 nIniIDs.resize( nbNodes );
564 for ( int i = 0; i < nbNodes; ++i )
566 const SMDS_MeshNode* offNode = offElem->GetNode( i );
567 nIniIDs[ i ] = offNode->GetID() + nShift;
569 elemType.Init( offElem, /*basicOnly=*/false );
570 const SMDS_MeshElement* iniElem = editor.AddElement( nIniIDs, elemType );
571 initMDS->SetMeshElementOnShape( iniElem, theSolidID );
572 iniElem->setIsMarked( offElem->isMarked() );
573 editor.ClearLastCreated();
578 //================================================================================
580 * \brief set elements of layers boundary to sub-meshes
581 * \param [in] theEOS - vlEdges of a sub-shape
582 * \param [inout] theMesh - the mesh
583 * \param [in] theN2E - map of node to vlEdge
584 * \param [in] theNShift - nb of nodes in the mesh before adding the offset mesh
586 //================================================================================
588 void setBnd2Sub( VLEdgesOnShape & theEOS,
589 SMESH_Mesh* theInitMesh,
590 SMESH_Mesh* theOffsMesh,
591 const TNode2VLEdge & theN2E,
592 const smIdType theNShift,
593 TIDSortedNodeSet& theNodesToCheckCoinc )
595 SMESHDS_Mesh* offsetMDS = theOffsMesh->GetMeshDS();
596 SMESHDS_Mesh* initMDS = theInitMesh->GetMeshDS();
598 SMESHDS_SubMesh* sm = offsetMDS->MeshElements( theEOS._offsetShape );
599 if ( !sm || ( sm->NbElements() + sm->NbNodes() == 0 ))
602 for ( SMDS_NodeIteratorPtr nIt = sm->GetNodes(); nIt->more(); )
604 const SMDS_MeshNode* offNode = nIt->next();
605 VLEdge* vlEd = theN2E( offNode );
606 const SMDS_MeshNode* nOffBnd = vlEd->_nodes.back().Node();
607 const SMDS_MeshNode* nIniBnd = initMDS->FindNode( nOffBnd->GetID() + theNShift );
608 initMDS->UnSetNodeOnShape( nIniBnd );
610 switch( theEOS._initShape.ShapeType() ) {
611 case TopAbs_FACE: initMDS->SetNodeOnFace( nIniBnd, TopoDS::Face( theEOS._initShape ),
612 vlEd->_uv[0], vlEd->_uv[1] );
614 case TopAbs_EDGE: initMDS->SetNodeOnEdge( nIniBnd, TopoDS::Edge( theEOS._initShape ),
617 case TopAbs_VERTEX: initMDS->SetNodeOnVertex( nIniBnd, TopoDS::Vertex( theEOS._initShape ));
623 std::vector< const SMDS_MeshNode* > iniNodes, iniNodesBnd;
624 std::vector< VLEdge*> edgesVec;
625 for ( SMDS_ElemIteratorPtr eIt = sm->GetElements(); eIt->more(); )
627 const SMDS_MeshElement* offElem = eIt->next();
628 if ( offElem->GetType() != SMDSAbs_Face &&
629 offElem->GetType() != SMDSAbs_Edge )
632 const int nbNodes = offElem->NbCornerNodes();
633 iniNodes.resize( nbNodes );
634 iniNodesBnd.resize( nbNodes );
635 for ( int i = 0; i < nbNodes; ++i )
637 const SMDS_MeshNode* nOff = offElem->GetNode( i );
638 const SMDS_MeshNode* nIni = initMDS->FindNode( nOff->GetID() + theNShift );
639 iniNodes[ i ] = nIni;
640 if ( !theN2E.IsBound( nOff ))
642 VLEdge* vlEd = theN2E( nOff );
643 const SMDS_MeshNode* nOffBnd = vlEd->_nodes.back().Node();
644 const SMDS_MeshNode* nIniBnd = initMDS->FindNode( nOffBnd->GetID() + theNShift );
645 iniNodesBnd[ i ] = nIniBnd;
647 if ( const SMDS_MeshElement* iniElem = initMDS->FindElement( iniNodes ))
649 initMDS->UnSetElementOnShape( iniElem );
650 initMDS->SetMeshElementOnShape( iniElem, theEOS._initShape );
652 // move the face to the top of prisms, on init mesh boundary
653 initMDS->ChangeElementNodes( iniElem, iniNodesBnd.data(), nbNodes );
655 if ( theEOS._toCheckCoinc )
656 theNodesToCheckCoinc.insert( iniNodesBnd.begin(), iniNodesBnd.end() );
662 //================================================================================
664 * \brief set elements of WOVL face to sub-meshes
665 * \param [in] theEOS - vlEdges of a sub-shape
666 * \param [inout] theMesh - the mesh
667 * \param [in] theN2E - map of node to vlEdge
668 * \param [in] theNShift - nb of nodes in the mesh before adding the offset mesh
670 //================================================================================
672 void setBnd2FVWL( VLEdgesOnShape & theEOS,
673 SMESH_Mesh* theInitMesh,
674 SMESH_Mesh* theOffsMesh,
675 //const TNode2VLEdge & theN2E,
676 const smIdType theNShift )
678 SMESHDS_Mesh* offsetMDS = theOffsMesh->GetMeshDS();
679 SMESHDS_Mesh* initMDS = theInitMesh->GetMeshDS();
681 SMESHDS_SubMesh* sm = offsetMDS->MeshElements( theEOS._offsetShape );
682 if ( !sm || ( sm->NbElements() + sm->NbNodes() == 0 ))
685 for ( SMDS_NodeIteratorPtr nIt = sm->GetNodes(); nIt->more(); )
687 const SMDS_MeshNode* offNode = nIt->next();
688 const SMDS_MeshNode* iniNode = initMDS->FindNode( offNode->GetID() + theNShift );
689 initMDS->UnSetNodeOnShape( iniNode );
691 switch( theEOS._initShape.ShapeType() ) {
692 case TopAbs_FACE: initMDS->SetNodeOnFace( iniNode, TopoDS::Face( theEOS._initShape ));
694 case TopAbs_EDGE: initMDS->SetNodeOnEdge( iniNode, TopoDS::Edge( theEOS._initShape ));
696 case TopAbs_VERTEX: initMDS->SetNodeOnVertex( iniNode, TopoDS::Vertex( theEOS._initShape ));
702 std::vector< const SMDS_MeshNode* > iniNodes;
703 for ( SMDS_ElemIteratorPtr eIt = sm->GetElements(); eIt->more(); )
705 const SMDS_MeshElement* offElem = eIt->next();
706 // if ( offElem->GetType() != SMDSAbs_Face )
709 int nbNodes = offElem->NbCornerNodes();
710 iniNodes.resize( nbNodes );
711 for ( int i = 0; i < nbNodes; ++i )
713 const SMDS_MeshNode* nOff = offElem->GetNode( i );
714 const SMDS_MeshNode* nIni = initMDS->FindNode( nOff->GetID() + theNShift );
715 iniNodes[ i ] = nIni;
717 if ( const SMDS_MeshElement* iniElem = initMDS->FindElement( iniNodes ))
719 initMDS->UnSetElementOnShape( iniElem );
720 initMDS->SetMeshElementOnShape( iniElem, theEOS._initShapeID );
722 for ( const SMDS_MeshNode* n : iniNodes )
723 if ( n->GetPosition()->GetDim() == 3 )
725 initMDS->UnSetNodeOnShape( n );
726 if ( theEOS._initShape.ShapeType() == TopAbs_FACE )
727 initMDS->SetNodeOnFace( n, theEOS._initShapeID );
729 initMDS->SetNodeOnEdge( n, theEOS._initShapeID );
738 //================================================================================
740 * \brief Create a temporary mesh used to hold elements of the offset shape
742 //================================================================================
744 SMESH_Mesh* StdMeshers_Cartesian_VL::ViscousBuilder::MakeOffsetMesh()
749 //================================================================================
753 //================================================================================
755 StdMeshers_Cartesian_VL::ViscousBuilder::ViscousBuilder( const StdMeshers_ViscousLayers* hyp,
756 const SMESH_Mesh & mainMesh,
757 const TopoDS_Shape & mainShape)
758 : _hyp( hyp ), _offsetMesh( new TmpMesh )
760 // find shapes with viscous layers
761 TopTools_IndexedMapOfShape faces;
762 TopExp::MapShapes( mainShape, TopAbs_FACE, faces );
763 const SMESHDS_Mesh* iniMDS = mainMesh.GetMeshDS();
765 if ( hyp->IsToIgnoreShapes() )
767 for ( int i = 1; i <= faces.Size(); ++i )
768 _shapesWVL.insert( iniMDS->ShapeToIndex( faces( i )));
770 for ( TGeomID& sID : hyp->GetBndShapes() )
771 _shapesWVL.erase( sID );
775 for ( TGeomID& sID : hyp->GetBndShapes() )
776 _shapesWVL.insert( sID );
779 for ( TGeomID fID : _shapesWVL )
781 const TopoDS_Shape & face = iniMDS->IndexToShape( fID );
782 for ( TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next() )
783 _shapesWVL.insert( iniMDS->ShapeToIndex( exp.Current() ));
784 for ( TopExp_Explorer exp( face, TopAbs_VERTEX ); exp.More(); exp.Next() )
785 _shapesWVL.insert( iniMDS->ShapeToIndex( exp.Current() ));
788 // find EDGE and adjacent FACEs w/o VL
789 for ( int i = 1; i <= faces.Size(); ++i )
791 const TopoDS_Shape& face = faces( i );
792 TGeomID fID = iniMDS->ShapeToIndex( face );
793 if ( _shapesWVL.count( fID ))
795 for ( TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next() )
796 _edge2facesWOVL[ iniMDS->ShapeToIndex( exp.Current() )].push_back( fID );
801 //================================================================================
805 //================================================================================
807 StdMeshers_Cartesian_VL::ViscousBuilder::~ViscousBuilder()
809 delete _offsetMesh; _offsetMesh = 0;
812 //================================================================================
814 * \brief Create an offset shape from a given one
815 * \param [in] theShape - input shape
816 * \param [in] theMesh - main mesh
817 * \param [out] theError - error description
818 * \return TopoDS_Shape - result offset shape
820 //================================================================================
823 StdMeshers_Cartesian_VL::ViscousBuilder::MakeOffsetShape(const TopoDS_Shape & theShape,
824 SMESH_Mesh & theMesh,
825 std::string & theError )
827 double offset = -_hyp->GetTotalThickness();
828 double tol = Precision::Confusion();
829 _makeOffset.Initialize( theShape, offset, tol, BRepOffset_Skin, /*Intersection=*/false,
830 /*selfInter=*/false, GeomAbs_Intersection );
831 // exclude inlet FACEs
832 SMESHDS_Mesh* meshDS = theMesh.GetMeshDS();
833 for ( TopExp_Explorer fEx( theShape, TopAbs_FACE ); fEx.More(); fEx.Next() )
835 TGeomID fID = meshDS->ShapeToIndex( fEx.Current() );
836 if ( !_shapesWVL.count( fID ))
837 _makeOffset.SetOffsetOnFace( TopoDS::Face( fEx.Current()), 0 );
840 _makeOffset.MakeOffsetShape();
841 if ( _makeOffset.IsDone() )
843 _offsetShape = _makeOffset.Shape();
844 SMESH_MesherHelper::WriteShape( _offsetShape );////
846 _offsetMesh->ShapeToMesh( _offsetShape );
847 _offsetMesh->GetSubMesh( _offsetShape )->DependsOn();
851 switch ( _makeOffset.Error() )
853 case BRepOffset_NoError:
854 theError = "OK. Offset performed successfully.";break;
855 case BRepOffset_BadNormalsOnGeometry:
856 theError = "Degenerated normal on input data.";break;
857 case BRepOffset_C0Geometry:
858 theError = "C0 continuity of input data.";break;
859 case BRepOffset_NullOffset:
860 theError = "Null offset of all faces.";break;
861 case BRepOffset_NotConnectedShell:
862 theError = "Incorrect set of faces to remove, the remaining shell is not connected.";break;
863 case BRepOffset_CannotTrimEdges:
864 theError = "Can not trim edges.";break;
865 case BRepOffset_CannotFuseVertices:
866 theError = "Can not fuse vertices.";break;
867 case BRepOffset_CannotExtentEdge:
868 theError = "Can not extent edge.";break;
870 theError = "operation not done.";
872 theError = "BRepOffset_MakeOffset error: " + theError;
874 return TopoDS_Shape();
877 //================================================================================
879 * \brief Return a sub-shape of the offset shape generated from a given initial sub-shape
881 //================================================================================
883 TopoDS_Shape StdMeshers_Cartesian_VL::ViscousBuilder::getOffsetSubShape( const TopoDS_Shape& S )
885 const TopTools_ListOfShape& newShapes = _makeOffset.Generated( S );
886 for ( const TopoDS_Shape& ns : newShapes )
887 if ( ns.ShapeType() == S.ShapeType() )
889 return TopoDS_Shape();
892 //================================================================================
894 * \brief Create prismatic mesh between _offsetShape and theShape
895 * \param [out] theMesh - mesh to fill in
896 * \param [in] theShape - initial shape
897 * \return bool - is Ok
899 //================================================================================
901 bool StdMeshers_Cartesian_VL::ViscousBuilder::MakeViscousLayers( SMESH_Mesh & theMesh,
902 const TopoDS_Shape & theShape )
904 SMESHDS_Mesh* offsetMDS = _offsetMesh->GetMeshDS();
905 SMESHDS_Mesh* initMDS = theMesh.GetMeshDS();
906 offsetMDS->SetAllCellsNotMarked();
908 // Compute heights of viscous layers
909 std::vector< double > vlH;
910 computeVLHeight( _hyp, vlH );
912 std::vector< VLEdgesOnShape > edgesOnShape;
913 edgesOnShape.reserve( offsetMDS->MaxShapeIndex() + 1 );
916 // loop on sub-shapes to project nodes from offset boundary to initial boundary
917 TopAbs_ShapeEnum types[3] = { TopAbs_VERTEX, TopAbs_EDGE, TopAbs_FACE };
918 for ( TopAbs_ShapeEnum shType : types )
920 TopTools_IndexedMapOfShape shapes;
921 TopExp::MapShapes( theShape, shType, shapes );
922 for ( int i = 1; i <= shapes.Size(); ++i )
924 edgesOnShape.resize( edgesOnShape.size() + 1 );
925 VLEdgesOnShape& EOS = edgesOnShape.back();
927 EOS._initShape = shapes( i );
928 EOS._offsetShape = getOffsetSubShape( EOS._initShape );
929 EOS._initShapeID = initMDS->ShapeToIndex( EOS._initShape );
930 EOS._hasVL = _shapesWVL.count( EOS._initShapeID );
931 EOS._toCheckCoinc = false;
935 // project boundary nodes of offset mesh to boundary of init mesh
936 // (new nodes are created in the offset mesh)
937 switch( EOS._offsetShape.ShapeType() ) {
940 EOS._edges.resize( 1 );
941 EOS._edges[0]._nodes.resize( 2 );
942 EOS._edges[0]._nodes[0] = SMESH_Algo::VertexNode( TopoDS::Vertex( EOS._offsetShape ),
944 gp_Pnt offP = BRep_Tool::Pnt( TopoDS::Vertex( EOS._initShape ));
945 EOS._edges[0]._nodes[1] = offsetMDS->AddNode( offP.X(), offP.Y(), offP.Z() );
946 //EOS._edges[0]._length = offP.Distance( EOS._edges[0]._nodes[0] );
947 n2e.Bind( EOS._edges[0]._nodes[0].Node(), & EOS._edges[0] );
952 projectToEdge( EOS, offsetMDS, n2e );
957 projectToFace( EOS, offsetMDS, n2e );
963 // create nodes of layers
964 if ( _hyp->GetNumberLayers() > 1 )
966 //if ( _shapesWVL.count( EOS._initShapeID ))
967 for ( size_t i = 0; i < EOS._edges.size(); ++i )
969 divideVLEdge( &EOS._edges[ i ], vlH, offsetMDS );
973 } // loop on shape types
976 bool prismsOk = true;
977 for ( size_t i = 0; i < edgesOnShape.size(); ++i )
979 VLEdgesOnShape& EOS = edgesOnShape[ i ];
980 if ( EOS._initShape.ShapeType() == TopAbs_FACE && EOS._hasVL )
982 if ( !makePrisms( EOS, _offsetMesh, n2e ))
988 // create faces on FACEs WOVL
989 for ( size_t i = 0; i < edgesOnShape.size(); ++i )
991 VLEdgesOnShape& EOS = edgesOnShape[ i ];
992 if ( EOS._initShape.ShapeType() == TopAbs_EDGE && EOS._hasVL )
994 auto e2f = _edge2facesWOVL.find( EOS._initShapeID );
995 if ( e2f != _edge2facesWOVL.end() && !e2f->second.empty() )
997 TopoDS_Shape f = initMDS->IndexToShape( e2f->second[0] );
998 TopoDS_Shape f2 = getOffsetSubShape( f );
999 //cout << e2f->second[0] << " OFF " << offsetMDS->ShapeToIndex( f2 ) << endl;
1000 makeFaces( EOS, _offsetMesh, n2e, offsetMDS->ShapeToIndex( f2 ) );
1006 // copy offset mesh to the main one
1007 initMDS->Modified();
1008 initMDS->CompactMesh();
1009 smIdType nShift = initMDS->NbNodes();
1010 TGeomID solidID = initMDS->ShapeToIndex( theShape );
1011 copyMesh( _offsetMesh, & theMesh, solidID );
1016 if ( SMESH_subMesh * sm = theMesh.GetSubMesh( theShape ))
1018 sm->GetComputeError() =
1019 SMESH_ComputeError::New( COMPERR_ALGO_FAILED,
1020 "Viscous layers construction error: bad mesh on offset geometry" );
1025 // set elements of layers boundary to sub-meshes
1026 TIDSortedNodeSet nodesToCheckCoinc;
1027 for ( size_t i = 0; i < edgesOnShape.size(); ++i )
1029 VLEdgesOnShape& EOS = edgesOnShape[ i ];
1031 setBnd2Sub( EOS, &theMesh, _offsetMesh, n2e, nShift, nodesToCheckCoinc );
1033 setBnd2FVWL( EOS, &theMesh, _offsetMesh, nShift );
1036 // merge coincident nodes
1037 SMESH_MeshEditor editor( &theMesh );
1038 SMESH_MeshEditor::TListOfListOfNodes nodesToMerge;
1039 editor.FindCoincidentNodes( nodesToCheckCoinc, vlH[0]/50., nodesToMerge, false );
1040 editor.MergeNodes( nodesToMerge );
1043 if ( !_hyp->GetGroupName().empty() )
1045 SMDS_MeshGroup* group = _hyp->CreateGroup( _hyp->GetGroupName(), theMesh, SMDSAbs_Volume );
1047 for ( SMDS_ElemIteratorPtr it = initMDS->elementsIterator(); it->more(); )
1049 const SMDS_MeshElement* el = it->next();
1050 if ( el->isMarked() )