// Author : Anthony Geay (CEA/DEN)
#include "MEDCouplingMemArray.txx"
-#include "MCAuto.hxx"
#include "BBTree.txx"
#include "GenMathFormulae.hxx"
using namespace MEDCoupling;
+template class MemArray<int>;
+template class MemArray<double>;
template class DataArrayTemplate<int>;
template class DataArrayTemplate<double>;
}
}
+int DataArray::EffectiveCircPerm(int nbOfShift, int nbOfTuples)
+{
+ if(nbOfTuples<=0)
+ throw INTERP_KERNEL::Exception("DataArray::EffectiveCircPerm : number of tuples is expected to be > 0 !");
+ if(nbOfShift>=0)
+ {
+ return nbOfShift%nbOfTuples;
+ }
+ else
+ {
+ int tmp(-nbOfShift);
+ tmp=tmp%nbOfTuples;
+ return nbOfTuples-tmp;
+ }
+}
+
std::size_t DataArray::getHeapMemorySizeWithoutChildren() const
{
std::size_t sz1=_name.capacity();
* is to delete this array using decrRef() as it is no more needed. The array
* does not contain any textual info on components.
* \throw If \a this->getNumberOfComponents() != 2.
+ * \sa fromCartToPolar
*/
DataArrayDouble *DataArrayDouble::fromPolarToCart() const
{
* on the third component is copied from \a this array. The caller
* is to delete this array using decrRef() as it is no more needed.
* \throw If \a this->getNumberOfComponents() != 3.
+ * \sa fromCartToCyl
*/
DataArrayDouble *DataArrayDouble::fromCylToCart() const
{
* on the third component is copied from \a this array. The caller
* is to delete this array using decrRef() as it is no more needed.
* \throw If \a this->getNumberOfComponents() != 3.
+ * \sa fromCartToSpher
*/
DataArrayDouble *DataArrayDouble::fromSpherToCart() const
{
return ret.retn();
}
+/*!
+ * This method returns a newly created array to be deallocated that contains the result of conversion from cartesian to polar.
+ * This method expects that \a this has exactly 2 components.
+ * \sa fromPolarToCart
+ */
+DataArrayDouble *DataArrayDouble::fromCartToPolar() const
+{
+ MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
+ checkAllocated();
+ int nbOfComp(getNumberOfComponents()),nbTuples(getNumberOfTuples());
+ if(nbOfComp!=2)
+ throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToPolar : must be an array with exactly 2 components !");
+ ret->alloc(nbTuples,2);
+ double *retPtr(ret->getPointer());
+ const double *ptr(begin());
+ for(int i=0;i<nbTuples;i++,ptr+=2,retPtr+=2)
+ {
+ retPtr[0]=sqrt(ptr[0]*ptr[0]+ptr[1]*ptr[1]);
+ retPtr[1]=atan2(ptr[1],ptr[0]);
+ }
+ return ret.retn();
+}
+
+/*!
+ * This method returns a newly created array to be deallocated that contains the result of conversion from cartesian to cylindrical.
+ * This method expects that \a this has exactly 3 components.
+ * \sa fromCylToCart
+ */
+DataArrayDouble *DataArrayDouble::fromCartToCyl() const
+{
+ MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
+ checkAllocated();
+ int nbOfComp(getNumberOfComponents()),nbTuples(getNumberOfTuples());
+ if(nbOfComp!=3)
+ throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCyl : must be an array with exactly 3 components !");
+ ret->alloc(nbTuples,3);
+ double *retPtr(ret->getPointer());
+ const double *ptr(begin());
+ for(int i=0;i<nbTuples;i++,ptr+=3,retPtr+=3)
+ {
+ retPtr[0]=sqrt(ptr[0]*ptr[0]+ptr[1]*ptr[1]);
+ retPtr[1]=atan2(ptr[1],ptr[0]);
+ retPtr[2]=ptr[2];
+ }
+ return ret.retn();
+}
+
+/*!
+ * This method returns a newly created array to be deallocated that contains the result of conversion from cartesian to spherical coordinates.
+ * \sa fromSpherToCart
+ */
+DataArrayDouble *DataArrayDouble::fromCartToSpher() const
+{
+ MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
+ checkAllocated();
+ int nbOfComp(getNumberOfComponents()),nbTuples(getNumberOfTuples());
+ if(nbOfComp!=3)
+ throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToSpher : must be an array with exactly 3 components !");
+ ret->alloc(nbTuples,3);
+ double *retPtr(ret->getPointer());
+ const double *ptr(begin());
+ for(int i=0;i<nbTuples;i++,ptr+=3,retPtr+=3)
+ {
+ retPtr[0]=sqrt(ptr[0]*ptr[0]+ptr[1]*ptr[1]+ptr[2]*ptr[2]);
+ retPtr[1]=acos(ptr[2]/retPtr[0]);
+ retPtr[2]=atan2(ptr[1],ptr[0]);
+ }
+ return ret.retn();
+}
+
+/*!
+ * 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.
+ * This method expects that \a this has exactly 3 components.
+ * \sa MEDCouplingFieldDouble::computeVectorFieldCyl
+ */
+DataArrayDouble *DataArrayDouble::fromCartToCylGiven(const DataArrayDouble *coords, const double center[3], const double vect[3]) const
+{
+ if(!coords)
+ throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCylGiven : input coords are NULL !");
+ MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
+ checkAllocated(); coords->checkAllocated();
+ int nbOfComp(getNumberOfComponents()),nbTuples(getNumberOfTuples());
+ if(nbOfComp!=3)
+ throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCylGiven : must be an array with exactly 3 components !");
+ if(coords->getNumberOfComponents()!=3)
+ throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCylGiven : coords array must have exactly 3 components !");
+ if(coords->getNumberOfTuples()!=nbTuples)
+ throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCylGiven : coords array must have the same number of tuples !");
+ ret->alloc(nbTuples,nbOfComp);
+ double magOfVect(sqrt(vect[0]*vect[0]+vect[1]*vect[1]+vect[2]*vect[2]));
+ if(magOfVect<1e-12)
+ throw INTERP_KERNEL::Exception("DataArrayDouble::fromCartToCylGiven : magnitude of vect is too low !");
+ double Ur[3],Uteta[3],Uz[3],*retPtr(ret->getPointer());
+ const double *coo(coords->begin()),*vectField(begin());
+ std::transform(vect,vect+3,Uz,std::bind2nd(std::multiplies<double>(),1./magOfVect));
+ for(int i=0;i<nbTuples;i++,vectField+=3,retPtr+=3,coo+=3)
+ {
+ std::transform(coo,coo+3,center,Ur,std::minus<double>());
+ 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];
+ double magOfTeta(sqrt(Uteta[0]*Uteta[0]+Uteta[1]*Uteta[1]+Uteta[2]*Uteta[2]));
+ std::transform(Uteta,Uteta+3,Uteta,std::bind2nd(std::multiplies<double>(),1./magOfTeta));
+ 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];
+ retPtr[0]=Ur[0]*vectField[0]+Ur[1]*vectField[1]+Ur[2]*vectField[2];
+ retPtr[1]=Uteta[0]*vectField[0]+Uteta[1]*vectField[1]+Uteta[2]*vectField[2];
+ retPtr[2]=Uz[0]*vectField[0]+Uz[1]*vectField[1]+Uz[2]*vectField[2];
+ }
+ ret->copyStringInfoFrom(*this);
+ return ret.retn();
+}
+
/*!
* Computes the doubly contracted product of every tensor defined by the tuple of \a this
* array contating 6 components.
declareAsNew();
}
+/*!
+ * \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.
+ */
+MCAuto<DataArrayDouble> DataArrayDouble::symmetry3DPlane(const double point[3], const double normalVector[3]) const
+{
+ checkAllocated();
+ if(getNumberOfComponents()!=3)
+ throw INTERP_KERNEL::Exception("DataArrayDouble::symmetry3DPlane : this is excepted to have 3 components !");
+ int nbTuples(getNumberOfTuples());
+ MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
+ ret->alloc(nbTuples,3);
+ Symmetry3DPlane(point,normalVector,nbTuples,begin(),ret->getPointer());
+ return ret;
+}
+
DataArrayDoubleIterator *DataArrayDouble::iterator()
{
return new DataArrayDoubleIterator(this);
}
}
+/*!
+ * Low static method that operates 3D rotation of 'nbNodes' 3D nodes whose coordinates are arranged in \a coordsIn
+ * around an axe ( \a center, \a vect) and with angle \a angle.
+ */
+void DataArrayDouble::Rotate3DAlg(const double *center, const double *vect, double angle, int nbNodes, const double *coordsIn, double *coordsOut)
+{
+ if(!center || !vect)
+ throw INTERP_KERNEL::Exception("DataArrayDouble::Rotate3DAlg : null vector in input !");
+ double sina(sin(angle));
+ double cosa(cos(angle));
+ double vectorNorm[3];
+ double matrix[9];
+ double matrixTmp[9];
+ double norm(sqrt(vect[0]*vect[0]+vect[1]*vect[1]+vect[2]*vect[2]));
+ if(norm<std::numeric_limits<double>::min())
+ throw INTERP_KERNEL::Exception("DataArrayDouble::Rotate3DAlg : magnitude of input vector is too close of 0. !");
+ std::transform(vect,vect+3,vectorNorm,std::bind2nd(std::multiplies<double>(),1/norm));
+ //rotation matrix computation
+ 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;
+ matrixTmp[0]=vectorNorm[0]*vectorNorm[0]; matrixTmp[1]=vectorNorm[0]*vectorNorm[1]; matrixTmp[2]=vectorNorm[0]*vectorNorm[2];
+ matrixTmp[3]=vectorNorm[1]*vectorNorm[0]; matrixTmp[4]=vectorNorm[1]*vectorNorm[1]; matrixTmp[5]=vectorNorm[1]*vectorNorm[2];
+ matrixTmp[6]=vectorNorm[2]*vectorNorm[0]; matrixTmp[7]=vectorNorm[2]*vectorNorm[1]; matrixTmp[8]=vectorNorm[2]*vectorNorm[2];
+ std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies<double>(),1-cosa));
+ std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
+ matrixTmp[0]=0.; matrixTmp[1]=-vectorNorm[2]; matrixTmp[2]=vectorNorm[1];
+ matrixTmp[3]=vectorNorm[2]; matrixTmp[4]=0.; matrixTmp[5]=-vectorNorm[0];
+ matrixTmp[6]=-vectorNorm[1]; matrixTmp[7]=vectorNorm[0]; matrixTmp[8]=0.;
+ std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies<double>(),sina));
+ std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
+ //rotation matrix computed.
+ double tmp[3];
+ for(int i=0; i<nbNodes; i++)
+ {
+ std::transform(coordsIn+i*3,coordsIn+(i+1)*3,center,tmp,std::minus<double>());
+ coordsOut[i*3]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+matrix[2]*tmp[2]+center[0];
+ coordsOut[i*3+1]=matrix[3]*tmp[0]+matrix[4]*tmp[1]+matrix[5]*tmp[2]+center[1];
+ coordsOut[i*3+2]=matrix[6]*tmp[0]+matrix[7]*tmp[1]+matrix[8]*tmp[2]+center[2];
+ }
+}
+
+void DataArrayDouble::Symmetry3DPlane(const double point[3], const double normalVector[3], int nbNodes, const double *coordsIn, double *coordsOut)
+{
+ double matrix[9],matrix2[9],matrix3[9];
+ double vect[3],crossVect[3];
+ INTERP_KERNEL::orthogonalVect3(normalVector,vect);
+ crossVect[0]=normalVector[1]*vect[2]-normalVector[2]*vect[1];
+ crossVect[1]=normalVector[2]*vect[0]-normalVector[0]*vect[2];
+ crossVect[2]=normalVector[0]*vect[1]-normalVector[1]*vect[0];
+ double nv(INTERP_KERNEL::norm<3>(vect)),ni(INTERP_KERNEL::norm<3>(normalVector)),nc(INTERP_KERNEL::norm<3>(crossVect));
+ matrix[0]=vect[0]/nv; matrix[1]=crossVect[0]/nc; matrix[2]=-normalVector[0]/ni;
+ matrix[3]=vect[1]/nv; matrix[4]=crossVect[1]/nc; matrix[5]=-normalVector[1]/ni;
+ matrix[6]=vect[2]/nv; matrix[7]=crossVect[2]/nc; matrix[8]=-normalVector[2]/ni;
+ matrix2[0]=vect[0]/nv; matrix2[1]=vect[1]/nv; matrix2[2]=vect[2]/nv;
+ matrix2[3]=crossVect[0]/nc; matrix2[4]=crossVect[1]/nc; matrix2[5]=crossVect[2]/nc;
+ matrix2[6]=normalVector[0]/ni; matrix2[7]=normalVector[1]/ni; matrix2[8]=normalVector[2]/ni;
+ for(int i=0;i<3;i++)
+ for(int j=0;j<3;j++)
+ {
+ double val(0.);
+ for(int k=0;k<3;k++)
+ val+=matrix[3*i+k]*matrix2[3*k+j];
+ matrix3[3*i+j]=val;
+ }
+ //rotation matrix computed.
+ double tmp[3];
+ for(int i=0; i<nbNodes; i++)
+ {
+ std::transform(coordsIn+i*3,coordsIn+(i+1)*3,point,tmp,std::minus<double>());
+ coordsOut[i*3]=matrix3[0]*tmp[0]+matrix3[1]*tmp[1]+matrix3[2]*tmp[2]+point[0];
+ coordsOut[i*3+1]=matrix3[3]*tmp[0]+matrix3[4]*tmp[1]+matrix3[5]*tmp[2]+point[1];
+ coordsOut[i*3+2]=matrix3[6]*tmp[0]+matrix3[7]*tmp[1]+matrix3[8]*tmp[2]+point[2];
+ }
+}
+
+void DataArrayDouble::GiveBaseForPlane(const double normalVector[3], double baseOfPlane[9])
+{
+ double vect[3],crossVect[3];
+ INTERP_KERNEL::orthogonalVect3(normalVector,vect);
+ crossVect[0]=normalVector[1]*vect[2]-normalVector[2]*vect[1];
+ crossVect[1]=normalVector[2]*vect[0]-normalVector[0]*vect[2];
+ crossVect[2]=normalVector[0]*vect[1]-normalVector[1]*vect[0];
+ double nv(INTERP_KERNEL::norm<3>(vect)),ni(INTERP_KERNEL::norm<3>(normalVector)),nc(INTERP_KERNEL::norm<3>(crossVect));
+ baseOfPlane[0]=vect[0]/nv; baseOfPlane[1]=vect[1]/nv; baseOfPlane[2]=vect[2]/nv;
+ baseOfPlane[3]=crossVect[0]/nc; baseOfPlane[4]=crossVect[1]/nc; baseOfPlane[5]=crossVect[2]/nc;
+ baseOfPlane[6]=normalVector[0]/ni; baseOfPlane[7]=normalVector[1]/ni; baseOfPlane[8]=normalVector[2]/ni;
+}
+
+/*!
+ * Low static method that operates 3D rotation of \a nbNodes 3D nodes whose coordinates are arranged in \a coords
+ * around the center point \a center and with angle \a angle.
+ */
+void DataArrayDouble::Rotate2DAlg(const double *center, double angle, int nbNodes, const double *coordsIn, double *coordsOut)
+{
+ double cosa=cos(angle);
+ double sina=sin(angle);
+ double matrix[4];
+ matrix[0]=cosa; matrix[1]=-sina; matrix[2]=sina; matrix[3]=cosa;
+ double tmp[2];
+ for(int i=0; i<nbNodes; i++)
+ {
+ std::transform(coordsIn+i*2,coordsIn+(i+1)*2,center,tmp,std::minus<double>());
+ coordsOut[i*2]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+center[0];
+ coordsOut[i*2+1]=matrix[2]*tmp[0]+matrix[3]*tmp[1]+center[1];
+ }
+}
+
DataArrayDoubleIterator::DataArrayDoubleIterator(DataArrayDouble *da):_da(da),_tuple_id(0),_nb_comp(0),_nb_tuple(0)
{
if(_da)
return ret.retn();
}
+/*!
+ * Elements of \a partOfThis are expected to be included in \a this.
+ * The returned array \a ret is so that this[ret]==partOfThis
+ *
+ * For example, if \a this array contents are [9,10,0,6,4,11,3,8] and if \a partOfThis contains [6,0,11,8]
+ * the return array will contain [3,2,5,7].
+ *
+ * \a this is expected to be a 1 compo allocated array.
+ * \param [in] partOfThis - A 1 compo allocated array
+ * \return - A newly allocated array to be dealed by caller having the same number of tuples than \a partOfThis.
+ * \throw if two same element is present twice in \a this
+ * \throw if an element in \a partOfThis is \b NOT in \a this.
+ */
+DataArrayInt *DataArrayInt::indicesOfSubPart(const DataArrayInt& partOfThis) const
+{
+ if(getNumberOfComponents()!=1 || partOfThis.getNumberOfComponents()!=1)
+ throw INTERP_KERNEL::Exception("DataArrayInt::indicesOfSubPart : this and input array must be one component array !");
+ checkAllocated(); partOfThis.checkAllocated();
+ int thisNbTuples(getNumberOfTuples()),nbTuples(partOfThis.getNumberOfTuples());
+ const int *thisPt(begin()),*pt(partOfThis.begin());
+ MCAuto<DataArrayInt> ret(DataArrayInt::New());
+ ret->alloc(nbTuples,1);
+ int *retPt(ret->getPointer());
+ std::map<int,int> m;
+ for(int i=0;i<thisNbTuples;i++,thisPt++)
+ m[*thisPt]=i;
+ if(m.size()!=thisNbTuples)
+ throw INTERP_KERNEL::Exception("DataArrayInt::indicesOfSubPart : some elements appears more than once !");
+ for(int i=0;i<nbTuples;i++,retPt++,pt++)
+ {
+ std::map<int,int>::const_iterator it(m.find(*pt));
+ if(it!=m.end())
+ *retPt=(*it).second;
+ else
+ {
+ std::ostringstream oss; oss << "DataArrayInt::indicesOfSubPart : At pos #" << i << " of input array value is " << *pt << " not in this !";
+ throw INTERP_KERNEL::Exception(oss.str());
+ }
+ }
+ return ret.retn();
+}
+
void DataArrayInt::aggregate(const DataArrayInt *other)
{
if(!other)
}
/*!
- * Checks if contents of \a this array are equal to that of an array filled with
+ * Checks if \a this array has the given size, and if its contents is equal to an array filled with
* iota(). This method is particularly useful for DataArrayInt instances that represent
- * a renumbering array to check the real need in renumbering. This method checks than \a this can be considered as an identity function
+ * a renumbering array, to check if there is a real need in renumbering.
+ * This method checks than \a this can be considered as an identity mapping
* of a set having \a sizeExpected elements into itself.
*
* \param [in] sizeExpected - The number of elements expected.
