From bad6ca6e34a5d40365f378ba31e927e8b32486c5 Mon Sep 17 00:00:00 2001 From: vbd Date: Thu, 13 Sep 2007 09:16:02 +0000 Subject: [PATCH] staffan : * restructured test suites * added IntersectorHexa --- src/INTERP_KERNEL/IntersectorHexa.cxx | 293 +++++++++++++++ src/INTERP_KERNEL/IntersectorHexa.hxx | 82 +++++ src/INTERP_KERNEL/Makefile.in | 2 +- src/INTERP_KERNEL/Test/HexaTests.hxx | 42 +++ .../Test/Interpolation3DTestSuite.hxx | 31 ++ src/INTERP_KERNEL/Test/MeshTestToolkit.cxx | 341 ++++++++++++++++++ src/INTERP_KERNEL/Test/MeshTestToolkit.hxx | 54 +++ .../Test/MultiElementTetraTests.hxx | 96 +++++ .../Test/SingleElementTetraTests.hxx | 121 +++++++ .../Test/TransformedTriangleIntersectTest.cxx | 15 + 10 files changed, 1076 insertions(+), 1 deletion(-) create mode 100644 src/INTERP_KERNEL/IntersectorHexa.cxx create mode 100644 src/INTERP_KERNEL/IntersectorHexa.hxx create mode 100644 src/INTERP_KERNEL/Test/HexaTests.hxx create mode 100644 src/INTERP_KERNEL/Test/Interpolation3DTestSuite.hxx create mode 100644 src/INTERP_KERNEL/Test/MeshTestToolkit.cxx create mode 100644 src/INTERP_KERNEL/Test/MeshTestToolkit.hxx create mode 100644 src/INTERP_KERNEL/Test/MultiElementTetraTests.hxx create mode 100644 src/INTERP_KERNEL/Test/SingleElementTetraTests.hxx diff --git a/src/INTERP_KERNEL/IntersectorHexa.cxx b/src/INTERP_KERNEL/IntersectorHexa.cxx new file mode 100644 index 000000000..b80ea655a --- /dev/null +++ b/src/INTERP_KERNEL/IntersectorHexa.cxx @@ -0,0 +1,293 @@ +#include "IntersectorHexa.hxx" + +#include "MEDMEM_define.hxx" +#include "MEDMEM_Mesh.hxx" + +#include "MeshUtils.hxx" + +#include "IntersectorTetra.hxx" + +using namespace MED_EN; +using namespace MEDMEM; + +namespace INTERP_UTILS +{ + + IntersectorHexa::IntersectorHexa(const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh, int targetCell, SplittingPolicy policy) + { + + // check type + const medGeometryElement targetType = targetMesh.getElementType(MED_CELL, targetCell); + assert(targetType == MED_HEXA8); + + const int numTetra = static_cast(policy); + + // pre-calculate nodes + calculateSubNodes(targetMesh, targetCell, policy); + + _tetra.reserve(numTetra); + _nodes.reserve(30); // we never have more than this + + switch(policy) + { + case PLANAR_FACE_5: + { + const int subZone[8] = { 1, 2, 3, 4, 5, 6, 7, 8 }; + fiveSplit(srcMesh, subZone); + } + break; + + case PLANAR_FACE_6: + { + const int subZone[8] = { 1, 2, 3, 4, 5, 6, 7, 8 }; + sixSplit(srcMesh, subZone); + } + break; + + case GENERAL_24: + { + calculateGeneral24Tetra(srcMesh); + } + break; + + case GENERAL_48: + { + calculateGeneral48Tetra(srcMesh); + } + break; + default: + assert(false); + } + } + + IntersectorHexa::~IntersectorHexa() + { + for(vector::iterator iter = _tetra.begin(); iter != _tetra.end(); ++iter) + { + delete *iter; + } + + // free potential sub-mesh nodes that have been allocated + vector::iterator iter = _nodes.begin() + 8; + while(iter != _nodes.end()) + { + delete[] *iter; + ++iter; + } + } + + void IntersectorHexa::fiveSplit(const MEDMEM::MESH& srcMesh, const int* const subZone) + { + static const int SPLIT_NODES_5[20] = + { + 0, 1, 5, 2, + 0, 4, 5, 7, + 0, 3, 7, 2, + 5, 6, 7, 2, + 0, 2, 5, 7 + }; + + for(int i = 0; i < 5; ++i) + { + const double* nodes[4]; + for(int j = 0; j < 4; ++j) + { + nodes[j] = getCoordsOfSubNode(subZone[ SPLIT_NODES_5[4*i+j] ]); + } + IntersectorTetra* t = new IntersectorTetra(srcMesh, nodes); + _tetra.push_back(t); + } + } + + void IntersectorHexa::sixSplit(const MEDMEM::MESH& srcMesh, const int* const subZone) + { + static const int SPLIT_NODES_6[24] = + { + 0, 1, 5, 6, + 0, 2, 1, 6, + 0, 5, 4, 6, + 0, 4, 7, 6, + 0, 3, 2, 6, + 0, 7, 3, 6 + }; + + for(int i = 0; i < 6; ++i) + { + const double* nodes[4]; + for(int j = 0; j < 4; ++j) + { + nodes[j] = getCoordsOfSubNode(subZone[ SPLIT_NODES_6[4*i+j] ]); + } + IntersectorTetra* t = new IntersectorTetra(srcMesh, nodes); + _tetra.push_back(t); + } + } + + void IntersectorHexa::calculateGeneral24Tetra(const MEDMEM::MESH& srcMesh) + { + // the two mesh nodes used in each tetrahedron + // the tetrahedra all have nodes (cellCenter, faceCenter, edgeNode1, edgeNode2) + static const int TETRA_EDGES[48] = + { + // face with center 9 + 1, 2, + 2, 6, + 6, 5, + 5, 1, + // face with center 10 + 1, 2, + 2, 3, + 3, 4, + 4, 1, + // face with center 11 + 1, 5, + 5, 8, + 8, 4, + 4, 1, + // face with center 12 + 2, 6, + 6, 7, + 7, 3, + 3, 2, + // face with center 13 + 6, 7, + 7, 8, + 8, 5, + 5, 6, + // face with center 14 + 3, 7, + 7, 8, + 8, 4, + 4, 3 + }; + + const double* nodes[4]; + + // get the cell center + nodes[0] = getCoordsOfSubNode(15); + + for(int faceCenterNode = 9; faceCenterNode <= 14; ++faceCenterNode) + { + // get the face center + nodes[1] = getCoordsOfSubNode(faceCenterNode); + for(int j = 0; j < 4; ++j) + { + const int row = 4*(faceCenterNode - 9) + j; + nodes[2] = getCoordsOfSubNode(TETRA_EDGES[2*row]); + nodes[3] = getCoordsOfSubNode(TETRA_EDGES[2*row + 1]); + + IntersectorTetra* t = new IntersectorTetra(srcMesh, nodes); + _tetra.push_back(t); + } + } + } + + void IntersectorHexa::calculateGeneral48Tetra(const MEDMEM::MESH& srcMesh) + { + // define 8 hexahedral subzones as in Grandy, p449 + // the values correspond to the nodes that correspond to nodes 1,2,3,4,5,6,7,8 in the subcell + // these nodes have correspondance 1 -> 0, 2 -> a, 3 -> e, 4 -> d, 5 -> b, 6 -> c, 7 -> g, 8 -> f + // with the Fig. 4.c in Grandy + static const int subZones[64] = + { + 1, 9, 22, 13, 10, 21, 27, 23, + 9, 2, 14, 22, 21, 11, 24, 27, + 13, 22, 17, 4, 23, 27, 26, 18, + 22, 14, 3, 17, 27, 24, 19, 26, + 10, 21, 27, 23, 5, 12, 25, 15, + 21, 11, 24, 27, 12, 6, 16, 25, + 23, 27, 26, 18, 15, 25, 20, 8, + 27, 24, 19, 26, 25, 16, 7, 20 + }; + + for(int i = 0; i < 8; ++i) + { + sixSplit(srcMesh, &subZones[8*i]); + } + } + + void IntersectorHexa::calculateSubNodes(const MEDMEM::MESH& targetMesh, int targetCell, SplittingPolicy policy) + { + // retrieve real mesh nodes + for(int node = 1; node <= 8 ; ++node) + { + // calculate only normal nodes + _nodes.push_back(getCoordsOfNode(node, targetCell, targetMesh)); + } + + // create sub-mesh nodes if needed + switch(policy) + { + case GENERAL_24: + { + static const int GENERAL_24_SUB_NODES[28] = + { + 1, 2, 5, 6, // sub-node 9 (face) + 1, 2, 3, 4, // sub-node 10 (face) + 1, 4, 5, 8, // sub-node 11 (face) + 2, 3, 6, 7, // sub-node 12 (face) + 5, 6, 7, 8, // sub-node 13 (face) + 3, 4, 7, 8, // sub-node 14 (face) + 10, 11, 12, 13 // sub-node 15 (cell) + }; + + for(int i = 0; i < 7; ++i) + { + double* barycenter = new double[3]; + calcBarycenter<4>(barycenter, &GENERAL_24_SUB_NODES[4*i]); + _nodes.push_back(barycenter); + } + } + break; + + case GENERAL_48: + { + static const int GENERAL_48_SUB_NODES[38] = + { + 1, 2, // sub-node 9 (edge) + 1, 5, // sub-node 10 (edge) + 2, 6, // sub-node 11 (edge) + 5, 6, // sub-node 12 (edge) + 1, 4, // sub-node 13 (edge) + 2, 3, // sub-node 14 (edge) + 5, 8, // sub-node 15 (edge) + 6, 7, // sub-node 16 (edge) + 3, 4, // sub-node 17 (edge) + 4, 8, // sub-node 18 (edge) + 3, 7, // sub-node 19 (edge) + 7, 8, // sub-node 20 (edge) + 9, 12, // sub-node 21 (face) + 13, 14,// sub-node 22 (face) + 10, 18,// sub-node 23 (face) + 11, 19,// sub-node 24 (face) + 15, 16,// sub-node 25 (face) + 17, 20,// sub-node 26 (face) + 21, 26 // sub-node 27 (cell) + }; + + for(int i = 0; i < 19; ++i) + { + double* barycenter = new double[3]; + calcBarycenter<2>(barycenter, &GENERAL_48_SUB_NODES[2*i]); + _nodes.push_back(barycenter); + } + } + break; + + default: + break; + } + } + + double IntersectorHexa::intersectSourceCell(int srcCell) + { + double volume = 0.0; + for(vector::iterator iter = _tetra.begin(); iter != _tetra.end(); ++iter) + { + volume += (*iter)->intersectSourceCell(srcCell); + } + return volume; + } + + +}; diff --git a/src/INTERP_KERNEL/IntersectorHexa.hxx b/src/INTERP_KERNEL/IntersectorHexa.hxx new file mode 100644 index 000000000..47efdca41 --- /dev/null +++ b/src/INTERP_KERNEL/IntersectorHexa.hxx @@ -0,0 +1,82 @@ +#ifndef __INTERSECTOR_HEXA_HXX__ +#define __INTERSECTOR_HEXA_HXX__ + +#include "MEDMEM_define.hxx" +#include "MEDMEM_Mesh.hxx" +#include "TargetIntersector.hxx" + + +namespace INTERP_UTILS +{ + + class IntersectorTetra; + + /** + * Class representing a hexahedron, which allows + * + * + */ + class IntersectorHexa : public TargetIntersector + { + + public: + + enum SplittingPolicy { PLANAR_FACE_5 = 5, PLANAR_FACE_6 = 6, GENERAL_24 = 24, GENERAL_48 = 48 }; + + IntersectorHexa(const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh, int targetCell, SplittingPolicy policy = GENERAL_48); + + ~IntersectorHexa(); + + virtual double intersectSourceCell(int srcCell); + + private: + + void fiveSplit(const MEDMEM::MESH& srcMesh, const int* const subZone); + + void sixSplit(const MEDMEM::MESH& srcMesh, const int* const subZone); + + void calculateGeneral24Tetra(const MEDMEM::MESH& srcMesh); + + void calculateGeneral48Tetra(const MEDMEM::MESH& srcMesh); + + void calculateSubNodes(const MEDMEM::MESH& targetMesh, int targetCell, SplittingPolicy policy); + + inline const double* getCoordsOfSubNode(int node); + + template + inline void calcBarycenter(double* barycenter, const int* const pts); + + vector _tetra; + + vector _nodes; + + }; + + + inline const double* IntersectorHexa::getCoordsOfSubNode(int node) + { + // replace at with [] for unsafe but faster access + return _nodes.at(node - 1); + } + + template + inline void IntersectorHexa::calcBarycenter(double* barycenter, const int* const pts) + { + barycenter[0] = barycenter[1] = barycenter[2] = 0.0; + for(int i = 0; i < n ; ++i) + { + const double* pt = getCoordsOfSubNode(pts[i]); + barycenter[0] += pt[0]; + barycenter[1] += pt[1]; + barycenter[2] += pt[2]; + } + + barycenter[0] /= n; + barycenter[1] /= n; + barycenter[2] /= n; + } + +}; + + +#endif diff --git a/src/INTERP_KERNEL/Makefile.in b/src/INTERP_KERNEL/Makefile.in index 1e21f76d9..cadadcff9 100644 --- a/src/INTERP_KERNEL/Makefile.in +++ b/src/INTERP_KERNEL/Makefile.in @@ -84,7 +84,7 @@ BIN_SRC = BIN_SERVER_IDL = BIN_CLIENT_IDL = -TEST_PROGS = testUnitTetra test_a_mesh +TEST_PROGS = testUnitTetra LDFLAGS+= -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome LDFLAGSFORBIN+= -L$(top_builddir)/lib@LIB_LOCATION_SUFFIX@/salome diff --git a/src/INTERP_KERNEL/Test/HexaTests.hxx b/src/INTERP_KERNEL/Test/HexaTests.hxx new file mode 100644 index 000000000..be7f75e5e --- /dev/null +++ b/src/INTERP_KERNEL/Test/HexaTests.hxx @@ -0,0 +1,42 @@ +#ifndef __HEXA_TESTS_HXX_ +#define __HEXA_TESTS_HXX_ + +#include "Interpolation3DTestSuite.hxx" + +namespace INTERP_UTILS +{ + class HexaTests : public Interpolation3DTestSuite + { + CPPUNIT_TEST_SUITE( HexaTests ); + + CPPUNIT_TEST( simpleHexaBox ); + CPPUNIT_TEST( reflexiveHexaBox ); + CPPUNIT_TEST( hexaBoxes ); + CPPUNIT_TEST( hexaBoxesMoved ); + + CPPUNIT_TEST_SUITE_END(); + + public: + void simpleHexaBox() + { + _testTools->intersectMeshes("meshes/BoxHexa1.med", "BoxHexa1", "meshes/BoxTetra2.med", "BoxTetra2", 65250, 1.0e-5); + } + + void reflexiveHexaBox() + { + _testTools->intersectMeshes("meshes/BoxHexa1.med", "BoxHexa1", "meshes/BoxHexa1.med", "BoxHexa1", 204000); + } + + void hexaBoxes() + { + _testTools->intersectMeshes("meshes/BoxHexa1.med", "BoxHexa1", "meshes/BoxHexa2.med", "BoxHexa2", 65250); + } + + void hexaBoxesMoved() + { + _testTools->intersectMeshes("meshes/MovedHexaBox1.med", "MovedHexaBox1", "meshes/MovedHexaBox2.med", "MovedHexaBox2", 65250); + } + + }; +} +#endif diff --git a/src/INTERP_KERNEL/Test/Interpolation3DTestSuite.hxx b/src/INTERP_KERNEL/Test/Interpolation3DTestSuite.hxx new file mode 100644 index 000000000..71dd4e163 --- /dev/null +++ b/src/INTERP_KERNEL/Test/Interpolation3DTestSuite.hxx @@ -0,0 +1,31 @@ +#ifndef __TU_INTERPOLATION_3D_TEST_SUITE_HXX__ +#define __TU_INTERPOLATION_3D_TEST_SUITE_HXX__ + +#include "MeshTestToolkit.hxx" + +#include + +namespace INTERP_UTILS +{ + + class Interpolation3DTestSuite : public CppUnit::TestFixture + { + + public: + void setUp() + { + _testTools = new MeshTestToolkit(); + } + + void tearDown() + { + delete _testTools; + } + + protected: + + MeshTestToolkit* _testTools; + + }; +} +#endif diff --git a/src/INTERP_KERNEL/Test/MeshTestToolkit.cxx b/src/INTERP_KERNEL/Test/MeshTestToolkit.cxx new file mode 100644 index 000000000..576b5e6a0 --- /dev/null +++ b/src/INTERP_KERNEL/Test/MeshTestToolkit.cxx @@ -0,0 +1,341 @@ +#include "MeshTestToolkit.hxx" +#include "MEDMEM_Mesh.hxx" + +#include +#include +#include +#include +#include + +#include "VectorUtils.hxx" + +#include "MEDMEM_Field.hxx" +#include "MEDMEM_Support.hxx" + +// levels : +// 1 - titles and volume results +// 2 - symmetry / diagonal results and intersection matrix output +// 3 - empty +// 4 - empty +// 5 - misc +#include "Log.hxx" + +#include + +//#define VOL_PREC 1.0e-6 + +using namespace MEDMEM; +using namespace std; +using namespace MED_EN; + +namespace INTERP_UTILS +{ + +double MeshTestToolkit::sumRow(const IntersectionMatrix& m, int i) const +{ + double vol = 0.0; + for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter) + { + if(iter->count(i) != 0.0) + { + map::const_iterator iter2 = iter->find(i); + vol += iter2->second; + } + } + return vol; +} + +double MeshTestToolkit::sumCol(const IntersectionMatrix& m, int i) const +{ + double vol = 0.0; + const std::map& col = m[i]; + for(map::const_iterator iter = col.begin() ; iter != col.end() ; ++iter) + { + vol += std::abs(iter->second); + } + return vol; +} + + +void MeshTestToolkit::getVolumes(MESH& mesh, const double*& tab) const +{ + SUPPORT *sup=new SUPPORT(&mesh,"dummy",MED_CELL); + FIELD* f=mesh.getVolume(sup); + tab = f->getValue(); + delete sup; +} + +double MeshTestToolkit::sumVolume(const IntersectionMatrix& m) const +{ + + vector volumes; + for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter) + { + for(map::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2) + { + volumes.push_back(iter2->second); + // vol += std::abs(iter2->second); + } + } + + // sum in ascending order to avoid rounding errors + + sort(volumes.begin(), volumes.end()); + const double vol = accumulate(volumes.begin(), volumes.end(), 0.0); + + return vol; +} + +bool MeshTestToolkit::testVolumes(const IntersectionMatrix& m, MESH& sMesh, MESH& tMesh) const +{ + bool ok = true; + + // source elements + const double* sVol = new double[sMesh.getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS)]; + getVolumes(sMesh, sVol); + + for(int i = 0; i < sMesh.getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS); ++i) + { + const double sum_row = sumRow(m, i+1); + if(!epsilonEqualRelative(sum_row, sVol[i], VOL_PREC)) + { + LOG(1, "Source volume inconsistent : vol of cell " << i << " = " << sVol[i] << " but the row sum is " << sum_row ); + ok = false; + } + LOG(1, "diff = " <::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2) + { + int j = iter2->first; + const double v1 = iter2->second; + //if(m2[j - 1].count(i+1) > 0) + // { + map theMap = m2.at(j-1); + const double v2 = theMap[i + 1]; + if(v1 != v2) + { + LOG(2, "V1( " << i << ", " << j << ") = " << v1 << " which is different from V2( " << j - 1 << ", " << i + 1 << ") = " << v2 << " | diff = " << v1 - v2 ); + if(!epsilonEqualRelative(v1, v2, VOL_PREC)) + { + LOG(2, "(" << i << ", " << j << ") fails"); + isSymmetric = false; + } + } + } + ++i; + } + if(!isSymmetric) + { + LOG(1, "*** matrices are not symmetric"); + } + return isSymmetric; +} + +bool MeshTestToolkit::testDiagonal(const IntersectionMatrix& m) const +{ + LOG(1, "Checking if matrix is diagonal" ); + int i = 1; + bool isDiagonal = true; + for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter) + { + for(map::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2) + { + int j = iter2->first; + const double vol = iter2->second; + if(vol != 0.0 && (i != j)) + { + LOG(2, "V( " << i - 1 << ", " << j << ") = " << vol << " which is not zero" ); + if(!epsilonEqual(vol, 0.0, VOL_PREC)) + { + LOG(2, "(" << i << ", " << j << ") fails"); + isDiagonal = false; + } + } + } + ++i; + } + if(!isDiagonal) + { + LOG(1, "*** matrix is not diagonal"); + } + return isDiagonal; +} + +void MeshTestToolkit::dumpIntersectionMatrix(const IntersectionMatrix& m) const +{ + int i = 0; + std::cout << "Intersection matrix is " << endl; + for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter) + { + for(map::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2) + { + + std::cout << "V(" << i << ", " << iter2->first << ") = " << iter2->second << endl; + + } + ++i; + } + std::cout << "Sum of volumes = " << sumVolume(m) << std::endl; +} + +void MeshTestToolkit::calcIntersectionMatrix(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, IntersectionMatrix& m) const +{ + const string dataDir = getenv("DATA_DIR"); + + LOG(1, std::endl << "=== -> intersecting src = " << mesh1 << ", target = " << mesh2 ); + + LOG(5, "Loading " << mesh1 << " from " << mesh1path); + MESH sMesh(MED_DRIVER, dataDir+mesh1path, mesh1); + + LOG(5, "Loading " << mesh2 << " from " << mesh2path); + MESH tMesh(MED_DRIVER, dataDir+mesh2path, mesh2); + + m = _interpolator->interpol_maillages(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 ); + +} + +void MeshTestToolkit::intersectMeshes(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, const double correctVol, const double prec, bool doubleTest) const +{ + LOG(1, std::endl << std::endl << "=============================" ); + + using std::string; + const string path1 = string(mesh1path) + string(mesh1); + const string path2 = string(mesh2path) + string(mesh2); + + const bool isTestReflexive = (path1.compare(path2) == 0); + + IntersectionMatrix matrix1; + calcIntersectionMatrix(mesh1path, mesh1, mesh2path, mesh2, matrix1); + +#if LOG_LEVEL >= 2 + dumpIntersectionMatrix(matrix1); +#endif + + std::cout.precision(16); + + const double vol1 = sumVolume(matrix1); + + if(!doubleTest) + { + LOG(1, "vol = " << vol1 <<" correctVol = " << correctVol ); + CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol1, prec * std::max(correctVol, vol1)); + + if(isTestReflexive) + { + CPPUNIT_ASSERT_EQUAL_MESSAGE("Reflexive test failed", true, testDiagonal(matrix1)); + } + } + else + { + + IntersectionMatrix matrix2; + calcIntersectionMatrix(mesh2path, mesh2, mesh1path, mesh1, matrix2); + +#if LOG_LEVEL >= 2 + dumpIntersectionMatrix(matrix2); +#endif + + const double vol2 = sumVolume(matrix2); + + LOG(1, "vol1 = " << vol1 << ", vol2 = " << vol2 << ", correctVol = " << correctVol ); + + CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol1, prec * std::max(vol1, correctVol)); + CPPUNIT_ASSERT_DOUBLES_EQUAL(correctVol, vol2, prec * std::max(vol2, correctVol)); + CPPUNIT_ASSERT_DOUBLES_EQUAL(vol1, vol2, prec * std::max(vol1, vol2)); + CPPUNIT_ASSERT_EQUAL_MESSAGE("Symmetry test failed", true, testSymmetric(matrix1, matrix2)); + } + +} + +std::pair MeshTestToolkit::countNumberOfMatrixEntries(const IntersectionMatrix& m) +{ + + int numElems = 0; + int numNonZero = 0; + for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter) + { + numElems += iter->size(); + for(map::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2) + { + if(!epsilonEqual(iter2->second, 0.0, VOL_PREC)) + { + ++numNonZero; + } + } + } + return std::make_pair(numElems, numNonZero); +} + + + +} diff --git a/src/INTERP_KERNEL/Test/MeshTestToolkit.hxx b/src/INTERP_KERNEL/Test/MeshTestToolkit.hxx new file mode 100644 index 000000000..8cb133301 --- /dev/null +++ b/src/INTERP_KERNEL/Test/MeshTestToolkit.hxx @@ -0,0 +1,54 @@ +#ifndef __TU_MESH_TEST_TOOLKIT_HXX__ +#define __TU_MESH_TEST_TOOLKIT_HXX__ + +#include "../Interpolation3D.hxx" + +#define ERR_TOL 1.0e-8 + +class MEDMEM::Interpolation3D; +using MEDMEM::Interpolation3D; +class MEDMEM::MESH; + +namespace INTERP_UTILS +{ + +class MeshTestToolkit +{ + +public: + MeshTestToolkit() : _interpolator(new Interpolation3D()) {} + + ~MeshTestToolkit() { delete _interpolator; } + + // 1.0e-5 here is due to limited precision of "correct" volumes calculated in Salome + void intersectMeshes(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, const double correctVol, const double prec = 1.0e-5, bool doubleTest = true) const; + + void dumpIntersectionMatrix(const IntersectionMatrix& m) const; + + double sumRow(const IntersectionMatrix& m, int i) const; + + double sumCol(const IntersectionMatrix& m, int i) const; + + void getVolumes( MEDMEM::MESH& mesh,const double*& tab) const; + + bool testVolumes(const IntersectionMatrix& m, MEDMEM::MESH& sMesh, MEDMEM::MESH& tMesh) const; + + double sumVolume(const IntersectionMatrix& m) const; + + bool areCompatitable( const IntersectionMatrix& m1, const IntersectionMatrix& m2) const; + + bool testSymmetric(const IntersectionMatrix& m1, const IntersectionMatrix& m2) const; + + bool testDiagonal(const IntersectionMatrix& m) const; + + void calcIntersectionMatrix(const char* mesh1path, const char* mesh1, const char* mesh2path, const char* mesh2, IntersectionMatrix& m) const; + + std::pair countNumberOfMatrixEntries(const IntersectionMatrix& m); + +protected: + + Interpolation3D* _interpolator; + +}; +} +#endif diff --git a/src/INTERP_KERNEL/Test/MultiElementTetraTests.hxx b/src/INTERP_KERNEL/Test/MultiElementTetraTests.hxx new file mode 100644 index 000000000..5812a2bfa --- /dev/null +++ b/src/INTERP_KERNEL/Test/MultiElementTetraTests.hxx @@ -0,0 +1,96 @@ +#ifndef __MULTI_ELEMENT_TETRA_TESTS_HXX_ +#define __MULTI_ELEMENT_TETRA_TESTS_HXX_ + +#include "Interpolation3DTestSuite.hxx" + +namespace INTERP_UTILS +{ + class MultiElementTetraTests : public Interpolation3DTestSuite + { + CPPUNIT_TEST_SUITE( MultiElementTetraTests ); + + CPPUNIT_TEST( dividedUnitTetraSimplerReflexive ); + CPPUNIT_TEST( dividedUnitTetraReflexive ); + 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( moderateBoxSmallReflexive ); + CPPUNIT_TEST( moderateBoxEvenSmallerReflexive ); + CPPUNIT_TEST( tinyBoxReflexive ); + + CPPUNIT_TEST_SUITE_END(); + public: + void dividedUnitTetraSimplerReflexive() + { + _testTools->intersectMeshes("meshes/DividedUnitTetraSimpler.med", "DividedUnitTetraSimpler", "meshes/DividedUnitTetraSimpler.med", "DividedUnitTetraSimpler", 0.1666667); + } + + void dividedUnitTetraReflexive() + { + _testTools->intersectMeshes("meshes/DividedUnitTetra.med", "DividedUnitTetra", "meshes/DividedUnitTetra.med", "DividedUnitTetra", 0.1666667); + } + + + void nudgedDividedUnitTetra() + { + _testTools->intersectMeshes("meshes/NudgedDividedUnitTetra.med", "NudgedDividedUnitTetra", "meshes/DividedUnitTetra.med", "DividedUnitTetra", 0.150191); + } + + void nudgedDividedUnitTetraSimpler() + { + _testTools->intersectMeshes("meshes/NudgedDividedUnitTetraSimpler.med", "NudgedDividedUnitTetraSimpler", "meshes/DividedUnitTetraSimpler.med", "DividedUnitTetraSimpler", 0.150191); + } + + void dividedGenTetra() + { + _testTools->intersectMeshes("meshes/DividedGenTetra1.med", "DividedGenTetra1", "meshes/DividedGenTetra2.med", "DividedGenTetra2", 0.546329); + } + + void boxReflexive() + { + _testTools->intersectMeshes("meshes/Box3.med", "Box3", "meshes/Box3.med", "Box3", 13.9954); + } + + void boxReflexiveModerate() + { + _testTools->intersectMeshes("meshes/Box1Moderate.med", "Box1Moderate", "meshes/Box1Moderate.med", "Box1Moderate", 1.0e6); + } + + void tetraBoxes() + { + _testTools->intersectMeshes("meshes/Box1.med", "Box1", "meshes/Box2.med", "Box2", 124.197); + } + + void moderateBoxes() + { + _testTools->intersectMeshes("meshes/Box1Moderate.med", "Box1Moderate", "meshes/Box2Moderate.med", "Box2Moderate", 376856); + } + + void moderateBoxesSmaller() + { + _testTools->intersectMeshes("meshes/BoxModSmall1.med", "BoxModSmall1", "meshes/BoxModSmall2.med", "BoxModSmall2", 321853); + } + + void moderateBoxSmallReflexive() + { + _testTools->intersectMeshes("meshes/BoxModSmall1.med", "BoxModSmall1", "meshes/BoxModSmall1.med", "BoxModSmall1", 1.44018e6); + } + + void moderateBoxEvenSmallerReflexive() + { + _testTools->intersectMeshes("meshes/BoxEvenSmaller1.med", "BoxEvenSmaller1", "meshes/BoxEvenSmaller1.med", "BoxEvenSmaller1", 1.44018e6); + } + + void tinyBoxReflexive() + { + _testTools->intersectMeshes("meshes/TinyBox.med", "TinyBox", "meshes/TinyBox.med", "TinyBox", 979200); + } + }; +} + +#endif diff --git a/src/INTERP_KERNEL/Test/SingleElementTetraTests.hxx b/src/INTERP_KERNEL/Test/SingleElementTetraTests.hxx new file mode 100644 index 000000000..2e8ff387a --- /dev/null +++ b/src/INTERP_KERNEL/Test/SingleElementTetraTests.hxx @@ -0,0 +1,121 @@ +#ifndef __SINGLE_ELEMENT_TETRA_TESTS_HXX_ +#define __SINGLE_ELEMENT_TETRA_TESTS_HXX_ + +#include "Interpolation3DTestSuite.hxx" + +namespace INTERP_UTILS +{ + class SingleElementTetraTests : public Interpolation3DTestSuite + { + CPPUNIT_TEST_SUITE( SingleElementTetraTests ); + + CPPUNIT_TEST( tetraReflexiveUnit ); + CPPUNIT_TEST( tetraReflexiveGeneral ); + CPPUNIT_TEST( tetraNudgedSimpler ); + CPPUNIT_TEST( tetraNudged ); + CPPUNIT_TEST( tetraCorner ); + CPPUNIT_TEST( tetraSimpleIncluded ); + CPPUNIT_TEST( tetraDegenEdge ); + CPPUNIT_TEST( tetraDegenFace ); + CPPUNIT_TEST( tetraDegenTranslatedInPlane ); + CPPUNIT_TEST( tetraComplexIncluded ); + CPPUNIT_TEST( tetraHalfstripOnly ); + CPPUNIT_TEST( tetraHalfstripOnly2 ); + CPPUNIT_TEST( tetraSimpleHalfstripOnly ); + CPPUNIT_TEST( generalTetra ); + CPPUNIT_TEST( trickyTetra1 ); + CPPUNIT_TEST( inconsistentTetra ); + + CPPUNIT_TEST_SUITE_END(); + + public: + + void tetraReflexiveUnit() + { + _testTools->intersectMeshes("meshes/UnitTetra.med", "UnitTetra", "meshes/UnitTetra.med", "UnitTetra", 1.0/6.0); + } + + void tetraReflexiveGeneral() + { + _testTools->intersectMeshes("meshes/GeneralTetra.med", "GeneralTetra", "meshes/GeneralTetra.med", "GeneralTetra", 0.428559); + } + + void tetraNudgedSimpler() + { + _testTools->intersectMeshes("meshes/UnitTetra.med", "UnitTetra", "meshes/NudgedSimpler.med", "NudgedSimpler", 0.152112); + } + + void tetraNudged() + { + _testTools->intersectMeshes("meshes/UnitTetra.med", "UnitTetra", "meshes/NudgedTetra.med", "NudgedTetra", 0.142896); + } + + void tetraCorner() + { + _testTools->intersectMeshes("meshes/UnitTetra.med", "UnitTetra", "meshes/CornerTetra.med", "CornerTetra", 0.0135435); + } + + void tetraSimpleIncluded() + { + _testTools->intersectMeshes("meshes/SimpleIncludedTetra.med", "SimpleIncludedTetra", "meshes/SimpleIncludingTetra.med", "SimpleIncludingTetra", 17.0156); + } + + void tetraDegenEdge() + { + _testTools->intersectMeshes("meshes/UnitTetraDegenT.med", "UnitTetraDegenT", "meshes/DegenEdgeXY.med", "DegenEdgeXY", 0.0); + } + + void tetraDegenFace() + { + _testTools->intersectMeshes("meshes/UnitTetraDegenT.med", "UnitTetraDegenT", "meshes/DegenFaceXYZ.med", "DegenFaceXYZ", 0.0); + } + + void tetraDegenTranslatedInPlane() + { + _testTools->intersectMeshes("meshes/UnitTetraDegenT.med", "UnitTetraDegenT", "meshes/DegenTranslatedInPlane.med", "DegenTranslatedInPlane", 0.0571667); + } + + void tetraComplexIncluded() + { + _testTools->intersectMeshes("meshes/ComplexIncludedTetra.med", "ComplexIncludedTetra", "meshes/ComplexIncludingTetra.med", "ComplexIncludingTetra", 17.0156); + } + + void tetraHalfstripOnly() + { + // NB this test is not completely significant : we should also verify that + // there are triangles on the element that give a non-zero volume + _testTools->intersectMeshes("meshes/HalfstripOnly.med", "HalfstripOnly", "meshes/UnitTetra.med", "UnitTetra", 0.0); + } + + void tetraHalfstripOnly2() + { + // NB this test is not completely significant : we should also verify that + // there are triangles on the element that give a non-zero volume + _testTools->intersectMeshes("meshes/HalfstripOnly2.med", "HalfstripOnly2", "meshes/UnitTetra.med", "UnitTetra", 0.0); + } + + void tetraSimpleHalfstripOnly() + { + // NB this test is not completely significant : we should also verify that + // there are triangles on the element that give a non-zero volume + _testTools->intersectMeshes("meshes/SimpleHalfstripOnly.med", "SimpleHalfstripOnly", "meshes/UnitTetra.med", "UnitTetra", 0.0); + } + + void generalTetra() + { + _testTools->intersectMeshes("meshes/GenTetra1.med", "GenTetra1", "meshes/GenTetra2.med", "GenTetra2", 4.91393); + } + + void trickyTetra1() + { + _testTools->intersectMeshes("meshes/UnitTetra.med", "UnitTetra", "meshes/TrickyTetra1.med", "TrickyTetra1", 0.0); + } + + void inconsistentTetra() + { + _testTools->intersectMeshes("meshes/LargeUnitTetra.med", "LargeUnitTetra", "meshes/LargeInconsistentTetra.med", "LargeInconsistent", 7.86231e7); + } + + }; +} +#endif diff --git a/src/INTERP_KERNEL/Test/TransformedTriangleIntersectTest.cxx b/src/INTERP_KERNEL/Test/TransformedTriangleIntersectTest.cxx index 6245d9fa8..d5e3f40c5 100644 --- a/src/INTERP_KERNEL/Test/TransformedTriangleIntersectTest.cxx +++ b/src/INTERP_KERNEL/Test/TransformedTriangleIntersectTest.cxx @@ -435,6 +435,21 @@ void TransformedTriangleIntersectTest::testTriangle3() // listed with yes in the tables above return true and // that the ones listed with no or not listed at all return false + +#ifdef OPTIMIZE + bool isZero[TT::NO_TRI_SEGMENT * TT::NO_DP]; + + for(TriSegment seg = TT::PQ ; seg < TT::NO_TRI_SEGMENT ; seg = TT::TriSegment(seg + 1)) + { + // check beforehand which double-products are zero + for(DoubleProduct dp = TT::C_YZ; dp < TT::NO_DP; dp = DoubleProduct(dp + 1)) + { + isZero[TT::NO_DP*int(seg) + int(dp)] = (tri->calcStableC(seg, dp) == 0.0); + } + } +#endif + + // corner in tetrahedron (3 possibilities) CPPUNIT_ASSERT_EQUAL(true , tri->testCornerInTetrahedron(TT::P)); CPPUNIT_ASSERT_EQUAL(false, tri->testCornerInTetrahedron(TT::Q)); -- 2.39.2