Commit 37a189ab authored by webmanue's avatar webmanue
Browse files

Merge branch '115-add-bar-element-2d-3d' into 'master'

Resolve "Add bar element (2D, 3D)"

Closes #115

See merge request mechanics-and-materials/ae108!104
parents 0b084769 77aaf5ce
Pipeline #138133 passed with stages
in 13 minutes and 58 seconds
// © 2022 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/elements/Bar.h"
\ No newline at end of file
......@@ -17,6 +17,7 @@ cmake_minimum_required(VERSION 3.12 FATAL_ERROR)
add_library(${PROJECT_NAME}-elements
AutomaticForcesTrait.cc
AutomaticStiffnessMatrixTrait.cc
Bar.cc
ComputeConsistentMassMatrixTrait.cc
ComputeEnergyTrait.cc
ComputeForcesTrait.cc
......
// © 2021, 2022 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.
#pragma once
#include "ae108/elements/ComputeEnergyTrait.h"
#include "ae108/elements/ComputeForcesTrait.h"
#include "ae108/elements/ComputeStiffnessMatrixTrait.h"
#include "ae108/elements/ElementBase.h"
#include "ae108/elements/tensor/as_vector.h"
namespace ae108 {
namespace elements {
template <class RealType_> struct BarProperties {
using real_type = RealType_;
real_type young_modulus;
real_type cross_section;
};
/**
* @brief Computes the stiffness matrix for a bar with the given
* axis and the given properties.
* @note see Dennis M. Kochmann, "Computational Mechanics I: Introduction to
* FEA Analysis", https://mm.ethz.ch/education, pp.77-81
*/
template <std::size_t Dimension_, class RealType_>
Eigen::Matrix<RealType_, 2 * Dimension_, 2 * Dimension_, Eigen::RowMajor>
bar_stiffness_matrix(const tensor::Tensor<RealType_, Dimension_> &axis,
const BarProperties<RealType_> &properties) noexcept {
static_assert(Dimension_ > 0);
/**
* @brief Computes the 1D stiffness matrix assuming an axis of unit length.
*/
constexpr auto stiffness_matrix_1D =
[](const BarProperties<RealType_> &properties) {
const auto factor = properties.cross_section * properties.young_modulus;
auto x = Eigen::Matrix<double, 2, 2>();
x << factor, -1. * factor, -1 * factor, factor;
return x;
};
/**
* @brief Computes the rotation matrix assuming an axis of unit length.
*/
constexpr auto rotation =
[](const tensor::Tensor<RealType_, Dimension_> &axis) {
auto x = Eigen::Matrix<RealType_, 2, 2 * Dimension_>::Zero().eval();
x.template block<1, Dimension_>(0, 0) = tensor::as_vector(&axis);
x.template block<1, Dimension_>(1, Dimension_) =
tensor::as_vector(&axis);
return x;
};
const auto K = stiffness_matrix_1D(properties);
const auto R = rotation(axis);
return std::pow(tensor::as_vector(&axis).norm(), -3.) *
(R.transpose() * K * R);
}
/**
* @tparam Dimension_ The dimension of physical space.
*/
template <std::size_t Dimension_, class ValueType_ = double,
class RealType_ = double>
struct Bar final
: ElementBase<Bar<Dimension_, ValueType_, RealType_>, std::size_t,
ValueType_, RealType_, 2, Dimension_> {
public:
/**
* @param matrix A symmetric stiffness matrix.
*/
explicit Bar(typename Bar::StiffnessMatrix matrix) noexcept
: stiffness_matrix_(std::move(matrix)) {}
const typename Bar::StiffnessMatrix &stiffness_matrix() const {
return stiffness_matrix_;
}
/**
* @brief The dimension of physical space.
*/
static constexpr typename Bar::size_type dimension() { return Dimension_; }
private:
typename Bar::StiffnessMatrix stiffness_matrix_;
};
template <std::size_t Dimension_, class ValueType_, class RealType_>
struct ComputeEnergyTrait<Bar<Dimension_, ValueType_, RealType_>> {
template <class Element>
typename Element::Energy
operator()(const Element &element,
const typename Element::NodalDisplacements &u,
const typename Element::Time &) const noexcept {
const auto v = tensor::as_vector(&u);
return typename Element::Energy{.5} * v.transpose() *
element.stiffness_matrix() * v;
}
};
template <std::size_t Dimension_, class ValueType_, class RealType_>
struct ComputeForcesTrait<Bar<Dimension_, ValueType_, RealType_>> {
template <class Element>
typename Element::Forces
operator()(const Element &element,
const typename Element::NodalDisplacements &u,
const typename Element::Time &) const noexcept {
typename Element::Forces forces;
tensor::as_vector(&forces) =
element.stiffness_matrix() * tensor::as_vector(&u);
return forces;
}
};
template <std::size_t Dimension_, class ValueType_, class RealType_>
struct ComputeStiffnessMatrixTrait<Bar<Dimension_, ValueType_, RealType_>> {
template <class Element>
typename Element::StiffnessMatrix
operator()(const Element &element,
const typename Element::NodalDisplacements &,
const typename Element::Time &) const noexcept {
return element.stiffness_matrix();
}
};
} // namespace elements
} // namespace ae108
\ No newline at end of file
// © 2021, 2022 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 "Element_Test.h"
#include "ae108/elements/Bar.h"
#include <gmock/gmock.h>
#include <numeric>
using testing::DoubleEq;
using testing::DoubleNear;
using testing::Test;
using testing::Types;
namespace ae108 {
namespace elements {
namespace {
template <class Element>
tensor::Tensor<double, Element::dimension()>
create_reference_element_axis() noexcept {
auto element_axis = tensor::Tensor<double, Element::dimension()>();
element_axis[0] = 1.;
return element_axis;
}
template <class Element>
tensor::Tensor<double, Element::dimension()>
create_rotated_and_stretched_element_axis() noexcept {
auto element_axis = tensor::Tensor<double, Element::dimension()>();
std::iota(element_axis.begin(), element_axis.end(), 3.);
return element_axis;
}
constexpr static BarProperties<double> properties{.3, .4};
template <class Element_> struct ReferenceConfiguration {
using Element = Element_;
static Element create_element() noexcept {
return Element(bar_stiffness_matrix<Element::dimension()>(
create_reference_element_axis<Element_>(), properties));
}
static typename Element::Time create_time() noexcept {
return typename Element::Time{0.};
}
};
template <class Element_> struct RotatedAndStretchedConfiguration {
using Element = Element_;
static Element create_element() noexcept {
return Element(bar_stiffness_matrix<Element::dimension()>(
create_rotated_and_stretched_element_axis<Element_>(), properties));
}
static typename Element::Time create_time() noexcept {
return typename Element::Time{0.};
}
};
using Configurations =
Types<ReferenceConfiguration<Bar<1>>, ReferenceConfiguration<Bar<2>>,
ReferenceConfiguration<Bar<3>>,
RotatedAndStretchedConfiguration<Bar<1>>,
RotatedAndStretchedConfiguration<Bar<2>>,
RotatedAndStretchedConfiguration<Bar<3>>>;
INSTANTIATE_TYPED_TEST_CASE_P(Bar_Test, Element_Test, Configurations);
struct Bar1D_Test : Test {
using Element = Bar<1>;
const Element element = ReferenceConfiguration<Element>::create_element();
};
TEST_F(Bar1D_Test, computes_energy_with_axial_displacement) {
const auto time = Element::Time{0.};
const Element::NodalDisplacements displacements = {{
{{0.}},
{{1.}},
}};
EXPECT_THAT(
element.computeEnergy(displacements, time),
DoubleEq(.5 * properties.cross_section * properties.young_modulus));
}
TEST_F(Bar1D_Test, computes_energy_with_double_axial_displacement) {
const auto time = Element::Time{0.};
const Element::NodalDisplacements displacements = {{
{{0.}},
{{2.}},
}};
EXPECT_THAT(
element.computeEnergy(displacements, time),
DoubleEq(.5 * 4. * properties.cross_section * properties.young_modulus));
}
struct Bar1D_rotated_Test : Test {
using Element = Bar<1>;
const Element element =
RotatedAndStretchedConfiguration<Element>::create_element();
const tensor::Tensor<double, 1> axis =
create_rotated_and_stretched_element_axis<Element>();
const double L = tensor::as_vector(&axis).norm();
};
TEST_F(Bar1D_rotated_Test, computes_energy_with_axial_displacement) {
const auto time = Element::Time{0.};
const Element::NodalDisplacements displacements = {{
{{0.}},
{{1.}},
}};
EXPECT_THAT(
element.computeEnergy(displacements, time),
DoubleEq(.5 * properties.cross_section * properties.young_modulus / L));
}
TEST_F(Bar1D_rotated_Test, computes_energy_with_double_axial_displacement) {
const auto time = Element::Time{0.};
const Element::NodalDisplacements displacements = {{
{{0.}},
{{2.}},
}};
EXPECT_THAT(element.computeEnergy(displacements, time),
DoubleEq(.5 * 4. * properties.cross_section *
properties.young_modulus / L));
}
struct Bar2D_Test : Test {
using Element = Bar<2>;
const Element element = ReferenceConfiguration<Element>::create_element();
};
TEST_F(Bar2D_Test, computes_energy_with_axial_displacement) {
const auto time = Element::Time{0.};
const Element::NodalDisplacements displacements = {{
{{0., 0.}},
{{1., 0.}},
}};
EXPECT_THAT(
element.computeEnergy(displacements, time),
DoubleEq(.5 * properties.cross_section * properties.young_modulus));
}
TEST_F(Bar2D_Test, computes_energy_with_double_axial_displacement) {
const auto time = Element::Time{0.};
const Element::NodalDisplacements displacements = {{
{{0., 0.}},
{{2., 0.}},
}};
EXPECT_THAT(
element.computeEnergy(displacements, time),
DoubleEq(.5 * 4. * properties.cross_section * properties.young_modulus));
}
TEST_F(Bar2D_Test, computes_energy_with_lateral_displacement) {
const auto time = Element::Time{0.};
const Element::NodalDisplacements displacements = {{
{{0., 0.}},
{{0., 1.}},
}};
EXPECT_THAT(element.computeEnergy(displacements, time), DoubleEq(0.));
}
TEST_F(Bar2D_Test,
computes_energy_with_combined_axial_and_lateral_displacement) {
const auto time = Element::Time{0.};
const Element::NodalDisplacements displacements = {{
{{0., 0.}},
{{1., 1.}},
}};
EXPECT_THAT(
element.computeEnergy(displacements, time),
DoubleEq(.5 * properties.cross_section * properties.young_modulus));
}
struct Bar3D_Test : Test {
using Element = Bar<3>;
const Element element = ReferenceConfiguration<Element>::create_element();
const tensor::Tensor<double, 3> axis =
create_reference_element_axis<Element>();
const double L = tensor::as_vector(&axis).norm();
};
TEST_F(Bar3D_Test, computes_energy_with_axial_displacement) {
const auto time = Element::Time{0.};
const Element::NodalDisplacements displacements = {{
{{0., 0., 0.}},
{{axis[0] / L, axis[1] / L, axis[2] / L}},
}};
EXPECT_THAT(
element.computeEnergy(displacements, time),
DoubleEq(.5 * properties.cross_section * properties.young_modulus));
}
TEST_F(Bar3D_Test, computes_energy_with_double_axial_displacement) {
const auto time = Element::Time{0.};
const Element::NodalDisplacements displacements = {{
{{0., 0., 0.}},
{{2., 0., 0.}},
}};
EXPECT_THAT(
element.computeEnergy(displacements, time),
DoubleEq(.5 * 4. * properties.cross_section * properties.young_modulus));
}
TEST_F(Bar3D_Test, computes_energy_with_lateral_displacement) {
const auto time = Element::Time{0.};
const Element::NodalDisplacements displacements = {{
{{0., 0., 0.}},
{{0., 1., 0.}},
}};
EXPECT_THAT(element.computeEnergy(displacements, time), DoubleEq(0.));
}
TEST_F(Bar3D_Test,
computes_energy_with_combined_axial_and_lateral_displacement) {
const auto time = Element::Time{0.};
const Element::NodalDisplacements displacements = {{
{{0., 0., 0.}},
{{1., 1., 0.}},
}};
EXPECT_THAT(
element.computeEnergy(displacements, time),
DoubleEq(.5 * properties.cross_section * properties.young_modulus));
}
struct Bar2D_rotated_Test : Test {
using Element = Bar<2>;
const Element element =
RotatedAndStretchedConfiguration<Element>::create_element();
const tensor::Tensor<double, 2> axis =
create_rotated_and_stretched_element_axis<Element>();
const double L = tensor::as_vector(&axis).norm();
};
TEST_F(Bar2D_rotated_Test, computes_energy_with_axial_displacement) {
const auto time = Element::Time{0.};
const Element::NodalDisplacements displacements = {{
{{0., 0.}},
{{axis[0] / L, axis[1] / L}},
}};
EXPECT_THAT(
element.computeEnergy(displacements, time),
DoubleEq(.5 * properties.cross_section * properties.young_modulus / L));
}
TEST_F(Bar2D_rotated_Test, computes_energy_with_double_axial_displacement) {
const auto time = Element::Time{0.};
const Element::NodalDisplacements displacements = {{
{{0., 0.}},
{{2. * axis[0] / L, 2. * axis[1] / L}},
}};
EXPECT_THAT(element.computeEnergy(displacements, time),
DoubleEq(.5 * 4. * properties.cross_section *
properties.young_modulus / L));
}
TEST_F(Bar2D_rotated_Test, computes_energy_with_lateral_displacement) {
const auto time = Element::Time{0.};
const Element::NodalDisplacements displacements = {{
{{0., 0.}},
{{axis[1] / L, -axis[0] / L}},
}};
EXPECT_THAT(element.computeEnergy(displacements, time),
DoubleNear(0., 1e-18));
}
TEST_F(Bar2D_rotated_Test,
computes_energy_with_combined_axial_and_lateral_displacement) {
const auto time = Element::Time{0.};
const Element::NodalDisplacements displacements = {{
{{0., 0.}},
{{(axis[0] + axis[1]) / L, (axis[1] - axis[0]) / L}},
}};
EXPECT_THAT(
element.computeEnergy(displacements, time),
DoubleEq(.5 * properties.cross_section * properties.young_modulus / L));
}
struct Bar3D_rotated_Test : Test {
using Element = Bar<3>;
const Element element =
RotatedAndStretchedConfiguration<Element>::create_element();
const tensor::Tensor<double, 3> axis =
create_rotated_and_stretched_element_axis<Element>();
const double L = tensor::as_vector(&axis).norm();
};
TEST_F(Bar3D_rotated_Test, computes_energy_with_axial_displacement) {
const auto time = Element::Time{0.};
const Element::NodalDisplacements displacements = {{
{{0., 0., 0.}},
{{axis[0] / L, axis[1] / L, axis[2] / L}},
}};
EXPECT_THAT(
element.computeEnergy(displacements, time),
DoubleEq(.5 * properties.cross_section * properties.young_modulus / L));
}
TEST_F(Bar3D_rotated_Test, computes_energy_with_double_axial_displacement) {
const auto time = Element::Time{0.};
const Element::NodalDisplacements displacements = {{
{{0., 0., 0.}},
{{2. * axis[0] / L, 2. * axis[1] / L, 2. * axis[2] / L}},
}};
EXPECT_THAT(element.computeEnergy(displacements, time),
DoubleEq(.5 * 4. * properties.cross_section *
properties.young_modulus / L));
}
TEST_F(Bar3D_rotated_Test, computes_energy_with_lateral_displacement) {
const auto time = Element::Time{0.};
const Element::NodalDisplacements displacements = {{
{{0., 0., 0.}},
{{-axis[1] / L, axis[0] / L, 0.}},
}};
EXPECT_THAT(element.computeEnergy(displacements, time),
DoubleNear(0., 1e-18));
}
TEST_F(Bar3D_rotated_Test,
computes_energy_with_combined_axial_and_lateral_displacement) {
const auto time = Element::Time{0.};
const Element::NodalDisplacements displacements = {{
{{0., 0., 0.}},
{{(axis[0] - axis[1]) / L, (axis[1] + axis[0]) / L, axis[2] / L}},
}};
EXPECT_THAT(
element.computeEnergy(displacements, time),
DoubleEq(.5 * properties.cross_section * properties.young_modulus / L));
}
} // namespace
} // namespace elements
} // namespace ae108
\ No newline at end of file
......@@ -15,6 +15,7 @@
cmake_minimum_required(VERSION 3.12 FATAL_ERROR)
add_executable(${PROJECT_NAME}-ElementsTests
Bar_Test.cc
CoreElement_Test.cc
ForceElement_Test.cc
TimoshenkoBeamElement_Test.cc
......
// © 2020, 2022 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/elements/Bar.h"
#include "ae108/assembly/Assembler.h"
#include "ae108/cpppetsc/Context.h"
#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 "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 Viewer = cpppetsc::Viewer<Policy>;
using BoundaryCondition = cpppetsc::MeshBoundaryCondition<Mesh>;
// In this example we will simulate six elements, A, B, C, D, E, F.
// Each of these elements has two vertices.
// 3---C---2
// | \ / |
// D E/F B
// | / \ |
// 0---A---1