problem with discontinuous flux at an internal interface

Asked by Nguyen Van Dang on 2013-03-04

I am working with a problem related to different materials and discontinuous flux at the interface.
Here, I am considering the Laplace equation over omega=[0,1]x[0,1]
Let u=u1 and D(x,y)=D1 on [0,1]x[0,1/2], u=u2 and D(x,y)=D2 on [0,1]x[1/2,1],
where u1 = x - 2*D2*y+1/2*D1+D2, u2 = x+D1*y
We can compute the flux at the interface y=1/2
Now I want to recover this solution by Fenics. To do this, I applied DirichletBC for extenal boundaries. Since the flux is discontinuous, I added kcoeff*dot(grad(u),n)*v*ds(0) to BilinearForm. However, this term did not effect on the results and the results were not correct. Can you please show me the way to fix it?
Thanks a lot,
Best regards,

My code is following:

ufl file
# Compile this form with FFC: ffc -l dolfin Poisson.ufl.
element = FiniteElement("Lagrange", triangle, 1)
u = TrialFunction(element)
v = TestFunction(element)
kcoeff = Coefficient(element);
n = FacetNormal("triangle");
f = Coefficient(element)
a = kcoeff*inner(grad(u), grad(v))*dx - kcoeff*dot(grad(u),n)*v*ds(0)
L = f*v*dx

# Compile this form with FFC: ffc -l dolfin FunctionSpace1Dtri.ufl
element = FiniteElement("DG", "triangle", 0)


#include <dolfin.h>
#include "Poisson.h"
#include "FunctionSpace1Dtri.h"

using namespace dolfin;

double Darray[2] = {3e-3,3e-3};

void ExactSol(Mesh mesh, FunctionSpace V, double*Darray,Function& exactsol)
  for (VertexIterator v(mesh);!v.end();++v)
    Point p = v->point();
    double val = 0.0;
    if (p.y()<=1.0/2.0) // u1
      val = p.x()-2.0*Darray[1]*p.y()+1.0/2.0*Darray[0]+Darray[1];
    else // u2
      val = p.x()+Darray[0]*p.y();

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

class Interface : public SubDomain
  bool inside(const Array<double>& x, bool on_boundary) const
    return fabs(x[1] - 0.5)<1e-10;

int main()
  // Create mesh
  UnitSquare mesh(8, 8);

  Poisson::FunctionSpace V(mesh);
  Function u0(V);
  ExactSol(mesh, V, Darray,u0);

  class DirichletBoundary : public SubDomain
    bool inside(const Array<double>& x, bool on_boundary) const
      return on_boundary;

  MeshFunction<unsigned int>sub_domains(mesh, mesh.topology().dim(),1);
  MeshFunction<unsigned int> interface_marker(mesh, mesh.topology().dim()-1,1);

  // Mark Omega0 as sub domain 0
  Omega0 omega0;
  omega0.mark(sub_domains, 0);

  Interface interface;
  interface.mark(interface_marker, 0);

  Function interface_fun(V); // for checking the interface
  for (FacetIterator fc(mesh);!fc.end();++fc)
     if (interface_marker.values()[fc->index()]==0)
       for (VertexIterator v(*fc);!v.end();++v)

  FunctionSpace1Dtri::FunctionSpace V1(mesh);
  Function kcoeff(V1);

  // set diffusion coefficients for two subdomains
  for (CellIterator c(mesh);!c.end();++c)


  DirichletBoundary dirichlet_boundary;
  DirichletBC bc0(V, u0, dirichlet_boundary);
  // Compute solution
  Function u(V);

  Poisson::BilinearForm a(V, V);
  a.kcoeff = kcoeff;
  //a.ds = interface_marker;
  Constant f(0.0);
  Poisson::LinearForm L(V);
  L.f = f;

  solve(a == L, u, bc0);

  Function error(V);
  for (VertexIterator v(mesh);!v.end();++v)

  return 0;

Question information

English Edit question
DOLFIN Edit question
No assignee Edit question
Solved by:
Nguyen Van Dang
Last query:
Last reply:
Nguyen Van Dang (dang-1032170) said : #1

I found the answer in UFL book.