GCC Code Coverage Report

Directory: ./
File: Solver/linearProblemLinearElasticity.hpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 3 7 42.9%
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:
13 //
14
15 #ifndef __LINEARPROBLEMLINEARELASTICITY_HPP__
16 #define __LINEARPROBLEMLINEARELASTICITY_HPP__
17
18 // System includes
19
20 // External includes
21
22 // Project includes
23 #include "Solver/linearProblem.hpp"
24 #include "FiniteElement/elementField.hpp"
25 #include "FiniteElement/elementVector.hpp"
26 #include "FiniteElement/elementMatrix.hpp"
27 #ifdef FELISCE_WITH_CVGRAPH
28 #include "Model/cvgMainSlave.hpp"
29 #endif
30
31 namespace felisce {
32
33 class LinearProblemLinearElasticity:
34 public LinearProblem {
35 public:
36 LinearProblemLinearElasticity();
37 8 ~LinearProblemLinearElasticity() override= default;;
38
39 void initialize(std::vector<GeometricMeshRegion::Pointer>& mesh, FelisceTransient::Pointer fstransient, MPI_Comm& comm, bool doUseSNES) override;
40 void initPerElementType(ElementType eltType, FlagMatrixRHS flagMatrixRHS = FlagMatrixRHS::matrix_and_rhs) override;
41 void initPerElementTypeBoundaryCondition(ElementType& eltType, FlagMatrixRHS flagMatrixRHS) override;
42 void computeElementArray(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint, felInt& iel, FlagMatrixRHS flagMatrixRHS = FlagMatrixRHS::matrix_and_rhs) override;
43 void computeElementArrayBoundaryCondition(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint, felInt& iel, FlagMatrixRHS flagMatrixRHS = FlagMatrixRHS::matrix_and_rhs) override;
44 124 void userFinalizeEssBCTransient() override{};
45 void updateCurrentVelocity();
46 void updateOldDisplacement();
47 void initializeDofBoundaryAndBD2VolMaps();
48 ✗ inline felInt iDisplacement() { return m_iDisplacement; };
49 void readStressInterface(std::map<int, std::map<int,std::map<int,int> > > externalMap, const PetscVector& externalSeqStress, felInt iUC0);
50 std::pair<double,double> computeTestQuantities();
51 void testQuantitiesComputer(felInt ielSupportDof);
52 ✗ void coupled(bool coupled) { m_coupled = coupled; };
53 double computeL2ScalarProduct(std::string l, std::string r);
54 void initPetscVecsForCoupling();
55 void applyEssentialBoundaryConditionDerivedProblem(int rank, FlagMatrixRHS flagMatrixRHS) override;
56 120 bool coupled() const {return m_coupled;};
57
58 /**
59 * @brief it computes the P1 normal field on the lateral surface
60 * This function compute the P1 outward normal starting from the discontinuos normal.
61 * 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
62 */
63 void computeNormalField(const std::vector<int> & labels, felInt idVecVariable);
64
65
66 /**
67 * @brief Used to compute the P1 outward normal
68 */
69 void normalFieldComputer(felInt ielSupportDof);
70
71 /**
72 * @brief Used to update fe used for the computation of the P1 outward normal
73 */
74 void updateFeNormVel(const std::vector<Point*>& elemPoint, const std::vector<int>& elemIdPoint);
75
76 /**
77 * @brief Used to normalized the P1 outward normal right after its computation
78 * It does a getValues from the sequential normalField and a setValue on the parallel normalField.
79 * It is then necessary to call an assembly on the parallel normalField after finished using it.
80 */
81 void normalizeComputer(felInt ielSupportDof);
82
83 /**
84 * @brief Virtual function to initialize element fields for the computation of the P1 outward normal
85 */
86 ✗ void initPerETNormVel() {};
87
88 ✗ inline DofBoundary& dofBD(std::size_t iBD) { return m_dofBD[iBD]; };
89
90 void sumOnBoundary(PetscVector& v, PetscVector& b, const DofBoundary& dofBD);
91
92 protected:
93 // Dof boundary utilities
94 std::vector<DofBoundary> m_dofBD = std::vector<DofBoundary>(1, DofBoundary());
95
96 protected:
97 CurrentFiniteElement* m_fe;
98 CurvilinearFiniteElement* m_curvFe;
99 felInt m_iDisplacement;
100
101 void initPerET();
102 void updateFe(const std::vector<Point*>& elemPoint,const std::vector<int>& /*elemIdPoint*/);
103 private:
104 void scalarProductComputer(felInt ielSupportDof);
105
106 public:
107 void computeBoundaryStress(const std::vector<int>& labels);
108 private:
109
110 void computeRHS(felInt iel);
111 //! Material properties
112 double m_mu;
113 double m_lambda;
114 double m_density;
115
116 double m_deltaT;
117 double m_gammaNM;
118 double m_betaNM;
119
120 std::vector<int> m_interfaceLabels;
121
122 bool m_coupled;
123 protected:
124 Tools::mapHolder<std::string,ElementField> m_elemFields;
125 //--------data for computeBoundaryStress--------
126 protected:
127 PetscVector m_computedStress, m_seqComputedStress;
128
129 std::vector<double> m_auxiliaryDoubles;
130 std::vector<std::string> m_auxiliaryStrings;
131 std::vector<int> m_auxiliaryInts;
132 private:
133 // Element fields (DOF_FIELD)
134 // std::vector, domain elements, displacement values
135 ElementField m_elemFieldDisplacementVolumeDofs;
136 // std::vector, boundary elements, stress values
137 ElementField m_elemFieldStress;
138 // temporary strcutures used in the loop
139 // matrix (d x d), grad(i,j) = d u^j / d x_i|_{current boundary dof}
140 UBlasMatrix m_grad;
141 // std::vector (d) stress = mu (grad n + grad^t n) +lambda tr(grad) n
142 UBlasVector m_stress;
143 ElementType m_volEltType,m_previousVolEltType;
144 void boundaryStressComputer(felInt ielSupportDof);
145 void initPerETBoundaryStress();
146 void updateFEBoundaryStress(const std::vector<Point*>& elemPoint, const std::vector<int>& elemIdPoint);
147
148 public:
149 void initFixedPointAcceleration();
150 /// This two methods implement aitken acceleration algorithm to
151 /// speed up fixed point iteration.
152 /// First you call accelerationPreStep
153 /// and you pass the input that you are about to use
154 /// Than you use it and you obtain a new input
155 /// you pass this input to accelerationPostStep
156 /// To obtain a smarter input based on that.
157 void accelerationPreStep(PetscVector& seqCurrentInput);
158 void accelerationPostStep(PetscVector& seqCurrentInput,felInt cIteration);
159 private:
160 /// Different available techniques
161 enum accelerationType {
162 FixedPoint = 0, /*!< Standard fixed point iterations. We need only x1 and only for computing norm of update */
163 Relaxation = 1, /*!< Relaxed fixed point iterations. We need only x1 */
164 Aitken = 2, /*!< Aitken acceleration. We need only x0,x1,fx0,fx1 */
165 IronsTuck = 3 /*!< Variant of Aitken acceleration. We need only x0,x1,fx0,fx1 */
166 };
167
168 accelerationType m_accMethod; ///< The acceleration method in use
169 double m_omegaAcceleration; ///< Accelaration parameter, computed during acceleration or setted
170
171
172 #ifdef FELISCE_WITH_CVGRAPH
173 public:
174 void sendData();
175 void readData() override;
176 private:
177 void assembleMassBoundaryAndInitKSP( std::size_t iConn ) override;
178 void cvgraphNaturalBC(felInt iel) override;
179 protected:
180 ElementField m_robinAux;
181
182 void massMatrixComputer(felInt ielSupportDof);
183 void initPerETMass();
184 void updateFE(const std::vector<Point*>& elemPoint, const std::vector<int>&);
185 #endif
186 };
187 }
188
189 #endif
190