src/INTERP_KERNEL/Geometric2DIntersector.txx

   1 // Copyright (C) 2007-2016  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // Author : Anthony Geay (CEA/DEN)
  20 #ifndef __GEOMETRIC2DINTERSECTOR_TXX__
  21 #define __GEOMETRIC2DINTERSECTOR_TXX__
  22
  23 #include "Geometric2DIntersector.hxx"
  24 #include "PlanarIntersectorP0P0.txx"
  25 #include "Planar2D1DIntersectorP0P0.txx"
  26 #include "PlanarIntersectorP0P1.txx"
  27 #include "PlanarIntersectorP1P0.txx"
  28 #include "PlanarIntersectorP1P1.txx"
  29 #include "PlanarIntersectorP1P0Bary.txx"
  30 #include "PlanarIntersectorP0P1Bary.txx"
  31 #include "CellModel.hxx"
  32
  33 #include "InterpKernelGeo2DQuadraticPolygon.hxx"
  34 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
  35 #include "InterpKernelGeo2DEdgeLin.hxx"
  36 #include "InterpKernelGeo2DNode.hxx"
  37
  38 #define GEO2D_INTERSECTOR    Geometric2DIntersector<MyMeshType,MyMatrix,InterpType>
  39 #define INTERSECTOR_TEMPLATE template<class MyMeshType, class MyMatrix, template <class MeshType, class TheMatrix, class ThisIntersector> class InterpType>
  40
  41 namespace INTERP_KERNEL
  42 {
  43   INTERSECTOR_TEMPLATE
  44   GEO2D_INTERSECTOR::Geometric2DIntersector(const MyMeshType& meshT, const MyMeshType& meshS,
  45                                             double dimCaracteristic, double md3DSurf, double minDot3DSurf, double medianPlane,
  46                                             double precision, int orientation):
  47     InterpType<MyMeshType,MyMatrix,GEO2D_INTERSECTOR >(meshT,meshS,dimCaracteristic, precision, md3DSurf, minDot3DSurf, medianPlane, true, orientation, 0),
  48     _precision(precision)
  49   {
  50   }
  51
  52   INTERSECTOR_TEMPLATE
  53   double GEO2D_INTERSECTOR::intersectGeometry(ConnType icellT,   ConnType icellS,
  54                                               ConnType nbNodesT, ConnType nbNodesS)
  55   {
  56     int orientation = 1;
  57     std::vector<double> CoordsT;
  58     std::vector<double> CoordsS;
  59     PlanarIntersector<MyMeshType,MyMatrix>::getRealCoordinates(icellT,icellS,nbNodesT,nbNodesS,CoordsT,CoordsS,orientation);
  60     NormalizedCellType tT=PlanarIntersector<MyMeshType,MyMatrix>::_meshT.getTypeOfElement(icellT);
  61     NormalizedCellType tS=PlanarIntersector<MyMeshType,MyMatrix>::_meshS.getTypeOfElement(icellS);
  62     QuadraticPolygon *p1=buildPolygonFrom(CoordsT,tT);
  63     QuadraticPolygon *p2=buildPolygonFrom(CoordsS,tS);
  64     double ret=p1->intersectWithAbs(*p2);
  65     delete p1; delete p2;
  66     return orientation*ret;
  67   }
  68
  69   INTERSECTOR_TEMPLATE
  70   double GEO2D_INTERSECTOR::intersectGeometry1D(ConnType icellT,   ConnType icellS,
  71                                                 ConnType nbNodesT, ConnType nbNodesS,
  72                                                 bool& isColinear)
  73   {
  74     int orientation = 1;
  75     std::vector<double> CoordsT;
  76     std::vector<double> CoordsS;
  77     PlanarIntersector<MyMeshType,MyMatrix>::getRealCoordinates(icellT,icellS,nbNodesT,nbNodesS,CoordsT,CoordsS,orientation);
  78     NormalizedCellType tT=PlanarIntersector<MyMeshType,MyMatrix>::_meshT.getTypeOfElement(icellT);
  79     NormalizedCellType tS=PlanarIntersector<MyMeshType,MyMatrix>::_meshS.getTypeOfElement(icellS);
  80     QuadraticPolygon *p1=buildPolygonFrom(CoordsT,tT);
  81     QuadraticPolygon *p2=buildPolygonOfOneEdgeFrom(CoordsS,tS);
  82     double ret=p1->intersectWithAbs1D(*p2, isColinear);
  83     delete p1; delete p2;
  84     return orientation*ret;
  85   }
  86
  87   INTERSECTOR_TEMPLATE
  88   double GEO2D_INTERSECTOR::intersectGeometryWithQuadrangle(const double             * quadrangle,
  89                                                             const std::vector<double>& sourceCoords,
  90                                                             bool                       isSourceQuad)
  91   {
  92     std::vector<Node *> nodes(4);
  93     nodes[0]=new Node(quadrangle[0],quadrangle[1]);
  94     nodes[1]=new Node(quadrangle[SPACEDIM],quadrangle[SPACEDIM+1]);
  95     nodes[2]=new Node(quadrangle[2*SPACEDIM],quadrangle[2*SPACEDIM+1]);
  96     nodes[3]=new Node(quadrangle[3*SPACEDIM],quadrangle[3*SPACEDIM+1]);
  97     int nbOfSourceNodes=sourceCoords.size()/SPACEDIM;
  98     std::vector<Node *> nodes2(nbOfSourceNodes);
  99     for(int i=0;i<nbOfSourceNodes;i++)
 100       nodes2[i]=new Node(sourceCoords[i*SPACEDIM],sourceCoords[i*SPACEDIM+1]);
 101     QuadraticPolygon *p1=QuadraticPolygon::BuildLinearPolygon(nodes);
 102     QuadraticPolygon *p2;
 103     if(!isSourceQuad)
 104       p2=QuadraticPolygon::BuildLinearPolygon(nodes2);
 105     else
 106       p2=QuadraticPolygon::BuildArcCirclePolygon(nodes2);
 107     double ret=p1->intersectWithAbs(*p2);
 108     delete p1; delete p2;
 109     return ret;
 110   }
 111
 112   INTERSECTOR_TEMPLATE
 113   double GEO2D_INTERSECTOR::intersectGeometryGeneral(const std::vector<double>& targetCoords,
 114                                                      const std::vector<double>& sourceCoords)
 115   {
 116     int nbOfTargetNodes=targetCoords.size()/SPACEDIM;
 117     std::vector<Node *> nodes(nbOfTargetNodes);
 118     for(int i=0;i<nbOfTargetNodes;i++)
 119       nodes[i]=new Node(targetCoords[i*SPACEDIM],targetCoords[i*SPACEDIM+1]);
 120     int nbOfSourceNodes=sourceCoords.size()/SPACEDIM;
 121     std::vector<Node *> nodes2(nbOfSourceNodes);
 122     for(int i=0;i<nbOfSourceNodes;i++)
 123       nodes2[i]=new Node(sourceCoords[i*SPACEDIM],sourceCoords[i*SPACEDIM+1]);
 124     QuadraticPolygon *p1=QuadraticPolygon::BuildLinearPolygon(nodes);
 125     QuadraticPolygon *p2=QuadraticPolygon::BuildLinearPolygon(nodes2);
 126     double ret=p1->intersectWithAbs(*p2);
 127     delete p1; delete p2;
 128     return ret;
 129   }
 130
 131   //================================================================================
 132   /*!
 133    * \brief Intersect a triangle and a polygon for P1P0 barycentric algorithm
 134    *  \param targetCell - list of coordinates of target polygon in full interlace
 135    *  \param targetCellQuadratic - specifies if target polygon is quadratic or not
 136    *  \param sourceTria - list of coordinates of source triangle
 137    *  \param res - coefficients a,b and c associated to nodes of sourceTria
 138    */
 139   //================================================================================
 140
 141   INTERSECTOR_TEMPLATE
 142   double GEO2D_INTERSECTOR::intersectGeoBary(const std::vector<double>& targetCell,
 143                                              bool                       targetCellQuadratic,
 144                                              const double *             sourceTria,
 145                                              std::vector<double>&       res)
 146   {
 147     std::vector<Node *> nodes(3);
 148     nodes[0]=new Node(sourceTria[0*SPACEDIM],sourceTria[0*SPACEDIM+1]);
 149     nodes[1]=new Node(sourceTria[1*SPACEDIM],sourceTria[1*SPACEDIM+1]);
 150     nodes[2]=new Node(sourceTria[2*SPACEDIM],sourceTria[2*SPACEDIM+1]);
 151     int nbOfTargetNodes=targetCell.size()/SPACEDIM;
 152     std::vector<Node *> nodes2(nbOfTargetNodes);
 153     for(int i=0;i<nbOfTargetNodes;i++)
 154       nodes2[i]=new Node(targetCell[i*SPACEDIM],targetCell[i*SPACEDIM+1]);
 155     QuadraticPolygon *p1=QuadraticPolygon::BuildLinearPolygon(nodes);
 156     QuadraticPolygon *p2;
 157     if(!targetCellQuadratic)
 158       p2=QuadraticPolygon::BuildLinearPolygon(nodes2);
 159     else
 160       p2=QuadraticPolygon::BuildArcCirclePolygon(nodes2);
 161     double barycenter[2];
 162     double ret=p1->intersectWithAbs(*p2,barycenter);
 163     delete p1; delete p2;
 164     if ( ret > std::numeric_limits<double>::min() )
 165     {
 166       std::vector<const double* > sourceCell(3);
 167       sourceCell[0] = &sourceTria[0];
 168       sourceCell[1] = &sourceTria[SPACEDIM];
 169       sourceCell[2] = &sourceTria[SPACEDIM*2];
 170       res.resize(3);
 171       barycentric_coords( sourceCell, barycenter, &res[0]);
 172       res[0] *= ret;
 173       res[1] *= ret;
 174       res[2] *= ret;
 175     }
 176     else
 177     {
 178       ret = 0;
 179     }
 180     return ret;
 181   }
 182
 183   INTERSECTOR_TEMPLATE
 184   QuadraticPolygon *GEO2D_INTERSECTOR::buildPolygonFrom(const std::vector<double>& coords, NormalizedCellType type)
 185   {
 186     int nbNodes=coords.size()/SPACEDIM;
 187     std::vector<Node *> nodes(nbNodes);
 188     for(int i=0;i<nbNodes;i++)
 189       nodes[i]=new Node(coords[i*SPACEDIM],coords[i*SPACEDIM+1]);
 190     if(!CellModel::GetCellModel(type).isQuadratic())
 191       return QuadraticPolygon::BuildLinearPolygon(nodes);
 192     else
 193       return QuadraticPolygon::BuildArcCirclePolygon(nodes);
 194   }
 195
 196   INTERSECTOR_TEMPLATE
 197   QuadraticPolygon *GEO2D_INTERSECTOR::buildPolygonOfOneEdgeFrom(const std::vector<double>& coords, NormalizedCellType type)
 198   {
 199     if(type==NORM_SEG2)
 200       {
 201         Node *node0=new Node(coords[0],coords[1]);
 202         Node *node1=new Node(coords[SPACEDIM],coords[SPACEDIM+1]);
 203         QuadraticPolygon *ret=new QuadraticPolygon;
 204         ret->pushBack(new EdgeLin(node0,node1));
 205         node0->decrRef(); node1->decrRef();
 206         return ret;
 207       }
 208     else if(type==NORM_SEG3)
 209       {
 210         Node *nodeBg=new Node(coords[0],coords[1]);
 211         Node *nodeEnd=new Node(coords[SPACEDIM],coords[SPACEDIM+1]);
 212         Node *nodeMiddle=new Node(coords[2*SPACEDIM],coords[2*SPACEDIM+1]);
 213         QuadraticPolygon *ret=new QuadraticPolygon;
 214         ret->pushBack(new EdgeArcCircle(nodeBg,nodeMiddle,nodeEnd));
 215         nodeBg->decrRef(); nodeEnd->decrRef(); nodeMiddle->decrRef();
 216         return ret;
 217       }
 218     else
 219       throw INTERP_KERNEL::Exception("buildPolygonOfOneEdgeFrom : trying to build such non close QuadraticPolygon with 1D type !");
 220   }
 221
 222   INTERSECTOR_TEMPLATE
 223   QuadraticPolygon *GEO2D_INTERSECTOR::buildPolygonAFrom(ConnType cell, int nbOfPoints, NormalizedCellType type)
 224   {
 225     const ConnType *startOfCellNodeConn=PlanarIntersector<MyMeshType,MyMatrix>::_connectT+OTT<ConnType,numPol>::conn2C(PlanarIntersector<MyMeshType,MyMatrix>::_connIndexT[OTT<ConnType,numPol>::ind2C(cell)]);
 226     std::vector<Node *> nodes(nbOfPoints);
 227     for(int i=0;i<nbOfPoints;i++)
 228       nodes[i]=new Node(PlanarIntersector<MyMeshType,MyMatrix>::_coordsT+OTT<ConnType,numPol>::coo2C(startOfCellNodeConn[i])*SPACEDIM);
 229     if(CellModel::GetCellModel(type).isQuadratic())
 230       return QuadraticPolygon::BuildLinearPolygon(nodes);
 231     else
 232       return QuadraticPolygon::BuildArcCirclePolygon(nodes);
 233   }
 234
 235   INTERSECTOR_TEMPLATE
 236   QuadraticPolygon *GEO2D_INTERSECTOR::buildPolygonBFrom(ConnType cell, int nbOfPoints, NormalizedCellType type)
 237   {
 238     const ConnType *startOfCellNodeConn=PlanarIntersector<MyMeshType,MyMatrix>::_connectS+OTT<ConnType,numPol>::conn2C(PlanarIntersector<MyMeshType,MyMatrix>::_connIndexS[OTT<ConnType,numPol>::ind2C(cell)]);
 239     std::vector<Node *> nodes(nbOfPoints);
 240     for(int i=0;i<nbOfPoints;i++)
 241       nodes[i]=new Node(PlanarIntersector<MyMeshType,MyMatrix>::_coordsS+OTT<ConnType,numPol>::coo2C(startOfCellNodeConn[i])*SPACEDIM);
 242     const CellModel& cm=CellModel::GetCellModel(type);
 243     if(!cm.isQuadratic())
 244       return QuadraticPolygon::BuildLinearPolygon(nodes);
 245     else
 246       return QuadraticPolygon::BuildArcCirclePolygon(nodes);
 247   }
 248 }
 249
 250 #endif