00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include "StdMeshers_RadialPrism_3D.hxx"
00030
00031 #include "StdMeshers_ProjectionUtils.hxx"
00032 #include "StdMeshers_NumberOfLayers.hxx"
00033 #include "StdMeshers_LayerDistribution.hxx"
00034 #include "StdMeshers_Prism_3D.hxx"
00035 #include "StdMeshers_Regular_1D.hxx"
00036
00037 #include "SMDS_MeshNode.hxx"
00038 #include "SMESHDS_SubMesh.hxx"
00039 #include "SMESH_Gen.hxx"
00040 #include "SMESH_Mesh.hxx"
00041 #include "SMESH_MesherHelper.hxx"
00042 #include "SMESH_subMesh.hxx"
00043
00044 #include "utilities.h"
00045
00046 #include <BRepAdaptor_Curve.hxx>
00047 #include <BRepBuilderAPI_MakeEdge.hxx>
00048 #include <BRepTools.hxx>
00049 #include <BRep_Tool.hxx>
00050 #include <TopExp_Explorer.hxx>
00051 #include <TopoDS.hxx>
00052 #include <TopoDS_Shell.hxx>
00053 #include <TopoDS_Solid.hxx>
00054 #include <TopTools_MapOfShape.hxx>
00055 #include <gp.hxx>
00056 #include <gp_Pnt.hxx>
00057
00058
00059 using namespace std;
00060
00061 #define RETURN_BAD_RESULT(msg) { MESSAGE(")-: Error: " << msg); return false; }
00062 #define gpXYZ(n) gp_XYZ(n->X(),n->Y(),n->Z())
00063
00064 typedef StdMeshers_ProjectionUtils TAssocTool;
00065
00066
00067
00068
00069
00070
00071 StdMeshers_RadialPrism_3D::StdMeshers_RadialPrism_3D(int hypId, int studyId, SMESH_Gen* gen)
00072 :SMESH_3D_Algo(hypId, studyId, gen)
00073 {
00074 _name = "RadialPrism_3D";
00075 _shapeType = (1 << TopAbs_SOLID);
00076
00077 _compatibleHypothesis.push_back("LayerDistribution");
00078 _compatibleHypothesis.push_back("NumberOfLayers");
00079 myNbLayerHypo = 0;
00080 myDistributionHypo = 0;
00081 }
00082
00083
00087
00088
00089 StdMeshers_RadialPrism_3D::~StdMeshers_RadialPrism_3D()
00090 {}
00091
00092
00093
00094
00095
00096
00097 bool StdMeshers_RadialPrism_3D::CheckHypothesis(SMESH_Mesh& aMesh,
00098 const TopoDS_Shape& aShape,
00099 SMESH_Hypothesis::Hypothesis_Status& aStatus)
00100 {
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110 myNbLayerHypo = 0;
00111 myDistributionHypo = 0;
00112
00113 list <const SMESHDS_Hypothesis * >::const_iterator itl;
00114
00115 const list <const SMESHDS_Hypothesis * >&hyps = GetUsedHypothesis(aMesh, aShape);
00116 if ( hyps.size() == 0 )
00117 {
00118 aStatus = SMESH_Hypothesis::HYP_MISSING;
00119 return false;
00120 }
00121
00122 if ( hyps.size() > 1 )
00123 {
00124 aStatus = SMESH_Hypothesis::HYP_ALREADY_EXIST;
00125 return false;
00126 }
00127
00128 const SMESHDS_Hypothesis *theHyp = hyps.front();
00129
00130 string hypName = theHyp->GetName();
00131
00132 if (hypName == "NumberOfLayers")
00133 {
00134 myNbLayerHypo = static_cast<const StdMeshers_NumberOfLayers *>(theHyp);
00135 aStatus = SMESH_Hypothesis::HYP_OK;
00136 return true;
00137 }
00138 if (hypName == "LayerDistribution")
00139 {
00140 myDistributionHypo = static_cast<const StdMeshers_LayerDistribution *>(theHyp);
00141 aStatus = SMESH_Hypothesis::HYP_OK;
00142 return true;
00143 }
00144 aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
00145 return true;
00146 }
00147
00148
00149
00150
00151
00152
00153 bool StdMeshers_RadialPrism_3D::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aShape)
00154 {
00155 TopExp_Explorer exp;
00156 SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
00157
00158 myHelper = new SMESH_MesherHelper( aMesh );
00159 myHelper->IsQuadraticSubMesh( aShape );
00160
00161 std::auto_ptr<SMESH_MesherHelper> helperDeleter( myHelper );
00162
00163
00164 TopoDS_Solid solid = TopoDS::Solid( aShape );
00165 TopoDS_Shell outerShell = BRepTools::OuterShell( solid );
00166 TopoDS_Shape innerShell;
00167 int nbShells = 0;
00168 for ( TopoDS_Iterator It (solid); It.More(); It.Next(), ++nbShells )
00169 if ( !outerShell.IsSame( It.Value() ))
00170 innerShell = It.Value();
00171 if ( nbShells != 2 )
00172 return error(COMPERR_BAD_SHAPE, SMESH_Comment("Must be 2 shells but not ")<<nbShells);
00173
00174
00175
00176
00177
00178 TAssocTool::TShapeShapeMap shape2ShapeMap;
00179 if ( !TAssocTool::FindSubShapeAssociation( outerShell, &aMesh,
00180 innerShell, &aMesh,
00181 shape2ShapeMap) )
00182 return error(COMPERR_BAD_SHAPE,"Topology of inner and outer shells seems different" );
00183
00184
00185
00186
00187
00188 TNode2ColumnMap node2columnMap;
00189 myLayerPositions.clear();
00190
00191 for ( exp.Init( outerShell, TopAbs_FACE ); exp.More(); exp.Next() )
00192 {
00193
00194 TopoDS_Face outFace = TopoDS::Face( exp.Current() );
00195 TopoDS_Face inFace;
00196 if ( !shape2ShapeMap.IsBound( outFace )) {
00197 return error(SMESH_Comment("Corresponding inner face not found for face #" )
00198 << meshDS->ShapeToIndex( outFace ));
00199 } else {
00200 inFace = TopoDS::Face( shape2ShapeMap( outFace ));
00201 }
00202
00203
00204 TNodeNodeMap nodeIn2OutMap;
00205 if ( ! TAssocTool::FindMatchingNodesOnFaces( inFace, &aMesh, outFace, &aMesh,
00206 shape2ShapeMap, nodeIn2OutMap ))
00207 return error(COMPERR_BAD_INPUT_MESH,SMESH_Comment("Mesh on faces #")
00208 << meshDS->ShapeToIndex( outFace ) << " and "
00209 << meshDS->ShapeToIndex( inFace ) << " seems different" );
00210
00211
00212
00213 SMDS_ElemIteratorPtr faceIt = meshDS->MeshElements( inFace )->GetElements();
00214 while ( faceIt->more() )
00215 {
00216 const SMDS_MeshElement* face = faceIt->next();
00217 if ( !face || face->GetType() != SMDSAbs_Face )
00218 continue;
00219 int nbNodes = face->NbNodes();
00220 if ( face->IsQuadratic() )
00221 nbNodes /= 2;
00222
00223
00224 vector< const TNodeColumn* > columns( nbNodes );
00225 for ( int i = 0; i < nbNodes; ++i )
00226 {
00227 const SMDS_MeshNode* nIn = face->GetNode( i );
00228 TNode2ColumnMap::iterator n_col = node2columnMap.find( nIn );
00229 if ( n_col != node2columnMap.end() ) {
00230 columns[ i ] = & n_col->second;
00231 }
00232 else {
00233 TNodeNodeMap::iterator nInOut = nodeIn2OutMap.find( nIn );
00234 if ( nInOut == nodeIn2OutMap.end() )
00235 RETURN_BAD_RESULT("No matching node for "<< nIn->GetID() <<
00236 " in face "<< face->GetID());
00237 columns[ i ] = makeNodeColumn( node2columnMap, nIn, nInOut->second );
00238 }
00239 }
00240
00241 StdMeshers_Prism_3D::AddPrisms( columns, myHelper );
00242 }
00243 }
00244
00245 return true;
00246 }
00247
00248
00256
00257
00258 TNodeColumn* StdMeshers_RadialPrism_3D::makeNodeColumn( TNode2ColumnMap& n2ColMap,
00259 const SMDS_MeshNode* outNode,
00260 const SMDS_MeshNode* inNode)
00261 {
00262 SMESHDS_Mesh * meshDS = myHelper->GetMeshDS();
00263 int shapeID = myHelper->GetSubShapeID();
00264
00265 if ( myLayerPositions.empty() ) {
00266 gp_Pnt pIn = gpXYZ( inNode ), pOut = gpXYZ( outNode );
00267 computeLayerPositions( pIn, pOut );
00268 }
00269 int nbSegments = myLayerPositions.size() + 1;
00270
00271 TNode2ColumnMap::iterator n_col =
00272 n2ColMap.insert( make_pair( outNode, TNodeColumn() )).first;
00273 TNodeColumn & column = n_col->second;
00274 column.resize( nbSegments + 1 );
00275 column.front() = outNode;
00276 column.back() = inNode;
00277
00278 gp_XYZ p1 = gpXYZ( outNode );
00279 gp_XYZ p2 = gpXYZ( inNode );
00280 for ( int z = 1; z < nbSegments; ++z )
00281 {
00282 double r = myLayerPositions[ z - 1 ];
00283 gp_XYZ p = ( 1 - r ) * p1 + r * p2;
00284 SMDS_MeshNode* n = meshDS->AddNode( p.X(), p.Y(), p.Z() );
00285 meshDS->SetNodeInVolume( n, shapeID );
00286 column[ z ] = n;
00287 }
00288
00289 return & column;
00290 }
00291
00292
00293
00298
00299
00300
00301 class TNodeDistributor: public StdMeshers_Regular_1D
00302 {
00303 list <const SMESHDS_Hypothesis *> myUsedHyps;
00304 public:
00305
00306 static TNodeDistributor* GetDistributor(SMESH_Mesh& aMesh)
00307 {
00308 const int myID = -1000;
00309 map < int, SMESH_1D_Algo * > & algoMap = aMesh.GetGen()->_map1D_Algo;
00310 map < int, SMESH_1D_Algo * >::iterator id_algo = algoMap.find( myID );
00311 if ( id_algo == algoMap.end() )
00312 return new TNodeDistributor( myID, 0, aMesh.GetGen() );
00313 return static_cast< TNodeDistributor* >( id_algo->second );
00314 }
00315
00316 bool Compute( vector< double > & positions,
00317 gp_Pnt pIn,
00318 gp_Pnt pOut,
00319 SMESH_Mesh& aMesh,
00320 const StdMeshers_LayerDistribution* hyp)
00321 {
00322 double len = pIn.Distance( pOut );
00323 if ( len <= DBL_MIN ) return error("Too close points of inner and outer shells");
00324
00325 if ( !hyp || !hyp->GetLayerDistribution() )
00326 return error( "Invalid LayerDistribution hypothesis");
00327 myUsedHyps.clear();
00328 myUsedHyps.push_back( hyp->GetLayerDistribution() );
00329
00330 TopoDS_Edge edge = BRepBuilderAPI_MakeEdge( pIn, pOut );
00331 SMESH_Hypothesis::Hypothesis_Status aStatus;
00332 if ( !StdMeshers_Regular_1D::CheckHypothesis( aMesh, edge, aStatus ))
00333 return error( "StdMeshers_Regular_1D::CheckHypothesis() failed "
00334 "with LayerDistribution hypothesis");
00335
00336 BRepAdaptor_Curve C3D(edge);
00337 double f = C3D.FirstParameter(), l = C3D.LastParameter();
00338 list< double > params;
00339 if ( !StdMeshers_Regular_1D::computeInternalParameters( aMesh, C3D, len, f, l, params, false ))
00340 return error("StdMeshers_Regular_1D failed to compute layers distribution");
00341
00342 positions.clear();
00343 positions.reserve( params.size() );
00344 for (list<double>::iterator itU = params.begin(); itU != params.end(); itU++)
00345 positions.push_back( *itU / len );
00346 return true;
00347 }
00348 protected:
00349
00350 TNodeDistributor( int hypId, int studyId, SMESH_Gen* gen)
00351 : StdMeshers_Regular_1D( hypId, studyId, gen)
00352 {
00353 }
00354
00355 virtual const list <const SMESHDS_Hypothesis *> &
00356 GetUsedHypothesis(SMESH_Mesh &, const TopoDS_Shape &, const bool)
00357 {
00358 return myUsedHyps;
00359 }
00360
00361 };
00362
00363
00368
00369
00370 bool StdMeshers_RadialPrism_3D::computeLayerPositions(const gp_Pnt& pIn,
00371 const gp_Pnt& pOut)
00372 {
00373 if ( myNbLayerHypo )
00374 {
00375 int nbSegments = myNbLayerHypo->GetNumberOfLayers();
00376 myLayerPositions.resize( nbSegments - 1 );
00377 for ( int z = 1; z < nbSegments; ++z )
00378 myLayerPositions[ z - 1 ] = double( z )/ double( nbSegments );
00379 return true;
00380 }
00381 if ( myDistributionHypo ) {
00382 SMESH_Mesh * mesh = myHelper->GetMesh();
00383 if ( !TNodeDistributor::GetDistributor(*mesh)->Compute( myLayerPositions, pIn, pOut,
00384 *mesh, myDistributionHypo ))
00385 {
00386 error( TNodeDistributor::GetDistributor(*mesh)->GetComputeError() );
00387 return false;
00388 }
00389 }
00390 RETURN_BAD_RESULT("Bad hypothesis");
00391 }
00392
00393
00394
00395
00396
00397
00398
00399 bool StdMeshers_RadialPrism_3D::Evaluate(SMESH_Mesh& aMesh,
00400 const TopoDS_Shape& aShape,
00401 MapShapeNbElems& aResMap)
00402 {
00403
00404 TopoDS_Solid solid = TopoDS::Solid( aShape );
00405 TopoDS_Shell outerShell = BRepTools::OuterShell( solid );
00406 TopoDS_Shape innerShell;
00407 int nbShells = 0;
00408 for ( TopoDS_Iterator It (solid); It.More(); It.Next(), ++nbShells )
00409 if ( !outerShell.IsSame( It.Value() ))
00410 innerShell = It.Value();
00411 if ( nbShells != 2 ) {
00412 std::vector<int> aResVec(SMDSEntity_Last);
00413 for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
00414 SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
00415 aResMap.insert(std::make_pair(sm,aResVec));
00416 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
00417 smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
00418 return false;
00419 }
00420
00421
00422 TAssocTool::TShapeShapeMap shape2ShapeMap;
00423 if ( !TAssocTool::FindSubShapeAssociation( outerShell, &aMesh,
00424 innerShell, &aMesh,
00425 shape2ShapeMap) ) {
00426 std::vector<int> aResVec(SMDSEntity_Last);
00427 for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
00428 SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
00429 aResMap.insert(std::make_pair(sm,aResVec));
00430 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
00431 smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
00432 return false;
00433 }
00434
00435
00436 int nb0d_Out=0, nb2d_3_Out=0, nb2d_4_Out=0;
00437
00438 for (TopExp_Explorer exp(outerShell, TopAbs_FACE); exp.More(); exp.Next()) {
00439
00440 SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
00441 MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
00442 std::vector<int> aVec = (*anIt).second;
00443 nb0d_Out += aVec[SMDSEntity_Node];
00444 nb2d_3_Out += Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
00445 nb2d_4_Out += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
00446 }
00447 int nb1d_Out = 0;
00448 TopTools_MapOfShape tmpMap;
00449 for (TopExp_Explorer exp(outerShell, TopAbs_EDGE); exp.More(); exp.Next()) {
00450 if( tmpMap.Contains( exp.Current() ) )
00451 continue;
00452 tmpMap.Add( exp.Current() );
00453 SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
00454 MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
00455 std::vector<int> aVec = (*anIt).second;
00456 nb0d_Out += aVec[SMDSEntity_Node];
00457 nb1d_Out += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
00458 }
00459 tmpMap.Clear();
00460 for (TopExp_Explorer exp(outerShell, TopAbs_VERTEX); exp.More(); exp.Next()) {
00461 if( tmpMap.Contains( exp.Current() ) )
00462 continue;
00463 tmpMap.Add( exp.Current() );
00464 nb0d_Out++;
00465 }
00466
00467
00468 int nb0d_In=0, nb2d_3_In=0, nb2d_4_In=0;
00469
00470 for (TopExp_Explorer exp(innerShell, TopAbs_FACE); exp.More(); exp.Next()) {
00471
00472 SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
00473 MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
00474 std::vector<int> aVec = (*anIt).second;
00475 nb0d_In += aVec[SMDSEntity_Node];
00476 nb2d_3_In += Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
00477 nb2d_4_In += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
00478 }
00479 int nb1d_In = 0;
00480 tmpMap.Clear();
00481 bool IsQuadratic = false;
00482 bool IsFirst = true;
00483 for (TopExp_Explorer exp(innerShell, TopAbs_EDGE); exp.More(); exp.Next()) {
00484 if( tmpMap.Contains( exp.Current() ) )
00485 continue;
00486 tmpMap.Add( exp.Current() );
00487 SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
00488 MapShapeNbElemsItr anIt = aResMap.find(aSubMesh);
00489 std::vector<int> aVec = (*anIt).second;
00490 nb0d_In += aVec[SMDSEntity_Node];
00491 nb1d_In += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
00492 if(IsFirst) {
00493 IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
00494 IsFirst = false;
00495 }
00496 }
00497 tmpMap.Clear();
00498 for (TopExp_Explorer exp(innerShell, TopAbs_VERTEX); exp.More(); exp.Next()) {
00499 if( tmpMap.Contains( exp.Current() ) )
00500 continue;
00501 tmpMap.Add( exp.Current() );
00502 nb0d_In++;
00503 }
00504
00505 bool IsOK = (nb0d_Out==nb0d_In) && (nb1d_Out==nb1d_In) &&
00506 (nb2d_3_Out==nb2d_3_In) && (nb2d_4_Out==nb2d_4_In);
00507 if(!IsOK) {
00508 std::vector<int> aResVec(SMDSEntity_Last);
00509 for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
00510 SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
00511 aResMap.insert(std::make_pair(sm,aResVec));
00512 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
00513 smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
00514 return false;
00515 }
00516
00517 int nbLayers = 0;
00518 if( myNbLayerHypo ) {
00519 nbLayers = myNbLayerHypo->GetNumberOfLayers();
00520 }
00521 if ( myDistributionHypo ) {
00522 if ( !myDistributionHypo->GetLayerDistribution() ) {
00523 std::vector<int> aResVec(SMDSEntity_Last);
00524 for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
00525 SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
00526 aResMap.insert(std::make_pair(sm,aResVec));
00527 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
00528 smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
00529 return false;
00530 }
00531 TopExp_Explorer exp(outerShell, TopAbs_VERTEX);
00532 TopoDS_Vertex Vout = TopoDS::Vertex(exp.Current());
00533 TopoDS_Vertex Vin = TopoDS::Vertex( shape2ShapeMap(Vout) );
00534 if ( myLayerPositions.empty() ) {
00535 gp_Pnt pIn = BRep_Tool::Pnt(Vin);
00536 gp_Pnt pOut = BRep_Tool::Pnt(Vout);
00537 computeLayerPositions( pIn, pOut );
00538 }
00539 nbLayers = myLayerPositions.size() + 1;
00540 }
00541
00542 std::vector<int> aResVec(SMDSEntity_Last);
00543 for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
00544 if(IsQuadratic) {
00545 aResVec[SMDSEntity_Quad_Penta] = nb2d_3_Out * nbLayers;
00546 aResVec[SMDSEntity_Quad_Hexa] = nb2d_4_Out * nbLayers;
00547 int nb1d = ( nb2d_3_Out*3 + nb2d_4_Out*4 ) / 2;
00548 aResVec[SMDSEntity_Node] = nb0d_Out * ( 2*nbLayers - 1 ) - nb1d * nbLayers;
00549 }
00550 else {
00551 aResVec[SMDSEntity_Node] = nb0d_Out * ( nbLayers - 1 );
00552 aResVec[SMDSEntity_Penta] = nb2d_3_Out * nbLayers;
00553 aResVec[SMDSEntity_Hexa] = nb2d_4_Out * nbLayers;
00554 }
00555 SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
00556 aResMap.insert(std::make_pair(sm,aResVec));
00557
00558 return true;
00559 }