Source code for mmgp.FE_utils

# -*- coding: utf-8 -*-
#
# This file is subject to the terms and conditions defined in
# file 'LICENSE.txt', which is part of this source code package.
#
#
import Muscat.Containers.Filters.FilterObjects as FilterObjects
import Muscat.Containers.MeshInspectionTools as MIT
import networkx
import numpy as np
from Muscat.Containers import MeshFieldOperations as MFO
from Muscat.Containers.Mesh import Mesh
from Muscat.FE import FETools as FT
from Muscat.FE.Fields import FEField as FF
import time
from typing import Tuple, Callable




def compute_FE_projection_operators(origin_mesh: Mesh, target_mesh: Mesh,
                                    verbose: bool = False) -> Tuple[np.ndarray, np.ndarray]:
    """Compute Finite Element Projection Operators.

    Args:
        origin_mesh (Mesh): The original mesh data.
        target_mesh (Mesh): The target mesh data.
        verbose (bool, optional): Whether to display verbose output. Defaults to False.

    Returns:
        tuple(np.ndarray, np.ndarray): A tuple containing two projection operators (projOperator, invProjOperator).
    """
    start = time.time()
    if verbose is True: # pragma: no cover
        print("Computing direct FE interpolation operator:")
    space, numberings, _, _ = FT.PrepareFEComputation(
        origin_mesh, numberOfComponents=1)
    input_FE_field = FF.FEField(
        name="dummy",
        mesh=origin_mesh,
        space=space,
        numbering=numberings[0])
    proj_operator = MFO.GetFieldTransferOp(
        input_FE_field, target_mesh.nodes, method="Interp/Clamp", verbose=verbose)[0]

    if verbose is True: # pragma: no cover
        print("Computing inverse FE interpolation operator:")
    space, numberings, _, _ = FT.PrepareFEComputation(
        target_mesh, numberOfComponents=1)
    input_FE_field = FF.FEField(
        name="dummy",
        mesh=target_mesh,
        space=space,
        numbering=numberings[0])
    inv_proj_operator = MFO.GetFieldTransferOp(
        input_FE_field, origin_mesh.nodes, method="Interp/Clamp", verbose=verbose)[0]

    if verbose is True: # pragma: no cover
        print(
            "Duration computation of FE interpolation operators: {:#.6g} s".format(
                time.time() -
                start))

    return proj_operator, inv_proj_operator





def compute_node_to_node_graph(
        in_mesh: Mesh, dimensionality: int = None, dist_func: Callable=None) -> networkx.Graph:
    '''Creates a networkx graph from the node connectivity on a Mesh through edges

    Parameters
    ----------
    in_mesh : Mesh
        input mesh
    dimensionality : int
        dimension of the elements considered to initalize the graph
    dist_func : func
        function applied to the lengh of the edges of the mesh, and attached of the
        corresponding edge of the graph of the mesh

    Returns
    -------
    networkx.Graph
        Element to element graph
    '''
    if dimensionality is None: # pragma: no cover
        dimensionality = in_mesh.GetDimensionality()

    if dist_func is None:
        def dist_func(x):
            return x

    el_filter = FilterObjects.ElementFilter(dimensionality=dimensionality)
    mesh = MIT.ExtractElementsByElementFilter(in_mesh, el_filter)

    node_connectivity, _ = MIT.ComputeNodeToNodeConnectivity(mesh)

    G = initialize_graph_points_from_mesh_points(in_mesh)
    edges = []
    for i in range(node_connectivity.shape[0]):
        for j in node_connectivity[i][node_connectivity[i] > i]:
            length = np.linalg.norm(in_mesh.nodes[i] - in_mesh.nodes[j])
            edges.append((i, j, dist_func(length)))
    G.add_weighted_edges_from(edges)

    return G





def initialize_graph_points_from_mesh_points(
        in_mesh: Mesh) -> networkx.Graph:
    '''Initializes a networkx graph with nodes consistant with the number of nodes of a Mesh.
    This enables further edge addition compatible with the connectivity of the elements of the Mesh.

    Parameters
    ----------
    in_mesh : Mesh
        input mesh

    Returns
    -------
    networkx.Graph
        initialized graph
    '''
    G = networkx.Graph()
    G.add_nodes_from(np.arange(in_mesh.GetNumberOfNodes()))
    return G