Commit c8cd584a authored by Naumann, Andreas's avatar Naumann, Andreas
Browse files

speicherfehler in robinbc, ausgabe bei zu kleinem element in 1d, compiler...

speicherfehler in robinbc, ausgabe bei zu kleinem element in 1d, compiler fehler in rosenbrock und itl_oemsolver.hh
parent 9a7aa655
#include "BunnyProject.h"
#include "AdaptStationary.h"
/// RHS function
class F : public AbstractFunction<double, WorldVector<double> >
{
public:
F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {}
/// Implementation of AbstractFunction::operator().
double operator()(const WorldVector<double>& x) const
{
return -2 * x[0];
}
};
/// boundary
class G : public AbstractFunction<double, WorldVector<double> >
{
public:
/// Implementation of AbstractFunction::operator().
double operator()(const WorldVector<double>& x) const
{
return 10000.0;
}
};
Bunnydemo::Bunnydemo():
bunny("bunny"),
ballCenter(NULL),
adaptInfo(NULL),
adapt(NULL),
matrixOperator(NULL),
rhsOperator(NULL)
{}
Bunnydemo::~Bunnydemo() {
if(ballCenter != NULL)
delete ballCenter;
if(adaptInfo != NULL)
delete adaptInfo;
if(adapt != NULL)
delete adapt;
if(matrixOperator != NULL)
delete matrixOperator;
if(rhsOperator != NULL)
delete rhsOperator;
}
void Bunnydemo::create(string& filename) {
// ===== init parameters =====
Parameters::init(false, filename);
// ===== create projection =====
ballCenter = new WorldVector< double > ();
ballCenter->set(0.0);
new BallProject(1, VOLUME_PROJECTION, *ballCenter, 1.0);
// ===== create and init the scalar problem =====
bunny.initialize(INIT_ALL);
// === create adapt info ===
adaptInfo = new AdaptInfo("bunny->adapt");
// === create adapt ===
adapt = new AdaptStationary("bunny->adapt", &bunny, adaptInfo);
// ===== create matrix operator =====
matrixOperator = new Operator(bunny.getFeSpace());
matrixOperator->addSecondOrderTerm(new Laplace_SOT);
bunny.addMatrixOperator(matrixOperator);
// ===== create rhs operator =====
rhsOperator = new Operator(bunny.getFeSpace());
int degree = bunny.getFeSpace()->getBasisFcts()->getDegree();
rhsOperator->addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree)));
bunny.addVectorOperator(rhsOperator);
}
int Bunnydemo::solve(SolutionInformation& info) {
int retCode = adapt->adapt();
info.dofVec = bunny.getSolution();
return retCode;
}
#ifndef BUNNYPROJECT_H
#define BUNNYPROJECT_H
#include "Project.h"
#include "FixVec.h"
#include "ProblemScal.h"
class Bunnydemo : public Project {
ProblemScal bunny;
WorldVector< double >* ballCenter;
AdaptInfo* adaptInfo;
AdaptStationary* adapt;
Operator* matrixOperator;
Operator* rhsOperator;
public:
Bunnydemo();
~Bunnydemo();
void create(string& filename);
int solve(SolutionInformation&);
};
#endif
#include "ElliptProject.h"
#include "AdaptStationary.h"
/// Dirichlet boundary function
class G : public AbstractFunction<double, WorldVector<double> >
{
public:
/// Implementation of AbstractFunction::operator().
double operator()(const WorldVector<double>& x) const
{
return exp(-10.0 * (x * x));
}
};
/// RHS function
class F : public AbstractFunction<double, WorldVector<double> >
{
public:
F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {}
/// Implementation of AbstractFunction::operator().
double operator()(const WorldVector<double>& x) const
{
int dow = x.getSize();
double r2 = (x * x);
double ux = exp(-10.0 * r2);
return -(400.0 * r2 - 20.0 * dow) * ux;
}
};
Elliptdemo::Elliptdemo():
ellipt("ellipt"),
adaptInfo(NULL),
adapt(NULL),
matrixOperator(NULL),
rhsOperator(NULL)
{}
Elliptdemo::~Elliptdemo() {
if(matrixOperator != NULL)
delete matrixOperator;
if(rhsOperator != NULL)
delete rhsOperator;
if(adapt != NULL)
delete adapt;
if(adaptInfo != NULL)
delete adaptInfo;
}
void Elliptdemo::create(std::string& filename) {
// ===== init parameters =====
Parameters::init(true, filename);
// ===== create and init the scalar problem =====
ellipt.initialize(INIT_ALL);
// === create adapt info ===
adaptInfo = new AdaptInfo("ellipt->adapt");
// === create adapt ===
adapt = new AdaptStationary("ellipt->adapt", &ellipt, adaptInfo);
// ===== create matrix operator =====
matrixOperator = new Operator(ellipt.getFeSpace());
matrixOperator->addSecondOrderTerm(new Laplace_SOT);
ellipt.addMatrixOperator(matrixOperator);
// ===== create rhs operator =====
int degree = ellipt.getFeSpace()->getBasisFcts()->getDegree();
rhsOperator = new Operator(ellipt.getFeSpace());
rhsOperator->addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree)));
ellipt.addVectorOperator(rhsOperator);
// ===== add boundary conditions =====
ellipt.addDirichletBC(1, new G);
}
int Elliptdemo::solve(SolutionInformation& info) {
int retCode = adapt->adapt();
info.dofVec = ellipt.getSolution();
return retCode;
}
#ifndef ELLIPTPROJECT_H
#define ELLIPTPROJECT_H
#include "Project.h"
#include "ProblemScal.h"
class Elliptdemo : public Project {
ProblemScal ellipt;
AdaptInfo* adaptInfo;
AdaptStationary* adapt;
Operator* matrixOperator;
Operator* rhsOperator;
public:
Elliptdemo();
~Elliptdemo();
void create(std::string& filename);
int solve(SolutionInformation& info);
};
#endif
#include "HeatProject.h"
#include "TimedObject.h"
/// Dirichlet boundary function
class G : public AbstractFunction<double, WorldVector<double> >,
public TimedObject
{
public:
/// Implementation of AbstractFunction::operator().
double operator()(const WorldVector<double>& x) const
{
return sin(M_PI * (*timePtr)) * exp(-10.0 * (x * x));
}
};
/// RHS function
class F : public AbstractFunction<double, WorldVector<double> >,
public TimedObject
{
public:
F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {}
/// Implementation of AbstractFunction::operator().
double operator()(const WorldVector<double>& x) const
{
int dim = x.getSize();
double r2 = x * x;
double ux = sin(M_PI * (*timePtr)) * exp(-10.0 * r2);
double ut = M_PI * cos(M_PI * (*timePtr)) * exp(-10.0 * r2);
return ut -(400.0 * r2 - 20.0 * dim) * ux;
}
};
Heatdemo::Heatdemo():
heatSpace("heat->space"),
heat(heatSpace),
adaptInfo(NULL),
adaptInfoInitial(NULL),
A(NULL),
C(NULL),
Fop(NULL)
{}
Heatdemo::~Heatdemo() {
if(adaptInfo != NULL)
delete adaptInfo;
if(adaptInfoInitial != NULL)
delete adaptInfoInitial;
if(A != NULL)
delete A;
if(C != NULL)
delete C;
if(Fop != NULL)
delete Fop;
}
void Heatdemo::create(string& filename) {
// ===== init parameters =====
Parameters::init(false, filename);
// Parameters::readArgv(argc, argv);
// ===== create and init stationary problem =====
heatSpace.initialize(INIT_ALL);
// ===== create instationary problem =====
heat.initialize(INIT_ALL);
// create adapt info
adaptInfo = new AdaptInfo("heat->adapt");
// create initial adapt info
adaptInfoInitial = new AdaptInfo("heat->initial->adapt");
// create instationary adapt
adaptInstat = new AdaptInstationary("heat->adapt",
&heatSpace,
adaptInfo,
&heat,
adaptInfoInitial);
// ===== create rhs functions =====
int degree = heatSpace.getFeSpace()->getBasisFcts()->getDegree();
F *rhsFct = new F(degree);
rhsFct->setTimePtr(heat.getRHSTimePtr());
// ===== create operators =====
double one = 1.0;
double zero = 0.0;
// create laplace
A = new Operator(heatSpace.getFeSpace());
A->addSecondOrderTerm(new Laplace_SOT);
A->setUhOld(heat.getOldSolution());
if (*(heat.getThetaPtr()) != 0.0)
heatSpace.addMatrixOperator(A, heat.getThetaPtr(), &one);
if (*(heat.getTheta1Ptr()) != 0.0)
heatSpace.addVectorOperator(A, heat.getTheta1Ptr(), &zero);
// create zero order operator
C = new Operator(heatSpace.getFeSpace());
C->addZeroOrderTerm(new Simple_ZOT);
C->setUhOld(heat.getOldSolution());
heatSpace.addMatrixOperator(C, heat.getTau1Ptr(), heat.getTau1Ptr());
heatSpace.addVectorOperator(C, heat.getTau1Ptr(), heat.getTau1Ptr());
// create RHS operator
Fop = new Operator(heatSpace.getFeSpace());
Fop->addZeroOrderTerm(new CoordsAtQP_ZOT(rhsFct));
heatSpace.addVectorOperator(Fop);
// ===== create boundary functions =====
G *boundaryFct = new G;
boundaryFct->setTimePtr(heat.getBoundaryTimePtr());
heat.setExactSolution(boundaryFct);
}
int Heatdemo::solve(SolutionInformation& info) {
int retCode = adaptInstat->adapt();
info.dofVec = heatSpace.getSolution();
return retCode;
}
#ifndef HEATPROJECT_H
#define HEATPROJECT_H
#include "Project.h"
#include "ProblemScal.h"
#include "ProblemInstat.h"
#include "HeatProject.hh"
class Heatdemo : public Project {
ProblemScal heatSpace;
Heat heat;
AdaptInfo* adaptInfo;
AdaptInfo* adaptInfoInitial;
AdaptInstationary* adaptInstat;
Operator* A;
Operator* C;
Operator* Fop;
public:
Heatdemo();
~Heatdemo();
void create(string& filename);
int solve(SolutionInformation& info);
};
#endif
// ===========================================================================
// ===== instationary problem ================================================
// ===========================================================================
/// Instationary problem
class Heat : public ProblemInstatScal
{
public:
Heat(ProblemScal &heatSpace)
: ProblemInstatScal("heat", heatSpace)
{
// init theta scheme
theta = -1.0;
GET_PARAMETER(0, name + "->theta", "%f", &theta);
TEST_EXIT(theta >= 0)("theta not set!\n");
if (theta == 0.0) {
WARNING("You are using the explicit Euler scheme\n");
WARNING("Use a sufficiently small time step size!!!\n");
}
MSG("theta = %f\n", theta);
theta1 = theta - 1;
}
// ===== ProblemInstatBase methods ===================================
/// set the time in all needed functions!
void setTime(AdaptInfo *adaptInfo)
{
rhsTime = adaptInfo->getTime() - (1 - theta) * adaptInfo->getTimestep();
boundaryTime = adaptInfo->getTime();
tau1 = 1.0 / adaptInfo->getTimestep();
}
void closeTimestep(AdaptInfo *adaptInfo)
{
ProblemInstatScal::closeTimestep(adaptInfo);
WAIT;
}
// ===== initial problem methods =====================================
/// Used by \ref problemInitial to solve the system of the initial problem
void solve(AdaptInfo *adaptInfo)
{
problemStat->getSolution()->interpol(exactSolution);
}
/// Used by \ref problemInitial to do error estimation for the initial problem.
void estimate(AdaptInfo *adaptInfo)
{
double errMax, errSum;
errSum = Error<double>::L2Err(*exactSolution,
*(problemStat->getSolution()), 0, &errMax, false);
adaptInfo->setEstSum(errSum, 0);
adaptInfo->setEstMax(errMax, 0);
}
// ===== setting methods ===============================================
/// Sets \ref exactSolution;
void setExactSolution(AbstractFunction<double, WorldVector<double> > *fct)
{
exactSolution = fct;
}
// ===== getting methods ===============================================
/// Returns pointer to \ref theta.
double *getThetaPtr()
{
return &theta;
}
/// Returns pointer to \ref theta1.
double *getTheta1Ptr()
{
return &theta1;
}
/// Returns pointer to \ref tau1
double *getTau1Ptr()
{
return &tau1;
}
/// Returns pointer to \ref rhsTime.
double *getRHSTimePtr()
{
return &rhsTime;
}
/// Returns pointer to \ref theta1.
double *getBoundaryTimePtr()
{
return &boundaryTime;
}
private:
/// Used for theta scheme.
double theta;
/// theta - 1
double theta1;
/// 1.0 / timestep
double tau1;
/// time for right hand side functions.
double rhsTime;
/// time for boundary functions.
double boundaryTime;
/// Pointer to boundary function. Needed for initial problem.
AbstractFunction<double, WorldVector<double> > *exactSolution;
};
#include "NeumannProject.h"
#include "AdaptInfo.h"
#include "AdaptStationary.h"
class N : public AbstractFunction<double, WorldVector<double> >
{
public:
double operator()(const WorldVector<double>& x) const
{
return 1.0;
}
};
/// Dirichlet boundary function
class G : public AbstractFunction<double, WorldVector<double> >
{
public:
/// Implementation of AbstractFunction::operator().
double operator()(const WorldVector<double>& x) const
{
return exp(-10.0 * (x * x));
}
};
/// RHS function
class F : public AbstractFunction<double, WorldVector<double> >
{
public:
F(int degree) : AbstractFunction<double, WorldVector<double> >(degree) {}
/// Implementation of AbstractFunction::operator().
double operator()(const WorldVector<double>& x) const
{
int dow = x.getSize();
double r2 = (x*x);
double ux = exp(-10.0*r2);
return -(400.0*r2 - 20.0*dow)*ux;
}
};
Neumanndemo::Neumanndemo():
neumann("neumann"),
adaptInfo(NULL),
adapt(NULL),
matrixOperator(NULL),
rhsOperator(NULL)
{}
Neumanndemo::~Neumanndemo() {
if(matrixOperator != NULL)
delete matrixOperator;
if(rhsOperator != NULL)
delete rhsOperator;
if(adapt != NULL)
delete adapt;
if(adaptInfo != NULL)
delete adaptInfo;
}
void Neumanndemo::create(string& filename) {
// ===== init parameters =====
Parameters::init(true, filename);
// ===== create and init the scalar problem =====
neumann.initialize(INIT_ALL);
// === create adapt info ===
adaptInfo = new AdaptInfo("neumann->adapt", 1);
// === create adapt ===
adapt = new AdaptStationary("neumann->adapt",
&neumann,
adaptInfo);
// ===== create matrix operator =====
matrixOperator = new Operator(neumann.getFeSpace());
matrixOperator->addSecondOrderTerm(new Laplace_SOT);
neumann.addMatrixOperator(matrixOperator);
// ===== create rhs operator =====
int degree = neumann.getFeSpace()->getBasisFcts()->getDegree();
rhsOperator = new Operator(neumann.getFeSpace());
rhsOperator->addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree)));
neumann.addVectorOperator(rhsOperator);
// ===== add boundary conditions =====
neumann.addNeumannBC(1, new N);
neumann.addDirichletBC(2, new G);
}
int Neumanndemo::solve(SolutionInformation& info) {
int retCode = adapt->adapt();
info.dofVec = neumann.getSolution();
return retCode;
}
#ifndef NEUMANNPROJECT_H
#define NEUMANNPROJECT_H