/**
 * \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 <cassert>
#include <vector>

#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
#include "MathLib/LinAlg/LinAlg.h"
#include "MathLib/LinAlg/MatrixVectorTraits.h"
#include "NumLib/DOF/DOFTableUtil.h"
#include "NumLib/Fem/InitShapeMatrices.h"
#include "ProcessLib/Deformation/GMatrix.h"

namespace ProcessLib
{
namespace SmallDeformation
{
struct MaterialForcesInterface
{
    virtual std::vector<double> const& getMaterialForces(
        std::vector<double> const& local_x,
        std::vector<double>& nodal_values) = 0;

    virtual ~MaterialForcesInterface() = default;
};

template <int DisplacementDim, typename ShapeFunction,
          typename ShapeMatricesType, typename NodalForceVectorType,
          typename NodalDisplacementVectorType, typename GradientVectorType,
          typename GradientMatrixType, typename IPData,
          typename IntegrationMethod>
std::vector<double> const& getMaterialForces(
    std::vector<double> const& local_x, std::vector<double>& nodal_values,
    IntegrationMethod const& _integration_method, IPData const& _ip_data,
    MeshLib::Element const& element, bool const is_axially_symmetric)
{
    unsigned const n_integration_points =
        _integration_method.getNumberOfPoints();

    nodal_values.clear();
    auto local_b = MathLib::createZeroedVector<NodalDisplacementVectorType>(
        nodal_values, ShapeFunction::NPOINTS * DisplacementDim);

    for (unsigned ip = 0; ip < n_integration_points; ip++)
    {
        auto const& sigma = _ip_data[ip].sigma;
        auto const& N = _ip_data[ip].N;
        auto const& dNdx = _ip_data[ip].dNdx;

        auto const& psi = _ip_data[ip].free_energy_density;

        auto const x_coord =
            NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
                element, _ip_data[ip].N);

        // For the 2D case the 33-component is needed (and the four entries
        // of the non-symmetric matrix); In 3d there are nine entries.
        GradientMatrixType G(
            DisplacementDim * DisplacementDim + (DisplacementDim == 2 ? 1 : 0),
            ShapeFunction::NPOINTS * DisplacementDim);
        Deformation::computeGMatrix<DisplacementDim, ShapeFunction::NPOINTS>(
            dNdx, G, is_axially_symmetric, N, x_coord);

        GradientVectorType const grad_u =
            G * Eigen::Map<NodalForceVectorType const>(
                    local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);

        GradientVectorType eshelby_stress;
        eshelby_stress.setZero(DisplacementDim * DisplacementDim +
                               (DisplacementDim == 2 ? 1 : 0));

        if (DisplacementDim == 3)
        {
            eshelby_stress[0] = eshelby_stress[DisplacementDim + 1] =
                eshelby_stress[8] = psi;

            eshelby_stress[0] -= sigma[0] * grad_u[0] +
                                 sigma[3] / std::sqrt(2) * grad_u[3] +
                                 sigma[5] / std::sqrt(2) * grad_u[6];

            eshelby_stress[1] -= sigma[3] / std::sqrt(2) * grad_u[0] +
                                 sigma[1] * grad_u[3] +
                                 sigma[4] / std::sqrt(2) * grad_u[6];

            eshelby_stress[2] -= sigma[5] / std::sqrt(2) * grad_u[0] +
                                 sigma[4] / std::sqrt(2) * grad_u[3] +
                                 sigma[2] * grad_u[6];

            eshelby_stress[3] -= sigma[0] * grad_u[1] +
                                 sigma[3] / std::sqrt(2) * grad_u[4] +
                                 sigma[5] / std::sqrt(2) * grad_u[7];

            eshelby_stress[4] -= sigma[3] / std::sqrt(2) * grad_u[1] +
                                 sigma[1] * grad_u[4] +
                                 sigma[4] / std::sqrt(2) * grad_u[7];

            eshelby_stress[5] -= sigma[5] / std::sqrt(2) * grad_u[1] +
                                 sigma[4] / std::sqrt(2) * grad_u[4] +
                                 sigma[2] * grad_u[7];

            eshelby_stress[6] -= sigma[0] * grad_u[2] +
                                 sigma[3] / std::sqrt(2) * grad_u[5] +
                                 sigma[5] / std::sqrt(2) * grad_u[8];

            eshelby_stress[7] -= sigma[3] / std::sqrt(2) * grad_u[2] +
                                 sigma[1] * grad_u[5] +
                                 sigma[4] / std::sqrt(2) * grad_u[8];

            eshelby_stress[8] -= sigma[5] / std::sqrt(2) * grad_u[2] +
                                 sigma[4] / std::sqrt(2) * grad_u[5] +
                                 sigma[2] * grad_u[8];
        }
        else if (DisplacementDim == 2)
        {
            eshelby_stress[0] = eshelby_stress[DisplacementDim + 1] =
                eshelby_stress[4] = psi;

            eshelby_stress[0] -=
                sigma[0] * grad_u[0] + sigma[3] / std::sqrt(2) * grad_u[2];

            eshelby_stress[1] -=
                sigma[3] / std::sqrt(2) * grad_u[0] + sigma[1] * grad_u[2];

            eshelby_stress[2] -=
                sigma[0] * grad_u[1] + sigma[3] / std::sqrt(2) * grad_u[3];

            eshelby_stress[3] -=
                sigma[3] / std::sqrt(2) * grad_u[1] + sigma[1] * grad_u[3];

            // for axial-symmetric case the following is not zero in general
            eshelby_stress[4] -= sigma[2] * grad_u[4];
        }
        else
        {
            OGS_FATAL(
                "Material forces not implemented for displacement "
                "dimension "
                "other than 2 and 3.");
        }

        auto const& w = _ip_data[ip].integration_weight;
        local_b += G.transpose() * eshelby_stress * w;
    }

    return nodal_values;
}

template <typename LocalAssemblerInterface>
void writeMaterialForces(
    std::unique_ptr<GlobalVector>& material_forces,
    std::vector<std::unique_ptr<LocalAssemblerInterface>> const&
        local_assemblers,
    NumLib::LocalToGlobalIndexMap const& local_to_global_index_map,
    GlobalVector const& x)
{
    DBUG("Compute material forces for small deformation process.");

    // Prepare local storage and fetch values.
    MathLib::LinAlg::setLocalAccessibleVector(x);

    if (!material_forces)
    {
        material_forces =
            MathLib::MatrixVectorTraits<GlobalVector>::newInstance(x);
    }

    MathLib::LinAlg::set(*material_forces, 0);

    GlobalExecutor::executeDereferenced(
        [](const std::size_t mesh_item_id,
           LocalAssemblerInterface& local_assembler,
           const NumLib::LocalToGlobalIndexMap& dof_table,
           GlobalVector const& x,
           GlobalVector& node_values) {
            auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
            std::vector<double> local_data;
            auto const local_x = x.get(indices);

            local_assembler.getMaterialForces(local_x, local_data);

            assert(local_data.size() == indices.size());
            node_values.add(indices, local_data);
        },
        local_assemblers, local_to_global_index_map, x, *material_forces);
    MathLib::LinAlg::finalizeAssembly(*material_forces);
}

}  // namespace SmallDeformation
}  // namespace ProcessLib