Suppose you are given a general nonlinear equation f(A,B) where A is a known vector-valued function and B is an unknown vector-valued function and want to solve for the latter variable throughout the domain in question. Also suppose that the equation does NOT involve derivatives, i.e. it’s purely algebraic. For our purposes also suppose that A and B are finite element functions defined on an appropriate function space, e.g. P1 elements.
I am wondering if there’s an efficient way to solve this in legacy fenics without calling a NonlinearVariationalProblem and requiring integration over the domain. I’ve read a bit about possibly writing a custom local NewtonSolver such that the nonlinear equation is solved at the integration points of the domain but haven’t successfully implemented that yet. I also wonder if it’s possible to simply extract the ufl functions into numpy arrays then solving using numpy/scipy infrastructure and then cast back into the function space for the solution.
I’ve attached a pseudocode MWE for hyperelasticity that’s unfinished but clarifies how this worflow would occur.
from dolfin import *
import ufl
from tqdm import tqdm
# ----------------------------------------------------------------
# Input parameters
# ----------------------------------------------------------------
# Material parameters
mu = 20.0e6 # Matrix shear modulus [Pa]
nu = 0.45 # Matrix Poisson's ratio
kappa = 2*mu*(1+nu)/(3*(1-2*nu)) # Matrix bulk modulus
lamda = 2*mu*nu/(1 - 2*nu) # Matrix first Lame Parameter
# Define the domain corners
x_min = 0.0
x_max = 1.0
y_min = 0.0
y_max = 0.05
# Define the number of divisions in the x and y directions
nx = 80
ny = 6
quad_degree = 4
# ----------------------------------------------------------------
# Read mesh
# ----------------------------------------------------------------
mesh = RectangleMesh(Point(x_min, y_min), Point(x_max, y_max), nx, ny)
dim = mesh.geometric_dimension()
mesh_coord = mesh.coordinates()
# ----------------------------------------------------------------
# Define function spaces
# ----------------------------------------------------------------
u_elem = VectorElement('CG', mesh.ufl_cell(), degree=2, dim=2)
V = FunctionSpace(mesh, u_elem)
# ----------------------------------------------------------------
# Define boundary conditions
# ----------------------------------------------------------------
left = CompiledSubDomain("near(x[0], mesh_xmin) && on_boundary", mesh_xmin = x_min)
right = CompiledSubDomain("near(x[0], mesh_xmax) && on_boundary", mesh_xmax = x_max)
# Dirichlet boundary
BC_left = DirichletBC(V, Constant((0.0,0.0)), left)
disp_right = Constant((0.0,0.0))
BC_right = DirichletBC(V,disp_right,right)
BC = [BC_left]
# ----------------------------------------------------------------
# Define stresses and strains
# ----------------------------------------------------------------
# Define test & trial spaces
u, eta, du = Function(V), TestFunction(V), TrialFunction(V)
# Kinematics
d = u.geometric_dimension()
I = Identity(d)
def F(u):
return I + grad(u)
def C(u):
return F(u).T*F(u)
def Ic(u):
return tr(C(u))
def J(u):
return det(F(u))
# Stored strain energy density (compressible neo-Hookean model)
psi = (((mu/2)*(Ic(u) - 3) - mu*ln(J(u)) + (lamda/2)*(ln(J(u)))**2))
dx = Measure("dx", domain=mesh)
Pi_int = psi*dx
residual= derivative(Pi_int,u,eta)
Jac = derivative(residual, u, du)
# Usual solver
problem = NonlinearVariationalProblem(residual, u, BC, Jac)
solver = NonlinearVariationalSolver(problem)
### ARBITRARY NONLINEAR EQUATION TO BE SOLVED ###
def NonlinearFxn(A,B):
# A and B are functions in the same finite element function space
return A^3 - inner(A,B)*B
for i in tqdm(range(0,20)):
disp_right.assign(Constant((i/10, 0.0)))
solver.solve()
### SOLVE NONLINEAR EQUATION NonlinearFxn(u,B)==0 HERE FOR B ###
B = Function(V)
MyNonlinearEqnSolve(NonlinearFxn(u,B)==0)