GCC Code Coverage Report

Directory: ./
File: Solver/linearProblemLung.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 0 421 0.0%
Branches: 0 852 0.0%
Line Branch Exec Source
1 // ______ ______ _ _ _____ ______
2 // | ____| ____| | (_)/ ____| | ____|
3 // | |__ | |__ | | _| (___ ___| |__
4 // | __| | __| | | | |\___ \ / __| __|
5 // | | | |____| |____| |____) | (__| |____
6 // |_| |______|______|_|_____/ \___|______|
7 // Finite Elements for Life Sciences and Engineering
8 //
9 // License: LGL2.1 License
10 // FELiScE default license: LICENSE in root folder
11 //
12 // Main authors:
13 //
14
15 // System includes
16
17 // External includes
18
19 // Project includes
20 #include "DegreeOfFreedom/boundaryCondition.hpp"
21 #include "Solver/linearProblemLung.hpp"
22
23 namespace felisce
24 {
25
26 ✗ void LinearProblemLung::computeElementArrayLung(Vec & divVel,Mat & divBasis, PetscVector& rhsMinus ,int currentLabel,PetscVector& prevVel, PetscVector& prevDis, double currPre,const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint, felInt& iel, FlagMatrixRHS flagMatrixRHS)
27 {
28 (void) divVel;
29 (void) divBasis;
30 (void) rhsMinus;
31 (void) currentLabel;
32 (void) prevVel;
33 (void) prevDis;
34 (void) currPre;
35 (void) elemPoint;
36 (void) elemIdPoint;
37 (void) iel;
38 (void) flagMatrixRHS;
39 }
40
41 /***********************************************************************************/
42 /***********************************************************************************/
43
44 ✗ void LinearProblemLung::computeElementArrayLungDisplacement(PetscVector& emphysemaCoeff,PetscVector& volReg,PetscVector& volEvol,PetscVector& divVel,PetscMatrix& divBasis, PetscVector& rhsMinus,PetscVector& localExpansion, PetscVector& pCouplingB , PetscVector& volumeB,int currentLabel,PetscVector& prevVel, PetscVector& prevDis, PetscVector& displReal, double currPre,const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint,
45 felInt& iel, FlagMatrixRHS flagMatrixRHS)
46 {
47 (void) emphysemaCoeff;
48 (void) volReg;
49 (void) volEvol;
50 (void) divVel;
51 (void) divBasis;
52 (void) rhsMinus;
53 (void) localExpansion;
54 (void) pCouplingB;
55 (void) volumeB;
56 (void) currentLabel;
57 (void) prevVel;
58 (void) prevDis;
59 (void) displReal;
60 (void) currPre;
61 (void) elemPoint;
62 (void) elemIdPoint;
63 (void) iel;
64 (void) flagMatrixRHS;
65 }
66
67 /***********************************************************************************/
68 /***********************************************************************************/
69
70 ✗ void LinearProblemLung::computeElementArrayLungPressure(PetscVector& emphysemaCoeff, PetscVector& volReg,PetscVector& volEvol,PetscVector& divVel,PetscMatrix& divBasis, PetscVector& rhsMinus,PetscVector& localExpansion,PetscVector& pCouplingB, PetscVector& volumeB ,int currentLabel,PetscVector& prevVel, PetscVector& prevDis, double currPre,const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint,
71 felInt& iel, FlagMatrixRHS flagMatrixRHS)
72 {
73 (void) emphysemaCoeff;
74 (void) volReg;
75 (void) volEvol;
76 (void) divVel;
77 (void) divBasis;
78 (void) rhsMinus;
79 (void) localExpansion;
80 (void) pCouplingB;
81 (void) volumeB;
82 (void) currentLabel;
83 (void) prevVel;
84 (void) prevDis;
85 (void) currPre;
86 (void) elemPoint;
87 (void) elemIdPoint;
88 (void) iel;
89 (void) flagMatrixRHS;
90 }
91
92 /***********************************************************************************/
93 /***********************************************************************************/
94
95 ✗ void LinearProblemLung::computeElementLungVolume(int currentLabel, PetscVector& volReg,PetscVector& displ,const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint,
96 felInt& iel, FlagMatrixRHS flagMatrixRHS)
97 {
98 (void) currentLabel;
99 (void) volReg;
100 (void) displ;
101 (void) elemPoint;
102 (void) elemIdPoint;
103 (void) iel;
104 (void) flagMatrixRHS;
105 }
106
107 /***********************************************************************************/
108 /***********************************************************************************/
109
110 ✗ void LinearProblemLung::assembleMatrixRHSlung(Vec & divVel,Mat & divBasis, PetscVector& rhsMinus,PetscVector& prevVel, PetscVector& prevDis, double currPre, int rank, FlagMatrixRHS flagMatrixRHS)
111 {
112 IGNORE_UNUSED_RANK;
113 ElementType eltType; // geometric element type in the mesh.
114 ✗ int numPointPerElt = 0; // number of points per geometric element.
115 ✗ int currentLabel = 0; // number of label domain.
116 ✗ felInt numEltPerLabel = 0; // number of element for one label and one eltType.
117
118 // Use to define a "global" numbering of element in the mesh.
119 felInt numElement[ GeometricMeshRegion::m_numTypesOfElement ];
120 ✗ for (int ityp=0; ityp<GeometricMeshRegion::m_numTypesOfElement; ityp++ ) {
121 ✗ eltType = (ElementType)ityp;
122 ✗ numElement[eltType] = 0;
123 }
124
125 // contains points of the current element.
126 ✗ std::vector<Point*> elemPoint;
127
128 // contains ids of point of the current element.
129 ✗ std::vector<felInt> elemIdPoint;
130
131 // use to get element number in support dof mesh ordering.
132 felInt ielSupportDof;
133
134 // Assembly loop.
135 // First loop on geometric type.
136 ✗ const std::vector<ElementType>& bagElementTypeDomain = m_meshLocal[m_currentMesh]->bagElementTypeDomain();
137 ✗ for (std::size_t i = 0; i < bagElementTypeDomain.size(); ++i) {
138 ✗ eltType = bagElementTypeDomain[i];
139
140 // Resize array.
141 ✗ numPointPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
142 ✗ elemPoint.resize(numPointPerElt, nullptr);
143 ✗ elemIdPoint.resize(numPointPerElt, 0);
144
145 // Define all current finite element use in the problem from data
146 // file cnfiguration and allocate 1 and elemVec (question: virtual ?).
147 ✗ defineFiniteElement(eltType);
148 ✗ initElementArray();
149
150 // allocate array use for assemble with petsc.
151 ✗ allocateArrayForAssembleMatrixRHS(flagMatrixRHS);
152
153 // Virtual function use in derived problem to allocate elemenField necessary.
154 ✗ initPerElementType(eltType, flagMatrixRHS);
155
156 // Used by user to add specific term (source term for example with elemField).
157 ✗ userElementInit();
158
159 // Second loop on region of the mesh.
160 ✗ for(auto itRef = m_meshLocal[m_currentMesh]->intRefToBegEndMaps[eltType].begin(); itRef != m_meshLocal[m_currentMesh]->intRefToBegEndMaps[eltType].end(); itRef++) {
161 ✗ currentLabel = itRef->first;
162 ✗ numEltPerLabel = itRef->second.second;
163
164 // By default this virtual defines variable m_currentLabel: (m_currentLabel=label).
165 // We can switch on label region with that and define some parameters.
166 ✗ initPerDomain(currentLabel, flagMatrixRHS);
167
168 // Third loop on element in the region with the type: eltType. ("real" loop on elements).
169 ✗ for ( felInt iel = 0; iel < numEltPerLabel; iel++) {
170 ✗ FEL_ASSERT(!m_elementVector.empty());
171
172 // if(!FelisceParam::instance(this->instanceIndex()).duplicateSupportDof) {
173 // Return each id of point of the element and coordinate in two arrays: elemPoint and elemIdPoint.
174 ✗ setElemPoint(eltType, numElement[eltType], elemPoint, elemIdPoint, &ielSupportDof);
175
176 // Clear elementary matrices
177 ✗ for (std::size_t j = 0, size = m_matrices.size(); j < size ; j++)
178 ✗ m_elementMat[j]->zero();
179
180 // Clear elementary std::vector.
181 ✗ for (std::size_t j = 0, size = m_vectors.size(); j < size ; j++)
182 ✗ m_elementVector[j]->zero();
183
184 // Function call in derived problem to compute specific operators of the problem (Heat, N-S,...).
185 ✗ if (FelisceParam::instance(this->instanceIndex()).NSequationFlag == 1 && FelisceParam::instance(this->instanceIndex()).characteristicMethod != 0) {
186 // Use method of characteristics for N-S
187 ✗ computeElementArrayCharact(elemPoint, elemIdPoint, ielSupportDof, eltType, numElement[eltType], flagMatrixRHS);
188 } else {
189 ✗ computeElementArrayLung(divVel,divBasis,rhsMinus,currentLabel,prevVel, prevDis, currPre, elemPoint, elemIdPoint, ielSupportDof, flagMatrixRHS);
190 }
191
192 // Compute specific term of users.
193 ✗ userElementCompute(elemPoint, elemIdPoint, ielSupportDof);
194
195 // Add values of elemMat in the global matrix: m_matrices[0].
196 ✗ setValueMatrixRHS(ielSupportDof, ielSupportDof, flagMatrixRHS);
197 // } else {
198 // m_duplicateSupportDofAssemblyLoop(rank, eltType, numElement[eltType], elemPoint, elemIdPoint, flagMatrixRHS);
199 // }
200
201 ✗ numElement[eltType]++;
202 }
203 }
204 // desallocate array use for assemble with petsc.
205 ✗ desallocateArrayForAssembleMatrixRHS(flagMatrixRHS);
206 }
207
208 // Assembly loop for LumpedModelBC in case of NS model with implicit implementation
209 ✗ if (FelisceParam::instance(this->instanceIndex()).lumpedModelBCLabel.size()>0 && FelisceParam::instance(this->instanceIndex()).model == "NS") {
210 ✗ for (std::size_t iLabel = 0; iLabel < FelisceParam::instance(this->instanceIndex()).lumpedModelBCLabel.size(); iLabel++) {
211 ✗ if(FelisceParam::instance(this->instanceIndex()).lumpedModelBCAlgo[iLabel] == 2) // use a Enum here... 08/13 VM
212 ✗ assembleMatrixImplLumpedModelBC(rank, iLabel);
213 }
214 }
215
216 ✗ if( (flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs) || (flagMatrixRHS == FlagMatrixRHS::only_matrix)) {
217 // real assembling of the matrix manage by PETsC.
218 // TODO Jeremie 01 02 2012: assemble multiple matrix
219 ✗ for (std::size_t i = 0; i < m_matrices.size(); i++)
220 ✗ m_matrices[i].assembly(MAT_FINAL_ASSEMBLY);
221 }
222
223 // call with high level of verbose to print matrix in matlab format.
224 ✗ writeMatrixForMatlab(m_verbosity);
225 ✗ if( (flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs) || (flagMatrixRHS == FlagMatrixRHS::only_rhs)) {
226 // real assembling of the right hand side (RHS).
227 ✗ for (std::size_t i = 0; i < m_vectors.size(); i++) {
228 ✗ m_vectors[i].assembly();
229 }
230 }
231 // call with high level of verbose to print right hand side in matlab format.
232 ✗ writeRHSForMatlab(m_verbosity);
233 }
234
235 /***********************************************************************************/
236 /***********************************************************************************/
237
238 ✗ void LinearProblemLung::assembleMatrixRHSlungDisplacement(PetscVector& emphysemaCoeff, PetscVector& volReg,PetscVector& volEvol,PetscVector& divVel,PetscMatrix& divBasis, PetscVector& rhsMinus, PetscVector& localExpansion, PetscVector& pCouplingB, PetscVector& volumeB,PetscVector& prevVel, PetscVector& prevDis, PetscVector& displReal, double currPre, int rank, FlagMatrixRHS flagMatrixRHS)
239 {
240 IGNORE_UNUSED_RANK;
241 ElementType eltType; // geometric element type in the mesh.
242 ✗ int numPointPerElt = 0; // number of points per geometric element.
243 ✗ int currentLabel = 0; // number of label domain.
244 ✗ felInt numEltPerLabel = 0; // number of element for one label and one eltType.
245
246 // Use to define a "global" numbering of element in the mesh.
247 felInt numElement[ GeometricMeshRegion::m_numTypesOfElement ];
248 ✗ for (int ityp=0; ityp<GeometricMeshRegion::m_numTypesOfElement; ityp++ ) {
249 ✗ eltType = (ElementType)ityp;
250 ✗ numElement[eltType] = 0;
251 }
252
253 // contains points of the current element.
254 ✗ std::vector<Point*> elemPoint;
255
256 // contains ids of point of the current element.
257 ✗ std::vector<felInt> elemIdPoint;
258
259 // use to get element number in support dof mesh ordering.
260 felInt ielSupportDof;
261
262 // Assembly loop.
263 // First loop on geometric type.
264 ✗ const std::vector<ElementType>& bagElementTypeDomain = m_meshLocal[m_currentMesh]->bagElementTypeDomain();
265 ✗ for (std::size_t i = 0; i < bagElementTypeDomain.size(); ++i) {
266 ✗ eltType = bagElementTypeDomain[i];
267
268 // resize array.
269 ✗ numPointPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
270 ✗ elemPoint.resize(numPointPerElt, nullptr);
271 ✗ elemIdPoint.resize(numPointPerElt, 0);
272
273 // define all current finite element use in the problem from data
274 // file cnfiguration and allocate 1 and elemVec (question: virtual ?).
275 ✗ defineFiniteElement(eltType);
276 ✗ initElementArray();
277
278 // allocate array use for assemble with petsc.
279 ✗ allocateArrayForAssembleMatrixRHS(flagMatrixRHS);
280
281 // virtual function use in derived problem to allocate elemenField necessary.
282 ✗ initPerElementType(eltType, flagMatrixRHS);
283
284 // used by user to add specific term (source term for example with elemField).
285 ✗ userElementInit();
286
287 // second loop on region of the mesh.
288 ✗ for(auto itRef = m_meshLocal[m_currentMesh]->intRefToBegEndMaps[eltType].begin(); itRef != m_meshLocal[m_currentMesh]->intRefToBegEndMaps[eltType].end(); itRef++) {
289 ✗ currentLabel = itRef->first;
290 ✗ numEltPerLabel = itRef->second.second;
291
292 // By default this virtual defines variable m_currentLabel: (m_currentLabel=label).
293 // We can switch on label region with that and define some parameters.
294 ✗ initPerDomain(currentLabel, flagMatrixRHS);
295
296 // Third loop on element in the region with the type: eltType. ("real" loop on elements).
297 ✗ for ( felInt iel = 0; iel < numEltPerLabel; iel++) {
298 ✗ FEL_ASSERT(!m_elementVector.empty());
299
300 // if(!FelisceParam::instance(this->instanceIndex()).duplicateSupportDof) {
301 // return each id of point of the element and coordinate in two arrays: elemPoint and elemIdPoint.
302 ✗ setElemPoint(eltType, numElement[eltType], elemPoint, elemIdPoint, &ielSupportDof);
303
304 // clear elementary matrices
305 ✗ for (std::size_t j = 0, size = m_matrices.size(); j < size ; j++)
306 ✗ m_elementMat[j]->zero();
307
308 // clear elementary std::vector.
309 ✗ for (std::size_t j = 0, size = m_vectors.size(); j < size ; j++)
310 ✗ m_elementVector[j]->zero();
311
312 // function call in derived problem to compute specific operators of the problem (Heat, N-S,...).
313 ✗ if (FelisceParam::instance(this->instanceIndex()).NSequationFlag == 1 && FelisceParam::instance(this->instanceIndex()).characteristicMethod != 0) {
314 // use method of characteristics for N-S
315 ✗ computeElementArrayCharact(elemPoint, elemIdPoint, ielSupportDof, eltType, numElement[eltType], flagMatrixRHS);
316 } else {
317 ✗ computeElementArrayLungDisplacement(emphysemaCoeff, volReg, volEvol,divVel,divBasis,rhsMinus,localExpansion,pCouplingB, volumeB,currentLabel,prevVel, prevDis,displReal, currPre, elemPoint, elemIdPoint, ielSupportDof, flagMatrixRHS);
318 }
319
320 // compute specific term of users.
321 ✗ userElementCompute(elemPoint, elemIdPoint, ielSupportDof);
322
323 // add values of elemMat in the global matrix: m_matrices[0].
324 ✗ setValueMatrixRHS(ielSupportDof, ielSupportDof, flagMatrixRHS);
325 // } else {
326 // m_duplicateSupportDofAssemblyLoop(rank, eltType, numElement[eltType], elemPoint, elemIdPoint, flagMatrixRHS);
327 // }
328
329 ✗ numElement[eltType]++;
330 }
331 }
332 // desallocate array use for assemble with petsc.
333 ✗ desallocateArrayForAssembleMatrixRHS(flagMatrixRHS);
334 }
335
336 // Assembly loop for LumpedModelBC in case of NS model with implicit implementation
337 ✗ if (FelisceParam::instance(this->instanceIndex()).lumpedModelBCLabel.size()>0 && FelisceParam::instance(this->instanceIndex()).model == "NS") {
338 ✗ for (std::size_t iLabel = 0; iLabel < FelisceParam::instance(this->instanceIndex()).lumpedModelBCLabel.size(); iLabel++) {
339 ✗ if(FelisceParam::instance(this->instanceIndex()).lumpedModelBCAlgo[iLabel] == 2) // use a Enum here... 08/13 VM
340 ✗ assembleMatrixImplLumpedModelBC(rank, iLabel);
341 }
342 }
343
344 ✗ if( (flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs) || (flagMatrixRHS == FlagMatrixRHS::only_matrix)) {
345 // real assembling of the matrix manage by PETsC.
346 // todo Jeremie 01 02 2012: assemble multiple matrix
347 ✗ for (std::size_t i = 0; i < m_matrices.size(); i++)
348 ✗ m_matrices[i].assembly(MAT_FINAL_ASSEMBLY);
349 }
350
351 // call with high level of verbose to print matrix in matlab format.
352 ✗ writeMatrixForMatlab(m_verbosity);
353 ✗ if( (flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs) || (flagMatrixRHS == FlagMatrixRHS::only_rhs)) {
354 // real assembling of the right hand side (RHS).
355 ✗ for (std::size_t i = 0; i < m_vectors.size(); i++) {
356 ✗ m_vectors[i].assembly();
357 }
358 }
359 // call with high level of verbose to print right hand side in matlab format.
360 ✗ writeRHSForMatlab(m_verbosity);
361 }
362
363 /***********************************************************************************/
364 /***********************************************************************************/
365
366 ✗ void LinearProblemLung::assembleMatrixRHSlungPressure(PetscVector& emphysemaCoeff,PetscVector& volReg,PetscVector& volEvol,PetscVector& divVel,PetscMatrix& divBasis, PetscVector& rhsMinus,PetscVector& localExpansion, PetscVector& pCouplingB, PetscVector& volumeB,PetscVector& prevVel, PetscVector& prevDis, double currPre, int rank, FlagMatrixRHS flagMatrixRHS)
367 {
368 IGNORE_UNUSED_RANK;
369 ElementType eltType; // geometric element type in the mesh.
370 ✗ int numPointPerElt = 0; // number of points per geometric element.
371 ✗ int currentLabel = 0; // number of label domain.
372 ✗ felInt numEltPerLabel = 0; // number of element for one label and one eltType.
373
374 // use to define a "global" numbering of element in the mesh.
375 felInt numElement[ GeometricMeshRegion::m_numTypesOfElement ];
376 ✗ for (int ityp=0; ityp<GeometricMeshRegion::m_numTypesOfElement; ityp++ ) {
377 ✗ eltType = (ElementType)ityp;
378 ✗ numElement[eltType] = 0;
379 }
380
381 // contains points of the current element.
382 ✗ std::vector<Point*> elemPoint;
383
384 // contains ids of point of the current element.
385 ✗ std::vector<felInt> elemIdPoint;
386
387 // use to get element number in support dof mesh ordering.
388 felInt ielSupportDof;
389
390 // Assembly loop.
391 // First loop on geometric type.
392 ✗ const std::vector<ElementType>& bagElementTypeDomain = m_meshLocal[m_currentMesh]->bagElementTypeDomain();
393 ✗ for (std::size_t i = 0; i < bagElementTypeDomain.size(); ++i) {
394 ✗ eltType = bagElementTypeDomain[i];
395
396 // resize array.
397 ✗ numPointPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
398 ✗ elemPoint.resize(numPointPerElt, nullptr);
399 ✗ elemIdPoint.resize(numPointPerElt, 0);
400
401 // define all current finite element use in the problem from data
402 // file cnfiguration and allocate 1 and elemVec (question: virtual ?).
403 ✗ defineFiniteElement(eltType);
404 ✗ initElementArray();
405
406 // allocate array use for assemble with petsc.
407 ✗ allocateArrayForAssembleMatrixRHS(flagMatrixRHS);
408
409 // virtual function use in derived problem to allocate elemenField necessary.
410 ✗ initPerElementType(eltType, flagMatrixRHS);
411
412 // used by user to add specific term (source term for example with elemField).
413 ✗ userElementInit();
414
415 // second loop on region of the mesh.
416 ✗ for(auto itRef = m_meshLocal[m_currentMesh]->intRefToBegEndMaps[eltType].begin(); itRef != m_meshLocal[m_currentMesh]->intRefToBegEndMaps[eltType].end(); itRef++) {
417 ✗ currentLabel = itRef->first;
418 ✗ numEltPerLabel = itRef->second.second;
419
420 // By default this virtual defines variable m_currentLabel: (m_currentLabel=label).
421 // We can switch on label region with that and define some parameters.
422 ✗ initPerDomain(currentLabel, flagMatrixRHS);
423
424 // Third loop on element in the region with the type: eltType. ("real" loop on elements).
425 ✗ for ( felInt iel = 0; iel < numEltPerLabel; iel++) {
426 ✗ FEL_ASSERT(!m_elementVector.empty());
427
428 // if(!FelisceParam::instance(this->instanceIndex()).duplicateSupportDof) {
429 // return each id of point of the element and coordinate in two arrays: elemPoint and elemIdPoint.
430 ✗ setElemPoint(eltType, numElement[eltType], elemPoint, elemIdPoint, &ielSupportDof);
431
432 // clear elementary matrices
433 ✗ for (std::size_t j = 0, size = m_matrices.size(); j < size ; j++)
434 ✗ m_elementMat[j]->zero();
435
436 // clear elementary std::vector.
437 ✗ for (std::size_t j = 0, size = m_vectors.size(); j < size ; j++)
438 ✗ m_elementVector[j]->zero();
439
440 // function call in derived problem to compute specific operators of the problem (Heat, N-S,...).
441 ✗ if (FelisceParam::instance(this->instanceIndex()).NSequationFlag == 1 && FelisceParam::instance(this->instanceIndex()).characteristicMethod != 0) {
442 // use method of characteristics for N-S
443 ✗ computeElementArrayCharact(elemPoint, elemIdPoint, ielSupportDof, eltType, numElement[eltType], flagMatrixRHS);
444 } else
445 ✗ computeElementArrayLungPressure(emphysemaCoeff, volReg, volEvol,divVel,divBasis,rhsMinus,localExpansion,pCouplingB,volumeB,currentLabel,prevVel, prevDis, currPre, elemPoint, elemIdPoint, ielSupportDof, flagMatrixRHS);
446
447 // compute specific term of users.
448 ✗ userElementCompute(elemPoint, elemIdPoint, ielSupportDof);
449
450 // add values of elemMat in the global matrix: m_matrices[0].
451 ✗ setValueMatrixRHS(ielSupportDof, ielSupportDof, flagMatrixRHS);
452 // } else {
453 // m_duplicateSupportDofAssemblyLoop(rank, eltType, numElement[eltType], elemPoint, elemIdPoint, flagMatrixRHS);
454 // }
455
456 ✗ numElement[eltType]++;
457 }
458 }
459 // desallocate array use for assemble with petsc.
460 ✗ desallocateArrayForAssembleMatrixRHS(flagMatrixRHS);
461 }
462
463 // Assembly loop for LumpedModelBC in case of NS model with implicit implementation
464 ✗ if (FelisceParam::instance(this->instanceIndex()).lumpedModelBCLabel.size()>0 && FelisceParam::instance(this->instanceIndex()).model == "NS") {
465 ✗ for (std::size_t iLabel = 0; iLabel < FelisceParam::instance(this->instanceIndex()).lumpedModelBCLabel.size(); iLabel++) {
466 ✗ if(FelisceParam::instance(this->instanceIndex()).lumpedModelBCAlgo[iLabel] == 2) // use a Enum here... 08/13 VM
467 ✗ assembleMatrixImplLumpedModelBC(rank, iLabel);
468 }
469 }
470
471 ✗ if( (flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs) || (flagMatrixRHS == FlagMatrixRHS::only_matrix)) {
472 // real assembling of the matrix manage by PETsC.
473 // todo Jeremie 01 02 2012: assemble multiple matrix
474 ✗ for (std::size_t i = 0; i < m_matrices.size(); i++)
475 ✗ m_matrices[i].assembly(MAT_FINAL_ASSEMBLY);
476 }
477
478 // call with high level of verbose to print matrix in matlab format.
479 ✗ writeMatrixForMatlab(m_verbosity);
480 ✗ if( (flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs) || (flagMatrixRHS == FlagMatrixRHS::only_rhs)) {
481 // real assembling of the right hand side (RHS).
482 ✗ for (std::size_t i = 0; i < m_vectors.size(); i++) {
483 ✗ m_vectors[i].assembly();
484 }
485 }
486 // call with high level of verbose to print right hand side in matlab format.
487 ✗ writeRHSForMatlab(m_verbosity);
488 }
489
490 /***********************************************************************************/
491 /***********************************************************************************/
492
493 ✗ void LinearProblemLung::computeVolumeLungRegions(PetscVector& volReg, PetscVector& displ, int rank, FlagMatrixRHS flagMatrixRHS)
494 {
495 IGNORE_UNUSED_RANK;
496 ElementType eltType; // geometric element type in the mesh.
497 ✗ int numPointPerElt = 0; // number of points per geometric element.
498 ✗ int currentLabel = 0; // number of label domain.
499 ✗ felInt numEltPerLabel = 0; // number of element for one label and one eltType.
500
501 // use to define a "global" numbering of element in the mesh.
502 felInt numElement[ GeometricMeshRegion::m_numTypesOfElement ];
503 ✗ for (int ityp=0; ityp<GeometricMeshRegion::m_numTypesOfElement; ityp++ ) {
504 ✗ eltType = (ElementType)ityp;
505 ✗ numElement[eltType] = 0;
506 }
507
508 // contains points of the current element.
509 ✗ std::vector<Point*> elemPoint;
510
511 // contains ids of point of the current element.
512 ✗ std::vector<felInt> elemIdPoint;
513
514 // use to get element number in support dof mesh ordering.
515 felInt ielSupportDof;
516
517 // allocateArrayForAssembleMatrixRHS(flagMatrixRHS);
518
519 // Assembly loop.
520 // First loop on geometric type.
521 ✗ const std::vector<ElementType>& bagElementTypeDomain = m_meshLocal[m_currentMesh]->bagElementTypeDomain();
522 ✗ for (std::size_t i = 0; i < bagElementTypeDomain.size(); ++i) {
523 ✗ eltType = bagElementTypeDomain[i];
524
525 // resize array.
526 ✗ numPointPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
527 ✗ elemPoint.resize(numPointPerElt, nullptr);
528 ✗ elemIdPoint.resize(numPointPerElt, 0);
529
530 // define all current finite element use in the problem from data
531 // file cnfiguration and allocate 1 and elemVec (question: virtual ?).
532 // defineFiniteElement(eltType);
533 // initElementArray();
534
535 // second loop on region of the mesh.
536 ✗ for(auto itRef = m_meshLocal[m_currentMesh]->intRefToBegEndMaps[eltType].begin(); itRef != m_meshLocal[m_currentMesh]->intRefToBegEndMaps[eltType].end(); itRef++) {
537 ✗ currentLabel = itRef->first;
538 ✗ numEltPerLabel = itRef->second.second;
539
540
541 // Third loop on element in the region with the type: eltType. ("real" loop on elements).
542 ✗ for ( felInt iel = 0; iel < numEltPerLabel; iel++) {
543 ✗ FEL_ASSERT(!m_elementVector.empty());
544
545 ✗ setElemPoint(eltType, numElement[eltType], elemPoint, elemIdPoint, &ielSupportDof);
546 //
547 ✗ computeElementLungVolume(currentLabel, volReg, displ, elemPoint, elemIdPoint, ielSupportDof, flagMatrixRHS);
548
549
550 ✗ numElement[eltType]++;
551 }
552 }
553 }
554 // desallocateArrayForAssembleMatrixRHS(flagMatrixRHS);
555 }
556
557 /***********************************************************************************/
558 /***********************************************************************************/
559
560 ✗ void LinearProblemLung::computeVolumeLungRegions2(PetscVector& volReg,PetscVector& displ,int rank, FlagMatrixRHS flagMatrixRHS)
561 {
562 IGNORE_UNUSED_RANK;
563 ElementType eltType; // geometric element type in the mesh.
564 ✗ int numPointPerElt = 0; // number of points per geometric element.
565 ✗ int currentLabel = 0; // number of label domain.
566 ✗ felInt numEltPerLabel = 0; // number of element for one label and one eltType.
567
568 // use to define a "global" numbering of element in the mesh.
569 felInt numElement[ GeometricMeshRegion::m_numTypesOfElement ];
570 ✗ for (int ityp=0; ityp<GeometricMeshRegion::m_numTypesOfElement; ityp++ ) {
571 ✗ eltType = (ElementType)ityp;
572 ✗ numElement[eltType] = 0;
573 }
574
575 // contains points of the current element.
576 ✗ std::vector<Point*> elemPoint;
577
578 // contains ids of point of the current element.
579 ✗ std::vector<felInt> elemIdPoint;
580
581 // use to get element number in support dof mesh ordering.
582 felInt ielSupportDof;
583
584 // Assembly loop.
585 // First loop on geometric type.
586 ✗ const std::vector<ElementType>& bagElementTypeDomain = m_meshLocal[m_currentMesh]->bagElementTypeDomain();
587 ✗ for (std::size_t i = 0; i < bagElementTypeDomain.size(); ++i) {
588 ✗ eltType = bagElementTypeDomain[i];
589
590 // resize array.
591 ✗ numPointPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
592 ✗ elemPoint.resize(numPointPerElt, nullptr);
593 ✗ elemIdPoint.resize(numPointPerElt, 0);
594
595 // define all current finite element use in the problem from data
596 // file cnfiguration and allocate 1 and elemVec (question: virtual ?).
597 ✗ defineFiniteElement(eltType);
598 ✗ initElementArray();
599
600 // allocate array use for assemble with petsc.
601 ✗ allocateArrayForAssembleMatrixRHS(flagMatrixRHS);
602
603 // virtual function use in derived problem to allocate elemenField necessary.
604 ✗ initPerElementType(eltType, flagMatrixRHS);
605
606 // used by user to add specific term (source term for example with elemField).
607 ✗ userElementInit();
608
609 // second loop on region of the mesh.
610 ✗ for(auto itRef = m_meshLocal[m_currentMesh]->intRefToBegEndMaps[eltType].begin(); itRef != m_meshLocal[m_currentMesh]->intRefToBegEndMaps[eltType].end(); itRef++) {
611 ✗ currentLabel = itRef->first;
612 ✗ numEltPerLabel = itRef->second.second;
613
614
615 // By default this virtual defines variable m_currentLabel: (m_currentLabel=label).
616 // We can switch on label region with that and define some parameters.
617 ✗ initPerDomain(currentLabel, flagMatrixRHS);
618
619 // Third loop on element in the region with the type: eltType. ("real" loop on elements).
620 ✗ for ( felInt iel = 0; iel < numEltPerLabel; iel++) {
621 ✗ FEL_ASSERT(!m_elementVector.empty());
622
623 // if(!FelisceParam::instance(this->instanceIndex()).duplicateSupportDof) {
624 // return each id of point of the element and coordinate in two arrays: elemPoint and elemIdPoint.
625 ✗ setElemPoint(eltType, numElement[eltType], elemPoint, elemIdPoint, &ielSupportDof);
626
627 // clear elementary matrices
628 ✗ for (std::size_t j = 0, size = m_matrices.size(); j < size ; j++)
629 ✗ m_elementMat[j]->zero();
630
631 // clear elementary std::vector.
632 ✗ for (std::size_t j = 0, size = m_vectors.size(); j < size ; j++)
633 ✗ m_elementVector[j]->zero();
634
635 ✗ computeElementLungVolume(currentLabel, volReg, displ, elemPoint, elemIdPoint, ielSupportDof, flagMatrixRHS);
636
637 // compute specific term of users.
638 ✗ userElementCompute(elemPoint, elemIdPoint, ielSupportDof);
639
640 // add values of elemMat in the global matrix: m_matrices[0].
641 ✗ setValueMatrixRHS(ielSupportDof, ielSupportDof, flagMatrixRHS);
642 // } else {
643 // m_duplicateSupportDofAssemblyLoop(rank, eltType, numElement[eltType], elemPoint, elemIdPoint, flagMatrixRHS);
644 // }
645
646 ✗ numElement[eltType]++;
647 }
648 }
649 // desallocate array use for assemble with petsc.
650 ✗ desallocateArrayForAssembleMatrixRHS(flagMatrixRHS);
651 }
652
653 // Assembly loop for LumpedModelBC in case of NS model with implicit implementation
654 ✗ if (FelisceParam::instance(this->instanceIndex()).lumpedModelBCLabel.size()>0 && FelisceParam::instance(this->instanceIndex()).model == "NS") {
655 ✗ for (std::size_t iLabel = 0; iLabel < FelisceParam::instance(this->instanceIndex()).lumpedModelBCLabel.size(); iLabel++) {
656 ✗ if(FelisceParam::instance(this->instanceIndex()).lumpedModelBCAlgo[iLabel] == 2) // use a Enum here... 08/13 VM
657 ✗ assembleMatrixImplLumpedModelBC(rank, iLabel);
658 }
659 }
660
661 ✗ if(flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs || flagMatrixRHS == FlagMatrixRHS::only_matrix) {
662 // real assembling of the matrix manage by PETsC.
663 // todo Jeremie 01 02 2012: assemble multiple matrix
664 ✗ for (std::size_t i = 0; i < m_matrices.size(); i++)
665 ✗ m_matrices[i].assembly(MAT_FINAL_ASSEMBLY);
666 }
667
668 // call with high level of verbose to print matrix in matlab format.
669 ✗ writeMatrixForMatlab(m_verbosity);
670 ✗ if(flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs || flagMatrixRHS == FlagMatrixRHS::only_rhs) {
671 // real assembling of the right hand side (RHS).
672 ✗ for (std::size_t i = 0; i < m_vectors.size(); i++) {
673 ✗ m_vectors[i].assembly();
674 }
675 }
676 // call with high level of verbose to print right hand side in matlab format.
677 ✗ writeRHSForMatlab(m_verbosity);
678 }
679
680 /***********************************************************************************/
681 /***********************************************************************************/
682
683 ✗ void LinearProblemLung::normalLungMesh(PetscVector& modelVec, std::string modeLung, PetscVector& normalVec, felInt idVecVariable){
684
685 ✗ this->initPetscVectors();
686
687 ✗ m_seqVecs.Init("normalField");
688 ✗ m_seqVecs.Init("measStar");
689 ✗ m_vecs.Init("normalField");
690 ✗ m_vecs.Init("measStar");
691
692 ✗ m_vecs.Get("normalField").duplicateFrom(modelVec);
693 ✗ m_vecs.Get("measStar").duplicateFrom(modelVec);
694
695 ✗ gatherVector(m_vecs.Get("normalField"),m_seqVecs.Get("normalField"));
696 ✗ gatherVector(m_vecs.Get("measStar"),m_seqVecs.Get("measStar"));
697
698 ✗ m_vecs.Get("normalField").set(0.0);
699 ✗ m_vecs.Get("measStar").set(0.0);
700
701 ✗ std::cout << "ENTER NORMAL COMPUTATION" << std::endl;
702
703 ✗ std::cout << "size lis bc = " << m_boundaryConditionList.size() << std::endl;
704
705 ✗ std::vector<felInt> m_nLabels;
706
707 ✗ if(modeLung=="displacement"){
708 ✗ for (std::size_t i=0; i < m_boundaryConditionList.numDirichletBoundaryCondition(); i++){
709 ✗ std::cout << "dirichlet i = " << i << std::endl;
710 ✗ BoundaryCondition* BCn = m_boundaryConditionList.Dirichlet(i);
711 ✗ for(auto it_labelNumber = BCn->listLabel().begin(); it_labelNumber != BCn->listLabel().end(); it_labelNumber++) {
712 ✗ m_nLabels.push_back(*it_labelNumber);
713 ✗ std::cout << "het h " << std::endl;
714 }
715 }
716 }
717
718 ✗ if(modeLung=="pressure"){
719 ✗ for (std::size_t i=0; i < m_boundaryConditionList.numNeumannNormalBoundaryCondition(); i++){
720 ✗ std::cout << "neumannNormal i = " << i << std::endl;
721 ✗ BoundaryCondition* BCn = m_boundaryConditionList.NeumannNormal(i);
722 ✗ for(auto it_labelNumber = BCn->listLabel().begin(); it_labelNumber != BCn->listLabel().end(); it_labelNumber++) {
723 ✗ m_nLabels.push_back(*it_labelNumber);
724 ✗ std::cout << "het h " << std::endl;
725 }
726 }
727 }
728
729 ✗ std::cout << "Gets to normal computation " << std::endl;
730
731 ✗ computeNormalField(m_nLabels, idVecVariable);
732
733 ✗ normalVec.duplicateFrom(m_seqVecs.Get("normalField"));
734 ✗ normalVec.copyFrom(m_seqVecs.Get("normalField"));
735 ✗ VecView(normalVec.toPetsc(),PETSC_VIEWER_STDOUT_WORLD);
736 }
737
738 /***********************************************************************************/
739 /***********************************************************************************/
740
741 ✗ void LinearProblemLung::applyBoundaryConditionMecha(double ratioM,PetscVector& volCurr,PetscVector& volFRC,PetscVector& volTLC, double dt, double coeffReaction ,const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint, felInt& iel, FlagMatrixRHS flagMatrixRHS)
742 {
743 IGNORE_UNUSED_IEL;
744 IGNORE_UNUSED_ELEM_POINT;
745 IGNORE_UNUSED_ELEM_ID_POINT;
746 ✗ BoundaryCondition* BC = nullptr;
747 ✗ CurvilinearFiniteElement* fe = nullptr;
748 ✗ int iUnknown = -1;
749 ✗ int iblock = -1;
750 ✗ int cptComp = -1;
751
752 ✗ FEL_ASSERT(!m_elementVectorBD.empty());
753
754 //initialization
755
756 //used for mechanically induced ventilation in the lung model
757 //add surface integral terms for diaphragm and chest rection
758
759 // double cl, cr;
760 double coeff;
761 int nbLobes;
762 ✗ volCurr.getSize(&nbLobes);
763 // PetscScalar Vprev[nbLobes];
764 ✗ PetscScalar Vcurr[nbLobes];
765 ✗ PetscScalar Vtlc[nbLobes];
766 ✗ PetscScalar Vfrc[nbLobes];
767 ✗ PetscInt ix[nbLobes];
768 ✗ for(int i=0;i<nbLobes;i++) {
769 ✗ ix[i]=i;
770 }
771 ✗ volCurr.getValues(nbLobes,ix,Vcurr);
772 ✗ volTLC.getValues(nbLobes,ix,Vtlc);
773 ✗ volFRC.getValues(nbLobes,ix,Vfrc);
774
775 ✗ for (std::size_t iNN = 0; iNN < m_boundaryConditionList.numNeumannNormalBoundaryCondition(); iNN++) {
776 ✗ BC = m_boundaryConditionList.NeumannNormal(iNN);
777 ✗ int iVariable = m_listVariable.getVariableIdList(BC->getVariable().physicalVariable());
778 ✗ iUnknown = BC->getUnknown();
779 ✗ fe = m_listCurvilinearFiniteElement[iVariable];
780 ✗ FEL_ASSERT(fe);
781 ✗ for(auto it_labelNumber = BC->listLabel().begin();
782 ✗ it_labelNumber != BC->listLabel().end(); it_labelNumber++) {
783 ✗ if(*it_labelNumber == m_currentLabel) {
784
785 ✗ int iL=m_currentLabel-1;
786
787 ✗ double ratioVol=0;
788 ✗ ratioVol=std::abs((Vcurr[iL]*1000000+Vfrc[iL])/Vtlc[iL]);
789 ✗ if(ratioVol<ratioM){
790 ✗ coeff=coeffReaction*(Vtlc[iL]-Vfrc[iL])*(Vcurr[iL]*1000000)/(std::abs(Vtlc[iL]-Vcurr[iL]*1000000-Vfrc[iL])); //conversion of the current volume into milliliters, unit in which Vfrc and Vtlc are provided
791 } else{
792 ✗ const double alph=Vtlc[iL]-Vfrc[iL];
793 ✗ const double Vc=ratioM*Vtlc[iL]-Vfrc[iL];
794 ✗ const double a=std::abs(coeffReaction*alph*alph/((alph-Vc)*(alph-Vc)));
795 ✗ const double b=coeffReaction*alph*Vc/(alph-Vc)-a*Vc;
796 ✗ coeff=a*Vcurr[iL]*1000000+b;
797 }
798
799 // cl=dt*dt*coeff;
800 // cr=-dt*coeff;
801 //dt useless, should be re;oved fro; the argument (Nicolas)
802 //following lines to avoid warning
803 ✗ double valTmp=0;
804 ✗ valTmp=dt;
805 ✗ dt=0;
806 ✗ dt=valTmp;
807
808 ✗ cptComp = 0;
809 ✗ if (flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs){
810 ✗ for(auto it_comp = BC->getComp().begin(); it_comp != BC->getComp().end(); ++it_comp) {
811 ✗ iblock = m_listUnknown.getBlockPosition(iUnknown,*it_comp);
812 ✗ m_elementMatBD[3]->setZerosLung(0,*fe,iblock,iblock,m_listVariable[iVariable].numComponent()); //to ensure the matrix is null before computation
813 ✗ m_elementMatBD[3]->phi_i_phi_j(coeff,*fe,iblock,iblock,m_listVariable[iVariable].numComponent());
814 ✗ cptComp++;
815 }
816 ✗ } else if (flagMatrixRHS == FlagMatrixRHS::only_rhs){
817 ✗ for(auto it_comp = BC->getComp().begin(); it_comp != BC->getComp().end(); ++it_comp) {
818 ✗ iblock = m_listUnknown.getBlockPosition(iUnknown,*it_comp);
819 ✗ m_elementMatBD[3]->setZerosLung(0,*fe,iblock,iblock,m_listVariable[iVariable].numComponent()); //to ensure the matrix is null before computation
820 ✗ m_elementMatBD[3]->phi_i_phi_j(coeff,*fe,iblock,iblock,m_listVariable[iVariable].numComponent());
821 ✗ cptComp++;
822 }
823 }
824 }
825 }
826 }
827 }
828
829 /***********************************************************************************/
830 /***********************************************************************************/
831
832 ✗ void LinearProblemLung::assembleMatrixRHSBoundaryConditionMecha(double ratioM, PetscVector& volCurr,PetscVector& volFRC,PetscVector& volTLC, double dt, double coeffReaction ,int rank, FlagMatrixRHS flagMatrixRHS, const int idBCforLinCombMethod)
833 {
834 IGNORE_UNUSED_RANK;
835 ElementType eltType; //geometric element type in the mesh.
836 ✗ int numPointPerElt = 0; //number of points per geometric element.
837 ✗ int currentLabel = 0; //number of label domain.
838 ✗ felInt numEltPerLabel = 0; //number of element for one label and one eltType.
839 // use to define a "global" numbering of element in the mesh.
840 felInt numElement[ GeometricMeshRegion::m_numTypesOfElement ];
841 ✗ for (int ityp=0; ityp<GeometricMeshRegion::m_numTypesOfElement; ityp++ ) {
842 ✗ eltType = (ElementType)ityp;
843 ✗ numElement[eltType] = 0;
844 }
845
846 // I need to define the CurrentFEWithBd to correctly compute the stress tensor
847 ✗ if(m_meshLocal[m_currentMesh]->statusFaces() == false) {
848 ✗ m_meshLocal[m_currentMesh]->buildFaces();
849 }
850 ✗ const std::vector<ElementType>& bagElementTypeDomain = m_meshLocal[m_currentMesh]->bagElementTypeDomain();
851 ✗ for (std::size_t ielt = 0; ielt < bagElementTypeDomain.size(); ++ielt) {
852 ✗ eltType = bagElementTypeDomain[ielt];
853 ✗ defineCurrentFiniteElementWithBd(eltType);
854 }
855
856 // contains points of the current element.
857 ✗ std::vector<Point*> elemPoint;
858
859 // contains ids of point of the current element.
860 ✗ std::vector<felInt> elemIdPoint;
861
862 // contains ids of all the support elements associated to a mesh element
863 ✗ std::vector<felInt> vectorIdSupport;
864
865 // use to get element number in support dof mesh ordering.
866 felInt ielSupportDof;
867
868 ✗ allocateVectorBoundaryConditionDerivedLinPb();
869
870 //First loop on geometric type.
871 ✗ const std::vector<ElementType>& bagElementTypeDomainBoundary = m_meshLocal[m_currentMesh]->bagElementTypeDomainBoundary();
872 ✗ for (std::size_t i = 0; i < bagElementTypeDomainBoundary.size(); ++i) {
873 ✗ eltType = bagElementTypeDomainBoundary[i];
874 //resize array.
875 ✗ numPointPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
876 ✗ elemPoint.resize(numPointPerElt, nullptr);
877 ✗ elemIdPoint.resize(numPointPerElt, 0);
878
879 //define all current finite element use in the problem from data
880 //file configuration and allocate elemMat and elemVec (question: virtual ?).
881 ✗ defineCurvilinearFiniteElement(eltType);
882 ✗ initElementArrayBD();
883
884 //allocate array use for assemble with petsc.
885 ✗ allocateArrayForAssembleMatrixRHSBD(flagMatrixRHS);
886
887
888 //Function (1), (2), (3) doing the same stuff for different use: allocate/initialize elementField
889 // (1) virtual function use in derived problem to allocate elemenField necessary.
890 ✗ initPerElementTypeBoundaryCondition(eltType, flagMatrixRHS);
891
892 // (2) use by user to add specific term (term source for example with elemField).
893 ✗ userElementInitNaturalBoundaryCondition();
894
895 // (3) allocate Elemfield and setValue if BC = constant in space
896 // if BC = function : setValue in users
897 ✗ allocateElemFieldBoundaryCondition(idBCforLinCombMethod);
898 //second loop on region of the mesh.
899 ✗ for(auto itRef = m_meshLocal[m_currentMesh]->intRefToBegEndMaps[eltType].begin(); itRef != m_meshLocal[m_currentMesh]->intRefToBegEndMaps[eltType].end(); itRef++) {
900 ✗ currentLabel = itRef->first;
901 ✗ numEltPerLabel = itRef->second.second;
902
903 //By default this virtual define variable m_currentLabel: (m_currentLabel=label).
904 //We can switch on label region with that and define some parameters.
905 //We can also fill the elemField allocated in initPerElementTypeBoundaryCondition()
906 ✗ initPerDomainBoundaryCondition(elemPoint, elemIdPoint, currentLabel, numEltPerLabel, &ielSupportDof, eltType, numElement[eltType], flagMatrixRHS); //it updates label
907 //Third loop on element in the region with the type: eltType. ("real" loop on elements).
908 ✗ for ( felInt iel = 0; iel < numEltPerLabel; iel++) {
909 // return each id of point of the element and coordinate in two arrays: elemPoint and elemIdPoint.
910 ✗ setElemPoint(eltType, numElement[eltType], elemPoint, elemIdPoint, vectorIdSupport);
911
912 // loop over all the support elements
913 ✗ for (std::size_t it = 0; it < vectorIdSupport.size(); it++) {
914 // get the id of the support
915 ✗ ielSupportDof = vectorIdSupport[it];
916 // clear elementary matrix.
917 // todo Jeremie 01 02 2012: assemble multiple matrix
918 ✗ for (std::size_t j = 0; j < m_matrices.size(); j++)
919 ✗ m_elementMatBD[j]->zero();
920
921 // clear elementary std::vector.
922 ✗ for (std::size_t j = 0; j < m_vectors.size(); j++)
923 ✗ m_elementVectorBD[j]->zero();
924
925 // function call in derived problem to compute specific operators of the problem (Heat, N-S,...).initialized in (1)
926 ✗ computeElementArrayBoundaryCondition(elemPoint, elemIdPoint, ielSupportDof,flagMatrixRHS); //it updates fe..
927 // compute specific term of users (Neumann transient for example) initialized in (2)
928 ✗ userElementComputeNaturalBoundaryCondition(elemPoint, elemIdPoint, ielSupportDof, currentLabel);
929 // apply bc initialized in (3) if constant or initialized in (1) or (2) if functions
930 ✗ applyBoundaryConditionMecha(ratioM,volCurr, volFRC, volTLC, dt, coeffReaction ,elemPoint, elemIdPoint, ielSupportDof, flagMatrixRHS);
931 // add values of elemMat in the global matrix: m_matrices[0].
932 ✗ setValueMatrixRHSBD(ielSupportDof, ielSupportDof, flagMatrixRHS);
933 }
934 ✗ numElement[eltType]++;
935 }
936 }
937 //allocate array use for assemble with petsc.
938 ✗ desallocateArrayForAssembleMatrixRHS(flagMatrixRHS);
939 }
940 ✗ if( (flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs) || (flagMatrixRHS == FlagMatrixRHS::only_matrix)) {
941 //real assembling of the matrix manage by PETsC.
942 ✗ for (std::size_t i = 0; i < m_matrices.size(); i++)
943 ✗ m_matrices[i].assembly(MAT_FINAL_ASSEMBLY);
944 }
945
946 //useless - to avoid warning linked to the fact prevDis is not used - to be corrected (Nicolas)
947 // VecAYPX(prevDis,0,prevDis);
948 ////////////////////////////////////////////
949
950 //call with high level of verbose to print matrix in matlab format.
951 ✗ writeMatrixForMatlab(m_verbosity);
952 ✗ if( (flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs) || (flagMatrixRHS == FlagMatrixRHS::only_rhs)) {
953 //real assembling of the right hand side (RHS).
954 ✗ for (std::size_t index = 0; index < m_vectors.size(); index++) {
955 ✗ m_vectors[index].assembly();
956 }
957 }
958
959 //call with high level of verbose to print right hand side in matlab format.
960 ✗ writeRHSForMatlab(m_verbosity);
961 }
962
963 /***********************************************************************************/
964 /***********************************************************************************/
965
966 ✗ void LinearProblemLung::applyBCMecha(double ratioM, PetscVector& volCurr,PetscVector& volFRC,PetscVector& volTLC, double dt, double coeffReaction, const int & applicationOption, int rank, FlagMatrixRHS flagMatrixRHSEss, FlagMatrixRHS flagMatrixRHSNat,
967 const int idBCforLinCombMethod, bool doPrintBC,
968 ApplyNaturalBoundaryConditionsType do_apply_natural) {
969 ✗ if (doPrintBC && m_verbosity > 2) {
970 ✗ PetscPrintf(MpiInfo::petscComm(), "\n/============== Information about boundary condition.===============/\n");
971 ✗ m_boundaryConditionList.print(m_verbosity);
972 }
973 // In some dynamic problems we might want to be able to skip natural conditions, hence this flag.
974 // Its role is absolutely not the same as the other flag addedBoundaryFlag: in a same problem
975 // FelisceParam::instance(this->instanceIndex()).addedBoundaryFlag is probably the same all the way whereas do_apply_natural is deemed to be
976 // true in the first time iteration and then false otherwise in some time schemes.
977 ✗ if (do_apply_natural == ApplyNaturalBoundaryConditions::yes) {
978 // Natural boundary condition (Neumann, Neumann Normal, Robin, RobinNormal, EmbedFSI, Natural LumpedModelBC, Backflow stabilization for now).
979 ✗ if (m_boundaryConditionList.NaturalBoundaryCondition()) {
980 // Assemble Matrix and/or RHS (depends on flagMatrixRHSNat)
981 ✗ m_elementMat[3]->zero();
982 ✗ assembleMatrixRHSBoundaryConditionMecha(ratioM, volCurr, volFRC, volTLC, dt, coeffReaction,rank, flagMatrixRHSNat, idBCforLinCombMethod);
983 }
984 }
985
986 // Essential boundary condition (Dirichlet, Essential lumpedModelBC).
987 ✗ if (m_boundaryConditionList.EssentialBoundaryCondition() || FelisceParam::instance(this->instanceIndex()).useEssDerivedBoundaryCondition)
988 ✗ applyEssentialBoundaryCondition(applicationOption, rank, flagMatrixRHSEss);
989 }
990
991 /***********************************************************************************/
992 /***********************************************************************************/
993
994 ✗ void LinearProblemLung::computeNormalField(const std::vector<int> & labels, felInt idVecVariable)
995 {
996 ✗ m_auxiliaryInts.resize(1);
997 ✗ m_auxiliaryInts[0]=idVecVariable;//for instance m_iVelocity;
998 /*!
999 For each dof d of the surface the sum: \f$\sum_{K} m(K)\vec n(K)\f$ is computed.
1000 - K are all the boundary elements that share the dof d.
1001 - m(K) is the measure of the element
1002 - n(K) is the P0 normal defined on the element
1003 */
1004 ✗ this->assemblyLoopBoundaryGeneral(&LinearProblemLung::normalFieldComputer,labels,&LinearProblemLung::initPerETNormVel,&LinearProblemLung::updateFeNormVel);
1005
1006 ✗ m_vecs.Get("normalField").assembly();
1007 ✗ m_vecs.Get("measStar").assembly();
1008
1009 /// Than the vector is divided by \f$\sum_{K} m(K)\f$ to obtain the P1 field.
1010 ✗ pointwiseDivide(m_vecs.Get("normalField"), m_vecs.Get("normalField"), m_vecs.Get("measStar"));
1011 ✗ gatherVector( m_vecs.Get("normalField"), m_seqVecs.Get("normalField") );
1012
1013 /// using normalizeComputer we normalize the result.
1014 ✗ this->assemblyLoopBoundaryGeneral(&LinearProblemLung::normalizeComputer,labels,&LinearProblemLung::initPerETNormVel,&LinearProblemLung::updateFeNormVel);
1015
1016 ✗ m_vecs.Get("normalField").assembly();
1017 ✗ gatherVector( m_vecs.Get("normalField"), m_seqVecs.Get("normalField") );
1018 }
1019
1020 /***********************************************************************************/
1021 /***********************************************************************************/
1022
1023 ✗ void LinearProblemLung::normalFieldComputer(felInt ielSupportDof)
1024 {
1025 double value;
1026 ✗ for ( std::size_t iSurfDof(0); iSurfDof < ( std::size_t ) m_listCurvilinearFiniteElement[m_auxiliaryInts[0]]->numDof(); iSurfDof++) {
1027 ✗ double feMeas(m_listCurvilinearFiniteElement[m_auxiliaryInts[0]]->measure());
1028 ✗ for ( std::size_t iCoor(0); iCoor < ( std::size_t ) m_listCurvilinearFiniteElement[m_auxiliaryInts[0]]->numCoor(); iCoor++) {
1029 felInt iDofGlob;
1030 ✗ dof().loc2glob( ielSupportDof, iSurfDof, m_auxiliaryInts[0], iCoor, iDofGlob);
1031 ✗ AOApplicationToPetsc(m_ao,1,&iDofGlob);
1032 ✗ value=0.0;
1033 ✗ for ( std::size_t iQuadBd(0.0) ; iQuadBd < ( std::size_t )m_listCurvilinearFiniteElement[m_auxiliaryInts[0]]->numQuadraturePoint(); iQuadBd++) {
1034 ✗ value += m_listCurvilinearFiniteElement[m_auxiliaryInts[0]]->weightMeas( iQuadBd ) * m_listCurvilinearFiniteElement[m_auxiliaryInts[0]]->normal[iQuadBd](iCoor);
1035 }
1036 ✗ m_vecs.Get("normalField").setValue(iDofGlob, value, ADD_VALUES);
1037 ✗ m_vecs.Get("measStar").setValue(iDofGlob, feMeas, ADD_VALUES);
1038 }
1039 }
1040 }
1041
1042 /***********************************************************************************/
1043 /***********************************************************************************/
1044
1045 ✗ void LinearProblemLung::updateFeNormVel(const std::vector<Point*>& elemPoint,const std::vector<int>& /*elemIdPoint*/)
1046 {
1047 ✗ m_listCurvilinearFiniteElement[m_auxiliaryInts[0]]->updateMeasNormal(0,elemPoint);
1048 }
1049
1050 /***********************************************************************************/
1051 /***********************************************************************************/
1052
1053 ✗ void LinearProblemLung::normalizeComputer(felInt ielSupportDof)
1054 {
1055 ✗ const int numComp = m_listCurvilinearFiniteElement[m_auxiliaryInts[0]]->numCoor( );
1056 ✗ const int numDof = m_listCurvilinearFiniteElement[m_auxiliaryInts[0]]->numDof( );
1057
1058 ✗ std::vector< double > norm(numDof,0.0);
1059 ✗ std::vector< std::vector< double > > NN(numDof, std::vector< double >(numComp, 0.0));
1060 ✗ for ( int idof = 0; idof < numDof; idof++ ) {
1061 ✗ for ( int icomp = 0; icomp < numComp; icomp++ ) {
1062 int iDofGlob;
1063 double value;
1064 ✗ dof().loc2glob( ielSupportDof, idof, m_auxiliaryInts[0], icomp, iDofGlob);
1065 ✗ AOApplicationToPetsc(m_ao, 1, &iDofGlob);
1066 ✗ m_seqVecs.Get("normalField").getValues(1, &iDofGlob, &value);
1067 ✗ norm[idof] += value*value;
1068 ✗ NN[idof][icomp] = value;
1069 }
1070 ✗ norm[idof] = std::sqrt( norm[idof] );
1071 ✗ for ( int icomp = 0; icomp < numComp; icomp++ )
1072 ✗ NN[idof][icomp] = NN[idof][icomp]/norm[idof];
1073 }
1074 ✗ for ( int idof = 0; idof < numDof; idof++ ) {
1075 ✗ for ( int icomp = 0; icomp < numComp; icomp++ ) {
1076 int iDofGlob;
1077 ✗ dof().loc2glob( ielSupportDof, idof, m_auxiliaryInts[0], icomp, iDofGlob);
1078 ✗ AOApplicationToPetsc(m_ao, 1, &iDofGlob);
1079 ✗ m_vecs.Get("normalField").setValue(iDofGlob, NN[idof][icomp], INSERT_VALUES);
1080 }
1081 }
1082 }
1083
1084 }
1085