14    [](
int p) { 
return p; };
 
   16    [](
int p) { 
return p; };
 
   18    [](
int p) { 
return p; };
 
   20    [](
int p) { 
return p; };
 
   22    [](
int p) { 
return p; };
 
   25    int *faces_nodes, 
int *p, 
double *
N, 
double *diffN, 
double *phi_f_e[4][3],
 
   26    double *diff_phi_f_e[4][3], 
int gdim,
 
   27    PetscErrorCode (*base_polynomials)(
int p, 
double s, 
double *diff_s,
 
   28                                       double *
L, 
double *diffL,
 
   35            "expected to return derivatives");
 
   38  for (
int ff = 0; ff < 4; ff++) {
 
   40        &faces_nodes[3 * ff], p[ff], 
N, diffN, phi_f_e[ff], diff_phi_f_e[ff],
 
   41        gdim, 4, base_polynomials);
 
 
   48    int *faces_nodes, 
int p, 
double *
N, 
double *diffN, 
double *phi_f_e[3],
 
   49    double *diff_phi_f_e[3], 
int gdim, 
int nb,
 
   50    PetscErrorCode (*base_polynomials)(
int p, 
double s, 
double *diff_s,
 
   51                                       double *
L, 
double *diffL,
 
   54  constexpr int face_edges_nodes[3][2] = {{0, 1}, {1, 2}, {2, 0}};
 
   55  constexpr int face_opposite_edges_node[] = {2, 0, 1};
 
   76    for (
int ee = 0; ee < 3; ee++) {
 
   77      const int n0 = faces_nodes[face_edges_nodes[ee][0]];
 
   78      const int n1 = faces_nodes[face_edges_nodes[ee][1]];
 
   79      t_diff_ksi[ee](
i) = t_node_diff_ksi[n1](
i) - t_node_diff_ksi[n0](
i);
 
   80      t_edge_cross[ee](
i) = levi_civita(
i, 
j, 
k) * t_node_diff_ksi[n0](
j) *
 
   81                            t_node_diff_ksi[n1](
k);
 
   84    for (
int ee = 0; ee < 3; ee++) {
 
   85      t_edge_cross[ee](0) = 1;
 
   86      t_edge_cross[ee](1) = 0;
 
   87      t_edge_cross[ee](2) = 0;
 
   92  for (
int ee = 0; ee != 3; ee++) {
 
   93    const int i0 = faces_nodes[face_edges_nodes[ee][0]];
 
   94    const int i1 = faces_nodes[face_edges_nodes[ee][1]];
 
   95    const int iO = faces_nodes[face_opposite_edges_node[ee]];
 
   96    auto t_psi_f_e = getFTensor1FromPtr<3>(phi_f_e[ee]);
 
   97    for (
int ii = 0; ii != nb * gdim; ii += nb) {
 
   98      const double n0 = 
N[ii + i0];
 
   99      const double n1 = 
N[ii + i1];
 
  100      const double lambda = 
N[ii + iO];
 
  101      const double ksi = n1 - n0;
 
  102      base_polynomials(p, ksi, NULL, psi_l, NULL, 3);
 
  105      for (
int l = 0; 
l <= p - 1; 
l++) {
 
  106        t_psi_f_e(
i) = (
lambda * t_psi_l) * t_edge_cross[ee](
i);
 
  121    double diff_psi_l[3 * (p + 1)];
 
  122    for (
int ee = 0; ee != 3; ee++) {
 
  123      const int i0 = faces_nodes[face_edges_nodes[ee][0]];
 
  124      const int i1 = faces_nodes[face_edges_nodes[ee][1]];
 
  125      const int iO = faces_nodes[face_opposite_edges_node[ee]];
 
  127      for (
int ii = 0; ii != nb * gdim; ii += nb) {
 
  128        const double n0 = 
N[ii + i0];
 
  129        const double n1 = 
N[ii + i1];
 
  130        const double lambda = 
N[ii + iO];
 
  131        const double ksi = n1 - n0;
 
  132        base_polynomials(p, ksi, &t_diff_ksi[ee](0), psi_l, diff_psi_l, 3);
 
  135            &diff_psi_l[0], &diff_psi_l[p + 1], &diff_psi_l[2 * p + 2]};
 
  136        for (
int l = 0; 
l <= p - 1; 
l++) {
 
  137          t_diff_phi_f_e(
i, 
j) =
 
  138              (t_node_diff_ksi[iO](
j) * t_psi_l + 
lambda * t_diff_psi_l(
j)) *
 
 
  152    int *faces_nodes, 
int *p, 
double *
N, 
double *diffN, 
double *phi_f[],
 
  153    double *diff_phi_f[], 
int gdim,
 
  154    PetscErrorCode (*base_polynomials)(
int p, 
double s, 
double *diff_s,
 
  155                                       double *
L, 
double *diffL,
 
  163            "expected to return derivatives");
 
  166  for (
int ff = 0; ff < 4; ff++) {
 
  168        &faces_nodes[3 * ff], p[ff], 
N, diffN, phi_f[ff], diff_phi_f[ff], gdim,
 
  169        4, base_polynomials);
 
 
  175    int *face_nodes, 
int p, 
double *
N, 
double *diffN, 
double *phi_f,
 
  176    double *diff_phi_f, 
int gdim, 
int nb,
 
  177    PetscErrorCode (*base_polynomials)(
int p, 
double s, 
double *diff_s,
 
  178                                       double *
L, 
double *diffL,
 
  188  const int vert_i = face_nodes[1];
 
  189  const int vert_j = face_nodes[2];
 
  190  const int i0 = face_nodes[0];
 
  205    t_diff_ksi0i(
i) = t_node_diff_ksi[vert_i](
i) - t_node_diff_ksi[i0](
i);
 
  206    t_diff_ksi0j(
i) = t_node_diff_ksi[vert_j](
i) - t_node_diff_ksi[i0](
i);
 
  207    t_cross(
i) = levi_civita(
i, 
j, 
k) * t_node_diff_ksi[vert_i](
j) *
 
  208                 t_node_diff_ksi[vert_j](
k);
 
  215  double psi_l[p + 1], diff_psi_l[3 * (p + 1)];
 
  216  double psi_m[p + 1], diff_psi_m[3 * (p + 1)];
 
  217  auto t_psi_f = getFTensor1FromPtr<3>(phi_f);
 
  219  for (
int ii = 0; ii < gdim; ii++) {
 
  221    const int node_shift = ii * nb;
 
  222    const double ni = 
N[node_shift + vert_i];
 
  223    const double nj = 
N[node_shift + vert_j];
 
  224    const double n0 = 
N[node_shift + i0];
 
  225    const double ksi0i = ni - n0;
 
  226    const double ksi0j = nj - n0;
 
  227    double beta_0ij = n0 * ni * nj;
 
  228    base_polynomials(p, ksi0i, NULL, psi_l, NULL, 3);
 
  229    base_polynomials(p, ksi0j, NULL, psi_m, NULL, 3);
 
  233    for (; oo <= p - 3; oo++) {
 
  235      for (
int l = 0; 
l <= oo; 
l++) {
 
  238          t_psi_f(
i) = (beta_0ij * t_psi_l * psi_m[
m]) * t_cross(
i);
 
  256    for (
int ii = 0; ii < gdim; ii++) {
 
  258      const int node_shift = ii * nb;
 
  259      const double ni = 
N[node_shift + vert_i];
 
  260      const double nj = 
N[node_shift + vert_j];
 
  261      const double n0 = 
N[node_shift + i0];
 
  262      const double ksi0i = ni - n0;
 
  263      const double ksi0j = nj - n0;
 
  264      double beta_0ij = n0 * ni * nj;
 
  266      t_diff_beta_0ij(
i) = (ni * nj) * t_node_diff_ksi[i0](
i) +
 
  267                           (n0 * nj) * t_node_diff_ksi[vert_i](
i) +
 
  268                           (n0 * ni) * t_node_diff_ksi[vert_j](
i);
 
  269      base_polynomials(p, ksi0i, &t_diff_ksi0i(0), psi_l, diff_psi_l, 3);
 
  270      base_polynomials(p, ksi0j, &t_diff_ksi0j(0), psi_m, diff_psi_m, 3);
 
  274      for (; oo <= p - 3; oo++) {
 
  277            &diff_psi_l[0], &diff_psi_l[p + 1], &diff_psi_l[2 * p + 2]};
 
  278        for (
int l = 0; 
l <= oo; 
l++) {
 
  283                    &diff_psi_m[
m], &diff_psi_m[p + 1 + 
m],
 
  284                    &diff_psi_m[2 * p + 2 + 
m]};
 
  285            t_diff_phi_f(
i, 
j) = ((t_psi_l * psi_m[
m]) * t_diff_beta_0ij(
j) +
 
  286                                  (beta_0ij * psi_m[
m]) * t_diff_psi_l(
j) +
 
  287                                  (beta_0ij * t_psi_l) * t_diff_psi_m(
j)) *
 
 
  308    int p, 
double *
N, 
double *diffN, 
double *phi_v_e[6],
 
  309    double *diff_phi_v_e[6], 
int gdim,
 
  310    PetscErrorCode (*base_polynomials)(
int p, 
double s, 
double *diff_s,
 
  311                                       double *
L, 
double *diffL,
 
  314  constexpr int edges_nodes[6][2] = {{0, 1}, {1, 2}, {2, 0},
 
  315                                     {0, 3}, {1, 3}, {2, 3}};
 
  343  double diff_psi_l[3 * (p + 1)];
 
  345  for (
int ee = 0; ee != 6; ee++) {
 
  347        t_coords[edges_nodes[ee][1]](
i) - t_coords[edges_nodes[ee][0]](
i);
 
  348    t_diff_ksi0i(
i) = t_node_diff_ksi[edges_nodes[ee][1]](
i) -
 
  349                      t_node_diff_ksi[edges_nodes[ee][0]](
i);
 
  350    auto t_phi_v_e = getFTensor1FromPtr<3>(phi_v_e[ee]);
 
  351    for (
int ii = 0; ii != gdim; ii++) {
 
  352      const int node_shift = ii * 4;
 
  353      const double ni = 
N[node_shift + edges_nodes[ee][1]];
 
  354      const double n0 = 
N[node_shift + edges_nodes[ee][0]];
 
  355      const double beta_e = ni * n0;
 
  356      const double ksi0i = ni - n0;
 
  357      base_polynomials(p, ksi0i, NULL, psi_l, NULL, 3);
 
  359      for (
int l = 0; 
l <= p - 2; 
l++) {
 
  360        t_phi_v_e(
i) = (beta_e * t_psi_l) * t_tou_e(
i);
 
  368    for (
int ee = 0; ee != 6; ee++) {
 
  370          t_coords[edges_nodes[ee][1]](
i) - t_coords[edges_nodes[ee][0]](
i);
 
  371      t_diff_ksi0i(
i) = t_node_diff_ksi[edges_nodes[ee][1]](
i) -
 
  372                        t_node_diff_ksi[edges_nodes[ee][0]](
i);
 
  374      for (
int ii = 0; ii != gdim; ii++) {
 
  375        const int node_shift = ii * 4;
 
  376        const double ni = 
N[node_shift + edges_nodes[ee][1]];
 
  377        const double n0 = 
N[node_shift + edges_nodes[ee][0]];
 
  378        const double beta_e = ni * n0;
 
  379        const double ksi0i = ni - n0;
 
  380        t_diff_beta_e(
i) = ni * t_node_diff_ksi[edges_nodes[ee][0]](
i) +
 
  381                           t_node_diff_ksi[edges_nodes[ee][1]](
i) * n0;
 
  382        base_polynomials(p, ksi0i, &t_diff_ksi0i(0), psi_l, diff_psi_l, 3);
 
  385            diff_psi_l, &diff_psi_l[p + 1], &diff_psi_l[2 * p + 2]);
 
  386        for (
int l = 0; 
l <= p - 2; 
l++) {
 
  387          t_diff_phi_v_e(
i, 
j) =
 
  388              (t_diff_beta_e(
j) * t_psi_l + beta_e * t_diff_psi_l(
j)) *
 
 
  402    int p, 
double *
N, 
double *diffN, 
double *phi_v_f[], 
double *diff_phi_v_f[],
 
  404    PetscErrorCode (*base_polynomials)(
int p, 
double s, 
double *diff_s,
 
  405                                       double *
L, 
double *diffL,
 
  408  constexpr int faces_nodes[4][3] = {
 
  409      {0, 1, 3}, {1, 2, 3}, {0, 3, 2}, {0, 2, 1}};
 
  432  for (
int ff = 0; ff != 4; ff++) {
 
  433    const int v0 = faces_nodes[ff][0];
 
  434    const int vi = faces_nodes[ff][1];
 
  435    const int vj = faces_nodes[ff][2];
 
  436    t_tau0i[ff](
i) = t_coords[vi](
i) - t_coords[v0](
i);
 
  437    t_tau0j[ff](
i) = t_coords[vj](
i) - t_coords[v0](
i);
 
  438    t_diff_ksi0i[ff](
i) = t_node_diff_ksi[vi](
i) - t_node_diff_ksi[v0](
i);
 
  439    t_diff_ksi0j[ff](
i) = t_node_diff_ksi[vj](
i) - t_node_diff_ksi[v0](
i);
 
  442  double psi_l[p + 1], psi_m[p + 1];
 
  443  double diff_psi_l[3 * (p + 1)], diff_psi_m[3 * (p + 1)];
 
  444  for (
int ff = 0; ff != 4; ff++) {
 
  445    const int v0 = faces_nodes[ff][0];
 
  446    const int vi = faces_nodes[ff][1];
 
  447    const int vj = faces_nodes[ff][2];
 
  448    auto t_phi_v_f = getFTensor1FromPtr<3>(phi_v_f[ff]);
 
  450    for (
int ii = 0; ii < gdim; ii++) {
 
  451      const int node_shift = 4 * ii;
 
  452      const double n0 = 
N[node_shift + v0];
 
  453      const double ni = 
N[node_shift + vi];
 
  454      const double nj = 
N[node_shift + vj];
 
  455      const double beta_f = n0 * ni * nj;
 
  457      t_diff_beta_f(
i) = (ni * nj) * t_node_diff_ksi[v0](
i) +
 
  458                         (n0 * nj) * t_node_diff_ksi[vi](
i) +
 
  459                         (n0 * ni) * t_node_diff_ksi[vj](
i);
 
  460      const double ksi0i = ni - n0;
 
  461      const double ksi0j = nj - n0;
 
  462      base_polynomials(p, ksi0i, &t_diff_ksi0i[ff](0), psi_l, diff_psi_l, 3);
 
  463      base_polynomials(p, ksi0j, &t_diff_ksi0j[ff](0), psi_m, diff_psi_m, 3);
 
  466      for (
int oo = 0; oo <= p - 3; oo++) {
 
  469            diff_psi_l, &diff_psi_l[p + 1], &diff_psi_l[2 * p + 2]);
 
  470        for (
int l = 0; 
l <= oo; 
l++) {
 
  474                diff_psi_m[
m], diff_psi_m[p + 1 + 
m],
 
  475                diff_psi_m[2 * p + 2 + 
m]};
 
  476            const double a = beta_f * t_psi_l * psi_m[
m];
 
  477            t_phi_v_f(
i) = 
a * t_tau0i[ff](
i);
 
  480            t_phi_v_f(
i) = 
a * t_tau0j[ff](
i);
 
  484            t_diff_a(
j) = (t_psi_l * psi_m[
m]) * t_diff_beta_f(
j) +
 
  485                          (beta_f * psi_m[
m]) * t_diff_psi_l(
j) +
 
  486                          (beta_f * t_psi_l) * t_diff_psi_m(
j);
 
  487            t_diff_phi_v_f(
i, 
j) = t_diff_a(
j) * t_tau0i[ff](
i);
 
  489            t_diff_phi_v_f(
i, 
j) = t_diff_a(
j) * t_tau0j[ff](
i);
 
  499                 "wrong order %d != %d", jj,
 
 
  508    int p, 
double *
N, 
double *diffN, 
double *phi_v, 
double *diff_phi_v,
 
  510    PetscErrorCode (*base_polynomials)(
int p, 
double s, 
double *diff_s,
 
  511                                       double *
L, 
double *diffL,
 
  534  t_diff_ksi0i(
i) = t_node_diff_ksi[1](
i) - t_node_diff_ksi[0](
i);
 
  535  t_diff_ksi0j(
i) = t_node_diff_ksi[2](
i) - t_node_diff_ksi[0](
i);
 
  536  t_diff_ksi0k(
i) = t_node_diff_ksi[3](
i) - t_node_diff_ksi[0](
i);
 
  539  double diff_psi_l[3 * (p + 1)];
 
  541  double diff_psi_m[3 * (p + 1)];
 
  543  double diff_psi_n[3 * (p + 1)];
 
  545  auto t_phi_v = getFTensor1FromPtr<3>(phi_v);
 
  549  for (
int ii = 0; ii < gdim; ii++) {
 
  550    const int node_shift = ii * 4;
 
  551    const double n0 = 
N[node_shift + 0];
 
  552    const double ni = 
N[node_shift + 1];
 
  553    const double nj = 
N[node_shift + 2];
 
  554    const double nk = 
N[node_shift + 3];
 
  555    const double ksi0i = ni - n0;
 
  556    const double ksi0j = nj - n0;
 
  557    const double ksi0k = nk - n0;
 
  558    const double beta_v = n0 * ni * nj * nk;
 
  559    t_diff_beta_v(
i) = (ni * nj * nk) * t_node_diff_ksi[0](
i) +
 
  560                       (n0 * nj * nk) * t_node_diff_ksi[1](
i) +
 
  561                       (n0 * ni * nk) * t_node_diff_ksi[2](
i) +
 
  562                       (n0 * ni * nj) * t_node_diff_ksi[3](
i);
 
  563    base_polynomials(p, ksi0i, &t_diff_ksi0i(0), psi_l, diff_psi_l, 3);
 
  564    base_polynomials(p, ksi0j, &t_diff_ksi0j(0), psi_m, diff_psi_m, 3);
 
  565    base_polynomials(p, ksi0k, &t_diff_ksi0k(0), psi_n, diff_psi_n, 3);
 
  570    for (
int oo = 0; oo <= p - 4; oo++) {
 
  573          diff_psi_l, &diff_psi_l[p + 1], &diff_psi_l[2 * p + 2]);
 
  574      for (
int l = 0; 
l <= oo; 
l++) {
 
  577            diff_psi_m, &diff_psi_m[p + 1], &diff_psi_m[2 * p + 2]);
 
  578        for (
int m = 0; (
l + 
m) <= oo; 
m++) {
 
  582                                                     diff_psi_n[p + 1 + 
n],
 
  583                                                     diff_psi_n[2 * p + 2 + 
n]);
 
  584            const double a = beta_v * t_psi_l * t_psi_m * psi_n[
n];
 
  597            t_diff_a(
j) = (t_psi_l * t_psi_m * psi_n[
n]) * t_diff_beta_v(
j) +
 
  598                          (beta_v * t_psi_m * psi_n[
n]) * t_diff_psi_l(
j) +
 
  599                          (beta_v * t_psi_l * psi_n[
n]) * t_diff_psi_m(
j) +
 
  600                          (beta_v * t_psi_l * t_psi_m) * t_diff_psi_n(
j);
 
  601            t_diff_phi_v(N0, 
j) = t_diff_a(
j);
 
  602            t_diff_phi_v(N1, 
j) = 0;
 
  603            t_diff_phi_v(N2, 
j) = 0;
 
  605            t_diff_phi_v(N0, 
j) = 0;
 
  606            t_diff_phi_v(N1, 
j) = t_diff_a(
j);
 
  607            t_diff_phi_v(N2, 
j) = 0;
 
  609            t_diff_phi_v(N0, 
j) = 0;
 
  610            t_diff_phi_v(N1, 
j) = 0;
 
  611            t_diff_phi_v(N2, 
j) = t_diff_a(
j);
 
  625               "wrong order %d != %d", jj,
 
 
  635                                         double *diffN, 
double *phi_f,
 
  636                                         double *diff_phi_f, 
int gdim, 
int nb) {
 
  637  const int face_edges_nodes[3][2] = {{0, 1}, {1, 2}, {2, 0}};
 
  638  const int face_opposite_edges_node[] = {2, 0, 1};
 
  646  FTensor::Tensor2<double, 3, 3> t_diff_cross(0., 0., 0., 0., 0., 0., 0., 0.,
 
  651  const int i0 = faces_nodes[0];
 
  652  const int i1 = faces_nodes[1];
 
  653  const int i2 = faces_nodes[2];
 
  654  const int o[] = {faces_nodes[face_opposite_edges_node[0]],
 
  655                   faces_nodes[face_opposite_edges_node[1]],
 
  656                   faces_nodes[face_opposite_edges_node[2]]};
 
  670    t_diff_cross(
i, 
j) = 0;
 
  671    for (
int ee = 0; ee != 3; ee++) {
 
  672      int ei0 = faces_nodes[face_edges_nodes[ee][0]];
 
  673      int ei1 = faces_nodes[face_edges_nodes[ee][1]];
 
  674      t_cross[ee](0) = t_node_diff_ksi[ei0](1) * t_node_diff_ksi[ei1](2) -
 
  675                       t_node_diff_ksi[ei0](2) * t_node_diff_ksi[ei1](1);
 
  676      t_cross[ee](1) = t_node_diff_ksi[ei0](2) * t_node_diff_ksi[ei1](0) -
 
  677                       t_node_diff_ksi[ei0](0) * t_node_diff_ksi[ei1](2);
 
  678      t_cross[ee](2) = t_node_diff_ksi[ei0](0) * t_node_diff_ksi[ei1](1) -
 
  679                       t_node_diff_ksi[ei0](1) * t_node_diff_ksi[ei1](0);
 
  681          diffN[3 * o[ee] + 0], diffN[3 * o[ee] + 1], diffN[3 * o[ee] + 2]);
 
  682      t_diff_cross(
i, 
j) += t_cross[ee](
i) * t_diff_o(
j);
 
  687    t_diff_n0_p_n1(
i) = t_node_diff_ksi[i0](
i) + t_node_diff_ksi[i1](
i);
 
  688    t_diff_n0_p_n1_p_n2(
i) = t_diff_n0_p_n1(
i) + t_node_diff_ksi[i2](
i);
 
  690    for (
int ee = 0; ee != 3; ee++) {
 
  699  boost::shared_ptr<FTensor::Tensor2<double *, 3, 3>> t_diff_phi_ptr;
 
  701    t_diff_phi_ptr = boost::shared_ptr<FTensor::Tensor2<double *, 3, 3>>(
 
  709  double fi[p + 1], diff_fi[3 * p + 3];
 
  710  double fj[p + 1], diff_fj[3 * p + 3];
 
  711  double tmp_fj[p + 1], tmp_diff_fj[3 * p + 3];
 
  712  for (
int ii = 0; ii != gdim; ii++) {
 
  713    const int shift = ii * nb;
 
  714    double n0 = 
N[shift + i0];
 
  715    double n1 = 
N[shift + i1];
 
  716    double n2 = 
N[shift + i2];
 
  717    double *diff_n1 = (diff_phi_f) ? &t_node_diff_ksi[i1](0) : NULL;
 
  718    double *diff_n0_p_n1 = (diff_phi_f) ? &t_diff_n0_p_n1(0) : NULL;
 
  720                              diff_phi_f ? diff_fi : NULL, 3);
 
  722    for (
int pp = 0; pp <= p; pp++) {
 
  723      double *diff_n2 = (diff_phi_f) ? &t_node_diff_ksi[i2](0) : NULL;
 
  724      double *diff_n0_p_n1_p_n2 = (diff_phi_f) ? &t_diff_n0_p_n1_p_n2(0) : NULL;
 
  726                                diff_n0_p_n1_p_n2, tmp_fj,
 
  727                                diff_phi_f ? tmp_diff_fj : NULL, 3);
 
  731        diff_fj[0 * (p + 1) + pp] = tmp_diff_fj[0 * (pp + 1) + pp];
 
  732        diff_fj[1 * (p + 1) + pp] = tmp_diff_fj[1 * (pp + 1) + pp];
 
  733        diff_fj[2 * (p + 1) + pp] = tmp_diff_fj[2 * (pp + 1) + pp];
 
  736    double no0 = 
N[shift + o[0]];
 
  737    double no1 = 
N[shift + o[1]];
 
  738    double no2 = 
N[shift + o[2]];
 
  740    base0(
i) = no0 * t_cross[0](
i) + no1 * t_cross[1](
i) + no2 * t_cross[2](
i);
 
  742    for (
int oo = 0; oo < p; oo++) {
 
  745                                              &diff_fi[2 * p + 2], 1);
 
  746      for (
int ll = 0; ll <= oo; ll++) {
 
  747        const int mm = oo - ll;
 
  749          const double a = t_fi * fj[mm];
 
  751          t_phi(
i) = 
a * base0(
i);
 
  754                &diff_fj[0 + mm], &diff_fj[p + 1 + mm],
 
  755                &diff_fj[2 * p + 2 + mm], 1);
 
  756            (*t_diff_phi_ptr)(
i, 
j) =
 
  757                a * t_diff_cross(
i, 
j) +
 
  758                (t_diff_fi(
j) * fj[mm] + t_fi * t_diff_fj(
j)) * base0(
i);
 
  770               "wrong number of base functions " 
  771               "jj!=NBFACETRI_DEMKOWICZ_HDIV(p) " 
 
  781    int p, 
double *
N, 
double *diffN, 
int p_f[], 
double *phi_f[4],
 
  782    double *diff_phi_f[4], 
double *phi_v, 
double *diff_phi_v, 
int gdim) {
 
  784  const int opposite_face_node[4] = {2, 0, 1, 3};
 
  786  const int znf[] = {0, 2, 3};
 
  802  for (
int ff = 0; ff != 4; ++ff) {
 
  803    t_m_node_diff_ksi[ff](
i) = -t_node_diff_ksi[ff](
i);
 
  815  for (
int ii = 0; ii != gdim; ii++) {
 
  816    const int shift = 4 * ii;
 
  818    for (
int ff = 0; ff != 3; ff++) {
 
  819      const int fff = znf[ff];
 
  820      const int iO = opposite_face_node[fff];
 
  821      const double nO = 
N[shift + iO];
 
  822      for (
int pp = 1; pp <= p; pp++) {
 
  824            pp, 2 * pp + 2, nO, 1 - nO, &t_node_diff_ksi[iO](0),
 
  825            &t_m_node_diff_ksi[iO](0), &fk(ff, 0), &diff_fk(ff, 0), 3);
 
  830    for (
int oo = 2; oo <= p; oo++) {
 
  831      for (
int k = 1; 
k != oo; 
k++) {
 
  838          int sp[] = {ii * 3 * nb_dofs + 3 * s, ii * 3 * nb_dofs + 3 * s,
 
  839                      ii * 3 * nb_dofs + 3 * s};
 
  842                                            &phi_f[znf[0]][sp[0] + 
HVEC1],
 
  843                                            &phi_f[znf[0]][sp[0] + 
HVEC2], 3),
 
  845                                            &phi_f[znf[1]][sp[1] + 
HVEC1],
 
  846                                            &phi_f[znf[1]][sp[1] + 
HVEC2], 3),
 
  848                                            &phi_f[znf[2]][sp[2] + 
HVEC1],
 
  849                                            &phi_f[znf[2]][sp[2] + 
HVEC2], 3)};
 
  850          int sdp[] = {ii * 9 * nb_dofs + 9 * s, ii * 9 * nb_dofs + 9 * s,
 
  851                       ii * 9 * nb_dofs + 9 * s};
 
  854                  &diff_phi_f[znf[0]][sdp[0] + 
HVEC0_0],
 
  855                  &diff_phi_f[znf[0]][sdp[0] + 
HVEC0_1],
 
  856                  &diff_phi_f[znf[0]][sdp[0] + 
HVEC0_2],
 
  857                  &diff_phi_f[znf[0]][sdp[0] + 
HVEC1_0],
 
  858                  &diff_phi_f[znf[0]][sdp[0] + 
HVEC1_1],
 
  859                  &diff_phi_f[znf[0]][sdp[0] + 
HVEC1_2],
 
  860                  &diff_phi_f[znf[0]][sdp[0] + 
HVEC2_0],
 
  861                  &diff_phi_f[znf[0]][sdp[0] + 
HVEC2_1],
 
  862                  &diff_phi_f[znf[0]][sdp[0] + 
HVEC2_2], 9),
 
  864                  &diff_phi_f[znf[1]][sdp[1] + 
HVEC0_0],
 
  865                  &diff_phi_f[znf[1]][sdp[1] + 
HVEC0_1],
 
  866                  &diff_phi_f[znf[1]][sdp[1] + 
HVEC0_2],
 
  867                  &diff_phi_f[znf[1]][sdp[1] + 
HVEC1_0],
 
  868                  &diff_phi_f[znf[1]][sdp[1] + 
HVEC1_1],
 
  869                  &diff_phi_f[znf[1]][sdp[1] + 
HVEC1_2],
 
  870                  &diff_phi_f[znf[1]][sdp[1] + 
HVEC2_0],
 
  871                  &diff_phi_f[znf[1]][sdp[1] + 
HVEC2_1],
 
  872                  &diff_phi_f[znf[1]][sdp[1] + 
HVEC2_2], 9),
 
  874                  &diff_phi_f[znf[2]][sdp[2] + 
HVEC0_0],
 
  875                  &diff_phi_f[znf[2]][sdp[2] + 
HVEC0_1],
 
  876                  &diff_phi_f[znf[2]][sdp[2] + 
HVEC0_2],
 
  877                  &diff_phi_f[znf[2]][sdp[2] + 
HVEC1_0],
 
  878                  &diff_phi_f[znf[2]][sdp[2] + 
HVEC1_1],
 
  879                  &diff_phi_f[znf[2]][sdp[2] + 
HVEC1_2],
 
  880                  &diff_phi_f[znf[2]][sdp[2] + 
HVEC2_0],
 
  881                  &diff_phi_f[znf[2]][sdp[2] + 
HVEC2_1],
 
  882                  &diff_phi_f[znf[2]][sdp[2] + 
HVEC2_2], 9)};
 
  884            for (
int ff = 0; ff != 3; ff++) {
 
  886                                                    diff_fk(ff, 1 * p + 
k - 1),
 
  887                                                    diff_fk(ff, 2 * p + 
k - 1));
 
  888              t_phi_v(
i) = fk(ff, 
k - 1) * t_phi_f[ff](
i);
 
  889              t_diff_phi_v(
i, 
j) = t_diff_fk(
j) * t_phi_f[ff](
i) +
 
  890                                   fk(ff, 
k - 1) * t_diff_phi_f[ff](
i, 
j);
 
  903               "wrong number of base functions " 
  904               "jj!=NBVOLUMETET_DEMKOWICZ_HDIV(p) " 
 
#define FTENSOR_INDEX(DIM, I)
PetscErrorCode IntegratedJacobi_polynomials(int p, double alpha, double x, double t, double *diff_x, double *diff_t, double *L, double *diffL, const int dim)
Calculate integrated Jacobi approximation basis.
PetscErrorCode Jacobi_polynomials(int p, double alpha, double x, double t, double *diff_x, double *diff_t, double *L, double *diffL, const int dim)
Calculate Jacobi approximation basis.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define NBVOLUMETET_DEMKOWICZ_HDIV(P)
#define NBVOLUMETET_AINSWORTH_FACE_HDIV(P)
#define NBFACETRI_DEMKOWICZ_HDIV(P)
#define NBFACETRI_AINSWORTH_FACE_HDIV(P)
#define NBVOLUMETET_AINSWORTH_VOLUME_HDIV(P)
#define NBFACETRI_AINSWORTH_EDGE_HDIV(P)
FTensor::Index< 'i', SPACE_DIM > i
const double n
refractive index of diffusive medium
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > getFTensor2HVecFromPtr< 3, 3 >(double *ptr)
MoFEMErrorCode Hdiv_Ainsworth_FaceBubbleShapeFunctions(int *faces_nodes, int *p, double *N, double *diffN, double *phi_f[], double *diff_phi_f[], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Face bubble functions by Ainsworth .
MoFEMErrorCode Hdiv_Demkowicz_Face_MBTET_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f, double *diff_phi_f, int gdim, int nb)
MoFEMErrorCode Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET(int *faces_nodes, int *p, double *N, double *diffN, double *phi_f_e[4][3], double *diff_phi_f_e[4][3], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Hdiv base functions, Edge-based face functions by Ainsworth .
MoFEMErrorCode Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f_e[3], double *diff_phi_f_e[3], int gdim, int nb, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Hdiv base functions, Edge-based face functions by Ainsworth .
MoFEMErrorCode Hdiv_Demkowicz_Interior_MBTET(int p, double *N, double *diffN, int p_face[], double *phi_f[4], double *diff_phi_f[4], double *phi_v, double *diff_phi_v, int gdim)
MoFEMErrorCode Hdiv_Ainsworth_FaceBasedVolumeShapeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v_f[], double *diff_phi_v_f[], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
MoFEMErrorCode Hdiv_Ainsworth_FaceBubbleShapeFunctions_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f, double *diff_phi_f, int gdim, int nb, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Face bubble functions by Ainsworth .
MoFEMErrorCode Hdiv_Ainsworth_EdgeBasedVolumeShapeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v_e[6], double *diff_phi_v_e[6], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Hdiv base function, Edge-based interior (volume) functions by Ainsworth .
MoFEMErrorCode Hdiv_Ainsworth_VolumeBubbleShapeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v, double *diff_phi_v, int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Interior bubble functions by Ainsworth .
FTensor::Index< 'm', 3 > m
static boost::function< int(int)> broken_nbvolumetet_edge_hdiv
static boost::function< int(int)> broken_nbvolumetet_face_hdiv
static boost::function< int(int)> broken_nbfacetri_face_hdiv
static boost::function< int(int)> broken_nbvolumetet_volume_hdiv
static boost::function< int(int)> broken_nbfacetri_edge_hdiv