1 // Copyright (C) 2007-2011 CEA/DEN, EDF R&D, 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.
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
23 // SMESH SMESH : implementaion of SMESH idl descriptions
24 // File : StdMeshers_RadialPrism_3D.cxx
26 // Created : Fri Oct 20 11:37:07 2006
27 // Author : Edward AGAPOV (eap)
29 #include "StdMeshers_RadialPrism_3D.hxx"
31 #include "StdMeshers_ProjectionUtils.hxx"
32 #include "StdMeshers_NumberOfLayers.hxx"
33 #include "StdMeshers_LayerDistribution.hxx"
34 #include "StdMeshers_Prism_3D.hxx"
35 #include "StdMeshers_Regular_1D.hxx"
37 #include "SMDS_MeshNode.hxx"
38 #include "SMESHDS_SubMesh.hxx"
39 #include "SMESH_Gen.hxx"
40 #include "SMESH_Mesh.hxx"
41 #include "SMESH_MesherHelper.hxx"
42 #include "SMESH_subMesh.hxx"
44 #include "utilities.h"
46 #include <BRepAdaptor_Curve.hxx>
47 #include <BRepBuilderAPI_MakeEdge.hxx>
48 #include <BRepTools.hxx>
49 #include <BRep_Tool.hxx>
50 #include <TopExp_Explorer.hxx>
52 #include <TopoDS_Shell.hxx>
53 #include <TopoDS_Solid.hxx>
54 #include <TopTools_MapOfShape.hxx>
61 #define RETURN_BAD_RESULT(msg) { MESSAGE(")-: Error: " << msg); return false; }
62 #define gpXYZ(n) gp_XYZ(n->X(),n->Y(),n->Z())
64 typedef StdMeshers_ProjectionUtils TAssocTool;
66 //=======================================================================
67 //function : StdMeshers_RadialPrism_3D
69 //=======================================================================
71 StdMeshers_RadialPrism_3D::StdMeshers_RadialPrism_3D(int hypId, int studyId, SMESH_Gen* gen)
72 :SMESH_3D_Algo(hypId, studyId, gen)
74 _name = "RadialPrism_3D";
75 _shapeType = (1 << TopAbs_SOLID); // 1 bit per shape type
77 _compatibleHypothesis.push_back("LayerDistribution");
78 _compatibleHypothesis.push_back("NumberOfLayers");
80 myDistributionHypo = 0;
83 //================================================================================
87 //================================================================================
89 StdMeshers_RadialPrism_3D::~StdMeshers_RadialPrism_3D()
92 //=======================================================================
93 //function : CheckHypothesis
95 //=======================================================================
97 bool StdMeshers_RadialPrism_3D::CheckHypothesis(SMESH_Mesh& aMesh,
98 const TopoDS_Shape& aShape,
99 SMESH_Hypothesis::Hypothesis_Status& aStatus)
101 // check aShape that must have 2 shells
103 if ( TAssocTool::Count( aShape, TopAbs_SOLID, 0 ) != 1 ||
104 TAssocTool::Count( aShape, TopAbs_SHELL, 0 ) != 2 )
106 aStatus = HYP_BAD_GEOMETRY;
111 myDistributionHypo = 0;
113 list <const SMESHDS_Hypothesis * >::const_iterator itl;
115 const list <const SMESHDS_Hypothesis * >&hyps = GetUsedHypothesis(aMesh, aShape);
116 if ( hyps.size() == 0 )
118 aStatus = SMESH_Hypothesis::HYP_MISSING;
119 return false; // can't work with no hypothesis
122 if ( hyps.size() > 1 )
124 aStatus = SMESH_Hypothesis::HYP_ALREADY_EXIST;
128 const SMESHDS_Hypothesis *theHyp = hyps.front();
130 string hypName = theHyp->GetName();
132 if (hypName == "NumberOfLayers")
134 myNbLayerHypo = static_cast<const StdMeshers_NumberOfLayers *>(theHyp);
135 aStatus = SMESH_Hypothesis::HYP_OK;
138 if (hypName == "LayerDistribution")
140 myDistributionHypo = static_cast<const StdMeshers_LayerDistribution *>(theHyp);
141 aStatus = SMESH_Hypothesis::HYP_OK;
144 aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
148 //=======================================================================
151 //=======================================================================
153 bool StdMeshers_RadialPrism_3D::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape)
156 SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
158 myHelper = new SMESH_MesherHelper( aMesh );
159 myHelper->IsQuadraticSubMesh( aShape );
160 // to delete helper at exit from Compute()
161 std::auto_ptr<SMESH_MesherHelper> helperDeleter( myHelper );
164 TopoDS_Solid solid = TopoDS::Solid( aShape );
165 TopoDS_Shell outerShell = BRepTools::OuterShell( solid );
166 TopoDS_Shape innerShell;
168 for ( TopoDS_Iterator It (solid); It.More(); It.Next(), ++nbShells )
169 if ( !outerShell.IsSame( It.Value() ))
170 innerShell = It.Value();
172 return error(COMPERR_BAD_SHAPE, SMESH_Comment("Must be 2 shells but not ")<<nbShells);
174 // ----------------------------------
175 // Associate subshapes of the shells
176 // ----------------------------------
178 TAssocTool::TShapeShapeMap shape2ShapeMap;
179 if ( !TAssocTool::FindSubShapeAssociation( outerShell, &aMesh,
182 return error(COMPERR_BAD_SHAPE,"Topology of inner and outer shells seems different" );
184 // ------------------
186 // ------------------
188 TNode2ColumnMap node2columnMap;
189 myLayerPositions.clear();
191 for ( exp.Init( outerShell, TopAbs_FACE ); exp.More(); exp.Next() )
193 // Corresponding subshapes
194 TopoDS_Face outFace = TopoDS::Face( exp.Current() );
196 if ( !shape2ShapeMap.IsBound( outFace )) {
197 return error(SMESH_Comment("Corresponding inner face not found for face #" )
198 << meshDS->ShapeToIndex( outFace ));
200 inFace = TopoDS::Face( shape2ShapeMap( outFace ));
203 // Find matching nodes of in and out faces
204 TNodeNodeMap nodeIn2OutMap;
205 if ( ! TAssocTool::FindMatchingNodesOnFaces( inFace, &aMesh, outFace, &aMesh,
206 shape2ShapeMap, nodeIn2OutMap ))
207 return error(COMPERR_BAD_INPUT_MESH,SMESH_Comment("Mesh on faces #")
208 << meshDS->ShapeToIndex( outFace ) << " and "
209 << meshDS->ShapeToIndex( inFace ) << " seems different" );
213 SMDS_ElemIteratorPtr faceIt = meshDS->MeshElements( inFace )->GetElements();
214 while ( faceIt->more() ) // loop on faces on inFace
216 const SMDS_MeshElement* face = faceIt->next();
217 if ( !face || face->GetType() != SMDSAbs_Face )
219 int nbNodes = face->NbNodes();
220 if ( face->IsQuadratic() )
223 // find node columns for each node
224 vector< const TNodeColumn* > columns( nbNodes );
225 for ( int i = 0; i < nbNodes; ++i )
227 const SMDS_MeshNode* nIn = face->GetNode( i );
228 TNode2ColumnMap::iterator n_col = node2columnMap.find( nIn );
229 if ( n_col != node2columnMap.end() ) {
230 columns[ i ] = & n_col->second;
233 TNodeNodeMap::iterator nInOut = nodeIn2OutMap.find( nIn );
234 if ( nInOut == nodeIn2OutMap.end() )
235 RETURN_BAD_RESULT("No matching node for "<< nIn->GetID() <<
236 " in face "<< face->GetID());
237 columns[ i ] = makeNodeColumn( node2columnMap, nIn, nInOut->second );
241 StdMeshers_Prism_3D::AddPrisms( columns, myHelper );
243 } // loop on faces of out shell
248 //================================================================================
250 * \brief Create a column of nodes from outNode to inNode
251 * \param n2ColMap - map of node columns to add a created column
252 * \param outNode - botton node of a column
253 * \param inNode - top node of a column
254 * \retval const TNodeColumn* - a new column pointer
256 //================================================================================
258 TNodeColumn* StdMeshers_RadialPrism_3D::makeNodeColumn( TNode2ColumnMap& n2ColMap,
259 const SMDS_MeshNode* outNode,
260 const SMDS_MeshNode* inNode)
262 SMESHDS_Mesh * meshDS = myHelper->GetMeshDS();
263 int shapeID = myHelper->GetSubShapeID();
265 if ( myLayerPositions.empty() ) {
266 gp_Pnt pIn = gpXYZ( inNode ), pOut = gpXYZ( outNode );
267 computeLayerPositions( pIn, pOut );
269 int nbSegments = myLayerPositions.size() + 1;
271 TNode2ColumnMap::iterator n_col =
272 n2ColMap.insert( make_pair( outNode, TNodeColumn() )).first;
273 TNodeColumn & column = n_col->second;
274 column.resize( nbSegments + 1 );
275 column.front() = outNode;
276 column.back() = inNode;
278 gp_XYZ p1 = gpXYZ( outNode );
279 gp_XYZ p2 = gpXYZ( inNode );
280 for ( int z = 1; z < nbSegments; ++z )
282 double r = myLayerPositions[ z - 1 ];
283 gp_XYZ p = ( 1 - r ) * p1 + r * p2;
284 SMDS_MeshNode* n = meshDS->AddNode( p.X(), p.Y(), p.Z() );
285 meshDS->SetNodeInVolume( n, shapeID );
292 //================================================================================
293 //================================================================================
295 * \brief Class computing layers distribution using data of
296 * StdMeshers_LayerDistribution hypothesis
298 //================================================================================
299 //================================================================================
301 class TNodeDistributor: public StdMeshers_Regular_1D
303 list <const SMESHDS_Hypothesis *> myUsedHyps;
305 // -----------------------------------------------------------------------------
306 static TNodeDistributor* GetDistributor(SMESH_Mesh& aMesh)
308 const int myID = -1000;
309 map < int, SMESH_1D_Algo * > & algoMap = aMesh.GetGen()->_map1D_Algo;
310 map < int, SMESH_1D_Algo * >::iterator id_algo = algoMap.find( myID );
311 if ( id_algo == algoMap.end() )
312 return new TNodeDistributor( myID, 0, aMesh.GetGen() );
313 return static_cast< TNodeDistributor* >( id_algo->second );
315 // -----------------------------------------------------------------------------
316 bool Compute( vector< double > & positions,
320 const StdMeshers_LayerDistribution* hyp)
322 double len = pIn.Distance( pOut );
323 if ( len <= DBL_MIN ) return error("Too close points of inner and outer shells");
325 if ( !hyp || !hyp->GetLayerDistribution() )
326 return error( "Invalid LayerDistribution hypothesis");
328 myUsedHyps.push_back( hyp->GetLayerDistribution() );
330 TopoDS_Edge edge = BRepBuilderAPI_MakeEdge( pIn, pOut );
331 SMESH_Hypothesis::Hypothesis_Status aStatus;
332 if ( !StdMeshers_Regular_1D::CheckHypothesis( aMesh, edge, aStatus ))
333 return error( "StdMeshers_Regular_1D::CheckHypothesis() failed "
334 "with LayerDistribution hypothesis");
336 BRepAdaptor_Curve C3D(edge);
337 double f = C3D.FirstParameter(), l = C3D.LastParameter();
338 list< double > params;
339 if ( !StdMeshers_Regular_1D::computeInternalParameters( aMesh, C3D, len, f, l, params, false ))
340 return error("StdMeshers_Regular_1D failed to compute layers distribution");
343 positions.reserve( params.size() );
344 for (list<double>::iterator itU = params.begin(); itU != params.end(); itU++)
345 positions.push_back( *itU / len );
349 // -----------------------------------------------------------------------------
350 TNodeDistributor( int hypId, int studyId, SMESH_Gen* gen)
351 : StdMeshers_Regular_1D( hypId, studyId, gen)
354 // -----------------------------------------------------------------------------
355 virtual const list <const SMESHDS_Hypothesis *> &
356 GetUsedHypothesis(SMESH_Mesh &, const TopoDS_Shape &, const bool)
360 // -----------------------------------------------------------------------------
363 //================================================================================
365 * \brief Compute positions of nodes between the internal and the external surfaces
366 * \retval bool - is a success
368 //================================================================================
370 bool StdMeshers_RadialPrism_3D::computeLayerPositions(const gp_Pnt& pIn,
375 int nbSegments = myNbLayerHypo->GetNumberOfLayers();
376 myLayerPositions.resize( nbSegments - 1 );
377 for ( int z = 1; z < nbSegments; ++z )
378 myLayerPositions[ z - 1 ] = double( z )/ double( nbSegments );
381 if ( myDistributionHypo ) {
382 SMESH_Mesh * mesh = myHelper->GetMesh();
383 if ( !TNodeDistributor::GetDistributor(*mesh)->Compute( myLayerPositions, pIn, pOut,
384 *mesh, myDistributionHypo ))
386 error( TNodeDistributor::GetDistributor(*mesh)->GetComputeError() );
390 RETURN_BAD_RESULT("Bad hypothesis");
394 //=======================================================================
395 //function : Evaluate
397 //=======================================================================
399 bool StdMeshers_RadialPrism_3D::Evaluate(SMESH_Mesh& aMesh,
400 const TopoDS_Shape& aShape,
401 MapShapeNbElems& aResMap)
404 TopoDS_Solid solid = TopoDS::Solid( aShape );
405 TopoDS_Shell outerShell = BRepTools::OuterShell( solid );
406 TopoDS_Shape innerShell;
408 for ( TopoDS_Iterator It (solid); It.More(); It.Next(), ++nbShells )
409 if ( !outerShell.IsSame( It.Value() ))
410 innerShell = It.Value();
411 if ( nbShells != 2 ) {
412 std::vector<int> aResVec(SMDSEntity_Last);
413 for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
414 SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
415 aResMap.insert(std::make_pair(sm,aResVec));
416 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
417 smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
421 // Associate subshapes of the shells
422 TAssocTool::TShapeShapeMap shape2ShapeMap;
423 if ( !TAssocTool::FindSubShapeAssociation( outerShell, &aMesh,
426 std::vector<int> aResVec(SMDSEntity_Last);
427 for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
428 SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
429 aResMap.insert(std::make_pair(sm,aResVec));
430 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
431 smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
435 // get info for outer shell
436 int nb0d_Out=0, nb2d_3_Out=0, nb2d_4_Out=0;
437 //TopTools_SequenceOfShape FacesOut;
438 for (TopExp_Explorer exp(outerShell, TopAbs_FACE); exp.More(); exp.Next()) {
439 //FacesOut.Append(exp.Current());
440 SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
441 MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
442 std::vector<int> aVec = (*anIt).second;
443 nb0d_Out += aVec[SMDSEntity_Node];
444 nb2d_3_Out += Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
445 nb2d_4_Out += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
448 TopTools_MapOfShape tmpMap;
449 for (TopExp_Explorer exp(outerShell, TopAbs_EDGE); exp.More(); exp.Next()) {
450 if( tmpMap.Contains( exp.Current() ) )
452 tmpMap.Add( exp.Current() );
453 SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
454 MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
455 std::vector<int> aVec = (*anIt).second;
456 nb0d_Out += aVec[SMDSEntity_Node];
457 nb1d_Out += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
460 for (TopExp_Explorer exp(outerShell, TopAbs_VERTEX); exp.More(); exp.Next()) {
461 if( tmpMap.Contains( exp.Current() ) )
463 tmpMap.Add( exp.Current() );
467 // get info for inner shell
468 int nb0d_In=0, nb2d_3_In=0, nb2d_4_In=0;
469 //TopTools_SequenceOfShape FacesIn;
470 for (TopExp_Explorer exp(innerShell, TopAbs_FACE); exp.More(); exp.Next()) {
471 //FacesIn.Append(exp.Current());
472 SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
473 MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
474 std::vector<int> aVec = (*anIt).second;
475 nb0d_In += aVec[SMDSEntity_Node];
476 nb2d_3_In += Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
477 nb2d_4_In += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
481 bool IsQuadratic = false;
483 for (TopExp_Explorer exp(innerShell, TopAbs_EDGE); exp.More(); exp.Next()) {
484 if( tmpMap.Contains( exp.Current() ) )
486 tmpMap.Add( exp.Current() );
487 SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
488 MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
489 std::vector<int> aVec = (*anIt).second;
490 nb0d_In += aVec[SMDSEntity_Node];
491 nb1d_In += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
493 IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
498 for (TopExp_Explorer exp(innerShell, TopAbs_VERTEX); exp.More(); exp.Next()) {
499 if( tmpMap.Contains( exp.Current() ) )
501 tmpMap.Add( exp.Current() );
505 bool IsOK = (nb0d_Out==nb0d_In) && (nb1d_Out==nb1d_In) &&
506 (nb2d_3_Out==nb2d_3_In) && (nb2d_4_Out==nb2d_4_In);
508 std::vector<int> aResVec(SMDSEntity_Last);
509 for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
510 SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
511 aResMap.insert(std::make_pair(sm,aResVec));
512 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
513 smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
518 if( myNbLayerHypo ) {
519 nbLayers = myNbLayerHypo->GetNumberOfLayers();
521 if ( myDistributionHypo ) {
522 if ( !myDistributionHypo->GetLayerDistribution() ) {
523 std::vector<int> aResVec(SMDSEntity_Last);
524 for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
525 SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
526 aResMap.insert(std::make_pair(sm,aResVec));
527 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
528 smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
531 TopExp_Explorer exp(outerShell, TopAbs_VERTEX);
532 TopoDS_Vertex Vout = TopoDS::Vertex(exp.Current());
533 TopoDS_Vertex Vin = TopoDS::Vertex( shape2ShapeMap(Vout) );
534 if ( myLayerPositions.empty() ) {
535 gp_Pnt pIn = BRep_Tool::Pnt(Vin);
536 gp_Pnt pOut = BRep_Tool::Pnt(Vout);
537 computeLayerPositions( pIn, pOut );
539 nbLayers = myLayerPositions.size() + 1;
542 std::vector<int> aResVec(SMDSEntity_Last);
543 for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
545 aResVec[SMDSEntity_Quad_Penta] = nb2d_3_Out * nbLayers;
546 aResVec[SMDSEntity_Quad_Hexa] = nb2d_4_Out * nbLayers;
547 int nb1d = ( nb2d_3_Out*3 + nb2d_4_Out*4 ) / 2;
548 aResVec[SMDSEntity_Node] = nb0d_Out * ( 2*nbLayers - 1 ) - nb1d * nbLayers;
551 aResVec[SMDSEntity_Node] = nb0d_Out * ( nbLayers - 1 );
552 aResVec[SMDSEntity_Penta] = nb2d_3_Out * nbLayers;
553 aResVec[SMDSEntity_Hexa] = nb2d_4_Out * nbLayers;
555 SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
556 aResMap.insert(std::make_pair(sm,aResVec));