Spatially varying parameter in a simple Heat Equation

Question

If we want to solve the heat equation with a spatially varying thermal diffusion coefficient. How can it be done in C++ in FEniCS?

By spatially varying, I mean that the diffusion coefficient is not constant throughout the mesh. It is changing with location. It's ok if the mesh is Gmsh or FEniCS mesh.

If anyone has done similar thing in past with any problem, please consider sharing a code snippet. That will be very helpful.

Thanks,

Comments

Are you looking for functionality which is not offered by Expression?

---

Hi nate,

I understand how an expression can be used for a parameter that varies in space. But can it be used for a problem like this:

Consider a 2D mesh divided in 2 layers with different material properties. Initially, I want the material properties to be constant in a layer. Thereafter, I want them to vary spatially in each layer. Can you please tell me how can this be done?

Answer 1

Given that you don't want to use Expression and the explanation of your problem is rather cryptic, here's a simple implementation of formulating a diffusion coefficient by projecting a piecewise function onto a FE function space.

from dolfin import *

mesh = UnitSquareMesh(32, 32)

# Set up left and right side subdomains
left_side = AutoSubDomain(lambda x, on_exterior: x[0] <= 0.5)
right_side = AutoSubDomain(lambda x, on_exterior: x[0] >= 0.5)

# Assign a CellFunction we'll use to construct the volume integration element
# Label the left 1 and the right 2
cf = CellFunction('size_t', mesh, 0)
left_side.mark(cf, 1)
right_side.mark(cf, 2)
dx = Measure('dx', domain=mesh, subdomain_data=cf)

# Diffusion coefficient space
V_material = FunctionSpace(mesh, 'DG', 0)
D_trial, D_test = TrialFunction(V_material), TestFunction(V_material)
D = Function(V_material)  # The diffusion coefficient

# Initially we project constant values
D_left = Constant(1.0)
D_right = Constant(10.0)
a = D_trial*D_test*dx
L = D_left*D_test*dx(1) + D_right*D_test*dx(2)
solve(a == L, D)

plot(D, interactive=True)

# Now we can project expressions
D_left = Expression('x[0] + 1.0')
D_right = Expression('2.0 - x[0]')
a = D_trial*D_test*dx
L = D_left*D_test*dx(1) + D_right*D_test*dx(2)
solve(a == L, D)

plot(D, interactive=True)

# Solve Poisson equation with the diffusion coefficient, D
V = FunctionSpace(mesh, 'CG', 1)
u, v = TrialFunction(V), TestFunction(V)
a = dot(D*grad(u), grad(v))*dx
L = Constant(1.0)*v*dx

soln = Function(V)
solve(a == L, soln, [DirichletBC(V, Constant(0.0), 'on_boundary')])
plot(soln, interactive=True)

If you don't need something this complicated, you could simply use:

from dolfin import *

mesh = UnitSquareMesh(32, 32)

# Set up left and right side subdomains
left_side = AutoSubDomain(lambda x, on_exterior: x[0] <= 0.5)
right_side = AutoSubDomain(lambda x, on_exterior: x[0] >= 0.5)

# Assign a CellFunction we'll use to construct the volume integration element
# Label the left 1 and the right 2
cf = CellFunction('size_t', mesh, 0)
left_side.mark(cf, 1)
right_side.mark(cf, 2)
dx = Measure('dx', domain=mesh, subdomain_data=cf)

D_left = Expression('x[0] + 1.0')
D_right = Expression('2.0 - x[0]')

V = FunctionSpace(mesh, 'CG', 1)
u, v = TrialFunction(V), TestFunction(V)
a = dot(D_left*grad(u), grad(v))*dx(1) + dot(D_right*grad(u), grad(v))*dx(2) 
L = Constant(1.0)*v*dx

soln = Function(V)
solve(a == L, soln, [DirichletBC(V, Constant(0.0), 'on_boundary')])
plot(soln, interactive=True)

Comments

Dear Nate,

Thank you so much for taking the time and effort to answer my question. Your answer is very helpful. I' ll try to create a C++ version of the same.

The problem I am solving is that of wave propagation and I want to use a Gmsh mostly. But if things don't work out, I can use a FENiCS mesh. So, in my python code, I was successfully able to assign different material properties like E (Youngs modulus), nu (Poisson Ratio), rho (Density) to different subdomains. I could do this even for a Gmsh mesh. However, I was facing trouble doing it in C++. So I was trying the following code, but it does not work for parallel. However, I think if I try to translate your approach to C++, it might work in parallel.

  auto V0 = std::make_shared<DG0::FunctionSpace>(mesh);

  MeshFunction<std::size_t> subdomains(mesh, mesh.topology().dim(), true);

  Omega0 Omega_0;  // This subdomain is defined as (x[1] <= 0.5) for a Gmsh unitsquare
  Omega_0.mark(subdomains, 0);

  Omega1 Omega_1;  // This subdomain is defined as (x[1] >= 0.5) for a Gmsh unitsquare
  Omega_1.mark(subdomains, 1);

  // Material parameters -->  mass density
  Function rho(V0);
  Array<double> rho_values(2);
  rho_values[0] = 0.5;
  rho_values[1] = 2;

  std:size_t cell_no, subdomain_no;
  for (CellIterator c(mesh); !c.end(); ++c)
  {
   cell_no = c->index();                                 // Cell number
   subdomain_no = subdomains.values()[cell_no];         
   // Number subdomain identifier that owns the cell number " cell_no " 
   rho.vector()->setitem(cell_no, rho_values[subdomain_no]);  

   std::cout << "" << std::endl;
   std::cout << "- sub_domain_number: " << subdomains.values()[cell_no] << std::endl;
   std::cout << "- cell_number: " << cell_no << std::endl;
   std::cout << "- rho_in_cell: " << rho.vector()->getitem(cell_no) << std::endl;
  }

  plot(mesh);
  interactive();
  plot(subdomains);
  interactive();  
  plot(rho);
  interactive();  

---

Hi Nate,

Your dx() approach in python worked in parallel successfully. Is there a way to do this in c++? I can write the expression and the solve function in addition to what I have noted below, but do you know how to express the dx = Measure('dx', domain=mesh, subdomain_data=cf) part in c++ and UFL? Is that statement necessary in c++ or can I just define the linear and bilinear form with dx(1) and dx(2) in UFL using python? Thanks for your help.

class left_side : public SubDomain
  {
    bool inside(const Array<double>& x, bool on_boundary) const
    { return(x[1] <= 0.5); }
  };

  class right_side : public SubDomain
  {
    bool inside(const Array<double>& x, bool on_boundary) const
    { return(x[1] >= 0.5); }
  };

  auto mesh = std::make_shared<UnitSquareMesh>(20, 20);
 // Mesh mesh("Square.xml.gz");

  MeshFunction<std::size_t> subdomains(mesh, mesh->topology().dim());
  subdomains = 0;
 //MeshFunction<std::size_t> subdomains(mesh, mesh.topology().dim(), 0);
 // CellFunction<std::size_t>subdomains(mesh, (dolfin::MPI::rank(mesh->mpi_comm()) == 0));

  left_side Left;
  Left.mark(subdomains, 1);
  right_side Right;
  Right.mark(subdomains, 2);

---

You should be able to use dx(idx) in ufl. Then for the C++ code look into Form::set_cell_domains() here and pass in your subdomains MeshFunction

---

Hi Nate,

I tried on of your approach in C++ and it worked fr a simple problem, but I got stuck where I actually need it - in my problem. I have posted the Question here.

https://groups.google.com/forum/#!topic/fenics-support/FNu7ZTXJwEY

I can post it in this forum if you want.

Please tell me if you have a solution to this.

Thanks a lot!!