I am solving a very simple linear elasticity problem, and would like to output the stresses as well as the displacement. To this end, I am using the project function to project the stress quantity to a TensorFunctionSpace, as advised in previous answers on the forum:
from __future__ import print_function
from dolfin import *
import time
import numpy
# Form compiler options
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["optimize"] = True
# Sub domain for clamp at left end
def left(x, on_boundary):
return x[0] < 0.001 #and on_boundary
# Sub domain for rotation at right end
def right(x, on_boundary):
return x[0] > 0.99 #and on_boundary
# Load mesh and define function space
mesh = UnitCubeMesh(15, 15, 15)
# Define function space
V = VectorFunctionSpace(mesh, "CG", 1)
# Test and trial functions
u = TrialFunction(V)
delta_u = TestFunction(V)
E = 10.0
nu = 0.0
mu = E / (2.0*(1.0 + nu))
lmbda = E*nu / ((1.0 + nu)*(1.0 - 2.0*nu))
# External forces (body and applied tractions
f = Constant((1.0, 0.0, 0.0))
# Set up boundary condition at left end
zero = Constant((0.0, 0.0, 0.0))
bc = DirichletBC(V, zero, left)
def eps(u):
return sym(grad(u))
# Stress tensor
def sigma(u):
return 2.0*mu*eps(u) + lmbda*tr(eps(u))*Identity(len(u))
# Forms
a = inner(sigma(u), eps(delta_u))*dx
L = inner(f, delta_u)*dx
# Time-stepping
u = Function(V)
xdmf_file = XDMFFile(mesh.mpi_comm(), "out.xdmf")
start = time.time()
print("Started solution")
solve(a == L, u,
bcs=bc,
solver_parameters={"linear_solver": "cg"},)
print("Finished solution in %f seconds"%(time.time()-start))
u.rename("u", "displacement")
xdmf_file.write(u, 0.0)
# Start projecting stress
start = time.time()
print("Started stress calculation")
sigma_fs = TensorFunctionSpace(mesh, "CG", 1)
# Following line is very inefficient
sigma_out = project(sigma(u), sigma_fs)
sigma_out.rename("sigma", "stress")
xdmf_file.write(sigma_out, 0.0)
print("Finished stress calculation in %f seconds"%(time.time()-start))
However, as you can see from the output, it takes a lot of time to do the projection. It appears that the time required for the stress computation increases exponentially with the number of elements:
Started solution
Solving linear variational problem.
Finished solution in 0.189142 seconds
Started stress calculation
Finished stress calculation in 17.490856 seconds
I understand that the solver is a black box, but it is often very necessary to me to have not just stress in the output, but also other fields of equal complexity. If I am to see data other than the solution vector in the output, I require a faster way of calculating---or storing---these quantities. I doubt that I am the only one who had such a problem, so I am guessing that you might have a better advice.
Thanks in advance.
