### issue with mesh refinement on a subdomain

232
views
0
6 months ago by
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):

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
Could it just be a critical time step issue?  When I set
		# numerical parameters
self.MAX_ITER = 2000
self.dt = 0.01/2.0​
there is no divergence, but I didn't look more carefully at the solution.
written 6 months ago by David Kamensky