1 // Copyright (C) 2007-2019 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::MemArray<mcIdType>;
46 template class MEDCoupling::MemArray<double>;
47 template class MEDCoupling::DataArrayTemplate<mcIdType>;
48 template class MEDCoupling::DataArrayTemplate<double>;
49 template class MEDCoupling::DataArrayTemplateClassic<mcIdType>;
50 template class MEDCoupling::DataArrayTemplateClassic<double>;
51 template class MEDCoupling::DataArrayTemplateFP<double>;
52 template class MEDCoupling::DataArrayIterator<double>;
53 template class MEDCoupling::DataArrayIterator<mcIdType>;
54 template class MEDCoupling::DataArrayDiscrete<Int32>;
55 template class MEDCoupling::DataArrayDiscreteSigned<Int32>;
56 template class MEDCoupling::DataArrayDiscrete<Int64>;
57 template class MEDCoupling::DataArrayDiscreteSigned<Int64>;
58 template class MEDCoupling::DataArrayTuple<mcIdType>;
59 template class MEDCoupling::DataArrayTuple<double>;
60 template class MEDCoupling::DataArrayTuple<float>;
62 template<mcIdType SPACEDIM>
63 void DataArrayDouble::findCommonTuplesAlg(const double *bbox, mcIdType nbNodes, mcIdType limitNodeId, double prec, DataArrayIdType *c, DataArrayIdType *cI) const
65 const double *coordsPtr=getConstPointer();
66 BBTreePts<SPACEDIM,mcIdType> myTree(bbox,0,0,nbNodes,prec);
67 std::vector<bool> isDone(nbNodes);
68 for(mcIdType i=0;i<nbNodes;i++)
72 std::vector<mcIdType> intersectingElems;
73 myTree.getElementsAroundPoint(coordsPtr+i*SPACEDIM,intersectingElems);
74 if(intersectingElems.size()>1)
76 std::vector<mcIdType> commonNodes;
77 for(std::vector<mcIdType>::const_iterator it=intersectingElems.begin();it!=intersectingElems.end();it++)
81 commonNodes.push_back(*it);
84 if(!commonNodes.empty())
86 cI->pushBackSilent(cI->back()+ToIdType(commonNodes.size())+1);
88 c->insertAtTheEnd(commonNodes.begin(),commonNodes.end());
95 template<mcIdType SPACEDIM>
96 void DataArrayDouble::FindTupleIdsNearTuplesAlg(const BBTreePts<SPACEDIM,mcIdType>& myTree, const double *pos, mcIdType nbOfTuples, double eps,
97 DataArrayIdType *c, DataArrayIdType *cI)
99 for(mcIdType i=0;i<nbOfTuples;i++)
101 std::vector<mcIdType> intersectingElems;
102 myTree.getElementsAroundPoint(pos+i*SPACEDIM,intersectingElems);
103 std::vector<mcIdType> commonNodes;
104 for(std::vector<mcIdType>::const_iterator it=intersectingElems.begin();it!=intersectingElems.end();it++)
105 commonNodes.push_back(*it);
106 cI->pushBackSilent(cI->back()+ToIdType(commonNodes.size()));
107 c->insertAtTheEnd(commonNodes.begin(),commonNodes.end());
111 template<mcIdType SPACEDIM>
112 void DataArrayDouble::FindClosestTupleIdAlg(const BBTreePts<SPACEDIM,mcIdType>& myTree, double dist, const double *pos, mcIdType nbOfTuples, const double *thisPt, mcIdType thisNbOfTuples, mcIdType *res)
114 double distOpt(dist);
115 const double *p(pos);
117 for(mcIdType i=0;i<nbOfTuples;i++,p+=SPACEDIM,r++)
122 double ret=myTree.getElementsAroundPoint2(p,distOpt,elem);
123 if(ret!=std::numeric_limits<double>::max())
125 distOpt=std::max(ret,1e-4);
130 { distOpt=2*distOpt; continue; }
135 mcIdType DataArray::EffectiveCircPerm(mcIdType nbOfShift, mcIdType nbOfTuples)
138 throw INTERP_KERNEL::Exception("DataArray::EffectiveCircPerm : number of tuples is expected to be > 0 !");
141 return nbOfShift%nbOfTuples;
145 mcIdType tmp(-nbOfShift);
147 return nbOfTuples-tmp;
151 std::size_t DataArray::getHeapMemorySizeWithoutChildren() const
153 std::size_t sz1=_name.capacity();
154 std::size_t sz2=_info_on_compo.capacity();
156 for(std::vector<std::string>::const_iterator it=_info_on_compo.begin();it!=_info_on_compo.end();it++)
157 sz3+=(*it).capacity();
161 std::vector<const BigMemoryObject *> DataArray::getDirectChildrenWithNull() const
163 return std::vector<const BigMemoryObject *>();
167 * Sets the attribute \a _name of \a this array.
168 * See \ref MEDCouplingArrayBasicsName "DataArrays infos" for more information.
169 * \param [in] name - new array name
171 void DataArray::setName(const std::string& name)
177 * Copies textual data from an \a other DataArray. The copied data are
178 * - the name attribute,
179 * - the information of components.
181 * For more information on these data see \ref MEDCouplingArrayBasicsName "DataArrays infos".
183 * \param [in] other - another instance of DataArray to copy the textual data from.
184 * \throw If number of components of \a this array differs from that of the \a other.
186 void DataArray::copyStringInfoFrom(const DataArray& other)
188 if(_info_on_compo.size()!=other._info_on_compo.size())
189 throw INTERP_KERNEL::Exception("Size of arrays mismatches on copyStringInfoFrom !");
191 _info_on_compo=other._info_on_compo;
194 void DataArray::copyPartOfStringInfoFrom(const DataArray& other, const std::vector<std::size_t>& compoIds)
196 std::size_t nbOfCompoOth=other.getNumberOfComponents();
197 std::size_t newNbOfCompo=compoIds.size();
198 for(std::size_t i=0;i<newNbOfCompo;i++)
199 if(compoIds[i]>=nbOfCompoOth || compoIds[i]<0)
201 std::ostringstream oss; oss << "Specified component id is out of range (" << compoIds[i] << ") compared with nb of actual components (" << nbOfCompoOth << ")";
202 throw INTERP_KERNEL::Exception(oss.str().c_str());
204 for(std::size_t i=0;i<newNbOfCompo;i++)
205 setInfoOnComponent(i,other.getInfoOnComponent(compoIds[i]));
208 void DataArray::copyPartOfStringInfoFrom2(const std::vector<std::size_t>& compoIds, const DataArray& other)
210 if(compoIds.size()!=other.getNumberOfComponents())
211 throw INTERP_KERNEL::Exception("Given compoIds has not the same size as number of components of given array !");
212 std::size_t partOfCompoToSet=compoIds.size();
213 std::size_t nbOfCompo=getNumberOfComponents();
214 for(std::size_t i=0;i<partOfCompoToSet;i++)
215 if(compoIds[i]>=nbOfCompo || compoIds[i]<0)
217 std::ostringstream oss; oss << "Specified component id is out of range (" << compoIds[i] << ") compared with nb of actual components (" << nbOfCompo << ")";
218 throw INTERP_KERNEL::Exception(oss.str().c_str());
220 for(std::size_t i=0;i<partOfCompoToSet;i++)
221 setInfoOnComponent(compoIds[i],other.getInfoOnComponent(i));
224 bool DataArray::areInfoEqualsIfNotWhy(const DataArray& other, std::string& reason) const
226 std::ostringstream oss;
227 if(_name!=other._name)
229 oss << "Names DataArray mismatch : this name=\"" << _name << " other name=\"" << other._name << "\" !";
233 if(_info_on_compo!=other._info_on_compo)
235 oss << "Components DataArray mismatch : \nThis components=";
236 for(std::vector<std::string>::const_iterator it=_info_on_compo.begin();it!=_info_on_compo.end();it++)
237 oss << "\"" << *it << "\",";
238 oss << "\nOther components=";
239 for(std::vector<std::string>::const_iterator it=other._info_on_compo.begin();it!=other._info_on_compo.end();it++)
240 oss << "\"" << *it << "\",";
248 * Compares textual information of \a this DataArray with that of an \a other one.
249 * The compared data are
250 * - the name attribute,
251 * - the information of components.
253 * For more information on these data see \ref MEDCouplingArrayBasicsName "DataArrays infos".
254 * \param [in] other - another instance of DataArray to compare the textual data of.
255 * \return bool - \a true if the textual information is same, \a false else.
257 bool DataArray::areInfoEquals(const DataArray& other) const
260 return areInfoEqualsIfNotWhy(other,tmp);
263 void DataArray::reprWithoutNameStream(std::ostream& stream) const
265 stream << "Number of components : "<< getNumberOfComponents() << "\n";
266 stream << "Info of these components : ";
267 for(std::vector<std::string>::const_iterator iter=_info_on_compo.begin();iter!=_info_on_compo.end();iter++)
268 stream << "\"" << *iter << "\" ";
272 std::string DataArray::cppRepr(const std::string& varName) const
274 std::ostringstream ret;
275 reprCppStream(varName,ret);
280 * Sets information on all components. To know more on format of this information
281 * see \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
282 * \param [in] info - a vector of strings.
283 * \throw If size of \a info differs from the number of components of \a this.
285 void DataArray::setInfoOnComponents(const std::vector<std::string>& info)
287 if(getNumberOfComponents()!=info.size())
289 std::ostringstream oss; oss << "DataArray::setInfoOnComponents : input is of size " << info.size() << " whereas number of components is equal to " << getNumberOfComponents() << " !";
290 throw INTERP_KERNEL::Exception(oss.str().c_str());
296 * This method is only a dispatcher towards DataArrayDouble::setPartOfValues3, DataArrayInt::setPartOfValues3, DataArrayChar::setPartOfValues3 depending on the true
297 * type of \a this and \a aBase.
299 * \throw If \a aBase and \a this do not have the same type.
301 * \sa DataArrayDouble::setPartOfValues3, DataArrayInt::setPartOfValues3, DataArrayChar::setPartOfValues3.
303 void DataArray::setPartOfValuesBase3(const DataArray *aBase, const mcIdType *bgTuples, const mcIdType *endTuples, mcIdType bgComp, mcIdType endComp, mcIdType stepComp, bool strictCompoCompare)
306 throw INTERP_KERNEL::Exception("DataArray::setPartOfValuesBase3 : input aBase object is NULL !");
307 DataArrayDouble *this1(dynamic_cast<DataArrayDouble *>(this));
308 DataArrayIdType *this2(dynamic_cast<DataArrayIdType *>(this));
309 DataArrayChar *this3(dynamic_cast<DataArrayChar *>(this));
310 const DataArrayDouble *a1(dynamic_cast<const DataArrayDouble *>(aBase));
311 const DataArrayIdType *a2(dynamic_cast<const DataArrayIdType *>(aBase));
312 const DataArrayChar *a3(dynamic_cast<const DataArrayChar *>(aBase));
315 this1->setPartOfValues3(a1,bgTuples,endTuples,bgComp,endComp,stepComp,strictCompoCompare);
320 this2->setPartOfValues3(a2,bgTuples,endTuples,bgComp,endComp,stepComp,strictCompoCompare);
325 this3->setPartOfValues3(a3,bgTuples,endTuples,bgComp,endComp,stepComp,strictCompoCompare);
328 throw INTERP_KERNEL::Exception("DataArray::setPartOfValuesBase3 : input aBase object and this do not have the same type !");
331 std::vector<std::string> DataArray::getVarsOnComponent() const
333 std::size_t nbOfCompo=_info_on_compo.size();
334 std::vector<std::string> ret(nbOfCompo);
335 for(std::size_t i=0;i<nbOfCompo;i++)
336 ret[i]=getVarOnComponent(i);
340 std::vector<std::string> DataArray::getUnitsOnComponent() const
342 std::size_t nbOfCompo=_info_on_compo.size();
343 std::vector<std::string> ret(nbOfCompo);
344 for(std::size_t i=0;i<nbOfCompo;i++)
345 ret[i]=getUnitOnComponent(i);
350 * Returns information on a component specified by an index.
351 * To know more on format of this information
352 * see \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
353 * \param [in] i - the index (zero based) of the component of interest.
354 * \return std::string - a string containing the information on \a i-th component.
355 * \throw If \a i is not a valid component index.
357 std::string DataArray::getInfoOnComponent(std::size_t i) const
359 if(i<_info_on_compo.size())
360 return _info_on_compo[i];
363 std::ostringstream oss; oss << "DataArray::getInfoOnComponent : Specified component id is out of range (" << i << ") compared with nb of actual components (" << _info_on_compo.size();
364 throw INTERP_KERNEL::Exception(oss.str().c_str());
369 * Returns the var part of the full information of the \a i-th component.
370 * For example, if \c getInfoOnComponent(0) returns "SIGXY [N/m^2]", then
371 * \c getVarOnComponent(0) returns "SIGXY".
372 * If a unit part of information is not detected by presence of
373 * two square brackets, then the full information is returned.
374 * To read more about the component information format, see
375 * \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
376 * \param [in] i - the index (zero based) of the component of interest.
377 * \return std::string - a string containing the var information, or the full info.
378 * \throw If \a i is not a valid component index.
380 std::string DataArray::getVarOnComponent(std::size_t i) const
382 if(i<_info_on_compo.size())
384 return GetVarNameFromInfo(_info_on_compo[i]);
388 std::ostringstream oss; oss << "DataArray::getVarOnComponent : Specified component id is out of range (" << i << ") compared with nb of actual components (" << _info_on_compo.size();
389 throw INTERP_KERNEL::Exception(oss.str().c_str());
394 * Returns the unit part of the full information of the \a i-th component.
395 * For example, if \c getInfoOnComponent(0) returns "SIGXY [ N/m^2]", then
396 * \c getUnitOnComponent(0) returns " N/m^2".
397 * If a unit part of information is not detected by presence of
398 * two square brackets, then an empty string is returned.
399 * To read more about the component information format, see
400 * \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
401 * \param [in] i - the index (zero based) of the component of interest.
402 * \return std::string - a string containing the unit information, if any, or "".
403 * \throw If \a i is not a valid component index.
405 std::string DataArray::getUnitOnComponent(std::size_t i) const
407 if(i<_info_on_compo.size())
409 return GetUnitFromInfo(_info_on_compo[i]);
413 std::ostringstream oss; oss << "DataArray::getUnitOnComponent : Specified component id is out of range (" << i << ") compared with nb of actual components (" << _info_on_compo.size();
414 throw INTERP_KERNEL::Exception(oss.str().c_str());
419 * Returns the var part of the full component information.
420 * For example, if \a info == "SIGXY [N/m^2]", then this method returns "SIGXY".
421 * If a unit part of information is not detected by presence of
422 * two square brackets, then the whole \a info is returned.
423 * To read more about the component information format, see
424 * \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
425 * \param [in] info - the full component information.
426 * \return std::string - a string containing only var information, or the \a info.
428 std::string DataArray::GetVarNameFromInfo(const std::string& info)
430 std::size_t p1=info.find_last_of('[');
431 std::size_t p2=info.find_last_of(']');
432 if(p1==std::string::npos || p2==std::string::npos)
437 return std::string();
438 std::size_t p3=info.find_last_not_of(' ',p1-1);
439 return info.substr(0,p3+1);
443 * Returns the unit part of the full component information.
444 * For example, if \a info == "SIGXY [ N/m^2]", then this method returns " N/m^2".
445 * If a unit part of information is not detected by presence of
446 * two square brackets, then an empty string is returned.
447 * To read more about the component information format, see
448 * \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
449 * \param [in] info - the full component information.
450 * \return std::string - a string containing only unit information, if any, or "".
452 std::string DataArray::GetUnitFromInfo(const std::string& info)
454 std::size_t p1=info.find_last_of('[');
455 std::size_t p2=info.find_last_of(']');
456 if(p1==std::string::npos || p2==std::string::npos)
457 return std::string();
459 return std::string();
460 return info.substr(p1+1,p2-p1-1);
464 * This method put in info format the result of the merge of \a var and \a unit.
465 * The standard format for that is "var [unit]".
466 * Inversely you can retrieve the var part or the unit part of info string using resp. GetVarNameFromInfo and GetUnitFromInfo.
468 std::string DataArray::BuildInfoFromVarAndUnit(const std::string& var, const std::string& unit)
470 std::ostringstream oss;
471 oss << var << " [" << unit << "]";
475 std::string DataArray::GetAxisTypeRepr(MEDCouplingAxisType at)
480 return std::string("AX_CART");
482 return std::string("AX_CYL");
484 return std::string("AX_SPHER");
486 throw INTERP_KERNEL::Exception("DataArray::GetAxisTypeRepr : unrecognized axis type enum !");
491 * Returns a new DataArray by concatenating all given arrays, so that (1) the number
492 * of tuples in the result array is a sum of the number of tuples of given arrays and (2)
493 * the number of component in the result array is same as that of each of given arrays.
494 * Info on components is copied from the first of the given arrays. Number of components
495 * in the given arrays must be the same.
496 * \param [in] arrs - a sequence of arrays to include in the result array. All arrays must have the same type.
497 * \return DataArray * - the new instance of DataArray (that can be either DataArrayInt, DataArrayDouble, DataArrayChar).
498 * The caller is to delete this result array using decrRef() as it is no more
500 * \throw If all arrays within \a arrs are NULL.
501 * \throw If all not null arrays in \a arrs have not the same type.
502 * \throw If getNumberOfComponents() of arrays within \a arrs.
504 DataArray *DataArray::Aggregate(const std::vector<const DataArray *>& arrs)
506 std::vector<const DataArray *> arr2;
507 for(std::vector<const DataArray *>::const_iterator it=arrs.begin();it!=arrs.end();it++)
511 throw INTERP_KERNEL::Exception("DataArray::Aggregate : only null instance in input vector !");
512 std::vector<const DataArrayDouble *> arrd;
513 std::vector<const DataArrayIdType *> arri;
514 std::vector<const DataArrayChar *> arrc;
515 for(std::vector<const DataArray *>::const_iterator it=arr2.begin();it!=arr2.end();it++)
517 const DataArrayDouble *a=dynamic_cast<const DataArrayDouble *>(*it);
519 { arrd.push_back(a); continue; }
520 const DataArrayIdType *b=dynamic_cast<const DataArrayIdType *>(*it);
522 { arri.push_back(b); continue; }
523 const DataArrayChar *c=dynamic_cast<const DataArrayChar *>(*it);
525 { arrc.push_back(c); continue; }
526 throw INTERP_KERNEL::Exception("DataArray::Aggregate : presence of not null instance in inuput that is not in [DataArrayDouble, DataArrayInt, DataArrayChar] !");
528 if(arr2.size()==arrd.size())
529 return DataArrayDouble::Aggregate(arrd);
530 if(arr2.size()==arri.size())
531 return DataArrayIdType::Aggregate(arri);
532 if(arr2.size()==arrc.size())
533 return DataArrayChar::Aggregate(arrc);
534 throw INTERP_KERNEL::Exception("DataArray::Aggregate : all input arrays must have the same type !");
538 * Sets information on a component specified by an index.
539 * To know more on format of this information
540 * see \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
541 * \warning Don't pass NULL as \a info!
542 * \param [in] i - the index (zero based) of the component of interest.
543 * \param [in] info - the string containing the information.
544 * \throw If \a i is not a valid component index.
546 void DataArray::setInfoOnComponent(std::size_t i, const std::string& info)
548 if(i<_info_on_compo.size())
549 _info_on_compo[i]=info;
552 std::ostringstream oss; oss << "DataArray::setInfoOnComponent : Specified component id is out of range (" << i << ") compared with nb of actual components (" << _info_on_compo.size();
553 throw INTERP_KERNEL::Exception(oss.str().c_str());
558 * Sets information on all components. This method can change number of components
559 * at certain conditions; if the conditions are not respected, an exception is thrown.
560 * The number of components can be changed in \a this only if \a this is not allocated.
561 * The condition of number of components must not be changed.
563 * To know more on format of the component information see
564 * \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
565 * \param [in] info - a vector of component infos.
566 * \throw If \a this->getNumberOfComponents() != \a info.size() && \a this->isAllocated()
568 void DataArray::setInfoAndChangeNbOfCompo(const std::vector<std::string>& info)
570 if(getNumberOfComponents()!=info.size())
576 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 !";
577 throw INTERP_KERNEL::Exception(oss.str().c_str());
584 void DataArray::checkNbOfTuples(mcIdType nbOfTuples, const std::string& msg) const
586 if(getNumberOfTuples()!=nbOfTuples)
588 std::ostringstream oss; oss << msg << " : mismatch number of tuples : expected " << nbOfTuples << " having " << getNumberOfTuples() << " !";
589 throw INTERP_KERNEL::Exception(oss.str().c_str());
593 void DataArray::checkNbOfComps(std::size_t nbOfCompo, const std::string& msg) const
595 if (getNumberOfComponents()!=nbOfCompo)
597 std::ostringstream oss; oss << msg << " : mismatch number of components : expected " << nbOfCompo << " having " << getNumberOfComponents() << " !";
598 throw INTERP_KERNEL::Exception(oss.str().c_str());
602 void DataArray::checkNbOfElems(mcIdType nbOfElems, const std::string& msg) const
604 if(getNbOfElems()!=nbOfElems)
606 std::ostringstream oss; oss << msg << " : mismatch number of elems : Expected " << nbOfElems << " having " << getNbOfElems() << " !";
607 throw INTERP_KERNEL::Exception(oss.str().c_str());
611 void DataArray::checkNbOfTuplesAndComp(const DataArray& other, const std::string& msg) const
613 if(getNumberOfTuples()!=other.getNumberOfTuples())
615 std::ostringstream oss; oss << msg << " : mismatch number of tuples : expected " << other.getNumberOfTuples() << " having " << getNumberOfTuples() << " !";
616 throw INTERP_KERNEL::Exception(oss.str().c_str());
618 if(getNumberOfComponents()!=other.getNumberOfComponents())
620 std::ostringstream oss; oss << msg << " : mismatch number of components : expected " << other.getNumberOfComponents() << " having " << getNumberOfComponents() << " !";
621 throw INTERP_KERNEL::Exception(oss.str().c_str());
625 void DataArray::checkNbOfTuplesAndComp(mcIdType nbOfTuples, std::size_t nbOfCompo, const std::string& msg) const
627 checkNbOfTuples(nbOfTuples,msg);
628 checkNbOfComps(nbOfCompo,msg);
632 * Simply this method checks that \b value is in [0,\b ref).
634 void DataArray::CheckValueInRange(mcIdType ref, mcIdType value, const std::string& msg)
636 if(value<0 || value>=ref)
638 std::ostringstream oss; oss << "DataArray::CheckValueInRange : " << msg << " ! Expected in range [0," << ref << "[ having " << value << " !";
639 throw INTERP_KERNEL::Exception(oss.str().c_str());
644 * This method checks that [\b start, \b end) is compliant with ref length \b value.
645 * typically start in [0,\b value) and end in [0,\b value). If value==start and start==end, it is supported.
647 void DataArray::CheckValueInRangeEx(mcIdType value, mcIdType start, mcIdType end, const std::string& msg)
649 if(start<0 || start>=value)
651 if(value!=start || end!=start)
653 std::ostringstream oss; oss << "DataArray::CheckValueInRangeEx : " << msg << " ! Expected start " << start << " of input range, in [0," << value << "[ !";
654 throw INTERP_KERNEL::Exception(oss.str().c_str());
657 if(end<0 || end>value)
659 std::ostringstream oss; oss << "DataArray::CheckValueInRangeEx : " << msg << " ! Expected end " << end << " of input range, in [0," << value << "] !";
660 throw INTERP_KERNEL::Exception(oss.str().c_str());
664 void DataArray::CheckClosingParInRange(mcIdType ref, mcIdType value, const std::string& msg)
666 if(value<0 || value>ref)
668 std::ostringstream oss; oss << "DataArray::CheckClosingParInRange : " << msg << " ! Expected input range in [0," << ref << "] having closing open parenthesis " << value << " !";
669 throw INTERP_KERNEL::Exception(oss.str().c_str());
674 * 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,
675 * typically it is a whole slice of tuples of DataArray or cells, nodes of a mesh...
677 * The input \a sliceId should be an id in [0, \a nbOfSlices) that specifies the slice of work.
679 * \param [in] start - the start of the input slice of the whole work to perform split into slices.
680 * \param [in] stop - the stop of the input slice of the whole work to perform split into slices.
681 * \param [in] step - the step (that can be <0) of the input slice of the whole work to perform split into slices.
682 * \param [in] sliceId - the slice id considered
683 * \param [in] nbOfSlices - the number of slices (typically the number of cores on which the work is expected to be sliced)
684 * \param [out] startSlice - the start of the slice considered
685 * \param [out] stopSlice - the stop of the slice consided
687 * \throw If \a step == 0
688 * \throw If \a nbOfSlices not > 0
689 * \throw If \a sliceId not in [0,nbOfSlices)
691 void DataArray::GetSlice(mcIdType start, mcIdType stop, mcIdType step, mcIdType sliceId, mcIdType nbOfSlices, mcIdType& startSlice, mcIdType& stopSlice)
693 DataArrayTools<mcIdType>::GetSlice(start, stop, step, sliceId, nbOfSlices, startSlice, stopSlice);
696 mcIdType DataArray::GetNumberOfItemGivenBES(mcIdType begin, mcIdType end, mcIdType step, const std::string& msg)
698 return DataArrayTools<mcIdType>::GetNumberOfItemGivenBES(begin, end, step, msg);
701 mcIdType DataArray::GetNumberOfItemGivenBESRelative(mcIdType begin, mcIdType end, mcIdType step, const std::string& msg)
703 return DataArrayTools<mcIdType>::GetNumberOfItemGivenBESRelative(begin, end, step, msg);
706 mcIdType DataArray::GetPosOfItemGivenBESRelativeNoThrow(mcIdType value, mcIdType begin, mcIdType end, mcIdType step)
708 return DataArrayTools<mcIdType>::GetPosOfItemGivenBESRelativeNoThrow(value, begin, end, step);
712 * Returns a new instance of DataArrayDouble. The caller is to delete this array
713 * using decrRef() as it is no more needed.
715 DataArrayDouble *DataArrayDouble::New()
717 return new DataArrayDouble;
721 * Returns the only one value in \a this, if and only if number of elements
722 * (nb of tuples * nb of components) is equal to 1, and that \a this is allocated.
723 * \return double - the sole value stored in \a this array.
724 * \throw If at least one of conditions stated above is not fulfilled.
726 double DataArrayDouble::doubleValue() const
730 if(getNbOfElems()==1)
732 return *getConstPointer();
735 throw INTERP_KERNEL::Exception("DataArrayDouble::doubleValue : DataArrayDouble instance is allocated but number of elements is not equal to 1 !");
738 throw INTERP_KERNEL::Exception("DataArrayDouble::doubleValue : DataArrayDouble instance is not allocated !");
742 * Returns a full copy of \a this. For more info on copying data arrays see
743 * \ref MEDCouplingArrayBasicsCopyDeep.
744 * \return DataArrayDouble * - a new instance of DataArrayDouble. The caller is to
745 * delete this array using decrRef() as it is no more needed.
747 DataArrayDouble *DataArrayDouble::deepCopy() const
749 return new DataArrayDouble(*this);
753 * Checks that \a this array is consistently **increasing** or **decreasing** in value,
754 * with at least absolute difference value of |\a eps| at each step.
755 * If not an exception is thrown.
756 * \param [in] increasing - if \a true, the array values should be increasing.
757 * \param [in] eps - minimal absolute difference between the neighbor values at which
758 * the values are considered different.
759 * \throw If sequence of values is not strictly monotonic in agreement with \a
761 * \throw If \a this->getNumberOfComponents() != 1.
762 * \throw If \a this is not allocated.
764 void DataArrayDouble::checkMonotonic(bool increasing, double eps) const
766 if(!isMonotonic(increasing,eps))
769 throw INTERP_KERNEL::Exception("DataArrayDouble::checkMonotonic : 'this' is not INCREASING monotonic !");
771 throw INTERP_KERNEL::Exception("DataArrayDouble::checkMonotonic : 'this' is not DECREASING monotonic !");
776 * Checks that \a this array is consistently **increasing** or **decreasing** in value,
777 * with at least absolute difference value of |\a eps| at each step.
778 * \param [in] increasing - if \a true, array values should be increasing.
779 * \param [in] eps - minimal absolute difference between the neighbor values at which
780 * the values are considered different.
781 * \return bool - \a true if values change in accordance with \a increasing arg.
782 * \throw If \a this->getNumberOfComponents() != 1.
783 * \throw If \a this is not allocated.
785 bool DataArrayDouble::isMonotonic(bool increasing, double eps) const
788 if(getNumberOfComponents()!=1)
789 throw INTERP_KERNEL::Exception("DataArrayDouble::isMonotonic : only supported with 'this' array with ONE component !");
790 mcIdType nbOfElements(getNumberOfTuples());
791 const double *ptr=getConstPointer();
795 double absEps=fabs(eps);
798 for(mcIdType i=1;i<nbOfElements;i++)
800 if(ptr[i]<(ref+absEps))
808 for(mcIdType i=1;i<nbOfElements;i++)
810 if(ptr[i]>(ref-absEps))
818 void DataArrayDouble::writeVTK(std::ostream& ofs, mcIdType indent, const std::string& nameInFile, DataArrayByte *byteArr) const
820 static const char SPACE[4]={' ',' ',' ',' '};
822 std::string idt(indent,' ');
824 ofs << idt << "<DataArray type=\"Float32\" Name=\"" << nameInFile << "\" NumberOfComponents=\"" << getNumberOfComponents() << "\"";
826 bool areAllEmpty(true);
827 for(std::vector<std::string>::const_iterator it=_info_on_compo.begin();it!=_info_on_compo.end();it++)
831 for(std::size_t i=0;i<_info_on_compo.size();i++)
832 ofs << " ComponentName" << i << "=\"" << _info_on_compo[i] << "\"";
836 ofs << " format=\"appended\" offset=\"" << byteArr->getNumberOfTuples() << "\">";
837 INTERP_KERNEL::AutoPtr<float> tmp(new float[getNbOfElems()]);
839 // to make Visual C++ happy : instead of std::copy(begin(),end(),(float *)tmp);
840 for(const double *src=begin();src!=end();src++,pt++)
842 const char *data(reinterpret_cast<const char *>((float *)tmp));
843 std::size_t sz(getNbOfElems()*sizeof(float));
844 byteArr->insertAtTheEnd(data,data+sz);
845 byteArr->insertAtTheEnd(SPACE,SPACE+4);
849 ofs << " RangeMin=\"" << getMinValueInArray() << "\" RangeMax=\"" << getMaxValueInArray() << "\" format=\"ascii\">\n" << idt;
850 std::copy(begin(),end(),std::ostream_iterator<double>(ofs," "));
852 ofs << std::endl << idt << "</DataArray>\n";
855 void DataArrayDouble::reprCppStream(const std::string& varName, std::ostream& stream) const
857 mcIdType nbTuples=getNumberOfTuples();
858 std::size_t nbComp=getNumberOfComponents();
859 const double *data(getConstPointer());
860 stream.precision(17);
861 stream << "DataArrayDouble *" << varName << "=DataArrayDouble::New();" << std::endl;
862 if(nbTuples*nbComp>=1)
864 stream << "const double " << varName << "Data[" << nbTuples*nbComp << "]={";
865 std::copy(data,data+nbTuples*nbComp-1,std::ostream_iterator<double>(stream,","));
866 stream << data[nbTuples*nbComp-1] << "};" << std::endl;
867 stream << varName << "->useArray(" << varName << "Data,false,CPP_DEALLOC," << nbTuples << "," << nbComp << ");" << std::endl;
870 stream << varName << "->alloc(" << nbTuples << "," << nbComp << ");" << std::endl;
871 stream << varName << "->setName(\"" << getName() << "\");" << std::endl;
875 * Method that gives a quick overvien of \a this for python.
877 void DataArrayDouble::reprQuickOverview(std::ostream& stream) const
879 static const std::size_t MAX_NB_OF_BYTE_IN_REPR=300;
880 stream << "DataArrayDouble C++ instance at " << this << ". ";
883 std::size_t nbOfCompo(_info_on_compo.size());
886 mcIdType nbOfTuples(getNumberOfTuples());
887 stream << "Number of tuples : " << nbOfTuples << ". Number of components : " << nbOfCompo << "." << std::endl;
888 reprQuickOverviewData(stream,MAX_NB_OF_BYTE_IN_REPR);
891 stream << "Number of components : 0.";
894 stream << "*** No data allocated ****";
897 void DataArrayDouble::reprQuickOverviewData(std::ostream& stream, std::size_t maxNbOfByteInRepr) const
899 const double *data=begin();
900 mcIdType nbOfTuples(getNumberOfTuples());
901 std::size_t nbOfCompo(_info_on_compo.size());
902 std::ostringstream oss2; oss2 << "[";
904 std::string oss2Str(oss2.str());
905 bool isFinished=true;
906 for(mcIdType i=0;i<nbOfTuples && isFinished;i++)
911 for(std::size_t j=0;j<nbOfCompo;j++,data++)
914 if(j!=nbOfCompo-1) oss2 << ", ";
920 if(i!=nbOfTuples-1) oss2 << ", ";
921 std::string oss3Str(oss2.str());
922 if(oss3Str.length()<maxNbOfByteInRepr)
934 * Equivalent to DataArrayDouble::isEqual except that if false the reason of
937 * \param [in] other the instance to be compared with \a this
938 * \param [in] prec the precision to compare numeric data of the arrays.
939 * \param [out] reason In case of inequality returns the reason.
940 * \sa DataArrayDouble::isEqual
942 bool DataArrayDouble::isEqualIfNotWhy(const DataArrayDouble& other, double prec, std::string& reason) const
944 if(!areInfoEqualsIfNotWhy(other,reason))
946 return _mem.isEqual(other._mem,prec,reason);
950 * Checks if \a this and another DataArrayDouble are fully equal. For more info see
951 * \ref MEDCouplingArrayBasicsCompare.
952 * \param [in] other - an instance of DataArrayDouble to compare with \a this one.
953 * \param [in] prec - precision value to compare numeric data of the arrays.
954 * \return bool - \a true if the two arrays are equal, \a false else.
956 bool DataArrayDouble::isEqual(const DataArrayDouble& other, double prec) const
959 return isEqualIfNotWhy(other,prec,tmp);
963 * Checks if values of \a this and another DataArrayDouble are equal. For more info see
964 * \ref MEDCouplingArrayBasicsCompare.
965 * \param [in] other - an instance of DataArrayDouble to compare with \a this one.
966 * \param [in] prec - precision value to compare numeric data of the arrays.
967 * \return bool - \a true if the values of two arrays are equal, \a false else.
969 bool DataArrayDouble::isEqualWithoutConsideringStr(const DataArrayDouble& other, double prec) const
972 return _mem.isEqual(other._mem,prec,tmp);
976 * This method checks that all tuples in \a other are in \a this.
977 * If true, the output param \a tupleIds contains the tuples ids of \a this that correspond to tupes in \a this.
978 * 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.
980 * \param [in] other - the array having the same number of components than \a this.
981 * \param [out] tupleIds - the tuple ids containing the same number of tuples than \a other has.
982 * \sa DataArrayDouble::findCommonTuples
984 bool DataArrayDouble::areIncludedInMe(const DataArrayDouble *other, double prec, DataArrayIdType *&tupleIds) const
987 throw INTERP_KERNEL::Exception("DataArrayDouble::areIncludedInMe : input array is NULL !");
988 checkAllocated(); other->checkAllocated();
989 if(getNumberOfComponents()!=other->getNumberOfComponents())
990 throw INTERP_KERNEL::Exception("DataArrayDouble::areIncludedInMe : the number of components does not match !");
991 MCAuto<DataArrayDouble> a=DataArrayDouble::Aggregate(this,other);
992 DataArrayIdType *c=0,*ci=0;
993 a->findCommonTuples(prec,getNumberOfTuples(),c,ci);
994 MCAuto<DataArrayIdType> cSafe(c),ciSafe(ci);
995 mcIdType newNbOfTuples=-1;
996 MCAuto<DataArrayIdType> ids=DataArrayIdType::ConvertIndexArrayToO2N(a->getNumberOfTuples(),c->begin(),ci->begin(),ci->end(),newNbOfTuples);
997 MCAuto<DataArrayIdType> ret1=ids->selectByTupleIdSafeSlice(getNumberOfTuples(),a->getNumberOfTuples(),1);
998 tupleIds=ret1.retn();
999 return newNbOfTuples==getNumberOfTuples();
1003 * Searches for tuples coincident within \a prec tolerance. Each tuple is considered
1004 * as coordinates of a point in getNumberOfComponents()-dimensional space. The
1005 * distance separating two points is computed with the infinite norm.
1007 * Indices of coincident tuples are stored in output arrays.
1008 * A pair of arrays (\a comm, \a commIndex) is called "Surjective Format 2".
1010 * This method is typically used by MEDCouplingPointSet::findCommonNodes() and
1011 * MEDCouplingUMesh::mergeNodes().
1012 * \param [in] prec - minimal absolute distance between two tuples (infinite norm) at which they are
1013 * considered not coincident.
1014 * \param [in] limitTupleId - limit tuple id. If all tuples within a group of coincident
1015 * tuples have id strictly lower than \a limitTupleId then they are not returned.
1016 * \param [out] comm - the array holding ids (== indices) of coincident tuples.
1017 * \a comm->getNumberOfComponents() == 1.
1018 * \a comm->getNumberOfTuples() == \a commIndex->back().
1019 * \param [out] commIndex - the array dividing all indices stored in \a comm into
1020 * groups of (indices of) coincident tuples. Its every value is a tuple
1021 * index where a next group of tuples begins. For example the second
1022 * group of tuples in \a comm is described by following range of indices:
1023 * [ \a commIndex[1], \a commIndex[2] ). \a commIndex->getNumberOfTuples()-1
1024 * gives the number of groups of coincident tuples.
1025 * \throw If \a this is not allocated.
1026 * \throw If the number of components is not in [1,2,3,4].
1028 * \if ENABLE_EXAMPLES
1029 * \ref cpp_mcdataarraydouble_findcommontuples "Here is a C++ example".
1031 * \ref py_mcdataarraydouble_findcommontuples "Here is a Python example".
1033 * \sa DataArrayInt::ConvertIndexArrayToO2N(), DataArrayDouble::areIncludedInMe
1035 void DataArrayDouble::findCommonTuples(double prec, mcIdType limitTupleId, DataArrayIdType *&comm, DataArrayIdType *&commIndex) const
1038 std::size_t nbOfCompo=getNumberOfComponents();
1039 if ((nbOfCompo<1) || (nbOfCompo>4)) //test before work
1040 throw INTERP_KERNEL::Exception("DataArrayDouble::findCommonTuples : Unexpected spacedim of coords. Must be 1, 2, 3 or 4.");
1042 mcIdType nbOfTuples(getNumberOfTuples());
1044 MCAuto<DataArrayIdType> c(DataArrayIdType::New()),cI(DataArrayIdType::New()); c->alloc(0,1); cI->pushBackSilent(0);
1048 findCommonTuplesAlg<4>(begin(),nbOfTuples,limitTupleId,prec,c,cI);
1051 findCommonTuplesAlg<3>(begin(),nbOfTuples,limitTupleId,prec,c,cI);
1054 findCommonTuplesAlg<2>(begin(),nbOfTuples,limitTupleId,prec,c,cI);
1057 findCommonTuplesAlg<1>(begin(),nbOfTuples,limitTupleId,prec,c,cI);
1060 throw INTERP_KERNEL::Exception("DataArrayDouble::findCommonTuples : nb of components managed are 1,2,3 and 4 ! not implemented for other number of components !");
1063 commIndex=cI.retn();
1067 * This methods returns the minimal distance between the two set of points \a this and \a other.
1068 * So \a this and \a other have to have the same number of components. If not an INTERP_KERNEL::Exception will be thrown.
1069 * This method works only if number of components of \a this (equal to those of \a other) is in 1, 2 or 3.
1071 * \param [out] thisTupleId the tuple id in \a this corresponding to the returned minimal distance
1072 * \param [out] otherTupleId the tuple id in \a other corresponding to the returned minimal distance
1073 * \return the minimal distance between the two set of points \a this and \a other.
1074 * \sa DataArrayDouble::findClosestTupleId
1076 double DataArrayDouble::minimalDistanceTo(const DataArrayDouble *other, mcIdType& thisTupleId, mcIdType& otherTupleId) const
1078 MCAuto<DataArrayIdType> part1=findClosestTupleId(other);
1079 std::size_t nbOfCompo=getNumberOfComponents();
1080 mcIdType otherNbTuples=other->getNumberOfTuples();
1081 const double *thisPt(begin()),*otherPt(other->begin());
1082 const mcIdType *part1Pt(part1->begin());
1083 double ret=std::numeric_limits<double>::max();
1084 for(mcIdType i=0;i<otherNbTuples;i++,part1Pt++,otherPt+=nbOfCompo)
1087 for(std::size_t j=0;j<nbOfCompo;j++)
1088 tmp+=(otherPt[j]-thisPt[nbOfCompo*(*part1Pt)+j])*(otherPt[j]-thisPt[nbOfCompo*(*part1Pt)+j]);
1090 { ret=tmp; thisTupleId=*part1Pt; otherTupleId=i; }
1096 * This methods returns for each tuple in \a other which tuple in \a this is the closest.
1097 * So \a this and \a other have to have the same number of components. If not an INTERP_KERNEL::Exception will be thrown.
1098 * This method works only if number of components of \a this (equal to those of \a other) is in 1, 2 or 3.
1100 * \return a newly allocated (new object to be dealt by the caller) DataArrayInt having \c other->getNumberOfTuples() tuples and one components.
1101 * \sa DataArrayDouble::minimalDistanceTo
1103 DataArrayIdType *DataArrayDouble::findClosestTupleId(const DataArrayDouble *other) const
1106 throw INTERP_KERNEL::Exception("DataArrayDouble::findClosestTupleId : other instance is NULL !");
1107 checkAllocated(); other->checkAllocated();
1108 std::size_t nbOfCompo(getNumberOfComponents());
1109 if(nbOfCompo!=other->getNumberOfComponents())
1111 std::ostringstream oss; oss << "DataArrayDouble::findClosestTupleId : number of components in this is " << nbOfCompo;
1112 oss << ", whereas number of components in other is " << other->getNumberOfComponents() << "! Should be equal !";
1113 throw INTERP_KERNEL::Exception(oss.str().c_str());
1115 mcIdType nbOfTuples(other->getNumberOfTuples());
1116 mcIdType thisNbOfTuples(getNumberOfTuples());
1117 MCAuto<DataArrayIdType> ret=DataArrayIdType::New(); ret->alloc(nbOfTuples,1);
1119 getMinMaxPerComponent(bounds);
1124 double xDelta(fabs(bounds[1]-bounds[0])),yDelta(fabs(bounds[3]-bounds[2])),zDelta(fabs(bounds[5]-bounds[4]));
1125 double delta=std::max(xDelta,yDelta); delta=std::max(delta,zDelta);
1126 double characSize=pow((delta*delta*delta)/((double)thisNbOfTuples),1./3.);
1127 BBTreePts<3,mcIdType> myTree(begin(),0,0,getNumberOfTuples(),characSize*1e-12);
1128 FindClosestTupleIdAlg<3>(myTree,3.*characSize*characSize,other->begin(),nbOfTuples,begin(),thisNbOfTuples,ret->getPointer());
1133 double xDelta(fabs(bounds[1]-bounds[0])),yDelta(fabs(bounds[3]-bounds[2]));
1134 double delta=std::max(xDelta,yDelta);
1135 double characSize=sqrt(delta/(double)thisNbOfTuples);
1136 BBTreePts<2,mcIdType> myTree(begin(),0,0,getNumberOfTuples(),characSize*1e-12);
1137 FindClosestTupleIdAlg<2>(myTree,2.*characSize*characSize,other->begin(),nbOfTuples,begin(),thisNbOfTuples,ret->getPointer());
1142 double characSize=fabs(bounds[1]-bounds[0])/FromIdType<double>(thisNbOfTuples);
1143 BBTreePts<1,mcIdType> myTree(begin(),0,0,getNumberOfTuples(),characSize*1e-12);
1144 FindClosestTupleIdAlg<1>(myTree,1.*characSize*characSize,other->begin(),nbOfTuples,begin(),thisNbOfTuples,ret->getPointer());
1148 throw INTERP_KERNEL::Exception("Unexpected spacedim of coords for findClosestTupleId. Must be 1, 2 or 3.");
1154 * This method expects that \a this and \a otherBBoxFrmt arrays are bounding box arrays ( as the output of MEDCouplingPointSet::getBoundingBoxForBBTree method ).
1155 * 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
1156 * how many bounding boxes in \a otherBBoxFrmt.
1157 * So, this method expects that \a this and \a otherBBoxFrmt have the same number of components.
1159 * \param [in] otherBBoxFrmt - It is an array .
1160 * \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.
1161 * \sa MEDCouplingPointSet::getBoundingBoxForBBTree
1162 * \throw If \a this and \a otherBBoxFrmt have not the same number of components.
1163 * \throw If \a this and \a otherBBoxFrmt number of components is not even (BBox format).
1165 DataArrayIdType *DataArrayDouble::computeNbOfInteractionsWith(const DataArrayDouble *otherBBoxFrmt, double eps) const
1168 throw INTERP_KERNEL::Exception("DataArrayDouble::computeNbOfInteractionsWith : input array is NULL !");
1169 if(!isAllocated() || !otherBBoxFrmt->isAllocated())
1170 throw INTERP_KERNEL::Exception("DataArrayDouble::computeNbOfInteractionsWith : this and input array must be allocated !");
1171 std::size_t nbOfComp(getNumberOfComponents());
1172 mcIdType nbOfTuples(getNumberOfTuples());
1173 if(nbOfComp!=otherBBoxFrmt->getNumberOfComponents())
1175 std::ostringstream oss; oss << "DataArrayDouble::computeNbOfInteractionsWith : this number of components (" << nbOfComp << ") must be equal to the number of components of input array (" << otherBBoxFrmt->getNumberOfComponents() << ") !";
1176 throw INTERP_KERNEL::Exception(oss.str().c_str());
1180 std::ostringstream oss; oss << "DataArrayDouble::computeNbOfInteractionsWith : Number of components (" << nbOfComp << ") is not even ! It should be to be compatible with bbox format !";
1181 throw INTERP_KERNEL::Exception(oss.str().c_str());
1183 MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(nbOfTuples,1);
1184 const double *thisBBPtr(begin());
1185 mcIdType *retPtr(ret->getPointer());
1190 BBTree<3,mcIdType> bbt(otherBBoxFrmt->begin(),0,0,otherBBoxFrmt->getNumberOfTuples(),eps);
1191 for(mcIdType i=0;i<nbOfTuples;i++,retPtr++,thisBBPtr+=nbOfComp)
1192 *retPtr=bbt.getNbOfIntersectingElems(thisBBPtr);
1197 BBTree<2,mcIdType> bbt(otherBBoxFrmt->begin(),0,0,otherBBoxFrmt->getNumberOfTuples(),eps);
1198 for(mcIdType i=0;i<nbOfTuples;i++,retPtr++,thisBBPtr+=nbOfComp)
1199 *retPtr=bbt.getNbOfIntersectingElems(thisBBPtr);
1204 BBTree<1,mcIdType> bbt(otherBBoxFrmt->begin(),0,0,otherBBoxFrmt->getNumberOfTuples(),eps);
1205 for(mcIdType i=0;i<nbOfTuples;i++,retPtr++,thisBBPtr+=nbOfComp)
1206 *retPtr=bbt.getNbOfIntersectingElems(thisBBPtr);
1210 throw INTERP_KERNEL::Exception("DataArrayDouble::computeNbOfInteractionsWith : space dimension supported are [1,2,3] !");
1217 * Returns a copy of \a this array by excluding coincident tuples. Each tuple is
1218 * considered as coordinates of a point in getNumberOfComponents()-dimensional
1219 * space. The distance between tuples is computed using norm2. If several tuples are
1220 * not far each from other than \a prec, only one of them remains in the result
1221 * array. The order of tuples in the result array is same as in \a this one except
1222 * that coincident tuples are excluded.
1223 * \param [in] prec - minimal absolute distance between two tuples at which they are
1224 * considered not coincident.
1225 * \param [in] limitTupleId - limit tuple id. If all tuples within a group of coincident
1226 * tuples have id strictly lower than \a limitTupleId then they are not excluded.
1227 * \return DataArrayDouble * - the new instance of DataArrayDouble that the caller
1228 * is to delete using decrRef() as it is no more needed.
1229 * \throw If \a this is not allocated.
1230 * \throw If the number of components is not in [1,2,3,4].
1232 * \if ENABLE_EXAMPLES
1233 * \ref py_mcdataarraydouble_getdifferentvalues "Here is a Python example".
1236 DataArrayDouble *DataArrayDouble::getDifferentValues(double prec, mcIdType limitTupleId) const
1239 DataArrayIdType *c0=0,*cI0=0;
1240 findCommonTuples(prec,limitTupleId,c0,cI0);
1241 MCAuto<DataArrayIdType> c(c0),cI(cI0);
1242 mcIdType newNbOfTuples=-1;
1243 MCAuto<DataArrayIdType> o2n=DataArrayIdType::ConvertIndexArrayToO2N(getNumberOfTuples(),c0->begin(),cI0->begin(),cI0->end(),newNbOfTuples);
1244 return renumberAndReduce(o2n->getConstPointer(),newNbOfTuples);
1248 * Copy all components in a specified order from another DataArrayDouble.
1249 * Both numerical and textual data is copied. The number of tuples in \a this and
1250 * the other array can be different.
1251 * \param [in] a - the array to copy data from.
1252 * \param [in] compoIds - sequence of zero based indices of components, data of which is
1254 * \throw If \a a is NULL.
1255 * \throw If \a compoIds.size() != \a a->getNumberOfComponents().
1256 * \throw If \a compoIds[i] < 0 or \a compoIds[i] > \a this->getNumberOfComponents().
1258 * \if ENABLE_EXAMPLES
1259 * \ref py_mcdataarraydouble_setselectedcomponents "Here is a Python example".
1262 void DataArrayDouble::setSelectedComponents(const DataArrayDouble *a, const std::vector<std::size_t>& compoIds)
1265 throw INTERP_KERNEL::Exception("DataArrayDouble::setSelectedComponents : input DataArrayDouble is NULL !");
1267 copyPartOfStringInfoFrom2(compoIds,*a);
1268 std::size_t partOfCompoSz=compoIds.size();
1269 std::size_t nbOfCompo=getNumberOfComponents();
1270 mcIdType nbOfTuples=std::min(getNumberOfTuples(),a->getNumberOfTuples());
1271 const double *ac=a->getConstPointer();
1272 double *nc=getPointer();
1273 for(mcIdType i=0;i<nbOfTuples;i++)
1274 for(std::size_t j=0;j<partOfCompoSz;j++,ac++)
1275 nc[nbOfCompo*i+compoIds[j]]=*ac;
1278 * Checks if 0.0 value is present in \a this array. If it is the case, an exception
1280 * \throw If zero is found in \a this array.
1282 void DataArrayDouble::checkNoNullValues() const
1284 const double *tmp=getConstPointer();
1285 mcIdType nbOfElems=getNbOfElems();
1286 const double *where=std::find(tmp,tmp+nbOfElems,0.);
1287 if(where!=tmp+nbOfElems)
1288 throw INTERP_KERNEL::Exception("A value 0.0 have been detected !");
1292 * Computes minimal and maximal value in each component. An output array is filled
1293 * with \c 2 * \a this->getNumberOfComponents() values, so the caller is to allocate
1294 * enough memory before calling this method.
1295 * \param [out] bounds - array of size at least 2 *\a this->getNumberOfComponents().
1296 * It is filled as follows:<br>
1297 * \a bounds[0] = \c min_of_component_0 <br>
1298 * \a bounds[1] = \c max_of_component_0 <br>
1299 * \a bounds[2] = \c min_of_component_1 <br>
1300 * \a bounds[3] = \c max_of_component_1 <br>
1303 void DataArrayDouble::getMinMaxPerComponent(double *bounds) const
1306 std::size_t dim=getNumberOfComponents();
1307 for (std::size_t idim=0; idim<dim; idim++)
1309 bounds[idim*2]=std::numeric_limits<double>::max();
1310 bounds[idim*2+1]=-std::numeric_limits<double>::max();
1312 const double *ptr=getConstPointer();
1313 mcIdType nbOfTuples=getNumberOfTuples();
1314 for(mcIdType i=0;i<nbOfTuples;i++)
1316 for(std::size_t idim=0;idim<dim;idim++)
1318 if(bounds[idim*2]>ptr[i*dim+idim])
1320 bounds[idim*2]=ptr[i*dim+idim];
1322 if(bounds[idim*2+1]<ptr[i*dim+idim])
1324 bounds[idim*2+1]=ptr[i*dim+idim];
1331 * This method retrieves a newly allocated DataArrayDouble instance having same number of tuples than \a this and twice number of components than \a this
1332 * to store both the min and max per component of each tuples.
1333 * \param [in] epsilon the width of the bbox (identical in each direction) - 0.0 by default
1335 * \return a newly created DataArrayDouble instance having \c this->getNumberOfTuples() tuples and 2 * \c this->getNumberOfComponent() components
1337 * \throw If \a this is not allocated yet.
1339 DataArrayDouble *DataArrayDouble::computeBBoxPerTuple(double epsilon) const
1342 const double *dataPtr=getConstPointer();
1343 std::size_t nbOfCompo=getNumberOfComponents();
1344 mcIdType nbTuples=getNumberOfTuples();
1345 MCAuto<DataArrayDouble> bbox=DataArrayDouble::New();
1346 bbox->alloc(nbTuples,2*nbOfCompo);
1347 double *bboxPtr=bbox->getPointer();
1348 for(mcIdType i=0;i<nbTuples;i++)
1350 for(std::size_t j=0;j<nbOfCompo;j++)
1352 bboxPtr[2*nbOfCompo*i+2*j]=dataPtr[nbOfCompo*i+j]-epsilon;
1353 bboxPtr[2*nbOfCompo*i+2*j+1]=dataPtr[nbOfCompo*i+j]+epsilon;
1360 * For each tuples **t** in \a other, this method retrieves tuples in \a this that are equal to **t**.
1361 * Two tuples are considered equal if the euclidian distance between the two tuples is lower than \a eps.
1363 * \param [in] other a DataArrayDouble having same number of components than \a this.
1364 * \param [in] eps absolute precision representing distance (using infinite norm) between 2 tuples behind which 2 tuples are considered equal.
1365 * \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.
1366 * \a cI allows to extract information in \a c.
1367 * \param [out] cI is an indirection array that allows to extract the data contained in \a c.
1369 * \throw In case of:
1370 * - \a this is not allocated
1371 * - \a other is not allocated or null
1372 * - \a this and \a other do not have the same number of components
1373 * - if number of components of \a this is not in [1,2,3]
1375 * \sa MEDCouplingPointSet::getNodeIdsNearPoints, DataArrayDouble::getDifferentValues
1377 void DataArrayDouble::computeTupleIdsNearTuples(const DataArrayDouble *other, double eps, DataArrayIdType *& c, DataArrayIdType *& cI) const
1380 throw INTERP_KERNEL::Exception("DataArrayDouble::computeTupleIdsNearTuples : input pointer other is null !");
1382 other->checkAllocated();
1383 std::size_t nbOfCompo=getNumberOfComponents();
1384 std::size_t otherNbOfCompo=other->getNumberOfComponents();
1385 if(nbOfCompo!=otherNbOfCompo)
1386 throw INTERP_KERNEL::Exception("DataArrayDouble::computeTupleIdsNearTuples : number of components should be equal between this and other !");
1387 mcIdType nbOfTuplesOther=other->getNumberOfTuples();
1388 mcIdType nbOfTuples=getNumberOfTuples();
1389 MCAuto<DataArrayIdType> cArr(DataArrayIdType::New()),cIArr(DataArrayIdType::New()); cArr->alloc(0,1); cIArr->pushBackSilent(0);
1394 BBTreePts<3,mcIdType> myTree(begin(),0,0,nbOfTuples,eps);
1395 FindTupleIdsNearTuplesAlg<3>(myTree,other->getConstPointer(),nbOfTuplesOther,eps,cArr,cIArr);
1400 BBTreePts<2,mcIdType> myTree(begin(),0,0,nbOfTuples,eps);
1401 FindTupleIdsNearTuplesAlg<2>(myTree,other->getConstPointer(),nbOfTuplesOther,eps,cArr,cIArr);
1406 BBTreePts<1,mcIdType> myTree(begin(),0,0,nbOfTuples,eps);
1407 FindTupleIdsNearTuplesAlg<1>(myTree,other->getConstPointer(),nbOfTuplesOther,eps,cArr,cIArr);
1411 throw INTERP_KERNEL::Exception("Unexpected spacedim of coords for computeTupleIdsNearTuples. Must be 1, 2 or 3.");
1413 c=cArr.retn(); cI=cIArr.retn();
1417 * 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
1418 * around origin of 'radius' 1.
1420 * \param [in] eps absolute epsilon. under that value of delta between max and min no scale is performed.
1422 void DataArrayDouble::recenterForMaxPrecision(double eps)
1425 std::size_t dim=getNumberOfComponents();
1426 std::vector<double> bounds(2*dim);
1427 getMinMaxPerComponent(&bounds[0]);
1428 for(std::size_t i=0;i<dim;i++)
1430 double delta=bounds[2*i+1]-bounds[2*i];
1431 double offset=(bounds[2*i]+bounds[2*i+1])/2.;
1433 applyLin(1./delta,-offset/delta,i);
1435 applyLin(1.,-offset,i);
1440 * Returns the maximal value and all its locations within \a this one-dimensional array.
1441 * \param [out] tupleIds - a new instance of DataArrayInt containing indices of
1442 * tuples holding the maximal value. The caller is to delete it using
1443 * decrRef() as it is no more needed.
1444 * \return double - the maximal value among all values of \a this array.
1445 * \throw If \a this->getNumberOfComponents() != 1
1446 * \throw If \a this->getNumberOfTuples() < 1
1448 double DataArrayDouble::getMaxValue2(DataArrayIdType*& tupleIds) const
1452 double ret=getMaxValue(tmp);
1453 tupleIds=findIdsInRange(ret,ret);
1458 * Returns the minimal value and all its locations within \a this one-dimensional array.
1459 * \param [out] tupleIds - a new instance of DataArrayInt containing indices of
1460 * tuples holding the minimal value. The caller is to delete it using
1461 * decrRef() as it is no more needed.
1462 * \return double - the minimal value among all values of \a this array.
1463 * \throw If \a this->getNumberOfComponents() != 1
1464 * \throw If \a this->getNumberOfTuples() < 1
1466 double DataArrayDouble::getMinValue2(DataArrayIdType*& tupleIds) const
1470 double ret=getMinValue(tmp);
1471 tupleIds=findIdsInRange(ret,ret);
1476 * 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.
1477 * This method only works for single component array.
1479 * \return a value in [ 0, \c this->getNumberOfTuples() )
1481 * \throw If \a this is not allocated
1484 mcIdType DataArrayDouble::count(double value, double eps) const
1488 if(getNumberOfComponents()!=1)
1489 throw INTERP_KERNEL::Exception("DataArrayDouble::count : must be applied on DataArrayDouble with only one component, you can call 'rearrange' method before !");
1490 const double *vals=begin();
1491 mcIdType nbOfTuples=getNumberOfTuples();
1492 for(mcIdType i=0;i<nbOfTuples;i++,vals++)
1493 if(fabs(*vals-value)<=eps)
1499 * Returns the average value of \a this one-dimensional array.
1500 * \return double - the average value over all values of \a this array.
1501 * \throw If \a this->getNumberOfComponents() != 1
1502 * \throw If \a this->getNumberOfTuples() < 1
1504 double DataArrayDouble::getAverageValue() const
1506 if(getNumberOfComponents()!=1)
1507 throw INTERP_KERNEL::Exception("DataArrayDouble::getAverageValue : must be applied on DataArrayDouble with only one component, you can call 'rearrange' method before !");
1508 mcIdType nbOfTuples(getNumberOfTuples());
1510 throw INTERP_KERNEL::Exception("DataArrayDouble::getAverageValue : array exists but number of tuples must be > 0 !");
1511 const double *vals=getConstPointer();
1512 double ret=std::accumulate(vals,vals+nbOfTuples,0.);
1513 return ret/FromIdType<double>(nbOfTuples);
1517 * Returns the Euclidean norm of the vector defined by \a this array.
1518 * \return double - the value of the Euclidean norm, i.e.
1519 * the square root of the inner product of vector.
1520 * \throw If \a this is not allocated.
1522 double DataArrayDouble::norm2() const
1526 std::size_t nbOfElems=getNbOfElems();
1527 const double *pt=getConstPointer();
1528 for(std::size_t i=0;i<nbOfElems;i++,pt++)
1534 * Returns the maximum norm of the vector defined by \a this array.
1535 * This method works even if the number of components is different from one.
1536 * If the number of elements in \a this is 0, -1. is returned.
1537 * \return double - the value of the maximum norm, i.e.
1538 * the maximal absolute value among values of \a this array (whatever its number of components).
1539 * \throw If \a this is not allocated.
1541 double DataArrayDouble::normMax() const
1545 std::size_t nbOfElems(getNbOfElems());
1546 const double *pt(getConstPointer());
1547 for(std::size_t i=0;i<nbOfElems;i++,pt++)
1549 double val(std::abs(*pt));
1557 * Returns the maximum norm of for each component of \a this array.
1558 * If the number of elements in \a this is 0, -1. is returned.
1559 * \param [out] res - pointer to an array of result values, of size at least \a
1560 * this->getNumberOfComponents(), that is to be allocated by the caller.
1561 * \throw If \a this is not allocated.
1563 void DataArrayDouble::normMaxPerComponent(double * res) const
1566 mcIdType nbOfTuples(getNumberOfTuples());
1567 std::size_t nbOfCompos(getNumberOfComponents());
1568 std::fill(res, res+nbOfCompos, -1.0);
1569 const double *pt(getConstPointer());
1570 for(mcIdType i=0;i<nbOfTuples;i++)
1571 for (std::size_t j=0; j<nbOfCompos; j++, pt++)
1573 double val(std::abs(*pt));
1581 * Returns the minimum norm (absolute value) of the vector defined by \a this array.
1582 * This method works even if the number of components is different from one.
1583 * If the number of elements in \a this is 0, std::numeric_limits<double>::max() is returned.
1584 * \return double - the value of the minimum norm, i.e.
1585 * the minimal absolute value among values of \a this array (whatever its number of components).
1586 * \throw If \a this is not allocated.
1588 double DataArrayDouble::normMin() const
1591 double ret(std::numeric_limits<double>::max());
1592 std::size_t nbOfElems(getNbOfElems());
1593 const double *pt(getConstPointer());
1594 for(std::size_t i=0;i<nbOfElems;i++,pt++)
1596 double val(std::abs(*pt));
1604 * Accumulates values of each component of \a this array.
1605 * \param [out] res - an array of length \a this->getNumberOfComponents(), allocated
1606 * by the caller, that is filled by this method with sum value for each
1608 * \throw If \a this is not allocated.
1610 void DataArrayDouble::accumulate(double *res) const
1613 const double *ptr=getConstPointer();
1614 mcIdType nbTuple(getNumberOfTuples());
1615 std::size_t nbComps(getNumberOfComponents());
1616 std::fill(res,res+nbComps,0.);
1617 for(mcIdType i=0;i<nbTuple;i++)
1618 std::transform(ptr+i*nbComps,ptr+(i+1)*nbComps,res,res,std::plus<double>());
1622 * This method returns the min distance from an external tuple defined by [ \a tupleBg , \a tupleEnd ) to \a this and
1623 * the first tuple in \a this that matches the returned distance. If there is no tuples in \a this an exception will be thrown.
1626 * \a this is expected to be allocated and expected to have a number of components equal to the distance from \a tupleBg to
1627 * \a tupleEnd. If not an exception will be thrown.
1629 * \param [in] tupleBg start pointer (included) of input external tuple
1630 * \param [in] tupleEnd end pointer (not included) of input external tuple
1631 * \param [out] tupleId the tuple id in \a this that matches the min of distance between \a this and input external tuple
1632 * \return the min distance.
1633 * \sa MEDCouplingUMesh::distanceToPoint
1635 double DataArrayDouble::distanceToTuple(const double *tupleBg, const double *tupleEnd, mcIdType& tupleId) const
1638 mcIdType nbTuple(getNumberOfTuples());
1639 std::size_t nbComps(getNumberOfComponents());
1640 if(nbComps!=(std::size_t)std::distance(tupleBg,tupleEnd))
1641 { 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()); }
1643 throw INTERP_KERNEL::Exception("DataArrayDouble::distanceToTuple : no tuple in this ! No distance to compute !");
1644 double ret0=std::numeric_limits<double>::max();
1646 const double *work=getConstPointer();
1647 for(mcIdType i=0;i<nbTuple;i++)
1650 for(std::size_t j=0;j<nbComps;j++,work++)
1651 val+=(*work-tupleBg[j])*((*work-tupleBg[j]));
1655 { ret0=val; tupleId=i; }
1661 * Accumulate values of the given component of \a this array.
1662 * \param [in] compId - the index of the component of interest.
1663 * \return double - a sum value of \a compId-th component.
1664 * \throw If \a this is not allocated.
1665 * \throw If \a the condition ( 0 <= \a compId < \a this->getNumberOfComponents() ) is
1668 double DataArrayDouble::accumulate(std::size_t compId) const
1671 const double *ptr=getConstPointer();
1672 mcIdType nbTuple(getNumberOfTuples());
1673 std::size_t nbComps(getNumberOfComponents());
1675 throw INTERP_KERNEL::Exception("DataArrayDouble::accumulate : Invalid compId specified : No such nb of components !");
1677 for(mcIdType i=0;i<nbTuple;i++)
1678 ret+=ptr[i*nbComps+compId];
1683 * This method accumulate using addition tuples in \a this using input index array [ \a bgOfIndex, \a endOfIndex ).
1684 * The returned array will have same number of components than \a this and number of tuples equal to
1685 * \c std::distance(bgOfIndex,endOfIndex) \b minus \b one.
1687 * The input index array is expected to be ascendingly sorted in which the all referenced ids should be in [0, \c this->getNumberOfTuples).
1688 * 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.
1690 * \param [in] bgOfIndex - begin (included) of the input index array.
1691 * \param [in] endOfIndex - end (excluded) of the input index array.
1692 * \return DataArrayDouble * - the new instance having the same number of components than \a this.
1694 * \throw If bgOfIndex or end is NULL.
1695 * \throw If input index array is not ascendingly sorted.
1696 * \throw If there is an id in [ \a bgOfIndex, \a endOfIndex ) not in [0, \c this->getNumberOfTuples).
1697 * \throw If std::distance(bgOfIndex,endOfIndex)==0.
1699 DataArrayDouble *DataArrayDouble::accumulatePerChunck(const mcIdType *bgOfIndex, const mcIdType *endOfIndex) const
1701 if(!bgOfIndex || !endOfIndex)
1702 throw INTERP_KERNEL::Exception("DataArrayDouble::accumulatePerChunck : input pointer NULL !");
1704 std::size_t nbCompo(getNumberOfComponents());
1705 mcIdType nbOfTuples(getNumberOfTuples());
1706 std::size_t sz=std::distance(bgOfIndex,endOfIndex);
1708 throw INTERP_KERNEL::Exception("DataArrayDouble::accumulatePerChunck : invalid size of input index array !");
1710 MCAuto<DataArrayDouble> ret=DataArrayDouble::New(); ret->alloc(sz,nbCompo);
1711 const mcIdType *w=bgOfIndex;
1712 if(*w<0 || *w>=nbOfTuples)
1713 throw INTERP_KERNEL::Exception("DataArrayDouble::accumulatePerChunck : The first element of the input index not in [0,nbOfTuples) !");
1714 const double *srcPt=begin()+(*w)*nbCompo;
1715 double *tmp=ret->getPointer();
1716 for(std::size_t i=0;i<sz;i++,tmp+=nbCompo,w++)
1718 std::fill(tmp,tmp+nbCompo,0.);
1721 for(mcIdType j=w[0];j<w[1];j++,srcPt+=nbCompo)
1723 if(j>=0 && j<nbOfTuples)
1724 std::transform(srcPt,srcPt+nbCompo,tmp,tmp,std::plus<double>());
1727 std::ostringstream oss; oss << "DataArrayDouble::accumulatePerChunck : At rank #" << i << " the input index array points to id " << j << " should be in [0," << nbOfTuples << ") !";
1728 throw INTERP_KERNEL::Exception(oss.str().c_str());
1734 std::ostringstream oss; oss << "DataArrayDouble::accumulatePerChunck : At rank #" << i << " the input index array is not in ascendingly sorted.";
1735 throw INTERP_KERNEL::Exception(oss.str().c_str());
1738 ret->copyStringInfoFrom(*this);
1743 * 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.
1744 * This method expects that \a this as only one component. The returned array will have \a this->getNumberOfTuples()+1 tuple with also one component.
1745 * The ith element of returned array is equal to the sum of elements in \a this with rank strictly lower than i.
1747 * \return DataArrayDouble - A newly built array containing cum sum of \a this.
1749 MCAuto<DataArrayDouble> DataArrayDouble::cumSum() const
1752 checkNbOfComps(1,"DataArrayDouble::cumSum : this is expected to be single component");
1753 mcIdType nbOfTuple(getNumberOfTuples());
1754 MCAuto<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(nbOfTuple+1,1);
1755 double *ptr(ret->getPointer());
1757 const double *thisPtr(begin());
1758 for(mcIdType i=0;i<nbOfTuple;i++)
1759 ptr[i+1]=ptr[i]+thisPtr[i];
1764 * Converts each 2D point defined by the tuple of \a this array from the Polar to the
1765 * Cartesian coordinate system. The two components of the tuple of \a this array are
1766 * considered to contain (1) radius and (2) angle of the point in the Polar CS.
1767 * \return DataArrayDouble * - the new instance of DataArrayDouble, whose each tuple
1768 * contains X and Y coordinates of the point in the Cartesian CS. The caller
1769 * is to delete this array using decrRef() as it is no more needed. The array
1770 * does not contain any textual info on components.
1771 * \throw If \a this->getNumberOfComponents() != 2.
1772 * \sa fromCartToPolar
1774 DataArrayDouble *DataArrayDouble::fromPolarToCart() const
1777 std::size_t nbOfComp(getNumberOfComponents());
1779 throw INTERP_KERNEL::Exception("DataArrayDouble::fromPolarToCart : must be an array with exactly 2 components !");
1780 mcIdType nbOfTuple(getNumberOfTuples());
1781 DataArrayDouble *ret(DataArrayDouble::New());
1782 ret->alloc(nbOfTuple,2);
1783 double *w(ret->getPointer());
1784 const double *wIn(getConstPointer());
1785 for(mcIdType i=0;i<nbOfTuple;i++,w+=2,wIn+=2)
1787 w[0]=wIn[0]*cos(wIn[1]);
1788 w[1]=wIn[0]*sin(wIn[1]);
1794 * Converts each 3D point defined by the tuple of \a this array from the Cylindrical to
1795 * the Cartesian coordinate system. The three components of the tuple of \a this array
1796 * are considered to contain (1) radius, (2) azimuth and (3) altitude of the point in
1797 * the Cylindrical CS.
1798 * \return DataArrayDouble * - the new instance of DataArrayDouble, whose each tuple
1799 * contains X, Y and Z coordinates of the point in the Cartesian CS. The info
1800 * on the third component is copied from \a this array. The caller
1801 * is to delete this array using decrRef() as it is no more needed.
1802 * \throw If \a this->getNumberOfComponents() != 3.
1805 DataArrayDouble *DataArrayDouble::fromCylToCart() const
1808 std::size_t nbOfComp(getNumberOfComponents());
1810 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCylToCart : must be an array with exactly 3 components !");
1811 mcIdType nbOfTuple(getNumberOfTuples());
1812 DataArrayDouble *ret(DataArrayDouble::New());
1813 ret->alloc(getNumberOfTuples(),3);
1814 double *w(ret->getPointer());
1815 const double *wIn(getConstPointer());
1816 for(mcIdType i=0;i<nbOfTuple;i++,w+=3,wIn+=3)
1818 w[0]=wIn[0]*cos(wIn[1]);
1819 w[1]=wIn[0]*sin(wIn[1]);
1822 ret->setInfoOnComponent(2,getInfoOnComponent(2));
1827 * Converts each 3D point defined by the tuple of \a this array from the Spherical to
1828 * the Cartesian coordinate system. The three components of the tuple of \a this array
1829 * are considered to contain (1) radius, (2) polar angle and (3) azimuthal angle of the
1830 * point in the Cylindrical CS.
1831 * \return DataArrayDouble * - the new instance of DataArrayDouble, whose each tuple
1832 * contains X, Y and Z coordinates of the point in the Cartesian CS. The info
1833 * on the third component is copied from \a this array. The caller
1834 * is to delete this array using decrRef() as it is no more needed.
1835 * \throw If \a this->getNumberOfComponents() != 3.
1836 * \sa fromCartToSpher
1838 DataArrayDouble *DataArrayDouble::fromSpherToCart() const
1841 std::size_t nbOfComp(getNumberOfComponents());
1843 throw INTERP_KERNEL::Exception("DataArrayDouble::fromSpherToCart : must be an array with exactly 3 components !");
1844 mcIdType nbOfTuple(getNumberOfTuples());
1845 DataArrayDouble *ret(DataArrayDouble::New());
1846 ret->alloc(getNumberOfTuples(),3);
1847 double *w(ret->getPointer());
1848 const double *wIn(getConstPointer());
1849 for(mcIdType i=0;i<nbOfTuple;i++,w+=3,wIn+=3)
1851 w[0]=wIn[0]*cos(wIn[2])*sin(wIn[1]);
1852 w[1]=wIn[0]*sin(wIn[2])*sin(wIn[1]);
1853 w[2]=wIn[0]*cos(wIn[1]);
1859 * 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.
1860 * 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.
1861 * If \a at equals to AX_CYL the returned array will be the result of operation cylindric to cartesian of \a this...
1863 * \param [in] atOfThis - The axis type of \a this.
1864 * \return DataArrayDouble * - the new instance of DataArrayDouble (that must be dealed by caller) containing the result of the cartesianizification of \a this.
1866 DataArrayDouble *DataArrayDouble::cartesianize(MEDCouplingAxisType atOfThis) const
1869 std::size_t nbOfComp(getNumberOfComponents());
1870 MCAuto<DataArrayDouble> ret;
1878 ret=fromCylToCart();
1883 ret=fromPolarToCart();
1887 throw INTERP_KERNEL::Exception("DataArrayDouble::cartesianize : For AX_CYL, number of components must be in [2,3] !");
1891 ret=fromSpherToCart();
1896 ret=fromPolarToCart();
1900 throw INTERP_KERNEL::Exception("DataArrayDouble::cartesianize : For AX_CYL, number of components must be in [2,3] !");
1902 throw INTERP_KERNEL::Exception("DataArrayDouble::cartesianize : not recognized axis type ! Only AX_CART, AX_CYL and AX_SPHER supported !");
1904 ret->copyStringInfoFrom(*this);
1909 * This method returns a newly created array to be deallocated that contains the result of conversion from cartesian to polar.
1910 * This method expects that \a this has exactly 2 components.
1911 * \sa fromPolarToCart
1913 DataArrayDouble *DataArrayDouble::fromCartToPolar() const
1915 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
1917 std::size_t nbOfComp(getNumberOfComponents());
1918 mcIdType nbTuples(getNumberOfTuples());
1920 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToPolar : must be an array with exactly 2 components !");
1921 ret->alloc(nbTuples,2);
1922 double *retPtr(ret->getPointer());
1923 const double *ptr(begin());
1924 for(mcIdType i=0;i<nbTuples;i++,ptr+=2,retPtr+=2)
1926 retPtr[0]=sqrt(ptr[0]*ptr[0]+ptr[1]*ptr[1]);
1927 retPtr[1]=atan2(ptr[1],ptr[0]);
1933 * This method returns a newly created array to be deallocated that contains the result of conversion from cartesian to cylindrical.
1934 * This method expects that \a this has exactly 3 components.
1937 DataArrayDouble *DataArrayDouble::fromCartToCyl() const
1939 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
1941 std::size_t nbOfComp(getNumberOfComponents());
1942 mcIdType nbTuples(getNumberOfTuples());
1944 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCyl : must be an array with exactly 3 components !");
1945 ret->alloc(nbTuples,3);
1946 double *retPtr(ret->getPointer());
1947 const double *ptr(begin());
1948 for(mcIdType i=0;i<nbTuples;i++,ptr+=3,retPtr+=3)
1950 retPtr[0]=sqrt(ptr[0]*ptr[0]+ptr[1]*ptr[1]);
1951 retPtr[1]=atan2(ptr[1],ptr[0]);
1958 * This method returns a newly created array to be deallocated that contains the result of conversion from cartesian to spherical coordinates.
1959 * \sa fromSpherToCart
1961 DataArrayDouble *DataArrayDouble::fromCartToSpher() const
1963 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
1965 std::size_t nbOfComp(getNumberOfComponents());
1966 mcIdType nbTuples(getNumberOfTuples());
1968 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToSpher : must be an array with exactly 3 components !");
1969 ret->alloc(nbTuples,3);
1970 double *retPtr(ret->getPointer());
1971 const double *ptr(begin());
1972 for(mcIdType i=0;i<nbTuples;i++,ptr+=3,retPtr+=3)
1974 retPtr[0]=sqrt(ptr[0]*ptr[0]+ptr[1]*ptr[1]+ptr[2]*ptr[2]);
1975 retPtr[1]=acos(ptr[2]/retPtr[0]);
1976 retPtr[2]=atan2(ptr[1],ptr[0]);
1982 * 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.
1983 * This method expects that \a this has exactly 3 components.
1984 * \sa MEDCouplingFieldDouble::computeVectorFieldCyl
1986 DataArrayDouble *DataArrayDouble::fromCartToCylGiven(const DataArrayDouble *coords, const double center[3], const double vect[3]) const
1989 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCylGiven : input coords are NULL !");
1990 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
1991 checkAllocated(); coords->checkAllocated();
1992 std::size_t nbOfComp(getNumberOfComponents());
1993 mcIdType nbTuples(getNumberOfTuples());
1995 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCylGiven : must be an array with exactly 3 components !");
1996 if(coords->getNumberOfComponents()!=3)
1997 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCylGiven : coords array must have exactly 3 components !");
1998 if(coords->getNumberOfTuples()!=nbTuples)
1999 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCylGiven : coords array must have the same number of tuples !");
2000 ret->alloc(nbTuples,nbOfComp);
2001 double magOfVect(sqrt(vect[0]*vect[0]+vect[1]*vect[1]+vect[2]*vect[2]));
2003 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCylGiven : magnitude of vect is too low !");
2004 double Ur[3],Uteta[3],Uz[3],*retPtr(ret->getPointer());
2005 const double *coo(coords->begin()),*vectField(begin());
2006 std::transform(vect,vect+3,Uz,std::bind2nd(std::multiplies<double>(),1./magOfVect));
2007 for(mcIdType i=0;i<nbTuples;i++,vectField+=3,retPtr+=3,coo+=3)
2009 std::transform(coo,coo+3,center,Ur,std::minus<double>());
2010 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];
2011 double magOfTeta(sqrt(Uteta[0]*Uteta[0]+Uteta[1]*Uteta[1]+Uteta[2]*Uteta[2]));
2012 std::transform(Uteta,Uteta+3,Uteta,std::bind2nd(std::multiplies<double>(),1./magOfTeta));
2013 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];
2014 retPtr[0]=Ur[0]*vectField[0]+Ur[1]*vectField[1]+Ur[2]*vectField[2];
2015 retPtr[1]=Uteta[0]*vectField[0]+Uteta[1]*vectField[1]+Uteta[2]*vectField[2];
2016 retPtr[2]=Uz[0]*vectField[0]+Uz[1]*vectField[1]+Uz[2]*vectField[2];
2018 ret->copyStringInfoFrom(*this);
2023 * Computes the doubly contracted product of every tensor defined by the tuple of \a this
2024 * array containing 6 components.
2025 * \return DataArrayDouble * - the new instance of DataArrayDouble, whose each tuple
2026 * is calculated from the tuple <em>(t)</em> of \a this array as follows:
2027 * \f$ t[0]^2+t[1]^2+t[2]^2+2*t[3]^2+2*t[4]^2+2*t[5]^2\f$.
2028 * The caller is to delete this result array using decrRef() as it is no more needed.
2029 * \throw If \a this->getNumberOfComponents() != 6.
2031 DataArrayDouble *DataArrayDouble::doublyContractedProduct() const
2034 std::size_t nbOfComp(getNumberOfComponents());
2036 throw INTERP_KERNEL::Exception("DataArrayDouble::doublyContractedProduct : must be an array with exactly 6 components !");
2037 DataArrayDouble *ret=DataArrayDouble::New();
2038 mcIdType nbOfTuple=getNumberOfTuples();
2039 ret->alloc(nbOfTuple,1);
2040 const double *src=getConstPointer();
2041 double *dest=ret->getPointer();
2042 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=6)
2043 *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];
2048 * Computes the determinant of every square matrix defined by the tuple of \a this
2049 * array, which contains either 4, 6 or 9 components. The case of 6 components
2050 * corresponds to that of the upper triangular matrix.
2051 * \return DataArrayDouble * - the new instance of DataArrayDouble, whose each tuple
2052 * is the determinant of matrix of the corresponding tuple of \a this array.
2053 * The caller is to delete this result array using decrRef() as it is no more
2055 * \throw If \a this->getNumberOfComponents() is not in [4,6,9].
2057 DataArrayDouble *DataArrayDouble::determinant() const
2060 DataArrayDouble *ret=DataArrayDouble::New();
2061 mcIdType nbOfTuple=getNumberOfTuples();
2062 ret->alloc(nbOfTuple,1);
2063 const double *src=getConstPointer();
2064 double *dest=ret->getPointer();
2065 switch(getNumberOfComponents())
2068 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=6)
2069 *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];
2072 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=4)
2073 *dest=src[0]*src[3]-src[1]*src[2];
2076 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=9)
2077 *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];
2081 throw INTERP_KERNEL::Exception("DataArrayDouble::determinant : Invalid number of components ! must be in 4,6,9 !");
2086 * Computes 3 eigenvalues of every upper triangular matrix defined by the tuple of
2087 * \a this array, which contains 6 components.
2088 * \return DataArrayDouble * - the new instance of DataArrayDouble containing 3
2089 * components, whose each tuple contains the eigenvalues of the matrix of
2090 * corresponding tuple of \a this array.
2091 * The caller is to delete this result array using decrRef() as it is no more
2093 * \throw If \a this->getNumberOfComponents() != 6.
2095 DataArrayDouble *DataArrayDouble::eigenValues() const
2098 std::size_t nbOfComp=getNumberOfComponents();
2100 throw INTERP_KERNEL::Exception("DataArrayDouble::eigenValues : must be an array with exactly 6 components !");
2101 DataArrayDouble *ret=DataArrayDouble::New();
2102 mcIdType nbOfTuple=getNumberOfTuples();
2103 ret->alloc(nbOfTuple,3);
2104 const double *src=getConstPointer();
2105 double *dest=ret->getPointer();
2106 for(mcIdType i=0;i<nbOfTuple;i++,dest+=3,src+=6)
2107 INTERP_KERNEL::computeEigenValues6(src,dest);
2112 * Computes 3 eigenvectors of every upper triangular matrix defined by the tuple of
2113 * \a this array, which contains 6 components.
2114 * \return DataArrayDouble * - the new instance of DataArrayDouble containing 9
2115 * components, whose each tuple contains 3 eigenvectors of the matrix of
2116 * corresponding tuple of \a this array.
2117 * The caller is to delete this result array using decrRef() as it is no more
2119 * \throw If \a this->getNumberOfComponents() != 6.
2121 DataArrayDouble *DataArrayDouble::eigenVectors() const
2124 std::size_t nbOfComp=getNumberOfComponents();
2126 throw INTERP_KERNEL::Exception("DataArrayDouble::eigenVectors : must be an array with exactly 6 components !");
2127 DataArrayDouble *ret=DataArrayDouble::New();
2128 mcIdType nbOfTuple=getNumberOfTuples();
2129 ret->alloc(nbOfTuple,9);
2130 const double *src=getConstPointer();
2131 double *dest=ret->getPointer();
2132 for(mcIdType i=0;i<nbOfTuple;i++,src+=6)
2135 INTERP_KERNEL::computeEigenValues6(src,tmp);
2136 for(mcIdType j=0;j<3;j++,dest+=3)
2137 INTERP_KERNEL::computeEigenVectorForEigenValue6(src,tmp[j],1e-12,dest);
2143 * Computes the inverse matrix of every matrix defined by the tuple of \a this
2144 * array, which contains either 4, 6 or 9 components. The case of 6 components
2145 * corresponds to that of the upper triangular matrix.
2146 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2147 * same number of components as \a this one, whose each tuple is the inverse
2148 * matrix of the matrix of corresponding tuple of \a this array.
2149 * The caller is to delete this result array using decrRef() as it is no more
2151 * \throw If \a this->getNumberOfComponents() is not in [4,6,9].
2153 DataArrayDouble *DataArrayDouble::inverse() const
2156 std::size_t nbOfComp=getNumberOfComponents();
2157 if(nbOfComp!=6 && nbOfComp!=9 && nbOfComp!=4)
2158 throw INTERP_KERNEL::Exception("DataArrayDouble::inversion : must be an array with 4,6 or 9 components !");
2159 DataArrayDouble *ret=DataArrayDouble::New();
2160 mcIdType nbOfTuple=getNumberOfTuples();
2161 ret->alloc(nbOfTuple,nbOfComp);
2162 const double *src=getConstPointer();
2163 double *dest=ret->getPointer();
2165 for(mcIdType i=0;i<nbOfTuple;i++,dest+=6,src+=6)
2167 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];
2168 dest[0]=(src[1]*src[2]-src[4]*src[4])/det;
2169 dest[1]=(src[0]*src[2]-src[5]*src[5])/det;
2170 dest[2]=(src[0]*src[1]-src[3]*src[3])/det;
2171 dest[3]=(src[5]*src[4]-src[3]*src[2])/det;
2172 dest[4]=(src[5]*src[3]-src[0]*src[4])/det;
2173 dest[5]=(src[3]*src[4]-src[1]*src[5])/det;
2175 else if(nbOfComp==4)
2176 for(mcIdType i=0;i<nbOfTuple;i++,dest+=4,src+=4)
2178 double det=src[0]*src[3]-src[1]*src[2];
2180 dest[1]=-src[1]/det;
2181 dest[2]=-src[2]/det;
2185 for(mcIdType i=0;i<nbOfTuple;i++,dest+=9,src+=9)
2187 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];
2188 dest[0]=(src[4]*src[8]-src[7]*src[5])/det;
2189 dest[1]=(src[7]*src[2]-src[1]*src[8])/det;
2190 dest[2]=(src[1]*src[5]-src[4]*src[2])/det;
2191 dest[3]=(src[6]*src[5]-src[3]*src[8])/det;
2192 dest[4]=(src[0]*src[8]-src[6]*src[2])/det;
2193 dest[5]=(src[2]*src[3]-src[0]*src[5])/det;
2194 dest[6]=(src[3]*src[7]-src[6]*src[4])/det;
2195 dest[7]=(src[6]*src[1]-src[0]*src[7])/det;
2196 dest[8]=(src[0]*src[4]-src[1]*src[3])/det;
2202 * Computes the trace of every matrix defined by the tuple of \a this
2203 * array, which contains either 4, 6 or 9 components. The case of 6 components
2204 * corresponds to that of the upper triangular matrix.
2205 * \return DataArrayDouble * - the new instance of DataArrayDouble containing
2206 * 1 component, whose each tuple is the trace of
2207 * the matrix of corresponding tuple of \a this array.
2208 * The caller is to delete this result array using decrRef() as it is no more
2210 * \throw If \a this->getNumberOfComponents() is not in [4,6,9].
2212 DataArrayDouble *DataArrayDouble::trace() const
2215 std::size_t nbOfComp=getNumberOfComponents();
2216 if(nbOfComp!=6 && nbOfComp!=9 && nbOfComp!=4)
2217 throw INTERP_KERNEL::Exception("DataArrayDouble::trace : must be an array with 4,6 or 9 components !");
2218 DataArrayDouble *ret=DataArrayDouble::New();
2219 mcIdType nbOfTuple=getNumberOfTuples();
2220 ret->alloc(nbOfTuple,1);
2221 const double *src=getConstPointer();
2222 double *dest=ret->getPointer();
2224 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=6)
2225 *dest=src[0]+src[1]+src[2];
2226 else if(nbOfComp==4)
2227 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=4)
2228 *dest=src[0]+src[3];
2230 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=9)
2231 *dest=src[0]+src[4]+src[8];
2236 * Computes the stress deviator tensor of every stress tensor defined by the tuple of
2237 * \a this array, which contains 6 components.
2238 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2239 * same number of components and tuples as \a this array.
2240 * The caller is to delete this result array using decrRef() as it is no more
2242 * \throw If \a this->getNumberOfComponents() != 6.
2244 DataArrayDouble *DataArrayDouble::deviator() const
2247 std::size_t nbOfComp=getNumberOfComponents();
2249 throw INTERP_KERNEL::Exception("DataArrayDouble::deviator : must be an array with exactly 6 components !");
2250 DataArrayDouble *ret=DataArrayDouble::New();
2251 mcIdType nbOfTuple=getNumberOfTuples();
2252 ret->alloc(nbOfTuple,6);
2253 const double *src=getConstPointer();
2254 double *dest=ret->getPointer();
2255 for(mcIdType i=0;i<nbOfTuple;i++,dest+=6,src+=6)
2257 double tr=(src[0]+src[1]+src[2])/3.;
2269 * Computes the magnitude of every vector defined by the tuple of
2271 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2272 * same number of tuples as \a this array and one component.
2273 * The caller is to delete this result array using decrRef() as it is no more
2275 * \throw If \a this is not allocated.
2277 DataArrayDouble *DataArrayDouble::magnitude() const
2280 std::size_t nbOfComp=getNumberOfComponents();
2281 DataArrayDouble *ret=DataArrayDouble::New();
2282 mcIdType nbOfTuple=getNumberOfTuples();
2283 ret->alloc(nbOfTuple,1);
2284 const double *src=getConstPointer();
2285 double *dest=ret->getPointer();
2286 for(mcIdType i=0;i<nbOfTuple;i++,dest++)
2289 for(std::size_t j=0;j<nbOfComp;j++,src++)
2297 * Computes the maximal value within every tuple of \a this array.
2298 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2299 * same number of tuples as \a this array and one component.
2300 * The caller is to delete this result array using decrRef() as it is no more
2302 * \throw If \a this is not allocated.
2303 * \sa DataArrayDouble::maxPerTupleWithCompoId
2305 DataArrayDouble *DataArrayDouble::maxPerTuple() const
2308 std::size_t nbOfComp(getNumberOfComponents());
2309 MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
2310 mcIdType nbOfTuple(getNumberOfTuples());
2311 ret->alloc(nbOfTuple,1);
2312 const double *src=getConstPointer();
2313 double *dest=ret->getPointer();
2314 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=nbOfComp)
2315 *dest=*std::max_element(src,src+nbOfComp);
2320 * Computes the maximal value within every tuple of \a this array and it returns the first component
2321 * id for each tuple that corresponds to the maximal value within the tuple.
2323 * \param [out] compoIdOfMaxPerTuple - the new new instance of DataArrayInt containing the
2324 * same number of tuples and only one component.
2325 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2326 * same number of tuples as \a this array and one component.
2327 * The caller is to delete this result array using decrRef() as it is no more
2329 * \throw If \a this is not allocated.
2330 * \sa DataArrayDouble::maxPerTuple
2332 DataArrayDouble *DataArrayDouble::maxPerTupleWithCompoId(DataArrayIdType* &compoIdOfMaxPerTuple) const
2335 std::size_t nbOfComp(getNumberOfComponents());
2336 MCAuto<DataArrayDouble> ret0=DataArrayDouble::New();
2337 MCAuto<DataArrayIdType> ret1=DataArrayIdType::New();
2338 mcIdType nbOfTuple=getNumberOfTuples();
2339 ret0->alloc(nbOfTuple,1); ret1->alloc(nbOfTuple,1);
2340 const double *src=getConstPointer();
2341 double *dest=ret0->getPointer(); mcIdType *dest1=ret1->getPointer();
2342 for(mcIdType i=0;i<nbOfTuple;i++,dest++,dest1++,src+=nbOfComp)
2344 const double *loc=std::max_element(src,src+nbOfComp);
2346 *dest1=ToIdType(std::distance(src,loc));
2348 compoIdOfMaxPerTuple=ret1.retn();
2353 * This method returns a newly allocated DataArrayDouble instance having one component and \c this->getNumberOfTuples() * \c this->getNumberOfTuples() tuples.
2354 * \n This returned array contains the euclidian distance for each tuple in \a this.
2355 * \n So the returned array can be seen as a dense symmetrical matrix whose diagonal elements are equal to 0.
2356 * \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)
2358 * \warning use this method with care because it can leads to big amount of consumed memory !
2360 * \return A newly allocated (huge) MEDCoupling::DataArrayDouble instance that the caller should deal with.
2362 * \throw If \a this is not allocated.
2364 * \sa DataArrayDouble::buildEuclidianDistanceDenseMatrixWith
2366 DataArrayDouble *DataArrayDouble::buildEuclidianDistanceDenseMatrix() const
2369 std::size_t nbOfComp(getNumberOfComponents());
2370 mcIdType nbOfTuples(getNumberOfTuples());
2371 const double *inData=getConstPointer();
2372 MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
2373 ret->alloc(nbOfTuples*nbOfTuples,1);
2374 double *outData=ret->getPointer();
2375 for(mcIdType i=0;i<nbOfTuples;i++)
2377 outData[i*nbOfTuples+i]=0.;
2378 for(mcIdType j=i+1;j<nbOfTuples;j++)
2381 for(std::size_t k=0;k<nbOfComp;k++)
2382 { double delta=inData[i*nbOfComp+k]-inData[j*nbOfComp+k]; dist+=delta*delta; }
2384 outData[i*nbOfTuples+j]=dist;
2385 outData[j*nbOfTuples+i]=dist;
2392 * This method returns a newly allocated DataArrayDouble instance having one component and \c this->getNumberOfTuples() * \c other->getNumberOfTuples() tuples.
2393 * \n This returned array contains the euclidian distance for each tuple in \a other with each tuple in \a this.
2394 * \n So the returned array can be seen as a dense rectangular matrix with \c other->getNumberOfTuples() rows and \c this->getNumberOfTuples() columns.
2395 * \n Output rectangular matrix is sorted along rows.
2396 * \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)
2398 * \warning use this method with care because it can leads to big amount of consumed memory !
2400 * \param [in] other DataArrayDouble instance having same number of components than \a this.
2401 * \return A newly allocated (huge) MEDCoupling::DataArrayDouble instance that the caller should deal with.
2403 * \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.
2405 * \sa DataArrayDouble::buildEuclidianDistanceDenseMatrix
2407 DataArrayDouble *DataArrayDouble::buildEuclidianDistanceDenseMatrixWith(const DataArrayDouble *other) const
2410 throw INTERP_KERNEL::Exception("DataArrayDouble::buildEuclidianDistanceDenseMatrixWith : input parameter is null !");
2412 other->checkAllocated();
2413 std::size_t nbOfComp(getNumberOfComponents());
2414 std::size_t otherNbOfComp(other->getNumberOfComponents());
2415 if(nbOfComp!=otherNbOfComp)
2417 std::ostringstream oss; oss << "DataArrayDouble::buildEuclidianDistanceDenseMatrixWith : this nb of compo=" << nbOfComp << " and other nb of compo=" << otherNbOfComp << ". It should match !";
2418 throw INTERP_KERNEL::Exception(oss.str().c_str());
2420 mcIdType nbOfTuples(getNumberOfTuples());
2421 mcIdType otherNbOfTuples(other->getNumberOfTuples());
2422 const double *inData=getConstPointer();
2423 const double *inDataOther=other->getConstPointer();
2424 MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
2425 ret->alloc(otherNbOfTuples*nbOfTuples,1);
2426 double *outData=ret->getPointer();
2427 for(mcIdType i=0;i<otherNbOfTuples;i++,inDataOther+=nbOfComp)
2429 for(mcIdType j=0;j<nbOfTuples;j++)
2432 for(std::size_t k=0;k<nbOfComp;k++)
2433 { double delta=inDataOther[k]-inData[j*nbOfComp+k]; dist+=delta*delta; }
2435 outData[i*nbOfTuples+j]=dist;
2442 * This method expects that \a this stores 3 tuples containing 2 components each.
2443 * Each of this tuples represent a point into 2D space.
2444 * 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).
2445 * 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[.
2447 * \throw If \a this is not allocated.
2448 * \throw If \a this has not 3 tuples of 2 components
2449 * \throw If tuples/points in \a this are aligned
2451 void DataArrayDouble::asArcOfCircle(double center[2], double& radius, double& ang) const
2454 INTERP_KERNEL::QuadraticPlanarPrecision arcPrec(1e-14);
2455 if(getNumberOfTuples()!=3 && getNumberOfComponents()!=2)
2456 throw INTERP_KERNEL::Exception("DataArrayDouble::asArcCircle : this method expects");
2457 const double *pt(begin());
2458 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]));
2460 INTERP_KERNEL::AutoCppPtr<INTERP_KERNEL::EdgeLin> e1(new INTERP_KERNEL::EdgeLin(n0,n2)),e2(new INTERP_KERNEL::EdgeLin(n2,n1));
2461 INTERP_KERNEL::SegSegIntersector inters(*e1,*e2);
2462 bool colinearity(inters.areColinears());
2464 throw INTERP_KERNEL::Exception("DataArrayDouble::asArcOfCircle : 3 points in this have been detected as colinear !");
2466 INTERP_KERNEL::AutoCppPtr<INTERP_KERNEL::EdgeArcCircle> ret(new INTERP_KERNEL::EdgeArcCircle(n0,n2,n1));
2467 const double *c(ret->getCenter());
2468 center[0]=c[0]; center[1]=c[1];
2469 radius=ret->getRadius();
2470 ang=ret->getAngle();
2474 * Sorts value within every tuple of \a this array.
2475 * \param [in] asc - if \a true, the values are sorted in ascending order, else,
2476 * in descending order.
2477 * \throw If \a this is not allocated.
2479 void DataArrayDouble::sortPerTuple(bool asc)
2482 double *pt=getPointer();
2483 mcIdType nbOfTuple(getNumberOfTuples());
2484 std::size_t nbOfComp(getNumberOfComponents());
2486 for(mcIdType i=0;i<nbOfTuple;i++,pt+=nbOfComp)
2487 std::sort(pt,pt+nbOfComp);
2489 for(mcIdType i=0;i<nbOfTuple;i++,pt+=nbOfComp)
2490 std::sort(pt,pt+nbOfComp,std::greater<double>());
2495 * Modify all elements of \a this array, so that
2496 * an element _x_ becomes \f$ numerator / x \f$.
2497 * \warning If an exception is thrown because of presence of 0.0 element in \a this
2498 * array, all elements processed before detection of the zero element remain
2500 * \param [in] numerator - the numerator used to modify array elements.
2501 * \throw If \a this is not allocated.
2502 * \throw If there is an element equal to 0.0 in \a this array.
2504 void DataArrayDouble::applyInv(double numerator)
2507 double *ptr=getPointer();
2508 std::size_t nbOfElems=getNbOfElems();
2509 for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2511 if(std::abs(*ptr)>std::numeric_limits<double>::min())
2513 *ptr=numerator/(*ptr);
2517 std::ostringstream oss; oss << "DataArrayDouble::applyInv : presence of null value in tuple #" << i/getNumberOfComponents() << " component #" << i%getNumberOfComponents();
2519 throw INTERP_KERNEL::Exception(oss.str().c_str());
2526 * Modify all elements of \a this array, so that
2527 * an element _x_ becomes <em> val ^ x </em>. Contrary to DataArrayInt::applyPow
2528 * all values in \a this have to be >= 0 if val is \b not integer.
2529 * \param [in] val - the value used to apply pow on all array elements.
2530 * \throw If \a this is not allocated.
2531 * \warning If an exception is thrown because of presence of 0 element in \a this
2532 * array and \a val is \b not integer, all elements processed before detection of the zero element remain
2535 void DataArrayDouble::applyPow(double val)
2538 double *ptr=getPointer();
2539 std::size_t nbOfElems=getNbOfElems();
2541 bool isInt=((double)val2)==val;
2544 for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2550 std::ostringstream oss; oss << "DataArrayDouble::applyPow (double) : At elem # " << i << " value is " << *ptr << " ! must be >=0. !";
2551 throw INTERP_KERNEL::Exception(oss.str().c_str());
2557 for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2558 *ptr=pow(*ptr,val2);
2564 * Modify all elements of \a this array, so that
2565 * an element _x_ becomes \f$ val ^ x \f$.
2566 * \param [in] val - the value used to apply pow on all array elements.
2567 * \throw If \a this is not allocated.
2568 * \throw If \a val < 0.
2569 * \warning If an exception is thrown because of presence of 0 element in \a this
2570 * array, all elements processed before detection of the zero element remain
2573 void DataArrayDouble::applyRPow(double val)
2577 throw INTERP_KERNEL::Exception("DataArrayDouble::applyRPow : the input value has to be >= 0 !");
2578 double *ptr=getPointer();
2579 std::size_t nbOfElems=getNbOfElems();
2580 for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2586 * Returns a new DataArrayDouble created from \a this one by applying \a
2587 * FunctionToEvaluate to every tuple of \a this array. Textual data is not copied.
2588 * For more info see \ref MEDCouplingArrayApplyFunc
2589 * \param [in] nbOfComp - number of components in the result array.
2590 * \param [in] func - the \a FunctionToEvaluate declared as
2591 * \c bool (*\a func)(\c const \c double *\a pos, \c double *\a res),
2592 * where \a pos points to the first component of a tuple of \a this array
2593 * and \a res points to the first component of a tuple of the result array.
2594 * Note that length (number of components) of \a pos can differ from
2596 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2597 * same number of tuples as \a this array.
2598 * The caller is to delete this result array using decrRef() as it is no more
2600 * \throw If \a this is not allocated.
2601 * \throw If \a func returns \a false.
2603 DataArrayDouble *DataArrayDouble::applyFunc(std::size_t nbOfComp, FunctionToEvaluate func) const
2606 DataArrayDouble *newArr=DataArrayDouble::New();
2607 mcIdType nbOfTuples(getNumberOfTuples());
2608 std::size_t oldNbOfComp(getNumberOfComponents());
2609 newArr->alloc(nbOfTuples,nbOfComp);
2610 const double *ptr=getConstPointer();
2611 double *ptrToFill=newArr->getPointer();
2612 for(mcIdType i=0;i<nbOfTuples;i++)
2614 if(!func(ptr+i*oldNbOfComp,ptrToFill+i*nbOfComp))
2616 std::ostringstream oss; oss << "For tuple # " << i << " with value (";
2617 std::copy(ptr+oldNbOfComp*i,ptr+oldNbOfComp*(i+1),std::ostream_iterator<double>(oss,", "));
2618 oss << ") : Evaluation of function failed !";
2620 throw INTERP_KERNEL::Exception(oss.str().c_str());
2627 * Returns a new DataArrayDouble created from \a this one by applying a function to every
2628 * tuple of \a this array. Textual data is not copied.
2629 * For more info see \ref MEDCouplingArrayApplyFunc1.
2630 * \param [in] nbOfComp - number of components in the result array.
2631 * \param [in] func - the expression defining how to transform a tuple of \a this array.
2632 * Supported expressions are described \ref MEDCouplingArrayApplyFuncExpr "here".
2633 * \param [in] isSafe - By default true. If true invalid operation (division by 0. acos of value > 1. ...) leads to a throw of an exception.
2634 * If false the computation is carried on without any notification. When false the evaluation is a little faster.
2635 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2636 * same number of tuples as \a this array and \a nbOfComp components.
2637 * The caller is to delete this result array using decrRef() as it is no more
2639 * \throw If \a this is not allocated.
2640 * \throw If computing \a func fails.
2642 DataArrayDouble *DataArrayDouble::applyFunc(std::size_t nbOfComp, const std::string& func, bool isSafe) const
2644 INTERP_KERNEL::ExprParser expr(func);
2646 std::set<std::string> vars;
2647 expr.getTrueSetOfVars(vars);
2648 std::vector<std::string> varsV(vars.begin(),vars.end());
2649 return applyFuncNamedCompo(nbOfComp,varsV,func,isSafe);
2653 * Returns a new DataArrayDouble created from \a this one by applying a function to every
2654 * tuple of \a this array. Textual data is not copied. This method works by tuples (whatever its size).
2655 * If \a this is a one component array, call applyFuncOnThis instead that performs the same work faster.
2657 * For more info see \ref MEDCouplingArrayApplyFunc0.
2658 * \param [in] func - the expression defining how to transform a tuple of \a this array.
2659 * Supported expressions are described \ref MEDCouplingArrayApplyFuncExpr "here".
2660 * \param [in] isSafe - By default true. If true invalid operation (division by 0. acos of value > 1. ...) leads to a throw of an exception.
2661 * If false the computation is carried on without any notification. When false the evaluation is a little faster.
2662 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2663 * same number of tuples and components as \a this array.
2664 * The caller is to delete this result array using decrRef() as it is no more
2666 * \sa applyFuncOnThis
2667 * \throw If \a this is not allocated.
2668 * \throw If computing \a func fails.
2670 DataArrayDouble *DataArrayDouble::applyFunc(const std::string& func, bool isSafe) const
2672 std::size_t nbOfComp(getNumberOfComponents());
2674 throw INTERP_KERNEL::Exception("DataArrayDouble::applyFunc : output number of component must be > 0 !");
2676 mcIdType nbOfTuples(getNumberOfTuples());
2677 MCAuto<DataArrayDouble> newArr(DataArrayDouble::New());
2678 newArr->alloc(nbOfTuples,nbOfComp);
2679 INTERP_KERNEL::ExprParser expr(func);
2681 std::set<std::string> vars;
2682 expr.getTrueSetOfVars(vars);
2685 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 : ";
2686 std::copy(vars.begin(),vars.end(),std::ostream_iterator<std::string>(oss," "));
2687 throw INTERP_KERNEL::Exception(oss.str().c_str());
2691 expr.prepareFastEvaluator();
2692 newArr->rearrange(1);
2693 newArr->fillWithValue(expr.evaluateDouble());
2694 newArr->rearrange(nbOfComp);
2695 return newArr.retn();
2697 std::vector<std::string> vars2(vars.begin(),vars.end());
2698 double buff,*ptrToFill(newArr->getPointer());
2699 const double *ptr(begin());
2700 std::vector<double> stck;
2701 expr.prepareExprEvaluationDouble(vars2,1,1,0,&buff,&buff+1);
2702 expr.prepareFastEvaluator();
2705 for(mcIdType i=0;i<nbOfTuples;i++)
2707 for(std::size_t iComp=0;iComp<nbOfComp;iComp++,ptr++,ptrToFill++)
2710 expr.evaluateDoubleInternal(stck);
2711 *ptrToFill=stck.back();
2718 for(mcIdType i=0;i<nbOfTuples;i++)
2720 for(std::size_t iComp=0;iComp<nbOfComp;iComp++,ptr++,ptrToFill++)
2725 expr.evaluateDoubleInternalSafe(stck);
2727 catch(INTERP_KERNEL::Exception& e)
2729 std::ostringstream oss; oss << "For tuple # " << i << " component # " << iComp << " with value (";
2731 oss << ") : Evaluation of function failed !" << e.what();
2732 throw INTERP_KERNEL::Exception(oss.str().c_str());
2734 *ptrToFill=stck.back();
2739 return newArr.retn();
2743 * This method is a non const method that modify the array in \a this.
2744 * This method only works on one component array. It means that function \a func must
2745 * contain at most one variable.
2746 * This method is a specialization of applyFunc method with one parameter on one component array.
2748 * \param [in] func - the expression defining how to transform a tuple of \a this array.
2749 * Supported expressions are described \ref MEDCouplingArrayApplyFuncExpr "here".
2750 * \param [in] isSafe - By default true. If true invalid operation (division by 0. acos of value > 1. ...) leads to a throw of an exception.
2751 * If false the computation is carried on without any notification. When false the evaluation is a little faster.
2755 void DataArrayDouble::applyFuncOnThis(const std::string& func, bool isSafe)
2757 std::size_t nbOfComp(getNumberOfComponents());
2759 throw INTERP_KERNEL::Exception("DataArrayDouble::applyFuncOnThis : output number of component must be > 0 !");
2761 mcIdType nbOfTuples(getNumberOfTuples());
2762 INTERP_KERNEL::ExprParser expr(func);
2764 std::set<std::string> vars;
2765 expr.getTrueSetOfVars(vars);
2768 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 : ";
2769 std::copy(vars.begin(),vars.end(),std::ostream_iterator<std::string>(oss," "));
2770 throw INTERP_KERNEL::Exception(oss.str().c_str());
2774 expr.prepareFastEvaluator();
2775 std::vector<std::string> compInfo(getInfoOnComponents());
2777 fillWithValue(expr.evaluateDouble());
2778 rearrange(nbOfComp);
2779 setInfoOnComponents(compInfo);
2782 std::vector<std::string> vars2(vars.begin(),vars.end());
2783 double buff,*ptrToFill(getPointer());
2784 const double *ptr(begin());
2785 std::vector<double> stck;
2786 expr.prepareExprEvaluationDouble(vars2,1,1,0,&buff,&buff+1);
2787 expr.prepareFastEvaluator();
2790 for(mcIdType i=0;i<nbOfTuples;i++)
2792 for(std::size_t iComp=0;iComp<nbOfComp;iComp++,ptr++,ptrToFill++)
2795 expr.evaluateDoubleInternal(stck);
2796 *ptrToFill=stck.back();
2803 for(mcIdType i=0;i<nbOfTuples;i++)
2805 for(std::size_t iComp=0;iComp<nbOfComp;iComp++,ptr++,ptrToFill++)
2810 expr.evaluateDoubleInternalSafe(stck);
2812 catch(INTERP_KERNEL::Exception& e)
2814 std::ostringstream oss; oss << "For tuple # " << i << " component # " << iComp << " with value (";
2816 oss << ") : Evaluation of function failed !" << e.what();
2817 throw INTERP_KERNEL::Exception(oss.str().c_str());
2819 *ptrToFill=stck.back();
2827 * Returns a new DataArrayDouble created from \a this one by applying a function to every
2828 * tuple of \a this array. Textual data is not copied.
2829 * For more info see \ref MEDCouplingArrayApplyFunc2.
2830 * \param [in] nbOfComp - number of components in the result array.
2831 * \param [in] func - the expression defining how to transform a tuple of \a this array.
2832 * Supported expressions are described \ref MEDCouplingArrayApplyFuncExpr "here".
2833 * \param [in] isSafe - By default true. If true invalid operation (division by 0. acos of value > 1. ...) leads to a throw of an exception.
2834 * If false the computation is carried on without any notification. When false the evaluation is a little faster.
2835 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2836 * same number of tuples as \a this array.
2837 * The caller is to delete this result array using decrRef() as it is no more
2839 * \throw If \a this is not allocated.
2840 * \throw If \a func contains vars that are not in \a this->getInfoOnComponent().
2841 * \throw If computing \a func fails.
2843 DataArrayDouble *DataArrayDouble::applyFuncCompo(std::size_t nbOfComp, const std::string& func, bool isSafe) const
2845 return applyFuncNamedCompo(nbOfComp,getVarsOnComponent(),func,isSafe);
2849 * Returns a new DataArrayDouble created from \a this one by applying a function to every
2850 * tuple of \a this array. Textual data is not copied.
2851 * For more info see \ref MEDCouplingArrayApplyFunc3.
2852 * \param [in] nbOfComp - number of components in the result array.
2853 * \param [in] varsOrder - sequence of vars defining their order.
2854 * \param [in] func - the expression defining how to transform a tuple of \a this array.
2855 * Supported expressions are described \ref MEDCouplingArrayApplyFuncExpr "here".
2856 * \param [in] isSafe - By default true. If true invalid operation (division by 0. acos of value > 1. ...) leads to a throw of an exception.
2857 * If false the computation is carried on without any notification. When false the evaluation is a little faster.
2858 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2859 * same number of tuples as \a this array.
2860 * The caller is to delete this result array using decrRef() as it is no more
2862 * \throw If \a this is not allocated.
2863 * \throw If \a func contains vars not in \a varsOrder.
2864 * \throw If computing \a func fails.
2866 DataArrayDouble *DataArrayDouble::applyFuncNamedCompo(std::size_t nbOfComp, const std::vector<std::string>& varsOrder, const std::string& func, bool isSafe) const
2869 throw INTERP_KERNEL::Exception("DataArrayDouble::applyFuncNamedCompo : output number of component must be > 0 !");
2870 std::vector<std::string> varsOrder2(varsOrder);
2871 std::size_t oldNbOfComp(getNumberOfComponents());
2872 for(std::size_t i=varsOrder.size();i<oldNbOfComp;i++)
2873 varsOrder2.push_back(std::string());
2875 mcIdType nbOfTuples(getNumberOfTuples());
2876 INTERP_KERNEL::ExprParser expr(func);
2878 std::set<std::string> vars;
2879 expr.getTrueSetOfVars(vars);
2880 if(vars.size()>oldNbOfComp)
2882 std::ostringstream oss; oss << "The field has " << oldNbOfComp << " components and there are ";
2883 oss << vars.size() << " variables : ";
2884 std::copy(vars.begin(),vars.end(),std::ostream_iterator<std::string>(oss," "));
2885 throw INTERP_KERNEL::Exception(oss.str().c_str());
2887 MCAuto<DataArrayDouble> newArr(DataArrayDouble::New());
2888 newArr->alloc(nbOfTuples,nbOfComp);
2889 INTERP_KERNEL::AutoPtr<double> buff(new double[oldNbOfComp]);
2890 double *buffPtr(buff),*ptrToFill;
2891 std::vector<double> stck;
2892 for(std::size_t iComp=0;iComp<nbOfComp;iComp++)
2894 expr.prepareExprEvaluationDouble(varsOrder2,(int)oldNbOfComp,(int)nbOfComp,(int)iComp,buffPtr,buffPtr+oldNbOfComp);
2895 expr.prepareFastEvaluator();
2896 const double *ptr(getConstPointer());
2897 ptrToFill=newArr->getPointer()+iComp;
2900 for(mcIdType i=0;i<nbOfTuples;i++,ptrToFill+=nbOfComp,ptr+=oldNbOfComp)
2902 std::copy(ptr,ptr+oldNbOfComp,buffPtr);
2903 expr.evaluateDoubleInternal(stck);
2904 *ptrToFill=stck.back();
2910 for(mcIdType i=0;i<nbOfTuples;i++,ptrToFill+=nbOfComp,ptr+=oldNbOfComp)
2912 std::copy(ptr,ptr+oldNbOfComp,buffPtr);
2915 expr.evaluateDoubleInternalSafe(stck);
2916 *ptrToFill=stck.back();
2919 catch(INTERP_KERNEL::Exception& e)
2921 std::ostringstream oss; oss << "For tuple # " << i << " with value (";
2922 std::copy(ptr+oldNbOfComp*i,ptr+oldNbOfComp*(i+1),std::ostream_iterator<double>(oss,", "));
2923 oss << ") : Evaluation of function failed !" << e.what();
2924 throw INTERP_KERNEL::Exception(oss.str().c_str());
2929 return newArr.retn();
2932 void DataArrayDouble::applyFuncFast32(const std::string& func)
2935 INTERP_KERNEL::ExprParser expr(func);
2937 char *funcStr=expr.compileX86();
2939 *((void **)&funcPtr)=funcStr;//he he...
2941 double *ptr=getPointer();
2942 std::size_t nbOfComp=getNumberOfComponents();
2943 mcIdType nbOfTuples=getNumberOfTuples();
2944 std::size_t nbOfElems=nbOfTuples*nbOfComp;
2945 for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2950 void DataArrayDouble::applyFuncFast64(const std::string& func)
2953 INTERP_KERNEL::ExprParser expr(func);
2955 char *funcStr=expr.compileX86_64();
2957 *((void **)&funcPtr)=funcStr;//he he...
2959 double *ptr=getPointer();
2960 std::size_t nbOfComp=getNumberOfComponents();
2961 mcIdType nbOfTuples=getNumberOfTuples();
2962 std::size_t nbOfElems=nbOfTuples*nbOfComp;
2963 for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2969 * \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.
2971 MCAuto<DataArrayDouble> DataArrayDouble::symmetry3DPlane(const double point[3], const double normalVector[3]) const
2974 if(getNumberOfComponents()!=3)
2975 throw INTERP_KERNEL::Exception("DataArrayDouble::symmetry3DPlane : this is excepted to have 3 components !");
2976 mcIdType nbTuples(getNumberOfTuples());
2977 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
2978 ret->alloc(nbTuples,3);
2979 Symmetry3DPlane(point,normalVector,nbTuples,begin(),ret->getPointer());
2983 DataArrayDoubleIterator *DataArrayDouble::iterator()
2985 return new DataArrayDoubleIterator(this);
2989 * Returns a new DataArrayInt containing indices of tuples of \a this one-dimensional
2990 * array whose values are within a given range. Textual data is not copied.
2991 * \param [in] vmin - a lowest acceptable value (included).
2992 * \param [in] vmax - a greatest acceptable value (included).
2993 * \return DataArrayInt * - the new instance of DataArrayInt.
2994 * The caller is to delete this result array using decrRef() as it is no more
2996 * \throw If \a this->getNumberOfComponents() != 1.
2998 * \sa DataArrayDouble::findIdsNotInRange
3000 * \if ENABLE_EXAMPLES
3001 * \ref cpp_mcdataarraydouble_getidsinrange "Here is a C++ example".<br>
3002 * \ref py_mcdataarraydouble_getidsinrange "Here is a Python example".
3005 DataArrayIdType *DataArrayDouble::findIdsInRange(double vmin, double vmax) const
3008 if(getNumberOfComponents()!=1)
3009 throw INTERP_KERNEL::Exception("DataArrayDouble::findIdsInRange : this must have exactly one component !");
3010 const double *cptr(begin());
3011 MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(0,1);
3012 mcIdType nbOfTuples(getNumberOfTuples());
3013 for(mcIdType i=0;i<nbOfTuples;i++,cptr++)
3014 if(*cptr>=vmin && *cptr<=vmax)
3015 ret->pushBackSilent(i);
3020 * Returns a new DataArrayInt containing indices of tuples of \a this one-dimensional
3021 * array whose values are not within a given range. Textual data is not copied.
3022 * \param [in] vmin - a lowest not acceptable value (excluded).
3023 * \param [in] vmax - a greatest not acceptable value (excluded).
3024 * \return DataArrayInt * - the new instance of DataArrayInt.
3025 * The caller is to delete this result array using decrRef() as it is no more
3027 * \throw If \a this->getNumberOfComponents() != 1.
3029 * \sa DataArrayDouble::findIdsInRange
3031 DataArrayIdType *DataArrayDouble::findIdsNotInRange(double vmin, double vmax) const
3034 if(getNumberOfComponents()!=1)
3035 throw INTERP_KERNEL::Exception("DataArrayDouble::findIdsNotInRange : this must have exactly one component !");
3036 const double *cptr(begin());
3037 MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(0,1);
3038 mcIdType nbOfTuples(getNumberOfTuples());
3039 for(mcIdType i=0;i<nbOfTuples;i++,cptr++)
3040 if(*cptr<vmin || *cptr>vmax)
3041 ret->pushBackSilent(i);
3046 * Returns a new DataArrayDouble by concatenating two given arrays, so that (1) the number
3047 * of tuples in the result array is a sum of the number of tuples of given arrays and (2)
3048 * the number of component in the result array is same as that of each of given arrays.
3049 * Info on components is copied from the first of the given arrays. Number of components
3050 * in the given arrays must be the same.
3051 * \param [in] a1 - an array to include in the result array.
3052 * \param [in] a2 - another array to include in the result array.
3053 * \return DataArrayDouble * - the new instance of DataArrayDouble.
3054 * The caller is to delete this result array using decrRef() as it is no more
3056 * \throw If both \a a1 and \a a2 are NULL.
3057 * \throw If \a a1->getNumberOfComponents() != \a a2->getNumberOfComponents().
3059 DataArrayDouble *DataArrayDouble::Aggregate(const DataArrayDouble *a1, const DataArrayDouble *a2)
3061 std::vector<const DataArrayDouble *> tmp(2);
3062 tmp[0]=a1; tmp[1]=a2;
3063 return Aggregate(tmp);
3067 * Returns a new DataArrayDouble by concatenating all 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 * If the number of non null of elements in \a arr is equal to one the returned object is a copy of it
3073 * not the object itself.
3074 * \param [in] arr - a sequence of arrays to include in the result array.
3075 * \return DataArrayDouble * - the new instance of DataArrayDouble.
3076 * The caller is to delete this result array using decrRef() as it is no more
3078 * \throw If all arrays within \a arr are NULL.
3079 * \throw If getNumberOfComponents() of arrays within \a arr.
3081 DataArrayDouble *DataArrayDouble::Aggregate(const std::vector<const DataArrayDouble *>& arr)
3083 std::vector<const DataArrayDouble *> a;
3084 for(std::vector<const DataArrayDouble *>::const_iterator it4=arr.begin();it4!=arr.end();it4++)
3088 throw INTERP_KERNEL::Exception("DataArrayDouble::Aggregate : input list must contain at least one NON EMPTY DataArrayDouble !");
3089 std::vector<const DataArrayDouble *>::const_iterator it=a.begin();
3090 std::size_t nbOfComp((*it)->getNumberOfComponents());
3091 mcIdType nbt=(*it++)->getNumberOfTuples();
3092 for(mcIdType i=1;it!=a.end();it++,i++)
3094 if((*it)->getNumberOfComponents()!=nbOfComp)
3095 throw INTERP_KERNEL::Exception("DataArrayDouble::Aggregate : Nb of components mismatch for array aggregation !");
3096 nbt+=(*it)->getNumberOfTuples();
3098 MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
3099 ret->alloc(nbt,nbOfComp);
3100 double *pt=ret->getPointer();
3101 for(it=a.begin();it!=a.end();it++)
3102 pt=std::copy((*it)->getConstPointer(),(*it)->getConstPointer()+(*it)->getNbOfElems(),pt);
3103 ret->copyStringInfoFrom(*(a[0]));
3108 * Returns a new DataArrayDouble containing a dot product of two given arrays, so that
3109 * the i-th tuple of the result array is a sum of products of j-th components of i-th
3110 * tuples of given arrays (\f$ a_i = \sum_{j=1}^n a1_j * a2_j \f$).
3111 * Info on components and name is copied from the first of the given arrays.
3112 * Number of tuples and components in the given arrays must be the same.
3113 * \param [in] a1 - a given array.
3114 * \param [in] a2 - another given array.
3115 * \return DataArrayDouble * - the new instance of DataArrayDouble.
3116 * The caller is to delete this result array using decrRef() as it is no more
3118 * \throw If either \a a1 or \a a2 is NULL.
3119 * \throw If any given array is not allocated.
3120 * \throw If \a a1->getNumberOfTuples() != \a a2->getNumberOfTuples()
3121 * \throw If \a a1->getNumberOfComponents() != \a a2->getNumberOfComponents()
3123 DataArrayDouble *DataArrayDouble::Dot(const DataArrayDouble *a1, const DataArrayDouble *a2)
3126 throw INTERP_KERNEL::Exception("DataArrayDouble::Dot : input DataArrayDouble instance is NULL !");
3127 a1->checkAllocated();
3128 a2->checkAllocated();
3129 std::size_t nbOfComp(a1->getNumberOfComponents());
3130 if(nbOfComp!=a2->getNumberOfComponents())
3131 throw INTERP_KERNEL::Exception("Nb of components mismatch for array Dot !");
3132 mcIdType nbOfTuple(a1->getNumberOfTuples());
3133 if(nbOfTuple!=a2->getNumberOfTuples())
3134 throw INTERP_KERNEL::Exception("Nb of tuples mismatch for array Dot !");
3135 DataArrayDouble *ret=DataArrayDouble::New();
3136 ret->alloc(nbOfTuple,1);
3137 double *retPtr=ret->getPointer();
3138 const double *a1Ptr=a1->begin(),*a2Ptr(a2->begin());
3139 for(mcIdType i=0;i<nbOfTuple;i++)
3142 for(std::size_t j=0;j<nbOfComp;j++)
3143 sum+=a1Ptr[i*nbOfComp+j]*a2Ptr[i*nbOfComp+j];
3146 ret->setInfoOnComponent(0,a1->getInfoOnComponent(0));
3147 ret->setName(a1->getName());
3152 * Returns a new DataArrayDouble containing a cross product of two given arrays, so that
3153 * the i-th tuple of the result array contains 3 components of a vector which is a cross
3154 * product of two vectors defined by the i-th tuples of given arrays.
3155 * Info on components is copied from the first of the given arrays.
3156 * Number of tuples in the given arrays must be the same.
3157 * Number of components in the given arrays must be 3.
3158 * \param [in] a1 - a given array.
3159 * \param [in] a2 - another given array.
3160 * \return DataArrayDouble * - the new instance of DataArrayDouble.
3161 * The caller is to delete this result array using decrRef() as it is no more
3163 * \throw If either \a a1 or \a a2 is NULL.
3164 * \throw If \a a1->getNumberOfTuples() != \a a2->getNumberOfTuples()
3165 * \throw If \a a1->getNumberOfComponents() != 3
3166 * \throw If \a a2->getNumberOfComponents() != 3
3168 DataArrayDouble *DataArrayDouble::CrossProduct(const DataArrayDouble *a1, const DataArrayDouble *a2)
3171 throw INTERP_KERNEL::Exception("DataArrayDouble::CrossProduct : input DataArrayDouble instance is NULL !");
3172 std::size_t nbOfComp(a1->getNumberOfComponents());
3173 if(nbOfComp!=a2->getNumberOfComponents())
3174 throw INTERP_KERNEL::Exception("Nb of components mismatch for array crossProduct !");
3176 throw INTERP_KERNEL::Exception("Nb of components must be equal to 3 for array crossProduct !");
3177 mcIdType nbOfTuple(a1->getNumberOfTuples());
3178 if(nbOfTuple!=a2->getNumberOfTuples())
3179 throw INTERP_KERNEL::Exception("Nb of tuples mismatch for array crossProduct !");
3180 DataArrayDouble *ret=DataArrayDouble::New();
3181 ret->alloc(nbOfTuple,3);
3182 double *retPtr=ret->getPointer();
3183 const double *a1Ptr(a1->begin()),*a2Ptr(a2->begin());
3184 for(mcIdType i=0;i<nbOfTuple;i++)
3186 retPtr[3*i]=a1Ptr[3*i+1]*a2Ptr[3*i+2]-a1Ptr[3*i+2]*a2Ptr[3*i+1];
3187 retPtr[3*i+1]=a1Ptr[3*i+2]*a2Ptr[3*i]-a1Ptr[3*i]*a2Ptr[3*i+2];
3188 retPtr[3*i+2]=a1Ptr[3*i]*a2Ptr[3*i+1]-a1Ptr[3*i+1]*a2Ptr[3*i];
3190 ret->copyStringInfoFrom(*a1);
3195 * Returns a new DataArrayDouble containing maximal values of two given arrays.
3196 * Info on components is copied from the first of the given arrays.
3197 * Number of tuples and components in the given arrays must be the same.
3198 * \param [in] a1 - an array to compare values with another one.
3199 * \param [in] a2 - another array to compare values with the first one.
3200 * \return DataArrayDouble * - the new instance of DataArrayDouble.
3201 * The caller is to delete this result array using decrRef() as it is no more
3203 * \throw If either \a a1 or \a a2 is NULL.
3204 * \throw If \a a1->getNumberOfTuples() != \a a2->getNumberOfTuples()
3205 * \throw If \a a1->getNumberOfComponents() != \a a2->getNumberOfComponents()
3207 DataArrayDouble *DataArrayDouble::Max(const DataArrayDouble *a1, const DataArrayDouble *a2)
3210 throw INTERP_KERNEL::Exception("DataArrayDouble::Max : input DataArrayDouble instance is NULL !");
3211 std::size_t nbOfComp(a1->getNumberOfComponents());
3212 if(nbOfComp!=a2->getNumberOfComponents())
3213 throw INTERP_KERNEL::Exception("Nb of components mismatch for array Max !");
3214 mcIdType nbOfTuple(a1->getNumberOfTuples());
3215 if(nbOfTuple!=a2->getNumberOfTuples())
3216 throw INTERP_KERNEL::Exception("Nb of tuples mismatch for array Max !");
3217 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
3218 ret->alloc(nbOfTuple,nbOfComp);
3219 double *retPtr(ret->getPointer());
3220 const double *a1Ptr(a1->begin()),*a2Ptr(a2->begin());
3221 std::size_t nbElem(nbOfTuple*nbOfComp);
3222 for(std::size_t i=0;i<nbElem;i++)
3223 retPtr[i]=std::max(a1Ptr[i],a2Ptr[i]);
3224 ret->copyStringInfoFrom(*a1);
3229 * Returns a new DataArrayDouble containing minimal values of two given arrays.
3230 * Info on components is copied from the first of the given arrays.
3231 * Number of tuples and components in the given arrays must be the same.
3232 * \param [in] a1 - an array to compare values with another one.
3233 * \param [in] a2 - another array to compare values with the first one.
3234 * \return DataArrayDouble * - the new instance of DataArrayDouble.
3235 * The caller is to delete this result array using decrRef() as it is no more
3237 * \throw If either \a a1 or \a a2 is NULL.
3238 * \throw If \a a1->getNumberOfTuples() != \a a2->getNumberOfTuples()
3239 * \throw If \a a1->getNumberOfComponents() != \a a2->getNumberOfComponents()
3241 DataArrayDouble *DataArrayDouble::Min(const DataArrayDouble *a1, const DataArrayDouble *a2)
3244 throw INTERP_KERNEL::Exception("DataArrayDouble::Min : input DataArrayDouble instance is NULL !");
3245 std::size_t nbOfComp(a1->getNumberOfComponents());
3246 if(nbOfComp!=a2->getNumberOfComponents())
3247 throw INTERP_KERNEL::Exception("Nb of components mismatch for array min !");
3248 mcIdType nbOfTuple(a1->getNumberOfTuples());
3249 if(nbOfTuple!=a2->getNumberOfTuples())
3250 throw INTERP_KERNEL::Exception("Nb of tuples mismatch for array min !");
3251 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
3252 ret->alloc(nbOfTuple,nbOfComp);
3253 double *retPtr(ret->getPointer());
3254 const double *a1Ptr(a1->begin()),*a2Ptr(a2->begin());
3255 std::size_t nbElem(nbOfTuple*nbOfComp);
3256 for(std::size_t i=0;i<nbElem;i++)
3257 retPtr[i]=std::min(a1Ptr[i],a2Ptr[i]);
3258 ret->copyStringInfoFrom(*a1);
3263 * Returns a new DataArrayDouble that is the result of pow of two given arrays. There are 3
3266 * \param [in] a1 - an array to pow up.
3267 * \param [in] a2 - another array to sum up.
3268 * \return DataArrayDouble * - the new instance of DataArrayDouble.
3269 * The caller is to delete this result array using decrRef() as it is no more
3271 * \throw If either \a a1 or \a a2 is NULL.
3272 * \throw If \a a1->getNumberOfTuples() != \a a2->getNumberOfTuples()
3273 * \throw If \a a1->getNumberOfComponents() != 1 or \a a2->getNumberOfComponents() != 1.
3274 * \throw If there is a negative value in \a a1.
3276 DataArrayDouble *DataArrayDouble::Pow(const DataArrayDouble *a1, const DataArrayDouble *a2)
3279 throw INTERP_KERNEL::Exception("DataArrayDouble::Pow : at least one of input instances is null !");
3280 mcIdType nbOfTuple=a1->getNumberOfTuples();
3281 mcIdType nbOfTuple2=a2->getNumberOfTuples();
3282 std::size_t nbOfComp=a1->getNumberOfComponents();
3283 std::size_t nbOfComp2=a2->getNumberOfComponents();
3284 if(nbOfTuple!=nbOfTuple2)
3285 throw INTERP_KERNEL::Exception("DataArrayDouble::Pow : number of tuples mismatches !");
3286 if(nbOfComp!=1 || nbOfComp2!=1)
3287 throw INTERP_KERNEL::Exception("DataArrayDouble::Pow : number of components of both arrays must be equal to 1 !");
3288 MCAuto<DataArrayDouble> ret=DataArrayDouble::New(); ret->alloc(nbOfTuple,1);
3289 const double *ptr1(a1->begin()),*ptr2(a2->begin());
3290 double *ptr=ret->getPointer();
3291 for(mcIdType i=0;i<nbOfTuple;i++,ptr1++,ptr2++,ptr++)
3295 *ptr=pow(*ptr1,*ptr2);
3299 std::ostringstream oss; oss << "DataArrayDouble::Pow : on tuple #" << i << " of a1 value is < 0 (" << *ptr1 << ") !";
3300 throw INTERP_KERNEL::Exception(oss.str().c_str());
3307 * Apply pow on values of another DataArrayDouble to values of \a this one.
3309 * \param [in] other - an array to pow to \a this one.
3310 * \throw If \a other is NULL.
3311 * \throw If \a this->getNumberOfTuples() != \a other->getNumberOfTuples()
3312 * \throw If \a this->getNumberOfComponents() != 1 or \a other->getNumberOfComponents() != 1
3313 * \throw If there is a negative value in \a this.
3315 void DataArrayDouble::powEqual(const DataArrayDouble *other)
3318 throw INTERP_KERNEL::Exception("DataArrayDouble::powEqual : input instance is null !");
3319 mcIdType nbOfTuple=getNumberOfTuples();
3320 mcIdType nbOfTuple2=other->getNumberOfTuples();
3321 std::size_t nbOfComp=getNumberOfComponents();
3322 std::size_t nbOfComp2=other->getNumberOfComponents();
3323 if(nbOfTuple!=nbOfTuple2)
3324 throw INTERP_KERNEL::Exception("DataArrayDouble::powEqual : number of tuples mismatches !");
3325 if(nbOfComp!=1 || nbOfComp2!=1)
3326 throw INTERP_KERNEL::Exception("DataArrayDouble::powEqual : number of components of both arrays must be equal to 1 !");
3327 double *ptr=getPointer();
3328 const double *ptrc=other->begin();
3329 for(mcIdType i=0;i<nbOfTuple;i++,ptrc++,ptr++)
3332 *ptr=pow(*ptr,*ptrc);
3335 std::ostringstream oss; oss << "DataArrayDouble::powEqual : on tuple #" << i << " of this value is < 0 (" << *ptr << ") !";
3336 throw INTERP_KERNEL::Exception(oss.str().c_str());
3343 * This method is \b NOT wrapped into python because it can be useful only for performance reasons in C++ context.
3344 * All values in \a this must be 0. or 1. within eps error. 0 means false, 1 means true.
3345 * If an another value than 0 or 1 appear (within eps precision) an INTERP_KERNEL::Exception will be thrown.
3347 * \throw if \a this is not allocated.
3348 * \throw if \a this has not exactly one component.
3350 std::vector<bool> DataArrayDouble::toVectorOfBool(double eps) const
3353 if(getNumberOfComponents()!=1)
3354 throw INTERP_KERNEL::Exception("DataArrayDouble::toVectorOfBool : must be applied on single component array !");
3355 mcIdType nbt(getNumberOfTuples());
3356 std::vector<bool> ret(nbt);
3357 const double *pt(begin());
3358 for(mcIdType i=0;i<nbt;i++)
3362 else if(fabs(pt[i]-1.)<eps)
3366 std::ostringstream oss; oss << "DataArrayDouble::toVectorOfBool : the tuple #" << i << " has value " << pt[i] << " is invalid ! must be 0. or 1. !";
3367 throw INTERP_KERNEL::Exception(oss.str().c_str());
3374 * Useless method for end user. Only for MPI/Corba/File serialsation for multi arrays class.
3377 void DataArrayDouble::getTinySerializationIntInformation(std::vector<mcIdType>& tinyInfo) const
3382 tinyInfo[0]=getNumberOfTuples();
3383 tinyInfo[1]=ToIdType(getNumberOfComponents());
3393 * Useless method for end user. Only for MPI/Corba/File serialsation for multi arrays class.
3396 void DataArrayDouble::getTinySerializationStrInformation(std::vector<std::string>& tinyInfo) const
3400 std::size_t nbOfCompo(getNumberOfComponents());
3401 tinyInfo.resize(nbOfCompo+1);
3402 tinyInfo[0]=getName();
3403 for(std::size_t i=0;i<nbOfCompo;i++)
3404 tinyInfo[i+1]=getInfoOnComponent(i);
3409 tinyInfo[0]=getName();
3414 * Useless method for end user. Only for MPI/Corba/File serialsation for multi arrays class.
3415 * This method returns if a feeding is needed.
3417 bool DataArrayDouble::resizeForUnserialization(const std::vector<mcIdType>& tinyInfoI)
3419 mcIdType nbOfTuple=tinyInfoI[0];
3420 mcIdType nbOfComp=tinyInfoI[1];
3421 if(nbOfTuple!=-1 || nbOfComp!=-1)
3423 alloc(nbOfTuple,nbOfComp);
3430 * Useless method for end user. Only for MPI/Corba/File serialsation for multi arrays class.
3432 void DataArrayDouble::finishUnserialization(const std::vector<mcIdType>& tinyInfoI, const std::vector<std::string>& tinyInfoS)
3434 setName(tinyInfoS[0]);
3437 std::size_t nbOfCompo(getNumberOfComponents());
3438 for(std::size_t i=0;i<nbOfCompo;i++)
3439 setInfoOnComponent(i,tinyInfoS[i+1]);
3444 * Low static method that operates 3D rotation of 'nbNodes' 3D nodes whose coordinates are arranged in \a coordsIn
3445 * around an axe ( \a center, \a vect) and with angle \a angle.
3447 void DataArrayDouble::Rotate3DAlg(const double *center, const double *vect, double angle, mcIdType nbNodes, const double *coordsIn, double *coordsOut)
3449 if(!center || !vect)
3450 throw INTERP_KERNEL::Exception("DataArrayDouble::Rotate3DAlg : null vector in input !");
3451 double sina(sin(angle));
3452 double cosa(cos(angle));
3453 double vectorNorm[3];
3455 double matrixTmp[9];
3456 double norm(sqrt(vect[0]*vect[0]+vect[1]*vect[1]+vect[2]*vect[2]));
3457 if(norm<std::numeric_limits<double>::min())
3458 throw INTERP_KERNEL::Exception("DataArrayDouble::Rotate3DAlg : magnitude of input vector is too close of 0. !");
3459 std::transform(vect,vect+3,vectorNorm,std::bind2nd(std::multiplies<double>(),1/norm));
3460 //rotation matrix computation
3461 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;
3462 matrixTmp[0]=vectorNorm[0]*vectorNorm[0]; matrixTmp[1]=vectorNorm[0]*vectorNorm[1]; matrixTmp[2]=vectorNorm[0]*vectorNorm[2];
3463 matrixTmp[3]=vectorNorm[1]*vectorNorm[0]; matrixTmp[4]=vectorNorm[1]*vectorNorm[1]; matrixTmp[5]=vectorNorm[1]*vectorNorm[2];
3464 matrixTmp[6]=vectorNorm[2]*vectorNorm[0]; matrixTmp[7]=vectorNorm[2]*vectorNorm[1]; matrixTmp[8]=vectorNorm[2]*vectorNorm[2];
3465 std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies<double>(),1-cosa));
3466 std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
3467 matrixTmp[0]=0.; matrixTmp[1]=-vectorNorm[2]; matrixTmp[2]=vectorNorm[1];
3468 matrixTmp[3]=vectorNorm[2]; matrixTmp[4]=0.; matrixTmp[5]=-vectorNorm[0];
3469 matrixTmp[6]=-vectorNorm[1]; matrixTmp[7]=vectorNorm[0]; matrixTmp[8]=0.;
3470 std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies<double>(),sina));
3471 std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
3472 //rotation matrix computed.
3474 for(mcIdType i=0; i<nbNodes; i++)
3476 std::transform(coordsIn+i*3,coordsIn+(i+1)*3,center,tmp,std::minus<double>());
3477 coordsOut[i*3]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+matrix[2]*tmp[2]+center[0];
3478 coordsOut[i*3+1]=matrix[3]*tmp[0]+matrix[4]*tmp[1]+matrix[5]*tmp[2]+center[1];
3479 coordsOut[i*3+2]=matrix[6]*tmp[0]+matrix[7]*tmp[1]+matrix[8]*tmp[2]+center[2];
3483 void DataArrayDouble::Symmetry3DPlane(const double point[3], const double normalVector[3], mcIdType nbNodes, const double *coordsIn, double *coordsOut)
3485 double matrix[9],matrix2[9],matrix3[9];
3486 double vect[3],crossVect[3];
3487 INTERP_KERNEL::orthogonalVect3(normalVector,vect);
3488 crossVect[0]=normalVector[1]*vect[2]-normalVector[2]*vect[1];
3489 crossVect[1]=normalVector[2]*vect[0]-normalVector[0]*vect[2];
3490 crossVect[2]=normalVector[0]*vect[1]-normalVector[1]*vect[0];
3491 double nv(INTERP_KERNEL::norm<3>(vect)),ni(INTERP_KERNEL::norm<3>(normalVector)),nc(INTERP_KERNEL::norm<3>(crossVect));
3492 matrix[0]=vect[0]/nv; matrix[1]=crossVect[0]/nc; matrix[2]=-normalVector[0]/ni;
3493 matrix[3]=vect[1]/nv; matrix[4]=crossVect[1]/nc; matrix[5]=-normalVector[1]/ni;
3494 matrix[6]=vect[2]/nv; matrix[7]=crossVect[2]/nc; matrix[8]=-normalVector[2]/ni;
3495 matrix2[0]=vect[0]/nv; matrix2[1]=vect[1]/nv; matrix2[2]=vect[2]/nv;
3496 matrix2[3]=crossVect[0]/nc; matrix2[4]=crossVect[1]/nc; matrix2[5]=crossVect[2]/nc;
3497 matrix2[6]=normalVector[0]/ni; matrix2[7]=normalVector[1]/ni; matrix2[8]=normalVector[2]/ni;
3498 for(mcIdType i=0;i<3;i++)
3499 for(mcIdType j=0;j<3;j++)
3502 for(mcIdType k=0;k<3;k++)
3503 val+=matrix[3*i+k]*matrix2[3*k+j];
3506 //rotation matrix computed.
3508 for(mcIdType i=0; i<nbNodes; i++)
3510 std::transform(coordsIn+i*3,coordsIn+(i+1)*3,point,tmp,std::minus<double>());
3511 coordsOut[i*3]=matrix3[0]*tmp[0]+matrix3[1]*tmp[1]+matrix3[2]*tmp[2]+point[0];
3512 coordsOut[i*3+1]=matrix3[3]*tmp[0]+matrix3[4]*tmp[1]+matrix3[5]*tmp[2]+point[1];
3513 coordsOut[i*3+2]=matrix3[6]*tmp[0]+matrix3[7]*tmp[1]+matrix3[8]*tmp[2]+point[2];
3517 void DataArrayDouble::GiveBaseForPlane(const double normalVector[3], double baseOfPlane[9])
3519 double vect[3],crossVect[3];
3520 INTERP_KERNEL::orthogonalVect3(normalVector,vect);
3521 crossVect[0]=normalVector[1]*vect[2]-normalVector[2]*vect[1];
3522 crossVect[1]=normalVector[2]*vect[0]-normalVector[0]*vect[2];
3523 crossVect[2]=normalVector[0]*vect[1]-normalVector[1]*vect[0];
3524 double nv(INTERP_KERNEL::norm<3>(vect)),ni(INTERP_KERNEL::norm<3>(normalVector)),nc(INTERP_KERNEL::norm<3>(crossVect));
3525 baseOfPlane[0]=vect[0]/nv; baseOfPlane[1]=vect[1]/nv; baseOfPlane[2]=vect[2]/nv;
3526 baseOfPlane[3]=crossVect[0]/nc; baseOfPlane[4]=crossVect[1]/nc; baseOfPlane[5]=crossVect[2]/nc;
3527 baseOfPlane[6]=normalVector[0]/ni; baseOfPlane[7]=normalVector[1]/ni; baseOfPlane[8]=normalVector[2]/ni;
3531 * \param [in] seg2 : coordinates of input seg2 expected to have spacedim==2
3532 * \param [in] tri3 : coordinates of input tri3 also expected to have spacedim==2
3533 * \param [out] coeffs : the result of integration normalized to 1. along \a seg2 inside tri3 sorted by the node id of \a tri3
3534 * \param [out] length : the length of seg2. That is too say the length of integration
3536 void DataArrayDouble::ComputeIntegralOfSeg2IntoTri3(const double seg2[4], const double tri3[6], double coeffs[3], double& length)
3538 length=INTERP_KERNEL::norme_vecteur(seg2,seg2+2);
3540 INTERP_KERNEL::mid_of_seg2(seg2,seg2+2,mid);
3541 INTERP_KERNEL::barycentric_coords<2>(tri3,mid,coeffs); // integral along seg2 is equal to value at the center of SEG2 !
3545 * Low static method that operates 3D rotation of \a nbNodes 3D nodes whose coordinates are arranged in \a coords
3546 * around the center point \a center and with angle \a angle.
3548 void DataArrayDouble::Rotate2DAlg(const double *center, double angle, mcIdType nbNodes, const double *coordsIn, double *coordsOut)
3550 double cosa=cos(angle);
3551 double sina=sin(angle);
3553 matrix[0]=cosa; matrix[1]=-sina; matrix[2]=sina; matrix[3]=cosa;
3555 for(mcIdType i=0; i<nbNodes; i++)
3557 std::transform(coordsIn+i*2,coordsIn+(i+1)*2,center,tmp,std::minus<double>());
3558 coordsOut[i*2]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+center[0];
3559 coordsOut[i*2+1]=matrix[2]*tmp[0]+matrix[3]*tmp[1]+center[1];
3563 DataArrayDoubleIterator::DataArrayDoubleIterator(DataArrayDouble *da):DataArrayIterator<double>(da)
3567 DataArrayDoubleTuple::DataArrayDoubleTuple(double *pt, std::size_t nbOfComp):DataArrayTuple<double>(pt,nbOfComp)
3572 std::string DataArrayDoubleTuple::repr() const
3574 std::ostringstream oss; oss.precision(17); oss << "(";
3575 for(std::size_t i=0;i<_nb_of_compo-1;i++)
3576 oss << _pt[i] << ", ";
3577 oss << _pt[_nb_of_compo-1] << ")";
3581 double DataArrayDoubleTuple::doubleValue() const
3583 return this->zeValue();
3587 * This method returns a newly allocated instance the caller should dealed with by a MEDCoupling::DataArrayDouble::decrRef.
3588 * This method performs \b no copy of data. The content is only referenced using MEDCoupling::DataArrayDouble::useArray with ownership set to \b false.
3589 * 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
3590 * \b nbOfCompo=1 and \bnbOfTuples==this->_nb_of_elem.
3592 DataArrayDouble *DataArrayDoubleTuple::buildDADouble(std::size_t nbOfTuples, std::size_t nbOfCompo) const
3594 return this->buildDA(nbOfTuples,nbOfCompo);
3598 * Returns a full copy of \a this. For more info on copying data arrays see
3599 * \ref MEDCouplingArrayBasicsCopyDeep.
3600 * \return DataArrayInt * - a new instance of DataArrayInt.
3602 DataArrayInt32 *DataArrayInt32::deepCopy() const
3604 return new DataArrayInt32(*this);
3607 DataArrayInt32Iterator *DataArrayInt32::iterator()
3609 return new DataArrayInt32Iterator(this);
3613 DataArrayInt32Iterator::DataArrayInt32Iterator(DataArrayInt32 *da):DataArrayIterator<Int32>(da)
3617 DataArrayInt32Tuple::DataArrayInt32Tuple(Int32 *pt, std::size_t nbOfComp):DataArrayTuple<Int32>(pt,nbOfComp)
3621 std::string DataArrayInt32Tuple::repr() const
3623 std::ostringstream oss; oss << "(";
3624 for(std::size_t i=0;i<_nb_of_compo-1;i++)
3625 oss << _pt[i] << ", ";
3626 oss << _pt[_nb_of_compo-1] << ")";
3630 Int32 DataArrayInt32Tuple::intValue() const
3632 return this->zeValue();
3636 * This method returns a newly allocated instance the caller should dealed with by a MEDCoupling::DataArrayInt::decrRef.
3637 * This method performs \b no copy of data. The content is only referenced using MEDCoupling::DataArrayInt::useArray with ownership set to \b false.
3638 * 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
3639 * \b nbOfCompo=1 and \bnbOfTuples==this->_nb_of_elem.
3641 DataArrayInt32 *DataArrayInt32Tuple::buildDAInt(std::size_t nbOfTuples, std::size_t nbOfCompo) const
3643 return this->buildDA(nbOfTuples,nbOfCompo);
3646 DataArrayInt64Iterator *DataArrayInt64::iterator()
3648 return new DataArrayInt64Iterator(this);
3652 DataArrayInt64Iterator::DataArrayInt64Iterator(DataArrayInt64 *da):DataArrayIterator<Int64>(da)
3656 DataArrayInt64Tuple::DataArrayInt64Tuple(Int64 *pt, std::size_t nbOfComp):DataArrayTuple<Int64>(pt,nbOfComp)
3660 std::string DataArrayInt64Tuple::repr() const
3662 std::ostringstream oss; oss << "(";
3663 for(std::size_t i=0;i<_nb_of_compo-1;i++)
3664 oss << _pt[i] << ", ";
3665 oss << _pt[_nb_of_compo-1] << ")";
3669 Int64 DataArrayInt64Tuple::intValue() const
3671 return this->zeValue();
3674 DataArrayInt64 *DataArrayInt64Tuple::buildDAInt(std::size_t nbOfTuples, std::size_t nbOfCompo) const
3676 return this->buildDA(nbOfTuples,nbOfCompo);
3680 DataArrayInt64 *DataArrayInt64::deepCopy() const
3682 return new DataArrayInt64(*this);