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