From: ageay Date: Fri, 3 Dec 2010 17:59:11 +0000 (+0000) Subject: *** empty log message *** X-Git-Tag: V6_main_FINAL~1144 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=9fb9cc885756716843012aa3182bd4e515304937;p=tools%2Fmedcoupling.git *** empty log message *** --- diff --git a/src/INTERP_KERNEL/CellModel.cxx b/src/INTERP_KERNEL/CellModel.cxx index 88bc7a820..add613149 100644 --- a/src/INTERP_KERNEL/CellModel.cxx +++ b/src/INTERP_KERNEL/CellModel.cxx @@ -108,22 +108,22 @@ namespace INTERP_KERNEL { case NORM_POINT0: { - _nb_of_pts=0; _nb_of_sons=0; _dim=0; _extruded_type=NORM_SEG2; + _nb_of_pts=0; _nb_of_sons=0; _dim=0; _extruded_type=NORM_SEG2; _is_simplex=true; } break; case NORM_SEG2: { - _nb_of_pts=2; _nb_of_sons=0; _dim=1; _extruded_type=NORM_QUAD4; _quadratic_type=NORM_SEG3; + _nb_of_pts=2; _nb_of_sons=0; _dim=1; _extruded_type=NORM_QUAD4; _quadratic_type=NORM_SEG3; _is_simplex=true; } break; case NORM_SEG3: { - _nb_of_pts=3; _nb_of_sons=0; _dim=1; _extruded_type=NORM_QUAD8; _linear_type=NORM_SEG2; + _nb_of_pts=3; _nb_of_sons=0; _dim=1; _extruded_type=NORM_QUAD8; _linear_type=NORM_SEG2; _is_simplex=false; } break; case NORM_TETRA4: { - _nb_of_pts=4; _nb_of_sons=4; _dim=3; _quadratic_type=NORM_TETRA10; + _nb_of_pts=4; _nb_of_sons=4; _dim=3; _quadratic_type=NORM_TETRA10; _is_simplex=true; _sons_type[0]=NORM_TRI3; _sons_type[1]=NORM_TRI3; _sons_type[2]=NORM_TRI3; _sons_type[3]=NORM_TRI3; _sons_con[0][0]=0; _sons_con[0][1]=1; _sons_con[0][2]=2; _nb_of_sons_con[0]=3; _sons_con[1][0]=0; _sons_con[1][1]=3; _sons_con[1][2]=1; _nb_of_sons_con[1]=3; @@ -133,7 +133,7 @@ namespace INTERP_KERNEL break; case NORM_HEXA8: { - _nb_of_pts=8; _nb_of_sons=6; _dim=3; _quadratic_type=NORM_HEXA20; + _nb_of_pts=8; _nb_of_sons=6; _dim=3; _quadratic_type=NORM_HEXA20; _is_simplex=false; _sons_type[0]=NORM_QUAD4; _sons_type[1]=NORM_QUAD4; _sons_type[2]=NORM_QUAD4; _sons_type[3]=NORM_QUAD4; _sons_type[4]=NORM_QUAD4; _sons_type[5]=NORM_QUAD4; _sons_con[0][0]=0; _sons_con[0][1]=1; _sons_con[0][2]=2; _sons_con[0][3]=3; _nb_of_sons_con[0]=4; _sons_con[1][0]=4; _sons_con[1][1]=7; _sons_con[1][2]=6; _sons_con[1][3]=5; _nb_of_sons_con[1]=4; @@ -145,7 +145,7 @@ namespace INTERP_KERNEL break; case NORM_QUAD4: { - _nb_of_pts=4; _nb_of_sons=4; _dim=2; _quadratic_type=NORM_QUAD8; + _nb_of_pts=4; _nb_of_sons=4; _dim=2; _quadratic_type=NORM_QUAD8; _is_simplex=false; _sons_type[0]=NORM_SEG2; _sons_type[1]=NORM_SEG2; _sons_type[2]=NORM_SEG2; _sons_type[3]=NORM_SEG2; _sons_con[0][0]=0; _sons_con[0][1]=1; _nb_of_sons_con[0]=2; _sons_con[1][0]=1; _sons_con[1][1]=2; _nb_of_sons_con[1]=2; @@ -155,7 +155,7 @@ namespace INTERP_KERNEL break; case NORM_TRI3: { - _nb_of_pts=3; _nb_of_sons=3; _dim=2; _quadratic_type=NORM_TRI6; + _nb_of_pts=3; _nb_of_sons=3; _dim=2; _quadratic_type=NORM_TRI6; _is_simplex=true; _sons_type[0]=NORM_SEG2; _sons_type[1]=NORM_SEG2; _sons_type[2]=NORM_SEG2; _sons_con[0][0]=0; _sons_con[0][1]=1; _nb_of_sons_con[0]=2; _sons_con[1][0]=1; _sons_con[1][1]=2; _nb_of_sons_con[1]=2; @@ -164,7 +164,7 @@ namespace INTERP_KERNEL break; case NORM_TRI6: { - _nb_of_pts=6; _nb_of_sons=3; _dim=2; _linear_type=NORM_TRI3; + _nb_of_pts=6; _nb_of_sons=3; _dim=2; _linear_type=NORM_TRI3; _is_simplex=false; _sons_type[0]=NORM_SEG3; _sons_type[1]=NORM_SEG3; _sons_type[2]=NORM_SEG3; _sons_con[0][0]=0; _sons_con[0][1]=1; _sons_con[0][2]=3; _nb_of_sons_con[0]=3; _sons_con[1][0]=1; _sons_con[1][1]=2; _sons_con[1][2]=4; _nb_of_sons_con[1]=3; @@ -173,7 +173,7 @@ namespace INTERP_KERNEL break; case NORM_QUAD8: { - _nb_of_pts=8; _nb_of_sons=4; _dim=2; _linear_type=NORM_QUAD4; + _nb_of_pts=8; _nb_of_sons=4; _dim=2; _linear_type=NORM_QUAD4; _is_simplex=false; _sons_type[0]=NORM_SEG3; _sons_type[1]=NORM_SEG3; _sons_type[2]=NORM_SEG3; _sons_type[3]=NORM_SEG3; _sons_con[0][0]=0; _sons_con[0][1]=1; _sons_con[0][2]=4; _nb_of_sons_con[0]=3; _sons_con[1][0]=1; _sons_con[1][1]=2; _sons_con[1][2]=5; _nb_of_sons_con[1]=3; @@ -183,7 +183,7 @@ namespace INTERP_KERNEL break; case NORM_PYRA5: { - _nb_of_pts=5; _nb_of_sons=5; _dim=3; _quadratic_type=NORM_PYRA13; + _nb_of_pts=5; _nb_of_sons=5; _dim=3; _quadratic_type=NORM_PYRA13; _is_simplex=false; _sons_type[0]=NORM_QUAD4; _sons_type[1]=NORM_TRI3; _sons_type[2]=NORM_TRI3; _sons_type[3]=NORM_TRI3; _sons_type[4]=NORM_TRI3; _sons_con[0][0]=0; _sons_con[0][1]=1; _sons_con[0][2]=2; _sons_con[0][3]=3; _nb_of_sons_con[0]=4; _sons_con[1][0]=0; _sons_con[1][1]=4; _sons_con[1][2]=1; _nb_of_sons_con[1]=3; @@ -194,7 +194,7 @@ namespace INTERP_KERNEL break; case NORM_PENTA6: { - _nb_of_pts=6; _nb_of_sons=5; _dim=3; _quadratic_type=NORM_PENTA15; + _nb_of_pts=6; _nb_of_sons=5; _dim=3; _quadratic_type=NORM_PENTA15; _is_simplex=false; _sons_type[0]=NORM_TRI3; _sons_type[1]=NORM_TRI3; _sons_type[2]=NORM_QUAD4; _sons_type[3]=NORM_QUAD4; _sons_type[4]=NORM_QUAD4; _sons_con[0][0]=0; _sons_con[0][1]=1; _sons_con[0][2]=2; _nb_of_sons_con[0]=3; _sons_con[1][0]=3; _sons_con[1][1]=5; _sons_con[1][2]=4; _nb_of_sons_con[1]=3; @@ -205,7 +205,7 @@ namespace INTERP_KERNEL break; case NORM_TETRA10: { - _nb_of_pts=10; _nb_of_sons=4; _dim=3; _linear_type=NORM_TETRA4; + _nb_of_pts=10; _nb_of_sons=4; _dim=3; _linear_type=NORM_TETRA4; _is_simplex=false; _sons_type[0]=NORM_TRI6; _sons_type[1]=NORM_TRI6; _sons_type[2]=NORM_TRI6; _sons_type[3]=NORM_TRI6; _sons_con[0][0]=0; _sons_con[0][1]=1; _sons_con[0][2]=2; _sons_con[0][3]=4; _sons_con[0][4]=5; _sons_con[0][5]=6; _nb_of_sons_con[0]=6; _sons_con[1][0]=0; _sons_con[1][1]=3; _sons_con[1][2]=1; _sons_con[1][3]=7; _sons_con[1][4]=8; _sons_con[1][5]=4; _nb_of_sons_con[1]=6; @@ -215,7 +215,7 @@ namespace INTERP_KERNEL break; case NORM_PYRA13: { - _nb_of_pts=13; _nb_of_sons=5; _dim=3; _linear_type=NORM_PYRA5; + _nb_of_pts=13; _nb_of_sons=5; _dim=3; _linear_type=NORM_PYRA5; _is_simplex=false; _sons_type[0]=NORM_QUAD8; _sons_type[1]=NORM_TRI6; _sons_type[2]=NORM_TRI6; _sons_type[3]=NORM_TRI6; _sons_type[4]=NORM_TRI6; _sons_con[0][0]=0; _sons_con[0][1]=1; _sons_con[0][2]=2; _sons_con[0][3]=3; _sons_con[0][4]=5; _sons_con[0][5]=6; _sons_con[0][6]=7; _sons_con[0][7]=8; _nb_of_sons_con[0]=8; _sons_con[1][0]=0; _sons_con[1][1]=4; _sons_con[1][2]=1; _sons_con[1][3]=9; _sons_con[1][4]=10; _sons_con[1][5]=5; _nb_of_sons_con[1]=6; @@ -226,7 +226,7 @@ namespace INTERP_KERNEL break; case NORM_PENTA15: { - _nb_of_pts=15; _nb_of_sons=5; _dim=3; _linear_type=NORM_PENTA6; + _nb_of_pts=15; _nb_of_sons=5; _dim=3; _linear_type=NORM_PENTA6; _is_simplex=false; _sons_type[0]=NORM_TRI6; _sons_type[1]=NORM_TRI6; _sons_type[2]=NORM_QUAD8; _sons_type[3]=NORM_QUAD8; _sons_type[4]=NORM_QUAD8; _sons_con[0][0]=0; _sons_con[0][1]=1; _sons_con[0][2]=2; _sons_con[0][3]=6; _sons_con[0][4]=7; _sons_con[0][5]=8; _nb_of_sons_con[0]=6; _sons_con[1][0]=3; _sons_con[1][1]=5; _sons_con[1][2]=4; _sons_con[1][3]=11; _sons_con[1][4]=10; _sons_con[1][5]=9; _nb_of_sons_con[1]=6; @@ -237,7 +237,7 @@ namespace INTERP_KERNEL break; case NORM_HEXA20: { - _nb_of_pts=20; _nb_of_sons=6; _dim=3; _linear_type=NORM_HEXA8; + _nb_of_pts=20; _nb_of_sons=6; _dim=3; _linear_type=NORM_HEXA8; _is_simplex=false; _sons_type[0]=NORM_QUAD8; _sons_type[1]=NORM_QUAD8; _sons_type[2]=NORM_QUAD8; _sons_type[3]=NORM_QUAD8; _sons_type[4]=NORM_QUAD8; _sons_type[5]=NORM_QUAD8; _sons_con[0][0]=0; _sons_con[0][1]=1; _sons_con[0][2]=2; _sons_con[0][3]=3; _sons_con[0][4]=8; _sons_con[0][5]=9; _sons_con[0][6]=10; _sons_con[0][7]=11; _nb_of_sons_con[0]=8; _sons_con[1][0]=4; _sons_con[1][1]=7; _sons_con[1][2]=6; _sons_con[1][3]=5; _sons_con[1][4]=15; _sons_con[1][5]=14; _sons_con[1][6]=13; _sons_con[1][7]=12; _nb_of_sons_con[1]=8; @@ -249,12 +249,12 @@ namespace INTERP_KERNEL break; case NORM_POLYGON: { - _nb_of_pts=0; _nb_of_sons=0; _dim=2; _dyn=true; _extruded_type=NORM_POLYHED; + _nb_of_pts=0; _nb_of_sons=0; _dim=2; _dyn=true; _extruded_type=NORM_POLYHED; _is_simplex=false; } break; case NORM_POLYHED: { - _nb_of_pts=0; _nb_of_sons=0; _dim=3; _dyn=true; + _nb_of_pts=0; _nb_of_sons=0; _dim=3; _dyn=true; _is_simplex=false; } break; case NORM_ERROR: diff --git a/src/INTERP_KERNEL/CellModel.hxx b/src/INTERP_KERNEL/CellModel.hxx index 23b735512..a7a0828aa 100644 --- a/src/INTERP_KERNEL/CellModel.hxx +++ b/src/INTERP_KERNEL/CellModel.hxx @@ -46,6 +46,7 @@ namespace INTERP_KERNEL INTERPKERNEL_EXPORT bool isQuadratic() const { return _quadratic; } INTERPKERNEL_EXPORT unsigned getDimension() const { return _dim; } INTERPKERNEL_EXPORT bool isCompatibleWith(NormalizedCellType type) const; + INTERPKERNEL_EXPORT bool isSimplex() const { return _is_simplex; } //! sonId is in C format. INTERPKERNEL_EXPORT const unsigned *getNodesConstituentTheSon(unsigned sonId) const { return _sons_con[sonId]; } INTERPKERNEL_EXPORT unsigned getNumberOfNodes() const { return _nb_of_pts; } @@ -63,6 +64,7 @@ namespace INTERP_KERNEL private: bool _dyn; bool _quadratic; + bool _is_simplex; unsigned _dim; unsigned _nb_of_pts; unsigned _nb_of_sons; diff --git a/src/INTERP_KERNEL/InterpKernelCellSimplify.cxx b/src/INTERP_KERNEL/InterpKernelCellSimplify.cxx new file mode 100644 index 000000000..57b668c28 --- /dev/null +++ b/src/INTERP_KERNEL/InterpKernelCellSimplify.cxx @@ -0,0 +1,452 @@ +// 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 +// + +#include "InterpKernelCellSimplify.hxx" +#include "CellModel.hxx" + +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace INTERP_KERNEL; + +/*! + * This method takes as input a cell with type 'type' and whose connectivity is defined by (conn,lgth) + * It retrieves the same cell with a potentially different type (in return) whose connectivity is defined by (retConn,retLgth) + * \b WARNING for optimization reason the arrays 'retConn' and 'conn' can overlapped ! + */ +INTERP_KERNEL::NormalizedCellType CellSimplify::simplifyDegeneratedCell(INTERP_KERNEL::NormalizedCellType type, const int *conn, int lgth, + int *retConn, int& retLgth) throw(INTERP_KERNEL::Exception) +{ + const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::getCellModel(type); + std::set c(conn,conn+lgth); + c.erase(-1); + bool isObviousNonDegeneratedCell=((int)c.size()==lgth); + if(cm.isQuadratic() || isObviousNonDegeneratedCell) + {//quadratic do nothing for the moment. + retLgth=lgth; + int *tmp=new int[lgth];//no direct std::copy ! overlapping of conn and retConn ! + std::copy(conn,conn+lgth,tmp); + std::copy(tmp,tmp+lgth,retConn); + delete [] tmp; + return type; + } + if(cm.getDimension()==2) + { + int *tmp=new int[lgth]; + tmp[0]=conn[0]; + int newPos=1; + for(int i=1;i faces; + for(unsigned j=0;j nodes(conn,conn+lgth); + nodes.erase(-1); + int nbOfNodes=nodes.size(); + int magicNumber=100*nbOfNodes+nbOfFaces; + switch(magicNumber) + { + case 806: + return tryToUnPolyHex8(conn,nbOfFaces,lgth,retConn,retLgth); + case 605: + return tryToUnPolyPenta6(conn,nbOfFaces,lgth,retConn,retLgth); + case 505: + return tryToUnPolyPyra5(conn,nbOfFaces,lgth,retConn,retLgth); + case 404: + return tryToUnPolyTetra4(conn,nbOfFaces,lgth,retConn,retLgth); + default: + retLgth=lgth; + std::copy(conn,conn+lgth,retConn); + return INTERP_KERNEL::NORM_POLYHED; + } +} + +bool CellSimplify::orientOppositeFace(const int *baseFace, int *retConn, const int *sideFace, int lgthBaseFace) +{ + std::vector tmp2; + std::set bases(baseFace,baseFace+lgthBaseFace); + std::set sides(sideFace,sideFace+4); + std::set_intersection(bases.begin(),bases.end(),sides.begin(),sides.end(),std::back_insert_iterator< std::vector >(tmp2)); + if(tmp2.size()!=2) + return false; + std::vector< std::pair > baseEdges(lgthBaseFace); + std::vector< std::pair > oppEdges(lgthBaseFace); + std::vector< std::pair > sideEdges(4); + for(int i=0;i(baseFace[i],baseFace[(i+1)%lgthBaseFace]); + oppEdges[i]=std::pair(retConn[i],retConn[(i+1)%lgthBaseFace]); + } + for(int i=0;i<4;i++) + sideEdges[i]=std::pair(sideFace[i],sideFace[(i+1)%4]); + std::vector< std::pair > tmp; + std::set< std::pair > baseEdgesS(baseEdges.begin(),baseEdges.end()); + std::set< std::pair > sideEdgesS(sideEdges.begin(),sideEdges.end()); + std::set_intersection(baseEdgesS.begin(),baseEdgesS.end(),sideEdgesS.begin(),sideEdgesS.end(),std::back_insert_iterator< std::vector< std::pair > >(tmp)); + if(tmp.empty()) + { + //reverse sideFace + for(int i=0;i<4;i++) + { + std::pair p=sideEdges[i]; + std::pair r(p.second,p.first); + sideEdges[i]=r; + } + //end reverse sideFace + std::set< std::pair > baseEdgesS(baseEdges.begin(),baseEdges.end()); + std::set< std::pair > sideEdgesS(sideEdges.begin(),sideEdges.end()); + std::set_intersection(baseEdgesS.begin(),baseEdgesS.end(),sideEdgesS.begin(),sideEdgesS.end(),std::back_insert_iterator< std::vector< std::pair > >(tmp)); + if(tmp.empty()) + return false; + } + if(tmp.size()!=1) + return false; + bool found=false; + std::pair pInOpp; + for(int i=0;i<4 && !found;i++) + {//finding the pair(edge) in sideFace that do not include any node of tmp[0] edge + found=(tmp[0].first!=sideEdges[i].first && tmp[0].first!=sideEdges[i].second && + tmp[0].second!=sideEdges[i].first && tmp[0].second!=sideEdges[i].second); + if(found) + {//found ! reverse it + pInOpp.first=sideEdges[i].second; + pInOpp.second=sideEdges[i].first; + } + } + if(!found) + return false; + int pos=std::distance(baseEdges.begin(),std::find(baseEdges.begin(),baseEdges.end(),tmp[0])); + std::vector< std::pair >::iterator it=std::find(oppEdges.begin(),oppEdges.end(),pInOpp); + if(it==oppEdges.end())//the opposite edge of side face is not found opposite face ... maybe problem of orientation of polyhedron + return false; + int pos2=std::distance(oppEdges.begin(),it); + int offset=pos-pos2; + if(offset<0) + offset+=lgthBaseFace; + //this is the end copy the result + int *tmp3=new int[lgthBaseFace]; + for(int i=0;i(),(int)INTERP_KERNEL::NORM_QUAD4))==conn+lgth+nbOfFaces) + {//6 faces are QUAD4. + int oppositeFace=-1; + std::set conn1(conn,conn+4); + for(int i=1;i<6 && oppositeFace<0;i++) + { + std::vector tmp; + std::set conn2(conn+5*i,conn+5*i+4); + std::set_intersection(conn1.begin(),conn1.end(),conn2.begin(),conn2.end(),std::back_insert_iterator< std::vector >(tmp)); + if(tmp.empty()) + oppositeFace=i; + } + if(oppositeFace>=1) + {//oppositeFace of face#0 found. + int tmp2[4]; + if(tryToArrangeOppositeFace(conn,lgth,4,conn,conn+5*oppositeFace,6,tmp2)) + { + std::copy(conn,conn+4,retConn); + std::copy(tmp2,tmp2+4,retConn+4); + retLgth=8; + return INTERP_KERNEL::NORM_HEXA8; + } + } + } + retLgth=lgth; + std::copy(conn,conn+lgth,retConn); + return INTERP_KERNEL::NORM_POLYHED; +} + +/*! + * Cell with 'conn' connectivity has been detected as a good candidate. Full check of this. If yes NORM_PENTA6 is returned. + * If fails a POLYHED is returned. + */ +INTERP_KERNEL::NormalizedCellType CellSimplify::tryToUnPolyPenta6(const int *conn, int nbOfFaces, int lgth, int *retConn, int& retLgth) +{ + int nbOfTriFace=std::count(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_TRI3); + int nbOfQuadFace=std::count(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_QUAD4); + if(nbOfTriFace==2 && nbOfQuadFace==3) + { + int tri3_0=std::distance(conn+lgth,std::find(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_TRI3)); + int tri3_1=std::distance(conn+lgth,std::find(conn+lgth+tri3_0+1,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_TRI3)); + const int *tri_0=0,*tri_1=0; + const int *w=conn; + for(int i=0;i<5;i++) + { + if(i==tri3_0) + tri_0=w; + if(i==tri3_1) + tri_1=w; + w=std::find(w,conn+lgth,-1); + w++; + } + std::vector tmp; + std::set conn1(tri_0,tri_0+3); + std::set conn2(tri_1,tri_1+3); + std::set_intersection(conn1.begin(),conn1.end(),conn2.begin(),conn2.end(),std::back_insert_iterator< std::vector >(tmp)); + if(tmp.empty()) + { + int tmp2[3]; + if(tryToArrangeOppositeFace(conn,lgth,3,tri_0,tri_1,5,tmp2)) + { + std::copy(conn,conn+4,retConn); + std::copy(tmp2,tmp2+3,retConn+3); + retLgth=6; + return INTERP_KERNEL::NORM_PENTA6; + } + } + } + retLgth=lgth; + std::copy(conn,conn+lgth,retConn); + return INTERP_KERNEL::NORM_POLYHED; +} + +/*! + * Cell with 'conn' connectivity has been detected as a good candidate. Full check of this. If yes NORM_PYRA5 is returned. + * If fails a POLYHED is returned. + */ +INTERP_KERNEL::NormalizedCellType CellSimplify::tryToUnPolyPyra5(const int *conn, int nbOfFaces, int lgth, int *retConn, int& retLgth) +{ + int nbOfTriFace=std::count(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_TRI3); + int nbOfQuadFace=std::count(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_QUAD4); + if(nbOfTriFace==4 && nbOfQuadFace==1) + { + int quad4_pos=std::distance(conn+lgth,std::find(conn+lgth,conn+lgth+nbOfFaces,(int)INTERP_KERNEL::NORM_QUAD4)); + const int *quad4=0; + const int *w=conn; + for(int i=0;i<5 && quad4==0;i++) + { + if(i==quad4_pos) + quad4=w; + w=std::find(w,conn+lgth,-1); + w++; + } + std::set quad4S(quad4,quad4+4); + w=conn; + bool ok=true; + int point=-1; + for(int i=0;i<5 && ok;i++) + { + if(i!=quad4_pos) + { + std::vector tmp; + std::set conn2(w,w+3); + std::set_intersection(conn2.begin(),conn2.end(),quad4S.begin(),quad4S.end(),std::back_insert_iterator< std::vector >(tmp)); + ok=tmp.size()==2; + tmp.clear(); + std::set_difference(conn2.begin(),conn2.end(),quad4S.begin(),quad4S.end(),std::back_insert_iterator< std::vector >(tmp)); + ok=ok && tmp.size()==1; + if(ok) + { + if(point>=0) + ok=point==tmp[0]; + else + point=tmp[0]; + } + } + w=std::find(w,conn+lgth,-1); + w++; + } + if(ok && point>=0) + { + std::copy(quad4,quad4+4,retConn); + retConn[4]=point; + retLgth=5; + return INTERP_KERNEL::NORM_PYRA5; + } + } + retLgth=lgth; + std::copy(conn,conn+lgth,retConn); + return INTERP_KERNEL::NORM_POLYHED; +} + +/*! + * Cell with 'conn' connectivity has been detected as a good candidate. Full check of this. If yes NORM_TETRA4 is returned. + * If fails a POLYHED is returned. + */ +INTERP_KERNEL::NormalizedCellType CellSimplify::tryToUnPolyTetra4(const int *conn, int nbOfFaces, int lgth, int *retConn, int& retLgth) +{ + if(std::find_if(conn+lgth,conn+lgth+nbOfFaces,std::bind2nd(std::not_equal_to(),(int)INTERP_KERNEL::NORM_TRI3))==conn+lgth+nbOfFaces) + { + std::set tribase(conn,conn+3); + int point=-1; + bool ok=true; + for(int i=1;i<4 && ok;i++) + { + std::vector tmp; + std::set conn2(conn+i*4,conn+4*i+3); + std::set_intersection(conn2.begin(),conn2.end(),tribase.begin(),tribase.end(),std::back_insert_iterator< std::vector >(tmp)); + ok=tmp.size()==2; + tmp.clear(); + std::set_difference(conn2.begin(),conn2.end(),tribase.begin(),tribase.end(),std::back_insert_iterator< std::vector >(tmp)); + ok=ok && tmp.size()==1; + if(ok) + { + if(point>=0) + ok=point==tmp[0]; + else + point=tmp[0]; + } + } + if(ok && point>=0) + { + std::copy(conn,conn+3,retConn); + retConn[3]=point; + retLgth=4; + return INTERP_KERNEL::NORM_TETRA4; + } + } + retLgth=lgth; + std::copy(conn,conn+lgth,retConn); + return INTERP_KERNEL::NORM_POLYHED; +} diff --git a/src/INTERP_KERNEL/InterpKernelCellSimplify.hxx b/src/INTERP_KERNEL/InterpKernelCellSimplify.hxx new file mode 100644 index 000000000..3c2ec3925 --- /dev/null +++ b/src/INTERP_KERNEL/InterpKernelCellSimplify.hxx @@ -0,0 +1,47 @@ +// 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 __INTERPKERNELCELLSIMPLIFY_HXX__ +#define __INTERPKERNELCELLSIMPLIFY_HXX__ + +#include "NormalizedUnstructuredMesh.hxx" +#include "InterpKernelException.hxx" + +namespace INTERP_KERNEL +{ + class CellSimplify + { + public: + static INTERP_KERNEL::NormalizedCellType simplifyDegeneratedCell(INTERP_KERNEL::NormalizedCellType type, const int *conn, int lgth, + int *retConn, int& retLgth) throw(INTERP_KERNEL::Exception); + static int *getFullPolyh3DCell(INTERP_KERNEL::NormalizedCellType type, const int *conn, int lgth, + int& retNbOfFaces, int& retLgth); + static INTERP_KERNEL::NormalizedCellType tryToUnPoly2D(const int *conn, int lgth, int *retConn, int& retLgth); + static INTERP_KERNEL::NormalizedCellType tryToUnPoly3D(const int *conn, int nbOfFaces, int lgth, int *retConn, int& retLgth); + static INTERP_KERNEL::NormalizedCellType tryToUnPolyHex8(const int *conn, int nbOfFaces, int lgth, int *retConn, int& retLgth); + static INTERP_KERNEL::NormalizedCellType tryToUnPolyPenta6(const int *conn, int nbOfFaces, int lgth, int *retConn, int& retLgth); + static INTERP_KERNEL::NormalizedCellType tryToUnPolyPyra5(const int *conn, int nbOfFaces, int lgth, int *retConn, int& retLgth); + static INTERP_KERNEL::NormalizedCellType tryToUnPolyTetra4(const int *conn, int nbOfFaces, int lgth, int *retConn, int& retLgth); + static bool tryToArrangeOppositeFace(const int *conn, int lgth, int lgthBaseFace, const int *baseFace, const int *oppFaceId, int nbOfFaces, int *retConnOfOppFace); + static bool isWellOriented(const int *baseFace, int *retConn, const int *sideFace, int lgthBaseFace); + static bool orientOppositeFace(const int *baseFace, int *retConn, const int *sideFace, int lgthBaseFace); + }; +} + +#endif diff --git a/src/INTERP_KERNEL/Makefile.am b/src/INTERP_KERNEL/Makefile.am index e3dfc56cb..53a6650f0 100644 --- a/src/INTERP_KERNEL/Makefile.am +++ b/src/INTERP_KERNEL/Makefile.am @@ -44,6 +44,7 @@ PointLocator2DIntersector.txx \ INTERPKERNELDefines.hxx \ InterpKernelMatrix.hxx \ InterpKernelMeshQuality.hxx \ +InterpKernelCellSimplify.hxx \ Interpolation.hxx \ Interpolation.txx \ Interpolation2D.hxx \ @@ -189,7 +190,8 @@ dist_libinterpkernel_la_SOURCES = \ UnitTetraIntersectionBary.cxx \ InterpolationOptions.cxx \ DirectedBoundingBox.cxx \ - InterpKernelMeshQuality.cxx + InterpKernelMeshQuality.cxx \ + InterpKernelCellSimplify.cxx libinterpkernel_la_CPPFLAGS=-I$(srcdir)/Geometric2D -I$(srcdir)/Bases diff --git a/src/INTERP_KERNEL/SplitterTetra.txx b/src/INTERP_KERNEL/SplitterTetra.txx index d1e084876..fd0d06a82 100644 --- a/src/INTERP_KERNEL/SplitterTetra.txx +++ b/src/INTERP_KERNEL/SplitterTetra.txx @@ -242,7 +242,7 @@ namespace INTERP_KERNEL { // get sons connectivity NormalizedCellType faceType; - int *faceNodes, nbFaceNodes; + int *faceNodes, nbFaceNodes=-1; if ( cellModelCell.isDynamic() ) { faceNodes=new int[nbOfNodes4Type]; diff --git a/src/INTERP_KERNEL/TriangulationIntersector.txx b/src/INTERP_KERNEL/TriangulationIntersector.txx index de8a83954..0c17ed0c9 100644 --- a/src/INTERP_KERNEL/TriangulationIntersector.txx +++ b/src/INTERP_KERNEL/TriangulationIntersector.txx @@ -191,8 +191,7 @@ namespace INTERP_KERNEL //Compute the intersection area double inter_area[SPACEDIM], total_area = 0.; - double total_barycenter[SPACEDIM]; - total_barycenter[0]=total_barycenter[1] = 0.; + double total_barycenter[SPACEDIM]={0.,0.}; const ConnType nbNodesT=targetCell.size()/SPACEDIM; for(ConnType iT = 1; iTdecrRef(); } +void MEDCouplingFieldDiscretization::renumberEntitiesFromN2OArr(const int *new2OldPtr, int new2OldSz, DataArrayDouble *arr, const char *msg) +{ + int nbOfComp=arr->getNumberOfComponents(); + DataArrayDouble *arrCpy=arr->deepCopy(); + const double *ptSrc=arrCpy->getConstPointer(); + arr->reAlloc(new2OldSz); + double *ptToFill=arr->getPointer(); + for(int i=0;idecrRef(); +} + MEDCouplingFieldDiscretization::~MEDCouplingFieldDiscretization() { } @@ -395,6 +410,11 @@ void MEDCouplingFieldDiscretizationP0::renumberValuesOnCells(const MEDCouplingMe renumberEntitiesFromO2NArr(old2New,arr,"Cell"); } +void MEDCouplingFieldDiscretizationP0::renumberValuesOnCellsR(const MEDCouplingMesh *mesh, const int *new2old, int newSz, DataArrayDouble *arr) const +{ + renumberEntitiesFromN2OArr(new2old,newSz,arr,"Cell"); +} + /*! * This method returns a submesh of 'mesh' instance constituting cell ids contained in array defined as an interval [start;end). * @param di is an array returned that specifies entity ids (here cells ids) in mesh 'mesh' of entity in returned submesh. @@ -531,6 +551,13 @@ void MEDCouplingFieldDiscretizationP1::renumberValuesOnCells(const MEDCouplingMe { } +/*! + * Nothing to do it's not a bug. + */ +void MEDCouplingFieldDiscretizationP1::renumberValuesOnCellsR(const MEDCouplingMesh *mesh, const int *new2old, int newSz, DataArrayDouble *arr) const +{ +} + /*! * This method returns a submesh of 'mesh' instance constituting cell ids contained in array defined as an interval [start;end). * @param di is an array returned that specifies entity ids (here nodes ids) in mesh 'mesh' of entity in returned submesh. @@ -881,6 +908,11 @@ void MEDCouplingFieldDiscretizationGauss::renumberValuesOnCells(const MEDCouplin throw INTERP_KERNEL::Exception("Not implemented yet !"); } +void MEDCouplingFieldDiscretizationGauss::renumberValuesOnCellsR(const MEDCouplingMesh *mesh, const int *new2old, int newSz, DataArrayDouble *arr) const +{ + throw INTERP_KERNEL::Exception("Number of cells has changed and becomes higher with some cells that have been split ! Unable to conserve the Gauss field !"); +} + void MEDCouplingFieldDiscretizationGauss::setGaussLocalizationOnType(const MEDCouplingMesh *m, INTERP_KERNEL::NormalizedCellType type, const std::vector& refCoo, const std::vector& gsCoo, const std::vector& wg) throw(INTERP_KERNEL::Exception) { @@ -1180,6 +1212,11 @@ void MEDCouplingFieldDiscretizationGaussNE::renumberValuesOnCells(const MEDCoupl throw INTERP_KERNEL::Exception("Not implemented yet !"); } +void MEDCouplingFieldDiscretizationGaussNE::renumberValuesOnCellsR(const MEDCouplingMesh *mesh, const int *new2old, int newSz, DataArrayDouble *arr) 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 316fa753d..ad629d3ed 100644 --- a/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx +++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx @@ -68,6 +68,7 @@ namespace ParaMEDMEM virtual MEDCouplingMesh *buildSubMeshData(const MEDCouplingMesh *mesh, const int *start, const int *end, DataArrayInt *&di) const = 0; virtual void renumberValuesOnNodes(const int *old2New, DataArrayDouble *arr) const = 0; virtual void renumberValuesOnCells(const MEDCouplingMesh *mesh, const int *old2New, DataArrayDouble *arr) const = 0; + virtual void renumberValuesOnCellsR(const MEDCouplingMesh *mesh, const int *new2old, int newSz, DataArrayDouble *arr) const = 0; virtual void getSerializationIntArray(DataArrayInt *& arr) const; virtual void getTinySerializationIntInformation(std::vector& tinyInfo) const; virtual void getTinySerializationDbleInformation(std::vector& tinyInfo) const; @@ -88,6 +89,7 @@ namespace ParaMEDMEM protected: MEDCouplingFieldDiscretization(); static void renumberEntitiesFromO2NArr(const int *old2NewPtr, DataArrayDouble *arr, const char *msg); + static void renumberEntitiesFromN2OArr(const int *new2OldPtr, int new2OldSz, DataArrayDouble *arr, const char *msg); protected: double _precision; static const double DFLT_PRECISION; @@ -113,6 +115,7 @@ namespace ParaMEDMEM void getValueOnPos(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, int i, int j, int k, double *res) const; void renumberValuesOnNodes(const int *old2New, DataArrayDouble *arr) const; void renumberValuesOnCells(const MEDCouplingMesh *mesh, const int *old2New, DataArrayDouble *arr) const; + void renumberValuesOnCellsR(const MEDCouplingMesh *mesh, const int *new2old, int newSz, DataArrayDouble *arr) const; MEDCouplingMesh *buildSubMeshData(const MEDCouplingMesh *mesh, const int *start, const int *end, DataArrayInt *&di) const; public: static const char REPR[]; @@ -140,6 +143,7 @@ namespace ParaMEDMEM MEDCouplingMesh *buildSubMeshData(const MEDCouplingMesh *mesh, const int *start, const int *end, DataArrayInt *&di) const; void renumberValuesOnNodes(const int *old2New, DataArrayDouble *arr) const; void renumberValuesOnCells(const MEDCouplingMesh *mesh, const int *old2New, DataArrayDouble *arr) const; + void renumberValuesOnCellsR(const MEDCouplingMesh *mesh, const int *new2old, int newSz, DataArrayDouble *arr) const; public: static const char REPR[]; static const TypeOfField TYPE; @@ -195,6 +199,7 @@ namespace ParaMEDMEM MEDCouplingMesh *buildSubMeshData(const MEDCouplingMesh *mesh, const int *start, const int *end, DataArrayInt *&di) const; void renumberValuesOnNodes(const int *old2New, DataArrayDouble *arr) const; void renumberValuesOnCells(const MEDCouplingMesh *mesh, const int *old2New, DataArrayDouble *arr) const; + void renumberValuesOnCellsR(const MEDCouplingMesh *mesh, const int *new2old, int newSz, 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, @@ -244,6 +249,7 @@ namespace ParaMEDMEM MEDCouplingMesh *buildSubMeshData(const MEDCouplingMesh *mesh, const int *start, const int *end, DataArrayInt *&di) const; void renumberValuesOnNodes(const int *old2New, DataArrayDouble *arr) const; void renumberValuesOnCells(const MEDCouplingMesh *mesh, const int *old2New, DataArrayDouble *arr) const; + void renumberValuesOnCellsR(const MEDCouplingMesh *mesh, const int *new2old, int newSz, DataArrayDouble *arr) const; protected: MEDCouplingFieldDiscretizationGaussNE(const MEDCouplingFieldDiscretizationGaussNE& other); public: diff --git a/src/MEDCoupling/MEDCouplingFieldDouble.cxx b/src/MEDCoupling/MEDCouplingFieldDouble.cxx index 2c62b3f44..60b5caab8 100644 --- a/src/MEDCoupling/MEDCouplingFieldDouble.cxx +++ b/src/MEDCoupling/MEDCouplingFieldDouble.cxx @@ -75,6 +75,17 @@ void MEDCouplingFieldDouble::copyTinyStringsFrom(const MEDCouplingFieldDouble *o } } +/*! + * Copy only times, order, iteration from other. The underlying mesh is not impacted by this method. + */ +void MEDCouplingFieldDouble::copyTinyAttrFrom(const MEDCouplingFieldDouble *other) throw(INTERP_KERNEL::Exception) +{ + if(other) + { + _time_discr->copyTinyAttrFrom(*other->_time_discr); + } +} + std::string MEDCouplingFieldDouble::simpleRepr() const { std::ostringstream ret; @@ -1072,6 +1083,27 @@ bool MEDCouplingFieldDouble::zipConnectivity(int compType) throw(INTERP_KERNEL:: return false; } +/*! + * This method applyies ParaMEDMEM::MEDCouplingUMesh::simplexize on 'this->_mesh'. + * The semantic of 'policy' is given in ParaMEDMEM::MEDCouplingUMesh::simplexize method. + */ +bool MEDCouplingFieldDouble::simplexize(int policy) throw(INTERP_KERNEL::Exception) +{ + int oldNbOfCells=_mesh->getNumberOfCells(); + MEDCouplingAutoRefCountObjectPtr meshC2(_mesh->deepCpy()); + MEDCouplingAutoRefCountObjectPtr arr=meshC2->simplexize(policy); + int newNbOfCells=meshC2->getNumberOfCells(); + if(oldNbOfCells==newNbOfCells) + return false; + std::vector arrays; + _time_discr->getArrays(arrays); + for(std::vector::const_iterator iter=arrays.begin();iter!=arrays.end();iter++) + if(*iter) + _type->renumberValuesOnCellsR(_mesh,arr->getConstPointer(),arr->getNbOfElems(),*iter); + setMesh(meshC2); + return true; +} + MEDCouplingFieldDouble *MEDCouplingFieldDouble::doublyContractedProduct() const throw(INTERP_KERNEL::Exception) { MEDCouplingTimeDiscretization *td=_time_discr->doublyContractedProduct(); diff --git a/src/MEDCoupling/MEDCouplingFieldDouble.hxx b/src/MEDCoupling/MEDCouplingFieldDouble.hxx index 48c0bace0..19f804c31 100644 --- a/src/MEDCoupling/MEDCouplingFieldDouble.hxx +++ b/src/MEDCoupling/MEDCouplingFieldDouble.hxx @@ -33,6 +33,7 @@ namespace ParaMEDMEM public: static MEDCouplingFieldDouble *New(TypeOfField type, TypeOfTimeDiscretization td=NO_TIME); void copyTinyStringsFrom(const MEDCouplingFieldDouble *other) throw(INTERP_KERNEL::Exception); + void copyTinyAttrFrom(const MEDCouplingFieldDouble *other) throw(INTERP_KERNEL::Exception); std::string simpleRepr() const; std::string advancedRepr() const; bool isEqual(const MEDCouplingField *other, double meshPrec, double valsPrec) const; @@ -54,6 +55,8 @@ namespace ParaMEDMEM void checkCoherency() const throw(INTERP_KERNEL::Exception); NatureOfField getNature() const { return _nature; } void setNature(NatureOfField nat) throw(INTERP_KERNEL::Exception); + void setTimeTolerance(double val) { _time_discr->setTimeTolerance(val); } + double getTimeTolerance() const { return _time_discr->getTimeTolerance(); } void setTime(double val, int iteration, int order) { _time_discr->setTime(val,iteration,order); } void setStartTime(double val, int iteration, int order) { _time_discr->setStartTime(val,iteration,order); } void setEndTime(double val, int iteration, int order) { _time_discr->setEndTime(val,iteration,order); } @@ -111,6 +114,7 @@ namespace ParaMEDMEM bool mergeNodes(double eps) throw(INTERP_KERNEL::Exception); bool zipCoords() throw(INTERP_KERNEL::Exception); bool zipConnectivity(int compType) throw(INTERP_KERNEL::Exception); + bool simplexize(int policy) throw(INTERP_KERNEL::Exception); MEDCouplingFieldDouble *doublyContractedProduct() const throw(INTERP_KERNEL::Exception); MEDCouplingFieldDouble *determinant() const throw(INTERP_KERNEL::Exception); MEDCouplingFieldDouble *eigenValues() const throw(INTERP_KERNEL::Exception); diff --git a/src/MEDCoupling/MEDCouplingMemArray.cxx b/src/MEDCoupling/MEDCouplingMemArray.cxx index 5365cb97c..1b4a20b49 100644 --- a/src/MEDCoupling/MEDCouplingMemArray.cxx +++ b/src/MEDCoupling/MEDCouplingMemArray.cxx @@ -86,7 +86,7 @@ bool DataArray::areInfoEquals(const DataArray& other) const void DataArray::reprWithoutNameStream(std::ostream& stream) const { - stream << "Nb of components : "<< getNumberOfComponents() << "\n"; + stream << "Number of components : "<< getNumberOfComponents() << "\n"; stream << "Info of these components : "; for(std::vector::const_iterator iter=_info_on_compo.begin();iter!=_info_on_compo.end();iter++) stream << "\"" << *iter << "\" "; @@ -120,6 +120,17 @@ DataArrayDouble *DataArrayDouble::New() return new DataArrayDouble; } +bool DataArrayDouble::isAllocated() const +{ + return getConstPointer()!=0; +} + +void DataArrayDouble::checkAllocated() const throw(INTERP_KERNEL::Exception) +{ + if(!isAllocated()) + throw INTERP_KERNEL::Exception("DataArrayDouble::checkAllocated : Array is defined but not allocated ! Call alloc or setValues method first !"); +} + DataArrayDouble *DataArrayDouble::deepCopy() const { return new DataArrayDouble(*this); @@ -144,18 +155,47 @@ void DataArrayDouble::alloc(int nbOfTuple, int nbOfCompo) declareAsNew(); } -void DataArrayDouble::fillWithZero() +void DataArrayDouble::fillWithZero() throw(INTERP_KERNEL::Exception) +{ + fillWithValue(0.); +} + +void DataArrayDouble::fillWithValue(double val) throw(INTERP_KERNEL::Exception) { - _mem.fillWithValue(0.); + if(!getPointer()) + throw INTERP_KERNEL::Exception("DataArrayDouble::fillWithValue : allocate first !"); + _mem.fillWithValue(val); declareAsNew(); } -void DataArrayDouble::fillWithValue(double val) +void DataArrayDouble::iota(double init) throw(INTERP_KERNEL::Exception) { - _mem.fillWithValue(val); + double *ptr=getPointer(); + if(!ptr) + throw INTERP_KERNEL::Exception("DataArrayDouble::iota : allocate first !"); + if(getNumberOfComponents()!=1) + throw INTERP_KERNEL::Exception("DataArrayDouble::iota : works only for arrays with only one component!"); + int ntuples=getNumberOfTuples(); + for(int i=0;ivmax) + return false; + return true; +} + std::string DataArrayDouble::repr() const { std::ostringstream ret; @@ -208,8 +248,9 @@ bool DataArrayDouble::isEqualWithoutConsideringStr(const DataArrayDouble& other, return _mem.isEqual(other._mem,prec); } -void DataArrayDouble::reAlloc(int nbOfTuples) +void DataArrayDouble::reAlloc(int nbOfTuples) throw(INTERP_KERNEL::Exception) { + checkAllocated(); _mem.reAlloc(_info_on_compo.size()*nbOfTuples); _nb_of_tuples=nbOfTuples; declareAsNew(); @@ -227,6 +268,26 @@ DataArrayInt *DataArrayDouble::convertToIntArr() const return ret; } +DataArrayDouble *DataArrayDouble::fromNoInterlace() const throw(INTERP_KERNEL::Exception) +{ + if(_mem.isNull()) + throw INTERP_KERNEL::Exception("DataArrayDouble::fromNoInterlace : Not defined array !"); + double *tab=_mem.fromNoInterlace(getNumberOfComponents()); + DataArrayDouble *ret=DataArrayDouble::New(); + ret->useArray(tab,true,CPP_DEALLOC,getNumberOfTuples(),getNumberOfComponents()); + return ret; +} + +DataArrayDouble *DataArrayDouble::toNoInterlace() const throw(INTERP_KERNEL::Exception) +{ + if(_mem.isNull()) + throw INTERP_KERNEL::Exception("DataArrayDouble::fromNoInterlace : Not defined array !"); + double *tab=_mem.toNoInterlace(getNumberOfComponents()); + DataArrayDouble *ret=DataArrayDouble::New(); + ret->useArray(tab,true,CPP_DEALLOC,getNumberOfTuples(),getNumberOfComponents()); + return ret; +} + /*! * This method does \b not change the number of tuples after this call. * Only a permutation is done. If a permutation reduction is needed substr, or selectByTupleId should be used. @@ -377,6 +438,7 @@ DataArrayDouble *DataArrayDouble::substr(int tupleIdBg, int tupleIdEnd) const th */ DataArrayDouble *DataArrayDouble::changeNbOfComponents(int newNbOfComp, double dftValue) const throw(INTERP_KERNEL::Exception) { + checkAllocated(); DataArrayDouble *ret=DataArrayDouble::New(); ret->alloc(getNumberOfTuples(),newNbOfComp); const double *oldc=getConstPointer(); @@ -401,6 +463,7 @@ DataArrayDouble *DataArrayDouble::changeNbOfComponents(int newNbOfComp, double d DataArrayDouble *DataArrayDouble::keepSelectedComponents(const std::vector& compoIds) const throw(INTERP_KERNEL::Exception) { + checkAllocated(); MEDCouplingAutoRefCountObjectPtr ret(DataArrayDouble::New()); int newNbOfCompo=compoIds.size(); int oldNbOfCompo=getNumberOfComponents(); @@ -514,8 +577,9 @@ double DataArrayDouble::getAverageValue() const throw(INTERP_KERNEL::Exception) return ret/nbOfTuples; } -void DataArrayDouble::accumulate(double *res) const +void DataArrayDouble::accumulate(double *res) const throw(INTERP_KERNEL::Exception) { + checkAllocated(); const double *ptr=getConstPointer(); int nbTuple=getNumberOfTuples(); int nbComps=getNumberOfComponents(); @@ -524,8 +588,9 @@ void DataArrayDouble::accumulate(double *res) const std::transform(ptr+i*nbComps,ptr+(i+1)*nbComps,res,res,std::plus()); } -double DataArrayDouble::accumulate(int compId) const +double DataArrayDouble::accumulate(int compId) const throw(INTERP_KERNEL::Exception) { + checkAllocated(); const double *ptr=getConstPointer(); int nbTuple=getNumberOfTuples(); int nbComps=getNumberOfComponents(); @@ -537,6 +602,63 @@ double DataArrayDouble::accumulate(int compId) const return ret; } +DataArrayDouble *DataArrayDouble::fromPolarToCart() const throw(INTERP_KERNEL::Exception) +{ + int nbOfComp=getNumberOfComponents(); + if(nbOfComp!=2) + throw INTERP_KERNEL::Exception("DataArrayDouble::fromPolarToCart : must be an array with exactly 2 components !"); + int nbOfTuple=getNumberOfTuples(); + DataArrayDouble *ret=DataArrayDouble::New(); + ret->alloc(nbOfTuple,2); + double *w=ret->getPointer(); + const double *wIn=getConstPointer(); + for(int i=0;ialloc(getNumberOfTuples(),3); + double *w=ret->getPointer(); + const double *wIn=getConstPointer(); + for(int i=0;isetInfoOnComponent(2,getInfoOnComponent(2).c_str()); + return ret; +} + +DataArrayDouble *DataArrayDouble::fromSpherToCart() const throw(INTERP_KERNEL::Exception) +{ + int nbOfComp=getNumberOfComponents(); + if(nbOfComp!=3) + throw INTERP_KERNEL::Exception("DataArrayDouble::fromSpherToCart : must be an array with exactly 3 components !"); + int nbOfTuple=getNumberOfTuples(); + DataArrayDouble *ret=DataArrayDouble::New(); + ret->alloc(getNumberOfTuples(),3); + double *w=ret->getPointer(); + const double *wIn=getConstPointer(); + for(int i=0;ialloc(nbOfTuple,1); @@ -708,6 +831,7 @@ DataArrayDouble *DataArrayDouble::deviator() const throw(INTERP_KERNEL::Exceptio DataArrayDouble *DataArrayDouble::magnitude() const throw(INTERP_KERNEL::Exception) { + checkAllocated(); int nbOfComp=getNumberOfComponents(); DataArrayDouble *ret=DataArrayDouble::New(); int nbOfTuple=getNumberOfTuples(); @@ -726,6 +850,7 @@ DataArrayDouble *DataArrayDouble::magnitude() const throw(INTERP_KERNEL::Excepti DataArrayDouble *DataArrayDouble::maxPerTuple() const throw(INTERP_KERNEL::Exception) { + checkAllocated(); int nbOfComp=getNumberOfComponents(); DataArrayDouble *ret=DataArrayDouble::New(); int nbOfTuple=getNumberOfTuples(); @@ -739,6 +864,7 @@ DataArrayDouble *DataArrayDouble::maxPerTuple() const throw(INTERP_KERNEL::Excep void DataArrayDouble::sortPerTuple(bool asc) throw(INTERP_KERNEL::Exception) { + checkAllocated(); double *pt=getPointer(); int nbOfTuple=getNumberOfTuples(); int nbOfComp=getNumberOfComponents(); @@ -751,8 +877,9 @@ void DataArrayDouble::sortPerTuple(bool asc) throw(INTERP_KERNEL::Exception) declareAsNew(); } -void DataArrayDouble::applyLin(double a, double b, int compoId) +void DataArrayDouble::applyLin(double a, double b, int compoId) throw(INTERP_KERNEL::Exception) { + checkAllocated(); double *ptr=getPointer()+compoId; int nbOfComp=getNumberOfComponents(); int nbOfTuple=getNumberOfTuples(); @@ -763,6 +890,7 @@ void DataArrayDouble::applyLin(double a, double b, int compoId) DataArrayDouble *DataArrayDouble::applyFunc(int nbOfComp, FunctionToEvaluate func) const throw(INTERP_KERNEL::Exception) { + checkAllocated(); DataArrayDouble *newArr=DataArrayDouble::New(); int nbOfTuples=getNumberOfTuples(); int oldNbOfComp=getNumberOfComponents(); @@ -785,6 +913,7 @@ DataArrayDouble *DataArrayDouble::applyFunc(int nbOfComp, FunctionToEvaluate fun DataArrayDouble *DataArrayDouble::applyFunc(int nbOfComp, const char *func) const throw(INTERP_KERNEL::Exception) { + checkAllocated(); INTERP_KERNEL::ExprParser expr(func); expr.parse(); std::set vars; @@ -825,6 +954,7 @@ DataArrayDouble *DataArrayDouble::applyFunc(int nbOfComp, const char *func) cons DataArrayDouble *DataArrayDouble::applyFunc(const char *func) const throw(INTERP_KERNEL::Exception) { + checkAllocated(); INTERP_KERNEL::ExprParser expr(func); expr.parse(); expr.prepareExprEvaluationVec(); @@ -853,8 +983,9 @@ DataArrayDouble *DataArrayDouble::applyFunc(const char *func) const throw(INTERP return newArr; } -void DataArrayDouble::applyFuncFast32(const char *func) +void DataArrayDouble::applyFuncFast32(const char *func) throw(INTERP_KERNEL::Exception) { + checkAllocated(); INTERP_KERNEL::ExprParser expr(func); expr.parse(); char *funcStr=expr.compileX86(); @@ -869,8 +1000,9 @@ void DataArrayDouble::applyFuncFast32(const char *func) declareAsNew(); } -void DataArrayDouble::applyFuncFast64(const char *func) +void DataArrayDouble::applyFuncFast64(const char *func) throw(INTERP_KERNEL::Exception) { + checkAllocated(); INTERP_KERNEL::ExprParser expr(func); expr.parse(); char *funcStr=expr.compileX86_64(); @@ -888,7 +1020,7 @@ void DataArrayDouble::applyFuncFast64(const char *func) DataArrayInt *DataArrayDouble::getIdsInRange(double vmin, double vmax) const throw(INTERP_KERNEL::Exception) { if(getNumberOfComponents()!=1) - throw INTERP_KERNEL::Exception("DataArrayDouble::getIdsInRange : the default array must have only one component !"); + throw INTERP_KERNEL::Exception("DataArrayDouble::getIdsInRange : this must have exactly one component !"); const double *cptr=getConstPointer(); std::vector res; int nbOfTuples=getNumberOfTuples(); @@ -918,6 +1050,8 @@ DataArrayDouble *DataArrayDouble::aggregate(const DataArrayDouble *a1, const Dat DataArrayDouble *DataArrayDouble::dot(const DataArrayDouble *a1, const DataArrayDouble *a2) throw(INTERP_KERNEL::Exception) { + a1->checkAllocated(); + a2->checkAllocated(); int nbOfComp=a1->getNumberOfComponents(); if(nbOfComp!=a2->getNumberOfComponents()) throw INTERP_KERNEL::Exception("Nb of components mismatch for array dot !"); @@ -1166,6 +1300,17 @@ DataArrayInt *DataArrayInt::New() return new DataArrayInt; } +bool DataArrayInt::isAllocated() const +{ + return getConstPointer()!=0; +} + +void DataArrayInt::checkAllocated() const throw(INTERP_KERNEL::Exception) +{ + if(!isAllocated()) + throw INTERP_KERNEL::Exception("DataArrayInt::checkAllocated : Array is defined but not allocated ! Call alloc or setValues method first !"); +} + DataArrayInt *DataArrayInt::deepCopy() const { return new DataArrayInt(*this); @@ -1202,6 +1347,19 @@ void DataArrayInt::fillWithValue(int val) declareAsNew(); } +void DataArrayInt::iota(int init) throw(INTERP_KERNEL::Exception) +{ + int *ptr=getPointer(); + if(!ptr) + throw INTERP_KERNEL::Exception("DataArrayInt::iota : allocate first !"); + if(getNumberOfComponents()!=1) + throw INTERP_KERNEL::Exception("DataArrayInt::iota : works only for arrays with only one component!"); + int ntuples=getNumberOfTuples(); + for(int i=0;iuseArray(tab,true,CPP_DEALLOC,getNumberOfTuples(),getNumberOfComponents()); + return ret; +} + +DataArrayInt *DataArrayInt::toNoInterlace() const throw(INTERP_KERNEL::Exception) +{ + if(_mem.isNull()) + throw INTERP_KERNEL::Exception("DataArrayInt::toNoInterlace : Not defined array !"); + int *tab=_mem.toNoInterlace(getNumberOfComponents()); + DataArrayInt *ret=DataArrayInt::New(); + ret->useArray(tab,true,CPP_DEALLOC,getNumberOfTuples(),getNumberOfComponents()); + return ret; +} + void DataArrayInt::renumberInPlace(const int *old2New) { int nbTuples=getNumberOfTuples(); @@ -1418,6 +1596,19 @@ bool DataArrayInt::isIdentity() const return true; } +bool DataArrayInt::isUniform(int val) const +{ + if(getNumberOfComponents()!=1) + throw INTERP_KERNEL::Exception("DataArrayInt::isUniform : must be applied on DataArrayInt with only one component !"); + int nbOfTuples=getNumberOfTuples(); + const int *w=getConstPointer(); + const int *end=w+nbOfTuples; + for(;w!=end;w++) + if(*w!=val) + return false; + return true; +} + DataArrayDouble *DataArrayInt::convertToDblArr() const { DataArrayDouble *ret=DataArrayDouble::New(); @@ -1465,6 +1656,7 @@ DataArrayInt *DataArrayInt::substr(int tupleIdBg, int tupleIdEnd) const throw(IN */ DataArrayInt *DataArrayInt::changeNbOfComponents(int newNbOfComp, int dftValue) const throw(INTERP_KERNEL::Exception) { + checkAllocated(); DataArrayInt *ret=DataArrayInt::New(); ret->alloc(getNumberOfTuples(),newNbOfComp); const int *oldc=getConstPointer(); @@ -1487,8 +1679,9 @@ DataArrayInt *DataArrayInt::changeNbOfComponents(int newNbOfComp, int dftValue) return ret; } -void DataArrayInt::reAlloc(int nbOfTuples) +void DataArrayInt::reAlloc(int nbOfTuples) throw(INTERP_KERNEL::Exception) { + checkAllocated(); _mem.reAlloc(_info_on_compo.size()*nbOfTuples); _nb_of_tuples=nbOfTuples; declareAsNew(); @@ -1496,6 +1689,7 @@ void DataArrayInt::reAlloc(int nbOfTuples) DataArrayInt *DataArrayInt::keepSelectedComponents(const std::vector& compoIds) const throw(INTERP_KERNEL::Exception) { + checkAllocated(); MEDCouplingAutoRefCountObjectPtr ret(DataArrayInt::New()); int newNbOfCompo=compoIds.size(); int oldNbOfCompo=getNumberOfComponents(); @@ -1536,6 +1730,39 @@ void DataArrayInt::setArrayIn(DataArrayInt *newArray, DataArrayInt* &arrayToSet) } } +DataArrayInt *DataArrayInt::getIdsEqual(int val) const throw(INTERP_KERNEL::Exception) +{ + if(getNumberOfComponents()!=1) + throw INTERP_KERNEL::Exception("DataArrayInt::getIdsEqual : the array must have only one component !"); + const int *cptr=getConstPointer(); + std::vector res; + int nbOfTuples=getNumberOfTuples(); + for(int i=0;ialloc(res.size(),1); + std::copy(res.begin(),res.end(),ret->getPointer()); + return ret; +} + +DataArrayInt *DataArrayInt::getIdsEqualList(const std::vector& vals) const throw(INTERP_KERNEL::Exception) +{ + if(getNumberOfComponents()!=1) + throw INTERP_KERNEL::Exception("DataArrayInt::getIdsEqualList : the array must have only one component !"); + std::set vals2(vals.begin(),vals.end()); + const int *cptr=getConstPointer(); + std::vector res; + int nbOfTuples=getNumberOfTuples(); + for(int i=0;ialloc(res.size(),1); + std::copy(res.begin(),res.end(),ret->getPointer()); + return ret; +} + DataArrayInt *DataArrayInt::aggregate(const DataArrayInt *a1, const DataArrayInt *a2, int offsetA2) { int nbOfComp=a1->getNumberOfComponents(); diff --git a/src/MEDCoupling/MEDCouplingMemArray.hxx b/src/MEDCoupling/MEDCouplingMemArray.hxx index 0f5d600d9..57d7d8a0a 100644 --- a/src/MEDCoupling/MEDCouplingMemArray.hxx +++ b/src/MEDCoupling/MEDCouplingMemArray.hxx @@ -54,6 +54,7 @@ namespace ParaMEDMEM public: MemArray():_nb_of_elem(-1),_ownership(false),_dealloc(CPP_DEALLOC) { } MemArray(const MemArray& other); + bool isNull() const { return _pointer.isNull(); } const T *getConstPointerLoc(int offset) const { return _pointer.getConstPointerLoc(offset); } const T *getConstPointer() const { return _pointer.getConstPointer(); } T *getPointer() const { return _pointer.getPointer(); } @@ -64,6 +65,8 @@ namespace ParaMEDMEM void repr(int sl, std::ostream& stream) const; void reprZip(int sl, std::ostream& stream) const; void fillWithValue(const T& val); + T *fromNoInterlace(int nbOfComp) const; + T *toNoInterlace(int nbOfComp) const; void alloc(int nbOfElements); void reAlloc(int newNbOfElements); void useArray(const T *array, bool ownership, DeallocType type, int nbOfElem); @@ -113,11 +116,15 @@ namespace ParaMEDMEM { public: MEDCOUPLING_EXPORT static DataArrayDouble *New(); + MEDCOUPLING_EXPORT bool isAllocated() const; + MEDCOUPLING_EXPORT void checkAllocated() const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT DataArrayDouble *deepCopy() const; MEDCOUPLING_EXPORT DataArrayDouble *performCpy(bool deepCpy) const; MEDCOUPLING_EXPORT void alloc(int nbOfTuple, int nbOfCompo); - MEDCOUPLING_EXPORT void fillWithZero(); - MEDCOUPLING_EXPORT void fillWithValue(double val); + MEDCOUPLING_EXPORT void fillWithZero() throw(INTERP_KERNEL::Exception); + MEDCOUPLING_EXPORT void fillWithValue(double val) throw(INTERP_KERNEL::Exception); + MEDCOUPLING_EXPORT void iota(double init=0.) throw(INTERP_KERNEL::Exception); + MEDCOUPLING_EXPORT bool isUniform(double val, double eps) const; MEDCOUPLING_EXPORT std::string repr() const; MEDCOUPLING_EXPORT std::string reprZip() const; MEDCOUPLING_EXPORT void reprStream(std::ostream& stream) const; @@ -127,8 +134,10 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT bool isEqual(const DataArrayDouble& other, double prec) const; MEDCOUPLING_EXPORT bool isEqualWithoutConsideringStr(const DataArrayDouble& other, double prec) const; //!alloc or useArray should have been called before. - MEDCOUPLING_EXPORT void reAlloc(int nbOfTuples); + MEDCOUPLING_EXPORT void reAlloc(int nbOfTuples) throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT DataArrayInt *convertToIntArr() const; + MEDCOUPLING_EXPORT DataArrayDouble *fromNoInterlace() const throw(INTERP_KERNEL::Exception); + MEDCOUPLING_EXPORT DataArrayDouble *toNoInterlace() const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT void renumberInPlace(const int *old2New); MEDCOUPLING_EXPORT void renumberInPlaceR(const int *new2Old); MEDCOUPLING_EXPORT DataArrayDouble *renumber(const int *old2New) const; @@ -153,8 +162,11 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT double getMaxValue2(DataArrayInt*& tupleIds) const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT double getMinValue2(DataArrayInt*& tupleIds) const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT double getAverageValue() const throw(INTERP_KERNEL::Exception); - MEDCOUPLING_EXPORT void accumulate(double *res) const; - MEDCOUPLING_EXPORT double accumulate(int compId) const; + MEDCOUPLING_EXPORT void accumulate(double *res) const throw(INTERP_KERNEL::Exception); + MEDCOUPLING_EXPORT double accumulate(int compId) const throw(INTERP_KERNEL::Exception); + MEDCOUPLING_EXPORT DataArrayDouble *fromPolarToCart() const throw(INTERP_KERNEL::Exception); + MEDCOUPLING_EXPORT DataArrayDouble *fromCylToCart() const throw(INTERP_KERNEL::Exception); + MEDCOUPLING_EXPORT DataArrayDouble *fromSpherToCart() const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT DataArrayDouble *doublyContractedProduct() const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT DataArrayDouble *determinant() const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT DataArrayDouble *eigenValues() const throw(INTERP_KERNEL::Exception); @@ -165,12 +177,12 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT DataArrayDouble *magnitude() const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT DataArrayDouble *maxPerTuple() const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT void sortPerTuple(bool asc) throw(INTERP_KERNEL::Exception); - MEDCOUPLING_EXPORT void applyLin(double a, double b, int compoId); + MEDCOUPLING_EXPORT void applyLin(double a, double b, int compoId) throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT DataArrayDouble *applyFunc(int nbOfComp, FunctionToEvaluate func) const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT DataArrayDouble *applyFunc(int nbOfComp, const char *func) const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT DataArrayDouble *applyFunc(const char *func) const throw(INTERP_KERNEL::Exception); - MEDCOUPLING_EXPORT void applyFuncFast32(const char *func); - MEDCOUPLING_EXPORT void applyFuncFast64(const char *func); + MEDCOUPLING_EXPORT void applyFuncFast32(const char *func) throw(INTERP_KERNEL::Exception); + MEDCOUPLING_EXPORT void applyFuncFast64(const char *func) throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT DataArrayInt *getIdsInRange(double vmin, double vmax) const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT static DataArrayDouble *aggregate(const DataArrayDouble *a1, const DataArrayDouble *a2) throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT static DataArrayDouble *dot(const DataArrayDouble *a1, const DataArrayDouble *a2) throw(INTERP_KERNEL::Exception); @@ -197,6 +209,8 @@ namespace ParaMEDMEM { public: MEDCOUPLING_EXPORT static DataArrayInt *New(); + MEDCOUPLING_EXPORT bool isAllocated() const; + MEDCOUPLING_EXPORT void checkAllocated() const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT DataArrayInt *deepCopy() const; MEDCOUPLING_EXPORT DataArrayInt *performCpy(bool deepCpy) const; MEDCOUPLING_EXPORT void alloc(int nbOfTuple, int nbOfCompo); @@ -204,6 +218,7 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT bool isEqualWithoutConsideringStr(const DataArrayInt& other) const; MEDCOUPLING_EXPORT void fillWithZero(); MEDCOUPLING_EXPORT void fillWithValue(int val); + MEDCOUPLING_EXPORT void iota(int init=0) throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT std::string repr() const; MEDCOUPLING_EXPORT std::string reprZip() const; MEDCOUPLING_EXPORT void reprStream(std::ostream& stream) const; @@ -214,8 +229,10 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT DataArrayInt *invertArrayO2N2N2O(int newNbOfElem) const; MEDCOUPLING_EXPORT DataArrayInt *invertArrayN2O2O2N(int oldNbOfElem) const; //!alloc or useArray should have been called before. - MEDCOUPLING_EXPORT void reAlloc(int nbOfTuples); + MEDCOUPLING_EXPORT void reAlloc(int nbOfTuples) throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT DataArrayDouble *convertToDblArr() const; + MEDCOUPLING_EXPORT DataArrayInt *fromNoInterlace() const throw(INTERP_KERNEL::Exception); + MEDCOUPLING_EXPORT DataArrayInt *toNoInterlace() const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT void renumberInPlace(const int *old2New); MEDCOUPLING_EXPORT void renumberInPlaceR(const int *new2Old); MEDCOUPLING_EXPORT DataArrayInt *renumber(const int *old2New) const; @@ -223,6 +240,7 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT DataArrayInt *renumberAndReduce(const int *old2NewBg, int newNbOfTuple) const; MEDCOUPLING_EXPORT DataArrayInt *selectByTupleId(const int *new2OldBg, const int *new2OldEnd) const; MEDCOUPLING_EXPORT bool isIdentity() const; + MEDCOUPLING_EXPORT bool isUniform(int val) const; MEDCOUPLING_EXPORT DataArrayInt *substr(int tupleIdBg, int tupleIdEnd=-1) const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT DataArrayInt *changeNbOfComponents(int newNbOfComp, int dftValue) const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT DataArrayInt *keepSelectedComponents(const std::vector& compoIds) const throw(INTERP_KERNEL::Exception); @@ -233,6 +251,8 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT int *getPointer() const { return _mem.getPointer(); } MEDCOUPLING_EXPORT static void setArrayIn(DataArrayInt *newArray, DataArrayInt* &arrayToSet); MEDCOUPLING_EXPORT const int *getConstPointer() const { return _mem.getConstPointer(); } + MEDCOUPLING_EXPORT DataArrayInt *getIdsEqual(int val) const throw(INTERP_KERNEL::Exception); + MEDCOUPLING_EXPORT DataArrayInt *getIdsEqualList(const std::vector& vals) const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT static DataArrayInt *aggregate(const DataArrayInt *a1, const DataArrayInt *a2, int offsetA2); MEDCOUPLING_EXPORT static DataArrayInt *makePartition(const std::vector& groups, int newNb, std::vector< std::vector >& fidsOfGroups); MEDCOUPLING_EXPORT void useArray(const int *array, bool ownership, DeallocType type, int nbOfTuple, int nbOfCompo); diff --git a/src/MEDCoupling/MEDCouplingMemArray.txx b/src/MEDCoupling/MEDCouplingMemArray.txx index 5a1a048d5..1c6e8cefc 100644 --- a/src/MEDCoupling/MEDCouplingMemArray.txx +++ b/src/MEDCoupling/MEDCouplingMemArray.txx @@ -180,6 +180,32 @@ namespace ParaMEDMEM T *pt=_pointer.getPointer(); std::fill(pt,pt+_nb_of_elem,val); } + + template + T *MemArray::fromNoInterlace(int nbOfComp) const + { + const T *pt=_pointer.getConstPointer(); + int nbOfTuples=_nb_of_elem/nbOfComp; + T *ret=new T[_nb_of_elem]; + T *w=ret; + for(int i=0;i + T *MemArray::toNoInterlace(int nbOfComp) const + { + const T *pt=_pointer.getConstPointer(); + int nbOfTuples=_nb_of_elem/nbOfComp; + T *ret=new T[_nb_of_elem]; + T *w=ret; + for(int i=0;i void MemArray::alloc(int nbOfElements) diff --git a/src/MEDCoupling/MEDCouplingMesh.hxx b/src/MEDCoupling/MEDCouplingMesh.hxx index 524c23017..75ed3f309 100644 --- a/src/MEDCoupling/MEDCouplingMesh.hxx +++ b/src/MEDCoupling/MEDCouplingMesh.hxx @@ -91,6 +91,7 @@ namespace ParaMEDMEM virtual MEDCouplingMesh *mergeMyselfWith(const MEDCouplingMesh *other) const = 0; virtual MEDCouplingMesh *buildPart(const int *start, const int *end) const = 0; virtual MEDCouplingMesh *buildPartAndReduceNodes(const int *start, const int *end, DataArrayInt*& arr) const = 0; + virtual DataArrayInt *simplexize(int policy) throw(INTERP_KERNEL::Exception) = 0; virtual bool areCompatibleForMerge(const MEDCouplingMesh *other) const; static MEDCouplingMesh *mergeMeshes(const MEDCouplingMesh *mesh1, const MEDCouplingMesh *mesh2); //serialisation-unserialization diff --git a/src/MEDCoupling/MEDCouplingPointSet.cxx b/src/MEDCoupling/MEDCouplingPointSet.cxx index 87964cbdf..e587b2978 100644 --- a/src/MEDCoupling/MEDCouplingPointSet.cxx +++ b/src/MEDCoupling/MEDCouplingPointSet.cxx @@ -209,7 +209,7 @@ void MEDCouplingPointSet::findCommonNodes(int limitNodeId, double prec, DataArra findCommonNodesAlg<1>(bbox,nbNodesOld,limitNodeId,prec,c,cI); break; default: - throw INTERP_KERNEL::Exception("Unexpected spacedim of coords. Must be 1,2 or 3."); + throw INTERP_KERNEL::Exception("Unexpected spacedim of coords. Must be 1, 2 or 3."); } commIndex->alloc(cI.size(),1); std::copy(cI.begin(),cI.end(),commIndex->getPointer()); @@ -217,6 +217,55 @@ void MEDCouplingPointSet::findCommonNodes(int limitNodeId, double prec, DataArra std::copy(c.begin(),c.end(),comm->getPointer()); } +std::vector MEDCouplingPointSet::getNodeIdsNearPoint(const double *pos, double eps) const throw(INTERP_KERNEL::Exception) +{ + std::vector c,cI; + getNodeIdsNearPoints(pos,1,eps,c,cI); + return c; +} + +/*! + * Given a point given by its position 'pos' this method finds the set of node ids that are a a distance lower than eps. + * Position 'pos' is expected to be of size getSpaceDimension(). If not the behabiour is not warranted. + * This method throws an exception if no coordiantes are set. + */ +void MEDCouplingPointSet::getNodeIdsNearPoints(const double *pos, int nbOfNodes, double eps, std::vector& c, std::vector& cI) const throw(INTERP_KERNEL::Exception) +{ + if(!_coords) + throw INTERP_KERNEL::Exception("MEDCouplingPointSet::getNodeIdsNearPoint : no coordiantes set !"); + const double *coordsPtr=_coords->getConstPointer(); + if(!coordsPtr) + throw INTERP_KERNEL::Exception("MEDCouplingPointSet::getNodeIdsNearPoint : coordiante array set but no data inside it !"); + int spaceDim=getSpaceDimension(); + int nbNodes=getNumberOfNodes(); + std::vector bbox(2*nbNodes*spaceDim); + for(int i=0;i ret; + c.clear(); + cI.resize(1); cI[0]=0; + switch(spaceDim) + { + case 3: + findNodeIdsNearPointAlg<3>(bbox,pos,nbOfNodes,eps,c,cI); + break; + case 2: + findNodeIdsNearPointAlg<2>(bbox,pos,nbOfNodes,eps,c,cI); + break; + case 1: + findNodeIdsNearPointAlg<1>(bbox,pos,nbOfNodes,eps,c,cI); + break; + default: + throw INTERP_KERNEL::Exception("Unexpected spacedim of coords for getNodeIdsNearPoint. Must be 1, 2 or 3."); + } +} + /*! * @param comm in param in the same format than one returned by findCommonNodes method. * @param commI in param in the same format than one returned by findCommonNodes method. @@ -228,33 +277,28 @@ DataArrayInt *MEDCouplingPointSet::buildNewNumberingFromCommonNodesFormat(const DataArrayInt *ret=DataArrayInt::New(); int nbNodesOld=getNumberOfNodes(); ret->alloc(nbNodesOld,1); - std::fill(ret->getPointer(),ret->getPointer()+nbNodesOld,-1); - int *retPtr=ret->getPointer(); - std::vector commRemain(comm->getConstPointer(),comm->getConstPointer()+comm->getNumberOfTuples()); - std::vector commIRemain(commIndex->getConstPointer(),commIndex->getConstPointer()+commIndex->getNumberOfTuples()); + int *pt=ret->getPointer(); + std::fill(pt,pt+nbNodesOld,-1); + int nbOfGrps=commIndex->getNumberOfTuples()-1; + const int *cIPtr=commIndex->getPointer(); + const int *cPtr=comm->getPointer(); + for(int i=0;i::const_iterator iNode2=commRemain.begin(); - iNode2!=commRemain.begin()+commIRemain[1];iNode2++) - retPtr[*iNode2]=newNb; - int delta=commIRemain[1]; - commRemain.erase(commRemain.begin(),commRemain.begin()+commIRemain[1]); - commIRemain.erase(commIRemain.begin()); - std::transform(commIRemain.begin(),commIRemain.end(),commIRemain.begin(),std::bind2nd(std::minus(),delta)); + if(pt[iNode]==-1) + pt[iNode]=newNb++; + else + { + int grpId=-(pt[iNode]+2); + for(int j=cIPtr[grpId];j getNodeIdsNearPoint(const double *pos, double eps) const throw(INTERP_KERNEL::Exception); + void getNodeIdsNearPoints(const double *pos, int nbOfNodes, double eps, std::vector& c, std::vector& cI) const throw(INTERP_KERNEL::Exception); void findCommonNodes(int limitNodeId, double prec, DataArrayInt *&comm, DataArrayInt *&commIndex) const; DataArrayInt *buildNewNumberingFromCommonNodesFormat(const DataArrayInt *comm, const DataArrayInt *commIndex, int& newNbOfNodes) const; @@ -101,6 +103,9 @@ namespace ParaMEDMEM template void findCommonNodesAlg(std::vector& bbox, int nbNodes, int limitNodeId, double prec, std::vector& c, std::vector& cI) const; + template + void findNodeIdsNearPointAlg(std::vector& bbox, const double *pos, int nbNodes, double eps, + std::vector& c, std::vector& cI) const; protected: DataArrayDouble *_coords; }; diff --git a/src/MEDCoupling/MEDCouplingPointSet.txx b/src/MEDCoupling/MEDCouplingPointSet.txx index bae195566..9374a1cb6 100644 --- a/src/MEDCoupling/MEDCouplingPointSet.txx +++ b/src/MEDCoupling/MEDCouplingPointSet.txx @@ -36,32 +36,63 @@ namespace ParaMEDMEM BBTree myTree(&bbox[0],0,0,nbNodes,-prec); double bb[2*SPACEDIM]; double prec2=prec*prec; + std::vector isDone(nbNodes); + for(int i=0;i intersectingElems; + myTree.getIntersectingElems(bb,intersectingElems); + if(intersectingElems.size()>1) + { + std::vector commonNodes; + for(std::vector::const_iterator it=intersectingElems.begin();it!=intersectingElems.end();it++) + if(*it!=i) + if(*it>=limitNodeId) + if(INTERP_KERNEL::distance2(coordsPtr+SPACEDIM*i,coordsPtr+SPACEDIM*(*it)) + void MEDCouplingPointSet::findNodeIdsNearPointAlg(std::vector& bbox, const double *pos, int nbNodes, double eps, + std::vector& c, std::vector& cI) const + { + const double *coordsPtr=_coords->getConstPointer(); + BBTree myTree(&bbox[0],0,0,getNumberOfNodes(),-eps); + double bb[2*SPACEDIM]; + double eps2=eps*eps; for(int i=0;i intersectingElems; myTree.getIntersectingElems(bb,intersectingElems); - if(intersectingElems.size()>1) - { - std::vector commonNodes; - for(std::vector::const_iterator it=intersectingElems.begin();it!=intersectingElems.end();it++) - if(*it!=i) - if(*it>=limitNodeId) - if(INTERP_KERNEL::distance2(coordsPtr+SPACEDIM*i,coordsPtr+SPACEDIM*(*it)) commonNodes; + for(std::vector::const_iterator it=intersectingElems.begin();it!=intersectingElems.end();it++) + if(INTERP_KERNEL::distance2(pos+SPACEDIM*i,coordsPtr+SPACEDIM*(*it)) @@ -563,6 +565,59 @@ void MEDCouplingUMesh::convertToPolyTypes(const std::vector& cellIdsToConve computeTypes(); } +/*! + * This method is the opposite of ParaMEDMEM::MEDCouplingUMesh::convertToPolyTypes method. + * The aim is to take all polygons or polyhedrons cell and to try to traduce them into classical cells. + * + */ +void MEDCouplingUMesh::unPolyze() +{ + checkFullyDefined(); + if(getMeshDimension()<=1) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::unPolyze works on umeshes with meshdim equals to 2 or 3 !"); + int nbOfCells=getNumberOfCells(); + if(nbOfCells<1) + return ; + int initMeshLgth=getMeshLength(); + int *conn=_nodal_connec->getPointer(); + int *index=_nodal_connec_index->getPointer(); + int posOfCurCell=0; + int newPos=0; + int lgthOfCurCell; + for(int i=0;ireAlloc(newPos); + computeTypes(); +} + /*! * Array returned is the correspondance new to old. * The maximum value stored in returned array is the number of nodes of 'this' minus 1 after call of this method. @@ -1294,7 +1349,13 @@ std::string MEDCouplingUMesh::simpleRepr() const ret << "Unstructured mesh with name : \"" << getName() << "\"\n"; ret << "Mesh dimension : " << _mesh_dim << "\nSpace dimension : "; if(_coords!=0) - ret << getSpaceDimension() << "\n"; + { + const int spaceDim=getSpaceDimension(); + ret << spaceDim << "\nInfo attached on space dimension : "; + for(int i=0;igetInfoOnComponent(i) << "\" "; + ret << "\n"; + } else ret << msg0 << "\n"; ret << "Number of nodes : "; @@ -1732,6 +1793,32 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::buildDirectionVectorField() const return ret; } +/*! + * This method checks that 'this' is a contiguous mesh. The user is expected to call this method on a mesh with meshdim==1. + * If not an exception will thrown. If this is an empty mesh with no cell an exception will be thrown too. + * No consideration of coordinate is done by this method. + * A 1D mesh is said contiguous if : a cell i with nodal connectivity (k,p) the cell i+1 the nodal connectivity should be (p,m) + * If not false is returned. In case that false is returned a call to ParaMEDMEM::MEDCouplingUMesh::mergeNodes could be usefull. + */ +bool MEDCouplingUMesh::isContiguous1D() const throw(INTERP_KERNEL::Exception) +{ + if(getMeshDimension()!=1) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::isContiguous1D : this method has a sense only for 1D mesh !"); + int nbCells=getNumberOfCells(); + if(nbCells<1) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::isContiguous1D : this method has a sense for non empty mesh !"); + const int *connI=_nodal_connec_index->getConstPointer(); + const int *conn=_nodal_connec->getConstPointer(); + int ref=conn[connI[0]+2]; + for(int i=1;icheckFullyDefined(); - if(getMeshDimension()!=2 || getSpaceDimension()!=3) - throw INTERP_KERNEL::Exception("Invalid 'this' for buildExtrudedMeshFromThis method : must be meshDim==2 and spaceDim ==3 !"); - if(mesh1D->getMeshDimension()!=1 || mesh1D->getSpaceDimension()!=3) - throw INTERP_KERNEL::Exception("Invalid 'mesh1D' for buildExtrudedMeshFromThis method : must be meshDim==1 and spaceDim ==3 !"); + if(!mesh1D->isContiguous1D()) + throw INTERP_KERNEL::Exception("buildExtrudedMeshFromThis : 1D mesh passed in parameter is not contiguous !"); + if(getSpaceDimension()!=mesh1D->getSpaceDimension()) + throw INTERP_KERNEL::Exception("Invalid call to buildExtrudedMeshFromThis this and mesh1D must have same dimension !"); + if((getMeshDimension()!=2 || getSpaceDimension()!=3) && (getMeshDimension()!=1 || getSpaceDimension()!=2)) + throw INTERP_KERNEL::Exception("Invalid 'this' for buildExtrudedMeshFromThis method : must be (meshDim==2 and spaceDim==3) or (meshDim==1 and spaceDim==2) !"); + if(mesh1D->getMeshDimension()!=1) + throw INTERP_KERNEL::Exception("Invalid 'mesh1D' for buildExtrudedMeshFromThis method : must be meshDim==1 !"); bool isQuad=false; if(isPresenceOfQuadratic()) { @@ -1957,6 +2045,11 @@ MEDCouplingUMesh *MEDCouplingUMesh::buildExtrudedMeshFromThis(const MEDCouplingU newCoords=fillExtCoordsUsingTranslation(mesh1D,isQuad); break; } + case 1: + { + newCoords=fillExtCoordsUsingTranslAndAutoRotation(mesh1D,isQuad); + break; + } default: throw INTERP_KERNEL::Exception("Not implemented extrusion policy : must be in (0) !"); } @@ -1976,13 +2069,14 @@ DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslation(const MEDCoupli { int oldNbOfNodes=getNumberOfNodes(); int nbOf1DCells=mesh1D->getNumberOfCells(); + int spaceDim=getSpaceDimension(); DataArrayDouble *ret=DataArrayDouble::New(); std::vector isQuads; int nbOfLevsInVec=isQuad?2*nbOf1DCells+1:nbOf1DCells+1; - ret->alloc(oldNbOfNodes*nbOfLevsInVec,3); + ret->alloc(oldNbOfNodes*nbOfLevsInVec,spaceDim); double *retPtr=ret->getPointer(); const double *coords=getCoords()->getConstPointer(); - double *work=std::copy(coords,coords+3*oldNbOfNodes,retPtr); + double *work=std::copy(coords,coords+spaceDim*oldNbOfNodes,retPtr); std::vector v; std::vector c; double vec[3]; @@ -1995,23 +2089,152 @@ DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslation(const MEDCoupli c.resize(0); mesh1D->getCoordinatesOfNode(v[isQuad?2:1],c); mesh1D->getCoordinatesOfNode(v[0],c); - std::transform(c.begin(),c.begin()+3,c.begin()+3,vec,std::minus()); + std::transform(c.begin(),c.begin()+spaceDim,c.begin()+spaceDim,vec,std::minus()); for(int j=0;j()); + work=std::transform(vec,vec+spaceDim,retPtr+spaceDim*(i*oldNbOfNodes+j),work,std::plus()); if(isQuad) { c.resize(0); mesh1D->getCoordinatesOfNode(v[1],c); mesh1D->getCoordinatesOfNode(v[0],c); - std::transform(c.begin(),c.begin()+3,c.begin()+3,vec,std::minus()); + std::transform(c.begin(),c.begin()+spaceDim,c.begin()+spaceDim,vec,std::minus()); for(int j=0;j()); + work=std::transform(vec,vec+spaceDim,retPtr+spaceDim*(i*oldNbOfNodes+j),work,std::plus()); } } ret->copyStringInfoFrom(*getCoords()); return ret; } +/*! + * This method incarnates the policy 1 for MEDCouplingUMesh::buildExtrudedMeshFromThis method. + * @param mesh1D is the input 1D mesh used for translation and automatic rotation computation. + * @return newCoords new coords filled by this method. + */ +DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation(const MEDCouplingUMesh *mesh1D, bool isQuad) const throw(INTERP_KERNEL::Exception) +{ + if(mesh1D->getSpaceDimension()==2) + return fillExtCoordsUsingTranslAndAutoRotation2D(mesh1D,isQuad); + if(mesh1D->getSpaceDimension()==3) + return fillExtCoordsUsingTranslAndAutoRotation3D(mesh1D,isQuad); + throw INTERP_KERNEL::Exception("Not implemented rotation and translation alg. for spacedim other than 2 and 3 !"); +} + +/*! + * This method incarnates the policy 1 for MEDCouplingUMesh::buildExtrudedMeshFromThis method. + * @param mesh1D is the input 1D mesh used for translation and automatic rotation computation. + * @return newCoords new coords filled by this method. + */ +DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation2D(const MEDCouplingUMesh *mesh1D, bool isQuad) const throw(INTERP_KERNEL::Exception) +{ + if(isQuad) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation2D : not implemented for quadratic cells !"); + int oldNbOfNodes=getNumberOfNodes(); + int nbOf1DCells=mesh1D->getNumberOfCells(); + if(nbOf1DCells<2) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation2D : impossible to detect any angle of rotation ! Change extrusion policy 1->0 !"); + DataArrayDouble *ret=DataArrayDouble::New(); + int nbOfLevsInVec=nbOf1DCells+1; + ret->alloc(oldNbOfNodes*nbOfLevsInVec,2); + double *retPtr=ret->getPointer(); + retPtr=std::copy(getCoords()->getConstPointer(),getCoords()->getConstPointer()+getCoords()->getNbOfElems(),retPtr); + MEDCouplingUMesh *tmp=MEDCouplingUMesh::New(); + DataArrayDouble *tmp2=getCoords()->deepCopy(); + tmp->setCoords(tmp2); + tmp2->decrRef(); + const double *coo1D=mesh1D->getCoords()->getConstPointer(); + const int *conn1D=mesh1D->getNodalConnectivity()->getConstPointer(); + const int *connI1D=mesh1D->getNodalConnectivityIndex()->getConstPointer(); + for(int i=1;itranslate(vec); + double tmp3[2],radius,alpha,alpha0; + const double *p0=i+1rotate(end,0,angle); + retPtr=std::copy(tmp2->getConstPointer(),tmp2->getConstPointer()+tmp2->getNbOfElems(),retPtr); + } + tmp->decrRef(); + return ret; +} + +/*! + * This method incarnates the policy 1 for MEDCouplingUMesh::buildExtrudedMeshFromThis method. + * @param mesh1D is the input 1D mesh used for translation and automatic rotation computation. + * @return newCoords new coords filled by this method. + */ +DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation3D(const MEDCouplingUMesh *mesh1D, bool isQuad) const throw(INTERP_KERNEL::Exception) +{ + if(isQuad) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation3D : not implemented for quadratic cells !"); + int oldNbOfNodes=getNumberOfNodes(); + int nbOf1DCells=mesh1D->getNumberOfCells(); + if(nbOf1DCells<2) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation3D : impossible to detect any angle of rotation ! Change extrusion policy 1->0 !"); + DataArrayDouble *ret=DataArrayDouble::New(); + int nbOfLevsInVec=nbOf1DCells+1; + ret->alloc(oldNbOfNodes*nbOfLevsInVec,3); + double *retPtr=ret->getPointer(); + retPtr=std::copy(getCoords()->getConstPointer(),getCoords()->getConstPointer()+getCoords()->getNbOfElems(),retPtr); + MEDCouplingUMesh *tmp=MEDCouplingUMesh::New(); + DataArrayDouble *tmp2=getCoords()->deepCopy(); + tmp->setCoords(tmp2); + tmp2->decrRef(); + const double *coo1D=mesh1D->getCoords()->getConstPointer(); + const int *conn1D=mesh1D->getNodalConnectivity()->getConstPointer(); + const int *connI1D=mesh1D->getNodalConnectivityIndex()->getConstPointer(); + for(int i=1;itranslate(vec); + double tmp3[2],radius,alpha,alpha0; + const double *p0=i+11.e-7) + { + vecPlane[0]/=norm; vecPlane[1]/=norm; vecPlane[2]/=norm; + double norm2=sqrt(vecPlane[0]*vecPlane[0]+vecPlane[1]*vecPlane[1]); + double vec2[2]={vecPlane[1]/norm2,-vecPlane[0]/norm2}; + double s2=norm2; + double c2=cos(asin(s2)); + double m[3][3]={ + {vec2[0]*vec2[0]*(1-c2)+c2, vec2[0]*vec2[1]*(1-c2), vec2[1]*s2}, + {vec2[0]*vec2[1]*(1-c2), vec2[1]*vec2[1]*(1-c2)+c2, -vec2[0]*s2}, + {-vec2[1]*s2, vec2[0]*s2, c2} + }; + double p0r[3]={m[0][0]*p0[0]+m[0][1]*p0[1]+m[0][2]*p0[2], m[1][0]*p0[0]+m[1][1]*p0[1]+m[1][2]*p0[2], m[2][0]*p0[0]+m[2][1]*p0[1]+m[2][2]*p0[2]}; + double p1r[3]={m[0][0]*p1[0]+m[0][1]*p1[1]+m[0][2]*p1[2], m[1][0]*p1[0]+m[1][1]*p1[1]+m[1][2]*p1[2], m[2][0]*p1[0]+m[2][1]*p1[1]+m[2][2]*p1[2]}; + double p2r[3]={m[0][0]*p2[0]+m[0][1]*p2[1]+m[0][2]*p2[2], m[1][0]*p2[0]+m[1][1]*p2[1]+m[1][2]*p2[2], m[2][0]*p2[0]+m[2][1]*p2[1]+m[2][2]*p2[2]}; + INTERP_KERNEL::EdgeArcCircle::getArcOfCirclePassingThru(p0r,p1r,p2r,tmp3,radius,alpha,alpha0); + double cosangle=i+1rotate(end,vecPlane,angle); + + } + retPtr=std::copy(tmp2->getConstPointer(),tmp2->getConstPointer()+tmp2->getNbOfElems(),retPtr); + } + tmp->decrRef(); + return ret; +} + /*! * This method is private because not easy to use for end user. This method is const contrary to * MEDCouplingUMesh::buildExtrudedMeshFromThis method because this->_coords are expected to contain @@ -2023,7 +2246,7 @@ MEDCouplingUMesh *MEDCouplingUMesh::buildExtrudedMeshFromThisLowLev(int nbOfNode int nbOf1DCells=getNumberOfNodes()/nbOfNodesOf1Lev-1; int nbOf2DCells=getNumberOfCells(); int nbOf3DCells=nbOf2DCells*nbOf1DCells; - MEDCouplingUMesh *ret=MEDCouplingUMesh::New("Extruded",3); + MEDCouplingUMesh *ret=MEDCouplingUMesh::New("Extruded",getMeshDimension()+1); const int *conn=_nodal_connec->getConstPointer(); const int *connI=_nodal_connec_index->getConstPointer(); DataArrayInt *newConn=DataArrayInt::New(); @@ -2156,6 +2379,198 @@ void MEDCouplingUMesh::convertQuadraticCellsToLinear() throw(INTERP_KERNEL::Exce setConnectivity(newConn,newConnI,false); } +/*! + * This methods modify this by converting each cells into simplex cell, that is too say triangle for meshdim==2 or tetra for meshdim==3. + * This cut into simplex is performed following the parameter 'policy'. This method so typically increases the number of cells of this. + * This method returns new2old array that specifies a each cell of 'this' after the call what was its id it comes. + * + * The semantic of 'policy' parameter : + * - 1 only QUAD4. For QUAD4 the cut is done along 0-2 diagonal for QUAD4 + * - 2 only QUAD4. For QUAD4 the cut is done along 1-3 diagonal for QUAD4 + */ +DataArrayInt *MEDCouplingUMesh::simplexize(int policy) throw(INTERP_KERNEL::Exception) +{ + switch(policy) + { + case 0: + return simplexizePol0(); + case 1: + return simplexizePol1(); + default: + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::simplexize : unrecognized policy ! Must be 0 or 1 !"); + } +} + +bool MEDCouplingUMesh::areOnlySimplexCells() const throw(INTERP_KERNEL::Exception) +{ + checkFullyDefined(); + if(getMeshDimension()<1) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::areOnlySimplexCells : only available with meshes having a meshdim >= 1 !"); + int nbCells=getNumberOfCells(); + const int *conn=_nodal_connec->getConstPointer(); + const int *connI=_nodal_connec_index->getConstPointer(); + for(int i=0;ialloc(nbOfCells+nbOfCutCells,1); + if(nbOfCutCells==0) + { + ret->iota(0); + return ret; + } + int *retPt=ret->getPointer(); + DataArrayInt *newConn=DataArrayInt::New(); + DataArrayInt *newConnI=DataArrayInt::New(); + newConnI->alloc(nbOfCells+nbOfCutCells+1,1); + newConn->alloc(getMeshLength()+3*nbOfCutCells,1); + int *pt=newConn->getPointer(); + int *ptI=newConnI->getPointer(); + ptI[0]=0; + const int *oldc=_nodal_connec->getConstPointer(); + const int *ci=_nodal_connec_index->getConstPointer(); + for(int i=0;idecrRef(); + _nodal_connec=newConn; + _nodal_connec_index->decrRef(); + _nodal_connec_index=newConnI; + computeTypes(); + updateTime(); + return ret; +} + +/*! + * This method implements policy 1 of virtual method ParaMEDMEM::MEDCouplingUMesh::simplexize. + */ +DataArrayInt *MEDCouplingUMesh::simplexizePol1() throw(INTERP_KERNEL::Exception) +{ + if(getMeshDimension()!=2) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::simplexizePol0 : this policy is only available for mesh with meshdim == 2 !"); + int nbOfCells=getNumberOfCells(); + DataArrayInt *ret=DataArrayInt::New(); + int nbOfCutCells=getNumberOfCellsWithType(INTERP_KERNEL::NORM_QUAD4); + ret->alloc(nbOfCells+nbOfCutCells,1); + if(nbOfCutCells==0) + { + ret->iota(0); + return ret; + } + int *retPt=ret->getPointer(); + DataArrayInt *newConn=DataArrayInt::New(); + DataArrayInt *newConnI=DataArrayInt::New(); + newConnI->alloc(nbOfCells+nbOfCutCells+1,1); + newConn->alloc(getMeshLength()+3*nbOfCutCells,1); + int *pt=newConn->getPointer(); + int *ptI=newConnI->getPointer(); + ptI[0]=0; + const int *oldc=_nodal_connec->getConstPointer(); + const int *ci=_nodal_connec_index->getConstPointer(); + for(int i=0;idecrRef(); + _nodal_connec=newConn; + _nodal_connec_index->decrRef(); + _nodal_connec_index=newConnI; + computeTypes(); + updateTime(); + return ret; +} + + +/*! + * This method converts all degenerated cells to simpler cells. For example a NORM_QUAD4 cell consituted from 2 same node id in its + * nodal connectivity will be transform to a NORM_TRI3 cell. + * This method works \b only \b on \b linear cells. + * This method works on nodes ids, that is to say a call to ParaMEDMEM::MEDCouplingUMesh::mergeNodes + * method could be usefull before calling this method in case of presence of several pair of nodes located on same position. + * This method throws an exception if 'this' is not fully defined (connectivity). + * This method throws an exception too if a "too" degenerated cell is detected. For example a NORM_TRI3 with 3 times the same node id. + */ +void MEDCouplingUMesh::convertDegeneratedCells() throw(INTERP_KERNEL::Exception) +{ + checkFullyDefined(); + if(getMeshDimension()<=1) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::convertDegeneratedCells works on umeshes with meshdim equals to 2 or 3 !"); + int nbOfCells=getNumberOfCells(); + if(nbOfCells<1) + return ; + int initMeshLgth=getMeshLength(); + int *conn=_nodal_connec->getPointer(); + int *index=_nodal_connec_index->getPointer(); + int posOfCurCell=0; + int newPos=0; + int lgthOfCurCell; + for(int i=0;ireAlloc(newPos); + computeTypes(); +} + /*! * This method checks that all or only polygons (depending 'polyOnly' parameter) 2D cells are correctly oriented relative to 'vec' vector. * The 'vec' vector has to have a non nul norm. diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx index eb989c9a6..df4ddd292 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh.hxx @@ -79,6 +79,7 @@ namespace ParaMEDMEM 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 void unPolyze(); MEDCOUPLING_EXPORT DataArrayInt *zipCoordsTraducer(); MEDCOUPLING_EXPORT DataArrayInt *zipConnectivityTraducer(int compType); MEDCOUPLING_EXPORT void getReverseNodalConnectivity(DataArrayInt *revNodal, DataArrayInt *revNodalIndx) const; @@ -98,6 +99,7 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT MEDCouplingFieldDouble *getMeasureFieldOnNode(bool isAbs) const; MEDCOUPLING_EXPORT MEDCouplingFieldDouble *buildOrthogonalField() const; MEDCOUPLING_EXPORT MEDCouplingFieldDouble *buildDirectionVectorField() const; + MEDCOUPLING_EXPORT bool isContiguous1D() const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT void project1D(const double *pt, const double *v, double eps, double *res) const; MEDCOUPLING_EXPORT int getCellContainingPoint(const double *pos, double eps) const; MEDCOUPLING_EXPORT void getCellsContainingPoint(const double *pos, double eps, std::vector& elts) const; @@ -108,6 +110,9 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT bool isFullyQuadratic() const; MEDCOUPLING_EXPORT bool isPresenceOfQuadratic() const; MEDCOUPLING_EXPORT void convertQuadraticCellsToLinear() throw(INTERP_KERNEL::Exception); + MEDCOUPLING_EXPORT DataArrayInt *simplexize(int policy) throw(INTERP_KERNEL::Exception); + MEDCOUPLING_EXPORT bool areOnlySimplexCells() const throw(INTERP_KERNEL::Exception); + MEDCOUPLING_EXPORT void convertDegeneratedCells() throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT void are2DCellsNotCorrectlyOriented(const double *vec, bool polyOnly, std::vector& cells) const throw(INTERP_KERNEL::Exception); 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); @@ -142,10 +147,15 @@ namespace ParaMEDMEM void checkFullyDefined() const throw(INTERP_KERNEL::Exception); void reprConnectivityOfThisLL(std::ostringstream& stream) const; //tools + DataArrayInt *simplexizePol0() throw(INTERP_KERNEL::Exception); + DataArrayInt *simplexizePol1() throw(INTERP_KERNEL::Exception); void renumberNodesInConn(const int *newNodeNumbers); void fillCellIdsToKeepFromNodeIds(const int *begin, const int *end, bool fullyIn, std::vector& cellIdsKept) const; MEDCouplingUMesh *buildExtrudedMeshFromThisLowLev(int nbOfNodesOf1Lev, bool isQuad) const; DataArrayDouble *fillExtCoordsUsingTranslation(const MEDCouplingUMesh *mesh1D, bool isQuad) const; + DataArrayDouble *fillExtCoordsUsingTranslAndAutoRotation(const MEDCouplingUMesh *mesh1D, bool isQuad) const throw(INTERP_KERNEL::Exception); + DataArrayDouble *fillExtCoordsUsingTranslAndAutoRotation2D(const MEDCouplingUMesh *mesh1D, bool isQuad) const throw(INTERP_KERNEL::Exception); + DataArrayDouble *fillExtCoordsUsingTranslAndAutoRotation3D(const MEDCouplingUMesh *mesh1D, bool isQuad) const throw(INTERP_KERNEL::Exception); template void findCommonCellsBase(int compType, std::vector& res, std::vector& resI) const; bool areCellsEqualInPool(const std::vector& candidates, int compType, std::vector& result) const; diff --git a/src/MEDCoupling/MEDCouplingUMeshDesc.cxx b/src/MEDCoupling/MEDCouplingUMeshDesc.cxx index f45cd5eba..e8629b7b8 100644 --- a/src/MEDCoupling/MEDCouplingUMeshDesc.cxx +++ b/src/MEDCoupling/MEDCouplingUMeshDesc.cxx @@ -361,6 +361,11 @@ MEDCouplingPointSet *MEDCouplingUMeshDesc::buildFacePartOfMySelfNode(const int * return 0; } +DataArrayInt *MEDCouplingUMeshDesc::simplexize(int policy) throw(INTERP_KERNEL::Exception) +{ + throw INTERP_KERNEL::Exception("MEDCouplingUMeshDesc::simplexize : Not implemented yet !"); +} + void MEDCouplingUMeshDesc::findBoundaryNodes(std::vector& nodes) const { //not implemented yet diff --git a/src/MEDCoupling/MEDCouplingUMeshDesc.hxx b/src/MEDCoupling/MEDCouplingUMeshDesc.hxx index 7c965a31c..f11fe94a4 100644 --- a/src/MEDCoupling/MEDCouplingUMeshDesc.hxx +++ b/src/MEDCoupling/MEDCouplingUMeshDesc.hxx @@ -66,6 +66,7 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT MEDCouplingPointSet *buildPartOfMySelf(const int *start, const int *end, bool keepCoords) const; MEDCOUPLING_EXPORT MEDCouplingPointSet *buildPartOfMySelfNode(const int *start, const int *end, bool fullyIn) const; MEDCOUPLING_EXPORT MEDCouplingPointSet *buildFacePartOfMySelfNode(const int *start, const int *end, bool fullyIn) const; + MEDCOUPLING_EXPORT DataArrayInt *simplexize(int policy) throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT void findBoundaryNodes(std::vector& nodes) const; MEDCOUPLING_EXPORT MEDCouplingPointSet *buildBoundaryMesh(bool keepCoords) const; MEDCOUPLING_EXPORT void renumberCells(const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception); diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx index 569560584..4982c2e7a 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx @@ -143,6 +143,23 @@ namespace ParaMEDMEM CPPUNIT_TEST( testDataArrayIntInvertO2NNO21 ); CPPUNIT_TEST( testKeepSetSelectedComponent1 ); CPPUNIT_TEST( testKeepSetSelectedComponent2 ); + CPPUNIT_TEST( testDAIGetIdsEqual1 ); + CPPUNIT_TEST( testDAIGetIdsEqualList1 ); + CPPUNIT_TEST( testDAFromNoInterlace1 ); + CPPUNIT_TEST( testDAToNoInterlace1 ); + CPPUNIT_TEST( testDAIsUniform1 ); + CPPUNIT_TEST( testDADFromPolarToCart1 ); + CPPUNIT_TEST( testDADFromCylToCart1 ); + CPPUNIT_TEST( testDADFromSpherToCart1 ); + CPPUNIT_TEST( testUnPolyze1 ); + CPPUNIT_TEST( testConvertDegeneratedCells1 ); + CPPUNIT_TEST( testGetNodeIdsNearPoints1 ); + CPPUNIT_TEST( testFieldCopyTinyAttrFrom1 ); + CPPUNIT_TEST( testExtrudedMesh5 ); + CPPUNIT_TEST( testExtrudedMesh6 ); + CPPUNIT_TEST( testExtrudedMesh7 ); + CPPUNIT_TEST( testSimplexize1 ); + CPPUNIT_TEST( testSimplexize2 ); //MEDCouplingBasicsTestInterp.cxx CPPUNIT_TEST( test2DInterpP0P0_1 ); CPPUNIT_TEST( test2DInterpP0P0PL_1 ); @@ -314,6 +331,23 @@ namespace ParaMEDMEM void testDataArrayIntInvertO2NNO21(); void testKeepSetSelectedComponent1(); void testKeepSetSelectedComponent2(); + void testDAIGetIdsEqual1(); + void testDAIGetIdsEqualList1(); + void testDAFromNoInterlace1(); + void testDAToNoInterlace1(); + void testDAIsUniform1(); + void testDADFromPolarToCart1(); + void testDADFromCylToCart1(); + void testDADFromSpherToCart1(); + void testUnPolyze1(); + void testConvertDegeneratedCells1(); + void testGetNodeIdsNearPoints1(); + void testFieldCopyTinyAttrFrom1(); + void testExtrudedMesh5(); + void testExtrudedMesh6(); + void testExtrudedMesh7(); + void testSimplexize1(); + void testSimplexize2(); //MEDCouplingBasicsTestInterp.cxx void test2DInterpP0P0_1(); void test2DInterpP0P0PL_1(); diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest2.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest2.cxx index 0c36bdbad..2b2f4de2a 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest2.cxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest2.cxx @@ -3005,3 +3005,616 @@ void MEDCouplingBasicsTest::testKeepSetSelectedComponent2() a1->decrRef(); m1->decrRef(); } + +void MEDCouplingBasicsTest::testDAIGetIdsEqual1() +{ + const int tab1[7]={5,-2,-4,-2,3,2,-2}; + DataArrayInt *da=DataArrayInt::New(); + da->alloc(7,1); + std::copy(tab1,tab1+7,da->getPointer()); + DataArrayInt *da2=da->getIdsEqual(-2); + CPPUNIT_ASSERT_EQUAL(3,da2->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(1,da2->getNumberOfComponents()); + const int expected1[3]={1,3,6}; + CPPUNIT_ASSERT(std::equal(expected1,expected1+3,da2->getConstPointer())); + da2->decrRef(); + da->decrRef(); +} + +void MEDCouplingBasicsTest::testDAIGetIdsEqualList1() +{ + const int tab1[7]={5,-2,-4,-2,3,2,-2}; + DataArrayInt *da=DataArrayInt::New(); + da->alloc(7,1); + std::copy(tab1,tab1+7,da->getPointer()); + const int tab2[3]={3,-2,0}; + std::vector tab2V(tab2,tab2+3); + DataArrayInt *da2=da->getIdsEqualList(tab2V); + CPPUNIT_ASSERT_EQUAL(4,da2->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(1,da2->getNumberOfComponents()); + const int expected1[4]={1,3,4,6}; + CPPUNIT_ASSERT(std::equal(expected1,expected1+4,da2->getConstPointer())); + da2->decrRef(); + da->decrRef(); +} + +void MEDCouplingBasicsTest::testDAFromNoInterlace1() +{ + const int tab1[15]={1,11,21,31,41,2,12,22,32,42,3,13,23,33,43}; + DataArrayInt *da=DataArrayInt::New(); + da->alloc(5,3); + std::copy(tab1,tab1+15,da->getPointer()); + DataArrayInt *da2=da->fromNoInterlace(); + const int expected1[15]={1,2,3,11,12,13,21,22,23,31,32,33,41,42,43}; + CPPUNIT_ASSERT_EQUAL(5,da2->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(3,da2->getNumberOfComponents());// it's not a bug. Avoid to have 1 million components ! + CPPUNIT_ASSERT(std::equal(expected1,expected1+15,da2->getConstPointer())); + DataArrayDouble *da3=da->convertToDblArr(); + DataArrayDouble *da4=da3->fromNoInterlace(); + CPPUNIT_ASSERT_EQUAL(5,da4->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(3,da4->getNumberOfComponents());// it's not a bug. Avoid to have 1 million components ! + for(int i=0;i<15;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL((double)expected1[i],da4->getIJ(0,i),1e-14); + da4->decrRef(); + da3->decrRef(); + da2->decrRef(); + da->decrRef(); +} + +void MEDCouplingBasicsTest::testDAToNoInterlace1() +{ + const int tab1[15]={1,2,3,11,12,13,21,22,23,31,32,33,41,42,43}; + DataArrayInt *da=DataArrayInt::New(); + da->alloc(5,3); + std::copy(tab1,tab1+15,da->getPointer()); + DataArrayInt *da2=da->toNoInterlace(); + const int expected1[15]={1,11,21,31,41,2,12,22,32,42,3,13,23,33,43}; + CPPUNIT_ASSERT_EQUAL(5,da2->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(3,da2->getNumberOfComponents());// it's not a bug. Avoid to have 1 million components ! + CPPUNIT_ASSERT(std::equal(expected1,expected1+15,da2->getConstPointer())); + DataArrayDouble *da3=da->convertToDblArr(); + DataArrayDouble *da4=da3->toNoInterlace(); + CPPUNIT_ASSERT_EQUAL(5,da4->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(3,da4->getNumberOfComponents());// it's not a bug. Avoid to have 1 million components ! + for(int i=0;i<15;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL((double)expected1[i],da4->getIJ(0,i),1e-14); + da4->decrRef(); + da3->decrRef(); + da2->decrRef(); + da->decrRef(); +} + +void MEDCouplingBasicsTest::testDAIsUniform1() +{ + const int tab1[5]={1,1,1,1,1}; + DataArrayInt *da=DataArrayInt::New(); + da->alloc(5,1); + std::copy(tab1,tab1+5,da->getPointer()); + CPPUNIT_ASSERT(da->isUniform(1)); + da->setIJ(2,0,2); + CPPUNIT_ASSERT(!da->isUniform(1)); + da->setIJ(2,0,1); + CPPUNIT_ASSERT(da->isUniform(1)); + DataArrayDouble *da2=da->convertToDblArr(); + CPPUNIT_ASSERT(da2->isUniform(1.,1e-12)); + da2->setIJ(1,0,1.+1.e-13); + CPPUNIT_ASSERT(da2->isUniform(1.,1e-12)); + da2->setIJ(1,0,1.+1.e-11); + CPPUNIT_ASSERT(!da2->isUniform(1.,1e-12)); + da2->decrRef(); + da->decrRef(); +} + +void MEDCouplingBasicsTest::testDADFromPolarToCart1() +{ + const double tab1[4]={2.,0.2,2.5,0.7}; + DataArrayDouble *da=DataArrayDouble::New(); + da->alloc(2,2); + std::copy(tab1,tab1+4,da->getPointer()); + DataArrayDouble *da2=da->fromPolarToCart(); + const double expected1[4]={1.9601331556824833,0.39733866159012243, 1.9121054682112213,1.6105442180942275}; + for(int i=0;i<4;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected1[i],da2->getIJ(0,i),1e-13); + da2->decrRef(); + da->decrRef(); +} + +void MEDCouplingBasicsTest::testDADFromCylToCart1() +{ + const double tab1[6]={2.,0.2,4.,2.5,0.7,9.}; + DataArrayDouble *da=DataArrayDouble::New(); + da->alloc(2,3); + std::copy(tab1,tab1+6,da->getPointer()); + DataArrayDouble *da2=da->fromCylToCart(); + const double expected1[6]={1.9601331556824833,0.39733866159012243,4., 1.9121054682112213,1.6105442180942275,9.}; + for(int i=0;i<6;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected1[i],da2->getIJ(0,i),1e-13); + da2->decrRef(); + da->decrRef(); +} + +void MEDCouplingBasicsTest::testDADFromSpherToCart1() +{ + const double tab1[6]={2.,0.2,0.3,2.5,0.7,0.8}; + DataArrayDouble *da=DataArrayDouble::New(); + da->alloc(2,3); + std::copy(tab1,tab1+6,da->getPointer()); + DataArrayDouble *da2=da->fromSpherToCart(); + const double expected1[6]={0.37959212195737485,0.11742160338765303,1.9601331556824833, 1.1220769624465328,1.1553337045129035,1.9121054682112213}; + for(int i=0;i<6;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected1[i],da2->getIJ(0,i),1e-13); + da2->decrRef(); + da->decrRef(); +} + +void MEDCouplingBasicsTest::testUnPolyze1() +{ + const int elts[8]={0,1,2,3,4,5,6,7}; + std::vector eltsV(elts,elts+8); + MEDCouplingUMesh *mesh=build3DTargetMesh_1(); + mesh->convertToPolyTypes(eltsV); + mesh->unPolyze(); + MEDCouplingUMesh *mesh2=build3DTargetMesh_1(); + mesh->checkCoherency(); + CPPUNIT_ASSERT(mesh->isEqual(mesh2,1e-12)); + mesh->convertToPolyTypes(eltsV); + CPPUNIT_ASSERT(!mesh->isEqual(mesh2,1e-12)); + mesh->getNodalConnectivity()->setIJ(0,6,10); + mesh->getNodalConnectivity()->setIJ(0,7,9); + mesh->getNodalConnectivity()->setIJ(0,8,12); + mesh->getNodalConnectivity()->setIJ(0,9,13); + mesh->unPolyze(); + CPPUNIT_ASSERT(mesh->isEqual(mesh2,1e-12)); + mesh->convertToPolyTypes(eltsV); + mesh->getNodalConnectivity()->setIJ(0,6,12); + mesh->getNodalConnectivity()->setIJ(0,7,13); + mesh->getNodalConnectivity()->setIJ(0,8,10); + mesh->getNodalConnectivity()->setIJ(0,9,9); + mesh->unPolyze(); + CPPUNIT_ASSERT(mesh->isEqual(mesh2,1e-12)); + mesh->convertToPolyTypes(eltsV); + mesh->getNodalConnectivity()->setIJ(0,6,12); + mesh->getNodalConnectivity()->setIJ(0,7,10); + mesh->getNodalConnectivity()->setIJ(0,8,13); + mesh->getNodalConnectivity()->setIJ(0,9,9); + mesh->unPolyze(); + CPPUNIT_ASSERT(!mesh->isEqual(mesh2,1e-12)); + mesh->decrRef(); + mesh2->decrRef(); + // Test for 2D mesh + mesh=build2DTargetMesh_1(); + mesh2=build2DTargetMesh_1(); + eltsV.resize(5); + mesh->convertToPolyTypes(eltsV); + CPPUNIT_ASSERT(!mesh->isEqual(mesh2,1e-12)); + mesh->unPolyze(); + CPPUNIT_ASSERT(mesh->isEqual(mesh2,1e-12)); + mesh->decrRef(); + mesh2->decrRef(); +} + +void MEDCouplingBasicsTest::testConvertDegeneratedCells1() +{ + MEDCouplingUMesh *mesh=build3DTargetMesh_1(); + int conn[32]={0,1,3,3,9,10,12,12, 0,1,3,4,9,9,9,9, 1,1,1,1,10,12,9,10, 10,11,12,9,1,1,1,1}; + mesh->allocateCells(4); + mesh->insertNextCell(INTERP_KERNEL::NORM_HEXA8,8,conn); + mesh->insertNextCell(INTERP_KERNEL::NORM_HEXA8,8,conn+8); + mesh->insertNextCell(INTERP_KERNEL::NORM_HEXA8,8,conn+16); + mesh->insertNextCell(INTERP_KERNEL::NORM_HEXA8,8,conn+24); + mesh->finishInsertingCells(); + mesh->checkCoherency(); + CPPUNIT_ASSERT_EQUAL(4,mesh->getNumberOfCells()); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_HEXA8,mesh->getTypeOfCell(0)); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_HEXA8,mesh->getTypeOfCell(1)); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_HEXA8,mesh->getTypeOfCell(2)); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_HEXA8,mesh->getTypeOfCell(3)); + MEDCouplingFieldDouble *f1=mesh->getMeasureField(true); + mesh->convertDegeneratedCells(); + mesh->checkCoherency(); + MEDCouplingFieldDouble *f2=mesh->getMeasureField(true); + CPPUNIT_ASSERT_EQUAL(4,mesh->getNumberOfCells()); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_PENTA6,mesh->getTypeOfCell(0)); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_PYRA5,mesh->getTypeOfCell(1)); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_TETRA4,mesh->getTypeOfCell(2)); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_PYRA5,mesh->getTypeOfCell(3)); + for(int i=0;i<4;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(f1->getArray()->getIJ(0,i),f2->getArray()->getIJ(0,i),1e-5); + f1->decrRef(); + f2->decrRef(); + mesh->decrRef(); +} + +void MEDCouplingBasicsTest::testGetNodeIdsNearPoints1() +{ + MEDCouplingUMesh *mesh=build2DTargetMesh_1(); + DataArrayDouble *coords=mesh->getCoords(); + DataArrayDouble *tmp=DataArrayDouble::New(); + tmp->alloc(3,2); + const double vals[6]={0.2,0.2,0.1,0.2,0.2,0.2}; + std::copy(vals,vals+6,tmp->getPointer()); + DataArrayDouble *tmp2=DataArrayDouble::aggregate(coords,tmp); + tmp->decrRef(); + mesh->setCoords(tmp2); + tmp2->decrRef(); + const double pts[6]={0.2,0.2,0.1,0.3,-0.3,0.7}; + std::vector c=mesh->getNodeIdsNearPoint(pts,1e-7); + CPPUNIT_ASSERT_EQUAL(3,(int)c.size()); + CPPUNIT_ASSERT_EQUAL(4,c[0]); + CPPUNIT_ASSERT_EQUAL(9,c[1]); + CPPUNIT_ASSERT_EQUAL(11,c[2]); + c.clear(); + std::vector cI; + mesh->getNodeIdsNearPoints(pts,3,1e-7,c,cI); + CPPUNIT_ASSERT_EQUAL(4,(int)cI.size()); + CPPUNIT_ASSERT_EQUAL(4,(int)c.size()); + CPPUNIT_ASSERT_EQUAL(4,c[0]); + CPPUNIT_ASSERT_EQUAL(9,c[1]); + CPPUNIT_ASSERT_EQUAL(11,c[2]); + CPPUNIT_ASSERT_EQUAL(6,c[3]); + CPPUNIT_ASSERT_EQUAL(0,cI[0]); + CPPUNIT_ASSERT_EQUAL(3,cI[1]); + CPPUNIT_ASSERT_EQUAL(3,cI[2]); + CPPUNIT_ASSERT_EQUAL(4,cI[3]); + mesh->decrRef(); +} + +void MEDCouplingBasicsTest::testFieldCopyTinyAttrFrom1() +{ + MEDCouplingFieldDouble *f1=MEDCouplingFieldDouble::New(ON_CELLS,ONE_TIME); + f1->setName("f1"); + f1->setTimeTolerance(1.e-5); + f1->setDescription("f1Desc"); + f1->setTime(1.23,4,5); + MEDCouplingFieldDouble *f2=MEDCouplingFieldDouble::New(ON_CELLS,ONE_TIME); + f2->setName("f2"); + f2->setDescription("f2Desc"); + f2->setTime(6.78,9,10); + f2->setTimeTolerance(4.556e-12); + // + int dt,it; + f1->copyTinyAttrFrom(f2); + CPPUNIT_ASSERT_DOUBLES_EQUAL(4.556e-12,f1->getTimeTolerance(),1e-24); + CPPUNIT_ASSERT_DOUBLES_EQUAL(6.78,f1->getTime(dt,it),1e-12); + CPPUNIT_ASSERT_EQUAL(9,dt); + CPPUNIT_ASSERT_EQUAL(10,it); + CPPUNIT_ASSERT(std::string(f1->getName())=="f1");//name unchanged + CPPUNIT_ASSERT(std::string(f1->getDescription())=="f1Desc");//description unchanged + f1->decrRef(); + f2->decrRef(); + // + f1=MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME); + f1->setName("f1"); + f1->setTimeTolerance(1.e-5); + f1->setDescription("f1Desc"); + f2=MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME); + f2->setName("f2"); + f2->setDescription("f2Desc"); + f2->setTimeTolerance(4.556e-12); + // + f1->copyTinyAttrFrom(f2); + CPPUNIT_ASSERT_DOUBLES_EQUAL(4.556e-12,f1->getTimeTolerance(),1e-24); + CPPUNIT_ASSERT(std::string(f1->getName())=="f1");//name unchanged + CPPUNIT_ASSERT(std::string(f1->getDescription())=="f1Desc");//description unchanged + f1->decrRef(); + f2->decrRef(); + // + f1=MEDCouplingFieldDouble::New(ON_CELLS,CONST_ON_TIME_INTERVAL); + f1->setName("f1"); + f1->setTimeTolerance(1.e-5); + f1->setDescription("f1Desc"); + f1->setTime(1.23,4,5); + f1->setEndTime(5.43,2,1); + f2=MEDCouplingFieldDouble::New(ON_CELLS,CONST_ON_TIME_INTERVAL); + f2->setName("f2"); + f2->setDescription("f2Desc"); + f2->setTimeTolerance(4.556e-12); + f2->setTime(6.78,9,10); + f2->setEndTime(10.98,7,6); + // + f1->copyTinyAttrFrom(f2); + CPPUNIT_ASSERT_DOUBLES_EQUAL(4.556e-12,f1->getTimeTolerance(),1e-24); + CPPUNIT_ASSERT(std::string(f1->getName())=="f1");//name unchanged + CPPUNIT_ASSERT(std::string(f1->getDescription())=="f1Desc");//description unchanged + CPPUNIT_ASSERT_DOUBLES_EQUAL(6.78,f1->getTime(dt,it),1e-12); + CPPUNIT_ASSERT_EQUAL(9,dt); + CPPUNIT_ASSERT_EQUAL(10,it); + CPPUNIT_ASSERT_DOUBLES_EQUAL(10.98,f1->getEndTime(dt,it),1e-12); + CPPUNIT_ASSERT_EQUAL(7,dt); + CPPUNIT_ASSERT_EQUAL(6,it); + f1->decrRef(); + f2->decrRef(); + // + f1=MEDCouplingFieldDouble::New(ON_CELLS,LINEAR_TIME); + f1->setName("f1"); + f1->setTimeTolerance(1.e-5); + f1->setDescription("f1Desc"); + f1->setTime(1.23,4,5); + f1->setEndTime(5.43,2,1); + f2=MEDCouplingFieldDouble::New(ON_CELLS,LINEAR_TIME); + f2->setName("f2"); + f2->setDescription("f2Desc"); + f2->setTimeTolerance(4.556e-12); + f2->setTime(6.78,9,10); + f2->setEndTime(10.98,7,6); + // + f1->copyTinyAttrFrom(f2); + CPPUNIT_ASSERT_DOUBLES_EQUAL(4.556e-12,f1->getTimeTolerance(),1e-24); + CPPUNIT_ASSERT(std::string(f1->getName())=="f1");//name unchanged + CPPUNIT_ASSERT(std::string(f1->getDescription())=="f1Desc");//description unchanged + CPPUNIT_ASSERT_DOUBLES_EQUAL(6.78,f1->getTime(dt,it),1e-12); + CPPUNIT_ASSERT_EQUAL(9,dt); + CPPUNIT_ASSERT_EQUAL(10,it); + CPPUNIT_ASSERT_DOUBLES_EQUAL(10.98,f1->getEndTime(dt,it),1e-12); + CPPUNIT_ASSERT_EQUAL(7,dt); + CPPUNIT_ASSERT_EQUAL(6,it); + f1->decrRef(); + f2->decrRef(); +} + +/*! + * 1D -> 2D extrusion with rotation + */ +void MEDCouplingBasicsTest::testExtrudedMesh5() +{ + const double coo1[4]={0.,1.,2.,3.5}; + DataArrayDouble *a=DataArrayDouble::New(); + a->alloc(4,1); + std::copy(coo1,coo1+4,a->getPointer()); + MEDCouplingCMesh *b=MEDCouplingCMesh::New(); + b->setCoordsAt(0,a); + MEDCouplingUMesh *c=b->buildUnstructured(); + CPPUNIT_ASSERT_EQUAL(1,c->getSpaceDimension()); + c->changeSpaceDimension(2); + // + DataArrayDouble *d=DataArrayDouble::New(); + d->alloc(13,1); + d->iota(); + MEDCouplingCMesh *e=MEDCouplingCMesh::New(); + e->setCoordsAt(0,d); + MEDCouplingUMesh *f=e->buildUnstructured(); + DataArrayDouble *g=f->getCoords()->applyFunc(2,"3.5*IVec+x/6*3.14159265359*JVec"); + DataArrayDouble *h=g->fromPolarToCart(); + f->setCoords(h); + MEDCouplingUMesh *i=c->buildExtrudedMeshFromThis(f,1); + CPPUNIT_ASSERT_EQUAL(52,i->getNumberOfNodes()); + bool tmp2; + int tmp3; + DataArrayInt *tmp=i->mergeNodes(1e-9,tmp2,tmp3); + CPPUNIT_ASSERT(tmp2); + CPPUNIT_ASSERT_EQUAL(37,tmp3); + tmp->decrRef(); + i->convertDegeneratedCells(); + i->checkCoherency(); + CPPUNIT_ASSERT_EQUAL(36,i->getNumberOfCells()); + CPPUNIT_ASSERT_EQUAL(37,i->getNumberOfNodes()); + CPPUNIT_ASSERT_EQUAL(12,i->getNumberOfCellsWithType(INTERP_KERNEL::NORM_TRI3)); + CPPUNIT_ASSERT_EQUAL(24,i->getNumberOfCellsWithType(INTERP_KERNEL::NORM_QUAD4)); + const double expected1[3]={0.25,0.75,2.0625}; + MEDCouplingFieldDouble *j=i->getMeasureField(true); + for(int i=0;i<12;i++) + for(int k=0;k<3;k++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected1[k],j->getIJ(0,i*3+k),1e-10); + const double expected2[72]={0.62200846792814113, 0.16666666666681595, 1.4513530918323276, 0.38888888888923495, 2.6293994326053212, 0.7045454545460802, 0.45534180126145435, 0.45534180126150181, 1.0624642029433926, 1.0624642029435025, 1.9248539780597826, 1.9248539780599816, 0.16666666666661334, 0.62200846792815856, 0.38888888888876294, 1.4513530918323678, 0.70454545454522521, 2.629399432605394, -0.16666666666674007, 0.62200846792812436, -0.38888888888906142, 1.4513530918322881, -0.70454545454576778, 2.6293994326052488, -0.45534180126154766, 0.45534180126140844, -1.0624642029436118, 1.0624642029432834, -1.9248539780601803, 1.9248539780595841, -0.62200846792817499, 0.1666666666665495, -1.451353091832408, 0.388888888888613, -2.6293994326054668, 0.70454545454495332, -0.62200846792810593, -0.16666666666680507, -1.451353091832247, -0.38888888888921297, -2.6293994326051746, -0.70454545454604123, -0.45534180126135926, -0.45534180126159562, -1.0624642029431723, -1.0624642029437235, -1.9248539780593836, -1.9248539780603811, -0.1666666666664828, -0.62200846792819242, -0.38888888888846079, -1.4513530918324489, -0.70454545454467987, -2.6293994326055397, 0.16666666666687083, -0.62200846792808862, 0.38888888888936374, -1.4513530918322073, 0.70454545454631357, -2.6293994326051022, 0.45534180126164348, -0.45534180126131207, 1.0624642029438327, -1.0624642029430627, 1.9248539780605791, -1.9248539780591853, 0.62200846792821063, -0.16666666666641802, 1.4513530918324888, -0.38888888888831086, 2.6293994326056125, -0.70454545454440853}; + DataArrayDouble *m=i->getBarycenterAndOwner(); + for(int i=0;i<72;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected2[i],m->getIJ(0,i),1e-10); + // + m->decrRef(); + j->decrRef(); + i->decrRef(); + h->decrRef(); + g->decrRef(); + f->decrRef(); + e->decrRef(); + d->decrRef(); + c->decrRef(); + b->decrRef(); + a->decrRef(); +} + +/*! + * 1D -> 2D extrusion without rotation + */ +void MEDCouplingBasicsTest::testExtrudedMesh6() +{ + const double coo1[4]={0.,1.,2.,3.5}; + DataArrayDouble *a=DataArrayDouble::New(); + a->alloc(4,1); + std::copy(coo1,coo1+4,a->getPointer()); + MEDCouplingCMesh *b=MEDCouplingCMesh::New(); + b->setCoordsAt(0,a); + MEDCouplingUMesh *c=b->buildUnstructured(); + CPPUNIT_ASSERT_EQUAL(1,c->getSpaceDimension()); + c->changeSpaceDimension(2); + // + DataArrayDouble *d=DataArrayDouble::New(); + d->alloc(5,1); + d->iota(); + MEDCouplingCMesh *e=MEDCouplingCMesh::New(); + e->setCoordsAt(0,d); + MEDCouplingUMesh *f=e->buildUnstructured(); + DataArrayDouble *d2=f->getCoords()->applyFunc("x*x/2"); + f->setCoords(d2); + f->changeSpaceDimension(2); + // + const double center[2]={0.,0.}; + f->rotate(center,0,M_PI/3); + MEDCouplingUMesh *g=c->buildExtrudedMeshFromThis(f,0); + g->checkCoherency(); + const double expected1[]={ 0.4330127018922193, 0.4330127018922193, 0.649519052838329, 1.2990381056766578, 1.299038105676658, 1.948557158514987, 2.1650635094610955, 2.1650635094610964, 3.2475952641916446, 3.031088913245533, 3.0310889132455352, 4.546633369868303 }; + MEDCouplingFieldDouble *f1=g->getMeasureField(true); + for(int i=0;i<12;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected1[i],f1->getIJ(0,i),1e-12); + + const double expected2[]={0.625, 0.21650635094610962, 1.625, 0.21650635094610959, 2.8750000000000004, 0.21650635094610965, 1.1250000000000002, 1.0825317547305482, 2.125, 1.0825317547305482, 3.3750000000000004, 1.0825317547305484, 2.125, 2.8145825622994254, 3.125, 2.8145825622994254, 4.375, 2.8145825622994254, 3.6250000000000009, 5.4126587736527414, 4.625, 5.4126587736527414, 5.875, 5.4126587736527414}; + DataArrayDouble *f2=g->getBarycenterAndOwner(); + for(int i=0;i<24;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected2[i],f2->getIJ(0,i),1e-12); + // + f1->decrRef(); + f2->decrRef(); + g->decrRef(); + f->decrRef(); + e->decrRef(); + d->decrRef(); + d2->decrRef(); + c->decrRef(); + b->decrRef(); + a->decrRef(); +} + +/*! + * 2D -> 3D extrusion with rotation + */ +void MEDCouplingBasicsTest::testExtrudedMesh7() +{ + const double coo1[4]={0.,1.,2.,3.5}; + DataArrayDouble *a=DataArrayDouble::New(); + a->alloc(4,1); + std::copy(coo1,coo1+4,a->getPointer()); + MEDCouplingCMesh *b=MEDCouplingCMesh::New(); + b->setCoordsAt(0,a); + MEDCouplingUMesh *c=b->buildUnstructured(); + CPPUNIT_ASSERT_EQUAL(1,c->getSpaceDimension()); + c->changeSpaceDimension(2); + // + DataArrayDouble *d=DataArrayDouble::New(); + d->alloc(13,1); + d->iota(); + MEDCouplingCMesh *e=MEDCouplingCMesh::New(); + e->setCoordsAt(0,d); + MEDCouplingUMesh *f=e->buildUnstructured(); + DataArrayDouble *g=f->getCoords()->applyFunc(2,"3.5*IVec+x/6*3.14159265359*JVec"); + DataArrayDouble *h=g->fromPolarToCart(); + f->setCoords(h); + MEDCouplingUMesh *i=c->buildExtrudedMeshFromThis(f,1); + CPPUNIT_ASSERT_EQUAL(52,i->getNumberOfNodes()); + bool tmp2; + int tmp3; + DataArrayInt *tmp=i->mergeNodes(1e-9,tmp2,tmp3); + CPPUNIT_ASSERT(tmp2); + CPPUNIT_ASSERT_EQUAL(37,tmp3); + tmp->decrRef(); + i->convertDegeneratedCells(); + const double vec1[3]={10.,0.,0.}; + i->translate(vec1); + DataArrayDouble *g2=h->applyFunc(3,"13.5/3.5*x*IVec+0*JVec+13.5/3.5*y*KVec"); + f->setCoords(g2); + i->changeSpaceDimension(3); + MEDCouplingUMesh *i3=i->buildExtrudedMeshFromThis(f,1); + MEDCouplingFieldDouble *f2=i3->getMeasureField(true); + tmp=i->mergeNodes(1e-9,tmp2,tmp3); + CPPUNIT_ASSERT(tmp2); + CPPUNIT_ASSERT_EQUAL(444,tmp3); + tmp->decrRef(); + const double expected1[36]={1.327751058489274, 4.2942574094314701, 13.024068164857139, 1.3069177251569044, 4.1484240761012954, 12.297505664866796, 1.270833333332571, 3.8958333333309674, 11.039062499993179, 1.2291666666659207, 3.6041666666644425, 9.585937499993932, 1.1930822748415895, 3.3515759238941376, 8.3274943351204556, 1.1722489415082769, 3.2057425905609289, 7.6009318351210622, 1.1722489415082862, 3.2057425905609884, 7.6009318351213713, 1.1930822748416161, 3.3515759238943001, 8.3274943351212727, 1.2291666666659564, 3.6041666666646734, 9.5859374999950777, 1.2708333333326081, 3.8958333333311868, 11.039062499994293, 1.3069177251569224, 4.1484240761014384, 12.297505664867627, 1.3277510584902354, 4.2942574094346071, 13.024068164866796}; + int kk=0; + for(int ii=0;ii<12;ii++) + for(int jj=0;jj<36;jj++,kk++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected1[jj],f2->getIJ(0,kk),1e-9); + // + f2->decrRef(); + i3->decrRef(); + g2->decrRef(); + i->decrRef(); + h->decrRef(); + g->decrRef(); + f->decrRef(); + e->decrRef(); + d->decrRef(); + c->decrRef(); + b->decrRef(); + a->decrRef(); +} + +void MEDCouplingBasicsTest::testSimplexize1() +{ + MEDCouplingUMesh *m=build3DSurfTargetMesh_1(); + std::vector v(1); + v[0]=3; + m->convertToPolyTypes(v); + DataArrayInt *da=m->simplexize(0); + CPPUNIT_ASSERT_EQUAL(7,da->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(1,da->getNumberOfComponents()); + const int expected2[7]={0,0,1,2,3,4,4}; + for(int i=0;i<7;i++) + CPPUNIT_ASSERT_EQUAL(expected2[i],da->getIJ(i,0)); + m->checkCoherency(); + CPPUNIT_ASSERT_EQUAL(7,m->getNumberOfCells()); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_TRI3,m->getTypeOfCell(0)); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_TRI3,m->getTypeOfCell(1)); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_TRI3,m->getTypeOfCell(2)); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_TRI3,m->getTypeOfCell(3)); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_POLYGON,m->getTypeOfCell(4)); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_TRI3,m->getTypeOfCell(5)); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_TRI3,m->getTypeOfCell(6)); + const double expected1[7]={0.125,0.125,0.125,0.125,0.25,0.125,0.125}; + MEDCouplingFieldDouble *f=m->getMeasureField(false); + for(int i=0;i<7;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected1[i]*sqrt(2.),f->getIJ(i,0),1e-10); + std::set types=m->getAllTypes(); + CPPUNIT_ASSERT_EQUAL(2,(int)types.size()); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_TRI3,*(types.begin())); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_POLYGON,*(++(types.begin()))); + f->decrRef(); + da->decrRef(); + m->decrRef(); + // + m=build3DSurfTargetMesh_1(); + v[0]=3; + m->convertToPolyTypes(v); + da=m->simplexize(1); + CPPUNIT_ASSERT_EQUAL(7,da->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(1,da->getNumberOfComponents()); + for(int i=0;i<7;i++) + CPPUNIT_ASSERT_EQUAL(expected2[i],da->getIJ(i,0)); + m->checkCoherency(); + types=m->getAllTypes(); + CPPUNIT_ASSERT_EQUAL(2,(int)types.size()); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_TRI3,*(types.begin())); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_POLYGON,*(++(types.begin()))); + CPPUNIT_ASSERT_EQUAL(7,m->getNumberOfCells()); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_TRI3,m->getTypeOfCell(0)); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_TRI3,m->getTypeOfCell(1)); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_TRI3,m->getTypeOfCell(2)); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_TRI3,m->getTypeOfCell(3)); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_POLYGON,m->getTypeOfCell(4)); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_TRI3,m->getTypeOfCell(5)); + CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_TRI3,m->getTypeOfCell(6)); + f=m->getMeasureField(false); + for(int i=0;i<7;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected1[i]*sqrt(2.),f->getIJ(i,0),1e-10); + f->decrRef(); + da->decrRef(); + m->decrRef(); +} + +void MEDCouplingBasicsTest::testSimplexize2() +{ + MEDCouplingUMesh *m=build3DSurfTargetMesh_1(); + std::vector v(1); + v[0]=3; + m->convertToPolyTypes(v); + MEDCouplingFieldDouble *f1=MEDCouplingFieldDouble::New(ON_CELLS,ONE_TIME); + f1->setMesh(m); + DataArrayDouble *arr=DataArrayDouble::New(); + const double arr1[10]={10.,110.,20.,120.,30.,130.,40.,140.,50.,150.}; + arr->alloc(5,2); + std::copy(arr1,arr1+10,arr->getPointer()); + f1->setArray(arr); + arr->decrRef(); + // + f1->checkCoherency(); + CPPUNIT_ASSERT(f1->simplexize(0)); + f1->checkCoherency(); + const double expected1[14]={10.,110.,10.,110.,20.,120.,30.,130.,40.,140.,50.,150.,50.,150.}; + for(int i=0;i<14;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected1[i],f1->getIJ(0,i),1e-10); + CPPUNIT_ASSERT(!f1->simplexize(0)); + for(int i=0;i<14;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected1[i],f1->getIJ(0,i),1e-10); + // + f1->decrRef(); + m->decrRef(); +} diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py index a8e328c08..25119adda 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py @@ -354,6 +354,8 @@ class MEDCouplingBasicsTest(unittest.TestCase): tab1=[0,4] tab2=[0,2,3] # + subMesh=mesh.buildPart(tab1) + self.assertTrue(isinstance(subMesh,MEDCouplingUMesh)) subMesh=mesh.buildPartOfMySelf(tab1,True); self.assertTrue(isinstance(subMesh,MEDCouplingUMesh)) name=subMesh.getName(); @@ -4087,6 +4089,9 @@ class MEDCouplingBasicsTest(unittest.TestCase): pos1=[5.,30.,2.] self.assertEqual(16,m.getCellContainingPoint(pos1,1e-12)); # + elems=m2.giveElemsInBoundingBox([3.5,6.,12.2,25.,0.,1.5],1e-7) + self.assertEqual([1, 2, 4, 5, 7, 8, 10, 11, 13, 14, 16, 17],elems) + # pt=[2.4,12.7,-3.4] m.scale(pt,3.7); m3=m.buildUnstructured(); @@ -4274,15 +4279,15 @@ class MEDCouplingBasicsTest(unittest.TestCase): -0.4012053603397987, 0.8423032781211455, -0.3599436712889738#eigenvect 2 ] for i in xrange(5): - self.assertAlmostEqual(expected1[0],f2.getIJ(i,0),1e-13); - self.assertAlmostEqual(expected1[1],f2.getIJ(i,1),1e-13); - self.assertAlmostEqual(expected1[2],f2.getIJ(i,2),1e-13); - self.assertAlmostEqual(expected1[3],f2.getIJ(i,3),1e-13); - self.assertAlmostEqual(expected1[4],f2.getIJ(i,4),1e-13); - self.assertAlmostEqual(expected1[5],f2.getIJ(i,5),1e-13); - self.assertAlmostEqual(expected1[6],f2.getIJ(i,6),1e-13); - self.assertAlmostEqual(expected1[7],f2.getIJ(i,7),1e-13); - self.assertAlmostEqual(expected1[8],f2.getIJ(i,8),1e-13); + self.assertAlmostEqual(expected1[0],f2.getIJ(i,0),13); + self.assertAlmostEqual(expected1[1],f2.getIJ(i,1),13); + self.assertAlmostEqual(expected1[2],f2.getIJ(i,2),13); + self.assertAlmostEqual(expected1[3],f2.getIJ(i,3),13); + self.assertAlmostEqual(expected1[4],f2.getIJ(i,4),13); + self.assertAlmostEqual(expected1[5],f2.getIJ(i,5),13); + self.assertAlmostEqual(expected1[6],f2.getIJ(i,6),13); + self.assertAlmostEqual(expected1[7],f2.getIJ(i,7),13); + self.assertAlmostEqual(expected1[8],f2.getIJ(i,8),13); pass # pass @@ -4678,6 +4683,488 @@ class MEDCouplingBasicsTest(unittest.TestCase): pass # pass + + def testDAIGetIdsEqual1(self): + tab1=[5,-2,-4,-2,3,2,-2]; + da=DataArrayInt.New(); + da.setValues(tab1,7,1); + da2=da.getIdsEqual(-2); + self.assertEqual(3,da2.getNumberOfTuples()); + self.assertEqual(1,da2.getNumberOfComponents()); + expected1=[1,3,6]; + self.assertEqual(expected1,da2.getValues()); + pass + + def testDAIGetIdsEqualList1(self): + tab1=[5,-2,-4,-2,3,2,-2]; + da=DataArrayInt.New(); + da.setValues(tab1,7,1); + da2=da.getIdsEqualList([3,-2,0]); + self.assertEqual(4,da2.getNumberOfTuples()); + self.assertEqual(1,da2.getNumberOfComponents()); + expected1=[1,3,4,6]; + self.assertEqual(expected1,da2.getValues()); + pass + + def testDAFromNoInterlace1(self): + tab1=[1,11,21,31,41,2,12,22,32,42,3,13,23,33,43] + da=DataArrayInt.New(); + da.setValues(tab1,5,3); + da2=da.fromNoInterlace(); + expected1=[1,2,3,11,12,13,21,22,23,31,32,33,41,42,43] + self.assertEqual(5,da2.getNumberOfTuples()); + self.assertEqual(3,da2.getNumberOfComponents());# it's not a bug. Avoid to have 1 million components ! + self.assertEqual(expected1,da2.getValues()); + da3=da.convertToDblArr(); + da4=da3.fromNoInterlace(); + self.assertEqual(5,da4.getNumberOfTuples()); + self.assertEqual(3,da4.getNumberOfComponents());# it's not a bug. Avoid to have 1 million components ! + for i in xrange(15): + self.assertAlmostEqual(expected1[i],da4.getIJ(0,i),14); + pass + pass + + def testDAToNoInterlace1(self): + tab1=[1,2,3,11,12,13,21,22,23,31,32,33,41,42,43] + da=DataArrayInt.New(); + da.setValues(tab1,5,3); + da2=da.toNoInterlace(); + expected1=[1,11,21,31,41,2,12,22,32,42,3,13,23,33,43] + self.assertEqual(5,da2.getNumberOfTuples()); + self.assertEqual(3,da2.getNumberOfComponents());# it's not a bug. Avoid to have 1 million components ! + self.assertEqual(expected1,da2.getValues()); + da3=da.convertToDblArr(); + da4=da3.toNoInterlace(); + self.assertEqual(5,da4.getNumberOfTuples()); + self.assertEqual(3,da4.getNumberOfComponents());# it's not a bug. Avoid to have 1 million components ! + for i in xrange(15): + self.assertAlmostEqual(expected1[i],da4.getIJ(0,i),14); + pass + pass + + def testDAIsUniform1(self): + tab1=[1,1,1,1,1] + da=DataArrayInt.New(); + da.setValues(tab1,5,1); + self.assertTrue(da.isUniform(1)); + da.setIJ(2,0,2); + self.assertTrue(not da.isUniform(1)); + da.setIJ(2,0,1); + self.assertTrue(da.isUniform(1)); + da2=da.convertToDblArr(); + self.assertTrue(da2.isUniform(1.,1.e-12)); + da2.setIJ(1,0,1.+1.e-13); + self.assertTrue(da2.isUniform(1.,1.e-12)); + da2.setIJ(1,0,1.+1.e-11); + self.assertTrue(not da2.isUniform(1.,1.e-12)); + pass + + def testDADFromPolarToCart1(self): + tab1=[2.,0.2,2.5,0.7] + da=DataArrayDouble.New(); + da.setValues(tab1,2,2); + da2=da.fromPolarToCart(); + expected1=[1.9601331556824833,0.39733866159012243, 1.9121054682112213,1.6105442180942275] + for i in xrange(4): + self.assertAlmostEqual(expected1[i],da2.getIJ(0,i),13); + pass + pass + + def testDADFromCylToCart1(self): + tab1=[2.,0.2,4.,2.5,0.7,9.] + da=DataArrayDouble.New(); + da.setValues(tab1,2,3); + da2=da.fromCylToCart(); + expected1=[1.9601331556824833,0.39733866159012243,4., 1.9121054682112213,1.6105442180942275,9.] + for i in xrange(6): + self.assertAlmostEqual(expected1[i],da2.getIJ(0,i),13); + pass + pass + + def testDADFromSpherToCart1(self): + tab1=[2.,0.2,0.3,2.5,0.7,0.8] + da=DataArrayDouble.New(); + da.setValues(tab1,2,3); + da2=da.fromSpherToCart(); + expected1=[0.37959212195737485,0.11742160338765303,1.9601331556824833, 1.1220769624465328,1.1553337045129035,1.9121054682112213] + for i in xrange(6): + self.assertAlmostEqual(expected1[i],da2.getIJ(0,i),13); + pass + pass + + def testUnPolyze1(self): + elts=[0,1,2,3,4,5,6,7] + eltsV=elts; + mesh=MEDCouplingDataForTest.build3DTargetMesh_1(); + mesh.convertToPolyTypes(eltsV); + mesh.unPolyze(); + mesh2=MEDCouplingDataForTest.build3DTargetMesh_1(); + mesh.checkCoherency(); + self.assertTrue(mesh.isEqual(mesh2,1e-12)); + mesh.convertToPolyTypes(eltsV); + self.assertTrue(not mesh.isEqual(mesh2,1e-12)); + mesh.getNodalConnectivity().setIJ(0,6,10); + mesh.getNodalConnectivity().setIJ(0,7,9); + mesh.getNodalConnectivity().setIJ(0,8,12); + mesh.getNodalConnectivity().setIJ(0,9,13); + mesh.unPolyze(); + self.assertTrue(mesh.isEqual(mesh2,1e-12)); + mesh.convertToPolyTypes(eltsV); + mesh.getNodalConnectivity().setIJ(0,6,12); + mesh.getNodalConnectivity().setIJ(0,7,13); + mesh.getNodalConnectivity().setIJ(0,8,10); + mesh.getNodalConnectivity().setIJ(0,9,9); + mesh.unPolyze(); + self.assertTrue(mesh.isEqual(mesh2,1e-12)); + mesh.convertToPolyTypes(eltsV); + mesh.getNodalConnectivity().setIJ(0,6,12); + mesh.getNodalConnectivity().setIJ(0,7,10); + mesh.getNodalConnectivity().setIJ(0,8,13); + mesh.getNodalConnectivity().setIJ(0,9,9); + mesh.unPolyze(); + self.assertTrue(not mesh.isEqual(mesh2,1e-12)); + # Test for 2D mesh + mesh=MEDCouplingDataForTest.build2DTargetMesh_1(); + mesh2=MEDCouplingDataForTest.build2DTargetMesh_1(); + eltsV=eltsV[:5]; + mesh.convertToPolyTypes(eltsV); + self.assertTrue(not mesh.isEqual(mesh2,1e-12)); + mesh.unPolyze(); + self.assertTrue(mesh.isEqual(mesh2,1e-12)); + pass + + def testConvertDegeneratedCells1(self): + mesh=MEDCouplingDataForTest.build3DTargetMesh_1(); + conn=[0,1,3,3,9,10,12,12, 0,1,3,4,9,9,9,9, 1,1,1,1,10,12,9,10, 10,11,12,9,1,1,1,1] + mesh.allocateCells(4); + mesh.insertNextCell(NORM_HEXA8,8,conn[0:8]) + mesh.insertNextCell(NORM_HEXA8,8,conn[8:16]) + mesh.insertNextCell(NORM_HEXA8,8,conn[16:24]) + mesh.insertNextCell(NORM_HEXA8,8,conn[24:32]) + mesh.finishInsertingCells(); + mesh.checkCoherency(); + self.assertEqual(4,mesh.getNumberOfCells()); + self.assertEqual(NORM_HEXA8,mesh.getTypeOfCell(0)); + self.assertEqual(NORM_HEXA8,mesh.getTypeOfCell(1)); + self.assertEqual(NORM_HEXA8,mesh.getTypeOfCell(2)); + self.assertEqual(NORM_HEXA8,mesh.getTypeOfCell(3)); + f1=mesh.getMeasureField(True); + mesh.convertDegeneratedCells(); + mesh.checkCoherency(); + f2=mesh.getMeasureField(True); + self.assertEqual(4,mesh.getNumberOfCells()); + self.assertEqual(NORM_PENTA6,mesh.getTypeOfCell(0)); + self.assertEqual(NORM_PYRA5,mesh.getTypeOfCell(1)); + self.assertEqual(NORM_TETRA4,mesh.getTypeOfCell(2)); + self.assertEqual(NORM_PYRA5,mesh.getTypeOfCell(3)); + for i in xrange(4): + self.assertAlmostEqual(f1.getArray().getIJ(0,i),f2.getArray().getIJ(0,i),5); + pass + pass + + def testGetNodeIdsNearPoints1(self): + mesh=MEDCouplingDataForTest.build2DTargetMesh_1(); + coords=mesh.getCoords(); + tmp=DataArrayDouble.New(); + vals=[0.2,0.2,0.1,0.2,0.2,0.2] + tmp.setValues(vals,3,2); + tmp2=DataArrayDouble.aggregate(coords,tmp); + mesh.setCoords(tmp2); + pts=[0.2,0.2,0.1,0.3,-0.3,0.7] + c=mesh.getNodeIdsNearPoint(pts,1e-7); + self.assertEqual([4,9,11],c); + c,cI=mesh.getNodeIdsNearPoints(pts,3,1e-7); + self.assertEqual([0,3,3,4],cI); + self.assertEqual([4,9,11,6],c); + pass + + def testFieldCopyTinyAttrFrom1(self): + f1=MEDCouplingFieldDouble.New(ON_CELLS,ONE_TIME); + f1.setName("f1"); + f1.setTimeTolerance(1.e-5); + f1.setDescription("f1Desc"); + f1.setTime(1.23,4,5); + f2=MEDCouplingFieldDouble.New(ON_CELLS,ONE_TIME); + f2.setName("f2"); + f2.setDescription("f2Desc"); + f2.setTime(6.78,9,10); + f2.setTimeTolerance(4.556e-12); + # + f1.copyTinyAttrFrom(f2); + self.assertAlmostEqual(4.556e-12,f1.getTimeTolerance(),24); + t,dt,it=f1.getTime() + self.assertAlmostEqual(6.78,t,12); + self.assertEqual(9,dt); + self.assertEqual(10,it); + self.assertTrue(f1.getName()=="f1");#name unchanged + self.assertTrue(f1.getDescription()=="f1Desc");#description unchanged + # + f1=MEDCouplingFieldDouble.New(ON_CELLS,NO_TIME); + f1.setName("f1"); + f1.setTimeTolerance(1.e-5); + f1.setDescription("f1Desc"); + f2=MEDCouplingFieldDouble.New(ON_CELLS,NO_TIME); + f2.setName("f2"); + f2.setDescription("f2Desc"); + f2.setTimeTolerance(4.556e-12); + # + f1.copyTinyAttrFrom(f2); + self.assertAlmostEqual(4.556e-12,f1.getTimeTolerance(),24); + self.assertTrue(f1.getName()=="f1");#name unchanged + self.assertTrue(f1.getDescription()=="f1Desc");#description unchanged + # + f1=MEDCouplingFieldDouble.New(ON_CELLS,CONST_ON_TIME_INTERVAL); + f1.setName("f1"); + f1.setTimeTolerance(1.e-5); + f1.setDescription("f1Desc"); + f1.setTime(1.23,4,5); + f1.setEndTime(5.43,2,1); + f2=MEDCouplingFieldDouble.New(ON_CELLS,CONST_ON_TIME_INTERVAL); + f2.setName("f2"); + f2.setDescription("f2Desc"); + f2.setTimeTolerance(4.556e-12); + f2.setTime(6.78,9,10); + f2.setEndTime(10.98,7,6); + # + f1.copyTinyAttrFrom(f2); + self.assertAlmostEqual(4.556e-12,f1.getTimeTolerance(),24); + self.assertTrue(f1.getName()=="f1");#name unchanged + self.assertTrue(f1.getDescription()=="f1Desc");#description unchanged + t,dt,it=f1.getTime() + self.assertAlmostEqual(6.78,t,12); + self.assertEqual(9,dt); + self.assertEqual(10,it); + t,dt,it=f1.getEndTime() + self.assertAlmostEqual(10.98,t,12); + self.assertEqual(7,dt); + self.assertEqual(6,it); + # + f1=MEDCouplingFieldDouble.New(ON_CELLS,LINEAR_TIME); + f1.setName("f1"); + f1.setTimeTolerance(1.e-5); + f1.setDescription("f1Desc"); + f1.setTime(1.23,4,5); + f1.setEndTime(5.43,2,1); + f2=MEDCouplingFieldDouble.New(ON_CELLS,LINEAR_TIME); + f2.setName("f2"); + f2.setDescription("f2Desc"); + f2.setTimeTolerance(4.556e-12); + f2.setTime(6.78,9,10); + f2.setEndTime(10.98,7,6); + # + f1.copyTinyAttrFrom(f2); + self.assertAlmostEqual(4.556e-12,f1.getTimeTolerance(),24); + self.assertTrue(f1.getName()=="f1");#name unchanged + self.assertTrue(f1.getDescription()=="f1Desc");#description unchanged + t,dt,it=f1.getTime() + self.assertAlmostEqual(6.78,t,12); + self.assertEqual(9,dt); + self.assertEqual(10,it); + t,dt,it=f1.getEndTime() + self.assertAlmostEqual(10.98,t,12); + self.assertEqual(7,dt); + self.assertEqual(6,it); + pass + + def testExtrudedMesh5(self): + coo1=[0.,1.,2.,3.5] + a=DataArrayDouble.New(); + a.setValues(coo1,4,1); + b=MEDCouplingCMesh.New(); + b.setCoordsAt(0,a); + c=b.buildUnstructured(); + self.assertEqual(1,c.getSpaceDimension()); + c.changeSpaceDimension(2); + # + d=DataArrayDouble.New(); + d.alloc(13,1); + d.iota(); + e=MEDCouplingCMesh.New(); + e.setCoordsAt(0,d); + f=e.buildUnstructured(); + g=f.getCoords().applyFunc(2,"3.5*IVec+x/6*3.14159265359*JVec"); + h=g.fromPolarToCart(); + f.setCoords(h); + i=c.buildExtrudedMeshFromThis(f,1); + self.assertEqual(52,i.getNumberOfNodes()); + tmp,tmp2,tmp3=i.mergeNodes(1e-9); + self.assertTrue(tmp2); + self.assertEqual(37,tmp3); + i.convertDegeneratedCells(); + i.checkCoherency(); + self.assertEqual(36,i.getNumberOfCells()); + self.assertEqual(37,i.getNumberOfNodes()); + self.assertEqual(12,i.getNumberOfCellsWithType(NORM_TRI3)); + self.assertEqual(24,i.getNumberOfCellsWithType(NORM_QUAD4)); + expected1=[0.25,0.75,2.0625] + j=i.getMeasureField(True); + for ii in xrange(12): + for k in xrange(3): + self.assertAlmostEqual(expected1[k],j.getIJ(0,ii*3+k),10); + pass + pass + expected2=[0.62200846792814113, 0.16666666666681595, 1.4513530918323276, 0.38888888888923495, 2.6293994326053212, 0.7045454545460802, 0.45534180126145435, 0.45534180126150181, 1.0624642029433926, 1.0624642029435025, 1.9248539780597826, 1.9248539780599816, 0.16666666666661334, 0.62200846792815856, 0.38888888888876294, 1.4513530918323678, 0.70454545454522521, 2.629399432605394, -0.16666666666674007, 0.62200846792812436, -0.38888888888906142, 1.4513530918322881, -0.70454545454576778, 2.6293994326052488, -0.45534180126154766, 0.45534180126140844, -1.0624642029436118, 1.0624642029432834, -1.9248539780601803, 1.9248539780595841, -0.62200846792817499, 0.1666666666665495, -1.451353091832408, 0.388888888888613, -2.6293994326054668, 0.70454545454495332, -0.62200846792810593, -0.16666666666680507, -1.451353091832247, -0.38888888888921297, -2.6293994326051746, -0.70454545454604123, -0.45534180126135926, -0.45534180126159562, -1.0624642029431723, -1.0624642029437235, -1.9248539780593836, -1.9248539780603811, -0.1666666666664828, -0.62200846792819242, -0.38888888888846079, -1.4513530918324489, -0.70454545454467987, -2.6293994326055397, 0.16666666666687083, -0.62200846792808862, 0.38888888888936374, -1.4513530918322073, 0.70454545454631357, -2.6293994326051022, 0.45534180126164348, -0.45534180126131207, 1.0624642029438327, -1.0624642029430627, 1.9248539780605791, -1.9248539780591853, 0.62200846792821063, -0.16666666666641802, 1.4513530918324888, -0.38888888888831086, 2.6293994326056125, -0.70454545454440853] + m=i.getBarycenterAndOwner(); + for i in xrange(72): + self.assertAlmostEqual(expected2[i],m.getIJ(0,i),10); + pass + # + pass + + def testExtrudedMesh6(self): + coo1=[0.,1.,2.,3.5] + a=DataArrayDouble.New(); + a.setValues(coo1,4,1); + b=MEDCouplingCMesh.New(); + b.setCoordsAt(0,a); + c=b.buildUnstructured(); + self.assertEqual(1,c.getSpaceDimension()); + c.changeSpaceDimension(2); + # + d=DataArrayDouble.New(); + d.alloc(5,1); + d.iota(); + e=MEDCouplingCMesh.New(); + e.setCoordsAt(0,d); + f=e.buildUnstructured(); + d2=f.getCoords().applyFunc("x*x/2"); + f.setCoords(d2); + f.changeSpaceDimension(2); + # + center=[0.,0.] + f.rotate(center,[],pi/3); + g=c.buildExtrudedMeshFromThis(f,0); + g.checkCoherency(); + expected1=[ 0.4330127018922193, 0.4330127018922193, 0.649519052838329, 1.2990381056766578, 1.299038105676658, 1.948557158514987, 2.1650635094610955, 2.1650635094610964, 3.2475952641916446, 3.031088913245533, 3.0310889132455352, 4.546633369868303 ] + f1=g.getMeasureField(True); + for i in xrange(12): + self.assertAlmostEqual(expected1[i],f1.getIJ(0,i),12); + pass + expected2=[0.625, 0.21650635094610962, 1.625, 0.21650635094610959, 2.8750000000000004, 0.21650635094610965, 1.1250000000000002, 1.0825317547305482, 2.125, 1.0825317547305482, 3.3750000000000004, 1.0825317547305484, 2.125, 2.8145825622994254, 3.125, 2.8145825622994254, 4.375, 2.8145825622994254, 3.6250000000000009, 5.4126587736527414, 4.625, 5.4126587736527414, 5.875, 5.4126587736527414] + f2=g.getBarycenterAndOwner(); + for i in xrange(24): + self.assertAlmostEqual(expected2[i],f2.getIJ(0,i),12); + pass + pass + + def testExtrudedMesh7(self): + coo1=[0.,1.,2.,3.5] + a=DataArrayDouble.New(); + a.setValues(coo1,4,1); + b=MEDCouplingCMesh.New(); + b.setCoordsAt(0,a); + c=b.buildUnstructured(); + self.assertEqual(1,c.getSpaceDimension()); + c.changeSpaceDimension(2); + # + d=DataArrayDouble.New(); + d.alloc(13,1); + d.iota(); + e=MEDCouplingCMesh.New(); + e.setCoordsAt(0,d); + f=e.buildUnstructured(); + g=f.getCoords().applyFunc(2,"3.5*IVec+x/6*3.14159265359*JVec"); + h=g.fromPolarToCart(); + f.setCoords(h); + i=c.buildExtrudedMeshFromThis(f,1); + self.assertEqual(52,i.getNumberOfNodes()); + tmp,tmp2,tmp3=i.mergeNodes(1e-9); + self.assertTrue(tmp2); + self.assertEqual(37,tmp3); + i.convertDegeneratedCells(); + vec1=[10.,0.,0.] + i.translate(vec1); + g2=h.applyFunc(3,"13.5/3.5*x*IVec+0*JVec+13.5/3.5*y*KVec"); + f.setCoords(g2); + i.changeSpaceDimension(3); + i3=i.buildExtrudedMeshFromThis(f,1); + f2=i3.getMeasureField(True); + tmp,tmp2,tmp3=i.mergeNodes(1e-9); + self.assertTrue(tmp2); + self.assertEqual(444,tmp3); + expected1=[1.327751058489274, 4.2942574094314701, 13.024068164857139, 1.3069177251569044, 4.1484240761012954, 12.297505664866796, 1.270833333332571, 3.8958333333309674, 11.039062499993179, 1.2291666666659207, 3.6041666666644425, 9.585937499993932, 1.1930822748415895, 3.3515759238941376, 8.3274943351204556, 1.1722489415082769, 3.2057425905609289, 7.6009318351210622, 1.1722489415082862, 3.2057425905609884, 7.6009318351213713, 1.1930822748416161, 3.3515759238943001, 8.3274943351212727, 1.2291666666659564, 3.6041666666646734, 9.5859374999950777, 1.2708333333326081, 3.8958333333311868, 11.039062499994293, 1.3069177251569224, 4.1484240761014384, 12.297505664867627, 1.3277510584902354, 4.2942574094346071, 13.024068164866796] + for ii in xrange(12): + for jj in xrange(36): + self.assertAlmostEqual(expected1[jj],f2.getIJ(0,ii*36+jj),9); + pass + # + pass + + def testSimplexize1(self): + m=MEDCouplingDataForTest.build3DSurfTargetMesh_1(); + m.convertToPolyTypes([3]); + da=m.simplexize(0); + self.assertEqual(7,da.getNumberOfTuples()); + self.assertEqual(1,da.getNumberOfComponents()); + expected2=[0,0,1,2,3,4,4] + for i in xrange(7): + self.assertEqual(expected2[i],da.getIJ(i,0)); + pass + m.checkCoherency(); + self.assertEqual(7,m.getNumberOfCells()); + self.assertEqual(NORM_TRI3,m.getTypeOfCell(0)); + self.assertEqual(NORM_TRI3,m.getTypeOfCell(1)); + self.assertEqual(NORM_TRI3,m.getTypeOfCell(2)); + self.assertEqual(NORM_TRI3,m.getTypeOfCell(3)); + self.assertEqual(NORM_POLYGON,m.getTypeOfCell(4)); + self.assertEqual(NORM_TRI3,m.getTypeOfCell(5)); + self.assertEqual(NORM_TRI3,m.getTypeOfCell(6)); + expected1=[0.125,0.125,0.125,0.125,0.25,0.125,0.125] + f=m.getMeasureField(False); + for i in xrange(7): + self.assertAlmostEqual(expected1[i]*sqrt(2.),f.getIJ(i,0),10); + pass + types=m.getAllTypes(); + self.assertEqual([NORM_TRI3,NORM_POLYGON],types); + # + m=MEDCouplingDataForTest.build3DSurfTargetMesh_1(); + m.convertToPolyTypes([3]); + da=m.simplexize(1); + self.assertEqual(7,da.getNumberOfTuples()); + self.assertEqual(1,da.getNumberOfComponents()); + for i in xrange(7): + self.assertEqual(expected2[i],da.getIJ(i,0)); + pass + m.checkCoherency(); + types=m.getAllTypes(); + self.assertEqual([NORM_TRI3,NORM_POLYGON],types); + self.assertEqual(7,m.getNumberOfCells()); + self.assertEqual(NORM_TRI3,m.getTypeOfCell(0)); + self.assertEqual(NORM_TRI3,m.getTypeOfCell(1)); + self.assertEqual(NORM_TRI3,m.getTypeOfCell(2)); + self.assertEqual(NORM_TRI3,m.getTypeOfCell(3)); + self.assertEqual(NORM_POLYGON,m.getTypeOfCell(4)); + self.assertEqual(NORM_TRI3,m.getTypeOfCell(5)); + self.assertEqual(NORM_TRI3,m.getTypeOfCell(6)); + f=m.getMeasureField(False); + for i in xrange(7): + self.assertAlmostEqual(expected1[i]*sqrt(2.),f.getIJ(i,0),10); + pass + pass + + def testSimplexize2(self): + m=MEDCouplingDataForTest.build3DSurfTargetMesh_1(); + m.convertToPolyTypes([3]); + f1=MEDCouplingFieldDouble.New(ON_CELLS,ONE_TIME); + f1.setMesh(m); + arr=DataArrayDouble.New(); + arr1=[10.,110.,20.,120.,30.,130.,40.,140.,50.,150.] + arr.setValues(arr1,5,2); + f1.setArray(arr); + # + f1.checkCoherency(); + self.assertTrue(f1.simplexize(0)); + f1.checkCoherency(); + expected1=[10.,110.,10.,110.,20.,120.,30.,130.,40.,140.,50.,150.,50.,150.] + for i in xrange(14): + self.assertAlmostEqual(expected1[i],f1.getIJ(0,i),10); + pass + self.assertTrue(not f1.simplexize(0)); + for i in xrange(14): + self.assertAlmostEqual(expected1[i],f1.getIJ(0,i),10); + pass + # + pass def setUp(self): pass diff --git a/src/MEDCoupling_Swig/libMEDCoupling_Swig.i b/src/MEDCoupling_Swig/libMEDCoupling_Swig.i index d6b35f447..65112f642 100644 --- a/src/MEDCoupling_Swig/libMEDCoupling_Swig.i +++ b/src/MEDCoupling_Swig/libMEDCoupling_Swig.i @@ -59,11 +59,6 @@ using namespace INTERP_KERNEL; %feature("autodoc", "1"); %feature("docstring"); -%newobject ParaMEDMEM::DataArrayDouble::New; -%newobject ParaMEDMEM::DataArrayInt::New; -%newobject ParaMEDMEM::DataArrayDouble::convertToIntArr; -%newobject ParaMEDMEM::DataArrayInt::convertToDblArr; -%newobject ParaMEDMEM::MEDCouplingUMesh::New; %newobject ParaMEDMEM::MEDCouplingField::buildMeasureField; %newobject ParaMEDMEM::MEDCouplingFieldDouble::New; %newobject ParaMEDMEM::MEDCouplingFieldDouble::mergeFields; @@ -91,9 +86,11 @@ using namespace INTERP_KERNEL; %newobject ParaMEDMEM::MEDCouplingFieldDouble::operator-; %newobject ParaMEDMEM::MEDCouplingFieldDouble::operator*; %newobject ParaMEDMEM::MEDCouplingFieldDouble::operator/; -%newobject ParaMEDMEM::MEDCouplingUMesh::clone; -%newobject ParaMEDMEM::DataArrayDouble::deepCopy; -%newobject ParaMEDMEM::DataArrayDouble::performCpy; +%newobject ParaMEDMEM::MEDCouplingFieldDouble::clone; +%newobject ParaMEDMEM::MEDCouplingFieldDouble::cloneWithMesh; +%newobject ParaMEDMEM::MEDCouplingFieldDouble::buildNewTimeReprFromThis; +%newobject ParaMEDMEM::DataArrayInt::New; +%newobject ParaMEDMEM::DataArrayInt::convertToDblArr; %newobject ParaMEDMEM::DataArrayInt::deepCopy; %newobject ParaMEDMEM::DataArrayInt::performCpy; %newobject ParaMEDMEM::DataArrayInt::substr; @@ -105,6 +102,15 @@ using namespace INTERP_KERNEL; %newobject ParaMEDMEM::DataArrayInt::renumberAndReduce; %newobject ParaMEDMEM::DataArrayInt::invertArrayO2N2N2O; %newobject ParaMEDMEM::DataArrayInt::invertArrayN2O2O2N; +%newobject ParaMEDMEM::DataArrayInt::getIdsEqual; +%newobject ParaMEDMEM::DataArrayInt::getIdsEqualList; +%newobject ParaMEDMEM::DataArrayInt::aggregate; +%newobject ParaMEDMEM::DataArrayInt::fromNoInterlace; +%newobject ParaMEDMEM::DataArrayInt::toNoInterlace; +%newobject ParaMEDMEM::DataArrayDouble::New; +%newobject ParaMEDMEM::DataArrayDouble::convertToIntArr; +%newobject ParaMEDMEM::DataArrayDouble::deepCopy; +%newobject ParaMEDMEM::DataArrayDouble::performCpy; %newobject ParaMEDMEM::DataArrayDouble::aggregate; %newobject ParaMEDMEM::DataArrayDouble::dot; %newobject ParaMEDMEM::DataArrayDouble::crossProduct; @@ -130,21 +136,25 @@ using namespace INTERP_KERNEL; %newobject ParaMEDMEM::DataArrayDouble::renumber; %newobject ParaMEDMEM::DataArrayDouble::renumberR; %newobject ParaMEDMEM::DataArrayDouble::renumberAndReduce; -%newobject ParaMEDMEM::MEDCouplingFieldDouble::clone; -%newobject ParaMEDMEM::MEDCouplingFieldDouble::cloneWithMesh; -%newobject ParaMEDMEM::MEDCouplingFieldDouble::buildNewTimeReprFromThis; +%newobject ParaMEDMEM::DataArrayDouble::fromNoInterlace; +%newobject ParaMEDMEM::DataArrayDouble::toNoInterlace; +%newobject ParaMEDMEM::DataArrayDouble::fromPolarToCart; +%newobject ParaMEDMEM::DataArrayDouble::fromCylToCart; +%newobject ParaMEDMEM::DataArrayDouble::fromSpherToCart; %newobject ParaMEDMEM::MEDCouplingMesh::getCoordinatesAndOwner; %newobject ParaMEDMEM::MEDCouplingMesh::getBarycenterAndOwner; %newobject ParaMEDMEM::MEDCouplingMesh::buildOrthogonalField; %newobject ParaMEDMEM::MEDCouplingMesh::getCellIdsFullyIncludedInNodeIds; -%newobject ParaMEDMEM::MEDCouplingMesh::buildPart; %newobject ParaMEDMEM::MEDCouplingMesh::mergeMyselfWith; %newobject ParaMEDMEM::MEDCouplingMesh::fillFromAnalytic; +%newobject ParaMEDMEM::MEDCouplingMesh::getMeasureField; +%newobject ParaMEDMEM::MEDCouplingMesh::simplexize; %newobject ParaMEDMEM::MEDCouplingPointSet::zipCoordsTraducer; +%newobject ParaMEDMEM::MEDCouplingPointSet::buildBoundaryMesh; +%newobject ParaMEDMEM::MEDCouplingUMesh::New; +%newobject ParaMEDMEM::MEDCouplingUMesh::clone; %newobject ParaMEDMEM::MEDCouplingUMesh::zipConnectivityTraducer; %newobject ParaMEDMEM::MEDCouplingUMesh::buildDescendingConnectivity; -%newobject ParaMEDMEM::MEDCouplingPointSet::buildBoundaryMesh; -%newobject ParaMEDMEM::MEDCouplingMesh::getMeasureField; %newobject ParaMEDMEM::MEDCouplingUMesh::buildExtrudedMeshFromThis; %newobject ParaMEDMEM::MEDCouplingUMesh::mergeUMeshes; %newobject ParaMEDMEM::MEDCouplingUMesh::buildNewNumberingFromCommNodesFrmt; @@ -227,20 +237,18 @@ namespace ParaMEDMEM virtual int getMeshDimension() const throw(INTERP_KERNEL::Exception) = 0; virtual DataArrayDouble *getCoordinatesAndOwner() const = 0; virtual DataArrayDouble *getBarycenterAndOwner() const = 0; + virtual int getNumberOfCellsWithType(INTERP_KERNEL::NormalizedCellType type) const = 0; virtual INTERP_KERNEL::NormalizedCellType getTypeOfCell(int cellId) const = 0; virtual std::string simpleRepr() const = 0; virtual std::string advancedRepr() const = 0; // tools - virtual void getBoundingBox(double *bbox) const = 0; virtual MEDCouplingFieldDouble *getMeasureField(bool isAbs) const = 0; virtual MEDCouplingFieldDouble *getMeasureFieldOnNode(bool isAbs) const = 0; virtual MEDCouplingFieldDouble *fillFromAnalytic(TypeOfField t, int nbOfComp, const char *func) const throw(INTERP_KERNEL::Exception); 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 MEDCouplingMesh *buildPart(const int *start, const int *end) const = 0; virtual MEDCouplingMesh *mergeMyselfWith(const MEDCouplingMesh *other) const throw(INTERP_KERNEL::Exception) = 0; virtual bool areCompatibleForMerge(const MEDCouplingMesh *other) const; + virtual DataArrayInt *simplexize(int policy) throw(INTERP_KERNEL::Exception); static MEDCouplingMesh *mergeMeshes(const MEDCouplingMesh *mesh1, const MEDCouplingMesh *mesh2); %extend { @@ -304,6 +312,61 @@ namespace ParaMEDMEM self->scale(p,factor); delete [] p; } + + PyObject *getBoundingBox() const throw(INTERP_KERNEL::Exception) + { + int spaceDim=self->getSpaceDimension(); + double *tmp=new double[2*spaceDim]; + self->getBoundingBox(tmp); + PyObject *ret=convertDblArrToPyListOfTuple(tmp,2,spaceDim); + delete [] tmp; + return ret; + } + + PyObject *buildPart(PyObject *li) const + { + int size; + int *tmp=convertPyToNewIntArr2(li,&size); + MEDCouplingMesh *ret=self->buildPart(tmp,tmp+size); + delete [] tmp; + return convertMesh(ret, SWIG_POINTER_OWN | 0 ); + } + + PyObject *buildPartAndReduceNodes(PyObject *li) const + { + int size; + int *tmp=convertPyToNewIntArr2(li,&size); + DataArrayInt *arr=0; + MEDCouplingMesh *ret=self->buildPartAndReduceNodes(tmp,tmp+size,arr); + PyObject *res = PyList_New(2); + PyObject *obj0=convertMesh(ret, SWIG_POINTER_OWN | 0 ); + PyObject *obj1=SWIG_NewPointerObj(SWIG_as_voidptr(arr),SWIGTYPE_p_ParaMEDMEM__DataArrayInt, SWIG_POINTER_OWN | 0 ); + PyList_SetItem(res,0,obj0); + PyList_SetItem(res,1,obj1); + return res; + } + + void translate(PyObject *vector) + { + int sz; + double *v=convertPyToNewDblArr2(vector,&sz); + self->translate(v); + delete [] v; + } + + void rotate(PyObject *center, PyObject *vector, double alpha) + { + int sz; + double *c=convertPyToNewDblArr2(center,&sz); + if(!c) + return ; + double *v=convertPyToNewDblArr2(vector,&sz); + if(!v) + { delete [] c; return ; } + self->rotate(c,v,alpha); + delete [] c; + delete [] v; + } } }; } @@ -323,20 +386,14 @@ namespace ParaMEDMEM void setCoords(DataArrayDouble *coords); DataArrayDouble *getCoordinatesAndOwner() const; bool areCoordsEqual(const MEDCouplingPointSet& other, double prec) const; - void getBoundingBox(double *bbox) const; void zipCoords(); double getCaracteristicDimension() const; - void translate(const double *vector); void changeSpaceDimension(int newSpaceDim, double dftVal=0.) throw(INTERP_KERNEL::Exception); void tryToShareSameCoords(const MEDCouplingPointSet& other, double epsilon) throw(INTERP_KERNEL::Exception); virtual void tryToShareSameCoordsPermute(const MEDCouplingPointSet& other, double epsilon) throw(INTERP_KERNEL::Exception) = 0; static DataArrayDouble *mergeNodesArray(const MEDCouplingPointSet *m1, const MEDCouplingPointSet *m2); static MEDCouplingPointSet *buildInstanceFromMeshType(MEDCouplingMeshType type); - virtual MEDCouplingPointSet *buildPartOfMySelf(const int *start, const int *end, bool keepCoords) const = 0; - virtual MEDCouplingPointSet *buildPartOfMySelfNode(const int *start, const int *end, bool fullyIn) const = 0; - virtual MEDCouplingPointSet *buildFacePartOfMySelfNode(const int *start, const int *end, bool fullyIn) const = 0; virtual MEDCouplingPointSet *buildBoundaryMesh(bool keepCoords) const = 0; - virtual void renumberNodes(const int *newNodeNumbers, int newNbOfNodes); virtual bool isEmptyMesh(const std::vector& tinyInfo) const = 0; //! size of returned tinyInfo must be always the same. void getTinySerializationInformation(std::vector& tinyInfo, std::vector& littleStrings) const; @@ -344,7 +401,6 @@ namespace ParaMEDMEM void serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const; void unserialization(const std::vector& tinyInfo, const DataArrayInt *a1, DataArrayDouble *a2, const std::vector& littleStrings); - virtual void giveElemsInBoundingBox(const double *bbox, double eps, std::vector& elems) = 0; virtual void giveElemsInBoundingBox(const INTERP_KERNEL::DirectedBoundingBox& bbox, double eps, std::vector& elems) = 0; virtual DataArrayInt *zipCoordsTraducer() = 0; %extend @@ -410,26 +466,6 @@ namespace ParaMEDMEM self->findBoundaryNodes(nodes); return convertIntArrToPyList2(nodes); } - void rotate(PyObject *center, PyObject *vector, double alpha) - { - int sz; - double *c=convertPyToNewDblArr2(center,&sz); - if(!c) - return ; - double *v=convertPyToNewDblArr2(vector,&sz); - if(!v) - { delete [] c; return ; } - self->rotate(c,v,alpha); - delete [] c; - delete [] v; - } - void translate(PyObject *vector) - { - int sz; - double *v=convertPyToNewDblArr2(vector,&sz); - self->translate(v); - delete [] v; - } void renumberNodes(PyObject *li, int newNbOfNodes) { int size; @@ -448,6 +484,65 @@ namespace ParaMEDMEM delete [] p; return convertIntArrToPyList2(nodes); } + PyObject *getNodeIdsNearPoint(PyObject *pt, double eps) const throw(INTERP_KERNEL::Exception) + { + int size; + double *pos=convertPyToNewDblArr2(pt,&size); + if(sizegetSpaceDimension()) + { + delete [] pos; + throw INTERP_KERNEL::Exception("getNodeIdsNearPoint : to tiny array ! must be at least of size SpaceDim !"); + } + std::vector tmp; + try + { + tmp=self->getNodeIdsNearPoint(pos,eps); + } + catch(INTERP_KERNEL::Exception& e) + { + delete [] pos; + throw e; + } + delete [] pos; + return convertIntArrToPyList2(tmp); + } + + PyObject *getNodeIdsNearPoints(PyObject *pt, int nbOfNodes, double eps) const throw(INTERP_KERNEL::Exception) + { + std::vector c,cI; + int size; + double *pos=convertPyToNewDblArr2(pt,&size); + if(sizegetSpaceDimension()*nbOfNodes) + { + delete [] pos; + throw INTERP_KERNEL::Exception("getNodeIdsNearPoints : to tiny array ! must be at least of size SpaceDim*nbOfNodes !"); + } + try + { + self->getNodeIdsNearPoints(pos,nbOfNodes,eps,c,cI); + } + catch(INTERP_KERNEL::Exception& e) + { + delete [] pos; + throw e; + } + delete [] pos; + PyObject *ret=PyTuple_New(2); + PyTuple_SetItem(ret,0,convertIntArrToPyList2(c)); + PyTuple_SetItem(ret,1,convertIntArrToPyList2(cI)); + return ret; + } + + PyObject *giveElemsInBoundingBox(PyObject *bbox, double eps) + { + std::vector elems; + int size; + double *tmp=convertPyToNewDblArr2(bbox,&size); + self->giveElemsInBoundingBox(tmp,eps,elems); + delete [] tmp; + return convertIntArrToPyList2(elems); + } + static void rotate2DAlg(PyObject *center, double angle, int nbNodes, PyObject *coords) { int sz; @@ -503,7 +598,10 @@ namespace ParaMEDMEM void orientCorrectlyPolyhedrons() throw(INTERP_KERNEL::Exception); bool isPresenceOfQuadratic() const; MEDCouplingFieldDouble *buildDirectionVectorField() const; + bool isContiguous1D() const throw(INTERP_KERNEL::Exception); void convertQuadraticCellsToLinear() throw(INTERP_KERNEL::Exception); + void convertDegeneratedCells() throw(INTERP_KERNEL::Exception); + bool areOnlySimplexCells() const throw(INTERP_KERNEL::Exception); MEDCouplingFieldDouble *getEdgeRatioField() const throw(INTERP_KERNEL::Exception); MEDCouplingFieldDouble *getAspectRatioField() const throw(INTERP_KERNEL::Exception); MEDCouplingFieldDouble *getWarpField() const throw(INTERP_KERNEL::Exception); @@ -701,6 +799,7 @@ namespace ParaMEDMEM } } void convertToPolyTypes(const std::vector& cellIdsToConvert); + void unPolyze(); MEDCouplingUMesh *buildExtrudedMeshFromThis(const MEDCouplingUMesh *mesh1D, int policy); static MEDCouplingUMesh *mergeUMeshes(const MEDCouplingUMesh *mesh1, const MEDCouplingUMesh *mesh2) throw(INTERP_KERNEL::Exception); }; @@ -915,6 +1014,24 @@ namespace ParaMEDMEM PyTuple_SetItem(ret,1,SWIG_NewPointerObj(SWIG_as_voidptr(tmp),SWIGTYPE_p_ParaMEDMEM__DataArrayInt, SWIG_POINTER_OWN | 0 )); return ret; } + + PyObject *accumulate() const throw(INTERP_KERNEL::Exception) + { + int sz=self->getNumberOfComponents(); + double *tmp=new double[sz]; + try + { + self->accumulate(tmp); + } + catch(INTERP_KERNEL::Exception& e) + { + delete [] tmp; + throw e; + } + PyObject *ret=convertDblArrToPyList(tmp,sz); + delete [] tmp; + return ret; + } DataArrayDouble *keepSelectedComponents(PyObject *li) const throw(INTERP_KERNEL::Exception) { @@ -1141,6 +1258,7 @@ namespace ParaMEDMEM public: static MEDCouplingFieldDouble *New(TypeOfField type, TypeOfTimeDiscretization td=NO_TIME); void copyTinyStringsFrom(const MEDCouplingFieldDouble *other) throw(INTERP_KERNEL::Exception); + void copyTinyAttrFrom(const MEDCouplingFieldDouble *other) throw(INTERP_KERNEL::Exception); std::string simpleRepr() const; std::string advancedRepr() const; MEDCouplingFieldDouble *clone(bool recDeepCpy) const; @@ -1163,12 +1281,15 @@ namespace ParaMEDMEM int getNumberOfValues() const throw(INTERP_KERNEL::Exception); NatureOfField getNature() const { return _nature; } void setNature(NatureOfField nat) throw(INTERP_KERNEL::Exception); + void setTimeTolerance(double val); + double getTimeTolerance() const; void updateTime(); void changeUnderlyingMesh(const MEDCouplingMesh *other, int levOfCheck, double prec) throw(INTERP_KERNEL::Exception); void substractInPlaceDM(const MEDCouplingFieldDouble *f, int levOfCheck, double prec) throw(INTERP_KERNEL::Exception); bool mergeNodes(double eps) throw(INTERP_KERNEL::Exception); bool zipCoords() throw(INTERP_KERNEL::Exception); bool zipConnectivity(int compType) throw(INTERP_KERNEL::Exception); + bool simplexize(int policy) throw(INTERP_KERNEL::Exception); MEDCouplingFieldDouble *doublyContractedProduct() const throw(INTERP_KERNEL::Exception); MEDCouplingFieldDouble *determinant() const throw(INTERP_KERNEL::Exception); MEDCouplingFieldDouble *eigenValues() const throw(INTERP_KERNEL::Exception); diff --git a/src/ParaMEDMEM/InterpolationMatrix.cxx b/src/ParaMEDMEM/InterpolationMatrix.cxx index c7db0371c..1022280b1 100644 --- a/src/ParaMEDMEM/InterpolationMatrix.cxx +++ b/src/ParaMEDMEM/InterpolationMatrix.cxx @@ -526,7 +526,7 @@ namespace ParaMEDMEM for(vector >::const_iterator iter3=(*iter).begin();iter3!=(*iter).end();iter3++) res[(*iter3).first]+=(*iter3).second; set procsSet; - int id; + int id=-1; const vector >& mapping=_mapping.getSendingIds(); for(vector >::const_iterator iter2=mapping.begin();iter2!=mapping.end();iter2++) {