Define DirichletBC for a Vector Function

mesh
23

Thanks very much. I have successfully change it to a function by using

x = ufl.SpatialCoordinate(mesh)
lambda_ = ufl.exp()
mu = ufl.exp()

But I also want to know whether I can use the form like

def lambda_(x):
    return

to define. I think for some situation, I can not write the expression of it directly.

24

Please make a proper mathematical statement of what kind of function you would like. In general i would suggest you to read through the tutorial i linked to above, as it covers most usecases

25

Hey dokken, sorry for bothering you again. I would have thought I have solved the problem with your help. But I meet problem again when I try to solve this linear elasticity problem:

image
image460×671 20.3 KB

I guess that the solution is
image
That solution satisfies the displacement boundary condition. Besides, we can calculate that
image
So, it also satisfies the balance equation and force boundary condition:
image

If this solution is right, the displacement at point(L,L) should be
image
But the result I get from FEniCSx is

array([ 0.01343142, -0.1       ])

Could you please help find whether there is any error in my calculation or code?

import numpy as np
from mpi4py import MPI
from ufl import Identity, Measure, TestFunction, TrialFunction, VectorElement, dot, dx, inner, grad, nabla_div, sym

from dolfinx.generation import RectangleMesh, UnitSquareMesh
from dolfinx.mesh import CellType, locate_entities_boundary
from dolfinx.fem import (Constant, DirichletBC, Function, LinearProblem, FunctionSpace, VectorFunctionSpace, 
                         locate_dofs_geometrical)

L = 1
pixel = 20

mesh = RectangleMesh(MPI.COMM_WORLD, np.array([[0, 0, 0], [L, L, 0]]), [pixel, pixel], cell_type=CellType.quadrilateral)
V = VectorFunctionSpace(mesh, ("CG", 1))

nu = 0.2
mu_m = 1
mu_f = 5*mu_m
import ufl
x = ufl.SpatialCoordinate(mesh)
mu = ufl.exp(mu_m + 4*x[0]*(1-x[0])*(mu_f - mu_m))
lambda_ = ufl.exp((mu_m + 4*x[0]*(1-x[0])*(mu_f - mu_m))*2*nu/(1-2*nu))

#define boundary conditions
Vy = V.sub(1).collapse()
uDy_high = Function(Vy)
with uDy_high.vector.localForm() as uDy_high_loc:
    uDy_high_loc.set(-0.1*L)
                                               
def highside(x):
    return np.isclose(x[1], L)
boundary_high = locate_dofs_geometrical((V.sub(1), Vy), highside)
bch = DirichletBC(uDy_high, boundary_high, V.sub(1))
                                               
uDy_low = Function(Vy)
with uDy_low.vector.localForm() as uDy_low_loc:
    uDy_low_loc.set(0)
                                               
def lowside(x):
    return np.isclose(x[1], 0)
boundary_low = locate_dofs_geometrical((V.sub(1), Vy), lowside)
bcl = DirichletBC(uDy_low, boundary_low, V.sub(1))

Vx = V.sub(0).collapse()
uDx = Function(Vx)
with uDx.vector.localForm() as uDx_loc:
    uDx_loc.set(0)
                                               
def origin(x):
    return np.logical_and(np.isclose(x[0], 0), np.isclose(x[1], 0))
boundary_origin = locate_dofs_geometrical((V.sub(0), Vx), origin)
bco = DirichletBC(uDx, boundary_origin, V.sub(0))

u_D = Function(V)
with u_D.vector.localForm() as loc:
    loc.set(0)
bco = DirichletBC(u_D, locate_dofs_geometrical(V, origin))

bcs = [bco, bcl, bch]
                                               
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)
from petsc4py import PETSc
f = Constant(mesh, PETSc.ScalarType((0, 0)))
T = Constant(mesh, PETSc.ScalarType((0, 0)))

ds = Measure("ds", domain=mesh)
a = inner(sigma(u), epsilon(v)) * dx
L = inner(f, v) * dx + inner(T, v) * ds

problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()                                        

import pyvista
import dolfinx.plot
# Start virtual framebuffer
pyvista.start_xvfb(wait=0.05)

# Create plotter and pyvista grid
p = pyvista.Plotter(title="Deflection", window_size=[800, 800])
topology, cell_types = dolfinx.plot.create_vtk_topology(mesh, mesh.topology.dim)
grid = pyvista.UnstructuredGrid(topology, cell_types, mesh.geometry.x)

# Attach vector values to grid and warp grid by vector
vals_2D = uh.compute_point_values().real 
vals = np.zeros((vals_2D.shape[0], 3))
vals[:,:2] = vals_2D
grid["u"] = vals
actor_0 = p.add_mesh(grid, style="wireframe", color="k")
warped = grid.warp_by_vector("u", factor=1)
actor_1 = p.add_mesh(warped, show_edges=True)
p.view_xy()
if not pyvista.OFF_SCREEN:
   p.show()
else:
    fig_array = p.screenshot(f"component.png") 

end = pixel + 1
vals_2D[end*end-1, :]

Thanks very much!

26

ufl.exp(x) is the exponential function e^x; therefore, in your code you have defined \mu(x) and \lambda(x) as the exponential of the functions presented in your problem statement, i.e.

\mu(x) = e^{1 + 4x(1-x)(4)} \\ \lambda(x) = e^{2\nu [1 + 4x(1-x)(4)] / (1-2\nu)}

Changing these to

mu = mu_m + 4*x[0]*(1-x[0])*(mu_f - mu_m)
lambda_ = (mu_m + 4*x[0]*(1-x[0])*(mu_f - mu_m))*2*nu/(1-2*nu)

reproduces the analytical solution.

1 Like
27

Oh my god, I think I have misunderstood the ufl.exp as expression, thanks very much!

1 Like