How to solve a nonlinear form? #287
Description
Hi, I am trying to solve a single step compressible neohookean problem. But I always get nans in the newton solver.
the full setup is here below. The print of debug steps show that the sol has all nan entries except the dirichlet boundary conditions.
The running command I am using is python simpleNeo.py -m beam-tet.mesh -o 1 -K 20 -mu 0.3 -l 0.0 -0.01 0.0, the beam-tet.mesh is from the data/
Thank you for the help.
import os
import mfem.ser as mfem
from mfem.ser import intArray
from os.path import expanduser, join, dirname
import numpy as np
def save_current_state(iteration, mesh_obj, displacement_gf, folder="results/"):
"""
Saves the current mesh and displacement solution to separate files.
The iteration number is used in place of time for the filename.
"""
os.makedirs(folder, exist_ok=True)
mesh_filename = os.path.join(folder, f"mesh_t{iteration}.vtk")
sol_filename = os.path.join(folder, f"sol_t{iteration}.npy")
mesh_obj.PrintVTK(mesh_filename)
np.save(sol_filename, np.array(displacement_gf.GetDataArray()))
print(f" Saved state for iteration {iteration} to '{folder}'")
def run(order=1,
meshfile="",
K=1.0,
mu=1.0,
dirichlet_selector=None,
neumann_selector=None,
load_vector=None):
mesh = mfem.Mesh(meshfile, 1, 1)
dim = mesh.Dimension()
assert mesh.SpaceDimension() == dim, "invalid mesh"
# --- Dynamic Boundary Condition Setup ---
# This section remains unchanged and correctly sets boundary attributes
if mesh.bdr_attributes.Size() > 0:
for i in range(mesh.GetNBE()):
mesh.SetBdrAttribute(i, 3) # Default to insulated
if dirichlet_selector is not None:
for i in range(mesh.GetNBE()):
vertices = mesh.GetBdrElement(i).GetVerticesArray()
center = np.mean([mesh.GetVertexArray(v_idx) for v_idx in vertices], axis=0)
if dirichlet_selector(center):
mesh.SetBdrAttribute(i, 1) # Dirichlet
if neumann_selector is not None:
for i in range(mesh.GetNBE()):
vertices = mesh.GetBdrElement(i).GetVerticesArray()
center = np.mean([mesh.GetVertexArray(v_idx) for v_idx in vertices], axis=0)
if neumann_selector(center):
mesh.SetBdrAttribute(i, 2) # Neumann
mesh.SetAttributes()
# --- End Dynamic BC Setup ---
fec = mfem.H1_FECollection(order, dim)
fespace = mfem.FiniteElementSpace(mesh, fec, dim)
print(f"\nNumber of unknowns: {fespace.GetTrueVSize()}")
# --- Setup RHS (Force Vector) for loaded Neumann BC ---
if load_vector is None:
load_vector = [0.0] * dim
if len(load_vector) != dim:
raise ValueError(f"Load vector dimension ({len(load_vector)}) must match mesh dimension ({dim})")
f_vec = mfem.Vector(load_vector)
f_coeff = mfem.VectorConstantCoefficient(f_vec)
b = mfem.LinearForm(fespace)
neumann_marker = mfem.intArray([0] * mesh.bdr_attributes.Max())
if neumann_marker.Size() > 1:
neumann_marker[1] = 1 # Mark attribute 2 for Neumann load
b.AddBoundaryIntegrator(mfem.VectorBoundaryLFIntegrator(f_coeff), neumann_marker)
b.Assemble()
# --- End RHS Setup ---
# --- Setup Neo-Hookean Operator ---
print(f"\nMaterial Properties:")
print(f" Bulk Modulus (K): {K}")
print(f" Shear Modulus (mu): {mu}")
# Define the displacement GridFunction, which is our unknown
x = mfem.GridFunction(fespace)
x.Assign(0.0) # Initial guess is zero displacement
# Define the essential (Dirichlet) boundary based on attribute 1
ess_bdr = mfem.intArray([0] * mesh.bdr_attributes.Max())
ess_bdr[0] = 1
# Define the nonlinear form H(x) for the hyperelastic model
H = mfem.NonlinearForm(fespace)
model = mfem.NeoHookeanModel(mu, K) # Note: mu is shear modulus, K is bulk modulus
H.AddDomainIntegrator(mfem.HyperelasticNLFIntegrator(model))
H.SetEssentialBC(ess_bdr)
# --- End Operator Setup ---
# --- Solve the Nonlinear System using Newton's Method ---
print("\nSolving the nonlinear system with Newton's method...")
# Create the Newton solver
newton_solver = mfem.NewtonSolver()
newton_solver.SetSolver(mfem.GMRESSolver()) # Use GMRES to solve the Jacobian system at each step
newton_solver.SetOperator(H)
newton_solver.SetPrintLevel(1) # Print Newton iteration details
newton_solver.SetRelTol(1e-8)
newton_solver.SetAbsTol(1e-8)
newton_solver.SetMaxIter(25)
# debug steps
print(np.array(b.GetDataArray()))
print(np.array(x.GetDataArray()))
sol = mfem.Vector(x.Size())
H.GetGradient(x).Mult(x, sol)
print(np.array(sol.GetDataArray()))
H.Mult(x, sol)
print(np.array(sol.GetDataArray()))
# end of debug steps
newton_solver.Mult(b, x)
# --- End Solve ---
# Save the final state
save_current_state(0, mesh, x)
print(f"\nSimulation complete. Final state saved.")
if __name__ == "__main__":
from mfem.common.arg_parser import ArgParser
parser = ArgParser(description='Compressible Neo-Hookean Elasticity Solver')
parser.add_argument('-m', '--mesh', default='beam-tet.mesh', action='store', type=str, help='Mesh file to use.')
parser.add_argument('-o', '--order', action='store', default=1, type=int, help="Finite element order (polynomial degree)")
parser.add_argument('-K', '--K', default=50.0, type=float, help="Bulk Modulus")
parser.add_argument('-mu', '--mu', default=20.0, type=float, help="Shear Modulus")
parser.add_argument('-l', '--load-vector', nargs='+', type=float, default=[0.0, 0.5, 0.0], help="Load vector components (e.g., fx fy fz)")
args = parser.parse_args()
parser.print_options(args)
# Ensure the mesh file path is correct
meshfile = expanduser(join(os.path.dirname(__file__), args.mesh))
# --- Define Boundary Selector Functions ---
def select_fixed_end(p):
return p[0] < 1e-4
def select_pull_end(p):
return p[0] > 8-1e-2 # The beam-tet mesh is 8 units long
# --- End Selector Functions ---
run(order=args.order,
meshfile=meshfile,
K=args.K,
mu=args.mu,
dirichlet_selector=select_fixed_end,
neumann_selector=select_pull_end,
load_vector=args.load_vector)