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<Int32>;
50 template class MEDCoupling::DataArrayTemplateClassic<Int64>;
51 template class MEDCoupling::DataArrayTemplateClassic<double>;
52 template class MEDCoupling::DataArrayTemplateFP<double>;
53 template class MEDCoupling::DataArrayIterator<double>;
54 template class MEDCoupling::DataArrayIterator<mcIdType>;
55 template class MEDCoupling::DataArrayDiscrete<Int32>;
56 template class MEDCoupling::DataArrayDiscreteSigned<Int32>;
57 template class MEDCoupling::DataArrayDiscrete<Int64>;
58 template class MEDCoupling::DataArrayDiscreteSigned<Int64>;
59 template class MEDCoupling::DataArrayTuple<mcIdType>;
60 template class MEDCoupling::DataArrayTuple<double>;
61 template class MEDCoupling::DataArrayTuple<float>;
63 template<mcIdType SPACEDIM>
64 void DataArrayDouble::findCommonTuplesAlg(const double *bbox, mcIdType nbNodes, mcIdType limitNodeId, double prec, DataArrayIdType *c, DataArrayIdType *cI) const
66 const double *coordsPtr=getConstPointer();
67 BBTreePts<SPACEDIM,mcIdType> myTree(bbox,0,0,nbNodes,prec);
68 std::vector<bool> isDone(nbNodes);
69 for(mcIdType i=0;i<nbNodes;i++)
73 std::vector<mcIdType> intersectingElems;
74 myTree.getElementsAroundPoint(coordsPtr+i*SPACEDIM,intersectingElems);
75 if(intersectingElems.size()>1)
77 std::vector<mcIdType> commonNodes;
78 for(std::vector<mcIdType>::const_iterator it=intersectingElems.begin();it!=intersectingElems.end();it++)
82 commonNodes.push_back(*it);
85 if(!commonNodes.empty())
87 cI->pushBackSilent(cI->back()+ToIdType(commonNodes.size())+1);
89 c->insertAtTheEnd(commonNodes.begin(),commonNodes.end());
96 template<mcIdType SPACEDIM>
97 void DataArrayDouble::FindTupleIdsNearTuplesAlg(const BBTreePts<SPACEDIM,mcIdType>& myTree, const double *pos, mcIdType nbOfTuples, double eps,
98 DataArrayIdType *c, DataArrayIdType *cI)
100 for(mcIdType i=0;i<nbOfTuples;i++)
102 std::vector<mcIdType> intersectingElems;
103 myTree.getElementsAroundPoint(pos+i*SPACEDIM,intersectingElems);
104 std::vector<mcIdType> commonNodes;
105 for(std::vector<mcIdType>::const_iterator it=intersectingElems.begin();it!=intersectingElems.end();it++)
106 commonNodes.push_back(*it);
107 cI->pushBackSilent(cI->back()+ToIdType(commonNodes.size()));
108 c->insertAtTheEnd(commonNodes.begin(),commonNodes.end());
112 template<mcIdType SPACEDIM>
113 void DataArrayDouble::FindClosestTupleIdAlg(const BBTreePts<SPACEDIM,mcIdType>& myTree, double dist, const double *pos, mcIdType nbOfTuples, const double *thisPt, mcIdType thisNbOfTuples, mcIdType *res)
115 double distOpt(dist);
116 const double *p(pos);
118 for(mcIdType i=0;i<nbOfTuples;i++,p+=SPACEDIM,r++)
123 double ret=myTree.getElementsAroundPoint2(p,distOpt,elem);
124 if(ret!=std::numeric_limits<double>::max())
126 distOpt=std::max(ret,1e-4);
131 { distOpt=2*distOpt; continue; }
136 mcIdType DataArray::EffectiveCircPerm(mcIdType nbOfShift, mcIdType nbOfTuples)
139 throw INTERP_KERNEL::Exception("DataArray::EffectiveCircPerm : number of tuples is expected to be > 0 !");
142 return nbOfShift%nbOfTuples;
146 mcIdType tmp(-nbOfShift);
148 return nbOfTuples-tmp;
152 std::size_t DataArray::getHeapMemorySizeWithoutChildren() const
154 std::size_t sz1=_name.capacity();
155 std::size_t sz2=_info_on_compo.capacity();
157 for(std::vector<std::string>::const_iterator it=_info_on_compo.begin();it!=_info_on_compo.end();it++)
158 sz3+=(*it).capacity();
162 std::vector<const BigMemoryObject *> DataArray::getDirectChildrenWithNull() const
164 return std::vector<const BigMemoryObject *>();
168 * Sets the attribute \a _name of \a this array.
169 * See \ref MEDCouplingArrayBasicsName "DataArrays infos" for more information.
170 * \param [in] name - new array name
172 void DataArray::setName(const std::string& name)
178 * Copies textual data from an \a other DataArray. The copied data are
179 * - the name attribute,
180 * - the information of components.
182 * For more information on these data see \ref MEDCouplingArrayBasicsName "DataArrays infos".
184 * \param [in] other - another instance of DataArray to copy the textual data from.
185 * \throw If number of components of \a this array differs from that of the \a other.
187 void DataArray::copyStringInfoFrom(const DataArray& other)
189 if(_info_on_compo.size()!=other._info_on_compo.size())
190 throw INTERP_KERNEL::Exception("Size of arrays mismatches on copyStringInfoFrom !");
192 _info_on_compo=other._info_on_compo;
195 void DataArray::copyPartOfStringInfoFrom(const DataArray& other, const std::vector<std::size_t>& compoIds)
197 std::size_t nbOfCompoOth=other.getNumberOfComponents();
198 std::size_t newNbOfCompo=compoIds.size();
199 for(std::size_t i=0;i<newNbOfCompo;i++)
200 if(compoIds[i]>=nbOfCompoOth || compoIds[i]<0)
202 std::ostringstream oss; oss << "Specified component id is out of range (" << compoIds[i] << ") compared with nb of actual components (" << nbOfCompoOth << ")";
203 throw INTERP_KERNEL::Exception(oss.str().c_str());
205 for(std::size_t i=0;i<newNbOfCompo;i++)
206 setInfoOnComponent(i,other.getInfoOnComponent(compoIds[i]));
209 void DataArray::copyPartOfStringInfoFrom2(const std::vector<std::size_t>& compoIds, const DataArray& other)
211 if(compoIds.size()!=other.getNumberOfComponents())
212 throw INTERP_KERNEL::Exception("Given compoIds has not the same size as number of components of given array !");
213 std::size_t partOfCompoToSet=compoIds.size();
214 std::size_t nbOfCompo=getNumberOfComponents();
215 for(std::size_t i=0;i<partOfCompoToSet;i++)
216 if(compoIds[i]>=nbOfCompo || compoIds[i]<0)
218 std::ostringstream oss; oss << "Specified component id is out of range (" << compoIds[i] << ") compared with nb of actual components (" << nbOfCompo << ")";
219 throw INTERP_KERNEL::Exception(oss.str().c_str());
221 for(std::size_t i=0;i<partOfCompoToSet;i++)
222 setInfoOnComponent(compoIds[i],other.getInfoOnComponent(i));
225 bool DataArray::areInfoEqualsIfNotWhy(const DataArray& other, std::string& reason) const
227 std::ostringstream oss;
228 if(_name!=other._name)
230 oss << "Names DataArray mismatch : this name=\"" << _name << " other name=\"" << other._name << "\" !";
234 if(_info_on_compo!=other._info_on_compo)
236 oss << "Components DataArray mismatch : \nThis components=";
237 for(std::vector<std::string>::const_iterator it=_info_on_compo.begin();it!=_info_on_compo.end();it++)
238 oss << "\"" << *it << "\",";
239 oss << "\nOther components=";
240 for(std::vector<std::string>::const_iterator it=other._info_on_compo.begin();it!=other._info_on_compo.end();it++)
241 oss << "\"" << *it << "\",";
249 * Compares textual information of \a this DataArray with that of an \a other one.
250 * The compared data are
251 * - the name attribute,
252 * - the information of components.
254 * For more information on these data see \ref MEDCouplingArrayBasicsName "DataArrays infos".
255 * \param [in] other - another instance of DataArray to compare the textual data of.
256 * \return bool - \a true if the textual information is same, \a false else.
258 bool DataArray::areInfoEquals(const DataArray& other) const
261 return areInfoEqualsIfNotWhy(other,tmp);
264 void DataArray::reprWithoutNameStream(std::ostream& stream) const
266 stream << "Number of components : "<< getNumberOfComponents() << "\n";
267 stream << "Info of these components : ";
268 for(std::vector<std::string>::const_iterator iter=_info_on_compo.begin();iter!=_info_on_compo.end();iter++)
269 stream << "\"" << *iter << "\" ";
273 std::string DataArray::cppRepr(const std::string& varName) const
275 std::ostringstream ret;
276 reprCppStream(varName,ret);
281 * Sets information on all components. To know more on format of this information
282 * see \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
283 * \param [in] info - a vector of strings.
284 * \throw If size of \a info differs from the number of components of \a this.
286 void DataArray::setInfoOnComponents(const std::vector<std::string>& info)
288 if(getNumberOfComponents()!=info.size())
290 std::ostringstream oss; oss << "DataArray::setInfoOnComponents : input is of size " << info.size() << " whereas number of components is equal to " << getNumberOfComponents() << " !";
291 throw INTERP_KERNEL::Exception(oss.str().c_str());
297 * This method is only a dispatcher towards DataArrayDouble::setPartOfValues3, DataArrayInt::setPartOfValues3, DataArrayChar::setPartOfValues3 depending on the true
298 * type of \a this and \a aBase.
300 * \throw If \a aBase and \a this do not have the same type.
302 * \sa DataArrayDouble::setPartOfValues3, DataArrayInt::setPartOfValues3, DataArrayChar::setPartOfValues3.
304 void DataArray::setPartOfValuesBase3(const DataArray *aBase, const mcIdType *bgTuples, const mcIdType *endTuples, mcIdType bgComp, mcIdType endComp, mcIdType stepComp, bool strictCompoCompare)
307 throw INTERP_KERNEL::Exception("DataArray::setPartOfValuesBase3 : input aBase object is NULL !");
308 DataArrayDouble *this1(dynamic_cast<DataArrayDouble *>(this));
309 DataArrayIdType *this2(dynamic_cast<DataArrayIdType *>(this));
310 DataArrayChar *this3(dynamic_cast<DataArrayChar *>(this));
311 const DataArrayDouble *a1(dynamic_cast<const DataArrayDouble *>(aBase));
312 const DataArrayIdType *a2(dynamic_cast<const DataArrayIdType *>(aBase));
313 const DataArrayChar *a3(dynamic_cast<const DataArrayChar *>(aBase));
316 this1->setPartOfValues3(a1,bgTuples,endTuples,bgComp,endComp,stepComp,strictCompoCompare);
321 this2->setPartOfValues3(a2,bgTuples,endTuples,bgComp,endComp,stepComp,strictCompoCompare);
326 this3->setPartOfValues3(a3,bgTuples,endTuples,bgComp,endComp,stepComp,strictCompoCompare);
329 throw INTERP_KERNEL::Exception("DataArray::setPartOfValuesBase3 : input aBase object and this do not have the same type !");
332 std::vector<std::string> DataArray::getVarsOnComponent() const
334 std::size_t nbOfCompo=_info_on_compo.size();
335 std::vector<std::string> ret(nbOfCompo);
336 for(std::size_t i=0;i<nbOfCompo;i++)
337 ret[i]=getVarOnComponent(i);
341 std::vector<std::string> DataArray::getUnitsOnComponent() const
343 std::size_t nbOfCompo=_info_on_compo.size();
344 std::vector<std::string> ret(nbOfCompo);
345 for(std::size_t i=0;i<nbOfCompo;i++)
346 ret[i]=getUnitOnComponent(i);
351 * Returns information on a component specified by an index.
352 * To know more on format of this information
353 * see \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
354 * \param [in] i - the index (zero based) of the component of interest.
355 * \return std::string - a string containing the information on \a i-th component.
356 * \throw If \a i is not a valid component index.
358 std::string DataArray::getInfoOnComponent(std::size_t i) const
360 if(i<_info_on_compo.size())
361 return _info_on_compo[i];
364 std::ostringstream oss; oss << "DataArray::getInfoOnComponent : Specified component id is out of range (" << i << ") compared with nb of actual components (" << _info_on_compo.size();
365 throw INTERP_KERNEL::Exception(oss.str().c_str());
370 * Returns the var part of the full information of the \a i-th component.
371 * For example, if \c getInfoOnComponent(0) returns "SIGXY [N/m^2]", then
372 * \c getVarOnComponent(0) returns "SIGXY".
373 * If a unit part of information is not detected by presence of
374 * two square brackets, then the full information is returned.
375 * To read more about the component information format, see
376 * \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
377 * \param [in] i - the index (zero based) of the component of interest.
378 * \return std::string - a string containing the var information, or the full info.
379 * \throw If \a i is not a valid component index.
381 std::string DataArray::getVarOnComponent(std::size_t i) const
383 if(i<_info_on_compo.size())
385 return GetVarNameFromInfo(_info_on_compo[i]);
389 std::ostringstream oss; oss << "DataArray::getVarOnComponent : Specified component id is out of range (" << i << ") compared with nb of actual components (" << _info_on_compo.size();
390 throw INTERP_KERNEL::Exception(oss.str().c_str());
395 * Returns the unit part of the full information of the \a i-th component.
396 * For example, if \c getInfoOnComponent(0) returns "SIGXY [ N/m^2]", then
397 * \c getUnitOnComponent(0) returns " N/m^2".
398 * If a unit part of information is not detected by presence of
399 * two square brackets, then an empty string is returned.
400 * To read more about the component information format, see
401 * \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
402 * \param [in] i - the index (zero based) of the component of interest.
403 * \return std::string - a string containing the unit information, if any, or "".
404 * \throw If \a i is not a valid component index.
406 std::string DataArray::getUnitOnComponent(std::size_t i) const
408 if(i<_info_on_compo.size())
410 return GetUnitFromInfo(_info_on_compo[i]);
414 std::ostringstream oss; oss << "DataArray::getUnitOnComponent : Specified component id is out of range (" << i << ") compared with nb of actual components (" << _info_on_compo.size();
415 throw INTERP_KERNEL::Exception(oss.str().c_str());
420 * Returns the var part of the full component information.
421 * For example, if \a info == "SIGXY [N/m^2]", then this method returns "SIGXY".
422 * If a unit part of information is not detected by presence of
423 * two square brackets, then the whole \a info is returned.
424 * To read more about the component information format, see
425 * \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
426 * \param [in] info - the full component information.
427 * \return std::string - a string containing only var information, or the \a info.
429 std::string DataArray::GetVarNameFromInfo(const std::string& info)
431 std::size_t p1=info.find_last_of('[');
432 std::size_t p2=info.find_last_of(']');
433 if(p1==std::string::npos || p2==std::string::npos)
438 return std::string();
439 std::size_t p3=info.find_last_not_of(' ',p1-1);
440 return info.substr(0,p3+1);
444 * Returns the unit part of the full component information.
445 * For example, if \a info == "SIGXY [ N/m^2]", then this method returns " N/m^2".
446 * If a unit part of information is not detected by presence of
447 * two square brackets, then an empty string is returned.
448 * To read more about the component information format, see
449 * \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
450 * \param [in] info - the full component information.
451 * \return std::string - a string containing only unit information, if any, or "".
453 std::string DataArray::GetUnitFromInfo(const std::string& info)
455 std::size_t p1=info.find_last_of('[');
456 std::size_t p2=info.find_last_of(']');
457 if(p1==std::string::npos || p2==std::string::npos)
458 return std::string();
460 return std::string();
461 return info.substr(p1+1,p2-p1-1);
465 * This method put in info format the result of the merge of \a var and \a unit.
466 * The standard format for that is "var [unit]".
467 * Inversely you can retrieve the var part or the unit part of info string using resp. GetVarNameFromInfo and GetUnitFromInfo.
469 std::string DataArray::BuildInfoFromVarAndUnit(const std::string& var, const std::string& unit)
471 std::ostringstream oss;
472 oss << var << " [" << unit << "]";
476 std::string DataArray::GetAxisTypeRepr(MEDCouplingAxisType at)
481 return std::string("AX_CART");
483 return std::string("AX_CYL");
485 return std::string("AX_SPHER");
487 throw INTERP_KERNEL::Exception("DataArray::GetAxisTypeRepr : unrecognized axis type enum !");
492 * Returns a new DataArray by concatenating all given arrays, so that (1) the number
493 * of tuples in the result array is a sum of the number of tuples of given arrays and (2)
494 * the number of component in the result array is same as that of each of given arrays.
495 * Info on components is copied from the first of the given arrays. Number of components
496 * in the given arrays must be the same.
497 * \param [in] arrs - a sequence of arrays to include in the result array. All arrays must have the same type.
498 * \return DataArray * - the new instance of DataArray (that can be either DataArrayInt, DataArrayDouble, DataArrayChar).
499 * The caller is to delete this result array using decrRef() as it is no more
501 * \throw If all arrays within \a arrs are NULL.
502 * \throw If all not null arrays in \a arrs have not the same type.
503 * \throw If getNumberOfComponents() of arrays within \a arrs.
505 DataArray *DataArray::Aggregate(const std::vector<const DataArray *>& arrs)
507 std::vector<const DataArray *> arr2;
508 for(std::vector<const DataArray *>::const_iterator it=arrs.begin();it!=arrs.end();it++)
512 throw INTERP_KERNEL::Exception("DataArray::Aggregate : only null instance in input vector !");
513 std::vector<const DataArrayDouble *> arrd;
514 std::vector<const DataArrayIdType *> arri;
515 std::vector<const DataArrayChar *> arrc;
516 for(std::vector<const DataArray *>::const_iterator it=arr2.begin();it!=arr2.end();it++)
518 const DataArrayDouble *a=dynamic_cast<const DataArrayDouble *>(*it);
520 { arrd.push_back(a); continue; }
521 const DataArrayIdType *b=dynamic_cast<const DataArrayIdType *>(*it);
523 { arri.push_back(b); continue; }
524 const DataArrayChar *c=dynamic_cast<const DataArrayChar *>(*it);
526 { arrc.push_back(c); continue; }
527 throw INTERP_KERNEL::Exception("DataArray::Aggregate : presence of not null instance in inuput that is not in [DataArrayDouble, DataArrayInt, DataArrayChar] !");
529 if(arr2.size()==arrd.size())
530 return DataArrayDouble::Aggregate(arrd);
531 if(arr2.size()==arri.size())
532 return DataArrayIdType::Aggregate(arri);
533 if(arr2.size()==arrc.size())
534 return DataArrayChar::Aggregate(arrc);
535 throw INTERP_KERNEL::Exception("DataArray::Aggregate : all input arrays must have the same type !");
539 * Sets information on a component specified by an index.
540 * To know more on format of this information
541 * see \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
542 * \warning Don't pass NULL as \a info!
543 * \param [in] i - the index (zero based) of the component of interest.
544 * \param [in] info - the string containing the information.
545 * \throw If \a i is not a valid component index.
547 void DataArray::setInfoOnComponent(std::size_t i, const std::string& info)
549 if(i<_info_on_compo.size())
550 _info_on_compo[i]=info;
553 std::ostringstream oss; oss << "DataArray::setInfoOnComponent : Specified component id is out of range (" << i << ") compared with nb of actual components (" << _info_on_compo.size();
554 throw INTERP_KERNEL::Exception(oss.str().c_str());
559 * Sets information on all components. This method can change number of components
560 * at certain conditions; if the conditions are not respected, an exception is thrown.
561 * The number of components can be changed in \a this only if \a this is not allocated.
562 * The condition of number of components must not be changed.
564 * To know more on format of the component information see
565 * \ref MEDCouplingArrayBasicsCompoName "DataArrays infos".
566 * \param [in] info - a vector of component infos.
567 * \throw If \a this->getNumberOfComponents() != \a info.size() && \a this->isAllocated()
569 void DataArray::setInfoAndChangeNbOfCompo(const std::vector<std::string>& info)
571 if(getNumberOfComponents()!=info.size())
577 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 !";
578 throw INTERP_KERNEL::Exception(oss.str().c_str());
585 void DataArray::checkNbOfTuples(mcIdType nbOfTuples, const std::string& msg) const
587 if(getNumberOfTuples()!=nbOfTuples)
589 std::ostringstream oss; oss << msg << " : mismatch number of tuples : expected " << nbOfTuples << " having " << getNumberOfTuples() << " !";
590 throw INTERP_KERNEL::Exception(oss.str().c_str());
594 void DataArray::checkNbOfComps(std::size_t nbOfCompo, const std::string& msg) const
596 if (getNumberOfComponents()!=nbOfCompo)
598 std::ostringstream oss; oss << msg << " : mismatch number of components : expected " << nbOfCompo << " having " << getNumberOfComponents() << " !";
599 throw INTERP_KERNEL::Exception(oss.str().c_str());
603 void DataArray::checkNbOfElems(mcIdType nbOfElems, const std::string& msg) const
605 if(getNbOfElems()!=nbOfElems)
607 std::ostringstream oss; oss << msg << " : mismatch number of elems : Expected " << nbOfElems << " having " << getNbOfElems() << " !";
608 throw INTERP_KERNEL::Exception(oss.str().c_str());
612 void DataArray::checkNbOfTuplesAndComp(const DataArray& other, const std::string& msg) const
614 if(getNumberOfTuples()!=other.getNumberOfTuples())
616 std::ostringstream oss; oss << msg << " : mismatch number of tuples : expected " << other.getNumberOfTuples() << " having " << getNumberOfTuples() << " !";
617 throw INTERP_KERNEL::Exception(oss.str().c_str());
619 if(getNumberOfComponents()!=other.getNumberOfComponents())
621 std::ostringstream oss; oss << msg << " : mismatch number of components : expected " << other.getNumberOfComponents() << " having " << getNumberOfComponents() << " !";
622 throw INTERP_KERNEL::Exception(oss.str().c_str());
626 void DataArray::checkNbOfTuplesAndComp(mcIdType nbOfTuples, std::size_t nbOfCompo, const std::string& msg) const
628 checkNbOfTuples(nbOfTuples,msg);
629 checkNbOfComps(nbOfCompo,msg);
633 * Simply this method checks that \b value is in [0,\b ref).
635 void DataArray::CheckValueInRange(mcIdType ref, mcIdType value, const std::string& msg)
637 if(value<0 || value>=ref)
639 std::ostringstream oss; oss << "DataArray::CheckValueInRange : " << msg << " ! Expected in range [0," << ref << "[ having " << value << " !";
640 throw INTERP_KERNEL::Exception(oss.str().c_str());
645 * This method checks that [\b start, \b end) is compliant with ref length \b value.
646 * typically start in [0,\b value) and end in [0,\b value). If value==start and start==end, it is supported.
648 void DataArray::CheckValueInRangeEx(mcIdType value, mcIdType start, mcIdType end, const std::string& msg)
650 if(start<0 || start>=value)
652 if(value!=start || end!=start)
654 std::ostringstream oss; oss << "DataArray::CheckValueInRangeEx : " << msg << " ! Expected start " << start << " of input range, in [0," << value << "[ !";
655 throw INTERP_KERNEL::Exception(oss.str().c_str());
658 if(end<0 || end>value)
660 std::ostringstream oss; oss << "DataArray::CheckValueInRangeEx : " << msg << " ! Expected end " << end << " of input range, in [0," << value << "] !";
661 throw INTERP_KERNEL::Exception(oss.str().c_str());
665 void DataArray::CheckClosingParInRange(mcIdType ref, mcIdType value, const std::string& msg)
667 if(value<0 || value>ref)
669 std::ostringstream oss; oss << "DataArray::CheckClosingParInRange : " << msg << " ! Expected input range in [0," << ref << "] having closing open parenthesis " << value << " !";
670 throw INTERP_KERNEL::Exception(oss.str().c_str());
675 * 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,
676 * typically it is a whole slice of tuples of DataArray or cells, nodes of a mesh...
678 * The input \a sliceId should be an id in [0, \a nbOfSlices) that specifies the slice of work.
680 * \param [in] start - the start of the input slice of the whole work to perform split into slices.
681 * \param [in] stop - the stop of the input slice of the whole work to perform split into slices.
682 * \param [in] step - the step (that can be <0) of the input slice of the whole work to perform split into slices.
683 * \param [in] sliceId - the slice id considered
684 * \param [in] nbOfSlices - the number of slices (typically the number of cores on which the work is expected to be sliced)
685 * \param [out] startSlice - the start of the slice considered
686 * \param [out] stopSlice - the stop of the slice consided
688 * \throw If \a step == 0
689 * \throw If \a nbOfSlices not > 0
690 * \throw If \a sliceId not in [0,nbOfSlices)
692 void DataArray::GetSlice(mcIdType start, mcIdType stop, mcIdType step, mcIdType sliceId, mcIdType nbOfSlices, mcIdType& startSlice, mcIdType& stopSlice)
694 DataArrayTools<mcIdType>::GetSlice(start, stop, step, sliceId, nbOfSlices, startSlice, stopSlice);
697 mcIdType DataArray::GetNumberOfItemGivenBES(mcIdType begin, mcIdType end, mcIdType step, const std::string& msg)
699 return DataArrayTools<mcIdType>::GetNumberOfItemGivenBES(begin, end, step, msg);
702 mcIdType DataArray::GetNumberOfItemGivenBESRelative(mcIdType begin, mcIdType end, mcIdType step, const std::string& msg)
704 return DataArrayTools<mcIdType>::GetNumberOfItemGivenBESRelative(begin, end, step, msg);
707 mcIdType DataArray::GetPosOfItemGivenBESRelativeNoThrow(mcIdType value, mcIdType begin, mcIdType end, mcIdType step)
709 return DataArrayTools<mcIdType>::GetPosOfItemGivenBESRelativeNoThrow(value, begin, end, step);
713 * Returns a new instance of DataArrayDouble. The caller is to delete this array
714 * using decrRef() as it is no more needed.
716 DataArrayDouble *DataArrayDouble::New()
718 return new DataArrayDouble;
722 * Returns the only one value in \a this, if and only if number of elements
723 * (nb of tuples * nb of components) is equal to 1, and that \a this is allocated.
724 * \return double - the sole value stored in \a this array.
725 * \throw If at least one of conditions stated above is not fulfilled.
727 double DataArrayDouble::doubleValue() const
731 if(getNbOfElems()==1)
733 return *getConstPointer();
736 throw INTERP_KERNEL::Exception("DataArrayDouble::doubleValue : DataArrayDouble instance is allocated but number of elements is not equal to 1 !");
739 throw INTERP_KERNEL::Exception("DataArrayDouble::doubleValue : DataArrayDouble instance is not allocated !");
743 * Returns a full copy of \a this. For more info on copying data arrays see
744 * \ref MEDCouplingArrayBasicsCopyDeep.
745 * \return DataArrayDouble * - a new instance of DataArrayDouble. The caller is to
746 * delete this array using decrRef() as it is no more needed.
748 DataArrayDouble *DataArrayDouble::deepCopy() const
750 return new DataArrayDouble(*this);
754 * Checks that \a this array is consistently **increasing** or **decreasing** in value,
755 * with at least absolute difference value of |\a eps| at each step.
756 * If not an exception is thrown.
757 * \param [in] increasing - if \a true, the array values should be increasing.
758 * \param [in] eps - minimal absolute difference between the neighbor values at which
759 * the values are considered different.
760 * \throw If sequence of values is not strictly monotonic in agreement with \a
762 * \throw If \a this->getNumberOfComponents() != 1.
763 * \throw If \a this is not allocated.
765 void DataArrayDouble::checkMonotonic(bool increasing, double eps) const
767 if(!isMonotonic(increasing,eps))
770 throw INTERP_KERNEL::Exception("DataArrayDouble::checkMonotonic : 'this' is not INCREASING monotonic !");
772 throw INTERP_KERNEL::Exception("DataArrayDouble::checkMonotonic : 'this' is not DECREASING monotonic !");
777 * Checks that \a this array is consistently **increasing** or **decreasing** in value,
778 * with at least absolute difference value of |\a eps| at each step.
779 * \param [in] increasing - if \a true, array values should be increasing.
780 * \param [in] eps - minimal absolute difference between the neighbor values at which
781 * the values are considered different.
782 * \return bool - \a true if values change in accordance with \a increasing arg.
783 * \throw If \a this->getNumberOfComponents() != 1.
784 * \throw If \a this is not allocated.
786 bool DataArrayDouble::isMonotonic(bool increasing, double eps) const
789 if(getNumberOfComponents()!=1)
790 throw INTERP_KERNEL::Exception("DataArrayDouble::isMonotonic : only supported with 'this' array with ONE component !");
791 mcIdType nbOfElements(getNumberOfTuples());
792 const double *ptr=getConstPointer();
796 double absEps=fabs(eps);
799 for(mcIdType i=1;i<nbOfElements;i++)
801 if(ptr[i]<(ref+absEps))
809 for(mcIdType i=1;i<nbOfElements;i++)
811 if(ptr[i]>(ref-absEps))
819 void DataArrayDouble::writeVTK(std::ostream& ofs, mcIdType indent, const std::string& nameInFile, DataArrayByte *byteArr) const
821 static const char SPACE[4]={' ',' ',' ',' '};
823 std::string idt(indent,' ');
825 ofs << idt << "<DataArray type=\"Float32\" Name=\"" << nameInFile << "\" NumberOfComponents=\"" << getNumberOfComponents() << "\"";
827 bool areAllEmpty(true);
828 for(std::vector<std::string>::const_iterator it=_info_on_compo.begin();it!=_info_on_compo.end();it++)
832 for(std::size_t i=0;i<_info_on_compo.size();i++)
833 ofs << " ComponentName" << i << "=\"" << _info_on_compo[i] << "\"";
837 ofs << " format=\"appended\" offset=\"" << byteArr->getNumberOfTuples() << "\">";
838 INTERP_KERNEL::AutoPtr<float> tmp(new float[getNbOfElems()]);
840 // to make Visual C++ happy : instead of std::copy(begin(),end(),(float *)tmp);
841 for(const double *src=begin();src!=end();src++,pt++)
843 const char *data(reinterpret_cast<const char *>((float *)tmp));
844 std::size_t sz(getNbOfElems()*sizeof(float));
845 byteArr->insertAtTheEnd(data,data+sz);
846 byteArr->insertAtTheEnd(SPACE,SPACE+4);
850 ofs << " RangeMin=\"" << getMinValueInArray() << "\" RangeMax=\"" << getMaxValueInArray() << "\" format=\"ascii\">\n" << idt;
851 std::copy(begin(),end(),std::ostream_iterator<double>(ofs," "));
853 ofs << std::endl << idt << "</DataArray>\n";
856 void DataArrayDouble::reprCppStream(const std::string& varName, std::ostream& stream) const
858 mcIdType nbTuples=getNumberOfTuples();
859 std::size_t nbComp=getNumberOfComponents();
860 const double *data(getConstPointer());
861 stream.precision(17);
862 stream << "DataArrayDouble *" << varName << "=DataArrayDouble::New();" << std::endl;
863 if(nbTuples*nbComp>=1)
865 stream << "const double " << varName << "Data[" << nbTuples*nbComp << "]={";
866 std::copy(data,data+nbTuples*nbComp-1,std::ostream_iterator<double>(stream,","));
867 stream << data[nbTuples*nbComp-1] << "};" << std::endl;
868 stream << varName << "->useArray(" << varName << "Data,false,CPP_DEALLOC," << nbTuples << "," << nbComp << ");" << std::endl;
871 stream << varName << "->alloc(" << nbTuples << "," << nbComp << ");" << std::endl;
872 stream << varName << "->setName(\"" << getName() << "\");" << std::endl;
876 * Method that gives a quick overvien of \a this for python.
878 void DataArrayDouble::reprQuickOverview(std::ostream& stream) const
880 static const std::size_t MAX_NB_OF_BYTE_IN_REPR=300;
881 stream << "DataArrayDouble C++ instance at " << this << ". ";
884 std::size_t nbOfCompo(_info_on_compo.size());
887 mcIdType nbOfTuples(getNumberOfTuples());
888 stream << "Number of tuples : " << nbOfTuples << ". Number of components : " << nbOfCompo << "." << std::endl;
889 reprQuickOverviewData(stream,MAX_NB_OF_BYTE_IN_REPR);
892 stream << "Number of components : 0.";
895 stream << "*** No data allocated ****";
898 void DataArrayDouble::reprQuickOverviewData(std::ostream& stream, std::size_t maxNbOfByteInRepr) const
900 const double *data=begin();
901 mcIdType nbOfTuples(getNumberOfTuples());
902 std::size_t nbOfCompo(_info_on_compo.size());
903 std::ostringstream oss2; oss2 << "[";
905 std::string oss2Str(oss2.str());
906 bool isFinished=true;
907 for(mcIdType i=0;i<nbOfTuples && isFinished;i++)
912 for(std::size_t j=0;j<nbOfCompo;j++,data++)
915 if(j!=nbOfCompo-1) oss2 << ", ";
921 if(i!=nbOfTuples-1) oss2 << ", ";
922 std::string oss3Str(oss2.str());
923 if(oss3Str.length()<maxNbOfByteInRepr)
935 * Equivalent to DataArrayDouble::isEqual except that if false the reason of
938 * \param [in] other the instance to be compared with \a this
939 * \param [in] prec the precision to compare numeric data of the arrays.
940 * \param [out] reason In case of inequality returns the reason.
941 * \sa DataArrayDouble::isEqual
943 bool DataArrayDouble::isEqualIfNotWhy(const DataArrayDouble& other, double prec, std::string& reason) const
945 if(!areInfoEqualsIfNotWhy(other,reason))
947 return _mem.isEqual(other._mem,prec,reason);
951 * Checks if \a this and another DataArrayDouble are fully equal. For more info see
952 * \ref MEDCouplingArrayBasicsCompare.
953 * \param [in] other - an instance of DataArrayDouble to compare with \a this one.
954 * \param [in] prec - precision value to compare numeric data of the arrays.
955 * \return bool - \a true if the two arrays are equal, \a false else.
957 bool DataArrayDouble::isEqual(const DataArrayDouble& other, double prec) const
960 return isEqualIfNotWhy(other,prec,tmp);
964 * Checks if values of \a this and another DataArrayDouble are equal. For more info see
965 * \ref MEDCouplingArrayBasicsCompare.
966 * \param [in] other - an instance of DataArrayDouble to compare with \a this one.
967 * \param [in] prec - precision value to compare numeric data of the arrays.
968 * \return bool - \a true if the values of two arrays are equal, \a false else.
970 bool DataArrayDouble::isEqualWithoutConsideringStr(const DataArrayDouble& other, double prec) const
973 return _mem.isEqual(other._mem,prec,tmp);
977 * This method checks that all tuples in \a other are in \a this.
978 * If true, the output param \a tupleIds contains the tuples ids of \a this that correspond to tupes in \a this.
979 * 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.
981 * \param [in] other - the array having the same number of components than \a this.
982 * \param [out] tupleIds - the tuple ids containing the same number of tuples than \a other has.
983 * \sa DataArrayDouble::findCommonTuples
985 bool DataArrayDouble::areIncludedInMe(const DataArrayDouble *other, double prec, DataArrayIdType *&tupleIds) const
988 throw INTERP_KERNEL::Exception("DataArrayDouble::areIncludedInMe : input array is NULL !");
989 checkAllocated(); other->checkAllocated();
990 if(getNumberOfComponents()!=other->getNumberOfComponents())
991 throw INTERP_KERNEL::Exception("DataArrayDouble::areIncludedInMe : the number of components does not match !");
992 MCAuto<DataArrayDouble> a=DataArrayDouble::Aggregate(this,other);
993 DataArrayIdType *c=0,*ci=0;
994 a->findCommonTuples(prec,getNumberOfTuples(),c,ci);
995 MCAuto<DataArrayIdType> cSafe(c),ciSafe(ci);
996 mcIdType newNbOfTuples=-1;
997 MCAuto<DataArrayIdType> ids=DataArrayIdType::ConvertIndexArrayToO2N(a->getNumberOfTuples(),c->begin(),ci->begin(),ci->end(),newNbOfTuples);
998 MCAuto<DataArrayIdType> ret1=ids->selectByTupleIdSafeSlice(getNumberOfTuples(),a->getNumberOfTuples(),1);
999 tupleIds=ret1.retn();
1000 return newNbOfTuples==getNumberOfTuples();
1004 * Searches for tuples coincident within \a prec tolerance. Each tuple is considered
1005 * as coordinates of a point in getNumberOfComponents()-dimensional space. The
1006 * distance separating two points is computed with the infinite norm.
1008 * Indices of coincident tuples are stored in output arrays.
1009 * A pair of arrays (\a comm, \a commIndex) is called "Surjective Format 2".
1011 * This method is typically used by MEDCouplingPointSet::findCommonNodes() and
1012 * MEDCouplingUMesh::mergeNodes().
1013 * \param [in] prec - minimal absolute distance between two tuples (infinite norm) at which they are
1014 * considered not coincident.
1015 * \param [in] limitTupleId - limit tuple id. If all tuples within a group of coincident
1016 * tuples have id strictly lower than \a limitTupleId then they are not returned.
1017 * \param [out] comm - the array holding ids (== indices) of coincident tuples.
1018 * \a comm->getNumberOfComponents() == 1.
1019 * \a comm->getNumberOfTuples() == \a commIndex->back().
1020 * \param [out] commIndex - the array dividing all indices stored in \a comm into
1021 * groups of (indices of) coincident tuples. Its every value is a tuple
1022 * index where a next group of tuples begins. For example the second
1023 * group of tuples in \a comm is described by following range of indices:
1024 * [ \a commIndex[1], \a commIndex[2] ). \a commIndex->getNumberOfTuples()-1
1025 * gives the number of groups of coincident tuples.
1026 * \throw If \a this is not allocated.
1027 * \throw If the number of components is not in [1,2,3,4].
1029 * \if ENABLE_EXAMPLES
1030 * \ref cpp_mcdataarraydouble_findcommontuples "Here is a C++ example".
1032 * \ref py_mcdataarraydouble_findcommontuples "Here is a Python example".
1034 * \sa DataArrayInt::ConvertIndexArrayToO2N(), DataArrayDouble::areIncludedInMe
1036 void DataArrayDouble::findCommonTuples(double prec, mcIdType limitTupleId, DataArrayIdType *&comm, DataArrayIdType *&commIndex) const
1039 std::size_t nbOfCompo=getNumberOfComponents();
1040 if ((nbOfCompo<1) || (nbOfCompo>4)) //test before work
1041 throw INTERP_KERNEL::Exception("DataArrayDouble::findCommonTuples : Unexpected spacedim of coords. Must be 1, 2, 3 or 4.");
1043 mcIdType nbOfTuples(getNumberOfTuples());
1045 MCAuto<DataArrayIdType> c(DataArrayIdType::New()),cI(DataArrayIdType::New()); c->alloc(0,1); cI->pushBackSilent(0);
1049 findCommonTuplesAlg<4>(begin(),nbOfTuples,limitTupleId,prec,c,cI);
1052 findCommonTuplesAlg<3>(begin(),nbOfTuples,limitTupleId,prec,c,cI);
1055 findCommonTuplesAlg<2>(begin(),nbOfTuples,limitTupleId,prec,c,cI);
1058 findCommonTuplesAlg<1>(begin(),nbOfTuples,limitTupleId,prec,c,cI);
1061 throw INTERP_KERNEL::Exception("DataArrayDouble::findCommonTuples : nb of components managed are 1,2,3 and 4 ! not implemented for other number of components !");
1064 commIndex=cI.retn();
1068 * This methods returns the minimal distance between the two set of points \a this and \a other.
1069 * So \a this and \a other have to have the same number of components. If not an INTERP_KERNEL::Exception will be thrown.
1070 * This method works only if number of components of \a this (equal to those of \a other) is in 1, 2 or 3.
1072 * \param [out] thisTupleId the tuple id in \a this corresponding to the returned minimal distance
1073 * \param [out] otherTupleId the tuple id in \a other corresponding to the returned minimal distance
1074 * \return the minimal distance between the two set of points \a this and \a other.
1075 * \sa DataArrayDouble::findClosestTupleId
1077 double DataArrayDouble::minimalDistanceTo(const DataArrayDouble *other, mcIdType& thisTupleId, mcIdType& otherTupleId) const
1079 MCAuto<DataArrayIdType> part1=findClosestTupleId(other);
1080 std::size_t nbOfCompo=getNumberOfComponents();
1081 mcIdType otherNbTuples=other->getNumberOfTuples();
1082 const double *thisPt(begin()),*otherPt(other->begin());
1083 const mcIdType *part1Pt(part1->begin());
1084 double ret=std::numeric_limits<double>::max();
1085 for(mcIdType i=0;i<otherNbTuples;i++,part1Pt++,otherPt+=nbOfCompo)
1088 for(std::size_t j=0;j<nbOfCompo;j++)
1089 tmp+=(otherPt[j]-thisPt[nbOfCompo*(*part1Pt)+j])*(otherPt[j]-thisPt[nbOfCompo*(*part1Pt)+j]);
1091 { ret=tmp; thisTupleId=*part1Pt; otherTupleId=i; }
1097 * This methods returns for each tuple in \a other which tuple in \a this is the closest.
1098 * So \a this and \a other have to have the same number of components. If not an INTERP_KERNEL::Exception will be thrown.
1099 * This method works only if number of components of \a this (equal to those of \a other) is in 1, 2 or 3.
1101 * \return a newly allocated (new object to be dealt by the caller) DataArrayInt having \c other->getNumberOfTuples() tuples and one components.
1102 * \sa DataArrayDouble::minimalDistanceTo
1104 DataArrayIdType *DataArrayDouble::findClosestTupleId(const DataArrayDouble *other) const
1107 throw INTERP_KERNEL::Exception("DataArrayDouble::findClosestTupleId : other instance is NULL !");
1108 checkAllocated(); other->checkAllocated();
1109 std::size_t nbOfCompo(getNumberOfComponents());
1110 if(nbOfCompo!=other->getNumberOfComponents())
1112 std::ostringstream oss; oss << "DataArrayDouble::findClosestTupleId : number of components in this is " << nbOfCompo;
1113 oss << ", whereas number of components in other is " << other->getNumberOfComponents() << "! Should be equal !";
1114 throw INTERP_KERNEL::Exception(oss.str().c_str());
1116 mcIdType nbOfTuples(other->getNumberOfTuples());
1117 mcIdType thisNbOfTuples(getNumberOfTuples());
1118 MCAuto<DataArrayIdType> ret=DataArrayIdType::New(); ret->alloc(nbOfTuples,1);
1120 getMinMaxPerComponent(bounds);
1125 double xDelta(fabs(bounds[1]-bounds[0])),yDelta(fabs(bounds[3]-bounds[2])),zDelta(fabs(bounds[5]-bounds[4]));
1126 double delta=std::max(xDelta,yDelta); delta=std::max(delta,zDelta);
1127 double characSize=pow((delta*delta*delta)/((double)thisNbOfTuples),1./3.);
1128 BBTreePts<3,mcIdType> myTree(begin(),0,0,getNumberOfTuples(),characSize*1e-12);
1129 FindClosestTupleIdAlg<3>(myTree,3.*characSize*characSize,other->begin(),nbOfTuples,begin(),thisNbOfTuples,ret->getPointer());
1134 double xDelta(fabs(bounds[1]-bounds[0])),yDelta(fabs(bounds[3]-bounds[2]));
1135 double delta=std::max(xDelta,yDelta);
1136 double characSize=sqrt(delta/(double)thisNbOfTuples);
1137 BBTreePts<2,mcIdType> myTree(begin(),0,0,getNumberOfTuples(),characSize*1e-12);
1138 FindClosestTupleIdAlg<2>(myTree,2.*characSize*characSize,other->begin(),nbOfTuples,begin(),thisNbOfTuples,ret->getPointer());
1143 double characSize=fabs(bounds[1]-bounds[0])/FromIdType<double>(thisNbOfTuples);
1144 BBTreePts<1,mcIdType> myTree(begin(),0,0,getNumberOfTuples(),characSize*1e-12);
1145 FindClosestTupleIdAlg<1>(myTree,1.*characSize*characSize,other->begin(),nbOfTuples,begin(),thisNbOfTuples,ret->getPointer());
1149 throw INTERP_KERNEL::Exception("Unexpected spacedim of coords for findClosestTupleId. Must be 1, 2 or 3.");
1155 * This method expects that \a this and \a otherBBoxFrmt arrays are bounding box arrays ( as the output of MEDCouplingPointSet::getBoundingBoxForBBTree method ).
1156 * 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
1157 * how many bounding boxes in \a otherBBoxFrmt.
1158 * So, this method expects that \a this and \a otherBBoxFrmt have the same number of components.
1160 * \param [in] otherBBoxFrmt - It is an array .
1161 * \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.
1162 * \sa MEDCouplingPointSet::getBoundingBoxForBBTree
1163 * \throw If \a this and \a otherBBoxFrmt have not the same number of components.
1164 * \throw If \a this and \a otherBBoxFrmt number of components is not even (BBox format).
1166 DataArrayIdType *DataArrayDouble::computeNbOfInteractionsWith(const DataArrayDouble *otherBBoxFrmt, double eps) const
1169 throw INTERP_KERNEL::Exception("DataArrayDouble::computeNbOfInteractionsWith : input array is NULL !");
1170 if(!isAllocated() || !otherBBoxFrmt->isAllocated())
1171 throw INTERP_KERNEL::Exception("DataArrayDouble::computeNbOfInteractionsWith : this and input array must be allocated !");
1172 std::size_t nbOfComp(getNumberOfComponents());
1173 mcIdType nbOfTuples(getNumberOfTuples());
1174 if(nbOfComp!=otherBBoxFrmt->getNumberOfComponents())
1176 std::ostringstream oss; oss << "DataArrayDouble::computeNbOfInteractionsWith : this number of components (" << nbOfComp << ") must be equal to the number of components of input array (" << otherBBoxFrmt->getNumberOfComponents() << ") !";
1177 throw INTERP_KERNEL::Exception(oss.str().c_str());
1181 std::ostringstream oss; oss << "DataArrayDouble::computeNbOfInteractionsWith : Number of components (" << nbOfComp << ") is not even ! It should be to be compatible with bbox format !";
1182 throw INTERP_KERNEL::Exception(oss.str().c_str());
1184 MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(nbOfTuples,1);
1185 const double *thisBBPtr(begin());
1186 mcIdType *retPtr(ret->getPointer());
1191 BBTree<3,mcIdType> bbt(otherBBoxFrmt->begin(),0,0,otherBBoxFrmt->getNumberOfTuples(),eps);
1192 for(mcIdType i=0;i<nbOfTuples;i++,retPtr++,thisBBPtr+=nbOfComp)
1193 *retPtr=bbt.getNbOfIntersectingElems(thisBBPtr);
1198 BBTree<2,mcIdType> bbt(otherBBoxFrmt->begin(),0,0,otherBBoxFrmt->getNumberOfTuples(),eps);
1199 for(mcIdType i=0;i<nbOfTuples;i++,retPtr++,thisBBPtr+=nbOfComp)
1200 *retPtr=bbt.getNbOfIntersectingElems(thisBBPtr);
1205 BBTree<1,mcIdType> bbt(otherBBoxFrmt->begin(),0,0,otherBBoxFrmt->getNumberOfTuples(),eps);
1206 for(mcIdType i=0;i<nbOfTuples;i++,retPtr++,thisBBPtr+=nbOfComp)
1207 *retPtr=bbt.getNbOfIntersectingElems(thisBBPtr);
1211 throw INTERP_KERNEL::Exception("DataArrayDouble::computeNbOfInteractionsWith : space dimension supported are [1,2,3] !");
1218 * Returns a copy of \a this array by excluding coincident tuples. Each tuple is
1219 * considered as coordinates of a point in getNumberOfComponents()-dimensional
1220 * space. The distance between tuples is computed using norm2. If several tuples are
1221 * not far each from other than \a prec, only one of them remains in the result
1222 * array. The order of tuples in the result array is same as in \a this one except
1223 * that coincident tuples are excluded.
1224 * \param [in] prec - minimal absolute distance between two tuples at which they are
1225 * considered not coincident.
1226 * \param [in] limitTupleId - limit tuple id. If all tuples within a group of coincident
1227 * tuples have id strictly lower than \a limitTupleId then they are not excluded.
1228 * \return DataArrayDouble * - the new instance of DataArrayDouble that the caller
1229 * is to delete using decrRef() as it is no more needed.
1230 * \throw If \a this is not allocated.
1231 * \throw If the number of components is not in [1,2,3,4].
1233 * \if ENABLE_EXAMPLES
1234 * \ref py_mcdataarraydouble_getdifferentvalues "Here is a Python example".
1237 DataArrayDouble *DataArrayDouble::getDifferentValues(double prec, mcIdType limitTupleId) const
1240 DataArrayIdType *c0=0,*cI0=0;
1241 findCommonTuples(prec,limitTupleId,c0,cI0);
1242 MCAuto<DataArrayIdType> c(c0),cI(cI0);
1243 mcIdType newNbOfTuples=-1;
1244 MCAuto<DataArrayIdType> o2n=DataArrayIdType::ConvertIndexArrayToO2N(getNumberOfTuples(),c0->begin(),cI0->begin(),cI0->end(),newNbOfTuples);
1245 return renumberAndReduce(o2n->getConstPointer(),newNbOfTuples);
1249 * Copy all components in a specified order from another DataArrayDouble.
1250 * Both numerical and textual data is copied. The number of tuples in \a this and
1251 * the other array can be different.
1252 * \param [in] a - the array to copy data from.
1253 * \param [in] compoIds - sequence of zero based indices of components, data of which is
1255 * \throw If \a a is NULL.
1256 * \throw If \a compoIds.size() != \a a->getNumberOfComponents().
1257 * \throw If \a compoIds[i] < 0 or \a compoIds[i] > \a this->getNumberOfComponents().
1259 * \if ENABLE_EXAMPLES
1260 * \ref py_mcdataarraydouble_setselectedcomponents "Here is a Python example".
1263 void DataArrayDouble::setSelectedComponents(const DataArrayDouble *a, const std::vector<std::size_t>& compoIds)
1266 throw INTERP_KERNEL::Exception("DataArrayDouble::setSelectedComponents : input DataArrayDouble is NULL !");
1268 copyPartOfStringInfoFrom2(compoIds,*a);
1269 std::size_t partOfCompoSz=compoIds.size();
1270 std::size_t nbOfCompo=getNumberOfComponents();
1271 mcIdType nbOfTuples=std::min(getNumberOfTuples(),a->getNumberOfTuples());
1272 const double *ac=a->getConstPointer();
1273 double *nc=getPointer();
1274 for(mcIdType i=0;i<nbOfTuples;i++)
1275 for(std::size_t j=0;j<partOfCompoSz;j++,ac++)
1276 nc[nbOfCompo*i+compoIds[j]]=*ac;
1279 * Checks if 0.0 value is present in \a this array. If it is the case, an exception
1281 * \throw If zero is found in \a this array.
1283 void DataArrayDouble::checkNoNullValues() const
1285 const double *tmp=getConstPointer();
1286 mcIdType nbOfElems=getNbOfElems();
1287 const double *where=std::find(tmp,tmp+nbOfElems,0.);
1288 if(where!=tmp+nbOfElems)
1289 throw INTERP_KERNEL::Exception("A value 0.0 have been detected !");
1293 * Computes minimal and maximal value in each component. An output array is filled
1294 * with \c 2 * \a this->getNumberOfComponents() values, so the caller is to allocate
1295 * enough memory before calling this method.
1296 * \param [out] bounds - array of size at least 2 *\a this->getNumberOfComponents().
1297 * It is filled as follows:<br>
1298 * \a bounds[0] = \c min_of_component_0 <br>
1299 * \a bounds[1] = \c max_of_component_0 <br>
1300 * \a bounds[2] = \c min_of_component_1 <br>
1301 * \a bounds[3] = \c max_of_component_1 <br>
1304 void DataArrayDouble::getMinMaxPerComponent(double *bounds) const
1307 std::size_t dim=getNumberOfComponents();
1308 for (std::size_t idim=0; idim<dim; idim++)
1310 bounds[idim*2]=std::numeric_limits<double>::max();
1311 bounds[idim*2+1]=-std::numeric_limits<double>::max();
1313 const double *ptr=getConstPointer();
1314 mcIdType nbOfTuples=getNumberOfTuples();
1315 for(mcIdType i=0;i<nbOfTuples;i++)
1317 for(std::size_t idim=0;idim<dim;idim++)
1319 if(bounds[idim*2]>ptr[i*dim+idim])
1321 bounds[idim*2]=ptr[i*dim+idim];
1323 if(bounds[idim*2+1]<ptr[i*dim+idim])
1325 bounds[idim*2+1]=ptr[i*dim+idim];
1332 * This method retrieves a newly allocated DataArrayDouble instance having same number of tuples than \a this and twice number of components than \a this
1333 * to store both the min and max per component of each tuples.
1334 * \param [in] epsilon the width of the bbox (identical in each direction) - 0.0 by default
1336 * \return a newly created DataArrayDouble instance having \c this->getNumberOfTuples() tuples and 2 * \c this->getNumberOfComponent() components
1338 * \throw If \a this is not allocated yet.
1340 DataArrayDouble *DataArrayDouble::computeBBoxPerTuple(double epsilon) const
1343 const double *dataPtr=getConstPointer();
1344 std::size_t nbOfCompo=getNumberOfComponents();
1345 mcIdType nbTuples=getNumberOfTuples();
1346 MCAuto<DataArrayDouble> bbox=DataArrayDouble::New();
1347 bbox->alloc(nbTuples,2*nbOfCompo);
1348 double *bboxPtr=bbox->getPointer();
1349 for(mcIdType i=0;i<nbTuples;i++)
1351 for(std::size_t j=0;j<nbOfCompo;j++)
1353 bboxPtr[2*nbOfCompo*i+2*j]=dataPtr[nbOfCompo*i+j]-epsilon;
1354 bboxPtr[2*nbOfCompo*i+2*j+1]=dataPtr[nbOfCompo*i+j]+epsilon;
1361 * For each tuples **t** in \a other, this method retrieves tuples in \a this that are equal to **t**.
1362 * Two tuples are considered equal if the euclidian distance between the two tuples is lower than \a eps.
1364 * \param [in] other a DataArrayDouble having same number of components than \a this.
1365 * \param [in] eps absolute precision representing distance (using infinite norm) between 2 tuples behind which 2 tuples are considered equal.
1366 * \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.
1367 * \a cI allows to extract information in \a c.
1368 * \param [out] cI is an indirection array that allows to extract the data contained in \a c.
1370 * \throw In case of:
1371 * - \a this is not allocated
1372 * - \a other is not allocated or null
1373 * - \a this and \a other do not have the same number of components
1374 * - if number of components of \a this is not in [1,2,3]
1376 * \sa MEDCouplingPointSet::getNodeIdsNearPoints, DataArrayDouble::getDifferentValues
1378 void DataArrayDouble::computeTupleIdsNearTuples(const DataArrayDouble *other, double eps, DataArrayIdType *& c, DataArrayIdType *& cI) const
1381 throw INTERP_KERNEL::Exception("DataArrayDouble::computeTupleIdsNearTuples : input pointer other is null !");
1383 other->checkAllocated();
1384 std::size_t nbOfCompo=getNumberOfComponents();
1385 std::size_t otherNbOfCompo=other->getNumberOfComponents();
1386 if(nbOfCompo!=otherNbOfCompo)
1387 throw INTERP_KERNEL::Exception("DataArrayDouble::computeTupleIdsNearTuples : number of components should be equal between this and other !");
1388 mcIdType nbOfTuplesOther=other->getNumberOfTuples();
1389 mcIdType nbOfTuples=getNumberOfTuples();
1390 MCAuto<DataArrayIdType> cArr(DataArrayIdType::New()),cIArr(DataArrayIdType::New()); cArr->alloc(0,1); cIArr->pushBackSilent(0);
1395 BBTreePts<3,mcIdType> myTree(begin(),0,0,nbOfTuples,eps);
1396 FindTupleIdsNearTuplesAlg<3>(myTree,other->getConstPointer(),nbOfTuplesOther,eps,cArr,cIArr);
1401 BBTreePts<2,mcIdType> myTree(begin(),0,0,nbOfTuples,eps);
1402 FindTupleIdsNearTuplesAlg<2>(myTree,other->getConstPointer(),nbOfTuplesOther,eps,cArr,cIArr);
1407 BBTreePts<1,mcIdType> myTree(begin(),0,0,nbOfTuples,eps);
1408 FindTupleIdsNearTuplesAlg<1>(myTree,other->getConstPointer(),nbOfTuplesOther,eps,cArr,cIArr);
1412 throw INTERP_KERNEL::Exception("Unexpected spacedim of coords for computeTupleIdsNearTuples. Must be 1, 2 or 3.");
1414 c=cArr.retn(); cI=cIArr.retn();
1418 * 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
1419 * around origin of 'radius' 1.
1421 * \param [in] eps absolute epsilon. under that value of delta between max and min no scale is performed.
1423 void DataArrayDouble::recenterForMaxPrecision(double eps)
1426 std::size_t dim=getNumberOfComponents();
1427 std::vector<double> bounds(2*dim);
1428 getMinMaxPerComponent(&bounds[0]);
1429 for(std::size_t i=0;i<dim;i++)
1431 double delta=bounds[2*i+1]-bounds[2*i];
1432 double offset=(bounds[2*i]+bounds[2*i+1])/2.;
1434 applyLin(1./delta,-offset/delta,i);
1436 applyLin(1.,-offset,i);
1441 * Returns the maximal value and all its locations within \a this one-dimensional array.
1442 * \param [out] tupleIds - a new instance of DataArrayInt containing indices of
1443 * tuples holding the maximal value. The caller is to delete it using
1444 * decrRef() as it is no more needed.
1445 * \return double - the maximal value among all values of \a this array.
1446 * \throw If \a this->getNumberOfComponents() != 1
1447 * \throw If \a this->getNumberOfTuples() < 1
1449 double DataArrayDouble::getMaxValue2(DataArrayIdType*& tupleIds) const
1453 double ret=getMaxValue(tmp);
1454 tupleIds=findIdsInRange(ret,ret);
1459 * Returns the minimal value and all its locations within \a this one-dimensional array.
1460 * \param [out] tupleIds - a new instance of DataArrayInt containing indices of
1461 * tuples holding the minimal value. The caller is to delete it using
1462 * decrRef() as it is no more needed.
1463 * \return double - the minimal value among all values of \a this array.
1464 * \throw If \a this->getNumberOfComponents() != 1
1465 * \throw If \a this->getNumberOfTuples() < 1
1467 double DataArrayDouble::getMinValue2(DataArrayIdType*& tupleIds) const
1471 double ret=getMinValue(tmp);
1472 tupleIds=findIdsInRange(ret,ret);
1477 * 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.
1478 * This method only works for single component array.
1480 * \return a value in [ 0, \c this->getNumberOfTuples() )
1482 * \throw If \a this is not allocated
1485 mcIdType DataArrayDouble::count(double value, double eps) const
1489 if(getNumberOfComponents()!=1)
1490 throw INTERP_KERNEL::Exception("DataArrayDouble::count : must be applied on DataArrayDouble with only one component, you can call 'rearrange' method before !");
1491 const double *vals=begin();
1492 mcIdType nbOfTuples=getNumberOfTuples();
1493 for(mcIdType i=0;i<nbOfTuples;i++,vals++)
1494 if(fabs(*vals-value)<=eps)
1500 * Returns the average value of \a this one-dimensional array.
1501 * \return double - the average value over all values of \a this array.
1502 * \throw If \a this->getNumberOfComponents() != 1
1503 * \throw If \a this->getNumberOfTuples() < 1
1505 double DataArrayDouble::getAverageValue() const
1507 if(getNumberOfComponents()!=1)
1508 throw INTERP_KERNEL::Exception("DataArrayDouble::getAverageValue : must be applied on DataArrayDouble with only one component, you can call 'rearrange' method before !");
1509 mcIdType nbOfTuples(getNumberOfTuples());
1511 throw INTERP_KERNEL::Exception("DataArrayDouble::getAverageValue : array exists but number of tuples must be > 0 !");
1512 const double *vals=getConstPointer();
1513 double ret=std::accumulate(vals,vals+nbOfTuples,0.);
1514 return ret/FromIdType<double>(nbOfTuples);
1518 * Returns the Euclidean norm of the vector defined by \a this array.
1519 * \return double - the value of the Euclidean norm, i.e.
1520 * the square root of the inner product of vector.
1521 * \throw If \a this is not allocated.
1523 double DataArrayDouble::norm2() const
1527 std::size_t nbOfElems=getNbOfElems();
1528 const double *pt=getConstPointer();
1529 for(std::size_t i=0;i<nbOfElems;i++,pt++)
1535 * Returns the maximum norm of the vector defined by \a this array.
1536 * This method works even if the number of components is different from one.
1537 * If the number of elements in \a this is 0, -1. is returned.
1538 * \return double - the value of the maximum norm, i.e.
1539 * the maximal absolute value among values of \a this array (whatever its number of components).
1540 * \throw If \a this is not allocated.
1542 double DataArrayDouble::normMax() const
1546 std::size_t nbOfElems(getNbOfElems());
1547 const double *pt(getConstPointer());
1548 for(std::size_t i=0;i<nbOfElems;i++,pt++)
1550 double val(std::abs(*pt));
1558 * Returns the maximum norm of for each component of \a this array.
1559 * If the number of elements in \a this is 0, -1. is returned.
1560 * \param [out] res - pointer to an array of result values, of size at least \a
1561 * this->getNumberOfComponents(), that is to be allocated by the caller.
1562 * \throw If \a this is not allocated.
1564 void DataArrayDouble::normMaxPerComponent(double * res) const
1567 mcIdType nbOfTuples(getNumberOfTuples());
1568 std::size_t nbOfCompos(getNumberOfComponents());
1569 std::fill(res, res+nbOfCompos, -1.0);
1570 const double *pt(getConstPointer());
1571 for(mcIdType i=0;i<nbOfTuples;i++)
1572 for (std::size_t j=0; j<nbOfCompos; j++, pt++)
1574 double val(std::abs(*pt));
1582 * Returns the minimum norm (absolute value) of the vector defined by \a this array.
1583 * This method works even if the number of components is different from one.
1584 * If the number of elements in \a this is 0, std::numeric_limits<double>::max() is returned.
1585 * \return double - the value of the minimum norm, i.e.
1586 * the minimal absolute value among values of \a this array (whatever its number of components).
1587 * \throw If \a this is not allocated.
1589 double DataArrayDouble::normMin() const
1592 double ret(std::numeric_limits<double>::max());
1593 std::size_t nbOfElems(getNbOfElems());
1594 const double *pt(getConstPointer());
1595 for(std::size_t i=0;i<nbOfElems;i++,pt++)
1597 double val(std::abs(*pt));
1605 * Accumulates values of each component of \a this array.
1606 * \param [out] res - an array of length \a this->getNumberOfComponents(), allocated
1607 * by the caller, that is filled by this method with sum value for each
1609 * \throw If \a this is not allocated.
1611 void DataArrayDouble::accumulate(double *res) const
1614 const double *ptr=getConstPointer();
1615 mcIdType nbTuple(getNumberOfTuples());
1616 std::size_t nbComps(getNumberOfComponents());
1617 std::fill(res,res+nbComps,0.);
1618 for(mcIdType i=0;i<nbTuple;i++)
1619 std::transform(ptr+i*nbComps,ptr+(i+1)*nbComps,res,res,std::plus<double>());
1623 * This method returns the min distance from an external tuple defined by [ \a tupleBg , \a tupleEnd ) to \a this and
1624 * the first tuple in \a this that matches the returned distance. If there is no tuples in \a this an exception will be thrown.
1627 * \a this is expected to be allocated and expected to have a number of components equal to the distance from \a tupleBg to
1628 * \a tupleEnd. If not an exception will be thrown.
1630 * \param [in] tupleBg start pointer (included) of input external tuple
1631 * \param [in] tupleEnd end pointer (not included) of input external tuple
1632 * \param [out] tupleId the tuple id in \a this that matches the min of distance between \a this and input external tuple
1633 * \return the min distance.
1634 * \sa MEDCouplingUMesh::distanceToPoint
1636 double DataArrayDouble::distanceToTuple(const double *tupleBg, const double *tupleEnd, mcIdType& tupleId) const
1639 mcIdType nbTuple(getNumberOfTuples());
1640 std::size_t nbComps(getNumberOfComponents());
1641 if(nbComps!=(std::size_t)std::distance(tupleBg,tupleEnd))
1642 { 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()); }
1644 throw INTERP_KERNEL::Exception("DataArrayDouble::distanceToTuple : no tuple in this ! No distance to compute !");
1645 double ret0=std::numeric_limits<double>::max();
1647 const double *work=getConstPointer();
1648 for(mcIdType i=0;i<nbTuple;i++)
1651 for(std::size_t j=0;j<nbComps;j++,work++)
1652 val+=(*work-tupleBg[j])*((*work-tupleBg[j]));
1656 { ret0=val; tupleId=i; }
1662 * Accumulate values of the given component of \a this array.
1663 * \param [in] compId - the index of the component of interest.
1664 * \return double - a sum value of \a compId-th component.
1665 * \throw If \a this is not allocated.
1666 * \throw If \a the condition ( 0 <= \a compId < \a this->getNumberOfComponents() ) is
1669 double DataArrayDouble::accumulate(std::size_t compId) const
1672 const double *ptr=getConstPointer();
1673 mcIdType nbTuple(getNumberOfTuples());
1674 std::size_t nbComps(getNumberOfComponents());
1676 throw INTERP_KERNEL::Exception("DataArrayDouble::accumulate : Invalid compId specified : No such nb of components !");
1678 for(mcIdType i=0;i<nbTuple;i++)
1679 ret+=ptr[i*nbComps+compId];
1684 * This method accumulate using addition tuples in \a this using input index array [ \a bgOfIndex, \a endOfIndex ).
1685 * The returned array will have same number of components than \a this and number of tuples equal to
1686 * \c std::distance(bgOfIndex,endOfIndex) \b minus \b one.
1688 * The input index array is expected to be ascendingly sorted in which the all referenced ids should be in [0, \c this->getNumberOfTuples).
1689 * 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.
1691 * \param [in] bgOfIndex - begin (included) of the input index array.
1692 * \param [in] endOfIndex - end (excluded) of the input index array.
1693 * \return DataArrayDouble * - the new instance having the same number of components than \a this.
1695 * \throw If bgOfIndex or end is NULL.
1696 * \throw If input index array is not ascendingly sorted.
1697 * \throw If there is an id in [ \a bgOfIndex, \a endOfIndex ) not in [0, \c this->getNumberOfTuples).
1698 * \throw If std::distance(bgOfIndex,endOfIndex)==0.
1700 DataArrayDouble *DataArrayDouble::accumulatePerChunck(const mcIdType *bgOfIndex, const mcIdType *endOfIndex) const
1702 if(!bgOfIndex || !endOfIndex)
1703 throw INTERP_KERNEL::Exception("DataArrayDouble::accumulatePerChunck : input pointer NULL !");
1705 std::size_t nbCompo(getNumberOfComponents());
1706 mcIdType nbOfTuples(getNumberOfTuples());
1707 std::size_t sz=std::distance(bgOfIndex,endOfIndex);
1709 throw INTERP_KERNEL::Exception("DataArrayDouble::accumulatePerChunck : invalid size of input index array !");
1711 MCAuto<DataArrayDouble> ret=DataArrayDouble::New(); ret->alloc(sz,nbCompo);
1712 const mcIdType *w=bgOfIndex;
1713 if(*w<0 || *w>=nbOfTuples)
1714 throw INTERP_KERNEL::Exception("DataArrayDouble::accumulatePerChunck : The first element of the input index not in [0,nbOfTuples) !");
1715 const double *srcPt=begin()+(*w)*nbCompo;
1716 double *tmp=ret->getPointer();
1717 for(std::size_t i=0;i<sz;i++,tmp+=nbCompo,w++)
1719 std::fill(tmp,tmp+nbCompo,0.);
1722 for(mcIdType j=w[0];j<w[1];j++,srcPt+=nbCompo)
1724 if(j>=0 && j<nbOfTuples)
1725 std::transform(srcPt,srcPt+nbCompo,tmp,tmp,std::plus<double>());
1728 std::ostringstream oss; oss << "DataArrayDouble::accumulatePerChunck : At rank #" << i << " the input index array points to id " << j << " should be in [0," << nbOfTuples << ") !";
1729 throw INTERP_KERNEL::Exception(oss.str().c_str());
1735 std::ostringstream oss; oss << "DataArrayDouble::accumulatePerChunck : At rank #" << i << " the input index array is not in ascendingly sorted.";
1736 throw INTERP_KERNEL::Exception(oss.str().c_str());
1739 ret->copyStringInfoFrom(*this);
1744 * 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.
1745 * This method expects that \a this as only one component. The returned array will have \a this->getNumberOfTuples()+1 tuple with also one component.
1746 * The ith element of returned array is equal to the sum of elements in \a this with rank strictly lower than i.
1748 * \return DataArrayDouble - A newly built array containing cum sum of \a this.
1750 MCAuto<DataArrayDouble> DataArrayDouble::cumSum() const
1753 checkNbOfComps(1,"DataArrayDouble::cumSum : this is expected to be single component");
1754 mcIdType nbOfTuple(getNumberOfTuples());
1755 MCAuto<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(nbOfTuple+1,1);
1756 double *ptr(ret->getPointer());
1758 const double *thisPtr(begin());
1759 for(mcIdType i=0;i<nbOfTuple;i++)
1760 ptr[i+1]=ptr[i]+thisPtr[i];
1765 * Converts each 2D point defined by the tuple of \a this array from the Polar to the
1766 * Cartesian coordinate system. The two components of the tuple of \a this array are
1767 * considered to contain (1) radius and (2) angle of the point in the Polar CS.
1768 * \return DataArrayDouble * - the new instance of DataArrayDouble, whose each tuple
1769 * contains X and Y coordinates of the point in the Cartesian CS. The caller
1770 * is to delete this array using decrRef() as it is no more needed. The array
1771 * does not contain any textual info on components.
1772 * \throw If \a this->getNumberOfComponents() != 2.
1773 * \sa fromCartToPolar
1775 DataArrayDouble *DataArrayDouble::fromPolarToCart() const
1778 std::size_t nbOfComp(getNumberOfComponents());
1780 throw INTERP_KERNEL::Exception("DataArrayDouble::fromPolarToCart : must be an array with exactly 2 components !");
1781 mcIdType nbOfTuple(getNumberOfTuples());
1782 DataArrayDouble *ret(DataArrayDouble::New());
1783 ret->alloc(nbOfTuple,2);
1784 double *w(ret->getPointer());
1785 const double *wIn(getConstPointer());
1786 for(mcIdType i=0;i<nbOfTuple;i++,w+=2,wIn+=2)
1788 w[0]=wIn[0]*cos(wIn[1]);
1789 w[1]=wIn[0]*sin(wIn[1]);
1795 * Converts each 3D point defined by the tuple of \a this array from the Cylindrical to
1796 * the Cartesian coordinate system. The three components of the tuple of \a this array
1797 * are considered to contain (1) radius, (2) azimuth and (3) altitude of the point in
1798 * the Cylindrical CS.
1799 * \return DataArrayDouble * - the new instance of DataArrayDouble, whose each tuple
1800 * contains X, Y and Z coordinates of the point in the Cartesian CS. The info
1801 * on the third component is copied from \a this array. The caller
1802 * is to delete this array using decrRef() as it is no more needed.
1803 * \throw If \a this->getNumberOfComponents() != 3.
1806 DataArrayDouble *DataArrayDouble::fromCylToCart() const
1809 std::size_t nbOfComp(getNumberOfComponents());
1811 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCylToCart : must be an array with exactly 3 components !");
1812 mcIdType nbOfTuple(getNumberOfTuples());
1813 DataArrayDouble *ret(DataArrayDouble::New());
1814 ret->alloc(getNumberOfTuples(),3);
1815 double *w(ret->getPointer());
1816 const double *wIn(getConstPointer());
1817 for(mcIdType i=0;i<nbOfTuple;i++,w+=3,wIn+=3)
1819 w[0]=wIn[0]*cos(wIn[1]);
1820 w[1]=wIn[0]*sin(wIn[1]);
1823 ret->setInfoOnComponent(2,getInfoOnComponent(2));
1828 * Converts each 3D point defined by the tuple of \a this array from the Spherical to
1829 * the Cartesian coordinate system. The three components of the tuple of \a this array
1830 * are considered to contain (1) radius, (2) polar angle and (3) azimuthal angle of the
1831 * point in the Cylindrical CS.
1832 * \return DataArrayDouble * - the new instance of DataArrayDouble, whose each tuple
1833 * contains X, Y and Z coordinates of the point in the Cartesian CS. The info
1834 * on the third component is copied from \a this array. The caller
1835 * is to delete this array using decrRef() as it is no more needed.
1836 * \throw If \a this->getNumberOfComponents() != 3.
1837 * \sa fromCartToSpher
1839 DataArrayDouble *DataArrayDouble::fromSpherToCart() const
1842 std::size_t nbOfComp(getNumberOfComponents());
1844 throw INTERP_KERNEL::Exception("DataArrayDouble::fromSpherToCart : must be an array with exactly 3 components !");
1845 mcIdType nbOfTuple(getNumberOfTuples());
1846 DataArrayDouble *ret(DataArrayDouble::New());
1847 ret->alloc(getNumberOfTuples(),3);
1848 double *w(ret->getPointer());
1849 const double *wIn(getConstPointer());
1850 for(mcIdType i=0;i<nbOfTuple;i++,w+=3,wIn+=3)
1852 w[0]=wIn[0]*cos(wIn[2])*sin(wIn[1]);
1853 w[1]=wIn[0]*sin(wIn[2])*sin(wIn[1]);
1854 w[2]=wIn[0]*cos(wIn[1]);
1860 * 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.
1861 * 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.
1862 * If \a at equals to AX_CYL the returned array will be the result of operation cylindric to cartesian of \a this...
1864 * \param [in] atOfThis - The axis type of \a this.
1865 * \return DataArrayDouble * - the new instance of DataArrayDouble (that must be dealed by caller) containing the result of the cartesianizification of \a this.
1867 DataArrayDouble *DataArrayDouble::cartesianize(MEDCouplingAxisType atOfThis) const
1870 std::size_t nbOfComp(getNumberOfComponents());
1871 MCAuto<DataArrayDouble> ret;
1879 ret=fromCylToCart();
1884 ret=fromPolarToCart();
1888 throw INTERP_KERNEL::Exception("DataArrayDouble::cartesianize : For AX_CYL, number of components must be in [2,3] !");
1892 ret=fromSpherToCart();
1897 ret=fromPolarToCart();
1901 throw INTERP_KERNEL::Exception("DataArrayDouble::cartesianize : For AX_CYL, number of components must be in [2,3] !");
1903 throw INTERP_KERNEL::Exception("DataArrayDouble::cartesianize : not recognized axis type ! Only AX_CART, AX_CYL and AX_SPHER supported !");
1905 ret->copyStringInfoFrom(*this);
1910 * This method returns a newly created array to be deallocated that contains the result of conversion from cartesian to polar.
1911 * This method expects that \a this has exactly 2 components.
1912 * \sa fromPolarToCart
1914 DataArrayDouble *DataArrayDouble::fromCartToPolar() const
1916 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
1918 std::size_t nbOfComp(getNumberOfComponents());
1919 mcIdType nbTuples(getNumberOfTuples());
1921 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToPolar : must be an array with exactly 2 components !");
1922 ret->alloc(nbTuples,2);
1923 double *retPtr(ret->getPointer());
1924 const double *ptr(begin());
1925 for(mcIdType i=0;i<nbTuples;i++,ptr+=2,retPtr+=2)
1927 retPtr[0]=sqrt(ptr[0]*ptr[0]+ptr[1]*ptr[1]);
1928 retPtr[1]=atan2(ptr[1],ptr[0]);
1934 * This method returns a newly created array to be deallocated that contains the result of conversion from cartesian to cylindrical.
1935 * This method expects that \a this has exactly 3 components.
1938 DataArrayDouble *DataArrayDouble::fromCartToCyl() const
1940 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
1942 std::size_t nbOfComp(getNumberOfComponents());
1943 mcIdType nbTuples(getNumberOfTuples());
1945 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCyl : must be an array with exactly 3 components !");
1946 ret->alloc(nbTuples,3);
1947 double *retPtr(ret->getPointer());
1948 const double *ptr(begin());
1949 for(mcIdType i=0;i<nbTuples;i++,ptr+=3,retPtr+=3)
1951 retPtr[0]=sqrt(ptr[0]*ptr[0]+ptr[1]*ptr[1]);
1952 retPtr[1]=atan2(ptr[1],ptr[0]);
1959 * This method returns a newly created array to be deallocated that contains the result of conversion from cartesian to spherical coordinates.
1960 * \sa fromSpherToCart
1962 DataArrayDouble *DataArrayDouble::fromCartToSpher() const
1964 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
1966 std::size_t nbOfComp(getNumberOfComponents());
1967 mcIdType nbTuples(getNumberOfTuples());
1969 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToSpher : must be an array with exactly 3 components !");
1970 ret->alloc(nbTuples,3);
1971 double *retPtr(ret->getPointer());
1972 const double *ptr(begin());
1973 for(mcIdType i=0;i<nbTuples;i++,ptr+=3,retPtr+=3)
1975 retPtr[0]=sqrt(ptr[0]*ptr[0]+ptr[1]*ptr[1]+ptr[2]*ptr[2]);
1976 retPtr[1]=acos(ptr[2]/retPtr[0]);
1977 retPtr[2]=atan2(ptr[1],ptr[0]);
1983 * 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.
1984 * This method expects that \a this has exactly 3 components.
1985 * \sa MEDCouplingFieldDouble::computeVectorFieldCyl
1987 DataArrayDouble *DataArrayDouble::fromCartToCylGiven(const DataArrayDouble *coords, const double center[3], const double vect[3]) const
1990 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCylGiven : input coords are NULL !");
1991 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
1992 checkAllocated(); coords->checkAllocated();
1993 std::size_t nbOfComp(getNumberOfComponents());
1994 mcIdType nbTuples(getNumberOfTuples());
1996 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCylGiven : must be an array with exactly 3 components !");
1997 if(coords->getNumberOfComponents()!=3)
1998 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCylGiven : coords array must have exactly 3 components !");
1999 if(coords->getNumberOfTuples()!=nbTuples)
2000 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCylGiven : coords array must have the same number of tuples !");
2001 ret->alloc(nbTuples,nbOfComp);
2002 double magOfVect(sqrt(vect[0]*vect[0]+vect[1]*vect[1]+vect[2]*vect[2]));
2004 throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCylGiven : magnitude of vect is too low !");
2005 double Ur[3],Uteta[3],Uz[3],*retPtr(ret->getPointer());
2006 const double *coo(coords->begin()),*vectField(begin());
2007 std::transform(vect,vect+3,Uz,std::bind2nd(std::multiplies<double>(),1./magOfVect));
2008 for(mcIdType i=0;i<nbTuples;i++,vectField+=3,retPtr+=3,coo+=3)
2010 std::transform(coo,coo+3,center,Ur,std::minus<double>());
2011 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];
2012 double magOfTeta(sqrt(Uteta[0]*Uteta[0]+Uteta[1]*Uteta[1]+Uteta[2]*Uteta[2]));
2013 std::transform(Uteta,Uteta+3,Uteta,std::bind2nd(std::multiplies<double>(),1./magOfTeta));
2014 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];
2015 retPtr[0]=Ur[0]*vectField[0]+Ur[1]*vectField[1]+Ur[2]*vectField[2];
2016 retPtr[1]=Uteta[0]*vectField[0]+Uteta[1]*vectField[1]+Uteta[2]*vectField[2];
2017 retPtr[2]=Uz[0]*vectField[0]+Uz[1]*vectField[1]+Uz[2]*vectField[2];
2019 ret->copyStringInfoFrom(*this);
2024 * Computes the doubly contracted product of every tensor defined by the tuple of \a this
2025 * array containing 6 components.
2026 * \return DataArrayDouble * - the new instance of DataArrayDouble, whose each tuple
2027 * is calculated from the tuple <em>(t)</em> of \a this array as follows:
2028 * \f$ t[0]^2+t[1]^2+t[2]^2+2*t[3]^2+2*t[4]^2+2*t[5]^2\f$.
2029 * The caller is to delete this result array using decrRef() as it is no more needed.
2030 * \throw If \a this->getNumberOfComponents() != 6.
2032 DataArrayDouble *DataArrayDouble::doublyContractedProduct() const
2035 std::size_t nbOfComp(getNumberOfComponents());
2037 throw INTERP_KERNEL::Exception("DataArrayDouble::doublyContractedProduct : must be an array with exactly 6 components !");
2038 DataArrayDouble *ret=DataArrayDouble::New();
2039 mcIdType nbOfTuple=getNumberOfTuples();
2040 ret->alloc(nbOfTuple,1);
2041 const double *src=getConstPointer();
2042 double *dest=ret->getPointer();
2043 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=6)
2044 *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];
2049 * Computes the determinant of every square matrix defined by the tuple of \a this
2050 * array, which contains either 4, 6 or 9 components. The case of 6 components
2051 * corresponds to that of the upper triangular matrix.
2052 * \return DataArrayDouble * - the new instance of DataArrayDouble, whose each tuple
2053 * is the determinant of matrix of the corresponding tuple of \a this array.
2054 * The caller is to delete this result array using decrRef() as it is no more
2056 * \throw If \a this->getNumberOfComponents() is not in [4,6,9].
2058 DataArrayDouble *DataArrayDouble::determinant() const
2061 DataArrayDouble *ret=DataArrayDouble::New();
2062 mcIdType nbOfTuple=getNumberOfTuples();
2063 ret->alloc(nbOfTuple,1);
2064 const double *src=getConstPointer();
2065 double *dest=ret->getPointer();
2066 switch(getNumberOfComponents())
2069 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=6)
2070 *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];
2073 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=4)
2074 *dest=src[0]*src[3]-src[1]*src[2];
2077 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=9)
2078 *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];
2082 throw INTERP_KERNEL::Exception("DataArrayDouble::determinant : Invalid number of components ! must be in 4,6,9 !");
2087 * Computes 3 eigenvalues of every upper triangular matrix defined by the tuple of
2088 * \a this array, which contains 6 components.
2089 * \return DataArrayDouble * - the new instance of DataArrayDouble containing 3
2090 * components, whose each tuple contains the eigenvalues of the matrix of
2091 * corresponding tuple of \a this array.
2092 * The caller is to delete this result array using decrRef() as it is no more
2094 * \throw If \a this->getNumberOfComponents() != 6.
2096 DataArrayDouble *DataArrayDouble::eigenValues() const
2099 std::size_t nbOfComp=getNumberOfComponents();
2101 throw INTERP_KERNEL::Exception("DataArrayDouble::eigenValues : must be an array with exactly 6 components !");
2102 DataArrayDouble *ret=DataArrayDouble::New();
2103 mcIdType nbOfTuple=getNumberOfTuples();
2104 ret->alloc(nbOfTuple,3);
2105 const double *src=getConstPointer();
2106 double *dest=ret->getPointer();
2107 for(mcIdType i=0;i<nbOfTuple;i++,dest+=3,src+=6)
2108 INTERP_KERNEL::computeEigenValues6(src,dest);
2113 * Computes 3 eigenvectors of every upper triangular matrix defined by the tuple of
2114 * \a this array, which contains 6 components.
2115 * \return DataArrayDouble * - the new instance of DataArrayDouble containing 9
2116 * components, whose each tuple contains 3 eigenvectors of the matrix of
2117 * corresponding tuple of \a this array.
2118 * The caller is to delete this result array using decrRef() as it is no more
2120 * \throw If \a this->getNumberOfComponents() != 6.
2122 DataArrayDouble *DataArrayDouble::eigenVectors() const
2125 std::size_t nbOfComp=getNumberOfComponents();
2127 throw INTERP_KERNEL::Exception("DataArrayDouble::eigenVectors : must be an array with exactly 6 components !");
2128 DataArrayDouble *ret=DataArrayDouble::New();
2129 mcIdType nbOfTuple=getNumberOfTuples();
2130 ret->alloc(nbOfTuple,9);
2131 const double *src=getConstPointer();
2132 double *dest=ret->getPointer();
2133 for(mcIdType i=0;i<nbOfTuple;i++,src+=6)
2136 INTERP_KERNEL::computeEigenValues6(src,tmp);
2137 for(mcIdType j=0;j<3;j++,dest+=3)
2138 INTERP_KERNEL::computeEigenVectorForEigenValue6(src,tmp[j],1e-12,dest);
2144 * Computes the inverse matrix of every matrix defined by the tuple of \a this
2145 * array, which contains either 4, 6 or 9 components. The case of 6 components
2146 * corresponds to that of the upper triangular matrix.
2147 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2148 * same number of components as \a this one, whose each tuple is the inverse
2149 * matrix of the matrix of corresponding tuple of \a this array.
2150 * The caller is to delete this result array using decrRef() as it is no more
2152 * \throw If \a this->getNumberOfComponents() is not in [4,6,9].
2154 DataArrayDouble *DataArrayDouble::inverse() const
2157 std::size_t nbOfComp=getNumberOfComponents();
2158 if(nbOfComp!=6 && nbOfComp!=9 && nbOfComp!=4)
2159 throw INTERP_KERNEL::Exception("DataArrayDouble::inversion : must be an array with 4,6 or 9 components !");
2160 DataArrayDouble *ret=DataArrayDouble::New();
2161 mcIdType nbOfTuple=getNumberOfTuples();
2162 ret->alloc(nbOfTuple,nbOfComp);
2163 const double *src=getConstPointer();
2164 double *dest=ret->getPointer();
2166 for(mcIdType i=0;i<nbOfTuple;i++,dest+=6,src+=6)
2168 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];
2169 dest[0]=(src[1]*src[2]-src[4]*src[4])/det;
2170 dest[1]=(src[0]*src[2]-src[5]*src[5])/det;
2171 dest[2]=(src[0]*src[1]-src[3]*src[3])/det;
2172 dest[3]=(src[5]*src[4]-src[3]*src[2])/det;
2173 dest[4]=(src[5]*src[3]-src[0]*src[4])/det;
2174 dest[5]=(src[3]*src[4]-src[1]*src[5])/det;
2176 else if(nbOfComp==4)
2177 for(mcIdType i=0;i<nbOfTuple;i++,dest+=4,src+=4)
2179 double det=src[0]*src[3]-src[1]*src[2];
2181 dest[1]=-src[1]/det;
2182 dest[2]=-src[2]/det;
2186 for(mcIdType i=0;i<nbOfTuple;i++,dest+=9,src+=9)
2188 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];
2189 dest[0]=(src[4]*src[8]-src[7]*src[5])/det;
2190 dest[1]=(src[7]*src[2]-src[1]*src[8])/det;
2191 dest[2]=(src[1]*src[5]-src[4]*src[2])/det;
2192 dest[3]=(src[6]*src[5]-src[3]*src[8])/det;
2193 dest[4]=(src[0]*src[8]-src[6]*src[2])/det;
2194 dest[5]=(src[2]*src[3]-src[0]*src[5])/det;
2195 dest[6]=(src[3]*src[7]-src[6]*src[4])/det;
2196 dest[7]=(src[6]*src[1]-src[0]*src[7])/det;
2197 dest[8]=(src[0]*src[4]-src[1]*src[3])/det;
2203 * Computes the trace of every matrix defined by the tuple of \a this
2204 * array, which contains either 4, 6 or 9 components. The case of 6 components
2205 * corresponds to that of the upper triangular matrix.
2206 * \return DataArrayDouble * - the new instance of DataArrayDouble containing
2207 * 1 component, whose each tuple is the trace of
2208 * the matrix of corresponding tuple of \a this array.
2209 * The caller is to delete this result array using decrRef() as it is no more
2211 * \throw If \a this->getNumberOfComponents() is not in [4,6,9].
2213 DataArrayDouble *DataArrayDouble::trace() const
2216 std::size_t nbOfComp=getNumberOfComponents();
2217 if(nbOfComp!=6 && nbOfComp!=9 && nbOfComp!=4)
2218 throw INTERP_KERNEL::Exception("DataArrayDouble::trace : must be an array with 4,6 or 9 components !");
2219 DataArrayDouble *ret=DataArrayDouble::New();
2220 mcIdType nbOfTuple=getNumberOfTuples();
2221 ret->alloc(nbOfTuple,1);
2222 const double *src=getConstPointer();
2223 double *dest=ret->getPointer();
2225 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=6)
2226 *dest=src[0]+src[1]+src[2];
2227 else if(nbOfComp==4)
2228 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=4)
2229 *dest=src[0]+src[3];
2231 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=9)
2232 *dest=src[0]+src[4]+src[8];
2237 * Computes the stress deviator tensor of every stress tensor defined by the tuple of
2238 * \a this array, which contains 6 components.
2239 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2240 * same number of components and tuples as \a this array.
2241 * The caller is to delete this result array using decrRef() as it is no more
2243 * \throw If \a this->getNumberOfComponents() != 6.
2245 DataArrayDouble *DataArrayDouble::deviator() const
2248 std::size_t nbOfComp=getNumberOfComponents();
2250 throw INTERP_KERNEL::Exception("DataArrayDouble::deviator : must be an array with exactly 6 components !");
2251 DataArrayDouble *ret=DataArrayDouble::New();
2252 mcIdType nbOfTuple=getNumberOfTuples();
2253 ret->alloc(nbOfTuple,6);
2254 const double *src=getConstPointer();
2255 double *dest=ret->getPointer();
2256 for(mcIdType i=0;i<nbOfTuple;i++,dest+=6,src+=6)
2258 double tr=(src[0]+src[1]+src[2])/3.;
2270 * Computes the magnitude of every vector defined by the tuple of
2272 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2273 * same number of tuples as \a this array and one component.
2274 * The caller is to delete this result array using decrRef() as it is no more
2276 * \throw If \a this is not allocated.
2278 DataArrayDouble *DataArrayDouble::magnitude() const
2281 std::size_t nbOfComp=getNumberOfComponents();
2282 DataArrayDouble *ret=DataArrayDouble::New();
2283 mcIdType nbOfTuple=getNumberOfTuples();
2284 ret->alloc(nbOfTuple,1);
2285 const double *src=getConstPointer();
2286 double *dest=ret->getPointer();
2287 for(mcIdType i=0;i<nbOfTuple;i++,dest++)
2290 for(std::size_t j=0;j<nbOfComp;j++,src++)
2298 * Computes the maximal value within every tuple of \a this array.
2299 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2300 * same number of tuples as \a this array and one component.
2301 * The caller is to delete this result array using decrRef() as it is no more
2303 * \throw If \a this is not allocated.
2304 * \sa DataArrayDouble::maxPerTupleWithCompoId
2306 DataArrayDouble *DataArrayDouble::maxPerTuple() const
2309 std::size_t nbOfComp(getNumberOfComponents());
2310 MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
2311 mcIdType nbOfTuple(getNumberOfTuples());
2312 ret->alloc(nbOfTuple,1);
2313 const double *src=getConstPointer();
2314 double *dest=ret->getPointer();
2315 for(mcIdType i=0;i<nbOfTuple;i++,dest++,src+=nbOfComp)
2316 *dest=*std::max_element(src,src+nbOfComp);
2321 * Computes the maximal value within every tuple of \a this array and it returns the first component
2322 * id for each tuple that corresponds to the maximal value within the tuple.
2324 * \param [out] compoIdOfMaxPerTuple - the new new instance of DataArrayInt containing the
2325 * same number of tuples and only one component.
2326 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2327 * same number of tuples as \a this array and one component.
2328 * The caller is to delete this result array using decrRef() as it is no more
2330 * \throw If \a this is not allocated.
2331 * \sa DataArrayDouble::maxPerTuple
2333 DataArrayDouble *DataArrayDouble::maxPerTupleWithCompoId(DataArrayIdType* &compoIdOfMaxPerTuple) const
2336 std::size_t nbOfComp(getNumberOfComponents());
2337 MCAuto<DataArrayDouble> ret0=DataArrayDouble::New();
2338 MCAuto<DataArrayIdType> ret1=DataArrayIdType::New();
2339 mcIdType nbOfTuple=getNumberOfTuples();
2340 ret0->alloc(nbOfTuple,1); ret1->alloc(nbOfTuple,1);
2341 const double *src=getConstPointer();
2342 double *dest=ret0->getPointer(); mcIdType *dest1=ret1->getPointer();
2343 for(mcIdType i=0;i<nbOfTuple;i++,dest++,dest1++,src+=nbOfComp)
2345 const double *loc=std::max_element(src,src+nbOfComp);
2347 *dest1=ToIdType(std::distance(src,loc));
2349 compoIdOfMaxPerTuple=ret1.retn();
2354 * This method returns a newly allocated DataArrayDouble instance having one component and \c this->getNumberOfTuples() * \c this->getNumberOfTuples() tuples.
2355 * \n This returned array contains the euclidian distance for each tuple in \a this.
2356 * \n So the returned array can be seen as a dense symmetrical matrix whose diagonal elements are equal to 0.
2357 * \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)
2359 * \warning use this method with care because it can leads to big amount of consumed memory !
2361 * \return A newly allocated (huge) MEDCoupling::DataArrayDouble instance that the caller should deal with.
2363 * \throw If \a this is not allocated.
2365 * \sa DataArrayDouble::buildEuclidianDistanceDenseMatrixWith
2367 DataArrayDouble *DataArrayDouble::buildEuclidianDistanceDenseMatrix() const
2370 std::size_t nbOfComp(getNumberOfComponents());
2371 mcIdType nbOfTuples(getNumberOfTuples());
2372 const double *inData=getConstPointer();
2373 MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
2374 ret->alloc(nbOfTuples*nbOfTuples,1);
2375 double *outData=ret->getPointer();
2376 for(mcIdType i=0;i<nbOfTuples;i++)
2378 outData[i*nbOfTuples+i]=0.;
2379 for(mcIdType j=i+1;j<nbOfTuples;j++)
2382 for(std::size_t k=0;k<nbOfComp;k++)
2383 { double delta=inData[i*nbOfComp+k]-inData[j*nbOfComp+k]; dist+=delta*delta; }
2385 outData[i*nbOfTuples+j]=dist;
2386 outData[j*nbOfTuples+i]=dist;
2393 * This method returns a newly allocated DataArrayDouble instance having one component and \c this->getNumberOfTuples() * \c other->getNumberOfTuples() tuples.
2394 * \n This returned array contains the euclidian distance for each tuple in \a other with each tuple in \a this.
2395 * \n So the returned array can be seen as a dense rectangular matrix with \c other->getNumberOfTuples() rows and \c this->getNumberOfTuples() columns.
2396 * \n Output rectangular matrix is sorted along rows.
2397 * \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)
2399 * \warning use this method with care because it can leads to big amount of consumed memory !
2401 * \param [in] other DataArrayDouble instance having same number of components than \a this.
2402 * \return A newly allocated (huge) MEDCoupling::DataArrayDouble instance that the caller should deal with.
2404 * \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.
2406 * \sa DataArrayDouble::buildEuclidianDistanceDenseMatrix
2408 DataArrayDouble *DataArrayDouble::buildEuclidianDistanceDenseMatrixWith(const DataArrayDouble *other) const
2411 throw INTERP_KERNEL::Exception("DataArrayDouble::buildEuclidianDistanceDenseMatrixWith : input parameter is null !");
2413 other->checkAllocated();
2414 std::size_t nbOfComp(getNumberOfComponents());
2415 std::size_t otherNbOfComp(other->getNumberOfComponents());
2416 if(nbOfComp!=otherNbOfComp)
2418 std::ostringstream oss; oss << "DataArrayDouble::buildEuclidianDistanceDenseMatrixWith : this nb of compo=" << nbOfComp << " and other nb of compo=" << otherNbOfComp << ". It should match !";
2419 throw INTERP_KERNEL::Exception(oss.str().c_str());
2421 mcIdType nbOfTuples(getNumberOfTuples());
2422 mcIdType otherNbOfTuples(other->getNumberOfTuples());
2423 const double *inData=getConstPointer();
2424 const double *inDataOther=other->getConstPointer();
2425 MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
2426 ret->alloc(otherNbOfTuples*nbOfTuples,1);
2427 double *outData=ret->getPointer();
2428 for(mcIdType i=0;i<otherNbOfTuples;i++,inDataOther+=nbOfComp)
2430 for(mcIdType j=0;j<nbOfTuples;j++)
2433 for(std::size_t k=0;k<nbOfComp;k++)
2434 { double delta=inDataOther[k]-inData[j*nbOfComp+k]; dist+=delta*delta; }
2436 outData[i*nbOfTuples+j]=dist;
2443 * This method expects that \a this stores 3 tuples containing 2 components each.
2444 * Each of this tuples represent a point into 2D space.
2445 * 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).
2446 * 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[.
2448 * \throw If \a this is not allocated.
2449 * \throw If \a this has not 3 tuples of 2 components
2450 * \throw If tuples/points in \a this are aligned
2452 void DataArrayDouble::asArcOfCircle(double center[2], double& radius, double& ang) const
2455 INTERP_KERNEL::QuadraticPlanarPrecision arcPrec(1e-14);
2456 if(getNumberOfTuples()!=3 && getNumberOfComponents()!=2)
2457 throw INTERP_KERNEL::Exception("DataArrayDouble::asArcCircle : this method expects");
2458 const double *pt(begin());
2459 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]));
2461 INTERP_KERNEL::AutoCppPtr<INTERP_KERNEL::EdgeLin> e1(new INTERP_KERNEL::EdgeLin(n0,n2)),e2(new INTERP_KERNEL::EdgeLin(n2,n1));
2462 INTERP_KERNEL::SegSegIntersector inters(*e1,*e2);
2463 bool colinearity(inters.areColinears());
2465 throw INTERP_KERNEL::Exception("DataArrayDouble::asArcOfCircle : 3 points in this have been detected as colinear !");
2467 INTERP_KERNEL::AutoCppPtr<INTERP_KERNEL::EdgeArcCircle> ret(new INTERP_KERNEL::EdgeArcCircle(n0,n2,n1));
2468 const double *c(ret->getCenter());
2469 center[0]=c[0]; center[1]=c[1];
2470 radius=ret->getRadius();
2471 ang=ret->getAngle();
2475 * Sorts value within every tuple of \a this array.
2476 * \param [in] asc - if \a true, the values are sorted in ascending order, else,
2477 * in descending order.
2478 * \throw If \a this is not allocated.
2480 void DataArrayDouble::sortPerTuple(bool asc)
2483 double *pt=getPointer();
2484 mcIdType nbOfTuple(getNumberOfTuples());
2485 std::size_t nbOfComp(getNumberOfComponents());
2487 for(mcIdType i=0;i<nbOfTuple;i++,pt+=nbOfComp)
2488 std::sort(pt,pt+nbOfComp);
2490 for(mcIdType i=0;i<nbOfTuple;i++,pt+=nbOfComp)
2491 std::sort(pt,pt+nbOfComp,std::greater<double>());
2496 * Modify all elements of \a this array, so that
2497 * an element _x_ becomes \f$ numerator / x \f$.
2498 * \warning If an exception is thrown because of presence of 0.0 element in \a this
2499 * array, all elements processed before detection of the zero element remain
2501 * \param [in] numerator - the numerator used to modify array elements.
2502 * \throw If \a this is not allocated.
2503 * \throw If there is an element equal to 0.0 in \a this array.
2505 void DataArrayDouble::applyInv(double numerator)
2508 double *ptr=getPointer();
2509 std::size_t nbOfElems=getNbOfElems();
2510 for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2512 if(std::abs(*ptr)>std::numeric_limits<double>::min())
2514 *ptr=numerator/(*ptr);
2518 std::ostringstream oss; oss << "DataArrayDouble::applyInv : presence of null value in tuple #" << i/getNumberOfComponents() << " component #" << i%getNumberOfComponents();
2520 throw INTERP_KERNEL::Exception(oss.str().c_str());
2527 * Modify all elements of \a this array, so that
2528 * an element _x_ becomes <em> val ^ x </em>. Contrary to DataArrayInt::applyPow
2529 * all values in \a this have to be >= 0 if val is \b not integer.
2530 * \param [in] val - the value used to apply pow on all array elements.
2531 * \throw If \a this is not allocated.
2532 * \warning If an exception is thrown because of presence of 0 element in \a this
2533 * array and \a val is \b not integer, all elements processed before detection of the zero element remain
2536 void DataArrayDouble::applyPow(double val)
2539 double *ptr=getPointer();
2540 std::size_t nbOfElems=getNbOfElems();
2542 bool isInt=((double)val2)==val;
2545 for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2551 std::ostringstream oss; oss << "DataArrayDouble::applyPow (double) : At elem # " << i << " value is " << *ptr << " ! must be >=0. !";
2552 throw INTERP_KERNEL::Exception(oss.str().c_str());
2558 for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2559 *ptr=pow(*ptr,val2);
2565 * Modify all elements of \a this array, so that
2566 * an element _x_ becomes \f$ val ^ x \f$.
2567 * \param [in] val - the value used to apply pow on all array elements.
2568 * \throw If \a this is not allocated.
2569 * \throw If \a val < 0.
2570 * \warning If an exception is thrown because of presence of 0 element in \a this
2571 * array, all elements processed before detection of the zero element remain
2574 void DataArrayDouble::applyRPow(double val)
2578 throw INTERP_KERNEL::Exception("DataArrayDouble::applyRPow : the input value has to be >= 0 !");
2579 double *ptr=getPointer();
2580 std::size_t nbOfElems=getNbOfElems();
2581 for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2587 * Returns a new DataArrayDouble created from \a this one by applying \a
2588 * FunctionToEvaluate to every tuple of \a this array. Textual data is not copied.
2589 * For more info see \ref MEDCouplingArrayApplyFunc
2590 * \param [in] nbOfComp - number of components in the result array.
2591 * \param [in] func - the \a FunctionToEvaluate declared as
2592 * \c bool (*\a func)(\c const \c double *\a pos, \c double *\a res),
2593 * where \a pos points to the first component of a tuple of \a this array
2594 * and \a res points to the first component of a tuple of the result array.
2595 * Note that length (number of components) of \a pos can differ from
2597 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2598 * same number of tuples as \a this array.
2599 * The caller is to delete this result array using decrRef() as it is no more
2601 * \throw If \a this is not allocated.
2602 * \throw If \a func returns \a false.
2604 DataArrayDouble *DataArrayDouble::applyFunc(std::size_t nbOfComp, FunctionToEvaluate func) const
2607 DataArrayDouble *newArr=DataArrayDouble::New();
2608 mcIdType nbOfTuples(getNumberOfTuples());
2609 std::size_t oldNbOfComp(getNumberOfComponents());
2610 newArr->alloc(nbOfTuples,nbOfComp);
2611 const double *ptr=getConstPointer();
2612 double *ptrToFill=newArr->getPointer();
2613 for(mcIdType i=0;i<nbOfTuples;i++)
2615 if(!func(ptr+i*oldNbOfComp,ptrToFill+i*nbOfComp))
2617 std::ostringstream oss; oss << "For tuple # " << i << " with value (";
2618 std::copy(ptr+oldNbOfComp*i,ptr+oldNbOfComp*(i+1),std::ostream_iterator<double>(oss,", "));
2619 oss << ") : Evaluation of function failed !";
2621 throw INTERP_KERNEL::Exception(oss.str().c_str());
2628 * Returns a new DataArrayDouble created from \a this one by applying a function to every
2629 * tuple of \a this array. Textual data is not copied.
2630 * For more info see \ref MEDCouplingArrayApplyFunc1.
2631 * \param [in] nbOfComp - number of components in the result array.
2632 * \param [in] func - the expression defining how to transform a tuple of \a this array.
2633 * Supported expressions are described \ref MEDCouplingArrayApplyFuncExpr "here".
2634 * \param [in] isSafe - By default true. If true invalid operation (division by 0. acos of value > 1. ...) leads to a throw of an exception.
2635 * If false the computation is carried on without any notification. When false the evaluation is a little faster.
2636 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2637 * same number of tuples as \a this array and \a nbOfComp components.
2638 * The caller is to delete this result array using decrRef() as it is no more
2640 * \throw If \a this is not allocated.
2641 * \throw If computing \a func fails.
2643 DataArrayDouble *DataArrayDouble::applyFunc(std::size_t nbOfComp, const std::string& func, bool isSafe) const
2645 INTERP_KERNEL::ExprParser expr(func);
2647 std::set<std::string> vars;
2648 expr.getTrueSetOfVars(vars);
2649 std::vector<std::string> varsV(vars.begin(),vars.end());
2650 return applyFuncNamedCompo(nbOfComp,varsV,func,isSafe);
2654 * Returns a new DataArrayDouble created from \a this one by applying a function to every
2655 * tuple of \a this array. Textual data is not copied. This method works by tuples (whatever its size).
2656 * If \a this is a one component array, call applyFuncOnThis instead that performs the same work faster.
2658 * For more info see \ref MEDCouplingArrayApplyFunc0.
2659 * \param [in] func - the expression defining how to transform a tuple of \a this array.
2660 * Supported expressions are described \ref MEDCouplingArrayApplyFuncExpr "here".
2661 * \param [in] isSafe - By default true. If true invalid operation (division by 0. acos of value > 1. ...) leads to a throw of an exception.
2662 * If false the computation is carried on without any notification. When false the evaluation is a little faster.
2663 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2664 * same number of tuples and components as \a this array.
2665 * The caller is to delete this result array using decrRef() as it is no more
2667 * \sa applyFuncOnThis
2668 * \throw If \a this is not allocated.
2669 * \throw If computing \a func fails.
2671 DataArrayDouble *DataArrayDouble::applyFunc(const std::string& func, bool isSafe) const
2673 std::size_t nbOfComp(getNumberOfComponents());
2675 throw INTERP_KERNEL::Exception("DataArrayDouble::applyFunc : output number of component must be > 0 !");
2677 mcIdType nbOfTuples(getNumberOfTuples());
2678 MCAuto<DataArrayDouble> newArr(DataArrayDouble::New());
2679 newArr->alloc(nbOfTuples,nbOfComp);
2680 INTERP_KERNEL::ExprParser expr(func);
2682 std::set<std::string> vars;
2683 expr.getTrueSetOfVars(vars);
2686 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 : ";
2687 std::copy(vars.begin(),vars.end(),std::ostream_iterator<std::string>(oss," "));
2688 throw INTERP_KERNEL::Exception(oss.str().c_str());
2692 expr.prepareFastEvaluator();
2693 newArr->rearrange(1);
2694 newArr->fillWithValue(expr.evaluateDouble());
2695 newArr->rearrange(nbOfComp);
2696 return newArr.retn();
2698 std::vector<std::string> vars2(vars.begin(),vars.end());
2699 double buff,*ptrToFill(newArr->getPointer());
2700 const double *ptr(begin());
2701 std::vector<double> stck;
2702 expr.prepareExprEvaluationDouble(vars2,1,1,0,&buff,&buff+1);
2703 expr.prepareFastEvaluator();
2706 for(mcIdType i=0;i<nbOfTuples;i++)
2708 for(std::size_t iComp=0;iComp<nbOfComp;iComp++,ptr++,ptrToFill++)
2711 expr.evaluateDoubleInternal(stck);
2712 *ptrToFill=stck.back();
2719 for(mcIdType i=0;i<nbOfTuples;i++)
2721 for(std::size_t iComp=0;iComp<nbOfComp;iComp++,ptr++,ptrToFill++)
2726 expr.evaluateDoubleInternalSafe(stck);
2728 catch(INTERP_KERNEL::Exception& e)
2730 std::ostringstream oss; oss << "For tuple # " << i << " component # " << iComp << " with value (";
2732 oss << ") : Evaluation of function failed !" << e.what();
2733 throw INTERP_KERNEL::Exception(oss.str().c_str());
2735 *ptrToFill=stck.back();
2740 return newArr.retn();
2744 * This method is a non const method that modify the array in \a this.
2745 * This method only works on one component array. It means that function \a func must
2746 * contain at most one variable.
2747 * This method is a specialization of applyFunc method with one parameter on one component array.
2749 * \param [in] func - the expression defining how to transform a tuple of \a this array.
2750 * Supported expressions are described \ref MEDCouplingArrayApplyFuncExpr "here".
2751 * \param [in] isSafe - By default true. If true invalid operation (division by 0. acos of value > 1. ...) leads to a throw of an exception.
2752 * If false the computation is carried on without any notification. When false the evaluation is a little faster.
2756 void DataArrayDouble::applyFuncOnThis(const std::string& func, bool isSafe)
2758 std::size_t nbOfComp(getNumberOfComponents());
2760 throw INTERP_KERNEL::Exception("DataArrayDouble::applyFuncOnThis : output number of component must be > 0 !");
2762 mcIdType nbOfTuples(getNumberOfTuples());
2763 INTERP_KERNEL::ExprParser expr(func);
2765 std::set<std::string> vars;
2766 expr.getTrueSetOfVars(vars);
2769 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 : ";
2770 std::copy(vars.begin(),vars.end(),std::ostream_iterator<std::string>(oss," "));
2771 throw INTERP_KERNEL::Exception(oss.str().c_str());
2775 expr.prepareFastEvaluator();
2776 std::vector<std::string> compInfo(getInfoOnComponents());
2778 fillWithValue(expr.evaluateDouble());
2779 rearrange(nbOfComp);
2780 setInfoOnComponents(compInfo);
2783 std::vector<std::string> vars2(vars.begin(),vars.end());
2784 double buff,*ptrToFill(getPointer());
2785 const double *ptr(begin());
2786 std::vector<double> stck;
2787 expr.prepareExprEvaluationDouble(vars2,1,1,0,&buff,&buff+1);
2788 expr.prepareFastEvaluator();
2791 for(mcIdType i=0;i<nbOfTuples;i++)
2793 for(std::size_t iComp=0;iComp<nbOfComp;iComp++,ptr++,ptrToFill++)
2796 expr.evaluateDoubleInternal(stck);
2797 *ptrToFill=stck.back();
2804 for(mcIdType i=0;i<nbOfTuples;i++)
2806 for(std::size_t iComp=0;iComp<nbOfComp;iComp++,ptr++,ptrToFill++)
2811 expr.evaluateDoubleInternalSafe(stck);
2813 catch(INTERP_KERNEL::Exception& e)
2815 std::ostringstream oss; oss << "For tuple # " << i << " component # " << iComp << " with value (";
2817 oss << ") : Evaluation of function failed !" << e.what();
2818 throw INTERP_KERNEL::Exception(oss.str().c_str());
2820 *ptrToFill=stck.back();
2828 * Returns a new DataArrayDouble created from \a this one by applying a function to every
2829 * tuple of \a this array. Textual data is not copied.
2830 * For more info see \ref MEDCouplingArrayApplyFunc2.
2831 * \param [in] nbOfComp - number of components in the result array.
2832 * \param [in] func - the expression defining how to transform a tuple of \a this array.
2833 * Supported expressions are described \ref MEDCouplingArrayApplyFuncExpr "here".
2834 * \param [in] isSafe - By default true. If true invalid operation (division by 0. acos of value > 1. ...) leads to a throw of an exception.
2835 * If false the computation is carried on without any notification. When false the evaluation is a little faster.
2836 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2837 * same number of tuples as \a this array.
2838 * The caller is to delete this result array using decrRef() as it is no more
2840 * \throw If \a this is not allocated.
2841 * \throw If \a func contains vars that are not in \a this->getInfoOnComponent().
2842 * \throw If computing \a func fails.
2844 DataArrayDouble *DataArrayDouble::applyFuncCompo(std::size_t nbOfComp, const std::string& func, bool isSafe) const
2846 return applyFuncNamedCompo(nbOfComp,getVarsOnComponent(),func,isSafe);
2850 * Returns a new DataArrayDouble created from \a this one by applying a function to every
2851 * tuple of \a this array. Textual data is not copied.
2852 * For more info see \ref MEDCouplingArrayApplyFunc3.
2853 * \param [in] nbOfComp - number of components in the result array.
2854 * \param [in] varsOrder - sequence of vars defining their order.
2855 * \param [in] func - the expression defining how to transform a tuple of \a this array.
2856 * Supported expressions are described \ref MEDCouplingArrayApplyFuncExpr "here".
2857 * \param [in] isSafe - By default true. If true invalid operation (division by 0. acos of value > 1. ...) leads to a throw of an exception.
2858 * If false the computation is carried on without any notification. When false the evaluation is a little faster.
2859 * \return DataArrayDouble * - the new instance of DataArrayDouble containing the
2860 * same number of tuples as \a this array.
2861 * The caller is to delete this result array using decrRef() as it is no more
2863 * \throw If \a this is not allocated.
2864 * \throw If \a func contains vars not in \a varsOrder.
2865 * \throw If computing \a func fails.
2867 DataArrayDouble *DataArrayDouble::applyFuncNamedCompo(std::size_t nbOfComp, const std::vector<std::string>& varsOrder, const std::string& func, bool isSafe) const
2870 throw INTERP_KERNEL::Exception("DataArrayDouble::applyFuncNamedCompo : output number of component must be > 0 !");
2871 std::vector<std::string> varsOrder2(varsOrder);
2872 std::size_t oldNbOfComp(getNumberOfComponents());
2873 for(std::size_t i=varsOrder.size();i<oldNbOfComp;i++)
2874 varsOrder2.push_back(std::string());
2876 mcIdType nbOfTuples(getNumberOfTuples());
2877 INTERP_KERNEL::ExprParser expr(func);
2879 std::set<std::string> vars;
2880 expr.getTrueSetOfVars(vars);
2881 if(vars.size()>oldNbOfComp)
2883 std::ostringstream oss; oss << "The field has " << oldNbOfComp << " components and there are ";
2884 oss << vars.size() << " variables : ";
2885 std::copy(vars.begin(),vars.end(),std::ostream_iterator<std::string>(oss," "));
2886 throw INTERP_KERNEL::Exception(oss.str().c_str());
2888 MCAuto<DataArrayDouble> newArr(DataArrayDouble::New());
2889 newArr->alloc(nbOfTuples,nbOfComp);
2890 INTERP_KERNEL::AutoPtr<double> buff(new double[oldNbOfComp]);
2891 double *buffPtr(buff),*ptrToFill;
2892 std::vector<double> stck;
2893 for(std::size_t iComp=0;iComp<nbOfComp;iComp++)
2895 expr.prepareExprEvaluationDouble(varsOrder2,(int)oldNbOfComp,(int)nbOfComp,(int)iComp,buffPtr,buffPtr+oldNbOfComp);
2896 expr.prepareFastEvaluator();
2897 const double *ptr(getConstPointer());
2898 ptrToFill=newArr->getPointer()+iComp;
2901 for(mcIdType i=0;i<nbOfTuples;i++,ptrToFill+=nbOfComp,ptr+=oldNbOfComp)
2903 std::copy(ptr,ptr+oldNbOfComp,buffPtr);
2904 expr.evaluateDoubleInternal(stck);
2905 *ptrToFill=stck.back();
2911 for(mcIdType i=0;i<nbOfTuples;i++,ptrToFill+=nbOfComp,ptr+=oldNbOfComp)
2913 std::copy(ptr,ptr+oldNbOfComp,buffPtr);
2916 expr.evaluateDoubleInternalSafe(stck);
2917 *ptrToFill=stck.back();
2920 catch(INTERP_KERNEL::Exception& e)
2922 std::ostringstream oss; oss << "For tuple # " << i << " with value (";
2923 std::copy(ptr+oldNbOfComp*i,ptr+oldNbOfComp*(i+1),std::ostream_iterator<double>(oss,", "));
2924 oss << ") : Evaluation of function failed !" << e.what();
2925 throw INTERP_KERNEL::Exception(oss.str().c_str());
2930 return newArr.retn();
2933 void DataArrayDouble::applyFuncFast32(const std::string& func)
2936 INTERP_KERNEL::ExprParser expr(func);
2938 char *funcStr=expr.compileX86();
2940 *((void **)&funcPtr)=funcStr;//he he...
2942 double *ptr=getPointer();
2943 std::size_t nbOfComp=getNumberOfComponents();
2944 mcIdType nbOfTuples=getNumberOfTuples();
2945 std::size_t nbOfElems=nbOfTuples*nbOfComp;
2946 for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2951 void DataArrayDouble::applyFuncFast64(const std::string& func)
2954 INTERP_KERNEL::ExprParser expr(func);
2956 char *funcStr=expr.compileX86_64();
2958 *((void **)&funcPtr)=funcStr;//he he...
2960 double *ptr=getPointer();
2961 std::size_t nbOfComp=getNumberOfComponents();
2962 mcIdType nbOfTuples=getNumberOfTuples();
2963 std::size_t nbOfElems=nbOfTuples*nbOfComp;
2964 for(std::size_t i=0;i<nbOfElems;i++,ptr++)
2970 * \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.
2972 MCAuto<DataArrayDouble> DataArrayDouble::symmetry3DPlane(const double point[3], const double normalVector[3]) const
2975 if(getNumberOfComponents()!=3)
2976 throw INTERP_KERNEL::Exception("DataArrayDouble::symmetry3DPlane : this is excepted to have 3 components !");
2977 mcIdType nbTuples(getNumberOfTuples());
2978 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
2979 ret->alloc(nbTuples,3);
2980 Symmetry3DPlane(point,normalVector,nbTuples,begin(),ret->getPointer());
2984 DataArrayDoubleIterator *DataArrayDouble::iterator()
2986 return new DataArrayDoubleIterator(this);
2990 * Returns a new DataArrayInt containing indices of tuples of \a this one-dimensional
2991 * array whose values are within a given range. Textual data is not copied.
2992 * \param [in] vmin - a lowest acceptable value (included).
2993 * \param [in] vmax - a greatest acceptable value (included).
2994 * \return DataArrayInt * - the new instance of DataArrayInt.
2995 * The caller is to delete this result array using decrRef() as it is no more
2997 * \throw If \a this->getNumberOfComponents() != 1.
2999 * \sa DataArrayDouble::findIdsNotInRange
3001 * \if ENABLE_EXAMPLES
3002 * \ref cpp_mcdataarraydouble_getidsinrange "Here is a C++ example".<br>
3003 * \ref py_mcdataarraydouble_getidsinrange "Here is a Python example".
3006 DataArrayIdType *DataArrayDouble::findIdsInRange(double vmin, double vmax) const
3009 if(getNumberOfComponents()!=1)
3010 throw INTERP_KERNEL::Exception("DataArrayDouble::findIdsInRange : this must have exactly one component !");
3011 const double *cptr(begin());
3012 MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(0,1);
3013 mcIdType nbOfTuples(getNumberOfTuples());
3014 for(mcIdType i=0;i<nbOfTuples;i++,cptr++)
3015 if(*cptr>=vmin && *cptr<=vmax)
3016 ret->pushBackSilent(i);
3021 * Returns a new DataArrayInt containing indices of tuples of \a this one-dimensional
3022 * array whose values are not within a given range. Textual data is not copied.
3023 * \param [in] vmin - a lowest not acceptable value (excluded).
3024 * \param [in] vmax - a greatest not acceptable value (excluded).
3025 * \return DataArrayInt * - the new instance of DataArrayInt.
3026 * The caller is to delete this result array using decrRef() as it is no more
3028 * \throw If \a this->getNumberOfComponents() != 1.
3030 * \sa DataArrayDouble::findIdsInRange
3032 DataArrayIdType *DataArrayDouble::findIdsNotInRange(double vmin, double vmax) const
3035 if(getNumberOfComponents()!=1)
3036 throw INTERP_KERNEL::Exception("DataArrayDouble::findIdsNotInRange : this must have exactly one component !");
3037 const double *cptr(begin());
3038 MCAuto<DataArrayIdType> ret(DataArrayIdType::New()); ret->alloc(0,1);
3039 mcIdType nbOfTuples(getNumberOfTuples());
3040 for(mcIdType i=0;i<nbOfTuples;i++,cptr++)
3041 if(*cptr<vmin || *cptr>vmax)
3042 ret->pushBackSilent(i);
3047 * Returns a new DataArrayDouble by concatenating two given arrays, so that (1) the number
3048 * of tuples in the result array is a sum of the number of tuples of given arrays and (2)
3049 * the number of component in the result array is same as that of each of given arrays.
3050 * Info on components is copied from the first of the given arrays. Number of components
3051 * in the given arrays must be the same.
3052 * \param [in] a1 - an array to include in the result array.
3053 * \param [in] a2 - another array to include in the result array.
3054 * \return DataArrayDouble * - the new instance of DataArrayDouble.
3055 * The caller is to delete this result array using decrRef() as it is no more
3057 * \throw If both \a a1 and \a a2 are NULL.
3058 * \throw If \a a1->getNumberOfComponents() != \a a2->getNumberOfComponents().
3060 DataArrayDouble *DataArrayDouble::Aggregate(const DataArrayDouble *a1, const DataArrayDouble *a2)
3062 std::vector<const DataArrayDouble *> tmp(2);
3063 tmp[0]=a1; tmp[1]=a2;
3064 return Aggregate(tmp);
3068 * Returns a new DataArrayDouble by concatenating all given arrays, so that (1) the number
3069 * of tuples in the result array is a sum of the number of tuples of given arrays and (2)
3070 * the number of component in the result array is same as that of each of given arrays.
3071 * Info on components is copied from the first of the given arrays. Number of components
3072 * in the given arrays must be the same.
3073 * If the number of non null of elements in \a arr is equal to one the returned object is a copy of it
3074 * not the object itself.
3075 * \param [in] arr - a sequence of arrays to include in the result array.
3076 * \return DataArrayDouble * - the new instance of DataArrayDouble.
3077 * The caller is to delete this result array using decrRef() as it is no more
3079 * \throw If all arrays within \a arr are NULL.
3080 * \throw If getNumberOfComponents() of arrays within \a arr.
3082 DataArrayDouble *DataArrayDouble::Aggregate(const std::vector<const DataArrayDouble *>& arr)
3084 std::vector<const DataArrayDouble *> a;
3085 for(std::vector<const DataArrayDouble *>::const_iterator it4=arr.begin();it4!=arr.end();it4++)
3089 throw INTERP_KERNEL::Exception("DataArrayDouble::Aggregate : input list must contain at least one NON EMPTY DataArrayDouble !");
3090 std::vector<const DataArrayDouble *>::const_iterator it=a.begin();
3091 std::size_t nbOfComp((*it)->getNumberOfComponents());
3092 mcIdType nbt=(*it++)->getNumberOfTuples();
3093 for(mcIdType i=1;it!=a.end();it++,i++)
3095 if((*it)->getNumberOfComponents()!=nbOfComp)
3096 throw INTERP_KERNEL::Exception("DataArrayDouble::Aggregate : Nb of components mismatch for array aggregation !");
3097 nbt+=(*it)->getNumberOfTuples();
3099 MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
3100 ret->alloc(nbt,nbOfComp);
3101 double *pt=ret->getPointer();
3102 for(it=a.begin();it!=a.end();it++)
3103 pt=std::copy((*it)->getConstPointer(),(*it)->getConstPointer()+(*it)->getNbOfElems(),pt);
3104 ret->copyStringInfoFrom(*(a[0]));
3109 * Returns a new DataArrayDouble containing a dot product of two given arrays, so that
3110 * the i-th tuple of the result array is a sum of products of j-th components of i-th
3111 * tuples of given arrays (\f$ a_i = \sum_{j=1}^n a1_j * a2_j \f$).
3112 * Info on components and name is copied from the first of the given arrays.
3113 * Number of tuples and components in the given arrays must be the same.
3114 * \param [in] a1 - a given array.
3115 * \param [in] a2 - another given array.
3116 * \return DataArrayDouble * - the new instance of DataArrayDouble.
3117 * The caller is to delete this result array using decrRef() as it is no more
3119 * \throw If either \a a1 or \a a2 is NULL.
3120 * \throw If any given array is not allocated.
3121 * \throw If \a a1->getNumberOfTuples() != \a a2->getNumberOfTuples()
3122 * \throw If \a a1->getNumberOfComponents() != \a a2->getNumberOfComponents()
3124 DataArrayDouble *DataArrayDouble::Dot(const DataArrayDouble *a1, const DataArrayDouble *a2)
3127 throw INTERP_KERNEL::Exception("DataArrayDouble::Dot : input DataArrayDouble instance is NULL !");
3128 a1->checkAllocated();
3129 a2->checkAllocated();
3130 std::size_t nbOfComp(a1->getNumberOfComponents());
3131 if(nbOfComp!=a2->getNumberOfComponents())
3132 throw INTERP_KERNEL::Exception("Nb of components mismatch for array Dot !");
3133 mcIdType nbOfTuple(a1->getNumberOfTuples());
3134 if(nbOfTuple!=a2->getNumberOfTuples())
3135 throw INTERP_KERNEL::Exception("Nb of tuples mismatch for array Dot !");
3136 DataArrayDouble *ret=DataArrayDouble::New();
3137 ret->alloc(nbOfTuple,1);
3138 double *retPtr=ret->getPointer();
3139 const double *a1Ptr=a1->begin(),*a2Ptr(a2->begin());
3140 for(mcIdType i=0;i<nbOfTuple;i++)
3143 for(std::size_t j=0;j<nbOfComp;j++)
3144 sum+=a1Ptr[i*nbOfComp+j]*a2Ptr[i*nbOfComp+j];
3147 ret->setInfoOnComponent(0,a1->getInfoOnComponent(0));
3148 ret->setName(a1->getName());
3153 * Returns a new DataArrayDouble containing a cross product of two given arrays, so that
3154 * the i-th tuple of the result array contains 3 components of a vector which is a cross
3155 * product of two vectors defined by the i-th tuples of given arrays.
3156 * Info on components is copied from the first of the given arrays.
3157 * Number of tuples in the given arrays must be the same.
3158 * Number of components in the given arrays must be 3.
3159 * \param [in] a1 - a given array.
3160 * \param [in] a2 - another given array.
3161 * \return DataArrayDouble * - the new instance of DataArrayDouble.
3162 * The caller is to delete this result array using decrRef() as it is no more
3164 * \throw If either \a a1 or \a a2 is NULL.
3165 * \throw If \a a1->getNumberOfTuples() != \a a2->getNumberOfTuples()
3166 * \throw If \a a1->getNumberOfComponents() != 3
3167 * \throw If \a a2->getNumberOfComponents() != 3
3169 DataArrayDouble *DataArrayDouble::CrossProduct(const DataArrayDouble *a1, const DataArrayDouble *a2)
3172 throw INTERP_KERNEL::Exception("DataArrayDouble::CrossProduct : input DataArrayDouble instance is NULL !");
3173 std::size_t nbOfComp(a1->getNumberOfComponents());
3174 if(nbOfComp!=a2->getNumberOfComponents())
3175 throw INTERP_KERNEL::Exception("Nb of components mismatch for array crossProduct !");
3177 throw INTERP_KERNEL::Exception("Nb of components must be equal to 3 for array crossProduct !");
3178 mcIdType nbOfTuple(a1->getNumberOfTuples());
3179 if(nbOfTuple!=a2->getNumberOfTuples())
3180 throw INTERP_KERNEL::Exception("Nb of tuples mismatch for array crossProduct !");
3181 DataArrayDouble *ret=DataArrayDouble::New();
3182 ret->alloc(nbOfTuple,3);
3183 double *retPtr=ret->getPointer();
3184 const double *a1Ptr(a1->begin()),*a2Ptr(a2->begin());
3185 for(mcIdType i=0;i<nbOfTuple;i++)
3187 retPtr[3*i]=a1Ptr[3*i+1]*a2Ptr[3*i+2]-a1Ptr[3*i+2]*a2Ptr[3*i+1];
3188 retPtr[3*i+1]=a1Ptr[3*i+2]*a2Ptr[3*i]-a1Ptr[3*i]*a2Ptr[3*i+2];
3189 retPtr[3*i+2]=a1Ptr[3*i]*a2Ptr[3*i+1]-a1Ptr[3*i+1]*a2Ptr[3*i];
3191 ret->copyStringInfoFrom(*a1);
3196 * Returns a new DataArrayDouble containing maximal values of two given arrays.
3197 * Info on components is copied from the first of the given arrays.
3198 * Number of tuples and components in the given arrays must be the same.
3199 * \param [in] a1 - an array to compare values with another one.
3200 * \param [in] a2 - another array to compare values with the first one.
3201 * \return DataArrayDouble * - the new instance of DataArrayDouble.
3202 * The caller is to delete this result array using decrRef() as it is no more
3204 * \throw If either \a a1 or \a a2 is NULL.
3205 * \throw If \a a1->getNumberOfTuples() != \a a2->getNumberOfTuples()
3206 * \throw If \a a1->getNumberOfComponents() != \a a2->getNumberOfComponents()
3208 DataArrayDouble *DataArrayDouble::Max(const DataArrayDouble *a1, const DataArrayDouble *a2)
3211 throw INTERP_KERNEL::Exception("DataArrayDouble::Max : input DataArrayDouble instance is NULL !");
3212 std::size_t nbOfComp(a1->getNumberOfComponents());
3213 if(nbOfComp!=a2->getNumberOfComponents())
3214 throw INTERP_KERNEL::Exception("Nb of components mismatch for array Max !");
3215 mcIdType nbOfTuple(a1->getNumberOfTuples());
3216 if(nbOfTuple!=a2->getNumberOfTuples())
3217 throw INTERP_KERNEL::Exception("Nb of tuples mismatch for array Max !");
3218 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
3219 ret->alloc(nbOfTuple,nbOfComp);
3220 double *retPtr(ret->getPointer());
3221 const double *a1Ptr(a1->begin()),*a2Ptr(a2->begin());
3222 std::size_t nbElem(nbOfTuple*nbOfComp);
3223 for(std::size_t i=0;i<nbElem;i++)
3224 retPtr[i]=std::max(a1Ptr[i],a2Ptr[i]);
3225 ret->copyStringInfoFrom(*a1);
3230 * Returns a new DataArrayDouble containing minimal values of two given arrays.
3231 * Info on components is copied from the first of the given arrays.
3232 * Number of tuples and components in the given arrays must be the same.
3233 * \param [in] a1 - an array to compare values with another one.
3234 * \param [in] a2 - another array to compare values with the first one.
3235 * \return DataArrayDouble * - the new instance of DataArrayDouble.
3236 * The caller is to delete this result array using decrRef() as it is no more
3238 * \throw If either \a a1 or \a a2 is NULL.
3239 * \throw If \a a1->getNumberOfTuples() != \a a2->getNumberOfTuples()
3240 * \throw If \a a1->getNumberOfComponents() != \a a2->getNumberOfComponents()
3242 DataArrayDouble *DataArrayDouble::Min(const DataArrayDouble *a1, const DataArrayDouble *a2)
3245 throw INTERP_KERNEL::Exception("DataArrayDouble::Min : input DataArrayDouble instance is NULL !");
3246 std::size_t nbOfComp(a1->getNumberOfComponents());
3247 if(nbOfComp!=a2->getNumberOfComponents())
3248 throw INTERP_KERNEL::Exception("Nb of components mismatch for array min !");
3249 mcIdType nbOfTuple(a1->getNumberOfTuples());
3250 if(nbOfTuple!=a2->getNumberOfTuples())
3251 throw INTERP_KERNEL::Exception("Nb of tuples mismatch for array min !");
3252 MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
3253 ret->alloc(nbOfTuple,nbOfComp);
3254 double *retPtr(ret->getPointer());
3255 const double *a1Ptr(a1->begin()),*a2Ptr(a2->begin());
3256 std::size_t nbElem(nbOfTuple*nbOfComp);
3257 for(std::size_t i=0;i<nbElem;i++)
3258 retPtr[i]=std::min(a1Ptr[i],a2Ptr[i]);
3259 ret->copyStringInfoFrom(*a1);
3264 * Returns a new DataArrayDouble that is the result of pow of two given arrays. There are 3
3267 * \param [in] a1 - an array to pow up.
3268 * \param [in] a2 - another array to sum up.
3269 * \return DataArrayDouble * - the new instance of DataArrayDouble.
3270 * The caller is to delete this result array using decrRef() as it is no more
3272 * \throw If either \a a1 or \a a2 is NULL.
3273 * \throw If \a a1->getNumberOfTuples() != \a a2->getNumberOfTuples()
3274 * \throw If \a a1->getNumberOfComponents() != 1 or \a a2->getNumberOfComponents() != 1.
3275 * \throw If there is a negative value in \a a1.
3277 DataArrayDouble *DataArrayDouble::Pow(const DataArrayDouble *a1, const DataArrayDouble *a2)
3280 throw INTERP_KERNEL::Exception("DataArrayDouble::Pow : at least one of input instances is null !");
3281 mcIdType nbOfTuple=a1->getNumberOfTuples();
3282 mcIdType nbOfTuple2=a2->getNumberOfTuples();
3283 std::size_t nbOfComp=a1->getNumberOfComponents();
3284 std::size_t nbOfComp2=a2->getNumberOfComponents();
3285 if(nbOfTuple!=nbOfTuple2)
3286 throw INTERP_KERNEL::Exception("DataArrayDouble::Pow : number of tuples mismatches !");
3287 if(nbOfComp!=1 || nbOfComp2!=1)
3288 throw INTERP_KERNEL::Exception("DataArrayDouble::Pow : number of components of both arrays must be equal to 1 !");
3289 MCAuto<DataArrayDouble> ret=DataArrayDouble::New(); ret->alloc(nbOfTuple,1);
3290 const double *ptr1(a1->begin()),*ptr2(a2->begin());
3291 double *ptr=ret->getPointer();
3292 for(mcIdType i=0;i<nbOfTuple;i++,ptr1++,ptr2++,ptr++)
3296 *ptr=pow(*ptr1,*ptr2);
3300 std::ostringstream oss; oss << "DataArrayDouble::Pow : on tuple #" << i << " of a1 value is < 0 (" << *ptr1 << ") !";
3301 throw INTERP_KERNEL::Exception(oss.str().c_str());
3308 * Apply pow on values of another DataArrayDouble to values of \a this one.
3310 * \param [in] other - an array to pow to \a this one.
3311 * \throw If \a other is NULL.
3312 * \throw If \a this->getNumberOfTuples() != \a other->getNumberOfTuples()
3313 * \throw If \a this->getNumberOfComponents() != 1 or \a other->getNumberOfComponents() != 1
3314 * \throw If there is a negative value in \a this.
3316 void DataArrayDouble::powEqual(const DataArrayDouble *other)
3319 throw INTERP_KERNEL::Exception("DataArrayDouble::powEqual : input instance is null !");
3320 mcIdType nbOfTuple=getNumberOfTuples();
3321 mcIdType nbOfTuple2=other->getNumberOfTuples();
3322 std::size_t nbOfComp=getNumberOfComponents();
3323 std::size_t nbOfComp2=other->getNumberOfComponents();
3324 if(nbOfTuple!=nbOfTuple2)
3325 throw INTERP_KERNEL::Exception("DataArrayDouble::powEqual : number of tuples mismatches !");
3326 if(nbOfComp!=1 || nbOfComp2!=1)
3327 throw INTERP_KERNEL::Exception("DataArrayDouble::powEqual : number of components of both arrays must be equal to 1 !");
3328 double *ptr=getPointer();
3329 const double *ptrc=other->begin();
3330 for(mcIdType i=0;i<nbOfTuple;i++,ptrc++,ptr++)
3333 *ptr=pow(*ptr,*ptrc);
3336 std::ostringstream oss; oss << "DataArrayDouble::powEqual : on tuple #" << i << " of this value is < 0 (" << *ptr << ") !";
3337 throw INTERP_KERNEL::Exception(oss.str().c_str());
3344 * This method is \b NOT wrapped into python because it can be useful only for performance reasons in C++ context.
3345 * All values in \a this must be 0. or 1. within eps error. 0 means false, 1 means true.
3346 * If an another value than 0 or 1 appear (within eps precision) an INTERP_KERNEL::Exception will be thrown.
3348 * \throw if \a this is not allocated.
3349 * \throw if \a this has not exactly one component.
3351 std::vector<bool> DataArrayDouble::toVectorOfBool(double eps) const
3354 if(getNumberOfComponents()!=1)
3355 throw INTERP_KERNEL::Exception("DataArrayDouble::toVectorOfBool : must be applied on single component array !");
3356 mcIdType nbt(getNumberOfTuples());
3357 std::vector<bool> ret(nbt);
3358 const double *pt(begin());
3359 for(mcIdType i=0;i<nbt;i++)
3363 else if(fabs(pt[i]-1.)<eps)
3367 std::ostringstream oss; oss << "DataArrayDouble::toVectorOfBool : the tuple #" << i << " has value " << pt[i] << " is invalid ! must be 0. or 1. !";
3368 throw INTERP_KERNEL::Exception(oss.str().c_str());
3375 * Useless method for end user. Only for MPI/Corba/File serialsation for multi arrays class.
3378 void DataArrayDouble::getTinySerializationIntInformation(std::vector<mcIdType>& tinyInfo) const
3383 tinyInfo[0]=getNumberOfTuples();
3384 tinyInfo[1]=ToIdType(getNumberOfComponents());
3394 * Useless method for end user. Only for MPI/Corba/File serialsation for multi arrays class.
3397 void DataArrayDouble::getTinySerializationStrInformation(std::vector<std::string>& tinyInfo) const
3401 std::size_t nbOfCompo(getNumberOfComponents());
3402 tinyInfo.resize(nbOfCompo+1);
3403 tinyInfo[0]=getName();
3404 for(std::size_t i=0;i<nbOfCompo;i++)
3405 tinyInfo[i+1]=getInfoOnComponent(i);
3410 tinyInfo[0]=getName();
3415 * Useless method for end user. Only for MPI/Corba/File serialsation for multi arrays class.
3416 * This method returns if a feeding is needed.
3418 bool DataArrayDouble::resizeForUnserialization(const std::vector<mcIdType>& tinyInfoI)
3420 mcIdType nbOfTuple=tinyInfoI[0];
3421 mcIdType nbOfComp=tinyInfoI[1];
3422 if(nbOfTuple!=-1 || nbOfComp!=-1)
3424 alloc(nbOfTuple,nbOfComp);
3431 * Useless method for end user. Only for MPI/Corba/File serialsation for multi arrays class.
3433 void DataArrayDouble::finishUnserialization(const std::vector<mcIdType>& tinyInfoI, const std::vector<std::string>& tinyInfoS)
3435 setName(tinyInfoS[0]);
3438 std::size_t nbOfCompo(getNumberOfComponents());
3439 for(std::size_t i=0;i<nbOfCompo;i++)
3440 setInfoOnComponent(i,tinyInfoS[i+1]);
3445 * Low static method that operates 3D rotation of 'nbNodes' 3D nodes whose coordinates are arranged in \a coordsIn
3446 * around an axe ( \a center, \a vect) and with angle \a angle.
3448 void DataArrayDouble::Rotate3DAlg(const double *center, const double *vect, double angle, mcIdType nbNodes, const double *coordsIn, double *coordsOut)
3450 if(!center || !vect)
3451 throw INTERP_KERNEL::Exception("DataArrayDouble::Rotate3DAlg : null vector in input !");
3452 double sina(sin(angle));
3453 double cosa(cos(angle));
3454 double vectorNorm[3];
3456 double matrixTmp[9];
3457 double norm(sqrt(vect[0]*vect[0]+vect[1]*vect[1]+vect[2]*vect[2]));
3458 if(norm<std::numeric_limits<double>::min())
3459 throw INTERP_KERNEL::Exception("DataArrayDouble::Rotate3DAlg : magnitude of input vector is too close of 0. !");
3460 std::transform(vect,vect+3,vectorNorm,std::bind2nd(std::multiplies<double>(),1/norm));
3461 //rotation matrix computation
3462 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;
3463 matrixTmp[0]=vectorNorm[0]*vectorNorm[0]; matrixTmp[1]=vectorNorm[0]*vectorNorm[1]; matrixTmp[2]=vectorNorm[0]*vectorNorm[2];
3464 matrixTmp[3]=vectorNorm[1]*vectorNorm[0]; matrixTmp[4]=vectorNorm[1]*vectorNorm[1]; matrixTmp[5]=vectorNorm[1]*vectorNorm[2];
3465 matrixTmp[6]=vectorNorm[2]*vectorNorm[0]; matrixTmp[7]=vectorNorm[2]*vectorNorm[1]; matrixTmp[8]=vectorNorm[2]*vectorNorm[2];
3466 std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies<double>(),1-cosa));
3467 std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
3468 matrixTmp[0]=0.; matrixTmp[1]=-vectorNorm[2]; matrixTmp[2]=vectorNorm[1];
3469 matrixTmp[3]=vectorNorm[2]; matrixTmp[4]=0.; matrixTmp[5]=-vectorNorm[0];
3470 matrixTmp[6]=-vectorNorm[1]; matrixTmp[7]=vectorNorm[0]; matrixTmp[8]=0.;
3471 std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies<double>(),sina));
3472 std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
3473 //rotation matrix computed.
3475 for(mcIdType i=0; i<nbNodes; i++)
3477 std::transform(coordsIn+i*3,coordsIn+(i+1)*3,center,tmp,std::minus<double>());
3478 coordsOut[i*3]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+matrix[2]*tmp[2]+center[0];
3479 coordsOut[i*3+1]=matrix[3]*tmp[0]+matrix[4]*tmp[1]+matrix[5]*tmp[2]+center[1];
3480 coordsOut[i*3+2]=matrix[6]*tmp[0]+matrix[7]*tmp[1]+matrix[8]*tmp[2]+center[2];
3484 void DataArrayDouble::Symmetry3DPlane(const double point[3], const double normalVector[3], mcIdType nbNodes, const double *coordsIn, double *coordsOut)
3486 double matrix[9],matrix2[9],matrix3[9];
3487 double vect[3],crossVect[3];
3488 INTERP_KERNEL::orthogonalVect3(normalVector,vect);
3489 crossVect[0]=normalVector[1]*vect[2]-normalVector[2]*vect[1];
3490 crossVect[1]=normalVector[2]*vect[0]-normalVector[0]*vect[2];
3491 crossVect[2]=normalVector[0]*vect[1]-normalVector[1]*vect[0];
3492 double nv(INTERP_KERNEL::norm<3>(vect)),ni(INTERP_KERNEL::norm<3>(normalVector)),nc(INTERP_KERNEL::norm<3>(crossVect));
3493 matrix[0]=vect[0]/nv; matrix[1]=crossVect[0]/nc; matrix[2]=-normalVector[0]/ni;
3494 matrix[3]=vect[1]/nv; matrix[4]=crossVect[1]/nc; matrix[5]=-normalVector[1]/ni;
3495 matrix[6]=vect[2]/nv; matrix[7]=crossVect[2]/nc; matrix[8]=-normalVector[2]/ni;
3496 matrix2[0]=vect[0]/nv; matrix2[1]=vect[1]/nv; matrix2[2]=vect[2]/nv;
3497 matrix2[3]=crossVect[0]/nc; matrix2[4]=crossVect[1]/nc; matrix2[5]=crossVect[2]/nc;
3498 matrix2[6]=normalVector[0]/ni; matrix2[7]=normalVector[1]/ni; matrix2[8]=normalVector[2]/ni;
3499 for(mcIdType i=0;i<3;i++)
3500 for(mcIdType j=0;j<3;j++)
3503 for(mcIdType k=0;k<3;k++)
3504 val+=matrix[3*i+k]*matrix2[3*k+j];
3507 //rotation matrix computed.
3509 for(mcIdType i=0; i<nbNodes; i++)
3511 std::transform(coordsIn+i*3,coordsIn+(i+1)*3,point,tmp,std::minus<double>());
3512 coordsOut[i*3]=matrix3[0]*tmp[0]+matrix3[1]*tmp[1]+matrix3[2]*tmp[2]+point[0];
3513 coordsOut[i*3+1]=matrix3[3]*tmp[0]+matrix3[4]*tmp[1]+matrix3[5]*tmp[2]+point[1];
3514 coordsOut[i*3+2]=matrix3[6]*tmp[0]+matrix3[7]*tmp[1]+matrix3[8]*tmp[2]+point[2];
3518 void DataArrayDouble::GiveBaseForPlane(const double normalVector[3], double baseOfPlane[9])
3520 double vect[3],crossVect[3];
3521 INTERP_KERNEL::orthogonalVect3(normalVector,vect);
3522 crossVect[0]=normalVector[1]*vect[2]-normalVector[2]*vect[1];
3523 crossVect[1]=normalVector[2]*vect[0]-normalVector[0]*vect[2];
3524 crossVect[2]=normalVector[0]*vect[1]-normalVector[1]*vect[0];
3525 double nv(INTERP_KERNEL::norm<3>(vect)),ni(INTERP_KERNEL::norm<3>(normalVector)),nc(INTERP_KERNEL::norm<3>(crossVect));
3526 baseOfPlane[0]=vect[0]/nv; baseOfPlane[1]=vect[1]/nv; baseOfPlane[2]=vect[2]/nv;
3527 baseOfPlane[3]=crossVect[0]/nc; baseOfPlane[4]=crossVect[1]/nc; baseOfPlane[5]=crossVect[2]/nc;
3528 baseOfPlane[6]=normalVector[0]/ni; baseOfPlane[7]=normalVector[1]/ni; baseOfPlane[8]=normalVector[2]/ni;
3532 * \param [in] seg2 : coordinates of input seg2 expected to have spacedim==2
3533 * \param [in] tri3 : coordinates of input tri3 also expected to have spacedim==2
3534 * \param [out] coeffs : the result of integration normalized to 1. along \a seg2 inside tri3 sorted by the node id of \a tri3
3535 * \param [out] length : the length of seg2. That is too say the length of integration
3537 void DataArrayDouble::ComputeIntegralOfSeg2IntoTri3(const double seg2[4], const double tri3[6], double coeffs[3], double& length)
3539 length=INTERP_KERNEL::norme_vecteur(seg2,seg2+2);
3541 INTERP_KERNEL::mid_of_seg2(seg2,seg2+2,mid);
3542 INTERP_KERNEL::barycentric_coords<2>(tri3,mid,coeffs); // integral along seg2 is equal to value at the center of SEG2 !
3546 * Low static method that operates 3D rotation of \a nbNodes 3D nodes whose coordinates are arranged in \a coords
3547 * around the center point \a center and with angle \a angle.
3549 void DataArrayDouble::Rotate2DAlg(const double *center, double angle, mcIdType nbNodes, const double *coordsIn, double *coordsOut)
3551 double cosa=cos(angle);
3552 double sina=sin(angle);
3554 matrix[0]=cosa; matrix[1]=-sina; matrix[2]=sina; matrix[3]=cosa;
3556 for(mcIdType i=0; i<nbNodes; i++)
3558 std::transform(coordsIn+i*2,coordsIn+(i+1)*2,center,tmp,std::minus<double>());
3559 coordsOut[i*2]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+center[0];
3560 coordsOut[i*2+1]=matrix[2]*tmp[0]+matrix[3]*tmp[1]+center[1];
3564 DataArrayDoubleIterator::DataArrayDoubleIterator(DataArrayDouble *da):DataArrayIterator<double>(da)
3568 DataArrayDoubleTuple::DataArrayDoubleTuple(double *pt, std::size_t nbOfComp):DataArrayTuple<double>(pt,nbOfComp)
3573 std::string DataArrayDoubleTuple::repr() const
3575 std::ostringstream oss; oss.precision(17); oss << "(";
3576 for(std::size_t i=0;i<_nb_of_compo-1;i++)
3577 oss << _pt[i] << ", ";
3578 oss << _pt[_nb_of_compo-1] << ")";
3582 double DataArrayDoubleTuple::doubleValue() const
3584 return this->zeValue();
3588 * This method returns a newly allocated instance the caller should dealed with by a MEDCoupling::DataArrayDouble::decrRef.
3589 * This method performs \b no copy of data. The content is only referenced using MEDCoupling::DataArrayDouble::useArray with ownership set to \b false.
3590 * 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
3591 * \b nbOfCompo=1 and \bnbOfTuples==this->_nb_of_elem.
3593 DataArrayDouble *DataArrayDoubleTuple::buildDADouble(std::size_t nbOfTuples, std::size_t nbOfCompo) const
3595 return this->buildDA(nbOfTuples,nbOfCompo);
3599 * Returns a full copy of \a this. For more info on copying data arrays see
3600 * \ref MEDCouplingArrayBasicsCopyDeep.
3601 * \return DataArrayInt * - a new instance of DataArrayInt.
3603 DataArrayInt32 *DataArrayInt32::deepCopy() const
3605 return new DataArrayInt32(*this);
3608 DataArrayInt32Iterator *DataArrayInt32::iterator()
3610 return new DataArrayInt32Iterator(this);
3614 DataArrayInt32Iterator::DataArrayInt32Iterator(DataArrayInt32 *da):DataArrayIterator<Int32>(da)
3618 DataArrayInt32Tuple::DataArrayInt32Tuple(Int32 *pt, std::size_t nbOfComp):DataArrayTuple<Int32>(pt,nbOfComp)
3622 std::string DataArrayInt32Tuple::repr() const
3624 std::ostringstream oss; oss << "(";
3625 for(std::size_t i=0;i<_nb_of_compo-1;i++)
3626 oss << _pt[i] << ", ";
3627 oss << _pt[_nb_of_compo-1] << ")";
3631 Int32 DataArrayInt32Tuple::intValue() const
3633 return this->zeValue();
3637 * This method returns a newly allocated instance the caller should dealed with by a MEDCoupling::DataArrayInt::decrRef.
3638 * This method performs \b no copy of data. The content is only referenced using MEDCoupling::DataArrayInt::useArray with ownership set to \b false.
3639 * 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
3640 * \b nbOfCompo=1 and \bnbOfTuples==this->_nb_of_elem.
3642 DataArrayInt32 *DataArrayInt32Tuple::buildDAInt(std::size_t nbOfTuples, std::size_t nbOfCompo) const
3644 return this->buildDA(nbOfTuples,nbOfCompo);
3647 DataArrayInt64Iterator *DataArrayInt64::iterator()
3649 return new DataArrayInt64Iterator(this);
3653 DataArrayInt64Iterator::DataArrayInt64Iterator(DataArrayInt64 *da):DataArrayIterator<Int64>(da)
3657 DataArrayInt64Tuple::DataArrayInt64Tuple(Int64 *pt, std::size_t nbOfComp):DataArrayTuple<Int64>(pt,nbOfComp)
3661 std::string DataArrayInt64Tuple::repr() const
3663 std::ostringstream oss; oss << "(";
3664 for(std::size_t i=0;i<_nb_of_compo-1;i++)
3665 oss << _pt[i] << ", ";
3666 oss << _pt[_nb_of_compo-1] << ")";
3670 Int64 DataArrayInt64Tuple::intValue() const
3672 return this->zeValue();
3675 DataArrayInt64 *DataArrayInt64Tuple::buildDAInt(std::size_t nbOfTuples, std::size_t nbOfCompo) const
3677 return this->buildDA(nbOfTuples,nbOfCompo);
3681 DataArrayInt64 *DataArrayInt64::deepCopy() const
3683 return new DataArrayInt64(*this);