How to handle a subdomain inside an interval mesh?

Question

Hi!

I'm currently trying to establish a model to determine the behavior of a hanging cable tied to both ends (with application in structural engineering), and I'd like to modify locally the characteristics of the cable (corresponding to an instrumentation added on the cable).

I'm trying to do this thanks to the tutorial Handling Domains with Different Materials but I'm having a hard time understanding how subdomains have to be implemented in that case... In order to understand that, I'm trying to solve a (very) simplified problem where the solution I want to have in the end is a piecewise linear function equal to x on [0,0.3]U[0.5,1] and 3*x on [0.3,0.5] (for example).

Here is my code:

First, create mesh and subdomain:

from dolfin import *
import sys, math, numpy

mesh = UnitIntervalMesh(100)

subdomains = FacetFunction('size_t', mesh, 1) #MeshFunction should be used?

class domain_1(SubDomain):
      def inside(self, x, on_boundary):
      return True if 0.3 <= x[0] <= 0.5 else False

subdomains.set_all(0)
subdomain1 = domain_1()
subdomain1.mark(subdomains, 1)

At this time, plot(subdomains) returns the mesh with the two domains without problem, then:

dx=Measure('dx')[subdomains]
ds=Measure('ds')[subdomains]

V=FunctionSpace(mesh,"DG",1)

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

f1=Expression('x[0]',cell=interval)
f2=Expression('3*x[0]',cell=interval)

s=Function(V)

n=FacetNormal(mesh)

a=dot(grad(u),grad(v))*dx(0)+dot(grad(u),grad(v))*dx(1)
L=dot(inner(grad(f1),n),v)*ds(0)+dot(inner(grad(f2),n),v)*ds(1)

tol = 1E-16
def left_boundary(x, on_boundary):
    return on_boundary and abs(x[0]) < tol

def right_boundary(x, on_boundary):
    return on_boundary and abs(x[0]-1) < tol

Gamma_0_delr0 = DirichletBC(V, Constant(0.0), left_boundary)
Gamma_1_delr0 = DirichletBC(V, Constant(2.6), right_boundary)
bcs_delr0 = [Gamma_0_delr0, Gamma_1_delr0] #don't know how to impose the correct BCs...

solve(a==L,s,bcs_delr0)

plot(s)

This code returns an empty plot window ("nan" on y axis); which I believe is due to a bad implementation of boundary conditions... I tried to give it the solution at both ends of the mesh, but I believe that there is a lack of information on the boundaries of the subdomain 1, and maybe problems linked to the discontinuity of the function in 0.3 and 0.5...

I know this problem may be very basic, and I may lack of some mathematical background, so any help concerning subdomains and BCs would be much appreciated!

Regards.

Comments

Hi, do you have a specific reason for solving the problem with discontinuous elements?

---

Hi,

Indeed the DG here wasn't useful, I kept it from another code. When using classic CG the output is the function y=2.6*x that matches the BC I asked but not the slope break. I guess there's something wrong in my implementation of the Neumann conditions that should tell that f'(0.3)=1 and f'(0.5)=3 at the ends of the subdomain.

I guess the right question to ask is how can I write the linear form L so as to impose such BC on the subdomain?

Thanks for the answer.

---

Hi, what are the boundary conditions that you want? You have L as if you wanted to set the slope but then you use DirichletBC which fixes the value. You can't have both at the same time. Further, with only Neuman(slope) bcs the problem is not well posed.

---

This problem actually should allow me to understand how to impose Neumann BC on the boundaries of a subdomain that is included in the main domain, here is the description of the real problem I want to solve:

I have a free hanging elastic cable anchored to both ends at different levels, on which is installed an element that locally constrain the cable so that it has to be linear on this part. If we consider only the cable without the added element, the problem can be solved and its solution is called the catenary (See chap.1 of Max Irvine, Cable Structures).

Currently, I have a working code that returns the space coordinates of the points of the cable (represented as a line) but now I'm struggling to add the constraint. What I'd like to do is to impose the nullity of the curvature of the cable between two definite points. The problem is I have to impose 2 Dirichlet BCs on the anchors of the cable and another (Neumann?) condition that concerns the whole constrained subdomain of the cable.

I'm trying to implement a solution with Lagrange multiplier to constrain the PDE but in the examples I saw the multiplier is always defined over the whole domain and not only a part of it, and at this point I don't know if the problem I seek to resolve can be implemented that way. As a solution to this I was wondering if I could use Subdomains to locally define a Lagrange multiplier in order to constrain the equations only on the relevant part of the cable.

Hoping this explanation was clear enough given that I still don't understand very well FEM...

I can also provide pieces of my code or try to compute a minimal code showing more precisely what I want to do if needed, any hint or link to examples of similar problems is welcomed.

Thank you for your concern.

---

Okay. Please scratch what I said before. I forgot that you fix the values at the endpoints and the other bcs were for the middle of the domain. Also, suppose that you fix the height of the cable at two points in the middle. Don't you get a flat line in between them simply by setting the force in the middle to zero? Something like

from dolfin import *

A, B = 0.2, 0.8
middle = AutoSubDomain(lambda x: A - DOLFIN_EPS < x[0] < B + DOLFIN_EPS)

mesh = UnitIntervalMesh(100)
cell_f = CellFunction('size_t', mesh, 0)
middle.mark(cell_f, 1)

V = FunctionSpace(mesh, 'CG', 1)
u = TrialFunction(V)
v = TestFunction(V)

dx = Measure('dx')[cell_f]
a = inner(grad(u), grad(v))*dx
L = inner(Constant(-1), v)*dx(0)

bc_l = DirichletBC(V, Constant(1), 'x[0] < DOLFIN_EPS')
bc_r = DirichletBC(V, Constant(3), 'x[0] > 1- DOLFIN_EPS')
bc_c0 = DirichletBC(V, Constant(1.75), 'near(x[0], %g)' % A, 'pointwise')
bc_c1 = DirichletBC(V, Constant(1.75), 'near(x[0], %g)' % B, 'pointwise')
bcs = [bc_l, bc_r, bc_c0, bc_c1]

uh = Function(V)
solve(a == L, uh, bcs) 

plot(uh)
interactive()

---

Hi MiroK,

Sorry for the the delay, I've been working a bit more on the problem since I think I didn't ask the right question in the first place. Your answer definitely returns the output I wanted, and I came to a similar code while trying to solve the minimal problem I wrote above. Still, this solution works indeed for this minimal problem but I realized that it wouldn't work with my real problem. (I have to go with Lagrange multipliers)

In fact, what I would like to do is constrain locally (i.e a little part of my interval mesh in this case) a PDE thanks to Lagrange multipliers, so that the solution in this part of the mesh is linear (whereas it can be anything else that fits the PDE outside of the domain).

As a consequence, I have a PDE defined for the whole mesh (F(r), r=(x,y,z) the coordinates of the points of the mesh) and I want to constrain this PDE on the subdomain with LM so as to impose the nullity of the second derivative of r: (F(r)-Lambda*r'') should do this.

By the time I posted my question I didn't understand how exactly I could define the Lagrange multiplier lambda only on the subdomain I wanted, I've been looking for this for a while but it seems that I can't define a function such as lambda only on a part on the whole domain. I think there are 2 ways to answer the problem:

1. Try to define lambda on the subdomain using the constrained_domain argument of (Vector)FunctionSpace. In this case we define a whole function space on the subdomain and take lambda in this space. I didn't manage to make it work and I think that this argument is mostly dedicated to periodic boundary conditions so I may misunderstand its use.

2. Define lambda on the same function space as Trial and Test functions, so on the whole domain. In this case I have to cope with the fact that lambda has to have a definition inside and outside of the subdomain (I know the equation for lambda inside the subdomain but not outside). I'm currently trying to make this solution work and I'll post a resolution of the above minimal problem carried with LM as soon as possible.

Anyway thank you very much for your time,

Best regards.

---

Okay, but for the sake of completeness, please also include the math formulation of the problem that you are solving.