GCC Code Coverage Report

Directory: ./
File: Solver/linearProblemFSINitscheXFEMSolid.hxx
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 301 453 66.4%
Branches: 235 538 43.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: B. Fabreges
13 //
14
15 // System includes
16
17 // External includes
18
19 // Project includes
20 #include "FiniteElement/elementVector.hpp"
21 #include "FiniteElement/elementMatrix.hpp"
22
23 namespace felisce
24 {
25 /*!
26 \brief Matrix and rhs for the elastic string solver.
27 */
28 template<class SolidPrb>
29 2 LinearProblemFSINitscheXFEMSolid<SolidPrb>::LinearProblemFSINitscheXFEMSolid():
30
2/4
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
 2 SolidPrb()
31 {
32 2 m_isSeqLastTimeSolAllocated = false;
33 2 m_isSeqLastTimeVelSolAllocated = false;
34 2 m_isSeqStabExplicitFluidExtrapolAllocated = false;
35 2 m_isSeqRNStrucExtrapolAllocated = false;
36
37 2 m_elemFieldTimeSchemeType = NULL;
38 2 m_timeSchemeType = 0;
39
40 2 m_leftSide = 0;
41 2 m_rightSide = 0;
42
43 2 m_assemblingOfSolidProblemTerms = true;
44
45 2 m_useCurrentFluidSolution = false;
46
47
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 m_seqLastTimeSol = PetscVector::null();
48
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 m_seqLastTimeVelSol = PetscVector::null();
49
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 m_seqRNStrucExtrapol = PetscVector::null();
50
51
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 m_seqVecL2Norm = PetscVector::null();
52
53 2 m_isStrucOperatorExplicitlyComputed = false;
54 2 }
55
56 template<class SolidPrb>
57 4 LinearProblemFSINitscheXFEMSolid<SolidPrb>::~LinearProblemFSINitscheXFEMSolid()
58 {
59 4 if(m_isSeqLastTimeSolAllocated)
60 4 m_seqLastTimeSol.destroy();
61
62
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 if(m_isSeqLastTimeVelSolAllocated)
63 ✗ m_seqLastTimeVelSol.destroy();
64
65
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 if(m_isSeqStabExplicitFluidExtrapolAllocated)
66
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
 8 for(std::size_t i=0; i<m_seqStabExplicitFluidExtrapol.size(); ++i)
67 4 m_seqStabExplicitFluidExtrapol[i].destroy();
68
69
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 if(m_isSeqRNStrucExtrapolAllocated)
70 ✗ m_seqRNStrucExtrapol.destroy();
71
72
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
 4 if(m_seqVecL2Norm.isNotNull())
73 ✗ m_seqVecL2Norm.destroy();
74
75
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 2 times.
 12 for(std::size_t iPt=0; iPt < m_subItfPts.size(); ++iPt) {
76
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 8 delete m_subItfPts[iPt];
77 }
78 }
79
80 template<class SolidPrb>
81 2 void LinearProblemFSINitscheXFEMSolid<SolidPrb>::initialize(std::vector<GeometricMeshRegion::Pointer>& mesh, FelisceTransient::Pointer fstransient, MPI_Comm& comm, bool doUseSNES)
82 {
83 // call the initialize function of the solid problem
84
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 SolidPrb::initialize(mesh, fstransient, comm, doUseSNES);
85
86 // add the fluid variable to the list of physical variables
87 // unknowns are already define in the solid problem
88
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 this->m_listVariable.addVariable(velocity, 2);
89
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 this->m_listVariable.addVariable(pressure, 1);
90 2 this->m_listVariable.print(this->verbosity());
91
92 // send the unknown to the fluid problem (should be already initialize) to add them as variable there
93 2 m_pLinFluid->addSolidVariable();
94
95 // get the id of the displacement variable
96 2 m_iDisplacement = this->m_listVariable.getVariableIdList(displacement);
97
98 // get the number of component of the displacement
99 2 m_numSolidComponent = this->m_listVariable[m_iDisplacement].numComponent();
100
101 // get the time scheme type from the data file
102 2 m_elemFieldTimeSchemeType = FelisceParam::instance().elementFieldDynamicValue("TimeSchemeType");
103
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if(m_elemFieldTimeSchemeType != NULL)
104 2 m_timeSchemeType = (felInt) m_elemFieldTimeSchemeType->constant(Comp1);
105
106
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if (m_timeSchemeType == 2) {
107 // Robin-Neumann scheme
108 // extrapolation order
109 ✗ m_betaRN.resize(2);
110 ✗ if (FelisceParam::instance().fsiRNscheme == 0){
111 ✗ m_betaRN[0] = 0.;
112 ✗ m_betaRN[1] = 0.;
113 ✗ } else if (FelisceParam::instance().fsiRNscheme == 1){
114 ✗ m_betaRN[0] = 1.;
115 ✗ m_betaRN[1] = 0.;
116 ✗ } else if (FelisceParam::instance().fsiRNscheme == 2){
117 ✗ m_betaRN[0] = 2.;
118 ✗ m_betaRN[1] = -1.;
119 }
120 }
121
122 // sub element
123 2 m_numVerPerSolidElt = SolidPrb::mesh()->domainDim() + 1;
124
125 // create the std::vector of points for the solid elements
126 2 m_subItfPts.resize(m_numVerPerSolidElt);
127
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 2 times.
 6 for(std::size_t iver=0; iver < m_subItfPts.size(); ++iver) {
128 8 m_subItfPts[iver] = new Point();
129 }
130
131 // compute solid coefficient
132 2 m_massVel = FelisceParam::instance().densitySolid * FelisceParam::instance().thickness / this->m_fstransient->timeStep;
133 2 }
134
135 template<class SolidPrb>
136 80 void LinearProblemFSINitscheXFEMSolid<SolidPrb>::initPerElementTypeBD(GeometricMeshRegion::ElementType eltType, FlagMatrixRHS flagMatrixRHS)
137 {
138 // call the function of the solid problem
139 80 SolidPrb::initPerElementTypeBD(eltType, flagMatrixRHS);
140
141 // m_dispLastSol contains some old displacement
142 80 m_dispLastSol.initialize(DOF_FIELD, *(this->m_feDisp), m_numSolidComponent);
143
144 // source term coming from the coupling with the fluid and Nitsche-XFEM terms
145 80 m_elemNitscheXFEMSolidRHS.initialize(QUAD_POINT_FIELD, *(this->m_feDisp), m_numSolidComponent);
146
147 // allocate elementField for the fluid
148 80 m_fluidDofSol.resize(2);
149 80 m_fluidDofSol[0].initialize(DOF_FIELD, *m_feFluid, m_feFluid->numCoor());
150 80 m_fluidDofSol[1].initialize(DOF_FIELD, *m_fePres, 1);
151
152 // integration points
153 80 m_intPoints.clear();
154 80 m_intPoints.resize(this->m_feDisp->numQuadraturePoint());
155 80 }
156
157 template<class SolidPrb>
158 720 void LinearProblemFSINitscheXFEMSolid<SolidPrb>::computeElementArrayBD(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint, felInt& iel, FlagMatrixRHS flagMatrixRHS)
159 {
160 // Assembling of the solid problem terms
161
2/2
✓ Branch 0 taken 360 times.
✓ Branch 1 taken 360 times.
 720 if(m_assemblingOfSolidProblemTerms) {
162
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 360 times.
 360 if(!this->m_feDisp->hasOriginalQuadPoint())
163 ✗ this->m_feDisp->updateBasisAndNormalContra(0, elemPoint);
164
165
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 360 times.
 360 if(m_isStrucOperatorExplicitlyComputed) {
166 // update finite element
167 ✗ this->m_feDisp->updateMeasNormalContra(0, elemPoint);
168
169 // MATRIX //
170 // time scheme
171 ✗ this->m_elementMatBD[0]->phi_i_phi_j(m_massVel, *(this->m_feDisp), 0, 0, m_numSolidComponent);
172
173 // RHS //
174 // time scheme
175 ✗ m_dispLastSol.setValue(m_seqRNStrucExtrapol, *(this->m_feDisp), iel, m_iDisplacement, this->m_ao, this->dof());
176 ✗ this->m_elementVectorBD[0]->source(m_massVel, *(this->m_feDisp), m_dispLastSol, 0, m_numSolidComponent);
177
178 } else {
179 360 SolidPrb::computeElementArrayBD(elemPoint, elemIdPoint, iel, flagMatrixRHS);
180 }
181 } else {
182 360 assembleNitscheXFEMTerms(elemPoint, elemIdPoint, iel, flagMatrixRHS);
183 }
184 720 }
185
186 template<class SolidPrb>
187 720 void LinearProblemFSINitscheXFEMSolid<SolidPrb>::computeElementArrayBD(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint, const std::vector<Point*>& elemNormal, const std::vector< std::vector<Point*> >& elemTangent, felInt& iel, FlagMatrixRHS flagMatrixRHS)
188 {
189 // Assembling of the solid problem terms
190
2/2
✓ Branch 0 taken 360 times.
✓ Branch 1 taken 360 times.
 720 if(m_assemblingOfSolidProblemTerms) {
191
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 360 times.
 360 if(!this->m_feDisp->hasOriginalQuadPoint())
192 ✗ this->m_feDisp->updateBasisAndNormalContra(0, elemPoint);
193
194
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 360 times.
 360 if(m_isStrucOperatorExplicitlyComputed) {
195 // update finite element
196 ✗ this->m_feDisp->updateMeasNormalContra(0, elemPoint);
197
198 // MATRIX //
199 // time scheme
200 ✗ this->m_elementMatBD[0]->phi_i_phi_j(m_massVel, *(this->m_feDisp), 0, 0, m_numSolidComponent);
201
202 // RHS //
203 // time scheme
204 ✗ m_dispLastSol.setValue(m_seqRNStrucExtrapol, *(this->m_feDisp), iel, m_iDisplacement, this->m_ao, this->dof());
205 ✗ this->m_elementVectorBD[0]->source(m_massVel, *(this->m_feDisp), m_dispLastSol, 0, m_numSolidComponent);
206 } else {
207 360 SolidPrb::computeElementArrayBD(elemPoint, elemIdPoint, elemNormal, elemTangent, iel, flagMatrixRHS);
208 }
209 } else {
210 360 assembleNitscheXFEMTerms(elemPoint, elemIdPoint, iel, flagMatrixRHS);
211 }
212 720 }
213
214
215
216 template<class SolidPrb>
217 720 void LinearProblemFSINitscheXFEMSolid<SolidPrb>::assembleNitscheXFEMTerms(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint, felInt& iel, FlagMatrixRHS flagMatrixRHS)
218 {
219 (void) elemIdPoint;
220 (void) flagMatrixRHS;
221
222 // Assembling of the terms coming from the coupling with the fluid and the Nitsche-XFEM formulation
223 720 std::vector< std::vector<felInt> > vectorIdSupport;
224
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 720 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 720 if((m_timeSchemeType == 3) || (m_timeSchemeType == 2) || (m_timeSchemeType == 4) )
225
1/2
✓ Branch 2 taken 720 times.
✗ Branch 3 not taken.
 720 vectorIdSupport.resize(m_seqStabExplicitFluidExtrapol.size());
226 else
227 ✗ vectorIdSupport.resize(1);
228
229 // loop over the sub structure to compute the integrals involving the fluid velocity and pressure
230
1/2
✓ Branch 1 taken 720 times.
✗ Branch 2 not taken.
 720 felInt numSubElement = m_duplicateSupportElements->getNumSubItfEltPerItfElt(iel);
231
1/2
✓ Branch 0 taken 720 times.
✗ Branch 1 not taken.
 720 if(numSubElement > 0) {
232 felInt subElemId;
233 felInt fluidElemId;
234
1/2
✓ Branch 2 taken 720 times.
✗ Branch 3 not taken.
 720 std::vector<double> tmpPt(3);
235
1/2
✓ Branch 2 taken 720 times.
✗ Branch 3 not taken.
 720 std::vector<double> itfNormal(3);
236
1/2
✓ Branch 2 taken 720 times.
✗ Branch 3 not taken.
 720 std::vector<double> itfNormalInit(3);
237
238 // get the normal and the initial normal to this sub element
239 720 GeometricMeshRegion::ElementType eltTypeStruc = this->mesh()->bagElementTypeDomain()[0];
240
1/2
✓ Branch 4 taken 720 times.
✗ Branch 5 not taken.
 720 this->mesh()->listElementNormal(eltTypeStruc, iel).getCoor(itfNormal);
241
1/2
✓ Branch 4 taken 720 times.
✗ Branch 5 not taken.
 720 this->mesh()->listElementNormalInit(eltTypeStruc, iel).getCoor(itfNormalInit);
242
243
2/2
✓ Branch 2 taken 2056 times.
✓ Branch 3 taken 720 times.
 2776 for(felInt iSubElt=0; iSubElt < numSubElement; ++iSubElt) {
244 // get the id of the sub element
245
1/2
✓ Branch 1 taken 2056 times.
✗ Branch 2 not taken.
 2056 subElemId = m_duplicateSupportElements->getSubItfEltIdxItf(iel, iSubElt);
246
247 // fill the std::vector with the points of the sub-element
248
2/2
✓ Branch 0 taken 4112 times.
✓ Branch 1 taken 2056 times.
 6168 for(std::size_t iPt=0; iPt < m_numVerPerSolidElt; ++iPt) {
249
1/2
✓ Branch 1 taken 4112 times.
✗ Branch 2 not taken.
 4112 m_duplicateSupportElements->getSubItfEltVerCrd(subElemId, iPt, tmpPt);
250 4112 *m_subItfPts[iPt] = Point(tmpPt[0], tmpPt[1], tmpPt[2]);
251 }
252
253 // update measure
254
2/4
✓ Branch 1 taken 2056 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2056 times.
✗ Branch 5 not taken.
 2056 this->m_feDisp->updateSubElementMeasNormal(0, elemPoint, m_subItfPts);
255
256
257 // get the coordinate of the fluid element
258
1/2
✓ Branch 1 taken 2056 times.
✗ Branch 2 not taken.
 2056 fluidElemId = m_duplicateSupportElements->getMshEltIdxOfSubItfElt(subElemId);
259
260 2056 GeometricMeshRegion::ElementType eltType = m_pLinFluid->mesh()->bagElementTypeDomain()[0];
261 2056 int numPointsPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
262
263
1/2
✓ Branch 2 taken 2056 times.
✗ Branch 3 not taken.
 2056 std::vector<Point*> fluidElemPts(numPointsPerElt);
264
1/2
✓ Branch 2 taken 2056 times.
✗ Branch 3 not taken.
 2056 std::vector<felInt> fluidElemPtsId(numPointsPerElt);
265
266
1/2
✓ Branch 3 taken 2056 times.
✗ Branch 4 not taken.
 2056 m_pLinFluid->mesh()->getOneElement(eltType, fluidElemId, fluidElemPtsId, 0);
267
2/2
✓ Branch 0 taken 6168 times.
✓ Branch 1 taken 2056 times.
 8224 for (int iPoint = 0; iPoint < numPointsPerElt; ++iPoint)
268 6168 fluidElemPts[iPoint] = &m_pLinFluid->mesh()->listPoints()[fluidElemPtsId[iPoint]];
269
270
271 // update the element for the fluid finite element
272
1/2
✓ Branch 1 taken 2056 times.
✗ Branch 2 not taken.
 2056 m_feFluid->updateFirstDeriv(0, fluidElemPts);
273
274 // compute the integration point for this element
275
1/2
✓ Branch 1 taken 2056 times.
✗ Branch 2 not taken.
 2056 computeIntegrationPoints(fluidElemPts);
276
277 // get all the support elements for the fluid
278
6/8
✓ Branch 1 taken 2056 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2056 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1028 times.
✓ Branch 6 taken 1028 times.
✓ Branch 7 taken 1028 times.
✓ Branch 8 taken 1028 times.
 2056 if(FelisceParam::instance().useProjectionForNitscheXFEM || m_useCurrentFluidSolution)
279
2/2
✓ Branch 1 taken 1028 times.
✓ Branch 2 taken 1028 times.
 2056 for(std::size_t i=0; i<vectorIdSupport.size(); ++i)
280
1/2
✓ Branch 3 taken 1028 times.
✗ Branch 4 not taken.
 1028 m_pLinFluid->supportDofUnknown(0).getIdElementSupport(fluidElemId, vectorIdSupport[i]);
281 else
282
2/2
✓ Branch 1 taken 2056 times.
✓ Branch 2 taken 1028 times.
 3084 for(std::size_t i=0; i<vectorIdSupport.size(); ++i)
283
1/2
✓ Branch 3 taken 2056 times.
✗ Branch 4 not taken.
 2056 m_pLinFluid->supportDofUnknownOld(i,0).getIdElementSupport(fluidElemId, vectorIdSupport[i]);
284
285 // compute the elementary matrix and rhs
286
1/6
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2056 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 2056 switch(m_timeSchemeType) {
287 ✗ case 0: {
288 // Robin-Neumann schemes (direct implementation) (0th, 1st and 2nd order extrapolated variants)
289 ✗ for(std::size_t iSup=m_leftSide; iSup<vectorIdSupport[0].size() - m_rightSide; ++iSup) {
290 // force of the fluid to the structure
291 ✗ updateFluidVel(vectorIdSupport[0][iSup], m_pLinFluid->sequentialSolution());
292
293 ✗ computeFluidNormalDerivativeJump(2.*iSup - 1., itfNormal, itfNormalInit);
294 ✗ this->m_elementVectorBD[0]->source(-1., *(this->m_feDisp), m_elemNitscheXFEMSolidRHS, 0, m_numSolidComponent);
295 }
296 ✗ break;
297 }
298
299 ✗ case 1: {
300 // Robin - Robin //
301 // Nitsche penalty term (derivative of the displacement)
302 // This integral is here because the penalty parameter depends of the fluid element
303 ✗ double penaltyTerm = FelisceParam::instance().viscosity*FelisceParam::instance().NitschePenaltyParam/m_feFluid->diameter();
304 ✗ double tau = this->m_fstransient->timeStep;
305 ✗ double coef = 2. - m_leftSide - m_rightSide;
306
307 ✗ this->m_elementMatBD[0]->phi_i_phi_j(coef*penaltyTerm/tau, *(this->m_feDisp), 0, 0, m_numSolidComponent);
308
309 ✗ m_dispLastSol.setValue(m_seqLastTimeSol, *(this->m_feDisp), iel, m_iDisplacement, this->m_ao, this->dof());
310 ✗ this->m_elementVectorBD[0]->source(coef*penaltyTerm/tau, *(this->m_feDisp), m_dispLastSol, 0, m_numSolidComponent);
311
312
313 ✗ for(std::size_t iSup=m_leftSide; iSup<vectorIdSupport[0].size() - m_rightSide; ++iSup) {
314 // get the velocity and pressure at the dof of the structure
315 // Nitsche penalty term
316 ✗ if(FelisceParam::instance().useProjectionForNitscheXFEM)
317 ✗ updateFluidVel(vectorIdSupport[0][iSup], m_pLinFluid->sequentialSolution());
318 else
319 ✗ updateFluidVel(vectorIdSupport[0][iSup], m_pLinFluid->sequentialSolutionOld(0));
320
321 ✗ computeFluidVelOnSubStructure(itfNormalInit);
322 ✗ this->m_elementVectorBD[0]->source(penaltyTerm, *(this->m_feDisp), m_elemNitscheXFEMSolidRHS, 0, m_numSolidComponent);
323
324 ✗ computeFluidNormalDerivativeJump(2.*iSup - 1., itfNormal, itfNormalInit);
325 ✗ this->m_elementVectorBD[0]->source(-1., *(this->m_feDisp), m_elemNitscheXFEMSolidRHS, 0, m_numSolidComponent);
326 }
327 ✗ break;
328 }
329
330 ✗ case 2: {
331 // Robin-Neumann schemes a la Stenberg (0th, 1st and 2nd order extrapolated variants)
332 ✗ double alpha = 0.5;
333 ✗ double RNParam = FelisceParam::instance().densitySolid * FelisceParam::instance().thickness / this->m_fstransient->timeStep;
334 ✗ double penaltyParam1 = (m_feFluid->diameter()) / (FelisceParam::instance().viscosity * FelisceParam::instance().NitschePenaltyParam + RNParam*m_feFluid->diameter());
335 ✗ double penaltyParam2 = (FelisceParam::instance().NitschePenaltyParam *FelisceParam::instance().viscosity ) / (FelisceParam::instance().viscosity *FelisceParam::instance().NitschePenaltyParam + RNParam*m_feFluid->diameter());
336
337 ✗ if (m_useCurrentFluidSolution) {
338 ✗ m_dispLastSol.setValue(m_seqRNStrucExtrapol, *(this->m_feDisp), iel, m_iDisplacement, this->m_ao, this->dof());
339 ✗ this->m_elementVectorBD[0]->source(-1.0 * penaltyParam2 * RNParam , *(this->m_feDisp), m_dispLastSol, 0, m_numSolidComponent);
340
341 // get the velocity and pressure at the dof of the structure
342 // Nitsche penalty term
343 ✗ for(std::size_t iSup=m_leftSide; iSup< vectorIdSupport[0].size() - m_rightSide; ++iSup) {
344 ✗ updateFluidVel(vectorIdSupport[0][iSup], m_pLinFluid->sequentialSolution());
345
346 ✗ computeFluidVelOnSubStructure(itfNormalInit);
347 ✗ this->m_elementVectorBD[0]->source(alpha * penaltyParam2 * RNParam, *(this->m_feDisp), m_elemNitscheXFEMSolidRHS, 0, m_numSolidComponent);
348
349 // force of the fluid to the structure
350 ✗ computeFluidNormalDerivativeJump(2.*iSup - 1., itfNormal, itfNormalInit);
351 ✗ this->m_elementVectorBD[0]->source(-1.0 * penaltyParam1 * RNParam , *(this->m_feDisp), m_elemNitscheXFEMSolidRHS, 0, m_numSolidComponent);
352 }
353 } else {
354 // get the velocity and pressure at the dof of the structure
355 // Nitsche penalty term
356 ✗ for(std::size_t i=0; i<m_seqStabExplicitFluidExtrapol.size(); ++i) {
357 ✗ for(std::size_t iSup=m_leftSide; iSup<vectorIdSupport[i].size() - m_rightSide; ++iSup) {
358 ✗ updateFluidVel(vectorIdSupport[i][iSup], m_seqStabExplicitFluidExtrapol[i], i);
359 ✗ computeFluidNormalDerivativeJump(2.*iSup - 1., itfNormal, itfNormalInit);
360 ✗ this->m_elementVectorBD[0]->source(-1.0 * penaltyParam2 * m_betaRN[i], *(this->m_feDisp), m_elemNitscheXFEMSolidRHS, 0, m_numSolidComponent);
361 }
362 }
363 }
364
365 ✗ break;
366 }
367
368 2056 case 3:{
369 // stabilized explicit //
370 // Nitsche penalty term (derivative of the displacement)
371 // This integral is here because the penalty parameter depends of the fluid element
372
3/6
✓ Branch 1 taken 2056 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2056 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2056 times.
✗ Branch 8 not taken.
 2056 double penaltyTerm = FelisceParam::instance().viscosity*FelisceParam::instance().NitschePenaltyParam/m_feFluid->diameter();
373 2056 double tau = this->m_fstransient->timeStep;
374 2056 double coef = 2. - m_leftSide - m_rightSide;
375
376 // Add the penalty terms in the solid if we are using the classic Nitsche's formulation
377
2/4
✓ Branch 1 taken 2056 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2056 times.
✗ Branch 4 not taken.
 2056 if(!FelisceParam::instance().usePenaltyFreeNitscheMethod) {
378 // displacement part
379
1/2
✓ Branch 3 taken 2056 times.
✗ Branch 4 not taken.
 2056 this->m_elementMatBD[0]->phi_i_phi_j(coef*penaltyTerm/tau, *(this->m_feDisp), 0, 0, m_numSolidComponent); // gamma/(h tau) \int d . w
380
1/2
✓ Branch 2 taken 2056 times.
✗ Branch 3 not taken.
 2056 m_dispLastSol.setValue(m_seqLastTimeSol, *(this->m_feDisp), iel, m_iDisplacement, this->m_ao, this->dof());
381
1/2
✓ Branch 3 taken 2056 times.
✗ Branch 4 not taken.
 2056 this->m_elementVectorBD[0]->source(coef*penaltyTerm/tau, *(this->m_feDisp), m_dispLastSol, 0, m_numSolidComponent); // gamma/(h tau) \int dold . w
382
383 // get the velocity and pressure at the dof of the structure
384 // Fluid Nitsche penalty term
385
2/2
✓ Branch 1 taken 3084 times.
✓ Branch 2 taken 2056 times.
 5140 for(std::size_t i=0; i<m_seqStabExplicitFluidExtrapol.size(); ++i) {
386
2/2
✓ Branch 2 taken 6106 times.
✓ Branch 3 taken 3084 times.
 9190 for(std::size_t iSup=m_leftSide; iSup<vectorIdSupport[i].size() - m_rightSide; ++iSup) {
387
1/2
✓ Branch 4 taken 6106 times.
✗ Branch 5 not taken.
 6106 updateFluidVel(vectorIdSupport[i][iSup], m_seqStabExplicitFluidExtrapol[i], i);
388
1/2
✓ Branch 1 taken 6106 times.
✗ Branch 2 not taken.
 6106 computeFluidVelOnSubStructure(itfNormalInit);
389
1/2
✓ Branch 3 taken 6106 times.
✗ Branch 4 not taken.
 6106 this->m_elementVectorBD[0]->source(penaltyTerm, *(this->m_feDisp), m_elemNitscheXFEMSolidRHS, 0, m_numSolidComponent); // gamma/(h ) \int uf^{n-1} . w
390 }
391 }
392 }
393
394 // Nitsche consistency term
395
2/2
✓ Branch 2 taken 4090 times.
✓ Branch 3 taken 2056 times.
 6146 for(std::size_t iSup=m_leftSide; iSup<vectorIdSupport[0].size() - m_rightSide; ++iSup) {
396
6/8
✓ Branch 1 taken 4090 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4090 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2056 times.
✓ Branch 6 taken 2034 times.
✓ Branch 7 taken 2056 times.
✓ Branch 8 taken 2034 times.
 4090 if(FelisceParam::instance().useProjectionForNitscheXFEM || m_useCurrentFluidSolution)
397
1/2
✓ Branch 4 taken 2056 times.
✗ Branch 5 not taken.
 2056 updateFluidVel(vectorIdSupport[0][iSup], m_pLinFluid->sequentialSolution());
398 else
399
1/2
✓ Branch 4 taken 2034 times.
✗ Branch 5 not taken.
 2034 updateFluidVel(vectorIdSupport[0][iSup], m_pLinFluid->sequentialSolutionOld(0));
400
401
1/2
✓ Branch 1 taken 4090 times.
✗ Branch 2 not taken.
 4090 computeFluidNormalDerivativeJump(2.*iSup - 1., itfNormal, itfNormalInit);
402
1/2
✓ Branch 3 taken 4090 times.
✗ Branch 4 not taken.
 4090 this->m_elementVectorBD[0]->source(-1., *(this->m_feDisp), m_elemNitscheXFEMSolidRHS, 0, m_numSolidComponent); // -int sigma(u^{n-1},p^{n-1})n . w
403 }
404 2056 break;
405 }
406
407
408 ✗ case 4:{
409 // Robin-Neumann semi-implict (0th, 1st and 2nd order extrapolated variants) //
410 // Nitsche penalty term (derivative of the displacement)
411 // This integral is here because the penalty parameter depends of the fluid element
412 ✗ double penaltyTerm = FelisceParam::instance().viscosity*FelisceParam::instance().NitschePenaltyParam/m_feFluid->diameter();
413 ✗ double coef = 2. - m_leftSide - m_rightSide;
414
415 ✗ if(!m_isStrucOperatorExplicitlyComputed) {
416 ✗ m_dispLastSol.setValue(this->sequentialSolution(), *(this->m_feDisp), iel, m_iDisplacement, this->m_ao, this->dof());
417 ✗ this->m_elementVectorBD[0]->source(-coef*penaltyTerm, *(this->m_feDisp), m_dispLastSol, 0, m_numSolidComponent);
418 } else {
419 // matrix part
420 ✗ this->m_elementMatBD[0]->phi_i_phi_j(coef * penaltyTerm, *(this->m_feDisp), 0, 0, m_numSolidComponent);
421
422 // All the extrapolation parts are added in the rhs later in the model.
423 // There are indeed store in a std::vector during the real solid step and reused here later.
424 }
425
426 // get the velocity and pressure at the dof of the structure
427 // Nitsche penalty term
428 ✗ for(std::size_t i=0; i<m_seqStabExplicitFluidExtrapol.size(); ++i) {
429 ✗ for(std::size_t iSup=m_leftSide; iSup<vectorIdSupport[i].size() - m_rightSide; ++iSup) {
430 ✗ updateFluidVel(vectorIdSupport[i][iSup], m_seqStabExplicitFluidExtrapol[i], i);
431 ✗ computeFluidVelOnSubStructure(itfNormalInit);
432 ✗ this->m_elementVectorBD[0]->source(penaltyTerm, *(this->m_feDisp), m_elemNitscheXFEMSolidRHS, 0, m_numSolidComponent);
433 }
434 }
435
436 // Nitsche consistency term
437 ✗ for(std::size_t iSup=m_leftSide; iSup<vectorIdSupport[0].size() - m_rightSide; ++iSup) {
438 ✗ if(FelisceParam::instance().useProjectionForNitscheXFEM || m_useCurrentFluidSolution)
439 ✗ updateFluidVel(vectorIdSupport[0][iSup], m_pLinFluid->sequentialSolution());
440 else
441 ✗ updateFluidVel(vectorIdSupport[0][iSup], m_pLinFluid->sequentialSolutionOld(0));
442
443 ✗ computeFluidNormalDerivativeJump(2.*iSup - 1., itfNormal, itfNormalInit);
444 ✗ this->m_elementVectorBD[0]->source(-1., *(this->m_feDisp), m_elemNitscheXFEMSolidRHS, 0, m_numSolidComponent);
445 }
446
447 ✗ break;
448 }
449
450 ✗ default:{
451 ✗ FEL_ERROR("Bad time splitting scheme");
452 ✗ break;
453 }
454 }
455 }
456 720 }
457 720 }
458
459 template<class SolidPrb>
460 2 void LinearProblemFSINitscheXFEMSolid<SolidPrb>::setDuplicateSupportObject(DuplicateSupportDof* object)
461 {
462 2 m_duplicateSupportElements = object;
463 2 }
464
465 template<class SolidPrb>
466 2 void LinearProblemFSINitscheXFEMSolid<SolidPrb>::setFluidLinPrb(LinearProblemFSINitscheXFEMFluid* pLinPrb)
467 {
468 2 m_pLinFluid = pLinPrb;
469 2 }
470
471 template<class SolidPrb>
472 2 void LinearProblemFSINitscheXFEMSolid<SolidPrb>::setFluidFiniteElement(CurrentFiniteElement* velFE, CurrentFiniteElement* preFE)
473 {
474 2 m_feFluid = velFE;
475 2 m_fePres = preFE;
476 2 }
477
478
479 template<class SolidPrb>
480 80 void LinearProblemFSINitscheXFEMSolid<SolidPrb>::isAssemblingOfSolidProblemTerms(bool isit)
481 {
482 80 m_assemblingOfSolidProblemTerms = isit;
483 80 }
484
485 template<class SolidPrb>
486 ✗ void LinearProblemFSINitscheXFEMSolid<SolidPrb>::isStrucOperatorExplicitlyComputed(bool isit)
487 {
488 ✗ m_isStrucOperatorExplicitlyComputed = isit;
489 }
490
491 template<class SolidPrb>
492 40 void LinearProblemFSINitscheXFEMSolid<SolidPrb>::setUserBoundaryCondition()
493 {
494
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 36 times.
 40 if(this->m_fstransient->iteration == 1) {
495 // resize the arrays
496
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 m_applyDirichletBC.resize(this->m_listUnknown.size());
497
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 m_idDofBCVertices.resize(this->m_listUnknown.size());
498
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 m_valBC.resize(this->m_listUnknown.size());
499
500 // get the id of the boundary of the structure
501 // we identify the points that are in only one edge. There should only be "boundaries" elements.
502
1/2
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 4 std::vector<felInt> numEdgePerVertices(this->m_supportDofUnknown[0].numSupportDof(), 0);
503 4 felInt cntDof = 0;
504 4 felInt numComp = 0;
505 4 felInt idof = 0;
506
507
2/2
✓ Branch 3 taken 72 times.
✓ Branch 4 taken 4 times.
 76 for(std::size_t ielSup=0; ielSup<this->m_supportDofUnknown[0].iEle().size()-1; ++ielSup) {
508
2/2
✓ Branch 2 taken 144 times.
✓ Branch 3 taken 72 times.
 216 for(felInt idSup=0; idSup<this->m_supportDofUnknown[0].getNumSupportDof(ielSup); ++idSup) {
509 144 ++numEdgePerVertices[this->m_supportDofUnknown[0].iSupportDof()[this->m_supportDofUnknown[0].iEle()[ielSup] + idSup]];
510 }
511 }
512
513
2/2
✓ Branch 1 taken 80 times.
✓ Branch 2 taken 4 times.
 84 for(std::size_t ivert=0; ivert<numEdgePerVertices.size(); ++ivert) {
514
2/2
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 64 times.
 80 if(numEdgePerVertices[ivert] == 1) {
515 // we got a point on the border of the structure
516 // add it only if it's not a front point (not attached, to change if one wants to modeled a rigid immersed wall)
517
3/4
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 8 times.
 16 if(!m_duplicateSupportElements->checkFrontPoint(ivert)) {
518 8 cntDof = 0;
519
2/2
✓ Branch 1 taken 12 times.
✓ Branch 2 taken 8 times.
 20 for(std::size_t iUnknown=0; iUnknown<this->m_listUnknown.size(); ++iUnknown) {
520 12 numComp = this->m_listVariable[this->m_listUnknown.idVariable(iUnknown)].numComponent();
521
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 12 times.
 28 for(felInt icoor=0; icoor<numComp; ++icoor) {
522 16 idof = cntDof + ivert * numComp + icoor;
523
1/2
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
 16 AOApplicationToPetsc(this->ao(), 1, &idof);
524
1/2
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
 16 m_idDofBCVertices[iUnknown].push_back(idof);
525 }
526
527 12 cntDof += this->supportDofUnknown(iUnknown).numSupportDof() * numComp;
528 }
529 }
530 }
531 }
532
533 // By default, apply an homogeneous Dirichlet BC to all unknown of the solid.
534
2/2
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 4 times.
 10 for(std::size_t iUnknown=0; iUnknown<this->m_listUnknown.size(); ++iUnknown) {
535
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 m_applyDirichletBC[iUnknown] = true;
536
1/2
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
 6 m_valBC[iUnknown].resize(m_idDofBCVertices[iUnknown].size(), 0.);
537 }
538 4 }
539
540
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 if(!m_idDofBCVertices.empty()) {
541 40 userSetDirichletBC();
542
543 40 bool isDirichletBC = false;
544
2/2
✓ Branch 1 taken 60 times.
✓ Branch 2 taken 40 times.
 100 for(std::size_t iUnknown=0; iUnknown<this->m_listUnknown.size(); ++iUnknown) {
545
1/2
✓ Branch 2 taken 60 times.
✗ Branch 3 not taken.
 60 if(m_applyDirichletBC[iUnknown]) {
546 // Enforce the Dirichlet boundary conditions
547 60 this->matrix(0).zeroRows(m_idDofBCVertices[iUnknown].size(), m_idDofBCVertices[iUnknown].data(), 1.);
548
549 60 this->vector().setValues(m_idDofBCVertices[iUnknown].size(), m_idDofBCVertices[iUnknown].data(), m_valBC[iUnknown].data(), INSERT_VALUES);
550
551 60 isDirichletBC = true;
552 }
553 }
554
555
1/2
✓ Branch 0 taken 40 times.
✗ Branch 1 not taken.
 40 if(isDirichletBC)
556 40 this->vector().assembly();
557 }
558 40 }
559
560
561 template<class SolidPrb>
562 ✗ void LinearProblemFSINitscheXFEMSolid<SolidPrb>::userSetDirichletBC()
563 {
564 // nothing to do, to override in the user solid problem
565 }
566
567
568 template<class SolidPrb>
569 10196 void LinearProblemFSINitscheXFEMSolid<SolidPrb>::updateFluidVel(felInt ielSupportDof, const PetscVector& fluidSol, felInt iold)
570 {
571 felInt idDof;
572
573
2/2
✓ Branch 2 taken 20392 times.
✓ Branch 3 taken 10196 times.
 30588 for(std::size_t iUnknown=0; iUnknown<m_pLinFluid->listUnknown().size(); ++iUnknown) {
574 20392 felInt idVar = m_pLinFluid->listUnknown().idVariable(iUnknown);
575
576
6/8
✓ Branch 1 taken 20392 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 20392 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 8224 times.
✓ Branch 6 taken 12168 times.
✓ Branch 7 taken 8224 times.
✓ Branch 8 taken 12168 times.
 20392 if(FelisceParam::instance().useProjectionForNitscheXFEM || m_useCurrentFluidSolution) {
577 // get the value of the velocity and the pressure at the dofs
578
2/2
✓ Branch 2 taken 24672 times.
✓ Branch 3 taken 8224 times.
 32896 for (felInt iSupport=0; iSupport < m_pLinFluid->supportDofUnknown(iUnknown).getNumSupportDof(ielSupportDof); iSupport++) {
579
2/2
✓ Branch 3 taken 37008 times.
✓ Branch 4 taken 24672 times.
 61680 for(std::size_t iComp=0; iComp<m_pLinFluid->listVariable()[idVar].numComponent(); ++iComp) {
580 // local to global id of the dof
581
1/2
✓ Branch 2 taken 37008 times.
✗ Branch 3 not taken.
 37008 m_pLinFluid->dof().loc2glob(ielSupportDof, iSupport, idVar, iComp, idDof);
582
583 // petsc ordering
584
1/2
✓ Branch 2 taken 37008 times.
✗ Branch 3 not taken.
 37008 AOApplicationToPetsc(m_pLinFluid->ao(), 1, &idDof);
585
586 // get values
587
2/4
✓ Branch 2 taken 37008 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 37008 times.
✗ Branch 6 not taken.
 37008 fluidSol.getValues( 1, &idDof, &m_fluidDofSol[iUnknown].val(iComp, iSupport));
588 }
589 }
590 } else {
591 // get the value of the velocity and the pressure at the dofs
592
2/2
✓ Branch 2 taken 36504 times.
✓ Branch 3 taken 12168 times.
 48672 for (felInt iSupport=0; iSupport < m_pLinFluid->supportDofUnknownOld(iold, iUnknown).getNumSupportDof(ielSupportDof); iSupport++) {
593
2/2
✓ Branch 3 taken 54756 times.
✓ Branch 4 taken 36504 times.
 91260 for(std::size_t iComp=0; iComp<m_pLinFluid->listVariable()[idVar].numComponent(); ++iComp) {
594 // local to global id of the dof
595
1/2
✓ Branch 2 taken 54756 times.
✗ Branch 3 not taken.
 54756 m_pLinFluid->dofOld(iold).loc2glob(ielSupportDof, iSupport, idVar, iComp, idDof);
596
597 // petsc ordering
598
1/2
✓ Branch 2 taken 54756 times.
✗ Branch 3 not taken.
 54756 AOApplicationToPetsc(m_pLinFluid->aoOld(iold), 1, &idDof);
599
600 // get values
601
2/4
✓ Branch 2 taken 54756 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 54756 times.
✗ Branch 6 not taken.
 54756 fluidSol.getValues( 1, &idDof, &m_fluidDofSol[iUnknown].val(iComp, iSupport));
602 }
603 }
604 }
605 }
606 10196 }
607
608
609 template<class SolidPrb>
610 2056 void LinearProblemFSINitscheXFEMSolid<SolidPrb>::computeIntegrationPoints(std::vector<Point*>& ptElem)
611 {
612 // compute the integration of the substructure in the reference element of the structure element.
613 2056 this->m_feDisp->update(0, m_subItfPts);
614 2056 this->m_feDisp->computeCurrentQuadraturePoint();
615
616 2056 TypeShape fluidCellShape = (TypeShape) m_feFluid->geoEle().shape().typeShape();
617
618 // mapping from the current fluid element to the fluid reference element
619
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2056 times.
 2056 if(fluidCellShape == Quadrilateral) {
620 ✗ for(int ig=0; ig<this->m_feDisp->numQuadraturePoint(); ++ig) {
621 ✗ for(int icomp=0; icomp<m_feFluid->numCoor(); ++icomp) {
622 ✗ m_intPoints[ig][icomp] = 2*(this->m_feDisp->currentQuadPoint[ig][icomp] - (*ptElem[0])[icomp])/((*ptElem[2])[icomp] - (*ptElem[0])[icomp]) - 1.;
623 }
624 }
625
1/2
✓ Branch 0 taken 2056 times.
✗ Branch 1 not taken.
 2056 } else if(fluidCellShape == Triangle) {
626 2056 double det = (ptElem[2]->y() - ptElem[0]->y())*(ptElem[1]->x() - ptElem[0]->x())
627 2056 - (ptElem[2]->x() - ptElem[0]->x())*(ptElem[1]->y() - ptElem[0]->y());
628
629
2/2
✓ Branch 1 taken 4112 times.
✓ Branch 2 taken 2056 times.
 6168 for(int ig=0; ig<this->m_feDisp->numQuadraturePoint(); ++ig) {
630 4112 m_intPoints[ig][0] = (ptElem[2]->y() - ptElem[0]->y())*(this->m_feDisp->currentQuadPoint[ig].x() - ptElem[0]->x())/det + (ptElem[0]->x() - ptElem[2]->x()) * (this->m_feDisp->currentQuadPoint[ig].y() - ptElem[0]->y())/det;
631
632 4112 m_intPoints[ig][1] = (ptElem[0]->y() - ptElem[1]->y())*(this->m_feDisp->currentQuadPoint[ig].x() - ptElem[0]->x())/det + (ptElem[1]->x() - ptElem[0]->x())*(this->m_feDisp->currentQuadPoint[ig].y() - ptElem[0]->y())/det;
633 }
634
635 // m_scalDer[0] = (ptElem[2]->y() - ptElem[1]->y())/det;
636 // m_scalDer[1] = (ptElem[1]->x() - ptElem[2]->x())/det;
637 } else {
638 ✗ FEL_ERROR("This shape type is not implemented yet");
639 }
640 2056 }
641
642 template<class SolidPrb>
643 6106 void LinearProblemFSINitscheXFEMSolid<SolidPrb>::computeFluidVelOnSubStructure(const std::vector<double>& strucNormal)
644 {
645 // Compute the velocity at the integration points of the sub structure
646 // These points are computed in the function computeIntegrationPoints
647
2/2
✓ Branch 1 taken 3053 times.
✓ Branch 2 taken 3053 times.
 6106 if(m_numSolidComponent < m_feFluid->numCoor()) {
648
2/2
✓ Branch 1 taken 3053 times.
✓ Branch 2 taken 3053 times.
 6106 for(int icomp=0; icomp<this->m_feDisp->numRefCoor(); ++icomp) {
649
2/2
✓ Branch 1 taken 6106 times.
✓ Branch 2 taken 3053 times.
 9159 for(int ig=0; ig<this->m_feDisp->numQuadraturePoint(); ++ig) {
650 6106 m_elemNitscheXFEMSolidRHS.val(icomp, ig) = 0;
651
2/2
✓ Branch 1 taken 18318 times.
✓ Branch 2 taken 6106 times.
 24424 for(int idof=0; idof<m_feFluid->numDof(); ++idof) {
652
2/2
✓ Branch 1 taken 36636 times.
✓ Branch 2 taken 18318 times.
 54954 for(int jcomp=0; jcomp<m_feFluid->numCoor(); ++jcomp) {
653 36636 m_elemNitscheXFEMSolidRHS.val(icomp, ig) += m_fluidDofSol[0].val(jcomp, idof) * m_feFluid->refEle().basisFunction().phi(idof, m_intPoints[ig]) * strucNormal[jcomp];
654 }
655 }
656 }
657 }
658 } else {
659
2/2
✓ Branch 0 taken 6106 times.
✓ Branch 1 taken 3053 times.
 9159 for(int icomp=0; icomp<m_numSolidComponent; ++icomp) {
660
2/2
✓ Branch 1 taken 12212 times.
✓ Branch 2 taken 6106 times.
 18318 for(int ig=0; ig<this->m_feDisp->numQuadraturePoint(); ++ig) {
661 12212 m_elemNitscheXFEMSolidRHS.val(icomp, ig) = 0;
662
2/2
✓ Branch 1 taken 36636 times.
✓ Branch 2 taken 12212 times.
 48848 for(int idof=0; idof<m_feFluid->numDof(); ++idof) {
663 36636 m_elemNitscheXFEMSolidRHS.val(icomp, ig) += m_fluidDofSol[0].val(icomp, idof) * m_feFluid->refEle().basisFunction().phi(idof, m_intPoints[ig]);
664 }
665 }
666 }
667 }
668 6106 }
669
670 template<class SolidPrb>
671 4090 void LinearProblemFSINitscheXFEMSolid<SolidPrb>::computeFluidNormalDerivativeJump(double sideCoeff, const std::vector<double>& strucNormal, const std::vector<double>& strucNormalInit)
672 {
673 // compute the tensor at integration points of the sub structure
674 // These points are computed in the function computeIntegrationPoints
675 4090 double tmpval1 = 0;
676 4090 double tmpval2 = 0;
677
678
2/2
✓ Branch 1 taken 2045 times.
✓ Branch 2 taken 2045 times.
 4090 if(m_numSolidComponent < m_feFluid->numCoor()) {
679
2/2
✓ Branch 1 taken 2045 times.
✓ Branch 2 taken 2045 times.
 4090 for(int icomp=0; icomp<this->m_feDisp->numRefCoor(); ++icomp) {
680
2/2
✓ Branch 1 taken 4090 times.
✓ Branch 2 taken 2045 times.
 6135 for(int ig=0; ig<this->m_feDisp->numQuadraturePoint(); ++ig) {
681 4090 m_elemNitscheXFEMSolidRHS.val(icomp, ig) = 0;
682
683 // velocity
684
2/2
✓ Branch 1 taken 8180 times.
✓ Branch 2 taken 4090 times.
 12270 for(int kcomp=0; kcomp<m_feFluid->numCoor(); ++kcomp) {
685 8180 tmpval1 = 0;
686
2/2
✓ Branch 1 taken 24540 times.
✓ Branch 2 taken 8180 times.
 32720 for(int idof=0; idof<m_feFluid->numDof(); ++idof) {
687 24540 tmpval2 = 0;
688
2/2
✓ Branch 1 taken 49080 times.
✓ Branch 2 taken 24540 times.
 73620 for(int jcomp=0; jcomp<m_feFluid->numCoor(); ++jcomp) {
689 49080 tmpval2 += m_fluidDofSol[0].val(jcomp, idof) * strucNormal[jcomp];
690 }
691 // tmpval1 += m_scalDer[kcomp] * m_feFluid->refEle().basisFunction().dPhi(idof, kcomp, m_intPoints[ig]) * tmpval2;
692 24540 tmpval1 += m_feFluid->dPhi[ig](kcomp, idof) * tmpval2;
693 }
694 8180 m_elemNitscheXFEMSolidRHS.val(icomp, ig) += FelisceParam::instance().viscosity * tmpval1 * sideCoeff * strucNormalInit[kcomp];
695 }
696
697 // pressure
698 4090 tmpval1 = 0.;
699
2/2
✓ Branch 1 taken 8180 times.
✓ Branch 2 taken 4090 times.
 12270 for(int kcomp=0; kcomp<m_feFluid->numCoor(); ++kcomp)
700 8180 tmpval1 += strucNormal[kcomp] * strucNormalInit[kcomp];
701
702
2/2
✓ Branch 1 taken 12270 times.
✓ Branch 2 taken 4090 times.
 16360 for(int idof=0; idof<m_feFluid->numDof(); ++idof) {
703 12270 m_elemNitscheXFEMSolidRHS.val(icomp, ig) -= sideCoeff * m_fluidDofSol[1].val(0, idof) * m_fePres->refEle().basisFunction().phi(idof, m_intPoints[ig]) * tmpval1;
704 }
705 }
706 }
707
708 // if we use the symmetric formulation for the stress
709
1/2
✓ Branch 1 taken 2045 times.
✗ Branch 2 not taken.
 2045 if(FelisceParam::instance().useSymmetricStress) {
710
2/2
✓ Branch 1 taken 2045 times.
✓ Branch 2 taken 2045 times.
 4090 for(int icomp=0; icomp<this->m_feDisp->numRefCoor(); ++icomp) {
711
2/2
✓ Branch 1 taken 4090 times.
✓ Branch 2 taken 2045 times.
 6135 for(int ig=0; ig<this->m_feDisp->numQuadraturePoint(); ++ig) {
712
2/2
✓ Branch 1 taken 8180 times.
✓ Branch 2 taken 4090 times.
 12270 for(int jcomp=0; jcomp<m_feFluid->numCoor(); ++jcomp) {
713 8180 tmpval1 = 0;
714
2/2
✓ Branch 1 taken 24540 times.
✓ Branch 2 taken 8180 times.
 32720 for(int idof=0; idof<m_feFluid->numDof(); ++idof) {
715 24540 tmpval2 = 0;
716
2/2
✓ Branch 1 taken 49080 times.
✓ Branch 2 taken 24540 times.
 73620 for(int kcomp=0; kcomp<m_feFluid->numCoor(); ++kcomp) {
717 49080 tmpval2 += m_fluidDofSol[0].val(kcomp, idof) * sideCoeff * strucNormal[kcomp];
718 }
719 // tmpval1 += m_scalDer[jcomp] * m_feFluid->refEle().basisFunction().dPhi(idof, jcomp, m_intPoints[ig]) * tmpval2;
720 24540 tmpval1 += m_feFluid->dPhi[ig](jcomp, idof) * tmpval2;
721 }
722 8180 m_elemNitscheXFEMSolidRHS.val(icomp, ig) += FelisceParam::instance().viscosity * tmpval1 * strucNormalInit[jcomp];
723 }
724 }
725 }
726 }
727 } else {
728
2/2
✓ Branch 0 taken 4090 times.
✓ Branch 1 taken 2045 times.
 6135 for(int icomp=0; icomp<m_numSolidComponent; ++icomp) {
729
2/2
✓ Branch 1 taken 8180 times.
✓ Branch 2 taken 4090 times.
 12270 for(int ig=0; ig<this->m_feDisp->numQuadraturePoint(); ++ig) {
730 8180 m_elemNitscheXFEMSolidRHS.val(icomp, ig) = 0;
731
732 // velocity
733
2/2
✓ Branch 1 taken 16360 times.
✓ Branch 2 taken 8180 times.
 24540 for(int jcomp=0; jcomp<m_feFluid->numCoor(); ++jcomp) {
734 16360 tmpval1 = 0;
735
2/2
✓ Branch 1 taken 49080 times.
✓ Branch 2 taken 16360 times.
 65440 for(int idof=0; idof<m_feFluid->numDof(); ++idof) {
736 49080 tmpval1 += m_fluidDofSol[0].val(icomp, idof) * m_feFluid->dPhi[ig](jcomp, idof);
737 }
738
739 16360 m_elemNitscheXFEMSolidRHS.val(icomp, ig) += FelisceParam::instance().viscosity * tmpval1 * sideCoeff * strucNormal[jcomp];
740 }
741
742 // pressure
743
2/2
✓ Branch 1 taken 24540 times.
✓ Branch 2 taken 8180 times.
 32720 for(int idof=0; idof<m_feFluid->numDof(); ++idof) {
744 24540 m_elemNitscheXFEMSolidRHS.val(icomp, ig) -= sideCoeff * m_fluidDofSol[1].val(0, idof) * m_fePres->refEle().basisFunction().phi(idof, m_intPoints[ig]) * strucNormal[icomp];
745 }
746 }
747 }
748
749 // if we use the symmetric formulation for the stress
750
1/2
✓ Branch 1 taken 2045 times.
✗ Branch 2 not taken.
 2045 if(FelisceParam::instance().useSymmetricStress) {
751
2/2
✓ Branch 0 taken 4090 times.
✓ Branch 1 taken 2045 times.
 6135 for(int icomp=0; icomp<m_numSolidComponent; ++icomp) {
752
2/2
✓ Branch 1 taken 8180 times.
✓ Branch 2 taken 4090 times.
 12270 for(int ig=0; ig<this->m_feDisp->numQuadraturePoint(); ++ig) {
753
2/2
✓ Branch 1 taken 24540 times.
✓ Branch 2 taken 8180 times.
 32720 for(int idof=0; idof<m_feFluid->numDof(); ++idof) {
754 24540 tmpval1 = 0;
755
2/2
✓ Branch 1 taken 49080 times.
✓ Branch 2 taken 24540 times.
 73620 for(int jcomp=0; jcomp<m_feFluid->numCoor(); ++jcomp) {
756 49080 tmpval1 += m_fluidDofSol[0].val(jcomp, idof) * strucNormal[jcomp];
757 }
758
759 24540 m_elemNitscheXFEMSolidRHS.val(icomp, ig) += FelisceParam::instance().viscosity * tmpval1 * sideCoeff * m_feFluid->dPhi[ig](icomp, idof);
760 }
761 }
762 }
763 }
764 }
765 4090 }
766
767 /*!
768 * Gather the last time step displacement into m_seqLastTimeSol.
769 * If m_seqLastTimeSol is null, it will be created and allocated during the gathering
770 * according to the parallel shape of parStrucLastSol
771 */
772 template<class SolidPrb>
773 20 void LinearProblemFSINitscheXFEMSolid<SolidPrb>::gatherSeqLastTimeSol(PetscVector& parStrucLastSol)
774 {
775 20 this->gatherVector(parStrucLastSol, m_seqLastTimeSol);
776 20 m_isSeqLastTimeSolAllocated = true;
777 20 }
778
779 /*!
780 * Gather the last time step displacement-velocity into m_seqLastTimeVelSol.
781 * If m_seqLastTimeVelSol is null, it will be created and allocated during the gathering
782 * according to the parallel shape of parStrucLastVelSol
783 */
784 template<class SolidPrb>
785 void LinearProblemFSINitscheXFEMSolid<SolidPrb>::gatherSeqLastTimeVelSol(PetscVector& parStrucLastVelSol)
786 {
787 this->gatherVector(parStrucLastVelSol, m_seqLastTimeVelSol);
788 m_isSeqLastTimeVelSolAllocated = true;
789 }
790
791 /*!
792 * Gather the extrapolation Robin-Neumann of the displacement into m_seqRNStrucExtrapol.
793 * If m_seqRNStrucExtrapol is null, it will be created and allocated during the gathering
794 * according to the parallel shape of parStrucExtrapol
795 */
796 template<class SolidPrb>
797 ✗ void LinearProblemFSINitscheXFEMSolid<SolidPrb>::gatherRNStructExtrapol(PetscVector& parRNStrucExtrapol) {
798 ✗ this->gatherVector(parRNStrucExtrapol, m_seqRNStrucExtrapol);
799 ✗ m_isSeqRNStrucExtrapolAllocated = true;
800 }
801
802 /*!
803 * Gather the extrapolation of the fluid velocity for the stabilized explicit splitting scheme
804 * into m_seqStabExplicitFluidExtrapol.
805 * If m_seqStabExplicitFluidExtrapol is null, it will be created and allocated during the gathering
806 * according to the parallel shape of parFluidExtrapol
807 */
808 template<class SolidPrb>
809 38 void LinearProblemFSINitscheXFEMSolid<SolidPrb>::gatherStabExplicitFluidExtrapol(std::vector<PetscVector>& parFluidExtrapol)
810 {
811
2/2
✓ Branch 1 taken 76 times.
✓ Branch 2 taken 38 times.
 114 for(std::size_t i=0; i<parFluidExtrapol.size(); ++i)
812 76 this->gatherVector(parFluidExtrapol[i], m_seqStabExplicitFluidExtrapol[i]);
813
814 38 m_isSeqStabExplicitFluidExtrapolAllocated = true;
815 38 }
816
817 template<class SolidPrb>
818 20 void LinearProblemFSINitscheXFEMSolid<SolidPrb>::gatherStabExplicitFluidExtrapol(PetscVector& parFluidExtrapol)
819 {
820 20 this->gatherVector(parFluidExtrapol, m_seqStabExplicitFluidExtrapol[0]);
821 20 m_isSeqStabExplicitFluidExtrapolAllocated = true;
822 20 }
823
824 template<class SolidPrb>
825 20 void LinearProblemFSINitscheXFEMSolid<SolidPrb>::setStabExplicitFluidExtrapol(PetscVector& seqFluidExtrapol)
826 {
827 20 m_seqStabExplicitFluidExtrapol[0].copyFrom(seqFluidExtrapol);
828 20 }
829
830 template<class SolidPrb>
831 42 void LinearProblemFSINitscheXFEMSolid<SolidPrb>::initStabExplicitFluidExtrapol(felInt order)
832 {
833
2/2
✓ Branch 0 taken 27 times.
✓ Branch 1 taken 15 times.
 42 if(m_isSeqStabExplicitFluidExtrapolAllocated) {
834
2/2
✓ Branch 1 taken 47 times.
✓ Branch 2 taken 27 times.
 74 for(std::size_t i=0; i<m_seqStabExplicitFluidExtrapol.size(); ++i)
835 47 m_seqStabExplicitFluidExtrapol[i].destroy();
836
837 27 m_isSeqStabExplicitFluidExtrapolAllocated = false;
838 }
839
840 42 m_seqStabExplicitFluidExtrapol.resize(order);
841
2/2
✓ Branch 1 taken 64 times.
✓ Branch 2 taken 42 times.
 106 for(std::size_t i=0; i<m_seqStabExplicitFluidExtrapol.size(); ++i)
842
1/2
✓ Branch 3 taken 64 times.
✗ Branch 4 not taken.
 64 m_seqStabExplicitFluidExtrapol[i] = PetscVector::null();
843 42 }
844
845 template<class SolidPrb>
846 11 void LinearProblemFSINitscheXFEMSolid<SolidPrb>::deleteDynamicData()
847 {
848
1/2
✓ Branch 0 taken 11 times.
✗ Branch 1 not taken.
 11 if(m_isSeqStabExplicitFluidExtrapolAllocated) {
849
2/2
✓ Branch 1 taken 11 times.
✓ Branch 2 taken 11 times.
 22 for(std::size_t i=0; i<m_seqStabExplicitFluidExtrapol.size(); ++i) {
850 11 m_seqStabExplicitFluidExtrapol[i].destroy();
851
1/2
✓ Branch 3 taken 11 times.
✗ Branch 4 not taken.
 11 m_seqStabExplicitFluidExtrapol[i] = PetscVector::null();
852 }
853 11 m_isSeqStabExplicitFluidExtrapolAllocated = false;
854 }
855 11 }
856
857 template<class SolidPrb>
858 ✗ double LinearProblemFSINitscheXFEMSolid<SolidPrb>::computeL2Norm(PetscVector& inputVec)
859 {
860 ✗ if(m_seqVecL2Norm.isNull())
861 ✗ m_seqVecL2Norm.duplicateFrom(this->sequentialSolution());
862
863 ✗ this->gatherVector(inputVec, m_seqVecL2Norm);
864
865 GeometricMeshRegion::ElementType eltType;
866 ✗ felInt numPointPerElt = 0;
867 ✗ felInt numEltPerLabel = 0;
868
869 ✗ std::vector<Point*> elemPoint;
870 ✗ std::vector<felInt> elemIdPoint;
871 ✗ std::vector <felInt> vectorIdSupport;
872
873 // use to identify global id of dof.
874 felInt idDof;
875 ✗ felInt cpt = 0;
876
877 // Id link element to its supportDof.
878 felInt ielSupportDof;
879
880 // return value
881 ✗ double l2Norm = 0.;
882 ✗ double tmpL2Norm = 0.;
883
884 // initialize id global of elements.
885 felInt numElement[GeometricMeshRegion::m_numTypesOfElement ];
886 ✗ for (int ityp=0; ityp<GeometricMeshRegion::m_numTypesOfElement; ityp++ ) {
887 ✗ eltType = (GeometricMeshRegion::ElementType)ityp;
888 ✗ numElement[eltType] = 0;
889 }
890
891 CurvilinearFiniteElement *fePtr;
892 ✗ std::vector<double> dofValue;
893 ✗ std::vector<felInt> idGlobalDof;
894
895 ✗ const std::vector<GeometricMeshRegion::ElementType>& bagElementTypeDomain = this->meshLocal()->bagElementTypeDomain();
896 ✗ for (std::size_t i = 0; i < bagElementTypeDomain.size(); ++i) {
897 ✗ eltType = bagElementTypeDomain[i];
898 ✗ numPointPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
899 ✗ elemPoint.resize(numPointPerElt, 0);
900 ✗ elemIdPoint.resize(numPointPerElt, 0);
901
902 // define the curvilinear finite element
903 ✗ this->defineCurvilinearFiniteElement(eltType);
904 ✗ this->initElementArrayBD();
905 ✗ fePtr = this->m_listCurvilinearFiniteElement[m_iDisplacement];
906
907 ✗ dofValue.resize(fePtr->numDof()*m_numSolidComponent, 0.);
908 ✗ idGlobalDof.resize(fePtr->numDof()*m_numSolidComponent, 0);
909
910 ✗ for(auto itRef = this->meshLocal()->intRefToBegEndMaps[eltType].begin(); itRef != this->meshLocal()->intRefToBegEndMaps[eltType].end(); itRef++) {
911 ✗ numEltPerLabel = itRef->second.second;
912
913 ✗ for (felInt iel = 0; iel < numEltPerLabel; iel++) {
914 ✗ this->setElemPoint(eltType, numElement[eltType], elemPoint, elemIdPoint, vectorIdSupport);
915 ✗ fePtr->updateMeas(0, elemPoint);
916
917 // loop over all the support elements
918 ✗ for (std::size_t it = 0; it < vectorIdSupport.size(); it++) {
919 // get the id of the support
920 ✗ ielSupportDof = vectorIdSupport[it];
921
922 ✗ cpt = 0;
923 ✗ for(felInt iDof=0; iDof<fePtr->numDof(); iDof++) {
924 ✗ for(felInt iComp = 0; iComp<m_numSolidComponent; iComp++) {
925 ✗ this->dof().loc2glob(ielSupportDof, iDof, m_iDisplacement, iComp, idDof);
926 ✗ idGlobalDof[cpt] = idDof;
927 ✗ ++cpt;
928 }
929 }
930
931 // mapping to obtain new global numbering. (I/O array: idGlobalDof).
932 ✗ AOApplicationToPetsc(this->m_ao, fePtr->numDof()*m_numSolidComponent, idGlobalDof.data());
933
934 // gets values of associate dofs.
935 ✗ m_seqVecL2Norm.getValues(fePtr->numDof()*m_numSolidComponent, idGlobalDof.data(), dofValue.data());
936
937 ✗ for(felInt ig=0; ig<fePtr->numQuadraturePoint(); ++ig) {
938 ✗ double tmpval = 0.;
939 ✗ for(felInt iComp=0; iComp<m_numSolidComponent; ++iComp) {
940 ✗ double uicomp = 0.;
941 ✗ for(felInt idof = 0; idof<fePtr->numDof(); ++idof) {
942 ✗ uicomp += fePtr->phi[ig](idof) * dofValue[idof*m_numSolidComponent + iComp] ;
943 }
944 ✗ tmpval += uicomp * uicomp;
945 }
946 ✗ tmpL2Norm += tmpval * fePtr->weightMeas(ig);
947 }
948 }
949 ✗ numElement[eltType]++;
950 }
951 }
952 }
953
954 ✗ MPI_Allreduce(&tmpL2Norm, &l2Norm, 1, MPI_DOUBLE, MPI_SUM, MpiInfo::petscComm());
955 ✗ return l2Norm;
956 }
957 }
958