diff --git a/ikarus/finiteelements/mechanics/kirchhoffloveshell.hh b/ikarus/finiteelements/mechanics/kirchhoffloveshell.hh index 15e7373f2..7aa9572bc 100644 --- a/ikarus/finiteelements/mechanics/kirchhoffloveshell.hh +++ b/ikarus/finiteelements/mechanics/kirchhoffloveshell.hh @@ -194,8 +194,9 @@ namespace Ikarus { } Dune::CachedLocalBasis> localBasis; - std::function(const Eigen::Vector&, const double&)> volumeLoad; - std::function(const Eigen::Vector&, const double&)> + std::function(const Dune::FieldVector&, const double&)> + volumeLoad; + std::function(const Dune::FieldVector&, const double&)> neumannBoundaryLoad; const BoundaryPatch* neumannBoundary; DefaultMembraneStrain membraneStrain; @@ -233,7 +234,7 @@ namespace Ikarus { KinematicVariables kin; using namespace Dune; using namespace Dune::DerivativeDirections; - const auto [X, Jd, Hd] = geo.impl().zeroFirstAndSecondDerivativeOfPosition(gpPos); + 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(); @@ -257,7 +258,7 @@ namespace Ikarus { 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))); + const auto BV = toVoigt(toEigen(geo_->impl().secondFundamentalForm(gpPos))); kin.kappaV = BV - bV; return kin; } @@ -268,35 +269,34 @@ namespace Ikarus { K.setZero(); using namespace Dune::DerivativeDirections; using namespace Dune; - const auto uFunction = getDisplacementFunction(par, dx); + const auto uFunction = displacementFunction(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 intElement = geo_->integrationElement(gp.position()) * gp.weight(); const auto [C, epsV, kappaV, jE, J, h, H, a3N, a3] - = computeMaterialAndStrains(gp.position(), gpIndex, geo, uFunction); + = 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 bopIMembrane + = membraneStrain.derivative(gp.position(), jE, Nd, *geo_, uFunction, localBasis, i); - Eigen::Matrix bopIBending = bopBending(jE, h, Nd, Ndd, i, a3N, a3); + 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); + 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 + 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; @@ -312,14 +312,13 @@ namespace Ikarus { force.setZero(); using namespace Dune::DerivativeDirections; using namespace Dune; - const auto uFunction = getDisplacementFunction(par, dx); + const auto uFunction = displacementFunction(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); + = 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; @@ -327,24 +326,24 @@ namespace Ikarus { 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); + = 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(); + += bopIMembrane.transpose() * membraneForces * geo_->integrationElement(gp.position()) * gp.weight(); force.template segment<3>(3 * i) - += bopIBending.transpose() * moments * geo.integrationElement(gp.position()) * gp.weight(); + += bopIBending.transpose() * moments * geo_->integrationElement(gp.position()) * gp.weight(); } } // External forces volume forces over the domain if (volumeLoad) { - const auto u = getDisplacementFunction(par, dx); + const auto u = displacementFunction(par, dx); for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) { - Eigen::Vector fext = volumeLoad(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 = uFunction.evaluateDerivative(gpIndex, wrt(coeff(i))); - force.template segment(worlddim * i) - -= udCi * fext * geo.integrationElement(gp.position()) * gp.weight(); + force.template segment(worldDim * i) + -= udCi * fext * geo_->integrationElement(gp.position()) * gp.weight(); } } } @@ -368,7 +367,7 @@ namespace Ikarus { // 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; + force.template segment(worldDim * i) -= udCi * neumannValue * intElement; } } } @@ -385,17 +384,17 @@ namespace Ikarus { for (const auto& [gpIndex, gp] : uFunction.viewOverIntegrationPoints()) { const auto [C, epsV, kappaV, j, J, h, H, a3N, a3] - = computeMaterialAndStrains(gp.position(), gpIndex, geo, uFunction); + = 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(); + 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 auto u = uFunction.evaluate(gpIndex); + const Eigen::Vector fExt = volumeLoad(geo_->global(gp.position()), lambda); energy -= u.dot(fExt) * geo_->integrationElement(gp.position()) * gp.weight(); } } diff --git a/ikarus/finiteelements/mechanics/linearelastic.hh b/ikarus/finiteelements/mechanics/linearelastic.hh index 01a749405..d79f47b7c 100644 --- a/ikarus/finiteelements/mechanics/linearelastic.hh +++ b/ikarus/finiteelements/mechanics/linearelastic.hh @@ -219,8 +219,9 @@ namespace Ikarus { } Dune::CachedLocalBasis> localBasis; - std::function(const Eigen::Vector&, const double&)> volumeLoad; - std::function(const Eigen::Vector&, const double&)> + std::function(const Dune::FieldVector&, const double&)> + volumeLoad; + std::function(const Dune::FieldVector&, const double&)> neumannBoundaryLoad; const BoundaryPatch* neumannBoundary; double emod_; @@ -252,7 +253,7 @@ namespace Ikarus { if (volumeLoad) { for (const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) { const auto uVal = uFunction.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(); } } @@ -308,7 +309,7 @@ namespace Ikarus { // External forces volume forces over the domain if (volumeLoad) { for (const auto& [gpIndex, gp] : uFunction.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 = uFunction.evaluateDerivative(gpIndex, wrt(coeff(i))); force.template segment(myDim * i) diff --git a/ikarus/finiteelements/mechanics/nonlinearelastic.hh b/ikarus/finiteelements/mechanics/nonlinearelastic.hh index 1acf42831..4b00c61de 100644 --- a/ikarus/finiteelements/mechanics/nonlinearelastic.hh +++ b/ikarus/finiteelements/mechanics/nonlinearelastic.hh @@ -260,8 +260,9 @@ namespace Ikarus { } Dune::CachedLocalBasis> localBasis; - std::function(const Eigen::Vector&, const double&)> volumeLoad; - std::function(const Eigen::Vector&, const double&)> + std::function(const Dune::FieldVector&, const double&)> + volumeLoad; + std::function(const Dune::FieldVector&, const double&)> neumannBoundaryLoad; const BoundaryPatch* neumannBoundary; Material mat; @@ -289,7 +290,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(); } } @@ -345,7 +346,7 @@ namespace Ikarus { if (volumeLoad) { for (const auto& [gpIndex, gp] : uFunction.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 = uFunction.evaluateDerivative(gpIndex, wrt(coeff(i))); force.template segment(myDim * i) -= udCi * fExt * intElement; diff --git a/tests/src/testklshell.cpp b/tests/src/testklshell.cpp index f77384eeb..06376b8fe 100644 --- a/tests/src/testklshell.cpp +++ b/tests/src/testklshell.cpp @@ -208,4 +208,5 @@ int main(int argc, char** argv) { TestSuite t("Kirchhoff-Love"); t.subTest(singleElementTest()); t.subTest(NonLinearKLShellLoadControlTR()); + return t.exit(); }