src/MEDCoupling/MEDCouplingGaussLocalization.cxx

   1 // Copyright (C) 2007-2016  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 "MEDCouplingUMesh.hxx"
  24 #include "CellModel.hxx"
  25
  26 #include <cmath>
  27 #include <numeric>
  28 #include <sstream>
  29 #include <iterator>
  30 #include <algorithm>
  31
  32 using namespace MEDCoupling;
  33
  34 MEDCouplingGaussLocalization::MEDCouplingGaussLocalization(INTERP_KERNEL::NormalizedCellType type, const std::vector<double>& refCoo,
  35                                                            const std::vector<double>& gsCoo, const std::vector<double>& w)
  36 try:_type(type),_ref_coord(refCoo),_gauss_coord(gsCoo),_weight(w)
  37 {
  38   checkConsistencyLight();
  39 }
  40 catch(INTERP_KERNEL::Exception& e)
  41 {
  42     _type=INTERP_KERNEL::NORM_ERROR;
  43     _ref_coord.clear();
  44     _gauss_coord.clear();
  45     _weight.clear();
  46     throw e;
  47 }
  48
  49 MEDCouplingGaussLocalization::MEDCouplingGaussLocalization(INTERP_KERNEL::NormalizedCellType typ)
  50 try:_type(typ)
  51 {
  52   INTERP_KERNEL::CellModel::GetCellModel(_type);
  53 }
  54 catch(INTERP_KERNEL::Exception& e)
  55 {
  56     _type=INTERP_KERNEL::NORM_ERROR;
  57     throw e;
  58 }
  59
  60 void MEDCouplingGaussLocalization::setType(INTERP_KERNEL::NormalizedCellType typ)
  61 {
  62   INTERP_KERNEL::CellModel::GetCellModel(typ);//throws if not found. This is a check
  63   _type=typ;
  64 }
  65
  66 void MEDCouplingGaussLocalization::checkConsistencyLight() const
  67 {
  68   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(_type);
  69   int nbNodes=cm.getNumberOfNodes();
  70   int dim=cm.getDimension();
  71   if(!cm.isDynamic())
  72     {
  73       if((int)_ref_coord.size()!=nbNodes*dim)
  74         {
  75           std::ostringstream oss; oss << "Invalid size of refCoo : expecting to be : " << nbNodes << " (nbNodePerCell) * " << dim << " (dim) !";
  76           throw INTERP_KERNEL::Exception(oss.str().c_str());
  77         }
  78     }
  79   if(_gauss_coord.size()!=dim*_weight.size())
  80     {
  81       std::ostringstream oss; oss << "Invalid gsCoo size and weight size : gsCoo.size() must be equal to _weight.size() * " << dim << " (dim) !";
  82       throw INTERP_KERNEL::Exception(oss.str().c_str());
  83     }
  84 }
  85
  86 int MEDCouplingGaussLocalization::getDimension() const
  87 {
  88   if(_weight.empty())
  89     return -1;
  90   return (int)_gauss_coord.size()/(int)_weight.size();
  91 }
  92
  93 int MEDCouplingGaussLocalization::getNumberOfPtsInRefCell() const
  94 {
  95   int dim=getDimension();
  96   if(dim==0)
  97     return -1;
  98   return (int)_ref_coord.size()/dim;
  99 }
 100
 101 std::string MEDCouplingGaussLocalization::getStringRepr() const
 102 {
 103   std::ostringstream oss;
 104   oss << "CellType : " << INTERP_KERNEL::CellModel::GetCellModel(_type).getRepr() << std::endl;
 105   oss << "Ref coords : "; std::copy(_ref_coord.begin(),_ref_coord.end(),std::ostream_iterator<double>(oss,", ")); oss << std::endl;
 106   oss << "Localization coords : "; std::copy(_gauss_coord.begin(),_gauss_coord.end(),std::ostream_iterator<double>(oss,", ")); oss << std::endl;
 107   oss << "Weight : "; std::copy(_weight.begin(),_weight.end(),std::ostream_iterator<double>(oss,", ")); oss << std::endl;
 108   return oss.str();
 109 }
 110
 111 std::size_t MEDCouplingGaussLocalization::getMemorySize() const
 112 {
 113   std::size_t ret=0;
 114   ret+=_ref_coord.capacity()*sizeof(double);
 115   ret+=_gauss_coord.capacity()*sizeof(double);
 116   ret+=_weight.capacity()*sizeof(double);
 117   return ret;
 118 }
 119
 120 bool MEDCouplingGaussLocalization::isEqual(const MEDCouplingGaussLocalization& other, double eps) const
 121 {
 122   if(_type!=other._type)
 123     return false;
 124   if(!AreAlmostEqual(_ref_coord,other._ref_coord,eps))
 125     return false;
 126   if(!AreAlmostEqual(_gauss_coord,other._gauss_coord,eps))
 127     return false;
 128   if(!AreAlmostEqual(_weight,other._weight,eps))
 129     return false;
 130   return true;
 131 }
 132
 133 double MEDCouplingGaussLocalization::getRefCoord(int ptIdInCell, int comp) const
 134 {
 135   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(_type);
 136   int nbNodes=cm.getNumberOfNodes();
 137   int dim=cm.getDimension();
 138   if(ptIdInCell<0 || ptIdInCell>=nbNodes)
 139     throw INTERP_KERNEL::Exception("ptIdInCell specified is invalid : must be in [0;nbNodesPerCell) !");
 140   if(comp<0 || comp>=dim)
 141     throw INTERP_KERNEL::Exception("comp specified is invalid : must be in [0:dimOfCell) !");
 142   return _ref_coord[ptIdInCell*dim+comp];
 143 }
 144
 145 double MEDCouplingGaussLocalization::getGaussCoord(int gaussPtIdInCell, int comp) const
 146 {
 147   int dim=checkCoherencyOfRequest(gaussPtIdInCell,comp);
 148   return _gauss_coord[gaussPtIdInCell*dim+comp];
 149 }
 150
 151 double MEDCouplingGaussLocalization::getWeight(int gaussPtIdInCell, double newVal) const
 152 {
 153   checkCoherencyOfRequest(gaussPtIdInCell,0);
 154   return _weight[gaussPtIdInCell];
 155 }
 156
 157 /*!
 158  * Completely useless method for end user. Only for CORBA MPI serialization/unserialization.
 159  * push at the end of tinyInfo its basic serialization info. The size of pushed data is always the same.
 160  * @param tinyInfo inout parameter.
 161  */
 162 void MEDCouplingGaussLocalization::pushTinySerializationIntInfo(std::vector<int>& tinyInfo) const
 163 {
 164   tinyInfo.push_back((int)_type);
 165   tinyInfo.push_back(getNumberOfPtsInRefCell());
 166   tinyInfo.push_back(getNumberOfGaussPt());
 167 }
 168
 169 /*!
 170  * Completely useless method for end user. Only for CORBA MPI serialization/unserialization.
 171  * push at the end of tinyInfo its basic serialization info. The size of pushed data is \b NOT always the same contrary to pushTinySerializationIntInfo.
 172  * @param tinyInfo inout parameter.
 173  */
 174 void MEDCouplingGaussLocalization::pushTinySerializationDblInfo(std::vector<double>& tinyInfo) const
 175 {
 176   tinyInfo.insert(tinyInfo.end(),_ref_coord.begin(),_ref_coord.end());
 177   tinyInfo.insert(tinyInfo.end(),_gauss_coord.begin(),_gauss_coord.end());
 178   tinyInfo.insert(tinyInfo.end(),_weight.begin(),_weight.end());
 179 }
 180
 181 /*!
 182  * This method operates the exact inverse operation than MEDCouplingGaussLocalization::pushTinySerializationDblInfo method. This is one of the last step of unserialization process.
 183  * This method should be called on an object resized by buildNewInstanceFromTinyInfo static method.
 184  * This method takes in argument a pointer 'vals' that point to the begin of double data pushed remotely by pushTinySerializationDblInfo method.
 185  * This method returns the pointer 'vals' with an offset of size what it has been read in this method.
 186  */
 187 const double *MEDCouplingGaussLocalization::fillWithValues(const double *vals)
 188 {
 189   const double *work=vals;
 190   std::copy(work,work+_ref_coord.size(),_ref_coord.begin());
 191   work+=_ref_coord.size();
 192   std::copy(work,work+_gauss_coord.size(),_gauss_coord.begin());
 193   work+=_gauss_coord.size();
 194   std::copy(work,work+_weight.size(),_weight.begin());
 195   work+=_weight.size();
 196   return work;
 197 }
 198
 199 /*!
 200  * Given points in \a ptsInRefCoo in the reference coordinates of \a this (in _ref_coord attribute)
 201  * this method computes their coordinates in real world for each cells in \a mesh.
 202  * So the returned array will contain nbCells* \a ptsInRefCoo->getNumberOfTuples() tuples and the number of component
 203  * will be equal to the dimension of \a this.
 204  *
 205  * This method ignores Gauss points in \a this and only those in \a ptsInRefCoo are considered here.
 206  */
 207 MCAuto<DataArrayDouble> MEDCouplingGaussLocalization::localizePtsInRefCooForEachCell(const DataArrayDouble *ptsInRefCoo, const MEDCouplingUMesh *mesh) const
 208 {
 209   if(!ptsInRefCoo || !mesh)
 210     throw INTERP_KERNEL::Exception("MEDCouplingGaussLocalization::localizePtsInRefCooForEachCell : null pointer !");
 211   ptsInRefCoo->checkAllocated();
 212   mesh->checkConsistencyLight();
 213   //
 214   int nbCells(mesh->getNumberOfCells());
 215   const double *coords(mesh->getCoords()->begin());
 216   const int *connI(mesh->getNodalConnectivityIndex()->begin()),*conn(mesh->getNodalConnectivity()->begin());
 217   //
 218   int nbPts(ptsInRefCoo->getNumberOfTuples());
 219   INTERP_KERNEL::NormalizedCellType typ(getType());
 220   int dim(INTERP_KERNEL::CellModel::GetCellModel(typ).getDimension());
 221   MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
 222   ret->alloc(nbPts*nbCells,dim);
 223   double *retPtr(ret->getPointer());
 224   if(dim!=ptsInRefCoo->getNumberOfComponents())
 225     throw INTERP_KERNEL::Exception("MEDCouplingGaussLocalization::localizePtsInRefCooForEachCell : number of components of input coo is not equal to dim of element !");
 226   const std::vector<double>& wg(getWeights());
 227   INTERP_KERNEL::GaussCoords calculator;
 228   calculator.addGaussInfo(typ,dim, ptsInRefCoo->begin(),nbPts,&_ref_coord[0],getNumberOfPtsInRefCell());
 229   //
 230   for(int i=0;i<nbCells;i++,retPtr+=nbPts*dim)
 231     calculator.calculateCoords(getType(),coords,dim,conn+connI[i]+1,retPtr);
 232   return ret;
 233 }
 234
 235 /*!
 236  * This method sets the comp_th component of ptIdInCell_th point coordinate of reference element of type this->_type.
 237  * @throw if not 0<=ptIdInCell<nbOfNodePerCell or if not 0<=comp<dim
 238  */
 239 void MEDCouplingGaussLocalization::setRefCoord(int ptIdInCell, int comp, double newVal)
 240 {
 241   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(_type);
 242   int nbNodes=cm.getNumberOfNodes();
 243   int dim=cm.getDimension();
 244   if(ptIdInCell<0 || ptIdInCell>=nbNodes)
 245     throw INTERP_KERNEL::Exception("ptIdInCell specified is invalid : must be in [0;nbNodesPerCell) !");
 246   if(comp<0 || comp>=dim)
 247     throw INTERP_KERNEL::Exception("comp specified is invalid : must be in [0:dimOfCell) !");
 248   _ref_coord[ptIdInCell*dim+comp]=newVal;
 249 }
 250
 251 void MEDCouplingGaussLocalization::setGaussCoord(int gaussPtIdInCell, int comp, double newVal)
 252 {
 253   int dim=checkCoherencyOfRequest(gaussPtIdInCell,comp);
 254   _gauss_coord[gaussPtIdInCell*dim+comp]=newVal;
 255 }
 256
 257 void MEDCouplingGaussLocalization::setWeight(int gaussPtIdInCell, double newVal)
 258 {
 259   checkCoherencyOfRequest(gaussPtIdInCell,0);
 260   _weight[gaussPtIdInCell]=newVal;
 261 }
 262
 263 void MEDCouplingGaussLocalization::setRefCoords(const std::vector<double>& refCoo)
 264 {
 265   _ref_coord=refCoo;
 266 }
 267
 268 void MEDCouplingGaussLocalization::setGaussCoords(const std::vector<double>& gsCoo)
 269 {
 270   _gauss_coord=gsCoo;
 271 }
 272
 273 void MEDCouplingGaussLocalization::setWeights(const std::vector<double>& w)
 274 {
 275   _weight=w;
 276 }
 277
 278 /*!
 279  * The format of 'tinyData' parameter is the same than pushed in method MEDCouplingGaussLocalization::pushTinySerializationIntInfo.
 280  */
 281 MEDCouplingGaussLocalization MEDCouplingGaussLocalization::BuildNewInstanceFromTinyInfo(int dim, const std::vector<int>& tinyData)
 282 {
 283   std::vector<double> v1(dim*tinyData[1]),v2(dim*tinyData[2]),v3(tinyData[2]);
 284   return MEDCouplingGaussLocalization((INTERP_KERNEL::NormalizedCellType)tinyData[0],v1,v2,v3);
 285 }
 286
 287 int MEDCouplingGaussLocalization::checkCoherencyOfRequest(int gaussPtIdInCell, int comp) const
 288 {
 289   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(_type);
 290   int dim=cm.getDimension();
 291   int nbGsPts=getNumberOfGaussPt();
 292   if(gaussPtIdInCell<0 || gaussPtIdInCell>=nbGsPts)
 293     throw INTERP_KERNEL::Exception("gaussPtIdInCell specified is invalid : must be in [0:nbGsPts) !");
 294   if(comp<0 || comp>=dim)
 295     throw INTERP_KERNEL::Exception("comp specified is invalid : must be in [0:dimOfCell) !");
 296   return dim;
 297 }
 298
 299 bool MEDCouplingGaussLocalization::AreAlmostEqual(const std::vector<double>& v1, const std::vector<double>& v2, double eps)
 300 {
 301   std::size_t sz=v1.size();
 302   if(sz!=v2.size())
 303     return false;
 304   std::vector<double> tmp(sz);
 305   std::transform(v1.begin(),v1.end(),v2.begin(),tmp.begin(),std::minus<double>());
 306   std::transform(tmp.begin(),tmp.end(),tmp.begin(),std::ptr_fun<double,double>(fabs));
 307   return *std::max_element(tmp.begin(),tmp.end())<eps;
 308 }