Hi all! I am trying to simulate an elastic wave propagating, starting from a point source and impacting on a circlular subdomain inside a square domain.

I have written the following code. The problem is that as you can see from the vizualization, the result does not seem to be correct. Do you have some recommendation? Thanks you very much in advance.

```
import mshr
from dolfin import *
from mshr import *
import matplotlib.pyplot as plt;
import os, sys, shutil
from os import path
solver = KrylovSolver("gmres")
solver.parameters["relative_tolerance"] = 5e-6
solver.parameters["maximum_iterations"] = 1000
solver.parameters["monitor_convergence"] = False
fig = plt.figure(figsize=(10,8), dpi=80)
X1, Y1 = -20,-20;
X3, Y3 = 20, 20
a_c = 2.5 #circle radius
x_c = 5
y_c = 5
domain = mshr.Rectangle(dolfin.Point(X1, Y1), dolfin.Point(X3, Y3))- mshr.Circle(dolfin.Point(x_c, y_c), a_c)
n_elem = 2**6
mesh = mshr.generate_mesh(domain, n_elem)
V=VectorFunctionSpace(mesh, "CG", 1)
# Time variables
Tfinal = 10
Ndt = 200
dt = Tfinal / Ndt
c = 1
def circle(x, on_boundary):
return on_boundary and near(pow(x[0],2) + pow(x[1],2) == a_c**2, DOLFIN_EPS)
# boundary conditions for the circle
bc_circle = DirichletBC(V, Constant((0.0,0.0)), circle)
# Material properties
E = 1
nu =1
rho = 1
lmbda = 1
mu = 2
def sigma(v):
return lmbda*div(u)*Identity(2) + 2*mu*(grad(u))
def epsilon(v):
return sym(grad(v))
# Previous and current solution
g_ini = Constant((0.0, 0.0))
u1= interpolate(g_ini, V)
u0= interpolate(g_ini, V)
# Variational problem at each time
u = TrialFunction(V)
v = TestFunction(V)
T = Constant((0,0))
a = inner(sigma(u), epsilon(v))*dx
L = 2*inner(u1,v)*dx-inner(u0,v)*dx + inner(T,v)*ds
u = Function(V);
t = 0;
while t <= Tfinal:
plt.cla()
bcs = [bc_circle]
A, b = assemble_system(a, L, bcs)
delta = PointSource(V, Point(-0.5, -1.5), sin(c * 3 * t))
delta.apply(b)
solver.solve(A, u.vector(), b)
u0.assign(u1)
u1.assign(u)
plot(u)
plt.pause(0.0001)
t += dt
```
```