Sorry for the long post, the code I am using is from an article to learn how to use fenics and paraview, but I hope this gives enough information to reproduce the issue
from dolfin import *
import numpy as np
import os
import dolfin
#-------------------------------------------------------------
# Parameters
#-------------------------------------------------------------
INFO = 20
set_log_level(INFO)
# Optimization options for the form compiler
#parameters["num_threads"] = 1
parameters["mesh_partitioner"] = "SCOTCH"
parameters["form_compiler"]["quadrature_degree"] = 2
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["log_level"] = INFO
parameters["allow_extrapolation"] = True
ffc_options = {"optimize": True, \
"eliminate_zeros": True, \
"precompute_basis_const": True, \
"precompute_ip_const": True, \
"quadrature_degree": 2}
# Some user parameters
user_par = Parameters("user")
user_par.add("bounds_xmin",-0.5)
user_par.add("bounds_xmax",0.5)
user_par.add("bounds_ymin",-0.5)
user_par.add("bounds_ymax",0.5)
user_par.add("bounds_zmin",0.)
user_par.add("bounds_zmax", 2.5)
user_par.add("fe_order_u",1)
user_par.add("fe_order_p",1)
user_par.add("gamma_min",0.)
user_par.add("gamma_max",.2)
user_par.add("gamma_nsteps",10)
user_par.add("mesh_ref",10)
user_par.add("save_dir","results")
user_par.add("output_type","pvd")
user_par.add("plot",True)
# Non linear solver parameters
solver_par = NonlinearVariationalSolver.default_parameters()
solver_par.rename("solver")
solver_par["symmetric"]=True
solver_par["linear_solver"]= 'mumps' #use "mumps" in parallel
solver_par["lu_solver"]["same_nonzero_pattern"] = True
solver_par["lu_solver"]["verbose"] = True
solver_par["newton_solver"]["maximum_iterations"] = 20
solver_par["newton_solver"]["relaxation_parameter"] = .8
solver_par["newton_solver"]["relative_tolerance"] = 1e-5
solver_par["newton_solver"]["absolute_tolerance"] = 1e-5
# add user parameters in the global parameter set
parameters.add(user_par)
parameters.add(solver_par)
# Parse parameters from command line
parameters.parse()
info(parameters,True)
user_par = parameters.user
#-------------------------------------------------------------
# Geometry
#-------------------------------------------------------------
# Create the geometry and the mesh
xmin,xmax = user_par.bounds_xmin,user_par.bounds_xmax
ymin,ymax = user_par.bounds_ymin,user_par.bounds_ymax
zmin,zmax = user_par.bounds_zmin,user_par.bounds_zmax
geom = Box(xmin,ymin,zmin,xmax,ymax,zmax)
mesh = Mesh(geom,user_par.mesh_ref)
#-------------------------------------------------------------
# Definition of function spaces
#-------------------------------------------------------------
# Create function space
P2 = VectorFunctionSpace(mesh, "CG", user_par.fe_order_u) # Space for displacement
P1 = FunctionSpace(mesh, "CG", user_par.fe_order_p) # Space for pressure
V = MixedFunctionSpace([P1,P2])
V_u = V.sub(1)
V_p = V.sub(0)
ndim = P2.cell().d
# Create functions to define the energy and store the results
up = Function(V)
(p,u)=split(up)
# Create test and trial functions for the variational formulation
dup = TrialFunction(V)
vq = TestFunction(V)
(q,v) = TestFunctions(V)
#-------------------------------------------------------------
# Boundary conditions
#-------------------------------------------------------------
# Mark boundary subdomains
xtol = mesh.hmin()/4.
class ClampedBoundary(SubDomain):
def inside(self, x, on_boundary):
return x[2]-zmin<xtol and on_boundary
# Define the boundary conditions
zero_vector = Constant((0.0,0.0,0.0))
bc1 = DirichletBC(V_u, zero_vector, ClampedBoundary())
bc_u = [bc1]
#-------------------------------------------------------------
# Define boundaries with surface tension
#-------------------------------------------------------------
# Define the part on the boundary where surface tension should be applied
class SurfaceBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary
# Mark facets where apply surface tension with 1
boundary_parts = FacetFunction("size_t", mesh, 0)
surface_boundary = SurfaceBoundary()
surface_boundary.mark(boundary_parts, 1)
# Redefine element of area to include informations about surface tension
ds = ds[boundary_parts]
#-------------------------------------------------------------
# Kinematics
#-------------------------------------------------------------
I = Identity(ndim) # Identity tensor
F = I + grad(u) # Deformation gradient
C = transpose(F)*F # Right Cauchy-Green tensor
E = 0.5*( C - I ) # Green-Lagrange tensor
# Invariants of deformation tensors
Ic = tr(C)
J = det(F)
# Normal and tangent vectors in the reference configuration
N = FacetNormal(mesh)
# Element of area transformation operator
NansonOp = transpose(cofac(F))
# surface element vector in the deformed configuration
deformed_N = dot(NansonOp,N)
# norm of the surface element vector in the current configuration
current_element_of_area = sqrt(dot(deformed_N,deformed_N))
#-------------------------------------------------------------
# Energy and variational formulation
#-------------------------------------------------------------
# Lame’s parameters
mu, lmbda = Constant(1.), Constant(1000.)
# Bulk energy (strain energy for an almost incompressible neo-Hookean model)
bulk_energy_density = mu*(Ic - ndim) -(mu+ p)*ln(J) - 1/(2*lmbda)*p**2
bulk_energy = bulk_energy_density*dx
# Surface energy
gamma=Expression("t",t=0.00)
surface_energy_density = gamma*current_element_of_area
surface_energy = surface_energy_density*ds(1)
# Total potential energy
potential_energy = bulk_energy + surface_energy
# First directional derivative of the potential energy (a linear form in the test function vq)
F=derivative(potential_energy,up,vq)
# First directional derivative of the potential energy (a bilinear form in the test function vq and the trial function dup)
dF=derivative(F,up,dup)
# Setup the variational problem
varproblem = NonlinearVariationalProblem(F, up, bc_u, J=dF,form_compiler_parameters=ffc_options)
#-------------------------------------------------------------
# # Set up the solver (Newton solver)
#-------------------------------------------------------------
solver = NonlinearVariationalSolver(varproblem)
solver.parameters.update(parameters.solver)
#-------------------------------------------------------------
# # Solve the problem
#-------------------------------------------------------------
# loading parameter (list of values for the surface tension)
gamma_list = np.linspace(user_par.gamma_min,user_par.gamma_max,user_par.gamma_nsteps) # list of values of surface tension for the simulations
# directory and files to save the results
save_dir = parameters.user.save_dir
file_u = File(save_dir+"/displacement."+parameters.user.output_type)
file_p = File(save_dir+"/pressure."+parameters.user.output_type)
# Solve with Newton solver for each value of the surface tension, using the previous solution as a starting point.
for t in gamma_list:
# update the value of the surface tension
gamma.t = t
# solve the nonlinear problem (using Newton solver)
solver.solve()
# Save solution to file (readable by Paraview or Visit)
(p,u) = up.split()
file_u << (u,t)
file_p << (p,t)
# Plot and save png image
if parameters.user.plot:
plot_u = plot(u, mode = "displacement",title="Displacement field gamma=%.4f"%t,elevate=25.0)
plot_u.write_png(save_dir+"/displacement_%.4f"%t)
# save the parameters to file
File(save_dir+"/parameters.xml") << parameters
# get timings and save to file
if MPI.process_number() == 0:
timings_str = timings().str("Timings")
text_file = open(save_dir+"/timings.txt", "w")
text_file.write(timings_str)