Normal Dirichlet BC

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):
    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

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[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)

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])

with, "normal.bp", [solution]) as vtx:

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)
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.

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],0*x[1],0*x[2]+1])

works, while

f = dolfinx.fem.Function(ufl.FiniteElement("BDM", mesh.ufl_cell(), 1)) 
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.

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