11 using namespace MoFEM;
14 int *faces_nodes,
int *p,
double *
N,
double *diffN,
double *phi_f_e[4][3],
15 double *diff_phi_f_e[4][3],
int gdim,
16 PetscErrorCode (*base_polynomials)(
int p,
double s,
double *diff_s,
17 double *
L,
double *diffL,
21 for (
int ff = 0; ff < 4; ff++) {
22 if (diff_phi_f_e != NULL) {
24 &faces_nodes[3 * ff], p[ff],
N, diffN, phi_f_e[ff], diff_phi_f_e[ff],
25 gdim, 4, base_polynomials);
29 &faces_nodes[3 * ff], p[ff],
N, diffN, phi_f_e[ff], NULL, gdim, 4,
38 int *faces_nodes,
int p,
double *
N,
double *diffN,
double *phi_f_e[3],
39 double *diff_phi_f_e[3],
int gdim,
int nb,
40 PetscErrorCode (*base_polynomials)(
int p,
double s,
double *diff_s,
41 double *
L,
double *diffL,
44 const int face_edges_nodes[3][2] = {{0, 1}, {1, 2}, {2, 0}};
45 const int face_opposite_edges_node[] = {2, 0, 1};
65 for (
int ee = 0; ee < 3; ee++) {
66 const int n0 = faces_nodes[face_edges_nodes[ee][0]];
67 const int n1 = faces_nodes[face_edges_nodes[ee][1]];
68 t_diff_ksi0i[ee](
i) = t_node_diff_ksi[n1](
i) - t_node_diff_ksi[n0](
i);
69 t_edge_cross[ee](0) = t_node_diff_ksi[n0](1) * t_node_diff_ksi[n1](2) -
70 t_node_diff_ksi[n0](2) * t_node_diff_ksi[n1](1);
71 t_edge_cross[ee](1) = t_node_diff_ksi[n0](2) * t_node_diff_ksi[n1](0) -
72 t_node_diff_ksi[n0](0) * t_node_diff_ksi[n1](2);
73 t_edge_cross[ee](2) = t_node_diff_ksi[n0](0) * t_node_diff_ksi[n1](1) -
74 t_node_diff_ksi[n0](1) * t_node_diff_ksi[n1](0);
77 for (
int ee = 0; ee < 3; ee++) {
78 t_edge_cross[ee](0) = 1;
79 t_edge_cross[ee](1) = 0;
80 t_edge_cross[ee](2) = 0;
83 double psi_l[p + 1], diff_psi_l[3 * (p + 1)];
84 boost::shared_ptr<FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3>>
87 for (
int ee = 0; ee != 3; ee++) {
88 const int i0 = faces_nodes[face_edges_nodes[ee][0]];
89 const int i1 = faces_nodes[face_edges_nodes[ee][1]];
90 const int iO = faces_nodes[face_opposite_edges_node[ee]];
94 t_diff_phi_f_e_ptr = boost::shared_ptr<
103 for (
int ii = 0; ii != gdim; ii++) {
104 const int node_shift = ii * nb;
105 const double n0 =
N[node_shift + i0];
106 const double n1 =
N[node_shift + i1];
107 const double lambda =
N[node_shift + iO];
108 const double ksi0i = n1 - n0;
110 ierr = base_polynomials(p, ksi0i, &t_diff_ksi0i[ee](0), psi_l,
114 ierr = base_polynomials(p, ksi0i, NULL, psi_l, NULL, 3);
119 &diff_psi_l[0], &diff_psi_l[p + 1], &diff_psi_l[2 * p + 2]);
120 for (
int l = 0;
l <= p - 1;
l++) {
121 t_psi_f_e(
i) =
lambda * t_psi_l * t_edge_cross[ee](
i);
122 if (t_diff_phi_f_e_ptr) {
123 (*t_diff_phi_f_e_ptr)(
i,
j) =
124 (t_node_diff_ksi[iO](
j) * t_psi_l +
lambda * t_diff_psi_l(
j)) *
127 ++(*t_diff_phi_f_e_ptr);
138 int *faces_nodes,
int *p,
double *
N,
double *diffN,
double *phi_f[],
139 double *diff_phi_f[],
int gdim,
140 PetscErrorCode (*base_polynomials)(
int p,
double s,
double *diff_s,
141 double *
L,
double *diffL,
145 for (
int ff = 0; ff < 4; ff++) {
147 if (diff_phi_f != NULL) {
148 diff = diff_phi_f[ff];
153 &faces_nodes[3 * ff], p[ff],
N, diffN, phi_f[ff], diff, gdim, 4,
161 int *face_nodes,
int p,
double *
N,
double *diffN,
double *phi_f,
162 double *diff_phi_f,
int gdim,
int nb,
163 PetscErrorCode (*base_polynomials)(
int p,
double s,
double *diff_s,
164 double *
L,
double *diffL,
173 const int vert_i = face_nodes[1];
174 const int vert_j = face_nodes[2];
175 const int i0 = face_nodes[0];
190 t_diff_ksi0i(
i) = t_node_diff_ksi[vert_i](
i) - t_node_diff_ksi[i0](
i);
191 t_diff_ksi0j(
i) = t_node_diff_ksi[vert_j](
i) - t_node_diff_ksi[i0](
i);
192 t_cross(0) = t_node_diff_ksi[vert_i](1) * t_node_diff_ksi[vert_j](2) -
193 t_node_diff_ksi[vert_i](2) * t_node_diff_ksi[vert_j](1);
194 t_cross(1) = t_node_diff_ksi[vert_i](2) * t_node_diff_ksi[vert_j](0) -
195 t_node_diff_ksi[vert_i](0) * t_node_diff_ksi[vert_j](2);
196 t_cross(2) = t_node_diff_ksi[vert_i](0) * t_node_diff_ksi[vert_j](1) -
197 t_node_diff_ksi[vert_i](1) * t_node_diff_ksi[vert_j](0);
204 double psi_l[p + 1], diff_psi_l[3 * (p + 1)];
205 double psi_m[p + 1], diff_psi_m[3 * (p + 1)];
211 boost::shared_ptr<FTensor::Tensor2<double *, 3, 3> > t_diff_phi_f_ptr;
213 t_diff_phi_f_ptr = boost::shared_ptr<FTensor::Tensor2<double *, 3, 3> >(
221 for (
int ii = 0; ii < gdim; ii++) {
223 int node_shift = ii * nb;
224 const double ni =
N[node_shift + vert_i];
225 const double nj =
N[node_shift + vert_j];
226 const double n0 =
N[node_shift + i0];
227 const double ksi0i = ni - n0;
228 const double ksi0j = nj - n0;
229 double beta_0ij = n0 * ni * nj;
231 t_diff_beta_0ij(
i) = (ni * nj) * t_node_diff_ksi[i0](
i) +
232 (n0 * nj) * t_node_diff_ksi[vert_i](
i) +
233 (n0 * ni) * t_node_diff_ksi[vert_j](
i);
234 ierr = base_polynomials(p, ksi0i, &t_diff_ksi0i(0), psi_l, diff_psi_l, 3);
236 ierr = base_polynomials(p, ksi0j, &t_diff_ksi0j(0), psi_m, diff_psi_m, 3);
239 ierr = base_polynomials(p, ksi0i, NULL, psi_l, NULL, 3);
241 ierr = base_polynomials(p, ksi0j, NULL, psi_m, NULL, 3);
247 for (; oo <= p - 3; oo++) {
250 &diff_psi_l[2 * p + 2], 1);
251 for (
int l = 0;
l <= oo;
l++) {
255 diff_psi_m[
m], diff_psi_m[p + 1 +
m], diff_psi_m[2 * p + 2 +
m]);
256 t_psi_f(
i) = (beta_0ij * t_psi_l * psi_m[
m]) * t_cross(
i);
259 (*t_diff_phi_f_ptr)(
i,
j) =
260 ((t_psi_l * psi_m[
m]) * t_diff_beta_0ij(
j) +
261 (beta_0ij * psi_m[
m]) * t_diff_psi_l(
j) +
262 (beta_0ij * t_psi_l) * t_diff_psi_m(
j)) *
264 ++(*t_diff_phi_f_ptr);
281 int p,
double *
N,
double *diffN,
double *phi_v_e[6],
282 double *diff_phi_v_e[6],
int gdim,
283 PetscErrorCode (*base_polynomials)(
int p,
double s,
double *diff_s,
284 double *
L,
double *diffL,
287 const int edges_nodes[6][2] = {{0, 1}, {1, 2}, {2, 0},
288 {0, 3}, {1, 3}, {2, 3}};
316 double diff_psi_l[3 * (p + 1)];
318 for (
int ee = 0; ee != 6; ee++) {
320 t_coords[edges_nodes[ee][1]](
i) - t_coords[edges_nodes[ee][0]](
i);
321 t_diff_ksi0i(
i) = t_node_diff_ksi[edges_nodes[ee][1]](
i) -
322 t_node_diff_ksi[edges_nodes[ee][0]](
i);
330 &diff_phi_v_e[ee][
HVEC2_2], 9);
331 for (
int ii = 0; ii != gdim; ii++) {
332 const int node_shift = ii * 4;
333 const double ni =
N[node_shift + edges_nodes[ee][1]];
334 const double n0 =
N[node_shift + edges_nodes[ee][0]];
335 const double beta_e = ni * n0;
336 const double ksi0i = ni - n0;
338 t_diff_beta_e(
i) = ni * t_node_diff_ksi[edges_nodes[ee][0]](
i) +
339 t_node_diff_ksi[edges_nodes[ee][1]](
i) * n0;
341 base_polynomials(p, ksi0i, &t_diff_ksi0i(0), psi_l, diff_psi_l, 3);
344 ierr = base_polynomials(p, ksi0i, NULL, psi_l, NULL, 3);
349 &diff_psi_l[0], &diff_psi_l[p + 1], &diff_psi_l[2 * p + 2], 1);
350 for (
int l = 0;
l <= p - 2;
l++) {
351 t_psi_v_e(
i) = (beta_e * t_psi_l) * t_tou_e(
i);
354 t_diff_phi_v_e(
i,
j) =
355 (t_diff_beta_e(
j) * t_psi_l + beta_e * t_diff_psi_l(
j)) *
369 int p,
double *
N,
double *diffN,
double *phi_v_f[],
double *diff_phi_v_f[],
371 PetscErrorCode (*base_polynomials)(
int p,
double s,
double *diff_s,
372 double *
L,
double *diffL,
375 const int faces_nodes[4][3] = {{0, 1, 3}, {1, 2, 3}, {0, 2, 3}, {0, 1, 2}};
398 for (
int ff = 0; ff != 4; ff++) {
399 const int v0 = faces_nodes[ff][0];
400 const int vi = faces_nodes[ff][1];
401 const int vj = faces_nodes[ff][2];
402 t_tau0i[ff](
i) = t_coords[vi](
i) - t_coords[v0](
i);
403 t_tau0j[ff](
i) = t_coords[vj](
i) - t_coords[v0](
i);
404 t_diff_ksi0i[ff](
i) = t_node_diff_ksi[vi](
i) - t_node_diff_ksi[v0](
i);
405 t_diff_ksi0j[ff](
i) = t_node_diff_ksi[vj](
i) - t_node_diff_ksi[v0](
i);
408 double psi_l[p + 1], psi_m[p + 1];
409 double diff_psi_l[3 * (p + 1)], diff_psi_m[3 * (p + 1)];
410 for (
int ff = 0; ff != 4; ff++) {
411 const int v0 = faces_nodes[ff][0];
412 const int vi = faces_nodes[ff][1];
413 const int vj = faces_nodes[ff][2];
421 &diff_phi_v_f[ff][
HVEC2_2], 9);
422 for (
int ii = 0; ii < gdim; ii++) {
423 const int node_shift = 4 * ii;
424 const double n0 =
N[node_shift + v0];
425 const double ni =
N[node_shift + vi];
426 const double nj =
N[node_shift + vj];
427 const double beta_f = n0 * ni * nj;
429 t_diff_beta_f(
i) = (ni * nj) * t_node_diff_ksi[v0](
i) +
430 (n0 * nj) * t_node_diff_ksi[vi](
i) +
431 (n0 * ni) * t_node_diff_ksi[vj](
i);
432 const double ksi0i = ni - n0;
433 const double ksi0j = nj - n0;
434 ierr = base_polynomials(p, ksi0i, &t_diff_ksi0i[ff](0), psi_l, diff_psi_l,
437 ierr = base_polynomials(p, ksi0j, &t_diff_ksi0j[ff](0), psi_m, diff_psi_m,
442 for (
int oo = 0; oo <= p - 3; oo++) {
445 diff_psi_l, &diff_psi_l[p + 1], &diff_psi_l[2 * p + 2], 1);
446 for (
int l = 0;
l <= oo;
l++) {
451 &diff_psi_m[
m], &diff_psi_m[p + 1 +
m],
452 &diff_psi_m[2 * p + 2 +
m], 1);
453 const double a = beta_f * t_psi_l * psi_m[
m];
454 t_phi_v_f(
i) =
a * t_tau0i[ff](
i);
457 t_phi_v_f(
i) =
a * t_tau0j[ff](
i);
461 t_diff_a(
j) = (t_psi_l * psi_m[
m]) * t_diff_beta_f(
j) +
462 (beta_f * psi_m[
m]) * t_diff_psi_l(
j) +
463 (beta_f * t_psi_l) * t_diff_psi_m(
j);
464 t_diff_phi_v_f(
i,
j) = t_diff_a(
j) * t_tau0i[ff](
i);
466 t_diff_phi_v_f(
i,
j) = t_diff_a(
j) * t_tau0j[ff](
i);
476 "wrong order %d != %d", jj,
485 int p,
double *
N,
double *diffN,
double *phi_v,
double *diff_phi_v,
487 PetscErrorCode (*base_polynomials)(
int p,
double s,
double *diff_s,
488 double *
L,
double *diffL,
511 t_diff_ksi0i(
i) = t_node_diff_ksi[1](
i) - t_node_diff_ksi[0](
i);
512 t_diff_ksi0j(
i) = t_node_diff_ksi[2](
i) - t_node_diff_ksi[0](
i);
513 t_diff_ksi0k(
i) = t_node_diff_ksi[3](
i) - t_node_diff_ksi[0](
i);
516 double diff_psi_l[3 * (p + 1)];
518 double diff_psi_m[3 * (p + 1)];
520 double diff_psi_n[3 * (p + 1)];
529 for (
int ii = 0; ii < gdim; ii++) {
530 const int node_shift = ii * 4;
531 const double n0 =
N[node_shift + 0];
532 const double ni =
N[node_shift + 1];
533 const double nj =
N[node_shift + 2];
534 const double nk =
N[node_shift + 3];
535 const double ksi0i = ni - n0;
536 const double ksi0j = nj - n0;
537 const double ksi0k = nk - n0;
538 const double beta_v = n0 * ni * nj * nk;
539 t_diff_beta_v(
i) = (ni * nj * nk) * t_node_diff_ksi[0](
i) +
540 (n0 * nj * nk) * t_node_diff_ksi[1](
i) +
541 (n0 * ni * nk) * t_node_diff_ksi[2](
i) +
542 (n0 * ni * nj) * t_node_diff_ksi[3](
i);
543 ierr = base_polynomials(p, ksi0i, &t_diff_ksi0i(0), psi_l, diff_psi_l, 3);
545 ierr = base_polynomials(p, ksi0j, &t_diff_ksi0j(0), psi_m, diff_psi_m, 3);
547 ierr = base_polynomials(p, ksi0k, &t_diff_ksi0k(0), psi_n, diff_psi_n, 3);
553 for (
int oo = 0; oo <= p - 4; oo++) {
556 &diff_psi_l[2 * p + 2], 1);
557 for (
int l = 0;
l <= oo;
l++) {
560 diff_psi_m, &diff_psi_m[p + 1], &diff_psi_m[2 * p + 2], 1);
561 for (
int m = 0; (
l +
m) <= oo;
m++) {
565 diff_psi_n[p + 1 +
n],
566 diff_psi_n[2 * p + 2 +
n]);
567 const double a = beta_v * t_psi_l * t_psi_m * psi_n[
n];
580 t_diff_a(
j) = (t_psi_l * t_psi_m * psi_n[
n]) * t_diff_beta_v(
j) +
581 (beta_v * t_psi_m * psi_n[
n]) * t_diff_psi_l(
j) +
582 (beta_v * t_psi_l * psi_n[
n]) * t_diff_psi_m(
j) +
583 (beta_v * t_psi_l * t_psi_m) * t_diff_psi_n(
j);
584 t_diff_phi_v(N0,
j) = t_diff_a(
j);
585 t_diff_phi_v(N1,
j) = 0;
586 t_diff_phi_v(N2,
j) = 0;
588 t_diff_phi_v(N0,
j) = 0;
589 t_diff_phi_v(N1,
j) = t_diff_a(
j);
590 t_diff_phi_v(N2,
j) = 0;
592 t_diff_phi_v(N0,
j) = 0;
593 t_diff_phi_v(N1,
j) = 0;
594 t_diff_phi_v(N2,
j) = t_diff_a(
j);
608 "wrong order %d != %d", jj,
618 double *diffN,
double *phi_f,
619 double *diff_phi_f,
int gdim,
int nb) {
620 const int face_edges_nodes[3][2] = {{0, 1}, {1, 2}, {2, 0}};
621 const int face_opposite_edges_node[] = {2, 0, 1};
629 FTensor::Tensor2<double, 3, 3> t_diff_cross(0.,0.,0., 0.,0.,0., 0.,0.,0.);
633 const int i0 = faces_nodes[0];
634 const int i1 = faces_nodes[1];
635 const int i2 = faces_nodes[2];
636 const int o[] = {faces_nodes[face_opposite_edges_node[0]],
637 faces_nodes[face_opposite_edges_node[1]],
638 faces_nodes[face_opposite_edges_node[2]]};
652 t_diff_cross(
i,
j) = 0;
653 for (
int ee = 0; ee != 3; ee++) {
654 int ei0 = faces_nodes[face_edges_nodes[ee][0]];
655 int ei1 = faces_nodes[face_edges_nodes[ee][1]];
656 t_cross[ee](0) = t_node_diff_ksi[ei0](1) * t_node_diff_ksi[ei1](2) -
657 t_node_diff_ksi[ei0](2) * t_node_diff_ksi[ei1](1);
658 t_cross[ee](1) = t_node_diff_ksi[ei0](2) * t_node_diff_ksi[ei1](0) -
659 t_node_diff_ksi[ei0](0) * t_node_diff_ksi[ei1](2);
660 t_cross[ee](2) = t_node_diff_ksi[ei0](0) * t_node_diff_ksi[ei1](1) -
661 t_node_diff_ksi[ei0](1) * t_node_diff_ksi[ei1](0);
663 diffN[3 * o[ee] + 0], diffN[3 * o[ee] + 1], diffN[3 * o[ee] + 2]);
664 t_diff_cross(
i,
j) += t_cross[ee](
i) * t_diff_o(
j);
669 t_diff_n0_p_n1(
i) = t_node_diff_ksi[i0](
i) + t_node_diff_ksi[i1](
i);
670 t_diff_n0_p_n1_p_n2(
i) = t_diff_n0_p_n1(
i) + t_node_diff_ksi[i2](
i);
672 for (
int ee = 0; ee != 3; ee++) {
681 boost::shared_ptr<FTensor::Tensor2<double *, 3, 3> > t_diff_phi_ptr;
683 t_diff_phi_ptr = boost::shared_ptr<FTensor::Tensor2<double *, 3, 3> >(
691 double fi[p + 1], diff_fi[3 * p + 3];
692 double fj[p + 1], diff_fj[3 * p + 3];
693 double tmp_fj[p + 1], tmp_diff_fj[3 * p + 3];
694 for (
int ii = 0; ii != gdim; ii++) {
695 const int shift = ii * nb;
696 double n0 =
N[shift + i0];
697 double n1 =
N[shift + i1];
698 double n2 =
N[shift + i2];
699 double *diff_n1 = (diff_phi_f) ? &t_node_diff_ksi[i1](0) : NULL;
700 double *diff_n0_p_n1 = (diff_phi_f) ? &t_diff_n0_p_n1(0) : NULL;
702 diff_phi_f ? diff_fi : NULL, 3);
704 for (
int pp = 0; pp <= p; pp++) {
705 double *diff_n2 = (diff_phi_f) ? &t_node_diff_ksi[i2](0) : NULL;
706 double *diff_n0_p_n1_p_n2 = (diff_phi_f) ? &t_diff_n0_p_n1_p_n2(0) : NULL;
708 diff_n0_p_n1_p_n2, tmp_fj,
709 diff_phi_f ? tmp_diff_fj : NULL, 3);
713 diff_fj[0 * (p + 1) + pp] = tmp_diff_fj[0 * (pp + 1) + pp];
714 diff_fj[1 * (p + 1) + pp] = tmp_diff_fj[1 * (pp + 1) + pp];
715 diff_fj[2 * (p + 1) + pp] = tmp_diff_fj[2 * (pp + 1) + pp];
718 double no0 =
N[shift + o[0]];
719 double no1 =
N[shift + o[1]];
720 double no2 =
N[shift + o[2]];
722 base0(
i) = no0 * t_cross[0](
i) + no1 * t_cross[1](
i) + no2 * t_cross[2](
i);
724 for (
int oo = 0; oo < p; oo++) {
727 &diff_fi[2 * p + 2], 1);
728 for (
int ll = 0; ll <= oo; ll++) {
729 const int mm = oo - ll;
731 const double a = t_fi * fj[mm];
733 t_phi(
i) =
a * base0(
i);
736 &diff_fj[0 + mm], &diff_fj[p + 1 + mm],
737 &diff_fj[2 * p + 2 + mm], 1);
738 (*t_diff_phi_ptr)(
i,
j) =
739 a * t_diff_cross(
i,
j) +
740 (t_diff_fi(
j) * fj[mm] + t_fi * t_diff_fj(
j)) * base0(
i);
752 "wrong number of base functions "
753 "jj!=NBFACETRI_DEMKOWICZ_HDIV(p) "
763 int p,
double *
N,
double *diffN,
int p_f[],
double *phi_f[4],
764 double *diff_phi_f[4],
double *phi_v,
double *diff_phi_v,
int gdim) {
766 const int opposite_face_node[4] = {2, 0, 1, 3};
768 const int znf[] = {0, 2, 3};
784 for (
int ff = 0; ff != 4; ++ff) {
785 t_m_node_diff_ksi[ff](
i) = -t_node_diff_ksi[ff](
i);
798 for (
int ii = 0; ii != gdim; ii++) {
799 const int shift = 4 * ii;
801 for (
int ff = 0; ff != 3; ff++) {
802 const int fff = znf[ff];
803 const int iO = opposite_face_node[fff];
804 const double nO =
N[shift + iO];
805 for (
int pp = 1; pp <= p; pp++) {
807 pp, 2 * pp + 2, nO, 1 - nO, &t_node_diff_ksi[iO](0),
808 &t_m_node_diff_ksi[iO](0), &fk(ff, 0), &diff_fk(ff, 0), 3);
813 for (
int oo = 2; oo <= p; oo++) {
814 for (
int k = 1;
k != oo;
k++) {
821 int sp[] = {ii * 3 * nb_dofs + 3 * s, ii * 3 * nb_dofs + 3 * s,
822 ii * 3 * nb_dofs + 3 * s};
825 &phi_f[znf[0]][sp[0] +
HVEC1],
826 &phi_f[znf[0]][sp[0] +
HVEC2], 3),
828 &phi_f[znf[1]][sp[1] +
HVEC1],
829 &phi_f[znf[1]][sp[1] +
HVEC2], 3),
831 &phi_f[znf[2]][sp[2] +
HVEC1],
832 &phi_f[znf[2]][sp[2] +
HVEC2], 3)};
833 int sdp[] = {ii * 9 * nb_dofs + 9 * s, ii * 9 * nb_dofs + 9 * s,
834 ii * 9 * nb_dofs + 9 * s};
837 &diff_phi_f[znf[0]][sdp[0] +
HVEC0_0],
838 &diff_phi_f[znf[0]][sdp[0] +
HVEC0_1],
839 &diff_phi_f[znf[0]][sdp[0] +
HVEC0_2],
840 &diff_phi_f[znf[0]][sdp[0] +
HVEC1_0],
841 &diff_phi_f[znf[0]][sdp[0] +
HVEC1_1],
842 &diff_phi_f[znf[0]][sdp[0] +
HVEC1_2],
843 &diff_phi_f[znf[0]][sdp[0] +
HVEC2_0],
844 &diff_phi_f[znf[0]][sdp[0] +
HVEC2_1],
845 &diff_phi_f[znf[0]][sdp[0] +
HVEC2_2], 9),
847 &diff_phi_f[znf[1]][sdp[1] +
HVEC0_0],
848 &diff_phi_f[znf[1]][sdp[1] +
HVEC0_1],
849 &diff_phi_f[znf[1]][sdp[1] +
HVEC0_2],
850 &diff_phi_f[znf[1]][sdp[1] +
HVEC1_0],
851 &diff_phi_f[znf[1]][sdp[1] +
HVEC1_1],
852 &diff_phi_f[znf[1]][sdp[1] +
HVEC1_2],
853 &diff_phi_f[znf[1]][sdp[1] +
HVEC2_0],
854 &diff_phi_f[znf[1]][sdp[1] +
HVEC2_1],
855 &diff_phi_f[znf[1]][sdp[1] +
HVEC2_2], 9),
857 &diff_phi_f[znf[2]][sdp[2] +
HVEC0_0],
858 &diff_phi_f[znf[2]][sdp[2] +
HVEC0_1],
859 &diff_phi_f[znf[2]][sdp[2] +
HVEC0_2],
860 &diff_phi_f[znf[2]][sdp[2] +
HVEC1_0],
861 &diff_phi_f[znf[2]][sdp[2] +
HVEC1_1],
862 &diff_phi_f[znf[2]][sdp[2] +
HVEC1_2],
863 &diff_phi_f[znf[2]][sdp[2] +
HVEC2_0],
864 &diff_phi_f[znf[2]][sdp[2] +
HVEC2_1],
865 &diff_phi_f[znf[2]][sdp[2] +
HVEC2_2], 9)};
867 for (
int ff = 0; ff != 3; ff++) {
869 diff_fk(ff, 1 * p +
k - 1),
870 diff_fk(ff, 2 * p +
k - 1));
871 t_phi_v(
i) = fk(ff,
k - 1) * t_phi_f[ff](
i);
872 t_diff_phi_v(
i,
j) = t_diff_fk(
j) * t_phi_f[ff](
i) +
873 fk(ff,
k - 1) * t_diff_phi_f[ff](
i,
j);
886 "wrong number of base functions "
887 "jj!=NBVOLUMETET_DEMKOWICZ_HDIV(p) "