13 using namespace MoFEM;
41 for (
int ii = 0; ii != 4; ++ii) {
42 t_diff_n[ii](
i) = t_diff_n_tmp(
i);
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);
51 int lp = p >= 2 ? p - 2 + 1 : 0;
59 for (
int ii = 0; ii != 3; ++ii)
62 for (
int gg = 0; gg != gdim; ++gg) {
64 const int node_shift = gg * 4;
66 for (
int ii = 0; ii != 3; ++ii) {
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];
73 double ksi_ii =
N[node_shift + ii + 1] -
N[node_shift + 0];
76 &*l_ii.data().begin(),
77 &*diff_l_ii.data().begin(), 3);
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);
93 const double n[] = {
N[node_shift + 0],
N[node_shift + 1],
N[node_shift + 2],
98 const int tab[4][4] = {
99 {1, 2, 3, 0}, {2, 3, 0, 1}, {3, 0, 1, 2}, {0, 1, 2, 3}};
101 t_bk_diff(
i,
j,
k) = 0;
102 for (
int ii = 0; ii != 3; ++ii) {
103 const int i0 = tab[ii][0];
104 const int i1 = tab[ii][1];
105 const int i2 = tab[ii][2];
106 const int i3 = tab[ii][3];
107 auto &t_diff_n_i0 = t_diff_n[i0];
108 auto &t_diff_n_i1 = t_diff_n[i1];
109 auto &t_diff_n_i2 = t_diff_n[i2];
110 auto &t_diff_n_i3 = t_diff_n[i3];
112 t_k(
i,
j) = t_diff_n_i3(
i) * t_diff_n_i3(
j);
113 const double b =
n[i0] *
n[i1] *
n[i2];
114 t_bk(
i,
j) += b * t_k(
i,
j);
116 t_diff_b(
i) = t_diff_n_i0(
i) *
n[i1] *
n[i2] +
117 t_diff_n_i1(
i) *
n[i0] *
n[i2] +
118 t_diff_n_i2(
i) *
n[i0] *
n[i1];
119 t_bk_diff(
i,
j,
k) += t_k(
i,
j) * t_diff_b(
k);
123 for (
int o = p - 2 + 1; o <= p - 2 + 1; ++o) {
125 for (
int ii = 0; ii <= o; ++ii)
126 for (
int jj = 0; (ii + jj) <= o; ++jj) {
128 const int kk = o - ii - jj;
130 auto get_diff_l = [&](
const int y,
const int i) {
134 auto get_diff2_l = [&](
const int y,
const int i) {
136 diff2_l[y](0,
i), diff2_l[y](1,
i), diff2_l[y](2,
i),
137 diff2_l[y](3,
i), diff2_l[y](4,
i), diff2_l[y](5,
i),
138 diff2_l[y](6,
i), diff2_l[y](7,
i), diff2_l[y](8,
i));
142 auto t_diff_i = get_diff_l(0, ii);
143 auto t_diff2_i = get_diff2_l(0, ii);
145 auto t_diff_j = get_diff_l(1, jj);
146 auto t_diff2_j = get_diff2_l(1, jj);
148 auto t_diff_k = get_diff_l(2, kk);
149 auto t_diff2_k = get_diff2_l(2, kk);
152 t_diff_l2(
i) = t_diff_i(
i) * l_j * l_k + t_diff_j(
i) * l_i * l_k +
153 t_diff_k(
i) * l_i * l_j;
156 t_diff2_i(
i,
j) * l_j * l_k + t_diff_i(
i) * t_diff_j(
j) * l_k +
157 t_diff_i(
i) * l_j * t_diff_k(
j) +
159 t_diff2_j(
i,
j) * l_i * l_k + t_diff_j(
i) * t_diff_i(
j) * l_k +
160 t_diff_j(
i) * l_i * t_diff_k(
j) +
162 t_diff2_k(
i,
j) * l_i * l_j + t_diff_k(
i) * t_diff_i(
j) * l_j +
163 t_diff_k(
i) * l_i * t_diff_j(
j);
165 for (
int dd = 0;
dd != 3; ++
dd) {
168 t_axial_diff(
i,
j) = 0;
169 for (
int mm = 0; mm != 3; ++mm)
170 t_axial_diff(
dd, mm) = t_diff_l2(mm);
177 t_curl_A_bK_diff(
i,
j,
k) = t_curl_A(
i,
m) * t_bk_diff(
m,
j,
k);
180 t_axial_diff2(
i,
j,
k) = 0;
181 for (
int mm = 0; mm != 3; ++mm)
182 for (
int nn = 0; nn != 3; ++nn)
183 t_axial_diff2(
dd, mm, nn) = t_diff2_l2(mm, nn);
185 t_A_diff2(
i,
j,
k,
f) =
188 t_curl_A_diff2(
i,
j,
k) =
191 t_curl_A_diff2_bK(
i,
j,
k) = t_curl_A_diff2(
i,
m,
k) * t_bk(
m,
j);
194 t_curl_A_diff2_bK(
i,
f,
m));
203 "Wrong number of base functions %d != %d", zz,