1 // Copyright (C) 2007-2021 CEA/DEN, EDF R&D
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License, or (at your option) any later version.
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 // Lesser General Public License for more details.
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // Author : Anthony Geay (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 * Returns the var part of the full component information.
428 * For example, if \a info == "SIGXY [N/m^2]", then this method returns "SIGXY".
429 * If a unit part of information is not detected by presence of
430 * two square brackets, then the whole \a info is returned.
431 * To read more about the component information format, see
432 * \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
433 * \param [in] info - the full component information.
434 * \return std::string - a string containing only var information, or the \a info.
436 std::string DataArray::GetVarNameFromInfo(const std::string& info)
438 std::size_t p1=info.find_last_of('[');
439 std::size_t p2=info.find_last_of(']');
440 if(p1==std::string::npos || p2==std::string::npos)
445 return std::string();
446 std::size_t p3=info.find_last_not_of(' ',p1-1);
447 return info.substr(0,p3+1);
451 * Returns the unit part of the full component information.
452 * For example, if \a info == "SIGXY [ N/m^2]", then this method returns " N/m^2".
453 * If a unit part of information is not detected by presence of
454 * two square brackets, then an empty string is returned.
455 * To read more about the component information format, see
456 * \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
457 * \param [in] info - the full component information.
458 * \return std::string - a string containing only unit information, if any, or "".
460 std::string DataArray::GetUnitFromInfo(const std::string& info)
462 std::size_t p1=info.find_last_of('[');
463 std::size_t p2=info.find_last_of(']');
464 if(p1==std::string::npos || p2==std::string::npos)
465 return std::string();
467 return std::string();
468 return info.substr(p1+1,p2-p1-1);
472 * This method put in info format the result of the merge of \a var and \a unit.
473 * The standard format for that is "var [unit]".
474 * Inversely you can retrieve the var part or the unit part of info string using resp. GetVarNameFromInfo and GetUnitFromInfo.
476 std::string DataArray::BuildInfoFromVarAndUnit(const std::string& var, const std::string& unit)
478 std::ostringstream oss;
479 oss << var << " [" << unit << "]";
483 std::string DataArray::GetAxisTypeRepr(MEDCouplingAxisType at)
488 return std::string("AX_CART");
490 return std::string("AX_CYL");
492 return std::string("AX_SPHER");
494 throw INTERP_KERNEL::Exception("DataArray::GetAxisTypeRepr : unrecognized axis type enum !");
499 * Returns a new DataArray by concatenating all given arrays, so that (1) the number
500 * of tuples in the result array is a sum of the number of tuples of given arrays and (2)
501 * the number of component in the result array is same as that of each of given arrays.
502 * Info on components is copied from the first of the given arrays. Number of components
503 * in the given arrays must be the same.
504 * \param [in] arrs - a sequence of arrays to include in the result array. All arrays must have the same type.
505 * \return DataArray * - the new instance of DataArray (that can be either DataArrayInt, DataArrayDouble, DataArrayChar).
506 * The caller is to delete this result array using decrRef() as it is no more
508 * \throw If all arrays within \a arrs are NULL.
509 * \throw If all not null arrays in \a arrs have not the same type.
510 * \throw If getNumberOfComponents() of arrays within \a arrs.
512 DataArray *DataArray::Aggregate(const std::vector<const DataArray *>& arrs)
514 std::vector<const DataArray *> arr2;
515 for(std::vector<const DataArray *>::const_iterator it=arrs.begin();it!=arrs.end();it++)
519 throw INTERP_KERNEL::Exception("DataArray::Aggregate : only null instance in input vector !");
520 std::vector<const DataArrayDouble *> arrd;
521 std::vector<const DataArrayIdType *> arri;
522 std::vector<const DataArrayChar *> arrc;
523 for(std::vector<const DataArray *>::const_iterator it=arr2.begin();it!=arr2.end();it++)
525 const DataArrayDouble *a=dynamic_cast<const DataArrayDouble *>(*it);
527 { arrd.push_back(a); continue; }
528 const DataArrayIdType *b=dynamic_cast<const DataArrayIdType *>(*it);
530 { arri.push_back(b); continue; }
531 const DataArrayChar *c=dynamic_cast<const DataArrayChar *>(*it);
533 { arrc.push_back(c); continue; }
534 throw INTERP_KERNEL::Exception("DataArray::Aggregate : presence of not null instance in inuput that is not in [DataArrayDouble, DataArrayInt, DataArrayChar] !");
536 if(arr2.size()==arrd.size())
537 return DataArrayDouble::Aggregate(arrd);
538 if(arr2.size()==arri.size())
539 return DataArrayIdType::Aggregate(arri);
540 if(arr2.size()==arrc.size())
541 return DataArrayChar::Aggregate(arrc);
542 throw INTERP_KERNEL::Exception("DataArray::Aggregate : all input arrays must have the same type !");
546 * Sets information on a component specified by an index.
547 * To know more on format of this information
548 * see \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
549 * \warning Don't pass NULL as \a info!
550 * \param [in] i - the index (zero based) of the component of interest.
551 * \param [in] info - the string containing the information.
552 * \throw If \a i is not a valid component index.
554 void DataArray::setInfoOnComponent(std::size_t i, const std::string& info)
556 if(i<_info_on_compo.size())
557 _info_on_compo[i]=info;
560 std::ostringstream oss; oss << "DataArray::setInfoOnComponent : Specified component id is out of range (" << i << ") compared with nb of actual components (" << _info_on_compo.size();
561 throw INTERP_KERNEL::Exception(oss.str().c_str());
566 * Sets information on all components. This method can change number of components
567 * at certain conditions; if the conditions are not respected, an exception is thrown.
568 * The number of components can be changed in \a this only if \a this is not allocated.
569 * The condition of number of components must not be changed.
571 * To know more on format of the component information see
572 * \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
573 * \param [in] info - a vector of component infos.
574 * \throw If \a this->getNumberOfComponents() != \a info.size() && \a this->isAllocated()
576 void DataArray::setInfoAndChangeNbOfCompo(const std::vector<std::string>& info)
578 if(getNumberOfComponents()!=info.size())
584 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 !";
585 throw INTERP_KERNEL::Exception(oss.str().c_str());
592 void DataArray::checkNbOfTuples(mcIdType nbOfTuples, const std::string& msg) const
594 if(getNumberOfTuples()!=nbOfTuples)
596 std::ostringstream oss; oss << msg << " : mismatch number of tuples : expected " << nbOfTuples << " having " << getNumberOfTuples() << " !";
597 throw INTERP_KERNEL::Exception(oss.str().c_str());
601 void DataArray::checkNbOfComps(std::size_t nbOfCompo, const std::string& msg) const
603 if (getNumberOfComponents()!=nbOfCompo)
605 std::ostringstream oss; oss << msg << " : mismatch number of components : expected " << nbOfCompo << " having " << getNumberOfComponents() << " !";
606 throw INTERP_KERNEL::Exception(oss.str().c_str());
610 void DataArray::checkNbOfElems(mcIdType nbOfElems, const std::string& msg) const
612 if(getNbOfElems()!=nbOfElems)
614 std::ostringstream oss; oss << msg << " : mismatch number of elems : Expected " << nbOfElems << " having " << getNbOfElems() << " !";
615 throw INTERP_KERNEL::Exception(oss.str().c_str());
619 void DataArray::checkNbOfTuplesAndComp(const DataArray& other, const std::string& msg) const
621 if(getNumberOfTuples()!=other.getNumberOfTuples())
623 std::ostringstream oss; oss << msg << " : mismatch number of tuples : expected " << other.getNumberOfTuples() << " having " << getNumberOfTuples() << " !";
624 throw INTERP_KERNEL::Exception(oss.str().c_str());
626 if(getNumberOfComponents()!=other.getNumberOfComponents())
628 std::ostringstream oss; oss << msg << " : mismatch number of components : expected " << other.getNumberOfComponents() << " having " << getNumberOfComponents() << " !";
629 throw INTERP_KERNEL::Exception(oss.str().c_str());
633 void DataArray::checkNbOfTuplesAndComp(mcIdType nbOfTuples, std::size_t nbOfCompo, const std::string& msg) const
635 checkNbOfTuples(nbOfTuples,msg);
636 checkNbOfComps(nbOfCompo,msg);
640 * Simply this method checks that \b value is in [0,\b ref).
642 void DataArray::CheckValueInRange(mcIdType ref, mcIdType value, const std::string& msg)
644 if(value<0 || value>=ref)
646 std::ostringstream oss; oss << "DataArray::CheckValueInRange : " << msg << " ! Expected in range [0," << ref << "[ having " << value << " !";
647 throw INTERP_KERNEL::Exception(oss.str().c_str());
652 * This method checks that [\b start, \b end) is compliant with ref length \b value.
653 * typically start in [0,\b value) and end in [0,\b value). If value==start and start==end, it is supported.
655 void DataArray::CheckValueInRangeEx(mcIdType value, mcIdType start, mcIdType end, const std::string& msg)
657 if(start<0 || start>=value)
659 if(value!=start || end!=start)
661 std::ostringstream oss; oss << "DataArray::CheckValueInRangeEx : " << msg << " ! Expected start " << start << " of input range, in [0," << value << "[ !";
662 throw INTERP_KERNEL::Exception(oss.str().c_str());
665 if(end<0 || end>value)
667 std::ostringstream oss; oss << "DataArray::CheckValueInRangeEx : " << msg << " ! Expected end " << end << " of input range, in [0," << value << "] !";
668 throw INTERP_KERNEL::Exception(oss.str().c_str());
672 void DataArray::CheckClosingParInRange(mcIdType ref, mcIdType value, const std::string& msg)
674 if(value<0 || value>ref)
676 std::ostringstream oss; oss << "DataArray::CheckClosingParInRange : " << msg << " ! Expected input range in [0," << ref << "] having closing open parenthesis " << value << " !";
677 throw INTERP_KERNEL::Exception(oss.str().c_str());
682 * 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,
683 * typically it is a whole slice of tuples of DataArray or cells, nodes of a mesh...
685 * The input \a sliceId should be an id in [0, \a nbOfSlices) that specifies the slice of work.
687 * \param [in] start - the start of the input slice of the whole work to perform split into slices.
688 * \param [in] stop - the stop of the input slice of the whole work to perform split into slices.
689 * \param [in] step - the step (that can be <0) of the input slice of the whole work to perform split into slices.
690 * \param [in] sliceId - the slice id considered
691 * \param [in] nbOfSlices - the number of slices (typically the number of cores on which the work is expected to be sliced)
692 * \param [out] startSlice - the start of the slice considered
693 * \param [out] stopSlice - the stop of the slice consided
695 * \throw If \a step == 0
696 * \throw If \a nbOfSlices not > 0
697 * \throw If \a sliceId not in [0,nbOfSlices)
699 void DataArray::GetSlice(mcIdType start, mcIdType stop, mcIdType step, mcIdType sliceId, mcIdType nbOfSlices, mcIdType& startSlice, mcIdType& stopSlice)
701 DataArrayTools<mcIdType>::GetSlice(start, stop, step, sliceId, nbOfSlices, startSlice, stopSlice);
704 mcIdType DataArray::GetNumberOfItemGivenBES(mcIdType begin, mcIdType end, mcIdType step, const std::string& msg)
706 return DataArrayTools<mcIdType>::GetNumberOfItemGivenBES(begin, end, step, msg);
709 mcIdType DataArray::GetNumberOfItemGivenBESRelative(mcIdType begin, mcIdType end, mcIdType step, const std::string& msg)
711 return DataArrayTools<mcIdType>::GetNumberOfItemGivenBESRelative(begin, end, step, msg);
714 mcIdType DataArray::GetPosOfItemGivenBESRelativeNoThrow(mcIdType value, mcIdType begin, mcIdType end, mcIdType step)
716 return DataArrayTools<mcIdType>::GetPosOfItemGivenBESRelativeNoThrow(value, begin, end, step);
720 * Returns a new instance of DataArrayDouble. The caller is to delete this array
721 * using decrRef() as it is no more needed.
723 DataArrayDouble *DataArrayDouble::New()
725 return new DataArrayDouble;
729 * Returns the only one value in \a this, if and only if number of elements
730 * (nb of tuples * nb of components) is equal to 1, and that \a this is allocated.
731 * \return double - the sole value stored in \a this array.
732 * \throw If at least one of conditions stated above is not fulfilled.
734 double DataArrayDouble::doubleValue() const
738 if(getNbOfElems()==1)
740 return *getConstPointer();
743 throw INTERP_KERNEL::Exception("DataArrayDouble::doubleValue : DataArrayDouble instance is allocated but number of elements is not equal to 1 !");
746 throw INTERP_KERNEL::Exception("DataArrayDouble::doubleValue : DataArrayDouble instance is not allocated !");
750 * Returns a full copy of \a this. For more info on copying data arrays see
751 * \ref MEDCouplingArrayBasicsCopyDeep.
752 * \return DataArrayDouble * - a new instance of DataArrayDouble. The caller is to
753 * delete this array using decrRef() as it is no more needed.
755 DataArrayDouble *DataArrayDouble::deepCopy() const
757 return new DataArrayDouble(*this);
761 * Checks that \a this array is consistently **increasing** or **decreasing** in value,
762 * with at least absolute difference value of |\a eps| at each step.
763 * If not an exception is thrown.
764 * \param [in] increasing - if \a true, the array values should be increasing.
765 * \param [in] eps - minimal absolute difference between the neighbor values at which
766 * the values are considered different.
767 * \throw If sequence of values is not strictly monotonic in agreement with \a
769 * \throw If \a this->getNumberOfComponents() != 1.
770 * \throw If \a this is not allocated.
772 void DataArrayDouble::checkMonotonic(bool increasing, double eps) const
774 if(!isMonotonic(increasing,eps))
777 throw INTERP_KERNEL::Exception("DataArrayDouble::checkMonotonic : 'this' is not INCREASING monotonic !");
779 throw INTERP_KERNEL::Exception("DataArrayDouble::checkMonotonic : 'this' is not DECREASING monotonic !");
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 * \param [in] increasing - if \a true, array values should be increasing.
787 * \param [in] eps - minimal absolute difference between the neighbor values at which
788 * the values are considered different.
789 * \return bool - \a true if values change in accordance with \a increasing arg.
790 * \throw If \a this->getNumberOfComponents() != 1.
791 * \throw If \a this is not allocated.
793 bool DataArrayDouble::isMonotonic(bool increasing, double eps) const
796 if(getNumberOfComponents()!=1)
797 throw INTERP_KERNEL::Exception("DataArrayDouble::isMonotonic : only supported with 'this' array with ONE component !");
798 mcIdType nbOfElements(getNumberOfTuples());
799 const double *ptr=getConstPointer();
803 double absEps=fabs(eps);
806 for(mcIdType i=1;i<nbOfElements;i++)
808 if(ptr[i]<(ref+absEps))
816 for(mcIdType i=1;i<nbOfElements;i++)
818 if(ptr[i]>(ref-absEps))
826 void DataArrayDouble::writeVTK(std::ostream& ofs, mcIdType indent, const std::string& nameInFile, DataArrayByte *byteArr) const
828 static const char SPACE[4]={' ',' ',' ',' '};
830 std::string idt(indent,' ');
832 ofs << idt << "<DataArray type=\"Float32\" Name=\"" << nameInFile << "\" NumberOfComponents=\"" << getNumberOfComponents() << "\"";
834 bool areAllEmpty(true);
835 for(std::vector<std::string>::const_iterator it=_info_on_compo.begin();it!=_info_on_compo.end();it++)
839 for(std::size_t i=0;i<_info_on_compo.size();i++)
840 ofs << " ComponentName" << i << "=\"" << _info_on_compo[i] << "\"";
844 ofs << " format=\"appended\" offset=\"" << byteArr->getNumberOfTuples() << "\">";
845 INTERP_KERNEL::AutoPtr<float> tmp(new float[getNbOfElems()]);
847 // to make Visual C++ happy : instead of std::copy(begin(),end(),(float *)tmp);
848 for(const double *src=begin();src!=end();src++,pt++)
850 const char *data(reinterpret_cast<const char *>((float *)tmp));
851 std::size_t sz(getNbOfElems()*sizeof(float));
852 byteArr->insertAtTheEnd(data,data+sz);
853 byteArr->insertAtTheEnd(SPACE,SPACE+4);
857 ofs << " RangeMin=\"" << getMinValueInArray() << "\" RangeMax=\"" << getMaxValueInArray() << "\" format=\"ascii\">\n" << idt;
858 std::copy(begin(),end(),std::ostream_iterator<double>(ofs," "));
860 ofs << std::endl << idt << "</DataArray>\n";
863 void DataArrayDouble::reprCppStream(const std::string& varName, std::ostream& stream) const
865 mcIdType nbTuples=getNumberOfTuples();
866 std::size_t nbComp=getNumberOfComponents();
867 const double *data(getConstPointer());
868 stream.precision(17);
869 stream << "DataArrayDouble *" << varName << "=DataArrayDouble::New();" << std::endl;
870 if(nbTuples*nbComp>=1)
872 stream << "const double " << varName << "Data[" << nbTuples*nbComp << "]={";
873 std::copy(data,data+nbTuples*nbComp-1,std::ostream_iterator<double>(stream,","));
874 stream << data[nbTuples*nbComp-1] << "};" << std::endl;
875 stream << varName << "->useArray(" << varName << "Data,false,CPP_DEALLOC," << nbTuples << "," << nbComp << ");" << std::endl;
878 stream << varName << "->alloc(" << nbTuples << "," << nbComp << ");" << std::endl;
879 stream << varName << "->setName(\"" << getName() << "\");" << std::endl;
883 * Method that gives a quick overvien of \a this for python.
885 void DataArrayDouble::reprQuickOverview(std::ostream& stream) const
887 static const std::size_t MAX_NB_OF_BYTE_IN_REPR=300;
888 stream << "DataArrayDouble C++ instance at " << this << ". ";
891 std::size_t nbOfCompo(_info_on_compo.size());
894 mcIdType nbOfTuples(getNumberOfTuples());
895 stream << "Number of tuples : " << nbOfTuples << ". Number of components : " << nbOfCompo << "." << std::endl;
896 reprQuickOverviewData(stream,MAX_NB_OF_BYTE_IN_REPR);
899 stream << "Number of components : 0.";
902 stream << "*** No data allocated ****";
905 void DataArrayDouble::reprQuickOverviewData(std::ostream& stream, std::size_t maxNbOfByteInRepr) const
907 const double *data=begin();
908 mcIdType nbOfTuples(getNumberOfTuples());
909 std::size_t nbOfCompo(_info_on_compo.size());
910 std::ostringstream oss2; oss2 << "[";
912 std::string oss2Str(oss2.str());
913 bool isFinished=true;
914 for(mcIdType i=0;i<nbOfTuples && isFinished;i++)
919 for(std::size_t j=0;j<nbOfCompo;j++,data++)
922 if(j!=nbOfCompo-1) oss2 << ", ";
928 if(i!=nbOfTuples-1) oss2 << ", ";
929 std::string oss3Str(oss2.str());
930 if(oss3Str.length()<maxNbOfByteInRepr)
942 * Equivalent to DataArrayDouble::isEqual except that if false the reason of
945 * \param [in] other the instance to be compared with \a this
946 * \param [in] prec the precision to compare numeric data of the arrays.
947 * \param [out] reason In case of inequality returns the reason.
948 * \sa DataArrayDouble::isEqual
950 bool DataArrayDouble::isEqualIfNotWhy(const DataArrayDouble& other, double prec, std::string& reason) const
952 if(!areInfoEqualsIfNotWhy(other,reason))
954 return _mem.isEqual(other._mem,prec,reason);
958 * Checks if \a this and another DataArrayDouble are fully equal. For more info see
959 * \ref MEDCouplingArrayBasicsCompare.
960 * \param [in] other - an instance of DataArrayDouble to compare with \a this one.
961 * \param [in] prec - precision value to compare numeric data of the arrays.
962 * \return bool - \a true if the two arrays are equal, \a false else.
964 bool DataArrayDouble::isEqual(const DataArrayDouble& other, double prec) const
967 return isEqualIfNotWhy(other,prec,tmp);
971 * Checks if values of \a this and another DataArrayDouble are equal. For more info see
972 * \ref MEDCouplingArrayBasicsCompare.
973 * \param [in] other - an instance of DataArrayDouble to compare with \a this one.
974 * \param [in] prec - precision value to compare numeric data of the arrays.
975 * \return bool - \a true if the values of two arrays are equal, \a false else.
977 bool DataArrayDouble::isEqualWithoutConsideringStr(const DataArrayDouble& other, double prec) const
980 return _mem.isEqual(other._mem,prec,tmp);
984 * This method checks that all tuples in \a other are in \a this.
985 * If true, the output param \a tupleIds contains the tuples ids of \a this that correspond to tupes in \a this.
986 * 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.
988 * \param [in] other - the array having the same number of components than \a this.
989 * \param [out] tupleIds - the tuple ids containing the same number of tuples than \a other has.
990 * \sa DataArrayDouble::findCommonTuples
992 bool DataArrayDouble::areIncludedInMe(const DataArrayDouble *other, double prec, DataArrayIdType *&tupleIds) const
995 throw INTERP_KERNEL::Exception("DataArrayDouble::areIncludedInMe : input array is NULL !");
996 checkAllocated(); other->checkAllocated();
997 if(getNumberOfComponents()!=other->getNumberOfComponents())
998 throw INTERP_KERNEL::Exception("DataArrayDouble::areIncludedInMe : the number of components does not match !");
999 MCAuto<DataArrayDouble> a=DataArrayDouble::Aggregate(this,other);
1000 DataArrayIdType *c=0,*ci=0;
1001 a->findCommonTuples(prec,getNumberOfTuples(),c,ci);
1002 MCAuto<DataArrayIdType> cSafe(c),ciSafe(ci);
1003 mcIdType newNbOfTuples=-1;
1004 MCAuto<DataArrayIdType> ids=DataArrayIdType::ConvertIndexArrayToO2N(a->getNumberOfTuples(),c->begin(),ci->begin(),ci->end(),newNbOfTuples);
1005 MCAuto<DataArrayIdType> ret1=ids->selectByTupleIdSafeSlice(getNumberOfTuples(),a->getNumberOfTuples(),1);
1006 tupleIds=ret1.retn();
1007 return newNbOfTuples==getNumberOfTuples();
1011 * Searches for tuples coincident within \a prec tolerance. Each tuple is considered
1012 * as coordinates of a point in getNumberOfComponents()-dimensional space. The
1013 * distance separating two points is computed with the infinite norm.
1015 * Indices of coincident tuples are stored in output arrays.
1016 * A pair of arrays (\a comm, \a commIndex) is called "Surjective Format 2".
1018 * This method is typically used by MEDCouplingPointSet::findCommonNodes() and
1019 * MEDCouplingUMesh::mergeNodes().
1020 * \param [in] prec - minimal absolute distance between two tuples (infinite norm) at which they are
1021 * considered not coincident.
1022 * \param [in] limitTupleId - limit tuple id. If all tuples within a group of coincident
1023 * tuples have id strictly lower than \a limitTupleId then they are not returned.
1024 * \param [out] comm - the array holding ids (== indices) of coincident tuples.
1025 * \a comm->getNumberOfComponents() == 1.
1026 * \a comm->getNumberOfTuples() == \a commIndex->back().
1027 * \param [out] commIndex - the array dividing all indices stored in \a comm into
1028 * groups of (indices of) coincident tuples. Its every value is a tuple
1029 * index where a next group of tuples begins. For example the second
1030 * group of tuples in \a comm is described by following range of indices:
1031 * [ \a commIndex[1], \a commIndex[2] ). \a commIndex->getNumberOfTuples()-1
1032 * gives the number of groups of coincident tuples.
1033 * \throw If \a this is not allocated.
1034 * \throw If the number of components is not in [1,2,3,4].
1036 * \if ENABLE_EXAMPLES
1037 * \ref cpp_mcdataarraydouble_findcommontuples "Here is a C++ example".
1039 * \ref py_mcdataarraydouble_findcommontuples "Here is a Python example".
1041 * \sa DataArrayInt::ConvertIndexArrayToO2N(), DataArrayDouble::areIncludedInMe
1043 void DataArrayDouble::findCommonTuples(double prec, mcIdType limitTupleId, DataArrayIdType *&comm, DataArrayIdType *&commIndex) const
1046 std::size_t nbOfCompo=getNumberOfComponents();
1047 if ((nbOfCompo<1) || (nbOfCompo>4)) //test before work
1048 throw INTERP_KERNEL::Exception("DataArrayDouble::findCommonTuples : Unexpected spacedim of coords. Must be 1, 2, 3 or 4.");
1050 mcIdType nbOfTuples(getNumberOfTuples());
1052 MCAuto<DataArrayIdType> c(DataArrayIdType::New()),cI(DataArrayIdType::New()); c->alloc(0,1); cI->pushBackSilent(0);
1056 findCommonTuplesAlg<4>(begin(),nbOfTuples,limitTupleId,prec,c,cI);
1059 findCommonTuplesAlg<3>(begin(),nbOfTuples,limitTupleId,prec,c,cI);
1062 findCommonTuplesAlg<2>(begin(),nbOfTuples,limitTupleId,prec,c,cI);
1065 findCommonTuplesAlg<1>(begin(),nbOfTuples,limitTupleId,prec,c,cI);
1068 throw INTERP_KERNEL::Exception("DataArrayDouble::findCommonTuples : nb of components managed are 1,2,3 and 4 ! not implemented for other number of components !");
1071 commIndex=cI.retn();
1075 * This methods returns the minimal distance between the two set of points \a this and \a other.
1076 * So \a this and \a other have to have the same number of components. If not an INTERP_KERNEL::Exception will be thrown.
1077 * This method works only if number of components of \a this (equal to those of \a other) is in 1, 2 or 3.
1079 * \param [out] thisTupleId the tuple id in \a this corresponding to the returned minimal distance
1080 * \param [out] otherTupleId the tuple id in \a other corresponding to the returned minimal distance
1081 * \return the minimal distance between the two set of points \a this and \a other.
1082 * \sa DataArrayDouble::findClosestTupleId
1084 double DataArrayDouble::minimalDistanceTo(const DataArrayDouble *other, mcIdType& thisTupleId, mcIdType& otherTupleId) const
1086 MCAuto<DataArrayIdType> part1=findClosestTupleId(other);
1087 std::size_t nbOfCompo=getNumberOfComponents();
1088 mcIdType otherNbTuples=other->getNumberOfTuples();
1089 const double *thisPt(begin()),*otherPt(other->begin());
1090 const mcIdType *part1Pt(part1->begin());
1091 double ret=std::numeric_limits<double>::max();
1092 for(mcIdType i=0;i<otherNbTuples;i++,part1Pt++,otherPt+=nbOfCompo)
1095 for(std::size_t j=0;j<nbOfCompo;j++)
1096 tmp+=(otherPt[j]-thisPt[nbOfCompo*(*part1Pt)+j])*(otherPt[j]-thisPt[nbOfCompo*(*part1Pt)+j]);
1098 { ret=tmp; thisTupleId=*part1Pt; otherTupleId=i; }
1104 * This methods returns for each tuple in \a other which tuple in \a this is the closest.
1105 * So \a this and \a other have to have the same number of components. If not an INTERP_KERNEL::Exception will be thrown.
1106 * This method works only if number of components of \a this (equal to those of \a other) is in 1, 2 or 3.
1108 * \return a newly allocated (new object to be dealt by the caller) DataArrayInt having \c other->getNumberOfTuples() tuples and one components.
1109 * \sa DataArrayDouble::minimalDistanceTo
1111 DataArrayIdType *DataArrayDouble::findClosestTupleId(const DataArrayDouble *other) const
1114 throw INTERP_KERNEL::Exception("DataArrayDouble::findClosestTupleId : other instance is NULL !");
1115 checkAllocated(); other->checkAllocated();
1116 std::size_t nbOfCompo(getNumberOfComponents());
1117 if(nbOfCompo!=other->getNumberOfComponents())
1119 std::ostringstream oss; oss << "DataArrayDouble::findClosestTupleId : number of components in this is " << nbOfCompo;
1120 oss << ", whereas number of components in other is " << other->getNumberOfComponents() << "! Should be equal !";
1121 throw INTERP_KERNEL::Exception(oss.str().c_str());
1123 mcIdType nbOfTuples(other->getNumberOfTuples());
1124 mcIdType thisNbOfTuples(getNumberOfTuples());
1125 MCAuto<DataArrayIdType> ret=DataArrayIdType::New(); ret->alloc(nbOfTuples,1);
1127 getMinMaxPerComponent(bounds);
1132 double xDelta(fabs(bounds[1]-bounds[0])),yDelta(fabs(bounds[3]-bounds[2])),zDelta(fabs(bounds[5]-bounds[4]));
1133 double delta=std::max(xDelta,yDelta); delta=std::max(delta,zDelta);
1134 double characSize=pow((delta*delta*delta)/((double)thisNbOfTuples),1./3.);
1135 BBTreePts<3,mcIdType> myTree(begin(),0,0,getNumberOfTuples(),characSize*1e-12);
1136 FindClosestTupleIdAlg<3>(myTree,3.*characSize*characSize,other->begin(),nbOfTuples,begin(),thisNbOfTuples,ret->getPointer());
1141 double xDelta(fabs(bounds[1]-bounds[0])),yDelta(fabs(bounds[3]-bounds[2]));
1142 double delta=std::max(xDelta,yDelta);
1143 double characSize=sqrt(delta/(double)thisNbOfTuples);
1144 BBTreePts<2,mcIdType> myTree(begin(),0,0,getNumberOfTuples(),characSize*1e-12);
1145 FindClosestTupleIdAlg<2>(myTree,2.*characSize*characSize,other->begin(),nbOfTuples,begin(),thisNbOfTuples,ret->getPointer());
1150 double characSize=fabs(bounds[1]-bounds[0])/FromIdType<double>(thisNbOfTuples);
1151 BBTreePts<1,mcIdType> myTree(begin(),0,0,getNumberOfTuples(),characSize*1e-12);
1152 FindClosestTupleIdAlg<1>(myTree,1.*characSize*characSize,other->begin(),nbOfTuples,begin(),thisNbOfTuples,ret->getPointer());
1156 throw INTERP_KERNEL::Exception("Unexpected spacedim of coords for findClosestTupleId. Must be 1, 2 or 3.");
1162 * This method expects that \a this and \a otherBBoxFrmt arrays are bounding box arrays ( as the output of MEDCouplingPointSet::getBoundingBoxForBBTree method ).
1163 * 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
1164 * how many bounding boxes in \a otherBBoxFrmt.
1165 * So, this method expects that \a this and \a otherBBoxFrmt have the same number of components.
1167 * \param [in] otherBBoxFrmt - It is an array .
1168 * \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.
1169 * \sa MEDCouplingPointSet::getBoundingBoxForBBTree
1170 * \throw If \a this and \a otherBBoxFrmt have not the same number of components.
1171 * \throw If \a this and \a otherBBoxFrmt number of components is not even (BBox format).
1173 DataArrayIdType *DataArrayDouble::computeNbOfInteractionsWith(const DataArrayDouble *otherBBoxFrmt, double eps) const
1176 throw INTERP_KERNEL::Exception("DataArrayDouble::computeNbOfInteractionsWith : input array is NULL !");
1177 if(!isAllocated() || !otherBBoxFrmt->isAllocated())
1178 throw INTERP_KERNEL::Exception("DataArrayDouble::computeNbOfInteractionsWith : this and input array must be allocated !");
1179 std::size_t nbOfComp(getNumberOfComponents());
1180 mcIdType nbOfTuples(getNumberOfTuples());
1181 if(nbOfComp!=otherBBoxFrmt->getNumberOfComponents())
1183 std::ostringstream oss; oss << "DataArrayDouble::computeNbOfInteractionsWith : this number of components (" << nbOfComp << ") must be equal to the number of components of input array (" << otherBBoxFrmt->getNumberOfComponents() << ") !";
1184 throw INTERP_KERNEL::Exception(oss.str().c_str());
1188 std::ostringstream oss; oss << "DataArrayDouble::computeNbOfInteractionsWith : Number of components (" << nbOfComp << ") is not even ! It should be to be compatible with bbox format !";
1189 throw INTERP_KERNEL::Exception(oss.str().c_str());
1191 MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(nbOfTuples,1);
1192 const double *thisBBPtr(begin());
1193 mcIdType *retPtr(ret->getPointer());
1198 BBTree<3,mcIdType> bbt(otherBBoxFrmt->begin(),0,0,otherBBoxFrmt->getNumberOfTuples(),eps);
1199 for(mcIdType i=0;i<nbOfTuples;i++,retPtr++,thisBBPtr+=nbOfComp)
1200 *retPtr=bbt.getNbOfIntersectingElems(thisBBPtr);
1205 BBTree<2,mcIdType> bbt(otherBBoxFrmt->begin(),0,0,otherBBoxFrmt->getNumberOfTuples(),eps);
1206 for(mcIdType i=0;i<nbOfTuples;i++,retPtr++,thisBBPtr+=nbOfComp)
1207 *retPtr=bbt.getNbOfIntersectingElems(thisBBPtr);
1212 BBTree<1,mcIdType> bbt(otherBBoxFrmt->begin(),0,0,otherBBoxFrmt->getNumberOfTuples(),eps);
1213 for(mcIdType i=0;i<nbOfTuples;i++,retPtr++,thisBBPtr+=nbOfComp)
1214 *retPtr=bbt.getNbOfIntersectingElems(thisBBPtr);
1218 throw INTERP_KERNEL::Exception("DataArrayDouble::computeNbOfInteractionsWith : space dimension supported are [1,2,3] !");
1225 * Returns a copy of \a this array by excluding coincident tuples. Each tuple is
1226 * considered as coordinates of a point in getNumberOfComponents()-dimensional
1227 * space. The distance between tuples is computed using norm2. If several tuples are
1228 * not far each from other than \a prec, only one of them remains in the result
1229 * array. The order of tuples in the result array is same as in \a this one except
1230 * that coincident tuples are excluded.
1231 * \param [in] prec - minimal absolute distance between two tuples at which they are
1232 * considered not coincident.
1233 * \param [in] limitTupleId - limit tuple id. If all tuples within a group of coincident
1234 * tuples have id strictly lower than \a limitTupleId then they are not excluded.
1235 * \return DataArrayDouble * - the new instance of DataArrayDouble that the caller
1236 * is to delete using decrRef() as it is no more needed.
1237 * \throw If \a this is not allocated.
1238 * \throw If the number of components is not in [1,2,3,4].
1240 * \if ENABLE_EXAMPLES
1241 * \ref py_mcdataarraydouble_getdifferentvalues "Here is a Python example".
1244 DataArrayDouble *DataArrayDouble::getDifferentValues(double prec, mcIdType limitTupleId) const
1247 DataArrayIdType *c0=0,*cI0=0;
1248 findCommonTuples(prec,limitTupleId,c0,cI0);
1249 MCAuto<DataArrayIdType> c(c0),cI(cI0);
1250 mcIdType newNbOfTuples=-1;
1251 MCAuto<DataArrayIdType> o2n=DataArrayIdType::ConvertIndexArrayToO2N(getNumberOfTuples(),c0->begin(),cI0->begin(),cI0->end(),newNbOfTuples);
1252 return renumberAndReduce(o2n->getConstPointer(),newNbOfTuples);
1256 * Copy all components in a specified order from another DataArrayDouble.
1257 * Both numerical and textual data is copied. The number of tuples in \a this and
1258 * the other array can be different.
1259 * \param [in] a - the array to copy data from.
1260 * \param [in] compoIds - sequence of zero based indices of components, data of which is
1262 * \throw If \a a is NULL.
1263 * \throw If \a compoIds.size() != \a a->getNumberOfComponents().
1264 * \throw If \a compoIds[i] < 0 or \a compoIds[i] > \a this->getNumberOfComponents().
1266 * \if ENABLE_EXAMPLES
1267 * \ref py_mcdataarraydouble_setselectedcomponents "Here is a Python example".
1270 void DataArrayDouble::setSelectedComponents(const DataArrayDouble *a, const std::vector<std::size_t>& compoIds)
1273 throw INTERP_KERNEL::Exception("DataArrayDouble::setSelectedComponents : input DataArrayDouble is NULL !");
1275 copyPartOfStringInfoFrom2(compoIds,*a);
1276 std::size_t partOfCompoSz=compoIds.size();
1277 std::size_t nbOfCompo=getNumberOfComponents();
1278 mcIdType nbOfTuples=std::min(getNumberOfTuples(),a->getNumberOfTuples());
1279 const double *ac=a->getConstPointer();
1280 double *nc=getPointer();
1281 for(mcIdType i=0;i<nbOfTuples;i++)
1282 for(std::size_t j=0;j<partOfCompoSz;j++,ac++)
1283 nc[nbOfCompo*i+compoIds[j]]=*ac;
1286 * Checks if 0.0 value is present in \a this array. If it is the case, an exception
1288 * \throw If zero is found in \a this array.
1290 void DataArrayDouble::checkNoNullValues() const
1292 const double *tmp=getConstPointer();
1293 mcIdType nbOfElems=getNbOfElems();
1294 const double *where=std::find(tmp,tmp+nbOfElems,0.);
1295 if(where!=tmp+nbOfElems)
1296 throw INTERP_KERNEL::Exception("A value 0.0 have been detected !");
1300 * Computes minimal and maximal value in each component. An output array is filled
1301 * with \c 2 * \a this->getNumberOfComponents() values, so the caller is to allocate
1302 * enough memory before calling this method.
1303 * \param [out] bounds - array of size at least 2 *\a this->getNumberOfComponents().
1304 * It is filled as follows:<br>
1305 * \a bounds[0] = \c min_of_component_0 <br>
1306 * \a bounds[1] = \c max_of_component_0 <br>
1307 * \a bounds[2] = \c min_of_component_1 <br>
1308 * \a bounds[3] = \c max_of_component_1 <br>
1311 void DataArrayDouble::getMinMaxPerComponent(double *bounds) const
1314 std::size_t dim=getNumberOfComponents();
1315 for (std::size_t idim=0; idim<dim; idim++)
1317 bounds[idim*2]=std::numeric_limits<double>::max();
1318 bounds[idim*2+1]=-std::numeric_limits<double>::max();
1320 const double *ptr=getConstPointer();
1321 mcIdType nbOfTuples=getNumberOfTuples();
1322 for(mcIdType i=0;i<nbOfTuples;i++)
1324 for(std::size_t idim=0;idim<dim;idim++)
1326 if(bounds[idim*2]>ptr[i*dim+idim])
1328 bounds[idim*2]=ptr[i*dim+idim];
1330 if(bounds[idim*2+1]<ptr[i*dim+idim])
1332 bounds[idim*2+1]=ptr[i*dim+idim];
1339 * This method retrieves a newly allocated DataArrayDouble instance having same number of tuples than \a this and twice number of components than \a this
1340 * to store both the min and max per component of each tuples.
1341 * \param [in] epsilon the width of the bbox (identical in each direction) - 0.0 by default
1343 * \return a newly created DataArrayDouble instance having \c this->getNumberOfTuples() tuples and 2 * \c this->getNumberOfComponent() components
1345 * \throw If \a this is not allocated yet.
1347 DataArrayDouble *DataArrayDouble::computeBBoxPerTuple(double epsilon) const
1350 const double *dataPtr=getConstPointer();
1351 std::size_t nbOfCompo=getNumberOfComponents();
1352 mcIdType nbTuples=getNumberOfTuples();
1353 MCAuto<DataArrayDouble> bbox=DataArrayDouble::New();
1354 bbox->alloc(nbTuples,2*nbOfCompo);
1355 double *bboxPtr=bbox->getPointer();
1356 for(mcIdType i=0;i<nbTuples;i++)
1358 for(std::size_t j=0;j<nbOfCompo;j++)
1360 bboxPtr[2*nbOfCompo*i+2*j]=dataPtr[nbOfCompo*i+j]-epsilon;
1361 bboxPtr[2*nbOfCompo*i+2*j+1]=dataPtr[nbOfCompo*i+j]+epsilon;
1368 * For each tuples **t** in \a other, this method retrieves tuples in \a this that are equal to **t**.
1369 * Two tuples are considered equal if the euclidian distance between the two tuples is lower than \a eps.
1371 * \param [in] other a DataArrayDouble having same number of components than \a this.
1372 * \param [in] eps absolute precision representing distance (using infinite norm) between 2 tuples behind which 2 tuples are considered equal.
1373 * \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.
1374 * \a cI allows to extract information in \a c.
1375 * \param [out] cI is an indirection array that allows to extract the data contained in \a c.
1377 * \throw In case of:
1378 * - \a this is not allocated
1379 * - \a other is not allocated or null
1380 * - \a this and \a other do not have the same number of components
1381 * - if number of components of \a this is not in [1,2,3]
1383 * \sa MEDCouplingPointSet::getNodeIdsNearPoints, DataArrayDouble::getDifferentValues
1385 void DataArrayDouble::computeTupleIdsNearTuples(const DataArrayDouble *other, double eps, DataArrayIdType *& c, DataArrayIdType *& cI) const
1388 throw INTERP_KERNEL::Exception("DataArrayDouble::computeTupleIdsNearTuples : input pointer other is null !");
1390 other->checkAllocated();
1391 std::size_t nbOfCompo=getNumberOfComponents();
1392 std::size_t otherNbOfCompo=other->getNumberOfComponents();
1393 if(nbOfCompo!=otherNbOfCompo)
1394 throw INTERP_KERNEL::Exception("DataArrayDouble::computeTupleIdsNearTuples : number of components should be equal between this and other !");
1395 mcIdType nbOfTuplesOther=other->getNumberOfTuples();
1396 mcIdType nbOfTuples=getNumberOfTuples();
1397 MCAuto<DataArrayIdType> cArr(DataArrayIdType::New()),cIArr(DataArrayIdType::New()); cArr->alloc(0,1); cIArr->pushBackSilent(0);
1402 BBTreePts<3,mcIdType> myTree(begin(),0,0,nbOfTuples,eps);
1403 FindTupleIdsNearTuplesAlg<3>(myTree,other->getConstPointer(),nbOfTuplesOther,eps,cArr,cIArr);
1408 BBTreePts<2,mcIdType> myTree(begin(),0,0,nbOfTuples,eps);
1409 FindTupleIdsNearTuplesAlg<2>(myTree,other->getConstPointer(),nbOfTuplesOther,eps,cArr,cIArr);
1414 BBTreePts<1,mcIdType> myTree(begin(),0,0,nbOfTuples,eps);
1415 FindTupleIdsNearTuplesAlg<1>(myTree,other->getConstPointer(),nbOfTuplesOther,eps,cArr,cIArr);
1419 throw INTERP_KERNEL::Exception("Unexpected spacedim of coords for computeTupleIdsNearTuples. Must be 1, 2 or 3.");
1421 c=cArr.retn(); cI=cIArr.retn();
1425 * 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
1426 * around origin of 'radius' 1.
1428 * \param [in] eps absolute epsilon. under that value of delta between max and min no scale is performed.
1430 void DataArrayDouble::recenterForMaxPrecision(double eps)
1433 std::size_t dim=getNumberOfComponents();
1434 std::vector<double> bounds(2*dim);
1435 getMinMaxPerComponent(&bounds[0]);
1436 for(std::size_t i=0;i<dim;i++)
1438 double delta=bounds[2*i+1]-bounds[2*i];
1439 double offset=(bounds[2*i]+bounds[2*i+1])/2.;
1441 applyLin(1./delta,-offset/delta,i);
1443 applyLin(1.,-offset,i);
1448 * Returns the maximal value and all its locations within \a this one-dimensional array.
1449 * \param [out] tupleIds - a new instance of DataArrayInt containing indices of
1450 * tuples holding the maximal value. The caller is to delete it using
1451 * decrRef() as it is no more needed.
1452 * \return double - the maximal value among all values of \a this array.
1453 * \throw If \a this->getNumberOfComponents() != 1
1454 * \throw If \a this->getNumberOfTuples() < 1
1456 double DataArrayDouble::getMaxValue2(DataArrayIdType*& tupleIds) const
1460 double ret=getMaxValue(tmp);
1461 tupleIds=findIdsInRange(ret,ret);
1466 * Returns the minimal value and all its locations within \a this one-dimensional array.
1467 * \param [out] tupleIds - a new instance of DataArrayInt containing indices of
1468 * tuples holding the minimal value. The caller is to delete it using
1469 * decrRef() as it is no more needed.
1470 * \return double - the minimal value among all values of \a this array.
1471 * \throw If \a this->getNumberOfComponents() != 1
1472 * \throw If \a this->getNumberOfTuples() < 1
1474 double DataArrayDouble::getMinValue2(DataArrayIdType*& tupleIds) const
1478 double ret=getMinValue(tmp);
1479 tupleIds=findIdsInRange(ret,ret);
1484 * 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.
1485 * This method only works for single component array.
1487 * \return a value in [ 0, \c this->getNumberOfTuples() )
1489 * \throw If \a this is not allocated
1492 mcIdType DataArrayDouble::count(double value, double eps) const
1496 if(getNumberOfComponents()!=1)
1497 throw INTERP_KERNEL::Exception("DataArrayDouble::count : must be applied on DataArrayDouble with only one component, you can call 'rearrange' method before !");
1498 const double *vals=begin();
1499 mcIdType nbOfTuples=getNumberOfTuples();
1500 for(mcIdType i=0;i<nbOfTuples;i++,vals++)
1501 if(fabs(*vals-value)<=eps)
1507 * Returns the average value of \a this one-dimensional array.
1508 * \return double - the average value over all values of \a this array.
1509 * \throw If \a this->getNumberOfComponents() != 1
1510 * \throw If \a this->getNumberOfTuples() < 1
1512 double DataArrayDouble::getAverageValue() const
1514 if(getNumberOfComponents()!=1)
1515 throw INTERP_KERNEL::Exception("DataArrayDouble::getAverageValue : must be applied on DataArrayDouble with only one component, you can call 'rearrange' method before !");
1516 mcIdType nbOfTuples(getNumberOfTuples());
1518 throw INTERP_KERNEL::Exception("DataArrayDouble::getAverageValue : array exists but number of tuples must be > 0 !");
1519 const double *vals=getConstPointer();
1520 double ret=std::accumulate(vals,vals+nbOfTuples,0.);
1521 return ret/FromIdType<double>(nbOfTuples);
1525 * Returns the Euclidean norm of the vector defined by \a this array.
1526 * \return double - the value of the Euclidean norm, i.e.
1527 * the square root of the inner product of vector.
1528 * \throw If \a this is not allocated.
1530 double DataArrayDouble::norm2() const
1534 std::size_t nbOfElems=getNbOfElems();
1535 const double *pt=getConstPointer();
1536 for(std::size_t i=0;i<nbOfElems;i++,pt++)
1542 * Returns the maximum norm of the vector defined by \a this array.
1543 * This method works even if the number of components is different from one.
1544 * If the number of elements in \a this is 0, -1. is returned.
1545 * \return double - the value of the maximum norm, i.e.
1546 * the maximal absolute value among values of \a this array (whatever its number of components).
1547 * \throw If \a this is not allocated.
1549 double DataArrayDouble::normMax() const
1553 std::size_t nbOfElems(getNbOfElems());
1554 const double *pt(getConstPointer());
1555 for(std::size_t i=0;i<nbOfElems;i++,pt++)
1557 double val(std::abs(*pt));
1565 * Returns the maximum norm of for each component of \a this array.
1566 * If the number of elements in \a this is 0, -1. is returned.
1567 * \param [out] res - pointer to an array of result values, of size at least \a
1568 * this->getNumberOfComponents(), that is to be allocated by the caller.
1569 * \throw If \a this is not allocated.
1571 void DataArrayDouble::normMaxPerComponent(double * res) const
1574 mcIdType nbOfTuples(getNumberOfTuples());
1575 std::size_t nbOfCompos(getNumberOfComponents());
1576 std::fill(res, res+nbOfCompos, -1.0);
1577 const double *pt(getConstPointer());
1578 for(mcIdType i=0;i<nbOfTuples;i++)
1579 for (std::size_t j=0; j<nbOfCompos; j++, pt++)
1581 double val(std::abs(*pt));
1589 * Returns the minimum norm (absolute value) of the vector defined by \a this array.
1590 * This method works even if the number of components is different from one.
1591 * If the number of elements in \a this is 0, std::numeric_limits<double>::max() is returned.
1592 * \return double - the value of the minimum norm, i.e.
1593 * the minimal absolute value among values of \a this array (whatever its number of components).
1594 * \throw If \a this is not allocated.
1596 double DataArrayDouble::normMin() const
1599 double ret(std::numeric_limits<double>::max());
1600 std::size_t nbOfElems(getNbOfElems());
1601 const double *pt(getConstPointer());
1602 for(std::size_t i=0;i<nbOfElems;i++,pt++)
1604 double val(std::abs(*pt));
1612 * Accumulates values of each component of \a this array.
1613 * \param [out] res - an array of length \a this->getNumberOfComponents(), allocated
1614 * by the caller, that is filled by this method with sum value for each
1616 * \throw If \a this is not allocated.
1618 void DataArrayDouble::accumulate(double *res) const
1621 const double *ptr=getConstPointer();
1622 mcIdType nbTuple(getNumberOfTuples());
1623 std::size_t nbComps(getNumberOfComponents());
1624 std::fill(res,res+nbComps,0.);
1625 for(mcIdType i=0;i<nbTuple;i++)
1626 std::transform(ptr+i*nbComps,ptr+(i+1)*nbComps,res,res,std::plus<double>());
1630 * This method returns the min distance from an external tuple defined by [ \a tupleBg , \a tupleEnd ) to \a this and
1631 * the first tuple in \a this that matches the returned distance. If there is no tuples in \a this an exception will be thrown.
1634 * \a this is expected to be allocated and expected to have a number of components equal to the distance from \a tupleBg to
1635 * \a tupleEnd. If not an exception will be thrown.
1637 * \param [in] tupleBg start pointer (included) of input external tuple
1638 * \param [in] tupleEnd end pointer (not included) of input external tuple
1639 * \param [out] tupleId the tuple id in \a this that matches the min of distance between \a this and input external tuple
1640 * \return the min distance.
1641 * \sa MEDCouplingUMesh::distanceToPoint
1643 double DataArrayDouble::distanceToTuple(const double *tupleBg, const double *tupleEnd, mcIdType& tupleId) const
1646 mcIdType nbTuple(getNumberOfTuples());
1647 std::size_t nbComps(getNumberOfComponents());
1648 if(nbComps!=(std::size_t)std::distance(tupleBg,tupleEnd))
1649 { 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()); }
1651 throw INTERP_KERNEL::Exception("DataArrayDouble::distanceToTuple : no tuple in this ! No distance to compute !");
1652 double ret0=std::numeric_limits<double>::max();
1654 const double *work=getConstPointer();
1655 for(mcIdType i=0;i<nbTuple;i++)
1658 for(std::size_t j=0;j<nbComps;j++,work++)
1659 val+=(*work-tupleBg[j])*((*work-tupleBg[j]));
1663 { ret0=val; tupleId=i; }
1669 * Accumulate values of the given component of \a this array.
1670 * \param [in] compId - the index of the component of interest.
1671 * \return double - a sum value of \a compId-th component.
1672 * \throw If \a this is not allocated.
1673 * \throw If \a the condition ( 0 <= \a compId < \a this->getNumberOfComponents() ) is
1676 double DataArrayDouble::accumulate(std::size_t compId) const
1679 const double *ptr=getConstPointer();
1680 mcIdType nbTuple(getNumberOfTuples());
1681 std::size_t nbComps(getNumberOfComponents());
1683 throw INTERP_KERNEL::Exception("DataArrayDouble::accumulate : Invalid compId specified : No such nb of components !");
1685 for(mcIdType i=0;i<nbTuple;i++)
1686 ret+=ptr[i*nbComps+compId];
1691 * This method accumulate using addition tuples in \a this using input index array [ \a bgOfIndex, \a endOfIndex ).
1692 * The returned array will have same number of components than \a this and number of tuples equal to
1693 * \c std::distance(bgOfIndex,endOfIndex) \b minus \b one.
1695 * The input index array is expected to be ascendingly sorted in which the all referenced ids should be in [0, \c this->getNumberOfTuples).
1696 * 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.
1698 * \param [in] bgOfIndex - begin (included) of the input index array.
1699 * \param [in] endOfIndex - end (excluded) of the input index array.
1700 * \return DataArrayDouble * - the new instance having the same number of components than \a this.
1702 * \throw If bgOfIndex or end is NULL.
1703 * \throw If input index array is not ascendingly sorted.
1704 * \throw If there is an id in [ \a bgOfIndex, \a endOfIndex ) not in [0, \c this->getNumberOfTuples).
1705 * \throw If std::distance(bgOfIndex,endOfIndex)==0.
1707 DataArrayDouble *DataArrayDouble::accumulatePerChunck(const mcIdType *bgOfIndex, const mcIdType *endOfIndex) const
1709 if(!bgOfIndex || !endOfIndex)
1710 throw INTERP_KERNEL::Exception("DataArrayDouble::accumulatePerChunck : input pointer NULL !");
1712 std::size_t nbCompo(getNumberOfComponents());
1713 mcIdType nbOfTuples(getNumberOfTuples());
1714 std::size_t sz=std::distance(bgOfIndex,endOfIndex);
1716 throw INTERP_KERNEL::Exception("DataArrayDouble::accumulatePerChunck : invalid size of input index array !");
1718 MCAuto<DataArrayDouble> ret=DataArrayDouble::New(); ret->alloc(sz,nbCompo);
1719 const mcIdType *w=bgOfIndex;
1720 if(*w<0 || *w>=nbOfTuples)
1721 throw INTERP_KERNEL::Exception("DataArrayDouble::accumulatePerChunck : The first element of the input index not in [0,nbOfTuples) !");
1722 const double *srcPt=begin()+(*w)*nbCompo;
1723 double *tmp=ret->getPointer();
1724 for(std::size_t i=0;i<sz;i++,tmp+=nbCompo,w++)
1726 std::fill(tmp,tmp+nbCompo,0.);
1729 for(mcIdType j=w[0];j<w[1];j++,srcPt+=nbCompo)
1731 if(j>=0 && j<nbOfTuples)
1732 std::transform(srcPt,srcPt+nbCompo,tmp,tmp,std::plus<double>());
1735 std::ostringstream oss; oss << "DataArrayDouble::accumulatePerChunck : At rank #" << i << " the input index array points to id " << j << " should be in [0," << nbOfTuples << ") !";
1736 throw INTERP_KERNEL::Exception(oss.str().c_str());
1742 std::ostringstream oss; oss << "DataArrayDouble::accumulatePerChunck : At rank #" << i << " the input index array is not in ascendingly sorted.";
1743 throw INTERP_KERNEL::Exception(oss.str().c_str());
1746 ret->copyStringInfoFrom(*this);
1751 * 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.
1752 * This method expects that \a this as only one component. The returned array will have \a this->getNumberOfTuples()+1 tuple with also one component.
1753 * The ith element of returned array is equal to the sum of elements in \a this with rank strictly lower than i.
1755 * \return DataArrayDouble - A newly built array containing cum sum of \a this.
1757 MCAuto<DataArrayDouble> DataArrayDouble::cumSum() const
1760 checkNbOfComps(1,"DataArrayDouble::cumSum : this is expected to be single component");
1761 mcIdType nbOfTuple(getNumberOfTuples());
1762 MCAuto<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(nbOfTuple+1,1);
1763 double *ptr(ret->getPointer());
1765 const double *thisPtr(begin());
1766 for(mcIdType i=0;i<nbOfTuple;i++)
1767 ptr[i+1]=ptr[i]+thisPtr[i];
1772 * Converts each 2D point defined by the tuple of \a this array from the Polar to the
1773 * Cartesian coordinate system. The two components of the tuple of \a this array are
1774 * considered to contain (1) radius and (2) angle of the point in the Polar CS.
1775 * \return DataArrayDouble * - the new instance of DataArrayDouble, whose each tuple
1776 * contains X and Y coordinates of the point in the Cartesian CS. The caller
1777 * is to delete this array using decrRef() as it is no more needed. The array
1778 * does not contain any textual info on components.
1779 * \throw If \a this->getNumberOfComponents() != 2.
1780 * \sa fromCartToPolar
1782 DataArrayDouble *DataArrayDouble::fromPolarToCart() const
1785 std::size_t nbOfComp(getNumberOfComponents());
1787 throw INTERP_KERNEL::Exception("DataArrayDouble::fromPolarToCart : must be an array with exactly 2 components !");
1788 mcIdType nbOfTuple(getNumberOfTuples());
1789 DataArrayDouble *ret(DataArrayDouble::New());
1790 ret->alloc(nbOfTuple,2);
1791 double *w(ret->getPointer());
1792 const double *wIn(getConstPointer());
1793 for(mcIdType i=0;i<nbOfTuple;i++,w+=2,wIn+=2)
1795 w[0]=wIn[0]*cos(wIn[1]);
1796 w[1]=wIn[0]*sin(wIn[1]);
1802 * Converts each 3D point defined by the tuple of \a this array from the Cylindrical to
1803 * the Cartesian coordinate system. The three components of the tuple of \a this array
1804 * are considered to contain (1) radius, (2) azimuth and (3) altitude of the point in
1805 * the Cylindrical CS.
1806 * \return DataArrayDouble * - the new instance of DataArrayDouble, whose each tuple
1807 * contains X, Y and Z coordinates of the point in the Cartesian CS. The info
1808 * on the third component is copied from \a this array. The caller
1809 * is to delete this array using decrRef() as it is no more needed.
1810 * \throw If \a this->getNumberOfComponents() != 3.
1813 DataArrayDouble *DataArrayDouble::fromCylToCart() const
1816 std::size_t nbOfComp(getNumberOfComponents());
1818 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCylToCart : must be an array with exactly 3 components !");
1819 mcIdType nbOfTuple(getNumberOfTuples());
1820 DataArrayDouble *ret(DataArrayDouble::New());
1821 ret->alloc(getNumberOfTuples(),3);
1822 double *w(ret->getPointer());
1823 const double *wIn(getConstPointer());
1824 for(mcIdType i=0;i<nbOfTuple;i++,w+=3,wIn+=3)
1826 w[0]=wIn[0]*cos(wIn[1]);
1827 w[1]=wIn[0]*sin(wIn[1]);
1830 ret->setInfoOnComponent(2,getInfoOnComponent(2));
1835 * Converts each 3D point defined by the tuple of \a this array from the Spherical to
1836 * the Cartesian coordinate system. The three components of the tuple of \a this array
1837 * are considered to contain (1) radius, (2) polar angle and (3) azimuthal angle of the
1838 * point in the Cylindrical CS.
1839 * \return DataArrayDouble * - the new instance of DataArrayDouble, whose each tuple
1840 * contains X, Y and Z coordinates of the point in the Cartesian CS. The info
1841 * on the third component is copied from \a this array. The caller
1842 * is to delete this array using decrRef() as it is no more needed.
1843 * \throw If \a this->getNumberOfComponents() != 3.
1844 * \sa fromCartToSpher
1846 DataArrayDouble *DataArrayDouble::fromSpherToCart() const
1849 std::size_t nbOfComp(getNumberOfComponents());
1851 throw INTERP_KERNEL::Exception("DataArrayDouble::fromSpherToCart : must be an array with exactly 3 components !");
1852 mcIdType nbOfTuple(getNumberOfTuples());
1853 DataArrayDouble *ret(DataArrayDouble::New());
1854 ret->alloc(getNumberOfTuples(),3);
1855 double *w(ret->getPointer());
1856 const double *wIn(getConstPointer());
1857 for(mcIdType i=0;i<nbOfTuple;i++,w+=3,wIn+=3)
1859 w[0]=wIn[0]*cos(wIn[2])*sin(wIn[1]);
1860 w[1]=wIn[0]*sin(wIn[2])*sin(wIn[1]);
1861 w[2]=wIn[0]*cos(wIn[1]);
1867 * 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.
1868 * 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.
1869 * If \a at equals to AX_CYL the returned array will be the result of operation cylindric to cartesian of \a this...
1871 * \param [in] atOfThis - The axis type of \a this.
1872 * \return DataArrayDouble * - the new instance of DataArrayDouble (that must be dealed by caller) containing the result of the cartesianizification of \a this.
1874 DataArrayDouble *DataArrayDouble::cartesianize(MEDCouplingAxisType atOfThis) const
1877 std::size_t nbOfComp(getNumberOfComponents());
1878 MCAuto<DataArrayDouble> ret;
1887 ret=fromCylToCart();
1892 ret=fromPolarToCart();
1896 throw INTERP_KERNEL::Exception("DataArrayDouble::cartesianize : For AX_CYL, number of components must be in [2,3] !");
1900 ret=fromSpherToCart();
1905 ret=fromPolarToCart();
1909 throw INTERP_KERNEL::Exception("DataArrayDouble::cartesianize : For AX_CYL, number of components must be in [2,3] !");
1911 throw INTERP_KERNEL::Exception("DataArrayDouble::cartesianize : not recognized axis type ! Only AX_CART, AX_CYL and AX_SPHER supported !");
1913 ret->copyStringInfoFrom(*this);
1918 * This method returns a newly created array to be deallocated that contains the result of conversion from cartesian to polar.
1919 * This method expects that \a this has exactly 2 components.
1920 * \sa fromPolarToCart
1922 DataArrayDouble *DataArrayDouble::fromCartToPolar() const
1924 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
1926 std::size_t nbOfComp(getNumberOfComponents());
1927 mcIdType nbTuples(getNumberOfTuples());
1929 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToPolar : must be an array with exactly 2 components !");
1930 ret->alloc(nbTuples,2);
1931 double *retPtr(ret->getPointer());
1932 const double *ptr(begin());
1933 for(mcIdType i=0;i<nbTuples;i++,ptr+=2,retPtr+=2)
1935 retPtr[0]=sqrt(ptr[0]*ptr[0]+ptr[1]*ptr[1]);
1936 retPtr[1]=atan2(ptr[1],ptr[0]);
1942 * This method returns a newly created array to be deallocated that contains the result of conversion from cartesian to cylindrical.
1943 * This method expects that \a this has exactly 3 components.
1946 DataArrayDouble *DataArrayDouble::fromCartToCyl() const
1948 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
1950 std::size_t nbOfComp(getNumberOfComponents());
1951 mcIdType nbTuples(getNumberOfTuples());
1953 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCyl : must be an array with exactly 3 components !");
1954 ret->alloc(nbTuples,3);
1955 double *retPtr(ret->getPointer());
1956 const double *ptr(begin());
1957 for(mcIdType i=0;i<nbTuples;i++,ptr+=3,retPtr+=3)
1959 retPtr[0]=sqrt(ptr[0]*ptr[0]+ptr[1]*ptr[1]);
1960 retPtr[1]=atan2(ptr[1],ptr[0]);
1967 * This method returns a newly created array to be deallocated that contains the result of conversion from cartesian to spherical coordinates.
1968 * \sa fromSpherToCart
1970 DataArrayDouble *DataArrayDouble::fromCartToSpher() const
1972 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
1974 std::size_t nbOfComp(getNumberOfComponents());
1975 mcIdType nbTuples(getNumberOfTuples());
1977 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToSpher : must be an array with exactly 3 components !");
1978 ret->alloc(nbTuples,3);
1979 double *retPtr(ret->getPointer());
1980 const double *ptr(begin());
1981 for(mcIdType i=0;i<nbTuples;i++,ptr+=3,retPtr+=3)
1983 retPtr[0]=sqrt(ptr[0]*ptr[0]+ptr[1]*ptr[1]+ptr[2]*ptr[2]);
1984 retPtr[1]=acos(ptr[2]/retPtr[0]);
1985 retPtr[2]=atan2(ptr[1],ptr[0]);
1991 * 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.
1992 * This method expects that \a this has exactly 3 components.
1993 * \sa MEDCouplingFieldDouble::computeVectorFieldCyl
1995 DataArrayDouble *DataArrayDouble::fromCartToCylGiven(const DataArrayDouble *coords, const double center[3], const double vect[3]) const
1998 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCylGiven : input coords are NULL !");
1999 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
2000 checkAllocated(); coords->checkAllocated();
2001 std::size_t nbOfComp(getNumberOfComponents());
2002 mcIdType nbTuples(getNumberOfTuples());
2004 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCylGiven : must be an array with exactly 3 components !");
2005 if(coords->getNumberOfComponents()!=3)
2006 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCylGiven : coords array must have exactly 3 components !");
2007 if(coords->getNumberOfTuples()!=nbTuples)
2008 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCylGiven : coords array must have the same number of tuples !");
2009 ret->alloc(nbTuples,nbOfComp);
2010 double magOfVect(sqrt(vect[0]*vect[0]+vect[1]*vect[1]+vect[2]*vect[2]));
2012 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCylGiven : magnitude of vect is too low !");
2013 double Ur[3],Uteta[3],Uz[3],*retPtr(ret->getPointer());
2014 const double *coo(coords->begin()),*vectField(begin());
2015 std::transform(vect,vect+3,Uz,std::bind(std::multiplies<double>(),std::placeholders::_1,1./magOfVect));
2016 for(mcIdType i=0;i<nbTuples;i++,vectField+=3,retPtr+=3,coo+=3)
2018 std::transform(coo,coo+3,center,Ur,std::minus<double>());
2019 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];
2020 double magOfTeta(sqrt(Uteta[0]*Uteta[0]+Uteta[1]*Uteta[1]+Uteta[2]*Uteta[2]));
2021 std::transform(Uteta,Uteta+3,Uteta,std::bind(std::multiplies<double>(),std::placeholders::_1,1./magOfTeta));
2022 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];
2023 retPtr[0]=Ur[0]*vectField[0]+Ur[1]*vectField[1]+Ur[2]*vectField[2];
2024 retPtr[1]=Uteta[0]*vectField[0]+Uteta[1]*vectField[1]+Uteta[2]*vectField[2];
2025 retPtr[2]=Uz[0]*vectField[0]+Uz[1]*vectField[1]+Uz[2]*vectField[2];
2027 ret->copyStringInfoFrom(*this);
2032 * Computes the doubly contracted product of every tensor defined by the tuple of \a this
2033 * array containing 6 components.
2034 * \return DataArrayDouble * - the new instance of DataArrayDouble, whose each tuple
2035 * is calculated from the tuple <em>(t)</em> of \a this array as follows:
2036 * \f$ t[0]^2+t[1]^2+t[2]^2+2*t[3]^2+2*t[4]^2+2*t[5]^2\f$.
2037 * The caller is to delete this result array using decrRef() as it is no more needed.
2038 * \throw If \a this->getNumberOfComponents() != 6.
2040 DataArrayDouble *DataArrayDouble::doublyContractedProduct() const
2043 std::size_t nbOfComp(getNumberOfComponents());
2045 throw INTERP_KERNEL::Exception("DataArrayDouble::doublyContractedProduct : must be an array with exactly 6 components !");
2046 DataArrayDouble *ret=DataArrayDouble::New();
2047 mcIdType nbOfTuple=getNumberOfTuples();
2048 ret->alloc(nbOfTuple,1);
2049 const double *src=getConstPointer();
2050 double *dest=ret->getPointer();
2051 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=6)
2052 *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];
2057 * Computes the determinant of every square matrix defined by the tuple of \a this
2058 * array, which contains either 4, 6 or 9 components. The case of 6 components
2059 * corresponds to that of the upper triangular matrix.
2060 * \return DataArrayDouble * - the new instance of DataArrayDouble, whose each tuple
2061 * is the determinant of matrix of the corresponding tuple of \a this array.
2062 * The caller is to delete this result array using decrRef() as it is no more
2064 * \throw If \a this->getNumberOfComponents() is not in [4,6,9].
2066 DataArrayDouble *DataArrayDouble::determinant() const
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 switch(getNumberOfComponents())
2077 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=6)
2078 *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];
2081 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=4)
2082 *dest=src[0]*src[3]-src[1]*src[2];
2085 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=9)
2086 *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];
2090 throw INTERP_KERNEL::Exception("DataArrayDouble::determinant : Invalid number of components ! must be in 4,6,9 !");
2095 * Computes 3 eigenvalues of every upper triangular matrix defined by the tuple of
2096 * \a this array, which contains 6 components. The 6 components of tuples are expected to be stored as follow :<br>
2097 * \a tuple[0] = \c matrix_XX <br>
2098 * \a tuple[1] = \c matrix_YY <br>
2099 * \a tuple[2] = \c matrix_ZZ <br>
2100 * \a tuple[3] = \c matrix_XY <br>
2101 * \a tuple[4] = \c matrix_YZ <br>
2102 * \a tuple[5] = \c matrix_XZ <br>
2103 * \return DataArrayDouble * - the new instance of DataArrayDouble containing 3
2104 * components, whose each tuple contains the eigenvalues of the matrix of
2105 * corresponding tuple of \a this array.
2106 * The caller is to delete this result array using decrRef() as it is no more
2108 * \throw If \a this->getNumberOfComponents() != 6.
2110 DataArrayDouble *DataArrayDouble::eigenValues() const
2113 std::size_t nbOfComp=getNumberOfComponents();
2115 throw INTERP_KERNEL::Exception("DataArrayDouble::eigenValues : must be an array with exactly 6 components !");
2116 DataArrayDouble *ret=DataArrayDouble::New();
2117 mcIdType nbOfTuple=getNumberOfTuples();
2118 ret->alloc(nbOfTuple,3);
2119 const double *src=getConstPointer();
2120 double *dest=ret->getPointer();
2121 for(mcIdType i=0;i<nbOfTuple;i++,dest+=3,src+=6)
2122 INTERP_KERNEL::computeEigenValues6(src,dest);
2127 * Computes 3 eigenvectors of every upper triangular matrix defined by the tuple of
2128 * \a this array, which contains 6 components. The 6 components of tuples are expected to be stored as follow :<br>
2129 * \a tuple[0] = \c matrix_XX <br>
2130 * \a tuple[1] = \c matrix_YY <br>
2131 * \a tuple[2] = \c matrix_ZZ <br>
2132 * \a tuple[3] = \c matrix_XY <br>
2133 * \a tuple[4] = \c matrix_YZ <br>
2134 * \a tuple[5] = \c matrix_XZ <br>
2135 * \return DataArrayDouble * - the new instance of DataArrayDouble containing 9
2136 * components, whose each tuple contains 3 eigenvectors of the matrix of
2137 * corresponding tuple of \a this array.
2138 * The caller is to delete this result array using decrRef() as it is no more
2140 * \throw If \a this->getNumberOfComponents() != 6.
2142 DataArrayDouble *DataArrayDouble::eigenVectors() const
2145 std::size_t nbOfComp=getNumberOfComponents();
2147 throw INTERP_KERNEL::Exception("DataArrayDouble::eigenVectors : must be an array with exactly 6 components !");
2148 DataArrayDouble *ret=DataArrayDouble::New();
2149 mcIdType nbOfTuple=getNumberOfTuples();
2150 ret->alloc(nbOfTuple,9);
2151 const double *src=getConstPointer();
2152 double *dest=ret->getPointer();
2153 for(mcIdType i=0;i<nbOfTuple;i++,src+=6)
2156 INTERP_KERNEL::computeEigenValues6(src,tmp);
2157 for(mcIdType j=0;j<3;j++,dest+=3)
2158 INTERP_KERNEL::computeEigenVectorForEigenValue6(src,tmp[j],1e-12,dest);
2164 * Computes the inverse matrix of every matrix defined by the tuple of \a this
2165 * array, which contains either 4, 6 or 9 components. The case of 6 components
2166 * corresponds to that of the upper triangular matrix.
2167 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2168 * same number of components as \a this one, whose each tuple is the inverse
2169 * matrix of the matrix of corresponding tuple of \a this array.
2170 * The caller is to delete this result array using decrRef() as it is no more
2172 * \throw If \a this->getNumberOfComponents() is not in [4,6,9].
2174 DataArrayDouble *DataArrayDouble::inverse() const
2177 std::size_t nbOfComp=getNumberOfComponents();
2178 if(nbOfComp!=6 && nbOfComp!=9 && nbOfComp!=4)
2179 throw INTERP_KERNEL::Exception("DataArrayDouble::inversion : must be an array with 4,6 or 9 components !");
2180 DataArrayDouble *ret=DataArrayDouble::New();
2181 mcIdType nbOfTuple=getNumberOfTuples();
2182 ret->alloc(nbOfTuple,nbOfComp);
2183 const double *src=getConstPointer();
2184 double *dest=ret->getPointer();
2186 for(mcIdType i=0;i<nbOfTuple;i++,dest+=6,src+=6)
2188 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];
2189 dest[0]=(src[1]*src[2]-src[4]*src[4])/det;
2190 dest[1]=(src[0]*src[2]-src[5]*src[5])/det;
2191 dest[2]=(src[0]*src[1]-src[3]*src[3])/det;
2192 dest[3]=(src[5]*src[4]-src[3]*src[2])/det;
2193 dest[4]=(src[5]*src[3]-src[0]*src[4])/det;
2194 dest[5]=(src[3]*src[4]-src[1]*src[5])/det;
2196 else if(nbOfComp==4)
2197 for(mcIdType i=0;i<nbOfTuple;i++,dest+=4,src+=4)
2199 double det=src[0]*src[3]-src[1]*src[2];
2201 dest[1]=-src[1]/det;
2202 dest[2]=-src[2]/det;
2206 for(mcIdType i=0;i<nbOfTuple;i++,dest+=9,src+=9)
2208 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];
2209 dest[0]=(src[4]*src[8]-src[7]*src[5])/det;
2210 dest[1]=(src[7]*src[2]-src[1]*src[8])/det;
2211 dest[2]=(src[1]*src[5]-src[4]*src[2])/det;
2212 dest[3]=(src[6]*src[5]-src[3]*src[8])/det;
2213 dest[4]=(src[0]*src[8]-src[6]*src[2])/det;
2214 dest[5]=(src[2]*src[3]-src[0]*src[5])/det;
2215 dest[6]=(src[3]*src[7]-src[6]*src[4])/det;
2216 dest[7]=(src[6]*src[1]-src[0]*src[7])/det;
2217 dest[8]=(src[0]*src[4]-src[1]*src[3])/det;
2223 * Computes the trace of every matrix defined by the tuple of \a this
2224 * array, which contains either 4, 6 or 9 components. The case of 6 components
2225 * corresponds to that of the upper triangular matrix.
2226 * \return DataArrayDouble * - the new instance of DataArrayDouble containing
2227 * 1 component, whose each tuple is the trace of
2228 * the matrix of corresponding tuple of \a this array.
2229 * The caller is to delete this result array using decrRef() as it is no more
2231 * \throw If \a this->getNumberOfComponents() is not in [4,6,9].
2233 DataArrayDouble *DataArrayDouble::trace() const
2236 std::size_t nbOfComp=getNumberOfComponents();
2237 if(nbOfComp!=6 && nbOfComp!=9 && nbOfComp!=4)
2238 throw INTERP_KERNEL::Exception("DataArrayDouble::trace : must be an array with 4,6 or 9 components !");
2239 DataArrayDouble *ret=DataArrayDouble::New();
2240 mcIdType nbOfTuple=getNumberOfTuples();
2241 ret->alloc(nbOfTuple,1);
2242 const double *src=getConstPointer();
2243 double *dest=ret->getPointer();
2245 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=6)
2246 *dest=src[0]+src[1]+src[2];
2247 else if(nbOfComp==4)
2248 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=4)
2249 *dest=src[0]+src[3];
2251 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=9)
2252 *dest=src[0]+src[4]+src[8];
2257 * Computes the stress deviator tensor of every stress tensor defined by the tuple of
2258 * \a this array, which contains 6 components.
2259 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2260 * same number of components and tuples as \a this array.
2261 * The caller is to delete this result array using decrRef() as it is no more
2263 * \throw If \a this->getNumberOfComponents() != 6.
2265 DataArrayDouble *DataArrayDouble::deviator() const
2268 std::size_t nbOfComp=getNumberOfComponents();
2270 throw INTERP_KERNEL::Exception("DataArrayDouble::deviator : must be an array with exactly 6 components !");
2271 DataArrayDouble *ret=DataArrayDouble::New();
2272 mcIdType nbOfTuple=getNumberOfTuples();
2273 ret->alloc(nbOfTuple,6);
2274 const double *src=getConstPointer();
2275 double *dest=ret->getPointer();
2276 for(mcIdType i=0;i<nbOfTuple;i++,dest+=6,src+=6)
2278 double tr=(src[0]+src[1]+src[2])/3.;
2290 * Computes the magnitude of every vector defined by the tuple of
2292 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2293 * same number of tuples as \a this array and one component.
2294 * The caller is to delete this result array using decrRef() as it is no more
2296 * \throw If \a this is not allocated.
2298 DataArrayDouble *DataArrayDouble::magnitude() const
2301 std::size_t nbOfComp=getNumberOfComponents();
2302 DataArrayDouble *ret=DataArrayDouble::New();
2303 mcIdType nbOfTuple=getNumberOfTuples();
2304 ret->alloc(nbOfTuple,1);
2305 const double *src=getConstPointer();
2306 double *dest=ret->getPointer();
2307 for(mcIdType i=0;i<nbOfTuple;i++,dest++)
2310 for(std::size_t j=0;j<nbOfComp;j++,src++)
2318 * Computes the maximal value within every tuple of \a this array.
2319 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2320 * same number of tuples as \a this array and one component.
2321 * The caller is to delete this result array using decrRef() as it is no more
2323 * \throw If \a this is not allocated.
2324 * \sa DataArrayDouble::maxPerTupleWithCompoId
2326 DataArrayDouble *DataArrayDouble::maxPerTuple() const
2329 std::size_t nbOfComp(getNumberOfComponents());
2330 MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
2331 mcIdType nbOfTuple(getNumberOfTuples());
2332 ret->alloc(nbOfTuple,1);
2333 const double *src=getConstPointer();
2334 double *dest=ret->getPointer();
2335 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=nbOfComp)
2336 *dest=*std::max_element(src,src+nbOfComp);
2341 * Computes the maximal value within every tuple of \a this array and it returns the first component
2342 * id for each tuple that corresponds to the maximal value within the tuple.
2344 * \param [out] compoIdOfMaxPerTuple - the new new instance of DataArrayInt containing the
2345 * same number of tuples and only one component.
2346 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2347 * same number of tuples as \a this array and one component.
2348 * The caller is to delete this result array using decrRef() as it is no more
2350 * \throw If \a this is not allocated.
2351 * \sa DataArrayDouble::maxPerTuple
2353 DataArrayDouble *DataArrayDouble::maxPerTupleWithCompoId(DataArrayIdType* &compoIdOfMaxPerTuple) const
2356 std::size_t nbOfComp(getNumberOfComponents());
2357 MCAuto<DataArrayDouble> ret0=DataArrayDouble::New();
2358 MCAuto<DataArrayIdType> ret1=DataArrayIdType::New();
2359 mcIdType nbOfTuple=getNumberOfTuples();
2360 ret0->alloc(nbOfTuple,1); ret1->alloc(nbOfTuple,1);
2361 const double *src=getConstPointer();
2362 double *dest=ret0->getPointer(); mcIdType *dest1=ret1->getPointer();
2363 for(mcIdType i=0;i<nbOfTuple;i++,dest++,dest1++,src+=nbOfComp)
2365 const double *loc=std::max_element(src,src+nbOfComp);
2367 *dest1=ToIdType(std::distance(src,loc));
2369 compoIdOfMaxPerTuple=ret1.retn();
2374 * This method returns a newly allocated DataArrayDouble instance having one component and \c this->getNumberOfTuples() * \c this->getNumberOfTuples() tuples.
2375 * \n This returned array contains the euclidian distance for each tuple in \a this.
2376 * \n So the returned array can be seen as a dense symmetrical matrix whose diagonal elements are equal to 0.
2377 * \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)
2379 * \warning use this method with care because it can leads to big amount of consumed memory !
2381 * \return A newly allocated (huge) MEDCoupling::DataArrayDouble instance that the caller should deal with.
2383 * \throw If \a this is not allocated.
2385 * \sa DataArrayDouble::buildEuclidianDistanceDenseMatrixWith
2387 DataArrayDouble *DataArrayDouble::buildEuclidianDistanceDenseMatrix() const
2390 std::size_t nbOfComp(getNumberOfComponents());
2391 mcIdType nbOfTuples(getNumberOfTuples());
2392 const double *inData=getConstPointer();
2393 MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
2394 ret->alloc(nbOfTuples*nbOfTuples,1);
2395 double *outData=ret->getPointer();
2396 for(mcIdType i=0;i<nbOfTuples;i++)
2398 outData[i*nbOfTuples+i]=0.;
2399 for(mcIdType j=i+1;j<nbOfTuples;j++)
2402 for(std::size_t k=0;k<nbOfComp;k++)
2403 { double delta=inData[i*nbOfComp+k]-inData[j*nbOfComp+k]; dist+=delta*delta; }
2405 outData[i*nbOfTuples+j]=dist;
2406 outData[j*nbOfTuples+i]=dist;
2413 * This method returns a newly allocated DataArrayDouble instance having one component and \c this->getNumberOfTuples() * \c other->getNumberOfTuples() tuples.
2414 * \n This returned array contains the euclidian distance for each tuple in \a other with each tuple in \a this.
2415 * \n So the returned array can be seen as a dense rectangular matrix with \c other->getNumberOfTuples() rows and \c this->getNumberOfTuples() columns.
2416 * \n Output rectangular matrix is sorted along rows.
2417 * \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)
2419 * \warning use this method with care because it can leads to big amount of consumed memory !
2421 * \param [in] other DataArrayDouble instance having same number of components than \a this.
2422 * \return A newly allocated (huge) MEDCoupling::DataArrayDouble instance that the caller should deal with.
2424 * \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.
2426 * \sa DataArrayDouble::buildEuclidianDistanceDenseMatrix
2428 DataArrayDouble *DataArrayDouble::buildEuclidianDistanceDenseMatrixWith(const DataArrayDouble *other) const
2431 throw INTERP_KERNEL::Exception("DataArrayDouble::buildEuclidianDistanceDenseMatrixWith : input parameter is null !");
2433 other->checkAllocated();
2434 std::size_t nbOfComp(getNumberOfComponents());
2435 std::size_t otherNbOfComp(other->getNumberOfComponents());
2436 if(nbOfComp!=otherNbOfComp)
2438 std::ostringstream oss; oss << "DataArrayDouble::buildEuclidianDistanceDenseMatrixWith : this nb of compo=" << nbOfComp << " and other nb of compo=" << otherNbOfComp << ". It should match !";
2439 throw INTERP_KERNEL::Exception(oss.str().c_str());
2441 mcIdType nbOfTuples(getNumberOfTuples());
2442 mcIdType otherNbOfTuples(other->getNumberOfTuples());
2443 const double *inData=getConstPointer();
2444 const double *inDataOther=other->getConstPointer();
2445 MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
2446 ret->alloc(otherNbOfTuples*nbOfTuples,1);
2447 double *outData=ret->getPointer();
2448 for(mcIdType i=0;i<otherNbOfTuples;i++,inDataOther+=nbOfComp)
2450 for(mcIdType j=0;j<nbOfTuples;j++)
2453 for(std::size_t k=0;k<nbOfComp;k++)
2454 { double delta=inDataOther[k]-inData[j*nbOfComp+k]; dist+=delta*delta; }
2456 outData[i*nbOfTuples+j]=dist;
2463 * This method expects that \a this stores 3 tuples containing 2 components each.
2464 * Each of this tuples represent a point into 2D space.
2465 * 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).
2466 * 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[.
2468 * \throw If \a this is not allocated.
2469 * \throw If \a this has not 3 tuples of 2 components
2470 * \throw If tuples/points in \a this are aligned
2472 void DataArrayDouble::asArcOfCircle(double center[2], double& radius, double& ang) const
2475 INTERP_KERNEL::QuadraticPlanarPrecision arcPrec(1e-14);
2476 if(getNumberOfTuples()!=3 && getNumberOfComponents()!=2)
2477 throw INTERP_KERNEL::Exception("DataArrayDouble::asArcCircle : this method expects");
2478 const double *pt(begin());
2479 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]));
2481 INTERP_KERNEL::AutoCppPtr<INTERP_KERNEL::EdgeLin> e1(new INTERP_KERNEL::EdgeLin(n0,n2)),e2(new INTERP_KERNEL::EdgeLin(n2,n1));
2482 INTERP_KERNEL::SegSegIntersector inters(*e1,*e2);
2483 bool colinearity(inters.areColinears());
2485 throw INTERP_KERNEL::Exception("DataArrayDouble::asArcOfCircle : 3 points in this have been detected as colinear !");
2487 INTERP_KERNEL::AutoCppPtr<INTERP_KERNEL::EdgeArcCircle> ret(new INTERP_KERNEL::EdgeArcCircle(n0,n2,n1));
2488 const double *c(ret->getCenter());
2489 center[0]=c[0]; center[1]=c[1];
2490 radius=ret->getRadius();
2491 ang=ret->getAngle();
2495 * Sorts value within every tuple of \a this array.
2496 * \param [in] asc - if \a true, the values are sorted in ascending order, else,
2497 * in descending order.
2498 * \throw If \a this is not allocated.
2500 void DataArrayDouble::sortPerTuple(bool asc)
2503 double *pt=getPointer();
2504 mcIdType nbOfTuple(getNumberOfTuples());
2505 std::size_t nbOfComp(getNumberOfComponents());
2507 for(mcIdType i=0;i<nbOfTuple;i++,pt+=nbOfComp)
2508 std::sort(pt,pt+nbOfComp);
2510 for(mcIdType i=0;i<nbOfTuple;i++,pt+=nbOfComp)
2511 std::sort(pt,pt+nbOfComp,std::greater<double>());
2516 * Modify all elements of \a this array, so that
2517 * an element _x_ becomes \f$ numerator / x \f$.
2518 * \warning If an exception is thrown because of presence of 0.0 element in \a this
2519 * array, all elements processed before detection of the zero element remain
2521 * \param [in] numerator - the numerator used to modify array elements.
2522 * \throw If \a this is not allocated.
2523 * \throw If there is an element equal to 0.0 in \a this array.
2525 void DataArrayDouble::applyInv(double numerator)
2528 double *ptr=getPointer();
2529 std::size_t nbOfElems=getNbOfElems();
2530 for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2532 if(std::abs(*ptr)>std::numeric_limits<double>::min())
2534 *ptr=numerator/(*ptr);
2538 std::ostringstream oss; oss << "DataArrayDouble::applyInv : presence of null value in tuple #" << i/getNumberOfComponents() << " component #" << i%getNumberOfComponents();
2540 throw INTERP_KERNEL::Exception(oss.str().c_str());
2547 * Modify all elements of \a this array, so that
2548 * an element _x_ becomes <em> val ^ x </em>. Contrary to DataArrayInt::applyPow
2549 * all values in \a this have to be >= 0 if val is \b not integer.
2550 * \param [in] val - the value used to apply pow on all array elements.
2551 * \throw If \a this is not allocated.
2552 * \warning If an exception is thrown because of presence of 0 element in \a this
2553 * array and \a val is \b not integer, all elements processed before detection of the zero element remain
2556 void DataArrayDouble::applyPow(double val)
2559 double *ptr=getPointer();
2560 std::size_t nbOfElems=getNbOfElems();
2562 bool isInt=((double)val2)==val;
2565 for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2571 std::ostringstream oss; oss << "DataArrayDouble::applyPow (double) : At elem # " << i << " value is " << *ptr << " ! must be >=0. !";
2572 throw INTERP_KERNEL::Exception(oss.str().c_str());
2578 for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2579 *ptr=pow(*ptr,val2);
2585 * Modify all elements of \a this array, so that
2586 * an element _x_ becomes \f$ val ^ x \f$.
2587 * \param [in] val - the value used to apply pow on all array elements.
2588 * \throw If \a this is not allocated.
2589 * \throw If \a val < 0.
2590 * \warning If an exception is thrown because of presence of 0 element in \a this
2591 * array, all elements processed before detection of the zero element remain
2594 void DataArrayDouble::applyRPow(double val)
2598 throw INTERP_KERNEL::Exception("DataArrayDouble::applyRPow : the input value has to be >= 0 !");
2599 double *ptr=getPointer();
2600 std::size_t nbOfElems=getNbOfElems();
2601 for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2607 * Returns a new DataArrayDouble created from \a this one by applying \a
2608 * FunctionToEvaluate to every tuple of \a this array. Textual data is not copied.
2609 * For more info see \ref MEDCouplingArrayApplyFunc
2610 * \param [in] nbOfComp - number of components in the result array.
2611 * \param [in] func - the \a FunctionToEvaluate declared as
2612 * \c bool (*\a func)(\c const \c double *\a pos, \c double *\a res),
2613 * where \a pos points to the first component of a tuple of \a this array
2614 * and \a res points to the first component of a tuple of the result array.
2615 * Note that length (number of components) of \a pos can differ from
2617 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2618 * same number of tuples as \a this array.
2619 * The caller is to delete this result array using decrRef() as it is no more
2621 * \throw If \a this is not allocated.
2622 * \throw If \a func returns \a false.
2624 DataArrayDouble *DataArrayDouble::applyFunc(std::size_t nbOfComp, FunctionToEvaluate func) const
2627 DataArrayDouble *newArr=DataArrayDouble::New();
2628 mcIdType nbOfTuples(getNumberOfTuples());
2629 std::size_t oldNbOfComp(getNumberOfComponents());
2630 newArr->alloc(nbOfTuples,nbOfComp);
2631 const double *ptr=getConstPointer();
2632 double *ptrToFill=newArr->getPointer();
2633 for(mcIdType i=0;i<nbOfTuples;i++)
2635 if(!func(ptr+i*oldNbOfComp,ptrToFill+i*nbOfComp))
2637 std::ostringstream oss; oss << "For tuple # " << i << " with value (";
2638 std::copy(ptr+oldNbOfComp*i,ptr+oldNbOfComp*(i+1),std::ostream_iterator<double>(oss,", "));
2639 oss << ") : Evaluation of function failed !";
2641 throw INTERP_KERNEL::Exception(oss.str().c_str());
2648 * Returns a new DataArrayDouble created from \a this one by applying a function to every
2649 * tuple of \a this array. Textual data is not copied.
2650 * For more info see \ref MEDCouplingArrayApplyFunc1.
2651 * \param [in] nbOfComp - number of components in the result array.
2652 * \param [in] func - the expression defining how to transform a tuple of \a this array.
2653 * Supported expressions are described \ref MEDCouplingArrayApplyFuncExpr "here".
2654 * \param [in] isSafe - By default true. If true invalid operation (division by 0. acos of value > 1. ...) leads to a throw of an exception.
2655 * If false the computation is carried on without any notification. When false the evaluation is a little faster.
2656 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2657 * same number of tuples as \a this array and \a nbOfComp components.
2658 * The caller is to delete this result array using decrRef() as it is no more
2660 * \throw If \a this is not allocated.
2661 * \throw If computing \a func fails.
2663 DataArrayDouble *DataArrayDouble::applyFunc(std::size_t nbOfComp, const std::string& func, bool isSafe) const
2665 INTERP_KERNEL::ExprParser expr(func);
2667 std::set<std::string> vars;
2668 expr.getTrueSetOfVars(vars);
2669 std::vector<std::string> varsV(vars.begin(),vars.end());
2670 return applyFuncNamedCompo(nbOfComp,varsV,func,isSafe);
2674 * Returns a new DataArrayDouble created from \a this one by applying a function to every
2675 * tuple of \a this array. Textual data is not copied. This method works by tuples (whatever its size).
2676 * If \a this is a one component array, call applyFuncOnThis instead that performs the same work faster.
2678 * For more info see \ref MEDCouplingArrayApplyFunc0.
2679 * \param [in] func - the expression defining how to transform a tuple of \a this array.
2680 * Supported expressions are described \ref MEDCouplingArrayApplyFuncExpr "here".
2681 * \param [in] isSafe - By default true. If true invalid operation (division by 0. acos of value > 1. ...) leads to a throw of an exception.
2682 * If false the computation is carried on without any notification. When false the evaluation is a little faster.
2683 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2684 * same number of tuples and components as \a this array.
2685 * The caller is to delete this result array using decrRef() as it is no more
2687 * \sa applyFuncOnThis
2688 * \throw If \a this is not allocated.
2689 * \throw If computing \a func fails.
2691 DataArrayDouble *DataArrayDouble::applyFunc(const std::string& func, bool isSafe) const
2693 std::size_t nbOfComp(getNumberOfComponents());
2695 throw INTERP_KERNEL::Exception("DataArrayDouble::applyFunc : output number of component must be > 0 !");
2697 mcIdType nbOfTuples(getNumberOfTuples());
2698 MCAuto<DataArrayDouble> newArr(DataArrayDouble::New());
2699 newArr->alloc(nbOfTuples,nbOfComp);
2700 INTERP_KERNEL::ExprParser expr(func);
2702 std::set<std::string> vars;
2703 expr.getTrueSetOfVars(vars);
2706 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 : ";
2707 std::copy(vars.begin(),vars.end(),std::ostream_iterator<std::string>(oss," "));
2708 throw INTERP_KERNEL::Exception(oss.str().c_str());
2712 expr.prepareFastEvaluator();
2713 newArr->rearrange(1);
2714 newArr->fillWithValue(expr.evaluateDouble());
2715 newArr->rearrange(nbOfComp);
2716 return newArr.retn();
2718 std::vector<std::string> vars2(vars.begin(),vars.end());
2719 double buff,*ptrToFill(newArr->getPointer());
2720 const double *ptr(begin());
2721 std::vector<double> stck;
2722 expr.prepareExprEvaluationDouble(vars2,1,1,0,&buff,&buff+1);
2723 expr.prepareFastEvaluator();
2726 for(mcIdType i=0;i<nbOfTuples;i++)
2728 for(std::size_t iComp=0;iComp<nbOfComp;iComp++,ptr++,ptrToFill++)
2731 expr.evaluateDoubleInternal(stck);
2732 *ptrToFill=stck.back();
2739 for(mcIdType i=0;i<nbOfTuples;i++)
2741 for(std::size_t iComp=0;iComp<nbOfComp;iComp++,ptr++,ptrToFill++)
2746 expr.evaluateDoubleInternalSafe(stck);
2748 catch(INTERP_KERNEL::Exception& e)
2750 std::ostringstream oss; oss << "For tuple # " << i << " component # " << iComp << " with value (";
2752 oss << ") : Evaluation of function failed !" << e.what();
2753 throw INTERP_KERNEL::Exception(oss.str().c_str());
2755 *ptrToFill=stck.back();
2760 return newArr.retn();
2764 * This method is a non const method that modify the array in \a this.
2765 * This method only works on one component array. It means that function \a func must
2766 * contain at most one variable.
2767 * This method is a specialization of applyFunc method with one parameter on one component array.
2769 * \param [in] func - the expression defining how to transform a tuple of \a this array.
2770 * Supported expressions are described \ref MEDCouplingArrayApplyFuncExpr "here".
2771 * \param [in] isSafe - By default true. If true invalid operation (division by 0. acos of value > 1. ...) leads to a throw of an exception.
2772 * If false the computation is carried on without any notification. When false the evaluation is a little faster.
2776 void DataArrayDouble::applyFuncOnThis(const std::string& func, bool isSafe)
2778 std::size_t nbOfComp(getNumberOfComponents());
2780 throw INTERP_KERNEL::Exception("DataArrayDouble::applyFuncOnThis : output number of component must be > 0 !");
2782 mcIdType nbOfTuples(getNumberOfTuples());
2783 INTERP_KERNEL::ExprParser expr(func);
2785 std::set<std::string> vars;
2786 expr.getTrueSetOfVars(vars);
2789 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 : ";
2790 std::copy(vars.begin(),vars.end(),std::ostream_iterator<std::string>(oss," "));
2791 throw INTERP_KERNEL::Exception(oss.str().c_str());
2795 expr.prepareFastEvaluator();
2796 std::vector<std::string> compInfo(getInfoOnComponents());
2798 fillWithValue(expr.evaluateDouble());
2799 rearrange(nbOfComp);
2800 setInfoOnComponents(compInfo);
2803 std::vector<std::string> vars2(vars.begin(),vars.end());
2804 double buff,*ptrToFill(getPointer());
2805 const double *ptr(begin());
2806 std::vector<double> stck;
2807 expr.prepareExprEvaluationDouble(vars2,1,1,0,&buff,&buff+1);
2808 expr.prepareFastEvaluator();
2811 for(mcIdType i=0;i<nbOfTuples;i++)
2813 for(std::size_t iComp=0;iComp<nbOfComp;iComp++,ptr++,ptrToFill++)
2816 expr.evaluateDoubleInternal(stck);
2817 *ptrToFill=stck.back();
2824 for(mcIdType i=0;i<nbOfTuples;i++)
2826 for(std::size_t iComp=0;iComp<nbOfComp;iComp++,ptr++,ptrToFill++)
2831 expr.evaluateDoubleInternalSafe(stck);
2833 catch(INTERP_KERNEL::Exception& e)
2835 std::ostringstream oss; oss << "For tuple # " << i << " component # " << iComp << " with value (";
2837 oss << ") : Evaluation of function failed !" << e.what();
2838 throw INTERP_KERNEL::Exception(oss.str().c_str());
2840 *ptrToFill=stck.back();
2848 * Returns a new DataArrayDouble created from \a this one by applying a function to every
2849 * tuple of \a this array. Textual data is not copied.
2850 * For more info see \ref MEDCouplingArrayApplyFunc2.
2851 * \param [in] nbOfComp - number of components in the result array.
2852 * \param [in] func - the expression defining how to transform a tuple of \a this array.
2853 * Supported expressions are described \ref MEDCouplingArrayApplyFuncExpr "here".
2854 * \param [in] isSafe - By default true. If true invalid operation (division by 0. acos of value > 1. ...) leads to a throw of an exception.
2855 * If false the computation is carried on without any notification. When false the evaluation is a little faster.
2856 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2857 * same number of tuples as \a this array.
2858 * The caller is to delete this result array using decrRef() as it is no more
2860 * \throw If \a this is not allocated.
2861 * \throw If \a func contains vars that are not in \a this->getInfoOnComponent().
2862 * \throw If computing \a func fails.
2864 DataArrayDouble *DataArrayDouble::applyFuncCompo(std::size_t nbOfComp, const std::string& func, bool isSafe) const
2866 return applyFuncNamedCompo(nbOfComp,getVarsOnComponent(),func,isSafe);
2870 * Returns a new DataArrayDouble created from \a this one by applying a function to every
2871 * tuple of \a this array. Textual data is not copied.
2872 * For more info see \ref MEDCouplingArrayApplyFunc3.
2873 * \param [in] nbOfComp - number of components in the result array.
2874 * \param [in] varsOrder - sequence of vars defining their order.
2875 * \param [in] func - the expression defining how to transform a tuple of \a this array.
2876 * Supported expressions are described \ref MEDCouplingArrayApplyFuncExpr "here".
2877 * \param [in] isSafe - By default true. If true invalid operation (division by 0. acos of value > 1. ...) leads to a throw of an exception.
2878 * If false the computation is carried on without any notification. When false the evaluation is a little faster.
2879 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2880 * same number of tuples as \a this array.
2881 * The caller is to delete this result array using decrRef() as it is no more
2883 * \throw If \a this is not allocated.
2884 * \throw If \a func contains vars not in \a varsOrder.
2885 * \throw If computing \a func fails.
2887 DataArrayDouble *DataArrayDouble::applyFuncNamedCompo(std::size_t nbOfComp, const std::vector<std::string>& varsOrder, const std::string& func, bool isSafe) const
2890 throw INTERP_KERNEL::Exception("DataArrayDouble::applyFuncNamedCompo : output number of component must be > 0 !");
2891 std::vector<std::string> varsOrder2(varsOrder);
2892 std::size_t oldNbOfComp(getNumberOfComponents());
2893 for(std::size_t i=varsOrder.size();i<oldNbOfComp;i++)
2894 varsOrder2.push_back(std::string());
2896 mcIdType nbOfTuples(getNumberOfTuples());
2897 INTERP_KERNEL::ExprParser expr(func);
2899 std::set<std::string> vars;
2900 expr.getTrueSetOfVars(vars);
2901 if(vars.size()>oldNbOfComp)
2903 std::ostringstream oss; oss << "The field has " << oldNbOfComp << " components and there are ";
2904 oss << vars.size() << " variables : ";
2905 std::copy(vars.begin(),vars.end(),std::ostream_iterator<std::string>(oss," "));
2906 throw INTERP_KERNEL::Exception(oss.str().c_str());
2908 MCAuto<DataArrayDouble> newArr(DataArrayDouble::New());
2909 newArr->alloc(nbOfTuples,nbOfComp);
2910 INTERP_KERNEL::AutoPtr<double> buff(new double[oldNbOfComp]);
2911 double *buffPtr(buff),*ptrToFill;
2912 std::vector<double> stck;
2913 for(std::size_t iComp=0;iComp<nbOfComp;iComp++)
2915 expr.prepareExprEvaluationDouble(varsOrder2,(int)oldNbOfComp,(int)nbOfComp,(int)iComp,buffPtr,buffPtr+oldNbOfComp);
2916 expr.prepareFastEvaluator();
2917 const double *ptr(getConstPointer());
2918 ptrToFill=newArr->getPointer()+iComp;
2921 for(mcIdType i=0;i<nbOfTuples;i++,ptrToFill+=nbOfComp,ptr+=oldNbOfComp)
2923 std::copy(ptr,ptr+oldNbOfComp,buffPtr);
2924 expr.evaluateDoubleInternal(stck);
2925 *ptrToFill=stck.back();
2931 for(mcIdType i=0;i<nbOfTuples;i++,ptrToFill+=nbOfComp,ptr+=oldNbOfComp)
2933 std::copy(ptr,ptr+oldNbOfComp,buffPtr);
2936 expr.evaluateDoubleInternalSafe(stck);
2937 *ptrToFill=stck.back();
2940 catch(INTERP_KERNEL::Exception& e)
2942 std::ostringstream oss; oss << "For tuple # " << i << " with value (";
2943 std::copy(ptr+oldNbOfComp*i,ptr+oldNbOfComp*(i+1),std::ostream_iterator<double>(oss,", "));
2944 oss << ") : Evaluation of function failed !" << e.what();
2945 throw INTERP_KERNEL::Exception(oss.str().c_str());
2950 return newArr.retn();
2953 void DataArrayDouble::applyFuncFast32(const std::string& func)
2956 INTERP_KERNEL::ExprParser expr(func);
2958 char *funcStr=expr.compileX86();
2960 *((void **)&funcPtr)=funcStr;//he he...
2962 double *ptr=getPointer();
2963 std::size_t nbOfComp=getNumberOfComponents();
2964 mcIdType nbOfTuples=getNumberOfTuples();
2965 std::size_t nbOfElems=nbOfTuples*nbOfComp;
2966 for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2971 void DataArrayDouble::applyFuncFast64(const std::string& func)
2974 INTERP_KERNEL::ExprParser expr(func);
2976 char *funcStr=expr.compileX86_64();
2978 *((void **)&funcPtr)=funcStr;//he he...
2980 double *ptr=getPointer();
2981 std::size_t nbOfComp=getNumberOfComponents();
2982 mcIdType nbOfTuples=getNumberOfTuples();
2983 std::size_t nbOfElems=nbOfTuples*nbOfComp;
2984 for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2990 * \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.
2992 MCAuto<DataArrayDouble> DataArrayDouble::symmetry3DPlane(const double point[3], const double normalVector[3]) const
2995 if(getNumberOfComponents()!=3)
2996 throw INTERP_KERNEL::Exception("DataArrayDouble::symmetry3DPlane : this is excepted to have 3 components !");
2997 mcIdType nbTuples(getNumberOfTuples());
2998 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
2999 ret->alloc(nbTuples,3);
3000 Symmetry3DPlane(point,normalVector,nbTuples,begin(),ret->getPointer());
3004 DataArrayDoubleIterator *DataArrayDouble::iterator()
3006 return new DataArrayDoubleIterator(this);
3010 * Returns a new DataArrayInt containing indices of tuples of \a this one-dimensional
3011 * array whose values are within a given range. Textual data is not copied.
3012 * \param [in] vmin - a lowest acceptable value (included).
3013 * \param [in] vmax - a greatest acceptable value (included).
3014 * \return DataArrayInt * - the new instance of DataArrayInt.
3015 * The caller is to delete this result array using decrRef() as it is no more
3017 * \throw If \a this->getNumberOfComponents() != 1.
3019 * \sa DataArrayDouble::findIdsNotInRange
3021 * \if ENABLE_EXAMPLES
3022 * \ref cpp_mcdataarraydouble_getidsinrange "Here is a C++ example".<br>
3023 * \ref py_mcdataarraydouble_getidsinrange "Here is a Python example".
3026 DataArrayIdType *DataArrayDouble::findIdsInRange(double vmin, double vmax) const
3029 if(getNumberOfComponents()!=1)
3030 throw INTERP_KERNEL::Exception("DataArrayDouble::findIdsInRange : this must have exactly one component !");
3031 const double *cptr(begin());
3032 MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(0,1);
3033 mcIdType nbOfTuples(getNumberOfTuples());
3034 for(mcIdType i=0;i<nbOfTuples;i++,cptr++)
3035 if(*cptr>=vmin && *cptr<=vmax)
3036 ret->pushBackSilent(i);
3041 * Returns a new DataArrayInt containing indices of tuples of \a this one-dimensional
3042 * array whose values are not within a given range. Textual data is not copied.
3043 * \param [in] vmin - a lowest not acceptable value (excluded).
3044 * \param [in] vmax - a greatest not acceptable value (excluded).
3045 * \return DataArrayInt * - the new instance of DataArrayInt.
3046 * The caller is to delete this result array using decrRef() as it is no more
3048 * \throw If \a this->getNumberOfComponents() != 1.
3050 * \sa DataArrayDouble::findIdsInRange
3052 DataArrayIdType *DataArrayDouble::findIdsNotInRange(double vmin, double vmax) const
3055 if(getNumberOfComponents()!=1)
3056 throw INTERP_KERNEL::Exception("DataArrayDouble::findIdsNotInRange : this must have exactly one component !");
3057 const double *cptr(begin());
3058 MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(0,1);
3059 mcIdType nbOfTuples(getNumberOfTuples());
3060 for(mcIdType i=0;i<nbOfTuples;i++,cptr++)
3061 if(*cptr<vmin || *cptr>vmax)
3062 ret->pushBackSilent(i);
3067 * Returns a new DataArrayDouble by concatenating two given arrays, so that (1) the number
3068 * of tuples in the result array is a sum of the number of tuples of given arrays and (2)
3069 * the number of component in the result array is same as that of each of given arrays.
3070 * Info on components is copied from the first of the given arrays. Number of components
3071 * in the given arrays must be the same.
3072 * \param [in] a1 - an array to include in the result array.
3073 * \param [in] a2 - another array to include in the result array.
3074 * \return DataArrayDouble * - the new instance of DataArrayDouble.
3075 * The caller is to delete this result array using decrRef() as it is no more
3077 * \throw If both \a a1 and \a a2 are NULL.
3078 * \throw If \a a1->getNumberOfComponents() != \a a2->getNumberOfComponents().
3080 DataArrayDouble *DataArrayDouble::Aggregate(const DataArrayDouble *a1, const DataArrayDouble *a2)
3082 std::vector<const DataArrayDouble *> tmp(2);
3083 tmp[0]=a1; tmp[1]=a2;
3084 return Aggregate(tmp);
3088 * Returns a new DataArrayDouble by concatenating all given arrays, so that (1) the number
3089 * of tuples in the result array is a sum of the number of tuples of given arrays and (2)
3090 * the number of component in the result array is same as that of each of given arrays.
3091 * Info on components is copied from the first of the given arrays. Number of components
3092 * in the given arrays must be the same.
3093 * If the number of non null of elements in \a arr is equal to one the returned object is a copy of it
3094 * not the object itself.
3095 * \param [in] arr - a sequence of arrays to include in the result array.
3096 * \return DataArrayDouble * - the new instance of DataArrayDouble.
3097 * The caller is to delete this result array using decrRef() as it is no more
3099 * \throw If all arrays within \a arr are NULL.
3100 * \throw If getNumberOfComponents() of arrays within \a arr.
3102 DataArrayDouble *DataArrayDouble::Aggregate(const std::vector<const DataArrayDouble *>& arr)
3104 std::vector<const DataArrayDouble *> a;
3105 for(std::vector<const DataArrayDouble *>::const_iterator it4=arr.begin();it4!=arr.end();it4++)
3109 throw INTERP_KERNEL::Exception("DataArrayDouble::Aggregate : input list must contain at least one NON EMPTY DataArrayDouble !");
3110 std::vector<const DataArrayDouble *>::const_iterator it=a.begin();
3111 std::size_t nbOfComp((*it)->getNumberOfComponents());
3112 mcIdType nbt=(*it++)->getNumberOfTuples();
3113 for(mcIdType i=1;it!=a.end();it++,i++)
3115 if((*it)->getNumberOfComponents()!=nbOfComp)
3116 throw INTERP_KERNEL::Exception("DataArrayDouble::Aggregate : Nb of components mismatch for array aggregation !");
3117 nbt+=(*it)->getNumberOfTuples();
3119 MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
3120 ret->alloc(nbt,nbOfComp);
3121 double *pt=ret->getPointer();
3122 for(it=a.begin();it!=a.end();it++)
3123 pt=std::copy((*it)->getConstPointer(),(*it)->getConstPointer()+(*it)->getNbOfElems(),pt);
3124 ret->copyStringInfoFrom(*(a[0]));
3129 * Returns a new DataArrayDouble containing a dot product of two given arrays, so that
3130 * the i-th tuple of the result array is a sum of products of j-th components of i-th
3131 * tuples of given arrays (\f$ a_i = \sum_{j=1}^n a1_j * a2_j \f$).
3132 * Info on components and name is copied from the first of the given arrays.
3133 * Number of tuples and components in the given arrays must be the same.
3134 * \param [in] a1 - a given array.
3135 * \param [in] a2 - another given array.
3136 * \return DataArrayDouble * - the new instance of DataArrayDouble.
3137 * The caller is to delete this result array using decrRef() as it is no more
3139 * \throw If either \a a1 or \a a2 is NULL.
3140 * \throw If any given array is not allocated.
3141 * \throw If \a a1->getNumberOfTuples() != \a a2->getNumberOfTuples()
3142 * \throw If \a a1->getNumberOfComponents() != \a a2->getNumberOfComponents()
3144 DataArrayDouble *DataArrayDouble::Dot(const DataArrayDouble *a1, const DataArrayDouble *a2)
3147 throw INTERP_KERNEL::Exception("DataArrayDouble::Dot : input DataArrayDouble instance is NULL !");
3148 a1->checkAllocated();
3149 a2->checkAllocated();
3150 std::size_t nbOfComp(a1->getNumberOfComponents());
3151 if(nbOfComp!=a2->getNumberOfComponents())
3152 throw INTERP_KERNEL::Exception("Nb of components mismatch for array Dot !");
3153 mcIdType nbOfTuple(a1->getNumberOfTuples());
3154 if(nbOfTuple!=a2->getNumberOfTuples())
3155 throw INTERP_KERNEL::Exception("Nb of tuples mismatch for array Dot !");
3156 DataArrayDouble *ret=DataArrayDouble::New();
3157 ret->alloc(nbOfTuple,1);
3158 double *retPtr=ret->getPointer();
3159 const double *a1Ptr=a1->begin(),*a2Ptr(a2->begin());
3160 for(mcIdType i=0;i<nbOfTuple;i++)
3163 for(std::size_t j=0;j<nbOfComp;j++)
3164 sum+=a1Ptr[i*nbOfComp+j]*a2Ptr[i*nbOfComp+j];
3167 ret->setInfoOnComponent(0,a1->getInfoOnComponent(0));
3168 ret->setName(a1->getName());
3173 * Returns a new DataArrayDouble containing a cross product of two given arrays, so that
3174 * the i-th tuple of the result array contains 3 components of a vector which is a cross
3175 * product of two vectors defined by the i-th tuples of given arrays.
3176 * Info on components is copied from the first of the given arrays.
3177 * Number of tuples in the given arrays must be the same.
3178 * Number of components in the given arrays must be 3.
3179 * \param [in] a1 - a given array.
3180 * \param [in] a2 - another given array.
3181 * \return DataArrayDouble * - the new instance of DataArrayDouble.
3182 * The caller is to delete this result array using decrRef() as it is no more
3184 * \throw If either \a a1 or \a a2 is NULL.
3185 * \throw If \a a1->getNumberOfTuples() != \a a2->getNumberOfTuples()
3186 * \throw If \a a1->getNumberOfComponents() != 3
3187 * \throw If \a a2->getNumberOfComponents() != 3
3189 DataArrayDouble *DataArrayDouble::CrossProduct(const DataArrayDouble *a1, const DataArrayDouble *a2)
3192 throw INTERP_KERNEL::Exception("DataArrayDouble::CrossProduct : input DataArrayDouble instance is NULL !");
3193 std::size_t nbOfComp(a1->getNumberOfComponents());
3194 if(nbOfComp!=a2->getNumberOfComponents())
3195 throw INTERP_KERNEL::Exception("Nb of components mismatch for array crossProduct !");
3197 throw INTERP_KERNEL::Exception("Nb of components must be equal to 3 for array crossProduct !");
3198 mcIdType nbOfTuple(a1->getNumberOfTuples());
3199 if(nbOfTuple!=a2->getNumberOfTuples())
3200 throw INTERP_KERNEL::Exception("Nb of tuples mismatch for array crossProduct !");
3201 DataArrayDouble *ret=DataArrayDouble::New();
3202 ret->alloc(nbOfTuple,3);
3203 double *retPtr=ret->getPointer();
3204 const double *a1Ptr(a1->begin()),*a2Ptr(a2->begin());
3205 for(mcIdType i=0;i<nbOfTuple;i++)
3207 retPtr[3*i]=a1Ptr[3*i+1]*a2Ptr[3*i+2]-a1Ptr[3*i+2]*a2Ptr[3*i+1];
3208 retPtr[3*i+1]=a1Ptr[3*i+2]*a2Ptr[3*i]-a1Ptr[3*i]*a2Ptr[3*i+2];
3209 retPtr[3*i+2]=a1Ptr[3*i]*a2Ptr[3*i+1]-a1Ptr[3*i+1]*a2Ptr[3*i];
3211 ret->copyStringInfoFrom(*a1);
3216 * Returns a new DataArrayDouble containing maximal values of two given arrays.
3217 * Info on components is copied from the first of the given arrays.
3218 * Number of tuples and components in the given arrays must be the same.
3219 * \param [in] a1 - an array to compare values with another one.
3220 * \param [in] a2 - another array to compare values with the first one.
3221 * \return DataArrayDouble * - the new instance of DataArrayDouble.
3222 * The caller is to delete this result array using decrRef() as it is no more
3224 * \throw If either \a a1 or \a a2 is NULL.
3225 * \throw If \a a1->getNumberOfTuples() != \a a2->getNumberOfTuples()
3226 * \throw If \a a1->getNumberOfComponents() != \a a2->getNumberOfComponents()
3228 DataArrayDouble *DataArrayDouble::Max(const DataArrayDouble *a1, const DataArrayDouble *a2)
3231 throw INTERP_KERNEL::Exception("DataArrayDouble::Max : input DataArrayDouble instance is NULL !");
3232 std::size_t nbOfComp(a1->getNumberOfComponents());
3233 if(nbOfComp!=a2->getNumberOfComponents())
3234 throw INTERP_KERNEL::Exception("Nb of components mismatch for array Max !");
3235 mcIdType nbOfTuple(a1->getNumberOfTuples());
3236 if(nbOfTuple!=a2->getNumberOfTuples())
3237 throw INTERP_KERNEL::Exception("Nb of tuples mismatch for array Max !");
3238 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
3239 ret->alloc(nbOfTuple,nbOfComp);
3240 double *retPtr(ret->getPointer());
3241 const double *a1Ptr(a1->begin()),*a2Ptr(a2->begin());
3242 std::size_t nbElem(nbOfTuple*nbOfComp);
3243 for(std::size_t i=0;i<nbElem;i++)
3244 retPtr[i]=std::max(a1Ptr[i],a2Ptr[i]);
3245 ret->copyStringInfoFrom(*a1);
3250 * Returns a new DataArrayDouble containing minimal values of two given arrays.
3251 * Info on components is copied from the first of the given arrays.
3252 * Number of tuples and components in the given arrays must be the same.
3253 * \param [in] a1 - an array to compare values with another one.
3254 * \param [in] a2 - another array to compare values with the first one.
3255 * \return DataArrayDouble * - the new instance of DataArrayDouble.
3256 * The caller is to delete this result array using decrRef() as it is no more
3258 * \throw If either \a a1 or \a a2 is NULL.
3259 * \throw If \a a1->getNumberOfTuples() != \a a2->getNumberOfTuples()
3260 * \throw If \a a1->getNumberOfComponents() != \a a2->getNumberOfComponents()
3262 DataArrayDouble *DataArrayDouble::Min(const DataArrayDouble *a1, const DataArrayDouble *a2)
3265 throw INTERP_KERNEL::Exception("DataArrayDouble::Min : input DataArrayDouble instance is NULL !");
3266 std::size_t nbOfComp(a1->getNumberOfComponents());
3267 if(nbOfComp!=a2->getNumberOfComponents())
3268 throw INTERP_KERNEL::Exception("Nb of components mismatch for array min !");
3269 mcIdType nbOfTuple(a1->getNumberOfTuples());
3270 if(nbOfTuple!=a2->getNumberOfTuples())
3271 throw INTERP_KERNEL::Exception("Nb of tuples mismatch for array min !");
3272 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
3273 ret->alloc(nbOfTuple,nbOfComp);
3274 double *retPtr(ret->getPointer());
3275 const double *a1Ptr(a1->begin()),*a2Ptr(a2->begin());
3276 std::size_t nbElem(nbOfTuple*nbOfComp);
3277 for(std::size_t i=0;i<nbElem;i++)
3278 retPtr[i]=std::min(a1Ptr[i],a2Ptr[i]);
3279 ret->copyStringInfoFrom(*a1);
3284 * Returns a new DataArrayDouble that is the result of pow of two given arrays. There are 3
3287 * \param [in] a1 - an array to pow up.
3288 * \param [in] a2 - another array to sum up.
3289 * \return DataArrayDouble * - the new instance of DataArrayDouble.
3290 * The caller is to delete this result array using decrRef() as it is no more
3292 * \throw If either \a a1 or \a a2 is NULL.
3293 * \throw If \a a1->getNumberOfTuples() != \a a2->getNumberOfTuples()
3294 * \throw If \a a1->getNumberOfComponents() != 1 or \a a2->getNumberOfComponents() != 1.
3295 * \throw If there is a negative value in \a a1.
3297 DataArrayDouble *DataArrayDouble::Pow(const DataArrayDouble *a1, const DataArrayDouble *a2)
3300 throw INTERP_KERNEL::Exception("DataArrayDouble::Pow : at least one of input instances is null !");
3301 mcIdType nbOfTuple=a1->getNumberOfTuples();
3302 mcIdType nbOfTuple2=a2->getNumberOfTuples();
3303 std::size_t nbOfComp=a1->getNumberOfComponents();
3304 std::size_t nbOfComp2=a2->getNumberOfComponents();
3305 if(nbOfTuple!=nbOfTuple2)
3306 throw INTERP_KERNEL::Exception("DataArrayDouble::Pow : number of tuples mismatches !");
3307 if(nbOfComp!=1 || nbOfComp2!=1)
3308 throw INTERP_KERNEL::Exception("DataArrayDouble::Pow : number of components of both arrays must be equal to 1 !");
3309 MCAuto<DataArrayDouble> ret=DataArrayDouble::New(); ret->alloc(nbOfTuple,1);
3310 const double *ptr1(a1->begin()),*ptr2(a2->begin());
3311 double *ptr=ret->getPointer();
3312 for(mcIdType i=0;i<nbOfTuple;i++,ptr1++,ptr2++,ptr++)
3316 *ptr=pow(*ptr1,*ptr2);
3320 std::ostringstream oss; oss << "DataArrayDouble::Pow : on tuple #" << i << " of a1 value is < 0 (" << *ptr1 << ") !";
3321 throw INTERP_KERNEL::Exception(oss.str().c_str());
3328 * Apply pow on values of another DataArrayDouble to values of \a this one.
3330 * \param [in] other - an array to pow to \a this one.
3331 * \throw If \a other is NULL.
3332 * \throw If \a this->getNumberOfTuples() != \a other->getNumberOfTuples()
3333 * \throw If \a this->getNumberOfComponents() != 1 or \a other->getNumberOfComponents() != 1
3334 * \throw If there is a negative value in \a this.
3336 void DataArrayDouble::powEqual(const DataArrayDouble *other)
3339 throw INTERP_KERNEL::Exception("DataArrayDouble::powEqual : input instance is null !");
3340 mcIdType nbOfTuple=getNumberOfTuples();
3341 mcIdType nbOfTuple2=other->getNumberOfTuples();
3342 std::size_t nbOfComp=getNumberOfComponents();
3343 std::size_t nbOfComp2=other->getNumberOfComponents();
3344 if(nbOfTuple!=nbOfTuple2)
3345 throw INTERP_KERNEL::Exception("DataArrayDouble::powEqual : number of tuples mismatches !");
3346 if(nbOfComp!=1 || nbOfComp2!=1)
3347 throw INTERP_KERNEL::Exception("DataArrayDouble::powEqual : number of components of both arrays must be equal to 1 !");
3348 double *ptr=getPointer();
3349 const double *ptrc=other->begin();
3350 for(mcIdType i=0;i<nbOfTuple;i++,ptrc++,ptr++)
3353 *ptr=pow(*ptr,*ptrc);
3356 std::ostringstream oss; oss << "DataArrayDouble::powEqual : on tuple #" << i << " of this value is < 0 (" << *ptr << ") !";
3357 throw INTERP_KERNEL::Exception(oss.str().c_str());
3364 * This method is \b NOT wrapped into python because it can be useful only for performance reasons in C++ context.
3365 * All values in \a this must be 0. or 1. within eps error. 0 means false, 1 means true.
3366 * If an another value than 0 or 1 appear (within eps precision) an INTERP_KERNEL::Exception will be thrown.
3368 * \throw if \a this is not allocated.
3369 * \throw if \a this has not exactly one component.
3371 std::vector<bool> DataArrayDouble::toVectorOfBool(double eps) const
3374 if(getNumberOfComponents()!=1)
3375 throw INTERP_KERNEL::Exception("DataArrayDouble::toVectorOfBool : must be applied on single component array !");
3376 mcIdType nbt(getNumberOfTuples());
3377 std::vector<bool> ret(nbt);
3378 const double *pt(begin());
3379 for(mcIdType i=0;i<nbt;i++)
3383 else if(fabs(pt[i]-1.)<eps)
3387 std::ostringstream oss; oss << "DataArrayDouble::toVectorOfBool : the tuple #" << i << " has value " << pt[i] << " is invalid ! must be 0. or 1. !";
3388 throw INTERP_KERNEL::Exception(oss.str().c_str());
3395 * Useless method for end user. Only for MPI/Corba/File serialsation for multi arrays class.
3398 void DataArrayDouble::getTinySerializationIntInformation(std::vector<mcIdType>& tinyInfo) const
3403 tinyInfo[0]=getNumberOfTuples();
3404 tinyInfo[1]=ToIdType(getNumberOfComponents());
3414 * Useless method for end user. Only for MPI/Corba/File serialsation for multi arrays class.
3417 void DataArrayDouble::getTinySerializationStrInformation(std::vector<std::string>& tinyInfo) const
3421 std::size_t nbOfCompo(getNumberOfComponents());
3422 tinyInfo.resize(nbOfCompo+1);
3423 tinyInfo[0]=getName();
3424 for(std::size_t i=0;i<nbOfCompo;i++)
3425 tinyInfo[i+1]=getInfoOnComponent(i);
3430 tinyInfo[0]=getName();
3435 * Useless method for end user. Only for MPI/Corba/File serialsation for multi arrays class.
3436 * This method returns if a feeding is needed.
3438 bool DataArrayDouble::resizeForUnserialization(const std::vector<mcIdType>& tinyInfoI)
3440 mcIdType nbOfTuple=tinyInfoI[0];
3441 mcIdType nbOfComp=tinyInfoI[1];
3442 if(nbOfTuple!=-1 || nbOfComp!=-1)
3444 alloc(nbOfTuple,nbOfComp);
3451 * Useless method for end user. Only for MPI/Corba/File serialsation for multi arrays class.
3453 void DataArrayDouble::finishUnserialization(const std::vector<mcIdType>& tinyInfoI, const std::vector<std::string>& tinyInfoS)
3455 setName(tinyInfoS[0]);
3458 std::size_t nbOfCompo(getNumberOfComponents());
3459 for(std::size_t i=0;i<nbOfCompo;i++)
3460 setInfoOnComponent(i,tinyInfoS[i+1]);
3465 * Low static method that operates 3D rotation of 'nbNodes' 3D nodes whose coordinates are arranged in \a coordsIn
3466 * around an axe ( \a center, \a vect) and with angle \a angle.
3468 void DataArrayDouble::Rotate3DAlg(const double *center, const double *vect, double angle, mcIdType nbNodes, const double *coordsIn, double *coordsOut)
3470 if(!center || !vect)
3471 throw INTERP_KERNEL::Exception("DataArrayDouble::Rotate3DAlg : null vector in input !");
3472 double sina(sin(angle));
3473 double cosa(cos(angle));
3474 double vectorNorm[3];
3476 double matrixTmp[9];
3477 double norm(sqrt(vect[0]*vect[0]+vect[1]*vect[1]+vect[2]*vect[2]));
3478 if(norm<std::numeric_limits<double>::min())
3479 throw INTERP_KERNEL::Exception("DataArrayDouble::Rotate3DAlg : magnitude of input vector is too close of 0. !");
3480 std::transform(vect,vect+3,vectorNorm,std::bind(std::multiplies<double>(),std::placeholders::_1,1/norm));
3481 //rotation matrix computation
3482 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;
3483 matrixTmp[0]=vectorNorm[0]*vectorNorm[0]; matrixTmp[1]=vectorNorm[0]*vectorNorm[1]; matrixTmp[2]=vectorNorm[0]*vectorNorm[2];
3484 matrixTmp[3]=vectorNorm[1]*vectorNorm[0]; matrixTmp[4]=vectorNorm[1]*vectorNorm[1]; matrixTmp[5]=vectorNorm[1]*vectorNorm[2];
3485 matrixTmp[6]=vectorNorm[2]*vectorNorm[0]; matrixTmp[7]=vectorNorm[2]*vectorNorm[1]; matrixTmp[8]=vectorNorm[2]*vectorNorm[2];
3486 std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind(std::multiplies<double>(),std::placeholders::_1,1-cosa));
3487 std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
3488 matrixTmp[0]=0.; matrixTmp[1]=-vectorNorm[2]; matrixTmp[2]=vectorNorm[1];
3489 matrixTmp[3]=vectorNorm[2]; matrixTmp[4]=0.; matrixTmp[5]=-vectorNorm[0];
3490 matrixTmp[6]=-vectorNorm[1]; matrixTmp[7]=vectorNorm[0]; matrixTmp[8]=0.;
3491 std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind(std::multiplies<double>(),std::placeholders::_1,sina));
3492 std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
3493 //rotation matrix computed.
3495 for(mcIdType i=0; i<nbNodes; i++)
3497 std::transform(coordsIn+i*3,coordsIn+(i+1)*3,center,tmp,std::minus<double>());
3498 coordsOut[i*3]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+matrix[2]*tmp[2]+center[0];
3499 coordsOut[i*3+1]=matrix[3]*tmp[0]+matrix[4]*tmp[1]+matrix[5]*tmp[2]+center[1];
3500 coordsOut[i*3+2]=matrix[6]*tmp[0]+matrix[7]*tmp[1]+matrix[8]*tmp[2]+center[2];
3504 void DataArrayDouble::Symmetry3DPlane(const double point[3], const double normalVector[3], mcIdType nbNodes, const double *coordsIn, double *coordsOut)
3506 double matrix[9],matrix2[9],matrix3[9];
3507 double vect[3],crossVect[3];
3508 INTERP_KERNEL::orthogonalVect3(normalVector,vect);
3509 crossVect[0]=normalVector[1]*vect[2]-normalVector[2]*vect[1];
3510 crossVect[1]=normalVector[2]*vect[0]-normalVector[0]*vect[2];
3511 crossVect[2]=normalVector[0]*vect[1]-normalVector[1]*vect[0];
3512 double nv(INTERP_KERNEL::norm<3>(vect)),ni(INTERP_KERNEL::norm<3>(normalVector)),nc(INTERP_KERNEL::norm<3>(crossVect));
3513 matrix[0]=vect[0]/nv; matrix[1]=crossVect[0]/nc; matrix[2]=-normalVector[0]/ni;
3514 matrix[3]=vect[1]/nv; matrix[4]=crossVect[1]/nc; matrix[5]=-normalVector[1]/ni;
3515 matrix[6]=vect[2]/nv; matrix[7]=crossVect[2]/nc; matrix[8]=-normalVector[2]/ni;
3516 matrix2[0]=vect[0]/nv; matrix2[1]=vect[1]/nv; matrix2[2]=vect[2]/nv;
3517 matrix2[3]=crossVect[0]/nc; matrix2[4]=crossVect[1]/nc; matrix2[5]=crossVect[2]/nc;
3518 matrix2[6]=normalVector[0]/ni; matrix2[7]=normalVector[1]/ni; matrix2[8]=normalVector[2]/ni;
3519 for(mcIdType i=0;i<3;i++)
3520 for(mcIdType j=0;j<3;j++)
3523 for(mcIdType k=0;k<3;k++)
3524 val+=matrix[3*i+k]*matrix2[3*k+j];
3527 //rotation matrix computed.
3529 for(mcIdType i=0; i<nbNodes; i++)
3531 std::transform(coordsIn+i*3,coordsIn+(i+1)*3,point,tmp,std::minus<double>());
3532 coordsOut[i*3]=matrix3[0]*tmp[0]+matrix3[1]*tmp[1]+matrix3[2]*tmp[2]+point[0];
3533 coordsOut[i*3+1]=matrix3[3]*tmp[0]+matrix3[4]*tmp[1]+matrix3[5]*tmp[2]+point[1];
3534 coordsOut[i*3+2]=matrix3[6]*tmp[0]+matrix3[7]*tmp[1]+matrix3[8]*tmp[2]+point[2];
3538 void DataArrayDouble::GiveBaseForPlane(const double normalVector[3], double baseOfPlane[9])
3540 double vect[3],crossVect[3];
3541 INTERP_KERNEL::orthogonalVect3(normalVector,vect);
3542 crossVect[0]=normalVector[1]*vect[2]-normalVector[2]*vect[1];
3543 crossVect[1]=normalVector[2]*vect[0]-normalVector[0]*vect[2];
3544 crossVect[2]=normalVector[0]*vect[1]-normalVector[1]*vect[0];
3545 double nv(INTERP_KERNEL::norm<3>(vect)),ni(INTERP_KERNEL::norm<3>(normalVector)),nc(INTERP_KERNEL::norm<3>(crossVect));
3546 baseOfPlane[0]=vect[0]/nv; baseOfPlane[1]=vect[1]/nv; baseOfPlane[2]=vect[2]/nv;
3547 baseOfPlane[3]=crossVect[0]/nc; baseOfPlane[4]=crossVect[1]/nc; baseOfPlane[5]=crossVect[2]/nc;
3548 baseOfPlane[6]=normalVector[0]/ni; baseOfPlane[7]=normalVector[1]/ni; baseOfPlane[8]=normalVector[2]/ni;
3552 * \param [in] seg2 : coordinates of input seg2 expected to have spacedim==2
3553 * \param [in] tri3 : coordinates of input tri3 also expected to have spacedim==2
3554 * \param [out] coeffs : the result of integration normalized to 1. along \a seg2 inside tri3 sorted by the node id of \a tri3
3555 * \param [out] length : the length of seg2. That is too say the length of integration
3557 void DataArrayDouble::ComputeIntegralOfSeg2IntoTri3(const double seg2[4], const double tri3[6], double coeffs[3], double& length)
3559 length=INTERP_KERNEL::norme_vecteur(seg2,seg2+2);
3561 INTERP_KERNEL::mid_of_seg2(seg2,seg2+2,mid);
3562 INTERP_KERNEL::barycentric_coords<2>(tri3,mid,coeffs); // integral along seg2 is equal to value at the center of SEG2 !
3566 * Low static method that operates 3D rotation of \a nbNodes 3D nodes whose coordinates are arranged in \a coords
3567 * around the center point \a center and with angle \a angle.
3569 void DataArrayDouble::Rotate2DAlg(const double *center, double angle, mcIdType nbNodes, const double *coordsIn, double *coordsOut)
3571 double cosa=cos(angle);
3572 double sina=sin(angle);
3574 matrix[0]=cosa; matrix[1]=-sina; matrix[2]=sina; matrix[3]=cosa;
3576 for(mcIdType i=0; i<nbNodes; i++)
3578 std::transform(coordsIn+i*2,coordsIn+(i+1)*2,center,tmp,std::minus<double>());
3579 coordsOut[i*2]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+center[0];
3580 coordsOut[i*2+1]=matrix[2]*tmp[0]+matrix[3]*tmp[1]+center[1];
3584 DataArrayDoubleIterator::DataArrayDoubleIterator(DataArrayDouble *da):DataArrayIterator<double>(da)
3588 DataArrayDoubleTuple::DataArrayDoubleTuple(double *pt, std::size_t nbOfComp):DataArrayTuple<double>(pt,nbOfComp)
3593 std::string DataArrayDoubleTuple::repr() const
3595 std::ostringstream oss; oss.precision(17); oss << "(";
3596 for(std::size_t i=0;i<_nb_of_compo-1;i++)
3597 oss << _pt[i] << ", ";
3598 oss << _pt[_nb_of_compo-1] << ")";
3602 double DataArrayDoubleTuple::doubleValue() const
3604 return this->zeValue();
3608 * This method returns a newly allocated instance the caller should dealed with by a MEDCoupling::DataArrayDouble::decrRef.
3609 * This method performs \b no copy of data. The content is only referenced using MEDCoupling::DataArrayDouble::useArray with ownership set to \b false.
3610 * 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
3611 * \b nbOfCompo=1 and \bnbOfTuples==this->_nb_of_elem.
3613 DataArrayDouble *DataArrayDoubleTuple::buildDADouble(std::size_t nbOfTuples, std::size_t nbOfCompo) const
3615 return this->buildDA(nbOfTuples,nbOfCompo);
3619 * Returns a full copy of \a this. For more info on copying data arrays see
3620 * \ref MEDCouplingArrayBasicsCopyDeep.
3621 * \return DataArrayInt * - a new instance of DataArrayInt.
3623 DataArrayInt32 *DataArrayInt32::deepCopy() const
3625 return new DataArrayInt32(*this);
3628 MCAuto<DataArrayInt64> DataArrayInt32::convertToInt64Arr() const
3630 this->checkAllocated();
3631 MCAuto<DataArrayInt64> ret(DataArrayInt64::New());
3632 ret->alloc(this->getNumberOfTuples(),this->getNumberOfComponents());
3633 ret->copyStringInfoFrom(*this);
3634 const std::int32_t *pt(this->begin());
3635 std::for_each(ret->getPointer(),ret->getPointer()+ret->getNbOfElems(),[&pt](std::int64_t& val) { val = std::int64_t(*pt++); });
3639 DataArrayInt32Iterator *DataArrayInt32::iterator()
3641 return new DataArrayInt32Iterator(this);
3645 DataArrayInt32Iterator::DataArrayInt32Iterator(DataArrayInt32 *da):DataArrayIterator<Int32>(da)
3649 DataArrayInt32Tuple::DataArrayInt32Tuple(Int32 *pt, std::size_t nbOfComp):DataArrayTuple<Int32>(pt,nbOfComp)
3653 std::string DataArrayInt32Tuple::repr() const
3655 std::ostringstream oss; oss << "(";
3656 for(std::size_t i=0;i<_nb_of_compo-1;i++)
3657 oss << _pt[i] << ", ";
3658 oss << _pt[_nb_of_compo-1] << ")";
3662 Int32 DataArrayInt32Tuple::intValue() const
3664 return this->zeValue();
3668 * This method returns a newly allocated instance the caller should dealed with by a MEDCoupling::DataArrayInt::decrRef.
3669 * This method performs \b no copy of data. The content is only referenced using MEDCoupling::DataArrayInt::useArray with ownership set to \b false.
3670 * 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
3671 * \b nbOfCompo=1 and \bnbOfTuples==this->_nb_of_elem.
3673 DataArrayInt32 *DataArrayInt32Tuple::buildDAInt(std::size_t nbOfTuples, std::size_t nbOfCompo) const
3675 return this->buildDA(nbOfTuples,nbOfCompo);
3678 MCAuto<DataArrayInt32> DataArrayInt64::convertToInt32Arr() const
3680 this->checkAllocated();
3681 MCAuto<DataArrayInt32> ret(DataArrayInt32::New());
3682 ret->alloc(this->getNumberOfTuples(),this->getNumberOfComponents());
3683 ret->copyStringInfoFrom(*this);
3684 const std::int64_t *pt(this->begin());
3685 std::for_each(ret->getPointer(),ret->getPointer()+ret->getNbOfElems(),[&pt](std::int32_t& val) { val = std::int32_t(*pt++); });
3689 DataArrayInt64Iterator *DataArrayInt64::iterator()
3691 return new DataArrayInt64Iterator(this);
3695 DataArrayInt64Iterator::DataArrayInt64Iterator(DataArrayInt64 *da):DataArrayIterator<Int64>(da)
3699 DataArrayInt64Tuple::DataArrayInt64Tuple(Int64 *pt, std::size_t nbOfComp):DataArrayTuple<Int64>(pt,nbOfComp)
3703 std::string DataArrayInt64Tuple::repr() const
3705 std::ostringstream oss; oss << "(";
3706 for(std::size_t i=0;i<_nb_of_compo-1;i++)
3707 oss << _pt[i] << ", ";
3708 oss << _pt[_nb_of_compo-1] << ")";
3712 Int64 DataArrayInt64Tuple::intValue() const
3714 return this->zeValue();
3717 DataArrayInt64 *DataArrayInt64Tuple::buildDAInt(std::size_t nbOfTuples, std::size_t nbOfCompo) const
3719 return this->buildDA(nbOfTuples,nbOfCompo);
3723 DataArrayInt64 *DataArrayInt64::deepCopy() const
3725 return new DataArrayInt64(*this);