GCC Code Coverage Report

Directory: ./
File: Solver/linearProblemNS.hpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 7 12 58.3%
Branches: 0 0 -%
Line Branch Exec Source
1 // ______ ______ _ _ _____ ______
2 // | ____| ____| | (_)/ ____| | ____|
3 // | |__ | |__ | | _| (___ ___| |__
4 // | __| | __| | | | |\___ \ / __| __|
5 // | | | |____| |____| |____) | (__| |____
6 // |_| |______|______|_|_____/ \___|______|
7 // Finite Elements for Life Sciences and Engineering
8 //
9 // License: LGL2.1 License
10 // FELiScE default license: LICENSE in root folder
11 //
12 // Main authors: J. Foulon & J-F. Gerbeau & V. Martin
13 //
14
15 #ifndef _LinearProblemNS_HPP
16 #define _LinearProblemNS_HPP
17
18 // System includes
19 #include <vector>
20
21 // External includes
22
23 // Project includes
24 #include "Core/felisceParam.hpp"
25 #include "FiniteElement/elementField.hpp"
26 #include "Geometry/Tools/locator.hpp"
27 #include "Solver/bdf.hpp"
28 #include "Solver/cardiacCycle.hpp"
29 #include "Solver/linearProblem.hpp"
30 #include "Solver/lumpedModelBC.hpp"
31 #include "Solver/RISModel.hpp"
32
33 /*!
34 \file linearProblemNS.hpp
35 \authors J. Foulon & J-F. Gerbeau & V. Martin
36 \date 05/01/2011
37 \brief Monolithic algorithm for the Navier-Stokes equations
38 */
39
40 // #define A_MIN3(a,b,c) ( (a) < (b) ? ((a)<(c) ? 0 : 2) : ((b)<(c) ? 1 : 2) )
41 // #define A_MAX3(a,b,c) ( (a) > (b) ? ((a)>(c) ? 0 : 2) : ((b)>(c) ? 1 : 2) )
42
43 namespace felisce {
44 enum class FlagMeshUpdated {PreviousMesh = 0, CurrentMesh = 1};
45
46 /*!
47 \class LinearProblemNS
48 \authors J. Foulon & J-F. Gerbeau & V. Martin
49 \date 05/01/2011
50 \brief Manage specific functions for Navier-Stokes problem
51 */
52 class LinearProblemNS:
53 public LinearProblem {
54
55 public:
56 LinearProblemNS();
57 ~LinearProblemNS() override;
58
59 void initialize(std::vector<GeometricMeshRegion::Pointer>& mesh, FelisceTransient::Pointer fstransient, MPI_Comm& comm, bool doUseSNES) override;
60
61 54 void initializeTimeScheme(Bdf* bdf) override { m_bdf = bdf; }
62
63 void initializeLumpedModelBC(LumpedModelBC* lumpedModelBC);
64
65 void initializeRISModel(RISModel* risModel);
66
67 void initializeCardiacCycle(CardiacCycle* cardiacCycle);
68
69 void InitializeDerivedAttributes() override;
70
71 void initPerElementType(ElementType eltType, FlagMatrixRHS flagMatrixRHS=FlagMatrixRHS::matrix_and_rhs) override;
72
73 void initPerElementTypeBoundaryCondition(ElementType& eltType, FlagMatrixRHS flagMatrixRHS = FlagMatrixRHS::matrix_and_rhs) override;
74
75 void allocateElemFieldBoundaryConditionDerivedLinPb(const int idBCforLinCombMethod=0) override;
76
77 void initPerDomainBoundaryCondition(std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint,int label,
78 felInt numEltPerLabel, felInt* ielSupportDof, ElementType& eltType,
79 felInt numElement_eltType, FlagMatrixRHS flagMatrixRHS = FlagMatrixRHS::matrix_and_rhs) override;
80
81 void computeElementArrayBoundaryCondition(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint, felInt& iel, FlagMatrixRHS flagMatrixRHS) override;
82
83
84 void applyEssentialBoundaryConditionDerivedProblem(int rank, FlagMatrixRHS flagMatrixRHS) override;
85 void applyNaturalBoundaryConditionLumpedModelBC(std::vector<Point*>& elemPoint,const std::vector<felInt>& elemIdPoint,
86 int label, felInt numEltPerLabel,
87 felInt* ielSupportDof, std::size_t iLumpedModelBC, ElementType& eltType,
88 felInt numElement_eltType, FlagMatrixRHS flagMatrixRHS = FlagMatrixRHS::matrix_and_rhs);
89 void computeElementArray(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint, felInt& iel, FlagMatrixRHS flagMatrixRHS = FlagMatrixRHS::matrix_and_rhs) override;
90 //!Assembly loop on boundary element of a certain label in implicit LumpedModelBC
91 void assembleMatrixImplLumpedModelBC(int rank, const std::size_t iLabel) override;
92 // LumpedModelBC: NS model with implicit implementation
93 void userChangePattern(int numProc, int rankProc) override;
94
95 //! Functions used in method of characteristics
96 //==============================================
97
98 void computeElementArrayCharact(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint,
99 felInt& iel, ElementType& eltType, felInt& ielGeo, FlagMatrixRHS flagMatrixRHS = FlagMatrixRHS::matrix_and_rhs) override;
100 //similar to functions of std::set value by vectors in ElementField class, but for a field of type QUAD_POINT_FIELD
101 //obliged to put here because mesh and some functions of LinearProblem class are needed
102 void setValueCharact(PetscVector& v, ElementField& elemFieldCharact, CurrentFiniteElement& fe, felInt iel,
103 int idVariable, const ElementType& eltType, felInt ielGeo, const double dt);
104 //backtracking the characteristic line a length dt
105 int charactLine(PetscVector& v, int idVariable, int numComp, const double dt, int numStep, int& effStep,
106 Point& pt, Point& localCoord, ElementType& eltType, felInt& ielGeo, Point& vel);
107 //find the next point while backtracking the characteristic line a length step from a point by a variant of RK4 scheme
108 int nextPt(PetscVector& v, int idVariable, int numComp, const double step,
109 Point& pt, Point& localCoord, ElementType& eltType, felInt& ielGeo, Point& vel);
110 //backtracking the characteristic line in one element
111 int charactLineElem(PetscVector& v, int idVariable, int numComp, double dtTol,
112 Point& pt, Point& localCoord, ElementType& eltType, felInt& ielGeo, Point& vel, double& dtLeft);
113 //std::vector interpolation
114 Point vecInterp(PetscVector& v, int idVariable, int numComp,
115 const Point& localCoord, ElementType& eltType, felInt ielGeo);
116 //===================================
117
118 void gatherVectorBeforeAssembleMatrixRHS() override;
119 void initExtrapol(PetscVector& V_1) override;
120 void updateExtrapol(PetscVector& V_1) override;
121 void addMatrixRHS() override;
122
123 void reviseSolution(PetscVector& alphai, std::vector<PetscVector>& stationnarySolutions);
124
125 // functions non linear solver
126 void evalDynamicResidual();
127 void addMatrixRHSNl(const FlagMatrixRHS flagMatrixRHS);
128 inline PetscVector& previousSolution() {
129 return m_solExtrapol;
130 }
131
132
133 /**
134 * @brief it computes the P1 normal field on the lateral surface
135 * This function compute the P1 outward normal starting from the discontinuos normal.
136 * It also computes, for each boundary dof, the measure of the all the boundary elements that share the dof it has to be called when you init the structure
137 */
138 void computeNormalField(const std::vector<int> & labels, felInt idVecVariable);
139
140
141 /**
142 * @brief Used to compute the P1 outward normal
143 */
144 void normalFieldComputer(felInt ielSupportDof);
145
146 /**
147 * @brief Used to update fe used for the computation of the P1 outward normal
148 */
149 void updateFeNormVel(const std::vector<Point*>& elemPoint, const std::vector<int>& elemIdPoint);
150
151 /**
152 * @brief Used to normalized the P1 outward normal right after its computation
153 * It does a getValues from the sequential normalField and a setValue on the parallel normalField.
154 * It is then necessary to call an assembly on the parallel normalField after finished using it.
155 */
156 void normalizeComputer(felInt ielSupportDof);
157
158 /**
159 * @brief Virtual function to initialize element fields for the computation of the P1 outward normal
160 */
161 24 void initPerETNormVel() {};
162
163
164 // Access functions
165 inline const LumpedModelBC* lumpedModelBC() const {
166 return m_lumpedModelBC;
167 }
168
169 ✗ virtual void updateVolume() {};
170
171 // assemble the terms for the face-oriented stabilization
172 void assembleFaceOrientedStabilization();
173
174 void computeBoundaryStress(PetscVector& vec);
175
176 320 inline felInt idVarVel() const { return m_iVelocity; }
177 inline felInt idVarPre() const { return m_iPressure; }
178 ✗ inline felInt idVarLag() const { return m_iLagrange; }
179 inline felInt idVarJap() const { return m_iJapanCon; }
180 16 inline felInt iUnknownVel() const { return m_iUnknownVel; }
181 308 inline felInt iUnknownPre() const { return m_iUnknownPre; }
182 ✗ inline felInt iUnknownLag() const { return m_iUnknownLag; }
183 inline felInt iUnknownJap() const { return m_iUnknownJap; }
184
185 bool m_hasLocalizationChanged = false;
186 inline bool hasLocalizationChanged() { return m_hasLocalizationChanged; }
187 void localizeItfQpInBackMesh(std::size_t itfMsh, std::size_t bckMsh, felInt itfVar, std::vector<std::vector<felInt>>& mapItfQpInBackMesh);
188 void evaluateFluidStructurePattern();
189
190 void computeSolidFluidMaps();
191 void updateStructureVel(CurBaseFiniteElement& fe, std::vector<felInt>& elemIdPoint);
192
193 void updateNormal(CurvilinearFiniteElement& fe, int idMsh, std::vector<felInt>& elemIdPoint);
194
195 void computeBHStabLagLag(const double coeff, const int iq);
196 void computeBHStabLagJap(const double coeff, const int iq);
197 void computeBHStabJapJap(const double coeff, const int iq);
198
199 void assembleBarbosaHughesStab();
200 void assembleJapanConstraintInterface();
201 void assembleJapanConstraintBoundary();
202 void assembleBrezziPitkarantaStab();
203
204 void computeLagrangeMultiplierOnInterface(FlagMatrixRHS flagMatrixRHS=FlagMatrixRHS::matrix_and_rhs);
205
206 void setInterfaceExchangedData(std::vector<double>& itfVel, std::vector<double>& itfForc);
207
208 void computeInterfaceStress();
209 void computeElementaryInterfaceStress(const std::vector<Point*>& elemPoint, std::vector<felInt>& elemIdPoint, felInt iel, ElementVector& elemVec);
210
211 void computeItfDisp(std::vector<double>& itfDisp);
212 void moveItfMeshes();
213
214 protected:
215 CurrentFiniteElement* m_feVel;
216 CurrentFiniteElement* m_fePres;
217 CurvilinearFiniteElement* m_curvFeVel;
218 CurvilinearFiniteElement* m_curvFePre;
219 CurvilinearFiniteElement* m_curvFeLag;
220 CurvilinearFiniteElement* m_curvFeStruc;
221 CurrentFiniteElementWithBd* m_feVelWithBd;
222 CurrentFiniteElementWithBd* m_fePreWithBd;
223 Variable* m_velocity;
224 Variable* m_pressure;
225 Variable* m_displacement;
226 felInt m_iVelocity; //variables
227 felInt m_iPressure;
228 felInt m_iLagrange;
229 felInt m_iJapanCon;
230 felInt m_iDisplacement;
231 felInt m_iUnknownVel; //unknowns
232 felInt m_iUnknownPre;
233 felInt m_iUnknownLag;
234 felInt m_iUnknownJap;
235
236 int m_velBlock;
237 int m_preBlock;
238 int m_lagBlock;
239 int m_japBlock;
240
241 std::vector<std::vector<felInt> > m_mapLagQpInVelMesh;
242 std::vector<std::vector<felInt> > m_mapLagQpInVelMeshSeq;
243 std::vector<std::vector<felInt> > m_mapJapQpInVelMesh;
244 std::vector<std::vector<felInt> > m_mapJapQpInVelMeshSeq;
245
246 std::vector<double> m_realDisp;
247 std::vector<double> m_fictDisp;
248
249 bool m_initializeMapPnt;
250 std::vector<felInt> m_mapPntRealToFict;
251
252 double m_viscosity;
253 double m_density;
254 double m_theta;
255 bool m_useSymmetricStress;
256 PetscVector m_seqBdfRHS;
257 PetscVector m_solExtrapol;
258
259 std::vector<int> m_auxiliaryInts;
260
261 //============== LumpedModelBC============
262 LumpedModelBC* m_lumpedModelBC;
263 //! Vector of elemField to apply LumpedModelBC boundary conditions.
264 std::vector<ElementField> m_elemFieldLumpedModelBC;
265 std::vector<double> PlumpedModelBC;
266
267 ElementField m_elemFieldRHS;
268
269 //Reduced steklov
270 bool m_forceConvAndStabComputation;
271
272 FlagMeshUpdated m_meshUpdateState = FlagMeshUpdated::PreviousMesh;
273
274 ✗ virtual void userComputeItfDisp() {};
275 ✗ virtual void userMoveFictMesh() {};
276
277 private:
278 // private function for face oriented stabilization
279 void addNewFaceOrientedContributor(felInt size, felInt idElt, std::vector<bool>& vec);
280 void updateFaceOrientedFEWithBd(CurrentFiniteElementWithBd* fevel, CurrentFiniteElementWithBd* fepre, std::vector<felInt>& idOfFaces, felInt numPoints, felInt idElt, felInt& idSup);
281
282 //! boolean which indicate if boundary condition for pressure is apply
283 ElementField m_elemFieldAdv;
284 ElementField m_elemFieldRHSbdf;
285 ElementField m_elemFieldCharact;
286 ElementField m_elemFieldTotalPressureFormulation;
287 ElementField m_elemFieldSistole;
288 ElementField m_elemFieldDiastole;
289
290 //! Element field for Newton interation
291 ElementField m_elemFieldVelSeq;
292 ElementField m_elemFieldPresSeq;
293 PetscVector m_betaSeq;
294
295 Bdf* m_bdf;
296 CardiacCycle* m_cardiacCycle;
297 RISModel* m_ris;
298
299 bool allocateSeqVec;
300 bool allocateSeqVecExt;
301
302 // Characteristic Method
303 Locator::Pointer m_pLocator = nullptr;
304
305 //ALE
306 ElementField m_elemFieldVelMesh;
307 PetscVector m_beta;
308 int numDofExtensionProblem;
309
310 std::vector<PetscInt> m_petscToGlobal1;
311 std::vector<PetscInt> m_petscToGlobal2;
312 std::vector<PetscScalar> m_auxvec;
313
314 // Fictitious + lagrange multiplier
315 std::vector<double> *m_structVel = nullptr;
316 std::vector<double> *m_intfForc = nullptr;
317 ElementField m_structVelQuad;
318 std::vector<Point> m_ficInterface;
319
320 public:
321 ElementField m_structNormalQuad;
322
323 private:
324 ElementField m_seqVelDofPts;
325 ElementField m_seqLagDofPts;
326 CSRMatrixPattern m_csrNoInterface;
327 bool m_isCsrNoInterfaceStored;
328
329 //Theta-Method
330 std::vector<ElementMatrix*> m_elMatRHSThetaMethod;
331 void computeRHSThetaMethod(felInt& iel);
332 ElementField m_elemFieldPres;
333
334 public:
335
336 // distinction in order to add lbg term on current mesh with only_rhs.
337 1904 void onPreviousMesh(){ m_meshUpdateState = FlagMeshUpdated::PreviousMesh;}
338 360 void onCurrentMesh(){ m_meshUpdateState = FlagMeshUpdated::CurrentMesh;}
339
340
341 #ifdef FELISCE_WITH_CVGRAPH
342 public:
343 // Read/Send
344 void sendData();
345 void readData() override;
346
347 private:
348 void cvgraphNaturalBC(felInt iel) override;
349 ElementField m_robinAux;
350 // Functions for mass matrix boundary computation
351 void initPerETMass();
352 void updateFE(const std::vector<Point*>& elemPoint, const std::vector<int>&);
353 void massMatrixComputer(felInt ielSupportDof);
354 void assembleMassBoundaryAndInitKSP( std::size_t iConn ) override;
355 #endif
356 };
357 }
358
359 #endif
360