Skip to content
Snippets Groups Projects
Commit 8a79f686 authored by Naumann, Andreas's avatar Naumann, Andreas
Browse files

added test data creation programs

parent a56db0a5
No related branches found
No related tags found
No related merge requests found
Showing
with 1143 additions and 7 deletions
...@@ -28,8 +28,7 @@ set(ENABLE_PARMETIS off) ...@@ -28,8 +28,7 @@ set(ENABLE_PARMETIS off)
option(ENABLE_UMFPACK "use umfpack solver" false) option(ENABLE_UMFPACK "use umfpack solver" false)
option(ENABLE_MKL "use the mkl" false) option(ENABLE_MKL "use the mkl" false)
SET(MKL_DIR "" CACHE PATH "MKL directory") SET(MKL_DIR "" CACHE PATH "MKL directory")
#option(ENABLE_TESTS "compile the tests" false) option(ENABLE_TESTS "compile the tests" false)
SET(ENABLE_TESTS false)
if(ENABLE_INTEL) if(ENABLE_INTEL)
Message("please set the icc manually") Message("please set the icc manually")
...@@ -233,6 +232,7 @@ if(ENABLE_UMFPACK) ...@@ -233,6 +232,7 @@ if(ENABLE_UMFPACK)
INSTALL(FILES ${LIB_DIR}/AMD/Lib/libamd.a INSTALL(FILES ${LIB_DIR}/AMD/Lib/libamd.a
DESTINATION lib/amdis/amd/) DESTINATION lib/amdis/amd/)
SET(RPM_DEPEND_STR "blas") SET(RPM_DEPEND_STR "blas")
LIST(APPEND AMDiS_LIBS amdis blas amd umfpack)
endif(ENABLE_UMFPACK) endif(ENABLE_UMFPACK)
if(ENABLE_MKL) if(ENABLE_MKL)
...@@ -254,7 +254,7 @@ include_directories(${SOURCE_DIR}) ...@@ -254,7 +254,7 @@ include_directories(${SOURCE_DIR})
add_library(amdis SHARED ${AMDIS_SRC} ${PARALLEL_DOMAIN_AMDIS_SRC}) add_library(amdis SHARED ${AMDIS_SRC} ${PARALLEL_DOMAIN_AMDIS_SRC})
add_library(compositeFEM SHARED ${COMPOSITE_FEM_SRC}) add_library(compositeFEM SHARED ${COMPOSITE_FEM_SRC})
target_link_libraries(compositeFEM amdis) target_link_libraries(compositeFEM amdis)
LIST(APPEND AMDiS_LIBS amdis boost_system boost_iostreams)
if(WIN32) if(WIN32)
SET(COMPILEFLAGS "${COMPILEFLAGS} -D_SCL_SECURE_NO_WARNINGS -D_CRT_SECURE_NO_WARNINGS") SET(COMPILEFLAGS "${COMPILEFLAGS} -D_SCL_SECURE_NO_WARNINGS -D_CRT_SECURE_NO_WARNINGS")
endif(WIN32) endif(WIN32)
...@@ -329,8 +329,6 @@ set(CPACK_RPM_PACKAGE_REQUIRES "boost-devel >= 1.42, ${RPM_DEPEND_STR}") ...@@ -329,8 +329,6 @@ set(CPACK_RPM_PACKAGE_REQUIRES "boost-devel >= 1.42, ${RPM_DEPEND_STR}")
set(CPACK_DEBIAN_PACKAGE_DEPENDS "libboost-dev (>= 1.42)") set(CPACK_DEBIAN_PACKAGE_DEPENDS "libboost-dev (>= 1.42)")
include(CPack) include(CPack)
if(ENABLE_TESTS) if(ENABLE_TESTS)
#ENABLE_TESTING() ENABLE_TESTING()
#add_test(demo_test run_test.sh) add_subdirectory(test)
#add_test(demo_test ${CMAKE_CTEST_COMMAND} --build-and-test "${CMAKE_SOURCE_DIR}/test/demo_test" "${CMAKE_BINARY_DIR}/test/demo_test" --build-generator "${CMAKE_GENERATOR}" --build-makeprogram "${CMAKE_MAKE_PROGRAM}" "${CMAKE_BINARY_DIR}/test/demo_test")
endif(ENABLE_TESTS) endif(ENABLE_TESTS)
project(tests)
add_subdirectory(datacreation)
project(testdatacreation)
message("source-dir: ${AMDiS_SOURCE_DIR}")
include_directories(${AMDiS_SOURCE_DIR})
file(GLOB PROJECTFILES src/*Project.cpp)
foreach(projectfile ${PROJECTFILES})
#create creatorname
get_filename_component(CppName ${projectfile}, NAME_WE)
string(REPLACE "Project" "" creatorsuffix ${CppName})
set(PROJECTINCLUDE "${CppName}.h")
configure_file(src/creator_base.cpp src/creator${creatorsuffix}.cpp @ONLY)
add_executable(creator${creatorsuffix} src/creator${creatorsuffix}.cpp src/ProjectList.cpp ${projectfile})
target_link_libraries(creator${creatorsuffix} ${AMDiS_LIBS})
endforeach(projectfile)
#include "BallProject.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 dim = x.getSize();
double r2 = x * x;
double ux = exp(-10.0*r2);
return -(400.0 * r2 - 20.0 * dim) * ux;
}
};
BallDemo::BallDemo():
ball("ball"),
ballCenter(NULL),
adaptInfo(NULL),
adapt(NULL),
matrixOperator(NULL),
rhsOperator(NULL)
{
}
BallDemo::~BallDemo() {
if(matrixOperator != NULL)
delete matrixOperator;
if(rhsOperator != NULL)
delete rhsOperator;
if(ballCenter != NULL)
delete ballCenter;
if(adapt != NULL)
delete adapt;
if(adaptInfo != NULL)
delete adaptInfo;
}
void BallDemo::create(std::string& filename) {
// ===== init parameters =====
Parameters::init(false, filename);
ballCenter = new WorldVector< double >();
ballCenter->set(0.0);
// ===== create projection =====
new AMDiS::BallProject(1,
BOUNDARY_PROJECTION,
*ballCenter,
1.0);
// ===== create and init the scalar problem =====
ball.initialize(INIT_ALL);
// === create adapt info ===
adaptInfo = new AdaptInfo("ball->adapt");
// === create adapt ===
adapt = new AdaptStationary("ball->adapt", &ball, adaptInfo);
// ===== create matrix operator =====
matrixOperator=new Operator(ball.getFeSpace());
matrixOperator->addSecondOrderTerm(new Laplace_SOT);
ball.addMatrixOperator(matrixOperator);
// ===== create rhs operator =====
int degree = ball.getFeSpace()->getBasisFcts()->getDegree();
rhsOperator = new Operator(ball.getFeSpace());
rhsOperator->addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree)));
ball.addVectorOperator(rhsOperator);
// ===== add boundary conditions =====
ball.addDirichletBC(1, new G);
}
int BallDemo::solve(SolutionInformation& solinfo) {
assert(adaptInfo!=NULL);
assert(adapt!=NULL);
assert(matrixOperator!=NULL);
assert(rhsOperator!=NULL);
int retCode = adapt->adapt();
solinfo.dofVec = ball.getSolution();
return retCode;
}
void createProjectList(ProjectList& list) {
list.clear();
Project* demo = new BallDemo();
ProjectInfo ballInfo(demo, "init/ball.dat.2d", "../testdata/balldata_2d");
list.push_back(ballInfo);
ballInfo = ProjectInfo(demo, "init/ball.dat.3d", "../testdata/balldata_3d");
list.push_back(ballInfo);
}
#ifndef BALLPROJECT_H
#define BALLPROJECT_H
#include "Project.h"
#include "FixVec.h"
#include "ProblemScal.h"
using namespace AMDiS;
class BallDemo : public Project {
ProblemScal ball;
WorldVector< double >* ballCenter;
AdaptInfo* adaptInfo;
AdaptStationary* adapt;
Operator* matrixOperator;
Operator* rhsOperator;
public:
BallDemo() ;
~BallDemo() ;
void create(std::string& ) ;
int solve(SolutionInformation&);
};
#endif
#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;
}
void createProjectList(ProjectList& list) {
}
#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;
}
void createProjectList(ProjectList& list) {
}
#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;
}
void createProjectList(ProjectList& list) {
}
#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;
}
void createProjectList(ProjectList& list) {
}
#ifndef NEUMANNPROJECT_H
#define NEUMANNPROJECT_H
#include "Project.h"
#include "ProblemScal.h"
class Neumanndemo : public Project {
ProblemScal neumann;
AdaptInfo* adaptInfo;
AdaptStationary* adapt;
Operator* matrixOperator;
Operator* rhsOperator;
public:
Neumanndemo();
~Neumanndemo();
void create(string filename);
int solve(SolutionInformation&);
};
#endif
#include "ParametricProject.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.0 * x[0];
}
};
Parametricdemo::Parametricdemo():
parametric("parametric"),
adaptInfo(NULL),
adapt(NULL),
matrixOperator(NULL),
rhsOperator(NULL)
{}
Parametricdemo::~Parametricdemo() {
if(matrixOperator != NULL)
delete matrixOperator;
if(rhsOperator != NULL)
delete rhsOperator;
if(adapt != NULL)
delete adapt;
if(adaptInfo != NULL)
delete adaptInfo;
}
void Parametricdemo::create(string filename) {
// ===== init parameters =====
Parameters::init(false, filename);
// ===== create and init the scalar problem =====
parametric.initialize(INIT_ALL);
// === create adapt info ===
adaptInfo = new AdaptInfo("parametric->adapt", 1);
// === create adapt ===
adapt = new AdaptStationary("parametric->adapt",
&parametric,
adaptInfo);
// ===== create matrix operator =====
matrixOperator = new Operator(parametric.getFeSpace());
matrixOperator->addSecondOrderTerm(new Laplace_SOT);
parametric.addMatrixOperator(matrixOperator);
// ===== create rhs operator =====
rhsOperator = new Operator( parametric.getFeSpace());
int degree = parametric.getFeSpace()->getBasisFcts()->getDegree();
rhsOperator->addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree)));
parametric.addVectorOperator(rhsOperator);
}
int Parametricdemo::solve(SolutionInformation& info) {
int retCode = adapt->adapt();
info.dofVec = parametric.getSolution();
return retCode;
}
void createProjectList(ProjectList& list) {
}
#ifndef PARAMETRICPROJECT_H
#define PARAMETRICPROJECT_H
#include "Project.h"
#include "ProblemScal.h"
#include "ParametricProject.hh"
class Parametricdemo : public Project {
ParametricSphere parametric;
AdaptInfo* adaptInfo;
AdaptStationary* adapt;
Operator* matrixOperator;
Operator* rhsOperator;
public:
Parametricdemo();
~Parametricdemo();
void create(string filename);
int solve(SolutionInformation&);
};
#endif
/** \brief
* Implementation of a scalar problem. In \ref buildBeforeCoarsen() parametric
* coordinates for the vertices are filled in a DOFVector. This DOFVector is
* used in \ref parametric to parametrize elements while mesh traversal.
*/
class ParametricSphere : public ProblemScal
{
public:
/// constructor
ParametricSphere(const char *name)
: ProblemScal(name),
parametric(NULL)
{}
/// destructor
~ParametricSphere()
{
for (int i = 0; i < Global::getGeo(WORLD); i++)
delete parametricCoords[i];
delete parametric;
}
/// initialization of the base class and creation of \ref parametric.
void initialize(Flag initFlag,
ProblemScal *adoptProblem = NULL,
Flag adoptFlag = INIT_NOTHING)
{
ProblemScal::initialize(initFlag, adoptProblem, adoptFlag);
// ===== create projection =====
WorldVector<double> ballCenter;
ballCenter.set(0.0);
new BallProject(1,
VOLUME_PROJECTION,
ballCenter,
1.0);
// ===== create coords vector =====
for (int i = 0; i < Global::getGeo(WORLD); i++)
parametricCoords[i] = new DOFVector<double>(this->getFeSpace(),
"parametric coords");
// ===== create parametric object =====
parametric = new ParametricFirstOrder(&parametricCoords);
// ===== enable parametric traverse =====
this->getMesh()->setParametric(parametric);
};
/** \brief
* Implementation of ProblemStatBase::buildBeforeCoarsen().
*/
void buildBeforeCoarsen(AdaptInfo *adaptInfo, Flag flag) {
FUNCNAME("ParametricSphere::buildAfterCoarsen()");
MSG("calculation of parametric coordinates\n");
int preDOFs = this->getFeSpace()->getAdmin()->getNumberOfPreDOFs(VERTEX);
int dim = this->getMesh()->getDim();
int dow = Global::getGeo(WORLD);
WorldVector<double> vertexCoords;
const DegreeOfFreedom **dof;
DegreeOfFreedom vertexIndex;
// ===== disable parametric traverse =====
this->getMesh()->setParametric(NULL);
TraverseStack stack;
ElInfo *elInfo = NULL;
elInfo = stack.traverseFirst(this->getMesh(), -1,
Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
while (elInfo) {
dof = elInfo->getElement()->getDof();
for (int i = 0; i < dim + 1; i++) {
vertexCoords = elInfo->getCoord(i);
Projection::getProjection(1)->project(vertexCoords);
for (int j = 0; j < dow; j++)
(*(parametricCoords[j]))[dof[i][preDOFs]] = vertexCoords[j];
}
elInfo = stack.traverseNext(elInfo);
}
// ===== enable parametric traverse =====
this->getMesh()->setParametric(parametric);
}
protected:
/// DOFVector storing parametric coordinates.
WorldVector<DOFVector<double>*> parametricCoords;
/// Parametrizes vertex coordinates while mesh traversal.
ParametricFirstOrder *parametric;
};
#include "PeriodicProject.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 0.0;
}
};
/// 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 dim = x.getSize();
double r2 = (x*x);
double ux = exp(-10.0*r2);
return -(400.0*r2 - 20.0*dim)*ux;
}
};
Periodicdemo::Periodicdemo():
periodic("periodic"),
adaptInfo(NULL),
adapt(NULL),
matrixOperator(NULL),
rhsOperator(NULL)
{}
Periodicdemo::~Periodicdemo() {
if(matrixOperator != NULL)
delete matrixOperator;
if(rhsOperator != NULL)
delete rhsOperator;
if(adapt != NULL)
delete adapt;
if(adaptInfo != NULL)
delete adaptInfo;
}
void Periodicdemo::create(string filename) {
// ===== init parameters =====
Parameters::init(false, filename);
// ===== create and init the scalar problem =====
periodic.initialize(INIT_ALL);
// === create adapt info ===
adaptInfo = new AdaptInfo("periodic->adapt", 1);
// === create adapt ===
adapt = new AdaptStationary("periodic->adapt",
&periodic,
adaptInfo);
// ===== create matrix operator =====
matrixOperator = new Operator(periodic.getFeSpace());
matrixOperator->addSecondOrderTerm(new Laplace_SOT);
periodic.addMatrixOperator(matrixOperator);
// ===== create rhs operator =====
int degree = periodic.getFeSpace()->getBasisFcts()->getDegree();
rhsOperator = new Operator(periodic.getFeSpace());
rhsOperator->addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree)));
periodic.addVectorOperator(rhsOperator);
// ===== add boundary conditions =====
periodic.addDirichletBC(1, new G);
// ===== add periodic conditions =====
periodic.addPeriodicBC(-1);
}
int Periodicdemo::solve(SolutionInformation& info) {
int retCode = adapt->adapt();
info.dofVec = periodic.getSolution();
return retCode;
}
void createProjectList(ProjectList& list) {
}
#ifndef PERIODICPROJECT_H
#define PERIODICPROJECT_H
#include "Project.h"
#include "ProblemScal.h"
class Periodicdemo : public Project {
ProblemScal periodic;
AdaptInfo* adaptInfo;
AdaptStationary* adapt;
Operator* matrixOperator;
Operator* rhsOperator;
public:
Periodicdemo();
~Periodicdemo();
void create(string filename);
int solve(SolutionInformation&);
};
#endif
#ifndef PROJECT_H
#define PROJECT_H
#include <string>
#include <list>
#include "DOFVector.h"
#include "SystemVector.h"
using namespace AMDiS;
using namespace std;
struct SolutionInformation {
SystemVector* sysVec;
DOFVector< double >* dofVec;
SolutionInformation():
sysVec(NULL),
dofVec(NULL)
{}
};
void write(SolutionInformation&, std::string filename);
bool compare(SolutionInformation&, std::string filename);
class Project {
public:
virtual void create(std::string& ) = 0;
virtual int solve(SolutionInformation&) = 0;
};
class ProjectInfo {
string initfilename;
Project* project;
string outputfilename;
public:
ProjectInfo(Project* project,
string initfilename,
string outputfilename):
initfilename(initfilename),
project(project),
outputfilename(outputfilename) {}
void create() {
assert(initfilename.length() > 0);
assert(project!=NULL);
project->create(initfilename);
}
int solve(SolutionInformation& info) {
assert(project!=NULL);
return project->solve(info);
}
inline string getFilename() const {
return outputfilename;
}
};
typedef list< ProjectInfo > ProjectList;
void createProjectList(ProjectList& );
#endif
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment