v0.15.0
Loading...
Searching...
No Matches
UserDataOperators.hpp
Go to the documentation of this file.
1/** \file UserDataOperators.hpp
2 * \brief User data operators for finite element computations
3 *
4 * This file contains a comprehensive collection of user data operators used in
5 * MoFEM finite element computations. These operators provide functionality for
6 * calculating, transforming, and manipulating field values at integration points.
7 *
8 * ## Main Categories of Operators:
9 *
10 * ### Field Value Calculation Operators
11 * - **Scalar Field Operators**: Calculate scalar field values at integration points
12 * - OpCalculateScalarFieldValues: General scalar field value calculation
13 * - OpCalculateScalarFieldValuesFromPetscVecImpl: PETSc vector-based scalar values
14 * - **Vector Field Operators**: Calculate vector field values at integration points
15 * - OpCalculateVectorFieldValues: General vector field value calculation
16 * - OpCalculateDivergenceVectorFieldValues: Divergence calculation for vector fields
17 * - **Tensor Field Operators**: Calculate tensor field values at integration points
18 * - OpCalculateTensor2FieldValues: Second-rank tensor field calculations
19 * - OpCalculateTensor2SymmetricFieldValues: Symmetric tensor field calculations
20 *
21 * ### Field Gradient Calculation Operators
22 * - **Scalar Gradient Operators**:
23 * - OpCalculateScalarFieldGradient: Gradient of scalar fields
24 * - OpCalculateScalarFieldHessian: Hessian (second derivatives) of scalar fields
25 * - **Vector Gradient Operators**:
26 * - OpCalculateVectorFieldGradient: Gradient of vector fields
27 *
28 * ### Matrix and Tensor Operations
29 * - **Matrix Manipulation**:
30 * - OpInvertMatrix: Matrix inversion at integration points
31 * - OpScaleMatrix: Element-wise matrix scaling operations
32 * - OpCalculateTraceFromMat: Trace calculation for matrices
33 * - **Tensor Operations**:
34 * - OpSymmetrizeTensor: Tensor symmetrization operations
35 * - OpTensorTimesSymmetricTensor: Tensor multiplication operations
36 *
37 * ### Coordinate and Geometry Operations
38 * - **Transformation Operators**: Handle coordinate transformations and geometric operations
39 * - **Integration Point Operations**: Manage data at Gauss integration points
40 *
41 * ## Key Features:
42 *
43 * ### Template-Based Design
44 * Most operators are implemented as template structures allowing:
45 * - Flexible data types (double, float, etc.)
46 * - Various tensor dimensions and field dimensions
47 * - Different storage layouts and allocators
48 * - Compile-time optimization for specific use cases
49 *
50 * ### PETSc Integration
51 * Many operators provide seamless integration with PETSc:
52 * - Direct access to PETSc vectors and matrices
53 * - Efficient data exchange between MoFEM and PETSc
54 * - Support for distributed computing environments
55 *
56 * ### Field Space Support
57 * Operators support various finite element field spaces:
58 * - H1 (continuous) spaces for scalar and vector fields
59 * - Hdiv spaces for flux-like vector fields
60 * - Hcurl spaces for curl-conforming vector fields
61 * - L2 spaces for discontinuous fields
62 *
63 * ### Memory Management
64 * Robust memory management through:
65 * - Smart pointer usage (boost::shared_ptr)
66 * - RAII principles for resource management
67 * - Exception safety and error handling
68 *
69 * ## Usage Patterns:
70 *
71 * 1. **Field Evaluation**: Calculate field values at integration points for
72 * post-processing, visualization, or coupling with other physics
73 *
74 * 2. **Gradient Calculations**: Compute spatial derivatives for stress analysis,
75 * heat conduction, fluid dynamics, and other gradient-dependent phenomena
76 *
77 * 3. **Matrix Operations**: Perform linear algebra operations needed in
78 * constitutive relations, coordinate transformations, and iterative solvers
79 *
80 * 4. **Data Extraction**: Extract computed results from finite element solutions
81 * for analysis, visualization, or coupling with external codes
82 *
83 * ## Thread Safety and Performance:
84 * - Operators are designed to be thread-safe when used properly
85 * - Optimized for performance in parallel finite element computations
86 * - Minimal memory allocations during computation loops
87 *
88 * @note These operators are typically used within the MoFEM pipeline system
89 * and are automatically integrated into finite element assembly processes.
90 *
91 * @see ForcesAndSourcesCore::UserDataOperator for base class documentation
92 * @see EntitiesFieldData for field data management
93 * @see PipelineManager for operator pipeline management
94 */
95
96#ifndef __USER_DATA_OPERATORS_HPP__
97 #define __USER_DATA_OPERATORS_HPP__
98
99namespace MoFEM {
100
101/** \name Get values at Gauss pts */
102
103/**@{*/
104
105/** \name Scalar values */
106
107/**@{*/
108
109/** \brief Scalar field values at integration points
110 *
111 */
112
113/**
114 * @brief Operator for calculating scalar field values at integration points
115 *
116 * This template structure calculates scalar field values at integration points
117 * and stores them in a provided data container. It supports different storage
118 * types through template parameters and can optionally work with PETSc vectors.
119 *
120 * @tparam T Data type for the scalar field values (e.g., double, float)
121 * @tparam A Allocator type for the ublas vector storage
122 */
123template <class T, class A>
126
127 /**
128 * @brief Constructor for scalar field values calculation operator
129 *
130 * @param field_name Name of the scalar field to calculate values for
131 * @param data_ptr Shared pointer to ublas vector for storing calculated values
132 * @param zero_type Entity type for zero-level entities (default: MBVERTEX)
133 */
135 const std::string field_name,
136 boost::shared_ptr<ublas::vector<T, A>> data_ptr,
137 const EntityType zero_type = MBVERTEX)
139 dataPtr(data_ptr), zeroType(zero_type) {
140 if (!dataPtr)
141 THROW_MESSAGE("Pointer is not set");
142 }
143
144 /**
145 * @brief Constructor with PETSc vector support
146 *
147 * @param field_name Name of the scalar field to calculate values for
148 * @param data_ptr Shared pointer to ublas vector for storing calculated values
149 * @param data_vec Smart PETSc vector object for additional data storage
150 * @param zero_type Entity type for zero-level entities (default: MBVERTEX)
151 */
153 const std::string field_name,
154 boost::shared_ptr<ublas::vector<T, A>> data_ptr,
155 SmartPetscObj<Vec> data_vec, const EntityType zero_type = MBVERTEX)
157 dataPtr(data_ptr), zeroType(zero_type), dataVec(data_vec) {
158 if (!dataPtr)
159 THROW_MESSAGE("Pointer is not set");
160 }
161
162 /**
163 * \brief calculate values of scalar field at integration points
164 * @param side side entity number
165 * @param type side entity type
166 * @param data entity data
167 * @return error code
168 */
169 MoFEMErrorCode doWork(int side, EntityType type,
171
172protected:
173 boost::shared_ptr<ublas::vector<T, A>> dataPtr;
177};
178
179/**
180 * \brief Specialization of member function
181 *
182 */
183template <class T, class A>
185 int side, EntityType type, EntitiesFieldData::EntData &data) {
187 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented for T = %s",
188 typeid(T).name() // boost::core::demangle(typeid(T).name()).c_str()
189 );
191}
192
193/**
194 * @brief Specialization for double precision scalar field values calculation
195 *
196 * This structure is a specialization of OpCalculateScalarFieldValues_General
197 * for double precision scalar fields using DoubleAllocator. It provides
198 * concrete implementation for calculating scalar field values at integration
199 * points and storing them in double precision format.
200 *
201 * @ingroup mofem_forces_and_sources_user_data_operators
202 */
204 : public OpCalculateScalarFieldValues_General<double, DoubleAllocator> {
205
207 double, DoubleAllocator>::OpCalculateScalarFieldValues_General;
208
209 /**
210 * \brief calculate values of scalar field at integration points
211 * @param side side entity number
212 * @param type side entity type
213 * @param data entity data
214 * @return error code
215 */
216 MoFEMErrorCode doWork(int side, EntityType type,
219 VectorDouble &vec = *dataPtr;
220 const size_t nb_gauss_pts = getGaussPts().size2();
221 if (type == zeroType || vec.size() != nb_gauss_pts) {
222 vec.resize(nb_gauss_pts, false);
223 vec.clear();
224 }
225
226 const size_t nb_dofs = data.getFieldData().size();
227
228 if (nb_dofs) {
229
230 if (dataVec.use_count()) {
231 dotVector.resize(nb_dofs, false);
232 const double *array;
233 CHKERR VecGetArrayRead(dataVec, &array);
234 const auto &local_indices = data.getLocalIndices();
235 for (int i = 0; i != local_indices.size(); ++i)
236 if (local_indices[i] != -1)
237 dotVector[i] = array[local_indices[i]];
238 else
239 dotVector[i] = 0;
240 CHKERR VecRestoreArrayRead(dataVec, &array);
241 data.getFieldData().swap(dotVector);
242 }
243
244 const size_t nb_base_functions = data.getN().size2();
245 auto base_function = data.getFTensor0N();
246 auto values_at_gauss_pts = getFTensor0FromVec(vec);
247 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
248 auto field_data = data.getFTensor0FieldData();
249 size_t bb = 0;
250 for (; bb != nb_dofs; ++bb) {
251 values_at_gauss_pts += field_data * base_function;
252 ++field_data;
253 ++base_function;
254 }
255 // It is possible to have more base functions than dofs
256 for (; bb < nb_base_functions; ++bb)
257 ++base_function;
258 ++values_at_gauss_pts;
259 }
260
261 if (dataVec.use_count()) {
262 data.getFieldData().swap(dotVector);
263 }
264 }
265
267 }
268};
269
270/**
271 * @brief Calculate scalar field values from PETSc vector at integration points
272 *
273 * This template structure extracts scalar field values from a PETSc vector
274 * and calculates them at integration points. It supports different data contexts
275 * through the template parameter and can handle various entity types.
276 *
277 * @tparam CTX PETSc data context type specifying how data is accessed
278 *
279 * @ingroup mofem_forces_and_sources_user_data_operators
280 */
281template <PetscData::DataContext CTX>
284
285 /**
286 * @brief Constructor for PETSc vector-based scalar field calculation
287 *
288 * @param field_name Name of the scalar field to extract values from
289 * @param data_ptr Shared pointer to VectorDouble for storing calculated values
290 * @param zero_at_type Entity type where values should be zeroed (default: MBVERTEX)
291 */
293 const std::string field_name, boost::shared_ptr<VectorDouble> data_ptr,
294 const EntityType zero_at_type = MBVERTEX)
297 dataPtr(data_ptr), zeroAtType(zero_at_type) {
298 if (!dataPtr)
299 THROW_MESSAGE("Pointer is not set");
300 }
301
302 MoFEMErrorCode doWork(int side, EntityType type,
305
306 const size_t nb_gauss_pts = getGaussPts().size2();
307
308 VectorDouble &vec = *dataPtr;
309 if (type == zeroAtType || vec.size() != nb_gauss_pts) {
310 vec.resize(nb_gauss_pts, false);
311 vec.clear();
312 }
313
314 auto &local_indices = data.getLocalIndices();
315 const size_t nb_dofs = local_indices.size();
316 if (nb_dofs) {
317
318 const double *array;
319
320 auto get_array = [&](const auto ctx, auto vec) {
322 #ifndef NDEBUG
323 if ((getFEMethod()->data_ctx & ctx).none()) {
324 MOFEM_LOG_CHANNEL("SELF");
325 MOFEM_LOG("SELF", Sev::error)
326 << "In this case field degrees of freedom are read from vector. "
327 "That usually happens when time solver is used, and acces to "
328 "first or second rates is needed. You probably not set ts_u, "
329 "ts_u_t, or ts_u_tt and associated data structure, i.e. "
330 "data_ctx to CTX_SET_X, CTX_SET_X_T, or CTX_SET_X_TT "
331 "respectively";
332 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Vector not set!");
333 }
334 #endif
335 CHKERR VecGetArrayRead(vec, &array);
337 };
338
339 auto restore_array = [&](auto vec) {
340 return VecRestoreArrayRead(vec, &array);
341 };
342
343 switch (CTX) {
345 CHKERR get_array(PetscData::CtxSetX, getFEMethod()->ts_u);
346 break;
348 CHKERR get_array(PetscData::CtxSetX_T, getFEMethod()->ts_u_t);
349 break;
351 CHKERR get_array(PetscData::CtxSetX_TT, getFEMethod()->ts_u_tt);
352 break;
353 default:
354 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
355 "That case is not implemented");
356 }
357
358 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
359 for (int i = 0; i != local_indices.size(); ++i)
360 if (local_indices[i] != -1)
361 dot_dofs_vector[i] = array[local_indices[i]];
362 else
363 dot_dofs_vector[i] = 0;
364
365 switch (CTX) {
367 CHKERR restore_array(getFEMethod()->ts_u);
368 break;
370 CHKERR restore_array(getFEMethod()->ts_u_t);
371 break;
373 CHKERR restore_array(getFEMethod()->ts_u_tt);
374 break;
375 default:
376 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
377 "That case is not implemented");
378 }
379
380 const size_t nb_base_functions = data.getN().size2();
381 auto base_function = data.getFTensor0N();
382 auto values_at_gauss_pts = getFTensor0FromVec(vec);
383
384 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
385 size_t bb = 0;
386 for (; bb != nb_dofs; ++bb) {
387 values_at_gauss_pts += dot_dofs_vector[bb] * base_function;
388 ++base_function;
389 }
390 // Number of dofs can be smaller than number of Tensor_Dim x base
391 // functions
392 for (; bb < nb_base_functions; ++bb)
393 ++base_function;
394 ++values_at_gauss_pts;
395 }
396 }
398 }
399
400private:
401 boost::shared_ptr<VectorDouble> dataPtr;
403};
404
409
410/**
411 * \deprecated Name inconsistent with other operators
412 *
413 */
415
416/**@}*/
417
418/** \name Vector field values at integration points */
419
420/**@{*/
421
422/**
423 * @brief Calculate field values for tensor field rank 1, i.e. vector field
424 *
425 * This template structure calculates vector field values at integration points
426 * and stores them in a matrix container. It supports various tensor dimensions,
427 * data types, and storage layouts through template parameters.
428 *
429 * @tparam Tensor_Dim Dimension of the vector field (e.g., 2 for 2D, 3 for 3D)
430 * @tparam T Data type for the field values (e.g., double, float)
431 * @tparam L Layout type for the ublas matrix storage
432 * @tparam A Allocator type for the ublas matrix storage
433 */
434template <int Tensor_Dim, class T, class L, class A>
437
438 /**
439 * @brief Constructor for vector field values calculation operator
440 *
441 * @param field_name Name of the vector field to calculate values for
442 * @param data_ptr Shared pointer to ublas matrix for storing calculated values
443 * @param zero_type Entity type for zero-level entities (default: MBVERTEX)
444 */
446 const std::string field_name,
447 boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
448 const EntityType zero_type = MBVERTEX)
451 dataPtr(data_ptr), zeroType(zero_type) {
452 if (!dataPtr)
453 THROW_MESSAGE("Pointer is not set");
454 }
455
456 /**
457 * \brief calculate values of vector field at integration points
458 * @param side side entity number
459 * @param type side entity type
460 * @param data entity data
461 * @return error code
462 */
463 MoFEMErrorCode doWork(int side, EntityType type,
465
466protected:
467 boost::shared_ptr<ublas::matrix<T, L, A>> dataPtr;
469};
470
471template <int Tensor_Dim, class T, class L, class A>
474 int side, EntityType type, EntitiesFieldData::EntData &data) {
476 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
477 "Not implemented for T = %s and dim = %d",
478 typeid(T).name(), // boost::core::demangle(typeid(T).name()),
479 Tensor_Dim);
481}
482
483/** \brief Calculate field values (template specialization) for tensor field
484 * rank 1, i.e. vector field
485 *
486 */
487template <int Tensor_Dim>
489 ublas::row_major, DoubleAllocator>
491
493 boost::shared_ptr<MatrixDouble> data_ptr,
494 const EntityType zero_type = MBVERTEX)
497 dataPtr(data_ptr), zeroType(zero_type) {
498 if (!dataPtr)
499 THROW_MESSAGE("Pointer is not set");
500 }
501
503 boost::shared_ptr<MatrixDouble> data_ptr,
504 SmartPetscObj<Vec> data_vec,
505 const EntityType zero_type = MBVERTEX)
508 dataPtr(data_ptr), dataVec(data_vec), zeroType(zero_type) {
509 if (!dataPtr)
510 THROW_MESSAGE("Pointer is not set");
511 }
512
513 MoFEMErrorCode doWork(int side, EntityType type,
515
516protected:
517 boost::shared_ptr<MatrixDouble> dataPtr;
521};
522
523/**
524 * \brief Member function specialization calculating values for tenor field rank
525 *
526 */
527template <int Tensor_Dim>
529 Tensor_Dim, double, ublas::row_major,
530 DoubleAllocator>::doWork(int side, EntityType type,
533
534 const size_t nb_gauss_pts = getGaussPts().size2();
535 auto &mat = *dataPtr;
536 if (type == zeroType || mat.size2() != nb_gauss_pts) {
537 mat.resize(Tensor_Dim, nb_gauss_pts, false);
538 mat.clear();
539 }
540
541 const size_t nb_dofs = data.getFieldData().size();
542 if (nb_dofs) {
543
544 if (dataVec.use_count()) {
545 dotVector.resize(nb_dofs, false);
546 const double *array;
547 CHKERR VecGetArrayRead(dataVec, &array);
548 const auto &local_indices = data.getLocalIndices();
549 for (int i = 0; i != local_indices.size(); ++i)
550 if (local_indices[i] != -1)
551 dotVector[i] = array[local_indices[i]];
552 else
553 dotVector[i] = 0;
554 CHKERR VecRestoreArrayRead(dataVec, &array);
555 data.getFieldData().swap(dotVector);
556 }
557
558 if (nb_gauss_pts) {
559 const size_t nb_base_functions = data.getN().size2();
560 auto base_function = data.getFTensor0N();
561 auto values_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
562 FTensor::Index<'I', Tensor_Dim> I;
563 const size_t size = nb_dofs / Tensor_Dim;
564 if (nb_dofs % Tensor_Dim) {
565 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
566 "Nb. of DOFs is inconsistent with Tensor_Dim");
567 }
568 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
569 auto field_data = data.getFTensor1FieldData<Tensor_Dim>();
570
571 #ifndef NDEBUG
572 if (field_data.l2() != field_data.l2()) {
573 MOFEM_LOG("SELF", Sev::error) << "field data: " << field_data;
574 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
575 "Wrong number in coefficients");
576 }
577 #endif
578
579 size_t bb = 0;
580 for (; bb != size; ++bb) {
581
582 #ifndef SINGULARITY
583 #ifndef NDEBUG
584 if (base_function != base_function) {
585 MOFEM_LOG("SELF", Sev::error) << "base function: " << base_function;
586 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
587 "Wrong number number in base functions");
588 }
589 #endif
590 #endif
591
592 values_at_gauss_pts(I) += field_data(I) * base_function;
593 ++field_data;
594 ++base_function;
595 }
596 // Number of dofs can be smaller than number of Tensor_Dim x base
597 // functions
598 for (; bb < nb_base_functions; ++bb)
599 ++base_function;
600 ++values_at_gauss_pts;
601 }
602 }
603
604 if (dataVec.use_count()) {
605 data.getFieldData().swap(dotVector);
606 }
607 }
609}
610
611/**
612 * @brief Specialization for double precision vector field values calculation
613 *
614 * This structure is a specialization of OpCalculateVectorFieldValues_General
615 * for double precision vector fields using row_major layout and DoubleAllocator.
616 * It provides a convenient interface for calculating vector field values at
617 * integration points.
618 *
619 * @tparam Tensor_Dim Dimension of the vector field (e.g., 2 for 2D, 3 for 3D)
620 *
621 * @ingroup mofem_forces_and_sources_user_data_operators
622 */
623template <int Tensor_Dim>
626 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
627
629 Tensor_Dim, double, ublas::row_major,
630 DoubleAllocator>::OpCalculateVectorFieldValues_General;
631};
632
633/**@}*/
634
635/** \name Vector field values at integration points */
636
637/**@{*/
638
639/**
640 * @brief Calculate divergence of vector field at integration points
641 *
642 * This template structure calculates the divergence of a vector field at
643 * integration points. It supports different coordinate systems and tensor
644 * dimensions, making it suitable for various physical applications.
645 *
646 * @tparam Tensor_Dim Dimension of the vector field (e.g., 2 for 2D, 3 for 3D)
647 * @tparam COORDINATE_SYSTEM Coordinate system type (default: CARTESIAN)
648 */
649template <int Tensor_Dim, CoordinateTypes COORDINATE_SYSTEM = CARTESIAN>
652
653 /**
654 * @brief Constructor for vector field divergence calculation operator
655 *
656 * @param field_name Name of the vector field to calculate divergence for
657 * @param data_ptr Shared pointer to VectorDouble for storing calculated divergence values
658 * @param zero_type Entity type for zero-level entities (default: MBVERTEX)
659 */
661 const std::string field_name, boost::shared_ptr<VectorDouble> data_ptr,
662 const EntityType zero_type = MBVERTEX)
665 dataPtr(data_ptr), zeroType(zero_type) {
666 if (!dataPtr)
667 THROW_MESSAGE("Pointer is not set");
668 }
669
670 MoFEMErrorCode doWork(int side, EntityType type,
673
674 // When we move to C++17 add if constexpr()
675 if constexpr (COORDINATE_SYSTEM == POLAR || COORDINATE_SYSTEM == SPHERICAL)
676 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
677 "%s coordiante not implemented",
678 CoordinateTypesNames[COORDINATE_SYSTEM]);
679
680 const size_t nb_gauss_pts = getGaussPts().size2();
681 auto &vec = *dataPtr;
682 if (type == zeroType) {
683 vec.resize(nb_gauss_pts, false);
684 vec.clear();
685 }
686
687 const size_t nb_dofs = data.getFieldData().size();
688 if (nb_dofs) {
689
690 if (nb_gauss_pts) {
691 const size_t nb_base_functions = data.getN().size2();
692 auto values_at_gauss_pts = getFTensor0FromVec(vec);
693 FTensor::Index<'I', Tensor_Dim> I;
694 const size_t size = nb_dofs / Tensor_Dim;
695 #ifndef NDEBUG
696 if (nb_dofs % Tensor_Dim) {
697 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
698 "Number of dofs should multiple of dimensions");
699 }
700 #endif
701
702 // When we move to C++17 add if constexpr()
703 if constexpr (COORDINATE_SYSTEM == CARTESIAN) {
704 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim>();
705 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
706 auto field_data = data.getFTensor1FieldData<Tensor_Dim>();
707 size_t bb = 0;
708 for (; bb != size; ++bb) {
709 values_at_gauss_pts += field_data(I) * diff_base_function(I);
710 ++field_data;
711 ++diff_base_function;
712 }
713 // Number of dofs can be smaller than number of Tensor_Dim x base
714 // functions
715 for (; bb < nb_base_functions; ++bb)
716 ++diff_base_function;
717 ++values_at_gauss_pts;
718 }
719 }
720
721 // When we move to C++17 add if constexpr()
722 if constexpr (COORDINATE_SYSTEM == CYLINDRICAL) {
723 auto t_coords = getFTensor1CoordsAtGaussPts();
724 auto values_at_gauss_pts = getFTensor0FromVec(vec);
725 auto base_function = data.getFTensor0N();
726 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim>();
727 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
728 auto field_data = data.getFTensor1FieldData<Tensor_Dim>();
729 size_t bb = 0;
730 for (; bb != size; ++bb) {
731 values_at_gauss_pts += field_data(I) * diff_base_function(I);
732 values_at_gauss_pts +=
733 base_function * (field_data(0) / t_coords(0));
734 ++field_data;
735 ++base_function;
736 ++diff_base_function;
737 }
738 // Number of dofs can be smaller than number of Tensor_Dim x base
739 // functions
740 for (; bb < nb_base_functions; ++bb) {
741 ++base_function;
742 ++diff_base_function;
743 }
744 ++values_at_gauss_pts;
745 ++t_coords;
746 }
747 }
748 }
749 }
751 }
752
753protected:
754 boost::shared_ptr<VectorDouble> dataPtr;
756};
757
758/** \brief Approximate field values for given petsc vector
759 *
760 * \note Look at PetscData to see what vectors could be extracted with that user
761 * data operator.
762 *
763 * \ingroup mofem_forces_and_sources_user_data_operators
764 */
765template <int Tensor_Dim, PetscData::DataContext CTX>
768
770 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
771 const EntityType zero_at_type = MBVERTEX, bool throw_error = true)
774 dataPtr(data_ptr), zeroAtType(zero_at_type), throwError(throw_error) {
775 if (!dataPtr)
776 THROW_MESSAGE("Pointer is not set");
777 }
778
779 MoFEMErrorCode doWork(int side, EntityType type,
782
783 auto &local_indices = data.getLocalIndices();
784 const size_t nb_dofs = local_indices.size();
785 const size_t nb_gauss_pts = getGaussPts().size2();
786
787 MatrixDouble &mat = *dataPtr;
788 if (type == zeroAtType) {
789 mat.resize(Tensor_Dim, nb_gauss_pts, false);
790 mat.clear();
791 }
792 if (!nb_dofs)
794
795 if (!throwError) {
796 if ((getFEMethod()->data_ctx & PetscData::Switches(CTX)).none()) {
798 }
799 }
800
801 const double *array;
802
803 auto get_array = [&](const auto ctx, auto vec) {
805 #ifndef NDEBUG
806 if ((getFEMethod()->data_ctx & ctx).none()) {
807 MOFEM_LOG_CHANNEL("SELF");
808 MOFEM_LOG("SELF", Sev::error)
809 << "In this case field degrees of freedom are read from vector. "
810 "That usually happens when time solver is used, and access to "
811 "first or second rates is needed. You probably not set ts_u, "
812 "ts_u_t, or ts_u_tt and associated data structure, i.e. "
813 "data_ctx to CTX_SET_X, CTX_SET_DX, CTX_SET_X_T, or "
814 "CTX_SET_X_TT respectively";
815 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Vector not set");
816 }
817 #endif
818 CHKERR VecGetArrayRead(vec, &array);
820 };
821
822 auto restore_array = [&](auto vec) {
823 return VecRestoreArrayRead(vec, &array);
824 };
825
826 switch (CTX) {
828 CHKERR get_array(PetscData::CtxSetX, getFEMethod()->ts_u);
829 break;
831 CHKERR get_array(PetscData::CtxSetX, getFEMethod()->dx);
832 break;
834 CHKERR get_array(PetscData::CtxSetX_T, getFEMethod()->ts_u_t);
835 break;
837 CHKERR get_array(PetscData::CtxSetX_TT, getFEMethod()->ts_u_tt);
838 break;
839 default:
840 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
841 "That case is not implemented");
842 }
843
844 dotVector.resize(local_indices.size());
845 for (int i = 0; i != local_indices.size(); ++i)
846 if (local_indices[i] != -1)
847 dotVector[i] = array[local_indices[i]];
848 else
849 dotVector[i] = 0;
850
851 switch (CTX) {
853 CHKERR restore_array(getFEMethod()->ts_u);
854 break;
856 CHKERR restore_array(getFEMethod()->dx);
857 break;
859 CHKERR restore_array(getFEMethod()->ts_u_t);
860 break;
862 CHKERR restore_array(getFEMethod()->ts_u_tt);
863 break;
864 default:
865 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
866 "That case is not implemented");
867 }
868
869 const size_t nb_base_functions = data.getN().size2();
870 auto base_function = data.getFTensor0N();
871 auto values_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
872
873 FTensor::Index<'I', Tensor_Dim> I;
874 const size_t size = nb_dofs / Tensor_Dim;
875 if (nb_dofs % Tensor_Dim) {
876 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
877 }
878 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
879 auto field_data = getFTensor1FromArray<Tensor_Dim, Tensor_Dim>(dotVector);
880 size_t bb = 0;
881 for (; bb != size; ++bb) {
882 values_at_gauss_pts(I) += field_data(I) * base_function;
883 ++field_data;
884 ++base_function;
885 }
886 // Number of dofs can be smaller than number of Tensor_Dim x base
887 // functions
888 for (; bb < nb_base_functions; ++bb)
889 ++base_function;
890 ++values_at_gauss_pts;
891 }
893 }
894
895protected:
896 boost::shared_ptr<MatrixDouble> dataPtr;
900};
901
902/** \brief Get rate of values at integration pts for tensor field
903 * rank 1, i.e. vector field
904 *
905 * \ingroup mofem_forces_and_sources_user_data_operators
906 */
907template <int Tensor_Dim>
911
912/** \brief Get second rate of values at integration pts for tensor
913 * field rank 1, i.e. vector field
914 *
915 * \ingroup mofem_forces_and_sources_user_data_operators
916 */
917template <int Tensor_Dim>
921
922/** \brief Get second time second update vector at integration pts for tensor
923 * field rank 1, i.e. vector field
924 *
925 * \ingroup mofem_forces_and_sources_user_data_operators
926 */
927template <int Tensor_Dim>
931
932/**@}*/
933
934/** \name Tensor field values at integration points */
935
936/**@{*/
937
938/** \brief Calculate field values for tenor field rank 2.
939 *
940 */
941template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
944
946 const std::string field_name,
947 boost::shared_ptr<ublas::matrix<T, L, A>> data_ptr,
948 const EntityType zero_type = MBVERTEX)
951 dataPtr(data_ptr), zeroType(zero_type) {
952 if (!dataPtr)
953 THROW_MESSAGE("Pointer is not set");
954 }
955
956 MoFEMErrorCode doWork(int side, EntityType type,
958
959protected:
960 boost::shared_ptr<ublas::matrix<T, L, A>> dataPtr;
962};
963
964template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
967 doWork(int side, EntityType type, EntitiesFieldData::EntData &data) {
969 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
970 "Not implemented for T = %s, dim0 = %d and dim1 = %d",
971 typeid(T).name(), // boost::core::demangle(typeid(T).name()),
972 Tensor_Dim0, Tensor_Dim1);
974}
975
976template <int Tensor_Dim0, int Tensor_Dim1>
977struct OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, double,
978 ublas::row_major, DoubleAllocator>
980
982 const std::string field_name,
983 boost::shared_ptr<
984 ublas::matrix<double, ublas::row_major, DoubleAllocator>>
985 data_ptr,
986 const EntityType zero_type = MBVERTEX)
989 dataPtr(data_ptr), zeroType(zero_type) {
990 if (!dataPtr)
991 THROW_MESSAGE("Pointer is not set");
992 }
993
995 const std::string field_name,
996 boost::shared_ptr<
997 ublas::matrix<double, ublas::row_major, DoubleAllocator>>
998 data_ptr,
999 SmartPetscObj<Vec> data_vec, const EntityType zero_type = MBVERTEX)
1002 dataPtr(data_ptr), zeroType(zero_type), dataVec(data_vec) {
1003 if (!dataPtr)
1004 THROW_MESSAGE("Pointer is not set");
1005 }
1006
1007 MoFEMErrorCode doWork(int side, EntityType type,
1009
1010protected:
1011 boost::shared_ptr<ublas::matrix<double, ublas::row_major, DoubleAllocator>>
1016};
1017
1018template <int Tensor_Dim0, int Tensor_Dim1>
1020 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1021 DoubleAllocator>::doWork(int side, EntityType type,
1024
1025 MatrixDouble &mat = *dataPtr;
1026 const size_t nb_gauss_pts = data.getN().size1();
1027 if (type == zeroType) {
1028 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
1029 mat.clear();
1030 }
1031
1032 const size_t nb_dofs = data.getFieldData().size();
1033
1034 if (dataVec.use_count()) {
1035 dotVector.resize(nb_dofs, false);
1036 const double *array;
1037 CHKERR VecGetArrayRead(dataVec, &array);
1038 const auto &local_indices = data.getLocalIndices();
1039 for (int i = 0; i != local_indices.size(); ++i)
1040 if (local_indices[i] != -1)
1041 dotVector[i] = array[local_indices[i]];
1042 else
1043 dotVector[i] = 0;
1044 CHKERR VecRestoreArrayRead(dataVec, &array);
1045 data.getFieldData().swap(dotVector);
1046 }
1047
1048 if (nb_dofs) {
1049 const size_t nb_base_functions = data.getN().size2();
1050 auto base_function = data.getFTensor0N();
1051 auto values_at_gauss_pts =
1052 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1053 FTensor::Index<'i', Tensor_Dim0> i;
1054 FTensor::Index<'j', Tensor_Dim1> j;
1055 const size_t size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
1056 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1057 auto field_data = data.getFTensor2FieldData<Tensor_Dim0, Tensor_Dim1>();
1058 size_t bb = 0;
1059 for (; bb != size; ++bb) {
1060 values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
1061 ++field_data;
1062 ++base_function;
1063 }
1064 for (; bb < nb_base_functions; ++bb)
1065 ++base_function;
1066 ++values_at_gauss_pts;
1067 }
1068
1069 if (dataVec.use_count()) {
1070 data.getFieldData().swap(dotVector);
1071 }
1072 }
1074}
1075
1076/** \brief Get values at integration pts for tensor field rank 2, i.e. matrix
1077 * field
1078 *
1079 * \ingroup mofem_forces_and_sources_user_data_operators
1080 */
1081template <int Tensor_Dim0, int Tensor_Dim1>
1084 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1085
1087 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1088 DoubleAllocator>::OpCalculateTensor2FieldValues_General;
1089};
1090
1091/** \brief Get time direvarive values at integration pts for tensor field rank
1092 * 2, i.e. matrix field
1093 *
1094 * \ingroup mofem_forces_and_sources_user_data_operators
1095 */
1096template <int Tensor_Dim0, int Tensor_Dim1>
1099
1101 boost::shared_ptr<MatrixDouble> data_ptr,
1102 const EntityType zero_at_type = MBVERTEX)
1105 dataPtr(data_ptr), zeroAtType(zero_at_type) {
1106 if (!dataPtr)
1107 THROW_MESSAGE("Pointer is not set");
1108 }
1109
1110 MoFEMErrorCode doWork(int side, EntityType type,
1113
1114 const size_t nb_gauss_pts = getGaussPts().size2();
1115 MatrixDouble &mat = *dataPtr;
1116 if (type == zeroAtType) {
1117 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
1118 mat.clear();
1119 }
1120 const auto &local_indices = data.getLocalIndices();
1121 const size_t nb_dofs = local_indices.size();
1122 if (nb_dofs) {
1123 dotVector.resize(nb_dofs, false);
1124 const double *array;
1125 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
1126 for (size_t i = 0; i != local_indices.size(); ++i)
1127 if (local_indices[i] != -1)
1128 dotVector[i] = array[local_indices[i]];
1129 else
1130 dotVector[i] = 0;
1131 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
1132
1133 const size_t nb_base_functions = data.getN().size2();
1134
1135 auto base_function = data.getFTensor0N();
1136 auto values_at_gauss_pts =
1137 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1138 FTensor::Index<'i', Tensor_Dim0> i;
1139 FTensor::Index<'j', Tensor_Dim1> j;
1140 const size_t size = nb_dofs / (Tensor_Dim0 * Tensor_Dim1);
1141 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1142 auto field_data =
1143 getFTensor2FromPtr<Tensor_Dim0, Tensor_Dim1>(&*dotVector.begin());
1144 size_t bb = 0;
1145 for (; bb != size; ++bb) {
1146 values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
1147 ++field_data;
1148 ++base_function;
1149 }
1150 for (; bb < nb_base_functions; ++bb)
1151 ++base_function;
1152 ++values_at_gauss_pts;
1153 }
1154 }
1156 }
1157
1158protected:
1159 boost::shared_ptr<MatrixDouble> dataPtr; ///< Data computed into this matrix
1160 EntityType zeroAtType; ///< Zero values at Gauss point at this type
1161 VectorDouble dotVector; ///< Keeps temporary values of time derivatives
1162};
1163
1164/**
1165 * @brief Calculate symmetric tensor field values at integration pts.
1166 *
1167 * @tparam Tensor_Dim
1168
1169 * \ingroup mofem_forces_and_sources_user_data_operators
1170 */
1171template <int Tensor_Dim>
1174
1176 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1177 const EntityType zero_type = MBEDGE, const int zero_side = 0)
1180 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1181 if (!dataPtr)
1182 THROW_MESSAGE("Pointer is not set");
1183 }
1184
1186 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1187 SmartPetscObj<Vec> data_vec, const EntityType zero_type = MBEDGE,
1188 const int zero_side = 0)
1191 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side),
1192 dataVec(data_vec) {
1193 if (!dataPtr)
1194 THROW_MESSAGE("Pointer is not set");
1195 }
1196
1197 MoFEMErrorCode doWork(int side, EntityType type,
1200 MatrixDouble &mat = *dataPtr;
1201 const int nb_gauss_pts = getGaussPts().size2();
1202 if (type == this->zeroType && side == zeroSide) {
1203 mat.resize((Tensor_Dim * (Tensor_Dim + 1)) / 2, nb_gauss_pts, false);
1204 mat.clear();
1205 }
1206 const int nb_dofs = data.getFieldData().size();
1207 if (!nb_dofs)
1209
1210 if (dataVec.use_count()) {
1211 dotVector.resize(nb_dofs, false);
1212 const double *array;
1213 CHKERR VecGetArrayRead(dataVec, &array);
1214 const auto &local_indices = data.getLocalIndices();
1215 for (int i = 0; i != local_indices.size(); ++i)
1216 if (local_indices[i] != -1)
1217 dotVector[i] = array[local_indices[i]];
1218 else
1219 dotVector[i] = 0;
1220 CHKERR VecRestoreArrayRead(dataVec, &array);
1221 data.getFieldData().swap(dotVector);
1222 }
1223
1224 const int nb_base_functions = data.getN().size2();
1225 auto base_function = data.getFTensor0N();
1226 auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
1227 FTensor::Index<'i', Tensor_Dim> i;
1228 FTensor::Index<'j', Tensor_Dim> j;
1229 const int size = nb_dofs / ((Tensor_Dim * (Tensor_Dim + 1)) / 2);
1230 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1231 auto field_data = data.getFTensor2SymmetricFieldData<Tensor_Dim>();
1232 int bb = 0;
1233 for (; bb != size; ++bb) {
1234 values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
1235 ++field_data;
1236 ++base_function;
1237 }
1238 for (; bb < nb_base_functions; ++bb)
1239 ++base_function;
1240 ++values_at_gauss_pts;
1241 }
1242
1243 if (dataVec.use_count()) {
1244 data.getFieldData().swap(dotVector);
1245 }
1246
1248 }
1249
1250protected:
1251 boost::shared_ptr<MatrixDouble> dataPtr;
1253 const int zeroSide;
1256};
1257
1258/**
1259 * @brief Calculate symmetric tensor field rates ant integratio pts.
1260 *
1261 * @tparam Tensor_Dim
1262 *
1263 * \ingroup mofem_forces_and_sources_user_data_operators
1264 */
1265template <int Tensor_Dim>
1268
1270 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1271 const EntityType zero_type = MBEDGE, const int zero_side = 0)
1274 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
1275 if (!dataPtr)
1276 THROW_MESSAGE("Pointer is not set");
1277 }
1278
1279 MoFEMErrorCode doWork(int side, EntityType type,
1282 const int nb_gauss_pts = getGaussPts().size2();
1283 MatrixDouble &mat = *dataPtr;
1284 constexpr auto symm_size = (Tensor_Dim * (Tensor_Dim + 1)) / 2;
1285 if (type == zeroType && side == zeroSide) {
1286 mat.resize(symm_size, nb_gauss_pts, false);
1287 mat.clear();
1288 }
1289 auto &local_indices = data.getLocalIndices();
1290 const int nb_dofs = local_indices.size();
1291 if (!nb_dofs)
1293
1294 #ifndef NDEBUG
1295 if ((getFEMethod()->data_ctx & PetscData::CtxSetX_T).none()) {
1296 MOFEM_LOG_CHANNEL("SELF");
1297 MOFEM_LOG("SELF", Sev::error)
1298 << "In this case field degrees of freedom are read from vector. "
1299 "That usually happens when time solver is used, and acces to "
1300 "first rates is needed. You probably not set "
1301 "ts_u_t and associated data structure data_ctx to CTX_SET_X_T "
1302 "respectively";
1303 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Vector not set!");
1304 }
1305 #endif
1306
1307 dotVector.resize(nb_dofs, false);
1308 const double *array;
1309 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
1310 for (int i = 0; i != local_indices.size(); ++i)
1311 if (local_indices[i] != -1)
1312 dotVector[i] = array[local_indices[i]];
1313 else
1314 dotVector[i] = 0;
1315 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
1316
1317 const int nb_base_functions = data.getN().size2();
1318
1319 auto base_function = data.getFTensor0N();
1320 auto values_at_gauss_pts = getFTensor2SymmetricFromMat<Tensor_Dim>(mat);
1321 FTensor::Index<'i', Tensor_Dim> i;
1322 FTensor::Index<'j', Tensor_Dim> j;
1323 const int size = nb_dofs / symm_size;
1324 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1325 auto field_data = getFTensorDotData<Tensor_Dim>();
1326 int bb = 0;
1327 for (; bb != size; ++bb) {
1328 values_at_gauss_pts(i, j) += field_data(i, j) * base_function;
1329 ++field_data;
1330 ++base_function;
1331 }
1332 for (; bb < nb_base_functions; ++bb)
1333 ++base_function;
1334 ++values_at_gauss_pts;
1335 }
1336
1338 }
1339
1340protected:
1341 boost::shared_ptr<MatrixDouble> dataPtr;
1343 const int zeroSide;
1345
1346 template <int Dim> inline auto getFTensorDotData() {
1347 static_assert(Dim || !Dim, "not implemented");
1348 }
1349};
1350
1351template <>
1352template <>
1353inline auto
1356 &dotVector[0], &dotVector[1], &dotVector[2], &dotVector[3], &dotVector[4],
1357 &dotVector[5]);
1358}
1359
1360template <>
1361template <>
1362inline auto
1365 &dotVector[0], &dotVector[1], &dotVector[2]);
1366}
1367
1368/**@}*/
1369
1370/** \name Gradients and Hessian of scalar fields at integration points */
1371
1372/**@{*/
1373
1374/**
1375 * \brief Evaluate field gradient values for scalar field, i.e. gradient is
1376 * tensor rank 1 (vector)
1377 *
1378 */
1379template <int Tensor_Dim, class T, class L, class A>
1381 : public OpCalculateVectorFieldValues_General<Tensor_Dim, T, L, A> {
1382
1384 Tensor_Dim, T, L, A>::OpCalculateVectorFieldValues_General;
1385};
1386
1387/** \brief Evaluate field gradient values for scalar field, i.e. gradient is
1388 * tensor rank 1 (vector), specialization
1389 *
1390 */
1391template <int Tensor_Dim>
1393 ublas::row_major, DoubleAllocator>
1395 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1396
1398 Tensor_Dim, double, ublas::row_major,
1399 DoubleAllocator>::OpCalculateVectorFieldValues_General;
1400
1401 /**
1402 * \brief calculate gradient values of scalar field at integration points
1403 * @param side side entity number
1404 * @param type side entity type
1405 * @param data entity data
1406 * @return error code
1407 */
1408 MoFEMErrorCode doWork(int side, EntityType type,
1410};
1411
1412/**
1413 * \brief Member function specialization calculating scalar field gradients for
1414 * tenor field rank 1
1415 *
1416 */
1417template <int Tensor_Dim>
1419 Tensor_Dim, double, ublas::row_major,
1420 DoubleAllocator>::doWork(int side, EntityType type,
1423
1424 const size_t nb_gauss_pts = this->getGaussPts().size2();
1425 auto &mat = *this->dataPtr;
1426 if (type == this->zeroType) {
1427 mat.resize(Tensor_Dim, nb_gauss_pts, false);
1428 mat.clear();
1429 }
1430
1431 const int nb_dofs = data.getFieldData().size();
1432 if (nb_dofs) {
1433
1434 if (nb_gauss_pts) {
1435 const int nb_base_functions = data.getN().size2();
1436 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim>();
1437 auto gradients_at_gauss_pts = getFTensor1FromMat<Tensor_Dim>(mat);
1438
1439 #ifndef NDEBUG
1440 if (nb_dofs > nb_base_functions)
1441 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1442 "Number of base functions inconsistent with number of DOFs "
1443 "(%d > %d)",
1444 nb_dofs, nb_base_functions);
1445
1446 if (data.getDiffN().size2() != nb_base_functions * Tensor_Dim)
1447 SETERRQ(
1448 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1449 "Number of base functions inconsistent with number of derivatives "
1450 "(%lu != %d)",
1451 data.getDiffN().size2(), nb_base_functions);
1452
1453 if (data.getDiffN().size1() != nb_gauss_pts)
1454 SETERRQ(
1455 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1456 "Number of base functions inconsistent with number of integration "
1457 "pts (%lu != %lu)",
1458 data.getDiffN().size2(), nb_gauss_pts);
1459
1460 #endif
1461
1462 FTensor::Index<'I', Tensor_Dim> I;
1463 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1464 auto field_data = data.getFTensor0FieldData();
1465 int bb = 0;
1466 for (; bb != nb_dofs; ++bb) {
1467 gradients_at_gauss_pts(I) += field_data * diff_base_function(I);
1468 ++field_data;
1469 ++diff_base_function;
1470 }
1471 // Number of dofs can be smaller than number of base functions
1472 for (; bb < nb_base_functions; ++bb)
1473 ++diff_base_function;
1474 ++gradients_at_gauss_pts;
1475 }
1476 }
1477 }
1478
1480}
1481
1482/** \brief Get field gradients at integration pts for scalar field rank 0, i.e.
1483 * vector field
1484 *
1485 * \ingroup mofem_forces_and_sources_user_data_operators
1486 */
1487template <int Tensor_Dim>
1490 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1492 Tensor_Dim, double, ublas::row_major,
1493 DoubleAllocator>::OpCalculateScalarFieldGradient_General;
1494};
1495
1496/** \brief Evaluate field gradient values for scalar field, i.e. gradient is
1497 * tensor rank 1 (vector), specialization
1498 *
1499 */
1500template <int Tensor_Dim>
1503 Tensor_Dim, double, ublas::row_major, DoubleAllocator> {
1504
1506 Tensor_Dim, double, ublas::row_major,
1507 DoubleAllocator>::OpCalculateVectorFieldValues_General;
1508
1509 /**
1510 * \brief calculate gradient values of scalar field at integration points
1511 * @param side side entity number
1512 * @param type side entity type
1513 * @param data entity data
1514 * @return error code
1515 */
1516 MoFEMErrorCode doWork(int side, EntityType type,
1518};
1519
1520template <int Tensor_Dim>
1522 int side, EntityType type, EntitiesFieldData::EntData &data) {
1524
1525 const size_t nb_gauss_pts = this->getGaussPts().size2();
1526 auto &mat = *this->dataPtr;
1527 if (type == this->zeroType) {
1528 mat.resize(Tensor_Dim * Tensor_Dim, nb_gauss_pts, false);
1529 mat.clear();
1530 }
1531
1532 const int nb_dofs = data.getFieldData().size();
1533 if (nb_dofs) {
1534
1535 if (nb_gauss_pts) {
1536 const int nb_base_functions = data.getN().size2();
1537
1538 auto &hessian_base = data.getN(BaseDerivatives::SecondDerivative);
1539 #ifndef NDEBUG
1540 if (hessian_base.size1() != nb_gauss_pts) {
1541 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1542 "Wrong number of integration pts (%ld != %d)",
1543 hessian_base.size1(), nb_gauss_pts);
1544 }
1545 if (hessian_base.size2() != nb_base_functions * Tensor_Dim * Tensor_Dim) {
1546 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1547 "Wrong number of base functions (%ld != %d)",
1548 hessian_base.size2() / (Tensor_Dim * Tensor_Dim),
1549 nb_base_functions);
1550 }
1551 if (hessian_base.size2() < nb_dofs * Tensor_Dim * Tensor_Dim) {
1552 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1553 "Wrong number of base functions (%ld < %d)",
1554 hessian_base.size2(), nb_dofs * Tensor_Dim * Tensor_Dim);
1555 }
1556 #endif
1557
1558 auto t_diff2_base_function = getFTensor2FromPtr<Tensor_Dim, Tensor_Dim>(
1559 &*hessian_base.data().begin());
1560
1561 auto hessian_at_gauss_pts =
1562 getFTensor2FromMat<Tensor_Dim, Tensor_Dim>(mat);
1563
1564 FTensor::Index<'I', Tensor_Dim> I;
1565 FTensor::Index<'J', Tensor_Dim> J;
1566 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1567 auto field_data = data.getFTensor0FieldData();
1568 int bb = 0;
1569 for (; bb != nb_dofs; ++bb) {
1570 hessian_at_gauss_pts(I, J) +=
1571 field_data * t_diff2_base_function(I, J);
1572 ++field_data;
1573 ++t_diff2_base_function;
1574 }
1575 // Number of dofs can be smaller than number of base functions
1576 for (; bb < nb_base_functions; ++bb) {
1577 ++t_diff2_base_function;
1578 }
1579
1580 ++hessian_at_gauss_pts;
1581 }
1582 }
1583 }
1584
1586}
1587
1588/**}*/
1589
1590/** \name Gradients and hessian of tensor fields at integration points */
1591
1592/**@{*/
1593
1594/**
1595 * \brief Evaluate field gradient values for vector field, i.e. gradient is
1596 * tensor rank 2
1597 *
1598 * \tparam Tensor_Dim0 Dimension of the vector field
1599 * \tparam Tensor_Dim1 Dimension of the gradient (usually spatial dimension)
1600 * \tparam S Storage order for the field data
1601 * \tparam T Data type
1602 * \tparam L Layout type
1603 * \tparam A Allocator type
1604 *
1605 */
1606template <int Tensor_Dim0, int Tensor_Dim1, int S, class T, class L, class A>
1608 : public OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T,
1609 L, A> {
1610
1612 Tensor_Dim0, Tensor_Dim1, T, L, A>::OpCalculateTensor2FieldValues_General;
1613};
1614
1615template <int Tensor_Dim0, int Tensor_Dim1, int S>
1617 Tensor_Dim0, Tensor_Dim1, S, double, ublas::row_major, DoubleAllocator>
1619 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1620
1622 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1623 DoubleAllocator>::OpCalculateTensor2FieldValues_General;
1624
1625 /**
1626 * \brief calculate values of vector field at integration points
1627 * @param side side entity number
1628 * @param type side entity type
1629 * @param data entity data
1630 * @return error code
1631 */
1632 MoFEMErrorCode doWork(int side, EntityType type,
1634};
1635
1636/**
1637 * \brief Member function specialization calculating vector field gradients for
1638 * tenor field rank 2
1639 *
1640 */
1641template <int Tensor_Dim0, int Tensor_Dim1, int S>
1643 Tensor_Dim0, Tensor_Dim1, S, double, ublas::row_major,
1644 DoubleAllocator>::doWork(int side, EntityType type,
1647 if (!this->dataPtr)
1648 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1649 "Data pointer not allocated");
1650
1651 const size_t nb_gauss_pts = this->getGaussPts().size2();
1652 auto &mat = *this->dataPtr;
1653 if (type == this->zeroType) {
1654 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
1655 mat.clear();
1656 }
1657
1658 if (nb_gauss_pts) {
1659 const size_t nb_dofs = data.getFieldData().size();
1660
1661 if (nb_dofs) {
1662
1663 if (this->dataVec.use_count()) {
1664 this->dotVector.resize(nb_dofs, false);
1665 const double *array;
1666 CHKERR VecGetArrayRead(this->dataVec, &array);
1667 const auto &local_indices = data.getLocalIndices();
1668 for (int i = 0; i != local_indices.size(); ++i)
1669 if (local_indices[i] != -1)
1670 this->dotVector[i] = array[local_indices[i]];
1671 else
1672 this->dotVector[i] = 0;
1673 CHKERR VecRestoreArrayRead(this->dataVec, &array);
1674 data.getFieldData().swap(this->dotVector);
1675 }
1676
1677 const int nb_base_functions = data.getN().size2();
1678 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
1679 auto gradients_at_gauss_pts =
1680 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1681 FTensor::Index<'I', Tensor_Dim0> I;
1682 FTensor::Index<'J', Tensor_Dim1> J;
1683 int size = nb_dofs / Tensor_Dim0;
1684 if (nb_dofs % Tensor_Dim0) {
1685 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1686 "Data inconsistency");
1687 }
1688 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1689 auto field_data = getFTensor1FromPtr<Tensor_Dim0, S>(
1690 data.getFieldData().data().data());
1691
1692 #ifndef NDEBUG
1693 if (field_data.l2() != field_data.l2()) {
1694 MOFEM_LOG("SELF", Sev::error) << "field data " << field_data;
1695 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1696 "Wrong number in coefficients");
1697 }
1698 #endif
1699
1700 int bb = 0;
1701 for (; bb < size; ++bb) {
1702 #ifndef SINGULARITY
1703 #ifndef NDEBUG
1704 if (diff_base_function.l2() != diff_base_function.l2()) {
1705 MOFEM_LOG("SELF", Sev::error)
1706 << "diff_base_function: " << diff_base_function;
1707 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1708 "Wrong number number in base functions");
1709 }
1710 #endif
1711 #endif
1712
1713 gradients_at_gauss_pts(I, J) += field_data(I) * diff_base_function(J);
1714 ++field_data;
1715 ++diff_base_function;
1716 }
1717 // Number of dofs can be smaller than number of Tensor_Dim0 x base
1718 // functions
1719 for (; bb != nb_base_functions; ++bb)
1720 ++diff_base_function;
1721 ++gradients_at_gauss_pts;
1722 }
1723
1724 if (this->dataVec.use_count()) {
1725 data.getFieldData().swap(this->dotVector);
1726 }
1727 }
1728 }
1730}
1731
1732/** \brief Get field gradients at integration pts for scalar field rank 0, i.e.
1733 * vector field
1734 *
1735 * \tparam Tensor_Dim0 Dimension of the vector field
1736 * \tparam Tensor_Dim1 Dimension of the gradient (usually spatial dimension)
1737 * \tparam S Stride in the storage of the vector field, default is Tensor_Dim0
1738 *
1739 * \ingroup mofem_forces_and_sources_user_data_operators
1740 */
1741template <int Tensor_Dim0, int Tensor_Dim1, int S = Tensor_Dim0>
1743 : public OpCalculateVectorFieldGradient_General<Tensor_Dim0, Tensor_Dim1, S,
1744 double, ublas::row_major,
1745 DoubleAllocator> {
1746
1748 Tensor_Dim0, Tensor_Dim1, S, double, ublas::row_major,
1749 DoubleAllocator>::OpCalculateVectorFieldGradient_General;
1750};
1751
1752/** \brief Get field gradients time derivative at integration pts for scalar
1753 * field rank 0, i.e. vector field
1754 *
1755 * \ingroup mofem_forces_and_sources_user_data_operators
1756 */
1757template <int Tensor_Dim0, int Tensor_Dim1>
1760
1762 boost::shared_ptr<MatrixDouble> data_ptr,
1763 const EntityType zero_at_type = MBVERTEX)
1766 dataPtr(data_ptr), zeroAtType(zero_at_type) {
1767 if (!dataPtr)
1768 THROW_MESSAGE("Pointer is not set");
1769 }
1770
1771 MoFEMErrorCode doWork(int side, EntityType type,
1774
1775 const auto &local_indices = data.getLocalIndices();
1776 const int nb_dofs = local_indices.size();
1777 const int nb_gauss_pts = this->getGaussPts().size2();
1778
1779 auto &mat = *this->dataPtr;
1780 if (type == this->zeroAtType) {
1781 mat.resize(Tensor_Dim0 * Tensor_Dim1, nb_gauss_pts, false);
1782 mat.clear();
1783 }
1784 if (!nb_dofs)
1786
1787 dotVector.resize(nb_dofs, false);
1788 const double *array;
1789 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
1790 for (int i = 0; i != local_indices.size(); ++i)
1791 if (local_indices[i] != -1)
1792 dotVector[i] = array[local_indices[i]];
1793 else
1794 dotVector[i] = 0;
1795 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
1796
1797 const int nb_base_functions = data.getN().size2();
1798 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
1799 auto gradients_at_gauss_pts =
1800 getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1801 FTensor::Index<'I', Tensor_Dim0> I;
1802 FTensor::Index<'J', Tensor_Dim1> J;
1803 int size = nb_dofs / Tensor_Dim0;
1804 if (nb_dofs % Tensor_Dim0) {
1805 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1806 }
1807
1808 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1809 auto field_data = getFTensor1FromPtr<Tensor_Dim0>(&*dotVector.begin());
1810 int bb = 0;
1811 for (; bb < size; ++bb) {
1812 gradients_at_gauss_pts(I, J) += field_data(I) * diff_base_function(J);
1813 ++field_data;
1814 ++diff_base_function;
1815 }
1816 // Number of dofs can be smaller than number of Tensor_Dim0 x base
1817 // functions
1818 for (; bb != nb_base_functions; ++bb)
1819 ++diff_base_function;
1820 ++gradients_at_gauss_pts;
1821 }
1823 }
1824
1825private:
1826 boost::shared_ptr<MatrixDouble> dataPtr; ///< Data computed into this matrix
1827 EntityType zeroAtType; ///< Zero values at Gauss point at this type
1828 VectorDouble dotVector; ///< Keeps temporary values of time derivatives
1829};
1830
1831/**
1832 * \brief Evaluate field gradient values for symmetric 2nd order tensor field,
1833 * i.e. gradient is tensor rank 3
1834 *
1835 */
1836template <int Tensor_Dim0, int Tensor_Dim1, class T, class L, class A>
1838 : public OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T,
1839 L, A> {
1840
1842 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1843 const EntityType zero_type = MBVERTEX)
1844 : OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, T, L,
1845 A>(field_name, data_ptr,
1846 zero_type) {}
1847};
1848
1849template <int Tensor_Dim0, int Tensor_Dim1>
1851 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator>
1853 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1854
1856 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1857 const EntityType zero_type = MBVERTEX)
1858 : OpCalculateTensor2FieldValues_General<Tensor_Dim0, Tensor_Dim1, double,
1859 ublas::row_major,
1861 field_name, data_ptr, zero_type) {}
1862
1863 /**
1864 * \brief calculate values of vector field at integration points
1865 * @param side side entity number
1866 * @param type side entity type
1867 * @param data entity data
1868 * @return error code
1869 */
1870 MoFEMErrorCode doWork(int side, EntityType type,
1872};
1873
1874/**
1875 * \brief Member function specialization calculating tensor field gradients for
1876 * symmetric tensor field rank 2
1877 *
1878 */
1879template <int Tensor_Dim0, int Tensor_Dim1>
1880MoFEMErrorCode OpCalculateTensor2SymmetricFieldGradient_General<
1881 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1882 DoubleAllocator>::doWork(int side, EntityType type,
1885 if (!this->dataPtr)
1886 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1887 "Data pointer not allocated");
1888
1889 const size_t nb_gauss_pts = this->getGaussPts().size2();
1890 constexpr size_t msize = (Tensor_Dim0 * (Tensor_Dim0 + 1)) / 2;
1891 auto &mat = *this->dataPtr;
1892 if (type == this->zeroType) {
1893 mat.resize(msize * Tensor_Dim1, nb_gauss_pts, false);
1894 mat.clear();
1895 }
1896
1897 if (nb_gauss_pts) {
1898 const size_t nb_dofs = data.getFieldData().size();
1899
1900 if (nb_dofs) {
1901
1902 const int nb_base_functions = data.getN().size2();
1903 auto diff_base_function = data.getFTensor1DiffN<Tensor_Dim1>();
1904 auto gradients_at_gauss_pts =
1905 getFTensor3DgFromMat<Tensor_Dim0, Tensor_Dim1>(mat);
1906 FTensor::Index<'I', Tensor_Dim0> I;
1907 FTensor::Index<'J', Tensor_Dim0> J;
1908 FTensor::Index<'K', Tensor_Dim1> K;
1909 int size = nb_dofs / msize;
1910 if (nb_dofs % msize) {
1911 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1912 "Data inconsistency");
1913 }
1914 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
1915 auto field_data = data.getFTensor2SymmetricFieldData<Tensor_Dim0>();
1916 int bb = 0;
1917 for (; bb < size; ++bb) {
1918 gradients_at_gauss_pts(I, J, K) +=
1919 field_data(I, J) * diff_base_function(K);
1920 ++field_data;
1921 ++diff_base_function;
1922 }
1923 // Number of dofs can be smaller than number of Tensor_Dim0 x base
1924 // functions
1925 for (; bb != nb_base_functions; ++bb)
1926 ++diff_base_function;
1927 ++gradients_at_gauss_pts;
1928 }
1929 }
1930 }
1932}
1933
1934/** \brief Get field gradients at integration pts for symmetric tensorial field
1935 * rank 2
1936 *
1937 * \ingroup mofem_forces_and_sources_user_data_operators
1938 */
1939template <int Tensor_Dim0, int Tensor_Dim1>
1942 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1943
1945 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
1946 const EntityType zero_type = MBVERTEX)
1948 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1949 DoubleAllocator>(field_name, data_ptr, zero_type) {}
1950};
1951
1952template <int Tensor_Dim0, int Tensor_Dim1>
1955 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator> {
1956
1958 Tensor_Dim0, Tensor_Dim1, double, ublas::row_major,
1959 DoubleAllocator>::OpCalculateTensor2FieldValues_General;
1960
1961 /**
1962 * \brief calculate values of vector field at integration points
1963 * @param side side entity number
1964 * @param type side entity type
1965 * @param data entity data
1966 * @return error code
1967 */
1968 MoFEMErrorCode doWork(int side, EntityType type,
1970};
1971
1972/**
1973 * \brief Member function specialization calculating vector field gradients for
1974 * tenor field rank 2
1975 *
1976 */
1977template <int Tensor_Dim0, int Tensor_Dim1>
1979 int side, EntityType type, EntitiesFieldData::EntData &data) {
1981 if (!this->dataPtr)
1982 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1983 "Data pointer not allocated");
1984
1985 const size_t nb_gauss_pts = this->getGaussPts().size2();
1986 auto &mat = *this->dataPtr;
1987 if (type == this->zeroType) {
1988 mat.resize(Tensor_Dim0 * Tensor_Dim1 * Tensor_Dim1, nb_gauss_pts, false);
1989 mat.clear();
1990 }
1991
1992 if (nb_gauss_pts) {
1993 const size_t nb_dofs = data.getFieldData().size();
1994
1995 if (nb_dofs) {
1996
1997 if (this->dataVec.use_count()) {
1998 this->dotVector.resize(nb_dofs, false);
1999 const double *array;
2000 CHKERR VecGetArrayRead(this->dataVec, &array);
2001 const auto &local_indices = data.getLocalIndices();
2002 for (int i = 0; i != local_indices.size(); ++i)
2003 if (local_indices[i] != -1)
2004 this->dotVector[i] = array[local_indices[i]];
2005 else
2006 this->dotVector[i] = 0;
2007 CHKERR VecRestoreArrayRead(this->dataVec, &array);
2008 data.getFieldData().swap(this->dotVector);
2009 }
2010
2011 const int nb_base_functions = data.getN().size2();
2012
2013 auto &hessian_base = data.getN(BaseDerivatives::SecondDerivative);
2014 #ifndef NDEBUG
2015 if (hessian_base.size1() != nb_gauss_pts) {
2016 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2017 "Wrong number of integration pts (%d != %d)",
2018 hessian_base.size1(), nb_gauss_pts);
2019 }
2020 if (hessian_base.size2() !=
2021 nb_base_functions * Tensor_Dim1 * Tensor_Dim1) {
2022 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2023 "Wrong number of base functions (%d != %d)",
2024 hessian_base.size2() / (Tensor_Dim1 * Tensor_Dim1),
2025 nb_base_functions);
2026 }
2027 if (hessian_base.size2() < nb_dofs * Tensor_Dim1 * Tensor_Dim1) {
2028 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2029 "Wrong number of base functions (%d < %d)",
2030 hessian_base.size2(), nb_dofs * Tensor_Dim1 * Tensor_Dim1);
2031 }
2032 #endif
2033
2034 auto t_diff2_base_function = getFTensor2FromPtr<Tensor_Dim1, Tensor_Dim1>(
2035 &*hessian_base.data().begin());
2036
2037 auto t_hessian_at_gauss_pts =
2038 getFTensor3FromMat<Tensor_Dim0, Tensor_Dim1, Tensor_Dim1>(mat);
2039
2040 FTensor::Index<'I', Tensor_Dim0> I;
2041 FTensor::Index<'J', Tensor_Dim1> J;
2042 FTensor::Index<'K', Tensor_Dim1> K;
2043
2044 int size = nb_dofs / Tensor_Dim0;
2045 #ifndef NDEBUG
2046 if (nb_dofs % Tensor_Dim0) {
2047 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2048 "Data inconsistency");
2049 }
2050 #endif
2051
2052 for (int gg = 0; gg < nb_gauss_pts; ++gg) {
2053 auto field_data = data.getFTensor1FieldData<Tensor_Dim0>();
2054 int bb = 0;
2055 for (; bb < size; ++bb) {
2056 t_hessian_at_gauss_pts(I, J, K) +=
2057 field_data(I) * t_diff2_base_function(J, K);
2058 ++field_data;
2059 ++t_diff2_base_function;
2060 }
2061 // Number of dofs can be smaller than number of Tensor_Dim0 x base
2062 // functions
2063 for (; bb != nb_base_functions; ++bb)
2064 ++t_diff2_base_function;
2065 ++t_hessian_at_gauss_pts;
2066 }
2067
2068 if (this->dataVec.use_count()) {
2069 data.getFieldData().swap(this->dotVector);
2070 }
2071 }
2072 }
2074}
2075
2076/**@}*/
2077
2078/** \name Transform tensors and vectors */
2079
2080/**@{*/
2081
2082/**
2083 * @brief Calculate \f$ \pmb\sigma_{ij} = \mathbf{D}_{ijkl} \pmb\varepsilon_{kl}
2084 * \f$
2085 *
2086 * @tparam DIM
2087 *
2088 * \ingroup mofem_forces_and_sources_user_data_operators
2089 */
2090template <int DIM_01, int DIM_23, int S = 0>
2093
2096
2097 /**
2098 * @deprecated Do not use this constructor
2099 */
2102 boost::shared_ptr<MatrixDouble> in_mat,
2103 boost::shared_ptr<MatrixDouble> out_mat,
2104 boost::shared_ptr<MatrixDouble> d_mat)
2105 : UserOp(field_name, OPROW), inMat(in_mat), outMat(out_mat), dMat(d_mat) {
2106 // Only is run for vertices
2107 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
2108 if (!inMat)
2109 THROW_MESSAGE("Pointer for in mat is null");
2110 if (!outMat)
2111 THROW_MESSAGE("Pointer for out mat is null");
2112 if (!dMat)
2113 THROW_MESSAGE("Pointer for tensor mat is null");
2114 }
2115
2116 OpTensorTimesSymmetricTensor(boost::shared_ptr<MatrixDouble> in_mat,
2117 boost::shared_ptr<MatrixDouble> out_mat,
2118 boost::shared_ptr<MatrixDouble> d_mat)
2119 : UserOp(NOSPACE, OPSPACE), inMat(in_mat), outMat(out_mat), dMat(d_mat) {
2120 // Only is run for vertices
2121 if (!inMat)
2122 THROW_MESSAGE("Pointer for in mat is null");
2123 if (!outMat)
2124 THROW_MESSAGE("Pointer for out mat is null");
2125 if (!dMat)
2126 THROW_MESSAGE("Pointer for tensor mat is null");
2127 }
2128
2129 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
2131 const size_t nb_gauss_pts = getGaussPts().size2();
2132 auto t_D = getFTensor4DdgFromMat<DIM_01, DIM_23, S>(*(dMat));
2133 auto t_in = getFTensor2SymmetricFromMat<DIM_01>(*(inMat));
2134 outMat->resize((DIM_23 * (DIM_23 + 1)) / 2, nb_gauss_pts, false);
2135 auto t_out = getFTensor2SymmetricFromMat<DIM_23>(*(outMat));
2136 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
2137 t_out(i, j) = t_D(i, j, k, l) * t_in(k, l);
2138 ++t_in;
2139 ++t_out;
2140 }
2142 }
2143
2144private:
2145 FTensor::Index<'i', DIM_01> i;
2146 FTensor::Index<'j', DIM_01> j;
2147 FTensor::Index<'k', DIM_23> k;
2148 FTensor::Index<'l', DIM_23> l;
2149
2150 boost::shared_ptr<MatrixDouble> inMat;
2151 boost::shared_ptr<MatrixDouble> outMat;
2152 boost::shared_ptr<MatrixDouble> dMat;
2153};
2154
2155/**
2156 * @brief Operator for symmetrizing tensor fields
2157 *
2158 * This template structure symmetrizes tensor fields by computing the symmetric
2159 * part of an input tensor and storing the result in an output matrix. It's
2160 * commonly used in mechanics where symmetric tensors (like stress and strain)
2161 * are required.
2162 *
2163 * @tparam DIM Dimension of the tensor field
2164 */
2165template <int DIM>
2167
2170
2171 /**
2172 * @deprecated Do not use this constructor
2173 */
2175 boost::shared_ptr<MatrixDouble> in_mat,
2176 boost::shared_ptr<MatrixDouble> out_mat)
2177 : UserOp(field_name, OPROW), inMat(in_mat), outMat(out_mat) {
2178 // Only is run for vertices
2179 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
2180 if (!inMat)
2181 THROW_MESSAGE("Pointer not set for in matrix");
2182 if (!outMat)
2183 THROW_MESSAGE("Pointer not set for in matrix");
2184 }
2185
2186 /**
2187 * @brief Constructor for tensor symmetrization operator
2188 *
2189 * @param in_mat Shared pointer to input matrix containing asymmetric tensor
2190 * @param out_mat Shared pointer to output matrix for storing symmetric tensor
2191 */
2192 OpSymmetrizeTensor(boost::shared_ptr<MatrixDouble> in_mat,
2193 boost::shared_ptr<MatrixDouble> out_mat)
2194 : UserOp(NOSPACE, OPSPACE), inMat(in_mat), outMat(out_mat) {
2195 // Only is run for vertices
2196 if (!inMat)
2197 THROW_MESSAGE("Pointer not set for in matrix");
2198 if (!outMat)
2199 THROW_MESSAGE("Pointer not set for in matrix");
2200 }
2201
2202 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
2204 const size_t nb_gauss_pts = getGaussPts().size2();
2205 auto t_in = getFTensor2FromMat<DIM, DIM>(*(inMat));
2206 outMat->resize((DIM * (DIM + 1)) / 2, nb_gauss_pts, false);
2207 auto t_out = getFTensor2SymmetricFromMat<DIM>(*(outMat));
2208 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
2209 t_out(i, j) = (t_in(i, j) || t_in(j, i)) / 2;
2210 ++t_in;
2211 ++t_out;
2212 }
2214 }
2215
2216private:
2219 boost::shared_ptr<MatrixDouble> inMat;
2220 boost::shared_ptr<MatrixDouble> outMat;
2221};
2222
2223/**
2224 * @brief Operator for scaling matrix values by a scalar factor
2225 *
2226 * This structure performs element-wise scaling of matrix data by multiplying
2227 * all elements by a scalar value. It's useful for applying material properties,
2228 * coordinate transformations, or other scaling operations in finite element
2229 * computations.
2230 */
2232
2235
2236 /**
2237 * @deprecated Do not use this constructor
2238 */
2239 DEPRECATED OpScaleMatrix(const std::string field_name, const double scale,
2240 boost::shared_ptr<MatrixDouble> in_mat,
2241 boost::shared_ptr<MatrixDouble> out_mat)
2242 : UserOp(field_name, OPROW), inMat(in_mat), outMat(out_mat) {
2243 scalePtr = boost::make_shared<double>(scale);
2244 // Only is run for vertices
2245 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
2246 if (!inMat)
2247 THROW_MESSAGE("Pointer not set for in matrix");
2248 if (!outMat)
2249 THROW_MESSAGE("Pointer not set for in matrix");
2250 }
2251
2252 /**
2253 * @brief Constructor for matrix scaling operator
2254 *
2255 * @param scale_ptr Shared pointer to scalar value for scaling matrix elements
2256 * @param in_mat Shared pointer to input matrix to be scaled
2257 * @param out_mat Shared pointer to output matrix for storing scaled results
2258 */
2259 OpScaleMatrix(boost::shared_ptr<double> scale_ptr,
2260 boost::shared_ptr<MatrixDouble> in_mat,
2261 boost::shared_ptr<MatrixDouble> out_mat)
2262 : UserOp(NOSPACE, OPSPACE), scalePtr(scale_ptr), inMat(in_mat),
2263 outMat(out_mat) {
2264 if (!scalePtr)
2265 THROW_MESSAGE("Pointer not set for scale");
2266 if (!inMat)
2267 THROW_MESSAGE("Pointer not set for in matrix");
2268 if (!outMat)
2269 THROW_MESSAGE("Pointer not set for in matrix");
2270 }
2271
2272 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
2274 outMat->resize(inMat->size1(), inMat->size2(), false);
2275 noalias(*outMat) = (*scalePtr) * (*inMat);
2277 }
2278
2279private:
2280 boost::shared_ptr<double> scalePtr;
2281 boost::shared_ptr<MatrixDouble> inMat;
2282 boost::shared_ptr<MatrixDouble> outMat;
2283};
2284
2285/**@}*/
2286
2287/** \name H-div/H-curls (Vectorial bases) values at integration points */
2288
2289/**@{*/
2290
2291/** \brief Get vector field for H-div approximation
2292 * \ingroup mofem_forces_and_sources_user_data_operators
2293 */
2294template <int Base_Dim, int Field_Dim, class T, class L, class A>
2296
2297/** \brief Get vector field for H-div approximation
2298 * \ingroup mofem_forces_and_sources_user_data_operators
2299 */
2300template <int Field_Dim>
2302 ublas::row_major, DoubleAllocator>
2304
2306 boost::shared_ptr<MatrixDouble> data_ptr,
2307 SmartPetscObj<Vec> data_vec,
2308 const EntityType zero_type = MBEDGE,
2309 const int zero_side = 0)
2312 dataPtr(data_ptr), dataVec(data_vec), zeroType(zero_type),
2313 zeroSide(zero_side) {
2314 if (!dataPtr)
2315 THROW_MESSAGE("Pointer is not set");
2316 }
2317
2319 boost::shared_ptr<MatrixDouble> data_ptr,
2320 const EntityType zero_type = MBEDGE,
2321 const int zero_side = 0)
2323 field_name, data_ptr, SmartPetscObj<Vec>(), zero_type, zero_side) {}
2324
2325 /**
2326 * \brief Calculate values of vector field at integration points
2327 * @param side side entity number
2328 * @param type side entity type
2329 * @param data entity data
2330 * @return error code
2331 */
2332 MoFEMErrorCode doWork(int side, EntityType type,
2334
2335private:
2336 boost::shared_ptr<MatrixDouble> dataPtr;
2339 const int zeroSide;
2341};
2342
2343template <int Field_Dim>
2345 3, Field_Dim, double, ublas::row_major,
2346 DoubleAllocator>::doWork(int side, EntityType type,
2349 const size_t nb_integration_points = this->getGaussPts().size2();
2350 if (type == zeroType && side == zeroSide) {
2351 dataPtr->resize(Field_Dim, nb_integration_points, false);
2352 dataPtr->clear();
2353 }
2354 const size_t nb_dofs = data.getFieldData().size();
2355 if (!nb_dofs)
2357
2358 if (dataVec.use_count()) {
2359 dotVector.resize(nb_dofs, false);
2360 const double *array;
2361 CHKERR VecGetArrayRead(dataVec, &array);
2362 const auto &local_indices = data.getLocalIndices();
2363 for (int i = 0; i != local_indices.size(); ++i)
2364 if (local_indices[i] != -1)
2365 dotVector[i] = array[local_indices[i]];
2366 else
2367 dotVector[i] = 0;
2368 CHKERR VecRestoreArrayRead(dataVec, &array);
2369 data.getFieldData().swap(dotVector);
2370 }
2371
2372 const size_t nb_base_functions = data.getN().size2() / 3;
2373 FTensor::Index<'i', Field_Dim> i;
2374 auto t_n_hdiv = data.getFTensor1N<3>();
2375 auto t_data = getFTensor1FromMat<Field_Dim>(*dataPtr);
2376 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2377 auto t_dof = data.getFTensor0FieldData();
2378 int bb = 0;
2379 for (; bb != nb_dofs; ++bb) {
2380 t_data(i) += t_n_hdiv(i) * t_dof;
2381 ++t_n_hdiv;
2382 ++t_dof;
2383 }
2384 for (; bb != nb_base_functions; ++bb)
2385 ++t_n_hdiv;
2386 ++t_data;
2387 }
2388
2389 if (dataVec.use_count()) {
2390 data.getFieldData().swap(dotVector);
2391 }
2393}
2394
2395/** \brief Get vector field for H-div approximation
2396 *
2397 * \ingroup mofem_forces_and_sources_user_data_operators
2398 */
2399template <int Base_Dim, int Field_Dim = Base_Dim>
2402 Base_Dim, Field_Dim, double, ublas::row_major, DoubleAllocator> {
2404 Base_Dim, Field_Dim, double, ublas::row_major,
2405 DoubleAllocator>::OpCalculateHVecVectorField_General;
2406};
2407
2408/** \brief Get vector field for H-div approximation
2409 * \ingroup mofem_forces_and_sources_user_data_operators
2410 */
2411template <int Base_Dim, int Field_Dim = Base_Dim>
2413
2414template <int Field_Dim>
2417
2419 boost::shared_ptr<MatrixDouble> data_ptr,
2420 const EntityType zero_type = MBEDGE,
2421 const int zero_side = 0)
2424 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2425 if (!dataPtr)
2426 THROW_MESSAGE("Pointer is not set");
2427 }
2428
2429 /**
2430 * \brief Calculate values of vector field at integration points
2431 * @param side side entity number
2432 * @param type side entity type
2433 * @param data entity data
2434 * @return error code
2435 */
2436 MoFEMErrorCode doWork(int side, EntityType type,
2438
2439private:
2440 boost::shared_ptr<MatrixDouble> dataPtr;
2442 const int zeroSide;
2443};
2444
2445template <int Field_Dim>
2447 int side, EntityType type, EntitiesFieldData::EntData &data) {
2449
2450 const size_t nb_integration_points = this->getGaussPts().size2();
2451 if (type == zeroType && side == zeroSide) {
2452 dataPtr->resize(Field_Dim, nb_integration_points, false);
2453 dataPtr->clear();
2454 }
2455
2456 auto &local_indices = data.getIndices();
2457 const size_t nb_dofs = local_indices.size();
2458 if (nb_dofs) {
2459
2460 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
2461 const double *array;
2462 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
2463 for (size_t i = 0; i != nb_dofs; ++i)
2464 if (local_indices[i] != -1)
2465 dot_dofs_vector[i] = array[local_indices[i]];
2466 else
2467 dot_dofs_vector[i] = 0;
2468 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
2469
2470 const size_t nb_base_functions = data.getN().size2() / 3;
2471 FTensor::Index<'i', Field_Dim> i;
2472 auto t_n_hdiv = data.getFTensor1N<3>();
2473 auto t_data = getFTensor1FromMat<Field_Dim>(*dataPtr);
2474 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2475 int bb = 0;
2476 for (; bb != nb_dofs; ++bb) {
2477 t_data(i) += t_n_hdiv(i) * dot_dofs_vector[bb];
2478 ++t_n_hdiv;
2479 }
2480 for (; bb != nb_base_functions; ++bb)
2481 ++t_n_hdiv;
2482 ++t_data;
2483 }
2484 }
2485
2487}
2488
2489/**
2490 * @brief Calculate divergence of vector field
2491 * @ingroup mofem_forces_and_sources_user_data_operators
2492 *
2493 * @tparam BASE_DIM
2494 * @tparam SPACE_DIM
2495 */
2496template <int BASE_DIM, int SPACE_DIM>
2499
2501 boost::shared_ptr<VectorDouble> data_ptr,
2502 const EntityType zero_type = MBEDGE,
2503 const int zero_side = 0)
2506 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2507 if (!dataPtr)
2508 THROW_MESSAGE("Pointer is not set");
2509 }
2510
2511 MoFEMErrorCode doWork(int side, EntityType type,
2514 const size_t nb_integration_points = getGaussPts().size2();
2515 if (type == zeroType && side == zeroSide) {
2516 dataPtr->resize(nb_integration_points, false);
2517 dataPtr->clear();
2518 }
2519 const size_t nb_dofs = data.getFieldData().size();
2520 if (!nb_dofs)
2522 const size_t nb_base_functions = data.getN().size2() / BASE_DIM;
2525 auto t_n_diff_hdiv = data.getFTensor2DiffN<BASE_DIM, SPACE_DIM>();
2526 auto t_data = getFTensor0FromVec(*dataPtr);
2527 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2528 auto t_dof = data.getFTensor0FieldData();
2529 int bb = 0;
2530 for (; bb != nb_dofs; ++bb) {
2531 t_data += t_dof * t_n_diff_hdiv(j, j);
2532 ++t_n_diff_hdiv;
2533 ++t_dof;
2534 }
2535 for (; bb != nb_base_functions; ++bb)
2536 ++t_n_diff_hdiv;
2537 ++t_data;
2538 }
2540 }
2541
2542private:
2543 boost::shared_ptr<VectorDouble> dataPtr;
2545 const int zeroSide;
2546};
2547
2548/**
2549 * @brief Calculate gradient of vector field
2550 * @ingroup mofem_forces_and_sources_user_data_operators
2551 *
2552 * @tparam BASE_DIM
2553 * @tparam SPACE_DIM
2554 */
2555template <int BASE_DIM, int SPACE_DIM>
2558
2560 boost::shared_ptr<MatrixDouble> data_ptr,
2561 const EntityType zero_type = MBEDGE,
2562 const int zero_side = 0)
2565 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2566 if (!dataPtr)
2567 THROW_MESSAGE("Pointer is not set");
2568 }
2569
2570 MoFEMErrorCode doWork(int side, EntityType type,
2573 const size_t nb_integration_points = getGaussPts().size2();
2574 if (type == zeroType && side == zeroSide) {
2575 dataPtr->resize(BASE_DIM * SPACE_DIM, nb_integration_points, false);
2576 dataPtr->clear();
2577 }
2578 const size_t nb_dofs = data.getFieldData().size();
2579 if (!nb_dofs)
2581 const size_t nb_base_functions = data.getN().size2() / BASE_DIM;
2584 auto t_base_diff = data.getFTensor2DiffN<BASE_DIM, SPACE_DIM>();
2585 auto t_data = getFTensor2FromMat<BASE_DIM, SPACE_DIM>(*dataPtr);
2586 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2587 auto t_dof = data.getFTensor0FieldData();
2588 int bb = 0;
2589 for (; bb != nb_dofs; ++bb) {
2590 t_data(i, j) += t_dof * t_base_diff(i, j);
2591 ++t_base_diff;
2592 ++t_dof;
2593 }
2594 for (; bb != nb_base_functions; ++bb)
2595 ++t_base_diff;
2596 ++t_data;
2597 }
2599 }
2600
2601private:
2602 boost::shared_ptr<MatrixDouble> dataPtr;
2604 const int zeroSide;
2605};
2606
2607/**
2608 * @brief Calculate gradient of tensor field
2609 * @ingroup mofem_forces_and_sources_user_data_operators
2610 *
2611 * @tparam BASE_DIM
2612 * @tparam FIELD_DIM
2613 * @tparam SPACE_DIM
2614 */
2615template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM>
2617
2618/**
2619 * @brief Specialisation for 3D tensor field gradient calculation
2620 * @ingroup mofem_forces_and_sources_user_data_operators
2621 *
2622 */
2623template <>
2626
2628 boost::shared_ptr<MatrixDouble> data_ptr,
2629 const EntityType zero_type = MBEDGE,
2630 const int zero_side = 0)
2633 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2634 if (!dataPtr)
2635 THROW_MESSAGE("Pointer is not set");
2636 }
2637
2638 MoFEMErrorCode doWork(int side, EntityType type,
2641
2642 const size_t nb_integration_points = getGaussPts().size2();
2643 if (type == zeroType && side == zeroSide) {
2644 dataPtr->resize(27, nb_integration_points,
2645 false);
2646 dataPtr->clear();
2647 }
2648
2649 #ifndef NDEBUG
2650 if (data.getFieldData().size() % 3) {
2651 SETERRQ(
2652 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2653 "Data inconsistency, nb_dofs %% COEFF_DIM != 0, that is %ld %% %d "
2654 "!= 0",
2655 data.getFieldData().size(), 3);
2656 }
2657 #endif
2658
2659 const auto nb_dofs = data.getFieldData().size() / 3;
2660 if (!nb_dofs)
2662
2663 const size_t nb_base_functions = data.getN().size2() / 3;
2664 FTensor::Index<'i', 3> i;
2665 FTensor::Index<'j', 3> j;
2666 FTensor::Index<'k', 3> k;
2667
2668 auto t_base_diff = data.getFTensor2DiffN<3, 3>();
2669 auto t_data = getFTensor3FromMat<3, 3, 3>(*dataPtr);
2670 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2671 auto t_dof = data.getFTensor1FieldData<3>();
2672 int bb = 0;
2673 for (; bb != nb_dofs; ++bb) {
2674 t_data(k, i, j) += t_base_diff(i, j) * t_dof(k);
2675 ++t_base_diff;
2676 ++t_dof;
2677 }
2678 for (; bb != nb_base_functions; ++bb)
2679 ++t_base_diff;
2680 ++t_data;
2681 }
2683 }
2684
2685private:
2686 boost::shared_ptr<MatrixDouble> dataPtr;
2688 const int zeroSide;
2689};
2690
2691template <>
2694
2696 boost::shared_ptr<MatrixDouble> data_ptr,
2697 const EntityType zero_type = MBEDGE,
2698 const int zero_side = 0)
2701 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2702 if (!dataPtr)
2703 THROW_MESSAGE("Pointer is not set");
2704 }
2705
2706 MoFEMErrorCode doWork(int side, EntityType type,
2709
2710 const size_t nb_integration_points = getGaussPts().size2();
2711 if (type == zeroType && side == zeroSide) {
2712 dataPtr->resize(27, nb_integration_points, false);
2713 dataPtr->clear();
2714 }
2715
2716 const auto nb_dofs = data.getFieldData().size();
2717 if (!nb_dofs)
2719
2720 const size_t nb_base_functions = data.getN().size2() / 9;
2721
2722 #ifndef NDEBUG
2723 if (data.getDiffN().size1() != nb_integration_points) {
2724 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2725 "Wrong number of integration pts (%ld != %ld)",
2726 data.getDiffN().size1(), nb_integration_points);
2727 }
2728 if (data.getDiffN().size2() != nb_base_functions * 27) {
2729 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2730 "Wrong number of base functions (%ld != %ld)",
2731 data.getDiffN().size2() / 27, nb_base_functions);
2732 }
2733 if (nb_base_functions < nb_dofs) {
2734 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2735 "Wrong number of base functions (%ld < %ld)", nb_base_functions,
2736 nb_dofs);
2737 }
2738 #endif
2739
2740 FTensor::Index<'i', 3> i;
2741 FTensor::Index<'j', 3> j;
2742 FTensor::Index<'k', 3> k;
2743
2744 auto t_base_diff = data.getFTensor3DiffN<3, 3, 3>();
2745 auto t_data = getFTensor3FromMat<3, 3, 3>(*dataPtr);
2746 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2747 auto t_dof = data.getFTensor0FieldData();
2748 int bb = 0;
2749 for (; bb != nb_dofs; ++bb) {
2750 t_data(k, i, j) += t_base_diff(k, i, j) * t_dof;
2751 ++t_base_diff;
2752 ++t_dof;
2753 }
2754 for (; bb != nb_base_functions; ++bb)
2755 ++t_base_diff;
2756 ++t_data;
2757 }
2758
2759
2761 }
2762
2763private:
2764 boost::shared_ptr<MatrixDouble> dataPtr;
2766 const int zeroSide;
2767};
2768
2769/**
2770 * @brief Calculate gradient of vector field
2771 * @ingroup mofem_forces_and_sources_user_data_operators
2772 *
2773 * @tparam BASE_DIM
2774 * @tparam SPACE_DIM
2775 */
2776template <int BASE_DIM, int SPACE_DIM>
2779
2781 boost::shared_ptr<MatrixDouble> data_ptr,
2782 const EntityType zero_type = MBEDGE,
2783 const int zero_side = 0)
2786 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2787 if (!dataPtr)
2788 THROW_MESSAGE("Pointer is not set");
2789 }
2790
2791 MoFEMErrorCode doWork(int side, EntityType type,
2794 const size_t nb_integration_points = getGaussPts().size2();
2795 if (type == zeroType && side == zeroSide) {
2796 dataPtr->resize(BASE_DIM * SPACE_DIM * SPACE_DIM, nb_integration_points,
2797 false);
2798 dataPtr->clear();
2799 }
2800 const size_t nb_dofs = data.getFieldData().size();
2801 if (!nb_dofs)
2803
2804 const int nb_base_functions = data.getN().size2() / BASE_DIM;
2805
2806 #ifndef NDEBUG
2807 auto &hessian_base = data.getN(BaseDerivatives::SecondDerivative);
2808 if (hessian_base.size1() != nb_integration_points) {
2809 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2810 "Wrong number of integration pts (%d != %d)",
2811 hessian_base.size1(), nb_integration_points);
2812 }
2813 if (hessian_base.size2() !=
2814 BASE_DIM * nb_base_functions * SPACE_DIM * SPACE_DIM) {
2815 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2816 "Wrong number of base functions (%d != %d)",
2817 hessian_base.size2() / (BASE_DIM * SPACE_DIM * SPACE_DIM),
2818 nb_base_functions);
2819 }
2820 if (hessian_base.size2() < BASE_DIM * nb_dofs * SPACE_DIM * SPACE_DIM) {
2821 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2822 "Wrong number of base functions (%d < %d)", hessian_base.size2(),
2823 BASE_DIM * nb_dofs * SPACE_DIM * SPACE_DIM);
2824 }
2825 #endif
2826
2830
2831 auto t_base_diff2 =
2833 auto t_data = getFTensor3FromMat<BASE_DIM, SPACE_DIM, SPACE_DIM>(*dataPtr);
2834 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2835 auto t_dof = data.getFTensor0FieldData();
2836 int bb = 0;
2837 for (; bb != nb_dofs; ++bb) {
2838 t_data(i, j, k) += t_dof * t_base_diff2(i, j, k);
2839
2840 ++t_base_diff2;
2841 ++t_dof;
2842 }
2843 for (; bb != nb_base_functions; ++bb)
2844 ++t_base_diff2;
2845 ++t_data;
2846 }
2848 }
2849
2850private:
2851 boost::shared_ptr<MatrixDouble> dataPtr;
2853 const int zeroSide;
2854};
2855
2856/**
2857 * @brief Calculate divergence of vector field dot
2858 * @ingroup mofem_forces_and_sources_user_data_operators
2859 *
2860 * @tparam Tensor_Dim dimension of space
2861 */
2862template <int Tensor_Dim1, int Tensor_Dim2>
2865
2867 boost::shared_ptr<VectorDouble> data_ptr,
2868 const EntityType zero_type = MBEDGE,
2869 const int zero_side = 0)
2872 dataPtr(data_ptr), zeroType(zero_type), zeroSide(zero_side) {
2873 if (!dataPtr)
2874 THROW_MESSAGE("Pointer is not set");
2875 }
2876
2877 MoFEMErrorCode doWork(int side, EntityType type,
2880 const size_t nb_integration_points = getGaussPts().size2();
2881 if (type == zeroType && side == zeroSide) {
2882 dataPtr->resize(nb_integration_points, false);
2883 dataPtr->clear();
2884 }
2885
2886 const auto &local_indices = data.getLocalIndices();
2887 const int nb_dofs = local_indices.size();
2888 if (nb_dofs) {
2889
2890 std::array<double, MAX_DOFS_ON_ENTITY> dot_dofs_vector;
2891 const double *array;
2892 CHKERR VecGetArrayRead(getFEMethod()->ts_u_t, &array);
2893 for (size_t i = 0; i != local_indices.size(); ++i)
2894 if (local_indices[i] != -1)
2895 dot_dofs_vector[i] = array[local_indices[i]];
2896 else
2897 dot_dofs_vector[i] = 0;
2898 CHKERR VecRestoreArrayRead(getFEMethod()->ts_u_t, &array);
2899
2900 const size_t nb_base_functions = data.getN().size2() / Tensor_Dim1;
2901 FTensor::Index<'i', Tensor_Dim1> i;
2902 auto t_n_diff_hdiv = data.getFTensor2DiffN<Tensor_Dim1, Tensor_Dim2>();
2903 auto t_data = getFTensor0FromVec(*dataPtr);
2904 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
2905 int bb = 0;
2906 for (; bb != nb_dofs; ++bb) {
2907 double div = 0;
2908 for (auto ii = 0; ii != Tensor_Dim2; ++ii)
2909 div += t_n_diff_hdiv(ii, ii);
2910 t_data += dot_dofs_vector[bb] * div;
2911 ++t_n_diff_hdiv;
2912 }
2913 for (; bb != nb_base_functions; ++bb)
2914 ++t_n_diff_hdiv;
2915 ++t_data;
2916 }
2917 }
2919 }
2920
2921private:
2922 boost::shared_ptr<VectorDouble> dataPtr;
2924 const int zeroSide;
2925};
2926
2927/**
2928 * @brief Calculate curl of vector field
2929 * @ingroup mofem_forces_and_sources_user_data_operators
2930 *
2931 * @tparam Base_Dim base function dimension
2932 * @tparam Space_Dim dimension of space
2933 * @tparam Hcurl field dimension
2934 */
2935template <int Base_Dim, int Space_Dim> struct OpCalculateHcurlVectorCurl;
2936
2937/**
2938 * @brief Calculate curl of vector field
2939 * @ingroup mofem_forces_and_sources_user_data_operators
2940 *
2941 * @tparam Base_Dim base function dimension
2942 * @tparam Space_Dim dimension of space
2943 * @tparam Hcurl field dimension
2944 */
2945template <>
2948 OpCalculateHcurlVectorCurl(const std::string field_name,
2949 boost::shared_ptr<MatrixDouble> data_ptr,
2950 const EntityType zero_type = MBEDGE,
2951 const int zero_side = 0);
2952 MoFEMErrorCode doWork(int side, EntityType type,
2954
2955private:
2956 boost::shared_ptr<MatrixDouble> dataPtr;
2958 const int zeroSide;
2959};
2960
2961/**
2962 * @brief Calculate curl of vector field
2963 * @ingroup mofem_forces_and_sources_user_data_operators
2964 *
2965 * @tparam Field_Dim dimension of field
2966 * @tparam Space_Dim dimension of space
2967 */
2968template <>
2971
2972 OpCalculateHcurlVectorCurl(const std::string field_name,
2973 boost::shared_ptr<MatrixDouble> data_ptr,
2974 const EntityType zero_type = MBVERTEX,
2975 const int zero_side = 0);
2976
2977 MoFEMErrorCode doWork(int side, EntityType type,
2979
2980private:
2981 boost::shared_ptr<MatrixDouble> dataPtr;
2983 const int zeroSide;
2984};
2985
2986/**
2987 * @brief Calculate curl of vector field
2988 * @ingroup mofem_forces_and_sources_user_data_operators
2989 *
2990 * @tparam Field_Dim dimension of field
2991 * @tparam Space_Dim dimension of space
2992 */
2993template <>
2996
2997 OpCalculateHcurlVectorCurl(const std::string field_name,
2998 boost::shared_ptr<MatrixDouble> data_ptr,
2999 const EntityType zero_type = MBVERTEX,
3000 const int zero_side = 0);
3001
3002 MoFEMErrorCode doWork(int side, EntityType type,
3004
3005private:
3006 boost::shared_ptr<MatrixDouble> dataPtr;
3008 const int zeroSide;
3009};
3010
3011/**
3012 * @brief Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl
3013 * \ingroup mofem_forces_and_sources_user_data_operators
3014 *
3015 * @tparam Tensor_Dim0 rank of the field
3016 * @tparam Tensor_Dim1 dimension of space
3017 */
3018template <int Tensor_Dim0, int Tensor_Dim1>
3021
3023 boost::shared_ptr<MatrixDouble> data_ptr,
3024 boost::shared_ptr<double> scale_ptr,
3026 const EntityType zero_type = MBEDGE,
3027 const int zero_side = 0)
3030 dataPtr(data_ptr), scalePtr(scale_ptr), dataVec(data_vec),
3031 zeroType(zero_type), zeroSide(zero_side) {
3032 if (!dataPtr)
3033 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Pointer is not set");
3034 }
3035
3037 boost::shared_ptr<MatrixDouble> data_ptr,
3038 const EntityType zero_type = MBEDGE,
3039 const int zero_side = 0)
3040 : OpCalculateHVecTensorField(field_name, data_ptr, nullptr,
3041 SmartPetscObj<Vec>(), zero_type, zero_side) {
3042 }
3043
3044 MoFEMErrorCode doWork(int side, EntityType type,
3047 const size_t nb_integration_points = getGaussPts().size2();
3048 if (type == zeroType && side == zeroSide) {
3049 dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points, false);
3050 dataPtr->clear();
3051 }
3052 const size_t nb_dofs = data.getFieldData().size();
3053 if (nb_dofs) {
3054
3055 if (dataVec.use_count()) {
3056 dotVector.resize(nb_dofs, false);
3057 const double *array;
3058 CHKERR VecGetArrayRead(dataVec, &array);
3059 const auto &local_indices = data.getLocalIndices();
3060 for (int i = 0; i != local_indices.size(); ++i)
3061 if (local_indices[i] != -1)
3062 dotVector[i] = array[local_indices[i]];
3063 else
3064 dotVector[i] = 0;
3065 CHKERR VecRestoreArrayRead(dataVec, &array);
3066 data.getFieldData().swap(dotVector);
3067 }
3068
3069 double scale = (scalePtr) ? *scalePtr : 1.0;
3070 const size_t nb_base_functions = data.getN().size2() / 3;
3071 FTensor::Index<'i', Tensor_Dim0> i;
3072 FTensor::Index<'j', Tensor_Dim1> j;
3073 auto t_n_hvec = data.getFTensor1N<3>();
3074 auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*dataPtr);
3075 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
3076 auto t_dof = data.getFTensor1FieldData<Tensor_Dim0>();
3077 size_t bb = 0;
3078 for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
3079 t_data(i, j) += (scale * t_dof(i)) * t_n_hvec(j);
3080 ++t_n_hvec;
3081 ++t_dof;
3082 }
3083 for (; bb < nb_base_functions; ++bb)
3084 ++t_n_hvec;
3085 ++t_data;
3086 }
3087
3088 if (dataVec.use_count()) {
3089 data.getFieldData().swap(dotVector);
3090 }
3091 }
3093 }
3094
3095private:
3096 boost::shared_ptr<MatrixDouble> dataPtr;
3097 boost::shared_ptr<double> scalePtr;
3100 const int zeroSide;
3101 VectorDouble dotVector; ///< Keeps temporary values of time derivatives
3102};
3103
3104/** \brief Get tensor field for H-div approximation
3105 * \ingroup mofem_forces_and_sources_user_data_operators
3106 *
3107 * \warning This operator is not tested
3108 */
3109template <int Tensor_Dim0, int Tensor_Dim1>
3112
3114
3116 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
3117 SmartPetscObj<Vec> data_vec, EntityType broken_type,
3118 boost::shared_ptr<Range> broken_range_ptr = nullptr,
3119 boost::shared_ptr<double> scale_ptr = nullptr,
3120 const EntityType zero_type = MBEDGE, const int zero_side = 0)
3122 dataPtr(data_ptr), dataVec(data_vec), brokenType(broken_type),
3123 brokenRangePtr(broken_range_ptr), zeroType(zero_type) {
3124 if (!dataPtr)
3125 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Pointer is not set");
3126 }
3127
3128 /**
3129 * \brief Calculate values of vector field at integration points
3130 * @param side side entity number
3131 * @param type side entity type
3132 * @param data entity data
3133 * @return error code
3134 */
3135 MoFEMErrorCode doWork(int side, EntityType type,
3137
3138private:
3139 boost::shared_ptr<MatrixDouble> dataPtr;
3141 EntityType brokenType;
3142 boost::shared_ptr<Range> brokenRangePtr;
3143 boost::shared_ptr<double> scalePtr;
3146};
3147
3148template <int Tensor_Dim0, int Tensor_Dim1>
3151 int side, EntityType type, EntitiesFieldData::EntData &data) {
3153 const size_t nb_integration_points = OP::getGaussPts().size2();
3154 if (type == zeroType) {
3155 dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points, false);
3156 dataPtr->clear();
3157 }
3158 const size_t nb_dofs = data.getFieldData().size();
3159 if (!nb_dofs)
3161
3162 if (dataVec.use_count()) {
3163 dotVector.resize(nb_dofs, false);
3164 const double *array;
3165 CHKERR VecGetArrayRead(dataVec, &array);
3166 const auto &local_indices = data.getLocalIndices();
3167 for (int i = 0; i != local_indices.size(); ++i)
3168 if (local_indices[i] != -1)
3169 dotVector[i] = array[local_indices[i]];
3170 else
3171 dotVector[i] = 0;
3172 CHKERR VecRestoreArrayRead(dataVec, &array);
3173 data.getFieldData().swap(dotVector);
3174 }
3175
3176 /**
3177 * @brief Get side face dofs
3178 *
3179 * Find which base functions on borken space have adjacent given entity type
3180 * and are in the range ptr if given.
3181 *
3182 */
3183 auto get_get_side_face_dofs = [&]() {
3184 auto fe_type = OP::getFEType();
3185
3186 BaseFunction::DofsSideMap &side_dof_map =
3187 data.getFieldEntities()[0]->getDofSideMap().at(fe_type);
3188 std::vector<int> side_face_dofs;
3189 side_face_dofs.reserve(data.getIndices().size() / Tensor_Dim0);
3190
3191 for (
3192
3193 auto it = side_dof_map.get<1>().begin();
3194 it != side_dof_map.get<1>().end(); ++it
3195
3196 ) {
3197 if ((Tensor_Dim0 * it->dof) >= data.getIndices().size()) {
3198 break;
3199 }
3200 if (it->type == brokenType) {
3201 if (brokenRangePtr) {
3202 auto ent = OP::getSideEntity(it->side, brokenType);
3203 if (brokenRangePtr->find(ent) != brokenRangePtr->end()) {
3204 side_face_dofs.push_back(it->dof);
3205 }
3206 } else {
3207 side_face_dofs.push_back(it->dof);
3208 }
3209 }
3210 }
3211
3212 return side_face_dofs;
3213 };
3214
3215 auto side_face_dofs = get_get_side_face_dofs();
3216
3217 FTensor::Index<'i', Tensor_Dim0> i;
3218 FTensor::Index<'j', Tensor_Dim1> j;
3219 auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*dataPtr);
3220 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
3221 for (auto b : side_face_dofs) {
3222 auto t_row_base = data.getFTensor1N<3>(gg, b);
3223 auto t_dof = getFTensor1FromPtr<Tensor_Dim0>(data.getFieldData().data() +
3224 b * Tensor_Dim0);
3225 t_data(i, j) += t_dof(i) * t_row_base(j);
3226 }
3227 ++t_data;
3228 }
3229 *dataPtr *= (scalePtr) ? *scalePtr : 1.0;
3230
3231 if (dataVec.use_count()) {
3232 data.getFieldData().swap(dotVector);
3233 }
3234
3236}
3237
3238/**
3239 * @brief Calculate tenor field using tensor base, i.e. Hdiv/Hcurl
3240 * \ingroup mofem_forces_and_sources_user_data_operators
3241 *
3242 * @tparam Tensor_Dim0 rank of the field
3243 * @tparam Tensor_Dim1 dimension of space
3244 */
3245template <int Tensor_Dim0, int Tensor_Dim1>
3248
3250 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
3251 boost::shared_ptr<double> scale_ptr,
3253 const EntityType zero_type = MBEDGE, const int zero_side = 0)
3256 dataPtr(data_ptr), scalePtr(scale_ptr), dataVec(data_vec),
3257 zeroType(zero_type), zeroSide(zero_side) {
3258 if (!dataPtr)
3259 THROW_MESSAGE("Pointer is not set");
3260 }
3261
3263 boost::shared_ptr<MatrixDouble> data_ptr,
3264 const EntityType zero_type = MBEDGE,
3265 const int zero_side = 0)
3266 : OpCalculateHTensorTensorField(field_name, data_ptr, nullptr,
3267 SmartPetscObj<Vec>(), zero_type,
3268 zero_side) {}
3269
3270 MoFEMErrorCode doWork(int side, EntityType type,
3273 const size_t nb_integration_points = getGaussPts().size2();
3274 if (type == zeroType && side == zeroSide) {
3275 dataPtr->resize(Tensor_Dim0 * Tensor_Dim1, nb_integration_points, false);
3276 dataPtr->clear();
3277 }
3278 const size_t nb_dofs = data.getFieldData().size();
3279 if (!nb_dofs)
3281
3282 if (dataVec.use_count()) {
3283 dotVector.resize(nb_dofs, false);
3284 const double *array;
3285 CHKERR VecGetArrayRead(dataVec, &array);
3286 const auto &local_indices = data.getLocalIndices();
3287 for (int i = 0; i != local_indices.size(); ++i)
3288 if (local_indices[i] != -1)
3289 dotVector[i] = array[local_indices[i]];
3290 else
3291 dotVector[i] = 0;
3292 CHKERR VecRestoreArrayRead(dataVec, &array);
3293 data.getFieldData().swap(dotVector);
3294 }
3295
3296 double scale = (scalePtr) ? *scalePtr : 1.0;
3297 const size_t nb_base_functions =
3298 data.getN().size2() / (Tensor_Dim0 * Tensor_Dim1);
3299 FTensor::Index<'i', Tensor_Dim0> i;
3300 FTensor::Index<'j', Tensor_Dim1> j;
3301 auto t_n_hten = data.getFTensor2N<Tensor_Dim0, Tensor_Dim1>();
3302 auto t_data = getFTensor2FromMat<Tensor_Dim0, Tensor_Dim1>(*dataPtr);
3303 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
3304 auto t_dof = data.getFTensor0FieldData();
3305 size_t bb = 0;
3306 for (; bb != nb_dofs; ++bb) {
3307 t_data(i, j) += (scale * t_dof) * t_n_hten(i, j);
3308 ++t_n_hten;
3309 ++t_dof;
3310 }
3311 for (; bb < nb_base_functions; ++bb)
3312 ++t_n_hten;
3313 ++t_data;
3314 }
3315
3316 if (dataVec.use_count()) {
3317 data.getFieldData().swap(dotVector);
3318 }
3319
3321 }
3322
3323private:
3324 boost::shared_ptr<MatrixDouble> dataPtr;
3325 boost::shared_ptr<double> scalePtr;
3328 const int zeroSide;
3329 VectorDouble dotVector; ///< Keeps temporary values of time derivatives
3330};
3331
3332/**
3333 * @brief Calculate divergence of tonsorial field using vectorial base
3334 * \ingroup mofem_forces_and_sources_user_data_operators
3335 *
3336 * @tparam Tensor_Dim0 rank of the field
3337 * @tparam Tensor_Dim1 dimension of space
3338 */
3339template <int Tensor_Dim0, int Tensor_Dim1,
3340 CoordinateTypes CoordSys = CARTESIAN>
3343
3345 boost::shared_ptr<MatrixDouble> data_ptr,
3346 SmartPetscObj<Vec> data_vec,
3347 const EntityType zero_type = MBEDGE,
3348 const int zero_side = 0)
3351 dataPtr(data_ptr), dataVec(data_vec), zeroType(zero_type),
3352 zeroSide(zero_side) {
3353 if (!dataPtr)
3354 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Pointer is not set");
3355 }
3356
3358 boost::shared_ptr<MatrixDouble> data_ptr,
3359 const EntityType zero_type = MBEDGE,
3360 const int zero_side = 0)
3362 field_name, data_ptr, SmartPetscObj<Vec>(), zero_type, zero_side) {}
3363
3364 MoFEMErrorCode doWork(int side, EntityType type,
3367 const size_t nb_integration_points = getGaussPts().size2();
3368 if (type == zeroType && side == zeroSide) {
3369 dataPtr->resize(Tensor_Dim0, nb_integration_points, false);
3370 dataPtr->clear();
3371 }
3372 const size_t nb_dofs = data.getFieldData().size();
3373 if (nb_dofs) {
3374
3375 if (dataVec.use_count()) {
3376 dotVector.resize(nb_dofs, false);
3377 const double *array;
3378 CHKERR VecGetArrayRead(dataVec, &array);
3379 const auto &local_indices = data.getLocalIndices();
3380 for (int i = 0; i != local_indices.size(); ++i)
3381 if (local_indices[i] != -1)
3382 dotVector[i] = array[local_indices[i]];
3383 else
3384 dotVector[i] = 0;
3385 CHKERR VecRestoreArrayRead(dataVec, &array);
3386 data.getFieldData().swap(dotVector);
3387 }
3388
3389 const size_t nb_base_functions = data.getN().size2() / 3;
3390 FTensor::Index<'i', Tensor_Dim0> i;
3391 FTensor::Index<'j', Tensor_Dim1> j;
3392 auto t_n_diff_hvec = data.getFTensor2DiffN<3, Tensor_Dim1>();
3393 auto t_data = getFTensor1FromMat<Tensor_Dim0>(*dataPtr);
3394 auto t_base = data.getFTensor1N<3>();
3395 auto t_coords = getFTensor1CoordsAtGaussPts();
3396 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
3397 auto t_dof = data.getFTensor1FieldData<Tensor_Dim0>();
3398 size_t bb = 0;
3399 for (; bb != nb_dofs / Tensor_Dim0; ++bb) {
3400 double div = t_n_diff_hvec(j, j);
3401 t_data(i) += t_dof(i) * div;
3402 if constexpr (CoordSys == CYLINDRICAL) {
3403 t_data(i) += t_base(0) * (t_dof(i) / t_coords(0));
3404 }
3405 ++t_n_diff_hvec;
3406 ++t_dof;
3407 ++t_base;
3408 }
3409 for (; bb < nb_base_functions; ++bb) {
3410 ++t_base;
3411 ++t_n_diff_hvec;
3412 }
3413 ++t_data;
3414 ++t_coords;
3415 }
3416
3417 if (dataVec.use_count()) {
3418 data.getFieldData().swap(dotVector);
3419 }
3420 }
3422 }
3423
3424private:
3425 boost::shared_ptr<MatrixDouble> dataPtr;
3428 const int zeroSide;
3429
3430 VectorDouble dotVector; ///< Keeps temporary values of time derivatives
3431};
3432
3433/**
3434 * @brief Calculate divergence of tonsorial field using vectorial base
3435 * \ingroup mofem_forces_and_sources_user_data_operators
3436 *
3437 * \warning This operator is not tested
3438 *
3439 * @tparam Tensor_Dim0 rank of the field
3440 * @tparam Tensor_Dim1 dimension of space
3441 */
3442template <int Tensor_Dim0, int Tensor_Dim1,
3443 CoordinateTypes CoordSys = CARTESIAN>
3446
3448
3450 const std::string field_name, boost::shared_ptr<MatrixDouble> data_ptr,
3451 SmartPetscObj<Vec> data_vec, EntityType broken_type,
3452 boost::shared_ptr<Range> broken_range_ptr = nullptr,
3453 boost::shared_ptr<double> scale_ptr = nullptr,
3454 const EntityType zero_type = MBEDGE)
3456 dataPtr(data_ptr), dataVec(data_vec), brokenType(broken_type),
3457 brokenRangePtr(broken_range_ptr), scalePtr(scale_ptr),
3458 zeroType(zero_type) {
3459 if (!dataPtr)
3460 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Pointer is not set");
3461 }
3462
3463 MoFEMErrorCode doWork(int side, EntityType type,
3466 const size_t nb_integration_points = getGaussPts().size2();
3467 if (type == zeroType && side == 0) {
3468 dataPtr->resize(Tensor_Dim0, nb_integration_points, false);
3469 dataPtr->clear();
3470 }
3471
3472 const size_t nb_dofs = data.getFieldData().size();
3473 if (nb_dofs) {
3474
3475 if (dataVec.use_count()) {
3476 dotVector.resize(nb_dofs, false);
3477 const double *array;
3478 CHKERR VecGetArrayRead(dataVec, &array);
3479 const auto &local_indices = data.getLocalIndices();
3480 for (int i = 0; i != local_indices.size(); ++i)
3481 if (local_indices[i] != -1)
3482 dotVector[i] = array[local_indices[i]];
3483 else
3484 dotVector[i] = 0;
3485 CHKERR VecRestoreArrayRead(dataVec, &array);
3486 data.getFieldData().swap(dotVector);
3487 }
3488
3489 /**
3490 * @brief Get side face dofs
3491 *
3492 * Find which base functions on borken space have adjacent given entity
3493 * type and are in the range ptr if given.
3494 *
3495 */
3496 auto get_get_side_face_dofs = [&]() {
3497 auto fe_type = OP::getFEType();
3498
3499 BaseFunction::DofsSideMap &side_dof_map =
3500 data.getFieldEntities()[0]->getDofSideMap().at(fe_type);
3501 std::vector<int> side_face_dofs;
3502 side_face_dofs.reserve(data.getIndices().size() / Tensor_Dim0);
3503
3504 for (
3505
3506 auto it = side_dof_map.get<1>().begin();
3507 it != side_dof_map.get<1>().end(); ++it
3508
3509 ) {
3510 if ((Tensor_Dim0 * it->dof) >= data.getIndices().size()) {
3511 break;
3512 }
3513 if (it->type == brokenType) {
3514 if (brokenRangePtr) {
3515 auto ent = OP::getSideEntity(it->side, brokenType);
3516 if (brokenRangePtr->find(ent) != brokenRangePtr->end()) {
3517 side_face_dofs.push_back(it->dof);
3518 }
3519 } else {
3520 side_face_dofs.push_back(it->dof);
3521 }
3522 }
3523 }
3524
3525 return side_face_dofs;
3526 };
3527
3528 auto side_face_dofs = get_get_side_face_dofs();
3529
3530 FTensor::Index<'i', Tensor_Dim0> i;
3531 FTensor::Index<'j', Tensor_Dim1> j;
3532 auto t_data = getFTensor1FromMat<Tensor_Dim0>(*dataPtr);
3533 auto t_coords = getFTensor1CoordsAtGaussPts();
3534 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
3535 for (auto b : side_face_dofs) {
3536 auto t_dof = getFTensor1FromPtr<Tensor_Dim0>(
3537 data.getFieldData().data() + b * Tensor_Dim0);
3538 auto t_base = data.getFTensor1N<3>(gg, b);
3539 auto t_diff_base = data.getFTensor2DiffN<3, Tensor_Dim1>(gg, b);
3540 double div = t_diff_base(j, j);
3541 t_data(i) += t_dof(i) * div;
3542 if constexpr (CoordSys == CYLINDRICAL) {
3543 t_data(i) += t_base(0) * (t_dof(i) / t_coords(0));
3544 }
3545 }
3546 ++t_data;
3547 ++t_coords;
3548 }
3549 }
3550
3551 if (dataVec.use_count()) {
3552 data.getFieldData().swap(dotVector);
3553 }
3554
3556 }
3557
3558private:
3559 boost::shared_ptr<MatrixDouble> dataPtr;
3561 EntityType brokenType;
3562 boost::shared_ptr<Range> brokenRangePtr;
3563 boost::shared_ptr<double> scalePtr;
3566};
3567
3568/**
3569 * @brief Calculate trace of vector (Hdiv/Hcurl) space
3570 *
3571 * @tparam Tensor_Dim
3572 * @tparam OpBase
3573 */
3574template <int Tensor_Dim, typename OpBase>
3576
3578 boost::shared_ptr<MatrixDouble> data_ptr,
3579 boost::shared_ptr<double> scale_ptr,
3580 const EntityType zero_type = MBEDGE,
3581 const int zero_side = 0)
3582 : OpBase(field_name, OpBase::OPROW), dataPtr(data_ptr),
3583 scalePtr(scale_ptr), zeroType(zero_type), zeroSide(zero_side) {
3584 if (!dataPtr)
3585 THROW_MESSAGE("Pointer is not set");
3586 }
3587
3589 boost::shared_ptr<MatrixDouble> data_ptr,
3590 const EntityType zero_type = MBEDGE,
3591 const int zero_side = 0)
3592 : OpCalculateHVecTensorTrace(field_name, data_ptr, nullptr, zero_type,
3593 zero_side) {}
3594
3595 MoFEMErrorCode doWork(int side, EntityType type,
3598 const size_t nb_integration_points = OpBase::getGaussPts().size2();
3599 if (type == zeroType && side == 0) {
3600 dataPtr->resize(Tensor_Dim, nb_integration_points, false);
3601 dataPtr->clear();
3602 }
3603 const size_t nb_dofs = data.getFieldData().size();
3604 if (nb_dofs) {
3605 double scale_val = (scalePtr) ? *scalePtr : 1.0;
3606 auto t_normal = OpBase::getFTensor1NormalsAtGaussPts();
3607 const size_t nb_base_functions = data.getN().size2() / 3;
3608 auto t_base = data.getFTensor1N<3>();
3609 auto t_data = getFTensor1FromMat<Tensor_Dim>(*dataPtr);
3610 for (size_t gg = 0; gg != nb_integration_points; ++gg) {
3611 FTensor::Tensor1<double, Tensor_Dim> t_normalized_normal;
3612 t_normalized_normal(j) = t_normal(j);
3613 t_normalized_normal.normalize();
3614 auto t_dof = data.getFTensor1FieldData<Tensor_Dim>();
3615 size_t bb = 0;
3616 for (; bb != nb_dofs / Tensor_Dim; ++bb) {
3617 t_data(i) +=
3618 (scale_val * t_dof(i)) * (t_base(j) * t_normalized_normal(j));
3619 ++t_base;
3620 ++t_dof;
3621 }
3622 for (; bb < nb_base_functions; ++bb) {
3623 ++t_base;
3624 }
3625 ++t_data;
3626 ++t_normal;
3627 }
3628 }
3630 }
3631
3632private:
3633 boost::shared_ptr<MatrixDouble> dataPtr;
3634 boost::shared_ptr<double> scalePtr;
3636 const int zeroSide;
3637 FTensor::Index<'i', Tensor_Dim> i;
3638 FTensor::Index<'j', Tensor_Dim> j;
3639};
3640
3641/**@}*/
3642
3643/** \name Other operators */
3644
3645/**@{*/
3646
3647/**@}*/
3648
3649/** \name Operators for faces */
3650
3651/**@{*/
3652
3653/** \brief Transform local reference derivatives of shape functions to global
3654derivatives
3655
3656\ingroup mofem_forces_and_sources_tri_element
3657
3658*/
3659template <int DIM, int DERIVATIVE = 1> struct OpSetInvJacSpaceForFaceImpl;
3660
3663
3665 boost::shared_ptr<MatrixDouble> inv_jac_ptr);
3666
3667protected:
3668
3669 /**
3670 * @brief Apply transformation to the input matrix
3671 *
3672 * @tparam D1 dimension of the derivative of base functions in input
3673 * @tparam D2 dimension of the derivative of base functions in output
3674 * @tparam J1 nb of rows in jacobian (= dimension of space)
3675 * @tparam J2 nb of columns in jacobian (= dimension of reference element)
3676 * @param diff_n
3677 * @return MoFEMErrorCode
3678 */
3679 template <int D1, int D2, int J1, int J2>
3682
3683 static_assert(D2 == J2, "Dimension of jacobian and dimension of <out> "
3684 "directive does not match");
3685
3686 size_t nb_functions = diff_n.size2() / D1;
3687 if (nb_functions) {
3688 size_t nb_gauss_pts = diff_n.size1();
3689
3690 #ifndef NDEBUG
3691 if (nb_gauss_pts != getGaussPts().size2())
3692 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3693 "Wrong number of Gauss Pts");
3694 if (diff_n.size2() % D1)
3695 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3696 "Number of directives of base functions and D1 dimension does "
3697 "not match");
3698 #endif
3699
3700 diffNinvJac.resize(diff_n.size1(), D2 * nb_functions, false);
3701
3702 FTensor::Index<'i', D2> i;
3703 FTensor::Index<'k', D1> k;
3704 auto t_diff_n = getFTensor1FromPtr<D2>(&*diffNinvJac.data().begin());
3705 auto t_diff_n_ref = getFTensor1FromPtr<D1>(&*diff_n.data().begin());
3706 auto t_inv_jac = getFTensor2FromMat<J1, J2>(*invJacPtr);
3707 for (size_t gg = 0; gg != nb_gauss_pts; ++gg, ++t_inv_jac) {
3708 for (size_t dd = 0; dd != nb_functions; ++dd) {
3709 t_diff_n(i) = t_inv_jac(k, i) * t_diff_n_ref(k);
3710 ++t_diff_n;
3711 ++t_diff_n_ref;
3712 }
3713 }
3714
3715 diff_n.swap(diffNinvJac);
3716 }
3718 }
3719
3720 boost::shared_ptr<MatrixDouble> invJacPtr;
3722};
3723
3724template <>
3727
3729
3730 MoFEMErrorCode doWork(int side, EntityType type,
3732};
3733
3734template <>
3736 : public OpSetInvJacSpaceForFaceImpl<2, 1> {
3737
3738 using OpSetInvJacSpaceForFaceImpl<2, 1>::OpSetInvJacSpaceForFaceImpl;
3739
3740 MoFEMErrorCode doWork(int side, EntityType type,
3742};
3743
3744template <>
3746 : public OpSetInvJacSpaceForFaceImpl<2, 1> {
3747
3748 using OpSetInvJacSpaceForFaceImpl<2, 1>::OpSetInvJacSpaceForFaceImpl;
3749
3750 MoFEMErrorCode doWork(int side, EntityType type,
3752};
3753
3754template <int DERIVARIVE = 1>
3756 : public OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE> {
3757 OpSetInvJacH1ForFace(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3758 : OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE>(H1, inv_jac_ptr) {}
3759};
3760
3761template <int DERIVARIVE = 1>
3763 : public OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE> {
3764 OpSetInvJacL2ForFace(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3765 : OpSetInvJacSpaceForFaceImpl<2, DERIVARIVE>(L2, inv_jac_ptr) {}
3766};
3767
3768template <int DERIVARIVE = 1>
3770 : public OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE> {
3772 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3773 : OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE>(H1, inv_jac_ptr) {}
3774};
3775
3776template <int DERIVARIVE = 1>
3778 : public OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE> {
3780 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3781 : OpSetInvJacSpaceForFaceImpl<3, DERIVARIVE>(L2, inv_jac_ptr) {}
3782};
3783
3784/**
3785 * \brief Transform local reference derivatives of shape function to
3786 global derivatives for face
3787
3788 * \ingroup mofem_forces_and_sources_tri_element
3789 */
3790template <int DIM> struct OpSetInvJacHcurlFaceImpl;
3791
3792template <>
3795
3796 OpSetInvJacHcurlFaceImpl(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3798 invJacPtr(inv_jac_ptr) {}
3799
3800 MoFEMErrorCode doWork(int side, EntityType type,
3802
3803protected:
3804 boost::shared_ptr<MatrixDouble> invJacPtr;
3806};
3807
3808template <>
3810 using OpSetInvJacHcurlFaceImpl<2>::OpSetInvJacHcurlFaceImpl;
3811 MoFEMErrorCode doWork(int side, EntityType type,
3813};
3814
3817
3818/**
3819 * @brief Make Hdiv space from Hcurl space in 2d
3820 * @ingroup mofem_forces_and_sources_tri_element
3821 */
3831
3832/** \brief Transform Hcurl base fluxes from reference element to physical
3833 * triangle
3834 */
3836
3837/** \brief Apply contravariant (Piola) transfer to Hdiv space on face
3838 *
3839 * Covariant Piola transformation
3840 \f[
3841 \psi_i|_t = J^{-1}_{ij}\hat{\psi}_j\\
3842 \left.\frac{\partial \psi_i}{\partial \xi_j}\right|_t
3843 = J^{-1}_{ik}\frac{\partial \hat{\psi}_k}{\partial \xi_j}
3844 \f]
3845
3846 */
3847template <>
3850
3852 boost::shared_ptr<MatrixDouble> inv_jac_ptr);
3853
3854 MoFEMErrorCode doWork(int side, EntityType type,
3856
3857private:
3858 boost::shared_ptr<MatrixDouble> invJacPtr;
3859
3862};
3863
3866
3867/** \brief Apply contravariant (Piola) transfer to Hdiv space on face
3868 *
3869 * \note Hdiv space is generated by Hcurl space in 2d.
3870 *
3871 * Contravariant Piola transformation
3872 * \f[
3873 * \psi_i|_t = \frac{1}{\textrm{det}(J)}J_{ij}\hat{\psi}_j\\
3874 * \left.\frac{\partial \psi_i}{\partial \xi_j}\right|_t
3875 * =
3876 * \frac{1}{\textrm{det}(J)}J_{ik}\frac{\partial \hat{\psi}_k}{\partial \xi_j}
3877 * \f]
3878 *
3879 * \ingroup mofem_forces_and_sources
3880 *
3881 */
3883
3884template <>
3887
3889 boost::shared_ptr<MatrixDouble> jac_ptr)
3891 jacPtr(jac_ptr) {}
3892
3893 MoFEMErrorCode doWork(int side, EntityType type,
3895
3896protected:
3897 boost::shared_ptr<MatrixDouble> jacPtr;
3900};
3901
3902template <>
3906 2>::OpSetContravariantPiolaTransformOnFace2DImpl;
3907
3908 MoFEMErrorCode doWork(int side, EntityType type,
3910};
3911
3916
3917/**@}*/
3918
3919/** \name Operators for edges */
3920
3921/**@{*/
3922
3932
3933/**
3934 * @deprecated Name is deprecated and this is added for backward compatibility
3935 */
3938
3939/**@}*/
3940
3941/** \name Operator for fat prisms */
3942
3943/**@{*/
3944
3945/**
3946 * @brief Operator for fat prism element updating integration weights in the
3947 * volume.
3948 *
3949 * Jacobian on the distorted element is nonconstant. This operator updates
3950 * integration weight on prism to take into account nonconstat jacobian.
3951 *
3952 * \f[
3953 * W_i = w_i \left( \frac{1}{2V} \left\| \frac{\partial \mathbf{x}}{\partial
3954 * \pmb\xi} \right\| \right)
3955 * \f]
3956 * where \f$w_i\f$ is integration weight at integration point \f$i\f$,
3957 * \f$\mathbf{x}\f$ is physical coordinate, and \f$\pmb\xi\f$ is reference
3958 * element coordinate.
3959 *
3960 */
3970
3971/** \brief Calculate inverse of jacobian for face element
3972
3973 It is assumed that face element is XY plane. Applied
3974 only for 2d problems.
3975
3976 FIXME Generalize function for arbitrary face orientation in 3d space
3977 FIXME Calculate to Jacobins for two faces
3978
3979 \ingroup mofem_forces_and_sources_prism_element
3980
3981*/
3984
3985 OpCalculateInvJacForFatPrism(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
3987 invJacPtr(inv_jac_ptr), invJac(*invJacPtr) {}
3988
3992 MoFEMErrorCode doWork(int side, EntityType type,
3994
3995private:
3996 const boost::shared_ptr<MatrixDouble> invJacPtr;
3998};
3999
4000/** \brief Transform local reference derivatives of shape functions to global
4001derivatives
4002
4003FIXME Generalize to curved shapes
4004FIXME Generalize to case that top and bottom face has different shape
4005
4006\ingroup mofem_forces_and_sources_prism_element
4007
4008*/
4011
4012 OpSetInvJacH1ForFatPrism(boost::shared_ptr<MatrixDouble> inv_jac_ptr)
4014 invJacPtr(inv_jac_ptr), invJac(*invJacPtr) {}
4015
4019
4020 MoFEMErrorCode doWork(int side, EntityType type,
4022
4023private:
4024 const boost::shared_ptr<MatrixDouble> invJacPtr;
4027};
4028
4029// Flat prism
4030
4031/** \brief Calculate inverse of jacobian for face element
4032
4033 It is assumed that face element is XY plane. Applied
4034 only for 2d problems.
4035
4036 FIXME Generalize function for arbitrary face orientation in 3d space
4037 FIXME Calculate to Jacobins for two faces
4038
4039 \ingroup mofem_forces_and_sources_prism_element
4040
4041*/
4054
4055/** \brief Transform local reference derivatives of shape functions to global
4056derivatives
4057
4058FIXME Generalize to curved shapes
4059FIXME Generalize to case that top and bottom face has different shape
4060
4061\ingroup mofem_forces_and_sources_prism_element
4062
4063*/
4078
4079/**@}*/
4080
4081/** \name Operation on matrices at integration points */
4082
4083/**@{*/
4084
4085/**
4086 * @brief Operator for inverting matrices at integration points
4087 *
4088 * This template structure computes the inverse of square matrices and their
4089 * determinants at integration points. It's commonly used in finite element
4090 * methods for coordinate transformations, constitutive relations, and other
4091 * matrix operations requiring matrix inversion.
4092 *
4093 * @tparam DIM Dimension of the square matrix to be inverted
4094 */
4095template <int DIM>
4097
4098 /**
4099 * @brief Constructor for matrix inversion operator
4100 *
4101 * @param in_ptr Shared pointer to input matrix to be inverted
4102 * @param det_ptr Shared pointer to vector for storing matrix determinants
4103 * @param out_ptr Shared pointer to output matrix for storing inverted matrices
4104 */
4105 OpInvertMatrix(boost::shared_ptr<MatrixDouble> in_ptr,
4106 boost::shared_ptr<VectorDouble> det_ptr,
4107 boost::shared_ptr<MatrixDouble> out_ptr)
4109 outPtr(out_ptr), detPtr(det_ptr) {}
4110
4111 MoFEMErrorCode doWork(int side, EntityType type,
4113
4114private:
4115 boost::shared_ptr<MatrixDouble> inPtr;
4116 boost::shared_ptr<MatrixDouble> outPtr;
4117 boost::shared_ptr<VectorDouble> detPtr;
4118};
4119
4120template <int DIM>
4124
4125 if (!inPtr)
4126 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
4127 "Pointer for inPtr matrix not allocated");
4128 if (!detPtr)
4129 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
4130 "Pointer for detPtr matrix not allocated");
4131
4132 const auto nb_rows = inPtr->size1();
4133 const auto nb_integration_pts = inPtr->size2();
4134
4135 // Calculate determinant
4136 {
4137 detPtr->resize(nb_integration_pts, false);
4138 auto t_in = getFTensor2FromMat<DIM, DIM>(*inPtr);
4139 auto t_det = getFTensor0FromVec(*detPtr);
4140 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
4141 t_det = determinantTensor(t_in);
4142 ++t_in;
4143 ++t_det;
4144 }
4145 }
4146
4147 // Invert jacobian
4148 if (outPtr) {
4149 outPtr->resize(nb_rows, nb_integration_pts, false);
4150 auto t_in = getFTensor2FromMat<DIM, DIM>(*inPtr);
4151 auto t_out = getFTensor2FromMat<DIM, DIM>(*outPtr);
4152 auto t_det = getFTensor0FromVec(*detPtr);
4153 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
4154 CHKERR invertTensor(t_in, t_det, t_out);
4155 ++t_in;
4156 ++t_out;
4157 ++t_det;
4158 }
4159 }
4160
4162}
4163
4164/**@}*/
4165
4166/**
4167 * @brief Operator for calculating the trace of matrices at integration points
4168 *
4169 * This template structure computes the trace (sum of diagonal elements) of
4170 * square matrices at integration points. The trace is commonly used in
4171 * mechanics for calculating volumetric strain, pressure, and other scalar
4172 * quantities derived from tensors.
4173 *
4174 * @tparam DIM Dimension of the square matrix
4175 *
4176 * @ingroup mofem_forces_and_sources
4177 */
4178template <int DIM>
4180
4181 /**
4182 * @brief Constructor for matrix trace calculation operator
4183 *
4184 * @param in_ptr Shared pointer to input matrix for trace calculation
4185 * @param out_ptr Shared pointer to output vector for storing calculated traces
4186 */
4187 OpCalculateTraceFromMat(boost::shared_ptr<MatrixDouble> in_ptr,
4188 boost::shared_ptr<VectorDouble> out_ptr)
4190 outPtr(out_ptr) {}
4191
4192 MoFEMErrorCode doWork(int side, EntityType type,
4194
4195private:
4197 boost::shared_ptr<MatrixDouble> inPtr;
4198 boost::shared_ptr<VectorDouble> outPtr;
4199};
4200
4201template <int DIM>
4206
4207 if (!inPtr)
4208 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
4209 "Pointer for inPtr matrix not allocated");
4210
4211 const auto nb_integration_pts = inPtr->size2();
4212 // Invert jacobian
4213 if (outPtr) {
4214 outPtr->resize(nb_integration_pts, false);
4215 auto t_in = getFTensor2FromMat<DIM, DIM>(*inPtr);
4216 auto t_out = getFTensor0FromVec(*outPtr);
4217
4218 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
4219 t_out = t_in(i, i);
4220 ++t_in;
4221 ++t_out;
4222 }
4223 }
4224
4226}
4227
4228/**@}*/
4229
4230/** \brief Calculates the trace of an input matrix
4231
4232\ingroup mofem_forces_and_sources
4233
4234*/
4235
4236template <int DIM>
4239
4240 OpCalculateTraceFromSymmMat(boost::shared_ptr<MatrixDouble> in_ptr,
4241 boost::shared_ptr<VectorDouble> out_ptr)
4243 outPtr(out_ptr) {}
4244
4245 MoFEMErrorCode doWork(int side, EntityType type,
4247
4248private:
4250 boost::shared_ptr<MatrixDouble> inPtr;
4251 boost::shared_ptr<VectorDouble> outPtr;
4252};
4253
4254template <int DIM>
4259
4260 if (!inPtr)
4261 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
4262 "Pointer for inPtr matrix not allocated");
4263
4264 const auto nb_integration_pts = inPtr->size2();
4265 // Invert jacobian
4266 if (outPtr) {
4267 outPtr->resize(nb_integration_pts, false);
4268 auto t_in = getFTensor2SymmetricFromMat<DIM>(*inPtr);
4269 auto t_out = getFTensor0FromVec(*outPtr);
4270
4271 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
4272 t_out = t_in(i, i);
4273 ++t_in;
4274 ++t_out;
4275 }
4276 }
4277
4279}
4280
4281} // namespace MoFEM
4282
4283#endif // __USER_DATA_OPERATORS_HPP__
4284
4285/**
4286 * \defgroup mofem_forces_and_sources_user_data_operators Users Operators
4287 *
4288 * \brief Classes and functions used to evaluate fields at integration pts,
4289 *jacobians, etc..
4290 *
4291 * \ingroup mofem_forces_and_sources
4292 **/
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr int SPACE_DIM
Tensor1< T, Tensor_Dim > normalize()
static const char *const CoordinateTypesNames[]
Coordinate system names.
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
FieldSpace
approximation spaces
Definition definitions.h:82
@ L2
field with C-1 continuity
Definition definitions.h:88
@ H1
continuous field
Definition definitions.h:85
@ NOSPACE
Definition definitions.h:83
@ HCURL
field with continuous tangents
Definition definitions.h:86
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
CoordinateTypes
Coodinate system.
@ CYLINDRICAL
@ POLAR
@ CARTESIAN
@ SPHERICAL
#define DEPRECATED
Definition definitions.h:17
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
constexpr int BASE_DIM
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'J', DIM1 > J
Definition level_set.cpp:30
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VecAllocator< double > DoubleAllocator
Definition Types.hpp:62
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
static MoFEMErrorCode invertTensor(FTensor::Tensor2< T1, DIM, DIM > &t, T2 &det, FTensor::Tensor2< T3, DIM, DIM > &inv_t)
static auto determinantTensor(FTensor::Tensor2< T, DIM, DIM > &t)
Calculate the determinant of a tensor of rank DIM.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
constexpr IntegrationType I
constexpr AssemblyType A
double scale
Definition plastic.cpp:123
constexpr auto field_name
multi_index_container< DofsSideMapData, indexed_by< ordered_non_unique< tag< TypeSide_mi_tag >, composite_key< DofsSideMapData, member< DofsSideMapData, EntityType, &DofsSideMapData::type >, member< DofsSideMapData, int, &DofsSideMapData::side > > >, ordered_unique< tag< EntDofIdx_mi_tag >, member< DofsSideMapData, int, &DofsSideMapData::dof > > > > DofsSideMap
Map entity stype and side to element/entity dof index.
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
const VectorFieldEntities & getFieldEntities() const
Get field entities (const version)
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
Get DOF values on entity.
FTensor::Tensor3< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 *Tensor_Dim2 >, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2 > getFTensor3Diff2N(FieldApproximationBase base)
Get second derivatives of base functions for Hvec space.
const VectorInt & getLocalIndices() const
Get local indices of degrees of freedom on entity.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1FieldData()
Return FTensor of rank 1, i.e. vector from field data coefficients.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FieldData()
Return scalar files as a FTensor of rank 0.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
FTensor::Tensor2_symmetric< FTensor::PackPtr< double *,(Tensor_Dim *(Tensor_Dim+1))/2 >, Tensor_Dim > getFTensor2SymmetricFieldData()
Return symmetric FTensor rank 2, i.e. matrix from field data coefficients.
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2FieldData()
Return FTensor rank 2, i.e. matrix from field data coefficients.
const VectorInt & getIndices() const
Get global indices of degrees of freedom on entity.
FTensor::Tensor3< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 *Tensor_Dim2 >, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2 > getFTensor3DiffN(FieldApproximationBase base)
Get derivatives of base functions for tonsorial Hdiv space.
EntityType getFEType() const
Get dimension of finite element.
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
EntityHandle getSideEntity(const int side_number, const EntityType type)
Get the side entity.
@ OPROW
operator doWork function is executed on FE rows
@ OPSPACE
operator do Work is execute on space data
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
structure to get information from mofem into EntitiesFieldData
Calculate divergence of tonsorial field using vectorial base.
OpCalculateBrokenHVecTensorDivergence(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, SmartPetscObj< Vec > data_vec, EntityType broken_type, boost::shared_ptr< Range > broken_range_ptr=nullptr, boost::shared_ptr< double > scale_ptr=nullptr, const EntityType zero_type=MBEDGE)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Get tensor field for H-div approximation.
boost::shared_ptr< MatrixDouble > dataPtr
OpCalculateBrokenHVecTensorField(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, SmartPetscObj< Vec > data_vec, EntityType broken_type, boost::shared_ptr< Range > broken_range_ptr=nullptr, boost::shared_ptr< double > scale_ptr=nullptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Calculate values of vector field at integration points.
Calculate divergence of vector field at integration points.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateDivergenceVectorFieldValues(const std::string field_name, boost::shared_ptr< VectorDouble > data_ptr, const EntityType zero_type=MBVERTEX)
Constructor for vector field divergence calculation operator.
Calculate tenor field using tensor base, i.e. Hdiv/Hcurl.
VectorDouble dotVector
Keeps temporary values of time derivatives.
OpCalculateHTensorTensorField(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, boost::shared_ptr< double > scale_ptr, SmartPetscObj< Vec > data_vec=SmartPetscObj< Vec >(), const EntityType zero_type=MBEDGE, const int zero_side=0)
boost::shared_ptr< MatrixDouble > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateHTensorTensorField(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Calculate divergence of tonsorial field using vectorial base.
boost::shared_ptr< MatrixDouble > dataPtr
VectorDouble dotVector
Keeps temporary values of time derivatives.
OpCalculateHVecTensorDivergence(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, SmartPetscObj< Vec > data_vec, const EntityType zero_type=MBEDGE, const int zero_side=0)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateHVecTensorDivergence(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
OpCalculateHVecTensorField(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, boost::shared_ptr< double > scale_ptr, SmartPetscObj< Vec > data_vec=SmartPetscObj< Vec >(), const EntityType zero_type=MBEDGE, const int zero_side=0)
boost::shared_ptr< MatrixDouble > dataPtr
OpCalculateHVecTensorField(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
VectorDouble dotVector
Keeps temporary values of time derivatives.
boost::shared_ptr< double > scalePtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateHVecTensorGradient(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
OpCalculateHVecTensorGradient(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Calculate gradient of tensor field.
Calculate trace of vector (Hdiv/Hcurl) space.
boost::shared_ptr< MatrixDouble > dataPtr
boost::shared_ptr< double > scalePtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpCalculateHVecTensorTrace(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
FTensor::Index< 'j', Tensor_Dim > j
FTensor::Index< 'i', Tensor_Dim > i
OpCalculateHVecTensorTrace(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, boost::shared_ptr< double > scale_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
OpCalculateHVecVectorFieldDot(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Get vector field for H-div approximation.
OpCalculateHVecVectorField_General(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, SmartPetscObj< Vec > data_vec, const EntityType zero_type=MBEDGE, const int zero_side=0)
OpCalculateHVecVectorField_General(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Get vector field for H-div approximation.
Get vector field for H-div approximation.
Calculate gradient of vector field.
boost::shared_ptr< MatrixDouble > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateHVecVectorGradient(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Calculate gradient of vector field.
boost::shared_ptr< MatrixDouble > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateHVecVectorHessian(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
boost::shared_ptr< MatrixDouble > dataPtr
boost::shared_ptr< MatrixDouble > dataPtr
boost::shared_ptr< MatrixDouble > dataPtr
Calculate curl of vector field.
Calculate divergence of vector field dot.
boost::shared_ptr< VectorDouble > dataPtr
OpCalculateHdivVectorDivergenceDot(const std::string field_name, boost::shared_ptr< VectorDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Calculate divergence of vector field.
boost::shared_ptr< VectorDouble > dataPtr
OpCalculateHdivVectorDivergence(const std::string field_name, boost::shared_ptr< VectorDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Calculate inverse of jacobian for face element.
OpCalculateInvJacForFatPrism(MatrixDouble &inv_jac)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateInvJacForFatPrism(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
const boost::shared_ptr< MatrixDouble > invJacPtr
Calculate inverse of jacobian for face element.
OpCalculateInvJacForFlatPrism(MatrixDouble &inv_jac_f3)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Evaluate field gradient values for scalar field, i.e. gradient is tensor rank 1 (vector)
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Evaluate field gradient values for scalar field, i.e. gradient is tensor rank 1 (vector),...
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate gradient values of scalar field at integration points
Calculate scalar field values from PETSc vector at integration points.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateScalarFieldValuesFromPetscVecImpl(const std::string field_name, boost::shared_ptr< VectorDouble > data_ptr, const EntityType zero_at_type=MBVERTEX)
Constructor for PETSc vector-based scalar field calculation.
Scalar field values at integration points.
OpCalculateScalarFieldValues_General(const std::string field_name, boost::shared_ptr< ublas::vector< T, A > > data_ptr, SmartPetscObj< Vec > data_vec, const EntityType zero_type=MBVERTEX)
Constructor with PETSc vector support.
boost::shared_ptr< ublas::vector< T, A > > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate values of scalar field at integration points
OpCalculateScalarFieldValues_General(const std::string field_name, boost::shared_ptr< ublas::vector< T, A > > data_ptr, const EntityType zero_type=MBVERTEX)
Constructor for scalar field values calculation operator.
Specialization for double precision scalar field values calculation.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate values of scalar field at integration points
Get time direvarive values at integration pts for tensor field rank 2, i.e. matrix field.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateTensor2FieldValuesDot(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_at_type=MBVERTEX)
EntityType zeroAtType
Zero values at Gauss point at this type.
VectorDouble dotVector
Keeps temporary values of time derivatives.
boost::shared_ptr< MatrixDouble > dataPtr
Data computed into this matrix.
OpCalculateTensor2FieldValues_General(const std::string field_name, boost::shared_ptr< ublas::matrix< double, ublas::row_major, DoubleAllocator > > data_ptr, SmartPetscObj< Vec > data_vec, const EntityType zero_type=MBVERTEX)
OpCalculateTensor2FieldValues_General(const std::string field_name, boost::shared_ptr< ublas::matrix< double, ublas::row_major, DoubleAllocator > > data_ptr, const EntityType zero_type=MBVERTEX)
Calculate field values for tenor field rank 2.
OpCalculateTensor2FieldValues_General(const std::string field_name, boost::shared_ptr< ublas::matrix< T, L, A > > data_ptr, const EntityType zero_type=MBVERTEX)
boost::shared_ptr< ublas::matrix< T, L, A > > dataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Get values at integration pts for tensor field rank 2, i.e. matrix field.
OpCalculateTensor2SymmetricFieldGradient_General(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBVERTEX)
Evaluate field gradient values for symmetric 2nd order tensor field, i.e. gradient is tensor rank 3.
OpCalculateTensor2SymmetricFieldGradient_General(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBVERTEX)
Get field gradients at integration pts for symmetric tensorial field rank 2.
OpCalculateTensor2SymmetricFieldGradient(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBVERTEX)
Calculate symmetric tensor field rates ant integratio pts.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateTensor2SymmetricFieldValuesDot(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
Calculate symmetric tensor field values at integration pts.
OpCalculateTensor2SymmetricFieldValues(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, SmartPetscObj< Vec > data_vec, const EntityType zero_type=MBEDGE, const int zero_side=0)
OpCalculateTensor2SymmetricFieldValues(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBEDGE, const int zero_side=0)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Operator for calculating the trace of matrices at integration points.
boost::shared_ptr< MatrixDouble > inPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< VectorDouble > outPtr
OpCalculateTraceFromMat(boost::shared_ptr< MatrixDouble > in_ptr, boost::shared_ptr< VectorDouble > out_ptr)
Constructor for matrix trace calculation operator.
Calculates the trace of an input matrix.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< VectorDouble > outPtr
boost::shared_ptr< MatrixDouble > inPtr
OpCalculateTraceFromSymmMat(boost::shared_ptr< MatrixDouble > in_ptr, boost::shared_ptr< VectorDouble > out_ptr)
Get field gradients time derivative at integration pts for scalar field rank 0, i....
boost::shared_ptr< MatrixDouble > dataPtr
Data computed into this matrix.
OpCalculateVectorFieldGradientDot(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_at_type=MBVERTEX)
EntityType zeroAtType
Zero values at Gauss point at this type.
VectorDouble dotVector
Keeps temporary values of time derivatives.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Evaluate field gradient values for vector field, i.e. gradient is tensor rank 2.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate values of vector field at integration points
Approximate field values for given petsc vector.
OpCalculateVectorFieldValuesFromPetscVecImpl(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_at_type=MBVERTEX, bool throw_error=true)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpCalculateVectorFieldValues_General(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, const EntityType zero_type=MBVERTEX)
OpCalculateVectorFieldValues_General(const std::string field_name, boost::shared_ptr< MatrixDouble > data_ptr, SmartPetscObj< Vec > data_vec, const EntityType zero_type=MBVERTEX)
Calculate field values for tensor field rank 1, i.e. vector field.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate values of vector field at integration points
OpCalculateVectorFieldValues_General(const std::string field_name, boost::shared_ptr< ublas::matrix< T, L, A > > data_ptr, const EntityType zero_type=MBVERTEX)
Constructor for vector field values calculation operator.
boost::shared_ptr< ublas::matrix< T, L, A > > dataPtr
Specialization for double precision vector field values calculation.
Operator for inverting matrices at integration points.
boost::shared_ptr< VectorDouble > detPtr
OpInvertMatrix(boost::shared_ptr< MatrixDouble > in_ptr, boost::shared_ptr< VectorDouble > det_ptr, boost::shared_ptr< MatrixDouble > out_ptr)
Constructor for matrix inversion operator.
boost::shared_ptr< MatrixDouble > outPtr
boost::shared_ptr< MatrixDouble > inPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Make Hdiv space from Hcurl space in 2d.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Operator for fat prism element updating integration weights in the volume.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Operator for scaling matrix values by a scalar factor.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< MatrixDouble > outMat
DEPRECATED OpScaleMatrix(const std::string field_name, const double scale, boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat)
boost::shared_ptr< double > scalePtr
OpScaleMatrix(boost::shared_ptr< double > scale_ptr, boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat)
Constructor for matrix scaling operator.
boost::shared_ptr< MatrixDouble > inMat
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpSetContravariantPiolaTransformOnEdge2D(const FieldSpace space=HCURL)
OpSetContravariantPiolaTransformOnFace2DImpl(boost::shared_ptr< MatrixDouble > jac_ptr)
Apply contravariant (Piola) transfer to Hdiv space on face.
Apply contravariant (Piola) transfer to Hdiv space on face.
Transform Hcurl base fluxes from reference element to physical triangle.
OpSetInvJacH1ForFaceEmbeddedIn3DSpace(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
OpSetInvJacH1ForFace(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
Transform local reference derivatives of shape functions to global derivatives.
const boost::shared_ptr< MatrixDouble > invJacPtr
OpSetInvJacH1ForFatPrism(MatrixDouble &inv_jac)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpSetInvJacH1ForFatPrism(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
Transform local reference derivatives of shape functions to global derivatives.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpSetInvJacH1ForFlatPrism(MatrixDouble &inv_jac_f3)
boost::shared_ptr< MatrixDouble > invJacPtr
OpSetInvJacHcurlFaceImpl(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
Transform local reference derivatives of shape function to global derivatives for face.
OpSetInvJacL2ForFaceEmbeddedIn3DSpace(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
OpSetInvJacL2ForFace(boost::shared_ptr< MatrixDouble > inv_jac_ptr)
Transform local reference derivatives of shape functions to global derivatives.
MoFEMErrorCode applyTransform(MatrixDouble &diff_n)
Apply transformation to the input matrix.
OpSetInvJacToScalarBasesBasic(FieldSpace space, boost::shared_ptr< MatrixDouble > inv_jac_ptr)
boost::shared_ptr< MatrixDouble > invJacPtr
Operator for symmetrizing tensor fields.
FTensor::Index< 'i', DIM > i
boost::shared_ptr< MatrixDouble > outMat
OpSymmetrizeTensor(boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat)
Constructor for tensor symmetrization operator.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
FTensor::Index< 'j', DIM > j
DEPRECATED OpSymmetrizeTensor(const std::string field_name, boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat)
boost::shared_ptr< MatrixDouble > inMat
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< MatrixDouble > dMat
OpTensorTimesSymmetricTensor(boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat, boost::shared_ptr< MatrixDouble > d_mat)
boost::shared_ptr< MatrixDouble > inMat
DEPRECATED OpTensorTimesSymmetricTensor(const std::string field_name, boost::shared_ptr< MatrixDouble > in_mat, boost::shared_ptr< MatrixDouble > out_mat, boost::shared_ptr< MatrixDouble > d_mat)
boost::shared_ptr< MatrixDouble > outMat
static constexpr Switches CtxSetX
Solution vector switch.
static constexpr Switches CtxSetX_TT
Second time derivative switch.
std::bitset< 8 > Switches
Bitset type for context switches.
static constexpr Switches CtxSetX_T
First time derivative switch.
@ CTX_SET_X_T
Time derivative X_t is set.
@ CTX_SET_DX
Solution increment DX is set.
@ CTX_SET_X
Solution vector X is set.
@ CTX_SET_X_TT
Second time derivative X_tt is set.
intrusive_ptr for managing petsc objects