Process.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
// SPDX-License-Identifier: BSD-3-Clause
 
#pragma once
 
#include <string>
#include <tuple>
 
#include "AbstractJacobianAssembler.h"
#include "MaterialLib/MPL/Medium.h"
#include "MathLib/LinAlg/GlobalMatrixVectorTypes.h"
#include "MeshLib/Utils/IntegrationPointWriter.h"
#include "NumLib/ODESolver/ODESystem.h"
#include "ParameterLib/Parameter.h"
#include "ProcessLib/BoundaryConditionAndSourceTerm/BoundaryConditionCollection.h"
#include "ProcessLib/BoundaryConditionAndSourceTerm/SourceTermCollection.h"
#include "ProcessLib/Output/CellAverageData.h"
#include "ProcessLib/Output/ExtrapolatorData.h"
#include "ProcessLib/Output/SecondaryVariable.h"
#include "ProcessVariable.h"
#include "SubmeshAssemblySupport.h"
#include "VectorMatrixAssembler.h"
#include "processlib_export.h"
 
namespace MeshLib
{
class Mesh;
}
 
namespace ProcessLib
{
 
class Process
    : public NumLib::ODESystem<  // TODO: later on use a simpler ODE system
          NumLib::ODESystemTag::FirstOrderImplicitQuasilinear,
          NumLib::NonlinearSolverTag::Newton>,
      public SubmeshAssemblySupport
{
public:
    PROCESSLIB_EXPORT static const std::string constant_one_parameter_name;
 
    Process(std::string name_,
            MeshLib::Mesh& mesh,
            std::unique_ptr<AbstractJacobianAssembler>&& jacobian_assembler,
            std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const&
                parameters,
            unsigned const integration_order,
            std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
                process_variables,
            SecondaryVariableCollection&& secondary_variables,
            const bool use_monolithic_scheme = true);

    void preTimestep(std::vector<GlobalVector*> const& x, const double t,
                     const double delta_t, const int process_id);

    void postTimestep(std::vector<GlobalVector*> const& x,
                      std::vector<GlobalVector*> const& x_prev, const double t,
                      const double delta_t, int const process_id);

    void postNonLinearSolver(std::vector<GlobalVector*> const& x,
                             std::vector<GlobalVector*> const& x_prev,
                             const double t, double const dt,
                             int const process_id);
 
    void preIteration(const unsigned iter, GlobalVector const& x) final;

    void computeSecondaryVariable(double const t,
                                  double const dt,
                                  std::vector<GlobalVector*> const& x,
                                  GlobalVector const& x_prev,
                                  int const process_id);
 
    NumLib::IterationResult postIteration(GlobalVector const& x) final;
 
    void initialize(
        std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const&
            media);
 
    void setInitialConditions(
        std::vector<GlobalVector*>& process_solutions,
        std::vector<GlobalVector*> const& process_solutions_prev,
        double const t,
        int const process_id);
 
    MathLib::MatrixSpecifications getMatrixSpecifications(
        const int process_id) const override;
 
    void updateDeactivatedSubdomains(double const time, const int process_id);
 
    virtual bool isMonolithicSchemeUsed() const
    {
        return _use_monolithic_scheme;
    }
 
    virtual void extrapolateIntegrationPointValuesToNodes(
        const double /*t*/,
        std::vector<GlobalVector*> const& /*integration_point_values_vectors*/,
        std::vector<GlobalVector*>& /*nodal_values_vectors*/)
    {
    }
 
    void preAssemble(const double t, double const dt,
                     GlobalVector const& x) final;
 
    void assemble(const double t, double const dt,
                  std::vector<GlobalVector*> const& x,
                  std::vector<GlobalVector*> const& x_prev,
                  int const process_id, GlobalMatrix& M, GlobalMatrix& K,
                  GlobalVector& b) final;
 
    void assembleWithJacobian(const double t, double const dt,
                              std::vector<GlobalVector*> const& x,
                              std::vector<GlobalVector*> const& x_prev,
                              int const process_id, GlobalVector& b,
                              GlobalMatrix& Jac) final;
 
    /* Computes data necessary for output.
     *
     * \pre Must be called before postTimestep() since the latter will overwrite
     * previous states with current ones possibly affecting the data computed by
     * this method.
     */
    void preOutput(const double t, double const dt,
                   std::vector<GlobalVector*> const& x,
                   std::vector<GlobalVector*> const& x_prev,
                   int const process_id);
 
    std::vector<NumLib::IndexValueVector<GlobalIndexType>> const*
    getKnownSolutions(double const t, GlobalVector const& x,
                      int const process_id) const final
    {
        return _boundary_conditions[process_id].getKnownSolutions(t, x);
    }
 
    virtual NumLib::LocalToGlobalIndexMap const& getDOFTable(
        const int /*process_id*/) const
    {
        return *_local_to_global_index_map;
    }
 
    MeshLib::Mesh& getMesh() const { return _mesh; }
 
    std::vector<std::vector<std::reference_wrapper<ProcessVariable>>> const&
    getProcessVariables() const
    {
        return _process_variables;
    }
 
    std::vector<std::reference_wrapper<ProcessVariable>> const&
    getProcessVariables(const int process_id) const
    {
        return _process_variables[process_id];
    }
 
    std::vector<std::size_t> const& getActiveElementIDs() const
    {
        return _ids_of_active_elements;
    }
 
    SecondaryVariableCollection const& getSecondaryVariables() const
    {
        return _secondary_variables;
    }
 
    std::vector<std::unique_ptr<MeshLib::IntegrationPointWriter>> const&
    getIntegrationPointWriters() const
    {
        return _integration_point_writer;
    }
 
    // Used as a call back for CalculateSurfaceFlux process.
 
    virtual Eigen::Vector3d getFlux(
        std::size_t /*element_id*/,
        MathLib::Point3d const& /*p*/,
        double const /*t*/,
        std::vector<GlobalVector*> const& /*x*/) const
    {
        return Eigen::Vector3d{};
    }
 
    virtual void solveReactionEquation(
        std::vector<GlobalVector*>& /*x*/,
        std::vector<GlobalVector*> const& /*x_prev*/, double const /*t*/,
        double const /*dt*/, NumLib::EquationSystem& /*ode_sys*/,
        int const /*process_id*/)
    {
    }
 
    bool requiresNormalization() const override { return false; }
 
protected:
    std::vector<NumLib::LocalToGlobalIndexMap const*> getDOFTables(
        int const number_of_processes) const;
 
    NumLib::Extrapolator& getExtrapolator() const
    {
        return _extrapolator_data.getExtrapolator();
    }
 
    NumLib::LocalToGlobalIndexMap const& getSingleComponentDOFTable() const
    {
        return _extrapolator_data.getDOFTable();
    }

    void initializeProcessBoundaryConditionsAndSourceTerms(
        const NumLib::LocalToGlobalIndexMap& dof_table, const int process_id,
        std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const&
            media);
 
private:
    virtual void initializeConcreteProcess(
        NumLib::LocalToGlobalIndexMap const& dof_table,
        MeshLib::Mesh const& mesh,
        unsigned const integration_order) = 0;

    virtual void initializeBoundaryConditions(
        std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const&
            media);
 
    virtual void setInitialConditionsConcreteProcess(
        std::vector<GlobalVector*>& /*x*/,
        double const /*t*/,
        int const /*process_id*/)
    {
    }
 
    virtual void preAssembleConcreteProcess(const double /*t*/,
                                            double const /*dt*/,
                                            GlobalVector const& /*x*/)
    {
    }
 
    virtual void assembleConcreteProcess(
        const double t, double const dt, std::vector<GlobalVector*> const& x,
        std::vector<GlobalVector*> const& x_prev, int const process_id,
        GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) = 0;
 
    virtual void assembleWithJacobianConcreteProcess(
        const double t, double const dt, std::vector<GlobalVector*> const& x,
        std::vector<GlobalVector*> const& x_prev, int const process_id,
        GlobalVector& b, GlobalMatrix& Jac) = 0;
 
    virtual void preTimestepConcreteProcess(
        std::vector<GlobalVector*> const& /*x*/,
        const double /*t*/,
        const double /*dt*/,
        const int /*process_id*/)
    {
    }
 
    virtual void postTimestepConcreteProcess(
        std::vector<GlobalVector*> const& /*x*/,
        std::vector<GlobalVector*> const& /*x_prev*/,
        const double /*t*/,
        const double /*dt*/,
        int const /*process_id*/)
    {
    }
 
    virtual void postNonLinearSolverConcreteProcess(
        std::vector<GlobalVector*> const& /*x*/,
        std::vector<GlobalVector*> const& /*x_prev*/, const double /*t*/,
        double const /*dt*/, int const /*process_id*/)
    {
    }
 
    virtual void preIterationConcreteProcess(const unsigned /*iter*/,
                                             GlobalVector const& /*x*/)
    {
    }
 
    virtual void computeSecondaryVariableConcrete(
        double const /*t*/,
        double const /*dt*/,
        std::vector<GlobalVector*> const& /*x*/,
        GlobalVector const& /*x_prev*/,
        int const /*process_id*/)
    {
    }
 
    virtual NumLib::IterationResult postIterationConcreteProcess(
        GlobalVector const& /*x*/)
    {
        return NumLib::IterationResult::SUCCESS;
    }
 
    virtual void preOutputConcreteProcess(
        const double /*t*/, double const /*dt*/,
        std::vector<GlobalVector*> const& /*x*/,
        std::vector<GlobalVector*> const& /*x_prev*/, int const /*process_id*/)
    {
    }
 
protected:
    virtual void constructDofTable();

    void constructMonolithicProcessDofTable();

    void constructDofTableOfSpecifiedProcessStaggeredScheme(
        const int specified_process_id);

    virtual std::tuple<NumLib::LocalToGlobalIndexMap*, bool>
    getDOFTableForExtrapolatorData() const;

    std::vector<GlobalIndexType>
    getIndicesOfResiduumWithoutInitialCompensation() const override;

    void setReleaseNodalForces(GlobalVector const* r_neq,
                               int const process_id) override;
 
private:
    void initializeExtrapolator();

    void computeSparsityPattern();
 
public:
    std::string const name;
 
protected:
    MeshLib::Mesh& _mesh;
    std::unique_ptr<MeshLib::MeshSubset const> _mesh_subset_all_nodes;
 
    std::unique_ptr<NumLib::LocalToGlobalIndexMap> _local_to_global_index_map;
 
    SecondaryVariableCollection _secondary_variables;
 
    CellAverageData cell_average_data_;
 
    // TODO (CL) switch to the parallel vector matrix assembler here, once
    // feature complete, or use the assembly mixin and remove these two members.
    std::unique_ptr<ProcessLib::AbstractJacobianAssembler> _jacobian_assembler;
    VectorMatrixAssembler _global_assembler;
 
    const bool _use_monolithic_scheme;

    unsigned const _integration_order;

    std::vector<std::unique_ptr<MeshLib::IntegrationPointWriter>>
        _integration_point_writer;
 
    GlobalSparsityPattern _sparsity_pattern;
 
protected:
    std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>
        _process_variables;

    std::vector<BoundaryConditionCollection> _boundary_conditions;
 
private:
    std::vector<SourceTermCollection> _source_term_collections;
 
    ExtrapolatorData _extrapolator_data;

    std::vector<std::size_t> _ids_of_active_elements;
};
 
}  // namespace ProcessLib