Calculate CGGT tonsorial bubble base.
23 {
25
26
27
30
36
38 {
40 &diffN[2]);
41 for (int ii = 0; ii != 4; ++ii) {
42 t_diff_n[ii](
i) = t_diff_n_tmp(
i);
43 ++t_diff_n_tmp;
44 }
45 }
46
48 for (int ii = 0; ii != 3; ++ii)
49 t_diff_ksi[ii](
i) = t_diff_n[ii + 1](
i) - t_diff_n[0](
i);
50
51 int lp =
p >= 2 ?
p - 2 + 1 : 0;
58
59 for (int ii = 0; ii != 3;++ii)
60 diff2_l[ii].clear();
61
62 for (int gg = 0; gg != gdim; ++gg) {
63
64 const int node_shift = gg * 4;
65
66 for (int ii = 0; ii != 3; ++ii) {
67
68 auto &t_diff_ksi_ii = t_diff_ksi[ii];
70 auto &diff_l_ii = diff_l[ii];
71 auto &diff2_l_ii = diff2_l[ii];
72
73 double ksi_ii =
N[node_shift + ii + 1] -
N[node_shift + 0];
74
76 &*l_ii.data().begin(),
77 &*diff_l_ii.data().begin(), 3);
78
79 for (
int l = 1;
l < lp; ++
l) {
81 const double b = ((
double)
l / ((
double)
l + 1));
82 for (int d0 = 0; d0 != 3; ++d0)
83 for (int d1 = 0; d1 != 3; ++d1) {
84 const int r = 3 * d0 + d1;
85 diff2_l_ii(r,
l + 1) =
a * (t_diff_ksi_ii(d0) * diff_l_ii(d1,
l) +
86 t_diff_ksi_ii(d1) * diff_l_ii(d0,
l) +
87 ksi_ii * diff2_l_ii(r,
l)) -
88 b * diff2_l_ii(r,
l - 1);
89 }
90 }
91
92 }
93
94 const double n[] = {
N[node_shift + 0],
N[node_shift + 1],
N[node_shift + 2],
96
99 const int tab[4][4] = {
100 {1, 2, 3, 0}, {2, 3, 0, 1}, {3, 0, 1, 2}, {0, 1, 2, 3}};
102 t_bk_diff(
i,
j,
k) = 0;
103 for (int ii = 0; ii != 3; ++ii) {
104 const int i0 = tab[ii][0];
105 const int i1 = tab[ii][1];
106 const int i2 = tab[ii][2];
107 const int i3 = tab[ii][3];
108 auto &t_diff_n_i0 = t_diff_n[i0];
109 auto &t_diff_n_i1 = t_diff_n[i1];
110 auto &t_diff_n_i2 = t_diff_n[i2];
111 auto &t_diff_n_i3 = t_diff_n[i3];
113 t_k(
i,
j) = t_diff_n_i3(
i) * t_diff_n_i3(
j);
114 const double b =
n[i0] *
n[i1] *
n[i2];
115 t_bk(
i,
j) +=
b * t_k(
i,
j);
117 t_diff_b(
i) = t_diff_n_i0(
i) *
n[i1] *
n[i2] +
118 t_diff_n_i1(
i) *
n[i0] *
n[i2] +
119 t_diff_n_i2(
i) *
n[i0] *
n[i1];
120 t_bk_diff(
i,
j,
k) += t_k(
i,
j) * t_diff_b(
k);
121 }
122
123 int zz = 0;
124 for (
int o =
p - 2 + 1;
o <=
p - 2 + 1; ++
o) {
125
126 for (
int ii = 0; ii <=
o; ++ii)
127 for (
int jj = 0; (ii + jj) <=
o; ++jj) {
128
129 const int kk =
o - ii - jj;
130
131 auto get_diff_l = [&](
const int y,
const int i) {
134 };
135 auto get_diff2_l = [&](
const int y,
const int i) {
137 diff2_l[y](0,
i), diff2_l[y](1,
i), diff2_l[y](2,
i),
138 diff2_l[y](3,
i), diff2_l[y](4,
i), diff2_l[y](5,
i),
139 diff2_l[y](6,
i), diff2_l[y](7,
i), diff2_l[y](8,
i));
140 };
141
143 auto t_diff_i = get_diff_l(0, ii);
144 auto t_diff2_i = get_diff2_l(0, ii);
146 auto t_diff_j = get_diff_l(1, jj);
147 auto t_diff2_j = get_diff2_l(1, jj);
149 auto t_diff_k = get_diff_l(2, kk);
150 auto t_diff2_k = get_diff2_l(2, kk);
151
153 t_diff_l2(
i) = t_diff_i(
i) * l_j * l_k + t_diff_j(
i) * l_i * l_k +
154 t_diff_k(
i) * l_i * l_j;
157 t_diff2_i(
i,
j) * l_j * l_k + t_diff_i(
i) * t_diff_j(
j) * l_k +
158 t_diff_i(
i) * l_j * t_diff_k(
j) +
159
160 t_diff2_j(
i,
j) * l_i * l_k + t_diff_j(
i) * t_diff_i(
j) * l_k +
161 t_diff_j(
i) * l_i * t_diff_k(
j) +
162
163 t_diff2_k(
i,
j) * l_i * l_j + t_diff_k(
i) * t_diff_i(
j) * l_j +
164 t_diff_k(
i) * l_i * t_diff_j(
j);
165
166 for (
int dd = 0;
dd != 3; ++
dd) {
167
169 t_axial_diff(
i,
j) = 0;
170 for (int mm = 0; mm != 3; ++mm)
171 t_axial_diff(
dd, mm) = t_diff_l2(mm);
172
178 t_curl_A_bK_diff(
i,
j,
k) = t_curl_A(
i,
m) * t_bk_diff(
m,
j,
k);
179
181 t_axial_diff2(
i,
j,
k) = 0;
182 for (int mm = 0; mm != 3; ++mm)
183 for (int nn = 0; nn != 3; ++nn)
184 t_axial_diff2(
dd, mm, nn) = t_diff2_l2(mm, nn);
186 t_A_diff2(
i,
j,
k, f) =
189 t_curl_A_diff2(
i,
j,
k) =
192 t_curl_A_diff2_bK(
i,
j,
k) = t_curl_A_diff2(
i,
m,
k) * t_bk(
m,
j);
193
195 t_curl_A_diff2_bK(
i, f,
m));
196
197 ++t_phi;
198 ++zz;
199 }
200 }
201 }
204 "Wrong number of base functions %d != %d", zz,
206
207 }
208
210}
#define NBVOLUMETET_CCG_BUBBLE(P)
Bubble function for CGG H div space.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ MOFEM_DATA_INCONSISTENCY
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
PetscErrorCode Legendre_polynomials(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Calculate Legendre approximation basis.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
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)
UBlasMatrix< double > MatrixDouble
UBlasVector< double > VectorDouble
