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

changed: new class for analytic solutions (credit runar)

Browse Source 
git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@786 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
akva 2011-02-09 10:17:47 +00:00 committed by Knut Morten Okstad
parent c88a28533d
commit 940d949764
22 changed files with 356 additions and 108 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
Apps/SIMLinEl2D.C
View File

@ -1,4 +1,4 @@
// $Id: SIMLinEl2D.C,v 1.18 2011-02-08 09:06:02 kmo Exp $
// $Id: SIMLinEl2D.C,v 1.19 2011-02-09 10:07:36 rho Exp $
//==============================================================================
//!
//! \file SIMLinEl2D.C
@ -118,7 +118,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 Hole(a,F0,nu);
asol = new AnaSol(NULL,NULL,NULL,new Hole(a,F0,nu));
std::cout <<"\nAnalytical solution: Hole a="<< a <<" F0="<< F0
<<" nu="<< nu << std::endl;
}
@ -127,7 +127,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 Lshape(a,F0,nu);
asol = new AnaSol(NULL,NULL,NULL,new Lshape(a,F0,nu));
std::cout <<"\nAnalytical solution: Lshape a="<< a <<" F0="<< F0
<<" nu="<< nu << std::endl;
}
@ -136,7 +136,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 CanTS(L,H,F0);
asol = new AnaSol(NULL,NULL,NULL,new CanTS(L,H,F0));
std::cout <<"\nAnalytical solution: CanTS L="<< L <<" H="<< H
<<" F0="<< F0 << std::endl;
}
@ -144,7 +144,7 @@ bool SIMLinEl2D::parse (char* keyWord, std::istream& is)
{
double H = atof(strtok(NULL," "));
double M0 = atof(strtok(NULL," "));
asol = new CanTM(H,M0);
asol = new AnaSol(NULL,NULL,NULL,new CanTM(H,M0));
std::cout <<"\nAnalytical solution: CanTM H="<< H
<<" M0="<< M0 << std::endl;
}
@ -154,7 +154,7 @@ 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 CurvedBeam(u0,a,b,E);
asol = new AnaSol(NULL,NULL,NULL,new CurvedBeam(u0,a,b,E));
std::cout <<"\nAnalytical solution: Curved Beam a="<< a <<" b="<< b
<<" u0="<< u0 <<" E="<< E << std::endl;
}
@ -164,11 +164,11 @@ bool SIMLinEl2D::parse (char* keyWord, std::istream& is)
// Define the analytical boundary traction field
int code = (cline = strtok(NULL," ")) ? atoi(cline) : 0;
if (code > 0 && asol)
if (code > 0 && asol->getVectorSecSol())
{
std::cout <<"Pressure code "<< code <<": Analytical traction"<< std::endl;
this->setPropertyType(code,Property::NEUMANN);
myTracs[code] = new TractionField(*asol);
myTracs[code] = new TractionField(*(asol->getVectorSecSol()));
}
}
@ -201,11 +201,11 @@ bool SIMLinEl2D::parse (char* keyWord, std::istream& is)
return false;
}
if (asol)
if (asol->getVectorSecSol())
{
std::cout <<"\tTraction on P"<< press.patch
<<" E"<< (int)press.lindx << std::endl;
myTracs[1+i] = new TractionField(*asol);
myTracs[1+i] = new TractionField(*asol->getVectorSecSol());
}
else
{

7
Apps/SIMLinEl2D.h
View File

@ -1,4 +1,4 @@
// $Id: SIMLinEl2D.h,v 1.9 2010-12-18 16:23:52 kmo Exp $
// $Id: SIMLinEl2D.h,v 1.10 2011-02-09 10:07:31 rho Exp $
//==============================================================================
//!
//! \file SIMLinEl2D.h
@ -16,6 +16,7 @@
#include "SIM2D.h"
#include "SIMenums.h"
#include "AnaSol.h"
/*!
@ -64,11 +65,11 @@ protected:
virtual bool initNeumann(size_t propInd);
//! \brief Returns the analytical solution function, if any.
virtual TensorFunc* getAnaSol() const { return asol; }
virtual AnaSol* getAnaSol() const { return asol; }
private:
MaterialVec mVec; //!< Material data
TensorFunc* asol; //!< Analytical stress field
AnaSol* asol; //!< Analytical stress field
};
#endif

18
Apps/SIMLinEl3D.C
View File

@ -1,4 +1,4 @@
// $Id: SIMLinEl3D.C,v 1.28 2011-02-08 09:06:02 kmo Exp $
// $Id: SIMLinEl3D.C,v 1.29 2011-02-09 10:08:02 rho Exp $
//==============================================================================
//!
//! \file SIMLinEl3D.C
@ -128,7 +128,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';
@ -232,7 +232,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 Hole(a,F0,nu,true);
asol = new AnaSol(NULL,NULL,NULL,new Hole(a,F0,nu,true));
std::cout <<"\nAnalytical solution: Hole a="<< a <<" F0="<< F0
<<" nu="<< nu << std::endl;
}
@ -241,7 +241,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 Lshape(a,F0,nu,true);
asol = new AnaSol(NULL,NULL,NULL,new Lshape(a,F0,nu,true));
std::cout <<"\nAnalytical solution: Lshape a="<< a <<" F0="<< F0
<<" nu="<< nu << std::endl;
}
@ -250,7 +250,7 @@ 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 CanTS(L,H,F0,true);
asol = new AnaSol(NULL,NULL,NULL,new CanTS(L,H,F0,true));
std::cout <<"\nAnalytical solution: CanTS L="<< L <<" H="<< H
<<" F0="<< F0 << std::endl;
}
@ -260,11 +260,11 @@ bool SIMLinEl3D::parse (char* keyWord, std::istream& is)
// Define the analytical boundary traction field
int code = (cline = strtok(NULL," ")) ? atoi(cline) : 0;
if (code > 0 && asol)
if (code > 0 && asol->getVectorSecSol())
{
std::cout <<"Pressure code "<< code <<": Analytical traction"<< std::endl;
this->setPropertyType(code,Property::NEUMANN);
myTracs[code] = new TractionField(*asol);
myTracs[code] = new TractionField(*(asol->getVectorSecSol()));
}
}
@ -297,11 +297,11 @@ bool SIMLinEl3D::parse (char* keyWord, std::istream& is)
return false;
}
if (asol)
if (asol && asol->getVectorSecSol())
{
std::cout <<"\tTraction on P"<< press.patch
<<" F"<< (int)press.lindx << std::endl;
myTracs[1+i] = new TractionField(*asol);
myTracs[1+i] = new TractionField(*asol->getVectorSecSol());
}
else
{

7
Apps/SIMLinEl3D.h
View File

@ -1,4 +1,4 @@
// $Id: SIMLinEl3D.h,v 1.11 2010-12-18 16:23:52 kmo Exp $
// $Id: SIMLinEl3D.h,v 1.12 2011-02-09 10:07:57 rho Exp $
//==============================================================================
//!
//! \file SIMLinEl3D.h
@ -16,6 +16,7 @@
#include "SIM3D.h"
#include "SIMenums.h"
#include "AnaSol.h"
/*!
@ -64,11 +65,11 @@ protected:
virtual bool initNeumann(size_t propInd);
//! \brief Returns the analytical solution function, if any.
virtual TensorFunc* getAnaSol() const { return asol; }
virtual AnaSol* getAnaSol() const { return asol; }
private:
MaterialVec mVec; //!< Material data
TensorFunc* asol; //!< Analytical stress field
AnaSol* asol; //!< Analytical solution fields
};
#endif

2
src/Integrands/AnalyticSolutionsStokes.h
View File

@ -14,7 +14,7 @@
#define _ANALYTIC_SOLUTIONS_STOKES_H
#include "AnaSol.h"
#include "Tensor.h"
/*!
\brief Analytic solution for Poiseuille flow

15
src/Integrands/NavierStokesG2MP.C
View File

@ -1,4 +1,4 @@
// $Id: NavierStokesG2MP.C,v 1.2 2010-10-05 09:26:12 kmo Exp $
// $Id: NavierStokesG2MP.C,v 1.3 2011-02-08 12:41:32 rho Exp $
//==============================================================================
//!
//! \file NavierStokesG2MP.C
@ -146,7 +146,6 @@ bool NavierStokesG2MP::evalInt(LocalIntegral*& elmInt,
// Stabilization parameters
if (mu < rho*U*h) {
// RUNAR
if (time.dt == 0.0)
delta1 = 0.5/sqrt(1.0/(h*h) + (U*U)/(h*h));
else
@ -161,11 +160,7 @@ bool NavierStokesG2MP::evalInt(LocalIntegral*& elmInt,
//delta2 *= 0.01*detJW;
delta1 *= 0.001*detJW;
delta2 *= 0.001*detJW;
// RUNAR
//delta1 = 0.0001*h*h*detJW;
//delta2 = 0.0001*h*h*detJW;
elmInt = &myMats;
if (eM) {
Matrix& EM = *eM;
@ -379,8 +374,7 @@ bool NavierStokesG2MP::evalInt(LocalIntegral*& elmInt,
}
}
myMats.rhsOnly = false;
return true;
return getIntegralResult(elmInt);
}
@ -439,7 +433,6 @@ bool NavierStokesG2MP::evalInt (LocalIntegral*& elmInt, double detJW,
// Stabilization parameters
if (mu < rho*U*h) {
// RUNAR
delta1 = 0.5/sqrt(1.0/(h*h) + (U*U)/(h*h));
delta2 = h;
}
@ -451,7 +444,6 @@ bool NavierStokesG2MP::evalInt (LocalIntegral*& elmInt, double detJW,
//delta1 = 0.0001*h*h*detJW;
//delta2 = 0.0001*h*h*detJW;
elmInt = &myMats;
if (eM) {
Matrix& EM = *eM;
@ -655,6 +647,5 @@ bool NavierStokesG2MP::evalInt (LocalIntegral*& elmInt, double detJW,
}
}
myMats.rhsOnly = false;
return true;
return getIntegralResult(elmInt);
}

16
src/Integrands/StabilizedNavierStokes.C
View File

@ -1,4 +1,4 @@
// $Id: StabilizedNavierStokes.C,v 1.4 2010-12-30 17:38:19 kmo Exp $
// $Id: StabilizedNavierStokes.C,v 1.5 2011-02-08 12:43:49 rho Exp $
//==============================================================================
//!
//! \file StabilizedNavierStokes.C
@ -37,8 +37,7 @@ bool StabilizedNavierStokes::evalInt(LocalIntegral*& elmInt, double detJW,
int k, l;
double div, laplace, conv;
Vector vel(3);
elmInt = &myMats;
if (eM) {
Matrix& EM = *eM;
Vector& EV = *eVs[0];
@ -106,7 +105,7 @@ bool StabilizedNavierStokes::evalInt(LocalIntegral*& elmInt, double detJW,
ES((i-1)*nf+k) += fb(k)*N(i);
}
return true;
return getIntegralResult(elmInt);
}
@ -122,7 +121,6 @@ bool StabilizedNavierStokes::evalInt (LocalIntegral*& elmInt, double detJW,
const real delta = 0.001*h*h*detJW;
elmInt = &myMats;
if (eM) {
Matrix& EM = *eM;
Vector& EV = *eVs[0];
@ -213,7 +211,7 @@ bool StabilizedNavierStokes::evalInt (LocalIntegral*& elmInt, double detJW,
ES(nf*i) -= delta*fb(k)*dNdX(i,k);
}
return true;
return getIntegralResult(elmInt);
}
@ -226,7 +224,6 @@ bool StabilizedNavierStokes::evalInt(LocalIntegral*& elmInt, double detJW,
double div, laplace, conv;
Vector vel(3);
elmInt = &myMats;
if (eM) {
Matrix& EM = *eM;
Vector& EV = *eVs[0];
@ -294,7 +291,7 @@ bool StabilizedNavierStokes::evalInt(LocalIntegral*& elmInt, double detJW,
ES((i-1)*nf+k) += fb(k)*N(i);
}
return true;
return getIntegralResult(elmInt);
}
@ -302,7 +299,6 @@ bool StabilizedNavierStokes::evalBou (LocalIntegral*& elmInt, double detJW,
const Vector& N, const Matrix& dNdX,
const Vec3& X, const Vec3& normal) const
{
elmInt = &myMats;
if (!eS || !tracFld)
{
std::cerr <<" *** StabilizedNavierStokes::evalBou: Zero pointers."<< std::endl;
@ -322,5 +318,5 @@ bool StabilizedNavierStokes::evalBou (LocalIntegral*& elmInt, double detJW,
for (short int d = 1; d <= nsd; d++)
ES(nf*(a-1)+d) += T[d-1]*N(a)*detJW;
return true;
return getIntegralResult(elmInt);
}

14
src/Integrands/StabilizedStokes.C
View File

@ -1,4 +1,4 @@
// $Id: StabilizedStokes.C,v 1.6 2010-12-30 17:38:19 kmo Exp $
// $Id: StabilizedStokes.C,v 1.7 2011-02-08 12:42:36 rho Exp $
//==============================================================================
//!
//! \file StabilizedStokes.C
@ -42,7 +42,6 @@ bool StabilizedStokes::evalInt (LocalIntegral*& elmInt, double detJW,
const int nf = nsd + 1;
elmInt = &myMats;
if (eM) {
Matrix& EM = *eM;
@ -90,8 +89,7 @@ bool StabilizedStokes::evalInt (LocalIntegral*& elmInt, double detJW,
ES((i-1)*nf+k) += fb(k)*N(i);
}
myMats.rhsOnly = false;
return true;
return getIntegralResult(elmInt);
}
@ -107,7 +105,6 @@ bool StabilizedStokes::evalInt (LocalIntegral*& elmInt, double detJW,
const int nf = nsd + 1;
const real delta = 0.001*h*h*detJW;
elmInt = &myMats;
if (eM) {
Matrix& EM = *eM;
@ -177,8 +174,7 @@ bool StabilizedStokes::evalInt (LocalIntegral*& elmInt, double detJW,
ES(nf*i) -= delta*fb(k)*dNdX(i,k);
}
myMats.rhsOnly = false;
return true;
return getIntegralResult(elmInt);
}
@ -192,7 +188,6 @@ bool StabilizedStokes::evalInt (LocalIntegral*& elmInt, double detJW,
const int nf = nsd + 1;
elmInt = &myMats;
if (eM) {
Matrix& EM = *eM;
@ -240,8 +235,7 @@ bool StabilizedStokes::evalInt (LocalIntegral*& elmInt, double detJW,
ES((i-1)*nf+k) += fb(k)*N(i);
}
myMats.rhsOnly = false;
return true;
return getIntegralResult(elmInt);
}

18
src/LinAlg/LinSolParams.C
View File

@ -1,4 +1,4 @@
// $Id: LinSolParams.C,v 1.5 2010-12-06 09:08:28 rho Exp $
// $Id: LinSolParams.C,v 1.6 2011-02-08 12:45:08 rho Exp $
//==============================================================================
//!
//! \file LinSolParams.C
@ -28,6 +28,7 @@ void LinSolParams::setDefault ()
package = MAT_SOLVER_PETSC;
levels = 0;
overlap = 0;
nullspc = NONE;
atol = 1.0e-6;
rtol = 1.0e-6;
@ -46,6 +47,7 @@ void LinSolParams::copy (const LinSolParams& spar)
package = spar.package;
levels = spar.levels;
overlap = spar.overlap;
nullspc = spar.nullspc;
atol = spar.atol;
rtol = spar.rtol;
@ -114,6 +116,16 @@ bool LinSolParams::read (std::istream& is, int nparam)
maxIts = atoi(++c);
}
else if (!strncasecmp(cline,"nullspace",6)) {
char* c = strchr(cline,'=');
if (!strncasecmp(c,"CONSTANT",8))
nullspc = CONSTANT;
else if (!strncasecmp(c,"RIGID_BODY",10))
nullspc = RIGID_BODY;
else
nullspc = NONE;
}
else
{
std::cerr <<" *** LinSolParams::read: Unknown keyword: "
@ -150,8 +162,10 @@ void LinSolParams::setParams(KSP& ksp) const
//PCMGSetLevels(pc,levels,PETSC_NULL);
//PCMGSetType(pc,PC_MG_MULTIPLICATIVE);
PCSetType(pc,prec.c_str());
if (overlap > 0)
if (overlap > 0) {
PCASMSetType(pc,PC_ASM_BASIC);
PCASMSetOverlap(pc,overlap);
}
PCFactorSetMatSolverPackage(pc,package.c_str());
PCSetFromOptions(pc);
PCSetUp(pc);

18
src/LinAlg/LinSolParams.h
View File

@ -1,4 +1,4 @@
// $Id: LinSolParams.h,v 1.4 2010-12-06 09:07:51 rho Exp $
// $Id: LinSolParams.h,v 1.5 2011-02-08 12:44:55 rho Exp $
//==============================================================================
//!
//! \file LinSolParams.h
@ -23,6 +23,9 @@
#endif
//! \brief Null-space for matrix operator
enum NullSpace { NONE, CONSTANT, RIGID_BODY };
/*!
\brief Class for linear solver parameters.
\details It contains information about solver method, preconditioner
@ -71,6 +74,12 @@ public:
//! \brief Set maximum number of iterations
virtual void setMaxIterations(int its) { maxIts = its; }
//! \brief Set number of overlap
virtual void setOverlap(int olap) { overlap = olap; }
//! \brief Set null-space of matrix
virtual void setNullSpace(NullSpace nspc) { nullspc = nspc; }
//! \brief Get linear solver method
virtual const char* getMethod() const { return method.c_str(); }
@ -92,6 +101,12 @@ public:
//! \brief Get maximum number of iterations
virtual int getMaxIterations() const { return maxIts; }
//! \brief Get number of overlaps
virtual int getOverlap() const { return overlap; }
//! \brief Get number of overlaps
virtual NullSpace getNullSpace() const { return nullspc; }
//! \brief Set linear solver parameters for KSP object
virtual void setParams(KSP& ksp) const;
@ -105,6 +120,7 @@ private:
PetscInt levels; // Number of levels of fill to use
PetscInt maxIts; // Maximum number of iterations
PetscInt overlap; // Number of overlaps
NullSpace nullspc; // Null-space for matrix
#endif // HAS_PETSC
};

33
src/LinAlg/PETScMatrix.C
View File

@ -1,4 +1,4 @@
// $Id: PETScMatrix.C,v 1.8 2010-12-06 09:19:50 rho Exp $
// $Id: PETScMatrix.C,v 1.9 2011-02-08 12:47:38 rho Exp $
//==============================================================================
//!
//! \file PETScMatrix.C
@ -63,7 +63,7 @@ size_t PETScVector::dim() const
{
PetscInt size;
size = VecGetLocalSize(x);
VecGetLocalSize(x,&size);
return size;
}
@ -152,11 +152,17 @@ real PETScVector::Linfnorm() const
PETScMatrix::PETScMatrix(const LinSolParams& spar) : solParams(spar)
{
// Create matrix object. By default the matrix type is AIJ
A = MatCreate(PETSC_COMM_WORLD);
MatCreate(PETSC_COMM_WORLD,&A);
// Create linear solver object.
KSPCreate(PETSC_COMM_WORLD,&ksp);
// Create null space if any
if (solParams.getNullSpace() == CONSTANT) {
MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nsp);
KSPSetNullSpace(ksp,nsp);
}
// Create eigenvalue solver object.
//EPSCreate(PETSC_COMM_WORLD,&eps);
}
@ -165,11 +171,18 @@ PETScMatrix::PETScMatrix(const LinSolParams& spar) : solParams(spar)
PETScMatrix::PETScMatrix (const PETScMatrix& B) : solParams(B.solParams)
{
// Duplicate matrix.
A = MatDuplicate(B.A);
//A = MatDuplicate(B.A);
MatDuplicate(B.A,MAT_COPY_VALUES,&A);
// Create linear solver object.
KSPCreate(PETSC_COMM_WORLD,&ksp);
// Create null space, if any
if (solParams.getNullSpace() == CONSTANT) {
MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nsp);
KSPSetNullSpace(ksp,nsp);
}
// Create eigenvalue solver object.
//EPSCreate(PETSC_COMM_WORLD,&eps);
}
@ -180,8 +193,12 @@ PETScMatrix::~PETScMatrix ()
// Deallocation of eigenvalue solver object.
//EPSDestroy(eps);
// Deallocation of null space
if (solParams.getNullSpace() == CONSTANT)
MatNullSpaceDestroy(nsp);
// Deallocation of linear solver object.
//KSPDestroy(ksp);
KSPDestroy(ksp);
// Deallocation of matrix object.
MatDestroy(A);
@ -516,8 +533,6 @@ bool PETScMatrix::multiply (const SystemVector& B, SystemVector& C)
bool PETScMatrix::solve (SystemVector& B, bool newLHS)
{
PC pc;
PETScVector* Bptr = dynamic_cast<PETScVector*>(&B);
if (!Bptr)
return false;
@ -526,10 +541,12 @@ bool PETScMatrix::solve (SystemVector& B, bool newLHS)
VecDuplicate(Bptr->getVector(),&x);
VecCopy(Bptr->getVector(),x);
// Has lefthand side changed?
if (newLHS)
KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);
else
KSPSetOperators(ksp,A,A,SAME_PRECONDITIONER);
solParams.setParams(ksp);
KSPSolve(ksp,x,Bptr->getVector());

3
src/LinAlg/PETScMatrix.h
View File

@ -1,4 +1,4 @@
// $Id: PETScMatrix.h,v 1.6 2010-12-06 09:17:56 rho Exp $
// $Id: PETScMatrix.h,v 1.7 2011-02-08 12:46:29 rho Exp $
//==============================================================================
//!
//! \file PETScMatrix.h
@ -208,6 +208,7 @@ public:
private:
Mat A; //!< Linear system matrix
KSP ksp; //!< Linear solver
MatNullSpace nsp; //!< Null-space of linear operator
const LinSolParams& solParams; //!< Linear solver parameters
//EPS eps; //!< Eigenvalue solver

4
src/LinAlg/SAM.h
View File

@ -1,4 +1,4 @@
// $Id: SAM.h,v 1.24 2011-02-05 18:07:51 kmo Exp $
// $Id: SAM.h,v 1.25 2011-02-08 12:45:46 rho Exp $
//==============================================================================
//!
//! \file SAM.h
@ -153,7 +153,7 @@ public:
//! \param[in] iel Identifier for the element to get the equation numbers for
//! \param[in] nedof Number of degrees of freedom in the element
//! (used for internal consistency checking, unless zero)
bool getElmEqns(IntVec& meen, int iel, int nedof = 0) const;
virtual bool getElmEqns(IntVec& meen, int iel, int nedof = 0) const;
//! \brief Finds the matrix of equation numbers for a node.
//! \param[out] mnen Matrix of node equation numbers

2
src/LinAlg/SparseMatrix.C
View File

@ -14,7 +14,7 @@
#include "SparseMatrix.h"
#include "SAM.h"
#if defined(HAS_SUPERLU_MT)
#include "pdsp_defs.h"
#include "superlu/pdsp_defs.h"
#elif defined(HAS_SUPERLU)
#include "slu_ddefs.h"
#endif

13
src/Makefile
View File

@ -1,4 +1,4 @@
# $Id: Makefile,v 1.48 2011-02-08 09:06:02 kmo Exp $
# $Id: Makefile,v 1.49 2011-02-08 13:10:06 rho Exp $
#===============================================================================
#
# File: Makefile
@ -25,12 +25,14 @@ CSRC = SIM/SIMinput.C SIM/SIMbase.C SIM/SIM3D.C SIM/SIM2D.C SIM/SIM1D.C \
ASM/ASMs2D.C ASM/ASMs2DLag.C ASM/ASMs2DSpec.C \
ASM/ASMs2Dmx.C ASM/ASMs2DmxLag.C \
ASM/ASMs1D.C ASM/ASMs1DLag.C ASM/ASMs1DSpec.C \
Integrands/AnalyticSolutions.C Integrands/LinearElasticity.C \
Integrands/Elasticity.C Integrands/Poisson.C Integrands/Stokes.C \
Integrands/AnalyticSolutions.C Integrands/AnalyticSolutionsStokes.C \
Integrands/Poisson.C \
Integrands/LinearElasticity.C Integrands/Elasticity.C Integrands/Stokes.C \
Integrands/StabilizedStokes.C Integrands/StabilizedNavierStokes.C \
Integrands/NavierStokesG2.C Integrands/NavierStokesG2CN.C \
Integrands/NavierStokesG2MP.C Integrands/NavierStokesG2GenTheta.C \
Integrands/ChorinVelPred.C Integrands/ChorinVelPredBE.C Integrands/ChorinVelPredBDF2.C \
Integrands/ChorinVelPred.C Integrands/ChorinVelPredBE.C \
Integrands/ChorinVelPredBDF2.C \
Integrands/ChorinPressCorr.C Integrands/ChorinVelCorr.C \
LinAlg/SystemMatrix.C LinAlg/DenseMatrix.C LinAlg/SPRMatrix.C \
LinAlg/SparseMatrix.C LinAlg/PETScMatrix.C LinAlg/LinSolParams.C \
@ -40,7 +42,8 @@ CSRC = SIM/SIMinput.C SIM/SIMbase.C SIM/SIM3D.C SIM/SIM2D.C SIM/SIM1D.C \
Utility/GaussQuadrature.C Utility/Legendre.C \
Utility/ElementBlock.C Utility/VTF.C \
Utility/Vec3Oper.C Utility/Tensor.C Utility/Function.C \
Utility/Utilities.C Utility/Profiler.C Utility/Functions.C
Utility/Utilities.C Utility/Profiler.C Utility/Functions.C \
Utility/AnaSol.C
OBJS = $(CSRC:.C=.o) $(FSRC:.f=.o)
HDRS = $(CSRC:.C=.h) SIM/SIMenums.h SIM/Property.h \

18
src/SIM/SIMbase.C
View File

@ -1,4 +1,4 @@
// $Id: SIMbase.C,v 1.49 2011-02-08 09:32:18 kmo Exp $
// $Id: SIMbase.C,v 1.50 2011-02-08 12:52:12 rho Exp $
//==============================================================================
//!
//! \file SIMbase.C
@ -268,11 +268,10 @@ bool SIMbase::preprocess (const std::vector<int>& ignoredPatches, bool fixDup)
// Initialize data structures for the algebraic system
#ifdef PARALLEL_PETSC
if (nProc > 1)
mySam = new SAMpatchPara(l2gn);
else
#endif
#else
mySam = new SAMpatch();
#endif
return mySam->init(myModel,nnod) && ok;
}
@ -425,7 +424,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++)
@ -600,7 +599,6 @@ bool SIMbase::solutionNorms (const TimeDomain& time, const Vectors& psol,
Matrix& eNorm, Vector& gNorm)
{
NormBase* norm = myProblem->getNormIntegrand(this->getAnaSol());
if (!norm) norm = myProblem->getNormIntegrandScal(this->getScalarSol());
if (!norm)
{
std::cerr <<" *** SIMbase::solutionNorms: No integrand."<< std::endl;
@ -609,9 +607,9 @@ bool SIMbase::solutionNorms (const TimeDomain& time, const Vectors& psol,
PROFILE1("Norm integration");
size_t nNorms = norm->getNoFields() + norm->getNoSecFields();
const Vector& primsol = psol.front();
size_t nNorms = norm->getNoFields();
size_t i, nel = mySam->getNoElms();
eNorm.resize(nNorms,nel,true);
gNorm.resize(nNorms,true);
@ -925,7 +923,7 @@ bool SIMbase::writeGlvS (const Vector& psol,
bool haveAsol = false;
bool scalarEq = myModel.empty() ? false : myModel.front()->getNoFields() == 1;
if (scalarEq) {
if (this->getScalarSol())
if (this->getAnaSol()->hasScalarSol())
haveAsol = true;
}
else
@ -993,9 +991,9 @@ bool SIMbase::writeGlvS (const Vector& psol,
for (j = 1; cit != grid->end_XYZ() && ok; j++, cit++)
{
if (scalarEq)
ok = myProblem->evalSolScal(solPt,*this->getScalarSol(),*cit);
ok = myProblem->evalSecSolScal(solPt,*this->getAnaSol()->getScalarSecSol(),*cit);
else
ok = myProblem->evalSol(solPt,*this->getAnaSol(),*cit);
ok = myProblem->evalSecSol(solPt,*this->getAnaSol()->getVectorSecSol(),*cit);
if (ok)
field.fillColumn(j,solPt);
}

22
src/SIM/SIMbase.h
View File

@ -1,4 +1,4 @@
// $Id: SIMbase.h,v 1.33 2011-02-08 09:32:18 kmo Exp $
// $Id: SIMbase.h,v 1.35 2011-02-08 15:51:16 rho Exp $
//==============================================================================
//!
//! \file SIMbase.h
@ -19,6 +19,7 @@
#include "TimeDomain.h"
#include "Property.h"
#include "Function.h"
#include "AnaSol.h"
#include <map>
class ASMbase;
@ -97,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,27 +170,27 @@ public:
//! \param[out] inf Infinity norms in each spatial direction
//! \param[out] ind Global index of the node corresponding to the inf-value
//! \return L2-norm of the solution vector
double solutionNorms(const Vector& x, double* inf, size_t* ind) const;
virtual double solutionNorms(const Vector& x, double* inf, size_t* ind) const;
//! \brief Integrates some solution norm quantities.
//! \details If an analytical solution is provided, norms of the exact
//! error in the solution are computed as well.
//! \param[in] time Parameters for nonlinear/time-dependent simulations.
//! \param[in] psol Global primary solution vector
//! \param[in] psol Global primary solution vectors
//! \param[out] eNorm Element-wise norm quantities
//! \param[out] gNorm Global norm quantities
bool solutionNorms(const TimeDomain& time, const Vectors& psol,
Matrix& eNorm, Vector& gNorm);
virtual bool solutionNorms(const TimeDomain& time, const Vectors& psol,
Matrix& eNorm, Vector& gNorm);
//! \brief Integrates some solution norm quantities.
//! \details If an analytical solution is provided, norms of the exact
//! error in the solution are computed as well.
//! \param[in] psol Global primary solution vector
//! \param[in] psol Global primary solution vectors
//! \param[out] eNorm Element-wise norm quantities
//! \param[out] gNorm Global norm quantities
//!
//! \details Use this version for linear/stationary problems only.
bool solutionNorms(const Vectors& psol, Matrix& eNorm, Vector& gNorm)
virtual bool solutionNorms(const Vectors& psol, Matrix& eNorm, Vector& gNorm)
{ return this->solutionNorms(TimeDomain(),psol,eNorm,gNorm); }
//! \brief Performs a generalized eigenvalue analysis of the assembled system.
@ -322,10 +323,7 @@ protected:
virtual bool initNeumann(size_t) { return true; }
//! \brief Returns the analytical solution field, if any.
virtual TensorFunc* getAnaSol() const { return 0; }
//! \brief Returns the scalar analytical solution field, if any.
virtual VecFunc* getScalarSol() const { return 0; }
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

39
src/Utility/AnaSol.C Normal file
View File

@ -0,0 +1,39 @@
//==============================================================================
//!
//! \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 Normal file
View File

@ -0,0 +1,72 @@
//==============================================================================
//!
//! \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

8
src/Utility/Function.h
View File

@ -1,4 +1,4 @@
// $Id: Function.h,v 1.6 2010-09-16 10:19:08 kmo Exp $
// $Id: Function.h,v 1.7 2011-02-08 15:07:19 rho Exp $
//==============================================================================
//!
//! \file Function.h
@ -71,8 +71,10 @@ namespace utl
class Vec3;
class Tensor;
class SymmTensor;
//! \brief Scalar-valued unary function of a scalar value.
typedef utl::Function<real,real> ScalarFunc;
@ -82,8 +84,8 @@ typedef utl::Function<Vec3,real> RealFunc;
//! \brief Vector-valued unary function of a spatial point.
typedef utl::Function<Vec3,Vec3> VecFunc;
//! \brief Tensor-valued unary function of a spatial point.
typedef utl::Function<Vec3,SymmTensor> TensorFunc;
//! \brief Symmetric tensor-valued unary function of a spatial point.
typedef utl::Function<Vec3,Tensor> TensorFunc;
/*!

56
src/Utility/Functions.C
View File

@ -1,4 +1,4 @@
// $Id: Functions.C,v 1.7 2011-01-05 12:45:14 kmo Exp $
// $Id: Functions.C,v 1.8 2011-02-08 12:55:52 rho Exp $
//==============================================================================
//!
//! \file Functions.C
@ -65,6 +65,27 @@ real LinearZFunc::evaluate (const Vec3& X) const
}
real QuadraticXFunc::evaluate(const Vec3& X) const
{
real val = 0.5*(a-b);
return max*(a-X.x)*(X.x-b)/(val*val);
}
real QuadraticYFunc::evaluate(const Vec3& X) const
{
real val = 0.5*(a-b);
return max*(a-X.y)*(X.y-b)/(val*val);
}
real QuadraticZFunc::evaluate(const Vec3& X) const
{
real val = 0.5*(a-b);
return max*(a-X.z)*(X.z-b)/(val*val);
}
real LinearRotZFunc::evaluate (const Vec3& _X) const
{
// Always return zero if the argument has no time component
@ -95,7 +116,8 @@ real StepXYFunc::evaluate (const Vec3& X) const
const RealFunc* utl::parseRealFunc (char* cline, real A)
{
// Check for spatial variation
int linear = 0;
int linear = 0;
int quadratic = 0;
if (cline)
if (strcmp(cline,"X") == 0)
linear = 1;
@ -109,6 +131,12 @@ const RealFunc* utl::parseRealFunc (char* cline, real A)
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)
@ -155,8 +183,28 @@ const RealFunc* utl::parseRealFunc (char* cline, real A)
}
cline = strtok(NULL," ");
}
else
std::cout << C;
else if (quadratic && (cline = strtok(NULL," ")))
{
real a = atof(cline);
real b = atof(strtok(NULL," "));
real val = 0.5*(a-b);
val *= val;
std::cout << A/val << "*(" << char('W' + quadratic) << "-" << a << ")*("
<< b << "-" << char('W' + quadratic) << ")";
switch(quadratic) {
case 1:
f = new QuadraticXFunc(A,a,b);
break;
case 2:
f = new QuadraticYFunc(A,a,b);
break;
case 3:
f = new QuadraticZFunc(A,a,b);
break;
}
cline = strtok(NULL," ");
}
//std::cout << C;
// Check for time variation
if (!cline) return f;

59
src/Utility/Functions.h
View File

@ -1,4 +1,4 @@
// $Id: Functions.h,v 1.6 2011-01-05 12:45:14 kmo Exp $
// $Id: Functions.h,v 1.7 2011-02-08 12:55:33 rho Exp $
//==============================================================================
//!
//! \file Functions.h
@ -194,6 +194,63 @@ protected:
};
/*!
\brief A scalar function, quadratic in \a x.
*/
class QuadraticXFunc : public RealFunc
{
real max; // Max value of function
real a, b; // Roots where function is \a 0
public:
//! \brief Constructor initializing the function parameters.
QuadraticXFunc(real MAX, real A, real B) { max = MAX; a = A; b = B; }
protected:
//! \brief Evaluates the quadratic function.
virtual real evaluate(const Vec3& X) const;
};
/*!
\brief A scalar function, quadratic in \a y.
*/
class QuadraticYFunc : public RealFunc
{
real max; // Max value of function
real a, b; // Roots where function is \a 0
public:
//! \brief Constructor initializing the function parameters.
QuadraticYFunc(real MAX, real A, real B) { max = MAX; a = A; b = B; }
protected:
//! \brief Evaluates the quadratic function.
virtual real evaluate(const Vec3& X) const;
};
/*!
\brief A scalar function, quadratic in \a x.
*/
class QuadraticZFunc : public RealFunc
{
real max; // Max value of function
real a, b; // Roots where function is \a 0
public:
//! \brief Constructor initializing the function parameters.
QuadraticZFunc(real MAX, real A, real B) { max = MAX; a = A; b = B; }
protected:
//! \brief Evaluates the quadratic function.
virtual real evaluate(const Vec3& X) const;
};
/*!
\brief A scalar function, defining a linear rotation about the global Z-axis.
\details The time component of the function argument multiplied with the