text: /**
   * \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 <limits>
  #include <memory>
  #include <vector>
  
  #include "LocalAssemblerInterface.h"
  #include "MaterialLib/PhysicalConstant.h"
  #include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
  #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
  #include "NumLib/Extrapolation/ExtrapolatableElement.h"
  #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
  #include "NumLib/Fem/InitShapeMatrices.h"
  #include "NumLib/Fem/ShapeMatrixPolicy.h"
  #include "ParameterLib/Parameter.h"
  #include "ProcessLib/Deformation/BMatrixPolicy.h"
  #include "ProcessLib/Deformation/GMatrixPolicy.h"
  #include "ProcessLib/Deformation/LinearBMatrix.h"
  #include "ProcessLib/LocalAssemblerInterface.h"
  #include "ProcessLib/LocalAssemblerTraits.h"
  #include "ProcessLib/Utils/SetOrGetIntegrationPointData.h"
  #include "SmallDeformationProcessData.h"
  
  namespace ProcessLib
  {
  namespace SmallDeformation
  {
  namespace MPL = MaterialPropertyLib;
  
  template <typename BMatricesType, typename ShapeMatricesType,
            int DisplacementDim>
  struct IntegrationPointData final
  {
      explicit IntegrationPointData(
          MaterialLib::Solids::MechanicsBase<DisplacementDim> const&
              solid_material)
          : solid_material(solid_material),
            material_state_variables(
                solid_material.createMaterialStateVariables())
      {
      }
  
      typename BMatricesType::KelvinVectorType sigma, sigma_prev;
      typename BMatricesType::KelvinVectorType eps, eps_prev;
      double free_energy_density = 0;
  
      MaterialLib::Solids::MechanicsBase<DisplacementDim> const& solid_material;
      std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
          DisplacementDim>::MaterialStateVariables>
          material_state_variables;
  
      double integration_weight;
      typename ShapeMatricesType::NodalRowVectorType N;
      typename ShapeMatricesType::GlobalDimNodalMatrixType dNdx;
  
      void pushBackState()
      {
          eps_prev = eps;
          sigma_prev = sigma;
          material_state_variables->pushBackState();
      }
  
      EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
  };
  
  /// 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 SmallDeformationLocalAssembler
      : public SmallDeformationLocalAssemblerInterface<DisplacementDim>
  {
  public:
      using ShapeMatricesType =
          ShapeMatrixPolicyType<ShapeFunction, DisplacementDim>;
      using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType;
      using NodalVectorType = typename ShapeMatricesType::NodalVectorType;
      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 GMatricesType = GMatrixPolicyType<ShapeFunction, DisplacementDim>;
      using GradientVectorType = typename GMatricesType::GradientVectorType;
      using GradientMatrixType = typename GMatricesType::GradientMatrixType;
      using IpData =
          IntegrationPointData<BMatricesType, ShapeMatricesType, DisplacementDim>;
  
      SmallDeformationLocalAssembler(SmallDeformationLocalAssembler const&) =
          delete;
      SmallDeformationLocalAssembler(SmallDeformationLocalAssembler&&) = delete;
  
      SmallDeformationLocalAssembler(
          MeshLib::Element const& e,
          std::size_t const /*local_matrix_size*/,
          bool const is_axially_symmetric,
          unsigned const integration_order,
          SmallDeformationProcessData<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());
  
          for (unsigned ip = 0; ip < n_integration_points; ip++)
          {
              _ip_data.emplace_back(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;
  
              static const int kelvin_vector_size =
                  MathLib::KelvinVector::kelvin_vector_dimensions(
                      DisplacementDim);
              // Initialize current time step values
              ip_data.sigma.setZero(kelvin_vector_size);
              ip_data.eps.setZero(kelvin_vector_size);
  
              // Previous time step values are not initialized and are set later.
              ip_data.sigma_prev.resize(kelvin_vector_size);
              ip_data.eps_prev.resize(kelvin_vector_size);
  
              _secondary_data.N[ip] = shape_matrices[ip].N;
          }
      }
  
      /// Returns number of read integration points.
      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")
          {
              if (_process_data.initial_stress != nullptr)
              {
                  OGS_FATAL(
                      "Setting initial conditions for stress from integration "
                      "point data and from a parameter '{:s}' is not possible "
                      "simultaneously.",
                      _process_data.initial_stress->name);
              }
              return setSigma(values);
          }
  
          return 0;
      }
  
      void initializeConcrete() override
      {
          unsigned const n_integration_points =
              _integration_method.getNumberOfPoints();
          for (unsigned ip = 0; ip < n_integration_points; ip++)
          {
              auto& ip_data = _ip_data[ip];
  
              /// Set initial stress from parameter.
              if (_process_data.initial_stress != nullptr)
              {
                  ParameterLib::SpatialPosition const x_position{
                      std::nullopt, _element.getID(), ip,
                      MathLib::Point3d(
                          NumLib::interpolateCoordinates<ShapeFunction,
                                                         ShapeMatricesType>(
                              _element, ip_data.N))};
  
                  ip_data.sigma =
                      MathLib::KelvinVector::symmetricTensorToKelvinVector<
                          DisplacementDim>((*_process_data.initial_stress)(
                          std::numeric_limits<
                              double>::quiet_NaN() /* time independent */,
                          x_position));
              }
  
              ip_data.pushBackState();
          }
      }
  
      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(
              "SmallDeformationLocalAssembler: assembly without jacobian is not "
              "implemented.");
      }
  
      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();
  
          MPL::VariableArray variables_prev;
          MPL::VariableArray variables;
          ParameterLib::SpatialPosition x_position;
          x_position.setElementID(_element.getID());
  
          auto const& b = _process_data.specific_body_force;
  
          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;
  
              typename ShapeMatricesType::template MatrixType<DisplacementDim,
                                                              displacement_size>
                  N_u_op = ShapeMatricesType::template MatrixType<
                      DisplacementDim,
                      displacement_size>::Zero(DisplacementDim,
                                               displacement_size);
              for (int i = 0; i < DisplacementDim; ++i)
              {
                  N_u_op
                      .template block<1, displacement_size / DisplacementDim>(
                          i, i * displacement_size / DisplacementDim)
                      .noalias() = N;
              }
  
              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& state = _ip_data[ip].material_state_variables;
  
              eps.noalias() =
                  B *
                  Eigen::Map<typename BMatricesType::NodalForceVectorType const>(
                      local_x.data(), ShapeFunction::NPOINTS * DisplacementDim);
  
              variables_prev[static_cast<int>(MPL::Variable::stress)]
                  .emplace<
                      MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
                      sigma_prev);
              variables_prev[static_cast<int>(MPL::Variable::mechanical_strain)]
                  .emplace<
                      MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
                      eps_prev);
  
              double const T_ref =
                  _process_data.reference_temperature
                      ? (*_process_data.reference_temperature)(t, x_position)[0]
                      : std::numeric_limits<double>::quiet_NaN();
  
              variables_prev[static_cast<int>(MPL::Variable::temperature)]
                  .emplace<double>(T_ref);
              variables[static_cast<int>(MPL::Variable::mechanical_strain)]
                  .emplace<
                      MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
                      eps);
              variables[static_cast<int>(MPL::Variable::temperature)]
                  .emplace<double>(T_ref);
  
              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.");
              }
  
              MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> C;
              std::tie(sigma, state, C) = std::move(*solution);
  
              auto const rho = _process_data.solid_density(t, x_position)[0];
              local_b.noalias() -=
                  (B.transpose() * sigma - N_u_op.transpose() * rho * b) * w;
              local_Jac.noalias() += B.transpose() * C * B * w;
          }
      }
  
      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();
          }
  
          ParameterLib::SpatialPosition x_position;
          x_position.setElementID(_element.getID());
  
          for (unsigned ip = 0; ip < n_integration_points; ip++)
          {
              x_position.setIntegrationPoint(ip);
  
              auto& ip_data = _ip_data[ip];
  
              // Update free energy density needed for material forces.
              ip_data.free_energy_density =
                  ip_data.solid_material.computeFreeEnergyDensity(
                      t, x_position, dt, ip_data.eps, ip_data.sigma,
                      *ip_data.material_state_variables);
          }
      }
  
      std::vector<double> const& getMaterialForces(
          std::vector<double> const& local_x,
          std::vector<double>& nodal_values) override
      {
          return ProcessLib::SmallDeformation::getMaterialForces<
              DisplacementDim, ShapeFunction, ShapeMatricesType,
              typename BMatricesType::NodalForceVectorType,
              NodalDisplacementVectorType, GradientVectorType,
              GradientMatrixType>(local_x, nodal_values, _integration_method,
                                  _ip_data, _element, _is_axially_symmetric);
      }
  
      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& 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::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);
      }
  
      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);
      }
  
      unsigned getNumberOfIntegrationPoints() const override
      {
          return _integration_method.getNumberOfPoints();
      }
  
      typename MaterialLib::Solids::MechanicsBase<
          DisplacementDim>::MaterialStateVariables const&
      getMaterialStateVariablesAt(unsigned integration_point) const override
      {
          return *_ip_data[integration_point].material_state_variables;
      }
  
      void computeSecondaryVariableConcrete(
          double const /*t*/, double const /*dt*/, Eigen::VectorXd const& /*x*/,
          Eigen::VectorXd const& /*x_dot*/) override
      {
          int const elem_id = _element.getID();
          ParameterLib::SpatialPosition x_position;
          x_position.setElementID(elem_id);
          unsigned const n_integration_points =
              _integration_method.getNumberOfPoints();
  
          auto sigma_sum = MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
              Eigen::Matrix<double, 3, 3>::Zero());
  
          for (unsigned ip = 0; ip < n_integration_points; ip++)
          {
              x_position.setIntegrationPoint(ip);
              auto const& sigma = _ip_data[ip].sigma;
              sigma_sum += sigma;
          }
  
          Eigen::Matrix<double, 3, 3, 0, 3, 3> const sigma_avg =
              MathLib::KelvinVector::kelvinVectorToTensor(sigma_sum) /
              n_integration_points;
  
          Eigen::SelfAdjointEigenSolver<Eigen::Matrix<double, 3, 3>> e_s(
              sigma_avg);
  
          Eigen::Map<Eigen::Vector3d>(
              &(*_process_data.principal_stress_values)[elem_id * 3], 3) =
              e_s.eigenvalues();
  
          auto eigen_vectors = e_s.eigenvectors();
  
          for (auto i = 0; i < 3; i++)
          {
              Eigen::Map<Eigen::Vector3d>(
                  &(*_process_data.principal_stress_vector[i])[elem_id * 3], 3) =
                  eigen_vectors.col(i);
          }
      }
  
  private:
      SmallDeformationProcessData<DisplacementDim>& _process_data;
  
      std::vector<IpData, Eigen::aligned_allocator<IpData>> _ip_data;
  
      IntegrationMethod _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 SmallDeformation
  }  // namespace ProcessLib
  