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

Changed: Unified the class GenAlphaMats into the class HHTMats.

Browse Source
This commit is contained in:
Knut Morten Okstad
2015-11-26 13:49:55 +01:00
parent 797b05ac3e
commit 974557f37f
4 changed files with 17 additions and 150 deletions
Download Patch File Download Diff File Expand all files Collapse all files

102
src/ASM/GenAlphaMats.C
View File

@@ -1,102 +0,0 @@
// $Id$
//==============================================================================
//!
//! \file GenAlphaMats.C
//!
//! \date Jul 04 2013
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Representation of the element matrices for a dynamic FEM problem.
//!
//==============================================================================
#include "GenAlphaMats.h"
GenAlphaMats::GenAlphaMats (double alpha, double a, double b) : NewmarkMats(a,b)
{
beta = 0.25*(1.0-alpha)*(1.0-alpha);
gamma = 0.5 - alpha;
}
const Matrix& GenAlphaMats::getNewtonMatrix () const
{
Matrix& N = const_cast<Matrix&>(A.front());
double alphaPlus1 = 1.5 - gamma; // = 1.0 + alpha
N = A[2];
N.multiply(alphaPlus1*(1.0 + alpha2*gamma/(beta*h)));
if (A.size() > 3)
N.add(A[3],alphaPlus1);
N.add(A[1],(alphaPlus1*gamma*alpha1 + 1.0/h)/(beta*h));
#if SP_DEBUG > 2
std::cout <<"\nElement mass matrix"<< A[1];
if (A.size() > 3)
std::cout <<"Material stiffness matrix"<< A[2]
<<"Geometric stiffness matrix"<< A[3];
else
std::cout <<"Tangent stiffness matrix"<< A[2];
std::cout <<"Resulting Newton matrix"<< A[0];
#endif
return A.front();
}
const Vector& GenAlphaMats::getRHSVector () const
{
int ia = vec.size() - 1; // index to element acceleration vector (a)
int iv = vec.size() - 2; // index to element velocity vector (v)
#if SP_DEBUG > 2
std::cout <<"GenAlphaMats::getRHSVector "
<< (isPredictor ? "predictor" : "corrector") << std::endl;
#endif
if (A.size() > 2 && b.size() > 1 && vec.size() > 2)
{
Vector& Fi = const_cast<Vector&>(b[1]);
// Find the actual inertia force from the dynamic equilibrium equation
// (Fi = Fext - Fs - Fd) and store it in the second RHS vector Fi=b[1]
Fi = b.front(); // Fi = Fext - Fs
if (alpha1 > 0.0 && iv > 0)
Fi.add(A[1]*vec[iv],-alpha1); // Fi -= alpha1*M*v
if (alpha2 > 0.0 && iv > 0)
Fi.add(A[2]*vec[iv],-alpha2); // Fi -= alpha2*K*v
#if SP_DEBUG > 2
std::cout <<"Element inertia vector"<< Fi;
#endif
}
if (A.size() > 2 && vec.size() > 2)
{
Vector& RHS = const_cast<Vector&>(b.front());
#if SP_DEBUG > 2
std::cout <<"\nElement velocity vector"<< vec[iv];
std::cout <<"Element acceleration vector"<< vec[ia];
std::cout <<"S_ext - S_int"<< b.front();
std::cout <<"S_i"<< A[1]*vec[ia];
#endif
// Find the right-hand-side vector
double alphaPlus1 = 1.5 - gamma;
RHS *= alphaPlus1;
if (alpha1 > 0.0)
RHS.add(A[1]*vec[iv], alpha1*(isPredictor ? alphaPlus1 : -alphaPlus1));
if (alpha2 > 0.0)
RHS.add(A[2]*vec[iv], alpha2*(isPredictor ? alphaPlus1 : -alphaPlus1));
RHS.add(A[1]*vec[ia], isPredictor ? 1.0 : -1.0);
}
#if SP_DEBUG > 2
std::cout <<"\nElement right-hand-side vector"<< b.front();
#endif
return b.front();
}

39
src/ASM/GenAlphaMats.h
View File

@@ -1,39 +0,0 @@
// $Id$
//==============================================================================
//!
//! \file GenAlphaMats.h
//!
//! \date Jul 04 2013
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Representation of the element matrices for a dynamic FEM problem.
//!
//==============================================================================
#ifndef _GEN_ALPHA_MATS_H
#define _GEN_ALPHA_MATS_H
#include "NewmarkMats.h"
/*!
\brief Class representing the element matrices for a dynamic FEM problem
based on generalized alpha time integration.
*/
class GenAlphaMats : public NewmarkMats
{
public:
//! \brief The constructor initializes the time integration parameters.
GenAlphaMats(double alpha, double a, double b);
//! \brief Empty destructor.
virtual ~GenAlphaMats() {}
//! \brief Returns the element-level Newton matrix.
virtual const Matrix& getNewtonMatrix() const;
//! \brief Returns the element-level right-hand-side vector.
virtual const Vector& getRHSVector() const;
};
#endif

21
src/ASM/HHTMats.C
View File

@@ -14,10 +14,11 @@
#include "HHTMats.h" #include "HHTMats.h"
HHTMats::HHTMats (double alpha, double a, double b) : NewmarkMats(a,b) HHTMats::HHTMats (double alpha, double a, double b, bool old) : NewmarkMats(a,b)
{ {
beta = 0.25*(1.0-alpha)*(1.0-alpha); beta = 0.25*(1.0-alpha)*(1.0-alpha);
gamma = 0.5 - alpha; gamma = 0.5 - alpha;
oldHHT = old;
} }
@@ -37,7 +38,7 @@ const Matrix& HHTMats::getNewtonMatrix () const
N.add(A[1],(alphaPlus1*alpha1*gamma + 1.0/h)/(beta*h)); N.add(A[1],(alphaPlus1*alpha1*gamma + 1.0/h)/(beta*h));
#if SP_DEBUG > 2 #if SP_DEBUG > 2
std::cout <<"\nDynElmMats::getNewtonMatrix"; std::cout <<"\nHHTMats::getNewtonMatrix";
std::cout <<"\nElement mass matrix"<< A[1]; std::cout <<"\nElement mass matrix"<< A[1];
if (A.size() > 3) if (A.size() > 3)
std::cout <<"Material stiffness matrix"<< A[2] std::cout <<"Material stiffness matrix"<< A[2]
@@ -63,7 +64,7 @@ const Vector& HHTMats::getRHSVector () const
int ia = vec.size() - 1; // index to element acceleration vector (a) int ia = vec.size() - 1; // index to element acceleration vector (a)
int iv = vec.size() - 2; // index to element velocity vector (v) int iv = vec.size() - 2; // index to element velocity vector (v)
if (isPredictor && b.size() > 1) if ((isPredictor || oldHHT) && b.size() > 1)
{ {
Vector& Fia = const_cast<Vector&>(b[1]); Vector& Fia = const_cast<Vector&>(b[1]);
@@ -76,11 +77,12 @@ const Vector& HHTMats::getRHSVector () const
Fia.add(A[1]*vec[iv],-alpha1); // Fia -= alpha1*M*v Fia.add(A[1]*vec[iv],-alpha1); // Fia -= alpha1*M*v
if (alpha2 > 0.0 && iv > 0) if (alpha2 > 0.0 && iv > 0)
Fia.add(A[2]*vec[iv],-alpha2); // Fia -= alpha2*K*v Fia.add(A[2]*vec[iv],-alpha2); // Fia -= alpha2*K*v
#if SP_DEBUG > 2
std::cout <<"Element inertia vector"<< Fia;
#endif
} }
#if SP_DEBUG > 2 #if SP_DEBUG > 2
if (isPredictor && b.size() > 1)
std::cout <<"Element inertia vector"<< b[1];
std::cout <<"Element velocity vector"<< vec[iv]; std::cout <<"Element velocity vector"<< vec[iv];
std::cout <<"Element acceleration vector"<< vec[ia]; std::cout <<"Element acceleration vector"<< vec[ia];
if (b.size() > 2) if (b.size() > 2)
@@ -97,7 +99,9 @@ const Vector& HHTMats::getRHSVector () const
// the predictor step force vector after the element assembly. // the predictor step force vector after the element assembly.
Vector& RHS = const_cast<Vector&>(b.front()); Vector& RHS = const_cast<Vector&>(b.front());
double alphaPlus1 = 1.5 - gamma; double alphaPlus1 = 1.5 - gamma;
if (isPredictor) if (oldHHT)
RHS *= alphaPlus1; // RHS = -(1+alphaH)*S_int
else if (isPredictor)
{ {
int pa = iv - 1; // index to predicted element acceleration vector (A) int pa = iv - 1; // index to predicted element acceleration vector (A)
A[1].multiply(vec[pa]-vec[ia],RHS); // RHS = M*(A-a) A[1].multiply(vec[pa]-vec[ia],RHS); // RHS = M*(A-a)
@@ -121,9 +125,10 @@ const Vector& HHTMats::getRHSVector () const
RHS.add(A[1]*vec[iv],alphaPlus1*alpha1); // RHS -= (1+alphaH)*alpha1*M*v RHS.add(A[1]*vec[iv],alphaPlus1*alpha1); // RHS -= (1+alphaH)*alpha1*M*v
if (alpha2 > 0.0) if (alpha2 > 0.0)
RHS.add(A[2]*vec[iv],alphaPlus1*alpha2); // RHS -= (1+alphaH)*alpha2*K*v RHS.add(A[2]*vec[iv],alphaPlus1*alpha2); // RHS -= (1+alphaH)*alpha2*K*v
if (oldHHT)
RHS.add(A[1]*vec[ia], isPredictor ? 1.0 : -1.0); // RHS (+/-)= M*a
#if SP_DEBUG > 2 #if SP_DEBUG > 2
std::cout <<"Element right-hand-side vector"<< b.front(); std::cout <<"\nElement right-hand-side vector"<< b.front();
#endif #endif
return b.front(); return b.front();

5
src/ASM/HHTMats.h
View File

@@ -26,7 +26,7 @@ class HHTMats : public NewmarkMats
{ {
public: public:
//! \brief The constructor initializes the time integration parameters. //! \brief The constructor initializes the time integration parameters.
HHTMats(double alpha, double a, double b); HHTMats(double alpha, double a, double b, bool old = false);
//! \brief Empty destructor. //! \brief Empty destructor.
virtual ~HHTMats() {} virtual ~HHTMats() {}
@@ -34,6 +34,9 @@ public:
virtual const Matrix& getNewtonMatrix() const; virtual const Matrix& getNewtonMatrix() const;
//! \brief Returns the element-level right-hand-side vector. //! \brief Returns the element-level right-hand-side vector.
virtual const Vector& getRHSVector() const; virtual const Vector& getRHSVector() const;
private:
bool oldHHT; //!< If \e true, used toghether with NewmarkNLSIM
}; };
#endif #endif