linearconstraintbc.C

Go to the documentation of this file.

/*
 *
 *                 #####    #####   ######  ######  ###   ###
 *               ##   ##  ##   ##  ##      ##      ## ### ##
 *              ##   ##  ##   ##  ####    ####    ##  #  ##
 *             ##   ##  ##   ##  ##      ##      ##     ##
 *            ##   ##  ##   ##  ##      ##      ##     ##
 *            #####    #####   ##      ######  ##     ##
 *
 *
 *             OOFEM : Object Oriented Finite Element Code
 *
 *               Copyright (C) 1993 - 2025   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 "linearconstraintbc.h"
#include "classfactory.h"
#include "masterdof.h"
#include "floatmatrix.h"
#include "sparsemtrx.h"
#include "unknownnumberingscheme.h"
#include "function.h"
#include "timestep.h"
#include "datastream.h"
#include "contextioerr.h"
#include "node.h"
#include "domain.h"
 
#ifdef _OPENMP
#include <omp.h>
#endif
 
namespace oofem {
REGISTER_BoundaryCondition(LinearConstraintBC);
 
LinearConstraintBC :: LinearConstraintBC(int n, Domain *d) : ActiveBoundaryCondition(n, d),
    md( new Node(0, domain) )
{
    // this is internal lagrange multiplier used to enforce the receiver constrain
    // this allocates a new equation related to this constraint
    this->md->appendDof( new MasterDof( this->md.get(), ( DofIDItem ) ( d->giveNextFreeDofID() ) ) );
    this->lhsType.clear();
    this->rhsType.clear();
}
 
 
LinearConstraintBC :: ~LinearConstraintBC()
{
}
 
 
void LinearConstraintBC :: initializeFrom(InputRecord &ir)
{
    ActiveBoundaryCondition :: initializeFrom(ir);
    rhsTf = 0;
 
    IR_GIVE_FIELD(ir, weights, _IFT_LinearConstraintBC_weights);
    IR_GIVE_FIELD(ir, rhs, _IFT_LinearConstraintBC_rhs);
    IR_GIVE_FIELD(ir, dofmans, _IFT_LinearConstraintBC_dofmans);
    IR_GIVE_FIELD(ir, dofs, _IFT_LinearConstraintBC_dofs);
    if ( weights.giveSize() != dofmans.giveSize() ) {
        throw ValueInputException(ir, _IFT_LinearConstraintBC_weights, "Size mismatch weights and dofmans");
    }
    IR_GIVE_OPTIONAL_FIELD(ir, weightsTf, _IFT_LinearConstraintBC_weightsfuncs);
    IR_GIVE_OPTIONAL_FIELD(ir, rhsTf, _IFT_LinearConstraintBC_rhsfuncs);
 
    IR_GIVE_FIELD(ir, lhsType, _IFT_LinearConstraintBC_lhstype);
    IR_GIVE_FIELD(ir, rhsType, _IFT_LinearConstraintBC_rhstype);
}
 
 
void LinearConstraintBC :: giveLocArray(const UnknownNumberingScheme &r_s,  IntArray &locr, int &lambda_eq)
{
    int size = this->weights.giveSize();
 
    locr.resize(size);
    // assemble location array
    for ( int _i = 1; _i <= size; _i++ ) {
        Dof *idof = this->domain->giveDofManager( this->dofmans.at(_i) )->giveDofWithID( this->dofs.at(_i) );
        locr.at(_i) = r_s.giveDofEquationNumber(idof);
    }
 
    lambda_eq = r_s.giveDofEquationNumber( *md->begin() );
}
 
 
void LinearConstraintBC :: assemble(SparseMtrx &answer, TimeStep *tStep,
                                    CharType type, const UnknownNumberingScheme &r_s,
                                    const UnknownNumberingScheme &c_s, double scale,
                                    void*lock)
{
    int size = this->weights.giveSize();
    IntArray lambdaeq(1);
    FloatMatrix contrib(size, 1), contribt;
    IntArray locr(size);
 
    if ( !this->lhsType.contains( ( int ) type ) ) {
        return;
    }
    this->giveLocArray( r_s, locr, lambdaeq.at(1) );
 
    if ( this->isImposed(tStep) ) {
        for ( int _i = 1; _i <= size; _i++ ) { // loop over dofs
            double factor = 1.;
            if ( weightsTf.giveSize() ) {
                factor = domain->giveFunction( weightsTf.at(_i) )->evaluateAtTime( tStep->giveIntrinsicTime() );
            }
            contrib.at(_i, 1) = this->weights.at(_i) * factor;
        }
 
        contribt.times(scale);
        contribt.beTranspositionOf(contrib);
#ifdef _OPENMP
        if (lock) omp_set_lock(static_cast<omp_lock_t*>(lock));
#endif
        answer.assemble(lambdaeq, locr, contribt);
        answer.assemble(locr, lambdaeq, contrib);
#ifdef _OPENMP
        if (lock) omp_unset_lock(static_cast<omp_lock_t*>(lock));
#endif
    } else {
        // the bc is not imposed at specific time step, however in order to make the equation system regular
        // we initialize the allocated equation to the following form 1*labmda = 0, forcing lagrange multiplier
        // of inactive condition to be zero.
        FloatMatrix help(1, 1);
        help.at(1, 1) = 1.0;
#ifdef _OPENMP
        if (lock) omp_set_lock(static_cast<omp_lock_t*>(lock));
#endif
        answer.assemble(lambdaeq, lambdaeq, help);
#ifdef _OPENMP
        if (lock) omp_unset_lock(static_cast<omp_lock_t*>(lock));
#endif
    }
}
 
void LinearConstraintBC :: assembleVector(FloatArray &answer, TimeStep *tStep,
                                          CharType type, ValueModeType mode,
                                          const UnknownNumberingScheme &s, 
                                          FloatArray *eNorms,
                                          void* lock)
{
    IntArray loc, lambdaeq(1);
    FloatArray vec(1);
    double factor = 1.;
 
    if ( !this->rhsType.contains( ( int ) type ) ) {
        return;
    }
 
    if ( type == InternalForcesVector ) {
        // compute true residual
        int size = this->weights.giveSize();
        Dof *mdof = *md->begin();
        Dof *idof;
 
        // assemble location array
        for ( int _i = 1; _i <= size; _i++ ) {
            factor = 1.;
            if ( weightsTf.giveSize() ) {
                factor = domain->giveFunction( weightsTf.at(_i) )->evaluateAtTime( tStep->giveIntrinsicTime() );
            }
            idof = this->domain->giveDofManager( this->dofmans.at(_i) )->giveDofWithID( this->dofs.at(_i) );
            if ( s.giveDofEquationNumber(idof) ) {
#ifdef _OPENMP
                if (lock) omp_set_lock(static_cast<omp_lock_t*>(lock));
#endif
                answer.at( s.giveDofEquationNumber(idof) ) += mdof->giveUnknown(mode, tStep) * this->weights.at(_i) * factor;
#ifdef _OPENMP
                if (lock) omp_unset_lock(static_cast<omp_lock_t*>(lock));
#endif
            }
            if ( s.giveDofEquationNumber( mdof ) ) {
#ifdef _OPENMP
                if (lock) omp_set_lock(static_cast<omp_lock_t*>(lock));
#endif
                answer.at( s.giveDofEquationNumber( mdof ) ) += idof->giveUnknown(mode, tStep) * this->weights.at(_i) * factor;
#ifdef _OPENMP
                if (lock) omp_unset_lock(static_cast<omp_lock_t*>(lock));
#endif
            }
        }
    } else if ( type == ExternalForcesVector ) {
        // use rhs value
 
        if ( rhsTf ) {
            factor = domain->giveFunction(rhsTf)->evaluateAtTime( tStep->giveIntrinsicTime() );
        }
        this->giveLocArray( s, loc, lambdaeq.at(1) );
        vec.at(1) = rhs * factor;
#ifdef _OPENMP
        if (lock) omp_set_lock(static_cast<omp_lock_t*>(lock));
#endif
        answer.assemble(vec, lambdaeq);
#ifdef _OPENMP
        if (lock) omp_unset_lock(static_cast<omp_lock_t*>(lock));
#endif
    }
}
 
void LinearConstraintBC :: giveLocationArrays(std :: vector< IntArray > &rows, std :: vector< IntArray > &cols, CharType type, const UnknownNumberingScheme &r_s, const UnknownNumberingScheme &c_s)
{
    rows.resize(3);
    cols.resize(3);
 
    IntArray loc, lambdaeq(1);
    this->giveLocArray( r_s, loc, lambdaeq.at(1) );
    // column block
    rows [ 0 ] = loc;
    cols [ 0 ] = lambdaeq;
    // row block
    cols [ 1 ] = loc;
    rows [ 1 ] = lambdaeq;
    // diagonal enry (some sparse mtrx implementation requaire this)
    rows [ 2 ] = lambdaeq;
    cols [ 2 ] = lambdaeq;
}
 
 
void
LinearConstraintBC :: saveContext(DataStream &stream, ContextMode mode)
{
    if ( mode & CM_Definition ) {
        contextIOResultType iores;
        if ( ( iores = weights.storeYourself(stream) ) != CIO_OK ) {
            THROW_CIOERR(iores);
        }
        if ( ( iores = weightsTf.storeYourself(stream) ) != CIO_OK ) {
            THROW_CIOERR(iores);
        }
        if ( ( iores = dofmans.storeYourself(stream) ) != CIO_OK ) {
            THROW_CIOERR(iores);
        }
        if ( ( iores = dofs.storeYourself(stream) ) != CIO_OK ) {
            THROW_CIOERR(iores);
        }
        if ( !stream.write(rhs) ) {
            THROW_CIOERR(CIO_IOERR);
        }
        if ( !stream.write(rhsTf) ) {
            THROW_CIOERR(CIO_IOERR);
        }
        if ( ( iores = lhsType.storeYourself(stream) ) != CIO_OK ) {
            THROW_CIOERR(iores);
        }
        if ( ( iores = rhsType.storeYourself(stream) ) != CIO_OK ) {
            THROW_CIOERR(iores);
        }
    }
 
    md->saveContext(stream, mode);
}
 
 
void
LinearConstraintBC :: restoreContext(DataStream &stream, ContextMode mode)
{
    if ( mode & CM_Definition ) {
        contextIOResultType iores;
        if ( ( iores = weights.restoreYourself(stream) ) != CIO_OK ) {
            THROW_CIOERR(iores);
        }
        if ( ( iores = weightsTf.restoreYourself(stream) ) != CIO_OK ) {
            THROW_CIOERR(iores);
        }
        if ( ( iores = dofmans.restoreYourself(stream) ) != CIO_OK ) {
            THROW_CIOERR(iores);
        }
        if ( ( iores = dofs.restoreYourself(stream) ) != CIO_OK ) {
            THROW_CIOERR(iores);
        }
        if ( !stream.read(rhs) ) {
            THROW_CIOERR(CIO_IOERR);
        }
        if ( !stream.read(rhsTf) ) {
            THROW_CIOERR(CIO_IOERR);
        }
        if ( ( iores = lhsType.restoreYourself(stream) ) != CIO_OK ) {
            THROW_CIOERR(iores);
        }
        if ( ( iores = rhsType.restoreYourself(stream) ) != CIO_OK ) {
            THROW_CIOERR(iores);
        }
    }
 
    md->restoreContext(stream, mode);
}
} //end of oofem namespace