From 6c20f278d994697d0b6c3732d0b45ddcf8f503f8 Mon Sep 17 00:00:00 2001 From: ageay Date: Mon, 16 Aug 2010 12:38:38 +0000 Subject: [PATCH] Addition of functionalities in Fields. --- src/INTERP_KERNEL/VolSurfFormulae.hxx | 137 +++++- src/INTERP_KERNEL/VolSurfUser.txx | 71 ++- .../MEDCouplingAutoRefCountObjectPtr.hxx | 49 ++ src/MEDCoupling/MEDCouplingCMesh.cxx | 64 +++ src/MEDCoupling/MEDCouplingCMesh.hxx | 7 + src/MEDCoupling/MEDCouplingExtrudedMesh.cxx | 30 ++ src/MEDCoupling/MEDCouplingExtrudedMesh.hxx | 5 + src/MEDCoupling/MEDCouplingField.cxx | 22 +- src/MEDCoupling/MEDCouplingField.hxx | 3 +- .../MEDCouplingFieldDiscretization.cxx | 223 ++++++++- .../MEDCouplingFieldDiscretization.hxx | 36 +- src/MEDCoupling/MEDCouplingFieldDouble.cxx | 238 ++++++++-- src/MEDCoupling/MEDCouplingFieldDouble.hxx | 13 +- .../MEDCouplingGaussLocalization.hxx | 4 +- src/MEDCoupling/MEDCouplingMemArray.cxx | 54 +++ src/MEDCoupling/MEDCouplingMemArray.hxx | 4 + src/MEDCoupling/MEDCouplingMesh.cxx | 74 ++- src/MEDCoupling/MEDCouplingMesh.hxx | 12 +- src/MEDCoupling/MEDCouplingNatureOfField.hxx | 4 +- .../MEDCouplingNormalizedCartesianMesh.hxx | 4 +- .../MEDCouplingNormalizedUnstructuredMesh.hxx | 4 +- src/MEDCoupling/MEDCouplingPointSet.cxx | 14 + src/MEDCoupling/MEDCouplingPointSet.hxx | 1 + .../MEDCouplingTimeDiscretization.cxx | 86 +++- .../MEDCouplingTimeDiscretization.hxx | 20 +- src/MEDCoupling/MEDCouplingUMesh.cxx | 179 ++++++- src/MEDCoupling/MEDCouplingUMesh.hxx | 18 +- src/MEDCoupling/MEDCouplingUMeshDesc.cxx | 19 + src/MEDCoupling/MEDCouplingUMeshDesc.hxx | 4 + src/MEDCoupling/Makefile.am | 4 +- .../Test/MEDCouplingBasicsTest.hxx | 16 + .../Test/MEDCouplingBasicsTest0.cxx | 71 +++ .../Test/MEDCouplingBasicsTest1.cxx | 4 +- .../Test/MEDCouplingBasicsTest2.cxx | 442 ++++++++++++++++++ src/MEDCoupling_Swig/MEDCouplingBasicsTest.py | 379 ++++++++++++++- .../MEDCouplingDataForTest.py | 55 +++ src/MEDCoupling_Swig/libMEDCoupling_Swig.i | 80 +++- src/ParaMEDMEM/ParaFIELD.cxx | 2 +- 38 files changed, 2272 insertions(+), 180 deletions(-) create mode 100644 src/MEDCoupling/MEDCouplingAutoRefCountObjectPtr.hxx diff --git a/src/INTERP_KERNEL/VolSurfFormulae.hxx b/src/INTERP_KERNEL/VolSurfFormulae.hxx index 09c8e54ad..f486a9a8d 100644 --- a/src/INTERP_KERNEL/VolSurfFormulae.hxx +++ b/src/INTERP_KERNEL/VolSurfFormulae.hxx @@ -22,7 +22,7 @@ #include "InterpolationUtils.hxx" -#include +#include namespace INTERP_KERNEL { @@ -36,7 +36,7 @@ namespace INTERP_KERNEL inline double calculateLgthForSeg2(const double *p1, const double *p2, int spaceDim) { if(spaceDim==1) - return fabs(*p2-*p1); + return *p2-*p1; else { double ret=0; @@ -415,7 +415,8 @@ namespace INTERP_KERNEL /*! * Calculate Volume for Generic Polyedron, even not convex one, WARNING !!! The polyedron's faces must be correctly ordered. - * 2nd API avoiding to create double** arrays. + * 2nd API avoiding to create double** arrays. The returned value could be negative if polyhedrons faces are not oriented with normal outside of the + * polyhedron */ template inline double calculateVolumeForPolyh2(const ConnType *connec, int lgth, const double *coords) @@ -443,6 +444,94 @@ namespace INTERP_KERNEL return volume/6.; } + /*! + * This method returns the area oriented vector of a polygon. This method is useful for normal computation without + * any troubles if several edges are colinears. + * @param res must be of size at least 3 to store the result. + */ + template + inline double areaVectorOfPolygon(const ConnType *connec, int lgth, const double *coords, double *res) + { + res[0]=0.; res[1]=0.; res[2]=0.; + for(int ptId=0;ptId::coo2C(connec[ptId]); + const double *pti1=coords+3*OTT::coo2C(connec[(ptId+1)%lgth]); + res[0]+=pti[1]*pti1[2]-pti[2]*pti1[1]; + res[1]+=pti[2]*pti1[0]-pti[0]*pti1[2]; + res[2]+=pti[0]*pti1[1]-pti[1]*pti1[0]; + } + } + + inline double integrationOverA3DLine(double u1, double v1, double u2, double v2, double A, double B, double C, double J, double K) + { + return (-u2*u2*u2*u2*J*J*J*B*B+12*u1*C*C*K-4*u2*u2*u2*A*A*K+3*u1*u1*u1*u1*J*J*A*B-6*C*C*u2*u2*J-4*C*u2*u2*u2*J*J*B-8*C*u2*u2*u2*J*A-3*u2*u2*u2*u2*J*A*A+8*u1*u1*u1*J*A*K*B-4*u2*K*K*K*B*B-3*u2*u2*u2*u2*J*J*A*B-12*C*u2*u2*J*K*B+6*u1*u1*J*K*K*B*B+6*u1*u1*A*K*K*B+4*u1*u1*u1*A*A*K+u1*u1*u1*u1*J*J*J*B*B-8*u2*u2*u2*J*A*K*B-12*C*C*u2*K+12*u1*u1*C*A*K-4*u2*u2*u2*J*J*K*B*B+6*u1*u1*C*C*J+4*u1*u1*u1*J*J*K*B*B+4*u1*u1*u1*C*J*J*B+8*u1*u1*u1*C*J*A+12*u1*u1*C*J*K*B+12*u1*C*K*K*B+4*u1*K*K*K*B*B-12*C*u2*u2*A*K+3*u1*u1*u1*u1*J*A*A-12*C*u2*K*K*B-6*u2*u2*J*K*K*B*B-6*u2*u2*A*K*K*B)/24.; + } + + template + inline void barycenterOfPolyhedron(const ConnType *connec, int lgth, const double *coords, double *res) + { + int nbOfFaces=std::count(connec,connec+lgth,-1)+1; + res[0]=0.; res[1]=0.; res[2]=0.; + const int *work=connec; + for(int i=0;i(work,nbOfNodesOfCurFace,coords,normal); + double normOfNormal=sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]); + normal[0]/=normOfNormal; normal[1]/=normOfNormal; normal[2]/=normOfNormal; + double u[2]={normal[1],-normal[0]}; + double s=sqrt(u[0]*u[0]+u[1]*u[1]); + double c=normal[2]; + if(fabs(s)>1e-12) + { + u[0]/=std::abs(s); u[1]/=std::abs(s); + } + else + { u[0]=1.; u[1]=0.; } + //C : high in plane of polyhedron face : always constant + double w=normal[0]*coords[3*OTT::coo2C(work[0])]+ + normal[1]*coords[3*OTT::coo2C(work[0])+1]+ + normal[2]*coords[3*OTT::coo2C(work[0])+2]; + // A,B,D,F,G,H,L,M,N coeffs of rotation matrix defined by (u,c,s) + double A=u[0]*u[0]*(1-c)+c; + double B=u[0]*u[1]*(1-c); + double D=u[1]*s; + double F=B; + double G=u[1]*u[1]*(1-c)+c; + double H=-u[0]*s; + double L=-u[1]*s; + double M=u[0]*s; + double N=c; + for(int j=0;j::coo2C(work[j]); + const double *p2=coords+3*OTT::coo2C(work[(j+1)%nbOfNodesOfCurFace]); + double u1=A*p1[0]+B*p1[1]+D*p1[2]; + double u2=A*p2[0]+B*p2[1]+D*p2[2]; + double v1=F*p1[0]+G*p1[1]+H*p1[2]; + double v2=F*p2[0]+G*p2[1]+H*p2[2]; + if(fabs(u2-u1)>1e-12) + { + double J=(v2-v1)/(u2-u1); + double K=(v1*u2-v2*u1)/(u2-u1); + double gx=integrationOverA3DLine(u1,v1,u2,v2,A,B,-w*D,J,K); + double gy=integrationOverA3DLine(u1,v1,u2,v2,F,G,-w*H,J,K); + double gz=integrationOverA3DLine(u1,v1,u2,v2,L,M,-w*N,J,K); + res[0]+=gx*normal[0]; + res[1]+=gy*normal[1]; + res[2]+=gz*normal[2]; + } + } + work=work2+1; + } + double vol=calculateVolumeForPolyh2(connec,lgth,coords); + res[0]/=vol; res[1]/=vol; res[2]/=vol; + } + // ============================================================================================================================================ // Calculate Volume for NON Generic Polyedron. Only polydrons with bary included in pts is supported by this method. Result is always positive. // ============================================================================================================================================ @@ -556,7 +645,7 @@ namespace INTERP_KERNEL res[0]=0.; res[1]=0.; for(int i=0;i::coo2C(connec[i])]*coords[2*OTT::coo2C(connec[(i+1)%lgth])+1]+ + double cp=coords[2*OTT::coo2C(connec[i])]*coords[2*OTT::coo2C(connec[(i+1)%lgth])+1]- coords[2*OTT::coo2C(connec[i])+1]*coords[2*OTT::coo2C(connec[(i+1)%lgth])]; area+=cp; res[0]+=cp*(coords[2*OTT::coo2C(connec[i])]+coords[2*OTT::coo2C(connec[(i+1)%lgth])]); @@ -566,23 +655,37 @@ namespace INTERP_KERNEL res[1]/=3.*area; } - template + template inline void computePolygonBarycenter3D(const ConnType *connec, int lgth, const double *coords, double *res) { - double area[3]={0.,0.,0.}; + double area[3]; + areaVectorOfPolygon(connec,lgth,coords,area); + double norm=sqrt(area[0]*area[0]+area[1]*area[1]+area[2]*area[2]); + area[0]/=norm; area[1]/=norm; area[2]/=norm; res[0]=0.; res[1]=0.; res[2]=0.; - for(int i=0;i::coo2C(connec[0])]+ + coords[3*OTT::coo2C(connec[i])]+ + coords[3*OTT::coo2C(connec[i+1])])/3.; + v[1]=(coords[3*OTT::coo2C(connec[0])+1]+ + coords[3*OTT::coo2C(connec[i])+1]+ + coords[3*OTT::coo2C(connec[i+1])+1])/3.; + v[2]=(coords[3*OTT::coo2C(connec[0])+2]+ + coords[3*OTT::coo2C(connec[i])+2]+ + coords[3*OTT::coo2C(connec[i+1])+2])/3.; + ConnType tmpConn[3]={connec[0],connec[i],connec[i+1]}; + areaVectorOfPolygon(tmpConn,3,coords,tmpArea); + double norm2=sqrt(tmpArea[0]*tmpArea[0]+tmpArea[1]*tmpArea[1]+tmpArea[2]*tmpArea[2]); + if(norm2>1e-12) + { + tmpArea[0]/=norm2; tmpArea[1]/=norm2; tmpArea[2]/=norm2; + double signOfArea=area[0]*tmpArea[0]+area[1]*tmpArea[1]+area[2]*tmpArea[2]; + res[0]+=signOfArea*norm2*v[0]/norm; res[1]+=signOfArea*norm2*v[1]/norm; res[2]+=signOfArea*norm2*v[2]/norm; + } + } } } diff --git a/src/INTERP_KERNEL/VolSurfUser.txx b/src/INTERP_KERNEL/VolSurfUser.txx index 1273d6950..5a0109a43 100644 --- a/src/INTERP_KERNEL/VolSurfUser.txx +++ b/src/INTERP_KERNEL/VolSurfUser.txx @@ -73,7 +73,7 @@ namespace INTERP_KERNEL { const double **pts=new const double *[lgth]; for(int inod=0;inod::coo2C(connec[inod]); + pts[inod] = coords+SPACEDIM*OTT::coo2C(connec[inod]); double val=INTERP_KERNEL::calculateAreaForPolyg(pts,lgth,SPACEDIM); delete [] pts; return val; @@ -175,23 +175,82 @@ namespace INTERP_KERNEL throw INTERP_KERNEL::Exception("Invalid spaceDim specified : must be 1, 2 or 3"); } - template + template void computeBarycenter(NormalizedCellType type, const ConnType *connec, int lgth, const double *coords, double *res) { switch(type) { + case NORM_SEG3: + case NORM_SEG2: + { + std::copy(coords+SPACEDIM*OTT::coo2C(connec[0]), + coords+SPACEDIM*OTT::coo2C(connec[0]+1),res); + std::transform(res,res+SPACEDIM,coords+SPACEDIM*OTT::coo2C(connec[1]),res,std::plus()); + std::transform(res,res+SPACEDIM,res,std::bind2nd(std::multiplies(),0.5)); + break; + } case NORM_TRI3: + case NORM_TRI6: + { + std::copy(coords+SPACEDIM*OTT::coo2C(connec[0]), + coords+SPACEDIM*OTT::coo2C(connec[0]+1),res); + std::transform(res,res+SPACEDIM,coords+SPACEDIM*OTT::coo2C(connec[1]),res,std::plus()); + std::transform(res,res+SPACEDIM,coords+SPACEDIM*OTT::coo2C(connec[2]),res,std::plus()); + std::transform(res,res+SPACEDIM,res,std::bind2nd(std::multiplies(),1./3.)); + break; + } case NORM_QUAD4: case NORM_POLYGON: { if(SPACEDIM==2) - computePolygonBarycenter2D(connec,lgth,coords,res); + computePolygonBarycenter2D(connec,lgth,coords,res); else if(SPACEDIM==3) - computePolygonBarycenter3D(connec,lgth,coords,res); + computePolygonBarycenter3D(connec,lgth,coords,res); else throw INTERP_KERNEL::Exception("Impossible spacedim linked to cell 2D Cell !"); break; } + case NORM_QUAD8: + { + if(SPACEDIM==2) + computePolygonBarycenter2D(connec,lgth/2,coords,res); + else if(SPACEDIM==3) + computePolygonBarycenter3D(connec,lgth/2,coords,res); + else + throw INTERP_KERNEL::Exception("Impossible spacedim linked to cell 2D Cell !"); + break; + } + case NORM_TETRA4: + { + res[0]=coords[3*OTT::coo2C(connec[0])]; + res[1]=coords[3*OTT::coo2C(connec[0])+1]; + res[2]=coords[3*OTT::coo2C(connec[0])+2]; + res[0]+=coords[3*OTT::coo2C(connec[1])]; + res[1]+=coords[3*OTT::coo2C(connec[1])+1]; + res[2]+=coords[3*OTT::coo2C(connec[1])+2]; + res[0]+=coords[3*OTT::coo2C(connec[2])]; + res[1]+=coords[3*OTT::coo2C(connec[2])+1]; + res[2]+=coords[3*OTT::coo2C(connec[2])+2]; + res[0]+=coords[3*OTT::coo2C(connec[3])]; + res[1]+=coords[3*OTT::coo2C(connec[3])+1]; + res[2]+=coords[3*OTT::coo2C(connec[3])+2]; + res[0]/=4.; res[1]/=4.; res[2]/=4.; + break; + } + case NORM_PYRA5: + { + double tmp[3]; + computePolygonBarycenter3D(connec,lgth,coords,tmp); + res[0]=(coords[3*OTT::coo2C(connec[4])]+3.*tmp[0])/4.; + res[1]=(coords[3*OTT::coo2C(connec[4])+1]+3.*tmp[1])/4.; + res[1]=(coords[3*OTT::coo2C(connec[4])+2]+3.*tmp[2])/4.; + break; + } + case NORM_POLYHED: + { + barycenterOfPolyhedron(connec,lgth,coords,res); + break; + } default: throw INTERP_KERNEL::Exception("Not recognized cell type to get Barycenter on it !"); } @@ -204,8 +263,8 @@ namespace INTERP_KERNEL return computeBarycenter(type,connec,lgth,coords,res); if(spaceDim==2) return computeBarycenter(type,connec,lgth,coords,res); - //if(spaceDim==1) - // return computeBarycenter(type,connec,lgth,coords,res); + if(spaceDim==1) + return computeBarycenter(type,connec,lgth,coords,res); throw INTERP_KERNEL::Exception("Invalid spaceDim specified for compute barycenter : must be 1, 2 or 3"); } } diff --git a/src/MEDCoupling/MEDCouplingAutoRefCountObjectPtr.hxx b/src/MEDCoupling/MEDCouplingAutoRefCountObjectPtr.hxx new file mode 100644 index 000000000..d91a8e16f --- /dev/null +++ b/src/MEDCoupling/MEDCouplingAutoRefCountObjectPtr.hxx @@ -0,0 +1,49 @@ +// Copyright (C) 2007-2010 CEA/DEN, EDF R&D +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// + +#ifndef __PARAMEDMEM_MEDCOUPLINGAUTOREFCOUNTOBJECTPTR_HXX__ +#define __PARAMEDMEM_MEDCOUPLINGAUTOREFCOUNTOBJECTPTR_HXX__ + +#include "MEDCouplingRefCountObject.hxx" + +namespace ParaMEDMEM +{ + template + class MEDCouplingAutoRefCountObjectPtr + { + public: + MEDCouplingAutoRefCountObjectPtr(const MEDCouplingAutoRefCountObjectPtr& other):_ptr(0) { referPtr(other._ptr); } + MEDCouplingAutoRefCountObjectPtr(T *ptr):_ptr(ptr) { } + ~MEDCouplingAutoRefCountObjectPtr() { destroyPtr(); } + MEDCouplingAutoRefCountObjectPtr &operator=(const MEDCouplingAutoRefCountObjectPtr& other) { destroyPtr(); referPtr(other._ptr); return *this; } + T *operator->() { return _ptr ; } + const T *operator->() const { return _ptr; } + T& operator*() { return *_ptr; } + const T& operator*() const { return *_ptr; } + operator T *() { return _ptr; } + operator const T *() const { return _ptr; } + private: + void referPtr(T *ptr) { _ptr=ptr; if(_ptr) _ptr->incrRef(); } + void destroyPtr() { if(_ptr) _ptr->decrRef(); } + private: + T *_ptr; + }; +} + +#endif diff --git a/src/MEDCoupling/MEDCouplingCMesh.cxx b/src/MEDCoupling/MEDCouplingCMesh.cxx index 452a3bf84..e75584697 100644 --- a/src/MEDCoupling/MEDCouplingCMesh.cxx +++ b/src/MEDCoupling/MEDCouplingCMesh.cxx @@ -31,6 +31,31 @@ MEDCouplingCMesh::MEDCouplingCMesh():_x_array(0),_y_array(0),_z_array(0) { } +MEDCouplingCMesh::MEDCouplingCMesh(const MEDCouplingCMesh& other, bool deepCpy) +{ + if(deepCpy) + { + if(other._x_array) + _x_array=_x_array->deepCopy(); + if(other._y_array) + _y_array=_y_array->deepCopy(); + if(other._z_array) + _z_array=_z_array->deepCopy(); + } + else + { + _x_array=other._x_array; + if(_x_array) + _x_array->incrRef(); + _y_array=other._y_array; + if(_y_array) + _y_array->incrRef(); + _z_array=other._z_array; + if(_z_array) + _z_array->incrRef(); + } +} + MEDCouplingCMesh::~MEDCouplingCMesh() { if(_x_array) @@ -46,6 +71,16 @@ MEDCouplingCMesh *MEDCouplingCMesh::New() return new MEDCouplingCMesh; } +MEDCouplingMesh *MEDCouplingCMesh::deepCpy() const +{ + return clone(true); +} + +MEDCouplingCMesh *MEDCouplingCMesh::clone(bool recDeepCpy) const +{ + return new MEDCouplingCMesh(*this,recDeepCpy); +} + void MEDCouplingCMesh::updateTime() { if(_x_array) @@ -56,6 +91,24 @@ void MEDCouplingCMesh::updateTime() updateTimeWith(*_z_array); } +/*! + * This method copyies all tiny strings from other (name and components name). + * @throw if other and this have not same mesh type. + */ +void MEDCouplingCMesh::copyTinyStringsFrom(const MEDCouplingMesh *other) throw(INTERP_KERNEL::Exception) +{ + const MEDCouplingCMesh *otherC=dynamic_cast(other); + if(!otherC) + throw INTERP_KERNEL::Exception("MEDCouplingCMesh::copyTinyStringsFrom : meshes have not same type !"); + MEDCouplingMesh::copyTinyStringsFrom(other); + if(_x_array && otherC->_x_array) + _x_array->copyStringInfoFrom(*otherC->_x_array); + if(_y_array && otherC->_y_array) + _y_array->copyStringInfoFrom(*otherC->_y_array); + if(_z_array && otherC->_z_array) + _z_array->copyStringInfoFrom(*otherC->_z_array); +} + bool MEDCouplingCMesh::isEqual(const MEDCouplingMesh *other, double prec) const { const MEDCouplingCMesh *otherC=dynamic_cast(other); @@ -64,6 +117,12 @@ bool MEDCouplingCMesh::isEqual(const MEDCouplingMesh *other, double prec) const return true; } +void MEDCouplingCMesh::checkDeepEquivalWith(const MEDCouplingMesh *other, int cellCompPol, double prec, + DataArrayInt *&cellCor, DataArrayInt *&nodeCor) const throw(INTERP_KERNEL::Exception) +{ + throw INTERP_KERNEL::Exception("Not implemented yet !"); +} + void MEDCouplingCMesh::checkCoherency() const throw(INTERP_KERNEL::Exception) { const char msg0[]="Invalid "; @@ -386,6 +445,11 @@ DataArrayDouble *MEDCouplingCMesh::getBarycenterAndOwner() const return ret; } +void MEDCouplingCMesh::renumberCells(const int *old2NewBg, const int *old2NewEnd, bool check) throw(INTERP_KERNEL::Exception) +{ + throw INTERP_KERNEL::Exception("Functionnality of renumbering cell not available for CMesh !"); +} + void MEDCouplingCMesh::getTinySerializationInformation(std::vector& tinyInfo, std::vector& littleStrings) const { throw INTERP_KERNEL::Exception("Not implemented yet !"); diff --git a/src/MEDCoupling/MEDCouplingCMesh.hxx b/src/MEDCoupling/MEDCouplingCMesh.hxx index e70dce3cd..eda2487c4 100644 --- a/src/MEDCoupling/MEDCouplingCMesh.hxx +++ b/src/MEDCoupling/MEDCouplingCMesh.hxx @@ -31,9 +31,14 @@ namespace ParaMEDMEM { public: static MEDCouplingCMesh *New(); + MEDCouplingMesh *deepCpy() const; + MEDCouplingCMesh *clone(bool recDeepCpy) const; void updateTime(); MEDCouplingMeshType getType() const { return CARTESIAN; } + void copyTinyStringsFrom(const MEDCouplingMesh *other) throw(INTERP_KERNEL::Exception); bool isEqual(const MEDCouplingMesh *other, double prec) const; + void checkDeepEquivalWith(const MEDCouplingMesh *other, int cellCompPol, double prec, + DataArrayInt *&cellCor, DataArrayInt *&nodeCor) const throw(INTERP_KERNEL::Exception); void checkCoherency() const throw(INTERP_KERNEL::Exception); bool isStructured() const; int getNumberOfCells() const; @@ -62,6 +67,7 @@ namespace ParaMEDMEM MEDCouplingMesh *mergeMyselfWith(const MEDCouplingMesh *other) const; DataArrayDouble *getCoordinatesAndOwner() const; DataArrayDouble *getBarycenterAndOwner() const; + void renumberCells(const int *old2NewBg, const int *old2NewEnd, bool check) throw(INTERP_KERNEL::Exception); //some useful methods void getSplitCellValues(int *res) const; void getSplitNodeValues(int *res) const; @@ -73,6 +79,7 @@ namespace ParaMEDMEM const std::vector& littleStrings); private: MEDCouplingCMesh(); + MEDCouplingCMesh(const MEDCouplingCMesh& other, bool deepCpy); ~MEDCouplingCMesh(); private: DataArrayDouble *_x_array; diff --git a/src/MEDCoupling/MEDCouplingExtrudedMesh.cxx b/src/MEDCoupling/MEDCouplingExtrudedMesh.cxx index bfa591c18..7b41a4173 100644 --- a/src/MEDCoupling/MEDCouplingExtrudedMesh.cxx +++ b/src/MEDCoupling/MEDCouplingExtrudedMesh.cxx @@ -61,6 +61,20 @@ MEDCouplingMeshType MEDCouplingExtrudedMesh::getType() const return EXTRUDED; } +/*! + * This method copyies all tiny strings from other (name and components name). + * @throw if other and this have not same mesh type. + */ +void MEDCouplingExtrudedMesh::copyTinyStringsFrom(const MEDCouplingMesh *other) throw(INTERP_KERNEL::Exception) +{ + const MEDCouplingExtrudedMesh *otherC=dynamic_cast(other); + if(!otherC) + throw INTERP_KERNEL::Exception("MEDCouplingExtrudedMesh::copyTinyStringsFrom : meshes have not same type !"); + MEDCouplingMesh::copyTinyStringsFrom(other); + _mesh2D->copyTinyStringsFrom(otherC->_mesh2D); + _mesh1D->copyTinyStringsFrom(otherC->_mesh1D); +} + MEDCouplingExtrudedMesh::MEDCouplingExtrudedMesh(const MEDCouplingUMesh *mesh3D, MEDCouplingUMesh *mesh2D, int cell2DId) throw(INTERP_KERNEL::Exception) try:_mesh2D(mesh2D),_mesh1D(MEDCouplingUMesh::New()),_mesh3D_ids(0),_cell_2D_id(cell2DId) { @@ -131,6 +145,11 @@ int MEDCouplingExtrudedMesh::getMeshDimension() const return 3; } +MEDCouplingMesh *MEDCouplingExtrudedMesh::deepCpy() const +{ + return clone(true); +} + MEDCouplingExtrudedMesh *MEDCouplingExtrudedMesh::clone(bool recDeepCpy) const { return new MEDCouplingExtrudedMesh(*this,recDeepCpy); @@ -154,6 +173,12 @@ bool MEDCouplingExtrudedMesh::isEqual(const MEDCouplingMesh *other, double prec) return true; } +void MEDCouplingExtrudedMesh::checkDeepEquivalWith(const MEDCouplingMesh *other, int cellCompPol, double prec, + DataArrayInt *&cellCor, DataArrayInt *&nodeCor) const throw(INTERP_KERNEL::Exception) +{ + throw INTERP_KERNEL::Exception("MEDCouplingExtrudedMesh::checkDeepEquivalWith : not implemented yet !"); +} + INTERP_KERNEL::NormalizedCellType MEDCouplingExtrudedMesh::getTypeOfCell(int cellId) const { const int *ids=_mesh3D_ids->getConstPointer(); @@ -249,6 +274,11 @@ void MEDCouplingExtrudedMesh::updateTime() } } +void MEDCouplingExtrudedMesh::renumberCells(const int *old2NewBg, const int *old2NewEnd, bool check) throw(INTERP_KERNEL::Exception) +{ + throw INTERP_KERNEL::Exception("Functionnality of renumbering cells unavailable for ExtrudedMesh"); +} + MEDCouplingUMesh *MEDCouplingExtrudedMesh::build3DUnstructuredMesh() const { MEDCouplingUMesh *ret=_mesh2D->buildExtrudedMeshFromThis(_mesh1D,0); diff --git a/src/MEDCoupling/MEDCouplingExtrudedMesh.hxx b/src/MEDCoupling/MEDCouplingExtrudedMesh.hxx index 566b23616..a92394fbe 100644 --- a/src/MEDCoupling/MEDCouplingExtrudedMesh.hxx +++ b/src/MEDCoupling/MEDCouplingExtrudedMesh.hxx @@ -38,13 +38,17 @@ namespace ParaMEDMEM static MEDCouplingExtrudedMesh *New(const MEDCouplingUMesh *mesh3D, MEDCouplingUMesh *mesh2D, int cell2DId) throw(INTERP_KERNEL::Exception); static MEDCouplingExtrudedMesh *New(); MEDCouplingMeshType getType() const; + void copyTinyStringsFrom(const MEDCouplingMesh *other) throw(INTERP_KERNEL::Exception); bool isStructured() const; int getNumberOfCells() const; int getNumberOfNodes() const; int getSpaceDimension() const; int getMeshDimension() const; + MEDCouplingMesh *deepCpy() const; MEDCouplingExtrudedMesh *clone(bool recDeepCpy) const; bool isEqual(const MEDCouplingMesh *other, double prec) const; + void checkDeepEquivalWith(const MEDCouplingMesh *other, int cellCompPol, double prec, + DataArrayInt *&cellCor, DataArrayInt *&nodeCor) const throw(INTERP_KERNEL::Exception); INTERP_KERNEL::NormalizedCellType getTypeOfCell(int cellId) const; int getNumberOfCellsWithType(INTERP_KERNEL::NormalizedCellType type) const; void getNodeIdsOfCell(int cellId, std::vector& conn) const; @@ -52,6 +56,7 @@ namespace ParaMEDMEM void checkCoherency() const throw (INTERP_KERNEL::Exception); void getBoundingBox(double *bbox) const; void updateTime(); + void renumberCells(const int *old2NewBg, const int *old2NewEnd, bool check) throw(INTERP_KERNEL::Exception); MEDCouplingUMesh *getMesh2D() const { return _mesh2D; } MEDCouplingUMesh *getMesh1D() const { return _mesh1D; } DataArrayInt *getMesh3DIds() const { return _mesh3D_ids; } diff --git a/src/MEDCoupling/MEDCouplingField.cxx b/src/MEDCoupling/MEDCouplingField.cxx index c8d23335e..b50f34fe5 100644 --- a/src/MEDCoupling/MEDCouplingField.cxx +++ b/src/MEDCoupling/MEDCouplingField.cxx @@ -40,19 +40,37 @@ bool MEDCouplingField::isEqual(const MEDCouplingField *other, double meshPrec, d return _mesh->isEqual(other->_mesh,meshPrec); } -bool MEDCouplingField::areCompatible(const MEDCouplingField *other) const +/*! + * This method states if 'this' and 'other' are compatibles each other before performing any treatment. + * This method is good for methods like : mergeFields. + * This method is not very demanding compared to areStrictlyCompatible that is better for operation on fields. + */ +bool MEDCouplingField::areCompatibleForMerge(const MEDCouplingField *other) const { if(!_type->isEqual(other->_type,1.)) return false; if(_mesh==other->_mesh) return true; - return _mesh->areCompatible(other->_mesh); + return _mesh->areCompatibleForMerge(other->_mesh); +} + +/*! + * This method is more strict than MEDCouplingField::areCompatible method. + * This method is used for operation on fields to operate a first check before attempting operation. + */ +bool MEDCouplingField::areStrictlyCompatible(const MEDCouplingField *other) const +{ + if(!_type->isEqual(other->_type,1.e-12)) + return false; + return _mesh==other->_mesh; } void MEDCouplingField::updateTime() { if(_mesh) updateTimeWith(*_mesh); + if(_type) + updateTimeWith(*_type); } TypeOfField MEDCouplingField::getTypeOfField() const diff --git a/src/MEDCoupling/MEDCouplingField.hxx b/src/MEDCoupling/MEDCouplingField.hxx index 4a8f72ecc..d7d45f29f 100644 --- a/src/MEDCoupling/MEDCouplingField.hxx +++ b/src/MEDCoupling/MEDCouplingField.hxx @@ -40,7 +40,8 @@ namespace ParaMEDMEM { public: virtual void checkCoherency() const throw(INTERP_KERNEL::Exception) = 0; - virtual bool areCompatible(const MEDCouplingField *other) const; + virtual bool areCompatibleForMerge(const MEDCouplingField *other) const; + virtual bool areStrictlyCompatible(const MEDCouplingField *other) const; virtual bool isEqual(const MEDCouplingField *other, double meshPrec, double valsPrec) const; void setMesh(const ParaMEDMEM::MEDCouplingMesh *mesh); const ParaMEDMEM::MEDCouplingMesh *getMesh() const { return _mesh; } diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx index 62b69dc0e..677f85c12 100644 --- a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx +++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx @@ -85,12 +85,90 @@ TypeOfField MEDCouplingFieldDiscretization::getTypeOfFieldFromStringRepr(const c throw INTERP_KERNEL::Exception("Representation does not match with any field discretization !"); } +/*! + * Excepted for MEDCouplingFieldDiscretizationPerCell no underlying TimeLabel object : nothing to do in generally. + */ +void MEDCouplingFieldDiscretization::updateTime() +{ +} + +/*! + * Computes normL1 of DataArrayDouble instance arr. + * @param res output parameter expected to be of size arr->getNumberOfComponents(); + * @throw when the field discretization fails on getWeighting fields (gauss points for example) + */ +void MEDCouplingFieldDiscretization::normL1(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, bool isWAbs, double *res) const throw(INTERP_KERNEL::Exception) +{ + MEDCouplingFieldDouble *vol=getWeightingField(mesh,isWAbs); + int nbOfCompo=arr->getNumberOfComponents(); + int nbOfElems=getNumberOfTuples(mesh); + std::fill(res,res+nbOfCompo,0.); + double totVol=0.; + const double *arrPtr=arr->getConstPointer(); + const double *volPtr=vol->getArray()->getConstPointer(); + for(int i=0;i(),1/totVol)); + vol->decrRef(); +} + +/*! + * Computes normL2 of DataArrayDouble instance arr. + * @param res output parameter expected to be of size arr->getNumberOfComponents(); + * @throw when the field discretization fails on getWeighting fields (gauss points for example) + */ +void MEDCouplingFieldDiscretization::normL2(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, bool isWAbs, double *res) const throw(INTERP_KERNEL::Exception) +{ + MEDCouplingFieldDouble *vol=getWeightingField(mesh,isWAbs); + int nbOfCompo=arr->getNumberOfComponents(); + int nbOfElems=getNumberOfTuples(mesh); + std::fill(res,res+nbOfCompo,0.); + double totVol=0.; + const double *arrPtr=arr->getConstPointer(); + const double *volPtr=vol->getArray()->getConstPointer(); + for(int i=0;i(),1/totVol)); + std::transform(res,res+nbOfCompo,res,std::ptr_fun(std::sqrt)); + vol->decrRef(); +} + +/*! + * Computes integral of DataArrayDouble instance arr. + * @param res output parameter expected to be of size arr->getNumberOfComponents(); + * @throw when the field discretization fails on getWeighting fields (gauss points for example) + */ +void MEDCouplingFieldDiscretization::integral(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, bool isWAbs, double *res) const throw(INTERP_KERNEL::Exception) +{ + MEDCouplingFieldDouble *vol=getWeightingField(mesh,isWAbs); + int nbOfCompo=arr->getNumberOfComponents(); + int nbOfElems=getNumberOfTuples(mesh); + std::fill(res,res+nbOfCompo,0.); + const double *arrPtr=arr->getConstPointer(); + const double *volPtr=vol->getArray()->getConstPointer(); + double *tmp=new double[nbOfCompo]; + for (int i=0;i(),volPtr[i])); + std::transform(tmp,tmp+nbOfCompo,res,res,std::plus()); + } + delete [] tmp; + vol->decrRef(); +} + void MEDCouplingFieldDiscretization::getSerializationIntArray(DataArrayInt *& arr) const { arr=0; } - /*! * Empty : Not a bug */ @@ -117,6 +195,14 @@ void MEDCouplingFieldDiscretization::finishUnserialization(const std::vectorgetNumberOfCells(); } +void MEDCouplingFieldDiscretizationP0::renumberArraysForCell(const MEDCouplingMesh *, const std::vector& arrays, + const int *old2NewBg, const int *old2NewEnd, bool check) throw(INTERP_KERNEL::Exception) +{ + const int *array=old2NewBg; + if(check) + array=DataArrayInt::checkAndPreparePermutation(old2NewBg,old2NewEnd); + for(std::vector::const_iterator it=arrays.begin();it!=arrays.end();it++) + { + if(*it) + (*it)->renumberInPlace(array); + } + if(check) + delete [] array; +} + DataArrayDouble *MEDCouplingFieldDiscretizationP0::getLocalizationOfDiscValues(const MEDCouplingMesh *mesh) const { return mesh->getBarycenterAndOwner(); @@ -244,7 +345,7 @@ void MEDCouplingFieldDiscretizationP0::getValueOnPos(const DataArrayDouble *arr, /*! * Nothing to do. It's not a bug. */ -void MEDCouplingFieldDiscretizationP0::renumberValuesOnNodes(const DataArrayInt *old2New, DataArrayDouble *arr) const +void MEDCouplingFieldDiscretizationP0::renumberValuesOnNodes(const int *, DataArrayDouble *) const { } @@ -284,6 +385,14 @@ bool MEDCouplingFieldDiscretizationP1::isEqual(const MEDCouplingFieldDiscretizat return otherC!=0; } +/*! + * Nothing to do here. + */ +void MEDCouplingFieldDiscretizationP1::renumberArraysForCell(const MEDCouplingMesh *, const std::vector& arrays, + const int *old2NewBg, const int *old2NewEnd, bool check) throw(INTERP_KERNEL::Exception) +{ +} + int MEDCouplingFieldDiscretizationP1::getNumberOfTuples(const MEDCouplingMesh *mesh) const { return mesh->getNumberOfNodes(); @@ -358,10 +467,9 @@ void MEDCouplingFieldDiscretizationP1::getValueOnPos(const DataArrayDouble *arr, arr->getTuple(id,res); } -void MEDCouplingFieldDiscretizationP1::renumberValuesOnNodes(const DataArrayInt *old2New, DataArrayDouble *arr) const +void MEDCouplingFieldDiscretizationP1::renumberValuesOnNodes(const int *old2NewPtr, DataArrayDouble *arr) const { - int oldNbOfElems=old2New->getNbOfElems(); - const int *old2NewPtr=old2New->getConstPointer(); + int oldNbOfElems=arr->getNumberOfTuples(); int nbOfComp=arr->getNumberOfComponents(); int newNbOfTuples=(*std::max_element(old2NewPtr,old2NewPtr+oldNbOfElems))+1; DataArrayDouble *arrCpy=arr->deepCopy(); @@ -437,6 +545,12 @@ MEDCouplingFieldDiscretizationPerCell::MEDCouplingFieldDiscretizationPerCell(con _discr_per_cell=arr->deepCopy(); } +void MEDCouplingFieldDiscretizationPerCell::updateTime() +{ + if(_discr_per_cell) + updateTimeWith(*_discr_per_cell); +} + void MEDCouplingFieldDiscretizationPerCell::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception) { if(!_discr_per_cell) @@ -458,6 +572,27 @@ bool MEDCouplingFieldDiscretizationPerCell::isEqual(const MEDCouplingFieldDiscre return _discr_per_cell->isEqual(*otherC->_discr_per_cell); } +/*! + * This method is typically the first step of renumbering. The impact on _discr_per_cell is necessary here. + * virtualy by this method. + */ +void MEDCouplingFieldDiscretizationPerCell::renumberCells(const int *old2NewBg, const int *old2NewEnd, bool check) throw(INTERP_KERNEL::Exception) +{ + int nbCells=_discr_per_cell->getNumberOfTuples(); + const int *array=old2NewBg; + if(std::distance(old2NewBg,old2NewEnd)!=nbCells) + throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationPerCell::renumberCells unexpected size of array !"); + if(check) + array=DataArrayInt::checkAndPreparePermutation(old2NewBg,old2NewEnd); + // + DataArrayInt *dpc=_discr_per_cell->renumber(array); + _discr_per_cell->decrRef(); + _discr_per_cell=dpc; + // + if(check) + delete [] (int *)array; +} + void MEDCouplingFieldDiscretizationPerCell::buildDiscrPerCellIfNecessary(const MEDCouplingMesh *m) { if(!_discr_per_cell) @@ -509,7 +644,7 @@ const char *MEDCouplingFieldDiscretizationGauss::getStringRepr() const return REPR; } -int MEDCouplingFieldDiscretizationGauss::getNumberOfTuples(const MEDCouplingMesh *mesh) const +int MEDCouplingFieldDiscretizationGauss::getNumberOfTuples(const MEDCouplingMesh *) const { int ret=0; const int *dcPtr=_discr_per_cell->getConstPointer(); @@ -519,6 +654,36 @@ int MEDCouplingFieldDiscretizationGauss::getNumberOfTuples(const MEDCouplingMesh return ret; } +void MEDCouplingFieldDiscretizationGauss::renumberArraysForCell(const MEDCouplingMesh *, const std::vector& arrays, + const int *old2NewBg, const int *old2NewEnd, bool check) throw(INTERP_KERNEL::Exception) +{ + const int *array=old2NewBg; + if(check) + array=DataArrayInt::checkAndPreparePermutation(old2NewBg,old2NewEnd); + int nbOfCells=_discr_per_cell->getNumberOfTuples(); + int nbOfTuples=getNumberOfTuples(0); + const int *dcPtr=_discr_per_cell->getConstPointer(); + int *array2=new int[nbOfTuples];//stores the final conversion array old2New to give to arrays in renumberInPlace. + int *array3=new int[nbOfCells];//store for each cell in present dcp array (already renumbered) the offset needed by each cell in new numbering. + array3[0]=0; + for(int i=1;i::const_iterator it=arrays.begin();it!=arrays.end();it++) + if(*it) + (*it)->renumberInPlace(array2); + delete [] array2; + if(check) + delete [] (int*)array; +} + DataArrayDouble *MEDCouplingFieldDiscretizationGauss::getLocalizationOfDiscValues(const MEDCouplingMesh *mesh) const { throw INTERP_KERNEL::Exception("Not implemented yet !"); @@ -655,9 +820,11 @@ MEDCouplingMesh *MEDCouplingFieldDiscretizationGauss::buildSubMeshData(const int throw INTERP_KERNEL::Exception("Not implemented yet !"); } -void MEDCouplingFieldDiscretizationGauss::renumberValuesOnNodes(const DataArrayInt *old2New, DataArrayDouble *arr) const +/*! + * No implementation needed ! + */ +void MEDCouplingFieldDiscretizationGauss::renumberValuesOnNodes(const int *, DataArrayDouble *) const { - throw INTERP_KERNEL::Exception("Not implemented yet !"); } void MEDCouplingFieldDiscretizationGauss::setGaussLocalizationOnType(const MEDCouplingMesh *m, INTERP_KERNEL::NormalizedCellType type, const std::vector& refCoo, @@ -855,6 +1022,40 @@ int MEDCouplingFieldDiscretizationGaussNE::getNumberOfTuples(const MEDCouplingMe return ret; } +void MEDCouplingFieldDiscretizationGaussNE::renumberArraysForCell(const MEDCouplingMesh *mesh, const std::vector& arrays, + const int *old2NewBg, const int *old2NewEnd, bool check) throw(INTERP_KERNEL::Exception) +{ + const int *array=old2NewBg; + if(check) + array=DataArrayInt::checkAndPreparePermutation(old2NewBg,old2NewEnd); + int nbOfCells=mesh->getNumberOfCells(); + int nbOfTuples=getNumberOfTuples(mesh); + int *array2=new int[nbOfTuples];//stores the final conversion array old2New to give to arrays in renumberInPlace. + int *array3=new int[nbOfCells];//store for each cell in after renumbering the offset needed by each cell in new numbering. + array3[0]=0; + for(int i=1;igetTypeOfCell(std::distance(array,std::find(array,array+nbOfCells,i-1))); + const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::getCellModel(type); + array3[i]=array3[i-1]+cm.getNumberOfNodes(); + } + int j=0; + for(int i=0;igetTypeOfCell(i); + const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::getCellModel(type); + for(int k=0;k::const_iterator it=arrays.begin();it!=arrays.end();it++) + if(*it) + (*it)->renumberInPlace(array2); + delete [] array2; + if(check) + delete [] (int*)array; +} + DataArrayDouble *MEDCouplingFieldDiscretizationGaussNE::getLocalizationOfDiscValues(const MEDCouplingMesh *mesh) const { throw INTERP_KERNEL::Exception("Not implemented yet !"); @@ -907,9 +1108,11 @@ MEDCouplingMesh *MEDCouplingFieldDiscretizationGaussNE::buildSubMeshData(const i throw INTERP_KERNEL::Exception("Not implemented yet !"); } -void MEDCouplingFieldDiscretizationGaussNE::renumberValuesOnNodes(const DataArrayInt *old2New, DataArrayDouble *arr) const +/*! + * No implementation needed ! + */ +void MEDCouplingFieldDiscretizationGaussNE::renumberValuesOnNodes(const int *, DataArrayDouble *) const { - throw INTERP_KERNEL::Exception("Not implemented yet !"); } MEDCouplingFieldDiscretizationGaussNE::MEDCouplingFieldDiscretizationGaussNE(const MEDCouplingFieldDiscretizationGaussNE& other):MEDCouplingFieldDiscretization(other) diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx index ae98cef79..ade9f999c 100644 --- a/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx +++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx @@ -17,15 +17,18 @@ // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // -#ifndef __MEDCOUPLINGFIELDDISCRETIZATION_HXX__ -#define __MEDCOUPLINGFIELDDISCRETIZATION_HXX__ +#ifndef __PARAMEDMEM_MEDCOUPLINGFIELDDISCRETIZATION_HXX__ +#define __PARAMEDMEM_MEDCOUPLINGFIELDDISCRETIZATION_HXX__ #include "MEDCoupling.hxx" #include "MEDCouplingRefCountObject.hxx" #include "InterpKernelException.hxx" +#include "MEDCouplingTimeLabel.hxx" #include "MEDCouplingNatureOfField.hxx" #include "MEDCouplingGaussLocalization.hxx" +#include + namespace ParaMEDMEM { class DataArrayInt; @@ -33,27 +36,34 @@ namespace ParaMEDMEM class DataArrayDouble; class MEDCouplingFieldDouble; - class MEDCOUPLING_EXPORT MEDCouplingFieldDiscretization + class MEDCOUPLING_EXPORT MEDCouplingFieldDiscretization : public TimeLabel { public: static MEDCouplingFieldDiscretization *New(TypeOfField type); double getPrecision() const { return _precision; } void setPrecision(double val) { _precision=val; } + void updateTime(); static TypeOfField getTypeOfFieldFromStringRepr(const char *repr) throw(INTERP_KERNEL::Exception); virtual TypeOfField getEnum() const = 0; virtual bool isEqual(const MEDCouplingFieldDiscretization *other, double eps) const = 0; virtual MEDCouplingFieldDiscretization *clone() const = 0; virtual const char *getStringRepr() const = 0; virtual int getNumberOfTuples(const MEDCouplingMesh *mesh) const = 0; + virtual void normL1(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, bool isWAbs, double *res) const throw(INTERP_KERNEL::Exception); + virtual void normL2(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, bool isWAbs, double *res) const throw(INTERP_KERNEL::Exception); + virtual void integral(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, bool isWAbs, double *res) const throw(INTERP_KERNEL::Exception); virtual DataArrayDouble *getLocalizationOfDiscValues(const MEDCouplingMesh *mesh) const = 0; virtual void checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception) = 0; + virtual void renumberCells(const int *old2NewBg, const int *old2NewEnd, bool check) throw(INTERP_KERNEL::Exception); + virtual void renumberArraysForCell(const MEDCouplingMesh *mesh, const std::vector& arrays, + const int *old2NewBg, const int *old2NewEnd, bool check) throw(INTERP_KERNEL::Exception) = 0; virtual double getIJK(const MEDCouplingMesh *mesh, const DataArrayDouble *da, int cellId, int nodeIdInCell, int compoId) const throw(INTERP_KERNEL::Exception); virtual void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception) = 0; virtual MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh, bool isAbs) const = 0; virtual void getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const = 0; virtual void getValueOnPos(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, int i, int j, int k, double *res) const = 0; virtual MEDCouplingMesh *buildSubMeshData(const int *start, const int *end, const MEDCouplingMesh *mesh, DataArrayInt *&di) const = 0; - virtual void renumberValuesOnNodes(const DataArrayInt *old2New, DataArrayDouble *arr) const = 0; + virtual void renumberValuesOnNodes(const int *old2New, DataArrayDouble *arr) const = 0; virtual void getSerializationIntArray(DataArrayInt *& arr) const; virtual void getTinySerializationIntInformation(std::vector& tinyInfo) const; virtual void getTinySerializationDbleInformation(std::vector& tinyInfo) const; @@ -86,13 +96,15 @@ namespace ParaMEDMEM const char *getStringRepr() const; bool isEqual(const MEDCouplingFieldDiscretization *other, double eps) const; int getNumberOfTuples(const MEDCouplingMesh *mesh) const; + void renumberArraysForCell(const MEDCouplingMesh *mesh, const std::vector& arrays, + const int *old2NewBg, const int *old2NewEnd, bool check) throw(INTERP_KERNEL::Exception); DataArrayDouble *getLocalizationOfDiscValues(const MEDCouplingMesh *mesh) const; void checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception); void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception); MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh, bool isAbs) const; void getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const; void getValueOnPos(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, int i, int j, int k, double *res) const; - void renumberValuesOnNodes(const DataArrayInt *old2New, DataArrayDouble *arr) const; + void renumberValuesOnNodes(const int *old2New, DataArrayDouble *arr) const; MEDCouplingMesh *buildSubMeshData(const int *start, const int *end, const MEDCouplingMesh *mesh, DataArrayInt *&di) const; public: static const char REPR[]; @@ -107,6 +119,8 @@ namespace ParaMEDMEM const char *getStringRepr() const; bool isEqual(const MEDCouplingFieldDiscretization *other, double eps) const; int getNumberOfTuples(const MEDCouplingMesh *mesh) const; + void renumberArraysForCell(const MEDCouplingMesh *mesh, const std::vector& arrays, + const int *old2NewBg, const int *old2NewEnd, bool check) throw(INTERP_KERNEL::Exception); DataArrayDouble *getLocalizationOfDiscValues(const MEDCouplingMesh *mesh) const; void checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception); void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception); @@ -114,7 +128,7 @@ namespace ParaMEDMEM void getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const; void getValueOnPos(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, int i, int j, int k, double *res) const; MEDCouplingMesh *buildSubMeshData(const int *start, const int *end, const MEDCouplingMesh *mesh, DataArrayInt *&di) const; - void renumberValuesOnNodes(const DataArrayInt *old2New, DataArrayDouble *arr) const; + void renumberValuesOnNodes(const int *old2New, DataArrayDouble *arr) const; static DataArrayInt *invertArrayO2N2N2O(const MEDCouplingMesh *mesh, const DataArrayInt *di); public: static const char REPR[]; @@ -131,8 +145,10 @@ namespace ParaMEDMEM MEDCouplingFieldDiscretizationPerCell(); MEDCouplingFieldDiscretizationPerCell(const MEDCouplingFieldDiscretizationPerCell& other); ~MEDCouplingFieldDiscretizationPerCell(); + void updateTime(); void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception); bool isEqual(const MEDCouplingFieldDiscretization *other, double eps) const; + void renumberCells(const int *old2NewBg, const int *old2NewEnd, bool check) throw(INTERP_KERNEL::Exception); protected: void buildDiscrPerCellIfNecessary(const MEDCouplingMesh *m); protected: @@ -148,6 +164,8 @@ namespace ParaMEDMEM MEDCouplingFieldDiscretization *clone() const; const char *getStringRepr() const; int getNumberOfTuples(const MEDCouplingMesh *mesh) const; + void renumberArraysForCell(const MEDCouplingMesh *mesh, const std::vector& arrays, + const int *old2NewBg, const int *old2NewEnd, bool check) throw(INTERP_KERNEL::Exception); DataArrayDouble *getLocalizationOfDiscValues(const MEDCouplingMesh *mesh) const; void checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception); void getTinySerializationIntInformation(std::vector& tinyInfo) const; @@ -161,7 +179,7 @@ namespace ParaMEDMEM void getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const; void getValueOnPos(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, int i, int j, int k, double *res) const; MEDCouplingMesh *buildSubMeshData(const int *start, const int *end, const MEDCouplingMesh *mesh, DataArrayInt *&di) const; - void renumberValuesOnNodes(const DataArrayInt *old2New, DataArrayDouble *arr) const; + void renumberValuesOnNodes(const int *old2New, DataArrayDouble *arr) const; void setGaussLocalizationOnType(const MEDCouplingMesh *m, INTERP_KERNEL::NormalizedCellType type, const std::vector& refCoo, const std::vector& gsCoo, const std::vector& wg) throw(INTERP_KERNEL::Exception); void setGaussLocalizationOnCells(const MEDCouplingMesh *m, const int *begin, const int *end, const std::vector& refCoo, @@ -197,6 +215,8 @@ namespace ParaMEDMEM const char *getStringRepr() const; bool isEqual(const MEDCouplingFieldDiscretization *other, double eps) const; int getNumberOfTuples(const MEDCouplingMesh *mesh) const; + void renumberArraysForCell(const MEDCouplingMesh *mesh, const std::vector& arrays, + const int *old2NewBg, const int *old2NewEnd, bool check) throw(INTERP_KERNEL::Exception); DataArrayDouble *getLocalizationOfDiscValues(const MEDCouplingMesh *mesh) const; void checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception); double getIJK(const MEDCouplingMesh *mesh, const DataArrayDouble *da, int cellId, int nodeIdInCell, int compoId) const throw(INTERP_KERNEL::Exception); @@ -205,7 +225,7 @@ namespace ParaMEDMEM void getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const; void getValueOnPos(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, int i, int j, int k, double *res) const; MEDCouplingMesh *buildSubMeshData(const int *start, const int *end, const MEDCouplingMesh *mesh, DataArrayInt *&di) const; - void renumberValuesOnNodes(const DataArrayInt *old2New, DataArrayDouble *arr) const; + void renumberValuesOnNodes(const int *old2New, DataArrayDouble *arr) const; protected: MEDCouplingFieldDiscretizationGaussNE(const MEDCouplingFieldDiscretizationGaussNE& other); public: diff --git a/src/MEDCoupling/MEDCouplingFieldDouble.cxx b/src/MEDCoupling/MEDCouplingFieldDouble.cxx index 7a2914fcd..846866906 100644 --- a/src/MEDCoupling/MEDCouplingFieldDouble.cxx +++ b/src/MEDCoupling/MEDCouplingFieldDouble.cxx @@ -21,6 +21,7 @@ #include "MEDCouplingPointSet.hxx" #include "MEDCouplingTimeDiscretization.hxx" #include "MEDCouplingFieldDiscretization.hxx" +#include "MEDCouplingAutoRefCountObjectPtr.hxx" #include #include @@ -61,9 +62,14 @@ bool MEDCouplingFieldDouble::isEqual(const MEDCouplingField *other, double meshP return true; } -bool MEDCouplingFieldDouble::areCompatible(const MEDCouplingField *other) const +/*! + * This method states if 'this' and 'other' are compatibles each other before performing any treatment. + * This method is good for methods like : mergeFields. + * This method is not very demanding compared to areStrictlyCompatible that is better for operation on fields. + */ +bool MEDCouplingFieldDouble::areCompatibleForMerge(const MEDCouplingField *other) const { - if(!MEDCouplingField::areCompatible(other)) + if(!MEDCouplingField::areCompatibleForMerge(other)) return false; const MEDCouplingFieldDouble *otherC=dynamic_cast(other); if(!otherC) @@ -75,20 +81,76 @@ bool MEDCouplingFieldDouble::areCompatible(const MEDCouplingField *other) const return true; } +/*! + * This method is more strict than MEDCouplingField::areCompatible method. + * This method is used for operation on fields to operate a first check before attempting operation. + */ +bool MEDCouplingFieldDouble::areStrictlyCompatible(const MEDCouplingField *other) const +{ + if(!MEDCouplingField::areStrictlyCompatible(other)) + return false; + const MEDCouplingFieldDouble *otherC=dynamic_cast(other); + if(!otherC) + return false; + if(_nature!=otherC->_nature) + return false; + if(!_time_discr->areStrictlyCompatible(otherC->_time_discr)) + return false; + return true; +} + bool MEDCouplingFieldDouble::areCompatibleForMul(const MEDCouplingField *other) const { - if(!MEDCouplingField::areCompatible(other)) + if(!MEDCouplingField::areStrictlyCompatible(other)) return false; const MEDCouplingFieldDouble *otherC=dynamic_cast(other); if(!otherC) return false; if(_nature!=otherC->_nature) return false; - if(!_time_discr->areCompatibleForMul(otherC->_time_discr)) + if(!_time_discr->areStrictlyCompatibleForMul(otherC->_time_discr)) return false; return true; } +/*! + * This method performs a clone of mesh and a renumbering of underlying cells of it. The number of cells remains the same. + * The values of field are impacted in consequence to have the same geometrical field. + */ +void MEDCouplingFieldDouble::renumberCells(const int *old2NewBg, const int *old2NewEnd, bool check) throw(INTERP_KERNEL::Exception) +{ + if(!_mesh) + throw INTERP_KERNEL::Exception("Expecting a defined mesh to be able to operate a renumbering !"); + MEDCouplingAutoRefCountObjectPtr m=_mesh->deepCpy(); + m->renumberCells(old2NewBg,old2NewEnd,check); + // + _type->renumberCells(old2NewBg,old2NewEnd,check); + std::vector arrays; + _time_discr->getArrays(arrays); + _type->renumberArraysForCell(_mesh,arrays,old2NewBg,old2NewEnd,check); + // + setMesh(m); + updateTime(); +} + +/*! + * This method performs a clone of mesh and a renumbering of underlying nodes of it. The number of nodes remains not compulsory the same as renumberCells method. + * The values of field are impacted in consequence to have the same geometrical field. + */ +void MEDCouplingFieldDouble::renumberNodes(const int *old2NewBg, const int *old2NewEnd) throw(INTERP_KERNEL::Exception) +{ + const MEDCouplingPointSet *meshC=dynamic_cast(_mesh); + if(!meshC) + throw INTERP_KERNEL::Exception("Invalid mesh to apply renumberNodes on it !"); + MEDCouplingAutoRefCountObjectPtr meshC2((MEDCouplingPointSet *)meshC->deepCpy()); + std::vector arrays; + _time_discr->getArrays(arrays); + for(std::vector::const_iterator iter=arrays.begin();iter!=arrays.end();iter++) + _type->renumberValuesOnNodes(old2NewBg,*iter); + meshC2->renumberNodes(old2NewBg,*std::max_element(old2NewBg,old2NewEnd)+1); + setMesh(meshC2); +} + TypeOfTimeDiscretization MEDCouplingFieldDouble::getTimeDiscretization() const { return _time_discr->getEnum(); @@ -130,6 +192,8 @@ double MEDCouplingFieldDouble::accumulate(int compId) const const double *ptr=getArray()->getConstPointer(); int nbTuple=getArray()->getNumberOfTuples(); int nbComps=getArray()->getNumberOfComponents(); + if(compId>=nbComps) + throw INTERP_KERNEL::Exception("Invalid compId specified : No such nb of components !"); double ret=0.; for(int i=0;i()); } +/*! + * Returns the normL1 of current field on compId component : + * \f[ + * \frac{\sum_{0 \leq i < nbOfEntity}|val[i]*Vol[i]|}{\sum_{0 \leq i < nbOfEntity}|Vol[i]|} + * \f] + * If compId>=nbOfComponent an exception is thrown. + */ +double MEDCouplingFieldDouble::normL1(int compId, bool isWAbs) const throw(INTERP_KERNEL::Exception) +{ + if(!_mesh) + throw INTERP_KERNEL::Exception("No mesh underlying this field to perform normL1"); + int nbComps=getArray()->getNumberOfComponents(); + if(compId>=nbComps) + throw INTERP_KERNEL::Exception("Invalid compId specified : No such nb of components !"); + double *res=new double[nbComps]; + try + { + _type->normL1(_mesh,getArray(),isWAbs,res); + } + catch(INTERP_KERNEL::Exception& e) + { + delete [] res; + throw e; + } + double ret=res[compId]; + delete [] res; + return ret; +} + +/*! + * Returns the normL1 of current field on each components : + * \f[ + * \frac{\sum_{0 \leq i < nbOfEntity}|val[i]*Vol[i]|}{\sum_{0 \leq i < nbOfEntity}|Vol[i]|} + * \f] + * The res is expected to be of size getNumberOfComponents(). + */ +void MEDCouplingFieldDouble::normL1(bool isWAbs, double *res) const throw(INTERP_KERNEL::Exception) +{ + if(!_mesh) + throw INTERP_KERNEL::Exception("No mesh underlying this field to perform normL1"); + _type->normL1(_mesh,getArray(),isWAbs,res); +} + +/*! + * Returns the normL2 of current field on compId component : + * \f[ + * \sqrt{\frac{\sum_{0 \leq i < nbOfEntity}|val[i]^{2}*Vol[i]|}{\sum_{0 \leq i < nbOfEntity}|Vol[i]|}} + * \f] + * If compId>=nbOfComponent an exception is thrown. + */ +double MEDCouplingFieldDouble::normL2(int compId, bool isWAbs) const throw(INTERP_KERNEL::Exception) +{ + if(!_mesh) + throw INTERP_KERNEL::Exception("No mesh underlying this field to perform normL1"); + int nbComps=getArray()->getNumberOfComponents(); + if(compId>=nbComps) + throw INTERP_KERNEL::Exception("Invalid compId specified : No such nb of components !"); + double *res=new double[nbComps]; + try + { + _type->normL2(_mesh,getArray(),isWAbs,res); + } + catch(INTERP_KERNEL::Exception& e) + { + delete [] res; + throw e; + } + double ret=res[compId]; + delete [] res; + return ret; +} + +/*! + * Returns the normL2 of current field on each components : + * \f[ + * \sqrt{\frac{\sum_{0 \leq i < nbOfEntity}|val[i]^{2}*Vol[i]|}{\sum_{0 \leq i < nbOfEntity}|Vol[i]|}} + * \f] + * The res is expected to be of size getNumberOfComponents(). + */ +void MEDCouplingFieldDouble::normL2(bool isWAbs, double *res) const throw(INTERP_KERNEL::Exception) +{ + if(!_mesh) + throw INTERP_KERNEL::Exception("No mesh underlying this field to perform normL1"); + _type->normL2(_mesh,getArray(),isWAbs,res); +} + /*! * Returns the accumulation (the sum) of comId_th component of each tuples weigthed by the field * returns by getWeightingField relative of the _type of field of default array. * This method is usefull to check the conservativity of interpolation method. */ -double MEDCouplingFieldDouble::measureAccumulate(int compId, bool isWAbs) const +double MEDCouplingFieldDouble::integral(int compId, bool isWAbs) const throw(INTERP_KERNEL::Exception) { if(!_mesh) - throw INTERP_KERNEL::Exception("No mesh underlying this field to perform measureAccumulate"); - MEDCouplingFieldDouble *weight=_type->getWeightingField(_mesh,isWAbs); - const double *ptr=weight->getArray()->getConstPointer(); - int nbOfValues=weight->getArray()->getNbOfElems(); - double ret=0.; - for (int i=0; idecrRef(); + throw INTERP_KERNEL::Exception("No mesh underlying this field to perform integral"); + int nbComps=getArray()->getNumberOfComponents(); + if(compId>=nbComps) + throw INTERP_KERNEL::Exception("Invalid compId specified : No such nb of components !"); + double *res=new double[nbComps]; + try + { + _type->integral(_mesh,getArray(),isWAbs,res); + } + catch(INTERP_KERNEL::Exception& e) + { + delete [] res; + throw e; + } + double ret=res[compId]; + delete [] res; return ret; } @@ -174,24 +332,11 @@ double MEDCouplingFieldDouble::measureAccumulate(int compId, bool isWAbs) const * returns by getWeightingField relative of the _type of field of default array. * This method is usefull to check the conservativity of interpolation method. */ -void MEDCouplingFieldDouble::measureAccumulate(bool isWAbs, double *res) const +void MEDCouplingFieldDouble::integral(bool isWAbs, double *res) const throw(INTERP_KERNEL::Exception) { if(!_mesh) - throw INTERP_KERNEL::Exception("No mesh underlying this field to perform measureAccumulate2"); - MEDCouplingFieldDouble *weight=_type->getWeightingField(_mesh,isWAbs); - const double *ptr=weight->getArray()->getConstPointer(); - int nbOfValues=weight->getArray()->getNbOfElems(); - int nbComps=getArray()->getNumberOfComponents(); - const double *vals=getArray()->getConstPointer(); - std::fill(res,res+nbComps,0.); - double *tmp=new double[nbComps]; - for (int i=0; i(),ptr[i])); - std::transform(tmp,tmp+nbComps,res,res,std::plus()); - } - weight->decrRef(); - delete [] tmp; + throw INTERP_KERNEL::Exception("No mesh underlying this field to perform integral2"); + _type->integral(_mesh,getArray(),isWAbs,res); } /*! @@ -398,36 +543,29 @@ void MEDCouplingFieldDouble::serialize(DataArrayInt *&dataInt, std::vector((MEDCouplingMesh *)(_mesh)); + const MEDCouplingPointSet *meshC=dynamic_cast(_mesh); if(!meshC) throw INTERP_KERNEL::Exception("Invalid mesh to apply mergeNodes on it !"); + MEDCouplingAutoRefCountObjectPtr meshC2((MEDCouplingPointSet *)meshC->deepCpy()); bool ret; - DataArrayInt *arr=meshC->mergeNodes(eps,ret); + MEDCouplingAutoRefCountObjectPtr arr=meshC2->mergeNodes(eps,ret); if(!ret)//no nodes have been merged. return ret; std::vector arrays; _time_discr->getArrays(arrays); - try - { - for(std::vector::const_iterator iter=arrays.begin();iter!=arrays.end();iter++) - _type->renumberValuesOnNodes(arr,*iter); - } - catch(INTERP_KERNEL::Exception& e) - { - arr->decrRef(); - throw e; - } - arr->decrRef(); + for(std::vector::const_iterator iter=arrays.begin();iter!=arrays.end();iter++) + _type->renumberValuesOnNodes(arr->getConstPointer(),*iter); + setMesh(meshC2); return true; } MEDCouplingFieldDouble *MEDCouplingFieldDouble::mergeFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) { - if(!f1->areCompatible(f2)) + if(!f1->areCompatibleForMerge(f2)) throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply mergeFields on them !"); const MEDCouplingMesh *m1=f1->getMesh(); const MEDCouplingMesh *m2=f2->getMesh(); @@ -443,7 +581,7 @@ MEDCouplingFieldDouble *MEDCouplingFieldDouble::mergeFields(const MEDCouplingFie MEDCouplingFieldDouble *MEDCouplingFieldDouble::addFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) { - if(!f1->areCompatible(f2)) + if(!f1->areStrictlyCompatible(f2)) throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply addFields on them !"); MEDCouplingTimeDiscretization *td=f1->_time_discr->add(f2->_time_discr); td->copyTinyAttrFrom(*f1->_time_discr); @@ -454,7 +592,7 @@ MEDCouplingFieldDouble *MEDCouplingFieldDouble::addFields(const MEDCouplingField const MEDCouplingFieldDouble &MEDCouplingFieldDouble::operator+=(const MEDCouplingFieldDouble& other) { - if(!areCompatible(&other)) + if(!areStrictlyCompatible(&other)) throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply += on them !"); _time_discr->addEqual(other._time_discr); return *this; @@ -462,7 +600,7 @@ const MEDCouplingFieldDouble &MEDCouplingFieldDouble::operator+=(const MEDCoupli MEDCouplingFieldDouble *MEDCouplingFieldDouble::substractFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) { - if(!f1->areCompatible(f2)) + if(!f1->areStrictlyCompatible(f2)) throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply substractFields on them !"); MEDCouplingTimeDiscretization *td=f1->_time_discr->substract(f2->_time_discr); td->copyTinyAttrFrom(*f1->_time_discr); @@ -473,7 +611,7 @@ MEDCouplingFieldDouble *MEDCouplingFieldDouble::substractFields(const MEDCouplin const MEDCouplingFieldDouble &MEDCouplingFieldDouble::operator-=(const MEDCouplingFieldDouble& other) { - if(!areCompatible(&other)) + if(!areStrictlyCompatible(&other)) throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply -= on them !"); _time_discr->substractEqual(other._time_discr); return *this; @@ -500,7 +638,7 @@ const MEDCouplingFieldDouble &MEDCouplingFieldDouble::operator*=(const MEDCoupli MEDCouplingFieldDouble *MEDCouplingFieldDouble::divideFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) { - if(!f1->areCompatible(f2)) + if(!f1->areStrictlyCompatible(f2)) throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply divideFields on them !"); MEDCouplingTimeDiscretization *td=f1->_time_discr->divide(f2->_time_discr); td->copyTinyAttrFrom(*f1->_time_discr); @@ -511,7 +649,7 @@ MEDCouplingFieldDouble *MEDCouplingFieldDouble::divideFields(const MEDCouplingFi const MEDCouplingFieldDouble &MEDCouplingFieldDouble::operator/=(const MEDCouplingFieldDouble& other) { - if(!areCompatible(&other)) + if(!areStrictlyCompatible(&other)) throw INTERP_KERNEL::Exception("Fields are not compatible ; unable to apply /= on them !"); _time_discr->divideEqual(other._time_discr); return *this; diff --git a/src/MEDCoupling/MEDCouplingFieldDouble.hxx b/src/MEDCoupling/MEDCouplingFieldDouble.hxx index 336885e0e..8e5bf0137 100644 --- a/src/MEDCoupling/MEDCouplingFieldDouble.hxx +++ b/src/MEDCoupling/MEDCouplingFieldDouble.hxx @@ -33,8 +33,11 @@ namespace ParaMEDMEM public: static MEDCouplingFieldDouble *New(TypeOfField type, TypeOfTimeDiscretization td=NO_TIME); bool isEqual(const MEDCouplingField *other, double meshPrec, double valsPrec) const; - bool areCompatible(const MEDCouplingField *other) const; + bool areCompatibleForMerge(const MEDCouplingField *other) const; + bool areStrictlyCompatible(const MEDCouplingField *other) const; bool areCompatibleForMul(const MEDCouplingField *other) const; + void renumberCells(const int *old2NewBg, const int *old2NewEnd, bool check) throw(INTERP_KERNEL::Exception); + void renumberNodes(const int *old2NewBg, const int *old2NewEnd) throw(INTERP_KERNEL::Exception); MEDCouplingFieldDouble *clone(bool recDeepCpy) const; MEDCouplingFieldDouble *buildNewTimeReprFromThis(TypeOfTimeDiscretization td, bool deepCpy) const; TypeOfTimeDiscretization getTimeDiscretization() const; @@ -55,8 +58,12 @@ namespace ParaMEDMEM DataArrayDouble *getEndArray() const { return _time_discr->getEndArray(); } double accumulate(int compId) const; void accumulate(double *res) const; - double measureAccumulate(int compId, bool isWAbs) const; - void measureAccumulate(bool isWAbs, double *res) const; + double normL1(int compId, bool isWAbs) const throw(INTERP_KERNEL::Exception); + void normL1(bool isWAbs, double *res) const throw(INTERP_KERNEL::Exception); + double normL2(int compId, bool isWAbs) const throw(INTERP_KERNEL::Exception); + void normL2(bool isWAbs, double *res) const throw(INTERP_KERNEL::Exception); + double integral(int compId, bool isWAbs) const throw(INTERP_KERNEL::Exception); + void integral(bool isWAbs, double *res) const throw(INTERP_KERNEL::Exception); void getValueOnPos(int i, int j, int k, double *res) const throw(INTERP_KERNEL::Exception); void getValueOn(const double *spaceLoc, double *res) const throw(INTERP_KERNEL::Exception); void getValueOn(const double *spaceLoc, double time, double *res) const throw(INTERP_KERNEL::Exception); diff --git a/src/MEDCoupling/MEDCouplingGaussLocalization.hxx b/src/MEDCoupling/MEDCouplingGaussLocalization.hxx index 708d2ae1f..05a840b33 100644 --- a/src/MEDCoupling/MEDCouplingGaussLocalization.hxx +++ b/src/MEDCoupling/MEDCouplingGaussLocalization.hxx @@ -17,8 +17,8 @@ // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // -#ifndef __MEDCOUPLINGGAUSSLOCALIZATION_HXX__ -#define __MEDCOUPLINGGAUSSLOCALIZATION_HXX__ +#ifndef __PARAMEDMEM_MEDCOUPLINGGAUSSLOCALIZATION_HXX__ +#define __PARAMEDMEM_MEDCOUPLINGGAUSSLOCALIZATION_HXX__ #include "NormalizedUnstructuredMesh.hxx" #include "InterpKernelException.hxx" diff --git a/src/MEDCoupling/MEDCouplingMemArray.cxx b/src/MEDCoupling/MEDCouplingMemArray.cxx index 84ade1a95..24c1bb171 100644 --- a/src/MEDCoupling/MEDCouplingMemArray.cxx +++ b/src/MEDCoupling/MEDCouplingMemArray.cxx @@ -101,6 +101,33 @@ DataArrayInt *DataArrayDouble::convertToIntArr() const return ret; } +void DataArrayDouble::renumberInPlace(const int *old2New) +{ + int nbTuples=getNumberOfTuples(); + int nbOfCompo=getNumberOfComponents(); + double *tmp=new double[nbTuples*nbOfCompo]; + const double *iptr=getConstPointer(); + for(int i=0;ialloc(nbTuples,nbOfCompo); + ret->copyStringInfoFrom(*this); + const double *iptr=getConstPointer(); + double *optr=ret->getPointer(); + for(int i=0;iend()) will be kept. @@ -375,6 +402,33 @@ void DataArrayInt::useArray(const int *array, bool ownership, DeallocType type, declareAsNew(); } +void DataArrayInt::renumberInPlace(const int *old2New) +{ + int nbTuples=getNumberOfTuples(); + int nbOfCompo=getNumberOfComponents(); + int *tmp=new int[nbTuples*nbOfCompo]; + const int *iptr=getConstPointer(); + for(int i=0;ialloc(nbTuples,nbOfCompo); + ret->copyStringInfoFrom(*this); + const int *iptr=getConstPointer(); + int *optr=ret->getPointer(); + for(int i=0;igetMeshDimension()) + throw INTERP_KERNEL::Exception("checkFastEquivalWith : Mesh dimensions are not equal !"); + if(getSpaceDimension()!=other->getSpaceDimension()) + throw INTERP_KERNEL::Exception("checkFastEquivalWith : Space dimensions are not equal !"); + if(getNumberOfCells()!=other->getNumberOfCells()) + throw INTERP_KERNEL::Exception("checkFastEquivalWith : number of cells are not equal !"); +} + +/*! + * This method is very poor and looks only if 'this' and 'other' are candidate for merge of fields lying repectively on them. + */ +bool MEDCouplingMesh::areCompatibleForMerge(const MEDCouplingMesh *other) const { if(getMeshDimension()!=other->getMeshDimension()) return false; @@ -80,6 +143,15 @@ MEDCouplingFieldDouble *MEDCouplingMesh::fillFromAnalytic(TypeOfField t, int nbO return ret; } +/*! + * This method copyies all tiny strings from other (name and components name). + * @throw if other and this have not same mesh type. + */ +void MEDCouplingMesh::copyTinyStringsFrom(const MEDCouplingMesh *other) throw(INTERP_KERNEL::Exception) +{ + _name=other->_name; +} + /*! * This method builds a field lying on 'this' with 'nbOfComp' components. * 'func' is a string that is the expression to evaluate. diff --git a/src/MEDCoupling/MEDCouplingMesh.hxx b/src/MEDCoupling/MEDCouplingMesh.hxx index 44f1a73d7..378d8e4de 100644 --- a/src/MEDCoupling/MEDCouplingMesh.hxx +++ b/src/MEDCoupling/MEDCouplingMesh.hxx @@ -47,8 +47,17 @@ namespace ParaMEDMEM public: void setName(const char *name) { _name=name; } const char *getName() const { return _name.c_str(); } + virtual MEDCouplingMesh *deepCpy() const = 0; virtual MEDCouplingMeshType getType() const = 0; + virtual void copyTinyStringsFrom(const MEDCouplingMesh *other) throw(INTERP_KERNEL::Exception); + // comparison methods virtual bool isEqual(const MEDCouplingMesh *other, double prec) const { return _name==other->_name; } + virtual void checkDeepEquivalWith(const MEDCouplingMesh *other, int cellCompPol, double prec, + DataArrayInt *&cellCor, DataArrayInt *&nodeCor) const throw(INTERP_KERNEL::Exception) = 0; + virtual void checkFastEquivalWith(const MEDCouplingMesh *other, double prec) const throw(INTERP_KERNEL::Exception); + void checkGeoEquivalWith(const MEDCouplingMesh *other, int levOfCheck, double prec, + DataArrayInt *&cellCor, DataArrayInt *&nodeCor) const throw(INTERP_KERNEL::Exception); + // virtual void checkCoherency() const throw(INTERP_KERNEL::Exception) = 0; virtual bool isStructured() const = 0; virtual int getNumberOfCells() const = 0; @@ -71,8 +80,9 @@ namespace ParaMEDMEM virtual MEDCouplingFieldDouble *buildOrthogonalField() const = 0; virtual void rotate(const double *center, const double *vector, double angle) = 0; virtual void translate(const double *vector) = 0; + virtual void renumberCells(const int *old2NewBg, const int *old2NewEnd, bool check) throw(INTERP_KERNEL::Exception) = 0; virtual MEDCouplingMesh *mergeMyselfWith(const MEDCouplingMesh *other) const = 0; - virtual bool areCompatible(const MEDCouplingMesh *other) const; + virtual bool areCompatibleForMerge(const MEDCouplingMesh *other) const; static MEDCouplingMesh *mergeMeshes(const MEDCouplingMesh *mesh1, const MEDCouplingMesh *mesh2); //serialisation-unserialization virtual void getTinySerializationInformation(std::vector& tinyInfo, std::vector& littleStrings) const = 0; diff --git a/src/MEDCoupling/MEDCouplingNatureOfField.hxx b/src/MEDCoupling/MEDCouplingNatureOfField.hxx index bdc1c78dd..26bb01893 100644 --- a/src/MEDCoupling/MEDCouplingNatureOfField.hxx +++ b/src/MEDCoupling/MEDCouplingNatureOfField.hxx @@ -17,8 +17,8 @@ // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // -#ifndef __MEDCOUPLINGNATUREOFFIELD_HXX__ -#define __MEDCOUPLINGNATUREOFFIELD_HXX__ +#ifndef __PARAMEDMEM_MEDCOUPLINGNATUREOFFIELD_HXX__ +#define __PARAMEDMEM_MEDCOUPLINGNATUREOFFIELD_HXX__ namespace ParaMEDMEM { diff --git a/src/MEDCoupling/MEDCouplingNormalizedCartesianMesh.hxx b/src/MEDCoupling/MEDCouplingNormalizedCartesianMesh.hxx index ebb12e01a..c2643f327 100644 --- a/src/MEDCoupling/MEDCouplingNormalizedCartesianMesh.hxx +++ b/src/MEDCoupling/MEDCouplingNormalizedCartesianMesh.hxx @@ -17,8 +17,8 @@ // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // -#ifndef __MEDCouplingNormalizedCartesianMesh_HXX__ -#define __MEDCouplingNormalizedCartesianMesh_HXX__ +#ifndef __PARAMEDMEM_MEDCOUPLINGNORMALIZEDCARTESIANMESH_HXX__ +#define __PARAMEDMEM_MEDCOUPLINGNORMALIZEDCARTESIANMESH_HXX__ #include "NormalizedUnstructuredMesh.hxx" diff --git a/src/MEDCoupling/MEDCouplingNormalizedUnstructuredMesh.hxx b/src/MEDCoupling/MEDCouplingNormalizedUnstructuredMesh.hxx index b3ebcbc7f..5ebd3ae19 100644 --- a/src/MEDCoupling/MEDCouplingNormalizedUnstructuredMesh.hxx +++ b/src/MEDCoupling/MEDCouplingNormalizedUnstructuredMesh.hxx @@ -17,8 +17,8 @@ // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // -#ifndef __MEDCOUPLINGNORMALIZEDUNSTRUCTUREDMESH_HXX__ -#define __MEDCOUPLINGNORMALIZEDUNSTRUCTUREDMESH_HXX__ +#ifndef __PARAMEDMEM_MEDCOUPLINGNORMALIZEDUNSTRUCTUREDMESH_HXX__ +#define __PARAMEDMEM_MEDCOUPLINGNORMALIZEDUNSTRUCTUREDMESH_HXX__ #include "NormalizedUnstructuredMesh.hxx" diff --git a/src/MEDCoupling/MEDCouplingPointSet.cxx b/src/MEDCoupling/MEDCouplingPointSet.cxx index 819584473..b9c07a98b 100644 --- a/src/MEDCoupling/MEDCouplingPointSet.cxx +++ b/src/MEDCoupling/MEDCouplingPointSet.cxx @@ -98,6 +98,20 @@ DataArrayDouble *MEDCouplingPointSet::getCoordinatesAndOwner() const return _coords; } +/*! + * This method copyies all tiny strings from other (name and components name). + * @throw if other and this have not same mesh type. + */ +void MEDCouplingPointSet::copyTinyStringsFrom(const MEDCouplingMesh *other) throw(INTERP_KERNEL::Exception) +{ + const MEDCouplingPointSet *otherC=dynamic_cast(other); + if(!otherC) + throw INTERP_KERNEL::Exception("MEDCouplingPointSet::copyTinyStringsFrom : meshes have not same type !"); + MEDCouplingMesh::copyTinyStringsFrom(other); + if(_coords && otherC->_coords) + _coords->copyStringInfoFrom(*otherC->_coords); +} + bool MEDCouplingPointSet::isEqual(const MEDCouplingMesh *other, double prec) const { const MEDCouplingPointSet *otherC=dynamic_cast(other); diff --git a/src/MEDCoupling/MEDCouplingPointSet.hxx b/src/MEDCoupling/MEDCouplingPointSet.hxx index 77bf87f04..85d565f55 100644 --- a/src/MEDCoupling/MEDCouplingPointSet.hxx +++ b/src/MEDCoupling/MEDCouplingPointSet.hxx @@ -49,6 +49,7 @@ namespace ParaMEDMEM void setCoords(DataArrayDouble *coords); DataArrayDouble *getCoords() const { return _coords; } DataArrayDouble *getCoordinatesAndOwner() const; + void copyTinyStringsFrom(const MEDCouplingMesh *other) throw(INTERP_KERNEL::Exception); bool isEqual(const MEDCouplingMesh *other, double prec) const; bool areCoordsEqual(const MEDCouplingPointSet& other, double prec) const; virtual DataArrayInt *mergeNodes(double precision, bool& areNodesMerged) = 0; diff --git a/src/MEDCoupling/MEDCouplingTimeDiscretization.cxx b/src/MEDCoupling/MEDCouplingTimeDiscretization.cxx index aa307459c..b8da1c812 100644 --- a/src/MEDCoupling/MEDCouplingTimeDiscretization.cxx +++ b/src/MEDCoupling/MEDCouplingTimeDiscretization.cxx @@ -86,7 +86,22 @@ bool MEDCouplingTimeDiscretization::areCompatible(const MEDCouplingTimeDiscretiz return true; } -bool MEDCouplingTimeDiscretization::areCompatibleForMul(const MEDCouplingTimeDiscretization *other) const +bool MEDCouplingTimeDiscretization::areStrictlyCompatible(const MEDCouplingTimeDiscretization *other) const +{ + if(std::fabs(_time_tolerance-other->_time_tolerance)>1.e-16) + return false; + if(_array==0 && other->_array==0) + return true; + if(_array==0 || other->_array==0) + return false; + if(_array->getNumberOfComponents()!=other->_array->getNumberOfComponents()) + return false; + if(_array->getNumberOfTuples()!=other->_array->getNumberOfTuples()) + return false; + return true; +} + +bool MEDCouplingTimeDiscretization::areStrictlyCompatibleForMul(const MEDCouplingTimeDiscretization *other) const { if(std::fabs(_time_tolerance-other->_time_tolerance)>1.e-16) return false; @@ -104,7 +119,7 @@ bool MEDCouplingTimeDiscretization::areCompatibleForMul(const MEDCouplingTimeDis bool MEDCouplingTimeDiscretization::isEqual(const MEDCouplingTimeDiscretization *other, double prec) const { - if(!areCompatible(other)) + if(!areStrictlyCompatible(other)) return false; if(_array==other->_array) return true; @@ -366,9 +381,17 @@ bool MEDCouplingNoTimeLabel::areCompatible(const MEDCouplingTimeDiscretization * return otherC!=0; } -bool MEDCouplingNoTimeLabel::areCompatibleForMul(const MEDCouplingTimeDiscretization *other) const +bool MEDCouplingNoTimeLabel::areStrictlyCompatible(const MEDCouplingTimeDiscretization *other) const { - if(!MEDCouplingTimeDiscretization::areCompatibleForMul(other)) + if(!MEDCouplingTimeDiscretization::areStrictlyCompatible(other)) + return false; + const MEDCouplingNoTimeLabel *otherC=dynamic_cast(other); + return otherC!=0; +} + +bool MEDCouplingNoTimeLabel::areStrictlyCompatibleForMul(const MEDCouplingTimeDiscretization *other) const +{ + if(!MEDCouplingTimeDiscretization::areStrictlyCompatibleForMul(other)) return false; const MEDCouplingNoTimeLabel *otherC=dynamic_cast(other); return otherC!=0; @@ -570,19 +593,23 @@ bool MEDCouplingWithTimeStep::areCompatible(const MEDCouplingTimeDiscretization if(!MEDCouplingTimeDiscretization::areCompatible(other)) return false; const MEDCouplingWithTimeStep *otherC=dynamic_cast(other); - if(!otherC) - return false; - return std::fabs(_time-otherC->_time)<_time_tolerance; + return otherC!=0; } -bool MEDCouplingWithTimeStep::areCompatibleForMul(const MEDCouplingTimeDiscretization *other) const +bool MEDCouplingWithTimeStep::areStrictlyCompatible(const MEDCouplingTimeDiscretization *other) const { - if(!MEDCouplingTimeDiscretization::areCompatibleForMul(other)) + if(!MEDCouplingTimeDiscretization::areStrictlyCompatible(other)) return false; const MEDCouplingWithTimeStep *otherC=dynamic_cast(other); - if(!otherC) + return otherC!=0; +} + +bool MEDCouplingWithTimeStep::areStrictlyCompatibleForMul(const MEDCouplingTimeDiscretization *other) const +{ + if(!MEDCouplingTimeDiscretization::areStrictlyCompatibleForMul(other)) return false; - return std::fabs(_time-otherC->_time)<_time_tolerance; + const MEDCouplingWithTimeStep *otherC=dynamic_cast(other); + return otherC!=0; } bool MEDCouplingWithTimeStep::isEqual(const MEDCouplingTimeDiscretization *other, double prec) const @@ -839,19 +866,23 @@ bool MEDCouplingConstOnTimeInterval::areCompatible(const MEDCouplingTimeDiscreti if(!MEDCouplingTimeDiscretization::areCompatible(other)) return false; const MEDCouplingConstOnTimeInterval *otherC=dynamic_cast(other); - if(!otherC) - return false; - return (std::fabs(_start_time-otherC->_start_time)<_time_tolerance && std::fabs(_end_time-otherC->_end_time)<_time_tolerance); + return otherC!=0; } -bool MEDCouplingConstOnTimeInterval::areCompatibleForMul(const MEDCouplingTimeDiscretization *other) const +bool MEDCouplingConstOnTimeInterval::areStrictlyCompatible(const MEDCouplingTimeDiscretization *other) const { - if(!MEDCouplingTimeDiscretization::areCompatible(other)) + if(!MEDCouplingTimeDiscretization::areStrictlyCompatible(other)) return false; const MEDCouplingConstOnTimeInterval *otherC=dynamic_cast(other); - if(!otherC) + return otherC!=0; +} + +bool MEDCouplingConstOnTimeInterval::areStrictlyCompatibleForMul(const MEDCouplingTimeDiscretization *other) const +{ + if(!MEDCouplingTimeDiscretization::areStrictlyCompatible(other)) return false; - return (std::fabs(_start_time-otherC->_start_time)<_time_tolerance && std::fabs(_end_time-otherC->_end_time)<_time_tolerance); + const MEDCouplingConstOnTimeInterval *otherC=dynamic_cast(other); + return otherC!=0; } bool MEDCouplingConstOnTimeInterval::isEqual(const MEDCouplingTimeDiscretization *other, double prec) const @@ -1041,6 +1072,13 @@ MEDCouplingTwoTimeSteps::MEDCouplingTwoTimeSteps(const MEDCouplingTwoTimeSteps& _end_array=0; } +void MEDCouplingTwoTimeSteps::updateTime() +{ + MEDCouplingTimeDiscretization::updateTime(); + if(_end_array) + updateTimeWith(*_end_array); +} + void MEDCouplingTwoTimeSteps::copyTinyAttrFrom(const MEDCouplingTimeDiscretization& other) { MEDCouplingTimeDiscretization::copyTinyAttrFrom(other); @@ -1259,9 +1297,17 @@ bool MEDCouplingLinearTime::areCompatible(const MEDCouplingTimeDiscretization *o return otherC!=0; } -bool MEDCouplingLinearTime::areCompatibleForMul(const MEDCouplingTimeDiscretization *other) const +bool MEDCouplingLinearTime::areStrictlyCompatible(const MEDCouplingTimeDiscretization *other) const +{ + if(!MEDCouplingTimeDiscretization::areStrictlyCompatible(other)) + return false; + const MEDCouplingLinearTime *otherC=dynamic_cast(other); + return otherC!=0; +} + +bool MEDCouplingLinearTime::areStrictlyCompatibleForMul(const MEDCouplingTimeDiscretization *other) const { - if(!MEDCouplingTimeDiscretization::areCompatibleForMul(other)) + if(!MEDCouplingTimeDiscretization::areStrictlyCompatibleForMul(other)) return false; const MEDCouplingLinearTime *otherC=dynamic_cast(other); return otherC!=0; diff --git a/src/MEDCoupling/MEDCouplingTimeDiscretization.hxx b/src/MEDCoupling/MEDCouplingTimeDiscretization.hxx index 6077b7616..bdf62a4ad 100644 --- a/src/MEDCoupling/MEDCouplingTimeDiscretization.hxx +++ b/src/MEDCoupling/MEDCouplingTimeDiscretization.hxx @@ -17,8 +17,8 @@ // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // -#ifndef __MEDCOUPLINGTIMEDISCRETIZATION_HXX__ -#define __MEDCOUPLINGTIMEDISCRETIZATION_HXX__ +#ifndef __PARAMEDMEM_MEDCOUPLINGTIMEDISCRETIZATION_HXX__ +#define __PARAMEDMEM_MEDCOUPLINGTIMEDISCRETIZATION_HXX__ #include "MEDCoupling.hxx" #include "MEDCouplingTimeLabel.hxx" @@ -43,7 +43,8 @@ namespace ParaMEDMEM virtual void copyTinyAttrFrom(const MEDCouplingTimeDiscretization& other); virtual void checkCoherency() const throw(INTERP_KERNEL::Exception); virtual bool areCompatible(const MEDCouplingTimeDiscretization *other) const; - virtual bool areCompatibleForMul(const MEDCouplingTimeDiscretization *other) const; + virtual bool areStrictlyCompatible(const MEDCouplingTimeDiscretization *other) const; + virtual bool areStrictlyCompatibleForMul(const MEDCouplingTimeDiscretization *other) const; virtual bool isEqual(const MEDCouplingTimeDiscretization *other, double prec) const; virtual MEDCouplingTimeDiscretization *buildNewTimeReprFromThis(const MEDCouplingTimeDiscretization *other, TypeOfTimeDiscretization type, bool deepCpy) const; @@ -116,7 +117,8 @@ namespace ParaMEDMEM void divideEqual(const MEDCouplingTimeDiscretization *other); bool isEqual(const MEDCouplingTimeDiscretization *other, double prec) const; bool areCompatible(const MEDCouplingTimeDiscretization *other) const; - bool areCompatibleForMul(const MEDCouplingTimeDiscretization *other) const; + bool areStrictlyCompatible(const MEDCouplingTimeDiscretization *other) const; + bool areStrictlyCompatibleForMul(const MEDCouplingTimeDiscretization *other) const; MEDCouplingTimeDiscretization *performCpy(bool deepCpy) const; void checkNoTimePresence() const throw(INTERP_KERNEL::Exception) { } void checkTimePresence(double time) const throw(INTERP_KERNEL::Exception); @@ -155,7 +157,8 @@ namespace ParaMEDMEM void divideEqual(const MEDCouplingTimeDiscretization *other); bool isEqual(const MEDCouplingTimeDiscretization *other, double prec) const; bool areCompatible(const MEDCouplingTimeDiscretization *other) const; - bool areCompatibleForMul(const MEDCouplingTimeDiscretization *other) const; + bool areStrictlyCompatible(const MEDCouplingTimeDiscretization *other) const; + bool areStrictlyCompatibleForMul(const MEDCouplingTimeDiscretization *other) const; void getTinySerializationIntInformation(std::vector& tinyInfo) const; void getTinySerializationDbleInformation(std::vector& tinyInfo) const; void finishUnserialization(const std::vector& tinyInfoI, const std::vector& tinyInfoD, const std::vector& tinyInfoS); @@ -191,7 +194,8 @@ namespace ParaMEDMEM void finishUnserialization(const std::vector& tinyInfoI, const std::vector& tinyInfoD, const std::vector& tinyInfoS); MEDCouplingTimeDiscretization *performCpy(bool deepCpy) const; bool areCompatible(const MEDCouplingTimeDiscretization *other) const; - bool areCompatibleForMul(const MEDCouplingTimeDiscretization *other) const; + bool areStrictlyCompatible(const MEDCouplingTimeDiscretization *other) const; + bool areStrictlyCompatibleForMul(const MEDCouplingTimeDiscretization *other) const; bool isEqual(const MEDCouplingTimeDiscretization *other, double prec) const; std::vector< const DataArrayDouble *> getArraysForTime(double time) const throw(INTERP_KERNEL::Exception); void getValueForTime(double time, const std::vector& vals, double *res) const; @@ -233,6 +237,7 @@ namespace ParaMEDMEM MEDCouplingTwoTimeSteps(); ~MEDCouplingTwoTimeSteps(); public: + void updateTime(); void copyTinyAttrFrom(const MEDCouplingTimeDiscretization& other); DataArrayDouble *getEndArray() const; void checkCoherency() const throw(INTERP_KERNEL::Exception); @@ -274,7 +279,8 @@ namespace ParaMEDMEM void checkCoherency() const throw(INTERP_KERNEL::Exception); MEDCouplingTimeDiscretization *performCpy(bool deepCpy) const; bool areCompatible(const MEDCouplingTimeDiscretization *other) const; - bool areCompatibleForMul(const MEDCouplingTimeDiscretization *other) const; + bool areStrictlyCompatible(const MEDCouplingTimeDiscretization *other) const; + bool areStrictlyCompatibleForMul(const MEDCouplingTimeDiscretization *other) const; void getValueForTime(double time, const std::vector& vals, double *res) const; void getValueOnTime(int eltId, double time, double *value) const throw(INTERP_KERNEL::Exception); void getValueOnDiscTime(int eltId, int dt, int it, double *value) const throw(INTERP_KERNEL::Exception); diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index cd201a5df..05754cec4 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -25,6 +25,7 @@ #include "PointLocatorAlgos.txx" #include "BBTree.txx" #include "DirectedBoundingBox.hxx" +#include "MEDCouplingAutoRefCountObjectPtr.hxx" #include #include @@ -51,6 +52,11 @@ MEDCouplingUMesh *MEDCouplingUMesh::New(const char *meshName, int meshDim) return ret; } +MEDCouplingMesh *MEDCouplingUMesh::deepCpy() const +{ + return clone(true); +} + MEDCouplingUMesh *MEDCouplingUMesh::clone(bool recDeepCpy) const { return new MEDCouplingUMesh(*this,recDeepCpy); @@ -87,6 +93,20 @@ void MEDCouplingUMesh::checkCoherency() const throw(INTERP_KERNEL::Exception) throw INTERP_KERNEL::Exception(message.str().c_str()); } } + if(_nodal_connec) + { + if(_nodal_connec->getNumberOfComponents()!=1) + throw INTERP_KERNEL::Exception("Nodal connectivity array is expected to be with number of components set to one !"); + if(_nodal_connec->getInfoOnComponent(0)!="") + throw INTERP_KERNEL::Exception("Nodal connectivity array is expected to have no info on its single component !"); + } + if(_nodal_connec_index) + { + if(_nodal_connec_index->getNumberOfComponents()!=1) + throw INTERP_KERNEL::Exception("Nodal connectivity index array is expected to be with number of components set to one !"); + if(_nodal_connec_index->getInfoOnComponent(0)!="") + throw INTERP_KERNEL::Exception("Nodal connectivity index array is expected to have no info on its single component !"); + } } void MEDCouplingUMesh::setMeshDimension(int meshDim) @@ -181,6 +201,70 @@ bool MEDCouplingUMesh::isEqual(const MEDCouplingMesh *other, double prec) const return true; } +/*! + * This method looks if 'this' and 'other' are geometrically equivalent that is to say if each cell in 'other' correspond to one cell and only one + * in 'this' is found regarding 'prec' parameter and 'cellCompPol' parameter. + * + * In case of success cellCor and nodeCor are informed both. + * @param cellCompPol values are described in MEDCouplingUMesh::zipConnectivityTraducer method. + * @param cellCor output array giving the correspondance of cells from 'other' to 'this'. + * @param nodeCor output array giving the correspondance of nodes from 'other' to 'this'. + */ +void MEDCouplingUMesh::checkDeepEquivalWith(const MEDCouplingMesh *other, int cellCompPol, double prec, + DataArrayInt *&cellCor, DataArrayInt *&nodeCor) const throw(INTERP_KERNEL::Exception) +{ + const MEDCouplingUMesh *otherC=dynamic_cast(other); + if(!otherC) + throw INTERP_KERNEL::Exception("checkDeepEquivalWith : Two meshes are not not unstructured !"); + checkFastEquivalWith(other,prec); + if(_types!=otherC->_types) + throw INTERP_KERNEL::Exception("checkDeepEquivalWith : Types are not equal !"); + MEDCouplingAutoRefCountObjectPtr m=mergeUMeshes(this,otherC); + bool areNodesMerged; + MEDCouplingAutoRefCountObjectPtr da=m->mergeNodes(prec,areNodesMerged); + if(!areNodesMerged) + throw INTERP_KERNEL::Exception("checkDeepEquivalWith : Nodes are incompatible ! "); + int maxId=*std::max_element(da->getConstPointer(),da->getConstPointer()+getNumberOfNodes()); + const int *pt=std::find_if(da->getConstPointer()+getNumberOfNodes(),da->getConstPointer()+da->getNbOfElems(),std::bind2nd(std::greater(),maxId)); + if(pt!=da->getConstPointer()+da->getNbOfElems()) + throw INTERP_KERNEL::Exception("checkDeepEquivalWith : some nodes in other are not in this !"); + nodeCor=DataArrayInt::New(); + nodeCor->alloc(otherC->getNumberOfNodes(),1); + std::copy(da->getConstPointer()+getNumberOfNodes(),da->getConstPointer()+da->getNbOfElems(),nodeCor->getPointer()); + // + da=m->zipConnectivityTraducer(cellCompPol); + maxId=*std::max_element(da->getConstPointer(),da->getConstPointer()+getNumberOfCells()); + pt=std::find_if(da->getConstPointer()+getNumberOfCells(),da->getConstPointer()+da->getNbOfElems(),std::bind2nd(std::greater(),maxId)); + if(pt!=da->getConstPointer()+da->getNbOfElems()) + { + nodeCor->decrRef(); nodeCor=0; + throw INTERP_KERNEL::Exception("checkDeepEquivalWith : some cells in other are not in this !"); + } + cellCor=DataArrayInt::New(); + cellCor->alloc(otherC->getNumberOfCells(),1); + std::copy(da->getConstPointer()+getNumberOfCells(),da->getConstPointer()+da->getNbOfElems(),cellCor->getPointer()); +} + +/*! + * This method checks fastly that 'this' and 'other' are equal. + */ +void MEDCouplingUMesh::checkFastEquivalWith(const MEDCouplingMesh *other, double prec) const throw(INTERP_KERNEL::Exception) +{ + const MEDCouplingUMesh *otherC=dynamic_cast(other); + if(!otherC) + throw INTERP_KERNEL::Exception("checkFastEquivalWith : Two meshes are not not unstructured !"); + MEDCouplingPointSet::checkFastEquivalWith(other,prec); + int nbOfCells=getNumberOfCells(); + if(nbOfCells<1) + return ; + bool status=true; + status&=areCellsFrom2MeshEqual(otherC,0,prec); + status&=areCellsFrom2MeshEqual(otherC,nbOfCells/2,prec); + status&=areCellsFrom2MeshEqual(otherC,nbOfCells-1,prec); + if(!status) + throw INTERP_KERNEL::Exception("checkFastEquivalWith : Two meshes are not equal because on 3 test cells some difference have been detected !"); +} + /*! * \b WARNING this method do the assumption that connectivity lies on the coordinates set. * For speed reasons no check of this will be done. @@ -438,33 +522,42 @@ DataArrayInt *MEDCouplingUMesh::zipCoordsTraducer() * This method stands if 'cell1' and 'cell2' are equals regarding 'compType' policy. * The semantic of 'compType' is specified in MEDCouplingUMesh::zipConnectivityTraducer method. */ -bool MEDCouplingUMesh::areCellsEquals(int cell1, int cell2, int compType) const +bool MEDCouplingUMesh::areCellsEqual(int cell1, int cell2, int compType) const { switch(compType) { case 0: - return areCellsEquals0(cell1,cell2); + return areCellsEqual0(cell1,cell2); case 1: - return areCellsEquals1(cell1,cell2); + return areCellsEqual1(cell1,cell2); case 2: - return areCellsEquals2(cell1,cell2); + return areCellsEqual2(cell1,cell2); } throw INTERP_KERNEL::Exception("Unknown comparison asked ! Must be in 0,1 or 2."); } -bool MEDCouplingUMesh::areCellsEquals0(int cell1, int cell2) const +/*! + * This method is the last step of the MEDCouplingUMesh::zipConnectivityTraducer with policy 0. + */ +bool MEDCouplingUMesh::areCellsEqual0(int cell1, int cell2) const { const int *conn=getNodalConnectivity()->getConstPointer(); const int *connI=getNodalConnectivityIndex()->getConstPointer(); return std::equal(conn+connI[cell1],conn+connI[cell1+1],conn+connI[cell2]); } -bool MEDCouplingUMesh::areCellsEquals1(int cell1, int cell2) const +/*! + * This method is the last step of the MEDCouplingUMesh::zipConnectivityTraducer with policy 1. + */ +bool MEDCouplingUMesh::areCellsEqual1(int cell1, int cell2) const { throw INTERP_KERNEL::Exception("Policy comparison, not implemented yet !"); } -bool MEDCouplingUMesh::areCellsEquals2(int cell1, int cell2) const +/*! + * This method is the last step of the MEDCouplingUMesh::zipConnectivityTraducer with policy 2. + */ +bool MEDCouplingUMesh::areCellsEqual2(int cell1, int cell2) const { const int *conn=getNodalConnectivity()->getConstPointer(); const int *connI=getNodalConnectivityIndex()->getConstPointer(); @@ -473,6 +566,33 @@ bool MEDCouplingUMesh::areCellsEquals2(int cell1, int cell2) const return s1==s2; } +/*! + * This method compares 2 cells coming from two unstructured meshes : 'this' and 'other'. + * This method compares 2 cells having the same id 'cellId' in 'this' and 'other'. + */ +bool MEDCouplingUMesh::areCellsFrom2MeshEqual(const MEDCouplingUMesh *other, int cellId, double prec) const +{ + if(getTypeOfCell(cellId)!=other->getTypeOfCell(cellId)) + return false; + std::vector c1,c2; + getNodeIdsOfCell(cellId,c1); + other->getNodeIdsOfCell(cellId,c2); + int sz=c1.size(); + if(sz!=c2.size()) + return false; + for(int i=0;i n1,n2; + getCoordinatesOfNode(c1[0],n1); + other->getCoordinatesOfNode(c2[0],n2); + std::transform(n1.begin(),n1.end(),n2.begin(),n1.begin(),std::minus()); + std::transform(n1.begin(),n1.end(),n1.end(),std::ptr_fun(fabs)); + if(*std::max_element(n1.begin(),n1.end())>prec) + return false; + } + return true; +} + /*! * This method find in candidate pool defined by 'candidates' the cells equal following the polycy 'compType'. * If any true is returned and the results will be put at the end of 'result' output parameter. If not false is returned @@ -481,7 +601,7 @@ bool MEDCouplingUMesh::areCellsEquals2(int cell1, int cell2) const * If in 'candidates' pool -1 value is considered as an empty value. * WARNING this method returns only ONE set of result ! */ -bool MEDCouplingUMesh::areCellsEqualsInPool(const std::vector& candidates, int compType, std::vector& result) const +bool MEDCouplingUMesh::areCellsEqualInPool(const std::vector& candidates, int compType, std::vector& result) const { std::set cand(candidates.begin(),candidates.end()); cand.erase(-1); @@ -494,7 +614,7 @@ bool MEDCouplingUMesh::areCellsEqualsInPool(const std::vector& candidates, std::set::const_iterator begin2=iter; begin2++; for(std::set::const_iterator iter2=begin2;iter2!=cand.end();iter2++) { - if(areCellsEquals(*iter,*iter2,compType)) + if(areCellsEqual(*iter,*iter2,compType)) { if(!ret) { @@ -553,7 +673,7 @@ void MEDCouplingUMesh::findCommonCellsBase(int compType, std::vector& res, for(std::vector::const_iterator iter=candidates1.begin();iter!=candidates1.end();iter++) if(!isFetched[*iter]) candidates.push_back(*iter); - if(areCellsEqualsInPool(candidates,compType,res)) + if(areCellsEqualInPool(candidates,compType,res)) { int pos=resI.back(); resI.push_back(res.size()); @@ -790,6 +910,8 @@ void MEDCouplingUMesh::findBoundaryNodes(std::vector& nodes) const * This method renumber 'this' using 'newNodeNumbers' array of size this->getNumberOfNodes. * newNbOfNodes specifies the *std::max_element(newNodeNumbers,newNodeNumbers+this->getNumberOfNodes()) * This value is asked because often known by the caller of this method. + * This method, contrary to MEDCouplingMesh::renumberCells does NOT conserve the number of nodes before and after. + * * @param newNodeNumbers array specifying the new numbering. * @param newNbOfNodes the new number of nodes. */ @@ -816,6 +938,9 @@ void MEDCouplingUMesh::renumberNodes(const int *newNodeNumbers, int newNbOfNodes * This method renumbers cells of 'this' using the array specified by [old2NewBg;old2NewEnd) * If std::distance(old2NewBg,old2NewEnd)!=this->getNumberOfCells() an INTERP_KERNEL::Exception will be thrown. * + * Contrary to MEDCouplingPointSet::renumberNodes, this method makes a permutation without any fuse of cell. + * After the call of this method the number of cells remains the same as before. + * * If 'check' equals true the method will check that any elements in [old2NewBg;old2NewEnd) is unique ; if not * an INTERP_KERNEL::Exception will be thrown. When 'check' equals true [old2NewBg;old2NewEnd) is not expected to * be strictly in [0;this->getNumberOfCells()). @@ -824,7 +949,7 @@ void MEDCouplingUMesh::renumberNodes(const int *newNodeNumbers, int newNbOfNodes * To avoid any throw of SIGSEGV when 'check' equals false, the elements in [old2NewBg;old2NewEnd) should be unique and * should be contained in[0;this->getNumberOfCells()). */ -void MEDCouplingUMesh::renumberCells(const int *old2NewBg, const int *old2NewEnd, bool check) +void MEDCouplingUMesh::renumberCells(const int *old2NewBg, const int *old2NewEnd, bool check) throw(INTERP_KERNEL::Exception) { int nbCells=getNumberOfCells(); const int *array=old2NewBg; @@ -994,6 +1119,8 @@ int MEDCouplingUMesh::getNumberOfCellsWithType(INTERP_KERNEL::NormalizedCellType /*! * Appends the nodal connectivity in 'conn' of cell with id 'cellId'. + * All elements added in conn can be used by MEDCouplingUMesh::getCoordinatesOfNode method. + * That is to say -1 separator is omitted in returned conn. */ void MEDCouplingUMesh::getNodeIdsOfCell(int cellId, std::vector& conn) const { @@ -1863,6 +1990,24 @@ void MEDCouplingUMesh::orientCorrectlyPolyhedrons() throw(INTERP_KERNEL::Excepti updateTime(); } +/*! + * This method has a sense for meshes with spaceDim==3 and meshDim==2. + * If it is not the case an exception will be thrown. + * This method is fast because the first cell of 'this' is used to compute the plane. + * @param vec output of size at least 3 used to store the normal vector (with norm equal to Area ) of searched plane. + * @param pos output of size at least 3 used to store a point owned of searched plane. + */ +void MEDCouplingUMesh::getFastMiddlePlaneOfThis(double *vec, double *pos) const throw(INTERP_KERNEL::Exception) +{ + if(getMeshDimension()!=2 || getSpaceDimension()!=3) + throw INTERP_KERNEL::Exception("Invalid mesh to apply getFastMiddlePlaneOfThis on it : must be meshDim==2 and spaceDim==3 !"); + const int *conn=_nodal_connec->getConstPointer(); + const int *connI=_nodal_connec_index->getConstPointer(); + const double *coordsPtr=_coords->getConstPointer(); + INTERP_KERNEL::areaVectorOfPolygon(conn+1,connI[1]-connI[0]-1,coordsPtr,vec); + std::copy(coordsPtr+3*conn[1],coordsPtr+3*conn[1]+3,pos); +} + /*! * This method aggregate the bbox of each cell and put it into bbox parameter. * @param bbox out parameter of size 2*spacedim*nbOfcells. @@ -2087,15 +2232,9 @@ DataArrayDouble *MEDCouplingUMesh::getBarycenterAndOwner() const const double *coor=_coords->getConstPointer(); for(int i=0;i()); - nbOfPts++; - } - ptToFill=std::transform(tmp,tmp+spaceDim,ptToFill,std::bind2nd(std::divides(),(double)nbOfPts)); + INTERP_KERNEL::NormalizedCellType type=(INTERP_KERNEL::NormalizedCellType)nodal[nodalI[i]]; + INTERP_KERNEL::computeBarycenter2(type,nodal+nodalI[i]+1,nodalI[i+1]-nodalI[i]-1,coor,spaceDim,ptToFill); + ptToFill+=spaceDim; } delete [] tmp; return ret; diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx index 1ae0fdc74..aabaf06ef 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh.hxx @@ -33,10 +33,14 @@ namespace ParaMEDMEM public: MEDCOUPLING_EXPORT static MEDCouplingUMesh *New(); MEDCOUPLING_EXPORT static MEDCouplingUMesh *New(const char *meshName, int meshDim); + MEDCOUPLING_EXPORT MEDCouplingMesh *deepCpy() const; MEDCOUPLING_EXPORT MEDCouplingUMesh *clone(bool recDeepCpy) const; MEDCOUPLING_EXPORT void updateTime(); MEDCOUPLING_EXPORT MEDCouplingMeshType getType() const { return UNSTRUCTURED; } MEDCOUPLING_EXPORT bool isEqual(const MEDCouplingMesh *other, double prec) const; + MEDCOUPLING_EXPORT void checkDeepEquivalWith(const MEDCouplingMesh *other, int cellCompPol, double prec, + DataArrayInt *&cellCor, DataArrayInt *&nodeCor) const throw(INTERP_KERNEL::Exception); + MEDCOUPLING_EXPORT void checkFastEquivalWith(const MEDCouplingMesh *other, double prec) const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT void checkCoherency() const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT void setMeshDimension(int meshDim); MEDCOUPLING_EXPORT void allocateCells(int nbOfCells); @@ -61,10 +65,11 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT void serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const; MEDCOUPLING_EXPORT void unserialization(const std::vector& tinyInfo, const DataArrayInt *a1, DataArrayDouble *a2, const std::vector& littleStrings); //tools - MEDCOUPLING_EXPORT bool areCellsEquals(int cell1, int cell2, int compType) const; - MEDCOUPLING_EXPORT bool areCellsEquals0(int cell1, int cell2) const; - MEDCOUPLING_EXPORT bool areCellsEquals1(int cell1, int cell2) const; - MEDCOUPLING_EXPORT bool areCellsEquals2(int cell1, int cell2) const; + MEDCOUPLING_EXPORT bool areCellsEqual(int cell1, int cell2, int compType) const; + MEDCOUPLING_EXPORT bool areCellsEqual0(int cell1, int cell2) const; + MEDCOUPLING_EXPORT bool areCellsEqual1(int cell1, int cell2) const; + MEDCOUPLING_EXPORT bool areCellsEqual2(int cell1, int cell2) const; + MEDCOUPLING_EXPORT bool areCellsFrom2MeshEqual(const MEDCouplingUMesh *other, int cellId, double prec) const; MEDCOUPLING_EXPORT void convertToPolyTypes(const std::vector& cellIdsToConvert); MEDCOUPLING_EXPORT DataArrayInt *zipCoordsTraducer(); MEDCOUPLING_EXPORT DataArrayInt *zipConnectivityTraducer(int compType); @@ -77,7 +82,7 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT void findBoundaryNodes(std::vector& nodes) const; MEDCOUPLING_EXPORT MEDCouplingPointSet *buildBoundaryMesh(bool keepCoords) const; MEDCOUPLING_EXPORT void renumberNodes(const int *newNodeNumbers, int newNbOfNodes); - MEDCOUPLING_EXPORT void renumberCells(const int *old2NewBg, const int *old2NewEnd, bool check); + MEDCOUPLING_EXPORT void renumberCells(const int *old2NewBg, const int *old2NewEnd, bool check) throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT void giveElemsInBoundingBox(const double *bbox, double eps, std::vector& elems); MEDCOUPLING_EXPORT void giveElemsInBoundingBox(const INTERP_KERNEL::DirectedBoundingBox& bbox, double eps, std::vector& elems); MEDCOUPLING_EXPORT MEDCouplingFieldDouble *getMeasureField(bool isAbs) const; @@ -97,6 +102,7 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT void orientCorrectly2DCells(const double *vec, bool polyOnly) throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT void arePolyhedronsNotCorrectlyOriented(std::vector& cells) const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT void orientCorrectlyPolyhedrons() throw(INTERP_KERNEL::Exception); + MEDCOUPLING_EXPORT void getFastMiddlePlaneOfThis(double *vec, double *pos) const throw(INTERP_KERNEL::Exception); //utilities for MED File RW MEDCOUPLING_EXPORT bool checkConsecutiveCellTypes() const; MEDCOUPLING_EXPORT DataArrayInt *rearrange2ConsecutiveCellTypes(); @@ -123,7 +129,7 @@ namespace ParaMEDMEM DataArrayDouble *fillExtCoordsUsingTranslation(const MEDCouplingUMesh *mesh1D, bool isQuad) const; template void findCommonCellsBase(int compType, std::vector& res, std::vector& resI) const; - bool areCellsEqualsInPool(const std::vector& candidates, int compType, std::vector& result) const; + bool areCellsEqualInPool(const std::vector& candidates, int compType, std::vector& result) const; MEDCouplingUMesh *buildPartOfMySelfKeepCoords(const int *start, const int *end) const; template void getCellsContainingPointsAlg(const double *coords, const double *pos, int nbOfPoints, diff --git a/src/MEDCoupling/MEDCouplingUMeshDesc.cxx b/src/MEDCoupling/MEDCouplingUMeshDesc.cxx index 77a8c7257..e00231fff 100644 --- a/src/MEDCoupling/MEDCouplingUMeshDesc.cxx +++ b/src/MEDCoupling/MEDCouplingUMeshDesc.cxx @@ -55,6 +55,14 @@ MEDCouplingUMeshDesc *MEDCouplingUMeshDesc::New(const char *meshName, int meshDi return ret; } +/*! + * not implemented + */ +MEDCouplingMesh *MEDCouplingUMeshDesc::deepCpy() const +{ + return 0; +} + void MEDCouplingUMeshDesc::checkCoherency() const throw(INTERP_KERNEL::Exception) { for(std::set::const_iterator iter=_types.begin();iter!=_types.end();iter++) @@ -68,6 +76,12 @@ void MEDCouplingUMeshDesc::checkCoherency() const throw(INTERP_KERNEL::Exception } } +void MEDCouplingUMeshDesc::checkDeepEquivalWith(const MEDCouplingMesh *other, int cellCompPol, double prec, + DataArrayInt *&cellCor, DataArrayInt *&nodeCor) const throw(INTERP_KERNEL::Exception) +{ + throw INTERP_KERNEL::Exception("MEDCouplingUMeshDesc::checkDeepEquivalWith : not implemented yet !"); +} + void MEDCouplingUMeshDesc::setMeshDimension(unsigned meshDim) { _mesh_dim=meshDim; @@ -331,6 +345,11 @@ MEDCouplingPointSet *MEDCouplingUMeshDesc::buildBoundaryMesh(bool keepCoords) co return 0; } +void MEDCouplingUMeshDesc::renumberCells(const int *old2NewBg, const int *old2NewEnd, bool check) throw(INTERP_KERNEL::Exception) +{ + throw INTERP_KERNEL::Exception("Available for UMesh desc but not implemented yet !"); +} + void MEDCouplingUMeshDesc::renumberNodes(const int *newNodeNumbers, int newNbOfNodes) { MEDCouplingPointSet::renumberNodes(newNodeNumbers,newNbOfNodes); diff --git a/src/MEDCoupling/MEDCouplingUMeshDesc.hxx b/src/MEDCoupling/MEDCouplingUMeshDesc.hxx index b69645c20..a2471245d 100644 --- a/src/MEDCoupling/MEDCouplingUMeshDesc.hxx +++ b/src/MEDCoupling/MEDCouplingUMeshDesc.hxx @@ -33,7 +33,10 @@ namespace ParaMEDMEM public: MEDCOUPLING_EXPORT static MEDCouplingUMeshDesc *New(); MEDCOUPLING_EXPORT static MEDCouplingUMeshDesc *New(const char *meshName, int meshDim); + MEDCOUPLING_EXPORT MEDCouplingMesh *deepCpy() const; MEDCOUPLING_EXPORT void checkCoherency() const throw(INTERP_KERNEL::Exception); + MEDCOUPLING_EXPORT void checkDeepEquivalWith(const MEDCouplingMesh *other, int cellCompPol, double prec, + DataArrayInt *&cellCor, DataArrayInt *&nodeCor) const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT void setMeshDimension(unsigned meshDim); MEDCOUPLING_EXPORT int getNumberOfCells() const; MEDCOUPLING_EXPORT int getNumberOfFaces() const; @@ -60,6 +63,7 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT MEDCouplingPointSet *buildFacePartOfMySelfNode(const int *start, const int *end, bool fullyIn) const; MEDCOUPLING_EXPORT void findBoundaryNodes(std::vector& nodes) const; MEDCOUPLING_EXPORT MEDCouplingPointSet *buildBoundaryMesh(bool keepCoords) const; + MEDCOUPLING_EXPORT void renumberCells(const int *old2NewBg, const int *old2NewEnd, bool check) throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT void renumberNodes(const int *newNodeNumbers, int newNbOfNodes); MEDCOUPLING_EXPORT MEDCouplingFieldDouble *getMeasureField(bool isAbs) const; MEDCOUPLING_EXPORT MEDCouplingFieldDouble *getMeasureFieldOnNode(bool isAbs) const; diff --git a/src/MEDCoupling/Makefile.am b/src/MEDCoupling/Makefile.am index f02342f05..ca9dd5295 100644 --- a/src/MEDCoupling/Makefile.am +++ b/src/MEDCoupling/Makefile.am @@ -36,7 +36,9 @@ MEDCouplingCMesh.hxx MEDCouplingTimeDiscretization.hxx MEDCouplingFieldDiscretization.hxx MEDCouplingPointSet.hxx MEDCouplingPointSet.txx \ MEDCouplingUMeshDesc.hxx MEDCouplingNatureOfField.hxx \ MEDCouplingNormalizedCartesianMesh.hxx MEDCouplingNormalizedCartesianMesh.txx \ -MEDCouplingRemapper.hxx MEDCouplingExtrudedMesh.hxx MEDCouplingGaussLocalization.hxx +MEDCouplingRemapper.hxx MEDCouplingExtrudedMesh.hxx MEDCouplingGaussLocalization.hxx \ +MEDCouplingAutoRefCountObjectPtr.hxx + # Libraries targets dist_libmedcoupling_la_SOURCES = \ diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx index 7fb1d01fe..aa2f7be0f 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx @@ -91,6 +91,12 @@ namespace ParaMEDMEM CPPUNIT_TEST( testGaussPointNEField1 ); CPPUNIT_TEST( testCellOrientation1 ); CPPUNIT_TEST( testCellOrientation2 ); + CPPUNIT_TEST( testPolyhedronBarycenter ); + CPPUNIT_TEST( testNormL12Integ1D ); + CPPUNIT_TEST( testAreaBary2D ); + CPPUNIT_TEST( testAreaBary3D ); + CPPUNIT_TEST( testRenumberCellsForFields ); + CPPUNIT_TEST( testRenumberNodesForFields ); //MEDCouplingBasicsTestInterp.cxx CPPUNIT_TEST( test2DInterpP0P0_1 ); CPPUNIT_TEST( test2DInterpP0P0PL_1 ); @@ -210,6 +216,12 @@ namespace ParaMEDMEM void testGaussPointNEField1(); void testCellOrientation1(); void testCellOrientation2(); + void testPolyhedronBarycenter(); + void testNormL12Integ1D(); + void testAreaBary2D(); + void testAreaBary3D(); + void testRenumberCellsForFields(); + void testRenumberNodesForFields(); //MEDCouplingBasicsTestInterp.cxx void test2DInterpP0P0_1(); void test2DInterpP0P0PL_1(); @@ -298,6 +310,10 @@ namespace ParaMEDMEM static MEDCouplingUMesh *build1DTargetMesh_2(); static MEDCouplingUMesh *build2DCurveSourceMesh_2(); static MEDCouplingUMesh *build2DCurveTargetMesh_2(); + static MEDCouplingUMesh *build1DTargetMesh_3(); + static MEDCouplingUMesh *build2DCurveTargetMesh_3(); + static MEDCouplingUMesh *build2DTargetMesh_3(); + static MEDCouplingUMesh *build3DTargetMesh_3(); static double sumAll(const std::vector< std::map >& matrix); }; } diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest0.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest0.cxx index 4d518bd94..26a5be5d0 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest0.cxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest0.cxx @@ -746,6 +746,77 @@ MEDCouplingUMesh *MEDCouplingBasicsTest::build2DCurveTargetMesh_2() return ret; } +MEDCouplingUMesh *MEDCouplingBasicsTest::build1DTargetMesh_3() +{ + MEDCouplingUMesh *ret=MEDCouplingUMesh::New("1DMesh_3",1); + ret->allocateCells(4); + int conn[10]={0,1,2, 3,4, 6,5,7 ,9,8}; + ret->insertNextCell(INTERP_KERNEL::NORM_SEG3,3,conn); + ret->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,conn+3); + ret->insertNextCell(INTERP_KERNEL::NORM_SEG3,3,conn+5); + ret->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,conn+8); + ret->finishInsertingCells(); + double coords[10]={0.5,1.,0.8,5.,5.21,0.5,1.1,0.7,5.,5.31}; + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(10,1); + std::copy(coords,coords+10,myCoords->getPointer()); + ret->setCoords(myCoords); + myCoords->decrRef(); + return ret; +} + +MEDCouplingUMesh *MEDCouplingBasicsTest::build2DCurveTargetMesh_3() +{ + MEDCouplingUMesh *ret=MEDCouplingUMesh::New("2DCurveMesh_3",1); + ret->allocateCells(4); + int conn[10]={0,1,2, 3,4, 6,5,7 ,9,8}; + ret->insertNextCell(INTERP_KERNEL::NORM_SEG3,3,conn); + ret->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,conn+3); + ret->insertNextCell(INTERP_KERNEL::NORM_SEG3,3,conn+5); + ret->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,conn+8); + ret->finishInsertingCells(); + double coords[20]={0.5,0.5,1.,1.,0.8,0.8,5.,5.,5.21,5.21,0.5,0.5,1.1,1.1,0.7,0.7,5.,5.,5.31,5.31}; + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(10,2); + std::copy(coords,coords+20,myCoords->getPointer()); + ret->setCoords(myCoords); + myCoords->decrRef(); + return ret; +} + +MEDCouplingUMesh *MEDCouplingBasicsTest::build2DTargetMesh_3() +{ + MEDCouplingUMesh *ret=MEDCouplingUMesh::New("2DMesh_3",2); + ret->allocateCells(10); + int conn[52]={ + 0,1,2, 0,1,3,4, 0,1,3,5,4, 0,1,2,6,7,8, 0,1,3,4,6,9,2,10, + 0,2,1, 0,4,3,1, 0,4,5,3,1, 0,2,1,8,7,6, 0,4,3,1,10,2,9,6 + }; + ret->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,conn); + ret->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+3); + ret->insertNextCell(INTERP_KERNEL::NORM_POLYGON,5,conn+7); + ret->insertNextCell(INTERP_KERNEL::NORM_TRI6,6,conn+12); + ret->insertNextCell(INTERP_KERNEL::NORM_QUAD8,8,conn+18); + ret->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,conn+26); + ret->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+29); + ret->insertNextCell(INTERP_KERNEL::NORM_POLYGON,5,conn+33); + ret->insertNextCell(INTERP_KERNEL::NORM_TRI6,6,conn+38); + ret->insertNextCell(INTERP_KERNEL::NORM_QUAD8,8,conn+44); + ret->finishInsertingCells(); + double coords[22]={0.,0.,1.,0.,0.5,1.,1.,1.,0.,1.,0.5,2.,0.5,0.,0.75,0.5,0.25,0.5,1.,0.5,0.,0.5}; + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(11,2); + std::copy(coords,coords+22,myCoords->getPointer()); + ret->setCoords(myCoords); + myCoords->decrRef(); + ret->checkCoherency(); + return ret; +} + +MEDCouplingUMesh *MEDCouplingBasicsTest::build3DTargetMesh_3() +{ +} + double MEDCouplingBasicsTest::sumAll(const std::vector< std::map >& matrix) { double ret=0.; diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest1.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest1.cxx index ef2f42354..ca6ab611c 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest1.cxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest1.cxx @@ -1445,7 +1445,7 @@ void MEDCouplingBasicsTest::testFillFromAnalytic() f1->accumulate(values4); CPPUNIT_ASSERT_DOUBLES_EQUAL(3.6,values4[0],1.e-12); CPPUNIT_ASSERT_DOUBLES_EQUAL(7.2,values4[1],1.e-12); - f1->measureAccumulate(true,values4); + f1->integral(true,values4); CPPUNIT_ASSERT_DOUBLES_EQUAL(0.5,values4[0],1.e-12); CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,values4[1],1.e-12); f1->decrRef(); @@ -1515,7 +1515,7 @@ void MEDCouplingBasicsTest::testFillFromAnalytic2() f1->accumulate(values4); CPPUNIT_ASSERT_DOUBLES_EQUAL(3.6,values4[0],1.e-12); CPPUNIT_ASSERT_DOUBLES_EQUAL(7.2,values4[1],1.e-12); - f1->measureAccumulate(true,values4); + f1->integral(true,values4); CPPUNIT_ASSERT_DOUBLES_EQUAL(0.5,values4[0],1.e-12); CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,values4[1],1.e-12); f1->decrRef(); diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest2.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest2.cxx index b3f96a287..046ef0337 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest2.cxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest2.cxx @@ -27,6 +27,7 @@ #include #include +#include using namespace ParaMEDMEM; @@ -223,6 +224,30 @@ void MEDCouplingBasicsTest::testCellOrientation2() res1.clear(); m5->arePolyhedronsNotCorrectlyOriented(res1); CPPUNIT_ASSERT(res1.empty()); + MEDCouplingFieldDouble *f3=m5->getMeasureField(false); + CPPUNIT_ASSERT_EQUAL(15,f3->getArray()->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(1,f3->getNumberOfComponents()); + const double *f3Ptr=f3->getArray()->getConstPointer(); + const double expected1[15]={ + 0.075,0.0375,0.0375,0.075,0.075, + 0.1125,0.05625,0.05625,0.1125,0.1125, + 0.0625,0.03125,0.03125,0.0625,0.0625 + }; + for(int i=0;i<15;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(std::abs(expected1[i]),f3Ptr[i],1e-12); + f3->decrRef(); + DataArrayDouble *f4=m5->getBarycenterAndOwner(); + CPPUNIT_ASSERT_EQUAL(15,f4->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(3,f4->getNumberOfComponents()); + const double *f4Ptr=f4->getConstPointer(); + const double expected2[45]={ + -0.05,-0.05,0.15, 0.3666666666666667,-0.13333333333333333,0.15, 0.53333333333333333,0.033333333333333333,0.15, -0.05,0.45,0.15, 0.45,0.45,0.15, + -0.05,-0.05,0.525, 0.3666666666666667,-0.13333333333333333,0.525, 0.53333333333333333,0.033333333333333333,0.525, -0.05,0.45,0.525, 0.45,0.45,0.525, + -0.05,-0.05,0.875, 0.3666666666666667,-0.13333333333333333,0.875, 0.53333333333333333,0.033333333333333333,0.875, -0.05,0.45,0.875, 0.45,0.45,0.875 + }; + for(int i=0;i<45;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected2[i],f4Ptr[i],1e-12); + f4->decrRef(); m5->decrRef(); m3->decrRef(); m4->decrRef(); @@ -230,3 +255,420 @@ void MEDCouplingBasicsTest::testCellOrientation2() f2->decrRef(); m2->decrRef(); } + +/*! + * This test check polyhedron true barycenter computation. + */ +void MEDCouplingBasicsTest::testPolyhedronBarycenter() +{ + /*double coords [] = { 0.241310763507,0.0504777305619,0.0682283524903, 0.252501053866,-0.0625176732937,0.137272639894, 0.152262663601,0.241816569527,0.133812556197, 0.18047750211,-0.0789949051358,0.339098173401, 0.151741971857,0.238885278571,0.137715037333, 0.242532155481,-0.0928169086456,0.0678043417367, 0.240941965335,-0.015461491464,0.0617186345825, 0.24127650112,0.0499427876717,0.0679634099148, -0.145828917428,0.206291632565,0.0310071927543, 0.0125651775307,0.266262085828,0.105228430543, -0.0994066533286,0.233224271238,0.0572213839567, -0.0951345338317,0.234819509426,0.0592126284538, 0.136580574205,-0.205486212579,0.0572866072014, 0.0637270784978,-0.168886355238,0.446614057077, 0.041337157151,-0.213402568198,0.372407095999, 0.0411601970268,-0.202387875756,0.411334979491, -0.108355701857,0.193636239335,0.204886756738, 0.00639779029829,0.155296981517,0.252585892979, 0.0262473111702,-0.112919732543,0.424286639249, -0.224103052733,-0.139430015438,-0.0122352295701, -0.0312760589481,-0.274272003594,0.0323959636568, -0.166663422532,-0.217754445175,0.00392109070364, -0.30586619777,-0.0475168041091,-0.0144585228182, -0.280881480586,0.135571293538,0.00623923647986, -0.25548538234,0.156819217766,0.0645277879769, -0.131567009284,0.184133752309,0.206021802753, -0.196204010965,0.151602971681,0.212974777736, -0.183713879463,0.0802946639531,0.260115662599, -0.244241178767,-0.0738873389604,0.144590565817, -0.155804057829,-0.164892720025,0.210613950558, -0.170950800428,-0.215099334026,0.00610122860092, -0.30552634869,-0.0490020791904,-0.0132786533145, 0.271831011884,0.15105657296,0.0230534827908, 0.281919192283,0.0898544306288,-0.0625201489143, 0.260240727276,-0.0120688706637,-0.0532316588626, 0.244947737722,0.0197984684293,0.0309341209233, 0.23439631578,0.229825279875,0.0508520585381, 0.160921316875,0.265078502128,0.121716560626, -0.315088694175,0.0747700471918,-0.245836615071, -0.327728781776,0.0857114674649,-0.239431905957, -0.308385460634,0.145142997084,-0.149886828433, 0.0488236045164,0.309462801914,0.0849169148265, -0.0244964803395,0.33145611751,-0.0476415818061, 0.0060567994229,0.32418412014,0.0367779543812, -0.0950221448063,0.236675326003,0.0572594453983, 0.248723023186,0.0886648784791,-0.176629430538, 0.116796984,0.256596599567,-0.292863523603, 0.118024552914,0.229154257843,-0.34233232501, 0.217507892549,-0.0417822335742,-0.176771782888, -0.224429321304,0.0125595300114,-0.362064725588, 0.0937301100955,-0.0500824832657,-0.299713548444, -0.244162220397,0.0383853931293,-0.389856984411, -0.0281989366102,0.097392811563,-0.458244577284, -0.385010847162,0.10122766194,-0.140052859922, -0.377936358012,0.110875172128,-0.176207095463, 0.244483045556,-0.0991073977045,0.0575134372934, 0.262605120167,-0.100243191645,-0.0495620806935, 0.240306880972,-0.136153701579,-0.114745281696, 0.215763176129,-0.0836766059189,-0.183249640616, 0.237870396603,-0.132449578286,-0.121598854639, -0.0637683083097,-0.27921020214,-0.149112321992, -0.0856211014977,-0.2973233473,-0.0446878139589, 0.104675342288,-0.0625908305324,-0.290346256534, 0.0248264249186,-0.247797708548,-0.165830884019, 0.0719302438309,-0.178468260473,-0.211432157345, 0.142871843159,-0.208769948542,0.0454101128246, 0.167803379307,-0.207851396623,-0.088802726124, 0.12868717152,-0.230920439715,0.00760508389036, -0.0372812069535,-0.286740286332,0.00963701291166 }; + + int connN [347] = { + //polyhedron 0 + 0, 1, 2, 3, 4, -1, 1, 5, 0, 6, 7, -1, 0, 7, 2, 8, 9, 10, 11, -1, 1, 5, 3, 12, 13, 14, 15, -1, 16, 9, 17, 2, 4, -1, 4, 3, 17, 13, 18, -1, 5, 6, 12, 19, 20, 21, -1, 6, 7, 19, 8, 22, 23, -1, 23, 24, 8, 10, -1, 25, 11, 16, 9, -1, 24, 26, 10, 25, 11, -1, 12, 14, 20, -1, 27, 28, 18, 29, 13, 15, -1, 14, 15, 20, 29, 21, 30, -1, 26, 27, 25, 18, 16, 17, -1, 22, 19, 31, 21, 30, -1, 22, 31, 23, 28, 24, 27, 26, -1, 31, 30, 28, 29, + //polyhedron 1 + 0, 7, 2, 8, 9, 10, 11, -1, 32, 0, 33, 7, 34, 35, -1, 32, 0, 36, 2, 37, -1, 35, 7, 38, 8, 39, 40, -1, 2, 37, 9, 41, -1, 40, 8, 42, 10, 43, 44, -1, 41, 9, 43, 11, 44, -1, 44, 11, 10, -1, 32, 33, 36, 45, 46, 47, -1, 33, 34, 45, 48, -1, 35, 34, 38, 48, 49, 50, -1, 41, 43, 37, 42, 36, 46, -1, 38, 39, 49, 51, -1, 39, 40, 51, 42, 52, 46, 47, -1, 45, 47, 48, 52, 50, -1, 52, 51, 50, 49, + //polyhedron 2 + 6, 7, 19, 8, 22, 23, -1, 6, 35, 7, -1, 6, 35, 19, 38, -1, 35, 7, 38, 8, 39, 40, -1, 53, 22, 54, 19, 39, 38, -1, 23, 53, 8, 54, 40, -1, 53, 22, 23, -1, 39, 54, 40, + //polyhedron 3 + 35, 34, 38, 48, 49, 50, -1, 6, 35, 5, 34, 55, 56, -1, 6, 35, 19, 38, -1, 34, 56, 48, 57, 58, 59, -1, 60, 61, 49, 21, 38, 19, -1, 62, 50, 58, 48, -1, 60, 63, 49, 64, 50, 62, -1, 5, 6, 12, 19, 20, 21, -1, 55, 5, 65, 12, -1, 66, 67, 57, 65, 56, 55, -1, 63, 66, 64, 57, 59, -1, 64, 62, 59, 58, -1, 60, 63, 61, 66, 68, 67, -1, 61, 68, 21, 20, -1, 67, 68, 65, 20, 12}; + + + //{ 0,113,212,255,347 }; + + double barys [] = { -0.0165220465527,-0.0190922868195,0.158882733414, 0.0287618656076,0.135874379934,-0.14601588119, -0.147128055553,0.0465995097041,-0.049391174453, -0.00142506732317,-0.0996953090351,-0.115159183132 };*/ + + + int connN[]={0,3,2,1, -1, 4,5,6,7, -1, 0,4,7,3, -1, 3,7,6,2, -1, 2,6,5,1, -1, 1,5,4,0}; + double coords[]={0.,0.,0., 1.,0.,0., 1.,1.,0., 0.,1.,0., 0.,0.,1., 1.,0.,1., 1.,1.,1., 0.,1.,1., 0.5, 0.5, 0.5}; + MEDCouplingUMesh *meshN=MEDCouplingUMesh::New(); + meshN->setName("ForBary"); + meshN->setMeshDimension(3); + meshN->allocateCells(4); + meshN->insertNextCell(INTERP_KERNEL::NORM_POLYHED,29,connN); + /*meshN->insertNextCell(INTERP_KERNEL::NORM_POLYHED,113,connN); + meshN->insertNextCell(INTERP_KERNEL::NORM_POLYHED,99,connN+113); + meshN->insertNextCell(INTERP_KERNEL::NORM_POLYHED,43,connN+212); + meshN->insertNextCell(INTERP_KERNEL::NORM_POLYHED,92,connN+255);*/ + meshN->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(9,3);//9,3 + std::copy(coords,coords+27,myCoords->getPointer()); + meshN->setCoords(myCoords); + myCoords->decrRef(); + meshN->checkCoherency(); + // + std::vector res1; + meshN->arePolyhedronsNotCorrectlyOriented(res1); + meshN->orientCorrectlyPolyhedrons(); + CPPUNIT_ASSERT(res1.empty()); + const double *ref,*daPtr; + DataArrayDouble *da=meshN->getBarycenterAndOwner(); + CPPUNIT_ASSERT_EQUAL(1,da->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(3,da->getNumberOfComponents()); + daPtr=da->getConstPointer(); + ref=meshN->getCoords()->getConstPointer()+24; + for(int i=0;i<3;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(ref[i],daPtr[i],1e-12); + da->decrRef(); + // + const double center[]={0.,0.,0.}; + const double vec[]={0.,2.78,0.}; + da=meshN->getBarycenterAndOwner(); + daPtr=da->getConstPointer(); + ref=meshN->getCoords()->getConstPointer()+24; + for(int i=0;i<3;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(ref[i],daPtr[i],1e-12); + da->decrRef(); + // + meshN->rotate(center,vec,M_PI/7.); + meshN->translate(vec); + da=meshN->getBarycenterAndOwner(); + daPtr=da->getConstPointer(); + ref=meshN->getCoords()->getConstPointer()+24; + for(int i=0;i<3;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(ref[i],daPtr[i],1e-12); + da->decrRef(); + // + const double center2[]={1.12,3.45,6.78}; + const double vec2[]={4.5,9.3,2.8}; + meshN->rotate(center2,vec2,M_E); + meshN->translate(vec2); + da=meshN->getBarycenterAndOwner(); + daPtr=da->getConstPointer(); + ref=meshN->getCoords()->getConstPointer()+24; + for(int i=0;i<3;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(ref[i],daPtr[i],1e-5); + da->decrRef(); + // + meshN->decrRef(); +} + +void MEDCouplingBasicsTest::testNormL12Integ1D() +{ + MEDCouplingUMesh *m1=build1DTargetMesh_3(); + MEDCouplingFieldDouble *f1=MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME); + f1->setMesh(m1); + DataArrayDouble *array=DataArrayDouble::New(); + array->alloc(m1->getNumberOfCells(),3); + const double arr[12]={-5.23,15.45,-25.56,6.67,-16.78,26.89,-7.91,17.23,-27.43,8.21,-18.63,28.72}; + std::copy(arr,arr+12,array->getPointer()); + f1->setArray(array); + array->decrRef(); + // + const double *ptr; + DataArrayDouble *f3=m1->getBarycenterAndOwner(); + CPPUNIT_ASSERT_EQUAL(4,f3->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(1,f3->getNumberOfComponents()); + double expected9[4]={0.75,5.105,0.8,5.155}; + ptr=f3->getConstPointer(); + for(int i=0;i<4;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected9[i],ptr[i],1e-12); + f3->decrRef(); + // + MEDCouplingFieldDouble *f2=m1->getMeasureField(false); + CPPUNIT_ASSERT_EQUAL(4,f2->getArray()->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(1,f2->getNumberOfComponents()); + double expected1[4]={0.5,0.21,-0.6,-0.31}; + ptr=f2->getArray()->getConstPointer(); + for(int i=0;i<4;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected1[i],ptr[i],1e-12); + f2->decrRef(); + double expected2[4]={0.5,0.21,0.6,0.31}; + f2=m1->getMeasureField(true); + ptr=f2->getArray()->getConstPointer(); + for(int i=0;i<4;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected2[i],ptr[i],1e-12); + f2->decrRef(); + //integral + double res[3]; + f1->integral(false,res); + double expected3[3]={0.9866,-0.3615,0.4217}; + for(int i=0;i<3;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected3[i],res[i],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected3[0],f1->integral(0,false),1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected3[1],f1->integral(1,false),1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected3[2],f1->integral(2,false),1e-12); + f1->integral(true,res); + double expected4[3]={-3.4152,8.7639,-14.6879}; + for(int i=0;i<3;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected4[i],res[i],1e-12); + //normL1 + f1->normL1(false,res); + double expected5[3]={16.377,24.3225,34.6715}; + for(int i=0;i<3;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected5[i],res[i],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected5[0],f1->normL1(0,false),1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected5[1],f1->normL1(1,false),1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected5[2],f1->normL1(2,false),1e-12); + f1->normL1(true,res); + double expected6[3]={6.97950617284,16.8901851851852,27.0296913580247}; + for(int i=0;i<3;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected6[i],res[i],1e-12); + //normL2 + f1->normL2(false,res); + double expected7[3]={13.3073310622,23.1556650736,33.8113075021}; + for(int i=0;i<3;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected7[i],res[i],1e-9); + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected7[0],f1->normL2(0,false),1e-9); + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected7[1],f1->normL2(1,false),1e-9); + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected7[2],f1->normL2(2,false),1e-9); + f1->normL2(true,res); + double expected8[3]={7.09091097945239,16.9275542960123,27.0532714641609}; + for(int i=0;i<3;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected8[i],res[i],1e-9); + // + f1->decrRef(); + m1->decrRef(); + // Testing with 2D Curve + m1=build2DCurveTargetMesh_3(); + f2=m1->getMeasureField(false); + CPPUNIT_ASSERT_EQUAL(4,f2->getArray()->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(1,f2->getNumberOfComponents()); + ptr=f2->getArray()->getConstPointer(); + for(int i=0;i<4;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(sqrt(2.)*expected2[i],ptr[i],1e-12); + f2->decrRef(); + f2=m1->getMeasureField(true); + CPPUNIT_ASSERT_EQUAL(4,f2->getArray()->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(1,f2->getNumberOfComponents()); + ptr=f2->getArray()->getConstPointer(); + for(int i=0;i<4;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected2[i]*sqrt(2.),ptr[i],1e-12); + f2->decrRef(); + //bary + f3=m1->getBarycenterAndOwner(); + CPPUNIT_ASSERT_EQUAL(4,f3->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(2,f3->getNumberOfComponents()); + double expected10[8]={0.75,0.75,5.105,5.105,0.8,0.8,5.155,5.155}; + ptr=f3->getConstPointer(); + for(int i=0;i<8;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected10[i],ptr[i],1e-12); + f3->decrRef(); + // + f1=MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME); + f1->setMesh(m1); + array=DataArrayDouble::New(); + array->alloc(m1->getNumberOfCells(),3); + std::copy(arr,arr+12,array->getPointer()); + f1->setArray(array); + array->decrRef(); + f1->integral(false,res); + for(int i=0;i<3;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(sqrt(2.)*expected4[i],res[i],1e-12); + f1->integral(true,res); + for(int i=0;i<3;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(sqrt(2.)*expected4[i],res[i],1e-12); + f1->normL1(false,res); + for(int i=0;i<3;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected6[i],res[i],1e-12); + f1->normL1(true,res); + for(int i=0;i<3;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected6[i],res[i],1e-12); + f1->normL2(false,res); + for(int i=0;i<3;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected8[i],res[i],1e-12); + f1->normL2(true,res); + for(int i=0;i<3;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected8[i],res[i],1e-12); + // + f1->decrRef(); + m1->decrRef(); +} + +void MEDCouplingBasicsTest::testAreaBary2D() +{ + MEDCouplingUMesh *m1=build2DTargetMesh_3(); + MEDCouplingFieldDouble *f1=m1->getMeasureField(false); + CPPUNIT_ASSERT_EQUAL(10,f1->getArray()->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(1,f1->getNumberOfComponents()); + double expected1[10]={-0.5,-1,-1.5,-0.5,-1, 0.5,1,1.5,0.5,1}; + const double *ptr=f1->getArray()->getConstPointer(); + for(int i=0;i<10;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected1[i],ptr[i],1e-12); + f1->decrRef(); + f1=m1->getMeasureField(true); + ptr=f1->getArray()->getConstPointer(); + for(int i=0;i<10;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(std::abs(expected1[i]),ptr[i],1e-12); + f1->decrRef(); + DataArrayDouble *f2=m1->getBarycenterAndOwner(); + CPPUNIT_ASSERT_EQUAL(10,f2->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(2,f2->getNumberOfComponents()); + double expected2[20]={ + 0.5,0.3333333333333333,0.5,0.5,0.5,0.77777777777777777,0.5,0.3333333333333333,0.5,0.5, + 0.5,0.3333333333333333,0.5,0.5,0.5,0.77777777777777777,0.5,0.3333333333333333,0.5,0.5, + }; + ptr=f2->getConstPointer(); + for(int i=0;i<20;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected2[i],ptr[i],1e-12); + f2->decrRef(); + m1->changeSpaceDimension(3); + f1=m1->getMeasureField(false); + CPPUNIT_ASSERT_EQUAL(10,f1->getArray()->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(1,f1->getNumberOfComponents()); + ptr=f1->getArray()->getConstPointer(); + for(int i=0;i<10;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(std::abs(expected1[i]),ptr[i],1e-12); + f1->decrRef(); + f2=m1->getBarycenterAndOwner(); + CPPUNIT_ASSERT_EQUAL(10,f2->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(3,f2->getNumberOfComponents()); + ptr=f2->getConstPointer(); + double expected3[30]={ + 0.5,0.3333333333333333,0.,0.5,0.5,0.,0.5,0.77777777777777777,0.,0.5,0.3333333333333333,0.,0.5,0.5,0., + 0.5,0.3333333333333333,0.,0.5,0.5,0.,0.5,0.77777777777777777,0.,0.5,0.3333333333333333,0.,0.5,0.5,0. + }; + for(int i=0;i<30;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected3[i],ptr[i],1e-12); + f2->decrRef(); + m1->decrRef(); +} + +void MEDCouplingBasicsTest::testAreaBary3D() +{ +} + +void MEDCouplingBasicsTest::testRenumberCellsForFields() +{ + MEDCouplingUMesh *m=build2DTargetMesh_1(); + MEDCouplingFieldDouble *f=MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME); + f->setMesh(m); + DataArrayDouble *arr=DataArrayDouble::New(); + int nbOfCells=m->getNumberOfCells(); + arr->alloc(nbOfCells,3); + f->setArray(arr); + arr->decrRef(); + const double values1[15]={7.,107.,10007.,8.,108.,10008.,9.,109.,10009.,10.,110.,10010.,11.,111.,10011.}; + std::copy(values1,values1+15,arr->getPointer()); + const int renumber1[5]={3,1,0,4,2}; + double res[3]; + const double loc[]={-0.05,-0.05, 0.55,-0.25, 0.55,0.15, -0.05,0.45, 0.45,0.45}; + for(int j=0;j<5;j++) + { + f->getValueOn(loc+2*j,res); + for(int i=0;i<3;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(values1[i+3*j],res[i],1e-12); + } + f->renumberCells(renumber1,renumber1+5,false); + const double *ptr=f->getArray()->getConstPointer(); + const double expected1[15]={9.,109.,10009.,8.,108.,10008.,11.,111.,10011.,7.,107.,10007.,10.,110.,10010.}; + for(int i=0;i<15;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected1[i],ptr[i],1e-12); + //check that fields remains the same geometrically + for(int j=0;j<5;j++) + { + f->getValueOn(loc+2*j,res); + for(int i=0;i<3;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(values1[i+3*j],res[i],1e-12); + } + f->decrRef(); + //On gauss + f=MEDCouplingFieldDouble::New(ON_GAUSS_PT,NO_TIME); + f->setMesh(m); + const double _a=0.446948490915965; + const double _b=0.091576213509771; + const double _p1=0.11169079483905; + const double _p2=0.0549758718227661; + const double refCoo1[6]={ 0.,0., 1.,0., 0.,1. }; + const double gsCoo1[12]={ 2*_b-1, 1-4*_b, 2*_b-1, 2.07*_b-1, 1-4*_b, + 2*_b-1, 1-4*_a, 2*_a-1, 2*_a-1, 1-4*_a, 2*_a-1, 2*_a-1 }; + const double wg1[6]={ 4*_p2, 4*_p2, 4*_p2, 4*_p1, 4*_p1, 4*_p1 }; + std::vector _refCoo1(refCoo1,refCoo1+6); + std::vector _gsCoo1(gsCoo1,gsCoo1+12); + std::vector _wg1(wg1,wg1+6); + f->setGaussLocalizationOnType(INTERP_KERNEL::NORM_TRI3,_refCoo1,_gsCoo1,_wg1); + const double refCoo2[8]={ 0.,0., 1.,0., 1.,1., 0.,1. }; + std::vector _refCoo2(refCoo2,refCoo2+8); + _gsCoo1.resize(4); _wg1.resize(2); + f->setGaussLocalizationOnType(INTERP_KERNEL::NORM_QUAD4,_refCoo2,_gsCoo1,_wg1); + arr=DataArrayDouble::New(); + arr->alloc(18,2); + const double values2[36]={1.,1001.,2.,1002., 11.,1011.,12.,1012.,13.,1013.,14.,1014.,15.,1015.,16.,1016., 21.,1021.,22.,1022.,23.,1023.,24.,1024.,25.,1025.,26.,1026., 31.,1031.,32.,1032., 41.,1041.,42.,1042.}; + std::copy(values2,values2+36,arr->getPointer()); + f->setArray(arr); + arr->decrRef(); + f->checkCoherency(); + MEDCouplingFieldDouble *fCpy=f->clone(true); + CPPUNIT_ASSERT(f->isEqual(fCpy,1e-12,1e-12)); + f->renumberCells(renumber1,renumber1+5,false); + CPPUNIT_ASSERT(!f->isEqual(fCpy,1e-12,1e-12)); + double expected2[36]={21.,1021.,22.,1022.,23.,1023.,24.,1024.,25.,1025.,26.,1026., 11.,1011.,12.,1012.,13.,1013.,14.,1014.,15.,1015.,16.,1016., 41.,1041.,42.,1042., 1.,1001.,2.,1002., 31.,1031.,32.,1032.}; + ptr=f->getArray()->getConstPointer(); + for(int i=0;i<36;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected2[i],ptr[i],1e-12); + const int renumber2[5]={2,1,4,0,3};//reverse renumber1 + f->renumberCells(renumber2,renumber2+5,false); + CPPUNIT_ASSERT(f->isEqual(fCpy,1e-12,1e-12)); + fCpy->decrRef(); + f->decrRef(); + //GaussNE + f=MEDCouplingFieldDouble::New(ON_GAUSS_NE,NO_TIME); + f->setMesh(m); + arr=DataArrayDouble::New(); + arr->alloc(18,2); + const double values3[36]={1.,1001.,2.,1002.,3.,1003.,4.,1004., 11.,1011.,12.,1012.,13.,1013., 21.,1021.,22.,1022.,23.,1023., 31.,1031.,32.,1032.,33.,1033.,34.,1034., 41.,1041.,42.,1042.,43.,1043.,44.,1044.}; + std::copy(values3,values3+36,arr->getPointer()); + f->setArray(arr); + arr->decrRef(); + f->checkCoherency(); + fCpy=f->clone(true); + CPPUNIT_ASSERT(f->isEqual(fCpy,1e-12,1e-12)); + f->renumberCells(renumber1,renumber1+5,false); + CPPUNIT_ASSERT(!f->isEqual(fCpy,1e-12,1e-12)); + double expected3[36]={21.,1021.,22.,1022.,23.,1023.,11.,1011.,12.,1012.,13.,1013.,41.,1041.,42.,1042.,43.,1043.,44.,1044.,1.,1001.,2.,1002.,3.,1003.,4.,1004.,31.,1031.,32.,1032.,33.,1033.,34.,1034.}; + ptr=f->getArray()->getConstPointer(); + for(int i=0;i<36;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected3[i],ptr[i],1e-12); + f->renumberCells(renumber2,renumber2+5,false);//perform reverse operation of renumbering to check that the resulting field is equal. + CPPUNIT_ASSERT(f->isEqual(fCpy,1e-12,1e-12)); + fCpy->decrRef(); + f->decrRef(); + // + m->decrRef(); +} + +void MEDCouplingBasicsTest::testRenumberNodesForFields() +{ + MEDCouplingUMesh *m=build2DTargetMesh_1(); + MEDCouplingFieldDouble *f=MEDCouplingFieldDouble::New(ON_NODES,NO_TIME); + f->setMesh(m); + DataArrayDouble *arr=DataArrayDouble::New(); + int nbOfNodes=m->getNumberOfNodes(); + arr->alloc(nbOfNodes,3); + f->setArray(arr); + arr->decrRef(); + const double values1[27]={7.,107.,10007.,8.,108.,10008.,9.,109.,10009.,10.,110.,10010.,11.,111.,10011.,12.,112.,10012.,13.,113.,10013.,14.,114.,10014.,15.,115.,10015.}; + std::copy(values1,values1+27,arr->getPointer()); + f->checkCoherency(); + const int renumber1[9]={0,4,1,3,5,2,6,7,8}; + double res[3]; + const double loc[]={0.5432,-0.2432, 0.5478,0.1528}; + const double expected1[6]={9.0272, 109.0272, 10009.0272, 11.4124,111.4124,10011.4124}; + for(int j=0;j<2;j++) + { + f->getValueOn(loc+2*j,res); + for(int i=0;i<3;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected1[i+3*j],res[i],1e-12); + } + MEDCouplingFieldDouble *fCpy=f->clone(true); + CPPUNIT_ASSERT(f->isEqual(fCpy,1e-12,1e-12)); + f->renumberNodes(renumber1,renumber1+9); + CPPUNIT_ASSERT(!f->isEqual(fCpy,1e-12,1e-12)); + for(int j=0;j<2;j++) + { + f->getValueOn(loc+2*j,res); + for(int i=0;i<3;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected1[i+3*j],res[i],1e-12); + } + const double expected2[27]={7.,107.,10007.,9.,109.,10009.,12.,112.,10012.,10.,110.,10010.,8.,108.,10008.,11.,111.,10011.,13.,113.,10013.,14.,114.,10014.,15.,115.,10015.}; + for(int i=0;i<27;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected2[i],f->getArray()->getConstPointer()[i],1e-12); + const int renumber2[9]={0,2,5,3,1,4,6,7,8};//reverse of renumber2 + f->renumberNodes(renumber2,renumber2+9); + CPPUNIT_ASSERT(f->isEqual(fCpy,1e-12,1e-12)); + fCpy->decrRef(); + // + m->decrRef(); + f->decrRef(); +} diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py index 9a06753e0..e57506296 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py @@ -20,7 +20,7 @@ from libMEDCoupling_Swig import * import unittest -from math import pi +from math import pi,e,sqrt from MEDCouplingDataForTest import MEDCouplingDataForTest class MEDCouplingBasicsTest(unittest.TestCase): @@ -762,7 +762,7 @@ class MEDCouplingBasicsTest(unittest.TestCase): arr=f.getArray(); arrPtr=arr.getValues(); for i in xrange(15): - self.assertAlmostEqual(expected2[rexpected1[i]],arrPtr[i],1e-16); + self.assertAlmostEqual(expected2[rexpected1[i]],arrPtr[i],15); pass pass @@ -980,7 +980,7 @@ class MEDCouplingBasicsTest(unittest.TestCase): values4=f1.accumulate(); self.assertTrue(abs(3.6-values4[0])<1.e-12); self.assertTrue(abs(7.2-values4[1])<1.e-12); - values4=f1.measureAccumulate(True); + values4=f1.integral(True); self.assertTrue(abs(0.5-values4[0])<1.e-12); self.assertTrue(abs(1.-values4[1])<1.e-12); # @@ -1042,7 +1042,7 @@ class MEDCouplingBasicsTest(unittest.TestCase): values4=f1.accumulate(); self.assertTrue(abs(3.6-values4[0])<1.e-12); self.assertTrue(abs(7.2-values4[1])<1.e-12); - values4=f1.measureAccumulate(True); + values4=f1.integral(True); self.assertTrue(abs(0.5-values4[0])<1.e-12); self.assertTrue(abs(1.-values4[1])<1.e-12); pass @@ -1870,8 +1870,379 @@ class MEDCouplingBasicsTest(unittest.TestCase): m5.orientCorrectlyPolyhedrons(); res1=m5.arePolyhedronsNotCorrectlyOriented(); self.assertTrue(len(res1)==0); + f3=m5.getMeasureField(False); + self.assertEqual(15,f3.getArray().getNumberOfTuples()); + self.assertEqual(1,f3.getNumberOfComponents()); + f3Ptr=f3.getArray().getValues(); + expected1=[0.075,0.0375,0.0375,0.075,0.075, 0.1125,0.05625,0.05625,0.1125,0.1125, 0.0625,0.03125,0.03125,0.0625,0.0625]; + for i in xrange(15): + self.assertTrue(abs(expected1[i]-f3Ptr[i])<1e-12); + pass + f4=m5.getBarycenterAndOwner(); + self.assertEqual(15,f4.getNumberOfTuples()); + self.assertEqual(3,f4.getNumberOfComponents()); + f4Ptr=f4.getValues(); + expected2=[-0.05,-0.05,0.15, 0.3666666666666667,-0.13333333333333333,0.15, 0.53333333333333333,0.033333333333333333,0.15, -0.05,0.45,0.15, 0.45,0.45,0.15,-0.05,-0.05,0.525, 0.3666666666666667,-0.13333333333333333,0.525, 0.53333333333333333,0.033333333333333333,0.525, -0.05,0.45,0.525, 0.45,0.45,0.525,-0.05,-0.05,0.875, 0.3666666666666667,-0.13333333333333333,0.875, 0.53333333333333333,0.033333333333333333,0.875, -0.05,0.45,0.875, 0.45,0.45,0.875]; + for i in xrange(45): + self.assertTrue(abs(expected2[i]-f4Ptr[i])<1e-12); + pass + pass + + def testPolyhedronBarycenter(self): + connN=[0,3,2,1, -1, 4,5,6,7, -1, 0,4,7,3, -1, 3,7,6,2, -1, 2,6,5,1, -1, 1,5,4,0]; + coords=[0.,0.,0., 1.,0.,0., 1.,1.,0., 0.,1.,0., 0.,0.,1., 1.,0.,1., 1.,1.,1., 0.,1.,1., 0.5, 0.5, 0.5]; + meshN=MEDCouplingUMesh.New(); + meshN.setName("ForBary"); + meshN.setMeshDimension(3); + meshN.allocateCells(4); + meshN.insertNextCell(NORM_POLYHED,29,connN[0:29]) + meshN.finishInsertingCells(); + myCoords=DataArrayDouble.New(); + myCoords.setValues(coords,9,3); + meshN.setCoords(myCoords); + meshN.checkCoherency(); + # + res1=meshN.arePolyhedronsNotCorrectlyOriented(); + meshN.orientCorrectlyPolyhedrons(); + self.assertTrue(len(res1)==0); + da=meshN.getBarycenterAndOwner(); + self.assertEqual(1,da.getNumberOfTuples()); + self.assertEqual(3,da.getNumberOfComponents()); + daPtr=da.getValues(); + ref=meshN.getCoords().getValues()[24:]; + for i in xrange(3): + self.assertTrue(abs(ref[i]-daPtr[i])<1e-12); + pass + # + center=[0.,0.,0.] + vec=[0.,2.78,0.] + da=meshN.getBarycenterAndOwner(); + daPtr=da.getValues(); + ref=meshN.getCoords().getValues()[24:]; + for i in xrange(3): + self.assertTrue(abs(ref[i]-daPtr[i])<1e-12); + pass + # + meshN.rotate(center,vec,pi/7.); + meshN.translate(vec); + da=meshN.getBarycenterAndOwner(); + daPtr=da.getValues(); + ref=meshN.getCoords().getValues()[24:]; + for i in xrange(3): + self.assertTrue(abs(ref[i]-daPtr[i])<1e-12); + pass + # + center2=[1.12,3.45,6.78] + vec2=[4.5,9.3,2.8] + meshN.rotate(center2,vec2,e); + meshN.translate(vec2); + da=meshN.getBarycenterAndOwner(); + daPtr=da.getValues(); + ref=meshN.getCoords().getValues()[24:]; + for i in xrange(3): + self.assertTrue(abs(ref[i]-daPtr[i])<1e-5); + pass + pass + + def testNormL12Integ1D(self): + m1=MEDCouplingDataForTest.build1DTargetMesh_3(); + f1=MEDCouplingFieldDouble.New(ON_CELLS,NO_TIME); + f1.setMesh(m1); + array=DataArrayDouble.New(); + arr=[-5.23,15.45,-25.56,6.67,-16.78,26.89,-7.91,17.23,-27.43,8.21,-18.63,28.72] + array.setValues(arr,m1.getNumberOfCells(),3); + f1.setArray(array); + # + f3=m1.getBarycenterAndOwner(); + self.assertEqual(4,f3.getNumberOfTuples()); + self.assertEqual(1,f3.getNumberOfComponents()); + expected9=[0.75,5.105,0.8,5.155] + ptr=f3.getValues(); + for i in xrange(4): + self.assertTrue(abs(expected9[i]-ptr[i])<1e-12); + pass + # + f2=m1.getMeasureField(False); + self.assertEqual(4,f2.getArray().getNumberOfTuples()); + self.assertEqual(1,f2.getNumberOfComponents()); + expected1=[0.5,0.21,-0.6,-0.31] + ptr=f2.getArray().getValues(); + for i in xrange(4): + self.assertTrue(abs(expected1[i]-ptr[i])<1e-12); + pass + expected2=[0.5,0.21,0.6,0.31] + f2=m1.getMeasureField(True); + ptr=f2.getArray().getValues(); + for i in xrange(4): + self.assertTrue(abs(expected2[i]-ptr[i])<1e-12); + pass + #integral + res=f1.integral(False); + expected3=[0.9866,-0.3615,0.4217] + for i in xrange(3): + self.assertTrue(abs(expected3[i]-res[i])<1e-12); + pass + self.assertTrue(abs(expected3[0]-f1.integral(0,False))<1e-12); + self.assertTrue(abs(expected3[1]-f1.integral(1,False))<1e-12); + self.assertTrue(abs(expected3[2]-f1.integral(2,False))<1e-12); + res=f1.integral(True); + expected4=[-3.4152,8.7639,-14.6879] + for i in xrange(3): + self.assertTrue(abs(expected4[i]-res[i])<1e-12); + pass + #normL1 + res=f1.normL1(False); + expected5=[16.377,24.3225,34.6715] + for i in xrange(3): + self.assertTrue(abs(expected5[i]-res[i])<1e-12); + pass + self.assertTrue(abs(expected5[0]-f1.normL1(0,False))<1e-12); + self.assertTrue(abs(expected5[1]-f1.normL1(1,False))<1e-12); + self.assertTrue(abs(expected5[2]-f1.normL1(2,False))<1e-12); + res=f1.normL1(True); + expected6=[6.97950617284,16.8901851851852,27.0296913580247] + for i in xrange(3): + self.assertTrue(abs(expected6[i]-res[i])<1e-12); + pass + #normL2 + res=f1.normL2(False); + expected7=[13.3073310622,23.1556650736,33.8113075021] + for i in xrange(3): + self.assertTrue(abs(expected7[i]-res[i])<1e-9); + pass + self.assertTrue(abs(expected7[0]-f1.normL2(0,False))<1e-9); + self.assertTrue(abs(expected7[1]-f1.normL2(1,False))<1e-9); + self.assertTrue(abs(expected7[2]-f1.normL2(2,False))<1e-9); + res=f1.normL2(True); + expected8=[7.09091097945239,16.9275542960123,27.0532714641609] + for i in xrange(3): + self.assertTrue(abs(expected8[i]-res[i])<1e-9); + pass + # + # Testing with 2D Curve + m1=MEDCouplingDataForTest.build2DCurveTargetMesh_3(); + f2=m1.getMeasureField(False); + self.assertEqual(4,f2.getArray().getNumberOfTuples()); + self.assertEqual(1,f2.getNumberOfComponents()); + ptr=f2.getArray().getValues(); + for i in xrange(4): + self.assertTrue(abs(sqrt(2.)*expected2[i]-ptr[i])<1e-12); + pass + f2=m1.getMeasureField(True); + self.assertEqual(4,f2.getArray().getNumberOfTuples()); + self.assertEqual(1,f2.getNumberOfComponents()); + ptr=f2.getArray().getValues(); + for i in xrange(4): + self.assertTrue(abs(expected2[i]*sqrt(2.)-ptr[i])<1e-12); + pass + #bary + f3=m1.getBarycenterAndOwner(); + self.assertEqual(4,f3.getNumberOfTuples()); + self.assertEqual(2,f3.getNumberOfComponents()); + expected10=[0.75,0.75,5.105,5.105,0.8,0.8,5.155,5.155] + ptr=f3.getValues(); + for i in xrange(8): + self.assertTrue(abs(expected10[i]-ptr[i])<1e-12); + pass + # + f1=MEDCouplingFieldDouble.New(ON_CELLS,NO_TIME); + f1.setMesh(m1); + array=DataArrayDouble.New(); + array.setValues(arr,m1.getNumberOfCells(),3); + f1.setArray(array); + res=f1.integral(False); + for i in xrange(3): + self.assertTrue(abs(sqrt(2.)*expected4[i]-res[i])<1e-12); + pass + res=f1.integral(True); + for i in xrange(3): + self.assertTrue(abs(sqrt(2.)*expected4[i]-res[i])<1e-12); + pass + res=f1.normL1(False); + for i in xrange(3): + self.assertTrue(abs(expected6[i]-res[i])<1e-12); + pass + res=f1.normL1(True); + for i in xrange(3): + self.assertTrue(abs(expected6[i]-res[i])<1e-12); + pass + res=f1.normL2(False); + for i in xrange(3): + self.assertTrue(abs(expected8[i]-res[i])<1e-12); + pass + res=f1.normL2(True); + for i in xrange(3): + self.assertTrue(abs(expected8[i]-res[i])<1e-12); + pass pass + def testAreaBary2D(self): + m1=MEDCouplingDataForTest.build2DTargetMesh_3(); + f1=m1.getMeasureField(False); + self.assertEqual(10,f1.getArray().getNumberOfTuples()); + self.assertEqual(1,f1.getNumberOfComponents()); + expected1=[-0.5,-1,-1.5,-0.5,-1, 0.5,1,1.5,0.5,1] + ptr=f1.getArray().getValues(); + for i in xrange(10): + self.assertTrue(abs(expected1[i]-ptr[i])<1e-12); + pass + f1=m1.getMeasureField(True); + ptr=f1.getArray().getValues(); + for i in xrange(10): + self.assertTrue(abs(abs(expected1[i])-ptr[i])<1e-12); + pass + f2=m1.getBarycenterAndOwner(); + self.assertEqual(10,f2.getNumberOfTuples()); + self.assertEqual(2,f2.getNumberOfComponents()); + expected2=[0.5,0.3333333333333333,0.5,0.5,0.5,0.77777777777777777,0.5,0.3333333333333333,0.5,0.5,0.5,0.3333333333333333,0.5,0.5,0.5,0.77777777777777777,0.5,0.3333333333333333,0.5,0.5] + ptr=f2.getValues(); + for i in xrange(20): + self.assertTrue(abs(expected2[i]-ptr[i])<1e-12); + pass + m1.changeSpaceDimension(3); + f1=m1.getMeasureField(False); + self.assertEqual(10,f1.getArray().getNumberOfTuples()); + self.assertEqual(1,f1.getNumberOfComponents()); + ptr=f1.getArray().getValues(); + for i in xrange(10): + self.assertTrue(abs(abs(expected1[i])-ptr[i])<1e-12); + pass + f2=m1.getBarycenterAndOwner(); + self.assertEqual(10,f2.getNumberOfTuples()); + self.assertEqual(3,f2.getNumberOfComponents()); + ptr=f2.getValues(); + expected3=[0.5,0.3333333333333333,0.,0.5,0.5,0.,0.5,0.77777777777777777,0.,0.5,0.3333333333333333,0.,0.5,0.5,0., 0.5,0.3333333333333333,0.,0.5,0.5,0.,0.5,0.77777777777777777,0.,0.5,0.3333333333333333,0.,0.5,0.5,0.] + for i in xrange(30): + self.assertTrue(abs(expected3[i]-ptr[i])<1e-12); + pass + pass + + def testRenumberCellsForFields(self): + m=MEDCouplingDataForTest.build2DTargetMesh_1(); + f=MEDCouplingFieldDouble.New(ON_CELLS,NO_TIME); + f.setMesh(m); + arr=DataArrayDouble.New(); + nbOfCells=m.getNumberOfCells(); + values1=[7.,107.,10007.,8.,108.,10008.,9.,109.,10009.,10.,110.,10010.,11.,111.,10011.] + arr.setValues(values1,nbOfCells,3); + f.setArray(arr); + renumber1=[3,1,0,4,2] + loc=[-0.05,-0.05, 0.55,-0.25, 0.55,0.15, -0.05,0.45, 0.45,0.45] + for j in xrange(5): + res=f.getValueOn(loc[2*j:2*j+2]); + for i in xrange(3): + self.assertTrue(abs(values1[i+3*j]-res[i])<1e-12); + pass + pass + f.renumberCells(renumber1,False); + ptr=f.getArray().getValues(); + expected1=[9.,109.,10009.,8.,108.,10008.,11.,111.,10011.,7.,107.,10007.,10.,110.,10010.] + for i in xrange(15): + self.assertTrue(abs(expected1[i]-ptr[i])<1e-12); + pass + #check that fields remains the same geometrically + for j in xrange(5): + res=f.getValueOn(loc[2*j:2*(j+1)]); + for i in xrange(3): + self.assertTrue(abs(values1[i+3*j]-res[i])<1e-12); + pass + pass + #On gauss + f=MEDCouplingFieldDouble.New(ON_GAUSS_PT,NO_TIME); + f.setMesh(m); + _a=0.446948490915965; + _b=0.091576213509771; + _p1=0.11169079483905; + _p2=0.0549758718227661; + refCoo1=[ 0.,0., 1.,0., 0.,1. ] + gsCoo1=[ 2*_b-1, 1-4*_b, 2*_b-1, 2.07*_b-1, 1-4*_b, 2*_b-1, 1-4*_a, 2*_a-1, 2*_a-1, 1-4*_a, 2*_a-1, 2*_a-1 ]; + wg1=[ 4*_p2, 4*_p2, 4*_p2, 4*_p1, 4*_p1, 4*_p1 ] + _refCoo1=refCoo1[0:6]; + _gsCoo1=gsCoo1[0:12]; + _wg1=wg1[0:6]; + f.setGaussLocalizationOnType(NORM_TRI3,_refCoo1,_gsCoo1,_wg1); + refCoo2=[ 0.,0., 1.,0., 1.,1., 0.,1. ] + _refCoo2=refCoo2[0:8]; + _gsCoo1=_gsCoo1[0:4] + _wg1=_wg1[0:2] + f.setGaussLocalizationOnType(NORM_QUAD4,_refCoo2,_gsCoo1,_wg1); + arr=DataArrayDouble.New(); + values2=[1.,1001.,2.,1002., 11.,1011.,12.,1012.,13.,1013.,14.,1014.,15.,1015.,16.,1016., 21.,1021.,22.,1022.,23.,1023.,24.,1024.,25.,1025.,26.,1026., 31.,1031.,32.,1032., 41.,1041.,42.,1042.] + arr.setValues(values2,18,2); + f.setArray(arr); + f.checkCoherency(); + fCpy=f.clone(True); + self.assertTrue(f.isEqual(fCpy,1e-12,1e-12)); + f.renumberCells(renumber1,False); + self.assertTrue(not f.isEqual(fCpy,1e-12,1e-12)); + expected2=[21.,1021.,22.,1022.,23.,1023.,24.,1024.,25.,1025.,26.,1026., 11.,1011.,12.,1012.,13.,1013.,14.,1014.,15.,1015.,16.,1016., 41.,1041.,42.,1042., 1.,1001.,2.,1002., 31.,1031.,32.,1032.] + ptr=f.getArray().getValues(); + for i in xrange(36): + self.assertTrue(abs(expected2[i]-ptr[i])<1e-12); + pass + renumber2=[2,1,4,0,3] + f.renumberCells(renumber2,False); + self.assertTrue(f.isEqual(fCpy,1e-12,1e-12)); + #GaussNE + f=MEDCouplingFieldDouble.New(ON_GAUSS_NE,NO_TIME); + f.setMesh(m); + arr=DataArrayDouble.New(); + values3=[1.,1001.,2.,1002.,3.,1003.,4.,1004., 11.,1011.,12.,1012.,13.,1013., 21.,1021.,22.,1022.,23.,1023., 31.,1031.,32.,1032.,33.,1033.,34.,1034., 41.,1041.,42.,1042.,43.,1043.,44.,1044.] + arr.setValues(values3,18,2); + f.setArray(arr); + f.checkCoherency(); + fCpy=f.clone(True); + self.assertTrue(f.isEqual(fCpy,1e-12,1e-12)); + f.renumberCells(renumber1,False); + self.assertTrue(not f.isEqual(fCpy,1e-12,1e-12)); + expected3=[21.,1021.,22.,1022.,23.,1023.,11.,1011.,12.,1012.,13.,1013.,41.,1041.,42.,1042.,43.,1043.,44.,1044.,1.,1001.,2.,1002.,3.,1003.,4.,1004.,31.,1031.,32.,1032.,33.,1033.,34.,1034.] + ptr=f.getArray().getValues(); + for i in xrange(36): + self.assertTrue(abs(expected3[i]-ptr[i])<1e-12); + pass + f.renumberCells(renumber2,False);#perform reverse operation of renumbering to check that the resulting field is equal. + self.assertTrue(f.isEqual(fCpy,1e-12,1e-12)); + # + pass + + def testRenumberNodesForFields(self): + m=MEDCouplingDataForTest.build2DTargetMesh_1(); + f=MEDCouplingFieldDouble.New(ON_NODES,NO_TIME); + f.setMesh(m); + arr=DataArrayDouble.New(); + nbOfNodes=m.getNumberOfNodes(); + values1=[7.,107.,10007.,8.,108.,10008.,9.,109.,10009.,10.,110.,10010.,11.,111.,10011.,12.,112.,10012.,13.,113.,10013.,14.,114.,10014.,15.,115.,10015.] + arr.setValues(values1,nbOfNodes,3); + f.setArray(arr); + f.checkCoherency(); + renumber1=[0,4,1,3,5,2,6,7,8] + loc=[0.5432,-0.2432, 0.5478,0.1528] + expected1=[9.0272, 109.0272, 10009.0272, 11.4124,111.4124,10011.4124] + for j in xrange(2): + res=f.getValueOn(loc[2*j:2*j+2]); + for i in xrange(3): + self.assertTrue(abs(expected1[i+3*j]-res[i])<1e-12); + pass + pass + fCpy=f.clone(True); + self.assertTrue(f.isEqual(fCpy,1e-12,1e-12)); + f.renumberNodes(renumber1); + self.assertTrue(not f.isEqual(fCpy,1e-12,1e-12)); + for j in xrange(2): + res=f.getValueOn(loc[2*j:2*j+2]); + for i in xrange(3): + self.assertTrue(abs(expected1[i+3*j]-res[i])<1e-12); + pass + pass + expected2=[7.,107.,10007.,9.,109.,10009.,12.,112.,10012.,10.,110.,10010.,8.,108.,10008.,11.,111.,10011.,13.,113.,10013.,14.,114.,10014.,15.,115.,10015.] + for i in xrange(27): + self.assertTrue(abs(expected2[i]-f.getArray().getValues()[i])<1e-12); + pass + renumber2=[0,2,5,3,1,4,6,7,8] + f.renumberNodes(renumber2); + self.assertTrue(f.isEqual(fCpy,1e-12,1e-12)); + pass + def setUp(self): pass pass diff --git a/src/MEDCoupling_Swig/MEDCouplingDataForTest.py b/src/MEDCoupling_Swig/MEDCouplingDataForTest.py index b9b8b0c23..b726e5618 100644 --- a/src/MEDCoupling_Swig/MEDCouplingDataForTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingDataForTest.py @@ -252,6 +252,58 @@ class MEDCouplingDataForTest: targetMesh.setCoords(myCoords); return targetMesh; + def build1DTargetMesh_3(cls): + ret=MEDCouplingUMesh.New("1DMesh_3",1); + ret.allocateCells(4); + conn=[0,1,2, 3,4, 6,5,7 ,9,8] + ret.insertNextCell(NORM_SEG3,3,conn[0:3]) + ret.insertNextCell(NORM_SEG2,2,conn[3:5]) + ret.insertNextCell(NORM_SEG3,3,conn[5:8]) + ret.insertNextCell(NORM_SEG2,2,conn[8:10]) + ret.finishInsertingCells(); + coords=[0.5,1.,0.8,5.,5.21,0.5,1.1,0.7,5.,5.31] + myCoords=DataArrayDouble.New(); + myCoords.setValues(coords,10,1); + ret.setCoords(myCoords); + return ret; + + def build2DCurveTargetMesh_3(cls): + ret=MEDCouplingUMesh.New("2DCurveMesh_3",1); + ret.allocateCells(4); + conn=[0,1,2, 3,4, 6,5,7 ,9,8] + ret.insertNextCell(NORM_SEG3,3,conn[0:3]) + ret.insertNextCell(NORM_SEG2,2,conn[3:5]) + ret.insertNextCell(NORM_SEG3,3,conn[5:8]) + ret.insertNextCell(NORM_SEG2,2,conn[8:10]) + ret.finishInsertingCells(); + coords=[0.5,0.5,1.,1.,0.8,0.8,5.,5.,5.21,5.21,0.5,0.5,1.1,1.1,0.7,0.7,5.,5.,5.31,5.31] + myCoords=DataArrayDouble.New(); + myCoords.setValues(coords,10,2); + ret.setCoords(myCoords); + return ret; + + def build2DTargetMesh_3(cls): + ret=MEDCouplingUMesh.New("2DMesh_3",2); + ret.allocateCells(10); + conn=[0,1,2, 0,1,3,4, 0,1,3,5,4, 0,1,2,6,7,8, 0,1,3,4,6,9,2,10, 0,2,1, 0,4,3,1, 0,4,5,3,1, 0,2,1,8,7,6, 0,4,3,1,10,2,9,6] + ret.insertNextCell(NORM_TRI3,3,conn[0:3]) + ret.insertNextCell(NORM_QUAD4,4,conn[3:7]) + ret.insertNextCell(NORM_POLYGON,5,conn[7:12]) + ret.insertNextCell(NORM_TRI6,6,conn[12:18]) + ret.insertNextCell(NORM_QUAD8,8,conn[18:26]) + ret.insertNextCell(NORM_TRI3,3,conn[26:29]) + ret.insertNextCell(NORM_QUAD4,4,conn[29:33]) + ret.insertNextCell(NORM_POLYGON,5,conn[33:38]) + ret.insertNextCell(NORM_TRI6,6,conn[38:44]) + ret.insertNextCell(NORM_QUAD8,8,conn[44:52]) + ret.finishInsertingCells(); + coords=[0.,0.,1.,0.,0.5,1.,1.,1.,0.,1.,0.5,2.,0.5,0.,0.75,0.5,0.25,0.5,1.,0.5,0.,0.5] + myCoords=DataArrayDouble.New(); + myCoords.setValues(coords,11,2); + ret.setCoords(myCoords); + ret.checkCoherency(); + return ret; + build2DTargetMesh_1=classmethod(build2DTargetMesh_1) build2DSourceMesh_1=classmethod(build2DSourceMesh_1) build3DTargetMesh_1=classmethod(build3DTargetMesh_1) @@ -262,4 +314,7 @@ class MEDCouplingDataForTest: build3DTargetMeshMergeNode_1=classmethod(build3DTargetMeshMergeNode_1) build2DTargetMeshMerged_1=classmethod(build2DTargetMeshMerged_1) build2DTargetMesh_2=classmethod(build2DTargetMesh_2) + build1DTargetMesh_3=classmethod(build1DTargetMesh_3) + build2DCurveTargetMesh_3=classmethod(build2DCurveTargetMesh_3) + build2DTargetMesh_3=classmethod(build2DTargetMesh_3) pass diff --git a/src/MEDCoupling_Swig/libMEDCoupling_Swig.i b/src/MEDCoupling_Swig/libMEDCoupling_Swig.i index 42397071b..94a49e7f1 100644 --- a/src/MEDCoupling_Swig/libMEDCoupling_Swig.i +++ b/src/MEDCoupling_Swig/libMEDCoupling_Swig.i @@ -76,6 +76,8 @@ using namespace INTERP_KERNEL; %newobject ParaMEDMEM::DataArrayDouble::substr; %newobject ParaMEDMEM::MEDCouplingFieldDouble::clone; %newobject ParaMEDMEM::MEDCouplingFieldDouble::buildNewTimeReprFromThis; +%newobject ParaMEDMEM::MEDCouplingMesh::getCoordinatesAndOwner; +%newobject ParaMEDMEM::MEDCouplingMesh::getBarycenterAndOwner; %newobject ParaMEDMEM::MEDCouplingMesh::buildOrthogonalField; %newobject ParaMEDMEM::MEDCouplingMesh::mergeMyselfWith; %newobject ParaMEDMEM::MEDCouplingMesh::fillFromAnalytic; @@ -169,8 +171,18 @@ namespace ParaMEDMEM virtual void rotate(const double *center, const double *vector, double angle) = 0; virtual void translate(const double *vector) = 0; virtual MEDCouplingMesh *mergeMyselfWith(const MEDCouplingMesh *other) const throw(INTERP_KERNEL::Exception) = 0; - virtual bool areCompatible(const MEDCouplingMesh *other) const; + virtual bool areCompatibleForMerge(const MEDCouplingMesh *other) const; static MEDCouplingMesh *mergeMeshes(const MEDCouplingMesh *mesh1, const MEDCouplingMesh *mesh2); + %extend + { + void renumberCells(PyObject *li, bool check) throw(INTERP_KERNEL::Exception) + { + int size; + int *tmp=convertPyToNewIntArr2(li,&size); + self->renumberCells(tmp,tmp+size,check); + delete [] tmp; + } + } }; } @@ -405,13 +417,6 @@ namespace ParaMEDMEM PyList_SetItem(res,1,SWIG_From_bool(ret1)); return res; } - void renumberCells(PyObject *li, bool check) - { - int size; - int *tmp=convertPyToNewIntArr2(li,&size); - self->renumberCells(tmp,tmp+size,check); - delete [] tmp; - } PyObject *checkButterflyCells() { std::vector cells; @@ -639,7 +644,7 @@ namespace ParaMEDMEM { public: virtual void checkCoherency() const throw(INTERP_KERNEL::Exception); - virtual bool areCompatible(const MEDCouplingField *other) const; + virtual bool areCompatibleForMerge(const MEDCouplingField *other) const; virtual bool isEqual(const MEDCouplingField *other, double meshPrec, double valsPrec) const; void setMesh(const ParaMEDMEM::MEDCouplingMesh *mesh); void setName(const char *name); @@ -727,7 +732,9 @@ namespace ParaMEDMEM void applyFunc(int nbOfComp, const char *func) throw(INTERP_KERNEL::Exception); void applyFunc(const char *func) throw(INTERP_KERNEL::Exception); double accumulate(int compId) const throw(INTERP_KERNEL::Exception); - double measureAccumulate(int compId, bool isWAbs) const throw(INTERP_KERNEL::Exception); + double integral(int compId, bool isWAbs) const throw(INTERP_KERNEL::Exception); + double normL1(int compId, bool isWAbs) const throw(INTERP_KERNEL::Exception); + double normL2(int compId, bool isWAbs) const throw(INTERP_KERNEL::Exception); static MEDCouplingFieldDouble *mergeFields(const MEDCouplingFieldDouble *f1, const MEDCouplingFieldDouble *f2) throw(INTERP_KERNEL::Exception); MEDCouplingFieldDouble *operator+(const MEDCouplingFieldDouble& other) const throw(INTERP_KERNEL::Exception); const MEDCouplingFieldDouble &operator+=(const MEDCouplingFieldDouble& other) throw(INTERP_KERNEL::Exception); @@ -834,15 +841,64 @@ namespace ParaMEDMEM delete [] tmp; return ret; } - PyObject *measureAccumulate(bool isWAbs) const + PyObject *integral(bool isWAbs) const + { + int sz=self->getNumberOfTuples(); + double *tmp=new double[sz]; + self->integral(isWAbs,tmp); + PyObject *ret=convertDblArrToPyList(tmp,sz); + delete [] tmp; + return ret; + } + PyObject *normL1(bool isWAbs) const throw(INTERP_KERNEL::Exception) + { + int sz=self->getNumberOfTuples(); + double *tmp=new double[sz]; + self->normL1(isWAbs,tmp); + PyObject *ret=convertDblArrToPyList(tmp,sz); + delete [] tmp; + return ret; + } + PyObject *normL2(bool isWAbs) const throw(INTERP_KERNEL::Exception) { int sz=self->getNumberOfTuples(); double *tmp=new double[sz]; - self->measureAccumulate(isWAbs,tmp); + self->normL2(isWAbs,tmp); PyObject *ret=convertDblArrToPyList(tmp,sz); delete [] tmp; return ret; } + + void renumberCells(PyObject *li, bool check) throw(INTERP_KERNEL::Exception) + { + int size; + int *tmp=convertPyToNewIntArr2(li,&size); + try + { + self->renumberCells(tmp,tmp+size,check); + } + catch(INTERP_KERNEL::Exception& e) + { + delete [] tmp; + throw e; + } + delete [] tmp; + } + void renumberNodes(PyObject *li) throw(INTERP_KERNEL::Exception) + { + int size; + int *tmp=convertPyToNewIntArr2(li,&size); + try + { + self->renumberNodes(tmp,tmp+size); + } + catch(INTERP_KERNEL::Exception& e) + { + delete [] tmp; + throw e; + } + delete [] tmp; + } } }; } diff --git a/src/ParaMEDMEM/ParaFIELD.cxx b/src/ParaMEDMEM/ParaFIELD.cxx index 6822a19fb..e0f0ba085 100644 --- a/src/ParaMEDMEM/ParaFIELD.cxx +++ b/src/ParaMEDMEM/ParaFIELD.cxx @@ -213,7 +213,7 @@ namespace ParaMEDMEM double ParaFIELD::getVolumeIntegral(int icomp, bool isWAbs) const { CommInterface comm_interface = _topology->getProcGroup()->getCommInterface(); - double integral=_field->measureAccumulate(icomp,isWAbs); + double integral=_field->integral(icomp,isWAbs); double total=0.; const MPI_Comm* comm = (dynamic_cast(_topology->getProcGroup()))->getComm(); comm_interface.allReduce(&integral, &total, 1, MPI_DOUBLE, MPI_SUM, *comm); -- 2.39.2