Moving neumann BC

dolfinx
1

Hi,
I am trying to implement a moving Neumann boundary condition. I created a cantilever beam MWE example here, where I manually update the node of the load through a custom function. The code runs, but in the result nothing happens. The output of my expression out contains non-zero values, but when I interpolate this to my BC (load_func) it is all zeros.

1: Is manually finding the nodes to apply a load to a good way to achieve a moving load, or is there a simpler method?
2: Why is my interpolated load empty?
Thank you!

import numpy as np
import ufl
import dolfinx
from dolfinx import mesh, fem
from dolfinx.io import XDMFFile
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
from petsc4py import PETSc

print(f"dolfinx version: {dolfinx.__version__}")

# Create a uniform 2D rectangular mesh
mesh_domain = mesh.create_rectangle(MPI.COMM_WORLD, ((0.0, 0.0), (1.0, .2)), (40, 5), mesh.CellType.triangle)
tdim = mesh_domain.topology.dim - 1

V = fem.functionspace(mesh_domain, ("P", 1, (2,)))

# Define trial and test functions
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# Define material properties
E = 210e9
nu = 0.3
lambda_ = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))

# Strain and stress
def epsilon(u):
    return ufl.sym(ufl.grad(u))

def sigma(u):
    return lambda_ * ufl.div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)

# Dirichlet BC (0 displacement on left boundary)
def boundary_points(x):   # return left boundary points
    return np.isclose(x[0], min(x[0]))

boundary_dofs = fem.locate_dofs_geometrical(V, boundary_points)

uD = fem.Function(V)
uD.interpolate(lambda x: np.zeros((2, x.shape[1])))     # zero values
bc = fem.dirichletbc(value=uD, dofs=boundary_dofs)

# Neumann BC (moving load)
class MovingLoad:
    """
    Load moves over top nodes, from ~0.5x to x.
    It uses the closest nodes.
    """
    def __init__(self):
        self.t = 0

    def eval(self, x):
        out = np.zeros((2, x.shape[1]))   # initialize output
        top_nodes = np.where(np.isclose(x[1], max(x[1])))     # get all top nodes
        x_coord = np.argmin(np.abs(x[0, top_nodes] - (0.5 + self.t / 20)))      # get node closest to 0.5+t/20
        out[1, top_nodes[0][x_coord]] = 10      # set constant load value for that specific node
        print(f"load pure {np.nonzero(out)}")
        # Out does contain nodes with the loading
        return out

load = MovingLoad()
V_load = fem.functionspace(mesh_domain, ("P", 1, (2,)))
load_func = fem.Function(V_load, name="load")
load_func.interpolate(load.eval)
# load_func is empty when interpolated
print(f"load interpolated {np.nonzero(load_func.x.array)}")

# Solving and plotting for steps
# weak form
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = -ufl.inner(load_func, v) * ufl.ds

problem = dolfinx.fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})

with XDMFFile(mesh_domain.comm, "MWE_results.xdmf", "w") as write_file:
    write_file.write_mesh(mesh_domain)
    for step in range(1, 11):
        load.t = step   # update load position
        load_func.interpolate(load.eval)

        uh = problem.solve()    # solve problem

        write_file.write_function(uh, step)   # write solution to file
        write_file.write_function(load_func, step)
2

You assumption about the second dimension of x being top nodes are wrong.

Consider the following mwe (that also works in parallel):

import numpy as np
import ufl
from mpi4py import MPI
import dolfinx
import dolfinx.fem.petsc


print(f"dolfinx version: {dolfinx.__version__}")

# Create a uniform 2D rectangular mesh
mesh_domain = dolfinx.mesh.create_rectangle(
    MPI.COMM_WORLD, ((0.0, 0.0), (1.0, 0.2)), (40, 5), dolfinx.mesh.CellType.triangle
)
tdim = mesh_domain.topology.dim - 1

V = dolfinx.fem.functionspace(mesh_domain, ("P", 1, (2,)))

# Define trial and test functions
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# Define material properties
E = 210e9
nu = 0.3
lambda_ = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))


# Strain and stress
def epsilon(u):
    return ufl.sym(ufl.grad(u))


def sigma(u):
    return lambda_ * ufl.div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)


# Dirichlet BC (0 displacement on left boundary)
def boundary_points(x):  # return left boundary points
    return np.isclose(x[0], 0.0)


boundary_dofs = dolfinx.fem.locate_dofs_geometrical(V, boundary_points)

uD = dolfinx.fem.Function(V)
uD.interpolate(lambda x: np.zeros((2, x.shape[1])))  # zero values
bc = dolfinx.fem.dirichletbc(value=uD, dofs=boundary_dofs)


def top_nodes(x):
    return np.isclose(x[1], 0.2)


# Neumann BC (moving load)
class MovingLoad:
    """
    Load moves over top nodes, from ~0.5x to x.
    It uses the closest nodes.
    """

    def __init__(self):
        self.t = 0

    def eval(self, x):
        out = np.zeros((2, x.shape[1]), dtype=dolfinx.default_scalar_type)  #
        top_loc = top_nodes(x)
        top_pos = np.flatnonzero(top_loc)

        # Compute distance from 0.5+t/20
        x_filtered = np.abs(x[0][top_loc] - (0.5 + self.t / 20))

        # Parallel safe min
        min_loc = np.min(x_filtered)
        min_glob = mesh_domain.comm.allreduce(min_loc, op=MPI.MIN)
        loc = np.flatnonzero(
            np.isclose(x_filtered, min_glob)
        )  # get all coordinates closest to 0.5+t/20

        out[1, top_pos[loc]] = 10  # set constant load value for that specific node
        return out


load = MovingLoad()
V_load = dolfinx.fem.functionspace(mesh_domain, ("P", 1, (2,)))
load_func = dolfinx.fem.Function(V_load, name="load")
load_func.interpolate(load.eval)
load_func.x.scatter_forward()


# Solving and plotting for steps
# weak form
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = -ufl.inner(load_func, v) * ufl.ds


problem = dolfinx.fem.petsc.LinearProblem(
    a,
    L,
    bcs=[bc],
    petsc_options={
        "ksp_type": "preonly",
        "pc_type": "lu",
        "pc_factor_mat_solver_type": "mumps",
    },
)

with dolfinx.io.XDMFFile(mesh_domain.comm, "MWE_results.xdmf", "w") as write_file:
    write_file.write_mesh(mesh_domain)
    for step in range(1, 11):
        load.t = step  # update load position
        load_func.interpolate(load.eval)
        load_func.x.scatter_forward()
        uh = problem.solve()  # solve problem

        write_file.write_function(uh, step)  # write solution to file
        write_file.write_function(load_func, step)
2 Likes
3

That works exactly how I want it, thank you dokken! I also appreciate the extra effort to make it work in parallel.