optimizing Intersector3D

diff --git a/src/INTERP_KERNEL/Intersector3D.cxx b/src/INTERP_KERNEL/Intersector3D.cxx
index 9fd4194423af15382ff9ee0fd4f405b2e5f5f81c..605d09a845a9927614b1839fb74ba478ce9760af 100644 (file)
--- a/src/INTERP_KERNEL/Intersector3D.cxx
+++ b/src/INTERP_KERNEL/Intersector3D.cxx
@@ -28,11 +28,15 @@ namespace MEDMEM
   }
 
   /*
-   * Destructor
+   * Destruct
    *
    */
   Intersector3D::~Intersector3D()
   {
+               for (map<medGeometryElement, CELLMODEL*>::iterator it=_cellmodel_map.begin();
+                       it!=_cellmodel_map.end();
+                       it++)
+                               delete it->second;
   }
 
   /*
@@ -124,23 +128,23 @@ namespace MEDMEM
   {
 
     // get cell model for the element
-    CELLMODEL cellModel(type);
+    CELLMODEL* cellModel=retrieveCellModel(type);
 
-    assert(cellModel.getDimension() == 3);
+    assert(cellModel->getDimension() == 3);
     assert(element >= 1);
     assert(element <= _srcMesh.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS));
 
     // try to avoid reallocations with push_back()
     // 2 * the number of faces should be plenty
-    triangles.reserve(2 * cellModel.getNumberOfConstituents(1));
+    triangles.reserve(2 * cellModel->getNumberOfConstituents(1));
 
     // loop over faces
-    double** transformedNodes = new double*[cellModel.getNumberOfNodes()];
+    double** transformedNodes = new double*[cellModel->getNumberOfNodes()];
     bool isOutside[8] = {true, true, true, true, true, true, true, true};
     bool isTargetOutside = false;
     
     // calculate the coordinates of the nodes
-    for(int i = 1; i <= cellModel.getNumberOfNodes() ; ++i)
+    for(int i = 1; i <= cellModel->getNumberOfNodes() ; ++i)
       {
        const double* node = getCoordsOfNode(i, element, _srcMesh);
        transformedNodes[i-1] = new double[3];
@@ -165,9 +169,9 @@ namespace MEDMEM
 
     if(!isTargetOutside)
       {
-       for(int i = 1 ; i <= cellModel.getNumberOfConstituents(1) ; ++i)
+       for(int i = 1 ; i <= cellModel->getNumberOfConstituents(1) ; ++i)
          {
-           const medGeometryElement faceType = cellModel.getConstituentType(1, i);
+           const medGeometryElement faceType = cellModel->getConstituentType(1, i);
            CELLMODEL faceModel(faceType);
        
            assert(faceModel.getDimension() == 2);
@@ -177,10 +181,10 @@ namespace MEDMEM
            // get the nodes of the face
            for(int j = 1; j <= faceModel.getNumberOfNodes(); ++j)
              {
-               nodes[j-1] = cellModel.getNodeConstituent(1, i, j) - 1;
+               nodes[j-1] = cellModel->getNodeConstituent(1, i, j) - 1;
                
                assert(nodes[j-1] +1 >= 0);
-               assert(nodes[j-1] +1 <= cellModel.getNumberOfNodes());
+               assert(nodes[j-1] +1 <= cellModel->getNumberOfNodes());
              }
            //       std::cout << "here3 : " << nodes[0] << "," <<  nodes[1] << ", " <<  nodes[2] << std::endl;
            // create transformed triangles from face
@@ -226,7 +230,7 @@ namespace MEDMEM
        ++filtered;
       }
 
-    for(int i = 0 ; i < cellModel.getNumberOfNodes() ; ++i)
+    for(int i = 0 ; i < cellModel->getNumberOfNodes() ; ++i)
       {
        delete[] transformedNodes[i];
       }
@@ -239,16 +243,16 @@ namespace MEDMEM
   void Intersector3D::triangulate(const medGeometryElement type, const int element, std::vector<TransformedTriangle>& triangles, const TetraAffineTransform& T) const
   {
     // get cell model for the element
-    CELLMODEL cellModel(type);
+    CELLMODEL* cellModel=retrieveCellModel(type);
 
     assert(element >= 1);
     assert(element <= _srcMesh.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS));
-    assert(cellModel.getDimension() == 3);
+    assert(cellModel->getDimension() == 3);
 
     // loop over faces
-    for(int i = 1 ; i <= cellModel.getNumberOfConstituents(1) ; ++i)
+    for(int i = 1 ; i <= cellModel->getNumberOfConstituents(1) ; ++i)
       {
-       medGeometryElement faceType = cellModel.getConstituentType(1, i);
+       medGeometryElement faceType = cellModel->getConstituentType(1, i);
        CELLMODEL faceModel(faceType);
        
        assert(faceModel.getDimension() == 2);
@@ -259,10 +263,10 @@ namespace MEDMEM
        for(int j = 1; j <= faceModel.getNumberOfNodes(); ++j)
          {
            // offset of node from cellIdx
-           int localNodeNumber = cellModel.getNodeConstituent(1, i, j);
+           int localNodeNumber = cellModel->getNodeConstituent(1, i, j);
 
            assert(localNodeNumber >= 1);
-           assert(localNodeNumber <= cellModel.getNumberOfNodes());
+           assert(localNodeNumber <= cellModel->getNumberOfNodes());
 
            const double* node = getCoordsOfNode(localNodeNumber, element, _srcMesh);
            
@@ -312,7 +316,19 @@ namespace MEDMEM
   }
     
 #endif 
-
+CELLMODEL* retrieveCellModel(const medGeometryElement type)
+{
+       map<medGeometryElement, CELLMODEL*>::iterator it = _cellmodel_map.find(type);
+       if (it!=_cellmodel_map.end())
+               return iter->second;
+       else
+       {
+               CELLMODEL* cellmodel=new CELLMODEL(type);
+               _cellmodel_map.insert(make_pair(type,cellmodel));
+               return cellmodel;
+       }
+}
+}
 
  
 };

diff --git a/src/INTERP_KERNEL/Intersector3D.hxx b/src/INTERP_KERNEL/Intersector3D.hxx
index dc8db25c6f2ff665615376920ac925c114f3318f..3dafafaac6992dde61d2d42ae315a80323feecad 100644 (file)
--- a/src/INTERP_KERNEL/Intersector3D.hxx
+++ b/src/INTERP_KERNEL/Intersector3D.hxx
@@ -42,7 +42,8 @@ namespace MEDMEM
 
     const MESH& _srcMesh;
     const MESH& _targetMesh;
-    
+    map<MED_EN::medGeometryElement,MEDMEM::CELLMODEL*> _cellmodel_map;
+    double _dp_accuracy_factor;
   };
 
 //   inline bool Intersector3D::isTriangleOutsideTetra(double* points) const

diff --git a/src/INTERP_KERNEL/MeshUtils.hxx b/src/INTERP_KERNEL/MeshUtils.hxx
index cf9d2e9b08850f3a7deebbe1e987e906bbf8dc68..6b84fdbdb9dc2f0ea33a1bee1c57b8af77568d02 100644 (file)
--- a/src/INTERP_KERNEL/MeshUtils.hxx
+++ b/src/INTERP_KERNEL/MeshUtils.hxx
@@ -53,9 +53,10 @@ namespace INTERP_UTILS
    */
   inline int getNumberOfNodesForType(const medGeometryElement type)
   {
+#if DEBUG
     assert(type > 300);
     assert(type < 400);
-    
+#endif    
     // int(type) = yxx where y is dimension of entity (here 3)
     // and xx is the number of nodes of the element
     return static_cast<int>(type) - 300;

diff --git a/src/INTERP_KERNEL/TransformedTriangle_math.cxx b/src/INTERP_KERNEL/TransformedTriangle_math.cxx
index ebbf768f39e970d45c78e7ab1f65309eb644ec14..c9928202ed7004f33a81cf937cc872bec1d2d7eb 100644 (file)
--- a/src/INTERP_KERNEL/TransformedTriangle_math.cxx
+++ b/src/INTERP_KERNEL/TransformedTriangle_math.cxx
@@ -49,7 +49,7 @@ namespace INTERP_UTILS
   const long double TransformedTriangle::MULT_PREC_F = 4.0 * TransformedTriangle::MACH_EPS;
 
   /// Threshold for resetting double and triple products to zero; ( F / f in Grandy )
-  const long double TransformedTriangle::THRESHOLD_F = 20.0;
+  const long double TransformedTriangle::THRESHOLD_F = 500.0;
 
   /// Threshold for what is considered a small enough angle to warrant correction of triple products by Grandy, [57]
   const double TransformedTriangle::TRIPLE_PRODUCT_ANGLE_THRESHOLD = 0.1;
@@ -374,7 +374,7 @@ namespace INTERP_UTILS
     
     //? is this more stable? -> no subtraction
     //    return asin( dotProd / ( lenNormal * lenEdgeVec ) ) + 3.141592625358979 / 2.0;
-    return 3.14159262535 - acos( dotProd / ( lenNormal * lenEdgeVec ) );
+    return atan(1.0)*4.0 - acos( dotProd / ( lenNormal * lenEdgeVec ) );
 
   }