# Normal Dirichlet BC

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:

``````# 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):
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**2 + x**2, 4), x < 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, n, 0))
# impose displacement according to the normal
u_boundary.vector.set(norm*displacement)

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

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:

``````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-x_c, x-y_c) /
np.sqrt((x-x_c)**2 + (x-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)
``````

For more complicated shapes that disks, you need to approximate the normal in the appropriate space, see for instance:

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

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?

This error happen generally if you don’t handle well interpolation data.
`f = dolfinx.fem.Function(V) `
`f.interpolate(lambda x: 2*x)`
Good:
`f = dolfinx.fem.Function(V) `
`f.interpolate(lambda x: 2*x)`
NB: V is a function space

1 Like

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

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

works, while

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

doesn’t

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.

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

f.interpolate(cond)
``````
2 Likes