*/
TriangleFaceKey(int node1, int node2, int node3)
{
- sort3Ints(_nodes, node1, node2, node3);
+ Sort3Ints(_nodes, node1, node2, node3);
_hashVal = ( _nodes[0] + _nodes[1] + _nodes[2] ) % 29;
}
return _hashVal;
}
- inline void sort3Ints(int* sorted, int node1, int node2, int node3);
+ inline static void Sort3Ints(int* sorted, int node1, int node2, int node3);
private:
/// global numbers of the three nodes, sorted in ascending order
* @param x2 second integer
* @param x3 third integer
*/
- inline void TriangleFaceKey::sort3Ints(int* sorted, int x1, int x2, int x3)
+ inline void TriangleFaceKey::Sort3Ints(int* sorted, int x1, int x2, int x3)
{
if(x1 < x2)
{
void clearVolumesCache();
private:
- // member functions
- inline void createAffineTransform(const double** corners);
inline void checkIsOutside(const double* pt, bool* isOutside, const double errTol = DEFAULT_ABS_TOL) const;
inline void checkIsStrictlyOutside(const double* pt, bool* isStrictlyOutside, const double errTol = DEFAULT_ABS_TOL) const;
inline void calculateNode(typename MyMeshType::MyConnType globalNodeNum);
HashMap< int , double* > _nodes;
/// HashMap relating triangular faces to calculated volume contributions, used for caching
- HashMap< TriangleFaceKey, double
-// #ifdef WIN32
-// , hash_compare<TriangleFaceKey,TriangleFaceKeyComparator>
-// #endif
- > _volumes;
+ HashMap< TriangleFaceKey, double > _volumes;
/// reference to the source mesh
const MyMeshType& _src_mesh;
static const double SPARSE_TRUNCATION_LIMIT;
};
- /**
- * Creates the affine transform _t from the corners of the tetrahedron. Used by the constructors
- *
- * @param corners double*[4] array containing pointers to four double[3] arrays with the
- * coordinates of the corners of the tetrahedron
- */
- template<class MyMeshType>
- inline void SplitterTetra<MyMeshType>::createAffineTransform(const double** corners)
- {
- // create AffineTransform from tetrahedron
- _t = new TetraAffineTransform( corners );
- }
-
/**
* Function used to filter out elements by checking if they belong to one of the halfspaces
* x <= 0, x >= 1, y <= 0, y >= 1, z <= 0, z >= 1, (indexed 0 - 7). The function updates an array of boolean variables
}
/**
+ * SplitterTetra class computes for a list of cell ids of a given mesh \a srcMesh (badly named) the intersection with a
+ * single TETRA4 cell given by \a tetraCorners (of length 4) and \a nodesId (of length 4 too). \a nodedIds is given only to establish
+ * if a partial computation of a triangle has already been performed (to increase performance).
+ *
+ * The \a srcMesh can contain polyhedron cells.
+ *
+ *
* Constructor creating object from the four corners of the tetrahedron.
*
* @param srcMesh mesh containing the source elements
_coords[6]=tetraCorners[2][0]; _coords[7]=tetraCorners[2][1]; _coords[8]=tetraCorners[2][2];
_coords[9]=tetraCorners[3][0]; _coords[10]=tetraCorners[3][1]; _coords[11]=tetraCorners[3][2];
// create the affine transform
- createAffineTransform(tetraCorners);
+ _t=new TetraAffineTransform(_coords);
}
/**
* with corners specified in pts. If the tetrahedron is degenerate or almost degenerate,
* construction succeeds, but the determinant of the transform is set to 0.
*
- * @param pts a 4x3 matrix containing 4 points (pts[0], ..., pts[3]) of 3 coordinates each
+ * @param pts a 4x3 matrix containing 4 points (P0X,P0Y,P0Z,P1X,P1Y,P1Z,...) of 3 coordinates each
*/
- TetraAffineTransform::TetraAffineTransform(const double** pts)
+ TetraAffineTransform::TetraAffineTransform(const double *pts)
{
LOG(2,"Creating transform from tetraeder : ");
- LOG(2, vToStr(pts[0]) << ", " << vToStr(pts[1]) << ", " << vToStr(pts[2]) << ", " << vToStr(pts[3]));
// three last points -> linear transform
for(int i = 0; i < 3 ; ++i)
for(int j = 0 ; j < 3 ; ++j)
{
// NB we insert columns, not rows
- _linear_transform[3*j + i] = (pts[i+1])[j] - (pts[0])[j];
+ _linear_transform[3*j + i] = pts[(i+1)*3+j] - pts[j];
}
}
// remember _linear_transform for the reverse transformation
memcpy( _back_linear_transform, _linear_transform, 9*sizeof(double));
- memcpy( _back_translation, pts[0], 3*sizeof(double));
+ memcpy( _back_translation, pts, 3*sizeof(double));
calculateDeterminant();
// and O is the position vector of the point that is mapped onto the origin
for(int i = 0 ; i < 3 ; ++i)
{
- _translation[i] = -(_linear_transform[3*i]*(pts[0])[0] + _linear_transform[3*i+1]*(pts[0])[1] + _linear_transform[3*i+2]*(pts[0])[2]) ;
+ _translation[i] = -(_linear_transform[3*i]*pts[0] + _linear_transform[3*i+1]*pts[1] + _linear_transform[3*i+2]*pts[2]) ;
}
// precalculate determinant (again after inversion of transform)
for(int i = 0; i < 4 ; ++i)
{
double v[3];
- apply(v, pts[i]);
+ apply(v, pts+3*i);
LOG(4, vToStr(v))
for(int j = 0; j < 3; ++j)
{
{
public:
- TetraAffineTransform(const double** pts);
+ TetraAffineTransform(const double *pts);
void apply(double* destPt, const double* srcPt) const;
void UnitTetraIntersectionBaryTest::test_TetraAffineTransform_reverseApply()
{
- double nodes[4][3] = { {-4.0, 9.0, 3.0 },
- {11.0, 0.0, 2.0 },
- { 0.0, 0.0, 0.0 },
- { 2.0, 1.0,10.0 }};
+ double nodes[12] = { -4.0, 9.0, 3.0,
+ 11.0, 0.0, 2.0,
+ 0.0, 0.0, 0.0,
+ 2.0, 1.0,10.0 };
// double pSrc[3] = { -4.0, 9.0, 3.0 };
double pSrc[3] = { 40., -20., 100. };
double pDest[] = {1,1,1};
- const double* n[4] = { &nodes[0][0], &nodes[1][0], &nodes[2][0], &nodes[3][0] };
- TetraAffineTransform a(&n[0]);
+ TetraAffineTransform a(nodes);
a.apply( pDest, pSrc );
a.reverseApply( pDest, pDest );
CPPUNIT_ASSERT_DOUBLES_EQUAL( pSrc[0], pDest[0], 1e-12);
sourceMesh->decrRef();
targetMesh->decrRef();
}
+
+void MEDCouplingRemapperTest::testBugNonRegression1()
+{
+ // source
+ DataArrayDouble *coordsSrc(DataArrayDouble::New());
+ const double coordsSrcData[18]={-6.25,3.6084391824351605,264.85199999999998,-6.25,3.6084391824351605,289.05200000000002,-6.2499999999999991,-3.6084391824351618,264.85199999999998,-6.2499999999999991,-3.6084391824351618,289.05200000000002,-1.7763568394002505e-15,4.4408920985006262e-15,264.85199999999998,-1.7763568394002505e-15,4.4408920985006262e-15,289.05200000000002};
+ coordsSrc->useArray(coordsSrcData,false,CPP_DEALLOC,6,3);
+ DataArrayInt *connSrc(DataArrayInt::New()),*connISrc(DataArrayInt::New());
+ const int connSrcData[7]={16,2,0,4,3,1,5};
+ connSrc->useArray(connSrcData,false,CPP_DEALLOC,7,1);
+ const int connISrcData[2]={0,7};
+ connISrc->useArray(connISrcData,false,CPP_DEALLOC,2,1);
+ MEDCouplingUMesh *srcMesh(MEDCouplingUMesh::New("source",3));
+ srcMesh->setCoords(coordsSrc);
+ srcMesh->setConnectivity(connSrc,connISrc,true);
+ coordsSrc->decrRef(); connSrc->decrRef(); connISrc->decrRef();
+ // target
+ DataArrayDouble *coordsTrg(DataArrayDouble::New());
+const double coordsTrgData[36]={-2,1.1547005383792521,264.85199999999998,-2,0.57735026918962618,264.85199999999998,-2.5,0.2886751345948132,264.85199999999998,-2.5,1.443375672974065,264.85199999999998,-3.0000000000000004,1.1547005383792526,264.85199999999998,-3.0000000000000004,0.57735026918962662,264.85199999999998,-2,1.1547005383792521,289.05200000000002,-2,0.57735026918962618,289.05200000000002,-2.5,0.2886751345948132,289.05200000000002,-2.5,1.443375672974065,289.05200000000002,-3.0000000000000004,1.1547005383792526,289.05200000000002,-3.0000000000000004,0.57735026918962662,289.05200000000002};
+ coordsTrg->useArray(coordsTrgData,false,CPP_DEALLOC,12,3);
+ DataArrayInt *connTrg=DataArrayInt::New();
+ const int connTrgData[44]={31,0,1,2,5,4,3,-1,7,6,9,10,11,8,-1,3,9,6,0,-1,4,10,9,3,-1,5,11,10,4,-1,2,8,11,5,-1,1,7,8,2,-1,0,6,7,1};
+ connTrg->useArray(connTrgData,false,CPP_DEALLOC,44,1);
+ DataArrayInt *connITrg=DataArrayInt::New();
+ const int connITrgData[2]={0,44};
+ connITrg->useArray(connITrgData,false,CPP_DEALLOC,2,1);
+ MEDCouplingUMesh *trgMesh=MEDCouplingUMesh::New("target",3);
+ trgMesh->setCoords(coordsTrg);
+ trgMesh->setConnectivity(connTrg,connITrg,true);
+ coordsTrg->decrRef(); connTrg->decrRef(); connITrg->decrRef();
+ // Go !
+ const double valExpected(20.957814771583468);
+ MEDCouplingRemapper remapper;
+ remapper.setPrecision(1e-12);
+ remapper.setIntersectionType(INTERP_KERNEL::Triangulation);
+ CPPUNIT_ASSERT_EQUAL(1,remapper.prepare(srcMesh,trgMesh,"P0P0"));
+ std::vector<std::map<int,double> > matrx(remapper.getCrudeMatrix());
+ CPPUNIT_ASSERT_EQUAL(1,(int)matrx.size());
+ CPPUNIT_ASSERT_EQUAL(1,(int)matrx[0].size());
+ std::cerr << std::endl << matrx[0][0] << std::endl;
+ //
+ srcMesh->decrRef(); trgMesh->decrRef();
+}
+
CPPUNIT_TEST( testExtruded2 );
CPPUNIT_TEST( testPrepareEx1 );
CPPUNIT_TEST( testPartialTransfer1 );
+ CPPUNIT_TEST( testBugNonRegression1 );
CPPUNIT_TEST_SUITE_END();
public:
void test2DInterpP0P0_1();
void testExtruded2();
void testPrepareEx1();
void testPartialTransfer1();
+ //
+ void testBugNonRegression1();
private:
static MEDCouplingUMesh *build1DTargetMesh_2();
static MEDCouplingUMesh *build2DTargetMesh_3();
