/**
 * \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 <algorithm>
#include <limits>
#include <memory>
#include <vector>

#include "Damage.h"
#include "IntegrationPointData.h"
#include "LocalAssemblerInterface.h"
#include "MaterialLib/SolidModels/Ehlers.h"
#include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
#include "MeshLib/findElementsWithinRadius.h"
#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
#include "NumLib/Fem/InitShapeMatrices.h"
#include "NumLib/Fem/ShapeMatrixPolicy.h"
#include "NumLib/Function/Interpolation.h"
#include "ParameterLib/Parameter.h"
#include "ProcessLib/Deformation/BMatrixPolicy.h"
#include "ProcessLib/Deformation/Divergence.h"
#include "ProcessLib/Deformation/LinearBMatrix.h"
#include "ProcessLib/LocalAssemblerTraits.h"
#include "ProcessLib/Utils/SetOrGetIntegrationPointData.h"
#include "SmallDeformationNonlocalProcessData.h"

namespace ProcessLib
{
namespace SmallDeformationNonlocal
{
namespace MPL = MaterialPropertyLib;

/// Used for the extrapolation of the integration point values. It is ordered
/// (and stored) by integration points.
template <typename ShapeMatrixType>
struct SecondaryData
{
    std::vector<ShapeMatrixType, Eigen::aligned_allocator<ShapeMatrixType>> N;
};

template <typename ShapeFunction, typename IntegrationMethod,
          int DisplacementDim>
class SmallDeformationNonlocalLocalAssembler
    : public SmallDeformationNonlocalLocalAssemblerInterface<DisplacementDim>
{
public:
    using ShapeMatricesType =
        ShapeMatrixPolicyType<ShapeFunction, DisplacementDim>;
    using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType;
    using NodalVectorType = typename ShapeMatricesType::NodalVectorType;
    using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType;
    using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
    using BMatricesType = BMatrixPolicyType<ShapeFunction, DisplacementDim>;

    using BMatrixType = typename BMatricesType::BMatrixType;
    using StiffnessMatrixType = typename BMatricesType::StiffnessMatrixType;
    using NodalForceVectorType = typename BMatricesType::NodalForceVectorType;
    using NodalDisplacementVectorType =
        typename BMatricesType::NodalForceVectorType;
    using IpData =
        IntegrationPointData<BMatricesType, ShapeMatricesType, DisplacementDim>;

    SmallDeformationNonlocalLocalAssembler(
        SmallDeformationNonlocalLocalAssembler const&) = delete;
    SmallDeformationNonlocalLocalAssembler(
        SmallDeformationNonlocalLocalAssembler&&) = delete;

    SmallDeformationNonlocalLocalAssembler(
        MeshLib::Element const& e,
        std::size_t const /*local_matrix_size*/,
        bool const is_axially_symmetric,
        unsigned const integration_order,
        SmallDeformationNonlocalProcessData<DisplacementDim>& process_data)
        : _process_data(process_data),
          _integration_method(integration_order),
          _element(e),
          _is_axially_symmetric(is_axially_symmetric)
    {
        unsigned const n_integration_points =
            _integration_method.getNumberOfPoints();

        _ip_data.reserve(n_integration_points);
        _secondary_data.N.resize(n_integration_points);

        auto const shape_matrices =
            NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
                                      DisplacementDim>(e, is_axially_symmetric,
                                                       _integration_method);

        auto& solid_material =
            MaterialLib::Solids::selectSolidConstitutiveRelation(
                _process_data.solid_materials,
                _process_data.material_ids,
                e.getID());
        auto* ehlers_solid_material = dynamic_cast<
            MaterialLib::Solids::Ehlers::SolidEhlers<DisplacementDim>*>(
            &solid_material);
        if (ehlers_solid_material == nullptr)
        {
            OGS_FATAL(
                "The SmallDeformationNonlocalLocal process supports only "
                "Ehlers material at the moment. For other materials the "
                "interface must be extended first.");
        }

        for (unsigned ip = 0; ip < n_integration_points; ip++)
        {
            _ip_data.emplace_back(*ehlers_solid_material);
            auto& ip_data = _ip_data[ip];
            auto const& sm = shape_matrices[ip];
            _ip_data[ip].integration_weight =
                _integration_method.getWeightedPoint(ip).getWeight() *
                sm.integralMeasure * sm.detJ;

            ip_data.N = sm.N;
            ip_data.dNdx = sm.dNdx;

            // Initialize current time step values
            ip_data.sigma.setZero(MathLib::KelvinVector::KelvinVectorDimensions<
                                  DisplacementDim>::value);
            ip_data.eps.setZero(MathLib::KelvinVector::KelvinVectorDimensions<
                                DisplacementDim>::value);

            // Previous time step values are not initialized and are set later.
            ip_data.sigma_prev.resize(
                MathLib::KelvinVector::KelvinVectorDimensions<
                    DisplacementDim>::value);
            ip_data.eps_prev.resize(
                MathLib::KelvinVector::KelvinVectorDimensions<
                    DisplacementDim>::value);

            _secondary_data.N[ip] = shape_matrices[ip].N;

            ip_data.coordinates = getSingleIntegrationPointCoordinates(ip);
        }
    }

    std::size_t setIPDataInitialConditions(std::string const& name,
                                    double const* values,
                                    int const integration_order) override
    {
        if (integration_order !=
            static_cast<int>(_integration_method.getIntegrationOrder()))
        {
            OGS_FATAL(
                "Setting integration point initial conditions; The integration "
                "order of the local assembler for element {:d} is different "
                "from the integration order in the initial condition.",
                _element.getID());
        }

        if (name == "sigma_ip")
        {
            return setSigma(values);
        }

        if (name == "kappa_d_ip")
        {
            return ProcessLib::setIntegrationPointScalarData(values, _ip_data,
                                                             &IpData::kappa_d);
        }

        return 0;
    }

    void setIPDataInitialConditionsFromCellData(
        std::string const& name, std::vector<double> const& value) override
    {
        if (name == "kappa_d_ip")
        {
            if (value.size() != 1)
            {
                OGS_FATAL(
                    "CellData for kappa_d initial conditions has wrong number "
                    "of components. 1 expected, got {:d}.",
                    value.size());
            }
            setKappaD(value[0]);
        }
    }

    double alpha_0(double const distance2) const
    {
        double const internal_length2 = _process_data.internal_length_squared;
        return (distance2 > internal_length2)
                   ? 0
                   : (1 - distance2 / (internal_length2)) *
                         (1 - distance2 / (internal_length2));
    }

    void nonlocal(
        std::size_t const /*mesh_item_id*/,
        std::vector<
            std::unique_ptr<SmallDeformationNonlocalLocalAssemblerInterface<
                DisplacementDim>>> const& local_assemblers) override
    {
        auto const search_element_ids = MeshLib::findElementsWithinRadius(
            _element, _process_data.internal_length_squared);

        unsigned const n_integration_points =
            _integration_method.getNumberOfPoints();

        std::vector<double> distances;  // Cache for ip-ip distances.
        //
        // For every integration point in this element collect the neighbouring
        // integration points falling in given radius (internal length) and
        // compute the alpha_kl weight.
        //
        for (unsigned k = 0; k < n_integration_points; k++)
        {
            //
            // Collect the integration points.
            //

            auto const& xyz = _ip_data[k].coordinates;

            // For all neighbors of element
            for (auto const search_element_id : search_element_ids)
            {
                auto const& la = local_assemblers[search_element_id];
                la->getIntegrationPointCoordinates(xyz, distances);
                for (int ip = 0; ip < static_cast<int>(distances.size()); ++ip)
                {
                    if (distances[ip] >= _process_data.internal_length_squared)
                    {
                        continue;
                    }
                    // save into current ip_k
                    _ip_data[k].non_local_assemblers.push_back(
                        {la->getIPDataPtr(ip),
                         std::numeric_limits<double>::quiet_NaN(),
                         distances[ip]});
                }
            }
            if (_ip_data[k].non_local_assemblers.size() == 0)
            {
                OGS_FATAL("no neighbours found!");
            }

            double a_k_sum_m = 0;
            for (auto const& tuple : _ip_data[k].non_local_assemblers)
            {
                double const distance2_m = tuple.distance2;

                auto const& w_m = tuple.ip_l_pointer->integration_weight;

                a_k_sum_m += w_m * alpha_0(distance2_m);
            }

            //
            // Calculate alpha_kl =
            //       alpha_0(|x_k - x_l|) / int_{m \in ip} alpha_0(|x_k - x_m|)
            //
            for (auto& tuple : _ip_data[k].non_local_assemblers)
            {
                double const distance2_l = tuple.distance2;
                double const a_kl = alpha_0(distance2_l) / a_k_sum_m;

                // Store the a_kl already multiplied with the integration
                // weight of that l integration point.
                auto const w_l = tuple.ip_l_pointer->integration_weight;
                tuple.alpha_kl_times_w_l = a_kl * w_l;
            }
        }
    }

    Eigen::Vector3d getSingleIntegrationPointCoordinates(
        int integration_point) const
    {
        auto const& N = _secondary_data.N[integration_point];

        Eigen::Vector3d xyz = Eigen::Vector3d::Zero();  // Resulting coordinates
        auto* nodes = _element.getNodes();
        for (int i = 0; i < N.size(); ++i)
        {
            Eigen::Map<Eigen::Vector3d const> node_coordinates{
                nodes[i]->getCoords(), 3};
            xyz += node_coordinates * N[i];
        }
        return xyz;
    }

    /// For each of the current element's integration points the squared
    /// distance from the current integration point is computed and stored in
    /// the given distances cache.
    void getIntegrationPointCoordinates(
        Eigen::Vector3d const& coords,
        std::vector<double>& distances) const override
    {
        unsigned const n_integration_points =
            _integration_method.getNumberOfPoints();

        distances.resize(n_integration_points);

        for (unsigned ip = 0; ip < n_integration_points; ip++)
        {
            auto const& xyz = _ip_data[ip].coordinates;
            distances[ip] = (xyz - coords).squaredNorm();
        }
    }

    void assemble(double const /*t*/, double const /*dt*/,
                  std::vector<double> const& /*local_x*/,
                  std::vector<double> const& /*local_xdot*/,
                  std::vector<double>& /*local_M_data*/,
                  std::vector<double>& /*local_K_data*/,
                  std::vector<double>& /*local_b_data*/) override
    {
        OGS_FATAL(
            "SmallDeformationNonlocalLocalAssembler: assembly without jacobian "
            "is not "
            "implemented.");
    }

    void preAssemble(double const t, double const dt,
                     std::vector<double> const& local_x) override
    {
        auto const n_integration_points =
            _integration_method.getNumberOfPoints();

        MPL::VariableArray variables;
        MPL::VariableArray variables_prev;
        ParameterLib::SpatialPosition x_position;
        x_position.setElementID(_element.getID());

        for (unsigned ip = 0; ip < n_integration_points; ip++)
        {
            x_position.setIntegrationPoint(ip);

            auto const& N = _ip_data[ip].N;
            auto const& dNdx = _ip_data[ip].dNdx;

            auto const x_coord =
                NumLib::interpolateXCoordinate<ShapeFunction,
                                               ShapeMatricesType>(_element, N);
            auto const B = LinearBMatrix::computeBMatrix<
                DisplacementDim, ShapeFunction::NPOINTS,
                typename BMatricesType::BMatrixType>(dNdx, N, x_coord,
                                                     _is_axially_symmetric);
            auto const& eps_prev = _ip_data[ip].eps_prev;
            auto const& sigma_prev = _ip_data[ip].sigma_prev;

            auto& eps = _ip_data[ip].eps;
            auto& sigma = _ip_data[ip].sigma;
            auto& C = _ip_data[ip].C;
            auto& state = _ip_data[ip].material_state_variables;
            double const& damage_prev = _ip_data[ip].damage_prev;

            eps.noalias() =
                B *
                Eigen::Map<typename BMatricesType::NodalForceVectorType const>(
                    local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);

            // sigma is for plastic part only.
            std::unique_ptr<
                MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>>
                new_C;
            std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
                DisplacementDim>::MaterialStateVariables>
                new_state;

            // Compute sigma_eff from damage total stress sigma
            using KelvinVectorType = typename BMatricesType::KelvinVectorType;
            KelvinVectorType const sigma_eff_prev =
                sigma_prev /
                (1. - damage_prev);  // damage_prev is in [0,1) range. See
                                     // calculateDamage() function.

            variables_prev[static_cast<int>(MPL::Variable::stress)]
                .emplace<
                    MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
                    sigma_eff_prev);
            variables_prev[static_cast<int>(MPL::Variable::mechanical_strain)]
                .emplace<
                    MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
                    eps_prev);
            variables_prev[static_cast<int>(MPL::Variable::temperature)]
                .emplace<double>(_process_data.reference_temperature);
            variables[static_cast<int>(MPL::Variable::mechanical_strain)]
                .emplace<
                    MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
                    eps);
            variables[static_cast<int>(MPL::Variable::temperature)]
                .emplace<double>(_process_data.reference_temperature);

            auto&& solution = _ip_data[ip].solid_material.integrateStress(
                variables_prev, variables, t, x_position, dt, *state);

            if (!solution)
            {
                OGS_FATAL("Computation of local constitutive relation failed.");
            }

            std::tie(sigma, state, C) = std::move(*solution);

            /// Compute only the local kappa_d.
            {
                auto const& ehlers_material =
                    static_cast<MaterialLib::Solids::Ehlers::SolidEhlers<
                        DisplacementDim> const&>(_ip_data[ip].solid_material);
                auto const damage_properties =
                    ehlers_material.evaluatedDamageProperties(t, x_position);
                auto const material_properties =
                    ehlers_material.evaluatedMaterialProperties(t, x_position);

                // Ehlers material state variables
                auto& state_vars =
                    static_cast<MaterialLib::Solids::Ehlers::StateVariables<
                        DisplacementDim>&>(
                        *_ip_data[ip].material_state_variables);

                double const eps_p_eff_diff =
                    state_vars.eps_p.eff - state_vars.eps_p_prev.eff;

                _ip_data[ip].kappa_d = calculateDamageKappaD<DisplacementDim>(
                    eps_p_eff_diff, sigma, _ip_data[ip].kappa_d_prev,
                    damage_properties.h_d, material_properties);

                if (!_ip_data[ip].active_self)
                {
                    _ip_data[ip].active_self |= _ip_data[ip].kappa_d > 0;
                    if (_ip_data[ip].active_self)
                    {
                        for (auto const& tuple :
                             _ip_data[ip].non_local_assemblers)
                        {
                            // Activate the integration point.
                            tuple.ip_l_pointer->activated = true;
                        }
                    }
                }
            }
        }
    }

    void assembleWithJacobian(double const t, double const /*dt*/,
                              std::vector<double> const& local_x,
                              std::vector<double> const& /*local_xdot*/,
                              const double /*dxdot_dx*/, const double /*dx_dx*/,
                              std::vector<double>& /*local_M_data*/,
                              std::vector<double>& /*local_K_data*/,
                              std::vector<double>& local_b_data,
                              std::vector<double>& local_Jac_data) override
    {
        auto const local_matrix_size = local_x.size();

        auto local_Jac = MathLib::createZeroedMatrix<StiffnessMatrixType>(
            local_Jac_data, local_matrix_size, local_matrix_size);

        auto local_b = MathLib::createZeroedVector<NodalDisplacementVectorType>(
            local_b_data, local_matrix_size);

        unsigned const n_integration_points =
            _integration_method.getNumberOfPoints();

        ParameterLib::SpatialPosition x_position;
        x_position.setElementID(_element.getID());

        // Non-local integration.
        for (unsigned ip = 0; ip < n_integration_points; ip++)
        {
            x_position.setIntegrationPoint(ip);
            auto const& w = _ip_data[ip].integration_weight;

            auto const& N = _ip_data[ip].N;
            auto const& dNdx = _ip_data[ip].dNdx;

            auto const x_coord =
                NumLib::interpolateXCoordinate<ShapeFunction,
                                               ShapeMatricesType>(_element, N);
            auto const B = LinearBMatrix::computeBMatrix<
                DisplacementDim, ShapeFunction::NPOINTS,
                typename BMatricesType::BMatrixType>(dNdx, N, x_coord,
                                                     _is_axially_symmetric);

            auto& sigma = _ip_data[ip].sigma;
            auto& C = _ip_data[ip].C;
            double& damage = _ip_data[ip].damage;

            {
                double nonlocal_kappa_d = 0;

                if (_ip_data[ip].active_self || _ip_data[ip].activated)
                {
                    for (auto const& tuple : _ip_data[ip].non_local_assemblers)
                    {
                        // Get local variable for the integration point l.
                        double const kappa_d_l = tuple.ip_l_pointer->kappa_d;
                        double const a_kl_times_w_l = tuple.alpha_kl_times_w_l;
                        nonlocal_kappa_d += a_kl_times_w_l * kappa_d_l;
                    }
                }

                auto const& ehlers_material =
                    static_cast<MaterialLib::Solids::Ehlers::SolidEhlers<
                        DisplacementDim> const&>(_ip_data[ip].solid_material);

                //
                // Overnonlocal formulation
                //
                // See (Di Luzio & Bazant 2005, IJSS) for details.
                // The implementation would go here and would be for a given
                // gamma_nonlocal:
                //
                // Update nonlocal damage with local damage (scaled with 1 -
                // \gamma_{nonlocal}) for the current integration point and the
                // nonlocal integral part.
                // nonlocal_kappa_d = (1. - gamma_nonlocal) * kappa_d +
                //                    gamma_nonlocal * nonlocal_kappa_d;

                nonlocal_kappa_d = std::max(0., nonlocal_kappa_d);

                // Update damage based on nonlocal kappa_d
                {
                    auto const damage_properties =
                        ehlers_material.evaluatedDamageProperties(t,
                                                                  x_position);
                    damage = calculateDamage(nonlocal_kappa_d,
                                             damage_properties.alpha_d,
                                             damage_properties.beta_d);
                    damage = std::max(0., damage);
                }
                sigma = sigma * (1. - damage);
            }

            local_b.noalias() -= B.transpose() * sigma * w;
            local_Jac.noalias() += B.transpose() * C * (1. - damage) * B * w;
        }
    }

    void initializeConcrete() override
    {
        unsigned const n_integration_points =
            _integration_method.getNumberOfPoints();

        for (unsigned ip = 0; ip < n_integration_points; ip++)
        {
            _ip_data[ip].pushBackState();
        }
    }

    void postTimestepConcrete(Eigen::VectorXd const& /*local_x*/,
                              double const /*t*/, double const /*dt*/) override
    {
        unsigned const n_integration_points =
            _integration_method.getNumberOfPoints();

        for (unsigned ip = 0; ip < n_integration_points; ip++)
        {
            _ip_data[ip].pushBackState();
        }
    }

    void computeCrackIntegral(std::size_t mesh_item_id,
                              NumLib::LocalToGlobalIndexMap const& dof_table,
                              GlobalVector const& x,
                              double& crack_volume) override
    {
        auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
        auto local_x = x.get(indices);

        auto u = Eigen::Map<typename BMatricesType::NodalForceVectorType const>(
            local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);

        int const n_integration_points =
            _integration_method.getNumberOfPoints();

        for (int ip = 0; ip < n_integration_points; ip++)
        {
            auto const& dNdx = _ip_data[ip].dNdx;
            auto const& d = _ip_data[ip].damage;
            auto const& w = _ip_data[ip].integration_weight;

            double const div_u =
                Deformation::divergence<DisplacementDim,
                                        ShapeFunction::NPOINTS>(u, dNdx);
            crack_volume += div_u * d * w;
        }
    }

    Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
        const unsigned integration_point) const override
    {
        auto const& N = _secondary_data.N[integration_point];

        // assumes N is stored contiguously in memory
        return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
    }

    std::vector<double> const& getNodalValues(
        std::vector<double>& nodal_values) const override
    {
        nodal_values.clear();
        auto local_b = MathLib::createZeroedVector<NodalDisplacementVectorType>(
            nodal_values, ShapeFunction::NPOINTS * DisplacementDim);

        unsigned const n_integration_points =
            _integration_method.getNumberOfPoints();

        for (unsigned ip = 0; ip < n_integration_points; ip++)
        {
            auto const& w = _ip_data[ip].integration_weight;

            auto const& N = _ip_data[ip].N;
            auto const& dNdx = _ip_data[ip].dNdx;

            auto const x_coord =
                NumLib::interpolateXCoordinate<ShapeFunction,
                                               ShapeMatricesType>(_element, N);
            auto const B = LinearBMatrix::computeBMatrix<
                DisplacementDim, ShapeFunction::NPOINTS,
                typename BMatricesType::BMatrixType>(dNdx, N, x_coord,
                                                     _is_axially_symmetric);
            auto& sigma = _ip_data[ip].sigma;

            local_b.noalias() += B.transpose() * sigma * w;
        }

        return nodal_values;
    }

    std::vector<double> const& getIntPtFreeEnergyDensity(
        const double /*t*/,
        std::vector<GlobalVector*> const& /*x*/,
        std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
        std::vector<double>& cache) const override
    {
        cache.clear();
        cache.reserve(_ip_data.size());

        transform(
            cbegin(_ip_data), cend(_ip_data), back_inserter(cache),
            [](auto const& ip_data) { return ip_data.free_energy_density; });

        return cache;
    }

    std::vector<double> const& getIntPtEpsPV(
        const double /*t*/,
        std::vector<GlobalVector*> const& /*x*/,
        std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
        std::vector<double>& cache) const override
    {
        cache.clear();
        cache.reserve(_ip_data.size());

        transform(cbegin(_ip_data), cend(_ip_data), back_inserter(cache),
                  [](auto const& ip_data) { return *ip_data.eps_p_V; });

        return cache;
    }

    std::vector<double> const& getIntPtEpsPDXX(
        const double /*t*/,
        std::vector<GlobalVector*> const& /*x*/,
        std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
        std::vector<double>& cache) const override
    {
        cache.clear();
        cache.reserve(_ip_data.size());

        transform(cbegin(_ip_data), cend(_ip_data), back_inserter(cache),
                  [](auto const& ip_data) { return *ip_data.eps_p_D_xx; });

        return cache;
    }

    std::vector<double> const& getIntPtSigma(
        const double /*t*/,
        std::vector<GlobalVector*> const& /*x*/,
        std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
        std::vector<double>& cache) const override
    {
        return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
            _ip_data, &IpData::sigma, cache);
    }

    std::vector<double> const& getIntPtEpsilon(
        const double /*t*/,
        std::vector<GlobalVector*> const& /*x*/,
        std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
        std::vector<double>& cache) const override
    {
        return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
            _ip_data, &IpData::eps, cache);
    }

    std::size_t setSigma(double const* values)
    {
        return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
            values, _ip_data, &IpData::sigma);
    }

    // TODO (naumov) This method is same as getIntPtSigma but for arguments and
    // the ordering of the cache_mat.
    // There should be only one.
    std::vector<double> getSigma() const override
    {
        return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
            _ip_data, &IpData::sigma);
    }

    void setKappaD(double value)
    {
        for (auto& ip_data : _ip_data)
        {
            ip_data.kappa_d = value;
        }
    }
    std::vector<double> getKappaD() const override
    {
        unsigned const n_integration_points =
            _integration_method.getNumberOfPoints();

        std::vector<double> result_values;
        result_values.resize(n_integration_points);

        for (unsigned ip = 0; ip < n_integration_points; ++ip)
        {
            result_values[ip] = _ip_data[ip].kappa_d;
        }

        return result_values;
    }

    std::vector<double> const& getIntPtDamage(
        const double /*t*/,
        std::vector<GlobalVector*> const& /*x*/,
        std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
        std::vector<double>& cache) const override
    {
        return ProcessLib::getIntegrationPointScalarData(
            _ip_data, &IpData::damage, cache);
    }

    unsigned getNumberOfIntegrationPoints() const override
    {
        return _integration_method.getNumberOfPoints();
    }

    typename MaterialLib::Solids::MechanicsBase<
        DisplacementDim>::MaterialStateVariables const&
    getMaterialStateVariablesAt(int const integration_point) const override
    {
        return *_ip_data[integration_point].material_state_variables;
    }

private:
    std::vector<double> const& getIntPtSigma(std::vector<double>& cache,
                                             std::size_t const component) const
    {
        cache.clear();
        cache.reserve(_ip_data.size());

        for (auto const& ip_data : _ip_data)
        {
            if (component < 3)
            {  // xx, yy, zz components
                cache.push_back(ip_data.sigma[component]);
            }
            else
            {  // mixed xy, yz, xz components
                cache.push_back(ip_data.sigma[component] / std::sqrt(2));
            }
        }

        return cache;
    }

    std::vector<double> const& getIntPtEpsilon(
        std::vector<double>& cache, std::size_t const component) const
    {
        cache.clear();
        cache.reserve(_ip_data.size());

        for (auto const& ip_data : _ip_data)
        {
            if (component < 3)  // xx, yy, zz components
                cache.push_back(ip_data.eps[component]);
            else  // mixed xy, yz, xz components
                cache.push_back(ip_data.eps[component] / std::sqrt(2));
        }

        return cache;
    }

    IntegrationPointDataNonlocalInterface*
    getIPDataPtr(int const ip) override
    {
        return &_ip_data[ip];
    }

private:
    SmallDeformationNonlocalProcessData<DisplacementDim>& _process_data;

    std::vector<IpData, Eigen::aligned_allocator<IpData>> _ip_data;

    IntegrationMethod const _integration_method;
    MeshLib::Element const& _element;
    SecondaryData<typename ShapeMatrices::ShapeType> _secondary_data;
    bool const _is_axially_symmetric;

    static const int displacement_size =
        ShapeFunction::NPOINTS * DisplacementDim;
};

}  // namespace SmallDeformationNonlocal
}  // namespace ProcessLib