GCC Code Coverage Report

Directory: ./
File: Solver/linearProblemHarmonicExtension.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 86 184 46.7%
Branches: 61 260 23.5%
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 "Core/felisceTools.hpp"
21 #include "Solver/linearProblemHarmonicExtension.hpp"
22 #include "InputOutput/io.hpp"
23 #include "DegreeOfFreedom/boundaryCondition.hpp"
24
25 namespace felisce {
26 30 LinearProblemHarmonicExtension::LinearProblemHarmonicExtension():
27 LinearProblem("Poisson equation"),
28
4/8
✓ Branch 2 taken 30 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 30 times.
✗ Branch 6 not taken.
✓ Branch 10 taken 30 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 30 times.
✗ Branch 14 not taken.
 30 m_buildTeporaryMatrix(false)
29 {
30
31 30 }
32
33
34
35 60 LinearProblemHarmonicExtension::~LinearProblemHarmonicExtension() {
36
2/2
✓ Branch 0 taken 24 times.
✓ Branch 1 taken 6 times.
 60 if(m_buildTeporaryMatrix)
37 48 m_matrix.destroy();
38 }
39
40 30 void LinearProblemHarmonicExtension::initialize(std::vector<GeometricMeshRegion::Pointer>& mesh, FelisceTransient::Pointer fstransient, MPI_Comm& comm, bool doUseSNES) {
41
1/2
✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
 30 LinearProblem::initialize(mesh,comm,doUseSNES);
42 30 m_fstransient = fstransient;
43
1/2
✓ Branch 2 taken 30 times.
✗ Branch 3 not taken.
 30 std::vector<PhysicalVariable> listVariable(3);
44
1/2
✓ Branch 2 taken 30 times.
✗ Branch 3 not taken.
 30 std::vector<std::size_t> listNumComp(3);
45
46 30 listVariable[0] = velocity;
47 30 listNumComp[0] = this->dimension();
48
49 30 listVariable[1] = pressure;
50 30 listNumComp[1] = 1;
51
52 30 listVariable[2] = displacement;
53 30 listNumComp[2] = this->dimension();
54
55 //define unknown of the linear system.
56
1/2
✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
 30 m_listUnknown.push_back(displacement);
57
1/2
✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
 30 this->definePhysicalVariable(listVariable,listNumComp);
58
1/2
✓ Branch 2 taken 30 times.
✗ Branch 3 not taken.
 30 m_iDisplacement = this->listVariable().getVariableIdList(displacement);
59
60
1/2
✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
 30 const auto& r_felisce_param_instance = FelisceParam::instance();
61
62 30 m_scaleCoeff = r_felisce_param_instance.scaleCoeff;
63
64
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 26 times.
 30 if(r_felisce_param_instance.stiffening){
65
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 initializeStiffnessVector();
66 }
67 30 }
68
69
70 564 void LinearProblemHarmonicExtension::initPerElementType(ElementType eltType, FlagMatrixRHS flagMatrixRHS) {
71 (void) eltType;
72 (void) flagMatrixRHS;
73 564 m_fe = m_listCurrentFiniteElement[m_iDisplacement];
74 564 }
75
76 5060330 void LinearProblemHarmonicExtension::computeElementArray(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint, felInt& iel, FlagMatrixRHS flagMatrixRHS) {
77 (void) elemIdPoint;
78 (void) elemPoint;
79 (void) iel;
80 (void) flagMatrixRHS;
81
82 5060330 const auto& r_felisce_param_instance = FelisceParam::instance();
83
84
2/2
✓ Branch 0 taken 79130 times.
✓ Branch 1 taken 4981200 times.
 5060330 if ( r_felisce_param_instance.useElasticExtension ) {
85 79130 m_fe->updateFirstDeriv(0, elemPoint);
86 // Elastic extension (à la Tezduyar).
87 // Note: ( coef = 1/J -> to stiff small elements )
88
89 79130 double mu = r_felisce_param_instance.mu_lame;
90 79130 double lambda = r_felisce_param_instance.lambda_lame;
91
92
6/6
✓ Branch 1 taken 43130 times.
✓ Branch 2 taken 36000 times.
✓ Branch 4 taken 22290 times.
✓ Branch 5 taken 20840 times.
✓ Branch 6 taken 22290 times.
✓ Branch 7 taken 56840 times.
 79130 if(stiffenElementIds.size() > 0 && stiffenElementIds[iel] == 1){
93
94
2/2
✓ Branch 1 taken 2229 times.
✓ Branch 2 taken 20061 times.
 22290 if(m_fstransient->iteration == 1){
95
96
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2229 times.
 2229 if(r_felisce_param_instance.skipStiffeningFirstIteration == 0) {
97 ✗ mu = mu /m_fe->measure() * r_felisce_param_instance.initStiff;
98 ✗ lambda = lambda /m_fe->measure() * r_felisce_param_instance.initStiff;
99 } else {
100 2229 mu = mu / m_fe->measure();
101 2229 lambda = lambda / m_fe->measure();
102 }
103 } else {
104 20061 mu = mu /m_fe->measure() * r_felisce_param_instance.finalStiff;
105 20061 lambda = lambda /m_fe->measure() * r_felisce_param_instance.finalStiff;
106 }
107 } else{
108 56840 mu = mu / m_fe->measure();
109 56840 lambda = lambda / m_fe->measure();
110 }
111
112 79130 m_elementMat[0]->eps_phi_i_eps_phi_j(2*mu, *m_fe, 0, 0, this->dimension());
113 79130 m_elementMat[0]->div_phi_j_div_phi_i(lambda, *m_fe, 0, 0, this->dimension());
114
115
2/2
✓ Branch 1 taken 251200 times.
✓ Branch 2 taken 4730000 times.
 4981200 } else if ( m_fstransient->iteration == 0 ) { // Check if this is really called? Because iteration begins at 1 for HarmonicExtension. (This-Boilevin, 2017).
116 251200 m_fe->updateFirstDeriv(0, elemPoint);
117 // Harmonic extension
118 251200 m_elementMat[0]->grad_phi_i_grad_phi_j(1., *m_fe, 0, 0, this->dimension());
119 }
120 5060330 }
121
122 24 void LinearProblemHarmonicExtension::copyMatrixRHS() {
123 24 m_matrix.duplicateFrom(matrix(0),MAT_COPY_VALUES);
124 24 m_matrix.assembly(MAT_FINAL_ASSEMBLY);
125 24 m_buildTeporaryMatrix = true;
126 24 }
127
128 540 void LinearProblemHarmonicExtension::addMatrixRHS() {
129
2/2
✓ Branch 1 taken 470 times.
✓ Branch 2 taken 70 times.
 540 if ( ! FelisceParam::instance().useElasticExtension )
130 470 matrix(0).axpy(1,m_matrix,SAME_NONZERO_PATTERN);
131 540 }
132
133 ✗ void LinearProblemHarmonicExtension::userFinalizeEssBCTransient() {}
134
135
136 //--------------------------------------------------------------------------------------
137 ✗ void LinearProblemHarmonicExtension::readMatchFile_ALE(ListVariable , const std::string& fileName) {
138 ✗ m_numDofBySupportToMatch = 0;
139 //Read match file.
140 FILE *pFile;
141 ✗ pFile = fopen(fileName.c_str(),"r");
142 ✗ if (pFile==nullptr) {
143 ✗ FEL_WARNING("Impossible to read match file: "+fileName+".");
144 } else {
145 //Read supportDof in match file.
146 ✗ if (! fscanf(pFile, "%d", &m_totNumNodeToMatch) )
147 ✗ FEL_ERROR("totalNumberOfNode in interface not scanned in file.");
148 ✗ if (!m_totNumNodeToMatch) {
149 ✗ FEL_ERROR("Your match file is empty.\n");
150 }
151 ✗ m_listOfNodeToMatch = new felInt[m_totNumNodeToMatch];
152 ✗ felInt id = 0;
153 ✗ for (felInt indc = 0; indc < m_totNumNodeToMatch; indc++) {
154 ✗ if (! fscanf(pFile, "%d", &id) )
155 ✗ FEL_ERROR("id of interface supportDof not scanned in file.");
156 ✗ m_listOfNodeToMatch[indc] = id - 1;
157 }
158
159 ✗ m_numDofBySupportToMatch = this->listVariable()[m_iDisplacement].numComponent();
160 ✗ m_totNumDofToMatch = m_numDofBySupportToMatch * m_totNumNodeToMatch;
161 ✗ m_vectorDisp.resize(m_totNumDofToMatch, 0.);
162 ✗ m_listOfDofToMatch.resize(m_totNumDofToMatch, 0.);
163
164 }
165 ✗ fclose(pFile);
166 }
167
168 ✗ void LinearProblemHarmonicExtension::identifyIdDofToMatch_ALE(Dof&)
169 {
170 ✗ std::vector<felInt> idDofBySupport;
171
172 ✗ for(felInt isd = 0; isd<m_totNumNodeToMatch; isd++) {
173 ✗ dof().supportDofToDof(m_listOfNodeToMatch[isd], idDofBySupport);
174
175 ✗ for (std::size_t iComp = 0; iComp < this->listVariable()[m_iDisplacement].numComponent(); iComp++) {
176 ✗ m_listOfDofToMatch[ isd * this->listVariable()[m_iDisplacement].numComponent() + iComp ] = idDofBySupport[ iComp ];
177 }
178 ✗ idDofBySupport.clear();
179 }
180 }
181
182
183 //Author : Alexandre This
184 //Date : 21 July 2016
185 //This function could/should maybe be renamed. Its purpose is not only to read displacement but is also used to read velocity.
186 //The use of this function to read velocity predates this comment and the accompanying commit.
187 //Note also that the function is not only reading the data but is applying a transform on it (the simplest one being a scale factor applied to the motion)
188
189 // Modification by Ludovic Boilevin-Kayl/Alexandre This to use the correct numerical scheme used by M3DISIM (05/2017)
190 ✗ void LinearProblemHarmonicExtension::readDataDisplacement(std::vector<IO::Pointer>& io, double ReadTime)
191 {
192 felInt idHE;
193
194 ✗ if(m_fstransient->iteration==1){
195 ✗ m_vectorDisp.resize(m_totNumDofToMatch, 0.);
196 ✗ m_oldVectorDisp.resize(m_totNumDofToMatch, 0.);
197 ✗ m_oldOldVectorDisp.resize(m_totNumDofToMatch, 0.);
198 }
199 //We need an initial displacement before being able to use the velocity
200 ✗ if(m_fstransient->iteration == 1 or FelisceParam::instance().updateMeshByVelocity == false ) {
201 ✗ idHE = 1; // displacement
202 }
203 else {
204 ✗ idHE = 0; // velocity
205 }
206
207 ✗ const int iMesh = m_listVariable[idHE].idMesh();
208
209 ✗ if(this->m_mesh[m_currentMesh]->domainDim() == 3) {
210 ✗ io[iMesh]->readVariable(idHE, ReadTime, &m_vectorDisp[0], m_vectorDisp.size());
211 ✗ } else if (this->m_mesh[m_currentMesh]->domainDim() == 2) {
212 ✗ std::vector<double> auxDisp(3*m_totNumNodeToMatch, 0.);
213 ✗ io[iMesh]->readVariable(idHE, ReadTime, auxDisp.data(), auxDisp.size());
214
215 ✗ for (int i = 0; i < m_totNumNodeToMatch; ++i) {
216 ✗ for (int j = 0; j < 2; ++j) {
217 ✗ m_vectorDisp[2*i+j] = auxDisp[3*i+j];
218 }
219 }
220 ✗ } else {
221 ✗ FEL_ERROR("LinearProblemHarmonicExtension::readDataDisplacement: The value of m_dimension is different from 2 or 3 and this case if not handled.\n");
222 }
223
224 ✗ std::transform(m_vectorDisp.begin(), m_vectorDisp.end(), m_vectorDisp.begin(), std::bind(std::multiplies<double>(),FelisceParam::instance().spaceUnit, std::placeholders::_1));
225
226 ✗ if (m_fstransient->iteration == 1 && FelisceParam::instance().updateMeshByVelocity ) { // for iteration = 1 and for use of velocity,
227 ✗ io[iMesh]->readVariable(0, ReadTime, &m_oldVectorDisp[0], m_oldVectorDisp.size());
228 ✗ std::transform(m_oldVectorDisp.begin(), m_oldVectorDisp.end(), m_oldVectorDisp.begin(), std::bind(std::multiplies<double>(),FelisceParam::instance().spaceUnit, std::placeholders::_1));
229 }
230
231 //If we are stuck at the same iteration (ie RIS model recomputing), we reuse the previous iteration
232 //Miguel said that we maybe "shouldn't re-do the readDataDisplacement anyway"
233 //Should something be done in the case we use the velocity ?
234 ✗ if(FelisceParam::instance().updateMeshByVelocity == false && m_oldVectorDispIteration == m_fstransient->iteration){
235 ✗ std::cout << "Resetting old values of displacement" << std::endl;
236 ✗ for (felInt im = 0; im < m_totNumDofToMatch; im++) {
237 ✗ m_vectorDisp[im] = m_oldVectorDisp[im];
238 ✗ m_oldVectorDisp[im] = m_oldOldVectorDisp[im];
239 }
240 }
241
242
243 ✗ if (FelisceParam::instance().readDisplFromFile || FelisceParam::instance().readPreStressedDisplFromFile || FelisceParam::instance().initALEDispByFile ) {
244 ✗ HarmExtSol().duplicateFrom(this->sequentialSolution());
245 ✗ HarmExtSol().set(0.);
246 felInt iPos;
247 double value;
248 ✗ for (felInt im = 0; im < m_totNumDofToMatch; im++) {
249 ✗ iPos = m_listOfDofToMatch[im];
250 ✗ AOApplicationToPetsc(m_ao,1,&iPos);
251
252 ✗ if(m_fstransient->iteration == 1){
253 ✗ value = m_vectorDisp[im] * m_scaleCoeff; // value is always a displacement
254 }
255 else{
256 ✗ if(FelisceParam::instance().useElasticExtension){ // if useElasticExtension = true, the value used is the incremental value = d^(n+1) - d^(n)
257 ✗ if(FelisceParam::instance().updateMeshByVelocity){
258 //value = (m_vectorDisp[im]*m_fstransient->timeStep) * m_scaleCoeff ; // --> Not the same numerical scheme as M3DISIM!
259 ✗ value = (m_vectorDisp[im] + m_oldVectorDisp[im] ) * ( m_fstransient->timeStep / 2 ) * m_scaleCoeff ; // Based on a Newmark-beta time scheme with gamma = 1/2 and beta = 1/4
260 }
261 else{
262 ✗ value = (m_vectorDisp[im] - m_oldVectorDisp[im]) * m_scaleCoeff ;
263 }
264 }
265 else{ // if useElasticExtension = false, the value used is the absolute value = d^(n+1)
266 ✗ value = m_vectorDisp[im]*m_scaleCoeff;
267 }
268 }
269 ✗ HarmExtSol().setValue(iPos, value, INSERT_VALUES); // matching between _vecdisp and node of label
270
271 ✗ if(FelisceParam::instance().useElasticExtension && !FelisceParam::instance().updateMeshByVelocity){
272 ✗ m_oldOldVectorDisp[im] = m_oldVectorDisp[im];
273 ✗ m_oldVectorDisp[im] = m_vectorDisp[im];
274 }
275
276 ✗ if(FelisceParam::instance().useElasticExtension && FelisceParam::instance().updateMeshByVelocity){
277 ✗ m_oldVectorDisp[im] = m_vectorDisp[im]; // stock the current value of velocity into m_oldVectorDisp which will be used as the previous velocity
278 }
279
280 }
281 ✗ HarmExtSol().assembly();
282 }
283 ✗ m_oldVectorDispIteration = m_fstransient->iteration;
284 }
285
286 570 void LinearProblemHarmonicExtension::finalizeEssBCTransientDerivedProblem()
287 {
288
289
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 570 times.
 570 if (FelisceParam::instance().readDisplFromFile ) {
290
291 ✗ if (m_fstransient->iteration > 0) {
292 ✗ for(std::size_t iBC=0; iBC < m_boundaryConditionList.numDirichletBoundaryCondition(); iBC++) {
293 ✗ BoundaryCondition& BC = * m_boundaryConditionList.Dirichlet(iBC);
294 //===============getListDofOfBC===============//
295 ✗ int iUnknown = BC.getUnknown();
296 ✗ int idVar = m_listUnknown.idVariable(iUnknown);
297 ✗ if(BC.typeValueBC() == EnsightFile) {
298 ✗ felInt idEltInSupportLocal = 0;
299 ✗ felInt idEltInSupportGlobal = 0;
300 felInt idDof;
301 ✗ double aux =0.;
302 ✗ BC.ValueBCInSupportDof().clear();
303 ✗ for (unsigned iSupportDofBC = 0; iSupportDofBC < BC.idEltAndIdSupport().size(); iSupportDofBC++) {
304 ✗ m_supportDofUnknownLocal[iUnknown].getIdElementSupport(BC.idEltAndIdSupport()[iSupportDofBC].first, idEltInSupportLocal);
305 ✗ ISLocalToGlobalMappingApply(m_mappingElem[m_currentMesh],1,&idEltInSupportLocal,&idEltInSupportGlobal);
306 ✗ for(auto it_comp = BC.getComp().begin(); it_comp != BC.getComp().end(); it_comp++) {
307 ✗ dof().loc2glob(idEltInSupportGlobal, BC.idEltAndIdSupport()[iSupportDofBC].second, idVar, *it_comp, idDof);
308 ✗ AOApplicationToPetsc(m_ao,1,&idDof);
309 ✗ HarmExtSol().getValues( 1, &idDof, &aux);
310 ✗ BC.ValueBCInSupportDof().push_back(aux);
311 }
312
313 }
314
315 }
316 }
317 }
318
319 }
320 570 }
321
322 4 void LinearProblemHarmonicExtension::initializeStiffnessVector()
323 {
324
325
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 4 std::cout << "Initializing stiffness std::vector" << std::endl;
326
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 4 std::vector<int> stiffeningLabels = FelisceParam::instance().stiffeningLabels;
327
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 4 std::cout << "Getting num of 3D elts" << std::endl;
328
1/2
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
 4 int nb3DElts = m_mesh[m_currentMesh]->getNumElement3D();
329
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 4 std::cout << "I got the number of 3D elts" << std::endl;
330
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 PetscPrintf(PETSC_COMM_WORLD, "\n\n *** Computing stiffening rule ***\n");
331
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 stiffenElementIds.resize(nb3DElts);
332
1/2
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
 4 std::fill(stiffenElementIds.begin(), stiffenElementIds.end(), 0);
333
334 felInt startIndex;
335 felInt numRefEle;
336
337
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 4 std::cout << "Starting to loop on the label list" << std::endl;
338
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 4 times.
 8 for(unsigned int j = 0; j < stiffeningLabels.size(); ++j){
339 4 startIndex = -1;
340 4 numRefEle = 0;
341 4 int ref = 0;
342
343 4 ref = stiffeningLabels[j];
344
3/6
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
 4 std::cout << "Getting the elements with label " << ref << std::endl;
345
1/2
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
 4 m_mesh[m_currentMesh]->getElementsIDPerRef(m_mesh[m_currentMesh]->Tetra4, ref, &startIndex, &numRefEle);
346
347
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 4 std::cout << "Filling the stiffenIds" << std::endl;
348
2/4
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 if(startIndex != -1 && numRefEle != 0){
349
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 PetscPrintf(PETSC_COMM_WORLD, "Found %d elements (starting at index %d) for reference %d. Activating stiffening for those elements.\n\n", numRefEle, startIndex, ref);
350
2/2
✓ Branch 0 taken 8916 times.
✓ Branch 1 taken 4 times.
 8920 for(int i = startIndex; i < startIndex + numRefEle; ++i){
351 8916 stiffenElementIds[i] = 1;
352 }
353 4 }
354 else{
355 ✗ PetscPrintf(PETSC_COMM_WORLD, "Error for the elements under the unknown label reference %d. Stiffening is std::set to 0 for them.\n\n", ref);
356 }
357 }
358
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 4 std::cout << "Done." << std::endl;
359 4 }
360
361 }
362