src/StdMeshers/StdMeshers_RadialPrism_3D.cxx

   1 //  SMESH SMESH : implementaion of SMESH idl descriptions
   2 //
   3 //  Copyright (C) 2003  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 //  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   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.
  10 //
  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.
  15 //
  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
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22 //
  23 //
  24 // File      : StdMeshers_RadialPrism_3D.cxx
  25 // Module    : SMESH
  26 // Created   : Fri Oct 20 11:37:07 2006
  27 // Author    : Edward AGAPOV (eap)
  28
  29
  30 #include "StdMeshers_RadialPrism_3D.hxx"
  31
  32 #include "StdMeshers_ProjectionUtils.hxx"
  33 #include "StdMeshers_NumberOfLayers.hxx"
  34 #include "StdMeshers_LayerDistribution.hxx"
  35 #include "StdMeshers_Prism_3D.hxx"
  36 #include "StdMeshers_Regular_1D.hxx"
  37
  38 #include "SMDS_MeshNode.hxx"
  39 #include "SMESHDS_SubMesh.hxx"
  40 #include "SMESH_Gen.hxx"
  41 #include "SMESH_Mesh.hxx"
  42 #include "SMESH_MeshEditor.hxx"
  43 #include "SMESH_MesherHelper.hxx"
  44 #include "SMESH_subMesh.hxx"
  45
  46 #include "utilities.h"
  47
  48 #include <TopoDS_Solid.hxx>
  49 #include <TopoDS_Shell.hxx>
  50 #include <BRepTools.hxx>
  51 #include <BRepAdaptor_Curve.hxx>
  52 #include <TopTools_ListIteratorOfListOfShape.hxx>
  53 #include <BRepBuilderAPI_MakeEdge.hxx>
  54 #include <gp.hxx>
  55 #include <gp_Pnt.hxx>
  56
  57
  58 using namespace std;
  59
  60 #define RETURN_BAD_RESULT(msg) { MESSAGE(msg); return false; }
  61 #define gpXYZ(n) gp_XYZ(n->X(),n->Y(),n->Z())
  62
  63 typedef StdMeshers_ProjectionUtils TAssocTool;
  64
  65 //=======================================================================
  66 //function : StdMeshers_RadialPrism_3D
  67 //purpose  :
  68 //=======================================================================
  69
  70 StdMeshers_RadialPrism_3D::StdMeshers_RadialPrism_3D(int hypId, int studyId, SMESH_Gen* gen)
  71   :SMESH_3D_Algo(hypId, studyId, gen)
  72 {
  73   _name = "RadialPrism_3D";
  74   _shapeType = (1 << TopAbs_SOLID);     // 1 bit per shape type
  75
  76   _compatibleHypothesis.push_back("LayerDistribution");
  77   _compatibleHypothesis.push_back("NumberOfLayers");
  78   myNbLayerHypo = 0;
  79   myDistributionHypo = 0;
  80 }
  81
  82 //================================================================================
  83 /*!
  84  * \brief Destructor
  85  */
  86 //================================================================================
  87
  88 StdMeshers_RadialPrism_3D::~StdMeshers_RadialPrism_3D()
  89 {}
  90
  91 //=======================================================================
  92 //function : CheckHypothesis
  93 //purpose  :
  94 //=======================================================================
  95
  96 bool StdMeshers_RadialPrism_3D::CheckHypothesis(SMESH_Mesh&                          aMesh,
  97                                                 const TopoDS_Shape&                  aShape,
  98                                                 SMESH_Hypothesis::Hypothesis_Status& aStatus)
  99 {
 100   // check aShape that must have 2 shells
 101   if ( TAssocTool::Count( aShape, TopAbs_SOLID, 0 ) != 1 ||
 102        TAssocTool::Count( aShape, TopAbs_SHELL, 0 ) != 2 )
 103   {
 104     aStatus = HYP_BAD_GEOMETRY;
 105     return false;
 106   }
 107
 108   myNbLayerHypo = 0;
 109   myDistributionHypo = 0;
 110
 111   list <const SMESHDS_Hypothesis * >::const_iterator itl;
 112
 113   const list <const SMESHDS_Hypothesis * >&hyps = GetUsedHypothesis(aMesh, aShape);
 114   if ( hyps.size() == 0 )
 115   {
 116     aStatus = SMESH_Hypothesis::HYP_MISSING;
 117     return false;  // can't work with no hypothesis
 118   }
 119
 120   if ( hyps.size() > 1 )
 121   {
 122     aStatus = SMESH_Hypothesis::HYP_ALREADY_EXIST;
 123     return false;
 124   }
 125
 126   const SMESHDS_Hypothesis *theHyp = hyps.front();
 127
 128   string hypName = theHyp->GetName();
 129
 130   if (hypName == "NumberOfLayers")
 131   {
 132     myNbLayerHypo = static_cast<const StdMeshers_NumberOfLayers *>(theHyp);
 133     aStatus = SMESH_Hypothesis::HYP_OK;
 134     return true;
 135   }
 136   if (hypName == "LayerDistribution")
 137   {
 138     myDistributionHypo = static_cast<const StdMeshers_LayerDistribution *>(theHyp);
 139     aStatus = SMESH_Hypothesis::HYP_OK;
 140     return true;
 141   }
 142   aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
 143   return true;
 144 }
 145
 146 //=======================================================================
 147 //function : Compute
 148 //purpose  :
 149 //=======================================================================
 150
 151 bool StdMeshers_RadialPrism_3D::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape)
 152 {
 153   TopExp_Explorer exp;
 154   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 155
 156   myHelper = new SMESH_MesherHelper( aMesh );
 157   myHelper->IsQuadraticSubMesh( aShape );
 158   // to delete helper at exit from Compute()
 159   std::auto_ptr<SMESH_MesherHelper> helperDeleter( myHelper );
 160
 161   // get 2 shells
 162   TopoDS_Solid solid = TopoDS::Solid( aShape );
 163   TopoDS_Shell outerShell = BRepTools::OuterShell( solid );
 164   TopoDS_Shape innerShell;
 165   int nbShells = 0;
 166   for ( TopoDS_Iterator It (solid); It.More(); It.Next(), ++nbShells )
 167     if ( !outerShell.IsSame( It.Value() ))
 168       innerShell = It.Value();
 169   if ( nbShells != 2 )
 170     RETURN_BAD_RESULT("Must be 2 shells");
 171
 172   // ----------------------------------
 173   // Associate subshapes of the shells
 174   // ----------------------------------
 175
 176   TAssocTool::TShapeShapeMap shape2ShapeMap;
 177   if ( !TAssocTool::FindSubShapeAssociation( outerShell, &aMesh,
 178                                              innerShell, &aMesh,
 179                                              shape2ShapeMap) )
 180     RETURN_BAD_RESULT("FindSubShapeAssociation failed");
 181
 182   // ------------------
 183   // Make mesh
 184   // ------------------
 185
 186   TNode2ColumnMap node2columnMap;
 187   myLayerPositions.clear();
 188
 189   for ( exp.Init( outerShell, TopAbs_FACE ); exp.More(); exp.Next() )
 190   {
 191     // Corresponding subshapes
 192     TopoDS_Face outFace = TopoDS::Face( exp.Current() );
 193     TopoDS_Face inFace;
 194     if ( !shape2ShapeMap.IsBound( outFace )) {
 195       RETURN_BAD_RESULT("Association not found for face " << meshDS->ShapeToIndex( outFace ));
 196     } else {
 197       inFace = TopoDS::Face( shape2ShapeMap( outFace ));
 198     }
 199
 200     // Find matching nodes of in and out faces
 201     TNodeNodeMap nodeIn2OutMap;
 202     if ( ! TAssocTool::FindMatchingNodesOnFaces( inFace, &aMesh, outFace, &aMesh,
 203                                                  shape2ShapeMap, nodeIn2OutMap ))
 204       RETURN_BAD_RESULT("Different mesh on corresponding out and in faces: "
 205                         << meshDS->ShapeToIndex( outFace ) << " and "
 206                         << meshDS->ShapeToIndex( inFace ));
 207     // Create volumes
 208
 209     SMDS_ElemIteratorPtr faceIt = meshDS->MeshElements( inFace )->GetElements();
 210     while ( faceIt->more() ) // loop on faces on inFace
 211     {
 212       const SMDS_MeshElement* face = faceIt->next();
 213       if ( !face || face->GetType() != SMDSAbs_Face )
 214         continue;
 215       int nbNodes = face->NbNodes();
 216       if ( face->IsQuadratic() )
 217         nbNodes /= 2;
 218
 219       // find node columns for each node
 220       vector< const TNodeColumn* > columns( nbNodes );
 221       for ( int i = 0; i < nbNodes; ++i )
 222       {
 223         const SMDS_MeshNode* n = face->GetNode( i );
 224         TNode2ColumnMap::iterator n_col = node2columnMap.find( n );
 225         if ( n_col != node2columnMap.end() )
 226           columns[ i ] = & n_col->second;
 227         else
 228           columns[ i ] = makeNodeColumn( node2columnMap, n, nodeIn2OutMap[ n ] );
 229       }
 230
 231       StdMeshers_Prism_3D::AddPrisms( columns, myHelper );
 232     }
 233   } // loop on faces of out shell
 234
 235   return true;
 236 }
 237
 238 //================================================================================
 239 /*!
 240  * \brief Create a column of nodes from outNode to inNode
 241   * \param n2ColMap - map of node columns to add a created column
 242   * \param outNode - botton node of a column
 243   * \param inNode - top node of a column
 244   * \retval const TNodeColumn* - a new column pointer
 245  */
 246 //================================================================================
 247
 248 TNodeColumn* StdMeshers_RadialPrism_3D::makeNodeColumn( TNode2ColumnMap&     n2ColMap,
 249                                                         const SMDS_MeshNode* outNode,
 250                                                         const SMDS_MeshNode* inNode)
 251 {
 252   SMESHDS_Mesh * meshDS = myHelper->GetMeshDS();
 253   int shapeID = myHelper->GetSubShapeID();
 254
 255   if ( myLayerPositions.empty() )
 256     computeLayerPositions( gpXYZ( inNode ), gpXYZ( outNode ));
 257
 258   int nbSegments = myLayerPositions.size() + 1;
 259
 260   TNode2ColumnMap::iterator n_col =
 261     n2ColMap.insert( make_pair( outNode, TNodeColumn() )).first;
 262   TNodeColumn & column = n_col->second;
 263   column.resize( nbSegments + 1 );
 264   column.front() = outNode;
 265   column.back() =  inNode;
 266
 267   gp_XYZ p1 = gpXYZ( outNode );
 268   gp_XYZ p2 = gpXYZ( inNode );
 269   for ( int z = 1; z < nbSegments; ++z )
 270   {
 271     double r = myLayerPositions[ z - 1 ];
 272     gp_XYZ p = ( 1 - r ) * p1 + r * p2;
 273     SMDS_MeshNode* n = meshDS->AddNode( p.X(), p.Y(), p.Z() );
 274     meshDS->SetNodeInVolume( n, shapeID );
 275     column[ z ] = n;
 276   }
 277
 278   return & column;
 279 }
 280
 281 //================================================================================
 282 //================================================================================
 283 /*!
 284  * \brief Class computing layers distribution using data of
 285  *        StdMeshers_LayerDistribution hypothesis
 286  */
 287 //================================================================================
 288 //================================================================================
 289
 290 class TNodeDistributor: private StdMeshers_Regular_1D
 291 {
 292   list <const SMESHDS_Hypothesis *> myUsedHyps;
 293 public:
 294   // -----------------------------------------------------------------------------
 295   static TNodeDistributor* GetDistributor(SMESH_Mesh& aMesh)
 296   {
 297     const int myID = -1000;
 298     map < int, SMESH_1D_Algo * > & algoMap = aMesh.GetGen()->_map1D_Algo;
 299     map < int, SMESH_1D_Algo * >::iterator id_algo = algoMap.find( myID );
 300     if ( id_algo == algoMap.end() )
 301       return new TNodeDistributor( myID, 0, aMesh.GetGen() );
 302     return static_cast< TNodeDistributor* >( id_algo->second );
 303   }
 304   // -----------------------------------------------------------------------------
 305   bool Compute( vector< double > &                  positions,
 306                 gp_Pnt                              pIn,
 307                 gp_Pnt                              pOut,
 308                 SMESH_Mesh&                         aMesh,
 309                 const StdMeshers_LayerDistribution* hyp)
 310   {
 311     double len = pIn.Distance( pOut );
 312     if ( len <= DBL_MIN ) RETURN_BAD_RESULT("Bad points");
 313
 314     if ( !hyp || !hyp->GetLayerDistribution() )
 315       RETURN_BAD_RESULT("Bad StdMeshers_LayerDistribution hypothesis");
 316     myUsedHyps.clear();
 317     myUsedHyps.push_back( hyp->GetLayerDistribution() );
 318
 319     TopoDS_Edge edge = BRepBuilderAPI_MakeEdge( pIn, pOut );
 320     SMESH_Hypothesis::Hypothesis_Status aStatus;
 321     if ( !StdMeshers_Regular_1D::CheckHypothesis( aMesh, edge, aStatus ))
 322       RETURN_BAD_RESULT("StdMeshers_Regular_1D::CheckHypothesis() failed with status "<<aStatus);
 323
 324     BRepAdaptor_Curve C3D(edge);
 325     double f = C3D.FirstParameter(), l = C3D.LastParameter();
 326     list< double > params;
 327     if ( !StdMeshers_Regular_1D::computeInternalParameters( C3D, len, f, l, params, false ))
 328       RETURN_BAD_RESULT("StdMeshers_Regular_1D::computeInternalParameters() failed");
 329
 330     positions.clear();
 331     positions.reserve( params.size() );
 332     for (list<double>::iterator itU = params.begin(); itU != params.end(); itU++)
 333       positions.push_back( *itU / len );
 334     return true;
 335   }
 336 protected:
 337   // -----------------------------------------------------------------------------
 338   TNodeDistributor( int hypId, int studyId, SMESH_Gen* gen)
 339     : StdMeshers_Regular_1D( hypId, studyId, gen)
 340   {
 341   }
 342   // -----------------------------------------------------------------------------
 343   virtual const list <const SMESHDS_Hypothesis *> &
 344     GetUsedHypothesis(SMESH_Mesh &, const TopoDS_Shape &, const bool)
 345   {
 346     return myUsedHyps;
 347   }
 348   // -----------------------------------------------------------------------------
 349 };
 350
 351 //================================================================================
 352 /*!
 353  * \brief Compute positions of nodes between the internal and the external surfaces
 354   * \retval bool - is a success
 355  */
 356 //================================================================================
 357
 358 bool StdMeshers_RadialPrism_3D::computeLayerPositions(gp_Pnt pIn, gp_Pnt pOut)
 359 {
 360   if ( myNbLayerHypo )
 361   {
 362     int nbSegments = myNbLayerHypo->GetNumberOfLayers();
 363     myLayerPositions.resize( nbSegments - 1 );
 364     for ( int z = 1; z < nbSegments; ++z )
 365       myLayerPositions[ z - 1 ] = double( z )/ double( nbSegments );
 366     return true;
 367   }
 368   if ( myDistributionHypo ) {
 369     SMESH_Mesh * mesh = myHelper->GetMesh();
 370     return TNodeDistributor::GetDistributor(*mesh)->Compute( myLayerPositions, pIn, pOut,
 371                                                              *mesh, myDistributionHypo );
 372   }
 373   RETURN_BAD_RESULT("Bad hypothesis");
 374 }