GCC Code Coverage Report

Directory: ./
File: FiniteElement/elementField.hpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 54 71 76.1%
Branches: 23 39 59.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-F. Gerbeau
13 //
14
15 #ifndef _ELEMFIELD_H_INCLUDED
16 #define _ELEMFIELD_H_INCLUDED
17
18 // System includes
19 #include <vector>
20
21 // External includes
22 #include "Core/NoThirdPartyWarning/Petsc/ao.hpp"
23 #include "Core/NoThirdPartyWarning/Petsc/vec.hpp"
24
25 // Project includes
26 #include "PETScInterface/petscVector.hpp"
27 #include "Core/felisce.hpp"
28 #include "FiniteElement/currentFiniteElement.hpp"
29 #include "FiniteElement/currentFiniteElementWithBd.hpp"
30 #include "FiniteElement/curvilinearFiniteElement.hpp"
31 #include "DegreeOfFreedom/variable.hpp"
32 #include "DegreeOfFreedom/dof.hpp"
33 #include "Geometry/geometricMeshRegion.hpp"
34
35 namespace felisce
36 {
37 ///@name felisce Globals
38 ///@{
39
40 ///@}
41 ///@name Type Definitions
42 ///@{
43
44 ///@}
45 ///@name Enum's
46 ///@{
47
48 ///@}
49 ///@name Functions
50 ///@{
51
52 ///@}
53 ///@name felisce Classes
54 ///@{
55
56 /*!
57 \class ElementField
58 \authors J-F. Gerbeau
59
60 \brief Class handling the values of fields at the element level
61
62 An ElementField contains the values of a scalar field, e.g. a diffusion coefficient or a source term, at the element level.
63 It can be either constant, defined on the DOF of the finite element or defined on the quadrature points.
64
65 let \f$(c_i)_{i=1..n}$ denote the values of the field, \f$n$ being the member m_dim.
66
67 - Type CONSTANT_FIELD: m_dim = 1, the value of the field is just a constant.
68
69 - Type DOF_FIELD: m_dim = number of DOF of the element. In the element operators the values on the quadrature points
70 will be computed according to \f$\sum_{i=1}^{n} c_i \phi_i(\hat x_{ig}) \f$, where $\f \hat x_{ig} \f$ denote the quadrature points
71 coordinates
72
73 - Type QUAD_POINT_FIELD: m_dim = number of quadrature points of the element. In the element operators the field
74 values will used directly in the quadrature points.
75
76 A ElementField can be initialized from a constant, a function of (x,y,z), a std::vector.
77
78 \warning Remind to call the setValue function of the ElementField for every new element,
79 after having updated the current finite element (when the value of the ElementField is element dependent)
80
81 \todo ElementVectorField and ElementTensorField
82 */
83 class ElementField
84 {
85 public:
86
87 ///@name Type Definitions
88 ///@{
89
90 ///@}
91 ///@name Life Cycle
92 ///@{
93
94 /// constructors
95 ElementField();
96
97 // Destructor
98 ~ElementField();
99
100 /// Copy constructor, to avoid warning due to user-declared destructor.
101 ✗ ElementField(const ElementField&) = default;
102
103 ///@}
104 ///@name Operations
105 ///@{
106
107 ///@}
108 ///@name Operators
109 ///@{
110
111 //intialize methods to suppress pointer declaration, 20/11/2011 Jeremie
112 void initialize(const ElementFieldType type=CONSTANT_FIELD, const int numComp=1);
113 void initialize(const ElementFieldType type,const CurrentFiniteElement& fe,const int numComp=1);
114 void initialize(const ElementFieldType type,const CurrentFiniteElementWithBd& fewbd,const int numComp=1);
115 void initialize(const ElementFieldType type,const CurvilinearFiniteElement& fec,const int numComp=1);
116
117 14 void update(const double value) {
118
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 14 times.
 14 FEL_ASSERT(m_type==CONSTANT_FIELD);
119 14 val(0,0) = value;
120 14 }
121
122 /// set value by functions:
123 template <class FunctorXYZ>
124 inline void setValue(const FunctorXYZ& functorXYZ, CurrentFiniteElement& fe, const int icomp=0);
125 template <class FunctorXYZ>
126 inline void setValueVec(const FunctorXYZ& functorXYZ, CurrentFiniteElement& fe, const int icomp);
127
128 template <class FunctorXYZ>
129 inline void setValue(const FunctorXYZ& functorXYZ, CurrentFiniteElementWithBd& fewbd, const int icomp=0);
130
131 template <class FunctorXYZ>
132 inline void setValue(const FunctorXYZ& functorXYZ, CurvilinearFiniteElement& fe, const int icomp=0);
133
134 template <class FunctorXYZ>
135 inline void setValueVec(const FunctorXYZ& functorXYZ, CurvilinearFiniteElement& fe, const int icomp);
136
137 template <class FunctorXYZT>
138 inline void setValueVec(const FunctorXYZT& functorXYZT, CurvilinearFiniteElement& fe, const double time, const int icomp);
139
140 template <class FunctorXYZT>
141 inline void setValue(const FunctorXYZT& functorXYZT, CurrentFiniteElement& fe, const double time,const int icomp=0);
142
143 template <class FunctorXYZT>
144 inline void setValue(const FunctorXYZT& functorXYZT, CurrentFiniteElementWithBd& fewbd, const double time, const int icomp=0);
145
146 template <class FunctorXYZT>
147 inline void setValue(const FunctorXYZT& functorXYZT, CurvilinearFiniteElement& fe, const double time, const int icomp=0);
148
149 template <class FunctorT> //Attention! this function it's hidden by the above one because of the default value of icomp, I don't know if it is possible to call it. (matteo feb 2014)
150 inline void setValueT(const FunctorT& functorT, CurvilinearFiniteElement& fe, const double time);
151
152 /// set value by vectors:
153
154 void setValue(const std::vector<double>& value);
155
156 //---any type of field, component-wise constant. (think for instance about gravity and m_elemFieldRHS)
157 void setValue(double value, felInt icomp);
158
159 // set value (all components)
160 // using AO application to the numbering
161 void setValue(const PetscVector& v, const CurBaseFiniteElement& fe, const felInt iel, const int idVariable, const AO& ao, const Dof& dof);
162
163 // using ISLocalToGlobalMapping application to the numbering (to BDDC Thin Structures)
164 // this fonction is the same that the previous setValue function exect that it work with an ISLocalToGlobalMapping and not a AO mapping (to BDDC Thin Structures)
165 //----------------------------------------------------------------------
166 void setValue(const PetscVector& v,const CurvilinearFiniteElement& fe, const felInt iel, const int idVariable,ISLocalToGlobalMapping matisMapDof, const Dof& dof);
167
168 // set value (select a specific component)
169 void setValue(const PetscVector& v, const CurvilinearFiniteElement& fe, const felInt iel, const int idVariable, const AO& ao, const Dof& dof, const int iComp);
170
171 // set value using a specific component of the petsc std::vector and setting the result in
172 // a specific component of the element field
173 // this is used to set each component of the elementfield using a different petscvector
174 void setValue(const PetscVector& v, const CurvilinearFiniteElement& fe, const felInt iel, const int idVariable, const AO& ao, const Dof& dof, const int iComp, const int efComp);
175
176 void setValue(const PetscVector& v, const CurBaseFiniteElement& fe, const felInt iel, const int idVariable, const Dof& dof);
177
178 void setValue(const PetscVector& v, const CurrentFiniteElementWithBd& fe, felInt ibd, const felInt iel, const int idVariable, const AO& ao, const Dof& dof);
179
180 // The idea of this function is to take a petscVec with more than one component, let's say the displacement, compute the scalar product
181 // with the normal in the element (we assume it is constant for now)
182 // and then save the result in a scalar element fields
183 void setValueFdotNormal(const PetscVector& v,const CurvilinearFiniteElement& fe, const felInt iel, const int idVariable, const AO& ao, const Dof& dof);
184
185 void setValueMatching(double* DofValues, felInt* DofIdentities, felInt DofDimension, const CurvilinearFiniteElement& fe, const felInt iel, const int idVariable, const Dof& dof);
186
187 ///set value by ElementFieldDynamicValue read in data file
188 void setValue(CurrentFiniteElement fe, ElementFieldDynamicValue &elementFieldDynamicValue, const double time);
189
190 ///@}
191 ///@name Access
192 ///@{
193
194 88215374 inline ElementFieldType type() const {
195 88215374 return m_type;
196 }
197
198 4112 inline int numComp() const {
199 4112 return m_numComp;
200 }
201
202 /// Return the type of finite element
203 1711 ElementFieldType getType() const {
204 1711 return m_type;
205 }
206
207 UBlasMatrix& get_val();
208
209 // There probably is a UBlas function to do this, but I was not able to find it.
210 UBlasVector valAsVec() const;
211
212 ///@}
213 ///@name Public member Variables
214 ///@{
215
216 /*!
217 val is the matrix containing the values of the field. It has numComp rows (e.g. numComp=1 for a scalar field,
218 or numComp=number of components of a std::vector field). The number of columns is 1 if type=CONSTANT_FIELD, nDof if type=DOF_FIELD,
219 and nQuadPoint if type=QUAD_POINT_FIELD
220 */
221 UBlasMatrix val; // val(icomp,j)
222 UBlasMatrix valTmp;
223
224 ///@}
225 ///@name Inquiry
226 ///@{
227
228 ///@}
229 ///@name Input and output
230 ///@{
231
232 void print(int verbose,std::ostream& c=std::cout);
233
234 ///@}
235 protected:
236 ///@name Protected static Member Variables
237 ///@{
238
239 ///@}
240 ///@name Protected member Variables
241 ///@{
242
243 ///@}
244 ///@name Protected Operators
245 ///@{
246
247 ///@}
248 ///@name Protected Operations
249 ///@{
250
251 ///@}
252 ///@name Protected Access
253 ///@{
254
255 ///@}
256 ///@name Protected Inquiry
257 ///@{
258
259 ///@}
260 ///@name Protected LifeCycle
261 ///@{
262
263 ///@}
264 private:
265 ///@name Private static Member Variables
266 ///@{
267
268 ///@}
269 ///@name Private member Variables
270 ///@{
271 ElementFieldType m_type;
272 int m_numComp;
273 int m_dim;
274
275 ///Arrays use with Petsc vectors.
276 ///Use global ordering of felisce and find global ordering for petsc.
277 int* m_findDofPositionInPetscOrdering;
278
279 ///@}
280 ///@name Private Operators
281 ///@{
282
283 ///@}
284 ///@name Private Operations
285 ///@{
286
287 void setShape(CurrentFiniteElement fe, const ElementFieldType type, const int numComp);
288
289 ///@}
290 ///@name Private Operations
291 ///@{
292
293 ///@}
294 ///@name Private Access
295 ///@{
296
297 ///@}
298 ///@name Private Inquiry
299 ///@{
300
301 ///@}
302 ///@name Private LifeCycle
303 ///@{
304
305 ///@}
306
307 }; // class
308
309 ///@}
310 ///@name Type Definitions
311 ///@{
312
313 //----------------------------------------------------------------------
314 // Implementation of functors' setValue:
315 //----------------------------------------------------------------------
316 template <class FunctorXYZ>
317 418068 inline void ElementField::setValue(const FunctorXYZ& functorXYZ, CurrentFiniteElement& fe,const int icomp) {
318
1/3
✓ Branch 0 taken 418068 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 418068 switch (m_type) {
319 418068 case QUAD_POINT_FIELD:
320
2/2
✓ Branch 1 taken 139356 times.
✓ Branch 2 taken 278712 times.
 418068 if ( !fe.hasQuadPtCoor() ) {
321 139356 fe.computeCurrentQuadraturePoint();
322 }
323
2/2
✓ Branch 1 taken 6271020 times.
✓ Branch 2 taken 418068 times.
 6689088 for(int ig=0; ig<fe.numQuadraturePoint(); ig++) {
324 6271020 val(icomp,ig) = functorXYZ(fe.currentQuadPoint[ig].x(),fe.currentQuadPoint[ig].y(),fe.currentQuadPoint[ig].z());
325 }
326 418068 break;
327 ✗ case CONSTANT_FIELD:
328 case DOF_FIELD:
329 ✗ FEL_ERROR("ElementField::setValue with a functor is not yet implemented for this type of element field");
330 ✗ break;
331 }
332 418068 }
333
334 template <class FunctorXYZ>
335 88438 inline void ElementField::setValueVec(const FunctorXYZ& functorXYZ, CurrentFiniteElement& fe,const int icomp) {
336
1/3
✓ Branch 0 taken 88438 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 88438 switch (m_type) {
337 88438 case QUAD_POINT_FIELD:
338
2/2
✓ Branch 1 taken 44149 times.
✓ Branch 2 taken 44289 times.
 88438 if ( !fe.hasQuadPtCoor() ) {
339 44149 fe.computeCurrentQuadraturePoint();
340 }
341
2/2
✓ Branch 1 taken 265816 times.
✓ Branch 2 taken 88438 times.
 354254 for(int ig=0; ig<fe.numQuadraturePoint(); ig++) {
342 265816 val(icomp,ig) = functorXYZ(icomp,fe.currentQuadPoint[ig].x(),fe.currentQuadPoint[ig].y(),fe.currentQuadPoint[ig].z());
343 }
344 88438 break;
345 ✗ case CONSTANT_FIELD:
346 case DOF_FIELD:
347 ✗ FEL_ERROR("ElementField::setValue with a functor is not yet implemented for this type of element field");
348 ✗ break;
349 }
350 88438 }
351
352
353 //----------------------------------------------------------------------
354
355 template <class FunctorXYZ>
356 inline void ElementField::setValue(const FunctorXYZ& functorXYZ, CurrentFiniteElementWithBd& fewbd,const int icomp) {
357 switch (m_type) {
358 case QUAD_POINT_FIELD:
359 if ( !fewbd.hasQuadPtCoor() ) {
360 fewbd.computeCurrentQuadraturePointInternAndBd();
361 }
362 for(int ig=0; ig<fewbd.numQuadraturePointInternAndBd(); ig++) {
363 val(icomp,ig) = functorXYZ(fewbd.currentQuadPoint[ig].x(),fewbd.currentQuadPoint[ig].y(),fewbd.currentQuadPoint[ig].z());
364 }
365 break;
366 case CONSTANT_FIELD:
367 case DOF_FIELD:
368 FEL_ERROR("ElementField::setValue with a functor is not yet implemented for this type of element field");
369 break;
370 }
371 }
372
373 //----------------------------------------------------------------------
374 template <class FunctorXYZ>
375 2400 inline void ElementField::setValue(const FunctorXYZ& functorXYZ, CurvilinearFiniteElement& fe,const int icomp) {
376
1/3
✓ Branch 0 taken 2400 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 2400 switch (m_type) {
377 2400 case QUAD_POINT_FIELD:
378
1/2
✓ Branch 1 taken 2400 times.
✗ Branch 2 not taken.
 2400 if ( !fe.hasQuadPtCoor() ) {
379 2400 fe.computeCurrentQuadraturePoint();
380 }
381
2/2
✓ Branch 1 taken 4800 times.
✓ Branch 2 taken 2400 times.
 7200 for(int ig=0; ig<fe.numQuadraturePoint(); ig++) {
382 4800 val(icomp,ig) = functorXYZ(fe.currentQuadPoint[ig].x(),fe.currentQuadPoint[ig].y(),fe.currentQuadPoint[ig].z());
383 }
384 2400 break;
385 ✗ case CONSTANT_FIELD:
386 case DOF_FIELD:
387 ✗ FEL_ERROR("ElementField::setValue with a functor is not yet implemented for this type of element field");
388 ✗ break;
389 }
390 2400 }
391
392 template <class FunctorXYZ>
393 2400 inline void ElementField::setValueVec(const FunctorXYZ& functorXYZ, CurvilinearFiniteElement& fe,const int icomp) {
394
1/3
✓ Branch 0 taken 2400 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 2400 switch (m_type) {
395 2400 case QUAD_POINT_FIELD:
396
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2400 times.
 2400 if ( !fe.hasQuadPtCoor() ) {
397 ✗ fe.computeCurrentQuadraturePoint();
398 }
399
2/2
✓ Branch 1 taken 4800 times.
✓ Branch 2 taken 2400 times.
 7200 for(int ig=0; ig<fe.numQuadraturePoint(); ig++) {
400 4800 val(icomp,ig) = functorXYZ(icomp,fe.currentQuadPoint[ig].x(),fe.currentQuadPoint[ig].y(),fe.currentQuadPoint[ig].z());
401 }
402 2400 break;
403 ✗ case CONSTANT_FIELD:
404 case DOF_FIELD:
405 ✗ FEL_ERROR("ElementField::setValueVec with a functor is not yet implemented for this type of element field");
406 ✗ break;
407 }
408 2400 }
409
410 template <class FunctorXYZT>
411 inline void ElementField::setValueVec(const FunctorXYZT& functorXYZT, CurvilinearFiniteElement& fe, const double time, const int icomp) {
412 switch (m_type) {
413 case QUAD_POINT_FIELD:
414 if ( !fe.hasQuadPtCoor() ) {
415 fe.computeCurrentQuadraturePoint();
416 }
417 for(int ig=0; ig<fe.numQuadraturePoint(); ig++) {
418 val(icomp,ig) = functorXYZT(icomp,fe.currentQuadPoint[ig].x(),fe.currentQuadPoint[ig].y(),fe.currentQuadPoint[ig].z(),time);
419 }
420 break;
421 case CONSTANT_FIELD:
422 case DOF_FIELD:
423 FEL_ERROR("ElementField::setValueVec with a functor is not yet implemented for this type of element field");
424 break;
425 }
426 }
427
428 //----------------------------------------------------------------------
429 template <class FunctorXYZT>
430 21708 inline void ElementField::setValue(const FunctorXYZT& functorXYZT, CurrentFiniteElement& fe, const double time,const int icomp) {
431
1/3
✓ Branch 0 taken 19046 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 21708 switch (m_type) {
432 21708 case QUAD_POINT_FIELD:
433
1/2
✓ Branch 1 taken 19046 times.
✗ Branch 2 not taken.
 21708 if ( !fe.hasQuadPtCoor() ) {
434 21708 fe.computeCurrentQuadraturePoint();
435 }
436
2/2
✓ Branch 1 taken 65124 times.
✓ Branch 2 taken 19046 times.
 102804 for(int ig=0; ig<fe.numQuadraturePoint(); ig++) {
437 81096 val(icomp,ig) = functorXYZT(fe.currentQuadPoint[ig].x(),fe.currentQuadPoint[ig].y(),fe.currentQuadPoint[ig].z(),time);
438 }
439 21708 break;
440 ✗ case CONSTANT_FIELD:
441 case DOF_FIELD:
442 ✗ FEL_ERROR("ElementField::setValue with a functor is not yet implemented for this type of element field");
443 ✗ break;
444 }
445 21708 }
446
447 //----------------------------------------------------------------------
448 template <class FunctorXYZT>
449 inline void ElementField::setValue(const FunctorXYZT& functorXYZT, CurrentFiniteElementWithBd& fewbd, const double time,const int icomp) {
450 switch (m_type) {
451 case QUAD_POINT_FIELD:
452 if ( !fewbd.hasQuadPtCoor() ) {
453 fewbd.computeCurrentQuadraturePointInternAndBd();
454 }
455 for(int ig=0; ig<fewbd.numQuadraturePointInternAndBd(); ig++) {
456 val(icomp,ig) = functorXYZT(fewbd.currentQuadPoint[ig].x(),fewbd.currentQuadPoint[ig].y(),fewbd.currentQuadPoint[ig].z(),time);
457 }
458 break;
459 case CONSTANT_FIELD:
460 case DOF_FIELD:
461 FEL_ERROR("ElementField::setValue with a functor is not yet implemented for this type of element field");
462 break;
463 }
464 }
465
466 //----------------------------------------------------------------------
467 template <class FunctorXYZT>
468 inline void ElementField::setValue(const FunctorXYZT& functorXYZT, CurvilinearFiniteElement& fe, const double time,const int icomp) {
469 switch (m_type) {
470 case QUAD_POINT_FIELD:
471 if ( !fe.hasQuadPtCoor() ) {
472 fe.computeCurrentQuadraturePoint();
473 }
474 for(int ig=0; ig<fe.numQuadraturePoint(); ig++) {
475 val(icomp,ig) = functorXYZT(fe.currentQuadPoint[ig].x(),fe.currentQuadPoint[ig].y(),fe.currentQuadPoint[ig].z(),time);
476 }
477 break;
478 case CONSTANT_FIELD:
479 case DOF_FIELD:
480 FEL_ERROR("ElementField::setValue with a functor is not yet implemented for this type of element field");
481 break;
482 }
483 }
484
485 //----------------------------------------------------------------------
486 template <class FunctorT>
487 inline void ElementField::setValueT(const FunctorT& functorT, CurvilinearFiniteElement& fe, const double time) {
488 switch (m_type) {
489 case QUAD_POINT_FIELD:
490 if ( !fe.hasQuadPtCoor() ) {
491 fe.computeCurrentQuadraturePoint();
492 }
493 for(int ig=0; ig<fe.numQuadraturePoint(); ig++)
494 val(0,ig) = functorT(time);
495 break;
496 case CONSTANT_FIELD:
497 case DOF_FIELD:
498 FEL_ERROR("ElementField::setValue with a functor is not yet implemented for this type of element field");
499 break;
500 }
501 }
502
503 //============================================================================================
504 /*
505 class ElementVectorField{
506 public:
507 enum Type {
508 CONSTANT_FIELD,
509 DOF_FIELD,
510 QUAD_POINT_FIELD
511 }
512 private:
513 Type m_type;
514 std::size_t m_numComp;
515 std::size_t m_size;
516 public:
517 UBlasVector val;
518 ElementVectorField(ElementFieldType type,const CurrentFiniteElement& fe,const int numComp);
519 void set(UBlasVector value){
520 FEL_ASSERT(m_type==CONSTANT_FIELD);
521 FEL_ASSERT(m_numComp == value.size());
522 val = value;
523 }
524 inline ElementVectorField::Type type() const {return m_type;}
525 template <class Templ_functor>
526 inline void update(const Templ_functor& func,const CurrentFiniteElement& fe,const int icomp);
527 void print(int verbose,std::ostream& c=std::cout);
528 }
529 // Implementation
530 template <class Templ_functor>
531 inline void ElementVectorField::update(const Templ_functor& func,const CurrentFiniteElement& fe,const int icomp=0){
532 switch (m_type) {
533 case QUAD_POINT_FIELD:
534 if ( !fe.hasQuadPtCoor() ) { fe.computeCurrentQuadraturePoint(); }
535 for(int ig=0;ig<fe.numQuadraturePoint();ig++){
536 val(icomp,ig) = func(fe.currentQuadPoint[ig].x(),fe.currentQuadPoint[ig].y(),fe.currentQuadPoint[ig].z());
537 }
538 break;
539 default:
540 FEL_ERROR("ElementField::update with a functor is not yet implemented for this type of element field");
541 break;
542 }
543 }
544
545 */
546 ///@}
547 } /* namespace felisce.*/
548
549 #endif
550