GCC Code Coverage Report

Directory: ./
File: Solver/eigenProblem.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 485 870 55.7%
Branches: 379 1142 33.2%
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 "Solver/eigenProblem.hpp"
21 #include "FiniteElement/elementMatrix.hpp"
22 #include "FiniteElement/elementField.hpp"
23 #include "Core/felisceParam.hpp"
24
25 namespace felisce
26 {
27 /*! Constructor
28 */
29
6/12
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 15 taken 1 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 1 times.
✗ Branch 19 not taken.
✓ Branch 21 taken 1 times.
✗ Branch 22 not taken.
✓ Branch 36 taken 1 times.
✗ Branch 37 not taken.
 1 EigenProblem::EigenProblem()
30 {
31
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 nullPetscVector = PetscVector::null();
32 1 }
33
34 /***********************************************************************************/
35 /***********************************************************************************/
36
37 2 EigenProblem::~EigenProblem()
38 {
39
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 2 if(m_seqSolAllocated) {
40 2 m_seqSol.zeroEntries();
41 2 m_seqSol.destroy();
42 }
43
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 2 if (m_initNumDofUnknown) {
44
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 2 delete [] m_numDofUnknown;
45
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 2 delete [] m_numDofLocalUnknown;
46 }
47
48
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 2 if (m_buildSystem ) {
49 2 m_sol.destroy();
50 }
51
52 #ifdef FELISCE_WITH_SLEPC
53 2 EPSDestroy(&m_eps);
54 #endif
55
56 }
57
58 /***********************************************************************************/
59 /***********************************************************************************/
60
61 1 void EigenProblem::initialize(const GeometricMeshRegion::Pointer& mesh, FelisceTransient::Pointer fstransient, MPI_Comm& comm)
62 {
63 1 m_fstransient = fstransient;
64 1 m_petscComm = comm;
65 1 m_verbose = FelisceParam::verbose();
66 1 m_mesh = mesh;
67
68 #ifdef FELISCE_WITH_SLEPC
69 1 EPSCreate(comm, &m_eps);
70 #endif
71 1 }
72
73 /***********************************************************************************/
74 /***********************************************************************************/
75
76 1 void EigenProblem::definePhysicalVariable(const std::vector<PhysicalVariable>& listVariable, const std::vector<std::size_t>& listNumComp)
77 {
78
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
 2 for (std::size_t iVar = 0; iVar < listVariable.size(); iVar++) {
79
1/2
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
 1 m_listVariable.addVariable(listVariable[iVar],listNumComp[iVar]);
80 }
81
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
 1 if (listVariable.size() != m_listVariable.size()) {
82 ✗ FEL_ERROR("Problem with your input file. Not Enought physical variable defined for this problem!");
83 }
84
85 //link between unknown of the linear system and variable define in the problem.
86
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 Variable variable;
87 1 bool unknown = false;
88
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
 2 for (std::size_t iUnknown = 0; iUnknown < m_listUnknown.size(); iUnknown++) {
89
5/6
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
✓ Branch 6 taken 1 times.
 2 for (std::size_t iVar = 0; iVar < listVariable.size() && unknown == false; iVar++) {
90
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 variable = m_listVariable[iVar];
91
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 if ( m_listUnknown[iUnknown] == variable.physicalVariable()) {
92 1 m_listUnknown.idVariable(iUnknown) = iVar;
93 1 m_listVariable.listIdUnknownOfVariable(iVar) = static_cast<int>(iUnknown);
94 1 unknown = true;
95 }
96 }
97 1 unknown = false;
98 }
99
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if (m_verbose) {
100
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 PetscPrintf(m_petscComm, "\n/========== Information about the linear system ==========/\n");
101
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 PetscPrintf(m_petscComm, "<%s>.\n", nameOfTheProblem.c_str());
102
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 PetscPrintf(m_petscComm, "Unknonwns are: ");
103
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
 2 for (std::size_t iUnkonwn = 0; iUnkonwn < m_listUnknown.size(); iUnkonwn++) {
104
2/4
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
 1 PetscPrintf(m_petscComm, "<%s> ", m_listVariable[m_listUnknown.idVariable(iUnkonwn)].name().c_str());
105 }
106
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 PetscPrintf(m_petscComm, "\nPhysical variables associate: ");
107 }
108
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_listVariable.print(m_verbose);
109
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_listUnknown.setDefaultMask(m_listVariable);
110
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_listUnknown.print(1);
111
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if (m_numDofUnknown == nullptr)
112
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
 1 m_numDofUnknown = new felInt[m_listUnknown.size()];
113
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if (m_numDofLocalUnknown == nullptr)
114
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
 1 m_numDofLocalUnknown = new felInt[m_listUnknown.size()];
115 1 m_initNumDofUnknown = true;
116
117 //todo Jeremie 01 02 2012: assemble multiple matrix
118
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_Matrix.resize(m_numberOfMatrix);
119 1 }
120
121 /***********************************************************************************/
122 /***********************************************************************************/
123
124 1 void EigenProblem::computeDof(int numProc, int idProc)
125 {
126 // felInt numElement = m_mesh->getNumDomainElement() + m_mesh->getNumBoundaryElement();
127
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
 2 for ( std::size_t iUnknown = 0; iUnknown < m_listUnknown.size(); iUnknown++ ) {
128
1/2
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
 1 SupportDofMesh supportDofMesh(m_listVariable[m_listUnknown.idVariable(iUnknown)], m_mesh);
129
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_supportDofUnknown.push_back(supportDofMesh);
130
131
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if (m_verbose) {
132
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 PetscPrintf(m_petscComm, "\n/=======Global support dof mesh associate to the unknown: <%s>.======/\n",
133
1/2
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
 2 m_listVariable[m_listUnknown.idVariable(iUnknown)].name().c_str());
134
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 supportDofMesh.print(m_verbose);
135 }
136 1 }
137
138
1/2
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
 1 std::vector<SupportDofMesh*> pSupportDofUnknown(m_supportDofUnknown.size(), nullptr);
139
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
 2 for(std::size_t iUnknown=0; iUnknown < m_supportDofUnknown.size(); ++iUnknown)
140 1 pSupportDofUnknown[iUnknown] = &m_supportDofUnknown[iUnknown];
141
142
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_dof.setDof(m_listUnknown, m_listVariable, pSupportDofUnknown);
143 //m_dof.buildPattern();
144 //Initialize pattern, random repartition on processors before call ParMetis.
145
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_dof.initializePattern(numProc, idProc);
146
147 1 m_numDof = m_dof.numDof();
148
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
 2 for ( std::size_t iUnknown = 0; iUnknown < m_listUnknown.size(); iUnknown++ ) {
149 1 m_numDofUnknown[iUnknown] = m_dof.numDofPerUnknown()[iUnknown];
150 }
151
152
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if (m_verbose) {
153
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 PetscPrintf(m_petscComm, "\n/========== Information about dof ==========/\n");
154 }
155
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_dof.print(m_verbose);
156 1 }
157
158 /***********************************************************************************/
159 /***********************************************************************************/
160
161 1 void EigenProblem::partitionDof( idx_t numPart, int rankPart, MPI_Comm comm)
162 {
163 /*!
164 Partitioning of the dof.
165 */
166
167 /*
168 Declaration of the parameters for ParMetis.
169 idx_t define the type for ParMetis int.
170 */
171
172 // Input parameters
173 1 idx_t num_flag = 0; // C-style array index.
174 1 idx_t ncon = 1; // number of weight per vertices of the graph.
175 1 idx_t wgtFlag = 0; // No weight in the graph.
176 idx_t option[3]; // verbose level and random seed for parmetis.
177 1 idx_t vertexdist[numPart + 1]; // Repartition of dof between the processes.
178
179 1 real_t ubvec = static_cast<real_t>(1.05); // Imbalance tolerance
180
181 // Output parameters
182 idx_t edgecut; // Number of edge in the graph that are cut by the partitioning.
183 1 int numRows = m_dof.pattern().numRows();
184
1/2
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 std::vector<idx_t> dofPartitionning(m_dof.pattern().numRows(), 0); // The partition of the dof.
185
186 /*
187 Initialisation of the parameters.
188 */
189
190 // vertexdist
191 // vertexdist[j] = sum(size in process i), for i=0,..,j-1.
192 // initialize the first value (always 0)
193 1 vertexdist[0] = 0;
194
195 // Gather all the size in vertexdist, starting at index 1 (because vertexdist[0] is 0)
196
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 MPI_Allgather(&numRows, 1, MPI_INT, &vertexdist[1], 1, MPI_INT, comm);
197
198 // Compute the partial sum
199
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 for (int i = 1; i < numPart; i ++)
200 ✗ vertexdist[i + 1] += vertexdist[i];
201
202 // option
203 1 option[0] = 1; // 0 = default values, 1 = user defined.
204 1 option[1] = 100; // verbose level.
205 1 option[2] = 100; // random seed.
206
207 {
208 // tpwgts
209 1 real_t tpwgts[numPart*ncon]; // See parmetis manual for description
210
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 for ( felInt i= 0; i < numPart*ncon; i++)
211 1 tpwgts[i] = static_cast<real_t>(1.)/numPart;
212
213 /*
214 Parmetis call
215 */
216
217 1 MPI_Comm comm_bis = PETSC_COMM_WORLD;
218
3/6
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
 1 ParMETIS_V3_PartKway(vertexdist, &m_dof.pattern().rowPointer(0), &m_dof.pattern().columnIndex(0), nullptr, nullptr, &wgtFlag, &num_flag, &ncon, &numPart, tpwgts, &ubvec, option, &edgecut, dofPartitionning.data(), &comm_bis);
219 1 }
220
221 /*
222 Gather the partitioning to all process
223 */
224
225
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 m_dofPart.resize(m_dof.numDof(), 0);
226
227 {
228 1 felInt recvcount[numPart];
229
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 for(int i=0; i<numPart; i++)
230 1 recvcount[i] = vertexdist[i+1] - vertexdist[i];
231
232
1/2
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 MPI_Allgatherv(dofPartitionning.data(), dofPartitionning.size(), MPI_INT, m_dofPart.data(), recvcount, vertexdist, MPI_INT, comm);
233 1 }
234
235 // Count the number of local dofs
236
2/2
✓ Branch 1 taken 4225 times.
✓ Branch 2 taken 1 times.
 4226 for (felInt i = 0; i < m_dof.numDof(); i++) {
237
1/2
✓ Branch 1 taken 4225 times.
✗ Branch 2 not taken.
 4225 if (m_dofPart[i] == rankPart)
238 4225 m_numDofLocal++;
239 }
240
241 /*
242 Build the matrix pattern
243 */
244
245
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 m_dof.buildPattern(rankPart, m_dofPart.data());
246
247 // Pattern operations
248
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
 1 if(m_dof.patternIsChanged()) {
249 ✗ m_dof.mergeLocalCompPattern(m_petscComm, numPart, rankPart, m_dofPart);
250
251 //free memory
252 ✗ m_dof.clearPatternC();
253 }
254
255 /*!
256 Global ordering: AO is an application between global mesh ordering and global ordering of Petsc
257 AO: globalordering to globalordering
258 the action is:
259 - Take the global index associated to the i-th proc
260 - Create an application to a new global index in order to have contiguous global indexes on the processor
261
262 both global ordering are on input (loc2glob and pordering)
263 AOCreateBasic(comm, m_numDofLocal, loc2Glob, pordering, &m_ao);
264
265 Example:
266 Proc 0:
267 glob2loc = [0 1 2 3 8]
268 pordering =[0 1 2 3 4]
269
270 Proc 1:
271 glob2loc = [4 5 6 7]
272 pordering =[5 6 7 8]
273 */
274 {
275 felInt rstart;
276 1 felInt cptDof = 0;
277
278 1 felInt loc2Glob[m_numDofLocal];
279
2/2
✓ Branch 0 taken 4225 times.
✓ Branch 1 taken 1 times.
 4226 for (felInt i = 0; i < m_numDofLocal; i++)
280 4225 loc2Glob[i] = 0;
281
282
2/2
✓ Branch 1 taken 4225 times.
✓ Branch 2 taken 1 times.
 4226 for (felInt i = 0; i < m_dof.numDof(); i++) {
283
1/2
✓ Branch 1 taken 4225 times.
✗ Branch 2 not taken.
 4225 if (m_dofPart[i] == rankPart) {
284 // cptDof (local) -> i (global)
285 4225 loc2Glob[cptDof] = i;
286 4225 cptDof++;
287 }
288 }
289
290 // inclusive sum (rank_0 = numdoflocal_0, rank_i = sum_j=0^rank_i numdoflocal_i)
291
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 MPI_Scan(&m_numDofLocal, &rstart, 1, MPIU_INT, MPI_SUM, comm);
292 1 rstart -= m_numDofLocal;
293
294 1 felInt pordering[m_numDofLocal];
295
2/2
✓ Branch 0 taken 4225 times.
✓ Branch 1 taken 1 times.
 4226 for (felInt i= 0; i < m_numDofLocal; i++) {
296 // pordering: so global indexes are contiguous on a processor; pordering: global->global
297 4225 pordering[i] = rstart + i;
298 }
299
300
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 AOCreateBasic(comm, m_numDofLocal, loc2Glob, pordering, &m_ao);
301 // AOView(m_ao,PETSC_VIEWER_STDOUT_WORLD);
302 1 }
303
304 /*!
305 Compute the element repartition.
306 An element is owned by a process k if most of its dof are on process k.
307 In case of a draw, the element belongs to the process with the smallest id.
308 */
309
310 // Declaration of some variables
311 1 felInt max = 0; // used to compute indice_max.
312 1 felInt indice_max = 0; // id of the process that owns the element.
313 felInt idDof; // Local number of a dof.
314 felInt ielSupportDof; // Global number of a dof.
315 1 felInt numElemsPerRef = 0; // Number of elements per label.
316 int idVar; // id of a physical variable.
317 1 felInt numEltTot = 0; // Total number of element (to be computed).
318 felInt numEltPerType[GeometricMeshRegion::m_numTypesOfElement]; // Total number of element per type (to be computed).
319 1 felInt numDofPerProcs[numPart]; // Number of dofs per process.
320 GeometricMeshRegion::ElementType eltType;
321
322 // Initialization
323
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 for ( int i = 0; i < numPart; i++)
324 1 numDofPerProcs[i] = 0;
325
326
2/2
✓ Branch 0 taken 25 times.
✓ Branch 1 taken 1 times.
 26 for (int ityp=0; ityp<GeometricMeshRegion::m_numTypesOfElement; ityp++)
327 25 numEltPerType[ityp] = 0;
328
329
2/4
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
 1 m_eltPart.resize(m_mesh->getNumElement(), 0);
330
331 // bag containing all the domain element type.
332 1 const std::vector<ElementType>& bagElementTypeDomain =m_mesh->bagElementTypeDomain();
333
334 // for each domain element type
335
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
 2 for (std::size_t i = 0; i < bagElementTypeDomain.size(); ++i) {
336 1 eltType = bagElementTypeDomain[i];
337
338 // for each label where the current element type can be found
339
2/2
✓ Branch 6 taken 1 times.
✓ Branch 7 taken 1 times.
 2 for(auto itRef =m_mesh->intRefToBegEndMaps[eltType].begin(); itRef !=m_mesh->intRefToBegEndMaps[eltType].end(); itRef++) {
340 1 numElemsPerRef = itRef->second.second;
341
342 // for each local element with the current label
343
2/2
✓ Branch 0 taken 8192 times.
✓ Branch 1 taken 1 times.
 8193 for (felInt iel = 0; iel < numElemsPerRef; iel++) {
344
345 // for each unknown of the linear system
346
2/2
✓ Branch 1 taken 8192 times.
✓ Branch 2 taken 8192 times.
 16384 for (std::size_t iUnknown = 0; iUnknown < m_listUnknown.size(); iUnknown++) {
347 // get the global id of the current element
348 8192 idVar = m_listUnknown.idVariable(iUnknown);
349
1/2
✓ Branch 2 taken 8192 times.
✗ Branch 3 not taken.
 8192 m_supportDofUnknown[iUnknown].getIdElementSupport(eltType, numEltPerType[eltType], ielSupportDof);
350
351 // for each support dof of the current element
352
2/2
✓ Branch 2 taken 24576 times.
✓ Branch 3 taken 8192 times.
 32768 for (felInt iSupport = 0; iSupport < m_supportDofUnknown[iUnknown].getNumSupportDof(ielSupportDof); iSupport++) {
353
354 // for each components of the current unknown
355
2/2
✓ Branch 2 taken 24576 times.
✓ Branch 3 taken 24576 times.
 49152 for (std::size_t iComp = 0; iComp < m_listVariable[idVar].numComponent(); iComp++) {
356 // get the global id of the dof and increased numDofPerProcs
357
1/2
✓ Branch 1 taken 24576 times.
✗ Branch 2 not taken.
 24576 m_dof.loc2glob(ielSupportDof, iSupport, idVar, iComp, idDof);
358 24576 ++numDofPerProcs[m_dofPart[idDof]];
359 }
360 }
361 }
362
363 // Find the rank of the process that owns the most dof in the current element.
364
2/2
✓ Branch 0 taken 8192 times.
✓ Branch 1 taken 8192 times.
 16384 for (int j = 0; j < numPart; j++) {
365
1/2
✓ Branch 0 taken 8192 times.
✗ Branch 1 not taken.
 8192 if (max < numDofPerProcs[j]) {
366 8192 max = numDofPerProcs[j];
367 8192 indice_max = j;
368 }
369 }
370
371 // Fill m_eltPart
372 8192 m_eltPart[numEltTot] = indice_max;
373
374 // Increment the counter
375 8192 numEltPerType[eltType]++;
376 8192 numEltTot++;
377
378 // Clear local variables
379 8192 max = 0;
380 8192 indice_max = 0;
381
2/2
✓ Branch 0 taken 8192 times.
✓ Branch 1 taken 8192 times.
 16384 for ( int j = 0; j < numPart; j++)
382 8192 numDofPerProcs[j] = 0;
383 }
384 }
385 }
386
387 // bag containing all the "boundary" element type.
388 1 const std::vector<ElementType>& bagElementTypeDomainBoundary =m_mesh->bagElementTypeDomainBoundary();
389
390 // Do the same thing but with the "boundaries" element this time
391
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
 2 for (std::size_t i = 0; i < bagElementTypeDomainBoundary.size(); ++i) {
392 1 eltType = bagElementTypeDomainBoundary[i];
393
394
2/2
✓ Branch 6 taken 4 times.
✓ Branch 7 taken 1 times.
 5 for(auto itRef =m_mesh->intRefToBegEndMaps[eltType].begin(); itRef !=m_mesh->intRefToBegEndMaps[eltType].end(); itRef++) {
395 4 numElemsPerRef = itRef->second.second;
396
397
2/2
✓ Branch 0 taken 256 times.
✓ Branch 1 taken 4 times.
 260 for (felInt iel = 0; iel < numElemsPerRef; iel++) {
398
399
2/2
✓ Branch 1 taken 256 times.
✓ Branch 2 taken 256 times.
 512 for (std::size_t iUnknown = 0; iUnknown < m_listUnknown.size(); iUnknown++) {
400 256 idVar = m_listUnknown.idVariable(iUnknown);
401
1/2
✓ Branch 2 taken 256 times.
✗ Branch 3 not taken.
 256 m_supportDofUnknown[iUnknown].getIdElementSupport(eltType, numEltPerType[eltType], ielSupportDof);
402
403
2/2
✓ Branch 2 taken 512 times.
✓ Branch 3 taken 256 times.
 768 for (int iSupport = 0; iSupport < m_supportDofUnknown[iUnknown].getNumSupportDof(ielSupportDof); iSupport++) {
404
405
2/2
✓ Branch 2 taken 512 times.
✓ Branch 3 taken 512 times.
 1024 for (std::size_t iComp = 0; iComp < m_listVariable[idVar].numComponent(); iComp++) {
406
1/2
✓ Branch 1 taken 512 times.
✗ Branch 2 not taken.
 512 m_dof.loc2glob(ielSupportDof,iSupport,idVar,iComp,idDof);
407 512 numDofPerProcs[ m_dofPart[idDof] ]++;
408 }
409 }
410 }
411
412
2/2
✓ Branch 0 taken 256 times.
✓ Branch 1 taken 256 times.
 512 for (int j = 0; j < numPart; j++) {
413
1/2
✓ Branch 0 taken 256 times.
✗ Branch 1 not taken.
 256 if (max < numDofPerProcs[j] ) {
414 256 max = numDofPerProcs[j];
415 256 indice_max = j;
416 }
417 }
418
419 256 m_eltPart[numEltTot] = indice_max;
420 256 numEltPerType[eltType]++;
421 256 numEltTot++;
422
423 256 max = 0;
424 256 indice_max = 0;
425
2/2
✓ Branch 0 taken 256 times.
✓ Branch 1 taken 256 times.
 512 for (int j = 0; j < numPart; j++)
426 256 numDofPerProcs[j] = 0;
427 }
428 }
429 }
430 2 }
431
432 /***********************************************************************************/
433 /***********************************************************************************/
434
435 ✗ void EigenProblem::getSupportCoordinate(felInt idElt,
436 felInt idSupport,
437 int iUnknown,
438 Point& pt)
439 {
440 ✗ const felInt id = m_supportDofUnknownLocal[iUnknown].iSupportDof()[m_supportDofUnknownLocal[iUnknown].iEle()[idElt]
441 ✗ + idSupport];
442 ✗ pt = m_supportDofUnknownLocal[iUnknown].listNode()[id];
443 }
444
445 /***********************************************************************************/
446 /***********************************************************************************/
447
448 1 void EigenProblem::cutMesh()
449 {
450 1 partitionDof(static_cast<idx_t>(MpiInfo::numProc()), MpiInfo::rankProc(), MpiInfo::petscComm());
451 1 setLocalMeshAndSupportDofMesh(MpiInfo::rankProc());
452 1 }
453
454 /***********************************************************************************/
455 /***********************************************************************************/
456
457 1 void EigenProblem::setLocalMeshAndSupportDofMesh(int rankPart) {
458 1 int idVar = -1;
459
460 // Set the local mesh and fill the local to global mapping for the support elements (loc2GlobSupportElem)
461 1 std::vector<felInt> loc2GlobElem;
462
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_meshLocal = felisce::make_shared<GeometricMeshRegion>();
463
1/2
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
 1 m_meshLocal->setLocalMesh(*m_mesh, m_eltPart, rankPart, loc2GlobElem);
464
465
466 // for each unknown, create the support dof and put them in m_supportDofUnknownLocal
467
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
 2 for ( std::size_t iUnknown = 0; iUnknown < m_listUnknown.size(); iUnknown++ ) {
468 1 idVar = m_listUnknown.idVariable(iUnknown);
469
1/2
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
 1 SupportDofMesh supportDofMesh(m_listVariable[idVar], m_meshLocal, loc2GlobElem, m_supportDofUnknown[iUnknown]);
470
471
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_supportDofUnknownLocal.push_back(supportDofMesh);
472 1 }
473
474 // Create the local to global mapping for the elements
475
1/2
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
 1 IS_LOCAL_TO_GLOBAL_MAPPING_MACRO(m_petscComm, 1, loc2GlobElem.size(), loc2GlobElem.data(), PETSC_COPY_VALUES, &m_mappingElem);
476
477 // print info on the local mesh and supportDofMeshes
478
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if (m_verbose) {
479 // print the local mesh
480
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 PetscPrintf(m_petscComm, "\n/================== Local mesh ================/\n");
481
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 m_meshLocal->print(m_verbose);
482
483 // print the local supportDofMesh for each unknown
484
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
 2 for(std::size_t iUnknown=0; iUnknown<m_listUnknown.size(); iUnknown++) {
485
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 PetscPrintf(m_petscComm,
486 "\n/=========== Local support dof mesh associate to the unknown: <%s> ===========/\n",
487
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 2 m_listVariable[idVar].name().c_str());
488
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 m_supportDofUnknownLocal[iUnknown].print(m_verbose);
489 }
490 }
491 1 }
492
493 /***********************************************************************************/
494 /***********************************************************************************/
495
496 1 void EigenProblem::allocateMatrix()
497 {
498 // Create the local to global mapping for support dofs.
499 1 felInt loc2Glob[m_numDofLocal];
500 1 felInt cptDof= 0;
501
2/2
✓ Branch 0 taken 4225 times.
✓ Branch 1 taken 1 times.
 4226 for ( felInt iLDof = 0; iLDof < m_numDofLocal; iLDof++)
502 4225 loc2Glob[iLDof] = 0;
503
504
2/2
✓ Branch 1 taken 4225 times.
✓ Branch 2 taken 1 times.
 4226 for ( felInt iDof = 0; iDof < m_dof.numDof(); iDof++) {
505
2/4
✓ Branch 2 taken 4225 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 4225 times.
✗ Branch 5 not taken.
 4225 if ( m_dofPart[iDof] == MpiInfo::rankProc() ) {
506 4225 loc2Glob[cptDof] = iDof;
507 4225 cptDof++;
508 }
509 }
510
511 // global to local mapping on nodes
512
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 IS_LOCAL_TO_GLOBAL_MAPPING_MACRO(m_petscComm, 1, m_numDofLocal, loc2Glob, PETSC_COPY_VALUES, &m_mappingNodes);
513 // ISLocalToGlobalMappingView( m_mappingNodes, PETSC_VIEWER_STDOUT_WORLD);
514
515
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
 1 if ( MpiInfo::numProc() > 1 ) {
516 ✗ felInt iCSR[m_dof.pattern().numRows()+1];
517 ✗ felInt jCSR[m_dof.pattern().numNonzeros()];
518
519 ✗ iCSR[0] = 0;
520 ✗ std::set<felInt> jCSRSorted;
521 felInt connect;
522 felInt globValue;
523 ✗ for ( felInt ii = 0; ii < m_numDofLocal; ii++) {
524 ✗ connect = m_dof.pattern().numNonzerosInRow(ii);
525 ✗ iCSR[ii+1] = iCSR[ii] + connect;
526 ✗ for ( felInt jj = 0; jj < connect; jj++) {
527 ✗ globValue = m_dof.pattern().columnIndex( m_dof.pattern().rowPointer(ii) + jj );
528 ✗ AOApplicationToPetsc(m_ao, 1, &globValue);
529 ✗ jCSRSorted.insert( globValue );
530 }
531 ✗ felInt kk = 0;
532 ✗ for(auto it_jCSR = jCSRSorted.begin(); it_jCSR != jCSRSorted.end(); it_jCSR++) {
533 ✗ jCSR[ iCSR[ii] + kk] = *it_jCSR;
534 ✗ kk++;
535 }
536 ✗ jCSRSorted.clear();
537 }
538 ✗ for (std::size_t ii = 0; ii < m_Matrix.size(); ii++) {
539 ✗ m_Matrix[ii].createAIJ(m_petscComm, m_numDofLocal, m_numDofLocal, m_dof.numDof(), m_dof.numDof(), 0, FELISCE_PETSC_NULLPTR, 0, FELISCE_PETSC_NULLPTR);
540 ✗ m_Matrix[ii].setFromOptions();
541 ✗ m_Matrix[ii].mpiAIJSetPreallocationCSR( iCSR, jCSR, nullptr);
542 ✗ m_Matrix[ii].setOption(MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
543 }
544 ✗ } else {
545 1 const std::size_t numRows = m_dof.pattern().numRows();
546 1 felInt nnz[ numRows ];
547
2/2
✓ Branch 0 taken 4225 times.
✓ Branch 1 taken 1 times.
 4226 for ( felInt idof = 0; idof < static_cast<felInt>(numRows) ; idof++) {
548
1/2
✓ Branch 2 taken 4225 times.
✗ Branch 3 not taken.
 4225 nnz[idof] = m_dof.pattern().numNonzerosInRow(idof);
549 }
550
551
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 1 times.
 3 for (std::size_t ii = 0; ii < m_Matrix.size(); ii++) {
552
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 m_Matrix[ii].createSeqAIJ(m_petscComm, numRows, numRows, 0, nnz);
553
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 m_Matrix[ii].setFromOptions();
554
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 m_Matrix[ii].setOption(MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
555 }
556 1 }
557
558 1 m_numDof = m_dof.numDof();
559
560 // Create the std::vector solution
561
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_sol.createMPI(m_petscComm, m_numDofLocal, m_numDof);
562
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_sol.set( 0.);
563
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 allocateSequentialSolution();
564 1 m_buildSystem = true;
565
566 // Get the local size of all process
567
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 m_localSizeOfVector.resize(MpiInfo::numProc());
568
3/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 1 times.
 2 for (int ii = 0; ii < MpiInfo::numProc(); ii++)
569 1 m_localSizeOfVector[ii] = 0;
570
571 1 felInt sizeLocal = 0;
572
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_sol.getLocalSize(&sizeLocal);
573
574
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 MPI_Allgather(&sizeLocal, 1, MPI_INT, &m_localSizeOfVector[0], 1, MPI_INT, m_petscComm);
575
576 // For each unknown, create a local to global mapping to std::unordered_map the local id of dof to their global id in the petsc ordering.
577 1 felInt idGlobalDof = 0;
578 1 felInt sumPreviousSizeLocal = 0;
579 1 std::vector<felInt> loc2GlobPerUnknown;
580 felInt sizeIdDof;
581
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 m_mappingIdLocalDofPerUnknownToIdGlobalDofPetsc.resize(m_listUnknown.size());
582
583
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
 1 for ( int iProc = 0; iProc < MpiInfo::rankProc(); iProc++)
584 ✗ sumPreviousSizeLocal += m_localSizeOfVector[iProc];
585
586
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
 2 for ( std::size_t iUnknown = 0; iUnknown < m_listUnknown.size(); iUnknown++ )
587 1 m_numDofLocalUnknown[iUnknown] = 0;
588
589
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
 2 for ( std::size_t iUnknown = 0; iUnknown < m_listUnknown.size(); iUnknown++ ) {
590
2/2
✓ Branch 0 taken 4225 times.
✓ Branch 1 taken 1 times.
 4226 for ( felInt iDof = 0; iDof < m_numDofUnknown[iUnknown]; iDof++) {
591
1/2
✓ Branch 1 taken 4225 times.
✗ Branch 2 not taken.
 4225 AOApplicationToPetsc(m_ao, 1, &idGlobalDof);
592
4/8
✓ Branch 0 taken 4225 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 4225 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 4225 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 4225 times.
✗ Branch 9 not taken.
 4225 if ( (sumPreviousSizeLocal <= idGlobalDof) && (idGlobalDof < m_localSizeOfVector[MpiInfo::rankProc()]+sumPreviousSizeLocal) ) {
593 4225 m_numDofLocalUnknown[iUnknown]++;
594
1/2
✓ Branch 1 taken 4225 times.
✗ Branch 2 not taken.
 4225 loc2GlobPerUnknown.push_back(idGlobalDof);
595 }
596 4225 idGlobalDof++;
597 }
598 1 sizeIdDof = loc2GlobPerUnknown.size();
599
600
1/2
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
 1 IS_LOCAL_TO_GLOBAL_MAPPING_MACRO(m_petscComm, 1, sizeIdDof, loc2GlobPerUnknown.data(), PETSC_COPY_VALUES,
601 &m_mappingIdLocalDofPerUnknownToIdGlobalDofPetsc[iUnknown] );
602 1 loc2GlobPerUnknown.clear();
603 }
604 2 }
605
606 /***********************************************************************************/
607 /***********************************************************************************/
608
609 1 void EigenProblem::allocateSequentialSolution()
610 {
611 1 m_seqSol.createSeq(PETSC_COMM_SELF,m_numDof);
612
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_seqSol.set(0.0);
613 1 m_seqSolAllocated=true;
614 1 }
615
616 /***********************************************************************************/
617 /***********************************************************************************/
618
619 1 void EigenProblem::assembleMatrix()
620 {
621 ElementType eltType; // geometric element type in the mesh.
622 1 int numPointPerElt = 0; // number of points per geometric element.
623 1 int currentLabel = 0; // number of label domain.
624 1 felInt numEltPerLabel = 0; // number of element for one label and one eltType.
625
626 // use to define a "global" numbering of element in the mesh.
627 felInt numElement[ GeometricMeshRegion::m_numTypesOfElement ];
628
2/2
✓ Branch 0 taken 25 times.
✓ Branch 1 taken 1 times.
 26 for (int ityp=0; ityp<GeometricMeshRegion::m_numTypesOfElement; ityp++ ) {
629 25 eltType = (ElementType)ityp;
630 25 numElement[eltType] = 0;
631 }
632
633 // contains points of the current element.
634 1 std::vector<Point*> elemPoint;
635 // contains ids of point of the current element.
636 1 std::vector<felInt> elemIdPoint;
637 // use to get element number in support dof mesh ordering.
638 felInt ielSupportDof;
639
640 //Assembly loop.
641 //First loop on geometric type.
642 1 const std::vector<ElementType>& bagElementTypeDomain = m_meshLocal->bagElementTypeDomain();
643
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
 2 for (std::size_t i = 0; i < bagElementTypeDomain.size(); ++i) {
644 1 eltType = bagElementTypeDomain[i];
645 //resize array.
646 1 numPointPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
647
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 elemPoint.resize(numPointPerElt, nullptr);
648
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 elemIdPoint.resize(numPointPerElt, 0);
649 //define all current finite element use in the problem from data
650 //file cnfiguration and allocate 1 and elemVec (question: virtual ?).
651
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 defineFiniteElement(eltType);
652 //allocate array use for assemble with petsc.
653
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 allocateArrayForAssembleMatrix();
654 //virtual function use in derived problem to allocate elemenField necessary.
655
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 initPerElementType();
656 //second loop on region of the mesh.
657
2/2
✓ Branch 6 taken 1 times.
✓ Branch 7 taken 1 times.
 2 for(auto itRef = m_meshLocal->intRefToBegEndMaps[eltType].begin(); itRef != m_meshLocal->intRefToBegEndMaps[eltType].end(); itRef++) {
658 1 currentLabel = itRef->first;
659 1 numEltPerLabel = itRef->second.second;
660 //By default this virtual defnie variable m_currentLabel: (m_currentLabel=label).
661 //We can switch on label region with that and define some parameters.
662
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 initPerDomain(currentLabel);
663 //Third loop on element in the region with the type: eltType. ("real" loop on elements).
664
2/2
✓ Branch 0 taken 8192 times.
✓ Branch 1 taken 1 times.
 8193 for ( felInt iel = 0; iel < numEltPerLabel; iel++) {
665 // return each id of point of the element and coordinate in two arrays: elemPoint and elemIdPoint.
666
1/2
✓ Branch 1 taken 8192 times.
✗ Branch 2 not taken.
 8192 setElemPoint(eltType, numElement[eltType], elemPoint, elemIdPoint, &ielSupportDof);
667 //todo Jeremie 01 02 2012: assemble multiple matrix
668
2/2
✓ Branch 1 taken 16384 times.
✓ Branch 2 taken 8192 times.
 24576 for (std::size_t j = 0; j < m_Matrix.size(); j++)
669
1/2
✓ Branch 3 taken 16384 times.
✗ Branch 4 not taken.
 16384 m_elementMat[j]->zero();
670
1/2
✓ Branch 1 taken 8192 times.
✗ Branch 2 not taken.
 8192 computeElementArray(elemPoint, elemIdPoint, ielSupportDof);
671 // add values of elemMat in the global matrix: m_Matrix[0].
672
1/2
✓ Branch 1 taken 8192 times.
✗ Branch 2 not taken.
 8192 setValueMatrix(&ielSupportDof);
673 8192 numElement[eltType]++;
674 }
675 }
676 //desallocate array use for assemble with petsc.
677
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 desallocateArrayForAssembleMatrix();
678 }
679
680 //real assembling of the matrix manage by PETsC.
681
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 1 times.
 3 for (std::size_t i = 0; i < m_Matrix.size(); i++) {
682
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 m_Matrix[i].assembly(MAT_FINAL_ASSEMBLY);
683 }
684
685 //call with high level of verbose to print matrix in matlab format.
686
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 writeMatrixForMatlab(m_verbose);
687 1 }
688
689 /***********************************************************************************/
690 /***********************************************************************************/
691
692 ✗ void EigenProblem::assembleMatrixBD()
693 {
694 ElementType eltType; //geometric element type in the mesh.
695 ✗ int numPointPerElt = 0; //number of points per geometric element.
696 ✗ int currentLabel = 0; //number of label domain.
697 ✗ felInt numEltPerLabel = 0; //number of element for one label and one eltType.
698
699 // use to define a "global" numbering of element in the mesh.
700 felInt numElement[ GeometricMeshRegion::m_numTypesOfElement ];
701 ✗ for (int ityp=0; ityp<GeometricMeshRegion::m_numTypesOfElement; ityp++ ) {
702 ✗ eltType = (ElementType)ityp;
703 ✗ numElement[eltType] = 0;
704 }
705
706 // contains points of the current element.
707 ✗ std::vector<Point*> elemPoint;
708
709 // contains ids of point of the current element.
710 ✗ std::vector<felInt> elemIdPoint;
711
712 // contains the ids of the support elements of a mesh element
713 ✗ std::vector<felInt> vectorIdSupport;
714
715 // use to get element number in support dof mesh ordering.
716 felInt ielSupportDof;
717
718 //Assembly loop.
719 //First loop on geometric type.
720 ✗ const std::vector<ElementType>& bagElementTypeDomain = m_meshLocal->bagElementTypeDomain();
721 ✗ for (std::size_t i = 0; i < bagElementTypeDomain.size(); ++i) {
722 ✗ eltType = bagElementTypeDomain[i];
723
724 //resize array.
725 ✗ numPointPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
726 ✗ elemPoint.resize(numPointPerElt, nullptr);
727 ✗ elemIdPoint.resize(numPointPerElt, 0);
728
729 //define all current finite element use in the problem from data
730 //file configuration and allocate elemMat and elemVec (question: virtual ?).
731 ✗ defineCurvilinearFiniteElement(eltType);
732
733 //allocate array use for assemble with petsc.
734 ✗ allocateArrayForAssembleMatrixBD();
735
736 //virtual function use in derived problem to allocate elemenField necessary.
737 ✗ initPerElementTypeBD();
738
739 // use by user to add specific term (term source for example with elemField).
740 ✗ userElementInitBD();
741
742 //second loop on region of the mesh.
743 ✗ for(auto itRef = m_meshLocal->intRefToBegEndMaps[eltType].begin(); itRef != m_meshLocal->intRefToBegEndMaps[eltType].end(); itRef++) {
744 ✗ currentLabel = itRef->first;
745 ✗ numEltPerLabel = itRef->second.second;
746
747 //By default this virtual defnie variable m_currentLabel: (m_currentLabel=label).
748 //We can switch on label region with that and define some parameters.
749 ✗ initPerDomainBD(currentLabel);
750
751 //Third loop on element in the region with the type: eltType. ("real" loop on elements).
752 ✗ for ( felInt iel = 0; iel < numEltPerLabel; iel++) {
753 // return each id of point of the element and coordinate in two arrays: elemPoint and elemIdPoint.
754 ✗ setElemPoint(eltType, numElement[eltType], elemPoint, elemIdPoint, vectorIdSupport);
755
756 // loop over all the support elements
757 ✗ for (std::size_t it = 0; it < vectorIdSupport.size(); it++) {
758 // get the id of the support
759 ✗ ielSupportDof = vectorIdSupport[it];
760
761 //todo Jeremie 01 02 2012: assemble multiple matrix
762 ✗ for (std::size_t j = 0; j < m_Matrix.size(); j++)
763 ✗ m_elementMatBD[j]->zero();
764
765 // function call in derived problem to compute specific operators of the problem (Heat, N-S,...).
766 ✗ computeElementArrayBD(elemPoint, elemIdPoint, ielSupportDof);
767
768 // add values of elemMat in the global matrix: m_Matrix[0].
769 ✗ setValueMatrixBD(&ielSupportDof);
770 }
771 ✗ numElement[eltType]++;
772 }
773 }
774 //allocate array use for assemble with petsc.
775 ✗ desallocateArrayForAssembleMatrix();
776 }
777
778 ✗ for (std::size_t i = 0; i < m_Matrix.size(); i++) {
779 ✗ m_Matrix[i].assembly(MAT_FINAL_ASSEMBLY);
780 }
781 //call with high level of verbose to print matrix in matlab format.
782 ✗ writeMatrixForMatlab(m_verbose);
783 }
784
785 /***********************************************************************************/
786 /***********************************************************************************/
787
788 1 void EigenProblem::initElementArray()
789 {
790 int idVar;
791
792
1/2
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
 1 std::vector<std::size_t> numberCmp(m_listUnknown.size());
793
1/2
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
 1 std::vector<const CurBaseFiniteElement*> finiteElt(m_listUnknown.size());
794
795
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
 2 for (std::size_t n=0; n<m_listUnknown.size(); n++){
796 1 idVar = m_listUnknown.idVariable(n);
797 1 numberCmp[n] = m_listVariable[idVar].numComponent();
798 1 finiteElt[n] = m_listCurrentFiniteElement[idVar];
799 }
800
801
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 1 times.
 3 for (std::size_t i = 0; i < m_Matrix.size(); i++)
802
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 m_elementMat.push_back(felisce::make_shared<ElementMatrix>(finiteElt, numberCmp, numberCmp));
803 1 }
804
805 /***********************************************************************************/
806 /***********************************************************************************/
807
808 ✗ void EigenProblem::initElementArrayBD()
809 {
810 int idVar;
811
812 ✗ std::vector<std::size_t> numberCmp(m_listUnknown.size());
813 ✗ std::vector<const CurBaseFiniteElement*> finiteElt(m_listUnknown.size());
814 ✗ for (std::size_t n=0; n<m_listUnknown.size(); n++){
815 ✗ idVar = m_listUnknown.idVariable(n);
816 ✗ numberCmp[n] = m_listVariable[idVar].numComponent();
817 ✗ finiteElt[n] = m_listCurvilinearFiniteElement[idVar];
818 }
819
820 ✗ for (std::size_t i = 0; i < m_Matrix.size(); i++)
821 ✗ m_elementMatBD.push_back(felisce::make_shared<ElementMatrix>(finiteElt, numberCmp, numberCmp));
822 }
823
824 /***********************************************************************************/
825 /***********************************************************************************/
826
827 1 void EigenProblem::defineFiniteElement(const ElementType& eltType)
828 {
829
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
 1 if (m_listCurrentFiniteElement.size() > 0) {
830 ✗ m_listCurrentFiniteElement.clear();
831 //todo Jeremie 01 02 2012: assemble multiple matrix
832 ✗ m_elementMat.clear();
833 }
834 const RefElement *refEle;
835
836 1 int typeOfFiniteElement = 0;
837 CurrentFiniteElement* fe;
838 1 const GeoElement *geoEle = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second;
839
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
 2 for ( std::size_t iVar = 0, size = m_listVariable.size(); iVar < size; iVar++) {
840 1 typeOfFiniteElement = m_listVariable[iVar].finiteElementType(); // 0 = linear, 1 = quadratic, 2 = bubble
841 1 refEle = geoEle->defineFiniteEle(eltType,typeOfFiniteElement, *m_meshLocal);
842
1/2
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 fe = new CurrentFiniteElement(*refEle,*geoEle,m_listVariable[iVar].getDegreeOfExactness());
843 1 m_listCurrentFiniteElement.add(fe);
844 }
845 //Element matrix and std::vector initialisation
846 1 initElementArray();
847 1 }
848
849 /***********************************************************************************/
850 /***********************************************************************************/
851
852 ✗ void EigenProblem::defineCurvilinearFiniteElement(const ElementType& eltType)
853 {
854 ✗ if (m_listCurvilinearFiniteElement.size() > 0) {
855 ✗ m_listCurvilinearFiniteElement.clear();
856 ✗ std::cout << m_Matrix.size();
857 ✗ m_elementMatBD.clear();
858 }
859 //Initialisation of Finite Element curvlinear problem
860 const RefElement *refEle;
861 const GeoElement *geoEle;
862 ✗ int typeOfFiniteElement = 0;
863 CurvilinearFiniteElement* fe;
864 ✗ geoEle = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second;
865 ✗ for ( std::size_t iVar = 0; iVar < m_listVariable.size(); iVar++) {
866 ✗ typeOfFiniteElement = m_listVariable[iVar].finiteElementType();
867 ✗ refEle = geoEle->defineFiniteEle(eltType,typeOfFiniteElement, *m_meshLocal);
868 ✗ fe = new CurvilinearFiniteElement(*refEle,*geoEle,m_listVariable[iVar].degreeOfExactness());
869 ✗ FEL_ASSERT(fe);
870 ✗ m_listCurvilinearFiniteElement.add(fe);
871 }
872 //Element matrix and std::vector initialisation
873 ✗ initElementArrayBD();
874 }
875
876 /***********************************************************************************/
877 /***********************************************************************************/
878
879 8192 void EigenProblem::setElemPoint(ElementType& eltType, felInt iel, std::vector<Point*>& elemPoint,
880 std::vector<felInt>& elemIdPoint, felInt* ielSupportDof)
881 {
882 // int iVar = 0; // id of the variable. ielSupport always the same with all variable.
883 8192 int iUnknown = 0;
884 8192 m_supportDofUnknownLocal[iUnknown].getIdElementSupport(eltType, iel, *ielSupportDof);
885 8192 ISLocalToGlobalMappingApply(m_mappingElem,1,ielSupportDof,ielSupportDof);
886
887 8192 int numPointsPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
888 8192 m_meshLocal->getOneElement(eltType,iel,elemIdPoint,0);
889
2/2
✓ Branch 0 taken 24576 times.
✓ Branch 1 taken 8192 times.
 32768 for (int iPoint = 0; iPoint < numPointsPerElt; iPoint++) {
890 24576 elemPoint[iPoint] = &m_meshLocal->listPoints()[elemIdPoint[iPoint]]; // => m_meshLocal not const
891 }
892 8192 }
893
894 /***********************************************************************************/
895 /***********************************************************************************/
896
897 ✗ void EigenProblem::setElemPoint(ElementType& eltType, felInt iel, std::vector<Point*>& elemPoint,
898 std::vector<felInt>& elemIdPoint, std::vector<felInt>& vectorSupport)
899 {
900 // int iVar = 0; // id of the variable. ielSupport always the same with all variable.
901 // We assume that the support elements are the same for all the unknown.
902 ✗ int iUnknown = 0;
903 ✗ m_supportDofUnknownLocal[iUnknown].getIdElementSupport(eltType, iel, vectorSupport);
904 ✗ ISLocalToGlobalMappingApply(m_mappingElem, vectorSupport.size(), vectorSupport.data(), vectorSupport.data());
905 ✗ int numPointsPerElt = GeometricMeshRegion::m_numPointsPerElt[eltType];
906 ✗ m_meshLocal->getOneElement(eltType, iel, elemIdPoint, 0);
907 ✗ for (int iPoint = 0; iPoint < numPointsPerElt; iPoint++) {
908 ✗ elemPoint[iPoint] = &m_meshLocal->listPoints()[elemIdPoint[iPoint]]; // => m_meshLocal not const
909 }
910 }
911
912 /***********************************************************************************/
913 /***********************************************************************************/
914
915 8192 void EigenProblem::setValueMatrix(const felInt* ielSupportDof)
916 {
917 // \todo: works (a priori) only for the square matrix
918 8192 const felInt numDofTotal = m_elementMat[0]->mat().size1();
919 8192 felInt node1 = 0;
920 8192 felInt node2 = 0;
921 8192 felInt cptGposLine = 0;
922 8192 felInt cptGposLine2 = 0;
923 8192 felInt sizeSupport = 0;
924 8192 felInt sizeSupport2 = 0;
925 8192 felInt cptBlock = 0;
926 8192 int idVar = -1;
927 8192 int idVar2 = -1;
928
2/2
✓ Branch 1 taken 8192 times.
✓ Branch 2 taken 8192 times.
 16384 for ( std::size_t iUnknown = 0; iUnknown < m_listUnknown.size(); iUnknown++) {
929 8192 idVar = m_listUnknown.idVariable(iUnknown);
930 8192 sizeSupport = m_listCurrentFiniteElement[idVar]->numDof();
931 // Warning: in elemMat the order is e.g. ( ux, uy, p)
932 // in the matrix the order is (ux1,uy1,ux2,uy2,...,p1,p2,..)
933
2/2
✓ Branch 2 taken 8192 times.
✓ Branch 3 taken 8192 times.
 16384 for (std::size_t iComp = 0; iComp < m_listVariable[idVar].numComponent(); iComp++) {
934 8192 cptBlock++;
935
2/2
✓ Branch 0 taken 24576 times.
✓ Branch 1 taken 8192 times.
 32768 for ( int iSupport = 0; iSupport < sizeSupport; iSupport++) {
936
1/2
✓ Branch 1 taken 24576 times.
✗ Branch 2 not taken.
 24576 m_dof.loc2glob(*ielSupportDof, iSupport, idVar, iComp, node1);
937 24576 m_GposLine[cptGposLine]= node1;
938 24576 cptGposLine2 = 0;
939
2/2
✓ Branch 1 taken 24576 times.
✓ Branch 2 taken 24576 times.
 49152 for ( std::size_t iUnknown2 = 0; iUnknown2 <m_listUnknown.size(); iUnknown2++) {
940 24576 idVar2 = m_listUnknown.idVariable(iUnknown2);
941 24576 sizeSupport2 = m_listCurrentFiniteElement[idVar2]->numDof();
942 // order matters: comp before support to manage ( ux, uy, p) order in mat elem to (ux1,uy1,ux2,uy2,...,p1,p2,..)!!!!
943
2/2
✓ Branch 2 taken 24576 times.
✓ Branch 3 taken 24576 times.
 49152 for (std::size_t jComp = 0; jComp <m_listVariable[idVar2].numComponent(); jComp++) {
944
2/2
✓ Branch 0 taken 73728 times.
✓ Branch 1 taken 24576 times.
 98304 for (felInt jSupport = 0; jSupport < sizeSupport2; jSupport++) {
945
1/2
✓ Branch 1 taken 73728 times.
✗ Branch 2 not taken.
 73728 int iConnect = m_dof.getNumGlobComp( iUnknown, iComp);
946
1/2
✓ Branch 1 taken 73728 times.
✗ Branch 2 not taken.
 73728 int jConnect = m_dof.getNumGlobComp( iUnknown2, jComp);
947
2/4
✓ Branch 2 taken 73728 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 73728 times.
✗ Branch 5 not taken.
 73728 if ( m_listUnknown.mask()(iConnect, jConnect) > 0) {
948
1/2
✓ Branch 1 taken 73728 times.
✗ Branch 2 not taken.
 73728 m_dof.loc2glob(*ielSupportDof,jSupport,idVar2,jComp,node2);
949 73728 m_GposColumn[cptGposLine2] = node2;
950
2/2
✓ Branch 1 taken 147456 times.
✓ Branch 2 taken 73728 times.
 221184 for (std::size_t i = 0; i < m_Matrix.size(); i++) {
951
1/2
✓ Branch 4 taken 147456 times.
✗ Branch 5 not taken.
 147456 m_valueMat[i][cptGposLine2 + cptGposLine*numDofTotal] = m_elementMat[i]->mat()(cptGposLine, cptGposLine2);
952 }
953 73728 cptGposLine2++;
954 }
955 }
956 }
957 }
958 24576 cptGposLine++;
959 }
960 }
961 }
962
963
1/2
✓ Branch 1 taken 8192 times.
✗ Branch 2 not taken.
 8192 AOApplicationToPetsc(m_ao,cptGposLine,m_GposLine);
964
1/2
✓ Branch 1 taken 8192 times.
✗ Branch 2 not taken.
 8192 AOApplicationToPetsc(m_ao,cptGposLine2,m_GposColumn);
965
2/2
✓ Branch 1 taken 16384 times.
✓ Branch 2 taken 8192 times.
 24576 for (std::size_t i = 0; i < m_Matrix.size(); i++) {
966
1/2
✓ Branch 3 taken 16384 times.
✗ Branch 4 not taken.
 16384 m_Matrix[i].setValues(numDofTotal,m_GposLine,numDofTotal,m_GposColumn,m_valueMat[i],ADD_VALUES);
967 }
968 8192 }
969
970 /***********************************************************************************/
971 /***********************************************************************************/
972
973 1 void EigenProblem::allocateArrayForAssembleMatrix()
974 {
975 1 felInt numDofTotal = 0;
976 1 int idFE = -1; //(idFE = idVar)
977
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
 2 for (std::size_t iUnknown = 0; iUnknown < m_listUnknown.size(); iUnknown++) {
978 1 idFE = m_listUnknown.idVariable(iUnknown);
979 1 numDofTotal += m_listCurrentFiniteElement[idFE]->numDof()*m_listVariable[idFE].numComponent();
980 }
981
982
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if (m_GposColumn == nullptr)
983
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 m_GposColumn = new felInt[numDofTotal];
984
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 if (m_GposLine == nullptr)
985
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 m_GposLine = new felInt[numDofTotal];
986
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 1 times.
 3 for (std::size_t i = 0; i < m_Matrix.size(); i++) {
987
2/4
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
 2 m_valueMat.push_back(new double[numDofTotal*numDofTotal]);
988 }
989 1 }
990
991 /***********************************************************************************/
992 /***********************************************************************************/
993
994 ✗ void EigenProblem::setValueMatrixBD(const felInt* ielSupportDof)
995 {
996 ✗ const felInt numDofTotal = m_elementMatBD[0]->mat().size1();
997
998 ✗ felInt node1 = 0;
999 ✗ felInt node2 = 0;
1000 ✗ felInt cptGposLine = 0;
1001 ✗ felInt cptGposLine2 = 0;
1002 ✗ felInt sizeSupport = 0;
1003 ✗ felInt sizeSupport2 = 0;
1004 ✗ felInt cptBlock = 0;
1005 ✗ int idVar = -1;
1006 ✗ int idVar2 = -1;
1007
1008 ✗ for ( std::size_t iUnknown = 0; iUnknown < m_listUnknown.size(); iUnknown++) {
1009 ✗ idVar = m_listUnknown.idVariable(iUnknown);
1010 ✗ sizeSupport = m_listCurvilinearFiniteElement[idVar]->numDof();
1011 // order matters: comp before support to manage ( ux, uy, p) order in mat elem to (ux1,uy1,ux2,uy2,...,p1,p2,..)!!!!
1012 ✗ for (std::size_t iComp = 0; iComp < m_listVariable[idVar].numComponent(); iComp++) {
1013 ✗ cptBlock++;
1014 ✗ for ( int iSupport = 0; iSupport < sizeSupport; iSupport++) {
1015 ✗ m_dof.loc2glob(*ielSupportDof, iSupport, idVar, iComp, node1);
1016 ✗ m_GposLine[cptGposLine]= node1;
1017 ✗ cptGposLine2 = 0;
1018 ✗ for ( std::size_t iUnknown2 = 0; iUnknown2 <m_listUnknown.size(); iUnknown2++) {
1019 ✗ idVar2 = m_listUnknown.idVariable(iUnknown2);
1020 ✗ sizeSupport2 = m_listCurvilinearFiniteElement[idVar2]->numDof();
1021 // order matters: comp before support to manage ( ux, uy, p) order in mat elem to (ux1,uy1,ux2,uy2,...,p1,p2,..)!!!!
1022 ✗ for (std::size_t jComp = 0; jComp <m_listVariable[idVar2].numComponent(); jComp++) {
1023 ✗ for ( int jSupport = 0; jSupport < sizeSupport2; jSupport++) {
1024 ✗ int iConnect = m_dof.getNumGlobComp( iUnknown, iComp);
1025 ✗ int jConnect = m_dof.getNumGlobComp( iUnknown2, jComp);
1026 ✗ if ( m_listUnknown.mask()(iConnect, jConnect) > 0) {
1027 ✗ m_dof.loc2glob(*ielSupportDof,jSupport,idVar2,jComp,node2);
1028 ✗ m_GposColumn[cptGposLine2] = node2;
1029 ✗ for (std::size_t i = 0; i < m_Matrix.size(); i++) {
1030
1031 ✗ m_valueMat[i][cptGposLine2 + cptGposLine*numDofTotal] = m_elementMatBD[i]->mat()(cptGposLine, cptGposLine2);
1032 }
1033 ✗ cptGposLine2++;
1034 }
1035 }
1036 }
1037 }
1038 ✗ cptGposLine++;
1039 }
1040 }
1041 }
1042
1043 ✗ AOApplicationToPetsc(m_ao,cptGposLine,m_GposLine);
1044 ✗ AOApplicationToPetsc(m_ao,cptGposLine2,m_GposColumn);
1045 ✗ for (std::size_t i = 0; i < m_Matrix.size(); i++) {
1046 ✗ m_Matrix[i].setValues(numDofTotal,m_GposLine,numDofTotal,m_GposColumn,m_valueMat[i],ADD_VALUES);
1047 }
1048 }
1049
1050 /***********************************************************************************/
1051 /***********************************************************************************/
1052
1053 ✗ void EigenProblem::allocateArrayForAssembleMatrixBD()
1054 {
1055 ✗ felInt numDofTotal = 0;
1056 ✗ for (std::size_t iFe = 0; iFe < m_listCurvilinearFiniteElement.size(); iFe++) {
1057 ✗ numDofTotal += m_listCurvilinearFiniteElement[iFe]->numDof()*m_listVariable[iFe].numComponent();
1058 }
1059 ✗ m_GposColumn = new felInt[numDofTotal];
1060 ✗ m_GposLine = new felInt[numDofTotal];
1061
1062 ✗ for (std::size_t i = 0; i < m_Matrix.size(); i++) {
1063 ✗ m_valueMat.push_back(new double[numDofTotal*numDofTotal]);
1064 }
1065 }
1066
1067 /***********************************************************************************/
1068 /***********************************************************************************/
1069
1070 1 void EigenProblem::desallocateArrayForAssembleMatrix()
1071 {
1072
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 1 times.
 3 for (std::size_t i = 0; i < m_Matrix.size(); i++) {
1073
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 delete [] m_valueMat[i];
1074 }
1075 1 m_valueMat.clear();
1076
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 delete [] m_GposLine;
1077 1 m_GposLine = nullptr;
1078
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 delete [] m_GposColumn;
1079 1 m_GposColumn = nullptr;
1080 1 }
1081
1082 /***********************************************************************************/
1083 /***********************************************************************************/
1084
1085 ✗ void EigenProblem::scaleMatrix(felReal coef, int i)
1086 {
1087 ✗ m_Matrix[i].scale(coef);
1088 }
1089
1090 /***********************************************************************************/
1091 /***********************************************************************************/
1092
1093 ✗ void EigenProblem::scaleMatrix(felReal coef, PetscMatrix& M)
1094 {
1095 ✗ M.scale(coef);
1096 }
1097
1098 /***********************************************************************************/
1099 /***********************************************************************************/
1100
1101 ✗ void EigenProblem::printSolution(int verbose) const
1102 {
1103 ✗ if (verbose > 19) {
1104 ✗ PetscPrintf(m_petscComm,"\n/================Solution of the assembly system Ax=b================/\n");
1105 ✗ m_sol.view();
1106 }
1107 }
1108
1109 /***********************************************************************************/
1110 /***********************************************************************************/
1111
1112 ✗ void EigenProblem::printMatrix(PetscMatrix& M)
1113 {
1114 ✗ M.view();
1115 }
1116
1117 /***********************************************************************************/
1118 /***********************************************************************************/
1119
1120 ✗ void EigenProblem::printVector(PetscVector& V)
1121 {
1122 ✗ V.view();
1123 }
1124
1125 /***********************************************************************************/
1126 /***********************************************************************************/
1127
1128 1 void EigenProblem::writeMatrixForMatlab(int verbose) const
1129 {
1130
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 if (verbose > 20) {
1131 ✗ PetscPrintf(m_petscComm,"\n/===========Matrix write in Matlab file matrix.m=============/\n");
1132 PetscViewer viewer;
1133 ✗ PetscViewerCreate(m_petscComm,&viewer);
1134 ✗ PetscViewerASCIIOpen(m_petscComm,"matrix.m",&viewer);
1135 ✗ PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
1136 ✗ m_Matrix[0].view(viewer);
1137 }
1138 1 }
1139
1140 /***********************************************************************************/
1141 /***********************************************************************************/
1142
1143 ✗ void EigenProblem::writeMatrixForMatlab(std::string const& fileName, PetscMatrix& matrix) const
1144 {
1145 ✗ PetscPrintf(m_petscComm,"\n/===========Matrix write in Matlab file fileName=============/\n");
1146 PetscViewer viewer;
1147 ✗ PetscViewerCreate(m_petscComm,&viewer);
1148 ✗ PetscViewerASCIIOpen(m_petscComm,fileName.c_str(),&viewer);
1149 ✗ PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
1150 ✗ matrix.view(viewer);
1151 }
1152
1153 /***********************************************************************************/
1154 /***********************************************************************************/
1155
1156 ✗ void EigenProblem::writeVectorForMatlab(std::string const& fileName, PetscVector& vector) const
1157 {
1158 ✗ PetscPrintf(m_petscComm,"\n/===========std::vector write in Matlab file std::vector.m=============/\n");
1159 PetscViewer viewer;
1160 ✗ PetscViewerCreate(m_petscComm,&viewer);
1161 ✗ PetscViewerASCIIOpen(m_petscComm,fileName.c_str(),&viewer);
1162 ✗ PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
1163 ✗ vector.view(viewer);
1164 }
1165
1166 /***********************************************************************************/
1167 /***********************************************************************************/
1168
1169 ✗ void EigenProblem::buildSolver()
1170 {
1171 #ifdef FELISCE_WITH_SLEPC
1172 ✗ EPSSetOperators(m_eps,m_Matrix[0].toPetsc(),FELISCE_PETSC_NULLPTR);
1173 ✗ EPSSetDimensions(m_eps,1,PETSC_DECIDE,PETSC_DECIDE);
1174 ✗ EPSSetFromOptions(m_eps);
1175 #endif
1176 }
1177
1178 /***********************************************************************************/
1179 /***********************************************************************************/
1180
1181 ✗ void EigenProblem::buildSVD()
1182 {
1183 #ifdef FELISCE_WITH_SLEPC
1184 // TODO: Replace with this when SLEPc has SVDSetOperators in all versions
1185 // SVDSetOperators(m_svd,m_Matrix[0].toPetsc(), FELISCE_PETSC_NULLPTR);
1186 #if PETSC_VERSION_LT(3, 19, 0)
1187 SVDSetOperator(m_svd,m_Matrix[0].toPetsc());
1188 #else
1189 // TODO: Replace with this when SLEPc has SVDSetOperators in all versions
1190 ✗ SVDSetOperators(m_svd,m_Matrix[0].toPetsc(), FELISCE_PETSC_NULLPTR);
1191 #endif
1192 ✗ SVDSetType(m_svd,SVDCROSS);
1193 ✗ SVDSetFromOptions(m_svd);
1194 #endif
1195 }
1196
1197 /***********************************************************************************/
1198 /***********************************************************************************/
1199
1200 ✗ void EigenProblem::solve()
1201 {
1202 #ifdef FELISCE_WITH_SLEPC
1203 PetscInt nconv1;
1204 PetscInt nconv2;
1205 ✗ PetscReal sigma_1 = 0.;
1206 ✗ PetscReal sigma_n = 0.;
1207 //First request a singular value from one end of the spectrum
1208 ✗ EPSSetWhichEigenpairs(m_eps,EPS_LARGEST_MAGNITUDE);
1209 ✗ if (m_verbose >= 10) {
1210 ✗ m_Matrix[0].view();
1211 ✗ EPSView(m_eps,PETSC_VIEWER_STDOUT_WORLD);
1212 }
1213
1214 ✗ EPSSolve(m_eps);
1215
1216 //Get number of converged singular values
1217 ✗ EPSGetConverged(m_eps,&nconv1);
1218 //Get converged singular values: largest singular value is stored in sigma_1.
1219 ✗ if (nconv1 > 0) {
1220 ✗ EPSGetEigenpair(m_eps,0,&sigma_1,FELISCE_PETSC_NULLPTR,FELISCE_PETSC_NULLPTR,FELISCE_PETSC_NULLPTR);
1221 } else {
1222 ✗ PetscPrintf(PETSC_COMM_WORLD," Unable to compute large singular value!\n\n");
1223 }
1224
1225 //Request a singular value from the other end of the spectrum
1226 ✗ EPSSetWhichEigenpairs(m_eps,EPS_SMALLEST_MAGNITUDE);
1227 ✗ EPSSolve(m_eps);
1228 //Get number of converged eigenpairs
1229 ✗ EPSGetConverged(m_eps,&nconv2);
1230 //Get converged singular values: smallest singular value is stored in sigma_n.
1231 ✗ if (nconv2 > 0) {
1232 ✗ EPSGetEigenpair(m_eps,0,&sigma_n,FELISCE_PETSC_NULLPTR,FELISCE_PETSC_NULLPTR,FELISCE_PETSC_NULLPTR);
1233 } else {
1234 ✗ PetscPrintf(PETSC_COMM_WORLD," Unable to compute small singular value!\n\n");
1235 }
1236
1237 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1238 Display solution
1239 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1240 ✗ if (nconv1 > 0 && nconv2 > 0) {
1241 ✗ PetscPrintf(PETSC_COMM_WORLD," Computed singular values: sigma_1=%6F, sigma_n=%6F\n",sigma_1,sigma_n);
1242 ✗ PetscPrintf(PETSC_COMM_WORLD," Estimated condition number: sigma_1/sigma_n=%6F\n\n",sigma_1/sigma_n);
1243 }
1244 #endif
1245 }
1246
1247 /***********************************************************************************/
1248 /***********************************************************************************/
1249
1250 ✗ void EigenProblem::solveSVD()
1251 {
1252 #ifdef FELISCE_WITH_SLEPC
1253 ✗ SVDSetDimensions(m_svd,1,PETSC_DECIDE,PETSC_DECIDE);
1254 PetscInt nconv1;
1255 ✗ PetscReal sigma_1 = 0.;
1256 ✗ SVDSetWhichSingularTriplets(m_svd,SVD_LARGEST);
1257 ✗ SVDSolve(m_svd);
1258 ✗ SVDGetConverged(m_svd,&nconv1);
1259 ✗ if (nconv1 > 0) {
1260 ✗ SVDGetSingularTriplet(m_svd,0,&sigma_1,nullptr,nullptr);
1261 } else {
1262 ✗ PetscPrintf(PETSC_COMM_WORLD," Unable to compute large singular value!\n\n");
1263 }
1264
1265 PetscInt nconv2;
1266 ✗ PetscReal sigma_n = 0.;
1267 ✗ SVDSetWhichSingularTriplets(m_svd,SVD_SMALLEST);
1268 ✗ SVDSolve(m_svd);
1269 ✗ SVDGetConverged(m_svd,&nconv2);
1270 ✗ if (nconv2 > 0) {
1271 ✗ SVDGetSingularTriplet(m_svd,0,&sigma_n,nullptr,nullptr);
1272 } else {
1273 ✗ PetscPrintf(PETSC_COMM_WORLD," Unable to compute small singular value!\n\n");
1274 }
1275
1276 #endif
1277 }
1278
1279 /***********************************************************************************/
1280 /***********************************************************************************/
1281
1282 1 void EigenProblem::clearMatrix(const int i) {
1283 1 m_Matrix[i].setUnfactored();
1284 1 m_Matrix[i].zeroEntries();
1285 1 }
1286
1287 /***********************************************************************************/
1288 /***********************************************************************************/
1289
1290 1 void EigenProblem::gatherSolution() {
1291
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 PetscVector vout;
1292
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_sol.scatterToAll(vout,INSERT_VALUES,SCATTER_FORWARD);
1293
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 m_seqSol.copyFrom(vout);
1294
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 vout.destroy();
1295 1 }
1296
1297 /***********************************************************************************/
1298 /***********************************************************************************/
1299
1300 175 void EigenProblem::fromVecToDoubleStar(double* solution, PetscVector& sol, int rank, int dimension, felInt size) {
1301 IGNORE_UNUSED_ARGUMENT(dimension);
1302
1303
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 175 times.
 175 if (FelisceParam::verbose() >= 40 )
1304 ✗ sol.view();
1305 felInt sizeVec;
1306
1/2
✓ Branch 1 taken 175 times.
✗ Branch 2 not taken.
 175 sol.getSize(&sizeVec);
1307
1308 // Check if the demanded double* output has the same dimension of the input Vec
1309
1/2
✓ Branch 0 taken 175 times.
✗ Branch 1 not taken.
 175 if (size == 0)
1310 175 size = sizeVec;
1311
1312
1/2
✓ Branch 1 taken 175 times.
✗ Branch 2 not taken.
 175 PetscVector xsol;
1313
1/2
✓ Branch 1 taken 175 times.
✗ Branch 2 not taken.
 175 sol.scatterToZero(xsol, INSERT_VALUES, SCATTER_FORWARD);
1314
1315 175 felInt iout[size];
1316
2/2
✓ Branch 0 taken 739375 times.
✓ Branch 1 taken 175 times.
 739550 for (felInt i=0; i<size; i++)
1317 739375 iout[i] = i;
1318
1319
1/2
✓ Branch 1 taken 175 times.
✗ Branch 2 not taken.
 175 AOApplicationToPetsc(m_ao, size, iout);
1320
1321
1/2
✓ Branch 0 taken 175 times.
✗ Branch 1 not taken.
 175 if(rank == 0) {
1322
1/2
✓ Branch 1 taken 175 times.
✗ Branch 2 not taken.
 175 xsol.getValues( size, iout, solution);
1323 }
1324
1/2
✓ Branch 1 taken 175 times.
✗ Branch 2 not taken.
 175 xsol.destroy();
1325
1/2
✓ Branch 1 taken 175 times.
✗ Branch 2 not taken.
 350 }
1326
1327 /***********************************************************************************/
1328 /***********************************************************************************/
1329
1330 ✗ void EigenProblem::fromDoubleStarToVec(double* solution, PetscVector& sol, felInt size) {
1331 felInt ia;
1332 ✗ for (felInt i=0; i< size; i++) {
1333 ✗ ia = i;
1334 ✗ AOApplicationToPetsc(m_ao,1,&ia);
1335 ✗ sol.setValues(1,&ia,&solution[i],INSERT_VALUES);
1336 }
1337
1338 ✗ sol.assembly();
1339
1340 ✗ if (FelisceParam::verbose() >= 40 )
1341 ✗ sol.view();
1342 }
1343
1344 /***********************************************************************************/
1345 /***********************************************************************************/
1346
1347 ✗ void EigenProblem::writeSolution(int rank, IO& io, double& time, int iteration)
1348 {
1349 // to be sure that m_seqSol is up to date
1350 ✗ gatherSolution();
1351
1352 ✗ felInt iout[m_numDof];
1353 ✗ felReal valueVec[m_numDof];
1354 ✗ for (felInt i=0; i<m_numDof; i++)
1355 ✗ iout[i] = i;
1356 ✗ AOApplicationToPetsc(m_ao,m_numDof,iout);
1357 ✗ m_seqSol.getValues(m_numDof,iout,valueVec);
1358
1359 ✗ int idVar = -1;
1360 ✗ int idVar2 = -1;
1361 // Write solution in ensight format
1362 ✗ if ( rank == 0) {
1363 double* solution;
1364 ✗ for ( std::size_t iUnknown = 0; iUnknown < m_listUnknown.size(); iUnknown++) {
1365 ✗ felInt numValueStart = 0;
1366 ✗ idVar = m_listUnknown.idVariable(iUnknown);
1367 ✗ for ( std::size_t iUnknown2 = 0; iUnknown2 < iUnknown; iUnknown2++) {
1368 ✗ idVar2 = m_listUnknown.idVariable(iUnknown2);
1369 ✗ numValueStart += m_supportDofUnknown[iUnknown2].listNode().size()*m_listVariable[idVar2].numComponent();
1370 }
1371 ✗ switch (m_listVariable[idVar].numComponent()) {
1372 ✗ case 1:
1373 ✗ solution = new double[m_supportDofUnknown[iUnknown].listNode().size()];
1374 ✗ for (felInt i = 0; i < static_cast<felInt>(m_supportDofUnknown[iUnknown].listNode().size()); i++) {
1375 ✗ solution[i] = valueVec[i+numValueStart];
1376 }
1377 ✗ io.writeSolution(rank, time, iteration, 0, m_listVariable[idVar].name(), solution, m_supportDofUnknown[iUnknown].listNode().size(), m_mesh->domainDim(),true);
1378 ✗ break;
1379 ✗ case 2:
1380 ✗ solution = new double[m_supportDofUnknown[iUnknown].listNode().size()*3];
1381 ✗ for (felInt i = 0; i < static_cast<felInt>(m_supportDofUnknown[iUnknown].listNode().size()); i++) {
1382 ✗ solution[i*3] = valueVec[i*2+numValueStart];
1383 ✗ solution[i*3+1] = valueVec[i*2+1+numValueStart];
1384 ✗ solution[i*3+2] = 0.;
1385 }
1386 ✗ io.writeSolution(rank, time, iteration, 1, m_listVariable[idVar].name(), solution, m_supportDofUnknown[iUnknown].listNode().size(), m_mesh->domainDim(),true);
1387 ✗ break;
1388 ✗ case 3:
1389 ✗ solution = new double[m_supportDofUnknown[iUnknown].listNode().size()*3];
1390 ✗ for (felInt i = 0; i < static_cast<felInt>(m_supportDofUnknown[iUnknown].listNode().size()); i++) {
1391 ✗ solution[i*3] = valueVec[i*3+numValueStart];
1392 ✗ solution[i*3+1] = valueVec[i*3+1+numValueStart];
1393 ✗ solution[i*3+2] = valueVec[i*3+2+numValueStart];
1394 }
1395 ✗ io.writeSolution(rank, time, iteration, 1, m_listVariable[idVar].name(), solution, m_supportDofUnknown[iUnknown].listNode().size(), m_mesh->domainDim(),true);
1396 ✗ break;
1397 ✗ default:
1398 ✗ FEL_ERROR("Problem with the dimension of your solution to write it.")
1399 ✗ break;
1400 }
1401 //delete [] solution;
1402 }
1403 }
1404 }
1405
1406 /***********************************************************************************/
1407 /***********************************************************************************/
1408
1409 175 void EigenProblem::writeEnsightVector(double* solValue, int idIter,
1410 std::string varName)
1411 {
1412 int rankProc;
1413
1/2
✓ Branch 1 taken 175 times.
✗ Branch 2 not taken.
 175 MPI_Comm_rank(m_petscComm,&rankProc);
1414
1415 175 std::string iteration;
1416
2/2
✓ Branch 0 taken 70 times.
✓ Branch 1 taken 105 times.
 175 if (idIter < 10)
1417
2/4
✓ Branch 1 taken 70 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 70 times.
✗ Branch 5 not taken.
 70 iteration = "0000" + std::to_string(idIter);
1418
1/2
✓ Branch 0 taken 105 times.
✗ Branch 1 not taken.
 105 else if (idIter < 100)
1419
2/4
✓ Branch 1 taken 105 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 105 times.
✗ Branch 5 not taken.
 105 iteration = "000" + std::to_string(idIter);
1420 ✗ else if (idIter < 1000)
1421 ✗ iteration = "00" + std::to_string(idIter);
1422 ✗ else if (idIter < 10000)
1423 ✗ iteration = "0" + std::to_string(idIter);
1424 ✗ else if (idIter < 100000)
1425 ✗ iteration = std::to_string(idIter);
1426
1427
6/12
✓ Branch 1 taken 175 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 175 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 175 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 175 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 175 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 175 times.
✗ Branch 17 not taken.
 350 std::string fileName = FelisceParam::instance().resultDir + "/" + varName + "." + iteration + ".scl";
1428 FILE * solFile;
1429
1/2
✓ Branch 0 taken 175 times.
✗ Branch 1 not taken.
 175 if (rankProc == 0) {
1430
1/2
✓ Branch 2 taken 175 times.
✗ Branch 3 not taken.
 175 solFile = fopen(fileName.c_str(),"w");
1431
1/2
✓ Branch 1 taken 175 times.
✗ Branch 2 not taken.
 175 fprintf(solFile, "Scalar per node\n");
1432 175 int count = 0;
1433
2/2
✓ Branch 0 taken 739375 times.
✓ Branch 1 taken 175 times.
 739550 for(int j=0; j < m_numDof; j++) {
1434
1/2
✓ Branch 1 taken 739375 times.
✗ Branch 2 not taken.
 739375 fprintf(solFile,"%12.5e", solValue[j]);
1435 739375 count++;
1436
2/2
✓ Branch 0 taken 123200 times.
✓ Branch 1 taken 616175 times.
 739375 if(count == 6) {
1437
1/2
✓ Branch 1 taken 123200 times.
✗ Branch 2 not taken.
 123200 fprintf(solFile,"\n");
1438 123200 count = 0;
1439 }
1440 }
1441
1/2
✓ Branch 1 taken 175 times.
✗ Branch 2 not taken.
 175 fprintf(solFile,"\n");
1442
1/2
✓ Branch 1 taken 175 times.
✗ Branch 2 not taken.
 175 fclose(solFile);
1443 }
1444 175 }
1445
1446 /***********************************************************************************/
1447 /***********************************************************************************/
1448
1449 1 void EigenProblem::writeEnsightCase(const int numIt,
1450 const double dt, std::string varName, double time0)
1451 {
1452 FILE * pFile;
1453
4/8
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
 2 std::string caseName = FelisceParam::instance().resultDir + "/" + varName + ".case";
1454
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 pFile = fopen(caseName.c_str(),"w");
1455
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 fprintf(pFile,"FORMAT\n");
1456
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 fprintf(pFile,"type: ensight\n");
1457
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 fprintf(pFile, "GEOMETRY\n");
1458
3/6
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
 1 std::string nameGeo = "model: 1 " + FelisceParam::instance().outputMesh[0] + "\n";
1459
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 fprintf(pFile,"%s", nameGeo.c_str());
1460
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 fprintf(pFile, "VARIABLE\n");
1461
4/8
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
 2 std::string secName = "scalar per node: 1 " + varName + " " + varName +".*****.scl\n";
1462
1463
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 fprintf(pFile, "%s", secName.c_str());
1464
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 fprintf(pFile, "TIME\n");
1465
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 fprintf(pFile, "time std::set: %d \n", 1);
1466
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 fprintf(pFile, "number of steps: %d \n", numIt);
1467
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 fprintf(pFile, "filename start number: %d\n", 0);
1468
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 fprintf(pFile, "filename increment: %d\n", 1);
1469
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 fprintf(pFile, "time values:\n");
1470
1471 1 int count=0;
1472
2/2
✓ Branch 0 taken 25 times.
✓ Branch 1 taken 1 times.
 26 for(int j=0; j < numIt; j++) {
1473 25 double valT = time0 + j*dt;
1474
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 fprintf(pFile,"%lf ",valT);
1475 25 count++;
1476
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 21 times.
 25 if(count == 6) {
1477
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 fprintf(pFile,"\n");
1478 4 count=0;
1479 }
1480 }
1481
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 fclose(pFile);
1482 1 }
1483
1484 /***********************************************************************************/
1485 /***********************************************************************************/
1486
1487 1 void EigenProblem::writeEnsightCase(const int numIt,
1488 const double dt,
1489 std::vector<std::string> varName,
1490 double time0)
1491 {
1492 FILE * pFile;
1493
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 std::string caseName = FelisceParam::instance().resultDir + "/"
1494
2/4
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
 3 + varName[0] + ".case";
1495
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 pFile = fopen(caseName.c_str(),"w");
1496
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 fprintf(pFile,"FORMAT\n");
1497
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 fprintf(pFile,"type: ensight\n");
1498
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 fprintf(pFile, "GEOMETRY\n");
1499
3/6
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
 1 std::string nameGeo = "model: 1 " + FelisceParam::instance().outputMesh[0] + "\n";
1500
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 fprintf(pFile,"%s", nameGeo.c_str());
1501
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 fprintf(pFile, "VARIABLE\n");
1502
1503 1 std::string secName;
1504
2/2
✓ Branch 1 taken 6 times.
✓ Branch 2 taken 1 times.
 7 for (std::size_t i=0; i<varName.size(); i++) {
1505
4/8
✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 6 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 6 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 6 times.
✗ Branch 13 not taken.
 6 secName = "scalar per node: 1 " + varName[i] + " " + varName[i] +".*****.scl\n";
1506
1/2
✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
 6 fprintf(pFile, "%s", secName.c_str());
1507 }
1508
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 fprintf(pFile, "TIME\n");
1509
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 fprintf(pFile, "time std::set: %d \n", 1);
1510
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 fprintf(pFile, "number of steps: %d \n", numIt);
1511
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 fprintf(pFile, "filename start number: %d\n", 0);
1512
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 fprintf(pFile, "filename increment: %d\n", 1);
1513
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 fprintf(pFile, "time values:\n");
1514
1515 1 int count=0;
1516
2/2
✓ Branch 0 taken 25 times.
✓ Branch 1 taken 1 times.
 26 for(int j=0; j < numIt; j++) {
1517 25 double valT = time0 + j*dt;
1518
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 25 fprintf(pFile,"%lf ",valT);
1519 25 count++;
1520
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 21 times.
 25 if(count == 6) {
1521
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 fprintf(pFile,"\n");
1522 4 count=0;
1523 }
1524 }
1525
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 fclose(pFile);
1526 1 }
1527
1528 /***********************************************************************************/
1529 /***********************************************************************************/
1530
1531 ✗ void EigenProblem::writeGeoForUnknown(int iUnknown)
1532 {
1533 // Intialization
1534 // =============
1535 ✗ std::string fileNameGeo, fileNameRoot, fileName;
1536 ElementType eltType;
1537 ✗ Ensight ens(FelisceParam::instance().inputDirectory,FelisceParam::instance().inputFile,FelisceParam::instance().inputMesh[0],FelisceParam::instance().outputMesh[0],FelisceParam::instance().meshDir,FelisceParam::instance().resultDir);
1538 ✗ std::string keyWord;
1539 int the_ref;
1540 ✗ short verbose = 0;
1541 FILE * pFile;
1542 // int nVpEl = 0;
1543 // felInt startindex = 0;
1544 ✗ felInt numElemsPerRef = 0;
1545 ✗ felInt ielSupportDof = 0;
1546 ✗ std::vector<felInt> elem;
1547 felInt numElement[ GeometricMeshRegion::m_numTypesOfElement ];
1548 ✗ for (int ityp=0; ityp<GeometricMeshRegion::m_numTypesOfElement; ityp++ ) {
1549 ✗ eltType = (ElementType)ityp;
1550 ✗ numElement[eltType] = 0;
1551 }
1552 ✗ const int idVar = m_listUnknown.idVariable(iUnknown);
1553 ✗ fileNameGeo = listVariable().listVariable()[idVar].name()+".geo";
1554 ✗ fileNameRoot = FelisceParam::instance().resultDir;
1555 ✗ fileName = fileNameRoot + fileNameGeo;
1556 ✗ if (verbose>0) std::cout << "Writing file "<< fileName << std::endl;
1557
1558 ✗ pFile = fopen (fileName.c_str(),"w");
1559 ✗ if (pFile==nullptr) {
1560 ✗ std::string command = "mkdir -p " + FelisceParam::instance().resultDir;
1561 ✗ int ierr = system(command.c_str());
1562 ✗ if (ierr>0) {
1563 ✗ FEL_ERROR("Impossible to write "+fileName+" geo mesh");
1564 }
1565 ✗ pFile = fopen (fileName.c_str(),"w");
1566 }
1567 //! Descriptions lines (All Geometry files must contain these first six lines)
1568 ✗ fprintf(pFile,"Geometry file\n");
1569 ✗ fprintf(pFile,"Geometry file\n");
1570 ✗ fprintf(pFile,"node id assign\n");
1571 ✗ fprintf(pFile,"element id assign\n");
1572 ✗ fprintf(pFile,"coordinates\n");
1573 ✗ fprintf(pFile,"%8d\n",(int)m_supportDofUnknown[iUnknown].listNode().size());
1574
1575 //! write the vertices
1576 ✗ if ( m_mesh->numCoor() == 3 ) {
1577 ✗ for(felInt iv=0; iv< static_cast<felInt>(m_supportDofUnknown[iUnknown].listNode().size()); iv++) {
1578 ✗ fprintf( pFile,"%12.5e%12.5e%12.5e\n", m_supportDofUnknown[iUnknown].listNode()[iv].x(),
1579 ✗ m_supportDofUnknown[iUnknown].listNode()[iv].y(),m_supportDofUnknown[iUnknown].listNode()[iv].z() );
1580 }
1581 ✗ } else if ( m_mesh->numCoor() == 2 ) {
1582 ✗ double zcoor = 0.;
1583 ✗ for(felInt iv=0; iv< static_cast<felInt>(m_supportDofUnknown[iUnknown].listNode().size()); iv++) {
1584 ✗ fprintf( pFile, "%12.5e%12.5e%12.5e\n", m_supportDofUnknown[iUnknown].listNode()[iv].x(), m_supportDofUnknown[iUnknown].listNode()[iv].y(), zcoor);
1585 }
1586 } else {
1587 ✗ std::ostringstream msg;
1588 msg << "Error in writing geo File " << fileNameGeo << ". "
1589 ✗ << "Expecting mesh with 2 or 3 coordinates (it is " << m_mesh->numCoor() << ")."
1590 ✗ << std::endl;
1591 ✗ FEL_ERROR(msg.str());
1592 }
1593
1594 //!write the elements per references.
1595 //------------
1596 ✗ felInt cptpart = 0;
1597 ✗ std::string DescriptionLine,felName;
1598 ✗ GeometricMeshRegion::IntRefToElementType_type& intRefToEnum = m_mesh->intRefToEnum();
1599 ✗ for(auto it_ref = intRefToEnum.begin();
1600 ✗ it_ref != intRefToEnum.end(); ++it_ref) {
1601 ✗ the_ref = (*it_ref).first;
1602
1603 ✗ for(auto it_elem = it_ref->second.begin();
1604 ✗ it_elem != it_ref->second.end(); ++it_elem) {
1605 ✗ eltType = *it_elem;
1606 ✗ if (m_mesh->flagFormatMesh() == GeometricMeshRegion::FormatMedit) {
1607 ✗ cptpart++;
1608 ✗ fprintf ( pFile, "part%8d\n",cptpart );
1609 ✗ if ( m_listVariable[idVar].finiteElementType() == 0 )
1610 ✗ felName = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].first;
1611 ✗ else if ( m_listVariable[idVar].finiteElementType() == 1 )
1612 ✗ felName = GeometricMeshRegion::eltEnumToFelNameGeoEle[GeometricMeshRegion::eltLinearToEltQuad[eltType]].first;
1613 else {
1614 ✗ if(m_mesh->domainDim() == GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].second->numCoor())
1615 ✗ felName = GeometricMeshRegion::eltEnumToFelNameGeoEle[GeometricMeshRegion::eltLinearToEltBubble[eltType]].first;
1616 else
1617 ✗ felName = GeometricMeshRegion::eltEnumToFelNameGeoEle[eltType].first;
1618 }
1619
1620 ✗ auto it_felName = FelisceParam::instance().felName2MapintRef2DescLine.find(felName);
1621 ✗ if ( it_felName == FelisceParam::instance().felName2MapintRef2DescLine.end() ) {
1622 //std::cout << "Warning: ----------->trying to access an unexisting element: " << felName << " in input file !!!!!!!" << std::endl;
1623 ✗ DescriptionLine = felName + "m_ref_" + std::to_string(the_ref);
1624 } else {
1625 ✗ auto it_ref2 = it_felName->second.find(the_ref);
1626 ✗ if ( it_ref2 == it_felName->second.end() ) {
1627 /*std::cout << "Warning: trying to access an unexisting reference: " << the_ref << "for element : "
1628 << felName << " in input file !!!!! " << std::endl << std::flush;*/
1629 ✗ DescriptionLine = felName + "m_ref_" + std::to_string(the_ref);
1630 } else {
1631 ✗ DescriptionLine = it_ref2->second;
1632 }
1633 }
1634 ✗ fprintf (pFile, "%s", DescriptionLine.c_str() );
1635 }
1636 ✗ fprintf(pFile,"\n");
1637 ✗ keyWord = ens.eltFelNameToEns6Pair()[felName].second;
1638 ✗ fprintf( pFile, "%s", keyWord.c_str() );
1639 ✗ fprintf(pFile,"\n");
1640 ✗ numElemsPerRef =m_mesh->intRefToBegEndMaps[eltType][the_ref].second;
1641 // startindex =m_mesh->intRefToBegEndMaps[eltType][the_ref].first;
1642 ✗ fprintf(pFile,"%8d\n", numElemsPerRef);
1643 // nVpEl = GeometricMeshRegion::m_numPointsPerElt[eltType];
1644 ✗ for ( felInt iel = 0; iel < numElemsPerRef; iel++) {
1645 ✗ m_supportDofUnknown[iUnknown].getIdElementSupport(eltType, numElement[eltType], ielSupportDof);
1646 ✗ elem.resize(m_supportDofUnknown[iUnknown].getNumSupportDof(ielSupportDof),0);
1647 ✗ for ( int iSupport = 0; iSupport < m_supportDofUnknown[iUnknown].getNumSupportDof(ielSupportDof); iSupport++) {
1648 ✗ elem[iSupport] = m_supportDofUnknown[iUnknown].iSupportDof()[ m_supportDofUnknown[iUnknown].iEle()[ielSupportDof] + iSupport];
1649 }
1650 ✗ for (std::size_t ivert = 0; ivert < elem.size(); ivert++) {
1651 ✗ fprintf( pFile, "%8d", elem[ivert]+1);
1652 }
1653 ✗ fprintf(pFile,"\n");
1654 ✗ numElement[eltType]++;
1655 }
1656 }
1657 }
1658 ✗ fclose(pFile);
1659 }
1660
1661 /***********************************************************************************/
1662 /***********************************************************************************/
1663
1664 ✗ void EigenProblem::findCoordinateWithIdDof(felInt i, Point& pt)
1665 {
1666 ✗ felInt inf = 0;
1667 ✗ felInt sup = 0;
1668 ✗ int idUnknown = -1;
1669 ✗ int idVar = -1;
1670 ✗ bool findUnknown = false;
1671 ✗ for (std::size_t iUnknown = 0; iUnknown < m_listUnknown.size() && findUnknown == false; iUnknown++) {
1672 ✗ if (iUnknown == 0) {
1673 ✗ inf = 0;
1674 ✗ sup = m_numDofUnknown[iUnknown];
1675 } else {
1676 ✗ inf += m_numDofUnknown[iUnknown-1];
1677 ✗ if (iUnknown == m_listUnknown.size()-1)
1678 ✗ sup = m_numDof;
1679 else
1680 ✗ sup += m_numDofUnknown[iUnknown];
1681 }
1682
1683 ✗ if (inf <= i && i < sup) {
1684 ✗ idUnknown = iUnknown;
1685 ✗ findUnknown = true;
1686 }
1687 }
1688 ✗ const felInt idDofPerUnknown = i - inf;
1689 ✗ felInt idSupportDof = -1;
1690 ✗ idVar = m_listUnknown.idVariable(idUnknown);
1691 ✗ idSupportDof = idDofPerUnknown/m_listVariable[idVar].numComponent();
1692 ✗ pt = m_supportDofUnknown[idUnknown].listNode()[idSupportDof];
1693 }
1694
1695 }
1696