# SPDX-License-Identifier: MIT
# Copyright (c) 2011–2026 Joris J.C. Remmers
from pyfem.util.BaseModule import BaseModule
from numpy import zeros, array
from pyfem.fem.Assembly import assembleInternalForce, assembleTangentStiffness
from pyfem.fem.Assembly import assembleExternalForce, prepare, commit
from math import sin,cos,exp
import sys
from pyfem.util.logger import getLogger
logger = getLogger()
#------------------------------------------------------------------------------
#
#------------------------------------------------------------------------------
[docs]
class NonlinearSolver( BaseModule ):
def __init__( self , props , globdat ):
self.tol = 1.0e-3
self.iterMax = 10
self.maxCycle = sys.maxsize
self.maxLam = 1.0e20
self.dtime = 1.0
self.loadFunc = "t"
self.loadCases= []
BaseModule.__init__( self , props )
if self.maxLam > 1.0e19 and self.maxCycle == sys.maxsize:
self.maxCycle = 5
globdat.lam = 0.0
globdat.solverStatus.dtime = self.dtime
self.loadfunc = eval ( "lambda t : " + str(self.loadFunc) )
if hasattr(self,"loadTable"):
self.maxCycle = len(self.loadTable)
loadTable = zeros(self.maxCycle+1)
loadTable[1:] = self.loadTable
self.loadTable = loadTable
#------------------------------------------------------------------------------
#
#------------------------------------------------------------------------------
[docs]
def run( self , props , globdat ):
stat = globdat.solverStatus
stat.increaseStep()
self.writeHeader( stat.cycle )
dofCount = len(globdat.dofs)
a = globdat.state
Da = globdat.Dstate
Da[:] = zeros( dofCount )
fint = zeros( dofCount )
self.setLoadAndConstraints( globdat )
prepare( props , globdat )
K,fint = assembleTangentStiffness( props, globdat )
error = 1.
self.setLoadAndConstraints( globdat )
fext = assembleExternalForce ( props, globdat )
logger.info(' Newton-Raphson............ : L2-norm residual')
while error > self.tol:
stat.iiter += 1
da = globdat.dofs.solve( K, fext - fint )
Da[:] += da[:]
a [:] += da[:]
K,fint = assembleTangentStiffness( props, globdat )
# note that the code is different from the one presented in the book, which
# is slightly shorter for the sake of clarity.
# In the case of a prescribed displacement, the external force is zero
# and hence its norm is zero. In that case, the norm of the residue is not
# divided by the norm of the external force.
norm = globdat.dofs.norm( fext )
if norm < 1.0e-16:
error = globdat.dofs.norm( fext-fint )
else:
error = globdat.dofs.norm( fext-fint ) / norm
logger.info(f' Iteration {stat.iiter:4d} ........... : {error:6.4e}')
globdat.dofs.setConstrainFactor( 0.0 )
if stat.iiter == self.iterMax:
raise RuntimeError('Newton-Raphson iterations did not converge!')
# Converged
logger.info(' Converged')
globdat.elements.commitHistory()
Da[:] = zeros( len(globdat.dofs) )
globdat.fint = fint
commit ( props , globdat )
if stat.cycle == self.maxCycle or globdat.lam > self.maxLam:
globdat.active = False
self.writeFooter( globdat )
#-------------------------------------------------------------------------------
#
#-------------------------------------------------------------------------------
[docs]
def setLoadAndConstraints( self , globdat ):
if hasattr(self,"loadTable"):
cycle = globdat.solverStatus.cycle
globdat.lam = self.loadTable[cycle]
globdat.dlam = self.loadTable[cycle]-self.loadTable[cycle-1]
globdat.dofs.setConstrainFactor( globdat.dlam )
globdat.solverStatus.lam = globdat.lam
else:
globdat.lam = self.loadfunc( globdat.solverStatus.time )
lam0 = self.loadfunc( globdat.solverStatus.time - globdat.solverStatus.dtime )
globdat.dlam = globdat.lam - lam0
globdat.dofs.setConstrainFactor( globdat.dlam )
globdat.solverStatus.lam = globdat.lam
logger.debug(' ---- main load -------------------------')
logger.debug(f' loadFactor : {globdat.lam:4.2f}')
logger.debug(f' incr. loadFactor : {globdat.dlam:4.2f}')
for loadCase in self.loadCases:
loadProps = getattr( self.myProps, loadCase )
loadfunc = eval ( "lambda t : " + str(loadProps.loadFunc) )
lam = loadfunc( globdat.solverStatus.time )
lam0 = loadfunc( globdat.solverStatus.time - globdat.solverStatus.dtime )
dlam = lam - lam0
globdat.dofs.setConstrainFactor( dlam , loadProps.nodeTable )
logger.debug(f' ---- {loadCase} ---------------------')
logger.debug(f' loadFactor : {lam:4.2f}')
logger.debug(f' incr. loadFactor : {dlam:4.2f}')