src/MEDCoupling/MEDCouplingGaussLocalization.cxx

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