Multiple Subdomains in C++

Question

I'm trying to apply the Different material tutorial to C++.
As a simple test, I'm trying to solve the Poisson equation times a scalar parameter k on a square mesh where the value of k is different depending on whether we are on the left or right half of the mesh.

In particular, I don't see how I could translate in C++ the python code

k_values = [1.5, 50]  # values of k in the two subdomains
for cell_no in range(len(subdomains.array())):
subdomain_no = subdomains.array()[cell_no]
k.vector()[cell_no] = k_values[subdomain_no]

Also, I am wondering whether it would be possible to somehow define two dx(0) and dx(1) so that the definition of the bilinear form in the .ufl file would look like

a = k0*inner(grad(u), grad(v))*dx(0) + k1*inner(grad(u), grad(v))*dx(1)

---

My incomplete and not working code so far is

# The bilinear form a(u, v) and linear form L(v) for
# Poisson's equation.
#
# Compile this form with FFC: ffc -l dolfin Poisson.ufl.

V = FiniteElement("Lagrange", triangle, 1)

u = TrialFunction(V)
v = TestFunction(V)
f = Coefficient(V)

k = Coefficient(V) 

a = k*inner(grad(u), grad(v))*dx
L = f*v*dx

---

// This program solves Poisson's equation
//
//     - div grad u(x, y) = f(x, y)
//
// on the [-1,1]x[-1,1] square with source f given by
//
//     f(x, y) = - 4*exp(x^2 + y^2) * (1 + x^2 + y^2)
//
// and Dirichlet boundary conditions given by
//
//     u(x, y) = exp(1)

#include <dolfin.h>
#include "Poisson.h"
#include <mshr.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];
    double dy = x[1];
    values[0] = - 4*exp(dx*dx + dy*dy)*(1 + dx*dx + dy*dy);
  }
};

// Sub domain for Dirichlet boundary condition
class DirichletBoundary : public SubDomain
{
  bool inside(const Array<double>& x, bool on_boundary) const
  {
    return ((x[0] < -1 + DOLFIN_EPS or x[0] > 1 - DOLFIN_EPS or x[1] < -1 + DOLFIN_EPS or x[1] > 1 - DOLFIN_EPS) and on_boundary);
  }
};

// Left half of mesh
class Omega0 : public SubDomain
{
  bool inside(const Array<double>& x, bool on_boundary) const
  {
    if (x[0] <= 1.0)
      return true;
    else
      return false;
  }
};

// Right half of mesh
class Omega1 : public SubDomain
{
  bool inside(const Array<double>& x, bool on_boundary) const
  {
    if (x[0] >= 1.0)
      return true;
    else
      return false;
  }
};


int main()
{
  auto circle = mshr::Rectangle(dolfin::Point(-1.0, 1.0, 0.0),dolfin::Point(1.0, -1.0, 0.0));
  auto mesh = mshr::generate_mesh(circle, 32);

  auto V = std::make_shared<Poisson::FunctionSpace>(mesh);

    Omega0 omega0;
  auto omega0_function = std::make_shared<MeshFunction<std::size_t>>(mesh, 2);
    omega0.mark(*omega0_function, 0);

    Omega1 omega1;
  auto omega1_function = std::make_shared<MeshFunction<std::size_t>>(mesh, 2);
    omega1.mark(*omega1_function, 1);

  // Define boundary condition
  auto u0 = std::make_shared<Constant>(exp(1.0));
  auto boundary = std::make_shared<DirichletBoundary>();
  DirichletBC bc(V, u0, boundary);

  // Define variational forms
  Poisson::BilinearForm a(V, V);
  Poisson::LinearForm L(V);

  auto f = std::make_shared<Source>();
  L.f = f;

  // Compute solution
  Function u(V);
  solve(a == L, u, bc);

  // Save solution in VTK format
  File file("poisson.pvd");
  file << u;

  // Plot solution
  plot(u);
  interactive();

  return 0;
}

Thanks for your help!

Answer 1

It is possible to use

a = k0*inner(grad(u), grad(v))*dx(0) + k1*inner(grad(u), grad(v))*dx(1)

You need to set the cell domains using Form::set_cell_domains() see here.

On a separate note, your code

k.vector()[cell_no] = k_values[subdomain_no]

is invalid. The local and global cell numbers do not line up with the DoFs of the problem. Consider projecting or interpolating a piecewise function onto the finite element space to formulate k instead.

Comments

Hi nate,
Thank you for your answer.

I'm trying to apply you answer on my code, but so far I have been running into some problems.

Forgive my lack of understanding, just to be sure, what you're suggesting is something of the form

MeshFunction<std::size_t> subdomains(mesh, mesh->topology().dim()); 
subdomains = 0;

Omega0 omega0;
omega0.mark(subdomains, 1);

Omega1 omega1;
omega1.mark(subdomains, 2);

// Define variational forms
Poisson::BilinearForm a(V, V);
 Poisson::LinearForm L(V);

a.set_cell_domains(subdomains);

That is, I define my 2 subdomains, mark them into the object subdomains
and the use your suggestion to send the information about the subdomains of the bilinear form a in the .ufl file, where

a = k0*inner(grad(u), grad(v))*dx(1) + k1*inner(grad(u), grad(v))*dx(2)

is that right?

I seem to have some problems with pointers and such though, since the error reads

main.cpp:89:31: error: no matching function for call to ‘Poisson::Form_a::set_cell_domains(dolfin::MeshFunction<long unsigned int>&)’
  a.set_cell_domains(subdomains);
                               ^
In file included from /usr/include/dolfin/fem/dolfin_fem.h:20:0,
                 from /usr/include/dolfin.h:31,
                 from main.cpp:13:
/usr/include/dolfin/fem/Form.h:310:10: note: candidate: void dolfin::Form::set_cell_domains(std::shared_ptr<const dolfin::MeshFunction<long unsigned int> >)
     void set_cell_domains(std::shared_ptr<const MeshFunction<std::size_t>>

---

looks like you need to dereference subdomains into a shared_ptr. Or just instantiate subdomains as a shared_ptr from the start via something like

auto subdomains = std::make_shared<MeshFunction<std::size_t>>(mesh, mesh->topology().dim());

---

Yes, something like that.
With a couple more modifications the code compiles, but I now have a segmentation fault when launching the application ^^

I'll look into it again tomorrow.

Thanks for your help!

---

Of course I forgot a reference in a pointer... or something like that.

The whole thing works now. For anyone interested the code is

# Compile this form with FFC: ffc -l dolfin Poisson.ufl.

V = FiniteElement("Lagrange", triangle, 1)

u = TrialFunction(V)
v = TestFunction(V)
f = Coefficient(V)

k1 = Coefficient(V)
k2 = Coefficient(V)

a = k1*inner(grad(u), grad(v))*dx(1) + k2*inner(grad(u), grad(v))*dx(2)
L = f*v*dx

---

// This program solves Poisson's equation
//
//     - div grad u(x, y) = f(x, y)
//
// on the [-1,1]x[-1,1] square with source f given by
//
//     f(x, y) = - 4*exp(x^2 + y^2) * (1 + x^2 + y^2)
//
// and Dirichlet boundary conditions given by
//
//     u(x, y) = exp(1)

#include <dolfin.h>
#include "Poisson.h"
#include <mshr.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];
    double dy = x[1];
    values[0] = - 4*exp(dx*dx + dy*dy)*(1 + dx*dx + dy*dy);
  }
};

// Sub domain for Dirichlet boundary condition
class DirichletBoundary : public SubDomain
{
  bool inside(const Array<double>& x, bool on_boundary) const
  {
    return ((x[0] < -1 + DOLFIN_EPS or x[0] > 1 - DOLFIN_EPS or x[1] < -1 + DOLFIN_EPS or x[1] > 1 - DOLFIN_EPS) and on_boundary);
  }
};

// Left half of mesh
class Omega0 : public SubDomain
{
  bool inside(const Array<double>& x, bool on_boundary) const
  {
    if (x[0] <= 0.0)
      return true;
    else
      return false;
  }
};

// Right half of mesh
class Omega1 : public SubDomain
{
  bool inside(const Array<double>& x, bool on_boundary) const
  {
    if (x[0] >= 0.0)
      return true;
    else
      return false;
  }
};

class Test : public Expression
{
  void eval(Array<double>& values, const Array<double>& x) const
  {
    values[0] = 5; //cos(x[0]);
  }
};

class Test2 : public Expression
{
  void eval(Array<double>& values, const Array<double>& x) const
  {
    values[0] = 10; //sin(x[0]);
  }
};

int main()
{
  auto rect = std::make_shared<mshr::Rectangle>(dolfin::Point(-1.0, 1.0, 0.0),dolfin::Point(1.0, -1.0, 0.0));
  auto mesh = mshr::generate_mesh(*rect, 100);

  auto V = std::make_shared<Poisson::FunctionSpace>(mesh);

  auto subdomains = std::make_shared<MeshFunction<std::size_t>>(mesh, mesh->topology().dim());
  *subdomains = 0;

  Omega0 omega0;
  omega0.mark(*subdomains, 1);

  Omega1 omega1;
  omega1.mark(*subdomains, 2);

  // Define boundary condition
  auto u0 = std::make_shared<Constant>(exp(1.0));
  auto boundary = std::make_shared<DirichletBoundary>();
  DirichletBC bc(V, u0, boundary);

  // Define variational forms
  Poisson::BilinearForm a(V, V);
  Poisson::LinearForm L(V);

  a.set_cell_domains(subdomains);

  auto f = std::make_shared<Source>();
  auto k11 = std::make_shared<Test>();
  auto k22 = std::make_shared<Test2>();
  a.k1 = k11;
  a.k2 = k22;
  L.f = f;

  // Compute solution
  Function u(V);
  solve(a == L, u, bc);

  // Save solution in VTK format
  File file("poisson.pvd");
  file << u;

  // Plot solution
  plot(u);
  interactive();

  return 0;
}