Hi
I've tried to solve elastic wave equation with RK timestepping:
from gryphon import ESDIRK
from dolfin import *
from mshr import *
d = 2
L = 20.0;
H = 5.0
geom = Rectangle(Point(0., 0.), Point(L, H))
hole = Circle(Point(L / 2., H / 2.), H / 4)
geom = geom - hole
mesh = Mesh()
mesh_generator = CSGCGALMeshGenerator2D()
mesh_generator.parameters["mesh_resolution"] = -1. #
mesh_generator.parameters["cell_size"] = 0.4 #
mesh_generator.generate(geom, mesh)
left = CompiledSubDomain("near(x[0], side) && on_boundary", side=0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side=20.0)
V = VectorFunctionSpace(mesh, "CG", 1)
# material properties
rho = 100.0
E, nu = 100.0, 0.3
mu, lmbd = Constant(E / (2 * (1 + nu))), Constant(E * nu / ((1 + nu) * (1 - 2 * nu)))
# variational formulation
u = TrialFunction(V)
v = TestFunction(V)
I = Identity(2)
def L(q):
return lmbd * tr(q) * I + mu * (q + q.T)
# mass matrix
m = rho * inner(u, v) * dx
# stiffness matrix
k = inner(L(grad(u)), grad(v)) * dx
# body forces
f = Expression(('0.0', '0.0'))
b = inner(f, v) * dx
rhs = -inner(L(grad(u)), grad(v)) * dx
# initial conditions
W = Function(V)
W = interpolate(Expression(('0', '0')), V)
# boundary conditions
gu = Expression(('sin(2*t)', '0.0'), t=0.0)
bcu1 = DirichletBC(V, gu, left)
bcu2 = DirichletBC(V, Expression(('0', '0')), right)
# timestepping parameters
T = [0,10]
print('Entering time step loop')
# Create the ESDIRK object.
obj = ESDIRK(T,W,rhs,bcs=[bcu1,bcu2])
# Call the solver which will do the actual calculation.
obj.solve()
interactive()
I get following error:
Invalid ranks 1 and 1 in product.
no matter if I use displacement or velocity-stress formulation. All of my problems are in VectorFunctionSpace. I think that's the reason.
I think I should invert the mass matrix (no forcing term) and multiply it by stifness matrix to obtain proper RHS.
