src/SMESH_I/SMESH_Measurements_i.cxx

   1 // Copyright (C) 2007-2021  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, or (at your option) any later version.
  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 //  File   : SMESH_Measurements_i.cxx
  23 //  Author : Pavel TELKOV, Open CASCADE S.A.S. (pavel.telkov@opencascade.com)
  24
  25 #ifdef WIN32
  26 #define NOMINMAX
  27 #endif
  28
  29 #include "SMESH_Measurements_i.hxx"
  30
  31 #include "SMDS_ElemIterator.hxx"
  32 #include "SMDS_Mesh.hxx"
  33 #include "SMDS_MeshElement.hxx"
  34 #include "SMDS_MeshNode.hxx"
  35 #include "SMESHDS_Mesh.hxx"
  36 #include "SMESH_Filter_i.hxx"
  37 #include "SMESH_Gen_i.hxx"
  38 #include "SMESH_MeshAlgos.hxx"
  39 #include "SMESH_PythonDump.hxx"
  40
  41 #include <cmath>
  42
  43 using namespace SMESH;
  44
  45 /**
  46  * this local function to avoid uninitialized fields
  47  */
  48 static void initMeasure( SMESH::Measure& theMeasure)
  49 {
  50
  51   theMeasure.minX = theMeasure.minY = theMeasure.minZ = 0.;
  52   theMeasure.maxX = theMeasure.maxY = theMeasure.maxZ = 0.;
  53   theMeasure.node1 = theMeasure.node2 = -1;
  54   theMeasure.elem1 = theMeasure.elem2 = -1;
  55   theMeasure.value = 0.;
  56 }
  57
  58 //=============================================================================
  59 /*!
  60  *  SMESH_Gen_i::CreateMeasurements
  61  *
  62  *  Create measurement instance
  63  */
  64 //=============================================================================
  65
  66 SMESH::Measurements_ptr SMESH_Gen_i::CreateMeasurements()
  67 {
  68   SMESH::Measurements_i* aMeasure = new SMESH::Measurements_i();
  69   SMESH::Measurements_var anObj = aMeasure->_this();
  70   return anObj._retn();
  71 }
  72
  73
  74 /*
  75   Class       : Measurements
  76   Description : make measure of mesh qunatities
  77 */
  78
  79 //=======================================================================
  80 // name    : Measurements_i
  81 // Purpose : Constructor
  82 //=======================================================================
  83 Measurements_i::Measurements_i()
  84 : SALOME::GenericObj_i( SMESH_Gen_i::GetPOA() )
  85 {
  86   //Base class Salome_GenericObject do it inmplicitly by overriding PortableServer::POA_ptr _default_POA() method
  87   //PortableServer::ObjectId_var anObjectId =
  88   //  SMESH_Gen_i::GetPOA()->activate_object( this );
  89 }
  90
  91 //=======================================================================
  92 // name    : ~Measurements_i
  93 // Purpose : Destructor
  94 //=======================================================================
  95 Measurements_i::~Measurements_i()
  96 {
  97   //TPythonDump()<<this<<".UnRegister()";
  98 }
  99
 100 static bool getNodeNodeDistance (SMESH::Measure& theMeasure,
 101                                  const SMDS_MeshNode* theNode1,
 102                                  const SMDS_MeshNode* theNode2 = 0)
 103 {
 104   double dist = 0., dd = 0.;
 105
 106   if (!theNode1)
 107     return false;
 108
 109   dd = theNode1->X(); if (theNode2) dd -= theNode2->X(); theMeasure.minX = dd; dd *= dd; dist += dd;
 110   dd = theNode1->Y(); if (theNode2) dd -= theNode2->Y(); theMeasure.minY = dd; dd *= dd; dist += dd;
 111   dd = theNode1->Z(); if (theNode2) dd -= theNode2->Z(); theMeasure.minZ = dd; dd *= dd; dist += dd;
 112
 113   if (dist < 0)
 114     return false;
 115
 116   theMeasure.value = sqrt(dist);
 117   theMeasure.node1 = theNode1->GetID();
 118   theMeasure.node2 = theNode2 ? theNode2->GetID() : 0;
 119
 120   return true;
 121 }
 122
 123 static bool getNodeElemDistance (SMESH::Measure&        theMeasure,
 124                                  const SMDS_MeshNode*   theNode,
 125                                  SMESH_ElementSearcher* theElemSearcher)
 126 {
 127   if ( !theNode || !theElemSearcher )
 128     return false;
 129
 130   const SMDS_MeshElement* closestElement = 0;
 131   gp_Pnt        point = SMESH_NodeXYZ( theNode );
 132   gp_Pnt closestPoint = theElemSearcher->Project( point, SMDSAbs_All, &closestElement );
 133
 134   if ( closestElement )
 135   {
 136     theMeasure.value = point.Distance( closestPoint );
 137     theMeasure.node1 = theNode->GetID();
 138     theMeasure.elem2 = closestElement->GetID();
 139     theMeasure.maxX  = closestPoint.X();
 140     theMeasure.maxY  = closestPoint.Y();
 141     theMeasure.maxZ  = closestPoint.Z();
 142     theMeasure.minX  = closestPoint.X() - point.X();
 143     theMeasure.minY  = closestPoint.Y() - point.Y();
 144     theMeasure.minZ  = closestPoint.Z() - point.Z();
 145   }
 146
 147   return closestElement;
 148 }
 149
 150 static SMESHDS_Mesh* getMesh(SMESH::SMESH_IDSource_ptr theSource)
 151 {
 152   if (!CORBA::is_nil( theSource ))
 153   {
 154     SMESH::SMESH_Mesh_var mesh = theSource->GetMesh();
 155     SMESH_Mesh_i* anImplPtr = DownCast<SMESH_Mesh_i*>( mesh );
 156     if (anImplPtr)
 157       return anImplPtr->GetImpl().GetMeshDS();
 158   }
 159   return 0;
 160 }
 161
 162 static bool isNodeType (SMESH::array_of_ElementType_var theTypes)
 163 {
 164   return theTypes->length() > 0 && theTypes[0] == SMESH::NODE;
 165 }
 166
 167 static double getNumericalValue(SMESH::SMESH_IDSource_ptr            theSource,
 168                                 SMESH::Controls::NumericalFunctorPtr theFunctor)
 169 {
 170   double value = 0;
 171
 172   if ( !CORBA::is_nil( theSource ) ) {
 173     const SMESHDS_Mesh* aMesh = getMesh( theSource );
 174     if ( aMesh ) {
 175       theFunctor->SetMesh( aMesh );
 176
 177       SMESH::smIdType_array_var anElementsId = theSource->GetIDs();
 178       for ( CORBA::ULong i = 0; i < anElementsId->length(); i++) {
 179         value += theFunctor->GetValue( anElementsId[i] );
 180       }
 181     }
 182   }
 183   return value;
 184 }
 185
 186 //=======================================================================
 187 // name    : MinDistance
 188 // Purpose : minimal distance between two given entities
 189 //=======================================================================
 190 SMESH::Measure Measurements_i::MinDistance
 191  (SMESH::SMESH_IDSource_ptr theSource1,
 192   SMESH::SMESH_IDSource_ptr theSource2)
 193 {
 194   SMESH::Measure aMeasure;
 195   initMeasure(aMeasure);
 196
 197   if (CORBA::is_nil( theSource1 ))
 198     return aMeasure;
 199
 200   // if second source is null, min distance from theSource1 to the origin is calculated
 201   bool isOrigin =  CORBA::is_nil( theSource2 );
 202
 203   // calculate minimal distance between two mesh entities
 204   SMESH::array_of_ElementType_var types1 = theSource1->GetTypes();
 205   SMESH::array_of_ElementType_var types2;
 206   if ( !isOrigin ) types2 = theSource2->GetTypes();
 207
 208   // here we assume that type of all IDs defined by first type in array
 209   bool isNode1 = isNodeType(types1);
 210   bool isNode2 = isOrigin || isNodeType(types2);
 211
 212   SMESH::smIdType_array_var aElementsId1 = theSource1->GetIDs();
 213   SMESH::smIdType_array_var aElementsId2;
 214
 215   // compute distance between two entities
 216   /* NOTE: currently only node-to-node case is implemented
 217    * all other cases will be implemented later
 218    * below IF should be replaced by complete switch
 219    * on mesh entities types
 220    */
 221   if (isNode1 && isNode2)
 222   {
 223     // node - node
 224     const SMESHDS_Mesh* aMesh1 = getMesh( theSource1 );
 225     const SMESHDS_Mesh* aMesh2 = isOrigin ? 0 : getMesh( theSource2 );
 226     if ( !isOrigin ) aElementsId2 = theSource2->GetIDs();
 227     const SMDS_MeshNode* theNode1 = aMesh1 ? aMesh1->FindNode( aElementsId1[0] ) : 0;
 228     const SMDS_MeshNode* theNode2 = aMesh2 ? aMesh2->FindNode( aElementsId2[0] ) : 0;
 229     getNodeNodeDistance( aMeasure, theNode1, theNode2 );
 230   }
 231   if (isNode1 && !isNode2 && aElementsId1->length() == 1 )
 232   {
 233     // node - elements
 234     SMESHDS_Mesh* aMesh1 = getMesh( theSource1 );
 235     SMESHDS_Mesh* aMesh2 = getMesh( theSource2 );
 236     if ( aMesh1 && aMesh2 )
 237     {
 238       const SMDS_MeshNode* aNode    = aMesh1->FindNode( aElementsId1[0] );
 239       SMDS_ElemIteratorPtr anElemIt = SMESH_Mesh_i::GetElements( theSource2, SMESH::ALL );
 240       std::unique_ptr< SMESH_ElementSearcher > aSearcher
 241         ( SMESH_MeshAlgos::GetElementSearcher( *aMesh2, anElemIt ));
 242       getNodeElemDistance( aMeasure, aNode, aSearcher.get() );
 243     }
 244   }
 245   else
 246   {
 247     // NOT_IMPLEMENTED
 248   }
 249
 250   return aMeasure;
 251 }
 252
 253 //=======================================================================
 254 // name    : enlargeBoundingBox
 255 // Purpose :
 256 //=======================================================================
 257 static void enlargeBoundingBox(const SMDS_MeshNode* theNode,
 258                                SMESH::Measure&      theMeasure)
 259 {
 260   if (!theNode)
 261     return;
 262   if ( theMeasure.node1 == -1 ) {
 263     // we use this attribute as a flag that it is the first node added to the bnd box
 264     theMeasure.minX = theMeasure.maxX = theNode->X();
 265     theMeasure.minY = theMeasure.maxY = theNode->Y();
 266     theMeasure.minZ = theMeasure.maxZ = theNode->Z();
 267     theMeasure.node1 = theNode->GetID();
 268   }
 269   else {
 270     theMeasure.minX = std::min( theMeasure.minX, theNode->X() );
 271     theMeasure.maxX = std::max( theMeasure.maxX, theNode->X() );
 272     theMeasure.minY = std::min( theMeasure.minY, theNode->Y() );
 273     theMeasure.maxY = std::max( theMeasure.maxY, theNode->Y() );
 274     theMeasure.minZ = std::min( theMeasure.minZ, theNode->Z() );
 275     theMeasure.maxZ = std::max( theMeasure.maxZ, theNode->Z() );
 276   }
 277 }
 278
 279 //=======================================================================
 280 // name    : enlargeBoundingBox
 281 // Purpose :
 282 //=======================================================================
 283 static void enlargeBoundingBox(const SMESH::SMESH_IDSource_ptr theObject,
 284                                SMESH::Measure&                 theMeasure)
 285 {
 286   if ( CORBA::is_nil( theObject ) )
 287     return;
 288   const SMESHDS_Mesh* aMesh = getMesh( theObject );
 289   if ( !aMesh )
 290     return;
 291
 292   if ( DownCast<SMESH_Mesh_i*>( theObject )) // theObject is mesh
 293   {
 294     for (SMDS_NodeIteratorPtr aNodeIter = aMesh->nodesIterator(); aNodeIter->more(); )
 295       enlargeBoundingBox( aNodeIter->next(), theMeasure);
 296   }
 297   else
 298   {
 299     SMESH::array_of_ElementType_var types = theObject->GetTypes();
 300     SMESH::smIdType_array_var aElementsId = theObject->GetIDs();
 301     // here we assume that type of all IDs defined by first type in array
 302     const bool isNode = isNodeType( types );
 303     for(int i = 0, n = aElementsId->length(); i < n; i++)
 304     {
 305       if (isNode)
 306         enlargeBoundingBox( aMesh->FindNode( aElementsId[i] ), theMeasure);
 307       else
 308       {
 309         if ( const SMDS_MeshElement* elem = aMesh->FindElement( aElementsId[i] ))
 310           for (SMDS_NodeIteratorPtr aNodeIter = elem->nodeIterator(); aNodeIter->more(); )
 311             enlargeBoundingBox( aNodeIter->next(), theMeasure);
 312       }
 313     }
 314   }
 315 }
 316
 317 //=======================================================================
 318 // name    : BoundingBox
 319 // Purpose : compute common bounding box of entities
 320 //=======================================================================
 321 SMESH::Measure Measurements_i::BoundingBox (const SMESH::ListOfIDSources& theSources)
 322 {
 323   SMESH::Measure aMeasure;
 324   initMeasure(aMeasure);
 325
 326   // calculate bounding box on sources
 327   for ( int i = 0, n = theSources.length(); i < n ; ++i )
 328     enlargeBoundingBox( theSources[i], aMeasure );
 329
 330   return aMeasure;
 331 }
 332
 333 //=======================================================================
 334 // name    : Length
 335 // Purpose : sum of length of 1D elements of the source
 336 //=======================================================================
 337 double Measurements_i::Length(SMESH::SMESH_IDSource_ptr theSource)
 338 {
 339   return getNumericalValue( theSource, SMESH::Controls::NumericalFunctorPtr(new SMESH::Controls::Length()) );
 340 }
 341
 342 //=======================================================================
 343 // name    : Area
 344 // Purpose : sum of area of 2D elements of the source
 345 //=======================================================================
 346 double Measurements_i::Area(SMESH::SMESH_IDSource_ptr theSource)
 347 {
 348   return getNumericalValue( theSource, SMESH::Controls::NumericalFunctorPtr(new SMESH::Controls::Area()) );
 349 }
 350
 351 //=======================================================================
 352 // name    : Volume
 353 // Purpose : sum of volume of 3D elements of the source
 354 //=======================================================================
 355 double Measurements_i::Volume(SMESH::SMESH_IDSource_ptr theSource)
 356 {
 357   return getNumericalValue( theSource, SMESH::Controls::NumericalFunctorPtr(new SMESH::Controls::Volume()) );
 358 }
 359
 360 //=======================================================================
 361 //function : GravityCenter
 362 //purpose  : return gravity center of the source: average coordinates of all nodes
 363 //=======================================================================
 364
 365 SMESH::PointStruct Measurements_i::GravityCenter(SMESH::SMESH_IDSource_ptr theSource)
 366 {
 367   SMESH::PointStruct grCenter = { 0.,0.,0. };
 368   const SMESHDS_Mesh* mesh = getMesh( theSource );
 369   if ( !mesh )
 370     return grCenter;
 371
 372   // unmark all nodes; visited nodes will be marked
 373   SMESH_MeshAlgos::MarkElems( mesh->nodesIterator(), /*isMarked=*/false );
 374
 375   gp_XYZ sumCoord( 0,0,0 );
 376   int nodeCount = 0;
 377
 378   SMDS_ElemIteratorPtr eIt = SMESH_Mesh_i::GetElements( theSource, SMESH::ALL );
 379   while ( eIt->more() )
 380   {
 381     const SMDS_MeshElement*   elem = eIt->next();
 382     for ( SMDS_NodeIteratorPtr nIt = elem->nodeIterator(); nIt->more(); )
 383     {
 384       const SMDS_MeshNode* n = nIt->next();
 385       if ( !n->isMarked() )
 386       {
 387         sumCoord += SMESH_NodeXYZ( n );
 388         ++nodeCount;
 389         n->setIsMarked( true );
 390       }
 391     }
 392   }
 393   sumCoord /= nodeCount;
 394
 395   grCenter.x = sumCoord.X();
 396   grCenter.y = sumCoord.Y();
 397   grCenter.z = sumCoord.Z();
 398
 399   return grCenter;
 400 }
 401
 402 //=======================================================================
 403 //function : Angle
 404 //purpose  : Return angle in radians defined by 3 points <(p1,p2,p3)
 405 //=======================================================================
 406
 407 CORBA::Double Measurements_i::Angle(const SMESH::PointStruct& p1,
 408                                     const SMESH::PointStruct& p2,
 409                                     const SMESH::PointStruct& p3 )
 410 {
 411   gp_Vec v1( p1.x - p2.x, p1.y - p2.y, p1.z - p2.z );
 412   gp_Vec v2( p3.x - p2.x, p3.y - p2.y, p3.z - p2.z );
 413
 414   double angle = -1;
 415
 416   try
 417   {
 418     angle = v1.Angle( v2 );
 419   }
 420   catch(...)
 421   {
 422   }
 423   if ( std::isnan( angle ))
 424     angle = -1;
 425
 426   return angle;
 427 }