matlabexportmodule.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 <algorithm>
#include <vector>
#include <string>
#include <iostream>
#include <sstream>
#include <iterator>
 
#include "matlabexportmodule.h"
#include "engngm.h"
#include "node.h"
#include "mathfem.h"
#include "gausspoint.h"
#include "weakperiodicbc.h"
#include "solutionbasedshapefunction.h"
#include "timestep.h"
#include "classfactory.h"
#include "set.h"
#include "unknownnumberingscheme.h"
#include "prescribedmean.h"
#include "feinterpol.h"
    
#ifdef __FM_MODULE
#include "fm/tr21stokes.h"
#include "fm/tet21stokes.h"
#include "fm/stokesflow.h"
#endif
 
#ifdef __SM_MODULE
#include "sm/Elements/nlstructuralelement.h"
#include "sm/EngineeringModels/structengngmodel.h"
#include "sm/prescribeddispslipbcdirichletrc.h"
#include "sm/prescribeddispslipbcneumannrc.h"
#include "sm/prescribeddispslipmultiple.h"
#include "sm/transversereinfconstraint.h"
#endif
 
 
namespace oofem {
 
REGISTER_ExportModule( MatlabExportModule )
 
MatlabExportModule :: MatlabExportModule(int n, EngngModel *e) : ExportModule(n, e), internalVarsToExport(), primaryVarsToExport()
{
    exportMesh = false;
    exportData = false;
    exportArea = false;
    exportSpecials = false;
    exportReactionForces = false;
    reactionForcesDofManList.clear();
    dataDofManList.clear();
    exportIntegrationPointFields = false;
    elList.clear();
    reactionForcesNodeSet = 0;
    dataNodeSet = 0;
    IPFieldsElSet = 0;
    noscaling = false;
}
 
 
MatlabExportModule :: ~MatlabExportModule()
{ }
 
 
void
MatlabExportModule :: initializeFrom(InputRecord &ir)
{
    ExportModule :: initializeFrom(ir);
 
    exportMesh = ir.hasField(_IFT_MatlabExportModule_mesh);
    
    exportData = ir.hasField(_IFT_MatlabExportModule_data);
    if ( exportData ) {
        IR_GIVE_OPTIONAL_FIELD(ir, this->dataNodeSet, _IFT_MatlabExportModule_DataNodeSet);
    }
    
    exportArea = ir.hasField(_IFT_MatlabExportModule_area);
    exportSpecials = ir.hasField(_IFT_MatlabExportModule_specials);
    exportHomogenizeIST = ir.hasField(_IFT_MatlabExportModule_homogenizeInternalVars);
 
 
    exportReactionForces = ir.hasField(_IFT_MatlabExportModule_ReactionForces);
    reactionForcesDofManList.resize(0);
    if ( exportReactionForces ) {
        IR_GIVE_OPTIONAL_FIELD(ir, reactionForcesDofManList, _IFT_MatlabExportModule_DofManList);
        this->reactionForcesNodeSet = 0;
        IR_GIVE_OPTIONAL_FIELD(ir, this->reactionForcesNodeSet, _IFT_MatlabExportModule_ReactionForcesNodeSet);
    }
 
    exportIntegrationPointFields = ir.hasField(_IFT_MatlabExportModule_IntegrationPoints);
    elList.resize(0);
    IR_GIVE_OPTIONAL_FIELD(ir, internalVarsToExport, _IFT_MatlabExportModule_internalVarsToExport);
    if ( exportIntegrationPointFields ) {
        IR_GIVE_OPTIONAL_FIELD(ir, elList, _IFT_MatlabExportModule_ElementList);
        IPFieldsElSet = 0;
        IR_GIVE_OPTIONAL_FIELD(ir, IPFieldsElSet, _IFT_MatlabExportModule_IPFieldsElSet);
    }
 
    noscaling = ir.hasField(_IFT_MatlabExportModule_noScaledHomogenization);
 
}
 
 
 
void
MatlabExportModule :: computeArea(TimeStep *tStep)
{
    Domain *domain = emodel->giveDomain(1);
 
    smax.clear();
    smin.clear();
 
    for (int i = 1; i <= domain->giveNumberOfSpatialDimensions(); i++) {
        smax.push_back(domain->giveDofManager(1)->giveCoordinate(i));
        smin.push_back(domain->giveDofManager(1)->giveCoordinate(i));
    }
 
    for ( int i = 0; i < domain->giveNumberOfDofManagers(); i++ ) {
        for (int j = 0; j < domain->giveNumberOfSpatialDimensions(); j++) {
            smax.at(j)=max(smax.at(j), domain->giveDofManager(i+1)->giveCoordinate(j+1));
            smin.at(j)=min(smin.at(j), domain->giveDofManager(i+1)->giveCoordinate(j+1));
        }
    }
 
    Area = 0;
    Volume = 0;
 
    if ( domain->giveNumberOfSpatialDimensions() == 2 ) {
        for ( auto &elem : domain->giveElements() ) {
            Area += elem->computeArea();
        }
    } else {
 
        for (size_t i = 0; i < partVolume.size(); i++ ) {
            partVolume.at(i) = 0.0;
        }
 
        for ( auto &elem : domain->giveElements() ) {
            //
            double v;
#ifdef __SM_MODULE
            if ( NLStructuralElement *e = dynamic_cast< NLStructuralElement *>( elem.get() ) ) {
                v = e->computeCurrentVolume(tStep);
            } else {
#endif
                v = elem->computeVolume();
#ifdef __SM_MODULE
            }
#endif
 
            std :: string eName ( elem->giveClassName() );
            int j = -1;
 
            //printf("%s\n", eName.c_str());
 
            for ( size_t k = 0; k < partName.size(); k++ ) {
                //printf("partName.at(%u) = %s\n", k, partName.at(k).c_str() );
                if ( eName.compare(partName.at(k)) == 0 ) {
                    j = k;
                    break;
                }
            }
 
            if ( j == -1 ) {
                partName.push_back( elem->giveClassName() );
                partVolume.push_back( v );
            } else {
                partVolume.at(j) += v;
            }
 
            Volume += v;
 
        }
    }
 
}
 
 
void
MatlabExportModule :: doOutput(TimeStep *tStep, bool forcedOutput)
{
    if ( !( testTimeStepOutput(tStep) || forcedOutput ) ) {
        return;
    }
 
 
    int nelem = this->elList.giveSize();
    if ( nelem == 0 ) { // no list given, export all elements
        this->elList.enumerate(this->emodel->giveDomain(1)->giveNumberOfElements());
    }
 
    FILE *FID;
    FID = giveOutputStream(tStep);
    Domain *domain  = emodel->giveDomain(1);
    ndim=domain->giveNumberOfSpatialDimensions();
 
    // Output header
    fprintf( FID, "%%%% OOFEM generated export file \n");
    fprintf( FID, "%% Output for time %f\n", tStep->giveTargetTime() );
 
 
    fprintf( FID, "function [mesh area data specials ReactionForces IntegrationPointFields]=%s\n\n", functionname.c_str() );
 
    if ( exportMesh ) {
        doOutputMesh(tStep, FID);
    } else {
        fprintf(FID, "\tmesh=[];\n");
    }
 
    if ( exportData ) {
        doOutputData(tStep, FID);
    } else {
        fprintf(FID, "\tdata=[];\n");
    }
 
    if ( exportArea ) {
        computeArea(tStep);
        fprintf(FID, "\tarea.xmax=%f;\n", smax.at(0));
        fprintf(FID, "\tarea.xmin=%f;\n", smin.at(0));
        fprintf(FID, "\tarea.ymax=%f;\n", smax.at(1));
        fprintf(FID, "\tarea.ymin=%f;\n", smin.at(1));
        if ( ndim == 2 ) {
            fprintf(FID, "\tarea.area=%f;\n", Area);
            fprintf(FID, "\tvolume=[];\n");
        } else {
            fprintf(FID, "\tarea.zmax=%f;\n", smax.at(2));
            fprintf(FID, "\tarea.zmin=%f;\n", smin.at(2));
            fprintf(FID, "\tarea.area=[];\n");
            fprintf(FID, "\tarea.volume=%f;\n", Volume);
            for (size_t i=0; i<this->partName.size(); i++) {
                fprintf(FID, "\tarea.volume_%s=%f;\n", partName.at(i).c_str(), partVolume.at(i));
            }
        }
    } else {
        fprintf(FID, "\tarea.area=[];\n");
        fprintf(FID, "\tarea.volume=[];\n");
    }
 
    if ( exportSpecials ) {
        if ( !exportArea ) {
            computeArea(tStep);
        }
 
        doOutputSpecials(tStep, FID);
    } else {
        fprintf(FID, "\tspecials=[];\n");
    }
 
    // Reaction forces
    if ( exportReactionForces ) {
        doOutputReactionForces(tStep, FID);
    } else {
        fprintf(FID, "\tReactionForces=[];\n");
    }
 
    // Internal variables in integration points
    if ( exportIntegrationPointFields ) {
        doOutputIntegrationPointFields(tStep, FID);
    } else {
        fprintf(FID, "\tIntegrationPointFields=[];\n");
    }
 
    // Homogenized quantities
    if ( exportHomogenizeIST ) {
        doOutputHomogenizeDofIDs(tStep, FID);
    }
 
    fprintf(FID, "\nend\n");
    fclose(FID);
}
 
 
void
MatlabExportModule :: doOutputMesh(TimeStep *tStep, FILE *FID)
{
    Domain *domain  = emodel->giveDomain(1);
 
    fprintf(FID, "\tmesh.p=[");
    for ( auto &dman : domain->giveDofManagers() ) {
        for ( int j = 1; j <= domain->giveNumberOfSpatialDimensions(); j++) {
            double c = dman->giveCoordinate(j);
            fprintf(FID, "%f, ", c);
        }
        fprintf(FID, "; ");
    }
 
    fprintf(FID, "]';\n");
 
    int numberOfDofMans=domain->giveElement(1)->giveNumberOfDofManagers();
 
    fprintf(FID, "\tmesh.t=[");
    for ( auto &elem : domain->giveElements() ) {
        if ( elem->giveNumberOfDofManagers() == numberOfDofMans ) {
            for ( int j = 1; j <= elem->giveNumberOfDofManagers(); j++ ) {
                fprintf( FID, "%d,", elem->giveDofManagerNumber(j) );
            }
        }
        fprintf(FID, ";");
    }
 
    fprintf(FID, "]';\n");
}
 
 
void
MatlabExportModule :: doOutputData(TimeStep *tStep, FILE *FID)
{
    Domain *domain  = emodel->giveDomain(1);
    std :: vector< int >DofIDList;
    std :: vector< int > :: iterator it;
    std :: vector< std :: vector< double > >valuesList;
    
    if ( this->dataNodeSet > 0 ) {
        // Export data based on node set
        Set *set = domain->giveSet( this->dataNodeSet );
        dataDofManList = set->giveNodeList();
        for (auto iDM : dataDofManList) {
            DofManager *dman = domain->giveDofManager(iDM);
            for ( Dof *thisDof: *dman ) {
                it = std :: find( DofIDList.begin(), DofIDList.end(), thisDof->giveDofID() );
 
                double value = thisDof->giveUnknown(VM_Total, tStep);
                if ( it == DofIDList.end() ) {
                    DofIDList.push_back( thisDof->giveDofID() );
                    valuesList.push_back({value});
                } else {
                    std::size_t pos = it - DofIDList.begin();
                    valuesList[pos].push_back(value);
                }
            }
        }
    } else {
        // Export data from all dofmanagers
        for ( auto &dman : domain->giveDofManagers() ) {
            for ( Dof *thisDof: *dman ) {
                it = std :: find( DofIDList.begin(), DofIDList.end(), thisDof->giveDofID() );
 
                double value = thisDof->giveUnknown(VM_Total, tStep);
                if ( it == DofIDList.end() ) {
                    DofIDList.push_back( thisDof->giveDofID() );
                    valuesList.push_back({value});
                } else {
                    std::size_t pos = it - DofIDList.begin();
                    valuesList[pos].push_back(value);
                }
            }
        }
    }
 
 
 
    fprintf(FID, "\tdata.DofIDs=[");
    for ( auto &dofid : DofIDList ) {
        fprintf( FID, "%d, ", dofid );
    }
 
    fprintf(FID, "];\n");
 
    for ( size_t i = 0; i < valuesList.size(); i++ ) {
        fprintf(FID, "\tdata.a{%lu}=[", static_cast< long unsigned int >(i) + 1);
        for ( double val: valuesList[i] ) {
            fprintf( FID, "%f,", val );
        }
 
        fprintf(FID, "];\n");
    }
 
}
 
 
void
MatlabExportModule :: doOutputSpecials(TimeStep *tStep,    FILE *FID)
{
//    FloatArray v_hat, GradPTemp, v_hatTemp;
 
    Domain *domain  = emodel->giveDomain(1);
/*
    v_hat.clear();

#if 0
    for ( auto &elem : domain->giveElements() ) {
#ifdef __FM_MODULE
 
        if ( Tr21Stokes *Tr = dynamic_cast< Tr21Stokes * >( elem.get() ) ) {
            Tr->giveIntegratedVelocity(v_hatTemp, tStep);
            v_hat.add(v_hatTemp);
        } else if ( Tet21Stokes *Tet = dynamic_cast< Tet21Stokes * >( elem.get() ) ) {
            Tet->giveIntegratedVelocity(v_hatTemp, tStep);
            v_hat.add(v_hatTemp);
        }
 
#endif
    }
#endif
 
    // Compute intrinsic area/volume
    double intrinsicSize = 1.0;
 
    std :: vector <double> V;
 
    for ( int i = 0; i < (int)smax.size(); i++ ) {
        intrinsicSize *= ( smax.at(i) - smin.at(i) );
    }
 
    for ( double vh: v_hat ) {
        V.push_back(vh);
    }
 
    fprintf(FID, "\tspecials.velocitymean=[");
    if (V.size()>0) {
        for (int i=0; i<ndim; i++) {
            fprintf(FID, "%e", V.at(i));
            if (i!=(ndim-1)) fprintf (FID, ", ");
        }
        fprintf(FID, "];\n");
    } else {
        fprintf(FID, "]; %% No velocities\n");
    }
 
    */
 
    // Output weak periodic boundary conditions
    unsigned int wpbccount = 1, sbsfcount = 1, mcount = 1, pdsdcount=1, pdsncount=1, trccount=1;
    [[maybe_unused]] unsigned int pdsmcount=1;
 
    for ( auto &gbc : domain->giveBcs() ) {
        WeakPeriodicBoundaryCondition *wpbc = dynamic_cast< WeakPeriodicBoundaryCondition * >( gbc.get() );
        if ( wpbc ) {
            for ( int j = 1; j <= wpbc->giveNumberOfInternalDofManagers(); j++ ) {
                fprintf(FID, "\tspecials.weakperiodic{%u}.descType=%u;\n", wpbccount, wpbc->giveBasisType() );
                fprintf(FID, "\tspecials.weakperiodic{%u}.coefficients=[", wpbccount);
                for ( Dof *dof: *wpbc->giveInternalDofManager(j) ) {
                    double X = dof->giveUnknown(VM_Total, tStep);
                    fprintf(FID, "%e\t", X);
                }
 
                fprintf(FID, "];\n");
                wpbccount++;
            }
        }
        SolutionbasedShapeFunction *sbsf = dynamic_cast< SolutionbasedShapeFunction *>( gbc.get());
        if (sbsf) {
            fprintf(FID, "\tspecials.solutionbasedsf{%u}.values=[", sbsfcount);
            for ( Dof *dof: *sbsf->giveInternalDofManager(1) ) {                  // Only one internal dof manager
                double X = dof->giveUnknown(VM_Total, tStep);
                fprintf(FID, "%e\t", X);
            }
            fprintf(FID, "];\n");
            sbsfcount++;
        }
        PrescribedMean *m = dynamic_cast<PrescribedMean *> ( gbc.get() );
        if (m) {
            fprintf(FID, "\tspecials.prescribedmean{%u}.value=[", mcount);
            for ( Dof *dof: *m->giveInternalDofManager(1)) {
                double X = dof->giveUnknown(VM_Total, tStep);
                fprintf(FID, "%e\t", X);
            }
            fprintf(FID, "];\n");
            mcount++;
        }
#ifdef __SM_MODULE
        PrescribedDispSlipBCDirichletRC *pdsd = dynamic_cast<PrescribedDispSlipBCDirichletRC *>( gbc.get() );
        if (pdsd) {
            FloatArray stress, bStress, rStress;
            pdsd->computeStress(stress, tStep);
            pdsd->computeTransferStress(bStress, tStep);
            pdsd->computeReinfStress(rStress, tStep);
            fprintf(FID, "\tspecials.prescribeddispslipbcdirichletrc{%u}.stress=[", pdsdcount);
            for ( auto i : stress ) {
                fprintf(FID, "%e\t", i);
            }
            fprintf(FID, "];\n");
            fprintf(FID, "\tspecials.prescribeddispslipbcdirichletrc{%u}.transferstress=[", pdsdcount);
            for ( auto i : bStress ) {
                fprintf(FID, "%e\t", i);
            }
            fprintf(FID, "];\n");
            fprintf(FID, "\tspecials.prescribeddispslipbcdirichletrc{%u}.reinfstress=[", pdsdcount);
            for ( auto i : rStress ) {
                fprintf(FID, "%e\t", i);
            }
            fprintf(FID, "];\n");
            pdsdcount++;
        }
        PrescribedDispSlipBCNeumannRC *pdsn = dynamic_cast<PrescribedDispSlipBCNeumannRC *>( gbc.get() );
        if (pdsn) {
            FloatArray stress, bStress, rStress;
            pdsn->computeStress(stress, tStep);
            pdsn->computeTransferStress(bStress, tStep);
            pdsn->computeReinfStress(rStress, tStep);
            fprintf(FID, "\tspecials.prescribeddispslipbcneumannrc{%u}.stress=[", pdsncount);
            for ( auto i : stress ) {
                fprintf(FID, "%e\t", i);
            }
            fprintf(FID, "];\n");
            fprintf(FID, "\tspecials.prescribeddispslipbcneumannrc{%u}.transferstress=[", pdsncount);
            for ( auto i : bStress ) {
                fprintf(FID, "%e\t", i);
            }
            fprintf(FID, "];\n");
            fprintf(FID, "\tspecials.prescribeddispslipbcneumannrc{%u}.reinfstress=[", pdsncount);
            for ( auto i : rStress ) {
                fprintf(FID, "%e\t", i);
            }
            fprintf(FID, "];\n");
            pdsncount++;
        }
        PrescribedDispSlipMultiple *pdsm = dynamic_cast<PrescribedDispSlipMultiple *>( gbc.get() );
        if (pdsm) {
            FloatArray stress, bStress, rStress;
            pdsm->computeStress(stress, tStep);
            pdsm->computeTransferStress(bStress, tStep);
            pdsm->computeReinfStress(rStress, tStep);
            fprintf(FID, "\tspecials.prescribeddispslipmultiple{%u}.stress=[", pdsncount);
            for ( auto i : stress ) {
                fprintf(FID, "%e\t", i);
            }
            fprintf(FID, "];\n");
            fprintf(FID, "\tspecials.prescribeddispslipmultiple{%u}.transferstress=[", pdsncount);
            for ( auto i : bStress ) {
                fprintf(FID, "%e\t", i);
            }
            fprintf(FID, "];\n");
            fprintf(FID, "\tspecials.prescribeddispslipmultiple{%u}.reinfstress=[", pdsncount);
            for ( auto i : rStress ) {
                fprintf(FID, "%e\t", i);
            }
            fprintf(FID, "];\n");
            pdsmcount++;
        }
        TransverseReinfConstraint *trc = dynamic_cast<TransverseReinfConstraint *> ( gbc.get() );
        if (trc) {
            FloatArray lambda;
            trc->computeField(lambda, tStep);
            fprintf(FID, "\tspecials.transversereinfconstraint{%u}.stress=[", trccount);
            for ( auto i : lambda ) {
                fprintf(FID, "%e\t", i);
            }
            fprintf(FID, "];\n");
            trccount++;
        }
#endif
    }
}
 
 
void
MatlabExportModule :: doOutputReactionForces(TimeStep *tStep,    FILE *FID)
{
 
    int domainIndex = 1;
    Domain *domain  = emodel->giveDomain( domainIndex );
 
    FloatArray reactions;
    IntArray dofManMap, dofidMap, eqnMap;
#ifdef __SM_MODULE
    StructuralEngngModel *strEngMod = dynamic_cast< StructuralEngngModel * >(emodel);
    if ( strEngMod ) {
        strEngMod->buildReactionTable(dofManMap, dofidMap, eqnMap, tStep, domainIndex);
        strEngMod->computeReaction(reactions, tStep, 1);
    } else
#endif
    {
        OOFEM_ERROR("Cannot export reaction forces - only implemented for structural problems.");
    }
 
    // Set the nodes and elements to export based on sets
    if ( this->reactionForcesNodeSet > 0 ) {
        Set *set = domain->giveSet( this->reactionForcesNodeSet );
        reactionForcesDofManList = set->giveNodeList();
    }
 
 
    int numDofManToExport = this->reactionForcesDofManList.giveSize();
    if ( numDofManToExport == 0 ) { // No dofMan's given - export every dMan with reaction forces
 
        for (int i = 1; i <= domain->giveNumberOfDofManagers(); i++) {
            if ( dofManMap.contains(i) ) {
                this->reactionForcesDofManList.followedBy(i);
            }
        }
        numDofManToExport = this->reactionForcesDofManList.giveSize();
    }
 
 
    // Output header
    fprintf( FID, "\n %%%% Export of reaction forces \n\n" );
 
    // Output the dofMan numbers that are exported
    fprintf( FID, "\tReactionForces.DofManNumbers = [" );
    for ( int i = 1; i <= numDofManToExport; i++ ) {
        fprintf( FID, "%i ", this->reactionForcesDofManList.at(i) );
    }
    fprintf( FID, "];\n" );
 
 
    // Define the reaction forces as a cell object
    fprintf( FID, "\tReactionForces.ReactionForces = cell(%i,1); \n", numDofManToExport );
    fprintf( FID, "\tReactionForces.DofIDs = cell(%i,1); \n", numDofManToExport );
 
 
    // Output the reaction forces for each dofMan. If a certain dof is not prescribed zero is exported.
    IntArray dofIDs;
    for ( int i = 1; i <= numDofManToExport; i++ ) {
        int dManNum = this->reactionForcesDofManList.at(i);
 
        fprintf(FID, "\tReactionForces.ReactionForces{%i} = [", i);
        if ( dofManMap.contains( dManNum ) ) {
 
            DofManager *dofMan = domain->giveDofManager( dManNum );
            dofIDs.clear();
 
            for ( Dof *dof: *dofMan ) {
                int num = dof->giveEquationNumber( EModelDefaultPrescribedEquationNumbering() );
                int pos = eqnMap.findFirstIndexOf( num );
                dofIDs.followedBy(dof->giveDofID());
                if ( pos > 0 ) {
                    fprintf(FID, "%e ", reactions.at(pos));
                } else {
                    fprintf( FID, "%e ", 0.0 ); // if not prescibed output zero
                }
            }
        }
        fprintf(FID, "];\n");
 
        // Output dof ID's
 
        fprintf( FID, "\tReactionForces.DofIDs{%i} = [", i);
        if ( dofManMap.contains( dManNum ) ) {
            for ( int id: dofIDs ) {
                fprintf( FID, "%i ", id );
            }
        }
        fprintf(FID, "];\n");
    }
            
    // Output the current load level (useful for CALM solver)
    double loadLevel = domain->giveEngngModel()->giveLoadLevel();
    fprintf( FID, "\tReactionForces.LoadLevel = [" );
    fprintf( FID, "%.9e", loadLevel);
    fprintf( FID, "];\n" );
}
 
 
void
MatlabExportModule :: doOutputIntegrationPointFields(TimeStep *tStep,    FILE *FID)
{
 
    int domainIndex = 1;
    Domain *domain  = emodel->giveDomain( domainIndex );
 
    // Output header
    fprintf( FID, "\n %%%% Export of internal variables in integration points \n\n" );
    fprintf( FID, "\n %% for interpretation of internal var. numbers see internalstatetype.h\n");
 
 
    int numVars = this->internalVarsToExport.giveSize();
    // Output the internalVarsToExport-list
    fprintf( FID, "\tIntegrationPointFields.InternalVarsToExport = [" );
    for ( int i = 1; i <= numVars; i++ ) {
        fprintf( FID, "%i ", this->internalVarsToExport.at(i) );
    }
    fprintf( FID, "];\n" );
 
 
 
 
    if ( this->IPFieldsElSet > 0 ) {
        Set *set = domain->giveSet( this->IPFieldsElSet );
        elList = set->giveElementList();
    }
 
 
    FloatArray valueArray;
 
    int nelem = this->elList.giveSize();
 
    fprintf( FID, "\tIntegrationPointFields.Elements = cell(%i,1); \n", nelem );
 
    for ( int ielem = 1; ielem <= nelem; ielem++ ) {
        Element *el = domain->giveElement( this->elList.at(ielem) );
        fprintf( FID, "\tIntegrationPointFields.Elements{%i}.elementNumber = %i; \n", ielem, el->giveNumber());
 
        int numIntRules = el->giveNumberOfIntegrationRules();
        fprintf( FID, "\tIntegrationPointFields.Elements{%i}.integrationRule = cell(%i,1); \n", ielem, numIntRules);
        for ( int i = 1; i <= numIntRules; i++ ) {
            IntegrationRule *iRule = el->giveIntegrationRule(i-1);
 
            fprintf( FID, "\tIntegrationPointFields.Elements{%i}.integrationRule{%i}.ip = cell(%i,1); \n ",
                     ielem, i, iRule->giveNumberOfIntegrationPoints() );
 
            // Loop over integration points
            for ( GaussPoint *ip: *iRule ) {
 
                double weight = ip->giveWeight();
 
                fprintf( FID, "\tIntegrationPointFields.Elements{%i}.integrationRule{%i}.ip{%i}.ipWeight = %e; \n ",
                         ielem, i, ip->giveNumber(), weight);
 
 
                // export Gauss point coordinates
                fprintf( FID, "\tIntegrationPointFields.Elements{%i}.integrationRule{%i}.ip{%i}.coords = [",
                         ielem, i, ip->giveNumber());
 
                FloatArray coords;
                el->computeGlobalCoordinates( coords, ip->giveNaturalCoordinates() );
                for ( int ic = 1; ic <= coords.giveSize(); ic++ ) {
                    fprintf( FID, "%e ", coords.at(ic) );
                }
                fprintf( FID, "]; \n" );
 
                //export volume around Gauss point
                fprintf( FID, "\tIntegrationPointFields.Elements{%i}.integrationRule{%i}.ip{%i}.volume = %e; \n ",
                         ielem, i, ip->giveNumber(), el->computeVolumeAround(ip));
 
                // export internal variables
                fprintf( FID, "\tIntegrationPointFields.Elements{%i}.integrationRule{%i}.ip{%i}.valArray = cell(%i,1); \n",
                         ielem, i, ip->giveNumber(), numVars);
 
                for ( int iv = 1; iv <= numVars; iv++ ) {
                    fprintf( FID, "\tIntegrationPointFields.Elements{%i}.integrationRule{%i}.ip{%i}.valArray{%i} = [",
                             ielem, i, ip->giveNumber(), iv);
                    InternalStateType vartype = ( InternalStateType ) this->internalVarsToExport.at(iv);
                    el->giveIPValue(valueArray, ip, vartype, tStep);
                    int nv = valueArray.giveSize();
                    for ( int ic = 1; ic <= nv; ic++ ) {
                        fprintf( FID, "%.6e ", valueArray.at(ic) );
                    }
                    fprintf( FID, "]; \n" );
                }
            }
 
        }
    }
 
}
 
 
void
MatlabExportModule :: initialize()
{ 
    ExportModule :: initialize();
}
 
 
void
MatlabExportModule :: terminate()
{ }
 
 
FILE *
MatlabExportModule :: giveOutputStream(TimeStep *tStep)
{
    FILE *answer;
 
    char fext[100];
    sprintf( fext, "_m%d_%d", this->number, tStep->giveNumber() );
 
    if ( this->testSubStepOutput() ) {
        // include tStep version in output file name
        if ( this->emodel->isParallel() && this->emodel->giveNumberOfProcesses() > 1 ) {
            sprintf( fext, "_%03d_m%d_%d_%d", emodel->giveRank(), this->number, tStep->giveNumber(), tStep->giveSubStepNumber() );
        } else {
            sprintf( fext, "_m%d_%d_%d", this->number, tStep->giveNumber(), tStep->giveSubStepNumber() );
        }
    } else   {
        if ( this->emodel->isParallel() && this->emodel->giveNumberOfProcesses() > 1 ) {
            sprintf( fext, "_%03d_m%d_%d", emodel->giveRank(), this->number, tStep->giveNumber() );
        } else {
            sprintf( fext, "_m%d_%d", this->number, tStep->giveNumber() );
        }
    }
 
    std :: string fileName;
 
    fileName = this->emodel->giveOutputBaseFileName();
 
    size_t foundDot;
    foundDot = fileName.rfind(".");
 
    // CARL
    while (foundDot != std :: string :: npos) {
        fileName.replace(foundDot, 1, "_");
        foundDot = fileName.rfind(".");
//        printf("%s\n", fileName.c_str());
    }
 
    fileName += fext;
 
    std :: string temp;
    temp = fileName;
    size_t backslash = temp.rfind("/");
 
    if (backslash != std :: string :: npos ) {
        functionname = temp.substr(backslash+1, std :: string :: npos);
    } else {
        functionname = temp;
    }
 
    fileName += ".m";
 
    if ( ( answer = fopen(fileName.c_str(), "w") ) == NULL ) {
        OOFEM_ERROR("failed to open file %s", fileName.c_str() );
    }
 
    return answer;
}
 
void
MatlabExportModule :: doOutputHomogenizeDofIDs(TimeStep *tStep, FILE *FID)
{
    std :: vector <FloatArray> HomQuantities;
    double vol = 0.0;
 
    // Initialize vector of arrays constaining homogenized quantities
    HomQuantities.resize(internalVarsToExport.giveSize());
 
    int nelem = this->elList.giveSize();
    for (int i = 1; i<=nelem; i++) {
        Element *e = this->emodel->giveDomain(1)->giveElement(elList.at(i));
        FEInterpolation *Interpolation = e->giveInterpolation();
 
        vol += e->computeVolumeAreaOrLength();
 
        for ( auto &gp: *e->giveDefaultIntegrationRulePtr() ) {
 
            for (int j=0; j<internalVarsToExport.giveSize(); j++) {
                FloatArray elementValues;
                e->giveIPValue(elementValues, gp, (InternalStateType) internalVarsToExport[j], tStep);
                double detJ=fabs(Interpolation->giveTransformationJacobian( gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(e)));
 
                elementValues.times(gp->giveWeight()*detJ);
                if (HomQuantities[j].giveSize() == 0) {
                    HomQuantities[j].resize(elementValues.giveSize());
                };
                HomQuantities[j].add(elementValues);
            }
        }
    }
 
    if (noscaling) vol = 1.0;
 
    for ( std :: size_t i = 0; i < HomQuantities.size(); i ++) {
        FloatArray &thisIS = HomQuantities[i];
        thisIS.times(1.0/vol);
        fprintf(FID, "\tspecials.%s = [", __InternalStateTypeToString ( (InternalStateType) internalVarsToExport[i] ) );
 
        for (int j = 0; j<thisIS.giveSize(); j++) {
            fprintf(FID, "%e", thisIS.at(j+1));
            if (j!=(thisIS.giveSize()-1) ) {
                fprintf(FID, ", ");
            }
        }
        fprintf(FID, "];\n");
    }
 
}
 
} // end namespace oofem