Boundary conditions definition dolfinx.fem.locate_dofs_geometrical

dolfinx
1

Hello,
hope you are good.
Im a student and new to fenicsx
And I have a few questions.

I want to set the x-component of the displacement to 0 on the left edge of a plate (2D geometry). For this, I’m using the script snippet below. But I got curious and wanted to print the contents of left_dofs_x.

here is what I got:

dtype=int32), array([505, 506, 516, 525, 533, 540, 546, 551, 555, 558, 560], dtype=int32)]

According to the documentation, locate_dofs_geometrical returns an array of degrees-of-freedom indices (local to the process) for degrees-of-freedom whose coordinate evaluates to True for the marker function.
If V is a list of two function spaces, then a 2-D array of shape (number of dofs, 2) is returned.

  • My first question is: don’t V.sub(0) and V_ux correspond to the same subspace?

  • My second question is: are the two columns (arrays) indices of degrees of freedom (local to the process) whose coordinate evaluates to True for the marker function left_edge, or does each array have a different meaning?
    I noticed that the indices of the first array are double the indices of the second array.

# Fix horizontal displacement on left edge (ux = 0)
V_ux, mapping = V.sub(0).collapse() # subspace of V for component ux

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

left_dofs_x = fem.locate_dofs_geometrical((V.sub(0), V_ux), left_edge)
print("left_dofs_x =", left_dofs_x)
u0_x = fem.Function(V_ux) 
bc_left_x = fem.dirichletbc(u0_x, left_dofs_x, V.sub(0))

Thanks in advance for taking your time to read this and for any reply

2

You’d have to supply a full working code for someone to reproduce your results, but to answer your questions:

They reference the same subspace, but are not the same representation of that space: V_ux ( the collapse’d one) is a stand-alone space whereas V.sub(0) is still part of V. That also explains

The coefficient numbers of V.sub(0) reference the dof in the V representation whereas those of V_ux reference the dofs in the stand-alone V_ux itself.

If your V is based on a mixed element with twice the same element, then that would explain why the coefficient numbers for V.sub(0) appear double that of V_ux

3

First of all, thanks for your reply.

i understood better the difference between V.sub(0) and V_ux

I don’t think that my V is based on a mixed element, well i’m not sure

Here is the full working code

### Linear elasticity - 2D Plate under pure elongation  

#------------- Import modules -------------------
import numpy as np
from ufl import sym, grad, Identity, tr, inner, Measure, TestFunction, TrialFunction, sqrt

from mpi4py import MPI

from dolfinx import fem, io
import dolfinx.fem.petsc
from dolfinx.mesh import create_rectangle, CellType, locate_entities_boundary, meshtags #MeshTags


#------------- Mesh and geometry --------------------
length, height = 5.0, 1.0
Nx, Ny = 50, 10                      # mesh refinement
domain = create_rectangle(
    MPI.COMM_WORLD,
    [np.array([0.0, 0.0]), np.array([length, height])],
    [Nx, Ny],
    cell_type=CellType.triangle,
)
## domain dimension
dim = domain.topology.dim
print(f"Mesh topology dimension d={dim}.")

#-------------- Function space (displacement field) and solution function -------------------
degree = 1 # degree of polynomial approximation for the displacement field (1=linear, 2=quadratic, 3=cubic)

shape = (dim,) # dim is the topological dimension of the domain (2 here)

V = fem.functionspace(domain, ("P", degree, shape))
u_sol = fem.Function(V, name="Displacement")


#-------------- Variational formulation -------------------

#-------------- Material parameters definition -------------------
# SI units: m, Pa, kg/m^3

E = fem.Constant(domain, 210e9)     # Young's modulus (Pa)
nu = fem.Constant(domain, 0.3)      # Poisson's ratio

lmbda = E * nu / (1 + nu) / (1 - 2 * nu)
mu = E / 2 / (1 + nu)



#-------------- Definition of strain and stress functions -------------------
def epsilon(v):
    return sym(grad(v))


def sigma(v):
    return lmbda * tr(epsilon(v)) * Identity(dim) + 2.0 * mu * epsilon(v)


# print("mu (UFL):\n", mu)
# print("epsilon (UFL):\n", epsilon(u_sol))
# print("sigma (UFL):\n", sigma(u_sol))

#-------------- Test and trial functions (real and virtual displacement fields) -------------------
u = TrialFunction(V)
v = TestFunction(V)

#-------------- Integration measure -------------------
dx = Measure("dx", domain=domain)

# There is no body force in this case so set f = 0
f  = fem.Constant(domain, np.array([0.0, 0.0]))

#-------------- Bilinear and linear forms -------------------
a = inner(sigma(u), epsilon(v)) * dx
L = inner(f, v) * dx  # Initialize L with the volume term f (here zero)


#-------------- Boundary conditions -------------------
# Numerically, we must remove rigid body modes (translations + rotations).
# Here we impose minimal Dirichlet constraints to obtain a unique solution:
#
# - Fix both components of displacement (ux = uy = 0) at the bottom-left corner (x = 0, y = 0) -> remove 2 translations
# - Fix horizontal displacement (ux = 0) on the left edge -> remove rotation


# Locate the degrees of freedom at the bottom-left corner 
def bottom_left_point(x):
    return np.isclose(x[0], 0.0) & np.isclose(x[1], 0.0)

bottom_left_point_dofs = fem.locate_dofs_geometrical(V, bottom_left_point)
u0_blp = fem.Function(V) # u0_blp is initialized to zero by default

u0_blp = np.zeros((2,), dtype=np.float64) # zero displacement vector at bottom-left corner
bc_blp = fem.dirichletbc(u0_blp, bottom_left_point_dofs, V) # Fix both components at bottom-left corner (x = 0, y = 0) 


# Fix horizontal displacement (ux = 0) on the left edge
V_ux, mapping = V.sub(0).collapse() # subspace of V for the ux component

def left_edge(x):
#    print("x shape: ", x.shape)
    return np.isclose(x[0], 0.0)

left_dofs_x = fem.locate_dofs_geometrical((V.sub(0), V_ux), left_edge)
print("left_dofs_x =", left_dofs_x)
u0_x = fem.Function(V_ux) # u0_x is initialized to zero by default
bc_left_x = fem.dirichletbc(u0_x, left_dofs_x, V.sub(0))


# Fix horizontal displacement (ux = epsilon*L) on the right edge
epsilon = 0.001  # prescribed uniaxial strain (0.1%)

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

right_dofs_x = fem.locate_dofs_geometrical((V.sub(0), V_ux), right_edge)
print("right_dofs_x =", right_dofs_x)

ulength_x = fem.Function(V_ux)

#Method 1
# or a fem.Expression
# ulength_x.interpolate(lambda x: epsilon * length * np.ones(x.shape[1], dtype=np.float64))
# print(f"ulength_x = {ulength_x} et le type est {type(ulength_x)}")

#Method 2
ulength_x.x.array[:] = epsilon * length

bc_right_x = fem.dirichletbc(ulength_x, right_dofs_x, V.sub(0))


# List of boundary conditions
bcs = [bc_blp, bc_left_x, bc_right_x]


#-------------- Assemble and solve linear system -------------------
petsc_options = {"ksp_type": "preonly", "pc_type": "lu"}

problem = fem.petsc.LinearProblem(a, L, u=u_sol, bcs=bcs,
                                  petsc_options=petsc_options)
problem.solve()

# ------------- Display displacement, strain, stress and von Mises -------------------
# spaces for tensorial / scalar fields
V_eps = fem.functionspace(domain, ("P", degree, (dim, dim)))   # strain (tensor)
V_sig = fem.functionspace(domain, ("P", degree, (dim, dim)))   # stress (tensor)
V_vm  = fem.functionspace(domain, ("P", degree, ()))          # Von Mises (scalar)

def project_field(expr, Vtarget, name_expr):
    """Interpolation via fem.Expression + interpolate"""
    interp = fem.Expression(expr, Vtarget.element.interpolation_points())
    F = fem.Function(Vtarget, name = name_expr)
    F.interpolate(interp)
    return F

# expressions for strain, stress and von Mises
eps_expr = sym(grad(u_sol)) #epsilon(u_sol)   
sig_expr = lmbda * tr(eps_expr) * Identity(dim) + 2.0 * mu * eps_expr
s_dev = sig_expr - (1.0/3.0) * tr(sig_expr) * Identity(dim)
vm_expr = sqrt(1.5 * inner(s_dev, s_dev))  

# project/interpolate
eps_proj = project_field(eps_expr, V_eps, "Déformation")
sig_proj = project_field(sig_expr, V_sig, "Contrainte")
vm_proj  = project_field(vm_expr, V_vm, "Von_Mises")

# .pvd output for visualization with Paraview

# 2
# with io.VTKFile(domain.comm, "results/linear_elasticity_2D_PureElongation/linear_elasticity_2D_PureElongation.pvd", "w") as vtk:
#     vtk.write_mesh(domain)
#     vtk.write_function(u_sol)
#     vtk.write_function(eps_proj)
#     vtk.write_function(sig_proj)
#     vtk.write_function(vm_proj)

# 1
vtk = io.VTKFile(domain.comm, "results/linear_elasticity_2D_PureElongation/linear_elasticity_2D_PureElongation.pvd", "w")
vtk.write_mesh(domain)
vtk.write_function(u_sol)
vtk.write_function(eps_proj)
vtk.write_function(sig_proj)
vtk.write_function(vm_proj)
vtk.close()

if domain.comm.rank == 0:
    print("Solution written! (open with ParaView).")

I used Jeremy Bleyer’s github
and Js Dokken website

Could you recommend (if any exist) any documents to better understand how FEniCSx works?
Thanks

4

Ah you’re right, that’s why a MWE was important :wink: Your V is a vector-valued, created by V = fem.functionspace(domain, ("P", degree, shape)) with shape=(2,). So it carries 2 degrees of freedom for each node, while V_ux is collapsed to only the x-components. Hence the doubling; dof n in V_ux corresponds to dof 2n in V (and, per my previous post, also to dof 2n in V.sub(0) )

The prominent resources are:

But those mainly show how to use FEniCS, not really how FEniCS works. The FEniCS book is outdated, but can still help you understand some core concepts https://launchpadlibrarian.net/83776282/fenics-book-2011-10-27-final.pdf

And then there are more technical resources, but that’s quite a big step https://orbilu.uni.lu/bitstream/10993/59327/1/dolfinx-paper-zenodo.10447666.pdf

1 Like
5

Thanks a lot M.Stein :grin: