Commit 0719b3b4 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Added boundary conditions to RosenbrockStationary, several small bugfixes and other small changes.

parent 3277e26c
...@@ -887,7 +887,7 @@ namespace AMDiS { ...@@ -887,7 +887,7 @@ namespace AMDiS {
TEST_EXIT(fct)("No function defined!\n"); TEST_EXIT(fct)("No function defined!\n");
int deg = 2 * vec1.getFeSpace()->getBasisFcts()->getDegree(); int deg = fct->getDegree();
Quadrature* quad = Quadrature* quad =
Quadrature::provideQuadrature(vec1.getFeSpace()->getMesh()->getDim(), deg); Quadrature::provideQuadrature(vec1.getFeSpace()->getMesh()->getDim(), deg);
FastQuadrature *fastQuad = FastQuadrature *fastQuad =
...@@ -925,7 +925,7 @@ namespace AMDiS { ...@@ -925,7 +925,7 @@ namespace AMDiS {
TEST_EXIT(fct)("No function defined!\n"); TEST_EXIT(fct)("No function defined!\n");
int deg = 2 * vec1.getFeSpace()->getBasisFcts()->getDegree(); int deg = fct->getDegree();
Quadrature* quad = Quadrature* quad =
Quadrature::provideQuadrature(vec1.getFeSpace()->getMesh()->getDim(), deg); Quadrature::provideQuadrature(vec1.getFeSpace()->getMesh()->getDim(), deg);
FastQuadrature *fastQuad = FastQuadrature *fastQuad =
......
...@@ -555,7 +555,7 @@ namespace AMDiS { ...@@ -555,7 +555,7 @@ namespace AMDiS {
*/ */
void interpol(AbstractFunction<T, WorldVector<double> > *fct); 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. /// Writes the data of the DOFVector to an output stream.
void serialize(std::ostream &out) void serialize(std::ostream &out)
......
...@@ -1232,6 +1232,21 @@ namespace AMDiS { ...@@ -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, void ProblemVec::addRobinBC(BoundaryType type, int row, int col,
AbstractFunction<double, WorldVector<double> > *n, AbstractFunction<double, WorldVector<double> > *n,
AbstractFunction<double, WorldVector<double> > *r) AbstractFunction<double, WorldVector<double> > *r)
......
...@@ -218,39 +218,48 @@ namespace AMDiS { ...@@ -218,39 +218,48 @@ namespace AMDiS {
/// Interpolates fct to \ref solution. /// Interpolates fct to \ref solution.
void interpolInitialSolution(std::vector<AbstractFunction<double, WorldVector<double> >*> *fct); 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, void addMatrixOperator(Operator *op, int i, int j,
double *factor = NULL, double *estFactor = NULL); 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, void addMatrixOperator(Operator &op, int i, int j,
double *factor = NULL, double *estFactor = NULL); double *factor = NULL, double *estFactor = NULL);
/// Adds an Operator to \ref rhs. /// Adds an operator to \ref rhs.
void addVectorOperator(Operator *op, int i, void addVectorOperator(Operator *op, int i,
double *factor = NULL, double *estFactor = NULL); double *factor = NULL, double *estFactor = NULL);
/// Adds an Operator to \ref rhs. /// Adds an operator to \ref rhs.
void addVectorOperator(Operator &op, int i, void addVectorOperator(Operator &op, int i,
double *factor = NULL, double *estFactor = NULL); 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, virtual void addDirichletBC(BoundaryType type, int row, int col,
AbstractFunction<double, WorldVector<double> > *b); 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, virtual void addDirichletBC(BoundaryType type, int row, int col,
DOFVector<double> *vec); 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, virtual void addNeumannBC(BoundaryType type, int row, int col,
AbstractFunction<double, WorldVector<double> > *n); 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, virtual void addRobinBC(BoundaryType type, int row, int col,
AbstractFunction<double, WorldVector<double> > *n, AbstractFunction<double, WorldVector<double> > *n,
AbstractFunction<double, WorldVector<double> > *r); AbstractFunction<double, WorldVector<double> > *r);
/// Adds periodic boundary conditions. /// Adds a periodic boundary condition.
virtual void addPeriodicBC(BoundaryType type, int row, int col); virtual void addPeriodicBC(BoundaryType type, int row, int col);
/** \brief /** \brief
......
...@@ -1412,6 +1412,7 @@ namespace AMDiS { ...@@ -1412,6 +1412,7 @@ namespace AMDiS {
Quadrature::quad_nd[3] = Quadrature::quad_3d; Quadrature::quad_nd[3] = Quadrature::quad_3d;
}; };
Quadrature* Quadrature::provideQuadrature(int dim_, int degree_) Quadrature* Quadrature::provideQuadrature(int dim_, int degree_)
{ {
FUNCNAME("Quadrature::provideQuadrature()"); FUNCNAME("Quadrature::provideQuadrature()");
...@@ -1439,12 +1440,14 @@ namespace AMDiS { ...@@ -1439,12 +1440,14 @@ namespace AMDiS {
return (quad_nd[dim_][degree_]); return (quad_nd[dim_][degree_]);
} }
Quadrature::~Quadrature() Quadrature::~Quadrature()
{ {
delete lambda; delete lambda;
delete[] w; delete [] w;
} }
double Quadrature::integrateStdSimplex(AbstractFunction<double, DimVec<double> > *f) double Quadrature::integrateStdSimplex(AbstractFunction<double, DimVec<double> > *f)
{ {
FUNCNAME("Quadrature::integrateStdSimplex()"); FUNCNAME("Quadrature::integrateStdSimplex()");
...@@ -1462,15 +1465,15 @@ namespace AMDiS { ...@@ -1462,15 +1465,15 @@ namespace AMDiS {
return result; return result;
} }
FastQuadrature* FastQuadrature::provideFastQuadrature(const BasisFunction* bas_fcts, FastQuadrature* FastQuadrature::provideFastQuadrature(const BasisFunction* bas_fcts,
const Quadrature& quad, const Quadrature& quad,
Flag init_flag) Flag init_flag)
{ {
FUNCNAME("FastQuadrature::FastQuadrature()"); FUNCNAME("FastQuadrature::provideFastQuuadrature()");
FastQuadrature *quad_fast; FastQuadrature *quad_fast;
#ifdef _OPENMP #ifdef _OPENMP
#pragma omp critical (provideFastQuad) #pragma omp critical (provideFastQuad)
#endif #endif
...@@ -1482,11 +1485,13 @@ namespace AMDiS { ...@@ -1482,11 +1485,13 @@ namespace AMDiS {
if ((*fast)->basisFunctions == bas_fcts && (*fast)->quadrature == &quad) if ((*fast)->basisFunctions == bas_fcts && (*fast)->quadrature == &quad)
break; 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; quad_fast = *fast;
} else { } else {
if (fast == fastQuadList.end()) { 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); fastQuadList.push_front(quad_fast);
...@@ -1502,6 +1507,7 @@ namespace AMDiS { ...@@ -1502,6 +1507,7 @@ namespace AMDiS {
return quad_fast; return quad_fast;
} }
void FastQuadrature::init(Flag init_flag) void FastQuadrature::init(Flag init_flag)
{ {
int dim = quadrature->getDim(); int dim = quadrature->getDim();
...@@ -1567,6 +1573,7 @@ namespace AMDiS { ...@@ -1567,6 +1573,7 @@ namespace AMDiS {
} }
} }
FastQuadrature::FastQuadrature(const FastQuadrature& fastQuad) FastQuadrature::FastQuadrature(const FastQuadrature& fastQuad)
{ {
FUNCNAME("FastQuadrature::FastQuadrature()"); FUNCNAME("FastQuadrature::FastQuadrature()");
...@@ -1609,6 +1616,7 @@ namespace AMDiS { ...@@ -1609,6 +1616,7 @@ namespace AMDiS {
} }
} }
FastQuadrature::~FastQuadrature() FastQuadrature::~FastQuadrature()
{ {
int nPoints = quadrature->getNumPoints(); int nPoints = quadrature->getNumPoints();
...@@ -1621,11 +1629,13 @@ namespace AMDiS { ...@@ -1621,11 +1629,13 @@ namespace AMDiS {
delete D2Phi; delete D2Phi;
} }
const double FastQuadrature::getSecDer(int q,int i ,int j, int m) const const double FastQuadrature::getSecDer(int q,int i ,int j, int m) const
{ {
return (D2Phi) ? (*D2Phi)[q][i][j][m] : 0.0; return (D2Phi) ? (*D2Phi)[q][i][j][m] : 0.0;
} }
const VectorOfFixVecs<DimMat<double> > *FastQuadrature::getSecDer(int q) const const VectorOfFixVecs<DimMat<double> > *FastQuadrature::getSecDer(int q) const
{ {
return D2Phi ? (&((*D2Phi)[q])) : NULL; return D2Phi ? (&((*D2Phi)[q])) : NULL;
......
...@@ -324,9 +324,7 @@ namespace AMDiS { ...@@ -324,9 +324,7 @@ namespace AMDiS {
* Constructs a FastQuadrature for the given Quadrature, BasisFunction, and * Constructs a FastQuadrature for the given Quadrature, BasisFunction, and
* flag. * flag.
*/ */
FastQuadrature(BasisFunction* basFcts, FastQuadrature(BasisFunction* basFcts, Quadrature* quad, Flag flag)
Quadrature* quad,
Flag flag)
: init_flag(flag), : init_flag(flag),
phi(NULL), phi(NULL),
grdPhi(NULL), grdPhi(NULL),
......
...@@ -22,14 +22,13 @@ namespace AMDiS { ...@@ -22,14 +22,13 @@ namespace AMDiS {
// create barycentric coords for each vertex of each side // create barycentric coords for each vertex of each side
const Element *refElement = Global::getReferenceElement(dim); 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 all element sides
for (int i = 0; i < dim + 1; i++) { for (int i = 0; i < dim + 1; i++) {
coords[i] = new VectorOfFixVecs<DimVec<double> >(dim, dim, DEFAULT_VALUE, coords[i] =
DimVec<double>(dim, new VectorOfFixVecs<DimVec<double> >(dim, dim, DEFAULT_VALUE,
DEFAULT_VALUE, DimVec<double>(dim, DEFAULT_VALUE, 0.0));
0.0));
// for each vertex of the side // for each vertex of the side
for (int k = 0; k < dim; k++) { for (int k = 0; k < dim; k++) {
int index = refElement->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim), i, k); int index = refElement->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim), i, k);
...@@ -78,11 +77,9 @@ namespace AMDiS { ...@@ -78,11 +77,9 @@ namespace AMDiS {
// for all element sides // for all element sides
for (int i = 0; i < dim + 1; i++) { for (int i = 0; i < dim + 1; i++) {
coords[i] = new VectorOfFixVecs<DimVec<double> >(dim, dim, DEFAULT_VALUE, coords[i] =
DimVec<double>(dim, new VectorOfFixVecs<DimVec<double> >(dim, dim, DEFAULT_VALUE,
DEFAULT_VALUE, DimVec<double>(dim, DEFAULT_VALUE, 0.0));
0.0)
);
// for each vertex of the side // for each vertex of the side
for (int k = 0; k < dim; k++) { for (int k = 0; k < dim; k++) {
int index = refElement->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim), i, k); int index = refElement->getVertexOfPosition(INDEX_OF_DIM(dim - 1, dim), i, k);
...@@ -230,7 +227,6 @@ namespace AMDiS { ...@@ -230,7 +227,6 @@ namespace AMDiS {
(*op)->getAssembler(omp_get_thread_num())->initElement(elInfo); (*op)->getAssembler(omp_get_thread_num())->initElement(elInfo);
for (int face = 0; face < dim + 1; face++) { for (int face = 0; face < dim + 1; face++) {
elInfo->getNormal(face, normal); elInfo->getNormal(face, normal);
Quadrature *quadrature = Quadrature *quadrature =
...@@ -251,12 +247,8 @@ namespace AMDiS { ...@@ -251,12 +247,8 @@ namespace AMDiS {
std::vector<double> f(nPoints, 0.0); std::vector<double> f(nPoints, 0.0);
if (robinOperators) if (robinOperators)
(*robinOperators)[face]->evalZeroOrder(nPoints, (*robinOperators)[face]->evalZeroOrder(nPoints, uhAtQp,
uhAtQp, NULL, NULL, &f[0], 1.0);
NULL,
NULL,
&f[0],
1.0);
std::vector<WorldVector<double> > grdUh(nPoints); std::vector<WorldVector<double> > grdUh(nPoints);
std::vector<WorldVector<double> > A_grdUh(nPoints); std::vector<WorldVector<double> > A_grdUh(nPoints);
......
...@@ -90,6 +90,7 @@ namespace AMDiS { ...@@ -90,6 +90,7 @@ namespace AMDiS {
VectorOfFixVecs<DimVec<double> >**coords; VectorOfFixVecs<DimVec<double> >**coords;
}; };
class NeumannBC : public RobinBC class NeumannBC : public RobinBC
{ {
public: public:
...@@ -99,6 +100,13 @@ namespace AMDiS { ...@@ -99,6 +100,13 @@ namespace AMDiS {
FiniteElemSpace *colFeSpace = NULL) FiniteElemSpace *colFeSpace = NULL)
: RobinBC(type, j, NULL, rowFeSpace, colFeSpace) : RobinBC(type, j, NULL, rowFeSpace, colFeSpace)
{} {}
NeumannBC(BoundaryType type,
DOFVectorBase<double> *j,
FiniteElemSpace *rowFeSpace,
FiniteElemSpace *colFeSpace = NULL)
: RobinBC(type, j, NULL, rowFeSpace, colFeSpace)
{}
}; };
} }
......
...@@ -22,7 +22,7 @@ namespace AMDiS { ...@@ -22,7 +22,7 @@ namespace AMDiS {
stageSolutions.resize(rm->getStages()); stageSolutions.resize(rm->getStages());
for (int i = 0; i < rm->getStages(); i++) for (int i = 0; i < rm->getStages(); i++)
stageSolutions[i] = new SystemVector(*solution); stageSolutions[i] = new SystemVector(*solution);
} }
...@@ -48,6 +48,16 @@ namespace AMDiS { ...@@ -48,6 +48,16 @@ namespace AMDiS {
*tmp *= rm->getA(i, j); *tmp *= rm->getA(i, j);
*stageSolution += *tmp; *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); timeRhsVec->set(0.0);
for (int j = 0; j < i; j++) { for (int j = 0; j < i; j++) {
...@@ -121,4 +131,17 @@ namespace AMDiS { ...@@ -121,4 +131,17 @@ namespace AMDiS {
ProblemVec::addVectorOperator(opRhs, row, &minusOne, &minusOne); 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);
}
} }
...@@ -29,6 +29,18 @@ ...@@ -29,6 +29,18 @@
namespace AMDiS { namespace AMDiS {
struct RosenbrockBoundary
{
AbstractFunction<double, WorldVector<double> > *fct;
DOFVector<double> *vec;
int row;
int col;
};
class RosenbrockStationary : public ProblemVec class RosenbrockStationary : public ProblemVec
{ {
public: public:
...@@ -95,6 +107,59 @@ namespace AMDiS { ...@@ -95,6 +107,59 @@ namespace AMDiS {
tauGamma = ptr; 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: protected:
void init(); void init();
...@@ -112,6 +177,8 @@ namespace AMDiS { ...@@ -112,6 +177,8 @@ namespace AMDiS {
double *tauPtr; double *tauPtr;
double *tauGamma; double *tauGamma;
std::vector<RosenbrockBoundary> boundaries;
}; };
} }
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment