# problem with discontinuous flux at an internal interface

Asked by Nguyen Van Dang on 2013-03-04

Hello,
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]
-div(D(x,y)*u(x,y))=0
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,
Dang

My code is following:

ufl file
Poisson.ufl
# 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)
L = f*v*dx

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

main.cpp

#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();
}
exactsol.vector()->setitem(v->index(),val);
}
}

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);
plot(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
interface_fun.vector()->zero();
for (FacetIterator fc(mesh);!fc.end();++fc)
{
if (interface_marker.values()[fc->index()]==0)
{
for (VertexIterator v(*fc);!v.end();++v)
{
interface_fun.vector()->setitem(v->index(),1);
}
}
}
plot(interface_fun);

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

// set diffusion coefficients for two subdomains
for (CellIterator c(mesh);!c.end();++c)
{
kcoeff.vector()->setitem(c->index(),Darray[sub_domains.values()[c->index()]]);

}
plot(kcoeff);

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);
plot(u);

Function error(V);
for (VertexIterator v(mesh);!v.end();++v)
{
error.vector()->setitem(v->index(),u0.vector()->getitem(v->index())-u.vector()->getitem(v->index()));
}
plot(error);

return 0;
}

## Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Nguyen Van Dang
Solved:
2013-03-06
Last query:
2013-03-06