text: /**
   * \copyright
   * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
   *            Distributed under a Modified BSD License.
   *              See accompanying file LICENSE.txt or
   *              http://www.opengeosys.org/project/license
   *
   */
  
  #include "HydroMechanicsProcess.h"
  
  #include <cassert>
  
  #include "MeshLib/Elements/Utils.h"
  #include "NumLib/DOF/ComputeSparsityPattern.h"
  #include "ProcessLib/HydroMechanics/CreateLocalAssemblers.h"
  #include "ProcessLib/Process.h"
  
  #include "HydroMechanicsFEM.h"
  #include "HydroMechanicsProcessData.h"
  
  namespace ProcessLib
  {
  namespace HydroMechanics
  {
  template <int DisplacementDim>
  HydroMechanicsProcess<DisplacementDim>::HydroMechanicsProcess(
      std::string name,
      MeshLib::Mesh& mesh,
      std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
      std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
      unsigned const integration_order,
      std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
          process_variables,
      HydroMechanicsProcessData<DisplacementDim>&& process_data,
      SecondaryVariableCollection&& secondary_variables,
      NumLib::NamedFunctionCaller&& named_function_caller,
      bool const use_monolithic_scheme)
      : Process(std::move(name), mesh, std::move(jacobian_assembler), parameters,
                integration_order, std::move(process_variables),
                std::move(secondary_variables), std::move(named_function_caller),
                use_monolithic_scheme),
        _process_data(std::move(process_data))
  {
      _nodal_forces = MeshLib::getOrCreateMeshProperty<double>(
          mesh, "NodalForces", MeshLib::MeshItemType::Node, DisplacementDim);
  
      _hydraulic_flow = MeshLib::getOrCreateMeshProperty<double>(
          mesh, "HydraulicFlow", MeshLib::MeshItemType::Node, 1);
  }
  
  template <int DisplacementDim>
  bool HydroMechanicsProcess<DisplacementDim>::isLinear() const
  {
      return false;
  }
  
  template <int DisplacementDim>
  MathLib::MatrixSpecifications
  HydroMechanicsProcess<DisplacementDim>::getMatrixSpecifications(
      const int process_id) const
  {
      // For the monolithic scheme or the M process (deformation) in the staggered
      // scheme.
      if (_use_monolithic_scheme || process_id == 1)
      {
          auto const& l = *_local_to_global_index_map;
          return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
                  &l.getGhostIndices(), &this->_sparsity_pattern};
      }
  
      // For staggered scheme and H process (pressure).
      auto const& l = *_local_to_global_index_map_with_base_nodes;
      return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
              &l.getGhostIndices(), &_sparsity_pattern_with_linear_element};
  }
  
  template <int DisplacementDim>
  void HydroMechanicsProcess<DisplacementDim>::constructDofTable()
  {
      // Create single component dof in every of the mesh's nodes.
      _mesh_subset_all_nodes =
          std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
      // Create single component dof in the mesh's base nodes.
      _base_nodes = MeshLib::getBaseNodes(_mesh.getElements());
      _mesh_subset_base_nodes =
          std::make_unique<MeshLib::MeshSubset>(_mesh, _base_nodes);
  
      // TODO move the two data members somewhere else.
      // for extrapolation of secondary variables of stress or strain
      std::vector<MeshLib::MeshSubset> all_mesh_subsets_single_component{
          *_mesh_subset_all_nodes};
      _local_to_global_index_map_single_component =
          std::make_unique<NumLib::LocalToGlobalIndexMap>(
              std::move(all_mesh_subsets_single_component),
              // by location order is needed for output
              NumLib::ComponentOrder::BY_LOCATION);
  
      if (_use_monolithic_scheme)
      {
          // For pressure, which is the first
          std::vector<MeshLib::MeshSubset> all_mesh_subsets{
              *_mesh_subset_base_nodes};
  
          // For displacement.
          const int monolithic_process_id = 0;
          std::generate_n(std::back_inserter(all_mesh_subsets),
                          getProcessVariables(monolithic_process_id)[1]
                              .get()
                              .getNumberOfComponents(),
                          [&]() { return *_mesh_subset_all_nodes; });
  
          std::vector<int> const vec_n_components{1, DisplacementDim};
          _local_to_global_index_map =
              std::make_unique<NumLib::LocalToGlobalIndexMap>(
                  std::move(all_mesh_subsets), vec_n_components,
                  NumLib::ComponentOrder::BY_LOCATION);
          assert(_local_to_global_index_map);
      }
      else
      {
          // For displacement equation.
          const int process_id = 1;
          std::vector<MeshLib::MeshSubset> all_mesh_subsets;
          std::generate_n(
              std::back_inserter(all_mesh_subsets),
              getProcessVariables(process_id)[0].get().getNumberOfComponents(),
              [&]() { return *_mesh_subset_all_nodes; });
  
          std::vector<int> const vec_n_components{DisplacementDim};
          _local_to_global_index_map =
              std::make_unique<NumLib::LocalToGlobalIndexMap>(
                  std::move(all_mesh_subsets), vec_n_components,
                  NumLib::ComponentOrder::BY_LOCATION);
  
          // For pressure equation.
          // Collect the mesh subsets with base nodes in a vector.
          std::vector<MeshLib::MeshSubset> all_mesh_subsets_base_nodes{
              *_mesh_subset_base_nodes};
          _local_to_global_index_map_with_base_nodes =
              std::make_unique<NumLib::LocalToGlobalIndexMap>(
                  std::move(all_mesh_subsets_base_nodes),
                  // by location order is needed for output
                  NumLib::ComponentOrder::BY_LOCATION);
  
          _sparsity_pattern_with_linear_element = NumLib::computeSparsityPattern(
              *_local_to_global_index_map_with_base_nodes, _mesh);
  
          assert(_local_to_global_index_map);
          assert(_local_to_global_index_map_with_base_nodes);
      }
  }
  
  template <int DisplacementDim>
  void HydroMechanicsProcess<DisplacementDim>::initializeConcreteProcess(
      NumLib::LocalToGlobalIndexMap const& dof_table,
      MeshLib::Mesh const& mesh,
      unsigned const integration_order)
  {
      const int mechanical_process_id = _use_monolithic_scheme ? 0 : 1;
      const int deformation_variable_id = _use_monolithic_scheme ? 1 : 0;
      ProcessLib::HydroMechanics::createLocalAssemblers<
          DisplacementDim, HydroMechanicsLocalAssembler>(
          mesh.getDimension(), mesh.getElements(), dof_table,
          // use displacement process variable to set shape function order
          getProcessVariables(mechanical_process_id)[deformation_variable_id]
              .get()
              .getShapeFunctionOrder(),
          _local_assemblers, mesh.isAxiallySymmetric(), integration_order,
          _process_data);
  
      _secondary_variables.addSecondaryVariable(
          "sigma_xx",
          makeExtrapolator(1, getExtrapolator(), _local_assemblers,
                           &LocalAssemblerInterface::getIntPtSigmaXX));
  
      _secondary_variables.addSecondaryVariable(
          "sigma_yy",
          makeExtrapolator(1, getExtrapolator(), _local_assemblers,
                           &LocalAssemblerInterface::getIntPtSigmaYY));
  
      _secondary_variables.addSecondaryVariable(
          "sigma_zz",
          makeExtrapolator(1, getExtrapolator(), _local_assemblers,
                           &LocalAssemblerInterface::getIntPtSigmaZZ));
  
      _secondary_variables.addSecondaryVariable(
          "sigma_xy",
          makeExtrapolator(1, getExtrapolator(), _local_assemblers,
                           &LocalAssemblerInterface::getIntPtSigmaXY));
  
      if (DisplacementDim == 3)
      {
          _secondary_variables.addSecondaryVariable(
              "sigma_xz",
              makeExtrapolator(1, getExtrapolator(), _local_assemblers,
                               &LocalAssemblerInterface::getIntPtSigmaXZ));
  
          _secondary_variables.addSecondaryVariable(
              "sigma_yz",
              makeExtrapolator(1, getExtrapolator(), _local_assemblers,
                               &LocalAssemblerInterface::getIntPtSigmaYZ));
      }
  
      _secondary_variables.addSecondaryVariable(
          "epsilon_xx",
          makeExtrapolator(1, getExtrapolator(), _local_assemblers,
                           &LocalAssemblerInterface::getIntPtEpsilonXX));
  
      _secondary_variables.addSecondaryVariable(
          "epsilon_yy",
          makeExtrapolator(1, getExtrapolator(), _local_assemblers,
                           &LocalAssemblerInterface::getIntPtEpsilonYY));
  
      _secondary_variables.addSecondaryVariable(
          "epsilon_zz",
          makeExtrapolator(1, getExtrapolator(), _local_assemblers,
                           &LocalAssemblerInterface::getIntPtEpsilonZZ));
  
      _secondary_variables.addSecondaryVariable(
          "epsilon_xy",
          makeExtrapolator(1, getExtrapolator(), _local_assemblers,
                           &LocalAssemblerInterface::getIntPtEpsilonXY));
  
      if (DisplacementDim == 3)
      {
          _secondary_variables.addSecondaryVariable(
              "epsilon_xz",
              makeExtrapolator(1, getExtrapolator(), _local_assemblers,
                               &LocalAssemblerInterface::getIntPtEpsilonXZ));
  
          _secondary_variables.addSecondaryVariable(
              "epsilon_yz",
              makeExtrapolator(1, getExtrapolator(), _local_assemblers,
                               &LocalAssemblerInterface::getIntPtEpsilonYZ));
      }
  
      _secondary_variables.addSecondaryVariable(
          "velocity",
          makeExtrapolator(mesh.getDimension(), getExtrapolator(),
                           _local_assemblers,
                           &LocalAssemblerInterface::getIntPtDarcyVelocity));
  
      _process_data.pressure_interpolated =
          MeshLib::getOrCreateMeshProperty<double>(
              const_cast<MeshLib::Mesh&>(mesh), "pressure_interpolated",
              MeshLib::MeshItemType::Node, 1);
  
      // Initialize local assemblers after all variables have been set.
      GlobalExecutor::executeMemberOnDereferenced(
          &LocalAssemblerInterface::initialize, _local_assemblers,
          *_local_to_global_index_map);
  }
  
  template <int DisplacementDim>
  void HydroMechanicsProcess<DisplacementDim>::initializeBoundaryConditions()
  {
      if (_use_monolithic_scheme)
      {
          const int process_id_of_hydromechancs = 0;
          initializeProcessBoundaryConditionsAndSourceTerms(
              *_local_to_global_index_map, process_id_of_hydromechancs);
          return;
      }
  
      // Staggered scheme:
      // for the equations of pressure
      const int hydraulic_process_id = 0;
      initializeProcessBoundaryConditionsAndSourceTerms(
          *_local_to_global_index_map_with_base_nodes, hydraulic_process_id);
  
      // for the equations of deformation.
      const int mechanical_process_id = 1;
      initializeProcessBoundaryConditionsAndSourceTerms(
          *_local_to_global_index_map, mechanical_process_id);
  }
  
  template <int DisplacementDim>
  void HydroMechanicsProcess<DisplacementDim>::assembleConcreteProcess(
      const double t, GlobalVector const& x, GlobalMatrix& M, GlobalMatrix& K,
      GlobalVector& b)
  {
      DBUG("Assemble the equations for HydroMechanics");
  
      // Note: This assembly function is for the Picard nonlinear solver. Since
      // only the Newton-Raphson method is employed to simulate coupled HM
      // processes in this class, this function is actually not used so far.
  
      std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
          dof_table = {std::ref(*_local_to_global_index_map)};
      // Call global assembler for each local assembly item.
      GlobalExecutor::executeMemberDereferenced(
          _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
          dof_table, t, x, M, K, b, _coupled_solutions);
  }
  
  template <int DisplacementDim>
  void HydroMechanicsProcess<DisplacementDim>::
      assembleWithJacobianConcreteProcess(const double t, GlobalVector const& x,
                                          GlobalVector const& xdot,
                                          const double dxdot_dx,
                                          const double dx_dx, GlobalMatrix& M,
                                          GlobalMatrix& K, GlobalVector& b,
                                          GlobalMatrix& Jac)
  {
      std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
          dof_tables;
      // For the monolithic scheme
      if (_use_monolithic_scheme)
      {
          DBUG(
              "Assemble the Jacobian of HydroMechanics for the monolithic"
              " scheme.");
          dof_tables.emplace_back(*_local_to_global_index_map);
      }
      else
      {
          // For the staggered scheme
          if (_coupled_solutions->process_id == 0)
          {
              DBUG(
                  "Assemble the Jacobian equations of liquid fluid process in "
                  "HydroMechanics for the staggered scheme.");
          }
          else
          {
              DBUG(
                  "Assemble the Jacobian equations of mechanical process in "
                  "HydroMechanics for the staggered scheme.");
          }
          dof_tables.emplace_back(*_local_to_global_index_map_with_base_nodes);
          dof_tables.emplace_back(*_local_to_global_index_map);
      }
  
      const int process_id =
          _use_monolithic_scheme ? 0 : _coupled_solutions->process_id;
      ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
  
      GlobalExecutor::executeSelectedMemberDereferenced(
          _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
          _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, x,
          xdot, dxdot_dx, dx_dx, M, K, b, Jac, _coupled_solutions);
  
      auto copyRhs = [&](int const variable_id, auto& output_vector) {
          if (_use_monolithic_scheme)
          {
              transformVariableFromGlobalVector(b, variable_id, dof_tables[0],
                                                output_vector,
                                                std::negate<double>());
          }
          else
          {
              transformVariableFromGlobalVector(
                  b, 0, dof_tables[_coupled_solutions->process_id], output_vector,
                  std::negate<double>());
          }
      };
      if (_use_monolithic_scheme || _coupled_solutions->process_id == 0)
      {
          copyRhs(0, *_hydraulic_flow);
      }
      if (_use_monolithic_scheme || _coupled_solutions->process_id == 1)
      {
          copyRhs(1, *_nodal_forces);
      }
  }
  
  template <int DisplacementDim>
  void HydroMechanicsProcess<DisplacementDim>::preTimestepConcreteProcess(
      GlobalVector const& x, double const t, double const dt,
      const int process_id)
  {
      DBUG("PreTimestep HydroMechanicsProcess.");
  
      _process_data.dt = dt;
      _process_data.t = t;
  
      if (hasMechanicalProcess(process_id))
      {
          ProcessLib::ProcessVariable const& pv =
              getProcessVariables(process_id)[0];
          GlobalExecutor::executeSelectedMemberOnDereferenced(
              &LocalAssemblerInterface::preTimestep, _local_assemblers,
              pv.getActiveElementIDs(), *_local_to_global_index_map, x, t,
              dt);
      }
  }
  
  template <int DisplacementDim>
  void HydroMechanicsProcess<DisplacementDim>::postTimestepConcreteProcess(
      GlobalVector const& x, const double /*t*/, const double /*delta_t*/,
      const int process_id)
  {
      DBUG("PostTimestep HydroMechanicsProcess.");
      ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
      GlobalExecutor::executeSelectedMemberOnDereferenced(
          &LocalAssemblerInterface::postTimestep, _local_assemblers,
          pv.getActiveElementIDs(), getDOFTable(process_id), x);
  }
  
  template <int DisplacementDim>
  void HydroMechanicsProcess<DisplacementDim>::postNonLinearSolverConcreteProcess(
      GlobalVector const& x, const double t, const int process_id)
  {
      if (!hasMechanicalProcess(process_id))
      {
          return;
      }
  
      DBUG("PostNonLinearSolver HydroMechanicsProcess.");
      // Calculate strain, stress or other internal variables of mechanics.
      ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
      GlobalExecutor::executeSelectedMemberOnDereferenced(
          &LocalAssemblerInterface::postNonLinearSolver, _local_assemblers,
          pv.getActiveElementIDs(), getDOFTable(process_id), x, t,
          _use_monolithic_scheme);
  }
  
  template <int DisplacementDim>
  void HydroMechanicsProcess<DisplacementDim>::computeSecondaryVariableConcrete(
      const double t, GlobalVector const& x, const int process_id)
  {
      DBUG("Compute the secondary variables for HydroMechanicsProcess.");
      ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
      GlobalExecutor::executeSelectedMemberOnDereferenced(
          &LocalAssemblerInterface::computeSecondaryVariable, _local_assemblers,
          pv.getActiveElementIDs(), getDOFTable(process_id), t, x,
          _coupled_solutions);
  }
  
  template <int DisplacementDim>
  std::tuple<NumLib::LocalToGlobalIndexMap*, bool>
  HydroMechanicsProcess<DisplacementDim>::getDOFTableForExtrapolatorData() const
  {
      const bool manage_storage = false;
      return std::make_tuple(_local_to_global_index_map_single_component.get(),
                             manage_storage);
  }
  
  template <int DisplacementDim>
  NumLib::LocalToGlobalIndexMap const&
  HydroMechanicsProcess<DisplacementDim>::getDOFTable(const int process_id) const
  {
      if (hasMechanicalProcess(process_id))
      {
          return *_local_to_global_index_map;
      }
  
      // For the equation of pressure
      return *_local_to_global_index_map_with_base_nodes;
  }
  
  template class HydroMechanicsProcess<2>;
  template class HydroMechanicsProcess<3>;
  
  }  // namespace HydroMechanics
  }  // namespace ProcessLib
  