From: vbd <vbd>
Date: Fri, 3 Aug 2007 10:37:10 +0000 (+0000)
Subject: staffan :
X-Git-Tag: trio_trio_coupling~82
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=d3cbf0078b502aedbd4679e08896211fd95f760b;p=tools%2Fmedcoupling.git

staffan :
	* fixed nasty bug involving the order of elements in a table
	used for triple product calculation
	* rewrote MeshElement::triangulate using CELLMODEL -> does not
	give the correct coordinates yet it seems
---

diff --git a/src/INTERP_KERNEL/Interpolation3D.cxx b/src/INTERP_KERNEL/Interpolation3D.cxx
index 049b4c9a6..8e538178d 100644
--- a/src/INTERP_KERNEL/Interpolation3D.cxx
+++ b/src/INTERP_KERNEL/Interpolation3D.cxx
@@ -298,6 +298,13 @@ namespace MEDMEM
 	T.dump();
 	std::cout << std::endl;
 
+	if(T.determinant() == 0.0)
+	  {
+	    // tetra is planar
+	    std::cout << "Planar tetra -- volume 0" << std::endl;
+	    return 0.0;
+	  }
+
 	// triangulate source element faces (assumed tetraeder for the time being)
 	// do nothing
 	vector<TransformedTriangle> triangles;
diff --git a/src/INTERP_KERNEL/MeshElement.cxx b/src/INTERP_KERNEL/MeshElement.cxx
index 620e36d86..f57e20f27 100644
--- a/src/INTERP_KERNEL/MeshElement.cxx
+++ b/src/INTERP_KERNEL/MeshElement.cxx
@@ -2,6 +2,7 @@
 
 #include "TetraAffineTransform.hxx"
 #include "TransformedTriangle.hxx"
+#include "MEDMEM_CellModel.hxx"
 
 namespace INTERP_UTILS
 {
@@ -110,6 +111,66 @@ namespace INTERP_UTILS
    * @param      T          affine transform that is applied to the nodes of the triangles
    */
   void MeshElement::triangulate(std::vector<TransformedTriangle>& triangles, const TetraAffineTransform& T) const
+  {
+    // get cell model for the element
+    CELLMODEL cellModel(_type);
+
+    assert(cellModel.getDimension() == 3);
+
+    // start index in connectivity array for cell
+    //    const int cellIdx = _mesh->getConnectivityIndex(MED_NODAL, MED_CELL)[_index] - 1;
+
+    //    assert(cellIdx >= 0);
+    //assert(cellIdx < _mesh->getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS));
+
+    // loop over faces
+    for(int i = 1 ; i <= cellModel.getNumberOfConstituents(1) ; ++i)
+      {
+	medGeometryElement faceType = cellModel.getConstituentType(1, i);
+	CELLMODEL faceModel(faceType);
+
+	assert(faceModel.getDimension() == 2);
+
+	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)
+	  {
+	    // offset of node from cellIdx
+	    int localNodeNumber = cellModel.getNodeConstituent(1, i, j);
+
+	    assert(localNodeNumber >= 1);
+	    assert(localNodeNumber <= cellModel.getNumberOfNodes());
+
+	    const double* node = getCoordsOfNode(localNodeNumber);
+	    
+	    // transform 
+	    //{ not totally efficient since we transform each node once per face
+	    T.apply(&transformedNodes[3*(j-1)], node);
+
+	  }
+
+	assert(faceType == MED_TRIA3);
+
+	// create triangles from face
+	switch(faceType)
+	  {
+	  case MED_TRIA3:
+	    triangles.push_back(TransformedTriangle(&transformedNodes[0], &transformedNodes[3], &transformedNodes[6]));
+	    break;
+	    
+	  default:
+	    std::cout << "Only elements with triangular faces are supported at the moment." << std::endl;
+	  }
+      }
+    
+
+  }
+  
+
+#if 0
+  void MeshElement::triangulate(std::vector<TransformedTriangle>& triangles, const TetraAffineTransform& T) const
   {
     std::cout << "Triangulating element " << getIndex() << std::endl;
     switch(_type)
@@ -156,7 +217,7 @@ namespace INTERP_UTILS
 	break;
       }
   }
-
+#endif
   int MeshElement::getIndex() const
   {
     return _index + 1;
diff --git a/src/INTERP_KERNEL/TetraAffineTransform.hxx b/src/INTERP_KERNEL/TetraAffineTransform.hxx
index d8e5e6c2e..c564438d0 100644
--- a/src/INTERP_KERNEL/TetraAffineTransform.hxx
+++ b/src/INTERP_KERNEL/TetraAffineTransform.hxx
@@ -32,8 +32,21 @@ namespace INTERP_UTILS
 		_linearTransform[3*j + i] = (pts[i+1])[j] - (pts[0])[j];
 	      }
 	  }
-#if 0
+
 	calculateDeterminant();
+
+	std::cout << "determinant before inverse = " << _determinant << std::endl;
+	
+	// check that tetra is non-planar -> determinant is not zero
+	// otherwise set _determinant to zero to signal caller that transformation did not work
+	if(epsilonEqual(_determinant, 0.0))
+	  {
+	    _determinant = 0.0;
+	    return;
+	  }
+	
+#if 0
+
 	//assert(_determinant != 0.0);
 	
 	if(_determinant < 0.0) 
@@ -62,12 +75,14 @@ namespace INTERP_UTILS
       // precalculate determinant (again after inversion of transform)
       calculateDeterminant();
       
+      // self-check
+      std::cout << "transform determinant is " << _determinant << std::endl;
       std::cout << "*Self-check : Applying transformation to original points ... ";
       for(int i = 0; i < 4 ; ++i)
 	{
 	  double v[3];
 	  apply(v, pts[i]);
-	  //	  std::cout << vToStr(v) << std::endl;
+	  std::cout << vToStr(v) << std::endl;
 	  for(int j = 0; j < 3; ++j)
 	    {
 	      assert(epsilonEqual(v[j], (3*i+j == 3 || 3*i+j == 7 || 3*i+j == 11 ) ? 1.0 : 0.0));
diff --git a/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx b/src/INTERP_KERNEL/TransformedTriangle_intersect.cxx
index 966c1aebd..f91241aaf 100644
--- 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 << "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);
       }
   }
 
@@ -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
@@ -690,9 +690,9 @@ namespace INTERP_UTILS
 	X, O, // OX
 	Y, O, // OY
 	Z, O, // OZ 
+	X, Y,  // XY
 	Y, Z, // YZ
 	Z, X, // ZX
-	X, Y  // XY
       };
 
     // Grandy, [16]
@@ -700,6 +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;
     return (t1*t2 <= 0.0) && (t1 - t2 != 0.0);
   }
 
@@ -771,11 +772,11 @@ namespace INTERP_UTILS
     //? is it always YZ here ?
     //? changed to XY !
     const double normal = calcStableC(PQ, C_XY) + calcStableC(QR, C_XY) + calcStableC(RP, C_XY);
-    std::cout << "surface above corner " << corner << " : " << "n = " << normal << ", t = [" <<  calcTByDevelopingRow(corner, 1, false) << ", "  << calcTByDevelopingRow(corner, 2, false) << ", " << calcTByDevelopingRow(corner, 3, false) << "] - stable : " << calcStableT(corner)  << std::endl;
+    // std::cout << "surface above corner " << corner << " : " << "n = " << normal << ", t = [" <<  calcTByDevelopingRow(corner, 1, false) << ", "  << calcTByDevelopingRow(corner, 2, false) << ", " << calcTByDevelopingRow(corner, 3, false) << "] - stable : " << calcStableT(corner)  << std::endl;
     //? we don't care here if the triple product is "invalid", that is, the triangle does not surround one of the
     // edges going out from the corner (Grandy [53])
-        return ( calcTByDevelopingRow(corner, 2, false) * normal ) >= 0.0;
-    //   return ( calcStableT(corner) * normal ) >= 0.0;
+    //   return ( calcTByDevelopingRow(corner, 2, false) * normal ) >= 0.0;
+       return ( calcStableT(corner) * normal ) >= 0.0;
   }
 
   /**
diff --git a/src/INTERP_KERNEL/TransformedTriangle_math.cxx b/src/INTERP_KERNEL/TransformedTriangle_math.cxx
index 0980b3b20..c513a2432 100644
--- a/src/INTERP_KERNEL/TransformedTriangle_math.cxx
+++ b/src/INTERP_KERNEL/TransformedTriangle_math.cxx
@@ -215,10 +215,10 @@ namespace INTERP_UTILS
     if(_isTripleProductsCalculated)
       return;
 
-    std::cout << "Precalculating triple products" << std::endl;
+    // std::cout << "Precalculating triple products" << std::endl;
     for(TetraCorner corner = O ; corner <= Z ; corner = TetraCorner(corner + 1))
       {
-	std::cout << "- Triple product for corner " << corner << std::endl;
+	// std::cout << "- Triple product for corner " << corner << std::endl;
 	bool isGoodRowFound = false;
 
 	// find edge / row to use
@@ -309,7 +309,7 @@ namespace INTERP_UTILS
 	  {
 	    // this value will not be used
 	    // we set it to whatever
-	    // std::cout << "Triple product not calculated for corner " << corner << std::endl;
+	    std::cout << "Triple product not calculated for corner " << corner << std::endl;
 	    _tripleProducts[corner] = -3.14159265;
 	    _validTP[corner] = false;