incrementallinearstatic.C

Go to the documentation of this file.

/*
 *
 *                 #####    #####   ######  ######  ###   ###
 *               ##   ##  ##   ##  ##      ##      ## ### ##
 *              ##   ##  ##   ##  ####    ####    ##  #  ##
 *             ##   ##  ##   ##  ##      ##      ##     ##
 *            ##   ##  ##   ##  ##      ##      ##     ##
 *            #####    #####   ##      ######  ##     ##
 *
 *
 *             OOFEM : Object Oriented Finite Element Code
 *
 *               Copyright (C) 1993 - 2013   Borek Patzak
 *
 *
 *
 *       Czech Technical University, Faculty of Civil Engineering,
 *   Department of Structural Mechanics, 166 29 Prague, Czech Republic
 *
 *  This library is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU Lesser General Public
 *  License as published by the Free Software Foundation; either
 *  version 2.1 of the License, or (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public
 *  License along with this library; if not, write to the Free Software
 *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 */

#include "../sm/EngineeringModels/incrementallinearstatic.h"
#include "timestep.h"
#include "dof.h"
#include "domain.h"
#include "sparsemtrx.h"
#include "dictionary.h"
#include "verbose.h"
#include "classfactory.h"
#include "datastream.h"
#include "contextioerr.h"
#include "dofmanager.h"
#include "activebc.h"
#include "unknownnumberingscheme.h"

/*
#include "set.h"
#include "element.h"
#include "node.h"
*/

#include "boundarycondition.h"

#include <vector>
#include <set>

namespace oofem {
REGISTER_EngngModel(IncrementalLinearStatic);

IncrementalLinearStatic :: IncrementalLinearStatic(int i, EngngModel *_master) : StructuralEngngModel(i, _master),
    loadVector(), internalLoadVector(), incrementOfDisplacementVector(), discreteTimes()
{
    ndomains = 1;
    endOfTimeOfInterest = 0;
}


IncrementalLinearStatic :: ~IncrementalLinearStatic()
{}


NumericalMethod *IncrementalLinearStatic :: giveNumericalMethod(MetaStep *mStep)

{
    if ( !nMethod ) {
        nMethod.reset( classFactory.createSparseLinSolver(solverType, this->giveDomain(1), this) );
        if ( !nMethod ) {
            OOFEM_ERROR("linear solver creation failed for lstype %d", solverType);
        }
    }

    return nMethod.get();
}


IRResultType IncrementalLinearStatic :: initializeFrom(InputRecord *ir)
{
    IRResultType result;

    IR_GIVE_OPTIONAL_FIELD(ir, discreteTimes, _IFT_IncrementalLinearStatic_prescribedtimes);
    if ( discreteTimes.giveSize() > 0 ) {
        numberOfSteps = discreteTimes.giveSize();
        endOfTimeOfInterest = discreteTimes.at( discreteTimes.giveSize() );
        fixedSteps = false;
    } else {
        deltaT = 1.0;
        IR_GIVE_OPTIONAL_FIELD(ir, deltaT, _IFT_IncrementalLinearStatic_deltat);
        IR_GIVE_FIELD(ir, numberOfSteps, _IFT_EngngModel_nsteps);
        endOfTimeOfInterest = deltaT * numberOfSteps;
        fixedSteps = true;
    }
    IR_GIVE_OPTIONAL_FIELD(ir, endOfTimeOfInterest, _IFT_IncrementalLinearStatic_endoftimeofinterest);

    int val = 0;
    IR_GIVE_OPTIONAL_FIELD(ir, val, _IFT_EngngModel_lstype);
    solverType = ( LinSystSolverType ) val;

    val = 0;
    IR_GIVE_OPTIONAL_FIELD(ir, val, _IFT_EngngModel_smtype);
    sparseMtrxType = ( SparseMtrxType ) val;


    suppressOutput = ir->hasField(_IFT_EngngModel_suppressOutput);

    if(suppressOutput) {
      printf("Suppressing output.\n");
    }
    else {

    if ( ( outputStream = fopen(this->dataOutputFileName.c_str(), "w") ) == NULL ) {
      OOFEM_ERROR("Can't open output file %s", this->dataOutputFileName.c_str());
    }

    fprintf(outputStream, "%s", PRG_HEADER);
    fprintf(outputStream, "\nStarting analysis on: %s\n", ctime(& this->startTime) );
    fprintf(outputStream, "%s\n", simulationDescription.c_str());
  }

    //StructuralEngngModel::initializeFrom (ir);
    return IRRT_OK;
}


double IncrementalLinearStatic :: giveDiscreteTime(int iStep)
{
    if ( this->fixedSteps ) {
        return this->deltaT * iStep;
    } else {
        if ( ( iStep > 0 ) && ( iStep <= discreteTimes.giveSize() ) ) {
            return ( discreteTimes.at(iStep) );
        }
    }

    OOFEM_ERROR("invalid iStep");
    return 0.0;
}


TimeStep *IncrementalLinearStatic :: giveNextStep()
{
    if ( !currentStep ) {
        currentStep.reset( new TimeStep(giveNumberOfTimeStepWhenIcApply(), this, 0, 0.,  this->giveDiscreteTime(1), 0) );
    }

    previousStep = std :: move(currentStep);
    double dt = this->giveDiscreteTime(previousStep->giveNumber()+1) - previousStep->giveTargetTime();
    currentStep.reset( new TimeStep(*previousStep, dt) );
    return currentStep.get();
}

void IncrementalLinearStatic :: solveYourself()
{
    this->giveDiscreteTime(1); // ensures numberOfSteps defined
    StructuralEngngModel :: solveYourself();
}

void IncrementalLinearStatic :: solveYourselfAt(TimeStep *tStep)
{
    Domain *d = this->giveDomain(1);
    // Creates system of governing eq's and solves them at given time step


    // >>> beginning PH    
    // The following piece of code updates assignment of boundary conditions to dofs
    // (this allows to have multiple boundary conditions assigned to one dof
    // which can be arbitrarily turned on and off in time)
    // Almost the entire section has been copied from domain.C
    std :: vector< std :: map< int, int > > dof_bc( d->giveNumberOfDofManagers() );

    for ( int i = 1; i <= d->giveNumberOfBoundaryConditions(); ++i ) {
      GeneralBoundaryCondition *gbc = d->giveBc(i);
      
      if ( gbc->isImposed(tStep) ){

  if ( gbc->giveSetNumber() > 0 ) { 
    // Loop over nodes in set and store the bc number in each dof.
    Set *set = d->giveSet( gbc->giveSetNumber() );
    ActiveBoundaryCondition *active_bc = dynamic_cast< ActiveBoundaryCondition * >(gbc);
    BoundaryCondition *bc = dynamic_cast< BoundaryCondition * >(gbc);
    if ( bc || ( active_bc && active_bc->requiresActiveDofs() ) ) {
      const IntArray &appliedDofs = gbc->giveDofIDs();
      const IntArray &nodes = set->giveNodeList();
      for ( int inode = 1; inode <= nodes.giveSize(); ++inode ) {
        for ( int idof = 1; idof <= appliedDofs.giveSize(); ++idof ) {
    
    if  ( dof_bc [ nodes.at(inode) - 1 ].find( appliedDofs.at(idof) ) == dof_bc [ nodes.at(inode) - 1 ].end() ) {
      // is empty
      dof_bc [ nodes.at(inode) - 1 ] [ appliedDofs.at(idof) ] = i;

      DofManager * dofman = d->giveDofManager( nodes.at(inode) );
      Dof * dof = dofman->giveDofWithID( appliedDofs.at(idof) );

      dof->setBcId(i);

    } else {
      // another bc has been already prescribed at this time step to this dof
      OOFEM_WARNING("More than one boundary condition assigned at time %f to node %d dof %d. Considering boundary condition %d", tStep->giveTargetTime(),  nodes.at(inode), appliedDofs.at(idof), dof_bc [ nodes.at(inode) - 1 ] [appliedDofs.at(idof)] );
      
      
    }
        }
      }
    }
  }
      }
    }
    
    // to get proper number of equations
    this->forceEquationNumbering();
    // <<< end PH



    // Initiates the total displacement to zero.
    if ( tStep->isTheFirstStep() ) {
        for ( auto &dofman : d->giveDofManagers() ) {
            for ( Dof *dof: *dofman ) {
                dof->updateUnknownsDictionary(tStep->givePreviousStep(), VM_Total, 0.);
                dof->updateUnknownsDictionary(tStep, VM_Total, 0.);
            }
        }

        for ( auto &bc : d->giveBcs() ) {
            ActiveBoundaryCondition *abc;

            if ( ( abc = dynamic_cast< ActiveBoundaryCondition * >(bc.get()) ) ) {
                int ndman = abc->giveNumberOfInternalDofManagers();
                for ( int i = 1; i <= ndman; i++ ) {
                    DofManager *dofman = abc->giveInternalDofManager(i);
                    for ( Dof *dof: *dofman ) {
                        dof->updateUnknownsDictionary(tStep->givePreviousStep(), VM_Total, 0.);
                        dof->updateUnknownsDictionary(tStep, VM_Total, 0.);
                    }
                }
            }
        }
    }

    // Apply dirichlet b.c's on total values
    for ( auto &dofman : d->giveDofManagers() ) {
        for ( Dof *dof: *dofman ) {
            double tot = dof->giveUnknown( VM_Total, tStep->givePreviousStep() );
            if ( dof->hasBc(tStep) ) {
                tot += dof->giveBcValue(VM_Incremental, tStep);
            }

            dof->updateUnknownsDictionary(tStep, VM_Total, tot);
        }
    }

    int neq = this->giveNumberOfDomainEquations( 1, EModelDefaultEquationNumbering() );

#ifdef VERBOSE
    OOFEM_LOG_RELEVANT("Solving [step number %8d, time %15e, equations %d]\n", tStep->giveNumber(), tStep->giveTargetTime(), neq);
#endif

    if ( neq == 0 ) { // Allows for fully prescribed/empty problems.
        return;
    }

    incrementOfDisplacementVector.resize(neq);
    incrementOfDisplacementVector.zero();

#ifdef VERBOSE
    OOFEM_LOG_INFO("Assembling load\n");
#endif
    // Assembling the element part of load vector
    internalLoadVector.resize(neq);
    internalLoadVector.zero();
    this->assembleVector( internalLoadVector, tStep, InternalForceAssembler(),
                         VM_Total, EModelDefaultEquationNumbering(), this->giveDomain(1) );

    loadVector.resize(neq);
    loadVector.zero();
    this->assembleVector( loadVector, tStep, ExternalForceAssembler(),
                         VM_Total, EModelDefaultEquationNumbering(), this->giveDomain(1) );

    loadVector.subtract(internalLoadVector);
    this->updateSharedDofManagers(loadVector, EModelDefaultEquationNumbering(), ReactionExchangeTag);


#ifdef VERBOSE
    OOFEM_LOG_INFO("Assembling stiffness matrix\n");
#endif
    stiffnessMatrix.reset( classFactory.createSparseMtrx(sparseMtrxType) );
    if ( !stiffnessMatrix ) {
        OOFEM_ERROR("sparse matrix creation failed");
    }

    stiffnessMatrix->buildInternalStructure( this, 1, EModelDefaultEquationNumbering() );
    stiffnessMatrix->zero();
    this->assemble( *stiffnessMatrix, tStep, TangentAssembler(TangentStiffness),
                   EModelDefaultEquationNumbering(), this->giveDomain(1) );

#ifdef VERBOSE
    OOFEM_LOG_INFO("Solving ...\n");
#endif
    this->giveNumericalMethod( this->giveCurrentMetaStep() );
    NM_Status s = nMethod->solve(*stiffnessMatrix, loadVector, incrementOfDisplacementVector);
    if ( !( s & NM_Success ) ) {
        OOFEM_ERROR("No success in solving system.");
    }
}


double IncrementalLinearStatic :: giveUnknownComponent(ValueModeType mode, TimeStep *tStep, Domain *d, Dof *dof)
{
    if ( this->requiresUnknownsDictionaryUpdate() ) {
        int hash = this->giveUnknownDictHashIndx(mode, tStep);
        if ( dof->giveUnknowns()->includes(hash) ) {
            return dof->giveUnknowns()->at(hash);
        } else {
            OOFEM_ERROR("Dof unknowns dictionary does not contain unknown of value mode (%s)", __ValueModeTypeToString(mode) );
        }
    } else {
        OOFEM_ERROR("Only the mode requiresUnknownsDictionaryUpdate() is supported");
    }

    return 0.;
}


int IncrementalLinearStatic :: giveUnknownDictHashIndx(ValueModeType mode, TimeStep *tStep)
{
    return tStep->giveNumber() % 2;
}


void IncrementalLinearStatic :: updateDofUnknownsDictionary(DofManager *inode, TimeStep *tStep)
{
    // update DOF unknowns dictionary, where
    // unknowns are hold instead of keeping them in global unknowns
    // vectors in engng instances
    // this is necessary, because during solution equation numbers for
    // particular DOFs may changed, and it is necessary to keep them
    // in DOF level.

    double val;
    for ( Dof *dof: *inode ) {
        // skip slave DOFs (only master (primary) DOFs have to be updated).
        if ( !dof->isPrimaryDof() ) {
            continue;
        }
        val = dof->giveUnknown(VM_Total, tStep);
        if ( !dof->hasBc(tStep) ) {
            val += this->incrementOfDisplacementVector.at( dof->__giveEquationNumber() );
        }

        dof->updateUnknownsDictionary(tStep->givePreviousStep(), VM_Total, val);
        dof->updateUnknownsDictionary(tStep, VM_Total, val);
    }
}


contextIOResultType IncrementalLinearStatic :: saveContext(DataStream &stream, ContextMode mode)
{
    contextIOResultType iores;

    if ( ( iores = StructuralEngngModel :: saveContext(stream, mode) ) != CIO_OK ) {
        THROW_CIOERR(iores);
    }

    return CIO_OK;
}


contextIOResultType IncrementalLinearStatic :: restoreContext(DataStream &stream, ContextMode mode)
{
    contextIOResultType iores;

    if ( ( iores = StructuralEngngModel :: restoreContext(stream, mode) ) != CIO_OK ) {
        THROW_CIOERR(iores);
    }

    return CIO_OK;
}
} // end namespace oofem