GCC Code Coverage Report

Directory: ./
File: Solver/eigenProblem.hpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 6 16 37.5%
Branches: 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 _EIGENPROBLEM_HPP
16 #define _EIGENPROBLEM_HPP
17
18 // System includes
19 #include <vector>
20 #include <set>
21 #include <unordered_map>
22 #include <ostream>
23
24 // External includes
25 #include <parmetis.h>
26 #include "Core/NoThirdPartyWarning/Mpi/mpi.hpp"
27 #include "Core/NoThirdPartyWarning/Petsc/mat.hpp"
28 #include "Core/NoThirdPartyWarning/Petsc/ao.hpp"
29 #include "Core/NoThirdPartyWarning/Petsc/ksp.hpp"
30 #include "Core/NoThirdPartyWarning/Petsc/pc.hpp"
31 #include "Core/NoThirdPartyWarning/Petsc/error.hpp"
32 #ifdef FELISCE_WITH_SLEPC
33 #include "Core/NoThirdPartyWarning/Slepc/slepceps.hpp"
34 #include "Core/NoThirdPartyWarning/Slepc/slepcsvd.hpp"
35 #endif
36
37 // Project includes
38 #include "Geometry/geometricMeshRegion.hpp"
39 #include "Geometry/geometricFaces.hpp"
40 #include "DegreeOfFreedom/listVariable.hpp"
41 #include "DegreeOfFreedom/supportDofMesh.hpp"
42 #include "DegreeOfFreedom/dof.hpp"
43 #include "DegreeOfFreedom/boundaryConditionList.hpp"
44 #include "DegreeOfFreedom/boundaryCondition.hpp"
45 #include "InputOutput/io.hpp"
46 #include "FiniteElement/elementMatrix.hpp"
47 #include "FiniteElement/elementVector.hpp"
48 #include "FiniteElement/elementField.hpp"
49 #include "FiniteElement/listCurrentFiniteElement.hpp"
50 #include "FiniteElement/currentFiniteElement.hpp"
51 #include "FiniteElement/listCurvilinearFiniteElement.hpp"
52 #include "FiniteElement/currentFiniteElementWithBd.hpp"
53 #include "Core/chrono.hpp"
54 #include "Core/felisceParam.hpp"
55 #include "Core/felisceTransient.hpp"
56 #include "PETScInterface/petscMatrix.hpp"
57 #include "PETScInterface/petscVector.hpp"
58
59 /*
60 WARNING !
61 EigenProblem class files need slepc library to work properly.
62 */
63
64 namespace felisce
65 {
66 /*!
67 \class EigenProblem
68 \authors E. Schenone
69 \date 10/12/2012
70 \brief Manage general functions of an eigenvalues problem.
71 */
72 class EigenProblem
73 {
74 public:
75 typedef GeometricMeshRegion::ElementType ElementType;
76
77 // Constructor
78 //============
79 EigenProblem();
80
81 virtual ~EigenProblem();
82
83 void initialize(const GeometricMeshRegion::Pointer& mesh, FelisceTransient::Pointer fstransient, MPI_Comm& comm);
84
85 1 void fixIdOfTheProblemSolver(const int i) {
86 1 m_idProblem = i;
87 1 }
88
89 //Define Physical Variable associate to the problem
90 //=================================================
91 void definePhysicalVariable(const std::vector<PhysicalVariable>& listVariable, const std::vector<std::size_t>& listNumComp);
92 // Compute Dof
93 //============
94 /*!
95 \brief compute degrees of freedom with support dof creation.
96 Build pattern of the matrix.
97 */
98 void computeDof(int numProc, int idProc);
99 //! for one element and a local suppordof id, find coordniate of support.
100 void getSupportCoordinate(felInt idElt, felInt idSupport, int iUnknown, Point& pt);
101 // Partition function
102 //===================
103
104 /*!
105 \brief After partitionning we obtain element repartition and we create mesh region and "region" support dof.
106 */
107 void cutMesh();
108 /*!
109 \brief use ParMetis to partition dof (line of the matrix).
110 */
111 void partitionDof(idx_t numPart, int rankPart, MPI_Comm comm);
112
113 // Build Mesh and SupportDofMesh local with partition arrays
114 //==========================================================
115 /*!
116 \brief Create the local mesh and local supportDofMesh-> Create the local to global mapping for the elements
117 One element is on rankPart if the majority of its dof is in the process.
118 All points of the mesh are duplicated.
119 \param[in] rankPart Rank of the current process.
120 */
121 void setLocalMeshAndSupportDofMesh(int rankPart);
122
123 // Allocate Matrix
124 //========================
125 /*!
126 \brief use the pattern define with dof to allocate memory for the matrix in CSR format.
127 */
128 void allocateMatrix();
129 void allocateSequentialSolution();
130
131 // Assemble Matrix
132 //================
133 //! Asssembly loop on domain element with current finite element. Ths function call following functions.
134 /*!
135 \brief elementary loop to build matrix.
136 This function uses specific functions which are define in derived problem classes.
137 */
138 void assembleMatrix();
139 void allocateArrayForAssembleMatrix();
140 void desallocateArrayForAssembleMatrix();
141 void setElemPoint(ElementType& eltType, felInt iel, std::vector<Point*>& elemPoint, std::vector<felInt>& elemIdPoint, felInt* ielSupportDof);
142 void setElemPoint(ElementType& eltType, felInt iel, std::vector<Point*>& elemPoint, std::vector<felInt>& elemIdPoint, std::vector<felInt>& vectorSupport);
143 /*!
144 \brief use elementary calculus to fill the global matrix.
145 */
146 void setValueMatrix(const felInt* ielSupportDof);
147 //! Functions use to compute elementary operator for current finite element.
148 //! Initialize element array and elemField object with currentFiniteElement object.
149 void defineFiniteElement(const ElementType& eltType);
150 void initElementArray();
151 ✗ virtual void initPerElementType() {}
152
153 1 virtual void initPerDomain(int label) {
154 1 m_currentLabel = label;
155 1 }
156
157
158 //! Update finite element with element vertices coordinates and compute operators.
159 ✗ virtual void computeElementArray(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint,
160 felInt& iel, FlagMatrixRHS flagMatrixRHS = FlagMatrixRHS::matrix_and_rhs) {
161 IGNORE_UNUSED_ELEM_POINT;
162 IGNORE_UNUSED_ELEM_ID_POINT;
163 IGNORE_UNUSED_IEL;
164 IGNORE_UNUSED_FLAG_MATRIX_RHS;
165 }
166 ✗ virtual void computeElementArrayBD(const std::vector<Point*>&, const std::vector<felInt>&, felInt&) {}
167
168 //! Asssembly loop on domain element with curvilinear finite element. This function calls following functions.
169 /*!
170 \brief elementary loop to build matrix.
171 This function uses specific functions which are define in derived problem classes.
172 */
173 void assembleMatrixBD();
174 //! Functions use to compute elementary operator for curvilinear finite element.
175 //! Initialize element array and elemField object with curvilinearFiniteElement object.
176 void defineCurvilinearFiniteElement(const ElementType& eltType);
177 void initElementArrayBD();
178 ✗ virtual void userElementInitBD() {}
179 ✗ virtual void initPerElementTypeBD() {}
180
181 //! Initiallize some paramater associate to a mesh region.
182 ✗ virtual void initPerDomainBD(int label) {
183 ✗ m_currentLabel = label;
184 }
185
186 //! Fill global matrix with elementary values.
187 /*!
188 \brief use elementary calculus to fill the global matrix.
189 */
190 void setValueMatrixBD(const felInt* ielSupportDof);
191 //!Allocate memory for assemble matrix and rhs.
192 void allocateArrayForAssembleMatrixBD();
193
194
195 //Vector and Matrix operators
196 //===========================
197
198 // m_Matrix[i] = coef*m_Matrix[i]
199 void scaleMatrix(felReal coef, int i=0);
200
201 // M = coef*M
202 void scaleMatrix(felReal coef, PetscMatrix& M);
203
204 // Print solution of the system solve.
205 void printSolution(int verbose = 0) const;
206
207 //Write in file matrix.m, matrix in matlab format.
208 void writeMatrixForMatlab(int verbose = 0) const;
209
210 //Write in file filename, matrix in matlab format.
211 void writeMatrixForMatlab(std::string const& fileName, PetscMatrix& matrix) const;
212
213 //Write in file filenam, std::vector in matlab format.
214 void writeVectorForMatlab(std::string const& fileName, PetscVector& vector) const;
215
216 // Print a Petsc Matrix.
217 void printMatrix(PetscMatrix& M);
218
219 // Print a Petsc Vector.
220 void printVector(PetscVector& V);
221
222 //Solve linear system
223 //===================
224
225 /*!
226 \brief solve the linear solver with Petsc methods.
227 */
228 virtual void buildSolver();
229 void solve();
230
231 virtual void buildSVD();
232 void solveSVD();
233
234
235 //Write solution
236 //==============
237 void fromVecToDoubleStar(double* solution, PetscVector& sol, int rank, int dimension, felInt size=0);
238 void fromDoubleStarToVec(double* solution, PetscVector& sol, felInt size);
239 void writeSolution(int rank, IO& io, double& time, int iteration);
240 void writeEnsightVector(double* solValue, int idIter, std::string varName);
241 void writeEnsightCase(const int numIt, const double dt, std::string varName, const double time0=0.);
242 void writeEnsightCase(const int numIt, const double dt, std::vector<std::string> varName, const double time0=0.);
243
244 /*!
245 \brief write mesh to read solution with specific mesh (P2, Q2,...) .
246 */
247 void writeGeoForUnknown(int iUnknown);
248
249 //void initTime(double& time,double& timeStep){ m_time = time; m_timeStep = timeStep;}
250 void clearMatrix(const int i=0);
251 void gatherSolution();
252
253 //Tools
254 //=====
255 ✗ virtual void readData(IO& io) {
256 IGNORE_UNUSED_ARGUMENT(io);
257 }
258
259 //! get coodinate of support dof associate to the dof.
260 void findCoordinateWithIdDof(felInt i, Point& pt);
261
262 // Access Functions
263 //=================
264 inline const GeometricMeshRegion::Pointer mesh() const {
265 return m_mesh;
266 }
267 inline GeometricMeshRegion::Pointer mesh() {
268 return m_mesh;
269 }
270
271 inline const GeometricMeshRegion::Pointer meshLocal() const {
272 return m_meshLocal;
273 }
274 inline GeometricMeshRegion::Pointer meshLocal() {
275 return m_meshLocal;
276 }
277
278 inline const std::vector<SupportDofMesh> & supportDofUnknownLocal() const {
279 return m_supportDofUnknownLocal;
280 }
281 inline std::vector<SupportDofMesh> & supportDofUnknownLocal() {
282 return m_supportDofUnknownLocal;
283 }
284
285 inline const ISLocalToGlobalMapping & mappingElem() const {
286 return m_mappingElem;
287 }
288 inline ISLocalToGlobalMapping & mappingElem() {
289 return m_mappingElem;
290 }
291
292 inline const ListVariable & listVariable() const {
293 return m_listVariable;
294 }
295 ✗ inline ListVariable & listVariable() {
296 ✗ return m_listVariable;
297 }
298
299 inline const ListUnknown & listUnknown() const {
300 return m_listUnknown;
301 }
302 inline ListUnknown & listUnknown() {
303 return m_listUnknown;
304 }
305
306 inline const Dof & dof() const {
307 return m_dof;
308 }
309 inline Dof & dof() {
310 return m_dof;
311 }
312
313 inline const AO & ao() const {
314 return m_ao;
315 }
316 inline AO & ao() {
317 return m_ao;
318 }
319
320 inline const felInt & numDof() const {
321 return m_numDof;
322 }
323 inline felInt & numDof() {
324 return m_numDof;
325 }
326
327 inline const felInt & numDofLocal() const {
328 return m_numDofLocal;
329 }
330 inline felInt & numDofLocal() {
331 return m_numDofLocal;
332 }
333
334
335 inline const felInt & numDofPerUnknown(felInt i) const {
336 return m_numDofUnknown[i];
337 }
338 inline felInt & numDofPerUnknown(felInt i) {
339 return m_numDofUnknown[i];
340 }
341
342 inline felInt & numDof(PhysicalVariable unknown) {
343 int iUnknown = m_listUnknown.getUnknownIdList(unknown);
344 return m_numDofUnknown[iUnknown];
345 }
346
347 inline felInt & numDofLocalPerUnknown(PhysicalVariable unknown) {
348 int iUnknown = m_listUnknown.getUnknownIdList(unknown);
349 return m_numDofLocalUnknown[iUnknown];
350 }
351
352 inline ISLocalToGlobalMapping & mappingDofLocalToDofGlobal(PhysicalVariable unknown) {
353 int iUnknown = m_listUnknown.getUnknownIdList(unknown);
354 return m_mappingIdLocalDofPerUnknownToIdGlobalDofPetsc[iUnknown];
355 }
356
357
358 inline const PetscVector& solution() const {
359 return m_sol;
360 }
361 inline PetscVector& solution() {
362 return m_sol;
363 }
364
365 inline const PetscMatrix& Matrix(const int i=0) const {
366 return m_Matrix[i];
367 }
368 inline PetscMatrix& Matrix(const int i=0) {
369 return m_Matrix[i];
370 }
371
372 inline const ListCurrentFiniteElement & listCurrentFiniteElement() const {
373 return m_listCurrentFiniteElement;
374 }
375 inline ListCurrentFiniteElement & listCurrentFiniteElement() {
376 return m_listCurrentFiniteElement;
377 }
378
379 inline const PetscVector& seqSolution() const {
380 return m_seqSol;
381 }
382 inline PetscVector& seqSolution() {
383 return m_seqSol;
384 }
385
386 inline int dimension() {
387 return m_mesh->domainDim();
388 }
389
390 protected:
391 MPI_Comm m_petscComm;
392 PetscVector nullPetscVector;
393 #ifdef FELISCE_WITH_SLEPC
394 EPS m_eps;
395 SVD m_svd;
396 #endif
397 FelisceTransient::Pointer m_fstransient = nullptr;
398 //! input mesh which define domain of work.
399 GeometricMeshRegion::Pointer m_mesh = nullptr;
400 int m_idProblem = 0;
401 int m_verbose = 0;
402 // Number of matrix use in formulation.
403 int m_numberOfMatrix = 1;
404 //! small region of the global domain, obtain with a decomposition of the degrees of freedom.
405 GeometricMeshRegion::Pointer m_meshLocal = nullptr;
406
407 //! global information about support of the degrees of freedom.
408 std::vector<SupportDofMesh> m_supportDofUnknown;
409 //! specific information about support of the degrees of freedom for the local region of work for a processor.
410 std::vector<SupportDofMesh> m_supportDofUnknownLocal;
411 //! object which contains unknown of the linear system.
412 ListUnknown m_listUnknown;
413 //! object which contains all variable in the problem.
414 ListVariable m_listVariable;
415 //! array which contains degree of freedom repartition (size = numDof, value = rank of the processor which keep in memory the dof).
416 std::vector<int> m_dofPart;
417 //! array which contains element repartition (size = numDof, value = rank of the processor which keep in memory the element).
418 std::vector<int> m_eltPart;
419 //! object which manage degrees of freedom.
420 Dof m_dof;
421 //! solution of the problem.
422 PetscVector m_sol;
423 //! the whole solution of the problem gathered in the processor
424 PetscVector m_seqSol;
425 //! tells whenever sequential solution is allocated or not
426 bool m_seqSolAllocated = 0;
427 //! pointer to external vectors (e.g. vectors defined in another linearProblem)
428 std::vector<PetscVector*> m_externalVec;
429 std::vector<AO> m_externalAO;
430 std::vector<Dof*> m_externalDof;
431 //! number of dof contains by the current processor.
432 felInt m_numDofLocal = 0;
433 //! number of dof in the problem.
434 felInt m_numDof = 0;
435 //! number of dof for unknown
436 felInt* m_numDofUnknown = nullptr;
437 //! number of dof local for unknown
438 felInt* m_numDofLocalUnknown = nullptr;
439 //! mapping between global problem numbering and matrix dof numbering (specific to Petsc utilisation).
440 AO m_ao;
441
442 //! mapping between local to global ordering of dof.
443 ISLocalToGlobalMapping m_mappingNodes;
444 //! mapping between local to global ordering of element.
445 ISLocalToGlobalMapping m_mappingElem;
446 //! mapping between local to global ordering of element.
447 std::vector<ISLocalToGlobalMapping> m_mappingIdLocalDofPerUnknownToIdGlobalDofPetsc;
448
449 //! list of finite element in the problem.
450 ListCurrentFiniteElement m_listCurrentFiniteElement;
451
452 //! list of curvilinear finite element in the problem.
453 ListCurvilinearFiniteElement m_listCurvilinearFiniteElement;
454
455 int m_currentLabel;
456
457 //!Element matrix.
458 std::vector<ElementMatrix::Pointer> m_elementMat;
459 //!Element matrix boundary element.
460 std::vector<ElementMatrix::Pointer> m_elementMatBD;
461 std::vector<PetscMatrix> m_Matrix;
462
463 //!Temporary arrays use to assemble matrix
464 felInt* m_GposLine = nullptr;
465 felInt* m_GposColumn = nullptr;
466 std::vector<double*> m_valueMat;
467 //!Temporary arrays use to assemble matrix
468 felInt* m_GposVec = nullptr;
469 double* _GblockVec = nullptr;
470 felInt* m_ia = nullptr;
471
472 std::vector<felInt> m_localSizeOfVector;
473 felInt m_shiftValue = 0;
474 bool m_initNumDofUnknown = false;
475 bool m_buildSystem = false;
476 bool m_buildSolver = false;
477
478 std::string nameOfTheProblem = "No-name";
479
480 void pushBackExternalVec(PetscVector& extVec) {
481 m_externalVec.push_back(&extVec);
482 }
483 void pushBackExternalAO(AO extAO) {
484 m_externalAO.push_back(extAO);
485 }
486 void pushBackExternalDof(Dof& extDof) {
487 m_externalDof.push_back(&extDof);
488 }
489 };
490 }
491
492 #endif
493