GCC Code Coverage Report

Directory: ./
File: Solver/eigenProblemALP.hpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 0 17 0.0%
Branches: 0 10 0.0%
Line Branch Exec Source
1 // ______ ______ _ _ _____ ______
2 // | ____| ____| | (_)/ ____| | ____|
3 // | |__ | |__ | | _| (___ ___| |__
4 // | __| | __| | | | |\___ \ / __| __|
5 // | | | |____| |____| |____) | (__| |____
6 // |_| |______|______|_|_____/ \___|______|
7 // Finite Elements for Life Sciences and Engineering
8 //
9 // License: LGL2.1 License
10 // FELiScE default license: LICENSE in root folder
11 //
12 // Main authors: E. Schenone
13 //
14
15 #ifndef _EIGENPROBLEMALP_HPP
16 #define _EIGENPROBLEMALP_HPP
17
18 // System includes
19
20 // External includes
21
22 // Project includes
23 #include "Core/felisceParam.hpp"
24 #include "Solver/eigenProblem.hpp"
25
26 /*
27 WARNING !
28 EigenProblemALP class files need slepc library to work properly.
29 */
30
31 namespace felisce
32 {
33 /*!
34 \class EigenProblemALP
35 \authors E. Schenone
36 \date 31/01/2013
37 \brief Manage functions of a ALP eigenvalues problem and ALP-ROM solver.
38 */
39
40 class EigenProblemALP:
41 public EigenProblem {
42 public:
43 // Constructor
44 //============
45 EigenProblemALP();
46 ~EigenProblemALP() override;
47
48 virtual void initialize(const GeometricMeshRegion::Pointer& mesh, FelisceTransient::Pointer fstransient, MPI_Comm& comm);
49
50 //Define Physical Variable associate to the problem
51 //=================================================
52 void initPerElementType() override;
53
54 // Assemble Matrix
55 //================
56 //! Update finite element with element vertices coordinates and compute operators.
57 void computeElementArray(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint, felInt& iel, FlagMatrixRHS flagMatrixRHS = FlagMatrixRHS::matrix_and_rhs) override;
58
59 //! Visualize modes as a variable
60 virtual void writeMode(const int iter=0);
61 virtual void writeEnsightSolution(const int iter);
62
63 //! Solve eigenproblem to build rom basis
64 void buildSolver() override;
65 virtual void solve();
66
67 // Reduced model functions
68 //========================
69 // Initialization
70 void readData(IO& io) override;
71 void getFiberDirection(felInt iel, int iUnknown, std::vector<double>& elemFiber);
72 void readBasis();
73 virtual void initializeROM();
74 virtual void initializeSystemSolver();
75 void readEnsightFile(double* vec, std::string varName, int idIter, felInt size);
76 void readEigenValueFromFile();
77 virtual void computeGamma();
78 void computeGammaExtrap();
79 void computeMatrixZeta();
80 virtual void computeMatrixM();
81 virtual void initializeSolution();
82 ✗ void setIntegrationMethod(const int methodModel) {
83 ✗ m_method = methodModel;
84 }
85 virtual void setHeteroTauClose(std::vector<double>& valueTauClose);
86 virtual void setFhNf0(std::vector<double>& valuef0);
87 virtual void setIapp(std::vector<double>& iapp, int& idS);
88 // Compute solution
89 int secondOrderGlobalIndex(int& i, int& j);
90 int thirdOrderGlobalIndex(int& i, int& j, int& k);
91 int fourthOrderGlobalIndex(int& i, int& j, int& k, int& h);
92 virtual void computeTensor();
93 virtual void checkSolution(PetscVector& solution, PetscVector& refSol, double & error);
94 virtual void checkBasis(int size=0);
95 //Update functions
96 virtual void updateBasis();
97 virtual void setPotential();
98 void MGS(int size=0);
99 void MGS(std::vector<PetscVector>& extraSpace);
100 virtual void updateTensor();
101 //virtual
102 void updateBeta();
103 virtual void updateBetaMonolithic();
104 void updateBetaEI();
105 void updateBetaBdf2();
106 void ConjGrad(double* A, double* b, double* x, double tol, int maxIter, int n);
107 void solveOrthComp(double* A, double* b, double* x, double tol, int maxIter);
108 virtual void updateEigenvalue();
109
110 // Projection functions from FE space to RO one and viceversa
111 void projectOnBasis(PetscVector& vectorDof, double* vectorBasis, bool flagMass=true, int size=0);
112 void projectOnDof(double* vectorBasis, PetscVector& vectorDof, int size);
113
114 // Write an ECG
115 ✗ virtual void initializeECG() {}
116 ✗ virtual void writeECG(int /*iteration*/) {}
117 ✗ virtual void updateEcgOperator() {}
118
119 ✗ virtual void readMeasureEcgFilter() {}
120
121 ✗ virtual void getExternalVec(double*){}
122 ✗ virtual void computeControlForFEM(PetscVector){}
123 ✗ virtual void readDataControl(){}
124 ✗ virtual bool applyFilter(){return false;}
125
126 // Print&Auxiliary functions
127 //==========================
128 void viewALP(double* vToPrint, int vSize, std::ofstream outStream);
129 void viewALP(double* vToPrint, int vSize, std::string fName);
130
131
132 // Acces functions
133 //================
134
135 PetscVector& basisFct(const int i=0) {
136 return m_basis[i];
137 }
138 double eigenValueFirst() {
139 return m_eigenValue[0];
140 }
141 double* eigenValue() {
142 return m_eigenValue;
143 }
144 double* beta() {
145 return m_beta;
146 }
147 double* mu() {
148 return m_mu;
149 }
150 double* gamma() {
151 return m_gamma;
152 }
153 double* eta() {
154 return m_eta;
155 }
156 int dimRomBasis() {
157 return m_dimRomBasis;
158 }
159 int tensor3size() {
160 return m_size3;
161 }
162 int tensor4size() {
163 return m_size4;
164 }
165
166 template <class Templ_functor> void evalFunctionOnDof(const Templ_functor& func, std::vector<double>& array);
167 template <class Templ_functor> void evalFunctionOnDof(const Templ_functor& func,double time, std::vector<double>& array);
168
169 protected:
170 int m_problem = -1;
171 int m_method = -1;
172 CurrentFiniteElement* m_fePotTransMemb = nullptr;
173 felInt m_ipotTransMemb;
174 PetscVector m_U_0;
175 PetscVector m_U_0_seq;
176 PetscVector m_W_0;
177 PetscVector m_W_0_seq;
178 PetscVector m_Ue_0;
179 PetscVector m_Ue_0_seq;
180 PetscVector m_FhNf0;
181 ElementField m_elemFieldU0;
182 ElementField m_elemFieldFhNf0;
183 PetscVector m_initPot;
184 double* m_eigenValue = nullptr;
185 std::vector<PetscVector> m_basis;
186 std::vector<PetscVector> m_orthComp;
187
188 double m_coefChi;
189 felInt m_dimRomBasis = 0;
190 bool m_useImprovedRec;
191 felInt m_dimOrthComp;
192
193 int m_idL = 0;
194 int m_idG = 0;
195 int m_idK = 0;
196 int m_idKie = 0;
197 int m_idGs = 0;
198 int m_idKs = 0;
199
200 int m_tensorOrder;
201 double** m_matrixM = nullptr;
202 std::vector<PetscVector> m_zeta;
203 double* m_tensorB = nullptr;
204 double* m_tensorE = nullptr;
205 double* m_tensorEs = nullptr;
206 double* m_tensorQ = nullptr;
207 double* m_tensorT = nullptr;
208 double* m_tensorTs = nullptr;
209 double** m_tensorU = nullptr;
210 double* m_tensorY = nullptr;
211 int m_size2;
212 int m_size3;
213 int m_size4;
214 int m_numSource;
215 double* m_beta = nullptr;
216 double* m_mu = nullptr;
217 double* m_xi = nullptr;
218 double* m_beta_n_1 = nullptr;
219 double* m_mu_n_1 = nullptr;
220 double* m_xi_n_1 = nullptr;
221 double* m_gamma = nullptr;
222 double* m_eta = nullptr;
223 double* m_gamma_extrap = nullptr;
224 double* m_eta_extrap = nullptr;
225 double** m_source = nullptr;
226 double* m_modeMean = nullptr;
227 double* m_fiber = nullptr;
228 bool m_solutionInitialized = false;
229
230 PetscMatrix m_matrixRom;
231 PetscVector m_rhsRom;
232 PetscVector m_solutionRom;
233 KSPInterface::Pointer m_kspRom = felisce::make_shared<KSPInterface>();
234 };
235
236 ✗ template <class Templ_functor> void EigenProblemALP::evalFunctionOnDof(const Templ_functor& func, std::vector<double>& array) {
237 Point pt;
238 ✗ array.resize(m_numDof);
239 felInt position;
240 ✗ for( felInt i = 0; i < m_numDof; i++) {
241 ✗ findCoordinateWithIdDof(i,pt);
242 ✗ position = i;
243 ✗ AOApplicationToPetsc(m_ao,1,&position);
244 ✗ array[position] = func(pt.x(),pt.y(),pt.z());
245 }
246 }
247
248 template <class Templ_functor> void EigenProblemALP::evalFunctionOnDof(const Templ_functor& func,double time, std::vector<double>& array) {
249 Point pt;
250 array.resize(m_numDof);
251 felInt position;
252 for( felInt i = 0; i < m_numDof; i++) {
253 findCoordinateWithIdDof(i,pt);
254 position = i;
255 AOApplicationToPetsc(m_ao,1,&position);
256 array[position] = func(pt.x(),pt.y(),pt.z(),time);
257 }
258 }
259
260 }
261
262 #endif
263