GCC Code Coverage Report

Directory: ./
File: FiniteElement/currentFiniteElement.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 156 295 52.9%
Branches: 103 365 28.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: J-F. Gerbeau
13 //
14
15 // System includes
16 #include <cmath>
17
18 // External includes
19
20 // Project includes
21 #include "FiniteElement/currentFiniteElement.hpp"
22 #include "Tools/math_utilities.hpp"
23
24 namespace felisce
25 {
26 53102 CurrentFiniteElement::CurrentFiniteElement(
27 const RefElement& refEle,
28 const GeoElement& geoEle,
29 const std::vector<DegreeOfExactness>& degOfExactness
30 53102 ) : CurBaseFiniteElement(refEle, geoEle, degOfExactness)
31 {
32
1/2
✓ Branch 1 taken 53102 times.
✗ Branch 2 not taken.
 53102 dPhi.resize(m_numQuadraturePoint);
33
1/2
✓ Branch 1 taken 53102 times.
✗ Branch 2 not taken.
 53102 m_jacobian.resize(m_numQuadraturePoint);
34
35
1/2
✓ Branch 1 taken 53102 times.
✗ Branch 2 not taken.
 53102 m_indexQuadPoint.resize(2);
36 53102 m_indexQuadPoint[0] = 0;
37 53102 m_indexQuadPoint[1] = m_numQuadraturePoint;
38
39
2/2
✓ Branch 2 taken 233314 times.
✓ Branch 3 taken 53102 times.
 286416 for(int ig=m_indexQuadPoint[0]; ig<m_indexQuadPoint[1]; ig++) {
40
1/2
✓ Branch 2 taken 233314 times.
✗ Branch 3 not taken.
 233314 m_jacobian[ig].resize(m_numCoor,m_numCoor);
41
1/2
✓ Branch 2 taken 233314 times.
✗ Branch 3 not taken.
 233314 dPhi[ig].resize(m_numRefCoor,m_numDof);
42 }
43 53102 }
44
45 /***********************************************************************************/
46 /***********************************************************************************/
47
48 13976 CurrentFiniteElement::CurrentFiniteElement(
49 const RefElement& refEle,
50 const GeoElement& geoEle,
51 const DegreeOfExactness& degOfExactness
52
2/4
✓ Branch 2 taken 13976 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 13976 times.
✗ Branch 6 not taken.
 13976 ) : CurrentFiniteElement(refEle, geoEle, std::vector<DegreeOfExactness>(1, degOfExactness))
53 {
54 13976 }
55
56 /***********************************************************************************/
57 /***********************************************************************************/
58
59 205972 CurrentFiniteElement::~CurrentFiniteElement()
60 {
61 138806 m_jacobian.clear();
62 138806 dPhi.clear();
63 205972 }
64
65 /***********************************************************************************/
66 /***********************************************************************************/
67
68 ✗ void CurrentFiniteElement::allocateOnDofsDataStructures()
69 {
70 //Experimental..it has been tested ONLY on P1 tetra, but it should be more general (matteo)
71 ✗ m_dPhiRefOnDofs.resize(m_numDof);
72 ✗ m_phiGeoOnDofs.resize(m_numDof);
73 ✗ m_dPhiGeoOnDofs.resize(m_numDof);
74 ✗ dPhiOnDofs.resize(m_numDof);
75 ✗ phiOnDofs.resize(m_numDof);
76 ✗ m_jacobianOnDofs.resize(m_numDof);
77 ✗ for ( int idof(0); idof<m_numDof; idof++ ) {
78 ✗ m_dPhiRefOnDofs[idof].resize(m_numRefCoor,m_numDof);
79 ✗ m_phiGeoOnDofs[idof].resize(m_numPoint);
80 ✗ m_dPhiGeoOnDofs[idof].resize(m_numRefCoor,m_numPoint);
81 ✗ phiOnDofs[idof].resize(m_numDof);
82 ✗ dPhiOnDofs[idof].resize(m_numRefCoor,m_numDof);
83 ✗ m_jacobianOnDofs[idof].resize(m_numCoor,m_numCoor);
84
85 ✗ QuadraturePoint point(0.,0.,0.,0.);
86 ✗ for (int pcoor=0; pcoor<m_numRefCoor; pcoor++ ) {
87 ✗ point.coor(pcoor)=m_refEle.nodeCoor(idof,pcoor);
88 }
89 ✗ for(int i=0; i<m_numDof; i++ ) {
90 ✗ phiOnDofs[idof](i) = m_refEle.basisFunction().phi(i, point);
91 ✗ for(int icoor=0; icoor<m_numRefCoor; icoor++ ) {
92 ✗ m_dPhiRefOnDofs[idof](icoor,i) = m_refEle.basisFunction().dPhi(i,icoor, point);
93 }
94 }
95 ✗ for(int k=0; k<m_numPoint; k++) {
96 ✗ m_phiGeoOnDofs[idof](k)= m_geoEle.basisFunction().phi(k, point);
97 ✗ for(int icoor=0; icoor<m_numRefCoor; icoor++) {
98 ✗ m_dPhiGeoOnDofs[idof](icoor,k) = m_geoEle.basisFunction().dPhi(k,icoor, point);
99 }
100 }
101 }
102 ✗ m_phiDPhiOnDofs=true;
103 }
104
105 /***********************************************************************************/
106 /***********************************************************************************/
107
108 42102782 void CurrentFiniteElement::computeJacobian(const int iBlockBd,const std::size_t numBlockBd)
109 {
110
3/6
✓ Branch 0 taken 42102782 times.
✗ Branch 1 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 42102782 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 42102782 times.
 42102782 FEL_ASSERT( numBlockBd > 0 && numBlockBd < m_indexQuadPoint.size() );
111 // derivatives of geo map:
112
2/2
✓ Branch 2 taken 208082119 times.
✓ Branch 3 taken 42102782 times.
 250184901 for ( int ig = m_indexQuadPoint[iBlockBd]; ig <m_indexQuadPoint[iBlockBd+numBlockBd]; ig++) {
113
1/2
✓ Branch 4 taken 208082119 times.
✗ Branch 5 not taken.
 208082119 m_jacobian[ig] = prod(m_dPhiGeo[ig],m_point);
114 }
115 42102782 }
116
117 /***********************************************************************************/
118 /***********************************************************************************/
119
120 ✗ void CurrentFiniteElement::computeJacobianOnDofs()
121 {
122 // derivatives of geo map:
123 ✗ for ( int idof(0); idof<m_numDof; ++idof) {
124 ✗ m_jacobianOnDofs[idof] = prod(m_dPhiGeoOnDofs[idof],m_point);
125 }
126 }
127
128 /***********************************************************************************/
129 /***********************************************************************************/
130
131 130177 void CurrentFiniteElement::computeMeas()
132 {
133
2/2
✓ Branch 0 taken 722153 times.
✓ Branch 1 taken 130177 times.
 852330 for (int ig = 0; ig <m_numQuadraturePoint; ig++) {
134 722153 m_meas(ig) = MathUtilities::Det(m_jacobian[ig]);
135 722153 weightMeas(ig) = m_meas(ig) * m_quadratureRule->weight(ig);
136 }
137 130177 }
138
139 /***********************************************************************************/
140 /***********************************************************************************/
141
142 41955997 void CurrentFiniteElement::computeMeasDer(const int iBlockBd,const std::size_t numBlockBd)
143 {
144
3/6
✓ Branch 0 taken 41955997 times.
✗ Branch 1 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 41955997 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 41955997 times.
 41955997 FEL_ASSERT( numBlockBd > 0 && numBlockBd < m_indexQuadPoint.size() );
145
1/2
✓ Branch 1 taken 41955997 times.
✗ Branch 2 not taken.
 41955997 UBlasMatrix inv_jacobian(m_numCoor,m_numCoor);
146 double det;
147
2/2
✓ Branch 2 taken 207310142 times.
✓ Branch 3 taken 41955997 times.
 249266139 for (int ig = m_indexQuadPoint[iBlockBd]; ig <m_indexQuadPoint[iBlockBd+numBlockBd]; ig++) {
148
1/2
✓ Branch 2 taken 207310142 times.
✗ Branch 3 not taken.
 207310142 MathUtilities::InvertMatrix(m_jacobian[ig], inv_jacobian, det);
149 //! on internal quad points only (boundary quad pts are set from CurvilinearFE):
150
3/4
✓ Branch 0 taken 207310142 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 196054166 times.
✓ Branch 3 taken 11255976 times.
 207310142 if ( ig >=0 && ig < m_numQuadraturePoint ) {
151
1/2
✓ Branch 1 taken 196054166 times.
✗ Branch 2 not taken.
 196054166 m_meas(ig) = det;
152
3/6
✓ Branch 1 taken 196054166 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 196054166 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 196054166 times.
✗ Branch 8 not taken.
 196054166 weightMeas(ig) = m_meas(ig) * m_quadratureRule->weight(ig);
153 }
154
3/6
✓ Branch 2 taken 207310142 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 207310142 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 207310142 times.
✗ Branch 10 not taken.
 207310142 noalias(dPhi[ig]) = prod(inv_jacobian, m_dPhiRef[ig]);
155 }
156 41955997 }
157
158 /***********************************************************************************/
159 /***********************************************************************************/
160
161 ✗ void CurrentFiniteElement::computeDerOnDofs()
162 {
163 ✗ UBlasMatrix inv_jacobian(m_numCoor,m_numCoor);
164 double det;
165 ✗ for (int idof(0); idof<m_numDof; ++idof) {
166 ✗ MathUtilities::InvertMatrix(m_jacobianOnDofs[idof], inv_jacobian, det);
167 ✗ noalias(dPhiOnDofs[idof]) = prod(inv_jacobian, dPhiOnDofs[idof]);
168 }
169 }
170
171 /***********************************************************************************/
172 /***********************************************************************************/
173
174 16608 void CurrentFiniteElement::computePiolaMeasDer(const int iBlockBd,const std::size_t numBlockBd, const std::vector<double>& normal)
175 {
176
3/6
✓ Branch 0 taken 16608 times.
✗ Branch 1 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 16608 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 16608 times.
 16608 FEL_ASSERT( numBlockBd > 0 && numBlockBd < m_indexQuadPoint.size() );
177
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 16608 times.
 16608 FEL_ASSERT( (felInt)normal.size() == m_numCoor);
178
179
1/2
✓ Branch 1 taken 16608 times.
✗ Branch 2 not taken.
 16608 UBlasMatrix inv_jacobian(m_numCoor,m_numCoor);
180
1/2
✓ Branch 1 taken 16608 times.
✗ Branch 2 not taken.
 16608 UBlasVector meas(m_numCoor);
181 double det;
182
183
2/2
✓ Branch 2 taken 49824 times.
✓ Branch 3 taken 16608 times.
 66432 for (int ig = m_indexQuadPoint[iBlockBd]; ig <m_indexQuadPoint[iBlockBd+numBlockBd]; ig++) {
184
2/2
✓ Branch 0 taken 99648 times.
✓ Branch 1 taken 49824 times.
 149472 for(felInt i=0; i<m_numCoor; ++i)
185
1/2
✓ Branch 2 taken 99648 times.
✗ Branch 3 not taken.
 99648 meas(i) = normal[i];
186
187
1/2
✓ Branch 2 taken 49824 times.
✗ Branch 3 not taken.
 49824 MathUtilities::InvertMatrix(m_jacobian[ig], inv_jacobian, det);
188
3/6
✓ Branch 2 taken 49824 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 49824 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 49824 times.
✗ Branch 10 not taken.
 49824 noalias(dPhi[ig]) = prod(inv_jacobian, m_dPhiRef[ig]);
189
190
2/4
✓ Branch 0 taken 49824 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 49824 times.
✗ Branch 3 not taken.
 49824 if ( ig >=0 && ig < m_numQuadraturePoint )
191
1/2
✓ Branch 1 taken 49824 times.
✗ Branch 2 not taken.
 49824 m_meas(ig) = det;
192
193 // Now, the norm of the inverse of the transpose of the jacobian times the normal vector
194
2/4
✓ Branch 1 taken 49824 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 49824 times.
✗ Branch 5 not taken.
 49824 meas = prod(inv_jacobian, meas);
195
196 //! On internal quad points only (boundary quad pts are set from CurvilinearFE):
197
2/4
✓ Branch 0 taken 49824 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 49824 times.
✗ Branch 3 not taken.
 49824 if ( ig >=0 && ig < m_numQuadraturePoint ) {
198
2/4
✓ Branch 1 taken 49824 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 49824 times.
✗ Branch 5 not taken.
 49824 m_meas(ig) *= norm_2(meas);
199
3/6
✓ Branch 1 taken 49824 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 49824 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 49824 times.
✗ Branch 8 not taken.
 49824 weightMeas(ig) = m_meas(ig) * m_quadratureRule->weight(ig);
200 }
201 }
202 16608 }
203
204 /***********************************************************************************/
205 /***********************************************************************************/
206
207 4583540 void CurrentFiniteElement::update(const int id,const std::vector<Point*>& point)
208 {
209 4583540 updatePoint(point);
210 4583540 m_hasMeas = false;
211 4583540 m_hasFirstDeriv = false;
212 4583540 m_hasSecondDeriv = false;
213 4583540 m_hasQuadPtCoor = false;
214 4583540 m_currentId = id;
215 4583540 }
216
217 /***********************************************************************************/
218 /***********************************************************************************/
219
220 110481 void CurrentFiniteElement::updateMeas(const int id,const std::vector<Point*>& point)
221 {
222 110481 updatePoint(point);
223 110481 m_hasMeas = true;
224 110481 m_hasFirstDeriv = false;
225 110481 m_hasSecondDeriv = false;
226 110481 m_hasQuadPtCoor = false;
227 110481 m_currentId = id;
228 110481 computeJacobian(0,m_indexQuadPoint.size()-1); // on all quad points
229 110481 computeMeas();
230 110481 }
231
232 /***********************************************************************************/
233 /***********************************************************************************/
234
235 19696 void CurrentFiniteElement::updateMeasQuadPt(const int id,const std::vector<Point*>& point)
236 {
237 19696 updatePoint(point);
238 19696 m_hasMeas = true;
239 19696 m_hasFirstDeriv = false;
240 19696 m_hasSecondDeriv = false;
241 19696 m_hasQuadPtCoor = false; // std::set to "true" in computeCurrentQuadraturePoint()
242 19696 m_currentId = id;
243 19696 computeJacobian(0,m_indexQuadPoint.size()-1); // on all quad points
244 19696 computeMeas();
245 19696 computeCurrentQuadraturePoint();
246 19696 }
247
248 /***********************************************************************************/
249 /***********************************************************************************/
250
251 41941773 void CurrentFiniteElement::updateFirstDeriv(const int id,const std::vector<Point*>& point)
252 {
253 41941773 updatePoint(point);
254 41941773 m_hasMeas = true;
255 41941773 m_hasFirstDeriv = true;
256 41941773 m_hasSecondDeriv = false;
257 41941773 m_hasQuadPtCoor = false;
258 41941773 m_currentId = id;
259 41941773 computeJacobian(0,m_indexQuadPoint.size()-1); // on all quad points
260 41941773 computeMeasDer(0,m_indexQuadPoint.size()-1);
261 41941773 }
262
263 /***********************************************************************************/
264 /***********************************************************************************/
265
266 ✗ void CurrentFiniteElement::updateFirstDerivOnDofs(const int id,const std::vector<Point*>& point)
267 {
268 ✗ updateFirstDeriv(id,point);
269 ✗ if ( m_phiDPhiOnDofs == false ) {
270 ✗ FEL_ERROR("Error: you should call allocateOnDofsDataStructures");
271 }
272 ✗ computeJacobianOnDofs();
273 ✗ computeDerOnDofs();
274 }
275
276 /***********************************************************************************/
277 /***********************************************************************************/
278
279 2720 void CurrentFiniteElement::updateBasisAndFirstDeriv(const int id, const std::vector<Point*>& point) {
280
1/2
✓ Branch 1 taken 2720 times.
✗ Branch 2 not taken.
 2720 updatePoint(point);
281 2720 m_hasMeas = true;
282 2720 m_hasFirstDeriv = true;
283 2720 m_hasSecondDeriv = false;
284 2720 m_hasQuadPtCoor = false;
285 2720 m_hasOriginalQuadPoint = true;
286 2720 m_currentId = id;
287
288
1/2
✓ Branch 2 taken 2720 times.
✗ Branch 3 not taken.
 2720 std::vector<Point> pts(m_numQuadraturePoint);
289
2/2
✓ Branch 0 taken 8160 times.
✓ Branch 1 taken 2720 times.
 10880 for(felInt ipt=0; ipt<m_numQuadraturePoint; ++ipt)
290
1/2
✓ Branch 1 taken 8160 times.
✗ Branch 2 not taken.
 8160 pts[ipt] = m_quadratureRule->quadraturePoint(ipt);
291
292 // recompute the basis function at the original quadrature point
293
1/2
✓ Branch 1 taken 2720 times.
✗ Branch 2 not taken.
 2720 updateBasisWithQuadPoint(pts);
294
295 // first deriv
296
1/2
✓ Branch 2 taken 2720 times.
✗ Branch 3 not taken.
 2720 computeJacobian(0, m_indexQuadPoint.size()-1);
297
298
1/2
✓ Branch 2 taken 2720 times.
✗ Branch 3 not taken.
 2720 computeMeasDer(0, m_indexQuadPoint.size()-1);
299 2720 }
300
301 /***********************************************************************************/
302 /***********************************************************************************/
303
304 28112 void CurrentFiniteElement::updateSubElementFirstDeriv(const int id, const std::vector<Point*>& ptElem, const std::vector<Point*>& ptSubElem, CurBaseFiniteElement* feSurf)
305 {
306
1/2
✓ Branch 1 taken 28112 times.
✗ Branch 2 not taken.
 28112 updatePoint(ptElem);
307 28112 m_hasMeas = true;
308 28112 m_hasFirstDeriv = true;
309 28112 m_hasSecondDeriv = false;
310 28112 m_hasQuadPtCoor = false;
311 28112 m_hasOriginalQuadPoint = false;
312 28112 m_currentId = id;
313
314 // find the coordinate of the sub element in the reference element of the element
315
1/2
✓ Branch 3 taken 28112 times.
✗ Branch 4 not taken.
 28112 std::vector<Point*> ptSubRef(ptSubElem.size());
316
1/5
✗ Branch 2 not taken.
✓ Branch 3 taken 28112 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 28112 switch(m_geoEle.shape().typeShape()) {
317
318 ✗ case Quadrilateral: {
319 // square
320 ✗ for(std::size_t ipt=0; ipt<ptSubElem.size(); ++ipt) {
321 ✗ ptSubRef[ipt] = new Point();
322 ✗ ptSubRef[ipt]->x() = 2.*(ptSubElem[ipt]->x() - ptElem[0]->x())/(ptElem[2]->x() - ptElem[0]->x()) - 1.;
323 ✗ ptSubRef[ipt]->y() = 2.*(ptSubElem[ipt]->y() - ptElem[0]->y())/(ptElem[2]->y() - ptElem[0]->y()) - 1.;
324 }
325 ✗ break;
326 }
327
328 28112 case Triangle: {
329 // triangle
330 // this value can not be equal to zero unless the triangle is flat.
331 28112 const double det = (ptElem[2]->y() - ptElem[0]->y())*(ptElem[1]->x() - ptElem[0]->x())
332 28112 - (ptElem[2]->x() - ptElem[0]->x())*(ptElem[1]->y() - ptElem[0]->y());
333
334
2/2
✓ Branch 1 taken 67728 times.
✓ Branch 2 taken 28112 times.
 95840 for(std::size_t ipt=0; ipt<ptSubElem.size(); ++ipt) {
335
1/2
✓ Branch 1 taken 67728 times.
✗ Branch 2 not taken.
 135456 ptSubRef[ipt] = new Point();
336
337 135456 ptSubRef[ipt]->x() = (ptElem[2]->y() - ptElem[0]->y())*(ptSubElem[ipt]->x() - ptElem[0]->x())/det
338 67728 + (ptElem[0]->x() - ptElem[2]->x())*(ptSubElem[ipt]->y() - ptElem[0]->y())/det;
339
340 135456 ptSubRef[ipt]->y() = (ptElem[0]->y() - ptElem[1]->y())*(ptSubElem[ipt]->x() - ptElem[0]->x())/det
341 67728 + (ptElem[1]->x() - ptElem[0]->x())*(ptSubElem[ipt]->y() - ptElem[0]->y())/det;
342 }
343 28112 break;
344 }
345
346 ✗ case Tetrahedron: {
347 // coordinates of the subElement in the reference Tetrahedra
348 // inversion of the map matrix
349 //std::cout << " Vertex Current Elem K: ";
350 //for( int it=0; it<ptElem.size(); ++it){
351 // ptElem[it]->print(10,std::cout);
352 // std::cout<<std::endl;
353 //}
354 //std::cout << " Vertex Current SubElem K_i: ";
355 //for( int it=0; it<ptSubElem.size(); ++it){
356 // ptSubElem[it]->print(10,std::cout);
357 // std::cout<<std::endl;
358 //}
359
360 ✗ UBlasMatrix FLU(m_numCoor,m_numCoor);
361 ✗ UBlasMatrix inv_FLU(m_numCoor,m_numCoor);
362 ✗ UBlasVector refPoint(m_numCoor);
363 double det;
364
365 // Filling the matrix
366 ✗ for(int indx = 0; indx < m_numCoor ; ++indx){
367 ✗ FLU(0,indx) = ptElem[indx+1]->x() - ptElem[0]->x();
368 ✗ FLU(1,indx) = ptElem[indx+1]->y() - ptElem[0]->y();
369 ✗ FLU(2,indx) = ptElem[indx+1]->z() - ptElem[0]->z();
370 }
371
372 ✗ MathUtilities::InvertMatrix(FLU, inv_FLU, det);
373
374 ✗ for(std::size_t ipt=0; ipt<ptSubElem.size(); ++ipt) {
375 // Filling the rhs
376 ✗ refPoint[0] = ptSubElem[ipt]->x() - ptElem[0]->x();
377 ✗ refPoint[1] = ptSubElem[ipt]->y() - ptElem[0]->y();
378 ✗ refPoint[2] = ptSubElem[ipt]->z() - ptElem[0]->z();
379
380 ✗ refPoint = prod(inv_FLU,refPoint);
381
382 ✗ ptSubRef[ipt] = new Point();
383 ✗ ptSubRef[ipt]->x() = refPoint[0];
384 ✗ ptSubRef[ipt]->y() = refPoint[1];
385 ✗ ptSubRef[ipt]->z() = refPoint[2];
386 }
387
388 ✗ break;
389 }
390
391 ✗ case NullShape:
392 case Node:
393 case Segment:
394 case Hexahedron:
395 case Prism:
396 case Pyramid: {
397 ✗ FEL_ERROR("The transformation is not implemented for this type of mesh cells");
398 break;
399 }
400 }
401
402 // create the current finite element for the sub element and update it
403 CurBaseFiniteElement* feSubElt;
404
2/2
✓ Branch 0 taken 11504 times.
✓ Branch 1 taken 16608 times.
 28112 if(feSurf == nullptr) {
405 // volume sub elem
406
2/4
✓ Branch 2 taken 11504 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 11504 times.
✗ Branch 6 not taken.
 11504 feSubElt = new CurrentFiniteElement(m_refEle, m_geoEle, (DegreeOfExactness) m_quadratureRule->degreeOfExactness());
407
1/2
✓ Branch 1 taken 11504 times.
✗ Branch 2 not taken.
 11504 ((CurrentFiniteElement*) feSubElt)->updateMeasQuadPt(id, ptSubRef);
408 } else {
409 // curvi sub elem
410
2/4
✓ Branch 5 taken 16608 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 16608 times.
✗ Branch 9 not taken.
 16608 feSubElt = new CurvilinearFiniteElement(feSurf->refEle(), feSurf->geoEle(), (DegreeOfExactness) feSurf->quadratureRule().degreeOfExactness());
411
1/2
✓ Branch 1 taken 16608 times.
✗ Branch 2 not taken.
 16608 ((CurvilinearFiniteElement*) feSubElt)->updateMeasQuadPt(id, ptSubRef);
412 }
413
414 // recompute the basis function at the new quadrature points for the element
415
1/2
✓ Branch 1 taken 28112 times.
✗ Branch 2 not taken.
 28112 updateBasisWithQuadPoint(feSubElt->currentQuadPoint);
416
417 // update the derivatives for the element (no need to recompute the jacobian, it's
418 // a constant if the mapping from the reference element to the current element is affine)
419
1/2
✓ Branch 2 taken 28112 times.
✗ Branch 3 not taken.
 28112 computeJacobian(0, m_indexQuadPoint.size()-1);
420
421
2/2
✓ Branch 0 taken 11504 times.
✓ Branch 1 taken 16608 times.
 28112 if(feSurf == nullptr){
422
423
1/2
✓ Branch 2 taken 11504 times.
✗ Branch 3 not taken.
 11504 computeMeasDer(0, m_indexQuadPoint.size()-1);
424 }
425 else {
426
1/2
✓ Branch 2 taken 16608 times.
✗ Branch 3 not taken.
 16608 std::vector<double> normal(m_numCoor);
427 double norm;
428
1/2
✓ Branch 0 taken 16608 times.
✗ Branch 1 not taken.
 16608 if(m_numCoor == 2) {
429 16608 normal[0] = ptSubRef[0]->y() - ptSubRef[1]->y();
430 16608 normal[1] = ptSubRef[1]->x() - ptSubRef[0]->x();
431 16608 norm = std::sqrt(normal[0]*normal[0] + normal[1]*normal[1]);
432 16608 normal[0] /= norm;
433 16608 normal[1] /= norm;
434 } else {
435 // normal in 3d
436 ✗ normal[0] = (ptSubRef[1]->y() - ptSubRef[0]->y())*(ptSubRef[2]->z() - ptSubRef[0]->z()) - ((ptSubRef[1]->z() - ptSubRef[0]->z())*(ptSubRef[2]->y() - ptSubRef[0]->y()));
437 ✗ normal[1] = (ptSubRef[1]->z() - ptSubRef[0]->z())*(ptSubRef[2]->x() - ptSubRef[0]->x()) - ((ptSubRef[1]->x() - ptSubRef[0]->x())*(ptSubRef[2]->z() - ptSubRef[0]->z()));
438 ✗ normal[2] = (ptSubRef[1]->x() - ptSubRef[0]->x())*(ptSubRef[2]->y() - ptSubRef[0]->y()) - ((ptSubRef[2]->x() - ptSubRef[0]->x())*(ptSubRef[1]->y() - ptSubRef[0]->y()));
439 ✗ norm = std::sqrt(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);
440 ✗ normal[0] /= norm;
441 ✗ normal[1] /= norm;
442 ✗ normal[2] /= norm;
443 }
444
445
1/2
✓ Branch 2 taken 16608 times.
✗ Branch 3 not taken.
 16608 computePiolaMeasDer(0, m_indexQuadPoint.size()-1, normal);
446 16608 }
447
448 // change the weight of the quadrature point
449
2/2
✓ Branch 1 taken 67728 times.
✓ Branch 2 taken 28112 times.
 95840 for(int ig=0; ig<feSubElt->numQuadraturePoint(); ++ig) {
450
3/6
✓ Branch 1 taken 67728 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 67728 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 67728 times.
✗ Branch 8 not taken.
 67728 weightMeas(ig) = m_meas(ig) * feSubElt->weightMeas(ig);
451 }
452
2/2
✓ Branch 1 taken 16608 times.
✓ Branch 2 taken 28112 times.
 44720 for(int ig=feSubElt->numQuadraturePoint(); ig<m_numQuadraturePoint; ++ig) {
453
1/2
✓ Branch 1 taken 16608 times.
✗ Branch 2 not taken.
 16608 weightMeas(ig) = 0.;
454 }
455
456 // delete
457
2/2
✓ Branch 1 taken 67728 times.
✓ Branch 2 taken 28112 times.
 95840 for(std::size_t ipt=0; ipt<ptSubRef.size(); ++ipt)
458
1/2
✓ Branch 1 taken 67728 times.
✗ Branch 2 not taken.
 67728 delete ptSubRef[ipt];
459 28112 ptSubRef.clear();
460
461
1/2
✓ Branch 0 taken 28112 times.
✗ Branch 1 not taken.
 28112 delete feSubElt;
462 28112 }
463
464 /***********************************************************************************/
465 /***********************************************************************************/
466
467 ✗ void CurrentFiniteElement::identifyLocBDDof(const CurvilinearFiniteElement& curvFe, std::vector<int>& map) const
468 {
469 // This function compare the point supporting dofs of a curvilinear finite element,
470 // with the point supporting dofs of "this" a current finite element.
471 // the return value "std::unordered_map" connects the local ids of the dofs on the curvFe to the loc ids of the dofs
472 // on the currentFiniteElement.
473 // You probably do not want to use this function and you would like to use the currentFiniteElementWithBD
474 // Which does similar things in a better way.
475 ✗ map.resize(curvFe.numDof());
476 Point p,pref,pv,pvref;
477 ✗ for (int iBD(0); iBD<curvFe.numDof();iBD++){
478 ✗ map[iBD]=-1;
479 ✗ for (int pcoor=0; pcoor<m_numRefCoor-1; pcoor++ ) {
480 ✗ pref.coor(pcoor)=curvFe.refEle().nodeCoor(iBD,pcoor);
481 }
482 ✗ curvFe.coorMap(p,pref);
483 ✗ for (int iV(0); iV<this->numDof();iV++){
484 ✗ for (int pcoor=0; pcoor<m_numRefCoor; pcoor++ ) {
485 ✗ pvref.coor(pcoor)=m_refEle.nodeCoor(iV,pcoor);
486 }
487 ✗ coorMap(pv,pvref);
488 ✗ if ( p == pv ) {
489 ✗ map[iBD]=iV;
490 }
491 }
492 ✗ if (map[iBD]==-1) {
493 ✗ FEL_ERROR("err");
494 }
495 }
496 }
497
498 /***********************************************************************************/
499 /***********************************************************************************/
500
501 ✗ void CurrentFiniteElement::print(int verbose,std::ostream& c) const
502 {
503 ✗ if(verbose) {
504 ✗ c << "CurrentFiniteElement: " << "\n";
505 ✗ m_printBase(verbose, c);
506 ✗ int theNumQuadPointTotal = m_indexQuadPoint.back();
507 ✗ if ( theNumQuadPointTotal != m_numQuadraturePoint ) {
508 ✗ c << " numQuadPoint: " << m_numQuadraturePoint << "\n";
509 ✗ c << " numQuadPointTotal(intern+boundary): " << theNumQuadPointTotal << "\n";
510 }
511
512 ✗ if( m_hasPoint && verbose > 1 ) c << " point = " << m_point << "\n";
513 ✗ if( m_hasMeas && verbose > 1 ) c << " meas = " << measure() << "\n";
514 ✗ if( verbose > 2 ) {
515 ✗ if ( theNumQuadPointTotal != m_numQuadraturePoint ) {
516 ✗ c << " The numQuadPoint=" << m_numQuadraturePoint << " first are internal. The others live on the boundary elements of the element.\n";
517 }
518 ✗ for (int ig = 0; ig <theNumQuadPointTotal; ig++ ) {
519 ✗ c << " phi[ig=" <<ig<< "] = " << phi[ig] << "\n";
520 }
521 }
522 ✗ if( verbose > 2 ) {
523 ✗ for (int ig = 0; ig <theNumQuadPointTotal; ig++ ) {
524 ✗ c << " phiGeo[ig=" <<ig<< "] = " << m_phiGeo[ig] << "\n";
525 }
526 }
527 ✗ if( m_hasFirstDeriv && verbose > 2 ) {
528 ✗ for (int ig = 0; ig <theNumQuadPointTotal; ig++ ) {
529 ✗ c << " dphi[ig=" <<ig<< "] = " << dPhi[ig] << "\n";
530 }
531 }
532 ✗ if( verbose > 2 ) {
533 ✗ for (int ig = 0; ig <theNumQuadPointTotal; ig++ ) {
534 ✗ c << " dphiGeo[ig=" <<ig<< "] = " << m_dPhiGeo[ig] << "\n";
535 }
536 }
537 ✗ if( verbose > 2 ) {
538 ✗ for (int ig = 0; ig <theNumQuadPointTotal; ig++ ) {
539 ✗ c << " dphiRef[ig=" <<ig<< "] = " << m_dPhiRef[ig] << "\n";
540 }
541 }
542 ✗ if( m_hasFirstDeriv && verbose > 2 ) {
543 ✗ for (int ig = 0; ig <theNumQuadPointTotal; ig++ ) {
544 ✗ c << " jacobian[ig=" <<ig<< "] = " << m_jacobian[ig] << "\n";
545 }
546 }
547 ✗ if( m_hasQuadPtCoor && verbose > 2 ) {
548 ✗ for (int ig = 0; ig <theNumQuadPointTotal; ig++ ) {
549 ✗ c << " currentQuadPoint[ig=" <<ig<< "] : ";
550 ✗ currentQuadPoint[ig].print(verbose,c);
551 }
552 }
553 ✗ if( m_phiDPhiOnDofs && verbose > 2 ) {
554 ✗ for (int idof = 0; idof <m_numDof; idof++ ) {
555 ✗ c << " phiOnDofs[idof=" <<idof<< "] : "<<phiOnDofs[idof]<<std::endl;
556 }
557 ✗ for (int idof = 0; idof <m_numDof; idof++ ) {
558 ✗ c << " dPhiOnDofs[idof=" <<idof<< "] : "<<dPhiOnDofs[idof]<<std::endl;
559 }
560 ✗ for (int idof = 0; idof <m_numDof; idof++ ) {
561 ✗ c << " jacobianOnDodfs[idof=" <<idof<< "] : "<<m_jacobianOnDofs[idof]<<std::endl;
562 }
563 ✗ for (int idof = 0; idof <m_numDof; idof++ ) {
564 ✗ c << " phiGeoOnDofs[idof=" <<idof<< "] : "<<m_phiGeoOnDofs[idof]<<std::endl;
565 }
566 ✗ for (int idof = 0; idof <m_numDof; idof++ ) {
567 ✗ c << " dPhiGeoOnDofs[idof=" <<idof<< "] : "<<m_dPhiGeoOnDofs[idof]<<std::endl;
568 }
569 ✗ for (int idof = 0; idof <m_numDof; idof++ ) {
570 ✗ c << " dPhiRefOnDofs[idof=" <<idof<< "] : "<<m_dPhiRefOnDofs[idof]<<std::endl;
571 }
572 }
573 ✗ if ( verbose > 3) {
574 ✗ m_quadratureRule->print(verbose);
575 }
576 ✗ std::cout << std::endl;
577 }
578 }
579 }
580