Constraint in RT VectorFunctionSpace

variational formulation
1

Dear FENICS community.

I am trying to enforce a Dirichlet boundary condition only on a component of a tensor element in a Raviart-Thomas tensor Space. Below is what I have tried:

from dolfin import *

L = 10
nx = 50
ny = 50

mesh = RectangleMesh(Point(0.0, 0.0), Point(L, L), nx, ny)
degree = 1

rt_tensor = VectorElement('RT', mesh.ufl_cell(), degree)
H = FunctionSpace(mesh, rt_tensor)

def boundary(x, on_boundary):
return on_boundary

zero_2d = Constant((0.0, 0.0))

bc_1 = DirichletBC(H.sub(0), zero_2d, boundary)
bc_2 = DirichletBC(H.sub(1), zero_2d, boundary)

What I have read:

https://fenicsproject.org/qa/5236/setting-boundary-condition-vector-component-dirichletbc/
https://fenicsproject.org/qa/12290/how-to-constraint-only-one-component-in-dg-method/

What I want:

For U(x,y) in H, with

U(x,y)=\begin{bmatrix} u_{11}(x,y) & u_{12}(x,y)\\ u_{21}(x,y) & u_{22}(x,y) \end{bmatrix},

I need to enforce only u_{12}(x,y)=u_{21}(x,y)=0 on some domain boundary \partial\Omega. But running H.sub(0).sub(i), i=0,... does not work. Moreover, running

bc_1 = DirichletBC(H.sub(0), zero_2d, boundary) # will enforce u_11=u_12=0

and

bc_2 = DirichletBC(H.sub(1), zero_2d, boundary) # will enforce u21=u_22=0

Any clue is appreciated.

Thanks!!

2

In the lowest-order Raviart-Thomas function space, the degrees of freedom set the normal components of the function on the facets of your mesh.
If you set the boundary conditions as you do above, the function will not vanish at the boundary, contrary to what I believe you mean with

it will only set the normal component of your function to zero. Consider e.g.:

from dolfin import *

mesh = UnitSquareMesh(5,5)
n = FacetNormal(mesh)

V = VectorFunctionSpace(mesh, 'RT', 1)

u = Expression((('x[0]+1','x[1]+1'), ('x[0]-1','x[1]-1')), degree=1)

g = interpolate(u, V)

bc = DirichletBC(V, Constant(((0,0), (0,0))), 'on_boundary')

print(assemble( (g.sub(0)**2 + g.sub(1)**2) * ds ), assemble( (inner(g.sub(0),n)**2 + inner(g.sub(1),n)**2) * ds ))

bc.apply(g.vector())

print(assemble( (g.sub(0)**2 + g.sub(1)**2) * ds ), assemble( (inner(g.sub(0),n)**2 + inner(g.sub(1),n)**2) * ds ))

I am not sure how to treat your original problem though. I came up with using a MixedElement with four components, but this may not suit your problem at all.

from dolfin import *

mesh = UnitSquareMesh(5,5)
x = SpatialCoordinate(mesh)
n = FacetNormal(mesh)

V_ele = FiniteElement('CG', mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, MixedElement([V_ele, V_ele, V_ele, V_ele]))

u = as_vector([ x[0]+1, x[1]+1, x[0]-1, x[1]-1 ])  
               
g = project(u, V)

bc1 = DirichletBC(V.sub(0), Constant(0.0), 'on_boundary')
bc2 = DirichletBC(V.sub(3), Constant(0.0), 'on_boundary')

print(assemble( g.sub(0)**2 * ds ), assemble( g.sub(1)**2 * ds ), assemble( g.sub(2)**2 * ds), assemble( g.sub(3)**2 * ds ))

bc1.apply(g.vector())
bc2.apply(g.vector())
      
print(assemble( g.sub(0)**2 * ds ), assemble( g.sub(1)**2 * ds ), assemble( g.sub(2)**2 * ds), assemble( g.sub(3)**2 * ds ))
3

Dear @volkerk, thank you for your response. Indeed I would like to have the normal component of the functions to be zero at the boundary. Unfortunately, the approach using four component for CG is no suitable for my case. I am trying to impose BC in the crossed components of the tensor space H = FunctionSpace(mesh, rt_tensor). This could be equivalent to impose BC in one of the component of a vector space RT, where

rt_ = FiniteElement('RT', mesh.ufl_cell(), degree)
RT=Functionspace(mesh, rt_)

But I am clueless there too.

4

Hi Jesus, if I understood correctly, then what you need to do is impossible with RT elements. This happens because they are not a component-wise element as doing a vector CG, but instead the dofs multiply vectors. This is what gives origin to “setting only the normal component”, so maybe you could try some Nitsche-like strategies for what you need. You cannot impose what you want even if you choose a boundary with the appropriate orientation, because u_12 and u_21 belong to different vectors.

What do you need exactly? In elasticity, the requirement of symmetry is imposed weakly, so that is another option you could consider.

Best!
NB

5

Dear Nicolas,

Thank you for taking the time to read this post. Indeed, I am trying to implement a Reissner-Mindlin plate following the approach from https://onlinelibrary.wiley.com/doi/epdf/10.1002/num.21698. The weak symmetry has been imposed successfully. But I’m having issues when imposing the boundary conditions different from the clamped case (hard-clamped). For instance, in the case of a simply supported plate we have that M_{11}=0 (or M_{22}=0) on \partial\Omega (see https://www.researchgate.net/publication/228574087_Finite_Elements_for_the_Reissner-Mindlin_Plate), but I don’t know how to implement this (or the case showed in the post). The closest I’ve come to imposing this boundary condition has been with bc_1 and bc_2, which is incorrect. Well, looking at my post I realize that what I really need is u_{11}(x,y)=u_{22}(x,y)=0.

6

Hi Jesus,

no worries, I’m glad I can be of assistance. I had a quick look to the second reference, and I believe some comments are in place:
(i) All of the hard-clamped type of BC are not attainable by RT elements. In such a context, take a look at Nitsche’s trick. I have always thought this answer is the best explanation available: https://scicomp.stackexchange.com/a/19913
(ii) u is the displacement, so the condition you mention (u_11=u_22=0) seems wrong, as u has only one index.
(iii) I believe your approach of setting things component-wise to be a bit limiting, as you would really need to study things too independently at each boundary, and of course doing anything but the square will be impossible :). Keeping that in mind, if what you really want is to impose M_11=M_22=0 at some portion of the domain, you could consider an RT tensor (given by two RT vectors) and then build it awkwardly as

M = as_tensor([[s1[0], s2[0]] , [s2[1], s1[1]]]

so that imposing a boundary condition on s2 could somehow help? I’m not sure, but it could be an explorable idea.

To conclude, if what you want is just an implementation of the problem, then go for Nitsche’s trick. If you are interested in convergence or other theoretical aspects, maybe go for something more oriented towards Lagrange multipliers. That would still be a difficult problem to analyze, and probably out of scope.

Best!
NB

7

Hi, Nicolas.

You’ve given me some very interesting ideas for the implementation. In response to your comments:

(i) I will definitely take a closer look at Nitsche’s trick for the implementation.

(ii) Indeed, the use of U(x,y) was only to exemplify in some way what I wanted to achieve.

(iii) That idea has been the one I have been trying to implement lately. The stress tensor is implemented creating an enriched element https://fenicsproject.org/qa/10359/enriched-spaces-for-peers-help-on-implementation/. This is achieved “akwardly” as
RT = VectorElement('RT', mesh.ufl_cell(), 1)
b_element = VectorElement('Bubble', mesh.ufl_cell(), 3)
H = MixedElement((RT, b_element))
rt_,b =rt, bubble= TrialFunctions(H)

and then

 `sigma =rt_+Curl(b)`

where Curl is written using as_tensor . With this approach, I am able to access to sigma_ij. However, as you said, this approach is very limited since, for example, I don’t know the exact solution at the boundaries of the simply supported case at each component of s in sigma= as_tensor([[s[0,0], s[0,1]] , [s[1,0], s1[1,1]]] for general R-M plates and the RT space doesn’t allow to enforce s[i,i]or s[j,j] trough DirichletBC.

8

Hahah its funny that you mention that post, as I posted it myself. I don’t think I can help you further, but I have found myself implementing mixed elasticity with the AFW elements (BDM(k+1) for stress, DG k the rest) which are more efficient (slightly) and more handeble. Anyway, hope you find a way to move on.

Best!

1 Like