1 // Copyright (C) 2007-2013 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.
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.hxx"
24 #include "MEDCouplingFieldDouble.hxx"
25 #include "MEDCouplingFieldDiscretization.hxx"
26 #include "MEDCouplingAutoRefCountObjectPtr.hxx"
34 using namespace ParaMEDMEM;
36 MEDCouplingMesh::MEDCouplingMesh():_time(0.),_iteration(-1),_order(-1)
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)
46 std::size_t MEDCouplingMesh::getHeapMemorySize() 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 throw(INTERP_KERNEL::Exception)
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 << "\" !";
103 bool MEDCouplingMesh::isEqual(const MEDCouplingMesh *other, double prec) const throw(INTERP_KERNEL::Exception)
106 return isEqualIfNotWhy(other,prec,tmp);
110 * This method checks geo equivalence between two meshes : 'this' and 'other'.
111 * If no exception is throw 'this' and 'other' are geometrically equivalent regarding 'levOfCheck' level.
112 * This method is typically used to change the mesh of a field "safely" depending the 'levOfCheck' level considered.
113 * 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].
115 * @param levOfCheck input that specifies the level of check specified. The possible values are listed below.
116 * @param prec input that specifies precision for double float data used for comparison in meshes.
117 * @param cellCor output array not always informed (depending 'levOfCheck' param) that gives the corresponding array for cells from 'other' to 'this'.
118 * @param nodeCor output array not always informed (depending 'levOfCheck' param) that gives the corresponding array for nodes from 'other' to 'this'.
120 * Possible values for levOfCheck :
121 * - 0 for strict equality. This is the strongest level. 'cellCor' and 'nodeCor' params are never informed.
122 * - 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)
123 * - 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)
124 * - 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)
125 * - 2 for deep 'equality' as 0 option except that no control is done on all strings in mesh.
127 void MEDCouplingMesh::checkGeoEquivalWith(const MEDCouplingMesh *other, int levOfCheck, double prec,
128 DataArrayInt *&cellCor, DataArrayInt *&nodeCor) const throw(INTERP_KERNEL::Exception)
138 if(!isEqual(other,prec))
139 throw INTERP_KERNEL::Exception("checkGeoFitWith : Meshes are not equal !");
146 checkDeepEquivalWith(other,levOfCheck-10,prec,cellCor,nodeCor);
153 checkDeepEquivalOnSameNodesWith(other,levOfCheck-20,prec,cellCor);
158 checkFastEquivalWith(other,prec);
163 if(!isEqualWithoutConsideringStr(other,prec))
164 throw INTERP_KERNEL::Exception("checkGeoFitWith : Meshes are not equal without considering strings !");
168 throw INTERP_KERNEL::Exception("checkGeoFitWith : Invalid levOfCheck specified ! Value must be in 0,1,2,10,11 or 12.");
173 * Given a nodeIds range ['partBg','partEnd'), this method returns the set of cell ids in ascendant order whose connectivity of
174 * these cells are fully included in the range. As a consequence the returned set of cell ids does \b not \b always fit the nodes in ['partBg','partEnd')
175 * This method returns the corresponding cells in a newly created array that the caller has the responsability.
177 DataArrayInt *MEDCouplingMesh::getCellIdsFullyIncludedInNodeIds(const int *partBg, const int *partEnd) const
179 std::vector<int> crest;
180 std::set<int> p(partBg,partEnd);
181 int nbOfCells=getNumberOfCells();
182 for(int i=0;i<nbOfCells;i++)
184 std::vector<int> conn;
185 getNodeIdsOfCell(i,conn);
187 for(std::vector<int>::const_iterator iter=conn.begin();iter!=conn.end() && cont;iter++)
188 if(p.find(*iter)==p.end())
193 DataArrayInt *ret=DataArrayInt::New();
194 ret->alloc((int)crest.size(),1);
195 std::copy(crest.begin(),crest.end(),ret->getPointer());
200 * This method checks fastly that 'this' and 'other' are equal. All common checks are done here.
202 void MEDCouplingMesh::checkFastEquivalWith(const MEDCouplingMesh *other, double prec) const throw(INTERP_KERNEL::Exception)
204 if(getMeshDimension()!=other->getMeshDimension())
205 throw INTERP_KERNEL::Exception("checkFastEquivalWith : Mesh dimensions are not equal !");
206 if(getSpaceDimension()!=other->getSpaceDimension())
207 throw INTERP_KERNEL::Exception("checkFastEquivalWith : Space dimensions are not equal !");
208 if(getNumberOfCells()!=other->getNumberOfCells())
209 throw INTERP_KERNEL::Exception("checkFastEquivalWith : number of cells are not equal !");
213 * This method is very poor and looks only if 'this' and 'other' are candidate for merge of fields lying repectively on them.
215 bool MEDCouplingMesh::areCompatibleForMerge(const MEDCouplingMesh *other) const
217 if(getMeshDimension()!=other->getMeshDimension())
219 if(getSpaceDimension()!=other->getSpaceDimension())
225 * This method builds a field lying on 'this' with 'nbOfComp' components.
226 * 'func' is a pointer that points to a function that takes 2 arrays in parameter and returns a boolean.
227 * The first array is a in-param of size this->getSpaceDimension and the second an out param of size 'nbOfComp'.
228 * The return field will have type specified by 't'. 't' is also used to determine where values of field will be
230 * Contrary to other fillFromAnalytic methods this method requests a C++ function pointer as input.
231 * The 'func' is a callback that takes as first parameter an input array of size 'this->getSpaceDimension()',
232 * the second parameter is a pointer on a valid zone of size at least equal to 'nbOfComp' values. And too finish
233 * the returned value is a boolean that is equal to False in case of invalid evaluation (log(0) for example...)
234 * @param t type of field returned and specifies where the evaluation of func will be done.
235 * @param nbOfComp number of components of returned field.
236 * @param func pointer to a function that should return false if the evaluation failed. (division by 0. for example)
237 * @return field with counter = 1.
239 MEDCouplingFieldDouble *MEDCouplingMesh::fillFromAnalytic(TypeOfField t, int nbOfComp, FunctionToEvaluate func) const
241 MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(t,ONE_TIME);
243 ret->fillFromAnalytic(nbOfComp,func);
244 ret->synchronizeTimeWithSupport();
249 * This method copyies all tiny strings from other (name and components name).
250 * @throw if other and this have not same mesh type.
252 void MEDCouplingMesh::copyTinyStringsFrom(const MEDCouplingMesh *other) throw(INTERP_KERNEL::Exception)
255 _description=other->_description;
256 _time_unit=other->_time_unit;
260 * This method copies all attributes that are \b NOT arrays in this.
261 * All tiny attributes not usefully for state of 'this' are ignored.
263 void MEDCouplingMesh::copyTinyInfoFrom(const MEDCouplingMesh *other) throw(INTERP_KERNEL::Exception)
265 copyTinyStringsFrom(other);
267 _iteration=other->_iteration;
268 _order=other->_order;
272 * This method builds a field lying on 'this' with 'nbOfComp' components.
273 * 'func' is a string that is the expression to evaluate.
274 * The return field will have type specified by 't'. 't' is also used to determine where values of field will be
276 * This method is equivalent to those taking a C++ function pointer except that here the 'func' is informed by
277 * an interpretable input string.
279 * The dynamic interpretor uses \b alphabetical \b order to assign the component id to the var name.
281 * - "2*x+z" func : x stands for component #0 and z stands for component #1 \b NOT #2 !
283 * Some var names are reserved and have special meaning. IVec stands for (1,0,0,...). JVec stands for (0,1,0...).
284 * KVec stands for (0,0,1,...)... These keywords allows too differentate the evaluation of output components each other.
286 * If 'nbOfComp' equals to 4 for example and that 'this->getSpaceDimension()' equals to 3.
288 * For the input tuple T = (1.,3.,7.) :
289 * - '2*x+z' will return (5.,5.,5.,5.)
290 * - '2*x+0*y+z' will return (9.,9.,9.,9.)
291 * - '2*x*IVec+(x+z)*LVec' will return (2.,0.,0.,4.)
292 * - '2*x*IVec+(y+z)*KVec' will return (2.,0.,10.,0.)
294 * @param t type of field returned and specifies where the evaluation of func will be done.
295 * @param nbOfComp number of components of returned field.
296 * @param func expression.
297 * @return field with counter = 1.
299 MEDCouplingFieldDouble *MEDCouplingMesh::fillFromAnalytic(TypeOfField t, int nbOfComp, const char *func) const
301 MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(t,ONE_TIME);
303 ret->fillFromAnalytic(nbOfComp,func);
304 ret->synchronizeTimeWithSupport();
309 * This method builds a field lying on 'this' with 'nbOfComp' components.
310 * 'func' is a string that is the expression to evaluate.
311 * The return field will have type specified by 't'. 't' is also used to determine where values of field will be
312 * evaluate. This method is different than MEDCouplingMesh::fillFromAnalytic, because the info on components are used here to determine vars pos in 'func'.
314 * @param t type of field returned and specifies where the evaluation of func will be done.
315 * @param nbOfComp number of components of returned field.
316 * @param func expression.
317 * @return field with counter = 1.
319 MEDCouplingFieldDouble *MEDCouplingMesh::fillFromAnalytic2(TypeOfField t, int nbOfComp, const char *func) const
321 MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(t,ONE_TIME);
323 ret->fillFromAnalytic2(nbOfComp,func);
324 ret->synchronizeTimeWithSupport();
329 * This method builds a field lying on 'this' with 'nbOfComp' components.
330 * 'func' is a string that is the expression to evaluate.
331 * The return field will have type specified by 't'. 't' is also used to determine where values of field will be
332 * evaluate. This method is different than MEDCouplingMesh::fillFromAnalytic, because 'varsOrder' specifies the pos to assign of vars in 'func'.
334 * @param t type of field returned and specifies where the evaluation of func will be done.
335 * @param nbOfComp number of components of returned field.
336 * @param func expression.
337 * @return field with counter = 1.
339 MEDCouplingFieldDouble *MEDCouplingMesh::fillFromAnalytic3(TypeOfField t, int nbOfComp, const std::vector<std::string>& varsOrder, const char *func) const
341 MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(t,ONE_TIME);
343 ret->fillFromAnalytic3(nbOfComp,varsOrder,func);
344 ret->synchronizeTimeWithSupport();
349 * retruns a newly created mesh with counter=1
350 * that is the union of \b mesh1 and \b mesh2 if possible. The cells of \b mesh2 will appear after cells of \b mesh1. Idem for nodes.
351 * The only contraint is that \b mesh1 an \b mesh2 have the same mesh types. If it is not the case please use the other API of MEDCouplingMesh::MergeMeshes,
352 * with input vector of meshes.
354 MEDCouplingMesh *MEDCouplingMesh::MergeMeshes(const MEDCouplingMesh *mesh1, const MEDCouplingMesh *mesh2) throw(INTERP_KERNEL::Exception)
357 throw INTERP_KERNEL::Exception("MEDCouplingMesh::MergeMeshes : first parameter is an empty mesh !");
359 throw INTERP_KERNEL::Exception("MEDCouplingMesh::MergeMeshes : second parameter is an empty mesh !");
360 return mesh1->mergeMyselfWith(mesh2);
364 * retruns a newly created mesh with counter=1
365 * that is the union of meshes if possible. The cells of \b meshes[1] will appear after cells of \b meshes[0]. Idem for nodes.
366 * This method performs a systematic conversion to unstructured meshes before performing aggregation contrary to the other ParaMEDMEM::MEDCouplingMesh::MergeMeshes with
367 * two parameters that work only on the same type of meshes. So here it is possible to mix different type of meshes.
369 MEDCouplingMesh *MEDCouplingMesh::MergeMeshes(std::vector<const MEDCouplingMesh *>& meshes) throw(INTERP_KERNEL::Exception)
371 std::vector< MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> > ms1(meshes.size());
372 std::vector< const MEDCouplingUMesh * > ms2(meshes.size());
373 for(std::size_t i=0;i<meshes.size();i++)
377 MEDCouplingUMesh *cur=meshes[i]->buildUnstructured();
378 ms1[i]=cur; ms2[i]=cur;
382 std::ostringstream oss; oss << "MEDCouplingMesh::MergeMeshes(std::vector<const MEDCouplingMesh *>& meshes) : mesh at pos #" << i << " of input vector of size " << meshes.size() << " is empty !";
383 throw INTERP_KERNEL::Exception(oss.str().c_str());
386 return MEDCouplingUMesh::MergeUMeshes(ms2);
390 * \param [in] type the geometric type for which the dimension is asked.
391 * \return the dimension associated to the input geometric type \a type.
393 * \throw if type is equal to \c INTERP_KERNEL::NORM_ERROR or to an unexisting geometric type.
395 int MEDCouplingMesh::GetDimensionOfGeometricType(INTERP_KERNEL::NormalizedCellType type) throw(INTERP_KERNEL::Exception)
397 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
398 return (int) cm.getDimension();
402 * \param [in] type the geometric type for which the representation is asked.
403 * \return the string representation corresponding to the input geometric type \a type.
405 * \throw if type is equal to \c INTERP_KERNEL::NORM_ERROR or to an unexisting geometric type.
407 const char *MEDCouplingMesh::GetReprOfGeometricType(INTERP_KERNEL::NormalizedCellType type) throw(INTERP_KERNEL::Exception)
409 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
413 void MEDCouplingMesh::getCellsContainingPoint(const double *pos, double eps, std::vector<int>& elts) const
415 int ret=getCellContainingPoint(pos,eps);
419 void MEDCouplingMesh::getCellsContainingPoints(const double *pos, int nbOfPoints, double eps, std::vector<int>& elts, std::vector<int>& eltsIndex) const
421 eltsIndex.resize(nbOfPoints+1);
424 int spaceDim=getSpaceDimension();
425 const double *work=pos;
426 for(int i=0;i<nbOfPoints;i++,work+=spaceDim)
428 int ret=getCellContainingPoint(work,eps);
432 eltsIndex[i+1]=eltsIndex[i]+1;
435 eltsIndex[i+1]=eltsIndex[i];
440 * This method writes a file in VTK format into file 'fileName'.
441 * An exception is thrown if the file is not writable.
443 void MEDCouplingMesh::writeVTK(const char *fileName) const throw(INTERP_KERNEL::Exception)
446 writeVTKAdvanced(fileName,cda,pda);
449 void MEDCouplingMesh::writeVTKAdvanced(const char *fileName, const std::string& cda, const std::string& pda) const throw(INTERP_KERNEL::Exception)
451 std::ofstream ofs(fileName);
452 ofs << "<VTKFile type=\"" << getVTKDataSetType() << "\" version=\"0.1\" byte_order=\"LittleEndian\">\n";
453 writeVTKLL(ofs,cda,pda);
454 ofs << "</VTKFile>\n";