v0.15.5
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, bool skip_calculation=false)
 
 ~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
 
bool skipCalulation
 

Detailed Description

CGG User Polynomial Base.

Examples
/home/lk58p/mofem_install/vanilla_dev_release/mofem-cephas/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,
bool  skip_calculation = false 
)

Definition at line 7 of file CGGUserPolynomialBase.cpp.

9 : TetPolynomialBase(), cachePhiPtr(cache_phi_otr),
10 skipCalulation(skip_calculation) {}
CGGUserPolynomialBase::cachePhiPtr
boost::shared_ptr< CachePhi > cachePhiPtr
Definition CGGUserPolynomialBase.hpp:21

◆ ~CGGUserPolynomialBase()

CGGUserPolynomialBase::~CGGUserPolynomialBase ( )
default

Member Function Documentation

◆ getValue()

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

Definition at line 20 of file CGGUserPolynomialBase.cpp.

21 {
22 MoFEMFunctionBeginHot;
23
24 this->cTx = ctx_ptr->getInterface<EntPolynomialBaseCtx>();
25
26 int nb_gauss_pts = pts.size2();
27 if (!nb_gauss_pts) {
28 MoFEMFunctionReturnHot(0);
29 }
30
31 if (pts.size1() < 3) {
32 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
33 "Wrong dimension of pts, should be at least 3 rows with "
34 "coordinates");
35 }
36
37 const auto base = this->cTx->bAse;
38 EntitiesFieldData &data = this->cTx->dAta;
39
40 switch (this->cTx->sPace) {
41 case HDIV:
42 CHKERR getValueHdivForCGGBubble(pts);
43 break;
44 case L2:
45 data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 4, false);
46 CHKERR Tools::shapeFunMBTET(
47 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
48 &pts(0, 0), &pts(1, 0), &pts(2, 0), nb_gauss_pts);
49 data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(4, 3, false);
50 std::copy(Tools::diffShapeFunMBTET.begin(), Tools::diffShapeFunMBTET.end(),
51 data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin());
52 this->cTx->basePolynomialsType0 = Legendre_polynomials;
53 CHKERR getValueL2AinsworthBase(pts);
54 break;
55 default:
56 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
57 }
58
59 MoFEMFunctionReturnHot(0);
60}
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:63

◆ getValueHdivForCGGBubble()

MoFEMErrorCode CGGUserPolynomialBase::getValueHdivForCGGBubble ( MatrixDouble &  pts)
private

number of integration points

Definition at line 63 of file CGGUserPolynomialBase.cpp.

63 {
64 MoFEMFunctionBegin;
65
66 constexpr bool calculate_diff = false;
67 constexpr bool debug = false;
68
69 // This should be used only in case USER_BASE is selected
70 if (this->cTx->bAse != USER_BASE) {
71 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
72 "Wrong base, should be USER_BASE");
73 }
74 // get access to data structures on element
75 EntitiesFieldData &data = this->cTx->dAta;
76 // Get approximation order on element. Note that bubble functions are only
77 // on tetrahedron.
78 const int order = data.dataOnEntities[MBTET][0].getOrder();
79 /// number of integration points
80 const int nb_gauss_pts = pts.size2();
81
82 auto check_cache = [this](int order, int nb_gauss_pts) -> bool {
83 if (cachePhiPtr) {
84 return cachePhiPtr->get<0>() == order &&
85 cachePhiPtr->get<1>() == nb_gauss_pts;
86 } else {
87 return false;
88 }
89 };
90
91 if (skipCalulation) {
92 auto &phi = data.dataOnEntities[MBTET][0].getN(USER_BASE);
93 const int nb_base_functions = NBVOLUMETET_CCG_BUBBLE(order);
94 phi.resize(nb_gauss_pts, 9 * nb_base_functions, false);
95 std::fill(phi.data().begin(), phi.data().end(), 0.0);
96 MoFEMFunctionReturnHot(0);
97 }
98
99 if (check_cache(order, nb_gauss_pts)) {
100 auto &mat = cachePhiPtr->get<2>();
101 auto &phi = data.dataOnEntities[MBTET][0].getN(USER_BASE);
102 phi.resize(mat.size1(), mat.size2(), false);
103 noalias(phi) = mat;
104 if constexpr (calculate_diff) {
105 auto &diff_phi = data.dataOnEntities[MBTET][0].getDiffN(USER_BASE);
106 auto &diff_mat = cachePhiPtr->get<3>();
107 diff_phi.resize(diff_mat.size1(), diff_mat.size2(), false);
108 noalias(diff_phi) = diff_mat;
109 }
110 } else {
111 // calculate shape functions, i.e. barycentric coordinates
112 shapeFun.resize(nb_gauss_pts, 4, false);
113 CHKERR ShapeMBTET(&*shapeFun.data().begin(), &pts(0, 0), &pts(1, 0),
114 &pts(2, 0), nb_gauss_pts);
115 // derivatives of shape functions
116 double diff_shape_fun[12];
117 CHKERR ShapeDiffMBTET(diff_shape_fun);
118
119 const int nb_base_functions = NBVOLUMETET_CCG_BUBBLE(order);
120 // get base functions and set size
121 MatrixDouble &phi = data.dataOnEntities[MBTET][0].getN(USER_BASE);
122 phi.resize(nb_gauss_pts, 9 * nb_base_functions, false);
123 // finally calculate base functions
124 FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_phi(
125 &phi(0, 0), &phi(0, 1), &phi(0, 2),
126
127 &phi(0, 3), &phi(0, 4), &phi(0, 5),
128
129 &phi(0, 6), &phi(0, 7), &phi(0, 8));
130 CHKERR CGG_BubbleBase_MBTET(order, &shapeFun(0, 0), diff_shape_fun, t_phi,
131 nb_gauss_pts);
132
133 auto imagingray_fd_directive = [&](auto &diff_phi) {
134 MoFEMFunctionBegin;
135
136 MatrixComplexDouble temp_phi;
137 temp_phi.resize(nb_gauss_pts, 9 * nb_base_functions, false);
138
139 MatrixComplexDouble diff_shape_fun;
140 diff_shape_fun.resize(4, 3, false);
141 diff_shape_fun(0, 0) = diffN_MBTET0x;
142 diff_shape_fun(0, 1) = diffN_MBTET0y;
143 diff_shape_fun(0, 2) = diffN_MBTET0z;
144 diff_shape_fun(1, 0) = diffN_MBTET1x;
145 diff_shape_fun(1, 1) = diffN_MBTET1y;
146 diff_shape_fun(1, 2) = diffN_MBTET1z;
147 diff_shape_fun(2, 0) = diffN_MBTET2x;
148 diff_shape_fun(2, 1) = diffN_MBTET2y;
149 diff_shape_fun(2, 2) = diffN_MBTET2z;
150 diff_shape_fun(3, 0) = diffN_MBTET3x;
151 diff_shape_fun(3, 1) = diffN_MBTET3y;
152 diff_shape_fun(3, 2) = diffN_MBTET3z;
153
154 MatrixComplexDouble shape_fun;
155 shape_fun.resize(nb_gauss_pts, 4, false);
156 for (int dd = 0; dd < 3; ++dd) {
157 constexpr double eps = 1e-16;
158
159 for (auto gg = 0; gg < nb_gauss_pts; ++gg) {
160 std::array<std::complex<double>, 3> ksi{pts(0, gg), pts(1, gg),
161 pts(2, gg)};
162
163 ksi[dd] += std::complex<double>(0, eps);
164
165 shape_fun(gg, 0) = N_MBTET0(ksi[0], ksi[1], ksi[2]);
166 shape_fun(gg, 1) = N_MBTET1(ksi[0], ksi[1], ksi[2]);
167 shape_fun(gg, 2) = N_MBTET2(ksi[0], ksi[1], ksi[2]);
168 shape_fun(gg, 3) = N_MBTET3(ksi[0], ksi[1], ksi[2]);
169 }
170
171 auto t_complex_phi_in =
172 FTensor::Tensor2<FTensor::PackPtr<std::complex<double> *, 9>, 3, 3>{
173 &temp_phi(0, 0), &temp_phi(0, 1), &temp_phi(0, 2),
174 &temp_phi(0, 3), &temp_phi(0, 4), &temp_phi(0, 5),
175 &temp_phi(0, 6), &temp_phi(0, 7), &temp_phi(0, 8)};
176
177 CHKERR CGG_BubbleBase_MBTET(order, shape_fun.data().data(),
178 diff_shape_fun.data().data(),
179 t_complex_phi_in, nb_gauss_pts);
180
181 auto t_complex_phi_out =
182 FTensor::Tensor2<FTensor::PackPtr<std::complex<double> *, 9>, 3, 3>{
183 &temp_phi(0, 0), &temp_phi(0, 1), &temp_phi(0, 2),
184 &temp_phi(0, 3), &temp_phi(0, 4), &temp_phi(0, 5),
185 &temp_phi(0, 6), &temp_phi(0, 7), &temp_phi(0, 8)};
186 auto t_diff_phi = getFTensor3FromPtr<3, 3, 3>(diff_phi.data().data());
187 for (auto gg = 0; gg < nb_gauss_pts; ++gg) {
188 for (auto bb = 0; bb << nb_base_functions; ++bb) {
189 t_diff_phi(0, 0, dd) = std::imag(t_complex_phi_out(0, 0)) / eps;
190 t_diff_phi(0, 1, dd) = std::imag(t_complex_phi_out(0, 1)) / eps;
191 t_diff_phi(0, 2, dd) = std::imag(t_complex_phi_out(0, 2)) / eps;
192 t_diff_phi(1, 0, dd) = std::imag(t_complex_phi_out(1, 0)) / eps;
193 t_diff_phi(1, 1, dd) = std::imag(t_complex_phi_out(1, 1)) / eps;
194 t_diff_phi(1, 2, dd) = std::imag(t_complex_phi_out(1, 2)) / eps;
195 t_diff_phi(2, 0, dd) = std::imag(t_complex_phi_out(2, 0)) / eps;
196 t_diff_phi(2, 1, dd) = std::imag(t_complex_phi_out(2, 1)) / eps;
197 t_diff_phi(2, 2, dd) = std::imag(t_complex_phi_out(2, 2)) / eps;
198 if (debug) {
199 MOFEM_LOG("SELF", Sev::noisy)
200 << "Gauss pt " << gg << " direction " << dd
201 << " diff_phi: " << t_diff_phi(0, 0, dd) << ", "
202 << t_diff_phi(0, 1, dd) << ", " << t_diff_phi(0, 2, dd)
203 << "; " << t_diff_phi(1, 0, dd) << ", "
204 << t_diff_phi(1, 1, dd) << ", " << t_diff_phi(1, 2, dd)
205 << "; " << t_diff_phi(2, 0, dd) << ", "
206 << t_diff_phi(2, 1, dd) << ", " << t_diff_phi(2, 2, dd);
207 }
208 ++t_complex_phi_out;
209 ++t_diff_phi;
210 }
211 }
212
213 if (debug) {
214 auto t_phi = FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3>{
215 &phi(0, 0), &phi(0, 1), &phi(0, 2), &phi(0, 3), &phi(0, 4),
216 &phi(0, 5), &phi(0, 6), &phi(0, 7), &phi(0, 8)};
217 auto t_complex_phi_out =
218 FTensor::Tensor2<FTensor::PackPtr<std::complex<double> *, 9>, 3,
219 3>{
220 &temp_phi(0, 0), &temp_phi(0, 1), &temp_phi(0, 2),
221 &temp_phi(0, 3), &temp_phi(0, 4), &temp_phi(0, 5),
222 &temp_phi(0, 6), &temp_phi(0, 7), &temp_phi(0, 8)};
223
224 for (auto gg = 0; gg < nb_gauss_pts; ++gg) {
225 FTensor::Index<'i', 3> i;
226 FTensor::Index<'j', 3> j;
227 FTensor::Tensor2<std::complex<double>, 3, 3> t_diff;
228 for (auto bb = 0; bb << nb_base_functions; ++bb) {
229 t_diff(i, j) = t_complex_phi_out(i, j) - t_phi(i, j);
230 auto nrm2_complex = t_diff(i, j) * t_diff(i, j);
231 if (std::abs(nrm2_complex) > 1e-12) {
232 MOFEM_LOG("SELF", Sev::warning)
233 << "Difference detected at gauss pt " << gg
234 << " for direction " << dd
235 << " nrm2 = " << std::abs(nrm2_complex);
236 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
237 "Difference detected between complex step "
238 "differentiation and real valued differentiation");
239 }
240 ++t_phi;
241 ++t_complex_phi_out;
242 }
243 }
244 }
245
246
247
248 }
249
250 MoFEMFunctionReturn(0);
251 };
252
253 MatrixDouble &diff_phi = data.dataOnEntities[MBTET][0].getDiffN(USER_BASE);
254 if constexpr (calculate_diff) {
255 diff_phi.resize(nb_gauss_pts, 27 * nb_base_functions, false);
256 CHKERR imagingray_fd_directive(diff_phi);
257 }
258
259 if (cachePhiPtr) {
260 cachePhiPtr->get<0>() = order;
261 cachePhiPtr->get<1>() = nb_gauss_pts;
262 cachePhiPtr->get<2>().resize(phi.size1(), phi.size2(), false);
263 noalias(cachePhiPtr->get<2>()) = phi;
264 if constexpr (calculate_diff) {
265 cachePhiPtr->get<3>().resize(diff_phi.size1(), diff_phi.size2(), false);
266 noalias(cachePhiPtr->get<3>()) = diff_phi;
267 }
268 }
269 }
270
271 MoFEMFunctionReturn(0);
272}
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:38
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 12 of file CGGUserPolynomialBase.cpp.

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

Member Data Documentation

◆ cachePhiPtr

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

Definition at line 21 of file CGGUserPolynomialBase.hpp.

◆ shapeFun

MatrixDouble CGGUserPolynomialBase::shapeFun
private

Definition at line 20 of file CGGUserPolynomialBase.hpp.

◆ skipCalulation

bool CGGUserPolynomialBase::skipCalulation
private

Definition at line 22 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