Function.h

Go to the documentation of this file.

// $Id$
//==============================================================================
//==============================================================================
 
#ifndef UTL_FUNCTION_H
#define UTL_FUNCTION_H
 
#include "matrixnd.h"
#include "Tensor.h"
#include "Vec3.h"
#include <cstddef>
 
 
namespace utl
{
  template<class Arg, class Result>
  class Function
  {
  protected:
    Function() {}
 
  public:
    virtual ~Function() {}
 
    virtual bool isZero() const { return false; }
    virtual bool isConstant() const { return true; }
 
  protected:
    virtual Result evaluate(const Arg& x) const = 0;
 
  public:
    Result operator()(const Arg& x) const { return this->evaluate(x); }
 
    typedef Arg    Input;  
    typedef Result Output; 
  };
 
 
  template<class Arg, class Result>
  class Function2
  {
  protected:
    Function2() {}
 
  public:
    virtual ~Function2() {}
 
    virtual bool isZero() const { return false; }
 
  protected:
    virtual Result evaluate(const Arg& x, const Arg& y) const = 0;
 
  public:
    Result operator()(const Arg& x, const Arg& y) const
    { return this->evaluate(x,y); }
 
    typedef Arg    Input;  
    typedef Result Output; 
  };
 
 
  template<class Result>
  class SpatialFunction : public Function<Vec3,Result>
  {
  protected:
    explicit SpatialFunction(const Result& val) : zero(val) {}
 
  public:
    virtual Result deriv(const Vec3&, int) const { return zero; }
    virtual Result dderiv(const Vec3&, int, int) const { return zero; }
 
  protected:
    virtual std::vector<Real> evalGradient(const Vec3&) const { return {}; }
    virtual std::vector<Real> evalHessian(const Vec3&) const { return {}; }
    virtual std::vector<Real> evalTimeDerivative(const Vec3&) const { return {}; }
 
    Result zero; 
  };
}
 
 
class ScalarFunc : public utl::Function<Real,Real>
{
protected:
  ScalarFunc() {}
 
public:
  Real eval(Real x) const { return this->evaluate(x); }
  virtual Real deriv(Real) const { return Real(0); }
};
 
 
class FunctionBase
{
protected:
  FunctionBase() : ncmp(1) {}
 
public:
  virtual ~FunctionBase() {}
 
  virtual unsigned char getType() const = 0;
 
  virtual std::vector<Real> getValue(const Vec3&) const = 0;
  virtual Real getScalarValue(const Vec3&) const = 0;
 
  size_t dim() const { return ncmp; }
 
  virtual bool initPatch(size_t) { return true; }
 
  virtual bool inDomain(const Vec3&) const { return true; }
 
  virtual void setParam(const std::string&, Real) {}
  void setParam(const std::string& name, const Vec3& value)
  {
    this->setParam(name + "x", value.x);
    this->setParam(name + "y", value.y);
    this->setParam(name + "z", value.z);
  }
 
protected:
  size_t ncmp; 
};
 
 
class RealFunc : public utl::SpatialFunction<Real>, public FunctionBase
{
protected:
  RealFunc() : utl::SpatialFunction<Real>(Real(0)) {}
 
public:
  unsigned char getType() const override { return 1; }
 
  std::vector<Real> getValue(const Vec3& X) const override
  {
    return std::vector<Real>(1,this->evaluate(X));
  }
 
  virtual Vec3 gradient(const Vec3& X) const
  {
    Vec3 result;
    for (size_t d = 1; d <= 3; ++d)
      result[d-1] = this->deriv(X,d);
 
    return result;
  }
 
  virtual SymmTensor hessian(const Vec3& X) const
  {
    SymmTensor result(3);
    for (size_t d1 = 1; d1 <= 3; ++d1)
      for (size_t d2 = d1; d2 <= 3; ++d2)
        result(d1,d2) = this->dderiv(X, d1, d2);
 
    return result;
  }
 
  Real getScalarValue(const Vec3& X) const override { return this->evaluate(X); }
 
  Real timeDerivative(const Vec3& X) const { return this->deriv(X,4); }
};
 
 
class VecFunc : public utl::SpatialFunction<Vec3>, public FunctionBase
{
protected:
  explicit VecFunc(size_t n = 3) : utl::SpatialFunction<Vec3>(Vec3())
  {
    ncmp = n;
  }
 
public:
  unsigned char getType() const override { return 2; }
 
  std::vector<Real> getValue(const Vec3& X) const override
  {
    return this->evaluate(X).vec(ncmp);
  }
 
  Real getScalarValue(const Vec3& X) const override
  {
    return this->evaluate(X).length();
  }
 
  Tensor gradient(const Vec3& X) const
  {
    Tensor result(ncmp);
    result = this->evalGradient(X);
    return result;
  }
 
  utl::matrix3d<Real> hessian(const Vec3& X) const
  {
    utl::matrix3d<Real> result(ncmp,ncmp,ncmp);
    result.fill(this->evalHessian(X).data());
    return result;
  }
 
  Vec3 timeDerivative(const Vec3& X) const { return this->evalTimeDerivative(X); }
};
 
 
class TractionFunc : public utl::Function2<Vec3,Vec3>
{
public:
  virtual bool isNormalPressure() const { return false; }
 
  virtual Vec3 timeDerivative(const Vec3&, const Vec3&) const { return Vec3(); }
};
 
#endif