Directory: | ./ |
---|---|
File: | Solver/linearProblem.hpp |
Date: | 2024-04-14 07:32:34 |
Exec | Total | Coverage | |
---|---|---|---|
Lines: | 88 | 134 | 65.7% |
Branches: | 11 | 24 | 45.8% |
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: J. Foulon & J-F. Gerbeau & V. Martin | ||
13 | // | ||
14 | |||
15 | #ifndef LINEARPROBLEM_HPP | ||
16 | #define LINEARPROBLEM_HPP | ||
17 | |||
18 | // System includes | ||
19 | #include <unordered_set> | ||
20 | #include <vector> | ||
21 | #include <numeric> | ||
22 | |||
23 | // External includes | ||
24 | #include <parmetis.h> | ||
25 | #include "Core/NoThirdPartyWarning/Petsc/ao.hpp" | ||
26 | #include "Core/NoThirdPartyWarning/Petsc/error.hpp" | ||
27 | |||
28 | // Project includes | ||
29 | #include "PETScInterface/petscVector.hpp" | ||
30 | #include "PETScInterface/petscMatrix.hpp" | ||
31 | #include "PETScInterface/KSPInterface.hpp" | ||
32 | #include "PETScInterface/SNESInterface.hpp" | ||
33 | #include "Core/felisce.hpp" | ||
34 | #include "Core/shared_pointers.hpp" | ||
35 | #include "Core/felisceTools.hpp" | ||
36 | #include "Core/mpiInfo.hpp" | ||
37 | #include "Core/felisceTransient.hpp" | ||
38 | #include "DegreeOfFreedom/boundaryConditionList.hpp" | ||
39 | #include "DegreeOfFreedom/listVariable.hpp" | ||
40 | #include "DegreeOfFreedom/listUnknown.hpp" | ||
41 | #include "DegreeOfFreedom/supportDofMesh.hpp" | ||
42 | #include "DegreeOfFreedom/dof.hpp" | ||
43 | #include "DegreeOfFreedom/csrmatrixpattern.hpp" | ||
44 | #include "DegreeOfFreedom/dofBoundary.hpp" | ||
45 | #include "FiniteElement/elementVector.hpp" | ||
46 | #include "FiniteElement/elementMatrix.hpp" | ||
47 | #include "FiniteElement/listCurrentFiniteElement.hpp" | ||
48 | #include "FiniteElement/listCurrentFiniteElementWithBd.hpp" | ||
49 | #include "FiniteElement/listCurvilinearFiniteElement.hpp" | ||
50 | #include "Geometry/geometricMeshRegion.hpp" | ||
51 | #include "Geometry/point.hpp" | ||
52 | #include "InputOutput/io.hpp" | ||
53 | #include "Solver/CVGraphInterface.hpp" | ||
54 | #ifdef FELISCE_WITH_CVGRAPH | ||
55 | #include "Model/cvgMainSlave.hpp" | ||
56 | #endif | ||
57 | |||
58 | namespace felisce | ||
59 | { | ||
60 | ///@name felisce Globals | ||
61 | ///@{ | ||
62 | |||
63 | ///@} | ||
64 | ///@name Enum's | ||
65 | ///@{ | ||
66 | |||
67 | enum class FlagMatrixRHS : int | ||
68 | { | ||
69 | matrix_and_rhs = 0, | ||
70 | only_matrix = 1, | ||
71 | only_rhs = 2 | ||
72 | }; | ||
73 | |||
74 | namespace DirichletApplicationOptions { | ||
75 | enum Values { | ||
76 | nonSymmetricPseudoElimination = 0, | ||
77 | symmetricPseudoElimination = 1, | ||
78 | penalization = 2 | ||
79 | }; | ||
80 | } | ||
81 | |||
82 | enum class ApplyNaturalBoundaryConditions { | ||
83 | no, | ||
84 | yes | ||
85 | }; | ||
86 | |||
87 | ///@} | ||
88 | ///@name Type Definitions | ||
89 | ///@{ | ||
90 | |||
91 | // Forward declarations. | ||
92 | class Model; | ||
93 | class FelisceParam; | ||
94 | class FelisceTransient; | ||
95 | // class IO; | ||
96 | class Bdf; | ||
97 | class LumpedModelBC; | ||
98 | class CardiacCycle; | ||
99 | class RISModel; | ||
100 | class ElementField; | ||
101 | class BoundaryCondition; | ||
102 | #ifdef FELISCE_WITH_CVGRAPH | ||
103 | class cvgMainSlave; | ||
104 | #endif | ||
105 | |||
106 | typedef ApplyNaturalBoundaryConditions ApplyNaturalBoundaryConditionsType; | ||
107 | |||
108 | ///@} | ||
109 | ///@name Functions | ||
110 | ///@{ | ||
111 | |||
112 | ///@} | ||
113 | ///@name felisce Classes | ||
114 | ///@{ | ||
115 | |||
116 | // SNES staff | ||
117 | class SnesContext { | ||
118 | public: | ||
119 | SnesContext()= default;; | ||
120 | void initialize(LinearProblem& linearProblem, int rankProc, int sizeProc, | ||
121 | ApplyNaturalBoundaryConditionsType do_apply_natural, | ||
122 | FlagMatrixRHS flag); | ||
123 | |||
124 | LinearProblem* linearProblem; | ||
125 | int rankProc; | ||
126 | int sizeProc; | ||
127 | ApplyNaturalBoundaryConditionsType do_apply_natural; | ||
128 | FlagMatrixRHS flag_matrix_rhs; | ||
129 | }; | ||
130 | |||
131 | /*! | ||
132 | \class LinearProblem | ||
133 | \authors J. Foulon & J-F. Gerbeau & V. Martin | ||
134 | \date 05/01/2011 | ||
135 | \brief Manage general functions of an EDP problem. | ||
136 | */ | ||
137 | class LinearProblem | ||
138 | { | ||
139 | public: | ||
140 | ///@name Type Definitions | ||
141 | ///@{ | ||
142 | |||
143 | typedef GeometricMeshRegion::ElementType ElementType; | ||
144 | |||
145 | /// Two types of PetsVec, to improve readibility | ||
146 | enum vectorType { sequential, ///< Sequential vectors | ||
147 | parallel ///< Parallel vectors | ||
148 | }; | ||
149 | |||
150 | /// Pointer definition of LinearProblem | ||
151 | FELISCE_CLASS_POINTER_DEFINITION(LinearProblem); | ||
152 | |||
153 | ///@} | ||
154 | ///@name Life Cycle | ||
155 | ///@{ | ||
156 | |||
157 | /*! | ||
158 | * \brief Constructor | ||
159 | * | ||
160 | * \param[in] name Name of the linear problem | ||
161 | * \param[in] numMat Number of matrices | ||
162 | * \param[in] numVec Number of rhs | ||
163 | */ | ||
164 | LinearProblem(const std::string name = "No name", const std::size_t numMat = 1, const std::size_t numVec = 1); | ||
165 | |||
166 | /*! | ||
167 | * \brief Destructor | ||
168 | */ | ||
169 | virtual ~LinearProblem(); | ||
170 | |||
171 | ///@} | ||
172 | ///@name Operators | ||
173 | ///@{ | ||
174 | |||
175 | ///@} | ||
176 | ///@name Operations | ||
177 | ///@{ | ||
178 | |||
179 | /*! | ||
180 | * \brief Initialization steps. Beware: all quantities may not be initialised there. | ||
181 | * | ||
182 | * \param[in] mesh The vector of meshes | ||
183 | * \param[in] fstransient TODO | ||
184 | * \param[in] comm TODO | ||
185 | * \param[in] doUseSNES True if you handling a non linear problem and are willing to use Petsc Newton-Raphson method. False in most of the current cases! | ||
186 | */ | ||
187 | virtual void initialize(std::vector<GeometricMeshRegion::Pointer>& /*mesh*/, FelisceTransient::Pointer /*fstransient*/, MPI_Comm& /*comm*/, bool /*doUseSNES*/) = 0; | ||
188 | |||
189 | /*! | ||
190 | * \brief Init few quantities. | ||
191 | * | ||
192 | * Init a handful of attributes. | ||
193 | * | ||
194 | * This method is usually called in the beginning of the virtual #initialize() defined just above (\todo | ||
195 | * This design is error prone; it should be reversed: a non virtual initialize which includes a call to | ||
196 | * a virtual method that can be overloaded. Hence there would be no risk of forgetting the call of present method...) | ||
197 | * | ||
198 | * \param[in] doUseSNES True if you handling a non linear problem and are willing to use Petsc Newton-Raphson method. | ||
199 | */ | ||
200 | void initialize(std::vector<GeometricMeshRegion::Pointer>& mesh, MPI_Comm& comm, bool doUseSNES); | ||
201 | |||
202 | /*! | ||
203 | * \brief Finalizes the linear problem | ||
204 | */ | ||
205 | 16 | virtual void finalize(std::vector<GeometricMeshRegion::Pointer>& /*mesh*/, FelisceTransient::Pointer /*fstransient*/, MPI_Comm& /*comm*/) {} | |
206 | |||
207 | /*! | ||
208 | * \brief Initialize some attributes at the very end of the #Model::initializeLinearProblem() method. | ||
209 | * | ||
210 | * Rationale: for Elasticity and HyperElasticity, I need to set the values of attributes but can't do it | ||
211 | * in #initialize() because some of the quantities I need (namely the boundary conditions and the dof) aren't | ||
212 | * defined yet. | ||
213 | */ | ||
214 | 184 | virtual void InitializeDerivedAttributes() {} | |
215 | |||
216 | /*! | ||
217 | * \brief Set the ID of the problem solver | ||
218 | */ | ||
219 | 513 | inline void fixIdOfTheProblemSolver(int const id) { m_identifier_solver = id; }; | |
220 | |||
221 | |||
222 | |||
223 | // Define Physical Variable associate to the problem | ||
224 | //================================================== | ||
225 | void definePhysicalVariable(const std::vector<PhysicalVariable>& listVariable, const std::vector<std::size_t>& listNumComp); | ||
226 | |||
227 | 72 | virtual void userAddOtherVariables(std::vector<PhysicalVariable>& /*listVariable*/, std::vector<std::size_t>& /*listNumComp*/) {} | |
228 | |||
229 | 70 | virtual void userAddOtherUnknowns(std::vector<PhysicalVariable>& /*listVariable*/, std::vector<std::size_t>& /*listNumComp*/) {} | |
230 | |||
231 | |||
232 | |||
233 | // Compute Dof | ||
234 | //============ | ||
235 | |||
236 | /** | ||
237 | * @brief Compute global degrees of freedom with global support dof creation. | ||
238 | * Build pattern of the matrix. | ||
239 | */ | ||
240 | void computeDof(int numProc, int idProc); | ||
241 | |||
242 | /** | ||
243 | * @brief Change pattern function | ||
244 | */ | ||
245 | 437 | virtual void userChangePattern(int /*numProc*/, int /*rankProc*/) {}; | |
246 | |||
247 | /** | ||
248 | * @brief Initialize supprtDof in derived problem | ||
249 | */ | ||
250 | 511 | virtual void initSupportDofDerivedProblem() {}; | |
251 | |||
252 | /// for one element and a local suppordof id, find coordinate of support. | ||
253 | void getSupportCoordinate(ElementType& eltType, felInt iel, felInt idSupport, int iUnknown, Point& pt); | ||
254 | |||
255 | /// for one element and a local suppordof id, find coordinate of support. | ||
256 | void getSupportCoordinate(felInt idElt, felInt idSupport, int iUnknown, Point& pt); | ||
257 | |||
258 | |||
259 | |||
260 | // Partition function | ||
261 | //=================== | ||
262 | |||
263 | /** | ||
264 | * @brief After partitionning we obtain element repartition and we create mesh region and "region" support dof. | ||
265 | */ | ||
266 | void cutMesh(); | ||
267 | |||
268 | /** | ||
269 | * @brief use ParMetis to partition dof (line of the matrix). | ||
270 | * The function call parmetis to compute the dof partition which is stored in m_dofPart. | ||
271 | * Next, the application ordering (AO) used to fo from the global felisce ordering to the global petsc ordering | ||
272 | * is created. | ||
273 | * Finaly, the element partition is created and is store in m_eltPart. | ||
274 | * @param[in] numPart Number of process. | ||
275 | * @param[in] rankPart Rank of the current process. | ||
276 | * @param[in] comm MPI communicator. | ||
277 | */ | ||
278 | void partitionDof(idx_t numPart, int rankPart, MPI_Comm comm); | ||
279 | // void partitionDof_Petsc(idx_t numPart, int rankPart, MPI_Comm comm); | ||
280 | |||
281 | |||
282 | |||
283 | // Build Mesh and SupportDofMesh local with partition arrays | ||
284 | //========================================================== | ||
285 | |||
286 | /** | ||
287 | * @brief Create the local mesh and local supportDofMesh-> Create the local to global mapping for the elements | ||
288 | * One element is on rankPart if the majority of its dof is in the process. | ||
289 | * All points of the mesh are duplicated. | ||
290 | * @param[in] rankPart Rank of the current process. | ||
291 | */ | ||
292 | void setLocalMeshAndSupportDofMesh(int rankPart); | ||
293 | |||
294 | |||
295 | |||
296 | // Allocate Matrix and RHS | ||
297 | //======================== | ||
298 | |||
299 | /** | ||
300 | * @brief use the pattern defined with dof to allocate memory for the matrix in CSR format. | ||
301 | * This function is indeed allocating the matrix, but it is also creating the rhs, solution, serial solution | ||
302 | * vectors. It creates two mappings: | ||
303 | * - m_mappingNodes which std::unordered_map the local id to the global id of dofs. | ||
304 | * - m_mappingIdLocalDofPerUnknownToIdGlobalDofPetsc which is actually a std::vector of mappings that maps the | ||
305 | * local id to the global id in the petsc ordering for each unknown. | ||
306 | */ | ||
307 | void allocateMatrix(); | ||
308 | |||
309 | /** | ||
310 | * @brief Allocate the serial vector that contains the solution of the linear system. | ||
311 | */ | ||
312 | void allocateSequentialSolution(); | ||
313 | |||
314 | |||
315 | |||
316 | // Assemble Matrix | ||
317 | //================ | ||
318 | |||
319 | /** | ||
320 | * @brief elementary loop to build matrix. | ||
321 | * This function uses specific functions which are defined in derived problem classes. | ||
322 | */ | ||
323 | void assembleMatrixRHS(int rank, FlagMatrixRHS flagMatrixRHS = FlagMatrixRHS::matrix_and_rhs); | ||
324 | |||
325 | /** | ||
326 | * @brief Initialize the currentFiniteElement list. | ||
327 | */ | ||
328 | void defineFiniteElement(const ElementType& eltType); | ||
329 | |||
330 | /** | ||
331 | * @brief Initialize the currentFiniteElement list. | ||
332 | */ | ||
333 | void defineFiniteElement(const std::vector<ElementType>& eltType); | ||
334 | |||
335 | /** | ||
336 | * @brief Initialize element array | ||
337 | */ | ||
338 | void initElementArray(const bool initTranspose=false); | ||
339 | |||
340 | /** | ||
341 | * @brief Initialize the element fields | ||
342 | */ | ||
343 | ✗ | virtual void initPerElementType(ElementType /*eltType*/, FlagMatrixRHS /*flagMatrixRHS*/ = FlagMatrixRHS::matrix_and_rhs) {}; | |
344 | |||
345 | /** | ||
346 | * @brief Use by user to add some terms like source term define with elemField. | ||
347 | */ | ||
348 | 7387 | virtual void userElementInit() {}; | |
349 | |||
350 | /** | ||
351 | * @brief Initialize some parameter associated to a mesh region. | ||
352 | */ | ||
353 | virtual void initPerDomain(int label, FlagMatrixRHS flagMatrixRHS = FlagMatrixRHS::matrix_and_rhs); | ||
354 | |||
355 | /** | ||
356 | * @brief Get vertices, coordinates and support element of element | ||
357 | */ | ||
358 | void setElemPoint(ElementType& eltType, felInt iel, std::vector<Point*>& elemPoint, std::vector<felInt>& elemIdPoint, std::vector<felInt>& vectorSupport); | ||
359 | |||
360 | /** | ||
361 | * @brief Get vertices, coordinates and support element of element | ||
362 | */ | ||
363 | void setElemPoint(ElementType& eltType, felInt iel, std::vector<Point*>& elemPoint, std::vector<felInt>& elemIdPoint, felInt* ielSupportDof); | ||
364 | |||
365 | /**. | ||
366 | * @brief Update finite element with element vertices coordinates and compute operators. | ||
367 | */ | ||
368 | ✗ | virtual void computeElementArray(const std::vector<Point*>& /*elemPoint*/, const std::vector<felInt>& /*elemIdPoint*/, felInt& /*iel*/, FlagMatrixRHS /*flagMatrixRHS*/ = FlagMatrixRHS::matrix_and_rhs) {}; | |
369 | |||
370 | /**. | ||
371 | * @brief Update finite element with element vertices coordinates and compute operators. | ||
372 | */ | ||
373 | ✗ | virtual void computeElementArray(const std::vector<Point*>& /*elemPoint*/, const std::vector<felInt>& /*elemIdPoint*/, felInt& /*iel1*/, felInt& /*iel2*/, ElementType& /*eltType*/, felInt& /*ielGeo*/, FlagMatrixRHS /*flagMatrixRHS*/ = FlagMatrixRHS::matrix_and_rhs) {}; | |
374 | |||
375 | /**. | ||
376 | * @brief Update finite element with element vertices coordinates and compute operators for characteristics. | ||
377 | */ | ||
378 | ✗ | virtual void computeElementArrayCharact(const std::vector<Point*>& /*elemPoint*/, const std::vector<felInt>& /*elemIdPoint*/, felInt& /*iel*/, ElementType& /*eltType*/, felInt& /*ielGeo*/, FlagMatrixRHS /*flagMatrixRHS*/ = FlagMatrixRHS::matrix_and_rhs) {}; | |
379 | |||
380 | /**. | ||
381 | * @brief Compute user defined operators. | ||
382 | */ | ||
383 | 22641489 | virtual void userElementCompute(const std::vector<Point*>& /*elemPoint*/, const std::vector<felInt>& /*elemIdPoint*/, felInt& /*iel*/, FlagMatrixRHS /*flagMatrixRHS*/ = FlagMatrixRHS::matrix_and_rhs) {}; | |
384 | |||
385 | /**. | ||
386 | * @brief Compute user defined operators. | ||
387 | */ | ||
388 | 141128 | virtual void userElementCompute(const std::vector<Point*>& /*elemPoint*/, const std::vector<felInt>& /*elemIdPoint*/, felInt& /*iel1*/, felInt& /*iel2*/, ElementType& /*eltType*/, felInt& /*ielGeo*/, FlagMatrixRHS /*flagMatrixRHS*/ = FlagMatrixRHS::matrix_and_rhs) {}; | |
389 | |||
390 | /**. | ||
391 | * @brief Compute user defined operators. | ||
392 | */ | ||
393 | ✗ | virtual void userElementCompute(const std::vector<Point*>& /*elemPoint*/, const std::vector<felInt>& /*elemIdPoint*/, felInt& /*iel*/, const std::vector<Point*>& /*elemPointVol*/, std::vector<felInt>& /*elemIdPointVol*/, felInt& /*ielVol*/, std::vector<felInt>& /*elemIdPointLocal*/, FlagMatrixRHS /*flagMatrixRHS*/ = FlagMatrixRHS::matrix_and_rhs) {}; | |
394 | |||
395 | /** | ||
396 | * @brief use elementary calculus to fill the global matrix. | ||
397 | */ | ||
398 | void setValueMatrixRHS(const std::vector<felInt>& ielSupportDof1, const std::vector<felInt>& ielSupportDof2, FlagMatrixRHS flagMatrixRHS = FlagMatrixRHS::matrix_and_rhs, bool buildMatT=false); | ||
399 | |||
400 | /** | ||
401 | * @brief use elementary calculus to fill the global matrix. | ||
402 | */ | ||
403 | void setValueMatrixRHS(const felInt ielSupportDof1, const felInt ielSupportDof2, FlagMatrixRHS flagMatrixRHS = FlagMatrixRHS::matrix_and_rhs, bool buildMatT=false); | ||
404 | |||
405 | /** | ||
406 | * @brief use elementary calculus to fill the global matrix. | ||
407 | */ | ||
408 | void setValueCustomMatrix(const std::vector<felInt>& ielSupportDof1, const std::vector<felInt>& ielSupportDof2, PetscMatrix& mat, bool buildMatT = false); | ||
409 | |||
410 | /** | ||
411 | * @brief use elementary calculus to fill the global matrix. | ||
412 | */ | ||
413 | void setValueCustomMatrix(const felInt ielSupportDof1, const felInt ielSupportDof2, PetscMatrix& mat, bool buildMatT = false); | ||
414 | |||
415 | /** | ||
416 | * @brief Function to set the elemVector[0] into a generic vector vec. | ||
417 | */ | ||
418 | void setValueCustomVector(const std::vector<felInt>& ielSupportDof1, InsertMode mode, PetscVector& vec); | ||
419 | |||
420 | /** | ||
421 | * @brief Function to set the elemVector[0] into a generic vector vec. | ||
422 | */ | ||
423 | void setValueCustomVector(const felInt ielSupportDof1, InsertMode mode, PetscVector& vec); | ||
424 | |||
425 | /** | ||
426 | * @brief Allocate memory for assemble matrix and rhs. | ||
427 | */ | ||
428 | void allocateArrayForAssembleMatrixRHS(FlagMatrixRHS flagMatrixRHS = FlagMatrixRHS::matrix_and_rhs); | ||
429 | |||
430 | /** | ||
431 | * @brief Desallocate memory for assemble matrix and rhs. | ||
432 | */ | ||
433 | void desallocateArrayForAssembleMatrixRHS(FlagMatrixRHS flagMatrixRHS = FlagMatrixRHS::matrix_and_rhs); | ||
434 | |||
435 | /** | ||
436 | * @brief Assembly loop on boundary element of a certain label in implicit LumpedModelBC | ||
437 | */ | ||
438 | ✗ | virtual void assembleMatrixImplLumpedModelBC(int /*rank*/, const std::size_t /*iLabel*/) {}; | |
439 | |||
440 | /** | ||
441 | * @brief Assembly loop on crack elements pre-defined | ||
442 | */ | ||
443 | 27861 | virtual void derivedAssembleMatrixCrackModel(int /*rank*/) {}; | |
444 | |||
445 | /** | ||
446 | * @brief elementary loop to build matrix. | ||
447 | * This function uses specific functions which are define in derived problem classes. | ||
448 | */ | ||
449 | void assembleMatrixRHSBD(int rank, FlagMatrixRHS flagMatrixRHS = FlagMatrixRHS::matrix_and_rhs); | ||
450 | |||
451 | /** | ||
452 | * @brief Initialize the curvilinearFiniteElement list. | ||
453 | */ | ||
454 | void defineCurvilinearFiniteElement(const ElementType& eltType); | ||
455 | |||
456 | /** | ||
457 | * @brief Initialize the curvilinearFiniteElement list. | ||
458 | */ | ||
459 | void defineCurvilinearFiniteElement(const std::vector<ElementType>& eltType); | ||
460 | |||
461 | /** | ||
462 | * @brief Initialize the currentFiniteElementWithBd list. | ||
463 | */ | ||
464 | void defineCurrentFiniteElementWithBd(const ElementType& eltType); | ||
465 | |||
466 | /** | ||
467 | * @brief Initialize the currentFiniteElementWithBd list. | ||
468 | */ | ||
469 | void defineCurrentFiniteElementWithBd(const std::vector<ElementType>& eltType); | ||
470 | |||
471 | /** | ||
472 | * @brief Initialize element array | ||
473 | */ | ||
474 | void initElementArrayBD(const bool initTranspose=false); | ||
475 | |||
476 | /** | ||
477 | * @brief Initialize the element fields | ||
478 | */ | ||
479 | ✗ | virtual void initPerElementTypeBD(ElementType /*eltType*/, FlagMatrixRHS /*flagMatrixRHS*/ = FlagMatrixRHS::matrix_and_rhs) {}; | |
480 | |||
481 | /** | ||
482 | * @brief Use by user to add some terms like source term define with elemField. | ||
483 | */ | ||
484 | 8748 | virtual void userElementInitBD() {}; | |
485 | |||
486 | /** | ||
487 | * @brief Initialize some parameter associated to a mesh region. | ||
488 | */ | ||
489 | virtual void initPerDomainBD(int label, FlagMatrixRHS flagMatrixRHS = FlagMatrixRHS::matrix_and_rhs); | ||
490 | |||
491 | /** | ||
492 | * @brief Get vertices, coordinates, normal, tangent and support element of element | ||
493 | */ | ||
494 | void setElemPointNormalTangent(ElementType& eltType, felInt iel,std::vector<Point*>& elemPoint, std::vector<Point*>& elemNormal, | ||
495 | std::vector< std::vector<Point*> >& elemTangent, std::vector<felInt>& elemIdPoint, std::vector<felInt>& vectorSupport); | ||
496 | |||
497 | /**. | ||
498 | * @brief Update finite element with element vertices coordinates and compute operators. | ||
499 | */ | ||
500 | ✗ | virtual void computeElementArrayBD(const std::vector<Point*>& /*elemPoint*/, const std::vector<felInt>& /*elemIdPoint*/, felInt& /*iel*/, FlagMatrixRHS /*flagMatrixRHS*/ = FlagMatrixRHS::matrix_and_rhs) {}; | |
501 | |||
502 | /**. | ||
503 | * @brief Update finite element with element vertices coordinates and compute operators. | ||
504 | */ | ||
505 | ✗ | virtual void computeElementArrayBD(const std::vector<Point*>& /*elemPoint*/, const std::vector<felInt>& /*elemIdPoint*/, const std::vector<Point*>& /*elemNormal*/, const std::vector<std::vector<Point*> >& /*elemTangent*/, felInt& /*iel*/, FlagMatrixRHS /*flagMatrixRHS*/ = FlagMatrixRHS::matrix_and_rhs) {}; | |
506 | |||
507 | /**. | ||
508 | * @brief Compute user defined operators. | ||
509 | */ | ||
510 | 200236 | virtual void userElementComputeBD(const std::vector<Point*>& /*elemPoint*/, const std::vector<felInt>& /*elemIdPoint*/, felInt& /*iel*/) {}; | |
511 | |||
512 | /** | ||
513 | * @brief use elementary calculus to fill the global matrix. | ||
514 | */ | ||
515 | void setValueMatrixRHSBD(const std::vector<felInt>& ielSupportDof1, const std::vector<felInt>& ielSupportDof2, FlagMatrixRHS flagMatrixRHS = FlagMatrixRHS::matrix_and_rhs, bool buildMatT=false); | ||
516 | |||
517 | /** | ||
518 | * @brief use elementary calculus to fill the global matrix. | ||
519 | */ | ||
520 | void setValueMatrixRHSBD(const felInt ielSupportDof1, const felInt ielSupportDof2, FlagMatrixRHS flagMatrixRHS = FlagMatrixRHS::matrix_and_rhs, bool buildMatT=false); | ||
521 | |||
522 | /** | ||
523 | * @brief Place the m_elementMatrixBD[0] into a generic matrix. | ||
524 | */ | ||
525 | void setValueCustomMatrixBD(const std::vector<felInt>& ielSupportDof1, const std::vector<felInt>& ielSupportDof2, PetscMatrix& mat, bool buildMatT=false); | ||
526 | |||
527 | /** | ||
528 | * @brief Place the m_elementMatrixBD[0] into a generic matrix. | ||
529 | */ | ||
530 | void setValueCustomMatrixBD(const felInt ielSupportDof1, const felInt ielSupportDof2, PetscMatrix& mat, bool buildMatT=false); | ||
531 | |||
532 | /** | ||
533 | * @brief Place the m_elementVectorBD[0] into a generic vec. | ||
534 | */ | ||
535 | void setValueCustomVectorBD(const std::vector<felInt>& ielSupportDof1, InsertMode mode, PetscVector& vec) ; | ||
536 | |||
537 | /** | ||
538 | * @brief Place the m_elementVectorBD[0] into a generic vec. | ||
539 | */ | ||
540 | void setValueCustomVectorBD(const felInt ielSupportDof1, InsertMode mode, PetscVector& vec ); | ||
541 | |||
542 | /** | ||
543 | * @brief Allocate memory for assemble matrix and rhs. | ||
544 | */ | ||
545 | void allocateArrayForAssembleMatrixRHSBD(FlagMatrixRHS flagMatrixRHS = FlagMatrixRHS::matrix_and_rhs); | ||
546 | |||
547 | /** | ||
548 | * @brief elementary loop to build matrix for "mixed" domain (volumes+surfaces or surfaces+edges). | ||
549 | * This function is used when the 3D (resp. 2D) domain has some surfaces (resp. edges) that are not boundaries of any volume (resp. surface), see for instance atria+ventricles. | ||
550 | * This function uses specific functions which are defined in derived problem classes. | ||
551 | */ | ||
552 | void assembleMatrixRHSCurrentAndCurvilinearElement(int rank, FlagMatrixRHS flagMatrixRHS = FlagMatrixRHS::matrix_and_rhs); | ||
553 | |||
554 | /** | ||
555 | * @brief This a quite general function to have a general assembly loop on the boundary. | ||
556 | */ | ||
557 | template<class theClass> | ||
558 | void assemblyLoopBoundaryGeneral( void (theClass::*functionOfTheLoop)( felInt ), | ||
559 | const std::vector<felInt> & labels, void (theClass::*initPerElementType)(), | ||
560 | void (theClass::*updateFunction)(const std::vector<Point*>&,const std::vector<int>&)); | ||
561 | |||
562 | |||
563 | // Assemble Boundary Conditions | ||
564 | //============================= | ||
565 | |||
566 | /*! | ||
567 | * \brief Apply boundary counditions | ||
568 | * | ||
569 | * \param[in] doPrintBC If true (and verbose > 2) print informations about boundary conditions. | ||
570 | * Rationale: in a Newton-Raphson method, we do not want to print these informations for each iteration. | ||
571 | * \param[in] do_apply_natural ApplyNaturalBoundaryConditions::yes in most cases. | ||
572 | * Rationale: in Newmark time scheme, we do not want to apply thee ones in the time iterations. | ||
573 | * | ||
574 | */ | ||
575 | void applyBC(const int & applicationOption, int rank, FlagMatrixRHS flagMatrixRHSEss = FlagMatrixRHS::matrix_and_rhs, FlagMatrixRHS flagMatrixRHSNat = FlagMatrixRHS::matrix_and_rhs, const int idBCforLinCombMethod=0, | ||
576 | bool doPrintBC = true, ApplyNaturalBoundaryConditionsType do_apply_natural = ApplyNaturalBoundaryConditions::yes); | ||
577 | |||
578 | /** | ||
579 | * @brief Assembly loop for boundary conditions on boundary element with curvilinear finite element. | ||
580 | */ | ||
581 | void assembleMatrixRHSNaturalBoundaryCondition(int rank, FlagMatrixRHS flagMatrixRHS=FlagMatrixRHS::matrix_and_rhs, const int idBCforLinCombMethod=0); | ||
582 | |||
583 | /** | ||
584 | * @brief Allocate vector boundary condition for derived problems | ||
585 | */ | ||
586 | 6104 | virtual void allocateVectorBoundaryConditionDerivedLinPb() {}; | |
587 | |||
588 | /** | ||
589 | * @brief Initialize the element fields | ||
590 | */ | ||
591 | 728 | virtual void initPerElementTypeBoundaryCondition(ElementType& /*eltType*/, FlagMatrixRHS /*flagMatrixRHS*/) {}; | |
592 | |||
593 | /** | ||
594 | * @brief Use by user to add some terms | ||
595 | */ | ||
596 | 5420 | virtual void userElementInitNaturalBoundaryCondition() {}; | |
597 | |||
598 | /** | ||
599 | * @brief Allocate vector m_elemFieldRobin, m_elemFieldRobinNormal, m_elemFieldNeumann, m_elemFieldNeumannNormal, m_elemFieldBackflowStab. | ||
600 | */ | ||
601 | void allocateElemFieldBoundaryCondition(const int idBCforLinCombMethod = 0); | ||
602 | |||
603 | /** | ||
604 | * @brief Overwrite value of m_elemFieldRobin, m_elemFieldRobinNormal, m_elemFieldNeumann and m_elemFieldNeumannNormal in derived linear problem. | ||
605 | */ | ||
606 | 3016 | virtual void allocateElemFieldBoundaryConditionDerivedLinPb(const int /*idBCforLinCombMethod*/) {}; | |
607 | |||
608 | /** | ||
609 | * @brief Initiallize some paramater associate to a mesh region. | ||
610 | */ | ||
611 | 5736 | virtual void initPerDomainBoundaryCondition(std::vector<Point*>& /*elemPoint*/, const std::vector<felInt>& /*elemIdPoint*/, int label, felInt /*numEltPerLabel*/, felInt* /*ielSupportDof*/, ElementType& /*eltType*/, felInt /*numElement_eltType*/, FlagMatrixRHS /*flagMatrixRHS*/ = FlagMatrixRHS::matrix_and_rhs) { m_currentLabel = label; }; | |
612 | |||
613 | /** | ||
614 | * @brief Update finite element with element vertices coordinates and compute operators (for curvilinear finite element). | ||
615 | */ | ||
616 | virtual void computeElementArrayBoundaryCondition(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint,felInt& iel, FlagMatrixRHS flagMatrixRHS = FlagMatrixRHS::matrix_and_rhs); | ||
617 | |||
618 | /** | ||
619 | * @brief Compute user defined boundary conditions. | ||
620 | */ | ||
621 | 345554 | virtual void userElementComputeNaturalBoundaryCondition(const std::vector<Point*>& /*elemPoint*/, const std::vector<felInt>& /*elemIdPoint*/, felInt& /*iel*/, int /*label*/) {}; | |
622 | |||
623 | /** | ||
624 | * @brief Compute user defined boundary conditions. | ||
625 | */ | ||
626 | ✗ | virtual void userElementComputeNaturalBoundaryCondition(const std::vector<Point*>& /*elemPoint*/, const std::vector<felInt>& /*elemIdPoint*/, felInt& /*ielSupportDof1*/, felInt& /*ielSupportDof2*/, felInt& /*ielGeoGlobal*/, int /*label*/) {}; | |
627 | |||
628 | /** | ||
629 | * @brief Apply boundary condition of type: Neumann, NeumannNormal | ||
630 | */ | ||
631 | void applyNaturalBoundaryCondition(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint, felInt& iel, FlagMatrixRHS flagMatrixRHS = FlagMatrixRHS::only_rhs); | ||
632 | |||
633 | /** | ||
634 | * @brief Compute and apply boundary conditions | ||
635 | */ | ||
636 | virtual void computeAndApplyElementNaturalBoundaryCondition(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint, felInt& ielSupportDof1, felInt& ielSupportDof2, felInt& ielGeoGlobal, FlagMatrixRHS flagMatrixRHS = FlagMatrixRHS::matrix_and_rhs); | ||
637 | |||
638 | /**. | ||
639 | * @brief Compute operators for surface model | ||
640 | */ | ||
641 | ✗ | virtual void computeElementArraySurfaceModel(const std::vector<Point*>& /*elemPoint*/, const std::vector<felInt>& /*elemIdPoint*/, felInt& /*iel*/, FlagMatrixRHS /*flagMatrixRHS*/ = FlagMatrixRHS::matrix_and_rhs) {}; | |
642 | |||
643 | /**. | ||
644 | * @brief Essential boundary condition : Dirichlet, Essential LumpedModelBC | ||
645 | */ | ||
646 | void applyEssentialBoundaryCondition(int applicationOption, int rank, FlagMatrixRHS flagMatrixRHS); | ||
647 | |||
648 | /**. | ||
649 | * @brief Essential boundary condition derived problem | ||
650 | */ | ||
651 | 23353 | virtual void applyEssentialBoundaryConditionDerivedProblem(int /*rank*/, FlagMatrixRHS /*flagMatrixRHS*/) {}; | |
652 | |||
653 | /** | ||
654 | * @brief Define boundary condition and allocate element fields | ||
655 | */ | ||
656 | virtual void defineBC(); | ||
657 | |||
658 | /** | ||
659 | * @brief determine the dof for each essential boundary condition | ||
660 | */ | ||
661 | void determineDofAssociateToLabel(); | ||
662 | |||
663 | /** | ||
664 | * @brief Initialize the array idEltAndIdSupport for the boundary condition BC. | ||
665 | * @param[in] BC The boundary condition | ||
666 | */ | ||
667 | void determineElementSuppDofAssociateToBC(BoundaryCondition* BC); | ||
668 | |||
669 | /** | ||
670 | * @brief Set the value of the constant BC (in time AND space) | ||
671 | */ | ||
672 | void finalizeEssBCConstantInT(); | ||
673 | |||
674 | /** | ||
675 | * @brief Set the value of the constant BC in users (in time AND space) | ||
676 | */ | ||
677 | 664 | virtual void userFinalizeEssBCConstantInT() {}; | |
678 | |||
679 | /** | ||
680 | * @brief Call the functions for the essential BC | ||
681 | */ | ||
682 | void finalizeEssBCTransient(); | ||
683 | |||
684 | /** | ||
685 | * @brief Call the derived problem functions for the essential BC | ||
686 | */ | ||
687 | 27280 | virtual void finalizeEssBCTransientDerivedProblem() {}; | |
688 | |||
689 | /** | ||
690 | * @brief Call the user functions for the essential BC | ||
691 | */ | ||
692 | 23292 | virtual void userFinalizeEssBCTransient() {}; | |
693 | |||
694 | /** | ||
695 | * @brief ??? | ||
696 | */ | ||
697 | void getListDofOfBC(const BoundaryCondition& BC, std::unordered_set<felInt>& idDofBC, std::unordered_map<int, double>& idBCAndValue); | ||
698 | |||
699 | /** | ||
700 | * @brief ??? | ||
701 | */ | ||
702 | void getRings(); | ||
703 | |||
704 | /** | ||
705 | * @brief ??? | ||
706 | */ | ||
707 | virtual void userGetBC(BoundaryCondition* &BCIn,BoundaryCondition* &BCOut, BoundaryCondition* &BCLat); | ||
708 | |||
709 | /** | ||
710 | * @brief ??? | ||
711 | */ | ||
712 | void getListDofOfSurfaceByBc(const BoundaryCondition& BC, std::set<felInt>& idDofBC, int rankProc ); | ||
713 | |||
714 | /** | ||
715 | * @brief ??? | ||
716 | */ | ||
717 | ✗ | virtual double userComputeValueNeumannNormalTransient(const int /*iNeumannNormalBC*/) { return 0; }; | |
718 | |||
719 | /** | ||
720 | * @brief ??? | ||
721 | */ | ||
722 | void setValueBoundaryCondition(BoundaryCondition* bc, PetscVector& vecOfValue); | ||
723 | |||
724 | /** | ||
725 | * @brief ??? | ||
726 | */ | ||
727 | ✗ | virtual void defineRobinBCVector(PetscVector& /*robinVec*/, const int /*iRobin*/) {}; | |
728 | |||
729 | |||
730 | |||
731 | //Solve linear system | ||
732 | //=================== | ||
733 | |||
734 | /*! | ||
735 | * \brief Post assemble operations. | ||
736 | * | ||
737 | * Performs some operations after assembling is done. For instance, in dynamic systems the assembled matrices | ||
738 | * are not necessarily those of the system, but rather intermediate ones required to build the matrices of the system. | ||
739 | * For instance, in dynamic elasticity a new stiffness matrix is assembled and the matrix and RHS of the | ||
740 | * system are created from this one and others stored elsewhere. | ||
741 | * | ||
742 | * \internal I know Model class provides a postAssemblingMatrixRHS() method. However, it is not enough in the | ||
743 | * case of hyperelasticity: as Petsc's SNES is in charge of the Newton Raphson's loop, I must respect its interface | ||
744 | * and assembling is done at Petsc's call of the function given to SNESSetFunction(). | ||
745 | * I deliberately chose a very different method name to avoid any ambiguity. | ||
746 | */ | ||
747 | 17159 | virtual void computeSystemAfterAssembling() {}; | |
748 | |||
749 | /*! | ||
750 | * \brief Pre assemble operations for iterative solve. | ||
751 | * | ||
752 | * Perform some operations before assembling the matrices (for instance clearing them out.) | ||
753 | * | ||
754 | * \internal I know Model class provides a preAssemblingMatrixRHS() method. However, it is not enough in the | ||
755 | * case of hyperelasticity: as Petsc's SNES is in charge of the Newton Raphson's loop, I must respect its interface | ||
756 | * and assembling is done at Petsc's call of the function given to SNESSetFunction(). | ||
757 | * I deliberately chose a very different method name to avoid any ambiguity. | ||
758 | */ | ||
759 | void preSNESAssemble(); | ||
760 | |||
761 | /*! | ||
762 | * \brief Additional operations for preSNESAssemble() specific to derived problem. | ||
763 | */ | ||
764 | 12107 | virtual void preSNESAssembleAdditional() {}; | |
765 | |||
766 | /** | ||
767 | * @brief solve the linear solver with Petsc methods. | ||
768 | */ | ||
769 | void buildSolver(); | ||
770 | |||
771 | /** | ||
772 | * @brief Solve a ksp linear problem. | ||
773 | */ | ||
774 | void solve(int rank, int size, bool new_matrix = true); | ||
775 | |||
776 | /** | ||
777 | * @brief Solve with ksp a system with the same matrix as in the previous solve | ||
778 | */ | ||
779 | void solveWithSameMatrix(); | ||
780 | |||
781 | /*! | ||
782 | * \brief Solve a snes linear system. | ||
783 | * | ||
784 | * A SNES linear system is an iterative system that uses Petsc's Newton-Raphson method. | ||
785 | * | ||
786 | * \param[in] do_apply_natural Whether natural boundaries conditions are applied. | ||
787 | * ApplyNaturalBoundaryConditions::yes in most cases, but might be no for instance in dynamic systems where | ||
788 | * these conditions appear explicitly only on the very first iteration. | ||
789 | * \param[in] flag_matrix_rhs Specify on which element the assembling is done. | ||
790 | */ | ||
791 | void iterativeSolve(int rank, int size, ApplyNaturalBoundaryConditionsType do_apply_natural = ApplyNaturalBoundaryConditions::yes, FlagMatrixRHS flag_matrix_rhs = FlagMatrixRHS::matrix_and_rhs); | ||
792 | |||
793 | /** | ||
794 | * @brief Initializes the solver before the non-linear solve. | ||
795 | */ | ||
796 | 17471 | virtual void initNonLinearIteration() {}; | |
797 | |||
798 | /** | ||
799 | * @brief Finalizes the solver before the non-linear solve. | ||
800 | */ | ||
801 | 11680 | virtual void endNonLinearIteration() {}; | |
802 | |||
803 | /** | ||
804 | * @brief Initializes the solver for the current time step | ||
805 | */ | ||
806 | 5571 | virtual void initIteration() {}; | |
807 | |||
808 | /** | ||
809 | * @brief Finalizes the solver for the current time step | ||
810 | */ | ||
811 | 4205 | virtual void endIteration() {}; | |
812 | |||
813 | /** | ||
814 | * @brief Create and copy matrix and rhs withou boundary conditions | ||
815 | */ | ||
816 | void createAndCopyMatrixRHSWithoutBC(FlagMatrixRHS flag = FlagMatrixRHS::matrix_and_rhs); | ||
817 | |||
818 | /** | ||
819 | * @brief Compute the residual | ||
820 | */ | ||
821 | void computeResidual(); | ||
822 | |||
823 | |||
824 | |||
825 | //Vector and Matrix operators | ||
826 | //=========================== | ||
827 | |||
828 | /** | ||
829 | * @brief Add matrix and rhs | ||
830 | */ | ||
831 | 140 | virtual inline void addMatrixRHS() { matrix(0).axpy(1,matrix(1),SAME_NONZERO_PATTERN); } | |
832 | |||
833 | /** | ||
834 | * @brief Scales the LHS matrix by a constant. | ||
835 | */ | ||
836 | ✗ | virtual void addScaleMatrix(double /*coef*/ = 1.0) {}; | |
837 | |||
838 | /** | ||
839 | * @brief Scales the RHS vector by a constant. | ||
840 | */ | ||
841 | ✗ | virtual void addScaleRHS(double /*coef*/ = 1.0) {}; | |
842 | |||
843 | ✗ | virtual inline void copyMatrixRHS() {}; | |
844 | |||
845 | ✗ | virtual inline void inverseMatrix() {}; | |
846 | |||
847 | ✗ | virtual inline void solveExplicit() {}; | |
848 | |||
849 | /*! | ||
850 | * \brief Clear Matrix and/or RHS | ||
851 | * | ||
852 | * \param[in] flagMatrixRHS If 0, clear both matrix and RHS. If 1, only matrix, if 2, only RHS | ||
853 | */ | ||
854 | void clearMatrixRHS(const FlagMatrixRHS flagMatrixRHS = FlagMatrixRHS::matrix_and_rhs); | ||
855 | |||
856 | /*! | ||
857 | * \brief Clear the ith matrix. | ||
858 | * \param[in] i Matrix index | ||
859 | */ | ||
860 | void clearMatrix(std::size_t i = 0); | ||
861 | |||
862 | /*! | ||
863 | * \brief Clear the ith std::vector. | ||
864 | * \param[in] i Vector index | ||
865 | */ | ||
866 | void clearVector(std::size_t i = 0); | ||
867 | |||
868 | /** | ||
869 | * @brief This method provides the solution vector | ||
870 | */ | ||
871 | 24561 | inline void gatherSolution() { gatherVec(m_sol, m_seqSol); } | |
872 | |||
873 | /** | ||
874 | * @brief This method provides the NL iteration evaluation state vector | ||
875 | */ | ||
876 | 21834 | inline void gatherEvaluationState() { gatherVec(m_evaluationState, m_seqEvaluationState); }; | |
877 | |||
878 | // We should stop using this function and prefer directly gatherVec (petscVector.hpp, petscVector_friend.hpp). | ||
879 | // This function is left to avoid compatibility issues. | ||
880 | 50684 | inline void gatherVector(PetscVector& v, PetscVector& seqV) { gatherVec(v,seqV); } // TODO remove | |
881 | |||
882 | 260 | virtual void gatherVectorBeforeAssembleMatrixRHS() {}; | |
883 | |||
884 | void getSolution(double* & solution, int sizeProc, int numProc); | ||
885 | |||
886 | void setSolution(double* solution, int sizeProc, int numProc); | ||
887 | |||
888 | void getSolutionUnknown(PhysicalVariable unknown, PetscVector& v); | ||
889 | |||
890 | void getVector(double* & value, PetscVector& seqVec, int sizeProc, int numProc); | ||
891 | |||
892 | /** | ||
893 | * @brief This method gets the coodinate of support dof associate to the dof. | ||
894 | */ | ||
895 | void findCoordinateWithIdDof(felInt i, Point& pt); | ||
896 | |||
897 | /*! | ||
898 | * \brief Evaluates function f(x,y,z) on dof | ||
899 | * \param[in] func Functor to compute the values | ||
900 | * \param[in] array Array of values | ||
901 | */ | ||
902 | template <class Templ_functor> | ||
903 | void evalFunctionOnDof(const Templ_functor& func, std::vector<double>& array); // TODO is linearProblem a good place for this function? | ||
904 | |||
905 | /*! | ||
906 | * \brief Evaluates function f(x,y,z) on dof | ||
907 | * \param[in] func Functor to compute the values | ||
908 | * \param[in] array Array of values | ||
909 | * \param[in] extraData Extra data | ||
910 | */ | ||
911 | template <class Templ_functor> | ||
912 | void evalFunctionOnDof(const Templ_functor& func, std::vector<double>& array, std::vector<double>& extraData); // TODO is linearProblem a good place for this function? | ||
913 | |||
914 | /*! | ||
915 | * \brief Evaluates function f(x,y,z) on dof | ||
916 | * \param[in] func Functor to compute the values | ||
917 | * \param[in] array Array of values | ||
918 | * \param[in] extraData Extra data | ||
919 | */ | ||
920 | template <class Templ_functor> | ||
921 | void evalFunctionOnDof(const Templ_functor& func, std::vector<int>& array, std::vector<double>& extraData); // TODO is linearProblem a good place for this function? | ||
922 | |||
923 | /*! | ||
924 | * \brief Evaluates function f(x,y,z,t) on dof | ||
925 | * \param[in] func Functor to compute the values | ||
926 | * \param[in] time Time | ||
927 | * \param[in] array Array of values | ||
928 | */ | ||
929 | template <class Templ_functor> | ||
930 | void evalFunctionOnDof(const Templ_functor& func,double time, std::vector<double>& array); // TODO is linearProblem a good place for this function? | ||
931 | |||
932 | /*! | ||
933 | * \brief Evaluates function f(x,y,z,t) on dof | ||
934 | * \param[in] func Functor to compute the values | ||
935 | * \param[in] array Array of values | ||
936 | * \param[in] extraData Extra data | ||
937 | * \param[in] time Time | ||
938 | */ | ||
939 | template <class Templ_functor> | ||
940 | void evalFunctionOnDof(const Templ_functor& func, std::vector<double>& array, std::vector<double>& extraData,double time); // TODO is linearProblem a good place for this function? | ||
941 | |||
942 | /*! | ||
943 | * \brief Evaluates function f(x,y,z,t) on dof | ||
944 | * \param[in] func Functor to compute the values | ||
945 | * \param[in] time Time | ||
946 | * \param[in] array Array of values | ||
947 | * \param[in] extraData Extra data | ||
948 | */ | ||
949 | template <class Templ_functor> | ||
950 | void evalFunctionOnDof(const Templ_functor& func, double time, std::vector<double>& array, std::vector<double>& extraData); // TODO is linearProblem a good place for this function? | ||
951 | |||
952 | /*! | ||
953 | * \brief Evaluates function f(x,y,z) on dof | ||
954 | * \param[in] func Functor to compute the values | ||
955 | * \param[in] iComp I-th component | ||
956 | * \param[in] array Array of values | ||
957 | */ | ||
958 | template <class Templ_functor> | ||
959 | void evalFunctionOnDof(const Templ_functor& func, int iComp, std::vector<double>& array); // TODO is linearProblem a good place for this function? | ||
960 | |||
961 | /*! | ||
962 | * \brief Evaluates function f(x,y,z) on dof | ||
963 | * \param[in] func Functor to compute the values | ||
964 | * \param[in] array Array of values | ||
965 | */ | ||
966 | template <class Templ_functor> | ||
967 | void evalFunctionWithCompOnDof(const Templ_functor& func, std::vector<double>& array); // TODO is linearProblem a good place for this function? | ||
968 | |||
969 | /*! | ||
970 | * \brief Set values in m_seqSol and m_sol | ||
971 | * \param[in] functorXYZ Functor to compute the values | ||
972 | * \param[in] iVariable Variable ID | ||
973 | * \param[in] iComponent Component ID | ||
974 | */ | ||
975 | template <class FunctorXYZ> | ||
976 | inline void set_U_0(const FunctorXYZ& functorXYZ, int iVariable, int iComponent); | ||
977 | |||
978 | /*! | ||
979 | * \brief Set values in m_seqSol and m_sol | ||
980 | * \param[in] functorXYZ Functor to compute the values | ||
981 | * \param[in] iVariable Variable ID | ||
982 | * \param[in] iComponent Component ID | ||
983 | */ | ||
984 | template <class FunctorXYZ> | ||
985 | inline void set_U_0_parallel(const FunctorXYZ& functorXYZ, int iUnknown, int iComponent); | ||
986 | |||
987 | |||
988 | |||
989 | // Methods for auxiliary vectors | ||
990 | //============================== | ||
991 | |||
992 | /** | ||
993 | * @brief Initialize vectors | ||
994 | */ | ||
995 | void initPetscVectors(); | ||
996 | |||
997 | |||
998 | |||
999 | // Virtual methods widely used in derived linear problems | ||
1000 | //======================================================= | ||
1001 | ✗ | virtual void updateExtrapol(PetscVector& /*V_1*/) {}; | |
1002 | |||
1003 | ✗ | virtual void initExtrapol(PetscVector& /*V_1*/) {}; | |
1004 | |||
1005 | ✗ | virtual void initializeTimeScheme(Bdf* /*bdf*/) {}; | |
1006 | |||
1007 | ✗ | virtual void readData(IO& /*io*/) {}; | |
1008 | |||
1009 | ✗ | virtual void readData(IO& /*io*/, double /*iteration*/) {}; | |
1010 | |||
1011 | ///@} | ||
1012 | ///@name Access | ||
1013 | ///@{ | ||
1014 | |||
1015 | /** | ||
1016 | * @brief Returns the dimension of the problem (constant version) | ||
1017 | */ | ||
1018 | 57439957 | inline int dimension() const { return m_dimesion; } | |
1019 | |||
1020 | /** | ||
1021 | * @brief Returns the verbosity level (constant version) | ||
1022 | */ | ||
1023 | 581694 | inline const int& verbosity() const { return m_verbosity; } | |
1024 | |||
1025 | /** | ||
1026 | * @brief Returns the list of variables (constant version) | ||
1027 | */ | ||
1028 | inline const ListVariable& listVariable() const { return m_listVariable; } | ||
1029 | |||
1030 | /** | ||
1031 | * @brief Returns the list of variables | ||
1032 | */ | ||
1033 | 177714 | inline ListVariable& listVariable() { return m_listVariable; } | |
1034 | |||
1035 | /** | ||
1036 | * @brief Returns the list of unknowns (constant version) | ||
1037 | */ | ||
1038 | inline const ListUnknown& listUnknown() const { return m_listUnknown; } | ||
1039 | |||
1040 | /** | ||
1041 | * @brief Returns the list of unknowns | ||
1042 | */ | ||
1043 | 63207 | inline ListUnknown& listUnknown() { return m_listUnknown; } | |
1044 | |||
1045 | /** | ||
1046 | * @brief Returns the i-th linear problem mesh (constant version) | ||
1047 | */ | ||
1048 | inline const GeometricMeshRegion::Pointer mesh(int id) const { return m_mesh[id]; } | ||
1049 | |||
1050 | /** | ||
1051 | * @brief Returns the i-th linear problem mesh | ||
1052 | */ | ||
1053 | 4 | inline GeometricMeshRegion::Pointer mesh(int id) { return m_mesh[id]; } | |
1054 | |||
1055 | /** | ||
1056 | * @brief Returns the m_currentMesh linear problem mesh (constant version) | ||
1057 | */ | ||
1058 | ✗ | inline const GeometricMeshRegion::Pointer mesh() const { return m_mesh[m_currentMesh]; } | |
1059 | |||
1060 | /** | ||
1061 | * @brief Returns the m_currentMesh linear problem mesh | ||
1062 | */ | ||
1063 | 69222 | inline GeometricMeshRegion::Pointer mesh() { return m_mesh[m_currentMesh]; } | |
1064 | |||
1065 | /** | ||
1066 | * @brief Returns the linear problem mesh vector (constant version) | ||
1067 | */ | ||
1068 | inline const std::vector<GeometricMeshRegion::Pointer> & listMesh() const { return m_mesh; } | ||
1069 | |||
1070 | /** | ||
1071 | * @brief Returns the linear problem mesh vector | ||
1072 | */ | ||
1073 | inline std::vector<GeometricMeshRegion::Pointer> & listMesh() { return m_mesh; } | ||
1074 | |||
1075 | /** | ||
1076 | * @brief Returns the linear problem local mesh (constant version) | ||
1077 | */ | ||
1078 | inline const GeometricMeshRegion::Pointer meshLocal(int id) const { return m_meshLocal[id]; } | ||
1079 | |||
1080 | /** | ||
1081 | * @brief Returns the linear problem local mesh | ||
1082 | */ | ||
1083 | 824 | inline GeometricMeshRegion::Pointer meshLocal(int id) { return m_meshLocal[id]; } | |
1084 | |||
1085 | /** | ||
1086 | * @brief Returns the m_currentMesh linear problem local mesh (constant version) | ||
1087 | */ | ||
1088 | inline const GeometricMeshRegion::Pointer meshLocal() const { return m_meshLocal[m_currentMesh]; } | ||
1089 | |||
1090 | /** | ||
1091 | * @brief Returns the m_currentMesh linear problem local mesh | ||
1092 | */ | ||
1093 | 1966 | inline GeometricMeshRegion::Pointer meshLocal() { return m_meshLocal[m_currentMesh]; } | |
1094 | |||
1095 | /** | ||
1096 | * @brief Returns the linear problem local mesh (constant version) | ||
1097 | */ | ||
1098 | inline const std::vector<GeometricMeshRegion::Pointer>& listMeshLocal() const { return m_meshLocal; } | ||
1099 | |||
1100 | /** | ||
1101 | * @brief Returns the linear problem local mesh | ||
1102 | */ | ||
1103 | inline std::vector<GeometricMeshRegion::Pointer>& listMeshLocal() { return m_meshLocal; } | ||
1104 | |||
1105 | /** | ||
1106 | * @brief Returns the DoF object associated to the unknown (constant version). | ||
1107 | */ | ||
1108 | inline const Dof& dof() const { return m_dof; } | ||
1109 | |||
1110 | /** | ||
1111 | * @brief Returns the DoF object associated to the unknown. | ||
1112 | */ | ||
1113 | 90715950 | inline Dof& dof() { return m_dof; } | |
1114 | |||
1115 | /** | ||
1116 | * @brief Returns the dof partition vector (constant version) | ||
1117 | */ | ||
1118 | inline const std::vector<int>& dofPartition() const { return m_dofPart; } | ||
1119 | |||
1120 | /** | ||
1121 | * @brief Returns the dof partition vector | ||
1122 | */ | ||
1123 | 8 | inline std::vector<int>& dofPartition(){ return m_dofPart; } | |
1124 | |||
1125 | /** | ||
1126 | * @brief It returns the number of dofs(constant version) | ||
1127 | */ | ||
1128 | inline const felInt& numDof() const { return m_numDof; } | ||
1129 | |||
1130 | /** | ||
1131 | * @brief It returns the number of dofs | ||
1132 | */ | ||
1133 | 113360 | inline felInt& numDof() { return m_numDof; } | |
1134 | |||
1135 | /** | ||
1136 | * @brief It returns the number of dofs (local-constant version). | ||
1137 | */ | ||
1138 | inline const felInt& numDofLocal() const { return m_numDofLocal; } | ||
1139 | |||
1140 | /** | ||
1141 | * @brief It returns the number of dofs (local). | ||
1142 | */ | ||
1143 | ✗ | inline felInt& numDofLocal() { return m_numDofLocal; } | |
1144 | |||
1145 | /** | ||
1146 | * @brief It returns the number of dofs associated to the unknown. (constant version) | ||
1147 | * \param[in] i The ID of the unknown | ||
1148 | */ | ||
1149 | inline const felInt& numDofPerUnknown(felInt i) const { return m_numDofUnknown[i]; } | ||
1150 | |||
1151 | /** | ||
1152 | * @brief It returns the number of dofs associated to the unknown. | ||
1153 | * \param[in] i The ID of the unknown | ||
1154 | */ | ||
1155 | ✗ | inline felInt& numDofPerUnknown(felInt i) { return m_numDofUnknown[i]; } | |
1156 | |||
1157 | /** | ||
1158 | * @brief It returns the number of dofs associated to the unknown. (constant version) | ||
1159 | * \param[in] i The ID of the unknown | ||
1160 | */ | ||
1161 | inline felInt numDofPerUnknown(PhysicalVariable unknown) const { return m_numDofUnknown[m_listUnknown.getUnknownIdList(unknown)]; } | ||
1162 | |||
1163 | /** | ||
1164 | * @brief It returns the number of dofs associated to the unknown. | ||
1165 | * \param[in] i The ID of the unknown | ||
1166 | */ | ||
1167 | ✗ | inline felInt& numDofPerUnknown(PhysicalVariable unknown) { return m_numDofUnknown[m_listUnknown.getUnknownIdList(unknown)]; } | |
1168 | |||
1169 | /** | ||
1170 | * @brief It returns the local number of dofs associated to the unknown. (constant version) | ||
1171 | * \param[in] i The ID of the unknown | ||
1172 | */ | ||
1173 | inline felInt numDofLocalPerUnknown(PhysicalVariable unknown) const { return m_numDofLocalUnknown[m_listUnknown.getUnknownIdList(unknown)]; }; | ||
1174 | |||
1175 | /** | ||
1176 | * @brief It returns the local number of dofs associated to the unknown. (constant version) | ||
1177 | * \param[in] i The ID of the unknown | ||
1178 | */ | ||
1179 | 1164 | inline felInt& numDofLocalPerUnknown(PhysicalVariable unknown) { return m_numDofLocalUnknown[m_listUnknown.getUnknownIdList(unknown)]; }; | |
1180 | |||
1181 | /** | ||
1182 | * @brief Returns the vector of supportDof (constant version) | ||
1183 | */ | ||
1184 | inline const std::vector<SupportDofMesh>& supportDofUnknown() const { return m_supportDofUnknown; } | ||
1185 | |||
1186 | /** | ||
1187 | * @brief Returns the vector of supportDof | ||
1188 | */ | ||
1189 | ✗ | inline std::vector<SupportDofMesh>& supportDofUnknown() { return m_supportDofUnknown; } | |
1190 | |||
1191 | /** | ||
1192 | * @brief Returns the supportDof related to the given unknown (constant version) | ||
1193 | * \param[in] iunknown The ID of the unknown | ||
1194 | */ | ||
1195 | ✗ | inline const SupportDofMesh& supportDofUnknown(int iunknown) const { return m_supportDofUnknown[iunknown]; } | |
1196 | |||
1197 | /** | ||
1198 | * @brief Returns the supportDof related to the given unknown (constant version) | ||
1199 | * \param[in] iunknown The ID of the unknown | ||
1200 | */ | ||
1201 | 42893792 | inline SupportDofMesh& supportDofUnknown(int iunknown) { return m_supportDofUnknown[iunknown]; } | |
1202 | |||
1203 | /** | ||
1204 | * @brief Returns the vector of supportDofLocal (constant version) | ||
1205 | */ | ||
1206 | inline const std::vector<SupportDofMesh>& supportDofUnknownLocal() const { return m_supportDofUnknownLocal; } | ||
1207 | |||
1208 | /** | ||
1209 | * @brief Returns the vector of supportDofLocal | ||
1210 | */ | ||
1211 | inline std::vector<SupportDofMesh>& supportDofUnknownLocal() { return m_supportDofUnknownLocal; } | ||
1212 | |||
1213 | /** | ||
1214 | * @brief Returns the supportDofLocal related to the given unknown (constant version) | ||
1215 | * \param[in] iunknown The ID of the unknown | ||
1216 | */ | ||
1217 | inline const SupportDofMesh& supportDofUnknownLocal(int iunknown) const { return m_supportDofUnknownLocal[iunknown]; } | ||
1218 | |||
1219 | /** | ||
1220 | * @brief Returns the supportDofLocal related to the given unknown (constant version) | ||
1221 | * \param[in] iunknown The ID of the unknown | ||
1222 | */ | ||
1223 | inline SupportDofMesh& supportDofUnknownLocal(int iunknown) { return m_supportDofUnknownLocal[iunknown]; } | ||
1224 | |||
1225 | /** | ||
1226 | * @brief Returns the flag to know if m_ao is allocated (constant version) | ||
1227 | */ | ||
1228 | inline bool isAOAllocated() const { return m_initAO; } | ||
1229 | |||
1230 | /** | ||
1231 | * @brief Returns the flag to know if m_ao is allocated | ||
1232 | */ | ||
1233 | 11 | inline bool& isAOAllocated() { return m_initAO; } | |
1234 | |||
1235 | /** | ||
1236 | * @brief Returns the AO object (constant version) | ||
1237 | */ | ||
1238 | inline const AO& ao() const { return m_ao; } | ||
1239 | |||
1240 | /** | ||
1241 | * @brief Returns the AO object | ||
1242 | */ | ||
1243 | 350467 | inline AO& ao() { return m_ao; } | |
1244 | |||
1245 | /** | ||
1246 | * @brief Returns the flag to know if m_mappingNodes is allocated (constant version) | ||
1247 | */ | ||
1248 | inline bool isMappingNodesAllocated() const { return m_initMappingNodes; } | ||
1249 | |||
1250 | /** | ||
1251 | * @brief Returns the flag to know if m_mappingNodes is allocated | ||
1252 | */ | ||
1253 | 11 | inline bool& isMappingNodesAllocated() { return m_initMappingNodes; } | |
1254 | |||
1255 | /** | ||
1256 | * @brief Returns the local to global mapping of the dof (constant version) | ||
1257 | */ | ||
1258 | inline const ISLocalToGlobalMapping& mappingNodes() const { return m_mappingNodes; } | ||
1259 | |||
1260 | /** | ||
1261 | * @brief Returns the local to global mapping of the dof | ||
1262 | */ | ||
1263 | 11 | inline ISLocalToGlobalMapping& mappingNodes() { return m_mappingNodes; } | |
1264 | |||
1265 | /** | ||
1266 | * @brief Returns the flag to know if m_initMappingElem is allocated (constant version) | ||
1267 | */ | ||
1268 | inline bool isMappingElemAllocated() const { return m_initMappingElem; } | ||
1269 | |||
1270 | /** | ||
1271 | * @brief Returns the flag to know if m_initMappingElem is allocated | ||
1272 | */ | ||
1273 | 11 | inline bool& isMappingElemAllocated() { return m_initMappingElem; } | |
1274 | |||
1275 | /** | ||
1276 | * @brief Returns the local to global mapping of the elements (constant version) | ||
1277 | */ | ||
1278 | inline const ISLocalToGlobalMapping& mappingElem() const { return m_mappingElem[m_currentMesh]; } | ||
1279 | |||
1280 | /** | ||
1281 | * @brief Returns the local to global mapping of the elements | ||
1282 | */ | ||
1283 | 11 | inline ISLocalToGlobalMapping& mappingElem() { return m_mappingElem[m_currentMesh]; } | |
1284 | |||
1285 | /** | ||
1286 | * @brief Returns the local to global mapping of the elements (constant version) | ||
1287 | */ | ||
1288 | inline const ISLocalToGlobalMapping& mappingElem(int id) const { return m_mappingElem[id]; } | ||
1289 | |||
1290 | /** | ||
1291 | * @brief Returns the local to global mapping of the elements | ||
1292 | */ | ||
1293 | ✗ | inline ISLocalToGlobalMapping& mappingElem(int id) { return m_mappingElem[id]; } | |
1294 | |||
1295 | /** | ||
1296 | * @brief Returns the flag to know if m_initMappingElemSupport is allocated (constant version) | ||
1297 | */ | ||
1298 | inline bool isMappingElemSupportAllocated() const { return m_initMappingElemSupport; } | ||
1299 | |||
1300 | /** | ||
1301 | * @brief Returns the flag to know if m_initMappingElemSupport is allocated | ||
1302 | */ | ||
1303 | 11 | inline bool& isMappingElemSupportAllocated() { return m_initMappingElemSupport; } | |
1304 | |||
1305 | /** | ||
1306 | * @brief Returns the local to global mapping of the supportDof elements for the given unknown (constant version) | ||
1307 | * \param[in] iunknown The ID of the unknown | ||
1308 | */ | ||
1309 | inline const ISLocalToGlobalMapping* mappingElemSupportPerUnknown(int iUnknown) const { return m_mappingElemSupportPerUnknown[iUnknown]; } | ||
1310 | |||
1311 | /** | ||
1312 | * @brief Returns the local to global mapping of the supportDof elements for the given unknown | ||
1313 | * \param[in] iunknown The ID of the unknown | ||
1314 | */ | ||
1315 | 22 | inline ISLocalToGlobalMapping* mappingElemSupportPerUnknown(int iUnknown) { return m_mappingElemSupportPerUnknown[iUnknown]; } | |
1316 | |||
1317 | /** | ||
1318 | * @brief Returns the flag to know if m_initMappingLocalFelisceToGlobalPetsc is allocated (constant version) | ||
1319 | */ | ||
1320 | inline bool isMappingLocalFelisceToGlobalPetscAllocated() const { return m_initMappingLocalFelisceToGlobalPetsc; } | ||
1321 | |||
1322 | /** | ||
1323 | * @brief Returns the flag to know if m_initMappingLocalFelisceToGlobalPetsc is allocated | ||
1324 | */ | ||
1325 | 11 | inline bool& isMappingLocalFelisceToGlobalPetscAllocated() { return m_initMappingLocalFelisceToGlobalPetsc; } | |
1326 | |||
1327 | /** | ||
1328 | * @brief Returns the local id dof per unknown to the global id in PETSc ordering (constant version) | ||
1329 | * \param[in] iunknown The ID of the unknown | ||
1330 | */ | ||
1331 | inline const ISLocalToGlobalMapping& mappingIdFelisceToPetsc(int iUnknown) const { return m_mappingIdLocalDofPerUnknownToIdGlobalDofPetsc[iUnknown]; } | ||
1332 | |||
1333 | /** | ||
1334 | * @brief Returns the local id dof per unknown to the global id in PETSc ordering | ||
1335 | * \param[in] iunknown The ID of the unknown | ||
1336 | */ | ||
1337 | 22 | inline ISLocalToGlobalMapping& mappingIdFelisceToPetsc(int iUnknown) { return m_mappingIdLocalDofPerUnknownToIdGlobalDofPetsc[iUnknown]; } | |
1338 | |||
1339 | /** | ||
1340 | * @brief Returns the local id dof per unknown to the global id in PETSc ordering (constant version) | ||
1341 | * \param[in] iunknown The ID of the unknown | ||
1342 | */ | ||
1343 | inline const ISLocalToGlobalMapping& mappingDofLocalToDofGlobal(PhysicalVariable unknown) const { return m_mappingIdLocalDofPerUnknownToIdGlobalDofPetsc[m_listUnknown.getUnknownIdList(unknown)]; } | ||
1344 | |||
1345 | /** | ||
1346 | * @brief Returns the local id dof per unknown to the global id in PETSc ordering | ||
1347 | * \param[in] iunknown The ID of the unknown | ||
1348 | */ | ||
1349 | 1044 | inline ISLocalToGlobalMapping& mappingDofLocalToDofGlobal(PhysicalVariable unknown) { return m_mappingIdLocalDofPerUnknownToIdGlobalDofPetsc[m_listUnknown.getUnknownIdList(unknown)]; } | |
1350 | |||
1351 | /** | ||
1352 | * @brief Returns the flag to know if the system has been built (constant version) | ||
1353 | */ | ||
1354 | inline bool areSolutionAndRHSAllocated() const { return m_buildSystem; } // TODO change name | ||
1355 | |||
1356 | /** | ||
1357 | * @brief Returns the flag to know if the system has been built | ||
1358 | */ | ||
1359 | 11 | inline bool& areSolutionAndRHSAllocated() { return m_buildSystem; } // TODO change name | |
1360 | |||
1361 | /** | ||
1362 | * @brief Returns the number of matrices allocated | ||
1363 | */ | ||
1364 | 37218 | inline std::size_t numberOfMatrices() const { return m_matrices.size(); } | |
1365 | |||
1366 | /** | ||
1367 | * @brief Returns the vector of matrices | ||
1368 | */ | ||
1369 | inline std::vector<PetscMatrix> & matrices() { return m_matrices; } | ||
1370 | |||
1371 | /** | ||
1372 | * @brief Returns the i-th matrix (constant version) | ||
1373 | * \param[in] i Matrix index number | ||
1374 | */ | ||
1375 | inline PetscMatrix const& matrix(const std::size_t index = 0) const { return m_matrices[index]; } | ||
1376 | |||
1377 | /** | ||
1378 | * @brief Returns the i-th matrix | ||
1379 | * \param[in] i Matrix index number | ||
1380 | */ | ||
1381 | 103277 | inline PetscMatrix& matrix(const std::size_t index = 0) { return m_matrices[index]; } | |
1382 | |||
1383 | /** | ||
1384 | * @brief Returns the number of rhs allocated | ||
1385 | */ | ||
1386 | 22 | inline std::size_t numberOfVectors() const { return m_vectors.size(); } | |
1387 | |||
1388 | /** | ||
1389 | * @brief Returns the vector of rhs | ||
1390 | */ | ||
1391 | inline std::vector<PetscVector> & vectors() { return m_vectors; } | ||
1392 | |||
1393 | /** | ||
1394 | * @brief Returns the i-th rhs (constant version) | ||
1395 | * \param[in] i Rhs index number | ||
1396 | */ | ||
1397 | 41822 | inline const PetscVector& vector(const std::size_t index = 0) const { return m_vectors[index]; } | |
1398 | |||
1399 | /** | ||
1400 | * @brief Returns the i-th rhs | ||
1401 | * \param[in] i Rhs index number | ||
1402 | */ | ||
1403 | 1872353 | inline PetscVector& vector(const std::size_t index = 0) { return m_vectors[index]; } | |
1404 | |||
1405 | /** | ||
1406 | * @brief Returns the local solution (constant version) | ||
1407 | */ | ||
1408 | inline const PetscVector& solution() const { return m_sol; } | ||
1409 | |||
1410 | /** | ||
1411 | * @brief Returns the local solution | ||
1412 | */ | ||
1413 | 30272 | inline PetscVector& solution() { return m_sol; } | |
1414 | |||
1415 | /** | ||
1416 | * @brief Returns the flag to know if the sequential solution is allocated | ||
1417 | */ | ||
1418 | inline bool isSequentialSolutionAllocated() const { return m_seqSolAllocated; } | ||
1419 | |||
1420 | /** | ||
1421 | * @brief Returns the flag to know if the sequential solution is allocated | ||
1422 | */ | ||
1423 | 11 | inline bool& isSequentialSolutionAllocated() { return m_seqSolAllocated; } | |
1424 | |||
1425 | /** | ||
1426 | * @brief Returns the sequential solution (constant version) | ||
1427 | */ | ||
1428 | inline const PetscVector& sequentialSolution() const { return m_seqSol; } | ||
1429 | |||
1430 | /** | ||
1431 | * @brief Returns the sequential solution | ||
1432 | */ | ||
1433 | 14584400 | inline PetscVector & sequentialSolution() { return m_seqSol; } | |
1434 | |||
1435 | /** | ||
1436 | * @brief Returns the linear problem local mesh (constant version) | ||
1437 | */ | ||
1438 | inline std::size_t currentMesh() const { return m_currentMesh; } | ||
1439 | |||
1440 | /** | ||
1441 | * @brief Returns the linear problem local mesh | ||
1442 | */ | ||
1443 | 195480 | inline std::size_t& currentMesh() { return m_currentMesh; } | |
1444 | |||
1445 | /** | ||
1446 | * @brief Returns the list of boundary conditions (constant version) | ||
1447 | */ | ||
1448 | inline const BoundaryConditionList& getBoundaryConditionList() const { return m_boundaryConditionList; } | ||
1449 | |||
1450 | /** | ||
1451 | * @brief Returns the list of boundary conditions | ||
1452 | */ | ||
1453 | 1910 | inline BoundaryConditionList& getBoundaryConditionList() { return m_boundaryConditionList; } | |
1454 | |||
1455 | /** | ||
1456 | * @brief Returns the KSPInterface object | ||
1457 | */ | ||
1458 | 2220 | inline KSPInterface& kspInterface() { return *m_KSPInterface; } | |
1459 | |||
1460 | /** | ||
1461 | * @brief Returns the SNESInterface object | ||
1462 | */ | ||
1463 | 45149 | inline SNESInterface& snesInterface() {return *m_pSNESInterface;} | |
1464 | |||
1465 | /** | ||
1466 | * @brief Returns m_evaluationState (constant version) | ||
1467 | */ | ||
1468 | 20911 | inline const PetscVector& evaluationState() const { return m_evaluationState; } | |
1469 | |||
1470 | /** | ||
1471 | * @brief Returns m_evaluationState | ||
1472 | */ | ||
1473 | 341626 | inline PetscVector& evaluationState() { return m_evaluationState; } | |
1474 | |||
1475 | /** | ||
1476 | * @brief Returns m_seqEvaluationState (constant version) | ||
1477 | */ | ||
1478 | inline const PetscVector & seqEvaluationState() const { return m_seqEvaluationState; } | ||
1479 | |||
1480 | /** | ||
1481 | * @brief Returns m_seqEvaluationState | ||
1482 | */ | ||
1483 | 978769 | inline PetscVector & seqEvaluationState() { return m_seqEvaluationState; } | |
1484 | |||
1485 | /** | ||
1486 | * @brief Returns the linearized flag (constant version) | ||
1487 | */ | ||
1488 | inline bool linearizedFlag() const { return m_linearizedFlag; } | ||
1489 | |||
1490 | /** | ||
1491 | * @brief Returns the linearized flag | ||
1492 | */ | ||
1493 | 1968 | inline bool& linearizedFlag() { return m_linearizedFlag; } | |
1494 | |||
1495 | /** | ||
1496 | * @brief Returns the matrix without boundary condition | ||
1497 | */ | ||
1498 | ✗ | inline PetscMatrix& matrixWithoutBC() { return m_matrixWithoutBC; } | |
1499 | |||
1500 | /** | ||
1501 | * @brief Returns the rhs without boundary condition | ||
1502 | */ | ||
1503 | ✗ | inline PetscVector& rhsWithoutBC() { return m_rhsWithoutBC; } | |
1504 | |||
1505 | /** | ||
1506 | * @brief Returns the residual | ||
1507 | */ | ||
1508 | inline const PetscVector& residual() const { return m_residual; } | ||
1509 | |||
1510 | /** | ||
1511 | * @brief Returns the sequential residual | ||
1512 | */ | ||
1513 | ✗ | inline const PetscVector& seqResidual() const { return m_seqResidual; } | |
1514 | |||
1515 | /** | ||
1516 | * @brief Returns the static term flag (constant version) | ||
1517 | */ | ||
1518 | inline bool computeStaticTerm() const { return m_computeStaticTerm; } | ||
1519 | |||
1520 | /** | ||
1521 | * @brief Returns the static term flag | ||
1522 | */ | ||
1523 | 23604 | inline bool& computeStaticTerm() { return m_computeStaticTerm; } | |
1524 | |||
1525 | /** | ||
1526 | * @brief Returns the i-th extenal vector | ||
1527 | * \param[in] i External vector index | ||
1528 | */ | ||
1529 | 13010046 | inline PetscVector& externalVec(std::size_t i = 0) { return *m_externalVec[i]; } | |
1530 | |||
1531 | /** | ||
1532 | * @brief Returns the i-th extenal vector ordering | ||
1533 | * \param[in] i External vector index | ||
1534 | */ | ||
1535 | ✗ | inline AO& externalAO(std::size_t i = 0) { return m_externalAO[i]; } | |
1536 | |||
1537 | /** | ||
1538 | * @brief Push back external vector | ||
1539 | * \param[in] extVec The vector | ||
1540 | */ | ||
1541 |
1/2✓ Branch 1 taken 779 times.
✗ Branch 2 not taken.
|
779 | inline void pushBackExternalVec(PetscVector& extVec) { m_externalVec.push_back(&extVec); } |
1542 | |||
1543 | /** | ||
1544 | * @brief Push back external AO object | ||
1545 | * \param[in] extAO The vector | ||
1546 | */ | ||
1547 | 62 | inline void pushBackExternalAO(AO extAO) { m_externalAO.push_back(extAO); } | |
1548 | |||
1549 | /** | ||
1550 | * @brief Push back external Dof object | ||
1551 | * \param[in] extDof The vector | ||
1552 | */ | ||
1553 |
1/2✓ Branch 1 taken 62 times.
✗ Branch 2 not taken.
|
62 | inline void pushBackExternalDof(Dof& extDof) { m_externalDof.push_back(&extDof); } |
1554 | |||
1555 | /** | ||
1556 | * @brief Gather auxiliary vectors by name | ||
1557 | * \param[in] name The vector name | ||
1558 | */ | ||
1559 |
4/8✓ Branch 2 taken 18 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 18 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 18 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 18 times.
✗ Branch 12 not taken.
|
18 | inline void gather(std::string name) { gatherVec( m_vecs.Get(name), m_seqVecs.Get(name)); } |
1560 | |||
1561 | /*! \brief Get an auxiliary vector by name | ||
1562 | * \param[in] type Sequential or parallel | ||
1563 | * \param[in] name String defining the name of the vector | ||
1564 | * \return the petscVector | ||
1565 | */ | ||
1566 | 1104 | inline PetscVector& get(vectorType type, std::string name) | |
1567 | { | ||
1568 |
3/4✓ Branch 0 taken 1052 times.
✓ Branch 1 taken 52 times.
✓ Branch 3 taken 1052 times.
✗ Branch 4 not taken.
|
1104 | if ( type == sequential ) { if ( m_seqVecs.count(name) > 0 ) return m_seqVecs[name]; } |
1569 |
2/4✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 52 times.
✗ Branch 4 not taken.
|
52 | else if ( type == parallel ) { if ( m_vecs.count(name) > 0 ) return m_vecs[name]; } |
1570 | |||
1571 | ✗ | FEL_ERROR( "The std::vector " + name + " that you are looking for was not initialized." ); | |
1572 | ✗ | return m_vecs[name]; | |
1573 | } | ||
1574 | |||
1575 | /** | ||
1576 | * @brief CVGraphInterface (constant version) | ||
1577 | */ | ||
1578 | inline const CVGraphInterface::Pointer cvgraphInterf() const { return m_cvgraphInterf; } | ||
1579 | |||
1580 | /** | ||
1581 | * @brief CVGraphInterface | ||
1582 | */ | ||
1583 | ✗ | inline CVGraphInterface::Pointer cvgraphInterf() { return m_cvgraphInterf; } | |
1584 | |||
1585 | /** | ||
1586 | * @brief Returns the name of the linear problem | ||
1587 | */ | ||
1588 | ✗ | inline std::string name() const { return m_name; } | |
1589 | |||
1590 | /** | ||
1591 | * @brief Returns the pointer the to model owning this linear problem (constant version) | ||
1592 | */ | ||
1593 | inline Model* model() const { return mOwnerModel; } | ||
1594 | |||
1595 | /** | ||
1596 | * @brief Returns the pointer the to model owning this linear problem | ||
1597 | */ | ||
1598 | 42555 | inline Model*& model() { return mOwnerModel; } | |
1599 | |||
1600 | /** | ||
1601 | * @brief Returns the index if the parameter class attached to the model owning this linear problem (constant version) | ||
1602 | */ | ||
1603 | std::size_t instanceIndex() const; // { return mOwnerModel->instanceIndex(); } | ||
1604 | |||
1605 | /** | ||
1606 | * @brief Returns the index if the parameter class attached to the model owning this linear problem | ||
1607 | */ | ||
1608 | std::size_t& instanceIndex(); // { return mOwnerModel->instanceIndex(); } | ||
1609 | |||
1610 | /** | ||
1611 | * @brief Returns the ID of the linear problem(constant version) | ||
1612 | */ | ||
1613 | inline std::size_t identifierProblem() const { return m_identifier_problem; } | ||
1614 | |||
1615 | /** | ||
1616 | * @brief Returns the ID of the linear problem | ||
1617 | */ | ||
1618 | 513 | inline std::size_t& identifierProblem() { return m_identifier_problem; } | |
1619 | |||
1620 | ///@} | ||
1621 | ///@name Inquiry | ||
1622 | ///@{ | ||
1623 | |||
1624 | ///@} | ||
1625 | ///@name Input and output | ||
1626 | ///@{ | ||
1627 | |||
1628 | |||
1629 | //Write on files | ||
1630 | //============== | ||
1631 | |||
1632 | /// Print solution of the system solve. | ||
1633 | void printSolution(int verbose = 0, std::ostream& outstr = std::cout) const; // TODO move in a more suitable class | ||
1634 | |||
1635 | /// Write in file matrix.m, matrix in matlab format. | ||
1636 | void writeMatrixForMatlab(int verbose = 0, const std::string& filename = "matrix") const; // TODO move in a more suitable class | ||
1637 | |||
1638 | /// Write in file RHS.m, RHS in matlab format. | ||
1639 | void writeRHSForMatlab(int verbose = 0, const std::string& filename = "rhs") const; // TODO move in a more suitable class | ||
1640 | |||
1641 | // Write in file sol.m, solution in matlab format. | ||
1642 | void writeSolForMatlab(int verbose = 0, const std::string& filename = "solution") const; // TODO move in a more suitable class | ||
1643 | |||
1644 | /// Write in file filename, matrix in matlab format. | ||
1645 | void writeMatrixForMatlab(std::string const& fileName, PetscMatrix& matrix) const; // TODO move in a more suitable class | ||
1646 | |||
1647 | /// Write in file filenam, std::vector in matlab format. | ||
1648 | void writeVectorForMatlab(std::string const& fileName, PetscVector& vector) const; // TODO move in a more suitable class | ||
1649 | |||
1650 | /// Write in file matrix.mb, matrix in matlab format. | ||
1651 | void writeMatrixForMatlabBinary(int verbose = 0, const std::string& filename = "matrix") const; // TODO move in a more suitable class | ||
1652 | |||
1653 | /// Write in file RHS.mb, RHS in matlab format. | ||
1654 | void writeRHSForMatlabBinary(int verbose = 0, const std::string& filename = "rhs") const; // TODO move in a more suitable class | ||
1655 | |||
1656 | /// Write in file sol.mb, solution in matlab format. | ||
1657 | void writeSolForMatlabBinary(int verbose = 0, const std::string& filename = "solution") const; // TODO move in a more suitable class | ||
1658 | |||
1659 | /// Write in file filename, matrix in matlab format. | ||
1660 | void writeMatrixForMatlabBinary(std::string const& fileName, PetscMatrix& matrix) const; // TODO move in a more suitable class | ||
1661 | |||
1662 | /// Write in file filenam, std::vector in matlab format. | ||
1663 | void writeVectorForMatlabBinary(std::string const& fileName, PetscVector& vector) const; // TODO move in a more suitable class | ||
1664 | |||
1665 | /// Write solution from vector. | ||
1666 | void writeSolutionFromVec(PetscVector& v, int rank, std::vector<IO::Pointer>& io, double& time, int iteration, std::string prefix); // TODO move in a more suitable class | ||
1667 | |||
1668 | |||
1669 | //Write solution | ||
1670 | //============== | ||
1671 | |||
1672 | virtual void writeSolution(int rank, std::vector<IO::Pointer>& io, double& time, int iteration); | ||
1673 | |||
1674 | void writeSolutionBackup(const ListVariable& listVar, int rank, std::vector<IO::Pointer>& io, double& time, int iteration); | ||
1675 | |||
1676 | void fromVecToDoubleStar(double* solution, const PetscVector& sol, int rank, int dimension, felInt size=0); | ||
1677 | |||
1678 | void fromDoubleStarToVec(double* solution, PetscVector* sol, felInt size); | ||
1679 | |||
1680 | virtual void writeSupportVariableBackup(std::string name, double* supportVariable, int size, int dimension, int rank, IO::Pointer io, double& time, int iteration, bool manageMemory=false); | ||
1681 | |||
1682 | |||
1683 | //Write mesh to read solution with specific mesh (P2, Q2,...) | ||
1684 | //=========================================================== | ||
1685 | |||
1686 | void writeGeoForUnknown(int iUnknown, std::string nameOfGeoFile=""); | ||
1687 | |||
1688 | void writeGeoForUnknownEnsightGold(int iUnknown, std::map<felInt, std::vector<felInt> >& refToListOfIds, std::string nameOfGeoFile=""); | ||
1689 | |||
1690 | ///@} | ||
1691 | ///@name Operation that should be declared as friends function | ||
1692 | ///@{ | ||
1693 | |||
1694 | /// Remove average on solution. | ||
1695 | void removeAverageFromSolution(PhysicalVariable unknown, int rankProc); | ||
1696 | |||
1697 | /// Remove average on solution. Every processors do the same work at the moment (in progress). | ||
1698 | void removeAverageFromSolutionCurv(PhysicalVariable unknown, int rankProc); | ||
1699 | |||
1700 | void removeMaxSol(PhysicalVariable unknown, double valMax, int rankProc); | ||
1701 | |||
1702 | //Compute mean quantities or flux on a given label lists. | ||
1703 | // If PhysicalVariable has more than ncoor components -> compute Flux (typically for the velocity) | ||
1704 | // If PhysicalVariable is a scalar (num component = 1) -> compute Mean (typically for the Pressure) | ||
1705 | // --Version specific for the solution | ||
1706 | void computeMeanQuantity(PhysicalVariable unknown, const std::vector<int>& label, std::vector<double>& MeanQuantity); | ||
1707 | |||
1708 | //Compute mean quantities or flux on a given label lists. | ||
1709 | // If PhysicalVariable has more than ncoor components -> compute Flux (typically for the velocity) | ||
1710 | // If PhysicalVariable is a scalar (num component = 1) -> compute Mean (typically for the Pressure) | ||
1711 | // --Generic version for any petscVec | ||
1712 | void computeMeanQuantity(PetscVector& seqVec, PhysicalVariable unknown, const std::vector<int>& label, std::vector<double>& meanQuantity); | ||
1713 | |||
1714 | //Compute mean quantities or flux on a given label lists. | ||
1715 | // If PhysicalVariable has more than ncoor components -> compute Flux (typically for the velocity) | ||
1716 | // If PhysicalVariable is a scalar (num component = 1) -> compute Mean (typically for the Pressure) | ||
1717 | // --Generic version for any petscVec | ||
1718 | void computeMeanExternalQuantity(PetscVector& seqVec, AO externalAO, PhysicalVariable unknown, const std::vector<int>& label, std::vector<double>& meanQuantity); | ||
1719 | |||
1720 | //Compute mean quantities for each component on a given label lists. | ||
1721 | void computeMeanQuantityDomain(PhysicalVariable unknown, const std::vector<int>& label, std::vector<double>& MeanQuantity, std::vector<double>& measure); | ||
1722 | |||
1723 | //Compute the measure of a boundary part | ||
1724 | void computeMeasure(const std::vector<int>& label, std::vector<double>& measure); | ||
1725 | |||
1726 | //Compute the measure and the averaged normal of a boundary part | ||
1727 | void computeMeasureNormal(const std::vector<int>& label, std::vector<double>& measure, std::vector<double>& normal); | ||
1728 | |||
1729 | // Compute L2 error between present solution unknown (u or p) and the reference solution fctExactSolOnDof, where the latter | ||
1730 | // is analytically computed on dof or read from a file (in this case, please sure that numDof(unknown) = size of fctExactSolOnDof - Cuc 10/10/13) | ||
1731 | void errorL2(PhysicalVariable unknown, std::vector<double>& fctExactSolOnDof, double& resultL2Norm); | ||
1732 | |||
1733 | // Computing \int_\Omega u \dot u i.e. norm(u,2)^2 (unknown = u) - Gregory 25/09/12 | ||
1734 | void computel2NormSquared(PhysicalVariable unknown, int rankProc, double& l2NormSquared); | ||
1735 | |||
1736 | // Computing \int_\Omega \nabla u : \nabla u where u = unknown | ||
1737 | // It's semi-norm(u,H1)^2 - Gregory 25/09/12 | ||
1738 | void computeSemiNormH1Squared(PhysicalVariable unknown, int rankProc, double& semiNormH1Squared); | ||
1739 | |||
1740 | // Computing \int_\Gamma \rho*0.5 (u \dot u) (u \dot n) where u = unknown - Gregory 25/09/12 | ||
1741 | void computePowerBDu(PhysicalVariable unknown, int rankProc, const std::vector<int>& label, double density, std::vector<double>& PowerBDu); | ||
1742 | |||
1743 | // Computing \int_\Gamma (grad u \dot n) \cdot w, where u = unknown and w is given | ||
1744 | void computeConvBoundary(PhysicalVariable unknown, PetscVector& testFunction, int rankProc, const std::vector<int>& label, std::vector<double>& convBD); | ||
1745 | |||
1746 | /* Given two std::vector functions, u and w, compute \int_{Gamma} \grad_u \cdot n \cdot w */ | ||
1747 | double computeViscBoundary(int label, PetscVector& W); | ||
1748 | |||
1749 | // Computing \int_\Gamma p (u \dot n) - Gregory 25/09/12 | ||
1750 | // Typically unknow1 = p and unknow2 = u | ||
1751 | void computePowerBDpu(PhysicalVariable unknown1, PhysicalVariable unknown2, int rankProc, const std::vector<int>& label, std::vector<double>& PowerBDpu); | ||
1752 | |||
1753 | // Computing \int_\Omega u(n) \dot u(n-1) - Gregory 25/09/12 | ||
1754 | void computeProdUn_Un_1(PhysicalVariable unknown, int rankProc, Bdf& bdf, double& prodUn_Un_1); | ||
1755 | |||
1756 | // Computing \int_\Omega (u(n)-u(n-1))^2 - Gregory 25/09/12 | ||
1757 | void computeDiffSquaredUn_Un_1(PhysicalVariable unknown, int rankProc, Bdf& bdf, double& diffSquareUn_Un_1); | ||
1758 | |||
1759 | // Computing numerical Temam's trick : \int_\Omega div(u(n-1)) \dot u(n) v | ||
1760 | void computeTemamTerm(PhysicalVariable unknown, int rankProc, Bdf& bdf, double& temamTerm); | ||
1761 | |||
1762 | void computeBoundaryStress(int coeff, int label, std::vector<double>& stress); | ||
1763 | |||
1764 | void computeDomainIntegralOfDivergence(PhysicalVariable unknown, const std::vector<int>& label, std::vector<double>& integralOfDivergence); | ||
1765 | |||
1766 | ///@} | ||
1767 | protected: | ||
1768 | ///@name Protected static Member Variables | ||
1769 | ///@{ | ||
1770 | |||
1771 | ///@} | ||
1772 | ///@name Protected member Variables | ||
1773 | ///@{ | ||
1774 | |||
1775 | /// Spatial dimesion of the problem | ||
1776 | int m_dimesion = -1; | ||
1777 | |||
1778 | /// Verbosity level considered | ||
1779 | int m_verbosity = 0; | ||
1780 | |||
1781 | /// Object to manage time splitting | ||
1782 | FelisceTransient::Pointer m_fstransient = nullptr; | ||
1783 | |||
1784 | /// Object which contains variables of the linear system. | ||
1785 | ListVariable m_listVariable; | ||
1786 | |||
1787 | /// Object which contains unknown of the linear system. | ||
1788 | ListUnknown m_listUnknown; | ||
1789 | |||
1790 | ///@} | ||
1791 | ///@name Protected member Variables related to the meshes | ||
1792 | ///@{ | ||
1793 | |||
1794 | /// Set of meshes used by unknowns | ||
1795 | std::vector<std::size_t> m_meshUnknown; | ||
1796 | |||
1797 | /// Vector of input meshes defining the domains of work. | ||
1798 | std::vector<GeometricMeshRegion::Pointer> m_mesh; | ||
1799 | |||
1800 | /// Vector of the partitioned meshes | ||
1801 | std::vector<GeometricMeshRegion::Pointer> m_meshLocal; | ||
1802 | |||
1803 | ///@} | ||
1804 | ///@name Protected member Variables related to Dof, SupportDof and their partition | ||
1805 | ///@{ | ||
1806 | |||
1807 | /// Object which manage degrees of freedom | ||
1808 | Dof m_dof; | ||
1809 | |||
1810 | /// Array which contains degree of freedom repartition (size = numDof, value = rank of the processor which keep in memory the dof). | ||
1811 | std::vector<int> m_dofPart; // TODO if it is realted to the dof. why is it here and not in the dof class? | ||
1812 | |||
1813 | /// Array which contains element repartition (size = numDof, value = rank of the processor which keep in memory the element). //size = numElem.. | ||
1814 | std::vector<std::vector<int> > m_eltPart; // TODO if it is realted to the dof. why is it here and not in the dof class? | ||
1815 | |||
1816 | /// Number of dof in the problem | ||
1817 | felInt m_numDof = 0; | ||
1818 | |||
1819 | /// Number of dof contained by the current processor | ||
1820 | felInt m_numDofLocal = 0; | ||
1821 | |||
1822 | /// Number of dof for unknown | ||
1823 | std::vector<felInt> m_numDofUnknown; | ||
1824 | |||
1825 | /// Number of dof local for unknown | ||
1826 | std::vector<felInt> m_numDofLocalUnknown; | ||
1827 | |||
1828 | /// Global information about support of the degrees of freedom | ||
1829 | std::vector<SupportDofMesh> m_supportDofUnknown; | ||
1830 | |||
1831 | /// Boolean to know if m_supportDofUnknownOriginal is allocated | ||
1832 | bool m_areOriginalSupportDofMeshAllocated = false; | ||
1833 | |||
1834 | /// Global information about the original support of the degrees of freedom | ||
1835 | std::vector<SupportDofMesh*> m_supportDofUnknownOriginal; | ||
1836 | |||
1837 | /// Local information about support of the degrees of freedom | ||
1838 | std::vector<SupportDofMesh> m_supportDofUnknownLocal; | ||
1839 | |||
1840 | /// Boolean to know if m_ao is allocated | ||
1841 | bool m_initAO = false; | ||
1842 | |||
1843 | /// Mapping between global problem numbering and matrix dof numbering (specific to Petsc utilisation). | ||
1844 | AO m_ao; // TODO if it is realted to the dof why is here and not in the dof class? | ||
1845 | |||
1846 | /// Boolean to know if m_mappingNodes is allocated | ||
1847 | bool m_initMappingNodes = false; | ||
1848 | |||
1849 | /// Mapping between local to global ordering of support dof. | ||
1850 | ISLocalToGlobalMapping m_mappingNodes; // TODO do we really need it? // TODO if it is realted to the dof. why is it here and not in the dof class? | ||
1851 | |||
1852 | /// Boolean to know if m_mappingElem is allocated | ||
1853 | bool m_initMappingElem = false; | ||
1854 | |||
1855 | /// Mapping between local to global ordering of element. | ||
1856 | std::vector<ISLocalToGlobalMapping> m_mappingElem; // TODO if it is realted to the mesh. why is it here and not in the mesh class? | ||
1857 | |||
1858 | /// Boolean to know if m_mappingElemSupportPerUnknown is allocated | ||
1859 | bool m_initMappingElemSupport = false; | ||
1860 | |||
1861 | /// Mapping between local to global ordering of element support | ||
1862 | std::vector<ISLocalToGlobalMapping*> m_mappingElemSupportPerUnknown; // TODO if it is realted to the supportDof. why is it here and not in the supportDof class? | ||
1863 | |||
1864 | /// Boolean to know if m_mappingIdLocalDofPerUnknownToIdGlobalDofPetsc is allocated | ||
1865 | bool m_initMappingLocalFelisceToGlobalPetsc = false; | ||
1866 | |||
1867 | /// Mapping between local to global ordering of element. | ||
1868 | std::vector<ISLocalToGlobalMapping> m_mappingIdLocalDofPerUnknownToIdGlobalDofPetsc; // TODO do we really need it? | ||
1869 | |||
1870 | ///@} | ||
1871 | ///@name Protected member Variables used for the linear system | ||
1872 | ///@{ | ||
1873 | |||
1874 | /// Boolean to know if the system (m_matrices, m_vectors...) is allocated | ||
1875 | bool m_buildSystem = false; | ||
1876 | |||
1877 | /// Matrices used in formulation: the first one is the matrix used to compute the solution | ||
1878 | std::vector<PetscMatrix> m_matrices; | ||
1879 | |||
1880 | /// Rhs used in formulation: the first one is the rhs used to compute the solution | ||
1881 | std::vector<PetscVector> m_vectors; | ||
1882 | |||
1883 | /// Solution of the problem | ||
1884 | PetscVector m_sol; | ||
1885 | |||
1886 | /// Boolean to know if the sequential solution is allocated | ||
1887 | bool m_seqSolAllocated = false; | ||
1888 | |||
1889 | /// Sequential solution of the problem | ||
1890 | PetscVector m_seqSol; | ||
1891 | |||
1892 | ///@} | ||
1893 | ///@name Protected member Variables for assembling the matrix | ||
1894 | ///@{ | ||
1895 | |||
1896 | /// ID of the current mesh used for assembling | ||
1897 | std::size_t m_currentMesh = 0; | ||
1898 | |||
1899 | /// Label of the current element | ||
1900 | int m_currentLabel; | ||
1901 | |||
1902 | /// List of finite element in the problem. | ||
1903 | ListCurrentFiniteElement m_listCurrentFiniteElement; | ||
1904 | |||
1905 | /// List of curvilinear finite element in the problem. | ||
1906 | ListCurvilinearFiniteElement m_listCurvilinearFiniteElement; | ||
1907 | |||
1908 | /// List of finite element with Bd in the problem. | ||
1909 | ListCurrentFiniteElementWithBd m_listCurrentFiniteElementWithBd; | ||
1910 | |||
1911 | /// Element matrix | ||
1912 | std::vector<ElementMatrix::Pointer> m_elementMat; | ||
1913 | |||
1914 | /// Element matrix transpose | ||
1915 | std::vector<ElementMatrix::Pointer> m_elementMatT; | ||
1916 | |||
1917 | /// Element vector | ||
1918 | std::vector<ElementVector::Pointer> m_elementVector; | ||
1919 | |||
1920 | /// Element matrix boundary element | ||
1921 | std::vector<ElementMatrix::Pointer> m_elementMatBD; | ||
1922 | |||
1923 | /// Element matrix boundary element transpose | ||
1924 | std::vector<ElementMatrix::Pointer> m_elementMatTBD; | ||
1925 | |||
1926 | /// Element vector boundary element | ||
1927 | std::vector<ElementVector::Pointer> m_elementVectorBD; | ||
1928 | |||
1929 | /// Temporary array: mapping from local to global position of matrix/vector lines | ||
1930 | felInt* m_GposLine = nullptr; | ||
1931 | |||
1932 | /// Temporary array: mapping from local to global position of matrix columns | ||
1933 | felInt* m_GposColumn = nullptr; | ||
1934 | |||
1935 | ///@} | ||
1936 | ///@name Protected member Variables for assembling the boundary condition | ||
1937 | ///@{ | ||
1938 | |||
1939 | /* | ||
1940 | what if we use | ||
1941 | std::unordered_map<typeOfBC,std::vector<ElementField> > m_mapElemFieldBC | ||
1942 | such that m_mapElemFieldBC[ neumann ] ~ m_elemFieldNeumann? | ||
1943 | in this way the code would be more readable and we can loop over them and generalize functions that, | ||
1944 | for now are huge cut and paste. | ||
1945 | */ | ||
1946 | |||
1947 | /// Object which contains all objects boundaryCondition to define boundaries conditions of the problem. | ||
1948 | BoundaryConditionList m_boundaryConditionList; | ||
1949 | |||
1950 | /// Vector of elemField to apply Neumann boundary condition. | ||
1951 | std::vector<ElementField> m_elemFieldNeumann; | ||
1952 | |||
1953 | /// Vector of elemField to apply Neumann Normal boundary condition. | ||
1954 | std::vector<ElementField> m_elemFieldNeumannNormal; | ||
1955 | |||
1956 | /// Vector of elemField to apply Robin boundary condition. | ||
1957 | std::vector<ElementField> m_elemFieldRobin; | ||
1958 | |||
1959 | /// Vector of elemField to apply Robin Normal boundary condition. | ||
1960 | std::vector<ElementField> m_elemFieldRobinNormal; | ||
1961 | |||
1962 | /// Vector of elemField to apply EmbedFSI boundary condition. | ||
1963 | std::vector<ElementField> m_elemFieldEmbedFSI; | ||
1964 | |||
1965 | /// Vector of elemFieldNormal to apply EmbedFSI boundary condition. | ||
1966 | std::vector<ElementField> m_elemFieldNormalForEmbedFSI; | ||
1967 | |||
1968 | /// Vector of elemFieldAlpha to apply EmbedFSI boundary condition. | ||
1969 | std::vector<ElementField> m_elemFieldAlphaForEmbedFSI; | ||
1970 | |||
1971 | /// Vector of elemField to apply backflow stabilization boundary condition. | ||
1972 | std::vector<ElementField> m_elemFieldBackflowStab; | ||
1973 | |||
1974 | /// Boolean to know if m_idDofRings and m_idDofRingsSeq have been already computed | ||
1975 | bool m_ringsComputed = false; | ||
1976 | |||
1977 | /// Id of the dofs belonging to the rings | ||
1978 | std::unordered_set<int> m_idDofRings; | ||
1979 | |||
1980 | // Id of the dofs belonging to the rings seq vector | ||
1981 | std::unordered_set<int> m_idDofRingsSeq; | ||
1982 | |||
1983 | ///@} | ||
1984 | ///@name Protected member Variables for solving the system of equation | ||
1985 | ///@{ | ||
1986 | |||
1987 | /// Interface to PETSc linear solver (direct or iterative) | ||
1988 | KSPInterface::Pointer m_KSPInterface = felisce::make_shared<KSPInterface>(); | ||
1989 | |||
1990 | /// Interface to PETSc non linear solver (direct or iterative) | ||
1991 | SNESInterface::Pointer m_pSNESInterface = felisce::make_shared<SNESInterface>(); | ||
1992 | |||
1993 | /// Id of the solver if multiple solvers (default value = 0). | ||
1994 | int m_identifier_solver = 0; | ||
1995 | |||
1996 | /// Evaluation state of the residual of the problem (only useful in SNES method) | ||
1997 | PetscVector m_evaluationState; | ||
1998 | |||
1999 | /// Sequential evaluation state of the residual of the problem (only useful in SNES method) | ||
2000 | PetscVector m_seqEvaluationState; | ||
2001 | |||
2002 | /// Context for SNES solver | ||
2003 | SnesContext m_snesContext; | ||
2004 | |||
2005 | /// The flag that determines if the problem is linearized (used for FSI) | ||
2006 | bool m_linearizedFlag = false; | ||
2007 | |||
2008 | /// Matrix without boundary condition | ||
2009 | PetscMatrix m_matrixWithoutBC; | ||
2010 | |||
2011 | /// RHS without boundary condition | ||
2012 | PetscVector m_rhsWithoutBC; | ||
2013 | |||
2014 | /// Parallel residual | ||
2015 | PetscVector m_residual; | ||
2016 | |||
2017 | /// Sequential residual | ||
2018 | PetscVector m_seqResidual; | ||
2019 | |||
2020 | /// The flag that determines whether compute or not the static term | ||
2021 | bool m_computeStaticTerm = false; | ||
2022 | |||
2023 | ///@} | ||
2024 | ///@name Protected member Variables that point to external vectors (e.g. vectors defined in another linearProblem) | ||
2025 | ///@{ | ||
2026 | |||
2027 | /// Pointer to external vectors | ||
2028 | std::vector<PetscVector*> m_externalVec; | ||
2029 | |||
2030 | /// Pointer to external vectors ordering | ||
2031 | std::vector<AO> m_externalAO; | ||
2032 | |||
2033 | /// Pointer to external vectors of dof objects | ||
2034 | std::vector<Dof*> m_externalDof; | ||
2035 | |||
2036 | ///@} | ||
2037 | ///@name Protected member Variables auxiliary maps of petscVector | ||
2038 | ///@{ | ||
2039 | |||
2040 | /// Map of petscVector | ||
2041 | Tools::mapHolder<std::string,PetscVector> m_vecs; | ||
2042 | |||
2043 | /// Map of sequential petscVector | ||
2044 | Tools::mapHolder<std::string,PetscVector> m_seqVecs; | ||
2045 | |||
2046 | ///@} | ||
2047 | ///@name Protected member Variables CVGraphInterface | ||
2048 | ///@{ | ||
2049 | |||
2050 | /// CVGraphInterface | ||
2051 | CVGraphInterface::Pointer m_cvgraphInterf; | ||
2052 | |||
2053 | ///@} | ||
2054 | ///@name Protected Operators | ||
2055 | ///@{ | ||
2056 | |||
2057 | ///@} | ||
2058 | ///@name Protected Operations | ||
2059 | ///@{ | ||
2060 | |||
2061 | /** | ||
2062 | * \brief Set the number of matrixes | ||
2063 | */ | ||
2064 | 62 | inline void setNumberOfMatrix(const std::size_t number) { m_matrices.resize(number); } | |
2065 | |||
2066 | /** | ||
2067 | * \brief Set the number of vectors. | ||
2068 | */ | ||
2069 | 16 | inline void setNumberOfVector(const std::size_t number) { m_vectors.resize(number); } | |
2070 | |||
2071 | ///@} | ||
2072 | ///@name Protected Access | ||
2073 | ///@{ | ||
2074 | |||
2075 | ///@} | ||
2076 | ///@name Protected Inquiry | ||
2077 | ///@{ | ||
2078 | |||
2079 | ///@} | ||
2080 | ///@name Protected LifeCycle | ||
2081 | ///@{ | ||
2082 | |||
2083 | ///@} | ||
2084 | private: | ||
2085 | ///@name Private static Member Variables | ||
2086 | ///@{ | ||
2087 | |||
2088 | ///@} | ||
2089 | ///@name Private member Variables | ||
2090 | ///@{ | ||
2091 | |||
2092 | /// Name of the current linear problem | ||
2093 | std::string m_name; | ||
2094 | |||
2095 | /// The Model owner of this linear problem | ||
2096 | Model* mOwnerModel; | ||
2097 | |||
2098 | /// Id of the problem if multiple problem (default value = 0). | ||
2099 | std::size_t m_identifier_problem = 0; | ||
2100 | |||
2101 | ///@} | ||
2102 | ///@name Private Operations | ||
2103 | ///@{ | ||
2104 | |||
2105 | /** | ||
2106 | * @brief Assembly loop when some support dof are duplicated | ||
2107 | * This function is called when some support dof are duplicated. It's only the loop over all the support | ||
2108 | * element. | ||
2109 | */ | ||
2110 | void m_duplicateSupportDofAssemblyLoop(int rank, ElementType& eltType, felInt ielGeo, std::vector<Point*>& elemPoint, std::vector<felInt>& elemIdPoint, FlagMatrixRHS flagMatrixRHS); | ||
2111 | |||
2112 | /** | ||
2113 | * @brief Assembly loop for boundary condition when some support dof are duplicated | ||
2114 | * This function is called when some support dof are duplicated. It's only the loop over all the support | ||
2115 | * element. | ||
2116 | */ | ||
2117 | void m_duplicateSupportDofAssemblyLoopBoundaryCondition(int rank, ElementType& eltType, felInt ielGeoByType, std::vector<Point*>& elemPoint, std::vector<felInt>& elemIdPoint, FlagMatrixRHS flagMatrixRHS); | ||
2118 | |||
2119 | /** | ||
2120 | * @brief Helper function for imposing essential boundary condition | ||
2121 | */ | ||
2122 | template<DirichletApplicationOptions::Values> | ||
2123 | void applyEssentialBoundaryConditionHelper(FlagMatrixRHS flagMatrixRHS); | ||
2124 | |||
2125 | /** | ||
2126 | * @brief Apply essential boundary condition via penalization | ||
2127 | */ | ||
2128 | void applyEssentialBoundaryConditionViaPenalizationOnBC(FlagMatrixRHS flagMatrixRHS, const BoundaryCondition* const BC); | ||
2129 | |||
2130 | /** | ||
2131 | * @brief Apply essential boundary condition via symmetric pseudo elimination | ||
2132 | */ | ||
2133 | void applyEssentialBoundaryConditionViaSymmetricPseudoEliminationOnBC(FlagMatrixRHS flagMatrixRHS, const BoundaryCondition* const BC, std::vector<felInt>& rows); | ||
2134 | |||
2135 | ///@} | ||
2136 | ///@name Private Access | ||
2137 | ///@{ | ||
2138 | |||
2139 | ///@} | ||
2140 | ///@name Private Inquiry | ||
2141 | ///@{ | ||
2142 | |||
2143 | ///@} | ||
2144 | #ifdef FELISCE_WITH_CVGRAPH | ||
2145 | // Many functions and variables used with CVGraph that should be moved into CVGraphInterface | ||
2146 | public: | ||
2147 | ///@name Operations | ||
2148 | ///@{ | ||
2149 | |||
2150 | void setValueMatrixBD(felInt ielSupportDof); | ||
2151 | void computeResidualOnDofs( PetscVector& residual, PetscVector& stressOut, PetscVector& seqStressOut, std::size_t iConn); | ||
2152 | void sendInitialCondition (std::size_t iConn); | ||
2153 | void applyOnlyCVGBoundaryCondition(); | ||
2154 | void assembleRHSOnlyCVGraph(); | ||
2155 | void prepareResidual(std::size_t iConn); | ||
2156 | void sumOnBoundaryCVGraph(PetscVector& v, PetscVector& b, const DofBoundary& dofBD); | ||
2157 | |||
2158 | virtual void cvgAddExternalResidualToRHS(); | ||
2159 | virtual void cvgraphFinalizeEssBC(); | ||
2160 | |||
2161 | ✗ | virtual void readData() {}; | |
2162 | ✗ | virtual void assembleMassBoundaryAndInitKSP( std::size_t /*iConn*/ ) { FEL_ERROR("You need to implement this function in this linear problem"); } | |
2163 | ✗ | virtual void cvgraphNaturalBC(felInt /*iel*/) {}; | |
2164 | |||
2165 | template<class theClass> | ||
2166 | void assembleMassBoundaryAndInitKSP( void (theClass::*functionOfTheLoop)( felInt ), | ||
2167 | const std::vector<felInt> & labels, | ||
2168 | void (theClass::*initPerElementType)(), | ||
2169 | void (theClass::*updateFunction)(const std::vector<Point*>&,const std::vector<int>&), | ||
2170 | std::size_t iBD); | ||
2171 | |||
2172 | template<class theClass> | ||
2173 | void assembleMassBoundaryOnly( void (theClass::*functionOfTheLoop)( felInt ), | ||
2174 | const std::vector<felInt> & labels, | ||
2175 | void (theClass::*initPerElementType)(), | ||
2176 | void (theClass::*updateFunction)(const std::vector<Point*>&,const std::vector<int>&), | ||
2177 | DofBoundary &dofBD, PetscMatrix&massMatrix); | ||
2178 | |||
2179 | ///@} | ||
2180 | ///@name Access | ||
2181 | ///@{ | ||
2182 | |||
2183 | ✗ | inline std::vector<DofBoundary>& dofBDVec() { return m_dofBD; } | |
2184 | ✗ | inline std::vector<PetscMatrix>& massBDVec() { return m_massBD; } | |
2185 | ✗ | inline std::vector<KSPInterface>& kspMassBDVec() { return m_kspMassBD; } | |
2186 | ✗ | inline DofBoundary& dofBD(std::size_t iBD) { return m_dofBD[iBD]; }; | |
2187 | ✗ | inline cvgMainSlave*& slave() { return m_slave; } | |
2188 | |||
2189 | ///@} | ||
2190 | protected: | ||
2191 | ///@name Protected member Variables | ||
2192 | ///@{ | ||
2193 | |||
2194 | std::vector<felInt> m_globPosColumn; | ||
2195 | std::vector<felInt> m_globPosRow; | ||
2196 | std::vector<double> m_matrixValues; | ||
2197 | DofBoundary* m_auxiliaryDofBD = nullptr; ///< an auxiliary pointer to a DofBoundary to be used in setValueMatrixBD() | ||
2198 | |||
2199 | std::vector<KSPInterface> m_kspMassBD; | ||
2200 | std::vector<PetscMatrix> m_massBD; | ||
2201 | PetscVector m_stressOut, m_seqStressOut; | ||
2202 | PetscMatrix m_auxiliaryMatrix; | ||
2203 | |||
2204 | std::string m_cvgDirichletVariable; | ||
2205 | |||
2206 | std::vector<DofBoundary> m_dofBD = std::vector<DofBoundary>(1, DofBoundary()); | ||
2207 | |||
2208 | ///@} | ||
2209 | private: | ||
2210 | ///@name Private member Variables | ||
2211 | ///@{ | ||
2212 | |||
2213 | cvgMainSlave* m_slave = nullptr; | ||
2214 | |||
2215 | ///@} | ||
2216 | ///@name Operations | ||
2217 | ///@{ | ||
2218 | |||
2219 | void sendZero(std::size_t iConn); | ||
2220 | |||
2221 | ///@} | ||
2222 | #endif | ||
2223 | }; | ||
2224 | ///@} | ||
2225 | ///@name Type Definitions | ||
2226 | ///@{ | ||
2227 | |||
2228 | ///@} | ||
2229 | } /* namespace felisce.*/ | ||
2230 | |||
2231 | |||
2232 | #include "linearProblem_template.hpp" | ||
2233 | |||
2234 | #endif // LINEAR_PROBLEM_HPP | ||
2235 |