v0.15.0
Loading...
Searching...
No Matches
Public Types | Public Member Functions | Private Member Functions | Private Attributes | List of all members
CGGUserPolynomialBase Struct Reference

CGG User Polynomial Base. More...

#include "users_modules/eshelbian_plasticity/src/CGGUserPolynomialBase.hpp"

Inheritance diagram for CGGUserPolynomialBase:
[legend]
Collaboration diagram for CGGUserPolynomialBase:
[legend]

Public Types

using CachePhi = boost::tuple< int, int, MatrixDouble, MatrixDouble >
 

Public Member Functions

 CGGUserPolynomialBase (boost::shared_ptr< CachePhi > cache_phi=nullptr)
 
 ~CGGUserPolynomialBase ()=default
 
MoFEMErrorCode query_interface (boost::typeindex::type_index type_index, BaseFunctionUnknownInterface **iface) const
 
MoFEMErrorCode getValue (MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
 

Private Member Functions

MoFEMErrorCode getValueHdivForCGGBubble (MatrixDouble &pts)
 

Private Attributes

MatrixDouble shapeFun
 
boost::shared_ptr< CachePhi > cachePhiPtr
 

Detailed Description

CGG User Polynomial Base.

Examples
mofem/users_modules/eshelbian_plasticity/src/impl/EshelbianPlasticity.cpp.

Definition at line 6 of file CGGUserPolynomialBase.hpp.

Member Typedef Documentation

◆ CachePhi

using CGGUserPolynomialBase::CachePhi = boost::tuple<int, int, MatrixDouble, MatrixDouble>

Definition at line 8 of file CGGUserPolynomialBase.hpp.

Constructor & Destructor Documentation

◆ CGGUserPolynomialBase()

CGGUserPolynomialBase::CGGUserPolynomialBase ( boost::shared_ptr< CachePhi >  cache_phi = nullptr)

Definition at line 7 of file CGGUserPolynomialBase.cpp.

9 : TetPolynomialBase(), cachePhiPtr(cache_phi_otr) {}
CGGUserPolynomialBase::cachePhiPtr
boost::shared_ptr< CachePhi > cachePhiPtr
Definition CGGUserPolynomialBase.hpp:20

◆ ~CGGUserPolynomialBase()

CGGUserPolynomialBase::~CGGUserPolynomialBase ( )
default

Member Function Documentation

◆ getValue()

MoFEMErrorCode CGGUserPolynomialBase::getValue ( MatrixDouble &  pts,
boost::shared_ptr< BaseFunctionCtx >  ctx_ptr 
)

Definition at line 19 of file CGGUserPolynomialBase.cpp.

20 {
21 MoFEMFunctionBeginHot;
22
23 this->cTx = ctx_ptr->getInterface<EntPolynomialBaseCtx>();
24
25 int nb_gauss_pts = pts.size2();
26 if (!nb_gauss_pts) {
27 MoFEMFunctionReturnHot(0);
28 }
29
30 if (pts.size1() < 3) {
31 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
32 "Wrong dimension of pts, should be at least 3 rows with "
33 "coordinates");
34 }
35
36 const auto base = this->cTx->bAse;
37 EntitiesFieldData &data = this->cTx->dAta;
38
39 switch (this->cTx->sPace) {
40 case HDIV:
41 CHKERR getValueHdivForCGGBubble(pts);
42 break;
43 case L2:
44 data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 4, false);
45 CHKERR Tools::shapeFunMBTET(
46 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
47 &pts(0, 0), &pts(1, 0), &pts(2, 0), nb_gauss_pts);
48 data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(4, 3, false);
49 std::copy(Tools::diffShapeFunMBTET.begin(), Tools::diffShapeFunMBTET.end(),
50 data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin());
51 this->cTx->basePolynomialsType0 = Legendre_polynomials;
52 CHKERR getValueL2AinsworthBase(pts);
53 break;
54 default:
55 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
56 }
57
58 MoFEMFunctionReturnHot(0);
59}
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
L2
@ L2
field with C-1 continuity
Definition definitions.h:88
HDIV
@ HDIV
field with continuous normal traction
Definition definitions.h:87
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:456
MoFEM::Legendre_polynomials
PetscErrorCode Legendre_polynomials(int p, double s, double *diff_s, double *L, double *diffL, int dim)
Calculate Legendre approximation basis.
Definition LegendrePolynomial.cpp:132
CGGUserPolynomialBase::getValueHdivForCGGBubble
MoFEMErrorCode getValueHdivForCGGBubble(MatrixDouble &pts)
Definition CGGUserPolynomialBase.cpp:62

◆ getValueHdivForCGGBubble()

MoFEMErrorCode CGGUserPolynomialBase::getValueHdivForCGGBubble ( MatrixDouble &  pts)
private

number of integration points

Definition at line 62 of file CGGUserPolynomialBase.cpp.

62 {
63 MoFEMFunctionBegin;
64
65 constexpr bool calculate_diff = false;
66 constexpr bool debug = false;
67
68 // This should be used only in case USER_BASE is selected
69 if (this->cTx->bAse != USER_BASE) {
70 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
71 "Wrong base, should be USER_BASE");
72 }
73 // get access to data structures on element
74 EntitiesFieldData &data = this->cTx->dAta;
75 // Get approximation order on element. Note that bubble functions are only
76 // on tetrahedron.
77 const int order = data.dataOnEntities[MBTET][0].getOrder();
78 /// number of integration points
79 const int nb_gauss_pts = pts.size2();
80
81 auto check_cache = [this](int order, int nb_gauss_pts) -> bool {
82 if (cachePhiPtr) {
83 return cachePhiPtr->get<0>() == order &&
84 cachePhiPtr->get<1>() == nb_gauss_pts;
85 } else {
86 return false;
87 }
88 };
89
90 if (check_cache(order, nb_gauss_pts)) {
91 auto &mat = cachePhiPtr->get<2>();
92 auto &phi = data.dataOnEntities[MBTET][0].getN(USER_BASE);
93 phi.resize(mat.size1(), mat.size2(), false);
94 noalias(phi) = mat;
95 if constexpr (calculate_diff) {
96 auto &diff_phi = data.dataOnEntities[MBTET][0].getDiffN(USER_BASE);
97 auto &diff_mat = cachePhiPtr->get<3>();
98 diff_phi.resize(diff_mat.size1(), diff_mat.size2(), false);
99 noalias(diff_phi) = diff_mat;
100 }
101 } else {
102 // calculate shape functions, i.e. barycentric coordinates
103 shapeFun.resize(nb_gauss_pts, 4, false);
104 CHKERR ShapeMBTET(&*shapeFun.data().begin(), &pts(0, 0), &pts(1, 0),
105 &pts(2, 0), nb_gauss_pts);
106 // derivatives of shape functions
107 double diff_shape_fun[12];
108 CHKERR ShapeDiffMBTET(diff_shape_fun);
109
110 const int nb_base_functions = NBVOLUMETET_CCG_BUBBLE(order);
111 // get base functions and set size
112 MatrixDouble &phi = data.dataOnEntities[MBTET][0].getN(USER_BASE);
113 phi.resize(nb_gauss_pts, 9 * nb_base_functions, false);
114 // finally calculate base functions
115 FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_phi(
116 &phi(0, 0), &phi(0, 1), &phi(0, 2),
117
118 &phi(0, 3), &phi(0, 4), &phi(0, 5),
119
120 &phi(0, 6), &phi(0, 7), &phi(0, 8));
121 CHKERR CGG_BubbleBase_MBTET(order, &shapeFun(0, 0), diff_shape_fun, t_phi,
122 nb_gauss_pts);
123
124 auto imagingray_fd_directive = [&](auto &diff_phi) {
125 MoFEMFunctionBegin;
126
127 MatrixComplexDouble temp_phi;
128 temp_phi.resize(nb_gauss_pts, 9 * nb_base_functions, false);
129
130 MatrixComplexDouble diff_shape_fun;
131 diff_shape_fun.resize(4, 3, false);
132 diff_shape_fun(0, 0) = diffN_MBTET0x;
133 diff_shape_fun(0, 1) = diffN_MBTET0y;
134 diff_shape_fun(0, 2) = diffN_MBTET0z;
135 diff_shape_fun(1, 0) = diffN_MBTET1x;
136 diff_shape_fun(1, 1) = diffN_MBTET1y;
137 diff_shape_fun(1, 2) = diffN_MBTET1z;
138 diff_shape_fun(2, 0) = diffN_MBTET2x;
139 diff_shape_fun(2, 1) = diffN_MBTET2y;
140 diff_shape_fun(2, 2) = diffN_MBTET2z;
141 diff_shape_fun(3, 0) = diffN_MBTET3x;
142 diff_shape_fun(3, 1) = diffN_MBTET3y;
143 diff_shape_fun(3, 2) = diffN_MBTET3z;
144
145 MatrixComplexDouble shape_fun;
146 shape_fun.resize(nb_gauss_pts, 4, false);
147 for (int dd = 0; dd < 3; ++dd) {
148 constexpr double eps = 1e-16;
149
150 for (auto gg = 0; gg < nb_gauss_pts; ++gg) {
151 std::array<std::complex<double>, 3> ksi{pts(0, gg), pts(1, gg),
152 pts(2, gg)};
153
154 ksi[dd] += std::complex<double>(0, eps);
155
156 shape_fun(gg, 0) = N_MBTET0(ksi[0], ksi[1], ksi[2]);
157 shape_fun(gg, 1) = N_MBTET1(ksi[0], ksi[1], ksi[2]);
158 shape_fun(gg, 2) = N_MBTET2(ksi[0], ksi[1], ksi[2]);
159 shape_fun(gg, 3) = N_MBTET3(ksi[0], ksi[1], ksi[2]);
160 }
161
162 auto t_complex_phi_in =
163 FTensor::Tensor2<FTensor::PackPtr<std::complex<double> *, 9>, 3, 3>{
164 &temp_phi(0, 0), &temp_phi(0, 1), &temp_phi(0, 2),
165 &temp_phi(0, 3), &temp_phi(0, 4), &temp_phi(0, 5),
166 &temp_phi(0, 6), &temp_phi(0, 7), &temp_phi(0, 8)};
167
168 CHKERR CGG_BubbleBase_MBTET(order, shape_fun.data().data(),
169 diff_shape_fun.data().data(),
170 t_complex_phi_in, nb_gauss_pts);
171
172 auto t_complex_phi_out =
173 FTensor::Tensor2<FTensor::PackPtr<std::complex<double> *, 9>, 3, 3>{
174 &temp_phi(0, 0), &temp_phi(0, 1), &temp_phi(0, 2),
175 &temp_phi(0, 3), &temp_phi(0, 4), &temp_phi(0, 5),
176 &temp_phi(0, 6), &temp_phi(0, 7), &temp_phi(0, 8)};
177 auto t_diff_phi = getFTensor3FromPtr<3, 3, 3>(diff_phi.data().data());
178 for (auto gg = 0; gg < nb_gauss_pts; ++gg) {
179 for (auto bb = 0; bb << nb_base_functions; ++bb) {
180 t_diff_phi(0, 0, dd) = std::imag(t_complex_phi_out(0, 0)) / eps;
181 t_diff_phi(0, 1, dd) = std::imag(t_complex_phi_out(0, 1)) / eps;
182 t_diff_phi(0, 2, dd) = std::imag(t_complex_phi_out(0, 2)) / eps;
183 t_diff_phi(1, 0, dd) = std::imag(t_complex_phi_out(1, 0)) / eps;
184 t_diff_phi(1, 1, dd) = std::imag(t_complex_phi_out(1, 1)) / eps;
185 t_diff_phi(1, 2, dd) = std::imag(t_complex_phi_out(1, 2)) / eps;
186 t_diff_phi(2, 0, dd) = std::imag(t_complex_phi_out(2, 0)) / eps;
187 t_diff_phi(2, 1, dd) = std::imag(t_complex_phi_out(2, 1)) / eps;
188 t_diff_phi(2, 2, dd) = std::imag(t_complex_phi_out(2, 2)) / eps;
189 if (debug) {
190 MOFEM_LOG("SELF", Sev::noisy)
191 << "Gauss pt " << gg << " direction " << dd
192 << " diff_phi: " << t_diff_phi(0, 0, dd) << ", "
193 << t_diff_phi(0, 1, dd) << ", " << t_diff_phi(0, 2, dd)
194 << "; " << t_diff_phi(1, 0, dd) << ", "
195 << t_diff_phi(1, 1, dd) << ", " << t_diff_phi(1, 2, dd)
196 << "; " << t_diff_phi(2, 0, dd) << ", "
197 << t_diff_phi(2, 1, dd) << ", " << t_diff_phi(2, 2, dd);
198 }
199 ++t_complex_phi_out;
200 ++t_diff_phi;
201 }
202 }
203
204 if (debug) {
205 auto t_phi = FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3>{
206 &phi(0, 0), &phi(0, 1), &phi(0, 2), &phi(0, 3), &phi(0, 4),
207 &phi(0, 5), &phi(0, 6), &phi(0, 7), &phi(0, 8)};
208 auto t_complex_phi_out =
209 FTensor::Tensor2<FTensor::PackPtr<std::complex<double> *, 9>, 3,
210 3>{
211 &temp_phi(0, 0), &temp_phi(0, 1), &temp_phi(0, 2),
212 &temp_phi(0, 3), &temp_phi(0, 4), &temp_phi(0, 5),
213 &temp_phi(0, 6), &temp_phi(0, 7), &temp_phi(0, 8)};
214
215 for (auto gg = 0; gg < nb_gauss_pts; ++gg) {
216 FTensor::Index<'i', 3> i;
217 FTensor::Index<'j', 3> j;
218 FTensor::Tensor2<std::complex<double>, 3, 3> t_diff;
219 for (auto bb = 0; bb << nb_base_functions; ++bb) {
220 t_diff(i, j) = t_complex_phi_out(i, j) - t_phi(i, j);
221 auto nrm2_complex = t_diff(i, j) * t_diff(i, j);
222 if (std::abs(nrm2_complex) > 1e-12) {
223 MOFEM_LOG("SELF", Sev::warning)
224 << "Difference detected at gauss pt " << gg
225 << " for direction " << dd
226 << " nrm2 = " << std::abs(nrm2_complex);
227 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
228 "Difference detected between complex step "
229 "differentiation and real valued differentiation");
230 }
231 ++t_phi;
232 ++t_complex_phi_out;
233 }
234 }
235 }
236
237
238
239 }
240
241 MoFEMFunctionReturn(0);
242 };
243
244 MatrixDouble &diff_phi = data.dataOnEntities[MBTET][0].getDiffN(USER_BASE);
245 if constexpr (calculate_diff) {
246 diff_phi.resize(nb_gauss_pts, 27 * nb_base_functions, false);
247 CHKERR imagingray_fd_directive(diff_phi);
248 }
249
250 if (cachePhiPtr) {
251 cachePhiPtr->get<0>() = order;
252 cachePhiPtr->get<1>() = nb_gauss_pts;
253 cachePhiPtr->get<2>().resize(phi.size1(), phi.size2(), false);
254 noalias(cachePhiPtr->get<2>()) = phi;
255 if constexpr (calculate_diff) {
256 cachePhiPtr->get<3>().resize(diff_phi.size1(), diff_phi.size2(), false);
257 noalias(cachePhiPtr->get<3>()) = diff_phi;
258 }
259 }
260 }
261
262 MoFEMFunctionReturn(0);
263}
NBVOLUMETET_CCG_BUBBLE
#define NBVOLUMETET_CCG_BUBBLE(P)
Bubble function for CGG H div space.
Definition CGGTonsorialBubbleBase.hpp:19
eps
static const double eps
Definition check_base_functions_derivatives_on_tet.cpp:11
USER_BASE
@ USER_BASE
user implemented approximation base
Definition definitions.h:68
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
order
constexpr int order
Definition dg_projection.cpp:19
debug
static const bool debug
Definition dm_create_subdm.cpp:12
diffN_MBTET0z
#define diffN_MBTET0z
derivative of tetrahedral shape function
Definition fem_tools.h:34
N_MBTET0
#define N_MBTET0(x, y, z)
tetrahedral shape function
Definition fem_tools.h:28
ShapeDiffMBTET
PetscErrorCode ShapeDiffMBTET(double *diffN)
calculate derivatives of shape functions
Definition fem_tools.c:319
N_MBTET1
#define N_MBTET1(x, y, z)
tetrahedral shape function
Definition fem_tools.h:29
diffN_MBTET2z
#define diffN_MBTET2z
derivative of tetrahedral shape function
Definition fem_tools.h:40
diffN_MBTET3y
#define diffN_MBTET3y
derivative of tetrahedral shape function
Definition fem_tools.h:42
diffN_MBTET0y
#define diffN_MBTET0y
derivative of tetrahedral shape function
Definition fem_tools.h:33
diffN_MBTET1z
#define diffN_MBTET1z
derivative of tetrahedral shape function
Definition fem_tools.h:37
diffN_MBTET3z
#define diffN_MBTET3z
derivative of tetrahedral shape function
Definition fem_tools.h:43
ShapeMBTET
PetscErrorCode ShapeMBTET(double *N, const double *G_X, const double *G_Y, const double *G_Z, int DIM)
calculate shape functions
Definition fem_tools.c:306
diffN_MBTET1y
#define diffN_MBTET1y
derivative of tetrahedral shape function
Definition fem_tools.h:36
diffN_MBTET2y
#define diffN_MBTET2y
derivative of tetrahedral shape function
Definition fem_tools.h:39
diffN_MBTET2x
#define diffN_MBTET2x
derivative of tetrahedral shape function
Definition fem_tools.h:38
N_MBTET2
#define N_MBTET2(x, y, z)
tetrahedral shape function
Definition fem_tools.h:30
diffN_MBTET3x
#define diffN_MBTET3x
derivative of tetrahedral shape function
Definition fem_tools.h:41
diffN_MBTET1x
#define diffN_MBTET1x
derivative of tetrahedral shape function
Definition fem_tools.h:35
diffN_MBTET0x
#define diffN_MBTET0x
derivative of tetrahedral shape function
Definition fem_tools.h:32
N_MBTET3
#define N_MBTET3(x, y, z)
tetrahedral shape function
Definition fem_tools.h:31
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition LogManager.hpp:316
i
FTensor::Index< 'i', SPACE_DIM > i
Definition hcurl_divergence_operator_2d.cpp:27
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:19
EshelbianPlasticity::CGG_BubbleBase_MBTET
MoFEMErrorCode CGG_BubbleBase_MBTET(const int p, const double *N, const double *diffN, FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > &phi, const int gdim)
Calculate CGGT tonsorial bubble base.
Definition CGGTonsorialBubbleBase.cpp:216
FTensor::dd
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition ddTensor0.hpp:33
MoFEM::Types::MatrixComplexDouble
UBlasMatrix< std::complex< double > > MatrixComplexDouble
Definition Types.hpp:78
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition Types.hpp:77
phi
static double phi
Definition poisson_2d_dis_galerkin.cpp:30

◆ query_interface()

MoFEMErrorCode CGGUserPolynomialBase::query_interface ( boost::typeindex::type_index  type_index,
BaseFunctionUnknownInterface **  iface 
) const

Definition at line 11 of file CGGUserPolynomialBase.cpp.

13 {
14 *iface = const_cast<CGGUserPolynomialBase *>(this);
15 return 0;
16}
CGGUserPolynomialBase
CGG User Polynomial Base.
Definition CGGUserPolynomialBase.hpp:6

Member Data Documentation

◆ cachePhiPtr

boost::shared_ptr<CachePhi> CGGUserPolynomialBase::cachePhiPtr
private

Definition at line 20 of file CGGUserPolynomialBase.hpp.

◆ shapeFun

MatrixDouble CGGUserPolynomialBase::shapeFun
private

Definition at line 19 of file CGGUserPolynomialBase.hpp.

The documentation for this struct was generated from the following files:
Generated by Doxygen 1.9.8 and hosted at University of Glasgow