src/MEDCoupling/MEDCouplingMesh.cxx

   1 // Copyright (C) 2007-2013  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 /*!
 104  * Checks if \a this and another MEDCouplingMesh are fully equal.
 105  *  \param [in] other - an instance of MEDCouplingMesh to compare with \a this one.
 106  *  \param [in] prec - precision value used to compare node coordinates.
 107  *  \return bool - \c true if the two meshes are equal, \c false else.
 108  */
 109 bool MEDCouplingMesh::isEqual(const MEDCouplingMesh *other, double prec) const throw(INTERP_KERNEL::Exception)
 110 {
 111   std::string tmp;
 112   return isEqualIfNotWhy(other,prec,tmp);
 113 }
 114
 115 /*!
 116  * This method checks geo equivalence between two meshes : \a this and \a other.
 117  * If no exception is thrown \a this and \a other are geometrically equivalent regarding \a levOfCheck level.
 118  * This method is typically used to change the mesh of a field "safely" depending the \a levOfCheck level considered.
 119  *
 120  * In case of success cell \c other[i] is equal to the cell \c this[cellCor[i]].
 121  * In case of success node \c other->getCoords()[i] is equal to the node \c this->getCoords()[nodeCor[i]].
 122  *
 123  * If \a cellCor is null (or Py_None) it means that for all #i cell in \a other is equal to cell # i in \a this.
 124  *
 125  * If \a nodeCor is null (or Py_None) it means that for all #i node in \a other is equal to node # i in \a this.
 126  *
 127  * So null (or Py_None) returned in \a cellCor and/or \a nodeCor means identity array. This is for optimization reason to avoid to build useless arrays
 128  * for some \a levOfCheck (for example 0).
 129  *
 130  * **Warning a not null output does not mean that it is not identity !**
 131  *
 132  * \param [in] other - the mesh to be compared with \a this.
 133  * \param [in] levOfCheck - input that specifies the level of check specified. The possible values are listed below.
 134  * \param [in] prec - input that specifies precision for double float data used for comparison in meshes.
 135  * \param [out] cellCor - output array not always informed (depending \a levOfCheck param) that gives the corresponding array for cells from \a other to \a this.
 136  * \param [out] nodeCor - output array not always informed (depending \a levOfCheck param) that gives the corresponding array for nodes from \a other to \a this.
 137  *
 138  * Possible values for levOfCheck :
 139  *   - 0 for strict equality. This is the strongest level. \a cellCor and \a nodeCor params are never informed.
 140  *   - 10,11,12 (10+x) for less strict equality. Two meshes are compared geometrically. In case of success \a cellCor and \a 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)
 141  *   - 20,21,22 (20+x), 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 \a 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)
 142  *   - 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)
 143  *   - 2 for deep 'equality' as 0 option except that no control is done on all strings in mesh.
 144  *
 145  * So the most strict level of check is 0 (equality). The least strict is 12. If the level of check 12 throws, the 2 meshes \a this and \a other are not similar enough
 146  * to be compared. An interpolation using MEDCouplingRemapper class should be then used.
 147  */
 148 void MEDCouplingMesh::checkGeoEquivalWith(const MEDCouplingMesh *other, int levOfCheck, double prec,
 149                                           DataArrayInt *&cellCor, DataArrayInt *&nodeCor) const throw(INTERP_KERNEL::Exception)
 150 {
 151   cellCor=0;
 152   nodeCor=0;
 153   if(this==other)
 154     return ;
 155   switch(levOfCheck)
 156     {
 157     case 0:
 158       {
 159         if(!isEqual(other,prec))
 160           throw INTERP_KERNEL::Exception("checkGeoFitWith : Meshes are not equal !");
 161         return ;
 162       }
 163     case 10:
 164     case 11:
 165     case 12:
 166       {
 167         checkDeepEquivalWith(other,levOfCheck-10,prec,cellCor,nodeCor);
 168         return ;
 169       }
 170     case 20:
 171     case 21:
 172     case 22:
 173       {
 174         checkDeepEquivalOnSameNodesWith(other,levOfCheck-20,prec,cellCor);
 175         return ;
 176       }
 177     case 1:
 178       {
 179         checkFastEquivalWith(other,prec);
 180         return;
 181       }
 182     case 2:
 183       {
 184         if(!isEqualWithoutConsideringStr(other,prec))
 185           throw INTERP_KERNEL::Exception("checkGeoFitWith : Meshes are not equal without considering strings !");
 186         return ;
 187       }
 188     default:
 189       throw INTERP_KERNEL::Exception("checkGeoFitWith : Invalid levOfCheck specified ! Value must be in 0,1,2,10,11 or 12.");
 190     }
 191 }
 192
 193 /*!
 194  * Finds cells whose all nodes are in a given array of node ids.
 195  *  \param [in] partBg - the array of node ids.
 196  *  \param [in] partEnd - end of \a partBg, i.e. a pointer to a (last+1)-th element
 197  *          of \a partBg.
 198  *  \return DataArrayInt * - a new instance of DataArrayInt holding ids of found
 199  *          cells. The caller is to delete this array using decrRef() as it is no
 200  *          more needed.
 201  */
 202 DataArrayInt *MEDCouplingMesh::getCellIdsFullyIncludedInNodeIds(const int *partBg, const int *partEnd) const
 203 {
 204   std::vector<int> crest;
 205   std::set<int> p(partBg,partEnd);
 206   int nbOfCells=getNumberOfCells();
 207   for(int i=0;i<nbOfCells;i++)
 208     {
 209       std::vector<int> conn;
 210       getNodeIdsOfCell(i,conn);
 211       bool cont=true;
 212       for(std::vector<int>::const_iterator iter=conn.begin();iter!=conn.end() && cont;iter++)
 213         if(p.find(*iter)==p.end())
 214           cont=false;
 215       if(cont)
 216         crest.push_back(i);
 217     }
 218   DataArrayInt *ret=DataArrayInt::New();
 219   ret->alloc((int)crest.size(),1);
 220   std::copy(crest.begin(),crest.end(),ret->getPointer());
 221   return ret;
 222 }
 223
 224 /*!
 225  * This method checks fastly that \a this and \a other are equal. All common checks are done here.
 226  */
 227 void MEDCouplingMesh::checkFastEquivalWith(const MEDCouplingMesh *other, double prec) const throw(INTERP_KERNEL::Exception)
 228 {
 229   if(getMeshDimension()!=other->getMeshDimension())
 230     throw INTERP_KERNEL::Exception("checkFastEquivalWith : Mesh dimensions are not equal !");
 231   if(getSpaceDimension()!=other->getSpaceDimension())
 232     throw INTERP_KERNEL::Exception("checkFastEquivalWith : Space dimensions are not equal !");
 233   if(getNumberOfCells()!=other->getNumberOfCells())
 234     throw INTERP_KERNEL::Exception("checkFastEquivalWith : number of cells are not equal !");
 235 }
 236
 237 /*!
 238  * This method is very poor and looks only if \a this and \a other are candidate for merge of fields lying repectively on them.
 239  */
 240 bool MEDCouplingMesh::areCompatibleForMerge(const MEDCouplingMesh *other) const
 241 {
 242   if(getMeshDimension()!=other->getMeshDimension())
 243     return false;
 244   if(getSpaceDimension()!=other->getSpaceDimension())
 245     return false;
 246   return true;
 247 }
 248
 249 /*!
 250  * This method is equivalent to MEDCouplingMesh::buildPart method except that here the cell ids are specified using slice \a beginCellIds \a endCellIds and \a stepCellIds.
 251  *
 252  * \sa MEDCouplingMesh::buildPart
 253  */
 254 MEDCouplingMesh *MEDCouplingMesh::buildPartRange(int beginCellIds, int endCellIds, int stepCellIds) const throw(INTERP_KERNEL::Exception)
 255 {
 256   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> cellIds=DataArrayInt::Range(beginCellIds,endCellIds,stepCellIds);
 257   return buildPart(cellIds->begin(),cellIds->end());
 258 }
 259
 260 /*!
 261  * This method is equivalent to MEDCouplingMesh::buildPartAndReduceNodes method except that here the cell ids are specified using slice \a beginCellIds \a endCellIds and \a stepCellIds.
 262  *
 263  * \sa MEDCouplingMesh::buildPartAndReduceNodes
 264  */
 265 MEDCouplingMesh *MEDCouplingMesh::buildPartRangeAndReduceNodes(int beginCellIds, int endCellIds, int stepCellIds, int& beginOut, int& endOut, int& stepOut, DataArrayInt*& arr) const throw(INTERP_KERNEL::Exception)
 266 {
 267   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> cellIds=DataArrayInt::Range(beginCellIds,endCellIds,stepCellIds);
 268   return buildPartAndReduceNodes(cellIds->begin(),cellIds->end(),arr);
 269 }
 270
 271 /*!
 272  * This method builds a field lying on \a this with 'nbOfComp' components.
 273  * 'func' is a pointer that points to a function that takes 2 arrays in parameter and returns a boolean.
 274  * The first array is a in-param of size this->getSpaceDimension and the second an out param of size 'nbOfComp'.
 275  * The return field will have type specified by 't'. 't' is also used to determine where values of field will be
 276  * evaluate.
 277  * Contrary to other fillFromAnalytic methods this method requests a C++ function pointer as input.
 278  * The 'func' is a callback that takes as first parameter an input array of size 'this->getSpaceDimension()',
 279  * the second parameter is a pointer on a valid zone of size at least equal to 'nbOfComp' values. And too finish
 280  * the returned value is a boolean that is equal to False in case of invalid evaluation (log(0) for example...)
 281  *
 282  * \param t type of field returned and specifies where the evaluation of func will be done.
 283  * \param nbOfComp number of components of returned field.
 284  * \param func pointer to a function that should return false if the evaluation failed. (division by 0. for example)
 285  * \return field with counter = 1.
 286  */
 287 MEDCouplingFieldDouble *MEDCouplingMesh::fillFromAnalytic(TypeOfField t, int nbOfComp, FunctionToEvaluate func) const
 288 {
 289   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(t,ONE_TIME);
 290   ret->setMesh(this);
 291   ret->fillFromAnalytic(nbOfComp,func);
 292   ret->synchronizeTimeWithSupport();
 293   return ret.retn();
 294 }
 295
 296 /*!
 297  * This method copyies all tiny strings from other (name and components name).
 298  * @throw if other and this have not same mesh type.
 299  */
 300 void MEDCouplingMesh::copyTinyStringsFrom(const MEDCouplingMesh *other) throw(INTERP_KERNEL::Exception)
 301 {
 302   _name=other->_name;
 303   _description=other->_description;
 304   _time_unit=other->_time_unit;
 305 }
 306
 307 /*!
 308  * This method copies all attributes that are \b NOT arrays in this.
 309  * All tiny attributes not usefully for state of \a this are ignored.
 310  */
 311 void MEDCouplingMesh::copyTinyInfoFrom(const MEDCouplingMesh *other) throw(INTERP_KERNEL::Exception)
 312 {
 313   copyTinyStringsFrom(other);
 314   _time=other->_time;
 315   _iteration=other->_iteration;
 316   _order=other->_order;
 317 }
 318
 319 /*!
 320  * \anchor mcmesh_fillFromAnalytic
 321  * Creates a new MEDCouplingFieldDouble of a given type, one time, with given number of
 322  * components, lying on \a this mesh, with contents got by applying a specified
 323  * function to coordinates of field location points (defined by the given field type).
 324  * For example, if \a t == ParaMEDMEM::ON_CELLS, the function is applied to cell
 325  * barycenters.<br>
 326  * For more info on supported expressions that can be used in the function, see \ref
 327  * MEDCouplingArrayApplyFuncExpr. The function can include arbitrary named variables
 328  * (e.g. "x","y" or "va44") to refer to components of point coordinates. Names of
 329  * variables are sorted in \b alphabetical \b order to associate a variable name with a
 330  * component. For example, in the expression "2*x+z", "x" stands for the component #0
 331  * and "z" stands for the component #1 (\b not #2)!<br>
 332  * In a general case, a value resulting from the function evaluation is assigned to all
 333  * components of the field. But there is a possibility to have its own expression for
 334  * each component within one function. For this purpose, there are predefined variable
 335  * names (IVec, JVec, KVec, LVec etc) each dedicated to a certain component (IVec, to
 336  * the component #0 etc). A factor of such a variable is added to the
 337  * corresponding component only.<br>
 338  * For example, \a nbOfComp == 4, \a this->getSpaceDimension() == 3, coordinates of a
 339  * point are (1.,3.,7.), then
 340  *   - "2*x + z"               produces (5.,5.,5.,5.)
 341  *   - "2*x + 0*y + z"         produces (9.,9.,9.,9.)
 342  *   - "2*x*IVec + (x+z)*LVec" produces (2.,0.,0.,4.)
 343  *   - "2*y*IVec + z*KVec + x" produces (7.,1.,1.,4.)
 344  *
 345  *  \param [in] t - the field type. It defines, apart from other things, points to
 346  *         coordinates of which the function is applied to get field values.
 347  *  \param [in] nbOfComp - the number of components in the result field.
 348  *  \param [in] func - a string defining the expression which is evaluated to get
 349  *         field values.
 350  *  \return MEDCouplingFieldDouble * - a new instance of MEDCouplingFieldDouble. The
 351  *         caller is to delete this field using decrRef() as it is no more needed.
 352  *  \throw If the nodal connectivity of cells is not defined.
 353  *  \throw If computing \a func fails.
 354  *
 355  *  \ref cpp_mcmesh_fillFromAnalytic "Here is a C++ example".<br>
 356  *  \ref  py_mcmesh_fillFromAnalytic "Here is a Python example".
 357  */
 358 MEDCouplingFieldDouble *MEDCouplingMesh::fillFromAnalytic(TypeOfField t, int nbOfComp, const char *func) const
 359 {
 360   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(t,ONE_TIME);
 361   ret->setMesh(this);
 362   ret->fillFromAnalytic(nbOfComp,func);
 363   ret->synchronizeTimeWithSupport();
 364   return ret.retn();
 365 }
 366
 367 /*!
 368  * Creates a new MEDCouplingFieldDouble of a given type, one time, with given number of
 369  * components, lying on \a this mesh, with contents got by applying a specified
 370  * function to coordinates of field location points (defined by the given field type).
 371  * For example, if \a t == ParaMEDMEM::ON_CELLS, the function is applied to cell
 372  * barycenters. This method differs from
 373  * \ref MEDCouplingMesh::fillFromAnalytic(TypeOfField t, int nbOfComp, const char *func) const "fillFromAnalytic()"
 374  * by the way how variable
 375  * names, used in the function, are associated with components of coordinates of field
 376  * location points; here, a variable name corresponding to a component is retrieved from
 377  * a corresponding node coordinates array (where it is set via
 378  * DataArrayDouble::setInfoOnComponent()).<br>
 379  * For more info on supported expressions that can be used in the function, see \ref
 380  * MEDCouplingArrayApplyFuncExpr. <br>
 381  * In a general case, a value resulting from the function evaluation is assigned to all
 382  * components of a field value. But there is a possibility to have its own expression for
 383  * each component within one function. For this purpose, there are predefined variable
 384  * names (IVec, JVec, KVec, LVec etc) each dedicated to a certain component (IVec, to
 385  * the component #0 etc). A factor of such a variable is added to the
 386  * corresponding component only.<br>
 387  * For example, \a nbOfComp == 4, \a this->getSpaceDimension() == 3, names of
 388  * spatial components are "x", "y" and "z", coordinates of a
 389  * point are (1.,3.,7.), then
 390  *   - "2*x + z"               produces (9.,9.,9.,9.)
 391  *   - "2*x*IVec + (x+z)*LVec" produces (2.,0.,0.,8.)
 392  *   - "2*y*IVec + z*KVec + x" produces (7.,1.,1.,8.)
 393  *
 394  *  \param [in] t - the field type. It defines, apart from other things, the points to
 395  *         coordinates of which the function is applied to get field values.
 396  *  \param [in] nbOfComp - the number of components in the result field.
 397  *  \param [in] func - a string defining the expression which is evaluated to get
 398  *         field values.
 399  *  \return MEDCouplingFieldDouble * - a new instance of MEDCouplingFieldDouble. The
 400  *         caller is to delete this field using decrRef() as it is no more needed.
 401  *  \throw If the node coordinates are not defined.
 402  *  \throw If the nodal connectivity of cells is not defined.
 403  *  \throw If computing \a func fails.
 404  *
 405  *  \ref cpp_mcmesh_fillFromAnalytic2 "Here is a C++ example".<br>
 406  *  \ref  py_mcmesh_fillFromAnalytic2 "Here is a Python example".
 407  */
 408 MEDCouplingFieldDouble *MEDCouplingMesh::fillFromAnalytic2(TypeOfField t, int nbOfComp, const char *func) const
 409 {
 410   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(t,ONE_TIME);
 411   ret->setMesh(this);
 412   ret->fillFromAnalytic2(nbOfComp,func);
 413   ret->synchronizeTimeWithSupport();
 414   return ret.retn();
 415 }
 416
 417 /*!
 418  * Creates a new MEDCouplingFieldDouble of a given type, one time, with given number of
 419  * components, lying on \a this mesh, with contents got by applying a specified
 420  * function to coordinates of field location points (defined by the given field type).
 421  * For example, if \a t == ParaMEDMEM::ON_CELLS, the function is applied to cell
 422  * barycenters. This method differs from \ref  \ref mcmesh_fillFromAnalytic
 423  * "fillFromAnalytic()" by the way how variable
 424  * names, used in the function, are associated with components of coordinates of field
 425  * location points; here, a component index of a variable is defined by a
 426  * rank of the variable within the input array \a varsOrder.<br>
 427  * For more info on supported expressions that can be used in the function, see \ref
 428  * MEDCouplingArrayApplyFuncExpr.
 429  * In a general case, a value resulting from the function evaluation is assigned to all
 430  * components of the field. But there is a possibility to have its own expression for
 431  * each component within one function. For this purpose, there are predefined variable
 432  * names (IVec, JVec, KVec, LVec etc) each dedicated to a certain component (IVec, to
 433  * the component #0 etc). A factor of such a variable is added to the
 434  * corresponding component only.<br>
 435  * For example, \a nbOfComp == 4, \a this->getSpaceDimension() == 3, names of
 436  * spatial components are given in \a varsOrder: ["x", "y","z"], coordinates of a
 437  * point are (1.,3.,7.), then
 438  *   - "2*x + z"               produces (9.,9.,9.,9.)
 439  *   - "2*x*IVec + (x+z)*LVec" produces (2.,0.,0.,8.)
 440  *   - "2*y*IVec + z*KVec + x" produces (7.,1.,1.,8.)
 441  *
 442  *  \param [in] t - the field type. It defines, apart from other things, the points to
 443  *         coordinates of which the function is applied to get field values.
 444  *  \param [in] nbOfComp - the number of components in the result field.
 445  *  \param [in] varsOrder - the vector defining names of variables used to refer to
 446  *         components of coordinates of field location points. A variable named
 447  *         varsOrder[0] refers to the component #0 etc.
 448  *  \param [in] func - a string defining the expression which is evaluated to get
 449  *         field values.
 450  *  \return MEDCouplingFieldDouble * - a new instance of MEDCouplingFieldDouble. The
 451  *         caller is to delete this field using decrRef() as it is no more needed.
 452  *  \throw If the node coordinates are not defined.
 453  *  \throw If the nodal connectivity of cells is not defined.
 454  *  \throw If computing \a func fails.
 455  *
 456  *  \ref cpp_mcmesh_fillFromAnalytic3 "Here is a C++ example".<br>
 457  *  \ref  py_mcmesh_fillFromAnalytic3 "Here is a Python example".
 458  */
 459 MEDCouplingFieldDouble *MEDCouplingMesh::fillFromAnalytic3(TypeOfField t, int nbOfComp, const std::vector<std::string>& varsOrder, const char *func) const
 460 {
 461   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(t,ONE_TIME);
 462   ret->setMesh(this);
 463   ret->fillFromAnalytic3(nbOfComp,varsOrder,func);
 464   ret->synchronizeTimeWithSupport();
 465   return ret.retn();
 466 }
 467
 468 /*!
 469  * Creates a new MEDCouplingMesh by concatenating two given meshes, if possible.
 470  * Cells and nodes of
 471  * the first mesh precede cells and nodes of the second mesh within the result mesh.
 472  * The meshes must be of the same mesh type, else, an exception is thrown. The method
 473  * MergeMeshes(), accepting a vector of input meshes, has no such a limitation.
 474  *  \param [in] mesh1 - the first mesh.
 475  *  \param [in] mesh2 - the second mesh.
 476  *  \return MEDCouplingMesh * - the result mesh. It is a new instance of
 477  *          MEDCouplingMesh. The caller is to delete this mesh using decrRef() as it
 478  *          is no more needed.
 479  *  \throw If the meshes are of different mesh type.
 480  */
 481 MEDCouplingMesh *MEDCouplingMesh::MergeMeshes(const MEDCouplingMesh *mesh1, const MEDCouplingMesh *mesh2) throw(INTERP_KERNEL::Exception)
 482 {
 483   if(!mesh1)
 484     throw INTERP_KERNEL::Exception("MEDCouplingMesh::MergeMeshes : first parameter is an empty mesh !");
 485   if(!mesh2)
 486     throw INTERP_KERNEL::Exception("MEDCouplingMesh::MergeMeshes : second parameter is an empty mesh !");
 487   return mesh1->mergeMyselfWith(mesh2);
 488 }
 489
 490 /*!
 491  * Creates a new MEDCouplingMesh by concatenating all given meshes, if possible.
 492  * Cells and nodes of
 493  * the *i*-th mesh precede cells and nodes of the (*i*+1)-th mesh within the result mesh.
 494  * This method performs a systematic conversion to unstructured meshes before
 495  * performing aggregation contrary to the other MergeMeshes()
 496  * with two parameters that works only on the same type of meshes. So here it is possible
 497  * to mix different type of meshes.
 498  *  \param [in] meshes - a vector of meshes to concatenate.
 499  *  \return MEDCouplingMesh * - the result mesh. It is a new instance of
 500  *          MEDCouplingUMesh. The caller is to delete this mesh using decrRef() as it
 501  *          is no more needed.
 502  *  \throw If \a meshes.size() == 0.
 503  *  \throw If \a size[ *i* ] == NULL.
 504  *  \throw If the coordinates is not set in none of the meshes.
 505  *  \throw If \a meshes[ *i* ]->getMeshDimension() < 0.
 506  *  \throw If the \a meshes are of different dimension (getMeshDimension()).
 507  */
 508 MEDCouplingMesh *MEDCouplingMesh::MergeMeshes(std::vector<const MEDCouplingMesh *>& meshes) throw(INTERP_KERNEL::Exception)
 509 {
 510   std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> > ms1(meshes.size());
 511   std::vector< const MEDCouplingUMesh * > ms2(meshes.size());
 512   for(std::size_t i=0;i<meshes.size();i++)
 513     {
 514       if(meshes[i])
 515         {
 516           MEDCouplingUMesh *cur=meshes[i]->buildUnstructured();
 517           ms1[i]=cur;  ms2[i]=cur;
 518         }
 519       else
 520         {
 521           std::ostringstream oss; oss << "MEDCouplingMesh::MergeMeshes(std::vector<const MEDCouplingMesh *>& meshes) : mesh at pos #" << i << " of input vector of size " << meshes.size() << " is empty !";
 522           throw INTERP_KERNEL::Exception(oss.str().c_str());
 523         }
 524     }
 525   return MEDCouplingUMesh::MergeUMeshes(ms2);
 526 }
 527
 528 /*!
 529  * For example if \a type is INTERP_KERNEL::NORM_TRI3 , INTERP_KERNEL::NORM_POLYGON is returned.
 530  * If \a type is INTERP_KERNEL::NORM_HEXA8 , INTERP_KERNEL::NORM_POLYHED is returned.
 531  *
 532  * \param [in] type the geometric type for which the corresponding dynamic type, is asked.
 533  * \return the corresponding dynamic type, able to store the input \a type.
 534  *
 535  * \throw if type is equal to \c INTERP_KERNEL::NORM_ERROR or to an unexisting geometric type.
 536  */
 537 INTERP_KERNEL::NormalizedCellType MEDCouplingMesh::GetCorrespondingPolyType(INTERP_KERNEL::NormalizedCellType type) throw(INTERP_KERNEL::Exception)
 538 {
 539   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
 540   return cm.getCorrespondingPolyType();
 541 }
 542
 543 /*!
 544  * \param [in] type the geometric type for which the number of nodes consituting it, is asked.
 545  * \return number of nodes consituting the input geometric type \a type.
 546  *
 547  * \throw if type is dynamic as \c INTERP_KERNEL::NORM_POLYHED , \c INTERP_KERNEL::NORM_POLYGON , \c INTERP_KERNEL::NORM_QPOLYG
 548  * \throw if type is equal to \c INTERP_KERNEL::NORM_ERROR or to an unexisting geometric type.
 549  */
 550 int MEDCouplingMesh::GetNumberOfNodesOfGeometricType(INTERP_KERNEL::NormalizedCellType type) throw(INTERP_KERNEL::Exception)
 551 {
 552   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
 553   if(cm.isDynamic())
 554     throw INTERP_KERNEL::Exception("MEDCouplingMesh::GetNumberOfNodesOfGeometricType : the input geometric type is dynamic ! Impossible to return a fixed number of nodes constituting it !");
 555   return (int) cm.getNumberOfNodes();
 556 }
 557
 558 /*!
 559  * \param [in] type the geometric type for which the status static/dynamic is asked.
 560  * \return true for static geometric type, false for dynamic geometric type.
 561  *
 562  * \throw if type is equal to \c INTERP_KERNEL::NORM_ERROR or to an unexisting geometric type.
 563  */
 564 bool MEDCouplingMesh::IsStaticGeometricType(INTERP_KERNEL::NormalizedCellType type) throw(INTERP_KERNEL::Exception)
 565 {
 566   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
 567   return !cm.isDynamic();
 568 }
 569
 570 bool MEDCouplingMesh::IsLinearGeometricType(INTERP_KERNEL::NormalizedCellType type) throw(INTERP_KERNEL::Exception)
 571 {
 572   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
 573   return !cm.isQuadratic();
 574 }
 575
 576 /*!
 577  * \param [in] type the geometric type for which the dimension is asked.
 578  * \return the dimension associated to the input geometric type \a type.
 579  *
 580  * \throw if type is equal to \c INTERP_KERNEL::NORM_ERROR or to an unexisting geometric type.
 581  */
 582 int MEDCouplingMesh::GetDimensionOfGeometricType(INTERP_KERNEL::NormalizedCellType type) throw(INTERP_KERNEL::Exception)
 583 {
 584   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
 585   return (int) cm.getDimension();
 586 }
 587
 588 /*!
 589  * \param [in] type the geometric type for which the representation is asked.
 590  * \return the string representation corresponding to the input geometric type \a type.
 591  *
 592  * \throw if type is equal to \c INTERP_KERNEL::NORM_ERROR or to an unexisting geometric type.
 593  */
 594 const char *MEDCouplingMesh::GetReprOfGeometricType(INTERP_KERNEL::NormalizedCellType type) throw(INTERP_KERNEL::Exception)
 595 {
 596   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
 597   return cm.getRepr();
 598 }
 599
 600 /*!
 601  * Finds cells in contact with a ball (i.e. a point with precision).
 602  * \warning This method is suitable if the caller intends to evaluate only one
 603  *          point, for more points getCellsContainingPoints() is recommended as it is
 604  *          faster.
 605  *  \param [in] pos - array of coordinates of the ball central point.
 606  *  \param [in] eps - ball radius.
 607  *  \param [in,out] elts - vector returning ids of the found cells. It is cleared
 608  *         before inserting ids.
 609  *
 610  *  \ref cpp_mcumesh_getCellsContainingPoint "Here is a C++ example".<br>
 611  *  \ref  py_mcumesh_getCellsContainingPoint "Here is a Python example".
 612  */
 613 void MEDCouplingMesh::getCellsContainingPoint(const double *pos, double eps, std::vector<int>& elts) const
 614 {
 615   int ret=getCellContainingPoint(pos,eps);
 616   elts.push_back(ret);
 617 }
 618
 619 /*!
 620  * Finds cells in contact with several balls (i.e. points with precision).
 621  * This method is an extension of getCellContainingPoint() and
 622  * getCellsContainingPoint() for the case of multiple points.
 623  *  \param [in] pos - an array of coordinates of points in full interlace mode :
 624  *         X0,Y0,Z0,X1,Y1,Z1,... Size of the array must be \a
 625  *         this->getSpaceDimension() * \a nbOfPoints
 626  *  \param [in] nbOfPoints - number of points to locate within \a this mesh.
 627  *  \param [in] eps - radius of balls (i.e. the precision).
 628  *  \param [in,out] elts - vector returning ids of found cells.
 629  *  \param [in,out] eltsIndex - an array, of length \a nbOfPoints + 1,
 630  *         dividing cell ids in \a elts into groups each referring to one
 631  *         point. Its every element (except the last one) is an index pointing to the
 632  *         first id of a group of cells. For example cells in contact with the *i*-th
 633  *         point are described by following range of indices:
 634  *         [ \a eltsIndex[ *i* ], \a eltsIndex[ *i*+1 ] ) and the cell ids are
 635  *         \a elts[ \a eltsIndex[ *i* ]], \a elts[ \a eltsIndex[ *i* ] + 1 ], ...
 636  *         Number of cells in contact with the *i*-th point is
 637  *         \a eltsIndex[ *i*+1 ] - \a eltsIndex[ *i* ].
 638  *
 639  *  \ref cpp_mcumesh_getCellsContainingPoints "Here is a C++ example".<br>
 640  *  \ref  py_mcumesh_getCellsContainingPoints "Here is a Python example".
 641  */
 642 void MEDCouplingMesh::getCellsContainingPoints(const double *pos, int nbOfPoints, double eps, std::vector<int>& elts, std::vector<int>& eltsIndex) const
 643 {
 644   eltsIndex.resize(nbOfPoints+1);
 645   eltsIndex[0]=0;
 646   elts.clear();
 647   int spaceDim=getSpaceDimension();
 648   const double *work=pos;
 649   for(int i=0;i<nbOfPoints;i++,work+=spaceDim)
 650     {
 651       int ret=getCellContainingPoint(work,eps);
 652       if(ret>=0)
 653         {
 654           elts.push_back(ret);
 655           eltsIndex[i+1]=eltsIndex[i]+1;
 656         }
 657       else
 658         eltsIndex[i+1]=eltsIndex[i];
 659     }
 660 }
 661
 662 /*!
 663  * Writes \a this mesh into a VTK format file named as specified.
 664  *  \param [in] fileName - the name of the file to write in.
 665  *  \throw If \a fileName is not a writable file.
 666  */
 667 void MEDCouplingMesh::writeVTK(const char *fileName) const throw(INTERP_KERNEL::Exception)
 668 {
 669   std::string cda,pda;
 670   writeVTKAdvanced(fileName,cda,pda);
 671 }
 672
 673 void MEDCouplingMesh::writeVTKAdvanced(const char *fileName, const std::string& cda, const std::string& pda) const throw(INTERP_KERNEL::Exception)
 674 {
 675   std::ofstream ofs(fileName);
 676   ofs << "<VTKFile type=\""  << getVTKDataSetType() << "\" version=\"0.1\" byte_order=\"LittleEndian\">\n";
 677   writeVTKLL(ofs,cda,pda);
 678   ofs << "</VTKFile>\n";
 679   ofs.close();
 680 }