diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index cb39260e74c93ce362e3f4f8ff2e182cb1f9362d..5ca9f63479e14909fd9abbae3d620ae5bb16e3cb 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -100,6 +100,7 @@ check-script:
.build-library: &build-library
stage: build
+ coverage: '/^lines:\s*\d+\.\d+\%/'
script:
- useradd developer
- mkdir build
@@ -115,7 +116,9 @@ check-script:
'
artifacts:
reports:
- cobertura: build/coverage.xml
+ coverage_report:
+ coverage_format: cobertura
+ path: build/coverage.xml
junit: build/*/test/gtest-results.xml
when: always
diff --git a/cpppetsc/src/include/ae108/cpppetsc/Mesh.h b/cpppetsc/src/include/ae108/cpppetsc/Mesh.h
index 01e343be68a4891e5c4d581dba1a5501de11b43e..470d21d124387eb4607f752e616c40aae728baf5 100644
--- a/cpppetsc/src/include/ae108/cpppetsc/Mesh.h
+++ b/cpppetsc/src/include/ae108/cpppetsc/Mesh.h
@@ -189,9 +189,15 @@ private:
static UniqueEntity<DM> createDM();
/**
- * @brief Adds the default uniform section to the mesh.
+ * @brief Creates an uniform section with the provided degrees of freedom.
*/
- void addSection(const size_type dofPerVertex, const size_type dofPerElement);
+ UniqueEntity<PetscSection> createSection(const size_type dofPerVertex,
+ const size_type dofPerElement) const;
+
+ /**
+ * @brief Sets the default section of the mesh.
+ */
+ void setSection(UniqueEntity<PetscSection> section);
/**
* @brief Adds the chart as described by the connectivity to the mesh.
@@ -237,6 +243,9 @@ private:
void assertCorrectBaseMesh(const vector_type &vector) const;
UniqueEntity<DM> _mesh;
+ UniqueEntity<PetscSF> _migration;
+ UniqueEntity<PetscSection> _section;
+ UniqueEntity<PetscSF> _globalToNatural;
size_type _totalNumberOfElements = 0;
size_type _totalNumberOfVertices = 0;
};
@@ -309,7 +318,7 @@ Mesh<Policy> Mesh<Policy>::fromConnectivity(const size_type dimension,
distributeMesh(&mesh);
- mesh.addSection(dofPerVertex, dofPerElement);
+ mesh.setSection(mesh.createSection(dofPerVertex, dofPerElement));
mesh.setGlobalToNaturalSF(mesh.createReorderingSF(
mesh.canonicalRowIndices(dofPerVertex, dofPerElement)));
@@ -386,8 +395,9 @@ UniqueEntity<PetscSF> Mesh<Policy>::createReorderingSF(
template <class Policy>
void Mesh<Policy>::setGlobalToNaturalSF(UniqueEntity<PetscSF> globalToNatural) {
+ _globalToNatural = std::move(globalToNatural);
Policy::handleError(
- DMPlexSetGlobalToNaturalSF(_mesh.get(), globalToNatural.release()));
+ DMPlexSetGlobalToNaturalSF(_mesh.get(), _globalToNatural.get()));
}
template <class Policy>
@@ -405,13 +415,17 @@ Mesh<Policy> Mesh<Policy>::cloneWithDofs(const size_type dofPerVertex,
auto migration = PetscSF();
Policy::handleError(DMPlexGetMigrationSF(_mesh.get(), &migration));
if (migration) {
- auto cloned = PetscSF();
+ mesh._migration = [&]() {
+ auto sf = PetscSF();
+ Policy::handleError(
+ PetscSFDuplicate(migration, PETSCSF_DUPLICATE_GRAPH, &sf));
+ return makeUniqueEntity<Policy>(sf);
+ }();
Policy::handleError(
- PetscSFDuplicate(migration, PETSCSF_DUPLICATE_GRAPH, &cloned));
- Policy::handleError(DMPlexSetMigrationSF(mesh._mesh.get(), cloned));
+ DMPlexSetMigrationSF(mesh._mesh.get(), mesh._migration.get()));
}
- mesh.addSection(dofPerVertex, dofPerElement);
+ mesh.setSection(mesh.createSection(dofPerVertex, dofPerElement));
mesh.setGlobalToNaturalSF(mesh.createReorderingSF(
mesh.canonicalRowIndices(dofPerVertex, dofPerElement)));
@@ -479,15 +493,17 @@ void Mesh<Policy>::addChart(const Container &elementVertexIDs,
}
template <class Policy>
-void Mesh<Policy>::addSection(const size_type dofPerVertex,
- const size_type dofPerElement) {
- auto section = PetscSection();
+UniqueEntity<PetscSection>
+Mesh<Policy>::createSection(const size_type dofPerVertex,
+ const size_type dofPerElement) const {
const std::array<size_type, 1> numberOfComponents = {{1}};
std::vector<size_type> numberOfDofsPerDim(2);
numberOfDofsPerDim.front() = dofPerVertex;
numberOfDofsPerDim.back() = dofPerElement;
Policy::handleError(DMSetNumFields(_mesh.get(), numberOfComponents.size()));
+
+ auto section = PetscSection();
Policy::handleError(DMPlexCreateSection(
_mesh.get(), nullptr /* label */, numberOfComponents.data(),
numberOfDofsPerDim.data(), 0 /* number of boundary conditions */,
@@ -495,11 +511,13 @@ void Mesh<Policy>::addSection(const size_type dofPerVertex,
nullptr /* boundary condition components */,
nullptr /* boundary condition points */, nullptr /* permutation */,
§ion));
+ return makeUniqueEntity<Policy>(section);
+}
- Policy::handleError(DMSetSection(_mesh.get(), section));
-
- // create default global section
- Policy::handleError(DMGetGlobalSection(_mesh.get(), §ion));
+template <class Policy>
+void Mesh<Policy>::setSection(UniqueEntity<PetscSection> section) {
+ _section = std::move(section);
+ Policy::handleError(DMSetSection(_mesh.get(), _section.get()));
}
template <class Policy> void Mesh<Policy>::distributeMesh(Mesh *const mesh) {
diff --git a/cppptest/src/include/ae108/cppptest/Matchers.h b/cppptest/src/include/ae108/cppptest/Matchers.h
index f481ae688163a3a9341d80d1b425f0f8fd8b72d4..f38cb00d0b00ae7f8faf07809389bb7beefc2e75 100644
--- a/cppptest/src/include/ae108/cppptest/Matchers.h
+++ b/cppptest/src/include/ae108/cppptest/Matchers.h
@@ -98,6 +98,31 @@ MATCHER_P3(ScalarEqIfLocal, row, col, reference,
almost_equal_to_reference(arg(row, col), reference, result_listener);
}
+/**
+ * @brief Check that the value at index row is almost equal to reference (if it
+ * is available locally).
+ */
+MATCHER_P3(ScalarNearIfLocal, row, reference, tolerance,
+ std::string(negation ? "not " : "") + "equal to " +
+ ::testing::PrintToString(reference) + " at (" +
+ ::testing::PrintToString(row) + ")") {
+ return !isLocal(arg, row) ||
+ near_to_reference(arg(row), reference, tolerance, result_listener);
+}
+
+/**
+ * @brief Check that the value at index (row, col) is almost equal to reference
+ * (if it is available locally).
+ */
+MATCHER_P4(ScalarNearIfLocal, row, col, reference, tolerance,
+ std::string(negation ? "not " : "") + "equal to " +
+ ::testing::PrintToString(reference) + " at (" +
+ ::testing::PrintToString(row) + ", " +
+ ::testing::PrintToString(col) + ")") {
+ return !isLocal(arg, row) || near_to_reference(arg(row, col), reference,
+ tolerance, result_listener);
+}
+
/**
* @brief Check that the value has the given address.
*/
@@ -136,4 +161,4 @@ MATCHER_P2(ScalarNear, reference, tolerance,
}
} // namespace
} // namespace cppptest
-} // namespace ae108
\ No newline at end of file
+} // namespace ae108
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 1500634e4d2e8a6794838b15554f4425fc7bdfbb..ab8a805099b5ed3e06e45d15e992f5b57ef71e10 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -112,4 +112,11 @@ target_link_libraries(${PROJECT_NAME}-DynamicSolver
PRIVATE ae108::solve
PRIVATE ae108::assembly
PRIVATE ae108::elements
+)
+
+add_executable(${PROJECT_NAME}-CantileverBeam CantileverBeam.cc)
+target_link_libraries(${PROJECT_NAME}-CantileverBeam
+ PRIVATE ae108::solve
+ PRIVATE ae108::assembly
+ PRIVATE ae108::elements
)
\ No newline at end of file
diff --git a/examples/CantileverBeam.cc b/examples/CantileverBeam.cc
new file mode 100644
index 0000000000000000000000000000000000000000..5c36975f884752865ce34eaa6e8f0a81d7ee1fb7
--- /dev/null
+++ b/examples/CantileverBeam.cc
@@ -0,0 +1,217 @@
+// © 2021 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/assembly/AssemblerGroup.h"
+#include "ae108/cpppetsc/Context.h"
+#include "ae108/cpppetsc/Mesh.h"
+#include "ae108/cpppetsc/ParallelComputePolicy.h"
+#include "ae108/cpppetsc/Vector.h"
+#include "ae108/elements/ForceElement.h"
+#include "ae108/elements/TimoshenkoBeamElement.h"
+#include "ae108/elements/tensor/as_vector.h"
+#include "ae108/solve/NonlinearSolver.h"
+
+using namespace ae108;
+
+using Policy = cpppetsc::ParallelComputePolicy;
+using Context = cpppetsc::Context<Policy>;
+using Mesh = cpppetsc::Mesh<Policy>;
+using Vector = cpppetsc::Vector<Policy>;
+using BoundaryCondition = cpppetsc::MeshBoundaryCondition<Mesh>;
+
+// In this example we will calculate the tip displacement of a cantilver beam
+// with a force on the free end.
+
+// 0---A---1/B
+
+// Let's specify the parameters of this mesh.
+
+constexpr auto number_of_vertices_per_element = Mesh::size_type{2};
+constexpr auto number_of_elements = Mesh::size_type{2};
+constexpr auto number_of_vertices = Mesh::size_type{2};
+constexpr auto coordinate_dimension = Mesh::size_type{3};
+constexpr auto dof_per_vertex = Mesh::size_type{6};
+constexpr auto dof_per_element = Mesh::size_type{0};
+
+// The connectivity specifies the vertex indices per Element.
+
+using Connectivity =
+ std::array<std::vector<Mesh::size_type>, number_of_elements>;
+const auto connectivity = Connectivity{{
+ {{0, 1}}, // vertices of element A
+ {{1}}, // vertex of tip load B
+}};
+
+// Vertices 1 and 2 are located at the following coordinates in physical
+// space.
+
+using VertexPositions =
+ std::array<std::array<Mesh::real_type, coordinate_dimension>,
+ number_of_vertices>;
+constexpr auto vertex_positions = VertexPositions{{
+ {{0., 0., 0.}},
+ {{1., 0., 0.}},
+}};
+
+// Now we are ready to select the Timoshenko beam element
+
+using Element =
+ elements::TimoshenkoBeamElement<coordinate_dimension, Vector::value_type,
+ Vector::real_type>;
+
+// The Timoshenko beam element comes with a number of geometrical and material
+// related properties. These are stored in the Properties struct
+
+using Properties =
+ elements::TimoshenkoBeamProperties<Mesh::real_type, coordinate_dimension>;
+
+// Let us define some parameters for the linear elastic beam of circular cross
+// section
+
+constexpr Mesh::real_type young_modulus = 1.;
+constexpr Mesh::real_type poisson_ratio = 0.3;
+constexpr Mesh::real_type shear_modulus =
+ young_modulus / (2 * (1 + poisson_ratio));
+constexpr Mesh::real_type shear_correction_factor =
+ (7 + 6 * poisson_ratio) / 6 / (1 + poisson_ratio); // Cowper (1966)
+constexpr Mesh::real_type radius = 0.05;
+constexpr Mesh::real_type area = radius * radius * M_PI;
+constexpr Mesh::real_type area_moment =
+ M_PI_4 * radius * radius * radius * radius;
+constexpr Mesh::real_type polar_moment =
+ M_PI_2 * radius * radius * radius * radius;
+
+// 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>;
+
+using ForceElement = elements::ForceElement<dof_per_vertex, Vector::value_type,
+ Vector::real_type>;
+using ForceAssembler =
+ assembly::Assembler<ForceElement, assembly::DefaultFeaturePlugins, Policy>;
+
+// Finally we combine the assemblers in a AssemblerGroup.
+
+using GroupAssembler = assembly::AssemblerGroup<Assembler, ForceAssembler>;
+
+// Our goals is to minimize the energy. This is done by the solver.
+
+using Solver = solve::NonlinearSolver<GroupAssembler>;
+
+int main(int argc, char **argv) {
+ // MPI/PETSc/cpppetsc must be initialized before using it.
+
+ const auto context = Context(&argc, &argv);
+
+ // Let's create the mesh and an assembler.
+
+ const auto mesh =
+ Mesh::fromConnectivity(coordinate_dimension, connectivity,
+ number_of_vertices, dof_per_vertex, 0);
+ auto assembler = GroupAssembler();
+ auto &element_assembler = assembler.get<0>();
+ auto &force_assembler = assembler.get<1>();
+
+ const Properties properties = {young_modulus,
+ shear_modulus,
+ shear_correction_factor,
+ shear_correction_factor,
+ area,
+ area_moment,
+ area_moment,
+ polar_moment};
+
+ // 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()) {
+ switch (element.index()) {
+ case 0: {
+ elements::tensor::Tensor<Mesh::real_type, coordinate_dimension>
+ element_axis;
+ elements::tensor::as_vector(&element_axis) =
+ elements::tensor::as_vector(
+ &vertex_positions.at(connectivity.at(element.index()).at(1))) -
+ elements::tensor::as_vector(
+ &vertex_positions.at(connectivity.at(element.index()).at(0)));
+
+ element_assembler.emplaceElement(
+ element, timoshenko_beam_stiffness_matrix(element_axis, properties));
+ break;
+ }
+ case 1: {
+ force_assembler.emplaceElement(
+ element, ForceElement::Force{{0., 1e-6, 0., 0., 0., 0.}});
+ break;
+ }
+ }
+ }
+ // 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 clamp the beam completely (i.e. all dofs) at (0,0)
+
+ std::vector<BoundaryCondition> boundary_conditions;
+ for (const auto &vertex : mesh.localVertices()) {
+ switch (vertex.index()) {
+ case 0: {
+ // 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.});
+ // The displacement in z direction is zero.
+ boundary_conditions.push_back({vertex, 2, 0.});
+ // The rotation around x axis is zero.
+ boundary_conditions.push_back({vertex, 3, 0.});
+ // The rotation around y axis is zero.
+ boundary_conditions.push_back({vertex, 4, 0.});
+ // The rotation around z axis is zero.
+ boundary_conditions.push_back({vertex, 5, 0.});
+ break;
+ }
+ }
+ }
+
+ // We are ready to minimize the energy.
+
+ auto result = solver.computeSolution(
+ boundary_conditions, Vector::fromGlobalMesh(mesh), time, &assembler);
+
+ // Let us now 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);
+
+ if (Policy::isPrimaryRank())
+ global_result.unwrap().print();
+}
diff --git a/solve/test/boundaryConditionsToEquations_Test.cc b/solve/test/boundaryConditionsToEquations_Test.cc
index bda8cbfa0576df759a3d1478fc7ba7d1f1a18f73..81e286a9d03937f8ab14a03cb28d7014752200a4 100644
--- a/solve/test/boundaryConditionsToEquations_Test.cc
+++ b/solve/test/boundaryConditionsToEquations_Test.cc
@@ -23,6 +23,7 @@
#include <gmock/gmock.h>
using ae108::cppptest::ScalarEqIfLocal;
+using ae108::cppptest::ScalarNearIfLocal;
using testing::DoubleEq;
using testing::Eq;
using testing::Ge;
@@ -176,10 +177,10 @@ TYPED_TEST(boundaryConditionsToEquations_Test,
const auto image = apply(result, displacements);
for (auto row = size_type{0}; row < image.unwrap().size(); ++row) {
- EXPECT_THAT(image.unwrap(), ScalarEqIfLocal(row, 0.));
+ EXPECT_THAT(image.unwrap(), ScalarNearIfLocal(row, 0., 1e-17));
}
}
} // namespace
} // namespace solve
-} // namespace ae108
\ No newline at end of file
+} // namespace ae108
diff --git a/tests/examples/cantilever_beam/definition.json b/tests/examples/cantilever_beam/definition.json
new file mode 100644
index 0000000000000000000000000000000000000000..c73a3fce5adc1f1ed2c4e18fd5aa621607f18823
--- /dev/null
+++ b/tests/examples/cantilever_beam/definition.json
@@ -0,0 +1,11 @@
+{
+ "executable": [
+ "examples",
+ "ae108-examples-CantileverBeam"
+ ],
+ "compare_stdout": "numeric",
+ "mpi_processes": [
+ 1,
+ 3
+ ]
+}
\ No newline at end of file
diff --git a/tests/examples/cantilever_beam/references/stdout.txt b/tests/examples/cantilever_beam/references/stdout.txt
new file mode 100644
index 0000000000000000000000000000000000000000..679b240d9604ac8ac4ebb17d382274cc103805bb
--- /dev/null
+++ b/tests/examples/cantilever_beam/references/stdout.txt
@@ -0,0 +1,14 @@
+Vec Object: 1 MPI processes
+ type: seq
+0.
+0.
+0.
+0.
+0.
+0.
+0.
+-0.0682796
+0.
+0.
+0.
+-0.101859