DOFTableUtil.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 <Eigen/Core>
#include <range/v3/algorithm/fill.hpp>
 
#include "MathLib/LinAlg/LinAlg.h"
#include "MathLib/LinAlg/LinAlgEnums.h"
#include "NumLib/DOF/LocalToGlobalIndexMap.h"
 
namespace NumLib
{
double getNonGhostNodalValue(GlobalVector const& x,
                             MeshLib::Location const& location,
                             NumLib::LocalToGlobalIndexMap const& dof_table,
                             std::size_t const global_component_id);

double getNodalValue(GlobalVector const& x, MeshLib::Location const& location,
                     NumLib::LocalToGlobalIndexMap const& dof_table,
                     std::size_t const global_component_id);

std::vector<GlobalIndexType> getIndices(
    std::size_t const mesh_item_id,
    NumLib::LocalToGlobalIndexMap const& dof_table);

std::vector<GlobalIndexType> getIndices(
    std::size_t const variable_id,
    std::size_t const component_id,
    std::size_t const bulk_element_id,
    NumLib::LocalToGlobalIndexMap const& dof_table);

LocalToGlobalIndexMap::RowColumnIndices getRowColumnIndices(
    std::size_t const id,
    NumLib::LocalToGlobalIndexMap const& dof_table,
    std::vector<GlobalIndexType>& indices);

double norm(GlobalVector const& x, unsigned const global_component,
            MathLib::VecNormType norm_type,
            LocalToGlobalIndexMap const& dof_table);

template <typename Functor>
void transformVariableFromGlobalVector(
    GlobalVector const& input_vector, int const variable_id,
    NumLib::LocalToGlobalIndexMap const& local_to_global_index_map,
    MeshLib::PropertyVector<double>& output_vector, Functor map_function)
{
    MathLib::LinAlg::setLocalAccessibleVector(input_vector);
 
    // We fill the output with zeroes, because filling with NaN
    // "breaks" visualization, e.g. of heat or mass flow rate
    // with Taylor-Hood elements. I.e., all lower order
    // properties would have NaN on the higher order nodes.
    ranges::fill(output_vector, 0);
 
    int const n_components =
        local_to_global_index_map.getNumberOfVariableComponents(variable_id);
    for (int component = 0; component < n_components; ++component)
    {
        auto const& mesh_subset =
            local_to_global_index_map.getMeshSubset(variable_id, component);
        for (auto const& l : MeshLib::views::meshLocations(
                 mesh_subset, MeshLib::MeshItemType::Node))
        {
            auto const node_id = l.item_id;
            auto const input_index = local_to_global_index_map.getGlobalIndex(
                l, variable_id, component);
            double const value = input_vector[input_index];
 
            output_vector.getComponent(node_id, component) =
                map_function(value);
        }
    }
}

template <typename Functor>
void transformVariableFromGlobalVector(
    GlobalVector const& input_vector_on_bulk_mesh, int const variable_id,
    NumLib::LocalToGlobalIndexMap const& local_to_global_index_map,
    MeshLib::PropertyVector<double>& output_vector_on_submesh,
    std::span<std::size_t const> const&
        map_submesh_node_id_to_bulk_mesh_node_id,
    Functor map_function)
{
    MathLib::LinAlg::setLocalAccessibleVector(input_vector_on_bulk_mesh);
 
    // We fill the output with zeroes, because filling with NaN
    // "breaks" visualization, e.g. of heat or mass flow rate
    // with Taylor-Hood elements. I.e., all lower order
    // properties would have NaN on the higher order nodes.
    ranges::fill(output_vector_on_submesh, 0);
 
    int const n_components =
        local_to_global_index_map.getNumberOfVariableComponents(variable_id);
    for (int component = 0; component < n_components; ++component)
    {
        auto const& mesh_subset =
            local_to_global_index_map.getMeshSubset(variable_id, component);
        auto const mesh_id = mesh_subset.getMeshID();
 
        for (std::size_t submesh_node_id = 0;
             submesh_node_id < map_submesh_node_id_to_bulk_mesh_node_id.size();
             ++submesh_node_id)
        {
            std::size_t const bulk_node_id =
                map_submesh_node_id_to_bulk_mesh_node_id[submesh_node_id];
            MeshLib::Location const l(mesh_id, MeshLib::MeshItemType::Node,
                                      bulk_node_id);
            auto const input_index = local_to_global_index_map.getGlobalIndex(
                l, variable_id, component);
 
            // On higher-order meshes some d.o.f.s might not be defined on all
            // nodes. We silently ignore the d.o.f.s we cannot find.
            [[unlikely]] if (input_index == NumLib::MeshComponentMap::nop)
            {
                continue;
            }
 
            double const value = input_vector_on_bulk_mesh[input_index];
 
            output_vector_on_submesh.getComponent(submesh_node_id, component) =
                map_function(value);
        }
    }
}
 
Eigen::VectorXd getLocalX(
    std::size_t const mesh_item_id,
    std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
    std::vector<GlobalVector*> const& x);
}  // namespace NumLib