Setting up a minimal end-to-end `ae108` example

An ae108 example wires together five things: a mesh, an element type (shape + material model), an assembler, boundary conditions, and a solver — then writes the result out. This walks through examples/Basic.cc, the simplest complete example: two quad elements sharing an edge, fixed on one end and pulled on the other.

1. Aliases & includes

Each example starts by pulling in the pieces it needs and aliasing them to short, domain-level names, so the rest of the file reads in terms of Mesh, Vector, EssentialBoundaryCondition, etc. instead of fully-qualified ae108::... names:

examples/Basic.cc - Aliases and Includes

#include "ae108/assembly/Assembler.h"
#include "ae108/assembly/DefaultFeaturePlugins.h"
#include "ae108/cpppetsc/BoundaryCondition.h"
#include "ae108/cpppetsc/Context.h"
#include "ae108/cpppetsc/Mesh.h"
#include "ae108/cpppetsc/Vector.h"
#include "ae108/cpppetsc/mpi.h"
#include "ae108/elements/CoreElement.h"
#include "ae108/elements/embedding/IsoparametricEmbedding.h"
#include "ae108/elements/integrator/IsoparametricIntegrator.h"
#include "ae108/elements/materialmodels/Hookean.h"
#include "ae108/elements/quadrature/Quadrature.h"
#include "ae108/elements/shape/Quad4.h"
#include "ae108/solve/NonlinearSolver.h"
#include <array>
#include <vector>

namespace cpppetsc = ae108::cpppetsc;
namespace elements = ae108::elements;
namespace solve = ae108::solve;
namespace assembly = ae108::assembly;

using Context = cpppetsc::Context;
using Mesh = cpppetsc::Mesh;
using Vector = cpppetsc::Vector;
using EssentialBoundaryCondition = cpppetsc::EssentialBoundaryCondition;
template <int Dim> using LabelConstraints = cpppetsc::LabelConstraints<Dim>;

Policy (here ParallelComputePolicy) selects the compute backend and is threaded through nearly every other type; Context must be constructed before anything else touches MPI/PETSc/cpppetsc.

2. Describe the geometry

Geometry is just vertex positions plus, per element, which vertex indices it uses. Basic.cc connects two quads (A and B) sharing the edge between vertices 1 and 2:

examples/Basic.cc - Mesh Setup

//  3 ------- 2 ------- 5
//  |         |         |
//  |    A    |    B    |
//  0---------1---------4

// Let's specify the parameters of this mesh.

constexpr auto number_of_vertices_per_element = Mesh::size_type{4};
constexpr auto number_of_elements = Mesh::size_type{2};
constexpr auto number_of_vertices = Mesh::size_type{6};
constexpr auto vertex_dimension = Mesh::size_type{2};
constexpr auto dof_per_vertex = Mesh::size_type{2};

// The connectivity specifies the vertex indices per Element.

using Connectivity =
    std::array<std::array<Mesh::size_type, number_of_vertices_per_element>,
               number_of_elements>;
constexpr auto connectivity = Connectivity{{
    {{0, 1, 2, 3}}, // vertices of element A
    {{1, 4, 5, 2}}, // vertices of element B
}};

// Each of the vertices is located at the following coordinates in physical
// space.

using VertexPositions =
    std::array<std::array<Mesh::real_type, vertex_dimension>,
               number_of_vertices>;
constexpr auto vertex_positions = VertexPositions{{
    {{0., 0.}},
    {{1., 0.}},
    {{1., 1.}},
    {{0., 1.}},
    {{2., 0.}},
    {{2., 1.}},
}};

3. Choose element behavior (the type-building chain)

An element type is assembled from smaller pieces, each answering one question:

Type	Answers
`MaterialModel` (`Hookean<vertex_dimension, ...>`)	How does the material respond to strain?
`Shape` (`Quad4`)	What are the reference-space shape functions?
`Embedding` (`IsoparametricEmbedding<Shape>`)	How does reference space map to physical space?
`Quadrature`	Which points/weights are used to integrate?
`Integrator`	Combines `Quadrature` + `Embedding` to integrate over physical space.
`Element` (`CoreElement<MaterialModel, Integrator, ...>`)	The element used by the assembler.

examples/Basic.cc - Type Aliases

using MaterialModel =
    materialmodels::Hookean<vertex_dimension, Vector::value_type,
                            Vector::real_type>;

// We'll choose the Quad4 shape functions (4 shape functions) that are defined
// in 2D reference space.

using Shape = shape::Quad4;

// Since also our physical space is 2D we may use the isoparametric embedding.

using Embedding = embedding::IsoparametricEmbedding<Shape>;

// Let's choose a simple quadrature rule in reference space.

using Quadrature = quadrature::Quadrature<quadrature::QuadratureType::Cube,
                                          vertex_dimension, 1 /* order */>;

// We use an "Integrator" to integrate in physical space. Of course, it depends
// on the selected quadrature rule and the embedding. In the case of the
// isoparametric embedding, this embedding is specified by the "Shape".

using Integrator =
    integrator::IsoparametricIntegrator<Shape, Quadrature, Vector::value_type,
                                        Vector::real_type>;

// Now we are ready to select our element type.

using Element = elements::CoreElement<MaterialModel, Integrator,
                                      Vector::value_type, Vector::real_type>;

4. Assembler

Basic.cc has only one element type, so a single Assembler is enough (no AssemblerGroup needed). DefaultFeaturePlugins supplies the default assembly features (energy/force/stiffness) that most elements support:

examples/Basic.cc - Assembler and Solver

using Assembler = assembly::Assembler<Element, assembly::DefaultFeaturePlugins>;

// Our goals is to minimize the energy. This is done by the solver.

using Solver = solve::NonlinearSolver<Assembler>;

5. Build the mesh and populate the assembler

examples/Basic.cc - Mesh Creation

  auto mesh = Mesh::fromConnectivity<cpppetsc::ElementType::Quadrilateral,
                                     vertex_dimension>(
      connectivity, vertex_positions, dof_per_vertex, 0);

examples/Basic.cc - Assembly

  for (const auto &element : mesh.localElements()) {
    auto element_index = element.localIndex();

    const auto cell_coords =
        mesh.closureCoordinates<vertex_dimension>(element_index);

    assembler.emplaceElement(
        element, Element{.model = model,
                         .integrator = Integrator(Embedding(
                             Embedding::Collection<Embedding::PhysicalPoint>{{
                                 cell_coords[0],
                                 cell_coords[1],
                                 cell_coords[2],
                                 cell_coords[3],
                             }}))});
  }

Iterating mesh.localElements() rather than indexing the connectivity directly is what makes this loop MPI-safe: under a distributed mesh, each rank only sees (and only needs to assemble) the elements it locally owns.

6. Boundary conditions

Boundary conditions are set up as their own three-step pipeline — tag mesh points, describe constrained dofs, then resolve to local dof indices — which is covered in full in boundaryconditions.md. Here the left edge (x=0) is fully fixed and the right edge (x=2) is pulled by 0.5 in x:

examples/Basic.cc - Define Boundary Conditions

  // Define label value for all elements that should be constrained
  constexpr auto left_bc_tag = 5;
  constexpr auto right_bc_tag = 6;

  // Associate a tag with a constrained dof at the nodes
  const auto essential_bcs = std::vector<EssentialBoundaryCondition>{
      {.label_value = left_bc_tag, .constrained_dofs{{0, 0.0}, {1, 0.0}}},
      {.label_value = right_bc_tag, .constrained_dofs{{0, 0.5}, {1, 0.0}}}};

  // Define the a condition on coordinates that associates a edges/faces
  //   with
  // label value.
  const auto functional_bc = std::vector<LabelConstraints<vertex_dimension>>{
      {.condition =
           [](std::span<const PetscReal> coords) {
             const auto tolerance = .1;
             const auto x_pos = coords[0];
             return (std::abs(x_pos - 0.) < tolerance);
           },
       .label_value = left_bc_tag},
      {.condition =
           [](std::span<const PetscReal> coords) {
             const auto tolerance = .1;

             const auto x_pos = coords[0];
             return (std::abs(x_pos - 2.) < tolerance);
           },
       .label_value = right_bc_tag},
  };

examples/Basic.cc - Apply Boundary Conditions

  // Assign a label value to the boundary elements
  mesh.markBoundaryCells(functional_bc);

  // Get the actual constrained dofs to be used in the solver
  const auto boundary_conditions =
      extractBoundaryDofs<Mesh, dof_per_vertex>(mesh, essential_bcs);

7. Solve

examples/Basic.cc - Compute Solution

  auto initial_guess = Vector::fromGlobalMesh(mesh);

  // We are ready to minimize the energy.
  const auto result = solver.computeSolution(
      boundary_conditions, std::move(initial_guess), time, &assembler);

8. Output

Basic.cc writes the (undeformed) mesh to VTU up front, then, after solving, both prints the result (reordered into canonical per-vertex-dof order) and writes the deformed result to VTU:

examples/Basic.cc - Write and Print Solution

  result.unwrap().print();

  // Write a vtu file with the result that can directly be visualized in
  // Paraview
  result.unwrap().writeAsVtu(mesh, "ae108-examples-basic.vtu");
  auto full_result_local =
      Vector::fromDistributedInCanonicalOrder(result, mesh);
  full_result_local.unwrap().writeAsAscii(mesh, "ae108-examples-basic.txt");
}

Minimal end-to-end skeleton

int main(int argc, char **argv) {
  const auto context = Context(&argc, &argv);          // 0. init MPI/PETSc

  // 1-2. geometry: Connectivity + VertexPositions      (section 2)
  // 3. element type chain: MaterialModel/Shape/Embedding/Quadrature/
  //    Integrator/Element                               (section 3)
  // 4. Assembler + Solver                                (section 4)

  auto mesh = Mesh::fromConnectivity<cpppetsc::ElementType::Quadrilateral,
                                     vertex_dimension>(
      connectivity, vertex_positions, dof_per_vertex, 0);
  auto assembler = Assembler();

  for (const auto &element : mesh.localElements()) {    // 5. populate assembler
    // emplaceElement(...) built from cell_coords via Integrator/Embedding
  }

  const auto solver = Solver(&mesh);

  mesh.markBoundaryCells(functional_bc);                 // 6. boundary conditions
  const auto boundary_conditions =
      extractBoundaryDofs<Mesh, dof_per_vertex>(mesh, essential_bcs);

  const auto result = solver.computeSolution(            // 7. solve
      boundary_conditions, Vector::fromGlobalMesh(mesh), /* time */ 0., &assembler);

  result.unwrap().writeAsVtu(mesh, "result.vtu");         // 8. output
}

See examples/Basic.cc for the full, runnable version.

Building & running it

The example is registered in examples/CMakeLists.txt as an executable target linking against the libraries used above:

examples/CMakeLists.txt - Target Definition

add_executable(${PROJECT_NAME}-Basic Basic.cc)
target_link_libraries(${PROJECT_NAME}-Basic
                      PRIVATE ae108::solve
                      PRIVATE ae108::assembly
                      PRIVATE ae108::elements
)

Build it as part of the project's normal CMake configure/build, then run the resulting ae108-examples-Basic binary — it prints the solved displacements on the primary rank and writes ae108-examples-basic-mesh.vtu and ae108-examples-basic.vtu, which can be opened in ParaView (or similar) to inspect the mesh and the solved deformation.

Setting up a minimal end-to-end ae108 example