Added Navier plate solution to simply supported plates. Moved Poisson analytical solution to separate file.

git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1203 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
kmo 2011-09-25 13:18:12 +00:00 committed by Knut Morten Okstad
parent e52597abb6
commit af2f3d991b
4 changed files with 528 additions and 392 deletions

View File

@ -7,7 +7,7 @@
//! //!
//! \author Knut Morten Okstad / SINTEF //! \author Knut Morten Okstad / SINTEF
//! //!
//! \brief Analytic solutions for some linear elasticity and Poisson problems. //! \brief Analytic solutions for linear elasticity problems.
//! //!
//============================================================================== //==============================================================================
@ -183,217 +183,150 @@ SymmTensor CurvedBeam::evaluate (const Vec3& X) const
/*! /*!
\class NavierPlate \class NavierPlate
Navier plate solution for a constant pressure load. Navier plate solution for a constant pressure load or
a point load at an arbitrary position.
Reference: S. E. Weberg, "Todimensjonal elastisitetsteori", 1975, page 62. Reference: S. E. Weberg, "Todimensjonal elastisitetsteori", 1975, pages 61-67.
*/ */
NavierPlate::NavierPlate (double a, double b, double t, double F, NavierPlate::NavierPlate (double a, double b, double t, double E, double Poiss,
double E, double Poiss) : pz(F), nu(Poiss) double P) : pz(P), nu(Poiss), type(0), inc(2)
{ {
alpha = M_PI/a; alpha = M_PI/a;
beta = M_PI/b; beta = M_PI/b;
// Calculate and print the maximum displacement (at x=a/2, y=b/2) // Calculate and print the maximum displacement (at the centre x=a/2, y=b/2)
double x = 0.5*a;
double y = 0.5*b;
double D = E*t*t*t/(12.0-12.0*nu*nu); double D = E*t*t*t/(12.0-12.0*nu*nu);
double w = 0.0; double w = 0.0;
for (int m = 1; m < 100; m += 2) for (int m = 1; m < 100; m += inc)
for (int n = 1; n < 100; n += 2) for (int n = 1; n < 100; n += inc)
{ {
double am = alpha*m; double am = alpha*m;
double bn = beta*n; double bn = beta*n;
double am2 = am*am; double abmn = am*am + bn*bn;
double bn2 = bn*bn; w += sin(am*x)*sin(bn*y) / (double(m)*double(n)*abmn*abmn);
double pzmn = 16.0*pz / (M_PI*M_PI*m*n*(am2+bn2)*(am2+bn2));
w += (pzmn/D) * sin(0.5*am*a)*sin(0.5*bn*b);
} }
std::cout <<"NavierPlate: wmax = "<< w << std::endl; w *= 16.0*pz / (D*M_PI*M_PI);
std::streamsize oldPrec = std::cout.precision(10);
std::cout <<"\nNavierPlate: w_max = "<< w << std::endl;
std::cout.precision(oldPrec);
} }
NavierPlate::NavierPlate (double a, double b, double t, double E, double Poiss,
double P, double xi_, double eta_, double c, double d)
: pz(P), nu(Poiss), type(2), inc(1)
{
alpha = M_PI/a;
beta = M_PI/b;
xi = xi_*a;
eta = eta_*b;
c2 = 0.5*c;
d2 = 0.5*d;
if (c == 0.0 || d == 0.0) type = 1;
if (xi_ == 0.5 && eta_ == 0.5) inc = 2;
// Calculate and print the displacement at the centre x=a/2, y=b/2
double x = 0.5*a;
double y = 0.5*b;
double D = E*t*t*t/(12.0-12.0*nu*nu);
double w = 0.0;
for (int m = 1; m <= 100; m += inc)
for (int n = 1; n <= 100; n += inc)
{
double am = alpha*m;
double bn = beta*n;
double abmn = am*am + bn*bn;
double pzmn = sin(am*xi)*sin(bn*eta) / (abmn*abmn);
if (type == 2)
pzmn *= sin(am*c2)*sin(bn*d2) / (double(m)*double(n));
w += pzmn*sin(am*x)*sin(bn*y);
}
if (type == 1)
w *= 4.0*pz / (D*a*b);
else
w *= 16.0*pz / (D*M_PI*M_PI);
std::streamsize oldPrec = std::cout.precision(10);
std::cout <<"\nNavierPlate: w_centre = "<< w << std::endl;
std::cout.precision(oldPrec);
}
inline double SIN(double x) { return sin(fmod(x,2.0*M_PI)); }
inline double COS(double x) { return cos(fmod(x,2.0*M_PI)); }
void NavierPlate::addTerms (std::vector<double>& M, void NavierPlate::addTerms (std::vector<double>& M,
double x, double y, int m, int n) const double x, double y, int m, int n) const
{ {
double am = alpha*m; double am = alpha*m;
double bn = beta*n; double bn = beta*n;
double am2 = am*am; double am2 = am*am;
double bn2 = bn*bn; double bn2 = bn*bn;
double pzmn = 16.0*pz / (M_PI*M_PI*m*n*(am2+bn2)*(am2+bn2)); double dmn = double(m) * double(n);
M[0] += pzmn*(am2 + nu*bn2)*sin(am*x)*sin(bn*y);
M[1] += pzmn*(bn2 + nu*am2)*sin(am*x)*sin(bn*y); double pzmn = 0.0;
M[2] += pzmn*(nu - 1.0)*am*bn*cos(am*x)*cos(bn*y); switch (type) {
case 0: // uniform pressure
pzmn = 1.0 / dmn;
break;
case 1: // concentrated point load
pzmn = SIN(am*xi)*SIN(bn*eta);
break;
case 2: // partial load
pzmn = SIN(am*xi)*SIN(bn*eta) * SIN(am*c2)*SIN(bn*d2) / dmn;
break;
}
pzmn /= (am2+bn2)*(am2+bn2);
M[0] += pzmn*(am2 + nu*bn2)*SIN(am*x)*SIN(bn*y);
M[1] += pzmn*(bn2 + nu*am2)*SIN(am*x)*SIN(bn*y);
M[2] += pzmn*(nu - 1.0)*am*bn*COS(am*x)*COS(bn*y);
} }
SymmTensor NavierPlate::evaluate (const Vec3& X) const SymmTensor NavierPlate::evaluate (const Vec3& X) const
{ {
const int max_mn = 100; const int max_mn = 100;
const double eps = 1.0e-16; const double eps = 1.0e-8;
SymmTensor M(2); SymmTensor M(2);
double prev = 0.0; double prev = 0.0;
for (int i = 1; i <= max_mn; i += 2) for (int i = 1; i <= max_mn; i += inc)
{ {
for (int m = 1; m < i; m += 2) for (int j = 1; j < i; j += inc)
this->addTerms(M,X.x,X.y,m,i); {
for (int n = 1; n < i; n += 2) this->addTerms(M,X.x,X.y,i,j);
this->addTerms(M,X.x,X.y,i,n); this->addTerms(M,X.x,X.y,j,i);
}
this->addTerms(M,X.x,X.y,i,i); this->addTerms(M,X.x,X.y,i,i);
double norm = M.L2norm(); double norm = M.L2norm();
if (fabs(norm-prev) < eps) break; #ifdef SP_DEBUG
prev = norm; if (i == 1) std::cout <<"\nNavierPlate, X = "<< X.x <<" "<< X.y << std::endl;
std::cout << i <<": "<< M(1,1) <<" "<< M(2,2) <<" "<< M(1,2)
<<" -> "<< norm <<" "<< fabs(norm-prev)/norm << std::endl;
#endif
if (i%2)
if (fabs(norm-prev) < eps*norm)
break;
else
prev = norm;
} }
if (type == 1) // concentrated load
M *= 4.0*pz * (alpha/M_PI)*(beta/M_PI);
else // uniform or partial pressure load
M *= 16.0*pz / (M_PI*M_PI);
return M; return M;
} }
Vec3 Square2D::evaluate (const Vec3& X) const
{
double x = X.x;
double y = X.y;
double pi = M_PI;
Vec3 temp;
temp.x = -pi*sin(pi*x)*(2.0-y);
temp.y = -cos(pi*x);
return temp;
}
double Square2DHeat::evaluate (const Vec3& X) const
{
double x = X.x;
double y = X.y;
double pi = M_PI;
return -pi*pi*cos(pi*x)*(2.0-y);
}
Vec3 LshapePoisson::evaluate (const Vec3& X) const
{
double x = X.x;
double y = X.y;
double R = hypot(x,y);
double pi = M_PI;
double Rp = pow(R,-1.0/3.0);
double th = x > 0.0 ? asin(y/R) : pi-asin(y/R);
double frac = 2.0/3.0;
Vec3 temp;
temp.x = frac*Rp*sin(th/3.0);
temp.y = -frac*Rp*cos(th/3.0);
return temp;
}
Vec3 SquareSinus::evaluate (const Vec3& X) const
{
double x = X.x;
double y = X.y;
double pi = M_PI;
Vec3 temp;
temp.x = 2*pi*cos(2*pi*x)*sin(2*pi*y);
temp.y = 2*pi*sin(2*pi*x)*cos(2*pi*y);
return temp;
}
double SquareSinusSource::evaluate (const Vec3& X) const
{
double x = X.x;
double y = X.y;
double pi = M_PI;
return -2*2*2*pi*pi*sin(2*pi*x)*sin(2*pi*y);
}
Vec3 PoissonInteriorLayer::evaluate (const Vec3& X) const
{
double x = X.x;
double y = X.y;
double pi = M_PI;
Vec3 results;
double root = sqrt((x-1.25)*(x-1.25)+(y+.25)*(y+.25));
double u = SLOPE*(root-pi/3);
results.x = -SLOPE*(x-1.25) / root / (1+u*u);
results.y = -SLOPE*(y+0.25) / root / (1+u*u);
return results;
}
double PoissonInteriorLayerSol::evaluate (const Vec3& X) const
{
// picked up from http://hpfem.org/hermes2d/doc/src/benchmarks.html
double x = X.x;
double y = X.y;
double pi = M_PI;
return atan(SLOPE*(sqrt((x-1.25)*(x-1.25)+(y+.25)*(y+.25))-pi/3));
}
double PoissonInteriorLayerSource::evaluate(const Vec3& X) const
{
// nabla of the u above.
double x = X.x;
double y = X.y;
/****** maple did this for me ************/
double ans = SLOPE * pow(pow(x - 0.125e1, 0.2e1) + pow(y + 0.25e0, 0.2e1), -0.3e1 / 0.2e1) * pow(0.2e1 * x - 0.250e1, 0.2e1) / (0.1e1 + SLOPE * SLOPE * pow(sqrt(pow(x - 0.125e1, 0.2e1) + pow(y + 0.25e0, 0.2e1)) - 0.3141592654e1 / 0.3e1, 0.2e1)) / 0.4e1 - 0.2e1 * SLOPE * pow(pow(x - 0.125e1, 0.2e1) + pow(y + 0.25e0, 0.2e1), -0.1e1 / 0.2e1) / (0.1e1 + SLOPE * SLOPE * pow(sqrt(pow(x - 0.125e1, 0.2e1) + pow(y + 0.25e0, 0.2e1)) - 0.3141592654e1 / 0.3e1, 0.2e1)) + pow(SLOPE, 0.3e1) / (pow(x - 0.125e1, 0.2e1) + pow(y + 0.25e0, 0.2e1)) * pow(0.2e1 * x - 0.250e1, 0.2e1) * pow(0.1e1 + SLOPE * SLOPE * pow(sqrt(pow(x - 0.125e1, 0.2e1) + pow(y + 0.25e0, 0.2e1)) - 0.3141592654e1 / 0.3e1, 0.2e1), -0.2e1) * (sqrt(pow(x - 0.125e1, 0.2e1) + pow(y + 0.25e0, 0.2e1)) - 0.3141592654e1 / 0.3e1) / 0.2e1 + SLOPE * pow(pow(x - 0.125e1, 0.2e1) + pow(y + 0.25e0, 0.2e1), -0.3e1 / 0.2e1) * pow(0.2e1 * y + 0.50e0, 0.2e1) / (0.1e1 + SLOPE * SLOPE * pow(sqrt(pow(x - 0.125e1, 0.2e1) + pow(y + 0.25e0, 0.2e1)) - 0.3141592654e1 / 0.3e1, 0.2e1)) / 0.4e1 + pow(SLOPE, 0.3e1) / (pow(x - 0.125e1, 0.2e1) + pow(y + 0.25e0, 0.2e1)) * pow(0.2e1 * y + 0.50e0, 0.2e1) * pow(0.1e1 + SLOPE * SLOPE * pow(sqrt(pow(x - 0.125e1, 0.2e1) + pow(y + 0.25e0, 0.2e1)) - 0.3141592654e1 / 0.3e1, 0.2e1), -0.2e1) * (sqrt(pow(x - 0.125e1, 0.2e1) + pow(y + 0.25e0, 0.2e1)) - 0.3141592654e1 / 0.3e1) / 0.2e1;
/***********************************************************************/
return ans;
}
Vec3 PoissonCube::evaluate (const Vec3& X) const
{
double x = X.x;
double y = X.y;
double z = X.z;
double pi = M_PI;
Vec3 temp;
temp.x = -pi*cos(x*pi)*sin(y*pi)*sin(z*pi);
temp.y = -pi*sin(x*pi)*cos(y*pi)*sin(z*pi);
temp.z = -pi*sin(x*pi)*sin(y*pi)*cos(z*pi);
return temp;
}
double PoissonCubeSource::evaluate (const Vec3& X) const
{
double x = X.x;
double y = X.y;
double z = X.z;
double pi = M_PI;
return 3.0*pi*pi*sin(x*pi)*sin(y*pi)*sin(z*pi);
}
Vec3 PoissonLine::evaluate (const Vec3& X) const
{
double x = X.x;
double pi = M_PI;
Vec3 temp;
temp.x = -(pi/L)*cos(x*pi/L);
return temp;
}
double PoissonLineSource::evaluate (const Vec3& X) const
{
double x = X.x;
double pi = M_PI;
return (pi*pi)/(L*L)*sin(pi*x/L);
}

View File

@ -7,7 +7,7 @@
//! //!
//! \author Knut Morten Okstad / SINTEF //! \author Knut Morten Okstad / SINTEF
//! //!
//! \brief Analytic solutions for some linear elasticity and Poisson problems. //! \brief Analytic solutions for linear elasticity problems.
//! //!
//============================================================================== //==============================================================================
@ -152,9 +152,11 @@ private:
class NavierPlate : public STensorFunc class NavierPlate : public STensorFunc
{ {
public: public:
//! \brief Constructor with some default parameters. //! \brief Constructor for plate with constant pressure load.
NavierPlate(double a = 1.0, double b = 1.0, double t = 0.1, double F = 1.0, NavierPlate(double a, double b, double t, double E, double Poiss, double P);
double E = 2.1e11, double Poiss = 0.3); //! \brief Constructor for plate with partial pressure or point load.
NavierPlate(double a, double b, double t, double E, double Poiss, double P,
double xi_, double eta_, double c = 0.0, double d = 0.0);
//! \brief Empty destructor. //! \brief Empty destructor.
virtual ~NavierPlate() {} virtual ~NavierPlate() {}
@ -170,218 +172,13 @@ private:
double beta; //!< pi/(plate width) double beta; //!< pi/(plate width)
double pz; //!< Load parameter double pz; //!< Load parameter
double nu; //!< Poisson's ratio double nu; //!< Poisson's ratio
};
char type; //!< Load type parameter
/*! double xi; //!< X-position of point/partial load
\brief Analytic solution for the Poisson equation on a square. double eta; //!< Y-position of point/partial load
*/ double c2; //!< Partial load extension in X-direction
double d2; //!< Partial load extension in Y-direction
class Square2D : public VecFunc int inc; //!< Increment in Fourier term summation (1 or 2)
{
public:
//! \brief Empty constructor.
Square2D(double = 1.0) {}
//! \brief Empty destructor.
virtual ~Square2D() {}
protected:
//! \brief Evaluates the analytic flux vector at the point \a X.
virtual Vec3 evaluate(const Vec3& X) const;
};
/*!
\brief Heat source for the Poisson equation on a square.
*/
class Square2DHeat : public RealFunc
{
public:
//! \brief Empty constructor.
Square2DHeat(double = 1.0) {}
//! \brief Empty destructor.
virtual ~Square2DHeat() {}
protected:
//! \brief Evaluates the heat field at the point \a X.
virtual double evaluate(const Vec3& X) const;
};
/*!
\brief Analytic solution for the L-shape problem for the Poisson equation.
*/
class LshapePoisson : public VecFunc
{
public:
//! \brief Empty constructor.
LshapePoisson() {}
//! \brief Empty destructor.
virtual ~LshapePoisson() {}
protected:
//! \brief Evaluates the analytic flux vector at the point \a X.
virtual Vec3 evaluate(const Vec3& X) const;
};
/*!
\brief Sinus-solution of the Poisson equation on a square.
*/
class SquareSinus : public VecFunc
{
public:
//! \brief Empty constructor.
SquareSinus() {}
//! \brief Empty destructor.
virtual ~SquareSinus() {}
protected:
//! \brief Evaluates the analytic flux vector at the point \a X.
virtual Vec3 evaluate(const Vec3& X) const;
};
/*!
\brief Sinus-solution of the Poisson equation on a square.
*/
class SquareSinusSource : public RealFunc
{
public:
//! \brief Empty constructor.
SquareSinusSource() {}
//! \brief Empty destructor.
virtual ~SquareSinusSource() {}
protected:
//! \brief Evaluates the analytic flux vector at the point \a X.
virtual double evaluate(const Vec3& X) const;
};
/*!
\brief Poisson problem with smooth solution that exhibits a steep interior layer
*/
class PoissonInteriorLayer : public VecFunc
{
public:
//! \brief Empty constructor.
PoissonInteriorLayer(double s = 60.0) : SLOPE(s) {}
//! \brief Empty destructor.
virtual ~PoissonInteriorLayer() {}
protected:
//! \brief Evaluates the heat field at the point \a X.
virtual Vec3 evaluate(const Vec3& X) const;
private:
double SLOPE; //!< layer SLOPE (larger SLOPE gives problems for adaptive solvers)
};
class PoissonInteriorLayerSol : public RealFunc
{
public:
//! \brief Empty constructor.
PoissonInteriorLayerSol(double s = 60.0) : SLOPE(s) {}
//! \brief Empty destructor.
virtual ~PoissonInteriorLayerSol() {}
protected:
//! \brief Evaluates the exact solution of the heat distribution at the point \a X.
virtual double evaluate(const Vec3& X) const;
private:
double SLOPE; //!< layer SLOPE (larger SLOPE gives problems for adaptive solvers)
};
class PoissonInteriorLayerSource : public RealFunc
{
public:
//! \brief default constructor.
PoissonInteriorLayerSource(double s = 60.0) : SLOPE(s) {}
//! \brief Empty destructor.
virtual ~PoissonInteriorLayerSource() {}
protected:
virtual double evaluate(const Vec3& X) const;
private:
double SLOPE; //!< layer SLOPE (larger SLOPE gives problems for adaptive solvers)
};
/*!
\brief Analytic solution for the Poisson equation on a cube.
*/
class PoissonCube : public VecFunc
{
public:
//! \brief Empty Constructor.
PoissonCube() {}
//! \brief Empty destructor.
virtual ~PoissonCube() {}
protected:
//! \brief Evaluates the analytic flux vector at the point \a X.
virtual Vec3 evaluate(const Vec3& X) const;
};
/*!
\brief Heat source for the Poisson equation on a cube.
*/
class PoissonCubeSource : public RealFunc
{
public:
//! \brief Empty constructor.
PoissonCubeSource() {}
//! \brief Empty destructor.
virtual ~PoissonCubeSource() {}
protected:
//! \brief Evaluates the heat field at the point \a X.
virtual double evaluate(const Vec3& X) const;
};
/*!
\brief Analytic solution for the Poisson equation on a line.
*/
class PoissonLine : public VecFunc
{
public:
//! \brief Empty constructor.
PoissonLine(double r = 1.0) : L(r) {}
//! \brief Empty destructor.
virtual ~PoissonLine() {}
protected:
//! \brief Evaluates the analytic flux vector at the point \a X.
virtual Vec3 evaluate(const Vec3& X) const;
private:
double L; //!< Length parameter
};
/*!
\brief Heat source for the Poisson equation on a line.
*/
class PoissonLineSource : public RealFunc
{
public:
//! \brief Empty constructor.
PoissonLineSource(double r = 1.0) : L(r) {}
//! \brief Empty destructor.
virtual ~PoissonLineSource() {}
protected:
//! \brief Evaluates the heat field at the point \a X.
virtual double evaluate(const Vec3& X) const;
private:
double L; //!< Length parameter
}; };
#endif #endif

View File

@ -0,0 +1,164 @@
// $Id$
//==============================================================================
//!
//! \file PoissonSolutions.C
//!
//! \date Jul 1 2009
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Analytic solutions for Poisson problems.
//!
//==============================================================================
#include "PoissonSolutions.h"
#include "Vec3.h"
Vec3 Square2D::evaluate (const Vec3& X) const
{
double x = X.x;
double y = X.y;
double pi = M_PI;
Vec3 temp;
temp.x = -pi*sin(pi*x)*(2.0-y);
temp.y = -cos(pi*x);
return temp;
}
double Square2DHeat::evaluate (const Vec3& X) const
{
double x = X.x;
double y = X.y;
double pi = M_PI;
return -pi*pi*cos(pi*x)*(2.0-y);
}
Vec3 LshapePoisson::evaluate (const Vec3& X) const
{
double x = X.x;
double y = X.y;
double R = hypot(x,y);
double pi = M_PI;
double Rp = pow(R,-1.0/3.0);
double th = x > 0.0 ? asin(y/R) : pi-asin(y/R);
double frac = 2.0/3.0;
Vec3 temp;
temp.x = frac*Rp*sin(th/3.0);
temp.y = -frac*Rp*cos(th/3.0);
return temp;
}
Vec3 SquareSinus::evaluate (const Vec3& X) const
{
double x = X.x;
double y = X.y;
double pi = M_PI;
Vec3 temp;
temp.x = 2*pi*cos(2*pi*x)*sin(2*pi*y);
temp.y = 2*pi*sin(2*pi*x)*cos(2*pi*y);
return temp;
}
double SquareSinusSource::evaluate (const Vec3& X) const
{
double x = X.x;
double y = X.y;
double pi = M_PI;
return -2*2*2*pi*pi*sin(2*pi*x)*sin(2*pi*y);
}
Vec3 PoissonInteriorLayer::evaluate (const Vec3& X) const
{
double x = X.x - 1.25;
double y = X.y + 0.25;
double pi = M_PI;
Vec3 result;
double root = hypot(x,y);
double u = SLOPE*(root-pi/3.0);
result.x = -SLOPE*x / root / (1.0+u*u);
result.y = -SLOPE*y / root / (1.0+u*u);
return result;
}
/*!
\class PoissonInteriorLayerSol
Picked up from http://hpfem.org/hermes2d/doc/src/benchmarks.html
*/
double PoissonInteriorLayerSol::evaluate (const Vec3& X) const
{
double x = X.x;
double y = X.y;
double pi = M_PI;
return atan(SLOPE*(hypot(x-1.25,y+0.25)-pi/3.0));
}
double PoissonInteriorLayerSource::evaluate (const Vec3& X) const
{
double x = X.x;
double y = X.y;
/****** maple did this for me (nabla of the u above) ************/
return SLOPE * pow(pow(x - 0.125e1, 0.2e1) + pow(y + 0.25e0, 0.2e1), -0.3e1 / 0.2e1) * pow(0.2e1 * x - 0.250e1, 0.2e1) / (0.1e1 + SLOPE * SLOPE * pow(sqrt(pow(x - 0.125e1, 0.2e1) + pow(y + 0.25e0, 0.2e1)) - 0.3141592654e1 / 0.3e1, 0.2e1)) / 0.4e1 - 0.2e1 * SLOPE * pow(pow(x - 0.125e1, 0.2e1) + pow(y + 0.25e0, 0.2e1), -0.1e1 / 0.2e1) / (0.1e1 + SLOPE * SLOPE * pow(sqrt(pow(x - 0.125e1, 0.2e1) + pow(y + 0.25e0, 0.2e1)) - 0.3141592654e1 / 0.3e1, 0.2e1)) + pow(SLOPE, 0.3e1) / (pow(x - 0.125e1, 0.2e1) + pow(y + 0.25e0, 0.2e1)) * pow(0.2e1 * x - 0.250e1, 0.2e1) * pow(0.1e1 + SLOPE * SLOPE * pow(sqrt(pow(x - 0.125e1, 0.2e1) + pow(y + 0.25e0, 0.2e1)) - 0.3141592654e1 / 0.3e1, 0.2e1), -0.2e1) * (sqrt(pow(x - 0.125e1, 0.2e1) + pow(y + 0.25e0, 0.2e1)) - 0.3141592654e1 / 0.3e1) / 0.2e1 + SLOPE * pow(pow(x - 0.125e1, 0.2e1) + pow(y + 0.25e0, 0.2e1), -0.3e1 / 0.2e1) * pow(0.2e1 * y + 0.50e0, 0.2e1) / (0.1e1 + SLOPE * SLOPE * pow(sqrt(pow(x - 0.125e1, 0.2e1) + pow(y + 0.25e0, 0.2e1)) - 0.3141592654e1 / 0.3e1, 0.2e1)) / 0.4e1 + pow(SLOPE, 0.3e1) / (pow(x - 0.125e1, 0.2e1) + pow(y + 0.25e0, 0.2e1)) * pow(0.2e1 * y + 0.50e0, 0.2e1) * pow(0.1e1 + SLOPE * SLOPE * pow(sqrt(pow(x - 0.125e1, 0.2e1) + pow(y + 0.25e0, 0.2e1)) - 0.3141592654e1 / 0.3e1, 0.2e1), -0.2e1) * (sqrt(pow(x - 0.125e1, 0.2e1) + pow(y + 0.25e0, 0.2e1)) - 0.3141592654e1 / 0.3e1) / 0.2e1;
/****************************************************************************/
}
Vec3 PoissonCube::evaluate (const Vec3& X) const
{
double x = X.x;
double y = X.y;
double z = X.z;
double pi = M_PI;
Vec3 temp;
temp.x = -pi*cos(x*pi)*sin(y*pi)*sin(z*pi);
temp.y = -pi*sin(x*pi)*cos(y*pi)*sin(z*pi);
temp.z = -pi*sin(x*pi)*sin(y*pi)*cos(z*pi);
return temp;
}
double PoissonCubeSource::evaluate (const Vec3& X) const
{
double x = X.x;
double y = X.y;
double z = X.z;
double pi = M_PI;
return 3.0*pi*pi*sin(x*pi)*sin(y*pi)*sin(z*pi);
}
Vec3 PoissonLine::evaluate (const Vec3& X) const
{
double x = X.x;
double pi = M_PI;
Vec3 temp;
temp.x = -(pi/L)*cos(x*pi/L);
return temp;
}
double PoissonLineSource::evaluate (const Vec3& X) const
{
double x = X.x;
double pi = M_PI;
return (pi*pi)/(L*L)*sin(pi*x/L);
}

View File

@ -0,0 +1,242 @@
// $Id$
//==============================================================================
//!
//! \file PoissonSolutions.h
//!
//! \date Jul 1 2009
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Analytic solutions for Poisson problems.
//!
//==============================================================================
#ifndef _POISSON_SOLUTIONS_H
#define _POISSON_SOLUTIONS_H
#include "Function.h"
/*!
\brief Analytic solution for the Poisson equation on a square domain.
*/
class Square2D : public VecFunc
{
public:
//! \brief Empty constructor.
Square2D(double = 1.0) {}
//! \brief Empty destructor.
virtual ~Square2D() {}
protected:
//! \brief Evaluates the analytic flux vector at the point \a X.
virtual Vec3 evaluate(const Vec3& X) const;
};
/*!
\brief Heat source for the Poisson equation on a square domain.
*/
class Square2DHeat : public RealFunc
{
public:
//! \brief Empty constructor.
Square2DHeat(double = 1.0) {}
//! \brief Empty destructor.
virtual ~Square2DHeat() {}
protected:
//! \brief Evaluates the heat field at the point \a X.
virtual double evaluate(const Vec3& X) const;
};
/*!
\brief Analytic solution for the Poisson equation on the L-shape domain.
*/
class LshapePoisson : public VecFunc
{
public:
//! \brief Empty constructor.
LshapePoisson() {}
//! \brief Empty destructor.
virtual ~LshapePoisson() {}
protected:
//! \brief Evaluates the analytic flux vector at the point \a X.
virtual Vec3 evaluate(const Vec3& X) const;
};
/*!
\brief Sinusoidal solution of the Poisson equation on a square domain.
*/
class SquareSinus : public VecFunc
{
public:
//! \brief Empty constructor.
SquareSinus() {}
//! \brief Empty destructor.
virtual ~SquareSinus() {}
protected:
//! \brief Evaluates the analytic flux vector at the point \a X.
virtual Vec3 evaluate(const Vec3& X) const;
};
/*!
\brief Heat source for the sinusoidal solution on a square domain.
*/
class SquareSinusSource : public RealFunc
{
public:
//! \brief Empty constructor.
SquareSinusSource() {}
//! \brief Empty destructor.
virtual ~SquareSinusSource() {}
protected:
//! \brief Evaluates the heat field at the point \a X.
virtual double evaluate(const Vec3& X) const;
};
/*!
\brief Poisson problem with smooth solution and a steep interior layer.
*/
class PoissonInteriorLayer : public VecFunc
{
public:
//! \brief Default constructor.
PoissonInteriorLayer(double s = 60.0) : SLOPE(s) {}
//! \brief Empty destructor.
virtual ~PoissonInteriorLayer() {}
protected:
//! \brief Evaluates the heat flux at the point \a X.
virtual Vec3 evaluate(const Vec3& X) const;
private:
double SLOPE; //!< layer SLOPE (large value gives problems in adaptive solver)
};
class PoissonInteriorLayerSol : public RealFunc
{
public:
//! \brief Default constructor.
PoissonInteriorLayerSol(double s = 60.0) : SLOPE(s) {}
//! \brief Empty destructor.
virtual ~PoissonInteriorLayerSol() {}
protected:
//! \brief Evaluates the exact temperature distribution at the point \a X.
virtual double evaluate(const Vec3& X) const;
private:
double SLOPE; //!< layer SLOPE
};
class PoissonInteriorLayerSource : public RealFunc
{
public:
//! \brief Default constructor.
PoissonInteriorLayerSource(double s = 60.0) : SLOPE(s) {}
//! \brief Empty destructor.
virtual ~PoissonInteriorLayerSource() {}
protected:
//! \brief Evaluates the heat field at the point \a X.
virtual double evaluate(const Vec3& X) const;
private:
double SLOPE; //!< layer SLOPE
};
/*!
\brief Analytic solution for the Poisson equation on a cube domain.
*/
class PoissonCube : public VecFunc
{
public:
//! \brief Empty Constructor.
PoissonCube() {}
//! \brief Empty destructor.
virtual ~PoissonCube() {}
protected:
//! \brief Evaluates the analytic flux vector at the point \a X.
virtual Vec3 evaluate(const Vec3& X) const;
};
/*!
\brief Heat source for the Poisson equation on a cube domain.
*/
class PoissonCubeSource : public RealFunc
{
public:
//! \brief Empty constructor.
PoissonCubeSource() {}
//! \brief Empty destructor.
virtual ~PoissonCubeSource() {}
protected:
//! \brief Evaluates the heat field at the point \a X.
virtual double evaluate(const Vec3& X) const;
};
/*!
\brief Analytic solution for the Poisson equation on a line domain.
*/
class PoissonLine : public VecFunc
{
public:
//! \brief Default constructor.
PoissonLine(double r = 1.0) : L(r) {}
//! \brief Empty destructor.
virtual ~PoissonLine() {}
protected:
//! \brief Evaluates the analytic flux vector at the point \a X.
virtual Vec3 evaluate(const Vec3& X) const;
private:
double L; //!< Length parameter
};
/*!
\brief Heat source for the Poisson equation on a line.
*/
class PoissonLineSource : public RealFunc
{
public:
//! \brief Default constructor.
PoissonLineSource(double r = 1.0) : L(r) {}
//! \brief Empty destructor.
virtual ~PoissonLineSource() {}
protected:
//! \brief Evaluates the heat field at the point \a X.
virtual double evaluate(const Vec3& X) const;
private:
double L; //!< Length parameter
};
#endif