Auto adaptive Poisson equation (C++)¶
Implementation¶
Running this demo requires the files: main.cpp
,
AdaptivePoisson.ufl
and CMakeLists.txt
.
Under construction.
#include <dolfin.h>
#include "AdaptivePoisson.h"
using namespace dolfin;
// Source term (right-hand side)
class Source : public Expression
{
void eval(Array<double>& values, const Array<double>& x) const
{
double dx = x[0] - 0.5;
double dy = x[1] - 0.5;
values[0] = 10*exp(-(dx*dx + dy*dy) / 0.02);
}
};
// Normal derivative (Neumann boundary condition)
class dUdN : public Expression
{
void eval(Array<double>& values, const Array<double>& x) const
{ values[0] = sin(5*x[0]); }
};
// Sub domain for Dirichlet boundary condition
class DirichletBoundary : public SubDomain
{
bool inside(const Array<double>& x, bool on_boundary) const
{ return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS; }
};
int main()
{
// Create mesh and define function space
auto mesh = std::make_shared<UnitSquareMesh>(8, 8);
auto V = std::make_shared<AdaptivePoisson::BilinearForm::TrialSpace>(mesh);
// Define boundary condition
auto u0 = std::make_shared<Constant>(0.0);
auto boundary = std::make_shared<DirichletBoundary>();
auto bc = std::make_shared<DirichletBC>(V, u0, boundary);
// Define variational forms
auto a = std::make_shared<AdaptivePoisson::BilinearForm>(V, V);
auto L = std::make_shared<AdaptivePoisson::LinearForm>(V);
auto f = std::make_shared<Source>();
auto g = std::make_shared<dUdN>();
L->f = f;
L->g = g;
// Define Function for solution
auto u = std::make_shared<Function>(V);
// Define goal functional (quantity of interest)
auto M = std::make_shared<AdaptivePoisson::GoalFunctional>(mesh);
// Define error tolerance
double tol = 1.e-5;
// Solve equation a = L with respect to u and the given boundary
// conditions, such that the estimated error (measured in M) is less
// than tol
std::vector<std::shared_ptr<const DirichletBC>> bcs({bc});
auto problem = std::make_shared<LinearVariationalProblem>(a, L, u, bcs);
AdaptiveLinearVariationalSolver solver(problem, M);
solver.parameters("error_control")("dual_variational_solver")["linear_solver"]
= "cg";
solver.parameters("error_control")("dual_variational_solver")["symmetric"]
= true;
solver.solve(tol);
solver.summary();
// Output solution to XDMF for Paraview
XDMFFile("initial_solution.xdmf").write(u->root_node());
XDMFFile("final_solution.xdmf").write(u->leaf_node());
return 0;
}