In dolfinx, how to apply a fixed displacement condition to a vertex?

dolfinx
1

I am solving the thermoelastic equation and applying a fixed displacement condition at the origin of a rectangle. Currently, I can only achieve this in dolfin, and the code is as follows：

from dolfin import *

mesh = RectangleMesh(Point(0., 0.), Point(1, 1), 10, 10)

class Gamma(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0, DOLFIN_EPS) and near(x[1], 0, DOLFIN_EPS)
Center = Gamma()

E = Constant(50e3)
nu = Constant(0.2)
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
alpha = Constant(1e-5)

def eps(v):
    return 0.5 * (grad(v) + grad(v).T)
def sigma(v, dT):
    return (lmbda*tr(eps(v))- alpha*(3*lmbda+2*mu)*dT)*Identity(2) + 2.0*mu*eps(v)

VU = VectorFunctionSpace(mesh, 'CG', 1)
dU = TrialFunction(VU)
U_ = TestFunction(VU)
Wint = inner(sigma(dU, 1000), eps(U_)) * dx
aM = lhs(Wint)
LM = rhs(Wint)

bcu = DirichletBC(VU, Constant((0, 0)), Center, method="pointwise")
U = Function(VU)
solve(aM == LM, U, bcu)

I’m not sure how to implement the same functionality in dolfinx. Here is the code snippet in dolfinx where I can only impose a fixed displacement condition on a specific boundary：

import numpy as np
from ufl import (tr, grad, Identity, TrialFunction, TestFunction, lhs, rhs, inner, dx)
from dolfinx import (mesh, io, default_scalar_type, fem)
from dolfinx.fem import functionspace, dirichletbc
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
from pathlib import Path
from dolfinx.mesh import create_unit_square

domain = create_unit_square(MPI.COMM_WORLD, 10, 10)

def left(x):
    return np.isclose(x[0], 0)

E = 50e3
nu = 0.2
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
alpha = 1e-5

def eps(v):
    return 0.5 * (grad(v) + grad(v).T)
def sigma(v, dT):
    return (lmbda*tr(eps(v))- alpha*(3*lmbda+2*mu)*dT)*Identity(2) + 2.0*mu*eps(v)

VU = functionspace(domain, ("CG", 1, (domain.geometry.dim, )))
dU = TrialFunction(VU)
U_ = TestFunction(VU)
Wint = inner(sigma(dU, 1000), eps(U_)) * dx
aM = lhs(Wint)
LM = rhs(Wint)

fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, left)
u_D = np.array([0, 0], dtype = default_scalar_type)
bc = dirichletbc(u_D, fem.locate_dofs_topological(VU, fdim, boundary_facets), VU)

# Compute solution
problem = LinearProblem(aM, LM, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
Uf = problem.solve()

I hope someone can help me modify this. Thank you.

2

Replace fdim with 0.
i.e. instead of boundary facets, you want the boundary_vertices.

3

I have made the modifications as you instructed, but maybe I should change the definition of “left” to “vertex.” I am not sure how to make this modification.

def left(x):
    return np.isclose(x[0], 0)

fdim = 0
boundary_vertices = mesh.locate_entities_boundary(domain, fdim, left)
print(boundary_vertices)
u_D = np.array([0, 0], dtype = default_scalar_type)
bc = dirichletbc(u_D, fem.locate_dofs_topological(VU, fdim, boundary_vertices), VU)
4

Currently you are looking for all vertices that have x=0.
You want something like

def center(x):
    return np.isclose(x[0], 0) & np.isclose(x[1], 0)

if you want the point (0,0)

5

Thank you for your help. I really didn’t know how to use isclose to define a point, now it’s working well.
I have another question. In the diagram, the red square represents a grid cell. If I want to define boundary conditions at the origin, does the origin have to be located at a vertex (as shown in picture2), or is it unnecessary (as shown in picture1)? I’m more inclined towards picture2.
a1a6728ef7473540bc3e5e0d8f79e2f

6

If it is not defined at a vertex, it is not a condition that you can enforce strongly (without using DOLFINx_MPC), as the condition then becomes an equation \sum_i c_i\phi_i(x) for all i with basis functions with support around x.

If possible, I advice you to align the mesh such that the condition is at the point of a degree of freedom.

7

Indeed, that’s correct. Thank you for your patient explanation.