src/StdMeshers/StdMeshers_RadialPrism_3D.cxx

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