GCC Code Coverage Report

Directory: ./
File: Solver/linearProblemStokesContinuation.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 283 288 98.3%
Branches: 223 370 60.3%
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: J. Foulon & J-F. Gerbeau & V. Martin
13 //
14
15 // System includes
16 #include <stack>
17
18 // External includes
19
20 // Project includes
21 #include "Solver/linearProblemStokesContinuation.hpp"
22 #include "Core/felisceTransient.hpp"
23 #include "FiniteElement/elementVector.hpp"
24 #include "FiniteElement/elementMatrix.hpp"
25 #include "InputOutput/io.hpp"
26 #include "DegreeOfFreedom/boundaryCondition.hpp"
27
28 /*!
29 \file linearProblemStokesContinuation.cpp
30 \date 17/03/2022
31 \brief Continuation method for Stokes' equation
32 */
33
34 namespace felisce {
35 4 LinearProblemStokesContinuation::LinearProblemStokesContinuation():
36
5/10
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
✓ Branch 10 taken 4 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 4 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 4 times.
✗ Branch 17 not taken.
 4 LinearProblem("Poisson continuation equation")
37 {
38 4 }
39
40 20 void LinearProblemStokesContinuation::readData(IO& io,double iteration) {
41 // Read velocity variable from the .case file
42 20 std::vector<double> potData;
43
1/2
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
 20 potData.resize(m_mesh[m_currentMesh]->numPoints()*3, 0.);
44
1/2
✓ Branch 3 taken 20 times.
✗ Branch 4 not taken.
 20 io.readVariable(0,iteration, potData.data(), potData.size());
45
46 20 std::vector<double> potDataVelocity;
47 20 int dim = this->dimension();
48
1/2
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
 20 potDataVelocity.resize(m_mesh[m_currentMesh]->numPoints()*dim, 0.);
49
2/2
✓ Branch 3 taken 2840 times.
✓ Branch 4 taken 20 times.
 2860 for (int i=0; i < m_mesh[m_currentMesh]->numPoints(); ++i) {
50
2/2
✓ Branch 0 taken 5680 times.
✓ Branch 1 taken 2840 times.
 8520 for(int icomp=0; icomp<dim; icomp++) {
51 5680 potDataVelocity[dim*i+icomp]=potData[3*i+icomp];
52 }
53 }
54
55
1/2
✓ Branch 2 taken 20 times.
✗ Branch 3 not taken.
 20 m_potData.duplicateFrom(this->sequentialSolution());
56
57
1/2
✓ Branch 3 taken 20 times.
✗ Branch 4 not taken.
 20 std::vector<PetscInt> ao_potData( potDataVelocity.size());
58
2/2
✓ Branch 1 taken 5680 times.
✓ Branch 2 taken 20 times.
 5700 for (size_t i=0; i < potDataVelocity.size(); ++i) {
59 5680 ao_potData[i]=i;
60 }
61
62
1/2
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
 20 AOApplicationToPetsc(this->ao(),potDataVelocity.size(),ao_potData.data());
63
2/4
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 20 times.
✗ Branch 8 not taken.
 20 VecSetValues(m_potData.toPetsc(),potDataVelocity.size(),ao_potData.data(), potDataVelocity.data(), INSERT_VALUES);
64 20 }
65
66 4 void LinearProblemStokesContinuation::initialize(std::vector<GeometricMeshRegion::Pointer>& mesh, FelisceTransient::Pointer fstransient, MPI_Comm& comm, bool doUseSNES) {
67
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 LinearProblem::initialize(mesh,comm, doUseSNES);
68 4 this->m_fstransient = fstransient;
69 4 std::vector<PhysicalVariable> listVariable;
70 4 std::vector<std::size_t> listNumComp;
71
72 // primal variables
73
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 listVariable.push_back(velocity);
74
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 listNumComp.push_back(this->dimension());
75
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 listVariable.push_back(pressure);
76
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 listNumComp.push_back(1);
77
78 // dual variables with dummy names
79
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 listVariable.push_back(velocityDiv);
80
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 listNumComp.push_back(this->dimension());
81
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 listVariable.push_back(potThorax);
82
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 listNumComp.push_back(1);
83
84 //define unknown of the linear system.
85
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 m_listUnknown.push_back(velocity);
86
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 m_listUnknown.push_back(pressure);
87
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 m_listUnknown.push_back(velocityDiv);
88
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 m_listUnknown.push_back(potThorax);
89
90
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 userAddOtherUnknowns(listVariable,listNumComp);
91
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 userAddOtherVariables(listVariable,listNumComp);
92
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 definePhysicalVariable(listVariable,listNumComp);
93
94
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 m_iVelocity = m_listVariable.getVariableIdList(velocity);
95
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 m_iPressure = m_listVariable.getVariableIdList(pressure);
96
97
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 m_viscosity = FelisceParam::instance().viscosity;
98
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 m_useSymmetricStress = FelisceParam::instance().useSymmetricStress;
99 4 }
100
101 40 void LinearProblemStokesContinuation::initPerElementType(ElementType eltType, FlagMatrixRHS flagMatrixRHS) {
102 IGNORE_UNUSED_ARGUMENT(eltType);
103 IGNORE_UNUSED_ARGUMENT(flagMatrixRHS);
104 40 m_feVel = m_listCurrentFiniteElement[m_iVelocity];
105 40 m_fePres = m_listCurrentFiniteElement[m_iPressure];
106
107 // define elemField to store the given velocity measurements and the source terms
108 40 m_dataTerm.initialize(DOF_FIELD, *m_feVel, this->dimension());
109 40 m_sourceTerm.initialize(DOF_FIELD, *m_feVel, this->dimension());
110 40 }
111
112
113 1210 void LinearProblemStokesContinuation::computeElementArray(const std::vector<Point*>& elemPoint, const std::vector<felInt>& /*elemIdPoint*/, felInt& iel, FlagMatrixRHS flagMatrixRHS) {
114 1210 m_feVel->updateFirstDeriv(0, elemPoint);
115 1210 m_fePres->updateFirstDeriv(0, elemPoint);
116
117 1210 double hK = m_feVel->diameter();
118
119 // stabilization parameters
120 1210 double gamma_M = FelisceParam::instance().referenceValue2; // data fidelity coefficient
121 1210 double gamma_p = hK*hK*FelisceParam::instance().referenceValue3; // Brezzi-Pitkaranta coefficient
122 1210 double gamma_p_star = FelisceParam::instance().referenceValue4; // dual pressure coefficient
123 1210 double gamma_div = 0.1;
124 1210 double gamma_u_star = 0.1;
125
126 // time derivative coefficient
127 1210 double rhof = FelisceParam::instance(this->instanceIndex()).density;
128 1210 double coef_time = rhof/m_fstransient->timeStep;
129
130 // defining the block matrix of the system
131
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1210 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 1210 if (flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs || flagMatrixRHS == FlagMatrixRHS::only_matrix) {
132 // stabilizing div term acting on u, \int \nabla \cdot u \nabla \cdot v
133 1210 m_elementMat[0]->div_phi_j_div_phi_i(gamma_div, *m_feVel, 0, 0, this->dimension()); // block (0,0)
134 // data term, int_\Omega u \cdot v; to this block we will also add the cip stabilization
135 1210 m_elementMat[0]->phi_i_phi_j(gamma_M, *m_feVel, 0, 0, this->dimension()); // block (0,0)
136
137 // stabilizing term acting on p, int_\Omega grad p : grad q
138 1210 m_elementMat[0]->grad_phi_i_grad_phi_j(gamma_p, *m_fePres, this->dimension(), this->dimension(), 1); // block (1,1)
139
140 // time discretization for z, int_\Omega v \cdot z
141 1210 m_elementMat[0]->phi_i_phi_j(coef_time, *m_feVel, 0, 1+this->dimension(), this->dimension()); // block (0,2)
142
143 // Stokes blocks for z and y, int_\Omega \epsilon v : \epsilon z or int_\Omega \nabla v : \nabla z
144
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1210 times.
 1210 if (m_useSymmetricStress)
145 ✗ m_elementMat[0]->eps_phi_i_eps_phi_j(2*m_viscosity, *m_feVel, 0, 1+this->dimension(), this->dimension()); // block (0,2)
146 else
147 1210 m_elementMat[0]->grad_phi_i_grad_phi_j(m_viscosity, *m_feVel, 0, 1+this->dimension(), this->dimension()); // block (0,2)
148 // \int y \nabla \cdot v and -\int q \nabla \cdot z
149 1210 m_elementMat[0]->psi_j_div_phi_i(1., *m_feVel, *m_fePres, 0, 1+2*this->dimension()); // block (0,3)
150 1210 m_elementMat[0]->psi_i_div_phi_j(-1., *m_feVel, *m_fePres, this->dimension(), 1+this->dimension()); // block (1,2)
151
152 // time discretization for u, int_\Omega u \cdot w
153 1210 m_elementMat[0]->phi_i_phi_j(coef_time, *m_feVel, 1+this->dimension(), 0, this->dimension()); // block (2,0)
154
155 // Stokes blocks for u and p, int_\Omega \epsilon u : \epsilon v or int_\Omega \nabla u : \nabla v
156
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1210 times.
 1210 if (m_useSymmetricStress)
157 ✗ m_elementMat[0]->eps_phi_i_eps_phi_j(2*m_viscosity, *m_feVel, 1+this->dimension(), 0, this->dimension()); // block (2,0)
158 else
159 1210 m_elementMat[0]->grad_phi_i_grad_phi_j(m_viscosity, *m_feVel, 1+this->dimension(), 0, this->dimension()); // block (2,0)
160 // -int p \nabla \cdot w and \int x \nabla \cdot u
161 1210 m_elementMat[0]->psi_j_div_phi_i(-1., *m_feVel, *m_fePres, 1+this->dimension(), this->dimension()); // block (2,1)
162 1210 m_elementMat[0]->psi_i_div_phi_j(1., *m_feVel, *m_fePres, 1+2*this->dimension(), 0); // block (3,0)
163
164 // stabilizing term acting on z, -int_\Omega grad z : grad w
165 1210 m_elementMat[0]->grad_phi_i_grad_phi_j(-gamma_u_star, *m_feVel, 1+this->dimension(), 1+this->dimension(), this->dimension()); // block (2,2)
166 // stabilizing term acting on y, -int_\Omega y * x
167 1210 m_elementMat[0]->phi_i_phi_j(-gamma_p_star, *m_fePres, 1+2*this->dimension(), 1+2*this->dimension(), 1); // block (3,3)
168 }
169 // defining the rhs of the system
170
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1210 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 1210 if (flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs || flagMatrixRHS == FlagMatrixRHS::only_rhs) {
171 //rhs data term, using the measurements
172 1210 m_dataTerm.setValue(m_potData, *m_feVel, iel, m_iVelocity, m_ao, dof());
173
174 // gamma_M int_\Omega u_M \cdot v
175 1210 m_elementVector[0]->source(gamma_M, *m_feVel, m_dataTerm, 0, this->dimension());
176
177 //rhs f term, zero
178
179 // rhs source term, the solution at the previous time step
180 1210 m_sourceTerm.setValue(this->sequentialSolution(), *m_feVel, iel, m_iVelocity, m_ao, dof());
181
182 // rhof/tau int_\Omega u^{n-1} \cdot w
183 1210 m_elementVector[0]->source(coef_time, *m_feVel, m_sourceTerm, 1+this->dimension(), this->dimension());
184 }
185 1210 }
186
187 4 void LinearProblemStokesContinuation::userChangePattern(int numProc, int rankProc) {
188 IGNORE_UNUSED_ARGUMENT(numProc);
189 IGNORE_UNUSED_ARGUMENT(rankProc);
190 // compute initial (uniform) repartition of dof on processes
191
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 felInt numDofByProc = m_numDof/MpiInfo::numProc();
192
2/4
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
 4 std::vector<felInt> numDofPerProc(MpiInfo::numProc());
193
3/4
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 16 times.
✓ Branch 4 taken 4 times.
 20 for(felInt i=0; i<MpiInfo::numProc(); ++i) {
194
3/4
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
✓ Branch 4 taken 12 times.
 16 if(i == MpiInfo::numProc() - 1)
195 4 numDofPerProc[i] = m_numDof - i*numDofByProc;
196 else
197 12 numDofPerProc[i] = numDofByProc;
198 }
199
200 4 felInt shiftDof = 0;
201
3/4
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 6 times.
✓ Branch 4 taken 4 times.
 10 for(felInt i=0; i<MpiInfo::rankProc(); ++i)
202 6 shiftDof += numDofPerProc[i];
203
204
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 std::vector<felInt> rankDof(m_numDof, -1);
205
3/4
✓ Branch 1 taken 856 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 852 times.
✓ Branch 5 taken 4 times.
 856 for(felInt i=0; i<numDofPerProc[MpiInfo::rankProc()]; ++i)
206
1/2
✓ Branch 1 taken 852 times.
✗ Branch 2 not taken.
 852 rankDof[i + shiftDof] = MpiInfo::rankProc();
207
208 // build the edges or faces depending on the dimension (for the global mesh)
209 // In this function, "faces" refers to either edges or faces.
210
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 if(dimension() == 2)
211
1/2
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
 4 m_mesh[m_currentMesh]->buildEdges();
212 else
213 ✗ m_mesh[m_currentMesh]->buildFaces();
214
215 // variables
216 GeometricMeshRegion::ElementType eltType, eltTypeOpp;
217 felInt idVar1, idVar2;
218 felInt node1, node2;
219 felInt numFacesPerElement, numEltPerLabel;
220 felInt ielSupport, jelSupport;
221 felInt ielCurrentGeo, ielOppositeGeo;
222 4 std::vector<felInt> vecSupport, vecSupportOpposite;
223
1/2
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 4 std::vector<felInt> numElement(m_mesh[m_currentMesh]->m_numTypesOfElement, 0);
224
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 std::vector< std::set<felInt> > nodesNeighborhood(m_numDof);
225
226
227 // zeroth loop on the unknown of the linear problem
228
2/2
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 4 times.
 20 for (size_t iUnknown1 = 0; iUnknown1 < m_listUnknown.size(); ++iUnknown1) {
229 16 idVar1 = m_listUnknown.idVariable(iUnknown1);
230
231 16 ielCurrentGeo = 0;
232
233 // first loop on element type
234
2/2
✓ Branch 4 taken 16 times.
✓ Branch 5 taken 16 times.
 32 for (size_t i=0; i<m_mesh[m_currentMesh]->bagElementTypeDomain().size(); ++i) {
235 16 eltType = m_mesh[m_currentMesh]->bagElementTypeDomain()[i];
236 16 const GeoElement* geoEle = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second;
237 16 numElement[eltType] = 0;
238 16 numFacesPerElement = geoEle->numBdEle();
239
240 // second loop on region of the mesh.
241
2/2
✓ Branch 10 taken 16 times.
✓ Branch 11 taken 16 times.
 32 for(GeometricMeshRegion::IntRefToBegEndIndex_type::const_iterator itRef =m_mesh[m_currentMesh]->intRefToBegEndMaps[eltType].begin(); itRef !=m_mesh[m_currentMesh]->intRefToBegEndMaps[eltType].end(); itRef++) {
242 16 numEltPerLabel = itRef->second.second;
243
244 // Third loop on element
245
2/2
✓ Branch 0 taken 3872 times.
✓ Branch 1 taken 16 times.
 3888 for (felInt iel = 0; iel < numEltPerLabel; iel++) {
246
1/2
✓ Branch 3 taken 3872 times.
✗ Branch 4 not taken.
 3872 m_supportDofUnknown[iUnknown1].getIdElementSupport(eltType, numElement[eltType], vecSupport);
247
248
2/2
✓ Branch 1 taken 3872 times.
✓ Branch 2 taken 3872 times.
 7744 for(size_t ielSup=0; ielSup<vecSupport.size(); ++ielSup) {
249 3872 ielSupport = vecSupport[ielSup];
250
251 // for all support dof in this support element
252
2/2
✓ Branch 2 taken 11616 times.
✓ Branch 3 taken 3872 times.
 15488 for (felInt iSup = 0; iSup < m_supportDofUnknown[iUnknown1].getNumSupportDof(ielSupport); ++iSup) {
253 // for all component of the unknown
254
2/2
✓ Branch 2 taken 17424 times.
✓ Branch 3 taken 11616 times.
 29040 for (size_t iComp = 0; iComp < m_listVariable[idVar1].numComponent(); iComp++) {
255 // get the global id of the support dof
256
1/2
✓ Branch 2 taken 17424 times.
✗ Branch 3 not taken.
 17424 dof().loc2glob(ielSupport, iSup, idVar1, iComp, node1);
257
258 // if this node is on this process
259
3/4
✓ Branch 2 taken 17424 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 4356 times.
✓ Branch 5 taken 13068 times.
 17424 if(rankDof[node1] == MpiInfo::rankProc()) {
260 // loop over the boundaries of the element
261
2/2
✓ Branch 0 taken 13068 times.
✓ Branch 1 taken 4356 times.
 17424 for(felInt iface=0; iface < numFacesPerElement; ++iface) {
262 // check if this face is an inner face or a boundary
263 13068 ielOppositeGeo = ielCurrentGeo;
264 13068 eltTypeOpp = eltType;
265
1/2
✓ Branch 3 taken 13068 times.
✗ Branch 4 not taken.
 13068 bool isAnInnerFace =m_mesh[m_currentMesh]->getAdjElement(eltTypeOpp, ielOppositeGeo, iface);
266
267
2/2
✓ Branch 0 taken 12348 times.
✓ Branch 1 taken 720 times.
 13068 if(isAnInnerFace) {
268 // for all unknown
269
2/2
✓ Branch 1 taken 49392 times.
✓ Branch 2 taken 12348 times.
 61740 for (size_t iUnknown2 = 0; iUnknown2 < m_listUnknown.size(); ++iUnknown2) {
270 49392 idVar2 = m_listUnknown.idVariable(iUnknown2);
271
272 // get the support element of the opposite element
273
1/2
✓ Branch 2 taken 49392 times.
✗ Branch 3 not taken.
 49392 m_supportDofUnknown[iUnknown2].getIdElementSupport(ielOppositeGeo, vecSupportOpposite);
274
275 // for all support element
276
2/2
✓ Branch 1 taken 49392 times.
✓ Branch 2 taken 49392 times.
 98784 for(size_t jelSup=0; jelSup<vecSupportOpposite.size(); ++jelSup) {
277 49392 jelSupport = vecSupportOpposite[jelSup];
278
279 // for all support dof in this support element
280
2/2
✓ Branch 2 taken 148176 times.
✓ Branch 3 taken 49392 times.
 197568 for (felInt jSup = 0; jSup < m_supportDofUnknown[iUnknown2].getNumSupportDof(jelSupport); ++jSup) {
281
282 // for all component of the second unknown
283
2/2
✓ Branch 2 taken 222264 times.
✓ Branch 3 taken 148176 times.
 370440 for (size_t jComp = 0; jComp < m_listVariable[idVar2].numComponent(); jComp++) {
284
285 // If the two current components are connected
286
1/2
✓ Branch 2 taken 222264 times.
✗ Branch 3 not taken.
 222264 felInt iConnect = dof().getNumGlobComp(iUnknown1, iComp);
287
1/2
✓ Branch 2 taken 222264 times.
✗ Branch 3 not taken.
 222264 felInt jConnect = dof().getNumGlobComp(iUnknown2, jComp);
288
2/4
✓ Branch 2 taken 222264 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 222264 times.
✗ Branch 5 not taken.
 222264 if (m_listUnknown.mask()(iConnect, jConnect) > 0) {
289 // get the global id of the second support dof
290
1/2
✓ Branch 2 taken 222264 times.
✗ Branch 3 not taken.
 222264 dof().loc2glob(jelSupport, jSup, idVar2, jComp, node2);
291
292 // remove diagonal term to define pattern and use it in Parmetis.
293 // if the two global ids are different, they are neighboors
294
2/2
✓ Branch 0 taken 214032 times.
✓ Branch 1 taken 8232 times.
 222264 if(node1 != node2)
295
1/2
✓ Branch 2 taken 214032 times.
✗ Branch 3 not taken.
 214032 nodesNeighborhood[node1].insert(node2);
296 }
297 }
298 }
299 }
300 }
301 }
302 }
303 }
304 }
305 }
306 }
307 3872 ++ielCurrentGeo;
308 3872 ++numElement[eltType];
309 }
310 }
311 }
312 }
313
314 // Store the pattern in CSR style
315 4 std::vector<felInt> iCSR, jCSR;
316 4 felInt dofSize = 0;
317 4 felInt cptDof = 0;
318 felInt pos;
319
320
1/2
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 4 iCSR.resize(dof().pattern().numRows() + 1, 0);
321
2/2
✓ Branch 1 taken 3408 times.
✓ Branch 2 taken 4 times.
 3412 for(size_t iNode=0; iNode<nodesNeighborhood.size(); ++iNode)
322 3408 dofSize += nodesNeighborhood[iNode].size();
323
324
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 jCSR.resize(dofSize, 0);
325
2/2
✓ Branch 1 taken 3408 times.
✓ Branch 2 taken 4 times.
 3412 for(size_t iNode=0; iNode<nodesNeighborhood.size(); ++iNode) {
326
3/4
✓ Branch 2 taken 3408 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 852 times.
✓ Branch 5 taken 2556 times.
 3408 if(rankDof[iNode] == MpiInfo::rankProc()) {
327 852 iCSR[cptDof + 1] = iCSR[cptDof] + nodesNeighborhood[iNode].size();
328 852 pos = 0;
329
2/2
✓ Branch 5 taken 55956 times.
✓ Branch 6 taken 852 times.
 56808 for(std::set<felInt>::iterator it=nodesNeighborhood[iNode].begin(); it != nodesNeighborhood[iNode].end(); ++it) {
330 55956 jCSR[iCSR[cptDof] + pos] = *it;
331 55956 ++pos;
332 }
333 852 ++cptDof;
334 }
335 }
336
337 // Now, call merge pattern
338
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 dof().mergeGlobalPattern(iCSR, jCSR);
339 //std::cout<<"After global merge"<<std::endl;
340 //dof().pattern().print(1, std::cout);
341 4 }
342
343 // compared to the function linearProblemNS::assembleFaceOrientedStabilization, here we use different blocks and only the jumps in the velocity
344 20 void LinearProblemStokesContinuation::assembleCIPStabilization() {
345 // definition of variables
346 20 std::pair<bool, GeometricMeshRegion::ElementType> adjElt;
347 GeometricMeshRegion::ElementType eltType; // Type of element
348
349 20 felInt ielCurrentLocalGeo = 0; // local geometric id of the current element
350 20 felInt ielCurrentGlobalGeo = 0; // global geometric id of the current element
351 20 felInt ielOppositeGlobalGeo = 0; // global geometric id of the opposite element
352 felInt numFacesPerType; // number of faces of an element of a given type
353 felInt numPointPerElt; // number of vertices by element
354 felInt ielCurrentGlobalSup; // global support element id of the current element
355 felInt ielOppositeGlobalSup; // global support element id of the opposite element
356
357 20 std::vector<felInt> idOfFacesCurrent; // ids of the current element edges
358 20 std::vector<felInt> idOfFacesOpposite; // ids of the opposite element edges
359 20 std::vector<felInt> currentElemIdPoint; // ids of the vertices of the current element
360 20 std::vector<felInt> oppositeElemIdPoint; // ids of the vertices of the opposite element
361 20 std::vector<Point*> currentElemPoint; // point coordinates of the current element vertices
362 20 std::vector<Point*> oppositeElemPoint; // point coordinates of the opposite element vertices
363
364 20 FlagMatrixRHS flag = FlagMatrixRHS::only_matrix; // flag to only assemble the matrix
365 felInt idFaceToDo;
366 20 bool allDone = false;
367
368
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 ElementField elemFieldAdvFace;
369
370 felInt rank;
371
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
372
373 // bag element type vector
374 20 const std::vector<ElementType>& bagElementTypeDomain = m_meshLocal[m_currentMesh]->bagElementTypeDomain();
375 20 eltType = bagElementTypeDomain[0];
376
377
378 // Finite Element With Bd for the opposite element
379 CurrentFiniteElementWithBd* ptmp;
380 CurrentFiniteElementWithBd* oppositeFEWithBd;
381
382
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 m_velocity = &m_listVariable[m_listVariable.getVariableIdList(velocity)];
383
384 20 const GeoElement *geoEle = m_mesh[m_currentMesh]->eltEnumToFelNameGeoEle[eltType].second;
385 20 felInt typeOfFEVel = m_velocity->finiteElementType();
386
387
1/2
✓ Branch 3 taken 20 times.
✗ Branch 4 not taken.
 20 const RefElement *refEleVel = geoEle->defineFiniteEle(eltType, typeOfFEVel, *m_mesh[m_currentMesh]);
388
2/4
✓ Branch 3 taken 20 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 20 times.
✗ Branch 7 not taken.
 20 oppositeFEWithBd = new CurrentFiniteElementWithBd(*refEleVel, *geoEle, m_velocity->degreeOfExactness(), m_velocity->degreeOfExactness());
389
390 // initializing variables
391 20 numPointPerElt = m_meshLocal[m_currentMesh]->m_numPointsPerElt[eltType];
392 20 numFacesPerType = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second->numBdEle();
393
394 // resize of vectors
395
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 idOfFacesCurrent.resize(numFacesPerType, -1);
396
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 idOfFacesOpposite.resize(numFacesPerType, -1);
397
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 currentElemIdPoint.resize(numPointPerElt, -1);
398
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 oppositeElemIdPoint.resize(numPointPerElt, -1);
399
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 currentElemPoint.resize(numPointPerElt, NULL);
400
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 oppositeElemPoint.resize(numPointPerElt, NULL);
401
402 // define finite element
403
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 defineFiniteElement(eltType);
404
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 initElementArray();
405
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 defineCurrentFiniteElementWithBd(eltType);
406
407 // allocate arrays for assembling the matrix
408
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 allocateArrayForAssembleMatrixRHS(flag);
409
410 // init variables
411
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 initPerElementType(eltType, flag);
412 20 CurrentFiniteElementWithBd* feWithBd = m_listCurrentFiniteElementWithBd[m_iVelocity];
413
414 20 CurrentFiniteElementWithBd* firstCurrentFe = feWithBd;
415 20 CurrentFiniteElementWithBd* firstOppositeFe = oppositeFEWithBd;
416
417 // get informations on the current element
418
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 setElemPoint(eltType, 0, currentElemPoint, currentElemIdPoint, &ielCurrentGlobalSup);
419
420 // update the finite elements
421
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 feWithBd->updateFirstDeriv(0, currentElemPoint);
422
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 feWithBd->updateBdMeasNormal();
423
424 // get the global id of the first geometric element
425
1/2
✓ Branch 2 taken 20 times.
✗ Branch 3 not taken.
 20 ISLocalToGlobalMappingApply(m_mappingElem[m_currentMesh], 1, &ielCurrentLocalGeo, &ielCurrentGlobalGeo);
426
427 // the map to remember what contribution have already been computed
428 20 std::map<felInt, std::vector<bool> > contribDone;
429
2/4
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
 20 addNewFaceOrientedContributor(numFacesPerType, ielCurrentGlobalGeo, contribDone[ielCurrentGlobalGeo]);
430
431 // build the stack to know what is the next element
432
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 std::stack<felInt> nextInStack;
433
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 nextInStack.push(ielCurrentGlobalGeo);
434
435 // get all the faces of the element
436
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 if(dimension() == 2)
437
1/2
✓ Branch 3 taken 20 times.
✗ Branch 4 not taken.
 20 m_mesh[m_currentMesh]->getAllEdgeOfElement(ielCurrentGlobalGeo, idOfFacesCurrent);
438 else
439 ✗ m_mesh[m_currentMesh]->getAllFaceOfElement(ielCurrentGlobalGeo, idOfFacesCurrent);
440
441 // loop over all the element (snake style)
442
2/2
✓ Branch 1 taken 2605 times.
✓ Branch 2 taken 20 times.
 2625 while(!nextInStack.empty()) {
443 // check the faces and use the first one that is an inner face and that has not been done yet.
444 2605 idFaceToDo = -1;
445 2605 allDone = true;
446
3/4
✓ Branch 1 taken 10420 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7815 times.
✓ Branch 5 taken 2605 times.
 10420 for(size_t iface=0; iface<contribDone[ielCurrentGlobalGeo].size(); ++iface) {
447
4/6
✓ Branch 1 taken 7815 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7815 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2890 times.
✓ Branch 8 taken 4925 times.
 7815 if(!contribDone[ielCurrentGlobalGeo][iface]) {
448 // This is an inner face, check if the contribution has already been computed or not
449
2/2
✓ Branch 0 taken 1865 times.
✓ Branch 1 taken 1025 times.
 2890 if(idFaceToDo == -1) {
450 1865 idFaceToDo = iface;
451 } else {
452 1025 allDone = false;
453 }
454 }
455 }
456
457 // update the stack
458
2/2
✓ Branch 0 taken 1020 times.
✓ Branch 1 taken 1585 times.
 2605 if(!allDone)
459
1/2
✓ Branch 1 taken 1020 times.
✗ Branch 2 not taken.
 1020 nextInStack.push(ielCurrentGlobalGeo);
460
461
6/6
✓ Branch 1 taken 2060 times.
✓ Branch 2 taken 545 times.
✓ Branch 3 taken 1040 times.
✓ Branch 4 taken 1020 times.
✓ Branch 5 taken 1040 times.
✓ Branch 6 taken 1565 times.
 2605 if(nextInStack.top() == ielCurrentGlobalGeo && allDone)
462 1040 nextInStack.pop();
463
464
465 // assemble terms
466
2/2
✓ Branch 0 taken 1865 times.
✓ Branch 1 taken 740 times.
 2605 if(idFaceToDo != -1) {
467 // get the opposite id of the element
468 1865 ielOppositeGlobalGeo = ielCurrentGlobalGeo;
469
1/2
✓ Branch 3 taken 1865 times.
✗ Branch 4 not taken.
 1865 adjElt =m_mesh[m_currentMesh]->getAdjElement(ielOppositeGlobalGeo, idFaceToDo);
470
471 // get the type of the opposite element
472 // eltTypeOpp = adjElt.second;
473
474 // update the opposite finite element and set ielOppositeGlobalSup
475
1/2
✓ Branch 1 taken 1865 times.
✗ Branch 2 not taken.
 1865 updateFaceOrientedFEWithBd(oppositeFEWithBd, idOfFacesOpposite, numPointPerElt, ielOppositeGlobalGeo, ielOppositeGlobalSup);
476
477 // find the local id of the edge in the opposite element
478 1865 felInt localIdFaceOpposite = -1;
479
2/2
✓ Branch 1 taken 5595 times.
✓ Branch 2 taken 1865 times.
 7460 for(size_t jface=0; jface<idOfFacesOpposite.size(); ++jface) {
480
2/2
✓ Branch 2 taken 1865 times.
✓ Branch 3 taken 3730 times.
 5595 if(idOfFacesCurrent[idFaceToDo] == idOfFacesOpposite[jface]) {
481 1865 localIdFaceOpposite = jface;
482 }
483 }
484
485 // compute coefficients
486
1/2
✓ Branch 1 taken 1865 times.
✗ Branch 2 not taken.
 1865 CurvilinearFiniteElement* curvFe = &feWithBd->bdEle(idFaceToDo);
487 //CurvilinearFiniteElement* curvFeOP = &oppositeFEWithBd->bdEle(localIdFaceOpposite);
488
1/2
✓ Branch 1 taken 1865 times.
✗ Branch 2 not taken.
 1865 double hK = curvFe->diameter(); //0.5 * (feWithBd->diameter() + oppositeFEWithBd->diameter());
489 //cout << hK - curvFeOP->diameter() << endl;
490
1/2
✓ Branch 1 taken 1865 times.
✗ Branch 2 not taken.
 1865 double gamma_u = FelisceParam::instance().referenceValue1; // cip coefficient
491 1865 double gamma = hK*gamma_u;
492
493 // We can now compute all the integrals.
494 // current element for phi_i and phi_j
495
1/2
✓ Branch 3 taken 1865 times.
✗ Branch 4 not taken.
 1865 m_elementMat[0]->zero();
496
1/2
✓ Branch 4 taken 1865 times.
✗ Branch 5 not taken.
 1865 m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n(gamma, *feWithBd, *feWithBd, idFaceToDo, idFaceToDo, 0, 0, this->dimension()); // block (0,0)
497
1/2
✓ Branch 1 taken 1865 times.
✗ Branch 2 not taken.
 1865 setValueMatrixRHS(ielCurrentGlobalSup, ielCurrentGlobalSup, flag);
498
499 // current element for phi_i and opposite element for phi_j
500
1/2
✓ Branch 3 taken 1865 times.
✗ Branch 4 not taken.
 1865 m_elementMat[0]->zero();
501
1/2
✓ Branch 4 taken 1865 times.
✗ Branch 5 not taken.
 1865 m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n(gamma, *oppositeFEWithBd, *feWithBd, localIdFaceOpposite, idFaceToDo, 0, 0, this->dimension()); // block (0,0)
502
1/2
✓ Branch 1 taken 1865 times.
✗ Branch 2 not taken.
 1865 setValueMatrixRHS(ielCurrentGlobalSup, ielOppositeGlobalSup, flag);
503
504 // opposite element for phi_i and current element for phi_j
505
1/2
✓ Branch 3 taken 1865 times.
✗ Branch 4 not taken.
 1865 m_elementMat[0]->zero();
506
1/2
✓ Branch 4 taken 1865 times.
✗ Branch 5 not taken.
 1865 m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n(gamma, *feWithBd, *oppositeFEWithBd, idFaceToDo, localIdFaceOpposite, 0, 0, this->dimension()); // block (0,0)
507
1/2
✓ Branch 1 taken 1865 times.
✗ Branch 2 not taken.
 1865 setValueMatrixRHS(ielOppositeGlobalSup, ielCurrentGlobalSup, flag);
508
509 // opposite element for phi_i and phi_j
510
1/2
✓ Branch 3 taken 1865 times.
✗ Branch 4 not taken.
 1865 m_elementMat[0]->zero();
511
1/2
✓ Branch 4 taken 1865 times.
✗ Branch 5 not taken.
 1865 m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n(gamma, *oppositeFEWithBd, *oppositeFEWithBd, localIdFaceOpposite, localIdFaceOpposite, 0, 0, this->dimension()); // block (0,0)
512
1/2
✓ Branch 1 taken 1865 times.
✗ Branch 2 not taken.
 1865 setValueMatrixRHS(ielOppositeGlobalSup, ielOppositeGlobalSup, flag);
513
514 // update the map to say that this face is done
515
2/4
✓ Branch 1 taken 1865 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1865 times.
✗ Branch 5 not taken.
 1865 contribDone[ielCurrentGlobalGeo][idFaceToDo] = true;
516
517 // check if the opposite element is on this proc
518
2/2
✓ Branch 2 taken 1550 times.
✓ Branch 3 taken 315 times.
 1865 if(m_eltPart[m_currentMesh][ielOppositeGlobalGeo] == rank) {
519
3/4
✓ Branch 2 taken 1550 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1185 times.
✓ Branch 6 taken 365 times.
 1550 if(contribDone.find(ielOppositeGlobalGeo) == contribDone.end()) {
520 // not found in the map, add it and initialize it.
521
2/4
✓ Branch 1 taken 1185 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1185 times.
✗ Branch 5 not taken.
 1185 addNewFaceOrientedContributor(numFacesPerType, ielOppositeGlobalGeo, contribDone[ielOppositeGlobalGeo]);
522 }
523
2/4
✓ Branch 1 taken 1550 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1550 times.
✗ Branch 5 not taken.
 1550 contribDone[ielOppositeGlobalGeo][localIdFaceOpposite] = true;
524
525 // update the finite element as the opposite finite element.
526 1550 ptmp = feWithBd;
527 1550 feWithBd = oppositeFEWithBd;
528 1550 oppositeFEWithBd = ptmp;
529
530 1550 ielCurrentGlobalGeo = ielOppositeGlobalGeo;
531 1550 ielCurrentGlobalSup = ielOppositeGlobalSup;
532
1/2
✓ Branch 1 taken 1550 times.
✗ Branch 2 not taken.
 1550 idOfFacesCurrent = idOfFacesOpposite;
533 }
534 else {
535 // not on this rank, take the next one in the stack
536
1/2
✓ Branch 1 taken 315 times.
✗ Branch 2 not taken.
 315 if(!nextInStack.empty()) {
537 315 ielCurrentGlobalGeo = nextInStack.top();
538
1/2
✓ Branch 1 taken 315 times.
✗ Branch 2 not taken.
 315 updateFaceOrientedFEWithBd(feWithBd, idOfFacesCurrent, numPointPerElt, ielCurrentGlobalGeo, ielCurrentGlobalSup);
539 }
540 }
541 }
542 else {
543 // All contribution have been already computed on this element
544 // take the next element in the stack
545
2/2
✓ Branch 1 taken 720 times.
✓ Branch 2 taken 20 times.
 740 if(!nextInStack.empty()) {
546 720 ielCurrentGlobalGeo = nextInStack.top();
547
1/2
✓ Branch 1 taken 720 times.
✗ Branch 2 not taken.
 720 updateFaceOrientedFEWithBd(feWithBd, idOfFacesCurrent, numPointPerElt, ielCurrentGlobalGeo, ielCurrentGlobalSup);
548 }
549 }
550 }
551
552 // desallocate array use for assemble with petsc.
553
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 desallocateArrayForAssembleMatrixRHS(flag);
554
555 // desallocate opposite finite elements
556 20 feWithBd = firstCurrentFe;
557
1/2
✓ Branch 0 taken 20 times.
✗ Branch 1 not taken.
 20 delete firstOppositeFe;
558 20 }
559
560 1205 void LinearProblemStokesContinuation::addNewFaceOrientedContributor(felInt size, felInt idElt, std::vector<bool>& vec) {
561 felInt ielTmp;
562 1205 std::pair<bool, GeometricMeshRegion::ElementType> adjElt;
563
564
1/2
✓ Branch 1 taken 1205 times.
✗ Branch 2 not taken.
 1205 vec.resize(size, false);
565
2/2
✓ Branch 0 taken 3615 times.
✓ Branch 1 taken 1205 times.
 4820 for(felInt iface=0; iface<size; ++iface) {
566 // check if this face is an inner face or a boundary
567 3615 ielTmp = idElt;
568
1/2
✓ Branch 3 taken 3615 times.
✗ Branch 4 not taken.
 3615 adjElt =m_mesh[m_currentMesh]->getAdjElement(ielTmp, iface);
569
570
2/2
✓ Branch 0 taken 200 times.
✓ Branch 1 taken 3415 times.
 3615 if(!adjElt.first) {
571 // not an inner face, set the contribution to done.
572
1/2
✓ Branch 1 taken 200 times.
✗ Branch 2 not taken.
 200 vec[iface] = true;
573 }
574 }
575 1205 }
576
577 2900 void LinearProblemStokesContinuation::updateFaceOrientedFEWithBd(CurrentFiniteElementWithBd* fe, std::vector<felInt>& idOfFaces, felInt numPoints, felInt idElt, felInt& idSup) {
578
1/2
✓ Branch 2 taken 2900 times.
✗ Branch 3 not taken.
 2900 std::vector<felInt> elemIdPoint(numPoints, -1);
579
1/2
✓ Branch 2 taken 2900 times.
✗ Branch 3 not taken.
 2900 std::vector<Point*> elemPoint(numPoints, NULL);
580
581 // get information on the opposite element
582 // same as the function "setElemPoint" but here ielOppositeGlobalGeo is global
583 // we assume that the supportDofMesh is the same for all unknown (like in setElemPoint)
584
1/2
✓ Branch 2 taken 2900 times.
✗ Branch 3 not taken.
 2900 m_supportDofUnknown[0].getIdElementSupport(idElt, idSup);
585
1/2
✓ Branch 3 taken 2900 times.
✗ Branch 4 not taken.
 2900 m_mesh[m_currentMesh]->getOneElement(idElt, elemIdPoint);
586
2/2
✓ Branch 0 taken 8700 times.
✓ Branch 1 taken 2900 times.
 11600 for (felInt iPoint=0; iPoint<numPoints; ++iPoint)
587 8700 elemPoint[iPoint] = &(m_mesh[m_currentMesh]->listPoints()[elemIdPoint[iPoint]]);
588
589 // update the finite elements
590
1/2
✓ Branch 1 taken 2900 times.
✗ Branch 2 not taken.
 2900 fe->updateFirstDeriv(0, elemPoint);
591
1/2
✓ Branch 1 taken 2900 times.
✗ Branch 2 not taken.
 2900 fe->updateBdMeasNormal();
592
593 // get all the faces of the element
594
1/2
✓ Branch 1 taken 2900 times.
✗ Branch 2 not taken.
 2900 if(dimension() == 2)
595
1/2
✓ Branch 3 taken 2900 times.
✗ Branch 4 not taken.
 2900 m_mesh[m_currentMesh]->getAllEdgeOfElement(idElt, idOfFaces);
596 else
597 ✗ m_mesh[m_currentMesh]->getAllFaceOfElement(idElt, idOfFaces);
598 2900 }
599 }
600