-// Copyright (C) 2007-2016 CEA/DEN, EDF R&D
+// Copyright (C) 2007-2019 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
// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
//
// Author : Anthony Geay (CEA/DEN)
-#ifndef __CURVEINTERSECTOR_TXX__
-#define __CURVEINTERSECTOR_TXX__
+
+#pragma once
#include "CurveIntersector.hxx"
#include "InterpolationUtils.hxx"
+#include "PointLocatorAlgos.txx"
#include <limits>
int ibox=0;
for(long icell=0; icell<nbelems; icell++)
{
- int nb_nodes_per_elem = conn_index[icell+1]-conn_index[icell];
+ ConnType nb_nodes_per_elem = conn_index[icell+1]-conn_index[icell];
//initializing bounding box limits
for(int idim=0; idim<SPACEDIM; idim++)
{
bbox[2*SPACEDIM*ibox+2*idim+1] = -std::numeric_limits<double>::max();
}
//updating the bounding box with each node of the element
- for (int j=0; j<nb_nodes_per_elem; j++)
+ for (ConnType j=0; j<nb_nodes_per_elem; j++)
{
const double* coord_node = coords +
SPACEDIM*OTT<ConnType,numPol>
}
/*!
- Computes the bouding box of a given element. iP in numPol mode.
+ Computes the bounding box of a given element. iP in numPol mode.
*/
template<class MyMeshType, class MyMatrix>
void CurveIntersector<MyMeshType,MyMatrix>::getElemBB (double* bb,
/*!
* \param [in] startOfSeg - input coming from intersectSegments or intersectSegmentsInternal
- * \param [in] endOfSeg - input coming from intersectSegments or intersectSegmentsInternal. Assume that endOfSeg>startOfSeg.
+ * \param [in] endOfSeg - input coming from intersectSegments or intersectSegmentsInternal. NO Assume about sort
* \param [in] pt - position of point that the method computes the bary coords for.
*/
template<class MyMeshType, class MyMatrix>
- bool CurveIntersector<MyMeshType,MyMatrix>::ComputeBaryCoordsOf(double startOfSeg, double endOfSeg, double pt, double& startPos, double& endPos)
+ void CurveIntersector<MyMeshType,MyMatrix>::ComputeBaryCoordsOf(double startOfSeg, double endOfSeg, double pt, double& startPos, double& endPos)
{
double deno(endOfSeg-startOfSeg);
- startPos=(endOfSeg-pt)/deno;
- endPos=1.-startPos;
- return startPos>=0. && endPos>=0.;
+ startPos = (endOfSeg-pt)/deno;
+ startPos = std::max(startPos,0.); startPos = std::min(startPos,1.);
+ endPos=1.-startPos;
}
-
+
/*! Readjusts a set of bounding boxes so that they are extended
in all dimensions for avoiding missing interesting intersections
void CurveIntersector<MyMeshType,MyMatrix>::adjustBoundingBoxes (std::vector<double>& bbox,
double adjustmentEpsAbs)
{
- long size = bbox.size()/(2*SPACEDIM);
- for (int i=0; i<size; i++)
+ std::size_t size = bbox.size()/(2*SPACEDIM);
+ for (std::size_t i=0; i<size; i++)
{
for(int idim=0; idim<SPACEDIM; idim++)
{
template<class MyMeshType, class MyMatrix>
bool CurveIntersector<MyMeshType,MyMatrix>::getRealTargetCoordinates(ConnType icellT, std::vector<double>& coordsT) const
{
- int nbNodesT(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)+1] - _connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]);
+ ConnType nbNodesT(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)+1] - _connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]);
coordsT.resize(SPACEDIM*nbNodesT);
for (ConnType iT=0; iT<nbNodesT; iT++)
{
template<class MyMeshType, class MyMatrix>
typename MyMeshType::MyConnType CurveIntersector<MyMeshType,MyMatrix>::getNodeIdOfTargetCellAt(ConnType icellT, ConnType nodeIdInCellT) const
{
- int nbNodesT(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)+1] - _connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]);
+ ConnType nbNodesT(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)+1] - _connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]);
if(nodeIdInCellT>=0 && nodeIdInCellT<nbNodesT)
return OTT<ConnType,numPol>::coo2C(_connectT[OTT<ConnType,numPol>::conn2C(_connIndexT[OTT<ConnType,numPol>::ind2C(icellT)]+nodeIdInCellT)]);
else
template<class MyMeshType, class MyMatrix>
bool CurveIntersector<MyMeshType,MyMatrix>::getRealSourceCoordinates(ConnType icellS, std::vector<double>& coordsS) const
{
- int nbNodesS = _connIndexS[OTT<ConnType,numPol>::ind2C(icellS)+1] - _connIndexS[OTT<ConnType,numPol>::ind2C(icellS)];
+ ConnType nbNodesS = _connIndexS[OTT<ConnType,numPol>::ind2C(icellS)+1] - _connIndexS[OTT<ConnType,numPol>::ind2C(icellS)];
coordsS.resize(SPACEDIM*nbNodesS);
for(ConnType iS=0; iS<nbNodesS; iS++)
{
template<class MyMeshType, class MyMatrix>
typename MyMeshType::MyConnType CurveIntersector<MyMeshType,MyMatrix>::getNodeIdOfSourceCellAt(ConnType icellS, ConnType nodeIdInCellS) const
{
- int nbNodesS(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)+1] - _connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]);
+ ConnType nbNodesS(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)+1] - _connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]);
if(nodeIdInCellS>=0 && nodeIdInCellS<nbNodesS)
return OTT<ConnType,numPol>::coo2C(_connectS[OTT<ConnType,numPol>::conn2C(_connIndexS[OTT<ConnType,numPol>::ind2C(icellS)]+nodeIdInCellS)]);
else
std::vector<TDualSegment>& segments)
{
// get coordinates of cell nodes
- int nbNodes;
- std::vector<double> ncoords;
- std::vector<int> nodeIds;
+ ConnType nbNodes;
+ std::vector<double> ncoords;
+ std::vector<ConnType> nodeIds;
{
const ConnType *connect = mesh.getConnectivityPtr();
const ConnType *connIndex = mesh.getConnectivityIndexPtr();
ncoords.resize(SPACEDIM*nbNodes);
nodeIds.resize(nbNodes);
- for(int i=0; i<nbNodes; i++)
+ for(ConnType i=0; i<nbNodes; i++)
for(int idim=0; idim<SPACEDIM; idim++)
{
nodeIds[i] = connect[OTT<ConnType,numPol>::conn2C(connIndex[OTT<ConnType,numPol>::ind2C(icell)]+i)];
// fill segments
segments.clear();
segments.reserve( 2*nbNodes );
- for(int i=0; i<nbNodes-1; i++)
+ for(ConnType i=0; i<nbNodes-1; i++)
{
segments.push_back(TDualSegment());
TDualSegment& seg1 = segments.back();
}
}
}
+
+ template<class MyMeshType, class MyMatrix>
+ bool CurveIntersector<MyMeshType,MyMatrix>::projectionThis(const double *coordsT, const double *coordsS, double& xs0, double& xs1, double& xt) const
+ {
+ enum { X=0, Y };
+ switch(SPACEDIM)
+ {
+ case 1:
+ {
+ xt = coordsT[0];
+ xs0 = coordsS[0]; xs1 = coordsS[1];
+ return true;
+ }
+ case 2:
+ {
+ const double *s0(coordsS),*s1(coordsS + 2);
+ double s01[2] = { s1[X]-s0[X], s1[Y]-s0[Y] }; // src segment direction
+ double sSize = sqrt( s01[X]*s01[X] + s01[Y]*s01[Y] ); // src segment size
+ if( sSize < this->_precision )
+ return false;
+ s01[X] /= sSize; s01[Y] /= sSize; // normalize s01
+ double t[2] = { coordsT[X]-s0[X], coordsT[Y]-s0[Y] };
+ xs0 = 0. ; xs1 = sSize; xt = s01[X]*t[X] + s01[Y]*t[Y];
+ double proj_t_on_s[2] = { s0[X]+xt*s01[X], s0[Y]+xt*s01[Y] };
+ double dist_t_s_vect[2] = { coordsT[X]-proj_t_on_s[X], coordsT[Y]-proj_t_on_s[Y] };
+ double dist_t_s = sqrt( dist_t_s_vect[X]*dist_t_s_vect[X]+dist_t_s_vect[Y]*dist_t_s_vect[Y] );
+ return dist_t_s < this->_precision;
+ }
+ default:
+ throw Exception("CurveIntersector::projectionThis : space dimension of mesh must be 1 or 2");
+ }
+ }
template<class MyMeshType, class MyMatrix>
bool CurveIntersector<MyMeshType,MyMatrix>::projectionThis(const double *coordsT, const double *coordsS,
double x1 = std::min( xt1, xs1 );
return ( x0 < x1 ) ? ( x1 - x0 ) : 0.;
}
+
+ template<class MyMeshType>
+ class DummyMyMeshType1D
+ {
+ public:
+ static const int MY_SPACEDIM=1;
+ static const int MY_MESHDIM=8;
+ typedef mcIdType MyConnType;
+ static const INTERP_KERNEL::NumberingPolicy My_numPol=MyMeshType::My_numPol;
+ // begin
+ // useless, but for windows compilation ...
+ const double *getCoordinatesPtr() const { return nullptr; }
+ const MyConnType *getConnectivityPtr() const { return nullptr; }
+ const MyConnType *getConnectivityIndexPtr() const { return nullptr; }
+ INTERP_KERNEL::NormalizedCellType getTypeOfElement(MyConnType) const { return (INTERP_KERNEL::NormalizedCellType)0; }
+ // end
+ };
+
+ /*!
+ * This method determines if a target point ( \a coordsT ) is in source seg2 contained in \a coordsS. To do so _precision attribute is used.
+ * If target point is in, \a xs0, \a xs1 and \a xt are set to 1D referential for a further barycentric computation.
+ * This method deals with SPACEDIM == 2 (see projectionThis).
+ */
+ template<class MyMeshType, class MyMatrix>
+ bool CurveIntersector<MyMeshType,MyMatrix>::isPtIncludedInSeg(const double *coordsT, const double *coordsS, double& xs0, double& xs1, double& xt) const
+ {
+ if(!projectionThis(coordsT,coordsS,xs0,xs1,xt))
+ return false;
+ constexpr ConnType TAB[2]={0,1};
+ const double coordsS_1D[2]={xs0,xs1};
+ const double *coordsT_1D(&xt);
+ return PointLocatorAlgos<DummyMyMeshType1D<MyMeshType>>::isElementContainsPoint(coordsT_1D,NORM_SEG2,coordsS_1D,TAB,2,this->_precision);
+ }
/*!
* \brief Return length of intersection of two segments
double xs0,xs1,xt0,xt1;
return intersectSegmentsInternal(coordsT,coordsS,xs0,xs1,xt0,xt1);
}
-
}
-
-#endif
