How can I create a plot of the reaction forces on a mesh?
import pyvista
from dolfinx import mesh, fem, plot, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import numpy as np
import gmsh
from petsc4py import PETSc
import time
from matplotlib import pyplot as plt
pyvista.OFF_SCREEN = True
t1 = time.time()
def create_mesh(comm: MPI.Comm, model: gmsh.model, name: str, filename: str, mode: str):
"""Create a DOLFINx from a Gmsh model and output to file.
Args:
comm: MPI communicator top create the mesh on.
model: Gmsh model.
name: Name (identifier) of the mesh to add.
filename: XDMF filename.
mode: XDMF file mode. "w" (write) or "a" (append).
"""
msh, ct, ft = io.gmshio.model_to_mesh(model, comm, rank=0, gdim=2)
msh.name = name
ct.name = f"{msh.name}_cells"
ft.name = f"{msh.name}_facets"
msh.topology.create_connectivity(1, 2)
return msh, ct, ft
with io.XDMFFile(msh.comm, filename, mode) as file:
msh.topology.create_connectivity(1, 2)
file.write_mesh(msh)
file.write_meshtags(
ct, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry"
)
file.write_meshtags(
ft, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry"
)
def gmsh_wall(model: gmsh.model, name: str) -> gmsh.model:
"""Create a Gmsh model of a polygon with rectangular holes.
Args:
model: Gmsh model to add the mesh to.
name: Name (identifier) of the mesh to add.
Returns:
Gmsh model with a polygon and rectangular holes mesh added.
"""
model.add(name)
model.setCurrent(name)
# Define the points for the polygon
polygon_points = [
model.occ.addPoint(0, 0, 0),
model.occ.addPoint(10, 0, 0),
model.occ.addPoint(10, 5, 0),
model.occ.addPoint(0, 5, 0)
]
# Create lines for the polygon
polygon_lines = [
model.occ.addLine(polygon_points[0], polygon_points[1]),
model.occ.addLine(polygon_points[1], polygon_points[2]),
model.occ.addLine(polygon_points[2], polygon_points[3]),
model.occ.addLine(polygon_points[3], polygon_points[0])
]
# Create a curve loop and plane surface for the polygon
polygon_loop = model.occ.addCurveLoop(polygon_lines)
# Define the points for the rectangular holes
hole_points_1 = [
model.occ.addPoint(1, 1, 0),
model.occ.addPoint(6, 1, 0),
model.occ.addPoint(6, 2, 0),
model.occ.addPoint(1, 2, 0)
]
hole_points_2 = [
model.occ.addPoint(3, 3, 0),
model.occ.addPoint(4, 3, 0),
model.occ.addPoint(4, 4, 0),
model.occ.addPoint(3, 4, 0)
]
# Create lines for the rectangular holes
hole_lines_1 = [
model.occ.addLine(hole_points_1[0], hole_points_1[1]),
model.occ.addLine(hole_points_1[1], hole_points_1[2]),
model.occ.addLine(hole_points_1[2], hole_points_1[3]),
model.occ.addLine(hole_points_1[3], hole_points_1[0])
]
hole_lines_2 = [
model.occ.addLine(hole_points_2[0], hole_points_2[1]),
model.occ.addLine(hole_points_2[1], hole_points_2[2]),
model.occ.addLine(hole_points_2[2], hole_points_2[3]),
model.occ.addLine(hole_points_2[3], hole_points_2[0])
]
# Create curve loops for the holes
hole_loop_1 = model.occ.addCurveLoop(hole_lines_1)
hole_loop_2 = model.occ.addCurveLoop(hole_lines_2)
# Subtract the holes from the polygon surface
polygon_surface = model.occ.addPlaneSurface([polygon_loop, hole_loop_1, hole_loop_2])
model.occ.synchronize()
# Add physical tags for the circumference of each hole
model.addPhysicalGroup(1, [hole_lines_1[2]], tag=2)
model.setPhysicalName(1, 2, "gd1")
model.addPhysicalGroup(1, [hole_lines_2[2]], tag=3)
model.setPhysicalName(1, 3, "gd2")
#add physical tags for the lower boundery
model.addPhysicalGroup(1, [polygon_lines[0]], tag=4)
model.setPhysicalName(1, 4, "floor")
model.addPhysicalGroup(2, [polygon_surface], tag=5)
model.setPhysicalName(2, 5, "wall")
model.mesh.setOrder(1)
model.mesh.generate(dim=2)
return model
gmsh.initialize()
#gmsh.option.setNumber("General.Terminal", 0)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.1)
# Create model
model = gmsh.model()
model = gmsh_wall(model, "wall")
model.setCurrent("wall")
msh, ct, ft = create_mesh(MPI.COMM_WORLD, model, "wall", "out_gmsh/mesh.xdmf", "w")
yieldstrength = 8273709.0 # Yield strength in N/m^2
E = fem.Constant(msh, 700 * yieldstrength) # Young's modulus in N/m^2
nu = fem.Constant(msh, 0.156) # Poisson's ratio
density = 1800 # Density in kg/m^3
thickness = fem.Constant(msh, 0.1) # Thickness in meters (10 cm)
g = 10 # Gravitational acceleration in m/s^2
k_spring = 4e3 / 0.004 # Spring stiffness in N/m
lmbda = E*nu/(1+nu)/(1-2*nu) # Lame's first parameter
mu = E/2/(1+nu) # Lame's second parameter
sig0 = fem.Constant(msh, 15.1e6) # yield strength in N/m^2
Et = E / 100.0 # tangent modulus
H = E * Et / (E - Et) # hardening modulus
# function space
V = fem.functionspace(msh, ("Lagrange", 1, (2, )))
du = ufl.TrialFunction(V)
u_ = ufl.TestFunction(V)
ds = ufl.Measure("ds", domain=msh, subdomain_data=ft)
# boundary conditions
class BoundaryCondition():
def __init__(self, type, marker, values):
self._type = type
if type == "Dirichlet":
facets = ft.find(marker)
dofs = fem.locate_dofs_topological(V, 1, facets)
#self._bc = fem.dirichletbc(values, dofs, V)
self._bc = fem.dirichletbc(values, dofs)
elif type == "Neumann":
self._bc = ufl.inner(values, u_) * ds(marker)
elif type == "Robin":
self._bc = values[0] * ufl.inner(du - values[1], u_)* ds(marker)
else:
raise TypeError("Unknown boundary condition: {0:s}".format(type))
@property
def bc(self):
return self._bc
@property
def type(self):
return self._type
# Define the Dirichlet boundary condition
u_bc = fem.Function(V)
boundary_conditions = [BoundaryCondition("Dirichlet", 4, u_bc),
BoundaryCondition("Dirichlet", 2, u_bc),
BoundaryCondition("Dirichlet", 3, u_bc),
#BoundaryCondition("Robin", 2, (k_spring, u_bc)),
#BoundaryCondition("Robin", 3, (k_spring, u_bc)),
]
bcs = []
F = 0
for condition in boundary_conditions:
if condition.type == "Dirichlet":
bcs.append(condition.bc)
else:
F += condition.bc
# Body force (weight)
f = fem.Constant(msh, default_scalar_type((0, -density * g)))
# Traction force
T = fem.Constant(msh, default_scalar_type((0, 0)))
# Define stress and strain
def epsilon(u):
return ufl.sym(ufl.grad(u))
def sigma(u):
return lmbda * ufl.tr(epsilon(u)) * ufl.Identity(msh.topology.dim) + 2 * mu * epsilon(u)
# Define variational problem
a = thickness * ufl.inner(sigma(du), epsilon(u_)) * ufl.dx #+ ufl.lhs(F)
L = thickness * ufl.inner(f, u_) * ufl.dx #+ ufl.rhs(F)
u = fem.Function(V, name="Displacement")
# Solve the problem
problem = fem.petsc.LinearProblem(a, L, bcs, u=u)
u_h = problem.solve()
# Compute the reaction forces on the Dirichlet boundary
residual = ufl.action(a, u) - L
v_reac = fem.Function(V)
virtual_work_form = fem.form(ufl.action(residual, v_reac))
def one(x):
values = np.zeros((1, x.shape[1]))
values[0] = 1.0
return values
def y(x):
values = np.zeros((1, x.shape[1]))
values[0] = x[1]
return values
u_bc.sub(1).interpolate(one)
fem.set_bc(v_reac.vector, bcs)
print(f"Vertical reaction Ry = {fem.assemble_scalar(virtual_work_form):.6f}")