rhs_function.h

Go to the documentation of this file.

 
#ifndef rhs_function_h
#define rhs_function_h
 
#include <deal.II/base/point.h>
#include <deal.II/base/function.h>
#include <deal.II/base/geometric_utilities.h>
#include <deal.II/base/tensor.h>
#include <cmath>
 
using namespace dealii;
 
template <int dim>
class RHSFunction : public Function<dim>
{
public:
    enum FunctionType {
        TYPE_X_Y,     
        TYPE_EXP_COS, 
        TYPE_ZERO,    
        TYPE_SIN_SIN, 
        TYPE_R_COS    
    };
 
    RHSFunction(FunctionType type)
        : Function<dim>(1), function_type(type) {}
 
    virtual double value(const Point<dim> & p,
                         const unsigned int component = 0) const override;
 
    virtual Tensor<1, dim> gradient(const Point<dim> & p,
                                    const unsigned int component = 0) const override;
 
private:
    FunctionType function_type;
 
    int R = 2;
};
 
// Out-of-line implementation of RHSFunction::value() (see class
// description for documentation). Dispatch is done with a chain of
// `if/else if` branches; the default branch returns 0.
template <int dim>
double RHSFunction<dim>::value(const Point<dim> &p, const unsigned int /*comp*/) const
{
    double return_value;
    if (function_type == TYPE_ZERO) {
        return 0;
    }
    else if (function_type == TYPE_SIN_SIN) {
        // Alternative right-hand side scaled with pi (commented out):
        //double pi = M_PI;
        //return 2 * pi * pi * std::sin(pi * p[0]) * std::sin(pi * p[1]);
        return 2 * std::sin(p[0]) * std::sin(p[1]);
    }
    else if (function_type == TYPE_X_Y) {
        return -10 * p[0] * p[1];
    }
    else if (function_type == TYPE_R_COS) {
        const auto r = p.norm();
        const auto theta = std::atan(p[1]/p[0]);
        return (r - 1/r - (1/r)*(1/R)*(1/R)) * std::cos(theta/R);
    }
    return_value = 0;
    return return_value;
}
 
 
// Out-of-line implementation of RHSFunction::gradient() (see class
// description for documentation). Returns the analytic gradient of the
// function selected by FunctionType; the TYPE_R_COS branch is an
// alternative test case using polar coordinates that we did not end up
// using for the final manufactured solution but is retained for reference.
template <int dim>
Tensor<1, dim> RHSFunction<dim>::gradient(const Point<dim> &p,
                                          const unsigned int /*comp*/) const
{
    Tensor<1, dim> return_values;
    return_values = 0;
    if (function_type == TYPE_ZERO) {
        return_values[0] = 0;
        return_values[1] = 0;
    }
    else if (function_type == TYPE_X_Y) {
        return_values[0] = -10 * p[1];
        return_values[1] = -10 * p[0];
    }
    else if (function_type == TYPE_SIN_SIN) {
        double pi = M_PI;
        return_values[0] = 2 * pi * pi * pi * std::cos(pi * p[0]) * std::sin(pi * p[1]);
        return_values[1] = 2 * pi * pi * pi * std::sin(pi * p[0]) * std::cos(pi * p[1]);
        return return_values;
    }
    else if (function_type == TYPE_R_COS) {
        // RHS gradient of an alternative test case in polar coordinates.
        // We ended up using a different test case in the final program.
        const auto r = p.norm();
        const auto theta = std::atan(p[1]/p[0]);
        return_values[0] = (1 + (1/r)*(1/r)*(1 + (1/R)*(1/R))) * std::cos(theta/R);
        return_values[1] = (r - 1/r - (1/r) * (1/R) * (1/R)) * (-1) * std::sin(theta/R);
    }
    return return_values;
}
 
 
//template class RHSFunction<2>;
//template class RHSFunction<3>;
 
#endif