From f81f44518f35c247efbefe85079ae61007598d46 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20D=C3=B6lz?= <48793870+jdoelz@users.noreply.github.com> Date: Thu, 23 Jul 2026 18:11:59 +0200 Subject: [PATCH 1/3] Add reaction diffusion implementation --- Bembel/LinearForm | 33 +----- Bembel/src/LinearForm/ReactionTerm.hpp | 77 ++++++++++++++ examples/CMakeLists.txt | 56 +--------- examples/GrayScott.cpp | 142 +++++++++++++++++++++++++ 4 files changed, 226 insertions(+), 82 deletions(-) create mode 100644 Bembel/src/LinearForm/ReactionTerm.hpp create mode 100644 examples/GrayScott.cpp diff --git a/Bembel/LinearForm b/Bembel/LinearForm index 836b657f0..d61266aaf 100644 --- a/Bembel/LinearForm +++ b/Bembel/LinearForm @@ -1,32 +1,9 @@ -// This file is part of Bembel, the higher order C++ boundary element library. -// -// Copyright (C) 2022 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_MODULE_ -#define BEMBEL_LINEARFORM_MODULE_ - -/** - * \ingroup Modules - * \defgroup LinearForm LinearForm - * \brief The LinearForm template class must be specialized for the assembly of - *the right hand side. - **/ - -#include "Quadrature" - -#include "src/LinearForm/LinearForm.hpp" - -#include "src/LinearForm/DirichletTrace.hpp" -#include "src/LinearForm/NeumannTrace.hpp" +#include "src/LinearForm/Basis2D.hpp" +#include "src/LinearForm/Basis3D.hpp" +#include "src/LinearForm/BasisIntersection.hpp" #include "src/LinearForm/DiscreteLinearForm.hpp" #include "src/LinearForm/RotatedTangentialTrace.hpp" #include "src/LinearForm/TangentialTrace.hpp" +#include "src/LinearForm/ReactionTerm.hpp" -#endif // BEMBEL_LINEARFORM_MODULE_ +#endif // BEMBEL_LINEARFORM_MODULE_ \ No newline at end of file diff --git a/Bembel/src/LinearForm/ReactionTerm.hpp b/Bembel/src/LinearForm/ReactionTerm.hpp new file mode 100644 index 000000000..6988d7be4 --- /dev/null +++ b/Bembel/src/LinearForm/ReactionTerm.hpp @@ -0,0 +1,77 @@ +// This file is part of Bembel, the higher order C++ boundary element library. +// +// Copyright (C) 2022 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_SRC_LINEARFORM_REACTIONTERM_HPP_ +#define BEMBEL_SRC_LINEARFORM_REACTIONTERM_HPP_ + +namespace Bembel { +/** + * \ingroup LinearForm + * \brief This class takes care of the + * assembly of the linear form of the reaction term, i.e., uv^2. + */ +template +void reactionLinearFrom(const AnsatzSpace &ansatz_space, + const Eigen::Matrix &u, + const Eigen::Matrix &v, + const Functor &functor, + Eigen::Matrix *reaction_lf, + int deg = 4) { + GaussSquare GS; + auto Q = GS[deg]; + SurfacePoint qp; + const auto longu = (ansatz_space.get_transformation_matrix() * u).eval(); + const auto longv = (ansatz_space.get_transformation_matrix() * v).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(); + const unsigned int polynomial_degree = super_space.get_polynomial_degree(); + const unsigned n_shape_fun = + (polynomial_degree + 1) * (polynomial_degree + 1); + // compute linear form of reaction term + Eigen::Matrix reaction_lf_long = + Eigen::Matrix::Zero( + n_shape_fun * number_of_elements, 1); + for (auto i = 0; i < Q.w_.size(); ++i) { + const Eigen::Matrix basis_evaluated_at_pt = + super_space.basis(Q.xi_.col(i)); + for (auto element = et.cpbegin(); element != et.cpend(); ++element) { + 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 h = element->get_h(); + const Scalar u_val = + longu.segment(n_shape_fun * element->id_, n_shape_fun).transpose() * + basis_evaluated_at_pt; + const Scalar v_val = + longv.segment(n_shape_fun * element->id_, n_shape_fun).transpose() * + basis_evaluated_at_pt; + reaction_lf_long.block(n_shape_fun * element->id_, 0, n_shape_fun, 1) += + x_kappa * Q.w_(i) * functor(u_val / h, v_val / h) * h * + basis_evaluated_at_pt; + } + } + reaction_lf->derived() = + ansatz_space.get_transformation_matrix().transpose() * reaction_lf_long; +} + +} // namespace Bembel +#endif // BEMBEL_SRC_LINEARFORM_REACTIONTERM_HPP_ \ No newline at end of file diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 2e6ac53c9..58633a930 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -1,62 +1,10 @@ - -# copy geometry files for testing -configure_file(${CMAKE_SOURCE_DIR}/geo/sphere.dat ${CMAKE_CURRENT_BINARY_DIR}/ COPYONLY) -configure_file(${CMAKE_SOURCE_DIR}/geo/torus.dat ${CMAKE_CURRENT_BINARY_DIR}/ COPYONLY) -configure_file(${CMAKE_SOURCE_DIR}/geo/cube_small.dat ${CMAKE_CURRENT_BINARY_DIR}/ COPYONLY) - -############################################################################### -# examples being run from CI -############################################################################### - set(CIFILES - AnsatzSpace - BlockClusterTree - Geometry - MassMatrix - LaplaceBeltrami - LaplaceBeltramiMultiGrid LaplaceSingleLayerFull LaplaceSingleLayerH2 - LaplaceHypersingularH2 LaplaceDoubleLayerH2 LaplaceAdjointDoubleLayerH2 LazyEigenSum + GrayScott HelmholtzSingleLayerFull HelmholtzSingleLayerH2 - HelmholtzDoubleLayerH2 - HelmholtzAdjointDoubleLayerH2 - HelmholtzHypersingularH2 - HomogenisedLaplaceSingleLayerFull - HomogenisedLaplaceSingleLayerH2 - MaxwellSingleLayerFull - MaxwellSingleLayerH2 - ) - -foreach(file IN LISTS CIFILES) - add_executable(${file} ${file}.cpp) - target_link_libraries(${file} Eigen3::Eigen) - add_test(${file}Example ${file}) - set_tests_properties(${file}Example PROPERTIES LABELS "E2Etests") -endforeach() - -enable_testing() - -add_custom_target(E2Etests) -add_dependencies(E2Etests ${CIFILES}) - -############################################################################### -# examples not being run from CI -############################################################################### - -set(NOCIFILES - Quadrature - VTKSurfaceExport - VTKDomainExport - FullLaplaceWorkflow - FullMaxwellWorkflow - ) - - foreach(file IN LISTS NOCIFILES) - add_executable(${file} ${file}.cpp) - target_link_libraries(${file} Eigen3::Eigen) -endforeach() + HelmholtzDoubleLayerH2 \ No newline at end of file diff --git a/examples/GrayScott.cpp b/examples/GrayScott.cpp new file mode 100644 index 000000000..9a66ee9c1 --- /dev/null +++ b/examples/GrayScott.cpp @@ -0,0 +1,142 @@ +// This file is part of Bembel, the higher order C++ boundary element library. +// +// Copyright (C) 2022 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 + +int main() { + using namespace Bembel; + using namespace Eigen; + IO::Stopwatch sw; + std::function init_u = [](const Vector3d &in) { + // small perturbation + return 0.6581 + 0.01 * (double)rand() / RAND_MAX - 0.005; + // return 0.6581 + 0.005; + // return 1.0; + }; + + std::function init_v = [](const Vector3d &in) { + // small perturbation + return 0.2279 + 0.01 * (double)rand() / RAND_MAX - 0.005; + // return 0.2279 - 0.005; + }; + + std::function init_cf = [](const Vector3d &in) { + return 1.0; + }; + + std::function reaction = + [](const double &u_val, const double &v_val) { + return u_val * v_val * v_val; + }; + + Eigen::VectorXd u_init; + Eigen::VectorXd v_init; + + Geometry geometry("sphere.dat"); + std::cout << "\n" << std::string(60, '=') << "\n"; + // Iterate over polynomial degree. + unsigned int polynomial_degree = 2; + // Iterate over refinement levels + unsigned int refinement_level = 4; + std::cout << "Degree " << polynomial_degree << " Level " << refinement_level + << "\n"; + sw.tic(); + // Build ansatz space + AnsatzSpace ansatz_space_lb( + geometry, refinement_level, polynomial_degree); + + AnsatzSpace ansatz_space_mass( + geometry, refinement_level, polynomial_degree); + // Gray Scott Model + // u_t = r_u * \laplacian u -uv^2 + f(1-u) + // v_y = r_v * \laplacian v +uv^2 - (f+k)v + double delta_t = 0.1; + double diffusion_rate_u = 0.01; + double diffusion_rate_v = + 0.001; // for sake of spots on the surface, please set the + // diffusion_rate_v to 0.0005 or less + double f = 0.1; + double k = 0.05; + // Build discrete operator for Laplace Beltrami and Identity + DiscreteLocalOperator disc_op_lb(ansatz_space_lb); + disc_op_lb.compute(); + DiscreteLocalOperator disc_op_mass(ansatz_space_mass); + disc_op_mass.compute(); + // Initial solutions of u and v + SparseLU, COLAMDOrdering> solver; + solver.compute(disc_op_mass.get_discrete_operator()); + DiscreteLinearForm, MassMatrixScalarCont> disc_lf( + ansatz_space_mass); + disc_lf.get_linear_form().set_function(init_u); + disc_lf.compute(); + u_init = solver.solve(disc_lf.get_discrete_linear_form()); + disc_lf.get_linear_form().set_function(init_v); + disc_lf.compute(); + v_init = solver.solve(disc_lf.get_discrete_linear_form()); + disc_lf.get_linear_form().set_function(init_cf); + disc_lf.compute(); + // Left hand side + const Eigen::SparseMatrix &stiff_mat = + disc_op_lb.get_discrete_operator(); + const Eigen::SparseMatrix &mass_mat = + disc_op_mass.get_discrete_operator(); + + Eigen::SimplicialLLT> solver_u; + Eigen::SimplicialLLT> solver_v; + solver_u.compute((1 + delta_t * f) * mass_mat + + delta_t * diffusion_rate_u * stiff_mat); + solver_v.compute((1 + (f + k) * delta_t) * mass_mat + + delta_t * diffusion_rate_v * stiff_mat); + sw.tic(); + Eigen::VectorXd u = u_init; + Eigen::VectorXd v = v_init; + // Solve the solutions of u and v at the end of the time + for (int i = 0; i < 20000; ++i) { + // Right hand side + Eigen::VectorXd reaction_term; + reactionLinearFrom(ansatz_space_lb, u, v, reaction, &reaction_term); + Eigen::Matrix rhs_u = + (delta_t * f) * disc_lf.get_discrete_linear_form() + + disc_op_mass.get_discrete_operator() * u - delta_t * reaction_term; + + Eigen::Matrix rhs_v = + disc_op_mass.get_discrete_operator() * v + delta_t * reaction_term; + u = solver_u.solve(rhs_u); + v = solver_v.solve(rhs_v); + } + auto difft = sw.toc(); + std::cout << "RD is solved in " << difft << " sec.\n"; + + // Visualization + FunctionEvaluator evaluator(ansatz_space_mass); + evaluator.set_function(u); + VTKSurfaceExport writer(geometry, 6); + std::function fun1 = + [&](int patch_number, const Eigen::Vector2d &reference_domain_point) { + auto retval = + evaluator.evaluateOnPatch(patch_number, reference_domain_point); + return double(retval(0, 0)); + }; + writer.addDataSet("u", fun1); + evaluator.set_function(v); + writer.addDataSet("v", fun1); + writer.writeToFile("sol.vtp"); + std::cout << "\n" << std::string(60, '=') << "\n"; + + return 0; +} \ No newline at end of file From 10d4e535058cf9c3d3b8bede1c14a274b38963ad Mon Sep 17 00:00:00 2001 From: jdoelz <48793870+jdoelz@users.noreply.github.com> Date: Fri, 24 Jul 2026 08:22:44 +0200 Subject: [PATCH 2/3] minor corrections and typos --- Bembel/LinearForm | 34 ++++++++++++--- Bembel/src/LinearForm/ReactionTerm.hpp | 7 +--- examples/CMakeLists.txt | 57 +++++++++++++++++++++++++- examples/GrayScott.cpp | 6 +-- 4 files changed, 87 insertions(+), 17 deletions(-) diff --git a/Bembel/LinearForm b/Bembel/LinearForm index d61266aaf..ffb76b4d4 100644 --- a/Bembel/LinearForm +++ b/Bembel/LinearForm @@ -1,9 +1,33 @@ -#include "src/LinearForm/Basis2D.hpp" -#include "src/LinearForm/Basis3D.hpp" -#include "src/LinearForm/BasisIntersection.hpp" +// This file is part of Bembel, the higher order C++ boundary element library. +// +// Copyright (C) 2022 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_MODULE_ +#define BEMBEL_LINEARFORM_MODULE_ + +/** + * \ingroup Modules + * \defgroup LinearForm LinearForm + * \brief The LinearForm template class must be specialized for the assembly of + *the right hand side. + **/ + +#include "Quadrature" + +#include "src/LinearForm/LinearForm.hpp" + +#include "src/LinearForm/DirichletTrace.hpp" +#include "src/LinearForm/NeumannTrace.hpp" #include "src/LinearForm/DiscreteLinearForm.hpp" +#include "src/LinearForm/ReactionTerm.hpp" #include "src/LinearForm/RotatedTangentialTrace.hpp" #include "src/LinearForm/TangentialTrace.hpp" -#include "src/LinearForm/ReactionTerm.hpp" -#endif // BEMBEL_LINEARFORM_MODULE_ \ No newline at end of file +#endif // BEMBEL_LINEARFORM_MODULE_ diff --git a/Bembel/src/LinearForm/ReactionTerm.hpp b/Bembel/src/LinearForm/ReactionTerm.hpp index 6988d7be4..72a5ea98d 100644 --- a/Bembel/src/LinearForm/ReactionTerm.hpp +++ b/Bembel/src/LinearForm/ReactionTerm.hpp @@ -19,7 +19,7 @@ namespace Bembel { * assembly of the linear form of the reaction term, i.e., uv^2. */ template -void reactionLinearFrom(const AnsatzSpace &ansatz_space, +void reactionLinearForm(const AnsatzSpace &ansatz_space, const Eigen::Matrix &u, const Eigen::Matrix &v, const Functor &functor, @@ -45,12 +45,7 @@ void reactionLinearFrom(const AnsatzSpace &ansatz_space, super_space.basis(Q.xi_.col(i)); for (auto element = et.cpbegin(); element != et.cpend(); ++element) { 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); diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 58633a930..c44091433 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -1,10 +1,63 @@ + +# copy geometry files for testing +configure_file(${CMAKE_SOURCE_DIR}/geo/sphere.dat ${CMAKE_CURRENT_BINARY_DIR}/ COPYONLY) +configure_file(${CMAKE_SOURCE_DIR}/geo/torus.dat ${CMAKE_CURRENT_BINARY_DIR}/ COPYONLY) +configure_file(${CMAKE_SOURCE_DIR}/geo/cube_small.dat ${CMAKE_CURRENT_BINARY_DIR}/ COPYONLY) + +############################################################################### +# examples being run from CI +############################################################################### + set(CIFILES + AnsatzSpace + BlockClusterTree + Geometry + MassMatrix + LaplaceBeltrami + LaplaceBeltramiMultiGrid LaplaceSingleLayerFull LaplaceSingleLayerH2 + LaplaceHypersingularH2 LaplaceDoubleLayerH2 LaplaceAdjointDoubleLayerH2 LazyEigenSum - GrayScott HelmholtzSingleLayerFull HelmholtzSingleLayerH2 - HelmholtzDoubleLayerH2 \ No newline at end of file + HelmholtzDoubleLayerH2 + HelmholtzAdjointDoubleLayerH2 + HelmholtzHypersingularH2 + HomogenisedLaplaceSingleLayerFull + HomogenisedLaplaceSingleLayerH2 + MaxwellSingleLayerFull + MaxwellSingleLayerH2 + ) + +foreach(file IN LISTS CIFILES) + add_executable(${file} ${file}.cpp) + target_link_libraries(${file} Eigen3::Eigen) + add_test(${file}Example ${file}) + set_tests_properties(${file}Example PROPERTIES LABELS "E2Etests") +endforeach() + +enable_testing() + +add_custom_target(E2Etests) +add_dependencies(E2Etests ${CIFILES}) + +############################################################################### +# examples not being run from CI +############################################################################### + +set(NOCIFILES + GrayScott + Quadrature + VTKSurfaceExport + VTKDomainExport + FullLaplaceWorkflow + FullMaxwellWorkflow + ) + + foreach(file IN LISTS NOCIFILES) + add_executable(${file} ${file}.cpp) + target_link_libraries(${file} Eigen3::Eigen) +endforeach() diff --git a/examples/GrayScott.cpp b/examples/GrayScott.cpp index 9a66ee9c1..66518a3b6 100644 --- a/examples/GrayScott.cpp +++ b/examples/GrayScott.cpp @@ -16,8 +16,6 @@ #include #include -#include - int main() { using namespace Bembel; using namespace Eigen; @@ -64,7 +62,7 @@ int main() { geometry, refinement_level, polynomial_degree); // Gray Scott Model // u_t = r_u * \laplacian u -uv^2 + f(1-u) - // v_y = r_v * \laplacian v +uv^2 - (f+k)v + // v_t = r_v * \laplacian v +uv^2 - (f+k)v double delta_t = 0.1; double diffusion_rate_u = 0.01; double diffusion_rate_v = @@ -109,7 +107,7 @@ int main() { for (int i = 0; i < 20000; ++i) { // Right hand side Eigen::VectorXd reaction_term; - reactionLinearFrom(ansatz_space_lb, u, v, reaction, &reaction_term); + reactionLinearForm(ansatz_space_lb, u, v, reaction, &reaction_term); Eigen::Matrix rhs_u = (delta_t * f) * disc_lf.get_discrete_linear_form() + disc_op_mass.get_discrete_operator() * u - delta_t * reaction_term; From ddeeaf97c65935d8e83cb301acc5a693df02e33f Mon Sep 17 00:00:00 2001 From: jdoelz <48793870+jdoelz@users.noreply.github.com> Date: Fri, 24 Jul 2026 08:47:36 +0200 Subject: [PATCH 3/3] moved GrayScott to unsupported --- Bembel/LinearForm | 1 - examples/CMakeLists.txt | 1 - unsupported/Bembel/GrayScott | 28 +++++++++++++++++++ .../Bembel/src/GrayScott}/ReactionTerm.hpp | 2 +- unsupported/examples/CMakeLists.txt | 1 + .../examples}/GrayScott.cpp | 19 +++++++------ 6 files changed, 40 insertions(+), 12 deletions(-) create mode 100644 unsupported/Bembel/GrayScott rename {Bembel/src/LinearForm => unsupported/Bembel/src/GrayScott}/ReactionTerm.hpp (98%) rename {examples => unsupported/examples}/GrayScott.cpp (88%) diff --git a/Bembel/LinearForm b/Bembel/LinearForm index ffb76b4d4..836b657f0 100644 --- a/Bembel/LinearForm +++ b/Bembel/LinearForm @@ -26,7 +26,6 @@ #include "src/LinearForm/DirichletTrace.hpp" #include "src/LinearForm/NeumannTrace.hpp" #include "src/LinearForm/DiscreteLinearForm.hpp" -#include "src/LinearForm/ReactionTerm.hpp" #include "src/LinearForm/RotatedTangentialTrace.hpp" #include "src/LinearForm/TangentialTrace.hpp" diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index c44091433..2e6ac53c9 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -49,7 +49,6 @@ add_dependencies(E2Etests ${CIFILES}) ############################################################################### set(NOCIFILES - GrayScott Quadrature VTKSurfaceExport VTKDomainExport diff --git a/unsupported/Bembel/GrayScott b/unsupported/Bembel/GrayScott new file mode 100644 index 000000000..1d8e39441 --- /dev/null +++ b/unsupported/Bembel/GrayScott @@ -0,0 +1,28 @@ +// This file is part of Bembel, the higher order C++ boundary element library. +// +// Copyright (C) 2022 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_GRAYSCOTT_MODULE_ +#define BEMBEL_GRAYSCOTT_MODULE_ + +/** + * \ingroup Modules + * \defgroup GrayScott GrayScott + * \brief Everything which is needed to solve the Gray-Scott model using + * Bembel. + **/ + +#include +#include +#include + +#include "src/GrayScott/ReactionTerm.hpp" + +#endif // BEMBEL_GRAYSCOTT_MODULE_ diff --git a/Bembel/src/LinearForm/ReactionTerm.hpp b/unsupported/Bembel/src/GrayScott/ReactionTerm.hpp similarity index 98% rename from Bembel/src/LinearForm/ReactionTerm.hpp rename to unsupported/Bembel/src/GrayScott/ReactionTerm.hpp index 72a5ea98d..870e43a96 100644 --- a/Bembel/src/LinearForm/ReactionTerm.hpp +++ b/unsupported/Bembel/src/GrayScott/ReactionTerm.hpp @@ -69,4 +69,4 @@ void reactionLinearForm(const AnsatzSpace &ansatz_space, } } // namespace Bembel -#endif // BEMBEL_SRC_LINEARFORM_REACTIONTERM_HPP_ \ No newline at end of file +#endif // BEMBEL_SRC_LINEARFORM_REACTIONTERM_HPP_ diff --git a/unsupported/examples/CMakeLists.txt b/unsupported/examples/CMakeLists.txt index 546a27c18..a8b76dc42 100644 --- a/unsupported/examples/CMakeLists.txt +++ b/unsupported/examples/CMakeLists.txt @@ -7,6 +7,7 @@ set(OTHER_NOCIFILES AugmentedEFIE AugmentedEFIEStabilized AugmentedEFIEBalun + GrayScott ) foreach(file IN LISTS OTHER_NOCIFILES) diff --git a/examples/GrayScott.cpp b/unsupported/examples/GrayScott.cpp similarity index 88% rename from examples/GrayScott.cpp rename to unsupported/examples/GrayScott.cpp index 66518a3b6..01185b3b7 100644 --- a/examples/GrayScott.cpp +++ b/unsupported/examples/GrayScott.cpp @@ -15,30 +15,31 @@ #include #include #include +#include int main() { using namespace Bembel; using namespace Eigen; IO::Stopwatch sw; - std::function init_u = [](const Vector3d &in) { + std::function init_u = [](const Vector3d& in) { // small perturbation return 0.6581 + 0.01 * (double)rand() / RAND_MAX - 0.005; // return 0.6581 + 0.005; // return 1.0; }; - std::function init_v = [](const Vector3d &in) { + std::function init_v = [](const Vector3d& in) { // small perturbation return 0.2279 + 0.01 * (double)rand() / RAND_MAX - 0.005; // return 0.2279 - 0.005; }; - std::function init_cf = [](const Vector3d &in) { + std::function init_cf = [](const Vector3d& in) { return 1.0; }; - std::function reaction = - [](const double &u_val, const double &v_val) { + std::function reaction = + [](const double& u_val, const double& v_val) { return u_val * v_val * v_val; }; @@ -89,9 +90,9 @@ int main() { disc_lf.get_linear_form().set_function(init_cf); disc_lf.compute(); // Left hand side - const Eigen::SparseMatrix &stiff_mat = + const Eigen::SparseMatrix& stiff_mat = disc_op_lb.get_discrete_operator(); - const Eigen::SparseMatrix &mass_mat = + const Eigen::SparseMatrix& mass_mat = disc_op_mass.get_discrete_operator(); Eigen::SimplicialLLT> solver_u; @@ -124,8 +125,8 @@ int main() { FunctionEvaluator evaluator(ansatz_space_mass); evaluator.set_function(u); VTKSurfaceExport writer(geometry, 6); - std::function fun1 = - [&](int patch_number, const Eigen::Vector2d &reference_domain_point) { + std::function fun1 = + [&](int patch_number, const Eigen::Vector2d& reference_domain_point) { auto retval = evaluator.evaluateOnPatch(patch_number, reference_domain_point); return double(retval(0, 0));