So, I tried to implement the Gonzalez scheme as follows, but at every time step, the residual is “not a number”… This might seem like a naïve/awkward situation, but I am not yet an advanced user of GetFEM, and I don’t know the correct way to update the variables in the time loop. The boundary conditions seem to be ok, on the other hand.
The situation is that of a 3D hyperelastic beam, characterized by the constitutive law given in the previous message, with a square section (of side 1) and axis aligned with the y-axis, of length L=4. The beam is clamped on the left end section (y=0) and loaded on the right end section (y=4) by a constant surface load parallel to the section, which is why I declared both f_old and f outside of the time loop (both are constant and equal to -10).
I made a grep in the GetFEM directory to try to find an example where analogous time loops are implemented, but every time I found “automated” things like Previous_Dot_u or md.shift_variables_for_time_integration. I would like to use less automated procedures, though, as described in the code below.
Is that possible? If yes, what am I missing?
Thanks for your help!
import numpy as np
import getfem as gf
def SetMeshRegionFromCondition(mesh,regionTag,condition):
pidCondition=np.compress(condition, range(mesh.nbpts()))
facesFromCondition=mesh.faces_from_pid(pidCondition)
mesh.set_region(regionTag,facesFromCondition)
gf.util('trace level', 1)
md = gf.Model("real")
E = 1e2 # Young's modulus [Pa]
nu = 0.49 # Poisson's ratio [-]
kappa = E/(3*(1-2*nu))
mu = E/(2*(1+nu))
rho = 500.0 # density
nbElems = 10
meshSize = 1/nbElems
fem_order = 2
# traction_force = 10 #Surfacic load on the Neumann boundary of the body in the vertical direction
# surfaceLoad=[0, 0, -traction_force]
# volumeLoad=[0, 0, 0]
m=gf.Mesh('regular simplices', np.arange(0,1+meshSize,meshSize),
np.arange(0,4+meshSize,meshSize),np.arange(0,1+meshSize,meshSize))
meshPoints=m.pts()
xCoord,yCoord,zCoord=meshPoints[0,:],meshPoints[1,:],meshPoints[2,:]
m.export_to_pos('Gonzalez.pos')
# Integration method
mim=gf.MeshIm(m, gf.Integ('IM_TETRAHEDRON(5)'))
# Create a MeshFem for u field of dimension 3 (i.e. a vector field)
mfu = gf.MeshFem(m,3) # Displacement
mfrhs = gf.MeshFem(m,3) # Right-hand side
# Choice of the finite element method
mfu.set_classical_fem(fem_order)
mfrhs.set_classical_fem(fem_order)
# mff=gf.MeshFem(mesh,1) # plot scalar Von-mises field
#Time related parameters
steps = 100
T_max = 1.0
dt = T_max/steps
#Dirichlet and Neumann conditions
#Boundaries definition
cleft=(abs(yCoord) < 1e-6)
cright=(abs(yCoord-4) < 1e-6)
#Boundaries creation
DIRICHLET_BOUNDARY = 1
NEUMANN_BOUNDARY = 2
SetMeshRegionFromCondition(mesh=m,regionTag=DIRICHLET_BOUNDARY,condition=cleft)
SetMeshRegionFromCondition(mesh=m,regionTag=NEUMANN_BOUNDARY,condition=cright)
md.add_fem_variable("u", mfu) # displacements
md.add_fem_variable("v", mfu) # velocities
md.add_fem_data("u_old", mfu)
md.add_fem_data("v_old", mfu)
# Initial conditions
U_0 = mfu.eval('[0, 0, 0]')
V_0 = 0.*U_0
md.set_variable('u_old', U_0)
md.set_variable('v_old', V_0)
f = mfrhs.eval('[0, 0, -10]')
f_old = mfrhs.eval('[0, 0, -10]')
md.add_initialized_data("K", E/(3*(1-2*nu)))
md.add_initialized_data("G", E/(2*(1+nu)))
md.add_macro("F", "Id(3)+Grad(u)")
md.add_macro("F_old", "Id(3)+Grad(u_old)")
md.add_macro("C", "F'*F") # "Right_Cauchy_Green(F)")
md.add_macro("C_old", "F_old'*F_old") # "Right_Cauchy_Green(F_old)")
md.add_macro("F_avg", "0.5*(F+F_old)")
md.add_macro("C_avg", "0.5*(C+C_old)")
md.add_macro("S(CC)", "K/2*log(Det(CC))*Inv(CC)+G*pow(Det(CC),-4/3)*(Id(3)-1/3*Trace(CC)*Inv(CC))")
md.add_macro("W(CC)", "K/8*sqr(log(Det(CC)))+G/2*(pow(Det(CC),-4/3)*(Trace(CC))-3)")
md.add_macro("SS", "S(C_avg)+(2*W(C)-2*W(C_old)-2*S(C_avg):(C-C_old))/Norm_sqr(C-C_old)*(C-C_old)")
# elasticity + inertia
md.add_initialized_data("rho", rho)
md.add_initialized_data("dt", dt)
md.add_linear_term(mim, "(u-u_old-dt/2*(v+v_old)).Test_u")
md.add_nonlinear_term(mim, "rho/dt*(v-v_old).Test_v + (F_avg*SS):Grad(Test_v)")
# boundary conditions
md.add_Dirichlet_condition_with_multipliers(mim, 'u', 0, DIRICHLET_BOUNDARY)
md.add_initialized_fem_data('f_old', mfrhs, f_old)
md.add_initialized_fem_data('f', mfrhs, f)
md.add_source_term_brick(mim, 'v', '(f+f_old)/2', NEUMANN_BOUNDARY)
# md.add_linear_term(mim, "((f+f_old)/2).Test_v", NEUMANN_BOUNDARY)
# Time loop
for timeStepIndex,timeStep in enumerate(np.arange(0.,T_max+dt,dt)):
print('Time step: %g' % timeStep)
md.solve('max_res', 1E-6, 'very noisy')
U = md.variable('u')
V = md.variable('v')
#Post processing of the solution
mfu.export_to_vtk('Displacement_%d.vtk' % timeStepIndex, mfu, U, 'Displacement')
# Update u_old and v_old
md.set_variable('u_old', U)
md.set_variable('v_old', V)