GCC Code Coverage Report

Directory: ./
File: FiniteElement/elementVector.hpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 17 17 100.0%
Branches: 2 4 50.0%
Line Branch Exec Source
1 // ______ ______ _ _ _____ ______
2 // | ____| ____| | (_)/ ____| | ____|
3 // | |__ | |__ | | _| (___ ___| |__
4 // | __| | __| | | | |\___ \ / __| __|
5 // | | | |____| |____| |____) | (__| |____
6 // |_| |______|______|_|_____/ \___|______|
7 // Finite Elements for Life Sciences and Engineering
8 //
9 // License: LGL2.1 License
10 // FELiScE default license: LICENSE in root folder
11 //
12 // Main authors: J. Foulon & J-F. Gerbeau & V. Martin & Miguel Fernandez
13 //
14
15 #ifndef _ELEMVEC_H_INCLUDED
16 #define _ELEMVEC_H_INCLUDED
17
18 // System includes
19 #include <vector>
20 #include <atomic>
21
22 // External includes
23
24 // Project includes
25 #include "Core/felisce.hpp"
26 #include "Core/shared_pointers.hpp"
27 #include "FiniteElement/elementField.hpp"
28 #include "FiniteElement/curBaseFiniteElement.hpp"
29
30 namespace felisce
31 {
32 ///@name felisce Globals
33 ///@{
34
35 ///@}
36 ///@name Type Definitions
37 ///@{
38
39 ///@}
40 ///@name Enum's
41 ///@{
42
43 ///@}
44 ///@name Functions
45 ///@{
46
47 ///@}
48 ///@name felisce Classes
49 ///@{
50 class ElementVector
51 {
52 public:
53 ///@name Type Definitions
54 ///@{
55
56 /// Pointer definition of ElementVector
57 FELISCE_CLASS_POINTER_DEFINITION(ElementVector);
58
59 ///@}
60 ///@name Life Cycle
61 ///@{
62
63 // default constructor
64 ElementVector() = default;
65
66 //! constructor for an arbitrary number of finite elements
67 ElementVector(const std::vector<const CurBaseFiniteElement*>& fe, const std::vector<std::size_t>& nbr);
68
69 // constructor for 1 finite element
70 FELISCE_DEPRECATED_MESSAGE("WARNING:: This constructor will be removed in the near future. Please use the main constructor which consider constant vector.")
71 ElementVector(const CurBaseFiniteElement& fe1, const std::size_t nbr1);
72
73 // constructor for 2 finite elements
74 FELISCE_DEPRECATED_MESSAGE("WARNING:: This constructor will be removed in the near future. Please use the main constructor which consider constant vector.")
75 ElementVector(const CurBaseFiniteElement& fe1, const std::size_t nbr1,
76 const CurBaseFiniteElement& fe2, const std::size_t nbr2);
77
78 // constructor for 3 finite elements
79 FELISCE_DEPRECATED_MESSAGE("WARNING:: This constructor will be removed in the near future. Please use the main constructor which consider constant vector.")
80 ElementVector(const CurBaseFiniteElement& fe1, const std::size_t nbr1,
81 const CurBaseFiniteElement& fe2, const std::size_t nbr2,
82 const CurBaseFiniteElement& fe3, const std::size_t nbr3);
83
84
85 // default destructor
86 324810 ~ElementVector() = default;
87
88 ///@}
89 ///@name Operators
90 ///@{
91
92 /*!
93 * \brief operator +=
94 *
95 * This operator acts upon the underlying UBlasVector and is is assumed all the other attributes are the same (it
96 * is checked by several assert in debug mode).
97 */
98 void operator+=(const ElementVector& vec);
99
100 /*!
101 * \brief operator -=
102 *
103 * This operator acts upon the underlying UBlasVector and is is assumed all the other attributes are the same (it
104 * is checked by several assert in debug mode).
105 */
106 void operator-=(const ElementVector& vec);
107
108 //! Operator multiplication by a scalar.
109 void operator*=(const double factor);
110
111 /**
112 * @brief Access operators
113 */
114 const double& operator () (const std::size_t i) const {
115 return m_vec[i];
116 }
117
118 const double& operator [] (const std::size_t i) const {
119 return m_vec[i];
120 }
121
122 double& operator () (const std::size_t i) {
123 return m_vec[i];
124 }
125
126 532 double& operator [] (const std::size_t i) {
127 532 return m_vec[i];
128 }
129
130 ///@}
131 ///@name Operations
132 ///@{
133
134 76637977 UBlasVectorRange vecBlock(const int i) {
135
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 76637977 times.
 76637977 FEL_ASSERT(i<static_cast<int>(numBlockRow()));
136
1/2
✓ Branch 5 taken 76637977 times.
✗ Branch 6 not taken.
 153275954 return UBlasVectorRange(m_vec,UBlasRange( m_firstRow[i], m_firstRow[i]+m_numRow[i]));
137 }
138
139 /**
140 * @brief Size
141 */
142 357 std::size_t size() const {
143 357 return m_vec.size();
144 }
145
146 /**
147 * @brief Sets to zero the main vector
148 */
149 28185015 void zero() {
150 28185015 m_vec.clear();
151 28185015 }
152
153 void source(const double coef,const CurBaseFiniteElement& fe,const ElementField& elField,const int iblock,const int numComp);
154
155 void source_lumped(const double coef, const CurvilinearFiniteElement& fe,const ElementField& elfield,const int iblock,const int numComp);
156
157 void sourceForComp(const double coef, const CurvilinearFiniteElement& fe,const ElementField& elfield,const int iblock,const int iComp);
158
159 void grad_f_phi_i(const double coef, const CurrentFiniteElement& fe1,const CurrentFiniteElement& fe2,const ElementField& elfield2,const int iblock=0);
160
161 void curl_f_phi_i(const double coef, const CurrentFiniteElement& fe1,const CurrentFiniteElement& fe2,const ElementField& elfield2,const int iblock=0);
162
163 void WGX_dot_phi_i(const CurrentFiniteElement& feWGX, const CurrentFiniteElement& feWSS, const ElementField& WSSField, const int iblock);
164
165 void WGY_dot_phi_i(const CurrentFiniteElement& feWGY, const CurrentFiniteElement& feWSS, const ElementField& WSSField, const int iblock);
166
167 void WGZ_dot_phi_i(const CurrentFiniteElement& feWGZ, const CurrentFiniteElement& feWSS, const ElementField& WSSField, const int iblock);
168
169 void stressX_dot_phi_i(const double viscosity, const CurrentFiniteElement& fePres, const CurrentFiniteElement& feStress,const CurrentFiniteElement& feVel,const ElementField& elfield_pressure,const ElementField& elfield_velocity,const int iblock=0);
170
171 void stressY_dot_phi_i(const double viscosity, const CurrentFiniteElement& fePres, const CurrentFiniteElement& feStress,const CurrentFiniteElement& feVel,const ElementField& elfield_pressure,const ElementField& elfield_velocity,const int iblock=0);
172
173 void stressZ_dot_phi_i(const double viscosity, const CurrentFiniteElement& fePres, const CurrentFiniteElement& feStress,const CurrentFiniteElement& feVel,const ElementField& elfield_pressure,const ElementField& elfield_velocity,const int iblock=0);
174
175 void grad_f_grad_phi_i(const double coef, const CurrentFiniteElement& fe,const ElementField& elfield,const int iblock=0);
176
177 void grad_f_grad_phi_i(const double coef,const CurvilinearFiniteElement& fe, const ElementField& elfield, const std::size_t iblock,const std::size_t num);
178
179 void f_grad_phi_i(const double coef, const CurrentFiniteElement& fe,const ElementField& elfield,const int iblock=0);
180
181 void div_f_phi_i(const double coef, const CurrentFiniteElement& fe1,const CurrentFiniteElement& fe2,const ElementField& elfield1,const int iblock=0);
182
183 void f_div_phi_i(const double coef, const CurrentFiniteElement& fe1,const CurrentFiniteElement& fe2,const ElementField& elfield2,const int iblock=0);
184
185 // \int_Gamma f v \cdot n (To compute /intGamma p v.n for instance) Attention ! is this the same of: f_n_dot_phi_i ? (different only for DOf case)
186 //----------------------------------------------------------------------
187 void f_phi_i_scalar_n(const double coef,const CurvilinearFiniteElement& fe,const ElementField& elfield,const int iblock=0);
188
189 // \int_Gamma f v \cdot J F^-T n (to evaluate follower pressure residual)
190 //----------------------------------------------------------------------
191 void f_phi_i_scalar_piola_n(const double coef,const CurvilinearFiniteElement& fe,const ElementField& elf,const ElementField& elfd,const int iblock=0);
192
193 /// contact
194 // /int_Gamma [disp \cdot n - gap ]_+ v \cdot n
195 //----------------------------------------------------------------------
196 void penalty_contact_phi_i_scalar_n(const double coef,const std::vector<double>& normal, const CurvilinearFiniteElement& fe,const ElementField& disp, const ElementField& gap);
197
198 /// penalization valves
199 // int_Gamma [u^(n,k-1) cdot n]_+ v
200 //----------------------------------------------------------------------
201 void penalty_phi_i_scalar_n(const double coef, CurvilinearFiniteElement& fe, const ElementField& vel, std::size_t iblock);
202
203 /// contact
204 // /int_Gamma [disp \cdot n - gap ]_+ v \cdot n
205 //----------------------------------------------------------------------
206 void nitsche_contact_phi_i_scalar_n(const double coef,const std::vector<double>& normal,
207 const CurvilinearFiniteElement& fe,const ElementField& disp, const ElementField& gap, const ElementField& lambda);
208
209 /// contact
210 // /int_Gamma [disp \cdot n - gap ]_+ v \cdot n
211 //----------------------------------------------------------------------
212 void nitsche_contact_phi_i_scalar_n(const double coef,
213 const CurvilinearFiniteElement& fe,const ElementField& disp, const ElementField& gap, const ElementField& normal, const ElementField& lambda);
214
215 // \int_Gamma [ p^n v \cdot n ] Attention ! is this the same of: f_phi_i_scalar_n ? (different only for DOf case)
216 //----------------------------------------------------------------------
217 void f_n_dot_phi_i(const double coef,const CurvilinearFiniteElement& fe1,const CurvilinearFiniteElement& fe2, const ElementField& elfield,const int iblock);
218
219 //----------------------------------------------------------------------
220 // ex. int_Sigma [ u^{n+1} * N * q ]
221 void f_dot_n_phi_i(const double coef,const CurvilinearFiniteElement& fe1,const CurvilinearFiniteElement& fe2, const ElementField& elfield,const int iblock);
222
223 // For Nitsche-XFEM formulation
224 void g_dot_n_f_phi_i(const double coef, const CurrentFiniteElement& fe, const ElementField& elt_f, const ElementField& elt_g, const ElementField& elt_n, felInt iblock, felInt numComp);
225
226 //----------------------------------------------------------------------
227 // operator \int_{Omega} Grad u^{n} dot Grad LS^{n} / norm(Grad LS^{n}) v (LS^{n} a level std::set of u^{n})
228 void grad_u_dot_grad_LSu_v(const double coef,const ElementField& u,const ElementField& LSu,const std::vector<double>& data,double avHeavU0,double avHeavMinusU0,double insidePar,double outsidePar,const CurvilinearFiniteElement& fe,std::size_t iblock);
229
230 /**
231 *Operator \f[\int_\Sigma \nabla v_1 \cdot v_2 \cdot \phi_i\f]
232 */
233 void grad_vec2_dot_vec1_dot_phi_i(const double coef, const ElementField& vec1, const ElementField& vec2, const CurrentFiniteElement& fe, const int iblock, const int numComp);
234
235 //----------------------------------------------------------------------
236 // Operator \int_\Sigma \eps vec2 \dot vec1 \dot \phi_i
237 void eps_vec2_dot_vec1_dot_phi_i(const double coef, const ElementField& vec1, const ElementField& vec2, const CurrentFiniteElement& fe, const int iblock, const int numComp);
238
239 //----------------------------------------------------------------------
240 // Operator \int_\Sigma (\grad \phi_i \dot vec1) \dot vec2
241 void grad_phi_i_dot_vec1_dot_vec2(const double coef, const ElementField& vec1, const ElementField& vec2, const CurrentFiniteElement& fe, const int iblock, const int numComp);
242
243 //----------------------------------------------------------------------
244 // Operator \int_\Sigma (\eps \phi_i \dot vec1) \dot vec2
245 void eps_phi_i_dot_vec1_dot_vec2(const double coef, const ElementField& vec1, const ElementField& vec2, const CurrentFiniteElement& fe, const int iblock, const int numComp);
246
247 //----------------------------------------------------------------------
248 // Operator \int_\Sigma (\eps vec2 vec1) \dot \esp \phi_i vec1
249 void eps_vec2_vec1_dot_eps_phi_i_vec1(const double coef, const ElementField& vec1, const ElementField& vec2, const CurrentFiniteElement& fe, const int iblock, const int numComp);
250
251 /**
252 * Operator \f[\int_\Sigma \nabla vec_2 vec_1 \cdot \nabla \phi_i vec1\f]
253 */
254 void grad_vec2_vec1_dot_grad_phi_i_vec1(const double coef, const ElementField& vec1, const ElementField& vec2, const CurrentFiniteElement& fe, const int iblock, const int numComp);
255
256 //----------------------------------------------------------------------
257 // Operators \f$ \int_{\Sigma} (\grad \phi_i \dot n) f \f$,
258 //! where $\Sigma$ is a boundary element. (For Nitsche's treatment...)
259 void grad_phi_i_dot_n_f(const double coef,const CurrentFiniteElementWithBd& fewbd,const ElementField& elfield,
260 const int iblockBd, const int numblockBd,const int iblock, const int numComp);
261
262 //----------------------------------------------------------------------
263 // Operators \f$ \int_{\Sigma} (\grad f \dot n) \phi_i \f$,
264 //! where $\Sigma$ is a boundary element. (For Nitsche's treatment...)
265 void grad_f_dot_n_phi_i(const double coef,const CurrentFiniteElementWithBd& fewbd, const ElementField& elfield,
266 const int iblockBd, const int numblockBd, const int iblock, const int numComp);
267 void grad_f_Transp_dot_n_phi_i(const double coef,const CurrentFiniteElementWithBd& fewbd,
268 const ElementField& elfield,
269 const int iblockBd, const int numblockBd,
270 const int iblock, const int numComp) ; // not tested yet, saverio 8/10/2012
271
272 //operator \int_{\Sigma} [ f^{n-1}\cdot n ] f^{n}\phi_i
273 void f_dot_n_f_phi_i(const double coef,const CurvilinearFiniteElement& fe,
274 const ElementField& elfield1,const ElementField& elfield2,
275 const int iblock, const int numComp);
276
277 //ATTENTION RANGEMENT
278 //computes std::vector with component int(element fe)(coef*div(elfield)*div(phi_i)) at position i
279 void div_u_div_phi_i(const double coef,const CurrentFiniteElement& fe, const ElementField& elfield, const int iblock);
280
281 //computes vector with component int(element fe)(elfield.phi_i) at position i
282 void f_scalar_phi_i(const double coef,const CurrentFiniteElement& fe, const ElementField& elfield,const int iblock);
283
284 //! Hyperelasticity case
285 void GradientBasedHyperElasticityVector_grad_phi_i(const double coef,
286 const CurrentFiniteElement& feDisplacement,
287 int quadraturePointIndex,
288 UBlasVector& rhsPart,
289 std::size_t iblock);
290
291 void sourceHyp1DCont(std::vector<double>& coef,const CurrentFiniteElement& fe1,const CurrentFiniteElement& fe2,
292 const ElementField& elfield1,const ElementField& elfield2,const int iblock=0,int iterm=0);
293
294 void sourceHyp1DMom(std::vector<double>& coef,const CurrentFiniteElement& fe1,const CurrentFiniteElement& fe2,
295 const ElementField& elfield1,const ElementField& elfield2,const int iblock=0,int iterm=0);
296
297 //----------------------------------------------------------------------
298 // term
299 // int_Sigma 1/2(|| w.n || - w.n ) (w.phi)
300 void abs_u_dot_n_u_dot_phi_i(const double coef,const ElementField& vel,const CurvilinearFiniteElement& feBd,
301 const int iblock, const int numComp);
302
303 void u_grad_u_phi_i(const double coef,const CurrentFiniteElement& fe,const ElementField& elField,const int iblock,const int numComp);
304
305 void u_grad_u_grad_phi_i(const double coef,const CurrentFiniteElement& fe,const ElementField& elField,const int iblock,const int numComp);
306
307 // /int_Gamma ( f \cdot n ) ( v \cdot n )
308 void f_dot_n_phi_i_dot_n(const double coef,const CurvilinearFiniteElement& fe,const ElementField& elfield,const int iblock);
309
310 // /int_Gamma ( f \times n ) ( v \times n )
311 //----------------------------------------------------------------------
312 void f_times_n_phi_i_times_n(const double coef,const CurvilinearFiniteElement& fe,const ElementField& elfield,const int iblock);
313
314 // /int_Gamma ( \Grad f n \times n ) ( v \times n )
315 //----------------------------------------------------------------------
316 void grad_f_dot_n_times_n_phi_i_times_n(const double coef,const CurrentFiniteElementWithBd& fewbd, const ElementField& elfield,
317 const int iblockBd, const int numblockBd, const int iblock, const int numComp);
318
319 void stab_supgAdvDiffCT(const double dt, const double stab1, const double density, const double visc,
320 const CurrentFiniteElement& feVel, const CurrentFiniteElement& fePres,
321 const ElementField& vel, const ElementField& pres, const ElementField& rhs);
322
323 void stab_supgProjCT(const double dt, const double stab1, const double density, const double visc, const CurrentFiniteElement& fe,
324 const ElementField& beta, const ElementField& vel, const ElementField& pres,
325 const ElementField& rhs);
326
327 void stab_supgRhsFS(const double dt, const double stabSUPG, const double density, const double viscosity,
328 const double c0, const double c1, const double c2,
329 const CurrentFiniteElement& feVel, const CurrentFiniteElement& fePre,
330 const ElementField& vel, const ElementField& pres, const ElementField& velRHS);
331
332 //----------------------------------------------------------------------
333 // qudratic term: alpha(x) T^2 (for scalars!)
334 void quadratic(const double coef, const CurrentFiniteElement& fe,const ElementField& elfield, const int iblock,const int numComp);
335
336 void sGrad_phi_i_tensor_sGrad_f(const double coef, const ElementField& f, const ElementField& normalField, const CurvilinearFiniteElement& fe, UBlasMatrix tensor);
337 void sGrad_phi_i_tensor_sGrad_f(const double coef, const ElementField& elemCoef, const ElementField& f, const ElementField& normalField, const CurvilinearFiniteElement& fe, UBlasMatrix tensor);
338
339 //---------------------------------------------------------------------------------------------------------------
340 // Some explanation about SUPG:
341 // noting T:time, L:length, M:mass
342 // WARNING: tau has a dimension of T*L^3*M^{-1} -> remember this when you use P1-P1
343 // It's advised to choose your units knowing this.
344 //---------------------------------------------------------------------------------------------------------------
345 void stab_supg(const double dt, const double stab1, const double density, const double visc,
346 const CurrentFiniteElement& fe,const ElementField& vel,const ElementField& rhs,
347 double& Re_elem,double& tau, const int type);
348
349 void transpirationSurfStabRHS( const double coef, const ElementField& elemCoef, const ElementField& normalField, const ElementField& force, const CurvilinearFiniteElement& fe, UBlasMatrix tensor);
350
351 // computes non-linear residual associated to \int_G(d) p n_K \cdot v for thin-walled solids
352 void followerPressure(const double coef, const ElementField& press, const ElementField& disp, const CurvilinearFiniteElement& fe);
353
354 void velocityConstraint(const double coef, const ElementField& disp, const CurvilinearFiniteElement& fe, const int iblock);
355
356 ///@}
357 ///@name Access
358 ///@{
359
360 675983746 UBlasVector& vec() {
361 675983746 return m_vec;
362 }
363
364 2 const UBlasVector& vec() const {
365 2 return m_vec;
366 }
367
368 76637985 std::size_t numBlockRow() const {
369 76637985 return m_numRow.size();
370 }
371
372 ///@}
373 ///@name Inquiry
374 ///@{
375
376 ///@}
377 ///@name Input and output
378 ///@{
379
380 void print(int verbose,std::ostream& c=std::cout);
381
382 ///@}
383 ///@name Friends
384 ///@{
385
386 ///@}
387 protected:
388 ///@name Protected static Member Variables
389 ///@{
390
391 ///@}
392 ///@name Protected member Variables
393 ///@{
394
395 ///@}
396 ///@name Protected Operators
397 ///@{
398
399 ///@}
400 ///@name Protected Operations
401 ///@{
402
403 ///@}
404 ///@name Protected Access
405 ///@{
406
407 ///@}
408 ///@name Protected Inquiry
409 ///@{
410
411 ///@}
412 ///@name Protected LifeCycle
413 ///@{
414
415 ///@}
416 private:
417 ///@name Private static Member Variables
418 ///@{
419
420 ///@}
421 ///@name Private member Variables
422 ///@{
423
424 UBlasVector m_vec; //!< The element vectorrows
425 std::vector<std::size_t> m_numRow; // m_numRow[i]=nb of rows in the i-th block row
426 std::vector<std::size_t> m_firstRow; //m_firstRow[i]=index of first row of i-th block
427
428 // TODO:: All temp classes should be replaced with on time execution (otherwise much memory consumed)
429 UBlasVector m_tmpVecDof;//!< temporary vector (size = numDof)
430 UBlasVector m_tmpVecComp;//!< temporary vector (size = numComp)
431 UBlasVector m_tmpVecCoor;//!< temporary vector (size = numCoor)
432 UBlasVector m_tmpVec2Coor;//!< temporary vector (size = numCoor)
433 UBlasMatrix m_tmpMatDof;
434
435 ///@}
436 ///@name Private Operators
437 ///@{
438
439 ///@}
440 ///@name Private Operations
441 ///@{
442
443 ///@}
444 ///@name Private Access
445 ///@{
446
447 ///@}
448 ///@name Private Inquiry
449 ///@{
450
451 ///@}
452 ///@name Private LifeCycle
453 ///@{
454
455 ///@}
456 };
457 ///@}
458 ///@name Type Definitions
459 ///@{
460
461 ///@}
462 } /* namespace felisce.*/
463
464 #endif
465