GCC Code Coverage Report

Directory: ./
File: Solver/linearProblemElasticCurvedBeam.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 126 132 95.5%
Branches: 116 218 53.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:
13 //
14
15 // System includes
16
17 // External includes
18
19 // Project includes
20 #include "Solver/linearProblemElasticCurvedBeam.hpp"
21 #include "FiniteElement/elementVector.hpp"
22 #include "FiniteElement/elementMatrix.hpp"
23
24 namespace felisce
25 {
26 5 LinearProblemElasticCurvedBeam::LinearProblemElasticCurvedBeam():
27
17/34
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 5 times.
✗ Branch 6 not taken.
✓ Branch 11 taken 5 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 5 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 5 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 5 times.
✗ Branch 21 not taken.
✓ Branch 23 taken 5 times.
✗ Branch 24 not taken.
✓ Branch 26 taken 5 times.
✗ Branch 27 not taken.
✓ Branch 29 taken 5 times.
✗ Branch 30 not taken.
✓ Branch 32 taken 5 times.
✗ Branch 33 not taken.
✓ Branch 35 taken 5 times.
✗ Branch 36 not taken.
✓ Branch 38 taken 5 times.
✗ Branch 39 not taken.
✓ Branch 41 taken 5 times.
✗ Branch 42 not taken.
✓ Branch 44 taken 5 times.
✗ Branch 45 not taken.
✓ Branch 47 taken 5 times.
✗ Branch 48 not taken.
✓ Branch 50 taken 5 times.
✗ Branch 51 not taken.
✓ Branch 53 taken 5 times.
✗ Branch 54 not taken.
 5 LinearProblem("Elastic Curved Beam")
28 {
29 5 m_createNormalTangent = true;
30 5 }
31
32
33 10 LinearProblemElasticCurvedBeam::~LinearProblemElasticCurvedBeam()
34 {
35
36 10 Elasticity.clear();
37
38 }
39
40
41 720 void LinearProblemElasticCurvedBeam::userElementComputeBD(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint,felInt& iel)
42 {
43 (void) elemPoint;
44 (void) elemIdPoint;
45 (void) iel;
46 720 }
47
48
49 5 void LinearProblemElasticCurvedBeam::initialize(std::vector<GeometricMeshRegion::Pointer>& mesh, FelisceTransient::Pointer fstransient, MPI_Comm& comm, bool doUseSNES)
50 {
51
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 LinearProblem::initialize(mesh,comm, doUseSNES);
52 5 m_fstransient = fstransient;
53
54
1/2
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
 5 std::vector<PhysicalVariable> listVariable(2);
55
1/2
✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
 5 std::vector<std::size_t> listNumComp(2);
56 5 listVariable[0] = displacement;
57 5 listNumComp[0] = 2;
58 5 listVariable[1] = rotation;
59 5 listNumComp[1] = 1;
60
61 //define unknown of the linear system.
62
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 m_listUnknown.push_back(displacement);
63
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 m_listUnknown.push_back(rotation);
64
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 this->definePhysicalVariable(listVariable, listNumComp);
65
66 //compute normal & tangent on node of the mesh
67
1/2
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
 5 m_mesh[m_currentMesh]->computeNormalTangent();
68 5 m_mesh[m_currentMesh]->createNormalTangent() = m_createNormalTangent;
69 //allocate matrix to build elementary elastic matrix
70
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 Elasticity.resize(2,2);
71
4/8
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 5 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 5 times.
✗ Branch 11 not taken.
 5 m_B0.resize(2,6); m_B1.resize(2,6); m_B2.resize(2,6); m_B.resize(2,6);
72
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 m_phiGeo_reinterpol.resize(2);
73
2/4
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
 5 m_F0.resize(2,2); m_F1.resize(2,2);
74
2/4
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
 5 m_F0_reinterpol.resize(2,2); m_F1_reinterpol.resize(2,2);
75
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 m_reinterpolation = FelisceParam::instance().reinterpolation;
76
77 5 }
78
79 160 void LinearProblemElasticCurvedBeam::initPerElementTypeBD(ElementType eltType, FlagMatrixRHS flagMatrixRHS)
80 {
81 (void) eltType;
82 (void) flagMatrixRHS;
83
84 160 m_iDisp = m_listVariable.getVariableIdList(displacement);
85 160 m_iRota = m_listVariable.getVariableIdList(rotation);
86 160 m_feDisp = m_listCurvilinearFiniteElement[m_iDisp];
87 160 m_feRota = m_listCurvilinearFiniteElement[m_iRota];
88 160 m_numPoint = m_feDisp->numPoint();
89 160 m_numCoor = m_feDisp->numCoor();
90
1/2
✓ Branch 2 taken 160 times.
✗ Branch 3 not taken.
 320 m_matVel = new ElementMatrix(std::vector<const CurBaseFiniteElement*>{m_feRota},
91
1/2
✓ Branch 4 taken 160 times.
✗ Branch 5 not taken.
 320 std::vector<std::size_t>{m_listVariable[m_iRota].numComponent()},
92
3/6
✓ Branch 4 taken 160 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 160 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 160 times.
✗ Branch 11 not taken.
 160 std::vector<std::size_t>{m_listVariable[m_iRota].numComponent()});
93
94 // define elemField to put term u_n_1.phi_j in RHS.
95 160 m_elemField.initialize(DOF_FIELD,*m_feDisp,m_listVariable[m_iDisp].numComponent() );
96
97 160 }
98
99
100 ✗ void LinearProblemElasticCurvedBeam::computeElementArrayBD(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint, felInt& iel, FlagMatrixRHS flagMatrixRHS) {
101 (void) elemIdPoint;
102 (void) elemPoint;
103 (void) iel;
104 (void) flagMatrixRHS;
105 }
106
107
108 1560 void LinearProblemElasticCurvedBeam::computeElementArrayBD(const std::vector<Point*>& elemPoint, const std::vector<felInt>& elemIdPoint, const std::vector<Point*>& elemNormal, const std::vector< std::vector<Point*> >& elemTangent, felInt& iel, FlagMatrixRHS flagMatrixRHS)
109 {
110 (void) elemIdPoint;
111 (void) elemPoint;
112 (void) iel;
113 (void) flagMatrixRHS;
114
115 // parameters
116
1/2
✓ Branch 1 taken 1560 times.
✗ Branch 2 not taken.
 1560 double rho = FelisceParam::instance().densitySolid;
117 1560 double tau = m_fstransient->timeStep;
118
119
1/2
✓ Branch 1 taken 1560 times.
✗ Branch 2 not taken.
 1560 double young = FelisceParam::instance().young;
120
1/2
✓ Branch 1 taken 1560 times.
✗ Branch 2 not taken.
 1560 double nu = FelisceParam::instance().poisson;
121
1/2
✓ Branch 1 taken 1560 times.
✗ Branch 2 not taken.
 1560 double depth = FelisceParam::instance().beam.depth;
122
1/2
✓ Branch 1 taken 1560 times.
✗ Branch 2 not taken.
 1560 double thickness = FelisceParam::instance().thickness;
123 1560 double massVel = rho*thickness/(tau*tau);
124
125 // update finite element shape functions and its first derivatives at quadrature points
126
1/2
✓ Branch 1 taken 1560 times.
✗ Branch 2 not taken.
 1560 m_feDisp->updateBeam(0, elemPoint, elemNormal, elemTangent, thickness);
127
1/2
✓ Branch 1 taken 1560 times.
✗ Branch 2 not taken.
 1560 m_feRota->updateBeam(0, elemPoint, elemNormal, elemTangent, thickness);
128
129
130 //creation of quadrature point in z
131 1560 const DegreeOfExactness degOfExactness_z = DegreeOfExactness_2;
132
1/2
✓ Branch 1 taken 1560 times.
✗ Branch 2 not taken.
 1560 m_quadratureRule_z = &listQuadratureRuleSegment.quadratureRuleByExactness(degOfExactness_z);
133 1560 m_numQuadraturePoint_z = m_quadratureRule_z->numQuadraturePoint();
134 1560 m_numQuadraturePoint_r = m_feDisp->numQuadraturePoint();
135
1/2
✓ Branch 1 taken 1560 times.
✗ Branch 2 not taken.
 1560 m_weightMeas_z.resize(m_numQuadraturePoint_z);
136
2/2
✓ Branch 0 taken 3120 times.
✓ Branch 1 taken 1560 times.
 4680 for (int iz = 0; iz<m_numQuadraturePoint_z; iz++)
137
1/2
✓ Branch 1 taken 3120 times.
✗ Branch 2 not taken.
 3120 m_weightMeas_z[iz] = m_quadratureRule_z->weight(iz);
138
139 // elementary elastic matrix
140
2/2
✓ Branch 1 taken 3120 times.
✓ Branch 2 taken 1560 times.
 4680 for(int ir=0; ir<m_feDisp->numQuadraturePoint(); ir++) {
141 3120 m_Grr0 = 0.; m_Grr1 = 0.; m_Grr2 = 0.;
142
1/2
✓ Branch 2 taken 3120 times.
✗ Branch 3 not taken.
 3120 m_dPhiGeo = m_feDisp->dPhiGeo(ir);
143
1/2
✓ Branch 2 taken 3120 times.
✗ Branch 3 not taken.
 3120 m_phiGeo = m_feDisp->phiGeo(ir);
144 // convariant basis linear part
145
2/4
✓ Branch 1 taken 3120 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3120 times.
✗ Branch 5 not taken.
 3120 m_F0 = m_feDisp->F0(ir);
146
2/4
✓ Branch 1 taken 3120 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3120 times.
✗ Branch 5 not taken.
 3120 m_F1 = m_feDisp->F1(ir);
147 // convariant component G_rr
148
1/2
✓ Branch 1 taken 3120 times.
✗ Branch 2 not taken.
 3120 m_feDisp->compute_Grr0_Grr1_Grr2(m_F0, m_F1);
149 3120 m_Grr0 = m_feDisp->Grr0();
150 3120 m_Grr1 = m_feDisp->Grr1();
151 3120 m_Grr2 = m_feDisp->Grr2();
152 // polynomial reinterpolation tying points in treatment of numeral locking
153
1/2
✓ Branch 0 taken 3120 times.
✗ Branch 1 not taken.
 3120 if (m_reinterpolation) {
154
1/2
✓ Branch 1 taken 3120 times.
✗ Branch 2 not taken.
 3120 m_reinterpoledPhiGeo();
155
1/2
✓ Branch 1 taken 3120 times.
✗ Branch 2 not taken.
 3120 m_compute_F0_F1_reinterpoled(elemNormal, elemPoint, thickness);
156 }
157
1/2
✓ Branch 1 taken 3120 times.
✗ Branch 2 not taken.
 3120 m_computeStiffnessMatrix(elemTangent, m_reinterpolation);
158
159
2/2
✓ Branch 0 taken 6240 times.
✓ Branch 1 taken 3120 times.
 9360 for(int iz=0; iz<m_numQuadraturePoint_z; iz++) {
160
1/2
✓ Branch 1 taken 6240 times.
✗ Branch 2 not taken.
 6240 Elasticity.clear();
161
1/2
✓ Branch 1 taken 6240 times.
✗ Branch 2 not taken.
 6240 double z = m_quadratureRule_z->quadraturePointCoor(iz,0);
162 6240 m_Grr = m_Grr0+z*m_Grr1+z*z*m_Grr2;
163 6240 m_Girr = 1./m_Grr;
164
5/10
✓ Branch 1 taken 6240 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6240 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6240 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 6240 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 6240 times.
✗ Branch 14 not taken.
 6240 m_B = m_B0 + z*m_B1 + z*z*m_B2;
165
1/2
✓ Branch 1 taken 6240 times.
✗ Branch 2 not taken.
 6240 m_weightMeas = std::sqrt(0.25*thickness*thickness*m_Grr) * m_feDisp->weightMeas[ir] * m_weightMeas_z[iz];
166
1/2
✓ Branch 1 taken 6240 times.
✗ Branch 2 not taken.
 6240 Elasticity(0,0) = young/(1-nu*nu)*m_Girr*m_Girr;
167
1/2
✓ Branch 1 taken 6240 times.
✗ Branch 2 not taken.
 6240 Elasticity(1,1) = 2.*(young/(1+nu))/(thickness*thickness)*m_Girr;
168
169
1/2
✓ Branch 1 taken 6240 times.
✗ Branch 2 not taken.
 6240 UBlasMatrix temp_type (2,6);
170
3/6
✓ Branch 1 taken 6240 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6240 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6240 times.
✗ Branch 8 not taken.
 6240 temp_type= prod(trans(m_B),Elasticity);
171
4/8
✓ Branch 1 taken 6240 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6240 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6240 times.
✗ Branch 8 not taken.
✓ Branch 13 taken 6240 times.
✗ Branch 14 not taken.
 6240 m_elementMatBD[0]->mat() += prod(temp_type, m_B) * m_weightMeas * depth;
172 6240 }
173 }
174
175 // elementary velocity matrix
176
1/2
✓ Branch 1 taken 1560 times.
✗ Branch 2 not taken.
 1560 m_feDisp->updateMeasQuadPt(0, elemPoint);
177
1/2
✓ Branch 2 taken 1560 times.
✗ Branch 3 not taken.
 1560 m_matVel->mat().clear();
178
1/2
✓ Branch 1 taken 1560 times.
✗ Branch 2 not taken.
 1560 m_matVel->phi_i_phi_j(1.,*m_feDisp,0,0,1);
179
3/6
✓ Branch 2 taken 1560 times.
✗ Branch 3 not taken.
✓ Branch 7 taken 1560 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1560 times.
✗ Branch 11 not taken.
 1560 m_elementMatBD[0]->matBlock(0,0) += massVel*m_matVel->mat();
180
3/6
✓ Branch 2 taken 1560 times.
✗ Branch 3 not taken.
✓ Branch 7 taken 1560 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1560 times.
✗ Branch 11 not taken.
 1560 m_elementMatBD[0]->matBlock(1,1) += massVel*m_matVel->mat();
181
182
183 //rho*epsilon/tau^2 * M X^{n-1}
184
1/2
✓ Branch 3 taken 1560 times.
✗ Branch 4 not taken.
 1560 m_elemField.setValue(externalVec(0), *m_feDisp, iel, m_iDisp, m_ao, dof());
185
6/12
✓ Branch 1 taken 1560 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1560 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 1560 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 1560 times.
✗ Branch 12 not taken.
✓ Branch 16 taken 1560 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 1560 times.
✗ Branch 20 not taken.
 1560 m_elementVectorBD[0]->vecBlock(0) += prod(massVel*m_matVel->mat(), row(m_elemField.get_val(),0));
186
6/12
✓ Branch 1 taken 1560 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1560 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 1560 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 1560 times.
✗ Branch 12 not taken.
✓ Branch 16 taken 1560 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 1560 times.
✗ Branch 20 not taken.
 1560 m_elementVectorBD[0]->vecBlock(1) += prod(massVel*m_matVel->mat(), row(m_elemField.get_val(),1));
187 1560 }
188
189
190 // to compute in each quadrature point basis functions reinterpoled
191 3120 void LinearProblemElasticCurvedBeam::m_reinterpoledPhiGeo() {
192
2/2
✓ Branch 0 taken 6240 times.
✓ Branch 1 taken 3120 times.
 9360 for (int iNode = 0; iNode<m_numPoint; iNode++)
193 6240 m_phiGeo_reinterpol(iNode) = 0.5;
194 3120 }
195
196
197 //----------------------------------------------------------------------
198 // to compute in a quadrature F0 & F1 reinterpoled
199 3120 void LinearProblemElasticCurvedBeam::m_compute_F0_F1_reinterpoled(const std::vector<Point*>& elemNormal, const std::vector<Point*>& elemPoint, double thickness) {
200
201 3120 int iRefcoor=0;
202 double sum1, sum2, sum3;
203 3120 m_F0_reinterpol.clear(); m_F1_reinterpol.clear();
204
205
2/2
✓ Branch 0 taken 6240 times.
✓ Branch 1 taken 3120 times.
 9360 for (int jcoor=0; jcoor<m_numCoor; jcoor++){
206 6240 sum1= 0.; sum2= 0.; sum3= 0.;
207
2/2
✓ Branch 0 taken 12480 times.
✓ Branch 1 taken 6240 times.
 18720 for (int iNode=0; iNode<m_numPoint; iNode++){
208 12480 sum1 += m_dPhiGeo(iRefcoor,iNode)*elemPoint[iNode]->coor(jcoor);
209 12480 sum2 += 0.5*thickness*m_phiGeo_reinterpol(iNode)*elemNormal[iNode]->coor(jcoor);
210 12480 sum3 += 0.5*thickness*m_dPhiGeo(iRefcoor,iNode)*elemNormal[iNode]->coor(jcoor);
211 }
212 6240 m_F0_reinterpol(0,jcoor) = sum1; m_F0_reinterpol(1,jcoor) = sum2; m_F1_reinterpol(0,jcoor) = sum3;
213 }
214 3120 }
215
216 3120 void LinearProblemElasticCurvedBeam::m_computeStiffnessMatrix(const std::vector< std::vector<Point*> >& elemTangent, bool reinterpolation) {
217 3120 double thickness = FelisceParam::instance().thickness;
218
219 3120 m_B0.clear(); m_B1.clear(); m_B2.clear();
220
2/2
✓ Branch 0 taken 6240 times.
✓ Branch 1 taken 3120 times.
 9360 for (int iNode=0; iNode<m_numPoint; iNode++){
221 6240 m_B0(0,iNode) = m_dPhiGeo(0,iNode)*m_F0(0,0);
222 6240 m_B1(0,iNode) = m_dPhiGeo(0,iNode)*m_F1(0,0);
223 6240 m_B0(0,iNode+m_numPoint) = m_dPhiGeo(0,iNode)*m_F0(0,1);
224 6240 m_B1(0,iNode+m_numPoint) = m_dPhiGeo(0,iNode)*m_F1(0,1);
225 6240 m_B1(0,iNode+2*m_numPoint) = 0.5*thickness*m_dPhiGeo(0,iNode)*(m_F0(0,0)*elemTangent[0][iNode]->coor(0)+m_F0(0,1)*elemTangent[0][iNode]->coor(1));
226 6240 m_B2(0,iNode+2*m_numPoint) = 0.5*thickness*m_dPhiGeo(0,iNode)*(m_F1(0,0)*elemTangent[0][iNode]->coor(0)+m_F1(0,1)*elemTangent[0][iNode]->coor(1));
227
228 //if reinterpolation B rz
229
1/2
✓ Branch 0 taken 6240 times.
✗ Branch 1 not taken.
 6240 if (reinterpolation) {
230 6240 m_B0(1,iNode) = m_dPhiGeo(0,iNode)*m_F0_reinterpol(1,0);
231 6240 m_B0(1,iNode+m_numPoint) = m_dPhiGeo(0,iNode)*m_F0_reinterpol(1,1);
232 6240 m_B0(1,iNode+2*m_numPoint) = 0.5*thickness*m_phiGeo_reinterpol(iNode)*(m_F0_reinterpol(0,0)*elemTangent[0][iNode]->coor(0)+m_F0_reinterpol(0,1)*elemTangent[0][iNode]->coor(1));
233 12480 m_B1(1,iNode+2*m_numPoint) = 0.5*thickness*(m_dPhiGeo(0,iNode)*(m_F0_reinterpol(1,0)*elemTangent[0][iNode]->coor(0)+m_F0_reinterpol(1,1)*elemTangent[0][iNode]->coor(1))+
234 6240 m_phiGeo_reinterpol(iNode)*(m_F1_reinterpol(0,0)*elemTangent[0][iNode]->coor(0)+m_F1_reinterpol(0,1)*elemTangent[0][iNode]->coor(1)));
235 }
236 else {
237 ✗ m_B0(1,iNode) = m_dPhiGeo(0,iNode)*m_F0(1,0);
238 ✗ m_B0(1,iNode+m_numPoint) = m_dPhiGeo(0,iNode)*m_F0(1,1);
239 ✗ m_B0(1,iNode+2*m_numPoint) = 0.5*thickness*m_phiGeo(iNode)*(m_F0(0,0)*elemTangent[0][iNode]->coor(0)+m_F0(0,1)*elemTangent[0][iNode]->coor(1));
240 ✗ m_B1(1,iNode+2*m_numPoint) = 0.5*thickness*(m_dPhiGeo(0,iNode)*(m_F0(1,0)*elemTangent[0][iNode]->coor(0)+m_F0(1,1)*elemTangent[0][iNode]->coor(1))+
241 ✗ m_phiGeo(iNode)*(m_F1(0,0)*elemTangent[0][iNode]->coor(0)+m_F1(0,1)*elemTangent[0][iNode]->coor(1)));
242 }
243 }
244 3120 }
245
246 }
247
248