Write 2nd order C++ -- hyperasticty example

dolfinx
1

I’m trying to print out a second order mesh using VTK in C++. I’ve been using the hyperelasticty demo as example. However, I keep getting an interpolation error when I interpolate the stress expresson onto the function I use for the output.

 Function element interpolation points has different shape to Expression interpolation points

these are my modification snippets

UFL

e = element("Lagrange", "hexahedron", 2, shape=(3,))
m_e = element("Lagrange", "hexahedron", 1, shape=(3,))

mesh = Mesh(m_e)
V = FunctionSpace(mesh, e)
...

sigma = (1 / J) * diff(psi, F) * F.T

main.cpp

    // Compute Cauchy stress. Construct appropriate Basix element for
    // stress.
    constexpr auto family = basix::element::family::P;
    auto cell_type
        = mesh::cell_type_to_basix_type(mesh->topology()->cell_type());
    constexpr int k = 2;
    constexpr bool discontinuous = true;
    basix::FiniteElement S_element = basix::create_element<U>(
        family, cell_type, k, basix::element::lagrange_variant::unset,
        basix::element::dpc_variant::unset, discontinuous);
    auto S = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace(
        mesh, S_element, std::vector<std::size_t>{3, 3}));
    auto sigma_expression = fem::create_expression<T, U>(
        *expression_hyperelasticity_sigma, {{"u", u}}, {});

    fem::Function sigma = fem::Function<T>(S);
    sigma.name = "cauchy_stress";
    sigma.interpolate(sigma_expression, *mesh);

    // Save solution in VTK format
    io::VTKFile file_s(mesh->comm(), "s.pvd", "w");
    std::vector<std::reference_wrapper<const fem::Function<T>>> sigma_vector = {sigma};
    file_s.write(sigma_vector, 0.0);

I’m using version 0.8.0. Can anybody help?

2

You would have to provide the whole ufl file and preferably a complete cpp file using a unit square (or unit cube), as it is rather hard to do a visual debug of the C++ interface, even for experienced developers.

3

sure here the code.

UFL

from basix.ufl import element
from ufl import (
    Coefficient,
    Constant,
    FunctionSpace,
    Identity,
    Mesh,
    TestFunction,
    TrialFunction,
    derivative,
    det,
    diff,
    ds,
    dx,
    grad,
    inner,
    ln,
    tr,
    variable,
)

# Function spaces
m_e = element("Lagrange", "tetrahedron", 1, shape=(3,))
e = element("Lagrange", "tetrahedron", 2, shape=(3,))
mesh = Mesh(m_e)
V = FunctionSpace(mesh, e)

# Trial and test functions
du = TrialFunction(V)  # Incremental displacement
v = TestFunction(V)  # Test function


# Functions
u = Coefficient(V)  # Displacement from previous iteration
B = Constant(mesh, shape=(3,))  # Body force per unit volume
T = Constant(mesh, shape=(3,))  # Traction force on the boundary

# Kinematics
d = len(u)
I = Identity(d)  # Identity tensor  # noqa: E741
F = variable(I + grad(u))  # Deformation gradient
C = F.T * F  # Right Cauchy-Green tensor

# Invariants of deformation tensors
Ic = tr(C)
J = det(F)

# Elasticity parameters
E = 10.0
nu = 0.3
mu = E / (2 * (1 + nu))
lmbda = E * nu / ((1 + nu) * (1 - 2 * nu))


# Stored strain energy density (compressible neo-Hookean model)
psi = (mu / 2) * (Ic - 3) - mu * ln(J) + (lmbda / 2) * (ln(J)) ** 2

# Total potential energy
Pi = psi * dx - inner(B, u) * dx - inner(T, u) * ds

# First variation of Pi (directional derivative about u in the direction
# of v)
F_form = derivative(Pi, u, v)

# Compute Jacobian of F
J_form = derivative(F_form, u, du)

# Compute Cauchy stress
sigma = (1 / J) * diff(psi, F) * F.T

forms = [F_form, J_form]
elements = [e]
expressions = [(sigma, [[0.25, 0.25, 0.25]])]

Main

// # Hyperelasticity
//
// Solve a compressible neo-Hookean model in 3D.


#include "hyperelasticity.h"
#include <basix/finite-element.h>
#include <climits>
#include <cmath>
#include <dolfinx.h>
#include <dolfinx/common/log.h>
#include <dolfinx/fem/assembler.h>
#include <dolfinx/fem/petsc.h>
#include <dolfinx/io/XDMFFile.h>
#include <dolfinx/la/Vector.h>
#include <dolfinx/mesh/Mesh.h>
#include <dolfinx/mesh/cell_types.h>
#include <dolfinx/nls/NewtonSolver.h>

using namespace dolfinx;
using T = PetscScalar;
using U = typename dolfinx::scalar_value_type_t<T>;

/// Hyperelastic problem class
class HyperElasticProblem
{
public:
  /// Constructor
  HyperElasticProblem(
      std::shared_ptr<fem::Form<T>> L, std::shared_ptr<fem::Form<T>> J,
      std::vector<std::shared_ptr<const fem::DirichletBC<T>>> bcs)
      : _l(L), _j(J), _bcs(bcs),
        _b(L->function_spaces()[0]->dofmap()->index_map,
           L->function_spaces()[0]->dofmap()->index_map_bs()),
        _matA(la::petsc::Matrix(fem::petsc::create_matrix(*J, "aij"), false))
  {
    auto map = L->function_spaces()[0]->dofmap()->index_map;
    const int bs = L->function_spaces()[0]->dofmap()->index_map_bs();
    std::int32_t size_local = bs * map->size_local();

    std::vector<PetscInt> ghosts(map->ghosts().begin(), map->ghosts().end());
    std::int64_t size_global = bs * map->size_global();
    VecCreateGhostBlockWithArray(map->comm(), bs, size_local, size_global,
                                 ghosts.size(), ghosts.data(),
                                 _b.array().data(), &_b_petsc);
  }

  /// Destructor
  virtual ~HyperElasticProblem()
  {
    if (_b_petsc)
      VecDestroy(&_b_petsc);
  }

  /// @brief  Form
  /// @return
  auto form()
  {
    return [](Vec x)
    {
      VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
      VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
    };
  }

  /// Compute F at current point x
  auto F()
  {
    return [&](const Vec x, Vec)
    {
      // Assemble b and update ghosts
      std::span<T> b(_b.mutable_array());
      std::fill(b.begin(), b.end(), 0.0);
      fem::assemble_vector<T>(b, *_l);
      VecGhostUpdateBegin(_b_petsc, ADD_VALUES, SCATTER_REVERSE);
      VecGhostUpdateEnd(_b_petsc, ADD_VALUES, SCATTER_REVERSE);

      // Set bcs
      Vec x_local;
      VecGhostGetLocalForm(x, &x_local);
      PetscInt n = 0;
      VecGetSize(x_local, &n);
      const T* array = nullptr;
      VecGetArrayRead(x_local, &array);
      fem::set_bc<T>(b, _bcs, std::span<const T>(array, n), -1.0);
      VecRestoreArrayRead(x, &array);
    };
  }

  /// Compute J = F' at current point x
  auto J()
  {
    return [&](const Vec, Mat A)
    {
      MatZeroEntries(A);
      fem::assemble_matrix(la::petsc::Matrix::set_block_fn(A, ADD_VALUES), *_j,
                           _bcs);
      MatAssemblyBegin(A, MAT_FLUSH_ASSEMBLY);
      MatAssemblyEnd(A, MAT_FLUSH_ASSEMBLY);
      fem::set_diagonal(la::petsc::Matrix::set_fn(A, INSERT_VALUES),
                        *_j->function_spaces()[0], _bcs);
      MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
    };
  }

  /// RHS vector
  Vec vector() { return _b_petsc; }

  /// Jacobian matrix
  Mat matrix() { return _matA.mat(); }

private:
  std::shared_ptr<fem::Form<T>> _l, _j;
  std::vector<std::shared_ptr<const fem::DirichletBC<T>>> _bcs;
  la::Vector<T> _b;
  Vec _b_petsc = nullptr;
  la::petsc::Matrix _matA;
};

int main(int argc, char* argv[])
{
  init_logging(argc, argv);
  PetscInitialize(&argc, &argv, nullptr, nullptr);

  // Set the logging thread name to show the process rank
  int mpi_rank;
  MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);
  std::string thread_name = "RANK " + std::to_string(mpi_rank);
  loguru::set_thread_name(thread_name.c_str());
  {
    // Create mesh and define function space
    auto mesh = std::make_shared<mesh::Mesh<U>>(mesh::create_box<U>(
        MPI_COMM_WORLD, {{{0.0, 0.0, 0.0}, {1.0, 1.0, 1.0}}}, {1, 1, 1},
        mesh::CellType::tetrahedron,
        mesh::create_cell_partitioner(mesh::GhostMode::none)));

    auto V = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace(
        functionspace_form_hyperelasticity_F_form, "u", mesh));

    auto B = std::make_shared<fem::Constant<T>>(std::vector<T>{0, 0, 0});
    auto traction = std::make_shared<fem::Constant<T>>(std::vector<T>{0, 0, 0});


    // Define solution function
    auto u = std::make_shared<fem::Function<T>>(V);
    auto a = std::make_shared<fem::Form<T>>(
        fem::create_form<T>(*form_hyperelasticity_J_form, {V, V}, {{"u", u}},
                            {{"B", B}, {"T", traction}}, {}));
    auto L = std::make_shared<fem::Form<T>>(
        fem::create_form<T>(*form_hyperelasticity_F_form, {V}, {{"u", u}},
                            {{"B", B}, {"T", traction}}, {}));

    auto u_rotation = std::make_shared<fem::Function<T>>(V);
    u_rotation->interpolate(
        [](auto x) -> std::pair<std::vector<T>, std::vector<std::size_t>>
        {
          constexpr U scale = 0.005;

          // Center of rotation
          constexpr U x1_c = 0.5;
          constexpr U x2_c = 0.5;

          // Large angle of rotation (60 degrees)
          constexpr U theta = 1.04719755;

          // New coordinates
          std::vector<U> fdata(3 * x.extent(1), 0.0);
          MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
              U, MDSPAN_IMPL_STANDARD_NAMESPACE::extents<
                     std::size_t, 3,
                     MDSPAN_IMPL_STANDARD_NAMESPACE::dynamic_extent>>
              f(fdata.data(), 3, x.extent(1));
          for (std::size_t p = 0; p < x.extent(1); ++p)
          {
            U x1 = x(1, p);
            U x2 = x(2, p);
            f(1, p) = scale
                      * (x1_c + (x1 - x1_c) * std::cos(theta)
                         - (x2 - x2_c) * std::sin(theta) - x1);
            f(2, p) = scale
                      * (x2_c + (x1 - x1_c) * std::sin(theta)
                         - (x2 - x2_c) * std::cos(theta) - x2);
          }

          return {std::move(fdata), {3, x.extent(1)}};
        });

    // Create Dirichlet boundary conditions
    auto bdofs_left = fem::locate_dofs_geometrical(
        *V,
        [](auto x)
        {
          constexpr U eps = 1.0e-6;
          std::vector<std::int8_t> marker(x.extent(1), false);
          for (std::size_t p = 0; p < x.extent(1); ++p)
          {
            if (std::abs(x(0, p)) < eps)
              marker[p] = true;
          }
          return marker;
        });
    auto bdofs_right = fem::locate_dofs_geometrical(
        *V,
        [](auto x)
        {
          constexpr U eps = 1.0e-6;
          std::vector<std::int8_t> marker(x.extent(1), false);
          for (std::size_t p = 0; p < x.extent(1); ++p)
          {
            if (std::abs(x(0, p) - 1) < eps)
              marker[p] = true;
          }
          return marker;
        });
    std::vector bcs = {
        std::make_shared<const fem::DirichletBC<T>>(std::vector<T>{0, 0, 0},
                                                    bdofs_left, V),
        std::make_shared<const fem::DirichletBC<T>>(u_rotation, bdofs_right)};

    HyperElasticProblem problem(L, a, bcs);
    nls::petsc::NewtonSolver newton_solver(mesh->comm());
    newton_solver.setF(problem.F(), problem.vector());
    newton_solver.setJ(problem.J(), problem.matrix());
    newton_solver.set_form(problem.form());
    newton_solver.rtol = 10 * std::numeric_limits<T>::epsilon();
    newton_solver.atol = 10 * std::numeric_limits<T>::epsilon();

    la::petsc::Vector _u(la::petsc::create_vector_wrap(*u->x()), false);
    auto [niter, success] = newton_solver.solve(_u.vec());
    std::cout << "Number of Newton iterations: " << niter << std::endl;

    // Compute Cauchy stress. Construct appropriate Basix element for
    // stress.
    constexpr auto family = basix::element::family::P;
    auto cell_type
        = mesh::cell_type_to_basix_type(mesh->topology()->cell_type());
    constexpr int k = 2;
    constexpr bool discontinuous = true;
    basix::FiniteElement S_element = basix::create_element<U>(
        family, cell_type, k, basix::element::lagrange_variant::unset,
        basix::element::dpc_variant::unset, discontinuous);
    auto S = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace(
        mesh, S_element, std::vector<std::size_t>{3, 3}));
    auto sigma_expression = fem::create_expression<T, U>(
        *expression_hyperelasticity_sigma, {{"u", u}}, {});

    auto sigma = fem::Function<T>(S);
    sigma.name = "cauchy_stress";
    sigma.interpolate(sigma_expression, *mesh);

    // Save solution in VTK format
    io::VTKFile file_u(mesh->comm(), "u.pvd", "w");
    file_u.write<T>({*u}, 0.0);

    // Save Cauchy stress in XDMF format
    io::XDMFFile file_sigma(mesh->comm(), "sigma.xdmf", "w");
    file_sigma.write_mesh(*mesh);
    file_sigma.write_function(sigma, 0.0);

    io::VTKFile file_s(mesh->comm(), "s.pvd", "w");
    file_u.write<T>({sigma}, 0.0);

  }

  PetscFinalize();

  return 0;
}
4

The problem is this.
YOu should let the expression be evaluated at the interpolation points of S_element, i.e.

interpolation_element = element("Lagrange", "tetrahedron", 2, shape=(3,3), discontinuous=True)
interpolation_points = []
for ip_i in interpolation_element._x:
    for ip_ij in ip_i:    
        interpolation_points.append(ip_ij)
import numpy as np
expressions = [(sigma, np.vstack(interpolation_points))]

@mscroggs @garth maybe we should make this interface a bit easier to use for the C++ users?

1 Like
5

thanks @dokken, it works now.