1 // Copyright (C) 2007-2023 CEA, EDF
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 (EDF R&D)
21 #include "MEDCouplingMemArray.txx"
24 #include "GenMathFormulae.hxx"
25 #include "InterpKernelAutoPtr.hxx"
26 #include "InterpKernelExprParser.hxx"
28 #include "InterpKernelAutoPtr.hxx"
29 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
30 #include "InterpKernelAutoPtr.hxx"
31 #include "InterpKernelGeo2DNode.hxx"
32 #include "InterpKernelGeo2DEdgeLin.hxx"
41 typedef double (*MYFUNCPTR)(double);
43 using namespace MEDCoupling;
45 template class MEDCOUPLING_EXPORT MEDCoupling::MemArray<mcIdType>;
46 template class MEDCOUPLING_EXPORT MEDCoupling::MemArray<double>;
47 template class MEDCOUPLING_EXPORT MEDCoupling::DataArrayTemplate<mcIdType>;
48 template class MEDCOUPLING_EXPORT MEDCoupling::DataArrayTemplate<double>;
49 template class MEDCOUPLING_EXPORT MEDCoupling::DataArrayTemplateClassic<Int32>;
50 template class MEDCOUPLING_EXPORT MEDCoupling::DataArrayTemplateClassic<Int64>;
51 template class MEDCOUPLING_EXPORT MEDCoupling::DataArrayTemplateClassic<double>;
52 template class MEDCOUPLING_EXPORT MEDCoupling::DataArrayTemplateFP<double>;
53 template class MEDCOUPLING_EXPORT MEDCoupling::DataArrayIterator<double>;
54 template class MEDCOUPLING_EXPORT MEDCoupling::DataArrayIterator<mcIdType>;
55 template class MEDCOUPLING_EXPORT MEDCoupling::DataArrayDiscrete<Int32>;
56 template class MEDCOUPLING_EXPORT MEDCoupling::DataArrayDiscreteSigned<Int32>;
57 template class MEDCOUPLING_EXPORT MEDCoupling::DataArrayDiscrete<Int64>;
58 template class MEDCOUPLING_EXPORT MEDCoupling::DataArrayDiscreteSigned<Int64>;
59 template class MEDCOUPLING_EXPORT MEDCoupling::DataArrayTuple<mcIdType>;
60 template class MEDCOUPLING_EXPORT MEDCoupling::DataArrayTuple<double>;
61 template class MEDCOUPLING_EXPORT MEDCoupling::DataArrayTuple<float>;
63 void MEDCoupling::DACheckNbOfTuplesAndComp(const DataArray *da, mcIdType nbOfTuples, std::size_t nbOfCompo, const std::string& msg)
66 throw INTERP_KERNEL::Exception("DACheckNbOfTuplesAndComp : null input object !");
67 da->checkNbOfTuplesAndComp(nbOfTuples,nbOfCompo,msg);
70 template<int SPACEDIM>
71 void DataArrayDouble::findCommonTuplesAlg(const double *bbox, mcIdType nbNodes, mcIdType limitNodeId, double prec, DataArrayIdType *c, DataArrayIdType *cI) const
73 const double *coordsPtr=getConstPointer();
74 BBTreePts<SPACEDIM,mcIdType> myTree(bbox,0,0,nbNodes,prec);
75 std::vector<bool> isDone(nbNodes);
76 for(mcIdType i=0;i<nbNodes;i++)
80 std::vector<mcIdType> intersectingElems;
81 myTree.getElementsAroundPoint(coordsPtr+i*SPACEDIM,intersectingElems);
82 if(intersectingElems.size()>1)
84 std::vector<mcIdType> commonNodes;
85 for(std::vector<mcIdType>::const_iterator it=intersectingElems.begin();it!=intersectingElems.end();it++)
89 commonNodes.push_back(*it);
92 if(!commonNodes.empty())
94 cI->pushBackSilent(cI->back()+ToIdType(commonNodes.size())+1);
96 c->insertAtTheEnd(commonNodes.begin(),commonNodes.end());
103 template<int SPACEDIM>
104 void DataArrayDouble::FindTupleIdsNearTuplesAlg(const BBTreePts<SPACEDIM,mcIdType>& myTree, const double *pos, mcIdType nbOfTuples, double eps,
105 DataArrayIdType *c, DataArrayIdType *cI)
107 for(mcIdType i=0;i<nbOfTuples;i++)
109 std::vector<mcIdType> intersectingElems;
110 myTree.getElementsAroundPoint(pos+i*SPACEDIM,intersectingElems);
111 std::vector<mcIdType> commonNodes;
112 for(std::vector<mcIdType>::const_iterator it=intersectingElems.begin();it!=intersectingElems.end();it++)
113 commonNodes.push_back(*it);
114 cI->pushBackSilent(cI->back()+ToIdType(commonNodes.size()));
115 c->insertAtTheEnd(commonNodes.begin(),commonNodes.end());
119 template<int SPACEDIM>
120 void DataArrayDouble::FindClosestTupleIdAlg(const BBTreePts<SPACEDIM,mcIdType>& myTree, double dist, const double *pos, mcIdType nbOfTuples, const double *thisPt, mcIdType thisNbOfTuples, mcIdType *res)
122 double distOpt = std::max(dist, std::numeric_limits<double>::epsilon());
123 const double *p(pos);
125 for(mcIdType i=0;i<nbOfTuples;i++,p+=SPACEDIM,r++)
130 double ret=myTree.getElementsAroundPoint2(p,distOpt,elem);
131 if(ret!=std::numeric_limits<double>::max())
133 distOpt=std::max(ret,1e-4);
138 { distOpt=2*distOpt; continue; }
143 mcIdType DataArray::EffectiveCircPerm(mcIdType nbOfShift, mcIdType nbOfTuples)
146 throw INTERP_KERNEL::Exception("DataArray::EffectiveCircPerm : number of tuples is expected to be > 0 !");
149 return nbOfShift%nbOfTuples;
153 mcIdType tmp(-nbOfShift);
155 return nbOfTuples-tmp;
159 std::size_t DataArray::getHeapMemorySizeWithoutChildren() const
161 std::size_t sz1=_name.capacity();
162 std::size_t sz2=_info_on_compo.capacity();
164 for(std::vector<std::string>::const_iterator it=_info_on_compo.begin();it!=_info_on_compo.end();it++)
165 sz3+=(*it).capacity();
169 std::vector<const BigMemoryObject *> DataArray::getDirectChildrenWithNull() const
171 return std::vector<const BigMemoryObject *>();
175 * Sets the attribute \a _name of \a this array.
176 * See \ref MEDCouplingArrayBasicsName "DataArrays infos" for more information.
177 * \param [in] name - new array name
179 void DataArray::setName(const std::string& name)
185 * Copies textual data from an \a other DataArray. The copied data are
186 * - the name attribute,
187 * - the information of components.
189 * For more information on these data see \ref MEDCouplingArrayBasicsName "DataArrays infos".
191 * \param [in] other - another instance of DataArray to copy the textual data from.
192 * \throw If number of components of \a this array differs from that of the \a other.
194 void DataArray::copyStringInfoFrom(const DataArray& other)
196 if(_info_on_compo.size()!=other._info_on_compo.size())
197 throw INTERP_KERNEL::Exception("Size of arrays mismatches on copyStringInfoFrom !");
199 _info_on_compo=other._info_on_compo;
202 void DataArray::copyPartOfStringInfoFrom(const DataArray& other, const std::vector<std::size_t>& compoIds)
204 std::size_t nbOfCompoOth=other.getNumberOfComponents();
205 std::size_t newNbOfCompo=compoIds.size();
206 for(std::size_t i=0;i<newNbOfCompo;i++)
207 if(compoIds[i]>=nbOfCompoOth)
209 std::ostringstream oss; oss << "Specified component id is out of range (" << compoIds[i] << ") compared with nb of actual components (" << nbOfCompoOth << ")";
210 throw INTERP_KERNEL::Exception(oss.str().c_str());
212 for(std::size_t i=0;i<newNbOfCompo;i++)
213 setInfoOnComponent(i,other.getInfoOnComponent(compoIds[i]));
216 void DataArray::copyPartOfStringInfoFrom2(const std::vector<std::size_t>& compoIds, const DataArray& other)
218 if(compoIds.size()!=other.getNumberOfComponents())
219 throw INTERP_KERNEL::Exception("Given compoIds has not the same size as number of components of given array !");
220 std::size_t partOfCompoToSet=compoIds.size();
221 std::size_t nbOfCompo=getNumberOfComponents();
222 for(std::size_t i=0;i<partOfCompoToSet;i++)
223 if(compoIds[i]>=nbOfCompo)
225 std::ostringstream oss; oss << "Specified component id is out of range (" << compoIds[i] << ") compared with nb of actual components (" << nbOfCompo << ")";
226 throw INTERP_KERNEL::Exception(oss.str().c_str());
228 for(std::size_t i=0;i<partOfCompoToSet;i++)
229 setInfoOnComponent(compoIds[i],other.getInfoOnComponent(i));
232 bool DataArray::areInfoEqualsIfNotWhy(const DataArray& other, std::string& reason) const
234 std::ostringstream oss;
235 if(_name!=other._name)
237 oss << "Names DataArray mismatch : this name=\"" << _name << " other name=\"" << other._name << "\" !";
241 if(_info_on_compo!=other._info_on_compo)
243 oss << "Components DataArray mismatch : \nThis components=";
244 for(std::vector<std::string>::const_iterator it=_info_on_compo.begin();it!=_info_on_compo.end();it++)
245 oss << "\"" << *it << "\",";
246 oss << "\nOther components=";
247 for(std::vector<std::string>::const_iterator it=other._info_on_compo.begin();it!=other._info_on_compo.end();it++)
248 oss << "\"" << *it << "\",";
256 * Compares textual information of \a this DataArray with that of an \a other one.
257 * The compared data are
258 * - the name attribute,
259 * - the information of components.
261 * For more information on these data see \ref MEDCouplingArrayBasicsName "DataArrays infos".
262 * \param [in] other - another instance of DataArray to compare the textual data of.
263 * \return bool - \a true if the textual information is same, \a false else.
265 bool DataArray::areInfoEquals(const DataArray& other) const
268 return areInfoEqualsIfNotWhy(other,tmp);
271 void DataArray::reprWithoutNameStream(std::ostream& stream) const
273 stream << "Number of components : "<< getNumberOfComponents() << "\n";
274 stream << "Info of these components : ";
275 for(std::vector<std::string>::const_iterator iter=_info_on_compo.begin();iter!=_info_on_compo.end();iter++)
276 stream << "\"" << *iter << "\" ";
280 std::string DataArray::cppRepr(const std::string& varName) const
282 std::ostringstream ret;
283 reprCppStream(varName,ret);
288 * Sets information on all components. To know more on format of this information
289 * see \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
290 * \param [in] info - a vector of strings.
291 * \throw If size of \a info differs from the number of components of \a this.
293 void DataArray::setInfoOnComponents(const std::vector<std::string>& info)
295 if(getNumberOfComponents()!=info.size())
297 std::ostringstream oss; oss << "DataArray::setInfoOnComponents : input is of size " << info.size() << " whereas number of components is equal to " << getNumberOfComponents() << " !";
298 throw INTERP_KERNEL::Exception(oss.str().c_str());
304 * This method is only a dispatcher towards DataArrayDouble::setPartOfValues3, DataArrayInt::setPartOfValues3, DataArrayChar::setPartOfValues3 depending on the true
305 * type of \a this and \a aBase.
307 * \throw If \a aBase and \a this do not have the same type.
309 * \sa DataArrayDouble::setPartOfValues3, DataArrayInt::setPartOfValues3, DataArrayChar::setPartOfValues3.
311 void DataArray::setPartOfValuesBase3(const DataArray *aBase, const mcIdType *bgTuples, const mcIdType *endTuples, mcIdType bgComp, mcIdType endComp, mcIdType stepComp, bool strictCompoCompare)
314 throw INTERP_KERNEL::Exception("DataArray::setPartOfValuesBase3 : input aBase object is NULL !");
315 DataArrayDouble *this1(dynamic_cast<DataArrayDouble *>(this));
316 DataArrayIdType *this2(dynamic_cast<DataArrayIdType *>(this));
317 DataArrayChar *this3(dynamic_cast<DataArrayChar *>(this));
318 const DataArrayDouble *a1(dynamic_cast<const DataArrayDouble *>(aBase));
319 const DataArrayIdType *a2(dynamic_cast<const DataArrayIdType *>(aBase));
320 const DataArrayChar *a3(dynamic_cast<const DataArrayChar *>(aBase));
323 this1->setPartOfValues3(a1,bgTuples,endTuples,bgComp,endComp,stepComp,strictCompoCompare);
328 this2->setPartOfValues3(a2,bgTuples,endTuples,bgComp,endComp,stepComp,strictCompoCompare);
333 this3->setPartOfValues3(a3,bgTuples,endTuples,bgComp,endComp,stepComp,strictCompoCompare);
336 throw INTERP_KERNEL::Exception("DataArray::setPartOfValuesBase3 : input aBase object and this do not have the same type !");
339 std::vector<std::string> DataArray::getVarsOnComponent() const
341 std::size_t nbOfCompo=_info_on_compo.size();
342 std::vector<std::string> ret(nbOfCompo);
343 for(std::size_t i=0;i<nbOfCompo;i++)
344 ret[i]=getVarOnComponent(i);
348 std::vector<std::string> DataArray::getUnitsOnComponent() const
350 std::size_t nbOfCompo=_info_on_compo.size();
351 std::vector<std::string> ret(nbOfCompo);
352 for(std::size_t i=0;i<nbOfCompo;i++)
353 ret[i]=getUnitOnComponent(i);
358 * Returns information on a component specified by an index.
359 * To know more on format of this information
360 * see \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
361 * \param [in] i - the index (zero based) of the component of interest.
362 * \return std::string - a string containing the information on \a i-th component.
363 * \throw If \a i is not a valid component index.
365 std::string DataArray::getInfoOnComponent(std::size_t i) const
367 if(i<_info_on_compo.size())
368 return _info_on_compo[i];
371 std::ostringstream oss; oss << "DataArray::getInfoOnComponent : Specified component id is out of range (" << i << ") compared with nb of actual components (" << _info_on_compo.size();
372 throw INTERP_KERNEL::Exception(oss.str().c_str());
377 * Returns the var part of the full information of the \a i-th component.
378 * For example, if \c getInfoOnComponent(0) returns "SIGXY [N/m^2]", then
379 * \c getVarOnComponent(0) returns "SIGXY".
380 * If a unit part of information is not detected by presence of
381 * two square brackets, then the full information is returned.
382 * To read more about the component information format, see
383 * \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
384 * \param [in] i - the index (zero based) of the component of interest.
385 * \return std::string - a string containing the var information, or the full info.
386 * \throw If \a i is not a valid component index.
388 std::string DataArray::getVarOnComponent(std::size_t i) const
390 if(i<_info_on_compo.size())
392 return GetVarNameFromInfo(_info_on_compo[i]);
396 std::ostringstream oss; oss << "DataArray::getVarOnComponent : Specified component id is out of range (" << i << ") compared with nb of actual components (" << _info_on_compo.size();
397 throw INTERP_KERNEL::Exception(oss.str().c_str());
402 * Returns the unit part of the full information of the \a i-th component.
403 * For example, if \c getInfoOnComponent(0) returns "SIGXY [ N/m^2]", then
404 * \c getUnitOnComponent(0) returns " N/m^2".
405 * If a unit part of information is not detected by presence of
406 * two square brackets, then an empty string is returned.
407 * To read more about the component information format, see
408 * \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
409 * \param [in] i - the index (zero based) of the component of interest.
410 * \return std::string - a string containing the unit information, if any, or "".
411 * \throw If \a i is not a valid component index.
413 std::string DataArray::getUnitOnComponent(std::size_t i) const
415 if(i<_info_on_compo.size())
417 return GetUnitFromInfo(_info_on_compo[i]);
421 std::ostringstream oss; oss << "DataArray::getUnitOnComponent : Specified component id is out of range (" << i << ") compared with nb of actual components (" << _info_on_compo.size();
422 throw INTERP_KERNEL::Exception(oss.str().c_str());
427 * Split \a st string into chunks of sz characters each.
428 * Size of input \a st must a be a multiple of \a sz. If not an exception is thrown.
430 std::vector<std::string> DataArray::SplitStringInChuncks(const std::string st, std::size_t sz)
432 std::size_t len = st.length();
433 std::size_t nbOfCompo(len/sz);
434 if( nbOfCompo*sz != len)
436 THROW_IK_EXCEPTION("DataArray::SplitStringInChuncks : Length of input string (" << len << ") is not equal to " << nbOfCompo << "*" << sz << " !");
438 std::vector<std::string> ret(nbOfCompo);
439 for(std::size_t i = 0 ; i < nbOfCompo ; ++i)
441 std::string part = st.substr(i*sz,sz);
442 std::size_t p3=part.find_last_not_of(" \t");
443 part = part.substr(0,p3+1);
450 * Returns the var part of the full component information.
451 * For example, if \a info == "SIGXY [N/m^2]", then this method returns "SIGXY".
452 * If a unit part of information is not detected by presence of
453 * two square brackets, then the whole \a info is returned.
454 * To read more about the component information format, see
455 * \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
456 * \param [in] info - the full component information.
457 * \return std::string - a string containing only var information, or the \a info.
459 std::string DataArray::GetVarNameFromInfo(const std::string& info)
461 std::size_t p1=info.find_last_of('[');
462 std::size_t p2=info.find_last_of(']');
463 if(p1==std::string::npos || p2==std::string::npos)
468 return std::string();
469 std::size_t p3=info.find_last_not_of(' ',p1-1);
470 return info.substr(0,p3+1);
474 * Returns the unit part of the full component information.
475 * For example, if \a info == "SIGXY [ N/m^2]", then this method returns " N/m^2".
476 * If a unit part of information is not detected by presence of
477 * two square brackets, then an empty string is returned.
478 * To read more about the component information format, see
479 * \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
480 * \param [in] info - the full component information.
481 * \return std::string - a string containing only unit information, if any, or "".
483 std::string DataArray::GetUnitFromInfo(const std::string& info)
485 std::size_t p1=info.find_last_of('[');
486 std::size_t p2=info.find_last_of(']');
487 if(p1==std::string::npos || p2==std::string::npos)
488 return std::string();
490 return std::string();
491 return info.substr(p1+1,p2-p1-1);
495 * This method put in info format the result of the merge of \a var and \a unit.
496 * The standard format for that is "var [unit]".
497 * Inversely you can retrieve the var part or the unit part of info string using resp. GetVarNameFromInfo and GetUnitFromInfo.
499 std::string DataArray::BuildInfoFromVarAndUnit(const std::string& var, const std::string& unit)
501 std::ostringstream oss;
502 oss << var << " [" << unit << "]";
506 std::string DataArray::GetAxisTypeRepr(MEDCouplingAxisType at)
511 return std::string("AX_CART");
513 return std::string("AX_CYL");
515 return std::string("AX_SPHER");
517 throw INTERP_KERNEL::Exception("DataArray::GetAxisTypeRepr : unrecognized axis type enum !");
522 * Returns a new DataArray by concatenating all given arrays, so that (1) the number
523 * of tuples in the result array is a sum of the number of tuples of given arrays and (2)
524 * the number of component in the result array is same as that of each of given arrays.
525 * Info on components is copied from the first of the given arrays. Number of components
526 * in the given arrays must be the same.
527 * \param [in] arrs - a sequence of arrays to include in the result array. All arrays must have the same type.
528 * \return DataArray * - the new instance of DataArray (that can be either DataArrayInt, DataArrayDouble, DataArrayChar).
529 * The caller is to delete this result array using decrRef() as it is no more
531 * \throw If all arrays within \a arrs are NULL.
532 * \throw If all not null arrays in \a arrs have not the same type.
533 * \throw If getNumberOfComponents() of arrays within \a arrs.
535 DataArray *DataArray::Aggregate(const std::vector<const DataArray *>& arrs)
537 std::vector<const DataArray *> arr2;
538 for(std::vector<const DataArray *>::const_iterator it=arrs.begin();it!=arrs.end();it++)
542 throw INTERP_KERNEL::Exception("DataArray::Aggregate : only null instance in input vector !");
543 std::vector<const DataArrayDouble *> arrd;
544 std::vector<const DataArrayIdType *> arri;
545 std::vector<const DataArrayChar *> arrc;
546 for(std::vector<const DataArray *>::const_iterator it=arr2.begin();it!=arr2.end();it++)
548 const DataArrayDouble *a=dynamic_cast<const DataArrayDouble *>(*it);
550 { arrd.push_back(a); continue; }
551 const DataArrayIdType *b=dynamic_cast<const DataArrayIdType *>(*it);
553 { arri.push_back(b); continue; }
554 const DataArrayChar *c=dynamic_cast<const DataArrayChar *>(*it);
556 { arrc.push_back(c); continue; }
557 throw INTERP_KERNEL::Exception("DataArray::Aggregate : presence of not null instance in inuput that is not in [DataArrayDouble, DataArrayInt, DataArrayChar] !");
559 if(arr2.size()==arrd.size())
560 return DataArrayDouble::Aggregate(arrd);
561 if(arr2.size()==arri.size())
562 return DataArrayIdType::Aggregate(arri);
563 if(arr2.size()==arrc.size())
564 return DataArrayChar::Aggregate(arrc);
565 throw INTERP_KERNEL::Exception("DataArray::Aggregate : all input arrays must have the same type !");
569 * Sets information on a component specified by an index.
570 * To know more on format of this information
571 * see \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
572 * \warning Don't pass NULL as \a info!
573 * \param [in] i - the index (zero based) of the component of interest.
574 * \param [in] info - the string containing the information.
575 * \throw If \a i is not a valid component index.
577 void DataArray::setInfoOnComponent(std::size_t i, const std::string& info)
579 if(i<_info_on_compo.size())
580 _info_on_compo[i]=info;
583 std::ostringstream oss; oss << "DataArray::setInfoOnComponent : Specified component id is out of range (" << i << ") compared with nb of actual components (" << _info_on_compo.size();
584 throw INTERP_KERNEL::Exception(oss.str().c_str());
589 * Sets information on all components. This method can change number of components
590 * at certain conditions; if the conditions are not respected, an exception is thrown.
591 * The number of components can be changed in \a this only if \a this is not allocated.
592 * The condition of number of components must not be changed.
594 * To know more on format of the component information see
595 * \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
596 * \param [in] info - a vector of component infos.
597 * \throw If \a this->getNumberOfComponents() != \a info.size() && \a this->isAllocated()
599 void DataArray::setInfoAndChangeNbOfCompo(const std::vector<std::string>& info)
601 if(getNumberOfComponents()!=info.size())
607 std::ostringstream oss; oss << "DataArray::setInfoAndChangeNbOfCompo : input is of size " << info.size() << " whereas number of components is equal to " << getNumberOfComponents() << " and this is already allocated !";
608 throw INTERP_KERNEL::Exception(oss.str().c_str());
615 void DataArray::checkNbOfTuples(mcIdType nbOfTuples, const std::string& msg) const
617 if(getNumberOfTuples()!=nbOfTuples)
619 std::ostringstream oss; oss << msg << " : mismatch number of tuples : expected " << nbOfTuples << " having " << getNumberOfTuples() << " !";
620 throw INTERP_KERNEL::Exception(oss.str().c_str());
624 void DataArray::checkNbOfComps(std::size_t nbOfCompo, const std::string& msg) const
626 if (getNumberOfComponents()!=nbOfCompo)
628 std::ostringstream oss; oss << msg << " : mismatch number of components : expected " << nbOfCompo << " having " << getNumberOfComponents() << " !";
629 throw INTERP_KERNEL::Exception(oss.str().c_str());
633 void DataArray::checkNbOfElems(mcIdType nbOfElems, const std::string& msg) const
635 if(getNbOfElems()!=nbOfElems)
637 std::ostringstream oss; oss << msg << " : mismatch number of elems : Expected " << nbOfElems << " having " << getNbOfElems() << " !";
638 throw INTERP_KERNEL::Exception(oss.str().c_str());
642 void DataArray::checkNbOfTuplesAndComp(const DataArray& other, const std::string& msg) const
644 if(getNumberOfTuples()!=other.getNumberOfTuples())
646 std::ostringstream oss; oss << msg << " : mismatch number of tuples : expected " << other.getNumberOfTuples() << " having " << getNumberOfTuples() << " !";
647 throw INTERP_KERNEL::Exception(oss.str().c_str());
649 if(getNumberOfComponents()!=other.getNumberOfComponents())
651 std::ostringstream oss; oss << msg << " : mismatch number of components : expected " << other.getNumberOfComponents() << " having " << getNumberOfComponents() << " !";
652 throw INTERP_KERNEL::Exception(oss.str().c_str());
656 void DataArray::checkNbOfTuplesAndComp(mcIdType nbOfTuples, std::size_t nbOfCompo, const std::string& msg) const
658 checkNbOfTuples(nbOfTuples,msg);
659 checkNbOfComps(nbOfCompo,msg);
663 * Simply this method checks that \b value is in [0,\b ref).
665 void DataArray::CheckValueInRange(mcIdType ref, mcIdType value, const std::string& msg)
667 if(value<0 || value>=ref)
669 std::ostringstream oss; oss << "DataArray::CheckValueInRange : " << msg << " ! Expected in range [0," << ref << "[ having " << value << " !";
670 throw INTERP_KERNEL::Exception(oss.str().c_str());
675 * This method checks that [\b start, \b end) is compliant with ref length \b value.
676 * typically start in [0,\b value) and end in [0,\b value). If value==start and start==end, it is supported.
678 void DataArray::CheckValueInRangeEx(mcIdType value, mcIdType start, mcIdType end, const std::string& msg)
680 if(start<0 || start>=value)
682 if(value!=start || end!=start)
684 std::ostringstream oss; oss << "DataArray::CheckValueInRangeEx : " << msg << " ! Expected start " << start << " of input range, in [0," << value << "[ !";
685 throw INTERP_KERNEL::Exception(oss.str().c_str());
688 if(end<0 || end>value)
690 std::ostringstream oss; oss << "DataArray::CheckValueInRangeEx : " << msg << " ! Expected end " << end << " of input range, in [0," << value << "] !";
691 throw INTERP_KERNEL::Exception(oss.str().c_str());
695 void DataArray::CheckClosingParInRange(mcIdType ref, mcIdType value, const std::string& msg)
697 if(value<0 || value>ref)
699 std::ostringstream oss; oss << "DataArray::CheckClosingParInRange : " << msg << " ! Expected input range in [0," << ref << "] having closing open parenthesis " << value << " !";
700 throw INTERP_KERNEL::Exception(oss.str().c_str());
705 * This method is useful to slice work among a pool of threads or processes. \a begin, \a end \a step is the input whole slice of work to perform,
706 * typically it is a whole slice of tuples of DataArray or cells, nodes of a mesh...
708 * The input \a sliceId should be an id in [0, \a nbOfSlices) that specifies the slice of work.
710 * \param [in] start - the start of the input slice of the whole work to perform split into slices.
711 * \param [in] stop - the stop of the input slice of the whole work to perform split into slices.
712 * \param [in] step - the step (that can be <0) of the input slice of the whole work to perform split into slices.
713 * \param [in] sliceId - the slice id considered
714 * \param [in] nbOfSlices - the number of slices (typically the number of cores on which the work is expected to be sliced)
715 * \param [out] startSlice - the start of the slice considered
716 * \param [out] stopSlice - the stop of the slice consided
718 * \throw If \a step == 0
719 * \throw If \a nbOfSlices not > 0
720 * \throw If \a sliceId not in [0,nbOfSlices)
722 void DataArray::GetSlice(mcIdType start, mcIdType stop, mcIdType step, mcIdType sliceId, mcIdType nbOfSlices, mcIdType& startSlice, mcIdType& stopSlice)
724 DataArrayTools<mcIdType>::GetSlice(start, stop, step, sliceId, nbOfSlices, startSlice, stopSlice);
727 mcIdType DataArray::GetNumberOfItemGivenBES(mcIdType begin, mcIdType end, mcIdType step, const std::string& msg)
729 return DataArrayTools<mcIdType>::GetNumberOfItemGivenBES(begin, end, step, msg);
732 mcIdType DataArray::GetNumberOfItemGivenBESRelative(mcIdType begin, mcIdType end, mcIdType step, const std::string& msg)
734 return DataArrayTools<mcIdType>::GetNumberOfItemGivenBESRelative(begin, end, step, msg);
737 mcIdType DataArray::GetPosOfItemGivenBESRelativeNoThrow(mcIdType value, mcIdType begin, mcIdType end, mcIdType step)
739 return DataArrayTools<mcIdType>::GetPosOfItemGivenBESRelativeNoThrow(value, begin, end, step);
743 * Returns a new instance of DataArrayDouble. The caller is to delete this array
744 * using decrRef() as it is no more needed.
746 DataArrayDouble *DataArrayDouble::New()
748 return new DataArrayDouble;
752 * Returns the only one value in \a this, if and only if number of elements
753 * (nb of tuples * nb of components) is equal to 1, and that \a this is allocated.
754 * \return double - the sole value stored in \a this array.
755 * \throw If at least one of conditions stated above is not fulfilled.
757 double DataArrayDouble::doubleValue() const
761 if(getNbOfElems()==1)
763 return *getConstPointer();
766 throw INTERP_KERNEL::Exception("DataArrayDouble::doubleValue : DataArrayDouble instance is allocated but number of elements is not equal to 1 !");
769 throw INTERP_KERNEL::Exception("DataArrayDouble::doubleValue : DataArrayDouble instance is not allocated !");
773 * Returns a full copy of \a this. For more info on copying data arrays see
774 * \ref MEDCouplingArrayBasicsCopyDeep.
775 * \return DataArrayDouble * - a new instance of DataArrayDouble. The caller is to
776 * delete this array using decrRef() as it is no more needed.
778 DataArrayDouble *DataArrayDouble::deepCopy() const
780 return new DataArrayDouble(*this);
784 * Checks that \a this array is consistently **increasing** or **decreasing** in value,
785 * with at least absolute difference value of |\a eps| at each step.
786 * If not an exception is thrown.
787 * \param [in] increasing - if \a true, the array values should be increasing.
788 * \param [in] eps - minimal absolute difference between the neighbor values at which
789 * the values are considered different.
790 * \throw If sequence of values is not strictly monotonic in agreement with \a
792 * \throw If \a this->getNumberOfComponents() != 1.
793 * \throw If \a this is not allocated.
795 void DataArrayDouble::checkMonotonic(bool increasing, double eps) const
797 if(!isMonotonic(increasing,eps))
800 throw INTERP_KERNEL::Exception("DataArrayDouble::checkMonotonic : 'this' is not INCREASING monotonic !");
802 throw INTERP_KERNEL::Exception("DataArrayDouble::checkMonotonic : 'this' is not DECREASING monotonic !");
807 * Checks that \a this array is consistently **increasing** or **decreasing** in value,
808 * with at least absolute difference value of |\a eps| at each step.
809 * \param [in] increasing - if \a true, array values should be increasing.
810 * \param [in] eps - minimal absolute difference between the neighbor values at which
811 * the values are considered different.
812 * \return bool - \a true if values change in accordance with \a increasing arg.
813 * \throw If \a this->getNumberOfComponents() != 1.
814 * \throw If \a this is not allocated.
816 bool DataArrayDouble::isMonotonic(bool increasing, double eps) const
819 if(getNumberOfComponents()!=1)
820 throw INTERP_KERNEL::Exception("DataArrayDouble::isMonotonic : only supported with 'this' array with ONE component !");
821 mcIdType nbOfElements(getNumberOfTuples());
822 const double *ptr=getConstPointer();
826 double absEps=fabs(eps);
829 for(mcIdType i=1;i<nbOfElements;i++)
831 if(ptr[i]<(ref+absEps))
839 for(mcIdType i=1;i<nbOfElements;i++)
841 if(ptr[i]>(ref-absEps))
849 void DataArrayDouble::writeVTK(std::ostream& ofs, mcIdType indent, const std::string& nameInFile, DataArrayByte *byteArr) const
851 static const char SPACE[4]={' ',' ',' ',' '};
853 std::string idt(indent,' ');
855 ofs << idt << "<DataArray type=\"Float32\" Name=\"" << nameInFile << "\" NumberOfComponents=\"" << getNumberOfComponents() << "\"";
857 bool areAllEmpty(true);
858 for(std::vector<std::string>::const_iterator it=_info_on_compo.begin();it!=_info_on_compo.end();it++)
862 for(std::size_t i=0;i<_info_on_compo.size();i++)
863 ofs << " ComponentName" << i << "=\"" << _info_on_compo[i] << "\"";
867 ofs << " format=\"appended\" offset=\"" << byteArr->getNumberOfTuples() << "\">";
868 INTERP_KERNEL::AutoPtr<float> tmp(new float[getNbOfElems()]);
870 // to make Visual C++ happy : instead of std::copy(begin(),end(),(float *)tmp);
871 for(const double *src=begin();src!=end();src++,pt++)
873 const char *data(reinterpret_cast<const char *>((float *)tmp));
874 std::size_t sz(getNbOfElems()*sizeof(float));
875 byteArr->insertAtTheEnd(data,data+sz);
876 byteArr->insertAtTheEnd(SPACE,SPACE+4);
880 ofs << " RangeMin=\"" << getMinValueInArray() << "\" RangeMax=\"" << getMaxValueInArray() << "\" format=\"ascii\">\n" << idt;
881 std::copy(begin(),end(),std::ostream_iterator<double>(ofs," "));
883 ofs << std::endl << idt << "</DataArray>\n";
886 void DataArrayDouble::reprCppStream(const std::string& varName, std::ostream& stream) const
888 mcIdType nbTuples=getNumberOfTuples();
889 std::size_t nbComp=getNumberOfComponents();
890 const double *data(getConstPointer());
891 stream.precision(17);
892 stream << "DataArrayDouble *" << varName << "=DataArrayDouble::New();" << std::endl;
893 if(nbTuples*nbComp>=1)
895 stream << "const double " << varName << "Data[" << nbTuples*nbComp << "]={";
896 std::copy(data,data+nbTuples*nbComp-1,std::ostream_iterator<double>(stream,","));
897 stream << data[nbTuples*nbComp-1] << "};" << std::endl;
898 stream << varName << "->useArray(" << varName << "Data,false,CPP_DEALLOC," << nbTuples << "," << nbComp << ");" << std::endl;
901 stream << varName << "->alloc(" << nbTuples << "," << nbComp << ");" << std::endl;
902 stream << varName << "->setName(\"" << getName() << "\");" << std::endl;
906 * Method that gives a quick overvien of \a this for python.
908 void DataArrayDouble::reprQuickOverview(std::ostream& stream) const
910 static const std::size_t MAX_NB_OF_BYTE_IN_REPR=300;
911 stream << "DataArrayDouble C++ instance at " << this << ". ";
914 std::size_t nbOfCompo(_info_on_compo.size());
917 mcIdType nbOfTuples(getNumberOfTuples());
918 stream << "Number of tuples : " << nbOfTuples << ". Number of components : " << nbOfCompo << "." << std::endl;
919 reprQuickOverviewData(stream,MAX_NB_OF_BYTE_IN_REPR);
922 stream << "Number of components : 0.";
925 stream << "*** No data allocated ****";
928 void DataArrayDouble::reprQuickOverviewData(std::ostream& stream, std::size_t maxNbOfByteInRepr) const
930 const double *data=begin();
931 mcIdType nbOfTuples(getNumberOfTuples());
932 std::size_t nbOfCompo(_info_on_compo.size());
933 std::ostringstream oss2; oss2 << "[";
935 std::string oss2Str(oss2.str());
936 bool isFinished=true;
937 for(mcIdType i=0;i<nbOfTuples && isFinished;i++)
942 for(std::size_t j=0;j<nbOfCompo;j++,data++)
945 if(j!=nbOfCompo-1) oss2 << ", ";
951 if(i!=nbOfTuples-1) oss2 << ", ";
952 std::string oss3Str(oss2.str());
953 if(oss3Str.length()<maxNbOfByteInRepr)
965 * Equivalent to DataArrayDouble::isEqual except that if false the reason of
968 * \param [in] other the instance to be compared with \a this
969 * \param [in] prec the precision to compare numeric data of the arrays.
970 * \param [out] reason In case of inequality returns the reason.
971 * \sa DataArrayDouble::isEqual
973 bool DataArrayDouble::isEqualIfNotWhy(const DataArrayDouble& other, double prec, std::string& reason) const
975 if(!areInfoEqualsIfNotWhy(other,reason))
977 return _mem.isEqual(other._mem,prec,reason);
981 * Checks if \a this and another DataArrayDouble are fully equal. For more info see
982 * \ref MEDCouplingArrayBasicsCompare.
983 * \param [in] other - an instance of DataArrayDouble to compare with \a this one.
984 * \param [in] prec - precision value to compare numeric data of the arrays.
985 * \return bool - \a true if the two arrays are equal, \a false else.
987 bool DataArrayDouble::isEqual(const DataArrayDouble& other, double prec) const
990 return isEqualIfNotWhy(other,prec,tmp);
994 * Checks if values of \a this and another DataArrayDouble are equal. For more info see
995 * \ref MEDCouplingArrayBasicsCompare.
996 * \param [in] other - an instance of DataArrayDouble to compare with \a this one.
997 * \param [in] prec - precision value to compare numeric data of the arrays.
998 * \return bool - \a true if the values of two arrays are equal, \a false else.
1000 bool DataArrayDouble::isEqualWithoutConsideringStr(const DataArrayDouble& other, double prec) const
1003 return _mem.isEqual(other._mem,prec,tmp);
1007 * This method checks that all tuples in \a other are in \a this.
1008 * If true, the output param \a tupleIds contains the tuples ids of \a this that correspond to tupes in \a this.
1009 * For each i in [ 0 , other->getNumberOfTuples() ) tuple #i in \a other is equal ( regarding input precision \a prec ) to tuple tupleIds[i] in \a this.
1011 * \param [in] other - the array having the same number of components than \a this.
1012 * \param [out] tupleIds - the tuple ids containing the same number of tuples than \a other has.
1013 * \sa DataArrayDouble::findCommonTuples
1015 bool DataArrayDouble::areIncludedInMe(const DataArrayDouble *other, double prec, DataArrayIdType *&tupleIds) const
1018 throw INTERP_KERNEL::Exception("DataArrayDouble::areIncludedInMe : input array is NULL !");
1019 checkAllocated(); other->checkAllocated();
1020 if(getNumberOfComponents()!=other->getNumberOfComponents())
1021 throw INTERP_KERNEL::Exception("DataArrayDouble::areIncludedInMe : the number of components does not match !");
1022 MCAuto<DataArrayDouble> a=DataArrayDouble::Aggregate(this,other);
1023 DataArrayIdType *c=0,*ci=0;
1024 a->findCommonTuples(prec,getNumberOfTuples(),c,ci);
1025 MCAuto<DataArrayIdType> cSafe(c),ciSafe(ci);
1026 mcIdType newNbOfTuples=-1;
1027 MCAuto<DataArrayIdType> ids=DataArrayIdType::ConvertIndexArrayToO2N(a->getNumberOfTuples(),c->begin(),ci->begin(),ci->end(),newNbOfTuples);
1028 MCAuto<DataArrayIdType> ret1=ids->selectByTupleIdSafeSlice(getNumberOfTuples(),a->getNumberOfTuples(),1);
1029 tupleIds=ret1.retn();
1030 return newNbOfTuples==getNumberOfTuples();
1034 * Searches for tuples coincident within \a prec tolerance. Each tuple is considered
1035 * as coordinates of a point in getNumberOfComponents()-dimensional space. The
1036 * distance separating two points is computed with the infinite norm.
1038 * Indices of coincident tuples are stored in output arrays.
1039 * A pair of arrays (\a comm, \a commIndex) is called "Surjective Format 2".
1041 * This method is typically used by MEDCouplingPointSet::findCommonNodes() and
1042 * MEDCouplingUMesh::mergeNodes().
1043 * \param [in] prec - minimal absolute distance between two tuples (infinite norm) at which they are
1044 * considered not coincident.
1045 * \param [in] limitTupleId - limit tuple id. If all tuples within a group of coincident
1046 * tuples have id strictly lower than \a limitTupleId then they are not returned.
1047 * \param [out] comm - the array holding ids (== indices) of coincident tuples.
1048 * \a comm->getNumberOfComponents() == 1.
1049 * \a comm->getNumberOfTuples() == \a commIndex->back().
1050 * \param [out] commIndex - the array dividing all indices stored in \a comm into
1051 * groups of (indices of) coincident tuples. Its every value is a tuple
1052 * index where a next group of tuples begins. For example the second
1053 * group of tuples in \a comm is described by following range of indices:
1054 * [ \a commIndex[1], \a commIndex[2] ). \a commIndex->getNumberOfTuples()-1
1055 * gives the number of groups of coincident tuples.
1056 * \throw If \a this is not allocated.
1057 * \throw If the number of components is not in [1,2,3,4].
1059 * \if ENABLE_EXAMPLES
1060 * \ref cpp_mcdataarraydouble_findcommontuples "Here is a C++ example".
1062 * \ref py_mcdataarraydouble_findcommontuples "Here is a Python example".
1064 * \sa DataArrayInt::ConvertIndexArrayToO2N(), DataArrayDouble::areIncludedInMe
1066 void DataArrayDouble::findCommonTuples(double prec, mcIdType limitTupleId, DataArrayIdType *&comm, DataArrayIdType *&commIndex) const
1069 std::size_t nbOfCompo=getNumberOfComponents();
1070 if ((nbOfCompo<1) || (nbOfCompo>4)) //test before work
1071 throw INTERP_KERNEL::Exception("DataArrayDouble::findCommonTuples : Unexpected spacedim of coords. Must be 1, 2, 3 or 4.");
1073 mcIdType nbOfTuples(getNumberOfTuples());
1075 MCAuto<DataArrayIdType> c(DataArrayIdType::New()),cI(DataArrayIdType::New()); c->alloc(0,1); cI->pushBackSilent(0);
1079 findCommonTuplesAlg<4>(begin(),nbOfTuples,limitTupleId,prec,c,cI);
1082 findCommonTuplesAlg<3>(begin(),nbOfTuples,limitTupleId,prec,c,cI);
1085 findCommonTuplesAlg<2>(begin(),nbOfTuples,limitTupleId,prec,c,cI);
1088 findCommonTuplesAlg<1>(begin(),nbOfTuples,limitTupleId,prec,c,cI);
1091 throw INTERP_KERNEL::Exception("DataArrayDouble::findCommonTuples : nb of components managed are 1,2,3 and 4 ! not implemented for other number of components !");
1094 commIndex=cI.retn();
1098 * This methods returns the minimal distance between the two set of points \a this and \a other.
1099 * So \a this and \a other have to have the same number of components. If not an INTERP_KERNEL::Exception will be thrown.
1100 * This method works only if number of components of \a this (equal to those of \a other) is in 1, 2 or 3.
1102 * \param [out] thisTupleId the tuple id in \a this corresponding to the returned minimal distance
1103 * \param [out] otherTupleId the tuple id in \a other corresponding to the returned minimal distance
1104 * \return the minimal distance between the two set of points \a this and \a other.
1105 * \sa DataArrayDouble::findClosestTupleId
1107 double DataArrayDouble::minimalDistanceTo(const DataArrayDouble *other, mcIdType& thisTupleId, mcIdType& otherTupleId) const
1109 MCAuto<DataArrayIdType> part1=findClosestTupleId(other);
1110 std::size_t nbOfCompo=getNumberOfComponents();
1111 mcIdType otherNbTuples=other->getNumberOfTuples();
1112 const double *thisPt(begin()),*otherPt(other->begin());
1113 const mcIdType *part1Pt(part1->begin());
1114 double ret=std::numeric_limits<double>::max();
1115 for(mcIdType i=0;i<otherNbTuples;i++,part1Pt++,otherPt+=nbOfCompo)
1118 for(std::size_t j=0;j<nbOfCompo;j++)
1119 tmp+=(otherPt[j]-thisPt[nbOfCompo*(*part1Pt)+j])*(otherPt[j]-thisPt[nbOfCompo*(*part1Pt)+j]);
1121 { ret=tmp; thisTupleId=*part1Pt; otherTupleId=i; }
1127 * This methods returns for each tuple in \a other which tuple in \a this is the closest.
1128 * So \a this and \a other have to have the same number of components. If not an INTERP_KERNEL::Exception will be thrown.
1129 * This method works only if number of components of \a this (equal to those of \a other) is in 1, 2 or 3.
1131 * \return a newly allocated (new object to be dealt by the caller) DataArrayInt having \c other->getNumberOfTuples() tuples and one components.
1132 * \sa DataArrayDouble::minimalDistanceTo
1134 DataArrayIdType *DataArrayDouble::findClosestTupleId(const DataArrayDouble *other) const
1137 throw INTERP_KERNEL::Exception("DataArrayDouble::findClosestTupleId : other instance is NULL !");
1138 checkAllocated(); other->checkAllocated();
1139 std::size_t nbOfCompo(getNumberOfComponents());
1140 if(nbOfCompo!=other->getNumberOfComponents())
1142 std::ostringstream oss; oss << "DataArrayDouble::findClosestTupleId : number of components in this is " << nbOfCompo;
1143 oss << ", whereas number of components in other is " << other->getNumberOfComponents() << "! Should be equal !";
1144 throw INTERP_KERNEL::Exception(oss.str().c_str());
1146 mcIdType nbOfTuples(other->getNumberOfTuples());
1147 mcIdType thisNbOfTuples(getNumberOfTuples());
1148 MCAuto<DataArrayIdType> ret=DataArrayIdType::New(); ret->alloc(nbOfTuples,1);
1150 getMinMaxPerComponent(bounds);
1155 double xDelta(fabs(bounds[1]-bounds[0])),yDelta(fabs(bounds[3]-bounds[2])),zDelta(fabs(bounds[5]-bounds[4]));
1156 double delta=std::max(xDelta,yDelta); delta=std::max(delta,zDelta);
1157 double characSize=pow((delta*delta*delta)/((double)thisNbOfTuples),1./3.);
1158 BBTreePts<3,mcIdType> myTree(begin(),0,0,getNumberOfTuples(),characSize*1e-12);
1159 FindClosestTupleIdAlg<3>(myTree,3.*characSize*characSize,other->begin(),nbOfTuples,begin(),thisNbOfTuples,ret->getPointer());
1164 double xDelta(fabs(bounds[1]-bounds[0])),yDelta(fabs(bounds[3]-bounds[2]));
1165 double delta=std::max(xDelta,yDelta);
1166 double characSize=sqrt(delta/(double)thisNbOfTuples);
1167 BBTreePts<2,mcIdType> myTree(begin(),0,0,getNumberOfTuples(),characSize*1e-12);
1168 FindClosestTupleIdAlg<2>(myTree,2.*characSize*characSize,other->begin(),nbOfTuples,begin(),thisNbOfTuples,ret->getPointer());
1173 double characSize=fabs(bounds[1]-bounds[0])/FromIdType<double>(thisNbOfTuples);
1174 BBTreePts<1,mcIdType> myTree(begin(),0,0,getNumberOfTuples(),characSize*1e-12);
1175 FindClosestTupleIdAlg<1>(myTree,1.*characSize*characSize,other->begin(),nbOfTuples,begin(),thisNbOfTuples,ret->getPointer());
1179 throw INTERP_KERNEL::Exception("Unexpected spacedim of coords for findClosestTupleId. Must be 1, 2 or 3.");
1185 * This method expects that \a this and \a otherBBoxFrmt arrays are bounding box arrays ( as the output of MEDCouplingPointSet::getBoundingBoxForBBTree method ).
1186 * This method will return a DataArrayInt array having the same number of tuples than \a this. This returned array tells for each cell in \a this
1187 * how many bounding boxes in \a otherBBoxFrmt.
1188 * So, this method expects that \a this and \a otherBBoxFrmt have the same number of components.
1190 * \param [in] otherBBoxFrmt - It is an array .
1191 * \param [in] eps - the absolute precision of the detection. when eps < 0 the bboxes are enlarged so more interactions are detected. Inversely when > 0 the bboxes are stretched.
1192 * \sa MEDCouplingPointSet::getBoundingBoxForBBTree
1193 * \throw If \a this and \a otherBBoxFrmt have not the same number of components.
1194 * \throw If \a this and \a otherBBoxFrmt number of components is not even (BBox format).
1196 DataArrayIdType *DataArrayDouble::computeNbOfInteractionsWith(const DataArrayDouble *otherBBoxFrmt, double eps) const
1199 throw INTERP_KERNEL::Exception("DataArrayDouble::computeNbOfInteractionsWith : input array is NULL !");
1200 if(!isAllocated() || !otherBBoxFrmt->isAllocated())
1201 throw INTERP_KERNEL::Exception("DataArrayDouble::computeNbOfInteractionsWith : this and input array must be allocated !");
1202 std::size_t nbOfComp(getNumberOfComponents());
1203 mcIdType nbOfTuples(getNumberOfTuples());
1204 if(nbOfComp!=otherBBoxFrmt->getNumberOfComponents())
1206 std::ostringstream oss; oss << "DataArrayDouble::computeNbOfInteractionsWith : this number of components (" << nbOfComp << ") must be equal to the number of components of input array (" << otherBBoxFrmt->getNumberOfComponents() << ") !";
1207 throw INTERP_KERNEL::Exception(oss.str().c_str());
1211 std::ostringstream oss; oss << "DataArrayDouble::computeNbOfInteractionsWith : Number of components (" << nbOfComp << ") is not even ! It should be to be compatible with bbox format !";
1212 throw INTERP_KERNEL::Exception(oss.str().c_str());
1214 MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(nbOfTuples,1);
1215 const double *thisBBPtr(begin());
1216 mcIdType *retPtr(ret->getPointer());
1221 BBTree<3,mcIdType> bbt(otherBBoxFrmt->begin(),0,0,otherBBoxFrmt->getNumberOfTuples(),eps);
1222 for(mcIdType i=0;i<nbOfTuples;i++,retPtr++,thisBBPtr+=nbOfComp)
1223 *retPtr=bbt.getNbOfIntersectingElems(thisBBPtr);
1228 BBTree<2,mcIdType> bbt(otherBBoxFrmt->begin(),0,0,otherBBoxFrmt->getNumberOfTuples(),eps);
1229 for(mcIdType i=0;i<nbOfTuples;i++,retPtr++,thisBBPtr+=nbOfComp)
1230 *retPtr=bbt.getNbOfIntersectingElems(thisBBPtr);
1235 BBTree<1,mcIdType> bbt(otherBBoxFrmt->begin(),0,0,otherBBoxFrmt->getNumberOfTuples(),eps);
1236 for(mcIdType i=0;i<nbOfTuples;i++,retPtr++,thisBBPtr+=nbOfComp)
1237 *retPtr=bbt.getNbOfIntersectingElems(thisBBPtr);
1241 throw INTERP_KERNEL::Exception("DataArrayDouble::computeNbOfInteractionsWith : space dimension supported are [1,2,3] !");
1248 * Returns a copy of \a this array by excluding coincident tuples. Each tuple is
1249 * considered as coordinates of a point in getNumberOfComponents()-dimensional
1250 * space. The distance between tuples is computed using norm2. If several tuples are
1251 * not far each from other than \a prec, only one of them remains in the result
1252 * array. The order of tuples in the result array is same as in \a this one except
1253 * that coincident tuples are excluded.
1254 * \param [in] prec - minimal absolute distance between two tuples at which they are
1255 * considered not coincident.
1256 * \param [in] limitTupleId - limit tuple id. If all tuples within a group of coincident
1257 * tuples have id strictly lower than \a limitTupleId then they are not excluded.
1258 * \return DataArrayDouble * - the new instance of DataArrayDouble that the caller
1259 * is to delete using decrRef() as it is no more needed.
1260 * \throw If \a this is not allocated.
1261 * \throw If the number of components is not in [1,2,3,4].
1263 * \if ENABLE_EXAMPLES
1264 * \ref py_mcdataarraydouble_getdifferentvalues "Here is a Python example".
1267 DataArrayDouble *DataArrayDouble::getDifferentValues(double prec, mcIdType limitTupleId) const
1270 DataArrayIdType *c0=0,*cI0=0;
1271 findCommonTuples(prec,limitTupleId,c0,cI0);
1272 MCAuto<DataArrayIdType> c(c0),cI(cI0);
1273 mcIdType newNbOfTuples=-1;
1274 MCAuto<DataArrayIdType> o2n=DataArrayIdType::ConvertIndexArrayToO2N(getNumberOfTuples(),c0->begin(),cI0->begin(),cI0->end(),newNbOfTuples);
1275 return renumberAndReduce(o2n->getConstPointer(),newNbOfTuples);
1279 * Copy all components in a specified order from another DataArrayDouble.
1280 * Both numerical and textual data is copied. The number of tuples in \a this and
1281 * the other array can be different.
1282 * \param [in] a - the array to copy data from.
1283 * \param [in] compoIds - sequence of zero based indices of components, data of which is
1285 * \throw If \a a is NULL.
1286 * \throw If \a compoIds.size() != \a a->getNumberOfComponents().
1287 * \throw If \a compoIds[i] < 0 or \a compoIds[i] > \a this->getNumberOfComponents().
1289 * \if ENABLE_EXAMPLES
1290 * \ref py_mcdataarraydouble_setselectedcomponents "Here is a Python example".
1293 void DataArrayDouble::setSelectedComponents(const DataArrayDouble *a, const std::vector<std::size_t>& compoIds)
1296 throw INTERP_KERNEL::Exception("DataArrayDouble::setSelectedComponents : input DataArrayDouble is NULL !");
1298 copyPartOfStringInfoFrom2(compoIds,*a);
1299 std::size_t partOfCompoSz=compoIds.size();
1300 std::size_t nbOfCompo=getNumberOfComponents();
1301 mcIdType nbOfTuples=std::min(getNumberOfTuples(),a->getNumberOfTuples());
1302 const double *ac=a->getConstPointer();
1303 double *nc=getPointer();
1304 for(mcIdType i=0;i<nbOfTuples;i++)
1305 for(std::size_t j=0;j<partOfCompoSz;j++,ac++)
1306 nc[nbOfCompo*i+compoIds[j]]=*ac;
1309 * Checks if 0.0 value is present in \a this array. If it is the case, an exception
1311 * \throw If zero is found in \a this array.
1313 void DataArrayDouble::checkNoNullValues() const
1315 const double *tmp=getConstPointer();
1316 mcIdType nbOfElems=getNbOfElems();
1317 const double *where=std::find(tmp,tmp+nbOfElems,0.);
1318 if(where!=tmp+nbOfElems)
1319 throw INTERP_KERNEL::Exception("A value 0.0 have been detected !");
1323 * Computes minimal and maximal value in each component. An output array is filled
1324 * with \c 2 * \a this->getNumberOfComponents() values, so the caller is to allocate
1325 * enough memory before calling this method.
1326 * \param [out] bounds - array of size at least 2 *\a this->getNumberOfComponents().
1327 * It is filled as follows:<br>
1328 * \a bounds[0] = \c min_of_component_0 <br>
1329 * \a bounds[1] = \c max_of_component_0 <br>
1330 * \a bounds[2] = \c min_of_component_1 <br>
1331 * \a bounds[3] = \c max_of_component_1 <br>
1334 void DataArrayDouble::getMinMaxPerComponent(double *bounds) const
1337 std::size_t dim=getNumberOfComponents();
1338 for (std::size_t idim=0; idim<dim; idim++)
1340 bounds[idim*2]=std::numeric_limits<double>::max();
1341 bounds[idim*2+1]=-std::numeric_limits<double>::max();
1343 const double *ptr=getConstPointer();
1344 mcIdType nbOfTuples=getNumberOfTuples();
1345 for(mcIdType i=0;i<nbOfTuples;i++)
1347 for(std::size_t idim=0;idim<dim;idim++)
1349 if(bounds[idim*2]>ptr[i*dim+idim])
1351 bounds[idim*2]=ptr[i*dim+idim];
1353 if(bounds[idim*2+1]<ptr[i*dim+idim])
1355 bounds[idim*2+1]=ptr[i*dim+idim];
1362 * This method retrieves a newly allocated DataArrayDouble instance having same number of tuples than \a this and twice number of components than \a this
1363 * to store both the min and max per component of each tuples.
1364 * \param [in] epsilon the width of the bbox (identical in each direction) - 0.0 by default
1366 * \return a newly created DataArrayDouble instance having \c this->getNumberOfTuples() tuples and 2 * \c this->getNumberOfComponent() components
1368 * \throw If \a this is not allocated yet.
1370 DataArrayDouble *DataArrayDouble::computeBBoxPerTuple(double epsilon) const
1373 const double *dataPtr=getConstPointer();
1374 std::size_t nbOfCompo=getNumberOfComponents();
1375 mcIdType nbTuples=getNumberOfTuples();
1376 MCAuto<DataArrayDouble> bbox=DataArrayDouble::New();
1377 bbox->alloc(nbTuples,2*nbOfCompo);
1378 double *bboxPtr=bbox->getPointer();
1379 for(mcIdType i=0;i<nbTuples;i++)
1381 for(std::size_t j=0;j<nbOfCompo;j++)
1383 bboxPtr[2*nbOfCompo*i+2*j]=dataPtr[nbOfCompo*i+j]-epsilon;
1384 bboxPtr[2*nbOfCompo*i+2*j+1]=dataPtr[nbOfCompo*i+j]+epsilon;
1391 * For each tuples **t** in \a other, this method retrieves tuples in \a this that are equal to **t**.
1392 * Two tuples are considered equal if the euclidian distance between the two tuples is lower than \a eps.
1394 * \param [in] other a DataArrayDouble having same number of components than \a this.
1395 * \param [in] eps absolute precision representing distance (using infinite norm) between 2 tuples behind which 2 tuples are considered equal.
1396 * \param [out] c will contain the set of tuple ids in \a this that are equal to to the tuple ids in \a other contiguously.
1397 * \a cI allows to extract information in \a c.
1398 * \param [out] cI is an indirection array that allows to extract the data contained in \a c.
1400 * \throw In case of:
1401 * - \a this is not allocated
1402 * - \a other is not allocated or null
1403 * - \a this and \a other do not have the same number of components
1404 * - if number of components of \a this is not in [1,2,3]
1406 * \sa MEDCouplingPointSet::getNodeIdsNearPoints, DataArrayDouble::getDifferentValues
1408 void DataArrayDouble::computeTupleIdsNearTuples(const DataArrayDouble *other, double eps, DataArrayIdType *& c, DataArrayIdType *& cI) const
1411 throw INTERP_KERNEL::Exception("DataArrayDouble::computeTupleIdsNearTuples : input pointer other is null !");
1413 other->checkAllocated();
1414 std::size_t nbOfCompo=getNumberOfComponents();
1415 std::size_t otherNbOfCompo=other->getNumberOfComponents();
1416 if(nbOfCompo!=otherNbOfCompo)
1417 throw INTERP_KERNEL::Exception("DataArrayDouble::computeTupleIdsNearTuples : number of components should be equal between this and other !");
1418 mcIdType nbOfTuplesOther=other->getNumberOfTuples();
1419 mcIdType nbOfTuples=getNumberOfTuples();
1420 MCAuto<DataArrayIdType> cArr(DataArrayIdType::New()),cIArr(DataArrayIdType::New()); cArr->alloc(0,1); cIArr->pushBackSilent(0);
1425 BBTreePts<3,mcIdType> myTree(begin(),0,0,nbOfTuples,eps);
1426 FindTupleIdsNearTuplesAlg<3>(myTree,other->getConstPointer(),nbOfTuplesOther,eps,cArr,cIArr);
1431 BBTreePts<2,mcIdType> myTree(begin(),0,0,nbOfTuples,eps);
1432 FindTupleIdsNearTuplesAlg<2>(myTree,other->getConstPointer(),nbOfTuplesOther,eps,cArr,cIArr);
1437 BBTreePts<1,mcIdType> myTree(begin(),0,0,nbOfTuples,eps);
1438 FindTupleIdsNearTuplesAlg<1>(myTree,other->getConstPointer(),nbOfTuplesOther,eps,cArr,cIArr);
1442 throw INTERP_KERNEL::Exception("Unexpected spacedim of coords for computeTupleIdsNearTuples. Must be 1, 2 or 3.");
1444 c=cArr.retn(); cI=cIArr.retn();
1448 * This method recenter tuples in \b this in order to be centered at the origin to benefit about the advantages of maximal precision to be around the box
1449 * around origin of 'radius' 1.
1451 * \param [in] eps absolute epsilon. under that value of delta between max and min no scale is performed.
1453 void DataArrayDouble::recenterForMaxPrecision(double eps)
1456 std::size_t dim=getNumberOfComponents();
1457 std::vector<double> bounds(2*dim);
1458 getMinMaxPerComponent(&bounds[0]);
1459 for(std::size_t i=0;i<dim;i++)
1461 double delta=bounds[2*i+1]-bounds[2*i];
1462 double offset=(bounds[2*i]+bounds[2*i+1])/2.;
1464 applyLin(1./delta,-offset/delta,i);
1466 applyLin(1.,-offset,i);
1471 * Returns the maximal value and all its locations within \a this one-dimensional array.
1472 * \param [out] tupleIds - a new instance of DataArrayInt containing indices of
1473 * tuples holding the maximal value. The caller is to delete it using
1474 * decrRef() as it is no more needed.
1475 * \return double - the maximal value among all values of \a this array.
1476 * \throw If \a this->getNumberOfComponents() != 1
1477 * \throw If \a this->getNumberOfTuples() < 1
1479 double DataArrayDouble::getMaxValue2(DataArrayIdType*& tupleIds) const
1483 double ret=getMaxValue(tmp);
1484 tupleIds=findIdsInRange(ret,ret);
1489 * Returns the minimal value and all its locations within \a this one-dimensional array.
1490 * \param [out] tupleIds - a new instance of DataArrayInt containing indices of
1491 * tuples holding the minimal value. The caller is to delete it using
1492 * decrRef() as it is no more needed.
1493 * \return double - the minimal value among all values of \a this array.
1494 * \throw If \a this->getNumberOfComponents() != 1
1495 * \throw If \a this->getNumberOfTuples() < 1
1497 double DataArrayDouble::getMinValue2(DataArrayIdType*& tupleIds) const
1501 double ret=getMinValue(tmp);
1502 tupleIds=findIdsInRange(ret,ret);
1507 * This method returns the number of values in \a this that are equals ( within an absolute precision of \a eps ) to input parameter \a value.
1508 * This method only works for single component array.
1510 * \return a value in [ 0, \c this->getNumberOfTuples() )
1512 * \throw If \a this is not allocated
1515 mcIdType DataArrayDouble::count(double value, double eps) const
1519 if(getNumberOfComponents()!=1)
1520 throw INTERP_KERNEL::Exception("DataArrayDouble::count : must be applied on DataArrayDouble with only one component, you can call 'rearrange' method before !");
1521 const double *vals=begin();
1522 mcIdType nbOfTuples=getNumberOfTuples();
1523 for(mcIdType i=0;i<nbOfTuples;i++,vals++)
1524 if(fabs(*vals-value)<=eps)
1530 * Returns the average value of \a this one-dimensional array.
1531 * \return double - the average value over all values of \a this array.
1532 * \throw If \a this->getNumberOfComponents() != 1
1533 * \throw If \a this->getNumberOfTuples() < 1
1535 double DataArrayDouble::getAverageValue() const
1537 if(getNumberOfComponents()!=1)
1538 throw INTERP_KERNEL::Exception("DataArrayDouble::getAverageValue : must be applied on DataArrayDouble with only one component, you can call 'rearrange' method before !");
1539 mcIdType nbOfTuples(getNumberOfTuples());
1541 throw INTERP_KERNEL::Exception("DataArrayDouble::getAverageValue : array exists but number of tuples must be > 0 !");
1542 const double *vals=getConstPointer();
1543 double ret=std::accumulate(vals,vals+nbOfTuples,0.);
1544 return ret/FromIdType<double>(nbOfTuples);
1548 * Returns the Euclidean norm of the vector defined by \a this array.
1549 * \return double - the value of the Euclidean norm, i.e.
1550 * the square root of the inner product of vector.
1551 * \throw If \a this is not allocated.
1553 double DataArrayDouble::norm2() const
1557 std::size_t nbOfElems=getNbOfElems();
1558 const double *pt=getConstPointer();
1559 for(std::size_t i=0;i<nbOfElems;i++,pt++)
1565 * Returns the maximum norm of the vector defined by \a this array.
1566 * This method works even if the number of components is different from one.
1567 * If the number of elements in \a this is 0, -1. is returned.
1568 * \return double - the value of the maximum norm, i.e.
1569 * the maximal absolute value among values of \a this array (whatever its number of components).
1570 * \throw If \a this is not allocated.
1572 double DataArrayDouble::normMax() const
1576 std::size_t nbOfElems(getNbOfElems());
1577 const double *pt(getConstPointer());
1578 for(std::size_t i=0;i<nbOfElems;i++,pt++)
1580 double val(std::abs(*pt));
1588 * Returns the maximum norm of for each component of \a this array.
1589 * If the number of elements in \a this is 0, -1. is returned.
1590 * \param [out] res - pointer to an array of result values, of size at least \a
1591 * this->getNumberOfComponents(), that is to be allocated by the caller.
1592 * \throw If \a this is not allocated.
1594 void DataArrayDouble::normMaxPerComponent(double * res) const
1597 mcIdType nbOfTuples(getNumberOfTuples());
1598 std::size_t nbOfCompos(getNumberOfComponents());
1599 std::fill(res, res+nbOfCompos, -1.0);
1600 const double *pt(getConstPointer());
1601 for(mcIdType i=0;i<nbOfTuples;i++)
1602 for (std::size_t j=0; j<nbOfCompos; j++, pt++)
1604 double val(std::abs(*pt));
1612 * Returns the minimum norm (absolute value) of the vector defined by \a this array.
1613 * This method works even if the number of components is different from one.
1614 * If the number of elements in \a this is 0, std::numeric_limits<double>::max() is returned.
1615 * \return double - the value of the minimum norm, i.e.
1616 * the minimal absolute value among values of \a this array (whatever its number of components).
1617 * \throw If \a this is not allocated.
1619 double DataArrayDouble::normMin() const
1622 double ret(std::numeric_limits<double>::max());
1623 std::size_t nbOfElems(getNbOfElems());
1624 const double *pt(getConstPointer());
1625 for(std::size_t i=0;i<nbOfElems;i++,pt++)
1627 double val(std::abs(*pt));
1635 * Accumulates values of each component of \a this array.
1636 * \param [out] res - an array of length \a this->getNumberOfComponents(), allocated
1637 * by the caller, that is filled by this method with sum value for each
1639 * \throw If \a this is not allocated.
1641 void DataArrayDouble::accumulate(double *res) const
1644 const double *ptr=getConstPointer();
1645 mcIdType nbTuple(getNumberOfTuples());
1646 std::size_t nbComps(getNumberOfComponents());
1647 std::fill(res,res+nbComps,0.);
1648 for(mcIdType i=0;i<nbTuple;i++)
1649 std::transform(ptr+i*nbComps,ptr+(i+1)*nbComps,res,res,std::plus<double>());
1653 * This method returns the min distance from an external tuple defined by [ \a tupleBg , \a tupleEnd ) to \a this and
1654 * the first tuple in \a this that matches the returned distance. If there is no tuples in \a this an exception will be thrown.
1657 * \a this is expected to be allocated and expected to have a number of components equal to the distance from \a tupleBg to
1658 * \a tupleEnd. If not an exception will be thrown.
1660 * \param [in] tupleBg start pointer (included) of input external tuple
1661 * \param [in] tupleEnd end pointer (not included) of input external tuple
1662 * \param [out] tupleId the tuple id in \a this that matches the min of distance between \a this and input external tuple
1663 * \return the min distance.
1664 * \sa MEDCouplingUMesh::distanceToPoint
1666 double DataArrayDouble::distanceToTuple(const double *tupleBg, const double *tupleEnd, mcIdType& tupleId) const
1669 mcIdType nbTuple(getNumberOfTuples());
1670 std::size_t nbComps(getNumberOfComponents());
1671 if(nbComps!=(std::size_t)std::distance(tupleBg,tupleEnd))
1672 { std::ostringstream oss; oss << "DataArrayDouble::distanceToTuple : size of input tuple is " << std::distance(tupleBg,tupleEnd) << " should be equal to the number of components in this : " << nbComps << " !"; throw INTERP_KERNEL::Exception(oss.str().c_str()); }
1674 throw INTERP_KERNEL::Exception("DataArrayDouble::distanceToTuple : no tuple in this ! No distance to compute !");
1675 double ret0=std::numeric_limits<double>::max();
1677 const double *work=getConstPointer();
1678 for(mcIdType i=0;i<nbTuple;i++)
1681 for(std::size_t j=0;j<nbComps;j++,work++)
1682 val+=(*work-tupleBg[j])*((*work-tupleBg[j]));
1686 { ret0=val; tupleId=i; }
1692 * Accumulate values of the given component of \a this array.
1693 * \param [in] compId - the index of the component of interest.
1694 * \return double - a sum value of \a compId-th component.
1695 * \throw If \a this is not allocated.
1696 * \throw If \a the condition ( 0 <= \a compId < \a this->getNumberOfComponents() ) is
1699 double DataArrayDouble::accumulate(std::size_t compId) const
1702 const double *ptr=getConstPointer();
1703 mcIdType nbTuple(getNumberOfTuples());
1704 std::size_t nbComps(getNumberOfComponents());
1706 throw INTERP_KERNEL::Exception("DataArrayDouble::accumulate : Invalid compId specified : No such nb of components !");
1708 for(mcIdType i=0;i<nbTuple;i++)
1709 ret+=ptr[i*nbComps+compId];
1714 * This method accumulate using addition tuples in \a this using input index array [ \a bgOfIndex, \a endOfIndex ).
1715 * The returned array will have same number of components than \a this and number of tuples equal to
1716 * \c std::distance(bgOfIndex,endOfIndex) \b minus \b one.
1718 * The input index array is expected to be ascendingly sorted in which the all referenced ids should be in [0, \c this->getNumberOfTuples).
1719 * This method is quite useful for users that need to put a field on cells to field on nodes on the same mesh without a need of conservation.
1721 * \param [in] bgOfIndex - begin (included) of the input index array.
1722 * \param [in] endOfIndex - end (excluded) of the input index array.
1723 * \return DataArrayDouble * - the new instance having the same number of components than \a this.
1725 * \throw If bgOfIndex or end is NULL.
1726 * \throw If input index array is not ascendingly sorted.
1727 * \throw If there is an id in [ \a bgOfIndex, \a endOfIndex ) not in [0, \c this->getNumberOfTuples).
1728 * \throw If std::distance(bgOfIndex,endOfIndex)==0.
1730 DataArrayDouble *DataArrayDouble::accumulatePerChunck(const mcIdType *bgOfIndex, const mcIdType *endOfIndex) const
1732 if(!bgOfIndex || !endOfIndex)
1733 throw INTERP_KERNEL::Exception("DataArrayDouble::accumulatePerChunck : input pointer NULL !");
1735 std::size_t nbCompo(getNumberOfComponents());
1736 mcIdType nbOfTuples(getNumberOfTuples());
1737 std::size_t sz=std::distance(bgOfIndex,endOfIndex);
1739 throw INTERP_KERNEL::Exception("DataArrayDouble::accumulatePerChunck : invalid size of input index array !");
1741 MCAuto<DataArrayDouble> ret=DataArrayDouble::New(); ret->alloc(sz,nbCompo);
1742 const mcIdType *w=bgOfIndex;
1743 if(*w<0 || *w>=nbOfTuples)
1744 throw INTERP_KERNEL::Exception("DataArrayDouble::accumulatePerChunck : The first element of the input index not in [0,nbOfTuples) !");
1745 const double *srcPt=begin()+(*w)*nbCompo;
1746 double *tmp=ret->getPointer();
1747 for(std::size_t i=0;i<sz;i++,tmp+=nbCompo,w++)
1749 std::fill(tmp,tmp+nbCompo,0.);
1752 for(mcIdType j=w[0];j<w[1];j++,srcPt+=nbCompo)
1754 if(j>=0 && j<nbOfTuples)
1755 std::transform(srcPt,srcPt+nbCompo,tmp,tmp,std::plus<double>());
1758 std::ostringstream oss; oss << "DataArrayDouble::accumulatePerChunck : At rank #" << i << " the input index array points to id " << j << " should be in [0," << nbOfTuples << ") !";
1759 throw INTERP_KERNEL::Exception(oss.str().c_str());
1765 std::ostringstream oss; oss << "DataArrayDouble::accumulatePerChunck : At rank #" << i << " the input index array is not in ascendingly sorted.";
1766 throw INTERP_KERNEL::Exception(oss.str().c_str());
1769 ret->copyStringInfoFrom(*this);
1774 * This method is close to numpy cumSum except that number of element is equal to \a this->getNumberOfTuples()+1. First element of DataArray returned is equal to 0.
1775 * This method expects that \a this as only one component. The returned array will have \a this->getNumberOfTuples()+1 tuple with also one component.
1776 * The ith element of returned array is equal to the sum of elements in \a this with rank strictly lower than i.
1778 * \return DataArrayDouble - A newly built array containing cum sum of \a this.
1780 MCAuto<DataArrayDouble> DataArrayDouble::cumSum() const
1783 checkNbOfComps(1,"DataArrayDouble::cumSum : this is expected to be single component");
1784 mcIdType nbOfTuple(getNumberOfTuples());
1785 MCAuto<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(nbOfTuple+1,1);
1786 double *ptr(ret->getPointer());
1788 const double *thisPtr(begin());
1789 for(mcIdType i=0;i<nbOfTuple;i++)
1790 ptr[i+1]=ptr[i]+thisPtr[i];
1795 * Converts each 2D point defined by the tuple of \a this array from the Polar to the
1796 * Cartesian coordinate system. The two components of the tuple of \a this array are
1797 * considered to contain (1) radius and (2) angle of the point in the Polar CS.
1798 * \return DataArrayDouble * - the new instance of DataArrayDouble, whose each tuple
1799 * contains X and Y coordinates of the point in the Cartesian CS. The caller
1800 * is to delete this array using decrRef() as it is no more needed. The array
1801 * does not contain any textual info on components.
1802 * \throw If \a this->getNumberOfComponents() != 2.
1803 * \sa fromCartToPolar
1805 DataArrayDouble *DataArrayDouble::fromPolarToCart() const
1808 std::size_t nbOfComp(getNumberOfComponents());
1810 throw INTERP_KERNEL::Exception("DataArrayDouble::fromPolarToCart : must be an array with exactly 2 components !");
1811 mcIdType nbOfTuple(getNumberOfTuples());
1812 DataArrayDouble *ret(DataArrayDouble::New());
1813 ret->alloc(nbOfTuple,2);
1814 double *w(ret->getPointer());
1815 const double *wIn(getConstPointer());
1816 for(mcIdType i=0;i<nbOfTuple;i++,w+=2,wIn+=2)
1818 w[0]=wIn[0]*cos(wIn[1]);
1819 w[1]=wIn[0]*sin(wIn[1]);
1825 * Converts each 3D point defined by the tuple of \a this array from the Cylindrical to
1826 * the Cartesian coordinate system. The three components of the tuple of \a this array
1827 * are considered to contain (1) radius, (2) azimuth and (3) altitude of the point in
1828 * the Cylindrical CS.
1829 * \return DataArrayDouble * - the new instance of DataArrayDouble, whose each tuple
1830 * contains X, Y and Z coordinates of the point in the Cartesian CS. The info
1831 * on the third component is copied from \a this array. The caller
1832 * is to delete this array using decrRef() as it is no more needed.
1833 * \throw If \a this->getNumberOfComponents() != 3.
1836 DataArrayDouble *DataArrayDouble::fromCylToCart() const
1839 std::size_t nbOfComp(getNumberOfComponents());
1841 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCylToCart : must be an array with exactly 3 components !");
1842 mcIdType nbOfTuple(getNumberOfTuples());
1843 DataArrayDouble *ret(DataArrayDouble::New());
1844 ret->alloc(getNumberOfTuples(),3);
1845 double *w(ret->getPointer());
1846 const double *wIn(getConstPointer());
1847 for(mcIdType i=0;i<nbOfTuple;i++,w+=3,wIn+=3)
1849 w[0]=wIn[0]*cos(wIn[1]);
1850 w[1]=wIn[0]*sin(wIn[1]);
1853 ret->setInfoOnComponent(2,getInfoOnComponent(2));
1858 * Converts each 3D point defined by the tuple of \a this array from the Spherical to
1859 * the Cartesian coordinate system. The three components of the tuple of \a this array
1860 * are considered to contain (1) radius, (2) polar angle and (3) azimuthal angle of the
1861 * point in the Cylindrical CS.
1862 * \return DataArrayDouble * - the new instance of DataArrayDouble, whose each tuple
1863 * contains X, Y and Z coordinates of the point in the Cartesian CS. The info
1864 * on the third component is copied from \a this array. The caller
1865 * is to delete this array using decrRef() as it is no more needed.
1866 * \throw If \a this->getNumberOfComponents() != 3.
1867 * \sa fromCartToSpher
1869 DataArrayDouble *DataArrayDouble::fromSpherToCart() const
1872 std::size_t nbOfComp(getNumberOfComponents());
1874 throw INTERP_KERNEL::Exception("DataArrayDouble::fromSpherToCart : must be an array with exactly 3 components !");
1875 mcIdType nbOfTuple(getNumberOfTuples());
1876 DataArrayDouble *ret(DataArrayDouble::New());
1877 ret->alloc(getNumberOfTuples(),3);
1878 double *w(ret->getPointer());
1879 const double *wIn(getConstPointer());
1880 for(mcIdType i=0;i<nbOfTuple;i++,w+=3,wIn+=3)
1882 w[0]=wIn[0]*cos(wIn[2])*sin(wIn[1]);
1883 w[1]=wIn[0]*sin(wIn[2])*sin(wIn[1]);
1884 w[2]=wIn[0]*cos(wIn[1]);
1890 * This method returns a new array containing the same number of tuples than \a this. To do this, this method needs \a at parameter to specify the convention of \a this.
1891 * All the tuples of the returned array will be in cartesian sense. So if \a at equals to AX_CART the returned array is basically a deep copy of \a this.
1892 * If \a at equals to AX_CYL the returned array will be the result of operation cylindric to cartesian of \a this...
1894 * \param [in] atOfThis - The axis type of \a this.
1895 * \return DataArrayDouble * - the new instance of DataArrayDouble (that must be dealed by caller) containing the result of the cartesianizification of \a this.
1897 DataArrayDouble *DataArrayDouble::cartesianize(MEDCouplingAxisType atOfThis) const
1900 std::size_t nbOfComp(getNumberOfComponents());
1901 MCAuto<DataArrayDouble> ret;
1910 ret=fromCylToCart();
1915 ret=fromPolarToCart();
1919 throw INTERP_KERNEL::Exception("DataArrayDouble::cartesianize : For AX_CYL, number of components must be in [2,3] !");
1923 ret=fromSpherToCart();
1928 ret=fromPolarToCart();
1932 throw INTERP_KERNEL::Exception("DataArrayDouble::cartesianize : For AX_CYL, number of components must be in [2,3] !");
1934 throw INTERP_KERNEL::Exception("DataArrayDouble::cartesianize : not recognized axis type ! Only AX_CART, AX_CYL and AX_SPHER supported !");
1936 ret->copyStringInfoFrom(*this);
1941 * This method returns a newly created array to be deallocated that contains the result of conversion from cartesian to polar.
1942 * This method expects that \a this has exactly 2 components.
1943 * \sa fromPolarToCart
1945 DataArrayDouble *DataArrayDouble::fromCartToPolar() const
1947 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
1949 std::size_t nbOfComp(getNumberOfComponents());
1950 mcIdType nbTuples(getNumberOfTuples());
1952 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToPolar : must be an array with exactly 2 components !");
1953 ret->alloc(nbTuples,2);
1954 double *retPtr(ret->getPointer());
1955 const double *ptr(begin());
1956 for(mcIdType i=0;i<nbTuples;i++,ptr+=2,retPtr+=2)
1958 retPtr[0]=sqrt(ptr[0]*ptr[0]+ptr[1]*ptr[1]);
1959 retPtr[1]=atan2(ptr[1],ptr[0]);
1965 * This method returns a newly created array to be deallocated that contains the result of conversion from cartesian to cylindrical.
1966 * This method expects that \a this has exactly 3 components.
1969 DataArrayDouble *DataArrayDouble::fromCartToCyl() const
1971 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
1973 std::size_t nbOfComp(getNumberOfComponents());
1974 mcIdType nbTuples(getNumberOfTuples());
1976 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCyl : must be an array with exactly 3 components !");
1977 ret->alloc(nbTuples,3);
1978 double *retPtr(ret->getPointer());
1979 const double *ptr(begin());
1980 for(mcIdType i=0;i<nbTuples;i++,ptr+=3,retPtr+=3)
1982 retPtr[0]=sqrt(ptr[0]*ptr[0]+ptr[1]*ptr[1]);
1983 retPtr[1]=atan2(ptr[1],ptr[0]);
1990 * This method returns a newly created array to be deallocated that contains the result of conversion from cartesian to spherical coordinates.
1991 * \sa fromSpherToCart
1993 DataArrayDouble *DataArrayDouble::fromCartToSpher() const
1995 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
1997 std::size_t nbOfComp(getNumberOfComponents());
1998 mcIdType nbTuples(getNumberOfTuples());
2000 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToSpher : must be an array with exactly 3 components !");
2001 ret->alloc(nbTuples,3);
2002 double *retPtr(ret->getPointer());
2003 const double *ptr(begin());
2004 for(mcIdType i=0;i<nbTuples;i++,ptr+=3,retPtr+=3)
2006 retPtr[0]=sqrt(ptr[0]*ptr[0]+ptr[1]*ptr[1]+ptr[2]*ptr[2]);
2007 retPtr[1]=acos(ptr[2]/retPtr[0]);
2008 retPtr[2]=atan2(ptr[1],ptr[0]);
2014 * This method returns a newly created array to be deallocated that contains the result of conversion from cartesian to cylindrical relative to the given \a center and a \a vector.
2015 * This method expects that \a this has exactly 3 components.
2016 * \sa MEDCouplingFieldDouble::computeVectorFieldCyl
2018 DataArrayDouble *DataArrayDouble::fromCartToCylGiven(const DataArrayDouble *coords, const double center[3], const double vect[3]) const
2021 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCylGiven : input coords are NULL !");
2022 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
2023 checkAllocated(); coords->checkAllocated();
2024 std::size_t nbOfComp(getNumberOfComponents());
2025 mcIdType nbTuples(getNumberOfTuples());
2027 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCylGiven : must be an array with exactly 3 components !");
2028 if(coords->getNumberOfComponents()!=3)
2029 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCylGiven : coords array must have exactly 3 components !");
2030 if(coords->getNumberOfTuples()!=nbTuples)
2031 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCylGiven : coords array must have the same number of tuples !");
2032 ret->alloc(nbTuples,nbOfComp);
2033 double magOfVect(sqrt(vect[0]*vect[0]+vect[1]*vect[1]+vect[2]*vect[2]));
2035 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCylGiven : magnitude of vect is too low !");
2036 double Ur[3],Uteta[3],Uz[3],*retPtr(ret->getPointer());
2037 const double *coo(coords->begin()),*vectField(begin());
2038 std::transform(vect,vect+3,Uz,std::bind(std::multiplies<double>(),std::placeholders::_1,1./magOfVect));
2039 for(mcIdType i=0;i<nbTuples;i++,vectField+=3,retPtr+=3,coo+=3)
2041 std::transform(coo,coo+3,center,Ur,std::minus<double>());
2042 Uteta[0]=Uz[1]*Ur[2]-Uz[2]*Ur[1]; Uteta[1]=Uz[2]*Ur[0]-Uz[0]*Ur[2]; Uteta[2]=Uz[0]*Ur[1]-Uz[1]*Ur[0];
2043 double magOfTeta(sqrt(Uteta[0]*Uteta[0]+Uteta[1]*Uteta[1]+Uteta[2]*Uteta[2]));
2044 std::transform(Uteta,Uteta+3,Uteta,std::bind(std::multiplies<double>(),std::placeholders::_1,1./magOfTeta));
2045 Ur[0]=Uteta[1]*Uz[2]-Uteta[2]*Uz[1]; Ur[1]=Uteta[2]*Uz[0]-Uteta[0]*Uz[2]; Ur[2]=Uteta[0]*Uz[1]-Uteta[1]*Uz[0];
2046 retPtr[0]=Ur[0]*vectField[0]+Ur[1]*vectField[1]+Ur[2]*vectField[2];
2047 retPtr[1]=Uteta[0]*vectField[0]+Uteta[1]*vectField[1]+Uteta[2]*vectField[2];
2048 retPtr[2]=Uz[0]*vectField[0]+Uz[1]*vectField[1]+Uz[2]*vectField[2];
2050 ret->copyStringInfoFrom(*this);
2055 * Computes the doubly contracted product of every tensor defined by the tuple of \a this
2056 * array containing 6 components.
2057 * \return DataArrayDouble * - the new instance of DataArrayDouble, whose each tuple
2058 * is calculated from the tuple <em>(t)</em> of \a this array as follows:
2059 * \f$ t[0]^2+t[1]^2+t[2]^2+2*t[3]^2+2*t[4]^2+2*t[5]^2\f$.
2060 * The caller is to delete this result array using decrRef() as it is no more needed.
2061 * \throw If \a this->getNumberOfComponents() != 6.
2063 DataArrayDouble *DataArrayDouble::doublyContractedProduct() const
2066 std::size_t nbOfComp(getNumberOfComponents());
2068 throw INTERP_KERNEL::Exception("DataArrayDouble::doublyContractedProduct : must be an array with exactly 6 components !");
2069 DataArrayDouble *ret=DataArrayDouble::New();
2070 mcIdType nbOfTuple=getNumberOfTuples();
2071 ret->alloc(nbOfTuple,1);
2072 const double *src=getConstPointer();
2073 double *dest=ret->getPointer();
2074 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=6)
2075 *dest=src[0]*src[0]+src[1]*src[1]+src[2]*src[2]+2.*src[3]*src[3]+2.*src[4]*src[4]+2.*src[5]*src[5];
2080 * Computes the determinant of every square matrix defined by the tuple of \a this
2081 * array, which contains either 4, 6 or 9 components. The case of 6 components
2082 * corresponds to that of the upper triangular matrix.
2083 * \return DataArrayDouble * - the new instance of DataArrayDouble, whose each tuple
2084 * is the determinant of matrix of the corresponding tuple of \a this array.
2085 * The caller is to delete this result array using decrRef() as it is no more
2087 * \throw If \a this->getNumberOfComponents() is not in [4,6,9].
2089 DataArrayDouble *DataArrayDouble::determinant() const
2092 DataArrayDouble *ret=DataArrayDouble::New();
2093 mcIdType nbOfTuple=getNumberOfTuples();
2094 ret->alloc(nbOfTuple,1);
2095 const double *src=getConstPointer();
2096 double *dest=ret->getPointer();
2097 switch(getNumberOfComponents())
2100 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=6)
2101 *dest=src[0]*src[1]*src[2]+2.*src[4]*src[5]*src[3]-src[0]*src[4]*src[4]-src[2]*src[3]*src[3]-src[1]*src[5]*src[5];
2104 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=4)
2105 *dest=src[0]*src[3]-src[1]*src[2];
2108 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=9)
2109 *dest=src[0]*src[4]*src[8]+src[1]*src[5]*src[6]+src[2]*src[3]*src[7]-src[0]*src[5]*src[7]-src[1]*src[3]*src[8]-src[2]*src[4]*src[6];
2113 throw INTERP_KERNEL::Exception("DataArrayDouble::determinant : Invalid number of components ! must be in 4,6,9 !");
2118 * Computes 3 eigenvalues of every upper triangular matrix defined by the tuple of
2119 * \a this array, which contains 6 components. The 6 components of tuples are expected to be stored as follow :<br>
2120 * \a tuple[0] = \c matrix_XX <br>
2121 * \a tuple[1] = \c matrix_YY <br>
2122 * \a tuple[2] = \c matrix_ZZ <br>
2123 * \a tuple[3] = \c matrix_XY <br>
2124 * \a tuple[4] = \c matrix_YZ <br>
2125 * \a tuple[5] = \c matrix_XZ <br>
2126 * \return DataArrayDouble * - the new instance of DataArrayDouble containing 3
2127 * components, whose each tuple contains the eigenvalues of the matrix of
2128 * corresponding tuple of \a this array.
2129 * The caller is to delete this result array using decrRef() as it is no more
2131 * \throw If \a this->getNumberOfComponents() != 6.
2133 DataArrayDouble *DataArrayDouble::eigenValues() const
2136 std::size_t nbOfComp=getNumberOfComponents();
2138 throw INTERP_KERNEL::Exception("DataArrayDouble::eigenValues : must be an array with exactly 6 components !");
2139 DataArrayDouble *ret=DataArrayDouble::New();
2140 mcIdType nbOfTuple=getNumberOfTuples();
2141 ret->alloc(nbOfTuple,3);
2142 const double *src=getConstPointer();
2143 double *dest=ret->getPointer();
2144 for(mcIdType i=0;i<nbOfTuple;i++,dest+=3,src+=6)
2145 INTERP_KERNEL::computeEigenValues6(src,dest);
2150 * Computes 3 eigenvectors of every upper triangular matrix defined by the tuple of
2151 * \a this array, which contains 6 components. The 6 components of tuples are expected to be stored as follow :<br>
2152 * \a tuple[0] = \c matrix_XX <br>
2153 * \a tuple[1] = \c matrix_YY <br>
2154 * \a tuple[2] = \c matrix_ZZ <br>
2155 * \a tuple[3] = \c matrix_XY <br>
2156 * \a tuple[4] = \c matrix_YZ <br>
2157 * \a tuple[5] = \c matrix_XZ <br>
2158 * \return DataArrayDouble * - the new instance of DataArrayDouble containing 9
2159 * components, whose each tuple contains 3 eigenvectors of the matrix of
2160 * corresponding tuple of \a this array.
2161 * The caller is to delete this result array using decrRef() as it is no more
2163 * \throw If \a this->getNumberOfComponents() != 6.
2165 DataArrayDouble *DataArrayDouble::eigenVectors() const
2168 std::size_t nbOfComp=getNumberOfComponents();
2170 throw INTERP_KERNEL::Exception("DataArrayDouble::eigenVectors : must be an array with exactly 6 components !");
2171 DataArrayDouble *ret=DataArrayDouble::New();
2172 mcIdType nbOfTuple=getNumberOfTuples();
2173 ret->alloc(nbOfTuple,9);
2174 const double *src=getConstPointer();
2175 double *dest=ret->getPointer();
2176 for(mcIdType i=0;i<nbOfTuple;i++,src+=6)
2179 INTERP_KERNEL::computeEigenValues6(src,tmp);
2180 for(mcIdType j=0;j<3;j++,dest+=3)
2181 INTERP_KERNEL::computeEigenVectorForEigenValue6(src,tmp[j],1e-12,dest);
2187 * Computes the inverse matrix of every matrix defined by the tuple of \a this
2188 * array, which contains either 4, 6 or 9 components. The case of 6 components
2189 * corresponds to that of the upper triangular matrix.
2190 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2191 * same number of components as \a this one, whose each tuple is the inverse
2192 * matrix of the matrix of corresponding tuple of \a this array.
2193 * The caller is to delete this result array using decrRef() as it is no more
2195 * \throw If \a this->getNumberOfComponents() is not in [4,6,9].
2197 DataArrayDouble *DataArrayDouble::inverse() const
2200 std::size_t nbOfComp=getNumberOfComponents();
2201 if(nbOfComp!=6 && nbOfComp!=9 && nbOfComp!=4)
2202 throw INTERP_KERNEL::Exception("DataArrayDouble::inversion : must be an array with 4,6 or 9 components !");
2203 DataArrayDouble *ret=DataArrayDouble::New();
2204 mcIdType nbOfTuple=getNumberOfTuples();
2205 ret->alloc(nbOfTuple,nbOfComp);
2206 const double *src=getConstPointer();
2207 double *dest=ret->getPointer();
2209 for(mcIdType i=0;i<nbOfTuple;i++,dest+=6,src+=6)
2211 double det=src[0]*src[1]*src[2]+2.*src[4]*src[5]*src[3]-src[0]*src[4]*src[4]-src[2]*src[3]*src[3]-src[1]*src[5]*src[5];
2212 dest[0]=(src[1]*src[2]-src[4]*src[4])/det;
2213 dest[1]=(src[0]*src[2]-src[5]*src[5])/det;
2214 dest[2]=(src[0]*src[1]-src[3]*src[3])/det;
2215 dest[3]=(src[5]*src[4]-src[3]*src[2])/det;
2216 dest[4]=(src[5]*src[3]-src[0]*src[4])/det;
2217 dest[5]=(src[3]*src[4]-src[1]*src[5])/det;
2219 else if(nbOfComp==4)
2220 for(mcIdType i=0;i<nbOfTuple;i++,dest+=4,src+=4)
2222 double det=src[0]*src[3]-src[1]*src[2];
2224 dest[1]=-src[1]/det;
2225 dest[2]=-src[2]/det;
2229 for(mcIdType i=0;i<nbOfTuple;i++,dest+=9,src+=9)
2231 double det=src[0]*src[4]*src[8]+src[1]*src[5]*src[6]+src[2]*src[3]*src[7]-src[0]*src[5]*src[7]-src[1]*src[3]*src[8]-src[2]*src[4]*src[6];
2232 dest[0]=(src[4]*src[8]-src[7]*src[5])/det;
2233 dest[1]=(src[7]*src[2]-src[1]*src[8])/det;
2234 dest[2]=(src[1]*src[5]-src[4]*src[2])/det;
2235 dest[3]=(src[6]*src[5]-src[3]*src[8])/det;
2236 dest[4]=(src[0]*src[8]-src[6]*src[2])/det;
2237 dest[5]=(src[2]*src[3]-src[0]*src[5])/det;
2238 dest[6]=(src[3]*src[7]-src[6]*src[4])/det;
2239 dest[7]=(src[6]*src[1]-src[0]*src[7])/det;
2240 dest[8]=(src[0]*src[4]-src[1]*src[3])/det;
2246 * Computes the trace of every matrix defined by the tuple of \a this
2247 * array, which contains either 4, 6 or 9 components. The case of 6 components
2248 * corresponds to that of the upper triangular matrix.
2249 * \return DataArrayDouble * - the new instance of DataArrayDouble containing
2250 * 1 component, whose each tuple is the trace of
2251 * the matrix of corresponding tuple of \a this array.
2252 * The caller is to delete this result array using decrRef() as it is no more
2254 * \throw If \a this->getNumberOfComponents() is not in [4,6,9].
2256 DataArrayDouble *DataArrayDouble::trace() const
2259 std::size_t nbOfComp=getNumberOfComponents();
2260 if(nbOfComp!=6 && nbOfComp!=9 && nbOfComp!=4)
2261 throw INTERP_KERNEL::Exception("DataArrayDouble::trace : must be an array with 4,6 or 9 components !");
2262 DataArrayDouble *ret=DataArrayDouble::New();
2263 mcIdType nbOfTuple=getNumberOfTuples();
2264 ret->alloc(nbOfTuple,1);
2265 const double *src=getConstPointer();
2266 double *dest=ret->getPointer();
2268 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=6)
2269 *dest=src[0]+src[1]+src[2];
2270 else if(nbOfComp==4)
2271 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=4)
2272 *dest=src[0]+src[3];
2274 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=9)
2275 *dest=src[0]+src[4]+src[8];
2280 * Computes the stress deviator tensor of every stress tensor defined by the tuple of
2281 * \a this array, which contains 6 components.
2282 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2283 * same number of components and tuples as \a this array.
2284 * The caller is to delete this result array using decrRef() as it is no more
2286 * \throw If \a this->getNumberOfComponents() != 6.
2288 DataArrayDouble *DataArrayDouble::deviator() const
2291 std::size_t nbOfComp=getNumberOfComponents();
2293 throw INTERP_KERNEL::Exception("DataArrayDouble::deviator : must be an array with exactly 6 components !");
2294 DataArrayDouble *ret=DataArrayDouble::New();
2295 mcIdType nbOfTuple=getNumberOfTuples();
2296 ret->alloc(nbOfTuple,6);
2297 const double *src=getConstPointer();
2298 double *dest=ret->getPointer();
2299 for(mcIdType i=0;i<nbOfTuple;i++,dest+=6,src+=6)
2301 double tr=(src[0]+src[1]+src[2])/3.;
2313 * Computes the magnitude of every vector defined by the tuple of
2315 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2316 * same number of tuples as \a this array and one component.
2317 * The caller is to delete this result array using decrRef() as it is no more
2319 * \throw If \a this is not allocated.
2321 DataArrayDouble *DataArrayDouble::magnitude() const
2324 std::size_t nbOfComp=getNumberOfComponents();
2325 DataArrayDouble *ret=DataArrayDouble::New();
2326 mcIdType nbOfTuple=getNumberOfTuples();
2327 ret->alloc(nbOfTuple,1);
2328 const double *src=getConstPointer();
2329 double *dest=ret->getPointer();
2330 for(mcIdType i=0;i<nbOfTuple;i++,dest++)
2333 for(std::size_t j=0;j<nbOfComp;j++,src++)
2340 DataArrayDouble *DataArrayDouble::operatePerTuple(std::function<double(const double *bg, const double *endd)> func) const
2343 std::size_t nbOfComp(getNumberOfComponents());
2344 MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
2345 mcIdType nbOfTuple(getNumberOfTuples());
2346 ret->alloc(nbOfTuple,1);
2347 const double *src=getConstPointer();
2348 double *dest=ret->getPointer();
2349 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=nbOfComp)
2350 *dest=func(src,src+nbOfComp);
2355 * Computes the maximal value within every tuple of \a this array.
2356 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2357 * same number of tuples as \a this array and one component.
2358 * The caller is to delete this result array using decrRef() as it is no more
2360 * \throw If \a this is not allocated.
2361 * \sa DataArrayDouble::maxPerTupleWithCompoId, DataArrayDouble::minPerTuple
2363 DataArrayDouble *DataArrayDouble::maxPerTuple() const
2365 return this->operatePerTuple([](const double *bg, const double *endd) { return *std::max_element(bg,endd); });
2369 * Computes the minimal value within every tuple of \a this array.
2370 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2371 * same number of tuples as \a this array and one component.
2372 * The caller is to delete this result array using decrRef() as it is no more
2374 * \throw If \a this is not allocated.
2375 * \sa DataArrayDouble::maxPerTuple
2377 DataArrayDouble *DataArrayDouble::minPerTuple() const
2379 return this->operatePerTuple([](const double *bg, const double *endd) { return *std::min_element(bg,endd); });
2383 * Computes the maximal value within every tuple of \a this array and it returns the first component
2384 * id for each tuple that corresponds to the maximal value within the tuple.
2386 * \param [out] compoIdOfMaxPerTuple - the new new instance of DataArrayInt containing the
2387 * same number of tuples and only one component.
2388 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2389 * same number of tuples as \a this array and one component.
2390 * The caller is to delete this result array using decrRef() as it is no more
2392 * \throw If \a this is not allocated.
2393 * \sa DataArrayDouble::maxPerTuple
2395 DataArrayDouble *DataArrayDouble::maxPerTupleWithCompoId(DataArrayIdType* &compoIdOfMaxPerTuple) const
2398 std::size_t nbOfComp(getNumberOfComponents());
2399 MCAuto<DataArrayDouble> ret0=DataArrayDouble::New();
2400 MCAuto<DataArrayIdType> ret1=DataArrayIdType::New();
2401 mcIdType nbOfTuple=getNumberOfTuples();
2402 ret0->alloc(nbOfTuple,1); ret1->alloc(nbOfTuple,1);
2403 const double *src=getConstPointer();
2404 double *dest=ret0->getPointer(); mcIdType *dest1=ret1->getPointer();
2405 for(mcIdType i=0;i<nbOfTuple;i++,dest++,dest1++,src+=nbOfComp)
2407 const double *loc=std::max_element(src,src+nbOfComp);
2409 *dest1=ToIdType(std::distance(src,loc));
2411 compoIdOfMaxPerTuple=ret1.retn();
2416 * This method returns a newly allocated DataArrayDouble instance having one component and \c this->getNumberOfTuples() * \c this->getNumberOfTuples() tuples.
2417 * \n This returned array contains the euclidian distance for each tuple in \a this.
2418 * \n So the returned array can be seen as a dense symmetrical matrix whose diagonal elements are equal to 0.
2419 * \n The returned array has only one component (and **not** \c this->getNumberOfTuples() components to avoid the useless memory consumption due to components info in returned DataArrayDouble)
2421 * \warning use this method with care because it can leads to big amount of consumed memory !
2423 * \return A newly allocated (huge) MEDCoupling::DataArrayDouble instance that the caller should deal with.
2425 * \throw If \a this is not allocated.
2427 * \sa DataArrayDouble::buildEuclidianDistanceDenseMatrixWith
2429 DataArrayDouble *DataArrayDouble::buildEuclidianDistanceDenseMatrix() const
2432 std::size_t nbOfComp(getNumberOfComponents());
2433 mcIdType nbOfTuples(getNumberOfTuples());
2434 const double *inData=getConstPointer();
2435 MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
2436 ret->alloc(nbOfTuples*nbOfTuples,1);
2437 double *outData=ret->getPointer();
2438 for(mcIdType i=0;i<nbOfTuples;i++)
2440 outData[i*nbOfTuples+i]=0.;
2441 for(mcIdType j=i+1;j<nbOfTuples;j++)
2444 for(std::size_t k=0;k<nbOfComp;k++)
2445 { double delta=inData[i*nbOfComp+k]-inData[j*nbOfComp+k]; dist+=delta*delta; }
2447 outData[i*nbOfTuples+j]=dist;
2448 outData[j*nbOfTuples+i]=dist;
2455 * This method returns a newly allocated DataArrayDouble instance having one component and \c this->getNumberOfTuples() * \c other->getNumberOfTuples() tuples.
2456 * \n This returned array contains the euclidian distance for each tuple in \a other with each tuple in \a this.
2457 * \n So the returned array can be seen as a dense rectangular matrix with \c other->getNumberOfTuples() rows and \c this->getNumberOfTuples() columns.
2458 * \n Output rectangular matrix is sorted along rows.
2459 * \n The returned array has only one component (and **not** \c this->getNumberOfTuples() components to avoid the useless memory consumption due to components info in returned DataArrayDouble)
2461 * \warning use this method with care because it can leads to big amount of consumed memory !
2463 * \param [in] other DataArrayDouble instance having same number of components than \a this.
2464 * \return A newly allocated (huge) MEDCoupling::DataArrayDouble instance that the caller should deal with.
2466 * \throw If \a this is not allocated, or if \a other is null or if \a other is not allocated, or if number of components of \a other and \a this differs.
2468 * \sa DataArrayDouble::buildEuclidianDistanceDenseMatrix
2470 DataArrayDouble *DataArrayDouble::buildEuclidianDistanceDenseMatrixWith(const DataArrayDouble *other) const
2473 throw INTERP_KERNEL::Exception("DataArrayDouble::buildEuclidianDistanceDenseMatrixWith : input parameter is null !");
2475 other->checkAllocated();
2476 std::size_t nbOfComp(getNumberOfComponents());
2477 std::size_t otherNbOfComp(other->getNumberOfComponents());
2478 if(nbOfComp!=otherNbOfComp)
2480 std::ostringstream oss; oss << "DataArrayDouble::buildEuclidianDistanceDenseMatrixWith : this nb of compo=" << nbOfComp << " and other nb of compo=" << otherNbOfComp << ". It should match !";
2481 throw INTERP_KERNEL::Exception(oss.str().c_str());
2483 mcIdType nbOfTuples(getNumberOfTuples());
2484 mcIdType otherNbOfTuples(other->getNumberOfTuples());
2485 const double *inData=getConstPointer();
2486 const double *inDataOther=other->getConstPointer();
2487 MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
2488 ret->alloc(otherNbOfTuples*nbOfTuples,1);
2489 double *outData=ret->getPointer();
2490 for(mcIdType i=0;i<otherNbOfTuples;i++,inDataOther+=nbOfComp)
2492 for(mcIdType j=0;j<nbOfTuples;j++)
2495 for(std::size_t k=0;k<nbOfComp;k++)
2496 { double delta=inDataOther[k]-inData[j*nbOfComp+k]; dist+=delta*delta; }
2498 outData[i*nbOfTuples+j]=dist;
2505 * This method expects that \a this stores 3 tuples containing 2 components each.
2506 * Each of this tuples represent a point into 2D space.
2507 * This method tries to find an arc of circle starting from first point (tuple) to 2nd and middle point (tuple) along 3nd and last point (tuple).
2508 * If such arc of circle exists, the corresponding center, radius of circle is returned. And additionnaly the length of arc expressed with an \a ang output variable in ]0,2*pi[.
2510 * \throw If \a this is not allocated.
2511 * \throw If \a this has not 3 tuples of 2 components
2512 * \throw If tuples/points in \a this are aligned
2514 void DataArrayDouble::asArcOfCircle(double center[2], double& radius, double& ang) const
2517 INTERP_KERNEL::QuadraticPlanarPrecision arcPrec(1e-14);
2518 if(getNumberOfTuples()!=3 && getNumberOfComponents()!=2)
2519 throw INTERP_KERNEL::Exception("DataArrayDouble::asArcCircle : this method expects");
2520 const double *pt(begin());
2521 MCAuto<INTERP_KERNEL::Node> n0(new INTERP_KERNEL::Node(pt[0],pt[1])),n1(new INTERP_KERNEL::Node(pt[2],pt[3])),n2(new INTERP_KERNEL::Node(pt[4],pt[5]));
2523 INTERP_KERNEL::AutoCppPtr<INTERP_KERNEL::EdgeLin> e1(new INTERP_KERNEL::EdgeLin(n0,n2)),e2(new INTERP_KERNEL::EdgeLin(n2,n1));
2524 INTERP_KERNEL::SegSegIntersector inters(*e1,*e2);
2525 bool colinearity(inters.areColinears());
2527 throw INTERP_KERNEL::Exception("DataArrayDouble::asArcOfCircle : 3 points in this have been detected as colinear !");
2529 INTERP_KERNEL::AutoCppPtr<INTERP_KERNEL::EdgeArcCircle> ret(new INTERP_KERNEL::EdgeArcCircle(n0,n2,n1));
2530 const double *c(ret->getCenter());
2531 center[0]=c[0]; center[1]=c[1];
2532 radius=ret->getRadius();
2533 ang=ret->getAngle();
2537 * Sorts value within every tuple of \a this array.
2538 * \param [in] asc - if \a true, the values are sorted in ascending order, else,
2539 * in descending order.
2540 * \throw If \a this is not allocated.
2542 void DataArrayDouble::sortPerTuple(bool asc)
2545 double *pt=getPointer();
2546 mcIdType nbOfTuple(getNumberOfTuples());
2547 std::size_t nbOfComp(getNumberOfComponents());
2549 for(mcIdType i=0;i<nbOfTuple;i++,pt+=nbOfComp)
2550 std::sort(pt,pt+nbOfComp);
2552 for(mcIdType i=0;i<nbOfTuple;i++,pt+=nbOfComp)
2553 std::sort(pt,pt+nbOfComp,std::greater<double>());
2558 * Modify all elements of \a this array, so that
2559 * an element _x_ becomes \f$ numerator / x \f$.
2560 * \warning If an exception is thrown because of presence of 0.0 element in \a this
2561 * array, all elements processed before detection of the zero element remain
2563 * \param [in] numerator - the numerator used to modify array elements.
2564 * \throw If \a this is not allocated.
2565 * \throw If there is an element equal to 0.0 in \a this array.
2567 void DataArrayDouble::applyInv(double numerator)
2570 double *ptr=getPointer();
2571 std::size_t nbOfElems=getNbOfElems();
2572 for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2574 if(std::abs(*ptr)>std::numeric_limits<double>::min())
2576 *ptr=numerator/(*ptr);
2580 std::ostringstream oss; oss << "DataArrayDouble::applyInv : presence of null value in tuple #" << i/getNumberOfComponents() << " component #" << i%getNumberOfComponents();
2582 throw INTERP_KERNEL::Exception(oss.str().c_str());
2589 * Modify all elements of \a this array, so that
2590 * an element _x_ becomes <em> val ^ x </em>. Contrary to DataArrayInt::applyPow
2591 * all values in \a this have to be >= 0 if val is \b not integer.
2592 * \param [in] val - the value used to apply pow on all array elements.
2593 * \throw If \a this is not allocated.
2594 * \warning If an exception is thrown because of presence of 0 element in \a this
2595 * array and \a val is \b not integer, all elements processed before detection of the zero element remain
2598 void DataArrayDouble::applyPow(double val)
2601 double *ptr=getPointer();
2602 std::size_t nbOfElems=getNbOfElems();
2604 bool isInt=((double)val2)==val;
2607 for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2613 std::ostringstream oss; oss << "DataArrayDouble::applyPow (double) : At elem # " << i << " value is " << *ptr << " ! must be >=0. !";
2614 throw INTERP_KERNEL::Exception(oss.str().c_str());
2620 for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2621 *ptr=pow(*ptr,val2);
2627 * Modify all elements of \a this array, so that
2628 * an element _x_ becomes \f$ val ^ x \f$.
2629 * \param [in] val - the value used to apply pow on all array elements.
2630 * \throw If \a this is not allocated.
2631 * \throw If \a val < 0.
2632 * \warning If an exception is thrown because of presence of 0 element in \a this
2633 * array, all elements processed before detection of the zero element remain
2636 void DataArrayDouble::applyRPow(double val)
2640 throw INTERP_KERNEL::Exception("DataArrayDouble::applyRPow : the input value has to be >= 0 !");
2641 double *ptr=getPointer();
2642 std::size_t nbOfElems=getNbOfElems();
2643 for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2649 * Returns a new DataArrayDouble created from \a this one by applying \a
2650 * FunctionToEvaluate to every tuple of \a this array. Textual data is not copied.
2651 * For more info see \ref MEDCouplingArrayApplyFunc
2652 * \param [in] nbOfComp - number of components in the result array.
2653 * \param [in] func - the \a FunctionToEvaluate declared as
2654 * \c bool (*\a func)(\c const \c double *\a pos, \c double *\a res),
2655 * where \a pos points to the first component of a tuple of \a this array
2656 * and \a res points to the first component of a tuple of the result array.
2657 * Note that length (number of components) of \a pos can differ from
2659 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2660 * same number of tuples as \a this array.
2661 * The caller is to delete this result array using decrRef() as it is no more
2663 * \throw If \a this is not allocated.
2664 * \throw If \a func returns \a false.
2666 DataArrayDouble *DataArrayDouble::applyFunc(std::size_t nbOfComp, FunctionToEvaluate func) const
2669 DataArrayDouble *newArr=DataArrayDouble::New();
2670 mcIdType nbOfTuples(getNumberOfTuples());
2671 std::size_t oldNbOfComp(getNumberOfComponents());
2672 newArr->alloc(nbOfTuples,nbOfComp);
2673 const double *ptr=getConstPointer();
2674 double *ptrToFill=newArr->getPointer();
2675 for(mcIdType i=0;i<nbOfTuples;i++)
2677 if(!func(ptr+i*oldNbOfComp,ptrToFill+i*nbOfComp))
2679 std::ostringstream oss; oss << "For tuple # " << i << " with value (";
2680 std::copy(ptr+oldNbOfComp*i,ptr+oldNbOfComp*(i+1),std::ostream_iterator<double>(oss,", "));
2681 oss << ") : Evaluation of function failed !";
2683 throw INTERP_KERNEL::Exception(oss.str().c_str());
2690 * Returns a new DataArrayDouble created from \a this one by applying a function to every
2691 * tuple of \a this array. Textual data is not copied.
2692 * For more info see \ref MEDCouplingArrayApplyFunc1.
2693 * \param [in] nbOfComp - number of components in the result array.
2694 * \param [in] func - the expression defining how to transform a tuple of \a this array.
2695 * Supported expressions are described \ref MEDCouplingArrayApplyFuncExpr "here".
2696 * \param [in] isSafe - By default true. If true invalid operation (division by 0. acos of value > 1. ...) leads to a throw of an exception.
2697 * If false the computation is carried on without any notification. When false the evaluation is a little faster.
2698 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2699 * same number of tuples as \a this array and \a nbOfComp components.
2700 * The caller is to delete this result array using decrRef() as it is no more
2702 * \throw If \a this is not allocated.
2703 * \throw If computing \a func fails.
2705 DataArrayDouble *DataArrayDouble::applyFunc(std::size_t nbOfComp, const std::string& func, bool isSafe) const
2707 INTERP_KERNEL::ExprParser expr(func);
2709 std::set<std::string> vars;
2710 expr.getTrueSetOfVars(vars);
2711 std::vector<std::string> varsV(vars.begin(),vars.end());
2712 return applyFuncNamedCompo(nbOfComp,varsV,func,isSafe);
2716 * Returns a new DataArrayDouble created from \a this one by applying a function to every
2717 * tuple of \a this array. Textual data is not copied. This method works by tuples (whatever its size).
2718 * If \a this is a one component array, call applyFuncOnThis instead that performs the same work faster.
2720 * For more info see \ref MEDCouplingArrayApplyFunc0.
2721 * \param [in] func - the expression defining how to transform a tuple of \a this array.
2722 * Supported expressions are described \ref MEDCouplingArrayApplyFuncExpr "here".
2723 * \param [in] isSafe - By default true. If true invalid operation (division by 0. acos of value > 1. ...) leads to a throw of an exception.
2724 * If false the computation is carried on without any notification. When false the evaluation is a little faster.
2725 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2726 * same number of tuples and components as \a this array.
2727 * The caller is to delete this result array using decrRef() as it is no more
2729 * \sa applyFuncOnThis
2730 * \throw If \a this is not allocated.
2731 * \throw If computing \a func fails.
2733 DataArrayDouble *DataArrayDouble::applyFunc(const std::string& func, bool isSafe) const
2735 std::size_t nbOfComp(getNumberOfComponents());
2737 throw INTERP_KERNEL::Exception("DataArrayDouble::applyFunc : output number of component must be > 0 !");
2739 mcIdType nbOfTuples(getNumberOfTuples());
2740 MCAuto<DataArrayDouble> newArr(DataArrayDouble::New());
2741 newArr->alloc(nbOfTuples,nbOfComp);
2742 INTERP_KERNEL::ExprParser expr(func);
2744 std::set<std::string> vars;
2745 expr.getTrueSetOfVars(vars);
2748 std::ostringstream oss; oss << "DataArrayDouble::applyFunc : this method works only with at most one var func expression ! If you need to map comps on variables please use applyFuncCompo or applyFuncNamedCompo instead ! Vars in expr are : ";
2749 std::copy(vars.begin(),vars.end(),std::ostream_iterator<std::string>(oss," "));
2750 throw INTERP_KERNEL::Exception(oss.str().c_str());
2754 expr.prepareFastEvaluator();
2755 newArr->rearrange(1);
2756 newArr->fillWithValue(expr.evaluateDouble());
2757 newArr->rearrange(nbOfComp);
2758 return newArr.retn();
2760 std::vector<std::string> vars2(vars.begin(),vars.end());
2761 double buff,*ptrToFill(newArr->getPointer());
2762 const double *ptr(begin());
2763 std::vector<double> stck;
2764 expr.prepareExprEvaluationDouble(vars2,1,1,0,&buff,&buff+1);
2765 expr.prepareFastEvaluator();
2768 for(mcIdType i=0;i<nbOfTuples;i++)
2770 for(std::size_t iComp=0;iComp<nbOfComp;iComp++,ptr++,ptrToFill++)
2773 expr.evaluateDoubleInternal(stck);
2774 *ptrToFill=stck.back();
2781 for(mcIdType i=0;i<nbOfTuples;i++)
2783 for(std::size_t iComp=0;iComp<nbOfComp;iComp++,ptr++,ptrToFill++)
2788 expr.evaluateDoubleInternalSafe(stck);
2790 catch(INTERP_KERNEL::Exception& e)
2792 std::ostringstream oss; oss << "For tuple # " << i << " component # " << iComp << " with value (";
2794 oss << ") : Evaluation of function failed !" << e.what();
2795 throw INTERP_KERNEL::Exception(oss.str().c_str());
2797 *ptrToFill=stck.back();
2802 return newArr.retn();
2806 * This method is a non const method that modify the array in \a this.
2807 * This method only works on one component array. It means that function \a func must
2808 * contain at most one variable.
2809 * This method is a specialization of applyFunc method with one parameter on one component array.
2811 * \param [in] func - the expression defining how to transform a tuple of \a this array.
2812 * Supported expressions are described \ref MEDCouplingArrayApplyFuncExpr "here".
2813 * \param [in] isSafe - By default true. If true invalid operation (division by 0. acos of value > 1. ...) leads to a throw of an exception.
2814 * If false the computation is carried on without any notification. When false the evaluation is a little faster.
2818 void DataArrayDouble::applyFuncOnThis(const std::string& func, bool isSafe)
2820 std::size_t nbOfComp(getNumberOfComponents());
2822 throw INTERP_KERNEL::Exception("DataArrayDouble::applyFuncOnThis : output number of component must be > 0 !");
2824 mcIdType nbOfTuples(getNumberOfTuples());
2825 INTERP_KERNEL::ExprParser expr(func);
2827 std::set<std::string> vars;
2828 expr.getTrueSetOfVars(vars);
2831 std::ostringstream oss; oss << "DataArrayDouble::applyFuncOnThis : this method works only with at most one var func expression ! If you need to map comps on variables please use applyFuncCompo or applyFuncNamedCompo instead ! Vars in expr are : ";
2832 std::copy(vars.begin(),vars.end(),std::ostream_iterator<std::string>(oss," "));
2833 throw INTERP_KERNEL::Exception(oss.str().c_str());
2837 expr.prepareFastEvaluator();
2838 std::vector<std::string> compInfo(getInfoOnComponents());
2840 fillWithValue(expr.evaluateDouble());
2841 rearrange(nbOfComp);
2842 setInfoOnComponents(compInfo);
2845 std::vector<std::string> vars2(vars.begin(),vars.end());
2846 double buff,*ptrToFill(getPointer());
2847 const double *ptr(begin());
2848 std::vector<double> stck;
2849 expr.prepareExprEvaluationDouble(vars2,1,1,0,&buff,&buff+1);
2850 expr.prepareFastEvaluator();
2853 for(mcIdType i=0;i<nbOfTuples;i++)
2855 for(std::size_t iComp=0;iComp<nbOfComp;iComp++,ptr++,ptrToFill++)
2858 expr.evaluateDoubleInternal(stck);
2859 *ptrToFill=stck.back();
2866 for(mcIdType i=0;i<nbOfTuples;i++)
2868 for(std::size_t iComp=0;iComp<nbOfComp;iComp++,ptr++,ptrToFill++)
2873 expr.evaluateDoubleInternalSafe(stck);
2875 catch(INTERP_KERNEL::Exception& e)
2877 std::ostringstream oss; oss << "For tuple # " << i << " component # " << iComp << " with value (";
2879 oss << ") : Evaluation of function failed !" << e.what();
2880 throw INTERP_KERNEL::Exception(oss.str().c_str());
2882 *ptrToFill=stck.back();
2890 * Returns a new DataArrayDouble created from \a this one by applying a function to every
2891 * tuple of \a this array. Textual data is not copied.
2892 * For more info see \ref MEDCouplingArrayApplyFunc2.
2893 * \param [in] nbOfComp - number of components in the result array.
2894 * \param [in] func - the expression defining how to transform a tuple of \a this array.
2895 * Supported expressions are described \ref MEDCouplingArrayApplyFuncExpr "here".
2896 * \param [in] isSafe - By default true. If true invalid operation (division by 0. acos of value > 1. ...) leads to a throw of an exception.
2897 * If false the computation is carried on without any notification. When false the evaluation is a little faster.
2898 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2899 * same number of tuples as \a this array.
2900 * The caller is to delete this result array using decrRef() as it is no more
2902 * \throw If \a this is not allocated.
2903 * \throw If \a func contains vars that are not in \a this->getInfoOnComponent().
2904 * \throw If computing \a func fails.
2906 DataArrayDouble *DataArrayDouble::applyFuncCompo(std::size_t nbOfComp, const std::string& func, bool isSafe) const
2908 return applyFuncNamedCompo(nbOfComp,getVarsOnComponent(),func,isSafe);
2912 * Returns a new DataArrayDouble created from \a this one by applying a function to every
2913 * tuple of \a this array. Textual data is not copied.
2914 * For more info see \ref MEDCouplingArrayApplyFunc3.
2915 * \param [in] nbOfComp - number of components in the result array.
2916 * \param [in] varsOrder - sequence of vars defining their order.
2917 * \param [in] func - the expression defining how to transform a tuple of \a this array.
2918 * Supported expressions are described \ref MEDCouplingArrayApplyFuncExpr "here".
2919 * \param [in] isSafe - By default true. If true invalid operation (division by 0. acos of value > 1. ...) leads to a throw of an exception.
2920 * If false the computation is carried on without any notification. When false the evaluation is a little faster.
2921 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2922 * same number of tuples as \a this array.
2923 * The caller is to delete this result array using decrRef() as it is no more
2925 * \throw If \a this is not allocated.
2926 * \throw If \a func contains vars not in \a varsOrder.
2927 * \throw If computing \a func fails.
2929 DataArrayDouble *DataArrayDouble::applyFuncNamedCompo(std::size_t nbOfComp, const std::vector<std::string>& varsOrder, const std::string& func, bool isSafe) const
2932 throw INTERP_KERNEL::Exception("DataArrayDouble::applyFuncNamedCompo : output number of component must be > 0 !");
2933 std::vector<std::string> varsOrder2(varsOrder);
2934 std::size_t oldNbOfComp(getNumberOfComponents());
2935 for(std::size_t i=varsOrder.size();i<oldNbOfComp;i++)
2936 varsOrder2.push_back(std::string());
2938 mcIdType nbOfTuples(getNumberOfTuples());
2939 INTERP_KERNEL::ExprParser expr(func);
2941 std::set<std::string> vars;
2942 expr.getTrueSetOfVars(vars);
2943 if(vars.size()>oldNbOfComp)
2945 std::ostringstream oss; oss << "The field has " << oldNbOfComp << " components and there are ";
2946 oss << vars.size() << " variables : ";
2947 std::copy(vars.begin(),vars.end(),std::ostream_iterator<std::string>(oss," "));
2948 throw INTERP_KERNEL::Exception(oss.str().c_str());
2950 MCAuto<DataArrayDouble> newArr(DataArrayDouble::New());
2951 newArr->alloc(nbOfTuples,nbOfComp);
2952 INTERP_KERNEL::AutoPtr<double> buff(new double[oldNbOfComp]);
2953 double *buffPtr(buff),*ptrToFill;
2954 std::vector<double> stck;
2955 for(std::size_t iComp=0;iComp<nbOfComp;iComp++)
2957 expr.prepareExprEvaluationDouble(varsOrder2,(int)oldNbOfComp,(int)nbOfComp,(int)iComp,buffPtr,buffPtr+oldNbOfComp);
2958 expr.prepareFastEvaluator();
2959 const double *ptr(getConstPointer());
2960 ptrToFill=newArr->getPointer()+iComp;
2963 for(mcIdType i=0;i<nbOfTuples;i++,ptrToFill+=nbOfComp,ptr+=oldNbOfComp)
2965 std::copy(ptr,ptr+oldNbOfComp,buffPtr);
2966 expr.evaluateDoubleInternal(stck);
2967 *ptrToFill=stck.back();
2973 for(mcIdType i=0;i<nbOfTuples;i++,ptrToFill+=nbOfComp,ptr+=oldNbOfComp)
2975 std::copy(ptr,ptr+oldNbOfComp,buffPtr);
2978 expr.evaluateDoubleInternalSafe(stck);
2979 *ptrToFill=stck.back();
2982 catch(INTERP_KERNEL::Exception& e)
2984 std::ostringstream oss; oss << "For tuple # " << i << " with value (";
2985 std::copy(ptr+oldNbOfComp*i,ptr+oldNbOfComp*(i+1),std::ostream_iterator<double>(oss,", "));
2986 oss << ") : Evaluation of function failed !" << e.what();
2987 throw INTERP_KERNEL::Exception(oss.str().c_str());
2992 return newArr.retn();
2995 void DataArrayDouble::applyFuncFast32(const std::string& func)
2998 INTERP_KERNEL::ExprParser expr(func);
3000 char *funcStr=expr.compileX86();
3002 *((void **)&funcPtr)=funcStr;//he he...
3004 double *ptr=getPointer();
3005 std::size_t nbOfComp=getNumberOfComponents();
3006 mcIdType nbOfTuples=getNumberOfTuples();
3007 std::size_t nbOfElems=nbOfTuples*nbOfComp;
3008 for(std::size_t i=0;i<nbOfElems;i++,ptr++)
3013 void DataArrayDouble::applyFuncFast64(const std::string& func)
3016 INTERP_KERNEL::ExprParser expr(func);
3018 char *funcStr=expr.compileX86_64();
3020 *((void **)&funcPtr)=funcStr;//he he...
3022 double *ptr=getPointer();
3023 std::size_t nbOfComp=getNumberOfComponents();
3024 mcIdType nbOfTuples=getNumberOfTuples();
3025 std::size_t nbOfElems=nbOfTuples*nbOfComp;
3026 for(std::size_t i=0;i<nbOfElems;i++,ptr++)
3032 * \return a new object that is the result of the symmetry along 3D plane defined by its normal vector \a normalVector and a point \a point.
3034 MCAuto<DataArrayDouble> DataArrayDouble::symmetry3DPlane(const double point[3], const double normalVector[3]) const
3037 if(getNumberOfComponents()!=3)
3038 throw INTERP_KERNEL::Exception("DataArrayDouble::symmetry3DPlane : this is excepted to have 3 components !");
3039 mcIdType nbTuples(getNumberOfTuples());
3040 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
3041 ret->alloc(nbTuples,3);
3042 Symmetry3DPlane(point,normalVector,nbTuples,begin(),ret->getPointer());
3046 DataArrayDoubleIterator *DataArrayDouble::iterator()
3048 return new DataArrayDoubleIterator(this);
3052 * Returns a new DataArrayInt containing indices of tuples of \a this one-dimensional
3053 * array whose values are within a given range. Textual data is not copied.
3054 * \param [in] vmin - a lowest acceptable value (included).
3055 * \param [in] vmax - a greatest acceptable value (included).
3056 * \return DataArrayInt * - the new instance of DataArrayInt.
3057 * The caller is to delete this result array using decrRef() as it is no more
3059 * \throw If \a this->getNumberOfComponents() != 1.
3061 * \sa DataArrayDouble::findIdsNotInRange
3063 * \if ENABLE_EXAMPLES
3064 * \ref cpp_mcdataarraydouble_getidsinrange "Here is a C++ example".<br>
3065 * \ref py_mcdataarraydouble_getidsinrange "Here is a Python example".
3068 DataArrayIdType *DataArrayDouble::findIdsInRange(double vmin, double vmax) const
3071 if(getNumberOfComponents()!=1)
3072 throw INTERP_KERNEL::Exception("DataArrayDouble::findIdsInRange : this must have exactly one component !");
3073 const double *cptr(begin());
3074 MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(0,1);
3075 mcIdType nbOfTuples(getNumberOfTuples());
3076 for(mcIdType i=0;i<nbOfTuples;i++,cptr++)
3077 if(*cptr>=vmin && *cptr<=vmax)
3078 ret->pushBackSilent(i);
3083 * Returns a new DataArrayInt containing indices of tuples of \a this one-dimensional
3084 * array whose values are not within a given range. Textual data is not copied.
3085 * \param [in] vmin - a lowest not acceptable value (excluded).
3086 * \param [in] vmax - a greatest not acceptable value (excluded).
3087 * \return DataArrayInt * - the new instance of DataArrayInt.
3088 * The caller is to delete this result array using decrRef() as it is no more
3090 * \throw If \a this->getNumberOfComponents() != 1.
3092 * \sa DataArrayDouble::findIdsInRange
3094 DataArrayIdType *DataArrayDouble::findIdsNotInRange(double vmin, double vmax) const
3097 if(getNumberOfComponents()!=1)
3098 throw INTERP_KERNEL::Exception("DataArrayDouble::findIdsNotInRange : this must have exactly one component !");
3099 const double *cptr(begin());
3100 MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(0,1);
3101 mcIdType nbOfTuples(getNumberOfTuples());
3102 for(mcIdType i=0;i<nbOfTuples;i++,cptr++)
3103 if(*cptr<vmin || *cptr>vmax)
3104 ret->pushBackSilent(i);
3109 * Returns a new DataArrayDouble by concatenating two given arrays, so that (1) the number
3110 * of tuples in the result array is a sum of the number of tuples of given arrays and (2)
3111 * the number of component in the result array is same as that of each of given arrays.
3112 * Info on components is copied from the first of the given arrays. Number of components
3113 * in the given arrays must be the same.
3114 * \param [in] a1 - an array to include in the result array.
3115 * \param [in] a2 - another array to include in the result array.
3116 * \return DataArrayDouble * - the new instance of DataArrayDouble.
3117 * The caller is to delete this result array using decrRef() as it is no more
3119 * \throw If both \a a1 and \a a2 are NULL.
3120 * \throw If \a a1->getNumberOfComponents() != \a a2->getNumberOfComponents().
3122 DataArrayDouble *DataArrayDouble::Aggregate(const DataArrayDouble *a1, const DataArrayDouble *a2)
3124 std::vector<const DataArrayDouble *> tmp(2);
3125 tmp[0]=a1; tmp[1]=a2;
3126 return Aggregate(tmp);
3130 * Returns a new DataArrayDouble by concatenating all given arrays, so that (1) the number
3131 * of tuples in the result array is a sum of the number of tuples of given arrays and (2)
3132 * the number of component in the result array is same as that of each of given arrays.
3133 * Info on components is copied from the first of the given arrays. Number of components
3134 * in the given arrays must be the same.
3135 * If the number of non null of elements in \a arr is equal to one the returned object is a copy of it
3136 * not the object itself.
3137 * \param [in] arr - a sequence of arrays to include in the result array.
3138 * \return DataArrayDouble * - the new instance of DataArrayDouble.
3139 * The caller is to delete this result array using decrRef() as it is no more
3141 * \throw If all arrays within \a arr are NULL.
3142 * \throw If getNumberOfComponents() of arrays within \a arr.
3144 DataArrayDouble *DataArrayDouble::Aggregate(const std::vector<const DataArrayDouble *>& arr)
3146 std::vector<const DataArrayDouble *> a;
3147 for(std::vector<const DataArrayDouble *>::const_iterator it4=arr.begin();it4!=arr.end();it4++)
3151 throw INTERP_KERNEL::Exception("DataArrayDouble::Aggregate : input list must contain at least one NON EMPTY DataArrayDouble !");
3152 std::vector<const DataArrayDouble *>::const_iterator it=a.begin();
3153 std::size_t nbOfComp((*it)->getNumberOfComponents());
3154 mcIdType nbt=(*it++)->getNumberOfTuples();
3155 for(mcIdType i=1;it!=a.end();it++,i++)
3157 if((*it)->getNumberOfComponents()!=nbOfComp)
3158 throw INTERP_KERNEL::Exception("DataArrayDouble::Aggregate : Nb of components mismatch for array aggregation !");
3159 nbt+=(*it)->getNumberOfTuples();
3161 MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
3162 ret->alloc(nbt,nbOfComp);
3163 double *pt=ret->getPointer();
3164 for(it=a.begin();it!=a.end();it++)
3165 pt=std::copy((*it)->getConstPointer(),(*it)->getConstPointer()+(*it)->getNbOfElems(),pt);
3166 ret->copyStringInfoFrom(*(a[0]));
3171 * Returns a new DataArrayDouble containing a dot product of two given arrays, so that
3172 * the i-th tuple of the result array is a sum of products of j-th components of i-th
3173 * tuples of given arrays (\f$ a_i = \sum_{j=1}^n a1_j * a2_j \f$).
3174 * Info on components and name is copied from the first of the given arrays.
3175 * Number of tuples and components in the given arrays must be the same.
3176 * \param [in] a1 - a given array.
3177 * \param [in] a2 - another given array.
3178 * \return DataArrayDouble * - the new instance of DataArrayDouble.
3179 * The caller is to delete this result array using decrRef() as it is no more
3181 * \throw If either \a a1 or \a a2 is NULL.
3182 * \throw If any given array is not allocated.
3183 * \throw If \a a1->getNumberOfTuples() != \a a2->getNumberOfTuples()
3184 * \throw If \a a1->getNumberOfComponents() != \a a2->getNumberOfComponents()
3186 DataArrayDouble *DataArrayDouble::Dot(const DataArrayDouble *a1, const DataArrayDouble *a2)
3189 throw INTERP_KERNEL::Exception("DataArrayDouble::Dot : input DataArrayDouble instance is NULL !");
3190 a1->checkAllocated();
3191 a2->checkAllocated();
3192 std::size_t nbOfComp(a1->getNumberOfComponents());
3193 if(nbOfComp!=a2->getNumberOfComponents())
3194 throw INTERP_KERNEL::Exception("Nb of components mismatch for array Dot !");
3195 mcIdType nbOfTuple(a1->getNumberOfTuples());
3196 if(nbOfTuple!=a2->getNumberOfTuples())
3197 throw INTERP_KERNEL::Exception("Nb of tuples mismatch for array Dot !");
3198 DataArrayDouble *ret=DataArrayDouble::New();
3199 ret->alloc(nbOfTuple,1);
3200 double *retPtr=ret->getPointer();
3201 const double *a1Ptr=a1->begin(),*a2Ptr(a2->begin());
3202 for(mcIdType i=0;i<nbOfTuple;i++)
3205 for(std::size_t j=0;j<nbOfComp;j++)
3206 sum+=a1Ptr[i*nbOfComp+j]*a2Ptr[i*nbOfComp+j];
3209 ret->setInfoOnComponent(0,a1->getInfoOnComponent(0));
3210 ret->setName(a1->getName());
3215 * Returns a new DataArrayDouble containing a cross product of two given arrays, so that
3216 * the i-th tuple of the result array contains 3 components of a vector which is a cross
3217 * product of two vectors defined by the i-th tuples of given arrays.
3218 * Info on components is copied from the first of the given arrays.
3219 * Number of tuples in the given arrays must be the same.
3220 * Number of components in the given arrays must be 3.
3221 * \param [in] a1 - a given array.
3222 * \param [in] a2 - another given array.
3223 * \return DataArrayDouble * - the new instance of DataArrayDouble.
3224 * The caller is to delete this result array using decrRef() as it is no more
3226 * \throw If either \a a1 or \a a2 is NULL.
3227 * \throw If \a a1->getNumberOfTuples() != \a a2->getNumberOfTuples()
3228 * \throw If \a a1->getNumberOfComponents() != 3
3229 * \throw If \a a2->getNumberOfComponents() != 3
3231 DataArrayDouble *DataArrayDouble::CrossProduct(const DataArrayDouble *a1, const DataArrayDouble *a2)
3234 throw INTERP_KERNEL::Exception("DataArrayDouble::CrossProduct : input DataArrayDouble instance is NULL !");
3235 std::size_t nbOfComp(a1->getNumberOfComponents());
3236 if(nbOfComp!=a2->getNumberOfComponents())
3237 throw INTERP_KERNEL::Exception("Nb of components mismatch for array crossProduct !");
3239 throw INTERP_KERNEL::Exception("Nb of components must be equal to 3 for array crossProduct !");
3240 mcIdType nbOfTuple(a1->getNumberOfTuples());
3241 if(nbOfTuple!=a2->getNumberOfTuples())
3242 throw INTERP_KERNEL::Exception("Nb of tuples mismatch for array crossProduct !");
3243 DataArrayDouble *ret=DataArrayDouble::New();
3244 ret->alloc(nbOfTuple,3);
3245 double *retPtr=ret->getPointer();
3246 const double *a1Ptr(a1->begin()),*a2Ptr(a2->begin());
3247 for(mcIdType i=0;i<nbOfTuple;i++)
3249 retPtr[3*i]=a1Ptr[3*i+1]*a2Ptr[3*i+2]-a1Ptr[3*i+2]*a2Ptr[3*i+1];
3250 retPtr[3*i+1]=a1Ptr[3*i+2]*a2Ptr[3*i]-a1Ptr[3*i]*a2Ptr[3*i+2];
3251 retPtr[3*i+2]=a1Ptr[3*i]*a2Ptr[3*i+1]-a1Ptr[3*i+1]*a2Ptr[3*i];
3253 ret->copyStringInfoFrom(*a1);
3258 * Returns a new DataArrayDouble containing maximal values of two given arrays.
3259 * Info on components is copied from the first of the given arrays.
3260 * Number of tuples and components in the given arrays must be the same.
3261 * \param [in] a1 - an array to compare values with another one.
3262 * \param [in] a2 - another array to compare values with the first one.
3263 * \return DataArrayDouble * - the new instance of DataArrayDouble.
3264 * The caller is to delete this result array using decrRef() as it is no more
3266 * \throw If either \a a1 or \a a2 is NULL.
3267 * \throw If \a a1->getNumberOfTuples() != \a a2->getNumberOfTuples()
3268 * \throw If \a a1->getNumberOfComponents() != \a a2->getNumberOfComponents()
3270 DataArrayDouble *DataArrayDouble::Max(const DataArrayDouble *a1, const DataArrayDouble *a2)
3273 throw INTERP_KERNEL::Exception("DataArrayDouble::Max : input DataArrayDouble instance is NULL !");
3274 std::size_t nbOfComp(a1->getNumberOfComponents());
3275 if(nbOfComp!=a2->getNumberOfComponents())
3276 throw INTERP_KERNEL::Exception("Nb of components mismatch for array Max !");
3277 mcIdType nbOfTuple(a1->getNumberOfTuples());
3278 if(nbOfTuple!=a2->getNumberOfTuples())
3279 throw INTERP_KERNEL::Exception("Nb of tuples mismatch for array Max !");
3280 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
3281 ret->alloc(nbOfTuple,nbOfComp);
3282 double *retPtr(ret->getPointer());
3283 const double *a1Ptr(a1->begin()),*a2Ptr(a2->begin());
3284 std::size_t nbElem(nbOfTuple*nbOfComp);
3285 for(std::size_t i=0;i<nbElem;i++)
3286 retPtr[i]=std::max(a1Ptr[i],a2Ptr[i]);
3287 ret->copyStringInfoFrom(*a1);
3292 * Returns a new DataArrayDouble containing minimal values of two given arrays.
3293 * Info on components is copied from the first of the given arrays.
3294 * Number of tuples and components in the given arrays must be the same.
3295 * \param [in] a1 - an array to compare values with another one.
3296 * \param [in] a2 - another array to compare values with the first one.
3297 * \return DataArrayDouble * - the new instance of DataArrayDouble.
3298 * The caller is to delete this result array using decrRef() as it is no more
3300 * \throw If either \a a1 or \a a2 is NULL.
3301 * \throw If \a a1->getNumberOfTuples() != \a a2->getNumberOfTuples()
3302 * \throw If \a a1->getNumberOfComponents() != \a a2->getNumberOfComponents()
3304 DataArrayDouble *DataArrayDouble::Min(const DataArrayDouble *a1, const DataArrayDouble *a2)
3307 throw INTERP_KERNEL::Exception("DataArrayDouble::Min : input DataArrayDouble instance is NULL !");
3308 std::size_t nbOfComp(a1->getNumberOfComponents());
3309 if(nbOfComp!=a2->getNumberOfComponents())
3310 throw INTERP_KERNEL::Exception("Nb of components mismatch for array min !");
3311 mcIdType nbOfTuple(a1->getNumberOfTuples());
3312 if(nbOfTuple!=a2->getNumberOfTuples())
3313 throw INTERP_KERNEL::Exception("Nb of tuples mismatch for array min !");
3314 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
3315 ret->alloc(nbOfTuple,nbOfComp);
3316 double *retPtr(ret->getPointer());
3317 const double *a1Ptr(a1->begin()),*a2Ptr(a2->begin());
3318 std::size_t nbElem(nbOfTuple*nbOfComp);
3319 for(std::size_t i=0;i<nbElem;i++)
3320 retPtr[i]=std::min(a1Ptr[i],a2Ptr[i]);
3321 ret->copyStringInfoFrom(*a1);
3326 * Returns a new DataArrayDouble that is the result of pow of two given arrays. There are 3
3329 * \param [in] a1 - an array to pow up.
3330 * \param [in] a2 - another array to sum up.
3331 * \return DataArrayDouble * - the new instance of DataArrayDouble.
3332 * The caller is to delete this result array using decrRef() as it is no more
3334 * \throw If either \a a1 or \a a2 is NULL.
3335 * \throw If \a a1->getNumberOfTuples() != \a a2->getNumberOfTuples()
3336 * \throw If \a a1->getNumberOfComponents() != 1 or \a a2->getNumberOfComponents() != 1.
3337 * \throw If there is a negative value in \a a1.
3339 DataArrayDouble *DataArrayDouble::Pow(const DataArrayDouble *a1, const DataArrayDouble *a2)
3342 throw INTERP_KERNEL::Exception("DataArrayDouble::Pow : at least one of input instances is null !");
3343 mcIdType nbOfTuple=a1->getNumberOfTuples();
3344 mcIdType nbOfTuple2=a2->getNumberOfTuples();
3345 std::size_t nbOfComp=a1->getNumberOfComponents();
3346 std::size_t nbOfComp2=a2->getNumberOfComponents();
3347 if(nbOfTuple!=nbOfTuple2)
3348 throw INTERP_KERNEL::Exception("DataArrayDouble::Pow : number of tuples mismatches !");
3349 if(nbOfComp!=1 || nbOfComp2!=1)
3350 throw INTERP_KERNEL::Exception("DataArrayDouble::Pow : number of components of both arrays must be equal to 1 !");
3351 MCAuto<DataArrayDouble> ret=DataArrayDouble::New(); ret->alloc(nbOfTuple,1);
3352 const double *ptr1(a1->begin()),*ptr2(a2->begin());
3353 double *ptr=ret->getPointer();
3354 for(mcIdType i=0;i<nbOfTuple;i++,ptr1++,ptr2++,ptr++)
3358 *ptr=pow(*ptr1,*ptr2);
3362 std::ostringstream oss; oss << "DataArrayDouble::Pow : on tuple #" << i << " of a1 value is < 0 (" << *ptr1 << ") !";
3363 throw INTERP_KERNEL::Exception(oss.str().c_str());
3370 * Apply pow on values of another DataArrayDouble to values of \a this one.
3372 * \param [in] other - an array to pow to \a this one.
3373 * \throw If \a other is NULL.
3374 * \throw If \a this->getNumberOfTuples() != \a other->getNumberOfTuples()
3375 * \throw If \a this->getNumberOfComponents() != 1 or \a other->getNumberOfComponents() != 1
3376 * \throw If there is a negative value in \a this.
3378 void DataArrayDouble::powEqual(const DataArrayDouble *other)
3381 throw INTERP_KERNEL::Exception("DataArrayDouble::powEqual : input instance is null !");
3382 mcIdType nbOfTuple=getNumberOfTuples();
3383 mcIdType nbOfTuple2=other->getNumberOfTuples();
3384 std::size_t nbOfComp=getNumberOfComponents();
3385 std::size_t nbOfComp2=other->getNumberOfComponents();
3386 if(nbOfTuple!=nbOfTuple2)
3387 throw INTERP_KERNEL::Exception("DataArrayDouble::powEqual : number of tuples mismatches !");
3388 if(nbOfComp!=1 || nbOfComp2!=1)
3389 throw INTERP_KERNEL::Exception("DataArrayDouble::powEqual : number of components of both arrays must be equal to 1 !");
3390 double *ptr=getPointer();
3391 const double *ptrc=other->begin();
3392 for(mcIdType i=0;i<nbOfTuple;i++,ptrc++,ptr++)
3395 *ptr=pow(*ptr,*ptrc);
3398 std::ostringstream oss; oss << "DataArrayDouble::powEqual : on tuple #" << i << " of this value is < 0 (" << *ptr << ") !";
3399 throw INTERP_KERNEL::Exception(oss.str().c_str());
3406 * This method is \b NOT wrapped into python because it can be useful only for performance reasons in C++ context.
3407 * All values in \a this must be 0. or 1. within eps error. 0 means false, 1 means true.
3408 * If an another value than 0 or 1 appear (within eps precision) an INTERP_KERNEL::Exception will be thrown.
3410 * \throw if \a this is not allocated.
3411 * \throw if \a this has not exactly one component.
3413 std::vector<bool> DataArrayDouble::toVectorOfBool(double eps) const
3416 if(getNumberOfComponents()!=1)
3417 throw INTERP_KERNEL::Exception("DataArrayDouble::toVectorOfBool : must be applied on single component array !");
3418 mcIdType nbt(getNumberOfTuples());
3419 std::vector<bool> ret(nbt);
3420 const double *pt(begin());
3421 for(mcIdType i=0;i<nbt;i++)
3425 else if(fabs(pt[i]-1.)<eps)
3429 std::ostringstream oss; oss << "DataArrayDouble::toVectorOfBool : the tuple #" << i << " has value " << pt[i] << " is invalid ! must be 0. or 1. !";
3430 throw INTERP_KERNEL::Exception(oss.str().c_str());
3437 * Useless method for end user. Only for MPI/Corba/File serialsation for multi arrays class.
3440 void DataArrayDouble::getTinySerializationIntInformation(std::vector<mcIdType>& tinyInfo) const
3445 tinyInfo[0]=getNumberOfTuples();
3446 tinyInfo[1]=ToIdType(getNumberOfComponents());
3456 * Useless method for end user. Only for MPI/Corba/File serialsation for multi arrays class.
3459 void DataArrayDouble::getTinySerializationStrInformation(std::vector<std::string>& tinyInfo) const
3463 std::size_t nbOfCompo(getNumberOfComponents());
3464 tinyInfo.resize(nbOfCompo+1);
3465 tinyInfo[0]=getName();
3466 for(std::size_t i=0;i<nbOfCompo;i++)
3467 tinyInfo[i+1]=getInfoOnComponent(i);
3472 tinyInfo[0]=getName();
3477 * Useless method for end user. Only for MPI/Corba/File serialsation for multi arrays class.
3478 * This method returns if a feeding is needed.
3480 bool DataArrayDouble::resizeForUnserialization(const std::vector<mcIdType>& tinyInfoI)
3482 mcIdType nbOfTuple=tinyInfoI[0];
3483 mcIdType nbOfComp=tinyInfoI[1];
3484 if(nbOfTuple!=-1 || nbOfComp!=-1)
3486 alloc(nbOfTuple,nbOfComp);
3493 * Useless method for end user. Only for MPI/Corba/File serialsation for multi arrays class.
3495 void DataArrayDouble::finishUnserialization(const std::vector<mcIdType>& tinyInfoI, const std::vector<std::string>& tinyInfoS)
3497 setName(tinyInfoS[0]);
3500 std::size_t nbOfCompo(getNumberOfComponents());
3501 for(std::size_t i=0;i<nbOfCompo;i++)
3502 setInfoOnComponent(i,tinyInfoS[i+1]);
3507 * Low static method that operates 3D rotation of 'nbNodes' 3D nodes whose coordinates are arranged in \a coordsIn
3508 * around an axe ( \a center, \a vect) and with angle \a angle.
3510 void DataArrayDouble::Rotate3DAlg(const double *center, const double *vect, double angle, mcIdType nbNodes, const double *coordsIn, double *coordsOut)
3512 if(!center || !vect)
3513 throw INTERP_KERNEL::Exception("DataArrayDouble::Rotate3DAlg : null vector in input !");
3514 double sina(sin(angle));
3515 double cosa(cos(angle));
3516 double vectorNorm[3];
3518 double matrixTmp[9];
3519 double norm(sqrt(vect[0]*vect[0]+vect[1]*vect[1]+vect[2]*vect[2]));
3520 if(norm<std::numeric_limits<double>::min())
3521 throw INTERP_KERNEL::Exception("DataArrayDouble::Rotate3DAlg : magnitude of input vector is too close of 0. !");
3522 std::transform(vect,vect+3,vectorNorm,std::bind(std::multiplies<double>(),std::placeholders::_1,1/norm));
3523 //rotation matrix computation
3524 matrix[0]=cosa; matrix[1]=0.; matrix[2]=0.; matrix[3]=0.; matrix[4]=cosa; matrix[5]=0.; matrix[6]=0.; matrix[7]=0.; matrix[8]=cosa;
3525 matrixTmp[0]=vectorNorm[0]*vectorNorm[0]; matrixTmp[1]=vectorNorm[0]*vectorNorm[1]; matrixTmp[2]=vectorNorm[0]*vectorNorm[2];
3526 matrixTmp[3]=vectorNorm[1]*vectorNorm[0]; matrixTmp[4]=vectorNorm[1]*vectorNorm[1]; matrixTmp[5]=vectorNorm[1]*vectorNorm[2];
3527 matrixTmp[6]=vectorNorm[2]*vectorNorm[0]; matrixTmp[7]=vectorNorm[2]*vectorNorm[1]; matrixTmp[8]=vectorNorm[2]*vectorNorm[2];
3528 std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind(std::multiplies<double>(),std::placeholders::_1,1-cosa));
3529 std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
3530 matrixTmp[0]=0.; matrixTmp[1]=-vectorNorm[2]; matrixTmp[2]=vectorNorm[1];
3531 matrixTmp[3]=vectorNorm[2]; matrixTmp[4]=0.; matrixTmp[5]=-vectorNorm[0];
3532 matrixTmp[6]=-vectorNorm[1]; matrixTmp[7]=vectorNorm[0]; matrixTmp[8]=0.;
3533 std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind(std::multiplies<double>(),std::placeholders::_1,sina));
3534 std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
3535 //rotation matrix computed.
3537 for(mcIdType i=0; i<nbNodes; i++)
3539 std::transform(coordsIn+i*3,coordsIn+(i+1)*3,center,tmp,std::minus<double>());
3540 coordsOut[i*3]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+matrix[2]*tmp[2]+center[0];
3541 coordsOut[i*3+1]=matrix[3]*tmp[0]+matrix[4]*tmp[1]+matrix[5]*tmp[2]+center[1];
3542 coordsOut[i*3+2]=matrix[6]*tmp[0]+matrix[7]*tmp[1]+matrix[8]*tmp[2]+center[2];
3546 void DataArrayDouble::Symmetry3DPlane(const double point[3], const double normalVector[3], mcIdType nbNodes, const double *coordsIn, double *coordsOut)
3548 double matrix[9],matrix2[9],matrix3[9];
3549 double vect[3],crossVect[3];
3550 INTERP_KERNEL::orthogonalVect3(normalVector,vect);
3551 crossVect[0]=normalVector[1]*vect[2]-normalVector[2]*vect[1];
3552 crossVect[1]=normalVector[2]*vect[0]-normalVector[0]*vect[2];
3553 crossVect[2]=normalVector[0]*vect[1]-normalVector[1]*vect[0];
3554 double nv(INTERP_KERNEL::norm<3>(vect)),ni(INTERP_KERNEL::norm<3>(normalVector)),nc(INTERP_KERNEL::norm<3>(crossVect));
3555 matrix[0]=vect[0]/nv; matrix[1]=crossVect[0]/nc; matrix[2]=-normalVector[0]/ni;
3556 matrix[3]=vect[1]/nv; matrix[4]=crossVect[1]/nc; matrix[5]=-normalVector[1]/ni;
3557 matrix[6]=vect[2]/nv; matrix[7]=crossVect[2]/nc; matrix[8]=-normalVector[2]/ni;
3558 matrix2[0]=vect[0]/nv; matrix2[1]=vect[1]/nv; matrix2[2]=vect[2]/nv;
3559 matrix2[3]=crossVect[0]/nc; matrix2[4]=crossVect[1]/nc; matrix2[5]=crossVect[2]/nc;
3560 matrix2[6]=normalVector[0]/ni; matrix2[7]=normalVector[1]/ni; matrix2[8]=normalVector[2]/ni;
3561 for(mcIdType i=0;i<3;i++)
3562 for(mcIdType j=0;j<3;j++)
3565 for(mcIdType k=0;k<3;k++)
3566 val+=matrix[3*i+k]*matrix2[3*k+j];
3569 //rotation matrix computed.
3571 for(mcIdType i=0; i<nbNodes; i++)
3573 std::transform(coordsIn+i*3,coordsIn+(i+1)*3,point,tmp,std::minus<double>());
3574 coordsOut[i*3]=matrix3[0]*tmp[0]+matrix3[1]*tmp[1]+matrix3[2]*tmp[2]+point[0];
3575 coordsOut[i*3+1]=matrix3[3]*tmp[0]+matrix3[4]*tmp[1]+matrix3[5]*tmp[2]+point[1];
3576 coordsOut[i*3+2]=matrix3[6]*tmp[0]+matrix3[7]*tmp[1]+matrix3[8]*tmp[2]+point[2];
3580 void DataArrayDouble::GiveBaseForPlane(const double normalVector[3], double baseOfPlane[9])
3582 double vect[3],crossVect[3];
3583 INTERP_KERNEL::orthogonalVect3(normalVector,vect);
3584 crossVect[0]=normalVector[1]*vect[2]-normalVector[2]*vect[1];
3585 crossVect[1]=normalVector[2]*vect[0]-normalVector[0]*vect[2];
3586 crossVect[2]=normalVector[0]*vect[1]-normalVector[1]*vect[0];
3587 double nv(INTERP_KERNEL::norm<3>(vect)),ni(INTERP_KERNEL::norm<3>(normalVector)),nc(INTERP_KERNEL::norm<3>(crossVect));
3588 baseOfPlane[0]=vect[0]/nv; baseOfPlane[1]=vect[1]/nv; baseOfPlane[2]=vect[2]/nv;
3589 baseOfPlane[3]=crossVect[0]/nc; baseOfPlane[4]=crossVect[1]/nc; baseOfPlane[5]=crossVect[2]/nc;
3590 baseOfPlane[6]=normalVector[0]/ni; baseOfPlane[7]=normalVector[1]/ni; baseOfPlane[8]=normalVector[2]/ni;
3594 * \param [in] seg2 : coordinates of input seg2 expected to have spacedim==2
3595 * \param [in] tri3 : coordinates of input tri3 also expected to have spacedim==2
3596 * \param [out] coeffs : the result of integration normalized to 1. along \a seg2 inside tri3 sorted by the node id of \a tri3
3597 * \param [out] length : the length of seg2. That is too say the length of integration
3599 void DataArrayDouble::ComputeIntegralOfSeg2IntoTri3(const double seg2[4], const double tri3[6], double coeffs[3], double& length)
3601 length=INTERP_KERNEL::norme_vecteur(seg2,seg2+2);
3603 INTERP_KERNEL::mid_of_seg2(seg2,seg2+2,mid);
3604 INTERP_KERNEL::barycentric_coords<2>(tri3,mid,coeffs); // integral along seg2 is equal to value at the center of SEG2 !
3608 * Low static method that operates 3D rotation of \a nbNodes 3D nodes whose coordinates are arranged in \a coords
3609 * around the center point \a center and with angle \a angle.
3611 void DataArrayDouble::Rotate2DAlg(const double *center, double angle, mcIdType nbNodes, const double *coordsIn, double *coordsOut)
3613 double cosa=cos(angle);
3614 double sina=sin(angle);
3616 matrix[0]=cosa; matrix[1]=-sina; matrix[2]=sina; matrix[3]=cosa;
3618 for(mcIdType i=0; i<nbNodes; i++)
3620 std::transform(coordsIn+i*2,coordsIn+(i+1)*2,center,tmp,std::minus<double>());
3621 coordsOut[i*2]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+center[0];
3622 coordsOut[i*2+1]=matrix[2]*tmp[0]+matrix[3]*tmp[1]+center[1];
3626 DataArrayDoubleIterator::DataArrayDoubleIterator(DataArrayDouble *da):DataArrayIterator<double>(da)
3630 DataArrayDoubleTuple::DataArrayDoubleTuple(double *pt, std::size_t nbOfComp):DataArrayTuple<double>(pt,nbOfComp)
3635 std::string DataArrayDoubleTuple::repr() const
3637 std::ostringstream oss; oss.precision(17); oss << "(";
3638 for(std::size_t i=0;i<_nb_of_compo-1;i++)
3639 oss << _pt[i] << ", ";
3640 oss << _pt[_nb_of_compo-1] << ")";
3644 double DataArrayDoubleTuple::doubleValue() const
3646 return this->zeValue();
3650 * This method returns a newly allocated instance the caller should dealed with by a MEDCoupling::DataArrayDouble::decrRef.
3651 * This method performs \b no copy of data. The content is only referenced using MEDCoupling::DataArrayDouble::useArray with ownership set to \b false.
3652 * This method throws an INTERP_KERNEL::Exception is it is impossible to match sizes of \b this that is too say \b nbOfCompo=this->_nb_of_elem and \bnbOfTuples==1 or
3653 * \b nbOfCompo=1 and \bnbOfTuples==this->_nb_of_elem.
3655 DataArrayDouble *DataArrayDoubleTuple::buildDADouble(std::size_t nbOfTuples, std::size_t nbOfCompo) const
3657 return this->buildDA(nbOfTuples,nbOfCompo);
3661 * Returns a full copy of \a this. For more info on copying data arrays see
3662 * \ref MEDCouplingArrayBasicsCopyDeep.
3663 * \return DataArrayInt * - a new instance of DataArrayInt.
3665 DataArrayInt32 *DataArrayInt32::deepCopy() const
3667 return new DataArrayInt32(*this);
3670 MCAuto<DataArrayInt64> DataArrayInt32::convertToInt64Arr() const
3672 this->checkAllocated();
3673 MCAuto<DataArrayInt64> ret(DataArrayInt64::New());
3674 ret->alloc(this->getNumberOfTuples(),this->getNumberOfComponents());
3675 ret->copyStringInfoFrom(*this);
3676 const std::int32_t *pt(this->begin());
3677 std::for_each(ret->getPointer(),ret->getPointer()+ret->getNbOfElems(),[&pt](std::int64_t& val) { val = std::int64_t(*pt++); });
3681 DataArrayInt32Iterator *DataArrayInt32::iterator()
3683 return new DataArrayInt32Iterator(this);
3687 DataArrayInt32Iterator::DataArrayInt32Iterator(DataArrayInt32 *da):DataArrayIterator<Int32>(da)
3691 DataArrayInt32Tuple::DataArrayInt32Tuple(Int32 *pt, std::size_t nbOfComp):DataArrayTuple<Int32>(pt,nbOfComp)
3695 std::string DataArrayInt32Tuple::repr() const
3697 std::ostringstream oss; oss << "(";
3698 for(std::size_t i=0;i<_nb_of_compo-1;i++)
3699 oss << _pt[i] << ", ";
3700 oss << _pt[_nb_of_compo-1] << ")";
3704 Int32 DataArrayInt32Tuple::intValue() const
3706 return this->zeValue();
3710 * This method returns a newly allocated instance the caller should dealed with by a MEDCoupling::DataArrayInt::decrRef.
3711 * This method performs \b no copy of data. The content is only referenced using MEDCoupling::DataArrayInt::useArray with ownership set to \b false.
3712 * This method throws an INTERP_KERNEL::Exception is it is impossible to match sizes of \b this that is too say \b nbOfCompo=this->_nb_of_elem and \bnbOfTuples==1 or
3713 * \b nbOfCompo=1 and \bnbOfTuples==this->_nb_of_elem.
3715 DataArrayInt32 *DataArrayInt32Tuple::buildDAInt(std::size_t nbOfTuples, std::size_t nbOfCompo) const
3717 return this->buildDA(nbOfTuples,nbOfCompo);
3720 MCAuto<DataArrayInt32> DataArrayInt64::convertToInt32Arr() const
3722 this->checkAllocated();
3723 MCAuto<DataArrayInt32> ret(DataArrayInt32::New());
3724 ret->alloc(this->getNumberOfTuples(),this->getNumberOfComponents());
3725 ret->copyStringInfoFrom(*this);
3726 const std::int64_t *pt(this->begin());
3727 std::for_each(ret->getPointer(),ret->getPointer()+ret->getNbOfElems(),[&pt](std::int32_t& val) { val = std::int32_t(*pt++); });
3731 DataArrayInt64Iterator *DataArrayInt64::iterator()
3733 return new DataArrayInt64Iterator(this);
3737 DataArrayInt64Iterator::DataArrayInt64Iterator(DataArrayInt64 *da):DataArrayIterator<Int64>(da)
3741 DataArrayInt64Tuple::DataArrayInt64Tuple(Int64 *pt, std::size_t nbOfComp):DataArrayTuple<Int64>(pt,nbOfComp)
3745 std::string DataArrayInt64Tuple::repr() const
3747 std::ostringstream oss; oss << "(";
3748 for(std::size_t i=0;i<_nb_of_compo-1;i++)
3749 oss << _pt[i] << ", ";
3750 oss << _pt[_nb_of_compo-1] << ")";
3754 Int64 DataArrayInt64Tuple::intValue() const
3756 return this->zeValue();
3759 DataArrayInt64 *DataArrayInt64Tuple::buildDAInt(std::size_t nbOfTuples, std::size_t nbOfCompo) const
3761 return this->buildDA(nbOfTuples,nbOfCompo);
3765 DataArrayInt64 *DataArrayInt64::deepCopy() const
3767 return new DataArrayInt64(*this);