diff --git a/CMakeLists.txt b/CMakeLists.txt
index 079e6716474f67a882806d0fbc8c8afc986340a5..3328aa26f30b02a53386b076a96108ec398712de 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -27,6 +27,7 @@ add_subdirectory(cppptest)
foreach(AE108_LIBRARY elements cpppetsc assembly solve cmdline)
add_subdirectory(${AE108_LIBRARY})
endforeach()
+add_subdirectory(examples)
include(CMakePackageConfigHelpers)
include(GNUInstallDirs)
diff --git a/examples/Basic.cc b/examples/Basic.cc
new file mode 100644
index 0000000000000000000000000000000000000000..01399e669d82725307b64d2a76282e272ffebe47
--- /dev/null
+++ b/examples/Basic.cc
@@ -0,0 +1,214 @@
+// © 2020 ETH Zurich, Mechanics and Materials Lab
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include "ae108/assembly/Assembler.h"
+#include "ae108/cpppetsc/Mesh.h"
+#include "ae108/cpppetsc/ParallelComputePolicy.h"
+#include "ae108/cpppetsc/Vector.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>
+
+using namespace ae108;
+
+using Policy = cpppetsc::ParallelComputePolicy;
+using Mesh = cpppetsc::Mesh<Policy>;
+using Vector = cpppetsc::Vector<Policy>;
+using BoundaryCondition = cpppetsc::MeshBoundaryCondition<Mesh>;
+
+namespace embedding = ae108::elements::embedding;
+namespace integrator = ae108::elements::integrator;
+namespace materialmodels = ae108::elements::materialmodels;
+namespace quadrature = ae108::elements::quadrature;
+namespace shape = ae108::elements::shape;
+
+// In this example we will simulate two elements, A and B.
+// Each of these elements has four vertices.
+
+// 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 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::value_type, dimension>, number_of_vertices>;
+constexpr auto vertex_positions = VertexPositions{{
+ {{0., 0.}},
+ {{1., 0.}},
+ {{1., 1.}},
+ {{0., 1.}},
+ {{2., 0.}},
+ {{2., 1.}},
+}};
+
+// Let's specify how the elements behave when deformed.
+
+// First, we select a Hookean (linear elastic) material model in 2D.
+
+using MaterialModel = materialmodels::Hookean<dimension>;
+
+// 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,
+ 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>;
+
+// Now we are ready to select our element type.
+
+using Element = elements::CoreElement<MaterialModel, Integrator>;
+
+// We will assemble e.g. energy using a collection of elements. This is done by
+// the assembler. (The list DefaultFeaturePlugins contain the features (e.g.
+// energy) that most elements support.)
+
+using Assembler =
+ assembly::Assembler<Element, assembly::DefaultFeaturePlugins, Policy>;
+
+// Our goals is to minimize the energy. This is done by the solver.
+
+using Solver = solve::NonlinearSolver<Assembler>;
+
+int main(int argc, char **argv) {
+ // PETSc must be initialized before using it.
+ Policy::handleError(PetscInitialize(&argc, &argv, NULL, NULL));
+
+ // We use a scope around our computation to make sure everything is cleaned up
+ // before we call PetscFinalize.
+ {
+ const auto mesh = Mesh::fromConnectivity(
+ dimension, connectivity, number_of_vertices, dof_per_vertex);
+ auto assembler = Assembler();
+
+ const auto model = MaterialModel(1.0, 0.);
+
+ // Depending on whether we use MPI, our mesh may be distributed and not all
+ // elements are present on this computational node.
+
+ // Let's add those elements that are "local" to the assembler.
+
+ for (const auto &element : mesh.localElements()) {
+ assembler.emplaceElement(
+ element, model,
+ Integrator(Embedding(Embedding::Collection<Embedding::PhysicalPoint>{{
+ vertex_positions.at(connectivity.at(element.index()).at(0)),
+ vertex_positions.at(connectivity.at(element.index()).at(1)),
+ vertex_positions.at(connectivity.at(element.index()).at(2)),
+ vertex_positions.at(connectivity.at(element.index()).at(3)),
+ }})));
+ }
+
+ // We need to create a solver. We do not use the time, so we can set it to
+ // zero.
+
+ const auto solver = Solver(&mesh);
+ const auto time = Element::Time{0.};
+
+ // Before we can produce meaningful results, we need to specify boundary
+ // conditions. Let's fix the nodes at x=0 and pull on the nodes at x=2.
+
+ std::vector<BoundaryCondition> boundary_conditions;
+ for (const auto &vertex : mesh.localVertices()) {
+ switch (vertex.index()) {
+ case 0:
+ case 3: {
+ // The displacement in x direction is zero.
+ boundary_conditions.push_back({vertex, 0, 0.});
+ // The displacement in y direction is zero.
+ boundary_conditions.push_back({vertex, 1, 0.});
+ break;
+ }
+ case 4:
+ case 5: {
+ // The displacement in x direction is .5.
+ boundary_conditions.push_back({vertex, 0, .5});
+ // The displacement in y direction is 0.
+ boundary_conditions.push_back({vertex, 1, 0.});
+ break;
+ }
+ }
+ }
+
+ // We are ready to minimize the energy.
+
+ const auto result = solver.computeSolution(
+ boundary_conditions, Vector::fromGlobalMesh(mesh), time, &assembler);
+
+ // As a final step, we print the result.
+ // We could use: result.unwrap().print();
+
+ // However, the output of this command changes when running the application
+ // in MPI-parallel mode. This happens because the vector "result" is
+ // distributed between the ranks in this case.
+
+ // So we gather all the results locally in the canonical data format first:
+ // [ vertex-0-dof-0, vertex-0-dof-1, vertex-1-dof-0, vertex-1-dof-1, ...].
+
+ const auto global_result =
+ Vector::fromDistributedInCanonicalOrder(result, mesh, dof_per_vertex);
+
+ // Then we print this global vector only on the primary rank.
+ // We expect the degrees of freedom in x direction to be between 0 and .5,
+ // and the degrees of freedom in y direction to be 0:
+ // Vertex 0: (0.0, 0.0)
+ // Vertex 1: (.25, 0.0)
+ // Vertex 2: (.25, 0.0)
+ // Vertex 3: (0.0, 0.0)
+ // Vertex 4: (0.5, 0.0)
+ // Vertex 5: (0.5, 0.0)
+
+ if (Policy::isPrimaryRank())
+ global_result.unwrap().print();
+ }
+ Policy::handleError(PetscFinalize());
+}
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..94d6106504373ba16e7a1798c5a1c618a4b7be12
--- /dev/null
+++ b/examples/CMakeLists.txt
@@ -0,0 +1,40 @@
+# © 2020 ETH Zurich, Mechanics and Materials Lab
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+cmake_minimum_required(VERSION 3.12 FATAL_ERROR)
+
+project(ae108-examples LANGUAGES C CXX)
+
+if (NOT TARGET ae108::elements)
+ find_package(ae108 REQUIRED)
+endif()
+
+add_executable(${PROJECT_NAME}-Basic Basic.cc)
+target_link_libraries(${PROJECT_NAME}-Basic
+ PRIVATE ae108::solve
+ PRIVATE ae108::assembly
+ PRIVATE ae108::elements
+)
+
+add_executable(${PROJECT_NAME}-MeshGeneration MeshGeneration.cc)
+target_link_libraries(${PROJECT_NAME}-MeshGeneration
+ PRIVATE ae108::solve
+ PRIVATE ae108::assembly
+ PRIVATE ae108::elements
+)
+
+add_executable(${PROJECT_NAME}-Output Output.cc)
+target_link_libraries(${PROJECT_NAME}-Output
+ PRIVATE ae108::cpppetsc
+)
diff --git a/examples/MeshGeneration.cc b/examples/MeshGeneration.cc
new file mode 100644
index 0000000000000000000000000000000000000000..0280241fad494b803c1811aad161f2fabf4d6c5e
--- /dev/null
+++ b/examples/MeshGeneration.cc
@@ -0,0 +1,101 @@
+// © 2020 ETH Zurich, Mechanics and Materials Lab
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include "ae108/cpppetsc/Mesh.h"
+#include "ae108/cpppetsc/ParallelComputePolicy.h"
+#include "ae108/elements/mesh/Mesh.h"
+#include "ae108/elements/mesh/generate_triangle_mesh.h"
+#include "ae108/elements/tensor/Tensor.h"
+#include <array>
+
+using namespace ae108;
+
+using Policy = cpppetsc::ParallelComputePolicy;
+using Mesh = cpppetsc::Mesh<Policy>;
+
+// In this example we will mesh a rectangle with triangles and distribute the
+// result with PETSc.
+
+// Let's specify the parameters of this mesh.
+
+constexpr auto dimension = Mesh::size_type{2};
+constexpr auto dof_per_vertex = Mesh::size_type{2};
+
+// We choose the rectangle to have x-length=3 and y-length=2.
+
+constexpr auto rectangle_size = elements::tensor::Tensor<double, 2>{{3., 2.}};
+
+// Also, we want to generate 2 * 4 triangles in x-direction, and 2 * 1 triangles
+// in y-direction.
+
+constexpr auto mesh_granularity =
+ elements::tensor::Tensor<std::size_t, 2>{{4, 1}};
+
+int main(int argc, char **argv) {
+ // PETSc must be initialized before using it.
+ Policy::handleError(PetscInitialize(&argc, &argv, NULL, NULL));
+
+ // We use a scope around our computation to make sure everything is cleaned up
+ // before we call PetscFinalize.
+ {
+ // First we generate the mesh locally. This mesh will look the same on every
+ // rank.
+
+ const auto generated_mesh = elements::mesh::generate_triangle_mesh(
+ rectangle_size, mesh_granularity);
+
+ // The number of positions in the mesh provides the number of vertices here.
+
+ const auto number_of_vertices = generated_mesh.number_of_positions();
+
+ // Let's create the distributed PETSc-mesh using the generated connectivity
+ // and print it to the console. We expect 2 * 4 2-cells (triangles) and
+ // 2 * (4 + 1) 0-cells (vertices). In parallel, some vertices are
+ // shared between ranks and we therefore expect a sum of at least 10
+ // 0-cells.
+
+ const auto mesh =
+ Mesh::fromConnectivity(dimension, generated_mesh.connectivity(),
+ number_of_vertices, dof_per_vertex);
+ mesh.print();
+
+ // Before we continue, we wait until all ranks are done printing.
+
+ Policy::handleError(PetscBarrier(nullptr));
+
+ // The PETSc mesh does not contain vertex positions, only the generated mesh
+ // does. Therefore, use the generated mesh to access vertex positions using
+ // an index. To test that this works, we print the position of the second
+ // vertex (index 1) on every rank that shares it. We expect this vertex to
+ // be at (0., 2.).
+
+ for (const auto &vertex : mesh.localVertices()) {
+ if (vertex.index() != 1) {
+ continue;
+ }
+
+ // We have found the vertex with index 1 locally.
+
+ const auto position = generated_mesh.position_of_vertex(1);
+ printf("The second vertex is at (%f, %f).\n", position.at(0),
+ position.at(1));
+ }
+
+ // If you choose enough processes (in my case 4 works), you will see that
+ // the vertex is shared by two ranks because more than one rank prints the
+ // message.
+ }
+
+ Policy::handleError(PetscFinalize());
+}
diff --git a/examples/Output.cc b/examples/Output.cc
new file mode 100644
index 0000000000000000000000000000000000000000..d39f14e4f8375cccab24fce195524117f4f53c39
--- /dev/null
+++ b/examples/Output.cc
@@ -0,0 +1,128 @@
+// © 2020 ETH Zurich, Mechanics and Materials Lab
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include "ae108/cpppetsc/Mesh.h"
+#include "ae108/cpppetsc/ParallelComputePolicy.h"
+#include "ae108/cpppetsc/Vector.h"
+#include "ae108/cpppetsc/Viewer.h"
+#include "ae108/cpppetsc/createVectorFromSource.h"
+#include "ae108/cpppetsc/setName.h"
+#include "ae108/cpppetsc/writeToViewer.h"
+#include <algorithm>
+#include <array>
+
+using namespace ae108;
+
+using Policy = cpppetsc::ParallelComputePolicy;
+using Mesh = cpppetsc::Mesh<Policy>;
+using Vector = cpppetsc::Vector<Policy>;
+using Viewer = cpppetsc::Viewer<Policy>;
+
+// In this example we will save data for two elements, A and B, to a HDF5 file.
+// Each of these elements has four vertices.
+
+// 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 topological_dimension = Mesh::TopologicalDimension{2};
+constexpr auto coordinate_dimension = Mesh::CoordinateDimension{3};
+constexpr auto dof_per_vertex = Mesh::size_type{1};
+constexpr auto dof_per_element = Mesh::size_type{0};
+
+// 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::value_type, coordinate_dimension>,
+ number_of_vertices>;
+constexpr auto vertex_positions = VertexPositions{{
+ {{0., 0., 0.}},
+ {{1., 0., 0.}},
+ {{1., 1., 0.}},
+ {{0., 1., 0.}},
+ {{2., 0., 0.}},
+ {{2., 1., 0.}},
+}};
+
+int main(int argc, char **argv) {
+ // PETSc must be initialized before using it.
+ Policy::handleError(PetscInitialize(&argc, &argv, NULL, NULL));
+
+ // We use a scope around our computation to make sure everything is cleaned up
+ // before we call PetscFinalize.
+ {
+ // First we create a mesh.
+ const auto mesh = Mesh::fromConnectivity(
+ topological_dimension, coordinate_dimension, connectivity,
+ number_of_vertices, dof_per_vertex, dof_per_element);
+
+ // Then we create a vector of coordinates from `vertex_positions`.
+ using DataSource = std::function<void(Mesh::size_type, Mesh::value_type *)>;
+ auto coordinates = cpppetsc::createVectorFromSource(
+ mesh, coordinate_dimension,
+ DataSource(
+ [&](const Mesh::size_type index, Mesh::value_type *const out) {
+ const auto &position = vertex_positions.at(index);
+ std::copy(position.begin(), position.end(), out);
+ }));
+ cpppetsc::setName("coordinates", &coordinates);
+
+ // Let's create a global vector and fill it with the vertex indices.
+ constexpr auto data_per_vertex = coordinate_dimension + 2;
+ auto data = cpppetsc::createVectorFromSource(
+ mesh, data_per_vertex,
+ DataSource(
+ [&](const Mesh::size_type index, Mesh::value_type *const out) {
+ std::fill_n(out, data_per_vertex,
+ static_cast<Mesh::value_type>(index));
+ }));
+ cpppetsc::setName("vertex_index_data", &data);
+
+ // Finally we write the results to "output.ae108".
+ // We start with writing the mesh.
+ auto viewer = Viewer::fromHdf5FilePath("output.ae108");
+ cpppetsc::writeToViewer(mesh, coordinates, &viewer);
+
+ // Then we write two vector fields of different dimension.
+ cpppetsc::writeToViewer(data, &viewer); // dimension 5
+ cpppetsc::writeToViewer(coordinates, &viewer); // dimension 3
+
+ // Note that Paraview supports both "point arrays" and "cell arrays"
+ cpppetsc::writeToViewer(
+ Mesh::vector_type::fromGlobalMesh(mesh.cloneWithDofs(1, 0)), &viewer);
+ cpppetsc::writeToViewer(
+ Mesh::vector_type::fromGlobalMesh(mesh.cloneWithDofs(0, 1)), &viewer);
+
+ fprintf(stderr,
+ "The data has been written to the file \"output.ae108\".\n");
+ }
+ Policy::handleError(PetscFinalize());
+}
\ No newline at end of file