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

Added proper external energy norm calculation

Browse Source 
git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1180 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
kmo 2011-09-22 14:08:59 +00:00 committed by Knut Morten Okstad
parent 2e5e41d7cc
commit 59de729ad2
2 changed files with 31 additions and 19 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

33
src/Integrands/Poisson.C
View File

@ -82,12 +82,16 @@ void Poisson::setMode (SIM::SolutionMode mode)
}
}
double Poisson::getHeat (const Vec3& X) const
{
return heatSrc ? (*heatSrc)(X) : 0.0;
}
double Poisson::getTraction (const Vec3& X, const Vec3& n) const
{
if (tracFld)
return (*tracFld)(X)*(n);
else
return 0;
return tracFld ? (*tracFld)(X)*n : 0.0;
}
@ -132,7 +136,7 @@ bool Poisson::evalInt (LocalIntegral*& elmInt, const FiniteElement& fe,
eM->multiply(fe.dNdX,CB,false,false,true); // EK += dNdX * CB
}
// Integrate heat source if defined
// Integrate heat source, if defined
if (eS && heatSrc)
eS->add(fe.N,(*heatSrc)(X)*fe.detJxW);
@ -306,7 +310,7 @@ bool PoissonNorm::evalInt (LocalIntegral*& elmInt, const FiniteElement& fe,
// Evaluate the constitutive matrix at this point
Matrix C;
if (!problem.formCmatrix(C,X,false)) return false;
if (!problem.formCmatrix(C,X,true)) return false;
// Evaluate the finite element heat flux field
Vector sigmah;
@ -314,15 +318,18 @@ bool PoissonNorm::evalInt (LocalIntegral*& elmInt, const FiniteElement& fe,
// Integrate the energy norm a(u^h,u^h)
pnorm[0] += sigmah.dot(C*sigmah)*fe.detJxW;
// Integrate the external energy (h,u^h)
pnorm[1] += problem.getHeat(X)*problem.evalSol(fe.N)*fe.detJxW;
if (anasol)
{
// Evaluate the analytical heat flux
Vector sigma((*anasol)(X).ptr(),sigmah.size());
// Integrate the energy norm a(u,u)
pnorm[1] += sigma.dot(C*sigma)*fe.detJxW;
pnorm[2] += sigma.dot(C*sigma)*fe.detJxW;
// Integrate the error in energy norm a(u-u^h,u-u^h)
sigma -= sigmah;
pnorm[2] += sigma.dot(C*sigma)*fe.detJxW;
pnorm[3] += sigma.dot(C*sigma)*fe.detJxW;
}
return true;
@ -335,12 +342,12 @@ bool PoissonNorm::evalBou (LocalIntegral*& elmInt, const FiniteElement& fe,
Poisson& problem = static_cast<Poisson&>(myProblem);
ElmNorm& pnorm = NormBase::getElmNormBuffer(elmInt);
// Evaluate the surface flux
double T = problem.getTraction(X,normal);
// Evaluate the displacement field
// Evaluate the surface heat flux
double t = problem.getTraction(X,normal);
// Evaluate the temperature field
double u = problem.evalSol(fe.N);
// Integrate the external energy
pnorm[1] += T*u*fe.detJxW;
// Integrate the external energy (t,u^h)
pnorm[1] += t*u*fe.detJxW;
return true;
}

17
src/Integrands/Poisson.h
View File

@ -44,14 +44,16 @@ public:
//! \brief Clears the integration point traction values.
void clearTracVal() { tracVal.clear(); }
//! \brief Defines the conductivity (constitutive property).
void setMaterial(double K) { kappa = K; }
//! \brief Defines the heat source.
void setSource(RealFunc* src) { heatSrc = src; }
//! \brief Evaluates the boundary traction field (if any) at specified point.
double getTraction(const Vec3& X, const Vec3& n) const;
//! \brief Defines the conductivity (constitutive property).
void setMaterial(double K) { kappa = K; }
//! \brief Evaluates the heat source (if any) at specified point.
double getHeat(const Vec3& X) const;
//! \brief Initializes current element for numerical integration.
//! \param[in] MNPC Matrix of nodal point correspondance for current element
@ -173,6 +175,12 @@ public:
//! \brief Empty destructor.
virtual ~PoissonNorm() {}
//! \brief Returns whether this norm has explicit boundary contributions.
virtual bool hasBoundaryTerms() const { return true; }
//! \brief Returns the number of norm quantities.
virtual size_t getNoFields() const { return anasol ? 4 : 2; }
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] fe Finite element data of current integration point
@ -188,9 +196,6 @@ public:
virtual bool evalBou(LocalIntegral*& elmInt, const FiniteElement& fe,
const Vec3& X, const Vec3& normal) const;
//! \brief Returns the number of norm quantities.
virtual size_t getNoFields() const { return anasol ? 3 : 1; }
private:
VecFunc* anasol; //!< Analytical heat flux
};