}
- _max_left=max_left;
- _min_right=min_right;
+ _max_left=max_left+1e-12;
+ _min_right=min_right-1e-12;
_left=new BBTree(bbs, &(new_elems_left[0]), level+1, new_elems_left.size());
_right=new BBTree(bbs, &(new_elems_right[0]), level+1, new_elems_right.size());
_right->getIntersectingElems(bb,elems);
}
+
+ /*! returns in \a elems the list of elements potentially containing the point pointed to by \a xx
+ \param xx pointer to query point coords
+ \param elems list of elements (given in 0-indexing) intersecting the bounding box
+ */
+
+ void getElementsAroundPoint(const double* xx, std::vector<int>& elems) const
+ {
+ // terminal node : return list of elements intersecting bb
+ if (_terminal)
+ {
+ for (int i=0; i<_nbelems; i++)
+ {
+ const double* const bb_ptr=_bb+_elems[i]*2*dim;
+ bool intersects = true;
+ for (int idim=0; idim<dim; idim++)
+ {
+ if (bb_ptr[idim*2]-xx[idim]>1e-12|| bb_ptr[idim*2+1]-xx[idim]<-1e-12)
+ intersects=false;
+ }
+ if (intersects)
+ {
+ elems.push_back(_elems[i]);
+ }
+ }
+ return;
+ }
+
+ //non terminal node
+ if (xx[_level%dim] < _min_right)
+ {
+ _left->getElementsAroundPoint(xx, elems);
+ return;
+ }
+ if (xx[_level%dim]> _max_left)
+ {
+ _right->getElementsAroundPoint(xx,elems);
+ return;
+ }
+ _left->getElementsAroundPoint(xx,elems);
+ _right->getElementsAroundPoint(xx,elems);
+ }
+
int size()
return _left->size()+_right->size();
}
};
-
#endif
* intersection matrix.
*
* The matrix is partially sparse : it is a vector of maps of integer - double pairs.
+ * It can also be an INTERP_KERNEL::Matrix object.
* The length of the vector is equal to the number of target elements - for each target element there is a map, regardless
* of whether the element intersects any source elements or not. But in the maps there are only entries for those source elements
* which have a non-zero intersection volume with the target element. The vector has indices running from
- * 0 to (#target elements - 1), meaning that the map for target element i is stored at index i - 1. In the maps, however,
+ * 0 to (nb target elements - 1), meaning that the map for target element i is stored at index i - 1. In the maps, however,
* the indexing is more natural : the intersection volume of the target element i with source element j is found at matrix[i-1][j].
*
* @param srcMesh 3-dimensional source mesh
* @param targetMesh 3-dimesional target mesh, containing only tetraedra
- * @param matrix vector of maps in which the result is stored
+ * @param result matrix in which the result is stored
*
*/
template<int SPACEDIM, int MESHDIM, class ConnType, NumberingPolicy numPol, class MatrixType, class MyMeshType>
}
}
+/*! Computes the triple product (XA^XB).XC (in 3D)*/
+inline double triple_product(const double* A, const double*B, const double*C, const double*X)
+{
+ double XA[3];
+ XA[0]=A[0]-X[0];
+ XA[1]=A[1]-X[1];
+ XA[2]=A[2]-X[2];
+ double XB[3];
+ XB[0]=B[0]-X[0];
+ XB[1]=B[1]-X[1];
+ XB[2]=B[2]-X[2];
+ double XC[3];
+ XC[0]=C[0]-X[0];
+ XC[1]=C[1]-X[1];
+ XC[2]=C[2]-X[2];
+
+ return
+ (XA[1]*XB[2]-XA[2]*XB[1])*XC[0]+
+ (XA[2]*XB[0]-XA[0]*XB[2])*XC[1]+
+ (XA[0]*XB[1]-XA[1]*XB[0])*XC[2];
+}
#endif
#include <map>
#include <ext/hash_map>
+#include "MEDMEM_define.hxx"
+#include "MEDMEM_CellModel.hxx"
+
using __gnu_cxx::hash_map;
namespace INTERP_KERNEL
inline void checkIsOutside(const double* pt, bool* isOutside) const;
inline void calculateNode(ConnType globalNodeNum);
inline void calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key);
+
/// disallow copying
IntersectorTetra(const IntersectorTetra& t);
/// reference to the source mesh
const NormalizedUnstructuredMesh<SPACEDIM,MESHDIM,ConnType,numPol,MyMeshType>& _srcMesh;
-
+
+
};
/**
#include "MeshUtils.hxx"
#include "VectorUtils.hxx"
#include "Log.hxx"
-// tony for the moment
-#include "MEDMEM_CellModel.hxx"
-#include "MEDMEM_define.hxx"
#include <cmath>
#include <cassert>
delete _t;
for(hash_map< int, double* >::iterator iter = _nodes.begin(); iter != _nodes.end() ; ++iter)
delete[] iter->second;
+
}
/**
////const medGeometryElement type = _srcMesh.getNumberOfNodesOfElement(element);
// get cell model for the element
- const CELLMODEL cellModel(type);
+ const CELLMODEL& cellModel= CELLMODEL_Map::retrieveCellModel(type);
// halfspace filtering
bool isOutside[8] = {true, true, true, true, true, true, true, true};
for(int i = 1 ; i <= cellModel.getNumberOfConstituents(1) ; ++i)
{
const medGeometryElement faceType = cellModel.getConstituentType(1, i);
- const CELLMODEL faceModel(faceType);
+ const CELLMODEL& faceModel= CELLMODEL_Map::retrieveCellModel(faceType);
assert(faceModel.getDimension() == 2);
Interpolation3D.hxx\
Interpolation2D.hxx\
Interpolation3DSurf.hxx\
-Remapper.hxx
-
+Remapper.hxx\
+PointLocator.hxx
# Libraries targets
TransformedTriangle_math.cxx\
BoundingBox.cxx\
TetraAffineTransform.cxx\
-Remapper.cxx
+Remapper.cxx \
+PointLocator.cxx
libinterpkernel_la_LDFLAGS= $(MED2_LIBS) $(HDF5_LIBS) ../MEDWrapper/V2_1/Core/libmed_V2_1.la $(STDLIB) ../MEDMEM/libmedmem.la -lutil -lm
ALL_FORTRAN_MODE
} NumberingPolicy;
+ class GenericMesh
+ {};
+
template<int SPACEDIM, int MESHDIM, class ConnType, NumberingPolicy numPol, class MyMeshType>
- class NormalizedUnstructuredMesh
+ class NormalizedUnstructuredMesh : public GenericMesh
{
public:
void getBoundingBox(double *boundingBox) const { asLeaf().getBoundingBox(boundingBox); }
--- /dev/null
+#include <list>
+#include "MEDMEM_Mesh.hxx"
+#include "MEDMEM_Exception.hxx"
+#include "MEDNormalizedUnstructuredMesh.txx"
+#include "PointLocatorAlgos.txx"
+#include "PointLocator.hxx"
+
+namespace INTERP_KERNEL {
+PointLocator::PointLocator(const MEDMEM::MESH& mesh)
+{
+ int meshdim=mesh.getMeshDimension();
+ int spacedim=mesh.getSpaceDimension();
+ if (meshdim != spacedim) throw MEDMEM::MEDEXCEPTION("Locator is not implemented for meshdim != spacedim");
+ switch (meshdim)
+ {
+ case 2:
+ _medmesh = new MEDNormalizedUnstructuredMesh<2,2> (&mesh);
+ _point_locator=new PointLocatorAlgos<2,2,int,ALL_FORTRAN_MODE,MEDNormalizedUnstructuredMesh<2,2> >(*(static_cast<MEDNormalizedUnstructuredMesh<2,2>* >(_medmesh)));
+ break;
+ case 3:
+ _medmesh = new MEDNormalizedUnstructuredMesh<3,3> (&mesh);
+ _point_locator=new PointLocatorAlgos<3,3,int,ALL_FORTRAN_MODE,MEDNormalizedUnstructuredMesh<3,3> >(*(static_cast<MEDNormalizedUnstructuredMesh<3,3>* >(_medmesh)));
+ break;
+ }
+}
+PointLocator::~PointLocator()
+{
+ delete _medmesh;
+ delete _point_locator;
+}
+
+ std::list<int> PointLocator::locate(const double*x)
+{
+ return _point_locator->locates(x);
+}
+}
--- /dev/null
+#ifndef _POINT_LOCATOR_HXX_
+#define _POINT_LOCATOR_HXX_
+#include <list>
+#include "NormalizedUnstructuredMesh.hxx"
+#include "MEDNormalizedUnstructuredMesh.hxx"
+#include "PointLocatorAlgos.txx"
+namespace INTERP_KERNEL
+{
+class PointLocator
+{
+public:
+ PointLocator(const MEDMEM::MESH& mesh);
+ virtual ~PointLocator();
+ std::list<int> locate(const double*x);
+
+private:
+ //MEDNormalizedUnstructuredMesh<2,2>* _medmesh;
+ //GenericPointLocatorAlgos*< 2,2,int,ALL_FORTRAN_MODE,MEDNormalizedUnstructuredMesh<2,2> >* _point_locator;
+
+ GenericMesh* _medmesh;
+ GenericPointLocatorAlgos* _point_locator;
+};
+}
+#endif
--- /dev/null
+#ifndef _POINT_LOCATOR_ALGOS_TXX_
+#define _POINT_LOCATOR_ALGOS_TXX_
+
+#include "MEDMEM_Exception.hxx"
+#include "InterpolationUtils.hxx"
+#include "MEDMEM_CellModel.hxx"
+#include "NormalizedUnstructuredMesh.hxx"
+#include "BBTree.txx"
+#include <list>
+namespace INTERP_KERNEL{
+
+ class GenericPointLocatorAlgos
+ {
+ public:
+ virtual ~GenericPointLocatorAlgos(){};
+ virtual std::list<int> locates(const double* x)=0;
+
+ };
+
+ template<int SPACEDIM, int MESHDIM, class ConnType, NumberingPolicy numPol, class MyMeshType> class PointLocatorAlgos: public GenericPointLocatorAlgos
+ {
+ private :
+ double* _bb;
+ BBTree<SPACEDIM>* _tree;
+ const NormalizedUnstructuredMesh<SPACEDIM,MESHDIM,ConnType,numPol,MyMeshType>& _mesh;
+
+ public:
+ PointLocatorAlgos<SPACEDIM,MESHDIM,ConnType,numPol,MyMeshType>(const NormalizedUnstructuredMesh<SPACEDIM,MESHDIM,ConnType,numPol,MyMeshType>& mesh):_mesh(mesh)
+ {
+ int nelem = _mesh.getNumberOfElements();
+ _bb = new double[SPACEDIM*2*nelem];
+ const ConnType* conn = _mesh.getConnectivityPtr();
+ const ConnType* conn_index = _mesh.getConnectivityIndexPtr();
+ const double* coords=_mesh.getCoordinatesPtr();
+ for (int i=0; i<nelem; i++)
+ {
+ for (int idim=0; idim<SPACEDIM; idim++)
+ {
+ _bb[2*(i*SPACEDIM+idim)]=HUGE;
+ _bb[2*(i*SPACEDIM+idim)+1]=-HUGE;
+ }
+ for (int index= conn_index[i]; index < conn_index[i+1];index++)
+ {
+ const double* coordelem = coords+OTT<ConnType,numPol>::ind2C(conn[OTT<ConnType,numPol>::ind2C(index)]);
+ for (int idim=0; idim<SPACEDIM;idim++)
+ {
+ _bb[2*(i*SPACEDIM+idim)]=(coordelem[idim]<_bb[2*(i*SPACEDIM+idim)])?coordelem[idim]:_bb[2*(i*SPACEDIM+idim)];
+ _bb[2*(i*SPACEDIM+idim)+1]=(coordelem[idim]>_bb[2*(i*SPACEDIM+idim)+1])?coordelem[idim]:_bb[2*(i*SPACEDIM+idim)+1];
+ }
+ }
+ }
+
+ _tree=new BBTree<SPACEDIM>(_bb,0,0,nelem);
+ }
+ ~PointLocatorAlgos<SPACEDIM,MESHDIM,ConnType,numPol,MyMeshType>()
+ {
+ delete[] _bb;
+ delete _tree;
+ }
+
+ std::list<int> locates(const double* x)
+ {
+ vector<int> candidates;
+ _tree->getElementsAroundPoint(x,candidates);
+ list<int> retlist;
+ for (int i=0; i< candidates.size(); i++)
+ {
+ if (elementContainsPoint(i,x))
+ retlist.push_back(OTT<ConnType,numPol>::indFC(i));
+ }
+ return retlist;
+ }
+
+ bool elementContainsPoint(int i, const double*x)
+ {
+ //as i is extracted from the BBTRee, it is already in C numbering
+ //it is not necessary to convert it from F to C
+
+ const double* coords= _mesh.getCoordinatesPtr();
+ const int* conn=_mesh.getConnectivityPtr();
+ const int* conn_index= _mesh.getConnectivityIndexPtr();
+ const int* conn_elem=conn+OTT<ConnType,numPol>::ind2C(conn_index[i]);
+
+
+ int nbnodes = conn_index[i+1]-conn_index[i];
+
+ // with dimension 2, it suffices to check all the edges
+ // and see if the sign of double products from the point
+ //is always the same.
+ // C
+ // / \
+ // / \
+ // Xo / \
+ // A-------B
+ //
+ //here XA^XC and XC^XB have different signs
+ //
+ if (SPACEDIM==2)
+ {
+ //in 2D, nbedges==nbnodes
+ int nbedges=nbnodes;
+ int sign[nbedges];
+ for (int iedge=0; iedge<nbedges; iedge++)
+ {
+ const double* A=coords+SPACEDIM*(OTT<ConnType,numPol>::ind2C(conn_elem[iedge]));
+ const double* B;
+ if (iedge+1< nbedges)
+ B=coords+SPACEDIM*(OTT<ConnType,numPol>::ind2C(conn_elem[iedge+1]));
+ else
+ B=coords+SPACEDIM*(OTT<ConnType,numPol>::ind2C(conn_elem[0]));
+
+ double a=mon_determinant(x, A, B);
+ if (a<-1e-12)
+ sign[iedge]=-1;
+ else if (a>1e-12)
+ sign[iedge]=1;
+ else
+ sign[iedge]=0;
+ }
+ return decide_from_sign(sign, nbedges);
+ }
+
+ if (SPACEDIM==3)
+ {
+ MED_EN::medGeometryElement elem_type;
+ switch (nbnodes) {
+ case 4 :
+ elem_type=MED_EN::MED_TETRA4;
+ break;
+ case 5:
+ elem_type=MED_EN::MED_PYRA5;
+ break;
+ case 6 :
+ elem_type=MED_EN::MED_PENTA6;
+ break;
+ case 8:
+ elem_type=MED_EN::MED_HEXA8;
+ break;
+ default:
+ throw MEDMEM::MEDEXCEPTION("PointLocatorAlgos : bad number of nodes in 3D locator");
+ }
+ const MEDMEM::CELLMODEL& model=MEDMEM::CELLMODEL_Map::retrieveCellModel(elem_type);
+ int nbfaces = model.getNumberOfConstituents(1);
+ int sign[nbfaces];
+ for (int iface=0; iface<nbfaces; iface++)
+ {
+ int* connface=model.getNodesConstituent(1,iface+1);
+ const double* AA=coords+SPACEDIM*(OTT<ConnType,numPol>::ind2C(conn_elem[connface[0]-1]));
+ const double* BB=coords+SPACEDIM*(OTT<ConnType,numPol>::ind2C(conn_elem[connface[1]-1]));
+ const double* CC=coords+SPACEDIM*(OTT<ConnType,numPol>::ind2C(conn_elem[connface[2]-1]));
+
+ double Vol=triple_product(AA,BB,CC,x);
+ if (Vol<-1e-12)
+ sign[iface]=-1;
+ else if (Vol>1e-12)
+ sign[iface]=1;
+ else
+ sign[iface]=0;
+ }
+ return decide_from_sign(sign, nbfaces);
+ }
+
+ }
+
+ bool decide_from_sign (const int* sign, int nbelem)
+ {
+ int min_sign =1;
+ int max_sign =-1;
+ for (int i=0; i<nbelem;i++)
+ {
+ min_sign=(sign[i]<min_sign)?sign[i]:min_sign;
+ max_sign=(sign[i]>max_sign)?sign[i]:max_sign;
+ }
+ if (min_sign==-1 && max_sign==1)
+ return false;
+ else
+ return true;
+ }
+ };
+}
+#endif
double bbox3[4]={5.0,6.0,7.0,9.0};
tree.getIntersectingElems(bbox3,elems);
CPPUNIT_ASSERT_EQUAL(2,(int)elems.size());
-
+ elems.clear();
+
+ double xx[2]={1.0,1.0};
+ tree.getElementsAroundPoint(xx,elems);
+ CPPUNIT_ASSERT_EQUAL(4,(int)elems.size());
+
delete[] bbox;
}
--- /dev/null
+#include "MEDMEM_Mesh.hxx"
+#include "MEDMEM_Meshing.hxx"
+
+MEDMEM::MESH* MEDMeshMaker(int dim, int nbedge, MED_EN::medGeometryElement type)
+{
+ MEDMEM::MESHING* mesh=new MEDMEM::MESHING();
+ mesh->setSpaceDimension(dim);
+ int nbnodes;
+ int nbelems;
+ switch (dim)
+ {
+ case 2:
+ nbnodes=(nbedge+1)*(nbedge+1);
+ if(type==MED_EN::MED_QUAD4)
+ nbelems=(nbedge*nbedge);
+ else
+ throw MEDMEM::MEDEXCEPTION("MEDMeshMaker: type not impletmented");
+ break;
+ case 3:
+ nbnodes=(nbedge+1)*(nbedge+1)*(nbedge+1);
+ if (type==MED_EN::MED_HEXA8)
+ nbelems= nbedge*nbedge*nbedge;
+ else
+ throw MEDMEM::MEDEXCEPTION("MEDMeshMaker: type not impletmented");
+ break;
+ }
+ double* coords = new double[dim*nbnodes];
+ int nz;
+ if (dim==2) nz =1; else nz=nbedge+1;
+ {
+ for (int ix=0; ix < nbedge+1; ix++)
+ for (int iy=0; iy<nbedge+1; iy++)
+ for (int iz=0; iz<nz;iz++)
+ {
+ int inode=(ix*(nbedge+1)*nz+iy*nz+iz);
+ coords[inode*dim]=double(ix)/double(nbedge);
+ coords[inode*dim+1]=double(iy)/double(nbedge);
+ if (dim==3)
+ coords[inode*dim+2]=double(iz)/double(nbedge);
+ }
+ }
+ mesh->setCoordinates(dim, nbnodes,coords,"CARTESIAN",MED_EN::MED_FULL_INTERLACE);
+
+ mesh->setNumberOfTypes(1,MED_EN::MED_CELL);
+ mesh->setTypes(&type,MED_EN::MED_CELL);
+ mesh->setNumberOfElements(&nbelems,MED_EN::MED_CELL);
+
+ int* conn = new int [nbelems*(type%100)];
+ if (dim==2)
+ {
+ for (int ix=0; ix<nbedge; ix++)
+ for (int iy=0; iy<nbedge; iy++)
+ {
+ int ielem=(ix*nbedge+iy);
+ conn [ielem*4]=ix*(nbedge+1)+iy+1;
+ conn [ielem*4+1]=ix*(nbedge+1)+iy+1+1;
+ conn [ielem*4+2]=(ix+1)*(nbedge+1)+iy+1+1;
+ conn [ielem*4+3]=(ix+1)*(nbedge+1)+iy+1;
+ }
+ }
+ if (dim==3)
+ {
+ for (int ix=0; ix<nbedge; ix++)
+ for (int iy=0; iy<nbedge; iy++)
+ for (int iz=0; iz<nbedge; iz++)
+ {
+ int ielem=(ix*nbedge*nbedge+iy*nbedge+iz);
+ conn [ielem*8]=ix*(nbedge+1)*(nbedge+1)+iy*(nbedge+1)+iz+1;
+ conn [ielem*8+1]=(ix+1)*(nbedge+1)*(nbedge+1)+iy*(nbedge+1)+iz+1;
+ conn [ielem*8+2]=(ix+1)*(nbedge+1)*(nbedge+1)+(iy+1)*(nbedge+1)+iz+1;
+ conn [ielem*8+3]=ix*(nbedge+1)*(nbedge+1)+(iy+1)*(nbedge+1)+iz+1;
+ conn [ielem*8+4]=ix*(nbedge+1)*(nbedge+1)+iy*(nbedge+1)+iz+1+1;
+ conn [ielem*8+5]=(ix+1)*(nbedge+1)*(nbedge+1)+iy*(nbedge+1)+iz+1+1;
+ conn [ielem*8+6]=(ix+1)*(nbedge+1)*(nbedge+1)+(iy+1)*(nbedge+1)+iz+1+1;
+ conn [ielem*8+7]=ix*(nbedge+1)*(nbedge+1)+(iy+1)*(nbedge+1)+iz+1+1;
+ }
+ }
+ mesh->setConnectivity(conn, MED_EN::MED_CELL,type);
+ mesh->setMeshDimension(dim);
+ return mesh;
+}
--- /dev/null
+#include "MEDMEM_define.hxx"
+namespace MEDMEM
+{
+ class MESH;
+}
+
+MEDMEM::MESH* MEDMeshMaker(int dim, int nbedge, MED_EN::medGeometryElement type);
TransformedTriangleIntersectTest.cxx \
BBTreeTest.cxx \
RemapperTest.cxx \
- SingleElementPlanarTests.cxx
+ SingleElementPlanarTests.cxx \
+ PointLocatorTest.cxx \
+ MEDMeshMaker.cxx
libInterpKernelTest_la_CPPFLAGS= @CPPUNIT_INCLUDES@ $(MED2_INCLUDES) $(HDF5_INCLUDES) \
-I$(srcdir)/.. -I$(srcdir)/../../MEDWrapper/V2_1/Core -I$(srcdir)/../../MEDMEM -DOPTIMIZE -DLOG_LEVEL=0
endif
AM_CPPFLAGS= $(libInterpKernelTest_la_CPPFLAGS)
+check-local:TestInterpKernel
+ ./TestInterpKernel
--- /dev/null
+#include "PointLocatorTest.hxx"
+#include "PointLocator.hxx"
+#include "MEDMeshMaker.hxx"
+
+#include <iostream>
+#include <list>
+
+namespace INTERP_TEST
+{
+
+
+ void PointLocatorTest::setUp()
+ {
+ }
+
+
+ void PointLocatorTest::tearDown()
+ {
+ }
+
+ /**
+ * Test that creates a tree in 2D and check that
+ * the results are correct in three
+ * cases :
+ * a non matching search
+ * a standard case
+ * a bbox overlapping the bboxes of the tree
+ */
+ void PointLocatorTest::test_PointLocator() {
+ MEDMEM::MESH* mesh2D= MEDMeshMaker(2,2,MED_EN::MED_QUAD4);
+ INTERP_KERNEL::PointLocator pl(*mesh2D) ;
+ double x[2]={0.0,0.0};
+ std::list<int> elems = pl.locate(x);
+ CPPUNIT_ASSERT_EQUAL(1,(int)elems.size());
+ CPPUNIT_ASSERT_EQUAL(1,(int)(*(elems.begin())));
+ elems.clear();
+
+ double x2[2]={0.25,0.25};
+ elems = pl.locate(x2);
+ CPPUNIT_ASSERT_EQUAL(1,(int)elems.size());
+ CPPUNIT_ASSERT_EQUAL(1,(int)(*(elems.begin())));
+ elems.clear();
+
+ double x3[2]={0.5,0.5};
+ elems = pl.locate(x3);
+ CPPUNIT_ASSERT_EQUAL(4,(int)elems.size());
+ elems.clear();
+
+ double x4[2]={-1.0,0.0};
+ elems = pl.locate(x4);
+ CPPUNIT_ASSERT_EQUAL(0,(int)elems.size());
+ elems.clear();
+
+ MEDMEM::MESH* mesh3D= MEDMeshMaker(3,2,MED_EN::MED_HEXA8);
+ INTERP_KERNEL::PointLocator pl3(*mesh3D);
+ double xx[3]={0.0,0.0,0.0};
+ elems = pl3.locate(xx);
+ CPPUNIT_ASSERT_EQUAL(1,(int)elems.size());
+ CPPUNIT_ASSERT_EQUAL(1,(int)*(elems.begin()));
+ elems.clear();
+
+ double xx2[3]={0.25,0.25,0.25};
+ elems = pl3.locate(xx2);
+ CPPUNIT_ASSERT_EQUAL(1,(int)elems.size());
+ CPPUNIT_ASSERT_EQUAL(1,(int)*(elems.begin()));
+ elems.clear();
+
+ double xx3[3]={0.5,0.5,0.5};
+ elems = pl3.locate(xx3);
+ CPPUNIT_ASSERT_EQUAL(8,(int)elems.size());
+ elems.clear();
+
+ double xx4[3]={-1.0,0.0,0.0};
+ elems = pl3.locate(x4);
+ CPPUNIT_ASSERT_EQUAL(0,(int)elems.size());
+ elems.clear();
+
+
+ }
+
+
+}
--- /dev/null
+#ifndef __TU_POINTLOCATOR_HXX__
+#define __TU_POINTLOCATOR_HXX__
+
+#include <cppunit/extensions/HelperMacros.h>
+#include "../PointLocator.hxx"
+
+namespace INTERP_TEST
+{
+
+ /**
+ * \brief Test suite testing some of the low level methods of TransformedTriangle.
+ *
+ */
+ class PointLocatorTest : public CppUnit::TestFixture
+ {
+
+ CPPUNIT_TEST_SUITE( PointLocatorTest );
+ CPPUNIT_TEST( test_PointLocator );
+ CPPUNIT_TEST_SUITE_END();
+
+
+ public:
+ void setUp();
+
+ void tearDown();
+
+ // tests
+ void test_PointLocator();
+
+ };
+
+
+
+
+}
+
+
+
+#endif
//target square is in reverse order as compared to initial square
- double source_integral=source_field.normL2(1,&source_areas);
+ double source_integral=source_field.normL2(1,&inv_source_areas);
double target_integral=target_field.normL2(1,&inv_target_areas);
CPPUNIT_ASSERT_DOUBLES_EQUAL(source_integral,target_integral,1e-10);
#include "SingleElementTetraTests.hxx"
#include "HexaTests.hxx"
#include "BBTreeTest.hxx"
+#include "PointLocatorTest.hxx"
#include "RemapperTest.hxx"
#include "MultiElement2DTests.hxx"
#include "SingleElementPlanarTests.hxx"
CPPUNIT_TEST_SUITE_REGISTRATION( INTERP_TEST::TransformedTriangleTest );
CPPUNIT_TEST_SUITE_REGISTRATION( BBTreeTest);
CPPUNIT_TEST_SUITE_REGISTRATION( RemapperTest);
+CPPUNIT_TEST_SUITE_REGISTRATION( PointLocatorTest);
CPPUNIT_TEST_SUITE_REGISTRATION( MultiElement2DTests );
CPPUNIT_TEST_SUITE_REGISTRATION( SingleElementPlanarTests );
