src/MEDCoupling/MEDCouplingMesh.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 // Author : Anthony Geay (CEA/DEN)
  20
  21 #include "MEDCouplingMesh.hxx"
  22 #include "MEDCouplingUMesh.hxx"
  23 #include "MEDCouplingMemArray.hxx"
  24 #include "MEDCouplingFieldDouble.hxx"
  25 #include "MEDCouplingFieldDiscretization.hxx"
  26 #include "MEDCouplingAutoRefCountObjectPtr.hxx"
  27
  28 #include <set>
  29 #include <cmath>
  30 #include <sstream>
  31 #include <fstream>
  32 #include <iterator>
  33
  34 using namespace ParaMEDMEM;
  35
  36 MEDCouplingMesh::MEDCouplingMesh():_time(0.),_iteration(-1),_order(-1)
  37 {
  38 }
  39
  40 MEDCouplingMesh::MEDCouplingMesh(const MEDCouplingMesh& other):_name(other._name),_description(other._description),
  41                                                                _time(other._time),_iteration(other._iteration),
  42                                                                _order(other._order),_time_unit(other._time_unit)
  43 {
  44 }
  45
  46 std::size_t MEDCouplingMesh::getHeapMemorySize() const
  47 {
  48   return _name.capacity()+_description.capacity()+_time_unit.capacity();
  49 }
  50
  51 /*!
  52  * This method is only for ParaMEDMEM in ParaFIELD constructor.
  53  */
  54 bool MEDCouplingMesh::isStructured() const
  55 {
  56   return getType()==CARTESIAN;
  57 }
  58
  59 bool MEDCouplingMesh::isEqualIfNotWhy(const MEDCouplingMesh *other, double prec, std::string& reason) const throw(INTERP_KERNEL::Exception)
  60 {
  61   if(!other)
  62     throw INTERP_KERNEL::Exception("MEDCouplingMesh::isEqualIfNotWhy : other instance is NULL !");
  63   std::ostringstream oss; oss.precision(15);
  64   if(_name!=other->_name)
  65     {
  66       oss << "Mesh names differ : this name = \"" << _name << "\" and other name = \"" << other->_name << "\" !";
  67       reason=oss.str();
  68       return false;
  69     }
  70   if(_description!=other->_description)
  71     {
  72       oss << "Mesh descriptions differ : this description = \"" << _description << "\" and other description = \"" << other->_description << "\" !";
  73       reason=oss.str();
  74       return false;
  75     }
  76   if(_iteration!=other->_iteration)
  77     {
  78       oss << "Mesh iterations differ : this iteration = \"" << _iteration << "\" and other iteration = \"" << other->_iteration << "\" !";
  79       reason=oss.str();
  80       return false;
  81     }
  82   if(_order!=other->_order)
  83     {
  84       oss << "Mesh orders differ : this order = \"" << _order << "\" and other order = \"" << other->_order << "\" !";
  85       reason=oss.str();
  86       return false;
  87     }
  88   if(_time_unit!=other->_time_unit)
  89     {
  90       oss << "Mesh time units differ : this time unit = \"" << _time_unit << "\" and other time unit = \"" << other->_time_unit << "\" !";
  91       reason=oss.str();
  92       return false;
  93     }
  94   if(fabs(_time-other->_time)>=1e-12)
  95     {
  96       oss << "Mesh times differ : this time = \"" << _time << "\" and other time = \"" << other->_time << "\" !";
  97       reason=oss.str();
  98       return false;
  99     }
 100   return true;
 101 }
 102
 103 bool MEDCouplingMesh::isEqual(const MEDCouplingMesh *other, double prec) const throw(INTERP_KERNEL::Exception)
 104 {
 105   std::string tmp;
 106   return isEqualIfNotWhy(other,prec,tmp);
 107 }
 108
 109 /*!
 110  * This method checks geo equivalence between two meshes : 'this' and 'other'.
 111  * If no exception is throw 'this' and 'other' are geometrically equivalent regarding 'levOfCheck' level.
 112  * This method is typically used to change the mesh of a field "safely" depending the 'levOfCheck' level considered.
 113  *
 114  * @param levOfCheck input that specifies the level of check specified. The possible values are listed below.
 115  * @param prec input that specifies precision for double float data used for comparison in meshes.
 116  * @param cellCor output array not always informed (depending 'levOfCheck' param) that gives the corresponding array for cells from 'other' to 'this'.
 117  * @param nodeCor output array not always informed (depending 'levOfCheck' param) that gives the corresponding array for nodes from 'other' to 'this'.
 118  *
 119  * Possible values for levOfCheck :
 120  *   - 0 for strict equality. This is the strongest level. 'cellCor' and 'nodeCor' params are never informed.
 121  *   - 10,11,12 for less strict equality. Two meshes are compared geometrically. In case of success 'cellCor' and 'nodeCor' are informed. Warning ! These equivalences are CPU/Mem costly. The 3 values correspond respectively to policy used for cell comparison (see MEDCouplingUMesh::zipConnectivityTraducer to have more details)
 122  *   - 20,21,22, for less strict equality. Two meshes are compared geometrically. The difference with the previous version is that nodes(coordinates) are expected to be the same between this and other. In case of success 'cellCor' is informed. Warning ! These equivalences are CPU/Mem costly. The 3 values correspond respectively to policy used for cell comparison (see MEDCouplingUMesh::zipConnectivityTraducer to have more details)
 123  *   - 1 for fast 'equality'. This is a lazy level. Just number of cells and number of nodes are considered here and 3 cells (begin,middle,end)
 124  *   - 2 for deep 'equality' as 0 option except that no control is done on all strings in mesh.
 125  */
 126 void MEDCouplingMesh::checkGeoEquivalWith(const MEDCouplingMesh *other, int levOfCheck, double prec,
 127                                           DataArrayInt *&cellCor, DataArrayInt *&nodeCor) const throw(INTERP_KERNEL::Exception)
 128 {
 129   cellCor=0;
 130   nodeCor=0;
 131   if(this==other)
 132     return ;
 133   switch(levOfCheck)
 134     {
 135     case 0:
 136       {
 137         if(!isEqual(other,prec))
 138           throw INTERP_KERNEL::Exception("checkGeoFitWith : Meshes are not equal !");
 139         return ;
 140       }
 141     case 10:
 142     case 11:
 143     case 12:
 144       {
 145         checkDeepEquivalWith(other,levOfCheck-10,prec,cellCor,nodeCor);
 146         return ;
 147       }
 148     case 20:
 149     case 21:
 150     case 22:
 151       {
 152         checkDeepEquivalOnSameNodesWith(other,levOfCheck-20,prec,cellCor);
 153         return ;
 154       }
 155     case 1:
 156       {
 157         checkFastEquivalWith(other,prec);
 158         return;
 159       }
 160     case 2:
 161       {
 162         if(!isEqualWithoutConsideringStr(other,prec))
 163           throw INTERP_KERNEL::Exception("checkGeoFitWith : Meshes are not equal without considering strings !");
 164         return ;
 165       }
 166     default:
 167       throw INTERP_KERNEL::Exception("checkGeoFitWith : Invalid levOfCheck specified ! Value must be in 0,1,2,10,11 or 12.");
 168     }
 169 }
 170
 171 /*!
 172  * Given a nodeIds range ['partBg','partEnd'), this method returns the set of cell ids in ascendant order whose connectivity of
 173  * these cells are fully included in the range. As a consequence the returned set of cell ids does \b not \b always fit the nodes in ['partBg','partEnd')
 174  * This method returns the corresponding cells in a newly created array that the caller has the responsability.
 175  */
 176 DataArrayInt *MEDCouplingMesh::getCellIdsFullyIncludedInNodeIds(const int *partBg, const int *partEnd) const
 177 {
 178   std::vector<int> crest;
 179   std::set<int> p(partBg,partEnd);
 180   int nbOfCells=getNumberOfCells();
 181   for(int i=0;i<nbOfCells;i++)
 182     {
 183       std::vector<int> conn;
 184       getNodeIdsOfCell(i,conn);
 185       bool cont=true;
 186       for(std::vector<int>::const_iterator iter=conn.begin();iter!=conn.end() && cont;iter++)
 187         if(p.find(*iter)==p.end())
 188           cont=false;
 189       if(cont)
 190         crest.push_back(i);
 191     }
 192   DataArrayInt *ret=DataArrayInt::New();
 193   ret->alloc((int)crest.size(),1);
 194   std::copy(crest.begin(),crest.end(),ret->getPointer());
 195   return ret;
 196 }
 197
 198 /*!
 199  * This method checks fastly that 'this' and 'other' are equal. All common checks are done here.
 200  */
 201 void MEDCouplingMesh::checkFastEquivalWith(const MEDCouplingMesh *other, double prec) const throw(INTERP_KERNEL::Exception)
 202 {
 203   if(getMeshDimension()!=other->getMeshDimension())
 204     throw INTERP_KERNEL::Exception("checkFastEquivalWith : Mesh dimensions are not equal !");
 205   if(getSpaceDimension()!=other->getSpaceDimension())
 206     throw INTERP_KERNEL::Exception("checkFastEquivalWith : Space dimensions are not equal !");
 207   if(getNumberOfCells()!=other->getNumberOfCells())
 208     throw INTERP_KERNEL::Exception("checkFastEquivalWith : number of cells are not equal !");
 209 }
 210
 211 /*!
 212  * This method is very poor and looks only if 'this' and 'other' are candidate for merge of fields lying repectively on them.
 213  */
 214 bool MEDCouplingMesh::areCompatibleForMerge(const MEDCouplingMesh *other) const
 215 {
 216   if(getMeshDimension()!=other->getMeshDimension())
 217     return false;
 218   if(getSpaceDimension()!=other->getSpaceDimension())
 219     return false;
 220   return true;
 221 }
 222
 223 /*!
 224  * This method builds a field lying on 'this' with 'nbOfComp' components.
 225  * 'func' is a pointer that points to a function that takes 2 arrays in parameter and returns a boolean.
 226  * The first array is a in-param of size this->getSpaceDimension and the second an out param of size 'nbOfComp'.
 227  * The return field will have type specified by 't'. 't' is also used to determine where values of field will be
 228  * evaluate.
 229  * Contrary to other fillFromAnalytic methods this method requests a C++ function pointer as input.
 230  * The 'func' is a callback that takes as first parameter an input array of size 'this->getSpaceDimension()',
 231  * the second parameter is a pointer on a valid zone of size at least equal to 'nbOfComp' values. And too finish
 232  * the returned value is a boolean that is equal to False in case of invalid evaluation (log(0) for example...)
 233  * @param t type of field returned and specifies where the evaluation of func will be done.
 234  * @param nbOfComp number of components of returned field.
 235  * @param func pointer to a function that should return false if the evaluation failed. (division by 0. for example)
 236  * @return field with counter = 1.
 237  */
 238 MEDCouplingFieldDouble *MEDCouplingMesh::fillFromAnalytic(TypeOfField t, int nbOfComp, FunctionToEvaluate func) const
 239 {
 240   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(t,ONE_TIME);
 241   ret->setMesh(this);
 242   ret->fillFromAnalytic(nbOfComp,func);
 243   ret->synchronizeTimeWithSupport();
 244   return ret.retn();
 245 }
 246
 247 /*!
 248  * This method copyies all tiny strings from other (name and components name).
 249  * @throw if other and this have not same mesh type.
 250  */
 251 void MEDCouplingMesh::copyTinyStringsFrom(const MEDCouplingMesh *other) throw(INTERP_KERNEL::Exception)
 252 {
 253   _name=other->_name;
 254   _description=other->_description;
 255   _time_unit=other->_time_unit;
 256 }
 257
 258 /*!
 259  * This method copies all attributes that are \b NOT arrays in this.
 260  * All tiny attributes not usefully for state of 'this' are ignored.
 261  */
 262 void MEDCouplingMesh::copyTinyInfoFrom(const MEDCouplingMesh *other) throw(INTERP_KERNEL::Exception)
 263 {
 264   copyTinyStringsFrom(other);
 265   _time=other->_time;
 266   _iteration=other->_iteration;
 267   _order=other->_order;
 268 }
 269
 270 /*!
 271  * This method builds a field lying on 'this' with 'nbOfComp' components.
 272  * 'func' is a string that is the expression to evaluate.
 273  * The return field will have type specified by 't'. 't' is also used to determine where values of field will be
 274  * evaluate.
 275  * This method is equivalent to those taking a C++ function pointer except that here the 'func' is informed by
 276  * an interpretable input string.
 277  *
 278  * The dynamic interpretor uses \b alphabetical \b order to assign the component id to the var name.
 279  * For example :
 280  * - "2*x+z" func : x stands for component #0 and z stands for component #1 \b NOT #2 !
 281  *
 282  * Some var names are reserved and have special meaning. IVec stands for (1,0,0,...). JVec stands for (0,1,0...).
 283  * KVec stands for (0,0,1,...)... These keywords allows too differentate the evaluation of output components each other.
 284  *
 285  * If 'nbOfComp' equals to 4 for example and that 'this->getSpaceDimension()' equals to 3.
 286  *
 287  * For the input tuple T = (1.,3.,7.) :
 288  *   - '2*x+z' will return (5.,5.,5.,5.)
 289  *   - '2*x+0*y+z' will return (9.,9.,9.,9.)
 290  *   - '2*x*IVec+(x+z)*LVec' will return (2.,0.,0.,4.)
 291  *   - '2*x*IVec+(y+z)*KVec' will return (2.,0.,10.,0.)
 292  *
 293  * @param t type of field returned and specifies where the evaluation of func will be done.
 294  * @param nbOfComp number of components of returned field.
 295  * @param func expression.
 296  * @return field with counter = 1.
 297  */
 298 MEDCouplingFieldDouble *MEDCouplingMesh::fillFromAnalytic(TypeOfField t, int nbOfComp, const char *func) const
 299 {
 300   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(t,ONE_TIME);
 301   ret->setMesh(this);
 302   ret->fillFromAnalytic(nbOfComp,func);
 303   ret->synchronizeTimeWithSupport();
 304   return ret.retn();
 305 }
 306
 307 /*!
 308  * This method builds a field lying on 'this' with 'nbOfComp' components.
 309  * 'func' is a string that is the expression to evaluate.
 310  * The return field will have type specified by 't'. 't' is also used to determine where values of field will be
 311  * evaluate. This method is different than MEDCouplingMesh::fillFromAnalytic, because the info on components are used here to determine vars pos in 'func'.
 312  *
 313  * @param t type of field returned and specifies where the evaluation of func will be done.
 314  * @param nbOfComp number of components of returned field.
 315  * @param func expression.
 316  * @return field with counter = 1.
 317  */
 318 MEDCouplingFieldDouble *MEDCouplingMesh::fillFromAnalytic2(TypeOfField t, int nbOfComp, const char *func) const
 319 {
 320   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(t,ONE_TIME);
 321   ret->setMesh(this);
 322   ret->fillFromAnalytic2(nbOfComp,func);
 323   ret->synchronizeTimeWithSupport();
 324   return ret.retn();
 325 }
 326
 327 /*!
 328  * This method builds a field lying on 'this' with 'nbOfComp' components.
 329  * 'func' is a string that is the expression to evaluate.
 330  * The return field will have type specified by 't'. 't' is also used to determine where values of field will be
 331  * evaluate. This method is different than MEDCouplingMesh::fillFromAnalytic, because 'varsOrder' specifies the pos to assign of vars in 'func'.
 332  *
 333  * @param t type of field returned and specifies where the evaluation of func will be done.
 334  * @param nbOfComp number of components of returned field.
 335  * @param func expression.
 336  * @return field with counter = 1.
 337  */
 338 MEDCouplingFieldDouble *MEDCouplingMesh::fillFromAnalytic3(TypeOfField t, int nbOfComp, const std::vector<std::string>& varsOrder, const char *func) const
 339 {
 340   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(t,ONE_TIME);
 341   ret->setMesh(this);
 342   ret->fillFromAnalytic3(nbOfComp,varsOrder,func);
 343   ret->synchronizeTimeWithSupport();
 344   return ret.retn();
 345 }
 346
 347 /*!
 348  * retruns a newly created mesh with counter=1
 349  * that is the union of \b mesh1 and \b mesh2 if possible. The cells of \b mesh2 will appear after cells of \b mesh1. Idem for nodes.
 350  * The only contraint is that \b mesh1 an \b mesh2 have the same mesh types. If it is not the case please use the other API of MEDCouplingMesh::MergeMeshes,
 351  * with input vector of meshes.
 352  */
 353 MEDCouplingMesh *MEDCouplingMesh::MergeMeshes(const MEDCouplingMesh *mesh1, const MEDCouplingMesh *mesh2) throw(INTERP_KERNEL::Exception)
 354 {
 355   if(!mesh1)
 356     throw INTERP_KERNEL::Exception("MEDCouplingMesh::MergeMeshes : first parameter is an empty mesh !");
 357   if(!mesh2)
 358     throw INTERP_KERNEL::Exception("MEDCouplingMesh::MergeMeshes : second parameter is an empty mesh !");
 359   return mesh1->mergeMyselfWith(mesh2);
 360 }
 361
 362 /*!
 363  * retruns a newly created mesh with counter=1
 364  * that is the union of meshes if possible. The cells of \b meshes[1] will appear after cells of \b meshes[0]. Idem for nodes.
 365  * This method performs a systematic conversion to unstructured meshes before performing aggregation contrary to the other ParaMEDMEM::MEDCouplingMesh::MergeMeshes with
 366  * two parameters that work only on the same type of meshes. So here it is possible to mix different type of meshes.
 367  */
 368 MEDCouplingMesh *MEDCouplingMesh::MergeMeshes(std::vector<const MEDCouplingMesh *>& meshes) throw(INTERP_KERNEL::Exception)
 369 {
 370   std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> > ms1(meshes.size());
 371   std::vector< const MEDCouplingUMesh * > ms2(meshes.size());
 372   for(std::size_t i=0;i<meshes.size();i++)
 373     {
 374       if(meshes[i])
 375         {
 376           MEDCouplingUMesh *cur=meshes[i]->buildUnstructured();
 377           ms1[i]=cur;  ms2[i]=cur;
 378         }
 379       else
 380         {
 381           std::ostringstream oss; oss << "MEDCouplingMesh::MergeMeshes(std::vector<const MEDCouplingMesh *>& meshes) : mesh at pos #" << i << " of input vector of size " << meshes.size() << " is empty !";
 382           throw INTERP_KERNEL::Exception(oss.str().c_str());
 383         }
 384     }
 385   return MEDCouplingUMesh::MergeUMeshes(ms2);
 386 }
 387
 388 /*!
 389  * \param [in] type the geometric type for which the dimension is asked.
 390  * \return the dimension associated to the input geometric type \a type.
 391  *
 392  * \throw if type is equal to \c INTERP_KERNEL::NORM_ERROR or to an unexisting geometric type.
 393  */
 394 int MEDCouplingMesh::GetDimensionOfGeometricType(INTERP_KERNEL::NormalizedCellType type) throw(INTERP_KERNEL::Exception)
 395 {
 396   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
 397   return (int) cm.getDimension();
 398 }
 399
 400 /*!
 401  * \param [in] type the geometric type for which the representation is asked.
 402  * \return the string representation corresponding to the input geometric type \a type.
 403  *
 404  * \throw if type is equal to \c INTERP_KERNEL::NORM_ERROR or to an unexisting geometric type.
 405  */
 406 const char *MEDCouplingMesh::GetReprOfGeometricType(INTERP_KERNEL::NormalizedCellType type) throw(INTERP_KERNEL::Exception)
 407 {
 408   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
 409   return cm.getRepr();
 410 }
 411
 412 void MEDCouplingMesh::getCellsContainingPoint(const double *pos, double eps, std::vector<int>& elts) const
 413 {
 414   int ret=getCellContainingPoint(pos,eps);
 415   elts.push_back(ret);
 416 }
 417
 418 void MEDCouplingMesh::getCellsContainingPoints(const double *pos, int nbOfPoints, double eps, std::vector<int>& elts, std::vector<int>& eltsIndex) const
 419 {
 420   eltsIndex.resize(nbOfPoints+1);
 421   eltsIndex[0]=0;
 422   elts.clear();
 423   int spaceDim=getSpaceDimension();
 424   const double *work=pos;
 425   for(int i=0;i<nbOfPoints;i++,work+=spaceDim)
 426     {
 427       int ret=getCellContainingPoint(work,eps);
 428       if(ret>=0)
 429         {
 430           elts.push_back(ret);
 431           eltsIndex[i+1]=eltsIndex[i]+1;
 432         }
 433       else
 434         eltsIndex[i+1]=eltsIndex[i];
 435     }
 436 }
 437
 438 /*!
 439  * This method writes a file in VTK format into file 'fileName'.
 440  * An exception is thrown if the file is not writable.
 441  */
 442 void MEDCouplingMesh::writeVTK(const char *fileName) const throw(INTERP_KERNEL::Exception)
 443 {
 444   std::string cda,pda;
 445   writeVTKAdvanced(fileName,cda,pda);
 446 }
 447
 448 void MEDCouplingMesh::writeVTKAdvanced(const char *fileName, const std::string& cda, const std::string& pda) const throw(INTERP_KERNEL::Exception)
 449 {
 450   std::ofstream ofs(fileName);
 451   ofs << "<VTKFile type=\""  << getVTKDataSetType() << "\" version=\"0.1\" byte_order=\"LittleEndian\">\n";
 452   writeVTKLL(ofs,cda,pda);
 453   ofs << "</VTKFile>\n";
 454   ofs.close();
 455 }