add three basic examples

- mesh generation, assembler and solver, output

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)

examples/Basic.cc

+// © 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());
+}

examples/CMakeLists.txt

+# © 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
+)

examples/MeshGeneration.cc

+// © 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());
+}

examples/Output.cc

+// © 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