<commit_before>from __future__ import absolute_import
import numpy as nm

import sfepy.linalg as la
from sfepy.discrete.integrals import Integral
from six.moves import range

def prepare_remap(indices, n_full):
    """
    Prepare vector for remapping range `[0, n_full]` to its subset given
    by `indices`.
    """
    remap = nm.empty((n_full,), dtype=nm.int32)
    remap.fill(-1)
    remap[indices] = nm.arange(indices.shape[0], dtype=nm.int32)

    return remap

def invert_remap(remap):
    """
    Return the inverse of `remap`, i.e. a mapping from a sub-range
    indices to a full range, see :func:`prepare_remap()`.
    """
    if remap is not None:
        inverse = nm.where(remap >= 0)[0].astype(nm.int32)

    else:
        inverse = None

    return inverse

def prepare_translate(old_indices, new_indices):
    """
    Prepare vector for translating `old_indices` to `new_indices`.

    Returns
    -------
    translate : array
        The translation vector. Then `new_ar = translate[old_ar]`.
    """
    old_indices = nm.asarray(old_indices)
    new_indices = nm.asarray(new_indices)

    translate = nm.zeros(old_indices.max() + 1, dtype=new_indices.dtype)
    translate[old_indices] = new_indices

    return translate

def compute_nodal_normals(nodes, region, field, return_imap=False):
    """
    Nodal normals are computed by simple averaging of element normals of
    elements every node is contained in.
    """
    dim = region.dim

    field.domain.create_surface_group(region)
    field.setup_surface_data(region)

    # Custom integral with quadrature points very close to facet vertices.
    coors = field.gel.surface_facet.coors
    centre = coors.sum(axis=0) / coors.shape[0]
    qp_coors = coors + 1e-8 * (centre - coors)
    # Unit normals -> weights = ones.
    qp_weights = nm.ones(qp_coors.shape[0], dtype=nm.float64)

    integral = Integral('aux', coors=qp_coors, weights=qp_weights)

    normals = nm.zeros((nodes.shape[0], dim), dtype=nm.float64)
    mask = nm.zeros((nodes.max() + 1,), dtype=nm.int32)
    imap = nm.empty_like(mask)
    imap.fill(nodes.shape[0]) # out-of-range index for normals.
    imap[nodes] = nm.arange(nodes.shape[0], dtype=nm.int32)

    cmap, _ = field.get_mapping(region, integral, 'surface')
    e_normals = cmap.normal[..., 0]

    sd = field.domain.surface_groups[region.name]
    econn = sd.get_connectivity()
    mask[econn] += 1

    # normals[imap[econn]] += e_normals
    im = imap[econn]
    for ii, en in enumerate(e_normals):
        normals[im[ii]] += en

    # All nodes must have a normal.
    if not nm.all(mask[nodes] > 0):
        raise ValueError('region %s has not complete faces!' % region.name)

    norm = la.norm_l2_along_axis(normals)[:, nm.newaxis]
    if (norm < 1e-15).any():
        raise ValueError('zero nodal normal! (a node in volume?)')
    normals /= norm

    if return_imap:
        return normals, imap

    else:<commit_msg>fix compute_nodal_normals() for higher order nodes<commit_after>from __future__ import absolute_import
import numpy as nm

import sfepy.linalg as la
from sfepy.discrete.integrals import Integral
from sfepy.discrete.fem.poly_spaces import PolySpace
from six.moves import range

def prepare_remap(indices, n_full):
    """
    Prepare vector for remapping range `[0, n_full]` to its subset given
    by `indices`.
    """
    remap = nm.empty((n_full,), dtype=nm.int32)
    remap.fill(-1)
    remap[indices] = nm.arange(indices.shape[0], dtype=nm.int32)

    return remap

def invert_remap(remap):
    """
    Return the inverse of `remap`, i.e. a mapping from a sub-range
    indices to a full range, see :func:`prepare_remap()`.
    """
    if remap is not None:
        inverse = nm.where(remap >= 0)[0].astype(nm.int32)

    else:
        inverse = None

    return inverse

def prepare_translate(old_indices, new_indices):
    """
    Prepare vector for translating `old_indices` to `new_indices`.

    Returns
    -------
    translate : array
        The translation vector. Then `new_ar = translate[old_ar]`.
    """
    old_indices = nm.asarray(old_indices)
    new_indices = nm.asarray(new_indices)

    translate = nm.zeros(old_indices.max() + 1, dtype=new_indices.dtype)
    translate[old_indices] = new_indices

    return translate

def compute_nodal_normals(nodes, region, field, return_imap=False):
    """
    Nodal normals are computed by simple averaging of element normals of
    elements every node is contained in.
    """
    dim = region.dim

    field.domain.create_surface_group(region)
    field.setup_surface_data(region)

    # Custom integral with quadrature points in nodes.
    ps = PolySpace.any_from_args('', field.gel.surface_facet,
                                 field.approx_order)
    qp_coors = ps.node_coors
    # Unit normals -> weights = ones.
    qp_weights = nm.ones(qp_coors.shape[0], dtype=nm.float64)

    integral = Integral('aux', coors=qp_coors, weights=qp_weights)

    normals = nm.zeros((nodes.shape[0], dim), dtype=nm.float64)
    mask = nm.zeros((nodes.max() + 1,), dtype=nm.int32)
    imap = nm.empty_like(mask)
    imap.fill(nodes.shape[0]) # out-of-range index for normals.
    imap[nodes] = nm.arange(nodes.shape[0], dtype=nm.int32)

    cmap, _ = field.get_mapping(region, integral, 'surface')
    e_normals = cmap.normal[..., 0]

    sd = field.surface_data[region.name]
    econn = sd.get_connectivity()
    mask[econn] += 1

    # normals[imap[econn]] += e_normals
    im = imap[econn]
    for ii, en in enumerate(e_normals):
        normals[im[ii]] += en

    # All nodes must have a normal.
    if not nm.all(mask[nodes] > 0):
        raise ValueError('region %s has not complete faces!' % region.name)

    norm = la.norm_l2_along_axis(normals)[:, nm.newaxis]
    if (norm < 1e-15).any():
        raise ValueError('zero nodal normal! (a node in volume?)')
    normals /= norm

    if return_imap:
        return normals, imap

    else:
{'fraction_non_alphanumeric': 0.12112533692722371, 'fraction_numerical': 0.011455525606469003, 'mean_word_length': 3.122916666666667, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2199955', 'n_tokens_mistral': 2098, 'n_tokens_neox': 1993, 'n_words': 497}