From 406e566afcfdb48e7c712fecce0707e18029d28c Mon Sep 17 00:00:00 2001 From: eap Date: Tue, 29 Sep 2009 06:34:30 +0000 Subject: [PATCH] 0020440: [CEA 349] P1P0 barycentric interpolators "Concerning engine intersectors TriangleIntersector,Geometric2DIntersector and ConvexIntersector have to be extended with and additional method:" + double intersectGeoBary(const std::vector& targetCell, bool targetCellQuadratic, const double *sourceCell, std::vector& res); --- src/INTERP_KERNEL/ConvexIntersector.hxx | 2 + src/INTERP_KERNEL/ConvexIntersector.txx | 79 +++++++++++-- src/INTERP_KERNEL/Geometric2DIntersector.hxx | 2 + src/INTERP_KERNEL/Geometric2DIntersector.txx | 93 +++++++++++++--- .../TriangulationIntersector.hxx | 3 + .../TriangulationIntersector.txx | 104 +++++++++++++++--- 6 files changed, 241 insertions(+), 42 deletions(-) diff --git a/src/INTERP_KERNEL/ConvexIntersector.hxx b/src/INTERP_KERNEL/ConvexIntersector.hxx index 094b4c76f..91e6d054b 100644 --- a/src/INTERP_KERNEL/ConvexIntersector.hxx +++ b/src/INTERP_KERNEL/ConvexIntersector.hxx @@ -23,6 +23,7 @@ #include "PlanarIntersectorP0P1.hxx" #include "PlanarIntersectorP1P0.hxx" #include "PlanarIntersectorP1P1.hxx" +#include "PlanarIntersectorP1P0Bary.hxx" namespace INTERP_KERNEL { @@ -41,6 +42,7 @@ namespace INTERP_KERNEL double intersectGeometry(ConnType icellT, ConnType icellS, ConnType nbNodesT, ConnType nbNodesS); double intersectGeometryWithQuadrangle(const double *quadrangle, const std::vector& sourceCoords, bool isSourceQuad); double intersectGeometryGeneral(const std::vector& targetCoords, const std::vector& sourceCoords); + double intersectGeoBary(const std::vector& targetCell, bool targetCellQuadratic, const double *sourceCell, std::vector& res); private : double _epsilon; }; diff --git a/src/INTERP_KERNEL/ConvexIntersector.txx b/src/INTERP_KERNEL/ConvexIntersector.txx index 70576f39f..317fd2393 100644 --- a/src/INTERP_KERNEL/ConvexIntersector.txx +++ b/src/INTERP_KERNEL/ConvexIntersector.txx @@ -24,18 +24,22 @@ #include "PlanarIntersectorP0P1.txx" #include "PlanarIntersectorP1P0.txx" #include "PlanarIntersectorP1P1.txx" +#include "PlanarIntersectorP1P0Bary.txx" #include "PolygonAlgorithms.txx" #include +#define CONVINTERSECTOR_TEMPLATE template class InterpType> +#define CONVEX_INTERSECTOR_ ConvexIntersector + namespace INTERP_KERNEL { - template class InterpType> - ConvexIntersector::ConvexIntersector(const MyMeshType& meshT, const MyMeshType& meshS, - double dimCaracteristic, double precision, double md3DSurf, - double medianPlane, bool doRotate , int oriantation, int printLevel) - :InterpType >(meshT,meshS,dimCaracteristic, precision, md3DSurf, medianPlane, doRotate, oriantation, printLevel), + CONVINTERSECTOR_TEMPLATE + CONVEX_INTERSECTOR_::ConvexIntersector(const MyMeshType& meshT, const MyMeshType& meshS, + double dimCaracteristic, double precision, double md3DSurf, + double medianPlane, bool doRotate , int oriantation, int printLevel) + :InterpType(meshT,meshS,dimCaracteristic, precision, md3DSurf, medianPlane, doRotate, oriantation, printLevel), _epsilon(precision*dimCaracteristic) { if(PlanarIntersector::_print_level >= 1) @@ -48,8 +52,9 @@ namespace INTERP_KERNEL } } - template class InterpType> - double ConvexIntersector::intersectGeometry(ConnType icellT, ConnType icellS, ConnType nbNodesT, ConnType nbNodesS) + CONVINTERSECTOR_TEMPLATE + double CONVEX_INTERSECTOR_::intersectGeometry(ConnType icellT, ConnType icellS, + ConnType nbNodesT, ConnType nbNodesS) { double result = 0; int orientation = 1; @@ -82,8 +87,10 @@ namespace INTERP_KERNEL return orientation*result; } - template class InterpType> - double ConvexIntersector::intersectGeometryWithQuadrangle(const double *quadrangle, const std::vector& sourceCoords, bool isSourceQuad) + CONVINTERSECTOR_TEMPLATE + double CONVEX_INTERSECTOR_::intersectGeometryWithQuadrangle(const double * quadrangle, + const std::vector& sourceCoords, + bool isSourceQuad) { double result = 0; int nbOfNodesS=sourceCoords.size()/SPACEDIM; @@ -112,8 +119,9 @@ namespace INTERP_KERNEL return result; } - template class InterpType> - double ConvexIntersector::intersectGeometryGeneral(const std::vector& targetCoords, const std::vector& sourceCoords) + CONVINTERSECTOR_TEMPLATE + double CONVEX_INTERSECTOR_::intersectGeometryGeneral(const std::vector& targetCoords, + const std::vector& sourceCoords) { double result = 0; int nbOfNodesS=sourceCoords.size()/SPACEDIM; @@ -131,6 +139,55 @@ namespace INTERP_KERNEL } return result; } + + //================================================================================ + /*! + * \brief Intersect a triangle and a polygon for P1P0 barycentric algorithm + * \param targetCell - list of coordinates of target polygon in full interlace + * \param targetCellQuadratic - specifies if target polygon is quadratic or not + * \param sourceTria - list of coordinates of source triangle + * \param res - coefficients a,b and c associated to nodes of sourceTria + */ + //================================================================================ + + CONVINTERSECTOR_TEMPLATE + double CONVEX_INTERSECTOR_::intersectGeoBary(const std::vector& targetCell, + bool targetCellQuadratic, + const double * sourceTria, + std::vector& res) + { + double area = 0; + double barycenter[SPACEDIM] = {0., 0.}; + int nbOfNodesT=targetCell.size()/SPACEDIM; + + /*** Compute the intersection area ***/ + INTERP_KERNEL::PolygonAlgorithms P(_epsilon, PlanarIntersector::_precision); + std::deque inter = P.intersectConvexPolygons(sourceTria, &targetCell[0], 3, nbOfNodesT); + double cross[SPACEDIM]; + int nb_inter =((int)inter.size())/SPACEDIM; + for(int i = 1; i(&inter[0],&inter[SPACEDIM*i],&inter[SPACEDIM*(i+1)],cross); + area += 0.5*norm(cross); + barycenter[0] += inter[SPACEDIM*i]; + barycenter[1] += inter[SPACEDIM*i+1]; + } + if ( area > std::numeric_limits::min() ) + { + barycenter[0] = ( barycenter[0] + inter[0] + inter[SPACEDIM*(nb_inter-1)] ) / nb_inter; + barycenter[1] = ( barycenter[1] + inter[1] + inter[SPACEDIM*(nb_inter-1)+1]) / nb_inter; + res.resize(3); + barycentric_coords( sourceTria, &barycenter[0], &res[0]); + res[0] *= area; + res[1] *= area; + res[2] *= area; + } + else + { + area = 0; + } + return area; + } } #endif diff --git a/src/INTERP_KERNEL/Geometric2DIntersector.hxx b/src/INTERP_KERNEL/Geometric2DIntersector.hxx index 257d6516b..303b24a80 100644 --- a/src/INTERP_KERNEL/Geometric2DIntersector.hxx +++ b/src/INTERP_KERNEL/Geometric2DIntersector.hxx @@ -23,6 +23,7 @@ #include "PlanarIntersectorP0P1.hxx" #include "PlanarIntersectorP1P0.hxx" #include "PlanarIntersectorP1P1.hxx" +#include "PlanarIntersectorP1P0Bary.hxx" namespace INTERP_KERNEL { @@ -42,6 +43,7 @@ namespace INTERP_KERNEL double intersectGeometry(ConnType icellT, ConnType icellS, ConnType nbNodesT, ConnType nbNodesS); double intersectGeometryWithQuadrangle(const double *quadrangle, const std::vector& sourceCoords, bool isSourceQuad); double intersectGeometryGeneral(const std::vector& targetCoords, const std::vector& sourceCoords); + double intersectGeoBary(const std::vector& targetCell, bool targetCellQuadratic, const double *sourceCell, std::vector& res); private: QuadraticPolygon *buildPolygonFrom(const std::vector& coords, NormalizedCellType type); QuadraticPolygon *buildPolygonAFrom(ConnType cell, int nbOfPoints, NormalizedCellType type); diff --git a/src/INTERP_KERNEL/Geometric2DIntersector.txx b/src/INTERP_KERNEL/Geometric2DIntersector.txx index 3592c35b2..44991812e 100644 --- a/src/INTERP_KERNEL/Geometric2DIntersector.txx +++ b/src/INTERP_KERNEL/Geometric2DIntersector.txx @@ -24,6 +24,7 @@ #include "PlanarIntersectorP0P1.txx" #include "PlanarIntersectorP1P0.txx" #include "PlanarIntersectorP1P1.txx" +#include "PlanarIntersectorP1P0Bary.txx" #include "CellModel.hxx" #include "QuadraticPolygon.hxx" @@ -31,18 +32,23 @@ #include "EdgeLin.hxx" #include "Node.hxx" +#define GEO2D_INTERSECTOR Geometric2DIntersector +#define INTERSECTOR_TEMPLATE template class InterpType> + namespace INTERP_KERNEL { - template class InterpType> - Geometric2DIntersector::Geometric2DIntersector(const MyMeshType& meshT, const MyMeshType& meshS, - double dimCaracteristic, double md3DSurf, double medianPlane, double precision, int orientation): - InterpType >(meshT,meshS,dimCaracteristic, precision, md3DSurf, medianPlane, true, orientation, 0) + INTERSECTOR_TEMPLATE + GEO2D_INTERSECTOR::Geometric2DIntersector(const MyMeshType& meshT, const MyMeshType& meshS, + double dimCaracteristic, double md3DSurf, double medianPlane, + double precision, int orientation): + InterpType(meshT,meshS,dimCaracteristic, precision, md3DSurf, medianPlane, true, orientation, 0) { QUADRATIC_PLANAR::_precision=dimCaracteristic*precision; } - template class InterpType> - double Geometric2DIntersector::intersectGeometry(ConnType icellT, ConnType icellS, ConnType nbNodesT, ConnType nbNodesS) + INTERSECTOR_TEMPLATE + double GEO2D_INTERSECTOR::intersectGeometry(ConnType icellT, ConnType icellS, + ConnType nbNodesT, ConnType nbNodesS) { int orientation = 1; std::vector CoordsT; @@ -57,8 +63,10 @@ namespace INTERP_KERNEL return ret; } - template class InterpType> - double Geometric2DIntersector::intersectGeometryWithQuadrangle(const double *quadrangle, const std::vector& sourceCoords, bool isSourceQuad) + INTERSECTOR_TEMPLATE + double GEO2D_INTERSECTOR::intersectGeometryWithQuadrangle(const double * quadrangle, + const std::vector& sourceCoords, + bool isSourceQuad) { std::vector nodes(4); nodes[0]=new Node(quadrangle[0],quadrangle[1]); @@ -80,8 +88,9 @@ namespace INTERP_KERNEL return ret; } - template class InterpType> - double Geometric2DIntersector::intersectGeometryGeneral(const std::vector& targetCoords, const std::vector& sourceCoords) + INTERSECTOR_TEMPLATE + double GEO2D_INTERSECTOR::intersectGeometryGeneral(const std::vector& targetCoords, + const std::vector& sourceCoords) { int nbOfTargetNodes=targetCoords.size()/SPACEDIM; std::vector nodes(nbOfTargetNodes); @@ -98,8 +107,60 @@ namespace INTERP_KERNEL return ret; } - template class InterpType> - QuadraticPolygon *Geometric2DIntersector::buildPolygonFrom(const std::vector& coords, NormalizedCellType type) + //================================================================================ + /*! + * \brief Intersect a triangle and a polygon for P1P0 barycentric algorithm + * \param targetCell - list of coordinates of target polygon in full interlace + * \param targetCellQuadratic - specifies if target polygon is quadratic or not + * \param sourceTria - list of coordinates of source triangle + * \param res - coefficients a,b and c associated to nodes of sourceTria + */ + //================================================================================ + + INTERSECTOR_TEMPLATE + double GEO2D_INTERSECTOR::intersectGeoBary(const std::vector& targetCell, + bool targetCellQuadratic, + const double * sourceTria, + std::vector& res) + { + std::vector nodes(3); + nodes[0]=new Node(sourceTria[0*SPACEDIM],sourceTria[0*SPACEDIM+1]); + nodes[1]=new Node(sourceTria[1*SPACEDIM],sourceTria[1*SPACEDIM+1]); + nodes[2]=new Node(sourceTria[2*SPACEDIM],sourceTria[2*SPACEDIM+1]); + int nbOfTargetNodes=targetCell.size()/SPACEDIM; + std::vector nodes2(nbOfTargetNodes); + for(int i=0;iintersectWith(*p2,barycenter); + delete p1; delete p2; + if ( ret > std::numeric_limits::min() ) + { + std::vector sourceCell(3); + sourceCell[0] = &sourceTria[0]; + sourceCell[1] = &sourceTria[SPACEDIM]; + sourceCell[2] = &sourceTria[SPACEDIM*2]; + res.resize(3); + barycentric_coords( sourceCell, barycenter, &res[0]); + res[0] *= ret; + res[1] *= ret; + res[2] *= ret; + } + else + { + ret = 0; + } + return ret; + } + + INTERSECTOR_TEMPLATE + QuadraticPolygon *GEO2D_INTERSECTOR::buildPolygonFrom(const std::vector& coords, NormalizedCellType type) { int nbNodes=coords.size()/SPACEDIM; std::vector nodes(nbNodes); @@ -111,8 +172,8 @@ namespace INTERP_KERNEL return QuadraticPolygon::buildArcCirclePolygon(nodes); } - template class InterpType> - QuadraticPolygon *Geometric2DIntersector::buildPolygonAFrom(ConnType cell, int nbOfPoints, NormalizedCellType type) + INTERSECTOR_TEMPLATE + QuadraticPolygon *GEO2D_INTERSECTOR::buildPolygonAFrom(ConnType cell, int nbOfPoints, NormalizedCellType type) { const ConnType *startOfCellNodeConn=PlanarIntersector::_connectT+OTT::conn2C(PlanarIntersector::_connIndexT[OTT::ind2C(cell)]); std::vector nodes(nbOfPoints); @@ -124,8 +185,8 @@ namespace INTERP_KERNEL return QuadraticPolygon::buildArcCirclePolygon(nodes); } - template class InterpType> - QuadraticPolygon *Geometric2DIntersector::buildPolygonBFrom(ConnType cell, int nbOfPoints, NormalizedCellType type) + INTERSECTOR_TEMPLATE + QuadraticPolygon *GEO2D_INTERSECTOR::buildPolygonBFrom(ConnType cell, int nbOfPoints, NormalizedCellType type) { const ConnType *startOfCellNodeConn=PlanarIntersector::_connectS+OTT::conn2C(PlanarIntersector::_connIndexS[OTT::ind2C(cell)]); std::vector nodes(nbOfPoints); diff --git a/src/INTERP_KERNEL/TriangulationIntersector.hxx b/src/INTERP_KERNEL/TriangulationIntersector.hxx index 2f93f9c40..f1c3dfd6f 100644 --- a/src/INTERP_KERNEL/TriangulationIntersector.hxx +++ b/src/INTERP_KERNEL/TriangulationIntersector.hxx @@ -23,6 +23,7 @@ #include "PlanarIntersectorP0P1.hxx" #include "PlanarIntersectorP1P0.hxx" #include "PlanarIntersectorP1P1.hxx" +#include "PlanarIntersectorP1P0Bary.hxx" namespace INTERP_KERNEL { @@ -40,6 +41,8 @@ namespace INTERP_KERNEL double intersectGeometry(ConnType icellT, ConnType icellS, ConnType nbNodesT, ConnType nbNodesS); double intersectGeometryWithQuadrangle(const double *quadrangle, const std::vector& sourceCoords, bool isSourceQuad); double intersectGeometryGeneral(const std::vector& targetCoords, const std::vector& sourceCoords); + double intersectGeoBary(const std::vector& targetCell, bool targetCellQuadratic, const double *sourceCell, std::vector& res); + }; } diff --git a/src/INTERP_KERNEL/TriangulationIntersector.txx b/src/INTERP_KERNEL/TriangulationIntersector.txx index ad7e91efb..f200553b7 100644 --- a/src/INTERP_KERNEL/TriangulationIntersector.txx +++ b/src/INTERP_KERNEL/TriangulationIntersector.txx @@ -23,19 +23,25 @@ #include "PlanarIntersectorP0P0.txx" #include "PlanarIntersectorP0P1.txx" #include "PlanarIntersectorP1P0.txx" +#include "PlanarIntersectorP1P0Bary.txx" #include "InterpolationUtils.hxx" #include "PlanarIntersector.hxx" #include +#define TRI_INTERSECTOR TriangulationIntersector +#define TRI_INTER_TEMPLATE template class InterpType> + namespace INTERP_KERNEL { - template class InterpType> - TriangulationIntersector::TriangulationIntersector(const MyMeshType& meshT, const MyMeshType& meshS, - double DimCaracteristic, double Precision, double md3DSurf, - double MedianPlane, int orientation, int PrintLevel) - :InterpType >(meshT,meshS,DimCaracteristic, Precision, md3DSurf, MedianPlane, true, orientation, PrintLevel) + TRI_INTER_TEMPLATE + TRI_INTERSECTOR::TriangulationIntersector(const MyMeshType& meshT, const MyMeshType& meshS, + double DimCaracteristic, double Precision, double md3DSurf, + double MedianPlane, int orientation, int PrintLevel) + :InterpType(meshT,meshS,DimCaracteristic, Precision, md3DSurf, + MedianPlane, true, orientation, PrintLevel) { if(PlanarIntersector::_print_level >= 1) { @@ -44,8 +50,9 @@ namespace INTERP_KERNEL } } - template class InterpType> - double TriangulationIntersector::intersectGeometry(ConnType icellT, ConnType icellS, ConnType nbNodesT, ConnType nbNodesS) + TRI_INTER_TEMPLATE + double TRI_INTERSECTOR::intersectGeometry(ConnType icellT, ConnType icellS, + ConnType nbNodesT, ConnType nbNodesS) { double result = 0.; int orientation = 1; @@ -89,8 +96,10 @@ namespace INTERP_KERNEL return orientation*result; } - template class InterpType> - double TriangulationIntersector::intersectGeometryWithQuadrangle(const double *quadrangle, const std::vector& sourceCoords, bool isSourceQuad) + TRI_INTER_TEMPLATE + double TRI_INTERSECTOR::intersectGeometryWithQuadrangle(const double * quadrangle, + const std::vector& sourceCoords, + bool isSourceQuad) { double result = 0.; ConnType nbNodesS=sourceCoords.size()/SPACEDIM; @@ -129,8 +138,9 @@ namespace INTERP_KERNEL return result; } - template class InterpType> - double TriangulationIntersector::intersectGeometryGeneral(const std::vector& targetCoords, const std::vector& sourceCoords) + TRI_INTER_TEMPLATE + double TRI_INTERSECTOR::intersectGeometryGeneral(const std::vector& targetCoords, + const std::vector& sourceCoords) { double result = 0.; ConnType nbNodesS=sourceCoords.size()/SPACEDIM; @@ -149,14 +159,78 @@ namespace INTERP_KERNEL ConnType nb_inter=((ConnType)inter.size())/2; if(nb_inter >3) inter=reconstruct_polygon(inter); for(ConnType i = 1; i(&inter[0],&inter[2*i],&inter[2*(i+1)],area); - result +=0.5*fabs(area[0]); - } + { + INTERP_KERNEL::crossprod<2>(&inter[0],&inter[2*i],&inter[2*(i+1)],area); + result +=0.5*fabs(area[0]); + } } } return result; } + + //================================================================================ + /*! + * \brief Intersect a triangle and a polygon for P1P0 barycentric algorithm + * \param targetCell - list of coordinates of target polygon in full interlace + * \param targetCellQuadratic - specifies if target polygon is quadratic or not + * \param sourceTria - list of coordinates of source triangle + * \param res - coefficients a,b and c associated to nodes of sourceTria + */ + //================================================================================ + + TRI_INTER_TEMPLATE + double TRI_INTERSECTOR::intersectGeoBary(const std::vector& targetCell, + bool targetCellQuadratic, + const double * sourceTria, + std::vector& res) + { + std::vector sourceCell(3); + sourceCell[0] = &sourceTria[0]; + sourceCell[1] = &sourceTria[SPACEDIM]; + sourceCell[2] = &sourceTria[SPACEDIM*2]; + + //Compute the intersection area + double inter_area[SPACEDIM], total_area = 0.; + double total_barycenter[SPACEDIM]; + total_barycenter[0]=total_barycenter[1] = 0.; + + const ConnType nbNodesT=targetCell.size()/SPACEDIM; + for(ConnType iT = 1; iT inter; + INTERP_KERNEL::intersec_de_triangle(&targetCell[0],&targetCell[SPACEDIM*iT],&targetCell[SPACEDIM*(iT+1)], + sourceCell[0], sourceCell[1], sourceCell[2], + inter, PlanarIntersector::_dim_caracteristic, + PlanarIntersector::_precision); + ConnType nb_inter=((ConnType)inter.size())/2; + if(nb_inter >3) inter=reconstruct_polygon(inter); + for(ConnType i = 1; i(&inter[0],&inter[2*i],&inter[2*(i+1)],inter_area); + inter_area[0] = 0.5 * fabs( inter_area[0] ); + total_area += inter_area[0]; + std::vector inter_bary=INTERP_KERNEL::bary_poly(inter); + total_barycenter[0] += inter_area[0] * inter_bary[0]; + total_barycenter[1] += inter_area[0] * inter_bary[1]; + } + } + if ( total_area > std::numeric_limits::min() ) + { + total_barycenter[0] /= total_area; + total_barycenter[1] /= total_area; + res.resize(3); + barycentric_coords( sourceCell, &total_barycenter[0], &res[0]); + res[0] *= total_area; + res[1] *= total_area; + res[2] *= total_area; + } + else + { + total_area = 0; + } + return total_area; + } + } #endif -- 2.39.2