flow_function.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2007 - 2018 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#ifndef dealii_flow_function_h
#define dealii_flow_function_h


#include <deal.II/base/config.h>

#include <deal.II/base/function.h>
#include <deal.II/base/point.h>
#include <deal.II/base/thread_management.h>

DEAL_II_NAMESPACE_OPEN

namespace Functions
{
  template <int dim>
  class FlowFunction : public Function<dim>
  {
  public:
    FlowFunction();

    virtual ~FlowFunction() override = default;

    void
    pressure_adjustment(double p);

    virtual void
    vector_values(const std::vector<Point<dim>> &   points,
                  std::vector<std::vector<double>> &values) const override = 0;
    virtual void
    vector_gradients(
      const std::vector<Point<dim>> &           points,
      std::vector<std::vector<Tensor<1, dim>>> &gradients) const override = 0;
    virtual void
    vector_laplacians(const std::vector<Point<dim>> &   points,
                      std::vector<std::vector<double>> &values) const = 0;

    virtual void
    vector_value(const Point<dim> &points,
                 Vector<double> &  value) const override;
    virtual double
    value(const Point<dim> & points,
          const unsigned int component) const override;
    virtual void
    vector_value_list(const std::vector<Point<dim>> &points,
                      std::vector<Vector<double>> &  values) const override;
    virtual void
    vector_gradient_list(
      const std::vector<Point<dim>> &           points,
      std::vector<std::vector<Tensor<1, dim>>> &gradients) const override;
    virtual void
    vector_laplacian_list(const std::vector<Point<dim>> &points,
                          std::vector<Vector<double>> &  values) const override;

    std::size_t
    memory_consumption() const;

  protected:
    double mean_pressure;

  private:
    mutable Threads::Mutex mutex;

    mutable std::vector<std::vector<double>> aux_values;

    mutable std::vector<std::vector<Tensor<1, dim>>> aux_gradients;
  };

  template <int dim>
  class PoisseuilleFlow : public FlowFunction<dim>
  {
  public:
    PoisseuilleFlow<dim>(const double r, const double Re);

    virtual ~PoisseuilleFlow() override = default;

    virtual void
    vector_values(const std::vector<Point<dim>> &   points,
                  std::vector<std::vector<double>> &values) const override;
    virtual void
    vector_gradients(
      const std::vector<Point<dim>> &           points,
      std::vector<std::vector<Tensor<1, dim>>> &gradients) const override;
    virtual void
    vector_laplacians(const std::vector<Point<dim>> &   points,
                      std::vector<std::vector<double>> &values) const override;

  private:
    const double radius;
    const double Reynolds;
  };


  template <int dim>
  class StokesCosine : public FlowFunction<dim>
  {
  public:
    StokesCosine(const double viscosity = 1., const double reaction = 0.);
    void
    set_parameters(const double viscosity, const double reaction);

    virtual ~StokesCosine() override = default;

    virtual void
    vector_values(const std::vector<Point<dim>> &   points,
                  std::vector<std::vector<double>> &values) const override;
    virtual void
    vector_gradients(
      const std::vector<Point<dim>> &           points,
      std::vector<std::vector<Tensor<1, dim>>> &gradients) const override;
    virtual void
    vector_laplacians(const std::vector<Point<dim>> &   points,
                      std::vector<std::vector<double>> &values) const override;

  private:
    double viscosity;
    double reaction;
  };


  class StokesLSingularity : public FlowFunction<2>
  {
  public:
    StokesLSingularity();

    virtual void
    vector_values(const std::vector<Point<2>> &     points,
                  std::vector<std::vector<double>> &values) const override;
    virtual void
    vector_gradients(
      const std::vector<Point<2>> &           points,
      std::vector<std::vector<Tensor<1, 2>>> &gradients) const override;
    virtual void
    vector_laplacians(const std::vector<Point<2>> &     points,
                      std::vector<std::vector<double>> &values) const override;

  private:
    double
    Psi(double phi) const;
    double
    Psi_1(double phi) const;
    double
    Psi_2(double phi) const;
    double
    Psi_3(double phi) const;
    double
    Psi_4(double phi) const;
    const double omega;
    static const double lambda;
    const double coslo;
    const double lp;
    const double lm;
  };

  class Kovasznay : public FlowFunction<2>
  {
  public:
    Kovasznay(const double Re, bool Stokes = false);

    virtual ~Kovasznay() override = default;

    virtual void
    vector_values(const std::vector<Point<2>> &     points,
                  std::vector<std::vector<double>> &values) const override;
    virtual void
    vector_gradients(
      const std::vector<Point<2>> &           points,
      std::vector<std::vector<Tensor<1, 2>>> &gradients) const override;
    virtual void
    vector_laplacians(const std::vector<Point<2>> &     points,
                      std::vector<std::vector<double>> &values) const override;

    double
    lambda() const;

  private:
    const double Reynolds;
    double       lbda;
    double       p_average;
    const bool   stokes;
  };

} // namespace Functions

DEAL_II_NAMESPACE_CLOSE

#endif