Line data    Source code

       1            0 : // Distributed under the MIT License.
       2              : // See LICENSE.txt for details.
       3              : 
       4              : #pragma once
       5              : 
       6              : #include <cstddef>
       7              : #include <pup.h>
       8              : #include <string>
       9              : #include <vector>
      10              : 
      11              : #include "DataStructures/Tensor/TypeAliases.hpp"
      12              : #include "Elliptic/BoundaryConditions/BoundaryCondition.hpp"
      13              : #include "Elliptic/BoundaryConditions/BoundaryConditionType.hpp"
      14              : #include "Options/Context.hpp"
      15              : #include "Options/String.hpp"
      16              : #include "Utilities/Gsl.hpp"
      17              : #include "Utilities/Serialization/CharmPupable.hpp"
      18              : #include "Utilities/TMPL.hpp"
      19              : 
      20              : /// \cond
      21              : class DataVector;
      22              : /// \endcond
      23              : 
      24              : namespace Poisson::BoundaryConditions {
      25              : 
      26              : /// Impose Robin boundary conditions \f$a u + b n_i \nabla^i u = c\f$. The
      27              : /// boundary condition is imposed as Neumann-type (i.e. on \f$n_i \nabla^i u\f$)
      28              : /// if \f$|b| > 0\f$ and as Dirichlet-type (i.e. on \f$u\f$) if \f$b = 0\f$.
      29              : template <size_t Dim>
      30            1 : class Robin : public elliptic::BoundaryConditions::BoundaryCondition<Dim> {
      31              :  private:
      32            0 :   using Base = elliptic::BoundaryConditions::BoundaryCondition<Dim>;
      33              : 
      34              :  public:
      35            0 :   static constexpr Options::String help =
      36              :       "Robin boundary conditions a * u + b * n_i grad(u)^i = c. The boundary "
      37              :       "condition is imposed as Neumann-type (i.e. on n_i grad(u)^i) if abs(b) "
      38              :       "> 0 and as Dirichlet-type (i.e. on u) if b = 0.";
      39              : 
      40            0 :   struct DirichletWeight {
      41            0 :     using type = double;
      42            0 :     static constexpr Options::String help = "The parameter 'a'";
      43              :   };
      44              : 
      45            0 :   struct NeumannWeight {
      46            0 :     using type = double;
      47            0 :     static constexpr Options::String help = "The parameter 'b'";
      48              :   };
      49              : 
      50            0 :   struct Constant {
      51            0 :     using type = double;
      52            0 :     static constexpr Options::String help = "The parameter 'c'";
      53              :   };
      54              : 
      55            0 :   using options = tmpl::list<DirichletWeight, NeumannWeight, Constant>;
      56              : 
      57            0 :   Robin() = default;
      58            0 :   Robin(const Robin&) = default;
      59            0 :   Robin& operator=(const Robin&) = default;
      60            0 :   Robin(Robin&&) = default;
      61            0 :   Robin& operator=(Robin&&) = default;
      62            0 :   ~Robin() = default;
      63              : 
      64              :   /// \cond
      65              :   explicit Robin(CkMigrateMessage* m) : Base(m) {}
      66              :   using PUP::able::register_constructor;
      67              :   WRAPPED_PUPable_decl_template(Robin);
      68              :   /// \endcond
      69              : 
      70            0 :   std::unique_ptr<domain::BoundaryConditions::BoundaryCondition> get_clone()
      71              :       const override {
      72              :     return std::make_unique<Robin>(*this);
      73              :   }
      74              : 
      75            0 :   Robin(double dirichlet_weight, double neumann_weight, double constant,
      76              :         const Options::Context& context = {});
      77              : 
      78            0 :   double dirichlet_weight() const { return dirichlet_weight_; }
      79            0 :   double neumann_weight() const { return neumann_weight_; }
      80            0 :   double constant() const { return constant_; }
      81              : 
      82            0 :   std::vector<elliptic::BoundaryConditionType> boundary_condition_types()
      83              :       const override {
      84              :     return {1, neumann_weight_ == 0.
      85              :                    ? elliptic::BoundaryConditionType::Dirichlet
      86              :                    : elliptic::BoundaryConditionType::Neumann};
      87              :   }
      88              : 
      89            0 :   using argument_tags = tmpl::list<>;
      90            0 :   using volume_tags = tmpl::list<>;
      91              : 
      92            0 :   void apply(gsl::not_null<Scalar<DataVector>*> field,
      93              :              gsl::not_null<Scalar<DataVector>*> n_dot_field_gradient,
      94              :              const tnsr::i<DataVector, Dim>& deriv_field) const;
      95              : 
      96            0 :   using argument_tags_linearized = tmpl::list<>;
      97            0 :   using volume_tags_linearized = tmpl::list<>;
      98              : 
      99            0 :   void apply_linearized(
     100              :       gsl::not_null<Scalar<DataVector>*> field_correction,
     101              :       gsl::not_null<Scalar<DataVector>*> n_dot_field_gradient_correction,
     102              :       const tnsr::i<DataVector, Dim>& deriv_field_correction) const;
     103              : 
     104            0 :   void pup(PUP::er& p) override;
     105              : 
     106              :  private:
     107            0 :   double dirichlet_weight_ = std::numeric_limits<double>::signaling_NaN();
     108            0 :   double neumann_weight_ = std::numeric_limits<double>::signaling_NaN();
     109            0 :   double constant_ = std::numeric_limits<double>::signaling_NaN();
     110              : };
     111              : 
     112              : template <size_t Dim>
     113            0 : bool operator==(const Robin<Dim>& lhs, const Robin<Dim>& rhs);
     114              : 
     115              : template <size_t Dim>
     116            0 : bool operator!=(const Robin<Dim>& lhs, const Robin<Dim>& rhs);
     117              : 
     118              : }  // namespace Poisson::BoundaryConditions