3D Neumann bc from height dependent function

In a 2d elasticity model i use a BC neuman compress condition which depends from a function like:

k=1.5
yl_max=564.0
yr_max=675.5430
left_function = Expression(("0.027*k*(yl_max-x[1])", 0), degree=1, k=k, yl_max=yl_max) #Compression in left Boundary
rigth_function = Expression(("-0.027*k*(yr_max-x[1])", 0), degree=1, k=k, yr_max=yr_max) #Compression in rigth Boundary
left = AutoSubDomain(lambda x: near(x[0],0.0,tol))
rigth = AutoSubDomain(lambda x: near(x[0],x_max,tol))
# Definition of Neumann boundary condition domain
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
left.mark(boundaries, 1)
rigth.mark(boundaries, 2)
ds = ds(subdomain_data=boundaries)

In this case, i know the max value for “y” in rigth and left bonudaries.
Now im created a 3D model, like a box whit topography at the top side like:


but, i have diferents z max value along the faces to define the compression funcion and i need to create a funcion depend from the diferent values. I have not found a corret way to do this. I hope you can help me.

Here’s a hacky code that should do what you want. Basically, it:

  1. Extracts the top surface of the mesh into a submesh
  2. Flattens the submesh onto the xy-plane
  3. Creates a function in the xy-plane that maps the height of the top surface
  4. Interpolates the height map using a UserExpression to create the pressure boundary condition

Warning: I have only tested this in serial, not in parallel.

from fenics import *

# Create a mesh with varying height around the edges
msh = UnitCubeMesh(2, 2, 2)
msh.coordinates()[:, 2] += 0.25*msh.coordinates()[:, 2]*(msh.coordinates()[:, 0] + msh.coordinates()[:, 1])

# Extract a submesh consisting of only the top surface
bmsh = BoundaryMesh(msh, "exterior")
mff = MeshFunction("size_t", bmsh, bmsh.topology().dim(), 1)
AutoSubDomain(lambda x: near(x[0], 0.0)).mark(mff, 0)
AutoSubDomain(lambda x: near(x[0], 1.0)).mark(mff, 0)
AutoSubDomain(lambda x: near(x[1], 0.0)).mark(mff, 0)
AutoSubDomain(lambda x: near(x[1], 1.0)).mark(mff, 0)
AutoSubDomain(lambda x: near(x[2], 0.0)).mark(mff, 0)
subbmsh = SubMesh(bmsh, mff, 1)

# Flatten the submesh of the top surface into the xy-plane
CG_top = FunctionSpace(subbmsh, "CG", 1)
VCG_top = VectorFunctionSpace(subbmsh, "CG", 1)
xy = interpolate(Expression(("x[0]", "x[1]", "0"), degree=1), VCG_top)
floormsh = create_mesh(xy)

# Create a function that maps the height of the top surface
CG_floor = FunctionSpace(floormsh, "CG", 1)
z = interpolate(Expression("x[2]", degree=1), CG_top)
z_floor = Function(CG_floor)
z_floor.vector()[:] = z.vector()[:]

# Create a function representing the desired boundary condition
CG_3D = FunctionSpace(msh, "CG", 1)
class Pressure(UserExpression):
    def __init__(self, height_map, **kwargs):
        super().__init__(**kwargs)
        self._height_map = height_map
    def eval(self, values, x):
        values[0] = self._height_map(x[0], x[1], 0.0) - x[2]
    def value_shape(self):
        return tuple()
pressure = interpolate(Pressure(z_floor), CG_3D)
with XDMFFile("pressure.xdmf") as xdmf:
    xdmf.write(pressure)