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

Some corrections in iteration norm evaluation for mixed methods. Added output of total reaction forces and norm (R,u)

Browse Source 
git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@913 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
kmo
2011-04-14 15:24:34 +00:00
committed by Knut Morten Okstad
parent fcaff5a760
commit 70c071c48e
18 changed files with 265 additions and 124 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
src/ASM/ASMbase.h
View File

@@ -99,7 +99,7 @@ public:
//! \brief Returns the number of parameter dimensions.
unsigned char getNoParamDim() const { return ndim; }
//! \brief Returns the number of solution fields.
virtual unsigned char getNoFields(int = 0) const { return nf; }
virtual unsigned char getNoFields(int b = 0) const { return b > 1 ? 0 : nf; }
//! \brief Returns local 1-based index of the node with given global number.
//! \details If the given node number is not present, 0 is returned.
@@ -113,6 +113,8 @@ public:
int getElmID(size_t iel) const;
//! \brief Returns the number of DOFs per node.
virtual unsigned char getNodalDOFs(size_t) const { return nf; }
//! \brief Returns which mixed field basis a node belongs to.
virtual unsigned char getNodalBasis(size_t) const { return 0; }
//! \brief Returns the global coordinates for the given node.
//! \param[in] inod 1-based node index local to current patch
virtual Vec3 getCoord(size_t inod) const = 0;

6
src/ASM/ASMs2Dmx.C
View File

@@ -74,6 +74,12 @@ unsigned char ASMs2Dmx::getNodalDOFs (size_t inod) const
}
unsigned char ASMs2Dmx::getNodalBasis (size_t inod) const
{
return inod <= nb1 ? 1 : 2;
}
void ASMs2Dmx::initMADOF (const int* sysMadof)
{
this->initMx(MLGN,sysMadof);

3
src/ASM/ASMs2Dmx.h
View File

@@ -61,6 +61,9 @@ public:
//! \brief Returns the number of DOFs per node.
//! \param[in] inod 1-based node index local to current patch
virtual unsigned char getNodalDOFs(size_t inod) const;
//! \brief Returns which mixed field basis a node belongs to.
//! \param[in] inod 1-based node index local to current patch
virtual unsigned char getNodalBasis(size_t inod) const;
//! \brief Initializes the patch level MADOF array for mixed problems.
virtual void initMADOF(const int* sysMadof);

6
src/ASM/ASMs2DmxLag.C
View File

@@ -68,6 +68,12 @@ unsigned char ASMs2DmxLag::getNodalDOFs (size_t inod) const
}
unsigned char ASMs2DmxLag::getNodalBasis (size_t inod) const
{
return inod <= nb1 ? 1 : 2;
}
void ASMs2DmxLag::initMADOF (const int* sysMadof)
{
this->initMx(MLGN,sysMadof);

6
src/ASM/ASMs2DmxLag.h
View File

@@ -34,9 +34,6 @@ public:
//! \brief Constructor creating an instance by reading the given input stream.
ASMs2DmxLag(std::istream& is, unsigned char n_s = 2,
unsigned char n_f1 = 2, unsigned char n_f2 = 1);
//! \brief Default constructor creating an empty patch.
ASMs2DmxLag(unsigned char n_s = 2,
unsigned char n_f1 = 2, unsigned char n_f2 = 1);
//! \brief Empty destructor.
virtual ~ASMs2DmxLag() {}
@@ -57,6 +54,9 @@ public:
//! \brief Returns the number of DOFs per node.
//! \param[in] inod 1-based node index local to current patch
virtual unsigned char getNodalDOFs(size_t inod) const;
//! \brief Returns which mixed field basis a node belongs to.
//! \param[in] inod 1-based node index local to current patch
virtual unsigned char getNodalBasis(size_t inod) const;
//! \brief Initializes the patch level MADOF array for mixed problems.
virtual void initMADOF(const int* sysMadof);

6
src/ASM/ASMs3Dmx.C
View File

@@ -74,6 +74,12 @@ unsigned char ASMs3Dmx::getNodalDOFs (size_t inod) const
}
unsigned char ASMs3Dmx::getNodalBasis (size_t inod) const
{
return inod <= nb1 ? 1 : 2;
}
void ASMs3Dmx::initMADOF (const int* sysMadof)
{
this->initMx(MLGN,sysMadof);

3
src/ASM/ASMs3Dmx.h
View File

@@ -61,6 +61,9 @@ public:
//! \brief Returns the number of DOFs per node.
//! \param[in] inod 1-based node index local to current patch
virtual unsigned char getNodalDOFs(size_t inod) const;
//! \brief Returns which mixed field basis a node belongs to.
//! \param[in] inod 1-based node index local to current patch
virtual unsigned char getNodalBasis(size_t inod) const;
//! \brief Initializes the patch level MADOF array for mixed problems.
virtual void initMADOF(const int* sysMadof);

6
src/ASM/ASMs3DmxLag.C
View File

@@ -68,6 +68,12 @@ unsigned char ASMs3DmxLag::getNodalDOFs (size_t inod) const
}
unsigned char ASMs3DmxLag::getNodalBasis (size_t inod) const
{
return inod <= nb1 ? 1 : 2;
}
void ASMs3DmxLag::initMADOF (const int* sysMadof)
{
this->initMx(MLGN,sysMadof);

6
src/ASM/ASMs3DmxLag.h
View File

@@ -34,9 +34,6 @@ public:
//! \brief Constructor creating an instance by reading the given input stream.
ASMs3DmxLag(std::istream& is, bool checkRHS = false,
unsigned char n_f1 = 3, unsigned char n_f2 = 1);
//! \brief Default constructor creating an empty patch.
ASMs3DmxLag(bool checkRHS = false,
unsigned char n_f1 = 3, unsigned char n_f2 = 1);
//! \brief Empty destructor.
virtual ~ASMs3DmxLag() {}
@@ -57,6 +54,9 @@ public:
//! \brief Returns the number of DOFs per node.
//! \param[in] inod 1-based node index local to current patch
virtual unsigned char getNodalDOFs(size_t inod) const;
//! \brief Returns which mixed field basis a node belongs to.
//! \param[in] inod 1-based node index local to current patch
virtual unsigned char getNodalBasis(size_t inod) const;
//! \brief Initializes the patch level MADOF array for mixed problems.
virtual void initMADOF(const int* sysMadof);

52
src/ASM/SAMpatch.C
View File

@@ -30,7 +30,8 @@ bool SAMpatch::init (const ASMVec& model, int numNod)
}
// Initialize the node/dof arrays (madof,msc) and compute ndof
this->initNodeDofs(model);
if (!this->initNodeDofs(model))
return false;
std::cout <<"\n\n >>> SAM model summary <<<"
<<"\nNumber of elements "<< nel
@@ -38,11 +39,13 @@ bool SAMpatch::init (const ASMVec& model, int numNod)
<<"\nNumber of dofs "<< ndof << std::endl;
// Initialize the element connectivity arrays (mpmnpc,mmnpc)
this->initElementConn(model);
if (!this->initElementConn(model))
return false;
// Initialize the constraint equation arrays (mpmceq,mmceq,ttcc)
this->initConstraintEqs(model);
if (nceq > 0)
if (!this->initConstraintEqs(model))
return false;
else if (nceq > 0)
std::cout <<"Number of constraints "<< nceq << std::endl;
// Initialize the dof-to-equation connectivity array (meqn)
@@ -52,21 +55,36 @@ bool SAMpatch::init (const ASMVec& model, int numNod)
}
void SAMpatch::initNodeDofs (const ASMVec& model)
bool SAMpatch::initNodeDofs (const ASMVec& model)
{
if (nnod < 1) return;
if (nnod < 1) return true;
int n;
size_t i, j;
// Check if we are doing mixed methods
for (i = 0; i < model.size() && nodeType.empty(); i++)
if (!model[i]->empty() && model[i]->getNoFields(2) > 0)
nodeType.resize(nnod,'D'); // Initialize all nodes as 'Displacement type'
// Initialize the array of accumulated DOFs for the nodes
madof = new int[nnod+1];
memset(madof,0,(nnod+1)*sizeof(int));
bool ierr = 0;
for (i = 0; i < model.size(); i++)
for (j = 0; j < model[i]->getNoNodes(); j++)
if ((n = model[i]->getNodeID(j+1)) > 0 && n <= nnod)
madof[n] = model[i]->getNodalDOFs(j+1);
{
if (madof[n] == 0)
madof[n] = model[i]->getNodalDOFs(j+1);
else if (madof[n] != model[i]->getNodalDOFs(j+1))
ierr++;
// Define the node type for mixed methods (used by norm evaluations)
if (!nodeType.empty() && model[i]->getNodalBasis(j+1) == 2)
nodeType[n-1] = 'P'; // This is a 'Pressure type' node (basis 2)
}
madof[0] = 1;
for (n = 0; n < nnod; n++)
@@ -91,12 +109,19 @@ void SAMpatch::initNodeDofs (const ASMVec& model)
if (nndof > 1) msc[idof1 ] *= bit->CY;
if (nndof > 2) msc[idof1+1] *= bit->CZ;
}
if (ierr == 0) return true;
std::cerr <<" *** SAMpatch::initNodeDOFs: Detected "<< ierr <<" nodes"
<<" nodes with conflicting number of DOFs in adjacent patches."
<< std::endl;
return false;
}
void SAMpatch::initElementConn (const ASMVec& model)
bool SAMpatch::initElementConn (const ASMVec& model)
{
if (nel < 1) return;
if (nel < 1) return true;
// Find the size of the element connectivity array
size_t i, j;
@@ -119,10 +144,12 @@ void SAMpatch::initElementConn (const ASMVec& model)
mmnpc[(mpmnpc[ip]++)-1] = model[j]->getNodeID(1+(*eit)[i]);
ip++;
}
return true;
}
void SAMpatch::initConstraintEqs (const ASMVec& model)
bool SAMpatch::initConstraintEqs (const ASMVec& model)
{
// Estimate the size of the constraint equation array
size_t j;
@@ -134,7 +161,8 @@ void SAMpatch::initConstraintEqs (const ASMVec& model)
// Initialize the constraint equation arrays
mpmceq = new int[nceq+1];
int ip = mpmceq[0] = 1;
if (nceq < 1) return;
if (nceq < 1) return true;
mmceq = new int[nmmceq];
ttcc = new real[nmmceq];
for (j = 0; j < model.size(); j++)
@@ -198,6 +226,8 @@ void SAMpatch::initConstraintEqs (const ASMVec& model)
// Reset the negative values in msc before calling SYSEQ
for (ip = 0; ip < ndof; ip++)
if (msc[ip] < 0) msc[ip] = 0;
return true;
}

8
src/ASM/SAMpatch.h
View File

@@ -1,4 +1,4 @@
// $Id: SAMpatch.h,v 1.5 2011-01-02 16:33:04 kmo Exp $
// $Id$
//==============================================================================
//!
//! \file SAMpatch.h
@@ -46,12 +46,12 @@ public:
protected:
//! \brief Initializes the nodal arrays \a MINEX, \a MADOF and \a MSC.
void initNodeDofs(const std::vector<ASMbase*>& model);
bool initNodeDofs(const std::vector<ASMbase*>& model);
//! \brief Initializes the element topology arrays \a MPMNPC and \a MMNPC.
void initElementConn(const std::vector<ASMbase*>& model);
bool initElementConn(const std::vector<ASMbase*>& model);
//! \brief Initializes the multi-point constraint arrays
//! \a MPMCEQ, \a MMCEQ and \a TTCC.
virtual void initConstraintEqs(const std::vector<ASMbase*>& model);
virtual bool initConstraintEqs(const std::vector<ASMbase*>& model);
};
#endif

73
src/ASM/SAMpatchPara.C
View File

@@ -13,13 +13,12 @@
#include "SAMpatchPara.h"
#include "PETScMatrix.h"
#include "LinAlgInit.h"
#include "ASMbase.h"
#include "MPC.h"
#include "LinAlgInit.h"
SAMpatchPara::SAMpatchPara(const IntVec& l2gn_mp)
SAMpatchPara::SAMpatchPara (const IntVec& l2gn_mp)
{
l2gn = l2gn_mp;
#ifdef PARALLEL_PETSC
@@ -94,10 +93,8 @@ bool SAMpatchPara::getNoDofCouplings (int ifirst, int ilast,
}
// Generate nnz for off-diagonal block
IntVec l2g;
Vector nnz;
l2g.resize(ndof);
nnz.resize(ndof);
IntVec l2g(ndof);
Vector nnz(ndof);
for (i = 0;i < ndof;i++) {
l2g[i] = meqn[i]-1;
nnz[i] = o_dofc[i].size();
@@ -130,12 +127,13 @@ bool SAMpatchPara::getNoDofCouplings (int ifirst, int ilast,
}
bool SAMpatchPara::initForAssembly (SystemVector& sysRHS, Vector* reactionForces) const
bool SAMpatchPara::initForAssembly (SystemVector& sysRHS,
Vector* reactionForces) const
{
sysRHS.redim(nleq);
sysRHS.init();
if (reactionForces)
reactionForces->resize(mpar[5],true);
reactionForces->resize(nspdof,true);
return nleq > 0 ? true : false;
}
@@ -243,9 +241,9 @@ bool SAMpatchPara::expandSolution (const SystemVector& solVec,
}
real SAMpatchPara::dot (const Vector& x, const Vector& y) const
real SAMpatchPara::dot (const Vector& x, const Vector& y, char dofType) const
{
real globVal = this->SAM::dot(x,y);
real globVal = this->SAM::dot(x,y,dofType);
#ifdef PARALLEL_PETSC
if (nProc > 1)
@@ -254,7 +252,7 @@ real SAMpatchPara::dot (const Vector& x, const Vector& y) const
for (size_t i = 0; i < ghostNodes.size(); i++)
{
int inod = ghostNodes[i];
if (madof[inod] - madof[inod-1] == mpar[17])
if (nodeType[inod-1] == dofType || dofType == 'A')
for (int j = madof[inod-1]; j < madof[inod]; j++)
locVal -= x(j)*y(j);
}
@@ -267,35 +265,44 @@ real SAMpatchPara::dot (const Vector& x, const Vector& y) const
}
real SAMpatchPara::normL2 (const Vector& x) const
real SAMpatchPara::normL2 (const Vector& x, char dofType) const
{
#ifdef PARALLEL_PETSC
if (nProc > 1 && nnodGlob > 1)
return this->norm2(x)/sqrt(mpar[17]*nnodGlob);
return this->norm2(x,dofType)/sqrt((madof[1]-madof[0])*nnodGlob);
// TODO,kmo: The above is not correct for mixed methods. We need to find the
// global number of DOFs of type dofType and use that in the denominator.
#endif
return this->SAM::normL2(x);
return this->SAM::normL2(x,dofType);
}
real SAMpatchPara::normInf (const Vector& x, size_t& comp) const
real SAMpatchPara::normInf (const Vector& x, size_t& comp, char dofType) const
{
real locmax = this->SAM::normInf(x,comp);
real locmax = this->SAM::normInf(x,comp,dofType);
#ifdef PARALLEL_PETSC
int nProc, myRank;
MPI_Comm_size(PETSC_COMM_WORLD,&nProc);
MPI_Comm_rank(PETSC_COMM_WORLD,&myRank);
if (nProc > 1)
{
int myRank;
MPI_Comm_rank(PETSC_COMM_WORLD,&myRank);
if (nProc > 1) {
Vector locval, globval;
locval.resize(2*nProc,0.0);
globval.resize(2*nProc,0.0);
int nndof = madof[1]-madof[0];
for (size_t i = 0; i < nodeType.size(); i++)
if (nodeType[i] == dofType)
{
nndof = madof[i+1]-madof[i];
break;
}
comp = meqn[(comp-1)*mpar[17]]/mpar[17]+1;
// TODO,kmo: Don't think this is correct in case of mixed methods
comp = meqn[(comp-1)*nndof]/nndof+1;
RealArray locval(2*nProc,0.0), globval(2*nProc,0.0);
locval[2*myRank] = locmax;
locval[2*myRank+1] = 1.0*comp;
MPI_Allreduce(&locval[0],&globval[0],2*nProc,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD);
// TODO,kmo: Is this calculation of comp correct? I doubth it...
for (int n = 0; n < nProc; n++)
if (globval[2*n] > locmax) {
locmax = globval[2*n];
@@ -307,7 +314,7 @@ real SAMpatchPara::normInf (const Vector& x, size_t& comp) const
}
void SAMpatchPara::initConstraintEqs (const std::vector<ASMbase*>& model)
bool SAMpatchPara::initConstraintEqs (const std::vector<ASMbase*>& model)
{
// TODO: Rewrite this calling the parent-class method, and then replacing
// the mpmceq array into the ndof-sized array used here (why is it needed?)
@@ -322,7 +329,8 @@ void SAMpatchPara::initConstraintEqs (const std::vector<ASMbase*>& model)
// Initialize the constraint equation arrays
mpmceq = new int[ndof];
memset(mpmceq,0,ndof*sizeof(int));
if (nceq < 1) return;
if (nceq < 1) return true;
mmceq = new int[nmmceq];
ttcc = new real[nmmceq];
int ip = 1;
@@ -386,6 +394,8 @@ void SAMpatchPara::initConstraintEqs (const std::vector<ASMbase*>& model)
// Reset the negative values in msc before calling SYSEQ
for (ip = 0; ip < ndof; ip++)
if (msc[ip] < 0) msc[ip] = 0;
return true;
}
@@ -492,14 +502,5 @@ bool SAMpatchPara::initSystemEquations ()
// Number of equations equals number of dofs
neq = ndof;
// Initialize the number of nodal DOFs
int inod, nndof;
mpar[16] = mpar[17] = madof[1]-madof[0];
for (inod = 1; inod < nnod && madof; inod++)
if ((nndof = madof[1]-madof[0]) < mpar[16])
mpar[16] = nndof;
else if (nndof > mpar[17])
mpar[17] = nndof;
return true;
}

22
src/ASM/SAMpatchPara.h
View File

@@ -1,4 +1,4 @@
// $Id: SAMpatchPara.h,v 1.4 2011-02-08 12:09:30 rho Exp $
// $Id$
//==============================================================================
//!
//! \file SAMpatchPara.h
@@ -49,7 +49,8 @@ public:
//! \param sysRHS The system right-hand-side load vector to be initialized
//! \param reactionForces Pointer to vector of nodal reaction forces
//! \return \e false if no free DOFs in the system, otherwise \e true
virtual bool initForAssembly(SystemVector& sysRHS, Vector* reactionForces = 0) const;
virtual bool initForAssembly(SystemVector& sysRHS,
Vector* reactionForces = 0) const;
//! \brief Adds element stiffness contributions to the system load vector.
//! \param sysRHS The right-hand-side system load vector
@@ -98,42 +99,43 @@ public:
virtual bool expandSolution(const SystemVector& solVec, Vector& displ) const;
//! \brief Computes the dot-product of two vectors.
//! \param[in] x First vector in dot-product.
//! \param[in] y Second vector in dot-product.
//! \param[in] x First vector in dot-product
//! \param[in] y Second vector in dot-product
//! \param[in] dofType Only consider nodes of this type (for mixed methods)
//! \return Value of dot-product
//!
//! \details Both vectors are defined over all nodes in the patch, i.e.
//! for parallel vectors the ghost entries are also included.
virtual real dot(const Vector& x, const Vector& y) const;
virtual real dot(const Vector& x, const Vector& y, char dofType = 'D') const;
//! \brief Computes the L2-norm of a vector.
//! \param[in] x Vector for norm computation
//! \param[in] dofType Only consider nodes of this type (for mixed methods)
//! \return Value of L2-norm
//!
//! \details The vector is defined over all nodes in the patch, i.e.
//! for parallel vectors the ghost entries are also included.
virtual real normL2(const Vector& x) const;
virtual real normL2(const Vector& x, char dofType = 'D') const;
//! \brief Computes the inf-norm of a vector.
//! \param[in] x Vector for norm computation
//! \param comp Local nodal DOF on input, index of the largest value on output
//! \param[in] dofType Only consider nodes of this type (for mixed methods)
//! \return Value of inf-norm
//!
//! \details The vector is defined over all nodes in the patch, i.e.
//! for parallel vectors the ghost entries are also included.
virtual real normInf(const Vector& x, size_t& comp) const;
virtual real normInf(const Vector& x, size_t& comp, char dofType = 'D') const;
protected:
//! \brief Initializes the multi-point constraint arrays
//! \a MPMCEQ, \a MMCEQ and \a TTCC.
virtual void initConstraintEqs(const std::vector<ASMbase*>& model);
virtual bool initConstraintEqs(const std::vector<ASMbase*>& model);
//! \brief Initializes the DOF-to-equation connectivity array \a MEQN.
virtual bool initSystemEquations();
private:
typedef std::set<int> IntSet; //!< General integer set
// Parameters for parallel computing
int nProc; //!< Number of processes
int nleq; //!< Number of equations for this processor

84
src/LinAlg/SAM.C
View File

@@ -63,7 +63,7 @@ void expand_(const real* solVec, const real* ttcc,
SAM::SAM () : nnod(mpar[0]), nel(mpar[1]), ndof(mpar[2]),
nceq(mpar[6]), neq(mpar[10]),
nspdof(mpar[5]), nceq(mpar[6]), neq(mpar[10]),
nmmnpc(mpar[14]), nmmceq(mpar[15])
{
// Initialize the parameters array to zero
@@ -166,6 +166,8 @@ void SAM::print (std::ostream& os) const
}
code[i] = '\0';
os <<'\t'<< code;
if ((size_t)n < nodeType.size())
os <<'\t'<< nodeType[n];
}
os << std::endl;
}
@@ -187,7 +189,6 @@ bool SAM::initSystemEquations ()
int i, j, idof;
int ndof1 = 0;
int ndof2 = 0;
int nspdof = 0;
int npdof = 0;
int nddof = 0;
for (i = 0; i < ndof; i++)
@@ -246,7 +247,6 @@ bool SAM::initSystemEquations ()
mpar[3] = ndof1;
mpar[4] = ndof2;
mpar[5] = nspdof;
mpar[7] = nspdof - nceq;
mpar[8] = npdof;
mpar[9] = nddof;
@@ -261,15 +261,6 @@ bool SAM::initSystemEquations ()
meqn[idof] = j++;
#endif
// Initialize the number of nodal DOFs
int inod, nndof;
mpar[16] = mpar[17] = madof[1]-madof[0];
for (inod = 1; inod < nnod && madof; inod++)
if ((nndof = madof[inod+1]-madof[inod]) < mpar[16])
mpar[16] = nndof;
else if (nndof > mpar[17])
mpar[17] = nndof;
if (ierr == 0) return true;
std::cerr <<"SAM::initSystemEquations: Failure "<< ierr << std::endl;
@@ -382,7 +373,7 @@ bool SAM::initForAssembly (SystemVector& sysRHS, Vector* reactionForces) const
sysRHS.redim(neq);
sysRHS.init();
if (reactionForces)
reactionForces->resize(mpar[5],true);
reactionForces->resize(nspdof,true);
return neq > 0 ? true : false;
}
@@ -655,16 +646,16 @@ bool SAM::expandVector (const real* solVec, Vector& dofVec) const
}
real SAM::dot (const Vector& x, const Vector& y) const
real SAM::dot (const Vector& x, const Vector& y, char dofType) const
{
if (mpar[16] == mpar[17]) // All nodes have the same number of DOFs
return x.dot(y);
if (nodeType.empty() || dofType == 'A')
return x.dot(y); // All nodes are of the same type, or consider all of them
// Consider the nodes with mpar[17] DOFs only
// Consider only the dofType nodes
int i, j, n = x.size() < y.size() ? x.size() : y.size();
real retVal = real(0);
for (i = 0; i < nnod; i++)
if (madof[i+1] - madof[i] == mpar[17])
if (nodeType[i] == dofType)
for (j = madof[i]-1; j < madof[i+1] && j < n; j++)
retVal += x[j]*y[j];
@@ -672,19 +663,19 @@ real SAM::dot (const Vector& x, const Vector& y) const
}
real SAM::normL2 (const Vector& x) const
real SAM::normL2 (const Vector& x, char dofType) const
{
if (x.empty())
return real(0);
else if (mpar[16] == mpar[17]) // All nodes have the same number of DOFs
else if (nodeType.empty()) // All nodes are of the same type
return x.norm2()/sqrt(x.size());
// Consider the nodes with mpar[17] DOFs only
int i, j, count = 0;
// Consider only the dofType nodes
int count, i, j, n = x.size();
real retVal = real(0);
for (i = 0; i < nnod; i++)
if (madof[i+1] - madof[i] == mpar[17])
for (j = madof[i]-1; j < madof[i+1] && j < (int)x.size(); j++)
for (count = i = 0; i < nnod; i++)
if (nodeType[i] == dofType)
for (j = madof[i]-1; j < madof[i+1] && j < n; j++)
{
retVal += x[j]*x[j];
count ++;
@@ -694,21 +685,28 @@ real SAM::normL2 (const Vector& x) const
}
real SAM::normInf (const Vector& x, size_t& comp) const
real SAM::normInf (const Vector& x, size_t& comp, char dofType) const
{
if (x.empty() || (int)comp > mpar[17])
if (x.empty() || comp < 1)
return real(0);
else if (mpar[16] == mpar[17]) // All nodes have the same number of DOFs
return x.normInf(--comp,mpar[17]);
else if (nodeType.empty())
{
// All nodes are of the same type with the same number of nodal DOFs
int nndof = madof[1] - madof[0];
if ((int)comp <= nndof)
return x.normInf(--comp,nndof);
else
return real(0);
}
// Consider the nodes with mpar[17] DOFs only
// Consider only the dofType nodes
int i, k = comp;
real retVal = real(0);
for (i = 0; i < nnod; i++)
if (madof[i+1] - madof[i] == mpar[17])
if (nodeType[i] == dofType)
{
int idof = madof[i] + k-2;
if (idof < (int)x.size())
if (idof >= 0 && idof < madof[i+1]-1)
if (fabs(x[idof]) > retVal)
{
comp = i+1;
@@ -723,9 +721,26 @@ real SAM::normInf (const Vector& x, size_t& comp) const
real SAM::normReact (const Vector& u, const Vector& rf) const
{
real retVal = real(0);
for (int i = 0; i < ndof; i++)
if (msc[i] < 0 && -msc[i] <= (int)rf.size())
retVal += 0.5*u[i]*rf(-msc[i]);
retVal += u[i]*rf(-msc[i]);
return 0.5*retVal;
}
real SAM::getReaction (int dir, const Vector& rf) const
{
real retVal = real(0);
if (dir > 0)
for (int i = 0; i < nnod; i++)
{
int idof = madof[i]+dir-2;
if (idof < madof[i+1]-1 && msc[idof] < 0 && -msc[idof] <= (int)rf.size())
retVal += rf(-msc[idof]);
}
return retVal;
}
@@ -733,7 +748,8 @@ real SAM::normReact (const Vector& u, const Vector& rf) const
bool SAM::getNodalReactions (int inod, const Vector& rf, Vector& nrf) const
{
if (inod < 1 || inod > nnod) return false;
if (inod < 1 || inod > nnod)
return false;
bool haveRF = false;
int ip = madof[inod-1]-1;

53
src/LinAlg/SAM.h
View File

@@ -23,7 +23,7 @@ class SystemVector;
/*!
\brief This class contains data and functions for the assembly of FE matrices.
\details The names and meanings of the data members of this class are
\details The names and meanings of (most of) the data members of this class
are adopted from Kolbein Bell's pionering work on the field.
See his reports on the SAM library for a thorough elaboration.
@@ -33,6 +33,10 @@ class SystemVector;
class SAM
{
protected:
typedef std::vector<int> IntVec; //!< General integer vector
typedef std::set<int> IntSet; //!< General integer set
public:
//! \brief The constructor initializes an empty object.
SAM();
@@ -51,12 +55,10 @@ public:
//! \brief Returns the number of equations (free DOFs) in the model.
virtual int getNoEquations() const { return neq; }
//! \brief Returns max number of dof couplings in the model.
//! \brief Returns max number of DOF couplings in the model.
int getMaxDofCouplings() const;
typedef std::vector<int> IntVec; //!< General integer vector
//! \brief Computes number of couplings for each dof in the system matrix.
//! \brief Computes number of couplings for each DOF in the system matrix.
//! \param[out] nnz Number of couplings (non-zeroes) for each DOF
//! \return \e false if number of couplings is not computed, otherwise \e true
bool getNoDofCouplings(IntVec& nnz) const;
@@ -72,8 +74,6 @@ public:
bool getDofCouplings(IntVec& irow, IntVec& jcol) const;
private:
typedef std::set<int> IntSet; //!< General integer set
//! \brief Finds the set of free DOFs coupled to each free DOF.
bool getDofCouplings(std::vector<IntSet>& dofc) const;
@@ -192,24 +192,38 @@ public:
bool expandVector(const Vector& solVec, Vector& displ) const;
//! \brief Computes the dot-product of two vectors of length NDOF.
virtual real dot(const Vector& x, const Vector& y) const;
//! \param[in] x The first vector of the dot-product
//! \param[in] y The second vector of the dot-product
//! \param[in] dofType Only consider nodes of this type (for mixed methods)
virtual real dot(const Vector& x, const Vector& y, char dofType = 'D') const;
//! \brief Computes the l2-norm of a vector of length NDOF.
real norm2(const Vector& x) const { return sqrt(this->dot(x,x)); }
//! \param[in] x The vector to compute the norm of
//! \param[in] dofType Only consider nodes of this type (for mixed methods)
real norm2(const Vector& x, char dofType = 'D') const
{ return sqrt(this->dot(x,x,dofType)); }
//! \brief Computes the L2-norm of a vector of length NDOF.
virtual real normL2(const Vector& x) const;
//! \param[in] x The vector to compute the norm of
//! \param[in] dofType Only consider nodes of this type (for mixed methods)
virtual real normL2(const Vector& x, char dofType = 'D') const;
//! \brief Computes the L_infinity-norm of a vector of length NDOF.
//! \param[in] x The vector to compute the norm of
//! \param comp Local nodal DOF on input, index of the largest value on output
virtual real normInf(const Vector& x, size_t& comp) const;
//! \brief Computes energy norm contributions from nodal reaction forces.
//! \param[in] dofType Only consider nodes of this type (for mixed methods)
virtual real normInf(const Vector& x, size_t& comp, char dofType = 'D') const;
//! \brief Computes the energy norm contributions from nodal reaction forces.
//! \param[in] u The (incremental) nodal displacement vector
//! \param[in] rf Compressed reaction force vector
//! \param[in] rf Compressed reaction force vector for the entire model
virtual real normReact(const Vector& u, const Vector& rf) const;
//! \brief Returns the total reaction force in the given coordinate direction.
//! \param[in] dir 1-based coordinate direction index
//! \param[in] rf Compressed reaction force vector for the entire model
virtual real getReaction(int dir, const Vector& rf) const;
//! \brief Returns a vector of reaction forces for a given node.
//! \param[in] inod Identifier for the node to get the reaction forces for
//! \param[in] rf Compressed reaction force vector for the entire model
//! \param[out] nrf Nodal reaction forces
//! \return \e true of the specified node has reaction forces
//! \return \e true if the specified node has reaction forces
bool getNodalReactions(int inod, const Vector& rf, Vector& nrf) const;
protected:
@@ -233,16 +247,23 @@ protected:
//! \return second = local DOF number in the range [1,NNDOF]
std::pair<int,int> getNodeAndLocalDof(int idof) const;
private:
int mpar[50]; //!< Matrix of parameters
protected:
// The following parameters are pointers to specific locations in MPAR
int& nnod; //!< Number of nodes
int& nel; //!< Number of elements
int& ndof; //!< Number of DOFs
int& nspdof; //!< Number of specified DOFs
int& nceq; //!< Number of constraint equations
int& neq; //!< Number of system equations
int& nmmnpc; //!< Number of elements in MMNPC
int& nmmceq; //!< Number of elements in MMCEQ
// The standard SAM arrays (see K. Bell's reports for detailed explanation)
// The standard SAM arrays (see K. Bell's reports for detailed explanation).
// We are using plane C-pointers for these items such that they more easily
// can be passed directly as arguments to FORTRAN subroutines.
int* mpmnpc; //!< Matrix of pointers to MNPCs in MMNPC
int* mmnpc; //!< Matrix of matrices of nodal point correspondances
int* madof; //!< Matrix of accumulated DOFs
@@ -253,7 +274,7 @@ protected:
int* minex; //!< Matrix of internal to external node numbers
int* meqn; //!< Matrix of equation numbers
int mpar[50]; //!< Matrix of parameters
std::vector<char> nodeType; //!< DOF classification when using mixed methods
friend class DenseMatrix;
friend class SPRMatrix;

10
src/SIM/NonLinSIM.C
View File

@@ -336,6 +336,9 @@ bool NonLinSIM::solutionNorms (const TimeDomain& time, const char* compName,
double dMax[nsd];
double normL2 = model->solutionNorms(solution.front(),dMax,iMax);
RealArray RF;
bool haveRFs = model->getCurrentReactions(RF,solution.front());
if (energyNorm)
if (!model->solutionNorms(time,solution,eNorm,gNorm))
gNorm.clear();
@@ -347,6 +350,13 @@ bool NonLinSIM::solutionNorms (const TimeDomain& time, const char* compName,
for (int d = 0; d < nsd; d++, D++)
std::cout <<"\n Max "<< D <<'-'<< compName
<<" : "<< dMax[d] <<" node "<< iMax[d];
if (haveRFs)
{
std::cout <<"\n Total reaction forces: Sum(R) =";
for (size_t i = 1; i < RF.size(); i++)
std::cout <<" "<< RF[i];
std::cout <<"\n "<< compName <<"*reactions: (R,u) = "<< RF.front();
}
if (gNorm.size() > 0)
{
std::cout <<"\n Energy norm: |u^h| = a(u^h,u^h)^0.5 : "<< gNorm(1);

30
src/SIM/SIMbase.C
View File

@@ -665,9 +665,9 @@ bool SIMbase::solveSystem (Vector& solution, int printSol,
void SIMbase::iterationNorms (const Vector& u, const Vector& r,
double& eNorm, double& rNorm, double& dNorm) const
{
eNorm = mySam->dot(r,u);
rNorm = mySam->norm2(r);
dNorm = mySam->norm2(u);
eNorm = mySam->dot(r,u,'A');
rNorm = mySam->norm2(r,'D');
dNorm = mySam->norm2(u,'D');
}
@@ -812,6 +812,30 @@ bool SIMbase::solutionNorms (const TimeDomain& time, const Vectors& psol,
}
/*!
The content of the output array \a RF is as follows:
\f[
RF[0] = \sum_{i=n}^{\rm nnod} \sum_{i=1}^{\rm nsd} f_n^i u_n^i
\quad\quad\mbox{(energy)}\f]
\f[
RF[i] = \sum_{n=1}^{\rm nnod} f_n^i \quad\forall\quad i=1,\ldots,{\rm nsd}
\f]
*/
bool SIMbase::getCurrentReactions (RealArray& RF, const Vector& psol) const
{
const Vector* reactionForces = myEqSys->getReactions();
if (!reactionForces) return false;
RF.resize(1+this->getNoSpaceDim());
RF.front() = 2.0*mySam->normReact(psol,*reactionForces);
for (size_t dir = 1; dir < RF.size(); dir++)
RF[dir] = mySam->getReaction(dir,*reactionForces);
return true;
}
bool SIMbase::systemModes (std::vector<Mode>& solution,
int nev, int ncv, int iop, double shift,
size_t iA, size_t iB)

11
src/SIM/SIMbase.h
View File

@@ -197,7 +197,7 @@ public:
//! \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 vectors
//! \param[in] psol Primary solution vectors
//! \param[out] eNorm Element-wise norm quantities
//! \param[out] gNorm Global norm quantities
virtual bool solutionNorms(const TimeDomain& time, const Vectors& psol,
@@ -206,7 +206,7 @@ public:
//! \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 vectors
//! \param[in] psol Primary solution vectors
//! \param[out] eNorm Element-wise norm quantities
//! \param[out] gNorm Global norm quantities
//!
@@ -214,6 +214,11 @@ public:
virtual bool solutionNorms(const Vectors& psol, Matrix& eNorm, Vector& gNorm)
{ return this->solutionNorms(TimeDomain(),psol,eNorm,gNorm); }
//! \brief Computes the total reaction forces in the model.
//! \param[out] RF Reaction force in each spatial direction + energy
//! \param[in] psol Primary solution vector
bool getCurrentReactions(RealArray& RF, const Vector& psol) const;
//! \brief Performs a generalized eigenvalue analysis of the assembled system.
//! \param[in] iop Which eigensolver method to use
//! \param[in] nev Number of eigenvalues/vector (see ARPack documentation)
@@ -272,7 +277,7 @@ public:
//! \brief Writes solution fields for a given load/time step to the VTF-file.
//! \details If an analytical solution is provided, the exact stress fields
//! are written to the VTF-file as well.
//! \param[in] psol Global primary solution vector
//! \param[in] psol Primary solution vector
//! \param[in] nViz Number of visualization points over each knot-span
//! \param[in] iStep Load/time step identifier
//! \param nBlock Running result block counter