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

Added: New class for generalized-alpha time integration.

Browse Source 
Changed: Calculate internal forces by integrating B^t*sigma also for
linear problems, instead of relying on K*d in NewmarkMats::getRHSVector.
This commit is contained in:
Knut Morten Okstad
2014-08-21 16:10:28 +02:00
parent 974557f37f
commit 0d2e30cde8
4 changed files with 291 additions and 20 deletions
Download Patch File Download Diff File Expand all files Collapse all files

30
src/ASM/NewmarkMats.C
View File

@@ -15,15 +15,26 @@
#include "NewmarkMats.h"
NewmarkMats::NewmarkMats (double a1, double a2, double b, double c)
NewmarkMats::NewmarkMats (double a1, double a2, double b, double c,
bool generalizedAlpha) : isPredictor(true), h(0.0)
{
alpha1 = a1;
alpha2 = a2;
if (generalizedAlpha)
{
alpha_m = b;
alpha_f = c;
double alpha = alpha_f - alpha_m;
beta = 0.25*(1.0-alpha)*(1.0-alpha);
gamma = 0.5 - alpha;
}
else
{
alpha_m = alpha_f = 1.0;
beta = b;
gamma = c;
isPredictor = true;
h = 0.0;
}
}
@@ -32,10 +43,8 @@ const Matrix& NewmarkMats::getNewtonMatrix () const
Matrix& N = const_cast<Matrix&>(A.front());
N = A[1];
if (alpha1 > 0.0)
N.multiply(1.0 + alpha1*gamma*h); // [N] = (1+alpha1*gamma*h)*[M]
N.add(A[2],(alpha2*gamma + beta*h)*h); // [N] += (alpha2*gamma+beta*h)*h*[K]
N.multiply(alpha_m + alpha_f*alpha1*gamma*h);
N.add(A[2],alpha_f*(alpha2*gamma + beta*h)*h);
#if SP_DEBUG > 2
std::cout <<"\nElement mass matrix"<< A[1];
@@ -56,13 +65,12 @@ const Vector& NewmarkMats::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 <<"\nf_ext"<< dF;
std::cout <<"\nf_ext - f_s"<< dF;
std::cout <<"f_i = M*a"<< A[1]*vec[ia];
if (alpha1 > 0.0)
std::cout <<"f_d1/alpha1 = M*v (alpha1="<< alpha1 <<")"<< A[1]*vec[iv];
if (alpha2 > 0.0)
std::cout <<"f_d2/alpha2 = K*v (alpha2="<< alpha2 <<")"<< A[2]*vec[iv];
std::cout <<"f_s = K*d"<< A[2]*vec.front();
#endif
dF.add(A[1]*vec[ia],-1.0); // dF = Fext - M*a
@@ -72,8 +80,6 @@ const Vector& NewmarkMats::getRHSVector () const
if (alpha2 > 0.0)
dF.add(A[2]*vec[iv],-alpha2); // dF -= alpha2*K*v
dF.add(A[2]*vec.front(),-1.0); // dF -= K*d
}
#if SP_DEBUG > 2

26
src/ASM/NewmarkMats.h
View File

@@ -26,11 +26,21 @@ class NewmarkMats : public ElmMats
{
public:
//! \brief The constructor initializes the time integration parameters.
NewmarkMats(double a1 = 0.0, double a2 = 0.0, double b = 0.0, double c = 0.0);
//! \param[in] a1 Mass-proportional damping coefficient
//! \param[in] a2 Stiffness-proportional damping coefficient
//! \param[in] b Time integration parameter
//! \param[in] c Time integration parameter
//! \param[in] generalizedAlpha If \e true, iterpret \a b and \a c as the
//! generalized alpha parameters \a alpha_m and \a alpha_f, respectively,
//! otherwise as \a beta and \a gamma
NewmarkMats(double a1 = 0.0, double a2 = 0.0, double b = 0.0, double c = 0.0,
bool generalizedAlpha = false);
//! \brief Empty destructor.
virtual ~NewmarkMats() {}
//! \brief Updates the time step size and the \a isPredictor flag.
//! \param[in] dt New time step size
//! \param[in] iter Newton-Raphson iteration counter
void setStepSize(double dt, int iter) { h = dt; isPredictor = iter == 0; }
//! \brief Returns the element-level Newton matrix.
@@ -39,13 +49,17 @@ public:
virtual const Vector& getRHSVector() const;
protected:
double alpha1; //!< Mass-proportional damping
double alpha2; //!< Stiffness-proportional damping
double beta; //!< Time integration parameter
double gamma; //!< Time integration parameter
bool isPredictor; //!< If \e true, we are in the predictor step
double h; //!< Time step size
bool isPredictor; //!< If \e true, we are in the predictor step
double alpha1; //!< Mass-proportional damping coefficient
double alpha2; //!< Stiffness-proportional damping coefficient
double beta; //!< Newmark time integration parameter
double gamma; //!< Newmark time integration parameter
private:
double alpha_m; //!< Generalized-alpha parameter
double alpha_f; //!< Generalized-alpha parameter
};
#endif

189
src/SIM/GenAlphaSIM.C Normal file
View File

@@ -0,0 +1,189 @@
// $Id$
//==============================================================================
//!
//! \file GenAlphaSIM.C
//!
//! \date Aug 21 2014
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Generalized-alpha driver for isogeometric dynamic FEM simulators.
//!
//==============================================================================
#include "GenAlphaSIM.h"
#include "SIMoutput.h"
#include "TimeStep.h"
#include "IFEM.h"
#include "tinyxml.h"
GenAlphaSIM::GenAlphaSIM (SIMbase& s) : NewmarkSIM(s), prevSol(3), tempSol(3)
{
// Default Newmark parameters (alpha = -0.1)
alpha_m = 1.0;
alpha_f = 0.9;
beta = 0.3025;
gamma = 0.6;
}
bool GenAlphaSIM::parse (const TiXmlElement* elem)
{
bool ok = this->NewmarkSIM::parse(elem);
if (!strcasecmp(elem->Value(),"newmarksolver"))
{
double alpha = -0.1;
const char* attr = elem->Attribute("alpha");
if (attr)
{
alpha = atof(attr);
alpha_f = alpha + 1.0;
alpha_m = 1.0;
}
else if ((attr = elem->Attribute("rho_inf")))
{
double rho = atof(attr);
alpha_f = 1.0/(1.0+rho);
alpha_m = (2.0-rho)/(1.0+rho);
alpha = alpha_f - alpha_m;
}
beta = 0.25*(1.0-alpha)*(1.0-alpha);
gamma = 0.5 - alpha;
}
return ok;
}
void GenAlphaSIM::initPrm ()
{
model.setIntegrationPrm(0,alpha1);
model.setIntegrationPrm(1,fabs(alpha2));
model.setIntegrationPrm(2,alpha_m);
model.setIntegrationPrm(3,alpha_f);
model.setIntegrationPrm(4,2.0);
}
bool GenAlphaSIM::initSol (size_t nSol)
{
this->MultiStepSIM::initSol(nSol);
size_t nDOFs = model.getNoDOFs();
for (size_t i = 0; i < 3; i++)
{
prevSol[i].resize(nDOFs,true);
tempSol[i].resize(nDOFs,true);
}
return true;
}
bool GenAlphaSIM::advanceStep (TimeStep& param, bool updateTime)
{
// Update displacement solutions between time steps
for (int n = solution.size()-3; n > 0; n--)
std::copy(solution[n-1].begin(),solution[n-1].end(),solution[n].begin());
// Update the previous solution
std::copy(solution.front().begin(),solution.front().end(),prevSol[0].begin());
if (solution.size() > 2)
{
size_t iA = solution.size() - 1;
size_t iV = solution.size() - 2;
std::copy(solution[iV].begin(),solution[iV].end(),prevSol[1].begin());
std::copy(solution[iA].begin(),solution[iA].end(),prevSol[2].begin());
}
return this->NewmarkSIM::advanceStep(param,updateTime);
}
bool GenAlphaSIM::predictStep (TimeStep& param)
{
if (!this->NewmarkSIM::predictStep(param))
return false;
size_t ipD = 0;
size_t ipA = solution.size() - 1;
size_t ipV = solution.size() - 2;
tempSol[0] = solution[ipD];
tempSol[1] = solution[ipV];
tempSol[2] = solution[ipA];
// Compute the intermediate solution
// {U}_(n+alpha) = (1-alpha)*{U}_n + alpha*{U}_(n+1)
solution[ipD].relax(alpha_f,prevSol[0]);
solution[ipV].relax(alpha_f,prevSol[1]);
solution[ipA].relax(alpha_m,prevSol[2]);
return true;
}
bool GenAlphaSIM::correctStep (TimeStep& param, bool converged)
{
if (solution.size() < 3)
{
std::cerr <<" *** GenAlphaSIM::correctStep: Too few solution vectors "
<< solution.size() << std::endl;
return false;
}
#ifdef SP_DEBUG
std::cout <<"\nGenAlphaSIM::correctStep(converged="
<< std::boolalpha << converged <<")";
#endif
size_t ipD = 0;
size_t ipA = solution.size() - 1;
size_t ipV = solution.size() - 2;
// Update current displacement, velocity and acceleration solutions
tempSol[0].add(linsol,beta*param.time.dt*param.time.dt);
tempSol[1].add(linsol,gamma*param.time.dt);
tempSol[2].add(linsol,1.0);
if (converged)
{
// Converged solution
solution[ipD] = tempSol[0];
solution[ipV] = tempSol[1];
solution[ipA] = tempSol[2];
}
else
{
// Compute new intermediate solution
// {U}_(n+alpha) = (1-alpha)*{U}_n + alpha*{U}_(n+1)
solution[ipD].relax(alpha_f,prevSol[0],tempSol[0]);
solution[ipV].relax(alpha_f,prevSol[1],tempSol[1]);
solution[ipA].relax(alpha_m,prevSol[2],tempSol[2]);
}
#if SP_DEBUG > 1
std::cout <<"\nCorrected displacement:"<< solution[ipD]
<<"Corrected velocity:"<< solution[ipV]
<<"Corrected acceleration:"<< solution[ipA];
#elif defined(SP_DEBUG)
if (converged)
std::cout <<"\nConverged displacement:"<< solution[ipD].max()
<<"\nConverged velocity:"<< solution[ipV].max()
<<"\nConverged acceleration:"<< solution[ipA].max() << std::endl;
#endif
model.updateRotations(linsol,alpha_f); //TODO look at this
return model.updateConfiguration(tempSol[0]);
}
void GenAlphaSIM::setSolution (const Vector& newSol, int idx)
{
solution[idx] = newSol;
if (idx == 0)
tempSol[0] = newSol;
else if (solution.size() > 2)
tempSol[solution.size()-3+idx] = newSol;
}

62
src/SIM/GenAlphaSIM.h Normal file
View File

@@ -0,0 +1,62 @@
// $Id$
//==============================================================================
//!
//! \file GenAlphaSIM.h
//!
//! \date Aug 21 2014
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Generalized-alpha driver for isogeometric dynamic FEM simulators.
//!
//==============================================================================
#ifndef _GEN_ALPHA_SIM_H
#define _GEN_ALPHA_SIM_H
#include "NewmarkSIM.h"
/*!
\brief Generalized-alpha driver for dynamic isogeometric FEM simulators.
*/
class GenAlphaSIM : public NewmarkSIM
{
public:
//! \brief The constructor initializes default solution parameters.
GenAlphaSIM(SIMbase& s);
//! \brief Empty destructor.
virtual ~GenAlphaSIM() {}
//! \brief Parses a data section from an XML document.
virtual bool parse(const TiXmlElement* elem);
//! \brief Initializes time integration parameters for the integrand.
virtual void initPrm();
//! \brief Initializes primary solution vectors.
//! \param[in] nSol Number of consequtive solutions stored
virtual bool initSol(size_t nSol = 3);
//! \brief Advances the time step one step forward.
//! \param param Time stepping parameters
//! \param[in] updateTime If \e false, the time parameters are not incremented
virtual bool advanceStep(TimeStep& param, bool updateTime = true);
//! \brief Modifies the current solution vector (used by sub-iterations only).
virtual void setSolution(const Vector& newSol, int idx);
protected:
//! \brief Calculates predicted velocities and accelerations.
virtual bool predictStep(TimeStep& param);
//! \brief Updates configuration variables (solution vector) in an iteration.
virtual bool correctStep(TimeStep& param, bool converged);
private:
double alpha_m; //!< Generalized-alpha parameter
double alpha_f; //!< Generalized-alpha parameter
Vectors prevSol; //!< Solution of last converged time step
Vectors tempSol; //!< Intermediate solution
};
#endif