Non-trivial boundary conditions for linear elasticity

General
1

Hello,
I am trying to solve a linear elasticity problem where I already have a manufactured solution, and I am trying to impose it as a Dirichlet boundary condition on the whole boundary. However I keep getting errors of all sorts. Please find attached below the code.

E = 1.0
nu = 0.3
mu = E / (2.0 * (1.0 + nu))
lambda_ = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))

def epsilon(u):
        return ufl.sym(ufl.grad(u))

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

def u_ex(x):
    return np.array(0.4 * ((x[0] - 0.5) ** 2 + (x[1] -0.5) ** 2), dtype=np.double)
    
def f(x):
    val = -0.8 * (lambda_ + 3*mu) * np.ones_like(x, dtype=np.double)
    return val[:-1, :]

def assemble_lhs_rhs(N=100):

    comm = MPI.COMM_WORLD

    config.mesh_params["degree"] = 1
    config.mesh_params["N"] = N
    
    cartesian_mesh_params = {"comm": MPI.COMM_SELF, **config.mesh_params}
    cart_mesh = CartesianMesh(**cartesian_mesh_params)
    domain = cart_mesh.mesh

    vec = True
    V = fem.VectorFunctionSpace(domain, ("CG", config.mesh_params["degree"]))
    W = fem.FunctionSpace(domain, ("CG", config.mesh_params["degree"]))

    def clamped_boundary(x):
        return np.ones_like(x[1], dtype=bool)

    fdim = domain.topology.dim - 1
    boundary_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, clamped_boundary)

    f_dolfinx = dolfinx.fem.Function(V)
    f_dolfinx.interpolate(f)
    u_dolfinx = dolfinx.fem.Function(W)
    u_dolfinx.interpolate(u_ex)

    print(u_dolfinx.dtype)

    bc0 = fem.dirichletbc(u_dolfinx, fem.locate_dofs_topological(V, fdim, boundary_facets), V.sub(0))
    bc1 = fem.dirichletbc(u_dolfinx, fem.locate_dofs_topological(V, fdim, boundary_facets), V.sub(1))
    bc = [bc0, bc1]

I get an error when trying to declare bc0 and bc1, either the error says that the boundary condition has no dtype or it has the wrong one… however, I am declaring dtype=np.double above, which is why I am not sure what’s going on under the hood. Any help would be greatly appreciated.

Thanks!

2

Your code is not complete. Could you provide a minimal working example so that someone can reproduce your error message and help you?

3

Here is a minimal working example. I hope that this is more clear than the initial code I shared:

import dolfinx, ufl
import numpy as np
from mpi4py import MPI

E = 1.0
nu = 0.3
mu = E / (2.0 * (1.0 + nu))
lambda_ = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))

def epsilon(u):
        return ufl.sym(ufl.grad(u))

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

def u_ex(x):
    return np.array(0.4 * ((x[0] - 0.5) ** 2 + (x[1] -0.5) ** 2), dtype=np.double)
    
def f(x):
    val = -0.8 * (lambda_ + 3*mu) * np.ones_like(x, dtype=np.double)
    return val[:-1, :]

def assemble_lhs_rhs(N=100):

    comm = MPI.COMM_WORLD
    domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, N, N, dolfinx.mesh.CellType.quadrilateral)

    V = dolfinx.fem.VectorFunctionSpace(domain, ("CG", 1))
    W = dolfinx.fem.FunctionSpace(domain, ("CG", 1))

    def clamped_boundary(x):
        return np.ones_like(x[1], dtype=bool)

    fdim = domain.topology.dim - 1
    boundary_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, clamped_boundary)

    f_dolfinx = dolfinx.fem.Function(V)
    f_dolfinx.interpolate(f)
    u_dolfinx = dolfinx.fem.Function(W)
    u_dolfinx.interpolate(u_ex)

    print(u_dolfinx.dtype)

    bc0 = dolfinx.fem.dirichletbc(u_dolfinx, dolfinx.fem.locate_dofs_topological(V, fdim, boundary_facets), V.sub(0))
    bc1 = dolfinx.fem.dirichletbc(u_dolfinx, dolfinx.fem.locate_dofs_topological(V, fdim, boundary_facets), V.sub(1))
    bc = [bc0, bc1]

if __name__ == "__main__":
      assemble_lhs_rhs(N=25)

Best,

Ferhat

4

Could you try,

bc0 = dolfinx.fem.dirichletbc(u_dolfinx, dolfinx.fem.locate_dofs_topological(V.sub(0), fdim, boundary_facets))
bc1 = dolfinx.fem.dirichletbc(u_dolfinx, dolfinx.fem.locate_dofs_topological(V.sub(1), fdim, boundary_facets))

though, I am not quite sure if this is what you intended.