Skip to content

Commit

Permalink
fixes
Browse files Browse the repository at this point in the history
  • Loading branch information
@rath3t
rath3t committed Jan 9, 2024
1 parent 7de7646 commit 1a4a214
Show file tree
Hide file tree
Showing 7 changed files with 78 additions and 51 deletions.
2 changes: 1 addition & 1 deletion .github/workflows/docsDryRun.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ on:
- '**.md'

jobs:
Build:
Build-Docs-Dry-Run:
runs-on: ubuntu-latest
container:
image: ikarusproject/ikarus-dev:latest
Expand Down
86 changes: 55 additions & 31 deletions ikarus/finiteelements/mechanics/kirchhoffloveshell.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,28 @@ namespace Ikarus {
static constexpr int membraneStrainSize = 3;
static constexpr int bendingStrainSize = 3;

/**
* \brief A structure representing kinematic variables.
*
* This structure holds various kinematic variables used in a mechanical analysis. It includes material tangent,
* membrane strain, bending strain, Jacobian matrices of deformed and reference geometries, Hessian matrices of
* deformed and reference geometries, the normal vector, and the normalized normal vector.
*
* \tparam ScalarType The scalar type for the matrix and vector elements.
*/
template <typename ScalarType = double>
struct KinematicVariables {
Eigen::Matrix<double, 3, 3> C; ///< material tangent
Eigen::Vector3<ScalarType> epsV; ///< membrane strain in Voigt notation
Eigen::Vector3<ScalarType> kappaV; ///< bending strain in Voigt notation
Eigen::Matrix<ScalarType, 2, 3> j; ///< Jacobian of the deformed geometry
Eigen::Matrix<double, 2, 3> J; ///< Jacobian of the reference geometry
Eigen::Matrix3<ScalarType> h; ///< Hessian of the deformed geometry
Eigen::Matrix3<double> H; ///< Hessian of the reference geometry
Eigen::Vector3<ScalarType> a3N; ///< Normal vector of the deformed geometry
Eigen::Vector3<ScalarType> a3; ///< normalized normal vector of the deformed geometry
};

/**
* @brief Constructor for the KirchhoffLoveShell class.
*
Expand Down Expand Up @@ -215,44 +237,46 @@ namespace Ikarus {
* \tparam geo The type of the geometry object.
* \tparam uFunction The type of the displacement field function.
*
* \return A tuple containing the material tangent, membrane strain, curvature variation,
* Jacobian matrix, metric tensor, Hessian matrix, second fundamental form,
* normalized normal vector, and normal vector at the given integration point.
* \return A tuple containing the material tangent, membrane strain, bending,
* Jacobian matrix of the reference position, Jacobian matrix of the current position, Hessian matrix of
* the current position, Hessian matrix of
* the reference position, normal vector, and normalized normal vector at the given
* integration point.
*/
auto computeMaterialAndStrains(const Dune::FieldVector<double, 2> &gpPos,
int gpIndex,
const Geometry &geo,
const auto &uFunction) const {
auto computeMaterialAndStrains(const Dune::FieldVector<double, 2> &gpPos, int gpIndex, const Geometry &geo,
const auto &uFunction) const {
using ScalarType = typename std::remove_cvref_t<decltype(uFunction)>::ctype;

KinematicVariables<ScalarType> kin;
using namespace Dune;
using namespace Dune::DerivativeDirections;
const auto [X, Jd, Hd] = geo.impl().zeroFirstAndSecondDerivativeOfPosition(gpPos);
const auto J = toEigen(Jd);
const auto H = toEigen(Hd);
const Eigen::Matrix<double, 2, 2> A = J*J.transpose();
Eigen::Matrix<double, 3, 3> G;
G.setZero();
G.block<2, 2>(0, 0) = A;
G(2, 2) = 1;
const auto [X, Jd, Hd] = geo.impl().zeroFirstAndSecondDerivativeOfPosition(gpPos);
kin.J = toEigen(Jd);
kin.H = toEigen(Hd);
const Eigen::Matrix<double, 2, 2> A = kin.J * kin.J.transpose();
Eigen::Matrix<double, 3, 3> G = Eigen::Matrix<double, 3, 3>::Zero();

G.block<2, 2>(0, 0) = A;
G(2, 2) = 1;
const Eigen::Matrix<double, 3, 3> GInv = G.inverse();
const auto C = materialTangent(GInv);

Eigen::Vector3<ScalarType> epsV=membraneStrain.value(gpPos,geo,uFunction);

const auto &Ndd = localBasis.evaluateSecondDerivatives(gpIndex);
const auto uasMatrix = Dune::viewAsEigenMatrixAsDynFixed(uFunction.coefficientsRef());
const auto hessianu = Ndd.transpose().template cast<ScalarType>()*uasMatrix;
const Eigen::Matrix3<ScalarType> h = H + hessianu;
const Eigen::Matrix<ScalarType, 3, 2> gradu = toEigen(
uFunction.evaluateDerivative(gpIndex, Dune::wrt(spatialAll), Dune::on(Dune::DerivativeDirections::referenceElement)));
const Eigen::Matrix<ScalarType, 2, 3> j = J + gradu.transpose();
const Eigen::Vector3<ScalarType> a3N = (j.row(0).cross(j.row(1)));
const Eigen::Vector3<ScalarType> a3 = a3N.normalized();
Eigen::Vector<ScalarType, 3> bV = h*a3;
kin.C = materialTangent(GInv);

kin.epsV = membraneStrain.value(gpPos, geo, uFunction);

const auto &Ndd = localBasis.evaluateSecondDerivatives(gpIndex);
const auto uasMatrix = Dune::viewAsEigenMatrixAsDynFixed(uFunction.coefficientsRef());
const auto hessianu = Ndd.transpose().template cast<ScalarType>() * uasMatrix;
kin.h = kin.H + hessianu;
const Eigen::Matrix<ScalarType, 3, 2> gradu = toEigen(uFunction.evaluateDerivative(
gpIndex, Dune::wrt(spatialAll), Dune::on(Dune::DerivativeDirections::referenceElement)));
kin.j = kin.J + gradu.transpose();
kin.a3N = (kin.j.row(0).cross(kin.j.row(1)));
kin.a3 = kin.a3N.normalized();
Eigen::Vector<ScalarType, 3> bV = kin.h * kin.a3;
bV(2) *= 2; // Voigt notation requires the two here
const auto BV = toVoigt(toEigen(geo.impl().secondFundamentalForm(gpPos)));
const auto kappaV = (BV - bV).eval();
return std::make_tuple(C, epsV, kappaV, j, J, h,H, a3N, a3);
kin.kappaV = BV - bV;
return kin;
}

template <typename ScalarType>
Expand Down
3 changes: 2 additions & 1 deletion tests/src/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,11 +18,12 @@ set(TEST_DEPENDING_ON_LOCALFEFUNCTIONS
testlinearelasticity
testnonlinearelasticityneohooke
testnonlinearelasticitysvk
testadaptivestepsizing
)

set(TEST_DEPENDING_ON_IGA testklshell testadaptivestepsizing)

set(TEST_NEED_MORE_TIME testadaptivestepsizing)
set(TEST_NEED_MORE_TIME)

foreach(programSourceFile ${programSourceFiles})
get_filename_component(programName ${programSourceFile} NAME_WLE)
Expand Down
18 changes: 11 additions & 7 deletions tests/src/checkfebyautodiff.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,14 +3,16 @@

#pragma once
#include <dune/common/test/testsuite.hh>
#include <ikarus/utils/basis.hh>
#include <ikarus/finiteelements/febases/autodifffe.hh>

#include <Eigen/Core>

#include <ikarus/finiteelements/febases/autodifffe.hh>
#include <ikarus/utils/basis.hh>

template <template <typename...> class FE, typename GridView, typename PreBasis, typename... ElementArgsType>
auto checkFEByAutoDiff(const GridView& gridView, const PreBasis& pb, const ElementArgsType&... eleArgs) {
auto basis = Ikarus::makeBasis(gridView, pb);
auto element = gridView.template begin<0>();
auto basis = Ikarus::makeBasis(gridView, pb);
auto element = gridView.template begin<0>();
auto localView = basis.flat().localView();
localView.bind(*element);
auto nDOF = localView.size();
Expand All @@ -21,7 +23,7 @@ auto checkFEByAutoDiff(const GridView& gridView, const PreBasis& pb, const Eleme
const std::string feClassName = Dune::className(fe);
Dune::TestSuite t("Check calculateScalarImpl() and calculateVectorImpl() by Automatic Differentiation for gridDim = "
+ feClassName);
using AutoDiffBasedFE = Ikarus::AutoDiffFE<decltype(fe),typename decltype(fe)::FERequirementType,false,true>;
using AutoDiffBasedFE = Ikarus::AutoDiffFE<decltype(fe), typename decltype(fe)::FERequirementType, false, true>;
AutoDiffBasedFE feAutoDiff(fe);

Eigen::VectorXd d;
Expand Down Expand Up @@ -49,12 +51,14 @@ auto checkFEByAutoDiff(const GridView& gridView, const PreBasis& pb, const Eleme
t.check(K.isApprox(KAutoDiff, tol),
"Mismatch between the stiffness matrices obtained from explicit implementation and the one based on "
"automatic differentiation for "
+ feClassName)<<"The difference is "<<(K-KAutoDiff);
+ feClassName)
<< "The difference is " << (K - KAutoDiff);

t.check(R.isApprox(RAutoDiff, tol),
"Mismatch between the residual vectors obtained from explicit implementation and the one based on "
"automatic differentiation for "
+ feClassName)<<"The difference is "<<(R-RAutoDiff);
+ feClassName)
<< "The difference is " << (R - RAutoDiff);

t.check(Dune::FloatCmp::eq(fe.calculateScalar(req), feAutoDiff.calculateScalar(req), tol),
"Mismatch between the energies obtained from explicit implementation and the one based on "
Expand Down
2 changes: 1 addition & 1 deletion tests/src/testadaptivestepsizing.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -92,7 +92,7 @@ auto KLShellAndAdaptiveStepSizing(const PathFollowingType& pft, const std::vecto
std::vector<ElementType> fes;

for (auto& element : elements(gridView))
fes.emplace_back(basis, element, E, nu, thickness, {}, &neumannBoundary, neumannBoundaryLoad);
fes.emplace_back(basis, element, E, nu, thickness, utils::LoadDefault{}, &neumannBoundary, neumannBoundaryLoad);

auto basisP = std::make_shared<const decltype(basis)>(basis);
DirichletValues dirichletValues(basisP->flat());
Expand Down
1 change: 0 additions & 1 deletion tests/src/testklshell.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -208,5 +208,4 @@ int main(int argc, char** argv) {
TestSuite t("Kirchhoff-Love");
t.subTest(singleElementTest());
t.subTest(NonLinearKLShellLoadControlTR());

}
17 changes: 8 additions & 9 deletions tests/src/testnonlinearelasticitysvk.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,20 +9,19 @@

using Dune::TestSuite;


template <typename Basis_, typename Material,typename FERequirements_ = Ikarus::FErequirements<>>
struct NonLinearElasticHelper : Ikarus::NonLinearElastic<Basis_, Material,FERequirements_, false> {
using Base = Ikarus::NonLinearElastic<Basis_, Material,FERequirements_, false>;
template <typename Basis_, typename Material, typename FERequirements_ = Ikarus::FErequirements<>>
struct NonLinearElasticHelper : Ikarus::NonLinearElastic<Basis_, Material, FERequirements_, false> {
using Base = Ikarus::NonLinearElastic<Basis_, Material, FERequirements_, false>;
using Base::Base;
using FlatBasis = typename Basis_::FlatBasis;

using LocalView = typename FlatBasis::LocalView;
using GridView = typename FlatBasis::GridView;

template <typename VolumeLoad = Ikarus::utils::LoadDefault, typename NeumannBoundaryLoad = Ikarus::utils::LoadDefault>
NonLinearElasticHelper(const Basis_& globalBasis, const typename LocalView::Element& element,const Material& mat, VolumeLoad p_volumeLoad = {},
const BoundaryPatch<GridView>* p_neumannBoundary = nullptr,
NeumannBoundaryLoad p_neumannBoundaryLoad = {})
NonLinearElasticHelper(const Basis_& globalBasis, const typename LocalView::Element& element, const Material& mat,
VolumeLoad p_volumeLoad = {}, const BoundaryPatch<GridView>* p_neumannBoundary = nullptr,
NeumannBoundaryLoad p_neumannBoundaryLoad = {})
: Base(globalBasis, element, mat, p_volumeLoad, p_neumannBoundary, p_neumannBoundaryLoad) {}
};

Expand Down Expand Up @@ -66,7 +65,7 @@ int main(int argc, char** argv) {
Dune::BitSetVector<1> neumannVertices(gridView.size(2), true);
BoundaryPatch neumannBoundary(gridView, neumannVertices);
t.subTest(checkFEByAutoDiff<NonLinearElasticHelper>(gridView, power<2>(lagrange<1>()), reducedMat, volumeLoad,
&neumannBoundary, neumannBoundaryLoad));
&neumannBoundary, neumannBoundaryLoad));
}

{
Expand All @@ -76,7 +75,7 @@ int main(int argc, char** argv) {
Dune::BitSetVector<1> neumannVertices(gridView.size(3), true);
BoundaryPatch neumannBoundary(gridView, neumannVertices);
t.subTest(checkFEByAutoDiff<NonLinearElasticHelper>(gridView, power<3>(lagrange<1>()), matSVK1, volumeLoad,
&neumannBoundary, neumannBoundaryLoad));
&neumannBoundary, neumannBoundaryLoad));
}

return t.exit();
Expand Down

0 comments on commit 1a4a214

Please sign in to comment.