GCC Code Coverage Report


Directory: ./
File: Solver/linearProblem_WriteMethods.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 180 425 42.4%
Branches: 157 780 20.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 & J-F. Gerbeau & V. Martin
13 //
14
15 // System includes
16
17 // External includes
18
19 // Project includes
20 #include "InputOutput/io.hpp"
21 #include "Solver/linearProblem.hpp"
22
23 namespace felisce
24 {
25 // Print solution of the system solve.
26 14901 void LinearProblem::printSolution(int verbose, std::ostream& outstr) const {
27 IGNORE_UNUSED_OUTSTR;
28
29
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 14901 times.
14901 if (verbose > 19) {
30 PetscPrintf(MpiInfo::petscComm(),"\n/================ Solution of the assembly system Ax=b ================/\n");
31 m_sol.view();
32 }
33 14901 }
34
35 //===============================================================
36 // MATLAB FUNCTIONS FOR DEBUG
37 // It is possible to save the objects also in binary format.
38 // The functions to do that are in the wrappers class: petscVector and petscMatrix.
39 // Also the use of this function is not always useful since you
40 // can directly call the function of the wrapper.
41 //
42 // Write in file matrix.m, matrix in matlab format.
43 33965 void LinearProblem::writeMatrixForMatlab(int verbose, const std::string& filename) const {
44
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 33965 times.
33965 if (verbose > 20) {
45 m_matrices[0].saveInMatlabFormat(MpiInfo::petscComm(),filename.c_str(),FelisceParam::instance().resultDir);
46 }
47 33965 }
48 //Write in file RHS.m, RHS in matlab format.
49 33965 void LinearProblem::writeRHSForMatlab(int verbose, const std::string& filename) const {
50
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 33965 times.
33965 if (verbose > 20) {
51 vector().saveInMatlabFormat(MpiInfo::petscComm(),filename,FelisceParam::instance().resultDir);
52 }
53 33965 }
54
55 //Write in file sol.m, RHS in matlab format.
56 void LinearProblem::writeSolForMatlab(int verbose, const std::string& filename) const {
57 if (verbose > 20) {
58 m_sol.saveInMatlabFormat(MpiInfo::petscComm(),filename,FelisceParam::instance().resultDir);
59 }
60 }
61
62 //Write in file filename, a matrix in matlab format.
63 void LinearProblem::writeMatrixForMatlab(std::string const& fileName, PetscMatrix& matrix) const {
64 matrix.saveInMatlabFormat(MpiInfo::petscComm(),fileName.c_str(),FelisceParam::instance().resultDir);
65 }
66
67 //Write in file filename, a std::vector in matlab format.
68 void LinearProblem::writeVectorForMatlab(std::string const& fileName, PetscVector& vector) const {
69 vector.saveInMatlabFormat(MpiInfo::petscComm(),fileName.c_str(),FelisceParam::instance().resultDir);
70 }
71
72 // Write in file matrix.mb, matrix in matlab format.
73 void LinearProblem::writeMatrixForMatlabBinary(int verbose, const std::string& filename) const {
74 if (verbose > 20) {
75 m_matrices[0].saveInBinaryFormat(MpiInfo::petscComm(),filename.c_str(),FelisceParam::instance().resultDir);
76 }
77 }
78 //Write in file RHS.m, RHS in matlab format.
79 void LinearProblem::writeRHSForMatlabBinary(int verbose, const std::string& filename) const {
80 if (verbose > 20) {
81 vector().saveInBinaryFormat(MpiInfo::petscComm(),filename,FelisceParam::instance().resultDir);
82 }
83 }
84
85 //Write in file sol.mb, RHS in matlab format.
86 void LinearProblem::writeSolForMatlabBinary(int verbose, const std::string& filename) const {
87 if (verbose > 20) {
88 m_sol.saveInBinaryFormat(MpiInfo::petscComm(),filename,FelisceParam::instance().resultDir);
89 }
90 }
91
92 //Write in file filename, a matrix in matlab format.
93 void LinearProblem::writeMatrixForMatlabBinary(std::string const& fileName, PetscMatrix& matrix) const {
94 matrix.saveInBinaryFormat(MpiInfo::petscComm(),fileName.c_str(),FelisceParam::instance().resultDir);
95 }
96
97 //Write in file filename, a std::vector in matlab format.
98 void LinearProblem::writeVectorForMatlabBinary(std::string const& fileName, PetscVector& vector) const {
99 vector.saveInBinaryFormat(MpiInfo::petscComm(),fileName.c_str(),FelisceParam::instance().resultDir);
100 }
101 //================================================================
102
103 12869 void LinearProblem::writeSolution(int rank, std::vector<IO::Pointer>& io, double& time, int iteration) {
104 // to be sure that m_seqSol is up to date
105
1/2
✓ Branch 1 taken 12869 times.
✗ Branch 2 not taken.
12869 gatherSolution();
106
2/2
✓ Branch 0 taken 9526 times.
✓ Branch 1 taken 3343 times.
12869 if ( rank != 0)
107 9526 return;
108
109
1/2
✓ Branch 2 taken 3343 times.
✗ Branch 3 not taken.
3343 std::vector<double> valueVec(m_numDof);
110
1/2
✓ Branch 1 taken 3343 times.
✗ Branch 2 not taken.
3343 m_seqSol.getAllValuesInAppOrdering(m_ao, valueVec);
111
112 // Write solution in ensight format
113 int idVar, idVar2, iMesh;
114 double* solution;
115
2/2
✓ Branch 1 taken 5164 times.
✓ Branch 2 taken 3343 times.
8507 for ( std::size_t iUnknown = 0; iUnknown < m_listUnknown.size(); iUnknown++) {
116 5164 felInt numValueStart = 0;
117 5164 idVar = m_listUnknown.idVariable(iUnknown);
118 5164 iMesh = m_listVariable[idVar].idMesh();
119
120
2/2
✓ Branch 0 taken 1839 times.
✓ Branch 1 taken 5164 times.
7003 for ( std::size_t iUnknown2 = 0; iUnknown2 < iUnknown; iUnknown2++) {
121 1839 idVar2 = m_listUnknown.idVariable(iUnknown2);
122
123 // the number of the previous suppDofs could be different from the number of nodes
124 //if NOT using FusionDof numSupportDof() == listNode().size()
125 //if DO using FusionDof numSupportDof() != listNode().size()
126 1839 numValueStart += m_supportDofUnknown[iUnknown2].numSupportDof()*m_listVariable[idVar2].numComponent();
127 }
128
129
3/4
✓ Branch 2 taken 1508 times.
✓ Branch 3 taken 956 times.
✓ Branch 4 taken 2700 times.
✗ Branch 5 not taken.
5164 switch (m_listVariable[idVar].numComponent()) {
130 1508 case 1:
131
2/4
✓ Branch 3 taken 1508 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 1508 times.
✗ Branch 8 not taken.
1508 solution = new double[m_supportDofUnknown[iUnknown].listNode().size()];
132
2/2
✓ Branch 3 taken 1931474 times.
✓ Branch 4 taken 1508 times.
1932982 for (felInt i = 0; i < static_cast<felInt>(m_supportDofUnknown[iUnknown].listNode().size()); i++) {
133
134
3/4
✓ Branch 1 taken 1931474 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 89610 times.
✓ Branch 4 taken 1841864 times.
1931474 if (FelisceParam::instance().FusionDof) {
135 89610 solution[i] = valueVec[m_supportDofUnknown[iUnknown].listPerm()[i]+numValueStart];
136 } else {
137 1841864 solution[i] = valueVec[i+numValueStart];
138 }
139 }
140
2/4
✓ Branch 8 taken 1508 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 1508 times.
✗ Branch 12 not taken.
1508 io[iMesh]->writeSolution(rank, time, iteration, 0, m_listVariable[idVar].name(), solution, m_supportDofUnknown[iUnknown].listNode().size(), dimension(), true);
141 1508 break;
142 956 case 2:
143
2/4
✓ Branch 3 taken 956 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 956 times.
✗ Branch 8 not taken.
956 solution = new double[m_supportDofUnknown[iUnknown].listNode().size()*3];
144
2/2
✓ Branch 3 taken 741795 times.
✓ Branch 4 taken 956 times.
742751 for (felInt i = 0; i < static_cast<felInt>(m_supportDofUnknown[iUnknown].listNode().size()); i++) {
145
146
3/4
✓ Branch 1 taken 741795 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 78390 times.
✓ Branch 4 taken 663405 times.
741795 if (FelisceParam::instance().FusionDof) {
147 78390 solution[i*3] = valueVec[m_supportDofUnknown[iUnknown].listPerm()[i]*2+numValueStart];
148 78390 solution[i*3+1] = valueVec[m_supportDofUnknown[iUnknown].listPerm()[i]*2+1+numValueStart];
149 78390 solution[i*3+2] = 0.;
150 } else {
151 663405 solution[i*3] = valueVec[i*2+numValueStart];
152 663405 solution[i*3+1] = valueVec[i*2+1+numValueStart];
153 663405 solution[i*3+2] = 0.;
154 }
155 }
156
2/4
✓ Branch 8 taken 956 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 956 times.
✗ Branch 12 not taken.
956 io[iMesh]->writeSolution(rank, time, iteration, 1, m_listVariable[idVar].name(), solution, m_supportDofUnknown[iUnknown].listNode().size(), dimension(),true);
157 956 break;
158 2700 case 3:
159
2/4
✓ Branch 3 taken 2700 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 2700 times.
✗ Branch 8 not taken.
2700 solution = new double[m_supportDofUnknown[iUnknown].listNode().size()*3];
160
2/2
✓ Branch 3 taken 2543217 times.
✓ Branch 4 taken 2700 times.
2545917 for (felInt i = 0; i < static_cast<felInt>(m_supportDofUnknown[iUnknown].listNode().size()); i++) {
161
162
3/4
✓ Branch 1 taken 2543217 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 22440 times.
✓ Branch 4 taken 2520777 times.
2543217 if (FelisceParam::instance().FusionDof) {
163 22440 solution[i*3] = valueVec[m_supportDofUnknown[iUnknown].listPerm()[i]*3+numValueStart];
164 22440 solution[i*3+1] = valueVec[m_supportDofUnknown[iUnknown].listPerm()[i]*3+1+numValueStart];
165 22440 solution[i*3+2] = valueVec[m_supportDofUnknown[iUnknown].listPerm()[i]*3+2+numValueStart];
166 } else {
167 2520777 solution[i*3] = valueVec[i*3+numValueStart];
168 2520777 solution[i*3+1] = valueVec[i*3+1+numValueStart];
169 2520777 solution[i*3+2] = valueVec[i*3+2+numValueStart];
170 }
171 }
172
2/4
✓ Branch 8 taken 2700 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 2700 times.
✗ Branch 12 not taken.
2700 io[iMesh]->writeSolution(rank, time, iteration, 1, m_listVariable[idVar].name(), solution, m_supportDofUnknown[iUnknown].listNode().size(), dimension(),true);
173 2700 break;
174 default:
175 FEL_ERROR("Problem with the dimension of your solution to write it.")
176 break;
177 }
178 }
179 // Do not need to delete here solution because it is deleted in EnsightCase::postProcess (manageMemory of io.writeSolution = true)
180 //delete [] solution;
181 3343 }
182
183 392 void LinearProblem::writeSolutionFromVec(PetscVector& v, int rank, std::vector<IO::Pointer>& io, double& time, int iteration, std::string prefix) {
184 // Create a sequential vector with the size of sol
185
1/2
✓ Branch 1 taken 392 times.
✗ Branch 2 not taken.
392 PetscVector xsol;
186
1/2
✓ Branch 1 taken 392 times.
✗ Branch 2 not taken.
392 v.scatterToZero(xsol, INSERT_VALUES, SCATTER_FORWARD);
187
2/2
✓ Branch 0 taken 294 times.
✓ Branch 1 taken 98 times.
392 if ( rank != 0) {
188 294 return;
189 }
190
191
1/2
✓ Branch 2 taken 98 times.
✗ Branch 3 not taken.
98 std::vector<double> valueVec( m_numDof );
192
1/2
✓ Branch 1 taken 98 times.
✗ Branch 2 not taken.
98 xsol.getAllValuesInAppOrdering( m_ao, valueVec );
193
194 98 int idVar = -1;
195 98 int idVar2 = -1;
196 98 int iMesh = -1;
197 // Write solution in ensight format
198 double* solution;
199
2/2
✓ Branch 1 taken 196 times.
✓ Branch 2 taken 98 times.
294 for ( std::size_t iUnknown = 0; iUnknown < m_listUnknown.size(); iUnknown++) {
200 196 felInt numValueStart = 0;
201 196 idVar = m_listUnknown.idVariable(iUnknown);
202 196 iMesh = m_listVariable[idVar].idMesh();
203
2/2
✓ Branch 0 taken 98 times.
✓ Branch 1 taken 196 times.
294 for ( std::size_t iUnknown2 = 0; iUnknown2 < iUnknown; iUnknown2++) {
204 98 idVar2 = m_listUnknown.idVariable(iUnknown2);
205 // the number of the previous suppDofs could be different from the number of nodes
206 //if NOT using FusionDof numSupportDof() == listNode().size()
207 //if DO using FusionDof numSupportDof() != listNode().size()
208 98 numValueStart += supportDofUnknown(iUnknown2).numSupportDof()*m_listVariable[idVar2].numComponent();
209 }
210
1/2
✓ Branch 1 taken 196 times.
✗ Branch 2 not taken.
196 std::string name=prefix;
211
212
3/4
✓ Branch 2 taken 98 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 90 times.
✗ Branch 5 not taken.
196 switch (m_listVariable[idVar].numComponent()) {
213
214 98 case 1:
215
2/4
✓ Branch 3 taken 98 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 98 times.
✗ Branch 8 not taken.
98 solution = new double[supportDofUnknown(iUnknown).listNode().size()];
216
2/2
✓ Branch 3 taken 1214462 times.
✓ Branch 4 taken 98 times.
1214560 for (felInt i = 0; i < static_cast<felInt>(supportDofUnknown(iUnknown).listNode().size()); i++) {
217
218
2/4
✓ Branch 1 taken 1214462 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1214462 times.
1214462 if (FelisceParam::instance().FusionDof) {
219 solution[i] = valueVec[supportDofUnknown(iUnknown).listPerm()[i]+numValueStart];
220 } else {
221 1214462 solution[i] = valueVec[i+numValueStart];
222 }
223 }
224
2/4
✓ Branch 2 taken 98 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 98 times.
✗ Branch 6 not taken.
98 name += std::string(m_listVariable[idVar].name());
225
1/2
✓ Branch 7 taken 98 times.
✗ Branch 8 not taken.
98 io[iMesh]->writeSolution(rank, time, iteration, 0,name, solution, supportDofUnknown(iUnknown).listNode().size(), dimension(),true);
226 98 break;
227
228 8 case 2:
229
2/4
✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 8 times.
✗ Branch 8 not taken.
8 solution = new double[supportDofUnknown(iUnknown).listNode().size()*3];
230
2/2
✓ Branch 3 taken 7616 times.
✓ Branch 4 taken 8 times.
7624 for (felInt i = 0; i < static_cast<felInt>(supportDofUnknown(iUnknown).listNode().size()); i++) {
231
232
2/4
✓ Branch 1 taken 7616 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 7616 times.
7616 if (FelisceParam::instance().FusionDof) {
233 solution[i*3] = valueVec[supportDofUnknown(iUnknown).listPerm()[i]*2+numValueStart];
234 solution[i*3+1] = valueVec[supportDofUnknown(iUnknown).listPerm()[i]*2+1+numValueStart];
235 solution[i*3+2] = 0.;
236 } else {
237 7616 solution[i*3] = valueVec[i*2+numValueStart];
238 7616 solution[i*3+1] = valueVec[i*2+1+numValueStart];
239 7616 solution[i*3+2] = 0.;
240 }
241 }
242
2/4
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 8 times.
✗ Branch 6 not taken.
8 name += std::string(m_listVariable[idVar].name());
243
1/2
✓ Branch 7 taken 8 times.
✗ Branch 8 not taken.
8 io[iMesh]->writeSolution(rank, time, iteration, 1,name, solution, supportDofUnknown(iUnknown).listNode().size(), dimension(),true);
244 8 break;
245
246 90 case 3:
247
2/4
✓ Branch 3 taken 90 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 90 times.
✗ Branch 8 not taken.
90 solution = new double[supportDofUnknown(iUnknown).listNode().size()*3];
248
2/2
✓ Branch 3 taken 1206846 times.
✓ Branch 4 taken 90 times.
1206936 for (felInt i = 0; i < static_cast<felInt>(supportDofUnknown(iUnknown).listNode().size()); i++) {
249
250
2/4
✓ Branch 1 taken 1206846 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1206846 times.
1206846 if (FelisceParam::instance().FusionDof) {
251 solution[i*3] = valueVec[supportDofUnknown(iUnknown).listPerm()[i]*3+numValueStart];
252 solution[i*3+1] = valueVec[supportDofUnknown(iUnknown).listPerm()[i]*3+1+numValueStart];
253 solution[i*3+2] = valueVec[supportDofUnknown(iUnknown).listPerm()[i]*3+2+numValueStart];
254 } else {
255 1206846 solution[i*3] = valueVec[i*3+numValueStart];
256 1206846 solution[i*3+1] = valueVec[i*3+1+numValueStart];
257 1206846 solution[i*3+2] = valueVec[i*3+2+numValueStart];
258 }
259 }
260
2/4
✓ Branch 2 taken 90 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 90 times.
✗ Branch 6 not taken.
90 name += std::string(m_listVariable[idVar].name());
261 // in write solution there is a delete [] solution
262
1/2
✓ Branch 7 taken 90 times.
✗ Branch 8 not taken.
90 io[iMesh]->writeSolution(rank, time, iteration, 1,name, solution, supportDofUnknown(iUnknown).listNode().size(), dimension(),true);
263 90 break;
264
265 default:
266 FEL_ERROR("Problem with the dimension of your solution to write it.")
267 break;
268 }
269 196 }
270
2/2
✓ Branch 2 taken 98 times.
✓ Branch 3 taken 294 times.
392 }
271
272 /*!
273 \brief write solution for backup in Ensight format.
274 */
275 void LinearProblem::writeSolutionBackup(const ListVariable& listVar, int rank, std::vector<IO::Pointer>& io, double& time, int iteration) {
276 gatherSolution();
277
278 std::vector<double> valueVec(m_numDof);
279 m_seqSol.getAllValuesInAppOrdering( m_ao, valueVec );
280
281 std::vector<std::size_t> printingVar;
282
283 int idVar = -1;
284 int idVar2 = -1;
285 int iMesh = -1;
286 // Write solution in ensight format
287 if ( rank == 0) {
288 double* solution;
289 for ( std::size_t iUnknown = 0; iUnknown < m_listUnknown.size(); iUnknown++) {
290 felInt numValueStart = 0;
291 idVar = m_listUnknown.idVariable(iUnknown);
292 iMesh = m_listVariable[idVar].idMesh();
293 for (std::size_t iListVar = 0; iListVar < listVar.size(); iListVar++) {
294 if ( m_listVariable[idVar].name() == listVar[iListVar].name() ) {
295 for ( std::size_t iUnknown2 = 0; iUnknown2 < iUnknown; iUnknown2++) {
296 idVar2 = m_listUnknown.idVariable(iUnknown2);
297 numValueStart += supportDofUnknown(iUnknown2).listNode().size()*m_listVariable[idVar2].numComponent();
298 }
299
300 switch (m_listVariable[idVar].numComponent()) {
301
302 case 1:
303 solution = new double[supportDofUnknown(iUnknown).listNode().size()];
304 for (felInt i = 0; i < static_cast<felInt>(supportDofUnknown(iUnknown).listNode().size()); i++) {
305 solution[i] = valueVec[i+numValueStart];
306 }
307 io[iMesh]->writeSolution(rank, time, iteration, 0, m_listVariable[idVar].name()+"_bkp", solution, supportDofUnknown(iUnknown).listNode().size(), dimension(),true);
308 break;
309
310 case 2:
311 solution = new double[supportDofUnknown(iUnknown).listNode().size()*3];
312 for (felInt i = 0; i < static_cast<felInt>(supportDofUnknown(iUnknown).listNode().size()); i++) {
313 solution[i*3] = valueVec[i*2+numValueStart];
314 solution[i*3+1] = valueVec[i*2+1+numValueStart];
315 solution[i*3+2] = 0.;
316 }
317 io[iMesh]->writeSolution(rank, time, iteration, 1, m_listVariable[idVar].name()+"_bkp", solution, supportDofUnknown(iUnknown).listNode().size(), dimension(),true);
318 break;
319
320 case 3:
321 solution = new double[supportDofUnknown(iUnknown).listNode().size()*3];
322 for (felInt i = 0; i < static_cast<felInt>(supportDofUnknown(iUnknown).listNode().size()); i++) {
323 solution[i*3] = valueVec[i*3+numValueStart];
324 solution[i*3+1] = valueVec[i*3+1+numValueStart];
325 solution[i*3+2] = valueVec[i*3+2+numValueStart];
326 }
327 io[iMesh]->writeSolution(rank, time, iteration, 1, m_listVariable[idVar].name()+"_bkp", solution, supportDofUnknown(iUnknown).listNode().size(), dimension(),true);
328 break;
329
330 default:
331 FEL_ERROR("Problem with the dimension of your solution to write it.")
332 break;
333 }
334 }
335 }
336 // Do not need to delete here solution because it is deleted in EnsightCase::postProcess (manageMemory of io.writeSolution = true)
337 //delete [] solution;
338 }
339 }
340 }
341
342 /*!
343 \brief move a Petsc std::vector to a double* in order to print it in an ensight file
344 */
345 void LinearProblem::fromVecToDoubleStar(double* solution, const PetscVector& sol, int rank, int dimension, felInt size) {
346 if (FelisceParam::verbose() >= 40 )
347 sol.view();
348 felInt sizeVec;
349 sol.getSize(&sizeVec);
350
351 // Check if the demanded double* outpuy has the same dimension of the input Vec
352 if (size == 0)
353 size = sizeVec;
354
355 PetscVector xsol;
356 sol.scatterToZero(xsol, INSERT_VALUES, SCATTER_FORWARD);
357 if(rank == 0) {
358 felInt iout[size];
359 felReal valueVec[size];
360 for (felInt i=0; i<size; i++)
361 iout[i] = i;
362
363 AOApplicationToPetsc(m_ao, size, iout);
364 xsol.getValues(size, iout, valueVec);
365
366 switch (dimension) {
367 case 1:
368 for(felInt i=0; i < size; i++) {
369 solution[i] = valueVec[i];
370 }
371 break;
372 case 2:
373 for(felInt i=0; i < size; i++) {
374 solution[i*3] = valueVec[i*2];
375 solution[i*3+1] = valueVec[i*2+1];
376 solution[i*3+2] = 0.;
377 }
378 break;
379 case 3:
380 for(felInt i=0; i < size; i++) {
381 solution[i*3] = valueVec[i*3];
382 solution[i*3+1] = valueVec[i*3+1];
383 solution[i*3+2] = valueVec[i*3+2];
384 }
385 break;
386 default:
387 FEL_ERROR("Problem with the dimension of your solution to write it.")
388 break;
389 }
390 }
391 }
392
393 /*!
394 \brief move a double* to a Petsc std::vector in order to read it from an ensight file
395 */
396 void LinearProblem::fromDoubleStarToVec(double* solution, PetscVector* sol, felInt size) {
397 sol->duplicateFrom(m_sol);
398 felInt ia;
399 for (felInt i=0; i< size; i++) {
400 ia = i;
401 AOApplicationToPetsc(m_ao,1,&ia);
402 sol->setValues(1,&ia,&solution[i],INSERT_VALUES);
403 }
404
405 sol->assembly();
406
407 if (FelisceParam::verbose() >= 40 )
408 sol->view();
409 }
410
411 /*!
412 \brief write support variables for backup in ensight files
413 Warning:
414 dimension = 1 -> sclar; dimensio = 2 -> 2size 3D vectors; dimension = 3 -> 3size 3D vectors
415 */
416 void LinearProblem::writeSupportVariableBackup(std::string name, double* supportVariable, int size, int dimension, int rank, IO::Pointer io, double& time, int iteration, bool manageMemory) {
417 io->writeSolution(rank, time, iteration, dimension-1, name, supportVariable, size, 1, manageMemory);
418 }
419
420 /*!
421 \brief write mesh to read solution with specific mesh (P2, Q2,...) .
422 */
423 22 void LinearProblem::writeGeoForUnknown(int iUnknown, std::string nameOfGeoFile) {
424 // Intialization
425 // =============
426 22 std::string fileNameGeo, fileNameRoot, fileName;
427 ElementType eltType;
428 22 std::string keyWord;
429 int the_ref;
430 22 short verbose = 0;
431 FILE * pFile;
432 22 felInt numElemsPerRef = 0;
433 22 felInt ielSupportDof = 0;
434 22 std::vector<felInt> elem;
435 22 std::vector<felInt> vectorIdSupport;
436
437 22 const int idVar = m_listUnknown.idVariable(iUnknown);
438 22 const int iMesh = m_listVariable[idVar].idMesh();
439
440
7/14
✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 22 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 22 times.
✗ Branch 8 not taken.
✓ Branch 11 taken 22 times.
✗ Branch 12 not taken.
✓ Branch 15 taken 22 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 22 times.
✗ Branch 19 not taken.
✓ Branch 21 taken 22 times.
✗ Branch 22 not taken.
22 Ensight ens(FelisceParam::instance().inputDirectory,FelisceParam::instance().inputFile,FelisceParam::instance().inputMesh[iMesh],FelisceParam::instance().outputMesh[iMesh],FelisceParam::instance().meshDir,FelisceParam::instance().resultDir);
441
442 felInt numElement[ GeometricMeshRegion::m_numTypesOfElement ];
443
2/2
✓ Branch 0 taken 550 times.
✓ Branch 1 taken 22 times.
572 for (int ityp=0; ityp<GeometricMeshRegion::m_numTypesOfElement; ityp++ ) {
444 550 eltType = (ElementType)ityp;
445 550 numElement[eltType] = 0;
446 }
447
448
449
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 22 times.
22 if(nameOfGeoFile.empty()) {
450 fileNameGeo = listVariable().listVariable()[idVar].name()+".geo";
451 fileNameRoot = FelisceParam::instance().resultDir;
452 fileName = fileNameRoot + fileNameGeo;
453 } else {
454
1/2
✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
22 fileName = nameOfGeoFile;
455 }
456
457
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 22 times.
22 if (verbose>0)
458 std::cout << "Writing file "<< fileName << std::endl;
459
460
1/2
✓ Branch 2 taken 22 times.
✗ Branch 3 not taken.
22 pFile = fopen (fileName.c_str(),"w");
461
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 20 times.
22 if (pFile==nullptr) {
462
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
2 std::string command = "mkdir -p " + FelisceParam::instance().resultDir;
463
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
2 int ierr = system(command.c_str());
464
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
2 if (ierr>0) {
465 FEL_ERROR("Impossible to write "+fileName+" geo mesh");
466 }
467
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
2 pFile = fopen (fileName.c_str(),"w");
468 2 }
469 //! Descriptions lines (All Geometry files must contain these first six lines)
470
1/2
✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
22 fprintf(pFile,"Geometry file\n");
471
1/2
✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
22 fprintf(pFile,"Geometry file\n");
472
1/2
✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
22 fprintf(pFile,"node id assign\n");
473
1/2
✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
22 fprintf(pFile,"element id assign\n");
474
1/2
✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
22 fprintf(pFile,"coordinates\n");
475
1/2
✓ Branch 4 taken 22 times.
✗ Branch 5 not taken.
22 fprintf(pFile,"%8d\n",(int)supportDofUnknown(iUnknown).listNode().size());
476
477 //! write the vertices
478
1/2
✗ Branch 3 not taken.
✓ Branch 4 taken 22 times.
22 if ( m_mesh[iMesh]->numCoor() == 3 ) {
479 for(felInt iv=0; iv< static_cast<felInt>(supportDofUnknown(iUnknown).listNode().size()); iv++) {
480 fprintf( pFile,"%12.5e%12.5e%12.5e\n", supportDofUnknown(iUnknown).listNode()[iv].x(),
481 supportDofUnknown(iUnknown).listNode()[iv].y(),supportDofUnknown(iUnknown).listNode()[iv].z() );
482 }
483
1/2
✓ Branch 3 taken 22 times.
✗ Branch 4 not taken.
22 } else if ( m_mesh[iMesh]->numCoor() == 2 ) {
484 22 double zcoor = 0.;
485
2/2
✓ Branch 3 taken 40644 times.
✓ Branch 4 taken 22 times.
40666 for(felInt iv=0; iv< static_cast<felInt>(supportDofUnknown(iUnknown).listNode().size()); iv++) {
486
1/2
✓ Branch 9 taken 40644 times.
✗ Branch 10 not taken.
40644 fprintf( pFile, "%12.5e%12.5e%12.5e\n", supportDofUnknown(iUnknown).listNode()[iv].x(), supportDofUnknown(iUnknown).listNode()[iv].y(), zcoor);
487 }
488 } else {
489 std::ostringstream msg;
490 msg << "Error in writing geo File " << fileNameGeo << ". "
491 << "Expecting mesh with 2 or 3 coordinates (it is " << m_mesh[iMesh]->numCoor() << ")."
492 << std::endl;
493 FEL_ERROR(msg.str());
494 }
495
496 //!write the elements per references.
497 //------------
498 22 felInt cptpart = 0;
499 22 std::string DescriptionLine,felName;
500 22 GeometricMeshRegion::IntRefToElementType_type& intRefToEnum = m_mesh[iMesh]->intRefToEnum();
501
2/2
✓ Branch 4 taken 110 times.
✓ Branch 5 taken 22 times.
132 for(auto it_ref = intRefToEnum.begin(); it_ref != intRefToEnum.end(); ++it_ref) {
502 110 the_ref = (*it_ref).first;
503
504
2/2
✓ Branch 6 taken 110 times.
✓ Branch 7 taken 110 times.
220 for(auto it_elem = it_ref->second.begin(); it_elem != it_ref->second.end(); ++it_elem) {
505 110 eltType = *it_elem;
506
507
1/2
✓ Branch 3 taken 110 times.
✗ Branch 4 not taken.
110 if (m_mesh[iMesh]->flagFormatMesh() == GeometricMeshRegion::FormatMedit) {
508 110 cptpart++;
509
1/2
✓ Branch 1 taken 110 times.
✗ Branch 2 not taken.
110 fprintf ( pFile, "part%8d\n",cptpart );
510
1/2
✓ Branch 2 taken 110 times.
✗ Branch 3 not taken.
110 if ( m_listVariable[idVar].finiteElementType() == 0 )
511
1/2
✓ Branch 2 taken 110 times.
✗ Branch 3 not taken.
110 felName = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].first;
512 else if ( m_listVariable[idVar].finiteElementType() == 1 )
513 felName = GeometricMeshRegion::eltEnumToFelNameGeoEle[GeometricMeshRegion::eltLinearToEltQuad[eltType]].first;
514 else {
515 if(m_mesh[iMesh]->domainDim() == GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second->numCoor())
516 felName = GeometricMeshRegion::eltEnumToFelNameGeoEle[GeometricMeshRegion::eltLinearToEltBubble[eltType]].first;
517 else
518 felName = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].first;
519 }
520
521
2/4
✓ Branch 1 taken 110 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 110 times.
✗ Branch 5 not taken.
110 auto it_felName = FelisceParam::instance().felName2MapintRef2DescLine.find(felName);
522
2/4
✓ Branch 1 taken 110 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 110 times.
✗ Branch 6 not taken.
110 if ( it_felName == FelisceParam::instance().felName2MapintRef2DescLine.end() ) {
523 //std::cout << "Warning: ----------->trying to access an unexisting element: " << felName << " in input file !!!!!!!" << std::endl;
524
3/6
✓ Branch 1 taken 110 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 110 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 110 times.
✗ Branch 8 not taken.
110 DescriptionLine = felName + "m_ref_" + std::to_string(the_ref);
525 } else {
526 auto it_ref2 = it_felName->second.find(the_ref);
527 if ( it_ref2 == it_felName->second.end() ) {
528 DescriptionLine = felName + "m_ref_" + std::to_string(the_ref);
529 } else {
530 DescriptionLine = it_ref2->second;
531 }
532 }
533
1/2
✓ Branch 2 taken 110 times.
✗ Branch 3 not taken.
110 fprintf (pFile, "%s", DescriptionLine.c_str() );
534 }
535
1/2
✓ Branch 1 taken 110 times.
✗ Branch 2 not taken.
110 fprintf(pFile,"\n");
536
2/4
✓ Branch 2 taken 110 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 110 times.
✗ Branch 6 not taken.
110 keyWord = ens.eltFelNameToEns6Pair()[felName].second;
537
1/2
✓ Branch 2 taken 110 times.
✗ Branch 3 not taken.
110 fprintf( pFile, "%s", keyWord.c_str() );
538
1/2
✓ Branch 1 taken 110 times.
✗ Branch 2 not taken.
110 fprintf(pFile,"\n");
539
540
1/2
✓ Branch 3 taken 110 times.
✗ Branch 4 not taken.
110 numElemsPerRef = m_mesh[iMesh]->intRefToBegEndMaps[eltType][the_ref].second;
541
2/4
✓ Branch 1 taken 110 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 110 times.
110 if(!FelisceParam::instance().duplicateSupportDof) {
542 fprintf(pFile,"%8d\n", numElemsPerRef);
543 } else {
544 110 felInt numSupportElem = 0;
545
2/2
✓ Branch 0 taken 79596 times.
✓ Branch 1 taken 110 times.
79706 for ( felInt iel = numElement[eltType]; iel < numElement[eltType] + numElemsPerRef; iel++) {
546
1/2
✓ Branch 2 taken 79596 times.
✗ Branch 3 not taken.
79596 supportDofUnknown(iUnknown).getIdElementSupport(eltType, iel, vectorIdSupport);
547 79596 numSupportElem += vectorIdSupport.size();
548 }
549
1/2
✓ Branch 1 taken 110 times.
✗ Branch 2 not taken.
110 fprintf(pFile,"%8d\n", numSupportElem);
550 }
551
552
2/2
✓ Branch 0 taken 79596 times.
✓ Branch 1 taken 110 times.
79706 for ( felInt iel = 0; iel < numElemsPerRef; iel++) {
553
1/2
✓ Branch 2 taken 79596 times.
✗ Branch 3 not taken.
79596 supportDofUnknown(iUnknown).getIdElementSupport(eltType, numElement[eltType], vectorIdSupport);
554
555
2/2
✓ Branch 1 taken 80420 times.
✓ Branch 2 taken 79596 times.
160016 for(std::size_t iSup=0; iSup<vectorIdSupport.size(); iSup++) {
556 80420 ielSupportDof = vectorIdSupport[iSup];
557
1/2
✓ Branch 3 taken 80420 times.
✗ Branch 4 not taken.
80420 elem.resize(supportDofUnknown(iUnknown).getNumSupportDof(ielSupportDof),0);
558
2/2
✓ Branch 2 taken 236904 times.
✓ Branch 3 taken 80420 times.
317324 for ( int iSupport = 0; iSupport < supportDofUnknown(iUnknown).getNumSupportDof(ielSupportDof); iSupport++) {
559 236904 elem[iSupport] = supportDofUnknown(iUnknown).iSupportDof()[ supportDofUnknown(iUnknown).iEle()[ielSupportDof] + iSupport];
560 }
561
2/2
✓ Branch 1 taken 236904 times.
✓ Branch 2 taken 80420 times.
317324 for (unsigned int ivert = 0; ivert < elem.size(); ivert++) {
562
1/2
✓ Branch 2 taken 236904 times.
✗ Branch 3 not taken.
236904 fprintf( pFile, "%8d", elem[ivert]+1);
563 }
564
1/2
✓ Branch 1 taken 80420 times.
✗ Branch 2 not taken.
80420 fprintf(pFile,"\n");
565 }
566 79596 numElement[eltType]++;
567 }
568 }
569 }
570
1/2
✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
22 fclose(pFile);
571 22 }
572
573 void LinearProblem::writeGeoForUnknownEnsightGold(int iUnknown, std::map<felInt, std::vector<felInt> >& refToListOfIds, std::string nameOfGeoFile) {
574 const int idVar = m_listUnknown.idVariable(iUnknown);
575 const int iMesh = m_listVariable[idVar].idMesh();
576
577 if (m_mesh[iMesh]->flagFormatMesh() == GeometricMeshRegion::FormatUndefined) {
578 felisce_error("Empty Data Structure!!! Have you read mesh file before?",__FILE__,__LINE__);
579 }
580
581 // init std::unordered_map
582 std::unordered_map<std::string, std::pair<ElementType, std::string> > eltFelNameToEnsPair;
583 eltFelNameToEnsPair["Nodes"] = std::make_pair( GeometricMeshRegion::Nod, "point" );
584 eltFelNameToEnsPair["Segments2"] = std::make_pair( GeometricMeshRegion::Seg2, "bar2" );
585 eltFelNameToEnsPair["Segments3"] = std::make_pair( GeometricMeshRegion::Seg3, "bar3" );
586 eltFelNameToEnsPair["Triangles3"] = std::make_pair( GeometricMeshRegion::Tria3, "tria3" );
587 eltFelNameToEnsPair["Triangles6"] = std::make_pair( GeometricMeshRegion::Tria6, "tria6" );
588 eltFelNameToEnsPair["Quadrangles4"] = std::make_pair( GeometricMeshRegion::Quad4, "quad4" );
589 eltFelNameToEnsPair["Quadrangles8"] = std::make_pair( GeometricMeshRegion::Quad8, "quad8" );
590 eltFelNameToEnsPair["Tetrahedra4"] = std::make_pair( GeometricMeshRegion::Tetra4, "tetra4" );
591 eltFelNameToEnsPair["Tetrahedra10"] = std::make_pair( GeometricMeshRegion::Tetra10,"tetra10" );
592 eltFelNameToEnsPair["Pyramids5"] = std::make_pair( GeometricMeshRegion::Pyram5, "pyramid5" );
593 eltFelNameToEnsPair["Pyramids13"] = std::make_pair( GeometricMeshRegion::Pyram13,"pyramid13" );
594 eltFelNameToEnsPair["Prisms6"] = std::make_pair( GeometricMeshRegion::Prism6, "penta6" );
595 eltFelNameToEnsPair["Prisms15"] = std::make_pair( GeometricMeshRegion::Prism15,"penta15" );
596 eltFelNameToEnsPair["Hexahedra8"] = std::make_pair( GeometricMeshRegion::Hexa8, "hexa8" );
597 eltFelNameToEnsPair["Hexahedra20"] = std::make_pair( GeometricMeshRegion::Hexa20, "hexa20" );
598
599 // clear std::unordered_map
600 refToListOfIds.clear();
601
602 // Intialization
603 // =============
604 std::string fileNameGeo, fileNameRoot, fileName;
605 std::string keyWord;
606 FILE * pFile;
607
608 if(nameOfGeoFile.empty()) {
609 fileNameGeo = listVariable().listVariable()[idVar].name()+".geo";;
610 fileNameRoot = FelisceParam::instance().resultDir;
611 fileName = fileNameRoot + fileNameGeo;
612 } else {
613 fileName = nameOfGeoFile;
614 }
615
616 pFile = fopen (fileName.c_str(),"w");
617 if (pFile==nullptr) {
618 std::string command = "mkdir -p " + FelisceParam::instance().resultDir;
619 int ierr = system(command.c_str());
620 if (ierr>0) {
621 FEL_ERROR("Impossible to write "+fileName+" geo mesh");
622 }
623 pFile = fopen (fileName.c_str(),"w");
624 }
625
626 // Descriptions lines (All Geometry files must contain these first six lines)
627 fprintf(pFile,"Geometry file\n");
628 fprintf(pFile,"Ensight Gold format\n");
629 fprintf(pFile,"node id given\n");
630 fprintf(pFile,"element id assign\n");
631
632 // list of type of element by labels
633 GeometricMeshRegion::IntRefToElementType_type& intRefToEnum = m_mesh[iMesh]->intRefToEnum();
634
635 felInt label;
636 bool(*comparison)(std::pair<felInt, Point >, std::pair<felInt, Point >) = Tools::pairComparison;
637
638 felInt numElemsPerRef = 0;
639 std::vector<felInt> vectorIdSupport;
640 felInt ielSupportDof;
641 felInt idSup;
642 Point pt;
643 std::vector<felInt> numElement(GeometricMeshRegion::m_numTypesOfElement, 0);
644 std::vector<felInt> numElement2(GeometricMeshRegion::m_numTypesOfElement, 0);
645 std::vector<felInt> numSupportElt(GeometricMeshRegion::m_numTypesOfElement, 0);
646
647 // loop over references
648 for (auto refIt = intRefToEnum.begin(); refIt != intRefToEnum.end(); ++refIt) {
649 label = refIt->first;
650 fprintf(pFile, "part\n");
651 fprintf(pFile, "%10d\n", label + 1);
652
653 std::string descriptionLine = "";
654 for (auto typeIt = refIt->second.begin(); typeIt != refIt->second.end(); ++typeIt)
655 descriptionLine += GeometricMeshRegion::eltEnumToFelNameGeoEle[*typeIt].first + "_";
656
657 descriptionLine += "ref_";
658 fprintf(pFile, "%s%d\n", descriptionLine.c_str(), label);
659 fprintf(pFile, "coordinates\n");
660
661 std::map<felInt, felInt> local2GlobalMapping;
662 std::set<std::pair<felInt, Point >, bool(*)(std::pair<felInt, Point >, std::pair<felInt, Point >)> listOfNodes(comparison);
663
664 // loop over the type of element having this reference
665 for (auto typeIt = refIt->second.begin(); typeIt != refIt->second.end(); ++typeIt) {
666 // fill the mappings
667 numElemsPerRef = m_mesh[iMesh]->intRefToBegEndMaps[*typeIt][label].second;
668 numSupportElt[*typeIt] = 0;
669
670 // loop over the element
671 for(felInt iel=0; iel<numElemsPerRef; ++iel) {
672 supportDofUnknown(iUnknown).getIdElementSupport(*typeIt, numElement[*typeIt], vectorIdSupport);
673
674 for(std::size_t iSup=0; iSup<vectorIdSupport.size(); ++iSup) {
675 ielSupportDof = vectorIdSupport[iSup];
676 for(int iSupport=0; iSupport<supportDofUnknown(iUnknown).getNumSupportDof(ielSupportDof); iSupport++) {
677 idSup = supportDofUnknown(iUnknown).iSupportDof()[supportDofUnknown(iUnknown).iEle()[ielSupportDof] + iSupport];
678 pt = supportDofUnknown(iUnknown).listNode()[idSup];
679 listOfNodes.insert(std::make_pair(idSup, pt));
680 }
681 ++numSupportElt[*typeIt];
682 }
683
684 ++numElement[*typeIt];
685 }
686 }
687
688 // print part and vertices
689 fprintf(pFile, "%10d\n", static_cast<int>(listOfNodes.size()));
690 felInt localCounter = 1;
691 for (auto nodeIt = listOfNodes.begin(); nodeIt != listOfNodes.end(); ++nodeIt) {
692 refToListOfIds[label].push_back(nodeIt->first);
693
694 fprintf(pFile, "%10d\n", nodeIt->first + 1);
695 local2GlobalMapping[nodeIt->first] = localCounter;
696 ++localCounter;
697 }
698
699 for(felInt icoor=0; icoor<3; ++icoor) {
700 for (auto nodeIt = listOfNodes.begin(); nodeIt != listOfNodes.end(); ++nodeIt) {
701 fprintf(pFile, "%12.5e\n", nodeIt->second[icoor]);
702 }
703 }
704
705 for (auto typeIt = refIt->second.begin(); typeIt != refIt->second.end(); ++typeIt) {
706 keyWord = eltFelNameToEnsPair[GeometricMeshRegion::eltEnumToFelNameGeoEle[*typeIt].first].second;
707 fprintf(pFile, "%s\n", keyWord.c_str());
708 fprintf(pFile, "%10d\n", numSupportElt[*typeIt]);
709
710 // loop over the element
711 numElemsPerRef = m_mesh[iMesh]->intRefToBegEndMaps[*typeIt][label].second;
712 for(felInt iel=0; iel<numElemsPerRef; ++iel) {
713 supportDofUnknown(iUnknown).getIdElementSupport(*typeIt, numElement2[*typeIt], vectorIdSupport);
714
715 for(std::size_t iSup=0; iSup<vectorIdSupport.size(); ++iSup) {
716 ielSupportDof = vectorIdSupport[iSup];
717 for(int iSupport=0; iSupport<supportDofUnknown(iUnknown).getNumSupportDof(ielSupportDof); iSupport++) {
718 fprintf(pFile, "%10d", local2GlobalMapping[supportDofUnknown(iUnknown).iSupportDof()[supportDofUnknown(iUnknown).iEle()[ielSupportDof] + iSupport]]);
719 }
720 fprintf(pFile, "\n");
721 }
722
723 ++numElement2[*typeIt];
724 }
725 }
726 }
727
728 fclose(pFile);
729 }
730 }//
731