GCC Code Coverage Report

Directory: ./
File: Model/bidomainThoraxModel.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 0 79 0.0%
Branches: 0 132 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: E. Schenone
13 //
14
15 // System includes
16
17 // External includes
18
19 // Project includes
20 #include "Model/bidomainThoraxModel.hpp"
21
22 namespace felisce {
23 ✗ BidomainThoraxModel::BidomainThoraxModel():
24 ✗ Model() {
25 ✗ for(std::size_t i=0; i<m_linearProblemBidomainThorax.size(); i++) {
26 ✗ m_linearProblemBidomainThorax[i] = nullptr;
27 }
28 ✗ m_name = "BidomainThorax";
29 }
30
31 ✗ BidomainThoraxModel::~BidomainThoraxModel() {
32 ✗ for(std::size_t i=0; i<m_linearProblemBidomainThorax.size(); i++) {
33 ✗ m_linearProblemBidomainThorax[i] = nullptr;
34 }
35 }
36
37 ✗ void BidomainThoraxModel::initializeDerivedModel() {
38 ✗ if ( FelisceParam::instance().withCVG ){
39 ✗ m_linearProblemBidomainThorax[0]->initPetscVectors();
40 }
41 }
42
43 ✗ void BidomainThoraxModel::initializeDerivedLinearProblem() {
44 ✗ for (std::size_t i=0; i<m_linearProblem.size(); i++) {
45 ✗ m_linearProblemBidomainThorax.push_back(dynamic_cast<LinearProblemBidomainThorax*>(m_linearProblem[i]));
46 }
47 }
48
49 ✗ void BidomainThoraxModel::preAssemblingMatrixRHS(std::size_t iProblem) {
50 ✗ if (m_fstransient->iteration == 0) {
51 ✗ m_linearProblemBidomainThorax[iProblem]->readData(*io());
52 }
53 }
54
55 // Pay attention on the call of this :
56 // call BidomainThoraxModel::writeSolution() function instead of (non-virtual) Model::writeSolution() function
57 ✗ void BidomainThoraxModel::writeSolution()
58 {
59 ✗ Model::writeSolution();
60 }
61
62 ✗ void BidomainThoraxModel::forward() {
63 //Write solution with ensight.
64 ✗ BidomainThoraxModel::writeSolution();
65 //Advance time step.
66 ✗ updateTime(FlagMatrixRHS::only_rhs);
67 ✗ printNewTimeIterationBanner();
68
69 ✗ preAssemblingMatrixRHS(0);
70
71 // Assembling of matrices at first iteration (because of needs some coefficients not defined at iteration 0).
72 ✗ if ( m_fstransient->iteration == 1 ) {
73 ✗ m_linearProblemBidomainThorax[0]->assembleMatrixRHS(MpiInfo::rankProc(), FlagMatrixRHS::matrix_and_rhs);
74 //Apply boundary conditions.
75 ✗ m_linearProblem[0]->finalizeEssBCConstantInT();
76 ✗ m_linearProblem[0]->finalizeEssBCTransient();
77 ✗ m_linearProblem[0]->applyBC(FelisceParam::instance().essentialBoundaryConditionsMethod, MpiInfo::rankProc(), FlagMatrixRHS::matrix_and_rhs, FlagMatrixRHS::matrix_and_rhs);
78 } else {
79 //Apply boundary conditions.
80 ✗ m_linearProblem[0]->finalizeEssBCConstantInT();
81 ✗ m_linearProblem[0]->finalizeEssBCTransient();
82 ✗ m_linearProblem[0]->applyBC(FelisceParam::instance().essentialBoundaryConditionsMethod, MpiInfo::rankProc(), FlagMatrixRHS::only_rhs, FlagMatrixRHS::only_rhs);
83 }
84
85 ✗ if (FelisceParam::instance().writeMatrixECG) {
86 ✗ m_linearProblemBidomainThorax[0]->calculateRHS();
87 }
88
89 // Apply CVGraph interface conditions
90 ✗ postAssemblingMatrixRHS(0);
91
92 //Solve the system _Matrix[0]*m_sol=_RHS.
93 ✗ m_linearProblemBidomainThorax[0]->solve(MpiInfo::rankProc(), MpiInfo::numProc());
94
95 // User-defined function to calculate residual = lines of transfer matrix
96 ✗ if (FelisceParam::instance().writeMatrixECG ) {
97 ✗ m_linearProblemBidomainThorax[0]->calculateRes();
98 }
99
100 }
101
102 #ifdef FELISCE_WITH_CVGRAPH
103 ✗ void BidomainThoraxModel::startIterationCVG() {
104 ✗ cvgMainSlave* slave=m_linearProblem[0]->slave();
105 ✗ if (slave==nullptr) {
106 ✗ FEL_ERROR("slave not initialized");
107 }
108 ✗ if( slave->newTimeStep() ) {
109 ✗ std::stringstream msg;
110 ✗ msg<<"Starting new time step " << m_fstransient->iteration <<" , preparing forward"<<std::endl;
111 ✗ PetscPrintf(MpiInfo::petscComm(),"%s",msg.str().c_str());
112
113 //Write solution with ensight.
114 ✗ BidomainThoraxModel::writeSolution();
115 //Advance time step.
116
117 ✗ if ( m_fstransient->iteration == 0 || m_fstransient->iteration == 1 ) {
118 //std::cout << "1.1 Slave new time step = " << m_fstransient->iteration << std::endl;
119 ✗ updateTime(FlagMatrixRHS::matrix_and_rhs);
120 ✗ printNewTimeIterationBanner();
121 //m_linearProblem[0]->clearMatrixRHS(FlagMatrixRHS::matrix_and_rhs);// Done in updateTime()
122 ✗ m_linearProblem[0]->assembleMatrixRHS(MpiInfo::rankProc(),FlagMatrixRHS::matrix_and_rhs);
123 ✗ m_linearProblem[0]->createAndCopyMatrixRHSWithoutBC();
124
125 ✗ m_linearProblem[0]->finalizeEssBCConstantInT();
126 //std::cout << "FinalizeEssBCTransient()...";
127 ✗ m_linearProblem[0]->finalizeEssBCTransient();
128 //std::cout << " Done FEBT." << std::endl;
129 //std::cout << "Apply BC...";
130 ✗ if (slave->thereIsAtLeastOneRobinCondition()){
131 //std::cout << " CVG Robin: we apply the robin condition in the MATRIX ONLY. The RHS parts are handled in cvgAddExternalResidualToRHS() and cvgAddRobinToRHS()." << std::endl;
132 ✗ m_linearProblem[0]->applyBC(FelisceParam::instance().essentialBoundaryConditionsMethod, MpiInfo::rankProc(), FlagMatrixRHS::matrix_and_rhs, FlagMatrixRHS::only_matrix,0, true, ApplyNaturalBoundaryConditions::yes);
133 ✗ if (m_fstransient->iteration == 1){
134 //std::cout << "Building the mass matrix at the CVG boundary" << std::endl;
135 ✗ m_linearProblemBidomainThorax[0]->initMassBoundaryForCVG();
136 }
137 }
138 else {
139 ✗ m_linearProblem[0]->applyBC(FelisceParam::instance().essentialBoundaryConditionsMethod, MpiInfo::rankProc(), FlagMatrixRHS::matrix_and_rhs, FlagMatrixRHS::matrix_and_rhs,0, true, ApplyNaturalBoundaryConditions::no);
140 }
141 }
142 else{
143 //std::cout << "1.2 Slave new time step = " << m_fstransient->iteration << std::endl;
144 ✗ updateTime(FlagMatrixRHS::only_rhs);
145 ✗ printNewTimeIterationBanner();
146
147 ✗ m_linearProblem[0]->finalizeEssBCConstantInT();
148 ✗ m_linearProblem[0]->finalizeEssBCTransient();
149 ✗ m_linearProblem[0]->applyBC(FelisceParam::instance().essentialBoundaryConditionsMethod, MpiInfo::rankProc(), FlagMatrixRHS::only_rhs, FlagMatrixRHS::only_rhs, 0, true, ApplyNaturalBoundaryConditions::no);
150
151 //m_linearProblemBidomainThorax[0]->addCurrentToRHS();
152 }
153
154 ✗ } else if( slave->redoTimeStep() ){
155 ✗ std::stringstream msg;
156 ✗ msg << "......................................" << std::endl;
157 ✗ msg << "Redo time step at time t= " << m_fstransient->time << std::endl;
158 ✗ PetscPrintf(MpiInfo::petscComm(),"%s",msg.str().c_str());
159
160
161 ✗ m_linearProblem[0]->clearMatrixRHS(FlagMatrixRHS::only_rhs);
162
163 // No call to assembleMatrixRHS. THIS IS DANGEROUS BUT PROVIDES A HUGE SPEEDUP
164 // This only works because the RHS without BC is 0. It is not the case if there are source terms for instance... Eliott
165 // I guess many things (e.g. loops on elements) are done in assembleMatrixRHS() even when FlagMatrixRHS is std::set to only_rhs and are therefore useless.
166 // m_linearProblem[0]->assembleMatrixRHS(MpiInfo::rankProc(),FlagMatrixRHS::only_rhs);
167
168 ✗ m_linearProblem[0]->finalizeEssBCConstantInT();
169 ✗ m_linearProblem[0]->finalizeEssBCTransient();
170 ✗ m_linearProblem[0]->applyBC(FelisceParam::instance().essentialBoundaryConditionsMethod, MpiInfo::rankProc(), FlagMatrixRHS::only_rhs, FlagMatrixRHS::only_rhs, 0, true, ApplyNaturalBoundaryConditions::no);
171 //m_linearProblemBidomainThorax[0]->addCurrentToRHS();
172
173 ✗ } else {
174 ✗ slave->msgInfo().print(1);
175 ✗ FEL_ERROR("Error: strange timeStatus");
176 }
177
178 ✗ if (slave->thereIsAtLeastOneRobinCondition()){
179 //std::cout << " CVG Robin: Add remaining term in RHS: beta_torso*massBD*u_torso" << std::endl;
180 ✗ m_linearProblemBidomainThorax[0]->cvgAddRobinToRHS();
181 }
182
183 //Solve the system _Matrix[0]*m_sol=_RHS.
184 //std::cout <<std::endl<<std::endl<< "####### SOLVE THE LINEAR SYSTEM ######" <<std::endl<< std::endl;
185 ✗ m_linearProblemBidomainThorax[0]->solve(MpiInfo::rankProc(), MpiInfo::numProc());
186
187
188 }
189 #endif
190
191
192 }
193