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