static const NumberingPolicy numPol=MyMeshType::My_numPol;
public:
ConvexIntersector(const MyMeshType& meshT, const MyMeshType& meshS,
- double dimCaracteristic, double precision, double medianPlane,
+ double dimCaracteristic, double precision, double md3DSurf, double medianPlane,
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);
{
template<class MyMeshType, class MyMatrix, template <class MeshType, class TheMatrix, class ThisIntersector> class InterpType>
ConvexIntersector<MyMeshType,MyMatrix,InterpType>::ConvexIntersector(const MyMeshType& meshT, const MyMeshType& meshS,
- double dimCaracteristic, double precision,
- double medianPlane, bool doRotate , int oriantation, int printLevel)
- :InterpType<MyMeshType,MyMatrix,ConvexIntersector<MyMeshType,MyMatrix,InterpType> >(meshT,meshS,dimCaracteristic, precision, medianPlane, doRotate, oriantation, printLevel),
+ double dimCaracteristic, double precision, double md3DSurf,
+ double medianPlane, bool doRotate , int oriantation, int printLevel)
+ :InterpType<MyMeshType,MyMatrix,ConvexIntersector<MyMeshType,MyMatrix,InterpType> >(meshT,meshS,dimCaracteristic, precision, md3DSurf, medianPlane, doRotate, oriantation, printLevel),
_epsilon(precision*dimCaracteristic)
{
if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 1)
static const NumberingPolicy numPol=MyMeshType::My_numPol;
public:
Geometric2DIntersector(const MyMeshType& meshT, const MyMeshType& meshS,
- double dimCaracteristic, double medianPlane, double precision, int orientation);
+ double dimCaracteristic, double md3DSurf, 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);
{
template<class MyMeshType, class MyMatrix, template <class MeshType, class TheMatrix, class ThisIntersector> class InterpType>
Geometric2DIntersector<MyMeshType,MyMatrix,InterpType>::Geometric2DIntersector(const MyMeshType& meshT, const MyMeshType& meshS,
- double dimCaracteristic, double medianPlane, double precision, int orientation):
- InterpType<MyMeshType,MyMatrix,Geometric2DIntersector<MyMeshType,MyMatrix,InterpType> >(meshT,meshS,dimCaracteristic, precision, medianPlane, true, orientation, 0)
+ double dimCaracteristic, double md3DSurf, double medianPlane, double precision, int orientation):
+ InterpType<MyMeshType,MyMatrix,Geometric2DIntersector<MyMeshType,MyMatrix,InterpType> >(meshT,meshS,dimCaracteristic, precision, md3DSurf, medianPlane, true, orientation, 0)
{
QUADRATIC_PLANAR::_precision=dimCaracteristic*precision;
}
Interpolation3DSurf(const InterpolationOptions& io);
void setOptions(double precision, int printLevel, double medianPlane,
IntersectionType intersectionType, bool doRotate, int orientation=0);
-
public:
- bool doRotate() const { return _do_rotate; }
- double medianPlane() const { return _median_plane; }
- double surf3DAdjustmentEps() const { return _surf_3D_adjustment_eps; }
- void setSurf3DAdjustmentEps(double val) { _surf_3D_adjustment_eps=val; }
template<class MyMeshType, class MyMatrixRow>
- void performAdjustmentOfBB(PlanarIntersector<MyMeshType,MyMatrixRow>* intersector, std::vector<double>& bbox) const
- { intersector->adjustBoundingBoxes(bbox,_surf_3D_adjustment_eps); }
- protected:
- bool _do_rotate;
- double _median_plane;
- double _surf_3D_adjustment_eps;
- static const double DFT_MEDIAN_PLANE;
- static const double DFT_SURF3D_ADJ_EPS;
+ void performAdjustmentOfBB(PlanarIntersector<MyMeshType,MyMatrixRow>* intersector, std::vector<double>& bbox) const
+ { intersector->adjustBoundingBoxes(bbox,InterpolationPlanar<Interpolation3DSurf>::getBoundingBoxAdjustment(),InterpolationPlanar<Interpolation3DSurf>::getBoundingBoxAdjustmentAbs()); }
};
}
namespace INTERP_KERNEL
{
- const double Interpolation3DSurf::DFT_MEDIAN_PLANE=0.5;
- const double Interpolation3DSurf::DFT_SURF3D_ADJ_EPS=1e-4;
-
- Interpolation3DSurf::Interpolation3DSurf():_do_rotate(true)
- ,_median_plane(DFT_MEDIAN_PLANE)
- ,_surf_3D_adjustment_eps(DFT_SURF3D_ADJ_EPS)
+ Interpolation3DSurf::Interpolation3DSurf()
{
}
Interpolation3DSurf::Interpolation3DSurf(const InterpolationOptions& io):InterpolationPlanar<Interpolation3DSurf>(io)
- ,_median_plane(io.getMedianPlane())
- ,_surf_3D_adjustment_eps(io.getBoundingBoxAdjustment())
{
}
IntersectionType intersectionType, bool doRotate, int orientation)
{
InterpolationPlanar<Interpolation3DSurf>::setOptions(precision,printLevel,intersectionType, orientation);
- _do_rotate=doRotate;
- _median_plane=medianPlane;
+ InterpolationPlanar<Interpolation3DSurf>::setDoRotate(doRotate);
+ InterpolationPlanar<Interpolation3DSurf>::setMedianPlane(medianPlane);
}
}
--- /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 "InterpolationOptions.hxx"
+
+const double INTERP_KERNEL::InterpolationOptions::DFT_MEDIAN_PLANE=0.5;
+
+const double INTERP_KERNEL::InterpolationOptions::DFT_SURF3D_ADJ_EPS=1.e-4;
+
+const double INTERP_KERNEL::InterpolationOptions::DFT_MAX_DIST_3DSURF_INTERSECT=-1.;
#ifndef __INTERPOLATIONOPTIONS_HXX__
#define __INTERPOLATIONOPTIONS_HXX__
-
-namespace INTERP_KERNEL {
+namespace INTERP_KERNEL
+{
typedef enum { Triangulation, Convex, Geometric2D } IntersectionType;
/// Type describing the different ways in which the hexahedron can be split into tetrahedra.
/// The PLANAR_* policies persume that each face is to be considered planar, while the general
/// policies make no such hypothesis. The integer at the end gives the number of tetrahedra
/// that result from the split.
typedef enum { PLANAR_FACE_5 = 5, PLANAR_FACE_6 = 6, GENERAL_24 = 24, GENERAL_48 = 48 } SplittingPolicy;
-
- class InterpolationOptions{
+ class InterpolationOptions
+ {
private :
int _print_level ;
IntersectionType _intersection_type;
double _precision;
double _median_plane ;
bool _do_rotate ;
+ //! this measure is relative to the caracteristic dimension
double _bounding_box_adjustment ;
+ //! this measure is absolute \b not relative to the cell size
+ double _bounding_box_adjustment_abs ;
+ double _max_distance_for_3Dsurf_intersect;
int _orientation ;
SplittingPolicy _splitting_policy ;
int getPrintLevel() const { return _print_level; }
void setPrintLevel(int pl) { _print_level=pl; }
- IntersectionType getIntersectionType() const { return InterpolationOptions::_intersection_type; }
- void setIntersectionType(IntersectionType it) { InterpolationOptions::_intersection_type=it; }
+ IntersectionType getIntersectionType() const { return _intersection_type; }
+ void setIntersectionType(IntersectionType it) { _intersection_type=it; }
- double getPrecision() const { return InterpolationOptions::_precision; }
- void setPrecision(double p) { InterpolationOptions::_precision=p; }
+ double getPrecision() const { return _precision; }
+ void setPrecision(double p) { _precision=p; }
- double getMedianPlane() const { return InterpolationOptions::_median_plane; }
- void setMedianPlane(double mp) { InterpolationOptions::_median_plane=mp; }
+ double getMedianPlane() const { return _median_plane; }
+ void setMedianPlane(double mp) { _median_plane=mp; }
- bool getDoRotate() const { return InterpolationOptions::_do_rotate; }
- void setDoRotate( bool dr) { InterpolationOptions::_do_rotate = dr; }
+ bool getDoRotate() const { return _do_rotate; }
+ void setDoRotate( bool dr) { _do_rotate = dr; }
- double getBoundingBoxAdjustment() const { return InterpolationOptions::_bounding_box_adjustment; }
- void setBoundingBoxAdjustment(double bba) { InterpolationOptions::_bounding_box_adjustment=bba; }
+ double getBoundingBoxAdjustment() const { return _bounding_box_adjustment; }
+ void setBoundingBoxAdjustment(double bba) { _bounding_box_adjustment=bba; }
+
+ double getBoundingBoxAdjustmentAbs() const { return _bounding_box_adjustment_abs; }
+ void setBoundingBoxAdjustmentAbs(double bba) { _bounding_box_adjustment_abs=bba; }
- int getOrientation() const { return InterpolationOptions::_orientation; }
- void setOrientation(int o) { InterpolationOptions::_orientation=o; }
+ double getMaxDistance3DSurfIntersect() const { return _max_distance_for_3Dsurf_intersect; }
+ void setMaxDistance3DSurfIntersect(double bba) { _max_distance_for_3Dsurf_intersect=bba; }
+
+ int getOrientation() const { return _orientation; }
+ void setOrientation(int o) { _orientation=o; }
SplittingPolicy getSplittingPolicy() const { return _splitting_policy; }
void setSplittingPolicy(SplittingPolicy sp) { _splitting_policy=sp; }
_print_level=0;
_intersection_type=Triangulation;
_precision=1e-12;
- _median_plane=0.5;
+ _median_plane=DFT_MEDIAN_PLANE;
_do_rotate=true;
- _bounding_box_adjustment=1e-4;
+ _bounding_box_adjustment=DFT_SURF3D_ADJ_EPS;
+ _bounding_box_adjustment_abs=0.;
+ _max_distance_for_3Dsurf_intersect=DFT_MAX_DIST_3DSURF_INTERSECT;
_orientation=0;
_splitting_policy=GENERAL_48;
}
+ private:
+ static const double DFT_MEDIAN_PLANE;
+ static const double DFT_SURF3D_ADJ_EPS;
+ static const double DFT_MAX_DIST_3DSURF_INTERSECT;
};
}
case Triangulation:
intersector=new TriangulationIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P0>(myMeshT,myMeshS,_dim_caracteristic,
InterpolationOptions::getPrecision(),
+ InterpolationOptions::getMaxDistance3DSurfIntersect(),
InterpolationOptions::getMedianPlane(),
InterpolationOptions::getOrientation(),
InterpolationOptions::getPrintLevel());
case Convex:
intersector=new ConvexIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P0>(myMeshT,myMeshS,_dim_caracteristic,
InterpolationOptions::getPrecision(),
+ InterpolationOptions::getMaxDistance3DSurfIntersect(),
InterpolationOptions::getDoRotate(),
InterpolationOptions::getMedianPlane(),
InterpolationOptions::getOrientation(),
break;
case Geometric2D:
intersector=new Geometric2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P0>(myMeshT, myMeshS, _dim_caracteristic,
+ InterpolationOptions::getMaxDistance3DSurfIntersect(),
InterpolationOptions::getMedianPlane(),
InterpolationOptions::getPrecision(),
InterpolationOptions::getOrientation());
case Triangulation:
intersector=new TriangulationIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P1>(myMeshT,myMeshS,_dim_caracteristic,
InterpolationOptions::getPrecision(),
+ InterpolationOptions::getMaxDistance3DSurfIntersect(),
InterpolationOptions::getMedianPlane(),
InterpolationOptions::getOrientation(),
InterpolationOptions::getPrintLevel());
case Convex:
intersector=new ConvexIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P1>(myMeshT,myMeshS,_dim_caracteristic,
InterpolationOptions::getPrecision(),
+ InterpolationOptions::getMaxDistance3DSurfIntersect(),
InterpolationOptions::getDoRotate(),
InterpolationOptions::getMedianPlane(),
InterpolationOptions::getOrientation(),
break;
case Geometric2D:
intersector=new Geometric2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP0P1>(myMeshT, myMeshS, _dim_caracteristic,
+ InterpolationOptions::getMaxDistance3DSurfIntersect(),
InterpolationOptions::getMedianPlane(),
InterpolationOptions::getPrecision(),
InterpolationOptions::getOrientation());
case Triangulation:
intersector=new TriangulationIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P0>(myMeshT,myMeshS,_dim_caracteristic,
InterpolationOptions::getPrecision(),
+ InterpolationOptions::getMaxDistance3DSurfIntersect(),
InterpolationOptions::getMedianPlane(),
InterpolationOptions::getOrientation(),
InterpolationOptions::getPrintLevel());
case Convex:
intersector=new ConvexIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P0>(myMeshT,myMeshS,_dim_caracteristic,
InterpolationOptions::getPrecision(),
+ InterpolationOptions::getMaxDistance3DSurfIntersect(),
InterpolationOptions::getDoRotate(),
InterpolationOptions::getMedianPlane(),
InterpolationOptions::getOrientation(),
break;
case Geometric2D:
intersector=new Geometric2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P0>(myMeshT, myMeshS, _dim_caracteristic,
+ InterpolationOptions::getMaxDistance3DSurfIntersect(),
InterpolationOptions::getMedianPlane(),
InterpolationOptions::getPrecision(),
InterpolationOptions::getOrientation());
case Triangulation:
intersector=new TriangulationIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P1>(myMeshT,myMeshS,_dim_caracteristic,
InterpolationOptions::getPrecision(),
+ InterpolationOptions::getMaxDistance3DSurfIntersect(),
InterpolationOptions::getMedianPlane(),
InterpolationOptions::getOrientation(),
InterpolationOptions::getPrintLevel());
case Convex:
intersector=new ConvexIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P1>(myMeshT,myMeshS,_dim_caracteristic,
InterpolationOptions::getPrecision(),
+ InterpolationOptions::getMaxDistance3DSurfIntersect(),
InterpolationOptions::getDoRotate(),
InterpolationOptions::getMedianPlane(),
InterpolationOptions::getOrientation(),
break;
case Geometric2D:
intersector=new Geometric2DIntersector<MyMeshType,MatrixType,PlanarIntersectorP1P1>(myMeshT, myMeshS, _dim_caracteristic,
+ InterpolationOptions::getMaxDistance3DSurfIntersect(),
InterpolationOptions::getMedianPlane(),
InterpolationOptions::getPrecision(),
InterpolationOptions::getOrientation());
BoundingBox.cxx \
TetraAffineTransform.cxx\
CellModel.cxx\
- UnitTetraIntersectionBary.cxx
+ UnitTetraIntersectionBary.cxx \
+ InterpolationOptions.cxx
libinterpkernel_la_CPPFLAGS=-I$(srcdir)/Geometric2D -I$(srcdir)/Bases
static const NumberingPolicy numPol=MyMeshType::My_numPol;
public:
//! \addtogroup InterpKerGrpIntPlan @{
- PlanarIntersector(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double medianPlane, bool doRotate, int orientation, int printLevel);
+ PlanarIntersector(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double md3DSurf, double medianPlane, bool doRotate, int orientation, int printLevel);
//! @}
virtual ~PlanarIntersector();
void createBoundingBoxes(const MyMeshType& mesh, std::vector<double>& bbox);
- void adjustBoundingBoxes(std::vector<double>& bbox, double Surf3DAdjustmentEps);
+ void adjustBoundingBoxes(std::vector<double>& bbox, double surf3DAdjustmentEps, double surf3DAdjustmentEpsAbs);
inline void getElemBB(double* bb, const MyMeshType& mesh, ConnType iP, ConnType nb_nodes);
protected :
int projectionThis(double *Coords_A, double *Coords_B, int nb_NodesA, int nb_NodesB);
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);
+ double getValueRegardingOption(double val) const;
+ static int projection(double *Coords_A, double *Coords_B,
+ int nb_NodesA, int nb_NodesB, double epsilon, double md3DSurf, double median_plane, bool do_rotate);
static void rotate3DTriangle( double* PP1, double*PP2, double*PP3,
TranslationRotationMatrix& rotation_matrix);
protected:
const MyMeshType& _meshT;
const MyMeshType& _meshS;
double _dim_caracteristic;
+ double _max_distance_3Dsurf_intersect;
double _precision;
double _median_plane;
bool _do_rotate;
namespace INTERP_KERNEL
{
template<class MyMeshType, class MyMatrix>
- PlanarIntersector<MyMeshType,MyMatrix>::PlanarIntersector(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double medianPlane, bool doRotate, int orientation, int printLevel):
+ PlanarIntersector<MyMeshType,MyMatrix>::PlanarIntersector(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double md3DSurf, double medianPlane, bool doRotate, int orientation, int printLevel):
_meshT(meshT),_meshS(meshS),
- _dim_caracteristic(dimCaracteristic),_precision(precision),_median_plane(medianPlane),
+ _dim_caracteristic(dimCaracteristic),_max_distance_3Dsurf_intersect(md3DSurf),_precision(precision),_median_plane(medianPlane),
_do_rotate(doRotate),_orientation(orientation),_print_level(printLevel)
{
_connectT=meshT.getConnectivityPtr();
\param bbox vector containing the bounding boxes
*/
template<class MyMeshType, class MyMatrix>
- void PlanarIntersector<MyMeshType,MyMatrix>::adjustBoundingBoxes(std::vector<double>& bbox, double Surf3DAdjustmentEps)
+ void PlanarIntersector<MyMeshType,MyMatrix>::adjustBoundingBoxes(std::vector<double>& bbox, double surf3DAdjustmentEps, double surf3DAdjustmentEpsAbs)
{
/* We build the segment tree for locating possible matching intersections*/
}
for(int idim=0; idim<SPACEDIM; idim++)
{
- bbox[i*2*SPACEDIM+2*idim ] -= Surf3DAdjustmentEps*max;
- bbox[i*2*SPACEDIM+2*idim+1] += Surf3DAdjustmentEps*max;
+ bbox[i*2*SPACEDIM+2*idim ] -= surf3DAdjustmentEps*max+surf3DAdjustmentEpsAbs;
+ bbox[i*2*SPACEDIM+2*idim+1] += surf3DAdjustmentEps*max+surf3DAdjustmentEpsAbs;
}
}
}
template<class MyMeshType, class MyMatrix>
void PlanarIntersector<MyMeshType,MyMatrix>::getRealCoordinates(ConnType icellT, ConnType icellS, ConnType nbNodesT, ConnType nbNodesS, std::vector<double>& coordsT, std::vector<double>& coordsS, int& orientation)
{
- /*double epsilon=_precision*_dim_caracteristic;
- coordsT.resize(SPACEDIM*nbNodesT);
- coordsS.resize(SPACEDIM*nbNodesS);
- int nb_dist_NodesT=nbNodesT;
- int nb_dist_NodesS=nbNodesS;
- int i_last = nbNodesT - 1;
- const double * Pi_last=_coordsT +_connectT[OTT<ConnType,numPol>::conn2C(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]+i_last)];
-
- for (int iT=0; iT<nbNodesT; iT++)
- {
- const double * PiT = _coordsT + SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectT[OTT<ConnType,numPol>::conn2C(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]+iT)]);
- if(distance2<SPACEDIM>(Pi_last, PiT)>epsilon)
- {
- for (int idim=0; idim<SPACEDIM; idim++)
- coordsT[SPACEDIM*iT+idim]=PiT[idim];
- i_last=iT; Pi_last = PiT;
- }
- else
- nb_dist_NodesT--;
- }
- coordsT.resize(nb_dist_NodesT*SPACEDIM);
- i_last = nbNodesS - 1;
- Pi_last=_coordsS + SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+i_last)]);
- for (int iS=0; iS<nbNodesS; iS++)
- {
- const double * PiS=_coordsS+SPACEDIM*OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+iS)]);
- if(distance2<SPACEDIM>(Pi_last, PiS)>epsilon)
- {
- for (int idim=0; idim<SPACEDIM; idim++)
- coordsS[SPACEDIM*iS+idim]=PiS[idim];
- i_last=iS; Pi_last = PiS;
- }
- else
- nb_dist_NodesS--;
- }
- coordsS.resize(nb_dist_NodesS*SPACEDIM);
- //project cells T and S on the median plane
- // and rotate the median plane
- if(SPACEDIM==3)
- orientation = projectionThis(&coordsT[0], &coordsS[0], nb_dist_NodesT, nb_dist_NodesS);
-
- //DEBUG PRINTS
- if(_print_level >= 3)
- {
- std::cout << std::endl << "Cell coordinates (possibly after projection)" << std::endl;
- std::cout << std::endl << "icellT= " << icellT << ", nb nodes T= " << nbNodesT << std::endl;
- for(int iT =0; iT< nbNodesT; iT++)
- {for (int idim=0; idim<SPACEDIM; idim++) std::cout << coordsT[SPACEDIM*iT+idim] << " "; std::cout << std::endl;}
- std::cout << std::endl << "icellS= " << icellS << ", nb nodes S= " << nbNodesS << std::endl;
- for(int iS =0; iS< nbNodesS; iS++)
- {for (int idim=0; idim<SPACEDIM; idim++) std::cout << coordsS[SPACEDIM*iS+idim]<< " "; std::cout << std::endl; }
- }*/
coordsT.resize(SPACEDIM*nbNodesT);
coordsS.resize(SPACEDIM*nbNodesS);
for (int idim=0; idim<SPACEDIM; idim++)
{for (int idim=0; idim<SPACEDIM; idim++) std::cout << coordsS[SPACEDIM*iS+idim]<< " "; std::cout << std::endl; }
}
}
+
+ /*!
+ * Filtering out zero surfaces and badly oriented surfaces
+ * _orientation = -1,0,1,2
+ * -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
+ * 2 : the absolute value of intersection is always taken into account
+ */
+ template<class MyMeshType, class MyMatrix>
+ double PlanarIntersector<MyMeshType,MyMatrix>::getValueRegardingOption(double val) const
+ {
+ if(_orientation==0)
+ return val;
+ if(_orientation==2)
+ return fabs(val);
+ if (( val > 0.0 && _orientation==1) || ( val < 0.0 && _orientation==-1 ))
+ return val;
+ return 0.;
+ }
template<class MyMeshType, class MyMatrix>
int PlanarIntersector<MyMeshType,MyMatrix>::projectionThis(double *Coords_A, double *Coords_B, int nb_NodesA, int nb_NodesB)
{
- return projection(Coords_A,Coords_B,nb_NodesA,nb_NodesB,_dim_caracteristic*_precision,_median_plane,_do_rotate);
+ return projection(Coords_A,Coords_B,nb_NodesA,nb_NodesB,_dim_caracteristic*_precision,_max_distance_3Dsurf_intersect,_median_plane,_do_rotate);
}
template<class MyMeshType, class MyMatrix>
int PlanarIntersector<MyMeshType,MyMatrix>::projection(double *Coords_A, double *Coords_B,
- int nb_NodesA, int nb_NodesB, double epsilon, double median_plane, bool do_rotate)
+ int nb_NodesA, int nb_NodesB, double epsilon, double md3DSurf, double median_plane, bool do_rotate)
{
double normal_A[3]={0,0,0};
double normal_B[3]={0,0,0};
normB = sqrt(dotprod<SPACEDIM>(normal_B,normal_B));
}
+ //fabien option
+ if(md3DSurf>0.)
+ {
+ double coords_GA[3];
+ for (int i=0;i<3;i++)
+ {
+ coords_GA[i]=0.;
+ for (int j=0;j<nb_NodesA;j++)
+ coords_GA[i]+=Coords_A[3*j+i];
+ coords_GA[i]/=nb_NodesA;
+ }
+ double G1[3],G2[3],G3[3];
+ for (int i=0;i<3;i++)
+ {
+ G1[i]=Coords_B[i]-coords_GA[i];
+ G2[i]=Coords_B[i+3]-coords_GA[i];
+ G3[i]=Coords_B[i+6]-coords_GA[i];
+ }
+ double prodvect[3];
+ prodvect[0]=G1[1]*G2[2]-G1[2]*G2[1];
+ prodvect[1]=G1[2]*G2[0]-G1[0]*G2[2];
+ prodvect[2]=G1[0]*G2[1]-G1[1]*G2[0];
+ double prodscal=prodvect[0]*G3[0]+prodvect[1]*G3[1]+prodvect[2]*G3[2];
+ if(fabs(prodscal)>md3DSurf)
+ return 0;
+ }
if(i_A2<nb_NodesA && i_B2<nb_NodesB)
{
//Build the normal of the median plane
typedef typename MyMeshType::MyConnType ConnType;
static const NumberingPolicy numPol=MyMeshType::My_numPol;
protected:
- PlanarIntersectorP0P0(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double medianPlane, bool doRotate, int orientation, int printLevel);
+ PlanarIntersectorP0P0(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double md3DSurf, double medianPlane, bool doRotate, int orientation, int printLevel);
public:
int getNumberOfRowsOfResMatrix() const;
int getNumberOfColsOfResMatrix() const;
{
template<class MyMeshType, class MyMatrix, class ConcreteP0P0Intersector>
PlanarIntersectorP0P0<MyMeshType,MyMatrix,ConcreteP0P0Intersector>::PlanarIntersectorP0P0(const MyMeshType& meshT, const MyMeshType& meshS,
- double dimCaracteristic, double precision, double medianPlane,
+ double dimCaracteristic, double precision, double md3DSurf, double medianPlane,
bool doRotate, int orientation, int printLevel):
- PlanarIntersector<MyMeshType,MyMatrix>(meshT,meshS,dimCaracteristic,precision,medianPlane,doRotate,orientation,printLevel)
+ PlanarIntersector<MyMeshType,MyMatrix>(meshT,meshS,dimCaracteristic,precision,md3DSurf,medianPlane,doRotate,orientation,printLevel)
{
}
int iS=*iter;
int nbNodesS=PlanarIntersector<MyMeshType,MyMatrix>::_connIndexS[iS+1]-PlanarIntersector<MyMeshType,MyMatrix>::_connIndexS[iS];
double surf=intersectGeometry(OTT<ConnType,numPol>::indFC(icellT),OTT<ConnType,numPol>::indFC(iS),nbNodesT,nbNodesS);
- //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 ))
+ surf=PlanarIntersector<MyMeshType,MyMatrix>::getValueRegardingOption(surf);
+ if(surf!=0.)
resRow.insert(std::make_pair(OTT<ConnType,numPol>::indFC(iS),surf));
}
}
typedef typename MyMeshType::MyConnType ConnType;
static const NumberingPolicy numPol=MyMeshType::My_numPol;
protected:
- PlanarIntersectorP0P1(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double medianPlane, bool doRotate, int orientation, int printLevel);
+ PlanarIntersectorP0P1(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double md3DSurf, double medianPlane, bool doRotate, int orientation, int printLevel);
public:
void intersectCells(ConnType icellT, const std::vector<ConnType>& icellsS, MyMatrix& res);
int getNumberOfRowsOfResMatrix() const;
{
template<class MyMeshType, class MyMatrix, class ConcreteP0P1Intersector>
PlanarIntersectorP0P1<MyMeshType,MyMatrix,ConcreteP0P1Intersector>::PlanarIntersectorP0P1(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)
+ double dimCaracteristic, double precision, double md3DSurf, double medianPlane,
+ bool doRotate, int orientation, int printLevel):
+ PlanarIntersector<MyMeshType,MyMatrix>(meshT,meshS,dimCaracteristic,precision,md3DSurf,medianPlane,doRotate,orientation,printLevel)
{
}
orientation=PlanarIntersector<MyMeshType,MyMatrix>::projectionThis(&sourceCellCoordsTmp[0],quadrangle,sourceCellCoords.size()/SPACEDIM,4);
NormalizedCellType tS=PlanarIntersector<MyMeshType,MyMatrix>::_meshS.getTypeOfElement(iS);
double surf=orientation*intersectGeometryWithQuadrangle(quadrangle,sourceCellCoordsTmp,CellModel::getCellModel(tS).isQuadratic());
- //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 ))
+ surf=PlanarIntersector<MyMeshType,MyMatrix>::getValueRegardingOption(surf);
+ if(surf!=0.)
{
typename MyMatrix::value_type::const_iterator iterRes=resRow.find(OTT<ConnType,numPol>::indFC(iS));
if(iterRes==resRow.end())
typedef typename MyMeshType::MyConnType ConnType;
static const NumberingPolicy numPol=MyMeshType::My_numPol;
protected:
- PlanarIntersectorP1P0(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double medianPlane, bool doRotate, int orientation, int printLevel);
+ PlanarIntersectorP1P0(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double md3DSurf, double medianPlane, bool doRotate, int orientation, int printLevel);
public:
void intersectCells(ConnType icellT, const std::vector<ConnType>& icellsS, MyMatrix& res);
int getNumberOfRowsOfResMatrix() const;
{
template<class MyMeshType, class MyMatrix, class ConcreteP1P0Intersector>
PlanarIntersectorP1P0<MyMeshType,MyMatrix,ConcreteP1P0Intersector>::PlanarIntersectorP1P0(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)
+ double dimCaracteristic, double precision, double md3DSurf, double medianPlane,
+ bool doRotate, int orientation, int printLevel):
+ PlanarIntersector<MyMeshType,MyMatrix>(meshT,meshS,dimCaracteristic,precision,md3DSurf,medianPlane,doRotate,orientation,printLevel)
{
}
if(SPACEDIM==3)
orientation=PlanarIntersector<MyMeshType,MyMatrix>::projectionThis(&targetCellCoordsTmp[0],quadrangle,targetCellCoords.size()/SPACEDIM,4);
double surf=orientation*intersectGeometryWithQuadrangle(quadrangle,targetCellCoordsTmp,isTargetQuad);
- if (( surf > 0.0 && PlanarIntersector<MyMeshType,MyMatrix>::_orientation >=0 ) || ( surf < 0.0 && PlanarIntersector<MyMeshType,MyMatrix>::_orientation <=0 ))
+ surf=PlanarIntersector<MyMeshType,MyMatrix>::getValueRegardingOption(surf);
+ if(surf!=0.)
{
typename MyMatrix::value_type::const_iterator iterRes=resRow.find(OTT<ConnType,numPol>::indFC(curNodeSInCmode));
if(iterRes==resRow.end())
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);
+ PlanarIntersectorP1P1(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double md3DSurf, double medianPlane, bool doRotate, int orientation, int printLevel);
public:
void intersectCells(ConnType icellT, const std::vector<ConnType>& icellsS, MyMatrix& res);
int getNumberOfRowsOfResMatrix() const;
{
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)
+ double dimCaracteristic, double precision, double md3DSurf, double medianPlane,
+ bool doRotate, int orientation, int printLevel):
+ PlanarIntersector<MyMeshType,MyMatrix>(meshT,meshS,dimCaracteristic,precision,md3DSurf,medianPlane,doRotate,orientation,printLevel)
{
}
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 ))
+ surf=PlanarIntersector<MyMeshType,MyMatrix>::getValueRegardingOption(surf);
+ if(surf!=0.)
{
typename MyMatrix::value_type::const_iterator iterRes=resRow.find(OTT<ConnType,numPol>::indFC(curNodeSInCmode));
if(iterRes==resRow.end())
static const NumberingPolicy numPol=MyMeshType::My_numPol;
public:
TriangulationIntersector(const MyMeshType& meshT, const MyMeshType& meshS,
- double dimCaracteristic, double precision, double medianPlane, int orientation, int printLevel);
+ double dimCaracteristic, double precision, double md3DSurf, 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);
{
template<class MyMeshType, class MyMatrix, template <class MeshType, class TheMatrix, class ThisIntersector> class InterpType>
TriangulationIntersector<MyMeshType,MyMatrix,InterpType>::TriangulationIntersector(const MyMeshType& meshT, const MyMeshType& meshS,
- double DimCaracteristic, double Precision,
+ double DimCaracteristic, double Precision, double md3DSurf,
double MedianPlane, int orientation, int PrintLevel)
- :InterpType<MyMeshType,MyMatrix,TriangulationIntersector<MyMeshType,MyMatrix,InterpType> >(meshT,meshS,DimCaracteristic, Precision, MedianPlane, true, orientation, PrintLevel)
+ :InterpType<MyMeshType,MyMatrix,TriangulationIntersector<MyMeshType,MyMatrix,InterpType> >(meshT,meshS,DimCaracteristic, Precision, md3DSurf, MedianPlane, true, orientation, PrintLevel)
{
if(PlanarIntersector<MyMeshType,MyMatrix>::_print_level >= 1)
{
_mesh->incrRef();
}
}
+
+
+MEDCouplingMesh *MEDCouplingField::buildSubMeshData(const int *start, const int *end, DataArrayInt *&di) const
+{
+ return _type->buildSubMeshData(start,end,_mesh,di);
+}
namespace ParaMEDMEM
{
+ class DataArrayInt;
class MEDCouplingMesh;
class MEDCouplingFieldDiscretization;
void setDescription(const char *desc) { _desc=desc; }
const char *getName() const { return _name.c_str(); }
TypeOfField getEntity() const;
+ MEDCouplingMesh *buildSubMeshData(const int *start, const int *end, DataArrayInt *&di) const;
protected:
void updateTime();
protected:
#include "MEDCouplingMesh.hxx"
#include "MEDCouplingFieldDouble.hxx"
+#include "MEDCouplingPointSet.hxx"
+
using namespace ParaMEDMEM;
const char MEDCouplingFieldDiscretizationP0::REPR[]="P0";
return mesh->getMeasureField();
}
+/*!
+ * This method returns a submesh of 'mesh' instance constituting cell ids contained in array defined as an interval [start;end).
+ * @ param di is an array returned that specifies entity ids (here cells ids) in mesh 'mesh' of entity in returned submesh.
+ * Example : The first cell id of returned mesh has the (*di)[0] id in 'mesh'
+ */
+MEDCouplingMesh *MEDCouplingFieldDiscretizationP0::buildSubMeshData(const int *start, const int *end, const MEDCouplingMesh *mesh, DataArrayInt *&di) const
+{
+ MEDCouplingPointSet* ret=((const MEDCouplingPointSet *) mesh)->buildPartOfMySelf(start,end,false);
+ di=DataArrayInt::New();
+ di->alloc(end-start,1);
+ int *pt=di->getPointer();
+ std::copy(start,end,pt);
+ return ret;
+}
+
TypeOfField MEDCouplingFieldDiscretizationP1::getEnum() const
{
return TYPE;
//Dual mesh to build
return 0;
}
+
+/*!
+ * This method invert array 'di' that is a conversion map from Old to New node numbering to New to Old node numbering.
+ */
+DataArrayInt *MEDCouplingFieldDiscretizationP1::invertArrayO2N2N2O(const MEDCouplingMesh *mesh, const DataArrayInt *di)
+{
+ DataArrayInt *ret=DataArrayInt::New();
+ ret->alloc(mesh->getNumberOfNodes(),1);
+ int nbOfOldNodes=di->getNumberOfTuples();
+ const int *old2New=di->getConstPointer();
+ int *pt=ret->getPointer();
+ for(int i=0;i!=nbOfOldNodes;i++)
+ if(old2New[i]!=-1)
+ pt[old2New[i]]=i;
+ return ret;
+}
+
+MEDCouplingMesh *MEDCouplingFieldDiscretizationP1::buildSubMeshData(const int *start, const int *end, const MEDCouplingMesh *mesh, DataArrayInt *&di) const
+{
+ MEDCouplingPointSet* ret=((const MEDCouplingPointSet *) mesh)->buildPartOfMySelf(start,end,true);
+ DataArrayInt *diInv=ret->zipCoordsTraducer();
+ di=invertArrayO2N2N2O(ret,diInv);
+ diInv->decrRef();
+ return ret;
+}
namespace ParaMEDMEM
{
+ class DataArrayInt;
class MEDCouplingMesh;
class DataArrayDouble;
class MEDCouplingFieldDouble;
virtual void checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception) = 0;
virtual void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception) = 0;
virtual MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh) const = 0;
+ virtual MEDCouplingMesh *buildSubMeshData(const int *start, const int *end, const MEDCouplingMesh *mesh, DataArrayInt *&di) const = 0;
};
class MEDCOUPLING_EXPORT MEDCouplingFieldDiscretizationP0 : public MEDCouplingFieldDiscretization
void checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception);
void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception);
MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh) const;
+ MEDCouplingMesh *buildSubMeshData(const int *start, const int *end, const MEDCouplingMesh *mesh, DataArrayInt *&di) const;
public:
static const char REPR[];
static const TypeOfField TYPE;
void checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception);
void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception);
MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh) const;
+ MEDCouplingMesh *buildSubMeshData(const int *start, const int *end, const MEDCouplingMesh *mesh, DataArrayInt *&di) const;
+ static DataArrayInt *invertArrayO2N2N2O(const MEDCouplingMesh *mesh, const DataArrayInt *di);
public:
static const char REPR[];
static const TypeOfField TYPE;
#include "MEDCouplingUMeshDesc.hxx"
#include "MEDCouplingMemArray.hxx"
+#include <cmath>
+#include <limits>
+#include <numeric>
+
using namespace ParaMEDMEM;
MEDCouplingPointSet::MEDCouplingPointSet():_coords(0)
}
}
+void MEDCouplingPointSet::zipCoords()
+{
+ checkFullyDefined();
+ DataArrayInt *traducer=zipCoordsTraducer();
+ traducer->decrRef();
+}
+
+void MEDCouplingPointSet::rotate(const double *center, const double *vector, double angle)
+{
+ int spaceDim=getSpaceDimension();
+ if(spaceDim==3)
+ rotate3D(center,vector,angle);
+ else if(spaceDim==2)
+ rotate2D(center,angle);
+ else
+ throw INTERP_KERNEL::Exception("MEDCouplingPointSet::rotate : invalid space dim for rotation must be 2 or 3");
+ _coords->declareAsNew();
+ updateTime();
+}
+
+void MEDCouplingPointSet::translate(const double *vector)
+{
+ double *coords=_coords->getPointer();
+ int nbNodes=getNumberOfNodes();
+ int dim=getSpaceDimension();
+ for(int i=0; i<nbNodes; i++)
+ for(int idim=0; idim<dim;idim++)
+ coords[i*dim+idim]+=vector[idim];
+ _coords->declareAsNew();
+ updateTime();
+}
+
MEDCouplingPointSet *MEDCouplingPointSet::buildInstanceFromMeshType(MEDCouplingMeshType type)
{
switch(type)
{
bool intersects = (bbtemp[idim*2]<bb2[idim*2+1])
&& (bb2[idim*2]<bbtemp[idim*2+1]) ;
- if (!intersects) return false;
+ if (!intersects)
+ return false;
}
return true;
}
+/*!
+ * 'This' is expected to be of spaceDim==3. Idem for 'center' and 'vect'
+ */
+void MEDCouplingPointSet::rotate3D(const double *center, const double *vect, double angle)
+{
+ double sina=sin(angle);
+ double cosa=cos(angle);
+ double vectorNorm[3];
+ double matrix[9];
+ double matrixTmp[9];
+ double norm=sqrt(vect[0]*vect[0]+vect[1]*vect[1]+vect[2]*vect[2]);
+ std::transform(vect,vect+3,vectorNorm,std::bind2nd(std::multiplies<double>(),1/norm));
+ //rotation matrix computation
+ matrix[0]=cosa; matrix[1]=0.; matrix[2]=0.; matrix[3]=0.; matrix[4]=cosa; matrix[5]=0.; matrix[6]=0.; matrix[7]=0.; matrix[8]=cosa;
+ matrixTmp[0]=vectorNorm[0]*vectorNorm[0]; matrixTmp[1]=vectorNorm[0]*vectorNorm[1]; matrixTmp[2]=vectorNorm[0]*vectorNorm[2];
+ matrixTmp[3]=vectorNorm[1]*vectorNorm[0]; matrixTmp[4]=vectorNorm[1]*vectorNorm[1]; matrixTmp[5]=vectorNorm[1]*vectorNorm[2];
+ matrixTmp[6]=vectorNorm[2]*vectorNorm[0]; matrixTmp[7]=vectorNorm[2]*vectorNorm[1]; matrixTmp[8]=vectorNorm[2]*vectorNorm[2];
+ std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies<double>(),1-cosa));
+ std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
+ matrixTmp[0]=0.; matrixTmp[1]=-vectorNorm[2]; matrixTmp[2]=vectorNorm[1];
+ matrixTmp[3]=vectorNorm[2]; matrixTmp[4]=0.; matrixTmp[5]=-vectorNorm[0];
+ matrixTmp[6]=-vectorNorm[1]; matrixTmp[7]=vectorNorm[0]; matrixTmp[8]=0.;
+ std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies<double>(),sina));
+ std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
+ //rotation matrix computed.
+ double *coords=_coords->getPointer();
+ int nbNodes=getNumberOfNodes();
+ double tmp[3];
+ for(int i=0; i<nbNodes; i++)
+ {
+ std::transform(coords+i*3,coords+(i+1)*3,center,tmp,std::minus<double>());
+ coords[i*3]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+matrix[2]*tmp[2]+center[0];
+ coords[i*3+1]=matrix[3]*tmp[0]+matrix[4]*tmp[1]+matrix[5]*tmp[2]+center[1];
+ coords[i*3+2]=matrix[6]*tmp[0]+matrix[7]*tmp[1]+matrix[8]*tmp[2]+center[2];
+ }
+}
+
+/*!
+ * 'This' is expected to be of spaceDim==2. Idem for 'center' and 'vect'
+ */
+void MEDCouplingPointSet::rotate2D(const double *center, double angle)
+{
+ double cosa=cos(angle);
+ double sina=sin(angle);
+ double matrix[4];
+ matrix[0]=cosa; matrix[1]=-sina; matrix[2]=sina; matrix[3]=cosa;
+ double *coords=_coords->getPointer();
+ int nbNodes=getNumberOfNodes();
+ double tmp[2];
+ for(int i=0; i<nbNodes; i++)
+ {
+ std::transform(coords+i*2,coords+(i+1)*2,center,tmp,std::minus<double>());
+ coords[i*2]=matrix[0]*tmp[0]+matrix[1]*tmp[1]+center[0];
+ coords[i*2+1]=matrix[2]*tmp[0]+matrix[3]*tmp[1]+center[1];
+ }
+}
DataArrayDouble *getCoords() const { return _coords; }
bool areCoordsEqual(const MEDCouplingPointSet& other, double prec) const;
void getBoundingBox(double *bbox) const;
+ void zipCoords();
+ void rotate(const double *center, const double *vector, double angle);
+ void translate(const double *vector);
static MEDCouplingPointSet *buildInstanceFromMeshType(MEDCouplingMeshType type);
virtual MEDCouplingPointSet *buildPartOfMySelf(const int *start, const int *end, bool keepCoords) const = 0;
+ virtual MEDCouplingPointSet *buildPartOfMySelfNode(const int *start, const int *end, bool fullyIn) const = 0;
+ virtual void renumberConnectivity(const int *newNodeNumbers) = 0;
//! size of returned tinyInfo must be always the same.
virtual void getTinySerializationInformation(std::vector<int>& tinyInfo, std::vector<std::string>& littleStrings) const;
virtual void resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector<std::string>& littleStrings);
virtual void unserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2,
const std::vector<std::string>& littleStrings);
virtual void giveElemsInBoundingBox(const double *bbox, double eps, std::vector<int>& elems) = 0;
+ virtual DataArrayInt *zipCoordsTraducer() = 0;
protected:
+ virtual void checkFullyDefined() const throw(INTERP_KERNEL::Exception) = 0;
static bool intersectsBoundingBox(const double* bb1, const double* bb2, int dim, double eps);
+ void rotate2D(const double *center, double angle);
+ void rotate3D(const double *center, const double *vect, double angle);
protected:
DataArrayDouble *_coords;
};
#include "VolSurfFormulae.hxx"
#include <sstream>
+#include <limits>
using namespace ParaMEDMEM;
}
}
-void MEDCouplingUMesh::zipCoords()
-{
- checkFullyDefined();
- DataArrayInt *traducer=zipCoordsTraducer();
- traducer->decrRef();
-}
-
struct MEDCouplingAccVisit
{
MEDCouplingAccVisit():_new_nb_of_nodes(0) { }
return ret;
}
+/*!
+ * Keeps from 'this' only cells which constituing point id is in the ids specified by 'start' and 'end'.
+ * The return newly allocated mesh will share the same coordinates as 'this'.
+ * Parameter 'fullyIn' specifies if a cell that has part of its nodes in ids array is kept or not.
+ * If 'fullyIn' is true only cells whose ids are \b fully contained in ['start','end') tab will be kept.
+ */
+MEDCouplingPointSet *MEDCouplingUMesh::buildPartOfMySelfNode(const int *start, const int *end, bool fullyIn) const
+{
+ std::set<int> fastFinder(start,end);
+ const int *conn=getNodalConnectivity()->getConstPointer();
+ const int *connIndex=getNodalConnectivityIndex()->getConstPointer();
+ int nbOfCells=getNumberOfCells();
+ std::vector<int> cellIdsKept;
+ for(int i=0;i<nbOfCells;i++)
+ {
+ std::set<int> connOfCell(conn+connIndex[i]+1,conn+connIndex[i+1]);
+ connOfCell.erase(-1);//polyhedron separator
+ int refLgth=std::min(connOfCell.size(),fastFinder.size());
+ std::set<int> locMerge;
+ std::insert_iterator< std::set<int> > it(locMerge,locMerge.begin());
+ std::set_intersection(connOfCell.begin(),connOfCell.end(),fastFinder.begin(),fastFinder.end(),it);
+ if(locMerge.size()==refLgth && fullyIn || locMerge.size()!=0 && !fullyIn)
+ cellIdsKept.push_back(i);
+ }
+ return buildPartOfMySelf(&cellIdsKept[0],&cellIdsKept[cellIdsKept.size()],true);
+}
+
+void MEDCouplingUMesh::renumberConnectivity(const int *newNodeNumbers)
+{
+ int *conn=getNodalConnectivity()->getPointer();
+ const int *connIndex=getNodalConnectivityIndex()->getConstPointer();
+ int nbOfCells=getNumberOfCells();
+ for(int i=0;i<nbOfCells;i++)
+ for(int iconn=connIndex[i]+1;iconn!=connIndex[i+1];iconn++)
+ {
+ int& node=conn[iconn];
+ if(node>=0)//avoid polyhedron separator
+ {
+ node=newNodeNumbers[node];
+ }
+ }
+}
+
/*!
* 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
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(node>=0)//avoid polyhedron separator
{
- 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] )
+ for (int idim=0; idim<dim; idim++)
{
- elem_bb[idim*2+1] = coords[node*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] ;
+ }
}
}
}
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);
+ area_vol[iel]=fabs(INTERP_KERNEL::calculateAreaForTria(coords+(dim_space*N1),
+ coords+(dim_space*N2),
+ coords+(dim_space*N3),
+ dim_space));
}
break ;
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) ;
+ area_vol[iel]=fabs(INTERP_KERNEL::calculateAreaForQuad(coords+dim_space*N1,
+ coords+dim_space*N2,
+ coords+dim_space*N3,
+ coords+dim_space*N4,
+ dim_space));
}
break ;
void unserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2,
const std::vector<std::string>& littleStrings);
//tools
- void zipCoords();
DataArrayInt *zipCoordsTraducer();
void getReverseNodalConnectivity(DataArrayInt *revNodal, DataArrayInt *revNodalIndx) const;
MEDCouplingPointSet *buildPartOfMySelf(const int *start, const int *end, bool keepCoords) const;
+ MEDCouplingPointSet *buildPartOfMySelfNode(const int *start, const int *end, bool fullyIn) const;
+ void renumberConnectivity(const int *newNodeNumbers);
void giveElemsInBoundingBox(const double *bbox, double eps, std::vector<int>& elems);
MEDCouplingFieldDouble *getMeasureField() const;
private:
#include "CellModel.hxx"
#include "MEDCouplingMemArray.hxx"
+#include <limits>
+
using namespace ParaMEDMEM;
MEDCouplingUMeshDesc::MEDCouplingUMeshDesc():_mesh_dim(-1),_desc_connec(0),_desc_connec_index(0),
return 0;
}
+MEDCouplingPointSet *MEDCouplingUMeshDesc::buildPartOfMySelfNode(const int *start, const int *end, bool fullyIn) const
+{
+ //not implemented yet
+ return 0;
+}
+
+void MEDCouplingUMeshDesc::renumberConnectivity(const int *newNodeNumbers)
+{
+}
+
MEDCouplingFieldDouble *MEDCouplingUMeshDesc::getMeasureField() const
{
//not implemented yet.
return 0;
}
+DataArrayInt *MEDCouplingUMeshDesc::zipCoordsTraducer()
+{
+ //not implemented yet.
+ return 0;
+}
+
void MEDCouplingUMeshDesc::computeTypes()
{
if(_desc_connec && _desc_connec_index)
_types.insert((INTERP_KERNEL::NormalizedCellType)conn[*pt]);
}
}
+
+void MEDCouplingUMeshDesc::checkFullyDefined() const throw(INTERP_KERNEL::Exception)
+{
+ if(!_desc_connec || !_desc_connec_index || !_nodal_connec_face || !_nodal_connec_face_index || !_coords)
+ throw INTERP_KERNEL::Exception("full connectivity and coordinates not set in unstructured mesh.");
+}
const std::vector<std::string>& littleStrings);
void giveElemsInBoundingBox(const double *bbox, double eps, std::vector<int>& elems);
MEDCouplingPointSet *buildPartOfMySelf(const int *start, const int *end, bool keepCoords) const;
+ MEDCouplingPointSet *buildPartOfMySelfNode(const int *start, const int *end, bool fullyIn) const;
+ void renumberConnectivity(const int *newNodeNumbers);
MEDCouplingFieldDouble *getMeasureField() const;
+ DataArrayInt *zipCoordsTraducer();
private:
MEDCouplingUMeshDesc();
~MEDCouplingUMeshDesc();
void computeTypes();
+ void checkFullyDefined() const throw(INTERP_KERNEL::Exception);
private:
unsigned _mesh_dim;
DataArrayInt *_desc_connec;
mesh->decrRef();
}
+void MEDCouplingBasicsTest::testBuildPartOfMySelfNode()
+{
+ MEDCouplingUMesh *mesh=build2DTargetMesh_1();
+ const int tab1[2]={5,7};
+ MEDCouplingPointSet *subMeshSimple=mesh->buildPartOfMySelfNode(tab1,tab1+2,true);
+ MEDCouplingUMesh *subMesh=dynamic_cast<MEDCouplingUMesh *>(subMeshSimple);
+ CPPUNIT_ASSERT(subMesh);
+ CPPUNIT_ASSERT_EQUAL(1,(int)subMesh->getAllTypes().size());
+ CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_QUAD4,*subMesh->getAllTypes().begin());
+ CPPUNIT_ASSERT_EQUAL(1,subMesh->getNumberOfCells());
+ CPPUNIT_ASSERT_EQUAL(5,subMesh->getNodalConnectivity()->getNbOfElems());
+ CPPUNIT_ASSERT_EQUAL(2,subMesh->getNodalConnectivityIndex()->getNbOfElems());
+ const int subConn[5]={4,7,8,5,4};
+ const int subConnIndex[3]={0,5};
+ CPPUNIT_ASSERT(std::equal(subConn,subConn+5,subMesh->getNodalConnectivity()->getPointer()));
+ CPPUNIT_ASSERT(std::equal(subConnIndex,subConnIndex+2,subMesh->getNodalConnectivityIndex()->getPointer()));
+ CPPUNIT_ASSERT(subMesh->getCoords()==mesh->getCoords());
+ subMeshSimple->decrRef();
+ //
+ subMeshSimple=mesh->buildPartOfMySelfNode(tab1,tab1+2,false);
+ subMesh=dynamic_cast<MEDCouplingUMesh *>(subMeshSimple);
+ CPPUNIT_ASSERT(subMesh);
+ CPPUNIT_ASSERT_EQUAL(2,(int)subMesh->getAllTypes().size());
+ CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_TRI3,*subMesh->getAllTypes().begin());
+ CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_QUAD4,*(++subMesh->getAllTypes().begin()));
+ CPPUNIT_ASSERT_EQUAL(3,subMesh->getNumberOfCells());
+ CPPUNIT_ASSERT_EQUAL(14,subMesh->getNodalConnectivity()->getNbOfElems());
+ CPPUNIT_ASSERT_EQUAL(4,subMesh->getNodalConnectivityIndex()->getNbOfElems());
+ const int subConn2[14]={3,4,5,2,4,6,7,4,3,4,7,8,5,4};
+ const int subConnIndex2[4]={0,4,9,14};
+ CPPUNIT_ASSERT(std::equal(subConn2,subConn2+14,subMesh->getNodalConnectivity()->getPointer()));
+ CPPUNIT_ASSERT(std::equal(subConnIndex2,subConnIndex2+4,subMesh->getNodalConnectivityIndex()->getPointer()));
+ CPPUNIT_ASSERT(subMesh->getCoords()==mesh->getCoords());
+ subMeshSimple->decrRef();
+ //testing the case where length of tab2 is greater than max number of node per cell.
+ const int tab2[7]={0,3,2,1,4,5,6};
+ subMeshSimple=mesh->buildPartOfMySelfNode(tab2,tab2+7,true);
+ subMesh=dynamic_cast<MEDCouplingUMesh *>(subMeshSimple);
+ CPPUNIT_ASSERT(subMesh);
+ CPPUNIT_ASSERT_EQUAL(2,(int)subMesh->getAllTypes().size());
+ CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_TRI3,*subMesh->getAllTypes().begin());
+ CPPUNIT_ASSERT_EQUAL(INTERP_KERNEL::NORM_QUAD4,*(++subMesh->getAllTypes().begin()));
+ CPPUNIT_ASSERT_EQUAL(3,subMesh->getNumberOfCells());
+ subMeshSimple->decrRef();
+ //
+ mesh->decrRef();
+}
+
void MEDCouplingBasicsTest::testZipCoords()
{
MEDCouplingUMesh *mesh=build2DTargetMesh_1();
targetMesh->decrRef();
}
+void MEDCouplingBasicsTest::test3DSurfInterpP1P1_1()
+{
+ MEDCouplingUMesh *sourceMesh=build3DSurfSourceMesh_1();
+ MEDCouplingUMesh *targetMesh=build3DSurfTargetMesh_2();
+ //
+ MEDCouplingNormalizedUnstructuredMesh<3,2> sourceWrapper(sourceMesh);
+ MEDCouplingNormalizedUnstructuredMesh<3,2> targetWrapper(targetMesh);
+ INTERP_KERNEL::Interpolation3DSurf 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*sqrt(2.),res[0][0],1.e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.05416666666666665*sqrt(2.),res[1][0],1.e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.02916666666666666*sqrt(2.),res[1][1],1.e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.08333333333333334*sqrt(2.),res[2][1],1.e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.05416666666666665*sqrt(2.),res[3][0],1.e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.02916666666666668*sqrt(2.),res[3][2],1.e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.1416666666666666*sqrt(2.),res[4][0],1.e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.02499999999999999*sqrt(2.),res[4][1],1.e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.02499999999999999*sqrt(2.),res[4][2],1.e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.09999999999999999*sqrt(2.),res[4][3],1.e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.02916666666666666*sqrt(2.),res[5][1],1.e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.09583333333333333*sqrt(2.),res[5][3],1.e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.08333333333333333*sqrt(2.),res[6][2],1.e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.02916666666666667*sqrt(2.),res[7][2],1.e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.09583333333333331*sqrt(2.),res[7][3],1.e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.04166666666666668*sqrt(2.),res[8][3],1.e-12);
+ res.clear();
+ }
+ //
+ sourceMesh->decrRef();
+ targetMesh->decrRef();
+}
+
+void MEDCouplingBasicsTest::test3DSurfInterpP0P0_2()
+{
+ MEDCouplingUMesh *sourceMesh=build3DSurfSourceMesh_1();
+ MEDCouplingUMesh *targetMesh=build3DSurfTargetMeshPerm_1();
+ //
+ MEDCouplingNormalizedUnstructuredMesh<3,2> sourceWrapper(sourceMesh);
+ MEDCouplingNormalizedUnstructuredMesh<3,2> targetWrapper(targetMesh);
+ INTERP_KERNEL::Interpolation3DSurf myInterpolator;
+ vector<map<int,double> > res;
+ myInterpolator.setPrecision(1e-12);
+ myInterpolator.setIntersectionType(INTERP_KERNEL::Triangulation);
+ {
+ myInterpolator.setOrientation(2);
+ myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,res,"P0P0");
+ CPPUNIT_ASSERT_EQUAL(5,(int)res.size());
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[0][0],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[0][1],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[1][0],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[2][0],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25*sqrt(2.),res[3][1],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[4][0],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[4][1],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.*sqrt(2.),sumAll(res),1e-12);
+ res.clear();
+ }
+ {
+ myInterpolator.setOrientation(0);
+ myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,res,"P0P0");
+ CPPUNIT_ASSERT_EQUAL(5,(int)res.size());
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[0][0],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[0][1],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[1][0],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.125*sqrt(2.),res[2][0],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25*sqrt(2.),res[3][1],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[4][0],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[4][1],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.75*sqrt(2.),sumAll(res),1e-12);
+ res.clear();
+ }
+ {
+ myInterpolator.setOrientation(1);
+ myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,res,"P0P0");
+ CPPUNIT_ASSERT_EQUAL(5,(int)res.size());
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[0][0],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[0][1],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[1][0],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25*sqrt(2.),res[3][1],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[4][0],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[4][1],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.875*sqrt(2.),sumAll(res),1e-12);
+ res.clear();
+ }
+ {
+ myInterpolator.setOrientation(-1);
+ myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,res,"P0P0");
+ CPPUNIT_ASSERT_EQUAL(5,(int)res.size());
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.125*sqrt(2.),res[2][0],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.125*sqrt(2.),sumAll(res),1e-12);
+ res.clear();
+ }
+ //clean up
+ sourceMesh->decrRef();
+ targetMesh->decrRef();
+}
+
+/*!
+ * Test of precision option implemented by Fabien that represents distance of "barycenter" to the other cell.
+ */
+void MEDCouplingBasicsTest::test3DSurfInterpP0P0_3()
+{
+ INTERP_KERNEL::Interpolation3DSurf myInterpolator;
+ vector<map<int,double> > res;
+ double vecTrans[3]={0.,0.,1.e-10};
+ double vec[3]={0.,-1.,0.};
+ double pt[3]={-0.3,-0.3,5.e-11};
+ const int N=32;
+ const double deltaA=M_PI/N;
+ myInterpolator.setPrecision(1e-12);
+ myInterpolator.setIntersectionType(INTERP_KERNEL::Triangulation);
+ myInterpolator.setMaxDistance3DSurfIntersect(1e-9);
+ for(int i=0;i<N;i++)
+ {
+ res.clear();
+ MEDCouplingUMesh *sourceMesh=build3DSurfSourceMesh_2();
+ sourceMesh->rotate(pt,vec,i*deltaA);
+ MEDCouplingUMesh *targetMesh=build3DSurfSourceMesh_2();
+ targetMesh->translate(vecTrans);
+ targetMesh->rotate(pt,vec,i*deltaA);
+ MEDCouplingNormalizedUnstructuredMesh<3,2> sourceWrapper(sourceMesh);
+ MEDCouplingNormalizedUnstructuredMesh<3,2> targetWrapper(targetMesh);
+ myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,res,"P0P0");
+ CPPUNIT_ASSERT_EQUAL(2,(int)res.size());
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,sumAll(res),1e-12);
+ sourceMesh->decrRef();
+ targetMesh->decrRef();
+ }
+ //
+ myInterpolator.setMaxDistance3DSurfIntersect(1e-11);
+ for(int i=0;i<N;i++)
+ {
+ res.clear();
+ MEDCouplingUMesh *sourceMesh=build3DSurfSourceMesh_2();
+ sourceMesh->rotate(pt,vec,i*deltaA);
+ MEDCouplingUMesh *targetMesh=build3DSurfSourceMesh_2();
+ targetMesh->translate(vecTrans);
+ targetMesh->rotate(pt,vec,i*deltaA);
+ MEDCouplingNormalizedUnstructuredMesh<3,2> sourceWrapper(sourceMesh);
+ MEDCouplingNormalizedUnstructuredMesh<3,2> targetWrapper(targetMesh);
+ myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,res,"P0P0");
+ CPPUNIT_ASSERT_EQUAL(2,(int)res.size());
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.,sumAll(res),1e-12);
+ sourceMesh->decrRef();
+ targetMesh->decrRef();
+ }
+ //
+ res.clear();
+ myInterpolator.setMaxDistance3DSurfIntersect(-1.);//unactivate fabien lookup
+ MEDCouplingUMesh *sourceMesh=build3DSurfSourceMesh_2();
+ MEDCouplingUMesh *targetMesh=build3DSurfSourceMesh_2();
+ targetMesh->translate(vecTrans);
+ myInterpolator.setBoundingBoxAdjustment(1e-11);
+ MEDCouplingNormalizedUnstructuredMesh<3,2> sourceWrapper0(sourceMesh);
+ MEDCouplingNormalizedUnstructuredMesh<3,2> targetWrapper0(targetMesh);
+ myInterpolator.interpolateMeshes(sourceWrapper0,targetWrapper0,res,"P0P0");
+ CPPUNIT_ASSERT_EQUAL(2,(int)res.size());
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.,sumAll(res),1e-12);
+ sourceMesh->decrRef();
+ targetMesh->decrRef();
+ //
+ res.clear();
+ sourceMesh=build3DSurfSourceMesh_2();
+ targetMesh=build3DSurfSourceMesh_2();
+ targetMesh->translate(vecTrans);
+ myInterpolator.setBoundingBoxAdjustment(1e-9);
+ MEDCouplingNormalizedUnstructuredMesh<3,2> sourceWrapper1(sourceMesh);
+ MEDCouplingNormalizedUnstructuredMesh<3,2> targetWrapper1(targetMesh);
+ myInterpolator.interpolateMeshes(sourceWrapper1,targetWrapper1,res,"P0P0");
+ CPPUNIT_ASSERT_EQUAL(2,(int)res.size());
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,sumAll(res),1e-12);
+ sourceMesh->decrRef();
+ targetMesh->decrRef();
+ //keeping the same bbox adj == 1.e-11 but trying rotation
+ res.clear();
+ sourceMesh=build3DSurfSourceMesh_2();
+ sourceMesh->rotate(pt,vec,M_PI/4.);
+ targetMesh=build3DSurfSourceMesh_2();
+ targetMesh->translate(vecTrans);
+ targetMesh->rotate(pt,vec,M_PI/4.);
+ myInterpolator.setBoundingBoxAdjustment(1e-11);
+ MEDCouplingNormalizedUnstructuredMesh<3,2> sourceWrapper2(sourceMesh);
+ MEDCouplingNormalizedUnstructuredMesh<3,2> targetWrapper2(targetMesh);
+ myInterpolator.interpolateMeshes(sourceWrapper2,targetWrapper2,res,"P0P0");
+ CPPUNIT_ASSERT_EQUAL(2,(int)res.size());
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,sumAll(res),1e-12);
+ sourceMesh->decrRef();
+ targetMesh->decrRef();
+}
+
void MEDCouplingBasicsTest::test3DInterpP0P0_1()
{
MEDCouplingUMesh *sourceMesh=build3DSourceMesh_1();
return sourceMesh;
}
+MEDCouplingUMesh *MEDCouplingBasicsTest::build3DSurfSourceMesh_2()
+{
+ double sourceCoords[12]={-0.3,-0.3,0., 0.7,-0.3,0., -0.3,0.7,0., 0.7,0.7,0.};
+ int sourceConn[6]={0,3,1,0,2,3};
+ MEDCouplingUMesh *sourceMesh=MEDCouplingUMesh::New();
+ sourceMesh->setMeshDimension(2);
+ sourceMesh->allocateCells(2);
+ sourceMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,sourceConn);
+ sourceMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,sourceConn+3);
+ sourceMesh->finishInsertingCells();
+ DataArrayDouble *myCoords=DataArrayDouble::New();
+ myCoords->alloc(4,3);
+ std::copy(sourceCoords,sourceCoords+12,myCoords->getPointer());
+ sourceMesh->setCoords(myCoords);
+ myCoords->decrRef();
+ return sourceMesh;
+}
+
MEDCouplingUMesh *MEDCouplingBasicsTest::build3DSurfTargetMesh_1()
{
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};
return targetMesh;
}
+/*!
+ * Idem build3DSurfTargetMesh_1 except that cell id 2 is not correctly numbered.
+ */
+MEDCouplingUMesh *MEDCouplingBasicsTest::build3DSurfTargetMeshPerm_1()
+{
+ 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[18]={0,3,4,1, 1,4,2, 4,2,5, 6,7,4,3, 7,8,5,4};
+ MEDCouplingUMesh *targetMesh=MEDCouplingUMesh::New();
+ targetMesh->setMeshDimension(2);
+ targetMesh->allocateCells(5);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,targetConn);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+4);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+7);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,targetConn+10);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,targetConn+14);
+ 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::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};
CPPUNIT_TEST( testDeepCopy );
CPPUNIT_TEST( testRevNodal );
CPPUNIT_TEST( testBuildPartOfMySelf );
+ CPPUNIT_TEST( testBuildPartOfMySelfNode );
CPPUNIT_TEST( testZipCoords );
CPPUNIT_TEST( testEqualMesh );
CPPUNIT_TEST( testEqualFieldDouble );
CPPUNIT_TEST( test3DSurfInterpP0P0_1 );
CPPUNIT_TEST( test3DSurfInterpP0P1_1 );
CPPUNIT_TEST( test3DSurfInterpP1P0_1 );
+ CPPUNIT_TEST( test3DSurfInterpP1P1_1 );
+ CPPUNIT_TEST( test3DSurfInterpP0P0_2 );
+ CPPUNIT_TEST( test3DSurfInterpP0P0_3 );
CPPUNIT_TEST( test3DInterpP0P0_1 );
CPPUNIT_TEST( test3DInterpP0P1_1 );
CPPUNIT_TEST( test3DInterpP1P0_1 );
void testDeepCopy();
void testRevNodal();
void testBuildPartOfMySelf();
+ void testBuildPartOfMySelfNode();
void testZipCoords();
void testEqualMesh();
void testEqualFieldDouble();
void test3DSurfInterpP0P0_1();
void test3DSurfInterpP0P1_1();
void test3DSurfInterpP1P0_1();
+ void test3DSurfInterpP1P1_1();
+ void test3DSurfInterpP0P0_2();
+ void test3DSurfInterpP0P0_3();
void test3DInterpP0P0_1();
void test3DInterpP0P1_1();
void test3DInterpP1P0_1();
MEDCouplingUMesh *build2DTargetMesh_1();
MEDCouplingUMesh *build2DTargetMesh_2();
MEDCouplingUMesh *build3DSurfSourceMesh_1();
+ MEDCouplingUMesh *build3DSurfSourceMesh_2();
MEDCouplingUMesh *build3DSurfTargetMesh_1();
+ MEDCouplingUMesh *build3DSurfTargetMeshPerm_1();
MEDCouplingUMesh *build3DSurfTargetMesh_2();
MEDCouplingUMesh *build3DSourceMesh_1();
MEDCouplingUMesh *build3DTargetMesh_1();
delete _local_field;
delete _icoco_field;
delete _union_group;
- }
+ }
+
+ void DEC::setNature(NatureOfField nature)
+ {
+ if(_local_field)
+ _local_field->getField()->setNature(nature);
+ }
+
/*! Attaches a local field to a DEC.
If the processor is on the receiving end of the DEC, the field
will be updated by a recvData() call.
localgroup=_source_group;
else
localgroup=_target_group;
- //delete _icoco_field;
+ //delete _icoco_field;
_icoco_field=new ICoCo::MEDField(*const_cast<ICoCo::TrioField* >(triofield), *localgroup);
attachLocalField(_icoco_field);
return;
{
public:
DEC():_local_field(0) { }
- DEC(ProcessorGroup& source_group, ProcessorGroup& target_group);
+ DEC(ProcessorGroup& source_group, ProcessorGroup& target_group);
+ void setNature(NatureOfField nature);
void attachLocalField( MEDCouplingFieldDouble* field);
void attachLocalField(const ParaFIELD* field);
void attachLocalField(const ICoCo::Field* field);
#include "ElementLocator.hxx"
#include "Topology.hxx"
#include "BlockTopology.hxx"
+#include "ParaFIELD.hxx"
#include "ParaMESH.hxx"
#include "ProcessorGroup.hxx"
#include "MPIProcessorGroup.hxx"
+#include "MEDCouplingFieldDouble.hxx"
#include <map>
#include <set>
namespace ParaMEDMEM
{
- ElementLocator::ElementLocator(const ParaMESH& sourceMesh,
+ ElementLocator::ElementLocator(const ParaFIELD& sourceField,
const ProcessorGroup& distant_group)
- : _local_para_mesh(sourceMesh),
- _local_cell_mesh(sourceMesh.getCellMesh()),
- _local_face_mesh(sourceMesh.getFaceMesh()),
- _local_group(*sourceMesh.getBlockTopology()->getProcGroup()),
+ : _local_para_field(sourceField),
+ _local_cell_mesh(sourceField.getSupport()->getCellMesh()),
+ _local_face_mesh(sourceField.getSupport()->getFaceMesh()),
+ _local_group(*sourceField.getSupport()->getBlockTopology()->getProcGroup()),
_distant_group(distant_group)
{
_union_group = _local_group.fuse(distant_group);
delete [] _domain_bounding_boxes;
}
+ const MPI_Comm *ElementLocator::getCommunicator() const
+ {
+ MPIProcessorGroup* group=static_cast<MPIProcessorGroup*> (_union_group);
+ return group->getComm();
+ }
+
// ==========================================================================
// Procedure for exchanging mesh between a distant proc and a local processor
// param idistantrank proc id on distant group
vector<int> elems;
double* distant_bb = _domain_bounding_boxes+rank*2*dim;
_local_cell_mesh->giveElemsInBoundingBox(distant_bb,getBoundingBoxAdjustment(),elems);
+
+ DataArrayInt *distant_ids_send;
+ MEDCouplingPointSet *send_mesh = (MEDCouplingPointSet *)_local_para_field.getField()->buildSubMeshData(&elems[0],&elems[elems.size()],distant_ids_send);
+ _exchangeMesh(send_mesh, distant_mesh, idistantrank, distant_ids_send, distant_ids);
+ distant_ids_send->decrRef();
+
+ if(send_mesh)
+ send_mesh->decrRef();
+#if 0
+ int* distant_ids_send=0;
//send_mesh contains null pointer if elems is empty
MEDCouplingPointSet* send_mesh = _local_cell_mesh->buildPartOfMySelf(&elems[0],&elems[elems.size()],false);
// Constituting an array containing the ids of the elements that are
// This array enables the correct redistribution of the data when the
// interpolated field is transmitted to the target array
- int* distant_ids_send=0;
if (elems.size()>0)
{
distant_ids_send = new int[elems.size()];
- int index=0;
- for (std::vector<int>::const_iterator iter = elems.begin(); iter!= elems.end(); iter++)
- {
- distant_ids_send[index]=*iter;
- index++;
- }
+ std::copy(elems.begin(),elems.end(),distant_ids_send);
}
_exchangeMesh(send_mesh, distant_mesh, idistantrank, distant_ids_send, distant_ids);
delete[] distant_ids_send;
if(send_mesh)
send_mesh->decrRef();
+#endif
}
void ElementLocator::exchangeMethod(const std::string& sourceMeth, int idistantrank, std::string& targetMeth)
void ElementLocator::_exchangeMesh( MEDCouplingPointSet* local_mesh,
MEDCouplingPointSet*& distant_mesh,
int iproc_distant,
- const int* distant_ids_send,
+ const DataArrayInt* distant_ids_send,
int*& distant_ids_recv)
{
CommInterface comm_interface=_union_group->getCommInterface();
//Getting tiny info of local mesh to allow the distant proc to initialize and allocate
//the transmitted mesh.
local_mesh->getTinySerializationInformation(tinyInfoLocal,tinyInfoLocalS);
- tinyInfoLocal.push_back(local_mesh->getNumberOfCells());
+ tinyInfoLocal.push_back(distant_ids_send->getNumberOfTuples());
tinyInfoDistant.resize(tinyInfoLocal.size());
std::fill(tinyInfoDistant.begin(),tinyInfoDistant.end(),0);
MPIProcessorGroup* group=static_cast<MPIProcessorGroup*> (_union_group);
else
distant_mesh_tmp->decrRef();
distant_ids_recv=new int[tinyInfoDistant.back()];
- comm_interface.sendRecv((void *)distant_ids_send,tinyInfoLocal.back(), MPI_INT,
+ comm_interface.sendRecv((void *)distant_ids_send->getConstPointer(),tinyInfoLocal.back(), MPI_INT,
iprocdistant_in_union, 1113,
distant_ids_recv,tinyInfoDistant.back(), MPI_INT,
iprocdistant_in_union,1113,
#include "InterpolationOptions.hxx"
+#include <mpi.h>
#include <vector>
#include <set>
namespace ParaMEDMEM
{
- class ParaMESH;
+ class ParaFIELD;
class ProcessorGroup;
class ParaSUPPORT;
class InterpolationMatrix;
class MEDCouplingPointSet;
+ class DataArrayInt;
class ElementLocator : public INTERP_KERNEL::InterpolationOptions
{
public:
- ElementLocator(const ParaMESH& sourceMesh, const ProcessorGroup& distant_group);
+ ElementLocator(const ParaFIELD& sourceField, const ProcessorGroup& distant_group);
virtual ~ElementLocator();
void exchangeMesh(int idistantrank,
MEDCouplingPointSet*& target_mesh,
int*& distant_ids);
void exchangeMethod(const std::string& sourceMeth, int idistantrank, std::string& targetMeth);
+ const std::vector<int>& getDistantProcIds() const { return _distant_proc_ids; }
+ const MPI_Comm *getCommunicator() const;
private:
- const ParaMESH& _local_para_mesh ;
+ void _computeBoundingBoxes();
+ bool _intersectsBoundingBox(int irank);
+ void _exchangeMesh(MEDCouplingPointSet* local_mesh, MEDCouplingPointSet*& distant_mesh,
+ int iproc_distant, const DataArrayInt* distant_ids_send,
+ int*& distant_ids_recv);
+ private:
+ const ParaFIELD& _local_para_field ;
MEDCouplingPointSet* _local_cell_mesh;
MEDCouplingPointSet* _local_face_mesh;
std::vector<MEDCouplingPointSet*> _distant_cell_meshes;
const ProcessorGroup& _local_group;
ProcessorGroup* _union_group;
std::vector<int> _distant_proc_ids;
-
- void _computeBoundingBoxes();
- bool _intersectsBoundingBox(int irank);
- void _exchangeMesh(MEDCouplingPointSet* local_mesh, MEDCouplingPointSet*& distant_mesh,
- int iproc_distant, const int* distant_ids_send,
- int*& distant_ids_recv);
};
}
--- /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 "GlobalizerMesh.hxx"
+#include "MEDCouplingFieldDouble.hxx"
+#include "CommInterface.hxx"
+
+using namespace std;
+
+namespace ParaMEDMEM
+{
+ GlobalizerMesh::GlobalizerMesh(const MPI_Comm *comm, MEDCouplingFieldDouble *localField):_comm(comm),_local_field(localField)
+ {
+ if(_local_field)
+ _local_field->incrRef();
+ }
+
+ GlobalizerMesh::~GlobalizerMesh()
+ {
+ if(_local_field)
+ _local_field->decrRef();
+ }
+
+ NatureOfField GlobalizerMesh::getLocalNature() const
+ {
+ return _local_field->getNature();
+ }
+
+ GlobalizerMeshWorkingSide::GlobalizerMeshWorkingSide(const MPI_Comm *comm, MEDCouplingFieldDouble *localField,
+ const std::string& distantMeth, const std::vector<int>& lazyProcs):GlobalizerMesh(comm,localField),_distant_method(distantMeth),_lazy_procs(lazyProcs)
+ {
+ }
+
+ GlobalizerMeshWorkingSide::~GlobalizerMeshWorkingSide()
+ {
+ }
+
+ const std::vector<int>& GlobalizerMeshWorkingSide::getProcIdsInInteraction() const
+ {
+ return _lazy_procs;
+ }
+
+ /*!
+ * connected with GlobalizerMeshLazySide::recvFromWorkingSide
+ */
+ void GlobalizerMeshWorkingSide::sendSumToLazySide(const std::vector< std::vector<int> >& distantLocEltIds, const std::vector< std::vector<double> >& partialSumRelToDistantIds)
+ {
+ int procId=0;
+ CommInterface comm;
+ for(vector<int>::const_iterator iter=_lazy_procs.begin();iter!=_lazy_procs.end();iter++,procId++)
+ {
+ const vector<int>& eltIds=distantLocEltIds[procId];
+ const vector<double>& valued=partialSumRelToDistantIds[procId];
+ int lgth=eltIds.size();
+ comm.send(&lgth,1,MPI_INT,*iter,1114,*_comm);
+ comm.send((void *)&eltIds[0],lgth,MPI_INT,*iter,1115,*_comm);
+ comm.send((void *)&valued[0],lgth,MPI_DOUBLE,*iter,1116,*_comm);
+ }
+ }
+
+ /*!
+ * connected with GlobalizerMeshLazySide::sendToWorkingSide
+ */
+ void GlobalizerMeshWorkingSide::recvSumFromLazySide(std::vector< std::vector<double> >& globalSumRelToDistantIds)
+ {
+ int procId=0;
+ CommInterface comm;
+ MPI_Status status;
+ for(vector<int>::const_iterator iter=_lazy_procs.begin();iter!=_lazy_procs.end();iter++,procId++)
+ {
+ std::vector<double>& vec=globalSumRelToDistantIds[procId];
+ comm.recv(&vec[0],vec.size(),MPI_DOUBLE,*iter,1117,*_comm,&status);
+ }
+ }
+
+ GlobalizerMeshLazySide::GlobalizerMeshLazySide(const MPI_Comm *comm, MEDCouplingFieldDouble *localField, const std::vector<int>& computeProcs):GlobalizerMesh(comm,localField),_compute_procs(computeProcs)
+ {
+ }
+
+ GlobalizerMeshLazySide::~GlobalizerMeshLazySide()
+ {
+ }
+
+ /*!
+ * connected with GlobalizerMeshWorkingSide::sendSumToLazySide
+ */
+ void GlobalizerMeshLazySide::recvFromWorkingSide()
+ {
+ _values_added.resize(_local_field->getNumberOfTuples());
+ int procId=0;
+ CommInterface comm;
+ _ids_per_working_proc.resize(_compute_procs.size());
+ MPI_Status status;
+ for(vector<int>::const_iterator iter=_compute_procs.begin();iter!=_compute_procs.end();iter++,procId++)
+ {
+ int lgth;
+ comm.recv(&lgth,1,MPI_INT,*iter,1114,*_comm,&status);
+ vector<int>& ids=_ids_per_working_proc[procId];
+ ids.resize(lgth);
+ vector<double> values(lgth);
+ comm.recv(&ids[0],lgth,MPI_INT,*iter,1115,*_comm,&status);
+ comm.recv(&values[0],lgth,MPI_DOUBLE,*iter,1116,*_comm,&status);
+ for(int i=0;i<lgth;i++)
+ _values_added[ids[i]]+=values[i];
+ }
+ }
+
+ /*!
+ * connected with GlobalizerMeshWorkingSide::recvSumFromLazySide
+ */
+ void GlobalizerMeshLazySide::sendToWorkingSide()
+ {
+ int procId=0;
+ CommInterface comm;
+ for(vector<int>::const_iterator iter=_compute_procs.begin();iter!=_compute_procs.end();iter++,procId++)
+ {
+ vector<int>& ids=_ids_per_working_proc[procId];
+ vector<double> valsToSend(ids.size());
+ vector<double>::iterator iter3=valsToSend.begin();
+ for(vector<int>::const_iterator iter2=ids.begin();iter2!=ids.end();iter2++,iter3++)
+ *iter3=_values_added[*iter2];
+ comm.send(&valsToSend[0],ids.size(),MPI_DOUBLE,*iter,1117,*_comm);
+ ids.clear();
+ }
+ _ids_per_working_proc.clear();
+ }
+}
+
--- /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 __GLOBALIZERMESH_HXX__
+#define __GLOBALIZERMESH_HXX__
+
+#include "MEDCouplingNatureOfField.hxx"
+
+#include <mpi.h>
+#include <string>
+#include <vector>
+
+namespace ParaMEDMEM
+{
+ class MEDCouplingFieldDouble;
+
+ class GlobalizerMesh
+ {
+ protected:
+ GlobalizerMesh(const MPI_Comm *comm, MEDCouplingFieldDouble *localField);
+ public:
+ NatureOfField getLocalNature() const;
+ virtual ~GlobalizerMesh();
+ protected:
+ const MPI_Comm *_comm;
+ MEDCouplingFieldDouble *_local_field;
+ };
+
+ class GlobalizerMeshWorkingSide : public GlobalizerMesh
+ {
+ public:
+ GlobalizerMeshWorkingSide(const MPI_Comm *comm, MEDCouplingFieldDouble *localField, const std::string& distantMeth, const std::vector<int>& lazyProcs);
+ ~GlobalizerMeshWorkingSide();
+ //
+ const std::vector<int>& getProcIdsInInteraction() const;
+ void sendSumToLazySide(const std::vector< std::vector<int> >& distantLocEltIds, const std::vector< std::vector<double> >& partialSumRelToDistantIds);
+ void recvSumFromLazySide(std::vector< std::vector<double> >& globalSumRelToDistantIds);
+ private:
+ std::string _distant_method;
+ std::vector<int> _lazy_procs;
+ };
+
+ class GlobalizerMeshLazySide : public GlobalizerMesh
+ {
+ public:
+ GlobalizerMeshLazySide(const MPI_Comm *comm, MEDCouplingFieldDouble *localField, const std::vector<int>& computeProcs);
+ ~GlobalizerMeshLazySide();
+ void recvFromWorkingSide();
+ void sendToWorkingSide();
+ private:
+ std::vector<int> _compute_procs;
+ std::vector<double> _values_added;
+ std::vector< std::vector<int> > _ids_per_working_proc;
+ };
+}
+
+#endif
for(int j=0;j<triofield._nodes_per_elem;j++)
{
conn[j]=(triofield._connectivity)[i*triofield._nodes_per_elem+j];
- }
+ }
_local_mesh->insertNextCell(elemtype,triofield._nodes_per_elem,conn);
}
delete[] conn;
// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
//
#include "ParaMESH.hxx"
+#include "ParaFIELD.hxx"
#include "ProcessorGroup.hxx"
#include "MxN_Mapping.hxx"
#include "InterpolationMatrix.hxx"
#include "MEDCouplingNormalizedUnstructuredMesh.txx"
#include "InterpolationOptions.hxx"
#include "NormalizedUnstructuredMesh.hxx"
+#include "GlobalizerMesh.hxx"
// class InterpolationMatrix
// This class enables the storage of an interpolation matrix Wij mapping
// param method interpolation method
// ====================================================================
- InterpolationMatrix::InterpolationMatrix(ParaMEDMEM::ParaMESH *source_support,
+ InterpolationMatrix::InterpolationMatrix(const ParaMEDMEM::ParaFIELD *source_field,
const ProcessorGroup& source_group,
const ProcessorGroup& target_group,
const DECOptions& dec_options,
const INTERP_KERNEL::InterpolationOptions& interp_options):
- _source_support(source_support->getCellMesh()),
+ _source_field(source_field),
+ _source_support(source_field->getSupport()->getCellMesh()),
_mapping(source_group, target_group, dec_options),
_source_group(source_group),
_target_group(target_group),
- _source_volume(0),
DECOptions(dec_options),
INTERP_KERNEL::InterpolationOptions(interp_options)
{
- int nbelems = _source_support->getNumberOfCells();
-
+ int nbelems = source_field->getField()->getNumberOfTuples();
_row_offsets.resize(nbelems+1);
- for (int i=0; i<nbelems+1; i++)
- {
- _row_offsets[i]=0;
- }
-
_coeffs.resize(nbelems);
+ _target_volume.resize(nbelems);
}
InterpolationMatrix::~InterpolationMatrix()
{
throw INTERP_KERNEL::Exception("local and distant meshes do not have the same space and mesh dimensions");
}
- std::string interpMethod(srcMeth);
- interpMethod+=targetMeth;
+ std::string interpMethod(targetMeth);
+ interpMethod+=srcMeth;
//creating the interpolator structure
vector<map<int,double> > surfaces;
int colSize=0;
source_wrapper.ReleaseTempArrays();
}
else if ( distant_support.getMeshDimension() == 3
- && distant_support.getSpaceDimension() ==3 )
+ && distant_support.getSpaceDimension() == 3 )
{
MEDCouplingNormalizedUnstructuredMesh<3,3> target_wrapper(distant_supportC);
MEDCouplingNormalizedUnstructuredMesh<3,3> source_wrapper(source_supportC);
}
int source_size=surfaces.size();
+ bool needTargetSurf=isSurfaceComputationNeeded(targetMeth);
- MEDCouplingFieldDouble *target_triangle_surf = distant_support.getMeasureField();
- MEDCouplingFieldDouble *source_triangle_surf = _source_support->getMeasureField();
-
- // Storing the source volumes
- _source_volume.resize(source_size);
- for (int i=0; i<source_size; i++)
- {
- _source_volume[i] = source_triangle_surf->getIJ(i,0);
- }
-
- source_triangle_surf->decrRef();
+ MEDCouplingFieldDouble *target_triangle_surf;
+ if(needTargetSurf)
+ target_triangle_surf = distant_support.getMeasureField();
//loop over the elements to build the interpolation
//matrix structures
- for (int ielem=0; ielem < surfaces.size(); ielem++)
+ for (int ielem=0; ielem < source_size; ielem++)
{
_row_offsets[ielem+1] += surfaces[ielem].size();
//_source_indices.push_back(make_pair(iproc_distant, ielem));
iter++)
{
//surface of the target triangle
- double surf = target_triangle_surf->getIJ(iter->first,0);
+ double surf;
+ if(needTargetSurf)
+ surf = target_triangle_surf->getIJ(iter->first,0);
//locating the (iproc, itriangle) pair in the list of columns
map<pair<int,int>,int >::iterator iter2 = _col_offsets.find(make_pair(iproc_distant,iter->first));
//(iproc, itriangle) is not registered in the list
//of distant elements
- col_id =_col_offsets.size()+1;
+ col_id =_col_offsets.size();
_col_offsets.insert(make_pair(make_pair(iproc_distant,iter->first),col_id));
_mapping.addElementFromSource(iproc_distant,
distant_elems[iter->first]);
- _target_volume.push_back(surf);
+ //target_volume.push_back(surf);
}
else
{
//ielem is the row,
//col_id is the number of the column
//iter->second is the value of the coefficient
-
+ if(needTargetSurf)
+ _target_volume[ielem].push_back(surf);
_coeffs[ielem].push_back(make_pair(col_id,iter->second));
}
}
- target_triangle_surf->decrRef();
+ if(needTargetSurf)
+ target_triangle_surf->decrRef();
+ }
+
+ void InterpolationMatrix::finishContributionW(GlobalizerMeshWorkingSide& globalizer)
+ {
+ NatureOfField nature=globalizer.getLocalNature();
+ switch(nature)
+ {
+ case ConservativeVolumic:
+ computeConservVolDenoW(globalizer);
+ break;
+ case Integral:
+ computeIntegralDenoW(globalizer);
+ break;
+ case IntegralGlobConstraint:
+ computeGlobConstraintDenoW(globalizer);
+ break;
+ default:
+ throw INTERP_KERNEL::Exception("Not recognized nature of field. Change nature of Field.");
+ break;
+ }
+ /*for(map<pair<int,int>,int>::iterator iter=_col_offsets.begin();iter!=_col_offsets.end();iter++)
+ if((*iter).second==9)
+ cout << (*iter).first.first << " ** " << (*iter).first.second << endl;
+ cout << "--> " << _col_offsets[make_pair(4,74)] << endl;
+ for(vector<vector<pair<int,double> > >::iterator iter3=_coeffs.begin();iter3!=_coeffs.end();iter3++)
+ for(vector<pair<int,double> >::iterator iter4=(*iter3).begin();iter4!=(*iter3).end();iter4++)
+ if((*iter4).first==569)
+ cout << " __ " << iter3-_coeffs.begin() << "___" << (*iter4).second << endl;
+ ostringstream st; st << "deno_" << _deno_multiply[0][0];
+ ofstream fs(st.str().c_str());
+ for(vector<vector<double> >::iterator iter1=_deno_multiply.begin();iter1!=_deno_multiply.end();iter1++)
+ {
+ for(vector<double>::iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
+ fs << *iter2 << " ";
+ fs << endl;
+ }*/
+ }
+
+ void InterpolationMatrix::finishContributionL(GlobalizerMeshLazySide& globalizer)
+ {
+ NatureOfField nature=globalizer.getLocalNature();
+ switch(nature)
+ {
+ case ConservativeVolumic:
+ computeConservVolDenoL(globalizer);
+ break;
+ case Integral:
+ computeIntegralDenoL(globalizer);
+ break;
+ case IntegralGlobConstraint:
+ computeGlobConstraintDenoL(globalizer);
+ break;
+ default:
+ throw INTERP_KERNEL::Exception("Not recognized nature of field. Change nature of Field.");
+ break;
+ }
+ }
+
+ void InterpolationMatrix::computeConservVolDenoW(GlobalizerMeshWorkingSide& globalizer)
+ {
+ computeGlobalRowSum(globalizer,_deno_multiply);
+ computeGlobalColSum(_deno_reverse_multiply);
}
+ void InterpolationMatrix::computeConservVolDenoL(GlobalizerMeshLazySide& globalizer)
+ {
+ globalizer.recvFromWorkingSide();
+ globalizer.sendToWorkingSide();
+ }
+
+ void InterpolationMatrix::computeIntegralDenoW(GlobalizerMeshWorkingSide& globalizer)
+ {
+ MEDCouplingFieldDouble *source_triangle_surf = _source_support->getMeasureField();
+ _deno_multiply.resize(_coeffs.size());
+ vector<vector<double> >::iterator iter6=_deno_multiply.begin();
+ const double *values=source_triangle_surf->getArray()->getConstPointer();
+ for(vector<vector<pair<int,double> > >::const_iterator iter4=_coeffs.begin();iter4!=_coeffs.end();iter4++,iter6++,values++)
+ {
+ (*iter6).resize((*iter4).size());
+ std::fill((*iter6).begin(),(*iter6).end(),*values);
+ }
+ source_triangle_surf->decrRef();
+ _deno_reverse_multiply=_target_volume;
+ }
+
+ /*!
+ * Nothing to do because surface computation is on working side.
+ */
+ void InterpolationMatrix::computeIntegralDenoL(GlobalizerMeshLazySide& globalizer)
+ {
+ }
+
+ void InterpolationMatrix::computeGlobConstraintDenoW(GlobalizerMeshWorkingSide& globalizer)
+ {
+ computeGlobalColSum(_deno_multiply);
+ computeGlobalRowSum(globalizer,_deno_reverse_multiply);
+ }
+
+ void InterpolationMatrix::computeGlobConstraintDenoL(GlobalizerMeshLazySide& globalizer)
+ {
+ globalizer.recvFromWorkingSide();
+ globalizer.sendToWorkingSide();
+ }
+
+ void InterpolationMatrix::computeGlobalRowSum(GlobalizerMeshWorkingSide& globalizer, std::vector<std::vector<double> >& denoStrorage)
+ {
+ //stores id in distant procs sorted by lazy procs connected with
+ vector< vector<int> > rowsPartialSumI;
+ //stores the corresponding values.
+ vector< vector<double> > rowsPartialSumD;
+ computeLocalRowSum(globalizer.getProcIdsInInteraction(),rowsPartialSumI,rowsPartialSumD);
+ globalizer.sendSumToLazySide(rowsPartialSumI,rowsPartialSumD);
+ globalizer.recvSumFromLazySide(rowsPartialSumD);
+ divideByGlobalRowSum(globalizer.getProcIdsInInteraction(),rowsPartialSumI,rowsPartialSumD,denoStrorage);
+ }
+
+ /*!
+ * @param distantProcs in parameter that indicates which lazy procs are concerned.
+ * @param resPerProcI out parameter that must be cleared before calling this method. The size of 1st dimension is equal to the size of 'distantProcs'.
+ * It contains the element ids (2nd dimension) of the corresponding lazy proc.
+ * @param resPerProcD out parameter with the same format than 'resPerProcI'. It contains corresponding sum values.
+ */
+ void InterpolationMatrix::computeLocalRowSum(const std::vector<int>& distantProcs, std::vector<std::vector<int> >& resPerProcI,
+ std::vector<std::vector<double> >& resPerProcD) const
+ {
+ resPerProcI.resize(distantProcs.size());
+ resPerProcD.resize(distantProcs.size());
+ std::vector<double> res(_col_offsets.size());
+ for(vector<vector<pair<int,double> > >::const_iterator iter=_coeffs.begin();iter!=_coeffs.end();iter++)
+ for(vector<pair<int,double> >::const_iterator iter3=(*iter).begin();iter3!=(*iter).end();iter3++)
+ res[(*iter3).first]+=(*iter3).second;
+ set<int> procsSet;
+ int id;
+ const vector<std::pair<int,int> >& mapping=_mapping.getSendingIds();
+ for(vector<std::pair<int,int> >::const_iterator iter2=mapping.begin();iter2!=mapping.end();iter2++)
+ {
+ std::pair<set<int>::iterator,bool> isIns=procsSet.insert((*iter2).first);
+ if(isIns.second)
+ id=std::find(distantProcs.begin(),distantProcs.end(),(*iter2).first)-distantProcs.begin();
+ resPerProcI[id].push_back((*iter2).second);
+ resPerProcD[id].push_back(res[iter2-mapping.begin()]);
+ }
+ /*
+ for(map<pair<int,int>, int >::const_iterator iter2=_col_offsets.begin();iter2!=_col_offsets.end();iter2++)
+ {
+ std::pair<set<int>::iterator,bool> isIns=procsSet.insert((*iter2).first.first);
+ if(isIns.second)
+ id=std::find(distantProcs.begin(),distantProcs.end(),(*iter2).first.first)-distantProcs.begin();
+ resPerProcI[id].push_back((*iter2).first.second);
+ resPerProcD[id].push_back(res[(*iter2).second]);
+ }*/
+ }
+
+ void InterpolationMatrix::divideByGlobalRowSum(const std::vector<int>& distantProcs, const std::vector<std::vector<int> >& resPerProcI,
+ const std::vector<std::vector<double> >& resPerProcD, std::vector<std::vector<double> >& deno)
+ {
+ map<int,double> fastSums;
+ int procId=0;
+ const vector<std::pair<int,int> >& mapping=_mapping.getSendingIds();
+ map< int, map<int,int> > distIdToLocId;
+ for(map< pair<int,int>,int >::iterator iter8=_col_offsets.begin();iter8!=_col_offsets.end();iter8++)
+ distIdToLocId[(*iter8).first.first][mapping[(*iter8).second].second]=(*iter8).first.second;
+
+ for(vector<int>::const_iterator iter1=distantProcs.begin();iter1!=distantProcs.end();iter1++,procId++)
+ {
+ const std::vector<int>& currentProcI=resPerProcI[procId];
+ const std::vector<double>& currentProcD=resPerProcD[procId];
+ vector<double>::const_iterator iter3=currentProcD.begin();
+ for(vector<int>::const_iterator iter2=currentProcI.begin();iter2!=currentProcI.end();iter2++,iter3++)
+ fastSums[_col_offsets[std::make_pair(*iter1,distIdToLocId[*iter1][*iter2])]]=*iter3;
+ }
+ deno.resize(_coeffs.size());
+ vector<vector<double> >::iterator iter6=deno.begin();
+ for(vector<vector<pair<int,double> > >::const_iterator iter4=_coeffs.begin();iter4!=_coeffs.end();iter4++,iter6++)
+ {
+ (*iter6).resize((*iter4).size());
+ vector<double>::iterator iter7=(*iter6).begin();
+ for(vector<pair<int,double> >::const_iterator iter5=(*iter4).begin();iter5!=(*iter4).end();iter5++,iter7++)
+ *iter7=fastSums[(*iter5).first];
+ }
+ }
+
+ void InterpolationMatrix::computeGlobalColSum(std::vector<std::vector<double> >& denoStrorage)
+ {
+ denoStrorage.resize(_coeffs.size());
+ vector<vector<double> >::iterator iter2=denoStrorage.begin();
+ for(vector<vector<pair<int,double> > >::const_iterator iter1=_coeffs.begin();iter1!=_coeffs.end();iter1++,iter2++)
+ {
+ (*iter2).resize((*iter1).size());
+ double sumOfCurrentRow=0.;
+ for(vector<pair<int,double> >::const_iterator iter3=(*iter1).begin();iter3!=(*iter1).end();iter3++)
+ sumOfCurrentRow+=(*iter3).second;
+ std::fill((*iter2).begin(),(*iter2).end(),sumOfCurrentRow);
+ }
+ }
// ==================================================================
// The call to this method updates the arrays on the target side
void InterpolationMatrix::prepare()
{
- int nbelems = _source_support->getNumberOfCells();
+ int nbelems = _source_field->getField()->getNumberOfTuples();
for (int ielem=0; ielem < nbelems; ielem++)
{
_row_offsets[ielem+1]+=_row_offsets[ielem];
- }
+ }
_mapping.prepareSendRecv();
}
//computing the matrix multiply on source side
if (_source_group.containsMyRank())
{
- int nbrows = _source_support->getNumberOfCells();
+ int nbrows = _coeffs.size();
// performing W.S
// W is the intersection matrix
{
int colid= _coeffs[irow][icol-_row_offsets[irow]].first;
double value = _coeffs[irow][icol-_row_offsets[irow]].second;
- target_value[(colid-1)*nbcomp+icomp]+=value*coeff_row;
+ double deno = _deno_multiply[irow][icol-_row_offsets[irow]];
+ target_value[colid*nbcomp+icomp]+=value*coeff_row/deno;
}
}
}
-
- // performing VT^(-1).(W.S)
- // where VT^(-1) is the inverse of the diagonal matrix containing
- // the volumes of target cells
-
- for (int i=0; i<_col_offsets.size();i++)
- {
- for (int icomp=0; icomp<nbcomp; icomp++)
- {
- target_value[i*nbcomp+icomp] /= _target_volume[i];
- }
- }
-
}
if (_target_group.containsMyRank())
void InterpolationMatrix::transposeMultiply(MEDCouplingFieldDouble& field) const
{
- // int nbcomp = field.getNumberOfComponents();
int nbcomp = field.getArray()->getNumberOfComponents();
vector<double> source_value(_col_offsets.size()* nbcomp,0.0);
_mapping.reverseSendRecv(&source_value[0],field);
//treatment of the transpose matrix multiply on the source side
if (_source_group.containsMyRank())
{
- int nbrows = _source_support->getNumberOfCells();
+ int nbrows = _coeffs.size();
double *array = field.getArray()->getPointer() ;
// Initialization
{
int colid = _coeffs[irow][icol-_row_offsets[irow]].first;
double value = _coeffs[irow][icol-_row_offsets[irow]].second;
-
+ double deno = _deno_reverse_multiply[irow][icol-_row_offsets[irow]];
for (int icomp=0; icomp<nbcomp; icomp++)
{
- double coeff_row = source_value[(colid-1)*nbcomp+icomp];
- array[irow*nbcomp+icomp] += value*coeff_row;
+ double coeff_row = source_value[colid*nbcomp+icomp];
+ array[irow*nbcomp+icomp] += value*coeff_row/deno;
}
}
}
-
- //performing VS^(-1).(WT.T)
- //VS^(-1) is the inverse of the diagonal matrix storing
- //volumes of the source cells
-
- for (int irow=0; irow<nbrows; irow++)
- {
- for (int icomp=0; icomp<nbcomp; icomp++)
- {
- array[irow*nbcomp+icomp] /= _source_volume[irow];
- }
- }
-
}
}
+
+ bool InterpolationMatrix::isSurfaceComputationNeeded(const std::string& method) const
+ {
+ return method=="P0";
+ }
}
namespace ParaMEDMEM
{
+ class GlobalizerMeshWorkingSide;
+ class GlobalizerMeshLazySide;
+
class InterpolationMatrix : public INTERP_KERNEL::InterpolationOptions,
public DECOptions
{
public:
- InterpolationMatrix(ParaMEDMEM::ParaMESH *source_support,
+ InterpolationMatrix(const ParaMEDMEM::ParaFIELD *source_field,
const ProcessorGroup& source_group,
const ProcessorGroup& target_group,
const DECOptions& dec_opt,
virtual ~InterpolationMatrix();
void addContribution(MEDCouplingPointSet& distant_support, int iproc_distant,
int* distant_elems, const std::string& srcMeth, const std::string& targetMeth);
+ void finishContributionW(GlobalizerMeshWorkingSide& globalizer);
+ void finishContributionL(GlobalizerMeshLazySide& globalizer);
void multiply(MEDCouplingFieldDouble& field) const;
void transposeMultiply(MEDCouplingFieldDouble& field)const;
void prepare();
int getNbRows() const { return _row_offsets.size(); }
MPIAccessDEC* getAccessDEC() { return _mapping.getAccessDEC(); }
private:
+ void computeConservVolDenoW(GlobalizerMeshWorkingSide& globalizer);
+ void computeIntegralDenoW(GlobalizerMeshWorkingSide& globalizer);
+ void computeGlobConstraintDenoW(GlobalizerMeshWorkingSide& globalizer);
+ void computeConservVolDenoL(GlobalizerMeshLazySide& globalizer);
+ void computeIntegralDenoL(GlobalizerMeshLazySide& globalizer);
+ void computeGlobConstraintDenoL(GlobalizerMeshLazySide& globalizer);
+
+ void computeLocalColSum(std::vector<double>& res) const;
+ void computeLocalRowSum(const std::vector<int>& distantProcs, std::vector<std::vector<int> >& resPerProcI,
+ std::vector<std::vector<double> >& resPerProcD) const;
+ void computeGlobalRowSum(GlobalizerMeshWorkingSide& globalizer, std::vector<std::vector<double> >& denoStrorage);
+ void computeGlobalColSum(std::vector<std::vector<double> >& denoStrorage);
+ void divideByGlobalRowSum(const std::vector<int>& distantProcs, const std::vector<std::vector<int> >& resPerProcI,
+ const std::vector<std::vector<double> >& resPerProcD, std::vector<std::vector<double> >& deno);
+ private:
+ bool isSurfaceComputationNeeded(const std::string& method) const;
+ private:
+ const ParaMEDMEM::ParaFIELD *_source_field;
std::vector<int> _row_offsets;
std::map<std::pair<int,int>, int > _col_offsets;
- MEDCouplingPointSet *_source_support;
+ MEDCouplingPointSet *_source_support;
MxN_Mapping _mapping;
const ProcessorGroup& _source_group;
const ProcessorGroup& _target_group;
- std::vector<double> _target_volume;
- std::vector<double> _source_volume;
+ std::vector< std::vector<double> > _target_volume;
std::vector<std::vector<std::pair<int,double> > > _coeffs;
+ std::vector<std::vector<double> > _deno_multiply;
+ std::vector<std::vector<double> > _deno_reverse_multiply;
};
}
#include "InterpolationMatrix.hxx"
#include "IntersectionDEC.hxx"
#include "ElementLocator.hxx"
-
-
+#include "GlobalizerMesh.hxx"
namespace ParaMEDMEM
{
*/
void IntersectionDEC::synchronize()
{
- ParaMEDMEM::ParaMESH* para_mesh = _local_field->getSupport();
- //cout <<"size of Interpolation Matrix"<<sizeof(InterpolationMatrix)<<endl;
- _interpolation_matrix = new InterpolationMatrix (para_mesh, *_source_group,*_target_group,*this,*this);
+ delete _interpolation_matrix;
+ _interpolation_matrix = new InterpolationMatrix (_local_field, *_source_group,*_target_group,*this,*this);
//setting up the communication DEC on both sides
if (_source_group->containsMyRank())
{
//locate the distant meshes
- ElementLocator locator(*para_mesh, *_target_group);
+ ElementLocator locator(*_local_field, *_target_group);
//transfering option from IntersectionDEC to ElementLocator
locator.setBoundingBoxAdjustment(getBoundingBoxAdjustment());
MEDCouplingPointSet* distant_mesh=0;
int* distant_ids=0;
+ std::string distantMeth;
for (int i=0; i<_target_group->size(); i++)
{
// int idistant_proc = (i+_source_group->myRank())%_target_group->size();
//gathers pieces of the target meshes that can intersect the local mesh
locator.exchangeMesh(idistant_proc,distant_mesh,distant_ids);
- std::string distantMeth;
- locator.exchangeMethod(_method,idistant_proc,distantMeth);
if (distant_mesh !=0)
{
+ locator.exchangeMethod(_method,idistant_proc,distantMeth);
//adds the contribution of the distant mesh on the local one
int idistant_proc_in_union=_union_group->translateRank(_target_group,idistant_proc);
std::cout <<"add contribution from proc "<<idistant_proc_in_union<<" to proc "<<_union_group->myRank()<<std::endl;
_interpolation_matrix->addContribution(*distant_mesh,idistant_proc_in_union,distant_ids,_method,distantMeth);
-
distant_mesh->decrRef();
delete[] distant_ids;
distant_mesh=0;
distant_ids=0;
}
- }
+ }
+ GlobalizerMeshWorkingSide globalizer(locator.getCommunicator(),_local_field->getField(),distantMeth,locator.getDistantProcIds());
+ _interpolation_matrix->finishContributionW(globalizer);
}
if (_target_group->containsMyRank())
{
- ElementLocator locator(*para_mesh, *_source_group);
+ ElementLocator locator(*_local_field, *_source_group);
//transfering option from IntersectionDEC to ElementLocator
locator.setBoundingBoxAdjustment(getBoundingBoxAdjustment());
//gathers pieces of the target meshes that can intersect the local mesh
locator.exchangeMesh(idistant_proc,distant_mesh,distant_ids);
std::cout << " Data sent from "<<_union_group->myRank()<<" to source proc "<< idistant_proc<<std::endl;
- std::string distantMeth;
- locator.exchangeMethod(_method,idistant_proc,distantMeth);
if (distant_mesh!=0)
{
+ std::string distantMeth;
+ locator.exchangeMethod(_method,idistant_proc,distantMeth);
distant_mesh->decrRef();
delete[] distant_ids;
distant_mesh=0;
distant_ids=0;
}
- }
+ }
+ GlobalizerMeshLazySide globalizer(locator.getCommunicator(),_local_field->getField(),locator.getDistantProcIds());
+ _interpolation_matrix->finishContributionL(globalizer);
}
_interpolation_matrix->prepare();
}
IntersectionDEC.hxx\
ExplicitCoincidentDEC.hxx\
ElementLocator.hxx\
+GlobalizerMesh.hxx\
ExplicitMapping.hxx\
ICoCoField.hxx \
ICoCoMEDField.hxx \
ExplicitCoincidentDEC.cxx\
IntersectionDEC.cxx\
ElementLocator.cxx\
+GlobalizerMesh.cxx\
MPIAccessDEC.cxx \
TimeInterpolator.cxx \
LinearTimeInterpolator.cxx\
void sendRecv(double* field, MEDCouplingFieldDouble& field) const ;
void reverseSendRecv(double* field, MEDCouplingFieldDouble& field) const ;
+ //
+ const std::vector<std::pair<int,int> >& getSendingIds() const { return _sending_ids; }//tmp
+
MPIAccessDEC* getAccessDEC(){ return _access_DEC; }
private :
ProcessorGroup* _union_group;
CPPUNIT_TEST(testIntersectionDEC_2D);
CPPUNIT_TEST(testIntersectionDEC_2DP0P1);
CPPUNIT_TEST(testIntersectionDEC_3D);
+ CPPUNIT_TEST(testIntersectionDECNonOverlapp_2D_P0P0);
+ CPPUNIT_TEST(testIntersectionDECNonOverlapp_2D_P0P1P1P0);
CPPUNIT_TEST(testSynchronousEqualIntersectionWithoutInterpNativeDEC_2D);
CPPUNIT_TEST(testSynchronousEqualIntersectionWithoutInterpDEC_2D);
void testIntersectionDEC_2D();
void testIntersectionDEC_2DP0P1();
void testIntersectionDEC_3D();
+ void testIntersectionDECNonOverlapp_2D_P0P0();
+ void testIntersectionDECNonOverlapp_2D_P0P1P1P0();
#ifdef MED_ENABLE_FVM
void testNonCoincidentDEC_2D();
void testNonCoincidentDEC_3D();
else
dec_emetteur.attachLocalField((ICoCo::Field*) &champ_recepteur);
+ dec_emetteur.setNature(ConservativeVolumic);
if(init)
dec_emetteur.synchronize();
// ParaMEDMEM::ParaSUPPORT* parasupport=new UnstructuredParaSUPPORT( support,*source_group);
ParaMEDMEM::ComponentTopology comptopo;
if(srcM=="P0")
- parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo);
+ {
+ parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo);
+ parafield->getField()->setNature(ConservativeVolumic);
+ }
else
parafield = new ParaFIELD(ON_NODES,NO_TIME,paramesh, comptopo);
int nb_local;
// ParaMEDMEM::ParaSUPPORT* parasupport=new UnstructuredParaSUPPORT(support,*target_group);
ParaMEDMEM::ComponentTopology comptopo;
if(targetM=="P0")
- parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo);
+ {
+ parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo);
+ parafield->getField()->setNature(ConservativeVolumic);
+ }
else
parafield = new ParaFIELD(ON_NODES,NO_TIME,paramesh, comptopo);
int nb_local;
// ParaMEDMEM::ParaSUPPORT* parasupport=new UnstructuredParaSUPPORT( support,*source_group);
ParaMEDMEM::ComponentTopology comptopo;
if(srcM=="P0")
- parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo);
+ {
+ parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo);
+ parafield->getField()->setNature(ConservativeVolumic);
+ }
else
parafield = new ParaFIELD(ON_NODES,NO_TIME,paramesh, comptopo);
int nb_local;
// ParaMEDMEM::ParaSUPPORT* parasupport=new UnstructuredParaSUPPORT(support,*target_group);
ParaMEDMEM::ComponentTopology comptopo;
if(targetM=="P0")
- parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo);
+ {
+ parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo);
+ parafield->getField()->setNature(ConservativeVolumic);
+ }
else
parafield = new ParaFIELD(ON_NODES,NO_TIME,paramesh, comptopo);
int nb_local;
testAsynchronousIntersectionDEC_2D(0.01,1,0.11,1,true,true,true,"P0","P0");
}
+void ParaMEDMEMTest::testIntersectionDECNonOverlapp_2D_P0P0()
+{
+ //
+ const double sourceCoordsAll[2][8]={{0.4,0.5,0.4,1.5,1.6,1.5,1.6,0.5},
+ {0.3,-0.5,1.6,-0.5,1.6,-1.5,0.3,-1.5}};
+ const double targetCoordsAll[3][16]={{0.7,1.45,0.7,1.65,0.9,1.65,0.9,1.45, 1.1,1.4,1.1,1.6,1.3,1.6,1.3,1.4},
+ {0.7,-0.6,0.7,0.7,0.9,0.7,0.9,-0.6, 1.1,-0.7,1.1,0.6,1.3,0.6,1.3,-0.7},
+ {0.7,-1.55,0.7,-1.35,0.9,-1.35,0.9,-1.55, 1.1,-1.65,1.1,-1.45,1.3,-1.45,1.3,-1.65}};
+ int conn4All[8]={0,1,2,3,4,5,6,7};
+ double targetResults[3][2]={{34.,34.},{38.333333333333336,42.666666666666664},{47.,47.}};
+ double targetResults2[3][2]={{0.28333333333333344,0.56666666666666687},{1.8564102564102569,2.0128205128205132},{1.0846153846153845,0.36153846153846159}};
+ double targetResults3[3][2]={{3.7777777777777781,7.5555555555555562},{24.511111111111113,26.355555555555558},{14.1,4.7}};
+ //
+ int size;
+ int rank;
+ MPI_Comm_size(MPI_COMM_WORLD,&size);
+ MPI_Comm_rank(MPI_COMM_WORLD,&rank);
+ //
+ if(size!=5)
+ return ;
+ int nproc_source = 2;
+ set<int> self_procs;
+ set<int> procs_source;
+ set<int> procs_target;
+
+ for (int i=0; i<nproc_source; i++)
+ procs_source.insert(i);
+ for (int i=nproc_source; i<size; i++)
+ procs_target.insert(i);
+ self_procs.insert(rank);
+ //
+ ParaMEDMEM::MEDCouplingUMesh *mesh=0;
+ ParaMEDMEM::ParaMESH *paramesh=0;
+ ParaMEDMEM::ParaFIELD* parafield=0;
+ ICoCo::Field* icocofield=0;
+ //
+ ParaMEDMEM::CommInterface interface;
+ //
+ ProcessorGroup* self_group = new ParaMEDMEM::MPIProcessorGroup(interface,self_procs);
+ ProcessorGroup* target_group = new ParaMEDMEM::MPIProcessorGroup(interface,procs_target);
+ ProcessorGroup* source_group = new ParaMEDMEM::MPIProcessorGroup(interface,procs_source);
+ //
+ MPI_Barrier(MPI_COMM_WORLD);
+ if(source_group->containsMyRank())
+ {
+ std::ostringstream stream; stream << "sourcemesh2D proc " << rank;
+ mesh=MEDCouplingUMesh::New(stream.str().c_str(),2);
+ mesh->allocateCells(2);
+ mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn4All);
+ mesh->finishInsertingCells();
+ DataArrayDouble *myCoords=DataArrayDouble::New();
+ myCoords->alloc(4,2);
+ const double *sourceCoords=sourceCoordsAll[rank];
+ std::copy(sourceCoords,sourceCoords+8,myCoords->getPointer());
+ mesh->setCoords(myCoords);
+ myCoords->decrRef();
+ paramesh=new ParaMESH(mesh,*source_group,"source mesh");
+ ParaMEDMEM::ComponentTopology comptopo;
+ parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo);
+ double *value=parafield->getField()->getArray()->getPointer();
+ value[0]=34+13*((double)rank);
+ }
+ else
+ {
+ std::ostringstream stream; stream << "targetmesh2D proc " << rank-nproc_source;
+ mesh=MEDCouplingUMesh::New(stream.str().c_str(),2);
+ mesh->allocateCells(2);
+ mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn4All);
+ mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn4All+4);
+ mesh->finishInsertingCells();
+ DataArrayDouble *myCoords=DataArrayDouble::New();
+ myCoords->alloc(8,2);
+ const double *targetCoords=targetCoordsAll[rank-nproc_source];
+ std::copy(targetCoords,targetCoords+16,myCoords->getPointer());
+ mesh->setCoords(myCoords);
+ myCoords->decrRef();
+ paramesh=new ParaMESH (mesh,*target_group,"target mesh");
+ ParaMEDMEM::ComponentTopology comptopo;
+ parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo);
+ }
+ //test 1 - Conservative volumic
+ ParaMEDMEM::IntersectionDEC dec(*source_group,*target_group);
+ parafield->getField()->setNature(ConservativeVolumic);
+ if (source_group->containsMyRank())
+ {
+ dec.setMethod("P0");
+ dec.attachLocalField(parafield);
+ dec.synchronize();
+ dec.setForcedRenormalization(false);
+ dec.sendData();
+ }
+ else
+ {
+ dec.setMethod("P0");
+ dec.attachLocalField(parafield);
+ dec.synchronize();
+ dec.setForcedRenormalization(false);
+ dec.recvData();
+ const double *res=parafield->getField()->getArray()->getConstPointer();
+ const double *expected=targetResults[rank-nproc_source];
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(expected[0],res[0],1e-13);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(expected[1],res[1],1e-13);
+ }
+ //test 2 - Integral
+ ParaMEDMEM::IntersectionDEC dec2(*source_group,*target_group);
+ parafield->getField()->setNature(Integral);
+ if (source_group->containsMyRank())
+ {
+ dec2.setMethod("P0");
+ dec2.attachLocalField(parafield);
+ dec2.synchronize();
+ dec2.setForcedRenormalization(false);
+ dec2.sendData();
+ }
+ else
+ {
+ dec2.setMethod("P0");
+ dec2.attachLocalField(parafield);
+ dec2.synchronize();
+ dec2.setForcedRenormalization(false);
+ dec2.recvData();
+ const double *res=parafield->getField()->getArray()->getConstPointer();
+ const double *expected=targetResults2[rank-nproc_source];
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(expected[0],res[0],1e-13);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(expected[1],res[1],1e-13);
+ }
+ //test 3 - Integral with global constraint
+ ParaMEDMEM::IntersectionDEC dec3(*source_group,*target_group);
+ parafield->getField()->setNature(IntegralGlobConstraint);
+ if (source_group->containsMyRank())
+ {
+ dec3.setMethod("P0");
+ dec3.attachLocalField(parafield);
+ dec3.synchronize();
+ dec3.setForcedRenormalization(false);
+ dec3.sendData();
+ }
+ else
+ {
+ dec3.setMethod("P0");
+ dec3.attachLocalField(parafield);
+ dec3.synchronize();
+ dec3.setForcedRenormalization(false);
+ dec3.recvData();
+ const double *res=parafield->getField()->getArray()->getConstPointer();
+ const double *expected=targetResults3[rank-nproc_source];
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(expected[0],res[0],1e-13);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(expected[1],res[1],1e-13);
+ }
+ //test 4 - Conservative volumic reversed
+ ParaMEDMEM::IntersectionDEC dec4(*source_group,*target_group);
+ parafield->getField()->setNature(ConservativeVolumic);
+ if (source_group->containsMyRank())
+ {
+ dec4.setMethod("P0");
+ dec4.attachLocalField(parafield);
+ dec4.synchronize();
+ dec4.setForcedRenormalization(false);
+ dec4.recvData();
+ const double *res=parafield->getField()->getArray()->getConstPointer();
+ CPPUNIT_ASSERT_EQUAL(1,parafield->getField()->getNumberOfTuples());
+ const double expected[]={37.8518518518519,43.5333333333333};
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(expected[rank],res[0],1e-13);
+ }
+ else
+ {
+ dec4.setMethod("P0");
+ dec4.attachLocalField(parafield);
+ dec4.synchronize();
+ dec4.setForcedRenormalization(false);
+ double *res=parafield->getField()->getArray()->getPointer();
+ const double *toSet=targetResults[rank-nproc_source];
+ res[0]=toSet[0];
+ res[1]=toSet[1];
+ dec4.sendData();
+ }
+ //test 5 - Integral reversed
+ ParaMEDMEM::IntersectionDEC dec5(*source_group,*target_group);
+ parafield->getField()->setNature(Integral);
+ if (source_group->containsMyRank())
+ {
+ dec5.setMethod("P0");
+ dec5.attachLocalField(parafield);
+ dec5.synchronize();
+ dec5.setForcedRenormalization(false);
+ dec5.recvData();
+ const double *res=parafield->getField()->getArray()->getConstPointer();
+ CPPUNIT_ASSERT_EQUAL(1,parafield->getField()->getNumberOfTuples());
+ const double expected[]={0.794600591715977,1.35631163708087};
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(expected[rank],res[0],1e-13);
+ }
+ else
+ {
+ dec5.setMethod("P0");
+ dec5.attachLocalField(parafield);
+ dec5.synchronize();
+ dec5.setForcedRenormalization(false);
+ double *res=parafield->getField()->getArray()->getPointer();
+ const double *toSet=targetResults2[rank-nproc_source];
+ res[0]=toSet[0];
+ res[1]=toSet[1];
+ dec5.sendData();
+ }
+ //test 6 - Integral with global constraint reversed
+ ParaMEDMEM::IntersectionDEC dec6(*source_group,*target_group);
+ parafield->getField()->setNature(IntegralGlobConstraint);
+ if (source_group->containsMyRank())
+ {
+ dec6.setMethod("P0");
+ dec6.attachLocalField(parafield);
+ dec6.synchronize();
+ dec6.setForcedRenormalization(false);
+ dec6.recvData();
+ const double *res=parafield->getField()->getArray()->getConstPointer();
+ CPPUNIT_ASSERT_EQUAL(1,parafield->getField()->getNumberOfTuples());
+ const double expected[]={36.4592592592593,44.5407407407407};
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(expected[rank],res[0],1e-13);
+ }
+ else
+ {
+ dec6.setMethod("P0");
+ dec6.attachLocalField(parafield);
+ dec6.synchronize();
+ dec6.setForcedRenormalization(false);
+ double *res=parafield->getField()->getArray()->getPointer();
+ const double *toSet=targetResults3[rank-nproc_source];
+ res[0]=toSet[0];
+ res[1]=toSet[1];
+ dec6.sendData();
+ }
+ //
+ delete parafield;
+ mesh->decrRef();
+ delete paramesh;
+ delete self_group;
+ delete target_group;
+ delete source_group;
+ //
+ MPI_Barrier(MPI_COMM_WORLD);
+}
+
+void ParaMEDMEMTest::testIntersectionDECNonOverlapp_2D_P0P1P1P0()
+{
+ int size;
+ int rank;
+ MPI_Comm_size(MPI_COMM_WORLD,&size);
+ MPI_Comm_rank(MPI_COMM_WORLD,&rank);
+ //
+ if(size!=5)
+ return ;
+ int nproc_source = 2;
+ set<int> self_procs;
+ set<int> procs_source;
+ set<int> procs_target;
+
+ for (int i=0; i<nproc_source; i++)
+ procs_source.insert(i);
+ for (int i=nproc_source; i<size; i++)
+ procs_target.insert(i);
+ self_procs.insert(rank);
+ //
+ ParaMEDMEM::MEDCouplingUMesh *mesh=0;
+ ParaMEDMEM::ParaMESH *paramesh=0;
+ ParaMEDMEM::ParaFIELD *parafieldP0=0,*parafieldP1=0;
+ ICoCo::Field* icocofield=0;
+ //
+ ParaMEDMEM::CommInterface interface;
+ //
+ ProcessorGroup* self_group = new ParaMEDMEM::MPIProcessorGroup(interface,self_procs);
+ ProcessorGroup* target_group = new ParaMEDMEM::MPIProcessorGroup(interface,procs_target);
+ ProcessorGroup* source_group = new ParaMEDMEM::MPIProcessorGroup(interface,procs_source);
+ //
+ MPI_Barrier(MPI_COMM_WORLD);
+ if(source_group->containsMyRank())
+ {
+ if(rank==0)
+ {
+ double coords[6]={-0.3,-0.3, 0.7,0.7, 0.7,-0.3};
+ int conn[3]={0,1,2};
+ int globalNode[3]={1,2,0};
+ mesh=MEDCouplingUMesh::New("Source mesh Proc0",2);
+ mesh->allocateCells(1);
+ mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,conn);
+ mesh->finishInsertingCells();
+ DataArrayDouble *myCoords=DataArrayDouble::New();
+ myCoords->alloc(3,2);
+ std::copy(coords,coords+6,myCoords->getPointer());
+ mesh->setCoords(myCoords);
+ myCoords->decrRef();
+ }
+ if(rank==1)
+ {
+ double coords[6]={-0.3,-0.3, -0.3,0.7, 0.7,0.7};
+ int conn[3]={0,1,2};
+ int globalNode[3]={1,3,2};
+ mesh=MEDCouplingUMesh::New("Source mesh Proc1",2);
+ mesh->allocateCells(1);
+ mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,conn);
+ mesh->finishInsertingCells();
+ DataArrayDouble *myCoords=DataArrayDouble::New();
+ myCoords->alloc(3,2);
+ std::copy(coords,coords+6,myCoords->getPointer());
+ mesh->setCoords(myCoords);
+ myCoords->decrRef();
+ }
+ paramesh=new ParaMESH(mesh,*source_group,"source mesh");
+ ParaMEDMEM::ComponentTopology comptopo;
+ parafieldP0 = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo);
+ parafieldP1 = new ParaFIELD(ON_NODES,NO_TIME,paramesh, comptopo);
+ double *valueP0=parafieldP0->getField()->getArray()->getPointer();
+ double *valueP1=parafieldP1->getField()->getArray()->getPointer();
+ parafieldP0->getField()->setNature(ConservativeVolumic);
+ parafieldP1->getField()->setNature(ConservativeVolumic);
+ if(rank==0)
+ {
+ valueP0[0]=31.;
+ valueP1[0]=34.; valueP1[1]=77.; valueP1[2]=53.;
+ }
+ if(rank==1)
+ {
+ valueP0[0]=47.;
+ valueP1[0]=34.; valueP1[1]=57.; valueP1[2]=77.;
+ }
+ }
+ else
+ {
+ if(rank==2)
+ {
+ double coords[10]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2 };
+ int conn[7]={0,3,4,1, 1,4,2};
+ int globalNode[5]={4,3,0,2,1};
+ mesh=MEDCouplingUMesh::New("Target mesh Proc2",2);
+ mesh->allocateCells(2);
+ mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn);
+ mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,conn+4);
+ mesh->finishInsertingCells();
+ DataArrayDouble *myCoords=DataArrayDouble::New();
+ myCoords->alloc(5,2);
+ std::copy(coords,coords+10,myCoords->getPointer());
+ mesh->setCoords(myCoords);
+ myCoords->decrRef();
+ }
+ if(rank==3)
+ {
+ double coords[6]={0.2,0.2, 0.7,-0.3, 0.7,0.2};
+ int conn[3]={0,2,1};
+ int globalNode[3]={1,0,5};
+ mesh=MEDCouplingUMesh::New("Target mesh Proc3",2);
+ mesh->allocateCells(1);
+ mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,conn);
+ mesh->finishInsertingCells();
+ DataArrayDouble *myCoords=DataArrayDouble::New();
+ myCoords->alloc(3,2);
+ std::copy(coords,coords+6,myCoords->getPointer());
+ mesh->setCoords(myCoords);
+ myCoords->decrRef();
+ }
+ if(rank==4)
+ {
+ double coords[12]={-0.3,0.2, -0.3,0.7, 0.2,0.7, 0.2,0.2, 0.7,0.7, 0.7,0.2};
+ int conn[8]={0,1,2,3, 3,2,4,5};
+ int globalNode[6]={2,6,7,1,8,5};
+ mesh=MEDCouplingUMesh::New("Target mesh Proc4",2);
+ mesh->allocateCells(2);
+ mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn);
+ mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+4);
+ mesh->finishInsertingCells();
+ DataArrayDouble *myCoords=DataArrayDouble::New();
+ myCoords->alloc(6,2);
+ std::copy(coords,coords+12,myCoords->getPointer());
+ mesh->setCoords(myCoords);
+ myCoords->decrRef();
+ }
+ ParaMEDMEM::ComponentTopology comptopo;
+ paramesh=new ParaMESH(mesh,*target_group,"target mesh");
+ parafieldP0 = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo);
+ parafieldP1 = new ParaFIELD(ON_NODES,NO_TIME,paramesh, comptopo);
+ parafieldP0->getField()->setNature(ConservativeVolumic);
+ parafieldP1->getField()->setNature(ConservativeVolumic);
+ }
+ // test 1 - P0 P1
+ ParaMEDMEM::IntersectionDEC dec(*source_group,*target_group);
+ if (source_group->containsMyRank())
+ {
+ dec.setMethod("P0");
+ dec.attachLocalField(parafieldP0);
+ dec.synchronize();
+ dec.setForcedRenormalization(false);
+ dec.sendData();
+ }
+ else
+ {
+ dec.setMethod("P1");
+ dec.attachLocalField(parafieldP1);
+ dec.synchronize();
+ dec.setForcedRenormalization(false);
+ dec.recvData();
+ /*const double *res=parafield->getField()->getArray()->getConstPointer();
+ const double *expected=targetResults[rank-nproc_source];
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(expected[0],res[0],1e-13);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(expected[1],res[1],1e-13);*/
+ }
+ //
+ delete parafieldP0;
+ delete parafieldP1;
+ mesh->decrRef();
+ delete paramesh;
+ delete self_group;
+ delete target_group;
+ delete source_group;
+ //
+ MPI_Barrier(MPI_COMM_WORLD);
+}
+
/*!
* Tests an asynchronous exchange between two codes
* one sends data with dtA as an interval, the max time being tmaxA
meshname<< "Mesh_2_"<< rank+1;
mesh=MEDLoader::ReadUMeshFromFile(strstream.str().c_str(),meshname.str().c_str(),0);
-
+
paramesh=new ParaMESH (mesh,*source_group,"source mesh");
// ParaMEDMEM::ParaSUPPORT* parasupport=new UnstructuredParaSUPPORT( support,*source_group);
ParaMEDMEM::ComponentTopology comptopo;
if(srcM=="P0")
- parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo);
+ {
+ parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo);
+ parafield->getField()->setNature(ConservativeVolumic);//InvertIntegral);//ConservativeVolumic);
+ }
else
parafield = new ParaFIELD(ON_NODES,NO_TIME,paramesh, comptopo);
meshname<< "Mesh_3_"<<rank-nproc_source+1;
mesh = MEDLoader::ReadUMeshFromFile(strstream.str().c_str(),meshname.str().c_str(),0);
-
+
paramesh=new ParaMESH (mesh,*target_group,"target mesh");
// ParaMEDMEM::ParaSUPPORT* parasupport=new UnstructuredParaSUPPORT(support,*target_group);
ParaMEDMEM::ComponentTopology comptopo;
if(targetM=="P0")
- parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo);
+ {
+ parafield = new ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo);
+ parafield->getField()->setNature(ConservativeVolumic);//InvertIntegral);//ConservativeVolumic);
+ }
else
parafield = new ParaFIELD(ON_NODES,NO_TIME,paramesh, comptopo);
paramesh=ParaMESH(mesh,source_group,"source mesh")
comptopo = ComponentTopology()
parafield = ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo)
+ parafield.getField().setNature(ConservativeVolumic)
nb_local=mesh.getNumberOfCells()
value = [1.0]*nb_local
parafield.getField().setValues(value)
paramesh=ParaMESH(mesh,target_group,"target mesh")
comptopo = ComponentTopology()
parafield = ParaFIELD(ON_CELLS,NO_TIME,paramesh, comptopo)
+ parafield.getField().setNature(ConservativeVolumic)
nb_local=mesh.getNumberOfCells()
value = [0.0]*nb_local
parafield.getField().setValues(value)
