I would like to define some of the parts of my problem in a separated function. To demonstrate what I want, I modified the Poisson demo and took the boundary conditions out of main() into construct_bc(). This way I get Segmentation Fault. How can I do it right?
#include <dolfin.h>
#include "Poisson.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;
}
};
DirichletBC construct_bc(Poisson::FunctionSpace const& V) {
Constant u0(0.0);
DirichletBoundary boundary;
return DirichletBC(V, u0, boundary);
}
int main()
{
// Create mesh and function space
UnitSquareMesh mesh(32, 32);
Poisson::FunctionSpace V(mesh);
// Define boundary condition
DirichletBC bc = construct_bc(V);
// Define variational forms
Poisson::BilinearForm a(V, V);
Poisson::LinearForm L(V);
Source f;
dUdN g;
L.f = f;
L.g = g;
// Compute solution
Function u(V);
solve(a == L, u, *bc);
return 0;
}