Skip to content

Commit

Permalink
add explicit derivatives of KL shell (#225)
Browse files Browse the repository at this point in the history
  • Loading branch information
rath3t committed Jan 10, 2024
1 parent 8837c21 commit 2dfc766
Show file tree
Hide file tree
Showing 25 changed files with 861 additions and 249 deletions.
80 changes: 80 additions & 0 deletions .github/workflows/docsDryRun.yml
Original file line number Diff line number Diff line change
@@ -0,0 +1,80 @@
# SPDX-FileCopyrightText: 2021-2024 The Ikarus Developers [email protected]
# SPDX-License-Identifier: CC0-1.0

name: Dry Run Documentation

on:
push:
paths-ignore:
- 'docs/**'
- '.github/workflows/ghpages.yml'
- '.github/workflows/createDockerContainer.yml'
- '**.md'

pull_request:
types: [opened]
branches:
- main
paths-ignore:
- 'docs/**'
- '.github/workflows/ghpages.yml'
- '.github/workflows/createDockerContainer.yml'
- '**.md'

jobs:
Build-Docs-Dry-Run:
runs-on: ubuntu-latest
container:
image: ikarusproject/ikarus-dev:latest
options: --memory-swap="20g" --memory="20g" --cpus="2" --user root
steps:
- uses: actions/checkout@v2
with:
path: 'repo'
- uses: actions/setup-python@v4
with:
python-version: 3.x
- name: Setup Mkdocs
run: |
pip install mkdocs
pip install git+https://${{ secrets.MKDOCS_TOKEN }}@github.com/squidfunk/mkdocs-material-insiders.git
pip install mkdocs-macros-plugin
pip install mkdocs-drawio-exporter
pip install mkdocs-glightbox
pip install mike
pip install mkdocs-bibtex
pip install pillow cairosvg
apt-get update
apt-get install libcairo2-dev libfreetype6-dev libffi-dev libjpeg-dev libpng-dev libz-dev -y
#wget https://github.com/jgraph/drawio-desktop/releases/download/v16.5.1/drawio-amd64-16.5.1.deb
#sudo apt-get install libayatana-appindicator3-1 libnotify4
#sudo dpkg -i drawio-amd64-16.5.1.deb
#sudo apt-get -y -f install
#sudo apt install libasound2 xvfb
- name: Build Doxygen
run: |
git clone https://github.com/jothepro/doxygen-awesome-css.git
cd doxygen-awesome-css
make install
apt-get update
apt-get install texlive-base -y
cd ..
cd repo
mkdir build_docs
cd build_docs
cmake .. -DBUILD_DOCS=1 -DDUNE_ENABLE_PYTHONBINDINGS=0
cmake --build . --target doxygen_ikarus
- name: Build Website
run: |
git clone https://github.com/ikarus-project/ikarus-project.github.io.git
cd ikarus-project.github.io
export MKDOCSOFFLINE=false
git fetch origin gh-pages --depth=1
cp -R ../repo/build_docs/ .
cp -R ../repo/ikarus/ .
git config --local user.email "${{ github.event.head_commit.author.email }}"
git config --local user.name "${{ github.event.head_commit.author.name }}"
cd build_docs/docs
mike deploy dev --config-file mkdocs.insiders.yml
4 changes: 2 additions & 2 deletions codespellignore
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
# SPDX-FileCopyrightText: 2021-2024 The Ikarus Developers [email protected]
# SPDX-License-Identifier: CC0-1.0

tE
te
nd
4 changes: 4 additions & 0 deletions docs/website/doxygen/Doxylocal
Original file line number Diff line number Diff line change
Expand Up @@ -112,3 +112,7 @@ DIR_GRAPH_MAX_DEPTH = 5
GENERATE_LEGEND = YES
DOT_MULTI_TARGETS = YES
UML_LOOK = NO
WARN_NO_PARAMDOC =YES
WARN_IF_INCOMPLETE_DOC=YES
WARN_IF_UNDOCUMENTED=YES
WARN_AS_ERROR=NO
2 changes: 1 addition & 1 deletion ikarus/finiteelements/mechanics/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,6 @@
add_subdirectory(materials)
# install headers
install(FILES nonlinearelastic.hh enhancedassumedstrains.hh linearelastic.hh materials.hh
kirchhoffloveshell.hh
kirchhoffloveshell.hh membranestrains.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/ikarus/finiteelements/mechanics
)
427 changes: 332 additions & 95 deletions ikarus/finiteelements/mechanics/kirchhoffloveshell.hh

Large diffs are not rendered by default.

13 changes: 6 additions & 7 deletions ikarus/finiteelements/mechanics/linearelastic.hh
Original file line number Diff line number Diff line change
Expand Up @@ -233,10 +233,10 @@ namespace Ikarus {
Dune::CachedLocalBasis<
std::remove_cvref_t<decltype(std::declval<LocalView>().tree().child(0).finiteElement().localBasis())>>
localBasis;
std::function<Eigen::Vector<double, Traits::worlddim>(const Eigen::Vector<double, Traits::worlddim>&,
std::function<Eigen::Vector<double, Traits::worlddim>(const Dune::FieldVector<double, Traits::worlddim>&,
const double&)>
volumeLoad;
std::function<Eigen::Vector<double, Traits::worlddim>(const Eigen::Vector<double, Traits::worlddim>&,
std::function<Eigen::Vector<double, Traits::worlddim>(const Dune::FieldVector<double, Traits::worlddim>&,
const double&)>
neumannBoundaryLoad;
const BoundaryPatch<GridView>* neumannBoundary;
Expand Down Expand Up @@ -270,7 +270,7 @@ namespace Ikarus {
if (volumeLoad) {
for (const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) {
const auto uVal = u.evaluate(gpIndex);
Eigen::Vector<double, Traits::worlddim> fext = volumeLoad(toEigen(geo.global(gp.position())), lambda);
Eigen::Vector<double, Traits::worlddim> fext = volumeLoad(geo.global(gp.position()), lambda);
energy -= uVal.dot(fext) * geo.integrationElement(gp.position()) * gp.weight();
}
}
Expand All @@ -294,7 +294,7 @@ namespace Ikarus {
const auto uVal = u.evaluate(quadPos);

// Value of the Neumann data at the current position
auto neumannValue = neumannBoundaryLoad(toEigen(intersection.geometry().global(curQuad.position())), lambda);
auto neumannValue = neumannBoundaryLoad(intersection.geometry().global(curQuad.position()), lambda);

energy -= neumannValue.dot(uVal) * curQuad.weight() * integrationElement;
}
Expand Down Expand Up @@ -327,7 +327,7 @@ namespace Ikarus {
if (volumeLoad) {
const auto u = getDisplacementFunction(par, dx);
for (const auto& [gpIndex, gp] : u.viewOverIntegrationPoints()) {
Eigen::Vector<double, Traits::worlddim> fext = volumeLoad(toEigen(geo.global(gp.position())), lambda);
Eigen::Vector<double, Traits::worlddim> fext = volumeLoad(geo.global(gp.position()), lambda);
for (size_t i = 0; i < numberOfNodes; ++i) {
const auto udCi = u.evaluateDerivative(gpIndex, wrt(coeff(i)));
force.template segment<myDim>(myDim * i)
Expand Down Expand Up @@ -357,8 +357,7 @@ namespace Ikarus {
const auto udCi = u.evaluateDerivative(quadPos, wrt(coeff(i)));

// Value of the Neumann data at the current position
auto neumannValue
= neumannBoundaryLoad(toEigen(intersection.geometry().global(curQuad.position())), lambda);
auto neumannValue = neumannBoundaryLoad(intersection.geometry().global(curQuad.position()), lambda);
force.template segment<myDim>(myDim * i) -= udCi * neumannValue * curQuad.weight() * integrationElement;
}
}
Expand Down
109 changes: 109 additions & 0 deletions ikarus/finiteelements/mechanics/membranestrains.hh
Original file line number Diff line number Diff line change
@@ -0,0 +1,109 @@
// SPDX-FileCopyrightText: 2021-2024 The Ikarus Developers [email protected]
// SPDX-License-Identifier: LGPL-3.0-or-later

/**
* \file membranestrains.hh
* \brief Implementation of membrane strain for shells
*/

#pragma once
#include <ranges>

#include <dune/common/fvector.hh>

#include <Eigen/Core>
namespace Ikarus {

struct DefaultMembraneStrain {
/**
* \brief Compute the strain vector at a given integration point.
*
* \param gpPos The position of the integration point.
* \param geo The geometry object providing the transposed Jacobian.
* \param uFunction The function representing the displacement field.
*
* \tparam Geometry The type of the geometry object.
*
* \return The strain vector at the given integration point.
*/
template <typename Geometry>
auto value(const Dune::FieldVector<double, 2>& gpPos, const Geometry& geo, const auto& uFunction) const
-> Eigen::Vector3<typename std::remove_cvref_t<decltype(uFunction)>::ctype> {
using ScalarType = typename std::remove_cvref_t<decltype(uFunction)>::ctype;
Eigen::Vector3<ScalarType> epsV;
const auto J = Dune::toEigen(geo.jacobianTransposed(gpPos));
using namespace Dune;
using namespace Dune::DerivativeDirections;
const Eigen::Matrix<ScalarType, 3, 2> gradu = toEigen(
uFunction.evaluateDerivative(gpPos, // Here the gpIndex could be passed
Dune::wrt(spatialAll), Dune::on(Dune::DerivativeDirections::referenceElement)));
const Eigen::Matrix<ScalarType, 2, 3> j = J + gradu.transpose();

epsV << J.row(0).dot(gradu.col(0)) + 0.5 * gradu.col(0).squaredNorm(),
J.row(1).dot(gradu.col(1)) + 0.5 * gradu.col(1).squaredNorm(),
j.row(0).dot(j.row(1)) - J.row(0).dot(J.row(1));
return epsV;
}

/**
* \brief Compute the strain-displacement matrix for a given node and integration point.
*
* \param gpPos The position of the integration point.
* \param jcur The Jacobian matrix.
* \param dNAtGp The derivative of the shape functions at the integration point.
* \param geo The geometry object of the finite element.
* \param uFunction The function representing the displacement field.
* \param localBasis The local basis object.
* \param node The FE node index.
*
* \tparam Geometry The type of the geometry object.
* \tparam ScalarType The scalar type for the matrix elements.
*
* \return The strain-displacement matrix for the given node and integration point.
*/
template <typename Geometry, typename ScalarType>
auto derivative(const Dune::FieldVector<double, 2>& gpPos, const Eigen::Matrix<ScalarType, 2, 3>& jcur,
const auto& dNAtGp, const Geometry& geo, const auto& uFunction, const auto& localBasis,
const int node) const {
Eigen::Matrix<ScalarType, 3, 3> bop;
bop.row(0) = jcur.row(0) * dNAtGp(node, 0);
bop.row(1) = jcur.row(1) * dNAtGp(node, 1);
bop.row(2) = jcur.row(0) * dNAtGp(node, 1) + jcur.row(1) * dNAtGp(node, 0);

return bop;
}

/**
* \brief Compute the second derivative of the membrane strain for a given node pair and integration point.
* \details This function computes the geometric tangent stiffness for a given node pair at a given integration
* point.
*
* \param gpPos The position of the integration point.
* \param dNAtGp The derivative of the shape functions at the integration point.
* \param geo The geometry object.
* \param uFunction The function representing the displacement field.
* \param localBasis The local basis object.
* \param S The strain vector.
* \param I The index of the first node.
* \param J The index of the second node.
*
* \tparam Geometry The type of the geometry object.
* \tparam ScalarType The scalar type for the matrix elements.
*
* \return The second derivative of the membrane strain.
*/
template <typename Geometry, typename ScalarType>
auto secondDerivative(const Dune::FieldVector<double, 2>& gpPos, const auto& dNAtGp, const Geometry& geo,
const auto& uFunction, const auto& localBasis, const Eigen::Vector3<ScalarType>& S, int I,
int J) const {
const auto& dN1i = dNAtGp(I, 0);
const auto& dN1j = dNAtGp(J, 0);
const auto& dN2i = dNAtGp(I, 1);
const auto& dN2j = dNAtGp(J, 1);
const ScalarType NS = dN1i * dN1j * S[0] + dN2i * dN2j * S[1] + (dN1i * dN2j + dN2i * dN1j) * S[2];
Eigen::Matrix<ScalarType, 3, 3> kg = Eigen::Matrix<double, 3, 3>::Identity() * NS;
return kg;
}
};

} // namespace Ikarus
14 changes: 6 additions & 8 deletions ikarus/finiteelements/mechanics/nonlinearelastic.hh
Original file line number Diff line number Diff line change
Expand Up @@ -275,10 +275,10 @@ namespace Ikarus {
Dune::CachedLocalBasis<
std::remove_cvref_t<decltype(std::declval<LocalView>().tree().child(0).finiteElement().localBasis())>>
localBasis;
std::function<Eigen::Vector<double, Traits::worlddim>(const Eigen::Vector<double, Traits::worlddim>&,
std::function<Eigen::Vector<double, Traits::worlddim>(const Dune::FieldVector<double, Traits::worlddim>&,
const double&)>
volumeLoad;
std::function<Eigen::Vector<double, Traits::worlddim>(const Eigen::Vector<double, Traits::worlddim>&,
std::function<Eigen::Vector<double, Traits::worlddim>(const Dune::FieldVector<double, Traits::worlddim>&,
const double&)>
neumannBoundaryLoad;
const BoundaryPatch<GridView>* neumannBoundary;
Expand Down Expand Up @@ -308,7 +308,7 @@ namespace Ikarus {
if (volumeLoad) {
for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) {
const auto u = uFunction.evaluate(gpIndex);
const Eigen::Vector<double, Traits::worlddim> fExt = volumeLoad(toEigen(geo.global(gp.position())), lambda);
const Eigen::Vector<double, Traits::worlddim> fExt = volumeLoad(geo.global(gp.position()), lambda);
energy -= u.dot(fExt) * geo.integrationElement(gp.position()) * gp.weight();
}
}
Expand All @@ -333,8 +333,7 @@ namespace Ikarus {
const auto u = uFunction.evaluate(quadPos);

// Value of the Neumann data at the current position
const auto neumannValue
= neumannBoundaryLoad(toEigen(intersection.geometry().global(curQuad.position())), lambda);
const auto neumannValue = neumannBoundaryLoad(intersection.geometry().global(curQuad.position()), lambda);
energy -= neumannValue.dot(u) * curQuad.weight() * intElement;
}
}
Expand Down Expand Up @@ -366,7 +365,7 @@ namespace Ikarus {
const auto u = displacementFunction(par, dx);
for (const auto& [gpIndex, gp] : u.viewOverIntegrationPoints()) {
const double intElement = geo.integrationElement(gp.position()) * gp.weight();
const Eigen::Vector<double, Traits::worlddim> fExt = volumeLoad(toEigen(geo.global(gp.position())), lambda);
const Eigen::Vector<double, Traits::worlddim> fExt = volumeLoad(geo.global(gp.position()), lambda);
for (size_t i = 0; i < numberOfNodes; ++i) {
const auto udCi = u.evaluateDerivative(gpIndex, wrt(coeff(i)));
force.template segment<myDim>(myDim * i) -= udCi * fExt * intElement;
Expand Down Expand Up @@ -395,8 +394,7 @@ namespace Ikarus {
const auto udCi = u.evaluateDerivative(quadPos, wrt(coeff(i)));

// Value of the Neumann data at the current position
const auto neumannValue
= neumannBoundaryLoad(toEigen(intersection.geometry().global(curQuad.position())), lambda);
const auto neumannValue = neumannBoundaryLoad(intersection.geometry().global(curQuad.position()), lambda);
force.template segment<myDim>(myDim * i) -= udCi * neumannValue * curQuad.weight() * intElement;
}
}
Expand Down
4 changes: 2 additions & 2 deletions ikarus/python/finiteelements/kirchhoffloveshell.hh
Original file line number Diff line number Diff line change
Expand Up @@ -51,8 +51,8 @@ namespace Ikarus::Python {
using Traits = typename KirchhoffLoveShell::Traits;
using FERequirements = typename KirchhoffLoveShell::FERequirementType;

using LoadFunction = std::function<Eigen::Vector<double, Traits::worlddim>(Eigen::Vector<double, Traits::worlddim>,
const double&)>;
using LoadFunction = std::function<Eigen::Vector<double, Traits::worlddim>(
Dune::FieldVector<double, Traits::worlddim>, const double&)>;
cls.def(pybind11::init([](const GlobalBasis& basis, const Element& element, double emod, double nu,
double thickness, const LoadFunction volumeLoad) {
return new KirchhoffLoveShell(basis, element, emod, nu, thickness, volumeLoad);
Expand Down
4 changes: 2 additions & 2 deletions ikarus/python/finiteelements/nonlinearelastic.hh
Original file line number Diff line number Diff line change
Expand Up @@ -39,8 +39,8 @@ namespace Ikarus::Python {
}),
pybind11::keep_alive<1, 2>(), pybind11::keep_alive<1, 3>());

using LoadFunction = std::function<Eigen::Vector<double, Traits::worlddim>(Eigen::Vector<double, Traits::worlddim>,
const double&)>;
using LoadFunction = std::function<Eigen::Vector<double, Traits::worlddim>(
Dune::FieldVector<double, Traits::worlddim>, const double&)>;

cls.def(pybind11::init(
[](const GlobalBasis& basis, const Element& element, const Material& mat,
Expand Down
4 changes: 2 additions & 2 deletions ikarus/python/finiteelements/registerelement.hh
Original file line number Diff line number Diff line change
Expand Up @@ -59,8 +59,8 @@ namespace Ikarus::Python {
}),
pybind11::keep_alive<1, 2>(), pybind11::keep_alive<1, 3>());

using LoadFunction = std::function<Eigen::Vector<double, Traits::worlddim>(Eigen::Vector<double, Traits::worlddim>,
const double&)>;
using LoadFunction = std::function<Eigen::Vector<double, Traits::worlddim>(
Dune::FieldVector<double, Traits::worlddim>, const double&)>;
if constexpr (defaultInitializers)
cls.def(
pybind11::init([](const GlobalBasis& basis, const Element& element, double emod, double nu,
Expand Down
1 change: 1 addition & 0 deletions ikarus/utils/functionsanitychecks.hh
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@
namespace Ikarus::utils {
namespace Impl {
/**
*@internal
* @brief Draws the function and returns the slope for a given function.
* @param functionName The name of the function.
* @param ftfunc The function to be evaluated.
Expand Down
14 changes: 14 additions & 0 deletions ikarus/utils/linearalgebrahelper.hh
Original file line number Diff line number Diff line change
Expand Up @@ -448,6 +448,20 @@ namespace Ikarus {
return vec;
}

/**
* \brief Create skew 3x3 matrix from 3d vector.
* \ingroup utils
* \tparam ScalarType The type of the coordinates in the vector.
* \param a The vector.
* \return The skew matrix.
*/
template <typename ScalarType>
Eigen::Matrix<ScalarType, 3, 3> skew(const Eigen::Vector<ScalarType, 3>& a) {
Eigen::Matrix<ScalarType, 3, 3> A;
A << 0, -a(2), a(1), a(2), 0, -a(0), -a(1), a(0), 0;
return A;
}

namespace Impl {
constexpr std::tuple<std::array<std::array<int, 2>, 1>, std::array<std::array<int, 2>, 3>,
std::array<std::array<int, 2>, 6>>
Expand Down
8 changes: 1 addition & 7 deletions python/ikarus/finite_elements/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -104,13 +104,7 @@ def KirchhoffLoveShell(
includes = ["ikarus/python/finiteelements/kirchhoffloveshell.hh"]

includes += ["ikarus/finiteelements/febases/autodifffe.hh"]
autodiffWrapper = "AutoDiffFE"
element_type = (
"Ikarus::"
+ autodiffWrapper
+ f"<Ikarus::KirchhoffLoveShell<{basis.cppTypeName},Ikarus::FERequirements<Eigen::Ref<Eigen::VectorXd>>,true>,"
f"Ikarus::FERequirements<Eigen::Ref<Eigen::VectorXd>>,true>"
)
element_type = f"Ikarus::KirchhoffLoveShell<{basis.cppTypeName},Ikarus::FERequirements<Eigen::Ref<Eigen::VectorXd>>,true>"

# else:
# element_type = "Ikarus::" + func.__name__ + f"<{basis.cppTypeName}, {material.cppTypeName} ,Ikarus::FERequirements<Eigen::Ref<Eigen::VectorXd>>,true>"
Expand Down
Loading

0 comments on commit 2dfc766

Please sign in to comment.