GCC Code Coverage Report

Directory: ./
File: Solver/linearProblemNitscheXFEM.hpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 0 3 0.0%
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: M. Fernandez & F.M. Gerosa
13 //
14
15 /*!
16 \file linearProblemNitscheXFEM.cpp
17 \authors M. Fernandez & F.M. Gerosa
18 \date 04/2018
19 \brief Solver for a fictitious domain and Nitsche-XFEM fluid formulation.
20 */
21
22 #ifndef _LinearProblemNitscheXFEM_HPP
23 #define _LinearProblemNitscheXFEM_HPP
24
25 // System includes
26 #include <vector>
27
28 // External includes
29
30 // Project includes
31 #include "Solver/linearProblem.hpp"
32 #include "Core/felisceParam.hpp"
33 #include "Core/felisceTransient.hpp"
34 #include "DegreeOfFreedom/duplicateSupportDof.hpp"
35 #include "FiniteElement/elementField.hpp"
36 #include "Solver/bdf.hpp"
37
38
39 namespace felisce {
40
41 // Class for the linear problem
42 class LinearProblemNitscheXFEM:
43 public LinearProblem
44 {
45
46 public:
47
48 // Constructor
49 LinearProblemNitscheXFEM();
50
51 // Destructor
52 ~LinearProblemNitscheXFEM();
53
54 // Pointer definition of LinearProblemNitscheXFEM
55 FELISCE_CLASS_POINTER_DEFINITION(LinearProblemNitscheXFEM);
56
57 // Initialize the linar problem (parameters, unknowns, ...)
58 void initialize(std::vector<GeometricMeshRegion::Pointer>& mesh, FelisceTransient::Pointer fstransient, MPI_Comm& comm, bool doUseSNES) override;
59
60 // Duplicate and initialize the support dof
61 void initSupportDofDerivedProblem() override;
62
63 // Initialize the bdf scheme (pointer to the one in the model)
64 ✗ void initializeTimeScheme(Bdf* bdf) override { m_bdf = bdf; }
65
66 // Get the current finite elements and initialize elemFields for the problem
67 void initPerElementType(ElementType eltType, FlagMatrixRHS flagMatrixRHS = FlagMatrixRHS::matrix_and_rhs) override;
68
69 // Get the current finite elements and initialize elemFields for Neumann BD
70 void initPerElementTypeBoundaryCondition(ElementType& eltType, FlagMatrixRHS flagMatrixRHS) override;
71
72 // Update matrix pattern
73 void userChangePattern(int numProc, int rankProc) override;
74
75 // Compute element array
76 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) override;
77
78 // Assemble front point tip elements
79 void assembleMatrixFrontPoints();
80
81 // Assemble DG terms on fictitious interface
82 void assembleMatrixTip();
83
84 // Assemble DG terms on tip elements faces
85 void assembleMatrixTipFaces();
86
87 // Assemble ghost penalty stabilisation
88 void assembleMatrixGhostPenalty();
89
90 // Assemble face oriented stabilisation
91 void assembleMatrixFaceOriented();
92
93 // Compute element array Neumann BD
94 void computeElementArrayBoundaryCondition(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint, felInt& ielSupportDof, FlagMatrixRHS flagMatrixRHS = FlagMatrixRHS::matrix_and_rhs) override;
95
96 // Compute and apply Neumann BD
97 void computeAndApplyElementNaturalBoundaryCondition(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint, felInt& ielSupportDof1, felInt& ielSupportDof2, felInt& ielGeo, FlagMatrixRHS flagMatrixRHS = FlagMatrixRHS::matrix_and_rhs) override;
98
99 // Assemble Darcy BD
100 void assembleMatrixBCDarcy();
101
102 // Get the pointer to the curvilinear finite element of the structure
103 void setStrucFiniteElement(CurvilinearFiniteElement* strucFE); // TODO D.C. remove this part once merged multimesh in master
104
105 // Add m_Matrix[1] to the current matrix (m_Matrix[0])
106 void addMatrixRHS() override;
107
108 // Get the reference to the object in the model duplicating the support dofs
109 void setDuplicateSupportObject(DuplicateSupportDof* object);
110
111 // Get a pointer to the structure linear problem
112 ✗ void setInterfaceMesh(GeometricMeshRegion* interfMesh) { m_interfMesh = interfMesh; } // TODO D.C. remove this part once merged multimesh in master
113
114 // Gather a vector to m_seqVelExtrapol
115 void gatherSeqVelExtrapol(std::vector<PetscVector>& parallelVec);
116
117 // Gather bdf->vector() to m_seqBdfRHS
118 void gatherSeqBdfRHS(std::vector<PetscVector>& parallelVec);
119
120 // Dynamics, copy original <-> current
121 void deleteDynamicData();
122
123 // Initialization objects previous time step
124 void initOldObjects();
125
126 // Update solution previous time step
127 void updateOldSolution(felInt countInitOld);
128
129 // Get solution previous time step
130 ✗ PetscVector const & solutionOld(felInt n) const { return m_solutionOld[n]; }
131
132 // Set the pointer to the structure velocity and stresses. (main_pvm)
133 void setInterfaceExchangedData(std::vector<double>& intfVelo, std::vector<double>& intfForc);
134
135 // Compute stresses on interface
136 void computeInterfaceStress();
137
138 // Print stresses
139 void printStructStress(int indexTime, double dt);
140
141 void printMeshPartition(int indexTime, double dt) const; // TODO remove
142
143 void writeDuplication(int rank, IO& io, double& time, int iteration) const;
144
145 void printSkipVolume(bool printSlipVolume);
146
147 protected:
148
149 CurrentFiniteElement* m_feVel; // volume velocity finite element
150 CurrentFiniteElement* m_fePres; // volume pressure finite element
151 CurrentFiniteElementWithBd* m_feVelWithBd; // volume finite element with edges
152 CurrentFiniteElementWithBd* m_fePreWithBd; // volume finite element with edges
153 CurvilinearFiniteElement *m_curvFeVel; // boundary velocity finite element
154 CurvilinearFiniteElement *m_curvFePress; // boundary pressure finite element
155 CurvilinearFiniteElement *m_curvFeDarcy; // finite element for pressure Darcy
156 CurvilinearFiniteElement *m_curvFeStruc; // finite element for the structure
157 Variable* m_velocity; // velocity
158 Variable* m_pressure; // pressure
159 Variable* m_presDarcy; // pressure darcy
160 felInt m_iVelocity; // id of the velocity
161 felInt m_iPressure; // id of the pressure
162 felInt m_iPreDarcy; // id of the pressure
163 felInt m_iUnknownVel; // id of the velocity unknown
164 felInt m_iUnknownPre; // id of the pressure unknown
165 felInt m_iUnknownDar; // id of the pressure darcy unknown
166 int m_velBlock; // index velocity block
167 int m_preBlock; // index pressure block
168 int m_darBlock; // index pressure darcy block
169 int m_useNSEquation; // flag to choose implementation method of advective term
170 double m_viscosity; // dynamic viscosity of the fluid
171 double m_density; // density of the fluid
172 double m_nitschePar; // Nitsche penalty parameter
173 bool m_useSymmetricStress; // expression of the Laplacian (epsilon(u) or grad(u))
174 bool m_useInterfaceStab; // interface flow stabilization flag
175
176 std::vector<PetscVector> m_seqVelExtrapol; // extrapolation of the velocity
177
178 // Duplicated support
179 DuplicateSupportDof *m_duplicateSupportElements; // manage the duplication of the support dofs
180
181 // Old supportDofUnknown
182 std::vector< std::vector<SupportDofMesh> > m_supportDofUnknownOld;
183
184 // Auxiliary arrays to store elements points
185 std::vector<Point*> m_mshPts;
186 std::vector<Point*> m_itfPts;
187 std::vector<Point*> m_subItfPts;
188 std::vector<Point> m_intPoints;
189
190 // Number of vertices per element
191 size_t m_numVerPerFluidElt;
192 size_t m_numVerPerSolidElt;
193
194 // Vectors to store previous time step info
195 std::vector<Dof> m_dofOld;
196 std::vector<AO> m_aoOld;
197
198 private:
199
200 // Compute elementary terms
201 void computeElementMatrixRHS();
202 void computeElementMatrixRHSPartDependingOnOldTimeStep(felInt iel, felInt ielOld, felInt idTimeStep);
203 void computeNitscheSigma_u_v(double velCoeff, double preCoeff);
204 void computeNitscheSigma_v_u(double velCoeff, double preCoeff);
205 void computeNitscheDpoint_Sigma_v(double velCoeff, double preCoeff);
206
207 // It contains all the Darcy-boundary integrals
208 void initPerElementTypeDarcy(ElementType& eltType, FlagMatrixRHS flagMatrixRHS);
209 void computeElementArrayDarcy(const std::vector<Point*>& elemPoint);
210 void computeAndApplyElementDarcy(const std::vector<Point*>& elemPoint, felInt ielSupportDof1, felInt ielSupportDof2, felInt ielGeoGlobal, FlagMatrixRHS flagMatrixRHS);
211 void elementComputeDarcyOnBoundary(const std::vector<Point*>& elemPoint);
212
213 // Update structure velocity
214 void updateStructureVel(felInt strucId);
215
216 // Update structure velocity on sub interface elements
217 void updateSubStructureVel(const std::vector<Point*>& ptElem, const std::vector<Point*>& ptSubElem);
218
219 // Compute elementary stress
220 void m_computeElementaryInterfaceStress(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint, const ElementType eltType, felInt iel, ElementVector& elemVec);
221
222 // Update quadrature point in substructure
223 void m_computeIntegrationPoints(std::vector<Point*>& ptElem); // TODO D.C. move to finite element class
224
225 // Evaluate structure velocity on sub element quadrature points
226 void m_computeFluidVelOnSubStructure();
227
228 // Evaluate the tensor at integration points of the sub structure
229 void m_computeFluidStress(double sideCoeff, const std::vector<double>& strucNormal);
230
231 // BDF time scheme
232 Bdf* m_bdf; // Bdf scheme, initialized in initializedTimeScheme
233 bool m_isSeqBdfRHSAllocated; // to know if m_seqBdfRHS has been allocated
234 std::vector<PetscVector> m_seqBdfRHS; // bdf part going to the rhs
235
236 // Sequential velocity
237 bool m_isSeqVelExtrapolAllocated; // to know if m_solExtrapol has been allocated
238
239 // Element fields
240 ElementField m_elemFieldAdv; // the extrapolated velocity
241 ElementField m_elemFieldAdvDup; // the extrapolated velocity on the duplicate elt
242 ElementField m_elemFieldAdvTmp;
243 ElementField m_elemFieldRHSbdf; // the bdf part that goes in the rhs
244 ElementField m_seqVelDofPts;
245 ElementField m_seqPreDofPts;
246 ElementField m_elemFieldNormal;
247 ElementField m_strucVelDofPts;
248 ElementField m_strucVelQuadPts;
249 ElementField m_ddotComp;
250 ElementField m_elemFieldSolidRHS;
251
252 // Interface mehs
253 GeometricMeshRegion* m_interfMesh; // structure linear problem
254
255 std::vector<PetscVector> m_solutionOld;
256 std::vector<PetscVector> m_sequentialSolutionOld;
257
258 // Communication with fluidToMaster
259 std::vector<double>* m_intfVelo;
260 std::vector<double>* m_intfForc;
261
262
263 double compAreaVolume(const std::vector<Point*>& elemPoint) const;
264 double m_totalVol = 0;
265 double m_skippedVol = 0;
266 double m_totalArea = 0;
267 double m_skippedArea = 0;
268 double m_tol = 0;
269 double m_eps = 1.e-8;
270 bool m_skipSmallVolume = true;
271 bool m_printSkipVolume = true;
272 };
273 }
274
275 #endif
276