AssemblyMixin.h

Go to the documentation of this file.

 
#pragma once
 
#include "MathLib/LinAlg/LinAlg.h"
#include "NumLib/DOF/GlobalMatrixProviders.h"
#include "ProcessLib/Assembly/ParallelVectorMatrixAssembler.h"
#include "ProcessLib/ProcessVariable.h"
#include "VectorMatrixAssembler.h"
 
namespace ProcessLib
{
struct SubmeshAssemblyData
{
    explicit SubmeshAssemblyData(
        MeshLib::Mesh const& mesh,
        std::vector<std::reference_wrapper<MeshLib::PropertyVector<double>>>&&
            residuum_vectors);
 
    MeshLib::PropertyVector<std::size_t> const& bulk_element_ids;
    MeshLib::PropertyVector<std::size_t> const& bulk_node_ids;
    std::vector<std::size_t> active_element_ids;
    std::vector<std::reference_wrapper<MeshLib::PropertyVector<double>>>
        residuum_vectors;
};
 
class AssemblyMixinBase
{
    enum class ActiveElementIDsState
    {
        UNINITIALIZED,
        HAS_DEACTIVATED_SUBDOMAINS,
        NO_DEACTIVATED_SUBDOMAINS
    };
 
protected:
    explicit AssemblyMixinBase(AbstractJacobianAssembler& jacobian_assembler)
        : pvma_{jacobian_assembler} {};
 
    void initializeAssemblyOnSubmeshes(
        const int process_id,
        MeshLib::Mesh& bulk_mesh,
        std::vector<std::reference_wrapper<MeshLib::Mesh>> const& submeshes,
        std::vector<std::string> const& residuum_names,
        std::vector<std::reference_wrapper<ProcessVariable>> const& pvs);
 
    void updateActiveElements(ProcessLib::ProcessVariable const& pv);
 
    static void copyResiduumVectorsToBulkMesh(
        GlobalVector const& rhs,
        NumLib::LocalToGlobalIndexMap const& local_to_global_index_map,
        std::vector<std::reference_wrapper<MeshLib::PropertyVector<double>>>
            residuum_vectors);
 
    static void copyResiduumVectorsToSubmesh(
        GlobalVector const& rhs,
        NumLib::LocalToGlobalIndexMap const& local_to_global_index_map,
        SubmeshAssemblyData const& sad);
 
private:
    void updateActiveElementsImpl(ProcessLib::ProcessVariable const& pv);
 
protected:
    std::vector<SubmeshAssemblyData> submesh_assembly_data_;
    std::vector<std::reference_wrapper<MeshLib::PropertyVector<double>>>
        residuum_vectors_bulk_;
 
    std::size_t b_submesh_id_ = 0;
 
    Assembly::ParallelVectorMatrixAssembler pvma_;
 
private:
    ActiveElementIDsState ids_state_ = ActiveElementIDsState::UNINITIALIZED;
};
 
template <typename Process>
class AssemblyMixin : private AssemblyMixinBase
{
    // Enforce correct use of CRTP, i.e., that Process is derived from
    // AssemblyMixin<Process>.
    using AssemblyMixinBase::AssemblyMixinBase;
    friend Process;
 
public:
    void initializeAssemblyOnSubmeshes(
        const int process_id,
        std::vector<std::reference_wrapper<MeshLib::Mesh>> const& submeshes,
        std::vector<std::string> const& residuum_names)
    {
        AssemblyMixinBase::initializeAssemblyOnSubmeshes(
            process_id, derived().getMesh(), submeshes, residuum_names,
            derived().getProcessVariables(process_id));
    }
 
    void updateActiveElements(const int process_id)
    {
        // convention: process variable 0 governs where assembly takes place
        // (active element IDs)
        ProcessLib::ProcessVariable const& pv =
            derived().getProcessVariables(process_id)[0];
 
        AssemblyMixinBase::updateActiveElements(pv);
    }
 
    // cppcheck-suppress functionStatic
    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*/)
    {
        /*
        DBUG("AssemblyMixin assemble(t={}, dt={}, process_id={}).", t, dt,
             process_id);
 
        assembleGeneric(&Assembly::ParallelVectorMatrixAssembler::assemble, t,
        dt, x, x_prev, process_id, M, K, b);
        */
        OGS_FATAL("Not yet implemented.");
    }
 
    void assembleWithJacobian(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,
                              GlobalMatrix& Jac)
    {
        DBUG("AssemblyMixin assembleWithJacobian(t={}, dt={}, process_id={}).",
             t, dt, process_id);
 
        assembleGeneric(
            &Assembly::ParallelVectorMatrixAssembler::assembleWithJacobian, t,
            dt, x, x_prev, process_id, M, K, b, Jac);
    }
 
private:
    Process& derived() { return static_cast<Process&>(*this); }
    Process const& derived() const
    {
        return static_cast<Process const&>(*this);
    }
 
    template <typename Method, typename... Jac>
    void assembleGeneric(Method global_assembler_method, 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, Jac&... jac_or_not_jac)
    {
        // TODO why not getDOFTables(x.size()); ?
        std::vector<NumLib::LocalToGlobalIndexMap const*> const dof_tables{
            derived()._local_to_global_index_map.get()};
 
        auto const& loc_asms = derived().local_assemblers_;
 
        if (!submesh_assembly_data_.empty())
        {
            auto& b_submesh = NumLib::GlobalVectorProvider::provider.getVector(
                b, b_submesh_id_);
 
            for (auto const& sad : submesh_assembly_data_)
            {
                b_submesh.setZero();
 
                (pvma_.*global_assembler_method)(
                    loc_asms, sad.active_element_ids, dof_tables, t, dt, x,
                    x_prev, process_id, M, K, b_submesh, jac_or_not_jac...);
 
                MathLib::LinAlg::axpy(b, 1.0, b_submesh);
 
                AssemblyMixinBase::copyResiduumVectorsToSubmesh(
                    b_submesh, *(dof_tables.front()), sad);
            }
 
            NumLib::GlobalVectorProvider::provider.releaseVector(b_submesh);
        }
        else
        {
            // convention: process variable 0 governs where assembly takes
            // place (active element IDs)
            ProcessLib::ProcessVariable const& pv =
                derived().getProcessVariables(process_id)[0];
 
            (pvma_.*global_assembler_method)(
                loc_asms, pv.getActiveElementIDs(), dof_tables, t, dt, x,
                x_prev, process_id, M, K, b, jac_or_not_jac...);
        }
 
        AssemblyMixinBase::copyResiduumVectorsToBulkMesh(
            b, *(dof_tables.front()), residuum_vectors_bulk_);
    }
};
 
}  // namespace ProcessLib