staffan :

	* correction of the MeshElement::getCoordsOfNode() which was
	totally wrong - it worked well with only one element because
	of the correspondance between the arrays connectivityIndex
	and connectivity that exists in this case

diff --git a/src/INTERP_KERNEL/Interpolation3D.cxx b/src/INTERP_KERNEL/Interpolation3D.cxx
index 8e538178deb2912ca42965cc79979bd35a2c5753..9355fbd27233ecb0a02fbe2bfc8bc27f4486ff2f 100644 (file)
--- a/src/INTERP_KERNEL/Interpolation3D.cxx
+++ b/src/INTERP_KERNEL/Interpolation3D.cxx
@@ -280,11 +280,16 @@ namespace MEDMEM
     double Interpolation3D::calculateIntersectionVolume(const MeshElement& srcElement, const MeshElement& targetElement)
       {
 
-       // std::cout << "Intersecting elems " << srcElement.getIndex() << " and " << targetElement.getIndex() << std::endl;
+       //std::cout << "Intersecting elems " << srcElement.getIndex() << " and " << targetElement.getIndex() << std::endl;
        // (a), (b) and (c) not yet implemented
 
        // (d) : without fine-level filtering (a) - (c) for the time being
-       
+
+       std::cout << "Source : ";
+       srcElement.dumpCoords();
+       std::cout << "Target : ";
+       targetElement.dumpCoords();
+
        // get array of points of target tetraeder
        const double* tetraCorners[4];
        for(int i = 0 ; i < 4 ; ++i)
@@ -294,10 +299,8 @@ namespace MEDMEM
 
        // create AffineTransform
        TetraAffineTransform T( tetraCorners );
-       std::cout << "Transform : " << std::endl;
-       T.dump();
-       std::cout << std::endl;
-
+       
+       // check if we have planar tetra element
        if(T.determinant() == 0.0)
          {
            // tetra is planar
@@ -305,6 +308,10 @@ namespace MEDMEM
            return 0.0;
          }
 
+       std::cout << "Transform : " << std::endl;
+       T.dump();
+       std::cout << std::endl;
+
        // triangulate source element faces (assumed tetraeder for the time being)
        // do nothing
        vector<TransformedTriangle> triangles;
@@ -328,7 +335,7 @@ namespace MEDMEM
 
        //? fault in article, Grandy, [8] : it is the determinant of the inverse transformation that 
        // should be used
-       return 1.0 / T.determinant() * volume ;
+       return std::abs(1.0 / T.determinant() * volume) ;
        
       }
 }

diff --git a/src/INTERP_KERNEL/MeshElement.cxx b/src/INTERP_KERNEL/MeshElement.cxx
index f57e20f272095127b94b097153bbf5af6e21cd6b..1501b53228dc7f2e0a6ff84ea052b8f4e731ec65 100644 (file)
--- a/src/INTERP_KERNEL/MeshElement.cxx
+++ b/src/INTERP_KERNEL/MeshElement.cxx
@@ -99,9 +99,12 @@ namespace INTERP_UTILS
    */
   const double* MeshElement::getCoordsOfNode(int node) const
   {
+    assert(node >= 1);
+    assert(node <= getNumberNodes());
     const int nodeOffset = node - 1;
     const int elemIdx = _mesh->getConnectivityIndex(MED_NODAL, MED_CELL)[_index] - 1;
-    return &(_mesh->getCoordinates(MED_FULL_INTERLACE)[elemIdx + 3*nodeOffset]);
+    const int connIdx = _mesh->getConnectivity(MED_FULL_INTERLACE, MED_NODAL, MED_CELL, MED_ALL_ELEMENTS)[elemIdx + nodeOffset] - 1;
+    return &(_mesh->getCoordinates(MED_FULL_INTERLACE)[3*connIdx]);
   }
   
   /**
@@ -130,9 +133,10 @@ namespace INTERP_UTILS
        CELLMODEL faceModel(faceType);
 
        assert(faceModel.getDimension() == 2);
+       assert(faceModel.getNumberOfNodes() == 3);
 
        double transformedNodes[3 * faceModel.getNumberOfNodes()];
-       const double* coords = _mesh->getCoordinates(MED_FULL_INTERLACE);
+       
 
        // loop over nodes of face
        for(int j = 1; j <= faceModel.getNumberOfNodes(); ++j)
@@ -149,11 +153,13 @@ namespace INTERP_UTILS
            //{ not totally efficient since we transform each node once per face
            T.apply(&transformedNodes[3*(j-1)], node);
 
+           std::cout << "Node " << localNodeNumber << " = " << vToStr(node) << " transformed to " << vToStr(&transformedNodes[3*(j-1)]) << std::endl;
+
          }
 
        assert(faceType == MED_TRIA3);
 
-       // create triangles from face
+       // create transformed triangles from face
        switch(faceType)
          {
          case MED_TRIA3:
@@ -218,11 +224,22 @@ namespace INTERP_UTILS
       }
   }
 #endif
+
   int MeshElement::getIndex() const
   {
     return _index + 1;
   }
   
+  void MeshElement::dumpCoords() const
+    {
+      std::cout << "Element " << _index + 1 << " has nodes " << std::endl;
+      for(int i = 1 ; i <= getNumberNodes() ; ++i)
+       {
+         std::cout << vToStr(getCoordsOfNode(i)) << ", ";
+       }
+      std::cout << std::endl;
+    }
+
 
   /// ElementBBoxOrder
   bool ElementBBoxOrder::operator()( MeshElement* elem1, MeshElement* elem2)

diff --git a/src/INTERP_KERNEL/MeshElement.hxx b/src/INTERP_KERNEL/MeshElement.hxx
index ab339da23d8e1a0e6ad0396f00afa11c1fb675d1..f58c36d0fae7962cad3db93f0552b21b21f624cf 100644 (file)
--- a/src/INTERP_KERNEL/MeshElement.hxx
+++ b/src/INTERP_KERNEL/MeshElement.hxx
@@ -97,6 +97,8 @@ namespace INTERP_UTILS
     
     int getIndex() const;
 
+    void dumpCoords() const;
+
   private:
     const int _index;
     const MEDMEM::MESH* _mesh;

diff --git a/src/INTERP_KERNEL/Test/Interpolation3DTest.cxx b/src/INTERP_KERNEL/Test/Interpolation3DTest.cxx
index c887c75d038027bd1f785e1b565e7087296600b6..3da899fd9da7c2c433e0bd15fc6e2a5fde90c414 100644 (file)
--- a/src/INTERP_KERNEL/Test/Interpolation3DTest.cxx
+++ b/src/INTERP_KERNEL/Test/Interpolation3DTest.cxx
@@ -105,8 +105,6 @@ void Interpolation3DTest::calcIntersectionMatrix(const char* mesh1path, const ch
 
   m = interpolator->interpol_maillages(sMesh, tMesh);
 
-  dumpIntersectionMatrix(m);
-
   std::cout << "Intersection calculation done. " << std::endl << std::endl;
 }
 

diff --git a/src/INTERP_KERNEL/TetraAffineTransform.hxx b/src/INTERP_KERNEL/TetraAffineTransform.hxx
index c564438d0ee04ee59341acb905b896df569a47d4..159b973acbda4de28ee9e504d345a92a96b44a7d 100644 (file)
--- a/src/INTERP_KERNEL/TetraAffineTransform.hxx
+++ b/src/INTERP_KERNEL/TetraAffineTransform.hxx
@@ -17,8 +17,8 @@ namespace INTERP_UTILS
     TetraAffineTransform(const double** pts)
     {
 
-      std::cout << "Creating transform from tetraeder : " << std::endl;
-      std::cout << vToStr(pts[0]) << ", " << vToStr(pts[1]) << ", " << vToStr(pts[2]) << ", " << vToStr(pts[3]) << ", " << std::endl;
+      // std::cout << "Creating transform from tetraeder : " << std::endl;
+      // std::cout << vToStr(pts[0]) << ", " << vToStr(pts[1]) << ", " << vToStr(pts[2]) << ", " << vToStr(pts[3]) << ", " << std::endl;
 
 #if 0
       do {

diff --git a/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx b/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx
index f91241aafd7bfb9e0851a1f5da7c6963ac6a4679..b601b7fe25e94eec82195fc7474ed95bfb37bd5f 100644 (file)
--- a/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx
+++ b/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx
@@ -109,12 +109,12 @@ namespace INTERP_UTILS
     // calculate point
     for(int i = 0; i < 3; ++i)
       {
-       std::cout << "tA = " << tA << " tB = " << tB << " alpha= " << alpha << std::endl;
+       // std::cout << "tA = " << tA << " tB = " << tB << " alpha= " << alpha << std::endl;
        pt[i] = (1 - alpha) * COORDS_TET_CORNER[3*corners[0] + i] + 
          alpha * COORDS_TET_CORNER[3*corners[1] + i];
-        std::cout << pt[i] << std::endl;
-        //assert(pt[i] >= 0.0);
-        //assert(pt[i] <= 1.0);
+       // std::cout << pt[i] << std::endl;
+        assert(pt[i] >= 0.0);
+        assert(pt[i] <= 1.0);
       }
   }
 
@@ -429,8 +429,8 @@ namespace INTERP_UTILS
     
     // special case : facet H = 0
     bool cond2 = (facet == NO_TET_FACET) ? testSegmentIntersectsHPlane(seg) : testSegmentIntersectsFacet(seg, facet);
-    std::cout << "Halfstrip tests (" << seg << ", " << edge << ") : " << (cVals[0]*cVals[1] < 0.0) << ", " << cond2 << ", " << (cVals[2]*cVals[3] > 0.0) << std::endl;
-    std::cout << "c2 = " << cVals[2] << ", c3 = " << cVals[3] << std::endl; 
+    // std::cout << "Halfstrip tests (" << seg << ", " << edge << ") : " << (cVals[0]*cVals[1] < 0.0) << ", " << cond2 << ", " << (cVals[2]*cVals[3] > 0.0) << std::endl;
+    // std::cout << "c2 = " << cVals[2] << ", c3 = " << cVals[3] << std::endl; 
   
     return (cVals[0]*cVals[1] < 0.0) && cond2 && (cVals[2]*cVals[3] > 0.0);
   }
@@ -657,7 +657,7 @@ namespace INTERP_UTILS
     const double cQR = calcStableC(QR, DoubleProduct(edge));
     const double cRP = calcStableC(RP, DoubleProduct(edge));
 
-    std::cout << "TriangleSurroundsEdge : edge = " << edge << " c = [" << cPQ << ", " << cQR << ", " << cRP << "]" << std::endl;
+    // std::cout << "TriangleSurroundsEdge : edge = " << edge << " c = [" << cPQ << ", " << cQR << ", " << cRP << "]" << std::endl;
 
     // if two or more c-values are zero we disallow x-edge intersection
     // Grandy, p.446
@@ -665,7 +665,7 @@ namespace INTERP_UTILS
     
     if(numZeros >= 2 ) 
       {
-       std::cout << "TriangleSurroundsEdge test fails due to too many 0 dp" << std::endl; 
+       // std::cout << "TriangleSurroundsEdge test fails due to too many 0 dp" << std::endl; 
       }
 
     return (cPQ*cQR >= 0.0) && (cQR*cRP >= 0.0) && (cRP*cPQ >= 0.0) && numZeros < 2;
@@ -700,7 +700,7 @@ namespace INTERP_UTILS
     const double t2 = calcStableT(TRIPLE_PRODUCTS[2*edge + 1]);
 
     //? should equality with zero use epsilon?
-    std::cout << "testEdgeIntersectsTriangle : t1 = " << t1 << " t2 = " << t2 << std::endl;
+    // std::cout << "testEdgeIntersectsTriangle : t1 = " << t1 << " t2 = " << t2 << std::endl;
     return (t1*t2 <= 0.0) && (t1 - t2 != 0.0);
   }
 
@@ -807,7 +807,7 @@ namespace INTERP_UTILS
 
     //? NB here we have no correction for precision - is this good?
     // Our authority Grandy says nothing
-    std::cout << "dp in triSurrRay for corner " << corner << " = [" << cPQ << ", " << cQR << ", " << cRP << "]" << std::endl;
+    // std::cout << "dp in triSurrRay for corner " << corner << " = [" << cPQ << ", " << cQR << ", " << cRP << "]" << std::endl;
     return ( cPQ*cQR > 0.0 ) && ( cPQ*cRP > 0.0 );
 
   }