From c15bbbd25183f443bd12d6ea11c7d14dc72bf1c0 Mon Sep 17 00:00:00 2001 From: rho Date: Wed, 23 Feb 2011 08:45:35 +0000 Subject: [PATCH] Added class for calculating boundary forces. Changed stress and strain computations. git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@817 e10b68d5-8a6e-419e-a041-bce267b0401d --- src/Integrands/ChorinVelPred.C | 232 +++++++++++++++++--------- src/Integrands/ChorinVelPred.h | 79 ++++++++- src/Integrands/Stokes.C | 288 ++++++++++++++++++++++++++------- src/Integrands/Stokes.h | 128 +++++++++++++++ 4 files changed, 580 insertions(+), 147 deletions(-) diff --git a/src/Integrands/ChorinVelPred.C b/src/Integrands/ChorinVelPred.C index 569198d6..62e58fb1 100644 --- a/src/Integrands/ChorinVelPred.C +++ b/src/Integrands/ChorinVelPred.C @@ -97,8 +97,19 @@ bool ChorinVelPred::initElementBou(const std::vector& MNPC) if (eS) eS->resize(nsd*nen,true); + int ierr = 0; + if (!eVs.empty() && !primsol.empty()) + if (ierr = utl::gather(MNPC,nsd,primsol[0],*eVs[0])) + std::cerr <<" *** ChorinVelPred::initElementBou: Detected " + << ierr <<" node numbers out of range."<< std::endl; + + if (!ePs.empty() && !psol.empty()) + if (ierr = utl::gather(MNPC,1,psol[0],*ePs[0])) + std::cerr <<" *** ChorinVelPred::initElementBou: Detected " + << ierr <<" node numbers out of range."<< std::endl; + myMats->withLHS = false; - return true; + return ierr == 0; } @@ -111,10 +122,20 @@ bool ChorinVelPred::initElementBou(const std::vector& MNPC1, if (eS) eS->resize(nsd*nen1,true); - myMats->withLHS = false; - return true; -} + int ierr = 0; + if (!eVs.empty() && !primsol.empty()) + if (ierr = utl::gather(MNPC1,nsd,primsol[0],*eVs[0])) + std::cerr <<" *** ChorinVelPred::initElementBou: Detected " + << ierr <<" node numbers out of range."<< std::endl; + + if (!ePs.empty() && !psol.empty()) + if (ierr = utl::gather(MNPC2,1,psol[0],*ePs[0])) + std::cerr <<" *** ChorinVelPred::initElementBou: Detected " + << ierr <<" node numbers out of range."<< std::endl; + myMats->withLHS = false; + return ierr == 0; +} NormBase* ChorinVelPred::getNormIntegrand (AnaSol* asol) const @@ -123,6 +144,12 @@ NormBase* ChorinVelPred::getNormIntegrand (AnaSol* asol) const } +NormBase* ChorinVelPred::getForceIntegrand (AnaSol* asol) const +{ + return new ChorinStokesForce(*const_cast(this),asol); +} + + bool ChorinVelPred::evalSol (Vector& s, const Vector& N, const Matrix& dNdX, const Vec3& X, const std::vector& MNPC) const @@ -193,7 +220,7 @@ bool ChorinVelPred::evalSol (Vector& s, const TensorFunc& asol, } -bool ChorinVelPred::strain(const Matrix& dNdX, SymmTensor& eps) const +bool ChorinVelPred::strain(const Matrix& dNdX, Tensor& eps) const { int i, k, l; @@ -203,27 +230,34 @@ bool ChorinVelPred::strain(const Matrix& dNdX, SymmTensor& eps) const return false; } + eps.zero(); if (!eVs.empty() && !eVs[0]->empty() && eps.dim() > 0) { Vector& EV = *(eVs[0]); for (i = 1;i <= dNdX.rows();i++) for (k = 1;k <= nsd;k++) - for (l = 1;l <= k;l++) - eps(k,l) = dNdX(i,l)*EV((i-1)*nsd+k); + dNdX(i,k)*EV((i-1)*nsd+l); + for (l = 1;l <= nsd;l++) + eps(k,l) += dNdX(i,l)*EV((i-1)*nsd+k); + } + else + return false; - eps *= mu; + if (formulation == Stokes::STRESS) { + eps += eps.transpose(); + eps *= 0.5; } return true; } -bool ChorinVelPred::stress(const Vector& N, const Matrix& dNdX, SymmTensor& sigma) const +bool ChorinVelPred::stress(const Vector& N, const Matrix& dNdX, Tensor& sigma) const { int i, k; // Compute strain if (!this->strain(dNdX,sigma)) return false; + sigma *= mu; // Compute pressure if (!ePs.empty() && !ePs[0]->empty() && sigma.dim() > 0) { @@ -231,7 +265,7 @@ bool ChorinVelPred::stress(const Vector& N, const Matrix& dNdX, SymmTensor& sigm Vector& EP = *(ePs[0]); for (i = 1;i <= N.size();i++) P += EP(i)*N(i); - + // Add pressure and strain contributions for (k = 1;k <= nsd;k++) sigma(k,k) -= P; @@ -243,7 +277,7 @@ bool ChorinVelPred::stress(const Vector& N, const Matrix& dNdX, SymmTensor& sigm bool ChorinVelPred::stress(const Vector& N1, const Vector& N2, const Matrix& dN1dX, const Matrix& dN2dX, - SymmTensor& sigma) const + Tensor& sigma) const { int i, k; @@ -337,7 +371,7 @@ bool ChorinStokesNorm::evalInt (LocalIntegral*& elmInt, double detJW, Vec3 Uh; Uh = 0.0; for (i = 1;i <= N.size();i++) for (k = 1;k <= nsd;k++) - Uh[k-1] += EV((i-1)*nsd + k)*N(i); + Uh[k-1] += EV((i-1)*nsd + k)*N(i); // L2-norm of velocity error U -= Uh; @@ -353,7 +387,7 @@ bool ChorinStokesNorm::evalInt (LocalIntegral*& elmInt, double detJW, Vec3 dPh; dPh = 0.0; for (i = 1;i <= N.size();i++) for (k = 1;k <= nsd;k++) - dPh[k-1] += EP(i)*dNdX(i,k); + dPh[k-1] += EP(i)*dNdX(i,k); // H1-seminorm of pressure dP -= dPh; @@ -368,44 +402,27 @@ bool ChorinStokesNorm::evalInt (LocalIntegral*& elmInt, double detJW, // Analytical velocity gradient Tensor gradU = (*anasol->getVectorSecSol())(X); - // Computed velocity gradient - Matrix gradUh(nsd,nsd); - gradUh.fill(0.0); - for (i = 1;i <= N.size();i++) - for (k = 1;k <= nsd;k++) - for (l = 1;l <= nsd;l++) - gradUh(k,l) += EV((i-1)*nsd + k)*dNdX(i,l); + // Numerical velocity gradient + Tensor gradUh(3); + gradUh.zero(); + problem.strain(dNdX,gradUh); - if (problem.getFormulation() == Stokes::LAPLACE) { - for (k = 1;k <= nsd;k++) - for (l = 1;l <= nsd;l++) { - // Energy norm of analytical solution - pnorm[ip] += mu*gradU(k,l)*gradU(k,l)*detJW; - // Energy norm of numerical solution - pnorm[ip+1] += mu*gradUh(k,l)*gradUh(k,l)*detJW; - // Energy norm of error - value = gradU(k,l)-gradUh(k,l); - pnorm[ip+2] += mu*value*value*detJW; - } + if (problem.getFormulation() == Stokes::STRESS) { + pnorm[ip++] += 2.0*mu*gradU.innerProd(gradU)*detJW; + pnorm[ip++] += 2.0*mu*gradUh.innerProd(gradUh)*detJW; + gradUh *= -1.0; + gradU += gradUh; + pnorm[ip++] += 2.0*mu*gradU.innerProd(gradU)*detJW; } else { - for (k = 1;k <= nsd;k++) - for (l = 1;l <= nsd;l++) { - // Strain of analytical solution - eps = 0.5*(gradU(k,l) + gradU(l,k)); - // Strain of computed solution - epsh = 0.5*(gradUh(k,l) + gradUh(l,k)); - // Energy norm of analytical solution - pnorm[ip] += mu*eps*eps*detJW; - // Energy norm of computed solution - pnorm[ip+1] += mu*epsh*epsh*detJW; - // Energy norm of error - eps -= epsh; - pnorm[ip+2] += mu*epsh*epsh*detJW; - } + pnorm[ip++] += mu*gradU.innerProd(gradU)*detJW; + pnorm[ip++] += mu*gradUh.innerProd(gradUh)*detJW; + gradUh *= -1.0; + gradU += gradUh; + pnorm[ip++] += mu*gradU.innerProd(gradU)*detJW; } } - + return true; } @@ -483,7 +500,7 @@ bool ChorinStokesNorm::evalInt(LocalIntegral*& elmInt, Vec3 Uh; Uh = 0.0; for (i = 1;i <= nbfU;i++) for (k = 1;k <= nsd;k++) - Uh[k-1] += EV((i-1)*nsd + k)*N1(i); + Uh[k-1] += EV((i-1)*nsd + k)*N1(i); // L2-norm of velocity error U -= Uh; @@ -499,7 +516,7 @@ bool ChorinStokesNorm::evalInt(LocalIntegral*& elmInt, Vec3 dPh; dPh = 0.0; for (i = 1;i <= nbfP;i++) for (k = 1;k <= nsd;k++) - dPh[k-1] += EP(i)*dN2dX(i,k); + dPh[k-1] += EP(i)*dN2dX(i,k); // H1-seminorm of pressure dP -= dPh; @@ -515,40 +532,22 @@ bool ChorinStokesNorm::evalInt(LocalIntegral*& elmInt, Tensor gradU = (*anasol->getVectorSecSol())(X); // Computed velocity gradient - Matrix gradUh(nsd,nsd); - gradUh.fill(0.0); - for (i = 1;i <= nbfU;i++) - for (k = 1;k <= nsd;k++) - for (l = 1;l <= nsd;l++) - gradUh(k,l) += EV((i-1)*nsd + k)*dN1dX(i,l); - - if (problem.getFormulation() == Stokes::LAPLACE) { - for (k = 1;k <= nsd;k++) - for (l = 1;l <= nsd;l++) { - // Energy norm of analytical solution - pnorm[ip] += mu*gradU(k,l)*gradU(k,l)*detJW; - // Energy norm of numerical solution - pnorm[ip+1] += mu*gradUh(k,l)*gradUh(k,l)*detJW; - // Energy norm of error - value = gradU(k,l)-gradUh(k,l); - pnorm[ip+2] += mu*value*value*detJW; - } + Tensor gradUh(nsd); + problem.strain(dN1dX,gradUh); + + if (problem.getFormulation() == Stokes::STRESS) { + pnorm[ip++] += 2.0*mu*gradU.innerProd(gradU)*detJW; + pnorm[ip++] += 2.0*mu*gradUh.innerProd(gradUh)*detJW; + gradUh *= -1.0; + gradU += gradUh; + pnorm[ip++] += 2.0*mu*gradU.innerProd(gradU)*detJW; } else { - for (k = 1;k <= nsd;k++) - for (l = 1;l <= nsd;l++) { - // Strain of analytical solution - eps = 0.5*(gradU(k,l) + gradU(l,k)); - // Strain of computed solution - epsh = 0.5*(gradUh(k,l) + gradUh(l,k)); - // Energy norm of analytical solution - pnorm[ip] += mu*eps*eps*detJW; - // Energy norm of computed solution - pnorm[ip+1] += mu*epsh*epsh*detJW; - // Energy norm of error - eps -= epsh; - pnorm[ip+2] += mu*epsh*epsh*detJW; - } + pnorm[ip++] += mu*gradU.innerProd(gradU)*detJW; + pnorm[ip++] += mu*gradUh.innerProd(gradUh)*detJW; + gradUh *= -1.0; + gradU += gradUh; + pnorm[ip++] += mu*gradU.innerProd(gradU)*detJW; } } @@ -556,3 +555,76 @@ bool ChorinStokesNorm::evalInt(LocalIntegral*& elmInt, } +bool ChorinStokesForce::initElementBou(const std::vector& 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() == Stokes::STRESS) + sigma += sigma.transpose(); + + // Analytical stress + sigma *= mu; + for (i = 1;i <= nsd;i++) + sigma(i,i) -= P; + + // Analytical traction + Vec3 t = sigma*normal; + + for (i = 0;i < nsd;i++) + pnorm[ip++] += t[i]*detJW; + + // Error in traction + t -= th; + pnorm[ip++] += t*t*detJW; + } + + return true; +} + + diff --git a/src/Integrands/ChorinVelPred.h b/src/Integrands/ChorinVelPred.h index 96670a72..135fc464 100644 --- a/src/Integrands/ChorinVelPred.h +++ b/src/Integrands/ChorinVelPred.h @@ -168,6 +168,12 @@ public: //! \param[in] asol Pointer to analytical solution field (optional) virtual NormBase* getNormIntegrand(AnaSol* asol = 0) const; + //! \brief Returns a pointer to an Integrand for boundary force evaluation. + //! \note The Integrand is allocated dynamically and has to be deleted + //! manually when leaving the scope of the pointer returned. + //! \param[in] asol Pointer to analytical solution field (optional) + virtual NormBase* getForceIntegrand(AnaSol* asol = 0) const; + //! \brief Accesses the velocity solution vectors of current patch Vector& getVelocity(int n = 0) { return this->getSolution(n); } @@ -192,17 +198,16 @@ public: //! \brief Get number of pressure solutions size_t getNoPressures() const { return psol.size(); } - protected: //! \brief Calculates viscous part of stress tensor at current point. //! \param[in] dNdX Basis function gradients at current point //! \param[out] eps Strain tensor at current point - virtual bool strain(const Matrix& dNdX, SymmTensor& eps) const; + virtual bool strain(const Matrix& dNdX, Tensor& eps) const; //! \brief Calculates the (Cauchy) stress tensor at current point //! \param[in] N Basis functions at current point //! \param[in] dNdX Basis function gradients at current point //! \param[out] sigma Strain tensor at current point - bool stress(const Vector& N, const Matrix& dNdX, SymmTensor& sigma) const; + bool stress(const Vector& N, const Matrix& dNdX, Tensor& sigma) const; //! \brief Calculates the (Cauchy) stress tensor at current point //! \param[in] N1 Velocity basis functions at current point @@ -210,8 +215,9 @@ public: //! \param[out] sigma Strain tensor at current point bool stress(const Vector& N1, const Vector& N2, const Matrix& dN1dX, const Matrix& dN2dX, - SymmTensor& sigma) const; - + Tensor& sigma) const; + + protected: bool incPressure; //!< Incremental pressure formulation bool mixedFEM; //!< Mixed FE formulation @@ -273,4 +279,67 @@ class ChorinStokesNorm : public StokesNorm }; +/*! + \brief Class representing for computing boundary force +*/ + +class ChorinStokesForce : public StokesForce +{ +public: + //! \brief The only constructor initializes its data members. + //! \param[in] p The Stokes problem to evaluate norms for + //! \param[in] a The analytical velocity and pressure fields (optional) + ChorinStokesForce(ChorinVelPred& p, AnaSol* a = 0) : StokesForce(p,a) {} + //! \brief Empty destructor. + virtual ~ChorinStokesForce() {} + + //! \brief Initializes current element for boundary integration. + //! \param[in] MNPC Matrix of nodal point correspondance for current element + virtual bool initElementBou(const std::vector& 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/Stokes.C b/src/Integrands/Stokes.C index da6b8fdb..1b4d7f06 100644 --- a/src/Integrands/Stokes.C +++ b/src/Integrands/Stokes.C @@ -70,8 +70,14 @@ 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 true; + return ierr == 0; } @@ -152,6 +158,12 @@ NormBase* Stokes::getNormIntegrand (AnaSol* asol) const } +NormBase* Stokes::getForceIntegrand (AnaSol* asol) const +{ + return new StokesForce(*const_cast(this),asol); +} + + bool Stokes::getIntegralResult (LocalIntegral*& elmInt) const { elmInt = myMats; @@ -159,15 +171,138 @@ bool Stokes::getIntegralResult (LocalIntegral*& elmInt) const } +bool Stokes::getBodyForce(const Vec3& X, Vector& f) const +{ + const Vec4* Y = dynamic_cast(&X); + if (!Y) return false; + + const double PI = 3.141592653589793238462; + const double x = X[0]; + const double y = X[1]; + const double t = Y->t; + + f.fill(0.0); + + + f(1) = rho*pow(sin(PI*x),2.0)*sin(2.0*PI*y)*cos(t) + + rho*PI*pow(sin(PI*x)*sin(2.0*PI*y)*sin(t),2.0)*sin(2.0*PI*x) - + 2.0*rho*PI*pow(sin(PI*x)*sin(PI*y)*sin(t),2.0)*sin(2.0*PI*x)*cos(2.0*PI*y) - + PI*sin(PI*x)*sin(PI*y)*sin(t) - + 2.0*pow(PI,2.0)*mu*sin(2.0*PI*y)*sin(t)*(cos(2*PI*x)-2*pow(sin(PI*x),2.0)); + + f(2) = -rho*sin(2.0*PI*x)*pow(sin(PI*y),2.0)*cos(t) - + 2.0*rho*PI*pow(sin(PI*x)*sin(PI*y)*sin(t),2.0)*cos(2.0*PI*x)*sin(2.0*PI*y) + + rho*PI*pow(sin(2.0*PI*x)*sin(PI*y)*sin(t),2.0)*sin(2.0*PI*y) + + PI*cos(PI*x)*cos(PI*y)*sin(t) - + 4*mu*pow(PI*sin(PI*y),2.0)*sin(2.0*PI*x)*sin(t) + + 2.0*mu*pow(PI,2.0)*sin(2.0*PI*x)*cos(2.0*PI*y)*sin(t); + + return true; +} + + +bool Stokes::strain(const Matrix& dNdX, Tensor& eps) const +{ + int i, k, l; + + const int nf = nsd+1; + + if (dNdX.cols() < nsd) { + std::cerr <<" *** Stokes::strain: Invalid dimension on dNdX " + << dNdX.rows() <<" "<< dNdX.cols() << std::endl; + return false; + } + + eps.zero(); + if (!eVs.empty() && !eVs[0]->empty() && eps.dim() > 0) { + Vector& EV = *(eVs[0]); + for (i = 1;i <= dNdX.rows();i++) + for (k = 1;k <= nsd;k++) + for (l = 1;l <= nsd;l++) + eps(k,l) += dNdX(i,l)*EV((i-1)*nf+k); + } + else + return false; + + if (formulation == Stokes::STRESS) { + eps += eps.transpose(); + eps *= 0.5; + } + + return true; +} + + +bool Stokes::stress(const Vector& N, const Matrix& dNdX, Tensor& sigma) const +{ + int i, k; + + const int nf = nsd+1; + + // Compute strain + if (!this->strain(dNdX,sigma)) + return false; + sigma *= mu; + + // Compute pressure + if (!eVs.empty() && !eVs[0]->empty() && sigma.dim() > 0) { + double P = 0.0; + Vector& EV = *(eVs[0]); + for (i = 1;i <= N.size();i++) + P += EV(i*nf)*N(i); + + // Add pressure and strain contributions + for (k = 1;k <= nsd;k++) + sigma(k,k) -= P; + } + + return true; +} + + +bool Stokes::stress(const Vector& N1, const Vector& N2, + const Matrix& dN1dX, const Matrix& dN2dX, + Tensor& sigma) const +{ + int i, k; + + const int nf = nsd+1; + + // Compute strain + if (!this->strain(dN1dX,sigma)) + return false; + + // Add pressure + if (!eVs.empty() && eVs[0]) { + double P = 0.0; + Vector& EV = *eVs[0]; + for (i = 1;i <= N2.size();i++) + P += EV(i*nf)*N2(i); + + // Add pressure to strain + for (k = 1;k <= nsd;k++) + sigma(k,k) -= P; + } + + return true; +} + + bool StokesNorm::initElement (const std::vector& 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 + const Vector& N, const Matrix& dNdX, + const Vec3& X) const { int i, k, l; double value, eps, epsh; @@ -221,7 +356,7 @@ bool StokesNorm::evalInt (LocalIntegral*& elmInt, double detJW, Vec3 Uh; Uh = 0.0; for (i = 1;i <= N.size();i++) for (k = 1;k <= nsd;k++) - Uh[k-1] += EV((i-1)*nf + k)*N(i); + Uh[k-1] += EV((i-1)*nf + k)*N(i); // L2-norm of velocity error U -= Uh; @@ -237,7 +372,7 @@ bool StokesNorm::evalInt (LocalIntegral*& elmInt, double detJW, Vec3 dPh; dPh = 0.0; for (i = 1;i <= N.size();i++) for (k = 1;k <= nsd;k++) - dPh[k-1] += EV(i*nf)*dNdX(i,k); + dPh[k-1] += EV(i*nf)*dNdX(i,k); // H1-seminorm of pressure dP -= dPh; @@ -251,73 +386,102 @@ bool StokesNorm::evalInt (LocalIntegral*& elmInt, double detJW, // Analytical velocity gradient Tensor gradU = (*anasol->getVectorSecSol())(X); + + // Numerical velocity gradient + Tensor gradUh(nsd); + problem.strain(dNdX,gradUh); - // Computed velocity gradient - Matrix gradUh(nsd,nsd); - for (i = 1;i <= N.size();i++) - for (k = 1;k <= nsd;k++) - for (l = 1;l <= nsd;l++) - gradUh(k,l) += EV((i-1)*nf + k)*dNdX(i,l); - - if (problem.getFormulation() == Stokes::LAPLACE) { - for (k = 1;k <= nsd;k++) - for (l = 1;l <= nsd;l++) { - // Energy norm of analytical solution - pnorm[ip++] += mu*gradU(k,l)*gradU(k,l)*detJW; - // Energy norm of numerical solution - pnorm[ip++] += mu*gradUh(k,l)*gradUh(k,l)*detJW; - // Energy norm of error - value = gradU(k,l)-gradUh(k,l); - pnorm[ip++] += mu*value*value*detJW; - } + if (problem.getFormulation() == Stokes::STRESS) { + pnorm[ip++] += 2.0*mu*gradU.innerProd(gradU)*detJW; + pnorm[ip++] += 2.0*mu*gradUh.innerProd(gradUh)*detJW; + gradUh *= -1.0; + gradU += gradUh; + pnorm[ip++] += 2.0*mu*gradU.innerProd(gradU)*detJW; } else { - for (k = 1;k <= nsd;k++) - for (l = 1;l <= nsd;l++) { - // Strain of analytical solution - eps = 0.5*(gradU(k,l) + gradU(l,k)); - // Strain of computed solution - epsh = 0.5*(gradUh(k,l) + gradUh(l,k)); - // Energy norm of analytical solution - pnorm[ip++] += mu*eps*eps*detJW; - // Energy norm of computed solution - pnorm[ip++] += mu*epsh*epsh*detJW; - // Energy norm of error - eps -= epsh; - pnorm[ip++] += mu*epsh*epsh*detJW; - } + pnorm[ip++] += mu*gradU.innerProd(gradU)*detJW; + pnorm[ip++] += mu*gradUh.innerProd(gradUh)*detJW; + gradUh *= -1.0; + gradU += gradUh; + pnorm[ip++] += mu*gradU.innerProd(gradU)*detJW; } + } + + return true; +} + + +bool StokesForce::initElementBou(const std::vector& 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() == Stokes::STRESS) + sigma += sigma.transpose(); + + // Analytical stress + sigma *= mu; + for (i = 1;i <= nsd;i++) + sigma(i,i) -= P; + + // Analytical traction + Vec3 t = sigma*normal; + + for (i = 0;i < nsd;i++) + pnorm[ip++] += t[i]*detJW; + + // Error in traction + t -= th; + pnorm[ip++] += t*t*detJW; } return true; } -bool Stokes::getBodyForce(const Vec3& X, Vector& f) const +size_t StokesForce::getNoFields() const { - const Vec4* Y = dynamic_cast(&X); - if (!Y) return false; - - const double PI = 3.141592653589793238462; - const double x = X[0]; - const double y = X[1]; - const double t = Y->t; - - f.fill(0.0); + size_t nsd = problem.getNoSpaceDim(); + + if (anasol) + return 2*nsd + 1; + else + return nsd; +} - f(1) = rho*pow(sin(PI*x),2.0)*sin(2.0*PI*y)*cos(t) + - rho*PI*pow(sin(PI*x)*sin(2.0*PI*y)*sin(t),2.0)*sin(2.0*PI*x) - - 2.0*rho*PI*pow(sin(PI*x)*sin(PI*y)*sin(t),2.0)*sin(2.0*PI*x)*cos(2.0*PI*y) - - PI*sin(PI*x)*sin(PI*y)*sin(t) - - 2.0*pow(PI,2.0)*mu*sin(2.0*PI*y)*sin(t)*(cos(2*PI*x)-2*pow(sin(PI*x),2.0)); - - f(2) = -rho*sin(2.0*PI*x)*pow(sin(PI*y),2.0)*cos(t) - - 2.0*rho*PI*pow(sin(PI*x)*sin(PI*y)*sin(t),2.0)*cos(2.0*PI*x)*sin(2.0*PI*y) + - rho*PI*pow(sin(2.0*PI*x)*sin(PI*y)*sin(t),2.0)*sin(2.0*PI*y) + - PI*cos(PI*x)*cos(PI*y)*sin(t) - - 4*mu*pow(PI*sin(PI*y),2.0)*sin(2.0*PI*x)*sin(t) + - 2.0*mu*pow(PI,2.0)*sin(2.0*PI*x)*cos(2.0*PI*y)*sin(t); - - return true; -} diff --git a/src/Integrands/Stokes.h b/src/Integrands/Stokes.h index 0bc51b47..6da47e7a 100644 --- a/src/Integrands/Stokes.h +++ b/src/Integrands/Stokes.h @@ -63,9 +63,29 @@ public: //! \brief Initializes current element for numerical integration. //! \param[in] MNPC Matrix of nodal point correspondance for current element virtual bool initElement(const std::vector& 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 @@ -106,6 +126,12 @@ public: //! \param[in] asol Pointer to analytical solution field (optional) virtual NormBase* getNormIntegrand(AnaSol* asol = 0) const; + //! \brief Returns a pointer to an Integrand for boundary force evaluation. + //! \note The Integrand is allocated dynamically and has to be deleted + //! manually when leaving the scope of the pointer returned. + //! \param[in] asol Pointer to analytical solution field (optional) + virtual NormBase* getForceIntegrand(AnaSol* asol = 0) const; + //! \brief Returns problem formulation type. SIM::Formulation getFormulation() const { return formulation; } @@ -131,6 +157,9 @@ public: //! \brief Returns a pointer to current element solution vector. virtual const Vector* getElementPressure() const { return 0; } + + //! \brief Returns the number of space dimensions + size_t getNoSpaceDim() const { return nsd; } //! \brief Returns the number of solution fields. virtual size_t getNoFields() const { return nsd; } @@ -149,6 +178,25 @@ public: virtual const char* getSecFieldLabel(size_t i, const char* prefix = 0) const { return 0; } + //! \brief Calculates viscous part of stress tensor at current point. + //! \param[in] dNdX Basis function gradients at current point + //! \param[out] eps Strain tensor at current point + virtual bool strain(const Matrix& dNdX, Tensor& eps) const; + + //! \brief Calculates the (Cauchy) stress tensor at current point + //! \param[in] N Basis functions at current point + //! \param[in] dNdX Basis function gradients at current point + //! \param[out] sigma Strain tensor at current point + bool stress(const Vector& N, const Matrix& dNdX, Tensor& sigma) const; + + //! \brief Calculates the (Cauchy) stress tensor at current point + //! \param[in] N1 Velocity basis functions at current point + //! \param[in] dNdX Basis function gradients at current point + //! \param[out] sigma Strain tensor at current point + bool stress(const Vector& N1, const Vector& N2, + const Matrix& dN1dX, const Matrix& dN2dX, + Tensor& sigma) const; + protected: //! \brief Utility used by the virtual \a evalInt and \a evalBou methods. //! \param elmInt Pointer to the integrated element quantities @@ -195,6 +243,9 @@ public: //! \brief Initializes current element for numerical integration. //! \param[in] MNPC Matrix of nodal point correspondance for current element virtual bool initElement(const std::vector& 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 @@ -218,4 +269,81 @@ protected: }; +/*! + \brief Class representing for computing boundary force +*/ + +class StokesForce : public NormBase +{ +public: + //! \brief The only constructor initializes its data members. + //! \param[in] p The Stokes problem to evaluate norms for + //! \param[in] a The analytical velocity and pressure fields (optional) + StokesForce(Stokes& p, AnaSol* a = 0) : problem(p), anasol(a) {} + //! \brief Empty destructor. + virtual ~StokesForce() {} + + //! \brief Initializes current element for boundary integration. + //! \param[in] MNPC Matrix of nodal point correspondance for current element + virtual bool initElementBou(const std::vector& 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