app-fvdgfem/formula.hh

Go to the documentation of this file.

#ifndef fvdgFormula_hh
#define fvdgFormula_hh

#include "basics/typedefs.hh"
#include "formula/formula.hh"

#include <cmath>

using concepts::Real;

// **************************************************************** Formulae **

class BoxSin : public concepts::Formula<Real> {
public:
  BoxSin() : pi_(3.14159265359) {}
  concepts::Formula<Real> * clone() const { return new BoxSin(); }
  Real operator() (const Real x, const Real t = 0.) const { return 0.; }
  Real operator() (const concepts::Real2d& x, const Real t = 0.) const {
    return 0.; 
  }
  Real operator() (const concepts::Real3d& x, const Real t = 0.0) const {
    return 0.25*std::sin(pi_*(2*x[0]+2*x[1]+2*x[2]-3-6*t));
  }
private:
    const Real pi_;
};

class InitTest : public concepts::Formula<Real> {
public:
  InitTest(Real a, Real b, Real c) : a_(a), b_(b), c_(c) {} 
  concepts::Formula<Real> * clone() const { return new InitTest(a_, b_, c_); }
  Real operator() (const Real x, const Real t = 0.) const { return 0.; }
  Real operator() (const concepts::Real2d& x, const Real t = 0.) const {
    return 0.; 
  }
  Real operator() (const concepts::Real3d& x, const Real t = 0.0) const {
    return a_*x[0]+b_*x[1]+c_*x[2];
  }
private:
  Real a_;
  Real b_;
  Real c_;
};

class TimeSquared : public concepts::Formula<Real> {
public:
  concepts::Formula<Real> * clone() const { return new TimeSquared; }
  Real operator() (const Real x, const Real t = 0.) const { return t*t; }
  Real operator() (const concepts::Real2d& x, const Real t = 0.) const {
    return t*t; 
  }
  Real operator() (const concepts::Real3d& x, const Real t = 0.0) const {
    return t*t;
  }
};

class HeatSolution : public concepts::Formula<Real> {
public:
  HeatSolution() : pi_(3.14159265359) {}
  concepts::Formula<Real> * clone() const { return new HeatSolution; }
  Real operator() (const Real x, const Real t = 0.) const { return 0.; }
  Real operator() (const concepts::Real2d& x, const Real t = 0.) const {
    return 0.; 
  }
  Real operator() (const concepts::Real3d& x, const Real t = 0.0) const {
    return std::pow(4*pi_*(t+0.1),-1.5)*std::exp(-x.l2_2()/(4*(t+0.1)));
  }
private:
  Real pi_;
};

class ADLoad : public concepts::Formula<Real> {
public:
  ADLoad(Real epsilon, Real nu) : epsilon_(epsilon), nu_(nu) {}
  concepts::Formula<Real> * clone() const { return new ADLoad(epsilon_, nu_); }
  Real operator() (const Real x, const Real t = 0.) const { return 0.; }
  Real operator() (const concepts::Real2d& x, const Real t = 0.) const {
    return 0.;
  }
  Real operator() (const concepts::Real3d& x, const Real t = 0.) const {
    const Real z = (x[0]+x[1]+x[2]+3-3*t)/nu_;
    const concepts::Real3d x_2( 1 - x[0]*x[0], 1 - x[1]*x[1], 1 - x[2]*x[2] );
    const concepts::Real3d x_2c(x_2[1]*x_2[2], x_2[0]*x_2[2], x_2[0]*x_2[1]);
    Real pre = 1.;
    for (uint i = 0; i < 3; ++i) {
      pre *= x_2[i];
    }
    // Part from time derivative
    const Real a = 3*pre*pre/(nu_*(1+z*z));
    // Part from flux
    Real b = 0;
    for (uint i = 0; i < 3; ++i) {
      b += 2.*x[i]*x_2c[i];
    }
    b *= std::atan(z);
    b -= pre*3./(2.*nu_*(1+z*z));
    b *= -2*pre*pre*pre*std::atan(z);
    // Part from Laplace
    Real c = 0;
    for (uint i = 0; i < 3; ++i) {
      c += 4.*x_2c[i];
    }
    c *= pre;
    for (uint i = 0; i < 3; ++i) {
      c -= 8.*x[i]*x[i]*x_2c[i]*x_2c[i];
    }
    c *= std::atan(z);
    for (uint i = 0; i < 3; ++i) {
      c += 8.*x[i]*x_2c[i]*pre/(nu_*(1+z*z));
    }
    c += pre*pre*6*z/(nu_*nu_*std::pow(1+z*z,2));
    c *= -epsilon_;
    return a + b + c;
  }
private:
  const Real epsilon_;
  const Real nu_;
};

class ADSolution : public concepts::Formula<Real> {
public:
  ADSolution(Real nu) : nu_(nu) {}
  concepts::Formula<Real> * clone() const { return new ADSolution(nu_); }
  Real operator() (const Real x, const Real t = 0.) const { return 0.; }
  Real operator() (const concepts::Real2d& x, const Real t = 0.) const {
    return 0.;
  }
  Real operator() (const concepts::Real3d& x, const Real t = 0.) const {
    const Real z = (x[0]+x[1]+x[2]+3-3*t)/nu_;
    return
      -std::pow((1-x[0]*x[0])*(1-x[1]*x[1])*(1-x[2]*x[2]),2)*std::atan(z);
  }
private:
  const Real nu_;
};


#endif // fvdgFormula_hh