SecondaryVariable.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 "BaseLib/Algorithm.h"
#include "NumLib/Extrapolation/ExtrapolatableElementCollection.h"
#include "NumLib/Extrapolation/Extrapolator.h"
#include "ProcessLib/Utils/TransposeInPlace.h"
 
namespace NumLib
{
class LocalToGlobalIndexMap;
}
 
namespace ProcessLib
{
struct SecondaryVariableFunctions final
{
    using Function = std::function<GlobalVector const&(
        const double t,
        std::vector<GlobalVector*> const& x,
        std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
        std::unique_ptr<GlobalVector>& result_cache)>;
 
    template <typename F1, typename F2>
    SecondaryVariableFunctions(const unsigned num_components_, F1&& eval_field_,
                               F2&& eval_residuals_)
        : num_components(num_components_),
          eval_field(std::forward<F1>(eval_field_)),
          eval_residuals(std::forward<F2>(eval_residuals_))
    {
        // Used to detect nasty implicit conversions.
        static_assert(
            std::is_same_v<
                GlobalVector const&,
                typename std::invoke_result_t<
                    F1, double const, std::vector<GlobalVector*> const&,
                    std::vector<NumLib::LocalToGlobalIndexMap const*> const&,
                    std::unique_ptr<GlobalVector>&>>,
            "The function eval_field_ does not return a const reference"
            " to a GlobalVector");
 
        static_assert(
            std::is_same_v<
                GlobalVector const&,
                typename std::invoke_result_t<
                    F2, double const, std::vector<GlobalVector*> const&,
                    std::vector<NumLib::LocalToGlobalIndexMap const*> const&,
                    std::unique_ptr<GlobalVector>&>>,
            "The function eval_residuals_ does not return a const reference"
            " to a GlobalVector");
    }
 
    template <typename F1>
    SecondaryVariableFunctions(const unsigned num_components_, F1&& eval_field_,
                               std::nullptr_t)
        : num_components(num_components_),
          eval_field(std::forward<F1>(eval_field_))
    {
        // Used to detect nasty implicit conversions.
        static_assert(
            std::is_same_v<
                GlobalVector const&,
                typename std::invoke_result_t<
                    F1, double const, std::vector<GlobalVector*> const&,
                    std::vector<NumLib::LocalToGlobalIndexMap const*> const&,
                    std::unique_ptr<GlobalVector>&>>,
            "The function eval_field_ does not return a const reference"
            " to a GlobalVector");
    }
 
    const unsigned num_components;  

    Function const eval_field;

    Function const eval_residuals;
};

struct SecondaryVariable final
{
    std::string const name;  

    SecondaryVariableFunctions fcts;
};

class SecondaryVariableCollection final
{
public:
    void addNameMapping(std::string const& internal_name,
                        std::string const& external_name);

    void addSecondaryVariable(std::string const& internal_name,
                              SecondaryVariableFunctions&& fcts);

    SecondaryVariable const& get(std::string const& external_name) const;
 
    std::map<std::string, std::string>::const_iterator begin() const;
    std::map<std::string, std::string>::const_iterator end() const;
 
private:
    std::map<std::string, std::string> _map_external_to_internal;

    std::map<std::string, SecondaryVariable> _configured_secondary_variables;
};

template <typename LocalAssemblerCollection>
SecondaryVariableFunctions makeExtrapolator(
    const unsigned num_components,
    NumLib::Extrapolator& extrapolator,
    LocalAssemblerCollection const& local_assemblers,
    typename NumLib::ExtrapolatableLocalAssemblerCollection<
        LocalAssemblerCollection>::IntegrationPointValuesMethod
        integration_point_values_method)
{
    auto const eval_field =
        [num_components, &extrapolator, &local_assemblers,
         integration_point_values_method](
            const double t,
            std::vector<GlobalVector*> const& x,
            std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
            std::unique_ptr<GlobalVector>& /*result_cache*/
            ) -> GlobalVector const&
    {
        auto const extrapolatables = NumLib::makeExtrapolatable(
            local_assemblers, integration_point_values_method);
        extrapolator.extrapolate(num_components, extrapolatables, t, x,
                                 dof_tables);
        return extrapolator.getNodalValues();
    };
 
    auto const eval_residuals =
        [num_components, &extrapolator, &local_assemblers,
         integration_point_values_method](
            const double t,
            std::vector<GlobalVector*> const& x,
            std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
            std::unique_ptr<GlobalVector>& /*result_cache*/
            ) -> GlobalVector const&
    {
        auto const extrapolatables = NumLib::makeExtrapolatable(
            local_assemblers, integration_point_values_method);
        extrapolator.calculateResiduals(num_components, extrapolatables, t, x,
                                        dof_tables);
        return extrapolator.getElementResiduals();
    };
    return {num_components, eval_field, eval_residuals};
}

template <typename LocalAssemblerCollection, typename IPDataAccessor>
SecondaryVariableFunctions makeExtrapolator2(
    const unsigned num_components,
    NumLib::Extrapolator& extrapolator,
    LocalAssemblerCollection const& local_assemblers,
    IPDataAccessor&& accessor)
{
    using LocalAssemblerInterface = std::remove_cvref_t<
        decltype(*std::declval<LocalAssemblerCollection>()[0])>;
    static_assert(std::is_invocable_r_v<std::vector<double>, IPDataAccessor,
                                        LocalAssemblerInterface const&>);
 
    if (num_components == 1)
    {
        auto method_wrapped =
            [accessor](
                LocalAssemblerInterface const& loc_asm, const double /*t*/,
                std::vector<GlobalVector*> const& /*x*/,
                std::vector<
                    NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
                std::vector<double>& cache) -> std::vector<double> const&
        {
            cache = accessor(loc_asm);
            return cache;
        };
 
        return makeExtrapolator(num_components, extrapolator, local_assemblers,
                                method_wrapped);
    }
 
    auto method_wrapped =
        [accessor, num_components](
            LocalAssemblerInterface const& loc_asm, const double /*t*/,
            std::vector<GlobalVector*> const& /*x*/,
            std::vector<
                NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
            std::vector<double>& cache) -> std::vector<double> const&
    {
        cache = accessor(loc_asm);
        transposeInPlace(cache, cache.size() / num_components);
        return cache;
    };
 
    return makeExtrapolator(num_components, extrapolator, local_assemblers,
                            method_wrapped);
}
 
}  // namespace ProcessLib