1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
|
#ifndef __BASESOLVER_H_INCLUDED__
#define __BASESOLVER_H_INCLUDED__
#include <deal.II/lac/trilinos_solver.h> //Trillinos direct sparse solver (parallel direct solver Paraklete comes with Trillinos 7.0 or higher)
#include "../commonIncludeFiles.h"
using namespace dealii;
using namespace std;
class BaseSolver{ //abstract class
public:
const int QUAD_POINTS_NUM_2D = 5;//number of quadrature points in each direction (QUAD_POINTS_NUM_2D^2 points per cell)
QGauss<2> quadrature_formula;
//constructor
BaseSolver():
quadrature_formula(QUAD_POINTS_NUM_2D),
dof_handler(triangulation) //initializing DoFHandler<2> member
{
cout << "entering BaseSolver constructor" << endl;
cout << "leaving BaseSolver constructor" << endl;
}; //constructor
virtual ~BaseSolver(){};
public:
static const bool debugOn = true;
//setting up linear solver
SolverControl solverControl;
TrilinosWrappers::SolverDirect linearSolver = TrilinosWrappers::SolverDirect(solverControl);//pass MUMPS or other parallel sover in 2nd argument line "Amesos_Mumps");
void setFlowParameter(cmdParametersStruct cmdParameters){
this->Reynolds = cmdParameters.Reynolds;
this->nu = 1 / Reynolds;
this->timeStep = cmdParameters.baseTimeStep / std::pow(2.0, cmdParameters.TimestepRefineCoeff);
};
//Flow parameters
double Reynolds;
double nu;
double timeStep;
const unsigned int dim = 2;
const FEValuesExtractors::Vector velocityExtractor = FEValuesExtractors::Vector(0);
const FEValuesExtractors::Scalar pressureExtractor = FEValuesExtractors::Scalar(dim);
ConstraintMatrix constraints;
virtual void setup_system();
virtual void setupBoundaries(SparseMatrixType& system_matrix);
virtual std::map<types::global_dof_index, double> getBoundaryValues(unsigned int curBoundaryIndicator)=0; //pure virtual f-n
//boundary variables
// std::map<types::global_dof_index, double> boundary_values;
std::vector<types::boundary_id> boundary_indicators;
// parallel::distributed::Triangulation<2,2> triangulation;
Triangulation<2,2> triangulation;
// ...variables for the sparsity pattern and values of the system matrix
// resulting from the discretization of the Laplace equation...
SparsityPattern sparsity_pattern;
SparseMatrixType system_matrix;
// ...and variables which will hold the right hand side and solution
// vectors.
Vector<double> solution;
Vector<double> system_rhs;
// FE_Q<2,2> feSystem;
DoFHandler<2> dof_handler;
FEValuesViews::Vector<2, 2> feValuesVelocity;
void setup_system_common(FiniteElement<2,2>& finiteElementSystem){ // } //of setup_system_common
TrilinosWrappers::SparseMatrix& convertDealIISparseToTrillinosMatrix(SparseMatrixType& system_matrix);
Vector<double>& solveLinearSystem(SparseMatrixType& system_matrix,const Vector<double>& system_rhs){
TrilinosWrappers::SparseMatrix& solverCompatibleSystem_matrix = convertDealIISparseToTrillinosMatrix(system_matrix);
// TrilinosWrappers::Vector* trilinosSolution;
Vector<double>* dealIIVector;
linearSolver.solve(solverCompatibleSystem_matrix, (*dealIIVector), system_rhs); //using TrilinosWrappers::SolverDirect::solve
// Vector<double> dealIIVector = new Vector<double>(trilinosSolution); //convert from Trilinos Vector to dealII
return (*dealIIVector);
};
}; // end of class
void BaseSolver::setupBoundaries(SparseMatrixType& system_matrix){
}
TrilinosWrappers::SparseMatrix& BaseSolver::convertDealIISparseToTrillinosMatrix(SparseMatrixType& system_matrix){
};
#endif
|
