Hi everyone.
I am using Sfepy for several weeks and i found it amazing.
I have some trouble to use it interactively:
I would like to run an elastodynamic calculation in interactive mode, but my code seems to give the steady state solution at each output
and i really do not know what is missing.
The configuration is a 2D plate subject to a pressure field on its top surface. I would like to see the dynamics of the deformation on short timescales.
Probably it is a basic mistake but i am really stuck with it.
Thank you very much for your help.
Regards.
Here is my code:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri May 7 15:40:38 2021
2D elastodynamics test-case:
A steel plate is subjet to a Gaussian pressure field centered on its top boundary
"""
from os.path import join as pjoin
import numpy as nm
from sfepy.discrete.fem import FEDomain, Field
from sfepy.mesh.mesh_generators import gen_block_mesh
from sfepy.discrete import (FieldVariable, Material, Integral,
Equation, Equations, Problem)
#from sfepy.mechanics.matcoefs import stiffness_from_lame
import sfepy.mechanics.matcoefs as mc
from sfepy.terms import Term
from sfepy.discrete.conditions import Conditions, EssentialBC
from sfepy.base.base import IndexedStruct
from sfepy.solvers.ls import ScipyDirect
from sfepy.solvers.nls import Newton
from sfepy.solvers.ts_solvers import SimpleTimeSteppingSolver
from sfepy.base.base import Struct
outputs_folder = "out_sfepy"
problem_name = "elastodynamics"
output_folder = pjoin(outputs_folder, problem_name)
output_format = "vtk"
#================================= Mesh ======================================
#----------- generation of a 2D block mesh:
Ndim=2
Ly = 2.0e-3 ; Lx = 10.0e-2
dims = [Lx, Ly]
shape = [1000, 20]
# Element size.
HX = Lx / (shape[0] - 1)
HY = Ly / (shape[1] - 1)
# generation of the mesh
mesh = gen_block_mesh(dims, shape, 0.5 * nm.array(dims),
name='user_block', verbose=False)
#=========================== Create a domain =================================
domain = FEDomain('domain', mesh)
#=============== Define volume "Omega" and boundaries "Gamma" ================
min_x, max_x = domain.get_mesh_bounding_box()[:, 0]
min_y, max_y = domain.get_mesh_bounding_box()[:, 1]
eps = 1e-8 * (max_x - min_x)
omega = domain.create_region('Omega', 'all')
gamma1 = domain.create_region('Gamma1','vertices in x < %.10f' % (min_x + eps),'facet')
gamma2 = domain.create_region('Gamma2','vertices in x > %.10f' % (max_x - eps),'facet')
gamma3 = domain.create_region('Gamma3','vertices in y > %.10f' % (max_y - eps),'facet')
gamma4 = domain.create_region('Gamma4','vertices in y < %.10f' % (min_y + eps),'facet')
#================= Define the finite element approximation ===================
# Using the field fu, we can define both the unknown variable  and the test variable .
field = Field.from_args('fu', nm.float64, 'vector', omega, approx_order=2)
u = FieldVariable('u', 'unknown', field)
du = FieldVariable('du', 'unknown', field)
ddu = FieldVariable('ddu', 'unknown', field)
v = FieldVariable('v', 'test', field, primary_var_name='u')
dv = FieldVariable('dv', 'test', field, primary_var_name='du')
ddv = FieldVariable('ddv', 'test', field, primary_var_name='ddu')
#============================ Material parameters ============================
# from E, nu :
E = 200.0e9 ; nu = 0.3
lam, mu = mc.lame_from_youngpoisson(E, nu, plane='strain')
#============================ For elastodynamics =============================
rho = 7800.0
# Longitudinal and shear wave propagation speeds.
cl = nm.sqrt((lam + 2.0 * mu) / rho)
cs = nm.sqrt(mu / rho)
#============================= Stiffness matrix ==============================
m = Material('m', D=mc.stiffness_from_lame(dim=2, lam=lam, mu=mu),values={'rho':rho})
#================================= Time ======================================
# Time-stepping parameters.
# Note: the Courant number C0 = dt * cl / H
dt = HX / cl # C0 = 1
print('dt=', dt)
timef = 1.5 * Lx / cl
#======================== Pressure distribution ==============================
def test_P(x, x0=0.0,R=1.0,P0=1.0):
#--- champ de pression sinusoidal:
#return P0*nm.sin(2.0*nm.pi*(x-x0)/R)
#--- pic de pression
return P0*nm.exp(-((x-x0)/R)**2.0)
def tractionLoad( ts, coors, mode='qp', equations=None, term=None, problem=None):
"""ts = TimeStepper, coor = coordinates of field nodes in region."""
if mode == 'qp':
val=nm.zeros((coors.shape[0], 3, 1),dtype=nm.float64)
for i in range(coors.shape[0]):
x=coors[i,0] #; y=coors[i,1]
pressure=test_P(x,0.5*dims[0], dims[0], 1.0e7)
val[i,:,0]=0.0, pressure, 0.0 # S11, S22, S12
return {'val' : val}
p = Material(name='load',kind='time-dependent',function=tractionLoad)
#=============================== Integration method ==========================
integral = Integral('i', order=3)
#========================== Definition of the "terms" ========================
t1 = Term.new('dw_lin_elastic(m.D, v, u)',integral, omega, m=m, v=v, u=u)
t2 = Term.new('dw_volume_dot(m.rho, ddv, ddu)',integral, omega, m=m, ddv=ddv, ddu=ddu)
#t3 = Term.new('dw_zero(dv, du)',integral, omega, m=m, dv=dv, du=du)
t4 = Term.new('dw_surface_ltr(load.val, v )',integral, gamma3, load=p,v=v)
eq = Equation('balance_of_forces in time',t1+t2+t4)
eqs = Equations([eq])
#============================ boundary conditions ============================
fix_u1 = EssentialBC('fix_u', gamma1, {'u.all' : 0.0})
fix_u2 = EssentialBC('fix_u', gamma2, {'u.all' : 0.0})
#=========================== Problem instance ================================
pb = Problem(problem_name,equations=eqs)
pb.set_bcs(ebcs=Conditions([fix_u1, fix_u2]))
pb.setup_output(output_dir=output_folder, output_format=output_format)
pb.update_materials()
#=============================== Solveur =====================================
ls = ScipyDirect({'use_presolve' : True})
nls_status = IndexedStruct()
nls = Newton({'i_max': 1, 'eps_a': 1e-6, 'eps_r' : 1e-6,}, lin_solver=ls, status=nls_status)
tss = SimpleTimeSteppingSolver({'t0' : 0.0, 't1' : 1.0e-9, 'n_step' : 10}, nls=nls, context=pb, verbose=True)
#================================ SOLVE ======================================
pb.set_solver(tss)
status = IndexedStruct()
vec = pb.solve(status=status)
print('Nonlinear solver status:\n', nls_status)
print('Stationary solver status:\n', status)
pb.save_state('linear_elasticity.vtk', vec)
