diff --git a/AMDiS/src/BasisFunction.h b/AMDiS/src/BasisFunction.h
index a4864bb1e040882627fb60ada4877f4a8ca81ed0..55f681632e7980e80d21342ee6ed3860a3f0e35c 100644
--- a/AMDiS/src/BasisFunction.h
+++ b/AMDiS/src/BasisFunction.h
@@ -181,7 +181,7 @@ namespace AMDiS {
virtual DimVec<double> *getCoords(int i) const = 0;
/** \brief
- * Returns a pointer to a const vector with interpolation coefficients of the
+ * Fills a vector with interpolation coefficients of the
* function f; if indices is a pointer to NULL, the coefficient for all
* basis functions are calculated and the i-th entry in the vector is the
* coefficient of the i-th basis function; if indices is non NULL, only the
@@ -197,16 +197,18 @@ namespace AMDiS {
* during mesh traversal.
* Must be implemented by sub classes.
*/
- virtual const double* interpol(const ElInfo *el_info,
- int n, const int *indices,
- AbstractFunction<double, WorldVector<double> > *f,
- double *coeff) = 0;
+ virtual void interpol(const ElInfo *el_info,
+ int n,
+ const int *indices,
+ AbstractFunction<double, WorldVector<double> > *f,
+ mtl::dense_vector<double> &coeff) const = 0;
/// WorldVector<double> valued interpol function.
- virtual const WorldVector<double>* interpol(const ElInfo *el_info, int no,
- const int *b_no,
- AbstractFunction<WorldVector<double>, WorldVector<double> > *f,
- WorldVector<double> *vec) = 0;
+ virtual void interpol(const ElInfo *el_info,
+ int no,
+ const int *b_no,
+ AbstractFunction<WorldVector<double>, WorldVector<double> > *f,
+ mtl::dense_vector<WorldVector<double> >& coeff) const = 0;
/// Returns the i-th local basis function
inline BasFctType *getPhi(int i) const
diff --git a/AMDiS/src/DOFVector.h b/AMDiS/src/DOFVector.h
index 56e219dbf34ce927399b4bc99f49bc560c549606..7e91601236c775ce5f5796c1d0d9fc9a08575668 100644
--- a/AMDiS/src/DOFVector.h
+++ b/AMDiS/src/DOFVector.h
@@ -572,9 +572,6 @@ namespace AMDiS {
/// Returns the average value of the DOFVector.
T average() const;
- /// Used by interpol while mesh traversal.
- void interpolFct(ElInfo* elinfo);
-
/// Prints \ref vec to stdout.
void print() const;
@@ -641,12 +638,6 @@ namespace AMDiS {
/// Specifies what operation should be performed after coarsening
CoarsenOperation coarsenOperation;
-
- /// Used by \ref interpol
- AbstractFunction<T, WorldVector<double> > *interFct;
-
- /// Used for mesh traversal
- static DOFVector<T> *traverseVector;
};
@@ -795,7 +786,7 @@ namespace AMDiS {
// y = a*x + y
template<typename T>
- void axpy(double a,const DOFVector<T>& x, DOFVector<T>& y);
+ void axpy(double a, const DOFVector<T>& x, DOFVector<T>& y);
// matrix vector product
template<typename T>
@@ -856,6 +847,12 @@ namespace AMDiS {
return vec->getUsedSize();
}
+ template<typename T>
+ inline int size(const DOFVector<T>& vec)
+ {
+ return vec.getUsedSize();
+ }
+
template<typename T>
inline void checkFeSpace(const FiniteElemSpace* feSpace, const std::vector<T>& vec)
{
diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh
index 7897d4d5a89ee1afc51a99871f7625805dda7e93..e749f3075b86b2360639db8367a79eee92b245b7 100644
--- a/AMDiS/src/DOFVector.hh
+++ b/AMDiS/src/DOFVector.hh
@@ -132,9 +132,6 @@ namespace AMDiS {
vec.clear();
}
- template<typename T>
- DOFVector<T> * DOFVector<T>::traverseVector = NULL;
-
template<typename T>
void DOFVectorBase<T>::addElementVector(T factor,
const ElementVector &elVec,
@@ -450,8 +447,6 @@ namespace AMDiS {
TEST_EXIT_DBG(fct)("No function to interpolate!\n");
- interFct = fct;
-
if (!this->getFeSpace()) {
MSG("no dof admin in vec %s, skipping interpolation\n",
this->getName().c_str());
@@ -476,36 +471,27 @@ namespace AMDiS {
return;
}
- traverseVector = this;
+ const BasisFunction *basFct = this->getFeSpace()->getBasisFcts();
+ const DOFAdmin* admin = this->getFeSpace()->getAdmin();
+ int nBasFcts = basFct->getNumber();
+ std::vector<DegreeOfFreedom> myLocalIndices(nBasFcts);
+ mtl::dense_vector<double> fctInterpolValues(nBasFcts);
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(this->getFeSpace()->getMesh(), -1,
Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
while (elInfo) {
- interpolFct(elInfo);
+ basFct->interpol(elInfo, 0, NULL, fct, fctInterpolValues);
+ basFct->getLocalIndices(const_cast<Element*>(elInfo->getElement()),
+ admin, myLocalIndices);
+ for (int i = 0; i < nBasFcts; i++)
+ vec[myLocalIndices[i]] = fctInterpolValues[i];
+
elInfo = stack.traverseNext(elInfo);
}
}
- template<typename T>
- void DOFVector<T>::interpolFct(ElInfo* elinfo)
- {
- const BasisFunction *basFct = traverseVector->getFeSpace()->getBasisFcts();
- const DOFAdmin* admin = traverseVector->getFeSpace()->getAdmin();
- const T *inter_val =
- const_cast<BasisFunction*>(basFct)->interpol(elinfo, 0, NULL,
- traverseVector->interFct, NULL);
-
- int nBasFcts = basFct->getNumber();
- std::vector<DegreeOfFreedom> myLocalIndices(nBasFcts);
- basFct->getLocalIndices(const_cast<Element*>(elinfo->getElement()),
- admin, myLocalIndices);
- for (int i = 0; i < nBasFcts; i++)
- (*traverseVector)[myLocalIndices[i]] = inter_val[i];
- }
-
-
template<typename T>
double DOFVector<T>::Int(Quadrature* q) const
{
@@ -1047,7 +1033,6 @@ namespace AMDiS {
this->nBasFcts = rhs.nBasFcts;
vec = rhs.vec;
this->elementVector.change_dim(this->nBasFcts);
- interFct = rhs.interFct;
coarsenOperation = rhs.coarsenOperation;
this->operators = rhs.operators;
this->operatorFactor = rhs.operatorFactor;
diff --git a/AMDiS/src/Lagrange.cc b/AMDiS/src/Lagrange.cc
index 06859ca1e5ded1493143d95c4cf861c8460d97e2..49d309125c98ad7098eb3cc2bc21e78d8e53e655 100644
--- a/AMDiS/src/Lagrange.cc
+++ b/AMDiS/src/Lagrange.cc
@@ -811,43 +811,21 @@ namespace AMDiS {
}
- const double* Lagrange::interpol(const ElInfo *el_info,
- int no, const int *b_no,
- AbstractFunction<double, WorldVector<double> > *f,
- double *vec)
+ void Lagrange::interpol(const ElInfo *el_info,
+ int no,
+ const int *b_no,
+ AbstractFunction<double, WorldVector<double> > *f,
+ mtl::dense_vector<double> &rvec) const
{
FUNCNAME("Lagrange::interpol()");
- static double* localVec = NULL;
- static int localVecSize = 0;
-
- double *rvec = NULL;
-
- if (vec) {
- rvec = vec;
- } else {
- if (localVec && nBasFcts > localVecSize) {
- delete [] localVec;
- localVec = new double[nBasFcts];
- }
- if (!localVec)
- localVec = new double[nBasFcts];
-
- localVecSize = nBasFcts;
- rvec = localVec;
- }
-
WorldVector<double> x;
el_info->testFlag(Mesh::FILL_COORDS);
if (b_no) {
- if (no <= 0 || no > getNumber()) {
- ERROR("something is wrong, doing nothing\n");
- rvec[0] = 0.0;
- return(const_cast<const double *>(rvec));
- }
-
+ TEST_EXIT_DBG(no >= 0 && no < getNumber())("Something is wrong!\n");
+
for (int i = 0; i < no; i++) {
if (b_no[i] < Global::getGeo(VERTEX, dim)) {
rvec[i] = (*f)(el_info->getCoord(b_no[i]));
@@ -865,23 +843,17 @@ namespace AMDiS {
rvec[i] = (*f)(x);
}
}
-
- return(const_cast<const double *>( rvec));
}
- const WorldVector<double>*
- Lagrange::interpol(const ElInfo *el_info,
- int no,
- const int *b_no,
- AbstractFunction<WorldVector<double>, WorldVector<double> > *f,
- WorldVector<double> *vec)
+ void Lagrange::interpol(const ElInfo *el_info,
+ int no,
+ const int *b_no,
+ AbstractFunction<WorldVector<double>, WorldVector<double> > *f,
+ mtl::dense_vector<WorldVector<double> > &rvec) const
{
- FUNCNAME("*Lagrange::interpol_d()");
+ FUNCNAME("*Lagrange::interpol()");
- static WorldVector<double> *inter;
- WorldVector<double> *rvec =
- vec ? vec : (inter ? inter : inter = new WorldVector<double>[getNumber()]);
WorldVector<double> x;
el_info->testFlag(Mesh::FILL_COORDS);
@@ -889,12 +861,8 @@ namespace AMDiS {
int vertices = Global::getGeo(VERTEX, dim);
if (b_no) {
- if (no <= 0 || no > getNumber()) {
- ERROR("something is wrong, doing nothing\n");
- for (int i = 0; i < dow; i++)
- rvec[0][i] = 0.0;
- return(const_cast<const WorldVector<double> *>( rvec));
- }
+ TEST_EXIT_DBG(no >= 0 && no < getNumber())("Something is wrong!\n");
+
for (int i = 0; i < no; i++) {
if (b_no[i] < Global::getGeo(VERTEX, dim)) {
rvec[i] = (*f)(el_info->getCoord(b_no[i]));
@@ -911,8 +879,6 @@ namespace AMDiS {
rvec[i] = (*f)(x);
}
}
-
- return(const_cast<const WorldVector<double> *>( rvec));
}
diff --git a/AMDiS/src/Lagrange.h b/AMDiS/src/Lagrange.h
index 79c5e38b4679b375087117243aa9981fc21f41fe..2e0dbf79031aea2aefda2bd53c2725ae62b2cce7 100644
--- a/AMDiS/src/Lagrange.h
+++ b/AMDiS/src/Lagrange.h
@@ -64,16 +64,15 @@ namespace AMDiS {
static Lagrange* getLagrange(int dim, int degree);
/// Implements BasisFunction::interpol
- const double *interpol(const ElInfo *, int, const int *,
- AbstractFunction<double, WorldVector<double> >*,
- double *);
+ void interpol(const ElInfo *, int, const int *,
+ AbstractFunction<double, WorldVector<double> >*,
+ mtl::dense_vector<double>&) const;
/// Implements BasisFunction::interpol
- const WorldVector<double> *interpol(const ElInfo *, int,
- const int *b_no,
- AbstractFunction<WorldVector<double>,
- WorldVector<double> >*,
- WorldVector<double> *);
+ void interpol(const ElInfo *, int,
+ const int *b_no,
+ AbstractFunction<WorldVector<double>, WorldVector<double> >*,
+ mtl::dense_vector<WorldVector<double> >&) const;
/// Returns the barycentric coordinates of the i-th basis function.
DimVec<double> *getCoords(int i) const;
diff --git a/AMDiS/src/ProblemInstat.h b/AMDiS/src/ProblemInstat.h
index 9528d0bdbf42ebb73849e76fc988ccb6600c1d73..6d4ff02a6859aff8d4db8a70abfc10940fdfd979 100644
--- a/AMDiS/src/ProblemInstat.h
+++ b/AMDiS/src/ProblemInstat.h
@@ -164,13 +164,13 @@ namespace AMDiS {
Flag adoptFlag = INIT_NOTHING);
/// Used in \ref initialize().
- void createUhOld();
+ virtual void createUhOld();
/// Implementation of ProblemTimeInterface::initTimestep().
- void initTimestep(AdaptInfo *adaptInfo);
+ virtual void initTimestep(AdaptInfo *adaptInfo);
/// Implementation of ProblemTimeInterface::closeTimestep().
- void closeTimestep(AdaptInfo *adaptInfo);
+ virtual void closeTimestep(AdaptInfo *adaptInfo);
/// Returns \ref problemStat.
inline ProblemStatSeq* getStatProblem()
@@ -190,11 +190,11 @@ namespace AMDiS {
}
/// Used by \ref problemInitial
- void transferInitialSolution(AdaptInfo *adaptInfo);
+ virtual void transferInitialSolution(AdaptInfo *adaptInfo);
- void serialize(std::ostream &out) {}
+ virtual void serialize(std::ostream &out) {}
- void deserialize(std::istream &in) {}
+ virtual void deserialize(std::istream &in) {}
protected:
/// Space problem solved in each timestep.
diff --git a/AMDiS/src/ProblemStat.cc b/AMDiS/src/ProblemStat.cc
index dfbb4a2dead232a1b4a39f19345d8d07b067287b..255e057e4ddba01b4d78d5693d3a498762dc963a 100644
--- a/AMDiS/src/ProblemStat.cc
+++ b/AMDiS/src/ProblemStat.cc
@@ -1065,7 +1065,7 @@ namespace AMDiS {
void ProblemStatSeq::dualAssemble(AdaptInfo *adaptInfo, Flag flag,
- bool asmMatrix, bool asmVector)
+ bool asmMatrix, bool asmVector)
{
FUNCNAME("ProblemStat::dualAssemble()");
@@ -1631,7 +1631,8 @@ namespace AMDiS {
void ProblemStatSeq::assembleOnOneMesh(FiniteElemSpace *feSpace,
Flag assembleFlag,
- DOFMatrix *matrix, DOFVector<double> *vector)
+ DOFMatrix *matrix,
+ DOFVector<double> *vector)
{
FUNCNAME("ProblemStat::assembleOnOnMesh()");