Formulation.h

// GmshDDM - Copyright (C) 2019-2022, A. Royer, C. Geuzaine, Université de Liège
//
// See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/gmsh/ddm/issues

#ifndef H_GMSHDDM_FORMULATION
#define H_GMSHDDM_FORMULATION

#include "InterfaceCompoundField.h"
#include "InterfaceField.h"
#include "SubdomainCompoundField.h"
#include "SubdomainField.h"

#include <gmshfem/Formulation.h>
#include <unordered_map>
#include <vector>

typedef struct _p_Mat *Mat;
typedef struct _p_Vec *Vec;
typedef struct _p_KSP *KSP;

namespace gmshddm
{


  namespace problem
  {

    template< class T_Scalar >
    int MatVectProductA(Mat A, Vec X, Vec Y);
    template< class T_Scalar >
    int MatVectProductI_A(Mat A, Vec X, Vec Y);
    template< class T_Scalar >
    int MatVectProductImpl(Mat A, Vec X, Vec Y, bool IA);
    template< class T_Scalar >
    int MatMatProductImpl(Mat A, Mat X, Mat Y, bool IA);
    template< class T_Scalar >
    int MatMatProductA(Mat A, Mat X, Mat Y);
    template< class T_Scalar >
    int MatMatProductI_A(Mat A, Mat X, Mat Y);

    template< class T_Scalar >
    class AbstractIterativeSolver;

    template< class T_Scalar >
    class Formulation
    {
     private:
      const std::string _name;
      std::vector< gmshfem::problem::Formulation< T_Scalar > * > _volume;
      std::vector< std::unordered_map< unsigned int, gmshfem::problem::Formulation< T_Scalar > * > > _surface;
      std::unordered_map< field::InterfaceFieldInterface< T_Scalar > *, field::InterfaceFieldInterface< T_Scalar > * > _interfaceFields;
      std::unordered_map< unsigned int, std::vector< unsigned long long > > _interfaceMapping;
      gmshfem::algebra::Vector< T_Scalar > _rhs;
      std::vector<gmshfem::algebra::Vector< T_Scalar >> _multi_rhs; // For many shots

      unsigned int _numberOfIterations;
      bool _physicalCommutator;
      bool _artificialCommutator;
      bool _IA;
      std::set< unsigned int > _physicalSourceTerms;
      std::set< unsigned int > _artificialSourceTerms;
      std::vector< std::set< unsigned int > > _numberedPhysicalSourceTerms;
      std::vector< bool > _activeShots;


      enum class sourceStatus {
        ACTIVE,
        INACTIVE,
        NOT_IN_NUMBERED_TERMS
      };
      sourceStatus indexStatusInShots(unsigned id) const;



     protected:
      void _runIterationOperations(const unsigned int iteration, const gmshfem::scalar::Precision< T_Scalar > rnorm);
      gmshfem::algebra::Vector< T_Scalar > _sortDofValues(const gmshfem::algebra::Vector< T_Scalar > &vector, const unsigned int fieldTag);
      void _fillG(const T_Scalar *array);
      void _fromPetscToInterfaceFields(Vec g);

      void _preProVolume();
      void _infoNumberOfDofs() const;
      void _assembleAllVolume();
      void _solveAllVolume();
      void _assembleAndSolveSurface();
      void _extractRHS();

      void _activateBilinear (const unsigned long long idom);
      void _deactivateBilinear (const unsigned long long idom);

      void _buildPetscRHS(Vec* RHS);
      void _buildPetscAllRHS(Mat* RHS);

     public:
      Formulation(const std::string &name, const std::vector< std::vector< unsigned int > > &topology);
      Formulation(const std::string &name, const std::vector< gmshfem::problem::Formulation< T_Scalar > * > &volume, const std::vector< std::unordered_map< unsigned int, gmshfem::problem::Formulation< T_Scalar > * > > &surface);
      ~Formulation();

      unsigned int size() const;
      unsigned int size(unsigned int i) const;

      const gmshfem::function::ScalarFunction< T_Scalar > physicalSource(const gmshfem::function::ScalarFunction< T_Scalar > &f);
      const gmshfem::function::VectorFunction< T_Scalar > physicalSource(const gmshfem::function::VectorFunction< T_Scalar > &f);
      void physicalSourceTerm(const unsigned int termTag);
      const gmshfem::function::ScalarFunction< T_Scalar > artificialSource(const gmshfem::function::ScalarFunction< T_Scalar > &f);
      const gmshfem::function::VectorFunction< T_Scalar > artificialSource(const gmshfem::function::VectorFunction< T_Scalar > &f);
      void artificialSourceTerm(const unsigned int termTag);
      void togglePhysicalAndArtificialSourceTerms();

      // Multi source options
      void setShotNumber(unsigned N);
      void enableShot(unsigned i);
      void disableShot(unsigned i);
      void toggleShot(unsigned i);
      bool isShotEnabled(unsigned i) const;
      void disableAllShots();
      void numberedPhysicalSourceTerm(const unsigned int id, const unsigned int termTag);

      // TMP
      void compareRHSexport();


      void addInterfaceField(field::InterfaceFieldInterface< T_Scalar > &interfaceField);
      void addInterfaceField(field::InterfaceFieldInterface< T_Scalar > &interfaceField, field::InterfaceFieldInterface< T_Scalar > &interfaceFieldMapping);
      field::InterfaceFieldInterface< T_Scalar > *getInterfaceField(const std::string &name) const;

      gmshfem::problem::Formulation< T_Scalar > &operator()(unsigned int i);
      gmshfem::problem::Formulation< T_Scalar > &operator()(unsigned int i, unsigned int j);

      gmshfem::common::Timer pre(bool mustAssemble = true);
      gmshfem::common::Timer assemble(int storeID = -1);

      gmshfem::common::Timer solve(const std::string &solver, const double tolerance = 1e-6, const int iterMax = 1000, const bool sameMatrixWithArtificialAndPhysicalSources = false, const bool skipFinalSolutionComputation = false);
      gmshfem::common::Timer solve(AbstractIterativeSolver<T_Scalar> &solver, const double tolerance = 1e-6, const int iterMax = 1000, const bool sameMatrixWithArtificialAndPhysicalSources = false, const bool skipFinalSolutionComputation = false);
      gmshfem::algebra::MatrixCCS< T_Scalar > computeMatrix();
      gmshfem::algebra::MatrixCCS< T_Scalar > computeMatrixBlock();

      void setSolutionIntoInterfaceFields(const std::vector< T_Scalar > &solution);
      const gmshfem::algebra::Vector< T_Scalar > &getRHS() const;
      gmshfem::algebra::Vector< T_Scalar > getGlobalRHS(); //Built on the fly, so not a reference. Modifies internally the Petsc representation of the local RHS.


      unsigned int numberOfIterations() const;

      friend int MatVectProductA< T_Scalar >(Mat A, Vec X, Vec Y);
      friend int MatVectProductI_A< T_Scalar >(Mat A, Vec X, Vec Y);
      friend int MatMatProductI_A< T_Scalar >(Mat A, Mat X, Mat Y);
      friend int MatMatProductA< T_Scalar >(Mat A, Mat X, Mat Y);
      friend int MatVectProductImpl< T_Scalar >(Mat A, Vec X, Vec Y, bool IA);
      friend int MatMatProductImpl< T_Scalar >(Mat A, Mat X, Mat Y, bool IA);


      // T_Scalar, PetscScalar, T_Interger
      template< class T_Scalar1, class T_Scalar2, class T_Integer >
      struct PetscInterface {

        static const T_Scalar1 *arrayInterface(const T_Scalar2 *array, const T_Integer size)
        {
          T_Scalar1 *ret = new T_Scalar1[size];
#pragma omp parallel for
          for(T_Integer i = 0; i < size; ++i) {
            ret[i] = T_Scalar1(array[i]);
          }
          return ret;
        }

        static void freeArray(const T_Scalar1 *array)
        {
          delete[] array;
        }
      };

    };


  } // namespace problem


} // namespace gmshddm

#endif // H_GMSHDDM_FORMULATION