swh:1:snp:f521c49ab17ef7db6ec70b2430e1ed203f50383f
Tip revision: 5dcf00c2a951f73adb11e19b3ec72f01452f2ee9 authored by Lars Bilke on 06 July 2021, 08:05:31 UTC
Moved LIS dependency from Conan to ExternalProject.
Moved LIS dependency from Conan to ExternalProject.
Tip revision: 5dcf00c
Interpolation.h
/**
* \file
* \copyright
* Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*
*/
#pragma once
#include <array>
#include <cassert>
#include "NumLib/Fem/CoordinatesMapping/NaturalNodeCoordinates.h"
#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
#include "NumLib/Fem/InitShapeMatrices.h"
namespace NumLib
{
namespace detail
{
//! \see ::NumLib::shapeFunctionInterpolate()
template <unsigned DOFOffset, typename NodalValues, typename ShapeMatrix>
void shapeFunctionInterpolate(const NodalValues& /*nodal_values*/,
const ShapeMatrix& /*shape_matrix_N*/)
{
}
//! \see ::NumLib::shapeFunctionInterpolate()
template <unsigned DOFOffset, typename NodalValues, typename ShapeMatrix,
typename... ScalarTypes>
void shapeFunctionInterpolate(const NodalValues& nodal_values,
const ShapeMatrix& shape_matrix_N,
double& interpolated_value,
ScalarTypes&... interpolated_values)
{
auto const num_nodes = shape_matrix_N.size();
double iv = 0.0;
for (auto n = decltype(num_nodes){0}; n < num_nodes; ++n)
{
iv += nodal_values[DOFOffset * num_nodes + n] * shape_matrix_N[n];
}
interpolated_value = iv;
shapeFunctionInterpolate<DOFOffset + 1>(nodal_values, shape_matrix_N,
interpolated_values...);
}
} // namespace detail
/*!
* Interpolates variables given at element nodes according to the given shape
* matrix.
*
* This function simply does the usual finite-element interpolation, i.e.
* multiplication of nodal values with the shape function.
*
* \param nodal_values vector of nodal values, ordered by component
* \param shape_matrix_N shape matrix of the point to which will be interpolated
* \param interpolated_value interpolated value of the first d.o.f. (output
* parameter) \param interpolated_values interpolated value of further d.o.f.
* (output parameter)
*
* \tparam NodalValues type of the container where nodal values are stored
* \tparam ShapeMatrix type of the shape matrix \f$N\f$.
* \tparam ScalarValues all of the types in this pack must currently be \c
* double.
*
* \note
* \c nodal_values have to be ordered by component and it is assumed that all
* passed d.o.f. are single-component and are interpolated using the same shape
* function.
*/
template <typename NodalValues, typename ShapeMatrix, typename... ScalarTypes>
void shapeFunctionInterpolate(const NodalValues& nodal_values,
const ShapeMatrix& shape_matrix_N,
double& interpolated_value,
ScalarTypes&... interpolated_values)
{
auto const num_nodal_dof = sizeof...(interpolated_values) + 1;
auto const num_nodes = shape_matrix_N.size();
assert(num_nodes * num_nodal_dof ==
static_cast<std::size_t>(nodal_values.size()));
(void)num_nodal_dof;
(void)num_nodes; // no warnings when not in debug build
detail::shapeFunctionInterpolate<0>(nodal_values, shape_matrix_N,
interpolated_value,
interpolated_values...);
}
/// Interpolates scalar \c node_values given in lower order element nodes (e.g.
/// the base nodes) to higher order element's nodes (e.g. quadratic nodes) and
/// writes the result into the global property vector.
///
/// The base nodes' values are copied. For each higher order node the shape
/// matrices are evaluated for the lower order element (the base nodes), and
/// used for the the scalar quantity interpolation.
template <typename LowerOrderShapeFunction, typename HigherOrderMeshElementType,
int GlobalDim, typename EigenMatrixType>
void interpolateToHigherOrderNodes(
MeshLib::Element const& element, bool const is_axially_symmetric,
Eigen::MatrixBase<EigenMatrixType> const& node_values,
MeshLib::PropertyVector<double>& interpolated_values_global_vector)
{
assert(dynamic_cast<HigherOrderMeshElementType const*>(&element));
assert(node_values.cols() == 1); // Scalar quantity only.
using SF = LowerOrderShapeFunction;
using ShapeMatricesType = ShapeMatrixPolicyType<SF, GlobalDim>;
int const number_base_nodes = element.getNumberOfBaseNodes();
int const number_all_nodes = element.getNumberOfNodes();
// Copy the values for linear nodes.
for (int n = 0; n < number_base_nodes; ++n)
{
std::size_t const global_index = element.getNodeIndex(n);
interpolated_values_global_vector[global_index] = node_values[n];
}
//
// Interpolate values for higher order nodes.
//
int const number_higher_order_nodes = number_all_nodes - number_base_nodes;
std::vector<MathLib::Point3d> higher_order_nodes;
higher_order_nodes.reserve(number_higher_order_nodes);
for (int n = 0; n < number_higher_order_nodes; ++n)
{
higher_order_nodes.emplace_back(
NaturalCoordinates<HigherOrderMeshElementType>::coordinates
[number_base_nodes + n]);
}
// Shape matrices evaluated at higher order nodes' coordinates.
auto const shape_matrices =
computeShapeMatrices<SF, ShapeMatricesType, GlobalDim,
ShapeMatrixType::N>(element, is_axially_symmetric,
higher_order_nodes);
for (int n = 0; n < number_higher_order_nodes; ++n)
{
std::size_t const global_index =
element.getNodeIndex(number_base_nodes + n);
interpolated_values_global_vector[global_index] =
shape_matrices[n].N * node_values;
}
}
} // namespace NumLib