### issue with mesh refinement on a subdomain

232

views

0

I am fairly new to Fenics. I am running a simple dynamic elastic bar and I need mesh refinement in certain regions. However, once I refine the mesh then everything blows up. I am highly confident that this is not a numerical thing and it is due to the mesh refinement. Below is a simplified version of the code:

```
from fenics import *
import numpy as np
import matplotlib
matplotlib.use('TkAgg')
from matplotlib import pylab as plt
import mshr as msh
import numpy as np
import scipy.linalg as scla
def clamped_boundary(x, on_boundary):
return on_boundary and x[0] < DOLFIN_EPS
def epsilon(u):
return 0.5*(nabla_grad(u) + nabla_grad(u).T)
def sigma(u,lambda_,mu,d):
return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)
class Wave:
def __init__( self ):
# physical parameters
self.L = 1
self.W = 0.2
self.w = self.L/100
self.mu = 1
self.rho = 1
self.delta = self.W/self.L
self.gamma = 0.4*self.delta**2
self.beta = 1.25
self.lambda_ = self.beta
self.g = self.gamma
# numerical parameters
self.MAX_ITER = 1000
self.dt = 0.01
def initiate_fem( self ):
# generating mesh
geometry = msh.Rectangle(Point(0.0,0.0), Point(self.L,self.W))
mesh = msh.generate_mesh(geometry,10)
cell_markers = CellFunction("bool", mesh)
cell_markers.set_all(False)
for cell in cells(mesh):
p = cell.midpoint()
if p.x() < 0.7 and p.x() > 0.3:
cell_markers[cell] = True
# refining mesh
# mesh = refine(mesh, cell_markers)
self.mesh = mesh
#define function space
self.V = VectorFunctionSpace(self.mesh, 'P', 1)
#define dirichlet boundary
self.bc = DirichletBC(self.V, Constant((0, 0)), clamped_boundary)
#define right hand side function
self.f = Constant((0, -self.rho*self.g))
self.T = Constant((0, 0))
#define functions
q = TrialFunction(self.V)
self.d = q.geometric_dimension()
v = TestFunction(self.V)
aq = inner(q,v)*dx
kq = -inner(sigma(q,self.lambda_,self.mu,self.d),epsilon(v))*dx
Lq = inner(self.f,v)*dx
Kq, bq = assemble_system( kq, Lq, self.bc )
self.Aq, dummy = assemble_system( aq, Lq, self.bc )
#define the mass and stiffness matrices
M = np.matrix( self.Aq.array() )
self.M_inv = np.linalg.inv(M)
self.K = np.matrix( Kq.array() )
#define the force term
c = np.matrix( bq.array() )
self.cp = np.transpose( c )
print(self.cp.shape)
def stormer_verlet( self ):
vtkfile = File('results/solution.pvd')
f = Function(self.V)
N = self.cp.shape[0]
q0 = np.zeros([N,1])
p0 = np.zeros([N,1])
for i in range(0,self.MAX_ITER):
print(i)
q0 = q0 + self.dt/2*self.M_inv*p0
p0 = p0 + self.dt*self.K*q0 + self.dt*self.cp
q0 = q0 + self.dt/2*self.M_inv*p0
f.vector().set_local( q0 )
if np.mod(i,10) == 0:
vtkfile << (f,i*self.dt)
```

The mesh refinement is happening in the commented line. Without mesh refinement the solution is perfectly fine. But once you add the mesh refinement, the solution blows up. Since the refinement is happening away from the boundary, I doubt that this is any numerical effect. The solution will be save in the folder "./results" and you can view it in paraview. And here is a script to run the code:

```
wave = Wave()
wave.initiate_fem()
wave.stormer_verlet()
```

I would very much appreciate it if someone could point out what the problem is. Thanks

Community: FEniCS Project

Please login to add an answer/comment or follow this question.

there is no divergence, but I didn't look more carefully at the solution.