1 // Copyright (C) 2007-2023 CEA/DEN, EDF R&D
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, or (at your option) any later version.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // Author : Anthony Geay (CEA/DEN)
21 #include "MEDCouplingMesh.hxx"
22 #include "MEDCouplingUMesh.hxx"
23 #include "MEDCouplingMemArray.txx"
24 #include "MEDCouplingFieldDouble.hxx"
25 #include "MEDCouplingFieldDiscretization.hxx"
34 using namespace MEDCoupling;
36 MEDCouplingMesh::MEDCouplingMesh():_time(0.),_iteration(-1),_order(-1)
40 MEDCouplingMesh::MEDCouplingMesh(const MEDCouplingMesh& other):RefCountObject(other),_name(other._name),_description(other._description),
41 _time(other._time),_iteration(other._iteration),
42 _order(other._order),_time_unit(other._time_unit)
46 std::size_t MEDCouplingMesh::getHeapMemorySizeWithoutChildren() const
48 return _name.capacity()+_description.capacity()+_time_unit.capacity();
52 * This method is only for ParaMEDMEM in ParaFIELD constructor.
54 bool MEDCouplingMesh::isStructured() const
56 return getType()==CARTESIAN;
59 bool MEDCouplingMesh::isEqualIfNotWhy(const MEDCouplingMesh *other, double prec, std::string& reason) const
62 throw INTERP_KERNEL::Exception("MEDCouplingMesh::isEqualIfNotWhy : other instance is NULL !");
63 std::ostringstream oss; oss.precision(15);
64 if(_name!=other->_name)
66 oss << "Mesh names differ : this name = \"" << _name << "\" and other name = \"" << other->_name << "\" !";
70 if(_description!=other->_description)
72 oss << "Mesh descriptions differ : this description = \"" << _description << "\" and other description = \"" << other->_description << "\" !";
76 if(_iteration!=other->_iteration)
78 oss << "Mesh iterations differ : this iteration = \"" << _iteration << "\" and other iteration = \"" << other->_iteration << "\" !";
82 if(_order!=other->_order)
84 oss << "Mesh orders differ : this order = \"" << _order << "\" and other order = \"" << other->_order << "\" !";
88 if(_time_unit!=other->_time_unit)
90 oss << "Mesh time units differ : this time unit = \"" << _time_unit << "\" and other time unit = \"" << other->_time_unit << "\" !";
94 if(fabs(_time-other->_time)>=1e-12)
96 oss << "Mesh times differ : this time = \"" << _time << "\" and other time = \"" << other->_time << "\" !";
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.
109 bool MEDCouplingMesh::isEqual(const MEDCouplingMesh *other, double prec) const
112 return isEqualIfNotWhy(other,prec,tmp);
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.
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]].
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.
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.
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).
130 * **Warning a not null output does not mean that it is not identity !**
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.
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.
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.
148 void MEDCouplingMesh::checkGeoEquivalWith(const MEDCouplingMesh *other, int levOfCheck, double prec,
149 DataArrayIdType *&cellCor, DataArrayIdType *&nodeCor) const
159 if(!isEqual(other,prec))
160 throw INTERP_KERNEL::Exception("checkGeoFitWith : Meshes are not equal !");
167 checkDeepEquivalWith(other,levOfCheck-10,prec,cellCor,nodeCor);
174 checkDeepEquivalOnSameNodesWith(other,levOfCheck-20,prec,cellCor);
179 checkFastEquivalWith(other,prec);
184 if(!isEqualWithoutConsideringStr(other,prec))
185 throw INTERP_KERNEL::Exception("checkGeoFitWith : Meshes are not equal without considering strings !");
189 throw INTERP_KERNEL::Exception("checkGeoFitWith : Invalid levOfCheck specified ! Value must be in 0,1,2,10,11 or 12.");
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
198 * \return DataArrayIdType * - a new instance of DataArrayIdType holding ids of found
199 * cells. The caller is to delete this array using decrRef() as it is no
202 DataArrayIdType *MEDCouplingMesh::getCellIdsFullyIncludedInNodeIds(const mcIdType *partBg, const mcIdType *partEnd) const
204 std::vector<mcIdType> crest;
205 std::set<mcIdType> p(partBg,partEnd);
206 mcIdType nbOfCells=getNumberOfCells();
207 for(mcIdType i=0;i<nbOfCells;i++)
209 std::vector<mcIdType> conn;
210 getNodeIdsOfCell(i,conn);
212 for(std::vector<mcIdType>::const_iterator iter=conn.begin();iter!=conn.end() && cont;iter++)
213 if(p.find(*iter)==p.end())
218 DataArrayIdType *ret=DataArrayIdType::New();
219 ret->alloc(crest.size(),1);
220 std::copy(crest.begin(),crest.end(),ret->getPointer());
225 * This method checks fastly that \a this and \a other are equal. All common checks are done here.
227 void MEDCouplingMesh::checkFastEquivalWith(const MEDCouplingMesh *other, double prec) const
230 throw INTERP_KERNEL::Exception("MEDCouplingMesh::checkFastEquivalWith : input mesh is null !");
231 if(getMeshDimension()!=other->getMeshDimension())
232 throw INTERP_KERNEL::Exception("checkFastEquivalWith : Mesh dimensions are not equal !");
233 if(getSpaceDimension()!=other->getSpaceDimension())
234 throw INTERP_KERNEL::Exception("checkFastEquivalWith : Space dimensions are not equal !");
235 if(getNumberOfCells()!=other->getNumberOfCells())
236 throw INTERP_KERNEL::Exception("checkFastEquivalWith : number of cells are not equal !");
240 * This method is very poor and looks only if \a this and \a other are candidate for merge of fields lying respectively on them.
242 bool MEDCouplingMesh::areCompatibleForMerge(const MEDCouplingMesh *other) const
245 throw INTERP_KERNEL::Exception("MEDCouplingMesh::areCompatibleForMerge : input mesh is null !");
246 if(getMeshDimension()!=other->getMeshDimension())
248 if(getSpaceDimension()!=other->getSpaceDimension())
254 * 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.
255 * \b WARNING , there is a big difference compared to MEDCouplingMesh::buildPart method.
256 * If the input range is equal all cells in \a this, \a this is returned !
258 * \return a new ref to be managed by the caller. Warning this ref can be equal to \a this if input slice is exactly equal to the whole cells in the same order.
260 * \sa MEDCouplingMesh::buildPart
262 MEDCouplingMesh *MEDCouplingMesh::buildPartRange(mcIdType beginCellIds, mcIdType endCellIds, mcIdType stepCellIds) const
264 if(beginCellIds==0 && endCellIds==getNumberOfCells() && stepCellIds==1)
266 MEDCouplingMesh *ret(const_cast<MEDCouplingMesh *>(this));
272 MCAuto<DataArrayIdType> cellIds=DataArrayIdType::Range(beginCellIds,endCellIds,stepCellIds);
273 return buildPart(cellIds->begin(),cellIds->end());
278 * 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.
280 * \sa MEDCouplingMesh::buildPartAndReduceNodes
282 MEDCouplingMesh *MEDCouplingMesh::buildPartRangeAndReduceNodes(mcIdType beginCellIds, mcIdType endCellIds, mcIdType stepCellIds, mcIdType& beginOut, mcIdType& endOut, mcIdType& stepOut, DataArrayIdType*& arr) const
284 MCAuto<DataArrayIdType> cellIds=DataArrayIdType::Range(beginCellIds,endCellIds,stepCellIds);
285 return buildPartAndReduceNodes(cellIds->begin(),cellIds->end(),arr);
289 * This method builds a field lying on \a this with 'nbOfComp' components.
290 * 'func' is a pointer that points to a function that takes 2 arrays in parameter and returns a boolean.
291 * The first array is a in-param of size this->getSpaceDimension and the second an out param of size 'nbOfComp'.
292 * The return field will have type specified by 't'. 't' is also used to determine where values of field will be
294 * Contrary to other fillFromAnalytic methods this method requests a C++ function pointer as input.
295 * The 'func' is a callback that takes as first parameter an input array of size 'this->getSpaceDimension()',
296 * the second parameter is a pointer on a valid zone of size at least equal to 'nbOfComp' values. And too finish
297 * the returned value is a boolean that is equal to False in case of invalid evaluation (log(0) for example...)
299 * \param t type of field returned and specifies where the evaluation of func will be done.
300 * \param nbOfComp number of components of returned field.
301 * \param func pointer to a function that should return false if the evaluation failed. (division by 0. for example)
302 * \return field with counter = 1.
304 MEDCouplingFieldDouble *MEDCouplingMesh::fillFromAnalytic(TypeOfField t, int nbOfComp, FunctionToEvaluate func) const
306 MCAuto<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(t,ONE_TIME);
308 ret->fillFromAnalytic(nbOfComp,func);
309 ret->synchronizeTimeWithSupport();
314 * This method copyies all tiny strings from other (name and components name).
315 * @throw if other and this have not same mesh type.
317 void MEDCouplingMesh::copyTinyStringsFrom(const MEDCouplingMesh *other)
320 throw INTERP_KERNEL::Exception("MEDCouplingMesh::copyTinyStringsFrom : input mesh is null !");
322 _description=other->_description;
323 _time_unit=other->_time_unit;
327 * This method copies all attributes that are \b NOT arrays in this.
328 * All tiny attributes not usefully for state of \a this are ignored.
330 void MEDCouplingMesh::copyTinyInfoFrom(const MEDCouplingMesh *other)
333 _iteration=other->_iteration;
334 _order=other->_order;
335 copyTinyStringsFrom(other);
339 * \anchor mcmesh_fillFromAnalytic
340 * Creates a new MEDCouplingFieldDouble of a given type, one time, with given number of
341 * components, lying on \a this mesh, with contents got by applying a specified
342 * function to coordinates of field location points (defined by the given field type).
343 * For example, if \a t == MEDCoupling::ON_CELLS, the function is applied to cell
345 * For more info on supported expressions that can be used in the function, see \ref
346 * MEDCouplingArrayApplyFuncExpr. The function can include arbitrary named variables
347 * (e.g. "x","y" or "va44") to refer to components of point coordinates. Names of
348 * variables are sorted in \b alphabetical \b order to associate a variable name with a
349 * component. For example, in the expression "2*x+z", "x" stands for the component #0
350 * and "z" stands for the component #1 (\b not #2)!<br>
351 * In a general case, a value resulting from the function evaluation is assigned to all
352 * components of the field. But there is a possibility to have its own expression for
353 * each component within one function. For this purpose, there are predefined variable
354 * names (IVec, JVec, KVec, LVec etc) each dedicated to a certain component (IVec, to
355 * the component #0 etc). A factor of such a variable is added to the
356 * corresponding component only.<br>
357 * For example, \a nbOfComp == 4, \a this->getSpaceDimension() == 3, coordinates of a
358 * point are (1.,3.,7.), then
359 * - "2*x + z" produces (5.,5.,5.,5.)
360 * - "2*x + 0*y + z" produces (9.,9.,9.,9.)
361 * - "2*x*IVec + (x+z)*LVec" produces (2.,0.,0.,4.)
362 * - "2*y*IVec + z*KVec + x" produces (7.,1.,1.,4.)
364 * \param [in] t - the field type. It defines, apart from other things, points to
365 * coordinates of which the function is applied to get field values.
366 * \param [in] nbOfComp - the number of components in the result field.
367 * \param [in] func - a string defining the expression which is evaluated to get
369 * \return MEDCouplingFieldDouble * - a new instance of MEDCouplingFieldDouble. The
370 * caller is to delete this field using decrRef() as it is no more needed.
371 * \throw If the nodal connectivity of cells is not defined.
372 * \throw If computing \a func fails.
374 * \if ENABLE_EXAMPLES
375 * \ref cpp_mcmesh_fillFromAnalytic "Here is a C++ example".<br>
376 * \ref py_mcmesh_fillFromAnalytic "Here is a Python example".
379 MEDCouplingFieldDouble *MEDCouplingMesh::fillFromAnalytic(TypeOfField t, int nbOfComp, const std::string& func) const
381 MCAuto<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(t,ONE_TIME);
383 ret->fillFromAnalytic(nbOfComp,func);
384 ret->synchronizeTimeWithSupport();
389 * Creates a new MEDCouplingFieldDouble of a given type, one time, with given number of
390 * components, lying on \a this mesh, with contents got by applying a specified
391 * function to coordinates of field location points (defined by the given field type).
392 * For example, if \a t == MEDCoupling::ON_CELLS, the function is applied to cell
393 * barycenters. This method differs from
394 * \ref MEDCouplingMesh::fillFromAnalytic(TypeOfField t, int nbOfComp, const std::string& func) const "fillFromAnalytic()"
395 * by the way how variable
396 * names, used in the function, are associated with components of coordinates of field
397 * location points; here, a variable name corresponding to a component is retrieved from
398 * a corresponding node coordinates array (where it is set via
399 * DataArrayDouble::setInfoOnComponent()).<br>
400 * For more info on supported expressions that can be used in the function, see \ref
401 * MEDCouplingArrayApplyFuncExpr. <br>
402 * In a general case, a value resulting from the function evaluation is assigned to all
403 * components of a field value. But there is a possibility to have its own expression for
404 * each component within one function. For this purpose, there are predefined variable
405 * names (IVec, JVec, KVec, LVec etc) each dedicated to a certain component (IVec, to
406 * the component #0 etc). A factor of such a variable is added to the
407 * corresponding component only.<br>
408 * For example, \a nbOfComp == 4, \a this->getSpaceDimension() == 3, names of
409 * spatial components are "x", "y" and "z", coordinates of a
410 * point are (1.,3.,7.), then
411 * - "2*x + z" produces (9.,9.,9.,9.)
412 * - "2*x*IVec + (x+z)*LVec" produces (2.,0.,0.,8.)
413 * - "2*y*IVec + z*KVec + x" produces (7.,1.,1.,8.)
415 * \param [in] t - the field type. It defines, apart from other things, the points to
416 * coordinates of which the function is applied to get field values.
417 * \param [in] nbOfComp - the number of components in the result field.
418 * \param [in] func - a string defining the expression which is evaluated to get
420 * \return MEDCouplingFieldDouble * - a new instance of MEDCouplingFieldDouble. The
421 * caller is to delete this field using decrRef() as it is no more needed.
422 * \throw If the node coordinates are not defined.
423 * \throw If the nodal connectivity of cells is not defined.
424 * \throw If computing \a func fails.
426 * \if ENABLE_EXAMPLES
427 * \ref cpp_mcmesh_fillFromAnalytic2 "Here is a C++ example".<br>
428 * \ref py_mcmesh_fillFromAnalytic2 "Here is a Python example".
431 MEDCouplingFieldDouble *MEDCouplingMesh::fillFromAnalyticCompo(TypeOfField t, int nbOfComp, const std::string& func) const
433 MCAuto<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(t,ONE_TIME);
435 ret->fillFromAnalyticCompo(nbOfComp,func);
436 ret->synchronizeTimeWithSupport();
441 * Creates a new MEDCouplingFieldDouble of a given type, one time, with given number of
442 * components, lying on \a this mesh, with contents got by applying a specified
443 * function to coordinates of field location points (defined by the given field type).
444 * For example, if \a t == MEDCoupling::ON_CELLS, the function is applied to cell
445 * barycenters. This method differs from \ref \ref mcmesh_fillFromAnalytic
446 * "fillFromAnalytic()" by the way how variable
447 * names, used in the function, are associated with components of coordinates of field
448 * location points; here, a component index of a variable is defined by a
449 * rank of the variable within the input array \a varsOrder.<br>
450 * For more info on supported expressions that can be used in the function, see \ref
451 * MEDCouplingArrayApplyFuncExpr.
452 * In a general case, a value resulting from the function evaluation is assigned to all
453 * components of the field. But there is a possibility to have its own expression for
454 * each component within one function. For this purpose, there are predefined variable
455 * names (IVec, JVec, KVec, LVec etc) each dedicated to a certain component (IVec, to
456 * the component #0 etc). A factor of such a variable is added to the
457 * corresponding component only.<br>
458 * For example, \a nbOfComp == 4, \a this->getSpaceDimension() == 3, names of
459 * spatial components are given in \a varsOrder: ["x", "y","z"], coordinates of a
460 * point are (1.,3.,7.), then
461 * - "2*x + z" produces (9.,9.,9.,9.)
462 * - "2*x*IVec + (x+z)*LVec" produces (2.,0.,0.,8.)
463 * - "2*y*IVec + z*KVec + x" produces (7.,1.,1.,8.)
465 * \param [in] t - the field type. It defines, apart from other things, the points to
466 * coordinates of which the function is applied to get field values.
467 * \param [in] nbOfComp - the number of components in the result field.
468 * \param [in] varsOrder - the vector defining names of variables used to refer to
469 * components of coordinates of field location points. A variable named
470 * varsOrder[0] refers to the component #0 etc.
471 * \param [in] func - a string defining the expression which is evaluated to get
473 * \return MEDCouplingFieldDouble * - a new instance of MEDCouplingFieldDouble. The
474 * caller is to delete this field using decrRef() as it is no more needed.
475 * \throw If the node coordinates are not defined.
476 * \throw If the nodal connectivity of cells is not defined.
477 * \throw If computing \a func fails.
479 * \if ENABLE_EXAMPLES
480 * \ref cpp_mcmesh_fillFromAnalytic3 "Here is a C++ example".<br>
481 * \ref py_mcmesh_fillFromAnalytic3 "Here is a Python example".
484 MEDCouplingFieldDouble *MEDCouplingMesh::fillFromAnalyticNamedCompo(TypeOfField t, int nbOfComp, const std::vector<std::string>& varsOrder, const std::string& func) const
486 MCAuto<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(t,ONE_TIME);
488 ret->fillFromAnalyticNamedCompo(nbOfComp,varsOrder,func);
489 ret->synchronizeTimeWithSupport();
494 * Creates a new MEDCouplingMesh by concatenating two given meshes, if possible.
496 * the first mesh precede cells and nodes of the second mesh within the result mesh.
497 * The meshes must be of the same mesh type, else, an exception is thrown. The method
498 * MergeMeshes(), accepting a vector of input meshes, has no such a limitation.
499 * \param [in] mesh1 - the first mesh.
500 * \param [in] mesh2 - the second mesh.
501 * \return MEDCouplingMesh * - the result mesh. It is a new instance of
502 * MEDCouplingMesh. The caller is to delete this mesh using decrRef() as it
504 * \throw If the meshes are of different mesh type.
506 MEDCouplingMesh *MEDCouplingMesh::MergeMeshes(const MEDCouplingMesh *mesh1, const MEDCouplingMesh *mesh2)
509 throw INTERP_KERNEL::Exception("MEDCouplingMesh::MergeMeshes : first parameter is an empty mesh !");
511 throw INTERP_KERNEL::Exception("MEDCouplingMesh::MergeMeshes : second parameter is an empty mesh !");
512 return mesh1->mergeMyselfWith(mesh2);
516 * Creates a new MEDCouplingMesh by concatenating all given meshes, if possible.
518 * the *i*-th mesh precede cells and nodes of the (*i*+1)-th mesh within the result mesh.
519 * This method performs a systematic conversion to unstructured meshes before
520 * performing aggregation contrary to the other MergeMeshes()
521 * with two parameters that works only on the same type of meshes. So here it is possible
522 * to mix different type of meshes.
523 * \param [in] meshes - a vector of meshes to concatenate.
524 * \return MEDCouplingMesh * - the result mesh. It is a new instance of
525 * MEDCouplingUMesh. The caller is to delete this mesh using decrRef() as it
527 * \throw If \a meshes.size() == 0.
528 * \throw If \a size[ *i* ] == NULL.
529 * \throw If the coordinates is not set in none of the meshes.
530 * \throw If \a meshes[ *i* ]->getMeshDimension() < 0.
531 * \throw If the \a meshes are of different dimension (getMeshDimension()).
533 MEDCouplingMesh *MEDCouplingMesh::MergeMeshes(std::vector<const MEDCouplingMesh *>& meshes)
535 std::vector< MCAuto<MEDCouplingUMesh> > ms1(meshes.size());
536 std::vector< const MEDCouplingUMesh * > ms2(meshes.size());
537 for(std::size_t i=0;i<meshes.size();i++)
541 MEDCouplingUMesh *cur=meshes[i]->buildUnstructured();
542 ms1[i]=cur; ms2[i]=cur;
546 std::ostringstream oss; oss << "MEDCouplingMesh::MergeMeshes(std::vector<const MEDCouplingMesh *>& meshes) : mesh at pos #" << i << " of input vector of size " << meshes.size() << " is empty !";
547 throw INTERP_KERNEL::Exception(oss.str().c_str());
550 return MEDCouplingUMesh::MergeUMeshes(ms2);
554 * For example if \a type is INTERP_KERNEL::NORM_TRI3 , INTERP_KERNEL::NORM_POLYGON is returned.
555 * If \a type is INTERP_KERNEL::NORM_HEXA8 , INTERP_KERNEL::NORM_POLYHED is returned.
557 * \param [in] type the geometric type for which the corresponding dynamic type, is asked.
558 * \return the corresponding dynamic type, able to store the input \a type.
560 * \throw if type is equal to \c INTERP_KERNEL::NORM_ERROR or to an unexisting geometric type.
562 INTERP_KERNEL::NormalizedCellType MEDCouplingMesh::GetCorrespondingPolyType(INTERP_KERNEL::NormalizedCellType type)
564 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
565 return cm.getCorrespondingPolyType();
569 * \param [in] type the geometric type for which the number of nodes consituting it, is asked.
570 * \return number of nodes consituting the input geometric type \a type.
572 * \throw if type is dynamic as \c INTERP_KERNEL::NORM_POLYHED , \c INTERP_KERNEL::NORM_POLYGON , \c INTERP_KERNEL::NORM_QPOLYG
573 * \throw if type is equal to \c INTERP_KERNEL::NORM_ERROR or to an unexisting geometric type.
575 mcIdType MEDCouplingMesh::GetNumberOfNodesOfGeometricType(INTERP_KERNEL::NormalizedCellType type)
577 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
579 throw INTERP_KERNEL::Exception("MEDCouplingMesh::GetNumberOfNodesOfGeometricType : the input geometric type is dynamic ! Impossible to return a fixed number of nodes constituting it !");
580 return ToIdType( cm.getNumberOfNodes());
584 * \param [in] type the geometric type for which the status static/dynamic is asked.
585 * \return true for static geometric type, false for dynamic geometric type.
587 * \throw if type is equal to \c INTERP_KERNEL::NORM_ERROR or to an unexisting geometric type.
589 bool MEDCouplingMesh::IsStaticGeometricType(INTERP_KERNEL::NormalizedCellType type)
591 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
592 return !cm.isDynamic();
595 bool MEDCouplingMesh::IsLinearGeometricType(INTERP_KERNEL::NormalizedCellType type)
597 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
598 return !cm.isQuadratic();
602 * \param [in] type the geometric type for which the dimension is asked.
603 * \return the dimension associated to the input geometric type \a type.
605 * \throw if type is equal to \c INTERP_KERNEL::NORM_ERROR or to an unexisting geometric type.
607 int MEDCouplingMesh::GetDimensionOfGeometricType(INTERP_KERNEL::NormalizedCellType type)
609 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
610 return (int) cm.getDimension();
614 * \param [in] type the geometric type for which the representation is asked.
615 * \return the string representation corresponding to the input geometric type \a type.
617 * \throw if type is equal to \c INTERP_KERNEL::NORM_ERROR or to an unexisting geometric type.
619 const char *MEDCouplingMesh::GetReprOfGeometricType(INTERP_KERNEL::NormalizedCellType type)
621 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
626 * Finds cells in contact with several balls (i.e. points with precision).
627 * This method is an extension of getCellContainingPoint() and
628 * getCellsContainingPoint() for the case of multiple points.
629 * \param [in] pos - an array of coordinates of points in full interlace mode :
630 * X0,Y0,Z0,X1,Y1,Z1,... Size of the array must be \a
631 * this->getSpaceDimension() * \a nbOfPoints
632 * \param [in] nbOfPoints - number of points to locate within \a this mesh.
633 * \param [in] eps - radius of balls (i.e. the precision).
634 * \param [out] elts - vector returning ids of found cells.
635 * \param [out] eltsIndex - an array, of length \a nbOfPoints + 1,
636 * dividing cell ids in \a elts into groups each referring to one
637 * point. Its every element (except the last one) is an index pointing to the
638 * first id of a group of cells. For example cells in contact with the *i*-th
639 * point are described by following range of indices:
640 * [ \a eltsIndex[ *i* ], \a eltsIndex[ *i*+1 ] ) and the cell ids are
641 * \a elts[ \a eltsIndex[ *i* ]], \a elts[ \a eltsIndex[ *i* ] + 1 ], ...
642 * Number of cells in contact with the *i*-th point is
643 * \a eltsIndex[ *i*+1 ] - \a eltsIndex[ *i* ].
645 * \if ENABLE_EXAMPLES
646 * \ref cpp_mcumesh_getCellsContainingPoints "Here is a C++ example".<br>
647 * \ref py_mcumesh_getCellsContainingPoints "Here is a Python example".
650 void MEDCouplingMesh::getCellsContainingPoints(const double *pos, mcIdType nbOfPoints, double eps, MCAuto<DataArrayIdType>& elts, MCAuto<DataArrayIdType>& eltsIndex) const
652 eltsIndex=DataArrayIdType::New(); elts=DataArrayIdType::New(); eltsIndex->alloc(nbOfPoints+1,1); eltsIndex->setIJ(0,0,0); elts->alloc(0,1);
653 mcIdType *eltsIndexPtr(eltsIndex->getPointer());
654 int spaceDim(getSpaceDimension());
655 const double *work(pos);
656 for(mcIdType i=0;i<nbOfPoints;i++,work+=spaceDim)
658 std::vector<mcIdType> ret;
659 getCellsContainingPoint(work,eps,ret);
660 elts->insertAtTheEnd(ret.begin(),ret.end());
661 eltsIndexPtr[i+1]=elts->getNumberOfTuples();
666 * Behaves like MEDCouplingMesh::getCellsContainingPoints for cells in \a this that are linear.
667 * For quadratic cells in \a this, this method behaves by just considering linear part of cells.
668 * This method is here only for backward compatibility (interpolation GaussPoints to GaussPoints).
670 * \sa MEDCouplingMesh::getCellsContainingPoints, MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss
672 void MEDCouplingMesh::getCellsContainingPointsLinearPartOnlyOnNonDynType(const double *pos, mcIdType nbOfPoints, double eps, MCAuto<DataArrayIdType>& elts, MCAuto<DataArrayIdType>& eltsIndex) const
674 this->getCellsContainingPoints(pos,nbOfPoints,eps,elts,eltsIndex);
678 * Method computing center of mass of whole mesh (\a this)
679 * \return DataArrayDouble * - a new instance of DataArrayDouble with one tuple of n components where n is space dimension
681 MCAuto<DataArrayDouble> MEDCouplingMesh::computeMeshCenterOfMass() const
683 MCAuto<DataArrayDouble> cellCenters( this->computeCellCenterOfMass() );
684 MCAuto<MEDCouplingFieldDouble> vol( this->getMeasureField(true) );
685 MCAuto<DataArrayDouble> volXCenter( DataArrayDouble::Multiply(cellCenters,vol->getArray()) );
686 MCAuto<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(1, this->getSpaceDimension());
687 volXCenter->accumulate( ret->getPointer() );
688 double volOfMesh(vol->accumulate(0));
689 ret->applyLin(1.0/volOfMesh,0.0);
694 * Writes \a this mesh into a VTK format file named as specified.
695 * \param [in] fileName - the name of the file to write in. If the extension is OK the fileName will be used directly.
696 * If extension is invalid or no extension the right extension will be appended.
697 * \return - the real fileName
698 * \throw If \a fileName is not a writable file.
699 * \sa getVTKFileNameOf
701 std::string MEDCouplingMesh::writeVTK(const std::string& fileName, bool isBinary) const
703 std::string ret(getVTKFileNameOf(fileName));
706 MCAuto<DataArrayByte> byteArr;
708 { byteArr=DataArrayByte::New(); byteArr->alloc(0,1); }
709 writeVTKAdvanced(ret,cda,pda,byteArr);
714 * This method takes in input a file name \a fileName and considering the VTK extension of \a this (depending on the type of \a this)
715 * returns a right file name. If the input \a fileName has a valid extension the returned string is equal to \a fileName.
717 * \sa getVTKFileExtension
719 std::string MEDCouplingMesh::getVTKFileNameOf(const std::string& fileName) const
722 std::string part0,part1;
723 SplitExtension(fileName,part0,part1);
724 std::string ext("."); ext+=getVTKFileExtension();
733 void MEDCouplingMesh::writeVTKAdvanced(const std::string& fileName, const std::string& cda, const std::string& pda, DataArrayByte *byteData) const
735 std::ofstream ofs(fileName.c_str());
736 ofs << "<VTKFile type=\"" << getVTKDataSetType() << "\" version=\"0.1\" byte_order=\"" << MEDCouplingByteOrderStr() << "\">\n";
737 writeVTKLL(ofs,cda,pda,byteData);
740 ofs << "<AppendedData encoding=\"raw\">\n_1234";
741 ofs << std::flush; ofs.close();
742 std::ofstream ofs2(fileName.c_str(),std::ios_base::binary | std::ios_base::app);
743 ofs2.write(byteData->begin(),byteData->getNbOfElems()); ofs2 << std::flush; ofs2.close();
744 std::ofstream ofs3(fileName.c_str(),std::ios_base::app); ofs3 << "\n</AppendedData>\n</VTKFile>\n"; ofs3.close();
748 ofs << "</VTKFile>\n";
753 void MEDCouplingMesh::SplitExtension(const std::string& fileName, std::string& baseName, std::string& extension)
755 std::size_t pos(fileName.find_last_of('.'));
756 if(pos==std::string::npos)
762 baseName=fileName.substr(0,pos);
763 extension=fileName.substr(pos);