vtktpsamodule.hpp

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
  Copyright 2025 NORCE AS
 
  This file is part of the Open Porous Media project (OPM).
 
  OPM is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 2 of the License, or
  (at your option) any later version.
 
  OPM 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 General Public License for more details.
 
  You should have received a copy of the GNU General Public License
  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 
  Consult the COPYING file in the top-level source directory of this
  module for the precise wording of the license and the list of
  copyright holders.
*/
#ifndef VTK_TPSA_MODULE_HPP
#define VTK_TPSA_MODULE_HPP
 
#include <dune/common/dynmatrix.hh>
#include <dune/common/fvector.hh>
 
#include <opm/common/utility/SymmTensor.hpp>
#include <opm/common/utility/VoigtArray.hpp>
#include <opm/material/densead/Math.hpp>
 
#include <opm/models/io/baseoutputmodule.hh>
#include <opm/models/io/vtkmultiwriter.hh>
#include <opm/models/io/vtktpsaparams.hpp>
#include <opm/models/tpsa/tpsabaseproperties.hpp>
#include <opm/models/utils/parametersystem.hpp>
#include <opm/models/utils/propertysystem.hh>
 
namespace Opm {
 
template <class TypeTag>
class VtkTpsaModule : public BaseOutputModule<TypeTag>
{
    using ParentType = BaseOutputModule<TypeTag>;
 
    using GridView = GetPropType<TypeTag, Properties::GridView>;
    using Simulator = GetPropType<TypeTag, Properties::Simulator>;
    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
 
    static constexpr auto vtkFormat = getPropValue<TypeTag, Properties::VtkOutputFormat>();
    using VtkMultiWriter = Opm::VtkMultiWriter<GridView, vtkFormat>;
 
    enum { enableMech = getPropValue<TypeTag, Properties::EnableMech>() };
 
    using BufferType = typename ParentType::BufferType;
    using ScalarBuffer = typename ParentType::ScalarBuffer;
    using VectorBuffer = typename ParentType::VectorBuffer;
    using TensorBuffer = typename ParentType::TensorBuffer;
 
    using SymTensor = SymmTensor<Scalar>;
    using Tensor = Dune::DynamicMatrix<Scalar>;
 
public:
    explicit VtkTpsaModule(const Simulator& simulator)
        : ParentType(simulator)
    {
        // Read runtime parameters
        if constexpr(enableMech) {
            params_.read();
        }
    }
 
    static void registerParameters()
    {
        if constexpr(enableMech) {
            VtkTpsaParams::registerParameters();
        }
    }
 
    void allocBuffers() override
    {
        // Check if vtk output is enabled
        if (!Parameters::Get<Parameters::EnableVtkOutput>()) {
            return;
        }
 
        // Allocate buffers
        if constexpr(enableMech) {
            // Displacement
            if (params_.displacementOutput_) {
                this->resizeVectorBuffer_(displacement_, BufferType::Dof);
            }
 
            // Rotation
            if (params_.rotationOutput_) {
                this->resizeVectorBuffer_(rotation_, BufferType::Dof);
            }
 
            // Solid pressure
            if (params_.solidPressureOutput_) {
                this->resizeScalarBuffer_(solidPres_, BufferType::Dof);
            }
 
            // Potential pressure force
            if (params_.potPresForceOutput_) {
                this->resizeScalarBuffer_(potPresForce_, BufferType::Dof);
            }
 
            // Stress
            if (params_.stressOutput_) {
                this->resizeTensorBuffer_(stress_, BufferType::Dof);
            }
 
            // Stress (without potential fracture contributions)
            if (params_.delStressOutput_) {
                this->resizeTensorBuffer_(delStress_, BufferType::Dof);
            }
 
            // Stress (without potential pressure/temperature/fracture forces)
            if (params_.linStressOutput_) {
                this->resizeTensorBuffer_(linStress_, BufferType::Dof);
            }
 
            // Strain
            if (params_.strainOutput_) {
                this->resizeTensorBuffer_(strain_, BufferType::Dof);
            }
        }
    }
 
    void processElement(const ElementContext& elemCtx) override
    {
        // Check if vtk output is enabled
        if (!Parameters::Get<Parameters::EnableVtkOutput>()) {
            return;
        }
 
        if constexpr(enableMech) {
            // Assign quantities from material state
            const auto& problem = elemCtx.problem();
            const auto& geoMechModel = problem.geoMechModel();
            const auto& faceStressInfo = geoMechModel.linearizer().getStressInfo();
            for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
                const unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
                const auto& materialState = geoMechModel.materialState(globalDofIdx, /*timeIdx=*/0);
 
                // Loop over x-, y- and z-dir (corresponding to dirIdx = 0, 1, 2)
                for (int dirIdx = 0; dirIdx < 3; ++dirIdx) {
                    // Displacement
                    if (params_.displacementOutput_) {
                        displacement_[globalDofIdx][dirIdx] = scalarValue(materialState.displacement(dirIdx));
                    }
 
                    // Rotation
                    if (params_.rotationOutput_) {
                        rotation_[globalDofIdx][dirIdx] = scalarValue(materialState.rotation(dirIdx));
                    }
                }
 
                // Solid pressure
                if (params_.solidPressureOutput_) {
                    solidPres_[globalDofIdx] = scalarValue(materialState.solidPressure());
                }
 
                // Potential pressure force
                if (params_.potPresForceOutput_) {
                    potPresForce_[globalDofIdx] =
                        geoMechModel.mechPotentialPressForce(globalDofIdx);
                }
 
                // Stress
                if (params_.stressOutput_ && !faceStressInfo.empty()) {
                    SymTensor stressSymTensor = geoMechModel.stress(globalDofIdx, true);
                    Tensor& stressTensor = stress_[globalDofIdx];
                    setTensorFromVoigt_(stressTensor, stressSymTensor);
                }
 
                // Stress (without potential fracture contributions)
                if (params_.delStressOutput_ && !faceStressInfo.empty()) {
                    SymTensor delStressSymTensor = geoMechModel.delstress(globalDofIdx);
                    Tensor& delStressTensor = delStress_[globalDofIdx];
                    setTensorFromVoigt_(delStressTensor, delStressSymTensor);
                }
 
                // Stress (without potential pressure/temperature/fracture forces)
                if (params_.linStressOutput_ && !faceStressInfo.empty()) {
                    SymTensor linStressSymTensor = geoMechModel.linstress(globalDofIdx);
                    Tensor& linStressTensor = linStress_[globalDofIdx];
                    setTensorFromVoigt_(linStressTensor, linStressSymTensor);
                }
 
                // Strain
                if (params_.strainOutput_ && !faceStressInfo.empty()) {
                    SymTensor strainSymTensor = geoMechModel.strain(globalDofIdx, true);
                    Tensor& strainTensor = strain_[globalDofIdx];
                    setTensorFromVoigt_(strainTensor, strainSymTensor);
                }
            }
        }
    }
 
    void commitBuffers(BaseOutputWriter& baseWriter) override
    {
        // Check if writer exists or mechanics output is enabled
        if (!dynamic_cast<VtkMultiWriter*>(&baseWriter)) {
            return;
        }
 
        if constexpr(enableMech) {
            // Displacement
            if (params_.displacementOutput_) {
                this->commitVectorBuffer_(baseWriter, "displacement", displacement_, BufferType::Dof);
            }
 
            // Rotation
            if (params_.rotationOutput_) {
                this->commitVectorBuffer_(baseWriter, "rotation", rotation_, BufferType::Dof);
            }
 
            // Solid pressure
            if (params_.solidPressureOutput_) {
                this->commitScalarBuffer_(baseWriter, "solid_pressure", solidPres_, BufferType::Dof);
            }
 
            if (params_.potPresForceOutput_) {
                this->commitScalarBuffer_(baseWriter,
                                          "potentialPressureForce",
                                          potPresForce_,
                                          BufferType::Dof);
            }
 
            if (params_.stressOutput_) {
                this->commitTensorBuffer_(baseWriter, "stress", stress_, BufferType::Dof);
            }
 
            if (params_.delStressOutput_) {
                this->commitTensorBuffer_(baseWriter, "delStress", delStress_, BufferType::Dof);
            }
 
            if (params_.linStressOutput_) {
                this->commitTensorBuffer_(baseWriter, "linStress", linStress_, BufferType::Dof);
            }
 
            if (params_.strainOutput_) {
                this->commitTensorBuffer_(baseWriter, "strain", strain_, BufferType::Dof);
            }
        }
    }
 
private:
    static void setTensorFromVoigt_(Tensor& tensor, const SymTensor& symTensor)
    {
        // Diagonal terms
        constexpr auto& ind = SymTensor::diag_indices;
        for (std::size_t i = 0; i < 3; ++i) {
            tensor[i][i] = symTensor[ind[i]];
        }
 
        // Off-diagonal terms
        tensor[0][1] = symTensor[VoigtIndex::XY];
        tensor[0][2] = symTensor[VoigtIndex::XZ];
        tensor[1][2] = symTensor[VoigtIndex::YZ];
        for (std::size_t i = 0; i < 3; ++i) {
            for (std::size_t j = 0; j < 3; ++j) {
                if (i > j) {
                    tensor[i][j] = tensor[j][i];
                }
            }
        }
    }
 
    VtkTpsaParams params_{};
 
    VectorBuffer displacement_{};
    VectorBuffer rotation_{};
 
    ScalarBuffer solidPres_{};
    ScalarBuffer potPresForce_{};
 
    TensorBuffer stress_{};
    TensorBuffer delStress_{};
    TensorBuffer linStress_{};
    TensorBuffer strain_{};
}; // class VtkTpsaModule
 
} // namespace Opm
 
#endif