src/MEDCoupling/MEDCouplingField.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 "MEDCouplingField.hxx"
  21 #include "MEDCouplingMesh.hxx"
  22 #include "MEDCouplingFieldDiscretization.hxx"
  23
  24 #include <sstream>
  25
  26 using namespace ParaMEDMEM;
  27
  28 bool MEDCouplingField::isEqualIfNotWhy(const MEDCouplingField *other, double meshPrec, double valsPrec, std::string& reason) const throw(INTERP_KERNEL::Exception)
  29 {
  30   if(!other)
  31     throw INTERP_KERNEL::Exception("MEDCouplingField::isEqualIfNotWhy : other instance is NULL !");
  32   std::ostringstream oss; oss.precision(15);
  33   if(_name!=other->_name)
  34     {
  35       oss << "Field names differ : this name = \"" << _name << "\" and other name = \"" << other->_name << "\" !";
  36       reason=oss.str();
  37       return false;
  38     }
  39   if(_desc!=other->_desc)
  40     {
  41       oss << "Field descriptions differ : this description = \"" << _desc << "\" and other description = \"" << other->_desc << "\" !";
  42       reason=oss.str();
  43       return false;
  44     }
  45   if(_nature!=other->_nature)
  46     {
  47       oss << "Field nature differ : this nature = \"" << MEDCouplingNatureOfField::getRepr(_nature) << "\" and other nature = \"" << MEDCouplingNatureOfField::getRepr(other->_nature) << "\" !";
  48       reason=oss.str();
  49       return false;
  50     }
  51   if(!_type->isEqualIfNotWhy(other->_type,valsPrec,reason))
  52     {
  53       reason.insert(0,"Spatial discretizations differ :");
  54       return false;
  55     }
  56   if(_mesh==0 && other->_mesh==0)
  57     return true;
  58   if(_mesh==0 || other->_mesh==0)
  59     {
  60       reason="Only one field between the two this and other has its underlying mesh defined !";
  61       return false;
  62     }
  63   if(_mesh==other->_mesh)
  64     return true;
  65   bool ret=_mesh->isEqualIfNotWhy(other->_mesh,meshPrec,reason);
  66   if(!ret)
  67     reason.insert(0,"Underlying meshes of fields differ for the following reason : ");
  68   return ret;
  69 }
  70
  71 bool MEDCouplingField::isEqual(const MEDCouplingField *other, double meshPrec, double valsPrec) const
  72 {
  73   std::string tmp;
  74   return isEqualIfNotWhy(other,meshPrec,valsPrec,tmp);
  75 }
  76
  77 bool MEDCouplingField::isEqualWithoutConsideringStr(const MEDCouplingField *other, double meshPrec, double valsPrec) const
  78 {
  79   if(!_type->isEqualWithoutConsideringStr(other->_type,valsPrec))
  80     return false;
  81   if(_nature!=other->_nature)
  82     return false;
  83   if(_mesh==0 && other->_mesh==0)
  84     return true;
  85   if(_mesh==0 || other->_mesh==0)
  86     return false;
  87   if(_mesh==other->_mesh)
  88     return true;
  89   return _mesh->isEqualWithoutConsideringStr(other->_mesh,meshPrec);
  90 }
  91
  92 /*!
  93  * This method states if 'this' and 'other' are compatibles each other before performing any treatment.
  94  * This method is good for methods like : mergeFields.
  95  * This method is not very demanding compared to areStrictlyCompatible that is better for operation on fields.
  96  */
  97 bool MEDCouplingField::areCompatibleForMerge(const MEDCouplingField *other) const
  98 {
  99   if(!_type->isEqual(other->_type,1.))
 100     return false;
 101   if(_nature!=other->_nature)
 102     return false;
 103   if(_mesh==other->_mesh)
 104     return true;
 105   return _mesh->areCompatibleForMerge(other->_mesh);
 106 }
 107
 108 /*!
 109  * This method is more strict than MEDCouplingField::areCompatibleForMerge method.
 110  * This method is used for operation on fields to operate a first check before attempting operation.
 111  */
 112 bool MEDCouplingField::areStrictlyCompatible(const MEDCouplingField *other) const
 113 {
 114   if(!_type->isEqual(other->_type,1.e-12))
 115     return false;
 116   if(_nature!=other->_nature)
 117     return false;
 118   return _mesh==other->_mesh;
 119 }
 120
 121 void MEDCouplingField::updateTime() const
 122 {
 123   if(_mesh)
 124     updateTimeWith(*_mesh);
 125   if(_type)
 126     updateTimeWith(*_type);
 127 }
 128
 129 TypeOfField MEDCouplingField::getTypeOfField() const
 130 {
 131   return _type->getEnum();
 132 }
 133
 134 /*!
 135  * This method returns the nature of field. This information is very important during interpolation process using ParaMEDMEM::MEDCouplingRemapper or ParaMEDMEM::InterpKernelDEC.
 136  * In other context than the two mentioned before this attribute of the field is not sensitive. This attribute is not store in MED file in MEDLoader.
 137  * More information of the semantic, and the consequence of this attribute in the result of the interpolation, is available \ref NatureOfField "here".
 138  */
 139 NatureOfField MEDCouplingField::getNature() const
 140 {
 141   return _nature;
 142 }
 143
 144 /*!
 145  * This method set the nature of field in \b this.This  information is very important during interpolation process using ParaMEDMEM::MEDCouplingRemapper or ParaMEDMEM::InterpKernelDEC.
 146  * In other context than the two mentioned before this attribute of the field is not sensitive. This attribute is not store in MED file in MEDLoader.
 147  * More information of the semantic, and the consequence of this attribute in the result of the interpolation, is available \ref TableNatureOfField "here".
 148  */
 149 void MEDCouplingField::setNature(NatureOfField nat) throw(INTERP_KERNEL::Exception)
 150 {
 151   _nature=nat;
 152 }
 153
 154 /*!
 155  * This method returns is case of success an instance of DataArrayDouble the user is in reponsability to deal with.
 156  * If 'this->_mesh' is not set an exception will be thrown.
 157  * For a field on node the array of coords will be returned. For a field on cell a ParaMEDMEM::DataArrayDouble instance
 158  * containing the barycenter of cells will be returned. And for a field on gauss point the explicit position of gauss points.
 159  */
 160 DataArrayDouble *MEDCouplingField::getLocalizationOfDiscr() const throw(INTERP_KERNEL::Exception)
 161 {
 162   if(!_mesh)
 163     throw INTERP_KERNEL::Exception("MEDCouplingField::getLocalizationOfDiscr : No mesh set !");
 164   return _type->getLocalizationOfDiscValues(_mesh);
 165 }
 166
 167 /*!
 168  * This method retrieves the measure field of 'this'. If no '_mesh' is defined an exception will be thrown.
 169  * Warning the retrieved field life cycle is the responsability of caller.
 170  */
 171 MEDCouplingFieldDouble *MEDCouplingField::buildMeasureField(bool isAbs) const throw(INTERP_KERNEL::Exception)
 172 {
 173   if(_mesh==0)
 174     throw INTERP_KERNEL::Exception("MEDCouplingField::getMeasureField : no mesh defined !!!");
 175   return _type->getMeasureField(_mesh,isAbs);
 176 }
 177
 178 void MEDCouplingField::setMesh(const MEDCouplingMesh *mesh)
 179 {
 180   if(mesh!=_mesh)
 181     {
 182       if(_mesh)
 183         _mesh->decrRef();
 184       _mesh=mesh;
 185       if(_mesh)
 186         {
 187           _mesh->incrRef();
 188           updateTimeWith(*_mesh);
 189         }
 190     }
 191 }
 192
 193 /*!
 194  * This method sets gauss localization by geometric type.
 195  * @param type geometric type on which the gauss localization will be set.
 196  * @param refCoo is the reference coordinates of the specified element. Its size has to be equal to nbOfNodesPerCell*dimOfType
 197  * @param gsCoo are the coordinates of Gauss points in reference element specified by 'refCoo'. Its size must be equal to wg.size()*dimOfType
 198  * @param wg are the weights on Gauss points. The size of this array is used to determine the number of Gauss point in the element.
 199  * @throw when size of 'RefCoo' is not valid regarding 'type' parameter, it throws too when the mesh is not set before or if it is not a field on Gauss points.
 200  */
 201 void MEDCouplingField::setGaussLocalizationOnType(INTERP_KERNEL::NormalizedCellType type, const std::vector<double>& refCoo,
 202                                                   const std::vector<double>& gsCoo, const std::vector<double>& wg) throw(INTERP_KERNEL::Exception)
 203 {
 204   if(!_mesh)
 205     throw INTERP_KERNEL::Exception("Mesh has to be set before calling setGaussLocalizationOnType method !");
 206   _type->setGaussLocalizationOnType(_mesh,type,refCoo,gsCoo,wg);
 207 }
 208
 209 /*!
 210  * This method sets on ids defined by [begin;end) their gauss localization. This method checks the coherency of cells ids in [begin;end) and 'refCoo' size.
 211  * If an incoherence appears an exception will be thrown and no seting will be performed.
 212  * An exception is thrown too if [begin,end) has a size lesser than 1.
 213  *
 214  * @param refCoo is the reference coordinates of the specified element. Its size has to be equal to nbOfNodesPerCell*dimOfType
 215  * @param gsCoo are the coordinates of Gauss points in reference element specified by 'refCoo'. Its size must be equal to wg.size()*dimOfType
 216  * @param wg are the weights on Gauss points. The size of this array is used to determine the number of Gauss point in the element.
 217  * @throw when size of 'RefCoo' is not valid regarding cells in [begin,end) parameters, it throws too when the mesh is not set before or if it is not a field on Gauss points.
 218  */
 219 void MEDCouplingField::setGaussLocalizationOnCells(const int *begin, const int *end, const std::vector<double>& refCoo,
 220                                                    const std::vector<double>& gsCoo, const std::vector<double>& wg) throw(INTERP_KERNEL::Exception)
 221 {
 222   if(!_mesh)
 223     throw INTERP_KERNEL::Exception("Mesh has to be set before calling setGaussLocalizationOnCells method !");
 224   _type->setGaussLocalizationOnCells(_mesh,begin,end,refCoo,gsCoo,wg);
 225 }
 226
 227 /*!
 228  * This method resets all Gauss loalizations if any.
 229  */
 230 void MEDCouplingField::clearGaussLocalizations()
 231 {
 232   if(!_mesh)
 233     throw INTERP_KERNEL::Exception("Mesh has to be set before calling clearGaussLocalizations method !");
 234   _type->clearGaussLocalizations();
 235 }
 236
 237 /*!
 238  * This method returns reference to the Gauss localization object corresponding to 'locId' id.
 239  * This method throws an exception if there is no mesh, invalid FieldDescription (different from Gauss) and if 'locId' is invalid because out of range given by
 240  * MEDCouplingField::getNbOfGaussLocalization method.
 241  * Warning this method is not const, so the returned object could be modified without any problem.
 242  */
 243 MEDCouplingGaussLocalization& MEDCouplingField::getGaussLocalization(int locId) throw(INTERP_KERNEL::Exception)
 244 {
 245   if(!_mesh)
 246     throw INTERP_KERNEL::Exception("Mesh has to be set before calling getGaussLocalization method !");
 247   return _type->getGaussLocalization(locId);
 248 }
 249
 250 /*!
 251  * This method returns reference to the Gauss localization object corresponding to 'locId' id.
 252  * This method throws an exception if there is no mesh, invalid FieldDescription (different from Gauss) and if several localization ids have been found
 253  * for a type.
 254  */
 255 int MEDCouplingField::getGaussLocalizationIdOfOneType(INTERP_KERNEL::NormalizedCellType type) const throw(INTERP_KERNEL::Exception)
 256 {
 257   if(!_mesh)
 258     throw INTERP_KERNEL::Exception("Mesh has to be set before calling getGaussLocalizationIdOfOneType method !");
 259   return _type->getGaussLocalizationIdOfOneType(type);
 260 }
 261
 262 /*!
 263  * This method returns number of Gauss localization available. Implicitely all ids in [0,getNbOfGaussLocalization()) is a valid Gauss localisation id.
 264  * This method throws an exception if there is no mesh, invalid FieldDescription (different from Gauss)
 265  */
 266 int MEDCouplingField::getNbOfGaussLocalization() const throw(INTERP_KERNEL::Exception)
 267 {
 268   if(!_mesh)
 269     throw INTERP_KERNEL::Exception("Mesh has to be set before calling getNbOfGaussLocalization method !");
 270   return _type->getNbOfGaussLocalization();
 271 }
 272
 273 /*!
 274  * This method returns an id of Gauss localization in [0,getNbOfGaussLocalization()) that corresponds to the localization of the cell specified by its cellId.
 275  * This methods throws an exception if there is no mesh, invalid FieldDescription (different from Gauss) or if at the cell with id 'cellId' in this->_mesh no
 276  * Gauss localization has been set.
 277  */
 278 int MEDCouplingField::getGaussLocalizationIdOfOneCell(int cellId) const throw(INTERP_KERNEL::Exception)
 279 {
 280   if(!_mesh)
 281     throw INTERP_KERNEL::Exception("Mesh has to be set before calling getGaussLocalizationIdOfOneCell method !");
 282   return _type->getGaussLocalizationIdOfOneCell(cellId);
 283 }
 284
 285 /*!
 286  * This method returns all cellIds that share the same Gauss localization given by 'locId' parameter (in range [0,getNbOfGaussLocalization()) ).
 287  * If no cells fit the Gauss localization given by 'locId' cellIds will be returned empty.
 288  * @param locId input that specifies the id of Gauss localization.
 289  * @param cellIds output parameter, that will contain the result if this method succeds. This parameter is systematically cleared when called.
 290  * @throw  if there is no mesh, invalid FieldDescription (different from Gauss) or if locId not in [0,getNbOfGaussLocalization())
 291  */
 292 void MEDCouplingField::getCellIdsHavingGaussLocalization(int locId, std::vector<int>& cellIds) const throw(INTERP_KERNEL::Exception)
 293 {
 294   cellIds.clear();
 295   if(!_mesh)
 296     throw INTERP_KERNEL::Exception("Mesh has to be set before calling getGaussLocalizationIdOfOneCell method !");
 297   _type->getCellIdsHavingGaussLocalization(locId,cellIds);
 298 }
 299
 300 /*!
 301  * This method returns reference to the Gauss localization object corresponding to 'locId' id.
 302  * This method throws an exception if there is no mesh, invalid FieldDescription (different from Gauss) and if 'locId' is invalid because out of range given by
 303  * MEDCouplingField::getNbOfGaussLocalization method.
 304  * Warning this method is const.
 305  */
 306 const MEDCouplingGaussLocalization& MEDCouplingField::getGaussLocalization(int locId) const throw(INTERP_KERNEL::Exception)
 307 {
 308   if(!_mesh)
 309     throw INTERP_KERNEL::Exception("Mesh has to be set before calling getGaussLocalization method !");
 310   return _type->getGaussLocalization(locId);
 311 }
 312
 313 MEDCouplingField::~MEDCouplingField()
 314 {
 315   if(_mesh)
 316     _mesh->decrRef();
 317 }
 318
 319 MEDCouplingField::MEDCouplingField(MEDCouplingFieldDiscretization *type, NatureOfField nature):_nature(nature),_mesh(0),_type(type)
 320 {
 321 }
 322
 323 MEDCouplingField::MEDCouplingField(TypeOfField type):_nature(NoNature),_mesh(0),_type(MEDCouplingFieldDiscretization::New(type))
 324 {
 325 }
 326
 327 MEDCouplingField::MEDCouplingField(const MEDCouplingField& other):RefCountObject(other),_name(other._name),_desc(other._desc),_nature(other._nature),
 328                                                                   _mesh(0),_type(other._type->clone())
 329 {
 330   if(other._mesh)
 331     {
 332       _mesh=other._mesh;
 333       _mesh->incrRef();
 334     }
 335 }
 336
 337 /*!
 338  * This method returns a submesh of 'mesh' instance constituting cell ids contained in array defined as an interval [start;end).
 339  * @param di is an array returned that specifies entity ids (nodes, cells ids...) in mesh 'mesh' of entity in returned submesh.
 340  */
 341 MEDCouplingMesh *MEDCouplingField::buildSubMeshData(const int *start, const int *end, DataArrayInt *&di) const
 342 {
 343   return _type->buildSubMeshData(_mesh,start,end,di);
 344 }
 345
 346 /*!
 347  * This method returns tuples ids implied by the mesh selection of the  cell ids contained in array defined as an interval [start;end).
 348  * \return a newly allocated DataArrayInt instance containing tuples ids.
 349  */
 350 DataArrayInt *MEDCouplingField::computeTupleIdsToSelectFromCellIds(const int *startCellIds, const int *endCellIds) const
 351 {
 352   return _type->computeTupleIdsToSelectFromCellIds(_mesh,startCellIds,endCellIds);
 353 }
 354
 355 /*!
 356  * This method returns number of tuples expected regarding its discretization and its _mesh attribute.
 357  * This method expected a not null _mesh instance. If null, an exception will be thrown.
 358  */
 359 int MEDCouplingField::getNumberOfTuplesExpected() const throw(INTERP_KERNEL::Exception)
 360 {
 361   if(_mesh)
 362     return _type->getNumberOfTuples(_mesh);
 363   else
 364     throw INTERP_KERNEL::Exception("MEDCouplingField::getNumberOfTuplesExpected : Empty mesh !");
 365 }
 366
 367 /*!
 368  * This method returns number of mesh placed expected regarding its discretization and its _mesh attribute.
 369  * This method expected a not null _mesh instance. If null, an exception will be thrown.
 370  */
 371 int MEDCouplingField::getNumberOfMeshPlacesExpected() const throw(INTERP_KERNEL::Exception)
 372 {
 373   if(_mesh)
 374     return _type->getNumberOfMeshPlaces(_mesh);
 375   else
 376     throw INTERP_KERNEL::Exception("MEDCouplingField::getNumberOfMeshPlacesExpected : Empty mesh !");
 377 }