Hello,
I wanted to implement a hyperelastic solid problem where the external load are applied in the form of localised forces on a Neumann boundary (say e.g on the left face of a rectangular domain). Based on the fenics tutorial and the post on Pointsources on a Nonlinear Problem I was able to get the problem working.
The problem is that the forces are not constant or known a priori, I have an expression for, that depends on the domain deformation :
\pmb{G} = | (det \pmb{F})\pmb{F}^{-T} \pmb{N} | \pmb{g}
where \pmb{F} is the gradient of the deformation, \pmb{g} are the known forces, and \pmb{N} is the normal vector.
The way I could think of to implement that was to apply a fe.project() for each concerned node, then to apply a PointSource, on the fe.NonlinearProblem.F() function. It’s probably a very nonoptimal way to do this, could you guide me to a better approach ? It takes about 4 times the time to compute compared to a constant forces case. This is the code, I hope it’s not too long :
from ufl import nabla_div, cofac
import numpy as np
import fenics as fe
def clamped_boundary(x, on_boundary):
return on_boundary and abs(x[1]) < tol
def neumann_boundary(x, on_boundary):
return on_boundary and ((abs(x[1] - H) < tol) or abs(abs(x[0]-x_center) - W / 2) < tol)
tol = 1E-14
# Geometry and material properties
.
.
.
# create Mesh
n_x_Direction = 10
n_y_Direction = 30
mesh = fe.RectangleMesh(fe.Point(x_center - W / 2, 0),
fe.Point(x_center + W / 2, H), n_x_Direction, n_y_Direction)
h = fe.Constant(H / n_y_Direction)
# create Function Space(s)
V = fe.VectorFunctionSpace(mesh, 'Lagrange', 2)
V_lagrangian = fe.FunctionSpace(mesh, 'Lagrange', 2)
fe.parameters["form_compiler"]["quadrature_degree"] = 4
fe.parameters["form_compiler"]["cpp_optimize"] = True
# Trial and Test Functions
du = fe.TrialFunction(V)
v = fe.TestFunction(V)
u_np1 = fe.Function(V)
# functions
u_n = fe.Function(V)
N_left = fe.Constant((-1.0, 0.0))
# Define boudnaries
coupling_boundary = fe.AutoSubDomain(neumann_boundary)
fixed_boundary = fe.AutoSubDomain(clamped_boundary)
# clamp the beam at the bottom
bc = fe.DirichletBC(V, fe.Constant((0, 0)), fixed_boundary)
class Problem(fe.NonlinearProblem):
def __init__(self, J, F, bc):
self.bilinear_form = J
self.linear_form = F
self.bc = bc
self.P = None
self.coords = None
self.surface_lagr_left = []
fe.NonlinearProblem.__init__(self)
def F(self, b, x):
fe.assemble(self.linear_form, tensor=b)
self.fill_surface_defs(surface_def_left, V_lagrangian)
for i in range(self.coords_left.shape[0]):
pt = fe.Point(
(self.coords_left[i, 0], self.coords_left[i, 1]))
for j in range(2):
pt_src = fe.PointSource(
V.sub(j), pt, -self.P_left[i, j]*self.surface_lagr_left[i])
pt_src.apply(b)
self.bc.apply(b, x)
def J(self, A, x):
fe.assemble(self.bilinear_form, tensor=A)
self.bc.apply(A)
def fill_coords(self, data):
if self.coords is None:
self.coords = data.copy()
self.coords_left_ids = (self.coords[:, 0] == (x_center-W/2))
self.coords_left = self.coords[self.coords_left_ids, :]
self.surface_lagr_left = [None]*self.coords_left.shape[0]
def fill_forces(self, data):
self.P = data.copy()
self.P_left = self.P[self.coords_left_ids, :]
def fill_surface_defs(self, surface_defs, subspace):
surface_def_left = surface_defs
for i in range(len(self.P_left)):
self.surface_lagr_left[i] = fe.project(surface_def_left, subspace)(
self.coords_left[i, 0], self.coords_left[i, 1])
I = fe.Identity(2)
F = I + fe.grad(u_np1) # Deformation gradient
C = F.T*F # Right Cauchy-Green tensor
J = fe.det(F) # Determinant of deformation fradient
# Left face -------------
n_left = fe.dot(cofac(F), N_left)
surface_def_left = fe.sqrt(fe.inner(n_left, n_left))
psi = (aa*fe.inner(F, F) + ee - dd*fe.ln(J))*fe.dx - \
fe.dot(fe.Constant((0.0, 0.0)), u_np1)*fe.dx
Form = fe.derivative(psi, u_np1, v)
Jac = fe.derivative(Form, u_np1, du)
problem = Problem(Jac, Form, bc)
solver = fe.NewtonSolver()
# Set Newton solver options
prm = solver.parameters
prm["relative_tolerance"] = 1e-6
prm["absolute_tolerance"] = 1e-6
t = 0.0
n = 0
.
.
.
for i in range(n_pseudo_time):
u_np1.assign(u_n)
problem.fill_coords(coordinates_)
problem.fill_forces(forces_[:, i])
solver.solve(problem, u_np1.vector())
n += 1
t += dt
Thank you in advance !