First try to move assembling from DimMat to dense2D taken from the very cool mtl.

AMDiS/src/AMDiS_fwd.h

@@ -47,6 +47,7 @@ namespace AMDiS {
   class ITL_BasePreconditioner;
   class LeafDataPeriodic;
   class LevelAdmin;
+  class MacroElement;
   class Marker;
   class Mesh; 
   class OEMSolver;
@@ -54,6 +55,7 @@ namespace AMDiS {
   class ProblemInstat;
   class ProblemIterationInterface;
   class ProblemVec;
+  class Projection;
   class PreconditionerScal;
   class Quadrature; 
   class RCNeighbourList;

AMDiS/src/DualTraverse.h

@@ -25,12 +25,10 @@
 #include "Traverse.h"
 #include "Flag.h"
 #include "MemoryManager.h"
+#include "AMDiS_fwd.h"
 
 namespace AMDiS {
 
-  class Mesh;
-  class ElInfo;
-
   /// Parallel traversal of two meshes. 
   class DualTraverse
   {
@@ -56,9 +54,7 @@ namespace AMDiS {
 		       ElInfo **elInfoSmall,
 		       ElInfo **elInfoLarge);
 
-    /** \brief
-     * Get next ElInfo combination
-     */
+    /// Get next ElInfo combination
     bool traverseNext(ElInfo **elInfoNext1,
 		      ElInfo **elInfoNext2,
 		      ElInfo **elInfoSmall,

AMDiS/src/ElInfo.cc

@@ -29,22 +29,21 @@ namespace AMDiS {
       neighbourCoord_(mesh_->getDim(), NO_INIT),
       oppVertex_(mesh_->getDim(), NO_INIT),
       grdLambda(mesh_->getDim(), NO_INIT),
-      subElemCoordsMat(NULL)
+      subElemCoordsMat(NULL),
+      subElemCoordsMat_mtl(2, 2)
   {
     projection_.set(NULL);
 
-    for (int i = 0; i < neighbourCoord_.getSize(); i++) {
+    for (int i = 0; i < neighbourCoord_.getSize(); i++)
       neighbourCoord_[i].init(mesh_->getDim());
-    }
 
     dimOfWorld = Global::getGeo(WORLD);
   }
 
   ElInfo::~ElInfo()
   {
-    if (subElemCoordsMat) {
+    if (subElemCoordsMat)
       DELETE subElemCoordsMat;
-    }
   }
 
 

AMDiS/src/ElInfo.h

@@ -22,21 +22,17 @@
 #ifndef AMDIS_ELINFO_H
 #define AMDIS_ELINFO_H
 
+#include <boost/numeric/mtl/mtl.hpp>
+
 #include "Flag.h"
 #include "Boundary.h"
 #include "Global.h"
 #include "FixVec.h"
 #include "Element.h"
+#include "AMDiS_fwd.h" 
 
 namespace AMDiS {
 
-  class MacroElement;
-  class Mesh;
-  class Element;
-  class BasisFunction;
-  class Projection;
-  template<typename ReturnType, typename ArgumentType> class AbstractFunction;
-
   /** \ingroup Traverse
    * \brief 
    * An ElInfo object holds informations wich are not stored in the corresponding
@@ -217,6 +213,10 @@ namespace AMDiS {
       return subElemCoordsMat;
     }
 
+    inline mtl::dense2D<double>& getSubElemCoordsMat_mtl4() {
+      return subElemCoordsMat_mtl;
+    }
+
     /** \} */ 
 
     /** \name setting methods
@@ -379,10 +379,18 @@ namespace AMDiS {
     virtual void getRefSimplexCoords(const BasisFunction *basisFcts,
 				     DimMat<double> *coords) const = 0;
 
+    virtual void getRefSimplexCoords(const BasisFunction *basisFcts,
+				     mtl::dense2D<double>& coords) const = 0;
+
     virtual void getSubElementCoords(const BasisFunction *basisFcts,
 				     int iChild,
 				     DimMat<double> *coords) const = 0;
 
+    virtual void getSubElementCoords(const BasisFunction *basisFcts,
+				     int iChild,
+				     mtl::dense2D<double>& coords) const = 0;
+
+
   protected:
     /// Pointer to the current mesh
     Mesh *mesh_;
@@ -483,6 +491,8 @@ namespace AMDiS {
      */
     DimMat<double> *subElemCoordsMat;
 
+    mtl::dense2D<double> subElemCoordsMat_mtl;
+
   public:
     /** \brief 
      * child_vertex[el_type][child][i] = father's local vertex index of new 

AMDiS/src/ElInfo2d.h

@@ -22,6 +22,8 @@
 #ifndef AMDIS_ELINFO2D_H
 #define AMDIS_ELINFO2D_H
 
+#include <boost/numeric/mtl/mtl.hpp>
+
 #include "ElInfo.h"
 #include "MemoryManager.h"
 
@@ -63,10 +65,18 @@ namespace AMDiS {
     void getRefSimplexCoords(const BasisFunction *basisFcts,
 			     DimMat<double> *coords) const;
 
+    void getRefSimplexCoords(const BasisFunction *basisFcts,
+			     mtl::dense2D<double>& coords) const {}
+
     void getSubElementCoords(const BasisFunction *basisFcts,
 			     int iChild,
 			     DimMat<double> *coords) const;
 
+    void getSubElementCoords(const BasisFunction *basisFcts,
+			     int iChild,
+			     mtl::dense2D<double>& coords) const {}
+
+
   protected:
     /// Temp vectors for function \ref calcGrdLambda.
     WorldVector<double> *e1, *e2, *normal;

AMDiS/src/ElInfo3d.h

@@ -22,9 +22,10 @@
 #ifndef AMDIS_ELINFO3D_H
 #define AMDIS_ELINFO3D_H
 
+#include <boost/numeric/mtl/mtl.hpp>
+
 #include "ElInfo.h"
 #include "MemoryManager.h"
-#include "OpenMP.h"
 
 namespace AMDiS {
 
@@ -93,10 +94,17 @@ namespace AMDiS {
     void getRefSimplexCoords(const BasisFunction *basisFcts,
 			     DimMat<double> *coords) const;
 
+    void getRefSimplexCoords(const BasisFunction *basisFcts,
+			     mtl::dense2D<double>& coords) const {}
+
     void getSubElementCoords(const BasisFunction *basisFcts,
 			     int iChild, 
 			     DimMat<double> *coords) const;
 
+    void getSubElementCoords(const BasisFunction *basisFcts,
+			     int iChild,
+			     mtl::dense2D<double>& coords) const {}
+
   protected:
     /// \ref el 's type. Is Filled automatically by the traversal routines.
     unsigned char el_type;

AMDiS/src/Mesh.cc

@@ -649,7 +649,7 @@ namespace AMDiS {
     default:
       ERROR_EXIT("invalid dim\n");
       return NULL;
-    };
+    }
   }
 
   bool Mesh::findElInfoAtPoint(const WorldVector<double>& xy,
@@ -690,9 +690,8 @@ namespace AMDiS {
 
     /* now, descend in tree to find leaf element at point */
     bool inside = findElementAtPointRecursive(mel_info, lambda, k, el_info);
-    for (int i = 0; i <= dim; i++) {
-      bary[i] = final_lambda[i];
-    }
+    for (int i = 0; i <= dim; i++)
+      bary[i] = final_lambda[i];   
   
     DELETE mel_info;
 
@@ -716,8 +715,6 @@ namespace AMDiS {
     return val;
   }
 
-
-
   bool Mesh::findElementAtPointRecursive(ElInfo *el_info,
 					 const DimVec<double>& lambda,
 					 int outside,
@@ -940,18 +937,16 @@ namespace AMDiS {
     node.serialize(out);
 
     // write admins
-    int i, size = static_cast<int>(admin.size());
+    int size = static_cast<int>(admin.size());
     out.write(reinterpret_cast<const char*>(&size), sizeof(int));
-    for (i = 0; i < size; i++) {
+    for (int i = 0; i < size; i++)
       admin[i]->serialize(out);
-    }
 
     // write macroElements
     size = static_cast<int>(macroElements.size());
     out.write(reinterpret_cast<const char*>(&size), sizeof(int));
-    for (i = 0; i < size; i++) {
+    for (int i = 0; i < size; i++)
       macroElements[i]->serialize(out);
-    }
 
     // write elementIndex
     out.write(reinterpret_cast<const char*>(&elementIndex), sizeof(int));
@@ -1020,9 +1015,9 @@ namespace AMDiS {
     in.read(reinterpret_cast<char*>(&size), sizeof(int));
     admin.resize(size, NULL);
     for (int i = 0; i < size; i++) {
-      if (!admin[i]) {
+      if (!admin[i])
 	admin[i] = NEW DOFAdmin(this);
-      }
+
       admin[i]->deserialize(in);
     }
 
@@ -1032,10 +1027,10 @@ namespace AMDiS {
     std::vector< std::vector<int> > neighbourIndices(size);
 
     for (int i = 0; i < static_cast<int>(macroElements.size()); i++) {
-      if (macroElements[i]) {
+      if (macroElements[i])
 	DELETE macroElements[i];
-      }
     }
+
     macroElements.resize(size);
     for (int i = 0; i < size; i++) {
       macroElements[i] = NEW MacroElement(dim);
@@ -1151,9 +1146,8 @@ namespace AMDiS {
       result += admin[i]->getUsedSize() * sizeof(DegreeOfFreedom);
     }
     
-    for (int i = 0; i < static_cast<int>(macroElements.size()); i++) {
-      result += macroElements[i]->calcMemoryUsage();
-    }
+    for (int i = 0; i < static_cast<int>(macroElements.size()); i++)
+      result += macroElements[i]->calcMemoryUsage();    
     
     return result;
   }

AMDiS/src/RefinementManager2d.cc

@@ -151,9 +151,8 @@ namespace AMDiS {
 
     Projection *projector = el_info->getProjection(0);
 
-    if (!projector || projector->getType() != VOLUME_PROJECTION) {
-      projector = el_info->getProjection(2);
-    }
+    if (!projector || projector->getType() != VOLUME_PROJECTION)
+      projector = el_info->getProjection(2);   
 
     if (el->getFirstChild() && projector && (!el->isNewCoordSet())) {
       WorldVector<double> *new_coord = NEW WorldVector<double>;
@@ -229,9 +228,8 @@ namespace AMDiS {
     /*  if there are functions to interpolate data to the finer grid, do so     */
     /****************************************************************************/
   
-    int iadmin;
     int nrAdmin = mesh->getNumberOfDOFAdmin();
-    for(iadmin = 0; iadmin < nrAdmin; iadmin++) {
+    for(int iadmin = 0; iadmin < nrAdmin; iadmin++) {
       std::list<DOFIndexedBase*>::iterator it;
       DOFAdmin* admin = const_cast<DOFAdmin*>(&mesh->getDOFAdmin(iadmin));
       std::list<DOFIndexedBase*>::iterator end = admin->endDOFIndexed();
@@ -358,15 +356,12 @@ namespace AMDiS {
   {
     FUNCNAME("RefinementManager2d::getRefinePatch()");
 
-    Triangle *el = dynamic_cast<Triangle*>(const_cast<Element*>( (*el_info)->getElement()));
-    int opp_vertex = 0;
-
     if ((*el_info)->getNeighbour(2) && (*el_info)->getOppVertex(2) != 2) {
       /****************************************************************************/
       /*  neighbour is not compatible devisible; refine neighbour first, store the*/
       /*  opp_vertex to traverse back to el                                       */
       /****************************************************************************/
-      opp_vertex = (*el_info)->getOppVertex(2);
+      int opp_vertex = (*el_info)->getOppVertex(2);
       
       ElInfo *neigh_info = stack->traverseNeighbour2d(*el_info, 2);
       neigh_info->getElement()->setMark(max(neigh_info->getElement()->getMark(), 1));
@@ -376,7 +371,9 @@ namespace AMDiS {
       /*  now go back to the original element and refine the patch                */
       /****************************************************************************/
       *el_info = neigh_info = stack->traverseNeighbour2d(neigh_info, opp_vertex);
-      TEST_EXIT_DBG(neigh_info->getElement() == el)("invalid traverse_neighbour1\n");
+      TEST_EXIT_DBG(neigh_info->getElement() == 
+		    dynamic_cast<Triangle*>(const_cast<Element*>((*el_info)->getElement())))
+	("invalid traverse_neighbour1\n");
     }
   
     if (refineList->setElement(1, (*el_info)->getNeighbour(2))) {
@@ -385,7 +382,7 @@ namespace AMDiS {
       *n_neigh = 2;
     }
 
-    return(0);
+    return 0;
   }
 
 }

AMDiS/src/Traverse.cc

@@ -1064,9 +1064,8 @@ namespace AMDiS {
 
     TEST_EXIT_DBG(traverse_mesh->getDim() == 3)("update only in 3d\n");
 
-    for (int i = stack_used; i > 0; i--) {
+    for (int i = stack_used; i > 0; i--)
       dynamic_cast<ElInfo3d*>(elinfo_stack[i])->update();
-    }
   }
 
   void TraverseStack::getCoordsInElem(const ElInfo *upperElInfo, 
@@ -1077,11 +1076,24 @@ namespace AMDiS {
 
     upperElInfo->getRefSimplexCoords(basisFcts, coords);
     
-    for (int i = 1; i <= levelDif; i++) {
+    for (int i = 1; i <= levelDif; i++)
+      upperElInfo->getSubElementCoords(basisFcts, 
+				       elinfo_stack[stack_used - levelDif + i]->getIChild(),
+				       coords);
+  }
+
+  void TraverseStack::getCoordsInElem(const ElInfo *upperElInfo, 
+				      const BasisFunction *basisFcts,
+				      mtl::dense2D<double>& coords)
+  {
+    int levelDif = elinfo_stack[stack_used]->getLevel() - upperElInfo->getLevel();
+
+    upperElInfo->getRefSimplexCoords(basisFcts, coords);
+    
+    for (int i = 1; i <= levelDif; i++)
       upperElInfo->getSubElementCoords(basisFcts, 
 				       elinfo_stack[stack_used - levelDif + i]->getIChild(),
 				       coords);
-    }
   }
 
 }

AMDiS/src/Traverse.h

@@ -29,21 +29,20 @@
 #ifndef AMDIS_TRAVERSE_H
 #define AMDIS_TRAVERSE_H
 
+#include <vector>
+#include <deque>
+#include <stack>
+
 #include "Flag.h"
 #include "Global.h"
 #include "ElInfo.h"
 #include "ElInfoStack.h"
 #include "MemoryManager.h"
 #include "OpenMP.h"
-#include <vector>
-#include <deque>
-#include <stack>
+#include "AMDiS_fwd.h"
 
 namespace AMDiS {
 
-  class MacroElement;
-  class Mesh;
-
   /** \ingroup Traverse
    * \brief
    * Mesh refinement and coarsening routines are examples of functions which 
@@ -59,9 +58,7 @@ namespace AMDiS {
   public:
     MEMORY_MANAGED(TraverseStack);
 
-    /** \brief
-     * Creates an empty TraverseStack
-     */
+    /// Creates an empty TraverseStack
     TraverseStack() 
       : traverse_mel(NULL),
         stack_size(0),
@@ -69,20 +66,16 @@ namespace AMDiS {
         save_stack_used(0),
         myThreadId_(0),
         maxThreads_(1)
-    {
-    }
+    {}
 
-    /** \brief
-     * Destructor
-     */
+    /// Destructor
     ~TraverseStack() 
     {
-      for (int i = 0; i < static_cast<int>(elinfo_stack.size()); i++) {
+      for (int i = 0; i < static_cast<int>(elinfo_stack.size()); i++)
 	DELETE elinfo_stack[i];
-      }
-      for (int i = 0; i < static_cast<int>(save_elinfo_stack.size()); i++) {
+
+      for (int i = 0; i < static_cast<int>(save_elinfo_stack.size()); i++) 
 	DELETE save_elinfo_stack[i];
-      }
     }
 
   public:
@@ -127,6 +120,9 @@ namespace AMDiS {
     void getCoordsInElem(const ElInfo *upperElInfo, const BasisFunction *basisFcts,
 			 DimMat<double> *coords);
 
+    void getCoordsInElem(const ElInfo *upperElInfo, const BasisFunction *basisFcts,
+			 mtl::dense2D<double>& coords);
+
     /// Is used for parallel mesh traverse.
     inline void setMyThreadId(int myThreadId) {
       myThreadId_ = myThreadId;