flow_function.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2007 - 2017 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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#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() = default;

    void pressure_adjustment(double p);

    virtual void vector_values (const std::vector<Point<dim> > &points,
                                std::vector<std::vector<double> > &values) const = 0;
    virtual void vector_gradients (const std::vector<Point<dim> >            &points,
                                   std::vector<std::vector<Tensor<1,dim> > > &gradients) const = 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;
    virtual double value (const Point<dim> &points, const unsigned int component) const;
    virtual void vector_value_list (const std::vector<Point<dim> > &points,
                                    std::vector<Vector<double> >   &values) const;
    virtual void vector_gradient_list (const std::vector<Point<dim> >            &points,
                                       std::vector<std::vector<Tensor<1,dim> > > &gradients) const;
    virtual void vector_laplacian_list (const std::vector<Point<dim> > &points,
                                        std::vector<Vector<double> >   &values) const;

    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() = default;

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

  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() = default;

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

  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;
    virtual void vector_gradients (const std::vector<Point<2> > &points,
                                   std::vector<std::vector<Tensor<1,2> > > &gradients) const;
    virtual void vector_laplacians (const std::vector<Point<2> > &points,
                                    std::vector<std::vector<double> >   &values) const;
  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() = default;

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

    double lambda () const;
  private:
    const double Reynolds;
    double lbda;
    double p_average;
    const bool stokes;
  };

}

DEAL_II_NAMESPACE_CLOSE

#endif