OilfieldToolsNet/IFEM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Added class for calculating boundary forces.

Browse Source 
Changed stress and strain computations.

git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@817 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
rho
2011-02-23 08:45:35 +00:00
committed by Knut Morten Okstad
parent 5579c2571f
commit c15bbbd251
4 changed files with 580 additions and 147 deletions
Download Patch File Download Diff File Expand all files Collapse all files

232
src/Integrands/ChorinVelPred.C
View File

@@ -97,8 +97,19 @@ bool ChorinVelPred::initElementBou(const std::vector<int>& MNPC)
if (eS) eS->resize(nsd*nen,true);
int ierr = 0;
if (!eVs.empty() && !primsol.empty())
if (ierr = utl::gather(MNPC,nsd,primsol[0],*eVs[0]))
std::cerr <<" *** ChorinVelPred::initElementBou: Detected "
<< ierr <<" node numbers out of range."<< std::endl;
if (!ePs.empty() && !psol.empty())
if (ierr = utl::gather(MNPC,1,psol[0],*ePs[0]))
std::cerr <<" *** ChorinVelPred::initElementBou: Detected "
<< ierr <<" node numbers out of range."<< std::endl;
myMats->withLHS = false;
return true;
return ierr == 0;
}
@@ -111,10 +122,20 @@ bool ChorinVelPred::initElementBou(const std::vector<int>& MNPC1,
if (eS) eS->resize(nsd*nen1,true);
myMats->withLHS = false;
return true;
}
int ierr = 0;
if (!eVs.empty() && !primsol.empty())
if (ierr = utl::gather(MNPC1,nsd,primsol[0],*eVs[0]))
std::cerr <<" *** ChorinVelPred::initElementBou: Detected "
<< ierr <<" node numbers out of range."<< std::endl;
if (!ePs.empty() && !psol.empty())
if (ierr = utl::gather(MNPC2,1,psol[0],*ePs[0]))
std::cerr <<" *** ChorinVelPred::initElementBou: Detected "
<< ierr <<" node numbers out of range."<< std::endl;
myMats->withLHS = false;
return ierr == 0;
}
NormBase* ChorinVelPred::getNormIntegrand (AnaSol* asol) const
@@ -123,6 +144,12 @@ NormBase* ChorinVelPred::getNormIntegrand (AnaSol* asol) const
}
NormBase* ChorinVelPred::getForceIntegrand (AnaSol* asol) const
{
return new ChorinStokesForce(*const_cast<ChorinVelPred*>(this),asol);
}
bool ChorinVelPred::evalSol (Vector& s, const Vector& N,
const Matrix& dNdX, const Vec3& X,
const std::vector<int>& MNPC) const
@@ -193,7 +220,7 @@ bool ChorinVelPred::evalSol (Vector& s, const TensorFunc& asol,
}
bool ChorinVelPred::strain(const Matrix& dNdX, SymmTensor& eps) const
bool ChorinVelPred::strain(const Matrix& dNdX, Tensor& eps) const
{
int i, k, l;
@@ -203,27 +230,34 @@ bool ChorinVelPred::strain(const Matrix& dNdX, SymmTensor& eps) const
return false;
}
eps.zero();
if (!eVs.empty() && !eVs[0]->empty() && eps.dim() > 0) {
Vector& EV = *(eVs[0]);
for (i = 1;i <= dNdX.rows();i++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= k;l++)
eps(k,l) = dNdX(i,l)*EV((i-1)*nsd+k); + dNdX(i,k)*EV((i-1)*nsd+l);
for (l = 1;l <= nsd;l++)
eps(k,l) += dNdX(i,l)*EV((i-1)*nsd+k);
}
else
return false;
eps *= mu;
if (formulation == Stokes::STRESS) {
eps += eps.transpose();
eps *= 0.5;
}
return true;
}
bool ChorinVelPred::stress(const Vector& N, const Matrix& dNdX, SymmTensor& sigma) const
bool ChorinVelPred::stress(const Vector& N, const Matrix& dNdX, Tensor& sigma) const
{
int i, k;
// Compute strain
if (!this->strain(dNdX,sigma))
return false;
sigma *= mu;
// Compute pressure
if (!ePs.empty() && !ePs[0]->empty() && sigma.dim() > 0) {
@@ -231,7 +265,7 @@ bool ChorinVelPred::stress(const Vector& N, const Matrix& dNdX, SymmTensor& sigm
Vector& EP = *(ePs[0]);
for (i = 1;i <= N.size();i++)
P += EP(i)*N(i);
// Add pressure and strain contributions
for (k = 1;k <= nsd;k++)
sigma(k,k) -= P;
@@ -243,7 +277,7 @@ bool ChorinVelPred::stress(const Vector& N, const Matrix& dNdX, SymmTensor& sigm
bool ChorinVelPred::stress(const Vector& N1, const Vector& N2,
const Matrix& dN1dX, const Matrix& dN2dX,
SymmTensor& sigma) const
Tensor& sigma) const
{
int i, k;
@@ -337,7 +371,7 @@ bool ChorinStokesNorm::evalInt (LocalIntegral*& elmInt, double detJW,
Vec3 Uh; Uh = 0.0;
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
Uh[k-1] += EV((i-1)*nsd + k)*N(i);
Uh[k-1] += EV((i-1)*nsd + k)*N(i);
// L2-norm of velocity error
U -= Uh;
@@ -353,7 +387,7 @@ bool ChorinStokesNorm::evalInt (LocalIntegral*& elmInt, double detJW,
Vec3 dPh; dPh = 0.0;
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
dPh[k-1] += EP(i)*dNdX(i,k);
dPh[k-1] += EP(i)*dNdX(i,k);
// H1-seminorm of pressure
dP -= dPh;
@@ -368,44 +402,27 @@ bool ChorinStokesNorm::evalInt (LocalIntegral*& elmInt, double detJW,
// Analytical velocity gradient
Tensor gradU = (*anasol->getVectorSecSol())(X);
// Computed velocity gradient
Matrix gradUh(nsd,nsd);
gradUh.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
gradUh(k,l) += EV((i-1)*nsd + k)*dNdX(i,l);
// Numerical velocity gradient
Tensor gradUh(3);
gradUh.zero();
problem.strain(dNdX,gradUh);
if (problem.getFormulation() == Stokes::LAPLACE) {
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++) {
// Energy norm of analytical solution
pnorm[ip] += mu*gradU(k,l)*gradU(k,l)*detJW;
// Energy norm of numerical solution
pnorm[ip+1] += mu*gradUh(k,l)*gradUh(k,l)*detJW;
// Energy norm of error
value = gradU(k,l)-gradUh(k,l);
pnorm[ip+2] += mu*value*value*detJW;
}
if (problem.getFormulation() == Stokes::STRESS) {
pnorm[ip++] += 2.0*mu*gradU.innerProd(gradU)*detJW;
pnorm[ip++] += 2.0*mu*gradUh.innerProd(gradUh)*detJW;
gradUh *= -1.0;
gradU += gradUh;
pnorm[ip++] += 2.0*mu*gradU.innerProd(gradU)*detJW;
}
else {
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++) {
// Strain of analytical solution
eps = 0.5*(gradU(k,l) + gradU(l,k));
// Strain of computed solution
epsh = 0.5*(gradUh(k,l) + gradUh(l,k));
// Energy norm of analytical solution
pnorm[ip] += mu*eps*eps*detJW;
// Energy norm of computed solution
pnorm[ip+1] += mu*epsh*epsh*detJW;
// Energy norm of error
eps -= epsh;
pnorm[ip+2] += mu*epsh*epsh*detJW;
}
pnorm[ip++] += mu*gradU.innerProd(gradU)*detJW;
pnorm[ip++] += mu*gradUh.innerProd(gradUh)*detJW;
gradUh *= -1.0;
gradU += gradUh;
pnorm[ip++] += mu*gradU.innerProd(gradU)*detJW;
}
}
return true;
}
@@ -483,7 +500,7 @@ bool ChorinStokesNorm::evalInt(LocalIntegral*& elmInt,
Vec3 Uh; Uh = 0.0;
for (i = 1;i <= nbfU;i++)
for (k = 1;k <= nsd;k++)
Uh[k-1] += EV((i-1)*nsd + k)*N1(i);
Uh[k-1] += EV((i-1)*nsd + k)*N1(i);
// L2-norm of velocity error
U -= Uh;
@@ -499,7 +516,7 @@ bool ChorinStokesNorm::evalInt(LocalIntegral*& elmInt,
Vec3 dPh; dPh = 0.0;
for (i = 1;i <= nbfP;i++)
for (k = 1;k <= nsd;k++)
dPh[k-1] += EP(i)*dN2dX(i,k);
dPh[k-1] += EP(i)*dN2dX(i,k);
// H1-seminorm of pressure
dP -= dPh;
@@ -515,40 +532,22 @@ bool ChorinStokesNorm::evalInt(LocalIntegral*& elmInt,
Tensor gradU = (*anasol->getVectorSecSol())(X);
// Computed velocity gradient
Matrix gradUh(nsd,nsd);
gradUh.fill(0.0);
for (i = 1;i <= nbfU;i++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
gradUh(k,l) += EV((i-1)*nsd + k)*dN1dX(i,l);
if (problem.getFormulation() == Stokes::LAPLACE) {
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++) {
// Energy norm of analytical solution
pnorm[ip] += mu*gradU(k,l)*gradU(k,l)*detJW;
// Energy norm of numerical solution
pnorm[ip+1] += mu*gradUh(k,l)*gradUh(k,l)*detJW;
// Energy norm of error
value = gradU(k,l)-gradUh(k,l);
pnorm[ip+2] += mu*value*value*detJW;
}
Tensor gradUh(nsd);
problem.strain(dN1dX,gradUh);
if (problem.getFormulation() == Stokes::STRESS) {
pnorm[ip++] += 2.0*mu*gradU.innerProd(gradU)*detJW;
pnorm[ip++] += 2.0*mu*gradUh.innerProd(gradUh)*detJW;
gradUh *= -1.0;
gradU += gradUh;
pnorm[ip++] += 2.0*mu*gradU.innerProd(gradU)*detJW;
}
else {
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++) {
// Strain of analytical solution
eps = 0.5*(gradU(k,l) + gradU(l,k));
// Strain of computed solution
epsh = 0.5*(gradUh(k,l) + gradUh(l,k));
// Energy norm of analytical solution
pnorm[ip] += mu*eps*eps*detJW;
// Energy norm of computed solution
pnorm[ip+1] += mu*epsh*epsh*detJW;
// Energy norm of error
eps -= epsh;
pnorm[ip+2] += mu*epsh*epsh*detJW;
}
pnorm[ip++] += mu*gradU.innerProd(gradU)*detJW;
pnorm[ip++] += mu*gradUh.innerProd(gradUh)*detJW;
gradUh *= -1.0;
gradU += gradUh;
pnorm[ip++] += mu*gradU.innerProd(gradU)*detJW;
}
}
@@ -556,3 +555,76 @@ bool ChorinStokesNorm::evalInt(LocalIntegral*& elmInt,
}
bool ChorinStokesForce::initElementBou(const std::vector<int>& MNPC)
{
return problem.initElementBou(MNPC);
}
bool ChorinStokesForce::initElementBou(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2,
size_t n1)
{
return problem.initElementBou(MNPC1,MNPC2,n1);
}
bool ChorinStokesForce::evalBou(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X, const Vec3& normal) const
{
const int nsd = dNdX.cols();
ElmNorm* eNorm = dynamic_cast<ElmNorm*>(elmInt);
if (!eNorm) {
std::cerr <<" *** StokesForce::evalBou: No force array."<< std::endl;
return false;
}
ChorinVelPred* prob = static_cast<ChorinVelPred*>(&problem);
double mu = prob->getViscosity(X);
// Numerical approximation of stress
Tensor sigmah(3);
prob->stress(N,dNdX,sigmah);
// Traction
Vec3 th = sigmah*normal;
// Numerical force term
ElmNorm& pnorm = *eNorm;
int i, ip = 0;
for (i = 0;i < nsd;i++)
pnorm[ip++] += th[i]*detJW;
// Analytical force term and error norm
if ( anasol && anasol->getScalarSol() && anasol->getVectorSecSol()) {
real P = (*anasol->getScalarSol())(X);
Tensor sigma = (*anasol->getVectorSecSol())(X);
// Symmetrice for stress formulation
if (prob->getFormulation() == Stokes::STRESS)
sigma += sigma.transpose();
// Analytical stress
sigma *= mu;
for (i = 1;i <= nsd;i++)
sigma(i,i) -= P;
// Analytical traction
Vec3 t = sigma*normal;
for (i = 0;i < nsd;i++)
pnorm[ip++] += t[i]*detJW;
// Error in traction
t -= th;
pnorm[ip++] += t*t*detJW;
}
return true;
}

79
src/Integrands/ChorinVelPred.h
View File

@@ -168,6 +168,12 @@ public:
//! \param[in] asol Pointer to analytical solution field (optional)
virtual NormBase* getNormIntegrand(AnaSol* asol = 0) const;
//! \brief Returns a pointer to an Integrand for boundary force evaluation.
//! \note The Integrand is allocated dynamically and has to be deleted
//! manually when leaving the scope of the pointer returned.
//! \param[in] asol Pointer to analytical solution field (optional)
virtual NormBase* getForceIntegrand(AnaSol* asol = 0) const;
//! \brief Accesses the velocity solution vectors of current patch
Vector& getVelocity(int n = 0) { return this->getSolution(n); }
@@ -192,17 +198,16 @@ public:
//! \brief Get number of pressure solutions
size_t getNoPressures() const { return psol.size(); }
protected:
//! \brief Calculates viscous part of stress tensor at current point.
//! \param[in] dNdX Basis function gradients at current point
//! \param[out] eps Strain tensor at current point
virtual bool strain(const Matrix& dNdX, SymmTensor& eps) const;
virtual bool strain(const Matrix& dNdX, Tensor& eps) const;
//! \brief Calculates the (Cauchy) stress tensor at current point
//! \param[in] N Basis functions at current point
//! \param[in] dNdX Basis function gradients at current point
//! \param[out] sigma Strain tensor at current point
bool stress(const Vector& N, const Matrix& dNdX, SymmTensor& sigma) const;
bool stress(const Vector& N, const Matrix& dNdX, Tensor& sigma) const;
//! \brief Calculates the (Cauchy) stress tensor at current point
//! \param[in] N1 Velocity basis functions at current point
@@ -210,8 +215,9 @@ public:
//! \param[out] sigma Strain tensor at current point
bool stress(const Vector& N1, const Vector& N2,
const Matrix& dN1dX, const Matrix& dN2dX,
SymmTensor& sigma) const;
Tensor& sigma) const;
protected:
bool incPressure; //!< Incremental pressure formulation
bool mixedFEM; //!< Mixed FE formulation
@@ -273,4 +279,67 @@ class ChorinStokesNorm : public StokesNorm
};
/*!
\brief Class representing for computing boundary force
*/
class ChorinStokesForce : public StokesForce
{
public:
//! \brief The only constructor initializes its data members.
//! \param[in] p The Stokes problem to evaluate norms for
//! \param[in] a The analytical velocity and pressure fields (optional)
ChorinStokesForce(ChorinVelPred& p, AnaSol* a = 0) : StokesForce(p,a) {}
//! \brief Empty destructor.
virtual ~ChorinStokesForce() {}
//! \brief Initializes current element for boundary integration.
//! \param[in] MNPC Matrix of nodal point correspondance for current element
virtual bool initElementBou(const std::vector<int>& MNPC);
//! \brief Initializes current element for boundary integration (mixed).
//! \param[in] MNPC1 Nodal point correspondance for the basis 1
//! \param[in] MNPC2 Nodal point correspondance for the basis 2
//! \param[in] n1 Number of nodes in basis 1 on this patch
virtual bool initElementBou(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2,
size_t n1);
//! \brief Evaluates the integrand at a boundary point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] time Parameters for nonlinear and time-dependent simulations
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] normal Boundary normal vector at current integration point
//!
//! \details The default implementation forwards to the stationary version.
//! Reimplement this method for time-dependent or non-linear problems.
virtual bool evalBou(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X, const Vec3& normal) const;
//! \brief Evaluates the integrand at a boundary point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] time Parameters for nonlinear and time-dependent simulations
//! \param[in] N1 Basis function values, field 1
//! \param[in] N2 Basis function values, field 2
//! \param[in] dN1dX Basis function gradients, field 1
//! \param[in] dN2dX Basis function gradients, field 2
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] normal Boundary normal vector at current integration point
//!
//! \details This interface is used for mixed formulations.
//! The default implementation forwards to the stationary version.
//! Reimplement this method for time-dependent or non-linear problems.
virtual bool evalBou(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N1, const Vector& N2,
const Matrix& dN1dX, const Matrix& dN2dX,
const Vec3& X, const Vec3& normal) const
{ return false; }
};
#endif

288
src/Integrands/Stokes.C
View File

@@ -70,8 +70,14 @@ bool Stokes::initElementBou (const std::vector<int>& MNPC)
{
eS->resize(nf*MNPC.size(),true);
int ierr = 0;
if (ierr = utl::gather(MNPC,nf,primsol[0],*eVs[0]))
std::cerr <<" *** Stokes::initElementBou: Detected "
<< ierr <<" node numbers out of range."<< std::endl;
myMats->withLHS = false;
return true;
return ierr == 0;
}
@@ -152,6 +158,12 @@ NormBase* Stokes::getNormIntegrand (AnaSol* asol) const
}
NormBase* Stokes::getForceIntegrand (AnaSol* asol) const
{
return new StokesForce(*const_cast<Stokes*>(this),asol);
}
bool Stokes::getIntegralResult (LocalIntegral*& elmInt) const
{
elmInt = myMats;
@@ -159,15 +171,138 @@ bool Stokes::getIntegralResult (LocalIntegral*& elmInt) const
}
bool Stokes::getBodyForce(const Vec3& X, Vector& f) const
{
const Vec4* Y = dynamic_cast<const Vec4*>(&X);
if (!Y) return false;
const double PI = 3.141592653589793238462;
const double x = X[0];
const double y = X[1];
const double t = Y->t;
f.fill(0.0);
f(1) = rho*pow(sin(PI*x),2.0)*sin(2.0*PI*y)*cos(t) +
rho*PI*pow(sin(PI*x)*sin(2.0*PI*y)*sin(t),2.0)*sin(2.0*PI*x) -
2.0*rho*PI*pow(sin(PI*x)*sin(PI*y)*sin(t),2.0)*sin(2.0*PI*x)*cos(2.0*PI*y) -
PI*sin(PI*x)*sin(PI*y)*sin(t) -
2.0*pow(PI,2.0)*mu*sin(2.0*PI*y)*sin(t)*(cos(2*PI*x)-2*pow(sin(PI*x),2.0));
f(2) = -rho*sin(2.0*PI*x)*pow(sin(PI*y),2.0)*cos(t) -
2.0*rho*PI*pow(sin(PI*x)*sin(PI*y)*sin(t),2.0)*cos(2.0*PI*x)*sin(2.0*PI*y) +
rho*PI*pow(sin(2.0*PI*x)*sin(PI*y)*sin(t),2.0)*sin(2.0*PI*y) +
PI*cos(PI*x)*cos(PI*y)*sin(t) -
4*mu*pow(PI*sin(PI*y),2.0)*sin(2.0*PI*x)*sin(t) +
2.0*mu*pow(PI,2.0)*sin(2.0*PI*x)*cos(2.0*PI*y)*sin(t);
return true;
}
bool Stokes::strain(const Matrix& dNdX, Tensor& eps) const
{
int i, k, l;
const int nf = nsd+1;
if (dNdX.cols() < nsd) {
std::cerr <<" *** Stokes::strain: Invalid dimension on dNdX "
<< dNdX.rows() <<" "<< dNdX.cols() << std::endl;
return false;
}
eps.zero();
if (!eVs.empty() && !eVs[0]->empty() && eps.dim() > 0) {
Vector& EV = *(eVs[0]);
for (i = 1;i <= dNdX.rows();i++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
eps(k,l) += dNdX(i,l)*EV((i-1)*nf+k);
}
else
return false;
if (formulation == Stokes::STRESS) {
eps += eps.transpose();
eps *= 0.5;
}
return true;
}
bool Stokes::stress(const Vector& N, const Matrix& dNdX, Tensor& sigma) const
{
int i, k;
const int nf = nsd+1;
// Compute strain
if (!this->strain(dNdX,sigma))
return false;
sigma *= mu;
// Compute pressure
if (!eVs.empty() && !eVs[0]->empty() && sigma.dim() > 0) {
double P = 0.0;
Vector& EV = *(eVs[0]);
for (i = 1;i <= N.size();i++)
P += EV(i*nf)*N(i);
// Add pressure and strain contributions
for (k = 1;k <= nsd;k++)
sigma(k,k) -= P;
}
return true;
}
bool Stokes::stress(const Vector& N1, const Vector& N2,
const Matrix& dN1dX, const Matrix& dN2dX,
Tensor& sigma) const
{
int i, k;
const int nf = nsd+1;
// Compute strain
if (!this->strain(dN1dX,sigma))
return false;
// Add pressure
if (!eVs.empty() && eVs[0]) {
double P = 0.0;
Vector& EV = *eVs[0];
for (i = 1;i <= N2.size();i++)
P += EV(i*nf)*N2(i);
// Add pressure to strain
for (k = 1;k <= nsd;k++)
sigma(k,k) -= P;
}
return true;
}
bool StokesNorm::initElement (const std::vector<int>& MNPC)
{
return problem.initElement(MNPC);
}
bool StokesNorm::initElementBou (const std::vector<int>& MNPC)
{
return problem.initElementBou(MNPC);
}
bool StokesNorm::evalInt (LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const
const Vector& N, const Matrix& dNdX,
const Vec3& X) const
{
int i, k, l;
double value, eps, epsh;
@@ -221,7 +356,7 @@ bool StokesNorm::evalInt (LocalIntegral*& elmInt, double detJW,
Vec3 Uh; Uh = 0.0;
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
Uh[k-1] += EV((i-1)*nf + k)*N(i);
Uh[k-1] += EV((i-1)*nf + k)*N(i);
// L2-norm of velocity error
U -= Uh;
@@ -237,7 +372,7 @@ bool StokesNorm::evalInt (LocalIntegral*& elmInt, double detJW,
Vec3 dPh; dPh = 0.0;
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
dPh[k-1] += EV(i*nf)*dNdX(i,k);
dPh[k-1] += EV(i*nf)*dNdX(i,k);
// H1-seminorm of pressure
dP -= dPh;
@@ -251,73 +386,102 @@ bool StokesNorm::evalInt (LocalIntegral*& elmInt, double detJW,
// Analytical velocity gradient
Tensor gradU = (*anasol->getVectorSecSol())(X);
// Numerical velocity gradient
Tensor gradUh(nsd);
problem.strain(dNdX,gradUh);
// Computed velocity gradient
Matrix gradUh(nsd,nsd);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
gradUh(k,l) += EV((i-1)*nf + k)*dNdX(i,l);
if (problem.getFormulation() == Stokes::LAPLACE) {
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++) {
// Energy norm of analytical solution
pnorm[ip++] += mu*gradU(k,l)*gradU(k,l)*detJW;
// Energy norm of numerical solution
pnorm[ip++] += mu*gradUh(k,l)*gradUh(k,l)*detJW;
// Energy norm of error
value = gradU(k,l)-gradUh(k,l);
pnorm[ip++] += mu*value*value*detJW;
}
if (problem.getFormulation() == Stokes::STRESS) {
pnorm[ip++] += 2.0*mu*gradU.innerProd(gradU)*detJW;
pnorm[ip++] += 2.0*mu*gradUh.innerProd(gradUh)*detJW;
gradUh *= -1.0;
gradU += gradUh;
pnorm[ip++] += 2.0*mu*gradU.innerProd(gradU)*detJW;
}
else {
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++) {
// Strain of analytical solution
eps = 0.5*(gradU(k,l) + gradU(l,k));
// Strain of computed solution
epsh = 0.5*(gradUh(k,l) + gradUh(l,k));
// Energy norm of analytical solution
pnorm[ip++] += mu*eps*eps*detJW;
// Energy norm of computed solution
pnorm[ip++] += mu*epsh*epsh*detJW;
// Energy norm of error
eps -= epsh;
pnorm[ip++] += mu*epsh*epsh*detJW;
}
pnorm[ip++] += mu*gradU.innerProd(gradU)*detJW;
pnorm[ip++] += mu*gradUh.innerProd(gradUh)*detJW;
gradUh *= -1.0;
gradU += gradUh;
pnorm[ip++] += mu*gradU.innerProd(gradU)*detJW;
}
}
return true;
}
bool StokesForce::initElementBou(const std::vector<int>& MNPC)
{
return problem.initElementBou(MNPC);
}
bool StokesForce::evalBou(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X, const Vec3& normal) const
{
const int nsd = dNdX.cols();
ElmNorm* eNorm = dynamic_cast<ElmNorm*>(elmInt);
if (!eNorm) {
std::cerr <<" *** StokesForce::evalBou: No force array."<< std::endl;
return false;
}
double mu = problem.getViscosity(X);
// Numerical approximation of stress
Tensor sigmah(3);
problem.stress(N,dNdX,sigmah);
// Traction
Vec3 th = sigmah*normal;
// Numerical force term
ElmNorm& pnorm = *eNorm;
int i, ip = 0;
for (i = 0;i < nsd;i++)
pnorm[ip++] += th[i]*detJW;
// Analytical force term and error norm
if (anasol->getScalarSol() && anasol->getVectorSecSol()) {
real P = (*anasol->getScalarSol())(X);
Tensor sigma = (*anasol->getVectorSecSol())(X);
// Symmetrice for stress formulation
if (problem.getFormulation() == Stokes::STRESS)
sigma += sigma.transpose();
// Analytical stress
sigma *= mu;
for (i = 1;i <= nsd;i++)
sigma(i,i) -= P;
// Analytical traction
Vec3 t = sigma*normal;
for (i = 0;i < nsd;i++)
pnorm[ip++] += t[i]*detJW;
// Error in traction
t -= th;
pnorm[ip++] += t*t*detJW;
}
return true;
}
bool Stokes::getBodyForce(const Vec3& X, Vector& f) const
size_t StokesForce::getNoFields() const
{
const Vec4* Y = dynamic_cast<const Vec4*>(&X);
if (!Y) return false;
const double PI = 3.141592653589793238462;
const double x = X[0];
const double y = X[1];
const double t = Y->t;
f.fill(0.0);
size_t nsd = problem.getNoSpaceDim();
if (anasol)
return 2*nsd + 1;
else
return nsd;
}
f(1) = rho*pow(sin(PI*x),2.0)*sin(2.0*PI*y)*cos(t) +
rho*PI*pow(sin(PI*x)*sin(2.0*PI*y)*sin(t),2.0)*sin(2.0*PI*x) -
2.0*rho*PI*pow(sin(PI*x)*sin(PI*y)*sin(t),2.0)*sin(2.0*PI*x)*cos(2.0*PI*y) -
PI*sin(PI*x)*sin(PI*y)*sin(t) -
2.0*pow(PI,2.0)*mu*sin(2.0*PI*y)*sin(t)*(cos(2*PI*x)-2*pow(sin(PI*x),2.0));
f(2) = -rho*sin(2.0*PI*x)*pow(sin(PI*y),2.0)*cos(t) -
2.0*rho*PI*pow(sin(PI*x)*sin(PI*y)*sin(t),2.0)*cos(2.0*PI*x)*sin(2.0*PI*y) +
rho*PI*pow(sin(2.0*PI*x)*sin(PI*y)*sin(t),2.0)*sin(2.0*PI*y) +
PI*cos(PI*x)*cos(PI*y)*sin(t) -
4*mu*pow(PI*sin(PI*y),2.0)*sin(2.0*PI*x)*sin(t) +
2.0*mu*pow(PI,2.0)*sin(2.0*PI*x)*cos(2.0*PI*y)*sin(t);
return true;
}

128
src/Integrands/Stokes.h
View File

@@ -63,9 +63,29 @@ public:
//! \brief Initializes current element for numerical integration.
//! \param[in] MNPC Matrix of nodal point correspondance for current element
virtual bool initElement(const std::vector<int>& MNPC);
//! \brief Initializes current element for numerical integration (mixed).
//! \param[in] MNPC1 Nodal point correspondance for the basis 1
//! \param[in] MNPC2 Nodal point correspondance for the basis 2
//! \param[in] n1 Number of nodes in basis 1 on this patch
virtual bool initElement(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2, size_t n1)
{
std::cerr <<" *** Integrand::initElement not implemented."<< std::endl;
return false;
}
//! \brief Initializes current element for boundary numerical integration.
//! \param[in] MNPC Matrix of nodal point correspondance for current element
virtual bool initElementBou(const std::vector<int>& MNPC);
//! \brief Initializes current element for boundary integration (mixed).
//! \param[in] MNPC1 Nodal point correspondance for the basis 1
//! \param[in] MNPC2 Nodal point correspondance for the basis 2
//! \param[in] n1 Number of nodes in basis 1 on this patch
virtual bool initElementBou(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2, size_t n1)
{
std::cerr <<" *** Integrand::initElementBou not implemented."<< std::endl;
return false;
}
//! \brief Evaluates the integrand at a boundary point.
//! \param elmInt The local integral object to receive the contributions
@@ -106,6 +126,12 @@ public:
//! \param[in] asol Pointer to analytical solution field (optional)
virtual NormBase* getNormIntegrand(AnaSol* asol = 0) const;
//! \brief Returns a pointer to an Integrand for boundary force evaluation.
//! \note The Integrand is allocated dynamically and has to be deleted
//! manually when leaving the scope of the pointer returned.
//! \param[in] asol Pointer to analytical solution field (optional)
virtual NormBase* getForceIntegrand(AnaSol* asol = 0) const;
//! \brief Returns problem formulation type.
SIM::Formulation getFormulation() const { return formulation; }
@@ -131,6 +157,9 @@ public:
//! \brief Returns a pointer to current element solution vector.
virtual const Vector* getElementPressure() const { return 0; }
//! \brief Returns the number of space dimensions
size_t getNoSpaceDim() const { return nsd; }
//! \brief Returns the number of solution fields.
virtual size_t getNoFields() const { return nsd; }
@@ -149,6 +178,25 @@ public:
virtual const char* getSecFieldLabel(size_t i, const char* prefix = 0) const
{ return 0; }
//! \brief Calculates viscous part of stress tensor at current point.
//! \param[in] dNdX Basis function gradients at current point
//! \param[out] eps Strain tensor at current point
virtual bool strain(const Matrix& dNdX, Tensor& eps) const;
//! \brief Calculates the (Cauchy) stress tensor at current point
//! \param[in] N Basis functions at current point
//! \param[in] dNdX Basis function gradients at current point
//! \param[out] sigma Strain tensor at current point
bool stress(const Vector& N, const Matrix& dNdX, Tensor& sigma) const;
//! \brief Calculates the (Cauchy) stress tensor at current point
//! \param[in] N1 Velocity basis functions at current point
//! \param[in] dNdX Basis function gradients at current point
//! \param[out] sigma Strain tensor at current point
bool stress(const Vector& N1, const Vector& N2,
const Matrix& dN1dX, const Matrix& dN2dX,
Tensor& sigma) const;
protected:
//! \brief Utility used by the virtual \a evalInt and \a evalBou methods.
//! \param elmInt Pointer to the integrated element quantities
@@ -195,6 +243,9 @@ public:
//! \brief Initializes current element for numerical integration.
//! \param[in] MNPC Matrix of nodal point correspondance for current element
virtual bool initElement(const std::vector<int>& MNPC);
//! \brief Initializes current element for boundary integration (mixed).
//! \param[in] MNPC Matrix of nodal point correspondance for current element
virtual bool initElementBou(const std::vector<int>& MNPC);
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
@@ -218,4 +269,81 @@ protected:
};
/*!
\brief Class representing for computing boundary force
*/
class StokesForce : public NormBase
{
public:
//! \brief The only constructor initializes its data members.
//! \param[in] p The Stokes problem to evaluate norms for
//! \param[in] a The analytical velocity and pressure fields (optional)
StokesForce(Stokes& p, AnaSol* a = 0) : problem(p), anasol(a) {}
//! \brief Empty destructor.
virtual ~StokesForce() {}
//! \brief Initializes current element for boundary integration.
//! \param[in] MNPC Matrix of nodal point correspondance for current element
virtual bool initElementBou(const std::vector<int>& MNPC);
//! \brief Initializes current element for boundary integration (mixed).
//! \param[in] MNPC1 Nodal point correspondance for the basis 1
//! \param[in] MNPC2 Nodal point correspondance for the basis 2
//! \param[in] n1 Number of nodes in basis 1 on this patch
virtual bool initElementBou(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2, size_t n1)
{ return false; }
//! \brief Evaluates the integrand at a boundary point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] time Parameters for nonlinear and time-dependent simulations
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] normal Boundary normal vector at current integration point
//!
//! \details The default implementation forwards to the stationary version.
//! Reimplement this method for time-dependent or non-linear problems.
virtual bool evalBou(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X, const Vec3& normal) const;
//! \brief Evaluates the integrand at a boundary point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] time Parameters for nonlinear and time-dependent simulations
//! \param[in] N1 Basis function values, field 1
//! \param[in] N2 Basis function values, field 2
//! \param[in] dN1dX Basis function gradients, field 1
//! \param[in] dN2dX Basis function gradients, field 2
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] normal Boundary normal vector at current integration point
//!
//! \details This interface is used for mixed formulations.
//! The default implementation forwards to the stationary version.
//! Reimplement this method for time-dependent or non-linear problems.
virtual bool evalBou(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N1, const Vector& N2,
const Matrix& dN1dX, const Matrix& dN2dX,
const Vec3& X, const Vec3& normal) const
{ return false; }
//! \brief Returns the number of primary norm quantities.
virtual size_t getNoFields() const;
//! \brief Returns the number of secondary norm quantities.
virtual size_t getNoSecFields() const { return 0; }
//! \brief Only boundary contributions here
virtual bool hasBoundaryTerms() const { return true; }
protected:
Stokes& problem; //!< The problem-specific data
AnaSol* anasol; //!< Analytical solution fields
};
#endif