This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Modelling internal BC

0 votes

Hi, I am new to FE and trying to model the following problem.

I have an elastic square with the following topology:

enter image description here

where blue and orange domains have different 'pressures', such that the force on the boundary between two domains is

enter image description here

Currently I define pressure as a discontinuous Lagrange element of degree 0, as in The FEniCS book Chapter 1.5.

My uff code looks like this:

V = VectorElement("Lagrange", triangle, 2)
W = FiniteElement("DG", triangle, 0)

u = TrialFunction(V)
v = TestFunction(V)
p = Coefficient(W)
E = Coefficient(W)
nu = Coefficient(W)

n = FacetNormal(triangle)


mu = E / (2.0*(1.0 + nu))

lmbda =  E*nu / ((1.0+nu) * (1.0 - 2.0*nu))

def sigma(v):
    return 2.0 * mu * sym(grad(v)) + lmbda*tr(sym(grad(v))) * Identity(2)

a = inner( sigma(u), sym(grad(v)) ) * dx
L = dot(p*n,v)*ds(0)

Since I want boundary integral to be calculated on every element face I define the Neumann BC as

class NeumannBoundary: public SubDomain{
    bool inside(const Array<double>& x, bool on_boundary) const{
        return true;
    }
};

I also use a simple Dirichlet BC:

class DirichletBoundary: public SubDomain{
    bool inside(const Array<double>& x, bool on_boundary) const{
        return (abs(x[0]) <= 10 && on_boundary);
    }
};

Where I prescribe u = (0,0)

When I run the solver, with zero 'pressure' in the outside (blue) material and 100 on the inside (orange) material, I get u=(0,0) everywhere.

If I set pressure of both regions to 100 I get the following solution:

enter image description here

Which is the same if I set pressure of the outer region to 100 and inner to 0.

So, I guess my boundary condition definition is wrong. Hence is the question: How to define internal surface conditions? (Sorry for the long intro)

asked Apr 27, 2016 by mdelmans FEniCS Novice (140 points)
...