Commit 2408a8d4 authored by webmanue's avatar webmanue
Browse files

Merge remote-tracking branch 'origin/master' into 98-update-reference-platform-to-ubuntu-jammy

parents 10ec14bb 0b3a29ac
Pipeline #135496 passed with stages
in 18 minutes and 42 seconds
......@@ -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
......
......@@ -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 */,
&section));
return makeUniqueEntity<Policy>(section);
}
Policy::handleError(DMSetSection(_mesh.get(), section));
// create default global section
Policy::handleError(DMGetGlobalSection(_mesh.get(), &section));
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) {
......
......@@ -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
......@@ -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
// © 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();
}
......@@ -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
{
"executable": [
"examples",
"ae108-examples-CantileverBeam"
],
"compare_stdout": "numeric",
"mpi_processes": [
1,
3
]
}
\ No newline at end of file
Vec Object: 1 MPI processes
type: seq
0.
0.
0.
0.
0.
0.
0.
-0.0682796
0.
0.
0.
-0.101859
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment