Select Git revision
Formulation.h
-
Boris Martin authoredBoris Martin authored
Formulation.h 7.16 KiB
// 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