GCC Code Coverage Report

Directory: ./
File: Solver/linearProblemFSINitscheXFEMFluid.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 1013 1358 74.6%
Branches: 926 2025 45.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/linearProblemFSINitscheXFEMFluid.hpp"
21 #include "Core/felisceTools.hpp"
22 #include "FiniteElement/elementVector.hpp"
23 #include "FiniteElement/elementMatrix.hpp"
24 #include "Geometry/neighFacesOfEdges.hpp"
25
26 namespace felisce {
27 2 LinearProblemFSINitscheXFEMFluid::LinearProblemFSINitscheXFEMFluid():
28 LinearProblem("Fictitious Domain Fluid Structure Interaction"),
29 2 m_viscosity(0.),
30 2 m_density(0.),
31
18/36
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 2 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 2 times.
✗ Branch 19 not taken.
✓ Branch 21 taken 2 times.
✗ Branch 22 not taken.
✓ Branch 24 taken 2 times.
✗ Branch 25 not taken.
✓ Branch 27 taken 2 times.
✗ Branch 28 not taken.
✓ Branch 30 taken 2 times.
✗ Branch 31 not taken.
✓ Branch 35 taken 2 times.
✗ Branch 36 not taken.
✓ Branch 38 taken 2 times.
✗ Branch 39 not taken.
✓ Branch 41 taken 2 times.
✗ Branch 42 not taken.
✓ Branch 44 taken 2 times.
✗ Branch 45 not taken.
✓ Branch 47 taken 2 times.
✗ Branch 48 not taken.
✓ Branch 50 taken 2 times.
✗ Branch 51 not taken.
✓ Branch 53 taken 2 times.
✗ Branch 54 not taken.
✓ Branch 60 taken 2 times.
✗ Branch 61 not taken.
✓ Branch 63 taken 2 times.
✗ Branch 64 not taken.
 2 m_bdf(nullptr) {
32 2 m_isSeqVelExtrapolAllocated = false;
33 2 m_isSeqBdfRHSAllocated = false;
34 2 m_isSeqOriginalSolAllocated = false;
35 2 m_isSeqRNFluidExtrapolAllocated = false;
36 2 m_feStruc = nullptr;
37 2 m_elemFieldTimeSchemeType = nullptr;
38 2 m_timeSchemeType = 0;
39 2 m_numOriginalDof = 0;
40 2 m_useCurrentFluidSolution = false;
41
42
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 m_seqVelExtrapol.resize(FelisceParam::instance().orderBdfNS);
43
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
 4 for(std::size_t i=0; i<m_seqVelExtrapol.size(); ++i)
44
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
 2 m_seqVelExtrapol[i] = PetscVector::null();
45
46
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 m_seqBdfRHS.resize(FelisceParam::instance().orderBdfNS);
47
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
 4 for(std::size_t i=0; i<m_seqBdfRHS.size(); ++i)
48
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
 2 m_seqBdfRHS[i] = PetscVector::null();
49 2 }
50
51 4 LinearProblemFSINitscheXFEMFluid::~LinearProblemFSINitscheXFEMFluid() {
52 4 m_bdf = nullptr;
53
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if (m_isSeqBdfRHSAllocated) {
54
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
 8 for(std::size_t i=0; i<m_seqBdfRHS.size(); ++i)
55 4 m_seqBdfRHS[i].destroy();
56 4 m_seqBdfRHS.clear();
57 }
58
59
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if (m_isSeqVelExtrapolAllocated) {
60
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
 8 for(std::size_t i=0; i<m_seqVelExtrapol.size(); ++i)
61 4 m_seqVelExtrapol[i].destroy();
62 4 m_seqVelExtrapol.clear();
63 }
64
65
2/2
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
 16 for(std::size_t iPt=0; iPt < m_mshPts.size(); ++iPt) {
66
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 12 delete m_mshPts[iPt];
67 }
68
69
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 2 times.
 12 for(std::size_t iPt=0; iPt < m_subItfPts.size(); ++iPt) {
70
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 8 delete m_subItfPts[iPt];
71 }
72
73
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if (m_isSeqOriginalSolAllocated)
74 4 m_tmpSeqOriginalSol.destroy();
75
76
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 2 times.
 12 for(std::size_t i=0; i<m_solutionOld.size(); ++i)
77 8 m_solutionOld[i].destroy();
78
79
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 2 times.
 12 for(std::size_t i=0; i<m_sequentialSolutionOld.size(); ++i)
80 8 m_sequentialSolutionOld[i].destroy();
81
82
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 2 times.
 12 for(std::size_t i=0; i<m_aoOld.size(); ++i)
83 8 AODestroy(&m_aoOld[i]);
84 }
85
86 2 void LinearProblemFSINitscheXFEMFluid::initialize(std::vector<GeometricMeshRegion::Pointer>& mesh, FelisceTransient::Pointer fstransient, MPI_Comm& comm, bool doUseSNES) {
87 // calling initialize in LinearProblem
88
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 LinearProblem::initialize(mesh,comm, doUseSNES);
89
90 // the problem
91 2 m_fstransient = fstransient;
92
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 setNumberOfMatrix(1);
93
94 // define variable of the linear system
95
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 std::vector<PhysicalVariable> listVariable(2);
96
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 std::vector<std::size_t> listNumComp(2);
97 2 listVariable[0] = velocity;
98 2 listNumComp[0] = 2; // because it only works for 2d in the fluid
99 2 listVariable[1] = pressure;
100 2 listNumComp[1] = 1;
101
102 //define unknown of the linear system.
103
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_listUnknown.push_back(velocity);
104
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_listUnknown.push_back(pressure);
105
106 // link between variable and unknowns
107
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 definePhysicalVariable(listVariable, listNumComp);
108
109 // data from the data file
110
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_viscosity = FelisceParam::instance().viscosity;
111
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_density = FelisceParam::instance().density;
112
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_useSymmetricStress = FelisceParam::instance().useSymmetricStress;
113
114 // get the time scheme type from the data file
115
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 m_elemFieldTimeSchemeType = FelisceParam::instance().elementFieldDynamicValue("TimeSchemeType");
116
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if(m_elemFieldTimeSchemeType != nullptr)
117
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_timeSchemeType = (felInt) m_elemFieldTimeSchemeType->constant(Comp1);
118
119
2/4
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
 2 if ((m_timeSchemeType == 0) || (m_timeSchemeType == 2)) {
120 // Robin-Neumann scheme
121 // extrapolation order
122 ✗ m_betaRN.resize(FelisceParam::instance().fsiRNscheme);
123 ✗ if (FelisceParam::instance().fsiRNscheme == 1){
124 ✗ m_betaRN[0] = 1.;
125 ✗ } else if (FelisceParam::instance().fsiRNscheme == 2){
126 ✗ m_betaRN[0] = 2.;
127 ✗ m_betaRN[1] = -1.;
128 }
129 }
130
131 // sub element
132 2 m_numVerPerFluidElt = m_mesh[m_currentMesh]->domainDim() + 1;
133 2 m_numVerPerSolidElt = m_mesh[m_currentMesh]->domainDim();
134
135 // create the std::vector of points for the fluid elements
136
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_mshPts.resize(m_numVerPerFluidElt);
137
2/2
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 2 times.
 8 for(std::size_t iver=0; iver < m_mshPts.size(); ++iver) {
138
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 12 m_mshPts[iver] = new Point();
139 }
140
141
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_itfPts.resize(m_numVerPerSolidElt);
142
143
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_subItfPts.resize(m_numVerPerSolidElt);
144
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 2 times.
 6 for(std::size_t iver=0; iver < m_subItfPts.size(); ++iver) {
145
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 8 m_subItfPts[iver] = new Point();
146 }
147 2 }
148
149
150 2 void LinearProblemFSINitscheXFEMFluid::addSolidVariable() {
151
2/2
✓ Branch 2 taken 3 times.
✓ Branch 3 taken 2 times.
 5 for(std::size_t iUnknown=0; iUnknown<m_pLinStruc->listUnknown().size(); ++iUnknown) {
152 3 m_listVariable.addVariable(m_pLinStruc->listUnknown()[iUnknown], m_pLinStruc->listVariable()[m_pLinStruc->listUnknown().idVariable(iUnknown)].numComponent());
153 }
154 2 }
155
156
157 160 void LinearProblemFSINitscheXFEMFluid::initPerElementType(ElementType eltType, FlagMatrixRHS flagMatrixRHS) {
158 IGNORE_UNUSED_ELT_TYPE;
159 IGNORE_UNUSED_FLAG_MATRIX_RHS;
160
161 // Get id of the variables
162 160 m_iVelocity = m_listVariable.getVariableIdList(velocity);
163 160 m_iPressure = m_listVariable.getVariableIdList(pressure);
164
165 // Get the finite elements
166 160 m_feVel = m_listCurrentFiniteElement[m_iVelocity];
167 160 m_fePres = m_listCurrentFiniteElement[m_iPressure];
168
169 // Get the variables of the problem
170 160 m_velocity = &m_listVariable[m_iVelocity];
171 160 m_pressure = &m_listVariable[m_iPressure];
172
173
174 // Initialize the element fields
175
2/4
✓ Branch 0 taken 160 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 160 times.
 160 if ((m_timeSchemeType == 0) || (m_timeSchemeType == 2)){
176 ✗ m_seqVelDofPts.resize(FelisceParam::instance().fsiRNscheme);
177 ✗ m_seqPreDofPts.resize(FelisceParam::instance().fsiRNscheme);
178 } else {
179 160 m_seqVelDofPts.resize(1);
180 160 m_seqPreDofPts.resize(1);
181 }
182 160 m_elemFieldAdv.initialize(DOF_FIELD, *m_feVel, dimension());
183 160 m_elemFieldAdvDup.initialize(DOF_FIELD, *m_feVel, dimension());
184 160 m_elemFieldAdvTmp.initialize(DOF_FIELD, *m_feVel, dimension());
185 160 m_elemFieldRHS.initialize(DOF_FIELD, *m_feVel, dimension());
186 160 m_elemFieldRHSbdf.initialize(DOF_FIELD, *m_feVel, dimension());
187 160 m_elemFieldTotalPressureFormulation.initialize(QUAD_POINT_FIELD, *m_feVel, dimension());
188 160 m_elemFieldNormal.initialize(CONSTANT_FIELD, m_feVel->numCoor());
189 160 m_elemFieldNormalInit.initialize(CONSTANT_FIELD, m_feVel->numCoor());
190 160 m_elemFieldAux.initialize(CONSTANT_FIELD, m_feVel->numCoor());
191
2/2
✓ Branch 1 taken 160 times.
✓ Branch 2 taken 160 times.
 320 for(std::size_t i=0; i<m_seqVelDofPts.size(); ++i){
192 160 m_seqVelDofPts[i].initialize(DOF_FIELD, *m_feVel, dimension());
193 160 m_seqPreDofPts[i].initialize(DOF_FIELD, *m_fePres, 1);
194 }
195
196 // solid
197 160 felInt idUnknownSolid = m_pLinStruc->listUnknown().getUnknownIdList(displacement);
198 160 felInt idVarSolid = m_pLinStruc->listUnknown().idVariable(idUnknownSolid);
199 160 m_strucVelQuadPts.initialize(QUAD_POINT_FIELD, *m_feVel, dimension());
200 160 m_strucVelDofPts.initialize(DOF_FIELD, *m_feStruc, m_pLinStruc->listVariable()[idVarSolid].numComponent());
201 160 m_ddotComp.initialize(QUAD_POINT_FIELD, *m_fePres, 1);
202 160 }
203
204
205 ✗ void LinearProblemFSINitscheXFEMFluid::initPerElementTypeBoundaryCondition(ElementType& eltType, FlagMatrixRHS flagMatrixRHS) {
206 IGNORE_UNUSED_ELT_TYPE;
207 IGNORE_UNUSED_FLAG_MATRIX_RHS;
208
209 // get the id of the variables
210 ✗ m_iVelocity = m_listVariable.getVariableIdList(velocity);
211 ✗ m_iPressure = m_listVariable.getVariableIdList(pressure);
212
213 // Get the curvilinear finite elements
214 ✗ m_curvFeVel = m_listCurvilinearFiniteElement[m_iVelocity];
215 ✗ m_curvFePress = m_listCurvilinearFiniteElement[m_iPressure];
216
217 // Get the variables of the problem
218 ✗ m_velocity = &m_listVariable[m_iVelocity];
219 ✗ m_pressure = &m_listVariable[m_iPressure];
220 }
221
222
223 // Virtual function in linearProblem. To allocate additional setValue. For Essential boundary conditions type
224 42 void LinearProblemFSINitscheXFEMFluid::finalizeEssBCTransientDerivedProblem() {
225
226 42 }
227
228
229
230 13 void LinearProblemFSINitscheXFEMFluid::userChangePattern(int numProc, int rankProc) {
231 // There may be integrals that involve the dof of the two duplicated support elements.
232 // The integrals for the ghost penalty modify also the pattern
233 // same for the face-oriented stabilization
234 (void) numProc;
235 (void) rankProc;
236
237 // definition of variables
238 felInt currentIelGeo, oppositeIelGeo;
239 felInt node1, node2;
240 felInt idVar1, idVar2;
241 felInt ielSupport, jelSupport;
242 13 std::vector<felInt> vecSupportCurrent, vecSupportOpposite;
243 bool isAnInnerFace;
244 bool computeBothSide;
245 bool addOppositeSupportElt;
246 bool isDomainElement;
247
248 GeometricMeshRegion::ElementType eltType, eltTypeOpp;
249 felInt numEltPerLabel;
250 felInt numElement[GeometricMeshRegion::m_numTypesOfElement];
251 felInt idByType;
252 felInt numFacesPerElement;
253
254 13 const std::vector<felInt>& iSupportDof = m_supportDofUnknown[0].iSupportDof();
255 13 const std::vector<felInt>& iEle = m_supportDofUnknown[0].iEle();
256 13 const std::vector<GeometricMeshRegion::ElementType>& bagDomain =m_mesh[m_currentMesh]->bagElementTypeDomain();
257
258 // the initial repartition of dof
259 13 felInt numDof = dof().numDof();
260
1/2
✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
 13 felInt numDofByProc = numDof/MpiInfo::numProc();
261 13 felInt numDofProc = 0;
262
1/2
✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
 13 felInt beginIdDofLocal = MpiInfo::rankProc()*numDofByProc;
263 felInt endIdDofLocal;
264
265
3/6
✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 13 times.
✗ Branch 7 not taken.
 13 if(MpiInfo::rankProc() == MpiInfo::numProc()-1) {
266
1/2
✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
 13 numDofProc = numDof - MpiInfo::rankProc()*numDofByProc;
267 13 endIdDofLocal = numDof;
268 } else {
269 ✗ numDofProc = numDofByProc;
270 ✗ endIdDofLocal = beginIdDofLocal + numDofProc;
271 }
272
273
274 // the complementary pattern
275 13 std::vector<felInt> iCSR, jCSR;
276
1/2
✓ Branch 2 taken 13 times.
✗ Branch 3 not taken.
 13 std::vector< std::set<felInt> > nodesNeighborhood(numDofProc);
277
278 // build the edges or faces depending on the dimension
279 // In this function, "faces" refers to either edges or faces.
280
3/10
✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 13 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 13 times.
✗ Branch 11 not taken.
 13 if(FelisceParam::instance().useGhostPenalty || FelisceParam::instance().NSStabType == 1) {
281
1/2
✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
 13 if(dimension() == 2)
282
1/2
✓ Branch 3 taken 13 times.
✗ Branch 4 not taken.
 13 m_mesh[m_currentMesh]->buildEdges();
283 else
284 ✗ m_mesh[m_currentMesh]->buildFaces();
285 }
286
287
288
289 // -------------------------------------------------------------------------- //
290 // Complementary pattern for : //
291 // - An interface with a tip inside the fluid domain
292 // - Robin-Neumann scheme (terms like u1 v2) //
293 // - ghost penalty without face-oriented stabilization (opposite element) //
294 // -------------------------------------------------------------------------- //
295 // loop over the unknowns of the problem
296
2/2
✓ Branch 1 taken 26 times.
✓ Branch 2 taken 13 times.
 39 for (std::size_t iUnknown1 = 0; iUnknown1 < m_listUnknown.size(); iUnknown1++) {
297 26 idVar1 = m_listUnknown.idVariable(iUnknown1);
298
299 // loop over all the intersected elements
300
1/2
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
 26 auto& idInterElt = m_duplicateSupportElements->getIntersectedMshElt();
301
2/2
✓ Branch 4 taken 970 times.
✓ Branch 5 taken 26 times.
 996 for (auto ielIt = idInterElt.begin(); ielIt != idInterElt.end(); ++ielIt) {
302 // get the geometric id of the element and the ids of its support elements
303 970 currentIelGeo = ielIt->first;
304
1/2
✓ Branch 3 taken 970 times.
✗ Branch 4 not taken.
 970 m_mesh[m_currentMesh]->getTypeElemFromIdElem(currentIelGeo, eltType, idByType);
305
306 970 isDomainElement = false;
307
2/2
✓ Branch 1 taken 970 times.
✓ Branch 2 taken 970 times.
 1940 for(std::size_t i=0; i<bagDomain.size(); ++i)
308
2/2
✓ Branch 1 taken 918 times.
✓ Branch 2 taken 52 times.
 970 if(bagDomain[i] == eltType)
309 918 isDomainElement = true;
310
311 970 numFacesPerElement = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second->numBdEle();
312
313
1/2
✓ Branch 2 taken 970 times.
✗ Branch 3 not taken.
 970 m_supportDofUnknown[iUnknown1].getIdElementSupport(currentIelGeo, vecSupportCurrent);
314
315 // loop over the support elements
316
2/2
✓ Branch 1 taken 1940 times.
✓ Branch 2 taken 970 times.
 2910 for(std::size_t ielSup=0; ielSup<vecSupportCurrent.size(); ++ielSup) {
317 1940 ielSupport = vecSupportCurrent[ielSup];
318
319 // for all support dof in this support element
320
2/2
✓ Branch 2 taken 5716 times.
✓ Branch 3 taken 1940 times.
 7656 for (felInt iSup = 0; iSup < m_supportDofUnknown[iUnknown1].getNumSupportDof(ielSupport); ++iSup) {
321 // for all component of the unknown
322
2/2
✓ Branch 2 taken 8574 times.
✓ Branch 3 taken 5716 times.
 14290 for (std::size_t iComp = 0; iComp < m_listVariable[idVar1].numComponent(); iComp++) {
323 // get the global id of the support dof
324
1/2
✓ Branch 2 taken 8574 times.
✗ Branch 3 not taken.
 8574 dof().loc2glob(ielSupport, iSup, idVar1, iComp, node1);
325
326 // If this support dof is owned by this process
327
2/4
✓ Branch 0 taken 8574 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 8574 times.
✗ Branch 3 not taken.
 8574 if(node1 >= beginIdDofLocal && node1 < endIdDofLocal) {
328 // for all unknown in the linear problem
329
2/2
✓ Branch 1 taken 17148 times.
✓ Branch 2 taken 8574 times.
 25722 for (std::size_t iUnknown2 = 0; iUnknown2 < m_listUnknown.size(); ++iUnknown2) {
330 17148 idVar2 = m_listUnknown.idVariable(iUnknown2);
331
332 // for all support element different than the ielSup
333 17148 std::size_t jelSup = 0;
334
2/2
✓ Branch 1 taken 34296 times.
✓ Branch 2 taken 17148 times.
 51444 while(jelSup < vecSupportCurrent.size()) {
335
2/2
✓ Branch 0 taken 17148 times.
✓ Branch 1 taken 17148 times.
 34296 if(jelSup != ielSup) {
336 17148 jelSupport = vecSupportCurrent[jelSup];
337
338 // for all support dof in this support element
339
2/2
✓ Branch 2 taken 50820 times.
✓ Branch 3 taken 17148 times.
 67968 for (felInt jSup = 0; jSup < supportDofUnknown(iUnknown2).getNumSupportDof(jelSupport); ++jSup) {
340 // for all component of the second unknown
341
2/2
✓ Branch 2 taken 76230 times.
✓ Branch 3 taken 50820 times.
 127050 for (std::size_t jComp = 0; jComp < m_listVariable[idVar2].numComponent(); jComp++) {
342 // If the two current components are connected
343
1/2
✓ Branch 2 taken 76230 times.
✗ Branch 3 not taken.
 76230 felInt iConnect = dof().getNumGlobComp(iUnknown1, iComp);
344
1/2
✓ Branch 2 taken 76230 times.
✗ Branch 3 not taken.
 76230 felInt jConnect = dof().getNumGlobComp(iUnknown2, jComp);
345
2/4
✓ Branch 2 taken 76230 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 76230 times.
✗ Branch 5 not taken.
 76230 if (m_listUnknown.mask()(iConnect, jConnect) > 0) {
346 // get the global id of the second support dof
347
1/2
✓ Branch 2 taken 76230 times.
✗ Branch 3 not taken.
 76230 dof().loc2glob(jelSupport, jSup, idVar2, jComp, node2);
348
2/2
✓ Branch 0 taken 76074 times.
✓ Branch 1 taken 156 times.
 76230 if(node1 != node2)
349
1/2
✓ Branch 2 taken 76074 times.
✗ Branch 3 not taken.
 76074 nodesNeighborhood[node1].insert(node2);
350 }
351 }
352 }
353 }
354 34296 ++jelSup;
355 }
356
357 // Ghost penalty but no face-oriented stabilization
358
5/12
✓ Branch 1 taken 17148 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 17148 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 17148 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 17148 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 17148 times.
 17148 if(FelisceParam::instance().useGhostPenalty && FelisceParam::instance().NSStabType != 1 && isDomainElement) {
359 // loop over the boundaries of the element
360 ✗ for(felInt iface=0; iface < numFacesPerElement; ++iface) {
361 // check if this face is an inner face or a boundary
362 ✗ oppositeIelGeo = currentIelGeo;
363 ✗ eltTypeOpp = eltType;
364 ✗ isAnInnerFace = m_mesh[m_currentMesh]->getAdjElement(eltTypeOpp, oppositeIelGeo, iface);
365
366 ✗ if(isAnInnerFace) {
367 // get the support element of the opposite element
368 ✗ supportDofUnknown(iUnknown2).getIdElementSupport(oppositeIelGeo, vecSupportOpposite);
369
370 // check if the opposite element is duplicated or not
371 ✗ if(vecSupportOpposite.size() > 1) {
372 // take the element on the same side
373 ✗ jelSupport = vecSupportOpposite[ielSup];
374 ✗ computeBothSide = false;
375 ✗ addOppositeSupportElt = true;
376 } else {
377 ✗ jelSupport = vecSupportOpposite[0];
378
379 // the opposite element is not duplicated, check if it is on the same side
380 ✗ std::size_t sizeCurrent = m_supportDofUnknown[0].getNumSupportDof(ielSupport);
381 ✗ std::size_t sizeOpposite = m_supportDofUnknown[0].getNumSupportDof(jelSupport);
382
383 ✗ felInt countSupportDofInCommon = 0;
384 ✗ for(std::size_t iSup1=0; iSup1<sizeCurrent; ++iSup1) {
385 ✗ for(std::size_t iSup2=0; iSup2<sizeOpposite; ++iSup2) {
386 ✗ if(iSupportDof[iEle[jelSupport] + iSup2] == iSupportDof[iEle[ielSupport] + iSup1]) {
387 ✗ ++countSupportDofInCommon;
388 }
389 }
390 }
391
392 ✗ if(countSupportDofInCommon > 1) {
393 // the elements are on the same side
394 ✗ computeBothSide = true;
395 ✗ addOppositeSupportElt = true;
396 } else {
397 ✗ computeBothSide = false;
398 ✗ addOppositeSupportElt = false;
399 }
400 }
401
402 ✗ if(addOppositeSupportElt) {
403 // for all support dof in this support element
404 ✗ for (felInt jSup = 0; jSup < supportDofUnknown(iUnknown2).getNumSupportDof(jelSupport); ++jSup) {
405 // for all component of the second unknown
406 ✗ for (std::size_t jComp = 0; jComp < m_listVariable[idVar2].numComponent(); jComp++) {
407 // If the two current components are connected
408 ✗ felInt iConnect = dof().getNumGlobComp(iUnknown1, iComp);
409 ✗ felInt jConnect = dof().getNumGlobComp(iUnknown2, jComp);
410 ✗ if (m_listUnknown.mask()(iConnect, jConnect) > 0) {
411 // get the global id of the second support dof
412 ✗ dof().loc2glob(jelSupport, jSup, idVar2, jComp, node2);
413 ✗ if(node1 != node2) {
414 ✗ nodesNeighborhood[node1].insert(node2);
415
416 ✗ if(computeBothSide)
417 ✗ nodesNeighborhood[node2].insert(node1);
418 }
419 }
420 }
421 }
422 }
423 }
424 }
425 }
426 }
427 }
428 }
429 }
430 }
431 }
432 }
433
434
435 // ----------------------------------------------------- //
436 // Complementary pattern for face-oriented stabilization //
437 // ----------------------------------------------------- //
438
2/4
✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 13 times.
✗ Branch 4 not taken.
 13 if(FelisceParam::instance().NSStabType == 1) {
439 // loop on the unknown of the linear problem
440
2/2
✓ Branch 1 taken 26 times.
✓ Branch 2 taken 13 times.
 39 for (std::size_t iUnknown1 = 0; iUnknown1 < m_listUnknown.size(); ++iUnknown1) {
441 26 idVar1 = m_listUnknown.idVariable(iUnknown1);
442 26 currentIelGeo = 0;
443
444 // first loop on element type
445
2/2
✓ Branch 4 taken 26 times.
✓ Branch 5 taken 26 times.
 52 for (std::size_t i=0; i<m_mesh[m_currentMesh]->bagElementTypeDomain().size(); ++i) {
446 26 eltType = m_mesh[m_currentMesh]->bagElementTypeDomain()[i];
447 26 const GeoElement* geoEle = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second;
448 26 numElement[eltType] = 0;
449 26 numFacesPerElement = geoEle->numBdEle();
450
451 // second loop on region of the mesh.
452
2/2
✓ Branch 8 taken 26 times.
✓ Branch 9 taken 26 times.
 52 for(auto itRef = m_mesh[m_currentMesh]->intRefToBegEndMaps[eltType].begin(); itRef !=m_mesh[m_currentMesh]->intRefToBegEndMaps[eltType].end(); itRef++) {
453 26 numEltPerLabel = itRef->second.second;
454
455 // Third loop on element
456
2/2
✓ Branch 0 taken 88972 times.
✓ Branch 1 taken 26 times.
 88998 for (felInt iel = 0; iel < numEltPerLabel; iel++) {
457
1/2
✓ Branch 2 taken 88972 times.
✗ Branch 3 not taken.
 88972 m_supportDofUnknown[iUnknown1].getIdElementSupport(eltType, numElement[eltType], vecSupportCurrent);
458
459
2/2
✓ Branch 1 taken 89890 times.
✓ Branch 2 taken 88972 times.
 178862 for(std::size_t ielSup=0; ielSup<vecSupportCurrent.size(); ++ielSup) {
460 89890 ielSupport = vecSupportCurrent[ielSup];
461
462 // for all support dof in this support element
463
2/2
✓ Branch 2 taken 269670 times.
✓ Branch 3 taken 89890 times.
 359560 for (felInt iSup = 0; iSup < m_supportDofUnknown[iUnknown1].getNumSupportDof(ielSupport); ++iSup) {
464 // for all component of the unknown
465
2/2
✓ Branch 2 taken 404505 times.
✓ Branch 3 taken 269670 times.
 674175 for (std::size_t iComp = 0; iComp < m_listVariable[idVar1].numComponent(); iComp++) {
466 // get the global id of the support dof
467
1/2
✓ Branch 2 taken 404505 times.
✗ Branch 3 not taken.
 404505 dof().loc2glob(ielSupport, iSup, idVar1, iComp, node1);
468
469 // if this node is on this process
470
2/4
✓ Branch 0 taken 404505 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 404505 times.
✗ Branch 3 not taken.
 404505 if(node1 >= beginIdDofLocal && node1 < endIdDofLocal) {
471 // loop over the boundaries of the element
472
2/2
✓ Branch 0 taken 1213515 times.
✓ Branch 1 taken 404505 times.
 1618020 for(felInt iface=0; iface < numFacesPerElement; ++iface) {
473 // check if this face is an inner face or a boundary
474 1213515 oppositeIelGeo = currentIelGeo;
475 1213515 eltTypeOpp = eltType;
476
1/2
✓ Branch 3 taken 1213515 times.
✗ Branch 4 not taken.
 1213515 isAnInnerFace = m_mesh[m_currentMesh]->getAdjElement(eltTypeOpp, oppositeIelGeo, iface);
477
478
2/2
✓ Branch 0 taken 1190349 times.
✓ Branch 1 taken 23166 times.
 1213515 if(isAnInnerFace) {
479 // for all unknown
480
2/2
✓ Branch 1 taken 2380698 times.
✓ Branch 2 taken 1190349 times.
 3571047 for (std::size_t iUnknown2 = 0; iUnknown2 < m_listUnknown.size(); ++iUnknown2) {
481 2380698 idVar2 = m_listUnknown.idVariable(iUnknown2);
482
483 // get the support element of the opposite element
484
1/2
✓ Branch 2 taken 2380698 times.
✗ Branch 3 not taken.
 2380698 supportDofUnknown(iUnknown2).getIdElementSupport(oppositeIelGeo, vecSupportOpposite);
485
486 // check if the opposite element and the current element are both duplicated
487
2/2
✓ Branch 2 taken 2354508 times.
✓ Branch 3 taken 26190 times.
 2380698 if(vecSupportCurrent.size() == vecSupportOpposite.size()) {
488 // take the element on the same side
489 2354508 jelSupport = vecSupportOpposite[ielSup];
490 2354508 addOppositeSupportElt = true;
491 } else {
492 26190 addOppositeSupportElt = false;
493 26190 jelSupport = -1;
494
6/6
✓ Branch 1 taken 34920 times.
✓ Branch 2 taken 21672 times.
✓ Branch 3 taken 30402 times.
✓ Branch 4 taken 4518 times.
✓ Branch 5 taken 30402 times.
✓ Branch 6 taken 26190 times.
 56592 for(std::size_t jelSup=0; jelSup < vecSupportOpposite.size() && !addOppositeSupportElt; ++jelSup) {
495 30402 jelSupport = vecSupportOpposite[jelSup];
496
497 // the opposite element is not duplicated, check if it is on the same side
498 30402 std::size_t sizeCurrent = m_supportDofUnknown[0].getNumSupportDof(ielSupport);
499 30402 std::size_t sizeOpposite = m_supportDofUnknown[0].getNumSupportDof(jelSupport);
500
501 30402 felInt countSupportDofInCommon = 0;
502
2/2
✓ Branch 0 taken 91206 times.
✓ Branch 1 taken 30402 times.
 121608 for(std::size_t iSup1=0; iSup1<sizeCurrent; ++iSup1) {
503
2/2
✓ Branch 0 taken 273618 times.
✓ Branch 1 taken 91206 times.
 364824 for(std::size_t iSup2=0; iSup2<sizeOpposite; ++iSup2) {
504
2/2
✓ Branch 4 taken 36324 times.
✓ Branch 5 taken 237294 times.
 273618 if(iSupportDof[iEle[jelSupport] + iSup2] == iSupportDof[iEle[ielSupport] + iSup1]) {
505 36324 ++countSupportDofInCommon;
506 }
507 }
508 }
509
510
2/2
✓ Branch 0 taken 17460 times.
✓ Branch 1 taken 12942 times.
 30402 if(countSupportDofInCommon > 1)
511 17460 addOppositeSupportElt = true;
512 }
513 }
514
515
2/2
✓ Branch 0 taken 2371968 times.
✓ Branch 1 taken 8730 times.
 2380698 if(addOppositeSupportElt) {
516 // for all support dof in this support element
517
2/2
✓ Branch 2 taken 7115904 times.
✓ Branch 3 taken 2371968 times.
 9487872 for (felInt jSup = 0; jSup < supportDofUnknown(iUnknown2).getNumSupportDof(jelSupport); ++jSup) {
518 // for all component of the second unknown
519
2/2
✓ Branch 2 taken 10673856 times.
✓ Branch 3 taken 7115904 times.
 17789760 for (std::size_t jComp = 0; jComp < m_listVariable[idVar2].numComponent(); jComp++) {
520 // If the two current components are connected
521
1/2
✓ Branch 2 taken 10673856 times.
✗ Branch 3 not taken.
 10673856 felInt iConnect = dof().getNumGlobComp(iUnknown1, iComp);
522
1/2
✓ Branch 2 taken 10673856 times.
✗ Branch 3 not taken.
 10673856 felInt jConnect = dof().getNumGlobComp(iUnknown2, jComp);
523
2/4
✓ Branch 2 taken 10673856 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 10673856 times.
✗ Branch 5 not taken.
 10673856 if (m_listUnknown.mask()(iConnect, jConnect) > 0) {
524 // get the global id of the second support dof
525
1/2
✓ Branch 2 taken 10673856 times.
✗ Branch 3 not taken.
 10673856 dof().loc2glob(jelSupport, jSup, idVar2, jComp, node2);
526
2/2
✓ Branch 0 taken 9883200 times.
✓ Branch 1 taken 790656 times.
 10673856 if(node1 != node2)
527
1/2
✓ Branch 2 taken 9883200 times.
✗ Branch 3 not taken.
 9883200 nodesNeighborhood[node1].insert(node2);
528 }
529 }
530 }
531 }
532 }
533 }
534 }
535 }
536 }
537 }
538 }
539 88972 ++currentIelGeo;
540 88972 ++numElement[eltType];
541 }
542 }
543 }
544 }
545 }
546
547
548 // ----------------------------------------------------------------- //
549 // build the complementary pattern and merge it with the current one //
550 // ----------------------------------------------------------------- //
551 13 felInt dofSize = 0;
552 felInt pos;
553
554
1/2
✓ Branch 4 taken 13 times.
✗ Branch 5 not taken.
 13 iCSR.resize(dof().pattern().numRows() + 1, 0);
555
2/2
✓ Branch 1 taken 72045 times.
✓ Branch 2 taken 13 times.
 72058 for(std::size_t iNode=0; iNode<nodesNeighborhood.size(); ++iNode)
556 72045 dofSize += nodesNeighborhood[iNode].size();
557
558
1/2
✓ Branch 0 taken 13 times.
✗ Branch 1 not taken.
 13 if(dofSize > 0) {
559
1/2
✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
 13 jCSR.resize(dofSize, 0);
560
2/2
✓ Branch 1 taken 72045 times.
✓ Branch 2 taken 13 times.
 72058 for(std::size_t iNode=0; iNode<nodesNeighborhood.size(); ++iNode) {
561 72045 iCSR[iNode + 1] = iCSR[iNode] + nodesNeighborhood[iNode].size();
562 72045 pos = 0;
563
2/2
✓ Branch 5 taken 2604318 times.
✓ Branch 6 taken 72045 times.
 2676363 for (auto it=nodesNeighborhood[iNode].begin(); it != nodesNeighborhood[iNode].end(); ++it) {
564 2604318 jCSR[iCSR[iNode] + pos] = *it;
565 2604318 ++pos;
566 }
567 }
568
569 // Now, call merge pattern
570
1/2
✓ Branch 2 taken 13 times.
✗ Branch 3 not taken.
 13 dof().mergeGlobalPattern(iCSR, jCSR);
571 }
572 13 }
573
574
575
576 141128 void LinearProblemFSINitscheXFEMFluid::computeElementArray(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint, felInt& iel1, felInt& iel2, ElementType& eltType, felInt& ielGeo, FlagMatrixRHS flagMatrixRHS) {
577 IGNORE_UNUSED_FLAG_MATRIX_RHS;
578 IGNORE_UNUSED_ELT_TYPE;
579 IGNORE_UNUSED_ELEM_ID_POINT;
580
581 141128 std::vector<felInt> strucElemIds;
582
583 // get the side of the support element
584 // here iel2 is refering to the solution because it will std::set the column in the global matrix
585 // (see the function setValueMatrixRHS in linearProblem.cpp)
586 sideOfInterface sideOfSupportElementForSolution;
587 double side;
588 felInt ielSupportDof;
589
1/2
✓ Branch 2 taken 141128 times.
✗ Branch 3 not taken.
 141128 m_supportDofUnknown[0].getIdElementSupport(ielGeo, ielSupportDof);
590
2/2
✓ Branch 0 taken 138296 times.
✓ Branch 1 taken 2832 times.
 141128 sideOfSupportElementForSolution = ielSupportDof == iel2 ? LEFT : RIGHT;
591
2/2
✓ Branch 0 taken 138296 times.
✓ Branch 1 taken 2832 times.
 141128 side = sideOfSupportElementForSolution == LEFT ? 1. : -1.;
592
593
2/2
✓ Branch 0 taken 138296 times.
✓ Branch 1 taken 2832 times.
 141128 if(iel1 == iel2) {
594 //////////////////////////////////
595 // VOLUME ELEMENT - SUB ELEMENT //
596 //////////////////////////////////
597
1/2
✓ Branch 1 taken 138296 times.
✗ Branch 2 not taken.
 138296 felInt numSubElement = m_duplicateSupportElements->getNumSubMshEltPerMshElt(ielGeo);
598
599
600 // -------------------------------------------
601 // Get the support element of the old solutions
602
1/2
✓ Branch 3 taken 138296 times.
✗ Branch 4 not taken.
 138296 std::vector<felInt> ielOld(m_solutionOld.size());
603
1/2
✓ Branch 3 taken 138296 times.
✗ Branch 4 not taken.
 138296 std::vector< std::vector<felInt> > oldSupportElement(m_solutionOld.size());
604
2/2
✓ Branch 1 taken 276592 times.
✓ Branch 2 taken 138296 times.
 414888 for(std::size_t i=0; i<oldSupportElement.size(); ++i) {
605
1/2
✓ Branch 4 taken 276592 times.
✗ Branch 5 not taken.
 276592 m_supportDofUnknownOld[i][0].getIdElementSupport(ielGeo, oldSupportElement[i]);
606
607 // std::set the id for the old solution
608
2/2
✓ Branch 2 taken 5552 times.
✓ Branch 3 taken 271040 times.
 276592 if(oldSupportElement[i].size() > 1) {
609 // this element is also duplicated for the old solution
610
2/2
✓ Branch 0 taken 2844 times.
✓ Branch 1 taken 2708 times.
 5552 ielOld[i] = sideOfSupportElementForSolution == LEFT ? oldSupportElement[i][0] : oldSupportElement[i][1];
611 } else {
612 // this element is NOT duplicated for the old solution
613 271040 ielOld[i] = oldSupportElement[i][0];
614 }
615 }
616
617
2/2
✓ Branch 0 taken 2832 times.
✓ Branch 1 taken 135464 times.
 138296 if(numSubElement > 0) {
618 // Volume sub elements
619 felInt subElemId;
620
1/2
✓ Branch 2 taken 2832 times.
✗ Branch 3 not taken.
 2832 std::vector<double> tmpPt(3);
621
622 // loop over the sub-elements
623
2/2
✓ Branch 0 taken 11216 times.
✓ Branch 1 taken 2832 times.
 14048 for(felInt iSubElt=0; iSubElt < numSubElement; ++iSubElt) {
624
1/2
✓ Branch 1 taken 11216 times.
✗ Branch 2 not taken.
 11216 subElemId = m_duplicateSupportElements->getSubMshEltIdxMsh(ielGeo, iSubElt);
625
626 // if the side of the sub-element is the same as the element
627
3/4
✓ Branch 1 taken 11216 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 5608 times.
✓ Branch 4 taken 5608 times.
 11216 if(m_duplicateSupportElements->getSubMshEltSide(ielGeo, iSubElt) == sideOfSupportElementForSolution) {
628 // --------------------------------------------------
629 // fill the std::vector with the points of the sub-element
630
2/2
✓ Branch 0 taken 16824 times.
✓ Branch 1 taken 5608 times.
 22432 for(std::size_t iPt=0; iPt < m_numVerPerFluidElt; ++iPt) {
631
1/2
✓ Branch 1 taken 16824 times.
✗ Branch 2 not taken.
 16824 m_duplicateSupportElements->getSubMshEltVerCrd(subElemId, iPt, tmpPt);
632 16824 *m_mshPts[iPt] = Point(tmpPt[0], tmpPt[1], tmpPt[2]);
633 }
634
635 // ---------------------
636 // Update finite element
637
1/2
✓ Branch 1 taken 5608 times.
✗ Branch 2 not taken.
 5608 m_feVel->updateSubElementFirstDeriv(0, elemPoint, m_mshPts);
638
1/2
✓ Branch 1 taken 5608 times.
✗ Branch 2 not taken.
 5608 m_fePres->updateSubElementFirstDeriv(0, elemPoint, m_mshPts);
639
640
641 // -----------------------
642 // build classic NS matrix
643
1/2
✓ Branch 1 taken 5608 times.
✗ Branch 2 not taken.
 5608 computeElementMatrixRHS();
644
3/4
✓ Branch 1 taken 11216 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 5608 times.
✓ Branch 4 taken 5608 times.
 11216 for(felInt i=0; i<FelisceParam::instance().orderBdfNS; ++i)
645
1/2
✓ Branch 2 taken 5608 times.
✗ Branch 3 not taken.
 5608 computeElementMatrixRHSPartDependingOnOldTimeStep(iel1, ielOld[i], i);
646 }
647 }
648 2832 } else {
649 // No volume sub elements
650
3/4
✓ Branch 1 taken 270928 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 135464 times.
✓ Branch 4 taken 135464 times.
 270928 for(felInt i=0; i<FelisceParam::instance().orderBdfNS; ++i) {
651
2/2
✓ Branch 2 taken 48 times.
✓ Branch 3 taken 135416 times.
 135464 if(oldSupportElement[i].size() > 1) {
652 // old element is duplicated, get the sub elements
653
1/2
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
 48 felInt numSubElementOld = m_duplicateSupportElements->getNumSubMshEltMshOld(ielGeo, i);
654
655 felInt subElemIdOld;
656
1/2
✓ Branch 2 taken 48 times.
✗ Branch 3 not taken.
 48 std::vector<double> tmpPtOld(3);
657 felInt ielOldTmp;
658
659
2/2
✓ Branch 0 taken 144 times.
✓ Branch 1 taken 48 times.
 192 for(felInt iSubElt=0; iSubElt < numSubElementOld; ++iSubElt) {
660
1/2
✓ Branch 1 taken 144 times.
✗ Branch 2 not taken.
 144 subElemIdOld = m_duplicateSupportElements->getSubMshEltIdxMshOld(ielGeo, iSubElt, i);
661
662 // --------------------------------------------------
663 // fill the std::vector with the points of the sub-element
664
2/2
✓ Branch 0 taken 432 times.
✓ Branch 1 taken 144 times.
 576 for(std::size_t iPt=0; iPt < m_numVerPerFluidElt; ++iPt) {
665
1/2
✓ Branch 1 taken 432 times.
✗ Branch 2 not taken.
 432 m_duplicateSupportElements->getSubMshEltVerCrdOld(subElemIdOld, iPt, tmpPtOld, i);
666 432 *m_mshPts[iPt] = Point(tmpPtOld[0], tmpPtOld[1], tmpPtOld[2]);
667 }
668
669 // ---------------------
670 // Update finite element
671
1/2
✓ Branch 1 taken 144 times.
✗ Branch 2 not taken.
 144 m_feVel->updateSubElementFirstDeriv(0, elemPoint, m_mshPts);
672
1/2
✓ Branch 1 taken 144 times.
✗ Branch 2 not taken.
 144 m_fePres->updateSubElementFirstDeriv(0, elemPoint, m_mshPts);
673
674
675 // --------------------------------------------------------
676 // Get the id of the support element owning the sub element
677
3/4
✓ Branch 1 taken 144 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 80 times.
✓ Branch 4 taken 64 times.
 144 ielOldTmp = m_duplicateSupportElements->getSubMshEltSideOld(subElemIdOld, i) == LEFT ? oldSupportElement[i][0] : oldSupportElement[i][1];
678
679
680 // -----------------------
681 // build classic NS matrix
682
1/2
✓ Branch 1 taken 144 times.
✗ Branch 2 not taken.
 144 computeElementMatrixRHSPartDependingOnOldTimeStep(iel1, ielOldTmp, i);
683 }
684 48 }
685 }
686
687 // Update finite element if needed
688
5/6
✓ Branch 1 taken 134104 times.
✓ Branch 2 taken 1360 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 134104 times.
✓ Branch 6 taken 1360 times.
✓ Branch 7 taken 134104 times.
 135464 if(!m_feVel->hasOriginalQuadPoint() || !m_fePres->hasOriginalQuadPoint()) {
689
1/2
✓ Branch 1 taken 1360 times.
✗ Branch 2 not taken.
 1360 m_feVel->updateBasisAndFirstDeriv(0, elemPoint);
690
1/2
✓ Branch 1 taken 1360 times.
✗ Branch 2 not taken.
 1360 m_fePres->updateBasisAndFirstDeriv(0, elemPoint);
691 } else {
692
1/2
✓ Branch 1 taken 134104 times.
✗ Branch 2 not taken.
 134104 m_feVel->updateFirstDeriv(0, elemPoint);
693
1/2
✓ Branch 1 taken 134104 times.
✗ Branch 2 not taken.
 134104 m_fePres->updateFirstDeriv(0, elemPoint);
694 }
695
696 // build classic NS matrix
697
1/2
✓ Branch 1 taken 135464 times.
✗ Branch 2 not taken.
 135464 computeElementMatrixRHS();
698
699 // build matrix and rhs depending on old time step if this element is not intersected
700 // by the old solutions
701
3/4
✓ Branch 1 taken 270928 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 135464 times.
✓ Branch 4 taken 135464 times.
 270928 for(felInt i=0; i<FelisceParam::instance().orderBdfNS; ++i) {
702
2/2
✓ Branch 2 taken 135416 times.
✓ Branch 3 taken 48 times.
 135464 if(oldSupportElement[i].size() == 1)
703
1/2
✓ Branch 2 taken 135416 times.
✗ Branch 3 not taken.
 135416 computeElementMatrixRHSPartDependingOnOldTimeStep(iel1, ielOld[i], i);
704 }
705 }
706
707
708
709 ////////////////////////////
710 // SUB STRUCTURE ELEMENTS //
711 ////////////////////////////
712
1/2
✓ Branch 1 taken 138296 times.
✗ Branch 2 not taken.
 138296 numSubElement = m_duplicateSupportElements->getNumSubItfEltPerMshElt(ielGeo);
713
2/2
✓ Branch 0 taken 2832 times.
✓ Branch 1 taken 135464 times.
 138296 if(numSubElement > 0) {
714 felInt subElemId;
715 felInt elemId;
716
1/2
✓ Branch 2 taken 2832 times.
✗ Branch 3 not taken.
 2832 std::vector<double> tmpPt(3);
717
1/2
✓ Branch 3 taken 2832 times.
✗ Branch 4 not taken.
 2832 std::vector<double> itfNormal(dimension(),0);
718
1/2
✓ Branch 3 taken 2832 times.
✗ Branch 4 not taken.
 2832 std::vector<double> itfNormalInit(dimension(),0);
719
720
2/2
✓ Branch 0 taken 4112 times.
✓ Branch 1 taken 2832 times.
 6944 for(felInt iSubElt=0; iSubElt < numSubElement; ++iSubElt) {
721 // get the id of the sub element
722
1/2
✓ Branch 1 taken 4112 times.
✗ Branch 2 not taken.
 4112 subElemId = m_duplicateSupportElements->getSubItfEltIdxMsh(ielGeo, iSubElt);
723
724 // get the id of the element owing this sub-element
725
1/2
✓ Branch 1 taken 4112 times.
✗ Branch 2 not taken.
 4112 elemId = m_duplicateSupportElements->getItfEltIdxOfSubItfElt(subElemId);
726
727 // get the normal to this sub element
728 4112 GeometricMeshRegion::ElementType eltTypeStruc = m_pLinStruc->mesh()->bagElementTypeDomain()[0];
729
1/2
✓ Branch 4 taken 4112 times.
✗ Branch 5 not taken.
 4112 m_pLinStruc->mesh()->listElementNormal(eltTypeStruc, elemId).getCoor(itfNormal);
730
1/2
✓ Branch 4 taken 4112 times.
✗ Branch 5 not taken.
 4112 m_pLinStruc->mesh()->listElementNormalInit(eltTypeStruc, elemId).getCoor(itfNormalInit);
731
732 // fill the std::vector with the points of the sub-element
733
2/2
✓ Branch 0 taken 8224 times.
✓ Branch 1 taken 4112 times.
 12336 for(std::size_t iPt=0; iPt < m_numVerPerSolidElt; ++iPt) {
734
1/2
✓ Branch 1 taken 8224 times.
✗ Branch 2 not taken.
 8224 m_duplicateSupportElements->getSubItfEltVerCrd(subElemId, iPt, tmpPt);
735 8224 *m_subItfPts[iPt] = Point(tmpPt[0], tmpPt[1], tmpPt[2]);
736 }
737
738 // Get the coordinates of the element owning this sub-element
739 4112 int numPointsPerElt = GeometricMeshRegion::m_numPointsPerElt[eltTypeStruc];
740
1/2
✓ Branch 2 taken 4112 times.
✗ Branch 3 not taken.
 4112 std::vector<felInt> strucElemPtsId(numPointsPerElt);
741
742
1/2
✓ Branch 3 taken 4112 times.
✗ Branch 4 not taken.
 4112 m_pLinStruc->mesh()->getOneElement(eltTypeStruc, elemId, strucElemPtsId, 0);
743
744
2/2
✓ Branch 0 taken 8224 times.
✓ Branch 1 taken 4112 times.
 12336 for (int iPoint = 0; iPoint < numPointsPerElt; ++iPoint)
745 8224 m_itfPts[iPoint] = &m_pLinStruc->mesh()->listPoint(strucElemPtsId[iPoint]);
746
747
748 // update the finite element
749
1/2
✓ Branch 1 taken 4112 times.
✗ Branch 2 not taken.
 4112 m_feVel->updateSubElementFirstDeriv(0, elemPoint, m_subItfPts, m_feStruc);
750
1/2
✓ Branch 1 taken 4112 times.
✗ Branch 2 not taken.
 4112 m_fePres->updateSubElementFirstDeriv(0, elemPoint, m_subItfPts, m_feStruc);
751
752 // Set element field
753
6/8
✓ Branch 1 taken 4112 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4112 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2056 times.
✓ Branch 6 taken 2056 times.
✓ Branch 7 taken 2056 times.
✓ Branch 8 taken 2056 times.
 4112 if(FelisceParam::instance().useProjectionForNitscheXFEM || m_useCurrentFluidSolution) {
754
2/2
✓ Branch 1 taken 2056 times.
✓ Branch 2 taken 2056 times.
 4112 for(std::size_t i=0; i<m_seqVelDofPts.size(); ++i){
755
2/4
✓ Branch 0 taken 2056 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2056 times.
 2056 if ((m_timeSchemeType == 0) || (m_timeSchemeType == 2)){
756 ✗ m_seqVelDofPts[i].setValue(m_seqRNFluidExtrapol[i], *m_feVel, iel2, m_iVelocity, m_ao, dof());
757 ✗ m_seqPreDofPts[i].setValue(m_seqRNFluidExtrapol[i], *m_fePres, iel2, m_iPressure, m_ao, dof());
758 }else{
759
1/2
✓ Branch 4 taken 2056 times.
✗ Branch 5 not taken.
 2056 m_seqVelDofPts[i].setValue(sequentialSolution(), *m_feVel, iel2, m_iVelocity, m_ao, dof());
760
1/2
✓ Branch 4 taken 2056 times.
✗ Branch 5 not taken.
 2056 m_seqPreDofPts[i].setValue(sequentialSolution(), *m_fePres, iel2, m_iPressure, m_ao, dof());
761 }
762 }
763 } else {
764
2/2
✓ Branch 1 taken 2056 times.
✓ Branch 2 taken 2056 times.
 4112 for(std::size_t i=0; i<m_seqVelDofPts.size(); ++i){
765
1/2
✓ Branch 6 taken 2056 times.
✗ Branch 7 not taken.
 2056 m_seqVelDofPts[i].setValue(m_sequentialSolutionOld[i], *m_feVel, ielOld[i], m_iVelocity, m_aoOld[i], m_dofOld[i]);
766
1/2
✓ Branch 6 taken 2056 times.
✗ Branch 7 not taken.
 2056 m_seqPreDofPts[i].setValue(m_sequentialSolutionOld[i], *m_fePres, ielOld[i], m_iPressure, m_aoOld[i], m_dofOld[i]);
767 }
768 }
769
770
1/2
✓ Branch 1 taken 4112 times.
✗ Branch 2 not taken.
 4112 m_elemFieldNormal.setValue(itfNormal);
771
1/2
✓ Branch 1 taken 4112 times.
✗ Branch 2 not taken.
 4112 m_elemFieldNormalInit.setValue(itfNormalInit);
772
773
774 // Inflow stabilization for NS
775
5/10
✓ Branch 1 taken 4112 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4112 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 4112 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 4112 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 4112 times.
✗ Branch 11 not taken.
 4112 if(FelisceParam::instance().NSequationFlag == 1 && FelisceParam::instance().useInterfaceFlowStab) {
776 // compute (\beta \cdot n)_- u_i v_i
777
2/4
✓ Branch 4 taken 4112 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4112 times.
✗ Branch 8 not taken.
 4112 m_elementMat[0]->f_dot_n_phi_i_phi_j(0.5 * side * FelisceParam::instance().density, *m_feVel, m_elemFieldAdv, m_elemFieldNormal, 0, 0, dimension(), 0);
778 }
779
780 // compute the elementary matrix and rhs
781
1/6
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 4112 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 4112 switch(m_timeSchemeType) {
782 ✗ case 0: {
783 // Robin-Neumann schemes (direct implementation) (0th, 1st and 2nd order extrapolated variants)
784
785 // MATRIX //
786 // \rho_s \epsilon / dt \int_S u_i . v_i
787 ✗ double penaltyParam = FelisceParam::instance().densitySolid * FelisceParam::instance().thickness / m_fstransient->timeStep;
788 ✗ m_elementMat[0]->phi_i_phi_j(penaltyParam, *m_feVel, 0, 0, dimension());
789
790
791 // RHS //
792 // update the value of the structure on that sub element (std::set m_strucVelQuadPts)
793 ✗ updateStructureVel(elemId, externalVec(0));
794 ✗ updateSubStructureVel(m_itfPts, m_subItfPts);
795
796 // \rho_s \epsilon / dt \int_S ( ddot^{n-1} + dt \pd ddot^{\star} ) . v_i
797 ✗ m_elementVector[0]->source(penaltyParam, *m_feVel, m_strucVelQuadPts, 0, dimension());
798
799 ✗ for(std::size_t i=0; i<m_seqVelDofPts.size(); ++i){
800 // velocity
801 ✗ if(m_useSymmetricStress)
802 ✗ m_elementVector[0]->eps_vec2_dot_vec1_dot_phi_i(-2. * side * m_viscosity * m_betaRN[i], m_elemFieldNormal, m_seqVelDofPts[i], *m_feVel, 0, dimension());
803 else
804 ✗ m_elementVector[0]->grad_vec2_dot_vec1_dot_phi_i(-side * m_viscosity * m_betaRN[i], m_elemFieldNormal, m_seqVelDofPts[i], *m_feVel, 0, dimension());
805
806 // pressure
807 ✗ for(felInt idim=0; idim<dimension(); ++idim)
808 ✗ m_elementVector[0]->source(side * m_elemFieldNormal.val(idim,0) * m_betaRN[i], *m_feVel, m_seqPreDofPts[i], idim, 1);
809 }
810
811 ✗ break;
812 }
813
814 ✗ case 1: {
815 // Robin-Robin
816 // MATRIX //
817 // Nitsche's penalty term, \gamma \nu / h u_i . v_i
818 ✗ double penaltyParam = m_viscosity*FelisceParam::instance().NitschePenaltyParam/m_feVel->diameter();
819 ✗ m_elementMat[0]->phi_i_phi_j(penaltyParam, *m_feVel, 0, 0, dimension());
820
821 // RHS //
822 // Nitsche's consistency term in the RHS //
823 // velocity
824 ✗ if(m_useSymmetricStress)
825 ✗ m_elementVector[0]->eps_vec2_dot_vec1_dot_phi_i(-2. * side * m_viscosity, m_elemFieldNormal, m_seqVelDofPts[0], *m_feVel, 0, dimension());
826 else
827 ✗ m_elementVector[0]->grad_vec2_dot_vec1_dot_phi_i(-side * m_viscosity, m_elemFieldNormal, m_seqVelDofPts[0], *m_feVel, 0, dimension());
828
829 // pressure
830 ✗ for(felInt idim=0; idim<dimension(); ++idim)
831 ✗ m_elementVector[0]->source(side * m_elemFieldNormal.val(idim,0), *m_feVel, m_seqPreDofPts[0], idim, 1);
832
833 // update the value of the solid velocity (std::set m_strucVelDofPts and m_strucVelQuadPts)
834 // m_strucVelDofPts.val.clear();
835 ✗ updateStructureVel(elemId, externalVec(0));
836
837 // m_strucVelQuadPts.val.clear();
838 ✗ updateSubStructureVel(m_itfPts, m_subItfPts);
839
840 // penalty term (RHS part)
841 ✗ m_elementVector[0]->source(penaltyParam, *m_feVel, m_strucVelQuadPts, 0, dimension());
842
843 ✗ break;
844 }
845
846 ✗ case 2: {
847 // Robin-Neumann schemes a la Stenberg (0th, 1st and 2nd order extrapolated variants)
848
849 // MATRIX //
850
851 ✗ double alpha = 0.5;
852 ✗ double RNParam = FelisceParam::instance().densitySolid * FelisceParam::instance().thickness / m_fstransient->timeStep;
853 ✗ double penaltyParam1 = (m_feVel->diameter()) / (m_viscosity*FelisceParam::instance().NitschePenaltyParam + RNParam*m_feVel->diameter());
854 ✗ double penaltyParam2 = (FelisceParam::instance().NitschePenaltyParam*m_viscosity ) / (m_viscosity*FelisceParam::instance().NitschePenaltyParam + RNParam*m_feVel->diameter());
855
856 // + u_i . v_i
857 ✗ m_elementMat[0]->phi_i_phi_j(alpha * penaltyParam2 * RNParam, *m_feVel, 0, 0, dimension());
858
859 // - \sigma(u_i, p_i)n_i . v_i
860 ✗ double coef = alpha * RNParam * penaltyParam1 + (1 - alpha);
861 ✗ computeNitscheSigma_u_v(coef * side * m_viscosity, coef * -side);
862
863 // - \sigma(0, -q_i)n_i . u_i (flip of sign because of the implementation of NS terms crossed with q)
864 ✗ computeNitscheSigma_v_u(0., alpha * penaltyParam1 * RNParam * -side);
865
866 // - \sigma(u_i, p_i)n_i . \sigma(0, -q_i)n_i (flip of sign because of the implementation of NS terms crossed with q)
867 ✗ computeStenbergSigma_u_Sigma_v( m_viscosity * alpha * penaltyParam1 , - alpha * penaltyParam1 );
868
869 // RHS //
870
871 // update the value of the structure on that sub element (std::set m_strucVelQuadPts)
872 ✗ updateStructureVel(elemId, externalVec(0));
873 ✗ updateSubStructureVel(m_itfPts, m_subItfPts);
874
875 // + (ddot^{n-1} + dt \pd ddot^{\star}) . v_i
876 ✗ m_elementVector[0]->source(alpha * penaltyParam2 * RNParam , *m_feVel, m_strucVelQuadPts, 0, dimension());
877
878 // - (ddot^{n-1} + dt \pd ddot^{\star}) . \sigma(0, -q_i) n_i (flip of sign because of the implementation of NS terms crossed with q)
879 ✗ computeNitscheDpoint_Sigma_v(0., -side * alpha * penaltyParam1 * RNParam);
880
881 // + \sigma( u^{\star}, p^{\star})n_i . v_i
882 ✗ for(std::size_t i=0; i<m_seqVelDofPts.size(); ++i){
883 // velocity
884 ✗ if(m_useSymmetricStress)
885 ✗ m_elementVector[0]->eps_vec2_dot_vec1_dot_phi_i(-2. * side * m_viscosity * m_betaRN[i] * alpha * penaltyParam2 , m_elemFieldNormal, m_seqVelDofPts[i], *m_feVel, 0, dimension());
886 else
887 ✗ m_elementVector[0]->grad_vec2_dot_vec1_dot_phi_i(-side * m_viscosity * m_betaRN[i] * alpha * penaltyParam2, m_elemFieldNormal, m_seqVelDofPts[i], *m_feVel, 0, dimension());
888
889 // pressure
890 ✗ for(felInt idim=0; idim<dimension(); ++idim)
891 ✗ m_elementVector[0]->source(side * m_elemFieldNormal.val(idim,0) * m_betaRN[i] * alpha * penaltyParam2, *m_feVel, m_seqPreDofPts[i], idim, 1);
892 }
893
894 // - \sigma( u^{\star}, p^{\star})n_i . \sigma(0, -q_i)n_i
895 ✗ for(std::size_t i=0; i<m_seqVelDofPts.size(); ++i){
896 ✗ std::vector<double> vec(dimension());
897 ✗ for(felInt idim=0; idim<dimension(); ++idim){
898 ✗ for(felInt j=0; j<dimension(); ++j){
899 ✗ vec[j] = 2 * m_elemFieldNormal.val(j,0) * m_elemFieldNormal.val(idim,0);
900 }
901 ✗ m_elemFieldAux.setValue(vec);
902 ✗ if(m_useSymmetricStress)
903 ✗ m_elementVector[0]->eps_vec2_dot_vec1_dot_phi_i(2. * m_viscosity * m_betaRN[i] * alpha * penaltyParam1 , m_elemFieldAux, m_seqVelDofPts[i], *m_feVel, dimension(), 1);
904 else
905 ✗ m_elementVector[0]->grad_vec2_dot_vec1_dot_phi_i( m_viscosity * m_betaRN[i] * alpha * penaltyParam1 , m_elemFieldNormal, m_seqVelDofPts[i], *m_feVel, dimension(), 1);
906 }
907 ✗ m_elementVector[0]->source(- m_betaRN[i] * alpha * penaltyParam1, *m_feVel, m_seqPreDofPts[i], dimension(), 1);
908 }
909
910
911 ✗ break;
912 }
913
914 4112 case 3: {
915 ///////////////////////////////
916 // stabilize explicit scheme //
917 ///////////////////////////////
918
919 // update the value of the solid velocity (std::set m_strucVelDofPts and m_strucVelQuadPts)
920
1/2
✓ Branch 2 taken 4112 times.
✗ Branch 3 not taken.
 4112 updateStructureVel(elemId, externalVec(0));
921
1/2
✓ Branch 1 taken 4112 times.
✗ Branch 2 not taken.
 4112 updateSubStructureVel(m_itfPts, m_subItfPts);
922
923
2/4
✓ Branch 1 taken 4112 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 4112 times.
 4112 if(FelisceParam::instance().usePenaltyFreeNitscheMethod) {
924 ✗ double pressureStabilization = 1./FelisceParam::instance().NitschePenaltyParam;
925
926 /*** MATRIX ***/
927 // Nitsche's consistency terms, \sigma(v_i, q_i)n_i . u_i - \sigma(u_i, p_i)n_i . v_i
928 ✗ computeNitscheSigma_u_v(side * m_viscosity, -side);
929 ✗ computeNitscheSigma_v_u(-side * m_viscosity, -side);
930
931 // stabilization of the coupling with the structure
932 // m_elementMat[0]->phi_i_phi_j(-pressureStabilization * m_feVel->diameter() / m_viscosity, *m_fePres, dimension(), dimension(), 1);
933 ✗ computeNitscheSigma_u_Sigma_v(pressureStabilization * m_feVel->diameter() / m_viscosity);
934
935 /*** RHS ***/
936 // -ddot \sigma(v_i, q_i) n_i //
937 ✗ computeNitscheDpoint_Sigma_v(-side * m_viscosity, -side);
938
939 // stabilization of the coupling with the structure
940 // m_elementVector[0]->source(-pressureStabilization * m_feVel->diameter() / m_viscosity, *m_fePres, m_seqPreDofPts[0], dimension(), 1);
941 ✗ computeNitscheSigma_uold_Sigma_v(pressureStabilization * m_feVel->diameter() / m_viscosity);
942 } else {
943 4112 double pressureStabilization = 1.;
944
2/4
✓ Branch 1 taken 4112 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4112 times.
✗ Branch 5 not taken.
 4112 double penaltyParam = m_viscosity*FelisceParam::instance().NitschePenaltyParam/m_feVel->diameter();
945
946 /*** MATRIX ***/
947 // Nitsche's penalty term, \gamma \mu / h u_i . v_i
948
1/2
✓ Branch 4 taken 4112 times.
✗ Branch 5 not taken.
 4112 m_elementMat[0]->phi_i_phi_j(penaltyParam, *m_feVel, 0, 0, dimension());
949
950 // Nitsche's consistency terms, - \sigma(v_i, -q_i)n_i . u_i
951
1/2
✓ Branch 1 taken 4112 times.
✗ Branch 2 not taken.
 4112 computeNitscheSigma_v_u(0., -side);
952
953 // stabilization of the coupling with the structure
954
1/2
✓ Branch 5 taken 4112 times.
✗ Branch 6 not taken.
 4112 m_elementMat[0]->phi_i_phi_j(-pressureStabilization/penaltyParam, *m_fePres, dimension(), dimension(), 1);
955
956 /*** RHS ***/
957 // penalty term, \gamma \mu / h \dot{d}^n v_i
958
1/2
✓ Branch 4 taken 4112 times.
✗ Branch 5 not taken.
 4112 m_elementVector[0]->source(penaltyParam, *m_feVel, m_strucVelQuadPts, 0, dimension());
959
960 // ddot \sigma(v_i, -q_i) n_i //
961
1/2
✓ Branch 1 taken 4112 times.
✗ Branch 2 not taken.
 4112 computeNitscheDpoint_Sigma_v(0., -side);
962
963 // stabilization of the coupling with the structure
964
1/2
✓ Branch 5 taken 4112 times.
✗ Branch 6 not taken.
 4112 m_elementVector[0]->source(-pressureStabilization/penaltyParam, *m_fePres, m_seqPreDofPts[0], dimension(), 1);
965
966
967 // RHS //
968 // Nitsche's consistency term in the RHS \sigma(u_i^{n-1}, p_i^{n-1}n_i . v_i
969 // velocity
970
1/2
✓ Branch 0 taken 4112 times.
✗ Branch 1 not taken.
 4112 if(m_useSymmetricStress)
971
1/2
✓ Branch 5 taken 4112 times.
✗ Branch 6 not taken.
 4112 m_elementVector[0]->eps_vec2_dot_vec1_dot_phi_i(-2. * side * m_viscosity, m_elemFieldNormal, m_seqVelDofPts[0], *m_feVel, 0, dimension());
972 else
973 ✗ m_elementVector[0]->grad_vec2_dot_vec1_dot_phi_i(-side * m_viscosity, m_elemFieldNormal, m_seqVelDofPts[0], *m_feVel, 0, dimension());
974
975 // pressure
976
2/2
✓ Branch 1 taken 8224 times.
✓ Branch 2 taken 4112 times.
 12336 for(felInt idim=0; idim<dimension(); ++idim)
977
2/4
✓ Branch 4 taken 8224 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 8224 times.
✗ Branch 8 not taken.
 8224 m_elementVector[0]->source(side * m_elemFieldNormal.val(idim,0), *m_feVel, m_seqPreDofPts[0], idim, 1);
978 }
979
980 4112 break;
981 }
982
983 ✗ case 4: {
984 // Robin-Neumann semi-implict schemes (0th, 1st and 2nd order extrapolated variants)
985 // MATRIX
986 // Nitsche's penalty term, \gamma \nu / h u_i . v_i
987 ✗ double pressureStabilization = 1.;
988 ✗ double penaltyParam = m_viscosity*FelisceParam::instance().NitschePenaltyParam/m_feVel->diameter();
989 ✗ m_elementMat[0]->phi_i_phi_j(penaltyParam, *m_feVel, 0, 0, dimension());
990
991 // Nitsche's consistency terms, -\sigma(u_i, p_i)n_i . v_i - \sigma(v_i, q_i)n_i . u_i
992 // u_i & v_i
993 ✗ computeNitscheSigma_v_u(0., -side);
994
995 // stabilization of the coupling with the structure
996 ✗ m_elementMat[0]->phi_i_phi_j(-pressureStabilization/penaltyParam, *m_fePres, dimension(), dimension(), 1);
997
998 // RHS //
999 // Nitsche's consistency term in the RHS //
1000 // velocity
1001 ✗ if(m_useSymmetricStress)
1002 ✗ m_elementVector[0]->eps_vec2_dot_vec1_dot_phi_i(-2. * side * m_viscosity, m_elemFieldNormal, m_seqVelDofPts[0], *m_feVel, 0, dimension());
1003 else
1004 ✗ m_elementVector[0]->grad_vec2_dot_vec1_dot_phi_i(-side * m_viscosity, m_elemFieldNormal, m_seqVelDofPts[0], *m_feVel, 0, dimension());
1005
1006 // pressure
1007 ✗ for(felInt idim=0; idim<dimension(); ++idim)
1008 ✗ m_elementVector[0]->source(side * m_elemFieldNormal.val(idim,0), *m_feVel, m_seqPreDofPts[0], idim, 1);
1009
1010
1011 // RHS //
1012 // Nitsche's consistency term in the RHS //
1013
1014 // update the value of the solid velocity (std::set m_strucVelDofPts and m_strucVelQuadPts)
1015 // m_strucVelDofPts.val.clear();
1016 ✗ updateStructureVel(elemId, m_pLinStruc->sequentialSolution());
1017
1018 // m_strucVelQuadPts.val.clear();
1019 ✗ updateSubStructureVel(m_itfPts, m_subItfPts);
1020
1021 // penalty term (RHS part)
1022 ✗ m_elementVector[0]->source(penaltyParam, *m_feVel, m_strucVelQuadPts, 0, dimension());
1023
1024 // ddot \sigma(v_i, q_i) n_i //
1025 // computeNitscheDpoint_Sigma_v(side * m_viscosity, -side);
1026 ✗ computeNitscheDpoint_Sigma_v(0., -side);
1027
1028 // stabilization of the coupling with the structure
1029 ✗ m_elementVector[0]->source(-pressureStabilization/penaltyParam, *m_fePres, m_seqPreDofPts[0], dimension(), 1);
1030 ✗ break;
1031 }
1032
1033 ✗ default: {
1034 ✗ FEL_ERROR("Bad Time splitting scheme");
1035 ✗ break;
1036 }
1037 }
1038 4112 }
1039 2832 }
1040 138296 } else {
1041
1042 // -------------------------------------------
1043 // Get the support element of the old solutions
1044
1/2
✓ Branch 3 taken 2832 times.
✗ Branch 4 not taken.
 2832 std::vector<felInt> ielOld(m_solutionOld.size());
1045
1/2
✓ Branch 3 taken 2832 times.
✗ Branch 4 not taken.
 2832 std::vector<felInt> ielOldDup(m_solutionOld.size());
1046
1/2
✓ Branch 3 taken 2832 times.
✗ Branch 4 not taken.
 2832 std::vector< std::vector<felInt> > oldSupportElement(m_solutionOld.size());
1047
2/2
✓ Branch 1 taken 5664 times.
✓ Branch 2 taken 2832 times.
 8496 for(std::size_t i=0; i<oldSupportElement.size(); ++i) {
1048
1/2
✓ Branch 4 taken 5664 times.
✗ Branch 5 not taken.
 5664 m_supportDofUnknownOld[i][0].getIdElementSupport(ielGeo, oldSupportElement[i]);
1049
1050 // std::set the id for the old solution
1051
2/2
✓ Branch 2 taken 5416 times.
✓ Branch 3 taken 248 times.
 5664 if(oldSupportElement[i].size() > 1) {
1052 // this element is also duplicated for the old solution
1053
2/2
✓ Branch 0 taken 2708 times.
✓ Branch 1 taken 2708 times.
 5416 if(sideOfSupportElementForSolution == LEFT) {
1054 2708 ielOld[i] = oldSupportElement[i][0];
1055 2708 ielOldDup[i] = oldSupportElement[i][1];
1056 } else {
1057 2708 ielOld[i] = oldSupportElement[i][1];
1058 2708 ielOldDup[i] = oldSupportElement[i][0];
1059 }
1060 } else {
1061 // this element is NOT duplicated for the old solution
1062 248 ielOld[i] = oldSupportElement[i][0];
1063 248 ielOldDup[i] = oldSupportElement[i][0];
1064 }
1065 }
1066
1067 ////////////////////////////
1068 // SUB STRUCTURE ELEMENTS //
1069 ////////////////////////////
1070
1/2
✓ Branch 1 taken 2832 times.
✗ Branch 2 not taken.
 2832 felInt numSubElement = m_duplicateSupportElements->getNumSubItfEltPerMshElt(ielGeo);
1071
1/2
✓ Branch 0 taken 2832 times.
✗ Branch 1 not taken.
 2832 if(numSubElement > 0) {
1072 felInt subElemId;
1073 felInt elemId;
1074
1/2
✓ Branch 2 taken 2832 times.
✗ Branch 3 not taken.
 2832 std::vector<double> tmpPt(3);
1075
1/2
✓ Branch 3 taken 2832 times.
✗ Branch 4 not taken.
 2832 std::vector<double> itfNormal(dimension(),0);
1076
1/2
✓ Branch 3 taken 2832 times.
✗ Branch 4 not taken.
 2832 std::vector<double> itfNormalInit(dimension(),0);
1077
1078 // std::set the extrapolation of velocity
1079
1/2
✓ Branch 1 taken 2832 times.
✗ Branch 2 not taken.
 2832 m_elemFieldAdv.val.clear();
1080
1/2
✓ Branch 1 taken 2832 times.
✗ Branch 2 not taken.
 2832 m_elemFieldAdvDup.val.clear();
1081
2/4
✓ Branch 1 taken 2832 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2832 times.
✗ Branch 4 not taken.
 2832 if(FelisceParam::instance().NSequationFlag == 1) {
1082
2/4
✓ Branch 1 taken 2832 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 2832 times.
 2832 if(FelisceParam::instance().useProjectionForNitscheXFEM) {
1083 ✗ for(felInt i=0; i<FelisceParam::instance().orderBdfNS; ++i) {
1084 ✗ m_elemFieldAdvTmp.setValue(m_seqVelExtrapol[i], *m_feVel, iel1, m_iVelocity, m_ao, dof());
1085 ✗ m_elemFieldAdv.val += m_elemFieldAdvTmp.val;
1086 ✗ m_elemFieldAdvTmp.setValue(m_seqVelExtrapol[i], *m_feVel, iel2, m_iVelocity, m_ao, dof());
1087 ✗ m_elemFieldAdvDup.val += m_elemFieldAdvTmp.val;
1088 }
1089 } else {
1090
3/4
✓ Branch 1 taken 5664 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2832 times.
✓ Branch 4 taken 2832 times.
 5664 for(felInt i=0; i<FelisceParam::instance().orderBdfNS; ++i) {
1091
1/2
✓ Branch 5 taken 2832 times.
✗ Branch 6 not taken.
 2832 m_elemFieldAdvTmp.setValue(m_seqVelExtrapol[i], *m_feVel, ielOld[i], m_iVelocity, m_aoOld[i], m_dofOld[i]);
1092
1/2
✓ Branch 1 taken 2832 times.
✗ Branch 2 not taken.
 2832 m_elemFieldAdv.val += m_elemFieldAdvTmp.val;
1093
1/2
✓ Branch 5 taken 2832 times.
✗ Branch 6 not taken.
 2832 m_elemFieldAdvTmp.setValue(m_seqVelExtrapol[i], *m_feVel, ielOldDup[i], m_iVelocity, m_aoOld[i], m_dofOld[i]);
1094
1/2
✓ Branch 1 taken 2832 times.
✗ Branch 2 not taken.
 2832 m_elemFieldAdvDup.val += m_elemFieldAdvTmp.val;
1095 }
1096 }
1097 }
1098
1099
1100
2/2
✓ Branch 0 taken 4112 times.
✓ Branch 1 taken 2832 times.
 6944 for(felInt iSubElt=0; iSubElt < numSubElement; ++iSubElt) {
1101 // get the id of the sub element
1102
1/2
✓ Branch 1 taken 4112 times.
✗ Branch 2 not taken.
 4112 subElemId = m_duplicateSupportElements->getSubItfEltIdxMsh(ielGeo, iSubElt);
1103
1104 // get the id of the element owing this sub element
1105
1/2
✓ Branch 1 taken 4112 times.
✗ Branch 2 not taken.
 4112 elemId = m_duplicateSupportElements->getItfEltIdxOfSubItfElt(subElemId);
1106
1107 // get the normal to this sub element
1108 4112 GeometricMeshRegion::ElementType eltTypeStruc = m_pLinStruc->mesh()->bagElementTypeDomain()[0];
1109
1/2
✓ Branch 4 taken 4112 times.
✗ Branch 5 not taken.
 4112 m_pLinStruc->mesh()->listElementNormal(eltTypeStruc, elemId).getCoor(itfNormal);
1110
1/2
✓ Branch 4 taken 4112 times.
✗ Branch 5 not taken.
 4112 m_pLinStruc->mesh()->listElementNormalInit(eltTypeStruc, elemId).getCoor(itfNormalInit);
1111
1112 // fill the std::vector with the points of the sub-element
1113
2/2
✓ Branch 0 taken 8224 times.
✓ Branch 1 taken 4112 times.
 12336 for(std::size_t iPt=0; iPt < m_numVerPerSolidElt; ++iPt) {
1114
1/2
✓ Branch 1 taken 8224 times.
✗ Branch 2 not taken.
 8224 m_duplicateSupportElements->getSubItfEltVerCrd(subElemId, iPt, tmpPt);
1115 8224 *m_subItfPts[iPt] = Point(tmpPt[0], tmpPt[1], tmpPt[2]);
1116 }
1117
1118 // Get the coordinates of the element owning this sub-element
1119 4112 int numPointsPerElt = GeometricMeshRegion::m_numPointsPerElt[eltTypeStruc];
1120
1/2
✓ Branch 2 taken 4112 times.
✗ Branch 3 not taken.
 4112 std::vector<felInt> strucElemPtsId(numPointsPerElt);
1121
1122
1/2
✓ Branch 3 taken 4112 times.
✗ Branch 4 not taken.
 4112 m_pLinStruc->mesh()->getOneElement(eltTypeStruc, elemId, strucElemPtsId, 0);
1123
1124
2/2
✓ Branch 0 taken 8224 times.
✓ Branch 1 taken 4112 times.
 12336 for (int iPoint = 0; iPoint < numPointsPerElt; ++iPoint)
1125 8224 m_itfPts[iPoint] = &m_pLinStruc->mesh()->listPoint(strucElemPtsId[iPoint]);
1126
1127
1128 // update the finite element
1129
1/2
✓ Branch 1 taken 4112 times.
✗ Branch 2 not taken.
 4112 m_feVel->updateSubElementFirstDeriv(0, elemPoint, m_subItfPts, m_feStruc);
1130
1/2
✓ Branch 1 taken 4112 times.
✗ Branch 2 not taken.
 4112 m_fePres->updateSubElementFirstDeriv(0, elemPoint, m_subItfPts, m_feStruc);
1131
1132 // Set element field
1133
6/8
✓ Branch 1 taken 4112 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4112 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2056 times.
✓ Branch 6 taken 2056 times.
✓ Branch 7 taken 2056 times.
✓ Branch 8 taken 2056 times.
 4112 if(FelisceParam::instance().useProjectionForNitscheXFEM || m_useCurrentFluidSolution > 0) {
1134
2/2
✓ Branch 1 taken 2056 times.
✓ Branch 2 taken 2056 times.
 4112 for(std::size_t i=0; i<m_seqVelDofPts.size(); ++i){
1135
2/4
✓ Branch 0 taken 2056 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2056 times.
 2056 if ((m_timeSchemeType == 0) || (m_timeSchemeType == 2)){
1136 ✗ m_seqVelDofPts[i].setValue(m_seqRNFluidExtrapol[i], *m_feVel, iel2, m_iVelocity, m_ao, dof());
1137 ✗ m_seqPreDofPts[i].setValue(m_seqRNFluidExtrapol[i], *m_fePres, iel2, m_iPressure, m_ao, dof());
1138 } else {
1139
1/2
✓ Branch 4 taken 2056 times.
✗ Branch 5 not taken.
 2056 m_seqVelDofPts[i].setValue(sequentialSolution(), *m_feVel, iel2, m_iVelocity, m_ao, dof());
1140
1/2
✓ Branch 4 taken 2056 times.
✗ Branch 5 not taken.
 2056 m_seqPreDofPts[i].setValue(sequentialSolution(), *m_fePres, iel2, m_iPressure, m_ao, dof());
1141 }
1142 }
1143 } else {
1144
2/2
✓ Branch 1 taken 2056 times.
✓ Branch 2 taken 2056 times.
 4112 for(std::size_t i=0; i<m_seqVelDofPts.size(); ++i){
1145
1/2
✓ Branch 6 taken 2056 times.
✗ Branch 7 not taken.
 2056 m_seqVelDofPts[i].setValue(m_sequentialSolutionOld[i], *m_feVel, ielOld[i], m_iVelocity, m_aoOld[i], m_dofOld[i]);
1146
1/2
✓ Branch 6 taken 2056 times.
✗ Branch 7 not taken.
 2056 m_seqPreDofPts[i].setValue(m_sequentialSolutionOld[i], *m_fePres, ielOld[i], m_iPressure, m_aoOld[i], m_dofOld[i]);
1147 }
1148 }
1149
1150
1/2
✓ Branch 1 taken 4112 times.
✗ Branch 2 not taken.
 4112 m_elemFieldNormal.setValue(itfNormal);
1151
1/2
✓ Branch 1 taken 4112 times.
✗ Branch 2 not taken.
 4112 m_elemFieldNormalInit.setValue(itfNormalInit);
1152
1153 // compute Inflow stabilization for NS equation
1154
5/10
✓ Branch 1 taken 4112 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4112 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 4112 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 4112 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 4112 times.
✗ Branch 11 not taken.
 4112 if(FelisceParam::instance().NSequationFlag == 1 && FelisceParam::instance().useInterfaceFlowStab) {
1155
2/4
✓ Branch 4 taken 4112 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4112 times.
✗ Branch 8 not taken.
 4112 m_elementMat[0]->f_dot_n_phi_i_phi_j(0.25 * side * FelisceParam::instance().density, *m_feVel, m_elemFieldAdv, m_elemFieldNormal, 0, 0, dimension(), 0);
1156
2/4
✓ Branch 4 taken 4112 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4112 times.
✗ Branch 8 not taken.
 4112 m_elementMat[0]->f_dot_n_phi_i_phi_j(0.25 * side * FelisceParam::instance().density, *m_feVel, m_elemFieldAdvDup, m_elemFieldNormal, 0, 0, dimension(), 0);
1157 }
1158
1159 // compute the elementary matrix and rhs
1160
1/6
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 4112 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 4112 switch(m_timeSchemeType) {
1161
1162 ✗ case 0: {
1163 // Robin-Neumann schemes (0th, 1st and 2nd order extrapolated variants)
1164 // MATRIX //
1165 ✗ computeNitscheSigma_u_v(-1 * side * m_viscosity, side);
1166
1167 // RHS //
1168 ✗ for(std::size_t i=0; i<m_seqVelDofPts.size(); ++i){
1169 // velocity
1170 ✗ if(m_useSymmetricStress)
1171 ✗ m_elementVector[0]->eps_vec2_dot_vec1_dot_phi_i(-2. * side * m_viscosity * m_betaRN[i], m_elemFieldNormal, m_seqVelDofPts[i], *m_feVel, 0, dimension());
1172 else
1173 ✗ m_elementVector[0]->grad_vec2_dot_vec1_dot_phi_i(-side * m_viscosity * m_betaRN[i], m_elemFieldNormal, m_seqVelDofPts[i], *m_feVel, 0, dimension());
1174
1175 // pressure
1176 ✗ for(felInt idim=0; idim<dimension(); ++idim)
1177 ✗ m_elementVector[0]->source(side * m_elemFieldNormal.val(idim,0) * m_betaRN[i], *m_feVel, m_seqPreDofPts[i], idim, 1);
1178 }
1179 ✗ break;
1180 }
1181
1182 ✗ case 1: {
1183 ✗ break;
1184 }
1185
1186 ✗ case 2: {
1187 // Robin-Neumann schemes a la Stenberg (0th, 1st and 2nd order extrapolated variants)
1188
1189 ✗ double alpha = 0.5;
1190 ✗ double RNParam = FelisceParam::instance().densitySolid * FelisceParam::instance().thickness / m_fstransient->timeStep;
1191 ✗ double penaltyParam1 = (m_feVel->diameter()) / (m_viscosity*FelisceParam::instance().NitschePenaltyParam + RNParam*m_feVel->diameter());
1192 ✗ double penaltyParam2 = (FelisceParam::instance().NitschePenaltyParam*m_viscosity ) / (m_viscosity*FelisceParam::instance().NitschePenaltyParam + RNParam*m_feVel->diameter());
1193
1194 // MATRIX //
1195
1196 // + \sigma(u_i, p_i)n_i . v_j
1197 ✗ computeNitscheSigma_u_v(-1 * side * m_viscosity * alpha * penaltyParam2, side * alpha * penaltyParam2);
1198
1199 // - \sigma(u_i, p_i)n_i . \sigma(0, -q_i)n_i
1200 ✗ computeStenbergSigma_u_Sigma_v( -m_viscosity * alpha * penaltyParam1 , alpha * penaltyParam1 );
1201
1202 // RHS //
1203
1204 // + \sigma( u^{\star}, p^{\star})n_i . v_i
1205 ✗ for(std::size_t i=0; i<m_seqVelDofPts.size(); ++i){
1206 // velocity
1207 ✗ if(m_useSymmetricStress)
1208 ✗ m_elementVector[0]->eps_vec2_dot_vec1_dot_phi_i(-2. * side * m_viscosity * m_betaRN[i] * alpha * penaltyParam2 , m_elemFieldNormal, m_seqVelDofPts[i], *m_feVel, 0, dimension());
1209 else
1210 ✗ m_elementVector[0]->grad_vec2_dot_vec1_dot_phi_i(-side * m_viscosity * m_betaRN[i] * alpha * penaltyParam2, m_elemFieldNormal, m_seqVelDofPts[i], *m_feVel, 0, dimension());
1211
1212 // pressure
1213 ✗ for(felInt idim=0; idim<dimension(); ++idim)
1214 ✗ m_elementVector[0]->source(side * m_elemFieldNormal.val(idim,0) * m_betaRN[i] * alpha * penaltyParam2, *m_feVel, m_seqPreDofPts[i], idim, 1);
1215 }
1216
1217 // - \sigma( u^{\star}, p^{\star})n_i . \sigma(0, -q_i)n_i
1218 ✗ for(std::size_t i=0; i<m_seqVelDofPts.size(); ++i){
1219 ✗ std::vector<double> vec(dimension());
1220 ✗ for(felInt idim=0; idim<dimension(); ++idim){
1221 ✗ for(felInt j=0; j<dimension(); ++j){
1222 ✗ vec[j] = 2 * m_elemFieldNormal.val(j,0) * m_elemFieldNormal.val(idim,0);
1223 }
1224 ✗ m_elemFieldAux.setValue(vec);
1225 ✗ if(m_useSymmetricStress)
1226 ✗ m_elementVector[0]->eps_vec2_dot_vec1_dot_phi_i(-2. * m_viscosity * m_betaRN[i] * alpha * penaltyParam1 , m_elemFieldAux, m_seqVelDofPts[i], *m_feVel, dimension(), 1);
1227 else
1228 ✗ m_elementVector[0]->grad_vec2_dot_vec1_dot_phi_i( -m_viscosity * m_betaRN[i] * alpha * penaltyParam1 , m_elemFieldNormal, m_seqVelDofPts[i], *m_feVel, dimension(), 1);
1229 }
1230 ✗ m_elementVector[0]->source(m_betaRN[i] * alpha * penaltyParam1, *m_feVel, m_seqPreDofPts[i], dimension(), 1);
1231 }
1232
1233 ✗ break;
1234 }
1235
1236 4112 case 3:{
1237 4112 break;
1238 }
1239
1240 ✗ case 4:{
1241 ✗ break;
1242 }
1243
1244 ✗ default: {
1245 ✗ FEL_ERROR("Bad Time splitting scheme");
1246 ✗ break;
1247 }
1248 }
1249 4112 }
1250 2832 }
1251 2832 }
1252 141128 }
1253
1254
1255
1256 40 void LinearProblemFSINitscheXFEMFluid::assembleMatrixFrontPoints() {
1257 // Assemble the matrix to take the crack point into account
1258 // Enforce continuity accros the fictitious interface
1259 40 felInt nbrVerFro = m_duplicateSupportElements->getNumFrontPoints();
1260
1261
1/2
✓ Branch 0 taken 40 times.
✗ Branch 1 not taken.
 40 if(nbrVerFro > 0) {
1262 40 std::vector<felInt> verFro;
1263 40 std::vector<felInt> verFro_temp; // generic front vertex std::vector for checking multi fronts within element
1264 felInt idElem;
1265 40 GeometricMeshRegion::ElementType eltType = m_mesh[m_currentMesh]->bagElementTypeDomain()[0];
1266
1267 40 felInt numPointPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
1268
1/2
✓ Branch 2 taken 40 times.
✗ Branch 3 not taken.
 40 std::vector<Point*> elemPoint(numPointPerElt);
1269
1/2
✓ Branch 2 taken 40 times.
✗ Branch 3 not taken.
 40 std::vector<felInt> elemIdPoint(numPointPerElt, 0);
1270
1271 40 std::vector<felInt> vectorIdSupport;
1272
1/2
✓ Branch 3 taken 40 times.
✗ Branch 4 not taken.
 40 std::vector<double> itfNormal(dimension(), 0.);
1273
1274 // define finite element
1275
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 defineFiniteElement(eltType);
1276
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 initElementArray();
1277
1278 // allocate arrays for assembling the matrix
1279 40 FlagMatrixRHS flag = FlagMatrixRHS::only_matrix;
1280
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 allocateArrayForAssembleMatrixRHS(flag);
1281
1282 // init variables
1283
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 initPerElementType(eltType, flag);
1284
1285
2/2
✓ Branch 0 taken 80 times.
✓ Branch 1 taken 40 times.
 120 for(felInt iVer=0; iVer<nbrVerFro; ++iVer) {
1286 // get the info on the front point
1287
1/2
✓ Branch 1 taken 80 times.
✗ Branch 2 not taken.
 80 m_duplicateSupportElements->getFrontPointInfo(iVer, verFro);
1288
1289 // get the fluid element in which the fictitious interface is
1290 80 idElem = verFro[1];
1291
1292 // get the support elements
1293
1/2
✓ Branch 1 taken 80 times.
✗ Branch 2 not taken.
 80 setElemPoint(eltType, idElem, elemPoint, elemIdPoint, vectorIdSupport);
1294
1295 // fill the std::vector with the points of the fictitious interface
1296 80 *m_subItfPts[0] = Point(m_pLinStruc->mesh()->listPoint(verFro[0]));
1297 80 *m_subItfPts[1] = Point(m_mesh[m_currentMesh]->listPoint(verFro[2]));
1298
1299 // now check if two fronts are in same element--then opposite vertex is other front)
1300
2/2
✓ Branch 0 taken 160 times.
✓ Branch 1 taken 80 times.
 240 for(felInt iVer2=0; iVer2<nbrVerFro; ++iVer2) {
1301
2/2
✓ Branch 0 taken 80 times.
✓ Branch 1 taken 80 times.
 160 if(iVer2!=iVer){
1302
1/2
✓ Branch 1 taken 80 times.
✗ Branch 2 not taken.
 80 m_duplicateSupportElements->getFrontPointInfo(iVer2, verFro_temp);
1303
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 80 times.
 80 if(verFro[1]==verFro_temp[1]){
1304 ✗ *m_subItfPts[1] = Point(m_pLinStruc->mesh()->listPoint(verFro[2]));
1305 }
1306 }
1307 }
1308 // special initial case for pre-determined crack point when nbrVerFro is 1 for one crack initially with the crack point being the same
1309 // NOTE THIS CURRENTLY ONLY WORKS FOR ONE CRACK IN A SOLID
1310
3/8
✓ Branch 1 taken 80 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 80 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 80 times.
 80 if(FelisceParam::instance().hasCrack && nbrVerFro==1) {
1311 ✗ *m_subItfPts[1] = Point(m_pLinStruc->mesh()->listPoint(verFro[2]));
1312 }
1313
1314
1315
1316 // update the finite element
1317
1/2
✓ Branch 1 taken 80 times.
✗ Branch 2 not taken.
 80 m_feVel->updateSubElementFirstDeriv(0, elemPoint, m_subItfPts, m_feStruc);
1318
1/2
✓ Branch 1 taken 80 times.
✗ Branch 2 not taken.
 80 m_fePres->updateSubElementFirstDeriv(0, elemPoint, m_subItfPts, m_feStruc);
1319
1320 // compute the normal
1321 80 itfNormal[0] = -(m_subItfPts[1]->y() - m_subItfPts[0]->y());
1322 80 itfNormal[1] = m_subItfPts[1]->x() - m_subItfPts[0]->x();
1323
1324 80 double norm = std::sqrt(itfNormal[0]*itfNormal[0] + itfNormal[1]*itfNormal[1]);
1325 80 itfNormal[0] /= norm;
1326 80 itfNormal[1] /= norm;
1327
1/2
✓ Branch 1 taken 80 times.
✗ Branch 2 not taken.
 80 m_elemFieldNormal.setValue(itfNormal);
1328
1329 // parameters
1330
2/4
✓ Branch 1 taken 80 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 80 times.
✗ Branch 5 not taken.
 80 double penaltyParam = m_viscosity*FelisceParam::instance().NitschePenaltyParam/m_feVel->diameter();
1331
1332 // loop over the support elements
1333
2/2
✓ Branch 1 taken 160 times.
✓ Branch 2 taken 80 times.
 240 for(std::size_t iSup1=0; iSup1 < vectorIdSupport.size(); ++iSup1) {
1334 160 felInt iel1 = vectorIdSupport[iSup1];
1335
2/2
✓ Branch 1 taken 320 times.
✓ Branch 2 taken 160 times.
 480 for(std::size_t iSup2=0; iSup2 < vectorIdSupport.size(); ++iSup2) {
1336 320 felInt iel2 = vectorIdSupport[iSup2];
1337
1338 // clear elementary matrix
1339
2/2
✓ Branch 1 taken 320 times.
✓ Branch 2 taken 320 times.
 640 for (std::size_t j = 0, size = numberOfMatrices(); j < size ; j++)
1340
1/2
✓ Branch 3 taken 320 times.
✗ Branch 4 not taken.
 320 m_elementMat[j]->zero();
1341
1342
1343 // get the side
1344 sideOfInterface sideOfSupportElementForSolution;
1345 double side;
1346
2/2
✓ Branch 0 taken 160 times.
✓ Branch 1 taken 160 times.
 320 sideOfSupportElementForSolution = iSup2 == 0 ? LEFT : RIGHT;
1347
2/2
✓ Branch 0 taken 160 times.
✓ Branch 1 taken 160 times.
 320 side = sideOfSupportElementForSolution == LEFT ? 1. : -1.;
1348
1349 // std::set the extrapolation of velocity
1350
5/10
✓ Branch 1 taken 320 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 320 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 320 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 320 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 320 times.
✗ Branch 11 not taken.
 320 if(FelisceParam::instance().NSequationFlag == 1 && FelisceParam::instance().useInterfaceFlowStab) {
1351 // Get the support element of the old solutions
1352
1/2
✓ Branch 3 taken 320 times.
✗ Branch 4 not taken.
 320 std::vector<felInt> ielOld(m_solutionOld.size());
1353
1/2
✓ Branch 3 taken 320 times.
✗ Branch 4 not taken.
 320 std::vector<felInt> ielOldDup(m_solutionOld.size());
1354
1/2
✓ Branch 3 taken 320 times.
✗ Branch 4 not taken.
 320 std::vector< std::vector<felInt> > oldSupportElement(m_solutionOld.size());
1355
2/2
✓ Branch 1 taken 640 times.
✓ Branch 2 taken 320 times.
 960 for(std::size_t i=0; i<oldSupportElement.size(); ++i) {
1356
1/2
✓ Branch 4 taken 640 times.
✗ Branch 5 not taken.
 640 m_supportDofUnknownOld[i][0].getIdElementSupport(idElem, oldSupportElement[i]);
1357
1358 // std::set the id for the old solution
1359
1/2
✓ Branch 2 taken 640 times.
✗ Branch 3 not taken.
 640 if(oldSupportElement[i].size() > 1) {
1360 // this element is also duplicated for the old solution
1361
2/2
✓ Branch 0 taken 320 times.
✓ Branch 1 taken 320 times.
 640 if(sideOfSupportElementForSolution == LEFT) {
1362 320 ielOld[i] = oldSupportElement[i][0];
1363 320 ielOldDup[i] = oldSupportElement[i][1];
1364 } else {
1365 320 ielOld[i] = oldSupportElement[i][1];
1366 320 ielOldDup[i] = oldSupportElement[i][0];
1367 }
1368 } else {
1369 // this element is NOT duplicated for the old solution
1370 ✗ ielOld[i] = oldSupportElement[i][0];
1371 ✗ ielOldDup[i] = oldSupportElement[i][0];
1372 }
1373 }
1374
1375
1376
1/2
✓ Branch 1 taken 320 times.
✗ Branch 2 not taken.
 320 m_elemFieldAdv.val.clear();
1377
1/2
✓ Branch 1 taken 320 times.
✗ Branch 2 not taken.
 320 m_elemFieldAdvDup.val.clear();
1378
2/4
✓ Branch 1 taken 320 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 320 times.
 320 if(FelisceParam::instance().useProjectionForNitscheXFEM) {
1379 ✗ for(felInt i=0; i<FelisceParam::instance().orderBdfNS; ++i) {
1380 ✗ m_elemFieldAdvTmp.setValue(m_seqVelExtrapol[i], *m_feVel, iel1, m_iVelocity, m_ao, dof());
1381 ✗ m_elemFieldAdv.val += m_elemFieldAdvTmp.val;
1382 ✗ m_elemFieldAdvTmp.setValue(m_seqVelExtrapol[i], *m_feVel, iel2, m_iVelocity, m_ao, dof());
1383 ✗ m_elemFieldAdvDup.val += m_elemFieldAdvTmp.val;
1384 }
1385 } else {
1386
3/4
✓ Branch 1 taken 640 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 320 times.
✓ Branch 4 taken 320 times.
 640 for(felInt i=0; i<FelisceParam::instance().orderBdfNS; ++i) {
1387
1/2
✓ Branch 5 taken 320 times.
✗ Branch 6 not taken.
 320 m_elemFieldAdvTmp.setValue(m_seqVelExtrapol[i], *m_feVel, ielOld[i], m_iVelocity, m_aoOld[i], m_dofOld[i]);
1388
1/2
✓ Branch 1 taken 320 times.
✗ Branch 2 not taken.
 320 m_elemFieldAdv.val += m_elemFieldAdvTmp.val;
1389
1/2
✓ Branch 5 taken 320 times.
✗ Branch 6 not taken.
 320 m_elemFieldAdvTmp.setValue(m_seqVelExtrapol[i], *m_feVel, ielOldDup[i], m_iVelocity, m_aoOld[i], m_dofOld[i]);
1390
1/2
✓ Branch 1 taken 320 times.
✗ Branch 2 not taken.
 320 m_elemFieldAdvDup.val += m_elemFieldAdvTmp.val;
1391 }
1392 }
1393 320 }
1394
1395
2/2
✓ Branch 0 taken 160 times.
✓ Branch 1 taken 160 times.
 320 if(iSup1 == iSup2) {
1396 // enforce the continuity
1397 // Nitsche's penalty term, \gamma \nu / h u_i . v_i
1398
1/2
✓ Branch 4 taken 160 times.
✗ Branch 5 not taken.
 160 m_elementMat[0]->phi_i_phi_j(penaltyParam, *m_feVel, 0, 0, dimension());
1399
1400 // Nitsche's consistency terms, -\sigma(u_i, p_i)n_i . v_i - \sigma(v_i, q_i)n_i . u_i
1401 // u_i & v_i
1402
1/2
✓ Branch 1 taken 160 times.
✗ Branch 2 not taken.
 160 computeNitscheSigma_u_v(0.5 * side * m_viscosity, -0.5 * side);
1403
1/2
✓ Branch 1 taken 160 times.
✗ Branch 2 not taken.
 160 computeNitscheSigma_v_u(0.5 * side * m_viscosity, -0.5 * side);
1404
1405 // stabilize inflow
1406
5/10
✓ Branch 1 taken 160 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 160 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 160 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 160 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 160 times.
✗ Branch 11 not taken.
 160 if(FelisceParam::instance().NSequationFlag == 1 && FelisceParam::instance().useInterfaceFlowStab){
1407 // (0.5 \rho_f \beta_1 \cdot n_2) (u_1 \cdot v_1)-(0.5 \rho_f \beta_2 \cdot n_2) (u_2 \cdot v_2)
1408
2/4
✓ Branch 4 taken 160 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 160 times.
✗ Branch 8 not taken.
 160 m_elementMat[0]->f_dot_n_phi_i_phi_j(0.5 * side * FelisceParam::instance().density, *m_feVel, m_elemFieldAdv, m_elemFieldNormal, 0, 0, dimension(), 0);
1409 }
1410 } else {
1411
1/2
✓ Branch 4 taken 160 times.
✗ Branch 5 not taken.
 160 m_elementMat[0]->phi_i_phi_j(-penaltyParam, *m_feVel, 0, 0, dimension());
1412
1413
1/2
✓ Branch 1 taken 160 times.
✗ Branch 2 not taken.
 160 computeNitscheSigma_u_v(-0.5 * side * m_viscosity, 0.5 * side);
1414
1/2
✓ Branch 1 taken 160 times.
✗ Branch 2 not taken.
 160 computeNitscheSigma_v_u( 0.5 * side * m_viscosity, -0.5 * side);
1415
1416 // stabilize inflow
1417
5/10
✓ Branch 1 taken 160 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 160 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 160 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 160 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 160 times.
✗ Branch 11 not taken.
 160 if(FelisceParam::instance().NSequationFlag == 1 && FelisceParam::instance().useInterfaceFlowStab) {
1418
2/4
✓ Branch 4 taken 160 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 160 times.
✗ Branch 8 not taken.
 160 m_elementMat[0]->f_dot_n_phi_i_phi_j(0.25 * side * FelisceParam::instance().density, *m_feVel, m_elemFieldAdv, m_elemFieldNormal, 0, 0, dimension(), 0);
1419
2/4
✓ Branch 4 taken 160 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 160 times.
✗ Branch 8 not taken.
 160 m_elementMat[0]->f_dot_n_phi_i_phi_j(0.25 * side * FelisceParam::instance().density, *m_feVel, m_elemFieldAdvDup, m_elemFieldNormal, 0, 0, dimension(), 0);
1420 }
1421 }
1422
1423 // add values of elemMat in the global matrix: matrix(0).
1424
1/2
✓ Branch 3 taken 320 times.
✗ Branch 4 not taken.
 320 setValueMatrixRHS(vectorIdSupport[iSup1], vectorIdSupport[iSup2], flag);
1425 }
1426 }
1427 }
1428
1429 // desallocate array use for assemble with petsc.
1430
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 desallocateArrayForAssembleMatrixRHS(flag);
1431 40 }
1432 40 }
1433
1434
1435
1436
1437 40 void LinearProblemFSINitscheXFEMFluid::assembleMatrixGhostPenalty() {
1438 // ------------- //
1439 // Ghost penalty //
1440 // ------------- //
1441
1442 // get the types of elements
1443 40 GeometricMeshRegion::ElementType eltType = m_meshLocal[m_currentMesh]->bagElementTypeDomain()[0];
1444 40 felInt numEdgePerType = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second->numEdge();
1445 40 felInt numPointPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];;
1446
1/2
✓ Branch 2 taken 40 times.
✗ Branch 3 not taken.
 40 std::vector<felInt> idOfFacesCurrent(numEdgePerType, -1);
1447
1/2
✓ Branch 2 taken 40 times.
✗ Branch 3 not taken.
 40 std::vector<felInt> idOfFacesOpposite(numEdgePerType, -1);
1448
1449 // define finite element
1450
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 defineFiniteElement(eltType);
1451
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 initElementArray();
1452
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 defineCurrentFiniteElementWithBd(eltType);
1453
1454 // allocate arrays for assembling the matrix
1455 40 FlagMatrixRHS flag = FlagMatrixRHS::only_matrix;
1456
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 allocateArrayForAssembleMatrixRHS(flag);
1457
1458 // init variables
1459
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 initPerElementType(eltType, flag);
1460 40 m_feVelWithBd = m_listCurrentFiniteElementWithBd[m_iVelocity];
1461
1462 // get global supportDofMesh informations
1463 40 std::vector<felInt>& iEle = m_supportDofUnknown[0].iEle();
1464 40 std::vector<felInt>& iSupportDof = m_supportDofUnknown[0].iSupportDof();
1465
1466 // element informations
1467 40 std::pair<bool, GeometricMeshRegion::ElementType> adjElt;
1468 40 std::vector<felInt> ielCurrentSupportElt;
1469 40 std::vector<felInt> ielOppositeSupportElt;
1470
1/2
✓ Branch 2 taken 40 times.
✗ Branch 3 not taken.
 40 std::vector<Point*> elemCurrentPoint(m_numVerPerFluidElt);
1471
1/2
✓ Branch 2 taken 40 times.
✗ Branch 3 not taken.
 40 std::vector<felInt> elemCurrentIdPoint(m_numVerPerFluidElt);
1472
1/2
✓ Branch 2 taken 40 times.
✗ Branch 3 not taken.
 40 std::vector<Point*> elemOppositePoint(m_numVerPerFluidElt);
1473
1/2
✓ Branch 2 taken 40 times.
✗ Branch 3 not taken.
 40 std::vector<felInt> elemOppositeIdPoint(m_numVerPerFluidElt);
1474 felInt currentIelGeo;
1475 felInt currentGlobalIelGeo;
1476 felInt oppositeGlobalIelGeo;
1477 felInt idCurrentSupElt;
1478 felInt idOppositeSupElt;
1479
1480 // loop over all the intersected elements
1481
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 auto& idInterElt = m_duplicateSupportElements->getIntersectedMshElt();
1482
6/8
✓ Branch 4 taken 1456 times.
✗ Branch 5 not taken.
✓ Branch 10 taken 1416 times.
✓ Branch 11 taken 40 times.
✓ Branch 12 taken 1456 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 1416 times.
✓ Branch 16 taken 40 times.
 1456 for (auto ielIt = idInterElt.begin(); ielIt != idInterElt.end() && ielIt->first < meshLocal()->numElements(eltType); ++ielIt) {
1483 // get the geometric id of the element and the ids of its support elements
1484 1416 currentIelGeo = ielIt->first;
1485
1/2
✓ Branch 2 taken 1416 times.
✗ Branch 3 not taken.
 1416 ISLocalToGlobalMappingApply(m_mappingElem[m_currentMesh], 1, &currentIelGeo, &currentGlobalIelGeo);
1486
1487 // get all the faces of the element
1488
1/2
✓ Branch 1 taken 1416 times.
✗ Branch 2 not taken.
 1416 if(dimension() == 2)
1489
1/2
✓ Branch 3 taken 1416 times.
✗ Branch 4 not taken.
 1416 m_mesh[m_currentMesh]->getAllEdgeOfElement(currentGlobalIelGeo, idOfFacesCurrent);
1490 else
1491 ✗ m_mesh[m_currentMesh]->getAllFaceOfElement(currentGlobalIelGeo, idOfFacesCurrent);
1492
1493 // get informations on the current element
1494
1/2
✓ Branch 1 taken 1416 times.
✗ Branch 2 not taken.
 1416 setElemPoint(eltType, currentIelGeo, elemCurrentPoint, elemCurrentIdPoint, ielCurrentSupportElt);
1495
1496 // for each support elements
1497
2/2
✓ Branch 1 taken 2832 times.
✓ Branch 2 taken 1416 times.
 4248 for(std::size_t iSup=0; iSup<ielCurrentSupportElt.size(); ++iSup) {
1498
2/2
✓ Branch 0 taken 1416 times.
✓ Branch 1 taken 1416 times.
 2832 sideOfInterface side = iSup == 0 ? LEFT : RIGHT;
1499 2832 idCurrentSupElt = ielCurrentSupportElt[iSup];
1500
1501 // loop over all the edges of the current support element
1502
2/2
✓ Branch 1 taken 8496 times.
✓ Branch 2 taken 2832 times.
 11328 for(std::size_t iface=0; iface < idOfFacesCurrent.size(); ++iface) {
1503 8496 bool computeGhostPenaltyOnThisEdge = false;
1504 8496 bool computeBothSide = false;
1505 8496 oppositeGlobalIelGeo = currentGlobalIelGeo;
1506
1507 // get the neighboor. If it exists, this is an inner edge
1508 // first = true if there is an adjacent element
1509 // second = the ElementType of the adjacent element if first is true
1510
1/2
✓ Branch 3 taken 8496 times.
✗ Branch 4 not taken.
 8496 adjElt = m_mesh[m_currentMesh]->getAdjElement(oppositeGlobalIelGeo, iface);
1511
1512
2/2
✓ Branch 0 taken 8336 times.
✓ Branch 1 taken 160 times.
 8496 if(adjElt.first) {
1513 // This edge is an inner edge
1514 // get information on the opposite element
1515
1/2
✓ Branch 2 taken 8336 times.
✗ Branch 3 not taken.
 8336 m_supportDofUnknown[0].getIdElementSupport(oppositeGlobalIelGeo, ielOppositeSupportElt);
1516
1/2
✓ Branch 3 taken 8336 times.
✗ Branch 4 not taken.
 8336 m_mesh[m_currentMesh]->getOneElement(oppositeGlobalIelGeo, elemOppositeIdPoint);
1517
2/2
✓ Branch 0 taken 25008 times.
✓ Branch 1 taken 8336 times.
 33344 for (felInt iPoint=0; iPoint<numPointPerElt; ++iPoint)
1518 25008 elemOppositePoint[iPoint] = &(m_mesh[m_currentMesh]->listPoints()[elemOppositeIdPoint[iPoint]]);
1519
1520
1521 8336 idOppositeSupElt = ielOppositeSupportElt[0];
1522
2/2
✓ Branch 1 taken 5344 times.
✓ Branch 2 taken 2992 times.
 8336 if(ielOppositeSupportElt.size() == 2) {
1523 // the other element is also duplicated, take the one on the same side
1524 5344 computeGhostPenaltyOnThisEdge = true;
1525
2/2
✓ Branch 0 taken 2672 times.
✓ Branch 1 taken 2672 times.
 5344 idOppositeSupElt = side == LEFT ? ielOppositeSupportElt[0] : ielOppositeSupportElt[1];
1526 } else {
1527 // the other element is not duplicated, nothing to do if it is not on the same side
1528 2992 std::size_t sizeCurrent = m_supportDofUnknown[0].getNumSupportDof(idCurrentSupElt);
1529 2992 std::size_t sizeOpposite = m_supportDofUnknown[0].getNumSupportDof(idOppositeSupElt);
1530
1531 2992 felInt countSupportDofInCommon = 0;
1532
2/2
✓ Branch 0 taken 8976 times.
✓ Branch 1 taken 2992 times.
 11968 for(std::size_t iSup1=0; iSup1<sizeCurrent; ++iSup1) {
1533
2/2
✓ Branch 0 taken 26928 times.
✓ Branch 1 taken 8976 times.
 35904 for(std::size_t iSup2=0; iSup2<sizeOpposite; ++iSup2) {
1534
2/2
✓ Branch 5 taken 3152 times.
✓ Branch 6 taken 23776 times.
 26928 if(iSupportDof[iEle[ielOppositeSupportElt[0]] + iSup2] == iSupportDof[iEle[idCurrentSupElt] + iSup1]) {
1535 3152 ++countSupportDofInCommon;
1536 }
1537 }
1538 }
1539
1540
2/2
✓ Branch 0 taken 1496 times.
✓ Branch 1 taken 1496 times.
 2992 if(countSupportDofInCommon > 1) {
1541 // the elements are on the same side
1542 1496 idOppositeSupElt = ielOppositeSupportElt[0];
1543 1496 computeGhostPenaltyOnThisEdge = true;
1544 1496 computeBothSide = true;
1545 }
1546 }
1547
1548 // add the contribution of this edge
1549
2/2
✓ Branch 0 taken 6840 times.
✓ Branch 1 taken 1496 times.
 8336 if(computeGhostPenaltyOnThisEdge) {
1550 // finite element with boundary for the opposite element
1551 6840 const GeoElement *geoEle = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second;
1552 6840 felInt typeOfFiniteElement = m_velocity->finiteElementType();
1553
1/2
✓ Branch 3 taken 6840 times.
✗ Branch 4 not taken.
 6840 const RefElement *refEle = geoEle->defineFiniteEle(eltType, typeOfFiniteElement, *m_meshLocal[m_currentMesh]);
1554
1/2
✓ Branch 3 taken 6840 times.
✗ Branch 4 not taken.
 6840 CurrentFiniteElementWithBd oppositeEltFEWithBd(*refEle, *geoEle, m_velocity->degreeOfExactness(), m_velocity->degreeOfExactness());
1555
1556
1557 // update the finite elements
1558
1/2
✓ Branch 1 taken 6840 times.
✗ Branch 2 not taken.
 6840 m_feVelWithBd->updateFirstDeriv(0, elemCurrentPoint);
1559
1/2
✓ Branch 1 taken 6840 times.
✗ Branch 2 not taken.
 6840 m_feVelWithBd->updateBdMeasNormal();
1560
1/2
✓ Branch 1 taken 6840 times.
✗ Branch 2 not taken.
 6840 oppositeEltFEWithBd.updateFirstDeriv(0, elemOppositePoint);
1561
1/2
✓ Branch 1 taken 6840 times.
✗ Branch 2 not taken.
 6840 oppositeEltFEWithBd.updateBdMeasNormal();
1562
1563
1564 // find the local id of the edge in the opposite element
1565 6840 felInt localIdFaceOpposite = -1;
1566
1/2
✓ Branch 1 taken 6840 times.
✗ Branch 2 not taken.
 6840 if(dimension() == 2)
1567
1/2
✓ Branch 3 taken 6840 times.
✗ Branch 4 not taken.
 6840 m_mesh[m_currentMesh]->getAllEdgeOfElement(eltType, oppositeGlobalIelGeo, idOfFacesOpposite);
1568 else
1569 ✗ m_mesh[m_currentMesh]->getAllFaceOfElement(eltType, oppositeGlobalIelGeo, idOfFacesOpposite);
1570
1571
2/2
✓ Branch 1 taken 20520 times.
✓ Branch 2 taken 6840 times.
 27360 for(std::size_t jface=0; jface<idOfFacesOpposite.size(); ++jface) {
1572
2/2
✓ Branch 2 taken 6840 times.
✓ Branch 3 taken 13680 times.
 20520 if(idOfFacesCurrent[iface] == idOfFacesOpposite[jface]) {
1573 6840 localIdFaceOpposite = jface;
1574 }
1575 }
1576
1577 // std::set coefficient of the ghost penalty
1578
3/6
✓ Branch 1 taken 6840 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6840 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6840 times.
✗ Branch 8 not taken.
 6840 double gammag = FelisceParam::instance().coefGhostPenalty * m_viscosity * m_feVelWithBd->bdEle(iface).diameter();
1579
1580 // ----------------------------------- //
1581 // current element for phi_i and phi_j //
1582 // ----------------------------------- //
1583
1/2
✓ Branch 3 taken 6840 times.
✗ Branch 4 not taken.
 6840 m_elementMat[0]->zero();
1584
1/2
✓ Branch 4 taken 6840 times.
✗ Branch 5 not taken.
 6840 m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n(gammag, *m_feVelWithBd, *m_feVelWithBd, iface, iface, 0, 0, dimension());
1585
1/2
✓ Branch 1 taken 6840 times.
✗ Branch 2 not taken.
 6840 setValueMatrixRHS(idCurrentSupElt, idCurrentSupElt, flag);
1586
1587 // -------------------------------------------------------- //
1588 // current element for phi_i and opposite element for phi_j //
1589 // -------------------------------------------------------- //
1590
1/2
✓ Branch 3 taken 6840 times.
✗ Branch 4 not taken.
 6840 m_elementMat[0]->zero();
1591
1/2
✓ Branch 4 taken 6840 times.
✗ Branch 5 not taken.
 6840 m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n(gammag, oppositeEltFEWithBd, *m_feVelWithBd, localIdFaceOpposite, iface, 0, 0, dimension());
1592
1/2
✓ Branch 1 taken 6840 times.
✗ Branch 2 not taken.
 6840 setValueMatrixRHS(idCurrentSupElt, idOppositeSupElt, flag);
1593
1594 // if the opposite element is not duplicated, it will not be handle in this loop
1595 // compute all the term in this case
1596
2/2
✓ Branch 0 taken 1496 times.
✓ Branch 1 taken 5344 times.
 6840 if(computeBothSide) {
1597 // -------------------------------------------------------- //
1598 // opposite element for phi_i and current element for phi_j //
1599 // -------------------------------------------------------- //
1600
1/2
✓ Branch 3 taken 1496 times.
✗ Branch 4 not taken.
 1496 m_elementMat[0]->zero();
1601
1/2
✓ Branch 4 taken 1496 times.
✗ Branch 5 not taken.
 1496 m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n(gammag, *m_feVelWithBd, oppositeEltFEWithBd, iface, localIdFaceOpposite, 0, 0, dimension());
1602
1/2
✓ Branch 1 taken 1496 times.
✗ Branch 2 not taken.
 1496 setValueMatrixRHS(idOppositeSupElt, idCurrentSupElt, flag);
1603
1604 // ------------------------------------ //
1605 // opposite element for phi_i and phi_j //
1606 // ------------------------------------ //
1607
1/2
✓ Branch 3 taken 1496 times.
✗ Branch 4 not taken.
 1496 m_elementMat[0]->zero();
1608
1/2
✓ Branch 4 taken 1496 times.
✗ Branch 5 not taken.
 1496 m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n(gammag, oppositeEltFEWithBd, oppositeEltFEWithBd, localIdFaceOpposite, localIdFaceOpposite, 0, 0, dimension());
1609
1/2
✓ Branch 1 taken 1496 times.
✗ Branch 2 not taken.
 1496 setValueMatrixRHS(idOppositeSupElt, idOppositeSupElt, flag);
1610 }
1611 6840 }
1612 }
1613 }
1614 }
1615 }
1616
1617 // desallocate array use for assemble with petsc.
1618
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 desallocateArrayForAssembleMatrixRHS(flag);
1619 40 }
1620
1621
1622
1623 40 void LinearProblemFSINitscheXFEMFluid::assembleMatrixFaceOriented() {
1624 // Use Face-oriented stabilization
1625 // E. Burman, M. A. Fernandez and P. Hansbo, Continuous interior penalty finite element method for
1626 // Oseen's equations, 2006.
1627 40 std::pair<bool, GeometricMeshRegion::ElementType> adjElt;
1628
1629 GeometricMeshRegion::ElementType eltType; // Type of element
1630 GeometricMeshRegion::ElementType eltTypeOpp; // Type of element
1631
1632 bool addOppositeSupportElt; // true if we need to compute the stabilization on this face
1633
1634 40 felInt ielCurrentLocalGeo = 0; // local geometric id of the current element
1635 40 felInt ielCurrentGlobalGeo = 0; // global geometric id of the current element
1636 40 felInt ielOppositeGlobalGeo = 0; // global geometric id of the opposite element
1637 felInt numFacesPerType; // number of faces of an element of a given type
1638 felInt numEltPerLabel; // number of element by reference
1639 felInt numPointPerElt; // number of vertices by element
1640 felInt currentSupport; // id of the current support element
1641 felInt oppositeSupport; // id of the opposite support element
1642 40 std::vector<felInt> ielCurrentGlobalSup; // local support element id of the current element
1643 40 std::vector<felInt> ielOppositeGlobalSup; // local support element id of the opposite element
1644
1645 40 std::vector<felInt> idOfFacesCurrent; // ids of the current element edges
1646 40 std::vector<felInt> idOfFacesOpposite; // ids of the opposite element edges
1647 40 std::vector<felInt> currentElemIdPoint; // ids of the vertices of the current element
1648 40 std::vector<felInt> oppositeElemIdPoint; // ids of the vertices of the opposite element
1649 40 std::vector<Point*> currentElemPoint; // point coordinates of the current element vertices
1650 40 std::vector<Point*> oppositeElemPoint; // point coordinates of the opposite element vertices
1651 40 std::vector<felInt> countElement; // count the element by type
1652
1653 40 FlagMatrixRHS flag = FlagMatrixRHS::only_matrix; // flag to only assemble the matrix
1654
1655
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 ElementField elemFieldAdvFace;
1656
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 ElementField elemFieldAdvFaceTmp;
1657
1658 40 const std::vector<felInt>& iSupportDof = m_supportDofUnknown[0].iSupportDof();
1659 40 const std::vector<felInt>& iEle = m_supportDofUnknown[0].iEle();
1660
1661 // resize std::vector
1662
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 countElement.resize(GeometricMeshRegion::m_numTypesOfElement, 0);
1663
1664 // volume bag
1665 40 const std::vector<ElementType>& bagElementTypeDomain = m_meshLocal[m_currentMesh]->bagElementTypeDomain();
1666
1667
1668 // finite element with bd for the opposite element
1669 40 std::unordered_map<GeometricMeshRegion::ElementType, CurrentFiniteElementWithBd*> oppositeEltFEWithBdVel;
1670 40 std::unordered_map<GeometricMeshRegion::ElementType, CurrentFiniteElementWithBd*> oppositeEltFEWithBdPre;
1671
1672
2/2
✓ Branch 1 taken 40 times.
✓ Branch 2 taken 40 times.
 80 for (std::size_t i=0; i<bagElementTypeDomain.size(); ++i) {
1673 40 eltType = bagElementTypeDomain[i];
1674
1675
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 m_velocity = &m_listVariable[m_listVariable.getVariableIdList(velocity)];
1676
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 m_pressure = &m_listVariable[m_listVariable.getVariableIdList(pressure)];
1677 40 const GeoElement *geoEle = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second;
1678
1679 40 felInt typeOfFEVel = m_velocity->finiteElementType();
1680
1/2
✓ Branch 3 taken 40 times.
✗ Branch 4 not taken.
 40 const RefElement *refEleVel = geoEle->defineFiniteEle(eltType, typeOfFEVel, *m_mesh[m_currentMesh]);
1681
3/6
✓ Branch 3 taken 40 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 40 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 40 times.
✗ Branch 10 not taken.
 40 oppositeEltFEWithBdVel[eltType] = new CurrentFiniteElementWithBd(*refEleVel, *geoEle, m_velocity->degreeOfExactness(), m_velocity->degreeOfExactness());
1682
1683 40 felInt typeOfFEPre = m_pressure->finiteElementType();
1684
1/2
✓ Branch 3 taken 40 times.
✗ Branch 4 not taken.
 40 const RefElement *refElePre = geoEle->defineFiniteEle(eltType, typeOfFEPre, *m_mesh[m_currentMesh]);
1685
3/6
✓ Branch 3 taken 40 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 40 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 40 times.
✗ Branch 10 not taken.
 40 oppositeEltFEWithBdPre[eltType] = new CurrentFiniteElementWithBd(*refElePre, *geoEle, m_pressure->degreeOfExactness(), m_pressure->degreeOfExactness());
1686 }
1687
1688
1689 // ----------------------------
1690 // First loop on geometric type
1691
2/2
✓ Branch 1 taken 40 times.
✓ Branch 2 taken 40 times.
 80 for (std::size_t itype=0; itype<bagElementTypeDomain.size(); ++itype) {
1692 // get element type and number of vertices and "faces"
1693 40 eltType = bagElementTypeDomain[itype];
1694 40 numPointPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
1695 40 numFacesPerType = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second->numBdEle();
1696
1697 // resize of vectors
1698
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 idOfFacesCurrent.resize(numFacesPerType, -1);
1699
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 idOfFacesOpposite.resize(numFacesPerType, -1);
1700
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 currentElemIdPoint.resize(numPointPerElt, -1);
1701
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 oppositeElemIdPoint.resize(numPointPerElt, -1);
1702
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 currentElemPoint.resize(numPointPerElt, nullptr);
1703
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 oppositeElemPoint.resize(numPointPerElt, nullptr);
1704
1705 // define finite element
1706
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 defineFiniteElement(eltType);
1707
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 initElementArray();
1708
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 defineCurrentFiniteElementWithBd(eltType);
1709
1710 // allocate arrays for assembling the matrix
1711
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 allocateArrayForAssembleMatrixRHS(flag);
1712
1713 // init variables
1714
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 initPerElementType(eltType, flag);
1715 40 m_feVelWithBd = m_listCurrentFiniteElementWithBd[m_iVelocity];
1716 40 m_fePreWithBd = m_listCurrentFiniteElementWithBd[m_iPressure];
1717
1718
1719 // ---------------------------------
1720 // second loop on region of the mesh
1721
2/2
✓ Branch 8 taken 40 times.
✓ Branch 9 taken 40 times.
 80 for(auto itRef = m_meshLocal[m_currentMesh]->intRefToBegEndMaps[eltType].begin(); itRef != m_meshLocal[m_currentMesh]->intRefToBegEndMaps[eltType].end(); itRef++) {
1722 // get the number of element with this reference
1723 40 numEltPerLabel = itRef->second.second;
1724
1725 // std::set the current label of LinearProblem to the current one of this loop
1726
1/2
✓ Branch 2 taken 40 times.
✗ Branch 3 not taken.
 40 initPerDomain(itRef->first, flag);
1727
1728 // ---------------------
1729 // third loop on element
1730
2/2
✓ Branch 0 taken 136880 times.
✓ Branch 1 taken 40 times.
 136920 for(felInt iel=0; iel<numEltPerLabel; ++iel) {
1731 // get the global id of the current geometric element
1732
1/2
✓ Branch 2 taken 136880 times.
✗ Branch 3 not taken.
 136880 ISLocalToGlobalMappingApply(m_mappingElem[m_currentMesh], 1, &ielCurrentLocalGeo, &ielCurrentGlobalGeo);
1733
1734 // get all the faces of the element
1735
1/2
✓ Branch 1 taken 136880 times.
✗ Branch 2 not taken.
 136880 if(dimension() == 2)
1736
1/2
✓ Branch 3 taken 136880 times.
✗ Branch 4 not taken.
 136880 m_mesh[m_currentMesh]->getAllEdgeOfElement(ielCurrentGlobalGeo, idOfFacesCurrent);
1737 else
1738 ✗ m_mesh[m_currentMesh]->getAllFaceOfElement(ielCurrentGlobalGeo, idOfFacesCurrent);
1739
1740 // get informations on the current element
1741
1/2
✓ Branch 2 taken 136880 times.
✗ Branch 3 not taken.
 136880 setElemPoint(eltType, countElement[eltType], currentElemPoint, currentElemIdPoint, ielCurrentGlobalSup);
1742
1743 // update the finite elements
1744
1/2
✓ Branch 1 taken 136880 times.
✗ Branch 2 not taken.
 136880 m_feVelWithBd->updateFirstDeriv(0, currentElemPoint);
1745
1/2
✓ Branch 1 taken 136880 times.
✗ Branch 2 not taken.
 136880 m_feVelWithBd->updateBdMeasNormal();
1746
1/2
✓ Branch 1 taken 136880 times.
✗ Branch 2 not taken.
 136880 m_fePreWithBd->updateFirstDeriv(0, currentElemPoint);
1747
1/2
✓ Branch 1 taken 136880 times.
✗ Branch 2 not taken.
 136880 m_fePreWithBd->updateBdMeasNormal();
1748
1749
1750 // loop over the support element of the current geometric element
1751
2/2
✓ Branch 1 taken 138296 times.
✓ Branch 2 taken 136880 times.
 275176 for(std::size_t ielSup=0; ielSup<ielCurrentGlobalSup.size(); ++ielSup) {
1752 138296 currentSupport = ielCurrentGlobalSup[ielSup];
1753
1754 // get the old support dof
1755
2/4
✓ Branch 2 taken 138296 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 138296 times.
✗ Branch 6 not taken.
 138296 std::vector< std::vector<felInt> > oldSupportElement(FelisceParam::instance().orderBdfNS);
1756
2/4
✓ Branch 2 taken 138296 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 138296 times.
✗ Branch 6 not taken.
 138296 std::vector<felInt> ielOld(FelisceParam::instance().orderBdfNS);
1757
1758
3/4
✓ Branch 1 taken 276592 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 138296 times.
✓ Branch 4 taken 138296 times.
 276592 for(felInt iextrap=0; iextrap<FelisceParam::instance().orderBdfNS; ++iextrap) {
1759 // get the old support element
1760
1/2
✓ Branch 4 taken 138296 times.
✗ Branch 5 not taken.
 138296 m_supportDofUnknownOld[iextrap][0].getIdElementSupport(ielCurrentGlobalGeo, oldSupportElement[iextrap]);
1761
1762 // std::set the id for the old solution
1763
2/2
✓ Branch 2 taken 2792 times.
✓ Branch 3 taken 135504 times.
 138296 if(oldSupportElement[iextrap].size() > 1)
1764 2792 ielOld[iextrap] = oldSupportElement[iextrap][ielSup];
1765 else
1766 135504 ielOld[iextrap] = oldSupportElement[iextrap][0];
1767 }
1768
1769 // compute coefficients
1770
1/2
✓ Branch 1 taken 138296 times.
✗ Branch 2 not taken.
 138296 m_elemFieldAdv.val.clear();
1771
2/4
✓ Branch 1 taken 138296 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 138296 times.
 138296 if(FelisceParam::instance().useProjectionForNitscheXFEM) {
1772 ✗ for(felInt iextrap=0; iextrap<FelisceParam::instance().orderBdfNS; ++iextrap) {
1773 ✗ m_elemFieldAdvTmp.setValue(m_seqVelExtrapol[iextrap], *m_feVel, currentSupport, m_iVelocity, m_ao, dof());
1774 ✗ m_elemFieldAdv.val += m_elemFieldAdvTmp.val;
1775 }
1776 } else {
1777
3/4
✓ Branch 1 taken 276592 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 138296 times.
✓ Branch 4 taken 138296 times.
 276592 for(felInt iextrap=0; iextrap<FelisceParam::instance().orderBdfNS; ++iextrap) {
1778
1/2
✓ Branch 5 taken 138296 times.
✗ Branch 6 not taken.
 138296 m_elemFieldAdvTmp.setValue(m_seqVelExtrapol[iextrap], *m_feVel, ielOld[iextrap], m_iVelocity, m_aoOld[iextrap], m_dofOld[iextrap]);
1779
1/2
✓ Branch 1 taken 138296 times.
✗ Branch 2 not taken.
 138296 m_elemFieldAdv.val += m_elemFieldAdvTmp.val;
1780 }
1781 }
1782
1783 138296 double normVelExtrap = 0;
1784
2/2
✓ Branch 1 taken 414888 times.
✓ Branch 2 taken 138296 times.
 553184 for(std::size_t jcol=0; jcol<m_elemFieldAdv.val.size2(); ++jcol)
1785
2/4
✓ Branch 1 taken 414888 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 414888 times.
✗ Branch 5 not taken.
 414888 normVelExtrap = std::max(normVelExtrap, norm_2(column(m_elemFieldAdv.val, jcol)));
1786
1787
1/2
✓ Branch 1 taken 138296 times.
✗ Branch 2 not taken.
 138296 double hK = m_feVelWithBd->diameter();
1788 138296 double ReK = m_density * normVelExtrap * hK / m_viscosity;
1789
1/2
✓ Branch 1 taken 138296 times.
✗ Branch 2 not taken.
 138296 double coeffVel = FelisceParam::instance().stabFOVel * std::min(1., ReK) * hK * hK;
1790
1791 138296 double gammaP = 0.;
1792
2/4
✓ Branch 1 taken 138296 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 138296 times.
 138296 if(FelisceParam::instance().NSequationFlag == 0) {
1793 ✗ gammaP = FelisceParam::instance().stabFOPre * hK * hK * hK / m_viscosity;
1794 } else {
1795
2/2
✓ Branch 0 taken 55048 times.
✓ Branch 1 taken 83248 times.
 138296 double xiPre = ReK < 1 ? m_density * hK * hK * hK / m_viscosity : hK * hK / normVelExtrap;
1796
1/2
✓ Branch 1 taken 138296 times.
✗ Branch 2 not taken.
 138296 gammaP = FelisceParam::instance().stabFOPre * xiPre;
1797 }
1798
1799 // --------------------------------------------------------
1800 // loop over all the faces of the current geometric element
1801
2/2
✓ Branch 1 taken 414888 times.
✓ Branch 2 taken 138296 times.
 553184 for(std::size_t iface=0; iface<idOfFacesCurrent.size(); ++iface) {
1802 // check if this face is an inner face or a boundary
1803 414888 ielOppositeGlobalGeo = ielCurrentGlobalGeo;
1804
1805 // first = true if there is an adjacent element
1806 // second = the ElementType of the adjacent element if first is true
1807
1/2
✓ Branch 3 taken 414888 times.
✗ Branch 4 not taken.
 414888 adjElt =m_mesh[m_currentMesh]->getAdjElement(ielOppositeGlobalGeo, iface);
1808
1809 // compute the stabilization only on inner faces
1810
2/2
✓ Branch 0 taken 406968 times.
✓ Branch 1 taken 7920 times.
 414888 if(adjElt.first) {
1811 // get the type of the opposite element
1812 406968 eltTypeOpp = adjElt.second;
1813
1814 // find the local id of the edge in the opposite element
1815 406968 felInt localIdFaceOpposite = -1;
1816
1/2
✓ Branch 1 taken 406968 times.
✗ Branch 2 not taken.
 406968 if(dimension() == 2)
1817
1/2
✓ Branch 3 taken 406968 times.
✗ Branch 4 not taken.
 406968 m_mesh[m_currentMesh]->getAllEdgeOfElement(ielOppositeGlobalGeo, idOfFacesOpposite);
1818 else
1819 ✗ m_mesh[m_currentMesh]->getAllFaceOfElement(ielOppositeGlobalGeo, idOfFacesOpposite);
1820
1821
2/2
✓ Branch 1 taken 1220904 times.
✓ Branch 2 taken 406968 times.
 1627872 for(std::size_t jface=0; jface<idOfFacesOpposite.size(); ++jface) {
1822
2/2
✓ Branch 2 taken 406968 times.
✓ Branch 3 taken 813936 times.
 1220904 if(idOfFacesCurrent[iface] == idOfFacesOpposite[jface]) {
1823 406968 localIdFaceOpposite = jface;
1824 }
1825 }
1826
1827 // get information on the opposite element
1828 // same as the function "setElemPoint" but here ielOppositeGlobalGeo is global
1829 // we assume that the supportDofMesh is the same for all unknown (like in setElemPoint)
1830
1/2
✓ Branch 2 taken 406968 times.
✗ Branch 3 not taken.
 406968 m_supportDofUnknown[0].getIdElementSupport(ielOppositeGlobalGeo, ielOppositeGlobalSup);
1831
1/2
✓ Branch 3 taken 406968 times.
✗ Branch 4 not taken.
 406968 m_mesh[m_currentMesh]->getOneElement(ielOppositeGlobalGeo, oppositeElemIdPoint);
1832
2/2
✓ Branch 0 taken 1220904 times.
✓ Branch 1 taken 406968 times.
 1627872 for (felInt iPoint=0; iPoint<numPointPerElt; ++iPoint)
1833 1220904 oppositeElemPoint[iPoint] = &(m_mesh[m_currentMesh]->listPoints()[oppositeElemIdPoint[iPoint]]);
1834
1835 // compare the number of support element
1836
2/2
✓ Branch 2 taken 402480 times.
✓ Branch 3 taken 4488 times.
 406968 if(ielCurrentGlobalSup.size() == ielOppositeGlobalSup.size()) {
1837 // both are duplicated or both are not. Take the element on the same side
1838 402480 oppositeSupport = ielOppositeGlobalSup[ielSup];
1839 402480 addOppositeSupportElt = true;
1840 } else {
1841 // One is duplicated, the other is not
1842 4488 addOppositeSupportElt = false;
1843
6/6
✓ Branch 1 taken 5984 times.
✓ Branch 2 taken 3712 times.
✓ Branch 3 taken 5208 times.
✓ Branch 4 taken 776 times.
✓ Branch 5 taken 5208 times.
✓ Branch 6 taken 4488 times.
 9696 for(std::size_t jelSup=0; jelSup < ielOppositeGlobalSup.size() && !addOppositeSupportElt; ++jelSup) {
1844 5208 oppositeSupport = ielOppositeGlobalSup[jelSup];
1845
1846 // check if this opposite support element is on the same side as the current support element
1847 5208 std::size_t sizeCurrent = m_supportDofUnknown[0].getNumSupportDof(currentSupport);
1848 5208 std::size_t sizeOpposite = m_supportDofUnknown[0].getNumSupportDof(oppositeSupport);
1849
1850 5208 felInt countSupportDofInCommon = 0;
1851
2/2
✓ Branch 0 taken 15624 times.
✓ Branch 1 taken 5208 times.
 20832 for(std::size_t iSup1=0; iSup1<sizeCurrent; ++iSup1) {
1852
2/2
✓ Branch 0 taken 46872 times.
✓ Branch 1 taken 15624 times.
 62496 for(std::size_t iSup2=0; iSup2<sizeOpposite; ++iSup2) {
1853
2/2
✓ Branch 4 taken 6224 times.
✓ Branch 5 taken 40648 times.
 46872 if(iSupportDof[iEle[oppositeSupport] + iSup2] == iSupportDof[iEle[currentSupport] + iSup1]) {
1854 6224 ++countSupportDofInCommon;
1855 }
1856 }
1857 }
1858
1859
2/2
✓ Branch 0 taken 2992 times.
✓ Branch 1 taken 2216 times.
 5208 if(countSupportDofInCommon > 1)
1860 2992 addOppositeSupportElt = true;
1861 }
1862 }
1863
1864 // if there really is an opposite support element
1865 // Even for an inner edge, it is possible that the current support element has no opposite
1866 // support element (if it is duplicated and the opposite element is not).
1867
2/2
✓ Branch 0 taken 405472 times.
✓ Branch 1 taken 1496 times.
 406968 if(addOppositeSupportElt) {
1868 // update the opposite finite elements
1869
2/4
✓ Branch 1 taken 405472 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 405472 times.
✗ Branch 5 not taken.
 405472 oppositeEltFEWithBdVel[eltTypeOpp]->updateFirstDeriv(0, oppositeElemPoint);
1870
2/4
✓ Branch 1 taken 405472 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 405472 times.
✗ Branch 5 not taken.
 405472 oppositeEltFEWithBdVel[eltTypeOpp]->updateBdMeasNormal();
1871
2/4
✓ Branch 1 taken 405472 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 405472 times.
✗ Branch 5 not taken.
 405472 oppositeEltFEWithBdPre[eltTypeOpp]->updateFirstDeriv(0, oppositeElemPoint);
1872
2/4
✓ Branch 1 taken 405472 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 405472 times.
✗ Branch 5 not taken.
 405472 oppositeEltFEWithBdPre[eltTypeOpp]->updateBdMeasNormal();
1873
1874
1875 // compute coefficient
1876
2/4
✓ Branch 2 taken 405472 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 405472 times.
✗ Branch 6 not taken.
 405472 elemFieldAdvFace.initialize(DOF_FIELD, m_feVelWithBd->bdEle(iface), dimension());
1877
2/4
✓ Branch 2 taken 405472 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 405472 times.
✗ Branch 6 not taken.
 405472 elemFieldAdvFaceTmp.initialize(DOF_FIELD, m_feVelWithBd->bdEle(iface), dimension());
1878
2/4
✓ Branch 1 taken 405472 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 405472 times.
 405472 if(FelisceParam::instance().useProjectionForNitscheXFEM) {
1879 ✗ for(felInt iextrap=0; iextrap<FelisceParam::instance().orderBdfNS; ++iextrap) {
1880 ✗ elemFieldAdvFaceTmp.setValue(m_seqVelExtrapol[iextrap], *m_feVelWithBd, iface, currentSupport, m_iVelocity, m_ao, dof());
1881 ✗ elemFieldAdvFace.val += elemFieldAdvFaceTmp.val;
1882 }
1883 } else {
1884
3/4
✓ Branch 1 taken 810944 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 405472 times.
✓ Branch 4 taken 405472 times.
 810944 for(felInt iextrap=0; iextrap<FelisceParam::instance().orderBdfNS; ++iextrap) {
1885
1/2
✓ Branch 5 taken 405472 times.
✗ Branch 6 not taken.
 405472 elemFieldAdvFaceTmp.setValue(m_seqVelExtrapol[iextrap], *m_feVelWithBd, iface, ielOld[iextrap], m_iVelocity, m_aoOld[iextrap], m_dofOld[iextrap]);
1886
1/2
✓ Branch 1 taken 405472 times.
✗ Branch 2 not taken.
 405472 elemFieldAdvFace.val += elemFieldAdvFaceTmp.val;
1887 }
1888 }
1889
1890
4/8
✓ Branch 1 taken 405472 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 405472 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 405472 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 405472 times.
✗ Branch 12 not taken.
 810944 UBlasVector normalVel = prod(trans(elemFieldAdvFace.val), m_feVelWithBd->bdEle(iface).normal[0]);
1891
1/2
✓ Branch 1 taken 405472 times.
✗ Branch 2 not taken.
 405472 double normVelBd = norm_inf(normalVel);
1892 405472 double gammaU = coeffVel * normVelBd;
1893
1894 // We can now compute all the integrals.
1895 // There is a minus sign for gammaP because there is one for q div(u).
1896 // current element for phi_i and phi_j
1897
1/2
✓ Branch 3 taken 405472 times.
✗ Branch 4 not taken.
 405472 m_elementMat[0]->zero();
1898
1/2
✓ Branch 4 taken 405472 times.
✗ Branch 5 not taken.
 405472 m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n(gammaU, *m_feVelWithBd, *m_feVelWithBd, iface, iface, 0, 0, dimension());
1899
1/2
✓ Branch 5 taken 405472 times.
✗ Branch 6 not taken.
 405472 m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n(-gammaP, *m_fePreWithBd, *m_fePreWithBd, iface, iface, dimension(), dimension(), 1);
1900 // m_elementMat[0]->grad_phi_j_grad_psi_i(-gammaP, *m_fePreWithBd, *m_fePreWithBd, iface, iface, dimension(), dimension(), 1);
1901
1/2
✓ Branch 1 taken 405472 times.
✗ Branch 2 not taken.
 405472 setValueMatrixRHS(currentSupport, currentSupport, flag);
1902
1903 // current element for phi_i and opposite element for phi_j
1904
1/2
✓ Branch 3 taken 405472 times.
✗ Branch 4 not taken.
 405472 m_elementMat[0]->zero();
1905
2/4
✓ Branch 4 taken 405472 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 405472 times.
✗ Branch 8 not taken.
 405472 m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n(gammaU, *oppositeEltFEWithBdVel[eltTypeOpp], *m_feVelWithBd, localIdFaceOpposite, iface, 0, 0, dimension());
1906
2/4
✓ Branch 5 taken 405472 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 405472 times.
✗ Branch 9 not taken.
 405472 m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n(-gammaP, *oppositeEltFEWithBdPre[eltTypeOpp], *m_fePreWithBd, localIdFaceOpposite, iface, dimension(), dimension(), 1);
1907 // m_elementMat[0]->grad_phi_j_grad_psi_i(gammaP, oppositeEltFEWithBdPre, *m_fePreWithBd, localIdFaceOpposite, iface, dimension(), dimension(), 1);
1908
1/2
✓ Branch 1 taken 405472 times.
✗ Branch 2 not taken.
 405472 setValueMatrixRHS(currentSupport, oppositeSupport, flag);
1909
1910 // current element for phi_i and opposite element for phi_j
1911
2/4
✓ Branch 4 taken 405472 times.
✗ Branch 5 not taken.
✓ Branch 10 taken 405472 times.
✗ Branch 11 not taken.
 405472 m_elementMat[0]->mat() = trans(m_elementMat[0]->mat());
1912
1/2
✓ Branch 1 taken 405472 times.
✗ Branch 2 not taken.
 405472 setValueMatrixRHS(oppositeSupport, currentSupport, flag);
1913
1914 // current element for phi_i and phi_j
1915
1/2
✓ Branch 3 taken 405472 times.
✗ Branch 4 not taken.
 405472 m_elementMat[0]->zero();
1916
3/6
✓ Branch 4 taken 405472 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 405472 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 405472 times.
✗ Branch 11 not taken.
 405472 m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n(gammaU, *oppositeEltFEWithBdVel[eltTypeOpp], *oppositeEltFEWithBdVel[eltTypeOpp], localIdFaceOpposite, localIdFaceOpposite, 0, 0, dimension());
1917
3/6
✓ Branch 5 taken 405472 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 405472 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 405472 times.
✗ Branch 12 not taken.
 405472 m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n(-gammaP, *oppositeEltFEWithBdPre[eltTypeOpp], *oppositeEltFEWithBdPre[eltTypeOpp], localIdFaceOpposite, localIdFaceOpposite, dimension(), dimension(), 1);
1918 // m_elementMat[0]->grad_phi_j_grad_psi_i(-gammaP, oppositeEltFEWithBdPre, oppositeEltFEWithBdPre, localIdFaceOpposite, localIdFaceOpposite, dimension(), dimension(), 1);
1919
1/2
✓ Branch 1 taken 405472 times.
✗ Branch 2 not taken.
 405472 setValueMatrixRHS(oppositeSupport, oppositeSupport, flag);
1920 405472 }
1921 }
1922 }
1923 138296 }
1924 136880 ++countElement[eltType];
1925 136880 ++ielCurrentLocalGeo;
1926 }
1927 }
1928 }
1929
1930 // desallocate array use for assemble with petsc.
1931
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 desallocateArrayForAssembleMatrixRHS(flag);
1932
1933 // desallocate opposite finite elements
1934
2/2
✓ Branch 1 taken 40 times.
✓ Branch 2 taken 40 times.
 80 for (std::size_t i=0; i<bagElementTypeDomain.size(); ++i) {
1935 40 eltType = bagElementTypeDomain[i];
1936
2/4
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 40 times.
✗ Branch 4 not taken.
 40 delete oppositeEltFEWithBdVel[eltType];
1937
2/4
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 40 times.
✗ Branch 4 not taken.
 40 delete oppositeEltFEWithBdPre[eltType];
1938 }
1939 40 }
1940
1941
1942
1943
1944 ✗ void LinearProblemFSINitscheXFEMFluid::computeElementArrayBoundaryCondition(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint, felInt& ielSupportDof, FlagMatrixRHS flagMatrixRHS) {
1945 IGNORE_UNUSED_FLAG_MATRIX_RHS;
1946 IGNORE_UNUSED_ELEM_ID_POINT;
1947
1948 // call to the function in linear problem that updates the curvilinear finite element
1949 //LinearProblem::computeElementArrayBoundaryCondition(elemPoint, elemIdPoint, ielSupportDof, flagMatrixRHS);
1950 ✗ if(m_curvFeVel->hasOriginalQuadPoint()) {
1951 ✗ m_curvFeVel->updateMeasNormalContra(0, elemPoint);
1952 ✗ m_curvFePress->updateMeasNormalContra(0, elemPoint);
1953 } else {
1954 ✗ m_curvFeVel->updateBasisAndNormalContra(0, elemPoint);
1955 ✗ m_curvFePress->updateBasisAndNormalContra(0, elemPoint);
1956 }
1957
1958 // we change the curvilinear finite element here
1959 // the integrals will automatically be computed on the physical part of the boundary
1960
1961 // get the global id of the geometric element
1962 ✗ felInt iel = -1;
1963 ✗ bool isfound = false;
1964 ✗ std::vector<felInt> tmpSupport;
1965 ✗ felInt start = ielSupportDof > m_mesh[m_currentMesh]->getNumElement() - 1 ? m_mesh[m_currentMesh]->getNumElement() - 1 : ielSupportDof;
1966 ✗ for(felInt i=start; i>0 && !isfound; --i) {
1967 // we start at ielSupportDof (should be close to the geometric id)
1968 ✗ m_supportDofUnknown[0].getIdElementSupport(i, tmpSupport);
1969 ✗ for(std::size_t j=0; j<tmpSupport.size(); ++j) {
1970 ✗ if(tmpSupport[j] == ielSupportDof) {
1971 ✗ iel = i;
1972 ✗ isfound = true;
1973 }
1974 }
1975 }
1976
1977 ✗ if(!isfound) {
1978 ✗ FEL_ERROR("geometric element not found, check the search function");
1979 }
1980
1981 // check if this element is intersected
1982 ✗ std::vector< std::vector< Point > > bndElem;
1983 ✗ m_duplicateSupportElements->getSubBndElt(iel, bndElem);
1984
1985 ✗ if(!bndElem.empty()) {
1986 // get the "side" of the current support element
1987 sideOfInterface side;
1988 felInt ielSupport;
1989 ✗ m_supportDofUnknown[0].getIdElementSupport(iel, ielSupport);
1990 ✗ side = ielSupport == ielSupportDof ? LEFT : RIGHT;
1991
1992 // get the side of the sub elements
1993 ✗ std::vector<sideOfInterface> bndElemSide;
1994 ✗ m_duplicateSupportElements->getSubBndEltSide(iel, bndElemSide);
1995
1996 // integration over the left or right sub element only
1997 ✗ for(std::size_t iSubElt=0; iSubElt<bndElemSide.size(); ++iSubElt) {
1998 ✗ if(bndElemSide[iSubElt] == side) {
1999 // get the points defining the sub element
2000 ✗ for(std::size_t iPt=0; iPt < m_numVerPerSolidElt; ++iPt) {
2001 ✗ *m_subItfPts[iPt] = bndElem[iSubElt][iPt];
2002 }
2003
2004 // update measure
2005 ✗ m_curvFeVel->updateSubElementMeasNormal(0, elemPoint, m_subItfPts);
2006 ✗ m_curvFePress->updateSubElementMeasNormal(0, elemPoint, m_subItfPts);
2007 }
2008 }
2009 }
2010 }
2011
2012
2013
2014 320 void LinearProblemFSINitscheXFEMFluid::computeNitscheSigma_u_v(double velCoeff, double preCoeff) {
2015 // \int_K (velCoeff \grad u vec - preCoeff p vec) \cdot v
2016
1/2
✓ Branch 0 taken 320 times.
✗ Branch 1 not taken.
 320 if(m_useSymmetricStress)
2017 320 m_elementMat[0]->phi_i_eps_phi_j_dot_vec(2.*velCoeff, m_elemFieldNormal, *m_feVel, 0, 0, dimension());
2018 else
2019 ✗ m_elementMat[0]->phi_i_grad_phi_j_dot_vec(velCoeff, m_elemFieldNormal, *m_feVel, 0, 0, dimension());
2020
2021 // same finite element for the pressure and the velocity
2022
2/2
✓ Branch 1 taken 640 times.
✓ Branch 2 taken 320 times.
 960 for(felInt idim=0; idim<dimension(); ++idim)
2023 640 m_elementMat[0]->phi_i_phi_j(preCoeff*m_elemFieldNormal.val(idim,0), *m_feVel, idim, dimension(), 1);
2024 320 }
2025
2026
2027 4432 void LinearProblemFSINitscheXFEMFluid::computeNitscheSigma_v_u(double velCoeff, double preCoeff) {
2028 // \int_K (velCoeff \grad v vec - preCoeff q vec) \cdot u
2029
1/2
✓ Branch 0 taken 4432 times.
✗ Branch 1 not taken.
 4432 if(m_useSymmetricStress)
2030 4432 m_elementMat[0]->eps_phi_i_dot_vec_phi_j(2.*velCoeff, m_elemFieldNormal, *m_feVel, 0, 0, dimension());
2031 else
2032 ✗ m_elementMat[0]->grad_phi_i_dot_vec_phi_j(velCoeff, m_elemFieldNormal, *m_feVel, 0, 0, dimension());
2033
2034 // same finite element for the pressure and velocity
2035
2/2
✓ Branch 1 taken 8864 times.
✓ Branch 2 taken 4432 times.
 13296 for(felInt idim=0; idim<dimension(); ++idim)
2036 8864 m_elementMat[0]->phi_i_phi_j(preCoeff*m_elemFieldNormal.val(idim,0), *m_feVel, dimension(), idim, 1);
2037 4432 }
2038
2039
2040 ✗ void LinearProblemFSINitscheXFEMFluid::computeNitscheSigma_u_Sigma_v(double coeff) {
2041 // same finite element for pressure and velocity
2042 // eps(u)n . eps(v)n, eps(u)n . qn, -pn . eps(v)n
2043 ✗ if(m_useSymmetricStress) {
2044 ✗ m_elementMat[0]->eps_phi_j_vec_dot_eps_phi_i_vec(4.*m_viscosity*m_viscosity*coeff, m_elemFieldNormal, *m_feVel, 0, 0, dimension());
2045
2046 ✗ for(felInt idim=0; idim<dimension(); ++idim) {
2047 ✗ m_elementMat[0]->phi_i_eps_phi_j_dot_vec(2.*m_viscosity*coeff*m_elemFieldNormal.val(idim,0), m_elemFieldNormal, *m_feVel, dimension(), idim, 1);
2048
2049 ✗ m_elementMat[0]->eps_phi_i_dot_vec_phi_j(-2.*m_viscosity*coeff*m_elemFieldNormal.val(idim,0), m_elemFieldNormal, *m_feVel, idim, dimension(), 1);
2050 }
2051 } else {
2052 ✗ m_elementMat[0]->grad_phi_j_vec_dot_grad_phi_i_vec(m_viscosity*m_viscosity*coeff, m_elemFieldNormal, *m_feVel, 0, 0, dimension());
2053
2054 ✗ for(felInt idim=0; idim<dimension(); ++idim) {
2055 ✗ m_elementMat[0]->phi_i_grad_phi_j_dot_vec(m_viscosity*coeff*m_elemFieldNormal.val(idim,0), m_elemFieldNormal, *m_feVel, dimension(), idim, 1);
2056
2057 ✗ m_elementMat[0]->grad_phi_i_dot_vec_phi_j(-m_viscosity*coeff*m_elemFieldNormal.val(idim,0), m_elemFieldNormal, *m_feVel, idim, dimension(), 1);
2058 }
2059 }
2060
2061 // -pn . qn = -pq
2062 ✗ m_elementMat[0]->phi_i_phi_j(-coeff, *m_fePres, dimension(), dimension(), 1);
2063 }
2064
2065
2066
2067
2068 ✗ void LinearProblemFSINitscheXFEMFluid::computeNitscheSigma_uold_Sigma_v(double coeff) {
2069 // same finite element for pressure and velocity
2070 // eps(uold)n . eps(v)n, eps(uold)n . qn, -poldn . eps(v)n
2071 ✗ if(m_useSymmetricStress) {
2072 ✗ m_elementVector[0]->eps_vec2_vec1_dot_eps_phi_i_vec1(4.*m_viscosity*m_viscosity*coeff, m_elemFieldNormal, m_seqVelDofPts[0], *m_feVel, 0, dimension());
2073
2074 ✗ for(felInt idim=0; idim<dimension(); ++idim) {
2075 ✗ m_elementVector[0]->eps_vec2_dot_vec1_dot_phi_i(2.*m_viscosity*coeff*m_elemFieldNormal.val(idim,0), m_elemFieldNormal, m_seqVelDofPts[0], *m_feVel, dimension(), 1);
2076
2077 ✗ m_elementVector[0]->eps_phi_i_dot_vec1_dot_vec2(-2.*m_viscosity*coeff*m_elemFieldNormal.val(idim,0), m_elemFieldNormal, m_seqVelDofPts[0], *m_feVel, idim, 1);
2078 }
2079 } else {
2080 ✗ m_elementVector[0]->grad_vec2_vec1_dot_grad_phi_i_vec1(m_viscosity*m_viscosity*coeff, m_elemFieldNormal, m_seqVelDofPts[0], *m_feVel, 0, dimension());
2081
2082 ✗ for(felInt idim=0; idim<dimension(); ++idim) {
2083 ✗ m_elementVector[0]->grad_vec2_dot_vec1_dot_phi_i(m_viscosity*coeff*m_elemFieldNormal.val(idim,0), m_elemFieldNormal, m_seqVelDofPts[0], *m_feVel, dimension(), 1);
2084
2085 ✗ m_elementVector[0]->grad_phi_i_dot_vec1_dot_vec2(-m_viscosity*coeff*m_elemFieldNormal.val(idim,0), m_elemFieldNormal, m_seqVelDofPts[0], *m_feVel, idim, 1);
2086 }
2087 }
2088
2089 // -poldn . qn = -pold q
2090 ✗ m_elementVector[0]->source(-coeff, *m_fePres, m_seqPreDofPts[0], dimension(), 1);
2091
2092 }
2093
2094
2095
2096 4112 void LinearProblemFSINitscheXFEMFluid::computeNitscheDpoint_Sigma_v(double velCoeff, double preCoeff) {
2097 // \int_K ddot \cdot (velCoeff \grad v vec - preCoeff q vec)
2098
1/2
✓ Branch 0 taken 4112 times.
✗ Branch 1 not taken.
 4112 if(m_useSymmetricStress)
2099 4112 m_elementVector[0]->eps_phi_i_dot_vec1_dot_vec2(2.*velCoeff, m_elemFieldNormal, m_strucVelQuadPts, *m_feVel, 0, dimension());
2100 else
2101 ✗ m_elementVector[0]->grad_phi_i_dot_vec1_dot_vec2(velCoeff, m_elemFieldNormal, m_strucVelQuadPts, *m_feVel, 0, dimension());
2102
2103
2/2
✓ Branch 1 taken 8224 times.
✓ Branch 2 taken 4112 times.
 12336 for(felInt idim=0; idim<dimension(); ++idim) {
2104
2/2
✓ Branch 1 taken 24672 times.
✓ Branch 2 taken 8224 times.
 32896 for(felInt ig=0; ig<m_feVel->numQuadraturePoint(); ++ig) {
2105 24672 m_ddotComp.val(0, ig) = m_strucVelQuadPts.val(idim, ig);
2106 }
2107 8224 m_elementVector[0]->source(preCoeff*m_elemFieldNormal.val(idim,0), *m_feVel, m_ddotComp, dimension(), 1);
2108 }
2109 4112 }
2110
2111 ✗ void LinearProblemFSINitscheXFEMFluid::computeStenbergSigma_u_Sigma_v(double velCoeff, double preCoeff) {
2112 // \int_K (velCoeff \grad u vec - preCoeff p vec) \cdot q vec
2113 // assuming same finite element for the pressure and the velocity
2114 ✗ std::vector<double> vec(dimension());
2115 ✗ for(felInt idim=0; idim<dimension(); ++idim){
2116 ✗ for(felInt i=0; i<dimension(); ++i){
2117 ✗ vec[i] = 2 * m_elemFieldNormal.val(i,0) * m_elemFieldNormal.val(idim,0);
2118 }
2119 ✗ m_elemFieldAux.setValue(vec);
2120 ✗ if(m_useSymmetricStress)
2121 ✗ m_elementMat[0]->phi_i_eps_phi_j_dot_vec(2. * velCoeff , m_elemFieldAux, *m_feVel, dimension(), idim, 1);
2122 else
2123 ✗ m_elementMat[0]->phi_i_grad_phi_j_dot_vec( velCoeff , m_elemFieldAux, *m_feVel, dimension(), idim, 1);
2124 }
2125 ✗ m_elementMat[0]->phi_i_phi_j(preCoeff, *m_feVel, dimension() , dimension(), 1);
2126 }
2127
2128
2129
2130
2131 2 void LinearProblemFSINitscheXFEMFluid::setStrucFiniteElement(CurvilinearFiniteElement* strucFE) {
2132 2 m_feStruc = strucFE;
2133 2 }
2134
2135 2 void LinearProblemFSINitscheXFEMFluid::setStrucLinPrb(LinearProblem* pLinPrb) {
2136 2 m_pLinStruc = pLinPrb;
2137 2 }
2138
2139
2140
2141
2142
2143 4112 void LinearProblemFSINitscheXFEMFluid::updateStructureVel(felInt strucId, PetscVector& strucVel) {
2144 felInt idSupportEltGlobal;
2145 felInt idDof;
2146
1/2
✓ Branch 2 taken 4112 times.
✗ Branch 3 not taken.
 4112 felInt idUnknown = m_pLinStruc->listUnknown().getUnknownIdList(displacement);
2147 4112 felInt idVar = m_pLinStruc->listUnknown().idVariable(idUnknown);
2148
2149 // get the global support element
2150 // there should be only one for the structure
2151
1/2
✓ Branch 2 taken 4112 times.
✗ Branch 3 not taken.
 4112 m_pLinStruc->supportDofUnknown(idUnknown).getIdElementSupport(strucId, idSupportEltGlobal);
2152
2153 // loop over all the support dof of this support element
2154
2/2
✓ Branch 2 taken 8224 times.
✓ Branch 3 taken 4112 times.
 12336 for (felInt iSupport = 0; iSupport < m_pLinStruc->supportDofUnknown(idUnknown).getNumSupportDof(idSupportEltGlobal); iSupport++) {
2155 // loop over the components
2156 // for(std::size_t iComp=0; iComp<m_pLinStruc->meshLocal()->numCoor(); ++iComp){
2157
2/2
✓ Branch 3 taken 12336 times.
✓ Branch 4 taken 8224 times.
 20560 for(std::size_t iComp=0; iComp<m_pLinStruc->listVariable()[idVar].numComponent(); ++iComp) {
2158 // local to global id of the dof
2159
1/2
✓ Branch 2 taken 12336 times.
✗ Branch 3 not taken.
 12336 m_pLinStruc->dof().loc2glob(idSupportEltGlobal, iSupport, idVar, iComp, idDof);
2160
2161 // petsc ordering
2162
1/2
✓ Branch 2 taken 12336 times.
✗ Branch 3 not taken.
 12336 AOApplicationToPetsc(m_pLinStruc->ao(), 1, &idDof);
2163
2164 // get values
2165
2/4
✓ Branch 1 taken 12336 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 12336 times.
✗ Branch 5 not taken.
 12336 strucVel.getValues(1, &idDof, &m_strucVelDofPts.val(iComp, iSupport));
2166 }
2167 }
2168 4112 }
2169
2170 4112 void LinearProblemFSINitscheXFEMFluid::updateSubStructureVel(const std::vector<Point*>& ptElem, const std::vector<Point*>& ptSubElem) {
2171
1/4
✗ Branch 1 not taken.
✓ Branch 2 taken 4112 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 4112 FEL_ASSERT(ptElem.size() == 2);
2172
2173 // compute the value of the displacement at the new integration points
2174 // get the coordinate of the integration point in the sub element
2175
1/2
✓ Branch 1 taken 4112 times.
✗ Branch 2 not taken.
 4112 m_feStruc->update(0, ptSubElem);
2176
1/2
✓ Branch 1 taken 4112 times.
✗ Branch 2 not taken.
 4112 m_feStruc->computeCurrentQuadraturePoint();
2177
2178 // compute the coordinate of these integration points in the reference element with the mapping of the
2179 // structure element (not the structure sub element)
2180
1/2
✓ Branch 3 taken 4112 times.
✗ Branch 4 not taken.
 4112 std::vector<Point> intPoints(m_feStruc->numQuadraturePoint());
2181 double dx, dy;
2182 4112 dx = std::abs(ptElem[1]->x() - ptElem[0]->x());
2183 4112 dy = std::abs(ptElem[1]->y() - ptElem[0]->y());
2184
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4112 times.
 4112 if(dx >= dy) {
2185 ✗ for(int ig=0; ig<m_feStruc->numQuadraturePoint(); ++ig) {
2186 ✗ intPoints[ig][0] = 2.*(m_feStruc->currentQuadPoint[ig].x() - ptElem[0]->x())/(ptElem[1]->x() - ptElem[0]->x()) - 1.;
2187 ✗ intPoints[ig][1] = 0.;
2188 }
2189 } else {
2190
2/2
✓ Branch 1 taken 8224 times.
✓ Branch 2 taken 4112 times.
 12336 for(int ig=0; ig<m_feStruc->numQuadraturePoint(); ++ig) {
2191 8224 intPoints[ig][0]= 2.*(m_feStruc->currentQuadPoint[ig].y() - ptElem[0]->y())/(ptElem[1]->y() - ptElem[0]->y()) - 1.;
2192 8224 intPoints[ig][1] = 0;
2193 }
2194 }
2195
2196 // update the structure finite element with the element (really needed ? we only want to access the basis
2197 // function directly to evaluate them at the computed integration points).
2198
1/2
✓ Branch 1 taken 4112 times.
✗ Branch 2 not taken.
 4112 m_feStruc->update(0, ptElem);
2199
2200
2201 // Now we can evaluate the basis function of the reference element at the integration points of
2202 // the sub element
2203
2/2
✓ Branch 1 taken 2056 times.
✓ Branch 2 taken 2056 times.
 4112 if(m_strucVelDofPts.numComp() == 1) {
2204 // The solid is a scalar function. The displacement is in the normal direction.
2205
2/2
✓ Branch 1 taken 4112 times.
✓ Branch 2 taken 2056 times.
 6168 for(int icomp=0; icomp<m_feVel->numCoor(); ++icomp) {
2206
2/2
✓ Branch 1 taken 8224 times.
✓ Branch 2 taken 4112 times.
 12336 for(int ig=0; ig<m_feStruc->numQuadraturePoint(); ++ig) {
2207
1/2
✓ Branch 1 taken 8224 times.
✗ Branch 2 not taken.
 8224 m_strucVelQuadPts.val(icomp, ig) = 0;
2208
2/2
✓ Branch 1 taken 16448 times.
✓ Branch 2 taken 8224 times.
 24672 for(int idof=0; idof<m_feStruc->numDof(); ++idof) {
2209
4/8
✓ Branch 1 taken 16448 times.
✗ Branch 2 not taken.
✓ Branch 7 taken 16448 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 16448 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 16448 times.
✗ Branch 14 not taken.
 16448 m_strucVelQuadPts.val(icomp, ig) += m_strucVelDofPts.val(0, idof) * m_feStruc->refEle().basisFunction().phi(idof, intPoints[ig]) * m_elemFieldNormalInit.val(icomp, 0);
2210 }
2211 }
2212 }
2213 } else {
2214 // The solid is a std::vector.
2215
2/2
✓ Branch 1 taken 4112 times.
✓ Branch 2 taken 2056 times.
 6168 for(int icomp=0; icomp<m_feVel->numCoor(); ++icomp) {
2216
2/2
✓ Branch 1 taken 8224 times.
✓ Branch 2 taken 4112 times.
 12336 for(int ig=0; ig<m_feStruc->numQuadraturePoint(); ++ig) {
2217
1/2
✓ Branch 1 taken 8224 times.
✗ Branch 2 not taken.
 8224 m_strucVelQuadPts.val(icomp, ig) = 0;
2218
2/2
✓ Branch 1 taken 16448 times.
✓ Branch 2 taken 8224 times.
 24672 for(int idof=0; idof<m_feStruc->numDof(); ++idof) {
2219
3/6
✓ Branch 1 taken 16448 times.
✗ Branch 2 not taken.
✓ Branch 7 taken 16448 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 16448 times.
✗ Branch 11 not taken.
 16448 m_strucVelQuadPts.val(icomp, ig) += m_strucVelDofPts.val(icomp, idof) * m_feStruc->refEle().basisFunction().phi(idof, intPoints[ig]);
2220 }
2221 }
2222 }
2223 }
2224 4112 }
2225
2226
2227 141072 void LinearProblemFSINitscheXFEMFluid::computeElementMatrixRHS() {
2228 // mass coefficient
2229 141072 double coef = m_density/m_fstransient->timeStep;
2230
2231 //================================
2232 // matrix
2233 //================================
2234
2235 // Laplacian operator
2236
1/2
✓ Branch 0 taken 141072 times.
✗ Branch 1 not taken.
 141072 if (m_useSymmetricStress) {
2237 // \mu \int \epsilon u^{n+1} : \epsilon v
2238 141072 m_elementMat[0]->eps_phi_i_eps_phi_j(2*m_viscosity,*m_feVel, 0, 0, dimension());
2239 } else {
2240 // \mu \int \nabla u^{n+1} : \nabla v
2241 ✗ m_elementMat[0]->grad_phi_i_grad_phi_j(m_viscosity, *m_feVel, 0, 0, dimension());
2242 }
2243
2244 // div and grad operator
2245 // \int q \nabla \cdot u^{n+1} and \int p \nabla \cdot v
2246 141072 m_elementMat[0]->psi_j_div_phi_i_and_psi_i_div_phi_j (*m_feVel, *m_fePres, 0, m_velocity->numComponent(), -1, -1);
2247
2248 // mass matrix from the time scheme
2249 // \int u^{n+1} \cdot v
2250 141072 m_elementMat[0]->phi_i_phi_j(m_bdf->coeffDeriv0()*coef, *m_feVel, 0, 0, dimension());
2251 141072 }
2252
2253
2254 141168 void LinearProblemFSINitscheXFEMFluid::computeElementMatrixRHSPartDependingOnOldTimeStep(felInt iel, felInt ielOld, felInt idTimeStep) {
2255 // --------------- //
2256 // convective term //
2257 // --------------- //
2258
1/3
✗ Branch 1 not taken.
✓ Branch 2 taken 141168 times.
✗ Branch 3 not taken.
 141168 switch(FelisceParam::instance().NSequationFlag) {
2259 ✗ case 0:
2260 // Stoke's model, no convective term
2261 ✗ break;
2262
2263 141168 case 1:
2264 // classic Navier Stokes convective term : \rho (u \cdot \grad) u
2265 // Natural condition: \sigma(u,p) \cdot n = 0
2266 // discretization : \rho (u^{n} \cdot \grad) u^{n+1}
2267 // DEFAULT METHOD
2268
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 141168 times.
 141168 if(FelisceParam::instance().useProjectionForNitscheXFEM)
2269 ✗ m_elemFieldAdv.setValue(m_seqVelExtrapol[idTimeStep], *m_feVel, iel, m_iVelocity, m_ao, dof());
2270 else
2271 141168 m_elemFieldAdv.setValue(m_seqVelExtrapol[idTimeStep], *m_feVel, ielOld, m_iVelocity, m_aoOld[idTimeStep], m_dofOld[idTimeStep]);
2272
2273 141168 m_elementMat[0]->u_grad_phi_j_phi_i(m_density, m_elemFieldAdv, *m_feVel, 0, 0, dimension());
2274
2275
2276 // Temam's operator to stabilize the convective term
2277 // Convective term in energy balance :
2278 // ... + \rho/2 \int_Gamma (u^n.n)|u^n+1|^2 - \rho/2 \int_Omega (\nabla . u^n) |u^n+1|^2 + ...
2279 // Add \rho/2 \int_Omega (\nabla . u^n) u^n+1 v to delete the second term.
2280 // it's consistent with the exact solution (divergence = 0)
2281 141168 m_elementMat[0]->div_u_phi_j_phi_i(0.5*m_density, m_elemFieldAdv, *m_feVel, 0, 0, dimension());
2282 141168 break;
2283
2284 ✗ default:
2285 ✗ FEL_ERROR("Error: Wrong or not implemented convection method, check the param NSequationFlag in the data file");
2286 }
2287
2288
2289
2290 //================================
2291 // RHS
2292 //================================
2293 // bdf
2294 // \frac{\rho}{\dt} \sum_{i=1}^k \alpha_i u_{n+1-i}
2295 // 'm_density' and not 'coef' because the time step is integrated into the bdf term
2296
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 141168 times.
 141168 if(FelisceParam::instance().useProjectionForNitscheXFEM)
2297 ✗ m_elemFieldRHSbdf.setValue(m_seqBdfRHS[idTimeStep], *m_feVel, iel, m_iVelocity, m_ao, dof());
2298 else
2299 141168 m_elemFieldRHSbdf.setValue(m_seqBdfRHS[idTimeStep], *m_feVel, ielOld, m_iVelocity, m_aoOld[idTimeStep], m_dofOld[idTimeStep]);
2300
2301
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 141168 times.
 141168 assert(!m_elementVector.empty());
2302 141168 m_elementVector[0]->source(m_density, *m_feVel, m_elemFieldRHSbdf, 0, dimension());
2303
2304
2305 // --------------------------------------------------------------------- //
2306 // SUPG Stabilization (works with a bdf 1 but wrong with a higher order) //
2307 // --------------------------------------------------------------------- //
2308
2/6
✗ Branch 2 not taken.
✓ Branch 3 taken 141168 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 141168 times.
✗ Branch 8 not taken.
 141168 if((m_velocity->finiteElementType() == m_pressure->finiteElementType()) || FelisceParam::instance().NSequationFlag == 1) {
2309
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 141168 times.
 141168 if(FelisceParam::instance().NSStabType == 0) {
2310 ✗ double stab1 = FelisceParam::instance().stabSUPG;
2311 ✗ double stab2 = FelisceParam::instance().stabdiv;
2312 ✗ int type = FelisceParam::instance().typeSUPG;
2313 double Re_elem, tau;
2314
2315 // compute advection velocity
2316 ✗ if(FelisceParam::instance().useProjectionForNitscheXFEM)
2317 ✗ m_elemFieldAdv.setValue(m_seqVelExtrapol[idTimeStep], *m_feVel, iel, m_iVelocity, m_ao, dof());
2318 else
2319 ✗ m_elemFieldAdv.setValue(m_seqVelExtrapol[idTimeStep], *m_feVel, ielOld, m_iVelocity, m_aoOld[idTimeStep], m_dofOld[idTimeStep]);
2320
2321 // Beware: if elemFieldRHS is modified in userNS, it wont be considered here...
2322 // (usually elemFieldRHS = 0 (no body force))
2323 ✗ m_elementMat[0]->stab_supg(m_fstransient->timeStep,
2324 ✗ stab1, stab2, m_density, m_viscosity, *m_feVel,
2325 ✗ m_elemFieldAdv, m_elemFieldRHS, *m_elementVector[0],
2326 Re_elem, tau, type);
2327 }
2328 }
2329 141168 }
2330
2331
2332
2333
2334
2335 /*!
2336 * Gather the extrapolation of the velocity into m_seqVelExtrapol.
2337 * If m_seqVelExtrapol is null, it will be created and allocated during the gathering
2338 * according to the parallel shape of parallelVec
2339 */
2340 20 void LinearProblemFSINitscheXFEMFluid::gatherSeqVelExtrapol(std::vector<PetscVector>& parallelVec) {
2341
2/2
✓ Branch 1 taken 20 times.
✓ Branch 2 taken 20 times.
 40 for(std::size_t i=0; i<parallelVec.size(); ++i)
2342 20 gatherVector(parallelVec[i], m_seqVelExtrapol[i]);
2343
2344 20 m_isSeqVelExtrapolAllocated = true;
2345 20 }
2346
2347 ✗ PetscVector& LinearProblemFSINitscheXFEMFluid::getVelExtrapol(felInt n) {
2348 ✗ return m_seqVelExtrapol[n];
2349 }
2350
2351
2352 /*!
2353 * Gather the rhs part of the time derivative of the velocity into m_seqBdfRHS.
2354 * If m_seqBdfRHS is null, it will be created and allocated during the gathering
2355 * according to the parallel shape of bdf->vector()
2356 */
2357 20 void LinearProblemFSINitscheXFEMFluid::gatherSeqBdfRHS(std::vector<PetscVector>& parallelVec) {
2358
2/2
✓ Branch 1 taken 20 times.
✓ Branch 2 taken 20 times.
 40 for(std::size_t i=0; i<parallelVec.size(); ++i)
2359 20 gatherVector(parallelVec[i], m_seqBdfRHS[i]);
2360
2361 20 m_isSeqBdfRHSAllocated = true;
2362 20 }
2363
2364
2365
2366 ✗ void LinearProblemFSINitscheXFEMFluid::addMatrixRHS() {
2367 ✗ matrix(0).axpy(1,matrix(1),SAME_NONZERO_PATTERN);
2368 }
2369
2370 /*!
2371 \brief Duplicate the support dof intersected by the structure
2372 */
2373 13 void LinearProblemFSINitscheXFEMFluid::initSupportDofDerivedProblem() {
2374 ///////////////////////////////
2375 // duplicate the support dof //
2376 ///////////////////////////////
2377 13 m_duplicateSupportElements->duplicateIntersectedSupportElements(m_supportDofUnknown);
2378 13 }
2379
2380
2381 /*!
2382 \brief Initialize the bdf for the solution defined on the non duplicated support dof mesh
2383 */
2384 2 void LinearProblemFSINitscheXFEMFluid::initOriginalBdf(felInt timeScheme, felInt extrapolOrder) {
2385 // create the std::vector solution with the original supportDofMesh
2386 // first, compute the number of dof in the original supportDofMesh
2387 2 felInt idVar = 0;
2388 2 felInt numSupportOfDof = 0;
2389 2 m_numOriginalDof = 0;
2390
2391
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 2 times.
 6 for (std::size_t iUnknown = 0; iUnknown < m_listUnknown.size(); ++iUnknown) {
2392 4 idVar = m_listUnknown.idVariable(iUnknown);
2393 4 numSupportOfDof = m_supportDofUnknownOriginal[iUnknown]->numSupportDof();
2394 4 m_numOriginalDof += numSupportOfDof * m_listVariable[idVar].numComponent();
2395 }
2396
2397 2 m_tmpSeqOriginalSol.create(PETSC_COMM_SELF);
2398 2 m_tmpSeqOriginalSol.setType(VECSEQ);
2399 2 m_tmpSeqOriginalSol.setSizes(PETSC_DECIDE, m_numOriginalDof);
2400
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_tmpSeqOriginalSol.set(0.);
2401
2402 // initialize the bdf storing all fluid solution with the original support dof mesh
2403 2 m_seqOriginalSolutions.defineOrder(FelisceParam::instance().orderBdfNS);
2404
2405
4/8
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
 2 if(FelisceParam::instance().orderBdfNS == 1 && timeScheme == 3 && extrapolOrder == 2) {
2406 2 m_seqOriginalSolutions.initialize(m_tmpSeqOriginalSol, m_tmpSeqOriginalSol);
2407 ✗ } else if ( ((timeScheme == 0) || (timeScheme == 2)) && extrapolOrder == 2 ) {
2408 ✗ m_seqOriginalSolutions.initialize(m_tmpSeqOriginalSol, m_tmpSeqOriginalSol);
2409 } else {
2410 ✗ m_seqOriginalSolutions.initialize(m_tmpSeqOriginalSol);
2411 }
2412 2 m_isSeqOriginalSolAllocated = true;
2413 2 }
2414
2415
2416 2 void LinearProblemFSINitscheXFEMFluid::setDuplicateSupportObject(DuplicateSupportDof* object) {
2417 2 m_duplicateSupportElements = object;
2418 2 }
2419
2420
2421 13 void LinearProblemFSINitscheXFEMFluid::copyCurrentSolToOriginalSol() {
2422
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 13 times.
 13 if(!m_isSeqOriginalSolAllocated) {
2423 ✗ FEL_ERROR("The bdf storing all the data with the original supportDofMesh should already be created at this point");
2424 }
2425
2426 // get the last original solution (will be deleted after the update anyway)
2427 // m_tmpSeqOriginalSol.copyFrom(m_seqOriginalSolutions.sol_n());
2428 // m_tmpSeqOriginalSol.set(0.);
2429 13 felInt numTotalDof = 0;
2430 13 felInt numTotalDofOriginal = 0;
2431
2432 // Copy only the value of the original support dofs
2433 // All the id of the duplicated points are at the end. The ids of the original points are thus from 0 to the
2434 // number of support dof in the original supportDofMesh
2435
2/2
✓ Branch 1 taken 26 times.
✓ Branch 2 taken 13 times.
 39 for (std::size_t iUnknown = 0; iUnknown < m_listUnknown.size(); ++iUnknown) {
2436 26 felInt idVar = m_listUnknown.idVariable(iUnknown);
2437 26 felInt numOriginalSupportDof = m_supportDofUnknownOriginal[iUnknown]->numSupportDof();
2438 26 felInt numOriginalDofOfUnknown = numOriginalSupportDof * m_listVariable[idVar].numComponent();
2439
2440
1/2
✓ Branch 2 taken 26 times.
✗ Branch 3 not taken.
 26 std::vector<felInt> ids(numOriginalDofOfUnknown);
2441
1/2
✓ Branch 2 taken 26 times.
✗ Branch 3 not taken.
 26 std::vector<double> values(numOriginalDofOfUnknown, 0.);
2442
2/2
✓ Branch 1 taken 70590 times.
✓ Branch 2 taken 26 times.
 70616 for(std::size_t i=0; i<ids.size(); ++i) {
2443 70590 ids[i] = i + numTotalDof;
2444 }
2445
2446
1/2
✓ Branch 4 taken 26 times.
✗ Branch 5 not taken.
 26 sequentialSolution().getValues(numOriginalDofOfUnknown, ids.data(), values.data());
2447
2448
2/2
✓ Branch 1 taken 70590 times.
✓ Branch 2 taken 26 times.
 70616 for(std::size_t i=0; i<ids.size(); ++i) {
2449 70590 ids[i] = i + numTotalDofOriginal;
2450 }
2451
2452
1/2
✓ Branch 3 taken 26 times.
✗ Branch 4 not taken.
 26 m_tmpSeqOriginalSol.setValues(numOriginalDofOfUnknown, ids.data(), values.data(), INSERT_VALUES);
2453
2454 26 numTotalDof += supportDofUnknown(iUnknown).numSupportDof() * m_listVariable[idVar].numComponent();
2455 26 numTotalDofOriginal += numOriginalDofOfUnknown;
2456 26 }
2457
2458 13 m_tmpSeqOriginalSol.assembly();
2459
2460 // update the bdf of the original support dof mesh
2461 13 m_seqOriginalSolutions.update(m_tmpSeqOriginalSol);
2462 13 }
2463
2464
2465 11 void LinearProblemFSINitscheXFEMFluid::copyOriginalToCurrent(PetscVector& vecSeqOriginal, PetscVector& vecSeqCurrent, PetscVector& vecParCurrent) {
2466
1/2
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
 11 PetscVector tmpSeqVec;
2467
2/4
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 11 times.
 11 if(vecSeqCurrent.isNull()) {
2468 ✗ tmpSeqVec.duplicateFrom(sequentialSolution());
2469 } else {
2470
1/2
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
 11 tmpSeqVec = vecSeqCurrent;
2471 }
2472
2473 // copy the serial solution
2474 11 felInt numTotalDof = 0;
2475 11 felInt numTotalDofOriginal = 0;
2476
2477
2/2
✓ Branch 1 taken 22 times.
✓ Branch 2 taken 11 times.
 33 for (std::size_t iUnknown = 0; iUnknown < m_listUnknown.size(); ++iUnknown) {
2478 22 felInt idVar = m_listUnknown.idVariable(iUnknown);
2479 22 felInt numComp = m_listVariable[idVar].numComponent();
2480
2481 22 felInt numOriginalSupportDof = m_supportDofUnknownOriginal[iUnknown]->numSupportDof();
2482 22 felInt numOriginalDofOfUnknown = numOriginalSupportDof * numComp;
2483
2484 22 felInt numSupportDof = supportDofUnknown(iUnknown).numSupportDof();
2485 22 felInt numDofOfUnknown = numSupportDof * numComp;
2486
2487
1/2
✓ Branch 2 taken 22 times.
✗ Branch 3 not taken.
 22 std::vector<felInt> ids(numDofOfUnknown, -1);
2488
1/2
✓ Branch 2 taken 22 times.
✗ Branch 3 not taken.
 22 std::vector<double> values(numDofOfUnknown, 0.);
2489
2490 // get the id of the dofs in the original std::vector (duplicate the id of the duplicate dof in the current)
2491
2/2
✓ Branch 0 taken 59730 times.
✓ Branch 1 taken 22 times.
 59752 for(felInt i=0; i<numOriginalDofOfUnknown; ++i) {
2492 59730 ids[i] = i + numTotalDofOriginal;
2493 }
2494
2495 22 felInt count = numOriginalDofOfUnknown;
2496
1/2
✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
 22 auto& idVerMsh = m_duplicateSupportElements->getIntersectedVerMsh();
2497
2/2
✓ Branch 4 taken 818 times.
✓ Branch 5 taken 22 times.
 840 for(auto it = idVerMsh.begin(); it != idVerMsh.end(); ++it) {
2498
2/2
✓ Branch 0 taken 1227 times.
✓ Branch 1 taken 818 times.
 2045 for(felInt i=0; i<numComp; ++i) {
2499 1227 ids[count++] = numTotalDofOriginal + it->first * numComp + i;
2500 }
2501 }
2502
2503 // get the values in the serial original std::vector (we get multiple times the duplicated ones)
2504
1/2
✓ Branch 3 taken 22 times.
✗ Branch 4 not taken.
 22 vecSeqOriginal.getValues(numDofOfUnknown, ids.data(), values.data());
2505
2506 // get the id in the current std::vector
2507
2/2
✓ Branch 1 taken 60957 times.
✓ Branch 2 taken 22 times.
 60979 for(std::size_t i=0; i<ids.size(); ++i) {
2508 60957 ids[i] = i + numTotalDof;
2509 }
2510
2511
1/2
✓ Branch 3 taken 22 times.
✗ Branch 4 not taken.
 22 tmpSeqVec.setValues(numDofOfUnknown, ids.data(), values.data(), INSERT_VALUES);
2512
2513 22 numTotalDof += numDofOfUnknown;
2514 22 numTotalDofOriginal += numOriginalDofOfUnknown;
2515 22 }
2516
2517
1/2
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
 11 tmpSeqVec.assembly();
2518
2519 // scatter the serial solution to the parallel one
2520
2/4
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 11 times.
✗ Branch 4 not taken.
 11 if(vecParCurrent.isNotNull()) {
2521
1/2
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
 11 tmpSeqVec.scatterToZeroNotCreatingVector(vecParCurrent, INSERT_VALUES, SCATTER_REVERSE);
2522 }
2523
2524
2/4
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 11 times.
 11 if(vecSeqCurrent.isNull()) {
2525 ✗ tmpSeqVec.destroy();
2526 }
2527 11 }
2528
2529
2530 11 void LinearProblemFSINitscheXFEMFluid::deleteDynamicData() {
2531 // here, we delete all the matrices, vectors and mapping of the linear problem.
2532 // everything will be rebuilt at the next iteration
2533 // This is almost a copy of the destructor of LinearProblem
2534
2535 // sequential solution
2536 11 bool& isSeqSolAlloc = isSequentialSolutionAllocated();
2537
1/2
✓ Branch 0 taken 11 times.
✗ Branch 1 not taken.
 11 if(isSeqSolAlloc) {
2538 11 sequentialSolution().destroy();
2539 11 isSeqSolAlloc = false;
2540 }
2541
2542 // element matrices and matrices
2543 11 m_elementMat.clear();
2544
2545 // element vectors
2546 11 m_elementVector.clear();
2547 11 m_elementMatBD.clear();
2548 11 m_elementVectorBD.clear();
2549
2550 // solution and rhs
2551 11 bool& areSolAndRhsAlloc = areSolutionAndRHSAllocated();
2552
1/2
✓ Branch 0 taken 11 times.
✗ Branch 1 not taken.
 11 if (areSolAndRhsAlloc) {
2553
2/2
✓ Branch 1 taken 11 times.
✓ Branch 2 taken 11 times.
 22 for (unsigned int i = 0; i < numberOfVectors(); i++)
2554 11 vector(i).destroy();
2555
2556 11 solution().destroy();
2557 11 areSolAndRhsAlloc = false;
2558 }
2559
2560 // finite elements
2561 11 m_listCurrentFiniteElement.clear();
2562 11 m_listCurvilinearFiniteElement.clear();
2563
2564 // mappings
2565 11 bool& isMapElemAlloc = isMappingElemAllocated();
2566
1/2
✓ Branch 0 taken 11 times.
✗ Branch 1 not taken.
 11 if(isMapElemAlloc) {
2567 11 ISLocalToGlobalMappingDestroy(&mappingElem());
2568 11 isMapElemAlloc = false;
2569 }
2570
2571 11 bool& isMapElemSupportAlloc = isMappingElemSupportAllocated();
2572
1/2
✓ Branch 0 taken 11 times.
✗ Branch 1 not taken.
 11 if(isMapElemSupportAlloc) {
2573
2/2
✓ Branch 1 taken 22 times.
✓ Branch 2 taken 11 times.
 33 for ( std::size_t iUnknown = 0; iUnknown < m_listUnknown.size(); iUnknown++ )
2574 22 ISLocalToGlobalMappingDestroy(mappingElemSupportPerUnknown(iUnknown));
2575
2576 11 isMapElemSupportAlloc = false;
2577 }
2578
2579 11 bool& isMapNodesAlloc = isMappingNodesAllocated();
2580
1/2
✓ Branch 0 taken 11 times.
✗ Branch 1 not taken.
 11 if(isMapNodesAlloc) {
2581 11 ISLocalToGlobalMappingDestroy(&mappingNodes());
2582 11 isMapNodesAlloc = false;
2583 }
2584
2585 11 bool& isMapFelToPetscAlloc = isMappingLocalFelisceToGlobalPetscAllocated();
2586
1/2
✓ Branch 0 taken 11 times.
✗ Branch 1 not taken.
 11 if(isMapFelToPetscAlloc) {
2587
2/2
✓ Branch 1 taken 22 times.
✓ Branch 2 taken 11 times.
 33 for ( std::size_t iUnknown = 0; iUnknown < m_listUnknown.size(); iUnknown++ )
2588 22 ISLocalToGlobalMappingDestroy(&mappingIdFelisceToPetsc(iUnknown));
2589
2590 11 isMapFelToPetscAlloc = false;
2591 }
2592
2593 // AO (copy the AO, it will be destroy with the old one)
2594 11 bool& isAOAlloc = isAOAllocated();
2595
1/2
✓ Branch 0 taken 11 times.
✗ Branch 1 not taken.
 11 if(isAOAlloc) {
2596 11 AODestroy(&ao());
2597 11 isAOAlloc = false;
2598 }
2599
2600
2601 // dofs (the pattern is useless for the old dof)
2602 11 dof().numDof() = 0;
2603 11 dof().numDofPerUnknown().clear();
2604 11 dof().clearPattern();
2605
2606
2607 // repartition of the dof
2608 11 m_dofPart.clear();
2609 11 m_eltPart.clear();
2610
2611
2612 // supportDofMesh (copy the old one)
2613 11 m_supportDofUnknown.clear();
2614 11 m_supportDofUnknownLocal.clear();
2615
2616
2617 // local mesh
2618 11 m_meshLocal[m_currentMesh]->domainDim() = GeometricMeshRegion::GeoMeshUndefined;
2619 11 m_meshLocal[m_currentMesh]->deleteElements();
2620
2621 // extrapolation and time rhs of this class
2622
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
 11 if(m_isSeqVelExtrapolAllocated) {
2623 ✗ for(std::size_t i=0; i<m_seqVelExtrapol.size(); ++i) {
2624 ✗ m_seqVelExtrapol[i].destroy();
2625 ✗ m_seqVelExtrapol[i] = PetscVector::null();
2626 }
2627 ✗ m_isSeqVelExtrapolAllocated = false;
2628 }
2629
2630
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
 11 if(m_isSeqBdfRHSAllocated) {
2631 ✗ for(std::size_t i=0; i<m_seqBdfRHS.size(); ++i) {
2632 ✗ m_seqBdfRHS[i].destroy();
2633 ✗ m_seqBdfRHS[i] = PetscVector::null();
2634 }
2635 ✗ m_isSeqBdfRHSAllocated = false;
2636 }
2637
2638
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
 11 if(m_isSeqRNFluidExtrapolAllocated) {
2639 ✗ for(std::size_t i=0; i<m_seqRNFluidExtrapol.size(); ++i) {
2640 ✗ m_seqRNFluidExtrapol[i].destroy();
2641 ✗ m_seqRNFluidExtrapol[i] = PetscVector::null();
2642 }
2643 ✗ m_isSeqRNFluidExtrapolAllocated = false;
2644 }
2645 11 }
2646
2647
2648 11 void LinearProblemFSINitscheXFEMFluid::copyOriginalSolToCurrentSol() {
2649
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
 11 if(!m_isSeqOriginalSolAllocated) {
2650 ✗ FEL_ERROR("The bdf storing all the data with the original supportDofMesh should already be created at this point");
2651 }
2652
2653 11 copyOriginalToCurrent(m_seqOriginalSolutions.sol_n(), sequentialSolution(), solution());
2654 11 }
2655
2656
2657 ✗ void LinearProblemFSINitscheXFEMFluid::copyOriginalBdfToCurrentBdf(felInt timeScheme, felInt extrapolOrder) {
2658 ✗ if(!m_isSeqOriginalSolAllocated) {
2659 ✗ FEL_ERROR("The bdf storing all the data with the original supportDofMesh should already be created at this point");
2660 }
2661
2662 // first, reinitialize the bdf with the solution to get vectors of right size
2663 ✗ felInt order = FelisceParam::instance().orderBdfNS;
2664
2665 // copy all bdf solution
2666 ✗ PetscVector zero = PetscVector::null();
2667 ✗ if(order == 1) {
2668 ✗ if(timeScheme == 3 && extrapolOrder == 2) {
2669 ✗ m_bdf->reinitialize(order, solution(), solution(), zero);
2670 ✗ copyOriginalToCurrent(m_seqOriginalSolutions.sol_n_1(), zero, m_bdf->sol_n_1());
2671 ✗ } else if ( ((timeScheme == 0) || (timeScheme == 2)) && extrapolOrder == 2) {
2672 ✗ m_bdf->reinitialize(order, solution(), solution(), zero);
2673 ✗ copyOriginalToCurrent(m_seqOriginalSolutions.sol_n_1(), zero, m_bdf->sol_n_1());
2674 }else {
2675 ✗ m_bdf->reinitialize(order, solution(), zero, zero);
2676 }
2677 ✗ } else if(order == 2) {
2678 ✗ m_bdf->reinitialize(order, solution(), solution(), zero);
2679 ✗ copyOriginalToCurrent(m_seqOriginalSolutions.sol_n_1(), zero, m_bdf->sol_n_1());
2680 } else {
2681 ✗ FEL_ERROR("Wrong bdf order");
2682 }
2683 }
2684
2685
2686
2687 2 void LinearProblemFSINitscheXFEMFluid::initOldObjects(felInt stabExplicitExtrapolOrder) {
2688 // compute how many old solution we need to store
2689
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 felInt order = std::max(stabExplicitExtrapolOrder, FelisceParam::instance().orderBdfNS);
2690
2691 // resize std::vector
2692
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_sequentialSolutionOld.resize(order);
2693
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_solutionOld.resize(order);
2694
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_aoOld.resize(order);
2695
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_supportDofUnknownOld.resize(order);
2696
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 m_dofOld.resize(order);
2697
2698
2699 // init old solution
2700
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 2 times.
 6 for(std::size_t i=0; i<m_sequentialSolutionOld.size(); ++i) {
2701
1/2
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
 4 m_sequentialSolutionOld[i].duplicateFrom(sequentialSolution());
2702
1/2
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
 4 m_sequentialSolutionOld[i].copyFrom(sequentialSolution());
2703 }
2704
2705
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 2 times.
 6 for(std::size_t i=0; i<m_solutionOld.size(); ++i) {
2706
1/2
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
 4 m_solutionOld[i].duplicateFrom(solution());
2707
1/2
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
 4 m_solutionOld[i].copyFrom(solution());
2708 }
2709
2710
2711 // init old AO
2712
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 2 times.
 6 for(std::size_t i=0; i<m_aoOld.size(); ++i) {
2713 4 m_aoOld[i] = ao();
2714
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 PetscObjectReference((PetscObject)ao());
2715 }
2716
2717
2718 // init old support dof
2719
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 2 times.
 6 for(std::size_t i=0; i<m_supportDofUnknownOld.size(); ++i)
2720
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 m_supportDofUnknownOld[i] = m_supportDofUnknown;
2721
2722
2723 // init old dof
2724
1/2
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
 2 std::vector<SupportDofMesh*> pSupportDofUnknown(m_supportDofUnknown.size(), nullptr);
2725
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 2 times.
 6 for(std::size_t i=0; i<m_dofOld.size(); ++i) {
2726
2/2
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 4 times.
 12 for(std::size_t iUnknown=0; iUnknown < m_supportDofUnknown.size(); ++iUnknown)
2727 8 pSupportDofUnknown[iUnknown] = &m_supportDofUnknownOld[i][iUnknown];
2728
2729
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 m_dofOld[i].setDof(m_listUnknown, m_listVariable, pSupportDofUnknown);
2730 4 m_dofOld[i].numDof() = dof().numDof();
2731
1/2
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
 4 m_dofOld[i].numDofPerUnknown() = dof().numDofPerUnknown();
2732 }
2733 2 }
2734
2735
2736 18 void LinearProblemFSINitscheXFEMFluid::updateOldSolution(felInt countInitOld) {
2737
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 2 times.
 18 if(countInitOld > 0) {
2738 // sequential solution
2739
1/2
✓ Branch 3 taken 16 times.
✗ Branch 4 not taken.
 16 m_sequentialSolutionOld[m_sequentialSolutionOld.size() - 1].destroy();
2740
2/2
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 16 times.
 32 for(std::size_t i=m_sequentialSolutionOld.size() - 1; i>0; --i)
2741
1/2
✓ Branch 3 taken 16 times.
✗ Branch 4 not taken.
 16 m_sequentialSolutionOld[i] = m_sequentialSolutionOld[i-1];
2742
2743
1/2
✓ Branch 3 taken 16 times.
✗ Branch 4 not taken.
 16 m_sequentialSolutionOld[0].duplicateFrom(sequentialSolution());
2744
1/2
✓ Branch 3 taken 16 times.
✗ Branch 4 not taken.
 16 m_sequentialSolutionOld[0].copyFrom(sequentialSolution());
2745
2746
2747 // parallel solution
2748
1/2
✓ Branch 3 taken 16 times.
✗ Branch 4 not taken.
 16 m_solutionOld[m_solutionOld.size() - 1].destroy();
2749
2/2
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 16 times.
 32 for(std::size_t i=m_solutionOld.size() - 1; i>0; --i)
2750
1/2
✓ Branch 3 taken 16 times.
✗ Branch 4 not taken.
 16 m_solutionOld[i] = m_solutionOld[i-1];
2751
2752
1/2
✓ Branch 3 taken 16 times.
✗ Branch 4 not taken.
 16 m_solutionOld[0].duplicateFrom(solution());
2753
1/2
✓ Branch 3 taken 16 times.
✗ Branch 4 not taken.
 16 m_solutionOld[0].copyFrom(solution());
2754
2755
2756 // AO
2757
1/2
✓ Branch 3 taken 16 times.
✗ Branch 4 not taken.
 16 AODestroy(&m_aoOld[m_aoOld.size() - 1]);
2758
2/2
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 16 times.
 32 for(std::size_t i=m_aoOld.size() - 1; i>0; --i)
2759 16 m_aoOld[i] = m_aoOld[i-1];
2760 16 m_aoOld[0] = ao();
2761
1/2
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
 16 PetscObjectReference((PetscObject)ao());
2762
2763
2764 // supportDofMesh (copy the old one)
2765
2/2
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 16 times.
 32 for(std::size_t i=m_supportDofUnknownOld.size() - 1; i>0; --i)
2766
1/2
✓ Branch 3 taken 16 times.
✗ Branch 4 not taken.
 16 m_supportDofUnknownOld[i] = m_supportDofUnknownOld[i-1];
2767
1/2
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
 16 m_supportDofUnknownOld[0] = m_supportDofUnknown;
2768
2769
2770 // dofs
2771
1/2
✓ Branch 3 taken 16 times.
✗ Branch 4 not taken.
 16 std::vector<SupportDofMesh*> pSupportDofUnknown(m_supportDofUnknown.size(), nullptr);
2772
2/2
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 16 times.
 32 for(std::size_t i=m_dofOld.size() - 1; i>0; --i) {
2773 16 m_dofOld[i].numDof() = m_dofOld[i-1].numDof();
2774
1/2
✓ Branch 5 taken 16 times.
✗ Branch 6 not taken.
 16 m_dofOld[i].numDofPerUnknown() = m_dofOld[i-1].numDofPerUnknown();
2775 }
2776
2777 16 m_dofOld[0].numDof() = dof().numDof();
2778
1/2
✓ Branch 5 taken 16 times.
✗ Branch 6 not taken.
 16 m_dofOld[0].numDofPerUnknown() = dof().numDofPerUnknown();
2779
2780
2781 // delete extrapolations
2782
2/2
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 16 times.
 32 for(std::size_t i=0; i<m_seqVelExtrapol.size(); ++i) {
2783
1/2
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
 16 m_seqVelExtrapol[i].destroy();
2784
2/4
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 16 times.
✗ Branch 6 not taken.
 16 m_seqVelExtrapol[i] = PetscVector::null();
2785 }
2786 16 m_isSeqVelExtrapolAllocated = false;
2787
2788 // delete rhs time
2789
2/2
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 16 times.
 32 for(std::size_t i=0; i<m_seqBdfRHS.size(); ++i) {
2790
1/2
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
 16 m_seqBdfRHS[i].destroy();
2791
2/4
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 16 times.
✗ Branch 6 not taken.
 16 m_seqBdfRHS[i] = PetscVector::null();
2792 }
2793 16 m_isSeqBdfRHSAllocated = false;
2794 16 } else {
2795 // copy the solutions
2796
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
 4 for(std::size_t i=m_sequentialSolutionOld.size() - 1; i>0; --i)
2797 2 m_sequentialSolutionOld[i].copyFrom(m_sequentialSolutionOld[i-1]);
2798 2 m_sequentialSolutionOld[0].copyFrom(sequentialSolution());
2799
2800
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
 4 for(std::size_t i=m_solutionOld.size() - 1; i>0; --i)
2801 2 m_solutionOld[i].copyFrom(m_solutionOld[i-1]);
2802 2 m_solutionOld[0].copyFrom(solution());
2803 }
2804 18 }
2805
2806
2807 /*!
2808 * Gather the extrapolation of the fluid velocity and pressure for the RN schemes
2809 * into m_seqStabExplicitFluidExtrapol.
2810 * If m_seqStabExplicitFluidExtrapol is null, it will be created and allocated during the gathering
2811 * according to the parallel shape of parFluidExtrapol
2812 */
2813 ✗ void LinearProblemFSINitscheXFEMFluid::gatherRNFluidExtrapol(std::vector<PetscVector>& parFluidExtrapol) {
2814 ✗ for(std::size_t i=0; i<parFluidExtrapol.size(); ++i)
2815 ✗ gatherVector(parFluidExtrapol[i], m_seqRNFluidExtrapol[i]);
2816
2817 ✗ m_isSeqRNFluidExtrapolAllocated = true;
2818 }
2819
2820 ✗ void LinearProblemFSINitscheXFEMFluid::gatherRNFluidExtrapol(PetscVector& parFluidExtrapol) {
2821 ✗ gatherVector(parFluidExtrapol, m_seqRNFluidExtrapol[0]);
2822 ✗ m_isSeqRNFluidExtrapolAllocated = true;
2823 }
2824
2825 ✗ void LinearProblemFSINitscheXFEMFluid::setRNFluidExtrapol(PetscVector& seqFluidExtrapol) {
2826 ✗ m_seqRNFluidExtrapol[0].copyFrom(seqFluidExtrapol);
2827 }
2828
2829 ✗ void LinearProblemFSINitscheXFEMFluid::initRNFluidExtrapol(felInt order) {
2830 ✗ if(m_isSeqRNFluidExtrapolAllocated) {
2831 ✗ for(std::size_t i=0; i<m_seqRNFluidExtrapol.size(); ++i)
2832 ✗ m_seqRNFluidExtrapol[i].destroy();
2833 }
2834
2835 ✗ m_seqRNFluidExtrapol.resize(order);
2836 ✗ for(std::size_t i=0; i<m_seqRNFluidExtrapol.size(); ++i)
2837 ✗ m_seqRNFluidExtrapol[i] = PetscVector::null();
2838 }
2839 }
2840