fe2fluidmaterial.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 "fm/Materials/fe2fluidmaterial.h"
#include "stokesflow.h"
#include "oofemtxtdatareader.h"
#include "domain.h"
#include "gausspoint.h"
#include "contextioerr.h"
#include "util.h"
#include "classfactory.h"
#include "dynamicinputrecord.h"
#include "mathfem.h"
 
#include <sstream>
 
//#define DEBUG_TANGENT
#define DEBUG_ERR ( 1e-6 )
 
namespace oofem {
REGISTER_Material(FE2FluidMaterial);
 
int FE2FluidMaterial :: n = 1;
 
void FE2FluidMaterialStatus :: setTimeStep(TimeStep *tStep)
{
    TimeStep *rveTStep = this->rve->giveCurrentStep(); // Should i create a new one if it is empty?
    rveTStep->setNumber( tStep->giveNumber() );
    rveTStep->setTime( tStep->giveTargetTime() );
    rveTStep->setTimeIncrement( tStep->giveTimeIncrement() );
}
 
FloatArrayF<6> FE2FluidMaterial :: computeDeviatoricStress3D(const FloatArrayF<6> &eps, GaussPoint *gp, TimeStep *tStep) const
{
    auto val = this->computeDeviatoricStress3D(eps, 0.0, gp, tStep);
    double r_vol = val.second + eps.at(1) + eps.at(2) + eps.at(3);
    if ( r_vol > 1e-9 ) {
        OOFEM_ERROR("RVE seems to be compressible;"
                    " extended macro-formulation which doesn't assume incompressibility is required");
    }
    return val.first;
}
 
std::pair<FloatArrayF<6>, double> FE2FluidMaterial :: computeDeviatoricStress3D(const FloatArrayF<6> &eps, double pressure, GaussPoint *gp, TimeStep *tStep) const
{
    FE2FluidMaterialStatus *ms = static_cast< FE2FluidMaterialStatus * >( this->giveStatus(gp) );
 
    ms->setTimeStep(tStep);
 
    MixedGradientPressureBC *bc = ms->giveBC();
 
    // Set input
    bc->setPrescribedDeviatoricGradientFromVoigt(eps);
    bc->setPrescribedPressure(pressure);
    // Solve subscale problem
    ms->giveRVE()->solveYourselfAt(ms->giveRVE()->giveCurrentStep());
 
    FloatArray stress_dev;
    double r_vol;
    bc->computeFields(stress_dev, r_vol, tStep);
 
    ms->letDeviatoricStressVectorBe(stress_dev);
    ms->letDeviatoricStrainRateVectorBe(eps);
    ms->letPressureBe(pressure);
    ms->markOldTangents(); // Mark this so that tangents are reevaluated if they are needed.
    // One could also just compute them here, but you don't actually need them if the problem has converged, so this method saves on that iteration.
    // Computing the tangents are often *more* expensive than computeFields, so this is well worth the time it saves
    return {stress_dev, r_vol};
}
 
FloatMatrixF<6,6> FE2FluidMaterial :: computeTangent3D(MatResponseMode mode, GaussPoint *gp, TimeStep *tStep) const
{
    FE2FluidMaterialStatus *ms = static_cast< FE2FluidMaterialStatus * >( this->giveStatus(gp) );
    ms->computeTangents(tStep);
    if ( mode != TangentStiffness ) {
        OOFEM_ERROR("Mode not implemented");
    }
    FloatMatrixF<6,6> x = ms->giveDeviatoricTangent();
    return x; //ms->giveDeviatoricTangent();
}
 
FluidDynamicMaterial::Tangents<6>
FE2FluidMaterial :: computeTangents3D(MatResponseMode mode, GaussPoint *gp, TimeStep *tStep) const
{
    //FloatMatrix &dsdd, FloatArray &dsdp, FloatArray &dedd, double &dedp,
    FE2FluidMaterialStatus *ms = static_cast< FE2FluidMaterialStatus * >( this->giveStatus(gp) );
    ms->computeTangents(tStep);
    if ( mode != TangentStiffness ) {
        OOFEM_ERROR("Mode not implemented");
    }
#if 0
    // Numerical ATS for debugging
    FloatMatrix numericalATS(6, 6);
    FloatArray dsig;
    FloatArray tempStrain(6);
 
    tempStrain.zero();
    FloatArray sig, strain, sigPert;
    double epspvol;
    computeDeviatoricStressVector(sig, epspvol, gp, tempStrain, 0., tStep);
    double h = 0.001; // Linear problem, size of this doesn't matter.
    for ( int k = 1; k <= 6; ++k ) {
        strain = tempStrain;
        strain.at(k) += h;
        double tmp = strain.at(1) + strain.at(2) + strain.at(3);
        strain.at(1) -= tmp/3.0;
        strain.at(2) -= tmp/3.0;
        strain.at(3) -= tmp/3.0;
        strain.printYourself();
        computeDeviatoricStressVector(sigPert, epspvol, gp, strain, 0., tStep);
        sigPert.printYourself();
        dsig.beDifferenceOf(sigPert, sig);
        numericalATS.setColumn(dsig, k);
    }
    numericalATS.times(1. / h);
 
    printf("Analytical deviatoric tangent = ");
    dsdd.printYourself();
    printf("Numerical deviatoric tangent = ");
    numericalATS.printYourself();
    numericalATS.subtract(dsdd);
    double norm = numericalATS.computeFrobeniusNorm();
    if ( norm > dsdd.computeFrobeniusNorm() * DEBUG_ERR && norm > 0.0 ) {
        OOFEM_ERROR("Error in deviatoric tangent");
    }
#endif
#if 0
    // Numerical ATS for debugging
    FloatArray strain(3);
    strain.zero();
    FloatArray sig, sigh;
    double epspvol, pressure = 0.0;
    double h = 1.00; // Linear problem, size of this doesn't matter.
    computeDeviatoricStressVector(sig, epspvol, gp, strain, pressure, tStep);
    computeDeviatoricStressVector(sigh, epspvol, gp, strain, pressure + h, tStep);
 
    FloatArray dsigh;
    dsigh.beDifferenceOf(sigh, sig);
    dsigh.times(1 / h);
 
    printf("Analytical deviatoric pressure tangent = ");
    dsdp.printYourself();
    printf("Numerical deviatoric pressure tangent = ");
    dsigh.printYourself();
    dsigh.subtract(dsdp);
    double norm = dsigh.computeNorm();
    if ( norm > dsdp.computeNorm() * DEBUG_ERR && norm > 0.0 ) {
        OOFEM_ERROR("Error in deviatoric pressure tangent");
    }
#endif
#if 0
    // Numerical ATS for debugging
    FloatArray tempStrain(3);
    tempStrain.zero();
    FloatArray sig, strain;
    double epspvol, epspvol11, epspvol22, epspvol12, pressure = 0.0;
    double h = 1.0; // Linear problem, size of this doesn't matter.
 
    computeDeviatoricStressVector(sig, epspvol, gp, tempStrain, pressure, tStep);
    strain = tempStrain;
    strain.at(1) += h;
    computeDeviatoricStressVector(sig, epspvol11, gp, strain, pressure, tStep);
    strain = tempStrain;
    strain.at(2) += h;
    computeDeviatoricStressVector(sig, epspvol22, gp, strain, pressure, tStep);
    strain = tempStrain;
    strain.at(3) += h;
    computeDeviatoricStressVector(sig, epspvol12, gp, strain, pressure, tStep);
 
    FloatArray dvol(3);
    dvol.at(1) = ( epspvol11 - epspvol ) / h;
    dvol.at(2) = ( epspvol22 - epspvol ) / h;
    dvol.at(3) = ( epspvol12 - epspvol ) / h;
    dvol.at(1) += 1.0;
    dvol.at(2) += 1.0;
 
    printf("Analytical volumetric deviatoric tangent = ");
    dedd.printYourself();
    printf("Numerical volumetric deviatoric tangent = ");
    dvol.printYourself();
    dvol.subtract(dedd);
    double norm = dvol.computeNorm();
    if ( norm > dedd.computeNorm() * DEBUG_ERR && norm > 0.0 ) {
        OOFEM_ERROR("Error in volumetric deviatoric tangent");
    }
#endif
#if 0
    // Numerical ATS for debugging
    FloatArray strain(3);
    strain.zero();
    FloatArray sig;
    double epspvol, epspvolh, pressure = 0.0;
    double h = 1.0; // Linear problem, size of this doesn't matter.
 
    computeDeviatoricStressVector(sig, epspvol, gp, strain, pressure, tStep);
    computeDeviatoricStressVector(sig, epspvolh, gp, strain, pressure + h, tStep);
 
    double dvol = -( epspvolh - epspvol ) / h;
 
    printf("Analytical volumetric pressure tangent = %e\n", dedp);
    printf("Numerical volumetric pressure tangent = %e\n", dvol);
 
    double norm = fabs(dvol - dedp);
    if ( norm > fabs(dedp) * DEBUG_ERR && norm > 0.0 ) {
        OOFEM_ERROR("Error in volumetric pressure tangent");
    }
#endif
    return {
        ms->giveDeviatoricTangent(),
        ms->giveDeviatoricPressureTangent(),
        ms->giveVolumetricDeviatoricTangent(),
        ms->giveVolumetricPressureTangent(),
    };
}
 
void FE2FluidMaterial :: initializeFrom(InputRecord &ir)
{
    FluidDynamicMaterial :: initializeFrom(ir);
    IR_GIVE_FIELD(ir, this->inputfile, _IFT_FE2FluidMaterial_fileName);
}
 
void FE2FluidMaterial :: giveInputRecord(DynamicInputRecord &input)
{
    FluidDynamicMaterial :: giveInputRecord(input);
    input.setField(this->inputfile, _IFT_FE2FluidMaterial_fileName);
}
 
 
int FE2FluidMaterial :: giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
{
    FE2FluidMaterialStatus *status = static_cast< FE2FluidMaterialStatus * >( this->giveStatus(gp) );
    if ( type == IST_VOFFraction ) {
        answer = FloatArray{status->giveVOFFraction()};
        return true;
    } else if ( type == IST_Pressure ) {
        answer = FloatArray{status->givePressure()};
        return true;
    } else if ( type == IST_Undefined ) { 
#if 0
        // Numerical ATS for debugging
        FloatArray strain(3);
        strain.zero();
        FloatArray sig;
        double epspvol, epspvolh, pressure = 0.0;
        double h = 1.0; // Linear problem, size of this doesn't matter.
 
        computeDeviatoricStressVector(sig, epspvol, gp, strain, pressure, tStep);
        computeDeviatoricStressVector(sig, epspvolh, gp, strain, pressure + h, tStep);
 
        double dvol = - ( epspvolh - epspvol ) / h;
 
        
        printf("Analytical volumetric pressure tangent = %f\n", status->giveVolumetricPressureTangent());
        printf("Numerical volumetric pressure tangent = %f\n", dvol);
 
        double norm = fabs(dvol - status->giveVolumetricPressureTangent());
        if ( norm > fabs(status->giveVolumetricPressureTangent()) * DEBUG_ERR && norm > 0.0 ) {
            OOFEM_ERROR("Error in volumetric pressure tangent");
        }
#endif
        answer = FloatArray{status->giveVolumetricPressureTangent()};
        return true;
    } else {
        return FluidDynamicMaterial :: giveIPValue(answer, gp, type, tStep);
    }
}
 
 
std::unique_ptr<MaterialStatus> FE2FluidMaterial :: CreateStatus(GaussPoint *gp) const
{
    int rank = -1;
    if ( this->domain->giveEngngModel()->isParallel() && this->domain->giveEngngModel()->giveNumberOfProcesses() > 1 ) {
        rank = this->domain->giveEngngModel()->giveRank();
    }
    return std::make_unique<FE2FluidMaterialStatus>(n++, rank, gp, this->inputfile);
}
 
int FE2FluidMaterial :: checkConsistency()
{
    return true;
}
 
FE2FluidMaterialStatus :: FE2FluidMaterialStatus(int n, int rank, GaussPoint *gp, const std :: string &inputfile) :
    FluidDynamicMaterialStatus(gp),
    voffraction(0.0),
    oldTangents(true)
{
    if ( !this->createRVE(n, rank ,gp, inputfile) ) {
        OOFEM_ERROR("Couldn't create RVE");
    }
}
 
// Uses an input file for now, should eventually create the RVE itself.
bool FE2FluidMaterialStatus :: createRVE(int n, int rank, GaussPoint *gp, const std :: string &inputfile)
{
    OOFEMTXTDataReader dr( inputfile.c_str() );
    this->rve = InstanciateProblem(dr, _processor, 0); // Everything but nrsolver is updated.
    dr.finish();
    this->rve->setProblemScale(microScale);
    this->rve->checkProblemConsistency();
    this->rve->initMetaStepAttributes( this->rve->giveMetaStep(1) );
    this->rve->giveNextStep(); // Makes sure there is a timestep (which we will modify before solving a step)
    this->rve->init();
 
    std :: ostringstream name;
    name << this->rve->giveOutputBaseFileName() << "-gp" << n;
    if ( rank >= 0 ) {
        name << "." << rank;
    }
 
    this->rve->letOutputBaseFileNameBe( name.str() );
 
    this->bc = dynamic_cast< MixedGradientPressureBC * >( this->rve->giveDomain(1)->giveBc(1) );
    if ( !this->bc ) {
        OOFEM_ERROR("RVE doesn't have necessary boundary condition; should have MixedGradientPressure as first b.c. (in first domain)");
    }
 
    return true;
}
 
void FE2FluidMaterialStatus :: printOutputAt(FILE *file, TimeStep *tStep) const
{
    FluidDynamicMaterialStatus :: printOutputAt(file, tStep);
}
 
void FE2FluidMaterialStatus :: updateYourself(TimeStep *tStep)
{
    double fluid_area = this->rve->giveDomain(1)->giveSize();
    double total_area = this->bc->domainSize();
    this->voffraction = fluid_area / total_area;
    FluidDynamicMaterialStatus :: updateYourself(tStep);
 
    this->rve->updateYourself(tStep);
    this->rve->terminate(tStep);
}
 
void FE2FluidMaterialStatus :: initTempStatus()
{
    FluidDynamicMaterialStatus :: initTempStatus();
}
 
void FE2FluidMaterialStatus :: saveContext(DataStream &stream, ContextMode mode)
{
    FluidDynamicMaterialStatus :: saveContext(stream, mode);
    this->rve->saveContext(stream, mode);
}
 
void FE2FluidMaterialStatus :: restoreContext(DataStream &stream, ContextMode mode)
{
    FluidDynamicMaterialStatus :: restoreContext(stream, mode);
    this->markOldTangents();
    this->rve->restoreContext(stream, mode);
}
 
void FE2FluidMaterialStatus :: markOldTangents() { this->oldTangents = true; }
 
void FE2FluidMaterialStatus :: computeTangents(TimeStep *tStep)
{
    if ( !tStep->isTheCurrentTimeStep() ) {
        OOFEM_ERROR("Only current timestep supported.");
    }
 
    if ( this->oldTangents ) {
        bc->computeTangents(this->giveDeviatoricTangent(),
                            this->giveDeviatoricPressureTangent(),
                            this->giveVolumetricDeviatoricTangent(),
                            this->giveVolumetricPressureTangent(),
                            tStep);
    }
 
    this->oldTangents = false;
}
 
double FE2FluidMaterial :: giveEffectiveViscosity(GaussPoint *gp, TimeStep *tStep) const
{
    FE2FluidMaterialStatus *status = static_cast< FE2FluidMaterialStatus * >( this->giveStatus(gp) );
    status->computeTangents(tStep);
    const auto &t = status->giveDeviatoricTangent();
    // Project against the normalized I_dev
    double v;
    if ( gp->giveMaterialMode() == _3dFlow ) {
        v = ( t.at(1, 1) * 2. - t.at(1, 2)    - t.at(1, 3) ) / 3. +
            ( -t.at(2, 1)    + t.at(2, 2) * 2. - t.at(2, 3) ) / 3. +
            ( -t.at(3, 1)    - t.at(3, 2)    + t.at(3, 3) * 2. ) / 3. +
            4. * ( t.at(4, 4) + t.at(5, 5) + t.at(6, 6) );
        v /= 16.;
    } else {
        v = ( t.at(1, 1) * 2. - t.at(1, 2) ) / 3. +
            ( -t.at(2, 1)    + t.at(2, 2) * 2. ) / 3. +
            4. * t.at(3, 3);
        v *= 9. / 56.;
    }
 
    return v;
}
} // end namespace oofem