In Integration measures — FEniCS Tutorial @ Sorbonne the tags where generated when creating a mesh. I use a simpler definition where tags might not be necessary, but are still needed for the definition of the marker.
I create my mesh and subdomain like this:
domain = create_rectangle(
MPI.COMM_WORLD, # parallel communicator
[np.array([0, 0]), np.array([length, height])],
[Nx, Ny],
cell_type=CellType.quadrilateral,
)
top_facets = mesh.locate_entities_boundary(domain, fdim, lambda x: np.logical_and(np.isclose(x[1], height), x[0] <= length / 2))
bottom_facets = mesh.locate_entities_boundary(domain, fdim, lambda x: np.logical_and(np.isclose(x[1], 0), x[0] <= length / 2))
# Combine top and bottom markers into a single measure
friction_marker = mesh.meshtags(domain, fdim, np.concatenate([top_facets, bottom_facets]), 1)
ds_friction = Measure("ds", domain=domain, subdomain_data=friction_marker)
Applying this in my bilinear form I gives 3 different results depending on using ds_friction(1), ds_friction(0) or ds_friction, but I don’t understand why they (or at least two of them) are not equivalent.
a = (
inner(sigma(u), epsilon(v)) * dx +
(1 / dt) * inner(u, v) * dx +
friction_coefficient * inner(u, v) * ds_friction(1)
)
This is my full code:
import numpy as np # type: ignore
from ufl import sym, grad, Identity, tr, inner, Measure, TestFunction, TrialFunction # type: ignore
from mpi4py import MPI # type: ignore
from dolfinx import fem, io, mesh
from dolfinx.mesh import create_rectangle, CellType
from dolfinx.fem.petsc import assemble_vector, assemble_matrix, create_vector, apply_lifting, set_bc
from dolfinx.fem import locate_dofs_geometrical, DirichletBC
from dolfinx.geometry import compute_collisions_points, compute_colliding_cells, bb_tree
from petsc4py import PETSc # type: ignore
import matplotlib.pyplot as plt # type: ignore
from dolfinx.fem import form
# Define temporal parameters
t = 0 # Start time
T = 5.0 # Final time
timesteps = 10
dt = T / timesteps # time step size
# Domain parameters
length, height = 10, 1.0
Nx, Ny = 50, 10
# Create mesh
domain = create_rectangle(
MPI.COMM_WORLD, # parallel communicator
[np.array([0, 0]), np.array([length, height])],
[Nx, Ny],
cell_type=CellType.quadrilateral,
)
dim = domain.topology.dim
print(f"Mesh topology dimension d={dim}.")
degree = 2
shape = (dim,) # this means we want a vector field of size `dim`
# Define the function space
V = fem.functionspace(domain, ("P", degree, shape))
u = TrialFunction(V)
v = TestFunction(V)
u_h = fem.Function(V, name="Displacement")
u_n = fem.Function(V)
# make a fem function to store the force
internal_force = fem.Function(V, name="InternalForce")
# Lamé parameters (for a given material)
E = 10.0 # Young's modulus
nu = 0.3 # Poisson's ratio
lambda_ = (E * nu) / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))
# Define strain and stress tensors
def epsilon(u):
return sym(grad(u))
def sigma(u):
return lambda_ * tr(epsilon(u)) * Identity(len(u)) + 2 * mu * epsilon(u)
# Mark the right boundary with ID 1
fdim = domain.topology.dim - 1
right_facets = mesh.locate_entities_boundary(domain, fdim, lambda x: np.isclose(x[0], length))
right_marker = mesh.meshtags(domain, fdim, right_facets, 1)
# Define the boundary measure using the right marker
ds_right = Measure("ds", domain=domain, subdomain_data=right_marker)
# mark the left boundary with ID 2
left_facets = mesh.locate_entities_boundary(domain, fdim, lambda x: np.isclose(x[0], 0, 0.5))
left_marker = mesh.meshtags(domain, fdim, left_facets, 2)
# Define the boundary measure using the left marker
ds_left = Measure("ds", domain=domain, subdomain_data=left_marker)
# Weak formulation
f_x_value = 1.0 # initial magnitude of the force in the x-direction
friction_coefficient = 0.5
u, v = TrialFunction(V), TestFunction(V)
# force f is a constant over the whole domain
f = fem.Constant(domain, np.array([f_x_value, 0.0]))
# f2 is only applied in ds
#f2 = fem.Constant(ds, np.array([f_x_value, 0.0]))
dx = Measure("dx", domain=domain)
#a = inner(sigma(u), epsilon(v)) * dx
# billinear form
# describes internal stress and strain
# applied to the whole domain
#L = inner(f, v) * ds(1)
# linear form
# describes external forces
# applied to the right boundary ds(1)
# define the friction boundary
# Locate top and bottom facets for friction
top_facets = mesh.locate_entities_boundary(domain, fdim, lambda x: np.logical_and(np.isclose(x[1], height), x[0] <= length / 2))
bottom_facets = mesh.locate_entities_boundary(domain, fdim, lambda x: np.logical_and(np.isclose(x[1], 0), x[0] <= length / 2))
# Combine top and bottom markers into a single measure
friction_marker = mesh.meshtags(domain, fdim, np.concatenate([top_facets, bottom_facets]), 1)
ds_friction = Measure("ds", domain=domain, subdomain_data=friction_marker)
def friction_boundary(x):
# return True if the point is on the friction boundary
# top and bottom facets
return np.logical_and(np.logical_and(np.isclose(x[1], height), x[0] <= length / 2), np.logical_and(np.isclose(x[1], 0), x[0] <= length / 2))
# Create a function for the boundary condition (only in y-direction)
u_bc = fem.Function(V)
with u_bc.x.petsc_vec.localForm() as loc:
loc.set(0.0)
# Apply boundary condition only on y-component (index 1) at the friction boundary
boundary_dofs = locate_dofs_geometrical((V.sub(1), V), friction_boundary)
bcs = [
fem.dirichletbc(u_bc, boundary_dofs, V)
]
a = (
inner(sigma(u), epsilon(v)) * dx +
(1 / dt) * inner(u, v) * dx +
friction_coefficient * inner(u, v) * ds_friction(2)
)
L = inner(f, v) * ds_right(1) + (1 / dt) * inner(u_n, v) * dx
# https://jsdokken.com/FEniCS23-tutorial/src/benefits_of_curved_meshes.html
bilinear_form = fem.form(a)
linear_form = fem.form(L)
A = assemble_matrix(bilinear_form, bcs=bcs)
A.assemble()
b = create_vector(linear_form)
solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
# define the problem
problem = fem.petsc.LinearProblem(
a, L, u=u_h, bcs=bcs,
petsc_options={"ksp_type": "preonly", "pc_type": "lu"}
)
vtk = io.VTKFile(domain.comm, "results/dynamic/dynamic_cylinder_fric.pvd", "w")
displacements = []
internal_forces = []
bb_tree = bb_tree(domain, domain.topology.dim)
points = np.array([[0,0.5,0], [10,0.5,0]], dtype=domain.geometry.x.dtype)
potential_colliding_cells = compute_collisions_points(bb_tree, points)
colliding_cells = compute_colliding_cells(domain, potential_colliding_cells, points)
print(colliding_cells)
for i in range(timesteps):
print(f"Time step {i}/{timesteps}")
t += dt
# Update the right hand side reusing the initial vector
with b.localForm() as loc_b:
loc_b.set(0)
assemble_vector(b, linear_form)
# Apply Dirichlet boundary condition to the vector
apply_lifting(b, [bilinear_form], [bcs])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, bcs)
# Solve linear problem
solver.solve(b, u_h.x.petsc_vec)
u_h.x.scatter_forward()
# Update solution at previous time step (u_n)
u_n.x.array[:] = u_h.x.array
# Calculate the internal force/stress
stress = sigma(u_h)
force_form = form(inner(stress, epsilon(u_h)) * ds_left(2))
force_integral = fem.assemble_scalar(force_form)
print(f"Internal force: {force_integral}")
internal_forces.append(force_integral)
# Assign computed force to the function 'internal_force'
internal_force.x.array[:] = force_integral
# Write the computed internal force to file
vtk.write_function(internal_force, t)
# Write solution to file
vtk.write_function(u_h, t)
points_on_proc = []
cells = []
for i, point in enumerate(points):
if len(colliding_cells.links(i)) > 0:
points_on_proc.append(point)
cells.append(colliding_cells.links(i)[0])
points_on_proc = np.array(points_on_proc, dtype=np.float64).reshape(-1, 3)
cells = np.array(cells, dtype=np.int32)
u_values = u_h.eval(points_on_proc, cells)
displacements.append(u_values)
vtk.close()