GCC Code Coverage Report

Directory: ./
File: Model/model.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 277 427 64.9%
Branches: 274 720 38.1%
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 & M. Fragu
13 //
14
15 // System includes
16 #include <ctime>
17
18 // External includes
19
20 // Project includes
21 #include "Core/mpiInfo.hpp"
22 #include "Model/model.hpp"
23 #include "Solver/CVGraphInterface.hpp"
24 #include "Solver/linearProblemHarmonicExtension.hpp"
25 #include "Tools/fe_utilities.hpp"
26
27 namespace felisce
28 {
29
30
4/8
✓ Branch 2 taken 467 times.
✗ Branch 3 not taken.
✓ Branch 13 taken 467 times.
✗ Branch 14 not taken.
✓ Branch 17 taken 467 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 467 times.
✗ Branch 21 not taken.
 467 Model::Model()
31 {
32 467 }
33
34 /***********************************************************************************/
35 /***********************************************************************************/
36
37 934 Model::~Model() {
38
39
2/2
✓ Branch 0 taken 461 times.
✓ Branch 1 taken 6 times.
 934 if (m_linearProblemIsInitialized) {
40
2/2
✓ Branch 2 taken 459 times.
✓ Branch 3 taken 2 times.
 922 if(FelisceParam::instance(this->instanceIndex()).duplicateSupportDof == false)
41 918 writeSolution();
42 }
43
44 // Prints chrono table
45 934 auto& r_instance = FelisceParam::instance(this->instanceIndex());
46
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 467 times.
 934 if (r_instance.chronoToScreen) {
47 ✗ r_instance.timer.ToScreen();
48 }
49
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 467 times.
 934 if (r_instance.chronoToFile) {
50 ✗ r_instance.timer.resultDir() = r_instance.resultDir;
51 ✗ r_instance.timer.ToFile();
52 }
53
54
2/2
✓ Branch 1 taken 513 times.
✓ Branch 2 taken 467 times.
 1960 for(std::size_t i=0; i<m_linearProblem.size(); i++) {
55
2/4
✓ Branch 1 taken 513 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 513 times.
✗ Branch 5 not taken.
 1026 if(m_linearProblem[i]) delete m_linearProblem[i];
56 }
57
58 #ifdef FELISCE_WITH_CVGRAPH
59 // necessary to destroy manually all petsc objects before the petsc finalize
60
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 467 times.
 934 if ( ! m_rhsWithForceTermWithoutBC.isNull() ) {
61 ✗ m_rhsWithForceTermWithoutBC.destroy();
62 }
63 #endif
64
2/2
✓ Branch 0 taken 461 times.
✓ Branch 1 taken 6 times.
 934 if (m_linearProblemIsInitialized) {
65 922 m_U_0.destroy();
66
67 922 PetscPrintf(MpiInfo::petscComm(), "\n==============================================================\n");
68 922 PetscPrintf(MpiInfo::petscComm(), "FELiScE %s : normal end\n",m_name.c_str());
69 // current date and time system
70 922 time_t now = time(nullptr);
71 922 char* date_time = ctime(&now);
72 922 PetscPrintf(MpiInfo::petscComm(), "%s",date_time);
73 922 PetscPrintf(MpiInfo::petscComm(), "==============================================================\n");
74
75
2/2
✓ Branch 0 taken 453 times.
✓ Branch 1 taken 8 times.
 922 if (m_callPetscInitialize)
76 906 PetscFinalize();
77 }
78
79
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 466 times.
 934 if (m_eigenProblemIsInitialized) {
80
81 2 PetscPrintf(MpiInfo::petscComm(), "\n==============================================================\n");
82 2 PetscPrintf(MpiInfo::petscComm(), "FELiScE %s : normal end\n",m_name.c_str());
83 // current date and time system
84 2 time_t now = time(nullptr);
85 2 char* date_time = ctime(&now);
86 2 PetscPrintf(MpiInfo::petscComm(), "%s",date_time);
87 2 PetscPrintf(MpiInfo::petscComm(), "==============================================================\n");
88
89 #ifdef FELISCE_WITH_SLEPC
90
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 2 if (m_callPetscInitialize)
91 2 SlepcFinalize();
92 #endif
93 }
94
95 934 m_ios.clear();
96 }
97
98 /***********************************************************************************/
99 /***********************************************************************************/
100
101 namespace { // anonymous
102 // Print which architecture is used
103 250 void printArchitecture(const MPI_Comm& comm) {
104
1/2
✓ Branch 1 taken 250 times.
✗ Branch 2 not taken.
 250 std::ostringstream oconv;
105
1/2
✓ Branch 1 taken 250 times.
✗ Branch 2 not taken.
 250 oconv << "FELiScE ";
106
1/2
✓ Branch 1 taken 250 times.
✗ Branch 2 not taken.
 250 oconv << 8 * sizeof(void*);
107
1/2
✓ Branch 1 taken 250 times.
✗ Branch 2 not taken.
 250 oconv << " bits ";
108
2/4
✓ Branch 1 taken 250 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 250 times.
✗ Branch 6 not taken.
 250 PetscPrintf(comm,"%s",oconv.str().c_str());
109 250 }
110 } // namespace anonymous
111
112 /***********************************************************************************/
113 /***********************************************************************************/
114
115 462 void Model::initializeModel(CommandLineOption& opt, FelisceTransient& fstransient, const MPI_Comm comm, const bool callPetscInitialize)
116 {
117 // Update instance index
118 462 m_instanceIndex = opt.instanceIndex();
119
1/2
✓ Branch 2 taken 462 times.
✗ Branch 3 not taken.
 462 auto& r_instance = FelisceParam::instance(this->instanceIndex());
120
121 // Get timer
122 462 auto& r_timer = r_instance.timer;
123
2/4
✓ Branch 2 taken 462 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 462 times.
✗ Branch 6 not taken.
 462 r_timer.Start("Model::initializeModel");
124
125
1/2
✓ Branch 1 taken 462 times.
✗ Branch 2 not taken.
 462 m_fstransient = felisce::make_shared<FelisceTransient>(fstransient);
126 462 m_callPetscInitialize = callPetscInitialize;
127
1/2
✓ Branch 2 taken 462 times.
✗ Branch 3 not taken.
 462 const std::string help = "Felisce\n\n";
128
129 //Initialize Petsc (MPI)
130 //======================
131
1/2
✓ Branch 0 taken 462 times.
✗ Branch 1 not taken.
 462 if (comm == MPI_COMM_WORLD) { // default communicator
132 #ifdef FELISCE_WITH_SLEPC
133
1/2
✓ Branch 4 taken 462 times.
✗ Branch 5 not taken.
 462 SlepcInitialize(&opt.argcPetsc(), &opt.argvPetsc(), FELISCE_PETSC_NULLPTR, help.c_str());
134 #else
135 PetscInitialize(&opt.argcPetsc(), &opt.argvPetsc(), FELISCE_PETSC_NULLPTR, help.c_str());
136 #endif
137
1/2
✓ Branch 1 taken 462 times.
✗ Branch 2 not taken.
 462 MpiInfo::petscComm()=MPI_COMM_WORLD;
138 } else { // user defined communicator (for use with CVGraph for instance...) This is a delicate part, think about before touching it!Saverio
139 // From petsc documentation:
140 // If you wish PETSc code to run ONLY on a subcommunicator of MPI_COMM_WORLD, create that communicator first and assign it to PETSC_COMM_WORLD BEFORE calling PetscInitialize().
141 // Thus if you are running a four process job and two processes will run PETSc and have PetscInitialize() and PetscFinalize() and two process will not, then do this.
142 // If ALL processes in the job are using PetscInitialize() and PetscFinalize() then you don't need to do this, even if different subcommunicators of the job are doing different things with PETSc.
143 ✗ PETSC_COMM_WORLD = comm;
144 ✗ if (callPetscInitialize) {
145 #ifdef FELISCE_WITH_SLEPC
146 ✗ SlepcInitialize(&opt.argcPetsc(),&opt.argvPetsc(), FELISCE_PETSC_NULLPTR, help.c_str());
147 #else
148 PetscInitialize(&opt.argcPetsc(),&opt.argvPetsc(), FELISCE_PETSC_NULLPTR, help.c_str());
149 #endif
150 }
151 ✗ MpiInfo::petscComm() = PETSC_COMM_WORLD;
152 }
153
154 int rankProc,numProc;
155
2/4
✓ Branch 1 taken 462 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 462 times.
✗ Branch 5 not taken.
 462 MPI_Comm_rank(MpiInfo::petscComm(),&rankProc);
156
2/4
✓ Branch 1 taken 462 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 462 times.
✗ Branch 5 not taken.
 462 MPI_Comm_size(MpiInfo::petscComm(),&numProc);
157
1/2
✓ Branch 1 taken 462 times.
✗ Branch 2 not taken.
 462 MpiInfo::rankProc()=rankProc;
158
1/2
✓ Branch 1 taken 462 times.
✗ Branch 2 not taken.
 462 MpiInfo::numProc()=numProc;
159
160 // Printing logo
161
3/4
✓ Branch 1 taken 462 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 128 times.
✓ Branch 4 taken 334 times.
 462 if (MpiInfo::rankProc() == 0) {
162 std::cout << " ______ ______ _ _ _____ ______ \n"
163 << " | ____| ____| | (_)/ ____| | ____| \n"
164 << " | |__ | |__ | | _| (___ ___| |__ \n"
165 << " | __| | __| | | | |\\___ \\ / __| __| \n"
166 << " | | | |____| |____| |____) | (__| |____ \n"
167 << " |_| |______|______|_|_____/ \\___|______| \n"
168
8/16
✓ Branch 1 taken 128 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 128 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 128 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 128 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 128 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 128 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 128 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 128 times.
✗ Branch 23 not taken.
 128 << " Finite Elements for Life Sciences and Engineering\n" << std::endl;
169 }
170
171 // current date and time system
172 462 time_t now = std::time(nullptr);
173 462 char* date_time = ctime(&now);
174
3/4
✓ Branch 1 taken 462 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 250 times.
✓ Branch 4 taken 212 times.
 462 if (FelisceParam::verbose()) {
175
2/4
✓ Branch 1 taken 250 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 250 times.
✗ Branch 5 not taken.
 250 PetscPrintf(MpiInfo::petscComm(), "\n================================================================\n");
176
177
1/2
✓ Branch 1 taken 250 times.
✗ Branch 2 not taken.
 250 printArchitecture(comm);
178
179
5/8
✓ Branch 1 taken 250 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 10 times.
✓ Branch 4 taken 240 times.
✓ Branch 6 taken 10 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 10 times.
✗ Branch 10 not taken.
 250 if(MpiInfo::numProc()==1) PetscPrintf(MpiInfo::petscComm(), "on 1 processor \n");
180
3/6
✓ Branch 1 taken 240 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 240 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 240 times.
✗ Branch 8 not taken.
 240 else PetscPrintf(MpiInfo::petscComm(), "on %d processors \n", MpiInfo::numProc());
181
2/4
✓ Branch 1 taken 250 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 250 times.
✗ Branch 5 not taken.
 250 PetscPrintf(MpiInfo::petscComm(), "%s", date_time);
182
2/4
✓ Branch 1 taken 250 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 250 times.
✗ Branch 5 not taken.
 250 PetscPrintf(MpiInfo::petscComm(), "=================================================================\n");
183 }
184
2/4
✓ Branch 1 taken 462 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 462 times.
✗ Branch 5 not taken.
 462 r_instance.print(FelisceParam::verbose());
185
186 //Read the mesh file
187 //==================
188
1/2
✓ Branch 1 taken 462 times.
✗ Branch 2 not taken.
 462 readMesh();
189
190 462 m_modelIsInitialized = true;
191
192
2/2
✓ Branch 1 taken 468 times.
✓ Branch 2 taken 462 times.
 930 for(std::size_t iio=0; iio<m_ios.size(); ++iio) {
193
1/2
✓ Branch 3 taken 468 times.
✗ Branch 4 not taken.
 468 m_ios[iio]->readInputFile();
194 }
195
196 //Error message if method of characteristics is used in parallel, to be removed when this method is implemented in parallel
197
3/4
✓ Branch 0 taken 446 times.
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 446 times.
 462 if (r_instance.NSequationFlag == 1 && r_instance.characteristicMethod != 0) {
198 ✗ if (MpiInfo::numProc() > 1)
199 ✗ FEL_ERROR("Method of characteristics has not been implemented yet in parallel!");
200 }
201
202 //Verify if this model use CVGraphInterface
203 //==================================
204
1/2
✓ Branch 1 taken 462 times.
✗ Branch 2 not taken.
 462 if ( !strcmp(r_instance.CVGraphInterface.c_str(),"default")) {
205 462 m_initializeCVGraphInterface = false;
206 } else {
207 ✗ m_initializeCVGraphInterface = true;
208 }
209
210
2/4
✓ Branch 2 taken 462 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 462 times.
✗ Branch 6 not taken.
 462 r_timer.Stop("Model::initializeModel");
211 462 }
212
213 /***********************************************************************************/
214 /***********************************************************************************/
215
216 456 void Model::initializeLinearProblem(const std::vector<LinearProblem*>& linearPb, bool doUseSNES)
217 {
218
1/2
✓ Branch 3 taken 456 times.
✗ Branch 4 not taken.
 456 std::vector<bool> aDoUseSNES(linearPb.size(), doUseSNES);
219
1/2
✓ Branch 1 taken 456 times.
✗ Branch 2 not taken.
 456 initializeLinearProblem(linearPb,aDoUseSNES);
220 456 }
221
222 /***********************************************************************************/
223 /***********************************************************************************/
224
225 461 void Model::initializeLinearProblem(const std::vector<LinearProblem*>& linearPb, const std::vector<bool>& doUseSNES)
226 {
227 461 auto& r_instance = FelisceParam::instance(this->instanceIndex());
228
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 461 times.
 461 FEL_CHECK( doUseSNES.size() == linearPb.size(),"Number of linear problems and doUseSNES differ");
229
230 // Check connections between variables and meshes
231 461 checkConnectionVariablesMeshs();
232
233 // Get timer
234 461 auto& r_timer = r_instance.timer;
235
236 // Initilize solution backup
237
2/4
✓ Branch 2 taken 461 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 461 times.
✗ Branch 6 not taken.
 461 r_timer.Start("Model::Initialize solution backup");
238 461 m_solutionBackup.resize(linearPb.size());
239
2/2
✓ Branch 1 taken 513 times.
✓ Branch 2 taken 461 times.
 974 for(std::size_t i=0; i<m_solutionBackup.size(); ++i) {
240 513 m_solutionBackup[i].initialize(r_instance.inputDirectory, r_instance.inputFile, r_instance.inputMesh,
241 513 r_instance.outputMesh, r_instance.meshDir, r_instance.resultDir, r_instance.prefixName);
242 }
243
2/4
✓ Branch 2 taken 461 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 461 times.
✗ Branch 6 not taken.
 461 r_timer.Stop("Model::Initialize solution backup");
244
245 // Copy all linear system in m_linearProblem
246
2/2
✓ Branch 1 taken 513 times.
✓ Branch 2 taken 461 times.
 974 for (std::size_t ipb = 0; ipb < linearPb.size(); ipb++) {
247 513 m_linearProblem.push_back(linearPb[ipb]);
248 513 linearPb[ipb]->identifierProblem() = ipb;
249 }
250
251 // Definition of specific linear problems in each model
252
2/4
✓ Branch 2 taken 461 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 461 times.
✗ Branch 6 not taken.
 461 r_timer.Start("Model::initializeDerivedLinearProblem");
253 461 initializeDerivedLinearProblem();
254
2/4
✓ Branch 2 taken 461 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 461 times.
✗ Branch 6 not taken.
 461 r_timer.Stop("Model::initializeDerivedLinearProblem");
255
256 // For each linear system,
257
2/2
✓ Branch 1 taken 513 times.
✓ Branch 2 taken 461 times.
 974 for (std::size_t ipb = 0; ipb < linearPb.size(); ipb++) {
258 // Set the model
259 513 m_linearProblem[ipb]->model() = this;
260
261 // Initialize it
262
3/6
✓ Branch 2 taken 513 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 513 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 513 times.
✗ Branch 9 not taken.
 513 r_timer.Start("LinearProblem" + std::to_string(ipb) + "::initialize");
263
1/2
✓ Branch 5 taken 513 times.
✗ Branch 6 not taken.
 513 m_linearProblem[ipb]->initialize(m_mesh, m_fstransient, MpiInfo::petscComm(), doUseSNES[ipb]);
264
3/6
✓ Branch 2 taken 513 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 513 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 513 times.
✗ Branch 9 not taken.
 513 r_timer.Stop("LinearProblem" + std::to_string(ipb) + "::initialize");
265
266
2/2
✓ Branch 2 taken 76 times.
✓ Branch 3 taken 437 times.
 513 if (r_instance.solver.size() < linearPb.size()) {
267 76 m_linearProblem[ipb]->fixIdOfTheProblemSolver(0);
268 } else {
269 437 m_linearProblem[ipb]->fixIdOfTheProblemSolver(ipb);
270 }
271
272 // And compute the dof.
273
3/6
✓ Branch 2 taken 513 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 513 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 513 times.
✗ Branch 9 not taken.
 513 r_timer.Start("LinearProblem" + std::to_string(ipb) + "::computeDof");
274 513 m_linearProblem[ipb]->computeDof(MpiInfo::numProc(), MpiInfo::rankProc());
275
3/6
✓ Branch 2 taken 513 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 513 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 513 times.
✗ Branch 9 not taken.
 513 r_timer.Stop("LinearProblem" + std::to_string(ipb) + "::computeDof");
276 }
277
278 // Define initial conditions
279
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 461 times.
 461 if (r_instance.hasInitialCondition) {
280 ✗ for (std::size_t ipb = 0; ipb < linearPb.size(); ipb++) {
281 ✗ for (std::size_t ii = 0; ii < r_instance.valueInitCond.size(); ii++) {
282 ✗ m_initialCondition.addVariable(*m_linearProblem[ipb]->listVariable().getVariable(r_instance.nameVariableInitCond[ii]));
283 }
284 ✗ m_initialCondition.print(FelisceParam::verbose(),std::cout);
285 }
286 }
287
288 // For each linear system
289
2/2
✓ Branch 1 taken 513 times.
✓ Branch 2 taken 461 times.
 974 for (std::size_t ipb = 0; ipb < linearPb.size(); ipb++) {
290 // userChangePattern: user function to change the pattern if needed in the specific application.
291
3/6
✓ Branch 2 taken 513 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 513 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 513 times.
✗ Branch 9 not taken.
 513 r_timer.Start("LinearProblem" + std::to_string(ipb) + "::userChangePattern");
292 513 m_linearProblem[ipb]->userChangePattern(MpiInfo::numProc(), MpiInfo::rankProc());
293
3/6
✓ Branch 2 taken 513 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 513 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 513 times.
✗ Branch 9 not taken.
 513 r_timer.Stop("LinearProblem" + std::to_string(ipb) + "::userChangePattern");
294
295 // cutMesh: partitioning of the dof and elements.
296
3/6
✓ Branch 2 taken 513 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 513 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 513 times.
✗ Branch 9 not taken.
 513 r_timer.Start("LinearProblem" + std::to_string(ipb) + "::cutMesh");
297 513 m_linearProblem[ipb]->cutMesh();
298
3/6
✓ Branch 2 taken 513 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 513 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 513 times.
✗ Branch 9 not taken.
 513 r_timer.Stop("LinearProblem" + std::to_string(ipb) + "::cutMesh");
299
300 // allocateMatrix: Allocate memory for the matrix _Matrix and std::vector _RHS in the linearProblem class
301
3/6
✓ Branch 2 taken 513 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 513 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 513 times.
✗ Branch 9 not taken.
 513 r_timer.Start("LinearProblem" + std::to_string(ipb) + "::allocateMatrix");
302 513 m_linearProblem[ipb]->allocateMatrix();
303
3/6
✓ Branch 2 taken 513 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 513 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 513 times.
✗ Branch 9 not taken.
 513 r_timer.Stop("LinearProblem" + std::to_string(ipb) + "::allocateMatrix");
304 }
305
306 // Specific initialization of the user model (useful after allocate Matrix of linear problem)
307
2/4
✓ Branch 2 taken 461 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 461 times.
✗ Branch 6 not taken.
 461 r_timer.Start("Model::initializeDerivedModel");
308 461 initializeDerivedModel();
309
2/4
✓ Branch 2 taken 461 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 461 times.
✗ Branch 6 not taken.
 461 r_timer.Stop("Model::initializeDerivedModel");
310
311
2/4
✓ Branch 2 taken 461 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 461 times.
✗ Branch 6 not taken.
 461 r_timer.Start("Model::setExternalVec");
312 // Make the connection between the different linear pb (to be defined in the derived class)
313 461 setExternalVec();
314
2/4
✓ Branch 2 taken 461 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 461 times.
✗ Branch 6 not taken.
 461 r_timer.Stop("Model::setExternalVec");
315
316
2/2
✓ Branch 1 taken 513 times.
✓ Branch 2 taken 461 times.
 974 for (std::size_t ipb = 0; ipb < linearPb.size(); ipb++) {
317 //Define boundary condition of the problem
318
3/6
✓ Branch 2 taken 513 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 513 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 513 times.
✗ Branch 9 not taken.
 513 r_timer.Start("LinearProblem" + std::to_string(ipb) + "::defineBC");
319 513 m_linearProblem[ipb]->defineBC();
320
3/6
✓ Branch 2 taken 513 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 513 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 513 times.
✗ Branch 9 not taken.
 513 r_timer.Stop("LinearProblem" + std::to_string(ipb) + "::defineBC");
321
322 //Compute boundary condition
323
3/6
✓ Branch 2 taken 513 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 513 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 513 times.
✗ Branch 9 not taken.
 513 r_timer.Start("LinearProblem" + std::to_string(ipb) + "::Compute boundary condition");
324 513 m_linearProblem[ipb]->determineDofAssociateToLabel();
325 513 m_linearProblem[ipb]->finalizeEssBCConstantInT();
326 513 m_linearProblem[ipb]->finalizeEssBCTransient();
327
3/6
✓ Branch 2 taken 513 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 513 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 513 times.
✗ Branch 9 not taken.
 513 r_timer.Stop("LinearProblem" + std::to_string(ipb) + "::Compute boundary condition");
328
329 //Apply specific operations before assembling
330
3/6
✓ Branch 2 taken 513 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 513 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 513 times.
✗ Branch 9 not taken.
 513 r_timer.Start("LinearProblem" + std::to_string(ipb) + "::preAssemblingMatrixRHS");
331 513 preAssemblingMatrixRHS(ipb);
332
3/6
✓ Branch 2 taken 513 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 513 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 513 times.
✗ Branch 9 not taken.
 513 r_timer.Stop("LinearProblem" + std::to_string(ipb) + "::preAssemblingMatrixRHS");
333
334 // Gather to write initial solution.
335 513 m_linearProblem[ipb]->gatherSolution();
336
337
3/6
✓ Branch 2 taken 513 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 513 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 513 times.
✗ Branch 9 not taken.
 513 r_timer.Start("LinearProblem" + std::to_string(ipb) + "::buildSolver");
338 513 m_linearProblem[ipb]->buildSolver();
339
3/6
✓ Branch 2 taken 513 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 513 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 513 times.
✗ Branch 9 not taken.
 513 r_timer.Stop("LinearProblem" + std::to_string(ipb) + "::buildSolver");
340 }
341
342 // Initialize object which manage output.
343
2/2
✓ Branch 1 taken 467 times.
✓ Branch 2 taken 461 times.
 928 for(std::size_t iio=0; iio<m_ios.size(); ++iio) {
344 467 m_ios[iio]->initializeOutput();
345 }
346
347
2/2
✓ Branch 1 taken 513 times.
✓ Branch 2 taken 461 times.
 974 for(std::size_t ipb = 0; ipb < m_linearProblem.size(); ipb++) {
348
3/6
✓ Branch 2 taken 513 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 513 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 513 times.
✗ Branch 9 not taken.
 513 r_timer.Start("LinearProblem" + std::to_string(ipb) + "::clearMatrixRHS");
349 513 m_linearProblem[ipb]->clearMatrixRHS();
350
3/6
✓ Branch 2 taken 513 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 513 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 513 times.
✗ Branch 9 not taken.
 513 r_timer.Stop("LinearProblem" + std::to_string(ipb) + "::clearMatrixRHS");
351
3/6
✓ Branch 2 taken 513 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 513 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 513 times.
✗ Branch 9 not taken.
 513 r_timer.Start("LinearProblem" + std::to_string(ipb) + "::InitializeDerivedAttributes");
352 513 m_linearProblem[ipb]->InitializeDerivedAttributes();
353
3/6
✓ Branch 2 taken 513 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 513 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 513 times.
✗ Branch 9 not taken.
 513 r_timer.Stop("LinearProblem" + std::to_string(ipb) + "::InitializeDerivedAttributes");
354 }
355
2/4
✓ Branch 2 taken 461 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 461 times.
✗ Branch 6 not taken.
 461 r_timer.Start("Model::setInitialCondition");
356 461 setInitialCondition();
357
2/4
✓ Branch 2 taken 461 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 461 times.
✗ Branch 6 not taken.
 461 r_timer.Stop("Model::setInitialCondition");
358
359 461 m_fstransient->iteration=0;
360
361 461 m_linearProblemIsInitialized = true;
362 461 }
363
364
365 /***********************************************************************************/
366 /***********************************************************************************/
367
368 461 void Model::checkConnectionVariablesMeshs()
369 {
370 461 auto& r_instance = FelisceParam::instance(this->instanceIndex());
371
372
2/2
✓ Branch 1 taken 455 times.
✓ Branch 2 taken 6 times.
 461 if( m_mesh.size() == 1 ) {
373 // If only one mesh provided -> every variables is defined on that mesh
374
1/2
✓ Branch 2 taken 455 times.
✗ Branch 3 not taken.
 455 r_instance.idMesh.resize(r_instance.physicalVariable.size(), 0);
375 }
376 else {
377
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 6 times.
 6 FEL_CHECK(r_instance.idMesh.size() == r_instance.physicalVariable.size(), "Error idmesh->size() != physicalVariable.size()");
378
2/2
✓ Branch 1 taken 19 times.
✓ Branch 2 taken 6 times.
 25 for(std::size_t i=0; i<r_instance.idMesh.size(); ++i){
379
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 19 times.
 19 FEL_CHECK(static_cast<std::size_t>(r_instance.idMesh[i]) < m_mesh.size(),"Error in idMesh mesh ID greater than m_mesh.size()");
380 }
381 }
382 461 }
383
384 /***********************************************************************************/
385 /***********************************************************************************/
386
387 393 void Model::initializeLinearProblem(LinearProblem* linearPb, bool doUseSNES)
388 {
389
1/2
✓ Branch 2 taken 393 times.
✗ Branch 3 not taken.
 393 std::vector<LinearProblem*> linPb(1, linearPb);
390
1/2
✓ Branch 1 taken 393 times.
✗ Branch 2 not taken.
 393 initializeLinearProblem(linPb, doUseSNES);
391 393 }
392
393 /***********************************************************************************/
394 /***********************************************************************************/
395
396 53 void Model::finalizeLinearProblem(const std::vector<LinearProblem*>& linearPb)
397 {
398 // For each linear system,
399
2/2
✓ Branch 1 taken 53 times.
✓ Branch 2 taken 53 times.
 106 for (std::size_t ipb = 0; ipb < linearPb.size(); ipb++) {
400
1/2
✓ Branch 4 taken 53 times.
✗ Branch 5 not taken.
 53 m_linearProblem[ipb]->finalize(m_mesh, m_fstransient, MpiInfo::petscComm());
401 }
402 53 }
403
404 /***********************************************************************************/
405 /***********************************************************************************/
406
407 53 void Model::finalizeLinearProblem(LinearProblem* linearPb)
408 {
409
1/2
✓ Branch 2 taken 53 times.
✗ Branch 3 not taken.
 53 std::vector<LinearProblem*> linPb(1, linearPb);
410
1/2
✓ Branch 1 taken 53 times.
✗ Branch 2 not taken.
 53 finalizeLinearProblem(linPb);
411 53 }
412
413 /***********************************************************************************/
414 /***********************************************************************************/
415
416 11420 void Model::updateTime(const FlagMatrixRHS flagMatrixRHS)
417 {
418
2/2
✓ Branch 1 taken 12232 times.
✓ Branch 2 taken 11420 times.
 23652 for(std::size_t i=0; i<m_linearProblem.size(); i++) {
419 12232 m_linearProblem[i]->clearMatrixRHS(flagMatrixRHS);
420 }
421 11420 m_fstransient->iteration++;
422 11420 m_fstransient->time +=m_fstransient->timeStep;
423 11420 }
424
425 /***********************************************************************************/
426 /***********************************************************************************/
427
428 24 void Model::SolveStaticProblem()
429 {
430 24 PetscPrintf(MpiInfo::petscComm(), "----------------------------------------------\n");
431 24 PetscPrintf(MpiInfo::petscComm(), "--------------- Static problem ---------------\n");
432 24 PetscPrintf(MpiInfo::petscComm(), "----------------------------------------------\n");
433
434 // Retrieve linear problem
435 24 auto p_linear_problem = m_linearProblem[0];
436
437 // Solve the system.
438 24 p_linear_problem->iterativeSolve(MpiInfo::rankProc(), MpiInfo::numProc());
439 24 }
440
441 /***********************************************************************************/
442 /***********************************************************************************/
443
444 344 void Model::SolveDynamicProblem()
445 {
446 344 PetscPrintf(MpiInfo::petscComm(), "----------------------------------------------\n");
447 344 PetscPrintf(MpiInfo::petscComm(), "--------------- Dynamic problem --------------\n");
448 344 PetscPrintf(MpiInfo::petscComm(), "----------------------------------------------\n");
449
450
2/2
✓ Branch 1 taken 10880 times.
✓ Branch 2 taken 344 times.
 11224 while (!hasFinished()) {
451 10880 forward();
452 }
453 344 }
454
455 /***********************************************************************************/
456 /***********************************************************************************/
457
458 11721 bool Model::hasFinished()
459 {
460 11721 const auto& r_instance = FelisceParam::instance(this->instanceIndex());
461
4/4
✓ Branch 1 taken 11673 times.
✓ Branch 2 taken 48 times.
✓ Branch 4 taken 322 times.
✓ Branch 5 taken 11351 times.
 11721 m_hasFinished = ( m_fstransient->time > r_instance.timeMax || m_fstransient->iteration >= r_instance.timeIterationMax );
462
463 11721 return m_hasFinished;
464 }
465
466 /***********************************************************************************/
467 /***********************************************************************************/
468
469 ✗ double Model::getTime() const
470 {
471 ✗ return m_fstransient->time;
472 }
473
474 /***********************************************************************************/
475 /***********************************************************************************/
476
477 11544 void Model::printNewTimeIterationBanner()
478 {
479 11544 PetscPrintf(MpiInfo::petscComm(),"\n----------------------------------------------\n");
480 11544 PetscPrintf(MpiInfo::petscComm(),"%s Iteration: %d, Time: %e \n", m_name.c_str(), m_fstransient->iteration, m_fstransient->time);
481 11544 PetscPrintf(MpiInfo::petscComm(),"----------------------------------------------\n");
482 11544 }
483
484 /***********************************************************************************/
485 /***********************************************************************************/
486
487 ✗ void Model::setTime(const double time)
488 {
489 ✗ m_fstransient->time = time;
490 }
491
492 /***********************************************************************************/
493 /***********************************************************************************/
494
495 462 void Model::readMesh()
496 {
497 462 auto& r_instance = FelisceParam::instance(this->instanceIndex());
498 462 const std::size_t nbMesh = r_instance.inputMesh.size();
499 462 m_mesh.resize(nbMesh);
500 462 m_ios.resize(nbMesh);
501
502
2/2
✓ Branch 0 taken 468 times.
✓ Branch 1 taken 462 times.
 930 for(std::size_t imesh=0; imesh<nbMesh; ++imesh) {
503 468 auto& p_io = m_ios[imesh];
504 468 auto& r_mesh = m_mesh[imesh];
505 936 p_io = felisce::make_shared<IO>(r_instance.inputDirectory, r_instance.inputFile, r_instance.inputMesh[imesh],
506 468 r_instance.outputMesh[imesh], r_instance.meshDir, r_instance.resultDir,
507 936 r_instance.prefixName[imesh]);
508 468 r_mesh = felisce::make_shared<GeometricMeshRegion>();
509
510 468 p_io->readMesh(*r_mesh, r_instance.spaceUnit);
511
512
2/2
✓ Branch 1 taken 250 times.
✓ Branch 2 taken 218 times.
 468 if (FelisceParam::verbose()) {
513 250 PetscPrintf(MpiInfo::petscComm(),"\n/========Information about the input mesh========/.\n");
514 250 PetscPrintf(MpiInfo::petscComm(),"Name of the input mesh: <%s>.\n", r_instance.inputMesh[imesh].c_str());
515 }
516
517
2/2
✓ Branch 1 taken 131 times.
✓ Branch 2 taken 337 times.
 468 if (MpiInfo::rankProc() == 0) {
518 131 r_mesh->print(FelisceParam::verbose());
519 }
520
521 468 int& num_coor = r_instance.numCoor;
522
2/2
✓ Branch 0 taken 462 times.
✓ Branch 1 taken 6 times.
 468 if (num_coor == 0) {
523 462 num_coor = r_mesh->numCoor();
524
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 6 times.
 6 } else if (num_coor != r_mesh->numCoor()) {
525 ✗ FEL_ERROR("Mesh dimensions are inconsistent. Combining 2D and 3D meshes is not supported.");
526 }
527
528 // Check if the mesh is inverted
529
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 468 times.
 468 if (r_instance.checkVolumeOrientation) {
530 ✗ if (FEUtilities::CheckVolumeIsInverted(*r_mesh)) {
531 ✗ FEL_ERROR("Mesh is inverted. Check connectivity");
532 }
533 }
534 }
535
536 462 m_meshIsRead = true;
537 462 }
538
539 /***********************************************************************************/
540 /***********************************************************************************/
541
542 126 void Model::writeMesh()
543 {
544 // Getting instance
545 126 const auto& r_instance = FelisceParam::instance(this->instanceIndex());
546
547
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 124 times.
 126 if(FelisceParam::verbose()>1) {
548 2 PetscPrintf(MpiInfo::petscComm(),"Write mesh\n");
549 }
550
551
1/2
✓ Branch 0 taken 126 times.
✗ Branch 1 not taken.
 126 if(!r_instance.duplicateSupportDof) {
552
2/2
✓ Branch 1 taken 127 times.
✓ Branch 2 taken 126 times.
 253 for(std::size_t imesh=0; imesh<m_mesh.size(); ++imesh)
553
2/4
✓ Branch 4 taken 127 times.
✗ Branch 5 not taken.
✓ Branch 9 taken 127 times.
✗ Branch 10 not taken.
 127 m_ios[imesh]->writeMesh(*m_mesh[imesh]);
554 } else {
555 // because writeGeoForUnknown write the mesh from the support element not from the mesh
556 // This way, the duplicated support elements are saved and the mesh has the same dof as the solutions
557 ✗ for(std::size_t ipb=0; ipb<m_linearProblem.size(); ++ipb) {
558 ✗ for(std::size_t iUnknown=0; iUnknown<m_linearProblem[ipb]->listUnknown().size(); ++iUnknown) {
559 ✗ m_linearProblem[ipb]->writeGeoForUnknown(iUnknown);
560 }
561 }
562 }
563
564 126 m_meshIsWritten = true;
565 126 }
566
567 /***********************************************************************************/
568 /***********************************************************************************/
569
570 ✗ void Model::buildEdges()
571 {
572 ✗ for(std::size_t imesh=0; imesh<m_mesh.size(); ++imesh)
573 ✗ m_mesh[imesh]->buildEdges();
574 }
575
576 /***********************************************************************************/
577 /***********************************************************************************/
578
579 ✗ void Model::buildFaces()
580 {
581 ✗ for(std::size_t imesh=0; imesh<m_mesh.size(); ++imesh)
582 ✗ m_mesh[imesh]->buildFaces();
583 }
584
585 /***********************************************************************************/
586 /***********************************************************************************/
587
588 11983 void Model::writeSolution()
589 {
590 11983 const auto& r_instance = FelisceParam::instance(this->instanceIndex());
591
592
2/2
✓ Branch 1 taken 3062 times.
✓ Branch 2 taken 8921 times.
 11983 if (MpiInfo::rankProc() == 0) {
593
2/2
✓ Branch 0 taken 125 times.
✓ Branch 1 taken 2937 times.
 3062 if (!m_meshIsWritten) {
594 125 writeMesh();
595 }
596 }
597
598 11983 if( (m_fstransient->iteration % r_instance.frequencyWriteSolution == 0)
599 ✗ || m_hasFinished
600
2/8
✗ Branch 0 not taken.
✓ Branch 1 taken 11983 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 11983 times.
✗ Branch 11 not taken.
 11983 || ( m_fstransient->iteration >= r_instance.intervalWriteSolution[0] && m_fstransient->iteration <= r_instance.intervalWriteSolution[1] )) {
601
602
2/2
✓ Branch 1 taken 108 times.
✓ Branch 2 taken 11875 times.
 11983 if(FelisceParam::verbose() > 1)
603 108 PetscPrintf(MpiInfo::petscComm(),"Write solutions\n");
604
605
2/2
✓ Branch 1 taken 12825 times.
✓ Branch 2 taken 11983 times.
 24808 for (std::size_t ipb = 0; ipb < m_linearProblem.size(); ipb++) {
606 12825 m_linearProblem[ipb]->writeSolution(MpiInfo::rankProc(), m_ios, m_fstransient->time, m_fstransient->iteration);
607 }
608
609
2/2
✓ Branch 1 taken 3062 times.
✓ Branch 2 taken 8921 times.
 11983 if (MpiInfo::rankProc() == 0) {
610
2/2
✓ Branch 1 taken 3066 times.
✓ Branch 2 taken 3062 times.
 6128 for(std::size_t iio=0; iio<m_ios.size(); ++iio) {
611
1/2
✓ Branch 2 taken 3066 times.
✗ Branch 3 not taken.
 3066 if ( m_ios[iio]->typeOutput == 1 ) // 1 : ensight{
612 3066 m_ios[iio]->postProcess(m_fstransient->time/*, !m_sameGeometricMeshForVariable*/);
613 }
614 }
615 }
616
617 // Warning. If using a bdf for time derivative resolution, you have to modify this call with your bdf method order.
618
2/2
✓ Branch 1 taken 12825 times.
✓ Branch 2 taken 11983 times.
 24808 for (std::size_t ipb = 0; ipb < m_linearProblem.size(); ipb++) {
619
1/2
✗ Branch 3 not taken.
✓ Branch 4 taken 12825 times.
 12825 if(m_solutionBackup[ipb].hasToBackup(m_fstransient->iteration, r_instance.frequencyBackup) ) {
620 ✗ if(FelisceParam::verbose() > 1)
621 ✗ PetscPrintf(MpiInfo::petscComm(),"Backup the solutions for linear problem %ld\n", static_cast<unsigned long>(ipb));
622
623 // Backup of linear problem variables
624 ✗ m_linearProblem[ipb]->writeSolutionBackup(m_initialCondition.listVariable(), MpiInfo::rankProc(), m_solutionBackup[ipb].listIO(), m_fstransient->time, m_fstransient->iteration);
625
626 // This part has to be improved (07/05/13, Elisa) TODO
627
628 // Backup of linear problem meshes
629 ✗ m_solutionBackup[ipb].backup(m_fstransient->time, m_mesh, MpiInfo::rankProc());
630 }
631 }
632 11983 }
633
634 /***********************************************************************************/
635 /***********************************************************************************/
636
637 ✗ void Model::ComputeStress(int coef, int label, std::vector<double>& stress)
638 {
639 ✗ m_linearProblem[0]->computeBoundaryStress(coef,label,stress);
640 }
641
642 /***********************************************************************************/
643 /***********************************************************************************/
644
645 ✗ void Model::addCallbackXYZ(const char *name, const CallbackXYZ& cb)
646 {
647 ✗ m_mapCallbackXYZ[name] = cb;
648 }
649
650 /***********************************************************************************/
651 /***********************************************************************************/
652
653 8 void Model::addCallbackXYZT(const char* name, const CallbackXYZT& cb)
654 {
655
3/6
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 8 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 8 times.
✗ Branch 9 not taken.
 8 m_mapCallbackXYZT[name] = cb;
656
657 //give the callback to each ElementFieldDynamicValue that want it
658 Component comp;
659 ElementFieldDynamicValue *efDV;
660 // iterate on all elementFieldDynamicValue contained in the std::unordered_map
661
1/2
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
 8 for (auto it = FelisceParam::instance(this->instanceIndex()).elementFieldDynamicValueMap.begin() ;
662
3/4
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 8 times.
✓ Branch 7 taken 8 times.
 16 it != FelisceParam::instance(this->instanceIndex()).elementFieldDynamicValueMap.end() ;
663 8 it++) {
664 8 efDV = it->second;
665 // iterate on all components
666
3/4
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 8 times.
 16 for (int icomp=0 ; icomp < efDV->numComp() ; ++icomp) {
667 8 comp = efDV->allComp[icomp];
668
669 // if the component has to be std::set with the function passed to
670 // Model::addCallbackXYZT, give it to the ElementFieldDynamicValue
671
11/24
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 8 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 8 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 4 times.
✓ Branch 14 taken 4 times.
✓ Branch 15 taken 8 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 8 times.
✗ Branch 19 not taken.
✓ Branch 21 taken 8 times.
✗ Branch 22 not taken.
✓ Branch 24 taken 4 times.
✓ Branch 25 taken 4 times.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
 8 if ( efDV->typeValueOfElementField(comp) == FROM_FUNCTION && efDV->functionName(comp) == static_cast<std::string>(name)) {
672
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 efDV->setCallbackXYZT(cb,comp);
673 }
674 }
675 }
676 8 }
677
678 /***********************************************************************************/
679 /***********************************************************************************/
680
681 ✗ CallbackXYZ& Model::callbackXYZ(const char* name)
682 {
683 ✗ auto it = m_mapCallbackXYZ.find(name);
684 ✗ auto end = m_mapCallbackXYZ.end();
685
686 ✗ if (it == end) {
687 ✗ const int n = m_mapCallbackXYZ.size();
688 ✗ std::ostringstream msg;
689 ✗ msg << "No such callback '" << name << "'. ";
690 ✗ msg << n << (n>1 ? " candidates" : " candidate") << (n>0 ? ": (" : ".");
691 ✗ std::size_t counter = 1;
692 ✗ for (it = m_mapCallbackXYZ.begin() ; it != end ; it++) {
693 ✗ ++counter;
694 ✗ msg << "'" << it->first << "'" << (counter == m_mapCallbackXYZ.size() ? ")." : ", ");
695 }
696 ✗ FEL_ERROR(msg.str().c_str());
697 }
698
699 ✗ return it->second;
700 }
701
702 /***********************************************************************************/
703 /***********************************************************************************/
704
705 ✗ CallbackXYZT& Model::callbackXYZT(const char* name)
706 {
707 ✗ auto it = m_mapCallbackXYZT.find(name);
708 ✗ auto end = m_mapCallbackXYZT.end();
709
710 ✗ if (it == end) {
711 ✗ int n = m_mapCallbackXYZT.size();
712 ✗ std::ostringstream msg;
713 ✗ msg << "No such callback '" << name << "'. ";
714 ✗ msg << n << (n>1 ? " candidates" : " candidate") << (n>0 ? ": (" : ".");
715 ✗ std::size_t counter = 1;
716 ✗ for (it = m_mapCallbackXYZT.begin() ; it != end; it++) {
717 ✗ ++counter;
718 ✗ msg << "'" << it->first << "'" << (counter == m_mapCallbackXYZT.size() ? ")." : ", ");
719 }
720 ✗ FEL_ERROR(msg.str().c_str());
721 }
722
723 ✗ return it->second;
724 }
725
726 /***********************************************************************************/
727 /***********************************************************************************/
728
729 ✗ double Model::evaluateCallback(const char* name, double x, double y, double z)
730 {
731 ✗ CallbackXYZ cb = callbackXYZ(name);
732 ✗ return cb(x,y,z);
733 }
734
735 /***********************************************************************************/
736 /***********************************************************************************/
737
738 ✗ double Model::evaluateCallback(const char* name, double x, double y, double z, double t)
739 {
740 ✗ CallbackXYZT cb = callbackXYZT(name);
741 ✗ return cb(x,y,z,t);
742 }
743
744 /***********************************************************************************/
745 /***********************************************************************************/
746
747 266 void Model::setInitialCondition()
748 {
749 266 const auto& r_instance = FelisceParam::instance(this->instanceIndex());
750
751
2/2
✓ Branch 1 taken 299 times.
✓ Branch 2 taken 266 times.
 565 for (std::size_t iProblem=0; iProblem<m_linearProblem.size(); iProblem++) {
752
1/2
✓ Branch 0 taken 299 times.
✗ Branch 1 not taken.
 299 if (!r_instance.restartSolution) {
753 299 m_U_0.duplicateFrom(m_linearProblem[iProblem]->vector());
754
1/2
✓ Branch 1 taken 299 times.
✗ Branch 2 not taken.
 299 m_U_0.set(0.0);
755
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 299 times.
 299 if (r_instance.hasInitialCondition) {
756
757 ✗ int idVar = -1;
758 double valueByDofU_0;
759 ✗ felInt aa = 0;
760
761 ✗ for (std::size_t iListVar = 0; iListVar < m_initialCondition.listVariable().size(); iListVar++) {
762 ✗ for ( std::size_t iUnknown = 0; iUnknown < m_linearProblem[iProblem]->listUnknown().size(); iUnknown++) {
763 ✗ idVar = m_linearProblem[iProblem]->listUnknown().idVariable(iUnknown);
764 ✗ if ( m_linearProblem[iProblem]->listVariable()[idVar].name() == m_initialCondition.listVariable()[iListVar].name() ) {
765 ✗ felInt numDofLpU = m_linearProblem[iProblem]->numDofLocalPerUnknown(m_linearProblem[iProblem]->listUnknown()[iUnknown]);
766 ✗ felInt idLocalValue[numDofLpU];
767 ✗ felInt idGlobalValue[numDofLpU];
768 ✗ for ( felInt ii = 0; ii < numDofLpU; ii++) {
769 ✗ idLocalValue[ii] = ii;
770 }
771 ✗ ISLocalToGlobalMappingApply(m_linearProblem[iProblem]->mappingDofLocalToDofGlobal(m_linearProblem[iProblem]->listUnknown()[iUnknown]),numDofLpU,&idLocalValue[0],&idGlobalValue[0]);
772 ✗ for (felInt ii = 0; ii < numDofLpU; ii++) {
773
774 // TODO this part of the function is not defined for a generic Model... maybe it should be moved or rewritten
775 ✗ if(r_instance.readDisplFromFile && r_instance.useALEformulation) {
776 ✗ aa = idGlobalValue[ii];
777 ✗ static_cast<LinearProblemHarmonicExtension*>(m_linearProblem[1])->HarmExtSol().getValues(1, &aa, &valueByDofU_0);
778 } else {
779 ✗ valueByDofU_0 = r_instance.valueInitCond[iListVar];
780 }
781 ✗ m_U_0.setValue(idGlobalValue[ii],valueByDofU_0,INSERT_VALUES);
782 }
783 }
784 }
785 }
786 }
787 299 m_U_0.assembly();
788 } else {
789 ✗ m_initialCondition.initializeFromFile(r_instance.restartSolutionDir, r_instance.prefixName[m_linearProblem[iProblem]->currentMesh()] + "_bkp.case");
790 ✗ userDefinedInitialConditions(iProblem);
791 }
792
793 //Initialize solution for solver.
794 299 m_linearProblem[iProblem]->solution().copyFrom(m_U_0);
795 }
796 266 }
797
798 /***********************************************************************************/
799 /***********************************************************************************/
800
801 ✗ void Model::writeMatrixForMatlab(const std::string& filename, unsigned int index) const
802 {
803 ✗ m_linearProblem[index]->writeMatrixForMatlab(std::numeric_limits<int>::max(), filename);
804 }
805
806 /***********************************************************************************/
807 /***********************************************************************************/
808
809 ✗ void Model::writeRHSForMatlab(const std::string& filename, unsigned int index) const
810 {
811 ✗ m_linearProblem[index]->writeRHSForMatlab(std::numeric_limits<int>::max(), filename);
812 }
813
814 /***********************************************************************************/
815 /***********************************************************************************/
816
817 #ifdef FELISCE_WITH_CVGRAPH
818
819 ✗ void Model::startIterationCVG()
820 {
821 ✗ cvgMainSlave* slave=m_linearProblem[0]->slave();
822 ✗ if (slave==nullptr) {
823 ✗ FEL_ERROR("slave not initialized");
824 }
825 ✗ if( slave->newTimeStep() ) {
826 ✗ std::stringstream msg;
827 ✗ msg<<"Starting a new time step"<<std::endl;
828 ✗ PetscPrintf(MpiInfo::petscComm(),"%s",msg.str().c_str());
829 ✗ cvgraphNewTimeStep();
830 ✗ m_rhsWithoutBC_zeroed=false;
831 ✗ } else if( slave->redoTimeStep() ){
832 ✗ std::stringstream msg;
833 ✗ msg << "......................................" << std::endl;
834 ✗ msg << "Redo time step at time t= " << m_fstransient->time << std::endl;
835 ✗ PetscPrintf(MpiInfo::petscComm(),"%s",msg.str().c_str());
836 ✗ rhsWithoutBC_restore();
837 ✗ m_linearProblem[0]->matrix(0).copyFrom(m_linearProblem[0]->matrixWithoutBC(),SAME_NONZERO_PATTERN);
838 ✗ m_linearProblem[0]->vector().copyFrom(m_linearProblem[0]->rhsWithoutBC());
839 ✗ } else if( slave->redoTimeStepOnlyCVGData() ) {
840 ✗ std::stringstream msg;
841 ✗ msg << "......................................" << std::endl;
842 ✗ msg << "Redo time step, with only cvg data, at time t= " << m_fstransient->time << std::endl;
843 ✗ PetscPrintf(MpiInfo::petscComm(),"%s",msg.str().c_str());
844 ✗ rhsWithoutBC_zeroEntries();
845 ✗ m_linearProblem[0]->matrix(0).copyFrom(m_linearProblem[0]->matrixWithoutBC(),SAME_NONZERO_PATTERN);
846 ✗ m_linearProblem[0]->vector().zeroEntries();
847 ✗ m_linearProblem[0]->applyBC(FelisceParam::instance(this->instanceIndex()).essentialBoundaryConditionsMethod, MpiInfo::rankProc(), FlagMatrixRHS::only_matrix);//We insert dirichlet and robin parts to the matrix.
848 ✗ } else {
849 ✗ slave->msgInfo().print(1);
850 ✗ FEL_ERROR("Error: unexpected timeStatus");
851 }
852 }
853
854 /***********************************************************************************/
855 /***********************************************************************************/
856
857 ✗ void Model::ImposeCVGraphInterfaceSolution(double* ImpValue)
858 {
859 ✗ m_cvgraphInterf->ImposeSolution(ImpValue);
860 }
861
862 /***********************************************************************************/
863 /***********************************************************************************/
864
865 ✗ void Model::ImposeCVGraphInterfaceSolutionAndComputedVariable(double* ImpValue, int idVariable)
866 {
867 ✗ m_cvgraphInterf->ImposeSolutionAndComputedVariable(ImpValue, idVariable);
868 }
869
870 /***********************************************************************************/
871 /***********************************************************************************/
872
873 ✗ void Model::ExtractCVGraphInterfaceSolution(double* valueDofInterface)
874 {
875 ✗ m_cvgraphInterf->ExtractSolution(m_linearProblem[0]->solution(), valueDofInterface);
876 }
877
878 /***********************************************************************************/
879 /***********************************************************************************/
880
881 ✗ void Model::ExtractCVGraphInterfaceSolution(double* valueDofInterface, int idVariable)
882 {
883 ✗ m_cvgraphInterf->ExtractSolution( m_linearProblem[0]->solution(), valueDofInterface, idVariable);
884 }
885
886 /***********************************************************************************/
887 /***********************************************************************************/
888
889 ✗ void Model::ExtractCVGraphInterfaceSolutionAndComputedVariable(double* valueDofInterface,int idVariable, std::vector<double>& interfaceComputedVariable)
890 {
891 ✗ m_cvgraphInterf->ExtractSolutionAndComputedVariable( m_linearProblem[0]->solution(), interfaceComputedVariable, valueDofInterface, idVariable);
892 }
893
894 /***********************************************************************************/
895 /***********************************************************************************/
896
897 ✗ void Model::StoreCVGraphInterfaceSolution(double* StoringValue)
898 {
899 ✗ m_cvgraphInterf->StoreSolution(StoringValue);
900 }
901
902 /***********************************************************************************/
903 /***********************************************************************************/
904
905 ✗ void Model::rhsWithoutBC_zeroEntries()
906 {
907 // If the m_rhsWithoutBC was already zeroed we do nothing
908 ✗ if ( !m_rhsWithoutBC_zeroed ) {
909
910 // We save a copy of it before erasing it.
911 // So that we can restore it later
912 ✗ if ( m_rhsWithForceTermWithoutBC.isNull() ) {
913 // In case it was not initialized we do it now.
914 ✗ m_rhsWithForceTermWithoutBC.duplicateFrom(m_linearProblem[0]->rhsWithoutBC());
915 }
916 ✗ m_rhsWithForceTermWithoutBC.copyFrom(m_linearProblem[0]->rhsWithoutBC());
917
918 // We std::set it to zero.
919 ✗ m_linearProblem[0]->rhsWithoutBC().zeroEntries();
920 ✗ m_rhsWithoutBC_zeroed=true;
921 }
922 }
923
924 /***********************************************************************************/
925 /***********************************************************************************/
926
927 ✗ void Model::rhsWithoutBC_restore()
928 {
929 // If it was zeroed, we restore it.
930 ✗ if (m_rhsWithoutBC_zeroed) {
931 ✗ m_linearProblem[0]->rhsWithoutBC().copyFrom(m_rhsWithForceTermWithoutBC);
932 }
933 }
934
935 #endif
936 }
937