GCC Code Coverage Report

Directory: ./
File: Solver/linearProblemPoissonContinuation.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 247 254 97.2%
Branches: 207 346 59.8%
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/linearProblemPoissonContinuation.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 linearProblemPoisson.cpp
30 \date 17/02/2022
31 \brief Continuation method for Poisson's equation
32 */
33 namespace felisce {
34 4 LinearProblemPoissonContinuation::LinearProblemPoissonContinuation():
35
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")
36 {
37 4 }
38
39 4 void LinearProblemPoissonContinuation::readData(IO& io,double iteration) {
40 // Read variable from file
41 4 std::vector<double> potData;
42
1/2
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 4 potData.resize(m_mesh[m_currentMesh]->numPoints(),0.);
43
1/2
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
 4 io.readVariable(0,iteration, potData.data(), potData.size());
44
45
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 m_potData.duplicateFrom(this->sequentialSolution());
46
47
1/2
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
 4 std::vector<PetscInt> ao_potData( potData.size());
48
2/2
✓ Branch 1 taken 568 times.
✓ Branch 2 taken 4 times.
 572 for (size_t i=0; i < potData.size(); ++i) {
49 568 ao_potData[i]=i;
50 }
51
52
1/2
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 4 AOApplicationToPetsc(this->ao(),potData.size(),ao_potData.data());
53
2/4
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
 4 VecSetValues(m_potData.toPetsc(),potData.size(),ao_potData.data(), potData.data(), INSERT_VALUES);
54 4 }
55
56 4 void LinearProblemPoissonContinuation::initialize(std::vector<GeometricMeshRegion::Pointer>& mesh, FelisceTransient::Pointer fstransient, MPI_Comm& comm, bool doUseSNES) {
57
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 LinearProblem::initialize(mesh,comm, doUseSNES);
58 4 this->m_fstransient = fstransient;
59
60
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 std::vector<PhysicalVariable> listVariable(2);
61
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 std::vector<size_t> listNumComp(2);
62 4 listVariable[0] = temperature; // primal variable, the discrete solution u_h
63 4 listNumComp[0] = 1;
64 4 listVariable[1] = potThorax; // dummy name for dual variable, the Langrange multiplier z_h
65 4 listNumComp[1] = 1;
66
67 //define unknown of the linear system.
68
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 m_listUnknown.push_back(temperature);
69
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 m_listUnknown.push_back(potThorax);
70
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 definePhysicalVariable(listVariable,listNumComp);
71
72
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 m_iTemperature = m_listVariable.getVariableIdList(temperature);
73
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 m_iPotThorax = m_listVariable.getVariableIdList(potThorax);
74 4 }
75
76 8 void LinearProblemPoissonContinuation::initPerElementType(ElementType /*eltType*/, FlagMatrixRHS /*flagMatrixRHS*/) {
77 8 m_feTemp = m_listCurrentFiniteElement[m_iTemperature];
78
79 // initialize elemField to store the data
80 8 m_elemField.initialize(DOF_FIELD, *m_feTemp, 1);
81 8 }
82
83 242 void LinearProblemPoissonContinuation::computeElementArray(const std::vector<Point*>& elemPoint, const std::vector<felInt>& /*elemIdPoint*/, felInt& iel, FlagMatrixRHS flagMatrixRHS) {
84 242 m_feTemp->updateFirstDeriv(0, elemPoint);
85
86 242 const double gamma_M = FelisceParam::instance(this->instanceIndex()).referenceValue2; // data fidelity coefficient
87 242 const int labelDataSubdomain = FelisceParam::instance(this->instanceIndex()).testCase; // label of the subdomain where data is given; 0 for the whole domain
88
89 // defining the block matrix of the system
90
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 242 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 242 if (flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs || flagMatrixRHS == FlagMatrixRHS::only_matrix) {
91 // block (0,0): gamma_M * int_\Omega_data u * v; in this block we will also add the cip stabilization
92
1/2
✓ Branch 0 taken 242 times.
✗ Branch 1 not taken.
 242 if (labelDataSubdomain == 0) {
93 242 m_elementMat[0]->phi_i_phi_j(gamma_M, *m_feTemp, 0, 0, 1);
94 }
95 ✗ else if (m_currentLabel == labelDataSubdomain) {
96 ✗ m_elementMat[0]->phi_i_phi_j(gamma_M, *m_feTemp, 0, 0, 1);
97 }
98
99 // block (0,1): int_\Omega grad v . grad z
100 242 m_elementMat[0]->grad_phi_i_grad_phi_j(1., *m_feTemp, 0, 1, 1);
101
102 // block (1,0): int_\Omega grad u . grad w
103 242 m_elementMat[0]->grad_phi_i_grad_phi_j(1., *m_feTemp, 1, 0, 1);
104
105 // block (1,1): -int_\Omega grad z . grad w
106 242 m_elementMat[0]->grad_phi_i_grad_phi_j(-1., *m_feTemp, 1, 1, 1);
107 }
108
109 // defining the rhs of the system
110
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 242 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 242 if (flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs || flagMatrixRHS == FlagMatrixRHS::only_rhs) {
111 // assuming that f is zero
112
113 // rhs data term
114 242 m_elemField.setValue(m_potData, *m_feTemp, iel, m_iTemperature, m_ao, dof());
115
1/2
✓ Branch 0 taken 242 times.
✗ Branch 1 not taken.
 242 if (labelDataSubdomain == 0) {
116 242 m_elementVector[0]->source(gamma_M, *m_feTemp, m_elemField, 0, 1);
117 }
118 ✗ else if (m_currentLabel == labelDataSubdomain) {
119 ✗ m_elementVector[0]->source(gamma_M, *m_feTemp, m_elemField, 0, 1);
120 }
121 }
122 242 }
123
124 4 void LinearProblemPoissonContinuation::userChangePattern(int numProc, int rankProc) {
125 IGNORE_UNUSED_ARGUMENT(numProc);
126 IGNORE_UNUSED_ARGUMENT(rankProc);
127 // compute initial (uniform) repartition of dof on processes
128
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 felInt numDofByProc = m_numDof/MpiInfo::numProc();
129
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());
130
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) {
131
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)
132 4 numDofPerProc[i] = m_numDof - i*numDofByProc;
133 else
134 12 numDofPerProc[i] = numDofByProc;
135 }
136
137 4 felInt shiftDof = 0;
138
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)
139 6 shiftDof += numDofPerProc[i];
140
141
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 std::vector<felInt> rankDof(m_numDof, -1);
142
3/4
✓ Branch 1 taken 288 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 284 times.
✓ Branch 5 taken 4 times.
 288 for(felInt i=0; i<numDofPerProc[MpiInfo::rankProc()]; ++i)
143
1/2
✓ Branch 1 taken 284 times.
✗ Branch 2 not taken.
 284 rankDof[i + shiftDof] = MpiInfo::rankProc();
144
145 // build the edges or faces depending on the dimension (for the global mesh)
146 // In this function, "faces" refers to either edges or faces.
147
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 if(dimension() == 2)
148
1/2
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
 4 m_mesh[m_currentMesh]->buildEdges();
149 else
150 ✗ m_mesh[m_currentMesh]->buildFaces();
151
152 // variables
153 GeometricMeshRegion::ElementType eltType, eltTypeOpp;
154 felInt idVar1, idVar2;
155 felInt node1, node2;
156 felInt numFacesPerElement, numEltPerLabel;
157 felInt ielSupport, jelSupport;
158 felInt ielCurrentGeo, ielOppositeGeo;
159 4 std::vector<felInt> vecSupport, vecSupportOpposite;
160
1/2
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 4 std::vector<felInt> numElement(m_mesh[m_currentMesh]->m_numTypesOfElement, 0);
161
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 std::vector< std::set<felInt> > nodesNeighborhood(m_numDof);
162
163 // zeroth loop on the unknown of the linear problem
164
2/2
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 4 times.
 12 for (size_t iUnknown1 = 0; iUnknown1 < m_listUnknown.size(); ++iUnknown1) {
165 8 idVar1 = m_listUnknown.idVariable(iUnknown1);
166
167 8 ielCurrentGeo = 0;
168
169 // first loop on element type
170
2/2
✓ Branch 4 taken 8 times.
✓ Branch 5 taken 8 times.
 16 for (size_t i=0; i<m_mesh[m_currentMesh]->bagElementTypeDomain().size(); ++i) {
171 8 eltType = m_mesh[m_currentMesh]->bagElementTypeDomain()[i];
172 8 const GeoElement* geoEle = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second;
173 8 numElement[eltType] = 0;
174 8 numFacesPerElement = geoEle->numBdEle();
175
176 // second loop on region of the mesh.
177
2/2
✓ Branch 10 taken 8 times.
✓ Branch 11 taken 8 times.
 16 for(GeometricMeshRegion::IntRefToBegEndIndex_type::const_iterator itRef =m_mesh[m_currentMesh]->intRefToBegEndMaps[eltType].begin(); itRef !=m_mesh[m_currentMesh]->intRefToBegEndMaps[eltType].end(); itRef++) {
178 8 numEltPerLabel = itRef->second.second;
179
180 // Third loop on element
181
2/2
✓ Branch 0 taken 1936 times.
✓ Branch 1 taken 8 times.
 1944 for (felInt iel = 0; iel < numEltPerLabel; iel++) {
182
1/2
✓ Branch 3 taken 1936 times.
✗ Branch 4 not taken.
 1936 m_supportDofUnknown[iUnknown1].getIdElementSupport(eltType, numElement[eltType], vecSupport);
183
184
2/2
✓ Branch 1 taken 1936 times.
✓ Branch 2 taken 1936 times.
 3872 for(size_t ielSup=0; ielSup<vecSupport.size(); ++ielSup) {
185 1936 ielSupport = vecSupport[ielSup];
186
187 // for all support dof in this support element
188
2/2
✓ Branch 2 taken 5808 times.
✓ Branch 3 taken 1936 times.
 7744 for (felInt iSup = 0; iSup < m_supportDofUnknown[iUnknown1].getNumSupportDof(ielSupport); ++iSup) {
189 // for all component of the unknown
190
2/2
✓ Branch 2 taken 5808 times.
✓ Branch 3 taken 5808 times.
 11616 for (size_t iComp = 0; iComp < m_listVariable[idVar1].numComponent(); iComp++) {
191 // get the global id of the support dof
192
1/2
✓ Branch 2 taken 5808 times.
✗ Branch 3 not taken.
 5808 dof().loc2glob(ielSupport, iSup, idVar1, iComp, node1);
193
194 // if this node is on this process
195
3/4
✓ Branch 2 taken 5808 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1452 times.
✓ Branch 5 taken 4356 times.
 5808 if(rankDof[node1] == MpiInfo::rankProc()) {
196 // loop over the boundaries of the element
197
2/2
✓ Branch 0 taken 4356 times.
✓ Branch 1 taken 1452 times.
 5808 for(felInt iface=0; iface < numFacesPerElement; ++iface) {
198 // check if this face is an inner face or a boundary
199 4356 ielOppositeGeo = ielCurrentGeo;
200 4356 eltTypeOpp = eltType;
201
1/2
✓ Branch 3 taken 4356 times.
✗ Branch 4 not taken.
 4356 bool isAnInnerFace =m_mesh[m_currentMesh]->getAdjElement(eltTypeOpp, ielOppositeGeo, iface);
202
203
2/2
✓ Branch 0 taken 4116 times.
✓ Branch 1 taken 240 times.
 4356 if(isAnInnerFace) {
204 // for all unknown
205
2/2
✓ Branch 1 taken 8232 times.
✓ Branch 2 taken 4116 times.
 12348 for (size_t iUnknown2 = 0; iUnknown2 < m_listUnknown.size(); ++iUnknown2) {
206 8232 idVar2 = m_listUnknown.idVariable(iUnknown2);
207
208 // get the support element of the opposite element
209
1/2
✓ Branch 2 taken 8232 times.
✗ Branch 3 not taken.
 8232 m_supportDofUnknown[iUnknown2].getIdElementSupport(ielOppositeGeo, vecSupportOpposite);
210
211 // for all support element
212
2/2
✓ Branch 1 taken 8232 times.
✓ Branch 2 taken 8232 times.
 16464 for(size_t jelSup=0; jelSup<vecSupportOpposite.size(); ++jelSup) {
213 8232 jelSupport = vecSupportOpposite[jelSup];
214
215 // for all support dof in this support element
216
2/2
✓ Branch 2 taken 24696 times.
✓ Branch 3 taken 8232 times.
 32928 for (felInt jSup = 0; jSup < m_supportDofUnknown[iUnknown2].getNumSupportDof(jelSupport); ++jSup) {
217
218 // for all component of the second unknown
219
2/2
✓ Branch 2 taken 24696 times.
✓ Branch 3 taken 24696 times.
 49392 for (size_t jComp = 0; jComp < m_listVariable[idVar2].numComponent(); jComp++) {
220
221 // If the two current components are connected
222
1/2
✓ Branch 2 taken 24696 times.
✗ Branch 3 not taken.
 24696 felInt iConnect = dof().getNumGlobComp(iUnknown1, iComp);
223
1/2
✓ Branch 2 taken 24696 times.
✗ Branch 3 not taken.
 24696 felInt jConnect = dof().getNumGlobComp(iUnknown2, jComp);
224
2/4
✓ Branch 2 taken 24696 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 24696 times.
✗ Branch 5 not taken.
 24696 if (m_listUnknown.mask()(iConnect, jConnect) > 0) {
225 // get the global id of the second support dof
226
1/2
✓ Branch 2 taken 24696 times.
✗ Branch 3 not taken.
 24696 dof().loc2glob(jelSupport, jSup, idVar2, jComp, node2);
227
228 // remove diagonal term to define pattern and use it in Parmetis.
229 // if the two global ids are different, they are neighboors
230
2/2
✓ Branch 0 taken 21952 times.
✓ Branch 1 taken 2744 times.
 24696 if(node1 != node2)
231
1/2
✓ Branch 2 taken 21952 times.
✗ Branch 3 not taken.
 21952 nodesNeighborhood[node1].insert(node2);
232 }
233 }
234 }
235 }
236 }
237 }
238 }
239 }
240 }
241 }
242 }
243 1936 ++ielCurrentGeo;
244 1936 ++numElement[eltType];
245 }
246 }
247 }
248 }
249
250 // Store the pattern in CSR style
251 4 std::vector<felInt> iCSR, jCSR;
252 4 felInt dofSize = 0;
253 4 felInt cptDof = 0;
254 felInt pos;
255
256
1/2
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 4 iCSR.resize(dof().pattern().numRows() + 1, 0);
257
2/2
✓ Branch 1 taken 1136 times.
✓ Branch 2 taken 4 times.
 1140 for(size_t iNode=0; iNode<nodesNeighborhood.size(); ++iNode)
258 1136 dofSize += nodesNeighborhood[iNode].size();
259
260
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 jCSR.resize(dofSize, 0);
261
2/2
✓ Branch 1 taken 1136 times.
✓ Branch 2 taken 4 times.
 1140 for(size_t iNode=0; iNode<nodesNeighborhood.size(); ++iNode) {
262
3/4
✓ Branch 2 taken 1136 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 284 times.
✓ Branch 5 taken 852 times.
 1136 if(rankDof[iNode] == MpiInfo::rankProc()) {
263 284 iCSR[cptDof + 1] = iCSR[cptDof] + nodesNeighborhood[iNode].size();
264 284 pos = 0;
265
2/2
✓ Branch 5 taken 6028 times.
✓ Branch 6 taken 284 times.
 6312 for(std::set<felInt>::iterator it=nodesNeighborhood[iNode].begin(); it != nodesNeighborhood[iNode].end(); ++it) {
266 6028 jCSR[iCSR[cptDof] + pos] = *it;
267 6028 ++pos;
268 }
269 284 ++cptDof;
270 }
271 }
272
273 // Now, call merge pattern
274
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 dof().mergeGlobalPattern(iCSR, jCSR);
275 4 }
276
277 4 void LinearProblemPoissonContinuation::assembleFaceOrientedStabilization() {
278 // definition of variables
279 4 std::pair<bool, GeometricMeshRegion::ElementType> adjElt;
280 GeometricMeshRegion::ElementType eltType; // Type of element
281
282 4 felInt ielCurrentLocalGeo = 0; // local geometric id of the current element
283 4 felInt ielCurrentGlobalGeo = 0; // global geometric id of the current element
284 4 felInt ielOppositeGlobalGeo = 0; // global geometric id of the opposite element
285 felInt numFacesPerType; // number of faces of an element of a given type
286 felInt numPointPerElt; // number of vertices by element
287 felInt ielCurrentGlobalSup; // global support element id of the current element
288 felInt ielOppositeGlobalSup; // global support element id of the opposite element
289
290 4 std::vector<felInt> idOfFacesCurrent; // ids of the current element edges
291 4 std::vector<felInt> idOfFacesOpposite; // ids of the opposite element edges
292 4 std::vector<felInt> currentElemIdPoint; // ids of the vertices of the current element
293 4 std::vector<felInt> oppositeElemIdPoint; // ids of the vertices of the opposite element
294 4 std::vector<Point*> currentElemPoint; // point coordinates of the current element vertices
295 4 std::vector<Point*> oppositeElemPoint; // point coordinates of the opposite element vertices
296
297 4 FlagMatrixRHS flag = FlagMatrixRHS::only_matrix; // flag to only assemble the matrix
298 felInt idFaceToDo;
299 4 bool allDone = false;
300
301
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 ElementField elemFieldAdvFace;
302
303 felInt rank;
304
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
305
306 // bag element type vector
307 4 const std::vector<ElementType>& bagElementTypeDomain = m_meshLocal[m_currentMesh]->bagElementTypeDomain();
308 4 eltType = bagElementTypeDomain[0];
309
310
311 // Finite Element With Bd for the opposite element
312 CurrentFiniteElementWithBd* ptmp;
313 CurrentFiniteElementWithBd* oppositeFEWithBd;
314
315
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 m_temperature = &m_listVariable[m_listVariable.getVariableIdList(temperature)];
316
317 4 const GeoElement *geoEle = m_mesh[m_currentMesh]->eltEnumToFelNameGeoEle[eltType].second;
318 4 felInt typeOfFEVel = m_temperature->finiteElementType();
319
320
1/2
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
 4 const RefElement *refEleVel = geoEle->defineFiniteEle(eltType, typeOfFEVel, *m_mesh[m_currentMesh]);
321
2/4
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 4 times.
✗ Branch 7 not taken.
 4 oppositeFEWithBd = new CurrentFiniteElementWithBd(*refEleVel, *geoEle, m_temperature->degreeOfExactness(), m_temperature->degreeOfExactness());
322
323 // initializing variables
324 4 numPointPerElt = m_meshLocal[m_currentMesh]->m_numPointsPerElt[eltType];
325 4 numFacesPerType = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second->numBdEle();
326
327 // resize of vectors
328
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 idOfFacesCurrent.resize(numFacesPerType, -1);
329
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 idOfFacesOpposite.resize(numFacesPerType, -1);
330
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 currentElemIdPoint.resize(numPointPerElt, -1);
331
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 oppositeElemIdPoint.resize(numPointPerElt, -1);
332
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 currentElemPoint.resize(numPointPerElt, NULL);
333
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 oppositeElemPoint.resize(numPointPerElt, NULL);
334
335 // define finite element
336
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 defineFiniteElement(eltType);
337
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 initElementArray();
338
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 defineCurrentFiniteElementWithBd(eltType);
339
340 // allocate arrays for assembling the matrix
341
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 allocateArrayForAssembleMatrixRHS(flag);
342
343 // init variables
344
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 initPerElementType(eltType, flag);
345 4 CurrentFiniteElementWithBd* feWithBd = m_listCurrentFiniteElementWithBd[m_iTemperature];
346
347 4 CurrentFiniteElementWithBd* firstCurrentFe = feWithBd;
348 4 CurrentFiniteElementWithBd* firstOppositeFe = oppositeFEWithBd;
349
350 // get informations on the current element
351
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 setElemPoint(eltType, 0, currentElemPoint, currentElemIdPoint, &ielCurrentGlobalSup);
352
353 // update the finite elements
354
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 feWithBd->updateFirstDeriv(0, currentElemPoint);
355
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 feWithBd->updateBdMeasNormal();
356
357 // get the global id of the first geometric element
358
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 ISLocalToGlobalMappingApply(m_mappingElem[m_currentMesh], 1, &ielCurrentLocalGeo, &ielCurrentGlobalGeo);
359
360 // the map to remember what contribution have already been computed
361 4 std::map<felInt, std::vector<bool> > contribDone;
362
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 4 addNewFaceOrientedContributor(numFacesPerType, ielCurrentGlobalGeo, contribDone[ielCurrentGlobalGeo]);
363
364 // build the stack to know what is the next element
365
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 std::stack<felInt> nextInStack;
366
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 nextInStack.push(ielCurrentGlobalGeo);
367
368 // get all the faces of the element
369
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 if(dimension() == 2)
370
1/2
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
 4 m_mesh[m_currentMesh]->getAllEdgeOfElement(ielCurrentGlobalGeo, idOfFacesCurrent);
371 else
372 ✗ m_mesh[m_currentMesh]->getAllFaceOfElement(ielCurrentGlobalGeo, idOfFacesCurrent);
373
374 // loop over all the element (snake style)
375
2/2
✓ Branch 1 taken 530 times.
✓ Branch 2 taken 4 times.
 534 while(!nextInStack.empty()) {
376 // check the faces and use the first one that is an inner face and that has not been done yet.
377 530 idFaceToDo = -1;
378 530 allDone = true;
379
3/4
✓ Branch 1 taken 2120 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1590 times.
✓ Branch 5 taken 530 times.
 2120 for(size_t iface=0; iface<contribDone[ielCurrentGlobalGeo].size(); ++iface) {
380
4/6
✓ Branch 1 taken 1590 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1590 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 577 times.
✓ Branch 8 taken 1013 times.
 1590 if(!contribDone[ielCurrentGlobalGeo][iface]) {
381 // This is an inner face, check if the contribution has already been computed or not
382
2/2
✓ Branch 0 taken 370 times.
✓ Branch 1 taken 207 times.
 577 if(idFaceToDo == -1) {
383 370 idFaceToDo = iface;
384 } else {
385 207 allDone = false;
386 }
387 }
388 }
389
390 // update the stack
391
2/2
✓ Branch 0 taken 206 times.
✓ Branch 1 taken 324 times.
 530 if(!allDone)
392
1/2
✓ Branch 1 taken 206 times.
✗ Branch 2 not taken.
 206 nextInStack.push(ielCurrentGlobalGeo);
393
394
6/6
✓ Branch 1 taken 416 times.
✓ Branch 2 taken 114 times.
✓ Branch 3 taken 210 times.
✓ Branch 4 taken 206 times.
✓ Branch 5 taken 210 times.
✓ Branch 6 taken 320 times.
 530 if(nextInStack.top() == ielCurrentGlobalGeo && allDone)
395 210 nextInStack.pop();
396
397
398 // assemble terms
399
2/2
✓ Branch 0 taken 370 times.
✓ Branch 1 taken 160 times.
 530 if(idFaceToDo != -1) {
400 // get the opposite id of the element
401 370 ielOppositeGlobalGeo = ielCurrentGlobalGeo;
402
1/2
✓ Branch 3 taken 370 times.
✗ Branch 4 not taken.
 370 adjElt =m_mesh[m_currentMesh]->getAdjElement(ielOppositeGlobalGeo, idFaceToDo);
403
404 // get the type of the opposite element
405 // eltTypeOpp = adjElt.second;
406
407 // update the opposite finite element and set ielOppositeGlobalSup
408
1/2
✓ Branch 1 taken 370 times.
✗ Branch 2 not taken.
 370 updateFaceOrientedFEWithBd(oppositeFEWithBd, idOfFacesOpposite, numPointPerElt, ielOppositeGlobalGeo, ielOppositeGlobalSup);
409
410 // find the local id of the edge in the opposite element
411 370 felInt localIdFaceOpposite = -1;
412
2/2
✓ Branch 1 taken 1110 times.
✓ Branch 2 taken 370 times.
 1480 for(size_t jface=0; jface<idOfFacesOpposite.size(); ++jface) {
413
2/2
✓ Branch 2 taken 370 times.
✓ Branch 3 taken 740 times.
 1110 if(idOfFacesCurrent[idFaceToDo] == idOfFacesOpposite[jface]) {
414 370 localIdFaceOpposite = jface;
415 }
416 }
417
418 // compute coefficients
419
1/2
✓ Branch 1 taken 370 times.
✗ Branch 2 not taken.
 370 const CurvilinearFiniteElement* curvFe = &feWithBd->bdEle(idFaceToDo);
420 //CurvilinearFiniteElement* curvFeOP = &oppositeFEWithBd->bdEle(localIdFaceOpposite);
421
1/2
✓ Branch 1 taken 370 times.
✗ Branch 2 not taken.
 370 double hK = curvFe->diameter(); //0.5 * (feWithBd->diameter() + oppositeFEWithBd->diameter());
422
2/4
✓ Branch 1 taken 370 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 370 times.
✗ Branch 5 not taken.
 370 const double gamma_u = FelisceParam::instance(this->instanceIndex()).referenceValue1; // cip stabilization parameter
423 370 double gamma = hK*gamma_u;
424
425 // We can now compute all the integrals.
426 // current element for phi_i and phi_j
427
1/2
✓ Branch 3 taken 370 times.
✗ Branch 4 not taken.
 370 m_elementMat[0]->zero();
428
1/2
✓ Branch 3 taken 370 times.
✗ Branch 4 not taken.
 370 m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n(gamma, *feWithBd, *feWithBd, idFaceToDo, idFaceToDo, 0, 0, 1); // block (0,0)
429
1/2
✓ Branch 1 taken 370 times.
✗ Branch 2 not taken.
 370 setValueMatrixRHS(ielCurrentGlobalSup, ielCurrentGlobalSup, flag);
430
431 // current element for phi_i and opposite element for phi_j
432
1/2
✓ Branch 3 taken 370 times.
✗ Branch 4 not taken.
 370 m_elementMat[0]->zero();
433
1/2
✓ Branch 3 taken 370 times.
✗ Branch 4 not taken.
 370 m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n(gamma, *oppositeFEWithBd, *feWithBd, localIdFaceOpposite, idFaceToDo, 0, 0, 1); // block (0,0)
434
1/2
✓ Branch 1 taken 370 times.
✗ Branch 2 not taken.
 370 setValueMatrixRHS(ielCurrentGlobalSup, ielOppositeGlobalSup, flag);
435
436 // opposite element for phi_i and current element for phi_j
437
1/2
✓ Branch 3 taken 370 times.
✗ Branch 4 not taken.
 370 m_elementMat[0]->zero();
438
1/2
✓ Branch 3 taken 370 times.
✗ Branch 4 not taken.
 370 m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n(gamma, *feWithBd, *oppositeFEWithBd, idFaceToDo, localIdFaceOpposite, 0, 0, 1); // block (0,0)
439
1/2
✓ Branch 1 taken 370 times.
✗ Branch 2 not taken.
 370 setValueMatrixRHS(ielOppositeGlobalSup, ielCurrentGlobalSup, flag);
440
441 // opposite element for phi_i and phi_j
442
1/2
✓ Branch 3 taken 370 times.
✗ Branch 4 not taken.
 370 m_elementMat[0]->zero();
443
1/2
✓ Branch 3 taken 370 times.
✗ Branch 4 not taken.
 370 m_elementMat[0]->grad_phi_j_dot_n_grad_psi_i_dot_n(gamma, *oppositeFEWithBd, *oppositeFEWithBd, localIdFaceOpposite, localIdFaceOpposite, 0, 0, 1); // block (0,0)
444
1/2
✓ Branch 1 taken 370 times.
✗ Branch 2 not taken.
 370 setValueMatrixRHS(ielOppositeGlobalSup, ielOppositeGlobalSup, flag);
445
446 // update the map to say that this face is done
447
2/4
✓ Branch 1 taken 370 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 370 times.
✗ Branch 5 not taken.
 370 contribDone[ielCurrentGlobalGeo][idFaceToDo] = true;
448
449 // check if the opposite element is on this proc
450
2/2
✓ Branch 2 taken 316 times.
✓ Branch 3 taken 54 times.
 370 if(m_eltPart[m_currentMesh][ielOppositeGlobalGeo] == rank) {
451
3/4
✓ Branch 2 taken 316 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 238 times.
✓ Branch 6 taken 78 times.
 316 if(contribDone.find(ielOppositeGlobalGeo) == contribDone.end()) {
452 // not found in the map, add it and initialize it.
453
2/4
✓ Branch 1 taken 238 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 238 times.
✗ Branch 5 not taken.
 238 addNewFaceOrientedContributor(numFacesPerType, ielOppositeGlobalGeo, contribDone[ielOppositeGlobalGeo]);
454 }
455
2/4
✓ Branch 1 taken 316 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 316 times.
✗ Branch 5 not taken.
 316 contribDone[ielOppositeGlobalGeo][localIdFaceOpposite] = true;
456
457 // update the finite element as the opposite finite element.
458 316 ptmp = feWithBd;
459 316 feWithBd = oppositeFEWithBd;
460 316 oppositeFEWithBd = ptmp;
461
462 316 ielCurrentGlobalGeo = ielOppositeGlobalGeo;
463 316 ielCurrentGlobalSup = ielOppositeGlobalSup;
464
1/2
✓ Branch 1 taken 316 times.
✗ Branch 2 not taken.
 316 idOfFacesCurrent = idOfFacesOpposite;
465 }
466 else {
467 // not on this rank, take the next one in the stack
468
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
 54 if(!nextInStack.empty()) {
469 54 ielCurrentGlobalGeo = nextInStack.top();
470
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
 54 updateFaceOrientedFEWithBd(feWithBd, idOfFacesCurrent, numPointPerElt, ielCurrentGlobalGeo, ielCurrentGlobalSup);
471 }
472 }
473 }
474 else {
475 // All contribution have been already computed on this element
476 // take the next element in the stack
477
2/2
✓ Branch 1 taken 156 times.
✓ Branch 2 taken 4 times.
 160 if(!nextInStack.empty()) {
478 156 ielCurrentGlobalGeo = nextInStack.top();
479
1/2
✓ Branch 1 taken 156 times.
✗ Branch 2 not taken.
 156 updateFaceOrientedFEWithBd(feWithBd, idOfFacesCurrent, numPointPerElt, ielCurrentGlobalGeo, ielCurrentGlobalSup);
480 }
481 }
482 }
483
484 // desallocate array use for assemble with petsc.
485
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 desallocateArrayForAssembleMatrixRHS(flag);
486
487 // desallocate opposite finite elements
488 4 feWithBd = firstCurrentFe;
489
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
 4 delete firstOppositeFe;
490 4 }
491
492 242 void LinearProblemPoissonContinuation::addNewFaceOrientedContributor(felInt size, felInt idElt, std::vector<bool>& vec) {
493 felInt ielTmp;
494 242 std::pair<bool, GeometricMeshRegion::ElementType> adjElt;
495
496
1/2
✓ Branch 1 taken 242 times.
✗ Branch 2 not taken.
 242 vec.resize(size, false);
497
2/2
✓ Branch 0 taken 726 times.
✓ Branch 1 taken 242 times.
 968 for(felInt iface=0; iface<size; ++iface) {
498 // check if this face is an inner face or a boundary
499 726 ielTmp = idElt;
500
1/2
✓ Branch 3 taken 726 times.
✗ Branch 4 not taken.
 726 adjElt = m_mesh[m_currentMesh]->getAdjElement(ielTmp, iface);
501
502
2/2
✓ Branch 0 taken 40 times.
✓ Branch 1 taken 686 times.
 726 if(!adjElt.first) {
503 // not an inner face, set the contribution to done.
504
1/2
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
 40 vec[iface] = true;
505 }
506 }
507 242 }
508
509 580 void LinearProblemPoissonContinuation::updateFaceOrientedFEWithBd(CurrentFiniteElementWithBd* fe, std::vector<felInt>& idOfFaces, felInt numPoints, felInt idElt, felInt& idSup) {
510
1/2
✓ Branch 2 taken 580 times.
✗ Branch 3 not taken.
 580 std::vector<felInt> elemIdPoint(numPoints, -1);
511
1/2
✓ Branch 2 taken 580 times.
✗ Branch 3 not taken.
 580 std::vector<Point*> elemPoint(numPoints, NULL);
512
513 // get information on the opposite element
514 // same as the function "setElemPoint" but here ielOppositeGlobalGeo is global
515 // we assume that the supportDofMesh is the same for all unknown (like in setElemPoint)
516
1/2
✓ Branch 2 taken 580 times.
✗ Branch 3 not taken.
 580 m_supportDofUnknown[0].getIdElementSupport(idElt, idSup);
517
1/2
✓ Branch 3 taken 580 times.
✗ Branch 4 not taken.
 580 m_mesh[m_currentMesh]->getOneElement(idElt, elemIdPoint);
518
2/2
✓ Branch 0 taken 1740 times.
✓ Branch 1 taken 580 times.
 2320 for (felInt iPoint=0; iPoint<numPoints; ++iPoint)
519 1740 elemPoint[iPoint] = &(m_mesh[m_currentMesh]->listPoints()[elemIdPoint[iPoint]]);
520
521 // update the finite elements
522
1/2
✓ Branch 1 taken 580 times.
✗ Branch 2 not taken.
 580 fe->updateFirstDeriv(0, elemPoint);
523
1/2
✓ Branch 1 taken 580 times.
✗ Branch 2 not taken.
 580 fe->updateBdMeasNormal();
524
525 // get all the faces of the element
526
1/2
✓ Branch 1 taken 580 times.
✗ Branch 2 not taken.
 580 if(dimension() == 2)
527
1/2
✓ Branch 3 taken 580 times.
✗ Branch 4 not taken.
 580 m_mesh[m_currentMesh]->getAllEdgeOfElement(idElt, idOfFaces);
528 else
529 ✗ m_mesh[m_currentMesh]->getAllFaceOfElement(idElt, idOfFaces);
530 580 }
531 }
532