From 1273dc58a7b618bcf106778e42ce6477acfa2f26 Mon Sep 17 00:00:00 2001 From: eap Date: Thu, 5 Jul 2012 09:39:10 +0000 Subject: [PATCH] 0021676: EDF 2283 NETGENPLUGIN: Improve Netgen 1D-2D-3D to generate pyramids in case where input 2D mesh includes quadrangles --- src/NETGENPlugin/NETGENPlugin_Mesher.cxx | 204 +++++++++++++++++++---- src/NETGENPlugin/NETGENPlugin_Mesher.hxx | 19 ++- 2 files changed, 183 insertions(+), 40 deletions(-) diff --git a/src/NETGENPlugin/NETGENPlugin_Mesher.cxx b/src/NETGENPlugin/NETGENPlugin_Mesher.cxx index 94a0447..c4320fc 100644 --- a/src/NETGENPlugin/NETGENPlugin_Mesher.cxx +++ b/src/NETGENPlugin/NETGENPlugin_Mesher.cxx @@ -43,28 +43,21 @@ #include #include #include -#include +#include -#include -#include +#include +#include #include #include -#include -#include #include -#include -#include #include -#include -#include #include #include #include #include #include #include -#include #include #include @@ -82,6 +75,9 @@ namespace netgen { extern volatile multithreadt multithread; } +#include +#include + using namespace nglib; using namespace std; @@ -211,7 +207,7 @@ void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_Hypothesis* hyp) mparams.secondorder = hyp->GetSecondOrder() ? 1 : 0; // quad-dominated surface meshing // only triangles are allowed for volumic mesh - if (!_isVolume) + //if (!_isVolume) mparams.quad = static_cast (hyp)->GetQuadAllowed() ? 1 : 0; _optimize = hyp->GetOptimize(); @@ -612,10 +608,11 @@ void NETGENPlugin_Mesher::RestrictLocalSize(netgen::Mesh& ngMesh, const gp_XYZ& */ //================================================================================ -bool NETGENPlugin_Mesher::fillNgMesh(const netgen::OCCGeometry& occgeom, +bool NETGENPlugin_Mesher::fillNgMesh(netgen::OCCGeometry& occgeom, netgen::Mesh& ngMesh, vector& nodeVec, - const list< SMESH_subMesh* > & meshedSM) + const list< SMESH_subMesh* > & meshedSM, + SMESH_ProxyMesh::Ptr proxyMesh) { TNode2IdMap nodeNgIdMap; for ( int i = 1; i < nodeVec.size(); ++i ) @@ -627,7 +624,7 @@ bool NETGENPlugin_Mesher::fillNgMesh(const netgen::OCCGeometry& occgeom, SMESH_MesherHelper helper (*_mesh); - int faceNgID = occgeom.fmap.Extent(); + int faceNgID = ngMesh.GetNFD(); list< SMESH_subMesh* >::const_iterator smIt, smEnd = meshedSM.end(); for ( smIt = meshedSM.begin(); smIt != smEnd; ++smIt ) @@ -636,7 +633,7 @@ bool NETGENPlugin_Mesher::fillNgMesh(const netgen::OCCGeometry& occgeom, if ( !visitedShapes.Add( sm->GetSubShape() )) continue; - SMESHDS_SubMesh * smDS = sm->GetSubMeshDS(); + const SMESHDS_SubMesh * smDS = sm->GetSubMeshDS(); if ( !smDS ) continue; switch ( sm->GetSubShape().ShapeType() ) @@ -799,19 +796,38 @@ bool NETGENPlugin_Mesher::fillNgMesh(const netgen::OCCGeometry& occgeom, // Find solids the geomFace bounds int solidID1 = 0, solidID2 = 0; - PShapeIteratorPtr solidIt = helper.GetAncestors( geomFace, *sm->GetFather(), TopAbs_SOLID); - while ( const TopoDS_Shape * solid = solidIt->next() ) + StdMeshers_QuadToTriaAdaptor* quadAdaptor = + dynamic_cast( proxyMesh.get() ); + if ( quadAdaptor ) { - int id = occgeom.somap.FindIndex ( *solid ); - if ( solidID1 && id != solidID1 ) solidID2 = id; - else solidID1 = id; + solidID1 = occgeom.somap.FindIndex( quadAdaptor->GetShape() ); + } + else + { + PShapeIteratorPtr solidIt = helper.GetAncestors( geomFace, *sm->GetFather(), TopAbs_SOLID); + while ( const TopoDS_Shape * solid = solidIt->next() ) + { + int id = occgeom.somap.FindIndex ( *solid ); + if ( solidID1 && id != solidID1 ) solidID2 = id; + else solidID1 = id; + } } - faceNgID++; - //_faceDescriptors[ faceNgID ].first = solidID1; - //_faceDescriptors[ faceNgID ].second = solidID2; // Add ng face descriptors of meshed faces + faceNgID++; ngMesh.AddFaceDescriptor (netgen::FaceDescriptor(faceNgID, solidID1, solidID2, 0)); + // if second oreder is required, even already meshed faces must be passed to NETGEN + int fID = occgeom.fmap.Add( geomFace ); + while ( fID < faceNgID ) // geomFace is already in occgeom.fmap, add a copy + fID = occgeom.fmap.Add( BRepBuilderAPI_Copy( geomFace, /*copyGeom=*/false )); + // Problem with the second order in a quadrangular mesh remains. + // 1) All quadrangles geberated by NETGEN are moved to an inexistent face + // by FillSMesh() (find AddFaceDescriptor) + // 2) Temporary triangles generated by StdMeshers_QuadToTriaAdaptor + // are on faces where quadrangles were. + // Due to these 2 points, wrong geom faces are used while conversion to qudratic + // of the mentioned above quadrangles and triangles + // Orient the face correctly in solidID1 (issue 0020206) bool reverse = false; if ( solidID1 ) { @@ -830,8 +846,11 @@ bool NETGENPlugin_Mesher::fillNgMesh(const netgen::OCCGeometry& occgeom, #ifdef DUMP_TRIANGLES cout << "SMESH face " << helper.GetMeshDS()->ShapeToIndex( geomFace ) - << " internal="<SetMeshElementOnShape(face, aFace); } - // create tetrahedra + // ------------------ + // Create tetrahedra + // ------------------ + for (i = 1; i <= nbVol; ++i) { const netgen::Element& elem = ngMesh.VolumeElement(i); @@ -1859,6 +1900,8 @@ bool NETGENPlugin_Mesher::Compute() return false; ngLib.setMesh(( Ng_Mesh*) ngMesh ); + ngMesh->ClearFaceDescriptors(); // we make descriptors our-self + if ( _simpleHyp ) { // Pass 1D simple parameters to NETGEN @@ -2063,14 +2106,38 @@ bool NETGENPlugin_Mesher::Compute() // --------------------- // generate volume mesh // --------------------- - // Fill ngMesh with nodes and faces of computed 2D submeshes - if ( !err && _isVolume && !meshedSM[ MeshDim_2D ].empty() ) + // Fill ngMesh with nodes and faces of computed 2D submeshes + if ( !err && _isVolume && ( !meshedSM[ MeshDim_2D ].empty() || mparams.quad )) { // load SMESH with computed segments and faces FillSMesh( occgeo, *ngMesh, initState, *_mesh, nodeVec, comment ); - // fill ng mesh + + // compute pyramids on quadrangles + SMESH_ProxyMesh::Ptr proxyMesh; + if ( _mesh->NbQuadrangles() > 0 ) + for ( int iS = 1; iS <= occgeo.somap.Extent(); ++iS ) + { + StdMeshers_QuadToTriaAdaptor* Adaptor = new StdMeshers_QuadToTriaAdaptor; + proxyMesh.reset( Adaptor ); + + int nbPyrams = _mesh->NbPyramids(); + Adaptor->Compute( *_mesh, occgeo.somap(iS) ); + if ( nbPyrams != _mesh->NbPyramids() ) + { + list< SMESH_subMesh* > quadFaceSM; + for (TopExp_Explorer face(occgeo.somap(iS), TopAbs_FACE); face.More(); face.Next()) + if ( Adaptor->GetProxySubMesh( face.Current() )) + { + quadFaceSM.push_back( _mesh->GetSubMesh( face.Current() )); + meshedSM[ MeshDim_2D ].remove( quadFaceSM.back() ); + } + fillNgMesh(occgeo, *ngMesh, nodeVec, quadFaceSM, proxyMesh); + } + } + // fill ngMesh with faces of sub-meshes err = ! ( fillNgMesh(occgeo, *ngMesh, nodeVec, meshedSM[ MeshDim_2D ])); initState = NETGENPlugin_ngMeshInfo(ngMesh); + //toPython( ngMesh, "/tmp/ngPython.py"); } if (!err && _isVolume) { @@ -2236,7 +2303,7 @@ bool NETGENPlugin_Mesher::Compute() for (int i = 1; i <= occgeo.somap.Extent(); i++) if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( occgeo.somap( i ))) { - bool smComputed = !sm->IsEmpty(); + bool smComputed = nbVol && !sm->IsEmpty(); if ( smComputed && internals.hasInternalVertexInSolid( sm->GetId() )) { int nbIntV = internals.getSolidsWithVertices().find( sm->GetId() )->second.size(); @@ -2247,7 +2314,7 @@ bool NETGENPlugin_Mesher::Compute() if ( !smComputed && ( !smError || smError->IsOK() )) { smError.reset( new SMESH_ComputeError( *error )); - if ( SMESH_Algo::GetMeshError( sm ) == SMESH_Algo::MEr_OK ) + if ( nbVol && SMESH_Algo::GetMeshError( sm ) == SMESH_Algo::MEr_OK ) smError->myName = COMPERR_WARNING; } pb3D = pb3D || ( smError && smError->IsKO() ); @@ -2543,6 +2610,75 @@ NETGENPlugin_Mesher::readErrors(const vector& nodeVec) return err; } +//================================================================================ +/*! + * \brief Write a python script creating an equivalent SALOME mesh. + * This is useful to see what mesh is passed as input for the next step of mesh + * generation (of mesh of higher dimension) + */ +//================================================================================ + +void NETGENPlugin_Mesher::toPython( const netgen::Mesh* ngMesh, + const std::string& pyFile) +{ + ofstream outfile(pyFile.c_str(), ios::out); + if ( !outfile ) return; + + outfile << "import smesh, SMESH" << endl + << "mesh = smesh.Mesh()" << endl << endl; + + using namespace netgen; + PointIndex pi; + for (pi = PointIndex::BASE; + pi < ngMesh->GetNP()+PointIndex::BASE; pi++) + { + outfile << "mesh.AddNode( "; + outfile << (*ngMesh)[pi](0) << ", "; + outfile << (*ngMesh)[pi](1) << ", "; + outfile << (*ngMesh)[pi](2) << ")" << endl; + } + + int nbDom = ngMesh->GetNDomains(); + for ( int i = 0; i < nbDom; ++i ) + outfile<< "grp" << i+1 << " = mesh.CreateEmptyGroup( SMESH.FACE, 'domain"<< i+1 << "')"<< endl; + + SurfaceElementIndex sei; + for (sei = 0; sei < ngMesh->GetNSE(); sei++) + { + outfile << "mesh.AddFace([ "; + Element2d sel = (*ngMesh)[sei]; + for (int j = 0; j < sel.GetNP(); j++) + outfile << sel[j] << ( j+1 < sel.GetNP() ? ", " : " ])"); + outfile << endl; + + if ((*ngMesh)[sei].GetIndex()) + { + if ( int dom1 = ngMesh->GetFaceDescriptor((*ngMesh)[sei].GetIndex ()).DomainIn()) + outfile << "grp"<< dom1 <<".Add([ " << (int)sei+1 << " ])" << endl; + if ( int dom2 = ngMesh->GetFaceDescriptor((*ngMesh)[sei].GetIndex ()).DomainOut()) + outfile << "grp"<< dom2 <<".Add([ " << (int)sei+1 << " ])" << endl; + } + } + + for (ElementIndex ei = 0; ei < ngMesh->GetNE(); ei++) + { + Element el = (*ngMesh)[ei]; + outfile << "mesh.AddVolume([ "; + for (int j = 0; j < el.GetNP(); j++) + outfile << el[j] << ( j+1 < el.GetNP() ? ", " : " ])"); + outfile << endl; + } + + for (int i = 1; i <= ngMesh->GetNSeg(); i++) + { + const Segment & seg = ngMesh->LineSegment (i); + outfile << "mesh.AddEdge([ " + << seg[0] << ", " + << seg[1] << " ])" << endl; + } + cout << "Write " << pyFile << endl; +} + //================================================================================ /*! * \brief Constructor of NETGENPlugin_ngMeshInfo diff --git a/src/NETGENPlugin/NETGENPlugin_Mesher.hxx b/src/NETGENPlugin/NETGENPlugin_Mesher.hxx index d1760a5..a525ea7 100644 --- a/src/NETGENPlugin/NETGENPlugin_Mesher.hxx +++ b/src/NETGENPlugin/NETGENPlugin_Mesher.hxx @@ -30,9 +30,11 @@ #define _NETGENPlugin_Mesher_HXX_ #include "NETGENPlugin_Defs.hxx" -#include "StdMeshers_FaceSide.hxx" -#include "SMDS_MeshElement.hxx" -#include "SMESH_Algo.hxx" + +#include +#include +#include +#include namespace nglib { #include @@ -103,16 +105,17 @@ class NETGENPLUGIN_EXPORT NETGENPlugin_Mesher static void RestrictLocalSize(netgen::Mesh& ngMesh, const gp_XYZ& p, const double size); static int FillSMesh(const netgen::OCCGeometry& occgeom, - const netgen::Mesh& ngMesh, + netgen::Mesh& ngMesh, const NETGENPlugin_ngMeshInfo& initState, SMESH_Mesh& sMesh, std::vector& nodeVec, SMESH_Comment& comment); - bool fillNgMesh(const netgen::OCCGeometry& occgeom, + bool fillNgMesh(netgen::OCCGeometry& occgeom, netgen::Mesh& ngMesh, std::vector& nodeVec, - const std::list< SMESH_subMesh* > & meshedSM); + const std::list< SMESH_subMesh* > & meshedSM, + SMESH_ProxyMesh::Ptr proxyMesh=SMESH_ProxyMesh::Ptr()); static void fixIntFaces(const netgen::OCCGeometry& occgeom, netgen::Mesh& ngMesh, @@ -134,6 +137,10 @@ class NETGENPLUGIN_EXPORT NETGENPlugin_Mesher static SMESH_ComputeErrorPtr readErrors(const std::vector< const SMDS_MeshNode* >& nodeVec); + + static void toPython( const netgen::Mesh* ngMesh, + const std::string& pyFile); // debug + private: SMESH_Mesh* _mesh; const TopoDS_Shape& _shape; -- 2.39.2