diff --git a/Bembel/LinearForm b/Bembel/LinearForm index 836b657f0..c7aa50c06 100644 --- a/Bembel/LinearForm +++ b/Bembel/LinearForm @@ -28,5 +28,6 @@ #include "src/LinearForm/DiscreteLinearForm.hpp" #include "src/LinearForm/RotatedTangentialTrace.hpp" #include "src/LinearForm/TangentialTrace.hpp" +#include "src/LinearForm/ReactionTerm.hpp" #endif // BEMBEL_LINEARFORM_MODULE_ diff --git a/Bembel/src/LinearForm/ReactionTerm.hpp b/Bembel/src/LinearForm/ReactionTerm.hpp new file mode 100644 index 000000000..6893495a9 --- /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_ diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 20f1b4a07..86eb1e554 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -20,6 +20,7 @@ set(CIFILES LaplaceDoubleLayerH2 LaplaceAdjointDoubleLayerH2 LazyEigenSum + GrayScott HelmholtzSingleLayerFull HelmholtzSingleLayerH2 HelmholtzDoubleLayerH2 diff --git a/examples/GrayScott.cpp b/examples/GrayScott.cpp new file mode 100644 index 000000000..0b1811099 --- /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; +}