--- /dev/null
+// Copyright (C) 2018 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, or (at your option) any later version.
+//
+// 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
+//
+// Author : Anthony Geay (EDF R&D)
+
+#ifndef __INTERPKERNELASSERT_HXX__
+#define __INTERPKERNELASSERT_HXX__
+
+#include "InterpKernelException.hxx"
+
+#include <sstream>
+
+#define IKAssert(a) { bool verdict(a); \
+ if(!verdict) { std::ostringstream osszz; osszz << "Assertion \"" << #a << "\" failed into " << __FILE__ << " at line " << __LINE__ << " !"; throw INTERP_KERNEL::Exception(osszz.str()); } }
+
+#define IKAssertMsg(a,b) { bool verdict(a); \
+ if(!verdict) { std::ostringstream osszz; osszz << "Assertion \"" << #a << "\" failed into " << __FILE__ << " at line " << __LINE__ << " with message \"" << b << "\" !"; throw INTERP_KERNEL::Exception(osszz.str()); } }
+
+#endif
void BoundingBox::toCompactData(double data[6]) const
{
- std::copy(_coords,_coords+6,data);
+ data[0]=_coords[XMIN];
+ data[1]=_coords[XMAX];
+ data[2]=_coords[YMIN];
+ data[3]=_coords[YMAX];
+ data[4]=_coords[ZMIN];
+ data[5]=_coords[ZMAX];
}
}
#include "BBTree.txx"
+#include "InterpKernelAssert.hxx"
+
namespace INTERP_KERNEL
{
/**
{
//#if __cplusplus >= 201103L
constexpr int SPACEDIM=MyMeshType::MY_SPACEDIM;
+ using ConnType=typename MyMeshType::MyConnType;
+ IKAssert(SPACEDIM==3);
if(InterpolationOptions::getIntersectionType() != PointLocator)
INTERP_KERNEL::Exception("Invalid 1D/0D intersection type specified : must be PointLocator.");
std::string methC ( InterpolationOptions::filterInterpolationMethod(method) );
- if(methC=="P1P1")
+ if(methC!="P1P1")
throw Exception("Invalid method chosen must be in \"P1P1\".");
-
- typedef typename MyMeshType::MyConnType ConnType;
+ const double epsilon(getPrecision());
// create MeshElement objects corresponding to each element of the two meshes
const unsigned long numSrcElems(srcMesh.getNumberOfElements()), numTargetElems(targetMesh.getNumberOfElements());
}
adjustBoundingBoxes(bboxes);
- const double *bboxPtr=0;
+ const double *bboxPtr(nullptr);
if(numSrcElems>0)
bboxPtr=&bboxes[0];
- BBTree<3,ConnType> tree(bboxPtr, srcElemIdx, 0, numSrcElems);
-
- // for each target element, get source elements with which to calculate intersection
- // - calculate intersection by calling intersectCells
+ BBTree<SPACEDIM,ConnType> tree(bboxPtr, srcElemIdx, 0, numSrcElems);
+ const ConnType *trgConnPtr(targetMesh.getConnectivityPtr()),*trgConnIPtr(targetMesh.getConnectivityIndexPtr());
+ const ConnType *srcConnPtr(srcMesh.getConnectivityPtr()),*srcConnIPtr(srcMesh.getConnectivityIndexPtr());
+ const double *trgCooPtr(targetMesh.getCoordinatesPtr()),*srcCooPtr(srcMesh.getCoordinatesPtr());
for(unsigned long i = 0; i < numTargetElems; ++i)
{
- std::vector<ConnType> intersectElems;
-
- //tree.getElementsAroundPoint(targetBox, intersectElems);
-
- //if ( !intersectElems.empty() )
- // intersector->intersectCells(targetIdx,intersectElems,result);
+ IKAssert(trgConnIPtr[i+1]==i+1 && trgConnIPtr[i]==i);
+ std::vector<ConnType> srcSegCondidates;
+ const double *trgCellPosition(trgCooPtr+SPACEDIM*trgConnPtr[i]);
+ typename MatrixType::value_type& resRow(result[i]);
+ tree.getElementsAroundPoint(trgCellPosition, srcSegCondidates);
+ for(auto srcSeg: srcSegCondidates)
+ {
+ IKAssertMsg(srcConnIPtr[srcSeg+1]==2*(srcSeg+1) && srcConnIPtr[srcSeg]==2*srcSeg,"Only implemented for linear 1D source");
+ double bc0(0.),bc1(0.);
+ ConnType srcNode0(srcConnPtr[2*srcSeg]),srcNode1(srcConnPtr[2*srcSeg+1]);
+ if(IsPointOn3DSeg(srcCooPtr+SPACEDIM*srcNode0,srcCooPtr+SPACEDIM*srcNode1,trgCellPosition,epsilon,bc0,bc1))
+ {
+ resRow.insert(std::make_pair(srcNode0,bc0));
+ resRow.insert(std::make_pair(srcNode1,bc1));
+ continue;
+ }
+ }
}
-
- // free allocated memory
delete [] srcElemIdx;
-
for(unsigned long i = 0 ; i < numSrcElems ; ++i)
- {
- delete srcElems[i];
- }
+ delete srcElems[i];
return srcMesh.getNumberOfNodes();
-
}
//#endif
}
/* calcul la surface d'un triangle */
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
- inline double Surf_Tri(const double* P_1,const double* P_2,const double* P_3)
+ inline double Surf_Tri(const double *P_1,const double *P_2,const double *P_3)
{
double A=(P_3[1]-P_1[1])*(P_2[0]-P_1[0])-(P_2[1]-P_1[1])*(P_3[0]-P_1[0]);
double Surface = 0.5*fabs(A);
//fonction qui calcul le determinant des vecteurs: P3P1 et P3P2
//(cf doc CGAL).
- inline double mon_determinant(const double* P_1,
- const double* P_2,
- const double* P_3)
+ inline double mon_determinant(const double *P_1,
+ const double *P_2,
+ const double *P_3)
{
double mon_det=(P_1[0]-P_3[0])*(P_2[1]-P_3[1])-(P_2[0]-P_3[0])*(P_1[1]-P_3[1]);
return mon_det;
{
double X=P_1[0]-P_2[0];
double Y=P_1[1]-P_2[1];
- double norme=sqrt(X*X+Y*Y);
- return norme;
+ return sqrt(X*X+Y*Y);
}
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
}
}
+ /*!
+ * Compute barycentric coords of point \a point relative to segment S [segStart,segStop] in 3D space.
+ * If point \a point is further from S than eps false is returned.
+ * If point \a point projection on S is outside S false is also returned.
+ * If point \a point is closer from S than eps and its projection inside S true is returned and \a bc0 and \a bc1 output parameter set.
+ */
+ inline bool IsPointOn3DSeg(const double segStart[3], const double segStop[3], const double point[3], double eps, double& bc0, double& bc1)
+ {
+ double AB[3]={segStop[0]-segStart[0],segStop[1]-segStart[1],segStop[2]-segStart[2]},AP[3]={point[0]-segStart[0],point[1]-segStart[1],point[2]-segStart[2]};
+ double l_AB(sqrt(AB[0]*AB[0]+AB[1]*AB[1]+AB[2]*AB[2]));
+ double AP_dot_AB((AP[0]*AB[0]+AP[1]*AB[1]+AP[2]*AB[2])/(l_AB*l_AB));
+ double projOfPOnAB[3]={segStart[0]+AP_dot_AB*AB[0],segStart[1]+AP_dot_AB*AB[1],segStart[2]+AP_dot_AB*AB[2]};
+ double V_dist_P_AB[3]={point[0]-projOfPOnAB[0],point[1]-projOfPOnAB[1],point[2]-projOfPOnAB[2]};
+ double dist_P_AB(sqrt(V_dist_P_AB[0]*V_dist_P_AB[0]+V_dist_P_AB[1]*V_dist_P_AB[1]+V_dist_P_AB[2]*V_dist_P_AB[2]));
+ if(dist_P_AB>=eps)
+ return false;//to far from segment [segStart,segStop]
+ if(AP_dot_AB<-eps || AP_dot_AB>1.+eps)
+ return false;
+ AP_dot_AB=std::max(AP_dot_AB,0.); AP_dot_AB=std::min(AP_dot_AB,1.);
+ bc0=1.-AP_dot_AB; bc1=AP_dot_AB;
+ return true;
+ }
+
/*!
* Calculate pseudo barycentric coordinates of a point p with respect to the quadrangle vertices.
* This method makes the assumption that:
#define __MEDCOUPLINGNORMALIZEDUNSTRUCTUREDMESH_TXX__
#include "MEDCouplingNormalizedUnstructuredMesh.hxx"
+#include "InterpKernelAssert.hxx"
#include "MEDCouplingUMesh.hxx"
#include "MEDCoupling1GTUMesh.hxx"
template<int SPACEDIM,int MESHDIM>
void MEDCouplingNormalizedUnstructuredMesh<SPACEDIM,MESHDIM>::prepare()
{
+ IKAssert(_mesh->getSpaceDimension()==SPACEDIM);
const MEDCoupling::MEDCouplingUMesh *m1(dynamic_cast<const MEDCoupling::MEDCouplingUMesh *>(_mesh));
if(m1)
{
sourceMesh->allocateCells(0);
sourceMesh->finishInsertingCells();
DataArrayDouble *myCoords=DataArrayDouble::New();
- myCoords->alloc(0,0);
+ myCoords->alloc(0,2);
sourceMesh->setCoords(myCoords);
myCoords->decrRef();
MEDCouplingUMesh *targetMesh=MEDCouplingUMesh::New();
