diff --git a/src/Integrands/AnalyticSolutionsStokes.C b/src/Integrands/AnalyticSolutionsStokes.C deleted file mode 100644 index 07fac333..00000000 --- a/src/Integrands/AnalyticSolutionsStokes.C +++ /dev/null @@ -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(&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(&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(&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(&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; -} diff --git a/src/Integrands/AnalyticSolutionsStokes.h b/src/Integrands/AnalyticSolutionsStokes.h deleted file mode 100644 index 9f9413a2..00000000 --- a/src/Integrands/AnalyticSolutionsStokes.h +++ /dev/null @@ -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 diff --git a/src/Integrands/ChorinPressCorr.C b/src/Integrands/ChorinPressCorr.C deleted file mode 100644 index 1b0d6cea..00000000 --- a/src/Integrands/ChorinPressCorr.C +++ /dev/null @@ -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& 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& MNPC1, - const std::vector& 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 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& 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& MNPC1, - const std::vector& 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& 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; -} diff --git a/src/Integrands/ChorinPressCorr.h b/src/Integrands/ChorinPressCorr.h deleted file mode 100644 index 5909ca10..00000000 --- a/src/Integrands/ChorinPressCorr.h +++ /dev/null @@ -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 - -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& 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& MNPC1, - const std::vector& 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& 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& MNPC1, - const std::vector& 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& 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& 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 eVs; //!< Pointer to element velocity vectors - std::vector ePs; //!< Pointer to element pressure vectors - - std::vector usol; //!< Velocity solution vectors for patch - - ElmMats* myMats; //!< Local element matrices - - VecFunc* tracFld; //!< Pointer to boundary traction field - mutable std::map tracVal; //!< Traction field point values -}; - -#endif diff --git a/src/Integrands/ChorinVelCorr.C b/src/Integrands/ChorinVelCorr.C deleted file mode 100644 index 6558d0bb..00000000 --- a/src/Integrands/ChorinVelCorr.C +++ /dev/null @@ -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& 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& MNPC1, - const std::vector& 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& 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& MNPC1, - const std::vector& 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& 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; -} diff --git a/src/Integrands/ChorinVelCorr.h b/src/Integrands/ChorinVelCorr.h deleted file mode 100644 index 6fd896dd..00000000 --- a/src/Integrands/ChorinVelCorr.h +++ /dev/null @@ -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 - -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& 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& MNPC1, - const std::vector& 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& 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& MNPC1, - const std::vector& 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& 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& MNPC1, - const std::vector& 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 eVs; //!< Pointer to element velocity vectors - std::vector ePs; //!< Pointer to element pressure vectors - - std::vector psol; //!< Pressure solution vectors for patch - - ElmMats* myMats; //!< Local element matrices - - TractionFunc* tracFld; //!< Pointer to boundary traction field - mutable std::map tracVal; //!< Traction field point values -}; - -#endif diff --git a/src/Integrands/ChorinVelPred.C b/src/Integrands/ChorinVelPred.C deleted file mode 100644 index cc1933c0..00000000 --- a/src/Integrands/ChorinVelPred.C +++ /dev/null @@ -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& 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& MNPC1, - const std::vector& 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& 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& MNPC1, - const std::vector& 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(this),asol); -} - - -NormBase* ChorinVelPred::getForceIntegrand (AnaSol* asol) const -{ - return new ChorinStokesForce(*const_cast(this),asol); -} - - -bool ChorinVelPred::evalSol (Vector& s, const Vector& N, - const Matrix& dNdX, const Vec3& X, - const std::vector& 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& MNPC1, - const std::vector& 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(elmInt); - if (!eNorm) { - std::cerr <<" *** ChorinStokesNorm::evalInt: No norm array."<< std::endl; - return false; - } - - // Find current element velocity vector - Vector* eV = const_cast(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(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(elmInt); - if (!eNorm) { - std::cerr <<" *** ChorinStokesNorm::evalInt: No norm array."<< std::endl; - return false; - } - - // Find current element velocity vector - Vector* eV = const_cast(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(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& MNPC) -{ - return problem.initElementBou(MNPC); -} - - -bool ChorinStokesForce::initElementBou(const std::vector& MNPC1, - const std::vector& 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(elmInt); - if (!eNorm) { - std::cerr <<" *** StokesForce::evalBou: No force array."<< std::endl; - return false; - } - - ChorinVelPred* prob = static_cast(&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; -} - - diff --git a/src/Integrands/ChorinVelPred.h b/src/Integrands/ChorinVelPred.h deleted file mode 100644 index 824e115a..00000000 --- a/src/Integrands/ChorinVelPred.h +++ /dev/null @@ -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 - -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& 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& MNPC1, - const std::vector& 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& 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& MNPC1, - const std::vector& 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& 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& MNPC1, - const std::vector& 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 ePs; //!< Pointers to element pressure vectors - std::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& 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& MNPC1, - const std::vector& 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 diff --git a/src/Integrands/ChorinVelPredBDF2.C b/src/Integrands/ChorinVelPredBDF2.C deleted file mode 100644 index 1651b9bf..00000000 --- a/src/Integrands/ChorinVelPredBDF2.C +++ /dev/null @@ -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); -} diff --git a/src/Integrands/ChorinVelPredBDF2.h b/src/Integrands/ChorinVelPredBDF2.h deleted file mode 100644 index 17e8eea1..00000000 --- a/src/Integrands/ChorinVelPredBDF2.h +++ /dev/null @@ -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 diff --git a/src/Integrands/ChorinVelPredBE.C b/src/Integrands/ChorinVelPredBE.C deleted file mode 100644 index 61eb6ee2..00000000 --- a/src/Integrands/ChorinVelPredBE.C +++ /dev/null @@ -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); -} - diff --git a/src/Integrands/ChorinVelPredBE.h b/src/Integrands/ChorinVelPredBE.h deleted file mode 100644 index 65ee0cda..00000000 --- a/src/Integrands/ChorinVelPredBE.h +++ /dev/null @@ -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 diff --git a/src/Integrands/NavierStokesG2.C b/src/Integrands/NavierStokesG2.C deleted file mode 100644 index f511478e..00000000 --- a/src/Integrands/NavierStokesG2.C +++ /dev/null @@ -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& 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; -} diff --git a/src/Integrands/NavierStokesG2.h b/src/Integrands/NavierStokesG2.h deleted file mode 100644 index 900137f1..00000000 --- a/src/Integrands/NavierStokesG2.h +++ /dev/null @@ -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& MNPC) const; - - //! \brief Returns which integrand to use - int getIntegrandType() const { return 2; } -}; - -#endif diff --git a/src/Integrands/NavierStokesG2CN.C b/src/Integrands/NavierStokesG2CN.C deleted file mode 100644 index d63c856d..00000000 --- a/src/Integrands/NavierStokesG2CN.C +++ /dev/null @@ -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); -} diff --git a/src/Integrands/NavierStokesG2CN.h b/src/Integrands/NavierStokesG2CN.h deleted file mode 100644 index c5d79f2d..00000000 --- a/src/Integrands/NavierStokesG2CN.h +++ /dev/null @@ -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 diff --git a/src/Integrands/NavierStokesG2GenTheta.C b/src/Integrands/NavierStokesG2GenTheta.C deleted file mode 100644 index fc393636..00000000 --- a/src/Integrands/NavierStokesG2GenTheta.C +++ /dev/null @@ -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); -} diff --git a/src/Integrands/NavierStokesG2GenTheta.h b/src/Integrands/NavierStokesG2GenTheta.h deleted file mode 100644 index 3b8107e3..00000000 --- a/src/Integrands/NavierStokesG2GenTheta.h +++ /dev/null @@ -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 diff --git a/src/Integrands/NavierStokesG2MP.C b/src/Integrands/NavierStokesG2MP.C deleted file mode 100644 index e8d2bff4..00000000 --- a/src/Integrands/NavierStokesG2MP.C +++ /dev/null @@ -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); -} diff --git a/src/Integrands/NavierStokesG2MP.h b/src/Integrands/NavierStokesG2MP.h deleted file mode 100644 index 7e8d528f..00000000 --- a/src/Integrands/NavierStokesG2MP.h +++ /dev/null @@ -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 diff --git a/src/Integrands/StabilizedNavierStokes.C b/src/Integrands/StabilizedNavierStokes.C deleted file mode 100644 index 6e3bca3e..00000000 --- a/src/Integrands/StabilizedNavierStokes.C +++ /dev/null @@ -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); -} diff --git a/src/Integrands/StabilizedNavierStokes.h b/src/Integrands/StabilizedNavierStokes.h deleted file mode 100644 index f2000a32..00000000 --- a/src/Integrands/StabilizedNavierStokes.h +++ /dev/null @@ -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 diff --git a/src/Integrands/StabilizedStokes.C b/src/Integrands/StabilizedStokes.C deleted file mode 100644 index 8a4c0a1e..00000000 --- a/src/Integrands/StabilizedStokes.C +++ /dev/null @@ -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); -} - - diff --git a/src/Integrands/StabilizedStokes.h b/src/Integrands/StabilizedStokes.h deleted file mode 100644 index 404f45f6..00000000 --- a/src/Integrands/StabilizedStokes.h +++ /dev/null @@ -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 diff --git a/src/Integrands/Stokes.C b/src/Integrands/Stokes.C deleted file mode 100644 index e6f23859..00000000 --- a/src/Integrands/Stokes.C +++ /dev/null @@ -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& 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& 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& 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(this),asol); -} - - -NormBase* Stokes::getForceIntegrand (AnaSol* asol) const -{ - return new StokesForce(*const_cast(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(&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& MNPC) -{ - return problem.initElement(MNPC); -} - - -bool StokesNorm::initElementBou (const std::vector& 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(elmInt); - if (!eNorm) { - std::cerr <<" *** StokesNorm::evalInt: No norm array."<< std::endl; - return false; - } - - // Find current element solution vector - Vector* eV = const_cast(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& 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(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; -} - - diff --git a/src/Integrands/Stokes.h b/src/Integrands/Stokes.h deleted file mode 100644 index 208605fa..00000000 --- a/src/Integrands/Stokes.h +++ /dev/null @@ -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 - -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& 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& MNPC1, - const std::vector& 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& 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& MNPC1, - const std::vector& 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& 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 eVs; //!< Pointers to element solution vectors - - ElmMats* myMats; //!< Local element matrices - - TractionFunc* tracFld; //!< Pointer to boundary traction field - - mutable std::map 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& 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& 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& 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& MNPC1, - const std::vector& 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