GCC Code Coverage Report

Directory: ./
File: Solver/linearProblemNSSimplifiedFSI.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 521 708 73.6%
Branches: 640 1793 35.7%
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 // System includes
16
17 // External includes
18
19 // Project includes
20 #include "Solver/linearProblemNSSimplifiedFSI.hpp"
21 #include "FiniteElement/elementVector.hpp"
22 #include "FiniteElement/elementMatrix.hpp"
23 #include "Core/felisceTools.hpp"
24 #include "DegreeOfFreedom/boundaryCondition.hpp"
25
26 namespace felisce {
27 /** \brief Constructor of the LinearProblemNSSimplifiedFSI class.
28
29 By default m_autoregulated and m_iopCoupling are setted to false.
30 Then actual data are read from FelisceParam.
31 ___________________________________________________________________
32 Several _PetscVectors_ are created (empty and with no structure: just the pointer)
33 - _Parallel_:
34 + displacement,
35 + normalField,
36 + measStar,
37 + meanGradUN,
38 + meanCurvature,
39 + koiterCoefficients,
40 + normalVelocity,
41 + maxEigVec (if exportPrinciplaDirectionsAndCurvature is true),
42 + minEigVec (if exportPrinciplaDirectionsAndCurvature is true),
43 + partition (if exportDofPartition is true)
44 ___________________________________________________________________
45 - _Sequential_:
46 + displacement,
47 + displacementOld,
48 + normalField,
49 + measStar,
50 + meanGradUN,
51 + meanCurvature,
52 + koiterCoefficients,
53 + normalVelocity,
54 + normalVelocityOld,
55 + displacementCurrent,
56 + maxEigVec (if exportPrincipalDirectionsAndCurvature is true),
57 + minEigVec (if exportPrincipalDirectionsAndCurvature is true),
58 + partition (if exportDofPartition is true)
59 */
60 12 LinearProblemNSSimplifiedFSI::LinearProblemNSSimplifiedFSI():LinearProblemNS(),
61 12 m_autoregulated(false),
62
17/34
✓ Branch 7 taken 12 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 12 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 12 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 12 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 12 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 12 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 12 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 12 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 12 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 12 times.
✗ Branch 35 not taken.
✓ Branch 37 taken 12 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 12 times.
✗ Branch 41 not taken.
✓ Branch 43 taken 12 times.
✗ Branch 44 not taken.
✓ Branch 46 taken 12 times.
✗ Branch 47 not taken.
✓ Branch 49 taken 12 times.
✗ Branch 50 not taken.
✓ Branch 55 taken 12 times.
✗ Branch 56 not taken.
✓ Branch 58 taken 12 times.
✗ Branch 59 not taken.
 12 m_iopCoupling(false)
63 {
64
3/4
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
✓ Branch 4 taken 8 times.
 12 m_nonLinear=FelisceParam::instance().nonLinearKoiterModel?1:0;
65
2/4
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 12 times.
 12 if ( FelisceParam::instance().zeroDisplacementAtBoarders )
66 ✗ FelisceParam::instance().useEssDerivedBoundaryCondition = true;
67
3/4
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
✓ Branch 4 taken 8 times.
 12 if ( FelisceParam::instance().autoregulated )
68 4 this->autoregulated(true);
69
70
2/4
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
 12 m_vecs.Init("displacement");
71
2/4
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
 12 m_vecs.Init("normalField");
72
2/4
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
 12 m_vecs.Init("measStar");
73
2/4
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
 12 m_vecs.Init("meanGradUN");
74
2/4
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
 12 m_vecs.Init("meanCurvature");
75
2/4
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
 12 m_vecs.Init("koiterCoefficients");
76
2/4
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
 12 m_vecs.Init("normalVelocity");
77
2/4
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 12 times.
 12 if ( FelisceParam::instance().exportDofPartition )
78 ✗ m_vecs.Init("partition");
79
80
2/4
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
 12 m_seqVecs.Init("displacement");
81
2/4
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
 12 m_seqVecs.Init("displacementOld");
82
2/4
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
 12 m_seqVecs.Init("normalField");
83
2/4
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
 12 m_seqVecs.Init("measStar");
84
2/4
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
 12 m_seqVecs.Init("meanGradUN");
85
2/4
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
 12 m_seqVecs.Init("meanCurvature");
86
2/4
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
 12 m_seqVecs.Init("koiterCoefficients");
87
2/4
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
 12 m_seqVecs.Init("normalVelocity");
88
2/4
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
 12 m_seqVecs.Init("normalVelocityOld");
89
2/4
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
 12 m_seqVecs.Init("displacementCurrent");
90
91
2/4
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 12 times.
 12 if ( FelisceParam::instance().exportDofPartition )
92 ✗ m_seqVecs.Init("partition");
93
94
3/4
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 4 times.
 12 if ( FelisceParam::instance().exportPrincipalDirectionsAndCurvature ) {
95
2/4
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 8 times.
✗ Branch 6 not taken.
 8 m_vecs.Init("maxEigVec");
96
2/4
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 8 times.
✗ Branch 6 not taken.
 8 m_vecs.Init("minEigVec");
97
2/4
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 8 times.
✗ Branch 6 not taken.
 8 m_seqVecs.Init("maxEigVec");
98
2/4
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 8 times.
✗ Branch 6 not taken.
 8 m_seqVecs.Init("minEigVec");
99 }
100 12 }
101
102 /** \brief initialize LinearProblem::m_dofBD.
103
104 It is important to call this function when we want to use the utilities in DofBoundary.
105 This class allows you to define a mapping between the volume numbering and a new indexing on
106 the boundary. This also let you define PetscVector only on the boundary.
107
108 TODO: use these class to redefine all the vectors of this class only on the boundary: this will save a lot of memory
109 */
110 void
111 4 LinearProblemNSSimplifiedFSI::initializeDofBoundaryAndBD2VolMaps() {
112 /*! The object LinearProblem::m_dofBD is initialized. */
113 4 m_dofBD[0/*iBD*/].initialize(this);
114 /*! The list of dofs on the boundary, for each variable and for each component are computed. */
115
2/2
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 4 times.
 12 for ( felInt iComp(0); iComp < this->dimension(); ++iComp )
116
1/2
✓ Branch 5 taken 8 times.
✗ Branch 6 not taken.
 8 this->dofBD(/*iBD*/0).buildListOfBoundaryPetscDofs(this, FelisceParam::instance().labelNSMesh, this->iUnknownVel() , iComp);
117
1/2
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
 4 this->dofBD(/*iBD*/0).buildListOfBoundaryPetscDofs(this, FelisceParam::instance().labelNSMesh, this->iUnknownPre() , /*iComp*/0);
118 4 }
119
120 /*! \brief \f$\vec u \cdot \vec n \f$ are reset to zero
121
122 These quantities are used in IOPCouplModel, use this function
123 to std::set them to zero.
124 */
125 void
126 12 LinearProblemNSSimplifiedFSI::resetInterfaceQuantities() {
127
3/6
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 12 times.
✗ Branch 9 not taken.
 12 m_seqVecs.Get("normalVelocity").set(0.0);
128
3/6
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 12 times.
✗ Branch 9 not taken.
 12 m_seqVecs.Get("normalVelocityOld").set(0.0);
129 12 }
130 /** \brief It is a call to handleLumpedModel in #m_autoregulation (Autoregulation class)
131
132 @param t current time step
133
134 The parameters of the windkessel model will be updated accordingly to
135 what is prescribed in the autoregulation class.
136 */
137 void
138 44 LinearProblemNSSimplifiedFSI::handleWindkessel(const double& t) {
139 44 m_autoregulation.handleLumpedModel( t, this->lumpedModelBC());
140 44 }
141
142 /*! \brief wrapping of #m_autoregulation function Autoregulation::control()
143
144 \param[in] t current time step
145 \return the value of the active tension at the current time step
146 */
147 double
148 181621 LinearProblemNSSimplifiedFSI::control(const double& t) {
149 181621 return m_autoregulation.control(t);
150 }
151
152 /*!
153 \brief compute flow and mean incoming pressure for the current time step
154 \param[in] time the current time step
155 \param[in] rank the id of the current processor
156
157 This function is used for post processing
158 */
159 void
160 44 LinearProblemNSSimplifiedFSI::updateControl(const double& time, const int& rank) {
161
3/6
✓ Branch 2 taken 44 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 44 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 44 times.
✗ Branch 9 not taken.
 88 std::vector<int> labelInlet = m_label["inflow"];
162
3/6
✓ Branch 2 taken 44 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 44 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 44 times.
✗ Branch 9 not taken.
 88 std::vector<int> labelOutlet = m_label["outflow"];
163 44 double P(0),Q(0);
164
1/2
✓ Branch 3 taken 44 times.
✗ Branch 4 not taken.
 44 std::vector<double> pout(labelOutlet.size(),0.0);
165
1/2
✓ Branch 1 taken 44 times.
✗ Branch 2 not taken.
 44 if (labelInlet.size() > 0)
166 {
167
1/2
✓ Branch 3 taken 44 times.
✗ Branch 4 not taken.
 44 std::vector<double> quantity(labelInlet.size(), 0.0);
168
169 //! The flow rate is computed at the inlet
170
1/2
✓ Branch 1 taken 44 times.
✗ Branch 2 not taken.
 44 this->computeMeanQuantity(velocity, labelInlet, quantity);
171 44 Q = -quantity[0]; //TODO to be changed in case of several inputs
172
173 //! The mean incoming pressure is computed at the inlet
174
1/2
✓ Branch 1 taken 44 times.
✗ Branch 2 not taken.
 44 this->computeMeanQuantity(pressure, labelInlet, quantity);
175 44 P = quantity[0] / Autoregulation::m_cycleDuration;
176 44 } else {
177 ✗ std::cout<<__FILE__<<":"<<__LINE__<<" inflow(s) label(s) not defined"<<std::endl;
178 }
179
1/2
✓ Branch 1 taken 44 times.
✗ Branch 2 not taken.
 44 if (labelOutlet.size() > 0 ) {
180 //! The mean Pressure is computed in all the outlets
181
1/2
✓ Branch 1 taken 44 times.
✗ Branch 2 not taken.
 44 this->computeMeanQuantity(pressure, labelOutlet, pout);
182 } else {
183 ✗ std::cout<<__FILE__<<":"<<__LINE__<<" outlet(s) label(s) not defined"<<std::endl;
184 }
185
1/2
✓ Branch 2 taken 44 times.
✗ Branch 3 not taken.
 44 if (labelOutlet.size()*labelInlet.size() > 0 ) {
186 /*! The values are saved in #m_autoregulation, that stores
187 the values and update the mean values over the cardiac cycle
188 */
189
2/4
✓ Branch 1 taken 44 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 44 times.
✗ Branch 5 not taken.
 44 m_autoregulation.update(time,Q,P,pout,rank);
190 }
191 44 }
192
193 /*!
194 \brief wrapping of Autoregulation::exportSummary
195 \param[in] rank the id of the current processor
196 */
197 void
198 4 LinearProblemNSSimplifiedFSI::exportControl(const int& rank) {
199 4 m_autoregulation.exportSummary(rank);
200 4 }
201
202 /*!
203 \brief it initializes #m_autoregulation
204 \param[in] nc a number of cardiac cycles
205 */
206 void
207 8 LinearProblemNSSimplifiedFSI::initializeAutoregulation( const int& nc) {
208
209
1/2
✗ Branch 4 not taken.
✓ Branch 5 taken 8 times.
 8 if (m_mesh[m_currentMesh]->bagElementTypeDomainBoundary().size() != 1 ){
210 ✗ std::cout<<"m_mesh.bagElementTypeDomainBoundary().size() = "<<m_mesh[m_currentMesh]->bagElementTypeDomainBoundary().size()<<std::endl;
211 ✗ FEL_ERROR("bag element type domain boundary is either to big or empty, pay attention to call this function after the initializeLinearProblem of the model");
212 }
213 //! the initialize method of #m_autoregulation is called
214
1/2
✓ Branch 7 taken 8 times.
✗ Branch 8 not taken.
 8 m_label = m_autoregulation.initialize(GeometricMeshRegion::eltEnumToFelNameGeoEle[m_mesh[m_currentMesh]->bagElementTypeDomainBoundary()[0]].first);
215 /*! the difference between the pressure that will be used for this simulation and the reference one
216 is computed for as many cardiac cycles as specified in the input
217 */
218 8 m_autoregulation.computeExtraPressure(nc);
219 /*!
220 The value of the control function is exported for post-processing reason
221 */
222 8 m_autoregulation.exportControlFunctions(nc);
223 8 }
224
225 440 void LinearProblemNSSimplifiedFSI::initPerET() {
226 440 m_curvFeVel = m_listCurvilinearFiniteElement[m_iVelocity];
227 440 m_curvFePre = m_listCurvilinearFiniteElement[m_iPressure];
228 440 }
229 415853 void LinearProblemNSSimplifiedFSI::updateFe(const std::vector<Point*>& elemPoint,const std::vector<int>& /*elemIdPoint*/) {
230 415853 m_curvFeVel ->updateMeasNormal(0, elemPoint);
231 415853 m_curvFePre ->updateMeasNormal(0,elemPoint);
232 415853 }
233
234
235 8 void LinearProblemNSSimplifiedFSI::computeCurvatures() {
236 8 assemblyLoopBoundary(& LinearProblemNSSimplifiedFSI::curvaturesComputer);
237
3/6
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 8 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 8 times.
✗ Branch 9 not taken.
 8 m_vecs.Get("meanCurvature").assembly();
238
3/6
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 8 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 8 times.
✗ Branch 9 not taken.
 8 m_vecs.Get("koiterCoefficients").assembly();
239
7/14
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 8 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 8 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 8 times.
✗ Branch 13 not taken.
✓ Branch 16 taken 8 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 8 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 8 times.
✗ Branch 23 not taken.
 8 pointwiseDivide(m_vecs.Get("meanCurvature"), m_vecs.Get("meanCurvature"), m_vecs.Get("measStar"));
240
7/14
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 8 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 8 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 8 times.
✗ Branch 13 not taken.
✓ Branch 16 taken 8 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 8 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 8 times.
✗ Branch 23 not taken.
 8 pointwiseDivide(m_vecs.Get("koiterCoefficients"), m_vecs.Get("koiterCoefficients"), m_vecs.Get("measStar"));
241
5/10
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 8 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 8 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 8 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 8 times.
✗ Branch 16 not taken.
 8 gatherVector( m_vecs.Get("meanCurvature"), m_seqVecs.Get("meanCurvature") );
242
5/10
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 8 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 8 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 8 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 8 times.
✗ Branch 16 not taken.
 8 gatherVector( m_vecs.Get("koiterCoefficients"), m_seqVecs.Get("koiterCoefficients") );
243
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 if ( FelisceParam::instance().exportPrincipalDirectionsAndCurvature ) {
244
3/6
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 8 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 8 times.
✗ Branch 9 not taken.
 8 m_vecs.Get("maxEigVec").assembly();
245
3/6
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 8 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 8 times.
✗ Branch 9 not taken.
 8 m_vecs.Get("minEigVec").assembly();
246
7/14
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 8 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 8 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 8 times.
✗ Branch 13 not taken.
✓ Branch 16 taken 8 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 8 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 8 times.
✗ Branch 23 not taken.
 8 pointwiseDivide(m_vecs.Get("maxEigVec"), m_vecs.Get("maxEigVec"), m_vecs.Get("measStar"));
247
7/14
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 8 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 8 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 8 times.
✗ Branch 13 not taken.
✓ Branch 16 taken 8 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 8 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 8 times.
✗ Branch 23 not taken.
 8 pointwiseDivide(m_vecs.Get("minEigVec"), m_vecs.Get("minEigVec"), m_vecs.Get("measStar"));
248
5/10
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 8 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 8 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 8 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 8 times.
✗ Branch 16 not taken.
 8 gatherVector( m_vecs.Get("maxEigVec"), m_seqVecs.Get("maxEigVec") );
249
5/10
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 8 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 8 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 8 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 8 times.
✗ Branch 16 not taken.
 8 gatherVector( m_vecs.Get("minEigVec"), m_seqVecs.Get("minEigVec") );
250 }
251 8 }
252
253 //! \brief Initialize data members of the structure
254 12 void LinearProblemNSSimplifiedFSI::initStructure() {
255 /*!
256 When the iopCoupling is used, the vectors
257 for the coupling and for the corresponding iterations
258 are initialized.
259
260 Depending on the coupling technique:
261 - currentIop (sequential and parallel)
262 - x1 (only sequential)
263 - x0,fx0,fx1 (only sequential and only for Aitken and IronsTuck
264 */
265
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 8 times.
 12 if ( m_iopCoupling ) {
266
2/4
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
 4 m_vecs.Init("currentIop");
267
2/4
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
 4 m_seqVecs.Init("currentIop");
268
269 4 this->initFixedPointAcceleration();
270 }
271
272 12 this->initPetscVectors();
273
274
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 12 times.
 12 if ( FelisceParam::instance().exportDofPartition ) {
275 ✗ PetscVector aus;
276 ✗ std::vector<double> dofPartInDouble;
277 ✗ for (std::size_t idDofPart(0); idDofPart< m_dofPart.size(); idDofPart++)
278 ✗ dofPartInDouble.push_back(m_dofPart[idDofPart]);
279 ✗ this->fromDoubleStarToVec(&dofPartInDouble[0], &aus, m_dofPart.size());
280 ✗ m_vecs.Get("partition").assembly();
281 ✗ m_vecs.Get("partition").copyFrom(aus);
282 ✗ gatherVector(m_vecs.Get("partition"),m_seqVecs.Get("partition"));
283 }
284 //! Structure data and other parameters are read from FelisceParam
285 12 m_deltaT = m_fstransient->timeStep;
286 12 m_widthS = FelisceParam::instance().widthFiber + FelisceParam::instance().widthShell;
287 12 m_poissonS = FelisceParam::instance().poissonShell;
288 12 m_youngS = FelisceParam::instance().youngShell;
289 12 m_betaS = m_youngS* FelisceParam::instance().widthShell/(1.0 - m_poissonS*m_poissonS);
290 12 m_densityS = FelisceParam::instance().densityStructure;
291 12 m_viscS = FelisceParam::instance().viscosityStructure;
292 12 m_referencePressure = FelisceParam::instance().outsidePressure;
293 12 m_alphaS = m_widthS * m_densityS / m_deltaT;
294 12 m_fiberTensPassive = FelisceParam::instance().preStressFibers;
295 12 m_fiberYoung = FelisceParam::instance().youngFibers;
296 12 m_fiberDensity = FelisceParam::instance().fiberDensity;
297
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 12 times.
 12 if ( FelisceParam::verbose() > 1 ) {
298 ✗ std::stringstream out;
299 ✗ out << " fibers parameters--> ";
300 ✗ out << " tension: "<<m_fiberTensPassive;
301 ✗ out << " Young: "<<m_fiberYoung<<std::endl;
302 ✗ PetscPrintf(MpiInfo::petscComm(),"%s",out.str().c_str());
303 }
304
305
2/2
✓ Branch 1 taken 12 times.
✓ Branch 2 taken 12 times.
 24 for (std::size_t iEmbedFSI=0; iEmbedFSI < m_boundaryConditionList.numEmbedFSIBoundaryCondition(); iEmbedFSI++ ) {
306 12 BoundaryCondition* BC = m_boundaryConditionList.EmbedFSI(iEmbedFSI);
307
2/2
✓ Branch 6 taken 64 times.
✓ Branch 7 taken 12 times.
 76 for(auto it_labelNumber = BC->listLabel().begin(); it_labelNumber != BC->listLabel().end(); it_labelNumber++) {
308
1/2
✓ Branch 2 taken 64 times.
✗ Branch 3 not taken.
 64 m_embedLabels.push_back(*it_labelNumber);
309 }
310 }
311
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 if(m_fiberDensity.size() == 0){
312 // No fiberDensity was provided in the data file
313 // We create a fiberDensity equal to 1. on all the labels where there is a BC of EmbedFSI type
314 12 m_fiberDensity.resize(m_embedLabels.size());
315
2/2
✓ Branch 1 taken 64 times.
✓ Branch 2 taken 12 times.
 76 for(std::size_t i=0; i<m_fiberDensity.size();i++){
316 64 m_fiberDensity[i] = 1.;
317 }
318 }
319
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 12 times.
 12 if(m_fiberDensity.size() != m_embedLabels.size()){
320 ✗ std::cout << "m_fiberDensity.size() = " << m_fiberDensity.size() << std::endl;
321 ✗ std::cout << "m_embedLabels = " << m_embedLabels.size() << std::endl;
322 ✗ std::cout << "Error in " << __FILE__ << ": the number of values provided in the data file for fiberDensity should be the same as the number of EmbedFSI labels" << std::endl;
323 ✗ exit(1);
324 }
325
2/2
✓ Branch 1 taken 64 times.
✓ Branch 2 taken 12 times.
 76 for(std::size_t itFiberDens=0; itFiberDens<m_fiberDensity.size(); itFiberDens++){
326 64 m_label2embedFSI[m_embedLabels[itFiberDens]] = itFiberDens;
327 }
328
329 12 this->computeNormalField(m_embedLabels, m_iVelocity);
330 /*!
331 In a 3D problem the curvatures are computed by calling computeCurvatures() method.\n
332 In 2D the #m_radius is std::set from FelisceParam
333 */
334
2/4
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
 12 switch (this->dimension()) {
335 ✗ case 1:
336 ✗ FEL_ERROR("this structural model is not supposed to work on points")
337 ✗ break;
338 4 case 2:
339 4 m_radius = FelisceParam::instance().radius;
340 4 break;
341 8 case 3:
342 8 computeCurvatures();
343 8 break;
344 ✗ default:
345 ✗ FEL_ERROR("strange dimension parameter");
346 }
347 12 }
348
349 /*! \brief It updates the normalVelocity and displacementCurrent.
350 */
351 216 void LinearProblemNSSimplifiedFSI::updateNormalVelocityAndCurrentDisplacement() {
352 //! normalVelocityOld (seq) = normalVelocity (seq), normalVelocity is reset to zero.
353
5/10
✓ Branch 2 taken 216 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 216 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 216 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 216 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 216 times.
✗ Branch 16 not taken.
 216 m_seqVecs.Get("normalVelocityOld").copyFrom(m_seqVecs.Get("normalVelocity"));
354
3/6
✓ Branch 2 taken 216 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 216 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 216 times.
✗ Branch 9 not taken.
 216 m_vecs.Get("normalVelocity").set(0.0); //TODO remove!
355 //! normalVelocity (para) = solution(para) scalar normalField(para), normalVelocity is gathered
356
5/10
✓ Branch 2 taken 216 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 216 times.
✗ Branch 6 not taken.
✓ Branch 10 taken 216 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 216 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 216 times.
✗ Branch 17 not taken.
 216 pointwiseMult(m_vecs.Get("normalVelocity"),this->solution(), m_vecs.Get("normalField"));
357
5/10
✓ Branch 2 taken 216 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 216 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 216 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 216 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 216 times.
✗ Branch 16 not taken.
 216 gatherVector(m_vecs.Get("normalVelocity"), m_seqVecs.Get("normalVelocity"));
358 /*! The results of this operation is still in the form of a std::vector: its components are sum up and stored in the pressure block
359 by using computeScalarNormalVelocity()
360 */
361 216 computeScalarNormalVelocity();
362 //! displacementCurrent = displacement + deltaT * normalVelocity (all sequential)
363
5/10
✓ Branch 2 taken 216 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 216 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 216 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 216 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 216 times.
✗ Branch 16 not taken.
 216 m_seqVecs.Get("displacementCurrent").copyFrom(m_seqVecs.Get("displacement"));
364
5/10
✓ Branch 2 taken 216 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 216 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 216 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 216 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 216 times.
✗ Branch 16 not taken.
 216 m_seqVecs.Get("displacementCurrent").axpy(m_deltaT,m_seqVecs.Get("normalVelocity"));
365 216 }
366 /*! \brief It updates the displacement
367 TODO change the name of this function
368 */
369 64 void LinearProblemNSSimplifiedFSI::updateStructure() {
370
1/2
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
 64 gatherSolution();
371 //! A P1 approximation of gradient of the velocity in the normal direction is computed
372
1/2
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
 64 this->computeMeanGradUN();
373 //! displacementOld=displacement (sequential)
374
5/10
✓ Branch 2 taken 64 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 64 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 64 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 64 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 64 times.
✗ Branch 16 not taken.
 64 m_seqVecs.Get("displacementOld").copyFrom(m_seqVecs.Get("displacement"));
375 //! It updates the normalVelocity and displacementCurrent.
376
1/2
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
 64 this->updateNormalVelocityAndCurrentDisplacement();
377 /*!
378 The variable displacement is computed:\n
379 displacement is defined as vectorial variable, like the velocity, but it is a scalar defined as the sum of its (fake) components, this is done to avoid complications in parallel
380 so we first do the pointwise multiplication between the solution and the normal field to obatin the normal velocity then we multiply it by the time step and we add it to the displacement.\n
381 After that the function computeScalarDisplacement() is called to store the scalar displacement in the pressure block and also to compute its minimum and its maximum for post-process.
382 TODO: it is redundant, I think we should use displacementCurrent.
383 */
384
5/10
✓ Branch 2 taken 64 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 64 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 64 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 64 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 64 times.
✗ Branch 16 not taken.
 64 m_vecs.Get("displacement").axpy(m_deltaT, m_vecs.Get("normalVelocity"));
385
5/10
✓ Branch 2 taken 64 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 64 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 64 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 64 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 64 times.
✗ Branch 16 not taken.
 64 gatherVector( m_vecs.Get("displacement"), m_seqVecs.Get("displacement"));
386 64 double maxD(0),minD(0), gMaxD(0), gMinD(0);
387
1/2
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
 64 computeScalarDisplacement(maxD,minD);
388
2/4
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 64 times.
✗ Branch 5 not taken.
 64 MPI_Reduce(&maxD, &gMaxD, 1, MPI_DOUBLE, MPI_MAX, 0, MpiInfo::petscComm());
389
2/4
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 64 times.
✗ Branch 5 not taken.
 64 MPI_Reduce(&minD, &gMinD, 1, MPI_DOUBLE, MPI_MIN, 0, MpiInfo::petscComm());
390
2/4
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 64 times.
✗ Branch 5 not taken.
 64 std::stringstream s1,s2;
391
3/6
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 64 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 64 times.
✗ Branch 8 not taken.
 64 s1 << " Maximum displacement [micron]: "<< gMaxD * 1.e4 <<std::endl;
392
3/6
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 64 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 64 times.
✗ Branch 8 not taken.
 64 s2 << " Minimum displacement [micron]: "<< gMinD * 1.e4 <<std::endl;
393
3/6
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 64 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 64 times.
✗ Branch 9 not taken.
 64 PetscPrintf(MpiInfo::petscComm(), "%s",s1.str().c_str());
394
3/6
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 64 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 64 times.
✗ Branch 9 not taken.
 64 PetscPrintf(MpiInfo::petscComm(), "%s",s2.str().c_str());
395 64 }
396
397 204 void LinearProblemNSSimplifiedFSI::userElementInitNaturalBoundaryCondition() {
398
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 204 times.
 204 if ( this->dimension() == 1 )
399 ✗ FEL_ERROR("this structural model is not supposed to work on points");
400
401 204 m_curvFeVel = m_listCurvilinearFiniteElement[m_iVelocity];
402 204 m_curvFePre = m_listCurvilinearFiniteElement[m_iPressure];
403
404 204 m_displacementEF .initialize( DOF_FIELD, *m_curvFePre, 1 ); //current displacement
405 204 m_displacementOldEF .initialize( DOF_FIELD, *m_curvFePre, 1 ); //old displacement
406 204 m_gradUNEF .initialize( DOF_FIELD, *m_curvFeVel, this->dimension() ); //gradient of the normal velocity
407 204 m_TgradUNEF .initialize( DOF_FIELD, *m_curvFeVel, this->dimension() ); //gradient of the normal velocity projected onto the tangent plane
408 204 m_iopEF .initialize( DOF_FIELD, *m_curvFePre, 1 );
409 204 m_normalFieldEF .initialize( DOF_FIELD, *m_curvFeVel, this->dimension() );
410
411
2/2
✓ Branch 1 taken 52 times.
✓ Branch 2 taken 152 times.
 204 if ( this->dimension() == 3 ) {
412 52 m_meanCurvEF .initialize( DOF_FIELD, *m_curvFeVel, 3 ); //1: mean curv; 2: gauss curv; 3: 2*m_meanCurv^2-m_gaussCurv;
413 52 m_koiterCoefficientsEF .initialize( DOF_FIELD, *m_curvFeVel, 3 ); //c_1,c_2,c_3
414 52 m_firstLapCoeffEF .initialize( DOF_FIELD, *m_curvFeVel, 1 );
415 52 m_secondLapCoeffEF .initialize( DOF_FIELD, *m_curvFeVel, 1 );
416 52 m_thirdLapCoeffEF .initialize( DOF_FIELD, *m_curvFeVel, 1 );
417 52 m_minLapFibersCoeffEF .initialize( DOF_FIELD, *m_curvFeVel, 1 );
418 52 m_maxLapFibersCoeffEF .initialize( DOF_FIELD, *m_curvFeVel, 1 );
419 52 m_etagradunEF .initialize(DOF_FIELD, *m_curvFeVel, m_curvFeVel->numCoor() );
420 }
421 204 }
422
423 202344 void LinearProblemNSSimplifiedFSI::userElementComputeNaturalBoundaryCondition( const std::vector<Point*>& elemPoint, const std::vector<felInt>& /*elemIdPoint*/,felInt& iel, int label) {
424 BoundaryCondition* BC;
425
2/2
✓ Branch 1 taken 202344 times.
✓ Branch 2 taken 202344 times.
 404688 for (std::size_t iEmbedFSI=0; iEmbedFSI<m_boundaryConditionList.numEmbedFSIBoundaryCondition(); iEmbedFSI++ ) { //loop on all the surface regions that has the EmbedFSI BC
426 202344 BC = m_boundaryConditionList.EmbedFSI(iEmbedFSI);//get the the BC
427
3/4
✓ Branch 7 taken 202344 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 196477 times.
✓ Branch 11 taken 5867 times.
 202344 if ( std::find( BC->listLabel().begin(), BC->listLabel().end(),label) != BC->listLabel().end() ) {
428 double dispValue, dispValueOld, springCoeff;
429 196477 m_curvFeVel ->updateMeas(0,elemPoint);
430 196477 m_curvFePre ->updateMeas(0,elemPoint);
431
3/6
✓ Branch 3 taken 196477 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 196477 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 196477 times.
✗ Branch 10 not taken.
 196477 m_normalFieldEF .setValue( m_seqVecs.Get("normalField") , *m_curvFeVel , iel, m_iVelocity, m_ao, dof());
432
3/6
✓ Branch 3 taken 196477 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 196477 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 196477 times.
✗ Branch 10 not taken.
 196477 m_displacementEF .setValue( m_seqVecs.Get("displacement") , *m_curvFePre, iel, m_iPressure, m_ao, dof());
433
3/6
✓ Branch 3 taken 196477 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 196477 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 196477 times.
✗ Branch 10 not taken.
 196477 m_displacementOldEF.setValue( m_seqVecs.Get("displacementOld") , *m_curvFePre, iel, m_iPressure, m_ao, dof());
434
3/6
✓ Branch 3 taken 196477 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 196477 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 196477 times.
✗ Branch 10 not taken.
 196477 m_gradUNEF .setValue( m_seqVecs.Get("meanGradUN") , *m_curvFeVel , iel, m_iVelocity, m_ao, dof());
435 196477 computeTangentComponent(); //m_TgradUNEF
436
2/2
✓ Branch 0 taken 3040 times.
✓ Branch 1 taken 193437 times.
 196477 if ( m_iopCoupling )
437
3/6
✓ Branch 3 taken 3040 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 3040 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 3040 times.
✗ Branch 10 not taken.
 3040 m_iopEF.setValue( m_seqVecs.Get("currentIop"), *m_curvFePre, iel, m_iPressure, m_ao, dof());
438
439
2/3
✓ Branch 1 taken 3040 times.
✓ Branch 2 taken 193437 times.
✗ Branch 3 not taken.
 196477 switch ( this->dimension() ) {
440 3040 case 2:
441 { //======================================================================//
442 // in 2D we implemented a much simpler model,
443 // basically a first order transpiration + a linear koiter model with no fibers
444 //
445 //================================ Variables declaration================//
446 /*
447 betaS = E*h/(1 - nu^2);
448 alphaS = rho h / deltaT;
449 */
450
2/2
✓ Branch 1 taken 6080 times.
✓ Branch 2 taken 3040 times.
 9120 for ( int iLocDof(0); iLocDof < m_curvFeVel->numDof(); iLocDof++) {
451 6080 dispValue = m_displacementEF.val(0,iLocDof);
452 6080 dispValueOld = m_displacementOldEF.val(0,iLocDof);
453
454 //===================================================================================
455 // \mathcal{J}_1
456 6080 springCoeff = m_betaS / ( std::pow(m_radius,2) );
457
458 6080 m_elemFieldAlphaForEmbedFSI[iEmbedFSI].val(0,iLocDof) = springCoeff * m_deltaT + m_alphaS + m_viscS;
459
2/2
✓ Branch 1 taken 12160 times.
✓ Branch 2 taken 6080 times.
 18240 for ( int iCoor(0); iCoor < m_curvFeVel->numCoor(); iCoor++) {
460 12160 m_elemFieldEmbedFSI[iEmbedFSI].val(iCoor, iLocDof) =
461 12160 m_normalFieldEF.val(iCoor,iLocDof) * (
462 12160 - ( springCoeff - m_alphaS/m_fstransient->timeStep ) * dispValue // (-J_1 + rho h /delta T^2)\eta^k
463 12160 - m_alphaS/m_fstransient->timeStep * dispValueOld
464 12160 - m_referencePressure // -p^{ref}
465 12160 - m_iopEF.val(0,iLocDof)
466 ) //normal direction: displacement equation
467 12160 - m_elemFieldAlphaForEmbedFSI[iEmbedFSI].val(0,iLocDof) * dispValue * m_TgradUNEF.val(iCoor,iLocDof); //tangent direction: transpiration
468 }
469 }
470
3/6
✓ Branch 4 taken 3040 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3040 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3040 times.
✗ Branch 11 not taken.
 3040 m_elemFieldNormalForEmbedFSI[iEmbedFSI].setValue( m_seqVecs.Get("normalField"), *m_curvFeVel, iel, m_iVelocity, m_ao, dof());
471 }
472 3040 break;
473 193437 case 3:
474 {
475
1/2
✓ Branch 1 taken 193437 times.
✗ Branch 2 not taken.
 193437 ElementField currentNormalVelocity;
476
1/2
✓ Branch 1 taken 193437 times.
✗ Branch 2 not taken.
 193437 currentNormalVelocity.initialize( DOF_FIELD, *m_curvFePre, 1 );
477
3/6
✓ Branch 3 taken 193437 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 193437 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 193437 times.
✗ Branch 10 not taken.
 193437 currentNormalVelocity.setValue(m_seqVecs.Get("normalVelocity") , *m_curvFePre, iel, m_iPressure, m_ao, dof());
478
479
2/4
✓ Branch 1 taken 193437 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 193437 times.
 193437 if ( FelisceParam::instance().switchSimplFSIOff ) {
480 ✗ for ( int iLocDof(0); iLocDof < m_curvFeVel->numDof(); iLocDof++) {
481 ✗ m_elemFieldAlphaForEmbedFSI[iEmbedFSI].val(0,iLocDof) = .1;
482 ✗ for ( int iCoor(0); iCoor < m_curvFeVel->numCoor(); iCoor++) {
483 ✗ m_elemFieldEmbedFSI[iEmbedFSI].val(iCoor, iLocDof) = m_normalFieldEF.val(iCoor,iLocDof) * (- m_referencePressure
484 ✗ -m_iopEF.val(0,iLocDof)
485 ✗ +m_elemFieldAlphaForEmbedFSI[iEmbedFSI].val(0,iLocDof)*currentNormalVelocity.val(0,iLocDof)
486 );
487 }
488 }
489 ✗ m_elemFieldNormalForEmbedFSI[iEmbedFSI].setValue( m_seqVecs.Get("normalField"), *m_curvFeVel, iel, m_iVelocity, m_ao, dof());
490 } else {
491 //================================ Variables declaration================//
492 double normGradDisp, cNormGradDisp, cSquaredNormGradDisp, normGradDispMin, normGradDispMax, extraTension;
493
1/2
✓ Branch 3 taken 193437 times.
✗ Branch 4 not taken.
 193437 m_curv.update(m_normalFieldEF, m_curvFeVel->covMetric[0], m_curvFeVel->m_covBasis[0]);
494 /*
495 betaS = E*h/(1 - nu^2);
496 alphaS = rho h / deltaT;
497 fiberTens = k_0(S);
498 fiberYoung = k_1;
499 */
500 // -k_0(\bar S) + k_0(S)
501
3/4
✓ Branch 0 taken 181621 times.
✓ Branch 1 taken 11816 times.
✓ Branch 4 taken 181621 times.
✗ Branch 5 not taken.
 193437 extraTension = - ( m_autoregulated ? control(m_fstransient->time): 0. );
502
503
504
505
1/2
✓ Branch 1 taken 193437 times.
✗ Branch 2 not taken.
 193437 double fiberTensDensity = (m_fiberTensPassive + extraTension) * m_fiberDensity[ m_label2embedFSI[label] ];
506
1/2
✓ Branch 1 taken 193437 times.
✗ Branch 2 not taken.
 193437 double fiberYoungDensity = m_fiberYoung * m_fiberDensity[m_label2embedFSI[label]];
507
2/2
✓ Branch 1 taken 580311 times.
✓ Branch 2 taken 193437 times.
 773748 for ( int iLocDof(0); iLocDof < m_curvFeVel->numDof(); iLocDof++) {
508
2/2
✓ Branch 1 taken 1740933 times.
✓ Branch 2 taken 580311 times.
 2321244 for ( int ausComp(0); ausComp < m_curvFeVel->numCoor(); ausComp++) {
509
3/6
✓ Branch 1 taken 1740933 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1740933 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1740933 times.
✗ Branch 8 not taken.
 1740933 m_etagradunEF.val(ausComp,iLocDof) = m_displacementEF.val(0,iLocDof)*m_gradUNEF.val(ausComp,iLocDof);
510 }
511 }
512
513 // read the coefficients of the model
514
3/6
✓ Branch 3 taken 193437 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 193437 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 193437 times.
✗ Branch 10 not taken.
 193437 m_koiterCoefficientsEF.setValue( m_seqVecs.Get("koiterCoefficients"), *m_curvFeVel, iel, m_iVelocity, m_ao, dof());
515
3/6
✓ Branch 3 taken 193437 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 193437 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 193437 times.
✗ Branch 10 not taken.
 193437 m_meanCurvEF.setValue( m_seqVecs.Get("meanCurvature"), *m_curvFeVel, iel, m_iVelocity, m_ao, dof());
516
517 // compute the squared norm of the gradient
518
2/4
✓ Branch 1 taken 193437 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 193437 times.
✗ Branch 5 not taken.
 193437 normGradDisp = Curvatures::mNormOfGradient(m_displacementEF, m_curv.invMetricTensor() );
519 // compute the squared norm of the gradient weigthed by C
520
2/4
✓ Branch 1 taken 193437 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 193437 times.
✗ Branch 5 not taken.
 193437 cNormGradDisp = Curvatures::mNormOfGradient(m_displacementEF, m_curv.cMatrix() );
521 // compute the squared norm of the gradient weigthed by C^2
522
2/4
✓ Branch 1 taken 193437 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 193437 times.
✗ Branch 5 not taken.
 193437 cSquaredNormGradDisp = Curvatures::mNormOfGradient(m_displacementEF, m_curv.cSquaredMatrix() );
523 // compute the squared norm of the gradient in the min direction
524
2/4
✓ Branch 1 taken 193437 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 193437 times.
✗ Branch 5 not taken.
 193437 normGradDispMin = Curvatures::mNormOfGradient(m_displacementEF, m_curv.minProjTensor() );
525 // compute the squared norm of the gradient in the max direction
526
2/4
✓ Branch 1 taken 193437 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 193437 times.
✗ Branch 5 not taken.
 193437 normGradDispMax = Curvatures::mNormOfGradient(m_displacementEF, m_curv.maxProjTensor());
527
528
2/2
✓ Branch 1 taken 580311 times.
✓ Branch 2 taken 193437 times.
 773748 for ( int iLocDof(0); iLocDof < m_curvFeVel->numDof(); iLocDof++) {
529
1/2
✓ Branch 1 taken 580311 times.
✗ Branch 2 not taken.
 580311 dispValue = m_displacementEF.val(0,iLocDof);
530
1/2
✓ Branch 1 taken 580311 times.
✗ Branch 2 not taken.
 580311 dispValueOld = m_displacementOldEF.val(0,iLocDof);
531
532 //===================================================================================
533 // \mathcal{J}_1
534 // cf. Articles_Aletti/SimplifiedFSImodel and autoregulation/Contents/nonlinearkoiter.pdf p.5 (discretization in time)
535 1740933 springCoeff = 0.5 * m_betaS * (
536
1/2
✓ Branch 1 taken 580311 times.
✗ Branch 2 not taken.
 580311 2.0 * m_koiterCoefficientsEF.val(0,iLocDof) // 2c_1
537
1/2
✓ Branch 1 taken 580311 times.
✗ Branch 2 not taken.
 580311 - 3.0 * m_koiterCoefficientsEF.val(1,iLocDof) * dispValue * m_nonLinear // -3c_2 \eta
538
1/2
✓ Branch 1 taken 580311 times.
✗ Branch 2 not taken.
 580311 + 4.0 * m_koiterCoefficientsEF.val(2,iLocDof) * std::pow(dispValue,2) * m_nonLinear // 4c_3 \eta^2
539
1/2
✓ Branch 1 taken 580311 times.
✗ Branch 2 not taken.
 580311 + 2.0 * m_poissonS * m_meanCurvEF.val(2,iLocDof) * normGradDisp * m_nonLinear // 2\nu(2\rho_1^2-\rho_2)||\nabla\eta||^2
540 580311 + (1-m_poissonS) * cSquaredNormGradDisp * m_nonLinear // (1-\nu)||\nabla\eta||^2_{C^2}
541 )
542 580311 + (m_curv.coeffFibersLinearMin() + m_curv.coeffFibersLinearMax())*fiberTensDensity // This two lines should be more different if fibers were not aligned with curv dir
543 580311 + (m_curv.coeffFibersLinearMin() + m_curv.coeffFibersLinearMax())*fiberYoungDensity // This two lines should be more different if fibers were not aligned with curv dir
544 580311 + 0.5*(m_curv.coeffFibersLinearMin() * normGradDispMin + m_curv.coeffFibersLinearMax() * normGradDispMin)*fiberYoungDensity //ok
545 580311 - 3./2.*(m_curv.coeffFibersQuadraticMin() + m_curv.coeffFibersCubicMax())*fiberYoungDensity*dispValue // ok
546 580311 + 0.5*(m_curv.coeffFibersQuadraticMin() + m_curv.coeffFibersCubicMax())*fiberYoungDensity*std::pow(dispValue,2); // ok
547 //TODO we still do not compute the mean value in the nodes of the mesh for fibers coefficients
548
549 //===================================================================================
550 // all these quantities called LapCoeff refers to $\tens J(\eta,\nabla\eta)$
551 // for $\tens I$
552 1160622 m_firstLapCoeffEF.val(0,iLocDof) =
553
1/2
✓ Branch 1 taken 580311 times.
✗ Branch 2 not taken.
 580311 - 4.0 * m_poissonS * m_meanCurvEF.val(0,iLocDof) * dispValue
554
1/2
✓ Branch 1 taken 580311 times.
✗ Branch 2 not taken.
 580311 + 2.0 * m_poissonS * m_meanCurvEF.val(2,iLocDof) * std::pow(dispValue,2)
555
1/2
✓ Branch 1 taken 580311 times.
✗ Branch 2 not taken.
 580311 + normGradDisp;
556
557 // for $\tens C$
558
1/2
✓ Branch 1 taken 580311 times.
✗ Branch 2 not taken.
 580311 m_secondLapCoeffEF.val(0,iLocDof) = - 2.0 * (1-m_poissonS) * dispValue;
559
560 // for $\tens C^2
561
1/2
✓ Branch 2 taken 580311 times.
✗ Branch 3 not taken.
 580311 m_thirdLapCoeffEF.val(0,iLocDof) = (1-m_poissonS) * std::pow(dispValue,2);
562
563 // for $\varrho_{min} \tens P_{min}$
564 1160622 m_minLapFibersCoeffEF.val(0,iLocDof) =
565 fiberTensDensity
566 580311 - m_curv.minCurv() * dispValue * fiberYoungDensity
567 580311 + 0.5 * std::pow( m_curv.minCurv() * dispValue, 2) * fiberYoungDensity
568
1/2
✓ Branch 1 taken 580311 times.
✗ Branch 2 not taken.
 580311 + 0.5 * normGradDispMin * fiberYoungDensity;
569
570 // for $\varrho_{max} \tens P_{max}$
571 1160622 m_maxLapFibersCoeffEF.val(0,iLocDof) =
572 fiberTensDensity
573 580311 - m_curv.maxCurv() * dispValue * fiberYoungDensity
574 580311 + 0.5 * std::pow( m_curv.maxCurv() * dispValue, 2) * fiberYoungDensity
575
1/2
✓ Branch 1 taken 580311 times.
✗ Branch 2 not taken.
 580311 + 0.5 * normGradDispMax * fiberYoungDensity;
576
577
1/2
✓ Branch 2 taken 580311 times.
✗ Branch 3 not taken.
 580311 m_elemFieldAlphaForEmbedFSI[iEmbedFSI].val(0,iLocDof) = springCoeff * m_deltaT + m_alphaS + m_viscS; // ok
578 //===================================================================================
579
2/2
✓ Branch 1 taken 1740933 times.
✓ Branch 2 taken 580311 times.
 2321244 for ( int iCoor(0); iCoor < m_curvFeVel->numCoor(); iCoor++) {
580 //TODO the terms \Delta t \nabla_nn u^k are missing
581
1/2
✓ Branch 1 taken 1740933 times.
✗ Branch 2 not taken.
 1740933 m_elemFieldEmbedFSI[iEmbedFSI].val(iCoor, iLocDof) = //in the paper this is called \vec l^{k,k-1}
582
1/2
✓ Branch 1 taken 1740933 times.
✗ Branch 2 not taken.
 1740933 m_normalFieldEF.val(iCoor,iLocDof) * (
583 1740933 - ( springCoeff - m_alphaS/m_fstransient->timeStep ) * dispValue // (-J_1 + rho h /delta T^2)\eta^k
584 1740933 - m_alphaS/m_fstransient->timeStep * dispValueOld // ok
585 1740933 - m_referencePressure // -p^{ref}
586
1/2
✓ Branch 1 taken 1740933 times.
✗ Branch 2 not taken.
 1740933 - m_iopEF.val(0,iLocDof)
587 1740933 + 0.5 * m_betaS * ( 1.0 - m_poissonS ) * cNormGradDisp * m_nonLinear // -j0
588
1/2
✓ Branch 1 taken 1740933 times.
✗ Branch 2 not taken.
 1740933 + 0.5 * m_betaS * 2.0 * m_poissonS * m_meanCurvEF.val(0,iLocDof) * normGradDisp * m_nonLinear //-jo
589 1740933 + normGradDispMin * ( 0.5 * m_curv.coeffFiberProportionMin() * m_curv.minCurv() ) * fiberYoungDensity //-jo
590 1740933 + normGradDispMax * ( 0.5 * m_curv.coeffFiberProportionMax() * m_curv.maxCurv() ) * fiberYoungDensity //-jo
591 1740933 - extraTension * ( m_curv.coeffFiberProportionMin() * m_curv.minCurv() + m_curv.coeffFiberProportionMax() * m_curv.maxCurv()) // -jo
592 ) //normal direction: displacement equation
593
2/4
✓ Branch 2 taken 1740933 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1740933 times.
✗ Branch 6 not taken.
 1740933 - m_elemFieldAlphaForEmbedFSI[iEmbedFSI].val(0,iLocDof) * dispValue * m_TgradUNEF.val(iCoor,iLocDof); //tangent direction: transpiration
594 }
595
596
597 }
598
3/6
✓ Branch 4 taken 193437 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 193437 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 193437 times.
✗ Branch 11 not taken.
 193437 m_elemFieldNormalForEmbedFSI[iEmbedFSI].setValue( m_seqVecs.Get("normalField"), *m_curvFeVel, iel, m_iVelocity, m_ao, dof());
599
600
601
3/4
✓ Branch 1 taken 193437 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 181621 times.
✓ Branch 4 taken 11816 times.
 193437 if ( FelisceParam::instance().nonLinearKoiterModel ) {
602 // first laplacian contribution
603
3/6
✓ Branch 3 taken 181621 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 181621 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 181621 times.
✗ Branch 10 not taken.
 181621 m_elementVectorBD[0]->sGrad_phi_i_tensor_sGrad_f( -FelisceParam::instance().penValueForEmbedFSI * 0.5 * m_betaS, m_firstLapCoeffEF, m_displacementEF, m_normalFieldEF, *m_curvFeVel, m_curv.invMetricTensor() );
604
3/6
✓ Branch 3 taken 181621 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 181621 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 181621 times.
✗ Branch 10 not taken.
 181621 m_elementMatBD[0]->sGrad_psi_j_tensor_sGrad_phi_i( FelisceParam::instance().penValueForEmbedFSI * 0.5 * m_betaS * m_deltaT, m_firstLapCoeffEF, m_normalFieldEF, *m_curvFeVel, m_curv.invMetricTensor() );
605
606 // second laplacian contribution
607
3/6
✓ Branch 3 taken 181621 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 181621 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 181621 times.
✗ Branch 10 not taken.
 181621 m_elementVectorBD[0]->sGrad_phi_i_tensor_sGrad_f( -FelisceParam::instance().penValueForEmbedFSI * 0.5 * m_betaS, m_secondLapCoeffEF, m_displacementEF, m_normalFieldEF, *m_curvFeVel, m_curv.cMatrix() );
608
3/6
✓ Branch 3 taken 181621 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 181621 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 181621 times.
✗ Branch 10 not taken.
 181621 m_elementMatBD[0]->sGrad_psi_j_tensor_sGrad_phi_i( FelisceParam::instance().penValueForEmbedFSI * 0.5 * m_betaS * m_deltaT, m_secondLapCoeffEF, m_normalFieldEF, *m_curvFeVel, m_curv.cMatrix() );
609
610 // third laplacian contribution
611
3/6
✓ Branch 3 taken 181621 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 181621 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 181621 times.
✗ Branch 10 not taken.
 181621 m_elementVectorBD[0]->sGrad_phi_i_tensor_sGrad_f( -FelisceParam::instance().penValueForEmbedFSI * 0.5 * m_betaS, m_thirdLapCoeffEF, m_displacementEF, m_normalFieldEF, *m_curvFeVel, m_curv.cSquaredMatrix() );
612
3/6
✓ Branch 3 taken 181621 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 181621 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 181621 times.
✗ Branch 10 not taken.
 181621 m_elementMatBD[0]->sGrad_psi_j_tensor_sGrad_phi_i( FelisceParam::instance().penValueForEmbedFSI * 0.5 * m_betaS * m_deltaT, m_thirdLapCoeffEF, m_normalFieldEF, *m_curvFeVel, m_curv.cSquaredMatrix() );
613 }
614
615 // laplacian terms due to fibers
616 // min direction
617
3/6
✓ Branch 3 taken 193437 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 193437 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 193437 times.
✗ Branch 10 not taken.
 193437 m_elementVectorBD[0]->sGrad_phi_i_tensor_sGrad_f( -FelisceParam::instance().penValueForEmbedFSI, m_minLapFibersCoeffEF, m_displacementEF, m_normalFieldEF, *m_curvFeVel, m_curv.minFiberTensor() );
618
3/6
✓ Branch 3 taken 193437 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 193437 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 193437 times.
✗ Branch 10 not taken.
 193437 m_elementMatBD[0]->sGrad_psi_j_tensor_sGrad_phi_i( FelisceParam::instance().penValueForEmbedFSI * m_deltaT, m_minLapFibersCoeffEF, m_normalFieldEF, *m_curvFeVel, m_curv.minFiberTensor());
619
620 // max direction
621
3/6
✓ Branch 3 taken 193437 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 193437 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 193437 times.
✗ Branch 10 not taken.
 193437 m_elementVectorBD[0]->sGrad_phi_i_tensor_sGrad_f( -FelisceParam::instance().penValueForEmbedFSI, m_maxLapFibersCoeffEF, m_displacementEF, m_normalFieldEF, *m_curvFeVel, m_curv.maxFiberTensor() );
622
3/6
✓ Branch 3 taken 193437 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 193437 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 193437 times.
✗ Branch 10 not taken.
 193437 m_elementMatBD[0]->sGrad_psi_j_tensor_sGrad_phi_i( FelisceParam::instance().penValueForEmbedFSI * m_deltaT, m_maxLapFibersCoeffEF, m_normalFieldEF, *m_curvFeVel, m_curv.maxFiberTensor());
623
624
1/2
✓ Branch 1 taken 193437 times.
✗ Branch 2 not taken.
 193437 double stabParam = FelisceParam::instance().tssParam;
625
3/6
✓ Branch 3 taken 193437 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 193437 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 193437 times.
✗ Branch 11 not taken.
 193437 m_elementMatBD[0]->transpirationSurfStab( FelisceParam::instance().penValueForEmbedFSI * stabParam, m_elemFieldAlphaForEmbedFSI[iEmbedFSI], m_normalFieldEF, *m_curvFeVel, m_curv.invMetricTensor());
626
2/4
✓ Branch 1 taken 193437 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 193437 times.
✗ Branch 4 not taken.
 193437 if ( FelisceParam::instance().tssConsistent )
627
3/6
✓ Branch 3 taken 193437 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 193437 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 193437 times.
✗ Branch 11 not taken.
 193437 m_elementVectorBD[0]->transpirationSurfStabRHS( - FelisceParam::instance().penValueForEmbedFSI * stabParam, m_elemFieldAlphaForEmbedFSI[iEmbedFSI], m_normalFieldEF, m_etagradunEF, *m_curvFeVel, m_curv.invMetricTensor());
628 }
629 193437 }
630 193437 break;
631 ✗ default:
632 ✗ FEL_ERROR("strange dimension parameter");
633 }
634 }
635 }
636 202344 }
637
638 196477 void LinearProblemNSSimplifiedFSI::computeTangentComponent() {
639
2/2
✓ Branch 1 taken 586391 times.
✓ Branch 2 taken 196477 times.
 782868 for ( int iLocDof(0); iLocDof < m_curvFeVel->numDof(); iLocDof++) {
640
2/2
✓ Branch 1 taken 1753093 times.
✓ Branch 2 taken 586391 times.
 2339484 for ( int iCoor(0); iCoor < m_curvFeVel->numCoor(); iCoor++) {
641 1753093 m_TgradUNEF.val(iCoor,iLocDof)=0.0;
642
2/2
✓ Branch 1 taken 5247119 times.
✓ Branch 2 taken 1753093 times.
 7000212 for ( int iCoorAus(0); iCoorAus < m_curvFeVel->numCoor(); iCoorAus++) {
643 // kronecker delta
644 5247119 double delta=0.0;
645
2/2
✓ Branch 0 taken 1753093 times.
✓ Branch 1 taken 3494026 times.
 5247119 if( iCoor==iCoorAus )
646 1753093 delta=1.0;
647 5247119 m_TgradUNEF.val(iCoor,iLocDof) += ( delta - m_normalFieldEF.val(iCoor,iLocDof) * m_normalFieldEF.val(iCoorAus,iLocDof) ) * m_gradUNEF.val( iCoorAus, iLocDof );
648 }
649 }
650 }
651 196477 }
652
653 64 void LinearProblemNSSimplifiedFSI::initPerETWBD() {
654 64 m_curvFeVel = m_listCurvilinearFiniteElement[m_iVelocity];
655 64 m_curWBDFeVel = m_listCurrentFiniteElementWithBd[m_iVelocity];
656 64 }
657
658 193677 void LinearProblemNSSimplifiedFSI::updateFeWBD(const std::vector<Point*>& elemPoint,const std::vector<felInt>& elemIdPoint) {
659 193677 std::vector<Point*> elemPointVol;
660 193677 std::vector<felInt> elemIdPointVol;
661
1/2
✓ Branch 1 taken 193677 times.
✗ Branch 2 not taken.
 193677 m_curvFeVel->updateMeasNormal(0, elemPoint);
662
1/2
✓ Branch 5 taken 193677 times.
✗ Branch 6 not taken.
 193677 m_mesh[m_currentMesh]->getElementFromBdElem(elemIdPoint, elemIdPointVol, /* not used */elemPointVol, /*idLocFace*/m_auxiliaryInts[0], /*idElemVol*/m_auxiliaryInts[1]);
663
1/2
✓ Branch 1 taken 193677 times.
✗ Branch 2 not taken.
 193677 m_curWBDFeVel->updateFirstDeriv(0, elemPointVol);
664
1/2
✓ Branch 1 taken 193677 times.
✗ Branch 2 not taken.
 193677 m_curWBDFeVel->updateBdMeasNormal();
665 193677 }
666 193677 void LinearProblemNSSimplifiedFSI::meanGradUNComputer(felInt ielSupportDof){
667 double value;
668
2/2
✓ Branch 1 taken 580791 times.
✓ Branch 2 taken 193677 times.
 774468 for ( std::size_t iCoor(0); iCoor < ( std::size_t ) m_curvFeVel->numCoor(); iCoor++) {
669 //it is a std::vector with dimension = num of quad points on the bd elem
670 580791 std::vector<double> aus;
671
1/2
✓ Branch 3 taken 580791 times.
✗ Branch 4 not taken.
 580791 this->computeGradUN( m_curWBDFeVel, m_auxiliaryInts[1], m_iVelocity, aus, m_auxiliaryInts[0], iCoor);
672
673
1/2
✓ Branch 1 taken 580791 times.
✗ Branch 2 not taken.
 580791 double feMeas( m_curvFeVel->measure() );
674 580791 value = aus[0] * feMeas; //TODO here you should integrate to be more general (for P1 it is ok..), in the case of P1 the gradient is constant over the element
675
2/2
✓ Branch 1 taken 1741893 times.
✓ Branch 2 taken 580791 times.
 2322684 for ( std::size_t iSurfDof(0); iSurfDof < ( std::size_t ) m_curvFeVel->numDof(); iSurfDof++) {
676 felInt iDofGlob;
677
1/2
✓ Branch 2 taken 1741893 times.
✗ Branch 3 not taken.
 1741893 dof().loc2glob( ielSupportDof, iSurfDof, m_iVelocity, iCoor, iDofGlob);
678
1/2
✓ Branch 1 taken 1741893 times.
✗ Branch 2 not taken.
 1741893 AOApplicationToPetsc(m_ao,1,&iDofGlob);
679
3/6
✓ Branch 2 taken 1741893 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1741893 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 1741893 times.
✗ Branch 9 not taken.
 1741893 m_vecs.Get("meanGradUN").setValue(iDofGlob, value, ADD_VALUES);
680 }
681 580791 }
682 193677 }
683 64 void LinearProblemNSSimplifiedFSI::computeMeanGradUN() {
684
3/6
✓ Branch 2 taken 64 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 64 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 64 times.
✗ Branch 9 not taken.
 64 m_vecs.Get("meanGradUN").set(0.0);
685
3/6
✓ Branch 2 taken 64 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 64 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 64 times.
✗ Branch 9 not taken.
 64 m_seqVecs.Get("meanGradUN").set(0.0);
686
1/2
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
 64 m_auxiliaryInts.resize(2);
687
2/2
✓ Branch 1 taken 52 times.
✓ Branch 2 taken 12 times.
 64 if ( this->dimension() == 3 )
688
1/2
✓ Branch 3 taken 52 times.
✗ Branch 4 not taken.
 52 m_mesh[m_currentMesh]->buildFaces();
689 else
690
1/2
✓ Branch 3 taken 12 times.
✗ Branch 4 not taken.
 12 m_mesh[m_currentMesh]->buildEdges();
691
692 //TODO save this list of labels!
693 64 std::vector<felInt> labels;
694
2/2
✓ Branch 1 taken 64 times.
✓ Branch 2 taken 64 times.
 128 for (std::size_t iEmbedFSI=0; iEmbedFSI < m_boundaryConditionList.numEmbedFSIBoundaryCondition(); iEmbedFSI++ ) {
695 64 BoundaryCondition* BC = m_boundaryConditionList.EmbedFSI(iEmbedFSI);
696
2/2
✓ Branch 6 taken 636 times.
✓ Branch 7 taken 64 times.
 700 for(auto it_labelNumber = BC->listLabel().begin(); it_labelNumber != BC->listLabel().end(); it_labelNumber++) {
697
1/2
✓ Branch 2 taken 636 times.
✗ Branch 3 not taken.
 636 labels.push_back(*it_labelNumber);
698 }
699 }
700
701
1/2
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
 64 this->assemblyLoopBoundaryGeneral(&LinearProblemNSSimplifiedFSI::meanGradUNComputer,labels,&LinearProblemNSSimplifiedFSI::initPerETWBD,&LinearProblemNSSimplifiedFSI::updateFeWBD);
702
3/6
✓ Branch 2 taken 64 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 64 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 64 times.
✗ Branch 9 not taken.
 64 m_vecs.Get("meanGradUN").assembly();
703
7/14
✓ Branch 2 taken 64 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 64 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 64 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 64 times.
✗ Branch 13 not taken.
✓ Branch 16 taken 64 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 64 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 64 times.
✗ Branch 23 not taken.
 64 pointwiseDivide(m_vecs.Get("meanGradUN"), m_vecs.Get("meanGradUN"), m_vecs.Get("measStar"));
704
5/10
✓ Branch 2 taken 64 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 64 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 64 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 64 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 64 times.
✗ Branch 16 not taken.
 64 gatherVector( m_vecs.Get("meanGradUN"), m_seqVecs.Get("meanGradUN") );
705 64 }
706
707 580791 void LinearProblemNSSimplifiedFSI::computeGradUN( CurrentFiniteElementWithBd* feVol, felInt iel, int idVariable, std::vector<double> &grad, int ibdFace, int icomp) {
708 580791 grad.clear();
709
1/2
✓ Branch 2 taken 580791 times.
✗ Branch 3 not taken.
 580791 grad.resize( feVol->numQuadPointBdEle(ibdFace), 0. );
710
711 580791 felInt numCoor = feVol->numCoor(); // dimension of the space, (3).
712 580791 felInt numDof = feVol->numDof(); // number of dof on the volume
713 felInt idDof; // the index of the Dof in the global numbering
714 std::size_t iQuadVol;
715
716
2/2
✓ Branch 1 taken 3483306 times.
✓ Branch 2 taken 580791 times.
 4064097 for( std::size_t ilocg = 0; ilocg < grad.size(); ilocg++) {
717 3483306 iQuadVol = ilocg + feVol->indexQuadPoint( ibdFace + 1 ); // face ibd starts at ibd+1 (0 for internal quad points)
718
719
2/2
✓ Branch 0 taken 13931784 times.
✓ Branch 1 taken 3483306 times.
 17415090 for ( int idof = 0; idof < numDof; idof++) {
720 double solValue;//the value
721
1/2
✓ Branch 2 taken 13931784 times.
✗ Branch 3 not taken.
 13931784 dof().loc2glob( iel, idof, idVariable, icomp, idDof);
722
1/2
✓ Branch 1 taken 13931784 times.
✗ Branch 2 not taken.
 13931784 AOApplicationToPetsc( m_ao, 1, &idDof);
723
1/2
✓ Branch 2 taken 13931784 times.
✗ Branch 3 not taken.
 13931784 this->sequentialSolution().getValues(1, &idDof, &solValue);
724
2/2
✓ Branch 0 taken 41791032 times.
✓ Branch 1 taken 13931784 times.
 55722816 for ( int iauscomp = 0; iauscomp < numCoor; iauscomp++) {
725
3/6
✓ Branch 2 taken 41791032 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 41791032 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 41791032 times.
✗ Branch 10 not taken.
 41791032 grad[ilocg] += solValue * feVol->dPhi[ iQuadVol ]( iauscomp, idof ) * feVol->bdEle( ibdFace ).normal[ ilocg ]( iauscomp );
726 }
727 }
728 }
729 580791 }
730
731 152 void LinearProblemNSSimplifiedFSI::applyEssentialBoundaryConditionDerivedProblem(int rank, FlagMatrixRHS flagMatrixRHS) {
732
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 152 times.
 152 if ( FelisceParam::instance().zeroDisplacementAtBoarders ) {
733 ✗ if (FelisceParam::verbose() > 2)
734 ✗ std::cout<<"["<<rank<<"]--Applying zero-Dirichlet BCs on the velocity at the interface between inlet outlet and surface."<<std::endl;
735 ✗ getRings();
736 ✗ double bigValue=FelisceParam::instance().penValueForEmbedFSI*1e10;
737
738 ✗ for(auto itDofBC = m_idDofRings.begin(); itDofBC != m_idDofRings.end(); ++itDofBC) {
739 ✗ if( (flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs) || (flagMatrixRHS == FlagMatrixRHS::only_matrix))
740 ✗ matrix(0).setValue(*itDofBC, *itDofBC, bigValue, INSERT_VALUES);
741 ✗ if( (flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs) || (flagMatrixRHS == FlagMatrixRHS::only_rhs))
742 ✗ vector().setValue(*itDofBC,0.0, INSERT_VALUES);
743 }
744 }
745 152 }
746
747 22419 void LinearProblemNSSimplifiedFSI::curvaturesComputer(felInt ielSupportDof) {
748
1/2
✓ Branch 1 taken 22419 times.
✗ Branch 2 not taken.
 22419 ElementField normal;
749
750
1/2
✓ Branch 2 taken 22419 times.
✗ Branch 3 not taken.
 22419 normal.initialize( DOF_FIELD, *m_curvFeVel, m_curvFeVel->numCoor() );
751
3/6
✓ Branch 3 taken 22419 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 22419 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 22419 times.
✗ Branch 10 not taken.
 22419 normal.setValue( m_seqVecs.Get("normalField"), *m_curvFeVel, ielSupportDof, m_iVelocity, m_ao, dof());
752
1/2
✓ Branch 3 taken 22419 times.
✗ Branch 4 not taken.
 22419 m_curv.update(normal, m_curvFeVel->covMetric[0], m_curvFeVel->m_covBasis[0]);
753
754 // The curvatures have been computed in the element, but we want to define them on the nodes, computing a mean over the elements sharing that node
755 // We use a std::vector field, where each component stores a different thing
756 double ausMeanCurvature,ausKoiterCoeff,currCompMinEigV,currCompMaxEigV;
757
2/2
✓ Branch 1 taken 67257 times.
✓ Branch 2 taken 22419 times.
 89676 for ( std::size_t iSurfDof(0); iSurfDof < ( std::size_t ) m_curvFeVel->numDof(); iSurfDof++) {
758
2/2
✓ Branch 1 taken 201771 times.
✓ Branch 2 taken 67257 times.
 269028 for ( std::size_t iCoor(0); iCoor < ( std::size_t ) m_curvFeVel->numCoor(); iCoor++) {
759 felInt iDofGlob;
760
1/2
✓ Branch 2 taken 201771 times.
✗ Branch 3 not taken.
 201771 dof().loc2glob( ielSupportDof, iSurfDof, m_iVelocity, iCoor, iDofGlob);
761
1/2
✓ Branch 1 taken 201771 times.
✗ Branch 2 not taken.
 201771 AOApplicationToPetsc(m_ao,1,&iDofGlob);
762 201771 ausMeanCurvature=0.0;
763 201771 ausKoiterCoeff=0.0;
764 201771 currCompMaxEigV=0;
765 201771 currCompMinEigV=0;
766
2/2
✓ Branch 1 taken 1210626 times.
✓ Branch 2 taken 201771 times.
 1412397 for ( std::size_t iQuadBd(0.0) ; iQuadBd < ( std::size_t )m_curvFeVel->numQuadraturePoint(); iQuadBd++) {
767
2/4
✓ Branch 1 taken 1210626 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1210626 times.
✗ Branch 4 not taken.
 1210626 if ( FelisceParam::instance().exportPrincipalDirectionsAndCurvature ) {
768
2/4
✓ Branch 1 taken 1210626 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1210626 times.
✗ Branch 5 not taken.
 1210626 currCompMaxEigV += m_curvFeVel->weightMeas( iQuadBd ) * m_curv.maxEigVCartComp(iCoor);
769
2/4
✓ Branch 1 taken 1210626 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1210626 times.
✗ Branch 5 not taken.
 1210626 currCompMinEigV += m_curvFeVel->weightMeas( iQuadBd ) * m_curv.minEigVCartComp(iCoor);
770 }
771
2/2
✓ Branch 0 taken 403542 times.
✓ Branch 1 taken 807084 times.
 1210626 if ( iCoor == 0 ) {
772
1/2
✓ Branch 1 taken 403542 times.
✗ Branch 2 not taken.
 403542 ausMeanCurvature += m_curvFeVel->weightMeas( iQuadBd ) * m_curv.meanCurv();// \rho_1
773
1/2
✓ Branch 1 taken 403542 times.
✗ Branch 2 not taken.
 403542 ausKoiterCoeff += m_curvFeVel->weightMeas( iQuadBd ) * m_curv.coeffKoiterLinear(m_poissonS);// C_1
774
2/2
✓ Branch 0 taken 403542 times.
✓ Branch 1 taken 403542 times.
 807084 } else if ( iCoor == 1 ) {
775
1/2
✓ Branch 1 taken 403542 times.
✗ Branch 2 not taken.
 403542 ausMeanCurvature += m_curvFeVel->weightMeas( iQuadBd ) * m_curv.gaussCurv();// \rho_2
776
1/2
✓ Branch 1 taken 403542 times.
✗ Branch 2 not taken.
 403542 ausKoiterCoeff += m_curvFeVel->weightMeas( iQuadBd ) * m_curv.coeffKoiterQuadratic(m_poissonS);// C_2
777
1/2
✓ Branch 0 taken 403542 times.
✗ Branch 1 not taken.
 403542 } else if ( iCoor == 2 ) {
778
1/2
✓ Branch 1 taken 403542 times.
✗ Branch 2 not taken.
 403542 ausMeanCurvature += m_curvFeVel->weightMeas( iQuadBd ) * m_curv.coeffKoiterMeanOfSquaredCurvatures();// 2*std::pow(m_meanCurv,2)-m_gaussCurv;
779
1/2
✓ Branch 1 taken 403542 times.
✗ Branch 2 not taken.
 403542 ausKoiterCoeff += m_curvFeVel->weightMeas( iQuadBd ) * m_curv.coeffKoiterCubic(m_poissonS);// C_3
780 }
781 }
782
3/6
✓ Branch 2 taken 201771 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 201771 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 201771 times.
✗ Branch 9 not taken.
 201771 m_vecs.Get("meanCurvature").setValue(iDofGlob, ausMeanCurvature, ADD_VALUES);
783
3/6
✓ Branch 2 taken 201771 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 201771 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 201771 times.
✗ Branch 9 not taken.
 201771 m_vecs.Get("koiterCoefficients").setValue(iDofGlob, ausKoiterCoeff, ADD_VALUES);
784
2/4
✓ Branch 1 taken 201771 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 201771 times.
✗ Branch 4 not taken.
 201771 if ( FelisceParam::instance().exportPrincipalDirectionsAndCurvature ) {
785
3/6
✓ Branch 2 taken 201771 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 201771 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 201771 times.
✗ Branch 9 not taken.
 201771 m_vecs.Get("maxEigVec").setValue(iDofGlob, currCompMaxEigV, ADD_VALUES);
786
3/6
✓ Branch 2 taken 201771 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 201771 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 201771 times.
✗ Branch 9 not taken.
 201771 m_vecs.Get("minEigVec").setValue(iDofGlob, currCompMinEigV, ADD_VALUES);
787 }
788 }
789 }
790 22419 }
791
792 /*!
793 \brief It compute the displacement as a scalar.
794
795 \param[out] max, the maximum of the displacement
796 \param[in] min, the minimum of the displacement
797
798 The computation is done via assemblyLoopBoundary() and scalarComputer().
799
800 The displacement is the integral over time of the scalar product of the velocity times the normal.\n
801 Therefore it is first stored as std::vector. The i-th component of such std::vector contain \f$ v_i n_i\f$.
802 This function sum the components and store them into the pressure block of displacement (both para and seq).\n
803 At the same tame also the max and min of the displacement are computed.
804 */
805 64 void LinearProblemNSSimplifiedFSI::computeScalarDisplacement(double& max, double& min) {
806 64 m_auxiliaryDoubles.clear();
807
1/2
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
 64 m_auxiliaryDoubles.push_back(-1.e30);
808
1/2
✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
 64 m_auxiliaryDoubles.push_back(1.e30);
809
810 64 assemblyLoopBoundary(&LinearProblemNSSimplifiedFSI::scalarComputer);
811
3/6
✓ Branch 2 taken 64 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 64 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 64 times.
✗ Branch 9 not taken.
 64 m_vecs.Get("displacement").assembly();
812
813
5/10
✓ Branch 2 taken 64 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 64 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 64 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 64 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 64 times.
✗ Branch 16 not taken.
 64 gatherVector(m_vecs.Get("displacement"),m_seqVecs.Get("displacement"));
814
815 64 max=m_auxiliaryDoubles[0];
816 64 min=m_auxiliaryDoubles[1];
817 64 }
818
819 /*!
820 \brief It compute the normalVelocity as a scalar
821
822 As the displacement also the normalVelocity is a scalar product of the velocity times the normal.\n
823 Therefore it is first stored as std::vector. The i-th component of such std::vector contain \f$ v_i n_i\f$.
824 This function sum the components and store them into the pressure block of normalVelocity (both para and seq).
825 */
826 216 void LinearProblemNSSimplifiedFSI::computeScalarNormalVelocity() {
827 216 assemblyLoopBoundary(&LinearProblemNSSimplifiedFSI::normalVelocityScalarComputer);
828
3/6
✓ Branch 2 taken 216 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 216 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 216 times.
✗ Branch 9 not taken.
 216 m_vecs.Get("normalVelocity").assembly();
829
5/10
✓ Branch 2 taken 216 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 216 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 216 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 216 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 216 times.
✗ Branch 16 not taken.
 216 gatherVector(m_vecs.Get("normalVelocity"),m_seqVecs.Get("normalVelocity"));
830 216 }
831
832 ✗ double LinearProblemNSSimplifiedFSI::computeL2Difference(std::string current, std::string old) {
833 ✗ if ( current == "zero" )
834 ✗ return std::sqrt( this->computeL2ScalarProduct(old,old) );
835 ✗ if ( old == "zero" )
836 ✗ return std::sqrt( this->computeL2ScalarProduct(current,current) );
837
838 ✗ m_auxiliaryStrings.clear();
839 ✗ m_auxiliaryStrings.push_back(current);
840 ✗ m_auxiliaryStrings.push_back(old);
841
842 ✗ m_auxiliaryDoubles.clear();
843 ✗ m_auxiliaryDoubles.push_back(0.0);
844 ✗ assemblyLoopBoundary(&LinearProblemNSSimplifiedFSI::normL2BoundaryDifferenceComputer);
845 ✗ double result(0);
846 ✗ MPI_Allreduce(&m_auxiliaryDoubles[0],&result,1,MPI_DOUBLE,MPI_SUM,MpiInfo::petscComm());
847 ✗ return std::sqrt(result);
848 }
849
850 152 std::pair<double,double> LinearProblemNSSimplifiedFSI::computeTestQuantities( bool lumped ) {
851 152 std::pair<double,double> resultP;
852 152 m_auxiliaryDoubles.clear();
853
1/2
✓ Branch 1 taken 152 times.
✗ Branch 2 not taken.
 152 m_auxiliaryDoubles.resize(4,0.0);
854
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 152 times.
 152 if ( lumped )
855 ✗ assemblyLoopBoundary(&LinearProblemNSSimplifiedFSI::testQuantitiesLumpedComputer);
856 else
857
1/2
✓ Branch 1 taken 152 times.
✗ Branch 2 not taken.
 152 assemblyLoopBoundary(&LinearProblemNSSimplifiedFSI::testQuantitiesComputer);
858
1/2
✓ Branch 2 taken 152 times.
✗ Branch 3 not taken.
 152 std::vector<double> result(4,0);
859
2/4
✓ Branch 1 taken 152 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 152 times.
✗ Branch 7 not taken.
 152 MPI_Allreduce(m_auxiliaryDoubles.data(),result.data(),4,MPI_DOUBLE,MPI_SUM,MpiInfo::petscComm());
860
861 using std::sqrt;
862 152 resultP.first = std::sqrt(result[0])/( std::sqrt(result[1]) + 1e-12);
863 152 resultP.second = std::sqrt(result[2])/( std::sqrt(result[3]) + 1e-12);
864 152 return resultP;
865 152 }
866
867 ✗ double LinearProblemNSSimplifiedFSI::computeL2ScalarProduct(std::string l, std::string r) {
868 ✗ m_auxiliaryStrings.clear();
869 ✗ m_auxiliaryStrings.push_back(l);
870 ✗ m_auxiliaryStrings.push_back(r);
871
872 ✗ m_auxiliaryDoubles.clear();
873 ✗ m_auxiliaryDoubles.push_back(0.0);
874 ✗ assemblyLoopBoundary(&LinearProblemNSSimplifiedFSI::scalarProductComputer);
875 ✗ double result(0);
876 ✗ MPI_Allreduce(&m_auxiliaryDoubles[0],&result,1,MPI_DOUBLE,MPI_SUM,MpiInfo::petscComm());
877 ✗ return result;
878 }
879
880 ✗ void LinearProblemNSSimplifiedFSI::normL2BoundaryDifferenceComputer(felInt ielSupportDof) {
881 double current,old,valuei,valuej;
882 ✗ for ( std::size_t iSurfDof(0); iSurfDof < ( std::size_t ) m_curvFePre->numDof(); iSurfDof++) {
883 felInt iDofGlobScalar;
884 ✗ dof().loc2glob( ielSupportDof, iSurfDof, m_iPressure, 0, iDofGlobScalar);
885 ✗ AOApplicationToPetsc(m_ao,1,&iDofGlobScalar);
886
887 ✗ m_seqVecs.Get( m_auxiliaryStrings[0] ).getValues(1,&iDofGlobScalar, &current);
888 ✗ m_seqVecs.Get( m_auxiliaryStrings[1] ).getValues(1,&iDofGlobScalar, &old);
889
890 ✗ valuei=current-old;
891 ✗ for ( std::size_t jSurfDof(0); jSurfDof < ( std::size_t ) m_curvFePre->numDof(); jSurfDof++) {
892 felInt jDofGlobScalar;
893 ✗ dof().loc2glob( ielSupportDof, jSurfDof, m_iPressure, 0, jDofGlobScalar);
894 ✗ AOApplicationToPetsc(m_ao,1,&jDofGlobScalar);
895 ✗ m_seqVecs.Get( m_auxiliaryStrings[0] ).getValues(1,&jDofGlobScalar, &current);
896 ✗ m_seqVecs.Get( m_auxiliaryStrings[1] ).getValues(1,&jDofGlobScalar, &old);
897
898 ✗ valuej=current-old;
899 ✗ for(int ig=0; ig < m_curvFePre->numQuadraturePoint(); ig++) {
900 ✗ m_auxiliaryDoubles[0] += valuei*valuej * m_curvFePre->phi[ig](iSurfDof) * m_curvFePre->phi[ig](jSurfDof) * m_curvFePre->weightMeas(ig);
901 }
902 }
903 }
904 }
905
906 void
907 3040 LinearProblemNSSimplifiedFSI::testQuantitiesComputer(felInt ielSupportDof) {
908
909 double uni,unoldi,ciopi,x1i;
910 double value0i,value1i,value2i,value3i;
911
912 double unj,unoldj,ciopj,x1j;
913 double value0j,value1j,value2j,value3j;
914
915
2/2
✓ Branch 1 taken 6080 times.
✓ Branch 2 taken 3040 times.
 9120 for ( std::size_t iSurfDof(0); iSurfDof < ( std::size_t ) m_curvFePre->numDof(); iSurfDof++) {
916 felInt iDofGlobScalar;
917
1/2
✓ Branch 2 taken 6080 times.
✗ Branch 3 not taken.
 6080 dof().loc2glob( ielSupportDof, iSurfDof, m_iPressure, 0, iDofGlobScalar);
918
1/2
✓ Branch 1 taken 6080 times.
✗ Branch 2 not taken.
 6080 AOApplicationToPetsc(m_ao,1,&iDofGlobScalar);
919
920
3/6
✓ Branch 2 taken 6080 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 6080 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 6080 times.
✗ Branch 9 not taken.
 6080 m_seqVecs.Get( "normalVelocity" ).getValues(1,&iDofGlobScalar, &uni);
921
3/6
✓ Branch 2 taken 6080 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 6080 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 6080 times.
✗ Branch 9 not taken.
 6080 m_seqVecs.Get( "normalVelocityOld" ).getValues(1,&iDofGlobScalar, &unoldi);
922
3/6
✓ Branch 2 taken 6080 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 6080 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 6080 times.
✗ Branch 9 not taken.
 6080 m_seqVecs.Get( "currentIop" ).getValues(1,&iDofGlobScalar, &ciopi);
923
3/6
✓ Branch 2 taken 6080 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 6080 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 6080 times.
✗ Branch 9 not taken.
 6080 m_seqVecs.Get( "x1" ).getValues(1,&iDofGlobScalar, &x1i);
924
925 6080 value0i = uni-unoldi;
926 6080 value1i = uni;
927 6080 value2i = ciopi-x1i;
928 6080 value3i = ciopi;
929
930
931
2/2
✓ Branch 1 taken 12160 times.
✓ Branch 2 taken 6080 times.
 18240 for ( std::size_t jSurfDof(0); jSurfDof < ( std::size_t ) m_curvFePre->numDof(); jSurfDof++) {
932 felInt jDofGlobScalar;
933
1/2
✓ Branch 2 taken 12160 times.
✗ Branch 3 not taken.
 12160 dof().loc2glob( ielSupportDof, jSurfDof, m_iPressure, 0, jDofGlobScalar);
934
1/2
✓ Branch 1 taken 12160 times.
✗ Branch 2 not taken.
 12160 AOApplicationToPetsc(m_ao,1,&jDofGlobScalar);
935
936
3/6
✓ Branch 2 taken 12160 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12160 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 12160 times.
✗ Branch 9 not taken.
 12160 m_seqVecs.Get( "normalVelocity" ).getValues(1,&jDofGlobScalar, &unj);
937
3/6
✓ Branch 2 taken 12160 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12160 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 12160 times.
✗ Branch 9 not taken.
 12160 m_seqVecs.Get( "normalVelocityOld" ).getValues(1,&jDofGlobScalar, &unoldj);
938
3/6
✓ Branch 2 taken 12160 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12160 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 12160 times.
✗ Branch 9 not taken.
 12160 m_seqVecs.Get( "currentIop" ).getValues(1,&jDofGlobScalar, &ciopj);
939
3/6
✓ Branch 2 taken 12160 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12160 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 12160 times.
✗ Branch 9 not taken.
 12160 m_seqVecs.Get( "x1" ).getValues(1,&jDofGlobScalar, &x1j);
940
941 12160 value0j = unj-unoldj;
942 12160 value1j = unj;
943 12160 value2j = ciopj-x1j;
944 12160 value3j = ciopj;
945
946
2/2
✓ Branch 1 taken 36480 times.
✓ Branch 2 taken 12160 times.
 48640 for(int ig=0; ig < m_curvFePre->numQuadraturePoint(); ig++) {
947
3/6
✓ Branch 2 taken 36480 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 36480 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 36480 times.
✗ Branch 10 not taken.
 36480 double w = m_curvFePre->phi[ig](iSurfDof) * m_curvFePre->phi[ig](jSurfDof) * m_curvFePre->weightMeas(ig);
948 36480 m_auxiliaryDoubles[0] += value0i*value0j * w;
949 36480 m_auxiliaryDoubles[1] += value1i*value1j * w;
950 36480 m_auxiliaryDoubles[2] += value2i*value2j * w;
951 36480 m_auxiliaryDoubles[3] += value3i*value3j * w;
952 }
953 }
954 }
955 3040 }
956
957 void
958 ✗ LinearProblemNSSimplifiedFSI::testQuantitiesLumpedComputer(felInt ielSupportDof) {
959
960 double uni,unoldi,ciopi,x1i;
961 double value0i,value1i,value2i,value3i;
962
963 ✗ for ( std::size_t iSurfDof(0); iSurfDof < ( std::size_t ) m_curvFePre->numDof(); iSurfDof++) {
964 felInt iDofGlobScalar;
965 ✗ dof().loc2glob( ielSupportDof, iSurfDof, m_iPressure, 0, iDofGlobScalar);
966 ✗ AOApplicationToPetsc(m_ao,1,&iDofGlobScalar);
967
968 ✗ m_seqVecs.Get( "normalVelocity" ).getValues(1,&iDofGlobScalar, &uni);
969 ✗ m_seqVecs.Get( "normalVelocityOld" ).getValues(1,&iDofGlobScalar, &unoldi);
970 ✗ m_seqVecs.Get( "currentIop" ).getValues(1,&iDofGlobScalar, &ciopi);
971 ✗ m_seqVecs.Get( "x1" ).getValues(1,&iDofGlobScalar, &x1i);
972
973 using std::pow;
974 ✗ value0i = std::pow( uni-unoldi ,2);
975 ✗ value1i = std::pow( uni ,2);
976 ✗ value2i = std::pow( ciopi-x1i ,2);
977 ✗ value3i = std::pow( ciopi ,2);
978
979 ✗ for(int ig=0; ig < m_curvFePre->numQuadraturePoint(); ig++) {
980 ✗ double w = std::pow( m_curvFePre->phi[ig](iSurfDof) ,2) * m_curvFePre->weightMeas(ig);
981 ✗ m_auxiliaryDoubles[0] += value0i * w;
982 ✗ m_auxiliaryDoubles[1] += value1i * w;
983 ✗ m_auxiliaryDoubles[2] += value2i * w;
984 ✗ m_auxiliaryDoubles[3] += value3i * w;
985 }
986 }
987 }
988
989 ✗ void LinearProblemNSSimplifiedFSI::scalarProductComputer(felInt ielSupportDof) {
990 double l,r;
991 ✗ for ( std::size_t iSurfDof(0); iSurfDof < ( std::size_t ) m_curvFePre->numDof(); iSurfDof++) {
992 felInt iDofGlobScalar;
993 ✗ dof().loc2glob( ielSupportDof, iSurfDof, m_iPressure, 0, iDofGlobScalar);
994 ✗ AOApplicationToPetsc(m_ao,1,&iDofGlobScalar);
995 ✗ m_seqVecs.Get( m_auxiliaryStrings[0] ).getValues(1,&iDofGlobScalar, &r);
996 ✗ for ( std::size_t jSurfDof(0); jSurfDof < ( std::size_t ) m_curvFePre->numDof(); jSurfDof++) {
997 felInt jDofGlobScalar;
998 ✗ dof().loc2glob( ielSupportDof, jSurfDof, m_iPressure, 0, jDofGlobScalar);
999 ✗ AOApplicationToPetsc(m_ao,1,&jDofGlobScalar);
1000 ✗ m_seqVecs.Get( m_auxiliaryStrings[1] ).getValues(1,&iDofGlobScalar, &l);
1001 ✗ for(int ig=0; ig < m_curvFePre->numQuadraturePoint(); ig++) {
1002 ✗ m_auxiliaryDoubles[0] += l*r * m_curvFePre->phi[ig](iSurfDof) * m_curvFePre->phi[ig](jSurfDof) * m_curvFePre->weightMeas(ig);
1003 }
1004 }
1005 }
1006 }
1007
1008 196717 void LinearProblemNSSimplifiedFSI::normalVelocityScalarComputer(felInt ielSupportDof) {
1009 double value,valueAus;
1010
2/2
✓ Branch 1 taken 586871 times.
✓ Branch 2 taken 196717 times.
 783588 for ( std::size_t iSurfDof(0); iSurfDof < ( std::size_t ) m_curvFeVel->numDof(); iSurfDof++) {
1011 586871 value=0.0;
1012 felInt iDofGlobScalar;
1013
1/2
✓ Branch 2 taken 586871 times.
✗ Branch 3 not taken.
 586871 dof().loc2glob( ielSupportDof, iSurfDof, m_iPressure, 0, iDofGlobScalar);
1014
1/2
✓ Branch 1 taken 586871 times.
✗ Branch 2 not taken.
 586871 AOApplicationToPetsc(m_ao,1,&iDofGlobScalar);
1015
2/2
✓ Branch 1 taken 1754053 times.
✓ Branch 2 taken 586871 times.
 2340924 for ( std::size_t iCoor(0); iCoor < ( std::size_t ) m_curvFeVel->numCoor(); iCoor++) {
1016 felInt iDofGlob;
1017
1/2
✓ Branch 2 taken 1754053 times.
✗ Branch 3 not taken.
 1754053 dof().loc2glob( ielSupportDof, iSurfDof, m_iVelocity, iCoor, iDofGlob);
1018
1/2
✓ Branch 1 taken 1754053 times.
✗ Branch 2 not taken.
 1754053 AOApplicationToPetsc(m_ao,1,&iDofGlob);
1019
3/6
✓ Branch 2 taken 1754053 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1754053 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 1754053 times.
✗ Branch 9 not taken.
 1754053 m_seqVecs.Get("normalVelocity").getValues(1,&iDofGlob,&valueAus);
1020 1754053 value += valueAus;
1021 }
1022
3/6
✓ Branch 2 taken 586871 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 586871 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 586871 times.
✗ Branch 9 not taken.
 586871 m_vecs.Get("normalVelocity").setValue(iDofGlobScalar, value, INSERT_VALUES);
1023 }
1024 196717 }
1025
1026 /*! \brief A function of the loop to sum up the fake components of the displacement and to compute its maximum and minimum.
1027 */
1028 193677 void LinearProblemNSSimplifiedFSI::scalarComputer(felInt ielSupportDof) {
1029 double value,valueAus;
1030
2/2
✓ Branch 1 taken 580791 times.
✓ Branch 2 taken 193677 times.
 774468 for ( std::size_t iSurfDof(0); iSurfDof < ( std::size_t ) m_curvFeVel->numDof(); iSurfDof++) {
1031 580791 value=0.0;
1032 felInt iDofGlobScalar;
1033
1/2
✓ Branch 2 taken 580791 times.
✗ Branch 3 not taken.
 580791 dof().loc2glob( ielSupportDof, iSurfDof, m_iPressure, 0, iDofGlobScalar);
1034
1/2
✓ Branch 1 taken 580791 times.
✗ Branch 2 not taken.
 580791 AOApplicationToPetsc(m_ao,1,&iDofGlobScalar);
1035
2/2
✓ Branch 1 taken 1741893 times.
✓ Branch 2 taken 580791 times.
 2322684 for ( std::size_t iCoor(0); iCoor < ( std::size_t ) m_curvFeVel->numCoor(); iCoor++) {
1036 felInt iDofGlob;
1037
1/2
✓ Branch 2 taken 1741893 times.
✗ Branch 3 not taken.
 1741893 dof().loc2glob( ielSupportDof, iSurfDof, m_iVelocity, iCoor, iDofGlob);
1038
1/2
✓ Branch 1 taken 1741893 times.
✗ Branch 2 not taken.
 1741893 AOApplicationToPetsc(m_ao,1,&iDofGlob);
1039
3/6
✓ Branch 2 taken 1741893 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1741893 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 1741893 times.
✗ Branch 9 not taken.
 1741893 m_seqVecs.Get("displacement").getValues(1,&iDofGlob,&valueAus);
1040 1741893 value += valueAus;
1041 }
1042
2/2
✓ Branch 1 taken 751 times.
✓ Branch 2 taken 580040 times.
 580791 if( value > m_auxiliaryDoubles[0] )
1043 751 m_auxiliaryDoubles[0] = value;
1044
2/2
✓ Branch 1 taken 1077 times.
✓ Branch 2 taken 579714 times.
 580791 if( value < m_auxiliaryDoubles[1] )
1045 1077 m_auxiliaryDoubles[1] = value;
1046
3/6
✓ Branch 2 taken 580791 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 580791 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 580791 times.
✗ Branch 9 not taken.
 580791 m_vecs.Get("displacement").setValue(iDofGlobScalar, value, INSERT_VALUES);
1047 }
1048 193677 }
1049 void
1050 152 LinearProblemNSSimplifiedFSI::readIOPInterface(std::map<int,std::map<int,int> > mapIOP, PetscVector& seqCurrentIOP, std::map<int,int> labelMap, felInt rank) {
1051
2/4
✓ Branch 1 taken 152 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 152 times.
 152 if ( FelisceParam::verbose() > 0 ) {
1052 ✗ std::stringstream out;
1053 ✗ out<<"reading data at the interface...";
1054 ✗ PetscPrintf(MpiInfo::petscComm(), "%s",out.str().c_str());
1055 }
1056 152 felInt iDofGlobScalarPressure(0);
1057 152 felInt iDofGlobScalarIOP(0);
1058 152 felInt iDofGlobScalarPressureAus(0);
1059 152 double valueAus(0.0);
1060
4/6
✓ Branch 1 taken 152 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 304 times.
✗ Branch 7 not taken.
✓ Branch 10 taken 152 times.
✓ Branch 11 taken 152 times.
 304 for(auto it_labelNumber = FelisceParam::instance().labelIOPMesh.begin(); it_labelNumber != FelisceParam::instance().labelIOPMesh.end(); it_labelNumber++) {
1061
1/2
✓ Branch 2 taken 152 times.
✗ Branch 3 not taken.
 152 felInt labelInThisProblem=labelMap[*it_labelNumber];
1062
2/4
✓ Branch 1 taken 152 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 152 times.
 152 if ( FelisceParam::verbose() > 3 ) {
1063 ✗ std::cout << std::endl;
1064 ✗ std::cout << "label IOP mesh: "<<*it_labelNumber<<" label NS mesh: "<<labelInThisProblem<<std::endl;
1065 }
1066
4/6
✓ Branch 2 taken 152 times.
✗ Branch 3 not taken.
✓ Branch 8 taken 12464 times.
✗ Branch 9 not taken.
✓ Branch 12 taken 12312 times.
✓ Branch 13 taken 152 times.
 12464 for(auto it_suppDofIop = mapIOP[*it_labelNumber].begin(); it_suppDofIop != mapIOP[*it_labelNumber].end(); it_suppDofIop++) {
1067 12312 iDofGlobScalarIOP = it_suppDofIop->second;
1068
2/4
✓ Branch 4 taken 12312 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 12312 times.
✗ Branch 8 not taken.
 12312 iDofGlobScalarPressure = m_dofBD[0/*iBD*/].getPetscDofs( this->dof().getNumGlobComp(m_iUnknownPre,0), labelInThisProblem, it_suppDofIop->first );
1069 12312 iDofGlobScalarPressureAus = iDofGlobScalarPressure;
1070
1/2
✓ Branch 1 taken 12312 times.
✗ Branch 2 not taken.
 12312 AOPetscToApplication(m_ao,/*counter*/1,&iDofGlobScalarPressureAus);
1071
2/2
✓ Branch 1 taken 3078 times.
✓ Branch 2 taken 9234 times.
 12312 if ( m_dofPart[iDofGlobScalarPressureAus] == rank ) {
1072
1/2
✓ Branch 1 taken 3078 times.
✗ Branch 2 not taken.
 3078 seqCurrentIOP.getValues(1,&iDofGlobScalarIOP, &valueAus);
1073
3/6
✓ Branch 2 taken 3078 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 3078 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 3078 times.
✗ Branch 9 not taken.
 3078 m_vecs.Get("currentIop").setValue(iDofGlobScalarPressure, valueAus, INSERT_VALUES);
1074 }
1075 }
1076 }
1077
3/6
✓ Branch 2 taken 152 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 152 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 152 times.
✗ Branch 9 not taken.
 152 m_vecs.Get("currentIop").assembly();
1078
5/10
✓ Branch 2 taken 152 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 152 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 152 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 152 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 152 times.
✗ Branch 16 not taken.
 152 gatherVector(m_vecs.Get("currentIop"), m_seqVecs.Get("currentIop"));
1079
3/6
✓ Branch 2 taken 152 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 152 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 152 times.
✗ Branch 9 not taken.
 152 m_seqVecs.Get("currentIop").assembly();
1080
1081
2/4
✓ Branch 1 taken 152 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 152 times.
 152 if ( FelisceParam::verbose() > 0 ) {
1082 ✗ std::stringstream out;
1083 ✗ out<<"done"<<std::endl;
1084 ✗ PetscPrintf(MpiInfo::petscComm(), "%s",out.str().c_str());
1085 }
1086 152 }
1087
1088 /***********************************************************************************/
1089 /***********************************************************************************/
1090
1091 4 void LinearProblemNSSimplifiedFSI::initFixedPointAcceleration()
1092 {
1093 4 m_accMethod = accelerationType(FelisceParam::instance(this->instanceIndex()).accelerationMethod);
1094 4 m_omegaAcceleration = FelisceParam::instance(this->instanceIndex()).omegaAcceleration;
1095
2/4
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
 4 m_seqVecs.Init("x1");
1096
1/3
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 4 switch ( m_accMethod ) {
1097 4 case FixedPoint:
1098 case Relaxation:
1099 4 break;
1100 ✗ case Aitken:
1101 case IronsTuck:
1102 ✗ m_seqVecs.Init("x0");
1103 ✗ m_seqVecs.Init("fx0");
1104 ✗ m_seqVecs.Init("fx1");
1105 ✗ break;
1106 }
1107 4 }
1108
1109 /***********************************************************************************/
1110 /***********************************************************************************/
1111
1112 152 void LinearProblemNSSimplifiedFSI::accelerationPreStep(PetscVector& seqCurrentInput)
1113 {
1114 // pre and post with respect to reading the data at the interface
1115
1116 // method explanation in accelerationPostStep
1117
1/3
✓ Branch 0 taken 152 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 152 switch(m_accMethod) {
1118 152 case FixedPoint:
1119 case Relaxation:
1120 152 break;
1121 ✗ case Aitken:
1122 case IronsTuck:
1123 // in these two cases we have to move forward x0 and x1
1124 ✗ m_seqVecs.Get("x0" ).copyFrom(m_seqVecs.Get("x1" ));
1125 ✗ m_seqVecs.Get("fx0").copyFrom(m_seqVecs.Get("fx1"));
1126 ✗ break;
1127 }
1128 // in any case we need to the save the current value of iop which was the last point
1129 // where we have have evaluated the functional F.
1130
3/6
✓ Branch 2 taken 152 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 152 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 152 times.
✗ Branch 9 not taken.
 152 m_seqVecs.Get("x1").copyFrom(seqCurrentInput);
1131 // The following evaluation of the system will
1132 // be saved in fx1 in the function acceleration
1133 152 }
1134
1135 /***********************************************************************************/
1136 /***********************************************************************************/
1137
1138 152 void LinearProblemNSSimplifiedFSI::accelerationPostStep( PetscVector& seqCurrentInput, int nbOfCurrentIteration )
1139 {
1140 // We are solving a fixed point problem.
1141
1142 // here an example for IOPcouplModel:
1143 // the input to the coupled system is currentiop
1144 // then you solve NS you compute the displacement
1145 // than you solve PF given the displacement and the current iop
1146 // the solution of PF at the boarder is the outcome of function.
1147 // P_interface = F ( currentIop )
1148 // an evaluation of F implies the solution of a NS pb + a PF pb.
1149 // P_interface will be modified by this function and reused for a new iteration.
1150
1151 // Four (4) different methods:
1152 // 0. fixed point iterations
1153 // default method: we simly iterate F over the same value.
1154 // new x1=fx1;
1155 // 1. relaxation:
1156 // new x1= w * fx1 + (1 - w) x1
1157 // w=1 ==> back to fixed point iterations
1158 // 2. Aitken acceleration
1159 // new x1= w * fx1 + (1 - w) x1, but w is a funcion of ( x0, fx0, x1, fx1 )
1160 // 3. Irons-Tuck: the same as Aitken, should be more "stable"
1161
1/4
✓ Branch 0 taken 152 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 152 switch(m_accMethod) {
1162 152 case FixedPoint:
1163 case Relaxation:
1164 // there is no need to save explicitly fx1 as done in aitken and irons tuck
1165 // since it is stored in currentIop and it will be overridden with the correct value.
1166 152 break;
1167 ✗ case Aitken:
1168 {
1169 // at this point we just read the results from the PF problem
1170 // this has been stored in the petscvec currentIop, since in a
1171 // fixed point approach it is already the new input (new x1 = currentIop = fx1).
1172 // in all the other approaches this is just the lates outcome
1173 // of the F.
1174 ✗ m_seqVecs.Get("fx1").copyFrom( seqCurrentInput );
1175
1176 ✗ if ( nbOfCurrentIteration == 1 /*first iteration: w=1*/ ) {
1177 ✗ m_omegaAcceleration = 1;
1178 } else {
1179 // we compute the value for the omega parameter
1180 ✗ m_seqVecs.Init("dx"); // map element initialization
1181 ✗ m_seqVecs.Get("dx").duplicateFrom(m_seqVecs.Get("x1")); // structure of dx copied from x1
1182 ✗ m_seqVecs.Get("dx").copyFrom(m_seqVecs.Get("x1")); // dx = x1
1183 ✗ m_seqVecs.Get("dx").axpy(-1.,m_seqVecs.Get("x0")); // dx = x1 - x0
1184
1185 ✗ m_seqVecs.Init("tmp"); // map element initialization
1186 ✗ m_seqVecs.Get("tmp").duplicateFrom(m_seqVecs.Get("fx0")); // structure of tmp copied from fx0
1187 ✗ m_seqVecs.Get("tmp").copyFrom(m_seqVecs.Get("fx0")); // tmp = fx0
1188 ✗ m_seqVecs.Get("tmp").axpy(-1.,m_seqVecs.Get("fx1")); // tmp = fx0 - fx1
1189 ✗ m_seqVecs.Get("tmp").axpy(1.,m_seqVecs.Get("dx")); // tmp = fx0 - fx1 + x1 - x0
1190
1191 ✗ const double omega = computeL2ScalarProduct("dx" ,"tmp");
1192 ✗ const double xxnorm = computeL2ScalarProduct("tmp","tmp");
1193 ✗ if ( Tools::equal( xxnorm, 0.0 ) )
1194 ✗ m_omegaAcceleration = 1;
1195 else
1196 ✗ m_omegaAcceleration = omega / xxnorm;
1197
1198 // we destroy those auxiliary vectors
1199 ✗ m_seqVecs.Get("tmp").destroy();
1200 ✗ m_seqVecs.erase("tmp");
1201 ✗ m_seqVecs.Get("dx").destroy();
1202 ✗ m_seqVecs.erase("dx");
1203
1204 }
1205 }
1206 ✗ break;
1207 ✗ case IronsTuck:
1208 {
1209 ✗ m_seqVecs.Get("fx1").copyFrom(seqCurrentInput);
1210
1211 ✗ if ( nbOfCurrentIteration == 1 /*first iteration: w=1*/ ) {
1212 ✗ m_omegaAcceleration = 1;
1213 } else {
1214
1215 // we compute the value for the omega parameter
1216 ✗ m_seqVecs.Init("d1"); // map element initialization
1217 ✗ m_seqVecs.Get("d1").duplicateFrom(m_seqVecs.Get("x1")); // structure of d1 copied from x1
1218 ✗ m_seqVecs.Get("d1").copyFrom(m_seqVecs.Get("x1")); // d1 = x1
1219 ✗ m_seqVecs.Get("d1").axpy(-1.,m_seqVecs.Get("fx1")); // d1 = x1 - fx1
1220
1221 ✗ m_seqVecs.Init("tmp"); // map element initialization
1222 ✗ m_seqVecs.Get("tmp").duplicateFrom(m_seqVecs.Get("x0")); // structure of tmp copied from x0
1223 ✗ m_seqVecs.Get("tmp").copyFrom(m_seqVecs.Get("x0")); // tmp = x0
1224 ✗ m_seqVecs.Get("tmp").axpy(-1.,m_seqVecs.Get("fx0")); // tmp = x0 - fx0
1225 ✗ m_seqVecs.Get("tmp").axpy(-1.,m_seqVecs.Get("d1")); // tmp = x0 - fx0 - (x1 - fx1 )
1226
1227 ✗ const double xxnorm ( computeL2ScalarProduct("tmp","tmp") );
1228 ✗ const double omega ( computeL2ScalarProduct("d1","tmp") );
1229 ✗ if ( Tools::equal( xxnorm, 0.0 ) ) {
1230 ✗ m_omegaAcceleration = 1;
1231 } else {
1232 ✗ m_omegaAcceleration = ( 1 + omega / xxnorm ) * m_omegaAcceleration;
1233 }
1234
1235 // we destroy those auxiliary vectors
1236 ✗ m_seqVecs.Get("tmp").destroy();
1237 ✗ m_seqVecs.Get("d1").destroy();
1238 ✗ m_seqVecs.erase("d1");
1239 ✗ m_seqVecs.erase("tmp");
1240 }
1241 }
1242 ✗ break;
1243 }
1244
1245
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 152 times.
 152 if ( FelisceParam::verbose() > 1 ) {
1246 ✗ std::stringstream msg;
1247 ✗ msg<<"----> using omega = "<<m_omegaAcceleration<<std::endl;
1248 ✗ PetscPrintf(MpiInfo::petscComm(),"%s",msg.str().c_str());
1249 }
1250
1251 //for the fixed point this step is not needed
1252
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 152 times.
 152 if ( m_accMethod != FixedPoint ) {
1253 ✗ std::stringstream msg;
1254 ✗ msg<<" updating fixed point input "<<std::endl;
1255 ✗ PetscPrintf(MpiInfo::petscComm(),"%s",msg.str().c_str());
1256 // update of currentIop, i.e. the new x1, the input for the next evaluation of F
1257 ✗ seqCurrentInput.axpby( 1. - m_omegaAcceleration, m_omegaAcceleration, m_seqVecs.Get("x1") );
1258 }
1259 152 }
1260
1261 /***********************************************************************************/
1262 /***********************************************************************************/
1263
1264 }
1265