src/MEDCoupling/MEDCouplingGaussLocalization.cxx

   1 // Copyright (C) 2007-2022  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // Author : Anthony Geay (CEA/DEN)
  20
  21 #include "MEDCouplingGaussLocalization.hxx"
  22 #include "InterpKernelGaussCoords.hxx"
  23 #include "MEDCoupling1GTUMesh.hxx"
  24 #include "MEDCouplingUMesh.hxx"
  25 #include "CellModel.hxx"
  26
  27 #include <cmath>
  28 #include <numeric>
  29 #include <sstream>
  30 #include <iterator>
  31 #include <algorithm>
  32
  33 using namespace MEDCoupling;
  34
  35 MEDCouplingGaussLocalization::MEDCouplingGaussLocalization(INTERP_KERNEL::NormalizedCellType type, const std::vector<double>& refCoo,
  36                                                            const std::vector<double>& gsCoo, const std::vector<double>& w)
  37 :_type(type),_ref_coord(refCoo),_gauss_coord(gsCoo),_weight(w)
  38 {
  39   // Will potentially throw (and then release memory for above objects _ref_coord etc ...)
  40   checkConsistencyLight();
  41 }
  42
  43 MEDCouplingGaussLocalization::MEDCouplingGaussLocalization(INTERP_KERNEL::NormalizedCellType typ)
  44 :_type(typ)
  45 {
  46   // Will potentially throw
  47   INTERP_KERNEL::CellModel::GetCellModel(_type);
  48 }
  49
  50 void MEDCouplingGaussLocalization::setType(INTERP_KERNEL::NormalizedCellType typ)
  51 {
  52   INTERP_KERNEL::CellModel::GetCellModel(typ);//throws if not found. This is a check
  53   _type=typ;
  54 }
  55
  56 void MEDCouplingGaussLocalization::checkConsistencyLight() const
  57 {
  58   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(_type);
  59   int nbNodes=cm.getNumberOfNodes();
  60   int dim=cm.getDimension();
  61   if(!cm.isDynamic())
  62     {
  63       if(ToIdType(_ref_coord.size())!=nbNodes*dim)
  64         {
  65           std::ostringstream oss; oss << "Invalid size of refCoo : expecting to be : " << nbNodes << " (nbNodePerCell) * " << dim << " (dim) !";
  66           throw INTERP_KERNEL::Exception(oss.str().c_str());
  67         }
  68     }
  69   if(_gauss_coord.size()!=dim*_weight.size())
  70     {
  71       std::ostringstream oss; oss << "Invalid gsCoo size and weight size : gsCoo.size() must be equal to _weight.size() * " << dim << " (dim) !";
  72       throw INTERP_KERNEL::Exception(oss.str().c_str());
  73     }
  74 }
  75
  76 int MEDCouplingGaussLocalization::getDimension() const
  77 {
  78   if(_weight.empty())
  79     THROW_IK_EXCEPTION("getDimension : weight is empty !");
  80   return (int)_gauss_coord.size()/(int)_weight.size();
  81 }
  82
  83 int MEDCouplingGaussLocalization::getNumberOfPtsInRefCell() const
  84 {
  85   if(_gauss_coord.empty())
  86   {
  87     if(!_weight.empty())
  88       THROW_IK_EXCEPTION("getNumberOfPtsInRefCell : gauss_coords are empty whereas weights are not empty !");
  89     const INTERP_KERNEL::CellModel& cm = INTERP_KERNEL::CellModel::GetCellModel(_type);
  90     return ((int)_ref_coord.size()) / ((int)cm.getDimension());
  91   }
  92   int dim( getDimension() );
  93   return (int)_ref_coord.size()/dim;
  94 }
  95
  96 std::string MEDCouplingGaussLocalization::getStringRepr() const
  97 {
  98   std::ostringstream oss;
  99   oss << "CellType : " << INTERP_KERNEL::CellModel::GetCellModel(_type).getRepr() << std::endl;
 100   oss << "Ref coords : "; std::copy(_ref_coord.begin(),_ref_coord.end(),std::ostream_iterator<double>(oss,", ")); oss << std::endl;
 101   oss << "Localization coords : "; std::copy(_gauss_coord.begin(),_gauss_coord.end(),std::ostream_iterator<double>(oss,", ")); oss << std::endl;
 102   oss << "Weight : "; std::copy(_weight.begin(),_weight.end(),std::ostream_iterator<double>(oss,", ")); oss << std::endl;
 103   return oss.str();
 104 }
 105
 106 std::size_t MEDCouplingGaussLocalization::getMemorySize() const
 107 {
 108   std::size_t ret=0;
 109   ret+=_ref_coord.capacity()*sizeof(double);
 110   ret+=_gauss_coord.capacity()*sizeof(double);
 111   ret+=_weight.capacity()*sizeof(double);
 112   return ret;
 113 }
 114
 115 bool MEDCouplingGaussLocalization::isEqual(const MEDCouplingGaussLocalization& other, double eps) const
 116 {
 117   if(_type!=other._type)
 118     return false;
 119   if(!AreAlmostEqual(_ref_coord,other._ref_coord,eps))
 120     return false;
 121   if(!AreAlmostEqual(_gauss_coord,other._gauss_coord,eps))
 122     return false;
 123   if(!AreAlmostEqual(_weight,other._weight,eps))
 124     return false;
 125   return true;
 126 }
 127
 128 double MEDCouplingGaussLocalization::getRefCoord(int ptIdInCell, int comp) const
 129 {
 130   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(_type);
 131   int nbNodes=cm.getNumberOfNodes();
 132   int dim=cm.getDimension();
 133   if(ptIdInCell<0 || ptIdInCell>=nbNodes)
 134     throw INTERP_KERNEL::Exception("ptIdInCell specified is invalid : must be in [0;nbNodesPerCell) !");
 135   if(comp<0 || comp>=dim)
 136     throw INTERP_KERNEL::Exception("comp specified is invalid : must be in [0:dimOfCell) !");
 137   return _ref_coord[ptIdInCell*dim+comp];
 138 }
 139
 140 double MEDCouplingGaussLocalization::getGaussCoord(int gaussPtIdInCell, int comp) const
 141 {
 142   int dim=checkCoherencyOfRequest(gaussPtIdInCell,comp);
 143   return _gauss_coord[gaussPtIdInCell*dim+comp];
 144 }
 145
 146 double MEDCouplingGaussLocalization::getWeight(int gaussPtIdInCell, double newVal) const
 147 {
 148   checkCoherencyOfRequest(gaussPtIdInCell,0);
 149   return _weight[gaussPtIdInCell];
 150 }
 151
 152 /*!
 153  * Completely useless method for end user. Only for CORBA MPI serialization/unserialization.
 154  * push at the end of tinyInfo its basic serialization info. The size of pushed data is always the same.
 155  * @param tinyInfo inout parameter.
 156  */
 157 void MEDCouplingGaussLocalization::pushTinySerializationIntInfo(std::vector<mcIdType>& tinyInfo) const
 158 {
 159   tinyInfo.push_back(ToIdType(_type));
 160   tinyInfo.push_back(getNumberOfPtsInRefCell());
 161   tinyInfo.push_back(getNumberOfGaussPt());
 162 }
 163
 164 /*!
 165  * Completely useless method for end user. Only for CORBA MPI serialization/unserialization.
 166  * push at the end of tinyInfo its basic serialization info. The size of pushed data is \b NOT always the same contrary to pushTinySerializationIntInfo.
 167  * @param tinyInfo inout parameter.
 168  */
 169 void MEDCouplingGaussLocalization::pushTinySerializationDblInfo(std::vector<double>& tinyInfo) const
 170 {
 171   tinyInfo.insert(tinyInfo.end(),_ref_coord.begin(),_ref_coord.end());
 172   tinyInfo.insert(tinyInfo.end(),_gauss_coord.begin(),_gauss_coord.end());
 173   tinyInfo.insert(tinyInfo.end(),_weight.begin(),_weight.end());
 174 }
 175
 176 /*!
 177  * This method operates the exact inverse operation than MEDCouplingGaussLocalization::pushTinySerializationDblInfo method. This is one of the last step of unserialization process.
 178  * This method should be called on an object resized by buildNewInstanceFromTinyInfo static method.
 179  * This method takes in argument a pointer 'vals' that point to the begin of double data pushed remotely by pushTinySerializationDblInfo method.
 180  * This method returns the pointer 'vals' with an offset of size what it has been read in this method.
 181  */
 182 const double *MEDCouplingGaussLocalization::fillWithValues(const double *vals)
 183 {
 184   const double *work=vals;
 185   std::copy(work,work+_ref_coord.size(),_ref_coord.begin());
 186   work+=_ref_coord.size();
 187   std::copy(work,work+_gauss_coord.size(),_gauss_coord.begin());
 188   work+=_gauss_coord.size();
 189   std::copy(work,work+_weight.size(),_weight.begin());
 190   work+=_weight.size();
 191   return work;
 192 }
 193
 194 /*!
 195  * Given points in \a ptsInRefCoo in the reference coordinates of \a this (in _ref_coord attribute)
 196  * this method computes their coordinates in real world for each cells in \a mesh.
 197  * So the returned array will contain nbCells* \a ptsInRefCoo->getNumberOfTuples() tuples and the number of component
 198  * will be equal to the dimension of \a this.
 199  *
 200  * This method ignores Gauss points in \a this and only those in \a ptsInRefCoo are considered here.
 201  */
 202 MCAuto<DataArrayDouble> MEDCouplingGaussLocalization::localizePtsInRefCooForEachCell(const DataArrayDouble *ptsInRefCoo, const MEDCouplingUMesh *mesh) const
 203 {
 204   if(!ptsInRefCoo || !mesh)
 205     throw INTERP_KERNEL::Exception("MEDCouplingGaussLocalization::localizePtsInRefCooForEachCell : null pointer !");
 206   ptsInRefCoo->checkAllocated();
 207   mesh->checkConsistencyLight();
 208   //
 209   mcIdType nbCells=mesh->getNumberOfCells();
 210   const double *coords(mesh->getCoords()->begin());
 211   const mcIdType *connI(mesh->getNodalConnectivityIndex()->begin()),*conn(mesh->getNodalConnectivity()->begin());
 212   //
 213   mcIdType nbPts(ptsInRefCoo->getNumberOfTuples());
 214   INTERP_KERNEL::NormalizedCellType typ(getType());
 215   int dim(INTERP_KERNEL::CellModel::GetCellModel(typ).getDimension()),outDim(mesh->getSpaceDimension());
 216   MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
 217   ret->alloc(nbPts*nbCells,outDim);
 218   double *retPtr(ret->getPointer());
 219   if(dim!=ToIdType(ptsInRefCoo->getNumberOfComponents()))
 220     throw INTERP_KERNEL::Exception("MEDCouplingGaussLocalization::localizePtsInRefCooForEachCell : number of components of input coo is not equal to dim of element !");
 221   INTERP_KERNEL::GaussCoords calculator;
 222   calculator.addGaussInfo(typ,dim, ptsInRefCoo->begin(),nbPts,&_ref_coord[0],getNumberOfPtsInRefCell());
 223   //
 224   for(mcIdType i=0;i<nbCells;i++,retPtr+=nbPts*outDim)
 225     calculator.calculateCoords(getType(),coords,outDim,conn+connI[i]+1,retPtr);
 226   return ret;
 227 }
 228
 229 /*!
 230  * This method returns an unstructured mesh that represents the reference cell.
 231  */
 232 MCAuto<MEDCouplingUMesh> MEDCouplingGaussLocalization::buildRefCell() const
 233 {
 234   MCAuto<DataArrayDouble> coo(DataArrayDouble::New());
 235   const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel(getType()));
 236   if(getDimension()!=ToIdType(cm.getDimension()))
 237     throw INTERP_KERNEL::Exception("BuildRefCell : dimension mistmatch !");
 238   coo->alloc(cm.getNumberOfNodes(),getDimension());
 239   std::copy(_ref_coord.begin(),_ref_coord.end(),coo->getPointer());
 240   MCAuto<MEDCoupling1SGTUMesh> ret(MEDCoupling1SGTUMesh::New("",getType()));
 241   ret->setCoords(coo);
 242   MCAuto<DataArrayIdType> conn(DataArrayIdType::New());
 243   conn->alloc(cm.getNumberOfNodes(),1);
 244   conn->iota();
 245   ret->setNodalConnectivity(conn);
 246   return MCAuto<MEDCouplingUMesh>(ret->buildUnstructured());
 247 }
 248
 249 /*!
 250  * This method returns an array containing for each Gauss Points in \a this, function values relative to the points of the
 251  * reference cell. Number of components of returned array is equal to the number of points of the reference cell.
 252  */
 253 MCAuto<DataArrayDouble> MEDCouplingGaussLocalization::getShapeFunctionValues() const
 254 {
 255   MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
 256   int nbGaussPt(getNumberOfGaussPt()),nbPtsRefCell(getNumberOfPtsInRefCell()),dim(getDimension());
 257   ret->alloc(nbGaussPt,nbPtsRefCell);
 258   double *retPtr(ret->getPointer());
 259   for(int iGaussPt = 0 ; iGaussPt < nbGaussPt ; ++iGaussPt)
 260   {
 261     std::vector<double> curGaussPt(_gauss_coord.begin()+iGaussPt*dim,_gauss_coord.begin()+(iGaussPt+1)*dim);
 262     INTERP_KERNEL::GaussInfo gi(_type,curGaussPt,1,_ref_coord,nbPtsRefCell);
 263     gi.initLocalInfo();
 264     const double *funcVal( gi.getFunctionValues(0) );
 265     std::copy(funcVal,funcVal+nbPtsRefCell,retPtr);
 266     retPtr += nbPtsRefCell;
 267   }
 268   return ret;
 269 }
 270
 271 MCAuto<DataArrayDouble> MEDCouplingGaussLocalization::getDerivativeOfShapeFunctionValues() const
 272 {
 273   MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
 274   int nbGaussPt(getNumberOfGaussPt()),nbPtsRefCell(getNumberOfPtsInRefCell()),dim(getDimension());
 275   ret->alloc(nbGaussPt,nbPtsRefCell*dim);
 276   double *retPtr(ret->getPointer());
 277   for(int iGaussPt = 0 ; iGaussPt < nbGaussPt ; ++iGaussPt)
 278   {
 279     std::vector<double> curGaussPt(_gauss_coord.begin()+iGaussPt*dim,_gauss_coord.begin()+(iGaussPt+1)*dim);
 280     INTERP_KERNEL::GaussInfo gi(_type,curGaussPt,1,_ref_coord,nbPtsRefCell);
 281     gi.initLocalInfo();
 282     const double *devOfFuncVal( gi.getDerivativeOfShapeFunctionAt(0) );
 283     std::copy(devOfFuncVal,devOfFuncVal+nbPtsRefCell*dim,retPtr);
 284     retPtr += nbPtsRefCell*dim;
 285   }
 286   return ret;
 287 }
 288
 289 /*!
 290  * This method sets the comp_th component of ptIdInCell_th point coordinate of reference element of type this->_type.
 291  * @throw if not 0<=ptIdInCell<nbOfNodePerCell or if not 0<=comp<dim
 292  */
 293 void MEDCouplingGaussLocalization::setRefCoord(int ptIdInCell, int comp, double newVal)
 294 {
 295   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(_type);
 296   int nbNodes=cm.getNumberOfNodes();
 297   int dim=cm.getDimension();
 298   if(ptIdInCell<0 || ptIdInCell>=nbNodes)
 299     throw INTERP_KERNEL::Exception("ptIdInCell specified is invalid : must be in [0;nbNodesPerCell) !");
 300   if(comp<0 || comp>=dim)
 301     throw INTERP_KERNEL::Exception("comp specified is invalid : must be in [0:dimOfCell) !");
 302   _ref_coord[ptIdInCell*dim+comp]=newVal;
 303 }
 304
 305 void MEDCouplingGaussLocalization::setGaussCoord(int gaussPtIdInCell, int comp, double newVal)
 306 {
 307   int dim=checkCoherencyOfRequest(gaussPtIdInCell,comp);
 308   _gauss_coord[gaussPtIdInCell*dim+comp]=newVal;
 309 }
 310
 311 void MEDCouplingGaussLocalization::setWeight(int gaussPtIdInCell, double newVal)
 312 {
 313   checkCoherencyOfRequest(gaussPtIdInCell,0);
 314   _weight[gaussPtIdInCell]=newVal;
 315 }
 316
 317 void MEDCouplingGaussLocalization::setRefCoords(const std::vector<double>& refCoo)
 318 {
 319   _ref_coord=refCoo;
 320 }
 321
 322 void MEDCouplingGaussLocalization::setGaussCoords(const std::vector<double>& gsCoo)
 323 {
 324   _gauss_coord=gsCoo;
 325 }
 326
 327 void MEDCouplingGaussLocalization::setWeights(const std::vector<double>& w)
 328 {
 329   _weight=w;
 330 }
 331
 332 /*!
 333  * The format of 'tinyData' parameter is the same than pushed in method MEDCouplingGaussLocalization::pushTinySerializationIntInfo.
 334  */
 335 MEDCouplingGaussLocalization MEDCouplingGaussLocalization::BuildNewInstanceFromTinyInfo(mcIdType dim, const std::vector<mcIdType>& tinyData)
 336 {
 337   std::vector<double> v1(dim*tinyData[1]),v2(dim*tinyData[2]),v3(tinyData[2]);
 338   return MEDCouplingGaussLocalization((INTERP_KERNEL::NormalizedCellType)tinyData[0],v1,v2,v3);
 339 }
 340
 341 int MEDCouplingGaussLocalization::checkCoherencyOfRequest(mcIdType gaussPtIdInCell, int comp) const
 342 {
 343   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(_type);
 344   int dim=cm.getDimension();
 345   mcIdType nbGsPts=getNumberOfGaussPt();
 346   if(gaussPtIdInCell<0 || gaussPtIdInCell>=nbGsPts)
 347     throw INTERP_KERNEL::Exception("gaussPtIdInCell specified is invalid : must be in [0:nbGsPts) !");
 348   if(comp<0 || comp>=dim)
 349     throw INTERP_KERNEL::Exception("comp specified is invalid : must be in [0:dimOfCell) !");
 350   return dim;
 351 }
 352
 353 bool MEDCouplingGaussLocalization::AreAlmostEqual(const std::vector<double>& v1, const std::vector<double>& v2, double eps)
 354 {
 355   std::size_t sz=v1.size();
 356   if(sz!=v2.size())
 357     return false;
 358   std::vector<double> tmp(sz);
 359   std::transform(v1.begin(),v1.end(),v2.begin(),tmp.begin(),std::minus<double>());
 360   std::transform(tmp.begin(),tmp.end(),tmp.begin(),[](double c){return fabs(c);});
 361   return *std::max_element(tmp.begin(),tmp.end())<eps;
 362 }
 363
 364 MCAuto<DataArrayDouble> MEDCouplingGaussLocalization::GetDefaultReferenceCoordinatesOf(INTERP_KERNEL::NormalizedCellType type)
 365 {
 366   std::vector<double> retCpp(INTERP_KERNEL::GaussInfo::GetDefaultReferenceCoordinatesOf(type));
 367   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
 368   auto nbDim(cm.getDimension());
 369   std::size_t sz(retCpp.size());
 370   MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
 371   if( sz%std::size_t(nbDim) != 0 )
 372     THROW_IK_EXCEPTION("GetDefaultReferenceCoordinatesOf : unexpected size of defaut array : " << sz << " % " << nbDim << " != 0 !");
 373   ret->alloc(sz/size_t(nbDim),nbDim);
 374   std::copy(retCpp.begin(),retCpp.end(),ret->getPointer());
 375   return ret;
 376 }