Explicit Newmark in non linear mechanical problem

1

Dear All,
this is one of my first code in Fenics. I tried to write down the virtual work principle of an inextensible elastica (1D) subjected to viscous loads where a time-varying rotation is applied to one of the ends. The problem is non-linear and time-dependent so I would like to use explicit newmark method (not at all sure it’s well implemented).

This is the error that the compiler return:

Traceback (most recent call last):
  File "secondo.py", line 127, in <module>
    - (px(ux_upd(ux_n, vx_n, ax_n), uy_upd(uy_n,vy_n,ay_n), vx_upd(vx_n,ax_n,ax), vy_upd(vy_n,ay_n,ay))*ux_  \
TypeError: ux_upd() takes 3 positional arguments but 4 were given

And this is the code I wrote following essentially two tutorial

*https://fenicsproject.org/olddocs/dolfin/2019.1.0/python/demos/hyperelasticity/demo_hyperelasticity.py.html#equation-first-variation (non linear)

*https://fenicsproject.org/docs/dolfin/latest/python/demos/elastodynamics/demo_elastodynamics.py.html (time integration part)

If you have any suggestion or comment please let me know.

    from dolfin import *

    # Time-stepping parameters
    T = 1.0 				# final Time
    num_steps = 100 		# number of time steps
    dt = Constant(T / num_steps)		# time step size

    # GN-beta method
    beta=0.5										

    #  Material constants
    b = 1.0						# base
    l = 100.0					# length
    E = 1000.0					# Young's modulus
    rho = 1.0					# density
    I = Constant((b**4)/12)		# inertia
    B = Constant(E*I)			# bending stifness

    # Constraint constants
    A = 1.0				# amplitude
    om = 1.0/(10.0*pi)	# angular frequency

    # Viscous constants
    xit = 1.0
    xin = 5.0

    # Mesh and function space
    mesh = UnitIntervalMesh(10) 			# s \in [0,1]
    P1 = FiniteElement('P', interval , 2) 	# http://www.femtable.org/ ?hermitiane?
    el4 = MixedElement([P1, P1, P1, P1])	# mixed formulation for the 4 unknown fields 
    U = FunctionSpace(mesh, el4)
    el2 = MixedElement([P1, P1])
    V = FunctionSpace(mesh, el2)

    # Define Trial&Test functions
    alL_ = TestFunction(U) 
    u_ = TestFunction(V)

    #split trial&test functions to access the components
    ax_, ay_, lambd_, LAMBD_ = split(alL_)
    ux_, uy_ = split(u_)

    # Define function for displacement, velocities, accelerations and lagrangian multipliers
    u = Function(V)
    v = Function(V)
    alL = Function(U)
    u_n = Function(V)
    v_n = Function(V)
    alL_n = Function(U)

    # split functions to access the components
    ux, uy = split(u)
    vx, vy = split(v)
    ax, ay, lambd, LAMBD = split(alL)
    ux_n, uy_n = split(u_n)
    vx_n, vy_n = split(v_n)
    ax_n, ay_n, lambd_n, LAMBD_n = split(alL_n)

    # Define boundary 
    def left(x, on_boundary):
    	return near(x[0], 0) and on_boundary

    ################ time discretization ######################
    # Newmark formulas for updating displacement and velocity
    # Displacement
    def ux_upd(ux_n, vx_n, ax_n):
    	return ux_n + dt*vx_n + dt**2*(ax_n)/2

    def uy_upd(uy_n, vy_n, ay_n):
    	return uy_n + dt*vy_n + dt**2*(ay_n)/2

    # Velocity
    def vx_upd(vx_n, ax_n, ax):
    	return vx_n + dt*(1 - beta)*ax_n +beta*dt*ax

    def vy_upd(vy_n, ay_n, ay):
    	return vy_n + dt*(1 - beta)*ay_n +beta*dt*ay
    ################ time discretization ######################

    # Define rotation, curvature and constraints
    def theta(ux, uy):
    	return atan(div(uy)/(1+div(ux)))

    def chi(ux, uy):
    	return div(theta(ux, uy))

    def g(ux, uy):
    	return (1+div(ux))**2 + (div(uy))**2 - 1

    hh = Expression("A*sin(om*t)", A = A, om = om, t=0, degree=2)
    def h(ux, uy, hh):
    	return theta(ux, uy) - hh

    # Define the non conservative forces
    def px(ux, uy, vx, vy):
    	div(uy)*xin*(-div(uy)*vx + vy*(div(ux) + 1)) - xit*(div(ux) + 1)*(div(uy)*vy + vx*(div(ux) + 1))

    def px(ux, uy, vx, vy):
    	-div(uy)*xit*(div(uy)*vy + vx*(div(ux) + 1)) - xin*(div(ux) + 1)*(-div(uy)*vx + vy*(div(ux) + 1))

    # Perorm the variation through the derivative function 
    def dchi(ux, uy):
    	return derivative(chi(ux, uy), u, u_)

    def dg(ux, uy):
    	return derivative(g(ux, uy), u, u_)

    def dh(ux, uy):
    	return derivative(h(ux, uy, hh), u, u_)
    	
    # Define variational problem
    F = + rho*(ax*ax_ + ay*ay_)*dx \
    	- (B*chi(ux_upd(ux_n,vx_n,ax_n), uy_upd(uy_n,vy_n,ay_n))*dchi(ux_upd(ux_n,vx_n,ax_n), uy_upd(uy_n,vy_n,ay_n)))*dx \
    	+ g(ux_upd(ux_n,vx_n,ax_n), uy_upd(uy_n,vy_n,ay_n))*lambd_*dx + lambd*dg(ux_upd(ux_n,vx_n,ax_n), uy_upd(uy_n,vy_n,ay_n))*dx \
    	+ h(ux_upd(ux_n,vx_n,ax_n,hh), uy_upd(uy_n,vy_n,ay_n))*LAMBD_*ds + LAMBD*dh(ux_upd(ux_n,vx_n,ax_n), uy_upd(uy_n,vy_n,ay_n))*ds \
    	- (px(ux_upd(ux_n, vx_n, ax_n), uy_upd(uy_n,vy_n,ay_n), vx_upd(vx_n,ax_n,ax), vy_upd(vy_n,ay_n,ay))*ux_  \
    		+ py(ux_upd(ux_n,vx_n,ax_n), uy_upd(uy_n,vy_n,ay_n), vx_upd(vx_n,ax_n,ax), vy_upd(vy_n,ay_n,ay))*uy_)*dx

    # solver
    t = 0
    for n in range(num_steps):

    	# Update current time
    	t += dt

    	# Update the hh
    	hh.t = t

    	# Solve variational problem for time step 
    	solve(F == 0, alL)
        
    	# Update previous solution
    	u_n.assign(u)
    	v_n.assign(v)
    	alL_n.assign(alL)


    	print('time')

Thanks in advance to all the community

2

This is not related to dolfin as far as I can tell. You are calling a function ux_upd with 4 arguments whereas in the definition above it takes only 3.

3

Thanks you @bhaveshshrimali for your fast answer!

Now the problem is different:

Traceback (most recent call last):
  File "secondo.py", line 128, in <module>
    - (py(ux_upd(ux_n,vx_n,ax_n), uy_upd(uy_n,vy_n,ay_n), vx_upd(vx_n,ax_n,ax), vy_upd(vy_n,ay_n,ay))*uy_)*dx
TypeError: unsupported operand type(s) for *: 'NoneType' and 'Indexed'

code:

from ufl import *
from dolfin import *
from fenics import *

# Time-stepping parameters
T = 1.0 				# final Time
num_steps = 100 		# number of time steps
dt = Constant(T / num_steps)		# time step size

# GN-beta method
beta=0.5										

#  Material constants
b = 1.0						# base
l = 100.0					# length
E = 1000.0					# Young's modulus
rho = 1.0					# density
I = Constant((b**4)/12)		# inertia
B = Constant(E*I)			# bending stifness

# Constraint constants
A = 1.0				# amplitude
om = 1.0/(10.0*pi)	# angular frequency

# Viscous constants
xit = 1.0
xin = 5.0

# Mesh and function space
mesh = UnitIntervalMesh(10) 			# s \in [0,1]
P1 = FiniteElement('P', interval , 2) 	# http://www.femtable.org/ ?hermitiane?
el4 = MixedElement([P1, P1, P1, P1])	# mixed formulation for the 4 unknown fields 
U = FunctionSpace(mesh, el4)
el2 = MixedElement([P1, P1])
V = FunctionSpace(mesh, el2)

# Define Trial&Test functions
alL_ = TestFunction(U) 
u_ = TestFunction(V)

#split trial&test functions to access the components
ax_, ay_, lambd_, LAMBD_ = split(alL_)
ux_, uy_ = split(u_)

# Define function for displacement, velocities, accelerations and lagrangian multipliers
u = Function(V)
v = Function(V)
alL = Function(U)
u_n = Function(V)
v_n = Function(V)
alL_n = Function(U)

#split functions to access the components
ux, uy = split(u)
vx, vy = split(v)
ax, ay, lambd, LAMBD = split(alL)
ux_n, uy_n = split(u_n)
vx_n, vy_n = split(v_n)
ax_n, ay_n, lambd_n, LAMBD_n = split(alL_n)

# Define boundary 
def left(x, on_boundary):
	return near(x[0], 0) and on_boundary

################ time discretization ######################
# Newmark formulas for updating displacement and velocity
# Displacement
def ux_upd(ux_n, vx_n, ax_n):
	return ux_n + dt*vx_n + dt**2*(ax_n)/2

def uy_upd(uy_n, vy_n, ay_n):
	return uy_n + dt*vy_n + dt**2*(ay_n)/2

# Velocity
def vx_upd(vx_n, ax_n, ax):
	return vx_n + dt*(1 - beta)*ax_n +beta*dt*ax

def vy_upd(vy_n, ay_n, ay):
	return vy_n + dt*(1 - beta)*ay_n +beta*dt*ay
################ time discretization ######################

# Define rotation,curvature, constraints
def theta(ux, uy):
	return atan(div(uy)/(1+div(ux)))

def chi(ux, uy):
	return div(theta(ux, uy))

def g(ux, uy):
	return (1+div(ux))**2 + (div(uy))**2 - 1

hh = Expression("A*sin(om*t)", A = A, om = om, t=0, degree=2)
def h(ux, uy, hh):
	return theta(ux, uy) - hh

# Define the non conservative forces
def px(ux, uy, vx, vy):
	div(uy)*xin*(-div(uy)*vx + vy*(div(ux) + 1)) - xit*(div(ux) + 1)*(div(uy)*vy + vx*(div(ux) + 1))

def px(ux, uy, vx, vy):
	-div(uy)*xit*(div(uy)*vy + vx*(div(ux) + 1)) - xin*(div(ux) + 1)*(-div(uy)*vx + vy*(div(ux) + 1))

# Variation through the derivative function
def dchi(ux, uy):
	return derivative(chi(ux, uy), u, u_)

def dg(ux, uy):
	return derivative(g(ux, uy), u, u_)

def dh(ux, uy):
	return derivative(h(ux, uy, hh), u, u_)

F = + rho*(ax*ax_ + ay*ay_)*dx \
	- (B*chi(ux_upd(ux_n,vx_n,ax_n), uy_upd(uy_n,vy_n,ay_n))*dchi(ux_upd(ux_n,vx_n,ax_n), uy_upd(uy_n,vy_n,ay_n)))*dx \
	+ g(ux_upd(ux_n,vx_n,ax_n), uy_upd(uy_n,vy_n,ay_n))*lambd_*dx + lambd*dg(ux_upd(ux_n,vx_n,ax_n), uy_upd(uy_n,vy_n,ay_n))*dx \
	+ h(ux_upd(ux_n,vx_n,ax_n), uy_upd(uy_n,vy_n,ay_n),hh)*LAMBD_*ds + LAMBD*dh(ux_upd(ux_n,vx_n,ax_n), uy_upd(uy_n,vy_n,ay_n))*ds \
	- (px(ux_upd(ux_n, vx_n, ax_n), uy_upd(uy_n,vy_n,ay_n), vx_upd(vx_n,ax_n,ax), vy_upd(vy_n,ay_n,ay))*ux_)*dx \
	- (py(ux_upd(ux_n,vx_n,ax_n), uy_upd(uy_n,vy_n,ay_n), vx_upd(vx_n,ax_n,ax), vy_upd(vy_n,ay_n,ay))*uy_)*dx



# solver
t = 0
for n in range(num_steps):

	# Update current time
	t += dt

	# Update the hh
	hh.t = t

	# Solve variational problem for time step 
	solve(F == 0, alL)

	# Save solution to file (VTK)
    # _ux, _uy = u.split()
    # vtkfile_ux << (_ux, t)
    # vtkfile_uy << (_uy, t)
    
	# Update previous solution
	u_n.assign(u)
	v_n.assign(v)
	alL_n.assign(alL)


	print('time')

any hints?
thanks again

4

You aren’t returning anything in py also there are several typos in your code (for instance you have defined px twice instead of py). Try to carefully formulate the problem and then attempt solving it, that way you will be able to pin down the source of error faster. For the above error, try

def px(ux, uy, vx, vy):
    return div(uy)*xin*(-div(uy)*vx + vy*(div(ux) + 1)) - xit*(div(ux) + 1)*(div(uy)*vy + vx*(div(ux) + 1))

def py(ux, uy, vx, vy):
    return -div(uy)*xit*(div(uy)*vy + vx*(div(ux) + 1)) - xin*(div(ux) + 1)*(-div(uy)*vx + vy*(div(ux) + 1))

Also try to keep the test functions and functions distinct from each other.

1 Like
5

@bhaveshshrimali thank you again, my bad.

Unfortunately, I correct this typo but still, the terminal returns the same error as before

Regarding this:

are you referring in general or the suggestion is specific to this part?
Because in this definition

ux, uy, vx and vy are functions. I used the test functions only in here:

and in these definitions (u_)

Is that correct?
Thanks again for your patience