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