typedef enum
{
+ NORM_POINT0 = 0,
NORM_SEG2 = 1,
NORM_SEG3 = 2,
NORM_TRI3 = 3,
void CellModel::buildUniqueInstance()
{
+ _map_of_unique_instance.insert(make_pair(NORM_POINT0,CellModel(NORM_POINT0)));
_map_of_unique_instance.insert(make_pair(NORM_SEG2,CellModel(NORM_SEG2)));
_map_of_unique_instance.insert(make_pair(NORM_SEG3,CellModel(NORM_SEG3)));
_map_of_unique_instance.insert(make_pair(NORM_TRI3,CellModel(NORM_TRI3)));
_dyn=false;
switch(type)
{
+ case NORM_POINT0:
+ {
+ _nb_of_pts=0; _nb_of_sons=0; _dim=0;
+ }
+ break;
case NORM_SEG2:
{
_nb_of_pts=2; _nb_of_sons=0; _dim=1;
#include "PlanarIntersectorP0P0.hxx"
#include "PlanarIntersectorP0P1.hxx"
#include "PlanarIntersectorP1P0.hxx"
-#include "InterpolationUtils.hxx"
+#include "PlanarIntersectorP1P1.hxx"
namespace INTERP_KERNEL
{
bool doRotate, int orientation, int printLevel);
double intersectGeometry(ConnType icellT, ConnType icellS, ConnType nbNodesT, ConnType nbNodesS);
double intersectGeometryWithQuadrangle(const double *quadrangle, const std::vector<double>& sourceCoords, bool isSourceQuad);
+ double intersectGeometryGeneral(const std::vector<double>& targetCoords, const std::vector<double>& sourceCoords);
private :
double _epsilon;
};
#include "PlanarIntersectorP0P0.txx"
#include "PlanarIntersectorP0P1.txx"
#include "PlanarIntersectorP1P0.txx"
+#include "PlanarIntersectorP1P1.txx"
#include "PolygonAlgorithms.txx"
return result;
}
+
+ template<class MyMeshType, class MyMatrix, template <class MeshType, class TheMatrix, class ThisIntersector> class InterpType>
+ double ConvexIntersector<MyMeshType,MyMatrix,InterpType>::intersectGeometryGeneral(const std::vector<double>& targetCoords, const std::vector<double>& sourceCoords)
+ {
+ double result = 0;
+ int nbOfNodesS=sourceCoords.size()/SPACEDIM;
+ int nbOfNodesT=targetCoords.size()/SPACEDIM;
+ /*** Compute the intersection area ***/
+ INTERP_KERNEL::PolygonAlgorithms<SPACEDIM> P(_epsilon, PlanarIntersector<MyMeshType,MyMatrix>::_precision);
+ std::deque<double> inter = P.intersectConvexPolygons(&targetCoords[0], &sourceCoords[0],
+ nbOfNodesT, nbOfNodesS);
+ double area[SPACEDIM];
+ int nb_inter =((int)inter.size())/SPACEDIM;
+ for(int i = 1; i<nb_inter-1; i++)
+ {
+ INTERP_KERNEL::crossprod<SPACEDIM>(&inter[0],&inter[SPACEDIM*i],&inter[SPACEDIM*(i+1)],area);
+ result +=0.5*norm<SPACEDIM>(area);
+ }
+ return result;
+ }
}
#endif
#include "PlanarIntersectorP0P0.hxx"
#include "PlanarIntersectorP0P1.hxx"
#include "PlanarIntersectorP1P0.hxx"
+#include "PlanarIntersectorP1P1.hxx"
namespace INTERP_KERNEL
{
double dimCaracteristic, double medianPlane, double precision, int orientation);
double intersectGeometry(ConnType icellT, ConnType icellS, ConnType nbNodesT, ConnType nbNodesS);
double intersectGeometryWithQuadrangle(const double *quadrangle, const std::vector<double>& sourceCoords, bool isSourceQuad);
+ double intersectGeometryGeneral(const std::vector<double>& targetCoords, const std::vector<double>& sourceCoords);
private:
QuadraticPolygon *buildPolygonFrom(const std::vector<double>& coords, NormalizedCellType type);
QuadraticPolygon *buildPolygonAFrom(ConnType cell, int nbOfPoints, NormalizedCellType type);
#include "PlanarIntersectorP0P0.txx"
#include "PlanarIntersectorP0P1.txx"
#include "PlanarIntersectorP1P0.txx"
+#include "PlanarIntersectorP1P1.txx"
#include "CellModel.hxx"
#include "QuadraticPolygon.hxx"
return ret;
}
+ template<class MyMeshType, class MyMatrix, template <class MeshType, class TheMatrix, class ThisIntersector> class InterpType>
+ double Geometric2DIntersector<MyMeshType,MyMatrix,InterpType>::intersectGeometryGeneral(const std::vector<double>& targetCoords, const std::vector<double>& sourceCoords)
+ {
+ int nbOfTargetNodes=targetCoords.size()/SPACEDIM;
+ std::vector<Node *> nodes(nbOfTargetNodes);
+ for(int i=0;i<nbOfTargetNodes;i++)
+ nodes[i]=new Node(targetCoords[i*SPACEDIM],targetCoords[i*SPACEDIM+1]);
+ int nbOfSourceNodes=sourceCoords.size()/SPACEDIM;
+ std::vector<Node *> nodes2(nbOfSourceNodes);
+ for(int i=0;i<nbOfSourceNodes;i++)
+ nodes2[i]=new Node(sourceCoords[i*SPACEDIM],sourceCoords[i*SPACEDIM+1]);
+ QuadraticPolygon *p1=QuadraticPolygon::buildLinearPolygon(nodes);
+ QuadraticPolygon *p2=QuadraticPolygon::buildLinearPolygon(nodes2);
+ double ret=p1->intersectWith(*p2);
+ delete p1; delete p2;
+ return ret;
+ }
+
template<class MyMeshType, class MyMatrix, template <class MeshType, class TheMatrix, class ThisIntersector> class InterpType>
QuadraticPolygon *Geometric2DIntersector<MyMeshType,MyMatrix,InterpType>::buildPolygonFrom(const std::vector<double>& coords, NormalizedCellType type)
{
}
Interpolation3DSurf::Interpolation3DSurf(const InterpolationOptions& io):InterpolationPlanar<Interpolation3DSurf>(io)
+ ,_median_plane(io.getMedianPlane())
+ ,_surf_3D_adjustment_eps(io.getBoundingBoxAdjustment())
{
}
double getPrecision() const { return InterpolationOptions::_precision; }
void setPrecision(double p) { InterpolationOptions::_precision=p; }
- double getMedianPlane() { return InterpolationOptions::_median_plane; }
+ double getMedianPlane() const { return InterpolationOptions::_median_plane; }
void setMedianPlane(double mp) { InterpolationOptions::_median_plane=mp; }
- bool getDoRotate() { return InterpolationOptions::_do_rotate; }
+ bool getDoRotate() const { return InterpolationOptions::_do_rotate; }
void setDoRotate( bool dr) { InterpolationOptions::_do_rotate = dr; }
- double getBoundingBoxAdjustment() { return InterpolationOptions::_bounding_box_adjustment; }
+ double getBoundingBoxAdjustment() const { return InterpolationOptions::_bounding_box_adjustment; }
void setBoundingBoxAdjustment(double bba) { InterpolationOptions::_bounding_box_adjustment=bba; }
- int getOrientation() { return InterpolationOptions::_orientation; }
+ int getOrientation() const { return InterpolationOptions::_orientation; }
void setOrientation(int o) { InterpolationOptions::_orientation=o; }
- SplittingPolicy getSplittingPolicy() { return _splitting_policy; }
+ SplittingPolicy getSplittingPolicy() const { return _splitting_policy; }
void setSplittingPolicy(SplittingPolicy sp) { _splitting_policy=sp; }
void init()
{
_print_level=0;
_intersection_type=Triangulation;
- _precision=1e-12;;
+ _precision=1e-12;
_median_plane=0.5;
_do_rotate=true;
- _bounding_box_adjustment=0.1;
+ _bounding_box_adjustment=1e-4;
_orientation=0;
_splitting_policy=GENERAL_48;
}
break;
}
}
+ else if(meth=="P1P1")
+ {
+ switch (InterpolationOptions::getIntersectionType())
+ {
+ case Triangulation:
+ intersector=new TriangulationIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P1>(myMeshT,myMeshS,_dim_caracteristic,
+ InterpolationOptions::getPrecision(),
+ InterpolationOptions::getMedianPlane(),
+ InterpolationOptions::getOrientation(),
+ InterpolationOptions::getPrintLevel());
+ break;
+ case Convex:
+ intersector=new ConvexIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P1>(myMeshT,myMeshS,_dim_caracteristic,
+ InterpolationOptions::getPrecision(),
+ InterpolationOptions::getDoRotate(),
+ InterpolationOptions::getMedianPlane(),
+ InterpolationOptions::getOrientation(),
+ InterpolationOptions::getPrintLevel());
+ break;
+ case Geometric2D:
+ intersector=new Geometric2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P1>(myMeshT, myMeshS, _dim_caracteristic,
+ InterpolationOptions::getMedianPlane(),
+ InterpolationOptions::getPrecision(),
+ InterpolationOptions::getOrientation());
+ break;
+ }
+ }
else
- throw INTERP_KERNEL::Exception("Invalid method specified ! Must be in : \"P0P0\" \"P0P1\" or \"P1P0\"");
+ throw INTERP_KERNEL::Exception("Invalid method specified ! Must be in : \"P0P0\" \"P0P1\" \"P1P0\" or \"P1P1\"");
/****************************************************************/
/* Create a search tree based on the bounding boxes */
/* Instanciate the intersector and initialise the result vector */
std::transform(tmp,tmp+SPACEDIM,quadOut+3*SPACEDIM,std::bind2nd(std::multiplies<double>(),0.5));
}
+ /*!
+ * This method builds a potentially non-convex polygon cell built with the first point of 'triIn' the barycenter of two edges starting or ending with
+ * the first point of 'triIn' and the barycenter of 'triIn'.
+ *
+ * @param triIn is a 6 doubles array in full interlace mode, that represents a triangle.
+ * @param quadOut is a 8 doubles array filled after the following call.
+ */
+ template<int SPACEDIM>
+ inline void fillDualCellOfPolyg(const double *polygIn, int nPtsPolygonIn, double *polygOut)
+ {
+ //1st point
+ std::copy(polygIn,polygIn+SPACEDIM,polygOut);
+ std::transform(polygIn,polygIn+SPACEDIM,polygIn+SPACEDIM,polygOut+SPACEDIM,std::plus<double>());
+ //2nd point
+ std::transform(polygOut+SPACEDIM,polygOut+2*SPACEDIM,polygOut+SPACEDIM,std::bind2nd(std::multiplies<double>(),0.5));
+ double tmp[SPACEDIM];
+ //
+ for(int i=0;i<nPtsPolygonIn-2;i++)
+ {
+ std::transform(polygIn,polygIn+SPACEDIM,polygIn+(i+2)*SPACEDIM,tmp,std::plus<double>());
+ std::transform(tmp,tmp+SPACEDIM,polygOut+(2*i+3)*SPACEDIM,std::bind2nd(std::multiplies<double>(),0.5));
+ std::transform(polygIn+(i+1)*SPACEDIM,polygIn+(i+2)*SPACEDIM,tmp,tmp,std::plus<double>());
+ std::transform(tmp,tmp+SPACEDIM,polygOut+(2*i+2)*SPACEDIM,std::bind2nd(std::multiplies<double>(),1./3.));
+ }
+ }
+
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
/* calcul les coordonnées du barycentre d'un polygone */
/* le vecteur en entrée est constitué des coordonnées */
PlanarIntersector.txx \
PlanarIntersectorP0P0.hxx \
PlanarIntersectorP0P0.txx \
+PlanarIntersectorP0P1.hxx \
+PlanarIntersectorP0P1.txx \
+PlanarIntersectorP1P0.hxx \
+PlanarIntersectorP1P0.txx \
+PlanarIntersectorP1P1.hxx \
+PlanarIntersectorP1P1.txx \
PolygonAlgorithms.hxx \
PolygonAlgorithms.txx \
PolyhedronIntersector.hxx \
#define __PLANARINTERSECTOR_HXX__
#include "TargetIntersector.hxx"
+#include "NormalizedUnstructuredMesh.hxx"
namespace INTERP_KERNEL
{
int projectionThis(double *Coords_A, double *Coords_B, int nb_NodesA, int nb_NodesB);
void getRealTargetCoordinates(ConnType icellT, std::vector<double>& coordsT);
void getRealSourceCoordinates(ConnType icellS, std::vector<double>& coordsS);
+ void getRealTargetCoordinatesPermute(ConnType icellT, int offset, std::vector<double>& coordsT);
+ void getRealSourceCoordinatesPermute(ConnType icellS, int offset, std::vector<double>& coordsS);
void getRealCoordinates(ConnType icellT, ConnType icellS, ConnType nbNodesT, ConnType nbNodesS, std::vector<double>& coordsT, std::vector<double>& coordsS, int& orientation);
static int projection(double *Coords_A, double *Coords_B,
int nb_NodesA, int nb_NodesB, double epsilon, double median_plane, bool do_rotate);
for(int idim=0; idim<SPACEDIM; idim++)
coordsS[SPACEDIM*iS+idim]=_coordsS[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+iS)])+idim];
}
+
+ /*!
+ * @param icellT id in target mesh in format of MyMeshType.
+ * @param offset is a value in C format that indicates the number of circular permutation.
+ * @param coordsT output val that stores coordinates of the target cell automatically resized to the right length.
+ */
+ template<class MyMeshType, class MyMatrix>
+ void PlanarIntersector<MyMeshType,MyMatrix>::getRealTargetCoordinatesPermute(ConnType icellT, int offset, std::vector<double>& coordsT)
+ {
+ int nbNodesT=_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)+1]-_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)];
+ coordsT.resize(SPACEDIM*nbNodesT);
+ for (ConnType iTTmp=0; iTTmp<nbNodesT; iTTmp++)
+ {
+ ConnType iT=(iTTmp+offset)%nbNodesT;
+ for(int idim=0; idim<SPACEDIM; idim++)
+ coordsT[SPACEDIM*iTTmp+idim]=_coordsT[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectT[OTT<ConnType,numPol>::conn2C(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]+iT)])+idim];
+ }
+ }
+
+ /*!
+ * @param icellS id in source mesh in format of MyMeshType.
+ * @param offset is a value in C format that indicates the number of circular permutation.
+ * @param coordsS output val that stores coordinates of the source cell automatically resized to the right length.
+ */
+ template<class MyMeshType, class MyMatrix>
+ void PlanarIntersector<MyMeshType,MyMatrix>::getRealSourceCoordinatesPermute(ConnType icellS, int offset, std::vector<double>& coordsS)
+ {
+ int nbNodesS=_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)+1]-_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)];
+ coordsS.resize(SPACEDIM*nbNodesS);
+ for (ConnType iSTmp=0; iSTmp<nbNodesS; iSTmp++)
+ {
+ ConnType iS=(iSTmp+offset)%nbNodesS;
+ for(int idim=0; idim<SPACEDIM; idim++)
+ coordsS[SPACEDIM*iSTmp+idim]=_coordsS[SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+iS)])+idim];
+ }
+ }
/*!
* @param icellT id in target mesh in format of MyMeshType.
std::vector<double> targetCellCoords;
int orientation=1;
PlanarIntersector<MyMeshType,MyMatrix>::getRealTargetCoordinates(OTT<ConnType,numPol>::indFC(icellT),targetCellCoords);
- NormalizedCellType tT=PlanarIntersector<MyMeshType,MyMatrix>::_meshT.getTypeOfElement(icellT);
+ NormalizedCellType tT=PlanarIntersector<MyMeshType,MyMatrix>::_meshT.getTypeOfElement(OTT<ConnType,numPol>::indFC(icellT));
bool isTargetQuad=CellModel::getCellModel(tT).isQuadratic();
typename MyMatrix::value_type& resRow=res[icellT];
for(typename std::vector<ConnType>::const_iterator iter=icellsS.begin();iter!=icellsS.end();iter++)
--- /dev/null
+// Copyright (C) 2007-2008 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 __PLANARINTERSECTORP1P1_HXX__
+#define __PLANARINTERSECTORP1P1_HXX__
+
+#include "PlanarIntersector.hxx"
+
+namespace INTERP_KERNEL
+{
+ template<class MyMeshType, class MyMatrix, class ConcreteP1P1Intersector>
+ class PlanarIntersectorP1P1 : public PlanarIntersector<MyMeshType,MyMatrix>
+ {
+ public:
+ static const int SPACEDIM=MyMeshType::MY_SPACEDIM;
+ static const int MESHDIM=MyMeshType::MY_MESHDIM;
+ typedef typename MyMeshType::MyConnType ConnType;
+ static const NumberingPolicy numPol=MyMeshType::My_numPol;
+ protected:
+ PlanarIntersectorP1P1(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double medianPlane, bool doRotate, int orientation, int printLevel);
+ public:
+ void intersectCells(ConnType icellT, const std::vector<ConnType>& icellsS, MyMatrix& res);
+ int getNumberOfRowsOfResMatrix() const;
+ int getNumberOfColsOfResMatrix() const;
+
+ double intersectGeometryGeneral(const std::vector<double>& targetCoords, const std::vector<double>& sourceCoords) { return asLeaf().intersectGeometryGeneral(targetCoords,sourceCoords); }
+ protected:
+ ConcreteP1P1Intersector& asLeaf() { return static_cast<ConcreteP1P1Intersector&>(*this); }
+ };
+}
+
+#endif
--- /dev/null
+// Copyright (C) 2007-2008 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 __PLANARINTERSECTORP1P1_TXX__
+#define __PLANARINTERSECTORP1P1_TXX__
+
+#include "PlanarIntersectorP1P1.hxx"
+#include "InterpolationUtils.hxx"
+#include "CellModel.hxx"
+
+namespace INTERP_KERNEL
+{
+ template<class MyMeshType, class MyMatrix, class ConcreteP1P1Intersector>
+ PlanarIntersectorP1P1<MyMeshType,MyMatrix,ConcreteP1P1Intersector>::PlanarIntersectorP1P1(const MyMeshType& meshT, const MyMeshType& meshS,
+ double dimCaracteristic, double precision, double medianPlane,
+ bool doRotate, int orientation, int printLevel):
+ PlanarIntersector<MyMeshType,MyMatrix>(meshT,meshS,dimCaracteristic,precision,medianPlane,doRotate,orientation,printLevel)
+ {
+ }
+
+ template<class MyMeshType, class MyMatrix, class ConcreteP1P1Intersector>
+ int PlanarIntersectorP1P1<MyMeshType,MyMatrix,ConcreteP1P1Intersector>::getNumberOfRowsOfResMatrix() const
+ {
+ return PlanarIntersector<MyMeshType,MyMatrix>::_meshT.getNumberOfNodes();
+ }
+
+ template<class MyMeshType, class MyMatrix, class ConcreteP1P1Intersector>
+ int PlanarIntersectorP1P1<MyMeshType,MyMatrix,ConcreteP1P1Intersector>::getNumberOfColsOfResMatrix() const
+ {
+ return PlanarIntersector<MyMeshType,MyMatrix>::_meshS.getNumberOfNodes();
+ }
+
+ /*!
+ * This methods split on the fly, into triangles in order to compute dual mesh of target cell (with icellT id in target mesh in C mode).
+ */
+ template<class MyMeshType, class MyMatrix, class ConcreteP1P1Intersector>
+ void PlanarIntersectorP1P1<MyMeshType,MyMatrix,ConcreteP1P1Intersector>::intersectCells(ConnType icellT, const std::vector<ConnType>& icellsS, MyMatrix& res)
+ {
+ int nbNodesT=PlanarIntersector<MyMeshType,MyMatrix>::_connIndexT[icellT+1]-PlanarIntersector<MyMeshType,MyMatrix>::_connIndexT[icellT];
+ int orientation=1;
+ const ConnType *startOfCellNodeConn=PlanarIntersector<MyMeshType,MyMatrix>::_connectT+OTT<ConnType,numPol>::conn2C(PlanarIntersector<MyMeshType,MyMatrix>::_connIndexT[icellT]);
+ std::vector<double> polygT;
+ PlanarIntersector<MyMeshType,MyMatrix>::getRealTargetCoordinates(OTT<ConnType,numPol>::indFC(icellT),polygT);
+ for(int nodeIdT=0;nodeIdT<nbNodesT;nodeIdT++)
+ {
+ ConnType curNodeTInCmode=OTT<ConnType,numPol>::coo2C(startOfCellNodeConn[nodeIdT]);
+ PlanarIntersector<MyMeshType,MyMatrix>::getRealTargetCoordinatesPermute(OTT<ConnType,numPol>::indFC(icellT),nodeIdT,polygT);
+ std::vector<double> polygDualT(SPACEDIM*2*(nbNodesT-1));
+ fillDualCellOfPolyg<SPACEDIM>(&polygT[0],polygT.size()/SPACEDIM,&polygDualT[0]);
+ typename MyMatrix::value_type& resRow=res[curNodeTInCmode];
+ for(typename std::vector<ConnType>::const_iterator iter=icellsS.begin();iter!=icellsS.end();iter++)
+ {
+ int iS=*iter;
+ int nbNodesS=PlanarIntersector<MyMeshType,MyMatrix>::_connIndexS[iS+1]-PlanarIntersector<MyMeshType,MyMatrix>::_connIndexS[iS];
+ const ConnType *startOfCellNodeConnS=PlanarIntersector<MyMeshType,MyMatrix>::_connectS+OTT<ConnType,numPol>::conn2C(PlanarIntersector<MyMeshType,MyMatrix>::_connIndexS[iS]);
+ for(int nodeIdS=0;nodeIdS<nbNodesS;nodeIdS++)
+ {
+ ConnType curNodeSInCmode=OTT<ConnType,numPol>::coo2C(startOfCellNodeConnS[nodeIdS]);
+ std::vector<double> polygS;
+ PlanarIntersector<MyMeshType,MyMatrix>::getRealSourceCoordinatesPermute(OTT<ConnType,numPol>::indFC(iS),nodeIdS,polygS);
+ std::vector<double> polygDualS(SPACEDIM*2*(nbNodesS-1));
+ fillDualCellOfPolyg<SPACEDIM>(&polygS[0],polygS.size()/SPACEDIM,&polygDualS[0]);
+ std::vector<double> polygDualTTmp(polygDualT);
+ if(SPACEDIM==3)
+ orientation=PlanarIntersector<MyMeshType,MyMatrix>::projectionThis(&polygDualS[0],&polygDualTTmp[0],polygDualS.size()/SPACEDIM,polygDualT.size()/SPACEDIM);
+ double surf=orientation*intersectGeometryGeneral(polygDualTTmp,polygDualS);
+ //filtering out zero surfaces and badly oriented surfaces
+ // _orientation = -1,0,1
+ // -1 : the intersection is taken into account if target and cells have different orientation
+ // 0 : the intersection is always taken into account
+ // 1 : the intersection is taken into account if target and cells have the same orientation
+ if (( surf > 0.0 && PlanarIntersector<MyMeshType,MyMatrix>::_orientation >=0 ) || ( surf < 0.0 && PlanarIntersector<MyMeshType,MyMatrix>::_orientation <=0 ))
+ {
+ typename MyMatrix::value_type::const_iterator iterRes=resRow.find(OTT<ConnType,numPol>::indFC(curNodeSInCmode));
+ if(iterRes==resRow.end())
+ resRow.insert(std::make_pair(OTT<ConnType,numPol>::indFC(curNodeSInCmode),surf));
+ else
+ {
+ double val=(*iterRes).second+surf;
+ resRow.erase(OTT<ConnType,numPol>::indFC(curNodeSInCmode));
+ resRow.insert(std::make_pair(OTT<ConnType,numPol>::indFC(curNodeSInCmode),val));
+ }
+ }
+ }
+ }
+ }
+ }
+}
+
+#endif
else
{
double val=(*iterRes).second+volume;
- resRow.erase(OTT<ConnType,numPol>::indFC(*iterCellS));
+ resRow.erase(OTT<ConnType,numPol>::indFC(sourceNode));
resRow.insert(std::make_pair(OTT<ConnType,numPol>::indFC(sourceNode),val));
}
}
const int tab[3][2]={{1,2},{3,2},{1,3}};
const int *curTab=tab[case1];
double pt0[3]; pt0[0]=(tmp[curTab[case2]][0]+tmp[0][0])/2.; pt0[1]=(tmp[curTab[case2]][1]+tmp[0][1])/2.; pt0[2]=(tmp[curTab[case2]][2]+tmp[0][2])/2.;
- double pt1[3]; pt1[0]=(tmp[0][0]+tmp[curTab[0]][0]+tmp[curTab[1]][0])/3.; pt1[1]=(tmp[0][1]+tmp[curTab[0]][1]+tmp[curTab[1]][1])/3.; pt1[1]=(tmp[0][2]+tmp[curTab[0]][2]+tmp[curTab[1]][2])/3.;
+ double pt1[3]; pt1[0]=(tmp[0][0]+tmp[curTab[0]][0]+tmp[curTab[1]][0])/3.; pt1[1]=(tmp[0][1]+tmp[curTab[0]][1]+tmp[curTab[1]][1])/3.; pt1[2]=(tmp[0][2]+tmp[curTab[0]][2]+tmp[curTab[1]][2])/3.;
double pt2[3]; pt2[0]=(tmp[0][0]+tmp[1][0]+tmp[2][0]+tmp[3][0])/4.; pt2[1]=(tmp[0][1]+tmp[1][1]+tmp[2][1]+tmp[3][1])/4.; pt2[2]=(tmp[0][2]+tmp[1][2]+tmp[2][2]+tmp[3][2])/4.;
std::copy(pt1,pt1+3,output+case2*3);
std::copy(pt0,pt0+3,output+(abs(case2-1))*3);
#include "PlanarIntersectorP0P0.hxx"
#include "PlanarIntersectorP0P1.hxx"
#include "PlanarIntersectorP1P0.hxx"
+#include "PlanarIntersectorP1P1.hxx"
namespace INTERP_KERNEL
{
double dimCaracteristic, double precision, double medianPlane, int orientation, int printLevel);
double intersectGeometry(ConnType icellT, ConnType icellS, ConnType nbNodesT, ConnType nbNodesS);
double intersectGeometryWithQuadrangle(const double *quadrangle, const std::vector<double>& sourceCoords, bool isSourceQuad);
+ double intersectGeometryGeneral(const std::vector<double>& targetCoords, const std::vector<double>& sourceCoords);
};
}
return result;
}
+
+ template<class MyMeshType, class MyMatrix, template <class MeshType, class TheMatrix, class ThisIntersector> class InterpType>
+ double TriangulationIntersector<MyMeshType,MyMatrix,InterpType>::intersectGeometryGeneral(const std::vector<double>& targetCoords, const std::vector<double>& sourceCoords)
+ {
+ double result = 0.;
+ ConnType nbNodesS=sourceCoords.size()/SPACEDIM;
+ ConnType nbNodesT=targetCoords.size()/SPACEDIM;
+ //Compute the intersection area
+ double area[SPACEDIM];
+ for(ConnType iT = 1; iT<nbNodesT-1; iT++)
+ {
+ for(ConnType iS = 1; iS<nbNodesS-1; iS++)
+ {
+ std::vector<double> inter;
+ INTERP_KERNEL::intersec_de_triangle(&targetCoords[0],&targetCoords[SPACEDIM*iT],&targetCoords[SPACEDIM*(iT+1)],
+ &sourceCoords[0],&sourceCoords[SPACEDIM*iS],&sourceCoords[SPACEDIM*(iS+1)],
+ inter, PlanarIntersector<MyMeshType,MyMatrix>::_dim_caracteristic,
+ PlanarIntersector<MyMeshType,MyMatrix>::_precision);
+ ConnType nb_inter=((ConnType)inter.size())/2;
+ if(nb_inter >3) inter=reconstruct_polygon(inter);
+ for(ConnType i = 1; i<nb_inter-1; i++)
+ {
+ INTERP_KERNEL::crossprod<2>(&inter[0],&inter[2*i],&inter[2*(i+1)],area);
+ result +=0.5*fabs(area[0]);
+ }
+ }
+ }
+ return result;
+ }
}
#endif
//
#include "MEDCouplingField.hxx"
#include "MEDCouplingMesh.hxx"
+#include "MEDCouplingFieldDiscretization.hxx"
using namespace ParaMEDMEM;
updateTimeWith(*_mesh);
}
+TypeOfField MEDCouplingField::getEntity() const
+{
+ return _type->getEnum();
+}
+
void MEDCouplingField::setMesh(MEDCouplingMesh *mesh)
{
if(mesh!=_mesh)
{
if(_mesh)
_mesh->decrRef();
+ delete _type;
+}
+
+MEDCouplingField::MEDCouplingField(TypeOfField type):_mesh(0),_type(MEDCouplingFieldDiscretization::New(type))
+{
}
MEDCouplingField::MEDCouplingField(const MEDCouplingField& other):_name(other._name),_desc(other._name),
- _time(other._time),_dt(other._dt),_it(other._it),
- _mesh(0),_type(other._type)
+ _mesh(0),_type(other._type->clone())
{
if(other._mesh)
{
namespace ParaMEDMEM
{
class MEDCouplingMesh;
+ class MEDCouplingFieldDiscretization;
class MEDCOUPLING_EXPORT MEDCouplingField : public RefCountObject
{
public:
virtual void checkCoherency() const throw(INTERP_KERNEL::Exception) = 0;
void setMesh(ParaMEDMEM::MEDCouplingMesh *mesh);
- void setTime(double val) { _time=val; }
- double getTime() const { return _time; }
- void setDtIt(int dt, int it) { _dt=dt; _it=it; }
- void getDtIt(int& dt, int& it) { dt=_dt; it=_it; }
ParaMEDMEM::MEDCouplingMesh *getMesh() const { return _mesh; }
void setName(const char *name) { _name=name; }
void setDescription(const char *desc) { _desc=desc; }
const char *getName() const { return _name.c_str(); }
- TypeOfField getEntity() const { return _type; }
+ TypeOfField getEntity() const;
protected:
void updateTime();
protected:
- MEDCouplingField(TypeOfField type):_time(0.),_dt(-1),_it(-1),_mesh(0),_type(type) { }
+ MEDCouplingField(TypeOfField type);
MEDCouplingField(const MEDCouplingField& other);
virtual ~MEDCouplingField();
protected:
std::string _name;
std::string _desc;
- double _time;
- int _dt;
- int _it;
MEDCouplingMesh *_mesh;
- const TypeOfField _type;
+ MEDCouplingFieldDiscretization *_type;
};
}
--- /dev/null
+// Copyright (C) 2007-2008 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 "MEDCouplingFieldDiscretization.hxx"
+#include "MEDCouplingMesh.hxx"
+#include "MEDCouplingFieldDouble.hxx"
+
+using namespace ParaMEDMEM;
+
+const char MEDCouplingFieldDiscretizationP0::REPR[]="P0";
+
+const TypeOfField MEDCouplingFieldDiscretizationP0::TYPE=ON_CELLS;
+
+const char MEDCouplingFieldDiscretizationP1::REPR[]="P1";
+
+const TypeOfField MEDCouplingFieldDiscretizationP1::TYPE=ON_NODES;
+
+MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretization::New(TypeOfField type)
+{
+ switch(type)
+ {
+ case MEDCouplingFieldDiscretizationP0::TYPE:
+ return new MEDCouplingFieldDiscretizationP0;
+ case MEDCouplingFieldDiscretizationP1::TYPE:
+ return new MEDCouplingFieldDiscretizationP1;
+ default:
+ throw INTERP_KERNEL::Exception("Choosen discretization is not implemented yet.");
+ }
+}
+
+TypeOfField MEDCouplingFieldDiscretizationP0::getEnum() const
+{
+ return TYPE;
+}
+
+MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretizationP0::clone() const
+{
+ return new MEDCouplingFieldDiscretizationP0;
+}
+
+const char *MEDCouplingFieldDiscretizationP0::getStringRepr() const
+{
+ return REPR;
+}
+
+int MEDCouplingFieldDiscretizationP0::getNumberOfTuples(const MEDCouplingMesh *mesh) const
+{
+ return mesh->getNumberOfCells();
+}
+
+void MEDCouplingFieldDiscretizationP0::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception)
+{
+ if(mesh->getNumberOfCells()!=da->getNumberOfTuples())
+ {
+ std::ostringstream message;
+ message << "Field on cells invalid because there are " << mesh->getNumberOfCells();
+ message << " cells in mesh and " << da->getNumberOfTuples() << " tuples in field !";
+ throw INTERP_KERNEL::Exception(message.str().c_str());
+ }
+}
+
+MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationP0::getWeightingField(const MEDCouplingMesh *mesh) const
+{
+ return mesh->getMeasureField();
+}
+
+TypeOfField MEDCouplingFieldDiscretizationP1::getEnum() const
+{
+ return TYPE;
+}
+
+MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretizationP1::clone() const
+{
+ return new MEDCouplingFieldDiscretizationP1;
+}
+
+const char *MEDCouplingFieldDiscretizationP1::getStringRepr() const
+{
+ return REPR;
+}
+
+int MEDCouplingFieldDiscretizationP1::getNumberOfTuples(const MEDCouplingMesh *mesh) const
+{
+ return mesh->getNumberOfNodes();
+}
+
+void MEDCouplingFieldDiscretizationP1::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception)
+{
+ if(mesh->getNumberOfNodes()!=da->getNumberOfTuples())
+ {
+ std::ostringstream message;
+ message << "Field on nodes invalid because there are " << mesh->getNumberOfNodes();
+ message << " cells in mesh and " << da->getNumberOfTuples() << " tuples in field !";
+ throw INTERP_KERNEL::Exception(message.str().c_str());
+ }
+}
+
+MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationP1::getWeightingField(const MEDCouplingMesh *mesh) const
+{
+ //not implemented yet.
+ //Dual mesh to build
+ return 0;
+}
--- /dev/null
+// Copyright (C) 2007-2008 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 __MEDCOUPLINGFIELDDISCRETIZATION_HXX__
+#define __MEDCOUPLINGFIELDDISCRETIZATION_HXX__
+
+#include "MEDCoupling.hxx"
+#include "RefCountObject.hxx"
+#include "InterpKernelException.hxx"
+
+namespace ParaMEDMEM
+{
+ class MEDCouplingMesh;
+ class DataArrayDouble;
+ class MEDCouplingFieldDouble;
+
+ class MEDCOUPLING_EXPORT MEDCouplingFieldDiscretization
+ {
+ public:
+ static MEDCouplingFieldDiscretization *New(TypeOfField type);
+ virtual TypeOfField getEnum() const = 0;
+ virtual MEDCouplingFieldDiscretization *clone() const = 0;
+ virtual const char *getStringRepr() const = 0;
+ virtual int getNumberOfTuples(const MEDCouplingMesh *mesh) const = 0;
+ virtual void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception) = 0;
+ virtual MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh) const = 0;
+ };
+
+ class MEDCOUPLING_EXPORT MEDCouplingFieldDiscretizationP0 : public MEDCouplingFieldDiscretization
+ {
+ public:
+ TypeOfField getEnum() const;
+ MEDCouplingFieldDiscretization *clone() const;
+ const char *getStringRepr() const;
+ int getNumberOfTuples(const MEDCouplingMesh *mesh) const;
+ void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception);
+ MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh) const;
+ public:
+ static const char REPR[];
+ static const TypeOfField TYPE;
+ };
+
+ class MEDCOUPLING_EXPORT MEDCouplingFieldDiscretizationP1 : public MEDCouplingFieldDiscretization
+ {
+ public:
+ TypeOfField getEnum() const;
+ MEDCouplingFieldDiscretization *clone() const;
+ const char *getStringRepr() const;
+ int getNumberOfTuples(const MEDCouplingMesh *mesh) const;
+ void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception);
+ MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh) const;
+ public:
+ static const char REPR[];
+ static const TypeOfField TYPE;
+ };
+}
+
+#endif
//
#include "MEDCouplingFieldDouble.hxx"
#include "MEDCouplingMesh.hxx"
+#include "MEDCouplingTimeDiscretization.hxx"
+#include "MEDCouplingFieldDiscretization.hxx"
#include <sstream>
using namespace ParaMEDMEM;
-MEDCouplingFieldDouble *MEDCouplingFieldDouble::New(TypeOfField type)
+MEDCouplingFieldDouble *MEDCouplingFieldDouble::New(TypeOfField type, TypeOfTimeDiscretization td)
{
- return new MEDCouplingFieldDouble(type);
+ return new MEDCouplingFieldDouble(type,td);
}
MEDCouplingFieldDouble *MEDCouplingFieldDouble::clone(bool recDeepCpy) const
return new MEDCouplingFieldDouble(*this,recDeepCpy);
}
-MEDCouplingFieldDouble::MEDCouplingFieldDouble(TypeOfField type):MEDCouplingField(type),_array(0)
+MEDCouplingFieldDouble::MEDCouplingFieldDouble(TypeOfField type, TypeOfTimeDiscretization td):MEDCouplingField(type),
+ _time_discr(MEDCouplingTimeDiscretization::New(td))
{
}
-MEDCouplingFieldDouble::MEDCouplingFieldDouble(const MEDCouplingFieldDouble& other, bool deepCpy):MEDCouplingField(other),_array(0)
+MEDCouplingFieldDouble::MEDCouplingFieldDouble(const MEDCouplingFieldDouble& other, bool deepCpy):MEDCouplingField(other),
+ _time_discr(other._time_discr->performCpy(deepCpy))
{
- if(other._array)
- _array=other._array->performCpy(deepCpy);
}
MEDCouplingFieldDouble::~MEDCouplingFieldDouble()
{
- if(_array)
- _array->decrRef();
+ delete _time_discr;
}
void MEDCouplingFieldDouble::checkCoherency() const throw(INTERP_KERNEL::Exception)
{
if(!_mesh)
throw INTERP_KERNEL::Exception("Field invalid because no mesh specified !");
- if(!_array)
+ if(!getArray())
throw INTERP_KERNEL::Exception("Field invalid because no values set !");
- if(_type==ON_CELLS)
- {
- if(_mesh->getNumberOfCells()!=_array->getNumberOfTuples())
- {
- std::ostringstream message;
- message << "Field on cells invalid because there are " << _mesh->getNumberOfCells();
- message << " cells in mesh and " << _array->getNumberOfTuples() << " tuples in field !";
- throw INTERP_KERNEL::Exception(message.str().c_str());
- }
- }
- else if(_type==ON_NODES)
- {
- if(_mesh->getNumberOfNodes()!=_array->getNumberOfTuples())
- {
- std::ostringstream message;
- message << "Field on nodes invalid because there are " << _mesh->getNumberOfNodes();
- message << " cells in mesh and " << _array->getNumberOfTuples() << " tuples in field !";
- throw INTERP_KERNEL::Exception(message.str().c_str());
- }
- }
- else
- throw INTERP_KERNEL::Exception("Field of undifined type !!!");
+ _type->checkCoherencyBetween(_mesh,getArray());
+}
+
+double MEDCouplingFieldDouble::accumulate(int compId) const
+{
+ const double *ptr=getArray()->getConstPointer();
+ int nbTuple=getArray()->getNumberOfTuples();
+ int nbComps=getArray()->getNumberOfComponents();
+ double ret=0.;
+ for(int i=0;i<nbTuple;i++)
+ ret+=ptr[i*nbComps+compId];
+ return ret;
+}
+
+double MEDCouplingFieldDouble::measureAccumulate(int compId) const
+{
+ if(!_mesh)
+ throw INTERP_KERNEL::Exception("No mesh underlying this field to perform measureAccumulate");
+ MEDCouplingFieldDouble *weight=_type->getWeightingField(_mesh);
+ const double *ptr=weight->getArray()->getConstPointer();
+ int nbOfValues=weight->getArray()->getNbOfElems();
+ double ret=0.;
+ for (int i=0; i<nbOfValues; i++)
+ ret+=getIJ(i,compId)*ptr[i];
+ weight->decrRef();
+ return ret;
+}
+
+void MEDCouplingFieldDouble::getValueOn(const double *spaceLoc, double *res) const throw(INTERP_KERNEL::Exception)
+{
+ _time_discr->checkNoTimePresence();
+}
+
+void MEDCouplingFieldDouble::getValueOn(const double *spaceLoc, double time, double *res) const throw(INTERP_KERNEL::Exception)
+{
+ _time_discr->checkTimePresence(time);
}
void MEDCouplingFieldDouble::applyLin(double a, double b, int compoId)
{
- double *ptr=_array->getPointer();
+ double *ptr=getArray()->getPointer();
ptr+=compoId;
- int nbOfComp=_array->getNumberOfComponents();
- int nbOfTuple=_array->getNumberOfTuples();
+ int nbOfComp=getArray()->getNumberOfComponents();
+ int nbOfTuple=getArray()->getNumberOfTuples();
for(int i=0;i<nbOfTuple;i++,ptr+=nbOfComp)
*ptr=a*(*ptr)+b;
}
int MEDCouplingFieldDouble::getNumberOfComponents() const
{
- return _array->getNumberOfComponents();
+ return getArray()->getNumberOfComponents();
}
int MEDCouplingFieldDouble::getNumberOfTuples() const throw(INTERP_KERNEL::Exception)
{
if(!_mesh)
throw INTERP_KERNEL::Exception("Impossible to retrieve number of tuples because no mesh specified !");
- if(_type==ON_CELLS)
- return _mesh->getNumberOfCells();
- else if(_type==ON_NODES)
- return _mesh->getNumberOfNodes();
- else
- throw INTERP_KERNEL::Exception("Impossible to retrieve number of tuples because type of entity not implemented yet !");
+ return _type->getNumberOfTuples(_mesh);
}
void MEDCouplingFieldDouble::updateTime()
{
MEDCouplingField::updateTime();
- if(_array)
- updateTimeWith(*_array);
+ if(getArray())
+ updateTimeWith(*getArray());
}
void MEDCouplingFieldDouble::setArray(DataArrayDouble *array)
{
- if(array!=_array)
- {
- if(_array)
- _array->decrRef();
- _array=array;
- if(_array)
- _array->incrRef();
- declareAsNew();
- }
+ _time_discr->setArray(array,this);
}
#include "MEDCoupling.hxx"
#include "MEDCouplingField.hxx"
+#include "MEDCouplingTimeDiscretization.hxx"
#include "MemArray.hxx"
namespace ParaMEDMEM
class MEDCOUPLING_EXPORT MEDCouplingFieldDouble : public MEDCouplingField
{
public:
- static MEDCouplingFieldDouble *New(TypeOfField type);
+ static MEDCouplingFieldDouble *New(TypeOfField type, TypeOfTimeDiscretization td=NO_TIME);
MEDCouplingFieldDouble *clone(bool recDeepCpy) const;
void checkCoherency() const throw(INTERP_KERNEL::Exception);
- double getIJ(int tupleId, int compoId) const { return _array->getIJ(tupleId,compoId); }
+ void setTime(double val, int dt, int it) { _time_discr->setTime(val,dt,it); }
+ double getTime(int& dt, int& it) const { return _time_discr->getTime(dt,it); }
+ double getIJ(int tupleId, int compoId) const { return getArray()->getIJ(tupleId,compoId); }
void setArray(DataArrayDouble *array);
- DataArrayDouble *getArray() const { return _array; }
+ DataArrayDouble *getArray() const { return _time_discr->getArray(); }
+ double accumulate(int compId) const;
+ double measureAccumulate(int compId) const;
+ void getValueOn(const double *spaceLoc, double *res) const throw(INTERP_KERNEL::Exception);
+ void getValueOn(const double *spaceLoc, double time, double *res) const throw(INTERP_KERNEL::Exception);
//! \b temporary
void applyLin(double a, double b, int compoId);
int getNumberOfComponents() const;
int getNumberOfTuples() const throw(INTERP_KERNEL::Exception);
void updateTime();
private:
- MEDCouplingFieldDouble(TypeOfField type);
+ MEDCouplingFieldDouble(TypeOfField type, TypeOfTimeDiscretization td);
MEDCouplingFieldDouble(const MEDCouplingFieldDouble& other, bool deepCpy);
~MEDCouplingFieldDouble();
private:
- DataArrayDouble *_array;
+ MEDCouplingTimeDiscretization *_time_discr;
};
}
namespace ParaMEDMEM
{
+ class MEDCouplingFieldDouble;
+
class MEDCOUPLING_EXPORT MEDCouplingMesh : public RefCountObject
{
public:
virtual int getNumberOfNodes() const = 0;
virtual int getSpaceDimension() const = 0;
virtual int getMeshDimension() const = 0;
+ // tools
+ virtual void getBoundingBox(double *bbox) const = 0;
+ virtual MEDCouplingFieldDouble *getMeasureField() const = 0;
protected:
MEDCouplingMesh() { }
MEDCouplingMesh(const MEDCouplingMesh& other):_name(other._name) { }
boundingBox[SPACEDIM+i]=-std::numeric_limits<double>::max();
}
ParaMEDMEM::DataArrayDouble *array=_mesh->getCoords();
- const double *ptr=array->getPointer();
+ const double *ptr=array->getConstPointer();
int nbOfPts=array->getNbOfElems()/SPACEDIM;
for(int j=0;j<SPACEDIM;j++)
{
const double *MEDCouplingNormalizedUnstructuredMesh<SPACEDIM,MESHDIM>::getCoordinatesPtr() const
{
ParaMEDMEM::DataArrayDouble *array=_mesh->getCoords();
- return array->getPointer();
+ return array->getConstPointer();
}
template<int SPACEDIM,int MESHDIM>
_conn_for_interp=new int[initialConnSize-nbOfCell];
_conn_index_for_interp=new int[nbOfCell+1];
_conn_index_for_interp[0]=0;
- const int *work_conn=_mesh->getNodalConnectivity()->getPointer()+1;
- const int *work_conn_index=_mesh->getNodalConnectivityIndex()->getPointer();
+ const int *work_conn=_mesh->getNodalConnectivity()->getConstPointer()+1;
+ const int *work_conn_index=_mesh->getNodalConnectivityIndex()->getConstPointer();
int *work_conn_for_interp=_conn_for_interp;
int *work_conn_index_for_interp=_conn_index_for_interp;
for(int i=0;i<nbOfCell;i++)
--- /dev/null
+// Copyright (C) 2007-2008 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 "MEDCouplingPointSet.hxx"
+#include "MemArray.hxx"
+
+using namespace ParaMEDMEM;
+
+MEDCouplingPointSet::MEDCouplingPointSet():_coords(0)
+{
+}
+
+MEDCouplingPointSet::MEDCouplingPointSet(const MEDCouplingPointSet& other, bool deepCpy):MEDCouplingMesh(other),_coords(0)
+{
+ if(other._coords)
+ _coords=other._coords->performCpy(deepCpy);
+}
+
+MEDCouplingPointSet::~MEDCouplingPointSet()
+{
+ if(_coords)
+ _coords->decrRef();
+}
+
+int MEDCouplingPointSet::getNumberOfNodes() const
+{
+ if(_coords)
+ return _coords->getNumberOfTuples();
+ else
+ throw INTERP_KERNEL::Exception("Unable to get number of nodes because no coordinates specified !");
+}
+
+int MEDCouplingPointSet::getSpaceDimension() const
+{
+ if(_coords)
+ return _coords->getNumberOfComponents();
+ else
+ throw INTERP_KERNEL::Exception("Unable to get space dimension because no coordinates specified !");
+}
+
+void MEDCouplingPointSet::updateTime()
+{
+ if(_coords)
+ {
+ updateTimeWith(*_coords);
+ }
+}
+
+bool MEDCouplingPointSet::isStructured() const
+{
+ return false;
+}
+
+void MEDCouplingPointSet::setCoords(DataArrayDouble *coords)
+{
+ if( coords != _coords )
+ {
+ if (_coords)
+ _coords->decrRef();
+ _coords=coords;
+ if(_coords)
+ _coords->incrRef();
+ declareAsNew();
+ }
+}
+
+bool MEDCouplingPointSet::areCoordsEqual(const MEDCouplingPointSet& other, double prec) const
+{
+ return _coords->isEqual(other._coords,prec);
+}
+
+void MEDCouplingPointSet::getBoundingBox(double *bbox) const
+{
+ int dim=getSpaceDimension();
+ for (int idim=0; idim<dim; idim++)
+ {
+ bbox[idim*2]=std::numeric_limits<double>::max();
+ bbox[idim*2+1]=-std::numeric_limits<double>::max();
+ }
+ const double *coords=_coords->getConstPointer();
+ int nbnodes=getNumberOfNodes();
+ for (int i=0; i<nbnodes; i++)
+ {
+ for (int idim=0; idim<dim;idim++)
+ {
+ if ( bbox[idim*2] > coords[i*dim+idim] )
+ {
+ bbox[idim*2] = coords[i*dim+idim] ;
+ }
+ if ( bbox[idim*2+1] < coords[i*dim+idim] )
+ {
+ bbox[idim*2+1] = coords[i*dim+idim] ;
+ }
+ }
+ }
+}
+
+// =============================================
+// Intersect Bounding Box given 2 Bounding Boxes
+// =============================================
+bool MEDCouplingPointSet::intersectsBoundingBox(const double* bb1, const double* bb2, int dim, double eps)
+{
+ double bbtemp[2*dim];
+ double deltamax=0.0;
+
+ for (int i=0; i< dim; i++)
+ {
+ double delta = bb1[2*i+1]-bb1[2*i];
+ if ( delta > deltamax )
+ {
+ deltamax = delta ;
+ }
+ }
+ for (int i=0; i<dim; i++)
+ {
+ bbtemp[i*2]=bb1[i*2]-deltamax*eps;
+ bbtemp[i*2+1]=bb1[i*2+1]+deltamax*eps;
+ }
+
+ for (int idim=0; idim < dim; idim++)
+ {
+ bool intersects = (bbtemp[idim*2]<bb2[idim*2+1])
+ && (bb2[idim*2]<bbtemp[idim*2+1]) ;
+ if (!intersects) return false;
+ }
+ return true;
+}
--- /dev/null
+// Copyright (C) 2007-2008 CEA/DEN, EDF R&D
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+#ifndef __PARAMEDMEM_MEDCOUPLINGPOINTSET_HXX__
+#define __PARAMEDMEM_MEDCOUPLINGPOINTSET_HXX__
+
+#include "MEDCoupling.hxx"
+#include "MEDCouplingMesh.hxx"
+
+#include <vector>
+
+namespace ParaMEDMEM
+{
+ class DataArrayInt;
+ class DataArrayDouble;
+
+ class MEDCOUPLING_EXPORT MEDCouplingPointSet : public MEDCouplingMesh
+ {
+ protected:
+ MEDCouplingPointSet();
+ MEDCouplingPointSet(const MEDCouplingPointSet& other, bool deepCpy);
+ ~MEDCouplingPointSet();
+ public:
+ void updateTime();
+ bool isStructured() const;
+ int getNumberOfNodes() const;
+ int getSpaceDimension() const;
+ void setCoords(DataArrayDouble *coords);
+ DataArrayDouble *getCoords() const { return _coords; }
+ bool areCoordsEqual(const MEDCouplingPointSet& other, double prec) const;
+ void getBoundingBox(double *bbox) const;
+ virtual MEDCouplingPointSet *buildPartOfMySelf(const int *start, const int *end, bool keepCoords) const = 0;
+ //! size of returned tinyInfo must be always the same.
+ virtual void getTinySerializationInformation(std::vector<int>& tinyInfo) const = 0;
+ virtual void resizeForSerialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2) = 0;
+ virtual void serialize(DataArrayInt *&a1, DataArrayDouble *&a2) = 0;
+ virtual MEDCouplingPointSet *buildObjectFromUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2) = 0;
+ virtual void giveElemsInBoundingBox(const double *bbox, double eps, std::vector<int>& elems) = 0;
+ protected:
+ static bool intersectsBoundingBox(const double* bb1, const double* bb2, int dim, double eps);
+ protected:
+ DataArrayDouble *_coords;
+ };
+}
+
+#endif
--- /dev/null
+// Copyright (C) 2007-2008 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 "MEDCouplingRMesh.hxx"
+#include "MemArray.hxx"
+
+using namespace ParaMEDMEM;
+
+MEDCouplingRMesh::MEDCouplingRMesh():_x_array(0),_y_array(0),_z_array(0)
+{
+}
+
+MEDCouplingRMesh::~MEDCouplingRMesh()
+{
+ if(_x_array)
+ _x_array->decrRef();
+ if(_y_array)
+ _y_array->decrRef();
+ if(_z_array)
+ _z_array->decrRef();
+}
+
+MEDCouplingRMesh *MEDCouplingRMesh::New()
+{
+ return new MEDCouplingRMesh;
+}
+
+void MEDCouplingRMesh::updateTime()
+{
+ if(_x_array)
+ updateTimeWith(*_x_array);
+ if(_y_array)
+ updateTimeWith(*_y_array);
+ if(_z_array)
+ updateTimeWith(*_z_array);
+}
+
+bool MEDCouplingRMesh::isEqual(const MEDCouplingMesh *other, double prec) const
+{
+ const MEDCouplingRMesh *otherC=dynamic_cast<const MEDCouplingRMesh *>(other);
+ if(!otherC)
+ return false;
+ return true;
+}
+
+void MEDCouplingRMesh::checkCoherency() const throw(INTERP_KERNEL::Exception)
+{
+ const char msg0[]="Invalid ";
+ const char msg1[]=" array ! Must contain more than 1 element.";
+ if(_x_array)
+ if(_x_array->getNbOfElems()<2)
+ {
+ std::ostringstream os; os << msg0 << 'X' << msg1;
+ throw INTERP_KERNEL::Exception(os.str().c_str());
+ }
+ if(_y_array)
+ if(_y_array->getNbOfElems()<2)
+ {
+ std::ostringstream os; os << msg0 << 'Y' << msg1;
+ throw INTERP_KERNEL::Exception(os.str().c_str());
+ }
+ if(_z_array)
+ if(_z_array->getNbOfElems()<2)
+ {
+ std::ostringstream os; os << msg0 << 'Z' << msg1;
+ throw INTERP_KERNEL::Exception(os.str().c_str());
+ }
+}
+
+bool MEDCouplingRMesh::isStructured() const
+{
+ return true;
+}
+
+int MEDCouplingRMesh::getNumberOfCells() const
+{
+ int ret=1;
+ if(_x_array)
+ ret*=_x_array->getNbOfElems()-1;
+ if(_y_array)
+ ret*=_y_array->getNbOfElems()-1;
+ if(_z_array)
+ ret*=_z_array->getNbOfElems()-1;
+ return ret;
+}
+
+int MEDCouplingRMesh::getNumberOfNodes() const
+{
+ int ret=1;
+ if(_x_array)
+ ret*=_x_array->getNbOfElems();
+ if(_y_array)
+ ret*=_y_array->getNbOfElems();
+ if(_z_array)
+ ret*=_z_array->getNbOfElems();
+ return ret;
+}
+
+int MEDCouplingRMesh::getSpaceDimension() const
+{
+ int ret=0;
+ if(_x_array)
+ ret++;
+ if(_y_array)
+ ret++;
+ if(_z_array)
+ ret++;
+ return ret;
+}
+
+int MEDCouplingRMesh::getMeshDimension() const
+{
+ int ret=0;
+ if(_x_array)
+ ret++;
+ if(_y_array)
+ ret++;
+ if(_z_array)
+ ret++;
+ return ret;
+}
+
+DataArrayDouble *MEDCouplingRMesh::getCoordsAt(int i) const throw(INTERP_KERNEL::Exception)
+{
+ switch(i)
+ {
+ case 0:
+ return _x_array;
+ case 1:
+ return _y_array;
+ case 2:
+ return _z_array;
+ default:
+ throw INTERP_KERNEL::Exception("Invalid rank specified must be 0 or 1 or 2.");
+ }
+}
+
+void MEDCouplingRMesh::setCoords(DataArrayDouble *coordsX, DataArrayDouble *coordsY, DataArrayDouble *coordsZ)
+{
+ if(_x_array)
+ _x_array->decrRef();
+ _x_array=coordsX;
+ if(_x_array)
+ _x_array->incrRef();
+ if(_y_array)
+ _y_array->decrRef();
+ _y_array=coordsY;
+ if(_y_array)
+ _y_array->incrRef();
+ if(_z_array)
+ _z_array->decrRef();
+ _z_array=coordsZ;
+ if(_z_array)
+ _z_array->incrRef();
+ declareAsNew();
+}
+
+void MEDCouplingRMesh::getBoundingBox(double *bbox) const
+{
+ //not implemented yet !
+}
+
+MEDCouplingFieldDouble *MEDCouplingRMesh::getMeasureField() const
+{
+ //not implemented yet !
+}
--- /dev/null
+// Copyright (C) 2007-2008 CEA/DEN, EDF R&D
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+#ifndef __PARAMEDMEM_MEDCOUPLINGRMESH_HXX__
+#define __PARAMEDMEM_MEDCOUPLINGRMESH_HXX__
+
+#include "MEDCoupling.hxx"
+#include "MEDCouplingMesh.hxx"
+
+namespace ParaMEDMEM
+{
+ class DataArrayDouble;
+
+ class MEDCouplingRMesh : public MEDCouplingMesh
+ {
+ public:
+ static MEDCouplingRMesh *New();
+ void updateTime();
+ bool isEqual(const MEDCouplingMesh *other, double prec) const;
+ void checkCoherency() const throw(INTERP_KERNEL::Exception);
+ bool isStructured() const;
+ int getNumberOfCells() const;
+ int getNumberOfNodes() const;
+ int getSpaceDimension() const;
+ int getMeshDimension() const;
+ DataArrayDouble *getCoordsAt(int i) const throw(INTERP_KERNEL::Exception);
+ void setCoords(DataArrayDouble *coordsX,
+ DataArrayDouble *coordsY=0,
+ DataArrayDouble *coordsZ=0);
+ // tools
+ void getBoundingBox(double *bbox) const;
+ MEDCouplingFieldDouble *getMeasureField() const;
+ private:
+ MEDCouplingRMesh();
+ ~MEDCouplingRMesh();
+ private:
+ DataArrayDouble *_x_array;
+ DataArrayDouble *_y_array;
+ DataArrayDouble *_z_array;
+ };
+}
+
+#endif
+++ /dev/null
-// Copyright (C) 2007-2008 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 "MEDCouplingSMesh.hxx"
-#include "MemArray.hxx"
-
-using namespace ParaMEDMEM;
-
-MEDCouplingSMesh::MEDCouplingSMesh():_x_array(0),_y_array(0),_z_array(0)
-{
-}
-
-MEDCouplingSMesh::~MEDCouplingSMesh()
-{
- if(_x_array)
- _x_array->decrRef();
- if(_y_array)
- _y_array->decrRef();
- if(_z_array)
- _z_array->decrRef();
-}
-
-MEDCouplingSMesh *MEDCouplingSMesh::New()
-{
- return new MEDCouplingSMesh;
-}
-
-void MEDCouplingSMesh::updateTime()
-{
- if(_x_array)
- updateTimeWith(*_x_array);
- if(_y_array)
- updateTimeWith(*_y_array);
- if(_z_array)
- updateTimeWith(*_z_array);
-}
-
-bool MEDCouplingSMesh::isEqual(const MEDCouplingMesh *other, double prec) const
-{
- const MEDCouplingSMesh *otherC=dynamic_cast<const MEDCouplingSMesh *>(other);
- if(!otherC)
- return false;
- return true;
-}
-
-void MEDCouplingSMesh::checkCoherency() const throw(INTERP_KERNEL::Exception)
-{
- const char msg0[]="Invalid ";
- const char msg1[]=" array ! Must contain more than 1 element.";
- if(_x_array)
- if(_x_array->getNbOfElems()<2)
- {
- std::ostringstream os; os << msg0 << 'X' << msg1;
- throw INTERP_KERNEL::Exception(os.str().c_str());
- }
- if(_y_array)
- if(_y_array->getNbOfElems()<2)
- {
- std::ostringstream os; os << msg0 << 'Y' << msg1;
- throw INTERP_KERNEL::Exception(os.str().c_str());
- }
- if(_z_array)
- if(_z_array->getNbOfElems()<2)
- {
- std::ostringstream os; os << msg0 << 'Z' << msg1;
- throw INTERP_KERNEL::Exception(os.str().c_str());
- }
-}
-
-bool MEDCouplingSMesh::isStructured() const
-{
- return true;
-}
-
-int MEDCouplingSMesh::getNumberOfCells() const
-{
- int ret=1;
- if(_x_array)
- ret*=_x_array->getNbOfElems()-1;
- if(_y_array)
- ret*=_y_array->getNbOfElems()-1;
- if(_z_array)
- ret*=_z_array->getNbOfElems()-1;
- return ret;
-}
-
-int MEDCouplingSMesh::getNumberOfNodes() const
-{
- int ret=1;
- if(_x_array)
- ret*=_x_array->getNbOfElems();
- if(_y_array)
- ret*=_y_array->getNbOfElems();
- if(_z_array)
- ret*=_z_array->getNbOfElems();
- return ret;
-}
-
-int MEDCouplingSMesh::getSpaceDimension() const
-{
- int ret=0;
- if(_x_array)
- ret++;
- if(_y_array)
- ret++;
- if(_z_array)
- ret++;
- return ret;
-}
-
-int MEDCouplingSMesh::getMeshDimension() const
-{
- int ret=0;
- if(_x_array)
- ret++;
- if(_y_array)
- ret++;
- if(_z_array)
- ret++;
- return ret;
-}
-
-DataArrayDouble *MEDCouplingSMesh::getCoordsAt(int i) const throw(INTERP_KERNEL::Exception)
-{
- switch(i)
- {
- case 0:
- return _x_array;
- case 1:
- return _y_array;
- case 2:
- return _z_array;
- default:
- throw INTERP_KERNEL::Exception("Invalid rank specified must be 0 or 1 or 2.");
- }
-}
-
-void MEDCouplingSMesh::setCoords(DataArrayDouble *coordsX, DataArrayDouble *coordsY, DataArrayDouble *coordsZ)
-{
- if(_x_array)
- _x_array->decrRef();
- _x_array=coordsX;
- if(_x_array)
- _x_array->incrRef();
- if(_y_array)
- _y_array->decrRef();
- _y_array=coordsY;
- if(_y_array)
- _y_array->incrRef();
- if(_z_array)
- _z_array->decrRef();
- _z_array=coordsZ;
- if(_z_array)
- _z_array->incrRef();
- declareAsNew();
-}
-
+++ /dev/null
-// Copyright (C) 2007-2008 CEA/DEN, EDF R&D
-//
-// This library is free software; you can redistribute it and/or
-// modify it under the terms of the GNU Lesser General Public
-// License as published by the Free Software Foundation; either
-// version 2.1 of the License.
-//
-// This library is distributed in the hope that it will be useful,
-// but WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-// Lesser General Public License for more details.
-//
-// You should have received a copy of the GNU Lesser General Public
-// License along with this library; if not, write to the Free Software
-// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
-//
-// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
-//
-#ifndef __PARAMEDMEM_MEDCOUPLINGSMESH_HXX__
-#define __PARAMEDMEM_MEDCOUPLINGSMESH_HXX__
-
-#include "MEDCouplingMesh.hxx"
-
-namespace ParaMEDMEM
-{
- class DataArrayDouble;
-
- class MEDCouplingSMesh : public MEDCouplingMesh
- {
- public:
- static MEDCouplingSMesh *New();
- void updateTime();
- bool isEqual(const MEDCouplingMesh *other, double prec) const;
- void checkCoherency() const throw(INTERP_KERNEL::Exception);
- bool isStructured() const;
- int getNumberOfCells() const;
- int getNumberOfNodes() const;
- int getSpaceDimension() const;
- int getMeshDimension() const;
- DataArrayDouble *getCoordsAt(int i) const throw(INTERP_KERNEL::Exception);
- void setCoords(DataArrayDouble *coordsX,
- DataArrayDouble *coordsY=0,
- DataArrayDouble *coordsZ=0);
- private:
- MEDCouplingSMesh();
- ~MEDCouplingSMesh();
- private:
- DataArrayDouble *_x_array;
- DataArrayDouble *_y_array;
- DataArrayDouble *_z_array;
- };
-}
-
-#endif
--- /dev/null
+// Copyright (C) 2007-2008 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 "MEDCouplingTimeDiscretization.hxx"
+#include "MemArray.hxx"
+
+#include <cmath>
+
+using namespace ParaMEDMEM;
+
+const char MEDCouplingNoTimeLabel::EXCEPTION_MSG[]="MEDCouplingNoTimeLabel::setTime : no time info attached.";
+
+const char MEDCouplingWithTimeStep::EXCEPTION_MSG[]="No data on this time.";
+
+const double MEDCouplingTimeDiscretization::TIME_TOLERANCE_DFT=1.e-12;
+
+MEDCouplingTimeDiscretization *MEDCouplingTimeDiscretization::New(TypeOfTimeDiscretization type)
+{
+ switch(type)
+ {
+ case MEDCouplingNoTimeLabel::DISCRETIZATION:
+ return new MEDCouplingNoTimeLabel;
+ case MEDCouplingWithTimeStep::DISCRETIZATION:
+ return new MEDCouplingWithTimeStep;
+ default:
+ throw INTERP_KERNEL::Exception("Time discretization not implemented yet");
+ }
+}
+
+MEDCouplingTimeDiscretization::MEDCouplingTimeDiscretization():_time_tolerance(TIME_TOLERANCE_DFT),_array(0)
+{
+}
+
+MEDCouplingTimeDiscretization::MEDCouplingTimeDiscretization(const MEDCouplingTimeDiscretization& other, bool deepCpy):_time_tolerance(other._time_tolerance)
+{
+ if(other._array)
+ _array=other._array->performCpy(deepCpy);
+ else
+ _array=0;
+}
+
+MEDCouplingTimeDiscretization::~MEDCouplingTimeDiscretization()
+{
+ if(_array)
+ _array->decrRef();
+}
+
+void MEDCouplingTimeDiscretization::setArray(DataArrayDouble *array, TimeLabel *owner)
+{
+ if(array!=_array)
+ {
+ if(_array)
+ _array->decrRef();
+ _array=array;
+ if(_array)
+ _array->incrRef();
+ owner->declareAsNew();
+ }
+}
+
+void MEDCouplingTimeDiscretization::setArrays(const std::vector<DataArrayDouble *>& arrays, TimeLabel *owner) throw(INTERP_KERNEL::Exception)
+{
+ if(arrays.size()!=1)
+ throw INTERP_KERNEL::Exception("MEDCouplingTimeDiscretization::setArrays : number of arrays must be one.");
+ setArray(arrays.back(),owner);
+}
+
+void MEDCouplingTimeDiscretization::getArrays(std::vector<DataArrayDouble *>& arrays) const
+{
+ arrays.resize(1);
+ arrays[0]=_array;
+}
+
+bool MEDCouplingTimeDiscretization::isBefore(const MEDCouplingTimeDiscretization *other) const throw(INTERP_KERNEL::Exception)
+{
+ int dt,it;
+ double time1=getEndTime(dt,it)-_time_tolerance;
+ double time2=other->getStartTime(dt,it)+other->getTimeTolerance();
+ return time1<=time2;
+}
+
+bool MEDCouplingTimeDiscretization::isStrictlyBefore(const MEDCouplingTimeDiscretization *other) const throw(INTERP_KERNEL::Exception)
+{
+ int dt,it;
+ double time1=getEndTime(dt,it)+_time_tolerance;
+ double time2=other->getStartTime(dt,it)-other->getTimeTolerance();
+ return time1<time2;
+}
+
+MEDCouplingNoTimeLabel::MEDCouplingNoTimeLabel()
+{
+}
+
+MEDCouplingNoTimeLabel::MEDCouplingNoTimeLabel(const MEDCouplingTimeDiscretization& other, bool deepCpy):MEDCouplingTimeDiscretization(other,deepCpy)
+{
+}
+
+MEDCouplingTimeDiscretization *MEDCouplingNoTimeLabel::performCpy(bool deepCpy) const
+{
+ return new MEDCouplingNoTimeLabel(*this,deepCpy);
+}
+
+void MEDCouplingNoTimeLabel::checkTimePresence(double time) const throw(INTERP_KERNEL::Exception)
+{
+ throw INTERP_KERNEL::Exception(EXCEPTION_MSG);
+}
+
+DataArrayDouble *MEDCouplingNoTimeLabel::getArrayOnTime(double time) const throw(INTERP_KERNEL::Exception)
+{
+ throw INTERP_KERNEL::Exception(EXCEPTION_MSG);
+}
+
+bool MEDCouplingNoTimeLabel::isBefore(const MEDCouplingTimeDiscretization *other) const throw(INTERP_KERNEL::Exception)
+{
+ throw INTERP_KERNEL::Exception(EXCEPTION_MSG);
+}
+
+bool MEDCouplingNoTimeLabel::isStrictlyBefore(const MEDCouplingTimeDiscretization *other) const throw(INTERP_KERNEL::Exception)
+{
+ throw INTERP_KERNEL::Exception(EXCEPTION_MSG);
+}
+
+double MEDCouplingNoTimeLabel::getStartTime(int& dt, int& it) const throw(INTERP_KERNEL::Exception)
+{
+ throw INTERP_KERNEL::Exception(EXCEPTION_MSG);
+}
+
+double MEDCouplingNoTimeLabel::getEndTime(int& dt, int& it) const throw(INTERP_KERNEL::Exception)
+{
+ throw INTERP_KERNEL::Exception(EXCEPTION_MSG);
+}
+
+void MEDCouplingNoTimeLabel::setStartTime(double time, int dt, int it) throw(INTERP_KERNEL::Exception)
+{
+ throw INTERP_KERNEL::Exception(EXCEPTION_MSG);
+}
+
+void MEDCouplingNoTimeLabel::setEndTime(double time, int dt, int it) throw(INTERP_KERNEL::Exception)
+{
+ throw INTERP_KERNEL::Exception(EXCEPTION_MSG);
+}
+
+void MEDCouplingNoTimeLabel::getValueOnTime(int eltId, double time, double *value) const throw(INTERP_KERNEL::Exception)
+{
+ throw INTERP_KERNEL::Exception(EXCEPTION_MSG);
+}
+
+void MEDCouplingNoTimeLabel::getValueOnDiscTime(int eltId, int dt, int it, double *value) const throw(INTERP_KERNEL::Exception)
+{
+ throw INTERP_KERNEL::Exception(EXCEPTION_MSG);
+}
+
+MEDCouplingWithTimeStep::MEDCouplingWithTimeStep(const MEDCouplingWithTimeStep& other, bool deepCpy):MEDCouplingTimeDiscretization(other,deepCpy),
+ _time(other._time),_dt(other._dt),_it(other._it)
+{
+}
+
+MEDCouplingWithTimeStep::MEDCouplingWithTimeStep():_time(0.),_dt(-1),_it(-1)
+{
+}
+
+MEDCouplingTimeDiscretization *MEDCouplingWithTimeStep::performCpy(bool deepCpy) const
+{
+ return new MEDCouplingWithTimeStep(*this,deepCpy);
+}
+
+void MEDCouplingWithTimeStep::checkNoTimePresence() const throw(INTERP_KERNEL::Exception)
+{
+ throw INTERP_KERNEL::Exception("No time specified on a field defined on one time");
+}
+
+void MEDCouplingWithTimeStep::checkTimePresence(double time) const throw(INTERP_KERNEL::Exception)
+{
+ if(std::fabs(time-_time)>_time_tolerance)
+ {
+ std::ostringstream stream;
+ stream << "The field is defined on time " << _time << " with eps=" << _time_tolerance << " and asking time = " << time << " !";
+ throw INTERP_KERNEL::Exception(stream.str().c_str());
+ }
+}
+
+DataArrayDouble *MEDCouplingWithTimeStep::getArrayOnTime(double time) const throw(INTERP_KERNEL::Exception)
+{
+ if(std::fabs(time-_time)<=_time_tolerance)
+ {
+ if(_array)
+ _array->incrRef();
+ return _array;
+ }
+ else
+ throw INTERP_KERNEL::Exception(EXCEPTION_MSG);
+}
+
+void MEDCouplingWithTimeStep::getValueOnTime(int eltId, double time, double *value) const throw(INTERP_KERNEL::Exception)
+{
+ if(std::fabs(time-_time)<=_time_tolerance)
+ if(_array)
+ _array->getTuple(eltId,value);
+ else
+ throw INTERP_KERNEL::Exception("No array existing.");
+ else
+ throw INTERP_KERNEL::Exception(EXCEPTION_MSG);
+}
+
+void MEDCouplingWithTimeStep::getValueOnDiscTime(int eltId, int dt, int it, double *value) const throw(INTERP_KERNEL::Exception)
+{
+ if(_dt==dt && _it==it)
+ if(_array)
+ _array->getTuple(eltId,value);
+ else
+ throw INTERP_KERNEL::Exception("No array existing.");
+ else
+ throw INTERP_KERNEL::Exception("No data on this discrete time.");
+}
+
+MEDCouplingTwoTimeSteps::MEDCouplingTwoTimeSteps():_start_time(0.),_end_time(0.),_start_dt(-1),_end_dt(-1),_start_it(-1),_end_it(-1),_end_array(0)
+{
+}
+
+MEDCouplingTwoTimeSteps::~MEDCouplingTwoTimeSteps()
+{
+ if(_end_array)
+ _end_array->decrRef();
+}
+
+void MEDCouplingTwoTimeSteps::checkNoTimePresence() const throw(INTERP_KERNEL::Exception)
+{
+ throw INTERP_KERNEL::Exception("The field presents a time to be specified in every access !");
+}
+
+void MEDCouplingTwoTimeSteps::checkTimePresence(double time) const throw(INTERP_KERNEL::Exception)
+{
+ if(time<_start_time-_time_tolerance || time>_end_time+_time_tolerance)
+ {
+ std::ostringstream stream;
+ stream << "The field is defined between times " << _start_time << " and " << _end_time << " with tolerance ";
+ stream << _time_tolerance << " and trying to access on time = " << time;
+ throw INTERP_KERNEL::Exception(stream.str().c_str());
+ }
+}
+
+void MEDCouplingTwoTimeSteps::getArrays(std::vector<DataArrayDouble *>& arrays) const
+{
+ arrays.resize(2);
+ arrays[0]=_array;
+ arrays[1]=_end_array;
+}
--- /dev/null
+// Copyright (C) 2007-2008 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 __MEDCOUPLINGTIMEDISCRETIZATION_HXX__
+#define __MEDCOUPLINGTIMEDISCRETIZATION_HXX__
+
+#include "MEDCoupling.hxx"
+#include "RefCountObject.hxx"
+#include "InterpKernelException.hxx"
+
+#include <vector>
+
+namespace ParaMEDMEM
+{
+ class DataArrayDouble;
+ class TimeLabel;
+
+ class MEDCOUPLING_EXPORT MEDCouplingTimeDiscretization
+ {
+ protected:
+ MEDCouplingTimeDiscretization();
+ MEDCouplingTimeDiscretization(const MEDCouplingTimeDiscretization& other, bool deepCpy);
+ public:
+ static MEDCouplingTimeDiscretization *New(TypeOfTimeDiscretization type);
+ virtual MEDCouplingTimeDiscretization *performCpy(bool deepCpy) const = 0;
+ void setTimeTolerance(double val);
+ double getTimeTolerance() const { return _time_tolerance; }
+ virtual void checkNoTimePresence() const throw(INTERP_KERNEL::Exception) = 0;
+ virtual void checkTimePresence(double time) const throw(INTERP_KERNEL::Exception) = 0;
+ virtual void setArray(DataArrayDouble *array, TimeLabel *owner);
+ virtual void setArrays(const std::vector<DataArrayDouble *>& arrays, TimeLabel *owner) throw(INTERP_KERNEL::Exception);
+ DataArrayDouble *getArray() const { return _array; }
+ virtual DataArrayDouble *getEndArray() const { return _array; }
+ //! Warning contrary to getArray method this method returns an object to deal with.
+ virtual DataArrayDouble *getArrayOnTime(double time) const throw(INTERP_KERNEL::Exception) = 0;
+ virtual void getArrays(std::vector<DataArrayDouble *>& arrays) const;
+ virtual bool isBefore(const MEDCouplingTimeDiscretization *other) const throw(INTERP_KERNEL::Exception);
+ virtual bool isStrictlyBefore(const MEDCouplingTimeDiscretization *other) const throw(INTERP_KERNEL::Exception);
+ double getTime(int& dt, int& it) const throw(INTERP_KERNEL::Exception) { return getStartTime(dt,it); }
+ virtual double getStartTime(int& dt, int& it) const throw(INTERP_KERNEL::Exception) = 0;
+ virtual double getEndTime(int& dt, int& it) const throw(INTERP_KERNEL::Exception) = 0;
+ void setTime(double time, int dt, int it) throw(INTERP_KERNEL::Exception) { setStartTime(time,dt,it); }
+ virtual void setStartTime(double time, int dt, int it) throw(INTERP_KERNEL::Exception) = 0;
+ virtual void setEndTime(double time, int dt, int it) throw(INTERP_KERNEL::Exception) = 0;
+ virtual void getValueOnTime(int eltId, double time, double *value) const throw(INTERP_KERNEL::Exception) = 0;
+ virtual void getValueOnDiscTime(int eltId, int dt, int it, double *value) const throw(INTERP_KERNEL::Exception) = 0;
+ virtual ~MEDCouplingTimeDiscretization();
+ protected:
+ double _time_tolerance;
+ DataArrayDouble *_array;
+ protected:
+ static const double TIME_TOLERANCE_DFT;
+ };
+
+ class MEDCOUPLING_EXPORT MEDCouplingNoTimeLabel : public MEDCouplingTimeDiscretization
+ {
+ public:
+ MEDCouplingNoTimeLabel();
+ MEDCouplingNoTimeLabel(const MEDCouplingTimeDiscretization& other, bool deepCpy);
+ MEDCouplingTimeDiscretization *performCpy(bool deepCpy) const;
+ void checkNoTimePresence() const throw(INTERP_KERNEL::Exception) { }
+ void checkTimePresence(double time) const throw(INTERP_KERNEL::Exception);
+ DataArrayDouble *getArrayOnTime(double time) const throw(INTERP_KERNEL::Exception);
+ bool isBefore(const MEDCouplingTimeDiscretization *other) const throw(INTERP_KERNEL::Exception);
+ bool isStrictlyBefore(const MEDCouplingTimeDiscretization *other) const throw(INTERP_KERNEL::Exception);
+ double getStartTime(int& dt, int& it) const throw(INTERP_KERNEL::Exception);
+ double getEndTime(int& dt, int& it) const throw(INTERP_KERNEL::Exception);
+ void setStartTime(double time, int dt, int it) throw(INTERP_KERNEL::Exception);
+ void setEndTime(double time, int dt, int it) throw(INTERP_KERNEL::Exception);
+ void getValueOnTime(int eltId, double time, double *value) const throw(INTERP_KERNEL::Exception);
+ void getValueOnDiscTime(int eltId, int dt, int it, double *value) const throw(INTERP_KERNEL::Exception);
+ public:
+ static const TypeOfTimeDiscretization DISCRETIZATION=NO_TIME;
+ private:
+ static const char EXCEPTION_MSG[];
+ };
+
+ class MEDCOUPLING_EXPORT MEDCouplingWithTimeStep : public MEDCouplingTimeDiscretization
+ {
+ protected:
+ MEDCouplingWithTimeStep(const MEDCouplingWithTimeStep& other, bool deepCpy);
+ public:
+ MEDCouplingWithTimeStep();
+ MEDCouplingTimeDiscretization *performCpy(bool deepCpy) const;
+ void checkNoTimePresence() const throw(INTERP_KERNEL::Exception);
+ void checkTimePresence(double time) const throw(INTERP_KERNEL::Exception);
+ void setStartTime(double time, int dt, int it) throw(INTERP_KERNEL::Exception) { _time=time; _dt=dt; _it=it; }
+ void setEndTime(double time, int dt, int it) throw(INTERP_KERNEL::Exception) { _time=time; _dt=dt; _it=it; }
+ double getStartTime(int& dt, int& it) const throw(INTERP_KERNEL::Exception) { dt=_dt; it=_it; return _time; }
+ double getEndTime(int& dt, int& it) const throw(INTERP_KERNEL::Exception) { dt=_dt; it=_it; return _time; }
+ DataArrayDouble *getArrayOnTime(double time) const throw(INTERP_KERNEL::Exception);
+ void getValueOnTime(int eltId, double time, double *value) const throw(INTERP_KERNEL::Exception);
+ void getValueOnDiscTime(int eltId, int dt, int it, double *value) const throw(INTERP_KERNEL::Exception);
+ public:
+ static const TypeOfTimeDiscretization DISCRETIZATION=ONE_TIME;
+ private:
+ static const char EXCEPTION_MSG[];
+ protected:
+ double _time;
+ int _dt;
+ int _it;
+ };
+
+ class MEDCOUPLING_EXPORT MEDCouplingTwoTimeSteps : public MEDCouplingTimeDiscretization
+ {
+ protected:
+ MEDCouplingTwoTimeSteps();
+ ~MEDCouplingTwoTimeSteps();
+ public:
+ void checkNoTimePresence() const throw(INTERP_KERNEL::Exception);
+ void checkTimePresence(double time) const throw(INTERP_KERNEL::Exception);
+ void getArrays(std::vector<DataArrayDouble *>& arrays) const;
+ DataArrayDouble *getEndArray() const { return _end_array; }
+ void setStartTime(double time, int dt, int it) throw(INTERP_KERNEL::Exception) { _start_time=time; _start_dt=dt; _start_it=it; }
+ void setEndTime(double time, int dt, int it) throw(INTERP_KERNEL::Exception) { _end_time=time; _end_dt=dt; _end_it=it; }
+ double getStartTime(int& dt, int& it) const throw(INTERP_KERNEL::Exception) { dt=_start_dt; it=_start_it; return _start_time; }
+ double getEndTime(int& dt, int& it) const throw(INTERP_KERNEL::Exception) { dt=_end_dt; it=_end_it; return _end_time; }
+ protected:
+ double _start_time;
+ double _end_time;
+ int _start_dt;
+ int _end_dt;
+ int _start_it;
+ int _end_it;
+ DataArrayDouble *_end_array;
+ };
+
+ class MEDCOUPLING_EXPORT MEDCouplingLinearTime : public MEDCouplingTwoTimeSteps
+ {
+ public:
+ static const TypeOfTimeDiscretization DISCRETIZATION=LINEAR_TIME;
+ };
+}
+
+#endif
// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
//
#include "MEDCouplingUMesh.hxx"
+#include "MEDCouplingFieldDouble.hxx"
#include "CellModel.hxx"
+#include "VolSurfFormulae.hxx"
#include <sstream>
void MEDCouplingUMesh::updateTime()
{
+ MEDCouplingPointSet::updateTime();
if(_nodal_connec)
{
updateTimeWith(*_nodal_connec);
{
updateTimeWith(*_nodal_connec_index);
}
- if(_coords)
- {
- updateTimeWith(*_coords);
- }
}
MEDCouplingUMesh::MEDCouplingUMesh():_iterator(-1),_mesh_dim(-1),
- _nodal_connec(0),_nodal_connec_index(0),_coords(0)
+ _nodal_connec(0),_nodal_connec_index(0)
{
}
declareAsNew();
}
-void MEDCouplingUMesh::setCoords(DataArrayDouble *coords)
-{
- if( coords != _coords )
- {
- if (_coords)
- _coords->decrRef();
- _coords=coords;
- if(_coords)
- _coords->incrRef();
- declareAsNew();
- }
-}
-
void MEDCouplingUMesh::insertNextCell(INTERP_KERNEL::NormalizedCellType type, int size, const int *nodalConnOfCell)
{
int *pt=_nodal_connec_index->getPointer();
void MEDCouplingUMesh::finishInsertingCells()
{
- int *pt=_nodal_connec_index->getPointer();
+ const int *pt=_nodal_connec_index->getConstPointer();
int idx=pt[_iterator];
_nodal_connec->reAlloc(idx);
return false;
if(_types!=otherC->_types)
return false;
- if(!_coords->isEqual(otherC->_coords,prec))
+ if(!areCoordsEqual(*otherC,prec))
return false;
return true;
}
int *revNodalIndxPtr=new int[nbOfNodes+1];
revNodalIndx->useArray(revNodalIndxPtr,true,CPP_DEALLOC,nbOfNodes+1,1);
std::fill(revNodalIndxPtr,revNodalIndxPtr+nbOfNodes+1,0);
- const int *conn=_nodal_connec->getPointer();
- const int *connIndex=_nodal_connec_index->getPointer();
+ const int *conn=_nodal_connec->getConstPointer();
+ const int *connIndex=_nodal_connec_index->getConstPointer();
int nbOfCells=getNumberOfCells();
int nbOfEltsInRevNodal=0;
for(int eltId=0;eltId<nbOfCells;eltId++)
std::fill(traducer,traducer+nbOfNodes,-1);
ret->useArray(traducer,true,CPP_DEALLOC,nbOfNodes,1);
int nbOfCells=getNumberOfCells();
- const int *connIndex=_nodal_connec_index->getPointer();
+ const int *connIndex=_nodal_connec_index->getConstPointer();
int *conn=_nodal_connec->getPointer();
for(int i=0;i<nbOfCells;i++)
for(int j=connIndex[i]+1;j<connIndex[i+1];j++)
conn[j]=traducer[conn[j]];
DataArrayDouble *newCoords=DataArrayDouble::New();
double *newCoordsPtr=new double[newNbOfNodes*spaceDim];
- const double *oldCoordsPtr=_coords->getPointer();
+ const double *oldCoordsPtr=_coords->getConstPointer();
newCoords->useArray(newCoordsPtr,true,CPP_DEALLOC,newNbOfNodes,spaceDim);
int *work=std::find_if(traducer,traducer+nbOfNodes,std::bind2nd(std::not_equal_to<int>(),-1));
for(;work!=traducer+nbOfNodes;work=std::find_if(work,traducer+nbOfNodes,std::bind2nd(std::not_equal_to<int>(),-1)))
return ret;
}
-MEDCouplingUMesh *MEDCouplingUMesh::buildPartOfMySelf(const int *start, const int *end, bool keepCoords) const
+MEDCouplingPointSet *MEDCouplingUMesh::buildPartOfMySelf(const int *start, const int *end, bool keepCoords) const
{
MEDCouplingUMesh *ret=buildPartOfMySelfKeepCoords(start,end);
if(!keepCoords)
return ret;
}
+/*!
+ * Given a boundary box 'bbox' returns elements 'elems' contained in this 'bbox'.
+ * Warning 'elems' is incremented during the call so if elems is not empty before call returned elements will be
+ * added in 'elems' parameter.
+ */
+void MEDCouplingUMesh::giveElemsInBoundingBox(const double *bbox, double eps, std::vector<int>& elems)
+{
+ int dim=getSpaceDimension();
+ double* elem_bb=new double[2*dim];
+ const int* conn = getNodalConnectivity()->getConstPointer();
+ const int* conn_index= getNodalConnectivityIndex()->getConstPointer();
+ const double* coords = getCoords()->getConstPointer();
+ int nbOfCells=getNumberOfCells();
+ for ( int ielem=0; ielem<nbOfCells;ielem++ )
+ {
+ for (int i=0; i<dim; i++)
+ {
+ elem_bb[i*2]=std::numeric_limits<double>::max();
+ elem_bb[i*2+1]=-std::numeric_limits<double>::max();
+ }
+
+ for (int inode=conn_index[ielem]+1; inode<conn_index[ielem+1]; inode++)//+1 due to offset of cell type.
+ {
+ int node= conn[inode];
+
+ for (int idim=0; idim<dim; idim++)
+ {
+ if ( coords[node*dim+idim] < elem_bb[idim*2] )
+ {
+ elem_bb[idim*2] = coords[node*dim+idim] ;
+ }
+ if ( coords[node*dim+idim] > elem_bb[idim*2+1] )
+ {
+ elem_bb[idim*2+1] = coords[node*dim+idim] ;
+ }
+ }
+ }
+ if (intersectsBoundingBox(elem_bb, bbox, dim, eps))
+ {
+ elems.push_back(ielem);
+ }
+ }
+ delete [] elem_bb;
+}
+
INTERP_KERNEL::NormalizedCellType MEDCouplingUMesh::getTypeOfCell(int cellId) const
{
- int *ptI=_nodal_connec_index->getPointer();
- int *pt=_nodal_connec->getPointer();
+ const int *ptI=_nodal_connec_index->getConstPointer();
+ const int *pt=_nodal_connec->getConstPointer();
return (INTERP_KERNEL::NormalizedCellType) pt[ptI[cellId]];
}
void MEDCouplingUMesh::setConnectivity(DataArrayInt *conn, DataArrayInt *connIndex, bool isComputingTypes)
{
- if(_nodal_connec!=conn)
- {
- if(_nodal_connec)
- _nodal_connec->decrRef();
- _nodal_connec=conn;
- if(_nodal_connec)
- _nodal_connec->incrRef();
- }
- if(_nodal_connec_index!=connIndex)
- {
- if(_nodal_connec_index)
- _nodal_connec_index->decrRef();
- _nodal_connec_index=connIndex;
- if(_nodal_connec_index)
- _nodal_connec_index->incrRef();
- }
+ DataArrayInt::setArrayIn(conn,_nodal_connec);
+ DataArrayInt::setArrayIn(connIndex,_nodal_connec_index);
if(isComputingTypes)
computeTypes();
}
-MEDCouplingUMesh::MEDCouplingUMesh(const MEDCouplingUMesh& other, bool deepCpy):MEDCouplingMesh(other),_iterator(-1),_mesh_dim(other._mesh_dim),
- _nodal_connec(0),_nodal_connec_index(0),_coords(0),
+MEDCouplingUMesh::MEDCouplingUMesh(const MEDCouplingUMesh& other, bool deepCpy):MEDCouplingPointSet(other,deepCpy),_iterator(-1),_mesh_dim(other._mesh_dim),
+ _nodal_connec(0),_nodal_connec_index(0),
_types(other._types)
{
if(other._nodal_connec)
_nodal_connec=other._nodal_connec->performCpy(deepCpy);
if(other._nodal_connec_index)
_nodal_connec_index=other._nodal_connec_index->performCpy(deepCpy);
- if(other._coords)
- _coords=other._coords->performCpy(deepCpy);
}
MEDCouplingUMesh::~MEDCouplingUMesh()
_nodal_connec->decrRef();
if(_nodal_connec_index)
_nodal_connec_index->decrRef();
- if(_coords)
- _coords->decrRef();
}
void MEDCouplingUMesh::computeTypes()
if(_nodal_connec && _nodal_connec_index)
{
_types.clear();
- const int *conn=_nodal_connec->getPointer();
- const int *connIndex=_nodal_connec_index->getPointer();
+ const int *conn=_nodal_connec->getConstPointer();
+ const int *connIndex=_nodal_connec_index->getConstPointer();
int nbOfElem=_nodal_connec_index->getNbOfElems()-1;
for(const int *pt=connIndex;pt!=connIndex+nbOfElem;pt++)
_types.insert((INTERP_KERNEL::NormalizedCellType)conn[*pt]);
throw INTERP_KERNEL::Exception("Reverse nodal connectivity computation requires full connectivity and coordinates set in unstructured mesh.");
}
-bool MEDCouplingUMesh::isStructured() const
-{
- return false;
-}
-
int MEDCouplingUMesh::getNumberOfCells() const
{
if(_nodal_connec_index)
else
return _iterator;
else
- throw INTERP_KERNEL::Exception("Unable to get number of cells because no coordinates specified !");
+ throw INTERP_KERNEL::Exception("Unable to get number of cells because no connectivity specified !");
}
-int MEDCouplingUMesh::getNumberOfNodes() const
+int MEDCouplingUMesh::getMeshLength() const
{
- if(_coords)
- return _coords->getNumberOfTuples();
- else
- throw INTERP_KERNEL::Exception("Unable to get number of nodes because no coordinates specified !");
+ return _nodal_connec->getNbOfElems();
}
-int MEDCouplingUMesh::getSpaceDimension() const
+void MEDCouplingUMesh::getTinySerializationInformation(std::vector<int>& tinyInfo) const
{
- if(_coords)
- return _coords->getNumberOfComponents();
- else
- throw INTERP_KERNEL::Exception("Unable to get space dimension because no coordinates specified !");
+ tinyInfo.resize(5);
+ tinyInfo[0] = getSpaceDimension();
+ tinyInfo[1] = getMeshDimension();
+ tinyInfo[2] = getNumberOfNodes();
+ tinyInfo[3] = getNumberOfCells();
+ tinyInfo[4] = getMeshLength();
}
-int MEDCouplingUMesh::getMeshLength() const
+/*!
+ * @param tinyInfo must be equal to the result given by getTinySerializationInformation method.
+ */
+void MEDCouplingUMesh::resizeForSerialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2)
{
- return _nodal_connec->getNbOfElems();
+ a1->alloc(tinyInfo[4]+tinyInfo[3]+1,1);
+ a2->alloc(tinyInfo[2],tinyInfo[0]);
+}
+
+void MEDCouplingUMesh::serialize(DataArrayInt *&a1, DataArrayDouble *&a2)
+{
+ a1=DataArrayInt::New();
+ a1->alloc(getMeshLength()+getNumberOfCells()+1,1);
+ int *ptA1=a1->getPointer();
+ const int *conn=getNodalConnectivity()->getConstPointer();
+ const int *index=getNodalConnectivityIndex()->getConstPointer();
+ ptA1=std::copy(index,index+getNumberOfCells()+1,ptA1);
+ std::copy(conn,conn+getMeshLength(),ptA1);
+ a2=getCoords();
+ a2->incrRef();
+}
+
+/*!
+ * @param tinyInfo must be equal to the result given by getTinySerializationInformation method.
+ */
+MEDCouplingPointSet *MEDCouplingUMesh::buildObjectFromUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2)
+{
+ MEDCouplingUMesh* meshing = MEDCouplingUMesh::New() ;
+ // Coordinates
+ meshing->setCoords(a2);
+ // Connectivity
+ const int *recvBuffer=a1->getConstPointer();
+ DataArrayInt* myConnecIndex=DataArrayInt::New();
+ myConnecIndex->alloc(tinyInfo[3]+1,1);
+ std::copy(recvBuffer,recvBuffer+tinyInfo[3]+1,myConnecIndex->getPointer());
+ DataArrayInt* myConnec=DataArrayInt::New();
+ myConnec->alloc(tinyInfo[4],1);
+ std::copy(recvBuffer+tinyInfo[3]+1,recvBuffer+tinyInfo[3]+1+tinyInfo[4],myConnec->getPointer());
+ meshing->setConnectivity(myConnec, myConnecIndex) ;
+ myConnec->decrRef();
+ myConnecIndex->decrRef();
+ //
+ meshing->setMeshDimension(tinyInfo[1]);
+ return meshing;
}
MEDCouplingUMesh *MEDCouplingUMesh::buildPartOfMySelfKeepCoords(const int *start, const int *end) const
int nbOfElemsRet=end-start;
int *connIndexRet=new int[nbOfElemsRet+1];
connIndexRet[0]=0;
- const int *conn=_nodal_connec->getPointer();
- const int *connIndex=_nodal_connec_index->getPointer();
+ const int *conn=_nodal_connec->getConstPointer();
+ const int *connIndex=_nodal_connec_index->getConstPointer();
int newNbring=0;
for(const int *work=start;work!=end;work++,newNbring++)
connIndexRet[newNbring+1]=connIndexRet[newNbring]+connIndex[*work+1]-connIndex[*work];
return ret;
}
+/*!
+ * brief returns the volumes of the cells underlying the field \a field
+ *
+ * For 2D geometries, the returned field contains the areas.
+ * For 3D geometries, the returned field contains the volumes.
+ *
+ * param field field on which cells the volumes are required
+ * return field containing the volumes, area or length depending the meshdimension.
+ */
+MEDCouplingFieldDouble *MEDCouplingUMesh::getMeasureField() const
+{
+ int ipt, type ;
+ int nbelem = getNumberOfCells() ;
+ int dim_mesh = getMeshDimension();
+ int dim_space = getSpaceDimension() ;
+ const double *coords = getCoords()->getConstPointer() ;
+ const int *connec = getNodalConnectivity()->getConstPointer() ;
+ const int *connec_index = getNodalConnectivityIndex()->getConstPointer() ;
+
+
+ MEDCouplingFieldDouble* field = MEDCouplingFieldDouble::New(ON_CELLS);
+ DataArrayDouble* array = DataArrayDouble::New() ;
+ array->alloc(nbelem, 1) ;
+ double *area_vol = array->getPointer() ;
+
+ switch (dim_mesh)
+ {
+ case 2: // getting the areas
+ for ( int iel=0 ; iel<nbelem ; iel++ )
+ {
+ ipt = connec_index[iel] ;
+ type = connec[ipt] ;
+
+ switch ( type )
+ {
+ case INTERP_KERNEL::NORM_TRI3 :
+ case INTERP_KERNEL::NORM_TRI6 :
+ {
+ int N1 = connec[ipt+1];
+ int N2 = connec[ipt+2];
+ int N3 = connec[ipt+3];
+
+ area_vol[iel]=INTERP_KERNEL::calculateAreaForTria(coords+(dim_space*N1),
+ coords+(dim_space*N2),
+ coords+(dim_space*N3),
+ dim_space);
+ }
+ break ;
+
+ case INTERP_KERNEL::NORM_QUAD4 :
+ case INTERP_KERNEL::NORM_QUAD8 :
+ {
+ int N1 = connec[ipt+1];
+ int N2 = connec[ipt+2];
+ int N3 = connec[ipt+3];
+ int N4 = connec[ipt+4];
+
+ area_vol[iel]=INTERP_KERNEL::calculateAreaForQuad(coords+dim_space*N1,
+ coords+dim_space*N2,
+ coords+dim_space*N3,
+ coords+dim_space*N4,
+ dim_space) ;
+ }
+ break ;
+
+ case INTERP_KERNEL::NORM_POLYGON :
+ {
+ // We must remember that the first item is the type. That's
+ // why we substract 1 to get the number of nodes of this polygon
+ int size = connec_index[iel+1] - connec_index[iel] - 1 ;
+
+ const double **pts = new const double *[size] ;
+
+ for ( int inod=0 ; inod<size ; inod++ )
+ {
+ // Remember the first item is the type
+ pts[inod] = coords+dim_space*connec[ipt+inod+1] ;
+ }
+
+ area_vol[iel]=INTERP_KERNEL::calculateAreaForPolyg(pts,size,dim_space);
+ delete [] pts;
+ }
+ break ;
+
+ default :
+ throw INTERP_KERNEL::Exception("Bad Support to get Areas on it !");
+
+ } // End of switch
+
+ } // End of the loop over the cells
+ break;
+ case 3: // getting the volumes
+ for ( int iel=0 ; iel<nbelem ; iel++ )
+ {
+ ipt = connec_index[iel] ;
+ type = connec[ipt] ;
+
+ switch ( type )
+ {
+ case INTERP_KERNEL::NORM_TETRA4 :
+ case INTERP_KERNEL::NORM_TETRA10 :
+ {
+ int N1 = connec[ipt+1];
+ int N2 = connec[ipt+2];
+ int N3 = connec[ipt+3];
+ int N4 = connec[ipt+4];
+
+ area_vol[iel]=INTERP_KERNEL::calculateVolumeForTetra(coords+dim_space*N1,
+ coords+dim_space*N2,
+ coords+dim_space*N3,
+ coords+dim_space*N4) ;
+ }
+ break ;
+
+ case INTERP_KERNEL::NORM_PYRA5 :
+ case INTERP_KERNEL::NORM_PYRA13 :
+ {
+ int N1 = connec[ipt+1];
+ int N2 = connec[ipt+2];
+ int N3 = connec[ipt+3];
+ int N4 = connec[ipt+4];
+ int N5 = connec[ipt+5];
+
+ area_vol[iel]=INTERP_KERNEL::calculateVolumeForPyra(coords+dim_space*N1,
+ coords+dim_space*N2,
+ coords+dim_space*N3,
+ coords+dim_space*N4,
+ coords+dim_space*N5) ;
+ }
+ break ;
+
+ case INTERP_KERNEL::NORM_PENTA6 :
+ case INTERP_KERNEL::NORM_PENTA15 :
+ {
+ int N1 = connec[ipt+1];
+ int N2 = connec[ipt+2];
+ int N3 = connec[ipt+3];
+ int N4 = connec[ipt+4];
+ int N5 = connec[ipt+5];
+ int N6 = connec[ipt+6];
+
+ area_vol[iel]=INTERP_KERNEL::calculateVolumeForPenta(coords+dim_space*N1,
+ coords+dim_space*N2,
+ coords+dim_space*N3,
+ coords+dim_space*N4,
+ coords+dim_space*N5,
+ coords+dim_space*N6) ;
+ }
+ break ;
+
+ case INTERP_KERNEL::NORM_HEXA8 :
+ case INTERP_KERNEL::NORM_HEXA20 :
+ {
+ int N1 = connec[ipt+1];
+ int N2 = connec[ipt+2];
+ int N3 = connec[ipt+3];
+ int N4 = connec[ipt+4];
+ int N5 = connec[ipt+5];
+ int N6 = connec[ipt+6];
+ int N7 = connec[ipt+7];
+ int N8 = connec[ipt+8];
+
+ area_vol[iel]=INTERP_KERNEL::calculateVolumeForHexa(coords+dim_space*N1,
+ coords+dim_space*N2,
+ coords+dim_space*N3,
+ coords+dim_space*N4,
+ coords+dim_space*N5,
+ coords+dim_space*N6,
+ coords+dim_space*N7,
+ coords+dim_space*N8) ;
+ }
+ break ;
+
+ case INTERP_KERNEL::NORM_POLYHED :
+ {
+ throw INTERP_KERNEL::Exception("Not yet implemented !");
+ }
+ break ;
+
+ default:
+ throw INTERP_KERNEL::Exception("Bad Support to get Volume on it !");
+ }
+ }
+ break;
+ default:
+ throw INTERP_KERNEL::Exception("interpolation is not available for this dimension");
+ }
+
+ field->setArray(array) ;
+ array->decrRef();
+ field->setMesh(const_cast<MEDCouplingUMesh *>(this));
+ return field ;
+}
#define __PARAMEDMEM_MEDCOUPLINGUMESH_HXX__
#include "MEDCoupling.hxx"
-#include "MEDCouplingMesh.hxx"
+#include "MEDCouplingPointSet.hxx"
#include "MemArray.hxx"
#include <set>
namespace ParaMEDMEM
{
- class MEDCOUPLING_EXPORT MEDCouplingUMesh : public MEDCouplingMesh
+ class MEDCOUPLING_EXPORT MEDCouplingUMesh : public MEDCouplingPointSet
{
public:
static MEDCouplingUMesh *New();
void checkCoherency() const throw(INTERP_KERNEL::Exception);
void setMeshDimension(unsigned meshDim);
void allocateCells(int nbOfCells);
- void setCoords(DataArrayDouble *coords);
- DataArrayDouble *getCoords() const { return _coords; }
void insertNextCell(INTERP_KERNEL::NormalizedCellType type, int size, const int *nodalConnOfCell);
void finishInsertingCells();
const std::set<INTERP_KERNEL::NormalizedCellType>& getAllTypes() const { return _types; }
DataArrayInt *getNodalConnectivityIndex() const { return _nodal_connec_index; }
INTERP_KERNEL::NormalizedCellType getTypeOfCell(int cellId) const;
int getNumberOfNodesInCell(int cellId) const;
- bool isStructured() const;
int getNumberOfCells() const;
- int getNumberOfNodes() const;
- int getSpaceDimension() const;
int getMeshDimension() const { return _mesh_dim; }
int getMeshLength() const;
+ //! size of returned tinyInfo must be always the same.
+ void getTinySerializationInformation(std::vector<int>& tinyInfo) const;
+ void resizeForSerialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2);
+ void serialize(DataArrayInt *&a1, DataArrayDouble *&a2);
+ MEDCouplingPointSet *buildObjectFromUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2);
//tools
void zipCoords();
DataArrayInt *zipCoordsTraducer();
void getReverseNodalConnectivity(DataArrayInt *revNodal, DataArrayInt *revNodalIndx) const;
- MEDCouplingUMesh *buildPartOfMySelf(const int *start, const int *end, bool keepCoords) const;
+ MEDCouplingPointSet *buildPartOfMySelf(const int *start, const int *end, bool keepCoords) const;
+ void giveElemsInBoundingBox(const double *bbox, double eps, std::vector<int>& elems);
+ MEDCouplingFieldDouble *getMeasureField() const;
private:
MEDCouplingUMesh();
MEDCouplingUMesh(const MEDCouplingUMesh& other, bool deepCpy);
unsigned _mesh_dim;
DataArrayInt *_nodal_connec;
DataArrayInt *_nodal_connec_index;
- DataArrayDouble *_coords;
std::set<INTERP_KERNEL::NormalizedCellType> _types;
private:
static const char PART_OF_NAME[];
--- /dev/null
+// Copyright (C) 2007-2008 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 "MEDCouplingUMeshDesc.hxx"
+#include "CellModel.hxx"
+#include "MemArray.hxx"
+
+using namespace ParaMEDMEM;
+
+MEDCouplingUMeshDesc::MEDCouplingUMeshDesc():_mesh_dim(-1),_desc_connec(0),_desc_connec_index(0),
+ _nodal_connec_face(0),_nodal_connec_face_index(0)
+{
+}
+
+MEDCouplingUMeshDesc::~MEDCouplingUMeshDesc()
+{
+ if(_desc_connec)
+ _desc_connec->decrRef();
+ if(_desc_connec_index)
+ _desc_connec_index->decrRef();
+ if(_nodal_connec_face)
+ _nodal_connec_face->decrRef();
+ if(_nodal_connec_face_index)
+ _nodal_connec_face_index->decrRef();
+}
+
+MEDCouplingUMeshDesc *MEDCouplingUMeshDesc::New()
+{
+ return new MEDCouplingUMeshDesc;
+}
+
+void MEDCouplingUMeshDesc::checkCoherency() const throw(INTERP_KERNEL::Exception)
+{
+ for(std::set<INTERP_KERNEL::NormalizedCellType>::const_iterator iter=_types.begin();iter!=_types.end();iter++)
+ {
+ if(INTERP_KERNEL::CellModel::getCellModel(*iter).getDimension()!=_mesh_dim)
+ {
+ std::ostringstream message;
+ message << "MeshDesc invalid because dimension is " << _mesh_dim << " and there is presence of cell(s) with type " << (*iter);
+ throw INTERP_KERNEL::Exception(message.str().c_str());
+ }
+ }
+}
+
+void MEDCouplingUMeshDesc::setMeshDimension(unsigned meshDim)
+{
+ _mesh_dim=meshDim;
+ declareAsNew();
+}
+
+int MEDCouplingUMeshDesc::getNumberOfCells() const
+{
+ if(_desc_connec_index)
+ return _desc_connec_index->getNumberOfTuples()-1;
+ else
+ throw INTERP_KERNEL::Exception("Unable to get number of cells because no connectivity specified !");
+}
+
+int MEDCouplingUMeshDesc::getNumberOfFaces() const
+{
+ if(_nodal_connec_face_index)
+ return _nodal_connec_face_index->getNumberOfTuples()-1;
+ else
+ throw INTERP_KERNEL::Exception("Unable to get number of faces because no connectivity specified !");
+}
+
+int MEDCouplingUMeshDesc::getCellMeshLength() const
+{
+ return _desc_connec->getNbOfElems();
+}
+
+int MEDCouplingUMeshDesc::getFaceMeshLength() const
+{
+ return _nodal_connec_face->getNbOfElems();
+}
+
+void MEDCouplingUMeshDesc::setConnectivity(DataArrayInt *descConn, DataArrayInt *descConnIndex, DataArrayInt *nodalFaceConn, DataArrayInt *nodalFaceConnIndx)
+{
+ DataArrayInt::setArrayIn(descConn,_desc_connec);
+ DataArrayInt::setArrayIn(descConnIndex,_desc_connec_index);
+ DataArrayInt::setArrayIn(nodalFaceConn,_nodal_connec_face);
+ DataArrayInt::setArrayIn(nodalFaceConnIndx,_nodal_connec_face_index);
+ computeTypes();
+}
+
+void MEDCouplingUMeshDesc::getTinySerializationInformation(std::vector<int>& tinyInfo) const
+{
+ tinyInfo.resize(7);
+ tinyInfo[0]=getSpaceDimension();
+ tinyInfo[1]=getMeshDimension();
+ tinyInfo[2]=getNumberOfNodes();
+ tinyInfo[3]=getNumberOfCells();
+ tinyInfo[4]=getCellMeshLength();
+ tinyInfo[5]=getNumberOfFaces();
+ tinyInfo[6]=getFaceMeshLength();
+}
+
+void MEDCouplingUMeshDesc::resizeForSerialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2)
+{
+ a1->alloc(tinyInfo[4]+tinyInfo[3]+1+tinyInfo[6]+tinyInfo[5]+1,1);
+ a2->alloc(tinyInfo[2],tinyInfo[0]);
+}
+
+void MEDCouplingUMeshDesc::serialize(DataArrayInt *&a1, DataArrayDouble *&a2)
+{
+ a1=DataArrayInt::New();
+ a1->alloc(getCellMeshLength()+getNumberOfCells()+1+getFaceMeshLength()+getNumberOfFaces()+1,1);
+ int *ptA1=a1->getPointer();
+ const int *descConn=_desc_connec->getConstPointer();
+ const int *descConnIndex=_desc_connec_index->getConstPointer();
+ const int *faceConn=_nodal_connec_face->getConstPointer();
+ const int *faceConnIndex=_nodal_connec_face_index->getConstPointer();
+ ptA1=std::copy(descConn,descConn+getCellMeshLength(),ptA1);
+ ptA1=std::copy(descConnIndex,descConnIndex+getNumberOfCells()+1,ptA1);
+ ptA1=std::copy(faceConn,faceConn+getFaceMeshLength(),ptA1);
+ std::copy(faceConnIndex,faceConnIndex+getNumberOfFaces()+1,ptA1);
+ a2=getCoords();
+ a2->incrRef();
+}
+
+MEDCouplingPointSet *MEDCouplingUMeshDesc::buildObjectFromUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2)
+{
+ MEDCouplingUMeshDesc *meshing=MEDCouplingUMeshDesc::New();
+ meshing->setCoords(a2);
+ const int *recvBuffer=a1->getConstPointer();
+ DataArrayInt *descConn=DataArrayInt::New();
+ descConn->alloc(tinyInfo[4],1);
+ std::copy(recvBuffer,recvBuffer+tinyInfo[4],descConn->getPointer());
+ DataArrayInt *descConnIndex=DataArrayInt::New();
+ descConnIndex->alloc(tinyInfo[3]+1,1);
+ std::copy(recvBuffer+tinyInfo[4],recvBuffer+tinyInfo[4]+tinyInfo[3]+1,descConnIndex->getPointer());
+ DataArrayInt *faceConn=DataArrayInt::New();
+ faceConn->alloc(tinyInfo[6],1);
+ std::copy(recvBuffer+tinyInfo[4]+tinyInfo[3]+1,recvBuffer+tinyInfo[4]+tinyInfo[3]+1+tinyInfo[6],faceConn->getPointer());
+ DataArrayInt *faceConnIndex=DataArrayInt::New();
+ faceConnIndex->alloc(tinyInfo[5]+1,1);
+ std::copy(recvBuffer+tinyInfo[4]+tinyInfo[3]+1+tinyInfo[6],
+ recvBuffer+tinyInfo[4]+tinyInfo[3]+1+tinyInfo[6]+tinyInfo[5]+1,faceConnIndex->getPointer());
+ meshing->setConnectivity(descConn,descConnIndex,faceConn,faceConnIndex);
+ descConn->decrRef();
+ descConnIndex->decrRef();
+ faceConn->decrRef();
+ faceConnIndex->decrRef();
+ meshing->setMeshDimension(tinyInfo[1]);
+ return meshing;
+}
+
+void MEDCouplingUMeshDesc::giveElemsInBoundingBox(const double *bbox, double eps, std::vector<int>& elems)
+{
+ int dim=getSpaceDimension();
+ double* elem_bb=new double[2*dim];
+ const int* conn = _desc_connec->getConstPointer();
+ const int* conn_index= _desc_connec_index->getConstPointer();
+ const int* face = _nodal_connec_face->getConstPointer();
+ const int* face_index= _nodal_connec_face_index->getConstPointer();
+ const double* coords = getCoords()->getConstPointer();
+ int nbOfCells=getNumberOfCells();
+ for ( int ielem=0; ielem<nbOfCells;ielem++ )
+ {
+ for (int i=0; i<dim; i++)
+ {
+ elem_bb[i*2]=std::numeric_limits<double>::max();
+ elem_bb[i*2+1]=-std::numeric_limits<double>::max();
+ }
+
+ for (int iface=conn_index[ielem]+1; iface<conn_index[ielem+1]; iface++)//+1 due to offset of cell type.
+ {
+ for(int inode=face_index[iface]+1;inode<face_index[iface+1];inode++)
+ {
+ int node=face[inode];
+ for (int idim=0; idim<dim; idim++)
+ {
+ if ( coords[node*dim+idim] < elem_bb[idim*2] )
+ {
+ elem_bb[idim*2] = coords[node*dim+idim] ;
+ }
+ if ( coords[node*dim+idim] > elem_bb[idim*2+1] )
+ {
+ elem_bb[idim*2+1] = coords[node*dim+idim] ;
+ }
+ }
+ }
+ }
+ if (intersectsBoundingBox(elem_bb, bbox, dim, eps))
+ {
+ elems.push_back(ielem);
+ }
+ }
+ delete [] elem_bb;
+}
+
+MEDCouplingPointSet *MEDCouplingUMeshDesc::buildPartOfMySelf(const int *start, const int *end, bool keepCoords) const
+{
+ //not implemented yet.
+ return 0;
+}
+
+MEDCouplingFieldDouble *MEDCouplingUMeshDesc::getMeasureField() const
+{
+ //not implemented yet.
+ return 0;
+}
+
+void MEDCouplingUMeshDesc::computeTypes()
+{
+ if(_desc_connec && _desc_connec_index)
+ {
+ _types.clear();
+ const int *conn=_desc_connec->getConstPointer();
+ const int *connIndex=_desc_connec_index->getConstPointer();
+ int nbOfElem=_desc_connec_index->getNbOfElems()-1;
+ for(const int *pt=connIndex;pt!=connIndex+nbOfElem;pt++)
+ _types.insert((INTERP_KERNEL::NormalizedCellType)conn[*pt]);
+ }
+}
--- /dev/null
+// Copyright (C) 2007-2008 CEA/DEN, EDF R&D
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+#ifndef __PARAMEDMEM_MEDCOUPLINGUMESHDESC_HXX__
+#define __PARAMEDMEM_MEDCOUPLINGUMESHDESC_HXX__
+
+#include "MEDCouplingPointSet.hxx"
+#include "MEDCoupling.hxx"
+#include "NormalizedUnstructuredMesh.hxx"
+
+#include <set>
+
+namespace ParaMEDMEM
+{
+ class MEDCOUPLING_EXPORT MEDCouplingUMeshDesc : public MEDCouplingPointSet
+ {
+ public:
+ static MEDCouplingUMeshDesc *New();
+ void checkCoherency() const throw(INTERP_KERNEL::Exception);
+ void setMeshDimension(unsigned meshDim);
+ int getNumberOfCells() const;
+ int getNumberOfFaces() const;
+ int getCellMeshLength() const;
+ int getFaceMeshLength() const;
+ int getMeshDimension() const { return _mesh_dim; }
+ void setConnectivity(DataArrayInt *descConn, DataArrayInt *descConnIndex, DataArrayInt *nodalFaceConn, DataArrayInt *nodalFaceConnIndx);
+ //tools to overload
+ void getTinySerializationInformation(std::vector<int>& tinyInfo) const;
+ void resizeForSerialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2);
+ void serialize(DataArrayInt *&a1, DataArrayDouble *&a2);
+ MEDCouplingPointSet *buildObjectFromUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2);
+ void giveElemsInBoundingBox(const double *bbox, double eps, std::vector<int>& elems);
+ MEDCouplingPointSet *buildPartOfMySelf(const int *start, const int *end, bool keepCoords) const;
+ MEDCouplingFieldDouble *getMeasureField() const;
+ private:
+ MEDCouplingUMeshDesc();
+ ~MEDCouplingUMeshDesc();
+ void computeTypes();
+ private:
+ unsigned _mesh_dim;
+ DataArrayInt *_desc_connec;
+ DataArrayInt *_desc_connec_index;
+ DataArrayInt *_nodal_connec_face;
+ DataArrayInt *_nodal_connec_face_index;
+ std::set<INTERP_KERNEL::NormalizedCellType> _types;
+ };
+}
+
+#endif
MEDCouplingFieldDouble.hxx MEDCouplingMesh.hxx MEDCouplingUMesh.hxx TimeLabel.hxx \
MEDCouplingField.hxx MEDCouplingNormalizedUnstructuredMesh.hxx MemArray.hxx \
MEDCouplingNormalizedUnstructuredMesh.txx MemArray.txx RefCountObject.hxx \
-MEDCouplingSMesh.hxx
+MEDCouplingRMesh.hxx MEDCouplingTimeDiscretization.hxx \
+MEDCouplingFieldDiscretization.hxx MEDCouplingPointSet.hxx \
+MEDCouplingUMeshDesc.hxx
# Libraries targets
dist_libmedcoupling_la_SOURCES = \
- MEDCouplingField.cxx MEDCouplingFieldDouble.cxx \
- MEDCouplingUMesh.cxx MemArray.cxx TimeLabel.cxx \
- MEDCouplingSMesh.cxx
+ MEDCouplingField.cxx MEDCouplingFieldDouble.cxx \
+ MEDCouplingUMesh.cxx MemArray.cxx TimeLabel.cxx \
+ MEDCouplingRMesh.cxx MEDCouplingTimeDiscretization.cxx \
+ MEDCouplingFieldDiscretization.cxx \
+ MEDCouplingPointSet.cxx MEDCouplingUMeshDesc.cxx
libmedcoupling_la_LDFLAGS=
EXTRA_DIST += \
MEDCouplingFieldDouble.hxx \
+ MEDCouplingTimeDiscretization.hxx \
+ MEDCouplingFieldDiscretization.hxx \
MEDCouplingMesh.hxx \
MEDCouplingUMesh.hxx \
- MEDCouplingSMesh.hxx \
+ MEDCouplingRMesh.hxx \
TimeLabel.hxx \
MEDCouplingField.hxx \
MEDCouplingNormalizedUnstructuredMesh.hxx \
MemArray.hxx \
MEDCouplingNormalizedUnstructuredMesh.txx \
MemArray.txx \
- RefCountObject.hxx
+ RefCountObject.hxx \
+ MEDCouplingPointSet.hxx \
+ MEDCouplingUMeshDesc.hxx
declareAsNew();
}
-void DataArrayDouble::useArray(double *array, bool ownership, DeallocType type, int nbOfTuple, int nbOfCompo)
+void DataArrayDouble::setArrayIn(DataArrayDouble *newArray, DataArrayDouble* &arrayToSet)
+{
+ if(newArray!=arrayToSet)
+ {
+ if(arrayToSet)
+ arrayToSet->decrRef();
+ arrayToSet=newArray;
+ if(arrayToSet)
+ arrayToSet->incrRef();
+ }
+}
+
+void DataArrayDouble::useArray(const double *array, bool ownership, DeallocType type, int nbOfTuple, int nbOfCompo)
{
_nb_of_tuples=nbOfTuple;
_info_on_compo.resize(nbOfCompo);
declareAsNew();
}
-void DataArrayInt::useArray(int *array, bool ownership, DeallocType type, int nbOfTuple, int nbOfCompo)
+void DataArrayInt::useArray(const int *array, bool ownership, DeallocType type, int nbOfTuple, int nbOfCompo)
{
_nb_of_tuples=nbOfTuple;
_info_on_compo.resize(nbOfCompo);
_nb_of_tuples=nbOfTuples;
declareAsNew();
}
+
+void DataArrayInt::setArrayIn(DataArrayInt *newArray, DataArrayInt* &arrayToSet)
+{
+ if(newArray!=arrayToSet)
+ {
+ if(arrayToSet)
+ arrayToSet->decrRef();
+ arrayToSet=newArray;
+ if(arrayToSet)
+ arrayToSet->incrRef();
+ }
+}
#include "MEDCoupling.hxx"
#include "RefCountObject.hxx"
+#include "InterpKernelException.hxx"
#include <string>
#include <vector>
namespace ParaMEDMEM
{
+ template<class T>
+ class MEDCouplingPointer
+ {
+ public:
+ MEDCouplingPointer():_internal(0),_external(0) { }
+ void null() { _internal=0; _external=0; }
+ bool isNull() const { return _internal==0 && _external==0; }
+ void setInternal(T *pointer);
+ void setExternal(const T *pointer);
+ const T *getConstPointer() const { if(_internal) return _internal; else return _external; }
+ const T *getConstPointerLoc(int offset) const { if(_internal) return _internal+offset; else return _external+offset; }
+ T *getPointer() const { if(_internal) return _internal; throw INTERP_KERNEL::Exception("Trying to write on an external pointer."); }
+ private:
+ T *_internal;
+ const T *_external;
+ };
+
template<class T>
class MemArray
{
public:
- MemArray():_nb_of_elem(-1),_ownership(false),_pointer(0),_dealloc(CPP_DEALLOC) { }
+ MemArray():_nb_of_elem(-1),_ownership(false),_dealloc(CPP_DEALLOC) { }
MemArray(const MemArray<T>& other);
- T *getPointer() const { return _pointer; }
+ const T *getConstPointerLoc(int offset) const { return _pointer.getConstPointerLoc(offset); }
+ const T *getConstPointer() const { return _pointer.getConstPointer(); }
+ T *getPointer() const { return _pointer.getPointer(); }
MemArray<T> &operator=(const MemArray<T>& other);
- T operator[](int id) const { return _pointer[id]; }
- T& operator[](int id) { return _pointer[id]; }
+ T operator[](int id) const { return _pointer.getConstPointer()[id]; }
+ T& operator[](int id) { return _pointer.getPointer()[id]; }
void alloc(int nbOfElements);
void reAlloc(int newNbOfElements);
- void useArray(void *array, bool ownership, DeallocType type, int nbOfElem);
+ void useArray(const T *array, bool ownership, DeallocType type, int nbOfElem);
void writeOnPlace(int id, T element0, const T *others, int sizeOfOthers);
~MemArray() { destroy(); }
private:
private:
int _nb_of_elem;
bool _ownership;
- T *_pointer;
+ MEDCouplingPointer<T> _pointer;
+ //T *_pointer;
DeallocType _dealloc;
};
bool isEqual(DataArrayDouble *other, double prec) const;
//!alloc or useArray should have been called before.
void reAlloc(int nbOfTuples);
+ void getTuple(int tupleId, double *res) const { std::copy(_mem.getConstPointerLoc(tupleId*_info_on_compo.size()),_mem.getConstPointerLoc((tupleId+1)*_info_on_compo.size()),res); }
double getIJ(int tupleId, int compoId) const { return _mem[tupleId*_info_on_compo.size()+compoId]; }
void setIJ(int tupleId, int compoId, double newVal) { _mem[tupleId*_info_on_compo.size()+compoId]=newVal; }
double *getPointer() const { return _mem.getPointer(); }
- void useArray(double *array, bool ownership, DeallocType type, int nbOfTuple, int nbOfCompo);
+ static void setArrayIn(DataArrayDouble *newArray, DataArrayDouble* &arrayToSet);
+ const double *getConstPointer() const { return _mem.getConstPointer(); }
+ void useArray(const double *array, bool ownership, DeallocType type, int nbOfTuple, int nbOfCompo);
void writeOnPlace(int id, double element0, const double *others, int sizeOfOthers) { _mem.writeOnPlace(id,element0,others,sizeOfOthers); }
//! nothing to do here because this class does not aggregate any TimeLabel instance.
void updateTime() { }
void alloc(int nbOfTuple, int nbOfCompo);
//!alloc or useArray should have been called before.
void reAlloc(int nbOfTuples);
+ void getTuple(int tupleId, int *res) const { std::copy(_mem.getConstPointerLoc(tupleId*_info_on_compo.size()),_mem.getConstPointerLoc((tupleId+1)*_info_on_compo.size()),res); }
int getIJ(int tupleId, int compoId) const { return _mem[tupleId*_info_on_compo.size()+compoId]; }
void setIJ(int tupleId, int compoId, int newVal) { _mem[tupleId*_info_on_compo.size()+compoId]=newVal; }
int *getPointer() const { return _mem.getPointer(); }
- void useArray(int *array, bool ownership, DeallocType type, int nbOfTuple, int nbOfCompo);
+ static void setArrayIn(DataArrayInt *newArray, DataArrayInt* &arrayToSet);
+ const int *getConstPointer() const { return _mem.getConstPointer(); }
+ void useArray(const int *array, bool ownership, DeallocType type, int nbOfTuple, int nbOfCompo);
void writeOnPlace(int id, int element0, const int *others, int sizeOfOthers) { _mem.writeOnPlace(id,element0,others,sizeOfOthers); }
//! nothing to do here because this class does not aggregate any TimeLabel instance.
void updateTime() { }
namespace ParaMEDMEM
{
template<class T>
- MemArray<T>::MemArray(const MemArray<T>& other):_nb_of_elem(-1),_ownership(false),_pointer(0),_dealloc(CPP_DEALLOC)
+ void MEDCouplingPointer<T>::setInternal(T *pointer)
{
- if(other._pointer)
+ _internal=pointer;
+ _external=0;
+ }
+
+ template<class T>
+ void MEDCouplingPointer<T>::setExternal(const T *pointer)
+ {
+ _external=pointer;
+ _internal=0;
+ }
+
+ template<class T>
+ MemArray<T>::MemArray(const MemArray<T>& other):_nb_of_elem(-1),_ownership(false),_dealloc(CPP_DEALLOC)
+ {
+ if(!other._pointer.isNull())
{
T *pointer=new T[other._nb_of_elem];
- std::copy(other._pointer,other._pointer+other._nb_of_elem,pointer);
+ std::copy(other._pointer.getConstPointer(),other._pointer.getConstPointer()+other._nb_of_elem,pointer);
useArray(pointer,true,CPP_DEALLOC,other._nb_of_elem);
}
}
template<class T>
- void MemArray<T>::useArray(void *array, bool ownership, DeallocType type, int nbOfElem)
+ void MemArray<T>::useArray(const T *array, bool ownership, DeallocType type, int nbOfElem)
{
_nb_of_elem=nbOfElem;
destroy();
- _pointer=(T *)array;
+ if(ownership)
+ _pointer.setInternal((T *)array);
+ else
+ _pointer.setExternal(array);
_ownership=ownership;
_dealloc=type;
}
{
if(id+sizeOfOthers>=_nb_of_elem)
reAlloc(2*_nb_of_elem+sizeOfOthers+1);
- _pointer[id]=element0;
- memcpy(_pointer+id+1,others,sizeOfOthers*sizeof(T));
+ T *pointer=_pointer.getPointer();
+ pointer[id]=element0;
+ std::copy(others,others+sizeOfOthers,pointer+id+1);
}
template<class T>
{
destroy();
_nb_of_elem=nbOfElements;
- _pointer=new T[_nb_of_elem];
+ _pointer.setInternal(new T[_nb_of_elem]);
_ownership=true;
_dealloc=CPP_DEALLOC;
}
void MemArray<T>::reAlloc(int newNbOfElements)
{
T *pointer=new T[newNbOfElements];
- memcpy(pointer,_pointer,std::min<int>(_nb_of_elem,newNbOfElements)*sizeof(int));
- destroyPointer(_pointer,_dealloc);
- _pointer=pointer;
+ std::copy(_pointer.getConstPointer(),_pointer.getConstPointer()+std::min<int>(_nb_of_elem,newNbOfElements),pointer);
+ if(_ownership)
+ destroyPointer((T *)_pointer.getConstPointer(),_dealloc);
+ _pointer.setInternal(pointer);
_nb_of_elem=newNbOfElements;
_ownership=true;
_dealloc=CPP_DEALLOC;
void MemArray<T>::destroy()
{
if(_ownership)
- destroyPointer(_pointer,_dealloc);
- _pointer=0;
+ destroyPointer((T *)_pointer.getConstPointer(),_dealloc);
+ _pointer.null();
_ownership=false;
}
MemArray<T> &MemArray<T>::operator=(const MemArray<T>& other)
{
alloc(other._nb_of_elem);
- memcpy(_pointer,other._pointer,_nb_of_elem*sizeof(T));
+ std::copy(other._pointer.getConstPointer(),other._pointer.getConstPointer()+_nb_of_elem,_pointer.getPointer());
return *this;
}
}
ON_NODES = 1
} TypeOfField;
+ typedef enum
+ {
+ NO_TIME = 4,
+ ONE_TIME = 5,
+ LINEAR_TIME = 6
+ } TypeOfTimeDiscretization;
+
class RefCountObject : public TimeLabel
{
protected:
#include "MemArray.hxx"
#include "Interpolation2D.txx"
#include "Interpolation3DSurf.txx"
+#include "Interpolation3D.txx"
#include "MEDCouplingNormalizedUnstructuredMesh.txx"
using namespace std;
using namespace ParaMEDMEM;
+void MEDCouplingBasicsTest::testArray()
+{
+ int tmp1[6]={7,6,5,4,3,2};
+ const int tmp2[3]={8,9,10};
+ {
+ MemArray<int> mem;
+ mem.useArray(tmp1,false,CPP_DEALLOC,6);
+ CPPUNIT_ASSERT(tmp1==mem.getConstPointer());
+ CPPUNIT_ASSERT_THROW(mem.getPointer(),INTERP_KERNEL::Exception);
+ CPPUNIT_ASSERT_THROW(mem[2]=7,INTERP_KERNEL::Exception);
+ CPPUNIT_ASSERT_THROW(mem.writeOnPlace(0,12,tmp2,3),INTERP_KERNEL::Exception);
+ mem.writeOnPlace(4,12,tmp2,3);
+ }
+ {
+ int *tmp3=new int[6];
+ std::copy(tmp1,tmp1+6,tmp3);
+ MemArray<int> mem2;
+ mem2.useArray(tmp3,true,CPP_DEALLOC,6);
+ CPPUNIT_ASSERT(tmp3==mem2.getConstPointer());
+ CPPUNIT_ASSERT(tmp3==mem2.getPointer());
+ CPPUNIT_ASSERT_EQUAL(5,mem2[2]);
+ mem2[2]=7;
+ CPPUNIT_ASSERT_EQUAL(7,mem2[2]);
+ mem2.writeOnPlace(0,12,tmp2,3);
+ CPPUNIT_ASSERT_EQUAL(9,mem2[2]);
+ CPPUNIT_ASSERT_EQUAL(12,mem2[0]);
+ mem2.writeOnPlace(4,12,tmp2,3);
+ }
+}
+
void MEDCouplingBasicsTest::testMesh()
{
const int nbOfCells=6;
mesh->decrRef();
}
+void MEDCouplingBasicsTest::testMeshPointsCloud()
+{
+ double targetCoords[27]={-0.3,-0.3,0.5, 0.2,-0.3,1., 0.7,-0.3,1.5, -0.3,0.2,0.5, 0.2,0.2,1., 0.7,0.2,1.5, -0.3,0.7,0.5, 0.2,0.7,1., 0.7,0.7,1.5};
+ int *targetConn=0;
+ MEDCouplingUMesh *targetMesh=MEDCouplingUMesh::New();
+ targetMesh->setMeshDimension(0);
+ targetMesh->allocateCells(8);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_POINT0,0,targetConn);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_POINT0,0,targetConn);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_POINT0,0,targetConn);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_POINT0,0,targetConn);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_POINT0,0,targetConn);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_POINT0,0,targetConn);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_POINT0,0,targetConn);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_POINT0,0,targetConn);
+ targetMesh->finishInsertingCells();
+ DataArrayDouble *myCoords=DataArrayDouble::New();
+ myCoords->alloc(9,3);
+ std::copy(targetCoords,targetCoords+27,myCoords->getPointer());
+ targetMesh->setCoords(myCoords);
+ myCoords->decrRef();
+ //
+ targetMesh->checkCoherency();
+ CPPUNIT_ASSERT_EQUAL(3,targetMesh->getSpaceDimension());
+ CPPUNIT_ASSERT_EQUAL(8,targetMesh->getNumberOfCells());
+ CPPUNIT_ASSERT_EQUAL(9,targetMesh->getNumberOfNodes());
+ CPPUNIT_ASSERT_EQUAL(0,targetMesh->getMeshDimension());
+ //
+ targetMesh->decrRef();
+}
+
void MEDCouplingBasicsTest::testDeepCopy()
{
DataArrayDouble *array=DataArrayDouble::New();
const int tab1[2]={0,4};
const int tab2[3]={0,2,3};
//
- MEDCouplingUMesh *subMesh=mesh->buildPartOfMySelf(tab1,tab1+2,true);
+ MEDCouplingPointSet *subMeshSimple=mesh->buildPartOfMySelf(tab1,tab1+2,true);
+ MEDCouplingUMesh *subMesh=dynamic_cast<MEDCouplingUMesh *>(subMeshSimple);
+ CPPUNIT_ASSERT(subMesh);
std::string name(subMesh->getName());
CPPUNIT_ASSERT_EQUAL(2,(int)mesh->getAllTypes().size());
CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_TRI3,*mesh->getAllTypes().begin());
CPPUNIT_ASSERT(std::equal(subConnIndex,subConnIndex+3,subMesh->getNodalConnectivityIndex()->getPointer()));
subMesh->decrRef();
//
- subMesh=mesh->buildPartOfMySelf(tab2,tab2+3,true);
+ subMeshSimple=mesh->buildPartOfMySelf(tab2,tab2+3,true);
+ subMesh=dynamic_cast<MEDCouplingUMesh *>(subMeshSimple);
+ CPPUNIT_ASSERT(subMesh);
name=subMesh->getName();
CPPUNIT_ASSERT_EQUAL(2,(int)subMesh->getAllTypes().size());
CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_TRI3,*subMesh->getAllTypes().begin());
oldCoords->decrRef();
//
const int tab1[2]={0,4};
- MEDCouplingUMesh *subMesh=mesh->buildPartOfMySelf(tab1,tab1+2,true);
+ MEDCouplingPointSet *subMeshPtSet=mesh->buildPartOfMySelf(tab1,tab1+2,true);
+ MEDCouplingUMesh *subMesh=dynamic_cast<MEDCouplingUMesh *>(subMeshPtSet);
+ CPPUNIT_ASSERT(subMesh);
DataArrayInt *traducer=subMesh->zipCoordsTraducer();
const int expectedTraducer[9]={0,1,-1,2,3,4,-1,5,6};
CPPUNIT_ASSERT(std::equal(expectedTraducer,expectedTraducer+9,traducer->getPointer()));
CPPUNIT_ASSERT(std::equal(subConnIndex,subConnIndex+3,subMesh->getNodalConnectivityIndex()->getPointer()));
subMesh->decrRef();
//
- subMesh=mesh->buildPartOfMySelf(tab1,tab1+2,false);
+ subMeshPtSet=mesh->buildPartOfMySelf(tab1,tab1+2,false);
+ subMesh=dynamic_cast<MEDCouplingUMesh *>(subMeshPtSet);
+ CPPUNIT_ASSERT(subMesh);
CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_QUAD4,*subMesh->getAllTypes().begin());
CPPUNIT_ASSERT_EQUAL(2,subMesh->getNumberOfCells());
CPPUNIT_ASSERT_EQUAL(7,subMesh->getNumberOfNodes());
targetMesh->decrRef();
}
+void MEDCouplingBasicsTest::test2DInterpP1P1_1()
+{
+ MEDCouplingUMesh *sourceMesh=build2DSourceMesh_1();
+ MEDCouplingUMesh *targetMesh=build2DTargetMesh_2();
+ //
+ MEDCouplingNormalizedUnstructuredMesh<2,2> sourceWrapper(sourceMesh);
+ MEDCouplingNormalizedUnstructuredMesh<2,2> targetWrapper(targetMesh);
+ INTERP_KERNEL::Interpolation2D myInterpolator;
+ vector<map<int,double> > res;
+ INTERP_KERNEL::IntersectionType types[2]={INTERP_KERNEL::Triangulation, INTERP_KERNEL::Geometric2D};
+ for(int i=0;i<2;i++)
+ {
+ myInterpolator.setPrecision(1e-12);
+ myInterpolator.setIntersectionType(types[i]);
+ myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,res,"P1P1");
+ CPPUNIT_ASSERT_EQUAL(9,(int)res.size());
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.08333333333333334,res[0][0],1.e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.05416666666666665,res[1][0],1.e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.02916666666666666,res[1][1],1.e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.08333333333333334,res[2][1],1.e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.05416666666666665,res[3][0],1.e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.02916666666666668,res[3][2],1.e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.1416666666666666,res[4][0],1.e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.02499999999999999,res[4][1],1.e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.02499999999999999,res[4][2],1.e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.09999999999999999,res[4][3],1.e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.02916666666666666,res[5][1],1.e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.09583333333333333,res[5][3],1.e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.08333333333333333,res[6][2],1.e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.02916666666666667,res[7][2],1.e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.09583333333333331,res[7][3],1.e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.04166666666666668,res[8][3],1.e-12);
+ res.clear();
+ }
+ //clean up
+ sourceMesh->decrRef();
+ targetMesh->decrRef();
+}
+
void MEDCouplingBasicsTest::test3DSurfInterpP0P0_1()
{
MEDCouplingUMesh *sourceMesh=build3DSurfSourceMesh_1();
void MEDCouplingBasicsTest::test3DInterpP0P0_1()
{
MEDCouplingUMesh *sourceMesh=build3DSourceMesh_1();
+ MEDCouplingUMesh *targetMesh=build3DTargetMesh_1();
+ //
+ MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(sourceMesh);
+ MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(targetMesh);
+ INTERP_KERNEL::Interpolation3D myInterpolator;
+ vector<map<int,double> > res;
+ myInterpolator.setPrecision(1e-12);
+ myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,res,"P0P0");
+ CPPUNIT_ASSERT_EQUAL(8,(int)res.size());
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(8.e6,sumAll(res),1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(20833.33333333333,res[0][0],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(41666.66666666667,res[0][6],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(20833.33333333333,res[0][7],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(20833.33333333333,res[0][8],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(20833.33333333333,res[0][10],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(41666.66666666667,res[1][2],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(166666.6666666667,res[1][7],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(166666.6666666667,res[1][8],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(166666.6666666667,res[2][0],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(20833.33333333333,res[2][5],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(145833.3333333333,res[2][6],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(20833.33333333333,res[2][9],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(20833.33333333333,res[2][11],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(395833.3333333333,res[3][0],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(145833.3333333333,res[3][2],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(20833.33333333331,res[3][3],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(166666.6666666667,res[3][5],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(395833.3333333333,res[3][8],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(20833.33333333333,res[4][1],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(20833.33333333333,res[4][4],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(145833.3333333333,res[4][6],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(20833.33333333333,res[4][9],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(166666.6666666667,res[4][10],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(145833.3333333333,res[5][2],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(20833.33333333331,res[5][3],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(166666.6666666667,res[5][4],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(395833.3333333333,res[5][7],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(395833.3333333333,res[5][10],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(166666.6666666667,res[6][1],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(250000,res[6][6],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(541666.6666666667,res[6][9],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(166666.6666666667,res[6][11],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(83333.33333333331,res[7][0],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(479166.6666666667,res[7][1],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(333333.3333333333,res[7][2],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(624999.9999999997,res[7][3],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(479166.6666666667,res[7][4],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(479166.6666666667,res[7][5],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(83333.33333333333,res[7][6],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(83333.33333333331,res[7][7],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(83333.33333333333,res[7][8],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(83333.33333333333,res[7][9],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(83333.33333333331,res[7][10],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(479166.6666666667,res[7][11],1e-7);
+ //clean up
+ sourceMesh->decrRef();
+ targetMesh->decrRef();
+}
+void MEDCouplingBasicsTest::test3DInterpP0P1_1()
+{
+ MEDCouplingUMesh *sourceMesh=build3DTargetMesh_1();
+ MEDCouplingUMesh *targetMesh=build3DSourceMesh_1();
+ //
+ MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(sourceMesh);
+ MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(targetMesh);
+ INTERP_KERNEL::Interpolation3D myInterpolator;
+ vector<map<int,double> > res;
+ myInterpolator.setPrecision(1e-12);
+ myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,res,"P0P1");
+ CPPUNIT_ASSERT_EQUAL(9,(int)res.size());
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(244444.4444444445,res[0][4],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(145833.3333333333,res[0][5],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(291666.6666666666,res[0][6],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(151388.8888888889,res[0][7],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(125000,res[1][0],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(140277.7777777778,res[1][1],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(119444.4444444444,res[1][2],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(151388.8888888889,res[1][3],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(119444.4444444444,res[1][4],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(151388.8888888889,res[1][5],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(26388.88888888889,res[1][6],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(348611.1111111111,res[2][6],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(151388.8888888888,res[2][7],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(244444.4444444444,res[3][2],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(145833.3333333334,res[3][3],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(291666.6666666666,res[3][6],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(151388.8888888889,res[3][7],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(536111.111111111,res[4][5],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(297222.2222222221,res[4][7],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(223611.1111111111,res[5][1],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(125000,res[5][3],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(125000,res[5][5],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(26388.88888888892,res[5][7],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(833333.333333333,res[6][7],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(536111.1111111109,res[7][3],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(297222.2222222221,res[7][7],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(11111.1111111111,res[8][1],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(11111.11111111111,res[8][2],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(166666.6666666666,res[8][3],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(11111.11111111111,res[8][4],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(166666.6666666667,res[8][5],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(166666.6666666667,res[8][6],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1466666.666666668,res[8][7],1e-7);
+ //clean up
+ sourceMesh->decrRef();
+ targetMesh->decrRef();
+}
+
+void MEDCouplingBasicsTest::test3DInterpP1P0_1()
+{
+ MEDCouplingUMesh *sourceMesh=build3DSourceMesh_1();
+ MEDCouplingUMesh *targetMesh=build3DTargetMesh_1();
+ //
+ MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(sourceMesh);
+ MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(targetMesh);
+ INTERP_KERNEL::Interpolation3D myInterpolator;
+ vector<map<int,double> > res;
+ myInterpolator.setPrecision(1e-12);
+ myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,res,"P1P0");
+ CPPUNIT_ASSERT_EQUAL(8,(int)res.size());
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(125000,res[0][1],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(140277.7777777778,res[1][1],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(223611.1111111111,res[1][5],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(11111.1111111111,res[1][8],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(119444.4444444444,res[2][1],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(244444.4444444445,res[2][3],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(11111.11111111111,res[2][8],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(151388.8888888889,res[3][1],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(145833.3333333333,res[3][3],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(125000,res[3][5],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(536111.1111111109,res[3][7],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(166666.6666666667,res[3][8],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(244444.4444444445,res[4][0],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(119444.4444444445,res[4][1],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(11111.11111111111,res[4][8],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(145833.3333333333,res[5][0],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(151388.8888888889,res[5][1],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(536111.1111111109,res[5][4],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(125000,res[5][5],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(166666.6666666666,res[5][8],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(291666.6666666666,res[6][0],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(26388.88888888889,res[6][1],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(348611.1111111112,res[6][2],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(291666.6666666667,res[6][3],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(166666.6666666666,res[6][8],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(151388.8888888889,res[7][0],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(151388.8888888889,res[7][2],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(151388.8888888889,res[7][3],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(297222.2222222221,res[7][4],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(26388.88888888892,res[7][5],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(833333.333333333,res[7][6],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(297222.2222222222,res[7][7],1e-7);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1466666.666666668,res[7][8],1e-7);
//clean up
sourceMesh->decrRef();
+ targetMesh->decrRef();
}
MEDCouplingUMesh *MEDCouplingBasicsTest::build2DSourceMesh_1()
return targetMesh;
}
+MEDCouplingUMesh *MEDCouplingBasicsTest::build2DTargetMesh_2()
+{
+ double targetCoords[18]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 };
+ int targetConn[24]={0,3,4, 0,4,1, 1,4,2, 4,5,2, 3,6,4, 6,7,4, 4,7,5, 7,8,5 };
+ MEDCouplingUMesh *targetMesh=MEDCouplingUMesh::New();
+ targetMesh->setMeshDimension(2);
+ targetMesh->allocateCells(8);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+3);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+6);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+9);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+12);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+15);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+18);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+21);
+ targetMesh->finishInsertingCells();
+ DataArrayDouble *myCoords=DataArrayDouble::New();
+ myCoords->alloc(9,2);
+ std::copy(targetCoords,targetCoords+18,myCoords->getPointer());
+ targetMesh->setCoords(myCoords);
+ myCoords->decrRef();
+ return targetMesh;
+}
+
MEDCouplingUMesh *MEDCouplingBasicsTest::build3DSurfSourceMesh_1()
{
double sourceCoords[12]={-0.3,-0.3,0.5, 0.7,-0.3,1.5, -0.3,0.7,0.5, 0.7,0.7,1.5};
return targetMesh;
}
+MEDCouplingUMesh *MEDCouplingBasicsTest::build3DSurfTargetMesh_2()
+{
+ double targetCoords[27]={-0.3,-0.3,0.5, 0.2,-0.3,1., 0.7,-0.3,1.5, -0.3,0.2,0.5, 0.2,0.2,1., 0.7,0.2,1.5, -0.3,0.7,0.5, 0.2,0.7,1., 0.7,0.7,1.5};
+ int targetConn[24]={0,3,4, 0,4,1, 1,4,2, 4,5,2, 3,6,4, 6,7,4, 4,7,5, 7,8,5 };
+ MEDCouplingUMesh *targetMesh=MEDCouplingUMesh::New();
+ targetMesh->setMeshDimension(2);
+ targetMesh->allocateCells(8);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+3);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+6);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+9);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+12);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+15);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+18);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+21);
+ targetMesh->finishInsertingCells();
+ DataArrayDouble *myCoords=DataArrayDouble::New();
+ myCoords->alloc(9,3);
+ std::copy(targetCoords,targetCoords+27,myCoords->getPointer());
+ targetMesh->setCoords(myCoords);
+ myCoords->decrRef();
+ return targetMesh;
+}
+
MEDCouplingUMesh *MEDCouplingBasicsTest::build3DSourceMesh_1()
{
double sourceCoords[27]={ 0.0, 0.0, 200.0, 0.0, 0.0, 0.0, 0.0, 200.0, 200.0, 0.0, 200.0, 0.0, 200.0, 0.0, 200.0,
sourceMesh->finishInsertingCells();
DataArrayDouble *myCoords=DataArrayDouble::New();
myCoords->alloc(9,3);
- std::copy(sourceCoords,sourceCoords+12,myCoords->getPointer());
+ std::copy(sourceCoords,sourceCoords+27,myCoords->getPointer());
sourceMesh->setCoords(myCoords);
myCoords->decrRef();
return sourceMesh;
}
+MEDCouplingUMesh *MEDCouplingBasicsTest::build3DTargetMesh_1()
+{
+ double targetCoords[81]={ 0., 0., 0., 50., 0., 0. , 200., 0., 0. , 0., 50., 0., 50., 50., 0. , 200., 50., 0., 0., 200., 0., 50., 200., 0. , 200., 200., 0. ,
+ 0., 0., 50., 50., 0., 50. , 200., 0., 50. , 0., 50., 50., 50., 50., 50. , 200., 50., 50., 0., 200., 50., 50., 200., 50. , 200., 200., 50. ,
+ 0., 0., 200., 50., 0., 200. , 200., 0., 200. , 0., 50., 200., 50., 50., 200. , 200., 50., 200., 0., 200., 200., 50., 200., 200. , 200., 200., 200. };
+ int targetConn[64]={0,1,4,3,9,10,13,12, 1,2,5,4,10,11,14,13, 3,4,7,6,12,13,16,15, 4,5,8,7,13,14,17,16,
+ 9,10,13,12,18,19,22,21, 10,11,14,13,19,20,23,22, 12,13,16,15,21,22,25,24, 13,14,17,16,22,23,26,25};
+ MEDCouplingUMesh *targetMesh=MEDCouplingUMesh::New();
+ targetMesh->setMeshDimension(3);
+ targetMesh->allocateCells(12);
+ for(int i=0;i<8;i++)
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_HEXA8,8,targetConn+8*i);
+ targetMesh->finishInsertingCells();
+ DataArrayDouble *myCoords=DataArrayDouble::New();
+ myCoords->alloc(27,3);
+ std::copy(targetCoords,targetCoords+81,myCoords->getPointer());
+ targetMesh->setCoords(myCoords);
+ myCoords->decrRef();
+ return targetMesh;
+}
+
double MEDCouplingBasicsTest::sumAll(const std::vector< std::map<int,double> >& matrix)
{
double ret=0.;
class MEDCouplingBasicsTest : public CppUnit::TestFixture
{
CPPUNIT_TEST_SUITE(MEDCouplingBasicsTest);
+ CPPUNIT_TEST( testArray );
CPPUNIT_TEST( testMesh );
+ CPPUNIT_TEST( testMeshPointsCloud );
CPPUNIT_TEST( testDeepCopy );
CPPUNIT_TEST( testRevNodal );
CPPUNIT_TEST( testBuildPartOfMySelf );
CPPUNIT_TEST( test2DInterpP0P0_1 );
CPPUNIT_TEST( test2DInterpP0P1_1 );
CPPUNIT_TEST( test2DInterpP1P0_1 );
+ CPPUNIT_TEST( test2DInterpP1P1_1 );
CPPUNIT_TEST( test3DSurfInterpP0P0_1 );
CPPUNIT_TEST( test3DSurfInterpP0P1_1 );
CPPUNIT_TEST( test3DSurfInterpP1P0_1 );
CPPUNIT_TEST( test3DInterpP0P0_1 );
+ CPPUNIT_TEST( test3DInterpP0P1_1 );
+ CPPUNIT_TEST( test3DInterpP1P0_1 );
CPPUNIT_TEST_SUITE_END();
public:
+ void testArray();
void testMesh();
+ void testMeshPointsCloud();
void testDeepCopy();
void testRevNodal();
void testBuildPartOfMySelf();
void test2DInterpP0P0_1();
void test2DInterpP0P1_1();
void test2DInterpP1P0_1();
+ void test2DInterpP1P1_1();
void test3DSurfInterpP0P0_1();
void test3DSurfInterpP0P1_1();
void test3DSurfInterpP1P0_1();
void test3DInterpP0P0_1();
+ void test3DInterpP0P1_1();
+ void test3DInterpP1P0_1();
private:
MEDCouplingUMesh *build2DSourceMesh_1();
MEDCouplingUMesh *build2DTargetMesh_1();
+ MEDCouplingUMesh *build2DTargetMesh_2();
MEDCouplingUMesh *build3DSurfSourceMesh_1();
MEDCouplingUMesh *build3DSurfTargetMesh_1();
+ MEDCouplingUMesh *build3DSurfTargetMesh_2();
MEDCouplingUMesh *build3DSourceMesh_1();
+ MEDCouplingUMesh *build3DTargetMesh_1();
double sumAll(const std::vector< std::map<int,double> >& matrix);
};
}
//
#include "BlockTopology.hxx"
#include "MemArray.hxx"
-#include "MEDCouplingSMesh.hxx"
+#include "MEDCouplingRMesh.hxx"
#include "CommInterface.hxx"
#include "ProcessorGroup.hxx"
#include "MPIProcessorGroup.hxx"
* instead of making the best choice with respect to the
* values of the different axes.
*/
- BlockTopology::BlockTopology(const ProcessorGroup& group, MEDCouplingSMesh *grid):
+ BlockTopology::BlockTopology(const ProcessorGroup& group, MEDCouplingRMesh *grid):
_proc_group(&group), _dimension(grid->getSpaceDimension()), _owns_processor_group(false)
{
vector <int> axis_length(_dimension);
namespace ParaMEDMEM
{
class ComponentTopology;
- class MEDCouplingSMesh;
+ class MEDCouplingRMesh;
typedef enum{Block,Cycle} CYCLE_TYPE;
{
public:
BlockTopology() { }
- BlockTopology(const ProcessorGroup& group, MEDCouplingSMesh *grid);
+ BlockTopology(const ProcessorGroup& group, MEDCouplingRMesh *grid);
BlockTopology(const BlockTopology& geom_topo, const ComponentTopology& comp_topo);
BlockTopology(const ProcessorGroup& group, int nb_elem);
virtual ~BlockTopology();
localgroup=_source_group;
else
localgroup=_target_group;
+ //delete _icoco_field;
_icoco_field=new ICoCo::MEDField(*const_cast<ICoCo::TrioField* >(triofield), *localgroup);
attachLocalField(_icoco_field);
return;
// between the distant ids and the ids of the local reconstruction
// ==========================================================================
void ElementLocator::exchangeMesh(int idistantrank,
- MEDCouplingUMesh*& distant_mesh,
+ MEDCouplingPointSet*& distant_mesh,
int*& distant_ids)
{
int dim = _local_cell_mesh->getSpaceDimension();
return;
}
- set <int> elems;
+ vector<int> elems;
double* distant_bb = _domain_bounding_boxes+rank*2*dim;
- double* elem_bb=new double[2*dim];
-
- //defining pointers to med
- const int* conn = _local_cell_mesh->getNodalConnectivity()->getPointer() ;
- const int* conn_index= _local_cell_mesh->getNodalConnectivityIndex()->getPointer();
- const double* coords = _local_cell_mesh->getCoords()->getPointer() ;
-
- for ( int ielem=0; ielem<_local_cell_mesh->getNumberOfCells() ; ielem++)
- {
- for (int i=0; i<dim; i++)
- {
- elem_bb[i*2]=std::numeric_limits<double>::max();
- elem_bb[i*2+1]=-std::numeric_limits<double>::max();
- }
-
- for (int inode=conn_index[ielem]+1; inode<conn_index[ielem+1]; inode++)//+1 due to offset of cell type.
- {
- int node= conn[inode];
-
- for (int idim=0; idim<dim; idim++)
- {
- if ( coords[node*dim+idim] < elem_bb[idim*2] )
- {
- elem_bb[idim*2] = coords[node*dim+idim] ;
- }
- if ( coords[node*dim+idim] > elem_bb[idim*2+1] )
- {
- elem_bb[idim*2+1] = coords[node*dim+idim] ;
- }
- }
- }
- if (_intersectsBoundingBox(elem_bb, distant_bb, dim))
- {
- elems.insert(ielem);
- }
- }
+ _local_cell_mesh->giveElemsInBoundingBox(distant_bb,getBoundingBoxAdjustment(),elems);
//send_mesh contains null pointer if elems is empty
- MEDCouplingUMesh* send_mesh = _meshFromElems(elems);
-
+ MEDCouplingPointSet* send_mesh = _local_cell_mesh->buildPartOfMySelf(&elems[0],&elems[elems.size()],false);
// Constituting an array containing the ids of the elements that are
// going to be sent to the distant subdomain.
// This array enables the correct redistribution of the data when the
{
distant_ids_send = new int[elems.size()];
int index=0;
- for (std::set<int>::const_iterator iter = elems.begin(); iter!= elems.end(); iter++)
+ for (std::vector<int>::const_iterator iter = elems.begin(); iter!= elems.end(); iter++)
{
distant_ids_send[index]=*iter;
index++;
}
_exchangeMesh(send_mesh, distant_mesh, idistantrank, distant_ids_send, distant_ids);
delete[] distant_ids_send;
- delete[] elem_bb;
if(send_mesh)
send_mesh->decrRef();
}
CommInterface comm_interface =_union_group->getCommInterface();
int dim = _local_cell_mesh->getSpaceDimension();
_domain_bounding_boxes = new double[2*dim*_union_group->size()];
- const double* coords = _local_cell_mesh->getCoords()->getPointer() ;
-
- int nbnodes = _local_cell_mesh->getNumberOfNodes();
double * minmax=new double [2*dim];
- for (int idim=0; idim<dim; idim++)
- {
- minmax[idim*2]=std::numeric_limits<double>::max();
- minmax[idim*2+1]=-std::numeric_limits<double>::max();
- }
-
- for (int i=0; i<nbnodes; i++)
- {
- for (int idim=0; idim<dim;idim++)
- {
- if ( minmax[idim*2] > coords[i*dim+idim] )
- {
- minmax[idim*2] = coords[i*dim+idim] ;
- }
- if ( minmax[idim*2+1] < coords[i*dim+idim] )
- {
- minmax[idim*2+1] = coords[i*dim+idim] ;
- }
- }
- }
+ _local_cell_mesh->getBoundingBox(minmax);
MPIProcessorGroup* group=static_cast<MPIProcessorGroup*> (_union_group);
const MPI_Comm* comm = group->getComm();
return true;
}
- // =============================================
- // Intersect Bounding Box given 2 Bounding Boxes
- // =============================================
- bool ElementLocator::_intersectsBoundingBox(double* bb1, double* bb2, int dim)
- {
- double bbtemp[2*dim];
- double deltamax=0.0;
- double adjustment_eps=getBoundingBoxAdjustment();
-
- for (int i=0; i< dim; i++)
- {
- double delta = bb1[2*i+1]-bb1[2*i];
- if ( delta > deltamax )
- {
- deltamax = delta ;
- }
- // deltamax = (delta>deltamax)?delta:deltamax;
- }
- for (int i=0; i<dim; i++)
- {
- bbtemp[i*2]=bb1[i*2]-deltamax*adjustment_eps;
- bbtemp[i*2+1]=bb1[i*2+1]+deltamax*adjustment_eps;
- }
-
- for (int idim=0; idim < dim; idim++)
- {
- bool intersects = (bbtemp[idim*2]<bb2[idim*2+1])
- && (bb2[idim*2]<bbtemp[idim*2+1]) ;
- if (!intersects) return false;
- }
- return true;
- }
-
-
// ======================
// Exchanging meshes data
// ======================
- void ElementLocator::_exchangeMesh( MEDCouplingUMesh* local_mesh,
- MEDCouplingUMesh*& distant_mesh,
+ void ElementLocator::_exchangeMesh( MEDCouplingPointSet* local_mesh,
+ MEDCouplingPointSet*& distant_mesh,
int iproc_distant,
const int* distant_ids_send,
int*& distant_ids_recv)
// First stage : exchanging sizes
// ------------------------------
-
- int* send_buffer = new int[5];
- int* recv_buffer = new int[5];
-
- //treatment for non-empty mesh
- int nbconn=0;
- int nbelems=0;
-
- if (local_mesh !=0)
- {
- nbelems = local_mesh->getNumberOfCells();
- nbconn = local_mesh->getMeshLength();
- send_buffer[0] = local_mesh->getSpaceDimension();
- send_buffer[1] = local_mesh->getMeshDimension();
- send_buffer[2] = local_mesh->getNumberOfNodes();
- send_buffer[3] = nbelems;
- send_buffer[4] = nbconn;
- }
- else
- {
- for (int i=0; i<5; i++)
- {
- send_buffer[i]=0;
- }
- }
-
+ vector<int> tinyInfoLocal,tinyInfoDistant;
+ local_mesh->getTinySerializationInformation(tinyInfoLocal);
+ tinyInfoLocal.push_back(local_mesh->getNumberOfCells());
+ tinyInfoDistant.resize(tinyInfoLocal.size());
+ std::fill(tinyInfoDistant.begin(),tinyInfoDistant.end(),0);
MPIProcessorGroup* group=static_cast<MPIProcessorGroup*> (_union_group);
- const MPI_Comm* comm=(group->getComm());
+ const MPI_Comm* comm=group->getComm();
MPI_Status status;
-
+
// iproc_distant is the number of proc in distant group
// it must be converted to union numbering before communication
int iprocdistant_in_union = group->translateRank(&_distant_group,
iproc_distant);
-
- comm_interface.sendRecv(send_buffer, 5, MPI_INT, iprocdistant_in_union, 1112,
- recv_buffer, 5, MPI_INT,iprocdistant_in_union,1112,
+
+ comm_interface.sendRecv(&tinyInfoLocal[0], tinyInfoLocal.size(), MPI_INT, iprocdistant_in_union, 1112,
+ &tinyInfoDistant[0], tinyInfoDistant.size(), MPI_INT,iprocdistant_in_union,1112,
*comm, &status);
-
- int distant_space_dim = recv_buffer[0];
- int distant_mesh_dim = recv_buffer[1];
- int distant_nnodes = recv_buffer[2];
- int distant_nb_elems = recv_buffer[3];
- int distant_nb_conn = recv_buffer[4];
-
- delete[] send_buffer;
- delete[] recv_buffer;
-
- // Second stage : exchanging connectivity buffers
- // ----------------------------------------------
-
- int nb_integers = nbconn + 2*nbelems + 1;
- send_buffer = new int[nb_integers];
- const int* conn = 0;
- const int* global_numbering=0;
- int* ptr_buffer = send_buffer;
-
- if (local_mesh != 0)
- {
- conn = local_mesh->getNodalConnectivity()->getPointer();
-
- global_numbering = local_mesh->getNodalConnectivityIndex()->getPointer() ;
-
- //copying the data in the integer buffer
-
- memcpy(ptr_buffer, global_numbering, (nbelems+1)*sizeof(int));
- ptr_buffer += nbelems+1;
- memcpy(ptr_buffer,conn, nbconn*sizeof(int));
- ptr_buffer += nbconn;
- memcpy(ptr_buffer, distant_ids_send, nbelems*sizeof(int));
- }
-
- // Preparing the recv buffers
- int nb_recv_integers = distant_nb_conn + 2*distant_nb_elems + 1 ;
- recv_buffer=new int[nb_recv_integers];
-
- // Exchanging integer buffer
- comm_interface.sendRecv(send_buffer, nb_integers, MPI_INT,
+ DataArrayInt *v1Local=0;
+ DataArrayDouble *v2Local=0;
+ DataArrayInt *v1Distant=DataArrayInt::New();
+ DataArrayDouble *v2Distant=DataArrayDouble::New();
+ local_mesh->serialize(v1Local,v2Local);
+ local_mesh->resizeForSerialization(tinyInfoDistant,v1Distant,v2Distant);
+ comm_interface.sendRecv(v1Local->getPointer(), v1Local->getNbOfElems(), MPI_INT,
iprocdistant_in_union, 1111,
- recv_buffer, nb_recv_integers, MPI_INT,
+ v1Distant->getPointer(), v1Distant->getNbOfElems(), MPI_INT,
iprocdistant_in_union,1111,
*comm, &status);
-
- if ( nb_integers>0 )
- {
- delete[] send_buffer;
- }
-
- // Third stage : exchanging coordinates
- // ------------------------------------
-
- int nb_recv_floats = distant_space_dim*distant_nnodes;
- int nb_send_floats = 0;
- double* coords=0;
-
- if ( local_mesh!=0 )
- {
- nb_send_floats = local_mesh->getSpaceDimension()
- * local_mesh->getNumberOfNodes();
- coords = local_mesh->getCoords()->getPointer();
- }
-
- DataArrayDouble* myCoords=DataArrayDouble::New();
- myCoords->alloc(distant_nnodes,distant_space_dim);
-
- comm_interface.sendRecv(coords, nb_send_floats, MPI_DOUBLE,
+ comm_interface.sendRecv(v2Local->getPointer(), v2Local->getNbOfElems(), MPI_DOUBLE,
iprocdistant_in_union, 1112,
- myCoords->getPointer(), nb_recv_floats, MPI_DOUBLE,
+ v2Distant->getPointer(), v2Distant->getNbOfElems(), MPI_DOUBLE,
iprocdistant_in_union, 1112,
- *group->getComm(), &status);
-
-
- // Reconstructing an image of the distant mesh locally
-
- if ( nb_recv_integers>0 && distant_space_dim !=0 )
- {
- MEDCouplingUMesh* meshing = MEDCouplingUMesh::New() ;
-
- // Coordinates
- meshing->setCoords(myCoords) ;
- myCoords->decrRef();
- // Connectivity
-
- int *work=recv_buffer;
- DataArrayInt* myConnecIndex=DataArrayInt::New();
- myConnecIndex->alloc(distant_nb_elems+1,1);
- memcpy(myConnecIndex->getPointer(), work, (distant_nb_elems+1)*sizeof(int));
- work += distant_nb_elems + 1 ;
-
- DataArrayInt* myConnec=DataArrayInt::New();
- myConnec->alloc(distant_nb_conn,1);
- memcpy(myConnec->getPointer(), work, (distant_nb_conn)*sizeof(int));
- work+=distant_nb_conn;
- meshing->setConnectivity(myConnec, myConnecIndex) ;
- myConnec->decrRef();
- myConnecIndex->decrRef();
-
- // correspondence between the distant ids and the ids of
- // the local reconstruction
-
- distant_ids_recv=new int [distant_nb_elems];
- for (int i=0; i<distant_nb_elems; i++)
- {
- distant_ids_recv[i]=*work++;
- }
-
- // Mesh dimension
- meshing->setMeshDimension(distant_mesh_dim);
-
- distant_mesh=meshing;
- delete[] recv_buffer;
- }
-
- }
-
-
- // ==============
- // _meshFromElems
- // ==============
-
- MEDCouplingUMesh* ElementLocator::_meshFromElems(set<int>& elems)
- {
- //returns null pointer if there are no elems in the mesh
- if ( elems.size()==0 ) return 0;
-
- // Defining pointers
- const int* conn_mesh =
- const_cast<int*> (_local_cell_mesh->getNodalConnectivity()->getPointer());
-
- const int* conn_index =
- const_cast<int*> (_local_cell_mesh->getNodalConnectivityIndex()->getPointer());
-
- const double* coords =
- const_cast<double*> ( _local_cell_mesh->getCoords()->getPointer());
-
- set<int> nodes;
- int nbconn=0;
- for (set<int>::const_iterator iter=elems.begin(); iter!=elems.end(); iter++)
- {
- // Conn_index : C-like Addresses
- for (int inode=conn_index[*iter]+1; inode<conn_index[*iter+1]; inode++)
- {
- nodes.insert(conn_mesh[inode]);
- nbconn++ ;
- }
- }
-
- map<int,int> big2small;
- int i=0;
- for (set<int>::const_iterator iter=nodes.begin(); iter!=nodes.end(); iter++)
- {
- big2small[*iter]=i;
- i++;
- }
-
- // Memory allocate
- DataArrayInt *conn = DataArrayInt::New() ;
- conn->alloc(nbconn+elems.size(),1) ;
- int *connPtr=conn->getPointer();
-
- DataArrayInt * connIndex = DataArrayInt::New() ;
- connIndex->alloc(elems.size()+1,1) ;
- int* connIndexPtr=connIndex->getPointer();
-
- DataArrayDouble *new_coords = DataArrayDouble::New() ;
- new_coords->alloc(nodes.size(), _local_cell_mesh->getSpaceDimension()) ;
- double *new_coords_ptr = new_coords->getPointer();
-
- // New connectivity table
- int index=0;
- int mainIndex=0;
- for (set<int>::const_iterator iter=elems.begin(); iter!=elems.end(); iter++,mainIndex++)
- {
- connIndexPtr[mainIndex]=index;
- connPtr[index++]=conn_mesh[conn_index[*iter]];
- for (int inode = conn_index[*iter]+1; inode < conn_index[*iter+1]; inode++)
- {
- connPtr[index]=big2small[conn_mesh[inode]] ; // C-like number
- index++;
- }
- }
- connIndexPtr[mainIndex]=index;
- // Coordinates
- index=0;
- for (set<int>::const_iterator iter=nodes.begin(); iter!=nodes.end(); iter++)
+ *comm, &status);
+ if(v1Distant->getNbOfElems()>0)
{
- int dim = _local_cell_mesh->getSpaceDimension();
- for (int i=0; i<dim;i++)
- {
- new_coords_ptr[index]=coords[(*iter)*dim+i];
- index++;
- }
+ distant_mesh=local_mesh->buildObjectFromUnserialization(tinyInfoDistant,v1Distant,v2Distant);
}
-
- // Initialize
- MEDCouplingUMesh* meshing = MEDCouplingUMesh::New() ;
- meshing->setCoords(new_coords) ;
- new_coords->decrRef();
- meshing->setConnectivity(conn, connIndex) ;
- conn->decrRef();
- connIndex->decrRef();
- meshing->setMeshDimension(_local_cell_mesh->getMeshDimension());
-
- return meshing;
+ distant_ids_recv=new int[tinyInfoDistant.back()];
+ comm_interface.sendRecv((void *)distant_ids_send,tinyInfoLocal.back(), MPI_INT,
+ iprocdistant_in_union, 1113,
+ distant_ids_recv,tinyInfoDistant.back(), MPI_INT,
+ iprocdistant_in_union,1113,
+ *comm, &status);
+ v1Local->decrRef();
+ v2Local->decrRef();
+ v1Distant->decrRef();
+ v2Distant->decrRef();
}
}
#define __ELEMENTLOCATOR_HXX__
#include "InterpolationOptions.hxx"
-#include "MEDCouplingUMesh.hxx"
#include <vector>
#include <set>
class ProcessorGroup;
class ParaSUPPORT;
class InterpolationMatrix;
-
+ class MEDCouplingPointSet;
class ElementLocator : public INTERP_KERNEL::InterpolationOptions
{
virtual ~ElementLocator();
void exchangeMesh(int idistantrank,
- MEDCouplingUMesh*& target_mesh,
+ MEDCouplingPointSet*& target_mesh,
int*& distant_ids);
void exchangeMethod(const std::string& sourceMeth, int idistantrank, std::string& targetMeth);
private:
const ParaMESH& _local_para_mesh ;
- MEDCouplingUMesh* _local_cell_mesh;
- MEDCouplingUMesh* _local_face_mesh;
- std::vector<MEDCouplingUMesh*> _distant_cell_meshes;
- std::vector<MEDCouplingUMesh*> _distant_face_meshes;
+ MEDCouplingPointSet* _local_cell_mesh;
+ MEDCouplingPointSet* _local_face_mesh;
+ std::vector<MEDCouplingPointSet*> _distant_cell_meshes;
+ std::vector<MEDCouplingPointSet*> _distant_face_meshes;
double* _domain_bounding_boxes;
const ProcessorGroup& _distant_group;
const ProcessorGroup& _local_group;
void _computeBoundingBoxes();
bool _intersectsBoundingBox(int irank);
- bool _intersectsBoundingBox(double* bb1, double* bb2, int dim);
- void _exchangeMesh(MEDCouplingUMesh* local_mesh, MEDCouplingUMesh*& distant_mesh,
+ void _exchangeMesh(MEDCouplingPointSet* local_mesh, MEDCouplingPointSet*& distant_mesh,
int iproc_distant, const int* distant_ids_send,
int*& distant_ids_recv);
- MEDCouplingUMesh* _meshFromElems(std::set<int>& elems);
};
}
//field on the sending end
if (triofield._field!=0)
{
- _field = new ParaMEDMEM::ParaFIELD(ParaMEDMEM::ON_CELLS,_mesh, *_comp_topology );
+ _field = new ParaMEDMEM::ParaFIELD(ParaMEDMEM::ON_CELLS,ParaMEDMEM::ONE_TIME,_mesh, *_comp_topology );
ParaMEDMEM::DataArrayDouble *fieldArr=_field->getField()->getArray();
_field->getField()->setName(triofield.getName().c_str());
- _field->getField()->setTime(triofield._time1);
- _field->getField()->setDtIt(0,triofield._itnumber);
+ _field->getField()->setTime(triofield._time1,0,triofield._itnumber);
for (int i =0; i<triofield._nb_elems; i++)
for (int j=0; j<triofield._nb_field_components; j++)
{
else
{
// _field = new ParaMEDMEM::ParaFIELD(ParaMEDMEM::ON_CELLS,_support, *_comp_topology );
- _field = new ParaMEDMEM::ParaFIELD(ParaMEDMEM::ON_CELLS,_mesh, *_comp_topology );
+ _field = new ParaMEDMEM::ParaFIELD(ParaMEDMEM::ON_CELLS,ParaMEDMEM::ONE_TIME,_mesh, *_comp_topology );
_field->getField()->setName(triofield.getName().c_str());
- _field->getField()->setTime(triofield._time1);
- _field->getField()->setDtIt(0,triofield._itnumber);
+ _field->getField()->setTime(triofield._time1,0,triofield._itnumber);
// the trio field points to the pointer inside the MED field
triofield._field=const_cast<double*> (_field->getField()->getArray()->getPointer());
for (int i=0; i<triofield._nb_field_components*triofield._nb_elems;i++)
#include "Interpolation2D.txx"
#include "Interpolation3DSurf.txx"
#include "Interpolation3D.txx"
+#include "MEDCouplingUMesh.hxx"
#include "MEDCouplingNormalizedUnstructuredMesh.txx"
#include "InterpolationOptions.hxx"
-#include "VolSurfFormulae.hxx"
#include "NormalizedUnstructuredMesh.hxx"
// class InterpolationMatrix
// the subdomain and the actual elem ids on the distant subdomain
// ======================================================================
- void InterpolationMatrix::addContribution ( MEDCouplingUMesh& distant_support,
+ void InterpolationMatrix::addContribution ( MEDCouplingPointSet& distant_support,
int iproc_distant,
int* distant_elems,
const std::string& srcMeth,
vector<map<int,double> > surfaces;
int colSize=0;
//computation of the intersection volumes between source and target elements
-
+ MEDCouplingUMesh *distant_supportC=dynamic_cast<MEDCouplingUMesh *>(&distant_support);
+ MEDCouplingUMesh *source_supportC=dynamic_cast<MEDCouplingUMesh *>(_source_support);
if ( distant_support.getMeshDimension() == 2
&& distant_support.getSpaceDimension() == 3 )
{
- MEDCouplingNormalizedUnstructuredMesh<3,2> target_wrapper(&distant_support);
- MEDCouplingNormalizedUnstructuredMesh<3,2> source_wrapper(_source_support);
+ MEDCouplingNormalizedUnstructuredMesh<3,2> target_wrapper(distant_supportC);
+ MEDCouplingNormalizedUnstructuredMesh<3,2> source_wrapper(source_supportC);
INTERP_KERNEL::Interpolation3DSurf interpolator (*this);
colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod.c_str());
else if ( distant_support.getMeshDimension() == 2
&& distant_support.getSpaceDimension() == 2)
{
- MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(&distant_support);
- MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(_source_support);
+ MEDCouplingNormalizedUnstructuredMesh<2,2> target_wrapper(distant_supportC);
+ MEDCouplingNormalizedUnstructuredMesh<2,2> source_wrapper(source_supportC);
INTERP_KERNEL::Interpolation2D interpolator (*this);
colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod.c_str());
else if ( distant_support.getMeshDimension() == 3
&& distant_support.getSpaceDimension() ==3 )
{
- MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(&distant_support);
- MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(_source_support);
+ MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(distant_supportC);
+ MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(source_supportC);
INTERP_KERNEL::Interpolation3D interpolator (*this);
colSize=interpolator.interpolateMeshes(target_wrapper,source_wrapper,surfaces,interpMethod.c_str());
int source_size=surfaces.size();
- MEDCouplingFieldDouble *target_triangle_surf =
- getSupportVolumes(&distant_support);
- MEDCouplingFieldDouble *source_triangle_surf =
- getSupportVolumes(_source_support) ;
+ MEDCouplingFieldDouble *target_triangle_surf = distant_support.getMeasureField();
+ MEDCouplingFieldDouble *source_triangle_surf = _source_support->getMeasureField();
// Storing the source volumes
_source_volume.resize(source_size);
}
}
-
- MEDCouplingFieldDouble* InterpolationMatrix::getSupportVolumes(MEDCouplingMesh * mesh)
- {
- if(!mesh->isStructured())
- return getSupportUnstructuredVolumes((MEDCouplingUMesh *)mesh);
- else
- throw INTERP_KERNEL::Exception("Not implemented yet !!!");
- }
-
- // ====================================================================
- // brief returns the volumes of the cells underlying the field \a field
-
- // For 2D geometries, the returned field contains the areas.
- // For 3D geometries, the returned field contains the volumes.
-
- // param field field on which cells the volumes are required
- // return field containing the volumes
- // ====================================================================
-
- MEDCouplingFieldDouble* InterpolationMatrix::getSupportUnstructuredVolumes(MEDCouplingUMesh * mesh)
- {
- int ipt, type ;
- int nbelem = mesh->getNumberOfCells() ;
- int dim_mesh = mesh->getMeshDimension();
- int dim_space = mesh->getSpaceDimension() ;
- double *coords = mesh->getCoords()->getPointer() ;
- int *connec = mesh->getNodalConnectivity()->getPointer() ;
- int *connec_index = mesh->getNodalConnectivityIndex()->getPointer() ;
-
-
- MEDCouplingFieldDouble* field = MEDCouplingFieldDouble::New(ON_CELLS);
- DataArrayDouble* array = DataArrayDouble::New() ;
- array->alloc(nbelem, 1) ;
- double *area_vol = array->getPointer() ;
-
- switch (dim_mesh)
- {
- case 2: // getting the areas
- for ( int iel=0 ; iel<nbelem ; iel++ )
- {
- ipt = connec_index[iel] ;
- type = connec[ipt] ;
-
- switch ( type )
- {
- case INTERP_KERNEL::NORM_TRI3 :
- case INTERP_KERNEL::NORM_TRI6 :
- {
- int N1 = connec[ipt+1];
- int N2 = connec[ipt+2];
- int N3 = connec[ipt+3];
-
- area_vol[iel]=INTERP_KERNEL::calculateAreaForTria(coords+(dim_space*N1),
- coords+(dim_space*N2),
- coords+(dim_space*N3),
- dim_space);
- }
- break ;
-
- case INTERP_KERNEL::NORM_QUAD4 :
- case INTERP_KERNEL::NORM_QUAD8 :
- {
- int N1 = connec[ipt+1];
- int N2 = connec[ipt+2];
- int N3 = connec[ipt+3];
- int N4 = connec[ipt+4];
-
- area_vol[iel]=INTERP_KERNEL::calculateAreaForQuad(coords+dim_space*N1,
- coords+dim_space*N2,
- coords+dim_space*N3,
- coords+dim_space*N4,
- dim_space) ;
- }
- break ;
-
- case INTERP_KERNEL::NORM_POLYGON :
- {
- // We must remember that the first item is the type. That's
- // why we substract 1 to get the number of nodes of this polygon
- int size = connec_index[iel+1] - connec_index[iel] - 1 ;
-
- double **pts = new double *[size] ;
-
- for ( int inod=0 ; inod<size ; inod++ )
- {
- // Remember the first item is the type
- pts[inod] = coords+dim_space*connec[ipt+inod+1] ;
- }
-
- area_vol[iel]=INTERP_KERNEL::calculateAreaForPolyg((const double **)pts,
- size,dim_space);
- delete [] pts;
- }
- break ;
-
- default :
- throw INTERP_KERNEL::Exception("Bad Support to get Areas on it !");
-
- } // End of switch
-
- } // End of the loop over the cells
- break;
- case 3: // getting the volumes
- for ( int iel=0 ; iel<nbelem ; iel++ )
- {
- ipt = connec_index[iel] ;
- type = connec[ipt] ;
-
- switch ( type )
- {
- case INTERP_KERNEL::NORM_TETRA4 :
- case INTERP_KERNEL::NORM_TETRA10 :
- {
- int N1 = connec[ipt+1];
- int N2 = connec[ipt+2];
- int N3 = connec[ipt+3];
- int N4 = connec[ipt+4];
-
- area_vol[iel]=INTERP_KERNEL::calculateVolumeForTetra(coords+dim_space*N1,
- coords+dim_space*N2,
- coords+dim_space*N3,
- coords+dim_space*N4) ;
- }
- break ;
-
- case INTERP_KERNEL::NORM_PYRA5 :
- case INTERP_KERNEL::NORM_PYRA13 :
- {
- int N1 = connec[ipt+1];
- int N2 = connec[ipt+2];
- int N3 = connec[ipt+3];
- int N4 = connec[ipt+4];
- int N5 = connec[ipt+5];
-
- area_vol[iel]=INTERP_KERNEL::calculateVolumeForPyra(coords+dim_space*N1,
- coords+dim_space*N2,
- coords+dim_space*N3,
- coords+dim_space*N4,
- coords+dim_space*N5) ;
- }
- break ;
-
- case INTERP_KERNEL::NORM_PENTA6 :
- case INTERP_KERNEL::NORM_PENTA15 :
- {
- int N1 = connec[ipt+1];
- int N2 = connec[ipt+2];
- int N3 = connec[ipt+3];
- int N4 = connec[ipt+4];
- int N5 = connec[ipt+5];
- int N6 = connec[ipt+6];
-
- area_vol[iel]=INTERP_KERNEL::calculateVolumeForPenta(coords+dim_space*N1,
- coords+dim_space*N2,
- coords+dim_space*N3,
- coords+dim_space*N4,
- coords+dim_space*N5,
- coords+dim_space*N6) ;
- }
- break ;
-
- case INTERP_KERNEL::NORM_HEXA8 :
- case INTERP_KERNEL::NORM_HEXA20 :
- {
- int N1 = connec[ipt+1];
- int N2 = connec[ipt+2];
- int N3 = connec[ipt+3];
- int N4 = connec[ipt+4];
- int N5 = connec[ipt+5];
- int N6 = connec[ipt+6];
- int N7 = connec[ipt+7];
- int N8 = connec[ipt+8];
-
- area_vol[iel]=INTERP_KERNEL::calculateVolumeForHexa(coords+dim_space*N1,
- coords+dim_space*N2,
- coords+dim_space*N3,
- coords+dim_space*N4,
- coords+dim_space*N5,
- coords+dim_space*N6,
- coords+dim_space*N7,
- coords+dim_space*N8) ;
- }
- break ;
-
- case INTERP_KERNEL::NORM_POLYHED :
- {
- throw INTERP_KERNEL::Exception("Not yet implemented !");
- }
- break ;
-
- default:
- throw INTERP_KERNEL::Exception("Bad Support to get Volume on it !");
- }
- }
- break;
- default:
- throw INTERP_KERNEL::Exception("interpolation is not available for this dimension");
- }
-
- field->setArray(array) ;
- array->decrRef();
- field->setMesh(mesh) ;
-
- return field ;
- }
}
virtual ~InterpolationMatrix();
- void addContribution(MEDCouplingUMesh& distant_support, int iproc_distant,
+ void addContribution(MEDCouplingPointSet& distant_support, int iproc_distant,
int* distant_elems, const std::string& srcMeth, const std::string& targetMeth);
void multiply(MEDCouplingFieldDouble& field) const;
void transposeMultiply(MEDCouplingFieldDouble& field)const;
void prepare();
- int getNbRows() const {return _row_offsets.size();}
- MPIAccessDEC* getAccessDEC(){return _mapping.getAccessDEC();}
-
- static MEDCouplingFieldDouble *getSupportVolumes(MEDCouplingMesh *field);
- static MEDCouplingFieldDouble *getSupportUnstructuredVolumes(MEDCouplingUMesh *field);
+ int getNbRows() const { return _row_offsets.size(); }
+ MPIAccessDEC* getAccessDEC() { return _mapping.getAccessDEC(); }
private:
std::vector<int> _row_offsets;
std::map<std::pair<int,int>, int > _col_offsets;
- MEDCouplingUMesh *_source_support;
+ MEDCouplingPointSet *_source_support;
MxN_Mapping _mapping;
const ProcessorGroup& _source_group;
//transfering option from IntersectionDEC to ElementLocator
locator.setBoundingBoxAdjustment(getBoundingBoxAdjustment());
- MEDCouplingUMesh* distant_mesh=0;
+ MEDCouplingPointSet* distant_mesh=0;
int* distant_ids=0;
for (int i=0; i<_target_group->size(); i++)
{
//transfering option from IntersectionDEC to ElementLocator
locator.setBoundingBoxAdjustment(getBoundingBoxAdjustment());
- MEDCouplingUMesh* distant_mesh=0;
+ MEDCouplingPointSet* distant_mesh=0;
int* distant_ids=0;
for (int i=0; i<_source_group->size(); i++)
{
ParaMEDMEM::MEDCouplingUMesh *mesh=readUMeshFromFileLev1(fileName,meshName,meshDimRelToMax,familiesToKeep,typesToKeep,meshDim);
if(typeOfOutField==ON_CELLS)
keepSpecifiedMeshDim<MEDFieldDoublePerCellType>(fieldPerCellType,meshDim);
- ParaMEDMEM::MEDCouplingFieldDouble *ret=ParaMEDMEM::MEDCouplingFieldDouble::New(typeOfOutField);
- ret->setDtIt(iteration,order);
+ ParaMEDMEM::MEDCouplingFieldDouble *ret=ParaMEDMEM::MEDCouplingFieldDouble::New(typeOfOutField,ONE_TIME);
ret->setName(fieldName);
- ret->setTime(time);
+ ret->setTime(time,iteration,order);
ret->setMesh(mesh);
mesh->decrRef();
ParaMEDMEM::DataArrayDouble *arr=buildArrayFromRawData(fieldPerCellType);
}
if(myRank==0)
writeMasterFile(fileName,fileNames,mesh->getCellMesh()->getName());
- writeUMesh(fileNames[myRank].c_str(),mesh->getCellMesh());
+ writeUMesh(fileNames[myRank].c_str(),dynamic_cast<MEDCouplingUMesh *>(mesh->getCellMesh()));
}
void MEDLoader::writeParaField(const char *fileName, const char *meshName, ParaMEDMEM::ParaFIELD *f)
\endverbatim
*/
- ParaFIELD::ParaFIELD(TypeOfField type, ParaMESH* para_support, const ComponentTopology& component_topology)
+ ParaFIELD::ParaFIELD(TypeOfField type, TypeOfTimeDiscretization td, ParaMESH* para_support, const ComponentTopology& component_topology)
:_support(para_support),
_component_topology(component_topology),_topology(0),
_field(0),
int nb_components = component_topology.nbLocalComponents();
if (nb_components!=0)
{
- _field=MEDCouplingFieldDouble::New(type);
+ _field=MEDCouplingFieldDouble::New(type,td);
_field->setMesh(_support->getCellMesh());
DataArrayDouble *array=DataArrayDouble::New();
array->alloc(_field->getNumberOfTuples(),nb_components);
_field->setName("Default ParaFIELD name");
_field->setDescription("Default ParaFIELD description");
- _field->setDtIt(0,0);
- _field->setTime(0.0);
}
/*! \brief Constructor creating the ParaFIELD
double ParaFIELD::getVolumeIntegral(int icomp) const
{
CommInterface comm_interface = _topology->getProcGroup()->getCommInterface();
- MEDCouplingFieldDouble *volume=InterpolationMatrix::getSupportVolumes(_field->getMesh());
- double *ptr=volume->getArray()->getPointer();
- int nbOfValues=volume->getArray()->getNbOfElems();
- double integral=0.;
- for (int i=0; i<nbOfValues; i++)
- integral+=_field->getIJ(i,icomp)*ptr[i];
-
- volume->decrRef();
-
+ double integral=_field->measureAccumulate(icomp);
double total=0.;
const MPI_Comm* comm = (dynamic_cast<const MPIProcessorGroup*>(_topology->getProcGroup()))->getComm();
comm_interface.allReduce(&integral, &total, 1, MPI_DOUBLE, MPI_SUM, *comm);
return total;
}
}
-
{
public:
- ParaFIELD(TypeOfField type, ParaMESH* mesh, const ComponentTopology& component_topology);
+ ParaFIELD(TypeOfField type, TypeOfTimeDiscretization td, ParaMESH* mesh, const ComponentTopology& component_topology);
ParaFIELD(MEDCouplingFieldDouble* field, const ProcessorGroup& group);
#include "Topology.hxx"
#include "BlockTopology.hxx"
#include "MemArray.hxx"
-#include "MEDCouplingSMesh.hxx"
+#include "MEDCouplingRMesh.hxx"
#include "InterpKernelUtilities.hxx"
#include <iostream>
namespace ParaMEDMEM
{
- ParaGRID::ParaGRID(MEDCouplingSMesh* global_grid, Topology* topology) throw(INTERP_KERNEL::Exception)
+ ParaGRID::ParaGRID(MEDCouplingRMesh* global_grid, Topology* topology) throw(INTERP_KERNEL::Exception)
{
_block_topology = dynamic_cast<BlockTopology*>(topology);
coordinates_names.push_back(array->getName());
coordinates_units.push_back(array->getInfoOnComponentAt(0));
}
- _grid=MEDCouplingSMesh::New();
+ _grid=MEDCouplingRMesh::New();
_grid->set(xyz_array, coordinates_names,coordinates_units);
_grid->setName(global_grid->getName());
_grid->setDescription(global_grid->getDescription());*/
{
class Topology;
class BlockTopology;
- class MEDCouplingSMesh;
+ class MEDCouplingRMesh;
class ParaGRID
{
public:
- ParaGRID(MEDCouplingSMesh* global_grid, Topology* topology) throw(INTERP_KERNEL::Exception);
+ ParaGRID(MEDCouplingRMesh* global_grid, Topology* topology) throw(INTERP_KERNEL::Exception);
BlockTopology * getBlockTopology() const { return _block_topology; }
virtual ~ParaGRID();
- MEDCouplingSMesh* getGrid() const { return _grid; }
+ MEDCouplingRMesh* getGrid() const { return _grid; }
private:
- MEDCouplingSMesh* _grid;
+ MEDCouplingRMesh* _grid;
// structured grid topology
ParaMEDMEM::BlockTopology* _block_topology;
// stores the x,y,z axes on the global grid
namespace ParaMEDMEM
{
- ParaMESH::ParaMESH( MEDCouplingUMesh *subdomain_mesh, MEDCouplingUMesh *subdomain_face,
+ ParaMESH::ParaMESH( MEDCouplingPointSet *subdomain_mesh, MEDCouplingPointSet *subdomain_face,
DataArrayInt *CorrespElt_local2global, DataArrayInt *CorrespFace_local2global,
DataArrayInt *CorrespNod_local2global, const ProcessorGroup& proc_group ):
_cell_mesh(subdomain_mesh),
CorrespNod_local2global->incrRef();
}
- ParaMESH::ParaMESH( MEDCouplingUMesh *mesh, const ProcessorGroup& proc_group, const std::string& name):
+ ParaMESH::ParaMESH( MEDCouplingPointSet *mesh, const ProcessorGroup& proc_group, const std::string& name):
_cell_mesh(mesh),
_face_mesh(0),
_my_domain_id(proc_group.myRank()),
#ifndef __PARAMESH_HXX__
#define __PARAMESH_HXX__
-#include "MEDCouplingUMesh.hxx"
+#include "MEDCouplingPointSet.hxx"
#include "ProcessorGroup.hxx"
+#include "MemArray.hxx"
#include <string>
#include <vector>
class ParaMESH
{
public:
- ParaMESH( MEDCouplingUMesh *subdomain_mesh,
- MEDCouplingUMesh *subdomain_face,
+ ParaMESH( MEDCouplingPointSet *subdomain_mesh,
+ MEDCouplingPointSet *subdomain_face,
DataArrayInt *CorrespElt_local2global,
DataArrayInt *CorrespFace_local2global,
DataArrayInt *CorrespNod_local2global,
const ProcessorGroup& proc_group ) ;
- ParaMESH( MEDCouplingUMesh *mesh,
+ ParaMESH( MEDCouplingPointSet *mesh,
const ProcessorGroup& proc_group, const std::string& name);
virtual ~ParaMESH();
Topology* getTopology() const { return _explicit_topology; }
bool isStructured() const { return _cell_mesh->isStructured(); }
- MEDCouplingUMesh *getCellMesh() const { return _cell_mesh; }
- MEDCouplingUMesh *getFaceMesh() const { return _face_mesh; }
+ MEDCouplingPointSet *getCellMesh() const { return _cell_mesh; }
+ MEDCouplingPointSet *getFaceMesh() const { return _face_mesh; }
BlockTopology* getBlockTopology() const { return _block_topology; }
const int* getGlobalNumberingNode() const { return _node_global->getPointer(); }
private:
//mesh object underlying the ParaMESH object
- MEDCouplingUMesh *_cell_mesh ;
- MEDCouplingUMesh *_face_mesh ;
+ MEDCouplingPointSet *_cell_mesh ;
+ MEDCouplingPointSet *_face_mesh ;
//id of the local grid
int _my_domain_id;
ParaMEDMEMTest_IntersectionDEC.cxx \
ParaMEDMEMTest_StructuredCoincidentDEC.cxx \
ParaMEDMEMTest_MEDLoader.cxx \
+ ParaMEDMEMTest_ICocoTrio.cxx \
MPIAccessDECTest.cxx \
test_AllToAllDEC.cxx \
test_AllToAllvDEC.cxx \
#endif
CPPUNIT_TEST(testStructuredCoincidentDEC);
CPPUNIT_TEST(testStructuredCoincidentDEC);
+ CPPUNIT_TEST(testICocoTrio1);
CPPUNIT_TEST(testMEDLoaderRead1);
CPPUNIT_TEST(testMEDLoaderPolygonRead);
CPPUNIT_TEST(testMEDLoaderPolyhedronRead);
void testAsynchronousSlowSourceIntersectionDEC_2D();
void testAsynchronousFastSourceIntersectionDEC_2D();
//
+ void testICocoTrio1();
+ //
void testMEDLoaderRead1();
void testMEDLoaderPolygonRead();
void testMEDLoaderPolyhedronRead();
--- /dev/null
+#include "ParaMEDMEMTest.hxx"
+#include <string>
+#include "CommInterface.hxx"
+#include "ProcessorGroup.hxx"
+#include "MPIProcessorGroup.hxx"
+#include "DEC.hxx"
+#include "IntersectionDEC.hxx"
+#include <set>
+#include <time.h>
+#include "ICoCoTrioField.hxx"
+#include <iostream>
+#include <assert.h>
+
+using namespace std;
+using namespace ParaMEDMEM;
+using namespace ICoCo;
+
+typedef enum {sync_and,sync_or} synctype;
+void synchronize_bool(bool& stop, synctype s)
+{
+ int my_stop;
+ int my_stop_temp = stop?1:0;
+
+ if (s==sync_and)
+ MPI_Allreduce(&my_stop_temp,&my_stop,1,MPI_INTEGER,MPI_MIN,MPI_COMM_WORLD);
+ else if (s==sync_or)
+ MPI_Allreduce(&my_stop_temp,&my_stop,1,MPI_INTEGER,MPI_MAX,MPI_COMM_WORLD);
+ stop =(my_stop==1);
+}
+
+void synchronize_dt(double& dt)
+{
+ double dttemp=dt;
+ MPI_Allreduce(&dttemp,&dt,1,MPI_DOUBLE,MPI_MIN,MPI_COMM_WORLD);
+}
+
+
+void affiche( const TrioField& field)
+{
+ cout <<field.getName()<<endl;
+ for (int ele=0;ele<field._nb_elems;ele++)
+ cout <<ele <<": "<<field._field[ele]<<endl;;
+
+}
+
+void remplit_coord(double* coords)
+{
+ coords[0*3+0]=0.;
+ coords[0*3+1]=0.;
+ coords[0*3+2]=0.;
+
+ coords[1*3+0]=1.;
+ coords[1*3+1]=0.;
+ coords[1*3+2]=0.;
+
+
+ coords[2*3+0]=0.;
+ coords[2*3+1]=0.;
+ coords[2*3+2]=1.;
+
+ coords[3*3+0]=1.;
+ coords[3*3+1]=0.;
+ coords[3*3+2]=1.;
+
+ for (int i=4;i<8;i++)
+ {
+ for (int d=0;d<3;d++)
+ coords[i*3+d]=coords[(i-4)*3+d];
+ coords[i*3+1]+=1e-5;
+ }
+
+}
+
+void init_quad(TrioField& champ_quad)
+{
+
+ champ_quad.setName("champ_quad");
+ champ_quad._space_dim=3;
+ champ_quad._mesh_dim=2;
+ champ_quad._nbnodes=8;
+ champ_quad._nodes_per_elem=4;
+ champ_quad._nb_elems=2;
+ champ_quad._itnumber=0;
+ champ_quad._time1=0;
+ champ_quad._time2=1;
+ champ_quad._nb_field_components=1;
+
+ champ_quad._coords=new double[champ_quad._nbnodes*champ_quad._space_dim];
+ //memcpy(afield._coords,sommets.addr(),champ_quad._nbnodes*champ_quad._space_dim*sizeof(double));
+
+ remplit_coord(champ_quad._coords);
+
+
+ champ_quad._connectivity=new int[champ_quad._nb_elems*champ_quad._nodes_per_elem];
+ champ_quad._connectivity[0*champ_quad._nodes_per_elem+0]=0;
+ champ_quad._connectivity[0*champ_quad._nodes_per_elem+1]=1;
+ champ_quad._connectivity[0*champ_quad._nodes_per_elem+2]=3;
+ champ_quad._connectivity[0*champ_quad._nodes_per_elem+3]=2;
+ champ_quad._connectivity[1*champ_quad._nodes_per_elem+0]=4;
+ champ_quad._connectivity[1*champ_quad._nodes_per_elem+1]=5;
+ champ_quad._connectivity[1*champ_quad._nodes_per_elem+2]=7;
+ champ_quad._connectivity[1*champ_quad._nodes_per_elem+3]=6;
+
+
+ champ_quad._has_field_ownership=false;
+ champ_quad._field=0;
+ //champ_quad._field=new double[champ_quad._nb_elems];
+ // assert(champ_quad._nb_field_components==1);
+}
+void init_triangle(TrioField& champ_triangle)
+{
+
+ champ_triangle.setName("champ_triangle");
+ champ_triangle._space_dim=3;
+ champ_triangle._mesh_dim=2;
+ champ_triangle._nbnodes=8;
+ champ_triangle._nodes_per_elem=3;
+ champ_triangle._nb_elems=4;
+ champ_triangle._itnumber=0;
+ champ_triangle._time1=0;
+ champ_triangle._time2=1;
+ champ_triangle._nb_field_components=1;
+
+ champ_triangle._coords=new double[champ_triangle._nbnodes*champ_triangle._space_dim];
+ //memcpy(afield._coords,sommets.addr(),champ_triangle._nbnodes*champ_triangle._space_dim*sizeof(double));
+ remplit_coord(champ_triangle._coords);
+
+ champ_triangle._connectivity=new int[champ_triangle._nb_elems*champ_triangle._nodes_per_elem];
+ champ_triangle._connectivity[0*champ_triangle._nodes_per_elem+0]=0;
+ champ_triangle._connectivity[0*champ_triangle._nodes_per_elem+1]=1;
+ champ_triangle._connectivity[0*champ_triangle._nodes_per_elem+2]=2;
+ champ_triangle._connectivity[1*champ_triangle._nodes_per_elem+0]=1;
+ champ_triangle._connectivity[1*champ_triangle._nodes_per_elem+1]=3;
+ champ_triangle._connectivity[1*champ_triangle._nodes_per_elem+2]=2;
+
+ champ_triangle._connectivity[2*champ_triangle._nodes_per_elem+0]=4;
+ champ_triangle._connectivity[2*champ_triangle._nodes_per_elem+1]=5;
+ champ_triangle._connectivity[2*champ_triangle._nodes_per_elem+2]=7;
+ champ_triangle._connectivity[3*champ_triangle._nodes_per_elem+0]=4;
+ champ_triangle._connectivity[3*champ_triangle._nodes_per_elem+1]=7;
+ champ_triangle._connectivity[3*champ_triangle._nodes_per_elem+2]=6;
+
+ champ_triangle._has_field_ownership=false;
+ // champ_triangle._field=new double[champ_triangle._nb_elems];
+ champ_triangle._field=0;
+}
+
+void ParaMEDMEMTest::testICocoTrio1()
+{
+ int size;
+ int rank;
+ MPI_Comm_size(MPI_COMM_WORLD,&size);
+ MPI_Comm_rank(MPI_COMM_WORLD,&rank);
+
+ //the test is meant to run on five processors
+ if (size !=2) return ;
+
+ CommInterface comm;
+ set<int> emetteur_ids;
+ set<int> recepteur_ids;
+ emetteur_ids.insert(0);
+ recepteur_ids.insert(1);
+
+ MPIProcessorGroup recepteur_group(comm,recepteur_ids);
+ MPIProcessorGroup emetteur_group(comm,emetteur_ids);
+
+
+ string cas;
+ if (recepteur_group.containsMyRank())
+ {
+ cas="recepteur";
+
+ }
+ else
+ cas="emetteur";
+
+ IntersectionDEC dec_emetteur(emetteur_group, recepteur_group);
+
+ TrioField champ_emetteur, champ_recepteur;
+
+ init_triangle(champ_emetteur);
+ //init_triangle(champ_emetteur);
+ init_quad(champ_recepteur);
+ //init_emetteur(champ_recepteur);
+
+ if (cas=="emetteur")
+ {
+ champ_emetteur._field=new double[champ_emetteur._nb_elems];
+ for (int ele=0;ele<champ_emetteur._nb_elems;ele++)
+ champ_emetteur._field[ele]=1;
+
+ champ_emetteur._has_field_ownership=true;
+ }
+
+
+ MPI_Barrier(MPI_COMM_WORLD);
+
+ clock_t clock0= clock ();
+ int compti=0;
+
+ bool init=true; // first time step ??
+ bool stop=false;
+ //boucle sur les pas de quads
+ while (!stop) {
+
+ compti++;
+ clock_t clocki= clock ();
+ cout << compti << " CLOCK " << (clocki-clock0)*1.e-6 << endl;
+ for (int non_unif=0;non_unif<2;non_unif++)
+ {
+ // if (champ_recepteur._field)
+ // delete [] champ_recepteur._field;
+ champ_recepteur._field=0;
+ // champ_recepteur._has_field_ownership=false;
+
+
+
+ if (cas=="emetteur")
+ if (non_unif)
+ champ_emetteur._field[0]=40;
+ bool ok=false; // Is the time interval successfully solved ?
+
+ // Loop on the time interval tries
+ if(1) {
+
+
+ if (cas=="emetteur")
+ dec_emetteur.attachLocalField((ICoCo::Field*) &champ_emetteur);
+ else
+ dec_emetteur.attachLocalField((ICoCo::Field*) &champ_recepteur);
+
+
+ if(init)
+ dec_emetteur.synchronize();
+ init=false;
+
+ if (cas=="emetteur") {
+ dec_emetteur.sendData();
+ affiche(champ_emetteur);
+ }
+ else if (cas=="recepteur") {
+ dec_emetteur.recvData();
+ affiche(champ_recepteur);
+ }
+ else
+ throw 0;
+ }
+ stop=true;
+ }
+}
+}
// ParaMEDMEM::ParaSUPPORT* parasupport=new UnstructuredParaSUPPORT( support,*source_group);
ParaMEDMEM::ComponentTopology comptopo;
if(srcM=="P0")
- parafield = new ParaFIELD(ON_CELLS,paramesh, comptopo);
+ parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo);
else
- parafield = new ParaFIELD(ON_NODES,paramesh, comptopo);
+ parafield = new ParaFIELD(ON_NODES,NO_TIME,paramesh, comptopo);
int nb_local;
if(srcM=="P0")
nb_local=mesh->getNumberOfCells();
// ParaMEDMEM::ParaSUPPORT* parasupport=new UnstructuredParaSUPPORT(support,*target_group);
ParaMEDMEM::ComponentTopology comptopo;
if(targetM=="P0")
- parafield = new ParaFIELD(ON_CELLS,paramesh, comptopo);
+ parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo);
else
- parafield = new ParaFIELD(ON_NODES,paramesh, comptopo);
+ parafield = new ParaFIELD(ON_NODES,NO_TIME,paramesh, comptopo);
int nb_local;
if(targetM=="P0")
nb_local=mesh->getNumberOfCells();
// ParaMEDMEM::ParaSUPPORT* parasupport=new UnstructuredParaSUPPORT( support,*source_group);
ParaMEDMEM::ComponentTopology comptopo;
if(srcM=="P0")
- parafield = new ParaFIELD(ON_CELLS,paramesh, comptopo);
+ parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo);
else
- parafield = new ParaFIELD(ON_NODES,paramesh, comptopo);
+ parafield = new ParaFIELD(ON_NODES,NO_TIME,paramesh, comptopo);
int nb_local;
if(srcM=="P0")
nb_local=mesh->getNumberOfCells();
// ParaMEDMEM::ParaSUPPORT* parasupport=new UnstructuredParaSUPPORT(support,*target_group);
ParaMEDMEM::ComponentTopology comptopo;
if(targetM=="P0")
- parafield = new ParaFIELD(ON_CELLS,paramesh, comptopo);
+ parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo);
else
- parafield = new ParaFIELD(ON_NODES,paramesh, comptopo);
+ parafield = new ParaFIELD(ON_NODES,NO_TIME,paramesh, comptopo);
int nb_local;
if(targetM=="P0")
nb_local=mesh->getNumberOfCells();
// ParaMEDMEM::ParaSUPPORT* parasupport=new UnstructuredParaSUPPORT( support,*source_group);
ParaMEDMEM::ComponentTopology comptopo;
if(srcM=="P0")
- parafield = new ParaFIELD(ON_CELLS,paramesh, comptopo);
+ parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo);
else
- parafield = new ParaFIELD(ON_NODES,paramesh, comptopo);
+ parafield = new ParaFIELD(ON_NODES,NO_TIME,paramesh, comptopo);
int nb_local;
if(srcM=="P0")
// ParaMEDMEM::ParaSUPPORT* parasupport=new UnstructuredParaSUPPORT(support,*target_group);
ParaMEDMEM::ComponentTopology comptopo;
if(targetM=="P0")
- parafield = new ParaFIELD(ON_CELLS,paramesh, comptopo);
+ parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo);
else
- parafield = new ParaFIELD(ON_NODES,paramesh, comptopo);
+ parafield = new ParaFIELD(ON_NODES,NO_TIME,paramesh, comptopo);
int nb_local;
if(targetM=="P0")
paramesh=new ParaMESH (mesh,source_group,"source mesh");
ParaMEDMEM::ComponentTopology comptopo(6);
- parafield = new ParaFIELD(ON_CELLS,paramesh, comptopo);
+ parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo);
int nb_local=mesh->getNumberOfCells();
const int* global_numbering = paramesh->getGlobalNumberingCell();
paramesh=new ParaMESH (mesh,self_group,"target mesh");
ParaMEDMEM::ComponentTopology comptopo(6, &target_group);
- parafield = new ParaFIELD(ON_CELLS,paramesh, comptopo);
+ parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo);
int nb_local=mesh->getNumberOfCells();
double *value=parafield->getField()->getArray()->getPointer();
paramesh=new ParaMESH (mesh,*source_group,"source mesh");
ParaMEDMEM::ComponentTopology comptopo;
- parafield = new ParaFIELD(ON_CELLS, paramesh, comptopo);
+ parafield = new ParaFIELD(ON_CELLS, NO_TIME, paramesh, comptopo);
int nb_local=mesh->getNumberOfCells();
double *value=parafield->getField()->getArray()->getPointer();
paramesh=new ParaMESH (mesh,*target_group,"target mesh");
ParaMEDMEM::ComponentTopology comptopo;
- parafield = new ParaFIELD(ON_CELLS,paramesh, comptopo);
+ parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo);
int nb_local=mesh->getNumberOfCells();
double *value=parafield->getField()->getArray()->getPointer();
