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

Some modifications regarding analytical solution fields, we need to distingush between Tensor fields and Symmetric Tensor fields, mostly because the Function classes are template-based.

Browse Source 
git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@798 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
kmo 2011-02-15 07:51:36 +00:00 committed by Knut Morten Okstad
parent aedb74ee5e
commit da2d5ad2ac
27 changed files with 497 additions and 460 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

50
Apps/FiniteDefElasticity/NonlinearElasticityTL.C
View File

@ -1,4 +1,4 @@
// $Id: NonlinearElasticityTL.C,v 1.4 2011-02-08 09:06:02 kmo Exp $
// $Id$
//==============================================================================
//!
//! \file NonlinearElasticityTL.C
@ -80,6 +80,54 @@ void NonlinearElasticityTL::setMode (SIM::SolutionMode mode)
}
bool NonlinearElasticityTL::evalInt (LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X) const
{
// Evaluate the deformation gradient, F, and the Green-Lagrange strains, E,
// and compute the nonlinear strain-displacement matrix B from dNdX and F
SymmTensor E(nsd), S(nsd);
if (!this->kinematics(dNdX,E))
return false;
// Evaluate the constitutive relation
double U = 0.0;
bool lHaveStrains = !E.isZero();
if (eKm || eKg || iS)
if (!this->constitutive(Cmat,S,U,E,X, (eKg || iS) && lHaveStrains))
return false;
if (eKm)
{
// Integrate the material stiffness matrix
CB.multiply(Cmat,Bmat).multiply(detJW); // CB = C*B*|J|*w
eKm->multiply(Bmat,CB,true,false,true); // EK += B^T * CB
}
if (eKg && lHaveStrains)
// Integrate the geometric stiffness matrix
this->formKG(*eKg,dNdX,S,detJW);
if (eM)
// Integrate the mass matrix
this->formMassMatrix(*eM,N,X,detJW);
if (iS && lHaveStrains)
{
// Integrate the internal forces
S *= -detJW;
if (!Bmat.multiply(S,*iS,true,true)) // ES -= B^T*S
return false;
}
if (eS)
// Integrate the load vector due to gravitation and other body forces
this->formBodyForce(*eS,N,X,detJW);
return this->getIntegralResult(elmInt);
}
bool NonlinearElasticityTL::evalBou (LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X, const Vec3& normal) const

17
Apps/FiniteDefElasticity/NonlinearElasticityTL.h
View File

@ -1,4 +1,4 @@
// $Id: NonlinearElasticityTL.h,v 1.2 2011-02-08 09:06:02 kmo Exp $
// $Id$
//==============================================================================
//!
//! \file NonlinearElasticityTL.h
@ -43,6 +43,16 @@ public:
//! \param[in] mode The solution mode to use
virtual void setMode(SIM::SolutionMode mode);
//! \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 a boundary point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] detJW Jacobian determinant times integration point weight
@ -69,8 +79,9 @@ protected:
virtual bool kinematics(const Matrix& dNdX, SymmTensor& E) const;
protected:
bool formB; //!< Flag determining whether we need to form the B-matrix
mutable Matrix F; //!< Deformation gradient
bool formB; //!< Flag determining whether we need to form the B-matrix
mutable Matrix F; //!< Deformation gradient
mutable Matrix CB; //!< Result of the matrix-matrix product C*B
};
#endif

4
Apps/FiniteDefElasticity/NonlinearElasticityUL.C
View File

@ -1,4 +1,4 @@
// $Id: NonlinearElasticityUL.C,v 1.5 2011-02-08 09:06:02 kmo Exp $
// $Id$
//==============================================================================
//!
//! \file NonlinearElasticityUL.C
@ -428,7 +428,7 @@ bool NonlinearElasticityUL::constitutive (Matrix& C, SymmTensor& sigma,
}
NormBase* NonlinearElasticityUL::getNormIntegrand (TensorFunc*) const
NormBase* NonlinearElasticityUL::getNormIntegrand (AnaSol*) const
{
return new ElasticityNormUL(*const_cast<NonlinearElasticityUL*>(this));
}

6
Apps/FiniteDefElasticity/NonlinearElasticityUL.h
View File

@ -1,4 +1,4 @@
// $Id: NonlinearElasticityUL.h,v 1.7 2011-02-08 09:06:02 kmo Exp $
// $Id$
//==============================================================================
//!
//! \file NonlinearElasticityUL.h
@ -89,7 +89,7 @@ public:
//! \note The Integrand object is allocated dynamically and has to be deleted
//! manually when leaving the scope of the pointer variable receiving the
//! returned pointer value.
virtual NormBase* getNormIntegrand(TensorFunc* = 0) const;
virtual NormBase* getNormIntegrand(AnaSol* = 0) const;
//! \brief Calculates some kinematic quantities at current point.
//! \param[in] dNdX Basis function gradients at current point
@ -143,7 +143,7 @@ class ElasticityNormUL : public ElasticityNorm
public:
//! \brief The only constructor initializes its data members.
//! \param[in] p The linear elasticity problem to evaluate norms for
ElasticityNormUL(NonlinearElasticityUL& p) : ElasticityNorm(p,0) {}
ElasticityNormUL(NonlinearElasticityUL& p) : ElasticityNorm(p) {}
//! \brief Empty destructor.
virtual ~ElasticityNormUL() {}

4
Apps/FiniteDefElasticity/NonlinearElasticityULMX.C
View File

@ -1,4 +1,4 @@
// $Id: NonlinearElasticityULMX.C,v 1.7 2011-02-08 09:06:02 kmo Exp $
// $Id$
//==============================================================================
//!
//! \file NonlinearElasticityULMX.C
@ -459,7 +459,7 @@ bool NonlinearElasticityULMX::finalizeElement (LocalIntegral*& elmInt)
}
NormBase* NonlinearElasticityULMX::getNormIntegrand (TensorFunc*) const
NormBase* NonlinearElasticityULMX::getNormIntegrand (AnaSol*) const
{
return new ElasticityNormULMX(*const_cast<NonlinearElasticityULMX*>(this));
}

4
Apps/FiniteDefElasticity/NonlinearElasticityULMX.h
View File

@ -1,4 +1,4 @@
// $Id: NonlinearElasticityULMX.h,v 1.3 2011-01-08 16:29:10 kmo Exp $
// $Id$
//==============================================================================
//!
//! \file NonlinearElasticityULMX.h
@ -89,7 +89,7 @@ public:
//! \note The Integrand object is allocated dynamically and has to be deleted
//! manually when leaving the scope of the pointer variable receiving the
//! returned pointer value.
virtual NormBase* getNormIntegrand(TensorFunc* = 0) const;
virtual NormBase* getNormIntegrand(AnaSol* = 0) const;
//! \brief Returns which integrand to be used.
virtual int getIntegrandType() const { return 4; }

40
Apps/SIMLinEl2D.C
View File

@ -1,4 +1,4 @@
// $Id: SIMLinEl2D.C,v 1.19 2011-02-09 10:07:36 rho Exp $
// $Id$
//==============================================================================
//!
//! \file SIMLinEl2D.C
@ -13,13 +13,14 @@
#include "SIMLinEl2D.h"
#include "LinearElasticity.h"
#include "NonlinearElasticity.h"
#include "FiniteDefElasticity/NonlinearElasticity.h"
#include "AnalyticSolutions.h"
#include "Functions.h"
#include "Utilities.h"
#include "Vec3Oper.h"
#include "Property.h"
#include "Tensor.h"
#include "AnaSol.h"
#include <string.h>
@ -34,14 +35,6 @@ SIMLinEl2D::SIMLinEl2D (int form, bool planeStress)
myProblem = new NonlinearElasticity(2,planeStress);
break;
}
asol = 0;
}
SIMLinEl2D::~SIMLinEl2D ()
{
if (asol) delete asol;
}
@ -55,7 +48,7 @@ bool SIMLinEl2D::parse (char* keyWord, std::istream& is)
{
double gx = atof(strtok(keyWord+7," "));
double gy = atof(strtok(NULL," "));
elp->setGravity(gx,gy,0.0);
elp->setGravity(gx,gy);
if (myPid == 0)
std::cout <<"\nGravitation vector: "
<< gx <<" "<< gy << std::endl;
@ -78,7 +71,7 @@ bool SIMLinEl2D::parse (char* keyWord, std::istream& is)
<< E <<" "<< nu <<" "<< rho << std::endl;
if (code == 0 || i == 0)
elp->setMaterial(E,nu,rho);
for (unsigned int j = 0; j < myProps.size() && code > 0; j++)
for (size_t j = 0; j < myProps.size() && code > 0; j++)
if (myProps[j].pindx == (size_t)code &&
myProps[j].pcode == Property::UNDEFINED)
{
@ -118,7 +111,7 @@ bool SIMLinEl2D::parse (char* keyWord, std::istream& is)
double a = atof(strtok(NULL," "));
double F0 = atof(strtok(NULL," "));
double nu = atof(strtok(NULL," "));
asol = new AnaSol(NULL,NULL,NULL,new Hole(a,F0,nu));
mySol = new AnaSol(new Hole(a,F0,nu));
std::cout <<"\nAnalytical solution: Hole a="<< a <<" F0="<< F0
<<" nu="<< nu << std::endl;
}
@ -127,7 +120,7 @@ bool SIMLinEl2D::parse (char* keyWord, std::istream& is)
double a = atof(strtok(NULL," "));
double F0 = atof(strtok(NULL," "));
double nu = atof(strtok(NULL," "));
asol = new AnaSol(NULL,NULL,NULL,new Lshape(a,F0,nu));
mySol = new AnaSol(new Lshape(a,F0,nu));
std::cout <<"\nAnalytical solution: Lshape a="<< a <<" F0="<< F0
<<" nu="<< nu << std::endl;
}
@ -136,7 +129,7 @@ bool SIMLinEl2D::parse (char* keyWord, std::istream& is)
double L = atof(strtok(NULL," "));
double H = atof(strtok(NULL," "));
double F0 = atof(strtok(NULL," "));
asol = new AnaSol(NULL,NULL,NULL,new CanTS(L,H,F0));
mySol = new AnaSol(new CanTS(L,H,F0));
std::cout <<"\nAnalytical solution: CanTS L="<< L <<" H="<< H
<<" F0="<< F0 << std::endl;
}
@ -144,7 +137,7 @@ bool SIMLinEl2D::parse (char* keyWord, std::istream& is)
{
double H = atof(strtok(NULL," "));
double M0 = atof(strtok(NULL," "));
asol = new AnaSol(NULL,NULL,NULL,new CanTM(H,M0));
mySol = new AnaSol(new CanTM(H,M0));
std::cout <<"\nAnalytical solution: CanTM H="<< H
<<" M0="<< M0 << std::endl;
}
@ -154,21 +147,24 @@ bool SIMLinEl2D::parse (char* keyWord, std::istream& is)
double b = atof(strtok(NULL," "));
double u0 = atof(strtok(NULL," "));
double E = atof(strtok(NULL," "));
asol = new AnaSol(NULL,NULL,NULL,new CurvedBeam(u0,a,b,E));
mySol = new AnaSol(new CurvedBeam(u0,a,b,E));
std::cout <<"\nAnalytical solution: Curved Beam a="<< a <<" b="<< b
<<" u0="<< u0 <<" E="<< E << std::endl;
}
else
{
std::cerr <<" ** SIMLinEl2D::parse: Unknown analytical solution "
<< cline << std::endl;
<< cline <<" (ignored)"<< std::endl;
return true;
}
// Define the analytical boundary traction field
int code = (cline = strtok(NULL," ")) ? atoi(cline) : 0;
if (code > 0 && asol->getVectorSecSol())
if (code > 0 && mySol->getStressSol())
{
std::cout <<"Pressure code "<< code <<": Analytical traction"<< std::endl;
this->setPropertyType(code,Property::NEUMANN);
myTracs[code] = new TractionField(*(asol->getVectorSecSol()));
myTracs[code] = new TractionField(*mySol->getStressSol());
}
}
@ -201,11 +197,11 @@ bool SIMLinEl2D::parse (char* keyWord, std::istream& is)
return false;
}
if (asol->getVectorSecSol())
if (mySol && mySol->getStressSol())
{
std::cout <<"\tTraction on P"<< press.patch
<<" E"<< (int)press.lindx << std::endl;
myTracs[1+i] = new TractionField(*asol->getVectorSecSol());
myTracs[1+i] = new TractionField(*mySol->getStressSol());
}
else
{

11
Apps/SIMLinEl2D.h
View File

@ -1,4 +1,4 @@
// $Id: SIMLinEl2D.h,v 1.10 2011-02-09 10:07:31 rho Exp $
// $Id$
//==============================================================================
//!
//! \file SIMLinEl2D.h
@ -16,7 +16,6 @@
#include "SIM2D.h"
#include "SIMenums.h"
#include "AnaSol.h"
/*!
@ -47,8 +46,8 @@ public:
//! \param[in] form Problem formulation option
//! \param[in] planeStress Plane stress/plane strain option
SIMLinEl2D(int form = SIM::LINEAR, bool planeStress = true);
//! \brief The destructor frees the dynamically allocated data.
virtual ~SIMLinEl2D();
//! \brief Empty destructor.
virtual ~SIMLinEl2D() {}
protected:
//! \brief Parses a data section from the input stream.
@ -64,12 +63,8 @@ protected:
//! \param[in] propInd Physical property index
virtual bool initNeumann(size_t propInd);
//! \brief Returns the analytical solution function, if any.
virtual AnaSol* getAnaSol() const { return asol; }
private:
MaterialVec mVec; //!< Material data
AnaSol* asol; //!< Analytical stress field
};
#endif

39
Apps/SIMLinEl3D.C
View File

@ -1,4 +1,4 @@
// $Id: SIMLinEl3D.C,v 1.29 2011-02-09 10:08:02 rho Exp $
// $Id$
//==============================================================================
//!
//! \file SIMLinEl3D.C
@ -13,13 +13,14 @@
#include "SIMLinEl3D.h"
#include "LinearElasticity.h"
#include "NonlinearElasticity.h"
#include "FiniteDefElasticity/NonlinearElasticity.h"
#include "AnalyticSolutions.h"
#include "Functions.h"
#include "Utilities.h"
#include "Vec3Oper.h"
#include "Property.h"
#include "Tensor.h"
#include "AnaSol.h"
#include <string.h>
@ -128,7 +129,7 @@ public:
T(1,i) = v1[i-1];
T(2,i) = v2[i-1];
T(3,i) = v3[i-1];
}
}
#ifdef PRINT_CS
s1 << v1 <<'\n';
s2 << v2 <<'\n';
@ -161,14 +162,6 @@ SIMLinEl3D::SIMLinEl3D (bool checkRHS, int form) : SIM3D(checkRHS)
myProblem = new NonlinearElasticity();
break;
}
asol = 0;
}
SIMLinEl3D::~SIMLinEl3D ()
{
if (asol) delete asol;
}
@ -208,8 +201,9 @@ bool SIMLinEl3D::parse (char* keyWord, std::istream& is)
<< E <<" "<< nu <<" "<< rho << std::endl;
if (code == 0 || i == 0)
elp->setMaterial(E,nu,rho);
for (unsigned int j = 0; j < myProps.size() && code > 0; j++)
if (myProps[j].pindx == code && myProps[j].pcode == Property::UNDEFINED)
for (size_t j = 0; j < myProps.size() && code > 0; j++)
if (myProps[j].pindx == (size_t)code &&
myProps[j].pcode == Property::UNDEFINED)
{
myProps[j].pindx = mVec.size();
myProps[j].pcode = Property::MATERIAL;
@ -232,7 +226,7 @@ bool SIMLinEl3D::parse (char* keyWord, std::istream& is)
double a = atof(strtok(NULL," "));
double F0 = atof(strtok(NULL," "));
double nu = atof(strtok(NULL," "));
asol = new AnaSol(NULL,NULL,NULL,new Hole(a,F0,nu,true));
mySol = new AnaSol(new Hole(a,F0,nu,true));
std::cout <<"\nAnalytical solution: Hole a="<< a <<" F0="<< F0
<<" nu="<< nu << std::endl;
}
@ -241,7 +235,7 @@ bool SIMLinEl3D::parse (char* keyWord, std::istream& is)
double a = atof(strtok(NULL," "));
double F0 = atof(strtok(NULL," "));
double nu = atof(strtok(NULL," "));
asol = new AnaSol(NULL,NULL,NULL,new Lshape(a,F0,nu,true));
mySol = new AnaSol(new Lshape(a,F0,nu,true));
std::cout <<"\nAnalytical solution: Lshape a="<< a <<" F0="<< F0
<<" nu="<< nu << std::endl;
}
@ -250,21 +244,24 @@ bool SIMLinEl3D::parse (char* keyWord, std::istream& is)
double L = atof(strtok(NULL," "));
double H = atof(strtok(NULL," "));
double F0 = atof(strtok(NULL," "));
asol = new AnaSol(NULL,NULL,NULL,new CanTS(L,H,F0,true));
mySol = new AnaSol(new CanTS(L,H,F0,true));
std::cout <<"\nAnalytical solution: CanTS L="<< L <<" H="<< H
<<" F0="<< F0 << std::endl;
}
else
{
std::cerr <<" ** SIMLinEl3D::parse: Unknown analytical solution "
<< cline << std::endl;
<< cline <<" (ignored)"<< std::endl;
return true;
}
// Define the analytical boundary traction field
int code = (cline = strtok(NULL," ")) ? atoi(cline) : 0;
if (code > 0 && asol->getVectorSecSol())
if (code > 0 && mySol->getStressSol())
{
std::cout <<"Pressure code "<< code <<": Analytical traction"<< std::endl;
this->setPropertyType(code,Property::NEUMANN);
myTracs[code] = new TractionField(*(asol->getVectorSecSol()));
myTracs[code] = new TractionField(*mySol->getStressSol());
}
}
@ -297,11 +294,11 @@ bool SIMLinEl3D::parse (char* keyWord, std::istream& is)
return false;
}
if (asol && asol->getVectorSecSol())
if (mySol && mySol->getStressSol())
{
std::cout <<"\tTraction on P"<< press.patch
<<" F"<< (int)press.lindx << std::endl;
myTracs[1+i] = new TractionField(*asol->getVectorSecSol());
myTracs[1+i] = new TractionField(*mySol->getStressSol());
}
else
{

11
Apps/SIMLinEl3D.h
View File

@ -1,4 +1,4 @@
// $Id: SIMLinEl3D.h,v 1.12 2011-02-09 10:07:57 rho Exp $
// $Id$
//==============================================================================
//!
//! \file SIMLinEl3D.h
@ -16,7 +16,6 @@
#include "SIM3D.h"
#include "SIMenums.h"
#include "AnaSol.h"
/*!
@ -47,8 +46,8 @@ public:
//! \param[in] checkRHS If \e true, ensure the model is in a right-hand system
//! \param[in] form Problem formulation option
SIMLinEl3D(bool checkRHS = false, int form = SIM::LINEAR);
//! \brief The destructor frees the dynamically allocated data.
virtual ~SIMLinEl3D();
//! \brief Empty destructor.
virtual ~SIMLinEl3D() {}
protected:
//! \brief Parses a data section from the input stream.
@ -64,12 +63,8 @@ protected:
//! \param[in] propInd Physical property index
virtual bool initNeumann(size_t propInd);
//! \brief Returns the analytical solution function, if any.
virtual AnaSol* getAnaSol() const { return asol; }
private:
MaterialVec mVec; //!< Material data
AnaSol* asol; //!< Analytical solution fields
};
#endif

94
src/ASM/IntegrandBase.h
View File

@ -1,4 +1,4 @@
// $Id: IntegrandBase.h,v 1.20 2011-02-08 12:10:25 rho Exp $
// $Id$
//==============================================================================
//!
//! \file IntegrandBase.h
@ -121,7 +121,6 @@ public:
//! the same element, provided the necessary integration point values are
//! stored internally in the object during the first integration loop.
virtual bool finalizeElement(LocalIntegral*&) { return true; }
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
@ -294,19 +293,6 @@ public:
return false;
}
//! \brief Evaluates the secondary solution at a result point.
//! \param[out] s The 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
virtual bool evalSecSol(Vector& s,
const Vector& N, const Matrix& dNdX,
const Vec3& X, const std::vector<int>& MNPC) const
{
return false;
}
//! \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
@ -325,44 +311,56 @@ public:
return this->evalSol(s,N1,dN1dX,X,MNPC1);
}
//! \brief Evaluates the analytical solution at an integration point.
//! \brief Evaluates the analytical primary solution at a result point.
//! \param[out] s The solution field values at current point
//! \param[in] asol The analytical solution field (vector field)
//! \param[in] X Cartesian coordinates of current point
virtual bool evalPrimSol(Vector& s, const VecFunc& asol, const Vec3& X) const
{
std::cerr <<" *** Integrand::evalPrimSol not implemented"<< std::endl;
return false;
}
//! \brief Evaluates the analytical secondary solution at a result point.
//! \param[out] s The solution field values at current point
//! \param[in] asol The analytical solution field (tensor field)
//! \param[in] X Cartesian coordinates of current point
virtual bool evalSol(Vector& s, const TensorFunc& asol, const Vec3& X) const
{
std::cerr <<" *** Integrand::evalSol (exact) not implemented"<< std::endl;
return false;
}
//! \brief Evaluates the analytical secondary solution at a result point.
//! \param[out] s The solution field values at current point
//! \param[in] asol The analytical solution field (symmetric tensor field)
//! \param[in] X Cartesian coordinates of current point
virtual bool evalSol(Vector& s, const STensorFunc& asol, const Vec3& X) const
{
std::cerr <<" *** Integrand::evalSol (exact) not implemented"<< std::endl;
return false;
}
//! \brief Evaluates the analytical primary solution at a result point.
//! \param[out] s The solution field value at current point
//! \param[in] asol The analytical solution field (scalar field)
//! \param[in] X Cartesian coordinates of current point
virtual bool evalPrimSol(real& s, const RealFunc& asol, const Vec3& X) const
{
std::cerr <<" *** Integrand::evalPrimSol not implemented"<< std::endl;
return false;
}
//! \brief Evaluates the analytical secondary solution at a result point.
//! \param[out] s The solution field values at current point
//! \param[in] asol The analytical solution field (vector field)
//! \param[in] X Cartesian coordinates of current point
virtual bool evalSol(Vector& s, const VecFunc& asol, const Vec3& X) const
{
std::cerr <<" *** Integrand::evalSol (exact) not implemented"<< std::endl;
return false;
}
//! \brief Evaluates the analytical secondary solution at an integration point.
//! \param[out] s The solution field values at current point
//! \param[in] asol The analytical solution field (tensor field)
//! \param[in] X Cartesian coordinates of current point
virtual bool evalSecSol(Vector& s, const TensorFunc& asol, const Vec3& X) const
{
return false;
}
//! \brief Evaluates the analytical scalar solution at an integration point.
//! \param[out] s The solution field values at current point
//! \param[in] asol The analytical solution field (vector field)
//! \param[in] X Cartesian coordinates of current point
virtual bool evalSolScal(real& s, const RealFunc& asol, const Vec3& X) const
{
std::cerr <<" *** Integrand::evalSol (exact) not implemented"<< std::endl;
return false;
}
//! \brief Evaluates the analytical secondary scalar solution at an integration point.
//! \param[out] s The solution field values at current point
//! \param[in] asol The analytical solution field (vector field)
//! \param[in] X Cartesian coordinates of current point
virtual bool evalSecSolScal(Vector& s, const VecFunc& asol, const Vec3& X) const
{
return false;
}
//! \brief Writes surface tractions/fluxes 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
@ -390,15 +388,9 @@ public:
}
//! \brief Returns a pointer to an Integrand for solution norm evaluation.
//! \param[in] asol Pointer to the analytical solution field (optional)
//! \param[in] asol Pointer to the analytical solution (optional)
virtual NormBase* getNormIntegrand(AnaSol* asol = 0) const { return 0; }
//! \brief Returns a pointer to an Integrand for solution norm evaluation.
//! \param[in] asol Pointer to the analytical solution field (optional)
//!
//! \details This version is used for scalar problems only.
virtual NormBase* getNormIntegrandScal(AnaSol* asol = 0) const { return 0; }
//! \brief Returns the number of solution vectors.
size_t getNoSolutions() const { return primsol.size(); }

100
src/Integrands/AnaSol.h Normal file
View File

@ -0,0 +1,100 @@
// $Id$
//==============================================================================
//!
//! \file AnaSol.h
//!
//! \date Dec 7 2010
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Analytical solution fields (primary and secondary).
//!
//==============================================================================
#ifndef _ANA_SOL_H
#define _ANA_SOL_H
#include "Function.h"
/*!
\brief Class for analytical solution fields (primary and secondary solution).
*/
class AnaSol
{
protected:
RealFunc* scalSol; //!< Primary scalar solution field
VecFunc* scalSecSol; //!< Secondary scalar solution field (gradient field)
VecFunc* vecSol; //!< Primary vector solution field
TensorFunc* vecSecSol; //!< Secondary solution field (vector gradient field)
STensorFunc* stressSol; //!< Secondary solution field (stress field)
public:
//! \brief Default constructor initializing all solution fields.
//! \param[in] s1 Primary scalar solution field
//! \param[in] s2 Secondary solution field, gradient
//! \param[in] v1 Primary vector solution field
//! \param[in] v2 Secondary solution field, gradient
//! \param[in] v3 Secondary solution field, symmetric stress tensor
//!
//! \details It is assumed that all the arguments are pointers to dynamically
//! allocation objects, as the class destructor will attempt to delete them.
AnaSol(RealFunc* s1 = 0, VecFunc* s2 = 0,
VecFunc* v1 = 0, TensorFunc* v2 = 0, STensorFunc* v3 = 0)
: scalSol(s1), scalSecSol(s2), vecSol(v1), vecSecSol(v2), stressSol(v3) {}
//! \brief Constructor initializing the symmetric stress tensor field only.
//! \param[in] sigma Symmetric stress tensor field
AnaSol(STensorFunc* sigma)
: scalSol(0), scalSecSol(0), vecSol(0), vecSecSol(0), stressSol(sigma) {}
//! \brief The destructor frees the analytical solution fields.
virtual ~AnaSol()
{
if (scalSol) delete scalSol;
if (scalSecSol) delete scalSecSol;
if (vecSol) delete vecSol;
if (vecSecSol) delete vecSecSol;
if (stressSol) delete stressSol;
}
//! \brief Checks whether a scalar solution is defined.
char hasScalarSol() const
{
if (scalSecSol)
return 2;
else if (scalSol)
return 1;
else
return 0;
}
//! \brief Checks whether a vector solution is defined.
char hasVectorSol() const
{
if (stressSol)
return 3;
else if (vecSecSol)
return 2;
else if (vecSol)
return 1;
else
return 0;
}
//! \brief Returns the scalar solution, if any.
RealFunc*getScalarSol() const { return scalSol; }
//! \brief Returns the secondary scalar solution, if any.
VecFunc* getScalarSecSol() const { return scalSecSol; }
//! \brief Returns the vector solution, if any.
VecFunc* getVectorSol() const { return vecSol; }
//! \brief Returns the secondary vector solution, if any.
TensorFunc* getVectorSecSol() const { return vecSecSol; }
//! \brief Returns the stress solution, if any.
STensorFunc* getStressSol() const { return stressSol; }
};
#endif

14
src/Integrands/AnalyticSolutions.C
View File

@ -1,4 +1,4 @@
// $Id: AnalyticSolutions.C,v 1.11 2011-02-08 12:19:52 rho Exp $
// $Id$
//==============================================================================
//!
//! \file AnalyticSolutions.C
@ -13,8 +13,6 @@
#include "AnalyticSolutions.h"
#include "Vec3.h"
#include "Vec3Oper.h"
#include <math.h>
/*!
@ -26,7 +24,7 @@
IJNME, 33, 1331-1364, 1992, page 1355 (eq. 36).
*/
Tensor Hole::evaluate (const Vec3& X) const
SymmTensor Hole::evaluate (const Vec3& X) const
{
double R = hypot(X.x,X.y);
double th = atan2(X.y,X.x);
@ -82,7 +80,7 @@ Lshape::Lshape (double r, double f, double P, bool use3D)
}
Tensor Lshape::evaluate (const Vec3& X) const
SymmTensor Lshape::evaluate (const Vec3& X) const
{
// Some constants (see Szabo & Babuska for an elaboration)
const double lambda = 0.544483737;
@ -112,7 +110,7 @@ Tensor Lshape::evaluate (const Vec3& X) const
}
Tensor CanTS::evaluate (const Vec3& X) const
SymmTensor CanTS::evaluate (const Vec3& X) const
{
double x = X.x/L;
double y = (is3D ? X.z : X.y)/H - 0.5;
@ -127,7 +125,7 @@ Tensor CanTS::evaluate (const Vec3& X) const
}
Tensor CanTM::evaluate (const Vec3& X) const
SymmTensor CanTM::evaluate (const Vec3& X) const
{
double y = (is3D ? X.z : X.y)/H - 0.5;
double I = H*H*H / 12.0;
@ -157,7 +155,7 @@ CurvedBeam::CurvedBeam (double u0, double Ri, double Ro, double E, bool use3D)
}
Tensor CurvedBeam::evaluate (const Vec3& X) const
SymmTensor CurvedBeam::evaluate (const Vec3& X) const
{
// Find polar coordinates
double r = hypot(X.x,X.y);

34
src/Integrands/AnalyticSolutions.h
View File

@ -1,4 +1,4 @@
// $Id: AnalyticSolutions.h,v 1.10 2011-02-08 12:19:37 rho Exp $
// $Id$
//==============================================================================
//!
//! \file AnalyticSolutions.h
@ -22,7 +22,7 @@
\brief Analytic solution for an infinite plate with a hole.
*/
class Hole : public TensorFunc
class Hole : public STensorFunc
{
public:
//! \brief Constructor with some default parameters.
@ -33,13 +33,13 @@ public:
protected:
//! \brief Evaluates the analytic stress tensor at the point \a x.
virtual Tensor evaluate(const Vec3& x) const;
virtual SymmTensor evaluate(const Vec3& x) const;
private:
bool is3D; //!< Flag telling whether to return a 3D stress tensor or not
double a; //!< Hole radius
double F0; //!< Load factor
double nu; //!< Poisson's ratio
bool is3D; //!< Flag telling whether to return a 3D stress tensor or not
};
@ -47,7 +47,7 @@ private:
\brief Analytic solution for the L-shaped domain.
*/
class Lshape : public TensorFunc
class Lshape : public STensorFunc
{
public:
//! \brief Constructor with some default parameters.
@ -57,13 +57,13 @@ public:
protected:
//! \brief Evaluates the analytic stress tensor at the point \a x.
virtual Tensor evaluate(const Vec3& x) const;
virtual SymmTensor evaluate(const Vec3& x) const;
private:
bool is3D; //!< Flag telling whether to return a 3D stress tensor or not
double a; //!< Length parameter
double F0; //!< Load factor
double nu; //!< Poisson's ratio
bool is3D; //!< Flag telling whether to return a 3D stress tensor or not
Tensor T; //!< Local-to-global stress transformation tensor
};
@ -73,7 +73,7 @@ private:
\brief Analytic solution for the cantilever beam with a tip shear load.
*/
class CanTS : public TensorFunc
class CanTS : public STensorFunc
{
public:
//! \brief Constructor with some default parameters.
@ -84,13 +84,13 @@ public:
protected:
//! \brief Evaluates the analytic stress tensor at the point \a x.
virtual Tensor evaluate(const Vec3& x) const;
virtual SymmTensor evaluate(const Vec3& x) const;
private:
bool is3D; //!< Flag telling whether to return a 3D stress tensor or not
double L; //!< Length
double H; //!< Height
double F0; //!< Load factor
bool is3D; //!< Flag telling whether to return a 3D stress tensor or not
};
@ -98,7 +98,7 @@ private:
\brief Analytic solution for the cantilever beam with a tip moment load.
*/
class CanTM : public TensorFunc
class CanTM : public STensorFunc
{
public:
//! \brief Constructor with some default parameters.
@ -109,12 +109,12 @@ public:
protected:
//! \brief Evaluates the analytic stress tensor at the point \a x.
virtual Tensor evaluate(const Vec3& x) const;
virtual SymmTensor evaluate(const Vec3& x) const;
private:
bool is3D; //!< Flag telling whether to return a 3D stress tensor or not
double H; //!< Height
double M0; //!< Load factor
bool is3D; //!< Flag telling whether to return a 3D stress tensor or not
};
@ -122,7 +122,7 @@ private:
\brief Analytic solution for the curved beam with end shear.
*/
class CurvedBeam : public TensorFunc
class CurvedBeam : public STensorFunc
{
public:
//! \brief Constructor with some default parameters.
@ -133,13 +133,13 @@ public:
protected:
//! \brief Evaluates the analytic stress tensor at the point \a x.
virtual Tensor evaluate(const Vec3& x) const;
virtual SymmTensor evaluate(const Vec3& x) const;
private:
bool is3D; //!< Flag telling whether to return a 3D stress tensor or not
double a; //!< Inner radius
double b; //!< Outer radius
double PN; //!< Load parameter
bool is3D; //!< Flag telling whether to return a 3D stress tensor or not
mutable Tensor T; //!< Local-to-global stress transformation tensor
};
@ -242,7 +242,7 @@ protected:
class PoissonLine : public VecFunc
{
public:
//! \brief Empty Constructor.
//! \brief Empty constructor.
PoissonLine(double r = 1.0) : L(r) {}
//! \brief Empty destructor.
virtual ~PoissonLine() {}

69
src/Integrands/AnalyticSolutionsStokes.C
View File

@ -1,3 +1,4 @@
// $Id$
//==============================================================================
//!
//! \file AnalyticSolutionsStokes.C
@ -12,39 +13,33 @@
#include "AnalyticSolutionsStokes.h"
#include "Vec3.h"
#include "Vec3Oper.h"
#include "Tensor.h"
#include <math.h>
// Define pi
#ifndef M_PI
#define PI 3.141592653589793
#else
#define PI M_PI
#endif
/*!
\class Poiseuille
Channel flow
Channel flow.
*/
Poiseuille::Poiseuille(double P, double diam, double len, double visc)
: Pin(P), D(diam), L(len), mu(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);
//vecGrad = new VelocityGrad(*this);
scalarSol = true;
vectorSol = true;
}
Poiseuille::~Poiseuille() {}
// Pressure solution
real Poiseuille::Pressure::evaluate(const Vec3& x) const
double Poiseuille::Pressure::evaluate(const Vec3& x) const
{
return params.Pin*(1.0-x[0]/params.L);
}
@ -53,15 +48,13 @@ real Poiseuille::Pressure::evaluate(const Vec3& x) const
// Velocity solution
Vec3 Poiseuille::Velocity::evaluate(const Vec3& x) const
{
Vec3 U;
U = 0.0;
double dPdx = -params.Pin/params.L;
double coeff = -dPdx/(2.0*params.mu);
U[0] = coeff*x[1]*(params.D-x[1]);
return U;
Vec3 u;
u[0] = coeff*x[1]*(params.D-x[1]);
return u;
}
@ -69,8 +62,6 @@ Vec3 Poiseuille::Velocity::evaluate(const Vec3& x) const
Vec3 Poiseuille::PressureGrad::evaluate(const Vec3& x) const
{
Vec3 dPdX;
dPdX = 0.0;
dPdX[0] = -params.Pin/params.L;
return dPdX;
@ -81,7 +72,6 @@ Vec3 Poiseuille::PressureGrad::evaluate(const Vec3& x) const
Tensor Poiseuille::VelocityGrad::evaluate(const Vec3& x) const
{
Tensor dUdX(3);
dUdX.zero();
dUdX(1,2) = params.Pin/(2.0*params.L*params.mu)*(1.0-2.0*x[1]);
return dUdX;
@ -91,31 +81,24 @@ Tensor Poiseuille::VelocityGrad::evaluate(const Vec3& x) const
/*!
\class TestSolution
Analytical test solution using trigonometrical functions
Analytical test solution using trigonometrical functions.
*/
TestSolution::TestSolution(double density, double visc)
: rho(density), mu(visc)
TestSolution::TestSolution(double dens, double visc)
: rho(dens), mu(visc)
{
scalSol = new Pressure();
vecSol = new Velocity();
scalSecSol = new PressureGrad();
vecSecSol = new VelocityGrad();
//vecGrad = new VelocityGrad(*this);
scalarSol = true;
vectorSol = true;
}
TestSolution::~TestSolution() {}
}
// Pressure
real TestSolution::Pressure::evaluate(const Vec3& x) const
double TestSolution::Pressure::evaluate(const Vec3& x) const
{
const Vec4* X = dynamic_cast<const Vec4*>(&x);
double t = X->t;
double t = X ? X->t : 0.0;
return cos(PI*x[0])*sin(PI*x[1])*sin(t);
}
@ -125,12 +108,11 @@ real TestSolution::Pressure::evaluate(const Vec3& x) const
Vec3 TestSolution::Velocity::evaluate(const Vec3& x) const
{
const Vec4* X = dynamic_cast<const Vec4*>(&x);
double t = X->t;
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);
u[2] = 0.0;
return u;
}
@ -140,12 +122,11 @@ Vec3 TestSolution::Velocity::evaluate(const Vec3& x) const
Vec3 TestSolution::PressureGrad::evaluate(const Vec3& x) const
{
const Vec4* X = dynamic_cast<const Vec4*>(&x);
double t = X->t;
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);
dPdX[2] = 0.0;
return dPdX;
}
@ -155,15 +136,13 @@ Vec3 TestSolution::PressureGrad::evaluate(const Vec3& x) const
Tensor TestSolution::VelocityGrad::evaluate(const Vec3& x) const
{
const Vec4* X = dynamic_cast<const Vec4*>(&x);
double t = X->t;
double t = X ? X->t : 0.0;
Tensor dUdX(3);
dUdX.zero();
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(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*PI*x[0])*sin(PI*x[1])*cos(PI*x[1])*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;
}

92
src/Integrands/AnalyticSolutionsStokes.h
View File

@ -1,3 +1,4 @@
// $Id$
//==============================================================================
//!
//! \file AnalyticSolutionsStokes.h
@ -16,140 +17,137 @@
#include "AnaSol.h"
#include "Tensor.h"
/*!
\brief Analytic solution for Poiseuille flow
\brief Analytic solution for Poiseuille flow.
*/
class Poiseuille : public AnaSol
{
public:
//! \brief Constructor for pressure inlet BC
Poiseuille(double P = 1.0, double diam = 1.0, double len = 1.0, double visc = 1.0);
//! \brief Destructor
virtual ~Poiseuille();
//! \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 value
double D; //!< Diameter of channel
double L; //!< Length of channel
double mu; //!< Fluid viscosity
double Pin; //!< Inlet pressure value
// Analytical pressure
//! \brief Analytical pressure field.
class Pressure : public RealFunc
{
private:
Poiseuille& params;
const Poiseuille& params;
protected:
real evaluate(const Vec3& x) const;
virtual double evaluate(const Vec3& x) const;
public:
Pressure(Poiseuille& p) : params(p) {}
Pressure(const Poiseuille& p) : params(p) {}
virtual ~Pressure() {}
};
// Analytical velocity
//! \brief Analytical velocity field.
class Velocity : public VecFunc
{
private:
Poiseuille& params;
const Poiseuille& params;
protected:
Vec3 evaluate(const Vec3& x) const;
virtual Vec3 evaluate(const Vec3& x) const;
public:
Velocity(Poiseuille& p) : params(p) {}
Velocity(const Poiseuille& p) : params(p) {}
virtual ~Velocity() {}
};
// Analytical pressure gradient
//! \brief Analytical pressure gradient field.
class PressureGrad : public VecFunc
{
private:
Poiseuille& params;
const Poiseuille& params;
protected:
Vec3 evaluate(const Vec3& x) const;
virtual Vec3 evaluate(const Vec3& x) const;
public:
PressureGrad(Poiseuille& p) : params(p) {}
PressureGrad(const Poiseuille& p) : params(p) {}
virtual ~PressureGrad() {}
};
// Analytical velocity gradient
//! \brief Analytical velocity gradient field.
class VelocityGrad : public TensorFunc
{
private:
Poiseuille& params;
const Poiseuille& params;
protected:
Tensor evaluate(const Vec3& x) const;
virtual Tensor evaluate(const Vec3& x) const;
public:
VelocityGrad(Poiseuille& p) : params(p) {}
VelocityGrad(const Poiseuille& p) : params(p) {}
virtual ~VelocityGrad() {}
};
};
/*!
\brief Artificial solution for Navier-Stokes equations
\brief Artificial solution for Navier-Stokes equations.
*/
class TestSolution : public AnaSol
{
public:
//! \brief Constructor
//! \param[in] rho Fluid density
//! \param[in] mu Fluid viscosity
TestSolution(double rho = 1.0, double mu = 1.0);
//! \brief Destructor
virtual ~TestSolution();
//! \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
double rho; //!< Fluid density
double mu; //!< Fluid viscosity
// Analytical pressure
//! \brief Analytical pressure field.
class Pressure : public RealFunc
{
protected:
real evaluate(const Vec3& x) const;
virtual double evaluate(const Vec3& x) const;
public:
Pressure() {}
virtual ~Pressure() {}
};
// Analytical velocity
//! \brief Analytical velocity field.
class Velocity : public VecFunc
{
protected:
Vec3 evaluate(const Vec3& x) const;
virtual Vec3 evaluate(const Vec3& x) const;
public:
Velocity() {}
virtual ~Velocity() {}
};
// Analytical pressure gradient
//! \brief Analytical pressure gradient field.
class PressureGrad : public VecFunc
{
protected:
Vec3 evaluate(const Vec3& x) const;
virtual Vec3 evaluate(const Vec3& x) const;
public:
PressureGrad() {}
virtual ~PressureGrad() {}
};
// Analytical velocity gradient
//! \brief Analytical velocity gradient field.
class VelocityGrad : public TensorFunc
{
protected:
Tensor evaluate(const Vec3& x) const;
virtual Tensor evaluate(const Vec3& x) const;
public:
VelocityGrad() {}

32
src/Integrands/Elasticity.C
View File

@ -1,4 +1,4 @@
// $Id: Elasticity.C,v 1.1 2011-02-08 09:06:02 kmo Exp $
// $Id$
//==============================================================================
//!
//! \file Elasticity.C
@ -17,19 +17,20 @@
#include "ElmNorm.h"
#include "Tensor.h"
#include "Vec3Oper.h"
#include "AnaSol.h"
#include "VTF.h"
Elasticity::Elasticity (unsigned short int n, bool ps)
Elasticity::Elasticity (unsigned short int n, bool ps) : nsd(n), planeStress(ps)
{
nsd = n;
planeStress = ps;
// Default material properties - typical values for steel (SI units)
Emod = 2.05e11;
nu = 0.29;
rho = 7.85e3;
// Default is zero gravity
g[0] = g[1] = g[2] = 0.0;
myMats = new ElmMats();
locSys = 0;
@ -161,7 +162,9 @@ bool Elasticity::initElement (const std::vector<int>& MNPC)
std::cerr <<" *** Elasticity::initElement: Detected "
<< ierr <<" node numbers out of range."<< std::endl;
myMats->withLHS = true;
if (myMats)
myMats->withLHS = true;
return ierr == 0;
}
@ -178,7 +181,9 @@ bool Elasticity::initElementBou (const std::vector<int>& MNPC)
std::cerr <<" *** Elasticity::initElement: Detected "
<< ierr <<" node numbers out of range."<< std::endl;
myMats->withLHS = false;
if (myMats)
myMats->withLHS = false;
return ierr == 0;
}
@ -585,9 +590,10 @@ bool Elasticity::evalSol (Vector& s, const Matrix& dNdX, const Vec3& X) const
}
bool Elasticity::evalSol (Vector& s, const TensorFunc& sol, const Vec3& X) const
bool Elasticity::evalSol (Vector& s, const STensorFunc& asol,
const Vec3& X) const
{
s = sol(X);
s = asol(X);
s.push_back(SymmTensor(s).vonMises());
return true;
}
@ -623,9 +629,13 @@ const char* Elasticity::getFieldLabel (size_t i, const char* prefix) const
}
NormBase* Elasticity::getNormIntegrand (TensorFunc* sol) const
NormBase* Elasticity::getNormIntegrand (AnaSol* asol) const
{
return new ElasticityNorm(*const_cast<Elasticity*>(this),sol);
if (asol)
return new ElasticityNorm(*const_cast<Elasticity*>(this),
asol->getStressSol());
else
return new ElasticityNorm(*const_cast<Elasticity*>(this));
}

17
src/Integrands/Elasticity.h
View File

@ -1,4 +1,4 @@
// $Id: Elasticity.h,v 1.1 2011-02-08 09:06:02 kmo Exp $
// $Id$
//==============================================================================
//!
//! \file Elasticity.h
@ -21,7 +21,6 @@
class LocalSystem;
class ElmMats;
class ElmNorm;
class Tensor;
class VTF;
@ -57,7 +56,7 @@ public:
void clearTracVal() { tracVal.clear(); }
//! \brief Defines the gravitation vector.
void setGravity(double gx, double gy, double gz)
void setGravity(double gx, double gy = 0.0, double gz = 0.0)
{ g[0] = gx; g[1] = gy; g[2] = gz; }
//! \brief Defines the body force field.
@ -98,9 +97,9 @@ public:
//! \brief Evaluates the analytical solution at an integration point.
//! \param[out] s The analytical stress values at current point
//! \param[in] sol The analytical solution field
//! \param[in] asol The analytical solution field
//! \param[in] X Cartesian coordinates of current point
virtual bool evalSol(Vector& s, const TensorFunc& sol, const Vec3& X) const;
virtual bool evalSol(Vector& s, const STensorFunc& asol, const Vec3& X) const;
//! \brief Evaluates the primary solution at a result point.
//! \param[in] N Basis function values at current point
@ -127,8 +126,8 @@ public:
//! \note The Integrand object is allocated dynamically and has to be deleted
//! manually when leaving the scope of the pointer variable receiving the
//! returned pointer value.
//! \param[in] sol Pointer to analytical solution field (optional)
virtual NormBase* getNormIntegrand(TensorFunc* sol = 0) const;
//! \param[in] asol Pointer to analytical solution fields (optional)
virtual NormBase* getNormIntegrand(AnaSol* asol = 0) const;
//! \brief Returns the number of secondary solution fields.
virtual size_t getNoFields() const;
@ -246,7 +245,7 @@ public:
//! \brief The only constructor initializes its data members.
//! \param[in] p The linear elasticity problem to evaluate norms for
//! \param[in] a The analytical stress field (optional)
ElasticityNorm(Elasticity& p, TensorFunc* a) : problem(p), anasol(a) {}
ElasticityNorm(Elasticity& p, STensorFunc* a = 0) : problem(p), anasol(a) {}
//! \brief Empty destructor.
virtual ~ElasticityNorm() {}
@ -297,7 +296,7 @@ protected:
Elasticity& problem; //!< The problem-specific data
private:
TensorFunc* anasol; //!< Analytical stress field
STensorFunc* anasol; //!< Analytical stress field
};
#endif

23
src/Integrands/Poisson.C
View File

@ -1,4 +1,4 @@
// $Id: Poisson.C,v 1.6 2010-12-27 18:29:01 kmo Exp $
// $Id$
//==============================================================================
//!
//! \file Poisson.C
@ -16,6 +16,7 @@
#include "ElmNorm.h"
#include "Tensor.h"
#include "Vec3Oper.h"
#include "AnaSol.h"
#include "VTF.h"
@ -86,7 +87,7 @@ bool Poisson::initElement (const std::vector<int>& MNPC)
int ierr = 0;
if (eV && !primsol.front().empty())
if (ierr = utl::gather(MNPC,1,primsol.front(),*eV))
if ((ierr = utl::gather(MNPC,1,primsol.front(),*eV)))
std::cerr <<" *** Poisson::initElement: Detected "
<< ierr <<" node numbers out of range."<< std::endl;
@ -199,8 +200,8 @@ bool Poisson::evalSol (Vector& q, const Vector&,
return false;
Vector Dtmp;
int ierr = 0;
if (ierr = utl::gather(MNPC,1,primsol.front(),Dtmp))
int ierr = utl::gather(MNPC,1,primsol.front(),Dtmp);
if (ierr > 0)
{
std::cerr <<" *** Poisson::evalSol: Detected "
<< ierr <<" node numbers out of range."<< std::endl;
@ -215,8 +216,7 @@ bool Poisson::evalSol (Vector& q, const Vector&,
}
bool Poisson::evalSol (Vector& q,
const Matrix& dNdX, const Vec3& X) const
bool Poisson::evalSol (Vector& q, const Matrix& dNdX, const Vec3& X) const
{
if (!eV || eV->empty())
{
@ -242,8 +242,7 @@ bool Poisson::evalSol (Vector& q,
}
bool Poisson::evalSolScal (Vector& s,
const VecFunc& asol, const Vec3& X) const
bool Poisson::evalSol (Vector& s, const VecFunc& asol, const Vec3& X) const
{
s = Vector(asol(X).ptr(),nsd);
return true;
@ -266,9 +265,13 @@ const char* Poisson::getFieldLabel (size_t i, const char* prefix) const
}
NormBase* Poisson::getNormIntegrandScal (VecFunc* asol) const
NormBase* Poisson::getNormIntegrand (AnaSol* asol) const
{
return new PoissonNorm(*const_cast<Poisson*>(this),asol);
if (asol)
return new PoissonNorm(*const_cast<Poisson*>(this),
asol->getScalarSecSol());
else
return new PoissonNorm(*const_cast<Poisson*>(this));
}

12
src/Integrands/Poisson.h
View File

@ -1,4 +1,4 @@
// $Id: Poisson.h,v 1.4 2010-12-27 18:29:02 kmo Exp $
// $Id$
//==============================================================================
//!
//! \file Poisson.h
@ -98,9 +98,9 @@ public:
//! \brief Evaluates the analytical solution at an integration point.
//! \param[out] s The analytical solution values at current point
//! \param[in] asol The analytical solution field
//! \param[in] asol The analytical solution field (heat flux)
//! \param[in] X Cartesian coordinates of current point
virtual bool evalSolScal(Vector& s, const VecFunc& asol, const Vec3& X) const;
virtual bool evalSol(Vector& s, const VecFunc& 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
@ -115,8 +115,8 @@ public:
//! \note The Integrand object is allocated dynamically and has to be deleted
//! manually when leaving the scope of the pointer variable receiving the
//! returned pointer value.
//! \param[in] asol Pointer to analytical solution field (optional)
virtual NormBase* getNormIntegrandScal(VecFunc* asol = 0) const;
//! \param[in] asol Pointer to analytical solution (optional)
virtual NormBase* getNormIntegrand(AnaSol* asol = 0) const;
//! \brief Returns the number of secondary solution fields.
virtual size_t getNoFields() const { return nsd; }
@ -173,7 +173,7 @@ public:
//! \brief The only constructor initializes its data members.
//! \param[in] p The Poisson problem to evaluate norms for
//! \param[in] a The analytical heat flux (optional)
PoissonNorm(Poisson& p, VecFunc* a) : problem(p), anasol(a) {}
PoissonNorm(Poisson& p, VecFunc* a = 0) : problem(p), anasol(a) {}
//! \brief Empty destructor.
virtual ~PoissonNorm() {}

74
src/SIM/SIMbase.C
View File

@ -1,4 +1,4 @@
// $Id: SIMbase.C,v 1.50 2011-02-08 12:52:12 rho Exp $
// $Id$
//==============================================================================
//!
//! \file SIMbase.C
@ -25,6 +25,7 @@
#include "LinSolParams.h"
#include "GlbNorm.h"
#include "ElmNorm.h"
#include "AnaSol.h"
#include "Tensor.h"
#include "Vec3.h"
#include "Vec3Oper.h"
@ -45,9 +46,10 @@ int SIMbase::num_threads_SLU = 1;
SIMbase::SIMbase ()
{
myProblem = 0;
mySol = 0;
myVtf = 0;
mySam = 0;
myEqSys = 0;
mySam = 0;
mySolParams = 0;
nGlPatches = 0;
@ -68,6 +70,7 @@ SIMbase::~SIMbase ()
#endif
if (myProblem) delete myProblem;
if (mySol) delete mySol;
if (myVtf) delete myVtf;
if (myEqSys) delete myEqSys;
if (mySam) delete mySam;
@ -424,7 +427,7 @@ bool SIMbase::assembleSystem (const TimeDomain& time, const Vectors& prevSol,
else
ok = false;
if (j == 0 && ok)
if (j == 0 && ok)
// All patches are referring to the same material, and we assume it has
// been initialized during input processing (thus no initMaterial call here)
for (i = 0; i < myModel.size() && ok; i++)
@ -540,19 +543,24 @@ bool SIMbase::solveSystem (Vector& solution, int printSol,
if (printSol < 1) return true;
// Compute and print solution norms
const int nsd = this->getNoSpaceDim();
size_t iMax[nsd];
double dMax[nsd];
double dNorm = this->solutionNorms(solution,dMax,iMax);
const size_t nf = myModel.front()->getNoFields(1);
size_t iMax[nf];
double dMax[nf];
double dNorm = this->solutionNorms(solution,dMax,iMax,nf);
if (myPid == 0)
{
std::cout <<"\n >>> Solution summary <<<\n"
<<"\nL2-norm : "<< dNorm;
char D = 'X';
for (int d = 0; d < nsd; d++, D++)
std::cout <<"\nMax "<< D <<'-'<< compName <<" : "
<< dMax[d] <<" node "<< iMax[d];
if (nf == 1)
std::cout <<"\nMax "<< compName <<" : "<< dMax[0] <<" node "<< iMax[0];
else
{
char D = 'X';
for (size_t d = 0; d < nf; d++, D++)
std::cout <<"\nMax "<< D <<'-'<< compName <<" : "
<< dMax[d] <<" node "<< iMax[d];
}
std::cout << std::endl;
}
@ -583,9 +591,12 @@ void SIMbase::iterationNorms (const Vector& u, const Vector& r,
}
double SIMbase::solutionNorms (const Vector& x, double* inf, size_t* ind) const
double SIMbase::solutionNorms (const Vector& x, double* inf,
size_t* ind, size_t nf) const
{
for (size_t d = 0; d < this->getNoSpaceDim(); d++)
if (nf == 0) nf = this->getNoSpaceDim();
for (size_t d = 0; d < nf; d++)
{
ind[d] = d+1;
inf[d] = mySam->normInf(x,ind[d]);
@ -598,7 +609,7 @@ double SIMbase::solutionNorms (const Vector& x, double* inf, size_t* ind) const
bool SIMbase::solutionNorms (const TimeDomain& time, const Vectors& psol,
Matrix& eNorm, Vector& gNorm)
{
NormBase* norm = myProblem->getNormIntegrand(this->getAnaSol());
NormBase* norm = myProblem->getNormIntegrand(mySol);
if (!norm)
{
std::cerr <<" *** SIMbase::solutionNorms: No integrand."<< std::endl;
@ -814,9 +825,9 @@ bool SIMbase::writeGlvBC (const int* nViz, int& nBlock) const
if (myModel[i]->empty()) continue; // skip empty patches
geomID++;
size_t nbc = myModel.front()->getNoFields(1);
RealArray flag(3,0.0);
size_t nbc = myModel[i]->getNoFields(1);
Matrix bc(nbc,myModel[i]->getNoNodes());
RealArray flag(3,0.0);
ASMbase::BCVec::const_iterator bit;
for (bit = myModel[i]->begin_BC(); bit != myModel[i]->end_BC(); bit++)
if ((node = myModel[i]->getNodeIndex(bit->node)))
@ -921,14 +932,12 @@ bool SIMbase::writeGlvS (const Vector& psol,
int geomID = 0, idBlock = 10;
std::vector<int> vID, dID[14], sID[14];
bool haveAsol = false;
bool scalarEq = myModel.empty() ? false : myModel.front()->getNoFields() == 1;
if (scalarEq) {
if (getAnaSol() && this->getAnaSol()->hasScalarSol())
haveAsol = true;
}
else
if (this->getAnaSol())
haveAsol = true;
bool scalarEq = myModel.front()->getNoFields() == 1;
if (mySol)
if (scalarEq)
haveAsol = mySol->hasScalarSol() > 1;
else
haveAsol = mySol->hasVectorSol() > 1;
for (i = 0; i < myModel.size(); i++)
{
@ -991,9 +1000,11 @@ bool SIMbase::writeGlvS (const Vector& psol,
for (j = 1; cit != grid->end_XYZ() && ok; j++, cit++)
{
if (scalarEq)
ok = myProblem->evalSecSolScal(solPt,*this->getAnaSol()->getScalarSecSol(),*cit);
ok = myProblem->evalSol(solPt,*mySol->getScalarSecSol(),*cit);
else if (mySol->hasVectorSol() == 3)
ok = myProblem->evalSol(solPt,*mySol->getStressSol(),*cit);
else
ok = myProblem->evalSecSol(solPt,*this->getAnaSol()->getVectorSecSol(),*cit);
ok = myProblem->evalSol(solPt,*mySol->getVectorSecSol(),*cit);
if (ok)
field.fillColumn(j,solPt);
}
@ -1201,7 +1212,6 @@ bool SIMbase::dumpSolution (const Vector& psol, std::ostream& os) const
Matrix field;
size_t i, j, k;
const int nsd = this->getNoSpaceDim();
for (i = 0; i < myModel.size(); i++)
{
@ -1211,15 +1221,19 @@ bool SIMbase::dumpSolution (const Vector& psol, std::ostream& os) const
os <<"\n# Patch: "<< i+1;
// Extract and write primary solution
size_t nf = myModel[i]->getNoFields(1);
Vector& patchSol = myProblem->getSolution();
myModel[i]->extractNodeVec(psol,patchSol);
for (int d = 1; d <= nsd; d++)
for (k = 1; k <= nf; k++)
{
os <<"# FE u_"<< char('w'+d);
os <<"# FE u";
if (nf > 1)
os <<'_'<< char('w'+k);
for (j = 1; j <= myModel[i]->getNoNodes(); j++)
{
std::pair<int,int> dofs = mySam->getNodeDOFs(j);
int idof = dofs.first+d-1;
int idof = dofs.first+k-1;
if (idof <= dofs.second)
os <<"\n"<< patchSol[idof-1];
}

14
src/SIM/SIMbase.h
View File

@ -1,4 +1,4 @@
// $Id: SIMbase.h,v 1.35 2011-02-08 15:51:16 rho Exp $
// $Id$
//==============================================================================
//!
//! \file SIMbase.h
@ -19,11 +19,11 @@
#include "TimeDomain.h"
#include "Property.h"
#include "Function.h"
#include "AnaSol.h"
#include <map>
class ASMbase;
class Integrand;
class AnaSol;
class VTF;
class SAMpatch;
class AlgEqSystem;
@ -98,7 +98,7 @@ public:
//! \brief Prints out problem-specific data to the given stream.
void printProblem(std::ostream& os) const;
//! \brief Returns the number of spatial dimensions in the model.
size_t getNoSpaceDim() const;
//! \brief Returns the model size in terms of number of DOFs.
@ -169,8 +169,10 @@ public:
//! \param[in] x Global primary solution vector
//! \param[out] inf Infinity norms in each spatial direction
//! \param[out] ind Global index of the node corresponding to the inf-value
//! \param[in] nf Number of components in the primary solution field
//! \return L2-norm of the solution vector
virtual double solutionNorms(const Vector& x, double* inf, size_t* ind) const;
virtual double solutionNorms(const Vector& x, double* inf,
size_t* ind, size_t nf = 0) const;
//! \brief Integrates some solution norm quantities.
//! \details If an analytical solution is provided, norms of the exact
@ -322,9 +324,6 @@ protected:
//! \brief Initializes for integration of Neumann terms for a given property.
virtual bool initNeumann(size_t) { return true; }
//! \brief Returns the analytical solution field, if any.
virtual AnaSol* getAnaSol() const { return 0; }
//! \brief Extract local solution vector(s) for a given patch.
//! \param[in] patch Pointer to the patch to extract solution vectors for
//! \param[in] sol Global primary solution vectors in DOF-order
@ -357,6 +356,7 @@ protected:
VecFuncMap myVectors; //!< Vector property fields
TracFuncMap myTracs; //!< Traction property fields
Integrand* myProblem; //!< Problem-specific data and methods
AnaSol* mySol; //!< Analytical/Exact solution
VTF* myVtf; //!< VTF-file for result visualization
// Parallel computing attributes

39
src/Utility/AnaSol.C
View File

@ -1,39 +0,0 @@
//==============================================================================
//!
//! \file AnaSol.h
//!
//! \date Des 7 2010
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Analytical solution fields (primary and secondary)
//!
//==============================================================================
#include "AnaSol.h"
AnaSol::AnaSol()
{
scalarSol = vectorSol = false;
}
AnaSol::AnaSol(RealFunc* sp, VecFunc* ss, VecFunc* vp, TensorFunc* vs)
: scalSol(sp), scalSecSol(ss), vecSol(vp), vecSecSol(vs)
{
scalarSol = vectorSol = false;
if (scalSol || scalSecSol) scalarSol = true;
if (vecSol || vecSecSol) vectorSol = true;
}
AnaSol::~AnaSol()
{
if (scalSol) delete scalSol;
if (scalSecSol) delete scalSecSol;
if (vecSol) delete vecSol;
if (vecSecSol) delete vecSecSol;
scalarSol = vectorSol = false;
}

72
src/Utility/AnaSol.h
View File

@ -1,72 +0,0 @@
//==============================================================================
//!
//! \file AnaSol.h
//!
//! \date Des 7 2010
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Analytical solution fields (primary and secondary)
//!
//==============================================================================
#ifndef ANA_SOL_H
#define ANA_SOL_H
#include "Function.h"
/*!
\brief Data type for analytical scalar solution fields (primary and secondary)
*/
class AnaSol
{
protected:
bool scalarSol; //!< If scalar solution field defined
bool vectorSol; //!< If a vector solution field is defined
// Scalar analytical solutions
RealFunc* scalSol; //<! Primary scalar solution field
VecFunc* scalSecSol; //<! Secondary scalar solution fields
// Vector solutions
VecFunc* vecSol; //<! Primary Vec3 solution fields
TensorFunc* vecSecSol; //<! Secondary vector solution fields
//GradientFunc* vecGrad; //<! Gradient of solution fields
public:
//! \brief Constructor
AnaSol();
//! \brief Constructor
//! \param[in] sp Primary scalar solution field
//! \param[in] ss Secondary scalar solution field
//! \param[in] vp Primary vector solution field
//! \param[in] vs Secondary vector solution field
AnaSol(RealFunc* sp, VecFunc* ss, VecFunc* vp, TensorFunc* vs);
//! \brief Destructor
virtual ~AnaSol();
//! \brief Returns \e true if a scalar solution is defined
virtual bool hasScalarSol() const { return scalarSol; }
//! \brief Returns \e true if a vector solution is defined
virtual bool hasVectorSol() const { return vectorSol; }
//! \brief Returns the scalar solution if any
virtual RealFunc* getScalarSol() const { return scalSol; }
//! \brief Returns the secondary scalar solution if any
virtual VecFunc* getScalarSecSol() const { return scalSecSol; }
//! \brief Returns the vector solution if any
virtual VecFunc* getVectorSol() const { return vecSol; }
//! \brief Returns the secondary vector solution if any
virtual TensorFunc* getVectorSecSol() const { return vecSecSol; }
};
#endif

9
src/Utility/Function.C
View File

@ -1,4 +1,4 @@
// $Id: Function.C,v 1.6 2010-10-14 19:10:55 kmo Exp $
// $Id$
//==============================================================================
//!
//! \file Function.C
@ -42,5 +42,10 @@ Vec3 PressureField::evaluate (const Vec3& x, const Vec3& n) const
Vec3 TractionField::evaluate (const Vec3& x, const Vec3& n) const
{
return stress(x) * n;
if (sigma) // symmetric tensor field
return (*sigma)(x) * n;
else if (sigmaN) // non-symmetric tensor field
return (*sigmaN)(x) * n;
else
return Vec3();
}

16
src/Utility/Function.h
View File

@ -1,4 +1,4 @@
// $Id: Function.h,v 1.7 2011-02-08 15:07:19 rho Exp $
// $Id$
//==============================================================================
//!
//! \file Function.h
@ -84,9 +84,12 @@ typedef utl::Function<Vec3,real> RealFunc;
//! \brief Vector-valued unary function of a spatial point.
typedef utl::Function<Vec3,Vec3> VecFunc;
//! \brief Symmetric tensor-valued unary function of a spatial point.
//! \brief Tensor-valued unary function of a spatial point.
typedef utl::Function<Vec3,Tensor> TensorFunc;
//! \brief Symmetric tensor-valued unary function of a spatial point.
typedef utl::Function<Vec3,SymmTensor> STensorFunc;
/*!
\brief Vector-valued binary function of a spatial point and normal vector.
@ -136,11 +139,14 @@ protected:
class TractionField : public TractionFunc
{
const TensorFunc& stress; //!< The tensor field to derive the traction from
const STensorFunc* sigma; //!< Symmetric tensor field to derive tractions from
const TensorFunc* sigmaN; //!< Tensor field to derive tractions from
public:
//! \brief Constructor initializing the tensor function reference.
TractionField(const TensorFunc& field) : stress(field) {}
//! \brief Constructor initializing the symmetric tensor function pointer.
TractionField(const STensorFunc& field) : sigma(&field), sigmaN(0) {}
//! \brief Constructor initializing the tensor function pointer.
TractionField(const TensorFunc& field) : sigma(0), sigmaN(&field) {}
//! \brief Empty destructor.
virtual ~TractionField() {}

60
src/Utility/Functions.C
View File

@ -1,4 +1,4 @@
// $Id: Functions.C,v 1.8 2011-02-08 12:55:52 rho Exp $
// $Id$
//==============================================================================
//!
//! \file Functions.C
@ -118,33 +118,34 @@ const RealFunc* utl::parseRealFunc (char* cline, real A)
// Check for spatial variation
int linear = 0;
int quadratic = 0;
if (cline)
if (strcmp(cline,"X") == 0)
linear = 1;
else if (strcmp(cline,"Y") == 0)
linear = 2;
else if (strcmp(cline,"Z") == 0)
linear = 3;
else if (strcmp(cline,"XrotZ") == 0)
linear = 4;
else if (strcmp(cline,"YrotZ") == 0)
linear = 5;
else if (strcmp(cline,"Tinit") == 0)
linear = 6;
else if (strcmp(cline,"quadX") == 0)
quadratic = 1;
else if (strcmp(cline,"quadY") == 0)
quadratic = 2;
else if (strcmp(cline,"quadZ") == 0)
quadratic = 1;
else if (strcmp(cline,"StepX") == 0)
linear = 7;
else if (strcmp(cline,"StepXY") == 0)
linear = 8;
if (!cline)
linear = -1;
else if (strcmp(cline,"X") == 0)
linear = 1;
else if (strcmp(cline,"Y") == 0)
linear = 2;
else if (strcmp(cline,"Z") == 0)
linear = 3;
else if (strcmp(cline,"XrotZ") == 0)
linear = 4;
else if (strcmp(cline,"YrotZ") == 0)
linear = 5;
else if (strcmp(cline,"Tinit") == 0)
linear = 6;
else if (strcmp(cline,"quadX") == 0)
quadratic = 1;
else if (strcmp(cline,"quadY") == 0)
quadratic = 2;
else if (strcmp(cline,"quadZ") == 0)
quadratic = 3;
else if (strcmp(cline,"StepX") == 0)
linear = 7;
else if (strcmp(cline,"StepXY") == 0)
linear = 8;
real C = A;
const RealFunc* f = 0;
if (linear && (cline = strtok(NULL," ")))
if (linear > 0 && (cline = strtok(NULL," ")))
{
C = real(1);
std::cout <<"("<< A <<"*";
@ -200,18 +201,19 @@ const RealFunc* utl::parseRealFunc (char* cline, real A)
break;
case 3:
f = new QuadraticZFunc(A,a,b);
break;
break;
}
cline = strtok(NULL," ");
}
//std::cout << C;
}
if (quadratic == 0)
std::cout << C;
// Check for time variation
if (!cline) return f;
const ScalarFunc* s = 0;
double freq = atof(cline);
if (cline = strtok(NULL," "))
if ((cline = strtok(NULL," ")))
{
std::cout <<" * sin("<< freq <<"*t + "<< cline <<")";
s = new SineFunc(C,freq,atof(cline));