--- /dev/null
+#include "IntersectorTetra.hxx"
+
+#include "TetraAffineTransform.hxx"
+#include "TransformedTriangle.hxx"
+#include "MeshUtils.hxx"
+#include "VectorUtils.hxx"
+#include "Log.hxx"
+
+#include <cmath>
+#include <cassert>
+#include <string>
+#include <sstream>
+
+using namespace MEDMEM;
+using namespace MED_EN;
+using namespace INTERP_UTILS;
+
+namespace MEDMEM
+{
+ /*
+ * Constructor
+ *
+ * @param srcMesh mesh containing the source elements
+ * @param targetMesh mesh containing the target elements
+ *
+ */
+ IntersectorTetra::IntersectorTetra(const MESH& srcMesh, const MESH& targetMesh, int targetCell)
+ : _srcMesh(srcMesh), filtered(0)
+ {
+ const medGeometryElement targetType = targetMesh.getElementType(MED_CELL, targetCell);
+
+ // maybe we should do something more civilized here
+ assert(targetType == MED_TETRA4);
+
+ // get array of points of target tetraeder
+ const double* tetraCorners[4];
+ for(int i = 0 ; i < 4 ; ++i)
+ {
+ tetraCorners[i] = getCoordsOfNode(i + 1, targetCell, targetMesh);
+ }
+
+ // create AffineTransform from tetrahedron
+ _t = new TetraAffineTransform( tetraCorners );
+
+ }
+
+ /*
+ * Destructor
+ *
+ */
+ IntersectorTetra::~IntersectorTetra()
+ {
+ delete _t;
+
+ for(hash_map< int, double*>::iterator iter = _nodes.begin(); iter != _nodes.end() ; ++iter)
+ {
+ delete[] iter->second;
+ }
+ }
+
+ /*
+ * Calculates the volume of intersection of an element in the source mesh and an element
+ * in the target mesh. The method is based on the algorithm of Grandy. It first calculates the transformation
+ * that takes the target tetrahedron into the unit tetrahedron. After that, the
+ * faces of the source element are triangulated and the calculated transformation is applied
+ * to each triangle. The algorithm of Grandy, implemented in INTERP_UTILS::TransformedTriangle is used
+ * to calculate the contribution to the volume from each triangle. The volume returned is the sum of these contributions
+ * divided by the determinant of the transformation.
+ *
+ * @pre The element in _targetMesh referenced by targetCell is of type MED_TETRA4.
+ * @param srcCell global number of the source element (1 <= srcCell < # source cells)
+ * @param targetCell global number of the target element (1 <= targetCell < # target cells) - this element must be a tetrahedron
+ */
+ double IntersectorTetra::intersectSourceCell(int element)
+ {
+
+ //{ could be done on outside
+ // check if we have planar tetra element
+ if(_t->determinant() == 0.0)
+ {
+ // tetra is planar
+ LOG(2, "Planar tetra -- volume 0");
+ return 0.0;
+ }
+
+ // get type of cell
+ const medGeometryElement type = _srcMesh.getElementType(MED_CELL, element);
+
+ // get cell model for the element
+ const CELLMODEL cellModel(type);
+
+ assert(cellModel.getDimension() == 3);
+ assert(element >= 1);
+ assert(element <= _srcMesh.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS));
+
+ // halfspace filtering
+ bool isOutside[8] = {true, true, true, true, true, true, true, true};
+ bool isTargetOutside = false;
+
+ // calculate the coordinates of the nodes
+ for(int i = 1; i <= cellModel.getNumberOfNodes() ; ++i)
+ {
+ // we could store mapping local -> global numbers too, but not sure it is worth it
+ const int globalNodeNum = getGlobalNumberOfNode(i, element, _srcMesh);
+ if(_nodes.find(globalNodeNum) == _nodes.end())
+ {
+ calculateNode(globalNodeNum);
+ }
+
+ checkIsOutside(_nodes[globalNodeNum], isOutside);
+
+ // local caching of globalNodeNum
+ // not sure this is efficient
+ // globalNodeNumbers[i] = globalNodeNum;
+ }
+
+ // halfspace filtering check
+ // NB : might not be beneficial for caching of triangles
+ for(int i = 0; i < 8; ++i)
+ {
+ if(isOutside[i])
+ {
+ isTargetOutside = true;
+ }
+ }
+
+ double totalVolume = 0.0;
+
+ if(!isTargetOutside)
+ {
+ for(int i = 1 ; i <= cellModel.getNumberOfConstituents(1) ; ++i)
+ {
+ const medGeometryElement faceType = cellModel.getConstituentType(1, i);
+ const CELLMODEL faceModel(faceType);
+
+ assert(faceModel.getDimension() == 2);
+
+ int faceNodes[faceModel.getNumberOfNodes()];
+
+ // get the nodes of the face
+ for(int j = 1; j <= faceModel.getNumberOfNodes(); ++j)
+ {
+ const int locNodeNum = cellModel.getNodeConstituent(1, i, j);
+ assert(locNodeNum >= 1);
+ assert(locNodeNum <= cellModel.getNumberOfNodes());
+
+ faceNodes[j-1] = getGlobalNumberOfNode(locNodeNum, element, _srcMesh);//globalNodeNumbers[locNodeNum];
+ }
+
+ switch(faceType)
+ {
+ case MED_TRIA3:
+ {
+ // create the face key
+ TriangleFaceKey key = TriangleFaceKey(faceNodes[0], faceNodes[1], faceNodes[2]);
+
+ // calculate the triangle if needed
+ if(_volumes.find(key) == _volumes.end())
+ {
+ TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[1]], _nodes[faceNodes[2]]);
+ calculateVolume(tri, key);
+ totalVolume += _volumes[key];
+ } else {
+ // count negative as face has reversed orientation
+ totalVolume -= _volumes[key];
+ }
+ }
+ break;
+
+ case MED_QUAD4:
+
+ // simple triangulation of faces along a diagonal :
+ //
+ // 2 ------ 3
+ // | / |
+ // | / |
+ // | / |
+ // | / |
+ // | / |
+ // | / |
+ // 1 ------ 4
+ //
+ //? not sure if this always works
+ {
+ // calculate the triangles if needed
+
+ // local nodes 1, 2, 3
+ TriangleFaceKey key1 = TriangleFaceKey(faceNodes[0], faceNodes[1], faceNodes[2]);
+ if(_volumes.find(key1) == _volumes.end())
+ {
+ TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[1]], _nodes[faceNodes[2]]);
+ calculateVolume(tri, key1);
+ totalVolume += _volumes[key1];
+ } else {
+ // count negative as face has reversed orientation
+ totalVolume -= _volumes[key1];
+ }
+
+ // local nodes 1, 3, 4
+ TriangleFaceKey key2 = TriangleFaceKey(faceNodes[0], faceNodes[2], faceNodes[3]);
+ if(_volumes.find(key2) == _volumes.end())
+ {
+ TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[2]], _nodes[faceNodes[3]]);
+ calculateVolume(tri, key2);
+ totalVolume += _volumes[key2];
+ }
+ else
+ {
+ // count negative as face has reversed orientation
+ totalVolume -= _volumes[key2];
+ }
+ }
+ break;
+
+ default:
+ std::cout << "+++ Error : Only elements with triangular and quadratilateral faces are supported at the moment." << std::endl;
+ assert(false);
+ }
+ }
+ }
+ else
+ {
+ ++filtered;
+ }
+
+ // reset if it is very small to keep the matrix sparse
+ // is this a good idea?
+ if(epsilonEqual(totalVolume, 0.0, 1.0e-11))
+ {
+ totalVolume = 0.0;
+ }
+
+ LOG(2, "Volume = " << totalVolume << ", det= " << _t->determinant());
+
+ // NB : fault in article, Grandy, [8] : it is the determinant of the inverse transformation
+ // that should be used (which is equivalent to dividing by the determinant)
+ return std::abs(1.0 / _t->determinant() * totalVolume) ;
+ }
+
+};
--- /dev/null
+#ifndef __INTERSECTOR_TETRA_HXX__
+#define __INTERSECTOR_TETRA_HXX__
+
+#include "MEDMEM_define.hxx"
+#include "MEDMEM_Mesh.hxx"
+
+#include "Intersector.hxx"
+#include <vector>
+#include <ext/hash_map>
+#include <functional>
+
+#include "TetraAffineTransform.hxx"
+#include "TransformedTriangle.hxx"
+
+using __gnu_cxx::hash_map;
+
+namespace INTERP_UTILS
+{
+ class TriangleFaceKey
+ {
+ public:
+ int _nodes[3];
+ int _hashVal;
+
+ TriangleFaceKey(int node1, int node2, int node3)
+ {
+ sort3Ints(_nodes, node1, node2, node3);
+ // assert(_nodes[0] < _nodes[1]);
+ //assert(_nodes[0] < _nodes[2]);
+ //assert(_nodes[1] < _nodes[2]);
+ _hashVal = ( _nodes[0] + _nodes[1] + _nodes[2] ) % 29;
+ }
+
+ bool operator==(const TriangleFaceKey& rhs) const
+ {
+ return _nodes[0] == rhs._nodes[0] && _nodes[1] == rhs._nodes[1] && _nodes[2] == rhs._nodes[2];
+ }
+
+ int hashVal() const
+ {
+ return _hashVal;
+ }
+
+ inline void sort3Ints(int* sorted, int node1, int node2, int node3);
+ };
+}
+
+namespace __gnu_cxx
+{
+ template<>
+ class hash<INTERP_UTILS::TriangleFaceKey>
+ {
+ public:
+ int operator()(const INTERP_UTILS::TriangleFaceKey& key) const
+ {
+ return key.hashVal();
+ }
+ };
+};
+
+namespace INTERP_UTILS
+{
+
+
+ /*
+ * Class calculating the volume of intersection of individual 3D elements.
+ *
+ */
+ class IntersectorTetra
+ {
+
+ public :
+ IntersectorTetra(const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh, int targetCell);
+
+ ~IntersectorTetra();
+
+ double intersectSourceCell(int srcCell);
+
+ mutable int filtered;
+
+ private:
+
+ TetraAffineTransform* _t;
+
+ hash_map< int, double*> _nodes;
+ hash_map< TriangleFaceKey, double > _volumes;
+
+ inline void checkIsOutside(const double* pt, bool* isOutside) const;
+ inline void calculateNode(int globalNodeNum);
+ inline void calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key);
+ // inline void sort3Ints(int* sorted, int node1, int node2, int node3);
+ //inline std::string createFaceKey(int node1, int node2, int node3);
+
+ const MEDMEM::MESH& _srcMesh;
+
+ };
+
+ inline void IntersectorTetra::checkIsOutside(const double* pt, bool* isOutside) const
+ {
+ isOutside[0] = isOutside[0] && (pt[0] <= 0.0);
+ isOutside[1] = isOutside[1] && (pt[0] >= 1.0);
+ isOutside[2] = isOutside[2] && (pt[1] <= 0.0);
+ isOutside[3] = isOutside[3] && (pt[1] >= 1.0);
+ isOutside[4] = isOutside[4] && (pt[2] <= 0.0);
+ isOutside[5] = isOutside[5] && (pt[2] >= 1.0);
+ isOutside[6] = isOutside[6] && (1.0 - pt[0] - pt[1] - pt[2] <= 0.0);
+ isOutside[7] = isOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] >= 1.0);
+ }
+
+ inline void IntersectorTetra::calculateNode(int globalNodeNum)
+ {
+ assert(globalNodeNum >= 0);
+ assert(globalNodeNum < _srcMesh.getNumberOfNodes());
+ const double* node = &(_srcMesh.getCoordinates(MED_EN::MED_FULL_INTERLACE)[3*globalNodeNum]);
+ double* transformedNode = new double[3];
+ assert(transformedNode != 0);
+
+ _t->apply(transformedNode, node);
+ _nodes[globalNodeNum] = transformedNode;
+ }
+
+ inline void IntersectorTetra::calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key)
+ {
+ const double vol = tri.calculateIntersectionVolume();
+ _volumes.insert(make_pair(key, vol));
+ }
+
+ inline void TriangleFaceKey::sort3Ints(int* sorted, int node1, int node2, int node3)
+ {
+ if(node1 < node2)
+ {
+ if(node1 < node3)
+ {
+ // node 1 is min
+ sorted[0] = node1;
+ sorted[1] = node2 < node3 ? node2 : node3;
+ sorted[2] = node2 < node3 ? node3 : node2;
+ }
+ else
+ {
+ // node3 , node1, node2
+ sorted[0] = node3;
+ sorted[1] = node1;
+ sorted[2] = node2;
+ }
+ }
+ else // x2 < x1
+ {
+ if(node2 < node3)
+ {
+ // node 2 is min
+ sorted[0] = node2;
+ sorted[1] = node1 < node3 ? node1 : node3;
+ sorted[2] = node1 < node3 ? node3 : node1;
+ }
+ else
+ {
+ // node 3, node 2, node 1
+ sorted[0] = node3;
+ sorted[1] = node2;
+ sorted[2] = node1;
+ }
+ }
+ }
+
+#if 0
+ inline std::string IntersectorTetra::createFaceKey(int node1, int node2, int node3)
+ {
+ int sorted[3];
+ sort3Ints(sorted, node1, node2, node3);
+
+ std::stringstream sstr;
+ sstr << node1 << "-" << node2 << "-" << node3;
+ return sstr.str();
+ }
+#endif
+
+};
+
+
+#endif
MeshUtils.hxx\
Intersector3D.hxx\
Log.hxx\
-TransformedTriangle_inline.hxx
+TransformedTriangle_inline.hxx\
+IntersectorTetra.hxx
# Libraries targets
Interpolation3D.cxx\
Intersector3D.cxx\
Interpolation3DSurf.cxx\
-Interpolation2D.cxx
+Interpolation2D.cxx\
+IntersectorTetra.cxx
# Executables targets
BIN =
CPPFLAGS+=$(BOOST_CPPFLAGS)
# optimization
-#CXXFLAGS+= -O1 -DOPTIMIZE
-#CPPFLAGS+= -O1 -DOPTIMIZE
+
CXXFLAGS+= -O2 -DOPTIMIZE_FILTER -DOPTIMIZE
CPPFLAGS+= -O2 -DOPTIMIZE_FILTER -DOPTIMIZE
#include "MEDMEM_define.hxx"
#include <cassert>
+#include <set>
+
using namespace MEDMEM;
using namespace MED_EN;
return static_cast<int>(type) - 300;
}
+ inline int getGlobalNumberOfNode(int node, int element, const MESH& mesh)
+ {
+ assert(node >= 1);
+ assert(node <= mesh.getNumberOfNodes());
+ const int nodeOffset = node - 1;
+ const int elemIdx = mesh.getConnectivityIndex(MED_NODAL, MED_CELL)[element - 1] - 1;
+ return mesh.getConnectivity(MED_FULL_INTERLACE, MED_NODAL, MED_CELL, MED_ALL_ELEMENTS)[elemIdx + nodeOffset] - 1;
+ }
+
};
m = interpolator->interpol_maillages(sMesh, tMesh);
- testVolumes(m, sMesh, tMesh);
+ // if reflexive, check volumes
+ if(strcmp(mesh1path,mesh2path) == 0)
+ {
+ const bool row_and_col_sums_ok = testVolumes(m, sMesh, tMesh);
+ CPPUNIT_ASSERT_EQUAL_MESSAGE("Row or column sums incorrect", true, row_and_col_sums_ok);
+ }
LOG(1, "Intersection calculation done. " << std::endl );
// multi - element
CPPUNIT_TEST( tetraComplexIncluded );
-#endif
+
CPPUNIT_TEST( dividedUnitTetraSimplerReflexive );
CPPUNIT_TEST( dividedUnitTetraReflexive );
- //#if 0
- //CPPUNIT_TEST( nudgedDividedUnitTetra );
- //CPPUNIT_TEST( nudgedDividedUnitTetraSimpler );
- //CPPUNIT_TEST( dividedGenTetra );
+
+ CPPUNIT_TEST( nudgedDividedUnitTetra );
+ CPPUNIT_TEST( nudgedDividedUnitTetraSimpler );
+ CPPUNIT_TEST( dividedGenTetra );
CPPUNIT_TEST( boxReflexive );
CPPUNIT_TEST( boxReflexiveModerate );
- //CPPUNIT_TEST( tetraBoxes );
- //CPPUNIT_TEST( moderateBoxes );
- //CPPUNIT_TEST( moderateBoxesSmaller );
+ CPPUNIT_TEST( tetraBoxes );
+ CPPUNIT_TEST( moderateBoxes );
+#endif
+ CPPUNIT_TEST( moderateBoxesSmaller );
CPPUNIT_TEST( moderateBoxSmallReflexive );
+#if 0
CPPUNIT_TEST( moderateBoxEvenSmallerReflexive );
CPPUNIT_TEST( tinyBoxReflexive );
- //CPPUNIT_TEST( simpleHexaBox );
- //#endif
+ CPPUNIT_TEST( simpleHexaBox );
+#endif
CPPUNIT_TEST_SUITE_END();
@COMMENCE@
# header files
-EXPORT_HEADERS = TestingUtils.hxx
+EXPORT_HEADERS = CppUnitTest.hxx \
+ TransformedTriangleTest.hxx \
+ TransformedTriangleIntersectTest.hxx \
+ Interpolation3DTest.hxx \
+ TestingUtils.hxx
# Libraries targets
# Executables targets
-BIN = PerfTest
+BIN = PerfTest TestInterpKernel
-BIN_SRC =
+BIN_SRC = CppUnitTest.cxx TransformedTriangleTest.cxx TransformedTriangleIntersectTest.cxx \
+Interpolation3DTest.cxx
BIN_CLIENT_IDL =
CPPFLAGS+=$(BOOST_CPPFLAGS) $(MED2_INCLUDES) $(HDF5_INCLUDES) -I$(top_srcdir)/src/INTERP_KERNEL
#include CppUnit
-#CXXFLAGS+= -I/usr/include/cppunit
-#CPPFLAGS+= -I/usr/include/cppunit
+CXXFLAGS+= -I/usr/include/cppunit
+CPPFLAGS+= -I/usr/include/cppunit
# for log
-CXXFLAGS+= -DLOG_LEVEL=0
-CPPFLAGS+= -DLOG_LEVEL=0
+CXXFLAGS+= -DLOG_LEVEL=3 #-DOPTIMIZE -O2
+CPPFLAGS+= -DLOG_LEVEL=3 #-DOPTIMIZE -O2
# for gcov
#CXXFLAGS+=-fprofile-arcs -ftest-coverage
#CPPFLAGS+=-fprofile-arcs -ftest-coverage
#for gprof
-CXXFLAGS+=-pg
-CPPFLAGS+=-pg
+#CXXFLAGS+=-pg
+#CPPFLAGS+=-pg
#LDFLAGS+=$(MED2_LIBS) $(HDF5_LIBS)
# change motivated by the bug KERNEL4778.
# for gcov
#LDFLAGS+= -lgcov
-LDFLAGS+= -pg
+#LDFLAGS+= -pg
#LDFLAGSFORBIN+=$(MED2_LIBS) $(HDF5_LIBS)
# change motivated by the bug KERNEL4778.
LIBS = @LIBS@ @CPPUNIT_LIBS@
-LDFLAGSFORBIN += $(LDFLAGS) -lm $(MED2_LIBS) $(HDF5_LIBS) \
- -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome -lmed_V2_1 -lmedmem \
- -linterpkernel
-#LDFLAGSFORBIN += $(LDFLAGS) -lm $(HDF5_LIBS) \
-# -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome \
-# -lcppunit -linterpkernel
+LDFLAGSFORBIN += $(LDFLAGS) -lm $(HDF5_LIBS) \
+ -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome \
+ -lcppunit -linterpkernel
-LDFLAGSFORBIN += -pg
+#LDFLAGSFORBIN += -pg
-UNIT_TEST_PROG = #PerfTest
+UNIT_TEST_PROG = PerfTest TestInterpKernel
@CONCLUDE@
calcIntersectionMatrix(mesh1path.c_str(), mesh1.c_str(), mesh2path.c_str(), mesh2.c_str(), m);
- //dumpIntersectionMatrix(m);
+ // dumpIntersectionMatrix(m);
return 0;
class TransformedTriangleTest;
class TransformedTriangleIntersectTest;
+
+
namespace INTERP_UTILS
{
};
+#include "TransformedTriangle_tables.hxx"
+
#endif
C_ZH, C_ZX, C_YZ, // OXY
C_XH, C_YH, C_ZH // XYZ
};
-#if 0
- template<TetraFacet facet, TriSegment seg>
- inline TransformedTriangle::DoubleProduct getDPForSegFacetIntersection()
- {
- return NO_DP;
- }
-
- template<>
- inline TransformedTriangle::DoubleProduct getDPForSegFacetIntersection<OYZ, PQ>
- {
- return C_XH;
- }
-
-#define DEF_DP_FOR_SEG_FACET(FACET,SEG,DP) template<> inline TransformedTriangle::DoubleProduct getDPForSegFacetIntersection<FACET,SEG> { return DP; }
- DEF_DP_FOR_SEG_FACET()
-
-#endif
// signs associated with entries in DP_FOR_SEGMENT_FACET_INTERSECTION
const double TransformedTriangle::SIGN_FOR_SEG_FACET_INTERSECTION[12] =
LOG(4, "tA = " << tA << " tB = " << tB << " alpha= " << alpha );
for(int i = 0; i < 3; ++i)
{
+
pt[i] = (1 - alpha) * COORDS_TET_CORNER[3*corners[0] + i] +
alpha * COORDS_TET_CORNER[3*corners[1] + i];
+#if 0
+ pt[i] = (1 - alpha) * getCoordinateForTetCorner<corners[0], i>() +
+ alpha * getCoordinateForTetCorner<corners[0], i>();
+#endif
LOG(6, pt[i] );
assert(pt[i] >= 0.0);
assert(pt[i] <= 1.0);
_isTripleProductsCalculated = true;
}
- /*
+ /**
* Calculates the angle between an edge of the tetrahedron and the triangle
*
* @param edge edge of the tetrahedron
