GCC Code Coverage Report

Directory: ./
File: FiniteElement/curvilinearFiniteElement.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 308 505 61.0%
Branches: 144 408 35.3%
Line Branch Exec Source
1 // ______ ______ _ _ _____ ______
2 // | ____| ____| | (_)/ ____| | ____|
3 // | |__ | |__ | | _| (___ ___| |__
4 // | __| | __| | | | |\___ \ / __| __|
5 // | | | |____| |____| |____) | (__| |____
6 // |_| |______|______|_|_____/ \___|______|
7 // Finite Elements for Life Sciences and Engineering
8 //
9 // License: LGL2.1 License
10 // FELiScE default license: LICENSE in root folder
11 //
12 // Main authors: J-F. Gerbeau
13 //
14
15 // System includes
16 #include <iomanip>
17 #include <array>
18
19 // External includes
20
21 // Project includes
22 #include "Core/felisceParam.hpp"
23 #include "Core/felisceTools.hpp"
24 #include "FiniteElement/curvilinearFiniteElement.hpp"
25 #include "Tools/math_utilities.hpp"
26
27 namespace felisce
28 {
29 117157 CurvilinearFiniteElement::CurvilinearFiniteElement(
30 const RefElement& refEle,
31 const GeoElement& geoEle,
32 const DegreeOfExactness& degOfExactness
33
1/2
✓ Branch 12 taken 117157 times.
✗ Branch 13 not taken.
 117157 ): CurBaseFiniteElement(refEle, geoEle, degOfExactness)
34 {
35 // Set the number of coordinates
36
1/2
✓ Branch 1 taken 117157 times.
✗ Branch 2 not taken.
 117157 const int numCoor = FelisceParam::instance().numCoor;
37
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 117149 times.
 117157 if (numCoor == 0) {
38 8 m_numCoor = m_geoEle.numCoor()+1;
39 } else { // We set manually. TODO: Implement line in 3D
40 117149 m_numCoor = numCoor;
41 }
42
43
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 117157 times.
 117157 FEL_ASSERT( m_numCoor <= 3 ); // curvilinear element not devised for 3D yet. VM
44
45
1/2
✓ Branch 2 taken 117157 times.
✗ Branch 3 not taken.
 117157 m_covBasis = std::vector<UBlasMatrix>(m_numQuadraturePoint);
46
1/2
✓ Branch 1 taken 117157 times.
✗ Branch 2 not taken.
 117157 m_F0.resize(m_numQuadraturePoint);
47
1/2
✓ Branch 1 taken 117157 times.
✗ Branch 2 not taken.
 117157 m_F1.resize(m_numQuadraturePoint);
48
1/2
✓ Branch 2 taken 117157 times.
✗ Branch 3 not taken.
 117157 m_covCompleteBasis = std::vector<UBlasMatrix>(m_numQuadraturePoint);
49
1/2
✓ Branch 2 taken 117157 times.
✗ Branch 3 not taken.
 117157 contravariantCompleteBasis = std::vector<UBlasMatrix>(m_numQuadraturePoint);
50
1/2
✓ Branch 2 taken 117157 times.
✗ Branch 3 not taken.
 117157 tangent = std::vector<UBlasMatrix>(m_numQuadraturePoint);
51
1/2
✓ Branch 2 taken 117157 times.
✗ Branch 3 not taken.
 117157 normal = std::vector<UBlasVector>(m_numQuadraturePoint);
52
1/2
✓ Branch 2 taken 117157 times.
✗ Branch 3 not taken.
 117157 covMetric = std::vector<UBlasMatrix>(m_numQuadraturePoint);
53
1/2
✓ Branch 2 taken 117157 times.
✗ Branch 3 not taken.
 117157 contravariantMetric = std::vector<UBlasMatrix>(m_numQuadraturePoint);
54
55
1/2
✓ Branch 1 taken 117157 times.
✗ Branch 2 not taken.
 117157 m_point.resize(m_numPoint,m_numCoor);
56
1/2
✓ Branch 1 taken 117157 times.
✗ Branch 2 not taken.
 117157 m_tangentVertice.resize(m_numCoor - 1);
57
2/2
✓ Branch 0 taken 165728 times.
✓ Branch 1 taken 117157 times.
 282885 for(int i=0 ; i<m_numCoor -1 ; i++)
58
1/2
✓ Branch 2 taken 165728 times.
✗ Branch 3 not taken.
 165728 m_tangentVertice[i].resize(m_numPoint,m_numCoor);
59
1/2
✓ Branch 1 taken 117157 times.
✗ Branch 2 not taken.
 117157 m_normalVertice.resize(m_numPoint,m_numCoor);
60
61
2/2
✓ Branch 0 taken 305235 times.
✓ Branch 1 taken 117157 times.
 422392 for(int ig=0; ig<m_numQuadraturePoint; ig++) {
62
1/2
✓ Branch 2 taken 305235 times.
✗ Branch 3 not taken.
 305235 m_covBasis[ig].resize(m_numRefCoor,m_numCoor);
63
1/2
✓ Branch 2 taken 305235 times.
✗ Branch 3 not taken.
 305235 m_F0[ig].resize(m_numCoor,m_numCoor);
64
1/2
✓ Branch 2 taken 305235 times.
✗ Branch 3 not taken.
 305235 m_F1[ig].resize(m_numCoor,m_numCoor);
65
1/2
✓ Branch 2 taken 305235 times.
✗ Branch 3 not taken.
 305235 m_covCompleteBasis[ig].resize(m_numCoor,m_numCoor);
66
1/2
✓ Branch 2 taken 305235 times.
✗ Branch 3 not taken.
 305235 contravariantCompleteBasis[ig].resize(m_numCoor,m_numCoor);
67
1/2
✓ Branch 2 taken 305235 times.
✗ Branch 3 not taken.
 305235 tangent[ig].resize(m_numCoor - 1,m_numCoor);
68
1/2
✓ Branch 2 taken 305235 times.
✗ Branch 3 not taken.
 305235 normal[ig].resize(m_numCoor);
69
1/2
✓ Branch 2 taken 305235 times.
✗ Branch 3 not taken.
 305235 covMetric[ig].resize(m_numRefCoor,m_numRefCoor);
70
1/2
✓ Branch 2 taken 305235 times.
✗ Branch 3 not taken.
 305235 contravariantMetric[ig].resize(m_numRefCoor,m_numRefCoor);
71 }
72 117157 }
73
74 /***********************************************************************************/
75 /***********************************************************************************/
76
77 462662 CurvilinearFiniteElement::~CurvilinearFiniteElement()
78 {
79
2/2
✓ Branch 0 taken 305235 times.
✓ Branch 1 taken 117157 times.
 844784 for(int i = 0; i < m_numQuadraturePoint; i++) {
80 610470 m_covBasis[i].clear();
81 610470 m_covCompleteBasis[i].clear();
82 610470 contravariantCompleteBasis[i].clear();
83 610470 covMetric[i].clear();
84 610470 contravariantMetric[i].clear();
85 610470 tangent[i].clear();
86 610470 normal[i].clear();
87 610470 m_F0[i].clear();
88 610470 m_F1[i].clear();
89 }
90 234314 m_covBasis.clear();
91 234314 m_covCompleteBasis.clear();
92 234314 contravariantCompleteBasis.clear();
93 234314 covMetric.clear();
94 234314 contravariantMetric.clear();
95 234314 tangent.clear();
96 234314 normal.clear();
97 234314 m_tangentVertice.clear();
98 234314 m_normalVertice.clear();
99 462662 }
100
101 /***********************************************************************************/
102 /***********************************************************************************/
103
104 7934853 void CurvilinearFiniteElement::computeCovariantBasis()
105 {
106 // Derivatives of geo map:
107
2/2
✓ Branch 0 taken 29757830 times.
✓ Branch 1 taken 7934853 times.
 37692683 for(int ig = 0; ig <m_numQuadraturePoint; ig++ ) {
108
2/4
✓ Branch 4 taken 29757830 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 29757830 times.
✗ Branch 8 not taken.
 29757830 noalias(m_covBasis[ig]) = prod(m_dPhiGeo[ig],m_point);
109 }
110 7934853 }
111
112 /***********************************************************************************/
113 /***********************************************************************************/
114
115 7934853 void CurvilinearFiniteElement::computeMeas()
116 {
117
2/2
✓ Branch 0 taken 29757830 times.
✓ Branch 1 taken 7934853 times.
 37692683 for(int ig = 0; ig <m_numQuadraturePoint; ig++ ) {
118
3/6
✓ Branch 4 taken 29757830 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 29757830 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 29757830 times.
✗ Branch 12 not taken.
 29757830 noalias(covMetric[ig]) = prod(m_covBasis[ig],trans(m_covBasis[ig]));
119 }
120
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 3730158 times.
✓ Branch 2 taken 4204695 times.
✗ Branch 3 not taken.
 7934853 switch (m_numRefCoor) {
121 ✗ case 0:
122 ✗ for(int ig = 0; ig <m_numQuadraturePoint; ig++ ) {
123 ✗ m_meas[ig] = 1.0;
124 ✗ weightMeas(ig) = 1.0;
125 }
126 ✗ break;
127 3730158 case 1:
128
2/2
✓ Branch 0 taken 7665926 times.
✓ Branch 1 taken 3730158 times.
 11396084 for(int ig = 0; ig <m_numQuadraturePoint; ig++ ) {
129 7665926 m_meas[ig] = std::sqrt(covMetric[ig](0,0));
130 7665926 weightMeas(ig) = m_meas(ig) * m_quadratureRule->weight(ig);
131 }
132 3730158 break;
133 4204695 case 2:
134
2/2
✓ Branch 0 taken 22091904 times.
✓ Branch 1 taken 4204695 times.
 26296599 for(int ig = 0; ig <m_numQuadraturePoint; ig++) {
135 22091904 m_meas[ig] = std::sqrt(covMetric[ig](0,0)*covMetric[ig](1,1) - covMetric[ig](0,1)*covMetric[ig](1,0));
136 22091904 weightMeas(ig) = m_meas(ig) * m_quadratureRule->weight(ig);
137 }
138 4204695 break;
139 ✗ default:
140 ✗ FEL_ERROR("CurvilinearFiniteElement::computeCovariantBasis: numRefCoor must be 0, 1 or 2");
141 break;
142 }
143 7934853 }
144
145 /***********************************************************************************/
146 /***********************************************************************************/
147
148 2760 void CurvilinearFiniteElement::computeContravariantMetric()
149 {
150
1/2
✓ Branch 1 taken 2760 times.
✗ Branch 2 not taken.
 2760 UBlasMatrix inv_cov_metric(m_numRefCoor,m_numRefCoor);
151 double det;
152
2/2
✓ Branch 0 taken 5520 times.
✓ Branch 1 taken 2760 times.
 8280 for(int ig = 0; ig <m_numQuadraturePoint; ig++ ) {
153
1/2
✓ Branch 2 taken 5520 times.
✗ Branch 3 not taken.
 5520 MathUtilities::InvertMatrix(covMetric[ig], inv_cov_metric, det);
154
2/4
✓ Branch 2 taken 5520 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 5520 times.
✗ Branch 6 not taken.
 5520 contravariantMetric[ig].assign(UBlasIdentityMatrix (m_numRefCoor));
155
2/4
✓ Branch 2 taken 5520 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 5520 times.
✗ Branch 7 not taken.
 5520 contravariantMetric[ig] = prod(inv_cov_metric, contravariantMetric[ig]);
156 }
157 2760 }
158
159 /***********************************************************************************/
160 /***********************************************************************************/
161
162 2760 void CurvilinearFiniteElement::computeContravariantBasis()
163 {
164
1/2
✓ Branch 1 taken 2760 times.
✗ Branch 2 not taken.
 2760 UBlasMatrix inv_cov_compete_basis(m_numCoor,m_numCoor);
165 double det;
166
2/2
✓ Branch 0 taken 5520 times.
✓ Branch 1 taken 2760 times.
 8280 for(int ig = 0; ig <m_numQuadraturePoint; ig++ ) {
167
2/2
✓ Branch 0 taken 5520 times.
✓ Branch 1 taken 5520 times.
 11040 for(int iRefcoor=0; iRefcoor<m_numRefCoor; iRefcoor++) {
168
2/2
✓ Branch 0 taken 11040 times.
✓ Branch 1 taken 5520 times.
 16560 for(int jcoor=0; jcoor<m_numCoor; jcoor++) {
169
2/4
✓ Branch 2 taken 11040 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 11040 times.
✗ Branch 7 not taken.
 11040 m_covCompleteBasis[ig](iRefcoor,jcoor) = m_covBasis[ig](iRefcoor,jcoor);
170 }
171 }
172
2/2
✓ Branch 0 taken 11040 times.
✓ Branch 1 taken 5520 times.
 16560 for(int jcoor=0; jcoor<m_numCoor; jcoor++) {
173
2/4
✓ Branch 2 taken 11040 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 11040 times.
✗ Branch 7 not taken.
 11040 m_covCompleteBasis[ig](m_numRefCoor,jcoor) = normal[ig](jcoor);
174 }
175
176
1/2
✓ Branch 2 taken 5520 times.
✗ Branch 3 not taken.
 5520 MathUtilities::InvertMatrix(m_covCompleteBasis[ig], inv_cov_compete_basis, det);
177
2/4
✓ Branch 2 taken 5520 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 5520 times.
✗ Branch 6 not taken.
 5520 contravariantCompleteBasis[ig].assign(UBlasIdentityMatrix(m_numCoor));
178
2/4
✓ Branch 2 taken 5520 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 5520 times.
✗ Branch 7 not taken.
 5520 contravariantCompleteBasis[ig] = prod(inv_cov_compete_basis, contravariantCompleteBasis[ig]);
179 }
180 2760 }
181
182 /***********************************************************************************/
183 /***********************************************************************************/
184
185 ✗ void CurvilinearFiniteElement::computeTangentNormal()
186 {
187 ✗ computeTangent();
188 ✗ computeNormal();
189 }
190
191 /***********************************************************************************/
192 /***********************************************************************************/
193
194 7005485 void CurvilinearFiniteElement::computeNormal()
195 {
196 double norm,n1,n2,n3;
197
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 3586100 times.
✓ Branch 2 taken 3419385 times.
✗ Branch 3 not taken.
 7005485 switch (m_numRefCoor) {
198 ✗ case 0:
199 ✗ break;
200 3586100 case 1:
201
2/2
✓ Branch 0 taken 7373400 times.
✓ Branch 1 taken 3586100 times.
 10959500 for(int ig = 0; ig <m_numQuadraturePoint; ig++) {
202 7373400 n1 = m_covBasis[ig](0,1);
203 7373400 n2 = - m_covBasis[ig](0,0);
204 7373400 norm = std::sqrt(n1*n1+n2*n2);
205 7373400 normal[ig](0) = m_sign*n1/norm;
206 7373400 normal[ig](1) = m_sign*n2/norm;
207 }
208 3586100 break;
209 3419385 case 2:
210
2/2
✓ Branch 0 taken 18464660 times.
✓ Branch 1 taken 3419385 times.
 21884045 for(int ig = 0; ig <m_numQuadraturePoint; ig++) {
211 18464660 n1 = m_covBasis[ig](0,1) * m_covBasis[ig](1,2) - m_covBasis[ig](0,2)*m_covBasis[ig](1,1);
212 18464660 n2 = m_covBasis[ig](0,2) * m_covBasis[ig](1,0) - m_covBasis[ig](0,0)*m_covBasis[ig](1,2);
213 18464660 n3 = m_covBasis[ig](0,0) * m_covBasis[ig](1,1) - m_covBasis[ig](0,1)*m_covBasis[ig](1,0);
214 18464660 norm = std::sqrt( n1 * n1 + n2 * n2 + n3 * n3 );
215 18464660 normal[ig](0) = m_sign*n1 / norm;
216 18464660 normal[ig](1) = m_sign*n2 / norm;
217 18464660 normal[ig](2) = m_sign*n3 / norm;
218 }
219 3419385 break;
220 ✗ default:
221 ✗ FEL_ERROR("CurvilinearFiniteElement::computeNormal: numRefCoor must be 0, 1 or 2");
222 break;
223 }
224 7005485 }
225
226 /***********************************************************************************/
227 /***********************************************************************************/
228
229 ✗ void CurvilinearFiniteElement::computeTangent()
230 {
231 double norm;
232 ✗ for(int ig = 0; ig <m_numQuadraturePoint; ig++ ) {
233 ✗ for(int irefcoor=0; irefcoor<m_numRefCoor; irefcoor++) {
234 ✗ norm = 0.0;
235 ✗ for(int jcoor=0; jcoor<m_numCoor; jcoor++) {
236 ✗ norm += m_covBasis[ig](irefcoor,jcoor)*m_covBasis[ig](irefcoor,jcoor);
237 }
238 ✗ norm = std::sqrt(norm);
239 ✗ for(int jcoor=0; jcoor<m_numCoor; jcoor++) {
240 ✗ tangent[ig](irefcoor,jcoor) = m_covBasis[ig](irefcoor,jcoor)/norm;
241 }
242 }
243 }
244 }
245
246 /***********************************************************************************/
247 /***********************************************************************************/
248
249 ✗ void CurvilinearFiniteElement::evaluateBasisFunctionsInPoint(const Point* p, std::vector<double>& values) const
250 {
251 // TODO: update these multiplications for being more efficient. Specially use bounded matrix
252 ✗ FEL_ASSERT(m_numCoor==3);
253 // 0. Computing x0
254 // For triangles is just m_point( 0, : );
255 // For quadrangles this is not true since, for csi = 0 we do not obtain a point, but we obtain
256 // the center of the quadrangle
257 ✗ UBlasBoundedVector<double, 3> x0;
258 ✗ for(std::size_t k = 0; k<3; k++) {
259 ✗ x0(k)=0;
260 ✗ for(int i(0); i<m_numDof; i++) {
261 ✗ x0(k)+=m_geoEle.basisFunction().phi(i, 0, 0, 0 )*m_point(i,k);
262 }
263 }
264
265 // 1. Inverse the mapping.
266 ✗ UBlasBoundedVector<double, 3> q;
267 ✗ for(std::size_t k = 0; k<3; k++) {
268 ✗ q(k) = p->coor(k) - x0(k);
269 }
270
271 ✗ const UBlasMatrix jacT= prod(m_dPhiGeo[0],m_point);
272
273 ✗ const UBlasMatrix jac = trans(jacT);
274 ✗ const UBlasVector qtilde = prod(jacT,q);
275
276 ✗ const UBlasBoundedMatrix<double, 2, 2> A = prod(jacT,jac);
277 double d;
278 ✗ UBlasBoundedMatrix<double, 2, 2> invA;
279 ✗ MathUtilities::InvertMatrix(A, invA, d);
280
281 // the point has been mapped from R^3 to the reference element.
282 ✗ const UBlasVector csi=prod(invA,qtilde);
283
284 // 2. evalute basis functions here
285 ✗ for(int i=0; i<m_numDof; i++ ) {
286 ✗ values[i]=m_refEle.basisFunction().phi(i, csi[0], csi[1], /*z*/ 0 );
287 }
288 }
289
290 /***********************************************************************************/
291 /***********************************************************************************/
292
293 540558 void CurvilinearFiniteElement::update(const int id, const std::vector<Point*>& point)
294 {
295 540558 updatePoint(point);
296 540558 m_hasMeas = false;
297 540558 m_hasNormal = false;
298 540558 m_hasTangent = false;
299 540558 m_hasQuadPtCoor = false;
300 540558 m_currentId = id;
301 540558 }
302
303 /***********************************************************************************/
304 /***********************************************************************************/
305
306 ✗ void CurvilinearFiniteElement::update(const int id, const UBlasMatrix& point, const std::vector<int>& ipt)
307 {
308 ✗ updatePoint(point, ipt);
309 ✗ m_hasMeas = false;
310 ✗ m_hasNormal = false;
311 ✗ m_hasTangent = false;
312 ✗ m_hasQuadPtCoor = false;
313 ✗ m_currentId = id;
314 }
315
316 /***********************************************************************************/
317 /***********************************************************************************/
318
319 430272 void CurvilinearFiniteElement::updateMeas(const int id, const std::vector<Point*>& point)
320 {
321 430272 updatePoint(point);
322 430272 m_hasMeas = true;
323 430272 m_hasNormal = false;
324 430272 m_hasTangent = false;
325 430272 m_hasQuadPtCoor = false;
326 430272 m_currentId = id;
327 430272 computeCovariantBasis();
328 430272 computeMeas();
329 430272 }
330
331 /***********************************************************************************/
332 /***********************************************************************************/
333
334 ✗ void CurvilinearFiniteElement::updateMeas(const int id, const UBlasMatrix& point, const std::vector<int>& ipt)
335 {
336 ✗ updatePoint(point, ipt);
337 ✗ m_hasMeas = true;
338 ✗ m_hasNormal = false;
339 ✗ m_hasTangent = false;
340 ✗ m_hasQuadPtCoor = false;
341 ✗ m_currentId = id;
342 ✗ computeCovariantBasis();
343 ✗ computeMeas();
344 }
345
346 /***********************************************************************************/
347 /***********************************************************************************/
348
349 499098 void CurvilinearFiniteElement::updateMeasQuadPt(const int id, const std::vector<Point*>& point)
350 {
351 499098 updatePoint(point);
352 499098 m_hasMeas = true;
353 499098 m_hasNormal = false;
354 499098 m_hasTangent = false;
355 499098 m_hasQuadPtCoor = true;
356 499098 m_currentId = id;
357 499098 computeCovariantBasis();
358 499098 computeMeas();
359 499098 computeCurrentQuadraturePoint();
360 499098 }
361
362 /***********************************************************************************/
363 /***********************************************************************************/
364
365 2894315 void CurvilinearFiniteElement::updateMeasNormal(const int id, const std::vector<Point*>& point)
366 {
367 2894315 updatePoint(point);
368 2894315 m_hasMeas = true;
369 2894315 m_hasNormal = true;
370 2894315 m_hasTangent = false;
371 2894315 m_hasQuadPtCoor = false;
372 2894315 m_currentId = id;
373 2894315 computeCovariantBasis();
374 2894315 computeMeas();
375 2894315 computeNormal();
376 2894315 }
377
378 /***********************************************************************************/
379 /***********************************************************************************/
380
381 26220 void CurvilinearFiniteElement::updateMeasNormalQuadPt(const int id, const std::vector<Point*>& point)
382 {
383 26220 updatePoint(point);
384 26220 m_hasMeas = true;
385 26220 m_hasNormal = true;
386 26220 m_hasTangent = false;
387 26220 m_hasQuadPtCoor = false;
388 26220 m_currentId = id;
389 26220 computeCovariantBasis();
390 26220 computeMeas();
391 26220 computeNormal();
392 26220 computeCurrentQuadraturePoint();
393 26220 }
394
395 /***********************************************************************************/
396 /***********************************************************************************/
397
398 4080132 void CurvilinearFiniteElement::updateMeasNormal(const int id, const UBlasMatrix& point, const std::vector<int>& ipt)
399 {
400 4080132 updatePoint(point, ipt);
401 4080132 m_hasMeas = true;
402 4080132 m_hasNormal = true;
403 4080132 m_hasTangent = false;
404 4080132 m_hasQuadPtCoor = false;
405 4080132 m_currentId = id;
406 4080132 computeCovariantBasis();
407 4080132 computeMeas();
408 4080132 computeNormal();
409 4080132 }
410
411 /***********************************************************************************/
412 /***********************************************************************************/
413
414 2760 void CurvilinearFiniteElement::updateMeasNormalContra(const int id, const std::vector<Point*>& point)
415 {
416 2760 updatePoint(point);
417 2760 m_hasMeas = true;
418 2760 m_hasNormal = true;
419 2760 m_hasTangent = false;
420 2760 m_hasQuadPtCoor = false;
421 2760 m_hasOriginalQuadPoint = true;
422 2760 m_currentId = id;
423 2760 computeCovariantBasis();
424 2760 computeMeas();
425 2760 computeNormal();
426 2760 computeContravariantMetric();
427 2760 computeContravariantBasis();
428 2760 }
429
430 /***********************************************************************************/
431 /***********************************************************************************/
432
433 ✗ void CurvilinearFiniteElement::updateMeasNormalContra(const int id, const UBlasMatrix& point, const std::vector<int>& ipt)
434 {
435 ✗ updatePoint(point, ipt);
436 ✗ m_hasMeas = true;
437 ✗ m_hasNormal = true;
438 ✗ m_hasTangent = false;
439 ✗ m_hasQuadPtCoor = false;
440 ✗ m_currentId = id;
441 ✗ computeCovariantBasis();
442 ✗ computeMeas();
443 ✗ computeNormal();
444 ✗ computeContravariantMetric();
445 ✗ computeContravariantBasis();
446 }
447
448 /***********************************************************************************/
449 /***********************************************************************************/
450
451 ✗ void CurvilinearFiniteElement::updateMeasNormalTangent(const int id, const std::vector<Point*>& point)
452 {
453 ✗ updatePoint(point);
454 ✗ m_hasMeas = true;
455 ✗ m_hasNormal = true;
456 ✗ m_hasTangent = true;
457 ✗ m_hasQuadPtCoor = false;
458 ✗ m_currentId = id;
459 ✗ computeCovariantBasis();
460 ✗ computeMeas();
461 ✗ computeTangentNormal();
462 }
463
464 /***********************************************************************************/
465 /***********************************************************************************/
466
467 ✗ void CurvilinearFiniteElement::updateMeasNormalTangent(const int id, const UBlasMatrix& point, const std::vector<int>& ipt)
468 {
469 ✗ updatePoint(point, ipt);
470 ✗ m_hasMeas = true;
471 ✗ m_hasNormal = true;
472 ✗ m_hasTangent = true;
473 ✗ m_hasQuadPtCoor = false;
474 ✗ m_currentId = id;
475 ✗ computeCovariantBasis();
476 ✗ computeMeas();
477 ✗ computeTangentNormal();
478 }
479
480 /***********************************************************************************/
481 /***********************************************************************************/
482
483 ✗ void CurvilinearFiniteElement::updateBasisAndNormalContra(const int id, const std::vector<Point*>& point)
484 {
485 ✗ updatePoint(point);
486 ✗ m_hasMeas = true;
487 ✗ m_hasNormal = true;
488 ✗ m_hasTangent = false;
489 ✗ m_hasQuadPtCoor = false;
490 ✗ m_hasOriginalQuadPoint = true;
491 ✗ m_currentId = id;
492
493 ✗ std::vector<Point> pts(m_numQuadraturePoint);
494 ✗ for(felInt ipt=0; ipt<m_numQuadraturePoint; ++ipt)
495 ✗ pts[ipt] = m_quadratureRule->quadraturePoint(ipt);
496
497 // recompute the basis function at the original quadrature point
498 ✗ updateBasisWithQuadPoint(pts);
499
500 // normal contra
501 ✗ computeCovariantBasis();
502 ✗ computeMeas();
503 ✗ computeNormal();
504 ✗ computeContravariantMetric();
505 ✗ computeContravariantBasis();
506 }
507
508 /***********************************************************************************/
509 /***********************************************************************************/
510
511 2056 void CurvilinearFiniteElement::updateSubElementMeasNormal(
512 const int id,
513 const std::vector<Point*>& ptElem,
514 const std::vector<Point*> ptSubElem
515 )
516 {
517
1/2
✓ Branch 1 taken 2056 times.
✗ Branch 2 not taken.
 2056 updatePoint(ptElem);
518 2056 m_hasMeas = true;
519 2056 m_hasNormal = true;
520 2056 m_hasTangent = false;
521 2056 m_hasQuadPtCoor = false;
522 2056 m_hasOriginalQuadPoint = false;
523 2056 m_currentId = id;
524
525 // find the coordinate of the sub element in the reference element of the element
526
1/2
✓ Branch 3 taken 2056 times.
✗ Branch 4 not taken.
 2056 std::vector<Point*> ptSubRef(ptSubElem.size());
527
1/4
✓ Branch 2 taken 2056 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 2056 switch(m_geoEle.shape().typeShape()) {
528 2056 case Segment: {
529
2/2
✓ Branch 1 taken 4112 times.
✓ Branch 2 taken 2056 times.
 6168 for(std::size_t ipt=0; ipt<ptSubElem.size(); ++ipt) {
530
1/2
✓ Branch 1 taken 4112 times.
✗ Branch 2 not taken.
 8224 ptSubRef[ipt] = new Point();
531 double dx, dy;
532 4112 dx = std::abs(ptElem[1]->x() - ptElem[0]->x());
533 4112 dy = std::abs(ptElem[1]->y() - ptElem[0]->y());
534
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4112 times.
 4112 if(dx >= dy)
535 ✗ ptSubRef[ipt]->x() = 2.*(ptSubElem[ipt]->x() - ptElem[0]->x())/(ptElem[1]->x() - ptElem[0]->x()) - 1;
536 else
537 4112 ptSubRef[ipt]->x() = 2.*(ptSubElem[ipt]->y() - ptElem[0]->y())/(ptElem[1]->y() - ptElem[0]->y()) - 1;
538
539 4112 ptSubRef[ipt]->y() = 0;
540 }
541 2056 break;
542 }
543 ✗ case NullShape:
544 case Node:
545 case Triangle: {
546 ✗ const double det1 = (ptElem[2]->y() - ptElem[0]->y())*(ptElem[1]->x() - ptElem[0]->x())
547 ✗ - (ptElem[2]->x() - ptElem[0]->x())*(ptElem[1]->y() - ptElem[0]->y());
548
549 ✗ const double det2 = (ptElem[2]->x() - ptElem[0]->x())*(ptElem[1]->z() - ptElem[0]->z())
550 ✗ - (ptElem[2]->z() - ptElem[0]->z())*(ptElem[1]->x() - ptElem[0]->x());
551
552 ✗ const double det3 = (ptElem[2]->z() - ptElem[0]->z())*(ptElem[1]->y() - ptElem[0]->y())
553 ✗ - (ptElem[2]->y() - ptElem[0]->y())*(ptElem[1]->z() - ptElem[0]->z());
554
555 ✗ if(std::abs(det1) >= std::abs(det2) && std::abs(det1) >= std::abs(det3)){
556 ✗ for(std::size_t ipt=0; ipt<ptSubElem.size(); ++ipt) {
557 ✗ ptSubRef[ipt] = new Point();
558
559 ✗ ptSubRef[ipt]->x() = (ptElem[2]->y() - ptElem[0]->y())*(ptSubElem[ipt]->x() - ptElem[0]->x())/det1
560 ✗ + (ptElem[0]->x() - ptElem[2]->x()) * (ptSubElem[ipt]->y() - ptElem[0]->y())/det1;
561
562 ✗ ptSubRef[ipt]->y() = (ptElem[0]->y() - ptElem[1]->y())*(ptSubElem[ipt]->x() - ptElem[0]->x())/det1
563 ✗ + (ptElem[1]->x() - ptElem[0]->x())*(ptSubElem[ipt]->y() - ptElem[0]->y())/det1;
564
565 ✗ ptSubRef[ipt]->z() = 0.0;
566 }
567 ✗ } else if(std::abs(det2) >= std::abs(det1) && std::abs(det2) >= std::abs(det3)){
568 ✗ for(std::size_t ipt=0; ipt<ptSubElem.size(); ++ipt) {
569 ✗ ptSubRef[ipt] = new Point();
570
571 ✗ ptSubRef[ipt]->x() = (ptElem[2]->x() - ptElem[0]->x())*(ptSubElem[ipt]->z() - ptElem[0]->z())/det2
572 ✗ + (ptElem[0]->z() - ptElem[2]->z()) * (ptSubElem[ipt]->x() - ptElem[0]->x())/det2;
573
574 ✗ ptSubRef[ipt]->y() = (ptElem[0]->x() - ptElem[1]->x())*(ptSubElem[ipt]->z() - ptElem[0]->z())/det2
575 ✗ + (ptElem[1]->z() - ptElem[0]->z())*(ptSubElem[ipt]->x() - ptElem[0]->x())/det2;
576
577 ✗ ptSubRef[ipt]->z() = 0.0;
578 }
579 } else {
580 ✗ for(std::size_t ipt=0; ipt<ptSubElem.size(); ++ipt) {
581 ✗ ptSubRef[ipt] = new Point();
582
583 ✗ ptSubRef[ipt]->x() = (ptElem[2]->z() - ptElem[0]->z())*(ptSubElem[ipt]->y() - ptElem[0]->y())/det3
584 ✗ + (ptElem[0]->y() - ptElem[2]->y()) * (ptSubElem[ipt]->z() - ptElem[0]->z())/det3;
585
586 ✗ ptSubRef[ipt]->y() = (ptElem[0]->z() - ptElem[1]->z())*(ptSubElem[ipt]->y() - ptElem[0]->y())/det3
587 ✗ + (ptElem[1]->y() - ptElem[0]->y())*(ptSubElem[ipt]->z() - ptElem[0]->z())/det3;
588
589 ✗ ptSubRef[ipt]->z() = 0.0;
590 }
591 }
592
593 ✗ break;
594 }
595 ✗ case Quadrilateral:
596 case Tetrahedron:
597 case Hexahedron:
598 case Prism:
599 case Pyramid: {
600 ✗ FEL_ERROR("The transformation is not implemented for this type of mesh cells");
601 break;
602 }
603 }
604
605 // create the curvilinear finite element for the sub element and update it
606
1/2
✓ Branch 2 taken 2056 times.
✗ Branch 3 not taken.
 4112 CurvilinearFiniteElement feSubElt(m_refEle, m_geoEle, (DegreeOfExactness) m_quadratureRule->degreeOfExactness());
607
608
1/2
✓ Branch 1 taken 2056 times.
✗ Branch 2 not taken.
 2056 feSubElt.updateMeasQuadPt(id, ptSubRef);
609
610 // delete sub element coordinates
611
2/2
✓ Branch 1 taken 4112 times.
✓ Branch 2 taken 2056 times.
 6168 for(std::size_t ipt=0; ipt<ptSubRef.size(); ++ipt)
612
1/2
✓ Branch 1 taken 4112 times.
✗ Branch 2 not taken.
 4112 delete ptSubRef[ipt];
613 2056 ptSubRef.clear();
614
615 // recompute the basis function at the new quadrature points for the element
616
1/2
✓ Branch 1 taken 2056 times.
✗ Branch 2 not taken.
 2056 updateBasisWithQuadPoint(feSubElt.currentQuadPoint);
617
618 // update the derivatives for the element
619
1/2
✓ Branch 1 taken 2056 times.
✗ Branch 2 not taken.
 2056 computeCovariantBasis();
620
1/2
✓ Branch 1 taken 2056 times.
✗ Branch 2 not taken.
 2056 computeMeas();
621
1/2
✓ Branch 1 taken 2056 times.
✗ Branch 2 not taken.
 2056 computeNormal();
622 // change the weight of the quadrature point
623
2/2
✓ Branch 1 taken 4112 times.
✓ Branch 2 taken 2056 times.
 6168 for(int ig=0; ig<feSubElt.numQuadraturePoint(); ++ig) {
624
3/6
✓ Branch 1 taken 4112 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4112 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4112 times.
✗ Branch 8 not taken.
 4112 weightMeas(ig) = m_meas(ig) * feSubElt.weightMeas(ig);
625 }
626
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2056 times.
 2056 for(int ig=feSubElt.numQuadraturePoint(); ig<m_numQuadraturePoint; ++ig) {
627 ✗ weightMeas(ig) = 0.0;
628 }
629 2056 }
630
631 /***********************************************************************************/
632 /***********************************************************************************/
633
634 ✗ void CurvilinearFiniteElement::print(int verbose,std::ostream& c) const
635 {
636 ✗ if(verbose) {
637 ✗ c << "CurvilinearFiniteElement: " << "\n";
638 ✗ m_printBase(verbose, c);
639
640 ✗ if( m_hasPoint && verbose > 1 ) c << " point = " << m_point << "\n";
641 ✗ if( m_hasMeas && verbose > 1 ) c << " meas = " << measure() << "\n";
642 ✗ if( verbose > 2 ) {
643 ✗ for(int ig = 0; ig <numQuadraturePoint(); ig++ ) {
644 ✗ c << " phi[ig=" <<ig<< "] = " << phi[ig] << "\n";
645 }
646 }
647 ✗ if( verbose > 2 ) {
648 ✗ for(int ig = 0; ig <numQuadraturePoint(); ig++ ) {
649 ✗ c << " phiGeo[ig=" <<ig<< "] = " << m_phiGeo[ig] << "\n";
650 }
651 }
652 ✗ if( verbose > 2 ) {
653 ✗ for(int ig = 0; ig <numQuadraturePoint(); ig++ ) {
654 ✗ c << " dphiGeo[ig=" <<ig<< "] = " << m_dPhiGeo[ig] << "\n";
655 }
656 }
657 ✗ if( verbose > 2 ) {
658 ✗ for(int ig = 0; ig <numQuadraturePoint(); ig++ ) {
659 ✗ c << " dphiRef[ig=" <<ig<< "] = " << m_dPhiRef[ig] << "\n";
660 }
661 }
662 ✗ if( m_hasNormal && verbose > 2 ) {
663 ✗ for(int ig = 0; ig <numQuadraturePoint(); ig++ ) {
664 ✗ c << " normal[ig=" <<ig<< "] = " << normal[ig] << "\n";
665 }
666 }
667 ✗ if( m_hasTangent && verbose > 2 ) {
668 ✗ for(int ig = 0; ig <numQuadraturePoint(); ig++ ) {
669 ✗ c << " tangent[ig=" <<ig<< "] = " << tangent[ig] << "\n";
670 }
671 }
672 ✗ if( m_hasNormal && verbose > 2 ) {
673 ✗ for(int ig = 0; ig <numQuadraturePoint(); ig++ ) {
674 ✗ c << " covBasis[ig=" <<ig<< "] = " << m_covBasis[ig] << "\n";
675 }
676 }
677 ✗ if( m_hasNormal && verbose > 2 ) {
678 ✗ for(int ig = 0; ig <numQuadraturePoint(); ig++ ) {
679 ✗ c << " covMetric[ig=" <<ig<< "] = " << covMetric[ig] << "\n";
680 }
681 }
682 ✗ if( m_hasQuadPtCoor && verbose > 2 ) {
683 ✗ for(int ig = 0; ig <numQuadraturePoint(); ig++ ) {
684 ✗ c << " contravariantMetric[ig=" <<ig<< "] = " << contravariantMetric[ig] << "\n";
685 }
686 }
687 ✗ if( m_hasQuadPtCoor && verbose > 2 ) {
688 ✗ for(int ig = 0; ig <numQuadraturePoint(); ig++ ) {
689 ✗ c << " currentQuadPoint[ig=" <<ig<< "] : ";
690 ✗ currentQuadPoint[ig].print(verbose,c);
691 }
692 }
693 ✗ std::cout << std::endl;
694 }
695 }
696
697 /***********************************************************************************/
698 /***********************************************************************************/
699
700 547178 void CurvilinearFiniteElement::updateNormalTangentMesh (const std::vector<Point*>& normalFE, const std::vector< std::vector<Point*> >& tangentFE)
701 {
702
2/2
✓ Branch 0 taken 1591016 times.
✓ Branch 1 taken 547178 times.
 2138194 for(int i=0; i<m_numPoint; i++) {
703
2/2
✓ Branch 0 taken 4661928 times.
✓ Branch 1 taken 1591016 times.
 6252944 for(int icoor=0; icoor<m_numCoor; icoor++) {
704 4661928 m_normalVertice(i,icoor) = normalFE[i]->coor(icoor);
705
2/2
✓ Branch 0 taken 9101616 times.
✓ Branch 1 taken 4661928 times.
 13763544 for(int j=0; j<m_numCoor - 1; j++)
706 9101616 m_tangentVertice[j](i,icoor) = tangentFE[j][i]->coor(icoor);
707 }
708 }
709 547178 m_hasNormal = true;
710 547178 m_hasTangent = true;
711 547178 }
712
713 /***********************************************************************************/
714 /***********************************************************************************/
715
716 547178 void CurvilinearFiniteElement::compute_weightMeas()
717 {
718
2/2
✓ Branch 0 taken 1591016 times.
✓ Branch 1 taken 547178 times.
 2138194 for(int ir = 0; ir<m_numQuadraturePoint; ir++)
719 1591016 weightMeas[ir] = m_quadratureRule->weight(ir);
720 547178 m_hasMeas = true;
721 547178 }
722
723 /***********************************************************************************/
724 /***********************************************************************************/
725
726 55560 void CurvilinearFiniteElement::compute_F0_F1_beam(const double thickness)
727 {
728 55560 int iRefcoor=0; // = r
729 double sum1, sum2, sum3;
730
2/2
✓ Branch 0 taken 111120 times.
✓ Branch 1 taken 55560 times.
 166680 for(int ir = 0; ir<m_numQuadraturePoint; ir++) {
731 111120 m_F1[ir].clear();
732
733
2/2
✓ Branch 0 taken 222240 times.
✓ Branch 1 taken 111120 times.
 333360 for(int jcoor=0; jcoor<m_numCoor; jcoor++){
734 222240 sum1= 0.0;
735 222240 sum2= 0.0;
736 222240 sum3= 0.0;
737
2/2
✓ Branch 0 taken 444480 times.
✓ Branch 1 taken 222240 times.
 666720 for(int iNode=0; iNode<m_numPoint; iNode++){
738 444480 sum1 += m_dPhiGeo[ir](iRefcoor,iNode)*m_point(iNode,jcoor);
739 444480 sum2 += 0.5*thickness*m_phiGeo[ir](iNode)*m_normalVertice(iNode,jcoor);
740 444480 sum3 += 0.5*thickness*m_dPhiGeo[ir](iRefcoor,iNode)*m_normalVertice(iNode,jcoor);
741 }
742 222240 m_F0[ir](0,jcoor) = sum1;
743 222240 m_F0[ir](1,jcoor) = sum2;
744 222240 m_F1[ir](0,jcoor) = sum3;
745 }
746 }
747 55560 }
748
749 /***********************************************************************************/
750 /***********************************************************************************/
751
752 491618 void CurvilinearFiniteElement::compute_F0_F1_thinShell(const double thickness)
753 {
754 std::array<double, 5> sum;
755
756
2/2
✓ Branch 0 taken 1479896 times.
✓ Branch 1 taken 491618 times.
 1971514 for(int ir = 0; ir<m_numQuadraturePoint; ir++) {
757
1/2
✓ Branch 2 taken 1479896 times.
✗ Branch 3 not taken.
 1479896 m_F1[ir].clear();
758
2/2
✓ Branch 0 taken 4439688 times.
✓ Branch 1 taken 1479896 times.
 5919584 for(int jcoor=0; jcoor<m_numCoor; jcoor++){
759
2/2
✓ Branch 0 taken 22198440 times.
✓ Branch 1 taken 4439688 times.
 26638128 for(int i=0 ; i<5; i++)
760 22198440 sum[i]= 0.0;
761
2/2
✓ Branch 0 taken 13379568 times.
✓ Branch 1 taken 4439688 times.
 17819256 for(int iNode=0; iNode<m_numPoint; iNode++){
762 // iRefcoor = 0 -> g_r // iRefcoor = 1 -> g_s
763
2/2
✓ Branch 0 taken 26759136 times.
✓ Branch 1 taken 13379568 times.
 40138704 for(int iRefcoor=0; iRefcoor<m_numRefCoor; iRefcoor++) {
764
2/4
✓ Branch 2 taken 26759136 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 26759136 times.
✗ Branch 6 not taken.
 26759136 sum[2*iRefcoor] += m_dPhiGeo[ir](iRefcoor,iNode)*m_point(iNode,jcoor);
765
2/4
✓ Branch 2 taken 26759136 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 26759136 times.
✗ Branch 6 not taken.
 26759136 sum[2*iRefcoor+1] += 0.5*thickness*m_dPhiGeo[ir](iRefcoor,iNode)*m_normalVertice(iNode,jcoor);
766 }
767 // g_z
768
2/4
✓ Branch 2 taken 13379568 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 13379568 times.
✗ Branch 6 not taken.
 13379568 sum[4] += 0.5*thickness*m_phiGeo[ir](iNode)*m_normalVertice(iNode,jcoor);
769 }
770
2/2
✓ Branch 0 taken 8879376 times.
✓ Branch 1 taken 4439688 times.
 13319064 for(int iRefcoor=0; iRefcoor<m_numRefCoor; iRefcoor++){
771
1/2
✓ Branch 3 taken 8879376 times.
✗ Branch 4 not taken.
 8879376 m_F0[ir](iRefcoor,jcoor) = sum[2*iRefcoor];
772
1/2
✓ Branch 3 taken 8879376 times.
✗ Branch 4 not taken.
 8879376 m_F1[ir](iRefcoor,jcoor) = sum[2*iRefcoor+1];
773 }
774
1/2
✓ Branch 3 taken 4439688 times.
✗ Branch 4 not taken.
 4439688 m_F0[ir](2,jcoor) = sum[4];
775 }
776 }
777 491618 }
778
779 /***********************************************************************************/
780 /***********************************************************************************/
781
782 795508 void CurvilinearFiniteElement::compute_Grr0_Grr1_Grr2(const UBlasMatrix& rF0, const UBlasMatrix& rF1)
783 {
784 795508 m_Grr0 = 0.0; m_Grr1 = 0.0; m_Grr2 = 0.0;
785
786
2/2
✓ Branch 0 taken 2330964 times.
✓ Branch 1 taken 795508 times.
 3126472 for(int jcoor = 0; jcoor<m_numCoor; jcoor++) {
787 2330964 m_Grr0 += rF0(0,jcoor)*rF0(0,jcoor);
788 2330964 m_Grr1 += 2.*rF0(0,jcoor)*rF1(0,jcoor);
789 2330964 m_Grr2 += rF1(0,jcoor)*rF1(0,jcoor);
790 }
791 795508 }
792
793 /***********************************************************************************/
794 /***********************************************************************************/
795
796 739948 void CurvilinearFiniteElement::compute_Gss0_Gss1_Gss2(const UBlasMatrix& rF0, const UBlasMatrix& rF1)
797 {
798 739948 m_Gss0 = 0.0; m_Gss1 = 0.0; m_Gss2 = 0.0;
799
800
2/2
✓ Branch 0 taken 2219844 times.
✓ Branch 1 taken 739948 times.
 2959792 for(int jcoor = 0; jcoor<m_numCoor; jcoor++) {
801 2219844 m_Gss0 += rF0(1,jcoor)*rF0(1,jcoor);
802 2219844 m_Gss1 += 2.*rF0(1,jcoor)*rF1(1,jcoor);
803 2219844 m_Gss2 += rF1(1,jcoor)*rF1(1,jcoor);
804 }
805 739948 }
806
807 /***********************************************************************************/
808 /***********************************************************************************/
809
810 739948 void CurvilinearFiniteElement::compute_Grs0_Grs1_Grs2(const UBlasMatrix& rF0, const UBlasMatrix& rF1)
811 {
812 739948 m_Grs0 = 0.0; m_Grs1 = 0.0; m_Grs2 = 0.0;
813
814
2/2
✓ Branch 0 taken 2219844 times.
✓ Branch 1 taken 739948 times.
 2959792 for(int jcoor = 0; jcoor<m_numCoor; jcoor++) {
815 2219844 m_Grs0 += rF0(0,jcoor)*rF0(1,jcoor);
816 2219844 m_Grs1 += rF0(0,jcoor)*rF1(1,jcoor)+rF0(1,jcoor)*rF1(0,jcoor);
817 2219844 m_Grs2 += rF1(0,jcoor)*rF1(1,jcoor);
818 }
819 739948 }
820
821 /***********************************************************************************/
822 /***********************************************************************************/
823
824 ✗ void CurvilinearFiniteElement::computeNormalThinShell()
825 {
826 double sum;
827
828 ✗ for(int ir = 0; ir<m_numQuadraturePoint; ir++) {
829 ✗ for(int jcoor=0; jcoor<m_numCoor; jcoor++){
830 ✗ sum= 0.0;
831 ✗ for(int iNode=0; iNode<m_numPoint; iNode++){
832 ✗ sum += m_phiGeo[ir](iNode)*m_normalVertice(iNode,jcoor);
833 }
834 ✗ normal[ir](jcoor) = sum;
835 }
836 }
837 }
838
839 /***********************************************************************************/
840 /***********************************************************************************/
841
842 55560 void CurvilinearFiniteElement::updateBeam (const int id, const std::vector<Point*>& point, const std::vector<Point*>& normalFE,const std::vector< std::vector<Point*> >& tangentFE, const double thickness)
843 {
844 55560 updatePoint(point);
845 55560 updateNormalTangentMesh (normalFE, tangentFE);
846 55560 m_hasNormal = true;
847 55560 m_hasTangent = true;
848 55560 m_hasQuadPtCoor = false;
849 55560 m_currentId = id;
850 55560 compute_weightMeas();
851 // compute convariant basis
852 55560 compute_F0_F1_beam(thickness);
853 55560 }
854
855 /***********************************************************************************/
856 /***********************************************************************************/
857
858 491618 void CurvilinearFiniteElement::updateThinShell(const int id, const std::vector<Point*>& point, const std::vector<Point*>& normalFE, const std::vector <std::vector<Point*> >& tangentFE, const double thickness)
859 {
860 491618 updatePoint(point);
861 491618 updateNormalTangentMesh (normalFE, tangentFE);
862 491618 m_hasNormal = true;
863 491618 m_hasTangent = true;
864 491618 m_hasQuadPtCoor = false;
865 491618 m_currentId = id;
866 491618 compute_weightMeas();
867 // compute covariant basis linear part
868 491618 compute_F0_F1_thinShell(thickness);
869 491618 }
870
871 /***********************************************************************************/
872 /***********************************************************************************/
873
874 ✗ bool CurvilinearFiniteElement::isInBallRange(double* c, double r)
875 {
876 //For each points of the element (3 for a triangle for example), compute the distance to the ball center.
877 //If the distance is shorter than R, the point is inside the ball
878 //We consider that the element is inside the ball if at least one element is part of the ball
879
880 ✗ for(std::size_t numPoint = 0; numPoint < m_point.size1(); numPoint++){
881 ✗ double distanceSquared = 0;
882 ✗ for(std::size_t numCoor = 0 ; numCoor < m_point.size2(); numCoor++){
883 ✗ const double coor = m_point(numPoint,numCoor);
884 ✗ const double ballCoor = c[numCoor];
885 //std::cout << coor << " - " << ballCoor << std::endl;
886 ✗ distanceSquared += (coor-ballCoor)*(coor-ballCoor);
887 }
888 ✗ if(distanceSquared < (r*r)){
889 //std::cout << "Found an element in ball range" << std::endl;
890 ✗ return true;
891 }
892 //std::cout << distanceSquared << std::endl;
893 //std::cout << std::endl;
894 }
895 //std::cout << "Did not found any element in ball range" << std::endl;
896 ✗ return false;
897 }
898
899 } // end class curvilinearFiniteElement
900