From b538f51c6cb192a3c4af3a715e76505dba1ccc0e Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Thu, 7 May 2009 06:21:02 +0000
Subject: [PATCH] First try to move assembling from DimMat to dense2D taken
from the very cool mtl.
---
AMDiS/src/AMDiS_fwd.h | 2 ++
AMDiS/src/DualTraverse.h | 8 ++------
AMDiS/src/ElInfo.cc | 9 ++++-----
AMDiS/src/ElInfo.h | 24 +++++++++++++++++-------
AMDiS/src/ElInfo2d.h | 10 ++++++++++
AMDiS/src/ElInfo3d.h | 10 +++++++++-
AMDiS/src/Mesh.cc | 30 ++++++++++++------------------
AMDiS/src/RefinementManager2d.cc | 19 ++++++++-----------
AMDiS/src/Traverse.cc | 20 ++++++++++++++++----
AMDiS/src/Traverse.h | 32 ++++++++++++++------------------
10 files changed, 94 insertions(+), 70 deletions(-)
diff --git a/AMDiS/src/AMDiS_fwd.h b/AMDiS/src/AMDiS_fwd.h
index ae8b776b..dbbadad6 100644
--- a/AMDiS/src/AMDiS_fwd.h
+++ b/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;
diff --git a/AMDiS/src/DualTraverse.h b/AMDiS/src/DualTraverse.h
index 91a93e55..920436a0 100644
--- a/AMDiS/src/DualTraverse.h
+++ b/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,
diff --git a/AMDiS/src/ElInfo.cc b/AMDiS/src/ElInfo.cc
index 84d14023..9bbbcac9 100644
--- a/AMDiS/src/ElInfo.cc
+++ b/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;
- }
}
diff --git a/AMDiS/src/ElInfo.h b/AMDiS/src/ElInfo.h
index 9841046b..ca786c39 100644
--- a/AMDiS/src/ElInfo.h
+++ b/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
diff --git a/AMDiS/src/ElInfo2d.h b/AMDiS/src/ElInfo2d.h
index e2ad7c9d..5e5cf616 100644
--- a/AMDiS/src/ElInfo2d.h
+++ b/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;
diff --git a/AMDiS/src/ElInfo3d.h b/AMDiS/src/ElInfo3d.h
index 48b92fc0..47c935bf 100644
--- a/AMDiS/src/ElInfo3d.h
+++ b/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;
diff --git a/AMDiS/src/Mesh.cc b/AMDiS/src/Mesh.cc
index 47d38e13..c0e2ce50 100644
--- a/AMDiS/src/Mesh.cc
+++ b/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;
}
diff --git a/AMDiS/src/RefinementManager2d.cc b/AMDiS/src/RefinementManager2d.cc
index 0699e562..ef9681c6 100644
--- a/AMDiS/src/RefinementManager2d.cc
+++ b/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;
}
}
diff --git a/AMDiS/src/Traverse.cc b/AMDiS/src/Traverse.cc
index 818f3599..c878898f 100644
--- a/AMDiS/src/Traverse.cc
+++ b/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);
- }
}
}
diff --git a/AMDiS/src/Traverse.h b/AMDiS/src/Traverse.h
index 7f4532a7..65d66ee6 100644
--- a/AMDiS/src/Traverse.h
+++ b/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;
--
GitLab