src/MEDCoupling/MEDCouplingGaussLocalization.cxx

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