diff --git a/src/Utility/Functions.C b/src/Utility/Functions.C index 848b23dd..989a7ca4 100644 --- a/src/Utility/Functions.C +++ b/src/Utility/Functions.C @@ -23,27 +23,41 @@ PressureField::PressureField (real p, int dir) : pdir(dir) } -real ConstTimeFunc::evaluate (const Vec3& x) const +real RampFunc::evaluate (const real& x) const { - const Vec4* X = dynamic_cast(&x); - return (*tfunc)(X ? X->t : real(0)); + return x < xmax ? fval*x/xmax : fval; } -real SpaceTimeFunc::evaluate (const Vec3& x) const +real DiracFunc::evaluate (const real& x) const { - const Vec4* X = dynamic_cast(&x); - return (*sfunc)(x) * (*tfunc)(X ? X->t : real(0)); + return fabs(x-xmax) < 1.0e-4 ? amp : real(0); } -real LinearTinitFunc::evaluate (const Vec3& x) const +real StepFunc::evaluate (const real& x) const { - const Vec4* X = dynamic_cast(&x); - if (X && X->t < Tinit) - return value*X->t/Tinit; - else - return value; + return x >= xmax ? amp : real(0); +} + + +real SineFunc::evaluate (const real& x) const +{ + return scale*sin(freq*x+phase); +} + + +real ConstTimeFunc::evaluate (const Vec3& X) const +{ + const Vec4* Xt = dynamic_cast(&X); + return (*tfunc)(Xt ? Xt->t : real(0)); +} + + +real SpaceTimeFunc::evaluate (const Vec3& X) const +{ + const Vec4* Xt = dynamic_cast(&X); + return (*sfunc)(X) * (*tfunc)(Xt ? Xt->t : real(0)); } @@ -65,61 +79,70 @@ real LinearZFunc::evaluate (const Vec3& X) const } -real QuadraticXFunc::evaluate(const Vec3& X) const +real QuadraticXFunc::evaluate (const Vec3& X) const { - real val = 0.5*(a-b); + real val = (a-b)/real(2); return max*(a-X.x)*(X.x-b)/(val*val); } -real QuadraticYFunc::evaluate(const Vec3& X) const +real QuadraticYFunc::evaluate (const Vec3& X) const { - real val = 0.5*(a-b); + real val = (a-b)/real(2); return max*(a-X.y)*(X.y-b)/(val*val); } -real QuadraticZFunc::evaluate(const Vec3& X) const +real QuadraticZFunc::evaluate (const Vec3& X) const { - real val = 0.5*(a-b); + real val = (a-b)/real(2); return max*(a-X.z)*(X.z-b)/(val*val); } -real LinearRotZFunc::evaluate (const Vec3& _X) const +real LinearRotZFunc::evaluate (const Vec3& X) const { // Always return zero if the argument has no time component - const Vec4* X = dynamic_cast(&_X); - if (!X) return real(0); + const Vec4* Xt = dynamic_cast(&X); + if (!Xt) return real(0); - real x = X->x - x0; - real y = X->y - y0; - real c = cos(A*X->t); - real s = sin(A*X->t); + real x = X.x - x0; + real y = X.y - y0; + real c = cos(A*Xt->t); + real s = sin(A*Xt->t); return rX ? x*c-y*s-x : x*s+y*c-y; } real StepXFunc::evaluate (const Vec3& X) const { - return X.x < x0 || X.x > x1 ? 0.0 : fv; + return X.x < x0 || X.x > x1 ? real(0) : fv; } - real StepXYFunc::evaluate (const Vec3& X) const { - return X.x < x0 || X.x > x1 || X.y < y0 || X.y > y1 ? 0.0 : fv; + return X.x < x0 || X.x > x1 || X.y < y0 || X.y > y1 ? real(0) : fv; } +/*! + The functions is assumed on the general form + \f[ f({\bf X},t) = A * g({\bf X}) * h(t) \f] + + The character string \a cline is assumed to contain first the definition + of the spatial function \a g( \b X ) and then the time function \a h(t). + Either of the two components may be omitted, for creating a space-function + constant in time, or a time function constant in space. +*/ + const RealFunc* utl::parseRealFunc (char* cline, real A) { // Check for spatial variation int linear = 0; int quadratic = 0; if (!cline) - linear = -1; + return 0; else if (strcmp(cline,"X") == 0) linear = 1; else if (strcmp(cline,"Y") == 0) @@ -130,18 +153,16 @@ const RealFunc* utl::parseRealFunc (char* cline, real A) linear = 4; else if (strcmp(cline,"YrotZ") == 0) linear = 5; - else if (strcmp(cline,"Tinit") == 0) + else if (strcmp(cline,"StepX") == 0) linear = 6; + else if (strcmp(cline,"StepXY") == 0) + linear = 7; else if (strcmp(cline,"quadX") == 0) quadratic = 1; else if (strcmp(cline,"quadY") == 0) quadratic = 2; else if (strcmp(cline,"quadZ") == 0) quadratic = 3; - else if (strcmp(cline,"StepX") == 0) - linear = 7; - else if (strcmp(cline,"StepXY") == 0) - linear = 8; real C = A; const RealFunc* f = 0; @@ -170,29 +191,25 @@ const RealFunc* utl::parseRealFunc (char* cline, real A) f = new LinearRotZFunc(false,A,atof(cline),atof(strtok(NULL," "))); break; case 6: - std::cout <<"RampT("<< cline <<"))"; - f = new LinearTinitFunc(A,atof(cline)); - break; - case 7: std::cout <<"StepX("<< cline <<"))"; f = new StepXFunc(A,atof(cline),atof(strtok(NULL," "))); break; - case 8: + case 7: std::cout <<"StepXY("<< cline <<"))"; f = new StepXYFunc(A,atof(cline),atof(strtok(NULL," "))); break; } cline = strtok(NULL," "); } - else if (quadratic && (cline = strtok(NULL," "))) + else if (quadratic > 0 && (cline = strtok(NULL," "))) { + C = real(1); real a = atof(cline); real b = atof(strtok(NULL," ")); - real val = 0.5*(a-b); - val *= val; - std::cout << A/val << "*(" << char('W' + quadratic) << "-" << a << ")*(" - << b << "-" << char('W' + quadratic) << ")"; - switch(quadratic) { + real val = (a-b)*(a-b)/real(4); + char var = 'W' + quadratic; + std::cout << A/val <<" * ("<< a <<"-"<< var <<")*("<< b <<"-"<< var <<")"; + switch (quadratic) { case 1: f = new QuadraticXFunc(A,a,b); break; @@ -205,23 +222,51 @@ const RealFunc* utl::parseRealFunc (char* cline, real A) } cline = strtok(NULL," "); } - if (quadratic == 0) + else // constant in space std::cout << C; // Check for time variation - if (!cline) return f; + if (!cline) return f; // constant in time const ScalarFunc* s = 0; - double freq = atof(cline); - if ((cline = strtok(NULL," "))) + if (strncmp(cline,"Ramp",4) == 0 || strcmp(cline,"Tinit") == 0) { - std::cout <<" * sin("<< freq <<"*t + "<< cline <<")"; - s = new SineFunc(C,freq,atof(cline)); + real xmax = atof(strtok(NULL," ")); + std::cout <<" * Ramp(t,"<< xmax <<")"; + s = new RampFunc(C,xmax); } - else + else if (strncmp(cline,"Dirac",5) == 0) { - std::cout <<" * "<< freq <<"*t"; - s = new LinearFunc(C*freq); + real xmax = atof(strtok(NULL," ")); + std::cout <<" * Dirac(t,"<< xmax <<")"; + s = new DiracFunc(C,xmax); + } + else if (strncmp(cline,"Step",4) == 0) + { + real xmax = atof(strtok(NULL," ")); + std::cout <<"Step(t,"<< xmax <<")"; + s = new StepFunc(C,xmax); + } + else if (strcmp(cline,"sin") == 0) + { + real freq = atof(strtok(NULL," ")); + if ((cline = strtok(NULL," "))) + { + real phase = atof(cline); + std::cout <<" * sin("<< freq <<"*t + "<< phase <<")"; + s = new SineFunc(C,freq,phase); + } + else + { + std::cout <<" * sin("<< freq <<"*t)"; + s = new SineFunc(C,freq); + } + } + else // linear in time + { + real scale = atof(cline); + std::cout <<" * "<< scale <<"*t"; + s = new LinearFunc(C*scale); } if (f) diff --git a/src/Utility/Functions.h b/src/Utility/Functions.h index 7ac8686e..c9c246b9 100644 --- a/src/Utility/Functions.h +++ b/src/Utility/Functions.h @@ -15,11 +15,10 @@ #define _FUNCTIONS_H #include "Function.h" -#include /*! - \brief A linear scalar function. + \brief A scalar-valued linear function. */ class LinearFunc : public ScalarFunc @@ -31,13 +30,70 @@ public: LinearFunc(real s = real(1)) : scale(s) {} protected: - //! \brief Evaluates the scalar function. + //! \brief Evaluates the function at \a x. virtual real evaluate(const real& x) const { return scale*x; } }; /*! - \brief A sinusoidal scalar function. + \brief A scalar-valued ramp function, linear up to \a xmax. +*/ + +class RampFunc : public ScalarFunc +{ + real fval; //!< Max function value + real xmax; //!< Function is linear from \a x = 0 to \a x = \a xmax + +public: + //! \brief Constructor initializing the function parameters. + RampFunc(real f = real(1), real x = real(1)) : fval(f), xmax(x) {} + +protected: + //! \brief Evaluates the function at \a x. + virtual real evaluate(const real& x) const; +}; + + +/*! + \brief A scalar-valued dirac function. +*/ + +class DiracFunc : public ScalarFunc +{ + real amp; //!< The amplitude of the dirac function + real xmax; //!< Associated \a x value + +public: + //! \brief Constructor initializing the function parameters. + DiracFunc(real a = real(1), real x = real(0)) : amp(a), xmax(x) {} + +protected: + //! \brief Evaluates the function at \a x. + virtual real evaluate(const real& x) const; +}; + + +/*! + \brief A scalar-valued step function. +*/ + +class StepFunc : public ScalarFunc +{ + real amp; //!< The amplitude of the step function + real xmax; //!< Associated \a x value + +public: + //! \brief Constructor initializing the function parameters. + StepFunc(real a, real x = real(0)) : amp(a), xmax(x) {} + +protected: + //! \brief Evaluates the function at \a x. + virtual real evaluate(const real& x) const; +}; + + +/*! + \brief A scalar-valued sinusoidal function. */ class SineFunc : public ScalarFunc @@ -52,13 +108,13 @@ public: : scale(s), freq(f), phase(p) {} protected: - //! \brief Evaluates the scalar function. - virtual real evaluate(const real& x) const { return scale*sin(freq*x+phase); } + //! \brief Evaluates the function at \a x. + virtual real evaluate(const real& x) const; }; /*! - \brief A scalar function, constant in space and time. + \brief A scalar-valued spatial function, constant in space and time. */ class ConstFunc : public RealFunc @@ -76,7 +132,7 @@ protected: /*! - \brief A scalar function, constant in space but varying in time. + \brief A scalar-valued spatial function, constant in space, varying in time. */ class ConstTimeFunc : public RealFunc @@ -91,12 +147,12 @@ public: protected: //! \brief Evaluates the time-varying function. - virtual real evaluate(const Vec3& x) const; + virtual real evaluate(const Vec3& X) const; }; /*! - \brief A scalar function, varying in space and time. + \brief A scalar-valued spatial function, varying in space and time. \details The function value is defined as a product between one space-dependent component and one time-dependent component. */ @@ -114,31 +170,12 @@ public: protected: //! \brief Evaluates the space-time function. - virtual real evaluate(const Vec3& x) const; + virtual real evaluate(const Vec3& X) const; }; /*! - \brief A scalar function, linear in \a t up to \a Tinit. -*/ - -class LinearTinitFunc : public RealFunc -{ - real value; //!< Max function value - real Tinit; //!< Function is linear from t = 0 to t = Tinit - -public: - //! \brief Constructor initializing the function parameters. - LinearTinitFunc(real value_, real Tinit_) : value(value_), Tinit(Tinit_) {} - -protected: - //! \brief Evaluates the linear function. - virtual real evaluate(const Vec3& x) const; -}; - - -/*! - \brief A scalar function, linear in \a x. + \brief A scalar-valued spatial function, linear in \a x. */ class LinearXFunc : public RealFunc @@ -157,7 +194,7 @@ protected: /*! - \brief A scalar function, linear in \a y. + \brief A scalar-valued spatial function, linear in \a y. */ class LinearYFunc : public RealFunc @@ -176,7 +213,7 @@ protected: /*! - \brief A scalar function, linear in \a z. + \brief A scalar-valued spatial function, linear in \a z. */ class LinearZFunc : public RealFunc @@ -195,7 +232,7 @@ protected: /*! - \brief A scalar function, quadratic in \a x. + \brief A scalar-valued spatial function, quadratic in \a x. */ class QuadraticXFunc : public RealFunc @@ -215,7 +252,7 @@ protected: /*! - \brief A scalar function, quadratic in \a y. + \brief A scalar-valued spatial function, quadratic in \a y. */ class QuadraticYFunc : public RealFunc @@ -235,7 +272,7 @@ protected: /*! - \brief A scalar function, quadratic in \a x. + \brief A scalar-valued spatial function, quadratic in \a z. */ class QuadraticZFunc : public RealFunc @@ -255,11 +292,11 @@ protected: /*! - \brief A scalar function, defining a linear rotation about the global Z-axis. + \brief A scalar-valued spatial function, defining a rotation about the Z-axis. \details The time component of the function argument multiplied with the function parameter \a A, is interpreted as the angle of rotation (in radians) - about the Z-axis passing through the point \a x0, \a y0. - The function then returns the translation in either \a x or \a y direction + about the global Z-axis passing through the point \a x0, \a y0. + The function then returns the translation in either X- or Y-direction (depending on the \a retX argument to the constructor) of the global point { \a X.x, \a X.y } corresponding to this rotation. \note If the function is passed a Vec3 object as argument (and not a Vec4), @@ -268,10 +305,10 @@ protected: class LinearRotZFunc : public RealFunc { - bool rX; //!< Flag telling whether to return the X (true) or Y component + bool rX; //!< Flag telling whether to return the X- (true) or Y-component real A; //!< Magnitude of the rotation - real x0; //!< Global x-coordinate of rotation centre - real y0; //!< Global y-coordinate of rotation centre + real x0; //!< Global X-coordinate of rotation centre + real y0; //!< Global Y-coordinate of rotation centre public: //! \brief Constructor initializing the function parameters. @@ -285,7 +322,7 @@ protected: /*! - \brief A scalar function, step in \a x. + \brief A scalar-valued spatial function, step in \a x. */ class StepXFunc : public RealFunc @@ -300,13 +337,13 @@ public: : fv(v), x0(X0), x1(X1) {} protected: - //! \brief Evaluates the linear function. + //! \brief Evaluates the step function. virtual real evaluate(const Vec3& X) const; }; /*! - \brief A scalar function, step in \a x and \a y. + \brief A scalar-valued spatial function, step in \a x and \a y. */ class StepXYFunc : public RealFunc @@ -325,15 +362,15 @@ public: : fv(v), x0(X0), y0(Y0), x1(X1), y1(Y1) {} protected: - //! \brief Evaluates the linear function. + //! \brief Evaluates the step function. virtual real evaluate(const Vec3& X) const; }; namespace utl { - //! \brief Creates a real function by parsing data from a character string. - const RealFunc* parseRealFunc(char* cline, real A = real(0)); + //! \brief Creates a scalar-valued function by parsing a character string. + const RealFunc* parseRealFunc(char* cline, real A = real(1)); } #endif