Next: Material Element Up: Engineering model Previous: Grid Representation - Domain Contents
class LinearStatic : public StructuralEngngModel
{
/*
This class implements LinearStatic Engineering problem.
Multiple loading works only if linear elastic material (such as isoLE) is used.
(Other non-linear materials kepp load history, so such multiple loading
will cause that next step will be assumed as new load increment,
not the total new load). Because they always copute real stresses acording
to reached strain state, they are not able to respond to linear analysis.
DESCRIPTION:
Solution of this problem is series of loading cases, maintained as sequence of
time-steps. This solution is in form of linear equation system Ax=b
TASK:
Creating Numerical method for solving Ax=b
Interfacing Numerical method to Elements
Managing time steps
*/
protected:
SparseMtrx* stiffnessMatrix;
FloatArray loadVector;
FloatArray displacementVector;
LinSystSolverType solverType;
SparseMtrxType sparseMtrxType;
/// Numerical method used to solve the problem
SparseLinearSystemNM *nMethod;
int initFlag;
#ifdef __PETSC_MODULE
Vec _loadVec, _dispVec;
#endif
public:
LinearStatic (int i, EngngModel* _master = NULL) ;
~LinearStatic ()
{delete stiffnessMatrix; if (nMethod) delete nMethod;}
// solving
void solveYourself ();
void solveYourselfAt (TimeStep *);
//int requiresNewLhs () {return 0;}
/**
Updates nodal values
(calls also this->updateDofUnknownsDictionary for updating dofs unknowns dictionaries
if model supports changes of static system). The element internal state update is also forced using
updateInternalState service.
*/
virtual void updateYourself (TimeStep*) ;
double giveUnknownComponent ( EquationID, ValueModeType, TimeStep*, Domain*, Dof*);
contextIOResultType saveContext (FILE* stream, void *obj = NULL) ;
contextIOResultType restoreContext (FILE* stream, void *obj = NULL);
void updateDomainLinks();
TimeStep* giveNextStep ();
NumericalMethod* giveNumericalMethod (TimeStep*);
//void printReactionForces (TimeStep *);
void terminate (TimeStep*);
IRResultType initializeFrom (InputRecord* ir);
// consistency check
virtual int checkConsistency (); // returns nonzero if o.k.
/** DOF printing routine. Called by DofManagers to print Dof specific part.
Dof class provides component printing routines, but emodel is responsible
for what will be printed at DOF level.
@param stream output stream
@param iDof dof to be processed
@param atTime solution step
*/
virtual void printDofOutputAt (FILE* stream, Dof* iDof, TimeStep* atTime);
// identification
const char* giveClassName () const { return "LinearStatic";}
classType giveClassID () const { return LinearStaticClass;}
fMode giveFormulation () { return TL; }
#ifdef __PARALLEL_MODE
/**
Determines the space necessary for send/receive buffer.
It uses related communication map pattern to determine the maximum size needed.
@param commMap communication map used to send/receive messages
@param buff communication buffer
@return upper bound of space needed
*/
int estimateMaxPackSize (IntArray& commMap, CommunicationBuffer& buff, int packUnpackType) ;
#endif
} ;
The listing of the LinearStatic class implementation follows. First, the giveNumericalMethod returns the numerical method, which will be used. The new instance is created if not defined according to nMethod attribute, which is initialized from input (see instanciateYourself). Note, that thanks to the general interface, defined by the NumericalMethod class and the component mapping approach, the further communication with the numerical method uses only the interface declared by the base NumerialMethod class and therefore there is no need to worry about particular numerical method details.
NumericalMethod* LinearStatic :: giveNumericalMethod (TimeStep* atTime)
// only one has reason for LinearStatic
// - SolutionOfLinearEquations
{
if (nMethod) return nMethod ;
SparseLinearSystemNM* nm;
if (solverType == ST_Direct) {
nm = (SparseLinearSystemNM*) new LDLTFactorization (1,this->giveDomain(1),this);
nMethod = nm;
return nm;
} else {
nm = (SparseLinearSystemNM*) new IMLSolver (1,this->giveDomain(1),this);
nMethod = nm;
return nm;
}
}
IRResultType
LinearStatic :: initializeFrom (InputRecord* ir)
{
// Required by IR_GIVE_FIELD macro
const char *__keyword, *__proc = "initializeFrom";
IRResultType result;
StructuralEngngModel::initializeFrom (ir);
int val = 0;
IR_GIVE_OPTIONAL_FIELD (ir, val, "lstype"); // Macro
solverType = (LinSystSolverType) val;
val = 0;
IR_GIVE_OPTIONAL_FIELD (ir, val, "smtype"); // Macro
sparseMtrxType = (SparseMtrxType) val;
return IRRT_OK;
}
double LinearStatic :: giveUnknownComponent (EquationID chc, ValueModeType mode,
TimeStep* tStep, Domain* d, Dof* dof)
// returns unknown quantity like displaacement, velocity of equation eq
// This function translates this request to numerical method language
{
int eq = dof->giveEquationNumber();
if (eq == 0) _error ("giveUnknownComponent: invalid equation number");
if (tStep != this->giveCurrentStep ()) {
_error ("giveUnknownComponent: unknown time step encountered");
return 0.;
}
if (chc != EID_MomentumBalance) {
_error ("giveUnknownComponent: Unknown is of undefined CharType for this problem");
return 0.;
}
switch (mode)
{
case VM_Total:
case VM_Incremental:
if (displacementVector.isNotEmpty()) return displacementVector.at(eq);
else return 0.;
// return nMethod-> giveUnknownComponent (LinearEquationSolution, eq);
default:
_error ("giveUnknownComponent: Unknown is of undefined type for this problem");
}
return 0.;
}
TimeStep* LinearStatic :: giveNextStep ()
{
int istep = this->giveNumberOfFirstStep();
//int mstep = 1;
StateCounterType counter = 1;
delete previousStep;
if (currentStep != NULL) {
istep = currentStep->giveNumber() + 1 ;
counter = currentStep->giveSolutionStateCounter() + 1;
}
previousStep = currentStep;
currentStep = new TimeStep (istep,this, 1, (double) istep, 0., counter);
// time and dt variables are set eq to 0 for staics - has no meaning
return currentStep;
}
void LinearStatic :: solveYourselfAt (TimeStep* tStep) {
//
// creates system of governing eq's and solves them at given time step
//
// first assemble problem at current time step
if (initFlag) {
#ifdef VERBOSE
OOFEM_LOG_INFO("Assembling stiffness matrix\n");
#endif
//
// first step assemble stiffness Matrix
//
/*
IntArray* mht = this -> GiveBanWidthVector ();
stiffnessMatrix = new Skyline ();
stiffnessMatrix -> checkSizeTowardsBanWidth (mht) ;
delete mht;
*/
stiffnessMatrix = ::CreateUsrDefSparseMtrx(sparseMtrxType); // new Skyline ();
if (stiffnessMatrix==NULL) _error ("solveYourselfAt: sparse matrix creation failed");
//stiffnessMatrix = new DynCompCol ();
//stiffnessMatrix = new CompCol ();
stiffnessMatrix->buildInternalStructure (this, 1, EID_MomentumBalance);
this -> assemble (stiffnessMatrix, tStep, EID_MomentumBalance, StiffnessMatrix, this->giveDomain(1));
//
// alocate space for displacementVector
//
//displacementVector = new FloatArray (this->giveNumberOfEquations());
displacementVector.resize (this->giveNumberOfEquations(EID_MomentumBalance));
displacementVector.zero();
#ifdef __PETSC_MODULE
if (solverType == ST_Petsc) {
this->givePetscContext(1, EID_MomentumBalance)->createVecGlobal (&_loadVec);
this->givePetscContext(1, EID_MomentumBalance)->createVecGlobal (&_dispVec);
}
#endif
initFlag = 0;
}
#ifdef VERBOSE
OOFEM_LOG_INFO("Assembling load\n");
#endif
#ifdef __PETSC_MODULE
// direct interface to PETSC
if (solverType == ST_Petsc) {
this->petsc_assembleVectorFromElements(_loadVec, tStep, EID_MomentumBalance,
ElementForceLoadVector, VM_Total,
this->giveDomain(1));
this->petsc_assembleVectorFromElements(_loadVec, tStep, EID_MomentumBalance,
ElementNonForceLoadVector, VM_Total,
this->giveDomain(1));
this->petsc_assembleVectorFromDofManagers(_loadVec, tStep, EID_MomentumBalance,
NodalLoadVector, VM_Total,
this->giveDomain(1)) ;
VecAssemblyBegin(_loadVec);
VecAssemblyEnd (_loadVec);
this->giveNumericalMethod(tStep);
#ifdef VERBOSE
OOFEM_LOG_INFO("Solving ...\n");
#endif
//nMethod -> solveYourselfAt(tStep);
PetscSolver *ps = dynamic_cast<PetscSolver*>(nMethod);
PetscSparseMtrx* psm = dynamic_cast<PetscSparseMtrx*>(stiffnessMatrix);
ps -> petsc_solve (psm, _loadVec, _dispVec);
this->givePetscContext(1, EID_MomentumBalance)->scatterG2N (_dispVec, &displacementVector, INSERT_VALUES);
} else
#endif
{
//
// assembling the element part of load vector
//
//loadVector = new FloatArray (this->giveNumberOfEquations());
loadVector.resize (this->giveNumberOfEquations(EID_MomentumBalance));
loadVector.zero();
this->assembleVectorFromElements(loadVector, tStep, EID_MomentumBalance,
ElementForceLoadVector, VM_Total, this->giveDomain(1));
this->assembleVectorFromElements(loadVector, tStep, EID_MomentumBalance,
ElementNonForceLoadVector, VM_Total, this->giveDomain(1));
//
// assembling the nodal part of load vector
//
this->assembleVectorFromDofManagers(loadVector, tStep, EID_MomentumBalance,
NodalLoadVector, VM_Total, this->giveDomain(1)) ;
//
// set-up numerical model
//
this->giveNumericalMethod(tStep);
/*
nMethod -> setSparseMtrxAsComponent ( LinearEquationLhs , stiffnessMatrix) ;
nMethod -> setFloatArrayAsComponent ( LinearEquationRhs , &loadVector) ;
nMethod -> setFloatArrayAsComponent ( LinearEquationSolution, &displacementVector) ;
*/
//
// call numerical model to solve arised problem
//
#ifdef VERBOSE
OOFEM_LOG_INFO("Solving ...\n");
#endif
//nMethod -> solveYourselfAt(tStep);
nMethod -> solve (stiffnessMatrix, &loadVector, &displacementVector);
}
tStep->incrementStateCounter(); // update solution state counter
//
// update nodes, elements, etc.
this->updateYourself(this->giveCurrentStep());
}
void LinearStatic :: updateYourself (TimeStep* stepN)
{
this->updateInternalState(stepN);
StructuralEngngModel::updateYourself(stepN);
}
contextIOResultType LinearStatic :: saveContext (FILE* stream, void *obj)
//
// saves state variable - displacement vector
//
{
contextIOResultType iores;
int closeFlag = 0;
if (stream==NULL) {
if (!this->giveContextFile(&stream,this->giveCurrentStep()->giveNumber(),
contextMode_write))
THROW_CIOERR(CIO_IOERR); // override
closeFlag = 1;
}
if ((iores = StructuralEngngModel :: saveContext (stream)) != CIO_OK)
THROW_CIOERR(iores);
if ((iores = displacementVector.storeYourself(stream)) != CIO_OK)
THROW_CIOERR(iores);
if (closeFlag) fclose (stream); // ensure consistent records
return CIO_OK;
}
contextIOResultType LinearStatic :: restoreContext (FILE* stream, void *obj)
//
// restore state variable - displacement vector
//
{
contextIOResultType iores;
int closeFlag = 0;
int istep = this->resolveCorrespondingStepNumber (obj);
if (stream == NULL) {
if (!this->giveContextFile(&stream, istep, contextMode_read))
THROW_CIOERR(CIO_IOERR); // override
closeFlag = 1;
}
if ((iores=StructuralEngngModel::restoreContext(stream,obj))!=CIO_OK)
THROW_CIOERR(iores);
if ((iores=displacementVector.restoreYourself(stream))!=CIO_OK)
THROW_CIOERR(iores);
if (closeFlag) fclose (stream); // ensure consistent records
return CIO_OK;
}
void
LinearStatic:: terminate(TimeStep* tStep)
{
StructuralEngngModel :: terminate (tStep);
this->printReactionForces (tStep, 1);
}
int
LinearStatic::checkConsistency ()
{
// check internal consistency
// if success returns nonzero
int i, nelem;
Element* ePtr;
StructuralElement* sePtr;
Domain* domain = this->giveDomain(1);
nelem = domain->giveNumberOfElements();
// check for proper element type
for (i=1; i<= nelem; i++) {
ePtr = domain->giveElement(i);
sePtr = DYNAMIC_CAST (StructuralElement, ePtr);
if (sePtr == NULL) {
printf ("Error: Element %d has no StructuralElement base\n",i);
return 0;
}
}
EngngModel :: checkConsistency ();
return 1;
}
void
LinearStatic::updateDomainLinks ()
{
EngngModel::updateDomainLinks();
this->giveNumericalMethod(giveCurrentStep())->setDomain (this->giveDomain(1));
}
Borek Patzak 2018-01-02