# Varying Pointsources applied on a Nonlinear problem

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"]["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)
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

``````

I was able to reduce the computation time by doing the projection on the whole domain at once, then evaluating it at each node : (See the modification on the `fill_durface_defs()` function

``````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"]["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 another function
surface_def_left_func = fe.Function(V_lagrangian)

# 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
surface_def_left_func.assign(
fe.project(surface_def_left, V_lagrangian))

for i in range(len(self.P_left)):
self.surface_lagr_left[i] = surface_def_left_func(
self.coords_left[i, 0], self.coords_left[i, 1])

I = fe.Identity(2)
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

``````

Is there still a way to further optimize this ?

Thank you!