GCC Code Coverage Report

Directory: ./
File: Solver/linearProblemHyperElasticity.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 271 311 87.1%
Branches: 172 374 46.0%
Line Branch Exec Source
1 // ______ ______ _ _ _____ ______
2 // | ____| ____| | (_)/ ____| | ____|
3 // | |__ | |__ | | _| (___ ___| |__
4 // | __| | __| | | | |\___ \ / __| __|
5 // | | | |____| |____| |____) | (__| |____
6 // |_| |______|______|_|_____/ \___|______|
7 // Finite Elements for Life Sciences and Engineering
8 //
9 // License: LGL2.1 License
10 // FELiScE default license: LICENSE in root folder
11 //
12 // Main authors: S. Gilles
13 //
14
15 // System includes
16
17 // External includes
18
19 // Project includes
20 #include "Solver/linearProblemHyperElasticity.hpp"
21 #include "DegreeOfFreedom/boundaryCondition.hpp"
22 #include "Core/felisceTransient.hpp"
23
24 namespace felisce
25 {
26
27 /**
28 * \brief This function is given to SNES to compute the function.
29 *
30 * However as a matter of fact the jacobian will also be computed here, as it is the way Felisce works
31 * (it hence avoids duplication or storage of intermediate matrixes)
32 *
33 * \param[in] snes SNES object. Not explicitly used but required by Petsc prototype.
34 * \param[in] evaluationState State at which to evaluate residual.
35 * \param[in] residual Vector to put residual. What is used as residual in LinearProblem is #_RHS.
36 * \param[in] context_as_void Optional user-defined function context. In our case a #SnesContext object,
37 * of which LinearProblem object is an attribute.
38 */
39 1044 PetscErrorCode snesFunctionImpl(SNES snes, Vec evaluationState, Vec residual, void* context_as_void) {
40 static_cast<void>(snes);
41 // static_cast<void>(residual); // residual as input argument should in fact be the same as linear_problem.RHS()
42
43
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1044 times.
 1044 FEL_ASSERT(context_as_void);
44
45 1044 SnesContext* context_ptr = static_cast<SnesContext*>(context_as_void);
46 1044 SnesContext& context = *context_ptr;
47 1044 LinearProblem* linear_problem = context.linearProblem;
48
49 {
50 // DEV TMP: not generic enough! (displacement is really related to my own peculiar problem...)
51
1/2
✓ Branch 1 taken 1044 times.
✗ Branch 2 not taken.
 1044 PetscInt NdofLocalPerUnknown = linear_problem->numDofLocalPerUnknown(/*PhysicalVariable::*/displacement);
52
53 1044 PetscVector& lp_evaluation_state = linear_problem->evaluationState();
54 1044 PetscVector& lp_rhs = linear_problem->vector();
55
56
1/2
✓ Branch 2 taken 1044 times.
✗ Branch 3 not taken.
 1044 std::vector<PetscInt> local_indexing(NdofLocalPerUnknown); // std::vector that contains integers from 0 to local size - 1
57 {
58 // std::iota(local_indexing.begin(), local_indexing.end(), 0); // std::iota only works with C++11!
59
60
2/2
✓ Branch 0 taken 2161602 times.
✓ Branch 1 taken 1044 times.
 2162646 for (int i = 0; i < NdofLocalPerUnknown; ++i) // substitute that also works in C++ 98 (arguably more easy to understand anyway)
61 2161602 local_indexing[i] = i;
62 }
63
64
1/2
✓ Branch 2 taken 1044 times.
✗ Branch 3 not taken.
 1044 std::vector<PetscInt> global_indexing(NdofLocalPerUnknown);
65
66
2/4
✓ Branch 2 taken 1044 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1044 times.
✗ Branch 6 not taken.
 2088 ISLocalToGlobalMappingApply(linear_problem->mappingDofLocalToDofGlobal(/*PhysicalVariable::*/displacement),
67 NdofLocalPerUnknown,
68 1044 local_indexing.data(),
69 global_indexing.data());
70
71 {
72 PetscScalar * evaluation_state_array;
73
1/2
✓ Branch 1 taken 1044 times.
✗ Branch 2 not taken.
 1044 VecGetArrayRead(evaluationState,(const PetscScalar**)&evaluation_state_array);
74
75 PetscScalar * residual_array;
76
1/2
✓ Branch 1 taken 1044 times.
✗ Branch 2 not taken.
 1044 VecGetArray(residual,&residual_array);
77
78
1/2
✓ Branch 2 taken 1044 times.
✗ Branch 3 not taken.
 1044 FEL_ASSERT_EQUAL(global_indexing.size(), NdofLocalPerUnknown);
79
80
1/2
✓ Branch 3 taken 1044 times.
✗ Branch 4 not taken.
 1044 lp_evaluation_state.setValues(global_indexing.size(),global_indexing.data(),evaluation_state_array,INSERT_VALUES);
81
1/2
✓ Branch 1 taken 1044 times.
✗ Branch 2 not taken.
 1044 lp_evaluation_state.assembly();
82
83
1/2
✓ Branch 3 taken 1044 times.
✗ Branch 4 not taken.
 1044 lp_rhs.setValues(global_indexing.size(),global_indexing.data(),residual_array,INSERT_VALUES);
84
1/2
✓ Branch 1 taken 1044 times.
✗ Branch 2 not taken.
 1044 lp_rhs.assembly();
85
86
1/2
✓ Branch 1 taken 1044 times.
✗ Branch 2 not taken.
 1044 VecRestoreArray(residual,&residual_array);
87
1/2
✓ Branch 1 taken 1044 times.
✗ Branch 2 not taken.
 1044 VecRestoreArrayRead(evaluationState,(const PetscScalar**)&evaluation_state_array);
88
89 }
90 1044 }
91
92 1044 linear_problem->preSNESAssemble();
93 1044 linear_problem->assembleMatrixRHS(context.rankProc, context.flag_matrix_rhs);
94 1044 linear_problem->computeSystemAfterAssembling();
95
96 //Apply boundary conditions.
97 1044 linear_problem->finalizeEssBCTransient();
98 {
99 enum { idBCforLinCombMethod = 0 };
100 1044 bool doPrintBC = false;
101 1044 linear_problem->applyBC(FelisceParam::instance().essentialBoundaryConditionsMethod,
102 context.rankProc,
103 FlagMatrixRHS::matrix_and_rhs,
104 FlagMatrixRHS::matrix_and_rhs,
105 idBCforLinCombMethod, doPrintBC,
106 context.do_apply_natural);
107 }
108
109 1044 VecCopy(linear_problem->vector().toPetsc(),residual);
110 //VecView(residual,PETSC_VIEWER_STDOUT_WORLD);
111 1044 PetscFunctionReturn(0);
112 }
113
114 /***********************************************************************************/
115 /***********************************************************************************/
116
117 /**
118 * \brief This function is supposed to be given to SNES to compute the jacobian.
119 *
120 * However the jacobian has already been computed within snesFunctionImpl, so this function is mostly a dummy one...
121 */
122 700 PetscErrorCode snesJacobianImpl(SNES snes, Vec evaluationState, Mat jacobian , Mat preconditioner, void* context_as_void) {
123 static_cast<void>(snes);
124 static_cast<void>(jacobian);
125 static_cast<void>(preconditioner);
126 static_cast<void>(evaluationState);
127
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 700 times.
 700 FEL_ASSERT(context_as_void);
128
129 700 PetscFunctionReturn(0);
130 }
131
132 /***********************************************************************************/
133 /***********************************************************************************/
134
135 24 LinearProblemHyperElasticity::LinearProblemHyperElasticity(SpatialForce force)
136 : LinearProblem("HyperElasticity equation"),
137
1/2
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
 24 m_force(force),
138 24 system_(HyperelasticityNS::static_),
139
8/16
✓ Branch 2 taken 24 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 24 times.
✗ Branch 6 not taken.
✓ Branch 10 taken 24 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 24 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 24 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 24 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 24 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 24 times.
✗ Branch 26 not taken.
 48 is_first_dynamic_iteration_done_(false)
140 {
141 24 }
142
143 /***********************************************************************************/
144 /***********************************************************************************/
145
146 96 LinearProblemHyperElasticity::~LinearProblemHyperElasticity()
147 = default;
148
149 /***********************************************************************************/
150 /***********************************************************************************/
151
152 24 void LinearProblemHyperElasticity::initialize(std::vector<GeometricMeshRegion::Pointer>& mesh, FelisceTransient::Pointer fstransient, MPI_Comm& comm, bool doUseSNES)
153 {
154 // Retrieve the instance
155
1/2
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
 24 const auto& r_instance = FelisceParam::instance();
156
157 using namespace HyperelasticityNS;
158
159
1/2
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
 24 LinearProblem::initialize(mesh, comm, doUseSNES);
160 24 m_fstransient = fstransient;
161
162
2/4
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 24 times.
✗ Branch 5 not taken.
 24 m_typeOfForceApplied = Private::HyperElasticity::AssignForceType(FelisceParam::instance().volumicOrSurfacicForce);
163
164
1/2
✓ Branch 2 taken 24 times.
✗ Branch 3 not taken.
 24 r_instance.hyperElasticLaw.hyperElasticLaws->InitParameters();
165
166
2/6
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 24 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 24 FEL_ASSERT(FelisceParam::instance().volumic_mass > 0.);
167
168 24 std::vector<PhysicalVariable> listVariable;
169 24 std::vector<std::size_t> listNumComp;
170
171
1/2
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
 24 listVariable.push_back(displacement);
172
1/2
✓ Branch 2 taken 24 times.
✗ Branch 3 not taken.
 24 listNumComp.push_back(this->dimension());
173
174
1/2
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
 24 m_listUnknown.push_back(displacement);
175
176
2/4
✓ Branch 2 taken 24 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 24 times.
 24 if (r_instance.hyperElasticLaw.pressureData->IsCompressible()) {
177 ✗ listVariable.push_back(pressure);
178 ✗ listNumComp.push_back(1);
179 ✗ m_listUnknown.push_back(pressure);
180 }
181
182
1/2
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
 24 definePhysicalVariable(listVariable,listNumComp);
183
1/2
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
 24 setNumberOfMatrix(MatrixNS::Nmatrices);
184
185
2/3
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 24 switch(r_instance.hyperElasticLaw.timeSchemeHyperelasticity) {
186 8 case none:
187 case midpoint:
188 8 break;
189 16 case half_sum:
190
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 setNumberOfVector(VectorNS::Nvector);
191 }
192 24 }
193
194 /***********************************************************************************/
195 /***********************************************************************************/
196
197 24 void LinearProblemHyperElasticity::InitializeDerivedAttributes()
198 {
199 // Retrieve the instance
200 24 const auto& r_instance = FelisceParam::instance();
201
202 using namespace HyperelasticityNS;
203
204 // Check here there are no Neumann or NeumannNormal boundary conditions in the case of the volumic force
205 // (it is not stricto sensu a limitation, so it could be lifted later, but at the moment it is a good way
206 // to make sure there is no mishap in the data parameter file, in which we could apply both volumic and
207 // surfacic force without being aware of it.)
208
1/2
✓ Branch 0 taken 24 times.
✗ Branch 1 not taken.
 24 if (m_typeOfForceApplied == Private::HyperElasticity::volumic) {
209
3/6
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 24 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 24 times.
 24 if (m_boundaryConditionList.numNeumannBoundaryCondition() || m_boundaryConditionList.numNeumannNormalBoundaryCondition()) {
210 ✗ FEL_ERROR("Neumann condition found along with a volumic force; it is not allowed at the moment. "
211 "Check there is no mishap in your data file.");
212 }
213 }
214
215 24 current_displacement_.duplicateFrom(this->solution());
216 24 current_velocity_.duplicateFrom(this->solution());
217 24 seq_current_displacement_.duplicateFrom(this->sequentialSolution());
218
219 24 current_velocity_.zeroEntries();
220 24 current_displacement_.zeroEntries();
221 24 seq_current_displacement_.zeroEntries();
222
223
2/3
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 24 switch(r_instance.hyperElasticLaw.timeSchemeHyperelasticity) {
224 8 case midpoint: {
225 8 midpoint_position_.duplicateFrom(this->solution());
226 8 seq_midpoint_position_.duplicateFrom(this->sequentialSolution());
227
228 8 midpoint_position_.zeroEntries();
229 8 seq_midpoint_position_.zeroEntries();
230 8 break;
231 }
232 16 case half_sum:
233 case none:
234 16 break;
235 }
236
237 24 m_snesContext.initialize(*this, MpiInfo::rankProc(), MpiInfo::numProc(), ApplyNaturalBoundaryConditions::no, FlagMatrixRHS::matrix_and_rhs);
238
239 24 snesInterface().setFunction(vector(), snesFunctionImpl, &m_snesContext);
240 24 snesInterface().setJacobian(matrix(0), matrix(0), snesJacobianImpl, &m_snesContext);
241 24 }
242
243 /***********************************************************************************/
244 /***********************************************************************************/
245
246 1044 void LinearProblemHyperElasticity::initPerElementType(ElementType eltType, FlagMatrixRHS flagMatrixRHS)
247 {
248 // Retrieve the instance
249 1044 const auto& r_instance = FelisceParam::instance();
250
251 1044 m_currentElementType = eltType;
252 1044 m_iDisplacement = m_listVariable.getVariableIdList(displacement);
253 1044 m_feDisplacement = m_listCurrentFiniteElement[m_iDisplacement];
254
255
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 1044 times.
 1044 if (r_instance.hyperElasticLaw.pressureData->IsCompressible()) {
256 ✗ r_instance.hyperElasticLaw.pressureData->Init(m_listVariable, m_listCurrentFiniteElement);
257
258 ✗ FEL_ASSERT(m_feDisplacement->numQuadraturePoint() == r_instance.hyperElasticLaw.pressureData->FiniteElementPressure()->numQuadraturePoint()
259 && "In its current state the code is unable to cope with different quadrature points for both variables.");
260 ✗ FEL_ASSERT(m_feDisplacement->numDof() == r_instance.hyperElasticLaw.pressureData->FiniteElementPressure()->numDof()
261 && "In its current state the code is unable to cope with different degres for both variables.");
262 }
263
264
2/4
✓ Branch 0 taken 1044 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1044 times.
✗ Branch 3 not taken.
 1044 if (flagMatrixRHS == FlagMatrixRHS::only_rhs || flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs)
265 1044 m_elemField.initialize(QUAD_POINT_FIELD, *m_feDisplacement, this->dimension());
266 1044 }
267
268 /***********************************************************************************/
269 /***********************************************************************************/
270
271 522000 void LinearProblemHyperElasticity::computeElementArray(
272 const std::vector<Point*>& elemPoint,
273 const std::vector<felInt>& elemIdPoint,
274 int& iel, FlagMatrixRHS flagMatrixRHS)
275 {
276 // Retrieve the instance
277
1/2
✓ Branch 1 taken 522000 times.
✗ Branch 2 not taken.
 522000 const auto& r_instance = FelisceParam::instance();
278
279 IGNORE_UNUSED_ELEM_ID_POINT;
280 using namespace HyperelasticityNS;
281
282
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 522000 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 522000 FEL_ASSERT(m_feDisplacement);
283 522000 CurrentFiniteElement& feDisplacement = *m_feDisplacement; // alias
284
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 522000 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 522000 FEL_ASSERT(flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs);
285
286
1/2
✓ Branch 1 taken 522000 times.
✗ Branch 2 not taken.
 522000 feDisplacement.updateFirstDeriv(0, elemPoint);
287
288
2/4
✓ Branch 2 taken 522000 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 522000 times.
 522000 if (r_instance.hyperElasticLaw.pressureData->IsCompressible())
289 ✗ r_instance.hyperElasticLaw.pressureData->FiniteElementPressure()->updateFirstDeriv(0, elemPoint);
290
291 522000 const int NquadraturePoint = feDisplacement.numQuadraturePoint();
292
293 522000 const std::size_t sizeMatrix = std::pow(this->dimension(), 2);
294
1/2
✓ Branch 1 taken 522000 times.
✗ Branch 2 not taken.
 522000 UBlasMatrix nonLinearPart(sizeMatrix, sizeMatrix);
295 522000 ElementVector& elementVector = *m_elementVector[0];
296
1/2
✓ Branch 1 taken 522000 times.
✗ Branch 2 not taken.
 522000 UBlasVector pressure;
297
298
2/3
✓ Branch 0 taken 42000 times.
✓ Branch 1 taken 480000 times.
✗ Branch 2 not taken.
 522000 switch (system_) {
299 42000 case HyperelasticityNS::static_: {
300
1/2
✓ Branch 3 taken 42000 times.
✗ Branch 4 not taken.
 42000 std::vector<UBlasVector> displacement_per_component = extractDisplacementPerComponent(iel, evaluationState(), seqEvaluationState());
301
302
2/4
✓ Branch 2 taken 42000 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 42000 times.
 42000 if (r_instance.hyperElasticLaw.pressureData->IsCompressible()) {
303 ✗ pressure = extractPressure(iel, evaluationState(), seqEvaluationState());
304 }
305
306 // DEV NOTE: nonLinearPart size could be defined once and for all. See when profiling whether it is
307 // a performance bottleneck.
308 42000 ElementMatrix& elementMatrix = *m_elementMat[0];
309
310 // At the moment all variables must share the same integration pattern, so no need to split calculation
311 // for displacement and pressure.
312
2/2
✓ Branch 0 taken 126000 times.
✓ Branch 1 taken 42000 times.
 168000 for(int quadraturePointIndex = 0; quadraturePointIndex < NquadraturePoint; ++quadraturePointIndex) {
313
1/2
✓ Branch 1 taken 126000 times.
✗ Branch 2 not taken.
 126000 UBlasMatrix linearPart;
314
1/2
✓ Branch 1 taken 126000 times.
✗ Branch 2 not taken.
 126000 UBlasVector rhsPart;
315
316 {
317 using namespace Private::HyperElasticity;
318
2/4
✓ Branch 1 taken 126000 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 126000 times.
✗ Branch 5 not taken.
 126000 calculateLinearAndNonLinearMatrixAndRhs(r_instance.hyperElasticLaw.hyperElasticLaws->GetHyperElasticLaw(this->m_currentLabel), feDisplacement,
319 displacement_per_component,
320 pressure,
321 quadraturePointIndex, linearPart,
322
1/2
✓ Branch 2 taken 126000 times.
✗ Branch 3 not taken.
 126000 nonLinearPart, rhsPart, r_instance.hyperElasticLaw.pressureData->IsCompressible());
323 }
324
325
2/4
✓ Branch 1 taken 126000 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 126000 times.
✗ Branch 5 not taken.
 126000 UBlasMatrix tangentMatrix = linearPart + nonLinearPart;
326
327
1/2
✓ Branch 1 taken 126000 times.
✗ Branch 2 not taken.
 126000 elementMatrix.grad_phi_i_GradientBasedHyperElastTensor_grad_phi_j(1., feDisplacement,
328 quadraturePointIndex,
329 tangentMatrix, 0, 0);
330
1/2
✓ Branch 1 taken 126000 times.
✗ Branch 2 not taken.
 126000 elementVector.GradientBasedHyperElasticityVector_grad_phi_i(1., feDisplacement,
331 quadraturePointIndex,
332 rhsPart, 0);
333 126000 }
334
335 // Set the source in the case of a volumic force. Surfacic ones are handled differently (through Neumann boundary condition)
336
1/2
✓ Branch 0 taken 42000 times.
✗ Branch 1 not taken.
 42000 if (m_typeOfForceApplied == Private::HyperElasticity::volumic) {
337 // Vector operators.
338
2/2
✓ Branch 1 taken 84000 times.
✓ Branch 2 taken 42000 times.
 126000 for (int iComp = 0; iComp < this->dimension(); ++iComp)
339
1/2
✓ Branch 1 taken 84000 times.
✗ Branch 2 not taken.
 84000 m_elemField.setValueVec(m_force, feDisplacement, iComp);
340
341
1/2
✓ Branch 2 taken 42000 times.
✗ Branch 3 not taken.
 42000 elementVector.source(1., feDisplacement, m_elemField, 0, this->dimension());
342 }
343
344 42000 break;
345 42000 }
346 480000 case HyperelasticityNS::dynamic_: {
347
348 // Note: the system matrix and rhs are NOT assembled here: they are calculated afterward within method
349 // computeSystemAfterAssembling() with the matrix calculated below and others calculated beforehand.
350
351 // Assemble the mass matrix in very first dynamic iteration. As it remains constant, it will be kept
352 // unchanged and uncleared throughout the rest of the program lifetime.
353
2/2
✓ Branch 0 taken 8000 times.
✓ Branch 1 taken 472000 times.
 480000 if (!is_first_dynamic_iteration_done_) {
354
1/4
✗ Branch 2 not taken.
✓ Branch 3 taken 8000 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 8000 FEL_ASSERT(m_elementMat[MatrixNS::mass_per_square_time]);
355 8000 ElementMatrix& elem_mass_matrix = *m_elementMat[MatrixNS::mass_per_square_time];
356
2/4
✓ Branch 1 taken 8000 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 8000 times.
✗ Branch 7 not taken.
 8000 elem_mass_matrix.phi_i_phi_j(2. * FelisceParam::instance().volumic_mass / std::pow(m_fstransient->timeStep, 2),
357 8000 feDisplacement, 0, 0, this->dimension());
358 }
359
360
361 // Assemble the new stiffness matrix.
362 480000 std::vector<UBlasVector> displacement_per_component;
363 480000 std::vector<UBlasVector> current_displacement_per_component;
364
365
2/4
✓ Branch 0 taken 240000 times.
✓ Branch 1 taken 240000 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 480000 switch(r_instance.hyperElasticLaw.timeSchemeHyperelasticity) {
366 240000 case midpoint:
367
1/2
✓ Branch 1 taken 240000 times.
✗ Branch 2 not taken.
 240000 displacement_per_component = extractDisplacementPerComponent(iel, midpoint_position_, seq_midpoint_position_);
368 240000 break;
369 240000 case half_sum: {
370
1/2
✓ Branch 3 taken 240000 times.
✗ Branch 4 not taken.
 240000 displacement_per_component = extractDisplacementPerComponent(iel, evaluationState(), seqEvaluationState());
371
372
2/2
✓ Branch 0 taken 4000 times.
✓ Branch 1 taken 236000 times.
 240000 if (!is_first_dynamic_iteration_done_)
373
1/2
✓ Branch 1 taken 4000 times.
✗ Branch 2 not taken.
 4000 current_displacement_per_component = extractDisplacementPerComponent(iel, current_displacement_, seq_current_displacement_);
374 240000 break;
375 }
376 ✗ case none:
377 ✗ break;
378 }
379
380
381
2/2
✓ Branch 0 taken 1440000 times.
✓ Branch 1 taken 480000 times.
 1920000 for (int quadraturePointIndex = 0; quadraturePointIndex < NquadraturePoint; ++quadraturePointIndex) {
382
1/2
✓ Branch 1 taken 1440000 times.
✗ Branch 2 not taken.
 1440000 UBlasMatrix linearPart;
383
1/2
✓ Branch 1 taken 1440000 times.
✗ Branch 2 not taken.
 1440000 UBlasVector rhsPart;
384
385 {
386 using namespace Private::HyperElasticity;
387
2/4
✓ Branch 1 taken 1440000 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1440000 times.
✗ Branch 5 not taken.
 1440000 calculateLinearAndNonLinearMatrixAndRhs(r_instance.hyperElasticLaw.hyperElasticLaws->GetHyperElasticLaw(this->m_currentLabel),
388 feDisplacement,
389 displacement_per_component,
390 pressure,
391 quadraturePointIndex,
392 linearPart,
393 nonLinearPart,
394 rhsPart,
395
1/2
✓ Branch 2 taken 1440000 times.
✗ Branch 3 not taken.
 1440000 r_instance.hyperElasticLaw.pressureData->IsCompressible()
396 );
397 }
398
399
2/4
✓ Branch 1 taken 1440000 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1440000 times.
✗ Branch 5 not taken.
 1440000 UBlasMatrix tangentMatrix = linearPart + nonLinearPart;
400
401
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1440000 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 1440000 FEL_ASSERT(flagMatrixRHS == FlagMatrixRHS::matrix_and_rhs);
402
403
2/4
✓ Branch 0 taken 720000 times.
✓ Branch 1 taken 720000 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 1440000 switch(r_instance.hyperElasticLaw.timeSchemeHyperelasticity) {
404 720000 case midpoint: {
405
1/2
✓ Branch 3 taken 720000 times.
✗ Branch 4 not taken.
 720000 m_elementMat[MatrixNS::new_stiffness]->grad_phi_i_GradientBasedHyperElastTensor_grad_phi_j(1., feDisplacement,
406 quadraturePointIndex,
407 tangentMatrix, 0, 0);
408
1/2
✓ Branch 1 taken 720000 times.
✗ Branch 2 not taken.
 720000 elementVector.GradientBasedHyperElasticityVector_grad_phi_i(1., feDisplacement,
409 quadraturePointIndex,
410 rhsPart, 0);
411 720000 break;
412 }
413 720000 case half_sum: {
414
1/2
✓ Branch 3 taken 720000 times.
✗ Branch 4 not taken.
 720000 m_elementMat[MatrixNS::new_stiffness]->grad_phi_i_GradientBasedHyperElastTensor_grad_phi_j(.5, feDisplacement,
415 quadraturePointIndex,
416 tangentMatrix, 0, 0);
417
1/2
✓ Branch 3 taken 720000 times.
✗ Branch 4 not taken.
 720000 m_elementVector[VectorNS::new_stiffness]->GradientBasedHyperElasticityVector_grad_phi_i(1., feDisplacement,
418 quadraturePointIndex,
419 rhsPart, 0);
420
421
422
423
2/2
✓ Branch 0 taken 12000 times.
✓ Branch 1 taken 708000 times.
 720000 if (!is_first_dynamic_iteration_done_) {
424 {
425 using namespace Private::HyperElasticity;
426
2/4
✓ Branch 1 taken 12000 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 12000 times.
✗ Branch 5 not taken.
 12000 calculateLinearAndNonLinearMatrixAndRhs(r_instance.hyperElasticLaw.hyperElasticLaws->GetHyperElasticLaw(this->m_currentLabel), feDisplacement,
427 current_displacement_per_component,
428 pressure,
429 quadraturePointIndex,
430 linearPart, nonLinearPart,
431 rhsPart,
432
1/2
✓ Branch 2 taken 12000 times.
✗ Branch 3 not taken.
 12000 r_instance.hyperElasticLaw.pressureData->IsCompressible()
433 );
434 }
435
436
2/4
✓ Branch 1 taken 12000 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 12000 times.
✗ Branch 5 not taken.
 12000 tangentMatrix = linearPart + nonLinearPart;
437
438
439
1/2
✓ Branch 3 taken 12000 times.
✗ Branch 4 not taken.
 12000 m_elementMat[MatrixNS::current_stiffness]->grad_phi_i_GradientBasedHyperElastTensor_grad_phi_j(.5, feDisplacement,
440 quadraturePointIndex,
441 tangentMatrix, 0, 0);
442
443
1/2
✓ Branch 3 taken 12000 times.
✗ Branch 4 not taken.
 12000 m_elementVector[VectorNS::current_stiffness]->GradientBasedHyperElasticityVector_grad_phi_i(1., feDisplacement,
444 quadraturePointIndex,
445 rhsPart, 0);
446 }
447 720000 break;
448 }
449 ✗ case none:
450 ✗ break;
451 }
452 1440000 }
453 480000 break;
454
455 480000 }
456
457 } // switch(system_)
458 522000 }
459
460 /***********************************************************************************/
461 /***********************************************************************************/
462
463 1044 void LinearProblemHyperElasticity::preSNESAssembleAdditional()
464 {
465 // Retrieve the instance
466 1044 const auto& r_instance = FelisceParam::instance();
467
468 using namespace HyperelasticityNS;
469
470 // clearMatrixRHS(); useless: already performed in preSNESAssemble()
471 1044 clearMatrix(MatrixNS::new_stiffness);
472
473
2/4
✓ Branch 0 taken 508 times.
✓ Branch 1 taken 536 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 1044 switch(r_instance.hyperElasticLaw.timeSchemeHyperelasticity) {
474 508 case midpoint: {
475
476 {
477 508 midpoint_position_.copyFrom(current_displacement_);
478 508 midpoint_position_.axpy(1., evaluationState());
479 508 midpoint_position_.scale(0.5);
480 }
481
482 {
483 508 gatherVector(midpoint_position_, seq_midpoint_position_);
484 }
485 508 break;
486 }
487 536 case half_sum: {
488 536 clearVector(VectorNS::new_stiffness);
489 536 break;
490 }
491 ✗ case none:
492 ✗ break;
493 }
494 1044 }
495
496 /***********************************************************************************/
497 /***********************************************************************************/
498
499 1044 void LinearProblemHyperElasticity::computeSystemAfterAssembling()
500 {
501 // Retrieve the instance
502 1044 const auto& r_instance = FelisceParam::instance();
503
504
2/2
✓ Branch 0 taken 960 times.
✓ Branch 1 taken 84 times.
 1044 if (system_ == HyperelasticityNS::dynamic_) {
505 using namespace HyperelasticityNS;
506
507 // In the very first dynamic iteration the mass matrices are assembled once and for all.
508
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 944 times.
 960 if (!is_first_dynamic_iteration_done_)
509 16 SetFirstDynamicIterationDone();
510
511 // Now all elements are settled to build the system.
512 {
513 // Tangent
514 960 matrix(MatrixNS::system).copyFrom(matrix(MatrixNS::mass_per_square_time), SAME_NONZERO_PATTERN);
515
516
2/4
✓ Branch 0 taken 480 times.
✓ Branch 1 taken 480 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 960 switch(r_instance.hyperElasticLaw.timeSchemeHyperelasticity) {
517 480 case midpoint: {
518 480 matrix(MatrixNS::system).axpy(0.5, matrix(MatrixNS::new_stiffness), SAME_NONZERO_PATTERN);
519 480 break;
520 }
521 480 case half_sum: {
522 480 matrix(MatrixNS::system).axpy(1., matrix(MatrixNS::new_stiffness), SAME_NONZERO_PATTERN);
523 480 matrix(MatrixNS::system).axpy(1., matrix(MatrixNS::current_stiffness), SAME_NONZERO_PATTERN);
524 480 break;
525 }
526 ✗ case none:
527 ✗ break;
528 }
529 }
530
531 {
532 // Residual
533
1/2
✓ Branch 1 taken 960 times.
✗ Branch 2 not taken.
 960 PetscVector vector;
534
1/2
✓ Branch 1 taken 960 times.
✗ Branch 2 not taken.
 960 vector.duplicateFrom(current_displacement_);
535 // < the contribution to RHS of current displacement term.
536
537
1/2
✓ Branch 2 taken 960 times.
✗ Branch 3 not taken.
 960 vector.copyFrom(evaluationState());
538
1/2
✓ Branch 1 taken 960 times.
✗ Branch 2 not taken.
 960 vector.axpy(-1., current_displacement_);
539
1/2
✓ Branch 2 taken 960 times.
✗ Branch 3 not taken.
 960 vector.axpy(-m_fstransient->timeStep, current_velocity_);
540
541
1/2
✓ Branch 1 taken 960 times.
✗ Branch 2 not taken.
 960 PetscVector dynamic_contribution_to_rhs;
542
1/2
✓ Branch 2 taken 960 times.
✗ Branch 3 not taken.
 960 dynamic_contribution_to_rhs.duplicateFrom(this->vector());
543
1/2
✓ Branch 2 taken 960 times.
✗ Branch 3 not taken.
 960 mult(matrix(MatrixNS::mass_per_square_time), vector, dynamic_contribution_to_rhs);
544
545
1/2
✓ Branch 2 taken 960 times.
✗ Branch 3 not taken.
 960 this->vector().axpy(1., dynamic_contribution_to_rhs);
546
547
548
2/4
✓ Branch 0 taken 480 times.
✓ Branch 1 taken 480 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 960 switch(r_instance.hyperElasticLaw.timeSchemeHyperelasticity) {
549 480 case midpoint:
550 480 break;
551 480 case half_sum: {
552
1/2
✓ Branch 3 taken 480 times.
✗ Branch 4 not taken.
 480 this->vector().axpy(1., this->vector(VectorNS::new_stiffness));
553
1/2
✓ Branch 3 taken 480 times.
✗ Branch 4 not taken.
 480 this->vector().axpy(1., this->vector(VectorNS::current_stiffness));
554 480 break;
555 }
556 ✗ case none:
557 ✗ break;
558 }
559 960 }
560 }
561 1044 }
562
563 /***********************************************************************************/
564 /***********************************************************************************/
565
566 526000 std::vector<UBlasVector> LinearProblemHyperElasticity::extractDisplacementPerComponent(
567 PetscInt iel,
568 PetscVector& mpi_vector,
569 PetscVector& seq_vector) const
570 {
571
1/2
✓ Branch 3 taken 526000 times.
✗ Branch 4 not taken.
 526000 std::vector<UBlasVector> ret(this->dimension());
572
573
2/2
✓ Branch 1 taken 1052000 times.
✓ Branch 2 taken 526000 times.
 1578000 for (int iComp = 0; iComp < this->dimension(); ++iComp) {
574
2/4
✓ Branch 1 taken 1052000 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 1052000 times.
✗ Branch 6 not taken.
 1052000 ret[iComp] = extractUnknownValue(m_iDisplacement, iComp, m_currentElementType, iel, mpi_vector, seq_vector);
575 }
576
577 526000 return ret;
578 }
579
580 /***********************************************************************************/
581 /***********************************************************************************/
582
583 ✗ inline UBlasVector LinearProblemHyperElasticity::extractPressure(
584 PetscInt iel,
585 PetscVector& mpi_vector,
586 PetscVector& seq_vector) const
587 {
588 // Retrieve the instance
589 ✗ const auto& r_instance = FelisceParam::instance();
590
591 ✗ FEL_ASSERT(r_instance.hyperElasticLaw.pressureData->IsCompressible());
592
593 ✗ return extractUnknownValue(r_instance.hyperElasticLaw.pressureData->IndexPressure(), 0, m_currentElementType, iel, mpi_vector, seq_vector);
594 }
595
596 /***********************************************************************************/
597 /***********************************************************************************/
598
599 1052000 UBlasVector LinearProblemHyperElasticity::extractUnknownValue(std::size_t index_unknown, std::size_t iComp, ElementType eltType,
600 std::size_t iel, PetscVector& mpi_vector, PetscVector& seq_vector) const
601 {
602 1052000 const SupportDofMesh& support = m_supportDofUnknown[index_unknown];
603
604 PetscInt id_element_support;
605
1/2
✓ Branch 1 taken 1052000 times.
✗ Branch 2 not taken.
 1052000 support.getIdElementSupport(eltType, iel, id_element_support);
606
607 1052000 const std::size_t Nlocal_support_dof = static_cast<std::size_t>(support.getNumSupportDof(id_element_support));
608
609
1/2
✓ Branch 1 taken 1052000 times.
✗ Branch 2 not taken.
 1052000 UBlasVector ret(Nlocal_support_dof);
610 PetscInt index_global_dof;
611
612
2/2
✓ Branch 0 taken 6312000 times.
✓ Branch 1 taken 1052000 times.
 7364000 for (std::size_t iSupport = 0; iSupport < Nlocal_support_dof; ++iSupport) {
613
1/2
✓ Branch 1 taken 6312000 times.
✗ Branch 2 not taken.
 6312000 m_dof.loc2glob(id_element_support, iSupport, index_unknown, iComp, index_global_dof);
614
615 PetscInt index_local_dof;
616
617 // m_mappingNodes is an empirical guess: it is the one that works not stupidly, ie returns a valid index if the dof is on current processor...
618 // DEV I'm not sure it will work with two variables, but the support that seems more logical failed utterly...
619
1/2
✓ Branch 1 taken 6312000 times.
✗ Branch 2 not taken.
 6312000 ISGlobalToLocalMappingApply(m_mappingNodes,
620 IS_GTOLM_MASK,
621 1, &index_global_dof,
622 FELISCE_PETSC_NULLPTR, &index_local_dof);
623
624 # ifndef NDEBUG
625 PetscInt localSize;
626
1/2
✓ Branch 1 taken 6312000 times.
✗ Branch 2 not taken.
 6312000 mpi_vector.getLocalSize(&localSize);
627
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6312000 times.
 6312000 if (index_local_dof >= localSize) {
628 int rank;
629 ✗ MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
630 ✗ std::cerr << '[' << rank << "] Dof index = " << index_local_dof
631 << " should be lower than local displacement size which "
632 ✗ "is " << localSize << std::endl;
633 ✗ FEL_ASSERT(false);
634 }
635 # endif // NDEBUG
636
637 // DEV \todo These indexes should be accessible locally with the use of ghost vertices...
638
2/2
✓ Branch 0 taken 103096 times.
✓ Branch 1 taken 6208904 times.
 6312000 if (index_local_dof == -1) {
639 // Case in which the values aren't stored on the current processor... I must fetch the global value!!!!
640 PetscScalar * seq_vector_array;
641
1/2
✓ Branch 1 taken 103096 times.
✗ Branch 2 not taken.
 103096 seq_vector.getArray(&seq_vector_array);
642
1/2
✓ Branch 1 taken 103096 times.
✗ Branch 2 not taken.
 103096 AOApplicationToPetsc(m_ao, 1, &index_global_dof);
643
1/2
✓ Branch 1 taken 103096 times.
✗ Branch 2 not taken.
 103096 ret(iSupport) = seq_vector_array[index_global_dof];
644
1/2
✓ Branch 1 taken 103096 times.
✗ Branch 2 not taken.
 103096 seq_vector.restoreArray(&seq_vector_array);
645 } else {
646 PetscScalar* mpi_vector_array;
647
1/2
✓ Branch 1 taken 6208904 times.
✗ Branch 2 not taken.
 6208904 mpi_vector.getArray(&mpi_vector_array);
648
1/2
✓ Branch 1 taken 6208904 times.
✗ Branch 2 not taken.
 6208904 ret(iSupport) = mpi_vector_array[index_local_dof];
649
1/2
✓ Branch 1 taken 6208904 times.
✗ Branch 2 not taken.
 6208904 mpi_vector.restoreArray(&mpi_vector_array);
650 }
651 }
652
653 2104000 return ret;
654 }
655
656 /***********************************************************************************/
657 /***********************************************************************************/
658
659 ✗ void LinearProblemHyperElasticity::computeElementArrayBoundaryCondition(
660 const std::vector<Point*>& elemPoint,
661 const std::vector<felInt>& elemIdPoint,
662 int& iel, FlagMatrixRHS flagMatrixRHS)
663 {
664 IGNORE_UNUSED_IEL;
665 IGNORE_UNUSED_ELEM_ID_POINT;
666
667 ✗ switch(flagMatrixRHS) {
668 ✗ case FlagMatrixRHS::only_rhs:
669 case FlagMatrixRHS::matrix_and_rhs: {
670 ✗ CurvilinearFiniteElement& feCurvDisplacement = *m_listCurvilinearFiniteElement[m_iDisplacement];
671 ✗ feCurvDisplacement.updateMeasNormal(0, elemPoint);
672
673 // Below: only when spatial function
674 ✗ if (m_boundaryConditionList.NeumannNormal(0)->typeValueBC() == FunctionS) {
675 {
676 ✗ const int component_index = 0;
677 ✗ m_elemFieldNeumannNormal[0].setValueVec(m_force, feCurvDisplacement, component_index);
678
679 // Important: in my first approach I added the following line, to mimic what was done in volumic force case:
680 //m_elemVecBD->source(1., feCurvDisplacement, m_elemFieldNeumannNormal[0], block, Ncomponent);
681 // It was wrong: the source was added more properly in applyNaturalBoundaryCondition()
682 }
683 }
684 ✗ break;
685 }
686 ✗ case FlagMatrixRHS::only_matrix:
687 ✗ break;
688 }
689 }
690
691 /***********************************************************************************/
692 /***********************************************************************************/
693
694 16 void LinearProblemHyperElasticity::GoToDynamic()
695 {
696 16 system_ = HyperelasticityNS::dynamic_;
697 16 UpdateDisplacementBetweenTimeStep();
698 16 }
699
700 /***********************************************************************************/
701 /***********************************************************************************/
702
703 320 void LinearProblemHyperElasticity::endIteration()
704 {
705 // Retrieve the instance
706 320 const auto& r_instance = FelisceParam::instance();
707
708 // First update velocities.
709 {
710
1/2
✓ Branch 1 taken 320 times.
✗ Branch 2 not taken.
 320 current_velocity_.scale(-1.);
711
1/2
✓ Branch 1 taken 320 times.
✗ Branch 2 not taken.
 320 PetscVector diff_displ;
712
1/2
✓ Branch 2 taken 320 times.
✗ Branch 3 not taken.
 320 diff_displ.duplicateFrom(this->solution());
713 {
714
1/2
✓ Branch 2 taken 320 times.
✗ Branch 3 not taken.
 320 diff_displ.copyFrom(this->solution());
715
1/2
✓ Branch 1 taken 320 times.
✗ Branch 2 not taken.
 320 diff_displ.axpy(-1., current_displacement_);
716
1/2
✓ Branch 2 taken 320 times.
✗ Branch 3 not taken.
 320 diff_displ.scale(2. / m_fstransient->timeStep);
717 }
718
1/2
✓ Branch 1 taken 320 times.
✗ Branch 2 not taken.
 320 current_velocity_.axpy(1., diff_displ);
719 320 }
720
721 // Then update displacement.
722 320 UpdateDisplacementBetweenTimeStep();
723
724 using namespace HyperelasticityNS;
725
726
2/4
✓ Branch 0 taken 160 times.
✓ Branch 1 taken 160 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 320 switch(r_instance.hyperElasticLaw.timeSchemeHyperelasticity) {
727 160 case midpoint:
728 160 break;
729 160 case half_sum:
730 160 matrix(MatrixNS::current_stiffness).copyFrom(matrix(MatrixNS::new_stiffness),SAME_NONZERO_PATTERN);
731 160 vector(VectorNS::current_stiffness).copyFrom(vector(VectorNS::new_stiffness));
732 160 break;
733 ✗ case none:
734 ✗ break;
735 }
736 320 }
737
738 /***********************************************************************************/
739 /***********************************************************************************/
740
741 16 inline void LinearProblemHyperElasticity::SetFirstDynamicIterationDone()
742 {
743
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 16 times.
 16 FEL_ASSERT(!is_first_dynamic_iteration_done_ && "Should only be called once!");
744 16 is_first_dynamic_iteration_done_ = true;
745 16 }
746
747 /***********************************************************************************/
748 /***********************************************************************************/
749
750 336 inline void LinearProblemHyperElasticity::UpdateDisplacementBetweenTimeStep()
751 {
752 336 current_displacement_.copyFrom(this->solution());
753 336 gatherSolution();
754 336 seq_current_displacement_.copyFrom(this->sequentialSolution());
755 336 }
756
757 }
758