Skip to content

Commit 0dfe755

Browse files
rath3tAlexanderMueller
rath3t
authored and
AlexanderMueller
committed
fix KL shell
1 parent 6078e73 commit 0dfe755

File tree

11 files changed

+49
-44
lines changed

11 files changed

+49
-44
lines changed

‎ikarus/finiteelements/mechanics/linearelastic.hh

Lines changed: 6 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -233,10 +233,10 @@ namespace Ikarus {
233233
Dune::CachedLocalBasis<
234234
std::remove_cvref_t<decltype(std::declval<LocalView>().tree().child(0).finiteElement().localBasis())>>
235235
localBasis;
236-
std::function<Eigen::Vector<double, Traits::worlddim>(const Eigen::Vector<double, Traits::worlddim>&,
236+
std::function<Eigen::Vector<double, Traits::worlddim>(const Dune::FieldVector<double, Traits::worlddim>&,
237237
const double&)>
238238
volumeLoad;
239-
std::function<Eigen::Vector<double, Traits::worlddim>(const Eigen::Vector<double, Traits::worlddim>&,
239+
std::function<Eigen::Vector<double, Traits::worlddim>(const Dune::FieldVector<double, Traits::worlddim>&,
240240
const double&)>
241241
neumannBoundaryLoad;
242242
const BoundaryPatch<GridView>* neumannBoundary;
@@ -270,7 +270,7 @@ namespace Ikarus {
270270
if (volumeLoad) {
271271
for (const auto& [gpIndex, gp] : eps.viewOverIntegrationPoints()) {
272272
const auto uVal = u.evaluate(gpIndex);
273-
Eigen::Vector<double, Traits::worlddim> fext = volumeLoad(toEigen(geo.global(gp.position())), lambda);
273+
Eigen::Vector<double, Traits::worlddim> fext = volumeLoad(geo.global(gp.position()), lambda);
274274
energy -= uVal.dot(fext) * geo.integrationElement(gp.position()) * gp.weight();
275275
}
276276
}
@@ -294,7 +294,7 @@ namespace Ikarus {
294294
const auto uVal = u.evaluate(quadPos);
295295

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

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

359359
// Value of the Neumann data at the current position
360-
auto neumannValue
361-
= neumannBoundaryLoad(toEigen(intersection.geometry().global(curQuad.position())), lambda);
360+
auto neumannValue = neumannBoundaryLoad(intersection.geometry().global(curQuad.position()), lambda);
362361
force.template segment<myDim>(myDim * i) -= udCi * neumannValue * curQuad.weight() * integrationElement;
363362
}
364363
}

‎ikarus/finiteelements/mechanics/membranestrains.hh

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -40,7 +40,8 @@ namespace Ikarus {
4040
const Eigen::Matrix<ScalarType, 2, 3> j = J + gradu.transpose();
4141

4242
epsV << J.row(0).dot(gradu.col(0)) + 0.5 * gradu.col(0).squaredNorm(),
43-
J.row(1).dot(gradu.col(1)) + 0.5 * gradu.col(1).squaredNorm(), j.row(0).dot(j.row(1));
43+
J.row(1).dot(gradu.col(1)) + 0.5 * gradu.col(1).squaredNorm(),
44+
j.row(0).dot(j.row(1)) - J.row(0).dot(J.row(1));
4445
return epsV;
4546
}
4647

‎ikarus/python/finiteelements/nonlinearelastic.hh

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -39,8 +39,8 @@ namespace Ikarus::Python {
3939
}),
4040
pybind11::keep_alive<1, 2>(), pybind11::keep_alive<1, 3>());
4141

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

4545
cls.def(pybind11::init(
4646
[](const GlobalBasis& basis, const Element& element, const Material& mat,

‎ikarus/python/finiteelements/registerelement.hh

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -59,8 +59,8 @@ namespace Ikarus::Python {
5959
}),
6060
pybind11::keep_alive<1, 2>(), pybind11::keep_alive<1, 3>());
6161

62-
using LoadFunction = std::function<Eigen::Vector<double, Traits::worlddim>(Eigen::Vector<double, Traits::worlddim>,
63-
const double&)>;
62+
using LoadFunction = std::function<Eigen::Vector<double, Traits::worlddim>(
63+
Dune::FieldVector<double, Traits::worlddim>, const double&)>;
6464
if constexpr (defaultInitializers)
6565
cls.def(
6666
pybind11::init([](const GlobalBasis& basis, const Element& element, double emod, double nu,

‎python/ikarus/finite_elements/__init__.py

Lines changed: 1 addition & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -104,13 +104,7 @@ def KirchhoffLoveShell(
104104
includes = ["ikarus/python/finiteelements/kirchhoffloveshell.hh"]
105105

106106
includes += ["ikarus/finiteelements/febases/autodifffe.hh"]
107-
autodiffWrapper = "AutoDiffFE"
108-
element_type = (
109-
"Ikarus::"
110-
+ autodiffWrapper
111-
+ f"<Ikarus::KirchhoffLoveShell<{basis.cppTypeName},Ikarus::FERequirements<Eigen::Ref<Eigen::VectorXd>>,true>,"
112-
f"Ikarus::FERequirements<Eigen::Ref<Eigen::VectorXd>>,true>"
113-
)
107+
element_type = f"Ikarus::KirchhoffLoveShell<{basis.cppTypeName},Ikarus::FERequirements<Eigen::Ref<Eigen::VectorXd>>,true>"
114108

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

‎tests/src/checkfebyautodiff.hh

Lines changed: 3 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -11,17 +11,16 @@
1111

1212
template <template <typename...> class FE, typename GridView, typename PreBasis, typename... ElementArgsType>
1313
auto checkFEByAutoDiff(const GridView& gridView, const PreBasis& pb, const ElementArgsType&... eleArgs) {
14-
auto basis = Ikarus::makeBasis(gridView, pb);
14+
auto basis = Ikarus::makeBasis(gridView, pb);
1515
Eigen::VectorXd d;
1616
d.setRandom(basis.flat().dimension());
1717
double lambda = 7.3;
1818

1919
auto req = Ikarus::FErequirements().addAffordance(Ikarus::AffordanceCollections::elastoStatics);
2020
req.insertGlobalSolution(Ikarus::FESolutions::displacement, d)
2121
.insertParameter(Ikarus::FEParameter::loadfactor, lambda);
22-
Dune::TestSuite t("Check calculateScalarImpl() and calculateVectorImpl() by Automatic Differentiation"
23-
);
24-
for(auto element : elements(gridView)) {
22+
Dune::TestSuite t("Check calculateScalarImpl() and calculateVectorImpl() by Automatic Differentiation");
23+
for (auto element : elements(gridView)) {
2524
auto localView = basis.flat().localView();
2625
localView.bind(element);
2726
auto nDOF = localView.size();
@@ -34,8 +33,6 @@ auto checkFEByAutoDiff(const GridView& gridView, const PreBasis& pb, const Eleme
3433
using AutoDiffBasedFE = Ikarus::AutoDiffFE<decltype(fe), typename decltype(fe)::FERequirementType, false, true>;
3534
AutoDiffBasedFE feAutoDiff(fe);
3635

37-
38-
3936
Eigen::MatrixXd K, KAutoDiff;
4037
K.setZero(nDOF, nDOF);
4138
KAutoDiff.setZero(nDOF, nDOF);

‎tests/src/testadaptivestepsizing.cpp

Lines changed: 15 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -3,9 +3,9 @@
33

44
#include <config.h>
55

6+
#include "checkfebyautodiff.hh"
67
#include "testcommon.hh"
78
#include "testhelpers.hh"
8-
#include "checkfebyautodiff.hh"
99

1010
#include <dune/common/test/testsuite.hh>
1111
#include <dune/functions/functionspacebases/boundarydofs.hh>
@@ -31,7 +31,7 @@
3131

3232
using Dune::TestSuite;
3333

34-
template <typename Basis_, typename FERequirements_ = Ikarus::FErequirements<>>
34+
template <typename Basis_, typename FERequirements_ = Ikarus::FERequirements<>>
3535
struct KirchhoffLoveShellHelper : Ikarus::KirchhoffLoveShell<Basis_, FERequirements_, false> {
3636
using Base = Ikarus::KirchhoffLoveShell<Basis_, FERequirements_, false>;
3737
using Base::Base;
@@ -112,8 +112,8 @@ auto KLShellAndAdaptiveStepSizing(const PathFollowingType& pft, const std::vecto
112112
for (auto& element : elements(gridView))
113113
fes.emplace_back(basis, element, E, nu, thickness, utils::LoadDefault{}, &neumannBoundary, neumannBoundaryLoad);
114114

115-
t.subTest(checkFEByAutoDiff<KirchhoffLoveShellHelper>(gridView, power<3>(nurbs()), E, nu, thickness, utils::LoadDefault{},
116-
&neumannBoundary, neumannBoundaryLoad));
115+
t.subTest(checkFEByAutoDiff<KirchhoffLoveShellHelper>(gridView, power<3>(nurbs()), E, nu, thickness,
116+
utils::LoadDefault{}, &neumannBoundary, neumannBoundaryLoad));
117117

118118
auto basisP = std::make_shared<const decltype(basis)>(basis);
119119
DirichletValues dirichletValues(basisP->flat());
@@ -166,7 +166,6 @@ auto KLShellAndAdaptiveStepSizing(const PathFollowingType& pft, const std::vecto
166166
auto nonLinOpFull
167167
= Ikarus::NonLinearOperator(functions(energyFunction, residualFunction, KFunction), parameter(d, lambda));
168168

169-
170169
auto nonLinOp = nonLinOpFull.template subOperator<1, 2>();
171170
auto linSolver = LinearSolver(SolverTypeTag::sd_SimplicialLDLT);
172171

@@ -219,13 +218,16 @@ auto KLShellAndAdaptiveStepSizing(const PathFollowingType& pft, const std::vecto
219218

220219
resetNonLinearOperatorParametersToZero(nonLinOp);
221220
const auto controlInfoWSS = crWSS.run();
222-
checkScalars(t, std::ranges::max(d), expectedResults[0][0], message1 + "<Max Displacement>");
223-
checkScalars(t, lambda, expectedResults[0][1], message1 + "<Lambda>");
221+
const double tolDisp = 1e-13;
222+
const double tolLoad = 1e-12;
223+
checkScalars(t, std::ranges::max(d), expectedResults[0][0], message1 + " <Max Displacement>", tolDisp);
224+
checkScalars(t, lambda, expectedResults[0][1], message1 + " <Lambda>", tolLoad);
224225
resetNonLinearOperatorParametersToZero(nonLinOp);
225226

226227
const auto controlInfoWoSS = crWoSS.run();
227-
checkScalars(t, std::ranges::max(d), expectedResults[1][0], message2 + "<Max Displacement>");
228-
checkScalars(t, lambda, expectedResults[1][1], message2 + "<Lambda>");
228+
229+
checkScalars(t, std::ranges::max(d), expectedResults[1][0], message2 + " <Max Displacement>", tolDisp);
230+
checkScalars(t, lambda, expectedResults[1][1], message2 + " <Lambda>", tolLoad);
229231
resetNonLinearOperatorParametersToZero(nonLinOp);
230232

231233
const int controlInfoWSSIterations
@@ -246,6 +248,9 @@ auto KLShellAndAdaptiveStepSizing(const PathFollowingType& pft, const std::vecto
246248
checkSolverInfos(t, expectedIterations[0], controlInfoWSS, loadSteps, message1);
247249
checkSolverInfos(t, expectedIterations[1], controlInfoWoSS, loadSteps, message2);
248250

251+
t.check(utils::checkGradient(nonLinOpFull, {.draw = false})) << "Check gradient failed";
252+
t.check(utils::checkHessian(nonLinOpFull, {.draw = false})) << "Check hessian failed";
253+
249254
return t;
250255
}
251256

@@ -264,7 +269,7 @@ int main(int argc, char** argv) {
264269
const std::vector<std::vector<double>> expectedResultsALC
265270
= {{0.1032139637288574, 0.0003103004514250302}, {0.162759603260405, 0.0007765975850229621}};
266271
const std::vector<std::vector<double>> expectedResultsLC
267-
= {{0.08741028329554587, 0.0002318693543601816}, {0.144353999993308, 6e-4}};
272+
= {{0.08741028329554552, 0.0002318693543601816}, {0.144353999993308, 6e-4}};
268273

269274
t.subTest(KLShellAndAdaptiveStepSizing(alc, expectedIterationsALC, expectedResultsALC, 3, 0.4));
270275
t.subTest(KLShellAndAdaptiveStepSizing(lc, expectedIterationsLC, expectedResultsLC, 2, 1e-4));

‎tests/src/testfeelement.hh

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -43,15 +43,15 @@ auto testFEElement(const PreBasis& preBasis, const std::string& elementName, con
4343
auto localView = flatBasis.localView();
4444

4545
auto volumeLoad = []<typename VectorType>([[maybe_unused]] const VectorType& globalCoord, auto& lamb) {
46-
VectorType fext;
46+
Eigen::Vector<typename VectorType::field_type, VectorType::dimension> fext;
4747
fext.setZero();
4848
fext[1] = 2 * lamb;
4949
fext[0] = lamb;
5050
return fext;
5151
};
5252

5353
auto neumannBoundaryLoad = []<typename VectorType>([[maybe_unused]] const VectorType& globalCoord, auto& lamb) {
54-
VectorType fext;
54+
Eigen::Vector<typename VectorType::field_type, VectorType::dimension> fext;
5555
fext.setZero();
5656
fext[1] = lamb / 40;
5757
fext[0] = 0;

‎tests/src/testhelpers.hh

Lines changed: 14 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -31,14 +31,23 @@ requires(std::convertible_to<Derived, Eigen::EigenBase<Derived> const&>and std::
3131
}
3232

3333
template <typename TestSuiteType, typename ScalarType>
34+
requires std::is_integral_v<ScalarType>
3435
void checkScalars(TestSuiteType& t, const ScalarType val, const ScalarType expectedVal,
3536
const std::string& messageIfFailed = "") {
3637
if constexpr (std::is_integral_v<ScalarType>)
37-
t.check(val == expectedVal) << std::setprecision(16) << "Incorrect Scalars:\t" << expectedVal << " Actual:\t" << val
38-
<< messageIfFailed;
39-
else
40-
t.check(Dune::FloatCmp::eq(val, expectedVal))
41-
<< std::setprecision(16) << "Incorrect Scalars:\t" << expectedVal << " Actual:\t" << val << messageIfFailed;
38+
t.check(val == expectedVal) << std::setprecision(16) << "Incorrect Scalar integer:\t" << expectedVal << " Actual:\t"
39+
<< val << messageIfFailed;
40+
}
41+
42+
template <typename TestSuiteType, typename ScalarType>
43+
requires(not std::is_integral_v<ScalarType>) void checkScalars(TestSuiteType& t, const ScalarType val,
44+
const ScalarType expectedVal,
45+
const std::string& messageIfFailed = "",
46+
double tol
47+
= Dune::FloatCmp::DefaultEpsilon<ScalarType>::value()) {
48+
t.check(Dune::FloatCmp::eq(val, expectedVal, tol))
49+
<< std::setprecision(16) << "Incorrect Scalar floating point:\t" << expectedVal << " Actual:\t" << val
50+
<< ". The used tolerance was " << tol << messageIfFailed;
4251
}
4352

4453
template <typename TestSuiteType, typename ControlInformation>

‎tests/src/testklshell.cpp

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -153,7 +153,7 @@ static auto NonLinearKLShellLoadControlTR() {
153153
return t;
154154
}
155155

156-
template <typename Basis_, typename FERequirements_ = Ikarus::FErequirements<>>
156+
template <typename Basis_, typename FERequirements_ = Ikarus::FERequirements<>>
157157
struct KirchhoffLoveShellHelper : Ikarus::KirchhoffLoveShell<Basis_, FERequirements_, false> {
158158
using Base = Ikarus::KirchhoffLoveShell<Basis_, FERequirements_, false>;
159159
using Base::Base;

‎tests/src/testnonlinearelasticitysvk.cpp

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -9,7 +9,7 @@
99

1010
using Dune::TestSuite;
1111

12-
template <typename Basis_, typename Material, typename FERequirements_ = Ikarus::FErequirements<>>
12+
template <typename Basis_, typename Material, typename FERequirements_ = Ikarus::FERequirements<>>
1313
struct NonLinearElasticHelper : Ikarus::NonLinearElastic<Basis_, Material, FERequirements_, false> {
1414
using Base = Ikarus::NonLinearElastic<Basis_, Material, FERequirements_, false>;
1515
using Base::Base;

0 commit comments

Comments
 (0)