coupledfieldselement.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
 */
 
#if 1
#include "sm/Elements/coupledfieldselement.h"
#include "sm/Materials/structuralms.h"
#include "sm/CrossSections/structuralcrosssection.h"
#include "sm/Elements/nlstructuralelement.h"
#include "node.h"
#include "material.h"
#include "gausspoint.h"
#include "gaussintegrationrule.h"
#include "floatmatrix.h"
#include "floatarray.h"
#include "intarray.h"
#include "domain.h"
#include "cltypes.h"
#include "mathfem.h"
#include "unknownnumberingscheme.h"
#include <cstdio>
 
namespace oofem {
CoupledFieldsElement :: CoupledFieldsElement(int i, Domain *aDomain) : NLStructuralElement(i, aDomain)
{
    nlGeo = 0;
}
 
 
void
CoupledFieldsElement :: computeLocationArrayOfDofIDs(const IntArray &dofIdArray, IntArray &answer)
{
    // Routine to extract compute the location array an element given an dofid array.
    answer.resize(0);
 
    int k = 0;
    for ( int i = 1; i <= numberOfDofMans; i++ ) {
        DofManager *dMan = this->giveDofManager(i);        
        for (int j = 1; j <= dofIdArray.giveSize(); j++ ) {   
            if ( dMan->hasDofID( (DofIDItem) dofIdArray.at(j) ) ) {
                Dof *d = dMan->giveDofWithID( dofIdArray.at(j) );
                answer.followedBy( k + d->giveNumber() );
                //answer.followedBy( k + j );
            }
        }
        k += dMan->giveNumberOfDofs( );
    }
}
 
 
void
CoupledFieldsElement :: computeVectorOfDofIDs(const IntArray &dofIdArray, ValueModeType valueMode, TimeStep *stepN, FloatArray &answer)
{
    // Routine to extract the solution vector for an element given an dofid array.
    // Size will be numberOfDofs and if a certain dofId does not exist a zero is used as value. 
    answer.resize( numberOfDofMans * dofIdArray.giveSize() ); // equal number of nodes for all fields
    answer.zero();
    int k = 1;
    for ( int i = 1; i <= numberOfDofMans; i++ ) {
        DofManager *dMan = this->giveDofManager(i);        
        for (int j = 1; j <= dofIdArray.giveSize(); j++ ) {   
            if ( dMan->hasDofID( (DofIDItem) dofIdArray.at(j) ) ) {
                Dof *d = dMan->giveDofWithID( dofIdArray.at(j) );
                answer.at(k) = d->giveUnknown(valueMode, stepN);
            }
            k++;
        }
    }
}
 
 
 
void
CoupledFieldsElement :: giveInternalForcesVectorGen(FloatArray &answer, TimeStep *tStep, int useUpdatedGpRecord, 
    void (*Nfunc)(GaussPoint*, FloatMatrix), void (*Bfunc)(GaussPoint*, FloatMatrix, int, int), //(GaussPoint*, FloatMatrix)
    void (*NStressFunc)(GaussPoint*, FloatArray), void (*BStressFunc)(GaussPoint*, FloatArray),
    double (*dVFunc)(GaussPoint*))
{
    // General implementation of internal forces that computes
    // f = sum_gp( N^T*GenStress_N + B^T*GenStress_B ) * dV
 
    FloatArray NStress, BStress, vGenStress, NS, BS;
    FloatMatrix N, B;
 
    for ( int j = 0; j < this->giveNumberOfIntegrationRules(); j++ ) {
        for ( auto &gp: this->giveIntegrationRule(j) ) {
            double dV  = this->computeVolumeAround(gp);
            
            // compute generalized stress measures
            if ( NStressFunc && Nfunc ) {
                Nfunc(gp, N);
                NStressFunc(gp, NStress);
                NS.beTProductOf(N, NStress);
                answer.add(dV, NS);
            }
 
            if ( BStressFunc && Bfunc ) {
                Bfunc(gp, B, 1, 3);
                BStressFunc(gp, BStress);
                BS.beTProductOf(B, BStress);
                answer.add(dV, BS);
            }
        }
    }
}
 
 
void
CoupledFieldsElement :: computeStiffnessMatrixGen(FloatMatrix &answer, MatResponseMode rMode, TimeStep *tStep, 
        void (*Nfunc)(GaussPoint*, FloatMatrix), 
        void (*Bfunc)(GaussPoint*, FloatMatrix),
        void (*NStiffness)(FloatMatrix, MatResponseMode, GaussPoint*, TimeStep*), 
        void (*BStiffness)(FloatMatrix, MatResponseMode, GaussPoint*, TimeStep*),
        double (*volumeAround)(GaussPoint*) )
{
    FloatMatrix B, DB, N, DN, D_B, D_N;
 
    bool matStiffSymmFlag = this->giveCrossSection()->isCharacteristicMtrxSymmetric(rMode);
    answer.resize(0,0);
 
    for ( auto &gp: this->giveIntegrationRule(0) ) {
        double dV = this->computeVolumeAround(gp);
 
 
        // compute int_V ( N^t * D_N * N )dV
        if ( NStiffness && Nfunc ) {
            Nfunc(gp, N);
            NStiffness(D_N, rMode, gp, tStep);
            DN.beProductOf(D_N, N);
            if ( matStiffSymmFlag ) {
                answer.plusProductSymmUpper(N, DN, dV);
            } else {
                answer.plusProductUnsym(N, DN, dV);
            }
        }
 
 
        // compute int_V ( B^t * D_B * B )dV
        if ( BStiffness && Bfunc ) {
            Bfunc(gp, B);
            BStiffness(D_B, rMode, gp, tStep);
            DB.beProductOf(D_B, B);
            if ( matStiffSymmFlag ) {
                answer.plusProductSymmUpper(B, DB, dV);
            } else {
                answer.plusProductUnsym(B, DB, dV);
            }
        }
    }
 
    if ( matStiffSymmFlag ) {
        answer.symmetrized();
    }
}
 
 
} // end namespace oofem
 
#endif