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

changed: move Stokes integrands to the Apps directory

Browse Source 
git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@836 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
akva
2011-02-23 11:43:33 +00:00
committed by Knut Morten Okstad
parent 437619c6e7
commit d51768ca23
26 changed files with 0 additions and 8324 deletions
Download Patch File Download Diff File Expand all files Collapse all files

148
src/Integrands/AnalyticSolutionsStokes.C
View File

@@ -1,148 +0,0 @@
// $Id$
//==============================================================================
//!
//! \file AnalyticSolutionsStokes.C
//!
//! \date Des 10 2010
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Analytic solutions for some Stokes and Navier-Stokes problems.
//!
//==============================================================================
#include "AnalyticSolutionsStokes.h"
#include "Vec3.h"
// Define pi
#ifndef M_PI
#define PI 3.141592653589793
#else
#define PI M_PI
#endif
/*!
\class Poiseuille
Channel flow.
*/
Poiseuille::Poiseuille(double P, double diam, double len, double visc)
: Pin(P), D(diam), L(len), mu(visc)
{
scalSol = new Pressure(*this);
vecSol = new Velocity(*this);
scalSecSol = new PressureGrad(*this);
vecSecSol = new VelocityGrad(*this);
}
// Pressure solution
double Poiseuille::Pressure::evaluate(const Vec3& x) const
{
return params.Pin*(1.0-x[0]/params.L);
}
// Velocity solution
Vec3 Poiseuille::Velocity::evaluate(const Vec3& x) const
{
double dPdx = -params.Pin/params.L;
double coeff = -dPdx/(2.0*params.mu);
Vec3 u;
u[0] = coeff*x[1]*(params.D-x[1]);
return u;
}
// Pressure gradient
Vec3 Poiseuille::PressureGrad::evaluate(const Vec3& x) const
{
Vec3 dPdX;
dPdX[0] = -params.Pin/params.L;
return dPdX;
}
// Velocity gradient
Tensor Poiseuille::VelocityGrad::evaluate(const Vec3& x) const
{
Tensor dUdX(3);
dUdX(1,2) = params.Pin/(2.0*params.L*params.mu)*(1.0-2.0*x[1]);
return dUdX;
}
/*!
\class TestSolution
Analytical test solution using trigonometrical functions.
*/
TestSolution::TestSolution(double dens, double visc)
: rho(dens), mu(visc)
{
scalSol = new Pressure();
vecSol = new Velocity();
scalSecSol = new PressureGrad();
vecSecSol = new VelocityGrad();
}
// Pressure
double TestSolution::Pressure::evaluate(const Vec3& x) const
{
const Vec4* X = dynamic_cast<const Vec4*>(&x);
double t = X ? X->t : 0.0;
return cos(PI*x[0])*sin(PI*x[1])*sin(t);
}
// Velocity
Vec3 TestSolution::Velocity::evaluate(const Vec3& x) const
{
const Vec4* X = dynamic_cast<const Vec4*>(&x);
double t = X ? X->t : 0.0;
Vec3 u;
u[0] = pow(sin(PI*x[0]),2.0)*sin(2.0*PI*x[1])*sin(t);
u[1] = -sin(2.0*PI*x[0])*pow(sin(PI*x[1]),2.0)*sin(t);
return u;
}
// Pressure gradient
Vec3 TestSolution::PressureGrad::evaluate(const Vec3& x) const
{
const Vec4* X = dynamic_cast<const Vec4*>(&x);
double t = X ? X->t : 0.0;
Vec3 dPdX;
dPdX[0] = -PI*sin(PI*x[0])*sin(PI*x[1])*sin(t);
dPdX[1] = PI*cos(PI*x[0])*cos(PI*x[1])*sin(t);
return dPdX;
}
// Velocity gradient
Tensor TestSolution::VelocityGrad::evaluate(const Vec3& x) const
{
const Vec4* X = dynamic_cast<const Vec4*>(&x);
double t = X ? X->t : 0.0;
Tensor dUdX(3);
dUdX(1,1) = 2.0*PI*sin(PI*x[0])*cos(PI*x[0])*sin(2.0*PI*x[1])*sin(t);
dUdX(1,2) = 2.0*PI*pow(sin(PI*x[0]),2.0)*cos(2.0*PI*x[1])*sin(t);
dUdX(2,1) = -2.0*PI*cos(2.0*PI*x[0])*pow(sin(PI*x[1]),2.0)*sin(t);
dUdX(2,2) = -2.0*PI*sin(2.0*PI*x[0])*sin(PI*x[1])*cos(PI*x[1])*sin(t);
return dUdX;
}

158
src/Integrands/AnalyticSolutionsStokes.h
View File

@@ -1,158 +0,0 @@
// $Id$
//==============================================================================
//!
//! \file AnalyticSolutionsStokes.h
//!
//! \date Des 10 2010
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Analytical solutions for some Stokes and Navier-Stokes problems.
//!
//==============================================================================
#ifndef _ANALYTIC_SOLUTIONS_STOKES_H
#define _ANALYTIC_SOLUTIONS_STOKES_H
#include "AnaSol.h"
#include "Tensor.h"
/*!
\brief Analytic solution for Poiseuille flow.
*/
class Poiseuille : public AnaSol
{
public:
//! \brief Constructor with some default parameters.
Poiseuille(double P = 1.0, double diam = 1.0, double len = 1.0,
double visc = 1.0);
//! \brief Empty destructor.
virtual ~Poiseuille() {}
private:
double D; //!< Diameter of channel
double L; //!< Length of channel
double mu; //!< Fluid viscosity
double Pin; //!< Inlet pressure value
//! \brief Analytical pressure field.
class Pressure : public RealFunc
{
const Poiseuille& params;
protected:
virtual double evaluate(const Vec3& x) const;
public:
Pressure(const Poiseuille& p) : params(p) {}
virtual ~Pressure() {}
};
//! \brief Analytical velocity field.
class Velocity : public VecFunc
{
const Poiseuille& params;
protected:
virtual Vec3 evaluate(const Vec3& x) const;
public:
Velocity(const Poiseuille& p) : params(p) {}
virtual ~Velocity() {}
};
//! \brief Analytical pressure gradient field.
class PressureGrad : public VecFunc
{
const Poiseuille& params;
protected:
virtual Vec3 evaluate(const Vec3& x) const;
public:
PressureGrad(const Poiseuille& p) : params(p) {}
virtual ~PressureGrad() {}
};
//! \brief Analytical velocity gradient field.
class VelocityGrad : public TensorFunc
{
const Poiseuille& params;
protected:
virtual Tensor evaluate(const Vec3& x) const;
public:
VelocityGrad(const Poiseuille& p) : params(p) {}
virtual ~VelocityGrad() {}
};
};
/*!
\brief Artificial solution for Navier-Stokes equations.
*/
class TestSolution : public AnaSol
{
public:
//! \brief Constructor with some default parameters.
//! \param[in] dens Fluid density
//! \param[in] visc Fluid viscosity
TestSolution(double dens = 1.0, double visc = 1.0);
//! \brief Empty destructor.
virtual ~TestSolution() {}
private:
double rho; //!< Fluid density
double mu; //!< Fluid viscosity
//! \brief Analytical pressure field.
class Pressure : public RealFunc
{
protected:
virtual double evaluate(const Vec3& x) const;
public:
Pressure() {}
virtual ~Pressure() {}
};
//! \brief Analytical velocity field.
class Velocity : public VecFunc
{
protected:
virtual Vec3 evaluate(const Vec3& x) const;
public:
Velocity() {}
virtual ~Velocity() {}
};
//! \brief Analytical pressure gradient field.
class PressureGrad : public VecFunc
{
protected:
virtual Vec3 evaluate(const Vec3& x) const;
public:
PressureGrad() {}
virtual ~PressureGrad() {}
};
//! \brief Analytical velocity gradient field.
class VelocityGrad : public TensorFunc
{
protected:
virtual Tensor evaluate(const Vec3& x) const;
public:
VelocityGrad() {}
virtual ~VelocityGrad() {}
};
};
#endif

330
src/Integrands/ChorinPressCorr.C
View File

@@ -1,330 +0,0 @@
//==============================================================================
//!
//! \file ChorinPressCorr.C
//!
//! \date Sep 30 2010
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Integrand implementations for pressure correction in Chorin's method
//!
//==============================================================================
#include "ChorinPressCorr.h"
#include "Utilities.h"
#include "Vec3Oper.h"
ChorinPressCorr::ChorinPressCorr(unsigned short int n, double B0,
bool incPress, bool mixed)
: nsd(n), Beta0(B0), incPressure(incPress), mixedFEM(mixed)
{
// Set default values for density and viscosity
rho = mu = 1.0;
// Initialize pointers to zero
tracFld = 0;
// Initialize element matrix and vector
myMats = new ElmMats();
myMats->A.resize(1);
myMats->b.resize(1);
myMats->rhsOnly = false;
eM = &myMats->A[0];
eS = &myMats->b[0];
// Current velocity and pressure solutions are needed
int nUsols = 1;
usol.resize(nUsols);
eVs.resize(nUsols);
for (int i = 0;i < nUsols;i++)
eVs[i] = new Vector();
int nPsols = 1;
primsol.resize(nPsols);
ePs.resize(nPsols);
for (int i = 0;i < nPsols;i++)
ePs[i] = new Vector();
}
ChorinPressCorr::~ChorinPressCorr()
{
int i;
// Delete element matrix/vector
if (myMats) delete myMats;
// Delete element solution vectors
if (!eVs.empty())
for (i = 0;i < eVs.size();i++)
if (eVs[i]) delete eVs[i];
if (!ePs.empty())
for (i = 0;i < ePs.size();i++)
if (ePs[i]) delete ePs[i];
}
bool ChorinPressCorr::initElement(const std::vector<int>& MNPC)
{
const size_t nen = MNPC.size();
if (eM) eM->resize(nen,nen,true);
if (eS) eS->resize(nen,true);
int ierr = 0;
if (!ePs.empty() && !primsol.empty())
for (int i = 0;i < ePs.size();i++)
if (ierr = utl::gather(MNPC,1,primsol[i],*ePs[i]))
std::cerr << "*** ChorinPressCorr::initElement: Detected "
<< ierr << " node numbers out of range." << std::endl;
if (!eVs.empty() && !usol.empty())
for (int i = 0;i < eVs.size();i++)
if (ierr = utl::gather(MNPC,nsd,usol[i],*eVs[i]))
std::cerr << "*** ChorinPressCorr::initElement: Detected "
<< ierr << " node numbers out of range." << std::endl;
myMats->withLHS = true;
return ierr == 0;
}
bool ChorinPressCorr::initElement(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2,
size_t n1)
{
// Only pressure degrees of freedom
const size_t nen2 = MNPC2.size();
if (eM) eM->resize(nen2,nen2,true);
if (eS) eS->resize(nen2,true);
// RUNAR
std::vector<int> MNPC = MNPC2;
for (int i = 0;i < MNPC.size();i++)
MNPC[i] -= n1;
// Extract element pressure vectors
int ierr = 0;
if (!ePs.empty() && !primsol.empty())
for (int i = 0;i < ePs.size();i++)
if (ierr = utl::gather(MNPC2,1,primsol[i],*ePs[i]))
std::cerr << "*** ChorinPressCorrMixed:initElement: Detected "
<< ierr << " node numbers out of range." << std::endl;
// Extract element velocity vectors
if (!eVs.empty() && !usol.empty())
for (int i = 0;i < eVs.size();i++)
if (ierr = utl::gather(MNPC1,nsd,usol[i],*eVs[i]))
std::cerr << "*** ChorinPressCorrMixed::initElement: Detected "
<< ierr << " node numbers out of range." << std::endl;
myMats->withLHS = true;
return ierr == 0;
}
bool ChorinPressCorr::initElementBou(const std::vector<int>& MNPC)
{
const size_t nen = MNPC.size();
if (eS) eS->resize(nsd*nen,true);
myMats->withLHS = false;
return true;
}
bool ChorinPressCorr::initElementBou(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2,
size_t n1)
{
// Only pressure degrees of freedom
const size_t nen2 = MNPC2.size();
if (eS) eS->resize(nen2,true);
// Extract the element level velocity vector
// Extract element pressure vectors
int ierr = 0;
if (!ePs.empty() && !primsol.empty())
if (ierr = utl::gather(MNPC2,1,primsol[0],*ePs[0]))
std::cerr << "*** ChorinPressCorrMixed:initElement: Detected "
<< ierr << " node numbers out of range." << std::endl;
// Extract element velocity vectors
if (eVs.empty() && !usol.empty())
if (ierr = utl::gather(MNPC1,nsd,usol[0],*eVs[0]))
std::cerr << "*** ChorinPressCorrMixed::initElement: Detected "
<< ierr << " node numbers out of range." << std::endl;
myMats->withLHS = false;
return ierr == 0;
}
bool ChorinPressCorr::evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const
{
int i, j, k;
double divU;
double laplace, rhoDtInv;
double B0;
// Time coefficient
rhoDtInv = 1.0/time.dt;
rhoDtInv *= rho*detJW;
// Time integration parameter
if (fabs(time.t - time.dt) < 1.0e-8)
B0 = 1.0;
else
B0 = Beta0;
// Divergence of velocity at integration point
divU = 0.0;
Vector& EV = *eVs[0];
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
divU += EV((i-1)*nsd+k)*dNdX(i,k);
// Assembling lhs terms
if (eM && !myMats->rhsOnly) {
Matrix& EM = *eM;
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
laplace = 0.0;
for (k = 1;k <= nsd;k++)
laplace += dNdX(i,k)*dNdX(j,k);
laplace *= detJW;
EM(j,i) += laplace;
}
}
// Assembling rhs terms
if (eS) {
Vector& ES = *eS;
for (i = 1;i <= N.size();i++)
ES(i) -= B0*rhoDtInv*divU*N(i);
}
return getIntegralResult(elmInt);
}
bool ChorinPressCorr::evalInt(LocalIntegral*& elmInt,
const TimeDomain& time,
double detJW,
const Vector& N1,
const Vector& N2,
const Matrix& dN1dX,
const Matrix& dN2dX,
const Vec3& X) const
{
int i, j, k;
double divU;
double laplace, rhoDtInv;
double B0;
// Time coefficient
rhoDtInv = 1.0/time.dt;
rhoDtInv *= rho*detJW;
// Time integration parameter
if (fabs(time.t - time.dt) < 1.0e-10)
B0 = 1.0;
else
B0 = Beta0;
// Divergence of velocity at integration point
divU = 0.0;
Vector& EV = *eVs[0];
for (i = 1;i <= N1.size();i++)
for(k = 1;k <= nsd;k++)
divU += EV((i-1)*nsd+k)*dN1dX(i,k);
// Assembling matrix terms
if (eM && !myMats->rhsOnly) {
Matrix& EM = *eM;
for (i = 1;i <= N2.size();i++)
for (j = 1;j <= N2.size();j++) {
laplace = 0.0;
for (k = 1;k <= nsd;k++)
laplace += dN2dX(i,k)*dN2dX(j,k);
laplace *= detJW;
EM(j,i) += laplace;
}
}
// Assembling rhs terms
if (eS) {
Vector& ES = *eS;
for (i = 1;i <= N2.size();i++)
ES(i) -= B0*rhoDtInv*divU*N2(i);
}
return getIntegralResult(elmInt);
}
bool ChorinPressCorr::evalBou(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X, const Vec3& normal) const
{
return getIntegralResult(elmInt);
}
bool ChorinPressCorr::evalSol(Vector& s,
const Vector& N, const Matrix& dNdX,
const Vec3& X, const std::vector<int>& MNPC) const
{
int i, k;
Vector P;
if (primsol.empty()) {
std::cerr << " *** ChorinPressCorr::evalSol: No primary solution."
<< std::endl;
return false;
}
// Element pressure vector
int ierr = 0;
if (ierr = utl::gather(MNPC,1,primsol[0],P)) {
std::cerr << " *** ChorinPressCorr::evalSol: Detected "
<< ierr << " node numbers out of range." << std::endl;
return false;
}
// Evaluate flux
s.resize(nsd,true);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
s(k) -= P(i)*dNdX(i,k);
return true;
}
bool ChorinPressCorr::evalSolScal(Vector& s,
const VecFunc& asol, const Vec3& X) const
{
s = Vector(asol(X).ptr(),nsd);
return true;
}
bool ChorinPressCorr::getIntegralResult (LocalIntegral*& elmInt) const
{
elmInt = myMats;
return elmInt ? true : false;
}

218
src/Integrands/ChorinPressCorr.h
View File

@@ -1,218 +0,0 @@
//==============================================================================
//!
//! \file ChorinPressCorr.h
//!
//! \date Sep 30 2010
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Integrand implementations for pressure correction in Chorin's method
//!
//==============================================================================
#ifndef _CHORIN_PRESS_CORR_H
#define _CHORIN_PRESS_CORR_H
#include "IntegrandBase.h"
#include "SIMenums.h"
#include "ElmMats.h"
#include "Vec3.h"
#include "TimeDomain.h"
#include <map>
class LocalSystem;
class VTF;
/*!
\brief Class representing the integrand of the pressure correction step
in Chorin's method.
\details This class implements the pressure correction step of Chorin's
method using NURBS based FEM with equal order elements for velocity and
pressure. The resulting system is a Poisson equation for the pressure
increment.
*/
class ChorinPressCorr : public Integrand
{
public:
//! \brief The default constructor initializes all pointers to zero.
//! \param[in] n Number of spatial dimensions
ChorinPressCorr(unsigned short int n = 3, double B0 = 1.0,
bool incPress = false, bool mixed = false);
//! \brief The destructor frees the dynamically allocated data objects.
virtual ~ChorinPressCorr();
//! \brief Defines the traction field to use in Neumann boundary conditions.
void setTraction(VecFunc* tf) { tracFld = tf; }
//! \brief Clears the integration point traction values.
void clearTracVal() { tracVal.clear(); }
//! \brief Defines fluid properties for current volume patch.
//! \param[in] Density Mass density
//! \param[in] Viscosity Dynamic viscosity
void setFluidProperties(double Density, double Viscosity)
{ rho = Density; mu = Viscosity; }
//! \brief Initializes current element for numerical integration.
//! \param[in] MNPC Matrix of nodal point correspondance for current element
bool initElement(const std::vector<int>& MNPC);
//! \brief Initializes current element for numerical integration.
//! \param[in] MNPC1 Nodal point correspondance for the velocity
//! \param[in] MNPC2 Nodal point correspondance for the pressure
//! \param[in] n1 Number of nodes in basis 1 on this patch
bool initElement(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2, size_t n1);
//! \brief Initializes current element for boundary integration.
//! \param[in] MNPC Matrix of nodal point correspondance for current element
bool initElementBou(const std::vector<int>& MNPC);
//! \brief Initializes current element for numerical boundary integration.
//! \param[in] MNPC1 Nodal point correspondance for the velocity
//! \param[in] MNPC2 Nodal point correspondance for the pressure
bool initElementBou(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2, size_t);
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const;
//! \brief Evaluates the mixed field problem integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N1 Basis function values, field velocity
//! \param[in] N2 Basis function values, field pressure
//! \param[in] dN1dX Basis function gradients, velocity
//! \param[in] dN2dX Basis function gradients, pressure
//! \param[in] X Cartesian coordinates of current integration point
bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N1, const Vector& N2,
const Matrix& dN1dX, const Matrix& dN2dX,
const Vec3& X) 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] 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
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] N2 Basis function values, pressure
//! \param[in] dN1dX Basis function gradients, pressure
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] normal Boundary normal vector at current integration point
//!
//! \details The boundary integral is the same as that of the parent class.
//! It does not depend on the pressure and volumetric change fields.
//! Thus, this call is forwarded to the single-field parent class method.
virtual bool evalBou(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector&, const Vector& N2,
const Matrix&, const Matrix& dN2dX,
const Vec3& X, const Vec3& normal) const
{
return this->evalBou(elmInt,time,detJW,N2,dN2dX,X,normal);
}
//! \brief Evaluates the secondary solution at a result point.
//! \param[out] s Array of solution field values at current point
//! \param[in] N Basis function values at current point
//! \param[in] dNdX Basis function gradients at current point
//! \param[in] X Cartesian coordinates of current point
//! \param[in] MNPC Nodal point correspondance for the basis function values
bool evalSol(Vector& s,
const Vector& N, const Matrix& dNdX,
const Vec3& X, const std::vector<int>& MNPC) const;
//! \brief Evaluates the secondary solution at a result point (mixed problem).
//! \param[out] s The solution field values at current point
//! \param[in] N2 Basis function values at current point, pressure
//! \param[in] dN2dX Basis function gradients at current point, pressure
//! \param[in] X Cartesian coordinates of current point
//! \param[in] MNPC2 Nodal point correspondance for the pressure
virtual bool evalSol(Vector& s,
const Vector&, const Vector& N2,
const Matrix&, const Matrix& dN2dX, const Vec3& X,
const std::vector<int>& MNPC2) const
{
return this->ChorinPressCorr::evalSol(s,N2,dN2dX,X,MNPC2);
}
//! \brief Evaluates the analytical solution at an integration point.
//! \param[out] s The analytical secondary solution values at current point
//! \param[in] asol The analytical solution field
//! \param[in] X Cartesian coordinates of current point
bool evalSolScal(Vector& s, const VecFunc& asol, const Vec3& X) const;
//! \brief Returns whether there are any traction values to write to VTF.
bool hasTractionValues() const { return !tracVal.empty(); }
//! \brief Returns whether an incremental pressure formulation is used.
bool incrementalPressure() const { return incPressure; }
//! \brief If a mixed FE formulation is used
bool mixedFormulation() const { return mixedFEM; }
//! \brief Accesses the velocity solution vectors of current patch
Vector& getVelocity(int n = 0) { return usol[n]; }
//! \brief Accesses the pressure solution vectors of current patch
Vector& getPressure(int n = 0) { return getSolution(n); }
//! \brief If the element matrix should be assembled
void assembleMatrix(bool newMat) { myMats->rhsOnly = !newMat; }
//! \brief Get number of velocity solutions
size_t getNoVelocities() const { return usol.size(); }
//! \brief Get number of pressure solutions
size_t getNoPressures() const { return this->getNoSolutions(); }
protected:
//! \brief Utility used by the virtual \a evalInt and \a evalBou methods.
//! \param elmInt Pointer to the integrated element quantities
bool getIntegralResult(LocalIntegral*& elmInt) const;
// Physical properties (constant)
double mu; //!< Dynamic viscosity
double rho; //!< Fluid density
unsigned short int nsd; //!< Number of space dimensions (1, 2, or 3)
bool incPressure; //!< Incremental pressure flag
bool mixedFEM; //!< Mixed FE formulation
double Beta0; //!< Time integration parameter
// Finite element quantities
Matrix* eM; //!< Pointer to element matrix
Vector* eS; //!< Pointer to element rhs vector
std::vector<Vector*> eVs; //!< Pointer to element velocity vectors
std::vector<Vector*> ePs; //!< Pointer to element pressure vectors
std::vector<Vector> usol; //!< Velocity solution vectors for patch
ElmMats* myMats; //!< Local element matrices
VecFunc* tracFld; //!< Pointer to boundary traction field
mutable std::map<Vec3,Vec3> tracVal; //!< Traction field point values
};
#endif

349
src/Integrands/ChorinVelCorr.C
View File

@@ -1,349 +0,0 @@
//==============================================================================
//!
//! \file ChorinVelCorr.C
//!
//! \date Sep 30 2010
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Integrand implementations for velocity correction in Chorin's method
//!
//==============================================================================
#include "ChorinVelCorr.h"
#include "Utilities.h"
#include "Vec3Oper.h"
ChorinVelCorr::ChorinVelCorr(unsigned short int n, double B0,
bool incPress, bool mixed)
: nsd(n), Beta0(B0), incPressure(incPress), mixedFEM(mixed)
{
// Set default values for density and viscosity
rho = mu = 1.0;
// Initialize pointers to zero
tracFld = 0;
// Initialize element matrix and vector
myMats = new ElmMats();
myMats->A.resize(1);
myMats->b.resize(1);
myMats->rhsOnly = false;
eM = &myMats->A[0];
eS = &myMats->b[0];
// Velocity, temporary velocity and pressure correction is needed
int nUsols = 1;
primsol.resize(nUsols);
eVs.resize(nUsols);
for (int i = 0;i < nUsols;i++)
eVs[i] = new Vector();
int nPsols = 1;
psol.resize(nPsols);
ePs.resize(nPsols);
for (int i = 0;i < nPsols;i++)
ePs[i] = new Vector();
}
ChorinVelCorr::~ChorinVelCorr()
{
int i;
// Delete element matrix/vector
if (myMats) delete myMats;
// Delete element solution vectors
if (!eVs.empty())
for (i = 0;i < eVs.size();i++)
if (eVs[i]) delete eVs[i];
if (!ePs.empty())
for (i = 0;i < ePs.size();i++)
if (ePs[i]) delete ePs[i];
}
bool ChorinVelCorr::initElement(const std::vector<int>& MNPC)
{
const size_t nen = MNPC.size();
if (eM) eM->resize(nsd*nen,nsd*nen,true);
if (eS) eS->resize(nsd*nen,true);
int ierr = 0;
if (!eVs.empty() && !primsol.empty())
for (int i = 0;i < eVs.size();i++)
if (ierr = utl::gather(MNPC,nsd,primsol[i],*eVs[i]))
std::cerr << "*** ChorinVelCorr::initElement: Detected "
<< ierr << " node numbers out of range." << std::endl;
if (!ePs.empty() && !psol.empty())
for (int i = 0;i < ePs.size();i++)
if (ierr = utl::gather(MNPC,1,psol[i],*ePs[i]))
std::cerr << "*** ChorinVelCorr::initElement: Detected "
<< ierr << " node numbers out of range." << std::endl;
myMats->withLHS = true;
return ierr == 0;
}
bool ChorinVelCorr::initElement(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2,
size_t n1)
{
// Only velocity degrees of freedom
const size_t nen1 = MNPC1.size();
if (eM) eM->resize(nsd*nen1,nsd*nen1,true);
if (eS) eS->resize(nsd*nen1,true);
// Extract element velocity vector
int ierr = 0;
if (!eVs.empty() && !primsol.empty())
for (int i = 0;i < eVs.size();i++)
if (ierr = utl::gather(MNPC1,nsd,primsol[i],*eVs[i]))
std::cerr << "*** ChorinVelCorrMixed::initElement: Detected "
<< ierr << " node numbers out of range." << std::endl;
if (!ePs.empty() && !psol.empty())
for (int i = 0;i < ePs.size();i++)
if (ierr = utl::gather(MNPC2,1,psol[i],*ePs[i]))
std::cerr << "*** ChorinVelCorrMixed::initElement: Detected "
<< ierr << " node numbers out of range." << std::endl;
myMats->withLHS = true;
return ierr == 0;
}
bool ChorinVelCorr::initElementBou(const std::vector<int>& MNPC)
{
const size_t nen = MNPC.size();
if (eS) eS->resize(nsd*nen,true);
myMats->withLHS = false;
return true;
}
bool ChorinVelCorr::initElementBou(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2,
size_t n1)
{
// Only velocity degrees of freedom
const size_t nen1 = MNPC1.size();
if (eS) eS->resize(nsd*nen1,true);
// Extract element velocity vector
int ierr = 0;
if (!eVs.empty() && !primsol.empty())
if (ierr = utl::gather(MNPC1,nsd,primsol[0],*eVs[0]))
std::cerr << "*** ChorinVelCorrMixed:initElementBou: Detected "
<< ierr << " node numbers out of range." << std::endl;
// Extract element pressure vector
if (!ePs.empty() && !psol.empty())
if (ierr = utl::gather(MNPC2,1,psol[0],*ePs[0]))
std::cerr << "*** ChorinVelCorr:initElementBou: Detected "
<< ierr << " node numbers out of range." << std::endl;
myMats->withLHS = false;
return ierr == 0;
}
bool ChorinVelCorr::evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const
{
int i, j, k;
double divU, P;
double B0;
double mass, rhoDtInv;
Vector U(nsd), dPdX(nsd);
// Time coefficient
rhoDtInv = 1.0/time.dt;
rhoDtInv *= rho*detJW;
// Time integration parameter
if (fabs(time.t-time.dt) < 1.0e-8)
B0 = 1.0;
else
B0 = Beta0;
// Velocity at integration point
U.fill(0.0);
Vector& EV = *eVs[0];
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
U(k) += EV((i-1)*nsd + k)*N(i);
// Pressure at integration point
Vector& EP = *ePs[0];
P = 0.0;
for (i = 1;i <= N.size();i++)
P += EP(i)*N(i);
// Pressure gradient at integration point
dPdX.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
dPdX(k) += EP(i)*dNdX(i,k);
// Assembling lhs terms
if (eM) {
Matrix& EM = *eM;
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
mass = N(i)*N(j)*B0*rhoDtInv;
for (k = 1;k <= nsd;k++)
EM((j-1)*nsd+k,(i-1)*nsd+k) += mass;
}
}
// Assembling rhs terms
if (eS) {
Vector& ES = *eS;
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
//ES((i-1)*nsd+k) += B0*rhoDtInv*U(k)*N(i) - dPdX(k)*N(i)*detJW;
ES((i-1)*nsd+k) += B0*rhoDtInv*U(k)*N(i) + P*dNdX(i,k)*detJW;
}
return getIntegralResult(elmInt);
}
bool ChorinVelCorr::evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N1, const Vector& N2,
const Matrix& dN1dX, const Matrix& dN2dX,
const Vec3& X) const
{
int i, j, k;
double divU, P;
double B0;
double mass, rhoDtInv;
Vector U(nsd), dPdX(nsd);
// Time coefficient
rhoDtInv = 1.0/time.dt;
rhoDtInv *= rho*detJW;
// Time integration parameter
if (fabs(time.t-time.dt) < 1.0e-10)
B0 = 1.0;
else
B0 = Beta0;
// Velocity at integration point
U.fill(0.0);
Vector& EV = *eVs[0];
for (i = 1;i <= N1.size();i++)
for (k = 1;k <= nsd;k++)
U(k) += EV((i-1)*nsd + k)*N1(i);
// Pressure at integration point
Vector& EP = *ePs[0];
P = 0.0;
for (i = 1;i <= N2.size();i++)
P += EP(i)*N2(i);
// Pressure gradient at integration point
dPdX.fill(0.0);
for (i = 1;i <= N2.size();i++)
for (k = 1;k <= nsd;k++)
dPdX(k) += EP(i)*dN2dX(i,k);
// Assembling lhs terms
if (eM) {
Matrix& EM = *eM;
for (i = 1;i <= N1.size();i++)
for (j = 1;j <= N1.size();j++) {
mass = N1(i)*N1(j)*B0*rhoDtInv;
for (k = 1;k <= nsd;k++)
EM((j-1)*nsd+k,(i-1)*nsd+k) += mass;
}
}
// Assembling rhs terms
if (eS) {
Vector& ES = *eS;
for (i = 1;i <= N1.size();i++)
for (k = 1;k <= nsd;k++)
//ES((i-1)*nsd+k) += B0*rhoDtInv*U(k)*N1(i) - dPdX(k)*N1(i)*detJW;
ES((i-1)*nsd+k) += B0*rhoDtInv*U(k)*N1(i) + P*dN1dX(i,k)*detJW;
}
return getIntegralResult(elmInt);
}
bool ChorinVelCorr::evalBou(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X, const Vec3& normal) const
{
if (incPressure)
return getIntegralResult(elmInt);
if (!tracFld) {
std::cerr <<" *** ChorinVelCorr::evalBou: No tractions."<< std::endl;
return false;
}
else if (!eS) {
std::cerr <<" *** ChorinVelCorr::evalBou: No load vector."<< std::endl;
return false;
}
// Evaluate the surface traction
Vec3 T = (*tracFld)(X,normal);
// Store the traction value for vizualization
if (!T.isZero()) tracVal[X] = T;
// Integrate the force vector
Vector& ES = *eS;
for (size_t a = 1; a <= N.size(); a++)
for (short int d = 1; d <= nsd; d++)
ES(nsd*(a-1)+d) -= T[d-1]*N(a)*detJW;
return getIntegralResult(elmInt);
}
bool ChorinVelCorr::evalSol(Vector& s,
const Vector& N, const Matrix& dNdX,
const Vec3& X, const std::vector<int>& MNPC) const
{
return false;
}
bool ChorinVelCorr::evalSolScal(Vector& s,
const VecFunc& asol, const Vec3& X) const
{
s = Vector(asol(X).ptr(),nsd);
return true;
}
bool ChorinVelCorr::getIntegralResult (LocalIntegral*& elmInt) const
{
elmInt = myMats;
return elmInt ? true : false;
}

226
src/Integrands/ChorinVelCorr.h
View File

@@ -1,226 +0,0 @@
//==============================================================================
//!
//! \file ChorinVelCorr.h
//!
//! \date Sep 30 2010
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Integrand implementations for pressure correction in Chorin's method
//!
//==============================================================================
#ifndef _CHORIN_VEL_CORR_H
#define _CHORIN_VEL_CORR_H
#include "IntegrandBase.h"
#include "SIMenums.h"
#include "ElmMats.h"
#include "Vec3.h"
#include "TimeDomain.h"
#include <map>
class LocalSystem;
class VTF;
/*!
\brief Class representing the integrand of the velocity correction step
in Chorin's method.
\details This class implements the velocity correction step of Chorin's
method using NURBS based FEM with equal order elements for velocity and
pressure. The resulting system is a mass equation for the velocity
increment.
*/
class ChorinVelCorr : public Integrand
{
public:
//! \brief The default constructor initializes all pointers to zero.
//! \param[in] n Number of spatial dimensions
ChorinVelCorr(unsigned short int n = 3, double B0 = 1.0,
bool incPress = false, bool mixed = false);
//! \brief The destructor frees the dynamically allocated data objects.
virtual ~ChorinVelCorr();
//! \brief Defines the traction field to use in Neumann boundary conditions.
void setTraction(TractionFunc* tf) { tracFld = tf; }
//! \brief Clears the integration point traction values.
void clearTracVal() { tracVal.clear(); }
//! \brief Defines fluid properties for current volume patch.
//! \param[in] Density Mass density
//! \param[in] Viscosity Dynamic viscosity
void setFluidProperties(double Density, double Viscosity)
{ rho = Density; mu = Viscosity; }
//! \brief Initializes current element for numerical integration.
//! \param[in] MNPC Matrix of nodal point correspondance for current element
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
bool initElement(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2, size_t n1);
//! \brief Initializes current element for boundary integration.
//! \param[in] MNPC Matrix of nodal point correspondance for current element
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
bool initElementBou(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2, size_t n1);
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const;
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] time Parameters for nonlinear and time-dependent simulations
//! \param[in] detJW Jacobian determinant times integration point weight
//! \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
//!
//! \details This interface is used for mixed formulations only.
//! The default implementation forwards to the stationary version.
//! Reimplement this method for time-dependent or non-linear problems.
bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N1, const Vector& N2,
const Matrix& dN1dX, const Matrix& dN2dX,
const Vec3& X) 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] 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
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.
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 this->evalBou(elmInt,time,detJW,N1,dN1dX,X,normal); }
//! \brief Evaluates the secondary solution at a result point.
//! \param[out] s Array of solution field values at current point
//! \param[in] N Basis function values at current point
//! \param[in] dNdX Basis function gradients at current point
//! \param[in] X Cartesian coordinates of current point
//! \param[in] MNPC Nodal point correspondance for the basis function values
bool evalSol(Vector& s, const Vector& N, const Matrix& dNdX,
const Vec3& X, const std::vector<int>& MNPC) const;
//! \brief Evaluates the secondary solution at a result point (mixed problem).
//! \param[out] s The solution field values at current point
//! \param[in] N1 Basis function values at current point, field 1
//! \param[in] N2 Basis function values at current point, field 2
//! \param[in] dN1dX Basis function gradients at current point, field 1
//! \param[in] dN2dX Basis function gradients at current point, field 2
//! \param[in] X Cartesian coordinates of current point
//! \param[in] MNPC1 Nodal point correspondance for the basis 1
//! \param[in] MNPC2 Nodal point correspondance for the basis 2
bool evalSol(Vector& s,
const Vector& N1, const Vector& N2,
const Matrix& dN1dX, const Matrix& dN2dX,
const Vec3& X,
const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2) const
{ return this->evalSol(s,N1,dN1dX,X,MNPC1); }
//! \brief Evaluates the analytical solution at an integration point.
//! \param[out] s The analytical secondary solution values at current point
//! \param[in] asol The analytical solution field
//! \param[in] X Cartesian coordinates of current point
virtual bool evalSolScal(Vector& s, const VecFunc& asol,
const Vec3& X) const;
//! \brief Returns whether there are any traction values to write to VTF.
bool hasTractionValues() const { return !tracVal.empty(); }
//! \brief Mixed FE formulation
bool mixedFormulation() const { return mixedFEM; }
//! \brief Accesses the velocity solution vectors of current patch
Vector& getVelocity(int n = 0) { return this->getSolution(n); }
//! \brief Accesses the pressure solution vectors of the current patch
Vector& getPressure(int n = 0) { return psol[n]; }
//! \brief If the element matrix should be assembled
void assembleMatrix(bool newMat) { myMats->rhsOnly = !newMat; }
//! \brief Get number of velocity solutions
size_t getNoVelocities() const { return this->getNoSolutions(); }
//! \brief Get number of pressure solutions
size_t getNoPressures() const { return psol.size(); }
protected:
//! \brief Utility used by the virtual \a evalInt and \a evalBou methods.
//! \param elmInt Pointer to the integrated element quantities
bool getIntegralResult(LocalIntegral*& elmInt) const;
// Physical properties (constant)
double mu; //!< Dynamic viscosity
double rho; //!< Fluid density
unsigned short int nsd; //!< Number of space dimensions (1, 2, or 3)
bool incPressure; //!< Incremental pressure formulation
bool mixedFEM; //!< Mixed FE formulation
double Beta0; //!< Time integration parameter
// Finite element quantities
Matrix* eM; //!< Pointer to element matrix
Vector* eS; //!< Pointer to element rhs vector
std::vector<Vector*> eVs; //!< Pointer to element velocity vectors
std::vector<Vector*> ePs; //!< Pointer to element pressure vectors
std::vector<Vector> psol; //!< Pressure solution vectors for patch
ElmMats* myMats; //!< Local element matrices
TractionFunc* tracFld; //!< Pointer to boundary traction field
mutable std::map<Vec3,Vec3> tracVal; //!< Traction field point values
};
#endif

630
src/Integrands/ChorinVelPred.C
View File

@@ -1,630 +0,0 @@
//$Id: ChorinVelPred.C,v 1.3 2011-02-08 12:32:28 rho Exp $
//==============================================================================
//!
//! \file ChorinVelPred.C
//!
//! \date Sep 29 2010
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Integrand implementations for velocity prediction in Chorin's method
//!
//==============================================================================
#include "ChorinVelPred.h"
#include "ElmNorm.h"
#include "Tensor.h"
#include "Vec3Oper.h"
#include "Utilities.h"
ChorinVelPred::ChorinVelPred(short int n, SIM::Formulation form,
int itg, bool incPress, bool mixed)
: Stokes(n,(SIM::Formulation)form,itg), incPressure(incPress), mixedFEM(mixed)
{
// Number of fields equals number of space dimensions
nf = nsd;
}
ChorinVelPred::~ChorinVelPred()
{
// Free memory for pressure element solution vectors
if (!ePs.empty())
for (int n = 0;n < ePs.size();n++)
if (ePs[n]) delete ePs[n];
}
bool ChorinVelPred::initElement(const std::vector<int>& MNPC)
{
const size_t nen = MNPC.size();
if (eM) eM->resize(nsd*nen,nsd*nen,true);
if (eS) eS->resize(nsd*nen,true);
int ierr = 0;
if (!eVs.empty() && !primsol.empty())
for (int i = 0;i < primsol.size();i++)
if (ierr = utl::gather(MNPC,nsd,primsol[i],*eVs[i]))
std::cerr <<" *** ChorinVelPred::initElement: Detected "
<< ierr <<" node numbers out of range."<< std::endl;
if (!ePs.empty() && !psol.empty())
for (int i = 0;i < psol.size();i++)
if (ierr = utl::gather(MNPC,1,psol[i],*ePs[i]))
std::cerr <<" *** ChorinVelPred::initElement: Detected "
<< ierr <<" node numbers out of range."<< std::endl;
myMats->withLHS = true;
return ierr == 0;
}
bool ChorinVelPred::initElement(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2,
size_t n1)
{
// Only velocity degrees of freedom
const size_t nen1 = MNPC1.size();
if (eM) eM->resize(nsd*nen1,nsd*nen1,true);
if (eS) eS->resize(nsd*nen1,true);
// Extract element velocity vectors
int ierr = 0;
if (!eVs.empty() && !primsol.empty())
for (int i = 0;i < primsol.size();i++)
if (ierr = utl::gather(MNPC1,nsd,primsol[i],*eVs[i]))
std::cerr <<" *** ChorinVelPredMixed::initElement: Detected "
<< ierr <<" node numbers out of range."<< std::endl;
// Extract element pressure vectors
if (!ePs.empty() && !psol.empty())
for (int i = 0;i < psol.size();i++)
if (ierr = utl::gather(MNPC2,1,psol[i],*ePs[i]))
std::cerr <<" *** ChorinVelPredMixed::initElement: Detected "
<< ierr <<" node numbers out of range."<< std::endl;
myMats->withLHS = true;
return ierr == 0;
}
bool ChorinVelPred::initElementBou(const std::vector<int>& MNPC)
{
const size_t nen = MNPC.size();
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 ierr == 0;
}
bool ChorinVelPred::initElementBou(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2,
size_t)
{
// Only velocity degrees of freedom
const size_t nen1 = MNPC1.size();
if (eS) eS->resize(nsd*nen1,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
{
return new ChorinStokesNorm(*const_cast<ChorinVelPred*>(this),asol);
}
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
{
int ierr = 0;
if (!eVs.empty() && eVs[0] && !primsol.front().empty())
if (ierr = utl::gather(MNPC,nsd,primsol.front(),*(eVs[0]))) {
std::cerr << " *** ChorinVelPred::evalSol: Detected "
<< ierr << " node numbers out of range." << std::endl;
return false;
}
if (!ePs.empty() && ePs[0] && !psol.front().empty())
if (ierr = utl::gather(MNPC,1,psol.front(),*(ePs[0]))) {
std::cerr << " *** ChorinVelPred::evalSol: Detected "
<< ierr << " node numbers out of range." << std::endl;
return false;
}
// Evaluate stress tensor, sigma
SymmTensor sigma(nsd);
if (!this->stress(N,dNdX,sigma))
return false;
s = sigma;
return true;
}
bool ChorinVelPred::evalSol(Vector& s,
const Vector& N1, const Vector& N2,
const Matrix& dN1dX, const Matrix& dN2dX,
const Vec3& X,
const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2) const
{
// Extract element velocity vector
int ierr = 0;
if (!eVs.empty() && eVs[0] && !primsol.front().empty())
if (ierr = utl::gather(MNPC1,nsd,primsol.front(),*(eVs[0]))) {
std::cerr << " *** ChorinVelPredMixed::evalSol: Detected "
<< ierr << " node numbers out of range." << std::endl;
return false;
}
// Extract element pressure vector
if (!ePs.empty() && ePs[0] && !psol.front().empty())
if (ierr = utl::gather(MNPC2,1,psol.front(),*(ePs[0]))) {
std::cerr << " *** ChorinVelPredMixed::evalSol: Detected "
<< ierr << " node numbers out of range." << std::endl;
return false;
}
// Evaluate stress tensor, sigma
SymmTensor sigma(nsd);
if (!this->stress(N1,N2,dN1dX,dN2dX,sigma))
return false;
s = sigma;
return true;
}
bool ChorinVelPred::evalSol (Vector& s, const TensorFunc& asol,
const Vec3& X) const
{
s = asol(X);
return true;
}
bool ChorinVelPred::strain(const Matrix& dNdX, Tensor& eps) const
{
int i, k, l;
if (dNdX.cols() < nsd) {
std::cerr <<" *** ChorinVelPred::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)*nsd+k);
}
else
return false;
if (formulation == SIM::STRESS) {
eps += eps.transpose();
eps *= 0.5;
}
return true;
}
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) {
double P = 0.0;
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;
}
return true;
}
bool ChorinVelPred::stress(const Vector& N1, const Vector& N2,
const Matrix& dN1dX, const Matrix& dN2dX,
Tensor& sigma) const
{
int i, k;
// Compute strain
if (!this->strain(dN1dX,sigma))
return false;
// Add pressure
if (!ePs.empty() && ePs[0]) {
double P = 0.0;
Vector& EP = *ePs[0];
for (i = 1;i <= N2.size();i++)
P += EP(i)*N2(i);
// Add pressure to strain
for (k = 1;k <= nsd;k++)
sigma(k,k) -= P;
}
return true;
}
bool ChorinStokesNorm::evalInt (LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const
{
int i, k, l;
double value, eps, epsh;
const int nsd = dNdX.cols();
const int nf = nsd+1;
const int nbf = N.size();
if (!anasol->hasScalarSol() && !anasol->hasVectorSol()) {
std::cerr <<" *** ChorinStokesNorm::evalInt: No analytical solution."<< std::endl;
return false;
}
ElmNorm* eNorm = dynamic_cast<ElmNorm*>(elmInt);
if (!eNorm) {
std::cerr <<" *** ChorinStokesNorm::evalInt: No norm array."<< std::endl;
return false;
}
// Find current element velocity vector
Vector* eV = const_cast<Vector*>(problem.getElementVelocity());
if (!eV) {
std::cerr <<" *** ChorinStokesNorm::evalInt: No velocity" << std::endl;
return false;
}
Vector& EV = *eV;
// Find current element pressure vector
Vector* eP = const_cast<Vector*>(problem.getElementPressure());
if (!eP) {
std::cerr <<" *** ChorinStokesNorm::evalInt: No pressure" << std::endl;
return false;
}
Vector& EP = *eP;
ElmNorm& pnorm = *eNorm;
int ip = 0;
// Pressure norms
if (anasol->getScalarSol()) {
// Analytical pressure
real P = (*anasol->getScalarSol())(X);
// Computed pressure
real Ph = 0.0;
for (i = 1;i <= nbf;i++)
Ph += EP(i)*N(i);
// Integral of computed pressure
pnorm[ip++] += Ph*detJW;
// Integral of basis functions
pnorm[ip++] += detJW;
// L2-norm of pressure error
P -= Ph;
pnorm[ip++] += P*P*detJW;
}
// Velocity norms
if (anasol->getVectorSol()) {
// Analytical velocity
Vec3 U = (*anasol->getVectorSol())(X);
// Computed velocity
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);
// L2-norm of velocity error
U -= Uh;
pnorm[ip++] += U*U*detJW;
}
// Pressure gradient norm
if (anasol->getScalarSecSol()) {
// Analytical pressure gradient
Vec3 dP = (*anasol->getScalarSecSol())(X);
// Computed pressure gradient
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);
// H1-seminorm of pressure
dP -= dPh;
pnorm[ip++] += dP*dP*detJW;
}
// Energy norms
if (anasol->getVectorSecSol()) {
// Viscosity
const double mu = problem.getViscosity(X);
// Analytical velocity gradient
Tensor gradU = (*anasol->getVectorSecSol())(X);
// Numerical velocity gradient
Tensor gradUh(3);
gradUh.zero();
problem.strain(dNdX,gradUh);
if (problem.getFormulation() == SIM::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 {
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 ChorinStokesNorm::evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N1, const Vector& N2,
const Matrix& dN1dX, const Matrix& dN2dX,
const Vec3& X) const
{
int i, k, l;
double value, eps, epsh;
const int nsd = dN1dX.cols();
const int nf = nsd+1;
const int nbfU = N1.size();
const int nbfP = N2.size();
if (!anasol->hasScalarSol() && !anasol->hasVectorSol()) {
std::cerr <<" *** ChorinStokesNorm::evalInt: No analytical solution."<< std::endl;
return false;
}
ElmNorm* eNorm = dynamic_cast<ElmNorm*>(elmInt);
if (!eNorm) {
std::cerr <<" *** ChorinStokesNorm::evalInt: No norm array."<< std::endl;
return false;
}
// Find current element velocity vector
Vector* eV = const_cast<Vector*>(problem.getElementVelocity());
if (!eV) {
std::cerr <<" *** ChorinStokesNorm::evalInt: No svelocity" << std::endl;
return false;
}
Vector& EV = *eV;
// Find current element pressure vector
Vector* eP = const_cast<Vector*>(problem.getElementPressure());
if (!eP) {
std::cerr <<" *** ChorinStokesNorm::evalInt: No pressure" << std::endl;
return false;
}
Vector& EP = *eP;
ElmNorm& pnorm = *eNorm;
int ip = 0;
// Pressure norms
if (anasol->getScalarSol()) {
// Analytical pressure
real P = (*anasol->getScalarSol())(X);
// Computed pressure
real Ph = 0.0;
for (i = 1;i <= nbfP;i++)
Ph += EP(i)*N2(i);
// Integral of computed pressure
pnorm[ip++] += Ph*detJW;
// Integral of basis functions
pnorm[ip++] += detJW;
// L2-norm of pressure error
P -= Ph;
pnorm[ip++] += P*P*detJW;
}
// Velocity norms
if (anasol->getVectorSol()) {
// Analytical velocity
Vec3 U = (*anasol->getVectorSol())(X);
// Computed velocity
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);
// L2-norm of velocity error
U -= Uh;
pnorm[ip++] += U*U*detJW;
}
// Pressure gradient norm
if (anasol->getScalarSecSol()) {
// Analytical pressure gradient
Vec3 dP = (*anasol->getScalarSecSol())(X);
// Computed pressure gradient
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);
// H1-seminorm of pressure
dP -= dPh;
pnorm[ip++] += dP*dP*detJW;
}
// Energy norms
if (anasol->getVectorSecSol()) {
// Viscosity
const double mu = problem.getViscosity(X);
// Analytical velocity gradient
Tensor gradU = (*anasol->getVectorSecSol())(X);
// Computed velocity gradient
Tensor gradUh(nsd);
problem.strain(dN1dX,gradUh);
if (problem.getFormulation() == SIM::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 {
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 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() == SIM::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;
}

346
src/Integrands/ChorinVelPred.h
View File

@@ -1,346 +0,0 @@
//==============================================================================
//!
//! \file ChorinVelPred.h
//!
//! \date Sep 29 2010
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Integrand implementation for velocity prediction in Chorin's method
//!
//==============================================================================
#ifndef _CHORIN_VEL_PRED_H
#define _CHORIN_VEL_PRED_H
#include "Stokes.h"
#include "ElmMats.h"
#include "Vec3.h"
#include "TimeDomain.h"
#include "SIMenums.h"
#include <map>
class LocalSystem;
class VTF;
/*!
\brief Class representing the integrand of the velocity prediction step
in Chorin's method.
\details This class implements the velocity prediction step of Chorin's
method using NURBS based FEM with equal order elements for velocity and
pressure.
*/
class ChorinVelPred : public Stokes
{
protected:
//! \brief The default constructor initializes all pointers to zero.
//! \param[in] n Number of spatial dimensions
//! \param[in] form The solution formulation to use (Laplace/Stress)
//! \param[in] itg The integrandtype to use
//! \param[in] incPress \e true if incremental pressure formulation
ChorinVelPred(short int n, SIM::Formulation form,
int itg, bool incPress, bool mixed = false);
public:
//! \brief The destructor frees dynamically allocated objects.
virtual ~ChorinVelPred();
//! \brief Initializes current element for numerical integration.
//! \param[in] MNPC Matrix of nodal point correspondance for current element
bool initElement(const std::vector<int>& MNPC);
//! \brief Initializes current element for numerical integration.
//! \param[in] MNPC1 Nodal point correspondance for the velocity
//! \param[in] MNPC2 Nodal point correspondance for the pressure
//! \param[in] n1 Number of nodes in velocity basis on this patch
bool initElement(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2, size_t n1);
//! \brief Initializes current element for boundary integration.
//! \param[in] MNPC Matrix of nodal point correspondance for current element
bool initElementBou(const std::vector<int>& MNPC);
//! \brief Initializes current element for numerical boundary integration.
//! \param[in] MNPC1 Nodal point correspondance for the velocity
//! \param[in] MNPC2 Nodal point correspondance for the pressure
bool initElementBou(const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2, size_t);
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] time Parameters for nonlinear and time-dependent simulations
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
//!
//! \details This interface is used when \a getIntegrandType returns 1.
virtual bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const
{ return false; }
//! \brief Evaluates the integrand at an interior 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] d2NdX2 Basis function second derivatives
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] h Characteristic element length
//!
//! \details This interface is used when \a getIntegrandType returns 2.
virtual bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Matrix3D& d2NdX2, const Vec3& X,
double h = 0.0) const
{ 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] 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
virtual bool evalBou(LocalIntegral*& elmInt, const TimeDomain& time,
double detJW, const Vector& N, const Matrix& dNdX,
const Vec3& X, const Vec3& normal) const
{ 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] N1 Basis function values, velocity
//! \param[in] dN1dX Basis function gradients, velocity
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] normal Boundary normal vector at current integration point
//!
//! \details The boundary integral is the same as that of the parent class.
//! It does not depend on the pressure and volumetric change fields.
//! Thus, this call is forwarded to the single-field parent class method.
bool evalBou(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N1, const Vector&,
const Matrix& dN1dX, const Matrix&,
const Vec3& X, const Vec3& normal) const
{
return this->evalBou(elmInt,time,detJW,N1,dN1dX,X,normal);
}
//! \brief Evaluates the secondary solution at current integration point.
//! \param[out] s The solution field values
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] MNPC Matrix of nodal point correspondance
virtual bool evalSol(Vector& s,
const Vector& N, const Matrix& dNdX,
const Vec3& X, const std::vector<int>& MNPC) const;
//! \brief Evaluates the secondary solution at a result point (mixed problem).
//! \param[out] s The solution field values at current point
//! \param[in] N1 Basis function values at current point, velocity
//! \param[in] N2 Basis function values at current point, pressure
//! \param[in] dN1dX Basis function gradients at current point, velocity
//! \param[in] dN2dX Basis function gradients at current point, pressure
//! \param[in] X Cartesian coordinates of current point
//! \param[in] MNPC1 Nodal point correspondance for the velocity
//! \param[in] MNPC2 Nodal point correspondance for the pressure
bool evalSol(Vector& s,
const Vector& N1, const Vector& N2,
const Matrix& dN1dX, const Matrix& dN2dX, const Vec3& X,
const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2) const;
//! \brief Evaluates the analytical secondary solution at the given point.
//! \param[out] s The solution field values
//! \param[in] asol The analytical solution field
//! \param[in] X Cartesian coordinates of current integration point
virtual bool evalSol(Vector& s, const TensorFunc& asol, const Vec3& X) const;
//! \brief Returns a pointer to an Integrand for solution norm 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* 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); }
//! \brief Accesses the pressure solution vectors of the current patch
Vector& getPressure(int n = 0) { return psol[n]; }
//! \brief Returns a pointer to current element solution vector.
const Vector* getElementVelocity() const { return eVs[0]; }
//! \brief Returns a pointer to current element solution vector.
const Vector* getElementPressure() const { return ePs[0]; }
//! \brief If an incremental pressure formulation is used
bool incrementalPressure () const { return incPressure; }
//! \brief If a mixed FE formulation is used
bool mixedFormulation() const { return mixedFEM; }
//! \brief Get number of velocity solutions
size_t getNoVelocities() const { return this->getNoSolutions(); }
//! \brief Get number of pressure solutions
size_t getNoPressures() const { return psol.size(); }
//! \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:
bool incPressure; //!< Incremental pressure formulation
bool mixedFEM; //!< Mixed FE formulation
// Finite element quantities
std::vector<Vector*> ePs; //!< Pointers to element pressure vectors
std::vector<Vector> psol; //!< Pressure solution vectors for this patch
};
/*!
\brief Class representing the integrand of Stokes norms for Chorin method.
*/
class ChorinStokesNorm : public StokesNorm
{
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)
ChorinStokesNorm(Stokes& p, AnaSol* a) : StokesNorm(p,a) {}
//! \brief Empty destructor.
virtual ~ChorinStokesNorm() {}
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const;
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] time Parameters for nonlinear and time-dependent simulations
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N1 Basis function values, velocity
//! \param[in] N2 Basis function values, pressure
//! \param[in] dN1dX Basis function gradients, velocity
//! \param[in] dN2dX Basis function gradients, pressure
//! \param[in] X Cartesian coordinates of current integration point
//!
//! \details This interface is used for mixed formulations only.
//! The default implementation forwards to the stationary version.
//! Reimplement this method for time-dependent or non-linear problems.
virtual bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N1, const Vector& N2,
const Matrix& dN1dX, const Matrix& dN2dX,
const Vec3& X) const;
//! \brief Returns the number of primary norm quantities.
virtual size_t getNoFields() const { return anasol ? 4 : 2; }
//! \brief Returns the number of secondary norm quantities.
virtual size_t getNoSecFields() const { return anasol ? 4 : 1; }
};
/*!
\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

568
src/Integrands/ChorinVelPredBDF2.C
View File

@@ -1,568 +0,0 @@
//==============================================================================
//!
//! \file ChorinVelPredBDF2.C
//!
//! \date Nov 23 2010
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief BDF2 implementations for velocity prediction in Chorin's method
//!
//==============================================================================
#include "ChorinVelPredBDF2.h"
#include "Tensor.h"
#include "Vec3Oper.h"
#include "Utilities.h"
ChorinVelPredBDF2::ChorinVelPredBDF2(short int n, SIM::Formulation form,
int itg, bool mixed)
: ChorinVelPredBE(n,form,itg,true,mixed)
{
int i;
// Resize velocity element solution vectors
int nUsols = 3;
eVs.resize(nUsols);
for (i = 0;i < nUsols;i++)
eVs[i] = new Vector();
primsol.resize(nUsols);
// Resize pressure element solution vectors
int nPsols = 3;
ePs.resize(nPsols);
for (i = 0;i < nPsols;i++)
ePs[i] = new Vector();
psol.resize(nPsols);
}
ChorinVelPredBDF2::~ChorinVelPredBDF2() {}
bool ChorinVelPredBDF2::evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const
{
int i, j;
int k, l;
int dof;
real dtInv, temp, mass;
real P, conv, diff, laplace;
real B0, B1, B2;
Vector U(nsd), UP(nsd), UPP(nsd), dPdX(nsd);
// Inverse of timestep size
if (time.dt > 0.0)
dtInv = 1.0/time.dt;
else
dtInv = 0.0;
// Element velocity vectors
Vector& EVP = *eVs[1];
Vector& EVPP = *eVs[2];
// Element pressure vector
Vector& EP = *ePs[1];
Vector& EPP = *ePs[2];
if (fabs(time.t-time.dt) < 1.0e-8) {
EVPP = EVP;
EPP = EP;
B0 = 1.0;
B1 = -1.0;
B2 = 0.0;
}
else {
B0 = 1.5;
B1 = -2.0;
B2 = 0.5;
}
// Velocity in integration point at previous timesteps
UP.fill(0.0); UPP.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++) {
dof = (i-1)*nsd + k;
UP(k) += EVP(dof)*N(i);
UPP(k) += EVPP(dof)*N(i);
}
// Convective velocity in integration point
U = UP;
U *= 2.0;
U -= UPP;
// Pressure in integration point
P = 0.0;
for (i = 1;i <= N.size();i++)
P += EP(i)*N(i);
// Pressure gradient in integration point
// dPdX.fill(0.0);
// for (i = 1;i <= N.size();i++)
// for (k = 1;k <= nsd;k++)
// dPdX(k) += (2.0*EP(i) - EPP(i))*dNdX(i,k);
// Loop over basis functions
// Lefthand size terms
if (eM) {
Matrix& EM = *eM;
// Time derivative terms
temp = rho*dtInv*detJW;
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
mass = B0*temp*N(i)*N(j);
for (k = 1;k <= nsd;k++)
EM((j-1)*nsd + k,(i-1)*nsd + k) += mass;
}
// Convection terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
// Sum convection for each direction
conv = 0.0;
for (k = 1;k <= nsd;k++)
conv += U(k)*dNdX(i,k);
conv *= rho*N(j)*detJW;
// Same convection contribution to each velocity component
for (k = 1;k <= nsd;k++)
EM((j-1)*nsd + k,(i-1)*nsd + k) += conv;
}
// Viscous terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
laplace = 0.0;
for (k = 1;k <= nsd;k++)
laplace += dNdX(i,k)*dNdX(j,k);
laplace *= mu*detJW;
// Same viscous contribution to each velocity component
for (k = 1;k <= nsd;k++)
EM((j-1)*nsd + k,(i-1)*nsd + k) += laplace;
// Extra terms for stress-strain formulation
if (formulation == SIM::STRESS)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
EM((j-1)*nsd + l,(i-1)*nsd + k) += mu*dNdX(i,k)*dNdX(j,l)*detJW;
}
}
// Righthand side terms
if (eS) {
Vector& ES = *eS;
Vector fb(g,nsd);
fb *= this->getMassDensity(X)*detJW;
// Vector fb(nsd);
// this->getBodyForce(X,fb);
// fb *= detJW;
for (i = 1;i <= N.size();i++) {
// Time derivative terms
mass = rho*dtInv*N(i)*detJW;
for (k = 1;k <= nsd;k++)
ES((i-1)*nsd + k) += mass*(-B1*UP(k)-B2*UPP(k));
// Body forces
for (k = 1;k <= nsd;k++)
ES((i-1)*nsd + k) += fb(k)*N(i);
// Pressure term
for (k = 1;k <= nsd;k++)
ES((i-1)*nsd + k) += P*dNdX(i,k)*detJW;
//ES((i-1)*nsd + k) -= dPdX(k)*N(i)*detJW;
}
}
return getIntegralResult(elmInt);
}
bool ChorinVelPredBDF2::evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Matrix3D& d2NdX2, const Vec3& X,
double h) const
{
int i, j;
int k, l, m;
int dof;
real P, Uabs, delta;
real dtInv, dUdt, temp;
real conv, convI, convJ, diff, laplace;
real B0, B1, B2;
Vector U(nsd), UP(nsd), UPP(nsd), dPdX(nsd);
Matrix dUdX(nsd,nsd);
Matrix3D d2UdX2(nsd,nsd,nsd);
// Inverse of timestep size
if (time.dt > 0.0)
dtInv = 1.0/time.dt;
else
dtInv = 0.0;
// Element velocity vectors
Vector& EVP = *eVs[1];
Vector& EVPP = *eVs[2];
// Element pressure vector
Vector& EP = *ePs[1];
Vector& EPP = *ePs[2];
if (fabs(time.t-time.dt) < 1.0e-8) {
EVPP = EVP;
EPP = EP;
B0 = 1.0;
B1 = -1.0;
B2 = 0.0;
}
else {
B0 = 1.5;
B1 = -2.0;
B2 = 0.5;
}
// Previous velocities in integration point
UP.fill(0.0); UPP.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++) {
dof = (i-1)*nsd + k;
UP(k) += EVP(dof)*N(i);
UPP(k) += EVPP(dof)*N(i);
}
// Convective velocity in integration point
U = UP;
U *= 2.0;
U -= UPP;
// Velocity magnitude
Uabs = U.norm2();
// Pressure in integration point
P = 0.0;
for (i = 1;i <= N.size();i++)
P += EP(i)*N(i);
// Velocity gradient in integration point
dUdX.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
dUdX(k,l) += dNdX(i,l)*EVP((i-1)*nsd + k);
// Pressure gradient in integration point
dPdX.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
dPdX(k) += dNdX(i,k)*EP(i);
// Hessian for velocity
d2UdX2.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
for (m = 1;m <= nsd;m++)
d2UdX2(k,m,l) += d2NdX2(i,m,l)*EVP((i-1)*nsd + k);
// Stabilization parameters
if (mu < rho*Uabs*h) {
if (time.dt == 0)
delta = 0.5/sqrt(1.0/(h*h) + (Uabs*Uabs)/(h*h));
else
delta = 0.5/sqrt(dtInv*dtInv + (Uabs*Uabs)/(h*h));
}
delta *= detJW;
// Loop over basis functions
// Lefthand size terms
if (eM) {
Matrix& EM = *eM;
// Time derivative terms
temp = rho*dtInv*detJW;
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
dUdt = B0*temp*N(i)*N(j);
for (k = 1;k <= nsd;k++)
EM((j-1)*nsd + k,(i-1)*nsd + k) += dUdt;
}
// Convection terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
// Sum convection for each direction
conv = 0.0;
for (k = 1;k <= nsd;k++)
conv += U(k)*dNdX(i,k);
conv *= rho*N(j)*detJW;
// Same convection contribution to each velocity component
for (k = 1;k <= nsd;k++)
EM((j-1)*nsd + k,(i-1)*nsd + k) += conv;
}
// Viscous terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
laplace = 0.0;
for (k = 1;k <= nsd;k++)
laplace += dNdX(i,k)*dNdX(j,k);
laplace *= mu*detJW;
// Same viscous contribution to each velocity component
for (k = 1;k <= nsd;k++)
EM((j-1)*nsd + k,(i-1)*nsd + k) += laplace;
// Extra terms for stress-strain formulation
if (formulation == SIM::STRESS)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
EM((j-1)*nsd + l,(i-1)*nsd + k) += mu*dNdX(i,k)*dNdX(j,l)*detJW;
}
// SUPG stabilization terms
for (i = 1;i <= N.size();i++) {
// Convection for solution basis functions
convI = 0.0;
for (k = 1;k <= nsd;k++)
convI += U(k)*dNdX(i,k);
convI *= rho;
for (j = 1;j <= N.size();j++) {
// Convection for test functions
convJ = 0.0;
for (k = 1;k <= nsd;k++)
convJ += U(k)*dNdX(j,k);
convJ *= rho*delta;
// Convection terms
for (k = 1;k <= nsd;k++) {
EM((j-1)*nsd + k,(i-1)*nsd + k) += convI*convJ;
}
// Viscous terms
diff = 0.0;
for (k = 1;k <= nsd;k++)
diff += d2NdX2(i,k,k);
diff *= mu;
for (k = 1;k <= nsd;k++)
EM((j-1)*nsd + k,(i-1)*nsd + k) -= diff*convJ;
}
}
}
// Righthand side terms
if (eS) {
Vector& ES = *eS;
Vector fb(g,nsd);
fb *= this->getMassDensity(X);
for (i = 1;i <= N.size();i++) {
// Time derivative terms
dUdt = rho*N(i)*dtInv*detJW;
for (k = 1;k <= nsd;k++)
ES((i-1)*nsd + k) += dUdt*(-B1*UP(k)-B2*UPP(k));
// Convection for solution basis functions
convI = 0.0;
for (k = 1;k <= nsd;k++)
convI += U(k)*dNdX(i,k);
convI *= rho*delta;
// Body forces
for (k = 1;k <= nsd;k++)
ES((i-1)*nsd + l) += fb(k)*(N(i)*detJW + convI);
// Pressure terms
for (k = 1;k <= nsd;k++)
ES((i-1)*nsd + k) += P*dNdX(i,k)*detJW - dPdX(k)*convI;
//ES((i-1)*nsd + k) -= dPdX(k)*N(i)*detJW;
}
}
return getIntegralResult(elmInt);
}
bool ChorinVelPredBDF2::evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N1, const Vector& N2,
const Matrix& dN1dX, const Matrix& dN2dX,
const Vec3& X) const
{
int i, j;
int k, l;
int dof;
real dtInv, temp, mass;
real P, conv, diff, laplace;
real B0, B1, B2;
Vector U(nsd), UP(nsd), UPP(nsd), dPdX(nsd);
// Inverse of timestep size
if (time.dt > 0.0)
dtInv = 1.0/time.dt;
else
dtInv = 0.0;
// Element velocity vectors
Vector& EVP = *eVs[1];
Vector& EVPP = *eVs[2];
// Element pressure vector
Vector& EP = *ePs[1];
Vector& EPP = *ePs[2];
if (fabs(time.t-time.dt) < 1.0e-8) {
EVPP = EVP;
EPP = EP;
B0 = 1.0;
B1 = -1.0;
B2 = 0.0;
}
else {
B0 = 1.5;
B1 = -2.0;
B2 = 0.5;
}
// Velocity in integration point at previous timesteps
UP.fill(0.0); UPP.fill(0.0);
for (i = 1;i <= N1.size();i++)
for (k = 1;k <= nsd;k++) {
dof = (i-1)*nsd + k;
UP(k) += EVP(dof)*N1(i);
UPP(k) += EVPP(dof)*N1(i);
}
// Convective velocity in integration point
U = UP;
U *= 2.0;
U -= UPP;
// Pressure in integration point
P = 0.0;
for (i = 1;i <= N2.size();i++)
P += EP(i)*N2(i);
// Loop over basis functions
// Lefthand size terms
if (eM) {
Matrix& EM = *eM;
// Time derivative terms
temp = rho*dtInv*detJW;
for (i = 1;i <= N1.size();i++)
for (j = 1;j <= N1.size();j++) {
mass = B0*temp*N1(i)*N1(j);
for (k = 1;k <= nsd;k++)
EM((j-1)*nsd + k,(i-1)*nsd + k) += mass;
}
// Convection terms
for (i = 1;i <= N1.size();i++)
for (j = 1;j <= N1.size();j++) {
// Sum convection for each direction
conv = 0.0;
for (k = 1;k <= nsd;k++)
conv += U(k)*dN1dX(i,k);
conv *= rho*N1(j)*detJW;
// Same convection contribution to each velocity component
for (k = 1;k <= nsd;k++)
EM((j-1)*nsd + k,(i-1)*nsd + k) += conv;
}
// Viscous terms
for (i = 1;i <= N1.size();i++)
for (j = 1;j <= N1.size();j++) {
laplace = 0.0;
for (k = 1;k <= nsd;k++)
laplace += dN1dX(i,k)*dN1dX(j,k);
laplace *= mu*detJW;
// Same viscous contribution to each velocity component
for (k = 1;k <= nsd;k++)
EM((j-1)*nsd + k,(i-1)*nsd + k) += laplace;
// Extra terms for stress-strain formulation
if (formulation == SIM::STRESS)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
EM((j-1)*nsd + l,(i-1)*nsd + k) += mu*dN1dX(i,k)*dN1dX(j,l)*detJW;
}
}
// Righthand side terms
if (eS) {
Vector& ES = *eS;
Vector fb(g,nsd);
fb *= this->getMassDensity(X)*detJW;
for (i = 1;i <= N1.size();i++) {
// Time derivative terms
mass = rho*dtInv*N1(i)*detJW;
for (k = 1;k <= nsd;k++)
ES((i-1)*nsd + k) += mass*(-B1*UP(k)-B2*UPP(k));
// Body forces
for (k = 1;k <= nsd;k++)
ES((i-1)*nsd + k) += fb(k)*N1(i);
// Pressure term
for (k = 1;k <= nsd;k++)
ES((i-1)*nsd + k) += P*dN1dX(i,k)*detJW;
//ES((i-1)*nsd + k) -= dPdX(k)*N(i)*detJW;
}
}
return getIntegralResult(elmInt);
}
bool ChorinVelPredBDF2::evalBou(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X, const Vec3& normal) const
{
if (!tracFld) {
std::cerr <<" *** ChorinVelPredBDF2::evalBou: No tractions."<< std::endl;
return false;
}
else if (!eS) {
std::cerr <<" *** ChorinVelPredBDF2::evalBou: No load vector."<< std::endl;
return false;
}
// Evaluate the surface traction
Vec3 T = (*tracFld)(X,normal);
// Store the traction value for vizualization
if (!T.isZero()) tracVal[X] = T;
// Integrate the force vector
Vector& ES = *eS;
for (size_t a = 1; a <= N.size(); a++)
for (short int d = 1; d <= nsd; d++)
ES(nsd*(a-1)+d) -= T[d-1]*N(a)*detJW;
return getIntegralResult(elmInt);
}

94
src/Integrands/ChorinVelPredBDF2.h
View File

@@ -1,94 +0,0 @@
//==============================================================================
//!
//! \file ChorinVelPredBDF2.h
//!
//! \date Nov 23 2010
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief BDF2 implementation for velocity prediction in Chorin's method
//!
//==============================================================================
#ifndef _CHORIN_VEL_PRED_BDF2_H
#define _CHORIN_VEL_PRED_BDF2_H
#include "ChorinVelPredBE.h"
/*!
\brief Class representing the integrand of the velocity prediction step
in Chorin's method.
\details This class implements the velocity prediction step of Chorin's
method using NURBS based FEM with equal order elements for velocity and
pressure. Semi-implicit BDF2 time discretization.
*/
class ChorinVelPredBDF2 : public ChorinVelPredBE
{
public:
//! \brief The default constructor initializes all pointers to zero.
//! \param[in] n Number of spatial dimensions
//! \param[in] form The solution formulation to use
//! \param[in] itg The integrandtype to use
ChorinVelPredBDF2(short int n, SIM::Formulation form,
int itg, bool mixed = false);
//! \brief The destructor frees dynamically allocated objects.
virtual ~ChorinVelPredBDF2();
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] time Parameters for nonlinear and time-dependent simulations
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
//!
//! \details This interface is used when \a getIntegrandType returns 1.
bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const;
//! \brief Evaluates the integrand at an interior 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] d2NdX2 Basis function second derivatives
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] h Characteristic element length
//!
//! \details This interface is used when \a getIntegrandType returns 2.
bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Matrix3D& d2NdX2, const Vec3& X,
double h = 0.0) const;
//! \brief Evaluates the mixed field problem integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N1 Basis function values, field velocity
//! \param[in] N2 Basis function values, field pressure
//! \param[in] dN1dX Basis function gradients, velocity
//! \param[in] dN2dX Basis function gradients, pressure
//! \param[in] X Cartesian coordinates of current integration point
bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N1, const Vector& N2,
const Matrix& dN1dX, const Matrix& dN2dX,
const Vec3& X) 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] 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
bool evalBou(LocalIntegral*& elmInt, const TimeDomain& time,
double detJW, const Vector& N, const Matrix& dNdX,
const Vec3& X, const Vec3& normal) const;
};
#endif

519
src/Integrands/ChorinVelPredBE.C
View File

@@ -1,519 +0,0 @@
//==============================================================================
//!
//! \file ChorinVelPredBE.C
//!
//! \date Sep 29 2010
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Integrand implementations for velocity prediction in Chorin's method
//!
//==============================================================================
#include "ChorinVelPredBE.h"
#include "Tensor.h"
#include "Vec3Oper.h"
#include "Utilities.h"
ChorinVelPredBE::ChorinVelPredBE(short int n, SIM::Formulation form,
int itg, bool incPress, bool mixed)
: ChorinVelPred(n,form,itg,incPress,mixed)
{
int i;
// Resize velocity element solution vectors
int nUsols = 2;
eVs.resize(nUsols);
for (i = 0;i < nUsols;i++)
eVs[i] = new Vector();
primsol.resize(nUsols);
// Resize pressure element solution vectors
int nPsols = 2;
ePs.resize(nPsols);
for (i = 0;i < nPsols;i++)
ePs[i] = new Vector();
psol.resize(nPsols);
}
ChorinVelPredBE::~ChorinVelPredBE() {}
bool ChorinVelPredBE::evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const
{
int i, j;
int k, l;
double dtInv, temp, mass, P;
double conv, diff, laplace;
Vector U(nsd), UP(nsd), dPdX(nsd);
// Inverse of timestep size
if (time.dt > 0.0)
dtInv = 1.0/time.dt;
else
dtInv = 0.0;
// Element velocity vectors
Vector& EV = *eVs[0];
Vector& EVP = *eVs[0];
// Velocity in integration point
U.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
U(k) += EV((i-1)*nsd + k)*N(i);
// Velocity in integration point at previous timestep
UP.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
UP(k) += EVP((i-1)*nsd + k)*N(i);
// Element pressure vector
Vector& EP = *ePs[0];
// Pressure in integration point
P = 0.0;
for (i = 1;i <= N.size();i++)
P += EP(i)*N(i);
// Pressure gradient in integration point
// dPdX.fill(0.0);
// for (i = 1;i <= N.size();i++)
// for (k = 1;k <= nsd;k++)
// dPdX(k) += EP(i)*dNdX(i,k);
// Loop over basis functions
// Lefthand size terms
if (eM) {
Matrix& EM = *eM;
// Time derivative terms
temp = rho*dtInv*detJW;
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
mass = temp*N(i)*N(j);
for (k = 1;k <= nsd;k++)
EM((j-1)*nsd + k,(i-1)*nsd + k) += mass;
}
// Convection terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
// Sum convection for each direction
conv = 0.0;
for (k = 1;k <= nsd;k++)
conv += U(k)*dNdX(i,k);
conv *= rho*N(j)*detJW;
// Same convection contribution to each velocity component
for (k = 1;k <= nsd;k++)
EM((j-1)*nsd + k,(i-1)*nsd + k) += conv;
}
// Viscous terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
laplace = 0.0;
for (k = 1;k <= nsd;k++)
laplace += dNdX(i,k)*dNdX(j,k);
laplace *= mu*detJW;
// Same viscous contribution to each velocity component
for (k = 1;k <= nsd;k++)
EM((j-1)*nsd + k,(i-1)*nsd + k) += laplace;
// Extra terms for stress-strain formulation
if (formulation == SIM::STRESS)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
EM((j-1)*nsd + l,(i-1)*nsd + k) += mu*dNdX(i,k)*dNdX(j,l)*detJW;
}
}
// Righthand side terms
if (eS) {
Vector& ES = *eS;
Vector fb(g,nsd);
fb *= this->getMassDensity(X)*detJW;
// Vector fb(nsd);
// this->getBodyForce(X,fb);
// fb *= detJW;
for (i = 1;i <= N.size();i++) {
// Time derivative terms
mass = rho*dtInv*N(i)*detJW;
for (k = 1;k <= nsd;k++)
ES((i-1)*nsd + k) += mass*UP(k);
// Body forces
for (k = 1;k <= nsd;k++)
ES((i-1)*nsd + k) += fb(k)*N(i);
}
// Only for incremental pressure
if (incPressure)
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
ES((i-1)*nsd + k) += P*dNdX(i,k)*detJW;
//ES((i-1)*nsd + k) -= dPdX(k)*N(i)*detJW;
}
return getIntegralResult(elmInt);
}
bool ChorinVelPredBE::evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Matrix3D& d2NdX2, const Vec3& X,
double h) const
{
int i, j;
int k, l, m;
double P, Uabs, delta;
double dtInv, dUdt, temp;
double conv, convI, convJ, diff, laplace;
Vector U(nsd), UP(nsd), dPdX(nsd), res(nsd);
Matrix dUdX(nsd,nsd);
Matrix3D d2UdX2(nsd,nsd,nsd);
// Inverse of timestep size
if (time.dt > 0.0)
dtInv = 1.0/time.dt;
else
dtInv = 0.0;
// Element velocity vectors
Vector& EV = *eVs[0];
Vector& EVP = *eVs[0];
// Velocity in integration point
U.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
U(k) += EV((i-1)*nsd + k)*N(i);
// Velocity magnitude
Uabs = U.norm2();
// Previous velocity in integration point
UP.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
UP(k) += EVP((i-1)*nsd + k)*N(i);
// Element pressure vector
Vector& EP = *ePs[0];
// Pressure in integration point
P = 0.0;
for (i = 1;i <= N.size();i++)
P += EP(i)*N(i);
// Velocity gradient in integration point
dUdX.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
dUdX(k,l) += dNdX(i,l)*EV((i-1)*nsd + k);
// Pressure gradient in integration point
dPdX.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
dPdX(k) += dNdX(i,k)*EP(i);
// Hessian for velocity
d2UdX2.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
for (m = 1;m <= nsd;m++)
d2UdX2(k,m,l) += d2NdX2(i,m,l)*EV((i-1)*nsd + k);
// Stabilization parameters
if (mu < rho*Uabs*h) {
if (time.dt == 0)
delta = 0.5/sqrt(1.0/(h*h) + (Uabs*Uabs)/(h*h));
else
delta = 0.5/sqrt(dtInv*dtInv + (Uabs*Uabs)/(h*h));
}
delta *= detJW;
// Loop over basis functions
// Lefthand size terms
if (eM) {
Matrix& EM = *eM;
// Time derivative terms
temp = rho*dtInv*detJW;
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
dUdt = temp*N(i)*N(j);
for (k = 1;k <= nsd;k++)
EM((j-1)*nsd + k,(i-1)*nsd + k) += dUdt;
}
// Convection terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
// Sum convection for each direction
conv = 0.0;
for (k = 1;k <= nsd;k++)
conv += U(k)*dNdX(i,k);
conv *= rho*N(j)*detJW;
// Same convection contribution to each velocity component
for (k = 1;k <= nsd;k++)
EM((j-1)*nsd + k,(i-1)*nsd + k) += conv;
}
// Viscous terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
laplace = 0.0;
for (k = 1;k <= nsd;k++)
laplace += dNdX(i,k)*dNdX(j,k);
laplace *= mu*detJW;
// Same viscous contribution to each velocity component
for (k = 1;k <= nsd;k++)
EM((j-1)*nsd + k,(i-1)*nsd + k) += laplace;
// Extra terms for stress-strain formulation
if (formulation == SIM::STRESS)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
EM((j-1)*nsd + l,(i-1)*nsd + k) += mu*dNdX(i,k)*dNdX(j,l)*detJW;
}
// SUPG stabilization terms
for (i = 1;i <= N.size();i++) {
// Convection for solution basis functions
convI = 0.0;
for (k = 1;k <= nsd;k++)
convI += U(k)*dNdX(i,k);
convI *= rho;
for (j = 1;j <= N.size();j++) {
// Convection for test functions
convJ = 0.0;
for (k = 1;k <= nsd;k++)
convJ += U(k)*dNdX(j,k);
convJ *= rho*delta;
// Convection terms
for (k = 1;k <= nsd;k++)
EM((j-1)*nsd + k,(i-1)*nsd + k) += convI*convJ;
// Viscous terms
diff = 0.0;
for (k = 1;k <= nsd;k++)
diff += d2NdX2(i,k,k);
diff *= mu;
for (k = 1;k <= nsd;k++)
EM((j-1)*nsd + k,(i-1)*nsd + k) -= diff*convJ;
}
}
}
// Righthand side terms
if (eS) {
Vector& ES = *eS;
Vector fb(g,nsd);
fb *= this->getMassDensity(X);
for (i = 1;i <= N.size();i++) {
// Time derivative terms
dUdt = rho*N(i)*dtInv*detJW;
for (k = 1;k <= nsd;k++)
ES((i-1)*nsd + k) += dUdt*UP(k);
// Body forces
for (k = 1;k <= nsd;k++)
ES((i-1)*nsd + l) += fb(k)*N(i);
// Convection for test functions
convI = 0.0;
for (k = 1;k <= nsd;k++)
convI += U(k)*dNdX(i,k);
convI *= rho*delta;
// Body force
for (k = 1;k <= nsd;k++)
ES((i-1)*nsd + l) += fb(k)*(N(i)*detJW + convI);
// Pressure term
if (incPressure)
for (k = 1;k <= nsd;k++)
ES((i-1)*nsd + k) += P*dNdX(i,k)*detJW - dPdX(k)*convI;
// ES((i-1)*nsd + k) -= dPdX(k)*(N(i)*detJW + convI);
}
}
return getIntegralResult(elmInt);
}
bool ChorinVelPredBE::evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N1, const Vector& N2,
const Matrix& dN1dX, const Matrix& dN2dX,
const Vec3& X) const
{
int i, j;
int k, l;
double dtInv, temp, mass, P;
double conv, diff, laplace;
Vector U(nsd), UP(nsd), dPdX(nsd);
// Inverse of timestep size
if (time.dt > 0.0)
dtInv = 1.0/time.dt;
else
dtInv = 0.0;
// Element velocity vectors
Vector& EV = *eVs[0];
Vector& EVP = *eVs[0];
// Velocity in integration point
U.fill(0.0);
for (i = 1;i <= N1.size();i++)
for (k = 1;k <= nsd;k++)
U(k) += EV((i-1)*nsd + k)*N1(i);
// Velocity in integration point at previous timestep
UP.fill(0.0);
for (i = 1;i <= N1.size();i++)
for (k = 1;k <= nsd;k++)
UP(k) += EVP((i-1)*nsd + k)*N1(i);
// Element pressure vector
Vector& EP = *ePs[0];
// Pressure in integration point
P = 0.0;
for (i = 1;i <= N2.size();i++)
P += EP(i)*N2(i);
// Loop over basis functions
// Lefthand size terms
if (eM) {
Matrix& EM = *eM;
// Time derivative terms
temp = rho*dtInv*detJW;
for (i = 1;i <= N1.size();i++)
for (j = 1;j <= N1.size();j++) {
mass = temp*N1(i)*N1(j);
for (k = 1;k <= nsd;k++)
EM((j-1)*nsd + k,(i-1)*nsd + k) += mass;
}
// Convection terms
for (i = 1;i <= N1.size();i++)
for (j = 1;j <= N1.size();j++) {
// Sum convection for each direction
conv = 0.0;
for (k = 1;k <= nsd;k++)
conv += U(k)*dN1dX(i,k);
conv *= rho*N1(j)*detJW;
// Same convection contribution to each velocity component
for (k = 1;k <= nsd;k++)
EM((j-1)*nsd + k,(i-1)*nsd + k) += conv;
}
// Viscous terms
for (i = 1;i <= N1.size();i++)
for (j = 1;j <= N1.size();j++) {
laplace = 0.0;
for (k = 1;k <= nsd;k++)
laplace += dN1dX(i,k)*dN1dX(j,k);
laplace *= mu*detJW;
// Same viscous contribution to each velocity component
for (k = 1;k <= nsd;k++)
EM((j-1)*nsd + k,(i-1)*nsd + k) += laplace;
// Extra terms for stress-strain formulation
if (formulation == SIM::STRESS)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
EM((j-1)*nsd + l,(i-1)*nsd + k) += mu*dN1dX(i,k)*dN1dX(j,l)*detJW;
}
}
// Righthand side terms
if (eS) {
Vector& ES = *eS;
Vector fb(g,nsd);
fb *= this->getMassDensity(X)*detJW;
for (i = 1;i <= N1.size();i++) {
// Time derivative terms
mass = rho*dtInv*N1(i)*detJW;
for (k = 1;k <= nsd;k++)
ES((i-1)*nsd + k) += mass*UP(k);
// Body forces
for (k = 1;k <= nsd;k++)
ES((i-1)*nsd + k) += fb(k)*N1(i);
}
// Only for incremental pressure
if (incPressure)
for (i = 1;i <= N1.size();i++)
for (k = 1;k <= nsd;k++)
ES((i-1)*nsd + k) += P*dN1dX(i,k)*detJW;
}
return getIntegralResult(elmInt);
}
bool ChorinVelPredBE::evalBou(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X, const Vec3& normal) const
{
if (incPressure) {
if (!tracFld) {
std::cerr <<" *** ChorinVelPredBE::evalBou: No tractions."<< std::endl;
return false;
}
else if (!eS) {
std::cerr <<" *** ChorinVelPredBE::evalBou: No load vector."<< std::endl;
return false;
}
// Evaluate the surface traction
Vec3 T = (*tracFld)(X,normal);
// Store the traction value for vizualization
if (!T.isZero()) tracVal[X] = T;
// Integrate the force vector
Vector& ES = *eS;
for (size_t a = 1; a <= N.size(); a++)
for (short int d = 1; d <= nsd; d++)
ES(nsd*(a-1)+d) -= T[d-1]*N(a)*detJW;
}
return getIntegralResult(elmInt);
}

96
src/Integrands/ChorinVelPredBE.h
View File

@@ -1,96 +0,0 @@
//==============================================================================
//!
//! \file ChorinVelPredBE.h
//!
//! \date Sep 29 2010
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Implicit Euler implementation for velocity prediction in Chorin's method
//!
//==============================================================================
#ifndef _CHORIN_VEL_PRED_BE_H
#define _CHORIN_VEL_PRED_BE_H
#include "ChorinVelPred.h"
/*!
\brief Class representing the integrand of the velocity prediction step
in Chorin's method.
\details This class implements the velocity prediction step of Chorin's
method using NURBS based FEM with equal order elements for velocity and
pressure. Semi-implicit Euler time discretization.
*/
class ChorinVelPredBE : public ChorinVelPred
{
public:
//! \brief The default constructor initializes all pointers to zero.
//! \param[in] n Number of spatial dimensions
//! \param[in] form The solution formulation to use
//! \param[in] itg The integrandtype to use
//! \param[in] incPress \e true if incremental pressure formulation
ChorinVelPredBE(short int n, SIM::Formulation form, int itg,
bool incPress = false, bool mixed = false);
//! \brief The destructor frees dynamically allocated objects.
virtual ~ChorinVelPredBE();
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] time Parameters for nonlinear and time-dependent simulations
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
//!
//! \details This interface is used when \a getIntegrandType returns 1.
bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const;
//! \brief Evaluates the integrand at an interior 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] d2NdX2 Basis function second derivatives
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] h Characteristic element length
//!
//! \details This interface is used when \a getIntegrandType returns 2.
bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Matrix3D& d2NdX2, const Vec3& X,
double h = 0.0) const;
//! \brief Evaluates the mixed field problem integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N1 Basis function values, field velocity
//! \param[in] N2 Basis function values, field pressure
//! \param[in] dN1dX Basis function gradients, velocity
//! \param[in] dN2dX Basis function gradients, pressure
//! \param[in] X Cartesian coordinates of current integration point
virtual bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N1, const Vector& N2,
const Matrix& dN1dX, const Matrix& dN2dX,
const Vec3& X) 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] 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
virtual bool evalBou(LocalIntegral*& elmInt, const TimeDomain& time,
double detJW, const Vector& N, const Matrix& dNdX,
const Vec3& X, const Vec3& normal) const;
};
#endif

624
src/Integrands/NavierStokesG2.C
View File

@@ -1,624 +0,0 @@
// $Id: NavierStokesG2.C,v 1.4 2011-02-08 12:38:23 rho Exp $
//==============================================================================
//!
//! \file NavierStokesG2.C
//!
//! \date May 28 2010
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Integrand implementations for G2 stabilized Navier-Stokes problems.
//!
//==============================================================================
#include "NavierStokesG2.h"
#include "TimeDomain.h"
#include "Utilities.h"
NavierStokesG2::NavierStokesG2(short int n, SIM::Formulation form, int itg)
: StabilizedStokes(n,form,itg)
{
mu = 1.0e-3;
rho = 1.0;
// Need current solution and solution at previous timestep
eVs.resize(2);
for (int i = 0;i < 2;i++)
eVs[i] = new Vector();
primsol.resize(2);
}
NavierStokesG2::~NavierStokesG2()
{
}
bool NavierStokesG2::evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Matrix3D& d2NdX2, const Vec3& X,
double h) const
{
int i, j;
int k, l, m;
real diff, divU, divI, divJ, convI, convJ, laplace, U, P;
real delta1, delta2;
real dtInv, dUdt;
Vector vel(3), pvel(3), conv(3), res(3), dPdX(3);
Matrix dUdX(3,3);
Matrix3D d2UdX2(3,3,3);
// Inverse of timestep size
if (time.dt > 0.0)
dtInv = 1.0/time.dt;
else
dtInv = 0.0;
Vector& EV = *eVs[0];
// Velocity in integration point
vel.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
vel(k) += EV((i-1)*nf + k)*N(i);
U = 0.0;
for (k = 1;k <= nsd;k++)
U += vel(k)*vel(k);
U = sqrt(U);
Vector& EVP = *eVs[1];
// Previous velocity at integration point
pvel.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
pvel(k) += EVP((i-1)*nf + k)*N(i);
// Pressure in integration point
P = 0.0;
for (i = 1;i <= N.size();i++)
P += EV(i*nf)*N(i);
// Velocity gradient in integration point
dUdX.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
dUdX(k,l) += dNdX(i,l)*EV((i-1)*nf+k);
// Pressure gradient in integration point
dPdX.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
dPdX(k) += dNdX(i,k)*EV(i*nf);
// Hessian for velocity
d2UdX2.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
for (m = 1;m <= nsd;m++)
d2UdX2(k,m,l) += d2NdX2(i,m,l)*EV((i-1)*nf + k);
// Stabilization parameters
if (mu < rho*U*h) {
if (time.dt == 0.0)
delta1 = 0.5/sqrt(1.0/(h*h) + (U*U)/(h*h));
else
delta1 = 0.5/sqrt(dtInv*dtInv + (U*U)/(h*h));
delta2 = h;
}
else
delta1 = delta2 = h*h;
delta1 *= 0.000001*detJW;
delta2 *= 0.000001*detJW;
if (eM) {
Matrix& EM = *eM;
// Time derivative term lhs
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
dUdt = rho*N(i)*N(j)*dtInv*detJW;
for (k = 1;k <= nsd;k++)
EM((j-1)*nf + k, (i-1)*nf + k) += dUdt;
}
// Convection terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
// Sum convection for each direction
convJ = 0.0;
for (k = 1;k <= nsd;k++)
convJ += vel(k)*dNdX(i,k);
convJ *= rho*N(j)*detJW;
for (k = 1;k <= nsd;k++) {
EM((j-1)*nf+k,(i-1)*nf+k) += convJ;
for (l = 1;l <= nsd;l++)
EM((j-1)*nf+l,(i-1)*nf+k) += rho*dUdX(l,k)*N(j)*N(i)*detJW;
}
}
// Viscous terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
laplace = 0.0;
for (k = 1;k <= nsd;k++)
laplace += dNdX(i,k)*dNdX(j,k);
laplace *= mu*detJW;
for (k = 1;k <= nsd;k++)
EM(nf*(j-1)+k,nf*(i-1)+k) += laplace;
}
if (formulation == SIM::STRESS)
for (i = 1; i <= N.size(); i++)
for (j = 1; j <= N.size(); j++)
for (k = 1; k <= nsd; k++)
for (l = 1; l <= nsd; l++)
EM(nf*(j-1)+l,nf*(i-1)+k) += mu*dNdX(i,k)*dNdX(j,l)*detJW;
// Pressure and divergence terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++)
for (k = 1;k <= nsd;k++) {
divI = N(j)*dNdX(i,k)*detJW;
EM(nf*(i-1)+k,nf*j) -= divI;
EM(nf*j,nf*(i-1)+k) -= divI;
}
// Stabilization terms
// Computation of non-linear residual
res.fill(0.0);
for (k = 1;k <= nsd;k++) {
for (l = 1;l <= nsd;l++)
res(k) += rho*vel(l)*dUdX(k,l) - mu*d2UdX2(k,l,l);
res(k) += dPdX(k);
res(k) *= delta1;
}
// STABILIZATION TERMS
for (i = 1;i <= N.size();i++) {
// Convection for solution basis functions
convI = 0.0;
for (k = 1;k <= nsd;k++)
convI += vel(k)*dNdX(i,k);
convI *= rho;
for (j = 1;j <= N.size();j++) {
// Convection for test functions
convJ = 0.0;
for (k = 1;k <= nsd;k++)
convJ += vel(k)*dNdX(j,k);
convJ *= delta1*rho;
//--- MOMENTUM STABILIZATION ---
// Convection terms
for (k = 1;k <= nsd;k++) {
EM((j-1)*nf + k,(i-1)*nf + k) += convI*convJ;
for (l = 1;l <= nsd;l++) {
EM((j-1)*nf + l,(i-1)*nf + k) += rho*res(l)*dNdX(j,k)*N(i);
EM((j-1)*nf + l,(i-1)*nf + k) += convJ*rho*dUdX(l,k)*N(i);
}
}
// Pressure terms
for (l = 1;l <= nsd;l++)
EM((j-1)*nf + l,i*nf) += dNdX(i,l)*convJ;
// Viscous terms
diff = 0.0;
for (k = 1;k <= nsd;k++)
diff += d2NdX2(i,k,k);
diff *= mu;
for (k = 1;k <= nsd;k++)
EM((j-1)*nf + k,(i-1)*nf + k) -= diff*convJ;
// Continuity-continuity stabilization
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
EM((j-1)*nf+l,(i-1)*nf+k) += delta2*dNdX(i,k)*dNdX(j,l);
//--- CONTINUITY STABILIZATION ---
// Convection terms
for (k = 1;k <= nsd;k++)
EM(j*nf,(i-1)*nf+k) += delta1*convI*dNdX(j,k);
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
EM(j*nf,(i-1)*nf+k) += delta1*rho*dUdX(l,k)*N(i)*dNdX(j,l);
// Pressure terms
laplace = 0.0;
for (k = 1;k <= nsd;k++)
laplace += dNdX(i,k)*dNdX(j,k);
laplace *= delta1;
EM(j*nf,i*nf) += laplace;
// Viscous terms
diff = 0.0;
for (k = 1;k <= nsd;k++)
diff += d2NdX2(i,k,k);
diff *= delta1*mu;
for (k = 1;k <= nsd;k++)
EM(j*nf,(i-1)*nf+k) -= diff*dNdX(j,k);
}
}
}
if (eS) {
// Integrate rhs vector
Vector& ES = *eS;
Vector fb(g,nsd);
fb *= this->getMassDensity(X)*detJW;
for (i = 1;i <= N.size();i++) {
// --- Residual of Galerkin formulation ---
// -- Momentum equations --
// Time derivative terms
for (k = 1;k <= nsd;k++)
ES((i-1)*nf + k) += rho*(pvel(k)-vel(k))*N(i)*dtInv*detJW;
// Body forces
for (k = 1;k <= nsd;k++)
ES((i-1)*nf+k) += fb(k)*N(i);
// Convection terms
for (k = 1;k <= nsd;k++) {
conv.fill(0.0);
for (m = 1;m <= nsd;m++)
conv(k) += vel(m)*dUdX(k,m);
ES((i-1)*nf + k) -= rho*conv(k)*N(i)*detJW;
}
// Pressure terms
for (k = 1;k <= nsd;k++)
ES((i-1)*nf + k) += P*dNdX(i,k)*detJW;
// Viscous terms
for (k = 1;k <= nsd;k++)
for (m = 1;m <= nsd;m++)
ES((i-1)*nf + k) -= mu*dUdX(k,m)*dNdX(i,m)*detJW;
// -- Continuity equation --
divU = 0.0;
for (m = 1;m <= nsd;m++)
divU += dUdX(m,m);
ES(i*nf) += divU*N(i)*detJW;
//--- MOMENTUM STABILIZATION---
// Convection for solution basis functions
convI = 0.0;
for (k = 1;k <= nsd;k++)
convI += vel(k)*dNdX(i,k);
convI *= rho;
// Convection terms
for (k = 1;k <= nsd;k++)
ES((i-1)*nf+k) -= res(k)*convI;
// Divergence of velocity field and velocity test function
divU = 0.0;
for (k = 1;k <= nsd;k++)
divU += dUdX(k,k);
// Continuity residual term
for (k = 1;k <= nsd;k++)
ES((i-1)*nf+k) -= delta2*divU*dNdX(i,k);
//--- CONTINUITY STABILIZATION ---
// Momentum residual term
for (k = 1;k <= nsd;k++)
ES(i*nf) -= res(k)*dNdX(i,k);
}
}
return getIntegralResult(elmInt);
}
bool NavierStokesG2::evalInt (LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Matrix3D& d2NdX2,
const Vec3& X, double h) const
{
int i, j;
int k, l, m;
real diff, divU, divI, divJ, convI, convJ, laplace, U, P;
real delta1, delta2;
Vector vel(3), conv(3), res(3), dPdX(3);
Matrix dUdX(3,3);
Matrix3D d2UdX2(3,3,3);
Vector& EV = *eVs[0];
// Velocity in integration point
vel.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
vel(k) += EV((i-1)*nf + k)*N(i);
U = 0.0;
for (k = 1;k <= nsd;k++)
U += vel(k)*vel(k);
U = sqrt(U);
// Pressure in integration point
P = 0.0;
for (i = 1;i <= N.size();i++)
P += EV(i*nf)*N(i);
// Velocity gradient in integration point
dUdX.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
dUdX(k,l) += dNdX(i,l)*EV((i-1)*nf+k);
// Pressure gradient in integration point
dPdX.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
dPdX(k) += dNdX(i,k)*EV(i*nf);
// Hessian for velocity
d2UdX2.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
for (m = 1;m <= nsd;m++)
d2UdX2(k,m,l) += d2NdX2(i,m,l)*EV((i-1)*nf + k);
// Stabilization parameters
if (mu < rho*U*h) {
delta1 = 0.5/sqrt(1.0/(h*h) + (U*U)/(h*h));
delta2 = h;
}
else
delta1 = delta2 = h*h;
delta1 *= detJW;
delta2 *= detJW;
if (eM) {
Matrix& EM = *eM;
// Convection terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
// Sum convection for each direction
convJ = 0.0;
for (k = 1;k <= nsd;k++)
convJ += vel(k)*dNdX(i,k);
convJ *= rho*N(j)*detJW;
for (k = 1;k <= nsd;k++) {
EM((j-1)*nf+k,(i-1)*nf+k) += convJ;
for (l = 1;l <= nsd;l++)
EM((j-1)*nf+l,(i-1)*nf+k) += rho*dUdX(l,k)*N(j)*N(i)*detJW;
}
}
// Viscous terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
laplace = 0.0;
for (k = 1;k <= nsd;k++)
laplace += dNdX(i,k)*dNdX(j,k);
laplace *= mu*detJW;
for (k = 1;k <= nsd;k++)
EM(nf*(j-1)+k,nf*(i-1)+k) += laplace;
}
if (formulation == SIM::STRESS)
for (i = 1; i <= N.size(); i++)
for (j = 1; j <= N.size(); j++)
for (k = 1; k <= nsd; k++)
for (l = 1; l <= nsd; l++)
EM(nf*(j-1)+l,nf*(i-1)+k) += mu*dNdX(i,k)*dNdX(j,l)*detJW;
// Pressure and divergence terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++)
for (k = 1;k <= nsd;k++) {
divI = N(j)*dNdX(i,k)*detJW;
EM(nf*(i-1)+k,nf*j) -= divI;
EM(nf*j,nf*(i-1)+k) -= divI;
}
// Stabilization terms
// Computation of non-linear residual
res.fill(0.0);
for (k = 1;k <= nsd;k++) {
for (l = 1;l <= nsd;l++)
res(k) += rho*vel(l)*dUdX(k,l) - mu*d2UdX2(k,l,l);
res(k) += dPdX(k);
res(k) *= delta1;
}
// STABILIZATION TERMS
for (i = 1;i <= N.size();i++) {
// Convection for solution basis functions
convI = 0.0;
for (k = 1;k <= nsd;k++)
convI += vel(k)*dNdX(i,k);
convI *= rho;
for (j = 1;j <= N.size();j++) {
// Convection for test functions
convJ = 0.0;
for (k = 1;k <= nsd;k++)
convJ += vel(k)*dNdX(j,k);
convJ *= delta1*rho;
//--- MOMENTUM STABILIZATION ---
// Convection terms
for (k = 1;k <= nsd;k++) {
EM((j-1)*nf + k,(i-1)*nf + k) += convI*convJ;
for (l = 1;l <= nsd;l++) {
EM((j-1)*nf + l,(i-1)*nf + k) += rho*res(l)*dNdX(j,k)*N(i);
EM((j-1)*nf + l,(i-1)*nf + k) += convJ*rho*dUdX(l,k)*N(i);
}
}
// Pressure terms
for (l = 1;l <= nsd;l++)
EM((j-1)*nf + l,i*nf) += dNdX(i,l)*convJ;
// Viscous terms
diff = 0.0;
for (k = 1;k <= nsd;k++)
diff += d2NdX2(i,k,k);
diff *= mu;
for (k = 1;k <= nsd;k++)
EM((j-1)*nf + k,(i-1)*nf + k) -= diff*convJ;
// Continuity-continuity stabilization
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
EM((j-1)*nf+l,(i-1)*nf+k) += delta2*dNdX(i,k)*dNdX(j,l);
//--- CONTINUITY STABILIZATION ---
// Convection terms
for (k = 1;k <= nsd;k++)
EM(j*nf,(i-1)*nf+k) += delta1*convI*dNdX(j,k);
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
EM(j*nf,(i-1)*nf+k) += delta1*rho*dUdX(l,k)*N(i)*dNdX(j,l);
// Pressure terms
laplace = 0.0;
for (k = 1;k <= nsd;k++)
laplace += dNdX(i,k)*dNdX(j,k);
laplace *= delta1;
EM(j*nf,i*nf) += laplace;
// Viscous terms
diff = 0.0;
for (k = 1;k <= nsd;k++)
diff += d2NdX2(i,k,k);
diff *= delta1*mu;
for (k = 1;k <= nsd;k++)
EM(j*nf,(i-1)*nf+k) -= diff*dNdX(j,k);
}
}
}
if (eS) {
// Integrate rhs vector
Vector& ES = *eS;
Vector fb(g,nsd);
fb *= this->getMassDensity(X)*detJW;
for (i = 1;i <= N.size();i++) {
// --- Residual of Galerkin formulation ---
// -- Momentum equations --
// Body forces
for (k = 1;k <= nsd;k++)
ES((i-1)*nf+k) += fb(k)*N(i);
// Convection terms
for (k = 1;k <= nsd;k++) {
conv.fill(0.0);
for (m = 1;m <= nsd;m++)
conv(k) += vel(m)*dUdX(k,m);
ES((i-1)*nf + k) -= rho*conv(k)*N(i)*detJW;
}
// Pressure terms
for (k = 1;k <= nsd;k++)
ES((i-1)*nf + k) += P*dNdX(i,k)*detJW;
// Viscous terms
for (k = 1;k <= nsd;k++)
for (m = 1;m <= nsd;m++)
ES((i-1)*nf + k) -= mu*dUdX(k,m)*dNdX(i,m)*detJW;
// -- Continuity equation --
divU = 0.0;
for (m = 1;m <= nsd;m++)
divU += dUdX(m,m);
ES(i*nf) += divU*N(i)*detJW;
//--- MOMENTUM STABILIZATION---
// Convection for solution basis functions
convI = 0.0;
for (k = 1;k <= nsd;k++)
convI += vel(k)*dNdX(i,k);
convI *= rho;
// Convection terms
for (k = 1;k <= nsd;k++)
ES((i-1)*nf+k) -= res(k)*convI;
// Divergence of velocity field and velocity test function
divU = 0.0;
for (k = 1;k <= nsd;k++)
divU += dUdX(k,k);
// Continuity residual term
for (k = 1;k <= nsd;k++)
ES((i-1)*nf+k) -= delta2*divU*dNdX(i,k);
//--- CONTINUITY STABILIZATION ---
// Momentum residual term
for (k = 1;k <= nsd;k++)
ES(i*nf) -= res(k)*dNdX(i,k);
}
}
return getIntegralResult(elmInt);
}
bool NavierStokesG2::evalSol (Vector& s, const Vector& N,
const Matrix& dNdX, const Vec3& X,
const std::vector<int>& MNPC) const
{
if (primsol.empty())
{
std::cerr <<" *** NavierStokesG2::evalSol: No primary solution."
<< std::endl;
return false;
}
int ierr = 0;
if (ierr = utl::gather(MNPC,nsd,primsol[0],*eVs[0]))
{
std::cerr <<" *** NavierStokesG2::evalSol: Detected "
<< ierr <<" node numbers out of range."<< std::endl;
return false;
}
return true;
}

136
src/Integrands/NavierStokesG2.h
View File

@@ -1,136 +0,0 @@
//==============================================================================
//!
//! \file NavierStokesG2.h
//!
//! \date May 28 2010
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Integrand implementations for Navier-Stokes equations with G2 stabilization.
//!
//==============================================================================
#ifndef _NAVIER_STOKES_G2_H
#define _NAVIER_STOKES_G2_H
#include "StabilizedStokes.h"
#include "SIMenums.h"
/*!
\brief Class representing the integrand of the pressure stabilized Navier-Stokes problem.
\details This class supports G2-stabilized Navier-Stokes formulations using equal
order elements for velocity and pressure.
*/
class NavierStokesG2 : public StabilizedStokes
{
public:
//! \brief The default constructor initializes all pointers to zero
//! \param[in] n Number of spatial dimensions
NavierStokesG2(short int n, SIM::Formulation form, int itg = 3);
//! \brief The destructor clears the static work arrays used internally.
virtual ~NavierStokesG2();
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] time Parameters for nonlinear and time-dependent simulations
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
//!
//! \details This interface is used when \a getIntegrandType returns 1.
virtual bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const
{ return false; }
//! \brief Evaluates the integrand at an interior 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] d2NdX2 Basis function second derivatives
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] h Characteristic element length
//!
//! \details This interface is used when \a getIntegrandType returns 2.
//! The default implementation forwards to the stationary version.
//! Reimplement this method for time-dependent or non-linear problems.
virtual bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Matrix3D& d2NdX2, const Vec3& X,
double h = 0.0) const;
//! \brief Evaluates the integrand at an interior 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] Navg Volume-averaged basis function values over the element
//! \param[in] X Cartesian coordinates of current integration point
//!
//! \details This interface is used when \a getIntegrandType returns 3.
//! The default implementation forwards to the stationary version.
//! Reimplement this method for time-dependent or non-linear problems.
virtual bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Vector& Navg, const Vec3& X) const
{ return true; }
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const
{ return false; }
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] d2NdX2 Basis function second derivatives
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] h characteristic element length
virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Matrix3D& d2NdX2,
const Vec3& X, double h = 0.0) const;
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] Navg Average value of basis function over element
//! \param[in] X Cartesian coordinates of current integration point
virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vector& Navg, const Vec3& X) const
{ return false; }
//! \brief Evaluates the secondary solution at current integration point.
//! \param[out] s The solution field values
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] MNPC Matrix of nodal point correspondance
virtual bool evalSol(Vector& s,
const Vector& N, const Matrix& dNdX,
const Vec3& X, const std::vector<int>& MNPC) const;
//! \brief Returns which integrand to use
int getIntegrandType() const { return 2; }
};
#endif

639
src/Integrands/NavierStokesG2CN.C
View File

@@ -1,639 +0,0 @@
// $Id: NavierStokesG2CN.C,v 1.3 2011-02-08 12:39:20 rho Exp $
//==============================================================================
//!
//! \file NavierStokesG2CN.C
//!
//! \date August 9 2010
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Integrand implementations for G2 stabilized Navier-Stokes problems.
//! \details Time integration by the Crank-Nicholson method.
//!
//==============================================================================
#include "NavierStokesG2CN.h"
#include "TimeDomain.h"
#include "ElmNorm.h"
NavierStokesG2CN::NavierStokesG2CN(short int n, SIM::Formulation form, int itg)
: NavierStokesG2(n,form,itg) {}
NavierStokesG2CN::~NavierStokesG2CN() {}
bool NavierStokesG2CN::evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Matrix3D& d2NdX2, const Vec3& X,
double h) const
{
int i, j;
int k, l, m;
real temp1, temp2;
real diff, divU, divI, divJ, convI, convJ, laplace, U, P;
real pconvI, pdiff, Pn;
real delta1, delta2;
real dtInv, dUdt;
Vector vel(3), conv(3), visc(3), res(3), dPdX(3);
Vector pvel(3), resp(3), dPdXn(3);
Matrix dUdX(3,3), dUdXn(3,3);
Matrix3D d2UdX2(3,3,3), d2UdX2n(3,3,3);
const int nf = nsd + 1;
// Inverse of timestep size
if (time.dt > 0.0)
dtInv = 1.0/time.dt;
else
dtInv = 0.0;
Vector& EV = *eVs[0];
// Velocity in integration point
vel.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
vel(k) += EV((i-1)*nf + k)*N(i);
U = 0.0;
for (k = 1;k <= nsd;k++)
U += vel(k)*vel(k);
U = sqrt(U);
Vector& EVP = *eVs[1];
// Previous velocity at integration point
pvel.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
pvel(k) += EVP((i-1)*nf + k)*N(i);
// Pressure in integration point
P = 0.0;
for (i = 1;i <= N.size();i++)
P += EV(i*nf)*N(i);
// Pressure in integration point at previous timestep
Pn = 0.0;
for (i = 1;i <= N.size();i++)
Pn += EVP(i*nf)*N(i);
// Velocity gradient in integration point
dUdX.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
dUdX(k,l) += dNdX(i,l)*EV((i-1)*nf+k);
// Velocity gradient at previous timestep
dUdXn.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
dUdXn(k,l) += dNdX(i,l)*EVP((i-1)*nf+k);
// Pressure gradient in integration point
dPdX.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
dPdX(k) += dNdX(i,k)*EV(i*nf);
// Pressure gradient in integration point at previous timestep
dPdXn.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
dPdXn(k) += dNdX(i,k)*EVP(i*nf);
// Hessian for velocity
d2UdX2.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
for (m = 1;m <= nsd;m++)
d2UdX2(k,m,l) += d2NdX2(i,m,l)*EV((i-1)*nf + k);
// Hessian for velocity at previous timestep
d2UdX2n.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
for (m = 1;m <= nsd;m++)
d2UdX2n(k,m,l) += d2NdX2(i,m,l)*EVP((i-1)*nf + k);
// Stabilization parameters
if (mu < rho*U*h) {
if (time.dt == 0.0)
delta1 = 0.5/sqrt(1.0/(h*h) + (U*U)/(h*h));
else
delta1 = 0.5/sqrt(dtInv*dtInv + (U*U)/(h*h));
delta2 = h;
}
else
delta1 = delta2 = h*h;
delta1 *= detJW;
delta2 *= detJW;
if (eM) {
Matrix& EM = *eM;
// Time derivative term lhs
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
dUdt = rho*N(i)*N(j)*dtInv*detJW;
for (k = 1;k <= nsd;k++)
EM((j-1)*nf + k, (i-1)*nf + k) += dUdt;
}
// Convection terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
// Sum convection for each direction
convJ = 0.0;
for (k = 1;k <= nsd;k++)
convJ += vel(k)*dNdX(i,k);
convJ *= 0.5*rho*N(j)*detJW;
temp1 = 0.5*rho*N(i)*N(j)*detJW;
for (k = 1;k <= nsd;k++) {
EM((j-1)*nf+k,(i-1)*nf+k) += convJ;
for (l = 1;l <= nsd;l++)
EM((j-1)*nf+l,(i-1)*nf+k) += temp1*dUdX(l,k);
}
}
// Viscous terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
laplace = 0.0;
for (k = 1;k <= nsd;k++)
laplace += dNdX(i,k)*dNdX(j,k);
laplace *= 0.5*mu*detJW;
for (k = 1;k <= nsd;k++)
EM(nf*(j-1)+k,nf*(i-1)+k) += laplace;
}
temp1 = 0.5*mu*detJW;
if (formulation == SIM::STRESS)
for (i = 1; i <= N.size(); i++)
for (j = 1; j <= N.size(); j++)
for (k = 1; k <= nsd; k++)
for (l = 1; l <= nsd; l++)
EM(nf*(j-1)+l,nf*(i-1)+k) += temp1*dNdX(i,k)*dNdX(j,l);
// Pressure and divergence terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++)
for (k = 1;k <= nsd;k++) {
divI = N(j)*dNdX(i,k)*detJW;
EM(nf*(i-1)+k,nf*j) -= 0.5*divI;
EM(nf*j,nf*(i-1)+k) -= divI;
}
// Stabilization terms
// Computation of non-linear residual
res.fill(0.0);
resp.fill(0.0);
for (k = 1;k <= nsd;k++) {
for (l = 1;l <= nsd;l++) {
res(k) += rho*vel(l)*dUdX(k,l) - mu*d2UdX2(k,l,l);
resp(k) += rho*pvel(l)*dUdXn(k,l) - mu*d2UdX2n(k,l,l);
}
res(k) += dPdX(k);
res(k) *= 0.5*delta1;
resp(k) += dPdXn(k);
resp(k) *= 0.5*delta1;
}
// STABILIZATION TERMS
for (i = 1;i <= N.size();i++) {
// Convection for solution basis functions
convI = 0.0;
for (k = 1;k <= nsd;k++)
convI += vel(k)*dNdX(i,k);
convI *= rho;
for (j = 1;j <= N.size();j++) {
// Convection for test functions
convJ = 0.0;
for (k = 1;k <= nsd;k++)
convJ += vel(k)*dNdX(j,k);
convJ *= delta1*rho;
//--- MOMENTUM STABILIZATION ---
// Convection terms
for (k = 1;k <= nsd;k++) {
EM((j-1)*nf + k,(i-1)*nf + k) += 0.5*convI*convJ;
for (l = 1;l <= nsd;l++) {
EM((j-1)*nf + l,(i-1)*nf + k) += rho*res(l)*dNdX(j,k)*N(i);
EM((j-1)*nf + l,(i-1)*nf + k) += 0.5*convJ*rho*dUdX(l,k)*N(i);
}
}
// Pressure terms
for (l = 1;l <= nsd;l++)
EM((j-1)*nf + l,i*nf) += 0.5*dNdX(i,l)*convJ;
// Viscous terms
diff = 0.0;
for (k = 1;k <= nsd;k++)
diff += d2NdX2(i,k,k);
diff *= 0.5*mu;
for (k = 1;k <= nsd;k++)
EM((j-1)*nf + k,(i-1)*nf + k) -= diff*convJ;
// Continuity-continuity stabilization
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
EM((j-1)*nf+l,(i-1)*nf+k) += delta2*dNdX(i,k)*dNdX(j,l);
//--- CONTINUITY STABILIZATION ---
// Convection terms
for (k = 1;k <= nsd;k++)
EM(j*nf,(i-1)*nf+k) += 0.5*delta1*convI*dNdX(j,k);
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
EM(j*nf,(i-1)*nf+k) += 0.5*delta1*rho*dUdX(l,k)*N(i)*dNdX(j,l);
// Pressure terms
laplace = 0.0;
for (k = 1;k <= nsd;k++)
laplace += dNdX(i,k)*dNdX(j,k);
laplace *= 0.5*delta1;
EM(j*nf,i*nf) += laplace;
// Viscous terms
diff = 0.0;
for (k = 1;k <= nsd;k++)
diff += d2NdX2(i,k,k);
diff *= 0.5*delta1*mu;
for (k = 1;k <= nsd;k++)
EM(j*nf,(i-1)*nf+k) -= diff*dNdX(j,k);
}
}
}
if (eS) {
// Integrate rhs vector
Vector& ES = *eS;
Vector fb(g,nsd);
fb *= this->getMassDensity(X)*detJW;
for (i = 1;i <= N.size();i++) {
// --- Residual of Galerkin formulation ---
// -- Momentum equations --
// Time derivative terms
for (k = 1;k <= nsd;k++)
ES((i-1)*nf + k) += rho*(pvel(k)-vel(k))*N(i)*dtInv*detJW;
// Body forces
for (k = 1;k <= nsd;k++)
ES((i-1)*nf+k) += fb(k)*N(i);
// Convection terms
for (k = 1;k <= nsd;k++) {
conv.fill(0.0);
for (m = 1;m <= nsd;m++)
conv(k) += vel(m)*dUdX(k,m) + pvel(m)*dUdXn(k,m);
ES((i-1)*nf + k) -= 0.5*rho*conv(k)*N(i)*detJW;
}
// Pressure terms
for (k = 1;k <= nsd;k++)
ES((i-1)*nf + k) += 0.5*(P+Pn)*dNdX(i,k)*detJW;
// Viscous terms
for (k = 1;k <= nsd;k++)
for (m = 1;m <= nsd;m++)
ES((i-1)*nf + k) -= 0.5*mu*(dUdX(k,m) + dUdXn(k,m))*dNdX(i,m)*detJW;
// -- Continuity equation --
divU = 0.0;
for (m = 1;m <= nsd;m++)
divU += dUdX(m,m);
ES(i*nf) += divU*N(i)*detJW;
//--- MOMENTUM STABILIZATION---
// Convection for solution basis functions
convI = 0.0;
for (k = 1;k <= nsd;k++)
convI += vel(k)*dNdX(i,k);
convI *= rho;
pconvI = 0.0;
for (k = 1;k <= nsd;k++)
pconvI += pvel(k)*dNdX(i,k);
pconvI *= rho;
// Convection terms
for (k = 1;k <= nsd;k++) {
ES((i-1)*nf+k) -= res(k)*convI;
ES((i-1)*nf+k) -= resp(k)*pconvI;
}
// Divergence of velocity field and velocity test function
divU = 0.0;
for (k = 1;k <= nsd;k++)
divU += dUdX(k,k);
// Continuity residual term
for (k = 1;k <= nsd;k++)
ES((i-1)*nf+k) -= delta2*divU*dNdX(i,k);
//--- CONTINUITY STABILIZATION ---
// Momentum residual term
for (k = 1;k <= nsd;k++)
ES(i*nf) -= (res(k)+resp(k))*dNdX(i,k);
}
}
return getIntegralResult(elmInt);
}
bool NavierStokesG2CN::evalInt (LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Matrix3D& d2NdX2,
const Vec3& X, double h) const
{
int i, j;
int k, l, m;
real diff, divU, divI, divJ, convI, convJ, laplace, U, P;
real delta1, delta2;
Vector vel(3), conv(3), res(3), dPdX(3);
Matrix dUdX(3,3);
Matrix3D d2UdX2(3,3,3);
const int nf = nsd + 1;
Vector& EV = *eVs[0];
// Velocity in integration point
vel.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
vel(k) += EV((i-1)*nf + k)*N(i);
U = 0.0;
for (k = 1;k <= nsd;k++)
U += vel(k)*vel(k);
U = sqrt(U);
// Pressure in integration point
P = 0.0;
for (i = 1;i <= N.size();i++)
P += EV(i*nf)*N(i);
// Velocity gradient in integration point
dUdX.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
dUdX(k,l) += dNdX(i,l)*EV((i-1)*nf+k);
// Pressure gradient in integration point
dPdX.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
dPdX(k) += dNdX(i,k)*EV(i*nf);
// Hessian for velocity
d2UdX2.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
for (m = 1;m <= nsd;m++)
d2UdX2(k,m,l) += d2NdX2(i,m,l)*EV((i-1)*nf + k);
// Stabilization parameters
if (mu < rho*U*h) {
delta1 = 0.5/sqrt(1.0/(h*h) + (U*U)/(h*h));
delta2 = h;
}
else
delta1 = delta2 = h*h;
delta1 *= detJW;
delta2 *= detJW;
if (eM) {
Matrix& EM = *eM;
// Convection terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
// Sum convection for each direction
convJ = 0.0;
for (k = 1;k <= nsd;k++)
convJ += vel(k)*dNdX(i,k);
convJ *= rho*N(j)*detJW;
for (k = 1;k <= nsd;k++) {
EM((j-1)*nf+k,(i-1)*nf+k) += convJ;
for (l = 1;l <= nsd;l++)
EM((j-1)*nf+l,(i-1)*nf+k) += rho*dUdX(l,k)*N(j)*N(i)*detJW;
}
}
// Viscous terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
laplace = 0.0;
for (k = 1;k <= nsd;k++)
laplace += dNdX(i,k)*dNdX(j,k);
laplace *= mu*detJW;
for (k = 1;k <= nsd;k++)
EM(nf*(j-1)+k,nf*(i-1)+k) += laplace;
}
if (formulation == SIM::STRESS)
for (i = 1; i <= N.size(); i++)
for (j = 1; j <= N.size(); j++)
for (k = 1; k <= nsd; k++)
for (l = 1; l <= nsd; l++)
EM(nf*(j-1)+l,nf*(i-1)+k) += mu*dNdX(i,k)*dNdX(j,l)*detJW;
// Pressure and divergence terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++)
for (k = 1;k <= nsd;k++) {
divI = N(j)*dNdX(i,k)*detJW;
EM(nf*(i-1)+k,nf*j) -= divI;
EM(nf*j,nf*(i-1)+k) -= divI;
}
// Stabilization terms
// Computation of non-linear residual
res.fill(0.0);
for (k = 1;k <= nsd;k++) {
for (l = 1;l <= nsd;l++)
res(k) += rho*vel(l)*dUdX(k,l) - mu*d2UdX2(k,l,l);
res(k) += dPdX(k);
res(k) *= delta1;
}
// STABILIZATION TERMS
for (i = 1;i <= N.size();i++) {
// Convection for solution basis functions
convI = 0.0;
for (k = 1;k <= nsd;k++)
convI += vel(k)*dNdX(i,k);
convI *= rho;
for (j = 1;j <= N.size();j++) {
// Convection for test functions
convJ = 0.0;
for (k = 1;k <= nsd;k++)
convJ += vel(k)*dNdX(j,k);
convJ *= delta1*rho;
//--- MOMENTUM STABILIZATION ---
// Convection terms
for (k = 1;k <= nsd;k++) {
EM((j-1)*nf + k,(i-1)*nf + k) += convI*convJ;
for (l = 1;l <= nsd;l++) {
// RUNAR
//EM((j-1)*nf + l,(i-1)*nf + k) += rho*res(l)*dNdX(j,k)*N(i);
EM((j-1)*nf + l,(i-1)*nf + k) += convJ*rho*dUdX(l,k)*N(i);
}
}
// Pressure terms
for (l = 1;l <= nsd;l++)
EM((j-1)*nf + l,i*nf) += dNdX(i,l)*convJ;
// Viscous terms
diff = 0.0;
for (k = 1;k <= nsd;k++)
diff += d2NdX2(i,k,k);
diff *= mu;
for (k = 1;k <= nsd;k++)
EM((j-1)*nf + k,(i-1)*nf + k) -= diff*convJ;
// Continuity-continuity stabilization
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
EM((j-1)*nf+l,(i-1)*nf+k) += delta2*dNdX(i,k)*dNdX(j,l);
//--- CONTINUITY STABILIZATION ---
// Convection terms
for (k = 1;k <= nsd;k++)
EM(j*nf,(i-1)*nf+k) += delta1*convI*dNdX(j,k);
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
EM(j*nf,(i-1)*nf+k) += delta1*rho*dUdX(l,k)*N(i)*dNdX(j,l);
// Pressure terms
laplace = 0.0;
for (k = 1;k <= nsd;k++)
laplace += dNdX(i,k)*dNdX(j,k);
laplace *= delta1;
EM(j*nf,i*nf) += laplace;
// Viscous terms
diff = 0.0;
for (k = 1;k <= nsd;k++)
diff += d2NdX2(i,k,k);
diff *= delta1*mu;
for (k = 1;k <= nsd;k++)
EM(j*nf,(i-1)*nf+k) -= diff*dNdX(j,k);
}
}
}
if (eS) {
// Integrate rhs vector
Vector& ES = *eS;
Vector fb(g,nsd);
fb *= this->getMassDensity(X)*detJW;
for (i = 1;i <= N.size();i++) {
// --- Residual of Galerkin formulation ---
// -- Momentum equations --
// Body forces
for (k = 1;k <= nsd;k++)
ES((i-1)*nf+k) += fb(k)*N(i);
// Convection terms
for (k = 1;k <= nsd;k++) {
conv.fill(0.0);
for (m = 1;m <= nsd;m++)
conv(k) += vel(m)*dUdX(k,m);
ES((i-1)*nf + k) -= rho*conv(k)*N(i)*detJW;
}
// Pressure terms
for (k = 1;k <= nsd;k++)
ES((i-1)*nf + k) += P*dNdX(i,k)*detJW;
// Viscous terms
for (k = 1;k <= nsd;k++)
for (m = 1;m <= nsd;m++)
ES((i-1)*nf + k) -= mu*dUdX(k,m)*dNdX(i,m)*detJW;
// -- Continuity equation --
divU = 0.0;
for (m = 1;m <= nsd;m++)
divU += dUdX(m,m);
ES(i*nf) += divU*N(i)*detJW;
//--- MOMENTUM STABILIZATION---
// Convection for solution basis functions
convI = 0.0;
for (k = 1;k <= nsd;k++)
convI += vel(k)*dNdX(i,k);
convI *= rho;
// Convection terms
for (k = 1;k <= nsd;k++)
ES((i-1)*nf+k) -= res(k)*convI;
// Divergence of velocity field and velocity test function
divU = 0.0;
for (k = 1;k <= nsd;k++)
divU += dUdX(k,k);
// Continuity residual term
for (k = 1;k <= nsd;k++)
ES((i-1)*nf+k) -= delta2*divU*dNdX(i,k);
//--- CONTINUITY STABILIZATION ---
// Momentum residual term
for (k = 1;k <= nsd;k++)
ES(i*nf) -= res(k)*dNdX(i,k);
}
}
return getIntegralResult(elmInt);
}

123
src/Integrands/NavierStokesG2CN.h
View File

@@ -1,123 +0,0 @@
//==============================================================================
//!
//! \file NavierStokesG2CN.h
//!
//! \date August 9 2010
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Integrand implementations for Navier-Stokes equations with G2 stabilization
//! and
//!
//==============================================================================
#ifndef _NAVIER_STOKES_G2_CN_H
#define _NAVIER_STOKES_G2_CN_H
#include "NavierStokesG2.h"
/*!
\brief Class representing the integrand of the pressure stabilized Navier-Stokes problem.
\details This class supports G2-stabilized Navier-Stokes formulations using equal
order elements for velocity and pressure. Time integration by the Crank-Nicholson method.
*/
class NavierStokesG2CN : public NavierStokesG2
{
public:
//! \brief The default constructor initializes all pointers to zero
//! \param[in] n Number of spatial dimensions
NavierStokesG2CN(short int n, SIM::Formulation form, int itg = 3);
//! \brief The destructor clears the static work arrays used internally.
virtual ~NavierStokesG2CN();
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] time Parameters for nonlinear and time-dependent simulations
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
//!
//! \details This interface is used when \a getIntegrandType returns 1.
virtual bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const
{ return false; }
//! \brief Evaluates the integrand at an interior 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] d2NdX2 Basis function second derivatives
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] h Characteristic element length
//!
//! \details This interface is used when \a getIntegrandType returns 2.
//! The default implementation forwards to the stationary version.
//! Reimplement this method for time-dependent or non-linear problems.
virtual bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Matrix3D& d2NdX2, const Vec3& X,
double h = 0.0) const;
//! \brief Evaluates the integrand at an interior 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] Navg Volume-averaged basis function values over the element
//! \param[in] X Cartesian coordinates of current integration point
//!
//! \details This interface is used when \a getIntegrandType returns 3.
//! The default implementation forwards to the stationary version.
//! Reimplement this method for time-dependent or non-linear problems.
virtual bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Vector& Navg, const Vec3& X) const
{ return true; }
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const
{ return false; }
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] d2NdX2 Basis function second derivatives
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] h characteristic element length
virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Matrix3D& d2NdX2,
const Vec3& X, double h = 0.0) const;
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] Navg Average value of basis function over element
//! \param[in] X Cartesian coordinates of current integration point
virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vector& Navg, const Vec3& X) const
{ return false; }
};
#endif

660
src/Integrands/NavierStokesG2GenTheta.C
View File

@@ -1,660 +0,0 @@
//$Id: NavierStokesG2GenTheta.C,v 1.3 2011-02-08 12:40:31 rho Exp $
//==============================================================================
//!
//! \file NavierStokesG2GenTheta.C
//!
//! \date August 3 2010
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Integrand implementations for G2 stabilized Navier-Stokes problems.
//! \details Time integration by the generalized theta method.
//!
//==============================================================================
#include "NavierStokesG2GenTheta.h"
#include "TimeDomain.h"
#include "ElmNorm.h"
NavierStokesG2GenTheta::NavierStokesG2GenTheta(short int n, SIM::Formulation form, int itg)
: NavierStokesG2(n,form,itg) {}
NavierStokesG2GenTheta::~NavierStokesG2GenTheta() {}
bool NavierStokesG2GenTheta::evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Matrix3D& d2NdX2, const Vec3& X,
double h) const
{
int i, j;
int k, l, m;
real temp1, temp2;
real diff, divU, divI, divJ, convI, convJ, laplace, U, P;
real delta1, delta2;
real dtInv, dUdt;
Vector vel(3), pvel(3), conv(3), convp(3), diffp(3), res(3), resp(3), dPdX(3);
Matrix dUdX(3,3), dUdXn(3,3);
Matrix3D d2UdX2(3,3,3), d2UdX2n(3,3,3);
// Time integration parameters
const real thetau = thetat;
const real thetap = thetam;
const int nf = nsd + 1;
// Inverse of timestep size
if (time.dt > 0.0)
dtInv = 1.0/time.dt;
else
dtInv = 0.0;
Vector& EV = *eVs[0];
// Velocity in integration point
vel.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
vel(k) += EV((i-1)*nf + k)*N(i);
U = 0.0;
for (k = 1;k <= nsd;k++)
U += vel(k)*vel(k);
U = sqrt(U);
Vector& EVP = *eVs[1];
// Previous velocity at integration point
pvel.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
pvel(k) += EVP((i-1)*nf + k)*N(i);
// Pressure in integration point
P = 0.0;
for (i = 1;i <= N.size();i++)
P += EV(i*nf)*N(i);
// Velocity gradient in integration point
dUdX.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
dUdX(k,l) += dNdX(i,l)*EV((i-1)*nf+k);
// Velocity gradient at previous timestep
dUdXn.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
dUdXn(k,l) += dNdX(i,l)*EVP((i-1)*nf+k);
// Pressure gradient in integration point
dPdX.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
dPdX(k) += dNdX(i,k)*EV(i*nf);
// Hessian for velocity
d2UdX2.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
for (m = 1;m <= nsd;m++)
d2UdX2(k,m,l) += d2NdX2(i,m,l)*EV((i-1)*nf + k);
// Hessian for velocity at previous timestep
d2UdX2n.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
for (m = 1;m <= nsd;m++)
d2UdX2n(k,m,l) += d2NdX2(i,m,l)*EVP((i-1)*nf + k);
// Stabilization parameters
if (mu < rho*U*h) {
if (time.dt == 0.0)
delta1 = 0.5/sqrt(1.0/(h*h) + (U*U)/(h*h));
else
delta1 = 0.5/sqrt(dtInv*dtInv + (U*U)/(h*h));
delta2 = h;
}
else
delta1 = delta2 = h*h;
delta1 *= 0.01*detJW;
delta2 *= 0.01*detJW;
// RUNAR
//delta1 = 0.0001*h*h*detJW;
//delta2 = 0.0001*h*h*detJW;
if (eM) {
Matrix& EM = *eM;
// Time derivative term lhs
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
dUdt = rho*N(i)*N(j)*dtInv*detJW;
for (k = 1;k <= nsd;k++)
EM((j-1)*nf + k, (i-1)*nf + k) += dUdt;
}
// Convection terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
// Sum convection for each direction
convJ = 0.0;
for (k = 1;k <= nsd;k++)
convJ += vel(k)*dNdX(i,k);
convJ *= thetau*rho*N(j)*detJW;
temp1 = thetau*rho*N(i)*N(j)*detJW;
for (k = 1;k <= nsd;k++) {
EM((j-1)*nf+k,(i-1)*nf+k) += convJ;
for (l = 1;l <= nsd;l++)
EM((j-1)*nf+l,(i-1)*nf+k) += temp1*dUdX(l,k);
}
}
// Viscous terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
laplace = 0.0;
for (k = 1;k <= nsd;k++)
laplace += dNdX(i,k)*dNdX(j,k);
laplace *= thetau*mu*detJW;
for (k = 1;k <= nsd;k++)
EM(nf*(j-1)+k,nf*(i-1)+k) += laplace;
}
temp1 = thetau*mu*detJW;
if (formulation == SIM::STRESS)
for (i = 1; i <= N.size(); i++)
for (j = 1; j <= N.size(); j++)
for (k = 1; k <= nsd; k++)
for (l = 1; l <= nsd; l++)
EM(nf*(j-1)+l,nf*(i-1)+k) += temp1*dNdX(i,k)*dNdX(j,l);
// Pressure and divergence terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++)
for (k = 1;k <= nsd;k++) {
divI = N(j)*dNdX(i,k)*detJW;
EM(nf*(i-1)+k,nf*j) -= thetap*divI;
EM(nf*j,nf*(i-1)+k) -= divI;
}
// Stabilization terms
// Computation of non-linear residual
res.fill(0.0);
resp.fill(0.0);
for (k = 1;k <= nsd;k++) {
for (l = 1;l <= nsd;l++) {
res(k) += rho*vel(l)*dUdX(k,l) - mu*d2UdX2(k,l,l);
resp(k) += rho*pvel(l)*dUdXn(k,l) - mu*d2UdX2n(k,l,l);
}
//res(k) += dPdX(k);
res(k) *= delta1;
resp(k) *= delta1;
}
// Convection at previous timestep
convp.fill(0.0);
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
convp(k) += rho*pvel(l)*dUdXn(k,l);
// Viscous terms at previous timestep
diffp.fill(0.0);
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
diffp(k) += mu*d2UdX2(k,l,l);
// STABILIZATION TERMS
for (i = 1;i <= N.size();i++) {
// Convection for solution basis functions
convI = 0.0;
for (k = 1;k <= nsd;k++)
convI += vel(k)*dNdX(i,k);
convI *= rho;
for (j = 1;j <= N.size();j++) {
// Convection for test functions
convJ = 0.0;
for (k = 1;k <= nsd;k++)
convJ += vel(k)*dNdX(j,k);
convJ *= delta1*rho;
//--- MOMENTUM STABILIZATION ---
// Convection terms
for (k = 1;k <= nsd;k++) {
EM((j-1)*nf + k,(i-1)*nf + k) += thetau*convI*convJ;
for (l = 1;l <= nsd;l++) {
EM((j-1)*nf + l,(i-1)*nf + k) += thetau*rho*res(l)*dNdX(j,k)*N(i);
EM((j-1)*nf + l,(i-1)*nf + k) += thetau*convJ*rho*dUdX(l,k)*N(i);
EM((j-1)*nf + l,(i-1)*nf + k) += delta1*alpha*thetap*convp(l)*rho*N(i)*dNdX(j,k);
}
}
// Pressure terms
for (l = 1;l <= nsd;l++) {
EM((j-1)*nf + l,i*nf) += thetap*dNdX(i,l)*convJ;
for (k = 1;k <= nsd;k++)
EM((j-1)*nf + l,(i-1)*nf + k) += delta1*thetap*dPdX(l)*rho*N(i)*dNdX(j,k);
}
// Viscous terms
diff = 0.0;
for (k = 1;k <= nsd;k++)
diff += d2NdX2(i,k,k);
diff *= thetau*mu;
for (k = 1;k <= nsd;k++) {
EM((j-1)*nf + k,(i-1)*nf + k) -= diff*convJ;
for (l = 1;l <= nsd;l++)
EM((j-1)*nf + l,(i-1)*nf + k) -= delta1*alpha*thetap*diffp(l)*rho*N(i)*dNdX(j,k);
}
// Continuity-continuity stabilization
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
EM((j-1)*nf+l,(i-1)*nf+k) += delta2*dNdX(i,k)*dNdX(j,l);
//--- CONTINUITY STABILIZATION ---
// Convection terms
for (k = 1;k <= nsd;k++)
EM(j*nf,(i-1)*nf+k) += thetau*delta1*convI*dNdX(j,k);
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
EM(j*nf,(i-1)*nf+k) += thetau*delta1*rho*dUdX(l,k)*N(i)*dNdX(j,l);
// Pressure terms
laplace = 0.0;
for (k = 1;k <= nsd;k++)
laplace += dNdX(i,k)*dNdX(j,k);
laplace *= thetap*delta1;
EM(j*nf,i*nf) += laplace;
// Viscous terms
diff = 0.0;
for (k = 1;k <= nsd;k++)
diff += d2NdX2(i,k,k);
diff *= thetau*delta1*mu;
for (k = 1;k <= nsd;k++)
EM(j*nf,(i-1)*nf+k) -= diff*dNdX(j,k);
}
}
}
if (eS) {
// Integrate rhs vector
Vector& ES = *eS;
Vector fb(g,nsd);
fb *= this->getMassDensity(X)*detJW;
for (i = 1;i <= N.size();i++) {
// --- Residual of Galerkin formulation ---
// -- Momentum equations --
// Time derivative terms
for (k = 1;k <= nsd;k++)
ES((i-1)*nf + k) += rho*(pvel(k)-vel(k))*N(i)*dtInv*detJW;
// Body forces
for (k = 1;k <= nsd;k++)
ES((i-1)*nf+k) += thetap*fb(k)*N(i);
// Convection terms
for (k = 1;k <= nsd;k++) {
conv.fill(0.0);
for (m = 1;m <= nsd;m++)
conv(k) += vel(m)*dUdX(k,m);
ES((i-1)*nf + k) -= thetau*rho*conv(k)*N(i)*detJW;
}
// Pressure terms
for (k = 1;k <= nsd;k++)
ES((i-1)*nf + k) += thetap*P*dNdX(i,k)*detJW;
// Viscous terms
for (k = 1;k <= nsd;k++)
for (m = 1;m <= nsd;m++)
ES((i-1)*nf + k) -= mu*dUdX(k,m)*dNdX(i,m)*detJW;
// -- Continuity equation --
divU = 0.0;
for (m = 1;m <= nsd;m++)
divU += dUdX(m,m);
ES(i*nf) += divU*N(i)*detJW;
//--- MOMENTUM STABILIZATION---
// Convection for solution basis functions
convI = 0.0;
for (k = 1;k <= nsd;k++)
convI += vel(k)*dNdX(i,k);
convI *= rho;
// Convection terms
for (k = 1;k <= nsd;k++) {
ES((i-1)*nf+k) -= thetau*res(k)*convI;
ES((i-1)*nf+k) -= alpha*thetap*convp(k)*N(i)*detJW;
ES((i-1)*nf+k) -= delta1*alpha*thetap*convp(k)*convI;
}
// Viscous terms
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++) {
ES((i-1)*nf+k) += alpha*thetap*mu*dUdXn(k,l)*dNdX(j,l)*detJW;
ES((i-1)*nf+k) += delta1*alpha*thetap*mu*diffp(k)*convI;
}
// Divergence of velocity field and velocity test function
divU = 0.0;
for (k = 1;k <= nsd;k++)
divU += dUdX(k,k);
// Continuity residual term
for (k = 1;k <= nsd;k++)
ES((i-1)*nf+k) -= delta2*divU*dNdX(i,k);
//--- CONTINUITY STABILIZATION ---
// Momentum residual term
for (k = 1;k <= nsd;k++) {
ES(i*nf) -= thetau*res(k)*dNdX(i,k);
ES(i*nf) -= alpha*thetap*resp(k)*dNdX(i,k);
}
}
}
return getIntegralResult(elmInt);
}
bool NavierStokesG2GenTheta::evalInt (LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Matrix3D& d2NdX2,
const Vec3& X, double h) const
{
int i, j;
int k, l, m;
real diff, divU, divI, divJ, convI, convJ, laplace, U, P;
real delta1, delta2;
Vector vel(3), conv(3), res(3), dPdX(3);
Matrix dUdX(3,3);
Matrix3D d2UdX2(3,3,3);
const int nf = nsd + 1;
Vector& EV = *eVs[0];
// Velocity in integration point
vel.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
vel(k) += EV((i-1)*nf + k)*N(i);
U = 0.0;
for (k = 1;k <= nsd;k++)
U += vel(k)*vel(k);
U = sqrt(U);
// Pressure in integration point
P = 0.0;
for (i = 1;i <= N.size();i++)
P += EV(i*nf)*N(i);
// Velocity gradient in integration point
dUdX.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
dUdX(k,l) += dNdX(i,l)*EV((i-1)*nf+k);
// Pressure gradient in integration point
dPdX.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
dPdX(k) += dNdX(i,k)*EV(i*nf);
// Hessian for velocity
d2UdX2.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
for (m = 1;m <= nsd;m++)
d2UdX2(k,m,l) += d2NdX2(i,m,l)*EV((i-1)*nf + k);
// Stabilization parameters
if (mu < rho*U*h) {
delta1 = 0.5/sqrt(1.0/(h*h) + (U*U)/(h*h));
delta2 = h;
}
else
delta1 = delta2 = h*h;
delta1 *= detJW;
delta2 *= detJW;
// RUNAR
//delta1 = 0.0001*h*h*detJW;
//delta2 = 0.0001*h*h*detJW;
if (eM) {
Matrix& EM = *eM;
// Convection terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
// Sum convection for each direction
convJ = 0.0;
for (k = 1;k <= nsd;k++)
convJ += vel(k)*dNdX(i,k);
convJ *= rho*N(j)*detJW;
for (k = 1;k <= nsd;k++) {
EM((j-1)*nf+k,(i-1)*nf+k) += convJ;
for (l = 1;l <= nsd;l++)
EM((j-1)*nf+l,(i-1)*nf+k) += rho*dUdX(l,k)*N(j)*N(i)*detJW;
}
}
// Viscous terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
laplace = 0.0;
for (k = 1;k <= nsd;k++)
laplace += dNdX(i,k)*dNdX(j,k);
laplace *= mu*detJW;
for (k = 1;k <= nsd;k++)
EM(nf*(j-1)+k,nf*(i-1)+k) += laplace;
}
if (formulation == SIM::STRESS)
for (i = 1; i <= N.size(); i++)
for (j = 1; j <= N.size(); j++)
for (k = 1; k <= nsd; k++)
for (l = 1; l <= nsd; l++)
EM(nf*(j-1)+l,nf*(i-1)+k) += mu*dNdX(i,k)*dNdX(j,l)*detJW;
// Pressure and divergence terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++)
for (k = 1;k <= nsd;k++) {
divI = N(j)*dNdX(i,k)*detJW;
EM(nf*(i-1)+k,nf*j) -= divI;
EM(nf*j,nf*(i-1)+k) -= divI;
}
// Stabilization terms
// Computation of non-linear residual
res.fill(0.0);
for (k = 1;k <= nsd;k++) {
for (l = 1;l <= nsd;l++)
res(k) += rho*vel(l)*dUdX(k,l) - mu*d2UdX2(k,l,l);
res(k) += dPdX(k);
res(k) *= delta1;
}
// STABILIZATION TERMS
for (i = 1;i <= N.size();i++) {
// Convection for solution basis functions
convI = 0.0;
for (k = 1;k <= nsd;k++)
convI += vel(k)*dNdX(i,k);
convI *= rho;
for (j = 1;j <= N.size();j++) {
// Convection for test functions
convJ = 0.0;
for (k = 1;k <= nsd;k++)
convJ += vel(k)*dNdX(j,k);
convJ *= delta1*rho;
//--- MOMENTUM STABILIZATION ---
// Convection terms
for (k = 1;k <= nsd;k++) {
EM((j-1)*nf + k,(i-1)*nf + k) += convI*convJ;
for (l = 1;l <= nsd;l++) {
EM((j-1)*nf + l,(i-1)*nf + k) += rho*res(l)*dNdX(j,k)*N(i);
EM((j-1)*nf + l,(i-1)*nf + k) += convJ*rho*dUdX(l,k)*N(i);
}
}
// Pressure terms
for (l = 1;l <= nsd;l++)
EM((j-1)*nf + l,i*nf) += dNdX(i,l)*convJ;
// Viscous terms
diff = 0.0;
for (k = 1;k <= nsd;k++)
diff += d2NdX2(i,k,k);
diff *= mu;
for (k = 1;k <= nsd;k++)
EM((j-1)*nf + k,(i-1)*nf + k) -= diff*convJ;
// Continuity-continuity stabilization
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
EM((j-1)*nf+l,(i-1)*nf+k) += delta2*dNdX(i,k)*dNdX(j,l);
//--- CONTINUITY STABILIZATION ---
// Convection terms
for (k = 1;k <= nsd;k++)
EM(j*nf,(i-1)*nf+k) += delta1*convI*dNdX(j,k);
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
EM(j*nf,(i-1)*nf+k) += delta1*rho*dUdX(l,k)*N(i)*dNdX(j,l);
// Pressure terms
laplace = 0.0;
for (k = 1;k <= nsd;k++)
laplace += dNdX(i,k)*dNdX(j,k);
laplace *= delta1;
EM(j*nf,i*nf) += laplace;
// Viscous terms
diff = 0.0;
for (k = 1;k <= nsd;k++)
diff += d2NdX2(i,k,k);
diff *= delta1*mu;
for (k = 1;k <= nsd;k++)
EM(j*nf,(i-1)*nf+k) -= diff*dNdX(j,k);
}
}
}
if (eS) {
// Integrate rhs vector
Vector& ES = *eS;
Vector fb(g,nsd);
fb *= this->getMassDensity(X)*detJW;
for (i = 1;i <= N.size();i++) {
// --- Residual of Galerkin formulation ---
// -- Momentum equations --
// Body forces
for (k = 1;k <= nsd;k++)
ES((i-1)*nf+k) += fb(k)*N(i);
// Convection terms
for (k = 1;k <= nsd;k++) {
conv.fill(0.0);
for (m = 1;m <= nsd;m++)
conv(k) += vel(m)*dUdX(k,m);
ES((i-1)*nf + k) -= rho*conv(k)*N(i)*detJW;
}
// Pressure terms
for (k = 1;k <= nsd;k++)
ES((i-1)*nf + k) += P*dNdX(i,k)*detJW;
// Viscous terms
for (k = 1;k <= nsd;k++)
for (m = 1;m <= nsd;m++)
ES((i-1)*nf + k) -= mu*dUdX(k,m)*dNdX(i,m)*detJW;
// -- Continuity equation --
divU = 0.0;
for (m = 1;m <= nsd;m++)
divU += dUdX(m,m);
ES(i*nf) += divU*N(i)*detJW;
//--- MOMENTUM STABILIZATION---
// Convection for solution basis functions
convI = 0.0;
for (k = 1;k <= nsd;k++)
convI += vel(k)*dNdX(i,k);
convI *= rho;
// Convection terms
for (k = 1;k <= nsd;k++)
ES((i-1)*nf+k) -= res(k)*convI;
// Divergence of velocity field and velocity test function
divU = 0.0;
for (k = 1;k <= nsd;k++)
divU += dUdX(k,k);
// Continuity residual term
for (k = 1;k <= nsd;k++)
ES((i-1)*nf+k) -= delta2*divU*dNdX(i,k);
//--- CONTINUITY STABILIZATION ---
// Momentum residual term
for (k = 1;k <= nsd;k++)
ES(i*nf) -= res(k)*dNdX(i,k);
}
}
return getIntegralResult(elmInt);
}

137
src/Integrands/NavierStokesG2GenTheta.h
View File

@@ -1,137 +0,0 @@
//==============================================================================
//!
//! \file NavierStokesG2GenTheta.h
//!
//! \date August 3 2010
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Integrand implementations for Navier-Stokes equations with G2 stabilization
//! and generalized theta method for time integration
//!
//==============================================================================
#ifndef _NAVIER_STOKES_G2_GEN_THETA_H
#define _NAVIER_STOKES_G2_GEN_THETA_H
#include "NavierStokesG2.h"
/*!
\brief Class representing the integrand of the pressure stabilized Navier-Stokes problem.
\details This class supports G2-stabilized Navier-Stokes formulations using equal
order elements for velocity and pressure. Time integration by the generalized theta method.
*/
class NavierStokesG2GenTheta : public NavierStokesG2
{
public:
//! \brief The default constructor initializes all pointers to zero
//! \param[in] n Number of spatial dimensions
NavierStokesG2GenTheta(short int n, SIM::Formulation form, int itg = 3);
//! \brief The destructor clears the static work arrays used internally.
virtual ~NavierStokesG2GenTheta();
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] time Parameters for nonlinear and time-dependent simulations
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
//!
//! \details This interface is used when \a getIntegrandType returns 1.
virtual bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const
{ return false; }
//! \brief Evaluates the integrand at an interior 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] d2NdX2 Basis function second derivatives
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] h Characteristic element length
//!
//! \details This interface is used when \a getIntegrandType returns 2.
//! The default implementation forwards to the stationary version.
//! Reimplement this method for time-dependent or non-linear problems.
virtual bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Matrix3D& d2NdX2, const Vec3& X,
double h = 0.0) const;
//! \brief Evaluates the integrand at an interior 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] Navg Volume-averaged basis function values over the element
//! \param[in] X Cartesian coordinates of current integration point
//!
//! \details This interface is used when \a getIntegrandType returns 3.
//! The default implementation forwards to the stationary version.
//! Reimplement this method for time-dependent or non-linear problems.
virtual bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Vector& Navg, const Vec3& X) const
{ return true; }
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const
{ return false; }
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] d2NdX2 Basis function second derivatives
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] h characteristic element length
virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Matrix3D& d2NdX2,
const Vec3& X, double h = 0.0) const;
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] Navg Average value of basis function over element
//! \param[in] X Cartesian coordinates of current integration point
virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vector& Navg, const Vec3& X) const
{ return false; }
private:
int step; // Current step (1, 2 or 3)
real dt; // Current (sub-)timestep size
//static const real theta = 0.292893219; // Time integration parameter
//static const real thetam = 0.414213562; // Time integration parameter
//static const real thetat = 0.171572875; // Time integration parameter
//static const real alpha = 0.585786438; // Time integration parameter
//static const real beta = 0.414213562; // Time integration parameter
static const real theta = 0.5;
static const real thetam = 0.5;
static const real thetat = 0.5;
static const real alpha = 1.0;
static const real beta = 1.0;
};
#endif

651
src/Integrands/NavierStokesG2MP.C
View File

@@ -1,651 +0,0 @@
// $Id: NavierStokesG2MP.C,v 1.3 2011-02-08 12:41:32 rho Exp $
//==============================================================================
//!
//! \file NavierStokesG2MP.C
//!
//! \date August 11 2010
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Integrand implementations for G2 stabilized Navier-Stokes problems.
//! \details Time integration by the Crank-Nicholson method.
//!
//==============================================================================
#include "NavierStokesG2MP.h"
#include "TimeDomain.h"
#include "ElmNorm.h"
NavierStokesG2MP::NavierStokesG2MP(short int n, SIM::Formulation form, int itg)
: NavierStokesG2(n,form,itg) {}
NavierStokesG2MP::~NavierStokesG2MP() {}
bool NavierStokesG2MP::evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Matrix3D& d2NdX2, const Vec3& X,
double h) const
{
int i, j;
int k, l, m;
real temp1, temp2;
real diff, divU, divI, convI, convJ, laplace, U, P;
real Pb, Pn;
real delta1, delta2;
real dtInv, dUdt;
Vector vel(3), bvel(3), conv(3), res(3), dPdX(3), dPdXn(3), dPdXb(3);
Vector pvel(3);
Matrix dUdX(3,3), dUdXn(3,3), dUdXb(3,3);
Matrix3D d2UdX2(3,3,3), d2UdX2n(3,3,3), d2UdX2b(3,3,3);
const int nf = nsd + 1;
// Inverse of timestep size
if (time.dt > 0.0)
dtInv = 1.0/time.dt;
else
dtInv = 0.0;
Vector& EV = *eVs[0];
// Velocity in integration point
vel.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
vel(k) += EV((i-1)*nf + k)*N(i);
U = 0.0;
for (k = 1;k <= nsd;k++)
U += vel(k)*vel(k);
U = sqrt(U);
Vector& EVP = *eVs[1];
// Previous velocity at integration point
pvel.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
pvel(k) += EVP((i-1)*nf + k)*N(i);
// Intermediate velocity at integration point
bvel.fill(0.0);
for (k = 1;k <= nsd;k++)
bvel(k) = 0.5*(vel(k)+pvel(k));
// Pressure in integration point
P = 0.0;
for (i = 1;i <= N.size();i++)
P += EV(i*nf)*N(i);
// Pressure in integration point at previous timestep
Pn = 0.0;
for (i = 1;i <= N.size();i++)
Pn += EVP(i*nf)*N(i);
// Pressure in integration point at intermediate timestep
Pb = 0.5*(P+Pn);
// Velocity gradient in integration point
dUdX.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
dUdX(k,l) += dNdX(i,l)*EV((i-1)*nf+k);
// Velocity gradient at previous timestep
dUdXn.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
dUdXn(k,l) += dNdX(i,l)*EVP((i-1)*nf+k);
// Velocity gradient at intermediate timestep
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
dUdXb(k,l) = 0.5*(dUdX(k,l)+dUdXn(k,l));
// Pressure gradient in integration point
dPdX.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
dPdX(k) += dNdX(i,k)*EV(i*nf);
// Pressure gradient at previous timestep
dPdXn.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
dPdXn(k) += dNdX(i,k)*EVP(i*nf);
// Pressure gradient at intermediate timestep
for (k = 1;k <= nsd;k++)
dPdXb(k) = 0.5*(dPdX(k) + dPdXn(k));
// Hessian for velocity
d2UdX2.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
for (m = 1;m <= nsd;m++)
d2UdX2(k,m,l) += d2NdX2(i,m,l)*EV((i-1)*nf + k);
// Hessian for velocity at previous timestep
d2UdX2n.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
for (m = 1;m <= nsd;m++)
d2UdX2n(k,m,l) += d2NdX2(i,m,l)*EVP((i-1)*nf + k);
// Hessian for velocity at intermediate timestep
for (k = 1;k <= nsd;k++)
for (m = 1;m <= nsd;m++)
for (l = 1;l <= nsd;l++)
d2UdX2b(k,m,l) = 0.5*(d2UdX2(k,m,l) + d2UdX2n(k,m,l));
// Stabilization parameters
if (mu < rho*U*h) {
if (time.dt == 0.0)
delta1 = 0.5/sqrt(1.0/(h*h) + (U*U)/(h*h));
else
delta1 = 0.5/sqrt(dtInv*dtInv + (U*U)/(h*h));
delta2 = h;
}
else
delta1 = delta2 = h*h;
// RUNAR
//delta1 *= 0.01*detJW;
//delta2 *= 0.01*detJW;
delta1 *= 0.001*detJW;
delta2 *= 0.001*detJW;
if (eM) {
Matrix& EM = *eM;
// Time derivative term lhs
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
dUdt = rho*N(i)*N(j)*dtInv*detJW;
for (k = 1;k <= nsd;k++)
EM((j-1)*nf + k, (i-1)*nf + k) += dUdt;
}
// Convection terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
// Sum convection for each direction
convJ = 0.0;
for (k = 1;k <= nsd;k++)
convJ += bvel(k)*dNdX(i,k);
convJ *= 0.5*rho*N(j)*detJW;
temp1 = 0.5*rho*N(i)*N(j)*detJW;
for (k = 1;k <= nsd;k++) {
EM((j-1)*nf+k,(i-1)*nf+k) += convJ;
for (l = 1;l <= nsd;l++)
EM((j-1)*nf+l,(i-1)*nf+k) += temp1*dUdXb(l,k);
}
}
// Viscous terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
laplace = 0.0;
for (k = 1;k <= nsd;k++)
laplace += dNdX(i,k)*dNdX(j,k);
laplace *= 0.5*mu*detJW;
for (k = 1;k <= nsd;k++)
EM(nf*(j-1)+k,nf*(i-1)+k) += laplace;
}
temp1 = 0.5*mu*detJW;
if (formulation == SIM::STRESS)
for (i = 1; i <= N.size(); i++)
for (j = 1; j <= N.size(); j++)
for (k = 1; k <= nsd; k++)
for (l = 1; l <= nsd; l++)
EM(nf*(j-1)+l,nf*(i-1)+k) += temp1*dNdX(i,k)*dNdX(j,l);
// Pressure and divergence terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++)
for (k = 1;k <= nsd;k++) {
divI = N(j)*dNdX(i,k)*detJW;
EM(nf*(i-1)+k,nf*j) -= 0.5*divI;
EM(nf*j,nf*(i-1)+k) -= 0.5*divI;
}
// Stabilization terms
// Computation of non-linear residual
res.fill(0.0);
for (k = 1;k <= nsd;k++) {
for (l = 1;l <= nsd;l++)
res(k) += rho*bvel(l)*dUdXb(k,l) - mu*d2UdX2b(k,l,l);
res(k) += dPdXb(k);
res(k) *= delta1;
}
// STABILIZATION TERMS
for (i = 1;i <= N.size();i++) {
// Convection for solution basis functions
convI = 0.0;
for (k = 1;k <= nsd;k++)
convI += bvel(k)*dNdX(i,k);
convI *= rho;
for (j = 1;j <= N.size();j++) {
// Convection for test functions
convJ = 0.0;
for (k = 1;k <= nsd;k++)
convJ += bvel(k)*dNdX(j,k);
convJ *= delta1*rho;
//--- MOMENTUM STABILIZATION ---
// Convection terms
for (k = 1;k <= nsd;k++) {
EM((j-1)*nf + k,(i-1)*nf + k) += 0.5*convI*convJ;
for (l = 1;l <= nsd;l++) {
EM((j-1)*nf + l,(i-1)*nf + k) += 0.5*rho*res(l)*dNdX(j,k)*N(i);
EM((j-1)*nf + l,(i-1)*nf + k) += 0.5*convJ*rho*dUdXb(l,k)*N(i);
}
}
// Pressure terms
for (l = 1;l <= nsd;l++)
EM((j-1)*nf + l,i*nf) += 0.5*dNdX(i,l)*convJ;
// Viscous terms
diff = 0.0;
for (k = 1;k <= nsd;k++)
diff += d2NdX2(i,k,k);
diff *= 0.5*mu;
for (k = 1;k <= nsd;k++)
EM((j-1)*nf + k,(i-1)*nf + k) -= diff*convJ;
// Continuity-continuity stabilization
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
EM((j-1)*nf+l,(i-1)*nf+k) -= 0.5*delta2*dNdX(i,k)*dNdX(j,l);
//--- CONTINUITY STABILIZATION ---
// Convection terms
for (k = 1;k <= nsd;k++)
EM(j*nf,(i-1)*nf+k) += 0.5*delta1*convI*dNdX(j,k);
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
EM(j*nf,(i-1)*nf+k) += 0.5*delta1*rho*dUdXb(l,k)*N(i)*dNdX(j,l);
// Pressure terms
laplace = 0.0;
for (k = 1;k <= nsd;k++)
laplace += dNdX(i,k)*dNdX(j,k);
laplace *= 0.5*delta1;
EM(j*nf,i*nf) += laplace;
// Viscous terms
diff = 0.0;
for (k = 1;k <= nsd;k++)
diff += d2NdX2(i,k,k);
diff *= 0.5*delta1*mu;
for (k = 1;k <= nsd;k++)
EM(j*nf,(i-1)*nf+k) -= diff*dNdX(j,k);
}
}
}
if (eS) {
// Integrate rhs vector
Vector& ES = *eS;
Vector fb(g,nsd);
fb *= this->getMassDensity(X)*detJW;
for (i = 1;i <= N.size();i++) {
// --- Residual of Galerkin formulation ---
// -- Momentum equations --
// Time derivative terms
for (k = 1;k <= nsd;k++)
ES((i-1)*nf + k) += rho*(pvel(k)-vel(k))*N(i)*dtInv*detJW;
// Body forces
for (k = 1;k <= nsd;k++)
ES((i-1)*nf+k) += fb(k)*N(i);
// Convection terms
for (k = 1;k <= nsd;k++) {
conv.fill(0.0);
for (m = 1;m <= nsd;m++)
conv(k) += bvel(m)*dUdXb(k,m);
ES((i-1)*nf + k) -= rho*conv(k)*N(i)*detJW;
}
// Pressure terms
for (k = 1;k <= nsd;k++)
ES((i-1)*nf + k) += Pb*dNdX(i,k)*detJW;
// Viscous terms
for (k = 1;k <= nsd;k++)
for (m = 1;m <= nsd;m++)
ES((i-1)*nf + k) -= mu*dUdXb(k,m)*dNdX(i,m)*detJW;
// -- Continuity equation --
divU = 0.0;
for (m = 1;m <= nsd;m++)
divU += dUdXb(m,m);
ES(i*nf) += divU*N(i)*detJW;
//--- MOMENTUM STABILIZATION---
// Convection for solution basis functions
convI = 0.0;
for (k = 1;k <= nsd;k++)
convI += vel(k)*dNdX(i,k);
convI *= rho;
// Convection terms
for (k = 1;k <= nsd;k++)
ES((i-1)*nf+k) -= res(k)*convI;
// Continuity residual term
for (k = 1;k <= nsd;k++)
ES((i-1)*nf+k) += delta2*divU*dNdX(i,k);
//--- CONTINUITY STABILIZATION ---
// Momentum residual term
for (k = 1;k <= nsd;k++)
ES(i*nf) -= res(k)*dNdX(i,k);
}
}
return getIntegralResult(elmInt);
}
bool NavierStokesG2MP::evalInt (LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Matrix3D& d2NdX2,
const Vec3& X, double h) const
{
int i, j;
int k, l, m;
real diff, divU, divI, divJ, convI, convJ, laplace, U, P;
real delta1, delta2;
Vector vel(3), conv(3), res(3), dPdX(3);
Matrix dUdX(3,3);
Matrix3D d2UdX2(3,3,3);
const int nf = nsd + 1;
Vector& EV = *eVs[0];
// Velocity in integration point
vel.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
vel(k) += EV((i-1)*nf + k)*N(i);
U = 0.0;
for (k = 1;k <= nsd;k++)
U += vel(k)*vel(k);
U = sqrt(U);
// Pressure in integration point
P = 0.0;
for (i = 1;i <= N.size();i++)
P += EV(i*nf)*N(i);
// Velocity gradient in integration point
dUdX.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
dUdX(k,l) += dNdX(i,l)*EV((i-1)*nf+k);
// Pressure gradient in integration point
dPdX.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
dPdX(k) += dNdX(i,k)*EV(i*nf);
// Hessian for velocity
d2UdX2.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
for (m = 1;m <= nsd;m++)
d2UdX2(k,m,l) += d2NdX2(i,m,l)*EV((i-1)*nf + k);
// Stabilization parameters
if (mu < rho*U*h) {
delta1 = 0.5/sqrt(1.0/(h*h) + (U*U)/(h*h));
delta2 = h;
}
else
delta1 = delta2 = h*h;
delta1 *= detJW;
delta2 *= detJW;
// RUNAR
//delta1 = 0.0001*h*h*detJW;
//delta2 = 0.0001*h*h*detJW;
if (eM) {
Matrix& EM = *eM;
// Convection terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
// Sum convection for each direction
convJ = 0.0;
for (k = 1;k <= nsd;k++)
convJ += vel(k)*dNdX(i,k);
convJ *= rho*N(j)*detJW;
for (k = 1;k <= nsd;k++) {
EM((j-1)*nf+k,(i-1)*nf+k) += convJ;
for (l = 1;l <= nsd;l++)
EM((j-1)*nf+l,(i-1)*nf+k) += rho*dUdX(l,k)*N(j)*N(i)*detJW;
}
}
// Viscous terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
laplace = 0.0;
for (k = 1;k <= nsd;k++)
laplace += dNdX(i,k)*dNdX(j,k);
laplace *= mu*detJW;
for (k = 1;k <= nsd;k++)
EM(nf*(j-1)+k,nf*(i-1)+k) += laplace;
}
if (formulation == SIM::STRESS)
for (i = 1; i <= N.size(); i++)
for (j = 1; j <= N.size(); j++)
for (k = 1; k <= nsd; k++)
for (l = 1; l <= nsd; l++)
EM(nf*(j-1)+l,nf*(i-1)+k) += mu*dNdX(i,k)*dNdX(j,l)*detJW;
// Pressure and divergence terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++)
for (k = 1;k <= nsd;k++) {
divI = N(j)*dNdX(i,k)*detJW;
EM(nf*(i-1)+k,nf*j) -= divI;
EM(nf*j,nf*(i-1)+k) -= divI;
}
// Stabilization terms
// Computation of non-linear residual
res.fill(0.0);
for (k = 1;k <= nsd;k++) {
for (l = 1;l <= nsd;l++)
res(k) += rho*vel(l)*dUdX(k,l) - mu*d2UdX2(k,l,l);
res(k) += dPdX(k);
res(k) *= delta1;
}
// STABILIZATION TERMS
for (i = 1;i <= N.size();i++) {
// Convection for solution basis functions
convI = 0.0;
for (k = 1;k <= nsd;k++)
convI += vel(k)*dNdX(i,k);
convI *= rho;
for (j = 1;j <= N.size();j++) {
// Convection for test functions
convJ = 0.0;
for (k = 1;k <= nsd;k++)
convJ += vel(k)*dNdX(j,k);
convJ *= delta1*rho;
//--- MOMENTUM STABILIZATION ---
// Convection terms
for (k = 1;k <= nsd;k++) {
EM((j-1)*nf + k,(i-1)*nf + k) += convI*convJ;
for (l = 1;l <= nsd;l++) {
// RUNAR
//EM((j-1)*nf + l,(i-1)*nf + k) += rho*res(l)*dNdX(j,k)*N(i);
EM((j-1)*nf + l,(i-1)*nf + k) += convJ*rho*dUdX(l,k)*N(i);
}
}
// Pressure terms
for (l = 1;l <= nsd;l++)
EM((j-1)*nf + l,i*nf) += dNdX(i,l)*convJ;
// Viscous terms
diff = 0.0;
for (k = 1;k <= nsd;k++)
diff += d2NdX2(i,k,k);
diff *= mu;
for (k = 1;k <= nsd;k++)
EM((j-1)*nf + k,(i-1)*nf + k) -= diff*convJ;
// Continuity-continuity stabilization
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
EM((j-1)*nf+l,(i-1)*nf+k) += delta2*dNdX(i,k)*dNdX(j,l);
//--- CONTINUITY STABILIZATION ---
// Convection terms
for (k = 1;k <= nsd;k++)
EM(j*nf,(i-1)*nf+k) += delta1*convI*dNdX(j,k);
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
EM(j*nf,(i-1)*nf+k) += delta1*rho*dUdX(l,k)*N(i)*dNdX(j,l);
// Pressure terms
laplace = 0.0;
for (k = 1;k <= nsd;k++)
laplace += dNdX(i,k)*dNdX(j,k);
laplace *= delta1;
EM(j*nf,i*nf) += laplace;
// Viscous terms
diff = 0.0;
for (k = 1;k <= nsd;k++)
diff += d2NdX2(i,k,k);
diff *= delta1*mu;
for (k = 1;k <= nsd;k++)
EM(j*nf,(i-1)*nf+k) -= diff*dNdX(j,k);
}
}
}
if (eS) {
// Integrate rhs vector
Vector& ES = *eS;
Vector fb(g,nsd);
fb *= this->getMassDensity(X)*detJW;
for (i = 1;i <= N.size();i++) {
// --- Residual of Galerkin formulation ---
// -- Momentum equations --
// Body forces
for (k = 1;k <= nsd;k++)
ES((i-1)*nf+k) += fb(k)*N(i);
// Convection terms
for (k = 1;k <= nsd;k++) {
conv.fill(0.0);
for (m = 1;m <= nsd;m++)
conv(k) += vel(m)*dUdX(k,m);
ES((i-1)*nf + k) -= rho*conv(k)*N(i)*detJW;
}
// Pressure terms
for (k = 1;k <= nsd;k++)
ES((i-1)*nf + k) += P*dNdX(i,k)*detJW;
// Viscous terms
for (k = 1;k <= nsd;k++)
for (m = 1;m <= nsd;m++)
ES((i-1)*nf + k) -= mu*dUdX(k,m)*dNdX(i,m)*detJW;
// -- Continuity equation --
divU = 0.0;
for (m = 1;m <= nsd;m++)
divU += dUdX(m,m);
ES(i*nf) += divU*N(i)*detJW;
//--- MOMENTUM STABILIZATION---
// Convection for solution basis functions
convI = 0.0;
for (k = 1;k <= nsd;k++)
convI += vel(k)*dNdX(i,k);
convI *= rho;
// Convection terms
for (k = 1;k <= nsd;k++)
ES((i-1)*nf+k) -= res(k)*convI;
// Divergence of velocity field and velocity test function
divU = 0.0;
for (k = 1;k <= nsd;k++)
divU += dUdX(k,k);
// Continuity residual term
for (k = 1;k <= nsd;k++)
ES((i-1)*nf+k) -= delta2*divU*dNdX(i,k);
//--- CONTINUITY STABILIZATION ---
// Momentum residual term
for (k = 1;k <= nsd;k++)
ES(i*nf) -= res(k)*dNdX(i,k);
}
}
return getIntegralResult(elmInt);
}

123
src/Integrands/NavierStokesG2MP.h
View File

@@ -1,123 +0,0 @@
//==============================================================================
//!
//! \file NavierStokesG2MP.h
//!
//! \date August 11 2010
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Integrand implementations for Navier-Stokes equations with G2 stabilization
//! and. Time integration by mid-point method
//!
//==============================================================================
#ifndef _NAVIER_STOKES_G2_MP_H
#define _NAVIER_STOKES_G2_MP_H
#include "NavierStokesG2.h"
/*!
\brief Class representing the integrand of the pressure stabilized Navier-Stokes problem.
\details This class supports G2-stabilized Navier-Stokes formulations using equal
order elements for velocity and pressure. Time integration by the mid-point method.
*/
class NavierStokesG2MP : public NavierStokesG2
{
public:
//! \brief The default constructor initializes all pointers to zero
//! \param[in] n Number of spatial dimensions
NavierStokesG2MP(short int n, SIM::Formulation form, int itg = 3);
//! \brief The destructor clears the static work arrays used internally.
virtual ~NavierStokesG2MP();
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] time Parameters for nonlinear and time-dependent simulations
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
//!
//! \details This interface is used when \a getIntegrandType returns 1.
virtual bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const
{ return false; }
//! \brief Evaluates the integrand at an interior 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] d2NdX2 Basis function second derivatives
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] h Characteristic element length
//!
//! \details This interface is used when \a getIntegrandType returns 2.
//! The default implementation forwards to the stationary version.
//! Reimplement this method for time-dependent or non-linear problems.
virtual bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Matrix3D& d2NdX2, const Vec3& X,
double h = 0.0) const;
//! \brief Evaluates the integrand at an interior 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] Navg Volume-averaged basis function values over the element
//! \param[in] X Cartesian coordinates of current integration point
//!
//! \details This interface is used when \a getIntegrandType returns 3.
//! The default implementation forwards to the stationary version.
//! Reimplement this method for time-dependent or non-linear problems.
virtual bool evalInt(LocalIntegral*& elmInt,
const TimeDomain& time, double detJW,
const Vector& N, const Matrix& dNdX,
const Vector& Navg, const Vec3& X) const
{ return true; }
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const
{ return false; }
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] d2NdX2 Basis function second derivatives
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] h characteristic element length
virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Matrix3D& d2NdX2,
const Vec3& X, double h = 0.0) const;
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] Navg Average value of basis function over element
//! \param[in] X Cartesian coordinates of current integration point
virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vector& Navg, const Vec3& X) const
{ return false; }
};
#endif

322
src/Integrands/StabilizedNavierStokes.C
View File

@@ -1,322 +0,0 @@
// $Id: StabilizedNavierStokes.C,v 1.5 2011-02-08 12:43:49 rho Exp $
//==============================================================================
//!
//! \file StabilizedNavierStokes.C
//!
//! \date Mar 17 2010
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Integrand implementations for stabilized Navier-Stokes problems.
//!
//==============================================================================
#include "StabilizedNavierStokes.h"
#include "Utilities.h"
#include "Vec3Oper.h"
StabilizedNavierStokes::StabilizedNavierStokes(short int n, SIM::Formulation form, int itg)
: StabilizedStokes(n,form,itg)
{
mu = 1.0e-3;
rho = 1.0;
}
StabilizedNavierStokes::~StabilizedNavierStokes()
{
}
bool StabilizedNavierStokes::evalInt(LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const
{
int i, j;
int k, l;
double div, laplace, conv;
Vector vel(3);
if (eM) {
Matrix& EM = *eM;
Vector& EV = *eVs[0];
// Velocity in integration point
vel.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
vel(k) += EV((i-1)*nf + k)*N(i);
// Convection terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
// Sum convection for each direction
conv = 0.0;
for (k = 1;k <= nsd;k++)
conv += vel(k)*dNdX(j,k);
conv *= rho*N(i)*detJW;
for (k = 1;k <= nsd;k++)
EM((i-1)*nf+k,(j-1)*nf+k) += conv;
}
// Viscous terms
for (i = 1; i <= N.size(); i++)
for (j = 1; j <= N.size(); j++) {
laplace = 0.0;
for (k = 1; k <= nsd; k++)
laplace += dNdX(i,k)*dNdX(j,k);
laplace *= mu*detJW;
for (k = 1; k <= nsd; k++)
EM(nf*(j-1)+k,nf*(i-1)+k) += laplace;
}
if (formulation == SIM::STRESS)
for (i = 1; i <= N.size(); i++)
for (j = 1; j <= N.size(); j++)
for (k = 1; k <= nsd; k++)
for (l = 1; l <= nsd; l++)
EM(nf*(j-1)+k,nf*(i-1)+l) += mu*dNdX(i,k)*dNdX(j,l)*detJW;
// Pressure and divergence terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++)
for (k = 1;k <= nsd;k++) {
div = N(j)*dNdX(i,k)*detJW;
EM(nf*(i-1)+k,nf*j) -= div;
EM(nf*j,nf*(i-1)+k) -= div;
}
// Pressure stabilization
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++)
EM(nf*j,nf*i) -= N(i)*N(j)*detJW;
}
if (eS) {
// Integrate body force vector
Vector& ES = *eS;
Vector fb(g,nsd);
fb *= this->getMassDensity(X)*detJW;
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
ES((i-1)*nf+k) += fb(k)*N(i);
}
return getIntegralResult(elmInt);
}
bool StabilizedNavierStokes::evalInt (LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Matrix3D& d2NdX2,
const Vec3& X, double h) const
{
int i, j;
int k, l;
double div, laplace, conv;
Vector vel(3);
const real delta = 0.001*h*h*detJW;
if (eM) {
Matrix& EM = *eM;
Vector& EV = *eVs[0];
// Velocity in integration point
vel.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
vel(k) += EV((i-1)*nf + k)*N(i);
// Convection terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
// Sum convection for each direction
conv = 0.0;
for (k = 1;k <= nsd;k++)
conv += vel(k)*dNdX(j,k);
conv *= rho*N(i)*detJW;
for (k = 1;k <= nsd;k++)
EM((i-1)*nf+k,(j-1)*nf+k) += conv;
}
// Viscous terms
for (i = 1; i <= N.size(); i++)
for (j = 1; j <= N.size(); j++) {
laplace = 0.0;
for (k = 1; k <= nsd; k++)
laplace += dNdX(i,k)*dNdX(j,k);
laplace *= mu*detJW;
for (k = 1; k <= nsd; k++)
EM(nf*(j-1)+k,nf*(i-1)+k) += laplace;
}
if (formulation == SIM::STRESS)
for (i = 1; i <= N.size(); i++)
for (j = 1; j <= N.size(); j++)
for (k = 1; k <= nsd; k++)
for (l = 1; l <= nsd; l++)
EM(nf*(j-1)+k,nf*(i-1)+l) += mu*dNdX(i,k)*dNdX(j,l)*detJW;
// Pressure and divergence terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++)
for (k = 1;k <= nsd;k++) {
div = N(j)*dNdX(i,k)*detJW;
EM(nf*(i-1)+k,nf*j) -= div;
EM(nf*j,nf*(i-1)+k) -= div;
}
// Pressure stabilization
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
for (k = 1;k <= nsd;k++) {
laplace = 0.0;
for (l = 1;l <= nsd;l++)
laplace += d2NdX2(i,l,l);
laplace *= delta*mu;
EM(nf*j,nf*(i-1) + k) += laplace*dNdX(j,k);
EM(nf*j,nf*i) -= delta*dNdX(i,k)*dNdX(j,k);
}
}
if (formulation == SIM::STRESS)
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
if (k != l)
EM(nf*j,nf*(i-1)+l) += delta*mu*d2NdX2(i,k,l)*dNdX(j,k);
}
if (eS) {
// Integrate body force vector
Vector& ES = *eS;
Vector fb(g,nsd);
fb *= this->getMassDensity(X)*detJW;
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
ES((i-1)*nf+k) += fb(k)*N(i);
// Pressure stabilization term
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
ES(nf*i) -= delta*fb(k)*dNdX(i,k);
}
return getIntegralResult(elmInt);
}
bool StabilizedNavierStokes::evalInt(LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vector& Navg, const Vec3& X) const
{
int i, j;
int k, l;
double div, laplace, conv;
Vector vel(3);
if (eM) {
Matrix& EM = *eM;
Vector& EV = *eVs[0];
// Velocity in integration point
vel.fill(0.0);
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
vel(k) += EV((i-1)*nf + k)*N(i);
//Convection terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
// Sum convection for each direction
conv = 0.0;
for (k = 1;k <= nsd;k++)
conv += vel(k)*dNdX(j,k);
conv *= rho*detJW;
for (k = 1;k <= nsd;k++)
EM((i-1)*nf+k,(j-1)*nf+k) += conv*N(i);
}
// Viscous terms
for (i = 1; i <= N.size(); i++)
for (j = 1; j <= N.size(); j++) {
laplace = 0.0;
for (k = 1; k <= nsd; k++)
laplace += dNdX(i,k)*dNdX(j,k);
laplace *= mu*detJW;
for (k = 1; k <= nsd; k++)
EM(nf*(j-1)+k,nf*(i-1)+k) += laplace;
}
if (formulation == SIM::STRESS)
for (i = 1; i <= N.size(); i++)
for (j = 1; j <= N.size(); j++)
for (k = 1; k <= nsd; k++)
for (l = 1; l <= nsd; l++)
EM(nf*(j-1)+k,nf*(i-1)+l) += mu*dNdX(i,k)*dNdX(j,l)*detJW;
// Pressure and divergence terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++)
for (k = 1;k <= nsd;k++) {
div = N(j)*dNdX(i,k)*detJW;
EM(nf*(i-1)+k,nf*j) -= div;
EM(nf*j,nf*(i-1)+k) -= div;
}
// Pressure stabilization
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++)
EM(nf*j,nf*i) -= (N(i)-Navg(i))*(N(j)-Navg(j))*detJW;
}
if (eS) {
// Integrate body force vector
Vector& ES = *eS;
Vector fb(g,nsd);
fb *= this->getMassDensity(X)*detJW;
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
ES((i-1)*nf+k) += fb(k)*N(i);
}
return getIntegralResult(elmInt);
}
bool StabilizedNavierStokes::evalBou (LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X, const Vec3& normal) const
{
if (!eS || !tracFld)
{
std::cerr <<" *** StabilizedNavierStokes::evalBou: Zero pointers."<< std::endl;
return false;
}
// Evaluate the surface traction
Vec3 T = (*tracFld)(X,normal);
// Store the traction value for vizualization
if (!T.isZero())
tracVal.insert(std::make_pair(X,T));
// Integrate the force vector
Vector& ES = *eS;
for (size_t a = 1; a <= N.size(); a++)
for (short int d = 1; d <= nsd; d++)
ES(nf*(a-1)+d) += T[d-1]*N(a)*detJW;
return getIntegralResult(elmInt);
}

82
src/Integrands/StabilizedNavierStokes.h
View File

@@ -1,82 +0,0 @@
//==============================================================================
//!
//! \file StabilizedNavierStokes.h
//!
//! \date Mar 17 2010
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Integrand implementations for pressure stabilized Navier-Stokes problems.
//!
//==============================================================================
#ifndef _STABILIZED_NAVIER_STOKES_H
#define _STABILIZED_NAVIER_STOKES_H
#include "StabilizedStokes.h"
/*!
\brief Class representing the integrand of the pressure stabilized Navier-Stokes problem.
\details This class supports pressure stabilized Navier-Stokes formulations using equal
order elements for velocity and pressure.
*/
class StabilizedNavierStokes : public StabilizedStokes
{
public:
//! \brief The default constructor initializes all pointers to zero
//! \param[in] n Number of spatial dimensions
StabilizedNavierStokes(short int n, SIM::Formulation form, int itg = 3);
//! \brief The destructor clears the static work arrays used internally.
virtual ~StabilizedNavierStokes();
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const;
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] d2NdX2 Basis function second derivatives
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] h characteristic element length
virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Matrix3D& d2NdX2,
const Vec3& X, double h = 0.0) const;
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] Navg Average value of basis function over element
//! \param[in] X Cartesian coordinates of current integration point
virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vector& Navg, const Vec3& X) 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] 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
virtual bool evalBou(LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X, const Vec3& normal) const;
//! \brief Returns which integrand to use
int getIntegrandType() const { return 3; }
};
#endif

241
src/Integrands/StabilizedStokes.C
View File

@@ -1,241 +0,0 @@
// $Id: StabilizedStokes.C,v 1.7 2011-02-08 12:42:36 rho Exp $
//==============================================================================
//!
//! \file StabilizedStokes.C
//!
//! \date Feb 11 2010
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Integrand implementations for stabilized Stokes problems.
//!
//==============================================================================
#include "StabilizedStokes.h"
#include "Utilities.h"
#include "ElmNorm.h"
#include "Tensor.h"
#include "Vec3Oper.h"
#include "VTF.h"
StabilizedStokes::StabilizedStokes (short int n,
SIM::Formulation form, int itg)
: Stokes(n,form,itg)
{
// Resize element solution vectors
eVs.resize(2);
for (int i = 0;i < 2;i++)
eVs[i] = new Vector();
primsol.resize(2);
}
bool StabilizedStokes::evalInt (LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const
{
size_t i, j;
short int k, l;
double div, laplace;
const int nf = nsd + 1;
if (eM) {
Matrix& EM = *eM;
// Viscous terms
for (i = 1; i <= N.size(); i++)
for (j = 1; j <= N.size(); j++) {
laplace = 0.0;
for (k = 1; k <= nsd; k++)
laplace += dNdX(i,k)*dNdX(j,k);
laplace *= mu*detJW;
for (k = 1; k <= nsd; k++)
EM(nf*(j-1)+k,nf*(i-1)+k) += laplace;
}
if (formulation == SIM::STRESS)
for (i = 1; i <= N.size(); i++)
for (j = 1; j <= N.size(); j++)
for (k = 1; k <= nsd; k++)
for (l = 1; l <= nsd; l++)
EM(nf*(j-1)+k,nf*(i-1)+l) += mu*dNdX(i,k)*dNdX(j,l)*detJW;
// Pressure and divergence terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++)
for (k = 1;k <= nsd;k++) {
div = N(j)*dNdX(i,k)*detJW;
EM(nf*(i-1)+k,nf*j) -= div;
EM(nf*j,nf*(i-1)+k) -= div;
}
// Pressure stabilization
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++)
EM(nf*j,nf*i) -= N(i)*N(j)*detJW;
}
if (eS) {
// Integrate body force vector
Vector& ES = *eS;
Vector fb(g,nsd);
fb *= this->getMassDensity(X)*detJW;
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
ES((i-1)*nf+k) += fb(k)*N(i);
}
return getIntegralResult(elmInt);
}
bool StabilizedStokes::evalInt (LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Matrix3D& d2NdX2,
const Vec3& X, double h) const
{
size_t i, j;
short int k, l;
double div, laplace;
const int nf = nsd + 1;
const real delta = 0.001*h*h*detJW;
if (eM) {
Matrix& EM = *eM;
// Viscous terms
for (i = 1; i <= N.size(); i++)
for (j = 1; j <= N.size(); j++) {
laplace = 0.0;
for (k = 1; k <= nsd; k++)
laplace += dNdX(i,k)*dNdX(j,k);
laplace *= mu*detJW;
for (k = 1; k <= nsd; k++)
EM(nf*(j-1)+k,nf*(i-1)+k) += laplace;
}
if (formulation == SIM::STRESS)
for (i = 1; i <= N.size(); i++)
for (j = 1; j <= N.size(); j++)
for (k = 1; k <= nsd; k++)
for (l = 1; l <= nsd; l++)
EM(nf*(j-1)+k,nf*(i-1)+l) += mu*dNdX(i,k)*dNdX(j,l)*detJW;
// Pressure and divergence terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++)
for (k = 1;k <= nsd;k++) {
div = N(j)*dNdX(i,k)*detJW;
EM(nf*(i-1)+k,nf*j) -= div;
EM(nf*j,nf*(i-1)+k) -= div;
}
// Pressure stabilization
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++) {
for (k = 1;k <= nsd;k++) {
laplace = 0.0;
for (l = 1;l <= nsd;l++)
laplace += d2NdX2(i,l,l);
laplace *= delta*mu;
EM(nf*j,nf*(i-1) + k) += laplace*dNdX(j,k);
EM(nf*j,nf*i) -= delta*dNdX(i,k)*dNdX(j,k);
}
}
if (formulation == SIM::STRESS)
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++)
for (k = 1;k <= nsd;k++)
for (l = 1;l <= nsd;l++)
if (k != l)
EM(nf*j,nf*(i-1)+l) += delta*mu*d2NdX2(i,k,l)*dNdX(j,k);
}
if (eS) {
// Integrate body force vector
Vector& ES = *eS;
Vector fb(g,nsd);
fb *= this->getMassDensity(X)*detJW;
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
ES((i-1)*nf+k) += fb(k)*N(i);
// Pressure stabilization term
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
ES(nf*i) -= delta*fb(k)*dNdX(i,k);
}
return getIntegralResult(elmInt);
}
bool StabilizedStokes::evalInt (LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vector& Navg, const Vec3& X) const
{
size_t i, j;
short int k, l;
double div, laplace;
const int nf = nsd + 1;
if (eM) {
Matrix& EM = *eM;
// Viscous terms
for (i = 1; i <= N.size(); i++)
for (j = 1; j <= N.size(); j++) {
laplace = 0.0;
for (k = 1; k <= nsd; k++)
laplace += dNdX(i,k)*dNdX(j,k);
laplace *= mu*detJW;
for (k = 1; k <= nsd; k++)
EM(nf*(j-1)+k,nf*(i-1)+k) += laplace;
}
if (formulation == SIM::STRESS)
for (i = 1; i <= N.size(); i++)
for (j = 1; j <= N.size(); j++)
for (k = 1; k <= nsd; k++)
for (l = 1; l <= nsd; l++)
EM(nf*(j-1)+k,nf*(i-1)+l) += mu*dNdX(i,k)*dNdX(j,l)*detJW;
// Pressure and divergence terms
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++)
for (k = 1;k <= nsd;k++) {
div = N(j)*dNdX(i,k)*detJW;
EM(nf*(i-1)+k,nf*j) -= div;
EM(nf*j,nf*(i-1)+k) -= div;
}
// Pressure stabilization
for (i = 1;i <= N.size();i++)
for (j = 1;j <= N.size();j++)
EM(nf*j,nf*i) -= (N(i)-Navg(i))*(N(j)-Navg(j))*detJW;
}
if (eS) {
// Integrate body force vector
Vector& ES = *eS;
Vector fb(g,nsd);
fb *= this->getMassDensity(X)*detJW;
for (i = 1;i <= N.size();i++)
for (k = 1;k <= nsd;k++)
ES((i-1)*nf+k) += fb(k)*N(i);
}
return getIntegralResult(elmInt);
}

72
src/Integrands/StabilizedStokes.h
View File

@@ -1,72 +0,0 @@
// $Id: StabilizedStokes.h,v 1.4 2010-12-06 08:49:59 rho Exp $
//==============================================================================
//!
//! \file StabilizedStokes.h
//!
//! \date Feb 11 2010
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Integrand implementations for pressure stabilized Stokes problems.
//!
//==============================================================================
#ifndef _STABILIZED_STOKES_H
#define _STABILIZED_STOKES_H
#include "Stokes.h"
/*!
\brief Class representing the integrand of the pressure-stabilized
Stokes problem.
\details This class supports pressure stabilized Stokes formulations
using equal order elements for velocity and pressure.
*/
class StabilizedStokes : public Stokes
{
public:
//! \brief The default constructor initializes all pointers to zero.
//! \param[in] n Number of spatial dimensions
//! \param[in] form The solution formulation to use
//! \param[in] itg The integrandtype to use
StabilizedStokes(short int n, SIM::Formulation form = SIM::LAPLACE,
int itg = 3);
//! \brief The destructor frees dynamically allocated objects.
virtual ~StabilizedStokes() {}
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const;
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] d2NdX2 Basis function second derivatives
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] h Characteristic element length
virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Matrix3D& d2NdX2,
const Vec3& X, double h = 0.0) const;
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] Navg Volume-averaged basis function values over the element
//! \param[in] X Cartesian coordinates of current integration point
virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vector& Navg, const Vec3& X) const;
};
#endif

487
src/Integrands/Stokes.C
View File

@@ -1,487 +0,0 @@
// $Id: Stokes.C,v 1.3 2011-02-08 12:23:54 rho Exp $
//==============================================================================
//!
//! \file Stokes.C
//!
//! \date Sep 29 2010
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Integrand implementations for Stokes problems.
//!
//==============================================================================
#include "Stokes.h"
#include "Utilities.h"
#include "ElmNorm.h"
#include "Tensor.h"
#include "Vec3Oper.h"
#include "VTF.h"
Stokes::Stokes (unsigned short int n, SIM::Formulation form, int itg_type)
{
nsd = n;
nf = n+1;
formulation = form;
integrandType = itg_type;
mu = rho = 1.0;
tracFld = 0;
myMats = new ElmMats();
myMats->resize(1,1);
myMats->rhsOnly = false;
eM = &myMats->A[0];
eS = &myMats->b[0];
}
Stokes::~Stokes ()
{
if (myMats) delete myMats;
for (int i = 0; i < 2; i++)
if (eVs[i]) delete eVs[i];
}
bool Stokes::initElement (const std::vector<int>& MNPC)
{
const size_t nen = MNPC.size();
eM->resize(nf*nen,nf*nen,true);
eS->resize(nf*nen,true);
int ierr = 0;
for (size_t i = 0; i < 2 && i < primsol.size() && eVs[i]; i++)
if (ierr = utl::gather(MNPC,nf,primsol[i],*eVs[i]))
std::cerr <<" *** Stokes::initElement: Detected "
<< ierr <<" node numbers out of range."<< std::endl;
myMats->withLHS = true;
return ierr == 0;
}
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 ierr == 0;
}
bool Stokes::evalBou (LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X, const Vec3& normal) const
{
if (!tracFld)
{
std::cerr <<" *** Stokes::evalBou: No tractions."<< std::endl;
return false;
}
// Evaluate the surface traction
Vec3 T = (*tracFld)(X,normal);
// Store the traction value for vizualization
if (!T.isZero()) tracVal[X] = T;
// Integrate the force vector
Vector& ES = *eS;
for (size_t a = 1; a <= N.size(); a++)
for (short int d = 1; d <= nsd; d++)
ES(nf*(a-1)+d) += T[d-1]*N(a)*detJW;
return getIntegralResult(elmInt);
}
bool Stokes::writeGlvT (VTF* vtf, int iStep, int& nBlock) const
{
if (tracVal.empty())
return true;
else if (!vtf)
return false;
// Write boundary tractions as discrete point vectors to the VTF-file
if (!vtf->writeVectors(tracVal,++nBlock))
return false;
return vtf->writeVblk(nBlock,"Tractions",1,iStep);
}
bool Stokes::evalSol (Vector& s, const Vector& N,
const Matrix& dNdX, const Vec3& X,
const std::vector<int>& MNPC) const
{
return false;
}
bool Stokes::evalSol (Vector& s, const TensorFunc& asol, const Vec3& X) const
{
s = asol(X);
return true;
}
const char* Stokes::getFieldLabel(size_t i, const char* prefix) const
{
static const char* s[3] = { "u_x","u_y","u_z" };
static std::string label;
if (prefix)
label = prefix + std::string(" ");
else
label.clear();
label += s[i];
return label.c_str();
}
NormBase* Stokes::getNormIntegrand (AnaSol* asol) const
{
return new StokesNorm(*const_cast<Stokes*>(this),asol);
}
NormBase* Stokes::getForceIntegrand (AnaSol* asol) const
{
return new StokesForce(*const_cast<Stokes*>(this),asol);
}
bool Stokes::getIntegralResult (LocalIntegral*& elmInt) const
{
elmInt = myMats;
return elmInt ? true : false;
}
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 == SIM::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
{
int i, k, l;
double value, eps, epsh;
const int nsd = dNdX.cols();
const int nf = nsd+1;
const int nbf = N.size();
if (!anasol->hasScalarSol() && !anasol->hasVectorSol()) {
std::cerr <<" *** StokesNorm::evalInt: No analytical solution."<< std::endl;
return false;
}
ElmNorm* eNorm = dynamic_cast<ElmNorm*>(elmInt);
if (!eNorm) {
std::cerr <<" *** StokesNorm::evalInt: No norm array."<< std::endl;
return false;
}
// Find current element solution vector
Vector* eV = const_cast<Vector*>(problem.getElementSolution());
if (!eV) {
std::cerr <<" *** StokesNorm::evalInt: No solution"<< std::endl;
return false;
}
Vector& EV = *eV;
ElmNorm& pnorm = *eNorm;
int ip = 0;
// Pressure norms
if (anasol->getScalarSol()) {
// Analytical pressure
real P = (*anasol->getScalarSol())(X);
// Computed pressure
real Ph = 0.0;
for (i = 1;i <= nbf;i++)
Ph += EV(i*nf)*N(i);
// L2-norm of pressure error
P -= Ph;
pnorm[ip++] += P*P*detJW;
}
// Velocity norms
if (anasol->getVectorSol()) {
// Analytical velocity
Vec3 U = (*anasol->getVectorSol())(X);
// Computed velocity
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);
// L2-norm of velocity error
U -= Uh;
pnorm[ip++] += U*U*detJW;
}
// Pressure gradient norm
if (anasol->getScalarSecSol()) {
// Analytical pressure gradient
Vec3 dP = (*anasol->getScalarSecSol())(X);
// Computed pressure gradient
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);
// H1-seminorm of pressure
dP -= dPh;
pnorm[ip++] += dP*dP*detJW;
}
// Energy norms
if (anasol->getVectorSecSol()) {
// Viscosity
const double mu = problem.getViscosity(X);
// Analytical velocity gradient
Tensor gradU = (*anasol->getVectorSecSol())(X);
// Numerical velocity gradient
Tensor gradUh(nsd);
problem.strain(dNdX,gradUh);
if (problem.getFormulation() == SIM::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 {
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() == SIM::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;
}
size_t StokesForce::getNoFields() const
{
size_t nsd = problem.getNoSpaceDim();
if (anasol)
return 2*nsd + 1;
else
return nsd;
}

345
src/Integrands/Stokes.h
View File

@@ -1,345 +0,0 @@
// $Id: Stokes.h,v 1.4 2011-02-08 12:21:31 rho Exp $
//==============================================================================
//!
//! \file Stokes.h
//!
//! \date Sep 29 2010
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Base class for Integrand implementation of Stokes problems.
//!
//==============================================================================
#ifndef _STOKES_H
#define _STOKES_H
#include "IntegrandBase.h"
#include "ElmMats.h"
#include "Vec3.h"
#include "AnaSol.h"
#include <map>
class VTF;
/*!
\brief Base class representing the integrand of Stokes problems.
\details Implements common features for Stokes/Navier-Stokes problems.
*/
class Stokes : public Integrand
{
public:
//! \brief The default constructor initializes all pointers to zero.
//! \param[in] n Number of spatial dimensions
//! \param[in] form The solution formulation to use
//! \param[in] itg_type The integrand type to use
Stokes(unsigned short int n,
SIM::Formulation form = SIM::LAPLACE,
int itg_type = 1);
//! \brief The destructor frees dynamically allocated objects.
virtual ~Stokes();
//! \brief Defines the traction field to use in Neumann boundary conditions.
void setTraction(TractionFunc* tf) { tracFld = tf; }
//! \brief Clears the integration point traction values.
void clearTracVal() { tracVal.clear(); }
//! \brief Defines the gravitation vector.
void setGravity(double gx, double gy, double gz)
{ g[0] = gx; g[1] = gy; g[2] = gz; }
//! \brief Defines fluid properties for current volume patch.
//! \param[in] Density Mass density
//! \param[in] Viscosity Dynamic viscosity
void setFluidProperties(double Density, double Viscosity)
{ rho = Density; mu = Viscosity; }
//! \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
//! \param[in] detJW Jacobian determinant times integration point weight
//! \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
virtual bool evalBou(LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X, const Vec3& normal) const;
//! \brief Evaluates the secondary solution at current integration point.
//! \param[out] s The solution field values
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] MNPC Matrix of nodal point correspondance
virtual bool evalSol(Vector& s,
const Vector& N, const Matrix& dNdX,
const Vec3& X, const std::vector<int>& MNPC) const;
//! \brief Evaluates the analytical secondary solution at the given point.
//! \param[out] s The solution field values
//! \param[in] asol The analytical solution field
//! \param[in] X Cartesian coordinates of current integration point
virtual bool evalSol(Vector& s, const TensorFunc& asol, const Vec3& X) const;
//! \brief Writes the surface tractions for a given time step to VTF-file.
//! \param vtf The VTF-file object to receive the tractions
//! \param[in] iStep Load/time step identifier
//! \param nBlock Running result block counter
bool writeGlvT(VTF* vtf, int iStep, int& nBlock) const;
//! \brief Returns a pointer to an Integrand for solution norm 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* 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; }
//! \brief Returns which integrand to use.
virtual int getIntegrandType() const { return integrandType; }
//! \brief Returns the mass density at current point.
virtual double getMassDensity(const Vec3&) const { return rho; }
//! \brief Returns the body force per unit mass
//! \param[in] X Cartesian coordinate of current integration point
//! \param[out] f Body force vector per unit mass
virtual bool getBodyForce(const Vec3& X, Vector& f) const;
//! \brief Returns viscosity at current point.
virtual double getViscosity(const Vec3&) const { return mu; }
//! \brief Returns a pointer to current element solution vector.
const Vector* getElementSolution() const { return eVs[0]; }
//! \brief Returns a pointer to current element solution vector.
virtual const Vector* getElementVelocity() const { return 0; }
//! \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; }
//! \brief Returns the number of secondary solution fields.
virtual size_t getNoSecFields() const { return nsd*nsd; }
//! \brief Returns the name of the primary solution field component.
//! \param[in] i Field component index
//! \param[in] prefix Name prefix for all components
virtual const char* getFieldLabel(size_t i, const char* prefix = 0) const;
//! \brief Returns the name of a secondary solution field component.
//! \param[in] i Field component index
//! \param[in] prefix Name prefix for all components
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
bool getIntegralResult(LocalIntegral*& elmInt) const;
// Problem parameters
unsigned short int nsd; //!< Number of space dimensions (1, 2 or, 3)
unsigned short int nf; //!< Number of primary field variables
SIM::Formulation formulation; //!< Problem formulation flag
int integrandType; //!< Integrand type flag
// Physical properties (constant)
double mu; //!< Dynamic viscosity
double rho; //!< Mass density
double g[3]; //!< Gravitation vector
// Finite element quantities
Matrix* eM; //!< Pointer to element coefficient matrix
Vector* eS; //!< Pointer to element right-hand-side vector
std::vector<Vector*> eVs; //!< Pointers to element solution vectors
ElmMats* myMats; //!< Local element matrices
TractionFunc* tracFld; //!< Pointer to boundary traction field
mutable std::map<Vec3,Vec3> tracVal; //!< Traction field point values
};
/*!
\brief Class representing the integrand of Stokes energy norms.
*/
class StokesNorm : 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)
StokesNorm(Stokes& p, AnaSol* a) : problem(p), anasol(a) {}
//! \brief Empty destructor.
virtual ~StokesNorm() {}
//! \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
//! \param[in] detJW Jacobian determinant times integration point weight
//! \param[in] N Basis function values
//! \param[in] dNdX Basis function gradients
//! \param[in] X Cartesian coordinates of current integration point
virtual bool evalInt(LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const;
//! \brief Returns the number of primary norm quantities.
virtual size_t getNoFields() const { return anasol ? 2 : 0; }
//! \brief Returns the number of secondary norm quantities.
virtual size_t getNoSecFields() const { return anasol ? 3 : 1; }
protected:
Stokes& problem; //!< The problem-specific data
AnaSol* anasol; //!< Analytical solution fields
};
/*!
\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