--- /dev/null
+#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<int>(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<IntersectorTetra*>::iterator iter = _tetra.begin(); iter != _tetra.end(); ++iter)
+ {
+ delete *iter;
+ }
+
+ // free potential sub-mesh nodes that have been allocated
+ vector<const double*>::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<IntersectorTetra*>::iterator iter = _tetra.begin(); iter != _tetra.end(); ++iter)
+ {
+ volume += (*iter)->intersectSourceCell(srcCell);
+ }
+ return volume;
+ }
+
+
+};
--- /dev/null
+#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<int n>
+ inline void calcBarycenter(double* barycenter, const int* const pts);
+
+ vector<IntersectorTetra*> _tetra;
+
+ vector<const double*> _nodes;
+
+ };
+
+
+ inline const double* IntersectorHexa::getCoordsOfSubNode(int node)
+ {
+ // replace at with [] for unsafe but faster access
+ return _nodes.at(node - 1);
+ }
+
+ template<int n>
+ 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
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
--- /dev/null
+#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
--- /dev/null
+#ifndef __TU_INTERPOLATION_3D_TEST_SUITE_HXX__
+#define __TU_INTERPOLATION_3D_TEST_SUITE_HXX__
+
+#include "MeshTestToolkit.hxx"
+
+#include <cppunit/extensions/HelperMacros.h>
+
+namespace INTERP_UTILS
+{
+
+ class Interpolation3DTestSuite : public CppUnit::TestFixture
+ {
+
+ public:
+ void setUp()
+ {
+ _testTools = new MeshTestToolkit();
+ }
+
+ void tearDown()
+ {
+ delete _testTools;
+ }
+
+ protected:
+
+ MeshTestToolkit* _testTools;
+
+ };
+}
+#endif
--- /dev/null
+#include "MeshTestToolkit.hxx"
+#include "MEDMEM_Mesh.hxx"
+
+#include <iostream>
+#include <map>
+#include <vector>
+#include <cmath>
+#include <algorithm>
+
+#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 <cppunit/extensions/HelperMacros.h>
+
+//#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<int, double>::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<int, double>& col = m[i];
+ for(map<int, double>::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<double>* f=mesh.getVolume(sup);
+ tab = f->getValue();
+ delete sup;
+}
+
+double MeshTestToolkit::sumVolume(const IntersectionMatrix& m) const
+{
+
+ vector<double> volumes;
+ for(IntersectionMatrix::const_iterator iter = m.begin() ; iter != m.end() ; ++iter)
+ {
+ for(map<int, double>::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 = " <<sum_row - sVol[i] );
+ }
+
+ // target elements
+ const double* tVol = new double[tMesh.getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS)];
+ getVolumes(tMesh, tVol);
+ for(int i = 0; i < tMesh.getNumberOfElements(MED_CELL,MED_ALL_ELEMENTS); ++i)
+ {
+ const double sum_col = sumCol(m, i);
+ if(!epsilonEqualRelative(sum_col, tVol[i], VOL_PREC))
+ {
+ LOG(1, "Target volume inconsistent : vol of cell " << i << " = " << tVol[i] << " but the col sum is " << sum_col);
+ ok = false;
+ }
+ LOG(1, "diff = " <<sum_col - tVol[i] );
+ }
+ delete[] sVol;
+ delete[] tVol;
+
+ return ok;
+}
+
+bool MeshTestToolkit::areCompatitable(const IntersectionMatrix& m1, const IntersectionMatrix& m2) const
+{
+ bool compatitable = true;
+ int i = 0;
+ for(IntersectionMatrix::const_iterator iter = m1.begin() ; iter != m1.end() ; ++iter)
+ {
+ for(map<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
+ {
+ int j = iter2->first;
+ if(m2.at(j-1).count(i+1) == 0)
+ {
+ if(!epsilonEqual(iter2->second, 0.0, VOL_PREC))
+ {
+ LOG(2, "V1( " << i << ", " << j << ") exists, but V2( " << j - 1 << ", " << i + 1 << ") " << " does not " );
+ LOG(2, "(" << i << ", " << j << ") fails");
+ compatitable = false;
+ }
+ }
+ }
+ ++i;
+ }
+ if(!compatitable)
+ {
+ LOG(1, "*** matrices are not compatitable");
+ }
+ return compatitable;
+}
+
+bool MeshTestToolkit::testSymmetric(const IntersectionMatrix& m1, const IntersectionMatrix& m2) const
+{
+
+ int i = 0;
+ bool isSymmetric = true;
+
+ LOG(1, "Checking symmetry src - target" );
+ isSymmetric = isSymmetric & areCompatitable(m1, m2) ;
+ LOG(1, "Checking symmetry target - src" );
+ isSymmetric = isSymmetric & areCompatitable(m2, m1);
+
+ for(IntersectionMatrix::const_iterator iter = m1.begin() ; iter != m1.end() ; ++iter)
+ {
+ for(map<int, double>::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<int, double> 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<int, double>::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<int, double>::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<int,int> 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<int, double>::const_iterator iter2 = iter->begin() ; iter2 != iter->end() ; ++iter2)
+ {
+ if(!epsilonEqual(iter2->second, 0.0, VOL_PREC))
+ {
+ ++numNonZero;
+ }
+ }
+ }
+ return std::make_pair(numElems, numNonZero);
+}
+
+
+
+}
--- /dev/null
+#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<int,int> countNumberOfMatrixEntries(const IntersectionMatrix& m);
+
+protected:
+
+ Interpolation3D* _interpolator;
+
+};
+}
+#endif
--- /dev/null
+#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
--- /dev/null
+#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
// 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));
