tet21ghostsolid.C

Go to the documentation of this file.

/*
 *
 *                 #####    #####   ######  ######  ###   ###
 *               ##   ##  ##   ##  ##      ##      ## ### ##
 *              ##   ##  ##   ##  ####    ####    ##  #  ##
 *             ##   ##  ##   ##  ##      ##      ##     ##
 *            ##   ##  ##   ##  ##      ##      ##     ##
 *            #####    #####   ##      ######  ##     ##
 *
 *
 *             OOFEM : Object Oriented Finite Element Code
 *
 *               Copyright (C) 1993 - 2014   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/Elements/tet21ghostsolid.h"
#include "../sm/Elements/nlstructuralelement.h"
#include "../sm/CrossSections/structuralcrosssection.h"
#include "fei3dtetquad.h"
#include "fei3dtetlin.h"
#include "classfactory.h"
#include "gausspoint.h"
#include "gaussintegrationrule.h"
#include "load.h"
#include "element.h"
#include "dofmanager.h"
#include "load.h"
#include "bodyload.h"
#include "boundaryload.h"
#include "neumannmomentload.h"
#include "dof.h"
#include "mathfem.h"

#include <iostream>

#ifdef __FM_MODULE
 #include "fluiddynamicmaterial.h"
 #include "fluidcrosssection.h"
#endif

namespace oofem {
#define FIXEDJ 0
#define FIXEDFinvT 0
#define FIXEDpressure 0
#define USEUNCOUPLED 1
#define USENUMTAN 0
#define TESTTANGENT 0

#define VELOCITYCOEFF 1.0

REGISTER_Element(tet21ghostsolid);

FEI3dTetQuad tet21ghostsolid :: interpolation;
FEI3dTetLin tet21ghostsolid :: interpolation_lin;

IntArray tet21ghostsolid :: momentum_ordering(30);
IntArray tet21ghostsolid :: conservation_ordering(4);
IntArray tet21ghostsolid :: ghostdisplacement_ordering(30);
IntArray tet21ghostsolid :: velocitydofsonside = {
    1, 2, 3, 8, 9, 10, 15, 16, 17, 22, 23, 24, 28, 29, 30, 34, 35, 36
};
IntArray tet21ghostsolid :: displacementdofsonside = {
    4, 5, 6, 11, 12, 13, 18, 19, 20, 25, 26, 27, 31, 32, 33, 37, 38, 39
};

tet21ghostsolid :: tet21ghostsolid(int n, Domain *aDomain) : NLStructuralElement(n, aDomain), SpatialLocalizerInterface(this)
{
    numberOfGaussPoints = 4;
    numberOfDofMans = 10;
    computeItransform = true;

    double nu = .25, E = 1;
    Dghost.resize(6, 6);
    Dghost.zero();
    Dghost.at(1, 1) = 1 - nu;
    Dghost.at(1, 2) = Dghost.at(1, 3) = nu;
    Dghost.at(2, 2) = 1 - nu;
    Dghost.at(2, 1) = Dghost.at(2, 3) = nu;
    Dghost.at(3, 3) = 1 - nu;
    Dghost.at(3, 1) = Dghost.at(3, 2) = nu;
    Dghost.at(4, 4) = Dghost.at(5, 5) = Dghost.at(6, 6) = .5 * ( 1 - 2 * nu );
    Dghost.times( E / ( 1 + nu ) / ( 1 - 2 * nu ) );

    conservation_ordering = {
        7, 14, 21, 28
    };

    for ( int i = 0, j = 1; i < 10; ++i ) {
        momentum_ordering(i * 3 + 0) = j++;
        momentum_ordering(i * 3 + 1) = j++;
        momentum_ordering(i * 3 + 2) = j++;
        if ( i <= 3 ) {
            j++;
        }
        j = j + 3;
    }

    for ( int i = 0, j = 1; i < 10; ++i ) {
        j = j + 3;
        ghostdisplacement_ordering(i * 3 + 0) = j++;
        ghostdisplacement_ordering(i * 3 + 1) = j++;
        ghostdisplacement_ordering(i * 3 + 2) = j++;
        if ( i <= 3 ) {
            j++;
        }
    }
}

FEInterpolation *
tet21ghostsolid :: giveInterpolation() const
{
    return & interpolation;
}

FEInterpolation *
tet21ghostsolid :: giveInterpolation(DofIDItem id) const
{
    if ( id == P_f ) {
        return & interpolation_lin;
    } else {
        return & interpolation;
    }
}

void
tet21ghostsolid :: computeGaussPoints()
// Sets up the array containing the four Gauss points of the receiver.
{
    if ( integrationRulesArray.size() == 0 ) {
        integrationRulesArray.resize(1);
        integrationRulesArray [ 0 ].reset( new GaussIntegrationRule(1, this, 1, 6) );
        this->giveCrossSection()->setupIntegrationPoints(* integrationRulesArray [ 0 ], numberOfGaussPoints, this);
    }
}

void
tet21ghostsolid :: computeNumericStiffnessMatrix(FloatMatrix &answer, MatResponseMode rMode, TimeStep *tStep)
{
    FloatArray a, aPert, intF, intFPert, DintF;
    double eps = 1e-9;

    answer.resize(64, 64);
    answer.zero();

    this->computeVectorOf(VM_Total, tStep, a);

    giveInternalForcesVectorGivenSolution(intF, tStep, 0, a);

    for ( int i = 1; i <= answer.giveNumberOfColumns(); i++ ) {
        aPert = a;
        aPert.at(i) = aPert.at(i) + eps;

        giveInternalForcesVectorGivenSolution(intFPert, tStep, 0, aPert);

        DintF.operator = ( intF - intFPert );
        DintF.times(-1 / eps);
        answer.setColumn(DintF, i);
    }
}

void
tet21ghostsolid :: computeNumericStiffnessMatrixDebug(FloatMatrix &answer, MatResponseMode rMode, TimeStep *tStep)
{
    FloatArray a, aPert, intF, intFPert, DintF;
    double eps = 1e-9;

    answer.resize(64, 64);
    answer.zero();

    this->computeVectorOf(VM_Total, tStep, a);

    giveInternalForcesVectorGivenSolutionDebug(intF, tStep, 0, a, true);

    for ( int i = 1; i <= answer.giveNumberOfColumns(); i++ ) {
        aPert = a;
        aPert.at(i) = aPert.at(i) + eps;

        giveInternalForcesVectorGivenSolutionDebug(intFPert, tStep, 0, aPert, false);

        DintF.operator = ( intF - intFPert );
        DintF.times(-1 / eps);
        answer.setColumn(DintF, i);
    }
}


void
tet21ghostsolid :: computeStressVector(FloatArray &answer, const FloatArray &strain, GaussPoint *gp, TimeStep *tStep)
{
    this->giveStructuralCrossSection()->giveRealStress_3d(answer, gp, strain, tStep);
}

void
tet21ghostsolid :: computeConstitutiveMatrixAt(FloatMatrix &answer, MatResponseMode rMode, GaussPoint *gp, TimeStep *tStep)
{
    this->giveStructuralCrossSection()->giveStiffnessMatrix_3d(answer, rMode, gp, tStep);
}

void
tet21ghostsolid :: computeStiffnessMatrix(FloatMatrix &answer, MatResponseMode rMode, TimeStep *tStep)
{
#if USENUMTAN == 1
    computeNumericStiffnessMatrix(answer, rMode, tStep);
    return;

#endif

#ifdef __FM_MODULE
    FluidDynamicMaterial *fluidMaterial = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial();
#endif

    FloatMatrix afDu, afDw, bfuDu, bfuDp, bfpDu, bfpDw, cfwDu;
    FloatMatrix Kf, G, Kx, Ed, EdB, dNx;
    FloatArray Nlin, dNv, a, a_prev, a_inc;
    FloatArray aVelocity, aPressure, aGhostDisplacement, aIncGhostDisplacement;

    afDu.resize(30, 30);
    afDu.zero();
    bfuDu.resize(30, 30); // Remove?
    bfuDu.zero();
    bfuDp.resize(30, 4);
    bfuDp.zero();
    bfpDu.resize(4, 30);
    bfpDu.zero();
    bfpDw.resize(4, 30);
    bfpDw.zero();
    cfwDu.resize(30, 30);
    cfwDu.zero();

    this->computeVectorOf(VM_Total, tStep, a);
    if ( !tStep->isTheFirstStep() ) {
        this->computeVectorOf(VM_Total, tStep->givePreviousStep(), a_prev);
    } else {
        a_prev.resize( a.giveSize() );
        a_prev.zero();
    }
    a_inc.operator = ( a - a_prev );

    aVelocity.beSubArrayOf(a, momentum_ordering);
    aPressure.beSubArrayOf(a, conservation_ordering);
    aGhostDisplacement.beSubArrayOf(a, ghostdisplacement_ordering);
    aIncGhostDisplacement.beSubArrayOf(a_inc, ghostdisplacement_ordering);

    FloatArray aTotal;
    aTotal.operator = ( aVelocity + aIncGhostDisplacement * VELOCITYCOEFF ); // Assume deltaT=1 gives that the increment is the velocity

    for ( auto &gp : *this->giveDefaultIntegrationRulePtr() ) {
        // gp = this->giveDefaultIntegrationRulePtr()->getIntegrationPoint(0);
        const FloatArray &lcoords = gp->giveNaturalCoordinates();
        double detJ = fabs( ( this->interpolation.giveTransformationJacobian( lcoords, FEIElementGeometryWrapper(this) ) ) );
        double weight = gp->giveWeight();

        this->interpolation.evaldNdx( dNx, lcoords, FEIElementGeometryWrapper(this) );
        this->interpolation_lin.evalN( Nlin, lcoords, FEIElementGeometryWrapper(this) );

        // Move to small deformations
        dNv.resize(30); // dNv = [dN1/dx dN1/dy dN1/dz dN2/dx dN2/dy dN2/dz ... dN10/dz]
        for ( int k = 0; k < dNx.giveNumberOfRows(); k++ ) {
            dNv.at(k * 3 + 1) = dNx.at(k + 1, 1);
            dNv.at(k * 3 + 2) = dNx.at(k + 1, 2);
            dNv.at(k * 3 + 3) = dNx.at(k + 1, 3);
        }

#ifdef __FM_MODULE
        gp->setMaterialMode(_3dFlow);
        fluidMaterial->giveDeviatoricStiffnessMatrix(Ed, TangentStiffness, gp, tStep);
        gp->setMaterialMode(_3dMat);
#else
        OOFEM_ERROR("Fluid module missing\n");
#endif

        if ( nlGeometry == 0 ) {
            FloatMatrix B;
            this->computeBmatrixAt(gp, B);

            // Fluid part
            EdB.beProductOf(Ed, B);
            Kf.plusProductSymmUpper(B, EdB, detJ * weight);
            Kf.symmetrized();

            // Ghost solid part
            Kx.plusProductSymmUpper(B, EdB, detJ * weight);

            // Incompressibility part
            G.plusDyadUnsym(dNv, Nlin, -detJ * weight);
        } else {
            int Voigt6 [ 3 ] [ 3 ] = { { 1, 6, 5 }, { 6, 2, 4 }, { 5, 4, 3 } };
            int Voigt9 [ 3 ] [ 3 ] = { { 1, 6, 5 }, { 9, 2, 4 }, { 8, 7, 3 } };

            // Compute commonly used matrices
            FloatArray Fa, Finva, FinvTa;
            FloatMatrix F, Finv, FinvT, BH, B;

            this->computeBHmatrixAt(gp, BH);
            this->computeBmatrixAt(gp, B);

            computeDeformationGradientVectorFromDispl(Fa, gp, tStep, aGhostDisplacement);
            F.beMatrixForm(Fa);
            Finv.beInverseOf(F);
            Finva.beVectorForm(Finv);
            FinvT.beTranspositionOf(Finv);
            FinvTa.beVectorForm(FinvT);

            FloatArray FinvTaB;
            FinvTaB.beTProductOf(BH, FinvTa);

            FloatMatrix NlinFinvTaB;
            NlinFinvTaB.beDyadicProductOf(Nlin, FinvTaB);

            double J = F.giveDeterminant();
#if FIXEDJ == 1
            J = 1;
#endif

            // Compute Cauchy stress
            FloatArray fluidCauchyArray, epsf;
            FloatMatrix fluidCauchyMatrix;
            double epsvol, pressure;

            epsf.beProductOf(B, aTotal);
            pressure = Nlin.dotProduct(aPressure);

#if FIXEDpressure == 1
            pressure = 1.0;
#endif

            gp->setMaterialMode(_3dFlow);
            fluidMaterial->computeDeviatoricStressVector(fluidCauchyArray, epsvol, gp, epsf, pressure, tStep);
            gp->setMaterialMode(_3dMat);
            fluidCauchyMatrix.beMatrixFormOfStress(fluidCauchyArray);

            // afDu terms ==========
            // 1st term ------------
            FloatArray fluidTwoPointStressArray, BTS;
            FloatMatrix fluidTwoPointStressMatrix, afuT1;

            fluidTwoPointStressMatrix.beProductOf(fluidCauchyMatrix, FinvT);
            fluidTwoPointStressArray.beVectorForm(fluidTwoPointStressMatrix);

            BTS.beTProductOf(BH, fluidTwoPointStressArray);

            afuT1.beDyadicProductOf(BTS, FinvTaB);
            afuT1.times(J * weight * detJ);
#if FIXEDJ == 1
            afuT1.times(0.0);
#endif
            afDu.add(afuT1);

            // 2nd term ------------ // Works
            FloatMatrix afuT2;
            FloatMatrix G, BTG;

            G.resize(9, 9);
            G.zero();

            for ( int i = 1; i <= 3; i++ ) {
                for ( int j = 1; j <= 3; j++ ) {
                    for ( int k = 1; k <= 3; k++ ) {
                        for ( int l = 1; l <= 3; l++ ) {
                            for ( int m = 1; m <= 3; m++ ) {
                                G.at(Voigt9 [ i - 1 ] [ j - 1 ], Voigt9 [ k - 1 ] [ l - 1 ]) +=
                                    Ed.at(Voigt6 [ i - 1 ] [ m - 1 ], Voigt6 [ k - 1 ] [ l - 1 ]) * FinvT.at(m, j);
                            }
                        }
                    }
                }
            }

            BTG.beTProductOf(BH, G);
            afuT2.beProductOf(BTG, BH);
            afuT2.times(J * weight * detJ);

            afDu.add(afuT2);

            // 3rd term ------------
            FloatMatrix afuT3;
            FloatMatrix K, BTK;

            K.resize(9, 9);
            K.zero();

            for ( int i = 1; i <= 3; i++ ) {
                for ( int j = 1; j <= 3; j++ ) {
                    for ( int k = 1; k <= 3; k++ ) {
                        for ( int l = 1; l <= 3; l++ ) {
                            for ( int m = 1; m <= 3; m++ ) {
                                K.at(Voigt9 [ i - 1 ] [ j - 1 ], Voigt9 [ m - 1 ] [ l - 1 ]) +=
                                    fluidCauchyArray.at(Voigt6 [ i - 1 ] [ k - 1 ]) * FinvT.at(k, l) * FinvT.at(m, j);
                            }
                        }
                    }
                }
            }

            BTK.beTProductOf(BH, K);
            afuT3.beProductOf(BTK, BH);
            afuT3.times(-J * weight * detJ);
#if FIXEDFinvT == 1
            afuT3.times(0.0);
#endif
            afDu.add(afuT3);

            // afDw terms ========== // Works
            afDw.add(afuT2);

            // bfuDu term ==========
            // 1st term ------------ // Works
            FloatMatrix bfuT1;
            bfuT1.beDyadicProductOf(FinvTaB, FinvTaB);
            bfuT1.times(-J * pressure * weight * detJ);

            bfuDu.add(bfuT1);

            // 2nd term ------------ // Works
            FloatMatrix C, bfuT2, BTC;
            C.resize(9, 9);
            C.zero();

            for ( int i = 1; i <= 3; i++ ) {
                for ( int j = 1; j <= 3; j++ ) {
                    for ( int k = 1; k <= 3; k++ ) {
                        for ( int l = 1; l <= 3; l++ ) {
                            C.at(Voigt9 [ i - 1 ] [ j - 1 ], Voigt9 [ k - 1 ] [ l - 1 ]) +=
                                FinvT.at(i, l) * FinvT.at(k, j);
                        }
                    }
                }
            }

            BTC.beTProductOf(BH, C);
            bfuT2.beProductOf(BTC, BH);
            bfuT2.times(pressure * J * weight * detJ);

            bfuDu.add(bfuT2);

            // bfuDp term ========== // Works
            FloatMatrix bfuT3;
            bfuT3.beDyadicProductOf(FinvTaB, Nlin);
            bfuT3.times(-J * detJ * weight);

            bfuDp.add(bfuT3);

            // bfpDu term ==========
            // 1st term ------------
            FloatMatrix bfpT1, A_IVm, Bvm;
            FloatArray A_IVa, Bv;

            Bv.beProductOf(BH, aTotal);
            Bvm.beMatrixForm(Bv);
            double FinvTaBv = FinvTa.dotProduct(Bv);

            A_IVm.resize(3, 3);
            A_IVm.zero();

            for ( int k = 1; k <= 3; k++ ) {
                for ( int l = 1; l <= 3; l++ ) {
                    for ( int i = 1; i <= 3; i++ ) {
                        for ( int j = 1; j <= 3; j++ ) {
                            A_IVm.at(k, l) += Bvm.at(i, j) *    FinvT.at(i, l) * FinvT.at(k, j);
                        }
                    }
                }
            }

            A_IVa.beVectorForm(A_IVm);

            bfpT1 = NlinFinvTaB;
            bfpT1.times(-J * detJ * weight * FinvTaBv);
#if FIXEDJ == 1
            bfpT1.times(0.0);
#endif

            bfpDu.add(bfpT1);

            // 2nd term ------------
            FloatMatrix bfpT2;
            FloatArray AiiB;

            AiiB.beTProductOf(BH, A_IVa);
            bfpT2.beDyadicProductOf(Nlin, AiiB);
            bfpT2.times(J * detJ * weight);
#if FIXEDFinvT == 1
            bfpT2.times(0.0);
#endif
            bfpDu.add(bfpT2);


            // 3rd term ------------ // Works
            FloatMatrix bfpT3;
            bfpT3 = NlinFinvTaB;
            bfpT3.times(-J * detJ * weight);
            bfpDu.add(bfpT3);

            // bfpDw term ========== // Works
            //            bfpT1 = NlinFinvTaB;
            //            bfpT1.times(-J*detJ*weight);
            bfpDw.add(bfpT3);

            // cfwDu term ========== // Works
            FloatMatrix BTEX, cfwT1;
            BTEX.beTProductOf(B, Dghost);
            cfwT1.beProductOf(BTEX, B);
            cfwT1.times(detJ * weight);

            cfwDu.add(cfwT1);
        }
    }

    //if (this->globalNumber == 373) { dJda.printYourself("dJda_analytical"); }
    FloatMatrix GT, temp;

    GT.beTranspositionOf(G);
    Kf.symmetrized();
    Kx.symmetrized();

    answer.resize(64, 64);
    answer.zero();
    temp.resize(64, 64);
    temp.zero();

    temp.assemble(afDu, ghostdisplacement_ordering, ghostdisplacement_ordering);
    temp.assemble(afDw, ghostdisplacement_ordering, momentum_ordering);
    temp.assemble(bfuDu, ghostdisplacement_ordering, ghostdisplacement_ordering);
    temp.assemble(bfuDp, ghostdisplacement_ordering, conservation_ordering);
    temp.assemble(bfpDu, conservation_ordering, ghostdisplacement_ordering);
    temp.assemble(bfpDw, conservation_ordering, momentum_ordering);
    temp.assemble(cfwDu, momentum_ordering, ghostdisplacement_ordering);

    if ( this->computeItransform ) {
        giveRowTransformationMatrix(tStep);
    }

    answer.beProductOf(Itransform, temp);

    // ******************* This is for development purposes only. Compare to numerical stiffness tangent
#if TESTTANGENT
    if ( ( this->globalNumber == 292 ) && ( tStep->giveNumber() > 0 ) ) {
        double MaxErr = 0.0;
        FloatMatrix numtan;
        int imax, jmax;

        temp.zero();
        temp.assemble(afDu, ghostdisplacement_ordering, ghostdisplacement_ordering);
        temp.assemble(afDw, ghostdisplacement_ordering, momentum_ordering);
        temp.assemble(bfuDu, ghostdisplacement_ordering, ghostdisplacement_ordering);
        temp.assemble(bfuDp, ghostdisplacement_ordering, conservation_ordering);
        temp.assemble(bfpDu, conservation_ordering, ghostdisplacement_ordering);
        temp.assemble(bfpDw, conservation_ordering, momentum_ordering);
        temp.assemble(cfwDu, momentum_ordering, ghostdisplacement_ordering);

        computeNumericStiffnessMatrixDebug(numtan, rMode, tStep);
        for ( int i = 1; i <= numtan.giveNumberOfRows(); i++ ) {
            for ( int j = 1; j <= numtan.giveNumberOfRows(); j++ ) {
                double difference = numtan.at(i, j) - temp.at(i, j);
                if ( fabs(difference) > fabs(MaxErr) ) {
                    MaxErr = difference;
                    imax = i;
                    jmax = j;
                }
                if ( fabs(numtan.at(i, j) < 1e-15) ) {
                    numtan.at(i, j) = 0.0;
                }
                if ( fabs(temp.at(i, j) < 1e-15) ) {
                    temp.at(i, j) = 0.0;
                }
            }
        }

        printf( "globalNumber: %u, i:%u, j:%u, MaxErr: %e, Numtan=%f, Analytic=%f\n", this->globalNumber, imax, jmax, MaxErr, numtan.at(imax, jmax), temp.at(imax, jmax) );

        numtan.printYourselfToFile("/tmp/numtan.txt", false);
        temp.printYourselfToFile("/tmp/analytic.txt", false);
        if ( fabs(MaxErr) > 1e-5 ) {
            printf("The error is actually quite large..\n");
        }
    }
#endif
    // *******************
    //computeNumericStiffnessMatrix(answer, rMode, tStep);
}

void
tet21ghostsolid :: computeLoadVector(FloatArray &answer, BodyLoad *load, CharType type, ValueModeType mode, TimeStep *tStep)
{
    // Compute displacements used to compute J

    answer.resize(64);
    answer.zero();

    if ( type != ExternalForcesVector ) {
        answer.resize(0);
        return;
    }

#ifdef __FM_MODULE
    FluidDynamicMaterial *mat = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial();
#endif
    FloatArray N, gVector, temparray(30), dNv, inc, a, a_prev, u, u_prev, vload;
    FloatMatrix dNx, G;

    load->computeComponentArrayAt(gVector, tStep, VM_Total);
    temparray.zero();

    vload.resize(4);
    vload.zero();

    this->giveUnknownData(a_prev, a, inc, tStep);
    u.beSubArrayOf(a, ghostdisplacement_ordering);
    u_prev.beSubArrayOf(a_prev, ghostdisplacement_ordering);

    for ( auto &gp : *this->integrationRulesArray [ 0 ] ) {
        const FloatArray &lcoords = gp->giveNaturalCoordinates();

        double detJ = fabs( this->interpolation.giveTransformationJacobian( lcoords, FEIElementGeometryWrapper(this) ) );
        double dA = detJ * gp->giveWeight();

        // Body load
        if ( gVector.giveSize() ) {
#ifdef __FM_MODULE
            double rho = mat->give('d', gp);
#else
            OOFEM_ERROR("Missing FM module");
            double rho = 1.0;
#endif

            if ( this->nlGeometry ) {
                FloatArray Fa, temp;
                FloatMatrix F, Finv, FinvT;
                computeDeformationGradientVectorFromDispl(Fa, gp, tStep, u);
                F.beMatrixForm(Fa);
                Finv.beInverseOf(F);
                FinvT.beTranspositionOf(Finv);

                dA = dA * F.giveDeterminant();

                temp.beProductOf(Finv, gVector);
                gVector = temp;
            }

            this->interpolation.evalN( N, lcoords, FEIElementGeometryWrapper(this) );

            for ( int j = 0; j < N.giveSize(); j++ ) {
                temparray(3 * j + 0) += N(j) * rho * gVector(0) * dA;
                temparray(3 * j + 1) += N(j) * rho * gVector(1) * dA;
                temparray(3 * j + 2) += N(j) * rho * gVector(2) * dA;
            }
        }

        // "load" from previous step
        this->interpolation.evaldNdx( dNx, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );
        this->interpolation_lin.evalN( N, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );

        dNv.resize(30);
        for ( int k = 0; k < dNx.giveNumberOfRows(); k++ ) {
            dNv.at(k * 3 + 1) = dNx.at(k + 1, 1);
            dNv.at(k * 3 + 2) = dNx.at(k + 1, 2);
            dNv.at(k * 3 + 3) = dNx.at(k + 1, 3);
        }

        G.plusDyadUnsym(N, dNv, dA);
        FloatMatrix GT;
        GT.beTranspositionOf(G);

        vload.plusProduct(GT, u_prev, -0.0);
        //vload.printYourself();
    }

    if ( this->computeItransform ) {
        giveRowTransformationMatrix(tStep);
    }

    FloatArray temp;

    temp.resize(64);
    temp.zero();
    temp.assemble(temparray, this->ghostdisplacement_ordering);
    temp.assemble(vload, this->conservation_ordering);

    answer.beProductOf(Itransform, temp);
}

void
tet21ghostsolid :: giveInternalForcesVector(FloatArray &answer, TimeStep *tStep, int useUpdatedGpRecord)
{
    FloatArray a;
    this->computeVectorOf(VM_Total, tStep, a);
    giveInternalForcesVectorGivenSolution(answer, tStep, useUpdatedGpRecord, a);
}

void
tet21ghostsolid :: giveInternalForcesVectorGivenSolution(FloatArray &answer, TimeStep *tStep, int useUpdatedGpRecord, FloatArray &a)
{
#ifdef __FM_MODULE
    FluidDynamicMaterial *fluidMaterial = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial();
#endif

    FloatMatrix Kf, G, Kx, B, Ed, dNx;
    FloatArray Strain, Stress, Nlin, dNv, a_inc, a_prev, aVelocity, aPressure, aGhostDisplacement, aIncGhostDisplacement, fluidStress, epsf, dpdivv;
    FloatArray momentum, conservation, auxstress, divv;
    double pressure, epsvol = 0.0;

    if ( !tStep->isTheFirstStep() ) {
        this->computeVectorOf(VM_Total, tStep->givePreviousStep(), a_prev);
    } else {
        a_prev.resize( a.giveSize() );
        a_prev.zero();
    }
    a_inc.operator = ( a - a_prev );

    aVelocity.beSubArrayOf(a, momentum_ordering);
    aPressure.beSubArrayOf(a, conservation_ordering);
    aGhostDisplacement.beSubArrayOf(a, ghostdisplacement_ordering);
    aIncGhostDisplacement.beSubArrayOf(a_inc, ghostdisplacement_ordering);

    FloatArray aTotal;
    aTotal.operator = ( aVelocity + aIncGhostDisplacement * VELOCITYCOEFF ); // Assume deltaT=1 gives that the increment is the velocity

    for ( auto &gp : *this->giveDefaultIntegrationRulePtr() ) {
        const FloatArray &lcoords = gp->giveNaturalCoordinates();
        double detJ = fabs( ( this->interpolation.giveTransformationJacobian( lcoords, FEIElementGeometryWrapper(this) ) ) );
        double weight = gp->giveWeight();

        this->interpolation.evaldNdx( dNx, lcoords, FEIElementGeometryWrapper(this) );
        this->interpolation_lin.evalN( Nlin, lcoords, FEIElementGeometryWrapper(this) );

        dNv.resize(30);
        for ( int k = 0; k < dNx.giveNumberOfRows(); k++ ) {
            dNv.at(k * 3 + 1) = dNx.at(k + 1, 1);
            dNv.at(k * 3 + 2) = dNx.at(k + 1, 2);
            dNv.at(k * 3 + 3) = dNx.at(k + 1, 3);
        }

        pressure = Nlin.dotProduct(aPressure);

        if ( nlGeometry == 0 ) {
            this->computeBmatrixAt(gp, B);

            epsf.beProductOf(B, aTotal);

            // Momentum equation
            gp->setMaterialMode(_3dFlow);
#ifdef __FM_MODULE
            fluidMaterial->computeDeviatoricStressVector(fluidStress, epsvol, gp, epsf, pressure, tStep);
#else
            OOFEM_ERROR("Missing FM module");
#endif
            gp->setMaterialMode(_3dMat);

            momentum.plusProduct(B, fluidStress, detJ * weight);
            momentum.add(-pressure * detJ * weight, dNv);

            // Conservation equation
            divv.beTProductOf(dNv, aTotal);
            divv.times(-detJ * weight);

            conservation.add(divv.at(1), Nlin);

            // Ghost solid part
            Strain.beProductOf(B, aGhostDisplacement);
            Stress.beProductOf(Dghost, Strain);
            auxstress.plusProduct(B, Stress, detJ * weight);
        } else {
            FloatMatrix BH;

            this->computeBHmatrixAt(gp, BH);
            this->computeBmatrixAt(gp, B);

            // Compute deformation gradient etc.
            FloatArray Fa, FinvTa;
            FloatMatrix F, Finv, FinvT, fluidStressMatrix;

            computeDeformationGradientVectorFromDispl(Fa, gp, tStep, aGhostDisplacement);
            F.beMatrixForm(Fa);
            double J = F.giveDeterminant();

            Finv.beInverseOf(F);
            FinvT.beTranspositionOf(Finv);
            FinvTa.beVectorForm(FinvT);

            epsf.beProductOf(B, aTotal);

            // Momentum equation --------------------------------------------------

            // Compute fluid cauchy stress
            FloatArray fluidCauchyArray;
            FloatMatrix fluidCauchyMatrix;

            gp->setMaterialMode(_3dFlow);
#ifdef __FM_MODULE
            fluidMaterial->computeDeviatoricStressVector(fluidCauchyArray, epsvol, gp, epsf, pressure, tStep);
#else
            OOFEM_ERROR("Missing FM module");
#endif
            gp->setMaterialMode(_3dMat);

            // Transform to 1st Piola-Kirshhoff
            fluidCauchyMatrix.beMatrixFormOfStress(fluidCauchyArray);
            fluidStressMatrix.beProductOf(fluidCauchyMatrix, FinvT);
            fluidStress.beVectorForm(fluidStressMatrix);

            // Add to equation
            momentum.plusProduct(BH, fluidStress, detJ * J * weight);

            // Pressure term in momentum equation
            FloatArray ptemp;
            ptemp.beTProductOf(BH, FinvTa);

            momentum.add(-pressure * detJ * weight * J, ptemp);

            // Conservation equation ---------------------------------------------
            FloatArray BvaTotal;
            double FinvTaBvaTotal;

            BvaTotal.beProductOf(BH, aTotal);
            FinvTaBvaTotal = FinvTa.dotProduct(BvaTotal);
            conservation.add(-J * detJ * weight * FinvTaBvaTotal, Nlin);

            // Ghost solid part --------------------------------------------------
            Strain.beProductOf(B, aGhostDisplacement);
            Stress.beProductOf(Dghost, Strain);
            auxstress.plusProduct(B, Stress, detJ * weight);
        }
    }
    answer.resize(64);
    answer.zero();

    FloatArray temp;

    temp.resize(64);
    temp.zero();

    temp.assemble(momentum, ghostdisplacement_ordering);
    temp.assemble(conservation, conservation_ordering);
    temp.assemble(auxstress, momentum_ordering);

    if ( this->computeItransform ) {
        giveRowTransformationMatrix(tStep);
    }

    answer.beProductOf(Itransform, temp);
}

void
tet21ghostsolid :: giveInternalForcesVectorGivenSolutionDebug(FloatArray &answer, TimeStep *tStep, int useUpdatedGpRecord, FloatArray &a, bool ExtraLogging)
{
#ifdef __FM_MODULE
    FluidDynamicMaterial *fluidMaterial = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial();
#endif

    FloatMatrix Kf, G, Kx, B, Ed, dNx;
    FloatArray Strain, Stress, Nlin, dNv, a_inc, a_prev, aVelocity, aPressure, aGhostDisplacement, aIncGhostDisplacement, fluidStress, epsf, dpdivv;
    FloatArray momentum, conservation, auxstress, divv;
    double pressure, epsvol = 0.0;

    if ( !tStep->isTheFirstStep() ) {
        this->computeVectorOf(VM_Total, tStep->givePreviousStep(), a_prev);
    } else {
        a_prev.resize( a.giveSize() );
        a_prev.zero();
    }
    a_inc.operator = ( a - a_prev );

    aVelocity.beSubArrayOf(a, momentum_ordering);
    aPressure.beSubArrayOf(a, conservation_ordering);
    aGhostDisplacement.beSubArrayOf(a, ghostdisplacement_ordering);
    aIncGhostDisplacement.beSubArrayOf(a_inc, ghostdisplacement_ordering);

    FloatArray aTotal;
    aTotal.operator = ( aVelocity + aIncGhostDisplacement * VELOCITYCOEFF ); // Assume deltaT=1 gives that the increment is the velocity

    for ( auto &gp : *this->giveDefaultIntegrationRulePtr() ) {
        // gp = this->giveDefaultIntegrationRulePtr()->getIntegrationPoint(0);
        const FloatArray &lcoords = gp->giveNaturalCoordinates();
        double detJ = fabs( ( this->interpolation.giveTransformationJacobian( lcoords, FEIElementGeometryWrapper(this) ) ) );
        double weight = gp->giveWeight();

        this->interpolation.evaldNdx( dNx, lcoords, FEIElementGeometryWrapper(this) );
        this->interpolation_lin.evalN( Nlin, lcoords, FEIElementGeometryWrapper(this) );

        dNv.resize(30);
        for ( int k = 0; k < dNx.giveNumberOfRows(); k++ ) {
            dNv.at(k * 3 + 1) = dNx.at(k + 1, 1);
            dNv.at(k * 3 + 2) = dNx.at(k + 1, 2);
            dNv.at(k * 3 + 3) = dNx.at(k + 1, 3);
        }

        pressure = Nlin.dotProduct(aPressure);
#if FIXEDpressure == 1
        pressure = 1.0;
#endif

        if ( nlGeometry == 0 ) {
            this->computeBmatrixAt(gp, B);

            epsf.beProductOf(B, aTotal);

            // Momentum equation
            gp->setMaterialMode(_3dFlow);
#ifdef __FM_MODULE
            fluidMaterial->computeDeviatoricStressVector(fluidStress, epsvol, gp, epsf, pressure, tStep);
#else
            OOFEM_ERROR("Missing FM module");
#endif
            gp->setMaterialMode(_3dMat);

            momentum.plusProduct(B, fluidStress, detJ * weight);
            momentum.add(-pressure * detJ * weight, dNv);

            // Conservation equation
            divv.beTProductOf(dNv, aTotal);
            divv.times(-detJ * weight);

            conservation.add(divv.at(1), Nlin);

            // Ghost solid part
            Strain.beProductOf(B, aGhostDisplacement);
            Stress.beProductOf(Dghost, Strain);
            auxstress.plusProduct(B, Stress, detJ * weight);
        } else {
            FloatMatrix BH;

            this->computeBHmatrixAt(gp, BH);
            this->computeBmatrixAt(gp, B);

            // Compute deformation gradient etc.
            FloatArray Fa, FinvTa;
            FloatMatrix F, Finv, FinvT, fluidStressMatrix;

            computeDeformationGradientVectorFromDispl(Fa, gp, tStep, aGhostDisplacement);
            F.beMatrixForm(Fa);
            double J = F.giveDeterminant();
#if FIXEDJ == 1
            J = 1;
#endif

#if FIXEDFinvT == 1
            FloatArray aTrue, dispTrue;
            this->computeVectorOf(VM_Total, tStep, aTrue);
            dispTrue.beSubArrayOf(aTrue, ghostdisplacement_ordering);
            computeDeformationGradientVectorFromDispl(Fa, gp, tStep, dispTrue);
            F.beMatrixForm(Fa);
#endif

            Finv.beInverseOf(F);
            FinvT.beTranspositionOf(Finv);
            FinvTa.beVectorForm(FinvT);

            epsf.beProductOf(B, aTotal);

            // Momentum equation --------------------------------------------------

            // Compute fluid cauchy stress
            FloatArray fluidCauchy;
            FloatMatrix fluidCauchyMatrix;

            gp->setMaterialMode(_3dFlow);
#ifdef __FM_MODULE
            fluidMaterial->computeDeviatoricStressVector(fluidCauchy, epsvol, gp, epsf, pressure, tStep);
#else
            OOFEM_ERROR("Missing FM module");
#endif
            gp->setMaterialMode(_3dMat);

            // Transform to 1st Piola-Kirshhoff
            fluidCauchyMatrix.beMatrixFormOfStress(fluidCauchy);
            fluidStressMatrix.beProductOf(fluidCauchyMatrix, FinvT);
            fluidStress.beVectorForm(fluidStressMatrix);

            /*            if ( (this->globalNumber==292) && (ExtraLogging)) {
             *              std :: cout << "gp@" << gp << " ";
             *              F.printYourself("F in NUMERIC");
             *              epsf.printYourself("epsf in NUMERIC");
             *              fluidCauchy.printYourself("fluidCauchy in NUMERIC");
             *          } */

            // Add to equation
            momentum.plusProduct(BH, fluidStress, detJ * J * weight);

            // Pressure term in momentum equation
            FloatArray ptemp;
            ptemp.beTProductOf(BH, FinvTa);
            momentum.add(-pressure * detJ * weight * J, ptemp);

            // Conservation equation ---------------------------------------------
            FloatArray BvaTotal;
            double FinvTaBvaTotal;

            BvaTotal.beProductOf(BH, aTotal);
            FinvTaBvaTotal = FinvTa.dotProduct(BvaTotal);
            conservation.add(-J * detJ * weight * FinvTaBvaTotal, Nlin);

            // Ghost solid part --------------------------------------------------
            Strain.beProductOf(B, aGhostDisplacement);
            Stress.beProductOf(Dghost, Strain);
            auxstress.plusProduct(B, Stress, detJ * weight);
        }
    }
    answer.resize(64);
    answer.zero();

    FloatArray temp;

    temp.resize(64);
    temp.zero();

    temp.assemble(momentum, ghostdisplacement_ordering);
    temp.assemble(conservation, conservation_ordering);
    temp.assemble(auxstress, momentum_ordering);

    if ( this->computeItransform ) {
        giveRowTransformationMatrix(tStep);
    }

    answer.beProductOf(Itransform, temp);
} // **************** DEBUG

void
tet21ghostsolid :: giveDofManDofIDMask(int inode, IntArray &answer) const
{
    // Returns the mask for node number inode of this element.

    if ( inode <= 4 ) {
        answer = {
            V_u, V_v, V_w, D_u, D_v, D_w, P_f
        };
    } else {
        answer = {
            V_u, V_v, V_w, D_u, D_v, D_w
        };
    }
}

void
tet21ghostsolid :: computeBmatrixAt(GaussPoint *gp, FloatMatrix &answer, int li, int ui)
// Returns the [6x30] strain-displacement matrix {B} of the receiver, eva-
// luated at gp.
// B matrix  -  6 rows : epsilon-X, epsilon-Y, epsilon-Z, gamma-YZ, gamma-ZX, gamma-XY  :
{
    FloatMatrix dnx;

    this->interpolation.evaldNdx( dnx, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );

    answer.resize(6, 30);
    answer.zero();

    for ( int i = 1; i <= 10; i++ ) {
        answer.at(1, 3 * i - 2) = dnx.at(i, 1);
        answer.at(2, 3 * i - 1) = dnx.at(i, 2);
        answer.at(3, 3 * i - 0) = dnx.at(i, 3);

        answer.at(4, 3 * i - 1) = dnx.at(i, 3);
        answer.at(4, 3 * i - 0) = dnx.at(i, 2);

        answer.at(5, 3 * i - 2) = dnx.at(i, 3);
        answer.at(5, 3 * i - 0) = dnx.at(i, 1);

        answer.at(6, 3 * i - 2) = dnx.at(i, 2);
        answer.at(6, 3 * i - 1) = dnx.at(i, 1);
    }
}

void
tet21ghostsolid :: computeDeformationGradientVectorAt(FloatArray &answer, FloatArray lcoord, TimeStep *tStep)
{
    FloatArray F, u;
    FloatMatrix dNdx, BH, Fmatrix, Finv;
    FEInterpolation *interpolation = this->giveInterpolation();

    // Fetch displacements
    this->computeVectorOf({ 1, 2, 3 }, VM_Total, tStep, u);

    // Compute dNdx in point
    interpolation->evaldNdx( dNdx, lcoord, FEIElementGeometryWrapper(this) );

    // Compute displacement gradient BH
    BH.resize(9, dNdx.giveNumberOfRows() * 3);
    BH.zero();

    for ( int i = 1; i <= dNdx.giveNumberOfRows(); i++ ) {
        BH.at(1, 3 * i - 2) = dNdx.at(i, 1);     // du/dx
        BH.at(2, 3 * i - 1) = dNdx.at(i, 2);     // dv/dy
        BH.at(3, 3 * i - 0) = dNdx.at(i, 3);     // dw/dz
        BH.at(4, 3 * i - 1) = dNdx.at(i, 3);     // dv/dz
        BH.at(7, 3 * i - 0) = dNdx.at(i, 2);     // dw/dy
        BH.at(5, 3 * i - 2) = dNdx.at(i, 3);     // du/dz
        BH.at(8, 3 * i - 0) = dNdx.at(i, 1);     // dw/dx
        BH.at(6, 3 * i - 2) = dNdx.at(i, 2);     // du/dy
        BH.at(9, 3 * i - 1) = dNdx.at(i, 1);     // dv/dx
    }

    // Finally, compute deformation gradient F=BH*u+I
    F.beProductOf(BH, u);
    F.at(1) += 1.0;
    F.at(2) += 1.0;
    F.at(3) += 1.0;

    answer = F;
}

void
tet21ghostsolid :: computeBHmatrixAt(GaussPoint *gp, FloatMatrix &answer)
{
    FloatMatrix dnx;

    this->interpolation.evaldNdx( dnx, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );

    answer.resize(9, 30);
    answer.zero();

    for ( int i = 1; i <= dnx.giveNumberOfRows(); i++ ) {
        answer.at(1, 3 * i - 2) = dnx.at(i, 1);     // du/dx
        answer.at(2, 3 * i - 1) = dnx.at(i, 2);     // dv/dy
        answer.at(3, 3 * i - 0) = dnx.at(i, 3);     // dw/dz
        answer.at(4, 3 * i - 1) = dnx.at(i, 3);     // dv/dz
        answer.at(7, 3 * i - 0) = dnx.at(i, 2);     // dw/dy
        answer.at(5, 3 * i - 2) = dnx.at(i, 3);     // du/dz
        answer.at(8, 3 * i - 0) = dnx.at(i, 1);     // dw/dx
        answer.at(6, 3 * i - 2) = dnx.at(i, 2);     // du/dy
        answer.at(9, 3 * i - 1) = dnx.at(i, 1);     // dv/dx
    }
}

void
tet21ghostsolid :: giveUnknownData(FloatArray &u_prev, FloatArray &u, FloatArray &inc, TimeStep *tStep)
{
    this->computeVectorOf(VM_Total, tStep, u);

    if ( !tStep->isTheFirstStep() ) {
        this->computeVectorOf(VM_Total, tStep->givePreviousStep(), u_prev);
        this->computeVectorOf(VM_Incremental, tStep, inc);
    } else {
        inc.resize( u.giveSize() );
        inc.zero();
        u_prev.operator = ( inc );
    }
}

void
tet21ghostsolid :: computeDeformationGradientVectorFromDispl(FloatArray &answer, GaussPoint *gp, TimeStep *tStep, FloatArray &u)
{
    // Computes the deformation gradient in the Voigt format at the Gauss point gp of
    // the receiver at time step tStep.
    // Order of components: 11, 22, 33, 23, 13, 12, 32, 31, 21 in the 3D.

    // Obtain the current displacement vector of the element and subtract initial displacements (if present)
    if ( initialDisplacements ) {
        u.subtract(* initialDisplacements);
    }

    // Displacement gradient H = du/dX
    FloatMatrix BH;
    this->computeBHmatrixAt(gp, BH);
    answer.beProductOf(BH, u);

    // Deformation gradient F = H + I
    MaterialMode matMode = gp->giveMaterialMode();
    if ( matMode == _3dMat || matMode == _PlaneStrain ) {
        answer.at(1) += 1.0;
        answer.at(2) += 1.0;
        answer.at(3) += 1.0;
    } else if ( matMode == _PlaneStress ) {
        answer.at(1) += 1.0;
        answer.at(2) += 1.0;
    } else if ( matMode == _1dMat ) {
        answer.at(1) += 1.0;
    } else {
        OOFEM_ERROR( "MaterialMode is not supported yet (%s)", __MaterialModeToString(matMode) );
    }
}

void
tet21ghostsolid :: computeDeformationGradientVector(FloatArray &answer, GaussPoint *gp, TimeStep *tStep)
{
    FloatArray u;
    this->computeVectorOf({ 1, 2, 3 }, VM_Total, tStep, u);
    computeDeformationGradientVectorFromDispl(answer, gp, tStep, u);
}

double
tet21ghostsolid :: computeVolumeAround(GaussPoint *gp)
// Returns the portion of the receiver which is attached to gp.
{
    double determinant = fabs( this->interpolation.giveTransformationJacobian( gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) ) );
    double weight      = gp->giveWeight();
    return ( this->computeVolume() );

    return ( determinant * weight );
}


int
tet21ghostsolid :: giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
{
    if ( type == IST_Velocity ) {
        FloatArray N, a;
        FloatMatrix Nmat;

        this->computeVectorOf({ V_u, V_v, V_w }, VM_Total, tStep, a);
        this->interpolation.evalN( N, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );

        Nmat.resize(3, N.giveSize() * 3);
        for ( int i = 1; i <= N.giveSize(); i++ ) {
            Nmat.at(1, 3 * i - 2) = N.at(i);
            Nmat.at(2, 3 * i - 1) = N.at(i);
            Nmat.at(3, 3 * i - 0) = N.at(i);
        }

        answer.beProductOf(Nmat, a);
        return 1;
    } else if ( type == IST_Pressure ) {
        FloatArray N, a;

        this->computeVectorOf({ P_f }, VM_Total, tStep, a);
        this->interpolation_lin.evalN( N, gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );

        answer.resize(1);
        answer.at(1) = N.dotProduct(a);
        return 1;
    } else {
        MaterialMode matmode = gp->giveMaterialMode();
        gp->setMaterialMode(_3dFlow);
        int r = StructuralElement :: giveIPValue(answer, gp, type, tStep);
        gp->setMaterialMode(matmode);
        return r;
    }
}

bool
tet21ghostsolid :: giveRowTransformationMatrix(TimeStep *tStep)
{
    // Create a transformation matrix that switch all rows/equations located in OmegaF but not on GammaInt, i.e where we do not have a no slip condition

    Itransform.resize(64, 64);
    Itransform.zero();

#if USEUNCOUPLED == 0
    Itransform.beUnitMatrix();
    return 0;

#endif

    int m = 1;

    int row1, row2;
    IntArray row1List, row2List;

    for ( int i = 1; i <= this->giveNumberOfDofManagers(); i++ ) {
        IntArray DofIDs;
        int firstIndex = m;
        this->giveDofManDofIDMask(i, DofIDs);
        // DofIDs.printYourself();

        for ( int j = 1; j <= DofIDs.giveSize(); j++ ) {
            row1 = m;
            row2 = m;

            if ( DofIDs.at(j) == V_u || DofIDs.at(j) == V_v || DofIDs.at(j) == V_w ) {
                bool doSwitch = !this->giveDofManager(i)->giveDofWithID( DofIDs.at(j) )->hasBc(tStep);

                if ( doSwitch ) {     // Boundary condition not set, make switch
                    row1 = m;

                    for ( int n = 1; n <= DofIDs.giveSize(); n++ ) {
                        if ( ( DofIDs.at(j) == V_u && DofIDs.at(n) == D_u ) ||
                             ( DofIDs.at(j) == V_v && DofIDs.at(n) == D_v ) ||
                             ( DofIDs.at(j) == V_w && DofIDs.at(n) == D_w ) ) {
                            row2 = firstIndex + n - 1;
                        }
                    }
                    row1List.resizeWithValues(row1List.giveSize() + 1);
                    row1List.at( row1List.giveSize() ) = row1;
                    row2List.resizeWithValues(row2List.giveSize() + 1);
                    row2List.at( row2List.giveSize() ) = row2;
                }
            }
            m++;
        }
    }

    // Create identity matrix
    Itransform.beUnitMatrix();

    // Create tranformation matrix by switching rows
    for ( int i = 1; i <= row1List.giveSize(); i++ ) {
        Itransform.at( row1List.at(i), row1List.at(i) ) = 0;
        Itransform.at( row2List.at(i), row2List.at(i) ) = 0;

        Itransform.at( row1List.at(i), row2List.at(i) ) = 1;
        Itransform.at( row2List.at(i), row1List.at(i) ) = 1;
    }

    computeItransform = false;

    return 0;
}

// Some extension Interfaces to follow:

Interface *tet21ghostsolid :: giveInterface(InterfaceType it)
{
    switch ( it ) {
    case NodalAveragingRecoveryModelInterfaceType:
        return static_cast< NodalAveragingRecoveryModelInterface * >(this);

    case SpatialLocalizerInterfaceType:
        return static_cast< SpatialLocalizerInterface * >(this);

    case EIPrimaryUnknownMapperInterfaceType:
        return static_cast< EIPrimaryUnknownMapperInterface * >(this);

    default:
        return StructuralElement :: giveInterface(it);
        //return FMElement :: giveInterface(it);
    }
}

void tet21ghostsolid :: EIPrimaryUnknownMI_computePrimaryUnknownVectorAtLocal(ValueModeType mode,
                                                                              TimeStep *tStep, const FloatArray &lcoords, FloatArray &answer)
{
    FloatArray n, n_lin;
    this->interpolation.evalN( n, lcoords, FEIElementGeometryWrapper(this) );
    this->interpolation_lin.evalN( n_lin, lcoords, FEIElementGeometryWrapper(this) );
    answer.resize(4);
    answer.zero();
    for ( int i = 1; i <= n.giveSize(); i++ ) {
        answer(0) += n.at(i) * this->giveNode(i)->giveDofWithID(V_u)->giveUnknown(mode, tStep);
        answer(1) += n.at(i) * this->giveNode(i)->giveDofWithID(V_v)->giveUnknown(mode, tStep);
        answer(2) += n.at(i) * this->giveNode(i)->giveDofWithID(V_w)->giveUnknown(mode, tStep);
    }

    for ( int i = 1; i <= n_lin.giveSize(); i++ ) {
        answer(3) += n_lin.at(i) * this->giveNode(i)->giveDofWithID(P_f)->giveUnknown(mode, tStep);
    }
}

void
tet21ghostsolid :: NodalAveragingRecoveryMI_computeNodalValue(FloatArray &answer, int node, InternalStateType type, TimeStep *tStep)
{
    if ( type == IST_Pressure ) {
        answer.resize(1);
        if ( node <= 4 ) {
            answer.at(1) = this->giveNode(node)->giveDofWithID(P_f)->giveUnknown(VM_Total, tStep);
        } else {
            IntArray eNodes;
            this->interpolation.computeLocalEdgeMapping(eNodes, node - 4);
            answer.at(1) = 0.5 * (
                this->giveNode( eNodes.at(1) )->giveDofWithID(P_f)->giveUnknown(VM_Total, tStep) +
                this->giveNode( eNodes.at(2) )->giveDofWithID(P_f)->giveUnknown(VM_Total, tStep) );
        }
    } else {
        answer.clear();
    }
}

void
tet21ghostsolid :: computeBoundarySurfaceLoadVector(FloatArray &answer, BoundaryLoad *load, int boundary, CharType type, ValueModeType mode, TimeStep *tStep, bool global)
{
    answer.clear();
    if ( type != ExternalForcesVector ) {
        return;
    }

    FEInterpolation *fei = this->giveInterpolation();
    if ( !fei ) {
        OOFEM_ERROR("No interpolator available");
    }

    FloatArray n_vec, f(18);
    FloatMatrix n, T;
    FloatArray force;
    int nsd = fei->giveNsd();

    f.zero();
    IntegrationRule *iRule = fei->giveBoundaryIntegrationRule(load->giveApproxOrder(), boundary);

    for ( GaussPoint *gp : *iRule ) {
        FloatArray lcoords = gp->giveNaturalCoordinates();
        if ( load->giveFormulationType() == Load :: FT_Entity ) {
            load->computeValueAt(force, tStep, lcoords, mode);
        } else {
            FloatArray gcoords, elcoords;
            this->interpolation.surfaceLocal2global( gcoords, boundary, lcoords, FEIElementGeometryWrapper(this) );
            this->interpolation.global2local( elcoords, gcoords, FEIElementGeometryWrapper(this) );
            NeumannMomentLoad *thisLoad = dynamic_cast< NeumannMomentLoad * >(load);
            if ( thisLoad != NULL ) {
                FloatArray temp;
                thisLoad->computeValueAtBoundary(temp, tStep, gcoords, VM_Total, this, boundary);

                FloatArray F;
                FloatMatrix Fm, Finv, FinvT;
                this->computeDeformationGradientVectorAt(F, elcoords, tStep);
                Fm.beMatrixForm(F);
                double J = Fm.giveDeterminant();
                Finv.beInverseOf(Fm);
                FinvT.beTranspositionOf(Finv);
                force.beTProductOf(Finv, temp);
                force.times(J);
            } else {
                load->computeValueAt(force, tStep, gcoords, VM_Total);
            }
        }

        // We always want the global values in the end, so we might as well compute them here directly:
        // transform force
        if ( load->giveCoordSystMode() == Load :: CST_Global ) {
            // then just keep it in global c.s
        } else {
            // transform from local boundary to element local c.s
            /*if ( this->computeLoadLSToLRotationMatrix(T, boundary, gp) ) {
             *  force.rotatedWith(T, 'n');
             * }*/
            // then to global c.s
            if ( this->computeLoadGToLRotationMtrx(T) ) {
                force.rotatedWith(T, 't');
            }
        }

        // Construct n-matrix
        fei->boundaryEvalN( n_vec, boundary, lcoords, FEIElementGeometryWrapper(this) );
        n.beNMatrixOf(n_vec, nsd);

        double thickness = 1.0; // Should be the circumference for axisymm-elements.
        double dV = thickness * gp->giveWeight() * fei->boundaryGiveTransformationJacobian( boundary, lcoords, FEIElementGeometryWrapper(this) );
        f.plusProduct(n, force, dV);
    }

    FloatArray temp;
    FloatMatrix Itrans;

    answer.resize(39);
    answer.zero();
    answer.assemble(f, this->velocitydofsonside);

    delete iRule;
}
}