GCC Code Coverage Report

Directory: ./
File: Model/FSINitscheXFEMModel.hxx
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 393 624 63.0%
Branches: 238 546 43.6%
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: B. Fabreges and M. A. Fernandez
13 //
14
15 // System includes
16
17 // External includes
18
19 // Project includes
20 #include "Geometry/geometricMeshRegion.hpp"
21
22 namespace felisce {
23 template<class SolidPrb>
24
2/4
✓ Branch 16 taken 2 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 2 times.
✗ Branch 20 not taken.
 2 FSINitscheXFEMModel<SolidPrb>::FSINitscheXFEMModel():Model()
25 {
26
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_name = "Nitsche XFEM FSI Model";
27 2 m_isParFluidVelAllocated = false;
28 2 m_isTimeSchemeFluidExtrapolAllocated = false;
29 2 m_isStructureVelAllocated = false;
30 2 m_isStructureVelOldAllocated = false;
31 2 m_isStructureVelOldOldAllocated = false;
32 2 m_isStructureVelOldStoppingCriteriaAllocated = false;
33 2 m_isStructureVelExtrapolAllocated = false;
34 2 m_isSeqStructureTimeRHSAllocated = false;
35 2 m_isSeqStructureVelExtrapolAllocated = false;
36 2 m_isStructureDispAllocated = false;
37 2 m_isStructureDispOldAllocated = false;
38 2 m_isStopVecAllocated = false;
39 2 m_isStructureNitscheRHSTermsOldAllocated = false;
40 2 m_isStructureTimeSchemeRHSTermsAllocated = false;
41 2 m_isStructureExtrapolatedNitscheTermsRHSAllocated = false;
42
43 2 m_elemFieldTimeSchemeType = NULL;
44 2 m_timeSchemeType = 3;
45 2 m_subIteration = 1;
46 2 m_stoppingCriteriaTol = 1e-3;
47 2 m_timeSchemeExtrapolOrder = 0;
48 2 m_indexTime = 0;
49 2 m_countToInitOldObjects = 0;
50 2 m_reinitializeFluidBdf = false;
51 2 m_feStruc = NULL;
52 2 m_feVel = NULL;
53 2 m_fePre = NULL;
54
55 2 m_RNExtrapolOrder = 0;
56 2 m_numStrucComp = 0;
57 2 }
58
59 template<class SolidPrb>
60 4 FSINitscheXFEMModel<SolidPrb>::~FSINitscheXFEMModel()
61 {
62 4 if(m_isParFluidVelAllocated) {
63
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
 8 for(std::size_t i=0; i<m_parFluidVel.size(); ++i)
64 4 m_parFluidVel[i].destroy();
65 4 m_parFluidVel.clear();
66 }
67
68
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if(m_isStructureVelAllocated)
69 4 m_structureVel.destroy();
70
71
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if(m_isStructureVelOldAllocated)
72 4 m_structureVelOld.destroy();
73
74
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if(m_isStructureVelOldOldAllocated)
75 4 m_structureVelOldOld.destroy();
76
77
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 if(m_isStructureVelOldStoppingCriteriaAllocated)
78 ✗ m_structureVelOldStoppingCriteria.destroy();
79
80
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if(m_isStructureVelExtrapolAllocated)
81 4 m_structureVelExtrapol.destroy();
82
83
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if(m_isSeqStructureTimeRHSAllocated)
84 4 m_seqStructureTimeRHS.destroy();
85
86
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if(m_isSeqStructureVelExtrapolAllocated)
87 4 m_seqStructureVelExtrapol.destroy();
88
89
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if(m_isStructureDispOldAllocated)
90 4 m_structureDispOld.destroy();
91
92
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 if(m_isStructureNitscheRHSTermsOldAllocated)
93 ✗ m_structureNitscheRHSTermsOld.destroy();
94
95
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 if(m_isStructureTimeSchemeRHSTermsAllocated)
96 ✗ m_structureTimeSchemeRHSTerms.destroy();
97
98
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 if(m_isStructureExtrapolatedNitscheTermsRHSAllocated)
99 ✗ m_structureExtrapolatedNitscheTermsRHS.destroy();
100
101
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 if(m_isStopVecAllocated)
102 ✗ m_stopVec.destroy();
103
104
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if(m_isTimeSchemeFluidExtrapolAllocated) {
105
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 2 times.
 12 for(std::size_t i=0; i<m_timeSchemeFluidExtrapol.size(); ++i)
106 8 m_timeSchemeFluidExtrapol[i].destroy();
107 4 m_timeSchemeFluidExtrapol.clear();
108 }
109
110
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if(m_feStruc != NULL)
111
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 delete m_feStruc;
112
113
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if(m_feVel != NULL)
114
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 delete m_feVel;
115
116
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if(m_fePre != NULL)
117
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 delete m_fePre;
118 }
119
120 /*!
121 * Intersect the meshes and std::set the informations in the two solvers
122 */
123 template<class SolidPrb>
124 2 void FSINitscheXFEMModel<SolidPrb>::initializeDerivedLinearProblem()
125 {
126
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
 2 FEL_CHECK(m_linearProblem.size() == 2, "Need two solver for Nitsche XFEM FSI Model");
127
128 // type conversion of the linear problems
129
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_pLinFluid = dynamic_cast<LinearProblemFSINitscheXFEMFluid*>(m_linearProblem[0]);
130
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_pLinStruc = dynamic_cast<LinearProblemFSINitscheXFEMSolid<SolidPrb>*>(m_linearProblem[1]);
131
132
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 if(FelisceParam::instance().duplicateSupportDof) {
133
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
 2 FEL_CHECK(m_mesh.size() > 1, "Need at least 2 meshes: one for the fluid and one for the structure");
134
135 // build the edges of the fluid mesh
136 2 m_mesh[0]->buildEdges();
137
138 // We intersect the two meshes here (because this function is called before computeDof)
139 // The meshes can be found in the linear problem
140
3/6
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
 2 m_intersectMeshes.setMeshesName("fluidMesh", "strucMesh");
141 2 m_intersectMeshes.setVerbose(FelisceParam::verbose());
142 2 m_intersectMeshes.initAndIntersectMeshes(*m_mesh[0], *m_mesh[1]);
143
144 // give access to the intersection to each linear problems
145 2 m_pLinFluid->setDuplicateSupportObject(&m_intersectMeshes);
146 2 m_pLinStruc->setDuplicateSupportObject(&m_intersectMeshes);
147
148 // give direct access to the other linear problem
149 2 m_pLinFluid->setStrucLinPrb(m_pLinStruc);
150 2 m_pLinStruc->setFluidLinPrb(m_pLinFluid);
151
152 // compute normals
153 2 m_mesh[1]->computeNormalTangent();
154 2 m_mesh[1]->computeElementNormal(m_mesh[1]->listElementNormals());
155 2 m_mesh[1]->computeElementNormal(m_mesh[1]->listElementNormalsInit());
156 }
157 2 }
158
159 /*!
160 * Initialize time scheme
161 */
162 template<class SolidPrb>
163 2 void FSINitscheXFEMModel<SolidPrb>::initializeDerivedModel()
164 {
165 // get the type of time discretization.
166 2 m_elemFieldTimeSchemeType = FelisceParam::instance().elementFieldDynamicValue("TimeSchemeType");
167
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if(m_elemFieldTimeSchemeType != NULL) {
168 2 m_timeSchemeType = (felInt) m_elemFieldTimeSchemeType->constant(Comp1);
169
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if (m_timeSchemeType == 0) {
170 // Robin-Neumann scheme (direct implementation)
171 // Robin-Neumann extrapolation order
172 ✗ m_RNExtrapolOrder = (felInt) FelisceParam::instance().fsiRNscheme;
173 ✗ m_timeSchemeExtrapolOrder = m_RNExtrapolOrder;
174
175 // resize the sequential std::vector in m_pLinFluid and the parallel one here
176 ✗ m_pLinFluid->initRNFluidExtrapol(m_timeSchemeExtrapolOrder);
177 ✗ m_timeSchemeFluidExtrapol.resize(m_timeSchemeExtrapolOrder);
178
179
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 } else if(m_timeSchemeType == 2) {
180 // Robin-Neumann scheme a la Stenberg
181 // Robin-Neumann extrapolation order
182 ✗ m_RNExtrapolOrder = (felInt) FelisceParam::instance().fsiRNscheme;
183 ✗ m_timeSchemeExtrapolOrder = m_RNExtrapolOrder;
184
185 // resize the sequential std::vector in m_pLinFluid and the parallel one here
186 ✗ m_pLinFluid->initRNFluidExtrapol(m_timeSchemeExtrapolOrder);
187 ✗ m_timeSchemeFluidExtrapol.resize(m_timeSchemeExtrapolOrder);
188
189 // resize the sequential std::vector in m_pLinStruc
190 ✗ m_pLinStruc->initStabExplicitFluidExtrapol(m_timeSchemeExtrapolOrder);
191
192
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 } else if(m_timeSchemeType == 3) {
193 // stabilized explicit
194 // number of sub iteration (correction + 1)
195
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 if(FelisceParam::instance().couplingSchemeStoppingCriteria == 0)
196 2 m_subIteration = (felInt) m_elemFieldTimeSchemeType->constant(Comp2);
197 ✗ else if(FelisceParam::instance().couplingSchemeStoppingCriteria == 1) {
198 ✗ m_stoppingCriteriaTol = (double) m_elemFieldTimeSchemeType->constant(Comp2);
199 ✗ m_stopVec.duplicateFrom(m_pLinStruc->solution());
200 ✗ m_isStopVecAllocated = true;
201 }
202
203 // extrapolation order of the initial fluid solution in structure solver
204 2 m_timeSchemeExtrapolOrder = (felInt) m_elemFieldTimeSchemeType->constant(Comp3);
205
206 // resize the sequential std::vector in m_pLinStruc and the parallel one here
207 2 m_pLinStruc->initStabExplicitFluidExtrapol(m_timeSchemeExtrapolOrder);
208 2 m_timeSchemeFluidExtrapol.resize(m_timeSchemeExtrapolOrder);
209
210 ✗ } else if(m_timeSchemeType == 4) {
211 // Robin-Neumann semi-implicit
212 // Robin-Neumann extrapolation order
213 ✗ m_RNExtrapolOrder = (felInt) FelisceParam::instance().fsiRNscheme;
214
215 // number of sub iterations in the (fluid)-(solid/inertia) problem
216 ✗ if(FelisceParam::instance().couplingSchemeStoppingCriteria == 0)
217 ✗ m_subIteration = (felInt) m_elemFieldTimeSchemeType->constant(Comp2);
218 ✗ else if(FelisceParam::instance().couplingSchemeStoppingCriteria == 1) {
219 ✗ m_stoppingCriteriaTol = (double) m_elemFieldTimeSchemeType->constant(Comp2);
220 ✗ m_stopVec.duplicateFrom(m_pLinStruc->solution());
221 ✗ m_isStopVecAllocated = true;
222 }
223
224 // extrapolation order of the initial fluid solution in structure solver
225 ✗ m_timeSchemeExtrapolOrder = (felInt) m_elemFieldTimeSchemeType->constant(Comp3);
226
227 // resize the sequential std::vector in m_pLinStruc and the parallel one here
228 ✗ m_pLinStruc->initStabExplicitFluidExtrapol(m_timeSchemeExtrapolOrder);
229 ✗ m_timeSchemeFluidExtrapol.resize(m_timeSchemeExtrapolOrder);
230 }
231 }
232
233
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
 2 if (FelisceParam::instance().orderBdfNS > 2)
234 ✗ FEL_ERROR("BDF not yet implemented for order greater than 2 with FDFSI.");
235
236
237 // bdf schemes
238 // fluid
239 2 m_bdfFluid.defineOrder(FelisceParam::instance().orderBdfNS);
240 2 m_pLinFluid->initializeTimeScheme(&m_bdfFluid);
241
242
243 // velocity of the structure
244 // parallel vectors
245 2 m_structureVel.duplicateFrom(m_pLinStruc->solution());
246
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_structureVel.set( 0.);
247 2 m_isStructureVelAllocated = true;
248
249 2 m_structureVelOld.duplicateFrom(m_pLinStruc->solution());
250
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_structureVelOld.set( 0.);
251 2 m_isStructureVelOldAllocated = true;
252
253 2 m_structureVelOldOld.duplicateFrom(m_pLinStruc->solution());
254
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_structureVelOldOld.set( 0.);
255 2 m_isStructureVelOldOldAllocated = true;
256
257
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
 2 if(FelisceParam::instance().couplingSchemeStoppingCriteria == 1) {
258 ✗ m_structureVelOldStoppingCriteria.duplicateFrom(m_pLinStruc->solution());
259 ✗ m_structureVelOldStoppingCriteria.set( 0.);
260 ✗ m_isStructureVelOldStoppingCriteriaAllocated = true;
261 }
262
263 2 m_structureVelExtrapol.duplicateFrom(m_pLinStruc->solution());
264
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_structureVelExtrapol.set( 0.);
265 2 m_isStructureVelExtrapolAllocated = true;
266
267 // sequential vectors
268 2 m_seqStructureTimeRHS.duplicateFrom(m_pLinStruc->sequentialSolution());
269
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_seqStructureTimeRHS.set( 0.0);
270 2 m_isSeqStructureTimeRHSAllocated = true;
271
272 2 m_seqStructureVelExtrapol.duplicateFrom(m_pLinStruc->sequentialSolution());
273
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_seqStructureVelExtrapol.set( 0.0);
274 2 m_isSeqStructureVelExtrapolAllocated = true;
275
276
277 // displacement of the structure
278 // parallel vectors
279 2 m_structureDispOld.duplicateFrom(m_pLinStruc->solution());
280
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_structureDispOld.set( 0.);
281 2 m_isStructureDispOldAllocated = true;
282
283
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 2 if(m_timeSchemeType == 4 && m_RNExtrapolOrder == 2) {
284 ✗ m_structureNitscheRHSTermsOld.duplicateFrom(m_pLinStruc->vector());
285 ✗ m_structureNitscheRHSTermsOld.set( 0.);
286 ✗ m_isStructureNitscheRHSTermsOldAllocated = true;
287 }
288
289
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if(m_timeSchemeType == 4) {
290 ✗ m_structureTimeSchemeRHSTerms.duplicateFrom(m_pLinStruc->vector());
291 ✗ m_structureTimeSchemeRHSTerms.set( 0.);
292 ✗ m_isStructureTimeSchemeRHSTermsAllocated = true;
293
294 ✗ m_structureExtrapolatedNitscheTermsRHS.duplicateFrom(m_pLinStruc->vector());
295 ✗ m_structureExtrapolatedNitscheTermsRHS.set( 0.);
296 ✗ m_isStructureExtrapolatedNitscheTermsRHSAllocated = true;
297 }
298
299
300 // Finite element for the linear problems
301 // We assume that there is only one type of volume element for the structure
302 2 const std::vector<GeometricMeshRegion::ElementType>& bagElementTypeDomainStruc = m_pLinStruc->meshLocal()->bagElementTypeDomain();
303
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
 2 FEL_ASSERT(bagElementTypeDomainStruc.size() == 1);
304
305 const RefElement *refEleStruc;
306 const GeoElement *geoEleStruc;
307
308 2 geoEleStruc = GeometricMeshRegion::eltEnumToFelNameGeoEle[bagElementTypeDomainStruc[0]].second;
309
310 2 int idDisplacement = m_pLinStruc->listVariable().getVariableIdList(displacement);
311 2 int typeOfFiniteElement = m_pLinStruc->listVariable()[idDisplacement].finiteElementType();
312
313
1/2
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 refEleStruc = geoEleStruc->defineFiniteEle(bagElementTypeDomainStruc[0], typeOfFiniteElement, *m_pLinStruc->meshLocal());
314
1/2
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
 2 m_feStruc = new CurvilinearFiniteElement(*refEleStruc, *geoEleStruc, m_pLinStruc->listVariable()[idDisplacement].degreeOfExactness());
315
316
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 FEL_ASSERT(m_feStruc);
317 2 m_pLinFluid->setStrucFiniteElement(m_feStruc);
318
319 // We assume that there is only one type of volume element for the fluid
320 2 const std::vector<GeometricMeshRegion::ElementType>& bagElementTypeDomainFluid = m_pLinFluid->meshLocal()->bagElementTypeDomain();
321
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
 2 FEL_ASSERT(bagElementTypeDomainFluid.size() == 1);
322
323 const RefElement *refEleVel;
324 const GeoElement *geoEleFluid;
325 const RefElement *refElePre;
326
327 2 geoEleFluid = GeometricMeshRegion::eltEnumToFelNameGeoEle[bagElementTypeDomainFluid[0]].second;
328
329 2 int idVel = m_pLinFluid->listVariable().getVariableIdList(velocity);
330 2 int idPre = m_pLinFluid->listVariable().getVariableIdList(pressure);
331 2 int typeOfFiniteElementVel = m_pLinFluid->listVariable()[idVel].finiteElementType();
332 2 int typeOfFiniteElementPre = m_pLinFluid->listVariable()[idPre].finiteElementType();
333
334
1/2
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 refEleVel = geoEleFluid->defineFiniteEle(bagElementTypeDomainFluid[0], typeOfFiniteElementVel, *m_pLinStruc->meshLocal());
335
1/2
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
 2 m_feVel = new CurrentFiniteElement(*refEleVel, *geoEleFluid, m_pLinFluid->listVariable()[idVel].degreeOfExactness());
336
337
1/2
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 refElePre = geoEleFluid->defineFiniteEle(bagElementTypeDomainFluid[0], typeOfFiniteElementPre, *m_pLinStruc->meshLocal());
338
1/2
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
 2 m_fePre = new CurrentFiniteElement(*refElePre, *geoEleFluid, m_pLinFluid->listVariable()[idPre].degreeOfExactness());
339
340
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 FEL_ASSERT(m_feVel);
341
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 FEL_ASSERT(m_fePre);
342 2 m_pLinStruc->setFluidFiniteElement(m_feVel, m_fePre);
343
344
345 // initialize the original bdf for the fluid
346 2 m_pLinFluid->initOriginalBdf(m_timeSchemeType, m_timeSchemeExtrapolOrder);
347
348 // Get the number of components in the displacement
349 2 felInt iUnknown = m_pLinStruc->listUnknown().getUnknownIdList(displacement);
350 2 felInt idVar = m_pLinStruc->listUnknown().idVariable(iUnknown);
351 2 m_numStrucComp = m_pLinStruc->listVariable()[idVar].numComponent();
352 2 }
353
354 template<class SolidPrb>
355 36 void FSINitscheXFEMModel<SolidPrb>::preAssemblingMatrixRHS(std::size_t iProblem)
356 {
357
2/2
✓ Branch 0 taken 18 times.
✓ Branch 1 taken 18 times.
 36 if(iProblem == 0) {
358 // Allocate the extrapolation of the velocity if necessary
359
1/2
✓ Branch 0 taken 18 times.
✗ Branch 1 not taken.
 18 if(!m_isParFluidVelAllocated) {
360 // There are orderBdfNS term in the velocity extrapolation
361 18 m_parFluidVel.resize(FelisceParam::instance().orderBdfNS);
362
363 // If we use a projection of the old solution, all the terms have the same size
364
5/6
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 16 times.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 16 times.
 18 if(FelisceParam::instance().useProjectionForNitscheXFEM || m_fstransient->iteration == 0) {
365
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
 4 for(std::size_t i=0; i<m_parFluidVel.size(); ++i)
366 2 m_parFluidVel[i].duplicateFrom(m_pLinFluid->vector());
367 } else {
368
2/2
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 16 times.
 32 for(std::size_t i=0; i<m_parFluidVel.size(); ++i)
369 16 m_parFluidVel[i].duplicateFrom(m_pLinFluid->solutionOld(i));
370 }
371
372 // std::set the extrapolation term to zero
373
2/2
✓ Branch 1 taken 18 times.
✓ Branch 2 taken 18 times.
 36 for(std::size_t i=0; i<m_parFluidVel.size(); ++i)
374
1/2
✓ Branch 2 taken 18 times.
✗ Branch 3 not taken.
 18 m_parFluidVel[i].set( 0.);
375
376 18 m_isParFluidVelAllocated = true;
377 }
378
379 // initialize fluid extrapolation in fluid linear problem
380
2/4
✓ Branch 0 taken 18 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 18 times.
 18 if((m_timeSchemeType == 0) || (m_timeSchemeType == 2)) {
381 ✗ if(!m_isTimeSchemeFluidExtrapolAllocated) {
382 ✗ if(FelisceParam::instance().useProjectionForNitscheXFEM || m_fstransient->iteration == 0) {
383 ✗ for(std::size_t i=0; i<m_timeSchemeFluidExtrapol.size(); ++i)
384 ✗ m_timeSchemeFluidExtrapol[i].duplicateFrom(m_pLinFluid->solution());
385 } else {
386 ✗ for(std::size_t i=0; i<m_timeSchemeFluidExtrapol.size(); ++i)
387 ✗ m_timeSchemeFluidExtrapol[i].duplicateFrom(m_pLinFluid->solutionOld(i));
388 }
389
390 ✗ for(std::size_t i=0; i<m_timeSchemeFluidExtrapol.size(); ++i)
391 ✗ m_timeSchemeFluidExtrapol[i].set(0.);
392
393 ✗ m_pLinFluid->initRNFluidExtrapol(m_timeSchemeExtrapolOrder);
394 ✗ m_pLinFluid->gatherRNFluidExtrapol(m_timeSchemeFluidExtrapol);
395
396 ✗ if(m_timeSchemeType == 0)
397 ✗ m_isTimeSchemeFluidExtrapolAllocated = true;
398 }
399 }
400 }
401
402
2/2
✓ Branch 0 taken 18 times.
✓ Branch 1 taken 18 times.
 36 if(iProblem == 1) {
403 // initialize fluid extrapolation in structure linear problem
404
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 18 if ((m_timeSchemeType == 3) || (m_timeSchemeType == 4)) {
405
1/2
✓ Branch 0 taken 18 times.
✗ Branch 1 not taken.
 18 if(!m_isTimeSchemeFluidExtrapolAllocated) {
406
5/6
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 16 times.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 16 times.
 18 if(FelisceParam::instance().useProjectionForNitscheXFEM || m_fstransient->iteration == 0) {
407
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 2 times.
 6 for(std::size_t i=0; i<m_timeSchemeFluidExtrapol.size(); ++i)
408 4 m_timeSchemeFluidExtrapol[i].duplicateFrom(m_pLinFluid->solution());
409 } else {
410
2/2
✓ Branch 1 taken 32 times.
✓ Branch 2 taken 16 times.
 48 for(std::size_t i=0; i<m_timeSchemeFluidExtrapol.size(); ++i)
411 32 m_timeSchemeFluidExtrapol[i].duplicateFrom(m_pLinFluid->solutionOld(i));
412 }
413
414
2/2
✓ Branch 1 taken 36 times.
✓ Branch 2 taken 18 times.
 54 for(std::size_t i=0; i<m_timeSchemeFluidExtrapol.size(); ++i)
415
1/2
✓ Branch 2 taken 36 times.
✗ Branch 3 not taken.
 36 m_timeSchemeFluidExtrapol[i].set(0.);
416
417 18 m_pLinStruc->initStabExplicitFluidExtrapol(m_timeSchemeExtrapolOrder);
418 18 m_pLinStruc->gatherStabExplicitFluidExtrapol(m_timeSchemeFluidExtrapol);
419 18 m_isTimeSchemeFluidExtrapolAllocated = true;
420 }
421
0/2
✗ Branch 0 not taken.
✗ Branch 1 not taken.
 18 } else if (m_timeSchemeType == 2) {
422 ✗ if(!m_isTimeSchemeFluidExtrapolAllocated) {
423 ✗ m_pLinStruc->initStabExplicitFluidExtrapol(m_timeSchemeExtrapolOrder);
424 ✗ m_pLinStruc->gatherStabExplicitFluidExtrapol(m_timeSchemeFluidExtrapol);
425 ✗ m_isTimeSchemeFluidExtrapolAllocated = true;
426 }
427 }
428 }
429 36 }
430
431 template<class SolidPrb>
432 ✗ void FSINitscheXFEMModel<SolidPrb>::postAssemblingMatrixRHS(std::size_t iProblem)
433 {
434 ✗ if(iProblem == 0) {
435 // m_Matrix = m_A + m_Matrix (add static matrix to the dynamic matrix to build
436 // complete matrix of the system
437 ✗ m_linearProblem[iProblem]->addMatrixRHS();
438 }
439 }
440
441 template<class SolidPrb>
442 2 void FSINitscheXFEMModel<SolidPrb>::setExternalVec()
443 {
444 // externalVec(0) in m_pLinFluid
445 2 m_pLinFluid->pushBackExternalVec(m_seqStructureVelExtrapol);
446
447 // externalVec(0) in m_pLinSolid
448 2 m_pLinStruc->pushBackExternalVec(m_seqStructureTimeRHS);
449 2 }
450
451 template<class SolidPrb>
452 20 void FSINitscheXFEMModel<SolidPrb>::forward()
453 {
454 // check solution on crack solid to determine crack spring constant -- if it exists in user
455 20 m_pLinStruc->userComputeCrackStress();
456
2/2
✓ Branch 1 taken 11 times.
✓ Branch 2 taken 9 times.
 20 if(m_intersectMeshes.hasIntersectionChange()) {
457 // reallocate everything for the fluid
458 // supportDofMeshes, duplication, and dofs
459 11 m_pLinFluid->computeDof(MpiInfo::numProc(), MpiInfo::rankProc());
460
461 // new pattern
462 11 m_pLinFluid->userChangePattern(MpiInfo::numProc(), MpiInfo::rankProc());
463
464 // cut the mesh
465 11 m_pLinFluid->cutMesh();
466
467 // allocate matrices and vectors
468 11 m_pLinFluid->allocateMatrix();
469
470 // determine dof for the boundary conditions
471 11 m_pLinFluid->determineDofAssociateToLabel();
472 11 m_pLinFluid->finalizeEssBCConstantInT();
473
474 // pre assembling to reallocate the extrapolations
475
2/2
✓ Branch 1 taken 22 times.
✓ Branch 2 taken 11 times.
 33 for (std::size_t ipb = 0; ipb < m_linearProblem.size(); ipb++) {
476 22 preAssemblingMatrixRHS(ipb);
477 }
478
479 // re build the solver (not possible to change the matrix size in the ksp without destroying it ?)
480 11 m_pLinFluid->buildSolver();
481
482 // Copy all the vectors of the fluid from the original support dof
483 11 m_pLinFluid->copyOriginalSolToCurrentSol();
484 }
485
486
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 18 times.
 20 if(m_fstransient->iteration == 0) {
487 // initialize the old objects
488 2 m_pLinFluid->initOldObjects(m_timeSchemeExtrapolOrder);
489
490 // Write solution for postprocessing (if required)
491 2 writeSolutionAndMeshes();
492 }
493
494 // Advance time step.
495 20 updateTime();
496
497 // Print time information
498 20 printNewTimeIterationBanner();
499
500 ///////////////////////
501 // Update Bdf scheme //
502 ///////////////////////
503
504
8/8
✓ Branch 1 taken 9 times.
✓ Branch 2 taken 11 times.
✓ Branch 3 taken 4 times.
✓ Branch 4 taken 5 times.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 2 times.
✓ Branch 8 taken 18 times.
✓ Branch 9 taken 2 times.
 20 if(m_intersectMeshes.hasIntersectionChange() || m_reinitializeFluidBdf || m_fstransient->iteration == 1) {
505 // initialize the bdf for the fluid
506
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 18 times.
 18 if(FelisceParam::instance().useProjectionForNitscheXFEM) {
507 ✗ m_pLinFluid->copyOriginalBdfToCurrentBdf(m_timeSchemeType, m_timeSchemeExtrapolOrder);
508 } else {
509
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 18 felInt order = std::max(m_timeSchemeExtrapolOrder, FelisceParam::instance().orderBdfNS);
510
1/2
✓ Branch 1 taken 18 times.
✗ Branch 2 not taken.
 18 PetscVector zero = PetscVector::null();
511
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 18 times.
 18 if(order == 1) {
512 ✗ m_bdfFluid.reinitialize(FelisceParam::instance().orderBdfNS, m_pLinFluid->solutionOld(0), zero, zero);
513
1/2
✓ Branch 0 taken 18 times.
✗ Branch 1 not taken.
 18 } else if(order == 2) {
514
2/4
✓ Branch 3 taken 18 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 18 times.
✗ Branch 7 not taken.
 18 m_bdfFluid.reinitialize(FelisceParam::instance().orderBdfNS, m_pLinFluid->solutionOld(1), m_pLinFluid->solutionOld(0), zero);
515 } else {
516 ✗ FEL_ERROR("This bdf order is not implemented for this problem");
517 }
518 18 }
519
520 // Compute the extrapolation of the velocity and initialize the one in the solver
521 18 m_bdfFluid.extrapolateByPart(m_parFluidVel);
522 18 m_pLinFluid->gatherSeqVelExtrapol(m_parFluidVel);
523
524 // Compute the rhs coming from the time scheme and initialize the one in the solver
525 18 m_bdfFluid.computeRHSTimeByPart(m_fstransient->timeStep, m_parFluidVel);
526 18 m_pLinFluid->gatherSeqBdfRHS(m_parFluidVel);
527
528 // std::set flag to false
529 18 m_reinitializeFluidBdf = false;
530 } else {
531 2 m_bdfFluid.update(m_pLinFluid->solution());
532
533 2 m_bdfFluid.extrapolateByPart(m_parFluidVel);
534 2 m_pLinFluid->gatherSeqVelExtrapol(m_parFluidVel);
535
536 2 m_bdfFluid.computeRHSTimeByPart(m_fstransient->timeStep, m_parFluidVel);
537 2 m_pLinFluid->gatherSeqBdfRHS(m_parFluidVel);
538 }
539
540 // structure
541 20 m_structureDispOld.axpby(2., -1., m_pLinStruc->solution());
542 20 m_pLinStruc->gatherVector(m_structureDispOld, m_seqStructureTimeRHS);
543 20 m_structureDispOld.copyFrom(m_pLinStruc->solution());
544 20 m_pLinStruc->gatherSeqLastTimeSol(m_structureDispOld);
545
546
547 // move the mesh if the displacement is a std::vector
548 // if( m_fstransient->iteration > 1 && m_numStrucComp == m_pLinFluid->dimension()) {
549 // m_pLinStruc->meshLocal()->moved() = false;
550 // m_pLinStruc->mesh()->moved() = false;
551 // m_pLinStruc->meshLocal()->moveMesh(m_interfaceDisplacement, 1);
552 // m_pLinStruc->mesh()->moveMesh(m_interfaceDisplacement, 1);
553 // }
554
555 // specific operation of the time schemes
556
1/6
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 20 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 20 switch(m_timeSchemeType) {
557 ✗ case 0: {
558 ✗ computeTimeSchemeExtrapol();
559 ✗ m_pLinFluid->gatherRNFluidExtrapol(m_timeSchemeFluidExtrapol);
560 ✗ m_pLinFluid->setUseCurrentFluidSolution(false);
561 ✗ m_pLinStruc->setUseCurrentFluidSolution(true);
562 ✗ solveProblemFluid();
563 ✗ solveProblemSolid();
564 ✗ break;
565 }
566
567 ✗ case 1: {
568 ✗ solveProblemSolid();
569 ✗ solveProblemFluid();
570 ✗ break;
571 }
572
573 ✗ case 2: {
574 // Check if we need to move the mesh to its initial position to assemble solid terms
575 ✗ bool moveTheMesh = m_fstransient->iteration > 1 && FelisceParam::instance().useDynamicStructure;
576
577 ✗ computeTimeSchemeExtrapol();
578
579 ✗ m_pLinFluid->gatherRNFluidExtrapol(m_timeSchemeFluidExtrapol);
580 ✗ m_pLinFluid->setUseCurrentFluidSolution(false);
581 ✗ solveProblemFluid();
582
583 ✗ m_pLinStruc->setUseCurrentFluidSolution(true);
584
585 // Assembly loop
586 // 1) build the term on the original mesh
587 ✗ if(moveTheMesh) {
588 ✗ m_pLinStruc->meshLocal()->moveMesh(m_interfaceDisplacement, 0);
589 ✗ m_pLinStruc->mesh()->moveMesh(m_interfaceDisplacement, 0);
590 }
591
592 ✗ m_pLinStruc->isAssemblingOfSolidProblemTerms(true);
593 ✗ m_pLinStruc->gatherRNStructExtrapol(m_structureVelExtrapol);
594 ✗ m_pLinStruc->assembleMatrixRHSBD(MpiInfo::rankProc());
595
596 // build the term with the fluid on the moved mesh
597 ✗ if(moveTheMesh) {
598 ✗ m_pLinStruc->meshLocal()->moveMesh(m_interfaceDisplacement, 1.0);
599 ✗ m_pLinStruc->mesh()->moveMesh(m_interfaceDisplacement, 1.0);
600 }
601
602 ✗ m_pLinStruc->isAssemblingOfSolidProblemTerms(false);
603 ✗ m_pLinStruc->initStabExplicitFluidExtrapol(1);
604 ✗ m_pLinStruc->assembleMatrixRHSBD(MpiInfo::rankProc());
605
606 ✗ m_pLinStruc->setUseCurrentFluidSolution(false);
607
608 ✗ m_pLinStruc->initStabExplicitFluidExtrapol(m_timeSchemeExtrapolOrder);
609 ✗ m_pLinStruc->gatherStabExplicitFluidExtrapol(m_timeSchemeFluidExtrapol);
610 ✗ m_pLinStruc->assembleMatrixRHSBD(MpiInfo::rankProc());
611
612 // Apply boundary conditions.
613 ✗ m_pLinStruc->setUserBoundaryCondition();
614
615 // Solve linear system.
616 ✗ PetscPrintf(PETSC_COMM_WORLD, "*** Structure solver ***\n");
617 ✗ m_pLinStruc->solve(MpiInfo::rankProc(), MpiInfo::numProc());
618
619 ////////////////////////////////////////////////////
620 // Compute the vec to link the two linear problem //
621 ////////////////////////////////////////////////////
622
623 ✗ computeTimeSchemeStructureVelExtrapol();
624
625 // gather solution
626 ✗ m_pLinStruc->gatherSolution();
627
628 ✗ break;
629 }
630
631 20 case 3: {
632 // initialize the fluid solution with the extrapolation
633 20 computeTimeSchemeExtrapol();
634 20 m_pLinStruc->gatherStabExplicitFluidExtrapol(m_timeSchemeFluidExtrapol);
635 20 m_pLinFluid->setUseCurrentFluidSolution(false);
636 20 m_pLinStruc->setUseCurrentFluidSolution(false);
637
638 20 felInt it = 0;
639 20 bool atConvergence = false;
640 20 double abstol = 1e-12;
641
2/2
✓ Branch 0 taken 40 times.
✓ Branch 1 taken 20 times.
 60 while(!atConvergence) {
642 // clear the matrices at each sub iteration
643 40 m_pLinFluid->clearMatrixRHS(FlagMatrixRHS::matrix_and_rhs);
644 40 m_pLinStruc->clearMatrixRHS(FlagMatrixRHS::matrix_and_rhs);
645
646 40 solveProblemSolid();
647 40 solveProblemFluid();
648
649 // use the fluid solution (not the extrapolation) for the next sub iterations
650
2/2
✓ Branch 0 taken 20 times.
✓ Branch 1 taken 20 times.
 40 if(it == 0) {
651 20 m_pLinStruc->initStabExplicitFluidExtrapol(1);
652 20 m_pLinStruc->gatherStabExplicitFluidExtrapol(m_pLinFluid->solution());
653
654 20 m_pLinFluid->setUseCurrentFluidSolution(true);
655 20 m_pLinStruc->setUseCurrentFluidSolution(true);
656 } else {
657 20 m_pLinStruc->setStabExplicitFluidExtrapol(m_pLinFluid->sequentialSolution());
658 }
659
660 // increment the number of iteration
661 40 ++it;
662
663 // compute the stopping criteria
664
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 if(FelisceParam::instance().couplingSchemeStoppingCriteria == 0) {
665 // fix number of iteration
666 40 atConvergence = it < m_subIteration ? false : true;
667 ✗ } else if(FelisceParam::instance().couplingSchemeStoppingCriteria == 1) {
668 // compute relative error on structure
669 ✗ m_stopVec.copyFrom(m_structureVel);
670 ✗ double rhs = std::sqrt(m_pLinStruc->computeL2Norm(m_stopVec));
671 ✗ m_stopVec.axpy(-1., m_structureVelOldStoppingCriteria);
672 ✗ double lhs = std::sqrt(m_pLinStruc->computeL2Norm(m_stopVec));
673
674 ✗ atConvergence = ((lhs < m_stoppingCriteriaTol * rhs) || (lhs < abstol)) && (it > 1) ? true : false;
675 ✗ PetscPrintf(PETSC_COMM_WORLD, "Relative Error after iteration %d : %12.5e --- %12.5e\n\n", it, lhs, m_stoppingCriteriaTol*rhs);
676 }
677 }
678 20 break;
679 }
680
681 ✗ case 4: {
682 ////////////////////////////////////////////////////
683 // fluid-(solid-inertia) problem //
684 ////////////////////////////////////////////////////
685
686 // gather the previous solid-velocity and solid-displacement for the (fluid)-(solid-inertia) problem
687 ✗ m_pLinStruc->gatherRNStructExtrapol(m_structureVelExtrapol);
688
689 // initialize the fluid solution with the extrapolation
690 ✗ computeTimeSchemeExtrapol();
691 ✗ m_pLinStruc->gatherStabExplicitFluidExtrapol(m_timeSchemeFluidExtrapol);
692
693 ✗ m_pLinStruc->setUseCurrentFluidSolution(false);
694 ✗ m_pLinFluid->setUseCurrentFluidSolution(false);
695
696 ✗ m_pLinStruc->isStrucOperatorExplicitlyComputed(true);
697
698 ✗ felInt it = 0;
699 ✗ bool atConvergence = false;
700 ✗ double abstol = 1e-12;
701 ✗ while(!atConvergence) {
702 // clear the matrices at each sub iteration
703 ✗ m_pLinStruc->clearMatrixRHS(FlagMatrixRHS::matrix_and_rhs);
704 ✗ m_pLinFluid->clearMatrixRHS(FlagMatrixRHS::matrix_and_rhs);
705
706 ✗ solveProblemSolid();
707 ✗ solveProblemFluid();
708
709 // use the fluid solution (not the extrapolation) for the next sub iterations
710 ✗ if(it == 0) {
711 ✗ m_pLinStruc->initStabExplicitFluidExtrapol(1);
712 ✗ m_pLinStruc->gatherStabExplicitFluidExtrapol(m_pLinFluid->solution());
713
714 ✗ m_pLinFluid->setUseCurrentFluidSolution(true);
715 ✗ m_pLinStruc->setUseCurrentFluidSolution(true);
716 } else {
717 ✗ m_pLinStruc->setStabExplicitFluidExtrapol(m_pLinFluid->sequentialSolution());
718 }
719
720 // increment the number of iteration
721 ✗ ++it;
722
723 // compute the stopping criteria
724 ✗ if(FelisceParam::instance().couplingSchemeStoppingCriteria == 0) {
725 // fix number of iteration
726 ✗ atConvergence = it < m_subIteration ? false : true;
727 ✗ } else if(FelisceParam::instance().couplingSchemeStoppingCriteria == 1) {
728 // compute relative error on structure
729 ✗ m_stopVec.copyFrom(m_structureVelExtrapol);
730 ✗ double rhs = std::sqrt(m_pLinStruc->computeL2Norm(m_stopVec));
731 ✗ m_stopVec.axpy(-1., m_structureVelOldStoppingCriteria);
732 ✗ double lhs = std::sqrt(m_pLinStruc->computeL2Norm(m_stopVec));
733
734 ✗ atConvergence = ((lhs < m_stoppingCriteriaTol * rhs) || (lhs < abstol)) && (it > 1) ? true : false;
735 ✗ PetscPrintf(PETSC_COMM_WORLD, "Relative Error after iteration %d : %12.5e --- %12.5e\n\n", it, lhs, m_stoppingCriteriaTol*rhs);
736 }
737 }
738
739
740 //////////////////////////////////
741 // solve the real solid problem //
742 //////////////////////////////////
743
744 ✗ m_pLinStruc->clearMatrixRHS(FlagMatrixRHS::matrix_and_rhs);
745 ✗ m_pLinStruc->isStrucOperatorExplicitlyComputed(false);
746 ✗ m_pLinStruc->setUseCurrentFluidSolution(true);
747 ✗ m_pLinStruc->initStabExplicitFluidExtrapol(1);
748 ✗ m_pLinStruc->gatherStabExplicitFluidExtrapol(m_pLinFluid->solution());
749
750 ✗ solveProblemSolid();
751
752 ✗ break;
753 }
754
755 ✗ default: {
756 ✗ solveProblemFluid();
757 ✗ solveProblemSolid();
758 ✗ break;
759 }
760 }
761
762 // post process //
763 // Write solution (to get the real solution define on the duplicated support dof, not on the original one)
764 20 writeSolutionAndMeshes();
765
766 // update the position of the interface and recompute the intersection
767
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 if(FelisceParam::instance().useDynamicStructure) {
768
769 20 updateInterfacePosition( m_pLinStruc->sequentialSolution() );
770 20 m_intersectMeshes.updateInterfacePosition( m_interfaceDisplacement );
771
772
1/2
✓ Branch 3 taken 20 times.
✗ Branch 4 not taken.
 20 m_mesh[1]->moveMesh(m_interfaceDisplacement, 1.0);
773 20 m_mesh[1]->computeElementNormal(m_mesh[1]->listElementNormals());
774
1/2
✓ Branch 3 taken 20 times.
✗ Branch 4 not taken.
 20 m_mesh[1]->moveMesh(m_interfaceDisplacement, 0.);
775
776 20 m_intersectMeshes.setIteration(m_fstransient->iteration);
777 20 m_intersectMeshes.intersectMeshes(*m_mesh[0], *m_mesh[1]);
778 }
779
780 // save the solution with the original support Dof and delete the matrices, vectors and mapping of the fluid
781
2/2
✓ Branch 1 taken 13 times.
✓ Branch 2 taken 7 times.
 20 if(m_intersectMeshes.hasIntersectionChange()) {
782 // copy for the fluid
783 13 m_pLinFluid->copyCurrentSolToOriginalSol();
784
785 // delete dynamic parts
786
2/2
✓ Branch 1 taken 11 times.
✓ Branch 2 taken 2 times.
 13 if(!hasFinished()) {
787 11 m_pLinFluid->updateOldSolution(1);
788 11 m_pLinFluid->deleteDynamicData();
789 11 m_pLinStruc->deleteDynamicData();
790
791
1/2
✓ Branch 0 taken 11 times.
✗ Branch 1 not taken.
 11 if(m_isTimeSchemeFluidExtrapolAllocated) {
792
2/2
✓ Branch 1 taken 22 times.
✓ Branch 2 taken 11 times.
 33 for(std::size_t i=0; i<m_timeSchemeFluidExtrapol.size(); ++i)
793 22 m_timeSchemeFluidExtrapol[i].destroy();
794
795 11 m_isTimeSchemeFluidExtrapolAllocated = false;
796 }
797
798
1/2
✓ Branch 0 taken 11 times.
✗ Branch 1 not taken.
 11 if(m_isParFluidVelAllocated) {
799
2/2
✓ Branch 1 taken 11 times.
✓ Branch 2 taken 11 times.
 22 for(std::size_t i=0; i<m_parFluidVel.size(); ++i)
800 11 m_parFluidVel[i].destroy();
801
802 11 m_isParFluidVelAllocated = false;
803 }
804 }
805
806 // std::set to true only to remember that the intersection has changed after this time step
807 // If the intersection is not changing at the next time step, we need to reinitialize the old objects
808 13 m_countToInitOldObjects = std::max(m_timeSchemeExtrapolOrder, FelisceParam::instance().orderBdfNS);
809 } else {
810 7 m_pLinFluid->updateOldSolution(m_countToInitOldObjects);
811
812
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 2 times.
 7 if(m_countToInitOldObjects > 0) {
813
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
 5 if(m_isTimeSchemeFluidExtrapolAllocated) {
814
2/2
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 5 times.
 15 for(std::size_t i=0; i<m_timeSchemeFluidExtrapol.size(); ++i)
815 10 m_timeSchemeFluidExtrapol[i].destroy();
816
817 5 m_isTimeSchemeFluidExtrapolAllocated = false;
818 }
819
820
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
 5 if(m_isParFluidVelAllocated) {
821
2/2
✓ Branch 1 taken 5 times.
✓ Branch 2 taken 5 times.
 10 for(std::size_t i=0; i<m_parFluidVel.size(); ++i)
822 5 m_parFluidVel[i].destroy();
823
824 5 m_isParFluidVelAllocated = false;
825 }
826
827 // pre assembling to reallocate the extrapolations
828
2/2
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 5 times.
 15 for (std::size_t ipb = 0; ipb < m_linearProblem.size(); ipb++) {
829 10 preAssemblingMatrixRHS(ipb);
830 }
831
832 5 m_countToInitOldObjects -= 1;
833 5 m_reinitializeFluidBdf = true;
834 } else {
835
2/6
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 2 if((m_timeSchemeType == 2) || (m_timeSchemeType == 3) || (m_timeSchemeType == 4)) {
836 // reset the right size for the extrapolation of the velocity
837 2 m_pLinStruc->initStabExplicitFluidExtrapol(m_timeSchemeExtrapolOrder);
838 }
839 }
840 }
841 20 }
842
843 template<class SolidPrb>
844 40 void FSINitscheXFEMModel<SolidPrb>::solveProblemFluid()
845 {
846 // Assembly loop for front points
847 40 m_pLinFluid->assembleMatrixFrontPoints();
848
849 // Assembly loop for the ghost penalty stabilization
850
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 if(FelisceParam::instance().useGhostPenalty)
851 40 m_pLinFluid->assembleMatrixGhostPenalty();
852
853 // Assembly loop for the face oriented stabilization
854
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 if(FelisceParam::instance().NSStabType == 1)
855 40 m_pLinFluid->assembleMatrixFaceOriented();
856
857 // Main assembly loop and assembly of the matrix
858 // The pattern is reduced to where values have been std::set till here, that is why this loop is called in last
859 40 m_pLinFluid->assembleMatrixRHS(MpiInfo::rankProc());
860
861 // Specific operations before solving the system
862 // postAssemblingMatrixRHS(0);
863
864 // Apply boundary conditions.
865 // m_pLinFluid->finalizeEssBCConstantInT();
866 40 m_pLinFluid->finalizeEssBCTransient();
867 40 m_pLinFluid->applyBC(FelisceParam::instance().essentialBoundaryConditionsMethod, MpiInfo::rankProc());
868
869 // Solve the linear system.
870 40 PetscPrintf(PETSC_COMM_WORLD, "*** Fluid solver ***\n");
871 40 m_pLinFluid->solve(MpiInfo::rankProc(), MpiInfo::numProc());
872
873 // Gather the solution to use it in the solid problem
874 40 m_pLinFluid->gatherSolution();
875 40 }
876
877 template<class SolidPrb>
878 40 void FSINitscheXFEMModel<SolidPrb>::solveProblemSolid()
879 {
880 // Check if we need to move the mesh to its initial position to assemble solid terms
881
3/4
✓ Branch 1 taken 36 times.
✓ Branch 2 taken 4 times.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
 40 bool moveTheMesh = m_fstransient->iteration > 1 && FelisceParam::instance().useDynamicStructure;
882
883 // Assembly loop
884 // 1) build the term on the original mesh
885
2/2
✓ Branch 0 taken 36 times.
✓ Branch 1 taken 4 times.
 40 if(moveTheMesh) {
886
1/2
✓ Branch 3 taken 36 times.
✗ Branch 4 not taken.
 36 m_pLinStruc->meshLocal()->moveMesh(m_interfaceDisplacement, 0);
887
1/2
✓ Branch 3 taken 36 times.
✗ Branch 4 not taken.
 36 m_pLinStruc->mesh()->moveMesh(m_interfaceDisplacement, 0);
888 }
889
890 40 m_pLinStruc->isAssemblingOfSolidProblemTerms(true);
891 40 m_pLinStruc->assembleMatrixRHSBD(MpiInfo::rankProc());
892
893
2/6
✗ Branch 0 not taken.
✓ Branch 1 taken 40 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 40 times.
 40 if(m_timeSchemeType == 4 && !m_pLinStruc->isStrucOperatorExplicitlyComputed()) {
894 ✗ m_structureTimeSchemeRHSTerms.copyFrom(m_pLinStruc->vector());
895 }
896
897 // build the term with the fluid on the moved mesh
898
2/2
✓ Branch 0 taken 36 times.
✓ Branch 1 taken 4 times.
 40 if(moveTheMesh) {
899
1/2
✓ Branch 3 taken 36 times.
✗ Branch 4 not taken.
 36 m_pLinStruc->meshLocal()->moveMesh(m_interfaceDisplacement, 1.0);
900
1/2
✓ Branch 3 taken 36 times.
✗ Branch 4 not taken.
 36 m_pLinStruc->mesh()->moveMesh(m_interfaceDisplacement, 1.0);
901 }
902
903 40 m_pLinStruc->isAssemblingOfSolidProblemTerms(false);
904 40 m_pLinStruc->assembleMatrixRHSBD(MpiInfo::rankProc());
905
906
2/6
✗ Branch 0 not taken.
✓ Branch 1 taken 40 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 40 times.
 40 if(m_timeSchemeType == 4 && !m_pLinStruc->isStrucOperatorExplicitlyComputed()) {
907 // update the vec
908 ✗ computeTimeSchemeStructureExtrapolRHS();
909
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 40 times.
 40 } else if(m_timeSchemeType == 4) {
910 // add the extrapolation to the rhs
911 ✗ m_pLinStruc->vector().axpy(-1., m_structureExtrapolatedNitscheTermsRHS);
912 }
913
914
915 // Apply boundary conditions.
916 40 m_pLinStruc->setUserBoundaryCondition();
917
918 // Solve linear system.
919 40 PetscPrintf(PETSC_COMM_WORLD, "*** Structure solver ***\n");
920
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 40 times.
 40 if(FelisceParam::instance().useODESolve)
921 ✗ m_pLinStruc->solve(MpiInfo::rankProc(), MpiInfo::numProc(), !m_pLinStruc->getUseCurrentFluidSolution());
922 else
923 40 m_linearProblem[1]->solve(MpiInfo::rankProc(), MpiInfo::numProc());
924
925
926 ////////////////////////////////////////////////////
927 // Compute the vec to link the two linear problem //
928 ////////////////////////////////////////////////////
929
930 40 computeTimeSchemeStructureVelExtrapol();
931
932 // gather solution
933 40 m_pLinStruc->gatherSolution();
934 40 }
935
936 template<class SolidPrb>
937 20 void FSINitscheXFEMModel<SolidPrb>::computeTimeSchemeExtrapol()
938 {
939
1/4
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 20 times.
✗ Branch 3 not taken.
 20 switch (m_timeSchemeExtrapolOrder) {
940 ✗ case 0:
941 ✗ break;
942 ✗ case 1:
943 ✗ m_timeSchemeFluidExtrapol[0].copyFrom(m_bdfFluid.sol_n());
944 ✗ break;
945 20 case 2:
946
2/4
✓ Branch 0 taken 20 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 20 times.
 20 if((m_timeSchemeType == 0) || (m_timeSchemeType == 2)) {
947 ✗ m_timeSchemeFluidExtrapol[0].copyFrom(m_bdfFluid.sol_n());
948 ✗ m_timeSchemeFluidExtrapol[1].copyFrom(m_bdfFluid.sol_n_1());
949 } else {
950
1/2
✓ Branch 2 taken 20 times.
✗ Branch 3 not taken.
 20 m_timeSchemeFluidExtrapol[0].set( 0.);
951
1/2
✓ Branch 2 taken 20 times.
✗ Branch 3 not taken.
 20 m_timeSchemeFluidExtrapol[1].set( 0.);
952 20 m_timeSchemeFluidExtrapol[0].axpy( 2.0, m_bdfFluid.sol_n());
953 20 m_timeSchemeFluidExtrapol[1].axpy( -1.0, m_bdfFluid.sol_n_1());
954 }
955 20 break;
956 ✗ default:
957 ✗ FEL_ERROR("unknown time scheme extrapolation order ");
958 ✗ break;
959 }
960 20 }
961
962 template<class SolidPrb>
963 40 void FSINitscheXFEMModel<SolidPrb>::computeTimeSchemeStructureVelExtrapol()
964 {
965
3/8
✓ Branch 0 taken 40 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 40 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 40 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 40 bool isRNScheme = (m_timeSchemeType == 0) || (m_timeSchemeType == 2) || (m_timeSchemeType == 4 && !m_pLinStruc->isStrucOperatorExplicitlyComputed());
966
967
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 40 times.
 40 if(isRNScheme) {
968 ✗ switch(m_RNExtrapolOrder) {
969 ✗ case 0:
970 ✗ break;
971 ✗ case 1:
972 ✗ m_structureVelOld.copyFrom(m_structureVel);
973 ✗ break;
974 ✗ case 2:
975 ✗ m_structureVelOldOld.copyFrom(m_structureVelOld);
976 ✗ m_structureVelOld.copyFrom(m_structureVel);
977 ✗ break;
978 ✗ default:
979 ✗ FEL_ERROR("unknown time scheme extrapolation order");
980 ✗ break;
981 }
982 }
983 else {
984
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 40 times.
 40 if(FelisceParam::instance().couplingSchemeStoppingCriteria == 1) {
985 ✗ if(m_timeSchemeType == 3)
986 ✗ m_structureVelOldStoppingCriteria.copyFrom(m_structureVel);
987 else
988 ✗ m_structureVelOldStoppingCriteria.copyFrom(m_structureVelExtrapol);
989 }
990 }
991
992 // compute the velocity of the structure
993 // dotX^n = ( X^n - X^(n-1) ) / tau
994 40 double tau = m_fstransient->timeStep;
995
2/8
✗ Branch 0 not taken.
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 40 times.
✗ Branch 8 not taken.
 40 if(m_timeSchemeType != 4 || (m_timeSchemeType == 4 && !m_pLinStruc->isStrucOperatorExplicitlyComputed())) {
996 40 m_structureVel.copyFrom(m_pLinStruc->solution());
997 40 m_structureVel.axpby(-1./tau, 1./tau, m_structureDispOld);
998 } else {
999 ✗ m_structureVelExtrapol.copyFrom(m_pLinStruc->solution());
1000 }
1001
1002
1003 // compute the extrapolation of the velocity
1004
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 40 times.
 40 if (isRNScheme) {
1005 ✗ switch (m_RNExtrapolOrder) {
1006 ✗ case 0:
1007 // ddot^{n-1} + dt \pd ddot^{\star} = ddot^{n-1}
1008 ✗ m_structureVelExtrapol.copyFrom(m_structureVel);
1009 ✗ break;
1010 ✗ case 1:
1011 // ddot^{n-1} + dt \pd ddot^{{\star} = 2*ddot^{n-1} - ddot^{n-2}
1012 ✗ m_structureVelExtrapol.copyFrom(m_structureVel);
1013 ✗ m_structureVelExtrapol.axpby(-1., 2. , m_structureVelOld);
1014 ✗ break;
1015 ✗ case 2:
1016 // ddot^{n-1} + dt \pd ddot^{{\star} = 3*ddot^{n-1} - 3*ddot^{n-2} + ddot^{n-3}
1017 ✗ m_structureVelExtrapol.copyFrom(m_structureVel);
1018 ✗ m_structureVelExtrapol.axpby(-3., 3., m_structureVelOld);
1019 ✗ m_structureVelExtrapol.axpy( 1. , m_structureVelOldOld);
1020 ✗ break;
1021 ✗ default:
1022 ✗ FEL_ERROR("unknown time scheme extrapolation order");
1023 ✗ break;
1024 }
1025
1/2
✓ Branch 0 taken 40 times.
✗ Branch 1 not taken.
 40 } else if(m_timeSchemeType == 3) {
1026 // No extrapolation is needed for scheme 3
1027 40 m_structureVelExtrapol.copyFrom(m_structureVel);
1028 }
1029
1030 // the serial std::vector is passed to the fluid problem
1031 40 m_pLinStruc->gatherVector(m_structureVelExtrapol, m_seqStructureVelExtrapol);
1032 40 }
1033
1034 template<class SolidPrb>
1035 ✗ void FSINitscheXFEMModel<SolidPrb>::computeTimeSchemeStructureExtrapolRHS()
1036 {
1037 ✗ switch(m_RNExtrapolOrder) {
1038 ✗ case 0:
1039 ✗ break;
1040
1041 ✗ case 1:
1042 ✗ m_structureExtrapolatedNitscheTermsRHS.copyFrom(m_pLinStruc->vector());
1043 ✗ m_structureExtrapolatedNitscheTermsRHS.axpby(-1., 1., m_structureTimeSchemeRHSTerms);
1044 ✗ break;
1045
1046 ✗ case 2:
1047 ✗ m_structureExtrapolatedNitscheTermsRHS.copyFrom(m_pLinStruc->vector());
1048 ✗ m_structureExtrapolatedNitscheTermsRHS.axpby(-1., 1., m_structureTimeSchemeRHSTerms);
1049 ✗ m_structureExtrapolatedNitscheTermsRHS.axpby(-1., 2. , m_structureNitscheRHSTermsOld);
1050
1051 ✗ m_structureNitscheRHSTermsOld.copyFrom(m_pLinStruc->vector());
1052 ✗ m_structureNitscheRHSTermsOld.axpby(-1., 1., m_structureTimeSchemeRHSTerms);
1053 ✗ break;
1054
1055 ✗ default:
1056 ✗ FEL_ERROR("unknown time scheme extrapolation order");
1057 ✗ break;
1058 }
1059 }
1060
1061 template<class SolidPrb>
1062 22 void FSINitscheXFEMModel<SolidPrb>::writeSolutionAndMeshes()
1063 {
1064 22 if( (m_fstransient->iteration % FelisceParam::instance().frequencyWriteSolution == 0)
1065 ✗ or m_hasFinished
1066
2/8
✗ Branch 0 not taken.
✓ Branch 1 taken 22 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 22 times.
✗ Branch 13 not taken.
 22 or ( m_fstransient->iteration >= FelisceParam::instance().intervalWriteSolution[0] && m_fstransient->iteration <= FelisceParam::instance().intervalWriteSolution[1] )) {
1067
1068 // write the mesh
1069
2/4
✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 22 times.
 22 if(FelisceParam::verbose()>1)
1070 ✗ PetscPrintf(MpiInfo::petscComm(),"Write meshes\n");
1071
1072 22 std::map<felInt, std::vector<felInt> > refToListOfIds;
1073 22 std::string fileName;
1074 22 std::string indexFile;
1075
1076
2/2
✓ Branch 1 taken 44 times.
✓ Branch 2 taken 22 times.
 66 for(std::size_t ipb=0; ipb<m_linearProblem.size(); ++ipb) {
1077
2/2
✓ Branch 4 taken 77 times.
✓ Branch 5 taken 44 times.
 121 for(std::size_t iUnknown=0; iUnknown<m_linearProblem[ipb]->listUnknown().size(); ++iUnknown) {
1078
1/2
✓ Branch 1 taken 77 times.
✗ Branch 2 not taken.
 77 std::stringstream oss;
1079
1/2
✓ Branch 1 taken 77 times.
✗ Branch 2 not taken.
 77 oss << m_indexTime;
1080
2/2
✓ Branch 0 taken 70 times.
✓ Branch 1 taken 7 times.
 77 if(m_indexTime < 10)
1081
2/4
✓ Branch 1 taken 70 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 70 times.
✗ Branch 5 not taken.
 70 indexFile = "0000" + oss.str();
1082
1/2
✓ Branch 0 taken 7 times.
✗ Branch 1 not taken.
 7 else if(m_indexTime < 100)
1083
2/4
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
 7 indexFile = "000" + oss.str();
1084 ✗ else if(m_indexTime < 1000)
1085 ✗ indexFile = "00" + oss.str();
1086 ✗ else if(m_indexTime < 10000)
1087 ✗ indexFile = "0" + oss.str();
1088 else
1089 ✗ indexFile = oss.str();
1090
1091 77 int idVar = m_linearProblem[ipb]->listUnknown().idVariable(iUnknown);
1092
6/12
✓ Branch 5 taken 77 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 77 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 77 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 77 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 77 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 77 times.
✗ Branch 21 not taken.
 77 fileName = FelisceParam::instance().resultDir + m_linearProblem[ipb]->listVariable().listVariable()[idVar].name()+".geo."+indexFile+".geo";
1093
1094
2/2
✓ Branch 0 taken 44 times.
✓ Branch 1 taken 33 times.
 77 if(ipb == 0) {
1095 // fluid, write mesh from the support dof to see the duplication
1096
2/2
✓ Branch 0 taken 22 times.
✓ Branch 1 taken 22 times.
 44 if(iUnknown == 0) {
1097 // write the mesh only once (with the velocity, the pressure mesh is the same)
1098
2/4
✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 22 times.
✗ Branch 4 not taken.
 22 if(FelisceParam::instance().outputFileFormat == 0)
1099
2/4
✓ Branch 2 taken 22 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 22 times.
✗ Branch 6 not taken.
 22 m_linearProblem[ipb]->writeGeoForUnknown(iUnknown, fileName);
1100 ✗ else if(FelisceParam::instance().outputFileFormat == 1) {
1101 ✗ m_linearProblem[ipb]->writeGeoForUnknownEnsightGold(iUnknown, refToListOfIds, fileName);
1102 ✗ m_ios[0]->ensight().setRefToListOfIds(refToListOfIds);
1103 }
1104 else
1105 ✗ FEL_ERROR("Unknown output file format");
1106 }
1107 } else {
1108 // structure, write mesh from the mesh points to see the displacement
1109
2/4
✓ Branch 3 taken 33 times.
✗ Branch 4 not taken.
✓ Branch 8 taken 33 times.
✗ Branch 9 not taken.
 33 m_ios[1]->writeMesh(*m_mesh[1], fileName);
1110 }
1111 }
1112 }
1113 22 m_meshIsWritten = true;
1114 22 ++m_indexTime;
1115
1116
1117 // write the solutions
1118
2/4
✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 22 times.
 22 if(FelisceParam::verbose() > 1)
1119 ✗ PetscPrintf(MpiInfo::petscComm(),"Write solutions\n");
1120
1121 //PetscVector aux;
1122 //aux.duplicateFrom(m_linearProblem[0]->solution());
1123 //aux.copyFrom(m_linearProblem[0]->solution());
1124 //if(m_fstransient->iteration > 0)
1125 // m_linearProblem[0]->solution().copyFrom(m_bdfFluid.sol_n());
1126
1127
2/2
✓ Branch 1 taken 44 times.
✓ Branch 2 taken 22 times.
 66 for (std::size_t ipb = 0; ipb < m_linearProblem.size(); ipb++) {
1128
2/4
✓ Branch 4 taken 44 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 44 times.
✗ Branch 8 not taken.
 44 m_linearProblem[ipb]->writeSolution(MpiInfo::rankProc(), m_ios, m_fstransient->time, m_fstransient->iteration);
1129 }
1130
1131 //m_linearProblem[0]->solution().copyFrom(aux);
1132 //m_linearProblem[0]->gatherSolution();
1133 //aux.destroy();
1134
1135
2/4
✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 22 times.
✗ Branch 4 not taken.
 22 if (MpiInfo::rankProc() == 0) {
1136
2/2
✓ Branch 1 taken 44 times.
✓ Branch 2 taken 22 times.
 66 for(std::size_t iio=0; iio<m_ios.size(); ++iio) {
1137
1/2
✓ Branch 2 taken 44 times.
✗ Branch 3 not taken.
 44 if(m_ios[iio]->typeOutput == 1 ) // 1 : ensight
1138
1/2
✓ Branch 4 taken 44 times.
✗ Branch 5 not taken.
 44 m_ios[iio]->postProcess(m_fstransient->time, 2);
1139 }
1140 }
1141 22 }
1142 22 }
1143
1144 template<class SolidPrb>
1145 20 void FSINitscheXFEMModel<SolidPrb>::updateInterfacePosition(PetscVector& disp)
1146 {
1147 20 const felInt dim = m_mesh[1]->numCoor();
1148 20 const felInt npt = m_mesh[1]->numPoints();
1149
1150
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 m_interfaceDisplacement.resize(dim * npt, 0.);
1151
1152 double* dispArray;
1153
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 disp.getArray(&dispArray);
1154
1155
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
 20 if( m_numStrucComp == 1 ) {
1156 10 auto& normals = m_mesh[1]->listNormals();
1157
1158 // displacement is scalar
1159
2/2
✓ Branch 0 taken 200 times.
✓ Branch 1 taken 10 times.
 210 for (felInt i = 0; i < npt; ++i)
1160
2/2
✓ Branch 0 taken 400 times.
✓ Branch 1 taken 200 times.
 600 for(felInt j = 0; j < dim; ++j)
1161 400 m_interfaceDisplacement[j + dim*i] = dispArray[i] * normals[i][j];
1162 } else {
1163 // displacement is a std::vector
1164
2/2
✓ Branch 0 taken 200 times.
✓ Branch 1 taken 10 times.
 210 for(felInt i = 0; i < npt; ++i)
1165
2/2
✓ Branch 0 taken 400 times.
✓ Branch 1 taken 200 times.
 600 for(felInt j = 0; j < dim; ++j)
1166 400 m_interfaceDisplacement[j + dim*i] = dispArray[j + dim*i];
1167 }
1168
1169
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 disp.restoreArray(&dispArray);
1170 20 }
1171 }
1172