From 43fbc25378a4d7e0ef2b72e3cc7e606d5dcdf416 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20D=C3=B6lz?= Date: Fri, 3 Jun 2022 10:27:46 +0200 Subject: [PATCH 01/20] Split storage of far field quadrature nodes from a single, large matrix to an std::vector of the size of the number of elements. This is to address cache based performance issues with larger refinement levels. Running the tests indicates that there is no significant change in performance. Logs below. Before changes: =============================================================================== Running tests... Test project /home/doelz/Dokumente/Codes/bembeldev_juergen/build Start 1: Spline 1/18 Test #1: Spline ............................ Passed 0.02 sec Start 2: GeometryImportAndEval 2/18 Test #2: GeometryImportAndEval ............. Passed 0.01 sec Start 3: SurfacePointUpdate 3/18 Test #3: SurfacePointUpdate ................ Passed 0.01 sec Start 4: Projector 4/18 Test #4: Projector ......................... Passed 0.02 sec Start 5: Glue 5/18 Test #5: Glue .............................. Passed 0.39 sec Start 6: FMMTransferMatrices 6/18 Test #6: FMMTransferMatrices ............... Passed 0.03 sec Start 7: FMMForwardTransformation 7/18 Test #7: FMMForwardTransformation .......... Passed 0.06 sec Start 8: FMMBackwardTransformation 8/18 Test #8: FMMBackwardTransformation ......... Passed 0.05 sec Start 9: DuffyTrick 9/18 Test #9: DuffyTrick ........................ Passed 10.56 sec Start 10: AnsatzSpaceExample 10/18 Test #10: AnsatzSpaceExample ................ Passed 0.01 sec Start 11: BlockClusterTreeExample 11/18 Test #11: BlockClusterTreeExample ........... Passed 12.52 sec Start 12: GeometryExample 12/18 Test #12: GeometryExample ................... Passed 0.01 sec Start 13: LaplaceSingleLayerFullExample 13/18 Test #13: LaplaceSingleLayerFullExample ..... Passed 21.49 sec Start 14: LaplaceSingleLayerH2Example 14/18 Test #14: LaplaceSingleLayerH2Example ....... Passed 19.89 sec Start 15: HelmholtzSingleLayerFullExample 15/18 Test #15: HelmholtzSingleLayerFullExample ... Passed 27.73 sec Start 16: HelmholtzSingleLayerH2Example 16/18 Test #16: HelmholtzSingleLayerH2Example ..... Passed 16.72 sec Start 17: MaxwellSingleLayerFullExample 17/18 Test #17: MaxwellSingleLayerFullExample ..... Passed 49.92 sec Start 18: MaxwellSingleLayerH2Example 18/18 Test #18: MaxwellSingleLayerH2Example ....... Passed 216.39 sec 100% tests passed, 0 tests failed out of 18 Total Test time (real) = 375.87 sec After changes: =============================================================================== Running tests... Test project /home/doelz/Dokumente/Codes/bembeldev_juergen/build Start 1: Spline 1/18 Test #1: Spline ............................ Passed 0.02 sec Start 2: GeometryImportAndEval 2/18 Test #2: GeometryImportAndEval ............. Passed 0.01 sec Start 3: SurfacePointUpdate 3/18 Test #3: SurfacePointUpdate ................ Passed 0.01 sec Start 4: Projector 4/18 Test #4: Projector ......................... Passed 0.02 sec Start 5: Glue 5/18 Test #5: Glue .............................. Passed 0.39 sec Start 6: FMMTransferMatrices 6/18 Test #6: FMMTransferMatrices ............... Passed 0.03 sec Start 7: FMMForwardTransformation 7/18 Test #7: FMMForwardTransformation .......... Passed 0.06 sec Start 8: FMMBackwardTransformation 8/18 Test #8: FMMBackwardTransformation ......... Passed 0.03 sec Start 9: DuffyTrick 9/18 Test #9: DuffyTrick ........................ Passed 10.42 sec Start 10: AnsatzSpaceExample 10/18 Test #10: AnsatzSpaceExample ................ Passed 0.02 sec Start 11: BlockClusterTreeExample 11/18 Test #11: BlockClusterTreeExample ........... Passed 12.42 sec Start 12: GeometryExample 12/18 Test #12: GeometryExample ................... Passed 0.01 sec Start 13: LaplaceSingleLayerFullExample 13/18 Test #13: LaplaceSingleLayerFullExample ..... Passed 17.39 sec Start 14: LaplaceSingleLayerH2Example 14/18 Test #14: LaplaceSingleLayerH2Example ....... Passed 19.76 sec Start 15: HelmholtzSingleLayerFullExample 15/18 Test #15: HelmholtzSingleLayerFullExample ... Passed 23.60 sec Start 16: HelmholtzSingleLayerH2Example 16/18 Test #16: HelmholtzSingleLayerH2Example ..... Passed 16.52 sec Start 17: MaxwellSingleLayerFullExample 17/18 Test #17: MaxwellSingleLayerFullExample ..... Passed 47.77 sec Start 18: MaxwellSingleLayerH2Example 18/18 Test #18: MaxwellSingleLayerH2Example ....... Passed 216.50 sec 100% tests passed, 0 tests failed out of 18 Total Test time (real) = 365.03 sec --- .../src/DuffyTrick/evaluateBilinearForm.hpp | 20 +++++++++++-------- .../DuffyTrick/farFieldQuadratureNodes.hpp | 13 +++++++----- Bembel/src/DuffyTrick/integrate0.hpp | 11 +++++----- Bembel/src/DuffyTrick/integrate1.hpp | 7 ++++--- Bembel/src/DuffyTrick/integrate2.hpp | 9 +++++---- Bembel/src/DuffyTrick/integrate3.hpp | 9 +++++---- Bembel/src/DuffyTrick/integrate4.hpp | 9 +++++---- Bembel/src/H2Matrix/H2Matrix.hpp | 3 ++- .../src/LinearOperator/DiscreteOperator.hpp | 11 +++++----- tests/DuffyTrick/test_integrate0.hpp | 3 ++- tests/DuffyTrick/test_integrate1.hpp | 3 ++- tests/DuffyTrick/test_integrate2.hpp | 3 ++- tests/DuffyTrick/test_integrate3.hpp | 2 +- tests/DuffyTrick/test_integrate4.hpp | 2 +- 14 files changed, 60 insertions(+), 45 deletions(-) diff --git a/Bembel/src/DuffyTrick/evaluateBilinearForm.hpp b/Bembel/src/DuffyTrick/evaluateBilinearForm.hpp index 8c99a12e1..9f125b164 100644 --- a/Bembel/src/DuffyTrick/evaluateBilinearForm.hpp +++ b/Bembel/src/DuffyTrick/evaluateBilinearForm.hpp @@ -19,7 +19,8 @@ template void evaluateBilinearForm( const LinearOperatorBase& linOp, const T& super_space, const ElementTreeNode& e1, const ElementTreeNode& e2, - const CubatureVector& GS, const Eigen::MatrixXd& ffield_qnodes, + const CubatureVector& GS, const Eigen::MatrixXd& ffield_qnodes1, + const Eigen::MatrixXd& ffield_qnodes2, Eigen::Matrix::Scalar, Eigen::Dynamic, Eigen::Dynamic>* intval) { ////////////////////////////////////////////////////////////////////////////// @@ -38,24 +39,27 @@ void evaluateBilinearForm( switch (cp(2)) { case 0: if (nfield_deg == ffield_deg) { - integrate0(linOp, super_space, e1, 0, e2, 0, ffield_qnodes, Q, intval); + integrate0(linOp, super_space, e1, 0, e2, 0, ffield_qnodes1, + ffield_qnodes2, Q, intval); return; } else { - integrate1(linOp, super_space, e1, 0, e2, 0, ffield_qnodes, Q, intval); + integrate1(linOp, super_space, e1, 0, e2, 0, ffield_qnodes1, + ffield_qnodes2, Q, intval); return; } case 1: assert(!"you should not have ended up here!"); case 2: - integrate2(linOp, super_space, e1, 0, e2, 0, ffield_qnodes, Q, intval); + integrate2(linOp, super_space, e1, 0, e2, 0, ffield_qnodes1, + ffield_qnodes2, Q, intval); return; case 3: - integrate3(linOp, super_space, e1, cp(0), e2, cp(1), ffield_qnodes, Q, - intval); + integrate3(linOp, super_space, e1, cp(0), e2, cp(1), ffield_qnodes1, + ffield_qnodes2, Q, intval); return; case 4: - integrate4(linOp, super_space, e1, cp(0), e2, cp(1), ffield_qnodes, Q, - intval); + integrate4(linOp, super_space, e1, cp(0), e2, cp(1), ffield_qnodes1, + ffield_qnodes2, Q, intval); return; default: assert(!"you should not have ended up here!"); diff --git a/Bembel/src/DuffyTrick/farFieldQuadratureNodes.hpp b/Bembel/src/DuffyTrick/farFieldQuadratureNodes.hpp index a976f3e93..3540ee0d0 100644 --- a/Bembel/src/DuffyTrick/farFieldQuadratureNodes.hpp +++ b/Bembel/src/DuffyTrick/farFieldQuadratureNodes.hpp @@ -17,26 +17,29 @@ namespace DuffyTrick { * is qNodes.col(k) = [xi, w, Chi(xi); dsChi(xi); dtChi(xi)]\in\Rbb^12 **/ template -Eigen::Matrix computeFfieldQnodes( +std::vector> computeFfieldQnodes( const T &super_space, const Cubature &Q) { - Eigen::Matrix ffield_qnodes; + std::vector> ffield_qnodes; int next = 0; // assume isotropic mesh width h! double h = (super_space.get_mesh().get_element_tree().cpbegin())->get_h(); auto nE = super_space.get_mesh().get_number_of_elements(); auto pbegin = super_space.get_mesh().get_element_tree().cpbegin(); auto pend = super_space.get_mesh().get_element_tree().cpend(); - ffield_qnodes.resize(12, nE * Q.xi_.cols()); + // ffield_qnodes.resize(12, nE * Q.xi_.cols()); + ffield_qnodes.reserve(nE); SurfacePoint surfpt; - for (auto it = pbegin; it != pend; ++it) + for (auto it = pbegin; it != pend; ++it) { + ffield_qnodes.emplace_back(12, Q.xi_.cols()); for (auto k = 0; k < Q.xi_.cols(); ++k) { // the quadrature weight is scaled by mesh width // this corresponds to a scaling of the basis functions // with respect to the L^2 norm! super_space.map2surface(*it, Q.xi_.col(k), h * Q.w_(k), &surfpt); - ffield_qnodes.col(it->id_ * Q.xi_.cols() + k) = surfpt; + ffield_qnodes[it->id_].col(k) = surfpt; } + } return ffield_qnodes; } } // namespace DuffyTrick diff --git a/Bembel/src/DuffyTrick/integrate0.hpp b/Bembel/src/DuffyTrick/integrate0.hpp index 2ef1f369e..f1146e968 100644 --- a/Bembel/src/DuffyTrick/integrate0.hpp +++ b/Bembel/src/DuffyTrick/integrate0.hpp @@ -19,21 +19,20 @@ namespace DuffyTrick { template void integrate0(const LinearOperatorBase &LinOp, const T &super_space, const ElementTreeNode &e1, int rot1, const ElementTreeNode &e2, - int rot2, const Eigen::MatrixXd &ffield_qnodes, - const Cubature &Q, + int rot2, const Eigen::MatrixXd &ffield_qnodes1, + const Eigen::MatrixXd &ffield_qnodes2, const Cubature &Q, Eigen::Matrix::Scalar, Eigen::Dynamic, Eigen::Dynamic> *intval) { intval->setZero(); for (auto i = 0; i < Q.w_.size(); ++i) for (auto j = 0; j < Q.w_.size(); ++j) - LinOp.evaluateIntegrand( - super_space, ffield_qnodes.col(e1.id_ * Q.w_.size() + i), - ffield_qnodes.col(e2.id_ * Q.w_.size() + j), intval); + LinOp.evaluateIntegrand(super_space, ffield_qnodes1.col(i), + ffield_qnodes2.col(j), intval); BEMBEL_UNUSED_(rot1); BEMBEL_UNUSED_(rot2); BEMBEL_UNUSED_(Q); return; } -} // namespace Duffy +} // namespace DuffyTrick } // namespace Bembel #endif diff --git a/Bembel/src/DuffyTrick/integrate1.hpp b/Bembel/src/DuffyTrick/integrate1.hpp index 2a865dfdf..9df1c055c 100644 --- a/Bembel/src/DuffyTrick/integrate1.hpp +++ b/Bembel/src/DuffyTrick/integrate1.hpp @@ -21,8 +21,8 @@ namespace DuffyTrick { template void integrate1(const LinearOperatorBase &LinOp, const T &super_space, const ElementTreeNode &e1, int rot1, const ElementTreeNode &e2, - int rot2, const Eigen::MatrixXd &ffield_qnodes, - const Cubature &Q, + int rot2, const Eigen::MatrixXd &ffield_qnodes1, + const Eigen::MatrixXd &ffield_qnodes2, const Cubature &Q, Eigen::Matrix::Scalar, Eigen::Dynamic, Eigen::Dynamic> *intval) { double h = e1.get_h(); @@ -38,7 +38,8 @@ void integrate1(const LinearOperatorBase &LinOp, const T &super_space, BEMBEL_UNUSED_(rot1); BEMBEL_UNUSED_(rot2); - BEMBEL_UNUSED_(ffield_qnodes); + BEMBEL_UNUSED_(ffield_qnodes1); + BEMBEL_UNUSED_(ffield_qnodes2); return; } } // namespace DuffyTrick diff --git a/Bembel/src/DuffyTrick/integrate2.hpp b/Bembel/src/DuffyTrick/integrate2.hpp index 2b88e988d..84272a8dc 100644 --- a/Bembel/src/DuffyTrick/integrate2.hpp +++ b/Bembel/src/DuffyTrick/integrate2.hpp @@ -24,8 +24,8 @@ namespace DuffyTrick { template void integrate2(const LinearOperatorBase &LinOp, const T &super_space, const ElementTreeNode &e1, int rot1, const ElementTreeNode &e2, - int rot2, const Eigen::MatrixXd &ffield_qnodes, - const Cubature &Q, + int rot2, const Eigen::MatrixXd &ffield_qnodes1, + const Eigen::MatrixXd &ffield_qnodes2, const Cubature &Q, Eigen::Matrix::Scalar, Eigen::Dynamic, Eigen::Dynamic> *intval) { intval->setZero(); @@ -58,10 +58,11 @@ void integrate2(const LinearOperatorBase &LinOp, const T &super_space, BEMBEL_UNUSED_(e2); BEMBEL_UNUSED_(rot1); BEMBEL_UNUSED_(rot2); - BEMBEL_UNUSED_(ffield_qnodes); + BEMBEL_UNUSED_(ffield_qnodes1); + BEMBEL_UNUSED_(ffield_qnodes2); return; } -} // namespace Duffy +} // namespace DuffyTrick } // namespace Bembel #endif diff --git a/Bembel/src/DuffyTrick/integrate3.hpp b/Bembel/src/DuffyTrick/integrate3.hpp index 09f7e527c..f323875a3 100644 --- a/Bembel/src/DuffyTrick/integrate3.hpp +++ b/Bembel/src/DuffyTrick/integrate3.hpp @@ -24,8 +24,8 @@ namespace DuffyTrick { template void integrate3(const LinearOperatorBase &LinOp, const T &super_space, const ElementTreeNode &e1, int rot1, const ElementTreeNode &e2, - int rot2, const Eigen::MatrixXd &ffield_qnodes, - const Cubature &Q, + int rot2, const Eigen::MatrixXd &ffield_qnodes1, + const Eigen::MatrixXd &ffield_qnodes2, const Cubature &Q, Eigen::Matrix::Scalar, Eigen::Dynamic, Eigen::Dynamic> *intval) { intval->setZero(); @@ -71,10 +71,11 @@ void integrate3(const LinearOperatorBase &LinOp, const T &super_space, LinOp.evaluateIntegrand(super_space, qp1, qp2, intval); } } - BEMBEL_UNUSED_(ffield_qnodes); + BEMBEL_UNUSED_(ffield_qnodes1); + BEMBEL_UNUSED_(ffield_qnodes2); return; } -} // namespace Duffy +} // namespace DuffyTrick } // namespace Bembel #endif diff --git a/Bembel/src/DuffyTrick/integrate4.hpp b/Bembel/src/DuffyTrick/integrate4.hpp index 53787ea9c..833499d28 100644 --- a/Bembel/src/DuffyTrick/integrate4.hpp +++ b/Bembel/src/DuffyTrick/integrate4.hpp @@ -24,8 +24,8 @@ namespace DuffyTrick { template void integrate4(const LinearOperatorBase &LinOp, const T &super_space, const ElementTreeNode &e1, int rot1, const ElementTreeNode &e2, - int rot2, const Eigen::MatrixXd &ffield_qnodes, - const Cubature &Q, + int rot2, const Eigen::MatrixXd &ffield_qnodes1, + const Eigen::MatrixXd &ffield_qnodes2, const Cubature &Q, Eigen::Matrix::Scalar, Eigen::Dynamic, Eigen::Dynamic> *intval) { intval->setZero(); @@ -55,10 +55,11 @@ void integrate4(const LinearOperatorBase &LinOp, const T &super_space, LinOp.evaluateIntegrand(super_space, qp6, qp4, intval); } } - BEMBEL_UNUSED_(ffield_qnodes); + BEMBEL_UNUSED_(ffield_qnodes1); + BEMBEL_UNUSED_(ffield_qnodes2); return; } -} // namespace Duffy +} // namespace DuffyTrick } // namespace Bembel #endif diff --git a/Bembel/src/H2Matrix/H2Matrix.hpp b/Bembel/src/H2Matrix/H2Matrix.hpp index d77360b17..77c539f51 100644 --- a/Bembel/src/H2Matrix/H2Matrix.hpp +++ b/Bembel/src/H2Matrix/H2Matrix.hpp @@ -183,7 +183,8 @@ class H2Matrix : public EigenBase> { // do integration Bembel::DuffyTrick::evaluateBilinearForm( linOp, super_space, element1, element2, GS, - ffield_qnodes, &intval); + ffield_qnodes[element1.id_], + ffield_qnodes[element2.id_], &intval); // insert into dense matrices of all block cluster trees for (int i = 0; i < vector_dimension; ++i) for (int j = 0; j < vector_dimension; ++j) diff --git a/Bembel/src/LinearOperator/DiscreteOperator.hpp b/Bembel/src/LinearOperator/DiscreteOperator.hpp index 67484b2c3..0b06d4186 100644 --- a/Bembel/src/LinearOperator/DiscreteOperator.hpp +++ b/Bembel/src/LinearOperator/DiscreteOperator.hpp @@ -29,8 +29,8 @@ struct DiscreteOperatorComputer< Eigen::Matrix *disc_op, const Derived &lin_op, const AnsatzSpace &ansatz_space) { GaussSquare GS; - const SuperSpace& super_space = ansatz_space.get_superspace(); - const ElementTree& element_tree = super_space.get_mesh().get_element_tree(); + const SuperSpace &super_space = ansatz_space.get_superspace(); + const ElementTree &element_tree = super_space.get_mesh().get_element_tree(); auto number_of_elements = element_tree.get_number_of_elements(); const auto vector_dimension = getFunctionSpaceVectorDimension::Form>(); @@ -57,9 +57,10 @@ struct DiscreteOperatorComputer< Eigen::Matrix intval( vector_dimension * polynomial_degree_plus_one_squared, vector_dimension * polynomial_degree_plus_one_squared); - DuffyTrick::evaluateBilinearForm(lin_op, super_space, *element1, - *element2, GS, ffield_qnodes, - &intval); + DuffyTrick::evaluateBilinearForm( + lin_op, super_space, *element1, *element2, GS, + ffield_qnodes[element1->id_], ffield_qnodes[element2->id_], + &intval); for (auto i = 0; i < vector_dimension; ++i) for (auto j = 0; j < vector_dimension; ++j) disc_op->block(polynomial_degree_plus_one_squared * diff --git a/tests/DuffyTrick/test_integrate0.hpp b/tests/DuffyTrick/test_integrate0.hpp index fa890292d..bb9ed5339 100644 --- a/tests/DuffyTrick/test_integrate0.hpp +++ b/tests/DuffyTrick/test_integrate0.hpp @@ -36,7 +36,8 @@ bool test_integrate0(const Bembel::AnsatzSpace &ansatz_space, intval.setZero(); error = 0; Bembel::DuffyTrick::integrate0(linOp, ansatz_space.get_superspace(), *it, - 0, *it2, 0, ffield_qnodes, Q, &intval); + 0, *it2, 0, ffield_qnodes[it->id_], + ffield_qnodes[it2->id_], Q, &intval); axis.col(0) << it->llc_(0), it->llc_(0) + h; axis.col(1) << it->llc_(1), it->llc_(1) + h; axis.col(2) << it2->llc_(0), it2->llc_(0) + h; diff --git a/tests/DuffyTrick/test_integrate1.hpp b/tests/DuffyTrick/test_integrate1.hpp index 1d40d675e..d0d0327dd 100644 --- a/tests/DuffyTrick/test_integrate1.hpp +++ b/tests/DuffyTrick/test_integrate1.hpp @@ -35,7 +35,8 @@ bool test_integrate1(const Bembel::AnsatzSpace &ansatz_space, intval.setZero(); error = 0; Bembel::DuffyTrick::integrate1(linOp, ansatz_space.get_superspace(), *it, - 0, *it2, 0, ffield_qnodes, Q, &intval); + 0, *it2, 0, ffield_qnodes, ffield_qnodes, + Q, &intval); axis.col(0) << it->llc_(0), it->llc_(0) + h; axis.col(1) << it->llc_(1), it->llc_(1) + h; axis.col(2) << it2->llc_(0), it2->llc_(0) + h; diff --git a/tests/DuffyTrick/test_integrate2.hpp b/tests/DuffyTrick/test_integrate2.hpp index 58716f689..fbebf7778 100644 --- a/tests/DuffyTrick/test_integrate2.hpp +++ b/tests/DuffyTrick/test_integrate2.hpp @@ -37,7 +37,8 @@ bool test_integrate2(const Bembel::AnsatzSpace &ansatz_space, error = 0; //////////////////////////////////////////////////////////////////////////// Bembel::DuffyTrick::integrate2(linOp, ansatz_space.get_superspace(), *it, 0, - *it, 0, ffield_qnodes, Q, &intval); + *it, 0, ffield_qnodes, ffield_qnodes, Q, + &intval); //////////////////////////////////////////////////////////////////////////// axis.col(0) << it->llc_(0), it->llc_(0) + h; axis.col(1) << it->llc_(1), it->llc_(1) + h; diff --git a/tests/DuffyTrick/test_integrate3.hpp b/tests/DuffyTrick/test_integrate3.hpp index d7b97934f..2060801db 100644 --- a/tests/DuffyTrick/test_integrate3.hpp +++ b/tests/DuffyTrick/test_integrate3.hpp @@ -43,7 +43,7 @@ bool test_integrate3(const Bembel::AnsatzSpace &ansatz_space, //////////////////////////////////////////////////////////////////////// Bembel::DuffyTrick::integrate3(linOp, ansatz_space.get_superspace(), *it, cp(0), *it2, cp(1), ffield_qnodes, - Q, &intval); + ffield_qnodes, Q, &intval); //////////////////////////////////////////////////////////////////////// axis.col(0) << it->llc_(0), it->llc_(0) + h; axis.col(1) << it->llc_(1), it->llc_(1) + h; diff --git a/tests/DuffyTrick/test_integrate4.hpp b/tests/DuffyTrick/test_integrate4.hpp index 3d4483a57..14af1628f 100644 --- a/tests/DuffyTrick/test_integrate4.hpp +++ b/tests/DuffyTrick/test_integrate4.hpp @@ -42,7 +42,7 @@ bool test_integrate4(const Bembel::AnsatzSpace &ansatz_space, //////////////////////////////////////////////////////////////////////// Bembel::DuffyTrick::integrate4(linOp, ansatz_space.get_superspace(), *it, cp(0), *it2, cp(1), ffield_qnodes, - Q, &intval); + ffield_qnodes, Q, &intval); //////////////////////////////////////////////////////////////////////// axis.col(0) << it->llc_(0), it->llc_(0) + h; axis.col(1) << it->llc_(1), it->llc_(1) + h; From 8500db2a0a517eb1b916cc660711a0056e4cad67 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20D=C3=B6lz?= Date: Sun, 10 Jul 2022 18:14:00 +0200 Subject: [PATCH 02/20] added DirichletTrace for DiscreteLinearFormOnReferenceDomain --- Bembel/LinearForm | 2 + .../DirichletTraceOnReferenceDomain.hpp | 72 ++++++++++++ .../DiscreteLinearFormOnReferenceDomain.hpp | 107 ++++++++++++++++++ 3 files changed, 181 insertions(+) create mode 100644 Bembel/src/LinearForm/DirichletTraceOnReferenceDomain.hpp create mode 100644 Bembel/src/LinearForm/DiscreteLinearFormOnReferenceDomain.hpp diff --git a/Bembel/LinearForm b/Bembel/LinearForm index ef7666871..cc8d942a7 100644 --- a/Bembel/LinearForm +++ b/Bembel/LinearForm @@ -12,8 +12,10 @@ #include "Quadrature" #include "src/LinearForm/LinearForm.hpp" +#include "src/LinearForm/DiscreteLinearFormOnReferenceDomain.hpp" #include "src/LinearForm/DirichletTrace.hpp" +#include "src/LinearForm/DirichletTraceOnReferenceDomain.hpp" #include "src/LinearForm/DiscreteLinearForm.hpp" #include "src/LinearForm/RotatedTangentialTrace.hpp" diff --git a/Bembel/src/LinearForm/DirichletTraceOnReferenceDomain.hpp b/Bembel/src/LinearForm/DirichletTraceOnReferenceDomain.hpp new file mode 100644 index 000000000..9f417ae9b --- /dev/null +++ b/Bembel/src/LinearForm/DirichletTraceOnReferenceDomain.hpp @@ -0,0 +1,72 @@ +// This file is part of Bembel, the higher order C++ boundary element library. +// It was written as part of a cooperation of J. Doelz, H. Harbrecht, S. Kurz, +// M. Multerer, S. Schoeps, and F. Wolf at Technische Universitaet Darmstadt, +// Universitaet Basel, and Universita della Svizzera italiana, Lugano. This +// source code is subject to the GNU General Public License version 3 and +// provided WITHOUT ANY WARRANTY, see for further +// information. +#ifndef BEMBEL_LINEARFORM_DIRICHLETTRACEONREFERENCEDOMAIN_H_ +#define BEMBEL_LINEARFORM_DIRICHLETTRACEONREFERENCEDOMAIN_H_ + +namespace Bembel { + +template +class DirichletTraceOnReferenceDomain; + +template +struct LinearFormTraits> { + typedef ScalarT Scalar; +}; + +/** + * \ingroup LinearForm + * \brief This class provides an implementation of the Dirichlet trace operator + * and a corresponding method to evaluate the linear form corresponding to the + * right hand side of the system via quadrature. + */ +template +class DirichletTraceOnReferenceDomain + : public LinearFormBase, Scalar> { + public: + DirichletTraceOnReferenceDomain() {} + void set_function( + const std::function &function) { + function_ = function; + } + template + void evaluateIntegrand_impl(const T &super_space, const SurfacePoint &p, + Eigen::Matrix *intval, + int patch) const { + auto polynomial_degree = super_space.get_polynomial_degree(); + auto polynomial_degree_plus_one_squared = + (polynomial_degree + 1) * (polynomial_degree + 1); + + // get evaluation points on unit square + auto s = p.segment<2>(0); + + // get quadrature weights + auto ws = p(2); + + // get points on geometry and tangential derivatives + auto x_f = p.segment<3>(3); + auto x_f_dx = p.segment<3>(6); + auto x_f_dy = p.segment<3>(9); + + // compute surface measures from tangential derivatives + auto x_kappa = x_f_dx.cross(x_f_dy).norm(); + + // integrand without basis functions + auto integrand = function_(patch, s) * x_kappa * ws; + + // multiply basis functions with integrand + super_space.addScaledBasis(intval, integrand, s); + + return; + }; + + private: + std::function function_; +}; +} // namespace Bembel + +#endif diff --git a/Bembel/src/LinearForm/DiscreteLinearFormOnReferenceDomain.hpp b/Bembel/src/LinearForm/DiscreteLinearFormOnReferenceDomain.hpp new file mode 100644 index 000000000..8f9a1536c --- /dev/null +++ b/Bembel/src/LinearForm/DiscreteLinearFormOnReferenceDomain.hpp @@ -0,0 +1,107 @@ +// This file is part of Bembel, the higher order C++ boundary element library. +// It was written as part of a cooperation of J. Doelz, H. Harbrecht, S. Kurz, +// M. Multerer, S. Schoeps, and F. Wolf at Technische Universitaet Darmstadt, +// Universitaet Basel, and Universita della Svizzera italiana, Lugano. This +// source code is subject to the GNU General Public License version 3 and +// provided WITHOUT ANY WARRANTY, see for further +// information. +#ifndef BEMBEL_LINEARFORM_DISCRETELINEARFORMONREFERENCEDOMAIN_H_ +#define BEMBEL_LINEARFORM_DISCRETELINEARFORMONREFERENCEDOMAIN_H_ + +namespace Bembel { +/** + * \ingroup LinearForm + * \brief This class, given a LinearForm with defined traits, takes care of the + * assembly of the right hand side. + */ +template +class DiscreteLinearFormOnReferenceDomain { + public: + ////////////////////////////////////////////////////////////////////////////// + // constructors + ////////////////////////////////////////////////////////////////////////////// + DiscreteLinearFormOnReferenceDomain() {} + DiscreteLinearFormOnReferenceDomain(const AnsatzSpace &ansatz_space) { + init_DiscreteLinearFormOnReferenceDomain(ansatz_space); + } + ////////////////////////////////////////////////////////////////////////////// + // init_DiscreteLinearFormOnReferenceDomain + ////////////////////////////////////////////////////////////////////////////// + void init_DiscreteLinearFormOnReferenceDomain( + const AnsatzSpace &ansatz_space) { + ansatz_space_ = ansatz_space; + deg_ = ansatz_space_.get_polynomial_degree() + 1; + return; + } + ////////////////////////////////////////////////////////////////////////////// + // compute + ////////////////////////////////////////////////////////////////////////////// + /** + * \todo Add inline commentary + **/ + void compute() { + GaussSquare GS; + const auto &Q = GS[deg_]; + SurfacePoint qp; + const auto &super_space = ansatz_space_.get_superspace(); + const auto &element_tree = super_space.get_mesh().get_element_tree(); + auto number_of_elements = element_tree.get_number_of_elements(); + auto polynomial_degree = super_space.get_polynomial_degree(); + auto polynomial_degree_plus_one_squared = + (polynomial_degree + 1) * (polynomial_degree + 1); + const auto function_space_dimension = + getFunctionSpaceVectorDimension::Form>(); + Eigen::Matrix::Scalar, Eigen::Dynamic, + function_space_dimension> + intval(polynomial_degree_plus_one_squared, function_space_dimension); + Eigen::Matrix::Scalar, Eigen::Dynamic, + function_space_dimension> + disc_lf_matrix(polynomial_degree_plus_one_squared * number_of_elements, + function_space_dimension); + disc_lf_matrix.setZero(); + for (auto element = element_tree.cpbegin(); element != element_tree.cpend(); + ++element) { + intval.setZero(); + for (auto i = 0; i < Q.w_.size(); ++i) { + super_space.map2surface(*element, Q.xi_.col(i), + element->get_h() * Q.w_(i), &qp); + lf_.evaluateIntegrand_impl(super_space, qp, &intval, element->patch_); + } + disc_lf_matrix.block(polynomial_degree_plus_one_squared * element->id_, 0, + polynomial_degree_plus_one_squared, + function_space_dimension) = intval; + } + disc_lf_ = + ansatz_space_.get_transformation_matrix().transpose() * + Eigen::Map::Scalar, + Eigen::Dynamic, 1>>( + disc_lf_matrix.data(), + disc_lf_matrix.rows() * disc_lf_matrix.cols()); + return; + } + ////////////////////////////////////////////////////////////////////////////// + // setter + ////////////////////////////////////////////////////////////////////////////// + void set_degree(const int °) { deg_ = deg; } + ////////////////////////////////////////////////////////////////////////////// + // getter + ////////////////////////////////////////////////////////////////////////////// + Derived &get_linear_form() { return lf_; } + const Eigen::Matrix::Scalar, + Eigen::Dynamic, 1> + &get_discrete_linear_form() const { + return disc_lf_; + } + ////////////////////////////////////////////////////////////////////////////// + // private member variables + ////////////////////////////////////////////////////////////////////////////// + private: + int deg_; + Derived lf_; + Eigen::Matrix::Scalar, Eigen::Dynamic, 1> + disc_lf_; + AnsatzSpace ansatz_space_; +}; + +} // namespace Bembel +#endif From 0899e122bde62042b45d022e3af61e3e32b4da41 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20D=C3=B6lz?= Date: Mon, 30 Jan 2023 20:48:54 +0100 Subject: [PATCH 03/20] use typedef SurfacePoint consistently throughout the code --- Bembel/Geometry | 4 ++-- Bembel/src/DuffyTrick/evaluateBilinearForm.hpp | 4 ++-- Bembel/src/DuffyTrick/farFieldQuadratureNodes.hpp | 13 ++++++------- Bembel/src/DuffyTrick/integrate0.hpp | 9 +++++---- Bembel/src/DuffyTrick/integrate1.hpp | 5 +++-- Bembel/src/DuffyTrick/integrate2.hpp | 5 +++-- Bembel/src/DuffyTrick/integrate3.hpp | 5 +++-- Bembel/src/DuffyTrick/integrate4.hpp | 5 +++-- Bembel/src/Geometry/Patch.hpp | 2 +- Bembel/src/Geometry/SurfacePoint.hpp | 2 +- tests/DuffyTrick/test_integrate1.hpp | 2 +- tests/DuffyTrick/test_integrate2.hpp | 2 +- tests/DuffyTrick/test_integrate3.hpp | 2 +- tests/DuffyTrick/test_integrate4.hpp | 2 +- tests/test_SurfacePointUpdate.cpp | 2 +- 15 files changed, 34 insertions(+), 30 deletions(-) diff --git a/Bembel/Geometry b/Bembel/Geometry index 6734770d0..aa18eb5ab 100644 --- a/Bembel/Geometry +++ b/Bembel/Geometry @@ -37,11 +37,11 @@ #include "Spline" +#include "src/Geometry/SurfacePoint.hpp" + #include "src/Geometry/Patch.hpp" #include "src/Geometry/PatchVector.hpp" #include "src/Geometry/GeometryIO.hpp" #include "src/Geometry/Geometry.hpp" -#include "src/Geometry/SurfacePoint.hpp" -#include "src/Geometry/Geometry.hpp" #endif // BEMBEL_GEOMETRY_MODULE_ diff --git a/Bembel/src/DuffyTrick/evaluateBilinearForm.hpp b/Bembel/src/DuffyTrick/evaluateBilinearForm.hpp index 92bffd7bd..e1d7f55cb 100644 --- a/Bembel/src/DuffyTrick/evaluateBilinearForm.hpp +++ b/Bembel/src/DuffyTrick/evaluateBilinearForm.hpp @@ -22,8 +22,8 @@ template void evaluateBilinearForm( const LinearOperatorBase& linOp, const T& super_space, const ElementTreeNode& e1, const ElementTreeNode& e2, - const CubatureVector& GS, const Eigen::MatrixXd& ffield_qnodes1, - const Eigen::MatrixXd& ffield_qnodes2, + const CubatureVector& GS, const std::vector& ffield_qnodes1, + const std::vector& ffield_qnodes2, Eigen::Matrix::Scalar, Eigen::Dynamic, Eigen::Dynamic>* intval) { ////////////////////////////////////////////////////////////////////////////// diff --git a/Bembel/src/DuffyTrick/farFieldQuadratureNodes.hpp b/Bembel/src/DuffyTrick/farFieldQuadratureNodes.hpp index 37661fbf3..41c05d379 100644 --- a/Bembel/src/DuffyTrick/farFieldQuadratureNodes.hpp +++ b/Bembel/src/DuffyTrick/farFieldQuadratureNodes.hpp @@ -20,27 +20,26 @@ namespace DuffyTrick { * is qNodes.col(k) = [xi, w, Chi(xi); dsChi(xi); dtChi(xi)]\in\Rbb^12 **/ template -std::vector> computeFfieldQnodes( - const T &super_space, const Cubature &Q) { - std::vector> ffield_qnodes; +std::vector> computeFfieldQnodes(const T &super_space, + const Cubature &Q) { + std::vector> ffield_qnodes; int next = 0; // assume isotropic mesh width h! double h = (super_space.get_mesh().get_element_tree().cpbegin())->get_h(); auto nE = super_space.get_mesh().get_number_of_elements(); auto pbegin = super_space.get_mesh().get_element_tree().cpbegin(); auto pend = super_space.get_mesh().get_element_tree().cpend(); - // ffield_qnodes.resize(12, nE * Q.xi_.cols()); ffield_qnodes.reserve(nE); SurfacePoint surfpt; for (auto it = pbegin; it != pend; ++it) { - ffield_qnodes.emplace_back(12, Q.xi_.cols()); + ffield_qnodes.emplace_back(Q.xi_.cols()); for (auto k = 0; k < Q.xi_.cols(); ++k) { // the quadrature weight is scaled by mesh width // this corresponds to a scaling of the basis functions // with respect to the L^2 norm! - super_space.map2surface(*it, Q.xi_.col(k), h * Q.w_(k), &surfpt); - ffield_qnodes[it->id_].col(k) = surfpt; + super_space.map2surface(*it, Q.xi_.col(k), h * Q.w_(k), + &ffield_qnodes[it->id_][k]); } } return ffield_qnodes; diff --git a/Bembel/src/DuffyTrick/integrate0.hpp b/Bembel/src/DuffyTrick/integrate0.hpp index e46e060e9..f095038d3 100644 --- a/Bembel/src/DuffyTrick/integrate0.hpp +++ b/Bembel/src/DuffyTrick/integrate0.hpp @@ -22,15 +22,16 @@ namespace DuffyTrick { template void integrate0(const LinearOperatorBase &LinOp, const T &super_space, const ElementTreeNode &e1, int rot1, const ElementTreeNode &e2, - int rot2, const Eigen::MatrixXd &ffield_qnodes1, - const Eigen::MatrixXd &ffield_qnodes2, const Cubature &Q, + int rot2, const std::vector &ffield_qnodes1, + const std::vector &ffield_qnodes2, + const Cubature &Q, Eigen::Matrix::Scalar, Eigen::Dynamic, Eigen::Dynamic> *intval) { intval->setZero(); for (auto i = 0; i < Q.w_.size(); ++i) for (auto j = 0; j < Q.w_.size(); ++j) - LinOp.evaluateIntegrand(super_space, ffield_qnodes1.col(i), - ffield_qnodes2.col(j), intval); + LinOp.evaluateIntegrand(super_space, ffield_qnodes1[i], ffield_qnodes2[j], + intval); BEMBEL_UNUSED_(rot1); BEMBEL_UNUSED_(rot2); BEMBEL_UNUSED_(Q); diff --git a/Bembel/src/DuffyTrick/integrate1.hpp b/Bembel/src/DuffyTrick/integrate1.hpp index f369d9e82..fa37711f4 100644 --- a/Bembel/src/DuffyTrick/integrate1.hpp +++ b/Bembel/src/DuffyTrick/integrate1.hpp @@ -24,8 +24,9 @@ namespace DuffyTrick { template void integrate1(const LinearOperatorBase &LinOp, const T &super_space, const ElementTreeNode &e1, int rot1, const ElementTreeNode &e2, - int rot2, const Eigen::MatrixXd &ffield_qnodes1, - const Eigen::MatrixXd &ffield_qnodes2, const Cubature &Q, + int rot2, const std::vector &ffield_qnodes1, + const std::vector &ffield_qnodes2, + const Cubature &Q, Eigen::Matrix::Scalar, Eigen::Dynamic, Eigen::Dynamic> *intval) { double h = e1.get_h(); diff --git a/Bembel/src/DuffyTrick/integrate2.hpp b/Bembel/src/DuffyTrick/integrate2.hpp index 73d1c3cd0..accf64bcb 100644 --- a/Bembel/src/DuffyTrick/integrate2.hpp +++ b/Bembel/src/DuffyTrick/integrate2.hpp @@ -27,8 +27,9 @@ namespace DuffyTrick { template void integrate2(const LinearOperatorBase &LinOp, const T &super_space, const ElementTreeNode &e1, int rot1, const ElementTreeNode &e2, - int rot2, const Eigen::MatrixXd &ffield_qnodes1, - const Eigen::MatrixXd &ffield_qnodes2, const Cubature &Q, + int rot2, const std::vector &ffield_qnodes1, + const std::vector &ffield_qnodes2, + const Cubature &Q, Eigen::Matrix::Scalar, Eigen::Dynamic, Eigen::Dynamic> *intval) { intval->setZero(); diff --git a/Bembel/src/DuffyTrick/integrate3.hpp b/Bembel/src/DuffyTrick/integrate3.hpp index 9eb839df2..cd0cc89d0 100644 --- a/Bembel/src/DuffyTrick/integrate3.hpp +++ b/Bembel/src/DuffyTrick/integrate3.hpp @@ -27,8 +27,9 @@ namespace DuffyTrick { template void integrate3(const LinearOperatorBase &LinOp, const T &super_space, const ElementTreeNode &e1, int rot1, const ElementTreeNode &e2, - int rot2, const Eigen::MatrixXd &ffield_qnodes1, - const Eigen::MatrixXd &ffield_qnodes2, const Cubature &Q, + int rot2, const std::vector &ffield_qnodes1, + const std::vector &ffield_qnodes2, + const Cubature &Q, Eigen::Matrix::Scalar, Eigen::Dynamic, Eigen::Dynamic> *intval) { intval->setZero(); diff --git a/Bembel/src/DuffyTrick/integrate4.hpp b/Bembel/src/DuffyTrick/integrate4.hpp index cfbd1968f..603fef5d1 100644 --- a/Bembel/src/DuffyTrick/integrate4.hpp +++ b/Bembel/src/DuffyTrick/integrate4.hpp @@ -27,8 +27,9 @@ namespace DuffyTrick { template void integrate4(const LinearOperatorBase &LinOp, const T &super_space, const ElementTreeNode &e1, int rot1, const ElementTreeNode &e2, - int rot2, const Eigen::MatrixXd &ffield_qnodes1, - const Eigen::MatrixXd &ffield_qnodes2, const Cubature &Q, + int rot2, const std::vector &ffield_qnodes1, + const std::vector &ffield_qnodes2, + const Cubature &Q, Eigen::Matrix::Scalar, Eigen::Dynamic, Eigen::Dynamic> *intval) { intval->setZero(); diff --git a/Bembel/src/Geometry/Patch.hpp b/Bembel/src/Geometry/Patch.hpp index c6cf4259a..ffd44df06 100644 --- a/Bembel/src/Geometry/Patch.hpp +++ b/Bembel/src/Geometry/Patch.hpp @@ -209,7 +209,7 @@ class Patch { // This is a combination of eval und evalJacobian, to avoid duplication of // work. See SurfacePoint.hpp - void updateSurfacePoint(Eigen::Matrix *srf_pt, + void updateSurfacePoint(SurfacePoint *srf_pt, const Eigen::Vector2d &ref_pt, double w, const Eigen::Vector2d &xi) const { const int x_location = diff --git a/Bembel/src/Geometry/SurfacePoint.hpp b/Bembel/src/Geometry/SurfacePoint.hpp index 2322be945..5aaf62cdf 100644 --- a/Bembel/src/Geometry/SurfacePoint.hpp +++ b/Bembel/src/Geometry/SurfacePoint.hpp @@ -38,6 +38,6 @@ * updateSurdacePoint method is specialized and should be used, since it avoids * redundant work. **/ -typedef Eigen::Matrix SurfacePoint; +typedef Eigen::Matrix SurfacePoint; #endif // BEMBEL_SRC_GEOMETRY_SURFACEPOINT_HPP_ diff --git a/tests/DuffyTrick/test_integrate1.hpp b/tests/DuffyTrick/test_integrate1.hpp index 756bc9f36..41ab4c339 100644 --- a/tests/DuffyTrick/test_integrate1.hpp +++ b/tests/DuffyTrick/test_integrate1.hpp @@ -20,7 +20,7 @@ bool test_integrate1(const Bembel::AnsatzSpace &ansatz_space, const Bembel::LinearOperatorBase &linOp) { Bembel::GaussSquare GS; auto Q = GS[maxqdeg]; - Eigen::MatrixXd ffield_qnodes(0, 0); + std::vector ffield_qnodes; Eigen::MatrixXd intval; Eigen::MatrixXd axis; intval.resize(1, 1); diff --git a/tests/DuffyTrick/test_integrate2.hpp b/tests/DuffyTrick/test_integrate2.hpp index ead39bed7..a4ecbd531 100644 --- a/tests/DuffyTrick/test_integrate2.hpp +++ b/tests/DuffyTrick/test_integrate2.hpp @@ -21,7 +21,7 @@ bool test_integrate2(const Bembel::AnsatzSpace &ansatz_space, Bembel::GaussSquare GS; auto Q = GS[maxqdeg]; - Eigen::MatrixXd ffield_qnodes(0, 0); + std::vector ffield_qnodes; Eigen::MatrixXd intval; Eigen::MatrixXd axis; intval.resize(1, 1); diff --git a/tests/DuffyTrick/test_integrate3.hpp b/tests/DuffyTrick/test_integrate3.hpp index 0a4602550..98e8e5f05 100644 --- a/tests/DuffyTrick/test_integrate3.hpp +++ b/tests/DuffyTrick/test_integrate3.hpp @@ -21,7 +21,7 @@ bool test_integrate3(const Bembel::AnsatzSpace &ansatz_space, Bembel::GaussSquare GS; auto Q = GS[maxqdeg]; - Eigen::MatrixXd ffield_qnodes(0, 0); + std::vector ffield_qnodes; Eigen::MatrixXd intval; Eigen::MatrixXd axis; intval.resize(1, 1); diff --git a/tests/DuffyTrick/test_integrate4.hpp b/tests/DuffyTrick/test_integrate4.hpp index b6f23a464..2b14f9c74 100644 --- a/tests/DuffyTrick/test_integrate4.hpp +++ b/tests/DuffyTrick/test_integrate4.hpp @@ -20,7 +20,7 @@ bool test_integrate4(const Bembel::AnsatzSpace &ansatz_space, const Bembel::LinearOperatorBase &linOp) { Bembel::GaussSquare GS; auto Q = GS[maxqdeg]; - Eigen::MatrixXd ffield_qnodes(0, 0); + std::vector ffield_qnodes; Eigen::MatrixXd intval; Eigen::MatrixXd axis; intval.resize(1, 1); diff --git a/tests/test_SurfacePointUpdate.cpp b/tests/test_SurfacePointUpdate.cpp index 06d84e4fd..166ce7ff7 100644 --- a/tests/test_SurfacePointUpdate.cpp +++ b/tests/test_SurfacePointUpdate.cpp @@ -24,7 +24,7 @@ int main() { for (auto x : Test::Constants::eq_points) { for (auto y : Test::Constants::eq_points) { auto pt = Eigen::Vector2d(x, y); - Eigen::Matrix srf_pt, srf_pt_ref; + SurfacePoint srf_pt, srf_pt_ref; srf_pt_ref.head(2) = pt; srf_pt_ref(2) = 3.1415; From 68d9ab3b7f0cbf92f12ce0cae7ed6fbcffb95fca Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20D=C3=B6lz?= Date: Mon, 30 Jan 2023 21:18:52 +0100 Subject: [PATCH 04/20] Typo --- Bembel/src/Geometry/SurfacePoint.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Bembel/src/Geometry/SurfacePoint.hpp b/Bembel/src/Geometry/SurfacePoint.hpp index 5aaf62cdf..2322be945 100644 --- a/Bembel/src/Geometry/SurfacePoint.hpp +++ b/Bembel/src/Geometry/SurfacePoint.hpp @@ -38,6 +38,6 @@ * updateSurdacePoint method is specialized and should be used, since it avoids * redundant work. **/ -typedef Eigen::Matrix SurfacePoint; +typedef Eigen::Matrix SurfacePoint; #endif // BEMBEL_SRC_GEOMETRY_SURFACEPOINT_HPP_ From e5b5c2d1f80a7b38502a3a8e27f7be21cecfbc77 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20D=C3=B6lz?= Date: Wed, 1 Feb 2023 19:52:53 +0100 Subject: [PATCH 05/20] Introduced typedef for improved compatiblity with older compilers --- Bembel/src/DuffyTrick/farFieldQuadratureNodes.hpp | 6 +++--- Bembel/src/DuffyTrick/integrate0.hpp | 5 ++--- Bembel/src/DuffyTrick/integrate1.hpp | 5 ++--- Bembel/src/DuffyTrick/integrate2.hpp | 5 ++--- Bembel/src/DuffyTrick/integrate3.hpp | 5 ++--- Bembel/src/DuffyTrick/integrate4.hpp | 5 ++--- Bembel/src/Geometry/SurfacePoint.hpp | 9 +++++++++ tests/DuffyTrick/test_integrate0.hpp | 2 +- tests/DuffyTrick/test_integrate1.hpp | 2 +- tests/DuffyTrick/test_integrate2.hpp | 2 +- tests/DuffyTrick/test_integrate3.hpp | 2 +- tests/DuffyTrick/test_integrate4.hpp | 2 +- 12 files changed, 27 insertions(+), 23 deletions(-) diff --git a/Bembel/src/DuffyTrick/farFieldQuadratureNodes.hpp b/Bembel/src/DuffyTrick/farFieldQuadratureNodes.hpp index 41c05d379..893781ee5 100644 --- a/Bembel/src/DuffyTrick/farFieldQuadratureNodes.hpp +++ b/Bembel/src/DuffyTrick/farFieldQuadratureNodes.hpp @@ -20,9 +20,9 @@ namespace DuffyTrick { * is qNodes.col(k) = [xi, w, Chi(xi); dsChi(xi); dtChi(xi)]\in\Rbb^12 **/ template -std::vector> computeFfieldQnodes(const T &super_space, - const Cubature &Q) { - std::vector> ffield_qnodes; +std::vector computeFfieldQnodes(const T &super_space, + const Cubature &Q) { + std::vector ffield_qnodes; int next = 0; // assume isotropic mesh width h! double h = (super_space.get_mesh().get_element_tree().cpbegin())->get_h(); diff --git a/Bembel/src/DuffyTrick/integrate0.hpp b/Bembel/src/DuffyTrick/integrate0.hpp index f095038d3..f5967aeff 100644 --- a/Bembel/src/DuffyTrick/integrate0.hpp +++ b/Bembel/src/DuffyTrick/integrate0.hpp @@ -22,9 +22,8 @@ namespace DuffyTrick { template void integrate0(const LinearOperatorBase &LinOp, const T &super_space, const ElementTreeNode &e1, int rot1, const ElementTreeNode &e2, - int rot2, const std::vector &ffield_qnodes1, - const std::vector &ffield_qnodes2, - const Cubature &Q, + int rot2, const ElementSurfacePoints &ffield_qnodes1, + const ElementSurfacePoints &ffield_qnodes2, const Cubature &Q, Eigen::Matrix::Scalar, Eigen::Dynamic, Eigen::Dynamic> *intval) { intval->setZero(); diff --git a/Bembel/src/DuffyTrick/integrate1.hpp b/Bembel/src/DuffyTrick/integrate1.hpp index fa37711f4..77a888240 100644 --- a/Bembel/src/DuffyTrick/integrate1.hpp +++ b/Bembel/src/DuffyTrick/integrate1.hpp @@ -24,9 +24,8 @@ namespace DuffyTrick { template void integrate1(const LinearOperatorBase &LinOp, const T &super_space, const ElementTreeNode &e1, int rot1, const ElementTreeNode &e2, - int rot2, const std::vector &ffield_qnodes1, - const std::vector &ffield_qnodes2, - const Cubature &Q, + int rot2, const ElementSurfacePoints &ffield_qnodes1, + const ElementSurfacePoints &ffield_qnodes2, const Cubature &Q, Eigen::Matrix::Scalar, Eigen::Dynamic, Eigen::Dynamic> *intval) { double h = e1.get_h(); diff --git a/Bembel/src/DuffyTrick/integrate2.hpp b/Bembel/src/DuffyTrick/integrate2.hpp index accf64bcb..8dc469b34 100644 --- a/Bembel/src/DuffyTrick/integrate2.hpp +++ b/Bembel/src/DuffyTrick/integrate2.hpp @@ -27,9 +27,8 @@ namespace DuffyTrick { template void integrate2(const LinearOperatorBase &LinOp, const T &super_space, const ElementTreeNode &e1, int rot1, const ElementTreeNode &e2, - int rot2, const std::vector &ffield_qnodes1, - const std::vector &ffield_qnodes2, - const Cubature &Q, + int rot2, const ElementSurfacePoints &ffield_qnodes1, + const ElementSurfacePoints &ffield_qnodes2, const Cubature &Q, Eigen::Matrix::Scalar, Eigen::Dynamic, Eigen::Dynamic> *intval) { intval->setZero(); diff --git a/Bembel/src/DuffyTrick/integrate3.hpp b/Bembel/src/DuffyTrick/integrate3.hpp index cd0cc89d0..8df2f4669 100644 --- a/Bembel/src/DuffyTrick/integrate3.hpp +++ b/Bembel/src/DuffyTrick/integrate3.hpp @@ -27,9 +27,8 @@ namespace DuffyTrick { template void integrate3(const LinearOperatorBase &LinOp, const T &super_space, const ElementTreeNode &e1, int rot1, const ElementTreeNode &e2, - int rot2, const std::vector &ffield_qnodes1, - const std::vector &ffield_qnodes2, - const Cubature &Q, + int rot2, const ElementSurfacePoints &ffield_qnodes1, + const ElementSurfacePoints &ffield_qnodes2, const Cubature &Q, Eigen::Matrix::Scalar, Eigen::Dynamic, Eigen::Dynamic> *intval) { intval->setZero(); diff --git a/Bembel/src/DuffyTrick/integrate4.hpp b/Bembel/src/DuffyTrick/integrate4.hpp index 603fef5d1..b3ca5f3f3 100644 --- a/Bembel/src/DuffyTrick/integrate4.hpp +++ b/Bembel/src/DuffyTrick/integrate4.hpp @@ -27,9 +27,8 @@ namespace DuffyTrick { template void integrate4(const LinearOperatorBase &LinOp, const T &super_space, const ElementTreeNode &e1, int rot1, const ElementTreeNode &e2, - int rot2, const std::vector &ffield_qnodes1, - const std::vector &ffield_qnodes2, - const Cubature &Q, + int rot2, const ElementSurfacePoints &ffield_qnodes1, + const ElementSurfacePoints &ffield_qnodes2, const Cubature &Q, Eigen::Matrix::Scalar, Eigen::Dynamic, Eigen::Dynamic> *intval) { intval->setZero(); diff --git a/Bembel/src/Geometry/SurfacePoint.hpp b/Bembel/src/Geometry/SurfacePoint.hpp index 2322be945..18492a640 100644 --- a/Bembel/src/Geometry/SurfacePoint.hpp +++ b/Bembel/src/Geometry/SurfacePoint.hpp @@ -11,6 +11,7 @@ // #ifndef BEMBEL_SRC_GEOMETRY_SURFACEPOINT_HPP_ #define BEMBEL_SRC_GEOMETRY_SURFACEPOINT_HPP_ +#include /** * \ingroup Geometry * \brief typedef of SurfacePoint @@ -40,4 +41,12 @@ **/ typedef Eigen::Matrix SurfacePoint; +/** + * \ingroup Geometry + * \brief typedef std::vector with aligned allocator of Eigen for + * compatibility with older compilers. + */ +typedef std::vector> + ElementSurfacePoints; + #endif // BEMBEL_SRC_GEOMETRY_SURFACEPOINT_HPP_ diff --git a/tests/DuffyTrick/test_integrate0.hpp b/tests/DuffyTrick/test_integrate0.hpp index 09a617bc8..a1de07961 100644 --- a/tests/DuffyTrick/test_integrate0.hpp +++ b/tests/DuffyTrick/test_integrate0.hpp @@ -20,7 +20,7 @@ bool test_integrate0(const Bembel::AnsatzSpace &ansatz_space, const Bembel::LinearOperatorBase &linOp) { Bembel::GaussSquare GS; auto Q = GS[maxqdeg]; - auto ffield_qnodes = + std::vector ffield_qnodes = Bembel::DuffyTrick::computeFfieldQnodes(ansatz_space.get_superspace(), Q); Eigen::MatrixXd intval; Eigen::MatrixXd axis; diff --git a/tests/DuffyTrick/test_integrate1.hpp b/tests/DuffyTrick/test_integrate1.hpp index 41ab4c339..362e4d45d 100644 --- a/tests/DuffyTrick/test_integrate1.hpp +++ b/tests/DuffyTrick/test_integrate1.hpp @@ -20,7 +20,7 @@ bool test_integrate1(const Bembel::AnsatzSpace &ansatz_space, const Bembel::LinearOperatorBase &linOp) { Bembel::GaussSquare GS; auto Q = GS[maxqdeg]; - std::vector ffield_qnodes; + ElementSurfacePoints ffield_qnodes; Eigen::MatrixXd intval; Eigen::MatrixXd axis; intval.resize(1, 1); diff --git a/tests/DuffyTrick/test_integrate2.hpp b/tests/DuffyTrick/test_integrate2.hpp index a4ecbd531..eeada5f82 100644 --- a/tests/DuffyTrick/test_integrate2.hpp +++ b/tests/DuffyTrick/test_integrate2.hpp @@ -21,7 +21,7 @@ bool test_integrate2(const Bembel::AnsatzSpace &ansatz_space, Bembel::GaussSquare GS; auto Q = GS[maxqdeg]; - std::vector ffield_qnodes; + ElementSurfacePoints ffield_qnodes; Eigen::MatrixXd intval; Eigen::MatrixXd axis; intval.resize(1, 1); diff --git a/tests/DuffyTrick/test_integrate3.hpp b/tests/DuffyTrick/test_integrate3.hpp index 98e8e5f05..8df2511d9 100644 --- a/tests/DuffyTrick/test_integrate3.hpp +++ b/tests/DuffyTrick/test_integrate3.hpp @@ -21,7 +21,7 @@ bool test_integrate3(const Bembel::AnsatzSpace &ansatz_space, Bembel::GaussSquare GS; auto Q = GS[maxqdeg]; - std::vector ffield_qnodes; + ElementSurfacePoints ffield_qnodes; Eigen::MatrixXd intval; Eigen::MatrixXd axis; intval.resize(1, 1); diff --git a/tests/DuffyTrick/test_integrate4.hpp b/tests/DuffyTrick/test_integrate4.hpp index 2b14f9c74..2d1b6dcd2 100644 --- a/tests/DuffyTrick/test_integrate4.hpp +++ b/tests/DuffyTrick/test_integrate4.hpp @@ -20,7 +20,7 @@ bool test_integrate4(const Bembel::AnsatzSpace &ansatz_space, const Bembel::LinearOperatorBase &linOp) { Bembel::GaussSquare GS; auto Q = GS[maxqdeg]; - std::vector ffield_qnodes; + ElementSurfacePoints ffield_qnodes; Eigen::MatrixXd intval; Eigen::MatrixXd axis; intval.resize(1, 1); From 00183498640ecbb231c377c8785d0244cdd824d2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20D=C3=B6lz?= Date: Tue, 7 Feb 2023 17:37:19 +0100 Subject: [PATCH 06/20] Replaced vector representation of SurfacePoint by a wrapper class --- .../src/DuffyTrick/evaluateBilinearForm.hpp | 4 +- Bembel/src/Geometry/SurfacePoint.hpp | 89 +++++++++++++++++-- Bembel/src/H2Matrix/H2Matrix.hpp | 2 +- tests/test_SurfacePointUpdate.cpp | 2 +- 4 files changed, 84 insertions(+), 13 deletions(-) diff --git a/Bembel/src/DuffyTrick/evaluateBilinearForm.hpp b/Bembel/src/DuffyTrick/evaluateBilinearForm.hpp index e1d7f55cb..5b3dd6f3e 100644 --- a/Bembel/src/DuffyTrick/evaluateBilinearForm.hpp +++ b/Bembel/src/DuffyTrick/evaluateBilinearForm.hpp @@ -22,8 +22,8 @@ template void evaluateBilinearForm( const LinearOperatorBase& linOp, const T& super_space, const ElementTreeNode& e1, const ElementTreeNode& e2, - const CubatureVector& GS, const std::vector& ffield_qnodes1, - const std::vector& ffield_qnodes2, + const CubatureVector& GS, const ElementSurfacePoints& ffield_qnodes1, + const ElementSurfacePoints& ffield_qnodes2, Eigen::Matrix::Scalar, Eigen::Dynamic, Eigen::Dynamic>* intval) { ////////////////////////////////////////////////////////////////////////////// diff --git a/Bembel/src/Geometry/SurfacePoint.hpp b/Bembel/src/Geometry/SurfacePoint.hpp index 18492a640..52b8389fb 100644 --- a/Bembel/src/Geometry/SurfacePoint.hpp +++ b/Bembel/src/Geometry/SurfacePoint.hpp @@ -11,10 +11,10 @@ // #ifndef BEMBEL_SRC_GEOMETRY_SURFACEPOINT_HPP_ #define BEMBEL_SRC_GEOMETRY_SURFACEPOINT_HPP_ -#include +// #include /** * \ingroup Geometry - * \brief typedef of SurfacePoint + * \brief Wrapper class for legacy SurfacePoint represented as an Eigen::Vector * * This typedef is essential for any evaluation of a bilinear form. It provides * all required geometry information, stored as follows: @@ -39,14 +39,85 @@ * updateSurdacePoint method is specialized and should be used, since it avoids * redundant work. **/ -typedef Eigen::Matrix SurfacePoint; +class SurfacePoint { + public: + /** + * \brief Supports legacy code. + */ + double& operator()(const int j) { return (data_(j)); } + /** + * \brief Supports legacy code. + */ + const double& operator()(const int j) const { return (data_(j)); } + /** + * \brief Supports legacy code. + */ + template + Eigen::Matrix segment(const int l) { + return data_.segment(l); + } + /** + * \brief Supports legacy code. + */ + template + const Eigen::Matrix segment(const int l) const { + return data_.segment(l); + } + /** + * \brief Supports legacy code. + */ + Eigen::Matrix segment(const int l, const int n) { + return data_.segment(l, n); + } + /** + * \brief Supports legacy code. + */ + const Eigen::Matrix segment(const int l, + const int n) const { + return data_.segment(l, n); + } + /** + * \brief Supports legacy code. + */ + Eigen::Matrix head(const int n) { + return data_.head(n); + } + /** + * \brief Supports legacy code. + */ + const Eigen::Matrix head(const int n) const { + return data_.head(n); + } + /** + * \brief Supports legacy code. + */ + Eigen::Matrix tail(const int n) { + return data_.tail(n); + } + /** + * \brief Supports legacy code. + */ + const Eigen::Matrix tail(const int n) const { + return data_.tail(n); + } + /** + * \brief Supports legacy code. + */ + Eigen::Matrix get_data() { return data_; } + /** + * \brief Supports legacy code. + */ + const Eigen::Matrix get_data() const { + return data_; + } + + // EIGEN_MAKE_ALIGNED_OPERATOR_NEW + + private: + Eigen::Matrix data_; +}; -/** - * \ingroup Geometry - * \brief typedef std::vector with aligned allocator of Eigen for - * compatibility with older compilers. - */ typedef std::vector> ElementSurfacePoints; -#endif // BEMBEL_SRC_GEOMETRY_SURFACEPOINT_HPP_ +#endif // BEMBEL_SRC_GEOMETRY_SURFACEPOINT_HPP_ \ No newline at end of file diff --git a/Bembel/src/H2Matrix/H2Matrix.hpp b/Bembel/src/H2Matrix/H2Matrix.hpp index 80288b3bb..590a91759 100644 --- a/Bembel/src/H2Matrix/H2Matrix.hpp +++ b/Bembel/src/H2Matrix/H2Matrix.hpp @@ -134,7 +134,7 @@ class H2Matrix : public EigenBase> { Bembel::GaussSquare GS; auto super_space = ansatz_space.get_superspace(); auto ffield_deg = linOp.get_FarfieldQuadratureDegree(polynomial_degree); - auto ffield_qnodes = + std::vector ffield_qnodes = Bembel::DuffyTrick::computeFfieldQnodes(super_space, GS[ffield_deg]); const int NumberOfFMMComponents = Bembel::LinearOperatorTraits::NumberOfFMMComponents; diff --git a/tests/test_SurfacePointUpdate.cpp b/tests/test_SurfacePointUpdate.cpp index 166ce7ff7..c6cacba40 100644 --- a/tests/test_SurfacePointUpdate.cpp +++ b/tests/test_SurfacePointUpdate.cpp @@ -38,7 +38,7 @@ int main() { geometry.get_geometry()[0].updateSurfacePoint(&srf_pt, pt, 3.1415, pt); - if ((srf_pt - srf_pt_ref).norm() > + if ((srf_pt.get_data() - srf_pt_ref.get_data()).norm() > Test::Constants::test_tolerance_geometry) return 1; } From 96523f42ef38dcb235b5b69880e4a67d5a9a2412 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20D=C3=B6lz?= Date: Wed, 8 Feb 2023 10:20:33 +0100 Subject: [PATCH 07/20] new surface point works now, although it still seems to be a bit slow --- .../src/AnsatzSpace/FunctionEvaluatorEval.hpp | 25 ++-- Bembel/src/AnsatzSpace/SuperSpace.hpp | 15 +-- Bembel/src/Geometry/Patch.hpp | 28 ++--- Bembel/src/Geometry/SurfacePoint.hpp | 111 +++++++----------- Bembel/src/Helmholtz/SingleLayerOperator.hpp | 50 ++------ Bembel/src/Helmholtz/SingleLayerPotential.hpp | 22 +--- Bembel/src/Identity/IdentityOperatorBase.hpp | 21 +--- Bembel/src/Laplace/SingleLayerOperator.hpp | 50 ++------ Bembel/src/Laplace/SingleLayerPotential.hpp | 22 +--- .../LaplaceBeltramiOperatorBase.hpp | 8 +- Bembel/src/LinearForm/DirichletTrace.hpp | 24 +--- .../src/LinearForm/RotatedTangentialTrace.hpp | 38 ++---- .../LinearOperator/Dummy/DummyOperator.hpp | 4 +- Bembel/src/Maxwell/SingleLayerOperator.hpp | 52 ++------ Bembel/src/Maxwell/SingleLayerPotential.hpp | 18 +-- Bembel/src/util/surfaceL2error.hpp | 23 +--- tests/test_SurfacePointUpdate.cpp | 29 ++--- 17 files changed, 141 insertions(+), 399 deletions(-) diff --git a/Bembel/src/AnsatzSpace/FunctionEvaluatorEval.hpp b/Bembel/src/AnsatzSpace/FunctionEvaluatorEval.hpp index a32bf2122..2d81f85c4 100644 --- a/Bembel/src/AnsatzSpace/FunctionEvaluatorEval.hpp +++ b/Bembel/src/AnsatzSpace/FunctionEvaluatorEval.hpp @@ -29,8 +29,7 @@ struct FunctionEvaluatorEval { Scalar, Eigen::Dynamic, getFunctionSpaceVectorDimension()> &coeff) const { - auto s = p.segment<2>(0); - return coeff.transpose() * super_space.basis(s) / element.get_h(); + return coeff.transpose() * super_space.basis(p.get_xi()) / element.get_h(); } }; @@ -47,15 +46,13 @@ struct FunctionEvaluatorEval { Scalar, Eigen::Dynamic, getFunctionSpaceVectorDimension()> &coeff) const { - auto s = p.segment<2>(0); - auto h = element.get_h(); - auto x_f_dx = p.segment<3>(6); - auto x_f_dy = p.segment<3>(9); + double h = element.get_h(); Eigen::Matrix::Scalar, Eigen::Dynamic, 1> - tangential_coefficients = coeff.transpose() * super_space.basis(s); - return (x_f_dx * tangential_coefficients(0) + - x_f_dy * tangential_coefficients(1)) / + tangential_coefficients = + coeff.transpose() * super_space.basis(p.get_xi()); + return (p.get_f_dx() * tangential_coefficients(0) + + p.get_f_dy() * tangential_coefficients(1)) / h; } @@ -67,14 +64,13 @@ struct FunctionEvaluatorEval { Scalar, Eigen::Dynamic, getFunctionSpaceVectorDimension()> &coeff) const { - auto s = p.segment<2>(0); - auto h = element.get_h(); + double h = element.get_h(); Eigen::Matrix::Scalar, Eigen::Dynamic, 1> - phiPhiVec_dx = super_space.basisDx(s); + phiPhiVec_dx = super_space.basisDx(p.get_xi()); Eigen::Matrix::Scalar, Eigen::Dynamic, 1> - phiPhiVec_dy = super_space.basisDy(s); + phiPhiVec_dy = super_space.basisDy(p.get_xi()); return (phiPhiVec_dx.dot(coeff.col(0)) + phiPhiVec_dy.dot(coeff.col(1))) / h / h; } @@ -93,8 +89,7 @@ struct FunctionEvaluatorEval { Scalar, Eigen::Dynamic, getFunctionSpaceVectorDimension()> &coeff) const { - auto s = p.segment<2>(0); - return coeff.transpose() * super_space.basis(s) / element.get_h(); + return coeff.transpose() * super_space.basis(p.get_xi()) / element.get_h(); } }; } // namespace Bembel diff --git a/Bembel/src/AnsatzSpace/SuperSpace.hpp b/Bembel/src/AnsatzSpace/SuperSpace.hpp index 2845883dd..a0be0ed00 100644 --- a/Bembel/src/AnsatzSpace/SuperSpace.hpp +++ b/Bembel/src/AnsatzSpace/SuperSpace.hpp @@ -143,19 +143,16 @@ struct SuperSpace { void addScaledSurfaceCurlInteraction( Eigen::Matrix* intval, Scalar w, const SurfacePoint& p1, const SurfacePoint& p2) const { - // surface measures - double kappa1 = p1.segment<3>(6).cross(p1.segment<3>(9)).norm(); - double kappa2 = p2.segment<3>(6).cross(p2.segment<3>(9)).norm(); // compute basis functions's surface curl. Each column of s_curl is a basis // function's surface curl at point s. Eigen::MatrixXd s_curl(3, polynomial_degree_plus_one_squared); - s_curl = (1.0 / kappa1) * - (-p1.segment<3>(6) * basisDy(p1.segment<2>(0)).transpose() + - p1.segment<3>(9) * basisDx(p1.segment<2>(0)).transpose()); + s_curl = (1.0 / p1.get_surface_measure()) * + (-p1.get_f_dx() * basisDy(p1.get_xi()).transpose() + + p1.get_f_dy() * basisDx(p1.get_xi()).transpose()); Eigen::MatrixXd t_curl(3, polynomial_degree_plus_one_squared); - t_curl = (1.0 / kappa2) * - (-p2.segment<3>(6) * basisDy(p2.segment<2>(0)).transpose() + - p2.segment<3>(9) * basisDx(p2.segment<2>(0)).transpose()); + t_curl = (1.0 / p2.get_surface_measure()) * + (-p2.get_f_dx() * basisDy(p2.get_xi()).transpose() + + p2.get_f_dy() * basisDx(p2.get_xi()).transpose()); // inner product of surface curls of any two basis functions for (int j = 0; j < polynomial_degree_plus_one_squared; ++j) for (int i = 0; i < polynomial_degree_plus_one_squared; ++i) diff --git a/Bembel/src/Geometry/Patch.hpp b/Bembel/src/Geometry/Patch.hpp index ffd44df06..0f1676f4a 100644 --- a/Bembel/src/Geometry/Patch.hpp +++ b/Bembel/src/Geometry/Patch.hpp @@ -209,9 +209,8 @@ class Patch { // This is a combination of eval und evalJacobian, to avoid duplication of // work. See SurfacePoint.hpp - void updateSurfacePoint(SurfacePoint *srf_pt, - const Eigen::Vector2d &ref_pt, double w, - const Eigen::Vector2d &xi) const { + void updateSurfacePoint(SurfacePoint *srf_pt, const Eigen::Vector2d &ref_pt, + double w, const Eigen::Vector2d &xi) const { const int x_location = Spl::FindLocationInKnotVector(ref_pt(0), unique_knots_x_); const int y_location = @@ -266,18 +265,17 @@ class Patch { const double bot = 1. / tmp[3]; const double botsqr = bot * bot; - (*srf_pt)(0) = xi(0); - (*srf_pt)(1) = xi(1); - (*srf_pt)(2) = w; - (*srf_pt)(3) = tmp[0] * bot; - (*srf_pt)(4) = tmp[1] * bot; - (*srf_pt)(5) = tmp[2] * bot; - (*srf_pt)(6) = (tmpDx[0] * tmp[3] - tmp[0] * tmpDx[3]) * botsqr; - (*srf_pt)(7) = (tmpDx[1] * tmp[3] - tmp[1] * tmpDx[3]) * botsqr; - (*srf_pt)(8) = (tmpDx[2] * tmp[3] - tmp[2] * tmpDx[3]) * botsqr; - (*srf_pt)(9) = (tmpDy[0] * tmp[3] - tmp[0] * tmpDy[3]) * botsqr; - (*srf_pt)(10) = (tmpDy[1] * tmp[3] - tmp[1] * tmpDy[3]) * botsqr; - (*srf_pt)(11) = (tmpDy[2] * tmp[3] - tmp[2] * tmpDy[3]) * botsqr; + Eigen::Vector3d tmpvec(tmp[0], tmp[1], tmp[2]); + Eigen::Vector3d tmpDxvec(tmpDx[0], tmpDx[1], tmpDx[2]); + Eigen::Vector3d tmpDyvec(tmpDy[0], tmpDy[1], tmpDy[2]); + + srf_pt->set_xi(xi); + srf_pt->set_w(w); + srf_pt->set_f(bot * tmpvec); + srf_pt->set_jacobian(botsqr * (Eigen::Matrix() + << tmpDxvec * tmp[3] - tmpvec * tmpDx[3], + tmpDyvec * tmp[3] - tmpvec * tmpDy[3]) + .finished()); delete[] buffer; return; } diff --git a/Bembel/src/Geometry/SurfacePoint.hpp b/Bembel/src/Geometry/SurfacePoint.hpp index 52b8389fb..a79916dd5 100644 --- a/Bembel/src/Geometry/SurfacePoint.hpp +++ b/Bembel/src/Geometry/SurfacePoint.hpp @@ -41,80 +41,55 @@ **/ class SurfacePoint { public: - /** - * \brief Supports legacy code. - */ - double& operator()(const int j) { return (data_(j)); } - /** - * \brief Supports legacy code. - */ - const double& operator()(const int j) const { return (data_(j)); } - /** - * \brief Supports legacy code. - */ - template - Eigen::Matrix segment(const int l) { - return data_.segment(l); - } - /** - * \brief Supports legacy code. - */ - template - const Eigen::Matrix segment(const int l) const { - return data_.segment(l); - } - /** - * \brief Supports legacy code. - */ - Eigen::Matrix segment(const int l, const int n) { - return data_.segment(l, n); - } - /** - * \brief Supports legacy code. - */ - const Eigen::Matrix segment(const int l, - const int n) const { - return data_.segment(l, n); - } - /** - * \brief Supports legacy code. - */ - Eigen::Matrix head(const int n) { - return data_.head(n); - } - /** - * \brief Supports legacy code. - */ - const Eigen::Matrix head(const int n) const { - return data_.head(n); - } - /** - * \brief Supports legacy code. - */ - Eigen::Matrix tail(const int n) { - return data_.tail(n); + // quadrature point + Eigen::Vector2d get_xi() { return xi_; } + const Eigen::Vector2d get_xi() const { return xi_; } + + // quadrature weight + double get_w() { return w_; } + const double get_w() const { return w_; } + + // geometry parametrization + Eigen::Vector3d get_f() { return f_; } + const Eigen::Vector3d get_f() const { return f_; } + + // geometry parametrization first derivatives + Eigen::Vector3d get_f_dx() { return jacobian_.col(0); } + const Eigen::Vector3d get_f_dx() const { return jacobian_.col(0); } + Eigen::Vector3d get_f_dy() { return jacobian_.col(1); } + const Eigen::Vector3d get_f_dy() const { return jacobian_.col(1); } + Eigen::Matrix get_jacobian() { return jacobian_; } + const Eigen::Matrix get_jacobian() const { return jacobian_; } + + // normal vectors + Eigen::Vector3d get_normal() { + return jacobian_.col(0).cross(jacobian_.col(1)); } - /** - * \brief Supports legacy code. - */ - const Eigen::Matrix tail(const int n) const { - return data_.tail(n); + const Eigen::Vector3d get_normal() const { + return jacobian_.col(0).cross(jacobian_.col(1)); } - /** - * \brief Supports legacy code. - */ - Eigen::Matrix get_data() { return data_; } - /** - * \brief Supports legacy code. - */ - const Eigen::Matrix get_data() const { - return data_; + Eigen::Vector3d get_unit_normal() { return get_normal().normalized(); } + const Eigen::Vector3d get_unit_normal() const { + return get_normal().normalized(); } - // EIGEN_MAKE_ALIGNED_OPERATOR_NEW + // surface measure + double get_surface_measure() { return get_normal().norm(); } + const double get_surface_measure() const { return get_normal().norm(); } + + // setter + void set_xi(const Eigen::Vector2d &xi) { xi_ = xi; } + void set_w(const double w) { w_ = w; } + void set_f(const Eigen::Vector3d &f) { f_ = f; } + void set_jacobian(const Eigen::Matrix &jacobian) { + jacobian_ = jacobian; + } private: - Eigen::Matrix data_; + Eigen::Vector2d xi_; + double w_; + Eigen::Vector3d f_; + Eigen::Matrix jacobian_; }; typedef std::vector> diff --git a/Bembel/src/Helmholtz/SingleLayerOperator.hpp b/Bembel/src/Helmholtz/SingleLayerOperator.hpp index 466866174..1dee1e590 100644 --- a/Bembel/src/Helmholtz/SingleLayerOperator.hpp +++ b/Bembel/src/Helmholtz/SingleLayerOperator.hpp @@ -43,64 +43,28 @@ class HelmholtzSingleLayerOperator Eigen::Matrix< typename LinearOperatorTraits::Scalar, Eigen::Dynamic, Eigen::Dynamic> *intval) const { - auto polynomial_degree = super_space.get_polynomial_degree(); - auto polynomial_degree_plus_one_squared = - (polynomial_degree + 1) * (polynomial_degree + 1); - - // get evaluation points on unit square - auto s = p1.segment<2>(0); - auto t = p2.segment<2>(0); - - // get quadrature weights - auto ws = p1(2); - auto wt = p2(2); - - // get points on geometry and tangential derivatives - auto x_f = p1.segment<3>(3); - auto x_f_dx = p1.segment<3>(6); - auto x_f_dy = p1.segment<3>(9); - auto y_f = p2.segment<3>(3); - auto y_f_dx = p2.segment<3>(6); - auto y_f_dy = p2.segment<3>(9); - - // compute surface measures from tangential derivatives - auto x_kappa = x_f_dx.cross(x_f_dy).norm(); - auto y_kappa = y_f_dx.cross(y_f_dy).norm(); - // integrand without basis functions - auto integrand = evaluateKernel(x_f, y_f) * x_kappa * y_kappa * ws * wt; + std::complex integrand = + evaluateKernel(p1.get_f(), p2.get_f()) * p1.get_surface_measure() * + p2.get_surface_measure() * p1.get_w() * p2.get_w(); // multiply basis functions with integrand and add to intval, this is an // efficient implementation of // (*intval) += super_space.BasisInteraction(s, t) * evaluateKernel(x_f, // y_f) // * x_kappa * y_kappa * ws * wt; - super_space.addScaledBasisInteraction(intval, integrand, s, t); + super_space.addScaledBasisInteraction(intval, integrand, p1.get_xi(), + p2.get_xi()); return; } Eigen::Matrix, 1, 1> evaluateFMMInterpolation_impl( const SurfacePoint &p1, const SurfacePoint &p2) const { - // get evaluation points on unit square - auto s = p1.segment<2>(0); - auto t = p2.segment<2>(0); - - // get points on geometry and tangential derivatives - auto x_f = p1.segment<3>(3); - auto x_f_dx = p1.segment<3>(6); - auto x_f_dy = p1.segment<3>(9); - auto y_f = p2.segment<3>(3); - auto y_f_dx = p2.segment<3>(6); - auto y_f_dy = p2.segment<3>(9); - - // compute surface measures from tangential derivatives - auto x_kappa = x_f_dx.cross(x_f_dy).norm(); - auto y_kappa = y_f_dx.cross(y_f_dy).norm(); - // interpolation Eigen::Matrix, 1, 1> intval; - intval(0) = evaluateKernel(x_f, y_f) * x_kappa * y_kappa; + intval(0) = evaluateKernel(p1.get_f(), p2.get_f()) * + p1.get_surface_measure() * p2.get_surface_measure(); return intval; } diff --git a/Bembel/src/Helmholtz/SingleLayerPotential.hpp b/Bembel/src/Helmholtz/SingleLayerPotential.hpp index 654570b10..e199464b1 100644 --- a/Bembel/src/Helmholtz/SingleLayerPotential.hpp +++ b/Bembel/src/Helmholtz/SingleLayerPotential.hpp @@ -41,29 +41,13 @@ class HelmholtzSingleLayerPotential const ElementTreeNode &element, const Eigen::Vector3d &point, const SurfacePoint &p) const { - // get evaluation points on unit square - auto s = p.segment<2>(0); - - // get quadrature weights - auto ws = p(2); - - // get points on geometry and tangential derivatives - auto x_f = p.segment<3>(3); - auto x_f_dx = p.segment<3>(6); - auto x_f_dy = p.segment<3>(9); - - // compute surface measures from tangential derivatives - auto x_kappa = x_f_dx.cross(x_f_dy).norm(); - // evaluate kernel - auto kernel = evaluateKernel(point, x_f); - + auto kernel = evaluateKernel(point, p.get_f()); // assemble Galerkin solution auto cauchy_value = fun_ev.evaluate(element, p); - // integrand without basis functions - auto integrand = kernel * cauchy_value * x_kappa * ws; - + auto integrand = + kernel * cauchy_value * p.get_surface_measure() * p.get_w(); return integrand; } diff --git a/Bembel/src/Identity/IdentityOperatorBase.hpp b/Bembel/src/Identity/IdentityOperatorBase.hpp index ab3b1d930..c486cc038 100644 --- a/Bembel/src/Identity/IdentityOperatorBase.hpp +++ b/Bembel/src/Identity/IdentityOperatorBase.hpp @@ -28,25 +28,10 @@ class IdentityOperatorBase : public LocalOperatorBase { void evaluateIntegrand_impl(const T &super_space, const SurfacePoint &p1, const SurfacePoint &p2, Eigen::MatrixXd *intval) const { - // get evaluation points on unit square - const auto s = p1.segment<2>(0); - - // get quadrature weights - const auto ws = p1(2); - - // get points on geometry and tangential derivatives - const auto &x_f = p1.segment<3>(3); - const auto &x_f_dx = p1.segment<3>(6); - const auto &x_f_dy = p1.segment<3>(9); - - // compute surface measures from tangential derivatives - const auto x_kappa = x_f_dx.cross(x_f_dy).norm(); - // integrand without basis functions - const auto integrand = x_kappa * ws; - - super_space.addScaledBasisInteraction(intval, integrand, s, s); - + const auto integrand = p1.get_surface_measure() * p1.get_w(); + super_space.addScaledBasisInteraction(intval, integrand, p1.get_xi(), + p1.get_xi()); return; } }; diff --git a/Bembel/src/Laplace/SingleLayerOperator.hpp b/Bembel/src/Laplace/SingleLayerOperator.hpp index 56dd1a2c2..96eda73cc 100644 --- a/Bembel/src/Laplace/SingleLayerOperator.hpp +++ b/Bembel/src/Laplace/SingleLayerOperator.hpp @@ -43,64 +43,28 @@ class LaplaceSingleLayerOperator Eigen::Matrix< typename LinearOperatorTraits::Scalar, Eigen::Dynamic, Eigen::Dynamic> *intval) const { - auto polynomial_degree = super_space.get_polynomial_degree(); - auto polynomial_degree_plus_one_squared = - (polynomial_degree + 1) * (polynomial_degree + 1); - - // get evaluation points on unit square - auto s = p1.segment<2>(0); - auto t = p2.segment<2>(0); - - // get quadrature weights - auto ws = p1(2); - auto wt = p2(2); - - // get points on geometry and tangential derivatives - auto x_f = p1.segment<3>(3); - auto x_f_dx = p1.segment<3>(6); - auto x_f_dy = p1.segment<3>(9); - auto y_f = p2.segment<3>(3); - auto y_f_dx = p2.segment<3>(6); - auto y_f_dy = p2.segment<3>(9); - - // compute surface measures from tangential derivatives - auto x_kappa = x_f_dx.cross(x_f_dy).norm(); - auto y_kappa = y_f_dx.cross(y_f_dy).norm(); - // integrand without basis functions - auto integrand = evaluateKernel(x_f, y_f) * x_kappa * y_kappa * ws * wt; + double integrand = evaluateKernel(p1.get_f(), p2.get_f()) * + p1.get_surface_measure() * p2.get_surface_measure() * + p1.get_w() * p2.get_w(); // multiply basis functions with integrand and add to intval, this is an // efficient implementation of // (*intval) += super_space.basisInteraction(s, t) * evaluateKernel(x_f, // y_f) // * x_kappa * y_kappa * ws * wt; - super_space.addScaledBasisInteraction(intval, integrand, s, t); + super_space.addScaledBasisInteraction(intval, integrand, p1.get_xi(), + p2.get_xi()); return; } Eigen::Matrix evaluateFMMInterpolation_impl( const SurfacePoint &p1, const SurfacePoint &p2) const { - // get evaluation points on unit square - auto s = p1.segment<2>(0); - auto t = p2.segment<2>(0); - - // get points on geometry and tangential derivatives - auto x_f = p1.segment<3>(3); - auto x_f_dx = p1.segment<3>(6); - auto x_f_dy = p1.segment<3>(9); - auto y_f = p2.segment<3>(3); - auto y_f_dx = p2.segment<3>(6); - auto y_f_dy = p2.segment<3>(9); - - // compute surface measures from tangential derivatives - auto x_kappa = x_f_dx.cross(x_f_dy).norm(); - auto y_kappa = y_f_dx.cross(y_f_dy).norm(); - // interpolation Eigen::Matrix intval; - intval(0) = evaluateKernel(x_f, y_f) * x_kappa * y_kappa; + intval(0) = evaluateKernel(p1.get_f(), p2.get_f()) * + p1.get_surface_measure() * p2.get_surface_measure(); return intval; } diff --git a/Bembel/src/Laplace/SingleLayerPotential.hpp b/Bembel/src/Laplace/SingleLayerPotential.hpp index 6b23d553f..ec4f27e43 100644 --- a/Bembel/src/Laplace/SingleLayerPotential.hpp +++ b/Bembel/src/Laplace/SingleLayerPotential.hpp @@ -41,29 +41,13 @@ class LaplaceSingleLayerPotential const ElementTreeNode &element, const Eigen::Vector3d &point, const SurfacePoint &p) const { - // get evaluation points on unit square - auto s = p.segment<2>(0); - - // get quadrature weights - auto ws = p(2); - - // get points on geometry and tangential derivatives - auto x_f = p.segment<3>(3); - auto x_f_dx = p.segment<3>(6); - auto x_f_dy = p.segment<3>(9); - - // compute surface measures from tangential derivatives - auto x_kappa = x_f_dx.cross(x_f_dy).norm(); - // evaluate kernel - auto kernel = evaluateKernel(point, x_f); - + auto kernel = evaluateKernel(point, p.get_f()); // assemble Galerkin solution auto cauchy_value = fun_ev.evaluate(element, p); - // integrand without basis functions - auto integrand = kernel * cauchy_value * x_kappa * ws; - + auto integrand = + kernel * cauchy_value * p.get_surface_measure() * p.get_w(); return integrand; } diff --git a/Bembel/src/LaplaceBeltrami/LaplaceBeltramiOperatorBase.hpp b/Bembel/src/LaplaceBeltrami/LaplaceBeltramiOperatorBase.hpp index ca96a1a51..613a7111b 100644 --- a/Bembel/src/LaplaceBeltrami/LaplaceBeltramiOperatorBase.hpp +++ b/Bembel/src/LaplaceBeltrami/LaplaceBeltramiOperatorBase.hpp @@ -30,15 +30,9 @@ class LaplaceBeltramiOperatorBase : public LocalOperatorBase { const unsigned int elements_per_direction = (1 << super_space.get_refinement_level()); const double h = 1. / ((double)(elements_per_direction)); - - // compute surface measures from tangential derivatives - double x_kappa = p1.segment<3>(6).cross(p1.segment<3>(9)).norm(); - // integrand without basis functions - double integrand = x_kappa * p1(2) / (h * h); - + double integrand = p1.get_surface_measure() * p1.get_w() / (h * h); super_space.addScaledSurfaceGradientInteraction(intval, integrand, p1, p2); - return; } }; diff --git a/Bembel/src/LinearForm/DirichletTrace.hpp b/Bembel/src/LinearForm/DirichletTrace.hpp index fd6e10e4a..fbc7c7186 100644 --- a/Bembel/src/LinearForm/DirichletTrace.hpp +++ b/Bembel/src/LinearForm/DirichletTrace.hpp @@ -38,30 +38,12 @@ class DirichletTrace : public LinearFormBase, Scalar> { void evaluateIntegrand_impl( const T &super_space, const SurfacePoint &p, Eigen::Matrix *intval) const { - auto polynomial_degree = super_space.get_polynomial_degree(); - auto polynomial_degree_plus_one_squared = - (polynomial_degree + 1) * (polynomial_degree + 1); - - // get evaluation points on unit square - auto s = p.segment<2>(0); - - // get quadrature weights - auto ws = p(2); - - // get points on geometry and tangential derivatives - auto x_f = p.segment<3>(3); - auto x_f_dx = p.segment<3>(6); - auto x_f_dy = p.segment<3>(9); - // compute surface measures from tangential derivatives - auto x_kappa = x_f_dx.cross(x_f_dy).norm(); - + auto x_kappa = p.get_f_dx().cross(p.get_f_dy()).norm(); // integrand without basis functions - auto integrand = function_(x_f) * x_kappa * ws; - + auto integrand = function_(p.get_f()) * p.get_surface_measure() * p.get_w(); // multiply basis functions with integrand - super_space.addScaledBasis(intval, integrand, s); - + super_space.addScaledBasis(intval, integrand, p.get_xi()); return; } diff --git a/Bembel/src/LinearForm/RotatedTangentialTrace.hpp b/Bembel/src/LinearForm/RotatedTangentialTrace.hpp index 69861501b..9ae7e082f 100644 --- a/Bembel/src/LinearForm/RotatedTangentialTrace.hpp +++ b/Bembel/src/LinearForm/RotatedTangentialTrace.hpp @@ -40,45 +40,23 @@ class RotatedTangentialTrace void evaluateIntegrand_impl( const T &super_space, const SurfacePoint &p, Eigen::Matrix *intval) const { - auto polynomial_degree = super_space.get_polynomial_degree(); - auto polynomial_degree_plus_one_squared = - (polynomial_degree + 1) * (polynomial_degree + 1); - - // get evaluation points on unit square - auto s = p.segment<2>(0); - - // get quadrature weights - auto ws = p(2); - - // get points on geometry and tangential derivatives - auto x_f = p.segment<3>(3); - auto x_f_dx = p.segment<3>(6); - auto x_f_dy = p.segment<3>(9); - - // compute surface measures from tangential derivatives - auto x_n = x_f_dx.cross(x_f_dy).normalized(); - // tangential component + quadrature weights // use n x f x n = f-n to avoid troubles with -flto flag in combination // of .cross() - auto fun_x_f = function_(x_f); - auto tangential_component = (fun_x_f - fun_x_f.dot(x_n) * x_n) * ws; + Eigen::Matrix fun_x_f = function_(p.get_f()); + Eigen::Vector3d x_n = p.get_unit_normal(); + Eigen::Matrix tangential_component = + (fun_x_f - fun_x_f.dot(x_n) * x_n) * p.get_w(); // extract tangential component - auto component_x = x_f_dx.dot(tangential_component); - auto component_y = x_f_dy.dot(tangential_component); + Eigen::Matrix components = + p.get_jacobian().transpose() * tangential_component; // evaluate shape functions - auto phiPhiVec = super_space.basis(s); - - // multiply basis functions with integrand - Eigen::Matrix phiPhiMat( - polynomial_degree_plus_one_squared, 2); - phiPhiMat.col(0) = component_x * phiPhiVec; - phiPhiMat.col(1) = component_y * phiPhiVec; + auto phiPhiVec = super_space.basis(p.get_xi()); // compute integrals - (*intval) += phiPhiMat; + (*intval) += phiPhiVec * components.transpose(); return; } diff --git a/Bembel/src/LinearOperator/Dummy/DummyOperator.hpp b/Bembel/src/LinearOperator/Dummy/DummyOperator.hpp index 2f8c2fb0d..8419f10a7 100644 --- a/Bembel/src/LinearOperator/Dummy/DummyOperator.hpp +++ b/Bembel/src/LinearOperator/Dummy/DummyOperator.hpp @@ -52,8 +52,8 @@ class DummyOperator : public LinearOperatorBase { const T &super_space, const SurfacePoint &p1, const SurfacePoint &p2, Eigen::Matrix::Scalar, Eigen::Dynamic, Eigen::Dynamic> *intval) const { - (*intval)(0, 0) += - test_func_(p1.segment(3, 2), p2.segment(3, 2)) * p1(2) * p2(2); + (*intval)(0, 0) += test_func_(p1.get_f().head<2>(), p2.get_f().head<2>()) * + p1.get_w() * p2.get_w(); return; } diff --git a/Bembel/src/Maxwell/SingleLayerOperator.hpp b/Bembel/src/Maxwell/SingleLayerOperator.hpp index 00220c765..62cdc15c6 100644 --- a/Bembel/src/Maxwell/SingleLayerOperator.hpp +++ b/Bembel/src/Maxwell/SingleLayerOperator.hpp @@ -43,70 +43,40 @@ class MaxwellSingleLayerOperator Eigen::Matrix< typename LinearOperatorTraits::Scalar, Eigen::Dynamic, Eigen::Dynamic> *intval) const { - auto polynomial_degree = super_space.get_polynomial_degree(); - auto polynomial_degree_plus_one_squared = - (polynomial_degree + 1) * (polynomial_degree + 1); - - // get evaluation points on unit square - auto s = p1.segment<2>(0); - auto t = p2.segment<2>(0); - - // get quadrature weights - auto ws = p1(2); - auto wt = p2(2); - - // get points on geometry and tangential derivatives - auto x_f = p1.segment<3>(3); - auto x_f_dx = p1.segment<3>(6); - auto x_f_dy = p1.segment<3>(9); - auto y_f = p2.segment<3>(3); - auto y_f_dx = p2.segment<3>(6); - auto y_f_dy = p2.segment<3>(9); - // compute surface measures from tangential derivatives auto h = 1. / (1 << super_space.get_refinement_level()); // h = 1 ./ (2^M) // integrand without basis functions - auto kernel_evaluation = evaluateKernel(x_f, y_f) * ws * wt; + auto kernel_evaluation = + evaluateKernel(p1.get_f(), p2.get_f()) * p1.get_w() * p2.get_w(); auto integrand_vector = kernel_evaluation; auto integrand_divergence = -kernel_evaluation / wavenumber2_ / h / h; // vector part: mulitply shape functions with integrand and add to buffer - super_space.addScaledVectorBasisInteraction(intval, integrand_vector, s, t, - x_f_dx, x_f_dy, y_f_dx, y_f_dy); + super_space.addScaledVectorBasisInteraction( + intval, integrand_vector, p1.get_xi(), p2.get_xi(), p1.get_f_dx(), + p1.get_f_dy(), p2.get_f_dx(), p2.get_f_dy()); // divergence part: multiply shape functions with integrand and add to // buffer super_space.addScaledVectorBasisDivergenceInteraction( - intval, integrand_divergence, s, t); + intval, integrand_divergence, p1.get_xi(), p2.get_xi()); return; } Eigen::Matrix, 4, 4> evaluateFMMInterpolation_impl( const SurfacePoint &p1, const SurfacePoint &p2) const { - // get evaluation points on unit square - auto s = p1.segment<2>(0); - auto t = p2.segment<2>(0); - - // get points on geometry and tangential derivatives - auto x_f = p1.segment<3>(3); - auto x_f_dx = p1.segment<3>(6); - auto x_f_dy = p1.segment<3>(9); - auto y_f = p2.segment<3>(3); - auto y_f_dx = p2.segment<3>(6); - auto y_f_dy = p2.segment<3>(9); - // evaluate kernel - auto kernel = evaluateKernel(x_f, y_f); + auto kernel = evaluateKernel(p1.get_f(), p2.get_f()); // interpolation Eigen::Matrix, 4, 4> intval; intval.setZero(); - intval(0, 0) = kernel * x_f_dx.dot(y_f_dx); - intval(0, 2) = kernel * x_f_dx.dot(y_f_dy); - intval(2, 0) = kernel * x_f_dy.dot(y_f_dx); - intval(2, 2) = kernel * x_f_dy.dot(y_f_dy); + intval(0, 0) = kernel * p1.get_f_dx().dot(p2.get_f_dx()); + intval(0, 2) = kernel * p1.get_f_dx().dot(p2.get_f_dy()); + intval(2, 0) = kernel * p1.get_f_dy().dot(p2.get_f_dx()); + intval(2, 2) = kernel * p1.get_f_dy().dot(p2.get_f_dy()); intval(1, 1) = -kernel / wavenumber2_; intval(1, 3) = -kernel / wavenumber2_; intval(3, 1) = -kernel / wavenumber2_; diff --git a/Bembel/src/Maxwell/SingleLayerPotential.hpp b/Bembel/src/Maxwell/SingleLayerPotential.hpp index 3b3b1ce5e..bac6e7f03 100644 --- a/Bembel/src/Maxwell/SingleLayerPotential.hpp +++ b/Bembel/src/Maxwell/SingleLayerPotential.hpp @@ -41,32 +41,20 @@ class MaxwellSingleLayerPotential const ElementTreeNode &element, const Eigen::Vector3d &point, const SurfacePoint &p) const { - // get evaluation points on unit square - auto s = p.segment<2>(0); - - // get quadrature weights - auto ws = p(2); - - // get points on geometry and tangential derivatives - auto x_f = p.segment<3>(3); - // compute surface measures from tangential derivatives auto h = element.get_h(); - // evaluate kernel - auto kernel = evaluateKernel(point, x_f); - auto kernel_gradient = evaluateKernelGrad(point, x_f); - + auto kernel = evaluateKernel(point, p.get_f()); + auto kernel_gradient = evaluateKernelGrad(point, p.get_f()); // assemble Galerkin solution auto scalar_part = fun_ev.evaluate(element, p); auto divergence_part = fun_ev.evaluateDiv(element, p); - // integrand without basis functions, note that the surface measure // disappears for the divergence // auto integrand = kernel * scalar_part * ws; auto integrand = (kernel * scalar_part + 1. / wavenumber2_ * kernel_gradient * divergence_part) * - ws; + p.get_w(); return integrand; } diff --git a/Bembel/src/util/surfaceL2error.hpp b/Bembel/src/util/surfaceL2error.hpp index 2164a5156..33b8ae1dc 100644 --- a/Bembel/src/util/surfaceL2error.hpp +++ b/Bembel/src/util/surfaceL2error.hpp @@ -23,8 +23,7 @@ double surfaceL2error(const AnsatzSpace &ansatz_space, GaussSquare GS; auto Q = GS[deg]; SurfacePoint qp; - const auto longvec = - (ansatz_space.get_transformation_matrix() * vec).eval(); + const auto longvec = (ansatz_space.get_transformation_matrix() * vec).eval(); const auto &super_space = ansatz_space.get_superspace(); const ElementTree &et = super_space.get_mesh().get_element_tree(); const unsigned int number_of_elements = et.get_number_of_elements(); @@ -34,24 +33,14 @@ double surfaceL2error(const AnsatzSpace &ansatz_space, for (auto element = et.cpbegin(); element != et.cpend(); ++element) { for (auto i = 0; i < Q.w_.size(); ++i) { super_space.map2surface(*element, Q.xi_.col(i), Q.w_(i), &qp); - // get evaluation points on unit square - const auto &s = qp.segment<2>(0); - // get quadrature weights - Scalar ws = qp(2); - // get points on geometry and tangential derivatives - const auto &x_f = qp.segment<3>(3); - const auto &x_f_dx = qp.segment<3>(6); - const auto &x_f_dy = qp.segment<3>(9); - const auto &normal = x_f_dx.cross(x_f_dy); - // compute surface measures from tangential derivatives - Scalar x_kappa = normal.norm(); // integrand without basis functions const Scalar val = longvec.segment(n_shape_fun * element->id_, n_shape_fun).transpose() * - super_space.basis(s); - retval += x_kappa * Q.w_(i) * element->get_h() * element->get_h() * - (functor(x_f) - val / element->get_h()) * - (functor(x_f) - val / element->get_h()); + super_space.basis(qp.get_xi()); + retval += qp.get_surface_measure() * Q.w_(i) * element->get_h() * + element->get_h() * + (functor(qp.get_f()) - val / element->get_h()) * + (functor(qp.get_f()) - val / element->get_h()); } } return sqrt(retval); diff --git a/tests/test_SurfacePointUpdate.cpp b/tests/test_SurfacePointUpdate.cpp index c6cacba40..87344588a 100644 --- a/tests/test_SurfacePointUpdate.cpp +++ b/tests/test_SurfacePointUpdate.cpp @@ -16,31 +16,26 @@ int main() { using namespace Bembel; + // initialize geometry Test::TestGeometryWriter::writeScreen(); - Bembel::Geometry geometry("test_Screen.dat"); assert(geometry.get_geometry().size() == 1); + // initialize tolerance + double tol = Test::Constants::test_tolerance_geometry; + for (auto x : Test::Constants::eq_points) { for (auto y : Test::Constants::eq_points) { auto pt = Eigen::Vector2d(x, y); - SurfacePoint srf_pt, srf_pt_ref; - - srf_pt_ref.head(2) = pt; - srf_pt_ref(2) = 3.1415; - srf_pt_ref.segment(3, 3) = geometry.get_geometry()[0].eval(pt); - auto dummy = geometry.get_geometry()[0].evalJacobian(pt); - srf_pt_ref.segment(6, 3) = dummy.col(0); - srf_pt_ref.segment(9, 3) = dummy.col(1); - - auto point = geometry.get_geometry()[0].eval(pt); - auto jacobian = geometry.get_geometry()[0].evalJacobian(pt); - + SurfacePoint srf_pt; geometry.get_geometry()[0].updateSurfacePoint(&srf_pt, pt, 3.1415, pt); - - if ((srf_pt.get_data() - srf_pt_ref.get_data()).norm() > - Test::Constants::test_tolerance_geometry) - return 1; + assert((srf_pt.get_xi() - pt).norm() < tol); + assert(std::abs(srf_pt.get_w() - 3.1415) < tol); + assert((srf_pt.get_f() - geometry.get_geometry()[0].eval(pt)).norm() < + tol); + assert( + (srf_pt.get_jacobian() - geometry.get_geometry()[0].evalJacobian(pt)) + .norm() < tol); } } return 0; From 9522ab0de0d72284b56ef6d5d5371dbb9d3cb782 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20D=C3=B6lz?= Date: Wed, 8 Feb 2023 10:44:43 +0100 Subject: [PATCH 08/20] fix compile error --- Bembel/src/DuffyTrick/evaluateBilinearForm.hpp | 4 ++-- Bembel/src/Geometry/SurfacePoint.hpp | 1 + Bembel/src/H2Matrix/H2Matrix.hpp | 2 +- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/Bembel/src/DuffyTrick/evaluateBilinearForm.hpp b/Bembel/src/DuffyTrick/evaluateBilinearForm.hpp index e1d7f55cb..5b3dd6f3e 100644 --- a/Bembel/src/DuffyTrick/evaluateBilinearForm.hpp +++ b/Bembel/src/DuffyTrick/evaluateBilinearForm.hpp @@ -22,8 +22,8 @@ template void evaluateBilinearForm( const LinearOperatorBase& linOp, const T& super_space, const ElementTreeNode& e1, const ElementTreeNode& e2, - const CubatureVector& GS, const std::vector& ffield_qnodes1, - const std::vector& ffield_qnodes2, + const CubatureVector& GS, const ElementSurfacePoints& ffield_qnodes1, + const ElementSurfacePoints& ffield_qnodes2, Eigen::Matrix::Scalar, Eigen::Dynamic, Eigen::Dynamic>* intval) { ////////////////////////////////////////////////////////////////////////////// diff --git a/Bembel/src/Geometry/SurfacePoint.hpp b/Bembel/src/Geometry/SurfacePoint.hpp index 18492a640..bbd7e0feb 100644 --- a/Bembel/src/Geometry/SurfacePoint.hpp +++ b/Bembel/src/Geometry/SurfacePoint.hpp @@ -50,3 +50,4 @@ typedef std::vector> ElementSurfacePoints; #endif // BEMBEL_SRC_GEOMETRY_SURFACEPOINT_HPP_ + diff --git a/Bembel/src/H2Matrix/H2Matrix.hpp b/Bembel/src/H2Matrix/H2Matrix.hpp index 80288b3bb..590a91759 100644 --- a/Bembel/src/H2Matrix/H2Matrix.hpp +++ b/Bembel/src/H2Matrix/H2Matrix.hpp @@ -134,7 +134,7 @@ class H2Matrix : public EigenBase> { Bembel::GaussSquare GS; auto super_space = ansatz_space.get_superspace(); auto ffield_deg = linOp.get_FarfieldQuadratureDegree(polynomial_degree); - auto ffield_qnodes = + std::vector ffield_qnodes = Bembel::DuffyTrick::computeFfieldQnodes(super_space, GS[ffield_deg]); const int NumberOfFMMComponents = Bembel::LinearOperatorTraits::NumberOfFMMComponents; From 2340a3b833b11184da6ae23a81453e5bf81edad7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20D=C3=B6lz?= Date: Sun, 12 Feb 2023 14:26:15 +0100 Subject: [PATCH 09/20] refactoring of shape functions, this made the code actually slower --- Bembel/Spline | 1 + Bembel/src/AnsatzSpace/Projector.hpp | 12 +- Bembel/src/Geometry/Patch.hpp | 154 ++++++++------------- Bembel/src/Geometry/SurfacePoint.hpp | 28 +--- Bembel/src/Spline/Bernstein.hpp | 181 ++++--------------------- Bembel/src/Spline/Bezierextraction.hpp | 28 +--- Bembel/src/Spline/Localize.hpp | 20 ++- Bembel/src/Spline/ShapeFunctions.hpp | 138 ++++++++++--------- tests/CMakeLists.txt | 1 + tests/test_SurfacePointUpdate.cpp | 6 +- 10 files changed, 188 insertions(+), 381 deletions(-) diff --git a/Bembel/Spline b/Bembel/Spline index 674c68fba..d286149ec 100644 --- a/Bembel/Spline +++ b/Bembel/Spline @@ -29,6 +29,7 @@ #include #include "src/util/Constants.hpp" +#include "src/util/Macros.hpp" #include "src/LinearOperator/DifferentialFormEnum.hpp" #include "src/Spline/DeBoor.hpp" diff --git a/Bembel/src/AnsatzSpace/Projector.hpp b/Bembel/src/AnsatzSpace/Projector.hpp index 63bc7770e..aaf68396c 100644 --- a/Bembel/src/AnsatzSpace/Projector.hpp +++ b/Bembel/src/AnsatzSpace/Projector.hpp @@ -149,18 +149,16 @@ inline _proj_info makeLocalProjectionTriplets( // Here, we suddenly use the degree for basisevaluation, i.e., // maximal_polynomial_degree-1. This is confusing, but correct and tested. { - double vals_y[Constants::MaxP + 1]; - double vals_x[Constants::MaxP + 1]; for (int iy = 0; iy < maximal_polynomial_degree; ++iy) { - Bembel::Basis::ShapeFunctionHandler::evalBasis( - maximal_polynomial_degree - 1, vals_y, mask[iy]); + Eigen::VectorXd vals_y = Bembel::Basis::ShapeFunctionHandler::evalBasis( + maximal_polynomial_degree - 1, mask[iy]); for (int ix = 0; ix < maximal_polynomial_degree; ++ix) { - Bembel::Basis::ShapeFunctionHandler::evalBasis( - maximal_polynomial_degree - 1, vals_x, mask[ix]); + Eigen::VectorXd vals_x = Bembel::Basis::ShapeFunctionHandler::evalBasis( + maximal_polynomial_degree - 1, mask[ix]); for (int jy = 0; jy < maximal_polynomial_degree; ++jy) { for (int jx = 0; jx < maximal_polynomial_degree; ++jx) { system(iy * masksize + ix, jy * maximal_polynomial_degree + jx) = - vals_x[jx] * vals_y[jy]; + vals_x(jx) * vals_y(jy); } } } diff --git a/Bembel/src/Geometry/Patch.hpp b/Bembel/src/Geometry/Patch.hpp index 0f1676f4a..a1535c43d 100644 --- a/Bembel/src/Geometry/Patch.hpp +++ b/Bembel/src/Geometry/Patch.hpp @@ -90,33 +90,23 @@ class Patch { Spl::Rescale(reference_point(1), unique_knots_y_[y_location], unique_knots_y_[y_location + 1]); - double *xbasis = new double[polynomial_degree_x_]; - double *ybasis = new double[polynomial_degree_y_]; - - Bembel::Basis::ShapeFunctionHandler::evalBasis(polynomial_degree_x_ - 1, - xbasis, scaledx); - Bembel::Basis::ShapeFunctionHandler::evalBasis(polynomial_degree_y_ - 1, - ybasis, scaledy); - - double tmp[4] = {0., 0., 0., 0.}; + Eigen::Vector4d tmp = Eigen::Vector4d::Zero(); + Eigen::VectorXd xbasis = Bembel::Basis::ShapeFunctionHandler::evalBasis( + polynomial_degree_x_ - 1, scaledx); + Eigen::VectorXd ybasis = Bembel::Basis::ShapeFunctionHandler::evalBasis( + polynomial_degree_y_ - 1, scaledy); for (int i = 0; i < polynomial_degree_x_; i++) { for (int j = 0; j < polynomial_degree_y_; j++) { - const double tpbasisval = xbasis[i] * ybasis[j]; + const double tpbasisval = xbasis(i) * ybasis(j); const int accs = 4 * (numy * (polynomial_degree_x_ * x_location + i) + polynomial_degree_y_ * y_location + j); -#pragma omp simd - for (int k = 0; k < 4; k++) tmp[k] += data_[accs + k] * tpbasisval; + for (int k = 0; k < 4; k++) + tmp(k) += data_[accs + k] * tpbasisval; } } - delete[] xbasis; - delete[] ybasis; - - Eigen::Vector3d out(tmp[0], tmp[1], tmp[2]); - // Rescaling by the NRBS weight, i.e. projection to 3D from 4D hom - - return out / tmp[3]; + return tmp.head<3>() / tmp[3]; } Eigen::Matrix evalJacobian( @@ -133,61 +123,41 @@ class Patch { Spl::Rescale(reference_point(1), unique_knots_y_[y_location], unique_knots_y_[y_location + 1]); - double *xbasis = new double[polynomial_degree_x_]; - double *ybasis = new double[polynomial_degree_y_]; - double *xbasisD = new double[polynomial_degree_x_]; - double *ybasisD = new double[polynomial_degree_y_]; - - Bembel::Basis::ShapeFunctionHandler::evalBasis(polynomial_degree_x_ - 1, - xbasis, scaledx); - Bembel::Basis::ShapeFunctionHandler::evalBasis(polynomial_degree_y_ - 1, - ybasis, scaledy); - Bembel::Basis::ShapeFunctionHandler::evalDerBasis(polynomial_degree_x_ - 1, - xbasisD, scaledx); - Bembel::Basis::ShapeFunctionHandler::evalDerBasis(polynomial_degree_y_ - 1, - ybasisD, scaledy); - - double tmp[4] = {0., 0., 0., 0.}; - double tmpDx[4] = {0., 0., 0., 0.}; - double tmpDy[4] = {0., 0., 0., 0.}; + Eigen::Vector4d tmp = Eigen::Vector4d::Zero(); + Eigen::Vector4d tmpDx = Eigen::Vector4d::Zero(); + Eigen::Vector4d tmpDy = Eigen::Vector4d::Zero(); + Eigen::VectorXd xbasis = Bembel::Basis::ShapeFunctionHandler::evalBasis( + polynomial_degree_x_ - 1, scaledx); + Eigen::VectorXd ybasis = Bembel::Basis::ShapeFunctionHandler::evalBasis( + polynomial_degree_y_ - 1, scaledy); + Eigen::VectorXd xbasisD = Bembel::Basis::ShapeFunctionHandler::evalDerBasis( + polynomial_degree_x_ - 1, scaledx); + Eigen::VectorXd ybasisD = Bembel::Basis::ShapeFunctionHandler::evalDerBasis( + polynomial_degree_y_ - 1, scaledy); for (int i = 0; i < polynomial_degree_x_; i++) { for (int j = 0; j < polynomial_degree_y_; j++) { - const double tpbasisval = xbasis[i] * ybasis[j]; - const double tpbasisvalDx = xbasisD[i] * ybasis[j]; - const double tpbasisvalDy = xbasis[i] * ybasisD[j]; + const double tpbasisval = xbasis(i) * ybasis(j); + const double tpbasisvalDx = xbasisD(i) * ybasis(j); + const double tpbasisvalDy = xbasis(i) * ybasisD(j); const int accs = 4 * (numy * (polynomial_degree_x_ * x_location + i) + polynomial_degree_y_ * y_location + j); // Here I add up the values of the basis functions in the dc // basis -#pragma omp simd for (int k = 0; k < 4; k++) { - tmp[k] += data_[accs + k] * tpbasisval; - tmpDx[k] += data_[accs + k] * tpbasisvalDx; - tmpDy[k] += data_[accs + k] * tpbasisvalDy; + tmp(k) += data_[accs + k] * tpbasisval; + tmpDx(k) += data_[accs + k] * tpbasisvalDx; + tmpDy(k) += data_[accs + k] * tpbasisvalDy; } } } - delete[] xbasis; - delete[] ybasis; - delete[] xbasisD; - delete[] ybasisD; - - Eigen::Matrix out; - - // Eigen::Vector3d out; - - double bot = 1. / (tmp[3] * tmp[3]); - -#pragma omp simd - for (int k = 0; k < 3; k++) { - out(k, 0) = (tmpDx[k] * tmp[3] - tmp[k] * tmpDx[3]) * bot; - out(k, 1) = (tmpDy[k] * tmp[3] - tmp[k] * tmpDy[3]) * bot; - } - - return out; + double botsqr = 1. / (tmp(3) * tmp(3)); + return botsqr * (Eigen::Matrix() + << tmpDx.head<3>() * tmp(3) - tmp.head<3>() * tmpDx(3), + tmpDy.head<3>() * tmp(3) - tmp.head<3>() * tmpDy(3)) + .finished(); } inline Eigen::Matrix evalNormal( @@ -221,43 +191,33 @@ class Patch { const double scaledy = Spl::Rescale(ref_pt(1), unique_knots_y_[y_location], unique_knots_y_[y_location + 1]); - // TODO(Felix) Do not use variable-length arrays in accordance to Google - // Style guide - double *buffer = - new double[2 * (polynomial_degree_x_ + polynomial_degree_y_) + 12]; - for (int i = 0; i < 12; ++i) buffer[i] = 0; - - double *tmp = buffer; - double *tmpDx = tmp + 4; - double *tmpDy = tmpDx + 4; - double *xbasis = tmpDy + 4; - double *ybasis = xbasis + polynomial_degree_x_; - double *xbasisD = ybasis + polynomial_degree_y_; - double *ybasisD = xbasisD + polynomial_degree_x_; - - Bembel::Basis::ShapeFunctionHandler::evalBasis(polynomial_degree_x_ - 1, - xbasis, scaledx); - Bembel::Basis::ShapeFunctionHandler::evalBasis(polynomial_degree_y_ - 1, - ybasis, scaledy); - Bembel::Basis::ShapeFunctionHandler::evalDerBasis(polynomial_degree_x_ - 1, - xbasisD, scaledx); - Bembel::Basis::ShapeFunctionHandler::evalDerBasis(polynomial_degree_y_ - 1, - ybasisD, scaledy); + Eigen::Vector4d tmp = Eigen::Vector4d::Zero(); + Eigen::Vector4d tmpDx = Eigen::Vector4d::Zero(); + Eigen::Vector4d tmpDy = Eigen::Vector4d::Zero(); + + Eigen::VectorXd xbasis = Bembel::Basis::ShapeFunctionHandler::evalBasis( + polynomial_degree_x_ - 1, scaledx); + Eigen::VectorXd ybasis = Bembel::Basis::ShapeFunctionHandler::evalBasis( + polynomial_degree_y_ - 1, scaledy); + Eigen::VectorXd xbasisD = Bembel::Basis::ShapeFunctionHandler::evalDerBasis( + polynomial_degree_x_ - 1, scaledx); + Eigen::VectorXd ybasisD = Bembel::Basis::ShapeFunctionHandler::evalDerBasis( + polynomial_degree_y_ - 1, scaledy); for (int i = 0; i < polynomial_degree_x_; ++i) { for (int j = 0; j < polynomial_degree_y_; ++j) { - const double tpbasisval = xbasis[i] * ybasis[j]; - const double tpbasisvalDx = xbasisD[i] * ybasis[j]; - const double tpbasisvalDy = xbasis[i] * ybasisD[j]; + const double tpbasisval = xbasis(i) * ybasis(j); + const double tpbasisvalDx = xbasisD(i) * ybasis(j); + const double tpbasisvalDy = xbasis(i) * ybasisD(j); const int accs = 4 * (numy * (polynomial_degree_x_ * x_location + i) + polynomial_degree_y_ * y_location + j); // Here I add up the values of the basis functions in the dc // basis for (int k = 0; k < 4; ++k) { - tmp[k] += data_[accs + k] * tpbasisval; - tmpDx[k] += data_[accs + k] * tpbasisvalDx; - tmpDy[k] += data_[accs + k] * tpbasisvalDy; + tmp(k) += data_[accs + k] * tpbasisval; + tmpDx(k) += data_[accs + k] * tpbasisvalDx; + tmpDy(k) += data_[accs + k] * tpbasisvalDy; } } } @@ -265,18 +225,14 @@ class Patch { const double bot = 1. / tmp[3]; const double botsqr = bot * bot; - Eigen::Vector3d tmpvec(tmp[0], tmp[1], tmp[2]); - Eigen::Vector3d tmpDxvec(tmpDx[0], tmpDx[1], tmpDx[2]); - Eigen::Vector3d tmpDyvec(tmpDy[0], tmpDy[1], tmpDy[2]); - srf_pt->set_xi(xi); srf_pt->set_w(w); - srf_pt->set_f(bot * tmpvec); - srf_pt->set_jacobian(botsqr * (Eigen::Matrix() - << tmpDxvec * tmp[3] - tmpvec * tmpDx[3], - tmpDyvec * tmp[3] - tmpvec * tmpDy[3]) - .finished()); - delete[] buffer; + srf_pt->set_f(bot * tmp.head<3>()); + srf_pt->set_jacobian( + botsqr * (Eigen::Matrix() + << tmpDx.head<3>() * tmp(3) - tmp.head<3>() * tmpDx(3), + tmpDy.head<3>() * tmp(3) - tmp.head<3>() * tmpDy(3)) + .finished()); return; } diff --git a/Bembel/src/Geometry/SurfacePoint.hpp b/Bembel/src/Geometry/SurfacePoint.hpp index 7ddafe0e4..6748d8942 100644 --- a/Bembel/src/Geometry/SurfacePoint.hpp +++ b/Bembel/src/Geometry/SurfacePoint.hpp @@ -14,30 +14,9 @@ // #include /** * \ingroup Geometry - * \brief Wrapper class for legacy SurfacePoint represented as an Eigen::Vector - * - * This typedef is essential for any evaluation of a bilinear form. It provides - * all required geometry information, stored as follows: - * (0) x-coordinate of the evaluation point in the parameter domain [0,1]^2 - * of the current element, i.e we map [0,1]^2->element->surface - * (1) y-coordinate of the evaluation point in the parameter domain [0,1]^2 - * of the current element, i.e we map [0,1]^2->element->surface - * (2) a quadrature weight. Can be left empty if not used as part of a - * quadrature. - * (3) x-coordinate of patch eval in space - * (4) y-coordinate of patch eval in space - * (5) z-coordinate of patch eval in space - * (6) x-component of derivative in x-dir - * (7) y-component of derivative in x-dir - * (8) z-component of derivative in x-dir - * (9) x-component of derivative in y-dir - * (10) y-component of derivative in y-dir - * (11) z-component of derivative in y-dir - * For application of the pull-back to the reference domain, one requires the - * jacobian of any point on the surface. Calling eval and evalJacobian of the - * Patch class introduces work that needs to be done twice. The - * updateSurdacePoint method is specialized and should be used, since it avoids - * redundant work. + * \brief This class stores quadrature and geometry information on quadrature + *points and provides various methods to compute differential geometric + *entities. **/ class SurfacePoint { public: @@ -96,4 +75,3 @@ typedef std::vector> ElementSurfacePoints; #endif // BEMBEL_SRC_GEOMETRY_SURFACEPOINT_HPP_ - diff --git a/Bembel/src/Spline/Bernstein.hpp b/Bembel/src/Spline/Bernstein.hpp index 800ae4ef3..bba57155e 100644 --- a/Bembel/src/Spline/Bernstein.hpp +++ b/Bembel/src/Spline/Bernstein.hpp @@ -49,68 +49,34 @@ inline constexpr double Bernstein(double evaluation_point) noexcept { //////////////////////////////////////////////////////////////////////////////// /// Hidden Classes //////////////////////////////////////////////////////////////////////////////// -template +template class HiddenBernsteinClass { public: - static inline T EvalCoefs(T *in, double evaluation_point) noexcept { - return in[N] * Bernstein(evaluation_point) + - HiddenBernsteinClass::EvalCoefs(in, evaluation_point); - } - static inline T EvalDerCoefs(T *in, double evaluation_point) noexcept { - return ( - (in[N + 1] - in[N]) * Bernstein(evaluation_point) + - HiddenBernsteinClass::EvalDerCoefs(in, evaluation_point)); - } - static inline void EvalBasisPEQ(T *in, double evaluation_point) noexcept { - in[N] += Bernstein(evaluation_point); - HiddenBernsteinClass::EvalBasisPEQ(in, evaluation_point); - return; - } - static inline void EvalDerBasisPEQ(T *in, double evaluation_point) noexcept { - in[N] += (P + 1) * (Bernstein(evaluation_point) - - Bernstein(evaluation_point)); - HiddenBernsteinClass::EvalDerBasisPEQ(in, evaluation_point); - return; - } - static inline void EvalBasis(T *in, double evaluation_point) noexcept { + static inline void EvalBasis(double *in, + const double evaluation_point) noexcept { in[N] = Bernstein(evaluation_point); - HiddenBernsteinClass::EvalBasis(in, evaluation_point); + HiddenBernsteinClass::EvalBasis(in, evaluation_point); return; } - static inline void EvalDerBasis(T *in, double evaluation_point) noexcept { + static inline void EvalDerBasis(double *in, + const double evaluation_point) noexcept { in[N] = (P + 1) * (Bernstein(evaluation_point) - Bernstein(evaluation_point)); - HiddenBernsteinClass::EvalDerBasis(in, evaluation_point); + HiddenBernsteinClass::EvalDerBasis(in, evaluation_point); return; } }; -template -class HiddenBernsteinClass { +template +class HiddenBernsteinClass<0, P> { public: - static inline T EvalCoefs(T *in, double evaluation_point) noexcept { - return in[0] * Bernstein<0, P>(evaluation_point); - } - static inline T EvalDerCoefs(T *in, double evaluation_point) noexcept { - // P needs to be passed lower to avoid infinite recursion - return (in[1] - in[0]) * Bernstein<0, P>(evaluation_point); - } - static inline void EvalBasisPEQ(T *in, double evaluation_point) noexcept { - in[0] += Bernstein<0, P>(evaluation_point); - return; - } - static inline void EvalDerBasisPEQ(T *in, double evaluation_point) noexcept { - // P needs to be passed lower to avoid infinite recursion - in[0] += (-P - 1) * Bernstein<0, P>(evaluation_point); - return; - } - - static inline void EvalBasis(T *in, double evaluation_point) noexcept { + static inline void EvalBasis(double *in, + const double evaluation_point) noexcept { in[0] = Bernstein<0, P>(evaluation_point); return; } - static inline void EvalDerBasis(T *in, double evaluation_point) noexcept { - // P needs to be passed lower to avoid infinite recursion + static inline void EvalDerBasis(double *in, + const double evaluation_point) noexcept { in[0] = (-P - 1) * Bernstein<0, P>(evaluation_point); return; } @@ -118,128 +84,39 @@ class HiddenBernsteinClass { // This specialization is needed to get a specialized recursion anchor for the // case P = 0. -template -class HiddenBernsteinClass { +template +class HiddenBernsteinClass<-1, P> { public: - static inline T EvalCoefs(T *in, double evaluation_point) noexcept { - (void)in; - (void)evaluation_point; - assert( - false && - "Pos.A This should not happen. Something is wrong with the recursion"); - }; - static inline T EvalDerCoefs(T *in, double evaluation_point) noexcept { - // P needs to be passed lower to avoid infinite recursion - (void)in; - (void)evaluation_point; - return 0; - }; - static inline void EvalBasis(T *in, double evaluation_point) noexcept { - (void)in; - (void)evaluation_point; + static inline void EvalBasis(double *in, + const double evaluation_point) noexcept { + BEMBEL_UNUSED_(in); + BEMBEL_UNUSED_(evaluation_point); assert( false && "Pos.C This should not happen. Something is wrong with the recursion"); }; - static inline void EvalDerBasis(T *in, double evaluation_point) noexcept { - (void)in; - (void)evaluation_point; - // P needs to be passed lower to avoid infinite recursion - return; - }; - static inline void EvalBasisPEQ(T *in, double evaluation_point) noexcept { - (void)in; - (void)evaluation_point; + static inline void EvalDerBasis(double *in, + const double evaluation_point) noexcept { + BEMBEL_UNUSED_(in); + BEMBEL_UNUSED_(evaluation_point); assert( false && "Pos.C This should not happen. Something is wrong with the recursion"); }; - static inline void EvalDerBasisPEQ(T *in, double evaluation_point) noexcept { - (void)in; - (void)evaluation_point; - // P needs to be passed lower to avoid infinite recursion - return; - }; }; //////////////////////////////////////////////////////////////////////////////// /// Evaluation Routines //////////////////////////////////////////////////////////////////////////////// -template -T EvalBernstein(T *in, double evaluation_point) noexcept { - return HiddenBernsteinClass::EvalCoefs(in, evaluation_point); -} - -template -void EvalBernstein(T *in, const std::vector &evaluation_points, - T *out) noexcept { - const int N = evaluation_points.size(); - for (int i = 0; i < N; i++) - out[i] = HiddenBernsteinClass::EvalCoefs(in, evaluation_points[i]); - return; -} - -template -std::vector EvalBernstein( - T *in, const std::vector &evaluation_points) noexcept { - const int N = evaluation_points.size(); - std::vector out(N); - for (int i = 0; i < N; i++) - out[i] = HiddenBernsteinClass::EvalCoefs(in, evaluation_points[i]); - return out; -} -template -void EvalBernsteinBasisPEQ(T *in, double evaluation_point) noexcept { - HiddenBernsteinClass::EvalBasisPEQ(in, evaluation_point); - return; -} - -template -void EvalBernsteinBasis(T *in, double evaluation_point) noexcept { - HiddenBernsteinClass::EvalBasis(in, evaluation_point); - return; -} -//////////////////////////////////////////////////////////////////////////////// -/// Evaluation of the Derivatives -//////////////////////////////////////////////////////////////////////////////// -template -T EvalBernsteinDer(T *in, double evaluation_point) noexcept { - return P * HiddenBernsteinClass::EvalDerCoefs( - in, evaluation_point); -} - -template -void EvalBernsteinDer(T *in, const std::vector &evaluation_points, - T *out) noexcept { - const int N = evaluation_points.size(); - for (int i = 0; i < N; i++) - out[i] = P * HiddenBernsteinClass::EvalDerCoefs( - in, evaluation_points[i]); - return; -} - -template -std::vector EvalBernsteinDer( - T *in, const std::vector &evaluation_points) noexcept { - const int N = evaluation_points.size(); - std::vector out(N); - for (int i = 0; i < N; i++) - out[i] = P * HiddenBernsteinClass::EvalDerCoefs( - in, evaluation_points[i]); - return out; -} - -template -void EvalBernsteinDerBasisPEQ(T *in, double evaluation_point) noexcept { - in[P] += P * Bernstein

(evaluation_point); - HiddenBernsteinClass::EvalDerBasisPEQ(in, evaluation_point); - return; +template +void EvalBernsteinBasis(double *in, const double evaluation_point) noexcept { + HiddenBernsteinClass::EvalBasis(in, evaluation_point); } -template -void EvalBernsteinDerBasis(T *in, double evaluation_point) noexcept { +template +void EvalBernsteinDerBasis(double *in, const double evaluation_point) noexcept { in[P] = P * Bernstein

(evaluation_point); - HiddenBernsteinClass::EvalDerBasis(in, evaluation_point); + HiddenBernsteinClass

::EvalDerBasis(in, evaluation_point); return; } diff --git a/Bembel/src/Spline/Bezierextraction.hpp b/Bembel/src/Spline/Bezierextraction.hpp index 9fb5086fe..e039ae5b1 100644 --- a/Bembel/src/Spline/Bezierextraction.hpp +++ b/Bembel/src/Spline/Bezierextraction.hpp @@ -53,24 +53,15 @@ Eigen::SparseMatrix MakeProjection(const std::vector &x_knots, Eigen::Matrix tempsolver = Eigen::Matrix::Zero(size_phi, size_phi); - double *coefficients_x = new double[polynomial_degree_x]; - double *coefficients_y = new double[polynomial_degree_y]; - - for (int i = 0; i < polynomial_degree_x; i++) coefficients_x[i] = 0; - - for (int i = 0; i < polynomial_degree_y; i++) coefficients_y[i] = 0; - - auto BezBasis = [](std::vector uniq, int deg, int pos, double pt, - double *coef) { + auto BezBasis = [](std::vector uniq, int deg, int pos, double pt) { std::div_t loc = std::div(pos, deg); if (pt > uniq[loc.quot + 1] || pt < uniq[loc.quot]) { return 0.; } else { - coef[loc.rem] = 1; - double out = Bembel::Basis::ShapeFunctionHandler::evalCoef( - deg - 1, coef, Rescale(pt, uniq[loc.quot], uniq[loc.quot + 1])); - coef[loc.rem] = 0; - return out; + Eigen::VectorXd unit = Eigen::VectorXd::Zero(deg); + unit(loc.rem) = 1.; + return unit.dot(Bembel::Basis::ShapeFunctionHandler::evalBasis( + deg - 1, Rescale(pt, uniq[loc.quot], uniq[loc.quot + 1]))); } }; @@ -81,18 +72,13 @@ Eigen::SparseMatrix MakeProjection(const std::vector &x_knots, for (int y = 0; y < size_phi_y; y++) { for (int x = 0; x < size_phi_x; x++) { tempsolver(y * size_phi_x + x, ix * size_phi_y + iy) = - BezBasis(x_unique_knots, polynomial_degree_x, ix, xpoint[x], - coefficients_x) * - BezBasis(y_unique_knots, polynomial_degree_y, iy, ypoint[y], - coefficients_y); + BezBasis(x_unique_knots, polynomial_degree_x, ix, xpoint[x]) * + BezBasis(y_unique_knots, polynomial_degree_y, iy, ypoint[y]); } } } } - delete[] coefficients_x; - delete[] coefficients_y; - Eigen::FullPivLU> lu_decomp(tempsolver); assert(lu_decomp.rank() == size_phi); /// The inverse matrix, for the solition of the linear system to come. diff --git a/Bembel/src/Spline/Localize.hpp b/Bembel/src/Spline/Localize.hpp index efa0cf49a..275095613 100644 --- a/Bembel/src/Spline/Localize.hpp +++ b/Bembel/src/Spline/Localize.hpp @@ -64,18 +64,14 @@ inline std::vector MakeInterpolationPoints( * \brief returns the coefficients to represent a function in the Bernstein * basis on [0,1]. **/ -inline Eigen::Matrix GetInterpolationMatrix( - int polynomial_degree, const std::vector &mask) { - Eigen::Matrix interpolationMatrix(polynomial_degree + 1, - polynomial_degree + 1); - - double val[Constants::MaxP + 1]; - for (int j = 0; j < polynomial_degree + 1; j++) { - Bembel::Basis::ShapeFunctionHandler::evalBasis(polynomial_degree, val, - mask[j]); - for (int i = 0; i < polynomial_degree + 1; i++) - interpolationMatrix(j, i) = val[i]; - } +inline Eigen::MatrixXd GetInterpolationMatrix(int polynomial_degree, + const std::vector &mask) { + Eigen::MatrixXd interpolationMatrix(polynomial_degree + 1, + polynomial_degree + 1); + + for (int j = 0; j < polynomial_degree + 1; j++) + interpolationMatrix.row(j) = Bembel::Basis::ShapeFunctionHandler::evalBasis( + polynomial_degree, mask[j]); return interpolationMatrix.inverse(); } diff --git a/Bembel/src/Spline/ShapeFunctions.hpp b/Bembel/src/Spline/ShapeFunctions.hpp index abee62cfc..989dff090 100644 --- a/Bembel/src/Spline/ShapeFunctions.hpp +++ b/Bembel/src/Spline/ShapeFunctions.hpp @@ -15,9 +15,6 @@ namespace Bembel { namespace Basis { -using funptr_doubleOut_doubleptrDoubleIn = double (*)(double*, double); -using funptr_voidOut_doubleptrDoubleIn = void (*)(double*, double); - /** * \ingroup Spline * \brief These routines implement a template recursion that allows to choose a @@ -27,79 +24,92 @@ using funptr_voidOut_doubleptrDoubleIn = void (*)(double*, double); template class PSpecificShapeFunctionHandler { public: - inline static double evalCoef(int p, double* ar, double x) { - return p == P ? Bembel::Basis::EvalBernstein(ar, x) - : PSpecificShapeFunctionHandler

::evalCoef(p, ar, x); - } - inline static double evalDerCoef(int p, double* ar, double x) { - return p == P ? Bembel::Basis::EvalBernsteinDer(ar, x) - : PSpecificShapeFunctionHandler

::evalDerCoef(p, ar, x); - } - inline static void evalBasis(int p, double* ar, double x) { - return p == P ? Bembel::Basis::EvalBernsteinBasis(ar, x) + /** + * \brief Evaluates the Bernstein basis of degree p at point x into ar. + */ + inline static void evalBasis(const int p, double* ar, const double x) { + return p == P ? Bembel::Basis::EvalBernsteinBasis

(ar, x) : PSpecificShapeFunctionHandler

::evalBasis(p, ar, x); } - inline static void evalDerBasis(int p, double* ar, double x) { + /** + * \brief Evaluates the Bernstein basis of degree p at point x into an + * Eigen::VectorXd. + * + * \attention For performance critical applications use the pointer version of + * this function to avoid slow memory allocation. + */ + inline static Eigen::Matrix evalBasis( + const int p, const double x) { + Eigen::VectorXd eval(p + 1); + evalBasis(p, eval.data(), x); + return eval; + } + /** + * \brief Evaluates the derivative of the Bernstein basis of degree p at point + * x into ar. + */ + inline static void evalDerBasis(const int p, double* ar, const double x) { return p == P - ? Bembel::Basis::EvalBernsteinDerBasis(ar, x) + ? Bembel::Basis::EvalBernsteinDerBasis

(ar, x) : PSpecificShapeFunctionHandler

::evalDerBasis(p, ar, x); } - inline static constexpr funptr_doubleOut_doubleptrDoubleIn ptrEvalCoef( - int p) { - return p == P ? &Bembel::Basis::EvalBernstein - : PSpecificShapeFunctionHandler

::ptrEvalCoef(p); - } - inline static constexpr funptr_doubleOut_doubleptrDoubleIn ptrEvalDerCoef( - int p) { - return p == P ? &Bembel::Basis::EvalBernsteinDer - : PSpecificShapeFunctionHandler

::ptrEvalDerCoef(p); - } - inline static constexpr funptr_voidOut_doubleptrDoubleIn ptrEvalBasis(int p) { - return p == P ? &Bembel::Basis::EvalBernsteinBasis - : PSpecificShapeFunctionHandler

::ptrEvalBasis(p); - } - inline static constexpr funptr_voidOut_doubleptrDoubleIn ptrEvalDerBasis( - int p) { - return p == P ? &Bembel::Basis::EvalBernsteinDerBasis - : PSpecificShapeFunctionHandler

::ptrEvalDerBasis(p); - } - inline static constexpr bool checkP(int p) { - static_assert(P > 0, "Polynomial degree must be larger than zero"); - return p <= Constants::MaxP; + /** + * \brief Evaluates the derivative of the Bernstein basis of degree p into an + * Eigen::VectorXd. + * + * \attention For performance critical applications use the pointer version of + * this function to avoid slow memory allocation. + */ + inline static Eigen::Matrix evalDerBasis( + const int p, const double x) { + Eigen::VectorXd eval(p + 1); + evalDerBasis(p, eval.data(), x); + return eval; } }; template <> class PSpecificShapeFunctionHandler<0> { public: - inline static double evalCoef(int p, double* ar, double x) { - return Bembel::Basis::EvalBernstein(ar, x); - } - inline static double evalDerCoef(int p, double* ar, double x) { - return Bembel::Basis::EvalBernsteinDer(ar, x); - } - inline static void evalBasis(int p, double* ar, double x) { - return Bembel::Basis::EvalBernsteinBasis(ar, x); - } - inline static void evalDerBasis(int p, double* ar, double x) { - return Bembel::Basis::EvalBernsteinDerBasis(ar, x); - } - inline static constexpr funptr_doubleOut_doubleptrDoubleIn ptrEvalCoef( - int p) { - return &Bembel::Basis::EvalBernstein; - } - inline static constexpr funptr_doubleOut_doubleptrDoubleIn ptrEvalDerCoef( - int p) { - return &Bembel::Basis::EvalBernsteinDer; - } - inline static constexpr funptr_voidOut_doubleptrDoubleIn ptrEvalBasis(int p) { - return &Bembel::Basis::EvalBernsteinBasis; - } - inline static constexpr funptr_voidOut_doubleptrDoubleIn ptrEvalDerBasis( - int p) { - return &Bembel::Basis::EvalBernsteinDerBasis; + /** + * \brief Evaluates the Bernstein basis of degree p at point x into ar. + */ + inline static void evalBasis(const int p, double* ar, const double x) { + Bembel::Basis::EvalBernsteinBasis<0>(ar, x); + return; + } + /** + * \brief Evaluates the Bernstein basis of degree p at point x into an + * Eigen::VectorXd. + * + * \attention For performance critical applications use the pointer version of + * this function to avoid slow memory allocation. + */ + inline static Eigen::VectorXd evalBasis(const int p, const double x) { + Eigen::VectorXd eval(p + 1); + evalBasis(p, eval.data(), x); + return eval; + } + /** + * \brief Evaluates the derivative of the Bernstein basis of degree p at point + * x into ar. + */ + inline static void evalDerBasis(const int p, double* ar, const double x) { + Bembel::Basis::EvalBernsteinDerBasis<0>(ar, x); + return; + } + /** + * \brief Evaluates the derivative of the Bernstein basis of degree p into an + * Eigen::VectorXd. + * + * \attention For performance critical applications use the pointer version of + * this function to avoid slow memory allocation. + */ + inline static Eigen::VectorXd evalDerBasis(const int p, const double x) { + Eigen::VectorXd eval(p + 1); + evalDerBasis(p, eval.data(), x); + return eval; } - inline static constexpr bool checkP(int p) { return Constants::MaxP >= 0; } }; using ShapeFunctionHandler = PSpecificShapeFunctionHandler; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 44c5fd986..4c6d838a3 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -1,5 +1,6 @@ set(UNITTESTS test_GeometryImportAndEval + test_Spline test_SurfacePointUpdate test_Projector test_Glue diff --git a/tests/test_SurfacePointUpdate.cpp b/tests/test_SurfacePointUpdate.cpp index 87344588a..4ab20ca09 100644 --- a/tests/test_SurfacePointUpdate.cpp +++ b/tests/test_SurfacePointUpdate.cpp @@ -26,13 +26,17 @@ int main() { for (auto x : Test::Constants::eq_points) { for (auto y : Test::Constants::eq_points) { - auto pt = Eigen::Vector2d(x, y); + Eigen::Vector2d pt = Eigen::Vector2d(x, y); SurfacePoint srf_pt; geometry.get_geometry()[0].updateSurfacePoint(&srf_pt, pt, 3.1415, pt); assert((srf_pt.get_xi() - pt).norm() < tol); assert(std::abs(srf_pt.get_w() - 3.1415) < tol); + assert((srf_pt.get_f().head<2>() - pt).norm() < tol); + assert(std::abs(srf_pt.get_f()(2) < tol)); assert((srf_pt.get_f() - geometry.get_geometry()[0].eval(pt)).norm() < tol); + assert((srf_pt.get_f_dx() - Eigen::Vector3d(1., 0., 0.)).norm() < tol); + assert((srf_pt.get_f_dy() - Eigen::Vector3d(0, 1., 0.)).norm() < tol); assert( (srf_pt.get_jacobian() - geometry.get_geometry()[0].evalJacobian(pt)) .norm() < tol); From a02bf44ca6538e70b4085d089dce7739ed42e678 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20D=C3=B6lz?= Date: Mon, 13 Feb 2023 11:29:43 +0100 Subject: [PATCH 10/20] refactoring surfacePointUpdate --- .vscode/settings.json | 6 +++ Bembel/src/Geometry/Patch.hpp | 58 ++++++++++++++++++---------- Bembel/src/Geometry/SurfacePoint.hpp | 12 ++++++ tests/test_SurfacePointUpdate.cpp | 5 --- 4 files changed, 56 insertions(+), 25 deletions(-) create mode 100644 .vscode/settings.json diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 000000000..bde692edd --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,6 @@ +{ + "files.associations": { + "quadrature": "cpp", + "spline": "cpp" + } +} \ No newline at end of file diff --git a/Bembel/src/Geometry/Patch.hpp b/Bembel/src/Geometry/Patch.hpp index a1535c43d..da8dbbe1f 100644 --- a/Bembel/src/Geometry/Patch.hpp +++ b/Bembel/src/Geometry/Patch.hpp @@ -101,8 +101,7 @@ class Patch { const double tpbasisval = xbasis(i) * ybasis(j); const int accs = 4 * (numy * (polynomial_degree_x_ * x_location + i) + polynomial_degree_y_ * y_location + j); - for (int k = 0; k < 4; k++) - tmp(k) += data_[accs + k] * tpbasisval; + for (int k = 0; k < 4; k++) tmp(k) += data_[accs + k] * tpbasisval; } } @@ -191,17 +190,33 @@ class Patch { const double scaledy = Spl::Rescale(ref_pt(1), unique_knots_y_[y_location], unique_knots_y_[y_location + 1]); - Eigen::Vector4d tmp = Eigen::Vector4d::Zero(); - Eigen::Vector4d tmpDx = Eigen::Vector4d::Zero(); - Eigen::Vector4d tmpDy = Eigen::Vector4d::Zero(); - - Eigen::VectorXd xbasis = Bembel::Basis::ShapeFunctionHandler::evalBasis( + // resize operations do nothing if the size does not change + srf_pt->get_buffer().resize(6); + srf_pt->get_buffer()[0].resize(4, 1); + srf_pt->get_buffer()[1].resize(4, 2); + srf_pt->get_buffer()[2].resize(polynomial_degree_x_, 1); + srf_pt->get_buffer()[3].resize(polynomial_degree_y_, 1); + srf_pt->get_buffer()[4].resize(polynomial_degree_x_, 1); + srf_pt->get_buffer()[5].resize(polynomial_degree_y_, 1); + + // improve readability + Eigen::MatrixXd &tmp = srf_pt->get_buffer()[0]; + Eigen::MatrixXd &tmpD = srf_pt->get_buffer()[1]; + Eigen::MatrixXd &xbasis = srf_pt->get_buffer()[2]; + Eigen::MatrixXd &ybasis = srf_pt->get_buffer()[3]; + Eigen::MatrixXd &xbasisD = srf_pt->get_buffer()[4]; + Eigen::MatrixXd &ybasisD = srf_pt->get_buffer()[5]; + + tmp.setZero(); + tmpD.setZero(); + + xbasis.col(0) = Bembel::Basis::ShapeFunctionHandler::evalBasis( polynomial_degree_x_ - 1, scaledx); - Eigen::VectorXd ybasis = Bembel::Basis::ShapeFunctionHandler::evalBasis( + ybasis.col(0) = Bembel::Basis::ShapeFunctionHandler::evalBasis( polynomial_degree_y_ - 1, scaledy); - Eigen::VectorXd xbasisD = Bembel::Basis::ShapeFunctionHandler::evalDerBasis( + xbasisD.col(0) = Bembel::Basis::ShapeFunctionHandler::evalDerBasis( polynomial_degree_x_ - 1, scaledx); - Eigen::VectorXd ybasisD = Bembel::Basis::ShapeFunctionHandler::evalDerBasis( + ybasisD.col(0) = Bembel::Basis::ShapeFunctionHandler::evalDerBasis( polynomial_degree_y_ - 1, scaledy); for (int i = 0; i < polynomial_degree_x_; ++i) { @@ -215,24 +230,27 @@ class Patch { // Here I add up the values of the basis functions in the dc // basis for (int k = 0; k < 4; ++k) { - tmp(k) += data_[accs + k] * tpbasisval; - tmpDx(k) += data_[accs + k] * tpbasisvalDx; - tmpDy(k) += data_[accs + k] * tpbasisvalDy; + tmp(k, 0) += data_[accs + k] * tpbasisval; + tmpD(k, 0) += data_[accs + k] * tpbasisvalDx; + tmpD(k, 1) += data_[accs + k] * tpbasisvalDy; } } } - const double bot = 1. / tmp[3]; + const double bot = 1. / tmp(3, 0); const double botsqr = bot * bot; srf_pt->set_xi(xi); srf_pt->set_w(w); - srf_pt->set_f(bot * tmp.head<3>()); - srf_pt->set_jacobian( - botsqr * (Eigen::Matrix() - << tmpDx.head<3>() * tmp(3) - tmp.head<3>() * tmpDx(3), - tmpDy.head<3>() * tmp(3) - tmp.head<3>() * tmpDy(3)) - .finished()); + srf_pt->set_f(bot * tmp.block<3, 1>(0, 0)); + // srf_pt->set_jacobian( + // botsqr * (Eigen::Matrix() + // << tmpDx.head<3>() * tmp(3) - tmp.head<3>() * tmpDx(3), + // tmpDy.head<3>() * tmp(3) - tmp.head<3>() * tmpDy(3)) + // .finished()); + srf_pt->set_jacobian(botsqr * + (tmpD.block<3, 2>(0, 0) * tmp(3, 0) - + tmp.block<3, 1>(0, 0) * tmpD.block<1, 2>(3, 0))); return; } diff --git a/Bembel/src/Geometry/SurfacePoint.hpp b/Bembel/src/Geometry/SurfacePoint.hpp index 6748d8942..f26294006 100644 --- a/Bembel/src/Geometry/SurfacePoint.hpp +++ b/Bembel/src/Geometry/SurfacePoint.hpp @@ -56,6 +56,16 @@ class SurfacePoint { double get_surface_measure() { return get_normal().norm(); } const double get_surface_measure() const { return get_normal().norm(); } + // buffer + std::vector> + &get_buffer() { + return buffer_; + } + const std::vector> + &get_buffer() const { + return buffer_; + } + // setter void set_xi(const Eigen::Vector2d &xi) { xi_ = xi; } void set_w(const double w) { w_ = w; } @@ -69,6 +79,8 @@ class SurfacePoint { double w_; Eigen::Vector3d f_; Eigen::Matrix jacobian_; + std::vector> + buffer_; }; typedef std::vector> diff --git a/tests/test_SurfacePointUpdate.cpp b/tests/test_SurfacePointUpdate.cpp index 4ab20ca09..3de6d9362 100644 --- a/tests/test_SurfacePointUpdate.cpp +++ b/tests/test_SurfacePointUpdate.cpp @@ -33,13 +33,8 @@ int main() { assert(std::abs(srf_pt.get_w() - 3.1415) < tol); assert((srf_pt.get_f().head<2>() - pt).norm() < tol); assert(std::abs(srf_pt.get_f()(2) < tol)); - assert((srf_pt.get_f() - geometry.get_geometry()[0].eval(pt)).norm() < - tol); assert((srf_pt.get_f_dx() - Eigen::Vector3d(1., 0., 0.)).norm() < tol); assert((srf_pt.get_f_dy() - Eigen::Vector3d(0, 1., 0.)).norm() < tol); - assert( - (srf_pt.get_jacobian() - geometry.get_geometry()[0].evalJacobian(pt)) - .norm() < tol); } } return 0; From 95db07c9ed59d9a82cfbf903f214f87a4fd705aa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20D=C3=B6lz?= Date: Mon, 13 Feb 2023 22:44:27 +0100 Subject: [PATCH 11/20] move memory allocation in updateSurfacePoint to memory allocation in SurfacePoint. Runtimes are comparable to master branch (that is, before start of refactoring) or quicker --- .gitignore | 5 +- .vscode/settings.json | 6 - Bembel/src/Geometry/Patch.hpp | 166 +++++++++++++++------------ Bembel/src/Geometry/SurfacePoint.hpp | 35 ++++-- Bembel/src/Spline/Basis.hpp | 20 ++-- Bembel/src/Spline/Bernstein.hpp | 8 ++ 6 files changed, 139 insertions(+), 101 deletions(-) delete mode 100644 .vscode/settings.json diff --git a/.gitignore b/.gitignore index 2aca7bcf9..a55321fc9 100644 --- a/.gitignore +++ b/.gitignore @@ -6,4 +6,7 @@ eigen3 *.out # Important for a nice gh-pages branch node_modules/ -package-lock.json \ No newline at end of file +package-lock.json +# dev +vtune +.vscode diff --git a/.vscode/settings.json b/.vscode/settings.json deleted file mode 100644 index bde692edd..000000000 --- a/.vscode/settings.json +++ /dev/null @@ -1,6 +0,0 @@ -{ - "files.associations": { - "quadrature": "cpp", - "spline": "cpp" - } -} \ No newline at end of file diff --git a/Bembel/src/Geometry/Patch.hpp b/Bembel/src/Geometry/Patch.hpp index da8dbbe1f..d1b7c0d61 100644 --- a/Bembel/src/Geometry/Patch.hpp +++ b/Bembel/src/Geometry/Patch.hpp @@ -90,22 +90,32 @@ class Patch { Spl::Rescale(reference_point(1), unique_knots_y_[y_location], unique_knots_y_[y_location + 1]); - Eigen::Vector4d tmp = Eigen::Vector4d::Zero(); - Eigen::VectorXd xbasis = Bembel::Basis::ShapeFunctionHandler::evalBasis( - polynomial_degree_x_ - 1, scaledx); - Eigen::VectorXd ybasis = Bembel::Basis::ShapeFunctionHandler::evalBasis( - polynomial_degree_y_ - 1, scaledy); + double *xbasis = new double[polynomial_degree_x_]; + double *ybasis = new double[polynomial_degree_y_]; + + Bembel::Basis::ShapeFunctionHandler::evalBasis(polynomial_degree_x_ - 1, + xbasis, scaledx); + Bembel::Basis::ShapeFunctionHandler::evalBasis(polynomial_degree_y_ - 1, + ybasis, scaledy); + + double tmp[4] = {0., 0., 0., 0.}; for (int i = 0; i < polynomial_degree_x_; i++) { for (int j = 0; j < polynomial_degree_y_; j++) { - const double tpbasisval = xbasis(i) * ybasis(j); + const double tpbasisval = xbasis[i] * ybasis[j]; const int accs = 4 * (numy * (polynomial_degree_x_ * x_location + i) + polynomial_degree_y_ * y_location + j); - for (int k = 0; k < 4; k++) tmp(k) += data_[accs + k] * tpbasisval; + for (int k = 0; k < 4; k++) tmp[k] += data_[accs + k] * tpbasisval; } } - return tmp.head<3>() / tmp[3]; + delete[] xbasis; + delete[] ybasis; + + Eigen::Vector3d out(tmp[0], tmp[1], tmp[2]); + // Rescaling by the NRBS weight, i.e. projection to 3D from 4D hom + + return out / tmp[3]; } Eigen::Matrix evalJacobian( @@ -122,41 +132,55 @@ class Patch { Spl::Rescale(reference_point(1), unique_knots_y_[y_location], unique_knots_y_[y_location + 1]); - Eigen::Vector4d tmp = Eigen::Vector4d::Zero(); - Eigen::Vector4d tmpDx = Eigen::Vector4d::Zero(); - Eigen::Vector4d tmpDy = Eigen::Vector4d::Zero(); - Eigen::VectorXd xbasis = Bembel::Basis::ShapeFunctionHandler::evalBasis( - polynomial_degree_x_ - 1, scaledx); - Eigen::VectorXd ybasis = Bembel::Basis::ShapeFunctionHandler::evalBasis( - polynomial_degree_y_ - 1, scaledy); - Eigen::VectorXd xbasisD = Bembel::Basis::ShapeFunctionHandler::evalDerBasis( - polynomial_degree_x_ - 1, scaledx); - Eigen::VectorXd ybasisD = Bembel::Basis::ShapeFunctionHandler::evalDerBasis( - polynomial_degree_y_ - 1, scaledy); + double *xbasis = new double[polynomial_degree_x_]; + double *ybasis = new double[polynomial_degree_y_]; + double *xbasisD = new double[polynomial_degree_x_]; + double *ybasisD = new double[polynomial_degree_y_]; + + Bembel::Basis::ShapeFunctionHandler::evalBasis(polynomial_degree_x_ - 1, + xbasis, scaledx); + Bembel::Basis::ShapeFunctionHandler::evalBasis(polynomial_degree_y_ - 1, + ybasis, scaledy); + Bembel::Basis::ShapeFunctionHandler::evalDerBasis(polynomial_degree_x_ - 1, + xbasisD, scaledx); + Bembel::Basis::ShapeFunctionHandler::evalDerBasis(polynomial_degree_y_ - 1, + ybasisD, scaledy); + + double tmp[4] = {0., 0., 0., 0.}; + double tmpDx[4] = {0., 0., 0., 0.}; + double tmpDy[4] = {0., 0., 0., 0.}; for (int i = 0; i < polynomial_degree_x_; i++) { for (int j = 0; j < polynomial_degree_y_; j++) { - const double tpbasisval = xbasis(i) * ybasis(j); - const double tpbasisvalDx = xbasisD(i) * ybasis(j); - const double tpbasisvalDy = xbasis(i) * ybasisD(j); + const double tpbasisval = xbasis[i] * ybasis[j]; + const double tpbasisvalDx = xbasisD[i] * ybasis[j]; + const double tpbasisvalDy = xbasis[i] * ybasisD[j]; const int accs = 4 * (numy * (polynomial_degree_x_ * x_location + i) + polynomial_degree_y_ * y_location + j); // Here I add up the values of the basis functions in the dc // basis for (int k = 0; k < 4; k++) { - tmp(k) += data_[accs + k] * tpbasisval; - tmpDx(k) += data_[accs + k] * tpbasisvalDx; - tmpDy(k) += data_[accs + k] * tpbasisvalDy; + tmp[k] += data_[accs + k] * tpbasisval; + tmpDx[k] += data_[accs + k] * tpbasisvalDx; + tmpDy[k] += data_[accs + k] * tpbasisvalDy; } } } - double botsqr = 1. / (tmp(3) * tmp(3)); - return botsqr * (Eigen::Matrix() - << tmpDx.head<3>() * tmp(3) - tmp.head<3>() * tmpDx(3), - tmpDy.head<3>() * tmp(3) - tmp.head<3>() * tmpDy(3)) - .finished(); + delete[] xbasis; + delete[] ybasis; + delete[] xbasisD; + delete[] ybasisD; + + double bot = 1. / (tmp[3] * tmp[3]); + Eigen::Matrix out; + for (int k = 0; k < 3; k++) { + out(k, 0) = (tmpDx[k] * tmp[3] - tmp[k] * tmpDx[3]) * bot; + out(k, 1) = (tmpDy[k] * tmp[3] - tmp[k] * tmpDy[3]) * bot; + } + + return out; } inline Eigen::Matrix evalNormal( @@ -190,67 +214,61 @@ class Patch { const double scaledy = Spl::Rescale(ref_pt(1), unique_knots_y_[y_location], unique_knots_y_[y_location + 1]); - // resize operations do nothing if the size does not change - srf_pt->get_buffer().resize(6); - srf_pt->get_buffer()[0].resize(4, 1); - srf_pt->get_buffer()[1].resize(4, 2); - srf_pt->get_buffer()[2].resize(polynomial_degree_x_, 1); - srf_pt->get_buffer()[3].resize(polynomial_degree_y_, 1); - srf_pt->get_buffer()[4].resize(polynomial_degree_x_, 1); - srf_pt->get_buffer()[5].resize(polynomial_degree_y_, 1); - - // improve readability - Eigen::MatrixXd &tmp = srf_pt->get_buffer()[0]; - Eigen::MatrixXd &tmpD = srf_pt->get_buffer()[1]; - Eigen::MatrixXd &xbasis = srf_pt->get_buffer()[2]; - Eigen::MatrixXd &ybasis = srf_pt->get_buffer()[3]; - Eigen::MatrixXd &xbasisD = srf_pt->get_buffer()[4]; - Eigen::MatrixXd &ybasisD = srf_pt->get_buffer()[5]; - - tmp.setZero(); - tmpD.setZero(); - - xbasis.col(0) = Bembel::Basis::ShapeFunctionHandler::evalBasis( - polynomial_degree_x_ - 1, scaledx); - ybasis.col(0) = Bembel::Basis::ShapeFunctionHandler::evalBasis( - polynomial_degree_y_ - 1, scaledy); - xbasisD.col(0) = Bembel::Basis::ShapeFunctionHandler::evalDerBasis( - polynomial_degree_x_ - 1, scaledx); - ybasisD.col(0) = Bembel::Basis::ShapeFunctionHandler::evalDerBasis( - polynomial_degree_y_ - 1, scaledy); + // allocating memory is expensive, use the buffer in SurfacePoint + srf_pt->allocate_buffer(2 * (polynomial_degree_x_ + polynomial_degree_y_) + + 12); + double *buffer = srf_pt->get_buffer(); + std::memset(buffer, 0, 12 * sizeof(double)); + + double *tmp = buffer; + double *tmpDx = tmp + 4; + double *tmpDy = tmpDx + 4; + double *xbasis = tmpDy + 4; + double *ybasis = xbasis + polynomial_degree_x_; + double *xbasisD = ybasis + polynomial_degree_y_; + double *ybasisD = xbasisD + polynomial_degree_x_; + + Bembel::Basis::ShapeFunctionHandler::evalBasis(polynomial_degree_x_ - 1, + xbasis, scaledx); + Bembel::Basis::ShapeFunctionHandler::evalBasis(polynomial_degree_y_ - 1, + ybasis, scaledy); + Bembel::Basis::ShapeFunctionHandler::evalDerBasis(polynomial_degree_x_ - 1, + xbasisD, scaledx); + Bembel::Basis::ShapeFunctionHandler::evalDerBasis(polynomial_degree_y_ - 1, + ybasisD, scaledy); for (int i = 0; i < polynomial_degree_x_; ++i) { for (int j = 0; j < polynomial_degree_y_; ++j) { - const double tpbasisval = xbasis(i) * ybasis(j); - const double tpbasisvalDx = xbasisD(i) * ybasis(j); - const double tpbasisvalDy = xbasis(i) * ybasisD(j); + const double tpbasisval = xbasis[i] * ybasis[j]; + const double tpbasisvalDx = xbasisD[i] * ybasis[j]; + const double tpbasisvalDy = xbasis[i] * ybasisD[j]; const int accs = 4 * (numy * (polynomial_degree_x_ * x_location + i) + polynomial_degree_y_ * y_location + j); // Here I add up the values of the basis functions in the dc // basis for (int k = 0; k < 4; ++k) { - tmp(k, 0) += data_[accs + k] * tpbasisval; - tmpD(k, 0) += data_[accs + k] * tpbasisvalDx; - tmpD(k, 1) += data_[accs + k] * tpbasisvalDy; + tmp[k] += data_[accs + k] * tpbasisval; + tmpDx[k] += data_[accs + k] * tpbasisvalDx; + tmpDy[k] += data_[accs + k] * tpbasisvalDy; } } } - const double bot = 1. / tmp(3, 0); + const double bot = 1. / tmp[3]; const double botsqr = bot * bot; + Eigen::Vector3d tmpvec(tmp[0], tmp[1], tmp[2]); + Eigen::Vector3d tmpDxvec(tmpDx[0], tmpDx[1], tmpDx[2]); + Eigen::Vector3d tmpDyvec(tmpDy[0], tmpDy[1], tmpDy[2]); + srf_pt->set_xi(xi); srf_pt->set_w(w); - srf_pt->set_f(bot * tmp.block<3, 1>(0, 0)); - // srf_pt->set_jacobian( - // botsqr * (Eigen::Matrix() - // << tmpDx.head<3>() * tmp(3) - tmp.head<3>() * tmpDx(3), - // tmpDy.head<3>() * tmp(3) - tmp.head<3>() * tmpDy(3)) - // .finished()); - srf_pt->set_jacobian(botsqr * - (tmpD.block<3, 2>(0, 0) * tmp(3, 0) - - tmp.block<3, 1>(0, 0) * tmpD.block<1, 2>(3, 0))); + srf_pt->set_f(bot * tmpvec); + srf_pt->set_jacobian(botsqr * (Eigen::Matrix() + << tmpDxvec * tmp[3] - tmpvec * tmpDx[3], + tmpDyvec * tmp[3] - tmpvec * tmpDy[3]) + .finished()); return; } diff --git a/Bembel/src/Geometry/SurfacePoint.hpp b/Bembel/src/Geometry/SurfacePoint.hpp index f26294006..d48257da2 100644 --- a/Bembel/src/Geometry/SurfacePoint.hpp +++ b/Bembel/src/Geometry/SurfacePoint.hpp @@ -20,6 +20,10 @@ **/ class SurfacePoint { public: + ~SurfacePoint() { + if (buffer_is_allocated()) delete[] buffer_; + } + // quadrature point Eigen::Vector2d get_xi() { return xi_; } const Eigen::Vector2d get_xi() const { return xi_; } @@ -57,14 +61,10 @@ class SurfacePoint { const double get_surface_measure() const { return get_normal().norm(); } // buffer - std::vector> - &get_buffer() { - return buffer_; - } - const std::vector> - &get_buffer() const { - return buffer_; - } + bool buffer_is_allocated() { return buffer_is_alloced_; } + const bool buffer_is_allocated() const { return buffer_is_alloced_; } + double *get_buffer() { return buffer_; } + const double *get_buffer() const { return buffer_; } // setter void set_xi(const Eigen::Vector2d &xi) { xi_ = xi; } @@ -73,14 +73,29 @@ class SurfacePoint { void set_jacobian(const Eigen::Matrix &jacobian) { jacobian_ = jacobian; } + void allocate_buffer(const int n) { + assert(n > 0); + if (n != buffer_size_) { + if (buffer_is_allocated()) { + delete[] buffer_; + buffer_is_alloced_ = false; + } + // TODO(Felix) Do not use variable-length arrays in accordance to Google + // Style guide + buffer_ = new double[n]; + buffer_size_ = n; + buffer_is_alloced_ = true; + } + } private: Eigen::Vector2d xi_; double w_; Eigen::Vector3d f_; Eigen::Matrix jacobian_; - std::vector> - buffer_; + bool buffer_is_alloced_ = false; + int buffer_size_ = 0; + double *buffer_; }; typedef std::vector> diff --git a/Bembel/src/Spline/Basis.hpp b/Bembel/src/Spline/Basis.hpp index 5e4b88d03..c86116e51 100644 --- a/Bembel/src/Spline/Basis.hpp +++ b/Bembel/src/Spline/Basis.hpp @@ -157,8 +157,8 @@ void Phi_times_Phi_(Eigen::Matrix *c, for (int iy = 0; iy < I; iy++) for (int ix = 0; ix < I; ix++) b[iy * I + ix] = X[ix] * Y[iy]; - for (int i = 0; i < (I * I); i++) - for (int j = 0; j < (I * I); j++) (*c)(i, j) += a[i] * b[j]; + for (int j = 0; j < (I * I); j++) + for (int i = 0; i < (I * I); i++) (*c)(i, j) += a[i] * b[j]; return; } @@ -187,14 +187,14 @@ void Div_Phi_times_Div_Phi_( phiphi_dx_

(&b_dx, 1., eta); phiphi_dy_

(&b_dy, 1., eta); - for (int i = 0; i < I2; ++i) - for (int j = 0; j < I2; ++j) (*c)(i, j) += a_dx[i] * b_dx[j]; - for (int i = 0; i < I2; ++i) - for (int j = 0; j < I2; ++j) (*c)(i, j + I2) += a_dx[i] * b_dy[j]; - for (int i = 0; i < I2; ++i) - for (int j = 0; j < I2; ++j) (*c)(i + I2, j) += a_dy[i] * b_dx[j]; - for (int i = 0; i < I2; ++i) - for (int j = 0; j < I2; ++j) (*c)(i + I2, j + I2) += a_dy[i] * b_dy[j]; + for (int j = 0; j < I2; ++j) + for (int i = 0; i < I2; ++i) (*c)(i, j) += a_dx[i] * b_dx[j]; + for (int j = 0; j < I2; ++j) + for (int i = 0; i < I2; ++i) (*c)(i, j + I2) += a_dx[i] * b_dy[j]; + for (int j = 0; j < I2; ++j) + for (int i = 0; i < I2; ++i) (*c)(i + I2, j) += a_dy[i] * b_dx[j]; + for (int j = 0; j < I2; ++j) + for (int i = 0; i < I2; ++i) (*c)(i + I2, j + I2) += a_dy[i] * b_dy[j]; return; } diff --git a/Bembel/src/Spline/Bernstein.hpp b/Bembel/src/Spline/Bernstein.hpp index bba57155e..0bc440ed1 100644 --- a/Bembel/src/Spline/Bernstein.hpp +++ b/Bembel/src/Spline/Bernstein.hpp @@ -52,12 +52,16 @@ inline constexpr double Bernstein(double evaluation_point) noexcept { template class HiddenBernsteinClass { public: + // think twice when replacing double* by something else. All other attempts + // were slower static inline void EvalBasis(double *in, const double evaluation_point) noexcept { in[N] = Bernstein(evaluation_point); HiddenBernsteinClass::EvalBasis(in, evaluation_point); return; } + // think twice when replacing double* by something else. All other attempts + // were slower static inline void EvalDerBasis(double *in, const double evaluation_point) noexcept { in[N] = (P + 1) * (Bernstein(evaluation_point) - @@ -108,11 +112,15 @@ class HiddenBernsteinClass<-1, P> { /// Evaluation Routines //////////////////////////////////////////////////////////////////////////////// +// think twice when replacing double* by something else. All other attempts +// were slower template void EvalBernsteinBasis(double *in, const double evaluation_point) noexcept { HiddenBernsteinClass::EvalBasis(in, evaluation_point); } +// think twice when replacing double* by something else. All other attempts +// were slower template void EvalBernsteinDerBasis(double *in, const double evaluation_point) noexcept { in[P] = P * Bernstein

(evaluation_point); From 210bdb7201c0303ec6a194104c6817cae8166374 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20D=C3=B6lz?= Date: Mon, 18 Mar 2024 15:28:38 +0100 Subject: [PATCH 12/20] add patch number to surface point --- Bembel/src/AnsatzSpace/FunctionEvaluator.hpp | 2 +- Bembel/src/AnsatzSpace/SuperSpace.hpp | 5 +++-- Bembel/src/Geometry/Patch.hpp | 6 ++++-- Bembel/src/Geometry/SurfacePoint.hpp | 6 ++++++ tests/test_SurfacePointUpdate.cpp | 2 +- 5 files changed, 15 insertions(+), 6 deletions(-) diff --git a/Bembel/src/AnsatzSpace/FunctionEvaluator.hpp b/Bembel/src/AnsatzSpace/FunctionEvaluator.hpp index 9f6f7316d..116e2ae19 100644 --- a/Bembel/src/AnsatzSpace/FunctionEvaluator.hpp +++ b/Bembel/src/AnsatzSpace/FunctionEvaluator.hpp @@ -100,7 +100,7 @@ class FunctionEvaluator { SurfacePoint sp; ansatz_space_.get_superspace().get_geometry()[patch].updateSurfacePoint( - &sp, ref_point, 1, element.mapToReferenceElement(ref_point)); + &sp, patch, ref_point, 1, element.mapToReferenceElement(ref_point)); return evaluate(element, sp); } Eigen::Matrix< diff --git a/Bembel/src/AnsatzSpace/SuperSpace.hpp b/Bembel/src/AnsatzSpace/SuperSpace.hpp index a0be0ed00..fb81e6403 100644 --- a/Bembel/src/AnsatzSpace/SuperSpace.hpp +++ b/Bembel/src/AnsatzSpace/SuperSpace.hpp @@ -107,7 +107,8 @@ struct SuperSpace { void map2surface(const ElementTreeNode& e, const Eigen::Vector2d& xi, double w, SurfacePoint* surf_pt) const { Eigen::Vector2d st = e.llc_ + e.get_h() * xi; - mesh_->get_geometry()[e.patch_].updateSurfacePoint(surf_pt, st, w, xi); + mesh_->get_geometry()[e.patch_].updateSurfacePoint(surf_pt, e.patch_, st, w, + xi); return; } ////////////////////////////////////////////////////////////////////////////// @@ -156,7 +157,7 @@ struct SuperSpace { // inner product of surface curls of any two basis functions for (int j = 0; j < polynomial_degree_plus_one_squared; ++j) for (int i = 0; i < polynomial_degree_plus_one_squared; ++i) - (*intval)(j * polynomial_degree_plus_one_squared + i) += + (*intval)(j* polynomial_degree_plus_one_squared + i) += w * s_curl.col(i).dot(t_curl.col(j)); } diff --git a/Bembel/src/Geometry/Patch.hpp b/Bembel/src/Geometry/Patch.hpp index d1b7c0d61..a2ed9363b 100644 --- a/Bembel/src/Geometry/Patch.hpp +++ b/Bembel/src/Geometry/Patch.hpp @@ -202,8 +202,9 @@ class Patch { // This is a combination of eval und evalJacobian, to avoid duplication of // work. See SurfacePoint.hpp - void updateSurfacePoint(SurfacePoint *srf_pt, const Eigen::Vector2d &ref_pt, - double w, const Eigen::Vector2d &xi) const { + void updateSurfacePoint(SurfacePoint *srf_pt, const int patch, + const Eigen::Vector2d &ref_pt, const double w, + const Eigen::Vector2d &xi) const { const int x_location = Spl::FindLocationInKnotVector(ref_pt(0), unique_knots_x_); const int y_location = @@ -262,6 +263,7 @@ class Patch { Eigen::Vector3d tmpDxvec(tmpDx[0], tmpDx[1], tmpDx[2]); Eigen::Vector3d tmpDyvec(tmpDy[0], tmpDy[1], tmpDy[2]); + srf_pt->set_patch(patch); srf_pt->set_xi(xi); srf_pt->set_w(w); srf_pt->set_f(bot * tmpvec); diff --git a/Bembel/src/Geometry/SurfacePoint.hpp b/Bembel/src/Geometry/SurfacePoint.hpp index d48257da2..d6f5d35ee 100644 --- a/Bembel/src/Geometry/SurfacePoint.hpp +++ b/Bembel/src/Geometry/SurfacePoint.hpp @@ -24,6 +24,10 @@ class SurfacePoint { if (buffer_is_allocated()) delete[] buffer_; } + // patch + int get_patch() { return patch_; } + const int get_patch() const { return patch_; } + // quadrature point Eigen::Vector2d get_xi() { return xi_; } const Eigen::Vector2d get_xi() const { return xi_; } @@ -67,6 +71,7 @@ class SurfacePoint { const double *get_buffer() const { return buffer_; } // setter + void set_patch(const int patch) { patch_ = patch; } void set_xi(const Eigen::Vector2d &xi) { xi_ = xi; } void set_w(const double w) { w_ = w; } void set_f(const Eigen::Vector3d &f) { f_ = f; } @@ -89,6 +94,7 @@ class SurfacePoint { } private: + int patch_; Eigen::Vector2d xi_; double w_; Eigen::Vector3d f_; diff --git a/tests/test_SurfacePointUpdate.cpp b/tests/test_SurfacePointUpdate.cpp index 3de6d9362..6e8c154d0 100644 --- a/tests/test_SurfacePointUpdate.cpp +++ b/tests/test_SurfacePointUpdate.cpp @@ -28,7 +28,7 @@ int main() { for (auto y : Test::Constants::eq_points) { Eigen::Vector2d pt = Eigen::Vector2d(x, y); SurfacePoint srf_pt; - geometry.get_geometry()[0].updateSurfacePoint(&srf_pt, pt, 3.1415, pt); + geometry.get_geometry()[0].updateSurfacePoint(&srf_pt, 0, pt, 3.1415, pt); assert((srf_pt.get_xi() - pt).norm() < tol); assert(std::abs(srf_pt.get_w() - 3.1415) < tol); assert((srf_pt.get_f().head<2>() - pt).norm() < tol); From 297586273be36b137bc20c20dadc932de429635c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20D=C3=B6lz?= Date: Mon, 18 Mar 2024 16:20:05 +0100 Subject: [PATCH 13/20] debugging --- Bembel/LinearForm | 1 - Bembel/src/LinearForm/DirichletTrace.hpp | 5 +- .../DirichletTraceOnReferenceDomain.hpp | 43 +++---- .../DiscreteLinearFormOnReferenceDomain.hpp | 107 ----------------- Bembel/src/LinearForm/NeumannTrace.hpp | 22 +--- Bembel/src/LinearForm/TangentialTrace.hpp | 25 +--- Bembel/src/util/surfaceL2error.hpp | 4 +- examples/CMakeLists.txt | 1 + .../DiscreteLinearFormOnReferenceDomain.cpp | 111 ++++++++++++++++++ 9 files changed, 137 insertions(+), 182 deletions(-) delete mode 100644 Bembel/src/LinearForm/DiscreteLinearFormOnReferenceDomain.hpp create mode 100644 examples/DiscreteLinearFormOnReferenceDomain.cpp diff --git a/Bembel/LinearForm b/Bembel/LinearForm index 7ea6cafa1..8e7229216 100644 --- a/Bembel/LinearForm +++ b/Bembel/LinearForm @@ -22,7 +22,6 @@ #include "Quadrature" #include "src/LinearForm/LinearForm.hpp" -#include "src/LinearForm/DiscreteLinearFormOnReferenceDomain.hpp" #include "src/LinearForm/DirichletTrace.hpp" #include "src/LinearForm/DirichletTraceOnReferenceDomain.hpp" diff --git a/Bembel/src/LinearForm/DirichletTrace.hpp b/Bembel/src/LinearForm/DirichletTrace.hpp index fbc7c7186..d8d33ccbe 100644 --- a/Bembel/src/LinearForm/DirichletTrace.hpp +++ b/Bembel/src/LinearForm/DirichletTrace.hpp @@ -38,10 +38,9 @@ class DirichletTrace : public LinearFormBase, Scalar> { void evaluateIntegrand_impl( const T &super_space, const SurfacePoint &p, Eigen::Matrix *intval) const { - // compute surface measures from tangential derivatives - auto x_kappa = p.get_f_dx().cross(p.get_f_dy()).norm(); + std::cout << p.get_xi() << std::endl; // integrand without basis functions - auto integrand = function_(p.get_f()) * p.get_surface_measure() * p.get_w(); + Scalar integrand = function_(p.get_f()) * p.get_surface_measure() * p.get_w(); // multiply basis functions with integrand super_space.addScaledBasis(intval, integrand, p.get_xi()); return; diff --git a/Bembel/src/LinearForm/DirichletTraceOnReferenceDomain.hpp b/Bembel/src/LinearForm/DirichletTraceOnReferenceDomain.hpp index 9f417ae9b..f612e0123 100644 --- a/Bembel/src/LinearForm/DirichletTraceOnReferenceDomain.hpp +++ b/Bembel/src/LinearForm/DirichletTraceOnReferenceDomain.hpp @@ -1,12 +1,15 @@ // This file is part of Bembel, the higher order C++ boundary element library. +// +// Copyright (C) 2024 see +// // It was written as part of a cooperation of J. Doelz, H. Harbrecht, S. Kurz, // M. Multerer, S. Schoeps, and F. Wolf at Technische Universitaet Darmstadt, // Universitaet Basel, and Universita della Svizzera italiana, Lugano. This // source code is subject to the GNU General Public License version 3 and // provided WITHOUT ANY WARRANTY, see for further // information. -#ifndef BEMBEL_LINEARFORM_DIRICHLETTRACEONREFERENCEDOMAIN_H_ -#define BEMBEL_LINEARFORM_DIRICHLETTRACEONREFERENCEDOMAIN_H_ +#ifndef BEMBEL_SRC_LINEARFORM_DIRICHLETTRACEONREFERENCEDOMAIN_H_ +#define BEMBEL_SRC_LINEARFORM_DIRICHLETTRACEONREFERENCEDOMAIN_H_ namespace Bembel { @@ -34,39 +37,21 @@ class DirichletTraceOnReferenceDomain function_ = function; } template - void evaluateIntegrand_impl(const T &super_space, const SurfacePoint &p, - Eigen::Matrix *intval, - int patch) const { - auto polynomial_degree = super_space.get_polynomial_degree(); - auto polynomial_degree_plus_one_squared = - (polynomial_degree + 1) * (polynomial_degree + 1); - - // get evaluation points on unit square - auto s = p.segment<2>(0); - - // get quadrature weights - auto ws = p(2); - - // get points on geometry and tangential derivatives - auto x_f = p.segment<3>(3); - auto x_f_dx = p.segment<3>(6); - auto x_f_dy = p.segment<3>(9); - - // compute surface measures from tangential derivatives - auto x_kappa = x_f_dx.cross(x_f_dy).norm(); - + void evaluateIntegrand_impl( + const T &super_space, const SurfacePoint &p, + Eigen::Matrix *intval) const { // integrand without basis functions - auto integrand = function_(patch, s) * x_kappa * ws; - + //std::cout << p.get_xi() << std::endl; + Scalar integrand = function_(p.get_patch(), p.get_xi()) * + p.get_surface_measure() * p.get_w(); // multiply basis functions with integrand - super_space.addScaledBasis(intval, integrand, s); - + super_space.addScaledBasis(intval, integrand, p.get_xi()); return; - }; + } private: std::function function_; }; } // namespace Bembel -#endif +#endif // BEMBEL_SRC_LINEARFORM_DIRICHLETTRACEONREFERENCEDOMAIN_H_ diff --git a/Bembel/src/LinearForm/DiscreteLinearFormOnReferenceDomain.hpp b/Bembel/src/LinearForm/DiscreteLinearFormOnReferenceDomain.hpp deleted file mode 100644 index 8f9a1536c..000000000 --- a/Bembel/src/LinearForm/DiscreteLinearFormOnReferenceDomain.hpp +++ /dev/null @@ -1,107 +0,0 @@ -// This file is part of Bembel, the higher order C++ boundary element library. -// It was written as part of a cooperation of J. Doelz, H. Harbrecht, S. Kurz, -// M. Multerer, S. Schoeps, and F. Wolf at Technische Universitaet Darmstadt, -// Universitaet Basel, and Universita della Svizzera italiana, Lugano. This -// source code is subject to the GNU General Public License version 3 and -// provided WITHOUT ANY WARRANTY, see for further -// information. -#ifndef BEMBEL_LINEARFORM_DISCRETELINEARFORMONREFERENCEDOMAIN_H_ -#define BEMBEL_LINEARFORM_DISCRETELINEARFORMONREFERENCEDOMAIN_H_ - -namespace Bembel { -/** - * \ingroup LinearForm - * \brief This class, given a LinearForm with defined traits, takes care of the - * assembly of the right hand side. - */ -template -class DiscreteLinearFormOnReferenceDomain { - public: - ////////////////////////////////////////////////////////////////////////////// - // constructors - ////////////////////////////////////////////////////////////////////////////// - DiscreteLinearFormOnReferenceDomain() {} - DiscreteLinearFormOnReferenceDomain(const AnsatzSpace &ansatz_space) { - init_DiscreteLinearFormOnReferenceDomain(ansatz_space); - } - ////////////////////////////////////////////////////////////////////////////// - // init_DiscreteLinearFormOnReferenceDomain - ////////////////////////////////////////////////////////////////////////////// - void init_DiscreteLinearFormOnReferenceDomain( - const AnsatzSpace &ansatz_space) { - ansatz_space_ = ansatz_space; - deg_ = ansatz_space_.get_polynomial_degree() + 1; - return; - } - ////////////////////////////////////////////////////////////////////////////// - // compute - ////////////////////////////////////////////////////////////////////////////// - /** - * \todo Add inline commentary - **/ - void compute() { - GaussSquare GS; - const auto &Q = GS[deg_]; - SurfacePoint qp; - const auto &super_space = ansatz_space_.get_superspace(); - const auto &element_tree = super_space.get_mesh().get_element_tree(); - auto number_of_elements = element_tree.get_number_of_elements(); - auto polynomial_degree = super_space.get_polynomial_degree(); - auto polynomial_degree_plus_one_squared = - (polynomial_degree + 1) * (polynomial_degree + 1); - const auto function_space_dimension = - getFunctionSpaceVectorDimension::Form>(); - Eigen::Matrix::Scalar, Eigen::Dynamic, - function_space_dimension> - intval(polynomial_degree_plus_one_squared, function_space_dimension); - Eigen::Matrix::Scalar, Eigen::Dynamic, - function_space_dimension> - disc_lf_matrix(polynomial_degree_plus_one_squared * number_of_elements, - function_space_dimension); - disc_lf_matrix.setZero(); - for (auto element = element_tree.cpbegin(); element != element_tree.cpend(); - ++element) { - intval.setZero(); - for (auto i = 0; i < Q.w_.size(); ++i) { - super_space.map2surface(*element, Q.xi_.col(i), - element->get_h() * Q.w_(i), &qp); - lf_.evaluateIntegrand_impl(super_space, qp, &intval, element->patch_); - } - disc_lf_matrix.block(polynomial_degree_plus_one_squared * element->id_, 0, - polynomial_degree_plus_one_squared, - function_space_dimension) = intval; - } - disc_lf_ = - ansatz_space_.get_transformation_matrix().transpose() * - Eigen::Map::Scalar, - Eigen::Dynamic, 1>>( - disc_lf_matrix.data(), - disc_lf_matrix.rows() * disc_lf_matrix.cols()); - return; - } - ////////////////////////////////////////////////////////////////////////////// - // setter - ////////////////////////////////////////////////////////////////////////////// - void set_degree(const int °) { deg_ = deg; } - ////////////////////////////////////////////////////////////////////////////// - // getter - ////////////////////////////////////////////////////////////////////////////// - Derived &get_linear_form() { return lf_; } - const Eigen::Matrix::Scalar, - Eigen::Dynamic, 1> - &get_discrete_linear_form() const { - return disc_lf_; - } - ////////////////////////////////////////////////////////////////////////////// - // private member variables - ////////////////////////////////////////////////////////////////////////////// - private: - int deg_; - Derived lf_; - Eigen::Matrix::Scalar, Eigen::Dynamic, 1> - disc_lf_; - AnsatzSpace ansatz_space_; -}; - -} // namespace Bembel -#endif diff --git a/Bembel/src/LinearForm/NeumannTrace.hpp b/Bembel/src/LinearForm/NeumannTrace.hpp index 1a14a2430..2798298ed 100644 --- a/Bembel/src/LinearForm/NeumannTrace.hpp +++ b/Bembel/src/LinearForm/NeumannTrace.hpp @@ -40,30 +40,12 @@ class NeumannTrace : public LinearFormBase, Scalar> { void evaluateIntegrand_impl( const T &super_space, const SurfacePoint &p, Eigen::Matrix *intval) const { - auto polynomial_degree = super_space.get_polynomial_degree(); - auto polynomial_degree_plus_one_squared = - (polynomial_degree + 1) * (polynomial_degree + 1); - - // get evaluation points on unit square - auto s = p.segment<2>(0); - - // get quadrature weights - auto ws = p(2); - - // get points on geometry and tangential derivatives - auto x_f = p.segment<3>(3); - auto x_f_dx = p.segment<3>(6); - auto x_f_dy = p.segment<3>(9); - - // compute (unnormalized) surface normal from tangential derivatives - auto x_f_n = x_f_dx.cross(x_f_dy); - // integrand without basis functions // dot: adjoint in first variable - auto integrand = x_f_n.dot(function_(x_f)) * ws; + auto integrand = p.get_normal().dot(function_(p.get_f)) * p.get_w(); // multiply basis functions with integrand - super_space.addScaledBasis(intval, integrand, s); + super_space.addScaledBasis(intval, integrand, p.get_xi()); return; } diff --git a/Bembel/src/LinearForm/TangentialTrace.hpp b/Bembel/src/LinearForm/TangentialTrace.hpp index fc79f9e6f..4a212ca9e 100644 --- a/Bembel/src/LinearForm/TangentialTrace.hpp +++ b/Bembel/src/LinearForm/TangentialTrace.hpp @@ -43,31 +43,18 @@ class TangentialTrace : public LinearFormBase, Scalar> { int polynomial_degree_plus_one_squared = (polynomial_degree + 1) * (polynomial_degree + 1); - // get evaluation points on unit square - Eigen::Vector2d s = p.segment<2>(0); - - // get quadrature weights - double ws = p(2); - - // get points on geometry and tangential derivatives - Eigen::Vector3d x_f = p.segment<3>(3); - Eigen::Vector3d x_f_dx = p.segment<3>(6); - Eigen::Vector3d x_f_dy = p.segment<3>(9); - - // compute surface measures from tangential derivatives - Eigen::Vector3d x_n = x_f_dx.cross(x_f_dy).normalized(); - // tangential component + quadrature weights - Eigen::Matrix fun_x_f = function_(x_f); - Eigen::Matrix tangential_component = fun_x_f.cross(x_n) * ws; + Eigen::Matrix fun_x_f = function_(p.get_f()); + Eigen::Matrix tangential_component = + fun_x_f.cross(p.get_unit_normal()) * p.get_w(); // extract tangential component - Scalar component_x = x_f_dx.dot(tangential_component); - Scalar component_y = x_f_dy.dot(tangential_component); + Scalar component_x = p.get_f_dx().dot(tangential_component); + Scalar component_y = p.get_f_dy().dot(tangential_component); // evaluate shape functions Eigen::Matrix phiPhiVec = - super_space.basis(s); + super_space.basis(p.get_xi()); // multiply basis functions with integrand Eigen::Matrix phiPhiMat( diff --git a/Bembel/src/util/surfaceL2error.hpp b/Bembel/src/util/surfaceL2error.hpp index 96f7bc653..3c9ae3012 100644 --- a/Bembel/src/util/surfaceL2error.hpp +++ b/Bembel/src/util/surfaceL2error.hpp @@ -32,9 +32,7 @@ double surfaceL2error(const AnsatzSpace &ansatz_space, for (auto i = 0; i < Q.w_.size(); ++i) { super_space.map2surface(*element, Q.xi_.col(i), Q.w_(i), &qp); // integrand without basis functions - const Scalar val = - longvec.segment(n_shape_fun * element->id_, n_shape_fun).transpose() * - super_space.basis(qp.get_xi()); + const Scalar val = fun_val.evaluate(*element, qp)(0); retval += qp.get_surface_measure() * Q.w_(i) * element->get_h() * element->get_h() * (functor(qp.get_f()) - val / element->get_h()) * diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index ca5a6d78f..4f0e8ebb2 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -11,6 +11,7 @@ configure_file(${CMAKE_SOURCE_DIR}/geo/cube_small.dat ${CMAKE_CURRENT_BINARY_DIR set(CIFILES AnsatzSpace BlockClusterTree + DiscreteLinearFormOnReferenceDomain Geometry MassMatrix LaplaceBeltrami diff --git a/examples/DiscreteLinearFormOnReferenceDomain.cpp b/examples/DiscreteLinearFormOnReferenceDomain.cpp new file mode 100644 index 000000000..4b1d447bf --- /dev/null +++ b/examples/DiscreteLinearFormOnReferenceDomain.cpp @@ -0,0 +1,111 @@ +// This file is part of Bembel, the higher order C++ boundary element library. +// +// Copyright (C) 2024 see +// +// It was written as part of a cooperation of J. Doelz, H. Harbrecht, S. Kurz, +// M. Multerer, S. Schoeps, and F. Wolf at Technische Universitaet Darmstadt, +// Universitaet Basel, and Universita della Svizzera italiana, Lugano. This +// source code is subject to the GNU General Public License version 3 and +// provided WITHOUT ANY WARRANTY, see for further +// information. + +#include +#include +#include +#include +#include +#include + +#include "Bembel/src/util/surfaceL2error.hpp" +#include "examples/Error.hpp" + +int main() { + using namespace Bembel; + using namespace Eigen; + + int polynomial_degree_max = 0; + int refinement_level_max = 1; + + std::function fun = [](const Vector3d &in) { + return in(0); + }; + Geometry geometry("sphere.dat"); + std::cout << "\n" << std::string(60, '=') << "\n"; + // Iterate over polynomial degree. + for (int polynomial_degree = 0; polynomial_degree < polynomial_degree_max + 1; + ++polynomial_degree) { + VectorXd error(refinement_level_max + 1); + // Iterate over refinement levels + for (int refinement_level = 0; refinement_level < refinement_level_max + 1; + ++refinement_level) { + std::cout << "Degree " << polynomial_degree << " Level " + << refinement_level << "\t\t"; + // Build ansatz space + AnsatzSpace ansatz_space(geometry, refinement_level, + polynomial_degree); + + // Set up and compute discrete operator + DiscreteLocalOperator disc_op(ansatz_space); + disc_op.compute(); + DiscreteLinearForm, MassMatrixScalarDisc> disc_lf( + ansatz_space); + disc_lf.get_linear_form().set_function(fun); + disc_lf.compute(); + + // Set up solver + SparseLU, COLAMDOrdering> solver; + solver.analyzePattern(disc_op.get_discrete_operator()); + solver.factorize(disc_op.get_discrete_operator()); + VectorXd x = solver.solve(disc_lf.get_discrete_linear_form()); + + // Set up function on reference domain + FunctionEvaluator evaluator(ansatz_space); + evaluator.set_function(x); + std::function fun_disc = + [&](int p, const Vector2d &x) { + return evaluator.evaluateOnPatch(p, x)(0); + }; + + // Set up linear form on reference domain + DiscreteLinearForm, + MassMatrixScalarDisc> + disc_lf2(ansatz_space); + disc_lf2.get_linear_form().set_function(fun_disc); + disc_lf2.compute(); + + VectorXd x2 = solver.solve(disc_lf2.get_discrete_linear_form()); + + // std::cout << std::endl; + // std::cout << x.head(10).transpose() << std::endl; + // std::cout << x2.head(10).transpose() << std::endl; + // std::cout << (x.head(10).array()/x2.head(10).array()).transpose() << + // std::endl; + std::cout << (disc_lf.get_discrete_linear_form() - + disc_lf2.get_discrete_linear_form()) + .norm() / + disc_lf.get_discrete_linear_form().norm() + << std::endl; + + VTKSurfaceExport writer(geometry, 5); + writer.addDataSet("x", ansatz_space, x); + writer.addDataSet("x2", ansatz_space, x2); + writer.addDataSet("fun_disc", fun_disc); + writer.addDataSet("fun", fun); + writer.writeToFile("gaga.vtp"); + + error(refinement_level) = surfaceL2error(ansatz_space, x, fun); + // std::cout << error(refinement_level) << std::endl; + } + + // estimate rate of convergence and check whether it is at least 90% of the + // expected value + // assert(checkRateOfConvergence(error.tail(2), polynomial_degree + 1, + // 0.9)); + + std::cout << std::endl; + } + // The VTKwriter sets up initial geomety information. + std::cout << "\n" << std::string(60, '=') << "\n"; + + return 0; +} From 138fe0d815e520b7c2b0797bdc77c305cebf0046 Mon Sep 17 00:00:00 2001 From: jdoelz <48793870+jdoelz@users.noreply.github.com> Date: Tue, 19 Mar 2024 08:57:32 +0100 Subject: [PATCH 14/20] fixed bug --- Bembel/src/Geometry/Patch.hpp | 1 + Bembel/src/Geometry/SurfacePoint.hpp | 7 +++- Bembel/src/LinearForm/DirichletTrace.hpp | 1 - .../DirichletTraceOnReferenceDomain.hpp | 3 +- examples/CMakeLists.txt | 2 +- ...pp => DirichletTraceOnReferenceDomain.cpp} | 38 ++++--------------- 6 files changed, 17 insertions(+), 35 deletions(-) rename examples/{DiscreteLinearFormOnReferenceDomain.cpp => DirichletTraceOnReferenceDomain.cpp} (70%) diff --git a/Bembel/src/Geometry/Patch.hpp b/Bembel/src/Geometry/Patch.hpp index 76f61d7c0..e163a7aec 100644 --- a/Bembel/src/Geometry/Patch.hpp +++ b/Bembel/src/Geometry/Patch.hpp @@ -336,6 +336,7 @@ class Patch { srf_pt->set_patch(patch); srf_pt->set_xi(xi); + srf_pt->set_xi_patch(ref_pt); srf_pt->set_w(w); srf_pt->set_f(bot * tmpvec); srf_pt->set_jacobian(botsqr * (Eigen::Matrix() diff --git a/Bembel/src/Geometry/SurfacePoint.hpp b/Bembel/src/Geometry/SurfacePoint.hpp index 4f061dd8b..7e124b45c 100644 --- a/Bembel/src/Geometry/SurfacePoint.hpp +++ b/Bembel/src/Geometry/SurfacePoint.hpp @@ -34,6 +34,10 @@ class SurfacePoint { Eigen::Vector2d get_xi() { return xi_; } const Eigen::Vector2d get_xi() const { return xi_; } + // quadrature point on patch + Eigen::Vector2d get_xi_patch() { return xi_patch_; } + const Eigen::Vector2d get_xi_patch() const { return xi_patch_; } + // quadrature weight double get_w() { return w_; } const double get_w() const { return w_; } @@ -75,6 +79,7 @@ class SurfacePoint { // setter void set_patch(const int patch) { patch_ = patch; } void set_xi(const Eigen::Vector2d &xi) { xi_ = xi; } + void set_xi_patch(const Eigen::Vector2d &xi_patch) { xi_patch_ = xi_patch; } void set_w(const double w) { w_ = w; } void set_f(const Eigen::Vector3d &f) { f_ = f; } void set_jacobian(const Eigen::Matrix &jacobian) { @@ -98,6 +103,7 @@ class SurfacePoint { private: int patch_; Eigen::Vector2d xi_; + Eigen::Vector2d xi_patch_; double w_; Eigen::Vector3d f_; Eigen::Matrix jacobian_; @@ -118,4 +124,3 @@ typedef std::vector> ElementSurfacePoints; #endif // BEMBEL_SRC_GEOMETRY_SURFACEPOINT_HPP_ - diff --git a/Bembel/src/LinearForm/DirichletTrace.hpp b/Bembel/src/LinearForm/DirichletTrace.hpp index d8d33ccbe..c40e4deae 100644 --- a/Bembel/src/LinearForm/DirichletTrace.hpp +++ b/Bembel/src/LinearForm/DirichletTrace.hpp @@ -38,7 +38,6 @@ class DirichletTrace : public LinearFormBase, Scalar> { void evaluateIntegrand_impl( const T &super_space, const SurfacePoint &p, Eigen::Matrix *intval) const { - std::cout << p.get_xi() << std::endl; // integrand without basis functions Scalar integrand = function_(p.get_f()) * p.get_surface_measure() * p.get_w(); // multiply basis functions with integrand diff --git a/Bembel/src/LinearForm/DirichletTraceOnReferenceDomain.hpp b/Bembel/src/LinearForm/DirichletTraceOnReferenceDomain.hpp index f612e0123..4a588b604 100644 --- a/Bembel/src/LinearForm/DirichletTraceOnReferenceDomain.hpp +++ b/Bembel/src/LinearForm/DirichletTraceOnReferenceDomain.hpp @@ -41,8 +41,7 @@ class DirichletTraceOnReferenceDomain const T &super_space, const SurfacePoint &p, Eigen::Matrix *intval) const { // integrand without basis functions - //std::cout << p.get_xi() << std::endl; - Scalar integrand = function_(p.get_patch(), p.get_xi()) * + Scalar integrand = function_(p.get_patch(), p.get_xi_patch()) * p.get_surface_measure() * p.get_w(); // multiply basis functions with integrand super_space.addScaledBasis(intval, integrand, p.get_xi()); diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 4f0e8ebb2..204fe2e33 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -11,7 +11,7 @@ configure_file(${CMAKE_SOURCE_DIR}/geo/cube_small.dat ${CMAKE_CURRENT_BINARY_DIR set(CIFILES AnsatzSpace BlockClusterTree - DiscreteLinearFormOnReferenceDomain + DirichletTraceOnReferenceDomain Geometry MassMatrix LaplaceBeltrami diff --git a/examples/DiscreteLinearFormOnReferenceDomain.cpp b/examples/DirichletTraceOnReferenceDomain.cpp similarity index 70% rename from examples/DiscreteLinearFormOnReferenceDomain.cpp rename to examples/DirichletTraceOnReferenceDomain.cpp index 4b1d447bf..395de90bb 100644 --- a/examples/DiscreteLinearFormOnReferenceDomain.cpp +++ b/examples/DirichletTraceOnReferenceDomain.cpp @@ -23,8 +23,8 @@ int main() { using namespace Bembel; using namespace Eigen; - int polynomial_degree_max = 0; - int refinement_level_max = 1; + int polynomial_degree_max = 3; + int refinement_level_max = 3; std::function fun = [](const Vector3d &in) { return in(0); @@ -39,7 +39,7 @@ int main() { for (int refinement_level = 0; refinement_level < refinement_level_max + 1; ++refinement_level) { std::cout << "Degree " << polynomial_degree << " Level " - << refinement_level << "\t\t"; + << refinement_level << std::endl; // Build ansatz space AnsatzSpace ansatz_space(geometry, refinement_level, polynomial_degree); @@ -73,35 +73,13 @@ int main() { disc_lf2.get_linear_form().set_function(fun_disc); disc_lf2.compute(); - VectorXd x2 = solver.solve(disc_lf2.get_discrete_linear_form()); - - // std::cout << std::endl; - // std::cout << x.head(10).transpose() << std::endl; - // std::cout << x2.head(10).transpose() << std::endl; - // std::cout << (x.head(10).array()/x2.head(10).array()).transpose() << - // std::endl; - std::cout << (disc_lf.get_discrete_linear_form() - - disc_lf2.get_discrete_linear_form()) - .norm() / - disc_lf.get_discrete_linear_form().norm() - << std::endl; - - VTKSurfaceExport writer(geometry, 5); - writer.addDataSet("x", ansatz_space, x); - writer.addDataSet("x2", ansatz_space, x2); - writer.addDataSet("fun_disc", fun_disc); - writer.addDataSet("fun", fun); - writer.writeToFile("gaga.vtp"); - - error(refinement_level) = surfaceL2error(ansatz_space, x, fun); - // std::cout << error(refinement_level) << std::endl; + // check correctness + assert(std::abs((disc_lf.get_discrete_linear_form() - + disc_lf2.get_discrete_linear_form()) + .norm() / + disc_lf.get_discrete_linear_form().norm()) < 1e-10); } - // estimate rate of convergence and check whether it is at least 90% of the - // expected value - // assert(checkRateOfConvergence(error.tail(2), polynomial_degree + 1, - // 0.9)); - std::cout << std::endl; } // The VTKwriter sets up initial geomety information. From 62e8959860e007e976d4506b7808381c71dfa103 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20D=C3=B6lz?= Date: Tue, 19 Mar 2024 17:41:57 +0100 Subject: [PATCH 15/20] remove most of the compile errors --- Bembel/src/Geometry/SurfacePoint.hpp | 9 ++- .../Helmholtz/AdjointDoubleLayerOperator.hpp | 58 ++------------ Bembel/src/Helmholtz/DoubleLayerOperator.hpp | 58 ++------------ Bembel/src/Helmholtz/DoubleLayerPotential.hpp | 17 +--- .../src/Helmholtz/HypersingularOperator.hpp | 62 ++++----------- .../SingleLayerOperator.hpp | 79 +++++-------------- .../SingleLayerPotential.hpp | 60 ++++++-------- Bembel/src/Laplace/DoubleLayerOperator.hpp | 58 ++------------ Bembel/src/Laplace/DoubleLayerPotential.hpp | 18 +---- .../Laplace/SingleLayerPotentialGradient.hpp | 19 +---- Bembel/src/LinearForm/NeumannTrace.hpp | 22 +----- Bembel/src/LinearForm/TangentialTrace.hpp | 26 ++---- tests/test_Bernstein.cpp | 29 ++++--- 13 files changed, 120 insertions(+), 395 deletions(-) diff --git a/Bembel/src/Geometry/SurfacePoint.hpp b/Bembel/src/Geometry/SurfacePoint.hpp index 57f5126df..7e124b45c 100644 --- a/Bembel/src/Geometry/SurfacePoint.hpp +++ b/Bembel/src/Geometry/SurfacePoint.hpp @@ -11,7 +11,9 @@ // #ifndef BEMBEL_SRC_GEOMETRY_SURFACEPOINT_HPP_ #define BEMBEL_SRC_GEOMETRY_SURFACEPOINT_HPP_ + #include + /** * \ingroup Geometry * \brief This class stores quadrature and geometry information on quadrature @@ -32,6 +34,10 @@ class SurfacePoint { Eigen::Vector2d get_xi() { return xi_; } const Eigen::Vector2d get_xi() const { return xi_; } + // quadrature point on patch + Eigen::Vector2d get_xi_patch() { return xi_patch_; } + const Eigen::Vector2d get_xi_patch() const { return xi_patch_; } + // quadrature weight double get_w() { return w_; } const double get_w() const { return w_; } @@ -73,6 +79,7 @@ class SurfacePoint { // setter void set_patch(const int patch) { patch_ = patch; } void set_xi(const Eigen::Vector2d &xi) { xi_ = xi; } + void set_xi_patch(const Eigen::Vector2d &xi_patch) { xi_patch_ = xi_patch; } void set_w(const double w) { w_ = w; } void set_f(const Eigen::Vector3d &f) { f_ = f; } void set_jacobian(const Eigen::Matrix &jacobian) { @@ -96,6 +103,7 @@ class SurfacePoint { private: int patch_; Eigen::Vector2d xi_; + Eigen::Vector2d xi_patch_; double w_; Eigen::Vector3d f_; Eigen::Matrix jacobian_; @@ -116,4 +124,3 @@ typedef std::vector> ElementSurfacePoints; #endif // BEMBEL_SRC_GEOMETRY_SURFACEPOINT_HPP_ - diff --git a/Bembel/src/Helmholtz/AdjointDoubleLayerOperator.hpp b/Bembel/src/Helmholtz/AdjointDoubleLayerOperator.hpp index c5d9d13ef..ee39bc2e8 100644 --- a/Bembel/src/Helmholtz/AdjointDoubleLayerOperator.hpp +++ b/Bembel/src/Helmholtz/AdjointDoubleLayerOperator.hpp @@ -43,66 +43,24 @@ class HelmholtzAdjointDoubleLayerOperator Eigen::Matrix::Scalar, Eigen::Dynamic, Eigen::Dynamic> *intval) const { - auto polynomial_degree = super_space.get_polynomial_degree(); - auto polynomial_degree_plus_one_squared = - (polynomial_degree + 1) * (polynomial_degree + 1); - - // get evaluation points on unit square - auto s = p1.segment<2>(0); - auto t = p2.segment<2>(0); - - // get quadrature weights - auto ws = p1(2); - auto wt = p2(2); - - // get points on geometry and tangential derivatives - auto x_f = p1.segment<3>(3); - auto x_f_dx = p1.segment<3>(6); - auto x_f_dy = p1.segment<3>(9); - auto y_f = p2.segment<3>(3); - auto y_f_dx = p2.segment<3>(6); - auto y_f_dy = p2.segment<3>(9); - - // compute surface measures from tangential derivatives - auto x_n = x_f_dx.cross(x_f_dy); - auto y_kappa = y_f_dx.cross(y_f_dy).norm(); - // integrand without basis functions - auto integrand = evaluateKernelGrad(x_f, x_n, y_f) * y_kappa * ws * wt; + auto integrand = + evaluateKernelGrad(p1.get_f(), p1.get_normal(), p2.get_f()) * + p2.get_surface_measure() * p1.get_w() * p2.get_w(); - // multiply basis functions with integrand and add to intval, this is an - // efficient implementation of - // (*intval) += super_space.BasisInteraction(s, t) * evaluateKernel(x_f, - // y_f) - // * x_kappa * y_kappa * ws * wt; - super_space.addScaledBasisInteraction(intval, integrand, s, t); + // multiply basis functions with integrand and add to intval + super_space.addScaledBasisInteraction(intval, integrand, p1.get_xi(), + p2.get_xi()); return; } Eigen::Matrix, 1, 1> evaluateFMMInterpolation_impl( const SurfacePoint &p1, const SurfacePoint &p2) const { - // get evaluation points on unit square - auto s = p1.segment<2>(0); - auto t = p2.segment<2>(0); - - // get points on geometry and tangential derivatives - auto x_f = p1.segment<3>(3); - auto x_f_dx = p1.segment<3>(6); - auto x_f_dy = p1.segment<3>(9); - auto y_f = p2.segment<3>(3); - auto y_f_dx = p2.segment<3>(6); - auto y_f_dy = p2.segment<3>(9); - - // compute surface measure in x and unnormalized normal in y from tangential - // derivatives - auto x_n = x_f_dx.cross(x_f_dy); - auto y_kappa = y_f_dx.cross(y_f_dy).norm(); - // interpolation Eigen::Matrix, 1, 1> intval; - intval(0) = evaluateKernelGrad(x_f, x_n, y_f) * y_kappa; - + intval(0) = evaluateKernelGrad(p1.get_f(), p1.get_normal(), p2.get_f()) * + p2.get_surface_measure(); return intval; } diff --git a/Bembel/src/Helmholtz/DoubleLayerOperator.hpp b/Bembel/src/Helmholtz/DoubleLayerOperator.hpp index 11407cbda..1abb7da8c 100644 --- a/Bembel/src/Helmholtz/DoubleLayerOperator.hpp +++ b/Bembel/src/Helmholtz/DoubleLayerOperator.hpp @@ -43,66 +43,24 @@ class HelmholtzDoubleLayerOperator Eigen::Matrix< typename LinearOperatorTraits::Scalar, Eigen::Dynamic, Eigen::Dynamic> *intval) const { - auto polynomial_degree = super_space.get_polynomial_degree(); - auto polynomial_degree_plus_one_squared = - (polynomial_degree + 1) * (polynomial_degree + 1); - - // get evaluation points on unit square - auto s = p1.segment<2>(0); - auto t = p2.segment<2>(0); - - // get quadrature weights - auto ws = p1(2); - auto wt = p2(2); - - // get points on geometry and tangential derivatives - auto x_f = p1.segment<3>(3); - auto x_f_dx = p1.segment<3>(6); - auto x_f_dy = p1.segment<3>(9); - auto y_f = p2.segment<3>(3); - auto y_f_dx = p2.segment<3>(6); - auto y_f_dy = p2.segment<3>(9); - - // compute surface measures from tangential derivatives - auto x_kappa = x_f_dx.cross(x_f_dy).norm(); - auto y_n = y_f_dx.cross(y_f_dy); - // integrand without basis functions - auto integrand = evaluateKernelGrad(x_f, y_f, y_n) * x_kappa * ws * wt; + auto integrand = + evaluateKernelGrad(p1.get_f(), p1.get_normal(), p2.get_f()) * + p1.get_surface_measure() * p1.get_w() * p2.get_w(); - // multiply basis functions with integrand and add to intval, this is an - // efficient implementation of - // (*intval) += super_space.BasisInteraction(s, t) * evaluateKernel(x_f, - // y_f) - // * x_kappa * y_kappa * ws * wt; - super_space.addScaledBasisInteraction(intval, integrand, s, t); + // multiply basis functions with integrand and add to intval + super_space.addScaledBasisInteraction(intval, integrand, p1.get_xi(), + p2.get_xi()); return; } Eigen::Matrix, 1, 1> evaluateFMMInterpolation_impl( const SurfacePoint &p1, const SurfacePoint &p2) const { - // get evaluation points on unit square - auto s = p1.segment<2>(0); - auto t = p2.segment<2>(0); - - // get points on geometry and tangential derivatives - auto x_f = p1.segment<3>(3); - auto x_f_dx = p1.segment<3>(6); - auto x_f_dy = p1.segment<3>(9); - auto y_f = p2.segment<3>(3); - auto y_f_dx = p2.segment<3>(6); - auto y_f_dy = p2.segment<3>(9); - - // compute surface measure in x and unnormalized normal in y from tangential - // derivatives - auto x_kappa = x_f_dx.cross(x_f_dy).norm(); - auto y_n = y_f_dx.cross(y_f_dy); - // interpolation Eigen::Matrix, 1, 1> intval; - intval(0) = evaluateKernelGrad(x_f, y_f, y_n) * x_kappa; - + intval(0) = evaluateKernelGrad(p1.get_f(), p1.get_normal(), p2.get_f()) * + p2.get_surface_measure(); return intval; } diff --git a/Bembel/src/Helmholtz/DoubleLayerPotential.hpp b/Bembel/src/Helmholtz/DoubleLayerPotential.hpp index 4d6ebccc6..44a14bcab 100644 --- a/Bembel/src/Helmholtz/DoubleLayerPotential.hpp +++ b/Bembel/src/Helmholtz/DoubleLayerPotential.hpp @@ -41,26 +41,13 @@ class HelmholtzDoubleLayerPotential const ElementTreeNode &element, const Eigen::Vector3d &point, const SurfacePoint &p) const { - // get evaluation points on unit square - auto s = p.segment<2>(0); - - // get quadrature weights - auto ws = p(2); - - // get points on geometry and tangential derivatives - auto x_f = p.segment<3>(3); - auto x_f_dx = p.segment<3>(6); - auto x_f_dy = p.segment<3>(9); - - // compute unnormalized normal from tangential derivatives - auto x_n = x_f_dx.cross(x_f_dy); - // assemble Galerkin solution auto cauchy_value = fun_ev.evaluate(element, p); // integrand without basis functions // dot: adjoint in first variable - auto integrand = evaluateKernelGrad(point, x_f, x_n) * cauchy_value * ws; + auto integrand = evaluateKernelGrad(point, p.get_f(), p.get_normal()) * + cauchy_value * p.get_w(); return integrand; } diff --git a/Bembel/src/Helmholtz/HypersingularOperator.hpp b/Bembel/src/Helmholtz/HypersingularOperator.hpp index a63d6392a..e072ce59b 100644 --- a/Bembel/src/Helmholtz/HypersingularOperator.hpp +++ b/Bembel/src/Helmholtz/HypersingularOperator.hpp @@ -41,46 +41,27 @@ class HelmholtzHypersingularOperator void evaluateIntegrand_impl(const T &super_space, const SurfacePoint &p1, const SurfacePoint &p2, Eigen::MatrixXcd *intval) const { - auto polynomial_degree = super_space.get_polynomial_degree(); - auto polynomial_degree_plus_one_squared = - (polynomial_degree + 1) * (polynomial_degree + 1); - - // get evaluation points on unit square - Eigen::Vector2d s = p1.segment<2>(0); - Eigen::Vector2d t = p2.segment<2>(0); - - // get quadrature weights - double ws = p1(2); - double wt = p2(2); - - // get points on geometry and tangential derivatives - Eigen::Vector3d x_f = p1.segment<3>(3); - Eigen::Vector3d x_f_dx = p1.segment<3>(6); - Eigen::Vector3d x_f_dy = p1.segment<3>(9); - Eigen::Vector3d y_f = p2.segment<3>(3); - Eigen::Vector3d y_f_dx = p2.segment<3>(6); - Eigen::Vector3d y_f_dy = p2.segment<3>(9); - // compute surface measures from tangential derivatives - Eigen::Vector3d x_n = x_f_dx.cross(x_f_dy); - Eigen::Vector3d y_n = y_f_dx.cross(y_f_dy); + Eigen::Vector3d x_n = p1.get_normal(); + Eigen::Vector3d y_n = p2.get_normal(); // compute h double h = 1. / (1 << super_space.get_refinement_level()); // h = 1 ./ (2^M) // evaluate kernel - std::complex kernel = evaluateKernel(x_f, y_f); + std::complex kernel = evaluateKernel(p1.get_f(), p2.get_f()); // integrand without basis functions std::complex integrandScalar = - -kernel * x_n.dot(y_n) * wavenumber2_ * ws * wt; + -kernel * x_n.dot(y_n) * wavenumber2_ * p1.get_w() * p2.get_w(); std::complex integrandCurl = - kernel * x_n.norm() * y_n.norm() * ws * wt / h / h; + kernel * x_n.norm() * y_n.norm() * p1.get_w() * p2.get_w() / h / h; // multiply basis functions with integrand and add to intval, this is an // efficient implementation of - super_space.addScaledBasisInteraction(intval, integrandScalar, s, t); + super_space.addScaledBasisInteraction(intval, integrandScalar, p1.get_xi(), + p2.get_xi()); super_space.addScaledSurfaceCurlInteraction(intval, integrandCurl, p1, p2); return; @@ -88,33 +69,18 @@ class HelmholtzHypersingularOperator Eigen::Matrix, 3, 3> evaluateFMMInterpolation_impl( const SurfacePoint &p1, const SurfacePoint &p2) const { - // get evaluation points on unit square - Eigen::Vector2d s = p1.segment<2>(0); - Eigen::Vector2d t = p2.segment<2>(0); - - // get points on geometry and tangential derivatives - Eigen::Vector3d x_f = p1.segment<3>(3); - Eigen::Vector3d x_f_dx = p1.segment<3>(6); - Eigen::Vector3d x_f_dy = p1.segment<3>(9); - Eigen::Vector3d y_f = p2.segment<3>(3); - Eigen::Vector3d y_f_dx = p2.segment<3>(6); - Eigen::Vector3d y_f_dy = p2.segment<3>(9); - - // compute surface measures from tangential derivatives - Eigen::Vector3d x_n = x_f_dx.cross(x_f_dy); - Eigen::Vector3d y_n = y_f_dx.cross(y_f_dy); - // evaluate kernel - std::complex kernel = evaluateKernel(x_f, y_f); + std::complex kernel = evaluateKernel(p1.get_f(), p2.get_f()); // interpolation Eigen::Matrix, 3, 3> intval; intval.setZero(); - intval(0, 0) = -kernel * wavenumber2_ * x_n.dot(y_n); - intval(1, 1) = kernel * x_f_dy.dot(y_f_dy); - intval(1, 2) = -kernel * x_f_dy.dot(y_f_dx); - intval(2, 1) = -kernel * x_f_dx.dot(y_f_dy); - intval(2, 2) = kernel * x_f_dx.dot(y_f_dx); + intval(0, 0) = + -kernel * wavenumber2_ * p1.get_normal().dot(p2.get_normal()); + intval(1, 1) = kernel * p1.get_f_dy().dot(p2.get_f_dy()); + intval(1, 2) = -kernel * p1.get_f_dy().dot(p2.get_f_dx()); + intval(2, 1) = -kernel * p1.get_f_dx().dot(p2.get_f_dy()); + intval(2, 2) = kernel * p1.get_f_dx().dot(p2.get_f_dx()); return intval; } diff --git a/Bembel/src/HomogenisedLaplace/SingleLayerOperator.hpp b/Bembel/src/HomogenisedLaplace/SingleLayerOperator.hpp index 0f528b4a8..760b53068 100644 --- a/Bembel/src/HomogenisedLaplace/SingleLayerOperator.hpp +++ b/Bembel/src/HomogenisedLaplace/SingleLayerOperator.hpp @@ -20,7 +20,7 @@ class HomogenisedLaplaceSingleLayerOperator; /** * \brief Specification of the LinerOperatorTraits for the Homogenised Laplace. */ -template<> +template <> struct LinearOperatorTraits { typedef Eigen::VectorXd EigenType; typedef Eigen::VectorXd::Scalar Scalar; @@ -36,8 +36,8 @@ struct LinearOperatorTraits { * \brief This class implements the specification of the integration for the * single layer operator for the homogenised Laplace. */ -class HomogenisedLaplaceSingleLayerOperator : public LinearOperatorBase< - HomogenisedLaplaceSingleLayerOperator> { +class HomogenisedLaplaceSingleLayerOperator + : public LinearOperatorBase { // implementation of the kernel evaluation, which may be based on the // information available from the superSpace public: @@ -47,74 +47,37 @@ class HomogenisedLaplaceSingleLayerOperator : public LinearOperatorBase< */ HomogenisedLaplaceSingleLayerOperator() { this->deg = getDegree(HomogenisedLaplaceSingleLayerOperator::precision); - this->cs = getCoefficients( - HomogenisedLaplaceSingleLayerOperator::precision); + this->cs = + getCoefficients(HomogenisedLaplaceSingleLayerOperator::precision); } - template - void evaluateIntegrand_impl(const T &super_space, const SurfacePoint &p1, - const SurfacePoint &p2, - Eigen::Matrix< - typename LinearOperatorTraits::Scalar, Eigen::Dynamic, Eigen::Dynamic> *intval) const { - auto polynomial_degree = super_space.get_polynomial_degree(); - auto polynomial_degree_plus_one_squared = (polynomial_degree + 1) - * (polynomial_degree + 1); - - // get evaluation points on unit square - auto s = p1.segment < 2 > (0); - auto t = p2.segment < 2 > (0); - - // get quadrature weights - auto ws = p1(2); - auto wt = p2(2); - - // get points on geometry and tangential derivatives - auto x_f = p1.segment < 3 > (3); - auto x_f_dx = p1.segment < 3 > (6); - auto x_f_dy = p1.segment < 3 > (9); - auto y_f = p2.segment < 3 > (3); - auto y_f_dx = p2.segment < 3 > (6); - auto y_f_dy = p2.segment < 3 > (9); - - // compute surface measures from tangential derivatives - auto x_kappa = x_f_dx.cross(x_f_dy).norm(); - auto y_kappa = y_f_dx.cross(y_f_dy).norm(); - + template + void evaluateIntegrand_impl( + const T &super_space, const SurfacePoint &p1, const SurfacePoint &p2, + Eigen::Matrix::Scalar, + Eigen::Dynamic, Eigen::Dynamic> *intval) const { // integrand without basis functions - auto integrand = evaluateKernel(x_f, y_f) * x_kappa * y_kappa * ws * wt; + auto integrand = evaluateKernel(p1.get_f(), p2.get_f()) * + p1.get_surface_measure() * p2.get_surface_measure() * + p1.get_w() * p2.get_w(); // multiply basis functions with integrand and add to intval, this is an // efficient implementation of // (*intval) += super_space.basisInteraction(s, t) // * evaluateKernel(x_f, y_f) * x_kappa * y_kappa * ws * wt; - super_space.addScaledBasisInteraction(intval, integrand, s, t); + super_space.addScaledBasisInteraction(intval, integrand, p1.get_f(), + p2.get_f()); return; } Eigen::Matrix evaluateFMMInterpolation_impl( const SurfacePoint &p1, const SurfacePoint &p2) const { - // get evaluation points on unit square - auto s = p1.segment < 2 > (0); - auto t = p2.segment < 2 > (0); - - // get points on geometry and tangential derivatives - auto x_f = p1.segment < 3 > (3); - auto x_f_dx = p1.segment < 3 > (6); - auto x_f_dy = p1.segment < 3 > (9); - auto y_f = p2.segment < 3 > (3); - auto y_f_dx = p2.segment < 3 > (6); - auto y_f_dy = p2.segment < 3 > (9); - - // compute surface measures from tangential derivatives - auto x_kappa = x_f_dx.cross(x_f_dy).norm(); - auto y_kappa = y_f_dx.cross(y_f_dy).norm(); - // interpolation Eigen::Matrix intval; - intval(0) = evaluateKernel(x_f, y_f) * x_kappa * y_kappa; - + intval(0) = evaluateKernel(p1.get_f(), p2.get_f()) * + p1.get_surface_measure() * p2.get_surface_measure(); return intval; } @@ -122,9 +85,9 @@ class HomogenisedLaplaceSingleLayerOperator : public LinearOperatorBase< * \brief Fundamental solution of the Homogenised Laplace problem */ double evaluateKernel(const Eigen::Vector3d &x, - const Eigen::Vector3d &y) const { - return k_mod(x - y) - + evaluate_solid_sphericals(x - y, this->cs, this->deg, false); + const Eigen::Vector3d &y) const { + return k_mod(x - y) + + evaluate_solid_sphericals(x - y, this->cs, this->deg, false); } /** diff --git a/Bembel/src/HomogenisedLaplace/SingleLayerPotential.hpp b/Bembel/src/HomogenisedLaplace/SingleLayerPotential.hpp index 983276ff2..b8c70d7bf 100644 --- a/Bembel/src/HomogenisedLaplace/SingleLayerPotential.hpp +++ b/Bembel/src/HomogenisedLaplace/SingleLayerPotential.hpp @@ -15,13 +15,13 @@ namespace Bembel { // forward declaration of class HomogenisedLaplaceSingleLayerPotential // in order to define traits -template +template class HomogenisedLaplaceSingleLayerPotential; /** * \brief Specification of the PotentialTraits for the Homogenised Laplace. */ -template +template struct PotentialTraits> { typedef Eigen::VectorXd::Scalar Scalar; static constexpr int OutputSpaceDimension = 1; @@ -32,51 +32,41 @@ struct PotentialTraits> { * \brief This class implements the specification of the integration for the * single layer potential for the homogenised Laplace. */ -template -class HomogenisedLaplaceSingleLayerPotential : public PotentialBase< - HomogenisedLaplaceSingleLayerPotential, LinOp> { +template +class HomogenisedLaplaceSingleLayerPotential + : public PotentialBase, + LinOp> { // implementation of the kernel evaluation, which may be based on the // information available from the superSpace public: /** - * \brief Constructs an object initialising the coefficients and the degree - * via the static variable HomogenisedLaplaceSingleLayerOperator::precision. - */ + * \brief Constructs an object initialising the coefficients and the degree + * via the static variable HomogenisedLaplaceSingleLayerOperator::precision. + */ HomogenisedLaplaceSingleLayerPotential() { - this->deg = getDegree( - HomogenisedLaplaceSingleLayerOperator::getPrecision()); - this->cs = getCoefficients( - HomogenisedLaplaceSingleLayerOperator::getPrecision()); + this->deg = + getDegree(HomogenisedLaplaceSingleLayerOperator::getPrecision()); + this->cs = + getCoefficients(HomogenisedLaplaceSingleLayerOperator::getPrecision()); } Eigen::Matrix< typename PotentialReturnScalar< - typename LinearOperatorTraits::Scalar, - double>::Scalar, 1, 1> evaluateIntegrand_impl( - const FunctionEvaluator &fun_ev, const ElementTreeNode &element, - const Eigen::Vector3d &point, const SurfacePoint &p) const { - // get evaluation points on unit square - auto s = p.segment < 2 > (0); - - // get quadrature weights - auto ws = p(2); - - // get points on geometry and tangential derivatives - auto x_f = p.segment < 3 > (3); - auto x_f_dx = p.segment < 3 > (6); - auto x_f_dy = p.segment < 3 > (9); - - // compute surface measures from tangential derivatives - auto x_kappa = x_f_dx.cross(x_f_dy).norm(); - + typename LinearOperatorTraits::Scalar, double>::Scalar, + 1, 1> + evaluateIntegrand_impl(const FunctionEvaluator &fun_ev, + const ElementTreeNode &element, + const Eigen::Vector3d &point, + const SurfacePoint &p) const { // evaluate kernel - auto kernel = evaluateKernel(point, x_f); + auto kernel = evaluateKernel(point, p.get_f()); // assemble Galerkin solution auto cauchy_value = fun_ev.evaluate(element, p); // integrand without basis functions - auto integrand = kernel * cauchy_value * x_kappa * ws; + auto integrand = + kernel * cauchy_value * p.get_surface_measure() * p.get_w(); return integrand; } @@ -85,9 +75,9 @@ class HomogenisedLaplaceSingleLayerPotential : public PotentialBase< * \brief Fundamental solution of the homogenised Laplace problem */ double evaluateKernel(const Eigen::Vector3d &x, - const Eigen::Vector3d &y) const { - return k_mod(x - y) - + evaluate_solid_sphericals(x - y, this->cs, this->deg, false); + const Eigen::Vector3d &y) const { + return k_mod(x - y) + + evaluate_solid_sphericals(x - y, this->cs, this->deg, false); } private: diff --git a/Bembel/src/Laplace/DoubleLayerOperator.hpp b/Bembel/src/Laplace/DoubleLayerOperator.hpp index 95260fb5f..ff1398ce8 100644 --- a/Bembel/src/Laplace/DoubleLayerOperator.hpp +++ b/Bembel/src/Laplace/DoubleLayerOperator.hpp @@ -43,68 +43,24 @@ class LaplaceDoubleLayerOperator Eigen::Matrix< typename LinearOperatorTraits::Scalar, Eigen::Dynamic, Eigen::Dynamic> *intval) const { - auto polynomial_degree = super_space.get_polynomial_degree(); - auto polynomial_degree_plus_one_squared = - (polynomial_degree + 1) * (polynomial_degree + 1); - - // get evaluation points on unit square - auto s = p1.segment<2>(0); - auto t = p2.segment<2>(0); - - // get quadrature weights - auto ws = p1(2); - auto wt = p2(2); - - // get points on geometry and tangential derivatives - auto x_f = p1.segment<3>(3); - auto x_f_dx = p1.segment<3>(6); - auto x_f_dy = p1.segment<3>(9); - auto y_f = p2.segment<3>(3); - auto y_f_dx = p2.segment<3>(6); - auto y_f_dy = p2.segment<3>(9); - - // compute surface measures from tangential derivatives - auto x_kappa = x_f_dx.cross(x_f_dy).norm(); - auto y_kappa = y_f_dx.cross(y_f_dy).norm(); - auto y_n = y_f_dx.cross(y_f_dy).normalized(); - // integrand without basis functions auto integrand = - evaluateKernelGrad(x_f, y_f, y_n) * y_kappa * x_kappa * ws * wt; + evaluateKernelGrad(p1.get_f(), p2.get_f(), p2.get_normal()) * + p1.get_surface_measure() * p1.get_w() * p2.get_w(); - // multiply basis functions with integrand and add to intval, this is an - // efficient implementation of - // (*intval) += super_space.BasisInteraction(s, t) * evaluateKernel(x_f, - // y_f) - // * x_kappa * y_kappa * ws * wt; - super_space.addScaledBasisInteraction(intval, integrand, s, t); + // multiply basis functions with integrand and add to intval + super_space.addScaledBasisInteraction(intval, integrand, p1.get_f(), + p2.get_f()); return; } Eigen::Matrix evaluateFMMInterpolation_impl( const SurfacePoint &p1, const SurfacePoint &p2) const { - // get evaluation points on unit square - auto s = p1.segment<2>(0); - auto t = p2.segment<2>(0); - - // get points on geometry and tangential derivatives - auto x_f = p1.segment<3>(3); - auto x_f_dx = p1.segment<3>(6); - auto x_f_dy = p1.segment<3>(9); - auto y_f = p2.segment<3>(3); - auto y_f_dx = p2.segment<3>(6); - auto y_f_dy = p2.segment<3>(9); - - // compute surface measure in x and unnormalized normal in y from tangential - // derivatives - auto x_kappa = x_f_dx.cross(x_f_dy).norm(); - auto y_n = y_f_dx.cross(y_f_dy); - // interpolation Eigen::Matrix intval; - intval(0) = evaluateKernelGrad(x_f, y_f, y_n) * x_kappa; - + intval(0) = evaluateKernelGrad(p1.get_f(), p2.get_f(), p2.get_normal()) * + p1.get_surface_measure(); return intval; } diff --git a/Bembel/src/Laplace/DoubleLayerPotential.hpp b/Bembel/src/Laplace/DoubleLayerPotential.hpp index c4d08cda3..a5d210679 100644 --- a/Bembel/src/Laplace/DoubleLayerPotential.hpp +++ b/Bembel/src/Laplace/DoubleLayerPotential.hpp @@ -41,26 +41,12 @@ class LaplaceDoubleLayerPotential const ElementTreeNode &element, const Eigen::Vector3d &point, const SurfacePoint &p) const { - // get evaluation points on unit square - auto s = p.segment<2>(0); - - // get quadrature weights - auto ws = p(2); - - // get points on geometry and tangential derivatives - auto x_f = p.segment<3>(3); - auto x_f_dx = p.segment<3>(6); - auto x_f_dy = p.segment<3>(9); - - // compute unnormalized normal from tangential derivatives - auto x_n = x_f_dx.cross(x_f_dy); - // assemble Galerkin solution auto cauchy_value = fun_ev.evaluate(element, p); // integrand without basis functions - auto integrand = - evaluateKernelGrad(point, x_f).dot(x_n) * cauchy_value * ws; + auto integrand = evaluateKernelGrad(point, p.get_f()).dot(p.get_normal()) * + cauchy_value * p.get_w(); return integrand; } diff --git a/Bembel/src/Laplace/SingleLayerPotentialGradient.hpp b/Bembel/src/Laplace/SingleLayerPotentialGradient.hpp index c04bca505..dac186379 100644 --- a/Bembel/src/Laplace/SingleLayerPotentialGradient.hpp +++ b/Bembel/src/Laplace/SingleLayerPotentialGradient.hpp @@ -41,28 +41,15 @@ class LaplaceSingleLayerPotentialGradient const ElementTreeNode &element, const Eigen::Vector3d &point, const SurfacePoint &p) const { - // get evaluation points on unit square - auto s = p.segment<2>(0); - - // get quadrature weights - auto ws = p(2); - - // get points on geometry and tangential derivatives - auto x_f = p.segment<3>(3); - auto x_f_dx = p.segment<3>(6); - auto x_f_dy = p.segment<3>(9); - - // compute surface measures from tangential derivatives - auto x_kappa = x_f_dx.cross(x_f_dy).norm(); - // evaluate kernel - auto kernel = evaluateKernelGrad(point, x_f); + auto kernel = evaluateKernelGrad(point, p.get_f()); // assemble Galerkin solution auto cauchy_value = fun_ev.evaluate(element, p); // integrand without basis functions - auto integrand = kernel * cauchy_value * x_kappa * ws; + auto integrand = + kernel * cauchy_value * p.get_surface_measure() * p.get_w(); return integrand; } diff --git a/Bembel/src/LinearForm/NeumannTrace.hpp b/Bembel/src/LinearForm/NeumannTrace.hpp index 1a14a2430..b4936eaf6 100644 --- a/Bembel/src/LinearForm/NeumannTrace.hpp +++ b/Bembel/src/LinearForm/NeumannTrace.hpp @@ -40,30 +40,12 @@ class NeumannTrace : public LinearFormBase, Scalar> { void evaluateIntegrand_impl( const T &super_space, const SurfacePoint &p, Eigen::Matrix *intval) const { - auto polynomial_degree = super_space.get_polynomial_degree(); - auto polynomial_degree_plus_one_squared = - (polynomial_degree + 1) * (polynomial_degree + 1); - - // get evaluation points on unit square - auto s = p.segment<2>(0); - - // get quadrature weights - auto ws = p(2); - - // get points on geometry and tangential derivatives - auto x_f = p.segment<3>(3); - auto x_f_dx = p.segment<3>(6); - auto x_f_dy = p.segment<3>(9); - - // compute (unnormalized) surface normal from tangential derivatives - auto x_f_n = x_f_dx.cross(x_f_dy); - // integrand without basis functions // dot: adjoint in first variable - auto integrand = x_f_n.dot(function_(x_f)) * ws; + auto integrand = p.get_normal().dot(function_(p.get_f())) * p.get_w(); // multiply basis functions with integrand - super_space.addScaledBasis(intval, integrand, s); + super_space.addScaledBasis(intval, integrand, p.get_xi()); return; } diff --git a/Bembel/src/LinearForm/TangentialTrace.hpp b/Bembel/src/LinearForm/TangentialTrace.hpp index fc79f9e6f..eb1901106 100644 --- a/Bembel/src/LinearForm/TangentialTrace.hpp +++ b/Bembel/src/LinearForm/TangentialTrace.hpp @@ -42,32 +42,18 @@ class TangentialTrace : public LinearFormBase, Scalar> { int polynomial_degree = super_space.get_polynomial_degree(); int polynomial_degree_plus_one_squared = (polynomial_degree + 1) * (polynomial_degree + 1); - - // get evaluation points on unit square - Eigen::Vector2d s = p.segment<2>(0); - - // get quadrature weights - double ws = p(2); - - // get points on geometry and tangential derivatives - Eigen::Vector3d x_f = p.segment<3>(3); - Eigen::Vector3d x_f_dx = p.segment<3>(6); - Eigen::Vector3d x_f_dy = p.segment<3>(9); - - // compute surface measures from tangential derivatives - Eigen::Vector3d x_n = x_f_dx.cross(x_f_dy).normalized(); - // tangential component + quadrature weights - Eigen::Matrix fun_x_f = function_(x_f); - Eigen::Matrix tangential_component = fun_x_f.cross(x_n) * ws; + Eigen::Matrix fun_x_f = function_(p.get_f()); + Eigen::Matrix tangential_component = + fun_x_f.cross(p.get_unit_normal()) * p.get_w(); // extract tangential component - Scalar component_x = x_f_dx.dot(tangential_component); - Scalar component_y = x_f_dy.dot(tangential_component); + Scalar component_x = p.get_f_dx().dot(tangential_component); + Scalar component_y = p.get_f_dy().dot(tangential_component); // evaluate shape functions Eigen::Matrix phiPhiVec = - super_space.basis(s); + super_space.basis(p.get_xi()); // multiply basis functions with integrand Eigen::Matrix phiPhiMat( diff --git a/tests/test_Bernstein.cpp b/tests/test_Bernstein.cpp index 0e876d41d..d2ef27ba7 100644 --- a/tests/test_Bernstein.cpp +++ b/tests/test_Bernstein.cpp @@ -15,25 +15,23 @@ int main() { using namespace Bembel; + using namespace Eigen; // We test the Bernstein with given control points against the deBoor code constexpr int P = Bembel::Constants::MaxP; - double* coefs = new double[P + 1]; - - for (int i = 0; i < P + 1; ++i) { - coefs[i] = (i + 1) / static_cast(P + 1); - } - - Eigen::Map coefs_vector(coefs, 1, P + 1); + VectorXd buffer(P + 1); + VectorXd coefs = + VectorXd::LinSpaced(1, P + 1, P + 1) / static_cast(P + 1); for (int p = 0; p <= Bembel::Constants::MaxP; ++p) { for (auto x : Test::Constants::eq_points) { - double result1 = Basis::ShapeFunctionHandler::evalCoef(p, coefs, x); + buffer.setZero(); + Basis::ShapeFunctionHandler::evalBasis(p, buffer.data(), x); + double result1 = buffer.dot(coefs); std::vector v = {x}; - double result2 = - Spl::DeBoor(Eigen::MatrixXd(coefs_vector.leftCols(p + 1)), - Spl::MakeBezierKnotVector(p + 1), v)(0); + double result2 = Spl::DeBoor(Eigen::MatrixXd(coefs.leftCols(p + 1)), + Spl::MakeBezierKnotVector(p + 1), v)(0); BEMBEL_TEST_IF(std::abs(result1 - result2) < Test::Constants::coefficient_accuracy); @@ -43,12 +41,13 @@ int main() { // Now, we do the same for the derivatives for (int p = 1; p <= Bembel::Constants::MaxP; ++p) { for (auto x : Test::Constants::eq_points) { - double result1 = Basis::ShapeFunctionHandler::evalDerCoef(p, coefs, x); + buffer.setZero(); + Basis::ShapeFunctionHandler::evalDerCoef(p, buffer, x); + double result1 = buffer.dot(coefs); std::vector v = {x}; - double result2 = - Spl::DeBoorDer(Eigen::MatrixXd(coefs_vector.leftCols(p + 1)), - Spl::MakeBezierKnotVector(p + 1), v)(0); + double result2 = Spl::DeBoorDer(Eigen::MatrixXd(coefs.leftCols(p + 1)), + Spl::MakeBezierKnotVector(p + 1), v)(0); BEMBEL_TEST_IF(std::abs(result1 - result2) < Test::Constants::coefficient_accuracy); From 61127dfd1b3741a614ccb7772da7f9343cf621ce Mon Sep 17 00:00:00 2001 From: jdoelz <48793870+jdoelz@users.noreply.github.com> Date: Wed, 20 Mar 2024 09:34:41 +0100 Subject: [PATCH 16/20] remove minor bugs caused by merging --- Bembel/src/Helmholtz/DoubleLayerOperator.hpp | 6 +++--- Bembel/src/HomogenisedLaplace/SingleLayerOperator.hpp | 9 +++------ Bembel/src/util/surfaceL2error.hpp | 4 ++-- tests/test_Bernstein.cpp | 9 ++++----- 4 files changed, 12 insertions(+), 16 deletions(-) diff --git a/Bembel/src/Helmholtz/DoubleLayerOperator.hpp b/Bembel/src/Helmholtz/DoubleLayerOperator.hpp index 1abb7da8c..485093b3f 100644 --- a/Bembel/src/Helmholtz/DoubleLayerOperator.hpp +++ b/Bembel/src/Helmholtz/DoubleLayerOperator.hpp @@ -45,7 +45,7 @@ class HelmholtzDoubleLayerOperator Eigen::Dynamic, Eigen::Dynamic> *intval) const { // integrand without basis functions auto integrand = - evaluateKernelGrad(p1.get_f(), p1.get_normal(), p2.get_f()) * + evaluateKernelGrad(p1.get_f(), p2.get_f(), p2.get_normal()) * p1.get_surface_measure() * p1.get_w() * p2.get_w(); // multiply basis functions with integrand and add to intval @@ -59,8 +59,8 @@ class HelmholtzDoubleLayerOperator const SurfacePoint &p1, const SurfacePoint &p2) const { // interpolation Eigen::Matrix, 1, 1> intval; - intval(0) = evaluateKernelGrad(p1.get_f(), p1.get_normal(), p2.get_f()) * - p2.get_surface_measure(); + intval(0) = evaluateKernelGrad(p1.get_f(), p2.get_f(), p2.get_normal()) * + p1.get_surface_measure(); return intval; } diff --git a/Bembel/src/HomogenisedLaplace/SingleLayerOperator.hpp b/Bembel/src/HomogenisedLaplace/SingleLayerOperator.hpp index 760b53068..d171166c9 100644 --- a/Bembel/src/HomogenisedLaplace/SingleLayerOperator.hpp +++ b/Bembel/src/HomogenisedLaplace/SingleLayerOperator.hpp @@ -62,12 +62,9 @@ class HomogenisedLaplaceSingleLayerOperator p1.get_surface_measure() * p2.get_surface_measure() * p1.get_w() * p2.get_w(); - // multiply basis functions with integrand and add to intval, this is an - // efficient implementation of - // (*intval) += super_space.basisInteraction(s, t) - // * evaluateKernel(x_f, y_f) * x_kappa * y_kappa * ws * wt; - super_space.addScaledBasisInteraction(intval, integrand, p1.get_f(), - p2.get_f()); + // multiply basis functions with integrand and add to intval + super_space.addScaledBasisInteraction(intval, integrand, p1.get_xi(), + p2.get_xi()); return; } diff --git a/Bembel/src/util/surfaceL2error.hpp b/Bembel/src/util/surfaceL2error.hpp index 3c9ae3012..e0f24c21e 100644 --- a/Bembel/src/util/surfaceL2error.hpp +++ b/Bembel/src/util/surfaceL2error.hpp @@ -35,8 +35,8 @@ double surfaceL2error(const AnsatzSpace &ansatz_space, const Scalar val = fun_val.evaluate(*element, qp)(0); retval += qp.get_surface_measure() * Q.w_(i) * element->get_h() * element->get_h() * - (functor(qp.get_f()) - val / element->get_h()) * - (functor(qp.get_f()) - val / element->get_h()); + (functor(qp.get_f()) - val) * + (functor(qp.get_f()) - val); } } return sqrt(retval); diff --git a/tests/test_Bernstein.cpp b/tests/test_Bernstein.cpp index d2ef27ba7..5b11a86f7 100644 --- a/tests/test_Bernstein.cpp +++ b/tests/test_Bernstein.cpp @@ -21,7 +21,7 @@ int main() { constexpr int P = Bembel::Constants::MaxP; VectorXd buffer(P + 1); VectorXd coefs = - VectorXd::LinSpaced(1, P + 1, P + 1) / static_cast(P + 1); + VectorXd::LinSpaced(P + 1, 1, P + 1) / static_cast(P + 1); for (int p = 0; p <= Bembel::Constants::MaxP; ++p) { for (auto x : Test::Constants::eq_points) { @@ -30,7 +30,7 @@ int main() { double result1 = buffer.dot(coefs); std::vector v = {x}; - double result2 = Spl::DeBoor(Eigen::MatrixXd(coefs.leftCols(p + 1)), + double result2 = Spl::DeBoor(MatrixXd(coefs.head(p + 1).transpose()), Spl::MakeBezierKnotVector(p + 1), v)(0); BEMBEL_TEST_IF(std::abs(result1 - result2) < @@ -42,11 +42,11 @@ int main() { for (int p = 1; p <= Bembel::Constants::MaxP; ++p) { for (auto x : Test::Constants::eq_points) { buffer.setZero(); - Basis::ShapeFunctionHandler::evalDerCoef(p, buffer, x); + Basis::ShapeFunctionHandler::evalDerBasis(p, buffer.data(), x); double result1 = buffer.dot(coefs); std::vector v = {x}; - double result2 = Spl::DeBoorDer(Eigen::MatrixXd(coefs.leftCols(p + 1)), + double result2 = Spl::DeBoorDer(MatrixXd(coefs.head(p + 1).transpose()), Spl::MakeBezierKnotVector(p + 1), v)(0); BEMBEL_TEST_IF(std::abs(result1 - result2) < @@ -54,6 +54,5 @@ int main() { } } - delete[] coefs; return 0; } From 380c1be2ef7b96a605420372f5cb231b14d337e0 Mon Sep 17 00:00:00 2001 From: jdoelz <48793870+jdoelz@users.noreply.github.com> Date: Wed, 20 Mar 2024 09:58:37 +0100 Subject: [PATCH 17/20] style guide --- Bembel/src/LinearForm/DirichletTrace.hpp | 3 ++- Bembel/src/LinearForm/DirichletTraceOnReferenceDomain.hpp | 6 +++--- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/Bembel/src/LinearForm/DirichletTrace.hpp b/Bembel/src/LinearForm/DirichletTrace.hpp index c40e4deae..43188c01b 100644 --- a/Bembel/src/LinearForm/DirichletTrace.hpp +++ b/Bembel/src/LinearForm/DirichletTrace.hpp @@ -39,7 +39,8 @@ class DirichletTrace : public LinearFormBase, Scalar> { const T &super_space, const SurfacePoint &p, Eigen::Matrix *intval) const { // integrand without basis functions - Scalar integrand = function_(p.get_f()) * p.get_surface_measure() * p.get_w(); + Scalar integrand = + function_(p.get_f()) * p.get_surface_measure() * p.get_w(); // multiply basis functions with integrand super_space.addScaledBasis(intval, integrand, p.get_xi()); return; diff --git a/Bembel/src/LinearForm/DirichletTraceOnReferenceDomain.hpp b/Bembel/src/LinearForm/DirichletTraceOnReferenceDomain.hpp index 4a588b604..96db5e22a 100644 --- a/Bembel/src/LinearForm/DirichletTraceOnReferenceDomain.hpp +++ b/Bembel/src/LinearForm/DirichletTraceOnReferenceDomain.hpp @@ -8,8 +8,8 @@ // source code is subject to the GNU General Public License version 3 and // provided WITHOUT ANY WARRANTY, see for further // information. -#ifndef BEMBEL_SRC_LINEARFORM_DIRICHLETTRACEONREFERENCEDOMAIN_H_ -#define BEMBEL_SRC_LINEARFORM_DIRICHLETTRACEONREFERENCEDOMAIN_H_ +#ifndef BEMBEL_SRC_LINEARFORM_DIRICHLETTRACEONREFERENCEDOMAIN_HPP_ +#define BEMBEL_SRC_LINEARFORM_DIRICHLETTRACEONREFERENCEDOMAIN_HPP_ namespace Bembel { @@ -53,4 +53,4 @@ class DirichletTraceOnReferenceDomain }; } // namespace Bembel -#endif // BEMBEL_SRC_LINEARFORM_DIRICHLETTRACEONREFERENCEDOMAIN_H_ +#endif // BEMBEL_SRC_LINEARFORM_DIRICHLETTRACEONREFERENCEDOMAIN_HPP_ From ed0b2e62b2b2f32b1d84a28082400d500d7abc6d Mon Sep 17 00:00:00 2001 From: jdoelz <48793870+jdoelz@users.noreply.github.com> Date: Fri, 20 Jun 2025 17:08:02 +0200 Subject: [PATCH 18/20] Update Bembel/src/Geometry/SurfacePoint.hpp Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- Bembel/src/Geometry/SurfacePoint.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Bembel/src/Geometry/SurfacePoint.hpp b/Bembel/src/Geometry/SurfacePoint.hpp index 7e124b45c..487727b3b 100644 --- a/Bembel/src/Geometry/SurfacePoint.hpp +++ b/Bembel/src/Geometry/SurfacePoint.hpp @@ -109,7 +109,7 @@ class SurfacePoint { Eigen::Matrix jacobian_; bool buffer_is_alloced_ = false; int buffer_size_ = 0; - double *buffer_; + double *buffer_ = nullptr; }; typedef std::vector> From 3eff3fcd82ddf708eaae67831fe41fc2036853e9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20D=C3=B6lz?= Date: Fri, 20 Jun 2025 17:24:46 +0200 Subject: [PATCH 19/20] fix compilation issues --- .../Laplace/AdjointDoubleLayerOperator.hpp | 58 +++---------------- Bembel/src/Laplace/DoubleLayerOperator.hpp | 4 +- Bembel/src/Laplace/HypersingularOperator.hpp | 48 ++------------- 3 files changed, 17 insertions(+), 93 deletions(-) diff --git a/Bembel/src/Laplace/AdjointDoubleLayerOperator.hpp b/Bembel/src/Laplace/AdjointDoubleLayerOperator.hpp index 332496778..0901dd225 100644 --- a/Bembel/src/Laplace/AdjointDoubleLayerOperator.hpp +++ b/Bembel/src/Laplace/AdjointDoubleLayerOperator.hpp @@ -43,65 +43,25 @@ class LaplaceAdjointDoubleLayerOperator Eigen::Matrix::Scalar, Eigen::Dynamic, Eigen::Dynamic> *intval) const { - auto polynomial_degree = super_space.get_polynomial_degree(); - auto polynomial_degree_plus_one_squared = - (polynomial_degree + 1) * (polynomial_degree + 1); - - // get evaluation points on unit square - auto s = p1.segment<2>(0); - auto t = p2.segment<2>(0); - - // get quadrature weights - auto ws = p1(2); - auto wt = p2(2); - - // get points on geometry and tangential derivatives - auto x_f = p1.segment<3>(3); - auto x_f_dx = p1.segment<3>(6); - auto x_f_dy = p1.segment<3>(9); - auto y_f = p2.segment<3>(3); - auto y_f_dx = p2.segment<3>(6); - auto y_f_dy = p2.segment<3>(9); - - // compute surface measures from tangential derivatives - auto x_n = x_f_dx.cross(x_f_dy); - auto y_kappa = y_f_dx.cross(y_f_dy).norm(); - // integrand without basis functions - auto integrand = evaluateKernelGrad(x_f, x_n, y_f) * y_kappa * ws * wt; + auto integrand = + evaluateKernelGrad(p1.get_f(), p1.get_normal(), p2.get_f()) * + p2.get_surface_measure() * p1.get_w() * p2.get_w(); + + // multiply basis functions with integrand and add to intval + super_space.addScaledBasisInteraction(intval, integrand, p1.get_xi(), + p2.get_xi()); - // multiply basis functions with integrand and add to intval, this is an - // efficient implementation of - // (*intval) += super_space.BasisInteraction(s, t) * evaluateKernel(x_f, - // y_f) - // * x_kappa * y_kappa * ws * wt; - super_space.addScaledBasisInteraction(intval, integrand, s, t); return; } Eigen::Matrix evaluateFMMInterpolation_impl( const SurfacePoint &p1, const SurfacePoint &p2) const { - // get evaluation points on unit square - auto s = p1.segment<2>(0); - auto t = p2.segment<2>(0); - - // get points on geometry and tangential derivatives - auto x_f = p1.segment<3>(3); - auto x_f_dx = p1.segment<3>(6); - auto x_f_dy = p1.segment<3>(9); - auto y_f = p2.segment<3>(3); - auto y_f_dx = p2.segment<3>(6); - auto y_f_dy = p2.segment<3>(9); - - // compute surface measure in x and unnormalized normal in y from tangential - // derivatives - auto x_n = x_f_dx.cross(x_f_dy); - auto y_kappa = y_f_dx.cross(y_f_dy).norm(); - // interpolation Eigen::Matrix intval; - intval(0) = evaluateKernelGrad(x_f, x_n, y_f) * y_kappa; + intval(0) = evaluateKernelGrad(p1.get_f(), p1.get_normal(), p2.get_f()) * + p2.get_surface_measure(); return intval; } diff --git a/Bembel/src/Laplace/DoubleLayerOperator.hpp b/Bembel/src/Laplace/DoubleLayerOperator.hpp index ff1398ce8..97cbef8c4 100644 --- a/Bembel/src/Laplace/DoubleLayerOperator.hpp +++ b/Bembel/src/Laplace/DoubleLayerOperator.hpp @@ -49,8 +49,8 @@ class LaplaceDoubleLayerOperator p1.get_surface_measure() * p1.get_w() * p2.get_w(); // multiply basis functions with integrand and add to intval - super_space.addScaledBasisInteraction(intval, integrand, p1.get_f(), - p2.get_f()); + super_space.addScaledBasisInteraction(intval, integrand, p1.get_xi(), + p2.get_xi()); return; } diff --git a/Bembel/src/Laplace/HypersingularOperator.hpp b/Bembel/src/Laplace/HypersingularOperator.hpp index 0a5918164..c0b82122c 100644 --- a/Bembel/src/Laplace/HypersingularOperator.hpp +++ b/Bembel/src/Laplace/HypersingularOperator.hpp @@ -42,36 +42,12 @@ class LaplaceHypersingularOperator Eigen::Matrix< typename LinearOperatorTraits::Scalar, Eigen::Dynamic, Eigen::Dynamic> *intval) const { - auto polynomial_degree = super_space.get_polynomial_degree(); - auto polynomial_degree_plus_one_squared = - (polynomial_degree + 1) * (polynomial_degree + 1); - - // get evaluation points on unit square - auto s = p1.segment<2>(0); - auto t = p2.segment<2>(0); - - // get quadrature weights - auto ws = p1(2); - auto wt = p2(2); - - // get points on geometry and tangential derivatives - auto x_f = p1.segment<3>(3); - auto x_f_dx = p1.segment<3>(6); - auto x_f_dy = p1.segment<3>(9); - auto y_f = p2.segment<3>(3); - auto y_f_dx = p2.segment<3>(6); - auto y_f_dy = p2.segment<3>(9); - - // compute surface measures from tangential derivatives - auto x_kappa = x_f_dx.cross(x_f_dy).norm(); - auto y_kappa = y_f_dx.cross(y_f_dy).norm(); - // compute h auto h = 1. / (1 << super_space.get_refinement_level()); // h = 1 ./ (2^M) // integrand without basis functions auto integrand = - evaluateKernel(x_f, y_f) * x_kappa * y_kappa * ws * wt / h / h; + evaluateKernel(p1.get_f(), p2.get_f()) * p1.get_surface_measure() * p2.get_surface_measure() * p1.get_w() * p2.get_w() / h / h; // multiply basis functions with integrand and add to intval, this is an // efficient implementation of @@ -82,28 +58,16 @@ class LaplaceHypersingularOperator Eigen::Matrix evaluateFMMInterpolation_impl( const SurfacePoint &p1, const SurfacePoint &p2) const { - // get evaluation points on unit square - auto s = p1.segment<2>(0); - auto t = p2.segment<2>(0); - - // get points on geometry and tangential derivatives - auto x_f = p1.segment<3>(3); - auto x_f_dx = p1.segment<3>(6); - auto x_f_dy = p1.segment<3>(9); - auto y_f = p2.segment<3>(3); - auto y_f_dx = p2.segment<3>(6); - auto y_f_dy = p2.segment<3>(9); - // evaluate kernel - auto kernel = evaluateKernel(x_f, y_f); + auto kernel = evaluateKernel(p1.get_f(), p2.get_f()); // interpolation Eigen::Matrix intval; intval.setZero(); - intval(0, 0) = kernel * x_f_dy.dot(y_f_dy); - intval(0, 1) = -kernel * x_f_dy.dot(y_f_dx); - intval(1, 0) = -kernel * x_f_dx.dot(y_f_dy); - intval(1, 1) = kernel * x_f_dx.dot(y_f_dx); + intval(0, 0) = kernel * p1.get_f_dy().dot(p2.get_f_dy()); + intval(0, 1) = -kernel * p1.get_f_dy().dot(p2.get_f_dx()); + intval(1, 0) = -kernel * p1.get_f_dx().dot(p2.get_f_dy()); + intval(1, 1) = kernel * p1.get_f_dx().dot(p2.get_f_dx()); return intval; } From f824e79f1946b92cfcfdd1dbac9b5b37467ad472 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20D=C3=B6lz?= Date: Fri, 20 Jun 2025 17:27:12 +0200 Subject: [PATCH 20/20] fixed style guide issue --- Bembel/src/Laplace/HypersingularOperator.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/Bembel/src/Laplace/HypersingularOperator.hpp b/Bembel/src/Laplace/HypersingularOperator.hpp index c0b82122c..b8ca27eb1 100644 --- a/Bembel/src/Laplace/HypersingularOperator.hpp +++ b/Bembel/src/Laplace/HypersingularOperator.hpp @@ -46,8 +46,9 @@ class LaplaceHypersingularOperator auto h = 1. / (1 << super_space.get_refinement_level()); // h = 1 ./ (2^M) // integrand without basis functions - auto integrand = - evaluateKernel(p1.get_f(), p2.get_f()) * p1.get_surface_measure() * p2.get_surface_measure() * p1.get_w() * p2.get_w() / h / h; + auto integrand = evaluateKernel(p1.get_f(), p2.get_f()) * + p1.get_surface_measure() * p2.get_surface_measure() * + p1.get_w() * p2.get_w() / h / h; // multiply basis functions with integrand and add to intval, this is an // efficient implementation of