diff --git a/.github/workflows/docsDryRun.yml b/.github/workflows/docsDryRun.yml new file mode 100644 index 000000000..dee844aa1 --- /dev/null +++ b/.github/workflows/docsDryRun.yml @@ -0,0 +1,80 @@ +# SPDX-FileCopyrightText: 2021-2024 The Ikarus Developers mueller@ibb.uni-stuttgart.de +# 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 diff --git a/codespellignore b/codespellignore index e9f217045..30b1d928d 100644 --- a/codespellignore +++ b/codespellignore @@ -1,4 +1,4 @@ # SPDX-FileCopyrightText: 2021-2024 The Ikarus Developers mueller@ibb.uni-stuttgart.de # SPDX-License-Identifier: CC0-1.0 - -tE +te +nd diff --git a/docs/website/doxygen/Doxylocal b/docs/website/doxygen/Doxylocal index 930733569..5973a52b9 100644 --- a/docs/website/doxygen/Doxylocal +++ b/docs/website/doxygen/Doxylocal @@ -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 diff --git a/ikarus/finiteelements/mechanics/CMakeLists.txt b/ikarus/finiteelements/mechanics/CMakeLists.txt index 1c87363bc..466295f21 100644 --- a/ikarus/finiteelements/mechanics/CMakeLists.txt +++ b/ikarus/finiteelements/mechanics/CMakeLists.txt @@ -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 ) diff --git a/ikarus/finiteelements/mechanics/kirchhoffloveshell.hh b/ikarus/finiteelements/mechanics/kirchhoffloveshell.hh index 9fa7e8638..d14cb040b 100644 --- a/ikarus/finiteelements/mechanics/kirchhoffloveshell.hh +++ b/ikarus/finiteelements/mechanics/kirchhoffloveshell.hh @@ -17,38 +17,11 @@ #include #include #include +#include #include namespace Ikarus { - /** - * @brief Helper function to calculate the energy for Kirchhoff-Love shell elements. - * - * This function calculates the energy for Kirchhoff-Love shell elements based on given strain and material - * properties. - * - * - * @tparam ScalarType The scalar type used for calculations. - * @param epsV The Green-Lagrange strains. - * @param Aconv Transformation matrix for strains. - * @param E Young's modulus. - * @param nu Poisson's ratio. - * @return The calculated energy. - */ - template - ScalarType energyHelper(const Eigen::Vector& epsV, const auto& Aconv, double E, double nu) { - const double lambda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu)); - const double mu = E / (2.0 * (1.0 + nu)); - const double lambdbar = 2.0 * lambda * mu / (lambda + 2.0 * mu); - Eigen::TensorFixedSize> moduli; - const auto AconvT = tensorView(Aconv, std::array({3, 3})); - moduli = lambdbar * dyadic(AconvT, AconvT).eval() + 2 * mu * symmetricFourthOrder(Aconv, Aconv); - - auto C = toVoigt(moduli); - auto C33 = C({0, 1, 5}, {0, 1, 5}).eval(); - return 0.5 * epsV.dot(C33 * epsV); - } - /** * @brief Kirchhoff-Love shell finite element class. * @@ -63,18 +36,42 @@ namespace Ikarus { template , bool useEigenRef = false> class KirchhoffLoveShell : public PowerBasisFE { public: - using Basis = Basis_; - using FlatBasis = typename Basis::FlatBasis; - using BasePowerFE = PowerBasisFE; // Handles globalIndices function - using FERequirementType = FERequirements_; - using ResultRequirementsType = ResultRequirements; - using LocalView = typename FlatBasis::LocalView; - using Element = typename LocalView::Element; - using Geometry = typename Element::Geometry; - using GridView = typename FlatBasis::GridView; - using Traits = TraitsFromLocalView; - static constexpr int myDim = Traits::mydim; - static constexpr int worlddim = Traits::worlddim; + using Basis = Basis_; + using FlatBasis = typename Basis::FlatBasis; + using BasePowerFE = PowerBasisFE; // Handles globalIndices function + using FERequirementType = FERequirements_; + using ResultRequirementsType = ResultRequirements; + using LocalView = typename FlatBasis::LocalView; + using Element = typename LocalView::Element; + using Geometry = typename Element::Geometry; + using GridView = typename FlatBasis::GridView; + using Traits = TraitsFromLocalView; + static constexpr int myDim = Traits::mydim; + static constexpr int worlddim = Traits::worlddim; + 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 + struct KinematicVariables { + Eigen::Matrix C; ///< material tangent + Eigen::Vector3 epsV; ///< membrane strain in Voigt notation + Eigen::Vector3 kappaV; ///< bending strain in Voigt notation + Eigen::Matrix j; ///< Jacobian of the deformed geometry + Eigen::Matrix J; ///< Jacobian of the reference geometry + Eigen::Matrix3 h; ///< Hessian of the deformed geometry + Eigen::Matrix3 H; ///< Hessian of the reference geometry + Eigen::Vector3 a3N; ///< Normal vector of the deformed geometry + Eigen::Vector3 a3; ///< normalized normal vector of the deformed geometry + }; /** * @brief Constructor for the KirchhoffLoveShell class. @@ -136,14 +133,16 @@ namespace Ikarus { * @tparam ScalarType The scalar type used for calculations. * @param par The FERequirements. * @param dx Optional additional displacement vector. - * @return A pair containing the displacement function and nodal displacements. + * @return The displacement function. */ - template + template auto getDisplacementFunction(const FERequirementType& par, const std::optional>& dx = std::nullopt) const { - const auto& d = par.getGlobalSolution(Ikarus::FESolutions::displacement); - + const auto& d = par.getGlobalSolution(FESolutions::displacement); + auto geo = std::make_shared::Entity::Geometry>( + this->localView().element().geometry()); Dune::BlockVector> disp(dispAtNodes.size()); + if (dx) for (auto i = 0U; i < disp.size(); ++i) for (auto k2 = 0U; k2 < worlddim; ++k2) @@ -154,10 +153,9 @@ namespace Ikarus { for (auto k2 = 0U; k2 < worlddim; ++k2) disp[i][k2] = d[this->localView().index(this->localView().tree().child(k2).localIndex(i))[0]]; - auto geo = std::make_shared::Entity::Geometry>( - this->localView().element().geometry()); Dune::StandardLocalFunction uFunction(localBasis, disp, geo); - return std::make_pair(uFunction, disp); + + return uFunction; } /** @@ -165,10 +163,31 @@ namespace Ikarus { * * Calculates the scalar value based on the given FERequirements. * - * @param par The FERequirements. + * @param req The FERequirements. * @return The calculated scalar value. */ - double calculateScalar(const FERequirementType& par) const { return calculateScalarImpl(par); } + double calculateScalar(const FERequirementType& req) const { return calculateScalarImpl(req); } + + /** + * @brief Calculate the vector associated with the given FERequirementType. + * + * @tparam ScalarType The scalar type for the calculation. + * @param req The FERequirementType object specifying the requirements for the calculation. + * @param force The vector to store the calculated result. + */ + void calculateVector(const FERequirementType& req, typename Traits::template VectorType<> force) const { + calculateVectorImpl(req, force); + } + /** + * @brief Calculate the matrix associated with the given FERequirementType. + * + * @tparam ScalarType The scalar type for the calculation. + * @param req The FERequirementType object specifying the requirements for the calculation. + * @param K The matrix to store the calculated result. + */ + void calculateMatrix(const FERequirementType& req, typename Traits::template MatrixType<> K) const { + calculateMatrixImpl(req, K); + } /** * @brief Calculate results at local coordinates. @@ -189,14 +208,15 @@ namespace Ikarus { Dune::CachedLocalBasis< std::remove_cvref_t().tree().child(0).finiteElement().localBasis())>> localBasis; - std::function(const Eigen::Vector&, + std::function(const Dune::FieldVector&, const double&)> volumeLoad; - std::function(const Eigen::Vector&, + std::function(const Dune::FieldVector&, const double&)> neumannBoundaryLoad; const BoundaryPatch* neumannBoundary; mutable Dune::BlockVector> dispAtNodes; + DefaultMembraneStrain membraneStrain; double emod_; double nu_; double thickness_; @@ -205,76 +225,211 @@ namespace Ikarus { protected: /** - * @brief Implementation to calculate the scalar value. + * \brief Compute material properties and strains at a given integration point. * - * Implementation to calculate the scalar value based on the given FERequirements and optional additional - * displacement. + * \param gpPos The position of the integration point. + * \param gpIndex The index of the integration point. + * \param geo The geometry object providing position and derivatives. + * \param uFunction The function representing the displacement field. * - * @tparam ScalarType The scalar type used for calculations. - * @param par The FERequirements. - * @param dx Optional additional displacement vector. - * @return The calculated scalar value. + * \tparam gpPos The type of the integration point position. + * \tparam gpIndex The type of the integration point index. + * \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, 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& gpPos, int gpIndex, const Geometry& geo, + const auto& uFunction) const { + using ScalarType = typename std::remove_cvref_t::ctype; + + KinematicVariables kin; + using namespace Dune; + using namespace Dune::DerivativeDirections; + const auto [X, Jd, Hd] = geo.impl().zeroFirstAndSecondDerivativeOfPosition(gpPos); + kin.J = toEigen(Jd); + kin.H = toEigen(Hd); + const Eigen::Matrix A = kin.J * kin.J.transpose(); + Eigen::Matrix G = Eigen::Matrix::Zero(); + + G.block<2, 2>(0, 0) = A; + G(2, 2) = 1; + const Eigen::Matrix GInv = G.inverse(); + 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() * uasMatrix; + kin.h = kin.H + hessianu; + const Eigen::Matrix 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 bV = kin.h * kin.a3; + bV(2) *= 2; // Voigt notation requires the two here + const auto BV = toVoigt(toEigen(geo.impl().secondFundamentalForm(gpPos))); + kin.kappaV = BV - bV; + return kin; + } + + template + void calculateMatrixImpl(const FERequirementType& par, typename Traits::template MatrixType K, + const std::optional>& dx = std::nullopt) const { + K.setZero(); + using namespace Dune::DerivativeDirections; + using namespace Dune; + const auto uFunction = getDisplacementFunction(par, dx); + const auto& lambda = par.getParameter(FEParameter::loadfactor); + const auto geo = this->localView().element().geometry(); + + for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) { + const auto intElement = geo.integrationElement(gp.position()) * gp.weight(); + const auto [C, epsV, kappaV, jE, J, h, H, a3N, a3] + = computeMaterialAndStrains(gp.position(), gpIndex, geo, uFunction); + const Eigen::Vector membraneForces = thickness_ * C * epsV; + const Eigen::Vector moments = Dune::power(thickness_, 3) / 12.0 * C * kappaV; + + const auto& Nd = localBasis.evaluateJacobian(gpIndex); + const auto& Ndd = localBasis.evaluateSecondDerivatives(gpIndex); + for (size_t i = 0; i < numberOfNodes; ++i) { + Eigen::Matrix bopIMembrane + = membraneStrain.derivative(gp.position(), jE, Nd, geo, uFunction, localBasis, i); + + Eigen::Matrix bopIBending = bopBending(jE, h, Nd, Ndd, i, a3N, a3); + for (size_t j = i; j < numberOfNodes; ++j) { + auto KBlock = K.template block(worlddim * i, worlddim * j); + Eigen::Matrix bopJMembrane + = membraneStrain.derivative(gp.position(), jE, Nd, geo, uFunction, localBasis, j); + Eigen::Matrix bopJBending = bopBending(jE, h, Nd, Ndd, j, a3N, a3); + KBlock += thickness_ * bopIMembrane.transpose() * C * bopJMembrane * intElement; + KBlock += Dune::power(thickness_, 3) / 12.0 * bopIBending.transpose() * C * bopJBending * intElement; + + Eigen::Matrix kgMembraneIJ + = membraneStrain.secondDerivative(gp.position(), Nd, geo, uFunction, localBasis, membraneForces, i, j); + Eigen::Matrix kgBendingIJ + = kgBending(jE, h, Nd, Ndd, a3N, a3, moments, i, j); + KBlock += kgMembraneIJ * intElement; + KBlock += kgBendingIJ * intElement; + } + } + } + K.template triangularView() = K.transpose(); + } + + template + void calculateVectorImpl(const FERequirementType& par, typename Traits::template VectorType force, + const std::optional>& dx = std::nullopt) const { + force.setZero(); + using namespace Dune::DerivativeDirections; + using namespace Dune; + const auto uFunction = getDisplacementFunction(par, dx); + const auto& lambda = par.getParameter(FEParameter::loadfactor); + const auto geo = this->localView().element().geometry(); + + // Internal forces + for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) { + const auto [C, epsV, kappaV, jE, J, h, H, a3N, a3] + = computeMaterialAndStrains(gp.position(), gpIndex, geo, uFunction); + const Eigen::Vector membraneForces = thickness_ * C * epsV; + const Eigen::Vector moments = Dune::power(thickness_, 3) / 12.0 * C * kappaV; + + const auto& Nd = localBasis.evaluateJacobian(gpIndex); + const auto& Ndd = localBasis.evaluateSecondDerivatives(gpIndex); + for (size_t i = 0; i < numberOfNodes; ++i) { + Eigen::Matrix bopIMembrane + = membraneStrain.derivative(gp.position(), jE, Nd, geo, uFunction, localBasis, i); + Eigen::Matrix bopIBending = bopBending(jE, h, Nd, Ndd, i, a3N, a3); + force.template segment<3>(3 * i) + += bopIMembrane.transpose() * membraneForces * geo.integrationElement(gp.position()) * gp.weight(); + force.template segment<3>(3 * i) + += bopIBending.transpose() * moments * geo.integrationElement(gp.position()) * gp.weight(); + } + } + + // External forces volume forces over the domain + if (volumeLoad) { + const auto u = getDisplacementFunction(par, dx); + for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) { + Eigen::Vector fext = volumeLoad(geo.global(gp.position()), lambda); + for (size_t i = 0; i < numberOfNodes; ++i) { + const auto udCi = uFunction.evaluateDerivative(gpIndex, wrt(coeff(i))); + force.template segment(worlddim * i) + -= udCi * fext * geo.integrationElement(gp.position()) * gp.weight(); + } + } + } + + if (not neumannBoundary and not neumannBoundaryLoad) return; + + // External forces, boundary forces, i.e. at the Neumann boundary + for (auto&& intersection : intersections(neumannBoundary->gridView(), this->localView().element())) { + if (not neumannBoundary or not neumannBoundary->contains(intersection)) continue; + + const auto& quadLine = QuadratureRules::rule(intersection.type(), order); + + for (const auto& curQuad : quadLine) { + // Local position of the quadrature point + const FieldVector& quadPos + = intersection.geometryInInside().global(curQuad.position()); + + const double intElement = intersection.geometry().integrationElement(curQuad.position()) * curQuad.weight(); + for (size_t i = 0; i < numberOfNodes; ++i) { + const auto udCi = uFunction.evaluateDerivative(quadPos, wrt(coeff(i))); + + // Value of the Neumann data at the current position + auto neumannValue = neumannBoundaryLoad(intersection.geometry().global(curQuad.position()), lambda); + force.template segment(worlddim * i) -= udCi * neumannValue * intElement; + } + } + } + } + template auto calculateScalarImpl(const FERequirementType& par, const std::optional>& dx = std::nullopt) const -> ScalarType { using namespace Dune::DerivativeDirections; using namespace Dune; - const auto [uFunction, uNodes] = getDisplacementFunction(par, dx); - const auto& lambda = par.getParameter(Ikarus::FEParameter::loadfactor); - const auto geo = this->localView().element().geometry(); - ScalarType energy = 0.0; - const auto uasMatrix = Dune::viewAsEigenMatrixAsDynFixed(uNodes); + const auto uFunction = getDisplacementFunction(par, dx); + const auto& lambda = par.getParameter(FEParameter::loadfactor); + const auto geo = this->localView().element().geometry(); + ScalarType energy = 0.0; for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) { - const auto [X, Jd, Hd] = geo.impl().zeroFirstAndSecondDerivativeOfPosition(gp.position()); - const auto J = toEigen(Jd); - const auto H = toEigen(Hd); - const Eigen::Matrix A = J * J.transpose(); - const Eigen::Matrix gradu = toEigen( - uFunction.evaluateDerivative(gpIndex, wrt(spatialAll), Dune::on(DerivativeDirections::referenceElement))); - const Eigen::Matrix j = J + gradu.transpose(); - - const auto& Ndd = localBasis.evaluateSecondDerivatives(gpIndex); - const auto h = H + Ndd.transpose().template cast() * uasMatrix; - const Eigen::Vector3 a3 = (j.row(0).cross(j.row(1))).normalized(); - Eigen::Vector bV = h * a3; - bV(2) *= 2; // Voigt notation requires the two here - - Eigen::Matrix G; - G.setZero(); - G.block<2, 2>(0, 0) = A; - G(2, 2) = 1; - const Eigen::Matrix GInv = G.inverse(); - - const auto epsV = toVoigt((0.5 * (j * j.transpose() - A)).eval()).eval(); - const auto BV = toVoigt(toEigen(geo.impl().secondFundamentalForm(gp.position()))); - const auto kappaV = (BV - bV).eval(); - const ScalarType membraneEnergy = thickness_ * energyHelper(epsV, GInv, emod_, nu_); - const ScalarType bendingEnergy = Dune::power(thickness_, 3) / 12.0 * energyHelper(kappaV, GInv, emod_, nu_); + const auto [C, epsV, kappaV, j, J, h, H, a3N, a3] + = computeMaterialAndStrains(gp.position(), gpIndex, geo, uFunction); + + const ScalarType membraneEnergy = 0.5 * thickness_ * epsV.dot(C * epsV); + const ScalarType bendingEnergy = 0.5 * Dune::power(thickness_, 3) / 12.0 * kappaV.dot(C * kappaV); energy += (membraneEnergy + bendingEnergy) * geo.integrationElement(gp.position()) * gp.weight(); } if (volumeLoad) { for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) { const auto u = uFunction.evaluate(gpIndex); - const Eigen::Vector fExt = volumeLoad(toEigen(geo.global(gp.position())), lambda); + const Eigen::Vector fExt = volumeLoad(geo.global(gp.position()), lambda); energy -= u.dot(fExt) * geo.integrationElement(gp.position()) * gp.weight(); } } - // line or surface loads, i.e., neumann boundary if (not neumannBoundary and not neumannBoundaryLoad) return energy; + // line or surface loads, i.e., neumann boundary const auto& element = this->localView().element(); for (auto&& intersection : intersections(neumannBoundary->gridView(), element)) { if (not neumannBoundary or not neumannBoundary->contains(intersection)) continue; - const auto& quadLine = Dune::QuadratureRules::rule(intersection.type(), order); + const auto& quadLine = QuadratureRules::rule(intersection.type(), order); for (const auto& curQuad : quadLine) { // Local position of the quadrature point - const Dune::FieldVector& quadPos + const FieldVector& quadPos = intersection.geometryInInside().global(curQuad.position()); const double intElement = intersection.geometry().integrationElement(curQuad.position()); @@ -283,13 +438,95 @@ 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; } } - return energy; } + + template + Eigen::Matrix kgBending(const Eigen::Matrix& jcur, + const Eigen::Matrix3& h, const auto& dN, const auto& ddN, + const Eigen::Vector3& a3N, + const Eigen::Vector3& a3, const Eigen::Vector3& S, + int I, int J) const { + Eigen::Matrix kg; + kg.setZero(); + + const auto& dN1i = dN(I, 0); + const auto& dN1j = dN(J, 0); + const auto& dN2i = dN(I, 1); + const auto& dN2j = dN(J, 1); + + const Eigen::Matrix P + = 1.0 / a3N.norm() * (Eigen::Matrix::Identity() - a3 * a3.transpose()); + + const auto a1dxI = Eigen::Matrix::Identity() + * dN1i; // the auto here allows the exploitation of the identity matrices, + // due to Eigen's expression templates + const auto a2dxI = Eigen::Matrix::Identity() * dN2i; + const auto a1dxJ = Eigen::Matrix::Identity() * dN1j; + const auto a2dxJ = Eigen::Matrix::Identity() * dN2j; + const auto a1 = jcur.row(0); + const auto a2 = jcur.row(1); + const Eigen::Matrix a3NdI = a1dxI.colwise().cross(a2) - a2dxI.colwise().cross(a1); + const Eigen::Matrix a3NdJ = a1dxJ.colwise().cross(a2) - a2dxJ.colwise().cross(a1); + Eigen::Matrix a3dI = P * a3NdI; + Eigen::Matrix a3dJ = P * a3NdJ; + for (int i = 0; i < 3; ++i) { + const auto a_albe = h.row(i).transpose(); + const auto& ddNI = ddN(I, i); + const auto& ddNJ = ddN(J, i); + Eigen::Vector3 vecd = P * a_albe; + + Eigen::Matrix a3Ndd + = 1.0 / a3N.squaredNorm() + * ((3 * a3 * a3.transpose() - Eigen::Matrix::Identity()) * (a3.dot(a_albe)) + - a_albe * a3.transpose() - a3 * a_albe.transpose()); + + Eigen::Matrix secondDerivativeDirectorIJ = skew(((dN2i * dN1j - dN1i * dN2j) * vecd).eval()); + kg -= (a3NdI.transpose() * a3Ndd * a3NdJ + secondDerivativeDirectorIJ + (ddNI * a3dJ + ddNJ * a3dI.transpose())) + * S[i] * (i == 2 ? 2 : 1); + } + + return kg; + } + + template + Eigen::Matrix bopBending(const Eigen::Matrix& jcur, + const Eigen::Matrix3& h, const auto& dN, const auto& ddN, + const int node, const Eigen::Vector3& a3N, + const Eigen::Vector3& a3) const { + const Eigen::Matrix a1dxI + = Eigen::Matrix::Identity() * dN(node, 0); // this should be double + // but the cross-product below complains otherwise + const Eigen::Matrix a2dxI = Eigen::Matrix::Identity() * dN(node, 1); + const auto a1 = jcur.row(0); + const auto a2 = jcur.row(1); + const Eigen::Matrix a3NdI + = a1dxI.colwise().cross(a2) - a2dxI.colwise().cross(a1); // minus needed since order has + // to be swapped to get column-wise cross product working + const Eigen::Matrix a3d1 + = 1.0 / a3N.norm() * (Eigen::Matrix::Identity() - a3 * a3.transpose()) * a3NdI; + + Eigen::Matrix bop = -(h * a3d1 + (a3 * ddN.row(node)).transpose()); + bop.row(2) *= 2; + + return bop; + } + + Eigen::Matrix materialTangent(const Eigen::Matrix& Aconv) const { + const double lambda = emod_ * nu_ / ((1.0 + nu_) * (1.0 - 2.0 * nu_)); + const double mu = emod_ / (2.0 * (1.0 + nu_)); + const double lambdbar = 2.0 * lambda * mu / (lambda + 2.0 * mu); + Eigen::TensorFixedSize> moduli; + const auto AconvT = tensorView(Aconv, std::array({3, 3})); + moduli = lambdbar * dyadic(AconvT, AconvT).eval() + 2.0 * mu * symmetricFourthOrder(Aconv, Aconv); + + auto C = toVoigt(moduli); + Eigen::Matrix C33 = C({0, 1, 5}, {0, 1, 5}); + return C33; + } }; } // namespace Ikarus diff --git a/ikarus/finiteelements/mechanics/linearelastic.hh b/ikarus/finiteelements/mechanics/linearelastic.hh index 33f0b7338..567ce62c3 100644 --- a/ikarus/finiteelements/mechanics/linearelastic.hh +++ b/ikarus/finiteelements/mechanics/linearelastic.hh @@ -233,10 +233,10 @@ namespace Ikarus { Dune::CachedLocalBasis< std::remove_cvref_t().tree().child(0).finiteElement().localBasis())>> localBasis; - std::function(const Eigen::Vector&, + std::function(const Dune::FieldVector&, const double&)> volumeLoad; - std::function(const Eigen::Vector&, + std::function(const Dune::FieldVector&, const double&)> neumannBoundaryLoad; const BoundaryPatch* neumannBoundary; @@ -270,7 +270,7 @@ namespace Ikarus { if (volumeLoad) { for (const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) { const auto uVal = u.evaluate(gpIndex); - Eigen::Vector fext = volumeLoad(toEigen(geo.global(gp.position())), lambda); + Eigen::Vector fext = volumeLoad(geo.global(gp.position()), lambda); energy -= uVal.dot(fext) * geo.integrationElement(gp.position()) * gp.weight(); } } @@ -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; } @@ -327,7 +327,7 @@ namespace Ikarus { if (volumeLoad) { const auto u = getDisplacementFunction(par, dx); for (const auto& [gpIndex, gp] : u.viewOverIntegrationPoints()) { - Eigen::Vector fext = volumeLoad(toEigen(geo.global(gp.position())), lambda); + Eigen::Vector 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 * i) @@ -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 * i) -= udCi * neumannValue * curQuad.weight() * integrationElement; } } diff --git a/ikarus/finiteelements/mechanics/membranestrains.hh b/ikarus/finiteelements/mechanics/membranestrains.hh new file mode 100644 index 000000000..891ca4e9d --- /dev/null +++ b/ikarus/finiteelements/mechanics/membranestrains.hh @@ -0,0 +1,109 @@ +// SPDX-FileCopyrightText: 2021-2024 The Ikarus Developers mueller@ibb.uni-stuttgart.de +// SPDX-License-Identifier: LGPL-3.0-or-later + +/** + * \file membranestrains.hh + * \brief Implementation of membrane strain for shells + */ + +#pragma once +#include + +#include + +#include +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 + auto value(const Dune::FieldVector& gpPos, const Geometry& geo, const auto& uFunction) const + -> Eigen::Vector3::ctype> { + using ScalarType = typename std::remove_cvref_t::ctype; + Eigen::Vector3 epsV; + const auto J = Dune::toEigen(geo.jacobianTransposed(gpPos)); + using namespace Dune; + using namespace Dune::DerivativeDirections; + const Eigen::Matrix gradu = toEigen( + uFunction.evaluateDerivative(gpPos, // Here the gpIndex could be passed + Dune::wrt(spatialAll), Dune::on(Dune::DerivativeDirections::referenceElement))); + const Eigen::Matrix 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 + auto derivative(const Dune::FieldVector& gpPos, const Eigen::Matrix& jcur, + const auto& dNAtGp, const Geometry& geo, const auto& uFunction, const auto& localBasis, + const int node) const { + Eigen::Matrix 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 + auto secondDerivative(const Dune::FieldVector& gpPos, const auto& dNAtGp, const Geometry& geo, + const auto& uFunction, const auto& localBasis, const Eigen::Vector3& 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 kg = Eigen::Matrix::Identity() * NS; + return kg; + } + }; + +} // namespace Ikarus diff --git a/ikarus/finiteelements/mechanics/nonlinearelastic.hh b/ikarus/finiteelements/mechanics/nonlinearelastic.hh index 250677acb..49b42abc1 100644 --- a/ikarus/finiteelements/mechanics/nonlinearelastic.hh +++ b/ikarus/finiteelements/mechanics/nonlinearelastic.hh @@ -275,10 +275,10 @@ namespace Ikarus { Dune::CachedLocalBasis< std::remove_cvref_t().tree().child(0).finiteElement().localBasis())>> localBasis; - std::function(const Eigen::Vector&, + std::function(const Dune::FieldVector&, const double&)> volumeLoad; - std::function(const Eigen::Vector&, + std::function(const Dune::FieldVector&, const double&)> neumannBoundaryLoad; const BoundaryPatch* neumannBoundary; @@ -308,7 +308,7 @@ namespace Ikarus { if (volumeLoad) { for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) { const auto u = uFunction.evaluate(gpIndex); - const Eigen::Vector fExt = volumeLoad(toEigen(geo.global(gp.position())), lambda); + const Eigen::Vector fExt = volumeLoad(geo.global(gp.position()), lambda); energy -= u.dot(fExt) * geo.integrationElement(gp.position()) * gp.weight(); } } @@ -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; } } @@ -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 fExt = volumeLoad(toEigen(geo.global(gp.position())), lambda); + const Eigen::Vector 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 * i) -= udCi * fExt * intElement; @@ -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 * i) -= udCi * neumannValue * curQuad.weight() * intElement; } } diff --git a/ikarus/python/finiteelements/kirchhoffloveshell.hh b/ikarus/python/finiteelements/kirchhoffloveshell.hh index 25dcda1c2..e81d941e8 100644 --- a/ikarus/python/finiteelements/kirchhoffloveshell.hh +++ b/ikarus/python/finiteelements/kirchhoffloveshell.hh @@ -51,8 +51,8 @@ namespace Ikarus::Python { using Traits = typename KirchhoffLoveShell::Traits; using FERequirements = typename KirchhoffLoveShell::FERequirementType; - using LoadFunction = std::function(Eigen::Vector, - const double&)>; + using LoadFunction = std::function( + Dune::FieldVector, 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); diff --git a/ikarus/python/finiteelements/nonlinearelastic.hh b/ikarus/python/finiteelements/nonlinearelastic.hh index 550014728..db6b58913 100644 --- a/ikarus/python/finiteelements/nonlinearelastic.hh +++ b/ikarus/python/finiteelements/nonlinearelastic.hh @@ -39,8 +39,8 @@ namespace Ikarus::Python { }), pybind11::keep_alive<1, 2>(), pybind11::keep_alive<1, 3>()); - using LoadFunction = std::function(Eigen::Vector, - const double&)>; + using LoadFunction = std::function( + Dune::FieldVector, const double&)>; cls.def(pybind11::init( [](const GlobalBasis& basis, const Element& element, const Material& mat, diff --git a/ikarus/python/finiteelements/registerelement.hh b/ikarus/python/finiteelements/registerelement.hh index b7d1d025d..c9e6cf033 100644 --- a/ikarus/python/finiteelements/registerelement.hh +++ b/ikarus/python/finiteelements/registerelement.hh @@ -59,8 +59,8 @@ namespace Ikarus::Python { }), pybind11::keep_alive<1, 2>(), pybind11::keep_alive<1, 3>()); - using LoadFunction = std::function(Eigen::Vector, - const double&)>; + using LoadFunction = std::function( + Dune::FieldVector, const double&)>; if constexpr (defaultInitializers) cls.def( pybind11::init([](const GlobalBasis& basis, const Element& element, double emod, double nu, diff --git a/ikarus/utils/functionsanitychecks.hh b/ikarus/utils/functionsanitychecks.hh index a348dcb8c..4f7e35888 100644 --- a/ikarus/utils/functionsanitychecks.hh +++ b/ikarus/utils/functionsanitychecks.hh @@ -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. diff --git a/ikarus/utils/linearalgebrahelper.hh b/ikarus/utils/linearalgebrahelper.hh index 9f914e694..1f0ed841b 100644 --- a/ikarus/utils/linearalgebrahelper.hh +++ b/ikarus/utils/linearalgebrahelper.hh @@ -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 + Eigen::Matrix skew(const Eigen::Vector& a) { + Eigen::Matrix A; + A << 0, -a(2), a(1), a(2), 0, -a(0), -a(1), a(0), 0; + return A; + } + namespace Impl { constexpr std::tuple, 1>, std::array, 3>, std::array, 6>> diff --git a/python/ikarus/finite_elements/__init__.py b/python/ikarus/finite_elements/__init__.py index f395d4eac..7072c1dbd 100644 --- a/python/ikarus/finite_elements/__init__.py +++ b/python/ikarus/finite_elements/__init__.py @@ -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">,true>," - f"Ikarus::FERequirements>,true>" - ) + element_type = f"Ikarus::KirchhoffLoveShell<{basis.cppTypeName},Ikarus::FERequirements>,true>" # else: # element_type = "Ikarus::" + func.__name__ + f"<{basis.cppTypeName}, {material.cppTypeName} ,Ikarus::FERequirements>,true>" diff --git a/tests/src/CMakeLists.txt b/tests/src/CMakeLists.txt index ecd914898..c55d20182 100644 --- a/tests/src/CMakeLists.txt +++ b/tests/src/CMakeLists.txt @@ -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) diff --git a/tests/src/checkfebyautodiff.hh b/tests/src/checkfebyautodiff.hh new file mode 100644 index 000000000..230617f42 --- /dev/null +++ b/tests/src/checkfebyautodiff.hh @@ -0,0 +1,69 @@ +// SPDX-FileCopyrightText: 2021-2024 The Ikarus Developers mueller@ibb.uni-stuttgart.de +// SPDX-License-Identifier: LGPL-3.0-or-later + +#pragma once +#include + +#include + +#include +#include + +template