# Normal Dirichlet BC

**URL:** <https://fenicsproject.discourse.group/t/normal-dirichlet-bc/10582>\
**Category:** dolfinx\
**Created:** [March 3, 2023, 2:27pm UTC](https://fenicsproject.discourse.group/t/normal-dirichlet-bc/10582 "2023-03-03T14:27:14Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![kara](https://avatars.discourse-cdn.com/v4/letter/k/e480ec/32.png) [@kara](https://fenicsproject.discourse.group/u/kara)\
**Post date:** [March 3, 2023, 2:27pm UTC](https://fenicsproject.discourse.group/t/normal-dirichlet-bc/10582/1 "2023-03-03T14:27:14Z")

</div>

Hello,  
I have created a disk mesh using gmsh and now I want to analyze its displacement under certain constraints. Specifically, I want to impose a displacement on the outer boundary of the disk using `fem.dirichletbc` in FEniCSx. However, I am unable to find any relevant information in the FEniCSx documentation. Can you guide me on how to implement this?  
Here is my code:

```auto
# gmsh model...
domain, cell_markers, facet_markers = gmshio.model_to_mesh(model, MPI.COMM_SELF, 0)
element = ("CG", 2)
V = fem.VectorFunctionSpace(domain, element)
# Params
young = 1000
poisson = 0.5
mu = 0.5*young/(1 + poisson)
rho = 1
lambda_ = young/(1 - poisson**2)
g = 9.81

def epsilon(u):
    return sym(grad(u)) 
def sigma(u):
    return lambda_ * nabla_div(u) * Identity(u.geometric_dimension()) + 2*mu*epsilon(u)

u = TrialFunction(V)
v = TestFunction(V)
source = fem.Constant(domain, ScalarType((0, 0, 0)))
lhs = inner(sigma(u), epsilon(v)) * dx
rhs = dot(source, v) * dx

def outer_frontier(x):
    return np.logical_and(np.isclose(x[0]**2 + x[1]**2, 4), x[1] < 2)

fdim = domain.topology.dim - 1
facets = mesh.locate_entities_boundary(domain, fdim, outer_frontier)
marked_values = np.full_like(facets, 1)
facet_tag = mesh.meshtags(domain, fdim, facets, marked_values)

b_dofs = fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.find(1))

u_boundary = fem.Function(V)
# Find normal to the facet of each cell and apply displacement
# Is there any solution like:
displacement = 0.01
n = ufl.FacetNormal(domain)
norm = ufl.as_vector((n[0], n[1], 0))
# impose displacement according to the normal
u_boundary.vector.set(norm*displacement) 

bc_normal = fem.dirichletbc(u_boundary, b_dofs, V)

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [March 3, 2023, 3:17pm UTC](https://fenicsproject.discourse.group/t/normal-dirichlet-bc/10582/2 "2023-03-03T15:17:30Z")

</div>

You are here mixing `ufl` objects and `petsc`/`numpy` data.

Note that the `FacetNormal` is not well defined at vertices (it can have multiple values).  
Thus, if you can compute the normal as an expression of the current coordinates (if you have a disk at `x_c, y_c`, the normal can be written explicitly.

Consider the following code:

```python
import numpy as np
import dolfinx
from mpi4py import MPI

x_c = 0.5
y_c = 0.5
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
V = dolfinx.fem.VectorFunctionSpace(mesh, ("Lagrange", 2))
u = dolfinx.fem.Function(V)

def n(x, displacement=0.2):
    normal = displacement*np.vstack([(x[0]-x_c, x[1]-y_c) /
                                     np.sqrt((x[0]-x_c)**2 + (x[1]-y_c)**2)])
    return normal

mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
cells = dolfinx.mesh.compute_incident_entities(
    mesh, boundary_facets, mesh.topology.dim-1, mesh.topology.dim)
u.interpolate(n, cells=cells)
u.x.scatter_forward()

boundary_dofs = dolfinx.fem.locate_dofs_topological(
    V, mesh.topology.dim-1, boundary_facets)
bc = dolfinx.fem.dirichletbc(u, boundary_dofs)

solution = dolfinx.fem.Function(V)
dolfinx.fem.petsc.set_bc(solution.vector, [bc])
solution.x.scatter_forward()

with dolfinx.io.VTXWriter(mesh.comm, "normal.bp", [solution]) as vtx:
    vtx.write(0)

```

 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/6/6e819e4c6b7c345bbee0a167b3f41d3dde08f00c.png)  
For more complicated shapes that disks, you need to approximate the normal in the appropriate space, see for instance:

> [@Normal vector of just one of the boundaries](https://fenicsproject.discourse.group/t/normal-vector-of-just-one-of-the-boundaries/7855/4):
>
> How are you going to transform this into a numpy array? I would suggest using a projection into a suitable function space, which you can do as: which shows you how to get a projection of the facet normal over a certain boundary marked with a mesh tag. You can of course adapt this to take in your stress and project in to a suitable space.

---

<div class="post-metadata">

**Author:** ![kara](https://avatars.discourse-cdn.com/v4/letter/k/e480ec/32.png) [@kara](https://fenicsproject.discourse.group/u/kara)\
**Post date:** [March 3, 2023, 3:31pm UTC](https://fenicsproject.discourse.group/t/normal-dirichlet-bc/10582/3 "2023-03-03T15:31:12Z")

</div>

Thank you too much, I’ll implement this and let you know at the begining of week

---

<div class="post-metadata">

**Author:** ![matteoicardi](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/matteoicardi/32/4768_2.png) [@matteoicardi](https://fenicsproject.discourse.group/u/matteoicardi)\
**Post date:** [March 30, 2023, 6:10pm UTC](https://fenicsproject.discourse.group/t/normal-dirichlet-bc/10582/4 "2023-03-30T18:10:17Z")

</div>

I am trying to do the same for BDM elements but I get an “Interpolation data has the wrong shape” error

I can’t understand how the interpolation work and how it depends on the function space

any help?

---

<div class="post-metadata">

**Author:** ![kara](https://avatars.discourse-cdn.com/v4/letter/k/e480ec/32.png) [@kara](https://fenicsproject.discourse.group/u/kara)\
**Post date:** [March 30, 2023, 7:14pm UTC](https://fenicsproject.discourse.group/t/normal-dirichlet-bc/10582/5 "2023-03-30T19:14:46Z")

</div>

This error happen generally if you don’t handle well interpolation data.  
Bad:  
`f = dolfinx.fem.Function(V) `  
`f.interpolate(lambda x: 2*x)`  
Good:  
`f = dolfinx.fem.Function(V) `  
`f.interpolate(lambda x: 2*x[0])`  
NB: V is a function space  
Please share a minimum code if this doesn’t help you.

---

<div class="post-metadata">

**Author:** ![matteoicardi](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/matteoicardi/32/4768_2.png) [@matteoicardi](https://fenicsproject.discourse.group/u/matteoicardi)\
**Post date:** [March 30, 2023, 8:05pm UTC](https://fenicsproject.discourse.group/t/normal-dirichlet-bc/10582/6 "2023-03-30T20:05:36Z")

</div>

Thanks kara,

my problem is extremely simple, I just want to set my BDM function to be a constant in the whole domain [0,0,1].  
If I use a generic vector space (Lagrangian) the interpolation is trivial, just passing a vector, for some reason for BDM that does not work. However your trick with a lambda function surprisingly works also for BDM. strange but true

minimal example

```auto
f = dolfinx.fem.Function(ufl.FiniteElement("BDM", mesh.ufl_cell(), 1)) 
f.interpolate(lambda x: [0*x[0],0*x[1],0*x[2]+1])

```

works, while

```auto
f = dolfinx.fem.Function(ufl.FiniteElement("BDM", mesh.ufl_cell(), 1)) 
f.interpolate( np.vstack([0,0,1]))

```

doesn’t

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [March 30, 2023, 8:13pm UTC](https://fenicsproject.discourse.group/t/normal-dirichlet-bc/10582/7 "2023-03-30T20:13:01Z")

</div>

> [@matteoicardi](#):
>
> `f.interpolate( np.vstack([0,0,1]))`

The input to interpolate should be a function that takes in a (3,num\_points) array of coordinates, and returns a (mesh.geometry.dim, num\_points) array, i.e.

```python
def cond(x):
    values=numpy.zeros((3,x.shape[1)), dtype=np.float64)
    values[2] = 1
    return values

f.interpolate(cond)

```