From 4bea06d35b9f0d0ce7bf9a063c264e67107d3e63 Mon Sep 17 00:00:00 2001 From: Thomas Oberbichler Date: Tue, 9 Jun 2020 21:39:26 +0200 Subject: [PATCH 01/19] Acquire GIL --- include/eqlib/LambdaConstraint.h | 1 + include/eqlib/LambdaObjective.h | 1 + 2 files changed, 2 insertions(+) diff --git a/include/eqlib/LambdaConstraint.h b/include/eqlib/LambdaConstraint.h index 29829285..56e245f7 100644 --- a/include/eqlib/LambdaConstraint.h +++ b/include/eqlib/LambdaConstraint.h @@ -33,6 +33,7 @@ class LambdaConstraint : public Constraint { public: // methods void compute(Ref rs, const std::vector>& gs, const std::vector>& hs) const override { + pybind11::gil_scoped_acquire acquire; m_compute(m_equations, m_variables, rs, gs, hs); } diff --git a/include/eqlib/LambdaObjective.h b/include/eqlib/LambdaObjective.h index e3607255..b196f512 100644 --- a/include/eqlib/LambdaObjective.h +++ b/include/eqlib/LambdaObjective.h @@ -30,6 +30,7 @@ class LambdaObjective : public Objective { public: // methods double compute(Ref g, Ref h) const override { + pybind11::gil_scoped_acquire acquire; return m_compute(m_variables, g, h); } From 90fa824f1aab8ae2d794bc43a9766725b23b7778 Mon Sep 17 00:00:00 2001 From: Thomas Oberbichler Date: Tue, 9 Jun 2020 21:40:29 +0200 Subject: [PATCH 02/19] Use symmetry --- include/eqlib/Problem.h | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/include/eqlib/Problem.h b/include/eqlib/Problem.h index 8079a256..6eab94f5 100644 --- a/include/eqlib/Problem.h +++ b/include/eqlib/Problem.h @@ -438,7 +438,11 @@ class Problem { index index = m_structure_hm.get_index(row.global, col.global); - data.hm_value(index) += h(row.local, col.local); + if (row.local < col.local) { + data.hm_value(index) += h(row.local, col.local); + } else { + data.hm_value(index) += h(col.local, row.local); + } } } From 7854875c5e665e82034096b0dd3f85780da4a09f Mon Sep 17 00:00:00 2001 From: Thomas Oberbichler Date: Tue, 9 Jun 2020 21:43:35 +0200 Subject: [PATCH 03/19] Fix density overflow --- include/eqlib/SparseStructure.h | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/include/eqlib/SparseStructure.h b/include/eqlib/SparseStructure.h index b6be2681..acfd84f8 100644 --- a/include/eqlib/SparseStructure.h +++ b/include/eqlib/SparseStructure.h @@ -78,13 +78,11 @@ class SparseStructure { double density() const noexcept { - const index size = rows() * cols(); - - if (size == 0) { + if (rows() == 0 || cols() == 0) { return 0; } - return (double)nb_nonzeros() / size; + return (double)nb_nonzeros() / rows() / cols(); } const std::vector& ia() const noexcept From 5d036771a8b79c1472974a2c0ea60ffc73d92a29 Mon Sep 17 00:00:00 2001 From: Thomas Oberbichler Date: Tue, 9 Jun 2020 21:49:33 +0200 Subject: [PATCH 04/19] Update hyperjet --- external_libraries/hyperjet/detail/fwd.h | 3 + external_libraries/hyperjet/detail/hyperjet.h | 100 +++++++++++++++++- external_libraries/hyperjet/detail/jet.h | 64 ++++++++++- external_libraries/hyperjet/detail/space.h | 63 +++++++++++ 4 files changed, 225 insertions(+), 5 deletions(-) diff --git a/external_libraries/hyperjet/detail/fwd.h b/external_libraries/hyperjet/detail/fwd.h index 3192380d..44e4bc0a 100644 --- a/external_libraries/hyperjet/detail/fwd.h +++ b/external_libraries/hyperjet/detail/fwd.h @@ -13,6 +13,7 @@ #include #include +#include #include #include @@ -32,6 +33,8 @@ HYPERJET_INLINE index length(const T& container) return static_cast(container.size()); } +constexpr index Dynamic = -1; + constexpr index init_size(const int size) { return size != -1 ? size : 0; diff --git a/external_libraries/hyperjet/detail/hyperjet.h b/external_libraries/hyperjet/detail/hyperjet.h index e17432d1..b018ed0d 100644 --- a/external_libraries/hyperjet/detail/hyperjet.h +++ b/external_libraries/hyperjet/detail/hyperjet.h @@ -256,10 +256,102 @@ class HyperJet { return result; } + auto backward(const std::vector>& xs) const + { + const index size = xs[0].size(); + + auto result = HyperJet::zero(size); + + result.f() = f(); + + backward_to(xs, result.g(), result.h(), true); + + return result; + } + + void backward_to(const std::vector>& xs, Eigen::Ref g, Eigen::Ref h, const bool full) const + { + const index size = g.size(); + + for (index i = 0; i < size; i++) { + for (index r = 0; r < this->size(); r++) { + assert(xs[r].size() == size); + g(i) += this->g(r) * xs[r].g(i); + } + + for (index j = i; j < size; j++) { + for (index r = 0; r < this->size(); r++) { + h(i, j) += this->g(r) * xs[r].h(i, j); + + for (index s = 0; s < this->size(); s++) { + h(i, j) += this->h(r, s) * xs[r].g(i) * xs[s].g(j); + } + } + } + } + + if (!full) { + return; + } + + for (index i = 0; i < size; i++) { + for (index j = 0; j < i; j++) { + h(i, j) = h(j, i); + } + } + } + + auto forward(const std::vector>& xs) const + { + const index size = length(xs); + + auto result = HyperJet::zero(size); + + result.f() = f(); + + forward_to(xs, result.g(), result.h(), true); + + return result; + } + + void forward_to(const std::vector>& xs, Eigen::Ref g, Eigen::Ref h, const bool full) const + { + const index size = g.size(); + + for (index i = 0; i < size; i++) { + g(i) += xs[i].g().dot(this->g()); + + for (index j = i; j < size; j++) { + for (index r = 0; r < this->size(); r++) { + for (index s = 0; s < this->size(); s++) { + h(i, j) += this->g(r) * xs[i].h(r, s) + this->h(r, s) * xs[i].g(r) * xs[j].g(r); + } + } + } + } + + if (!full) { + return; + } + + for (index i = 0; i < size; i++) { + for (index j = 0; j < i; j++) { + h(i, j) = h(j, i); + } + } + } + + std::string to_string() const + { + std::stringstream ss; + ss << *this; + return ss.str(); + } + public: // operators friend std::ostream& operator<<(std::ostream& out, const HyperJet& value) { - out << "HyperJet<" << value.m_f << ">"; + out << value.m_f << "hj"; return out; } @@ -778,7 +870,7 @@ class HyperJet { .def("__len__", &Type::size) .def("__pow__", &Type::pow) .def("__pow__", &Type::pow) - // .def("__repr__", &Type::to_string) + .def("__repr__", &Type::to_string) .def("abs", &Type::abs) .def("acos", &Type::acos) .def("arccos", &Type::acos) @@ -787,8 +879,12 @@ class HyperJet { .def("arctan2", &Type::atan2) .def("asin", &Type::asin) .def("atan", &Type::atan) + .def("backward", &Type::backward, "xs"_a) + .def("backward_to", &Type::backward_to, "xs"_a, "g"_a, "h"_a, "full"_a=true) .def("cos", &Type::cos) .def("enlarge", py::overload_cast(&Type::enlarge, py::const_), "left"_a = 0, "right"_a = 0) + .def("forward", &Type::forward, "xs"_a) + .def("forward_to", &Type::forward_to, "xs"_a, "g"_a, "h"_a, "full"_a=true) .def("sin", &Type::sin) .def("sqrt", &Type::sqrt) .def("tan", &Type::tan) diff --git a/external_libraries/hyperjet/detail/jet.h b/external_libraries/hyperjet/detail/jet.h index b17c3e96..da8bf374 100644 --- a/external_libraries/hyperjet/detail/jet.h +++ b/external_libraries/hyperjet/detail/jet.h @@ -50,7 +50,7 @@ class Jet { HYPERJET_INLINE static Type empty() { assert(TSize != -1); - return Type(Scalar(), Vector(TSize), Matrix(TSize, TSize)); + return Type(Scalar(), Vector(TSize)); } HYPERJET_INLINE static Type empty(const index size) @@ -189,10 +189,64 @@ class Jet { return result; } + auto backward(const std::vector>& xs) const + { + const index size = xs[0].size(); + + auto result = Jet::zero(size); + + result.f() = f(); + + backward_to(xs, result.g()); + + return result; + } + + void backward_to(const std::vector>& xs, Eigen::Ref g) const + { + const index size = g.size(); + + for (index i = 0; i < size; i++) { + for (index r = 0; r < this->size(); r++) { + assert(xs[r].size() == size); + g(i) += this->g(r) * xs[r].g(i); + } + } + } + + auto forward(const std::vector>& xs) const + { + const index size = length(xs); + + auto result = Jet::zero(size); + + result.f() = f(); + + forward_to(xs, result.g()); + + return result; + } + + void forward_to(const std::vector>& xs, Eigen::Ref g) const + { + const index size = g.size(); + + for (index i = 0; i < size; i++) { + g(i) = xs[i].g().dot(this->g()); + } + } + + std::string to_string() const + { + std::stringstream ss; + ss << *this; + return ss.str(); + } + public: // operators friend std::ostream& operator<<(std::ostream& out, const Jet& value) { - out << "Jet<" << value.m_f << ">"; + out << value.m_f << "j"; return out; } @@ -677,7 +731,7 @@ class Jet { .def("__len__", &Type::size) .def("__pow__", &Type::pow) .def("__pow__", &Type::pow) - // .def("__repr__", &Type::to_string) + .def("__repr__", &Type::to_string) .def("abs", &Type::abs) .def("acos", &Type::acos) .def("arccos", &Type::acos) @@ -686,8 +740,12 @@ class Jet { .def("arctan2", &Type::atan2) .def("asin", &Type::asin) .def("atan", &Type::atan) + .def("backward", &Type::backward, "xs"_a) + .def("backward_to", &Type::backward_to, "xs"_a, "g"_a) .def("cos", &Type::cos) .def("enlarge", py::overload_cast(&Type::enlarge, py::const_), "left"_a = 0, "right"_a = 0) + .def("forward", &Type::forward, "xs"_a) + .def("forward_to", &Type::forward_to, "xs"_a, "g"_a) .def("sin", &Type::sin) .def("sqrt", &Type::sqrt) .def("tan", &Type::tan) diff --git a/external_libraries/hyperjet/detail/space.h b/external_libraries/hyperjet/detail/space.h index 4a85f41f..3e488e13 100644 --- a/external_libraries/hyperjet/detail/space.h +++ b/external_libraries/hyperjet/detail/space.h @@ -26,6 +26,13 @@ struct Space<0, TScalar, TSize> { template using Matrix = Eigen::Matrix; + static Scalar empty() + { + static_assert(-1 <= TSize); + + return Scalar(); + } + static Scalar constant(TScalar value) { static_assert(-1 <= TSize); @@ -40,6 +47,20 @@ struct Space<0, TScalar, TSize> { return value; } + template + HYPERJET_INLINE static Eigen::Matrix variables(Eigen::Matrix value) + { + static_assert(-1 <= TSize); + + Eigen::Matrix result; + + for (index i = 0; i < TSize; i++) { + result(i) = variable(TOffset + i, value(i)); + } + + return result; + } + static TScalar f(const Scalar& variable) { return variable; @@ -84,6 +105,13 @@ struct Space<1, TScalar, TSize> { template using Matrix = Eigen::Matrix; + static Scalar empty() + { + static_assert(-1 <= TSize); + + return Scalar::empty(); + } + static Scalar constant(TScalar value) { static_assert(-1 <= TSize); @@ -100,6 +128,20 @@ struct Space<1, TScalar, TSize> { return result; } + template + HYPERJET_INLINE static Eigen::Matrix variables(Eigen::Matrix value) + { + static_assert(-1 <= TSize); + + Eigen::Matrix result; + + for (index i = 0; i < TSize; i++) { + result(i) = variable(TOffset + i, value(i)); + } + + return result; + } + static TScalar f(const Scalar& variable) { return variable.f(); @@ -150,6 +192,13 @@ struct Space<2, TScalar, TSize> { template using Matrix = Eigen::Matrix; + static Scalar empty() + { + static_assert(-1 <= TSize); + + return Scalar::empty(); + } + static Scalar constant(TScalar value) { static_assert(-1 <= TSize); @@ -166,6 +215,20 @@ struct Space<2, TScalar, TSize> { return result; } + template + HYPERJET_INLINE static Eigen::Matrix variables(Eigen::Matrix value) + { + static_assert(-1 <= TSize); + + Eigen::Matrix result; + + for (index i = 0; i < TSize; i++) { + result(i) = variable(TOffset + i, value(i)); + } + + return result; + } + static TScalar f(const Scalar& variable) { return variable.f(); From 60e3ae2c28dd989531f89c9749a5aa8423d1e999 Mon Sep 17 00:00:00 2001 From: Thomas Oberbichler Date: Tue, 9 Jun 2020 21:50:13 +0200 Subject: [PATCH 05/19] Add membrane --- include/eqlib/objectives/IgaMembrane3PAD.h | 179 +++++++++++++++++++++ 1 file changed, 179 insertions(+) create mode 100644 include/eqlib/objectives/IgaMembrane3PAD.h diff --git a/include/eqlib/objectives/IgaMembrane3PAD.h b/include/eqlib/objectives/IgaMembrane3PAD.h new file mode 100644 index 00000000..f3504572 --- /dev/null +++ b/include/eqlib/objectives/IgaMembrane3PAD.h @@ -0,0 +1,179 @@ +#pragma once + +#include "IgaUtilities.h" + +#include "../Node.h" +#include "../Objective.h" + +#include + +namespace eqlib { + +class IgaMembrane3PAD : public Objective { +private: // types + using Type = IgaMembrane3PAD; + + struct Data { + Matrix shape_functions; + Eigen::Matrix ref_a; + Eigen::Matrix3d transformation_matrix; + double weight; + }; + +private: // variables + std::vector> m_nodes; + Eigen::Matrix3d m_dm; + std::vector m_data; + +public: // constructor + IgaMembrane3PAD( + const std::vector>& nodes, + const double thickness, + const double youngs_modulus, + const double poissons_ratio) + : m_nodes(nodes) + { + m_variables.reserve(length(nodes) * 3); + + for (const auto node : nodes) { + m_variables.push_back(node->x()); + m_variables.push_back(node->y()); + m_variables.push_back(node->z()); + } + + m_dm << 1, poissons_ratio, 0, poissons_ratio, 1, 0, 0, 0, (1 - poissons_ratio) / 2; + m_dm *= youngs_modulus * thickness / (1 - std::pow(poissons_ratio, 2)); + } + +public: // methods + index add(const Matrix shape_functions, const double weight) + { + using namespace Eigen; + using namespace eqlib::iga_utilities; + + const Vector3d ref_a1 = evaluate_ref_geometry(m_nodes, shape_functions.row(1)); + const Vector3d ref_a2 = evaluate_ref_geometry(m_nodes, shape_functions.row(2)); + + const Vector3d ref_a(ref_a1.dot(ref_a1), ref_a2.dot(ref_a2), ref_a1.dot(ref_a2)); + + const Vector3d e1 = ref_a1.normalized(); + const Vector3d e2 = (ref_a2 - ref_a2.dot(e1) * e1).normalized(); + + const double det = ref_a(0) * ref_a(1) - ref_a(2) * ref_a(2); + + const Vector3d g_ab_con(ref_a(1) / det, ref_a(0) / det, -ref_a(2) / det); + + const Vector3d g_con1 = g_ab_con[0] * ref_a1 + g_ab_con[2] * ref_a2; + const Vector3d g_con2 = g_ab_con[2] * ref_a1 + g_ab_con[1] * ref_a2; + + const double eg11 = e1.dot(g_con1); + const double eg12 = e1.dot(g_con2); + const double eg21 = e2.dot(g_con1); + const double eg22 = e2.dot(g_con2); + + Matrix3d transformation_matrix; + transformation_matrix << eg11 * eg11, eg12 * eg12, 2 * eg11 * eg12, eg21 * eg21, eg22 * eg22, 2 * eg21 * eg22, 2 * eg11 * eg21, 2 * eg12 * eg22, 2 * (eg11 * eg22 + eg12 * eg21); + + m_data.emplace_back({shape_functions, ref_a, transformation_matrix, weight}); + + return m_data.size() - 1; + } + + template + double compute(Ref g, Ref h) const + { + using namespace eqlib::iga_utilities; + using Space = hyperjet::Space; + + static_assert(0 <= TOrder && TOrder <= 2); + + const index nb_nodes = length(m_nodes); + + Eigen::Matrix locations(nb_nodes, 3); + + for (index i = 0; i < nb_nodes; i++) { + locations.row(i) = m_nodes[i]->act_location(); + } + + double f = 0; + + if constexpr (TOrder > 0) { + g.setZero(); + } + + if constexpr (TOrder > 1) { + h.setZero(); + } + + for (const auto& [shape_functions, ref_a, transformation_matrix, weight] : m_data) { + const auto act_a1 = Space::variables<0, 3>(shape_functions.row(1) * locations); + const auto act_a2 = Space::variables<3, 3>(shape_functions.row(2) * locations); + + const Space::Vector<3> act_a(act_a1.dot(act_a1), act_a2.dot(act_a2), act_a1.dot(act_a2)); + + const auto eps = transformation_matrix * (act_a - ref_a).transpose() * 0.5; + + const auto result = eps.dot(m_dm * eps) * weight; + + const index a = shape_functions.cols() * 3; + + f += Space::f(result); + + if constexpr (TOrder > 0) { + for (index r = 0; r < a; r++) { + const index rd = r % 3; + const index ri = r / 3; + + g(0 + r) += result.g(0 + rd) * shape_functions(1, ri) + result.g(3 + rd) * shape_functions(2, ri); + } + } + + if constexpr (TOrder > 1) { + for (index r = 0; r < a; r++) { + const index rd = r % 3; + const index ri = r / 3; + + for (index s = r; s < a; s++) { + const index sd = s % 3; + const index si = s / 3; + + h(0 + r, 0 + s) += result.h(0 + rd, 0 + sd) * shape_functions(1, ri) * shape_functions(1, si) + + result.h(3 + rd, 0 + sd) * shape_functions(2, ri) * shape_functions(1, si) + + result.h(0 + rd, 3 + sd) * shape_functions(1, ri) * shape_functions(2, si) + + result.h(3 + rd, 3 + sd) * shape_functions(2, ri) * shape_functions(2, si); + } + } + } + } + + return f; + } + + double compute(Ref g, Ref h) const override + { + if (g.size() == 0) { + return compute<0>(g, h); + } else if (h.size() == 0) { + return compute<1>(g, h); + } else { + return compute<2>(g, h); + } + } + +public: // python + template + static void register_python(TModule& m) + { + namespace py = pybind11; + using namespace pybind11::literals; + + using Holder = Pointer; + using Base = Objective; + + py::class_(m, "IgaMembrane3PAD") + .def(py::init>, double, double, double>(), "nodes"_a, "thickness"_a, "youngs_modulus"_a, "poissons_ratio"_a) + .def("add", &Type::add, "shape_functions"_a, "weight"_a); + } +}; // class Point + +} // namespace eqlib \ No newline at end of file From c8ef12d9311579208c91ba5dacc9362f4935f514 Mon Sep 17 00:00:00 2001 From: Thomas Oberbichler Date: Tue, 9 Jun 2020 21:50:30 +0200 Subject: [PATCH 06/19] Format --- include/eqlib/objectives/IgaPointDistance.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/eqlib/objectives/IgaPointDistance.h b/include/eqlib/objectives/IgaPointDistance.h index 88c84d04..f1fcb4fa 100644 --- a/include/eqlib/objectives/IgaPointDistance.h +++ b/include/eqlib/objectives/IgaPointDistance.h @@ -70,7 +70,7 @@ class IgaPointDistance : public Objective { const Vector3D delta = point_b - point_a; - f += delta.squaredNorm() * weight / 2; + f += 0.5 * weight * delta.squaredNorm(); if constexpr (TOrder > 0) { index value_index = 0; From 5c84fcb888e736381c3fd5de15d9b7f8e14c4141 Mon Sep 17 00:00:00 2001 From: Thomas Oberbichler Date: Tue, 9 Jun 2020 21:50:57 +0200 Subject: [PATCH 07/19] Use backward mode --- include/eqlib/objectives/IgaPointDistanceAD.h | 92 +++++++++++++++++-- 1 file changed, 82 insertions(+), 10 deletions(-) diff --git a/include/eqlib/objectives/IgaPointDistanceAD.h b/include/eqlib/objectives/IgaPointDistanceAD.h index ca48d23e..40a19be9 100644 --- a/include/eqlib/objectives/IgaPointDistanceAD.h +++ b/include/eqlib/objectives/IgaPointDistanceAD.h @@ -28,8 +28,8 @@ class IgaPointDistanceAD : public Objective { IgaPointDistanceAD( std::vector> nodes_a, std::vector> nodes_b) - : m_nodes_a(nodes_a) - , m_nodes_b(nodes_b) + : m_nodes_a(nodes_a) + , m_nodes_b(nodes_b) { m_variables.reserve(length(nodes_a) * 3 + length(nodes_b) * 3); @@ -57,24 +57,96 @@ class IgaPointDistanceAD : public Objective { double compute(Ref g, Ref h) const { using namespace eqlib::iga_utilities; - using Space = hyperjet::Space<2, double, -1>; + using Space = hyperjet::Space; static_assert(0 <= TOrder && TOrder <= 2); - auto result = Space::Scalar::zero(nb_variables()); + double f = 0; + + if constexpr (TOrder > 0) { + g.setZero(); + } + + if constexpr (TOrder > 1) { + h.setZero(); + } for (const auto& [shape_functions_a, shape_functions_b, weight] : m_data) { - const auto x_a = evaluate_act_geometry_hj_a(m_nodes_a, shape_functions_a.row(0), nb_variables()); - const auto x_b = evaluate_act_geometry_hj_b(m_nodes_b, shape_functions_b.row(0), nb_variables()); + const auto x_a = Space::variables<0, 3>(evaluate_act_geometry(m_nodes_a, shape_functions_a.row(0))); + const auto x_b = Space::variables<3, 3>(evaluate_act_geometry(m_nodes_b, shape_functions_b.row(0))); const auto delta = x_a - x_b; - result += delta.squaredNorm() * weight; + const auto result = 0.5 * weight * delta.squaredNorm(); + + f += Space::f(result); + + const index a = shape_functions_a.cols() * 3; + const index b = shape_functions_b.cols() * 3; + + if constexpr (TOrder > 0) { + for (index r = 0; r < a; r++) { + const index rd = r % 3; + const index ri = r / 3; + + g(0 + r) = result.g(0 + rd) * shape_functions_a(0, ri); + } + + for (index r = 0; r < b; r++) { + const index rd = r % 3; + const index ri = r / 3; + + g(a + r) = result.g(3 + rd) * shape_functions_b(0, ri); + } + } + + if constexpr (TOrder > 1) { + for (index r = 0; r < a; r++) { + const index rd = r % 3; + const index ri = r / 3; + + for (index s = r; s < a; s++) { + const index sd = s % 3; + const index si = s / 3; + + if (rd != sd) { + continue; + } + + h(0 + r, 0 + s) = result.h(0 + rd, 0 + sd) * shape_functions_a(0, ri) * shape_functions_a(0, si); + } + + for (index s = 0; s < b; s++) { + const index sd = s % 3; + const index si = s / 3; + + if (rd != sd) { + continue; + } + + h(0 + r, a + s) = result.h(0 + rd, 3 + sd) * shape_functions_a(0, ri) * shape_functions_b(0, si); + } + } + + for (index r = 0; r < b; r++) { + const index rd = r % 3; + const index ri = r / 3; + + for (index s = r; s < b; s++) { + const index sd = s % 3; + const index si = s / 3; + + if (rd != sd) { + continue; + } + + h(a + r, a + s) = result.h(3 + rd, 3 + sd) * shape_functions_b(0, ri) * shape_functions_b(0, si); + } + } + } } - g = result.g() / 2; - h = result.h() / 2; - return result.f() / 2; + return f; } double compute(Ref g, Ref h) const override From 0dc7c1987ea6e79f3e6a54c7be53b83648307a94 Mon Sep 17 00:00:00 2001 From: Thomas Oberbichler Date: Tue, 9 Jun 2020 21:51:06 +0200 Subject: [PATCH 08/19] Format --- include/eqlib/objectives/IgaPointLocation.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/eqlib/objectives/IgaPointLocation.h b/include/eqlib/objectives/IgaPointLocation.h index 32748c9e..123a1952 100644 --- a/include/eqlib/objectives/IgaPointLocation.h +++ b/include/eqlib/objectives/IgaPointLocation.h @@ -56,7 +56,7 @@ class IgaPointLocation : public Objective { const Vector3D delta = act_x - target; - f += delta.dot(delta) * weight / 2; + f += 0.5 * weight * delta.dot(delta); if constexpr (TOrder > 0) { for (index i = 0; i < length(m_nodes); i++) { From 165294dc10e5994f2d5b0f8fcac5b1d3c5904d5c Mon Sep 17 00:00:00 2001 From: Thomas Oberbichler Date: Tue, 9 Jun 2020 21:51:24 +0200 Subject: [PATCH 09/19] Use backward mode --- .../eqlib/objectives/IgaRotationCouplingAD.h | 100 +++++++++++++++--- 1 file changed, 86 insertions(+), 14 deletions(-) diff --git a/include/eqlib/objectives/IgaRotationCouplingAD.h b/include/eqlib/objectives/IgaRotationCouplingAD.h index 234fa685..be770ac1 100644 --- a/include/eqlib/objectives/IgaRotationCouplingAD.h +++ b/include/eqlib/objectives/IgaRotationCouplingAD.h @@ -72,21 +72,22 @@ class IgaRotationCouplingAD : public Objective { double compute(Ref g, Ref h) const { using namespace eqlib::iga_utilities; - using Space = hyperjet::Space<2, double, -1>; + using Space = hyperjet::Space; + using std::asin; + using std::pow; static_assert(0 <= TOrder && TOrder <= 2); - auto result = Space::Scalar::zero(nb_variables()); + double f = 0; for (const auto& [shape_functions_a, shape_functions_b, ref_a3_a, ref_a3_b, axis, weight] : m_data) { - const auto a1_a = evaluate_act_geometry_hj_a(m_nodes_a, shape_functions_a.row(1), nb_variables()); - const auto a2_a = evaluate_act_geometry_hj_a(m_nodes_a, shape_functions_a.row(2), nb_variables()); + const auto a1_a = Space::variables<0, 3>(evaluate_act_geometry(m_nodes_a, shape_functions_a.row(1))); + const auto a2_a = Space::variables<3, 3>(evaluate_act_geometry(m_nodes_a, shape_functions_a.row(2))); - const auto a3_a = a1_a.cross(a2_a).normalized(); - - const auto a1_b = evaluate_act_geometry_hj_b(m_nodes_b, shape_functions_b.row(1), nb_variables()); - const auto a2_b = evaluate_act_geometry_hj_b(m_nodes_b, shape_functions_b.row(2), nb_variables()); + const auto a1_b = Space::variables<6, 3>(evaluate_act_geometry(m_nodes_b, shape_functions_b.row(1))); + const auto a2_b = Space::variables<9, 3>(evaluate_act_geometry(m_nodes_b, shape_functions_b.row(2))); + const auto a3_a = a1_a.cross(a2_a).normalized(); const auto a3_b = a1_b.cross(a2_b).normalized(); const auto w_a = a3_a - ref_a3_a.transpose(); @@ -95,17 +96,88 @@ class IgaRotationCouplingAD : public Objective { const auto omega_a = ref_a3_a.cross(w_a); const auto omega_b = ref_a3_b.cross(w_b); - const auto angle_a = omega_a.dot(axis).asin(); - const auto angle_b = omega_b.dot(axis).asin(); + const auto angle_a = asin(omega_a.dot(axis)); + const auto angle_b = asin(omega_b.dot(axis)); const auto angular_difference = angle_a - angle_b; - result += angular_difference.pow(2) * weight; + const auto result = 0.5 * weight * pow(angular_difference, 2); + + const index a = shape_functions_a.cols() * 3; + const index b = shape_functions_b.cols() * 3; + + if constexpr (TOrder > 0) { + for (index r = 0; r < a; r++) { + const index rd = r % 3; + const index ri = r / 3; + + g(0 + r) = result.g(0 + rd) * shape_functions_a(1, ri) + result.g(3 + rd) * shape_functions_a(2, ri); + } + + for (index r = 0; r < b; r++) { + const index rd = r % 3; + const index ri = r / 3; + + g(a + r) = result.g(6 + rd) * shape_functions_b(1, ri) + result.g(9 + rd) * shape_functions_b(2, ri); + } + } + + if constexpr (TOrder > 1) { + for (index r = 0; r < a; r++) { + const index rd = r % 3; + const index ri = r / 3; + + for (index s = 0; s < a; s++) { + const index sd = s % 3; + const index si = s / 3; + + h(0 + r, 0 + s) = result.h(0 + rd, 0 + sd) * shape_functions_a(1, ri) * shape_functions_a(1, si) + + result.h(3 + rd, 0 + sd) * shape_functions_a(2, ri) * shape_functions_a(1, si) + + result.h(0 + rd, 3 + sd) * shape_functions_a(1, ri) * shape_functions_a(2, si) + + result.h(3 + rd, 3 + sd) * shape_functions_a(2, ri) * shape_functions_a(2, si); + } + + for (index s = 0; s < b; s++) { + const index sd = s % 3; + const index si = s / 3; + + h(0 + r, a + s) = result.h(0 + rd, 6 + sd) * shape_functions_a(1, ri) * shape_functions_b(1, si) + + result.h(3 + rd, 6 + sd) * shape_functions_a(2, ri) * shape_functions_b(1, si) + + result.h(0 + rd, 9 + sd) * shape_functions_a(1, ri) * shape_functions_b(2, si) + + result.h(3 + rd, 9 + sd) * shape_functions_a(2, ri) * shape_functions_b(2, si); + } + } + + for (index r = 0; r < b; r++) { + const index rd = r % 3; + const index ri = r / 3; + + for (index s = 0; s < a; s++) { + const index sd = s % 3; + const index si = s / 3; + + h(a + r, 0 + s) = result.h(6 + rd, 0 + sd) * shape_functions_b(1, ri) * shape_functions_a(1, si) + + result.h(9 + rd, 0 + sd) * shape_functions_b(2, ri) * shape_functions_a(1, si) + + result.h(6 + rd, 3 + sd) * shape_functions_b(1, ri) * shape_functions_a(2, si) + + result.h(9 + rd, 3 + sd) * shape_functions_b(2, ri) * shape_functions_a(2, si); + } + + for (index s = 0; s < b; s++) { + const index sd = s % 3; + const index si = s / 3; + + h(a + r, a + s) = result.h(6 + rd, 6 + sd) * shape_functions_b(1, ri) * shape_functions_b(1, si) + + result.h(9 + rd, 6 + sd) * shape_functions_b(2, ri) * shape_functions_b(1, si) + + result.h(6 + rd, 9 + sd) * shape_functions_b(1, ri) * shape_functions_b(2, si) + + result.h(9 + rd, 9 + sd) * shape_functions_b(2, ri) * shape_functions_b(2, si); + } + } + } + + f += Space::f(result); } - g = result.g() / 2; - h = result.h() / 2; - return result.f() / 2; + return f; } double compute(Ref g, Ref h) const override From 247f7167f2c9652c61016db11a869809b7ebbe85 Mon Sep 17 00:00:00 2001 From: Thomas Oberbichler Date: Tue, 9 Jun 2020 21:52:26 +0200 Subject: [PATCH 10/19] Register membrane --- src/Module.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/Module.cpp b/src/Module.cpp index 2f018a4f..aa3ae665 100644 --- a/src/Module.cpp +++ b/src/Module.cpp @@ -25,6 +25,7 @@ #include #include +#include #include #include #include @@ -112,6 +113,9 @@ PYBIND11_MODULE(eqlib, m) // objectives: IgaNormalDistance eqlib::IgaNormalDistanceAD::register_python(m); + // objectives: IgaMembrane3PAD + eqlib::IgaMembrane3PAD::register_python(m); + // objectives: IgaPointDistance eqlib::IgaPointDistance::register_python(m); From 46a6d082dba9324c2bba07cd160c2770fd624173 Mon Sep 17 00:00:00 2001 From: Thomas Oberbichler Date: Tue, 9 Jun 2020 21:52:31 +0200 Subject: [PATCH 11/19] Use backward mode --- .../eqlib/objectives/IgaNormalDistanceAD.h | 86 ++++++++++++++++--- 1 file changed, 72 insertions(+), 14 deletions(-) diff --git a/include/eqlib/objectives/IgaNormalDistanceAD.h b/include/eqlib/objectives/IgaNormalDistanceAD.h index 7d5c9f42..678cacaf 100644 --- a/include/eqlib/objectives/IgaNormalDistanceAD.h +++ b/include/eqlib/objectives/IgaNormalDistanceAD.h @@ -57,33 +57,91 @@ class IgaNormalDistanceAD : public Objective { double compute(Ref g, Ref h) const { using namespace eqlib::iga_utilities; - using Space = hyperjet::Space<2, double, -1>; - using Scalar3hj = Space::Scalar; - using Vector3hj = Space::Vector<3>; + using Space = hyperjet::Space; static_assert(0 <= TOrder && TOrder <= 2); - auto result = Scalar3hj::zero(nb_variables()); + double f = 0; for (const auto& [shape_functions_a, shape_functions_b, weight] : m_data) { - const auto a1_a = evaluate_act_geometry_hj_a(m_nodes_a, shape_functions_a.row(1), nb_variables()); - const auto a2_a = evaluate_act_geometry_hj_a(m_nodes_a, shape_functions_a.row(2), nb_variables()); + const auto a1_a = Space::variables<0, 3>(evaluate_act_geometry(m_nodes_a, shape_functions_a.row(1))); + const auto a2_a = Space::variables<3, 3>(evaluate_act_geometry(m_nodes_a, shape_functions_a.row(2))); + + const auto a1_b = Space::variables<6, 3>(evaluate_act_geometry(m_nodes_b, shape_functions_b.row(1))); + const auto a2_b = Space::variables<9, 3>(evaluate_act_geometry(m_nodes_b, shape_functions_b.row(2))); const auto a3_a = a1_a.cross(a2_a).normalized(); - - const auto a1_b = evaluate_act_geometry_hj_b(m_nodes_b, shape_functions_b.row(1), nb_variables()); - const auto a2_b = evaluate_act_geometry_hj_b(m_nodes_b, shape_functions_b.row(2), nb_variables()); - const auto a3_b = a1_b.cross(a2_b).normalized(); const auto delta = a3_b - a3_a; - result += delta.dot(delta) * weight; + const auto result = 0.5 * weight * delta.dot(delta); + + const index a = shape_functions_a.cols() * 3; + const index b = shape_functions_b.cols() * 3; + + if constexpr(TOrder > 0) { + for (index r = 0; r < a; r++) { + const index rd = r % 3; + const index ri = r / 3; + + g(0 + r) = result.g(0 + rd) * shape_functions_a(1, ri) + result.g(3 + rd) * shape_functions_a(2, ri); + } + + for (index r = 0; r < b; r++) { + const index rd = r % 3; + const index ri = r / 3; + + g(a + r) = result.g(6 + rd) * shape_functions_b(1, ri) + result.g(9 + rd) * shape_functions_b(2, ri); + } + } + + if constexpr(TOrder > 1) { + for (index r = 0; r < a; r++) { + const index rd = r % 3; + const index ri = r / 3; + + for (index s = r; s < a; s++) { + const index sd = s % 3; + const index si = s / 3; + + h(0 + r, 0 + s) = result.h(0 + rd, 0 + sd) * shape_functions_a(1, ri) * shape_functions_a(1, si) + + result.h(3 + rd, 0 + sd) * shape_functions_a(2, ri) * shape_functions_a(1, si) + + result.h(0 + rd, 3 + sd) * shape_functions_a(1, ri) * shape_functions_a(2, si) + + result.h(3 + rd, 3 + sd) * shape_functions_a(2, ri) * shape_functions_a(2, si); + } + + for (index s = 0; s < b; s++) { + const index sd = s % 3; + const index si = s / 3; + + h(0 + r, a + s) = result.h(0 + rd, 6 + sd) * shape_functions_a(1, ri) * shape_functions_b(1, si) + + result.h(3 + rd, 6 + sd) * shape_functions_a(2, ri) * shape_functions_b(1, si) + + result.h(0 + rd, 9 + sd) * shape_functions_a(1, ri) * shape_functions_b(2, si) + + result.h(3 + rd, 9 + sd) * shape_functions_a(2, ri) * shape_functions_b(2, si); + } + } + + for (index r = 0; r < b; r++) { + const index rd = r % 3; + const index ri = r / 3; + + for (index s = r; s < b; s++) { + const index sd = s % 3; + const index si = s / 3; + + h(a + r, a + s) = result.h(6 + rd, 6 + sd) * shape_functions_b(1, ri) * shape_functions_b(1, si) + + result.h(9 + rd, 6 + sd) * shape_functions_b(2, ri) * shape_functions_b(1, si) + + result.h(6 + rd, 9 + sd) * shape_functions_b(1, ri) * shape_functions_b(2, si) + + result.h(9 + rd, 9 + sd) * shape_functions_b(2, ri) * shape_functions_b(2, si); + } + } + } + + f += Space::f(result); } - g = result.g() / 2; - h = result.h() / 2; - return result.f() / 2; + return f; } double compute(Ref g, Ref h) const override From 30a2e7722d83af69b5e1ad8defd90260f3b5fe0a Mon Sep 17 00:00:00 2001 From: Thomas Oberbichler Date: Tue, 9 Jun 2020 21:54:44 +0200 Subject: [PATCH 12/19] Check upper triangle --- tests/test_iga_normal_distance_ad.py | 7 +++---- tests/test_iga_point_distance.py | 3 ++- tests/test_iga_point_distance_ad.py | 3 ++- tests/test_iga_rotation_copling_ad.py | 3 ++- 4 files changed, 9 insertions(+), 7 deletions(-) diff --git a/tests/test_iga_normal_distance_ad.py b/tests/test_iga_normal_distance_ad.py index 0e3e82f7..2dd45c2e 100644 --- a/tests/test_iga_normal_distance_ad.py +++ b/tests/test_iga_normal_distance_ad.py @@ -1,7 +1,6 @@ -import eqlib as eq - import pytest - +import eqlib as eq +import numpy as np from numpy.testing import assert_almost_equal if __name__ == '__main__': @@ -64,4 +63,4 @@ def test_element(element): assert_almost_equal(f, DATA['exp_f']) assert_almost_equal(g, DATA['exp_g']) - assert_almost_equal(h, DATA['exp_h']) + assert_almost_equal(np.triu(h), np.triu(DATA['exp_h'])) diff --git a/tests/test_iga_point_distance.py b/tests/test_iga_point_distance.py index db2c13bb..e7062d16 100644 --- a/tests/test_iga_point_distance.py +++ b/tests/test_iga_point_distance.py @@ -1,5 +1,6 @@ import pytest import eqlib as eq +import numpy as np from numpy.testing import assert_almost_equal if __name__ == '__main__': @@ -52,4 +53,4 @@ def test_element(element): assert_almost_equal(f, DATA['exp_f']) assert_almost_equal(g, DATA['exp_g']) - assert_almost_equal(h, DATA['exp_h']) + assert_almost_equal(np.triu(h), np.triu(DATA['exp_h'])) diff --git a/tests/test_iga_point_distance_ad.py b/tests/test_iga_point_distance_ad.py index a1afd332..880088e5 100644 --- a/tests/test_iga_point_distance_ad.py +++ b/tests/test_iga_point_distance_ad.py @@ -1,5 +1,6 @@ import pytest import eqlib as eq +import numpy as np from numpy.testing import assert_almost_equal if __name__ == '__main__': @@ -52,4 +53,4 @@ def test_element(element): assert_almost_equal(f, DATA['exp_f']) assert_almost_equal(g, DATA['exp_g']) - assert_almost_equal(h, DATA['exp_h']) + assert_almost_equal(np.triu(h), np.triu(DATA['exp_h'])) diff --git a/tests/test_iga_rotation_copling_ad.py b/tests/test_iga_rotation_copling_ad.py index db5675ba..717c72b7 100644 --- a/tests/test_iga_rotation_copling_ad.py +++ b/tests/test_iga_rotation_copling_ad.py @@ -1,5 +1,6 @@ import pytest import eqlib as eq +import numpy as np from numpy.testing import assert_almost_equal if __name__ == '__main__': @@ -64,4 +65,4 @@ def test_element(element): assert_almost_equal(f, DATA['exp_f']) assert_almost_equal(g, DATA['exp_g']) - assert_almost_equal(h, DATA['exp_h']) + assert_almost_equal(np.triu(h), np.triu(DATA['exp_h'])) From 766b3029d75aebb9569ac5f4ba7d0f67065d7367 Mon Sep 17 00:00:00 2001 From: Thomas Oberbichler Date: Tue, 9 Jun 2020 22:06:03 +0200 Subject: [PATCH 13/19] Fix duplicate name --- external_libraries/hyperjet/detail/space.h | 24 +++++++++++----------- 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/external_libraries/hyperjet/detail/space.h b/external_libraries/hyperjet/detail/space.h index 3e488e13..75e24346 100644 --- a/external_libraries/hyperjet/detail/space.h +++ b/external_libraries/hyperjet/detail/space.h @@ -47,14 +47,14 @@ struct Space<0, TScalar, TSize> { return value; } - template - HYPERJET_INLINE static Eigen::Matrix variables(Eigen::Matrix value) + template + HYPERJET_INLINE static Eigen::Matrix variables(Eigen::Matrix value) { static_assert(-1 <= TSize); - Eigen::Matrix result; + Eigen::Matrix result; - for (index i = 0; i < TSize; i++) { + for (index i = 0; i < TDerivedSize; i++) { result(i) = variable(TOffset + i, value(i)); } @@ -128,14 +128,14 @@ struct Space<1, TScalar, TSize> { return result; } - template - HYPERJET_INLINE static Eigen::Matrix variables(Eigen::Matrix value) + template + HYPERJET_INLINE static Eigen::Matrix variables(Eigen::Matrix value) { static_assert(-1 <= TSize); - Eigen::Matrix result; + Eigen::Matrix result; - for (index i = 0; i < TSize; i++) { + for (index i = 0; i < TDerivedSize; i++) { result(i) = variable(TOffset + i, value(i)); } @@ -215,14 +215,14 @@ struct Space<2, TScalar, TSize> { return result; } - template - HYPERJET_INLINE static Eigen::Matrix variables(Eigen::Matrix value) + template + HYPERJET_INLINE static Eigen::Matrix variables(Eigen::Matrix value) { static_assert(-1 <= TSize); - Eigen::Matrix result; + Eigen::Matrix result; - for (index i = 0; i < TSize; i++) { + for (index i = 0; i < TDerivedSize; i++) { result(i) = variable(TOffset + i, value(i)); } From 6d5be8a42d8794fc5314bb3e92efc7710a20aa1f Mon Sep 17 00:00:00 2001 From: Thomas Oberbichler Date: Tue, 9 Jun 2020 22:11:56 +0200 Subject: [PATCH 14/19] Fix template --- include/eqlib/objectives/IgaMembrane3PAD.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/eqlib/objectives/IgaMembrane3PAD.h b/include/eqlib/objectives/IgaMembrane3PAD.h index f3504572..db2faf0d 100644 --- a/include/eqlib/objectives/IgaMembrane3PAD.h +++ b/include/eqlib/objectives/IgaMembrane3PAD.h @@ -109,7 +109,7 @@ class IgaMembrane3PAD : public Objective { const auto act_a1 = Space::variables<0, 3>(shape_functions.row(1) * locations); const auto act_a2 = Space::variables<3, 3>(shape_functions.row(2) * locations); - const Space::Vector<3> act_a(act_a1.dot(act_a1), act_a2.dot(act_a2), act_a1.dot(act_a2)); + const Space::template Vector<3> act_a(act_a1.dot(act_a1), act_a2.dot(act_a2), act_a1.dot(act_a2)); const auto eps = transformation_matrix * (act_a - ref_a).transpose() * 0.5; From d075873cce66dc5601b2d73faf884b91757b3808 Mon Sep 17 00:00:00 2001 From: Thomas Oberbichler Date: Tue, 9 Jun 2020 22:15:27 +0200 Subject: [PATCH 15/19] Fix template --- include/eqlib/objectives/IgaMembrane3PAD.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/eqlib/objectives/IgaMembrane3PAD.h b/include/eqlib/objectives/IgaMembrane3PAD.h index db2faf0d..8dc13cf0 100644 --- a/include/eqlib/objectives/IgaMembrane3PAD.h +++ b/include/eqlib/objectives/IgaMembrane3PAD.h @@ -109,7 +109,7 @@ class IgaMembrane3PAD : public Objective { const auto act_a1 = Space::variables<0, 3>(shape_functions.row(1) * locations); const auto act_a2 = Space::variables<3, 3>(shape_functions.row(2) * locations); - const Space::template Vector<3> act_a(act_a1.dot(act_a1), act_a2.dot(act_a2), act_a1.dot(act_a2)); + const typename Space::Vector<3> act_a(act_a1.dot(act_a1), act_a2.dot(act_a2), act_a1.dot(act_a2)); const auto eps = transformation_matrix * (act_a - ref_a).transpose() * 0.5; From 5ea801d4a700b3f8cfe7f47eab6b082613ac1b28 Mon Sep 17 00:00:00 2001 From: Thomas Oberbichler Date: Tue, 9 Jun 2020 22:16:41 +0200 Subject: [PATCH 16/19] Fix template --- include/eqlib/objectives/IgaMembrane3PAD.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/eqlib/objectives/IgaMembrane3PAD.h b/include/eqlib/objectives/IgaMembrane3PAD.h index 8dc13cf0..758d6642 100644 --- a/include/eqlib/objectives/IgaMembrane3PAD.h +++ b/include/eqlib/objectives/IgaMembrane3PAD.h @@ -106,8 +106,8 @@ class IgaMembrane3PAD : public Objective { } for (const auto& [shape_functions, ref_a, transformation_matrix, weight] : m_data) { - const auto act_a1 = Space::variables<0, 3>(shape_functions.row(1) * locations); - const auto act_a2 = Space::variables<3, 3>(shape_functions.row(2) * locations); + const auto act_a1 = Space::template variables<0, 3>(shape_functions.row(1) * locations); + const auto act_a2 = Space::template variables<3, 3>(shape_functions.row(2) * locations); const typename Space::Vector<3> act_a(act_a1.dot(act_a1), act_a2.dot(act_a2), act_a1.dot(act_a2)); From 02ca9aa651c7bbc2d5bd5bc11a0d0f7c7cb78a37 Mon Sep 17 00:00:00 2001 From: Thomas Oberbichler Date: Tue, 9 Jun 2020 22:19:11 +0200 Subject: [PATCH 17/19] Fix template --- include/eqlib/objectives/IgaMembrane3PAD.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/eqlib/objectives/IgaMembrane3PAD.h b/include/eqlib/objectives/IgaMembrane3PAD.h index 758d6642..eabba8fa 100644 --- a/include/eqlib/objectives/IgaMembrane3PAD.h +++ b/include/eqlib/objectives/IgaMembrane3PAD.h @@ -109,7 +109,7 @@ class IgaMembrane3PAD : public Objective { const auto act_a1 = Space::template variables<0, 3>(shape_functions.row(1) * locations); const auto act_a2 = Space::template variables<3, 3>(shape_functions.row(2) * locations); - const typename Space::Vector<3> act_a(act_a1.dot(act_a1), act_a2.dot(act_a2), act_a1.dot(act_a2)); + const typename Space::template Vector<3> act_a(act_a1.dot(act_a1), act_a2.dot(act_a2), act_a1.dot(act_a2)); const auto eps = transformation_matrix * (act_a - ref_a).transpose() * 0.5; From c815e89a269427900b392f3f7450535a303e2814 Mon Sep 17 00:00:00 2001 From: Thomas Oberbichler Date: Tue, 9 Jun 2020 22:21:48 +0200 Subject: [PATCH 18/19] Fix template --- include/eqlib/objectives/IgaNormalDistanceAD.h | 8 ++++---- include/eqlib/objectives/IgaPointDistanceAD.h | 4 ++-- include/eqlib/objectives/IgaRotationCouplingAD.h | 8 ++++---- 3 files changed, 10 insertions(+), 10 deletions(-) diff --git a/include/eqlib/objectives/IgaNormalDistanceAD.h b/include/eqlib/objectives/IgaNormalDistanceAD.h index 678cacaf..38568a39 100644 --- a/include/eqlib/objectives/IgaNormalDistanceAD.h +++ b/include/eqlib/objectives/IgaNormalDistanceAD.h @@ -64,11 +64,11 @@ class IgaNormalDistanceAD : public Objective { double f = 0; for (const auto& [shape_functions_a, shape_functions_b, weight] : m_data) { - const auto a1_a = Space::variables<0, 3>(evaluate_act_geometry(m_nodes_a, shape_functions_a.row(1))); - const auto a2_a = Space::variables<3, 3>(evaluate_act_geometry(m_nodes_a, shape_functions_a.row(2))); + const auto a1_a = Space::template variables<0, 3>(evaluate_act_geometry(m_nodes_a, shape_functions_a.row(1))); + const auto a2_a = Space::template variables<3, 3>(evaluate_act_geometry(m_nodes_a, shape_functions_a.row(2))); - const auto a1_b = Space::variables<6, 3>(evaluate_act_geometry(m_nodes_b, shape_functions_b.row(1))); - const auto a2_b = Space::variables<9, 3>(evaluate_act_geometry(m_nodes_b, shape_functions_b.row(2))); + const auto a1_b = Space::template variables<6, 3>(evaluate_act_geometry(m_nodes_b, shape_functions_b.row(1))); + const auto a2_b = Space::template variables<9, 3>(evaluate_act_geometry(m_nodes_b, shape_functions_b.row(2))); const auto a3_a = a1_a.cross(a2_a).normalized(); const auto a3_b = a1_b.cross(a2_b).normalized(); diff --git a/include/eqlib/objectives/IgaPointDistanceAD.h b/include/eqlib/objectives/IgaPointDistanceAD.h index 40a19be9..32e1bc6c 100644 --- a/include/eqlib/objectives/IgaPointDistanceAD.h +++ b/include/eqlib/objectives/IgaPointDistanceAD.h @@ -72,8 +72,8 @@ class IgaPointDistanceAD : public Objective { } for (const auto& [shape_functions_a, shape_functions_b, weight] : m_data) { - const auto x_a = Space::variables<0, 3>(evaluate_act_geometry(m_nodes_a, shape_functions_a.row(0))); - const auto x_b = Space::variables<3, 3>(evaluate_act_geometry(m_nodes_b, shape_functions_b.row(0))); + const auto x_a = Space::template variables<0, 3>(evaluate_act_geometry(m_nodes_a, shape_functions_a.row(0))); + const auto x_b = Space::template variables<3, 3>(evaluate_act_geometry(m_nodes_b, shape_functions_b.row(0))); const auto delta = x_a - x_b; diff --git a/include/eqlib/objectives/IgaRotationCouplingAD.h b/include/eqlib/objectives/IgaRotationCouplingAD.h index be770ac1..cac872ed 100644 --- a/include/eqlib/objectives/IgaRotationCouplingAD.h +++ b/include/eqlib/objectives/IgaRotationCouplingAD.h @@ -81,11 +81,11 @@ class IgaRotationCouplingAD : public Objective { double f = 0; for (const auto& [shape_functions_a, shape_functions_b, ref_a3_a, ref_a3_b, axis, weight] : m_data) { - const auto a1_a = Space::variables<0, 3>(evaluate_act_geometry(m_nodes_a, shape_functions_a.row(1))); - const auto a2_a = Space::variables<3, 3>(evaluate_act_geometry(m_nodes_a, shape_functions_a.row(2))); + const auto a1_a = Space::template variables<0, 3>(evaluate_act_geometry(m_nodes_a, shape_functions_a.row(1))); + const auto a2_a = Space::template variables<3, 3>(evaluate_act_geometry(m_nodes_a, shape_functions_a.row(2))); - const auto a1_b = Space::variables<6, 3>(evaluate_act_geometry(m_nodes_b, shape_functions_b.row(1))); - const auto a2_b = Space::variables<9, 3>(evaluate_act_geometry(m_nodes_b, shape_functions_b.row(2))); + const auto a1_b = Space::template variables<6, 3>(evaluate_act_geometry(m_nodes_b, shape_functions_b.row(1))); + const auto a2_b = Space::template variables<9, 3>(evaluate_act_geometry(m_nodes_b, shape_functions_b.row(2))); const auto a3_a = a1_a.cross(a2_a).normalized(); const auto a3_b = a1_b.cross(a2_b).normalized(); From 4f049fba8899d0f7cce351b9a0351ec8f80698a2 Mon Sep 17 00:00:00 2001 From: Thomas Oberbichler Date: Tue, 9 Jun 2020 22:34:55 +0200 Subject: [PATCH 19/19] Fix factor --- include/eqlib/objectives/IgaMembrane3PAD.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/eqlib/objectives/IgaMembrane3PAD.h b/include/eqlib/objectives/IgaMembrane3PAD.h index eabba8fa..5f31db8c 100644 --- a/include/eqlib/objectives/IgaMembrane3PAD.h +++ b/include/eqlib/objectives/IgaMembrane3PAD.h @@ -113,7 +113,7 @@ class IgaMembrane3PAD : public Objective { const auto eps = transformation_matrix * (act_a - ref_a).transpose() * 0.5; - const auto result = eps.dot(m_dm * eps) * weight; + const auto result = 0.5 * weight * eps.dot(m_dm * eps); const index a = shape_functions.cols() * 3;