issue with mesh refinement on a subdomain

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
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 = 1
		self.rho = 1 = self.W/self.L
		self.gamma = 0.4***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)
		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.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 )


	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):
			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()

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  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »