GCC Code Coverage Report

Directory: ./
File: Solver/linearProblemElasticity.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 201 235 85.5%
Branches: 154 330 46.7%
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 #include <numeric>
17 #include <cmath>
18
19 // External includes
20 #include "Core/NoThirdPartyWarning/Petsc/mat.hpp"
21
22 // Project includes
23 #include "linearProblemElasticity.hpp"
24 #include "Core/felisceTools.hpp"
25 #include "Core/felisceTransient.hpp"
26 #include "Geometry/geometricMeshRegion.hpp"
27 #include "FiniteElement/elementMatrix.hpp"
28 #include "FiniteElement/elementVector.hpp"
29 #include "DegreeOfFreedom/boundaryCondition.hpp"
30
31 namespace felisce
32 {
33
34 namespace { // anonymous
35
36 enum KinematicParameter2D {
37 planeStress = 1,
38 planeStrain = 2
39 };
40
41 /*!
42 * \brief Matrix defined so that strain = Gradient2Strain x gradient(displacement)
43 *
44 * No definition provided on purpose: this function is intended to be specialized
45 */
46 template<int DimensionT>
47 UBlasMatrix Gradient2Strain();
48
49
50 /*!
51 * \brief Engineering elasticity tensor
52 *
53 *
54 * \param[in] E Young modulus
55 * \param[in] nu Poisson coefficient
56 */
57 UBlasMatrix EngineeringElastTensor_2D(double E, double nu, int kinematic_param);
58
59 UBlasMatrix EngineeringElastTensor_3D(double E, double nu);
60
61 //! Initialize operator GradientBasedElastTensor
62 UBlasMatrix InitializeGradientBasedElastTensor(int dimension);
63
64 //! Assign the type of force applied from the integer read in the data file
65 ElasticityNS::ForceType AssignForceType(int data_parameter);
66
67
68 class ErrorInvalidVolumicMass : public std::exception {
69 public:
70 //! Constructor
71 ErrorInvalidVolumicMass();
72
73 //! Destructor.
74 ✗ ~ErrorInvalidVolumicMass() noexcept override = default;
75
76 //! Copy constructor, to avoid warning due to user-declared destructor.
77 ErrorInvalidVolumicMass(const ErrorInvalidVolumicMass&) = default;
78
79 //! what() overload
80 const char* what() const noexcept override;
81
82 private:
83
84 //! Message to be displayed by what()
85 std::string m_message;
86 };
87
88 } // namespace anonymous
89
90 36 LinearProblemElasticityDynamic::LinearProblemElasticityDynamic(SpatialForce force)
91 : LinearProblem("Dynamic elasticity equation"),
92 36 m_feDisplacement(nullptr),
93
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 36 m_force(force),
94 36 transient_parameters_(nullptr),
95
6/12
✓ Branch 2 taken 36 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 36 times.
✗ Branch 6 not taken.
✓ Branch 10 taken 36 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 36 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 36 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 36 times.
✗ Branch 20 not taken.
 72 system_(ElasticityNS::static_)
96 {
97
98 36 }
99
100 /***********************************************************************************/
101 /***********************************************************************************/
102
103 96 LinearProblemElasticityDynamic::~LinearProblemElasticityDynamic()
104 = default;
105
106 /***********************************************************************************/
107 /***********************************************************************************/
108
109 36 void LinearProblemElasticityDynamic::initialize(std::vector<GeometricMeshRegion::Pointer>& mesh, FelisceTransient::Pointer fstransient, MPI_Comm& comm, bool doUseSNES)
110 {
111
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 36 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 36 FEL_ASSERT(doUseSNES && "Should be true!!!");
112
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 36 LinearProblem::initialize(mesh, comm, doUseSNES);
113
114 36 transient_parameters_ = fstransient;
115
2/4
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
 36 m_typeOfForceApplied = AssignForceType(FelisceParam::instance().volumicOrSurfacicForce);
116
117
1/2
✓ Branch 2 taken 36 times.
✗ Branch 3 not taken.
 36 std::vector<PhysicalVariable> listVariable(1);
118
1/2
✓ Branch 2 taken 36 times.
✗ Branch 3 not taken.
 36 std::vector<std::size_t> listNumComp(1);
119 36 listVariable[0] = displacement;
120 36 listNumComp[0] = this->dimension();
121
122
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 36 m_listUnknown.push_back(displacement);
123
124
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 36 definePhysicalVariable(listVariable,listNumComp); // Must be called after m_listUnknown is initialized
125
126
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 36 setNumberOfMatrix(3);
127
128
2/4
✓ Branch 2 taken 36 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 36 times.
✗ Branch 7 not taken.
 36 m_gradientBasedElastTensor = std::move(InitializeGradientBasedElastTensor(this->dimension()));
129
130
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 36 const double volumic_mass = FelisceParam::instance().volumic_mass;
131
132
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 36 times.
 36 if (volumic_mass < 0.)
133 ✗ throw ErrorInvalidVolumicMass();
134
135
1/4
✗ Branch 1 not taken.
✓ Branch 2 taken 36 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 36 FEL_ASSERT(fstransient->timeStep > 1.e-10);
136 36 }
137
138 /***********************************************************************************/
139 /***********************************************************************************/
140
141 37 void LinearProblemElasticityDynamic::initPerElementType(ElementType eltType, FlagMatrixRHS flagMatrixRHS)
142 {
143 IGNORE_UNUSED_ELT_TYPE;
144 IGNORE_UNUSED_FLAG_MATRIX_RHS;
145 37 m_iDisplacement = m_listVariable.getVariableIdList(displacement);
146 37 m_feDisplacement = m_listCurrentFiniteElement[m_iDisplacement];
147
148 37 m_elemField.initialize(QUAD_POINT_FIELD, *m_feDisplacement, this->dimension());
149 37 }
150
151 /***********************************************************************************/
152 /***********************************************************************************/
153
154 36 void LinearProblemElasticityDynamic::InitializeDerivedAttributes()
155 {
156 // Check here there are no Neumann or NeumannNormal boundary conditions in the case of the volumic force
157 // (it is not stricto sensu a limitation, so it could be lifted later, but at the moment it is a good way
158 // to make sure there is no mishap in the data parameter file, in which we could apply both volumic and
159 // surfacic force without being aware of it.)
160
2/2
✓ Branch 0 taken 32 times.
✓ Branch 1 taken 4 times.
 36 if (m_typeOfForceApplied == ElasticityNS::volumic) {
161
3/6
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 32 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 32 times.
 32 if (m_boundaryConditionList.numNeumannBoundaryCondition() || m_boundaryConditionList.numNeumannNormalBoundaryCondition()) {
162 ✗ FEL_ERROR("Neumann condition found along with a volumic force; it is not allowed at the moment. "
163 "Check there is no mishap in your data file.");
164 }
165 }
166
167 36 current_displacement_.duplicateFrom(this->solution());
168 36 current_velocity_.duplicateFrom(this->solution());
169
170 36 current_velocity_.zeroEntries();
171 36 current_displacement_.zeroEntries();
172 36 }
173
174 /***********************************************************************************/
175 /***********************************************************************************/
176
177 4209 void LinearProblemElasticityDynamic::computeElementArray(const std::vector<Point*>& elemPoint,
178 const std::vector<felInt>& elemIdPoint,
179 int& iel, FlagMatrixRHS flagMatrixRHS) {
180 IGNORE_UNUSED_FLAG_MATRIX_RHS;
181 IGNORE_UNUSED_IEL;
182 IGNORE_UNUSED_ELEM_ID_POINT;
183
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4209 times.
 4209 FEL_ASSERT(m_feDisplacement);
184
185
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 4209 times.
 4209 FEL_ASSERT(this->numberOfMatrices() == 3);
186
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 4209 times.
 4209 FEL_ASSERT(static_cast<std::size_t>(this->numberOfMatrices()) == m_elementMat.size());
187
188 4209 CurrentFiniteElement& feDisplacement = *m_feDisplacement; // alias
189
190 4209 feDisplacement.updateFirstDeriv(0, elemPoint);
191
192
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 4209 times.
 4209 FEL_ASSERT(transient_parameters_);
193 4209 const FelisceTransient& transient_param = *transient_parameters_;
194
195 // First iteration is static case.
196
2/3
✓ Branch 0 taken 4149 times.
✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
 4209 switch(system_) {
197 4149 case ElasticityNS::static_: {
198
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4149 times.
 4149 FEL_ASSERT(transient_param.iteration >= 0);
199
1/2
✓ Branch 0 taken 4149 times.
✗ Branch 1 not taken.
 4149 if (transient_param.iteration == 0) {
200 // Matrix operators.
201 4149 m_elementMat[0]->grad_phi_i_GradientBasedElastTensor_grad_phi_j(1.0, feDisplacement, m_gradientBasedElastTensor, 0, 0);
202
203 // Set the source in the case of a volumic force. Surfacic ones are handled differently (through Neumann boundary condition)
204
2/2
✓ Branch 0 taken 2149 times.
✓ Branch 1 taken 2000 times.
 4149 if (m_typeOfForceApplied == ElasticityNS::volumic) {
205 // Vector operators.
206
2/2
✓ Branch 1 taken 4438 times.
✓ Branch 2 taken 2149 times.
 6587 for (int iComp = 0; iComp < this->dimension(); ++iComp)
207 4438 m_elemField.setValueVec(m_force, feDisplacement, iComp);
208
209 2149 m_elementVector[0]->source(1., feDisplacement, m_elemField, 0, this->dimension());
210 }
211 }
212 4149 break;
213 }
214 60 case ElasticityNS::dynamic_: {
215 // Matrix operators.
216
1/4
✗ Branch 2 not taken.
✓ Branch 3 taken 60 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 60 FEL_ASSERT(m_elementMat[0]);
217 60 ElementMatrix& complete_matrix = *m_elementMat[0];
218
219
1/2
✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
 60 ElementMatrix mass_matrix;
220
1/2
✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
 60 ElementMatrix half_stiffness_matrix; // stiffness matrix multiplied by 0.5
221
222 // Calculate mass and stiffness element matrix
223 {
224
225
1/2
✓ Branch 3 taken 60 times.
✗ Branch 4 not taken.
 60 mass_matrix.duplicateFrom(*m_elementMat[0]);
226
1/2
✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
 60 half_stiffness_matrix.duplicateFrom(mass_matrix);
227
228
1/2
✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
 60 half_stiffness_matrix.grad_phi_i_GradientBasedElastTensor_grad_phi_j(0.5, feDisplacement, m_gradientBasedElastTensor, 0, 0);
229
1/2
✓ Branch 2 taken 60 times.
✗ Branch 3 not taken.
 60 mass_matrix.phi_i_phi_j(1., feDisplacement, 0, 0, this->dimension());
230 }
231
232 // Calculate coefficients used in Newmark scheme.
233
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 60 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 60 FEL_ASSERT(transient_param.timeStep > 1.e-10);
234 60 const double inv_time_step = 1. / transient_param.timeStep;
235
1/2
✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
 60 const double velocity_mass_coefficient = 2. * FelisceParam::instance().volumic_mass * inv_time_step;
236 60 const double displacement_mass_coefficient = velocity_mass_coefficient * inv_time_step;
237
238 // Build the system matrix
239
1/2
✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
 60 complete_matrix = mass_matrix;
240
1/2
✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
 60 complete_matrix *= displacement_mass_coefficient;
241
1/2
✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
 60 complete_matrix += half_stiffness_matrix;
242
243 // Now we must take into account the terms related to previous iteration.
244 // We must for that fill two other matrixes: one for the displacement from previous iteration and the other
245 // for the velocity. Once assembled, these matrixes will be multiplied by the vectors from previous iteration
246 // and the RHS will hence be evaluated.
247 60 ElementMatrix& current_displacement_matrix = *m_elementMat[ElasticityNS::current_displacement];
248
249
1/2
✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
 60 current_displacement_matrix = mass_matrix;
250
1/2
✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
 60 current_displacement_matrix *= displacement_mass_coefficient;
251
1/2
✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
 60 current_displacement_matrix -= half_stiffness_matrix;
252
253 60 ElementMatrix& current_velocity_matrix = *m_elementMat[ElasticityNS::current_velocity];
254
1/2
✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
 60 current_velocity_matrix = mass_matrix;
255
1/2
✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
 60 current_velocity_matrix *= velocity_mass_coefficient;
256 60 break;
257 60 }
258 } // switch
259 4209 }
260
261 /***********************************************************************************/
262 /***********************************************************************************/
263
264 120 void LinearProblemElasticityDynamic::ComputeRHS()
265 {
266
1/2
✓ Branch 1 taken 120 times.
✗ Branch 2 not taken.
 120 PetscVector buf;
267
1/2
✓ Branch 1 taken 120 times.
✗ Branch 2 not taken.
 120 buf.duplicateFrom(current_displacement_);
268
269
1/2
✓ Branch 2 taken 120 times.
✗ Branch 3 not taken.
 120 mult(matrix(ElasticityNS::current_displacement), current_displacement_, buf);
270
1/2
✓ Branch 3 taken 120 times.
✗ Branch 4 not taken.
 120 multAdd(matrix(ElasticityNS::current_velocity), current_velocity_, buf, vector());
271 120 }
272
273 /***********************************************************************************/
274 /***********************************************************************************/
275
276 ✗ void LinearProblemElasticityDynamic::computeElementArrayBoundaryCondition(const std::vector<Point*>& elemPoint,
277 const std::vector<felInt>& elemIdPoint,
278 int& iel, FlagMatrixRHS flagMatrixRHS)
279 {
280 IGNORE_UNUSED_FLAG_MATRIX_RHS;
281 IGNORE_UNUSED_IEL;
282 IGNORE_UNUSED_ELEM_ID_POINT;
283
284 ✗ CurvilinearFiniteElement& feCurvDisplacement = *m_listCurvilinearFiniteElement[m_iDisplacement];
285 ✗ feCurvDisplacement.updateMeasNormal(0, elemPoint);
286
287 // Below: only when spatial function
288 ✗ const BoundaryCondition& neumann_condition = *m_boundaryConditionList.Neumann(0);
289
290 ✗ if (neumann_condition.typeValueBC() == FunctionS) {
291 ✗ const auto& comp_list = neumann_condition.getComp();
292
293 ✗ for (auto component_index : comp_list) { // more efficient.
294 ✗ m_elemFieldNeumann[0].setValueVec(m_force, feCurvDisplacement, component_index);
295 }
296 // Important: in my first approach I added the following line, to mimic what was done in volumic force case:
297 //m_elemVecBD->source(1., feCurvDisplacement, m_elemFieldNeumannNormal[0], block, Ncomponent);
298 // It was wrong: the source was added more properly in applyNaturalBoundaryCondition()
299 }
300 }
301
302 /***********************************************************************************/
303 /***********************************************************************************/
304
305 120 void LinearProblemElasticityDynamic::endIteration()
306 {
307 PetscInt size;
308
1/2
✓ Branch 1 taken 120 times.
✗ Branch 2 not taken.
 120 current_displacement_.getLocalSize(&size);
309
310
1/2
✓ Branch 1 taken 120 times.
✗ Branch 2 not taken.
 120 PetscInt NdofLocalPerUnknown = numDofLocalPerUnknown(displacement);
311
312
1/2
✓ Branch 2 taken 120 times.
✗ Branch 3 not taken.
 120 std::vector<PetscInt> local_indexing(NdofLocalPerUnknown); // std::vector that contains integers from 0 to local size - 1
313
314 // std::iota(local_indexing.begin(), local_indexing.end(), 0);
315
316
2/2
✓ Branch 0 taken 17010 times.
✓ Branch 1 taken 120 times.
 17130 for (int i = 0; i < NdofLocalPerUnknown; ++i) // substitute that also works in C++ 98 (arguably more easy to understand anyway)
317 17010 local_indexing[i] = i;
318
319 120 double vel_min = 1.e20;
320 120 double vel_max = -1.e20;
321 120 double disp_current_min = 1.e20;
322 120 double disp_current_max = -1.e20;
323 120 double disp_new_min = 1.e20;
324 120 double disp_new_max = -1.e20;
325
326 {
327 PetscScalar * solution;
328
1/2
✓ Branch 2 taken 120 times.
✗ Branch 3 not taken.
 120 this->solution().getArray(&solution);
329
330 PetscScalar * displacement_prev;
331
1/2
✓ Branch 1 taken 120 times.
✗ Branch 2 not taken.
 120 current_displacement_.getArray(&displacement_prev);
332
333 PetscScalar * my_velocity;
334
1/2
✓ Branch 1 taken 120 times.
✗ Branch 2 not taken.
 120 current_velocity_.getArray(&my_velocity);
335
336
2/2
✓ Branch 0 taken 17010 times.
✓ Branch 1 taken 120 times.
 17130 for (PetscInt i = 0; i < size; ++i) {
337 17010 my_velocity[i] *= -1.;
338 17010 my_velocity[i] += 2. / transient_parameters_->timeStep * (solution[i] - displacement_prev[i]);
339 17010 vel_min = std::min(vel_min, my_velocity[i]);
340 17010 vel_max = std::max(vel_max, my_velocity[i]);
341 17010 disp_current_min = std::min(disp_current_min, displacement_prev[i]);
342 17010 disp_current_max = std::max(disp_current_max, displacement_prev[i]);
343 }
344
345
1/2
✓ Branch 2 taken 120 times.
✗ Branch 3 not taken.
 120 this->solution().restoreArray(&solution);
346
1/2
✓ Branch 1 taken 120 times.
✗ Branch 2 not taken.
 120 current_displacement_.restoreArray(&displacement_prev);
347
1/2
✓ Branch 1 taken 120 times.
✗ Branch 2 not taken.
 120 current_velocity_.restoreArray(&my_velocity);
348 }
349
350
5/10
✓ Branch 1 taken 120 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 120 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 120 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 120 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 120 times.
✗ Branch 14 not taken.
 120 std::cout << "Displ prev [min, max] = [" << disp_current_min << ", " << disp_current_max << "]\n";
351
5/10
✓ Branch 1 taken 120 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 120 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 120 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 120 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 120 times.
✗ Branch 14 not taken.
 120 std::cout << "Velocities [min, max] = [" << vel_min << ", " << vel_max << "]\n";
352
353 // Make sure to do that after the velocity calculation!
354
1/2
✓ Branch 1 taken 120 times.
✗ Branch 2 not taken.
 120 UpdateCurrentDisplacement();
355
356 {
357 PetscScalar * displacement_new;
358
1/2
✓ Branch 1 taken 120 times.
✗ Branch 2 not taken.
 120 current_displacement_.getArray(&displacement_new);
359
360
2/2
✓ Branch 0 taken 17010 times.
✓ Branch 1 taken 120 times.
 17130 for (PetscInt i = 0u; i < size; ++i) {
361 17010 disp_new_min = std::min(disp_new_min, displacement_new[i]);
362 17010 disp_new_max = std::max(disp_new_max, displacement_new[i]);
363 }
364
365
1/2
✓ Branch 1 taken 120 times.
✗ Branch 2 not taken.
 120 current_displacement_.restoreArray(&displacement_new);
366
367 }
368
369
5/10
✓ Branch 1 taken 120 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 120 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 120 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 120 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 120 times.
✗ Branch 14 not taken.
 120 std::cout << "Displ new [min, max] = [" << disp_new_min << ", " << disp_new_max << "]\n";
370 120 }
371
372 /***********************************************************************************/
373 /***********************************************************************************/
374
375 156 void LinearProblemElasticityDynamic::UpdateCurrentDisplacement() {
376 156 current_displacement_.copyFrom(this->solution());
377 156 }
378
379 /***********************************************************************************/
380 /***********************************************************************************/
381
382 4 void LinearProblemElasticityDynamic::GoToDynamic() {
383 4 system_ = ElasticityNS::dynamic_;
384 4 }
385
386 /***********************************************************************************/
387 /***********************************************************************************/
388
389 namespace { // anonymous
390
391 template<>
392 20 UBlasMatrix Gradient2Strain<3>() {
393 20 UBlasMatrix ret(6, 9);
394
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 ret.clear();
395
396
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 ret(0, 0) = 1.;
397
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 ret(1, 4) = 1.;
398
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 ret(2, 8) = 1.;
399
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 ret(3, 1) = 1.;
400
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 ret(3, 3) = 1.;
401
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 ret(4, 5) = 1.;
402
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 ret(4, 7) = 1.;
403
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 ret(5, 2) = 1.;
404
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 ret(5, 6) = 1.;
405
406 20 return ret;
407 }
408
409 /***********************************************************************************/
410 /***********************************************************************************/
411
412 template<>
413 16 UBlasMatrix Gradient2Strain<2>() {
414 16 UBlasMatrix ret(3, 4);
415
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 ret.clear();
416
417
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 ret(0, 0) = 1.;
418
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 ret(1, 3) = 1.;
419
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 ret(2, 1) = 1.;
420
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 ret(2, 2) = 1.;
421
422 16 return ret;
423 }
424
425 /***********************************************************************************/
426 /***********************************************************************************/
427
428 20 UBlasMatrix EngineeringElastTensor_3D(const double E, const double nu)
429 {
430
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 UBlasMatrix ret(6, 6);
431
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 ret.clear();
432
433
3/6
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 20 times.
✗ Branch 8 not taken.
 20 ret(0, 0) = ret(1, 1) = ret(2, 2) = 1.;
434
435 20 const double inv = 1. / (1. - nu);
436
3/6
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 20 times.
✗ Branch 8 not taken.
 20 ret(3, 3) = ret(4, 4) = ret(5, 5) = 0.5 * (1. - 2. * nu) * inv;
437
438 {
439 20 const double non_diag_value = nu * inv;
440
441
2/2
✓ Branch 0 taken 60 times.
✓ Branch 1 taken 20 times.
 80 for (std::size_t i = 0u; i < 3u; ++i) {
442
2/2
✓ Branch 0 taken 60 times.
✓ Branch 1 taken 60 times.
 120 for (std::size_t j = i + 1u; j < 3u; ++j)
443
2/4
✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 60 times.
✗ Branch 5 not taken.
 60 ret(i, j) = ret(j, i) = non_diag_value;
444 }
445 }
446
447 20 const double factor = E * (1. - nu) / ((1. + nu) * (1. - 2. * nu));
448
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 20 ret *= factor;
449
450 40 return ret;
451 }
452
453 /***********************************************************************************/
454 /***********************************************************************************/
455
456 16 UBlasMatrix EngineeringElastTensor_2D(const double E, const double nu, const int kinematic_param)
457 {
458 16 UBlasMatrix ret(3, 3);
459
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 ret.clear();
460
461
1/3
✗ Branch 0 not taken.
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 switch(kinematic_param) {
462 ✗ case planeStress: {
463 ✗ ret(0, 0) = 1.;
464 ✗ ret(1, 1) = 1.;
465 ✗ ret(2, 2) = 0.5 * (1. - nu);
466 ✗ ret(0, 1) = nu;
467 ✗ ret(1, 0) = nu;
468
469 ✗ const double factor = E / (1. - std::pow(nu, 2));
470 ✗ ret *= factor;
471
472 ✗ return ret;
473 }
474 16 case planeStrain: {
475 16 const double inv = 1. / (1. - nu);
476
477
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 ret(0, 0) = 1.;
478
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 ret(1, 1) = 1.;
479
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 ret(2, 2) = 0.5 * inv * (1. - 2. * nu);
480
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 ret(0, 1) = nu * inv;
481
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 ret(1, 0) = nu * inv;
482
483 16 const double factor = E * (1. - nu) / ((1. + nu) * (1. - 2 * nu));
484
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 ret *= factor;
485
486 16 return ret;
487 }
488 }
489
490 ✗ FEL_ERROR("SHOULD NEVER BE REACHED!!!");
491 ✗ return ret; // to avoid warning
492 }
493
494 /***********************************************************************************/
495 /***********************************************************************************/
496
497 36 UBlasMatrix InitializeGradientBasedElastTensor(const int dimension)
498 {
499
2/4
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 36 times.
✗ Branch 5 not taken.
 36 UBlasMatrix engineering_elast_tensor, gradient_2_strain;
500
501 // Use Poisson coefficient and Young modulus read in the input parameter file
502
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 36 const double nu = FelisceParam::instance().poisson;
503
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 36 const double E = FelisceParam::instance().young;
504
1/2
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
 36 const int kinematic_parameter = FelisceParam::instance().planeStressStrain;
505
506
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 36 if (kinematic_parameter != planeStrain && kinematic_parameter != planeStress)
507 ✗ FEL_ERROR("Must be 1 (plane stress) or 2 (plane strain)");
508
509
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 36 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 36 FEL_ASSERT(std::fabs(nu -1.) > 1.e-10);
510
511
2/3
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 36 switch(dimension) {
512 16 case 2: {
513
2/4
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 16 times.
✗ Branch 6 not taken.
 16 engineering_elast_tensor = std::move(EngineeringElastTensor_2D(E, nu, kinematic_parameter));
514
2/4
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 16 times.
✗ Branch 6 not taken.
 16 gradient_2_strain = std::move(Gradient2Strain<2>());
515 16 break;
516 }
517 20 case 3: {
518
2/4
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 20 times.
✗ Branch 6 not taken.
 20 engineering_elast_tensor = std::move(EngineeringElastTensor_3D(E, nu));
519
2/4
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 20 times.
✗ Branch 6 not taken.
 20 gradient_2_strain = std::move(Gradient2Strain<3>());
520 20 break;
521 }
522 ✗ default: {
523 ✗ FEL_ERROR("Dimension must be 2 or 3");
524 ✗ exit(1);
525 }
526 }
527
528 // Intermediate step as Boost ublas gets a lousy syntax...
529
2/4
✓ Branch 1 taken 36 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 36 times.
✗ Branch 6 not taken.
 36 UBlasMatrix product1 = std::move(boost::numeric::ublas::prod(engineering_elast_tensor, gradient_2_strain));
530
531
3/6
✓ Branch 2 taken 36 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 36 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 36 times.
✗ Branch 10 not taken.
 108 return std::move(prod(boost::numeric::ublas::trans(gradient_2_strain), std::move(product1)));
532 36 }
533
534 /***********************************************************************************/
535 /***********************************************************************************/
536
537 36 ElasticityNS::ForceType AssignForceType(const int data_parameter) {
538
539 // Check whether a surfacic or volumic force has been chosen
540
2/3
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 36 switch (data_parameter) {
541 4 case static_cast<int>(ElasticityNS::surfacic):
542 4 return ElasticityNS::surfacic;
543 32 case static_cast<int>(ElasticityNS::volumic):
544 32 return ElasticityNS::volumic;
545 ✗ default:
546 ✗ FEL_ERROR("Force applied should be either volumic or surfacic; make sure [solid] volumicOrSurfacicForce "
547 "is correctly filled in your data file and encompass an acceptable value");
548 }
549
550 ✗ FEL_ERROR("Should never be reached; if so bug to correct in the function!!!");
551 ✗ return ElasticityNS::volumic; // to avoid warning about missing warning
552 }
553
554 /***********************************************************************************/
555 /***********************************************************************************/
556
557 ✗ ErrorInvalidVolumicMass::ErrorInvalidVolumicMass() {
558 ✗ m_message = std::move("No volumic mass read in input parameter file.");
559 }
560
561 /***********************************************************************************/
562 /***********************************************************************************/
563
564 ✗ const char* ErrorInvalidVolumicMass::what() const noexcept {
565 ✗ return m_message.c_str();
566 }
567
568 } // namespace anonymous
569
570 } // namespace felisce
571