fluidselftransportmodel.hh

#ifndef __DUNE_ACFEM_MODELS_MODULES_FLUIDSELFTRANSPORTMODEL_HH__
#define __DUNE_ACFEM_MODELS_MODULES_FLUIDSELFTRANSPORTMODEL_HH__
 
#include <dune/fem/function/localfunction/const.hh>
 
#include "../../expressions/expressionoperations.hh"
#include "../../expressions/terminal.hh"
 
#include "../modelbase.hh"
 
namespace Dune {
 
  namespace ACFem::PDEModel {
 
    template<class FunctionSpace>
    class FluidSelfTransportModel
      : public ModelBase<FunctionSpace>
      , public Expressions::SelfExpression<FluidSelfTransportModel<FunctionSpace> >
      , public TypedValueExpression
    {
      using ThisType = FluidSelfTransportModel;
      using BaseType = ModelBase<FunctionSpace>;
     public:
      using typename BaseType::RangeFieldType;
      using typename BaseType::DomainType;
      using typename BaseType::RangeType;
      using typename BaseType::JacobianRangeType;
      using BaseType::dimRange;
      using BaseType::dimDomain;
 
      enum {
        dimWorld = dimDomain
      };
 
      static_assert(dimDomain == dimRange && dimDomain == dimWorld,
                    "This is meant for dimensionworld vector-fields, "
                    "like fluids, deformations etc.");
 
      // Interface methods that need to be reimplemented
 
      FluidSelfTransportModel(const std::string& name = "")
        : name_(name == "" ? "(-U_iU_jD_iPhi_j+Bndry)" : name)
      {}
 
      std::string name() const
      {
        return name_;
      }
 
      template<class Intersection>
      auto classifyBoundary(const Intersection& intersection)
      {
        // true is correct as higher-level code has to take care of
        // the Dirichlet-(non-Dirichlet) splitting of the boundary.
        return std::make_pair(true, std::bitset<dimRange>());
      }
 
      JacobianRangeType flux(const RangeType& value) const
      {
        // We need to provide the tensor-product of value with itself.
        JacobianRangeType flux(0);
        for (int i = 0; i < dimWorld; ++i) {
          flux[i][i] = value[i] * value[i];
          for (int j = i+1; j < dimWorld; ++j) {
            flux[i][j] = flux[j][i] = - (value[i] * value[j]);
          }
        }
        return flux;
      }
 
      template<class Point>
      JacobianRangeType linearizedFlux(const RangeType& uBar, const RangeType& value) const
      {
        // actually the sum of tensor products ...
        JacobianRangeType flux(0);
        for (int i = 0; i < dimWorld; ++i) {
          flux[i][i] = 2.0 * uBar[i] * value[i];
          for (int j = i+1; j < dimWorld; ++j) {
            flux[i][j] = flux[j][i] = - (uBar[i] * value[j] + uBar[j] * value[i]);
          }
        }
        return flux;
      }
 
      RangeType fluxDivergence(const RangeType& value, const JacobianRangeType& jacobian) const
      {
        // I think this is correct then ....
        RangeType result;
        jacobian.mv(value, result);
        return result;
      }
 
      RangeType robinFlux(const DomainType& unitOuterNormal,
                          const RangeType& value) const
      {
        RangeType result = value;
        return result *= (value * unitOuterNormal); // scalar product
      }
 
      RangeType linearizedRobinFlux(const RangeType& uBar,
                                    const DomainType& unitOuterNormal,
                                    const RangeType& value) const
      {
        RangeType result;
        RangeFieldType flowBar(0), flow(0);
        for (int i = 0; i < dimWorld; ++i) {
          flowBar += unitOuterNormal[i] * uBar[i];
          flow += unitOuterNormal[i] * value[i];
        }
        for (int i = 0; i < dimWorld; ++i) {
          result[i] = flow * uBar[i] + flowBar * value[i];
        }
        return result;
      }
 
     protected:
      std::string name_;
    };
 
 
    template<class Object>
    static inline auto
    fluidSelfTransportModel(const Object& object, const std::string& name = "")
    {
      return expressionClosure(FluidSelfTransportModel<typename Object::FunctionSpaceType>(name));
    }
 
 
 
 
  } // namespace ACFem::PDEModel
 
  namespace ACFem {
 
    using PDEModel::fluidSelfTransportModel;
 
  }
 
}  //Namespace Dune
 
 
#endif // __DUNE_ACFEM_MODELS_MODULES_FLUIDSELFTRANSPORTMODEL_HH__