diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc
index 7abbcce4c1a827970983fed79561192a3dd4a43b..593b86031afc5d705c913fa5b81bdd88bbecf10d 100644
--- a/AMDiS/src/DOFVector.cc
+++ b/AMDiS/src/DOFVector.cc
@@ -887,7 +887,7 @@ namespace AMDiS {
TEST_EXIT(fct)("No function defined!\n");
- int deg = 2 * vec1.getFeSpace()->getBasisFcts()->getDegree();
+ int deg = fct->getDegree();
Quadrature* quad =
Quadrature::provideQuadrature(vec1.getFeSpace()->getMesh()->getDim(), deg);
FastQuadrature *fastQuad =
@@ -925,7 +925,7 @@ namespace AMDiS {
TEST_EXIT(fct)("No function defined!\n");
- int deg = 2 * vec1.getFeSpace()->getBasisFcts()->getDegree();
+ int deg = fct->getDegree();
Quadrature* quad =
Quadrature::provideQuadrature(vec1.getFeSpace()->getMesh()->getDim(), deg);
FastQuadrature *fastQuad =
diff --git a/AMDiS/src/DOFVector.h b/AMDiS/src/DOFVector.h
index 3841085dfa4c26e87107a538840f9a21a1b381b5..39c2d737dfc975cdcd24654708923507b57b5fd3 100644
--- a/AMDiS/src/DOFVector.h
+++ b/AMDiS/src/DOFVector.h
@@ -555,7 +555,7 @@ namespace AMDiS {
*/
void interpol(AbstractFunction<T, WorldVector<double> > *fct);
- void interpol(DOFVector<T> *v, double factor=1.0);
+ void interpol(DOFVector<T> *v, double factor = 1.0);
/// Writes the data of the DOFVector to an output stream.
void serialize(std::ostream &out)
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index 7824e8c208b5035596f5d581745a5876cd7ebce8..03baff9180f780ebab88d81092f576254172da9c 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -1232,6 +1232,21 @@ namespace AMDiS {
}
+ void ProblemVec::addNeumannBC(BoundaryType type, int row, int col,
+ DOFVector<double> *n)
+ {
+ FUNCNAME("ProblemVec::addNeumannBC()");
+
+ boundaryConditionSet = true;
+
+ NeumannBC *neumann =
+ new NeumannBC(type, n, componentSpaces[row], componentSpaces[col]);
+
+ if (rhs)
+ rhs->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(neumann);
+ }
+
+
void ProblemVec::addRobinBC(BoundaryType type, int row, int col,
AbstractFunction<double, WorldVector<double> > *n,
AbstractFunction<double, WorldVector<double> > *r)
diff --git a/AMDiS/src/ProblemVec.h b/AMDiS/src/ProblemVec.h
index bcf7c8e99ee8dece2150e0173fea6110d1bd1b7c..567e7d4dc4323ba716de00c73751b22018bc31f4 100644
--- a/AMDiS/src/ProblemVec.h
+++ b/AMDiS/src/ProblemVec.h
@@ -218,39 +218,48 @@ namespace AMDiS {
/// Interpolates fct to \ref solution.
void interpolInitialSolution(std::vector<AbstractFunction<double, WorldVector<double> >*> *fct);
- /// Adds an Operator to \ref A.
+ /// Adds an operator to \ref A.
void addMatrixOperator(Operator *op, int i, int j,
double *factor = NULL, double *estFactor = NULL);
- /// Adds an Operator to \ref A.
+ /// Adds an operator to \ref A.
void addMatrixOperator(Operator &op, int i, int j,
double *factor = NULL, double *estFactor = NULL);
- /// Adds an Operator to \ref rhs.
+ /// Adds an operator to \ref rhs.
void addVectorOperator(Operator *op, int i,
double *factor = NULL, double *estFactor = NULL);
- /// Adds an Operator to \ref rhs.
+ /// Adds an operator to \ref rhs.
void addVectorOperator(Operator &op, int i,
double *factor = NULL, double *estFactor = NULL);
- /// Adds dirichlet boundary conditions.
+ /// Adds a Dirichlet boundary condition, where the rhs is given by an
+ /// abstract function.
virtual void addDirichletBC(BoundaryType type, int row, int col,
AbstractFunction<double, WorldVector<double> > *b);
+ /// Adds a Dirichlet boundary condition, where the rhs is given by a DOF
+ /// vector.
virtual void addDirichletBC(BoundaryType type, int row, int col,
DOFVector<double> *vec);
- /// Adds neumann boundary conditions.
+ /// Adds a Neumann boundary condition, where the rhs is given by an
+ /// abstract function.
virtual void addNeumannBC(BoundaryType type, int row, int col,
AbstractFunction<double, WorldVector<double> > *n);
- /// Adds robin boundary conditions.
+ /// Adds a Neumann boundary condition, where the rhs is given by an DOF
+ /// vector.
+ virtual void addNeumannBC(BoundaryType type, int row, int col,
+ DOFVector<double> *n);
+
+ /// Adds Robin boundary condition.
virtual void addRobinBC(BoundaryType type, int row, int col,
AbstractFunction<double, WorldVector<double> > *n,
AbstractFunction<double, WorldVector<double> > *r);
- /// Adds periodic boundary conditions.
+ /// Adds a periodic boundary condition.
virtual void addPeriodicBC(BoundaryType type, int row, int col);
/** \brief
diff --git a/AMDiS/src/Quadrature.cc b/AMDiS/src/Quadrature.cc
index 38c4b8a2fdbc439867666049833e1d3d19baade7..e47c5e7cd388bf09b480f9b2202a01fad60da7ad 100644
--- a/AMDiS/src/Quadrature.cc
+++ b/AMDiS/src/Quadrature.cc
@@ -1412,6 +1412,7 @@ namespace AMDiS {
Quadrature::quad_nd[3] = Quadrature::quad_3d;
};
+
Quadrature* Quadrature::provideQuadrature(int dim_, int degree_)
{
FUNCNAME("Quadrature::provideQuadrature()");
@@ -1439,12 +1440,14 @@ namespace AMDiS {
return (quad_nd[dim_][degree_]);
}
+
Quadrature::~Quadrature()
{
delete lambda;
- delete[] w;
+ delete [] w;
}
+
double Quadrature::integrateStdSimplex(AbstractFunction<double, DimVec<double> > *f)
{
FUNCNAME("Quadrature::integrateStdSimplex()");
@@ -1462,15 +1465,15 @@ namespace AMDiS {
return result;
}
+
FastQuadrature* FastQuadrature::provideFastQuadrature(const BasisFunction* bas_fcts,
const Quadrature& quad,
Flag init_flag)
{
- FUNCNAME("FastQuadrature::FastQuadrature()");
+ FUNCNAME("FastQuadrature::provideFastQuuadrature()");
FastQuadrature *quad_fast;
-
#ifdef _OPENMP
#pragma omp critical (provideFastQuad)
#endif
@@ -1482,11 +1485,13 @@ namespace AMDiS {
if ((*fast)->basisFunctions == bas_fcts && (*fast)->quadrature == &quad)
break;
- if ((fast != fastQuadList.end()) && (((*fast)->init_flag & init_flag) == init_flag)) {
+ if (fast != fastQuadList.end() && ((*fast)->init_flag & init_flag) == init_flag) {
quad_fast = *fast;
} else {
if (fast == fastQuadList.end()) {
- quad_fast = new FastQuadrature(const_cast<BasisFunction*>(bas_fcts), const_cast<Quadrature*>(&quad), 0);
+ quad_fast =
+ new FastQuadrature(const_cast<BasisFunction*>(bas_fcts),
+ const_cast<Quadrature*>(&quad), 0);
fastQuadList.push_front(quad_fast);
@@ -1502,6 +1507,7 @@ namespace AMDiS {
return quad_fast;
}
+
void FastQuadrature::init(Flag init_flag)
{
int dim = quadrature->getDim();
@@ -1567,6 +1573,7 @@ namespace AMDiS {
}
}
+
FastQuadrature::FastQuadrature(const FastQuadrature& fastQuad)
{
FUNCNAME("FastQuadrature::FastQuadrature()");
@@ -1609,6 +1616,7 @@ namespace AMDiS {
}
}
+
FastQuadrature::~FastQuadrature()
{
int nPoints = quadrature->getNumPoints();
@@ -1621,11 +1629,13 @@ namespace AMDiS {
delete D2Phi;
}
+
const double FastQuadrature::getSecDer(int q,int i ,int j, int m) const
{
return (D2Phi) ? (*D2Phi)[q][i][j][m] : 0.0;
}
+
const VectorOfFixVecs<DimMat<double> > *FastQuadrature::getSecDer(int q) const
{
return D2Phi ? (&((*D2Phi)[q])) : NULL;
diff --git a/AMDiS/src/Quadrature.h b/AMDiS/src/Quadrature.h
index 9be2028caf02a43f2fe55607ab399f81a3c8048c..aa7037db7009e34ff8e440b040864f3ac5fbca4e 100644
--- a/AMDiS/src/Quadrature.h
+++ b/AMDiS/src/Quadrature.h
@@ -324,9 +324,7 @@ namespace AMDiS {
* Constructs a FastQuadrature for the given Quadrature, BasisFunction, and
* flag.
*/
- FastQuadrature(BasisFunction* basFcts,
- Quadrature* quad,
- Flag flag)
+ FastQuadrature(BasisFunction* basFcts, Quadrature* quad, Flag flag)
: init_flag(flag),
phi(NULL),
grdPhi(NULL),
diff --git a/AMDiS/src/RobinBC.cc b/AMDiS/src/RobinBC.cc
index bbac5f8437760cc0f2b80c2907eba7579e9d4a86..3b5111ef5a20331403989c826fa741937f335df4 100644
--- a/AMDiS/src/RobinBC.cc
+++ b/AMDiS/src/RobinBC.cc
@@ -22,14 +22,13 @@ namespace AMDiS {
// create barycentric coords for each vertex of each side
const Element *refElement = Global::getReferenceElement(dim);
- coords = new VectorOfFixVecs<DimVec<double> >*[dim+1];
+ coords = new VectorOfFixVecs<DimVec<double> >*[dim + 1];
// for all element sides
for (int i = 0; i < dim + 1; i++) {
- coords[i] = new VectorOfFixVecs<DimVec<double> >(dim, dim, DEFAULT_VALUE,
- DimVec<double>(dim,
- DEFAULT_VALUE,
- 0.0));
+ coords[i] =
+ new VectorOfFixVecs<DimVec<double> >(dim, dim, DEFAULT_VALUE,
+ DimVec<double>(dim, DEFAULT_VALUE, 0.0));
// for each vertex of the side
for (int k = 0; k < dim; k++) {
int index = refElement->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim), i, k);
@@ -78,11 +77,9 @@ namespace AMDiS {
// for all element sides
for (int i = 0; i < dim + 1; i++) {
- coords[i] = new VectorOfFixVecs<DimVec<double> >(dim, dim, DEFAULT_VALUE,
- DimVec<double>(dim,
- DEFAULT_VALUE,
- 0.0)
- );
+ coords[i] =
+ new VectorOfFixVecs<DimVec<double> >(dim, dim, DEFAULT_VALUE,
+ DimVec<double>(dim, DEFAULT_VALUE, 0.0));
// for each vertex of the side
for (int k = 0; k < dim; k++) {
int index = refElement->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim), i, k);
@@ -230,7 +227,6 @@ namespace AMDiS {
(*op)->getAssembler(omp_get_thread_num())->initElement(elInfo);
for (int face = 0; face < dim + 1; face++) {
-
elInfo->getNormal(face, normal);
Quadrature *quadrature =
@@ -251,12 +247,8 @@ namespace AMDiS {
std::vector<double> f(nPoints, 0.0);
if (robinOperators)
- (*robinOperators)[face]->evalZeroOrder(nPoints,
- uhAtQp,
- NULL,
- NULL,
- &f[0],
- 1.0);
+ (*robinOperators)[face]->evalZeroOrder(nPoints, uhAtQp,
+ NULL, NULL, &f[0], 1.0);
std::vector<WorldVector<double> > grdUh(nPoints);
std::vector<WorldVector<double> > A_grdUh(nPoints);
diff --git a/AMDiS/src/RobinBC.h b/AMDiS/src/RobinBC.h
index 1ba10e485437c9bd56920c2d298cd1693f231e3a..cc4ec1e16ce02bdb58976dee54e5dbf88285356c 100644
--- a/AMDiS/src/RobinBC.h
+++ b/AMDiS/src/RobinBC.h
@@ -90,6 +90,7 @@ namespace AMDiS {
VectorOfFixVecs<DimVec<double> >**coords;
};
+
class NeumannBC : public RobinBC
{
public:
@@ -99,6 +100,13 @@ namespace AMDiS {
FiniteElemSpace *colFeSpace = NULL)
: RobinBC(type, j, NULL, rowFeSpace, colFeSpace)
{}
+
+ NeumannBC(BoundaryType type,
+ DOFVectorBase<double> *j,
+ FiniteElemSpace *rowFeSpace,
+ FiniteElemSpace *colFeSpace = NULL)
+ : RobinBC(type, j, NULL, rowFeSpace, colFeSpace)
+ {}
};
}
diff --git a/AMDiS/src/time/RosenbrockStationary.cc b/AMDiS/src/time/RosenbrockStationary.cc
index d3bc410091c6411a05bdf840df972aa3b070c866..e1abde9e084a2cde2c76d483f0a529c19b470e36 100644
--- a/AMDiS/src/time/RosenbrockStationary.cc
+++ b/AMDiS/src/time/RosenbrockStationary.cc
@@ -22,7 +22,7 @@ namespace AMDiS {
stageSolutions.resize(rm->getStages());
for (int i = 0; i < rm->getStages(); i++)
- stageSolutions[i] = new SystemVector(*solution);
+ stageSolutions[i] = new SystemVector(*solution);
}
@@ -48,6 +48,16 @@ namespace AMDiS {
*tmp *= rm->getA(i, j);
*stageSolution += *tmp;
}
+
+ clock_t first = clock();
+
+ for (unsigned int j = 0; j < boundaries.size(); j++) {
+ boundaries[j].vec->interpol(boundaries[j].fct);
+ *(boundaries[j].vec) -= *(stageSolution->getDOFVector(boundaries[j].row));
+ }
+
+ MSG("Solution %.5f seconds\n", TIME_USED(first, clock()));
+
timeRhsVec->set(0.0);
for (int j = 0; j < i; j++) {
@@ -121,4 +131,17 @@ namespace AMDiS {
ProblemVec::addVectorOperator(opRhs, row, &minusOne, &minusOne);
}
+
+ void RosenbrockStationary::addDirichletBC(BoundaryType type, int row, int col,
+ AbstractFunction<double, WorldVector<double> > *fct)
+ {
+ FUNCNAME("RosenbrockStationary::addDirichletBC()");
+
+ DOFVector<double>* vec = new DOFVector<double>(componentSpaces[row], "vec");
+ RosenbrockBoundary bound = {fct, vec, row, col};
+ boundaries.push_back(bound);
+
+ ProblemVec::addDirichletBC(type, row, col, vec);
+ }
+
}
diff --git a/AMDiS/src/time/RosenbrockStationary.h b/AMDiS/src/time/RosenbrockStationary.h
index 02e4585816fade349d41538d8d620258ded23e21..d2016203455208b34d191b67753474f4b3694bc3 100644
--- a/AMDiS/src/time/RosenbrockStationary.h
+++ b/AMDiS/src/time/RosenbrockStationary.h
@@ -29,6 +29,18 @@
namespace AMDiS {
+ struct RosenbrockBoundary
+ {
+ AbstractFunction<double, WorldVector<double> > *fct;
+
+ DOFVector<double> *vec;
+
+ int row;
+
+ int col;
+ };
+
+
class RosenbrockStationary : public ProblemVec
{
public:
@@ -95,6 +107,59 @@ namespace AMDiS {
tauGamma = ptr;
}
+ /// Adds a Dirichlet boundary condition, where the rhs is given by an
+ /// abstract function.
+ void addDirichletBC(BoundaryType type, int row, int col,
+ AbstractFunction<double, WorldVector<double> > *b);
+
+ /// Adds a Dirichlet boundary condition, where the rhs is given by a DOF
+ /// vector.
+ void addDirichletBC(BoundaryType type, int row, int col,
+ DOFVector<double> *vec)
+ {
+ FUNCNAME("RosenbrockStationary::addDirichletBC()");
+
+ ERROR_EXIT("Not yet supported!\n");
+ }
+
+ /// Adds a Neumann boundary condition, where the rhs is given by an
+ /// abstract function.
+ void addNeumannBC(BoundaryType type, int row, int col,
+ AbstractFunction<double, WorldVector<double> > *n)
+ {
+ FUNCNAME("RosenbrockStationary::addNeumannBC()");
+
+ ERROR_EXIT("Not yet supported!\n");
+ }
+
+ /// Adds a Neumann boundary condition, where the rhs is given by an DOF
+ /// vector.
+ void addNeumannBC(BoundaryType type, int row, int col,
+ DOFVector<double> *n)
+ {
+ FUNCNAME("RosenbrockStationary::addNeumannBC()");
+
+ ERROR_EXIT("Not yet supported!\n");
+ }
+
+ /// Adds Robin boundary condition.
+ void addRobinBC(BoundaryType type, int row, int col,
+ AbstractFunction<double, WorldVector<double> > *n,
+ AbstractFunction<double, WorldVector<double> > *r)
+ {
+ FUNCNAME("RosenbrockStationary::addRobinBC()");
+
+ ERROR_EXIT("Not yet supported!\n");
+ }
+
+ /// Adds a periodic boundary condition.
+ void addPeriodicBC(BoundaryType type, int row, int col)
+ {
+ FUNCNAME("RosenbrockStationary::addPeriodicBC()");
+
+ ERROR_EXIT("Not yet supported!\n");
+ }
+
protected:
void init();
@@ -112,6 +177,8 @@ namespace AMDiS {
double *tauPtr;
double *tauGamma;
+
+ std::vector<RosenbrockBoundary> boundaries;
};
}