To make the problem more concrete, consider the following trivial 2-d elastic problem defined on a UnitSquareMesh composed of two triangles where the diagonal side appears as "/". All nodes are fixed, except at the right bottom corner where U=[1.0, 0.0] is prescribed. The yield stress is set to an extremely high value. The total elastic energy should be M=0.375. The upper left element is stress free due to BC, and in the lower right one we have s=[1, 0, -0.5].
However using the FEniCS Solid Mechanics interface, I found that M=0. Following the UFL file
degree = 1
elementA = VectorElement("Lagrange", triangle, degree)
elementT = VectorElement("Quadrature", triangle, degree, 9)
elementS = VectorElement("Quadrature", triangle, degree, 3)
u = Coefficient(elementA)
v = TestFunction(elementA)
du = TrialFunction(elementA)
f = Coefficient(elementA)
t = Coefficient(elementT)
s = Coefficient(elementS)
E = Constant(triangle)
nu = Constant(triangle)
def eps(u):
return as_vector([u[i].dx(i) for i in range(2)] + [u[i].dx(j) + u[j].dx(i) for i, j in [(0, 1)]])
def sigma(s):
return as_matrix([[s[0], s[2]], [s[2], s[1]]])
def eps_e(s):
return as_vector([(1-nu**2)/E*s[0]-nu/E*(1+nu)*s[1],
(1-nu**2)/E*s[1]-nu/E*(1+nu)*s[0],
2*(1+nu)*s[2]/E])
def tangent(t):
return as_matrix([[t[i*3 + j] for j in range(3)] for i in range(3)])
a = inner(eps(v), dot(tangent(t), eps(du)))*dx
L = inner(s, eps(v))*dx - inner(f, v)*dx
M = 0.5*inner(s, eps_e(s))*dx
forms = [a, L, M]
and the C++ code
#include <dolfin.h>
#include <FenicsSolidMechanics.h>
#include "Plas2D.h"
using namespace dolfin;
// Displacement at the right bottom corner
class DirichletCondition : public Expression
{
public:
DirichletCondition() : Expression(2) {}
void eval(Array<double>& values, const Array<double>& x) const
{
values[0] = 1.0;
values[1] = 0.0;
}
};
// Subdomain for Dirichlet boundary condition
class LeftTop : public SubDomain
{
bool inside(const Array<double>& x, bool on_boundary) const
{ return x[0] < 0.5 || x[1] > 0.5 ; }
};
// Subdomain for Dirichlet boundary condition
class RightBottom : public SubDomain
{
bool inside(const Array<double>& x, bool on_boundary) const
{ return x[0] > 0.5 && x[1] < 0.5 ; }
};
int main()
{
// Create mesh
UnitSquareMesh mesh(1, 1);
// Young's modulus and Poisson's ratio
const double E = 1.0;
const double nu = 0.0;
const double hardening_parameter = 0.0;
const double yield_stress = 100;
Constant E_(E);
Constant nu_(nu);
// Time parameter
double t = 0.0;
double dt = 0.0;
unsigned int steps = 0;
// Source term, RHS
Constant f(0.0, 0.0);
// Function spaces
Plas2D::FunctionSpace V(mesh);
// Extract elements for stress and tangent
std::shared_ptr<const FiniteElement> element_t;
std::shared_ptr<const FiniteElement> element_s;
Plas2D::CoefficientSpace_t Vt(mesh);
element_t = Vt.element();
Plas2D::CoefficientSpace_s Vs(mesh);
element_s = Vs.element();
// Create boundary condition
DirichletCondition U;
Constant zero(0.0, 0.0);
LeftTop lefttop;
RightBottom rightbottom;
DirichletBC bc1(V, zero, lefttop, "pointwise");
DirichletBC bc2(V, U, rightbottom, "pointwise");
std::vector<const DirichletBC*> bcs;
bcs.push_back(&bc1);
bcs.push_back(&bc2);
// Solution function
Function u(V);
// Object of class von Mises
const fenicssolid::VonMises J2(E, nu, yield_stress, hardening_parameter);
// Constituive update
std::shared_ptr<fenicssolid::ConstitutiveUpdate>
constitutive_update(new fenicssolid::ConstitutiveUpdate(u, *element_s, *Vs.dofmap(), J2));
// Create forms and attach functions
fenicssolid::QuadratureFunction tangent(mesh, element_t,
constitutive_update,
constitutive_update->w_tangent());
Plas2D::Form_a a(V, V);
a.t = tangent;
Plas2D::Form_L L(V);
L.f = f;
fenicssolid::QuadratureFunction stress(mesh, element_s,
constitutive_update->w_stress());
L.s = stress;
// Create PlasticityProblem
fenicssolid::PlasticityProblem nonlinear_problem(a, L, u, tangent,
stress, bcs, J2);
// Create nonlinear solver and set parameters
dolfin::NewtonSolver nonlinear_solver;
nonlinear_solver.parameters["convergence_criterion"] = "residual";
nonlinear_solver.parameters["maximum_iterations"] = 50;
nonlinear_solver.parameters["relative_tolerance"] = 1e-6;
nonlinear_solver.parameters["absolute_tolerance"] = 1e-15;
nonlinear_solver.parameters["linear_solver"] = "superlu_dist";
// nonlinear_solver.parameters["preconditioner"] = "hypre_amg";
// File names for output
File file1("output/disp.xdmf");
// Functional
Plas2D::Form_M M(mesh);
M.s = stress; M.E = E_; M.nu = nu_;
double M_value = 0.0;
// Load-disp info
unsigned int step = 0;
while (step <= steps)
{
std::cout << "step begin: " << step << std::endl;
std::cout << "time: " << t << std::endl;
// Solve non-linear problem
nonlinear_solver.solve(nonlinear_problem, *u.vector());
// Update variables
constitutive_update->update_history();
// Write output to files
file1 << u;
// Gtheta
M_value = assemble(M);
info("M = %.3e\n", M_value);
t += dt;
step++;
}
return 0;
}
Note that the displacement field found with this C++ code gives the correct result. However the total elastic energy functional is incorrectly evaluated. Through an output of the stress in the ConstitutiveUpdate.cpp
_w_stress[num_ip_per_cell*0 + ip] = trial_stress(0);
_w_stress[num_ip_per_cell*1 + ip] = trial_stress(1);
_w_stress[num_ip_per_cell*2 + ip] = trial_stress(3);
std::cout << "stress(0): " << trial_stress(0) << std::endl;
std::cout << "stress(1): " << trial_stress(1) << std::endl;
std::cout << "stress(3): " << trial_stress(3) << std::endl;
I found
step begin: 0
time: 0
stress(0): 0
stress(1): 0
stress(3): 0
stress(0): 0
stress(1): 0
stress(3): 0
Newton iteration 0: r (abs) = 1.000e+00 (tol = 1.000e-15) r (rel) = 1.000e+00 (tol = 1.000e-06)
stress(0): 1
stress(1): 0
stress(3): -0.5
stress(0): 0
stress(1): 0
stress(3): 0
Newton iteration 1: r (abs) = 0.000e+00 (tol = 1.000e-15) r (rel) = 0.000e+00 (tol = 1.000e-06)
Newton solver finished in 1 iterations and 1 linear solver iterations.
M = 0.000e+00
The stress is well calculated in both elements. However the non-zero stress values are overwritten by the zero ones, producing thus a zero elastic energy.
