#include <thirdorderstepsol.h>

Inheritance diagram for ThirdOrderIntegrator:

Collaboration diagram for ThirdOrderIntegrator:

Public Member Functions
	ThirdOrderIntegrator (const doublereal dT, const doublereal dSolutionTol, const integer iMaxIt, const bool bmod_res_test)

virtual	~ThirdOrderIntegrator (void)

virtual void	Residual (VectorHandler *pRes) const

virtual void	Jacobian (MatrixHandler *pJac) const

virtual void	Update (const VectorHandler *pSol) const

virtual doublereal	Advance (Solver pS, const doublereal TStep, const doublereal dAlph, const StepChange StType, std::deque< MyVectorHandler > &qX, std::deque< MyVectorHandler * > &qXPrime, MyVectorHandler const pX, MyVectorHandler const pXPrime, integer &EffIter, doublereal &Err, doublereal &SolErr)

Public Member Functions inherited from ImplicitStepIntegrator
	ImplicitStepIntegrator (const integer MaxIt, const doublereal dT, const doublereal dSolutionTol, const integer stp, const integer sts, const bool bmod_res_test)

virtual	~ImplicitStepIntegrator (void)

virtual void	EvalProd (doublereal Tau, const VectorHandler &f0, const VectorHandler &w, VectorHandler &z) const

virtual doublereal	TestScale (const NonlinearSolverTest *pTest, doublereal &dCoef) const

Public Member Functions inherited from StepIntegrator
	StepIntegrator (const integer MaxIt, const doublereal dT, const doublereal dSolutionTol, const integer stp, const integer sts)

virtual	~StepIntegrator (void)

void	SetDataManager (DataManager *pDatMan)

virtual integer	GetIntegratorNumPreviousStates (void) const

virtual integer	GetIntegratorNumUnknownStates (void) const

virtual integer	GetIntegratorMaxIters (void) const

virtual doublereal	GetIntegratorDTol (void) const

virtual doublereal	GetIntegratorDSolTol (void) const

virtual void	OutputTypes (const bool fpred)

virtual void	SetDriveHandler (const DriveHandler *pDH)

Public Member Functions inherited from NonlinearProblem
virtual	~NonlinearProblem (void)

Protected Member Functions
void	PredictDof_for_AfterPredict (const int DCount, const DofOrder::Order Order, const VectorHandler *const pSol=0) const

void	RealPredictDof (const int DCount, const DofOrder::Order Order, const VectorHandler *const pSol=0) const

void	UpdateDof (const int DCount, const DofOrder::Order Order, const VectorHandler *const pSol=0) const

virtual void	SetCoef (doublereal dT, doublereal dAlpha, enum StepChange NewStep)=0

virtual void	Predict (void)

Protected Member Functions inherited from StepIntegrator
template<class T >
void	UpdateLoop (const T const t, void(T::pUpd)(const int DCount, const DofOrder::Order Order, const VectorHandler const pSol) const, const VectorHandler const pSol=0) const

Protected Attributes
VectorHandler *	pXPrev

VectorHandler *	pXPrimePrev

doublereal	dT

doublereal	rho

doublereal	theta

doublereal	w0

doublereal	w1

doublereal	w2

doublereal	jx11

doublereal	jx12

doublereal	jx21

doublereal	jx22

doublereal	m0

doublereal	m1

doublereal	n0

doublereal	n1

bool	bAdvanceCalledFirstTime

MyVectorHandler	Restmp

std::vector< bool >	EqIsAlgebraic

std::vector< bool >	EqIsDifferential

SpMapMatrixHandler	Jacxi_xp

SpMapMatrixHandler	Jacxi_x

SpMapMatrixHandler	Jac_xp

SpMapMatrixHandler	Jac_x

MatrixHandler *	pJacxi_xp

MatrixHandler *	pJacxi_x

MatrixHandler *	pJac_xp

MatrixHandler *	pJac_x

Protected Attributes inherited from ImplicitStepIntegrator
VectorHandler *	pXCurr

VectorHandler *	pXPrimeCurr

bool	bModResTest

Protected Attributes inherited from StepIntegrator
DataManager *	pDM

const DataManager::DofVecType *	pDofs

bool	outputPred

integer	MaxIters

doublereal	dTol

doublereal	dSolTol

integer	steps

integer	unkstates

Additional Inherited Members
Public Types inherited from StepIntegrator
enum	{ DIFFERENTIAL = 0, ALGEBRAIC = 1 }

enum	StepChange { NEWSTEP, REPEATSTEP }

Detailed Description

Definition at line 40 of file thirdorderstepsol.h.

Constructor & Destructor Documentation

ThirdOrderIntegrator::ThirdOrderIntegrator	(	const doublereal	dT,
		const doublereal	dSolutionTol,
		const integer	iMaxIt,
		const bool	bmod_res_test
	)

Definition at line 52 of file thirdorderstepsol.cc.

References Jac_x, Jac_xp, Jacxi_x, Jacxi_xp, pJac_x, pJac_xp, pJacxi_x, and pJacxi_xp.

 : ImplicitStepIntegrator(iMaxIt, dT, dSolutionTol, 1, 2, bmod_res_test),
 pXPrev(0),
 pXPrimePrev(0),
 bAdvanceCalledFirstTime(true)
 {
         pJacxi_xp = &Jacxi_xp;
         pJacxi_x  = &Jacxi_x;
         pJac_xp   = &Jac_xp;
         pJac_x    = &Jac_x;
 };

ThirdOrderIntegrator::~ThirdOrderIntegrator ( void )

virtual

Definition at line 67 of file thirdorderstepsol.cc.

References NO_OP.

                                            {
         NO_OP;
 };

Member Function Documentation

doublereal ThirdOrderIntegrator::Advance	(	Solver *	pS,
		const doublereal	TStep,
		const doublereal	dAlph,
		const StepChange	StType,
		std::deque< MyVectorHandler * > &	qX,
		std::deque< MyVectorHandler * > &	qXPrime,
		MyVectorHandler *const	pX,
		MyVectorHandler *const	pXPrime,
		integer &	EffIter,
		doublereal &	Err,
		doublereal &	SolErr
	)

virtual

Implements StepIntegrator.

Definition at line 72 of file thirdorderstepsol.cc.

References DofOrder::ALGEBRAIC, ASSERT, bAdvanceCalledFirstTime, StepIntegrator::dSolTol, StepIntegrator::dTol, EqIsAlgebraic, EqIsDifferential, DataManager::iGetNumDofs(), Jac_x, Jac_xp, Jacxi_x, Jacxi_xp, DataManager::LinkToSolution(), StepIntegrator::MaxIters, StepIntegrator::pDM, StepIntegrator::pDofs, Solver::pGetNonlinearSolver(), Predict(), ImplicitStepIntegrator::pXCurr, pXPrev, ImplicitStepIntegrator::pXPrimeCurr, pXPrimePrev, MyVectorHandler::Resize(), SpMapMatrixHandler::Resize(), Restmp, SetCoef(), and NonlinearSolver::Solve().

 {
         if (bAdvanceCalledFirstTime) {
                 integer n = pDM->iGetNumDofs();
                 Restmp.Resize(n);
                 EqIsAlgebraic.resize(n);
                 EqIsDifferential.resize(n);
                 DataManager::DofIterator_const CurrDof = pDofs->begin();
                 for (int iCntm1 = 0; iCntm1 < n;
                         iCntm1++, ++CurrDof)
                 {
                         ASSERT(CurrDof != pDofs->end());
                         EqIsAlgebraic[iCntm1] = (
                                 CurrDof->EqOrder == DofOrder::ALGEBRAIC);
                         EqIsDifferential[iCntm1] = (!EqIsAlgebraic[iCntm1]);
                 }
                 CurrDof = pDofs->begin();
                 Jacxi_xp.Resize(n, n);
                 Jacxi_x.Resize(n, n);
                 Jac_xp.Resize(n, n);
                 Jac_x.Resize(n, n);
                 bAdvanceCalledFirstTime = false;
         }
         pXCurr  = pX;
         pXPrev  = qX[0];
 
         pXPrimeCurr  = pXPrime;
         pXPrimePrev  = qXPrime[0];
         
         SetCoef(TStep, dAlph, StType);  
 
         Predict();
         pDM->LinkToSolution(*pXCurr, *pXPrimeCurr);
         
         Err = 0.;        
         pS->pGetNonlinearSolver()->Solve(this, pS, MaxIters, dTol, 
                         EffIter, Err, dSolTol, SolErr);
         
         return Err;
 };

Here is the call graph for this function:

void ThirdOrderIntegrator::Jacobian ( MatrixHandler * pJac ) const

virtual

Implements NonlinearProblem.

Definition at line 258 of file thirdorderstepsol.cc.

References ASSERT, DataManager::AssJac(), DataManager::AssRes(), MyVectorHandler::Attach(), DEBUGCOUTFNAME, DataManager::dGetTime(), dT, EqIsDifferential, SpMapMatrixHandler::FakeThirdOrderMulAndSumWithShift(), DataManager::iGetNumDofs(), Jac_x, Jac_xp, Jacxi_x, Jacxi_xp, jx11, jx12, jx21, jx22, DataManager::LinkToSolution(), SpMapMatrixHandler::MulAndSumWithShift(), VectorHandler::pdGetVec(), StepIntegrator::pDM, pJac_x, pJac_xp, pJacxi_x, pJacxi_xp, ImplicitStepIntegrator::pXCurr, ImplicitStepIntegrator::pXPrimeCurr, MatrixHandler::Reset(), Restmp, DataManager::SetTime(), theta, and DataManager::Update().

 {
         DEBUGCOUTFNAME("ThirdOrderIntegrator::Jacobian");
         ASSERT(pDM != NULL);
         
         integer iNumDofs = pDM->iGetNumDofs();
 
         MyVectorHandler state, stateder;
         pJac->Reset();
 #if 0
         pJacxi_x->Reset();
         pJacxi_xp->Reset();
         pJac_x->Reset();
         pJac_xp->Reset();
 #endif
         /* theta*dT */
         state.Attach(iNumDofs, pXCurr->pdGetVec()+iNumDofs);
         stateder.Attach(iNumDofs, pXPrimeCurr->pdGetVec()+iNumDofs);
         pDM->SetTime(pDM->dGetTime() + theta*dT);
         pDM->LinkToSolution(state, stateder);
         pDM->Update();
 #warning This trick does not work with Coulomb friction (and other elements too)
         /*
          The reason is that 
                 a) Colomb friction residual (and other elements) can throw
                    to signal the need of a new jacobian
                 b) Have internal states that depend on the order of residual
                    calls
          However, this is needed for all the elements (almost all)
          that compute something during residual and use it 
          later during Jacobain computation
         */
         pDM->AssRes(Restmp, 1.);
         pDM->AssJac(*pJacxi_x, 1.);
         pDM->AssJac(*pJacxi_xp, 0.);
         
         /* dT */
         pDM->SetTime(pDM->dGetTime() - theta*dT);
         pDM->LinkToSolution(*pXCurr, *pXPrimeCurr);
         pDM->Update();
         pDM->AssRes(Restmp, 1.);
         pDM->AssJac(*pJac_x, 1.);
         pDM->AssJac(*pJac_xp, 0.);
         
         
         /* Attenzione: a differenza di quanto riportato a p. 16,
          * "Unconditionally stable multistep integration of ordinary
          * differential and differential-algebraic equations with
          * controlled algorithmic dissipation for multibody dynamic
          * applications"
          * qui il tempo finale e' in cima, il tempo theta in basso
          */ 
          
         /* 2,2 */
         // doublereal J22_x = (1.+3.*rho)/(6.*rho*(1.+rho))*dT;
         Jacxi_x.MulAndSumWithShift(*pJac, jx22, iNumDofs, iNumDofs);
         Jacxi_xp.FakeThirdOrderMulAndSumWithShift(*pJac, EqIsDifferential, 1. - jx22, iNumDofs, iNumDofs);
         
         /* 2,1 */
         // doublereal J21_x = -1./(6.*rho*(1.+rho)*(1.+rho))*dT;
         Jacxi_x.MulAndSumWithShift(*pJac, jx21, iNumDofs, 0);
         Jacxi_xp.FakeThirdOrderMulAndSumWithShift(*pJac, EqIsDifferential, -jx21, iNumDofs, 0);
         
         /* 1,2 */
         // doublereal J12_x = (1.+rho)*(1.+rho)/(6.*rho)*dT;
         Jac_x.MulAndSumWithShift(*pJac, jx12, 0, iNumDofs);
         Jac_xp.FakeThirdOrderMulAndSumWithShift(*pJac, EqIsDifferential, -jx12, 0, iNumDofs);
         
         /* 1,1 */
         // doublereal J11_x = (2.*rho-1.)/(6.*rho)*dT;
         Jac_x.MulAndSumWithShift(*pJac, jx11, 0, 0);
         Jac_xp.FakeThirdOrderMulAndSumWithShift(*pJac, EqIsDifferential, 1. -jx11, 0, 0);
         
         return;
 };

Here is the call graph for this function:

void ThirdOrderIntegrator::Predict ( void )

protectedvirtual

Definition at line 198 of file thirdorderstepsol.cc.

References DataManager::AfterPredict(), ASSERT, DEBUGCOUTFNAME, DataManager::LinkToSolution(), StepIntegrator::pDM, StepIntegrator::pDofs, PredictDof_for_AfterPredict(), ImplicitStepIntegrator::pXCurr, ImplicitStepIntegrator::pXPrimeCurr, RealPredictDof(), and StepIntegrator::UpdateLoop().

Referenced by Advance().

 {
         DEBUGCOUTFNAME("ThirdOrderIntegrator::Predict");
         ASSERT(pDM != NULL);
 
 #ifdef USE_SCHUR
         SchurDataManager* pSDM;
         if ((pSDM = dynamic_cast<SchurDataManager*> (pDM)) != 0) {
                 silent_cerr("Warning: ThirdOrderIntegrator currently is "
                         << "untested with the parallel solver" << std::endl);
         }
 #endif // USE_SCHUR
 
         DataManager::DofIterator_const CurrDof = pDofs->begin();
         
         /* 
          * Linear combination of previous step state and derivative 
          * etc....
         */
         /*
          * Predict for AfterPredict
         */
         UpdateLoop(this, &ThirdOrderIntegrator::PredictDof_for_AfterPredict);
         pDM->LinkToSolution(*pXCurr, *pXPrimeCurr);
         pDM->AfterPredict();
 
         /*
          * Vero Predict
          */
         UpdateLoop(this, &ThirdOrderIntegrator::RealPredictDof);
         return;
 };

Here is the call graph for this function:

void ThirdOrderIntegrator::PredictDof_for_AfterPredict	(	const int	DCount,
		const DofOrder::Order	Order,
		const VectorHandler *const	pSol = `0`
	)		const

protected

Definition at line 124 of file thirdorderstepsol.cc.

References DofOrder::ALGEBRAIC, DofOrder::DIFFERENTIAL, dT, MBDYN_EXCEPT_ARGS, VectorHandler::PutCoef(), ImplicitStepIntegrator::pXCurr, pXPrev, ImplicitStepIntegrator::pXPrimeCurr, and pXPrimePrev.

Referenced by Predict().

                                                        {
         if (Order == DofOrder::DIFFERENTIAL) {
                 doublereal dXnm1 = pXPrev->operator()(DCount);
                 doublereal dXPnm1 = 
                         pXPrimePrev->operator()(DCount);
                         
                 doublereal dXn = dXnm1+dXPnm1*dT;
                 doublereal dXPn = dXPnm1;
         
                 pXPrimeCurr->PutCoef(DCount, dXPn);
                 pXCurr->PutCoef(DCount, dXn);
                 
         } else if (Order == DofOrder::ALGEBRAIC) {
                 doublereal dXnm1 = pXPrev->operator()(DCount);
                 //doublereal dXInm1 = pXPrimePrev->operator()(DCount);
                 doublereal dXn = dXnm1;
                 doublereal dXIn = dXnm1*dT;
                 
                 pXCurr->PutCoef(DCount, dXn);
                 pXPrimeCurr->PutCoef(DCount, dXIn);
 
         } else {
                 silent_cerr("ThirdOrderIntegrator::PredictDof_for_AfterPredict:"
                         << "unknown order for dof " 
                         << DCount<< std::endl);
                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
         }
 }

Here is the call graph for this function:

void ThirdOrderIntegrator::RealPredictDof	(	const int	DCount,
		const DofOrder::Order	Order,
		const VectorHandler *const	pSol = `0`
	)		const

protected

Definition at line 154 of file thirdorderstepsol.cc.

References DofOrder::ALGEBRAIC, DofOrder::DIFFERENTIAL, dT, DataManager::iGetNumDofs(), VectorHandler::IncCoef(), MBDYN_EXCEPT_ARGS, StepIntegrator::pDM, VectorHandler::PutCoef(), ImplicitStepIntegrator::pXCurr, pXPrev, ImplicitStepIntegrator::pXPrimeCurr, pXPrimePrev, and theta.

Referenced by Predict().

                                                        {
         integer iNumDofs = pDM->iGetNumDofs();
         //simple copy of predicted state
         pXPrimeCurr->PutCoef(DCount+iNumDofs,
                 pXPrimeCurr->operator()(DCount));
         pXCurr->PutCoef(DCount+iNumDofs, pXCurr->operator()(DCount));
         if (Order == DofOrder::DIFFERENTIAL) {
                 //doublereal dXPnm1 = 
                 //      pXPrimePrev->operator()(DCount);
                 
                 /* tempo theta*/
 //              doublereal dXn = dXPnm1*(theta-1.)*dT;
 //              pXCurr->IncCoef(DCount+iNumDofs, dXn);
 //              pXCurr->PutCoef(DCount+iNumDofs,
 //                      m0*pXCurr->operator()(DCount)+m1*pXPrev->operator()(DCount)
 //                      +dT*(n0*pXPrimeCurr->operator()(DCount)+
 //                              n1*pXPrimePrev->operator()(DCount)));
                 pXCurr->IncCoef(DCount+iNumDofs,
                         pXPrimePrev->operator()(DCount)*theta*dT);
         } else if (Order == DofOrder::ALGEBRAIC) {
                 //doublereal dXnm1 = pXPrev->operator()(DCount);
                 
                 /* tempo theta*/
 //              doublereal dXIn = dXnm1*(theta-1.)*dT;
 //              pXPrimeCurr->IncCoef(DCount+iNumDofs, dXIn);
 //              pXPrimeCurr->PutCoef(DCount+iNumDofs,
 //                      m0*pXPrimeCurr->operator()(DCount)+m1*pXPrimePrev->operator()(DCount)
 //                      +dT*(n0*pXCurr->operator()(DCount)+
 //                              n1*pXPrev->operator()(DCount)));
 //              pXCurr->PutCoef(DCount+iNumDofs,
 //                      pXPrev->operator()(DCount));
                 pXPrimeCurr->IncCoef(DCount+iNumDofs,
                         pXPrev->operator()(DCount)*theta*dT);
         } else {
                 silent_cerr("ThirdOrderIntegrator::RealPredictDof: "
                         << "unknown order for dof " 
                         << DCount<< std::endl);
                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
         }
 }

Here is the call graph for this function:

void ThirdOrderIntegrator::Residual ( VectorHandler * pRes ) const

virtual

Implements NonlinearProblem.

Definition at line 231 of file thirdorderstepsol.cc.

References ASSERT, DataManager::AssRes(), MyVectorHandler::Attach(), DEBUGCOUTFNAME, DataManager::dGetTime(), dT, DataManager::iGetNumDofs(), DataManager::LinkToSolution(), VectorHandler::pdGetVec(), StepIntegrator::pDM, ImplicitStepIntegrator::pXCurr, ImplicitStepIntegrator::pXPrimeCurr, DataManager::SetTime(), theta, and DataManager::Update().

 {
         DEBUGCOUTFNAME("ThirdOrderIntegrator::Residual");
         ASSERT(pDM != NULL);
         
         integer iNumDofs = pDM->iGetNumDofs();
 
         MyVectorHandler state, stateder, res;
         
         /* theta*dT */
         state.Attach(iNumDofs, pXCurr->pdGetVec()+iNumDofs);
         stateder.Attach(iNumDofs, pXPrimeCurr->pdGetVec()+iNumDofs);
         res.Attach(iNumDofs, pRes->pdGetVec()+iNumDofs);
         pDM->SetTime(pDM->dGetTime() + theta*dT);
         pDM->LinkToSolution(state, stateder);
         pDM->Update();
         pDM->AssRes(res, 1.);
         
         /* dT */
         pDM->SetTime(pDM->dGetTime() - theta*dT);
         pDM->LinkToSolution(*pXCurr, *pXPrimeCurr);
         pDM->Update();
         pDM->AssRes(*pRes, 1.);
 
         return;
 };

Here is the call graph for this function:

virtual void ThirdOrderIntegrator::SetCoef	(	doublereal	dT,
		doublereal	dAlpha,
		enum StepChange	NewStep
	)

protectedpure virtual

Implemented in AdHocThirdOrderIntegrator, and TunableThirdOrderIntegrator.

Referenced by Advance().

void ThirdOrderIntegrator::Update ( const VectorHandler * pSol ) const

virtual

Implements NonlinearProblem.

Definition at line 364 of file thirdorderstepsol.cc.

References ASSERT, DEBUGCOUTFNAME, StepIntegrator::pDM, DataManager::Update(), UpdateDof(), and StepIntegrator::UpdateLoop().

 {
         DEBUGCOUTFNAME("ThirdOrderIntegrator::Predict");
         ASSERT(pDM != NULL);
 
 #ifdef USE_SCHUR
         SchurDataManager* pSDM;
         if ((pSDM = dynamic_cast<SchurDataManager*> (pDM)) != 0) {
                 silent_cerr("Warning: ThirdOrderIntegrator is untested "
                         << "with the parallel solver" << std::endl);
         }
 #endif // USE_SCHUR
 
         UpdateLoop(this,&ThirdOrderIntegrator::UpdateDof,pSol); 
         pDM->Update();
         return;
 };

Here is the call graph for this function:

void ThirdOrderIntegrator::UpdateDof	(	const int	DCount,
		const DofOrder::Order	Order,
		const VectorHandler *const	pSol = `0`
	)		const

protected

Definition at line 334 of file thirdorderstepsol.cc.

References DofOrder::ALGEBRAIC, DofOrder::DIFFERENTIAL, dT, DataManager::iGetNumDofs(), VectorHandler::IncCoef(), m0, MBDYN_EXCEPT_ARGS, n0, StepIntegrator::pDM, ImplicitStepIntegrator::pXCurr, ImplicitStepIntegrator::pXPrimeCurr, w0, and w1.

Referenced by Update().

                                                {
         integer iNumDofs = pDM->iGetNumDofs();
         doublereal dxp = pSol->operator()(DCount);
         doublereal dxp_xi = pSol->operator()(DCount+iNumDofs);
         if (Order == DofOrder::DIFFERENTIAL) {
                 
                 pXPrimeCurr->IncCoef(DCount, dxp);
                 pXPrimeCurr->IncCoef(DCount+iNumDofs, dxp_xi);
                 
                 pXCurr->IncCoef(DCount, dT*(w1*dxp_xi+w0*dxp));
                 pXCurr->IncCoef(DCount+iNumDofs, 
                         dT*(m0*w1*dxp_xi+(m0*w0+n0)*dxp));
         
         } else if (Order == DofOrder::ALGEBRAIC) {
                 pXCurr->IncCoef(DCount, dxp);
                 pXCurr->IncCoef(DCount+iNumDofs, dxp_xi);
                 
                 pXPrimeCurr->IncCoef(DCount, dT*(w1*dxp_xi+w0*dxp));
                 pXPrimeCurr->IncCoef(DCount+iNumDofs, 
                         dT*(m0*w1*dxp_xi+(m0*w0+n0)*dxp));
         } else {
                 silent_cerr("unknown order for dof " 
                         << DCount<< std::endl);
                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
         }
 };