37      boost::multi_index::indexed_by<
 
   39          boost::multi_index::hashed_unique<
 
   44                  member<BaseCacheItem, int, &BaseCacheItem::order>,
 
   45                  member<BaseCacheItem, int, &BaseCacheItem::nb_gauss_pts>>>>
 
   51      boost::multi_index::indexed_by<
 
   53          boost::multi_index::hashed_unique<
 
   58                  member<HDivBaseCacheItem, int, &HDivBaseCacheItem::order>,
 
   61                  member<HDivBaseCacheItem, int, &HDivBaseCacheItem::n0>,
 
   62                  member<HDivBaseCacheItem, int, &HDivBaseCacheItem::n1>,
 
   63                  member<HDivBaseCacheItem, int, &HDivBaseCacheItem::n2>>>>
 
   67  static std::array<std::map<const void *, BaseCacheMI>, 
LASTBASE>
 
   69  static std::array<std::map<const void *, BaseCacheMI>, 
LASTBASE>
 
   71  static std::array<std::map<const void *, HDivBaseFaceCacheMI>, 
LASTBASE>
 
   73  static std::array<std::map<const void *, BaseCacheMI>, 
LASTBASE>
 
 
   77std::array<std::map<const void *, TetBaseCache::BaseCacheMI>, 
LASTBASE>
 
   79std::array<std::map<const void *, TetBaseCache::BaseCacheMI>, 
LASTBASE>
 
   81std::array<std::map<const void *, TetBaseCache::HDivBaseFaceCacheMI>, 
LASTBASE>
 
   83std::array<std::map<const void *, TetBaseCache::BaseCacheMI>, 
LASTBASE>
 
   87TetPolynomialBase::query_interface(boost::typeindex::type_index type_index,
 
 
  100    auto erase = [&](
auto cache) {
 
  101      if (cache.find(
vPtr) != cache.end())
 
  105    for (
auto b = 0; b != 
LASTBASE; ++b) {
 
 
  136  PetscErrorCode (*base_polynomials)(
int p, 
double s, 
double *diff_s, 
double *
L,
 
  137                                     double *diffL, 
const int dim) =
 
  140  int nb_gauss_pts = pts.size2();
 
  142  int sense[6], 
order[6];
 
  148    double *h1_edge_n[6], *diff_h1_egde_n[6];
 
  149    for (
int ee = 0; ee != 6; ++ee) {
 
  152                "data inconsistency");
 
  157      data.
dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, nb_dofs,
 
  159      data.
dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
 
  170        h1_edge_n, diff_h1_egde_n, nb_gauss_pts, base_polynomials);
 
  172    for (
int ee = 0; ee != 6; ++ee) {
 
  174      data.
dataOnEntities[MBEDGE][ee].getDiffN(base).resize(0, 0, 
false);
 
  183    double *h1_face_n[4], *diff_h1_face_n[4];
 
  184    for (
int ff = 0; ff != 4; ++ff) {
 
  187                "data inconsistency");
 
  191      data.
dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, nb_dofs,
 
  193      data.
dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
 
  210        h1_face_n, diff_h1_face_n, nb_gauss_pts, base_polynomials);
 
  213    for (
int ff = 0; ff != 4; ++ff) {
 
  215      data.
dataOnEntities[MBTRI][ff].getDiffN(base).resize(0, 0, 
false);
 
  223    data.
dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, nb_vol_dofs,
 
  226                                                        3 * nb_vol_dofs, 
false);
 
  233        nb_gauss_pts, base_polynomials);
 
 
  248  const int nb_gauss_pts = pts.size2();
 
  251      (
unsigned int)nb_gauss_pts)
 
  253            "Base functions or nodes has wrong number of integration points " 
  259    auto &ptr = data.getBBAlphaIndicesSharedPtr(
field_name);
 
  273    auto &ptr = data.getBBDiffNSharedPtr(
field_name);
 
  279  auto get_alpha_by_name_ptr =
 
  281         const std::string &
field_name) -> boost::shared_ptr<MatrixInt> & {
 
  282    return data.getBBAlphaIndicesSharedPtr(
field_name);
 
  285  auto get_base_by_name_ptr =
 
  287         const std::string &
field_name) -> boost::shared_ptr<MatrixDouble> & {
 
  291  auto get_diff_base_by_name_ptr =
 
  293         const std::string &
field_name) -> boost::shared_ptr<MatrixDouble> & {
 
  297  auto get_alpha_by_order_ptr =
 
  298      [](
auto &data, 
const size_t o) -> boost::shared_ptr<MatrixInt> & {
 
  299    return data.getBBAlphaIndicesByOrderSharedPtr(o);
 
  302  auto get_base_by_order_ptr =
 
  303      [](
auto &data, 
const size_t o) -> boost::shared_ptr<MatrixDouble> & {
 
  304    return data.getBBNByOrderSharedPtr(o);
 
  307  auto get_diff_base_by_order_ptr =
 
  308      [](
auto &data, 
const size_t o) -> boost::shared_ptr<MatrixDouble> & {
 
  309    return data.getBBDiffNByOrderSharedPtr(o);
 
  313  auto &vertex_alpha = get_alpha(vert_ent_data);
 
  314  vertex_alpha.resize(4, 4, 
false);
 
  315  vertex_alpha.clear();
 
  316  for (
int n = 0; 
n != 4; ++
n)
 
  319  auto &vert_get_n = get_base(vert_ent_data);
 
  320  auto &vert_get_diff_n = get_diff_base(vert_ent_data);
 
  321  vert_get_n.resize(nb_gauss_pts, 4, 
false);
 
  322  vert_get_diff_n.resize(nb_gauss_pts, 12, 
false);
 
  324      1, 
lambda.size1(), vertex_alpha.size1(), &vertex_alpha(0, 0),
 
  326      &vert_get_diff_n(0, 0));
 
  327  for (
int n = 0; 
n != 4; ++
n) {
 
  328    const double f = boost::math::factorial<double>(
 
  330    for (
int g = 0; 
g != nb_gauss_pts; ++
g) {
 
  331      vert_get_n(
g, 
n) *= f;
 
  332      for (
int d = 0; d != 3; ++d)
 
  333        vert_get_diff_n(
g, 3 * 
n + d) *= f;
 
  341              "Wrong size of ent data");
 
  343    constexpr int edges_nodes[6][2] = {{0, 1}, {1, 2}, {2, 0},
 
  344                                       {0, 3}, {1, 3}, {2, 3}};
 
  345    for (
int ee = 0; ee != 6; ++ee) {
 
  347      const int sense = ent_data.getSense();
 
  350                "Sense of the edge unknown");
 
  351      const int order = ent_data.getOrder();
 
  355        if (get_alpha_by_order_ptr(ent_data, 
order)) {
 
  357              get_alpha_by_order_ptr(ent_data, 
order);
 
  359              get_base_by_order_ptr(ent_data, 
order);
 
  360          get_diff_base_by_name_ptr(ent_data, 
field_name) =
 
  361              get_diff_base_by_order_ptr(ent_data, 
order);
 
  363          auto &get_n = get_base(ent_data);
 
  364          auto &get_diff_n = get_diff_base(ent_data);
 
  365          get_n.resize(nb_gauss_pts, nb_dofs, 
false);
 
  366          get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, 
false);
 
  369          edge_alpha.resize(nb_dofs, 4, 
false);
 
  373            for (
int i = 0; 
i != edge_alpha.size1(); ++
i) {
 
  374              int a = edge_alpha(
i, edges_nodes[ee][0]);
 
  375              edge_alpha(
i, edges_nodes[ee][0]) =
 
  376                  edge_alpha(
i, edges_nodes[ee][1]);
 
  377              edge_alpha(
i, edges_nodes[ee][1]) = 
a;
 
  381              order, 
lambda.size1(), edge_alpha.size1(), &edge_alpha(0, 0),
 
  385          get_alpha_by_order_ptr(ent_data, 
order) =
 
  387          get_base_by_order_ptr(ent_data, 
order) =
 
  389          get_diff_base_by_order_ptr(ent_data, 
order) =
 
  390              get_diff_base_by_name_ptr(ent_data, 
field_name);
 
  395    for (
int ee = 0; ee != 6; ++ee) {
 
  397      ent_data.getBBAlphaIndicesSharedPtr(
field_name).reset();
 
  398      auto &get_n = get_base(ent_data);
 
  399      auto &get_diff_n = get_diff_base(ent_data);
 
  400      get_n.resize(nb_gauss_pts, 0, 
false);
 
  401      get_diff_n.resize(nb_gauss_pts, 0, 
false);
 
  409              "Wrong size of ent data");
 
  415    for (
int ff = 0; ff != 4; ++ff) {
 
  417      const int order = ent_data.getOrder();
 
  421        if (get_alpha_by_order_ptr(ent_data, 
order)) {
 
  423              get_alpha_by_order_ptr(ent_data, 
order);
 
  425              get_base_by_order_ptr(ent_data, 
order);
 
  426          get_diff_base_by_name_ptr(ent_data, 
field_name) =
 
  427              get_diff_base_by_order_ptr(ent_data, 
order);
 
  430          auto &get_n = get_base(ent_data);
 
  431          auto &get_diff_n = get_diff_base(ent_data);
 
  432          get_n.resize(nb_gauss_pts, nb_dofs, 
false);
 
  433          get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, 
false);
 
  435          auto &face_alpha = get_alpha(ent_data);
 
  436          face_alpha.resize(nb_dofs, 4, 
false);
 
  440          senseFaceAlpha.resize(face_alpha.size1(), face_alpha.size2(), 
false);
 
  442          constexpr int tri_nodes[4][3] = {
 
  443              {0, 1, 3}, {1, 2, 3}, {0, 2, 3}, {0, 1, 2}};
 
  444          for (
int d = 0; d != nb_dofs; ++d)
 
  445            for (
int n = 0; 
n != 3; ++
n)
 
  447                  face_alpha(d, tri_nodes[ff][
n]);
 
  450              order, 
lambda.size1(), face_alpha.size1(), &face_alpha(0, 0),
 
  454          get_alpha_by_order_ptr(ent_data, 
order) =
 
  456          get_base_by_order_ptr(ent_data, 
order) =
 
  458          get_diff_base_by_order_ptr(ent_data, 
order) =
 
  459              get_diff_base_by_name_ptr(ent_data, 
field_name);
 
  464    for (
int ff = 0; ff != 4; ++ff) {
 
  466      ent_data.getBBAlphaIndicesSharedPtr(
field_name).reset();
 
  467      auto &get_n = get_base(ent_data);
 
  468      auto &get_diff_n = get_diff_base(ent_data);
 
  469      get_n.resize(nb_gauss_pts, 0, 
false);
 
  470      get_diff_n.resize(nb_gauss_pts, 0, 
false);
 
  477              "Wrong size ent of ent data");
 
  480    const int order = ent_data.getOrder();
 
  482    if (get_alpha_by_order_ptr(ent_data, 
order)) {
 
  484          get_alpha_by_order_ptr(ent_data, 
order);
 
  486          get_base_by_order_ptr(ent_data, 
order);
 
  487      get_diff_base_by_name_ptr(ent_data, 
field_name) =
 
  488          get_diff_base_by_order_ptr(ent_data, 
order);
 
  491      auto &get_n = get_base(ent_data);
 
  492      auto &get_diff_n = get_diff_base(ent_data);
 
  493      get_n.resize(nb_gauss_pts, nb_dofs, 
false);
 
  494      get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, 
false);
 
  496        auto &tet_alpha = get_alpha(ent_data);
 
  497        tet_alpha.resize(nb_dofs, 4, 
false);
 
  501            order, 
lambda.size1(), tet_alpha.size1(), &tet_alpha(0, 0),
 
  505        get_alpha_by_order_ptr(ent_data, 
order) =
 
  507        get_base_by_order_ptr(ent_data, 
order) =
 
  509        get_diff_base_by_order_ptr(ent_data, 
order) =
 
  510            get_diff_base_by_name_ptr(ent_data, 
field_name);
 
  515    ent_data.getBBAlphaIndicesSharedPtr(
field_name).reset();
 
  516    auto &get_n = get_base(ent_data);
 
  517    auto &get_diff_n = get_diff_base(ent_data);
 
  518    get_n.resize(nb_gauss_pts, 0, 
false);
 
  519    get_diff_n.resize(nb_gauss_pts, 0, 
false);
 
 
  549  PetscErrorCode (*base_polynomials)(
int p, 
double s, 
double *diff_s, 
double *
L,
 
  550                                     double *diffL, 
const int dim) =
 
  553  int nb_gauss_pts = pts.size2();
 
  557  data.
dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, nb_dofs, 
false);
 
  558  data.
dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 3 * nb_dofs,
 
  574  auto interior_cache_ptr = get_interior_cache(base);
 
  576  if (interior_cache_ptr) {
 
  578        interior_cache_ptr->find(boost::make_tuple(volume_order, nb_gauss_pts));
 
  579    if (it != interior_cache_ptr->end()) {
 
  581      noalias(data.
dataOnEntities[MBTET][0].getDiffN(base)) = it->diffN;
 
  592      nb_gauss_pts, base_polynomials);
 
  594  if (interior_cache_ptr) {
 
  595    auto p = interior_cache_ptr->emplace(
 
 
  610  const int nb_gauss_pts = pts.size2();
 
  613      (
unsigned int)nb_gauss_pts)
 
  615             "Base functions or nodes has wrong number of integration points " 
  621    auto &ptr = data.getBBAlphaIndicesSharedPtr(
field_name);
 
  635    auto &ptr = data.getBBDiffNSharedPtr(
field_name);
 
  641  auto get_alpha_by_name_ptr =
 
  643         const std::string &
field_name) -> boost::shared_ptr<MatrixInt> & {
 
  644    return data.getBBAlphaIndicesSharedPtr(
field_name);
 
  647  auto get_base_by_name_ptr =
 
  649         const std::string &
field_name) -> boost::shared_ptr<MatrixDouble> & {
 
  653  auto get_diff_base_by_name_ptr =
 
  655         const std::string &
field_name) -> boost::shared_ptr<MatrixDouble> & {
 
  659  auto get_alpha_by_order_ptr =
 
  660      [](
auto &data, 
const size_t o) -> boost::shared_ptr<MatrixInt> & {
 
  661    return data.getBBAlphaIndicesByOrderSharedPtr(o);
 
  664  auto get_base_by_order_ptr =
 
  665      [](
auto &data, 
const size_t o) -> boost::shared_ptr<MatrixDouble> & {
 
  666    return data.getBBNByOrderSharedPtr(o);
 
  669  auto get_diff_base_by_order_ptr =
 
  670      [](
auto &data, 
const size_t o) -> boost::shared_ptr<MatrixDouble> & {
 
  671    return data.getBBDiffNByOrderSharedPtr(o);
 
  677              "Wrong size ent of ent data");
 
  680    const int order = ent_data.getOrder();
 
  683    if (get_alpha_by_order_ptr(ent_data, 
order)) {
 
  685          get_alpha_by_order_ptr(ent_data, 
order);
 
  687          get_base_by_order_ptr(ent_data, 
order);
 
  688      get_diff_base_by_name_ptr(ent_data, 
field_name) =
 
  689          get_diff_base_by_order_ptr(ent_data, 
order);
 
  692      auto &get_n = get_base(ent_data);
 
  693      auto &get_diff_n = get_diff_base(ent_data);
 
  694      get_n.resize(nb_gauss_pts, nb_dofs, 
false);
 
  695      get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, 
false);
 
  703                    "Inconsistent number of DOFs");
 
  705          auto &tri_alpha = get_alpha(ent_data);
 
  716                    "Inconsistent number of DOFs");
 
  718          auto &tet_alpha = get_alpha(ent_data);
 
  719          tet_alpha.resize(nb_dofs, 4, 
false);
 
  725            std::array<int *, 6> tet_edge_ptr{
 
  733                                                           tet_edge_ptr.data());
 
  736              std::array<int *, 6> tet_face_ptr{
 
  746                  face_n.data(), tet_face_ptr.data());
 
  756              order, 
lambda.size1(), tet_alpha.size1(), &tet_alpha(0, 0),
 
  760          get_alpha_by_order_ptr(ent_data, 
order) =
 
  762          get_base_by_order_ptr(ent_data, 
order) =
 
  764          get_diff_base_by_order_ptr(ent_data, 
order) =
 
  765              get_diff_base_by_name_ptr(ent_data, 
field_name);
 
  771    ent_data.getBBAlphaIndicesSharedPtr(
field_name).reset();
 
  772    auto &get_n = get_base(ent_data);
 
  773    auto &get_diff_n = get_diff_base(ent_data);
 
  774    get_n.resize(nb_gauss_pts, 0, 
false);
 
  775    get_diff_n.resize(nb_gauss_pts, 0, 
false);
 
 
  787    int volume_order, std::array<int, 4> &faces_order,
 
  788    std::array<int, 3 * 4> &faces_nodes,
 
  790    boost::function<
int(
int)> broken_nbfacetri_edge_hdiv,
 
  791    boost::function<
int(
int)> broken_nbfacetri_face_hdiv,
 
  792    boost::function<
int(
int)> broken_nbvolumetet_edge_hdiv,
 
  793    boost::function<
int(
int)> broken_nbvolumetet_face_hdiv,
 
  794    boost::function<
int(
int)> broken_nbvolumetet_volume_hdiv
 
  799  PetscErrorCode (*base_polynomials)(
int p, 
double s, 
double *diff_s, 
double *
L,
 
  800                                     double *diffL, 
const int dim) =
 
  803  int nb_gauss_pts = pts.size2();
 
  807  double *phi_f_e[4][3];
 
  809  double *diff_phi_f_e[4][3];
 
  810  double *diff_phi_f[4];
 
  817  for (
int ff = 0; ff != 4; ++ff) {
 
  819        broken_nbfacetri_edge_hdiv(faces_order[ff]));
 
  821    for (
int ee = 0; ee < 3; ee++) {
 
  822      N_face_edge(ff, ee).resize(nb_gauss_pts, 3 * face_edge_dofs, 
false);
 
  823      diffN_face_edge(ff, ee).resize(nb_gauss_pts, 9 * face_edge_dofs, 
false);
 
  824      phi_f_e[ff][ee] = &*
N_face_edge(ff, ee).data().begin();
 
  828        broken_nbfacetri_face_hdiv(faces_order[ff]));
 
  829    N_face_bubble[ff].resize(nb_gauss_pts, 3 * face_bubble_dofs, 
false);
 
  835  constexpr int nb_nodes_on_tet = 4;
 
  837  for (
int ff = 0; ff < 4; ff++) {
 
  839        &faces_nodes[3 * ff], broken_nbfacetri_edge_hdiv(faces_order[ff]),
 
  840        &*shape_functions.data().begin(), &*diff_shape_functions.data().begin(),
 
  841        phi_f_e[ff], diff_phi_f_e[ff], nb_gauss_pts, nb_nodes_on_tet,
 
  845  for (
int ff = 0; ff < 4; ff++) {
 
  847        &faces_nodes[3 * ff], broken_nbfacetri_face_hdiv(faces_order[ff]),
 
  848        &*shape_functions.data().begin(), &*diff_shape_functions.data().begin(),
 
  849        phi_f[ff], diff_phi_f[ff], nb_gauss_pts, nb_nodes_on_tet,
 
  858  double *diff_phi_v_e[6];
 
  859  double *diff_phi_v_f[4];
 
  863      broken_nbvolumetet_edge_hdiv(volume_order));
 
  866  for (
int ee = 0; ee != 6; ++ee) {
 
  867    N_volume_edge[ee].resize(nb_gauss_pts, 3 * volume_edge_dofs, 
false);
 
  872  if (volume_edge_dofs)
 
  874        broken_nbvolumetet_edge_hdiv(volume_order),
 
  875        &*shape_functions.data().begin(), &*diff_shape_functions.data().begin(),
 
  876        phi_v_e, diff_phi_v_e, nb_gauss_pts, base_polynomials);
 
  879      broken_nbvolumetet_face_hdiv(volume_order));
 
  882  for (
int ff = 0; ff != 4; ++ff) {
 
  883    N_volume_face[ff].resize(nb_gauss_pts, 3 * volume_face_dofs, 
false);
 
  888  if (volume_face_dofs)
 
  890        broken_nbvolumetet_face_hdiv(volume_order),
 
  891        &*shape_functions.data().begin(), &*diff_shape_functions.data().begin(),
 
  892        phi_v_f, diff_phi_v_f, nb_gauss_pts, base_polynomials);
 
  895      broken_nbvolumetet_volume_hdiv(volume_order));
 
  900  if (volume_bubble_dofs)
 
  902        broken_nbvolumetet_volume_hdiv(volume_order),
 
  903        &*shape_functions.data().begin(), &*diff_shape_functions.data().begin(),
 
  904        phi_v, diff_phi_v, nb_gauss_pts, base_polynomials);
 
 
  912  std::array<int, 4> faces_order;
 
  913  std::array<int, 4 * 3> faces_nodes;
 
  920            faces_nodes.begin());
 
  921  for (
int ff = 0; ff != 4; ff++) {
 
  928      faces_order, faces_nodes,
 
  946  int nb_gauss_pts = pts.size2();
 
  966      getFTensor1FromPtr<3>(&*(
N_face_edge(0, 0).data().begin())),
 
  967      getFTensor1FromPtr<3>(&*(
N_face_edge(0, 1).data().begin())),
 
  968      getFTensor1FromPtr<3>(&*(
N_face_edge(0, 2).data().begin())),
 
  969      getFTensor1FromPtr<3>(&*(
N_face_edge(1, 0).data().begin())),
 
  970      getFTensor1FromPtr<3>(&*(
N_face_edge(1, 1).data().begin())),
 
  971      getFTensor1FromPtr<3>(&*(
N_face_edge(1, 2).data().begin())),
 
  972      getFTensor1FromPtr<3>(&*(
N_face_edge(2, 0).data().begin())),
 
  973      getFTensor1FromPtr<3>(&*(
N_face_edge(2, 1).data().begin())),
 
  974      getFTensor1FromPtr<3>(&*(
N_face_edge(2, 2).data().begin())),
 
  975      getFTensor1FromPtr<3>(&*(
N_face_edge(3, 0).data().begin())),
 
  976      getFTensor1FromPtr<3>(&*(
N_face_edge(3, 1).data().begin())),
 
  977      getFTensor1FromPtr<3>(&*(
N_face_edge(3, 2).data().begin()))};
 
  992  for (
int ff = 0; ff != 4; ff++) {
 
  993    int face_order = faces_order[ff];
 
 1000                                                     3 * face_dofs, 
false);
 
 1001    data.
dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
 
 1002                                                         9 * face_dofs, 
false);
 
 1006      double *diff_base_ptr =
 
 1008      auto t_base = getFTensor1FromPtr<3>(base_ptr);
 
 1011      auto max_face_order =
 
 1012          std::max(face_order,
 
 1015          std::max(max_face_order,
 
 1018      for (
int gg = 0; gg != nb_gauss_pts; gg++) {
 
 1019        for (
int oo = 0; oo != max_face_order; oo++) {
 
 1025              for (
int ee = 0; ee != 3; ++ee) {
 
 1026                t_base(
i) = t_base_f_e[ff * 3 + ee](
i);
 
 1028                ++t_base_f_e[ff * 3 + ee];
 
 1030              for (
int ee = 0; ee != 3; ++ee) {
 
 1031                t_diff_base(
i, 
j) = t_diff_base_f_e[ff * 3 + ee](
i, 
j);
 
 1033                ++t_diff_base_f_e[ff * 3 + ee];
 
 1041              t_base(
i) = t_base_f_f[ff](
i);
 
 1044              t_diff_base(
i, 
j) = t_diff_base_f_f[ff](
i, 
j);
 
 1046              ++t_diff_base_f_f[ff];
 
 1061  data.
dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 3 * volume_dofs,
 
 1063  data.
dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
 
 1064                                                      9 * volume_dofs, 
false);
 
 1068    double *diff_base_ptr =
 
 1070    auto t_base = getFTensor1FromPtr<3>(base_ptr);
 
 1104    auto t_base_v = getFTensor1FromPtr<3>(base_ptr);
 
 1107    auto max_volume_order = std::max(
 
 1110    max_volume_order = std::max(
 
 1113    max_volume_order = std::max(
 
 1117    for (
int gg = 0; gg != nb_gauss_pts; gg++) {
 
 1118      for (
int oo = 0; oo < max_volume_order; oo++) {
 
 1125            for (
int ee = 0; ee < 6; ee++) {
 
 1126              t_base(
i) = t_base_v_e[ee](
i);
 
 1129              t_diff_base(
i, 
j) = t_diff_base_v_e[ee](
i, 
j);
 
 1131              ++t_diff_base_v_e[ee];
 
 1140            for (
int ff = 0; ff < 4; ff++) {
 
 1141              t_base(
i) = t_base_v_f[ff](
i);
 
 1144              t_diff_base(
i, 
j) = t_diff_base_v_f[ff](
i, 
j);
 
 1146              ++t_diff_base_v_f[ff];
 
 1155            t_base(
i) = t_base_v(
i);
 
 1158            t_diff_base(
i, 
j) = t_diff_base_v(
i, 
j);
 
 
 1185  int nb_gauss_pts = pts.size2();
 
 1191  int nb_dofs_volume =
 
 1199  int nb_dofs = 4 * nb_dofs_face + nb_dofs_volume;
 
 1200  data.
dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
 
 1202  data.
dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 9 * nb_dofs,
 
 1217  auto interior_cache_ptr = get_interior_cache(base);
 
 1219  if (interior_cache_ptr) {
 
 1221        interior_cache_ptr->find(boost::make_tuple(volume_order, nb_gauss_pts));
 
 1222    if (it != interior_cache_ptr->end()) {
 
 1224      noalias(data.
dataOnEntities[MBTET][0].getDiffN(base)) = it->diffN;
 
 1229  std::array<int, 4 * 3> faces_nodes = {0, 1, 3, 1, 2, 3, 0, 3, 2, 0, 2, 1};
 
 1230  std::array<int, 4> faces_order{volume_order, volume_order, volume_order,
 
 1235      faces_order, faces_nodes,
 
 1245  auto *base_ptr = &*data.
dataOnEntities[MBTET][0].getN(base).data().begin();
 
 1246  auto *diff_base_ptr =
 
 1248  auto t_base = getFTensor1FromPtr<3>(base_ptr);
 
 1253      getFTensor1FromPtr<3>(&*(
N_face_edge(0, 0).data().begin())),
 
 1254      getFTensor1FromPtr<3>(&*(
N_face_edge(0, 1).data().begin())),
 
 1255      getFTensor1FromPtr<3>(&*(
N_face_edge(0, 2).data().begin())),
 
 1256      getFTensor1FromPtr<3>(&*(
N_face_edge(1, 0).data().begin())),
 
 1257      getFTensor1FromPtr<3>(&*(
N_face_edge(1, 1).data().begin())),
 
 1258      getFTensor1FromPtr<3>(&*(
N_face_edge(1, 2).data().begin())),
 
 1259      getFTensor1FromPtr<3>(&*(
N_face_edge(2, 0).data().begin())),
 
 1260      getFTensor1FromPtr<3>(&*(
N_face_edge(2, 1).data().begin())),
 
 1261      getFTensor1FromPtr<3>(&*(
N_face_edge(2, 2).data().begin())),
 
 1262      getFTensor1FromPtr<3>(&*(
N_face_edge(3, 0).data().begin())),
 
 1263      getFTensor1FromPtr<3>(&*(
N_face_edge(3, 1).data().begin())),
 
 1264      getFTensor1FromPtr<3>(&*(
N_face_edge(3, 2).data().begin()))};
 
 1284      getFTensor1FromPtr<3>(&*(
N_face_bubble[3].data().begin()))};
 
 1322  auto t_base_v = getFTensor1FromPtr<3>(base_vol_ptr);
 
 1326  int count_dofs_face = 0;
 
 1327  int count_dofs_volume = 0;
 
 1329  auto max_volume_order =
 
 1330      std::max(volume_order,
 
 1333      std::max(max_volume_order,
 
 1336      std::max(max_volume_order,
 
 1339      std::max(max_volume_order,
 
 1341  max_volume_order = std::max(
 
 1345  for (
int gg = 0; gg != nb_gauss_pts; gg++) {
 
 1346    for (
int oo = 0; oo < max_volume_order; oo++) {
 
 1352          for (
auto ff = 0; ff != 4; ++ff) {
 
 1353            for (
int ee = 0; ee != 3; ++ee) {
 
 1354              t_base(
i) = t_base_f_e[ff * 3 + ee](
i);
 
 1356              ++t_base_f_e[ff * 3 + ee];
 
 1360            for (
int ee = 0; ee != 3; ++ee) {
 
 1361              t_diff_base(
i, 
j) = t_diff_base_f_e[ff * 3 + ee](
i, 
j);
 
 1363              ++t_diff_base_f_e[ff * 3 + ee];
 
 1372          for (
auto ff = 0; ff != 4; ++ff) {
 
 1373            t_base(
i) = t_base_f_f[ff](
i);
 
 1376            t_diff_base(
i, 
j) = t_diff_base_f_f[ff](
i, 
j);
 
 1378            ++t_diff_base_f_f[ff];
 
 1388          for (
int ee = 0; ee < 6; ++ee) {
 
 1389            t_base(
i) = t_base_v_e[ee](
i);
 
 1392            t_diff_base(
i, 
j) = t_diff_base_v_e[ee](
i, 
j);
 
 1394            ++t_diff_base_v_e[ee];
 
 1396            ++count_dofs_volume;
 
 1404          for (
int ff = 0; ff < 4; ff++) {
 
 1405            t_base(
i) = t_base_v_f[ff](
i);
 
 1408            t_diff_base(
i, 
j) = t_diff_base_v_f[ff](
i, 
j);
 
 1410            ++t_diff_base_v_f[ff];
 
 1412            ++count_dofs_volume;
 
 1421          t_base(
i) = t_base_v(
i);
 
 1424          t_diff_base(
i, 
j) = t_diff_base_v(
i, 
j);
 
 1428          ++count_dofs_volume;
 
 1434  if (nb_dofs != count_dofs / nb_gauss_pts) {
 
 1436    MOFEM_LOG(
"SELF", Sev::error) << 
"Nb dofs face: " << 4 * nb_dofs_face
 
 1437                                  << 
" -> " << count_dofs_face / nb_gauss_pts;
 
 1438    MOFEM_LOG(
"SELF", Sev::error) << 
"Nb dofs volume: " << nb_dofs_volume
 
 1439                                  << 
" -> " << count_dofs_volume / nb_gauss_pts;
 
 1441             "Number of dofs %d is different than expected %d",
 
 1442             count_dofs / nb_gauss_pts, nb_dofs);
 
 1446  if (interior_cache_ptr) {
 
 1447    auto p = interior_cache_ptr->emplace(
 
 
 1463             "This should be used only with DEMKOWICZ_JACOBI_BASE " 
 1467  int nb_gauss_pts = pts.size2();
 
 1473  double *diff_phi_f[4];
 
 1485  auto face_cache_ptr = get_face_cache_ptr();
 
 1488  for (
int ff = 0; ff != 4; ff++) {
 
 1490    int order = volume_order > face_order ? volume_order : face_order;
 
 1498    if (face_cache_ptr) {
 
 1499      auto it = face_cache_ptr->find(boost::make_tuple(
 
 1501          face_order, nb_gauss_pts,
 
 1506      if (it != face_cache_ptr->end()) {
 
 1508        noalias(data.
dataOnEntities[MBTRI][ff].getDiffN(base)) = it->diffN;
 
 1514    phi_f[ff] = &*data.
dataOnEntities[MBTRI][ff].getN(base).data().begin();
 
 1521        &*data.
dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
 
 1522        phi_f[ff], diff_phi_f[ff], nb_gauss_pts, 4);
 
 1523    if (face_cache_ptr) {
 
 1525          face_order, nb_gauss_pts, data.
facesNodes(ff, 0),
 
 1543  auto interior_cache_ptr = get_interior_cache();
 
 1552    for (
int v = 0; 
v != 1; ++
v) {
 
 1553      if (interior_cache_ptr) {
 
 1554        auto it = interior_cache_ptr->find(
 
 1555            boost::make_tuple(volume_order, nb_gauss_pts));
 
 1556        if (it != interior_cache_ptr->end()) {
 
 1558          noalias(data.
dataOnEntities[MBTET][0].getDiffN(base)) = it->diffN;
 
 1563      double *phi_v = &*data.
dataOnEntities[MBTET][0].getN(base).data().begin();
 
 1564      double *diff_phi_v =
 
 1568          volume_order, &data.
dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
 
 1569          &data.
dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), p_f, phi_f,
 
 1570          diff_phi_f, phi_v, diff_phi_v, nb_gauss_pts);
 
 1571      if (interior_cache_ptr) {
 
 1572        auto p = interior_cache_ptr->emplace(
 
 1581  for (
int ff = 0; ff != 4; ff++) {
 
 
 1600             "This should be used only with DEMKOWICZ_JACOBI_BASE " 
 1604  int nb_gauss_pts = pts.size2();
 
 1609  int nb_dofs = 4 * nb_dofs_face + nb_dofs_volume;
 
 1610  data.
dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
 
 1612  data.
dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 9 * nb_dofs,
 
 1630  auto interior_cache_ptr = get_interior_cache();
 
 1632  if (interior_cache_ptr) {
 
 1634        interior_cache_ptr->find(boost::make_tuple(volume_order, nb_gauss_pts));
 
 1635    if (it != interior_cache_ptr->end()) {
 
 1637      noalias(data.
dataOnEntities[MBTET][0].getDiffN(base)) = it->diffN;
 
 1642  std::array<MatrixDouble, 4> face_base_fun{
 
 1647  std::array<MatrixDouble, 4> face_diff_base{
 
 1653  int faces_nodes[4][3] = {{0, 1, 3}, {1, 2, 3}, {0, 3, 2}, {0, 2, 1}};
 
 1655  std::array<int, 4> p_f{volume_order, volume_order, volume_order,
 
 1657  std::array<double *, 4> phi_f{
 
 1658      &*face_base_fun[0].data().begin(), &*face_base_fun[1].data().begin(),
 
 1659      &*face_base_fun[2].data().begin(), &*face_base_fun[3].data().begin()};
 
 1660  std::array<double *, 4> diff_phi_f{
 
 1661      &*face_diff_base[0].data().begin(), &*face_diff_base[1].data().begin(),
 
 1662      &*face_diff_base[2].data().begin(), &*face_diff_base[3].data().begin()};
 
 1665  for (
int ff = 0; ff != 4; ff++) {
 
 1668        faces_nodes[ff], p_f[ff],
 
 1670        &*data.
dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
 
 1671        phi_f[ff], diff_phi_f[ff], nb_gauss_pts, 4);
 
 1674  MatrixDouble vol_bases(nb_gauss_pts, 3 * nb_dofs_volume);
 
 1675  MatrixDouble vol_diff_bases(nb_gauss_pts, 9 * nb_dofs_volume);
 
 1676  auto *phi_v = &*vol_bases.data().begin();
 
 1677  auto *diff_phi_v = &*vol_diff_bases.data().begin();
 
 1679      volume_order, &data.
dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
 
 1680      &data.
dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), p_f.data(),
 
 1681      phi_f.data(), diff_phi_f.data(), phi_v, diff_phi_v, nb_gauss_pts);
 
 1685      getFTensor1FromPtr<3>(phi_f[0]), getFTensor1FromPtr<3>(phi_f[1]),
 
 1686      getFTensor1FromPtr<3>(phi_f[2]), getFTensor1FromPtr<3>(phi_f[3])};
 
 1695      getFTensor1FromPtr<3>(&*vol_bases.data().begin());
 
 1699  auto t_base = getFTensor1FromPtr<3>(
 
 1707  for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
 
 1708    for (
int oo = 0; oo < volume_order; oo++) {
 
 1712        for (
auto ff = 0; ff != 4; ++ff) {
 
 1713          t_base(
i) = t_base_v_f[ff](
i);
 
 1716          t_diff_base(
i, 
j) = t_diff_base_v_f[ff](
i, 
j);
 
 1718          ++t_diff_base_v_f[ff];
 
 1724        t_base(
i) = t_base_v(
i);
 
 1727        t_diff_base(
i, 
j) = t_diff_base_v(
i, 
j);
 
 1734  if (interior_cache_ptr) {
 
 1735    auto p = interior_cache_ptr->emplace(
 
 
 1784  PetscErrorCode (*base_polynomials)(
int p, 
double s, 
double *diff_s, 
double *
L,
 
 1785                                     double *diffL, 
const int dim) =
 
 1788  int nb_gauss_pts = pts.size2();
 
 1792    int sense[6], 
order[6];
 
 1796    double *hcurl_edge_n[6], *diff_hcurl_edge_n[6];
 
 1797    for (
int ee = 0; ee != 6; ee++) {
 
 1800                "data inconsistency");
 
 1807                                                        3 * nb_dofs, 
false);
 
 1808      data.
dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
 
 1809                                                            9 * nb_dofs, 
false);
 
 1812      diff_hcurl_edge_n[ee] =
 
 1818        &*data.
dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
 
 1819        hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts, base_polynomials);
 
 1821    for (
int ee = 0; ee != 6; ee++) {
 
 1822      data.
dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, 
false);
 
 1823      data.
dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
 
 1835    double *hcurl_base_n[4], *diff_hcurl_base_n[4];
 
 1836    for (
int ff = 0; ff != 4; ff++) {
 
 1839                "data inconsistency");
 
 1844                                                       3 * nb_dofs, 
false);
 
 1845      data.
dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
 
 1846                                                           9 * nb_dofs, 
false);
 
 1849      diff_hcurl_base_n[ff] =
 
 1861        &*data.
dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
 
 1862        hcurl_base_n, diff_hcurl_base_n, nb_gauss_pts, base_polynomials);
 
 1864    for (
int ff = 0; ff != 4; ff++) {
 
 1865      data.
dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, 0, 
false);
 
 1866      data.
dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts, 0,
 
 1877                                                    3 * nb_vol_dofs, 
false);
 
 1878    data.
dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
 
 1879                                                        9 * nb_vol_dofs, 
false);
 
 1883        &*data.
dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
 
 1886        nb_gauss_pts, base_polynomials);
 
 1889    data.
dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 0, 
false);
 
 1890    data.
dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 0, 
false);
 
 
 1904             "This should be used only with DEMKOWICZ_JACOBI_BASE " 
 1909  int nb_gauss_pts = pts.size2();
 
 1913    int sense[6], 
order[6];
 
 1916               "wrong size of data structure, expected space for six edges " 
 1920    double *hcurl_edge_n[6], *diff_hcurl_edge_n[6];
 
 1921    for (
int ee = 0; ee != 6; ee++) {
 
 1924                "orintation of edges is not set");
 
 1931                                                        3 * nb_dofs, 
false);
 
 1932      data.
dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
 
 1933                                                            9 * nb_dofs, 
false);
 
 1936      diff_hcurl_edge_n[ee] =
 
 1942        &*data.
dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
 
 1943        hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts);
 
 1947    for (
int ee = 0; ee != 6; ee++) {
 
 1948      data.
dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, 
false);
 
 1949      data.
dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
 
 1960               "data structure for storing face h-curl base have wrong size " 
 1961               "should be four but is %ld",
 
 1964    double *hcurl_base_n[4], *diff_hcurl_base_n[4];
 
 1965    for (
int ff = 0; ff != 4; ff++) {
 
 1968                "orintation of face is not set");
 
 1973                                                       3 * nb_dofs, 
false);
 
 1974      data.
dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
 
 1975                                                           9 * nb_dofs, 
false);
 
 1978      diff_hcurl_base_n[ff] =
 
 1983              "data inconsistency, should be four faces");
 
 1987              "data inconsistency, should be three nodes on face");
 
 1992        &*data.
dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
 
 1993        hcurl_base_n, diff_hcurl_base_n, nb_gauss_pts);
 
 1997    for (
int ff = 0; ff != 4; ff++) {
 
 1998      data.
dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, 0, 
false);
 
 1999      data.
dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts, 0,
 
 2009                                                    3 * nb_vol_dofs, 
false);
 
 2010    data.
dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
 
 2011                                                        9 * nb_vol_dofs, 
false);
 
 2015        &*data.
dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
 
 2022    data.
dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 0, 
false);
 
 2023    data.
dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 0, 
false);
 
 
 2049                            boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
 
 2054  int nb_gauss_pts = pts.size2();
 
 2058  if (pts.size1() < 3)
 
 2061        "Wrong dimension of pts, should be at least 3 rows with coordinates");
 
 2067      data.
dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 4,
 
 2071          &pts(0, 0), &pts(1, 0), &pts(2, 0), nb_gauss_pts);
 
 2077        (
unsigned int)nb_gauss_pts) {
 
 2079               "Base functions or nodes has wrong number of integration points " 
 2083    data.
dataOnEntities[MBVERTEX][0].getDiffN(base).resize(4, 3, 
false);
 
 
 2137  switch (continuity) {
 
 2152            "Unknown (or not implemented) continuity");
 
 
 2167  auto set_ainsworth = [&dofs_side_map]() {
 
 2170    dofs_side_map.clear();
 
 2178        for (
auto ff = 0; ff != 4; ++ff) {
 
 2179          for (
int ee = 0; ee != 3; ++ee) {
 
 2189        for (
auto ff = 0; ff != 4; ++ff) {
 
 2198        for (
int ee = 0; ee < 6; ee++) {
 
 2206        for (
int ff = 0; ff < 4; ff++) {
 
 2222  auto set_demkowicz = [&dofs_side_map]() {
 
 2225    dofs_side_map.clear();
 
 2233        for (
auto ff = 0; ff != 4; ++ff) {
 
 2249  switch (continuity) {
 
 2264            "Unknown (or not implemented) continuity");
 
 
 2270template <
typename T>
 
 2272  auto it = cache.find(ptr);
 
 2273  if (it != cache.end()) {
 
 2276        << 
"Cache off " << cache_name << 
": " << it->second.size();
 
 2282        << 
"Cache on " << cache_name;
 
 
 2292                        std::string(
"hDivBaseFace") +
 
 
 2297bool TetPolynomialBase::switchCacheBaseInterior<HDIV>(
 
 2300                        std::string(
"hdivBaseInterior") +
 
 
 2305bool TetPolynomialBase::switchCacheBrokenBaseInterior<HDIV>(
 
 2308                        std::string(
"hdivBrokenBaseInterior") +
 
 
 2314                                                std::vector<void *> 
v) {
 
 2315  for (
auto fe_ptr : 
v) {
 
 2316    if (!TetPolynomialBase::switchCacheBaseFace<HDIV>(base, fe_ptr)) {
 
 2317      TetPolynomialBase::switchCacheBaseFace<HDIV>(base, fe_ptr);
 
 2319    if (!TetPolynomialBase::switchCacheBaseInterior<HDIV>(base, fe_ptr)) {
 
 2320      TetPolynomialBase::switchCacheBaseInterior<HDIV>(base, fe_ptr);
 
 2322    if (!TetPolynomialBase::switchCacheBrokenBaseInterior<HDIV>(base, fe_ptr)) {
 
 2323      TetPolynomialBase::switchCacheBrokenBaseInterior<HDIV>(base, fe_ptr);
 
 
 2330                                                 std::vector<void *> 
v) {
 
 2331  for (
auto fe_ptr : 
v) {
 
 2332    if (TetPolynomialBase::switchCacheBaseFace<HDIV>(base, fe_ptr)) {
 
 2333      TetPolynomialBase::switchCacheBaseFace<HDIV>(base, fe_ptr);
 
 2335    if (TetPolynomialBase::switchCacheBaseInterior<HDIV>(base, fe_ptr)) {
 
 2336      TetPolynomialBase::switchCacheBaseInterior<HDIV>(base, fe_ptr);
 
 2338    if (TetPolynomialBase::switchCacheBrokenBaseInterior<HDIV>(base, fe_ptr)) {
 
 2339      TetPolynomialBase::switchCacheBrokenBaseInterior<HDIV>(base, fe_ptr);
 
 
 2345void TetPolynomialBase::switchCacheBaseOn<HDIV>(std::vector<void *> 
v) {
 
 2346  for (
auto b = 0; b != 
LASTBASE; ++b) {
 
 
 2352void TetPolynomialBase::switchCacheBaseOff<HDIV>(std::vector<void *> 
v) {
 
 2353  for (
auto b = 0; b != 
LASTBASE; ++b) {
 
 
 2362                        std::string(
"hdivBaseInterior") +
 
 
 2368                                              std::vector<void *> 
v) {
 
 2369  for (
auto fe_ptr : 
v) {
 
 2370    if (!TetPolynomialBase::switchCacheBaseInterior<L2>(base, fe_ptr)) {
 
 2371      TetPolynomialBase::switchCacheBaseInterior<L2>(base, fe_ptr);
 
 
 2378                                               std::vector<void *> 
v) {
 
 2379  for (
auto fe_ptr : 
v) {
 
 2380    if (TetPolynomialBase::switchCacheBaseInterior<L2>(base, fe_ptr)) {
 
 2381      TetPolynomialBase::switchCacheBaseInterior<L2>(base, fe_ptr);
 
 
 2387void TetPolynomialBase::switchCacheBaseOn<L2>(std::vector<void *> 
v) {
 
 2388  for (
auto b = 0; b != 
LASTBASE; ++b) {
 
 
 2394void TetPolynomialBase::switchCacheBaseOff<L2>(std::vector<void *> 
v) {
 
 2395  for (
auto b = 0; b != 
LASTBASE; ++b) {
 
 
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
#define FTENSOR_INDEX(DIM, I)
auto tetCacheSwitch(const void *ptr, T &cache, std::string cache_name)
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
FieldSpace
approximation spaces
@ L2
field with C-1 continuity
@ HCURL
field with continuous tangents
@ HDIV
field with continuous normal traction
FieldContinuity
Field continuity.
@ CONTINUOUS
Regular field.
@ DISCONTINUOUS
Broken continuity (No effect on L2 space)
static const char *const FieldSpaceNames[]
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
static const char *const ApproximationBaseNames[]
#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 MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
#define NBVOLUMETET_AINSWORTH_EDGE_HDIV(P)
#define NBEDGE_DEMKOWICZ_HCURL(P)
#define NBVOLUMETET_H1(P)
Number of base functions on tetrahedron for H1 space.
#define NBVOLUMETET_AINSWORTH_HCURL(P)
#define NBFACETRI_AINSWORTH_HCURL(P)
#define NBVOLUMETET_DEMKOWICZ_HDIV(P)
PetscErrorCode H1_EdgeShapeFunctions_MBTET(int *sense, int *p, double *N, double *diffN, double *edgeN[], double *diff_edgeN[], int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
#define NBVOLUMETET_AINSWORTH_FACE_HDIV(P)
#define NBEDGE_H1(P)
Number of base function on edge for H1 space.
#define NBFACETRI_DEMKOWICZ_HDIV(P)
PetscErrorCode L2_Ainsworth_ShapeFunctions_MBTET(int p, double *N, double *diffN, double *L2N, double *diff_L2N, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Get base functions on tetrahedron for L2 space.
#define NBVOLUMETET_DEMKOWICZ_HCURL(P)
#define NBFACETRI_DEMKOWICZ_HCURL(P)
#define NBFACETRI_AINSWORTH_FACE_HDIV(P)
#define NBVOLUMETET_AINSWORTH_VOLUME_HDIV(P)
#define NBFACETRI_AINSWORTH_EDGE_HDIV(P)
PetscErrorCode H1_VolumeShapeFunctions_MBTET(int p, double *N, double *diffN, double *volumeN, double *diff_volumeN, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
#define NBEDGE_AINSWORTH_HCURL(P)
PetscErrorCode H1_FaceShapeFunctions_MBTET(int *faces_nodes, int *p, double *N, double *diffN, double *faceN[], double *diff_faceN[], int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
#define NBVOLUMETET_L2(P)
Number of base functions on tetrahedron for L2 space.
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium
FTensor::Index< 'j', 3 > j
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
UBlasMatrix< int > MatrixInt
implementation of Data Operators for Forces and Sources
FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > getFTensor2HVecFromPtr< 3, 3 >(double *ptr)
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_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 Hcurl_Ainsworth_VolumeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v, double *diff_phi_v, int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
H-curl volume base functions.
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 Hcurl_Ainsworth_EdgeBaseFunctions_MBTET(int *sense, int *p, double *N, double *diffN, double *edgeN[], double *diff_edgeN[], int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Edge based H-curl base functions on tetrahedral.
MoFEMErrorCode Hcurl_Demkowicz_FaceBaseFunctions_MBTET(int *faces_nodes, int *p, double *n, double *diff_n, double *phi[], double *diff_phi[], int nb_integration_pts)
Face base interior function.
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 Hcurl_Demkowicz_VolumeBaseFunctions_MBTET(int p, double *n, double *diff_n, double *phi, double *diff_phi, int nb_integration_pts)
Volume base interior function.
MoFEMErrorCode Hcurl_Demkowicz_EdgeBaseFunctions_MBTET(int *sense, int *p, double *n, double *diff_n, double *phi[], double *diff_phi[], int nb_integration_pts)
Edge based H-curl base functions on tetrahedral.
MoFEMErrorCode Hcurl_Ainsworth_FaceFunctions_MBTET(int *face_nodes, int *p, double *N, double *diffN, double *phi_f[4], double *diff_phi_f[4], int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Face H-curl functions.
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 .
constexpr auto field_name
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
multi_index_container< DofsSideMapData, indexed_by< ordered_non_unique< tag< TypeSide_mi_tag >, composite_key< DofsSideMapData, member< DofsSideMapData, EntityType, &DofsSideMapData::type >, member< DofsSideMapData, int, &DofsSideMapData::side > > >, ordered_unique< tag< EntDofIdx_mi_tag >, member< DofsSideMapData, int, &DofsSideMapData::dof > > > > DofsSideMap
Map entity stype and side to element/entity dof index.
static MoFEMErrorCode generateIndicesTriTet(const int N[], int *alpha[])
static MoFEMErrorCode baseFunctionsTet(const int N, const int gdim, const int n_alpha, const int *alpha, const double *lambda, const double *grad_lambda, double *base, double *grad_base)
static MoFEMErrorCode generateIndicesEdgeTet(const int N[], int *alpha[])
static MoFEMErrorCode generateIndicesVertexTet(const int N, int *alpha)
static MoFEMErrorCode generateIndicesTetTet(const int N, int *alpha)
Class used to pass element data to calculate base functions on tet,triangle,edge.
const FieldApproximationBase copyNodeBase
const std::string fieldName
PetscErrorCode(* basePolynomialsType0)(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
const FieldContinuity spaceContinuity
const FieldApproximationBase bAse
data structure for finite element entity
MatrixInt facesNodes
nodes on finite element faces
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
static constexpr int maxBrokenDofsOrder
Maximum order for broken space DOFs.
Calculate base functions on tetrahedral.
ublas::matrix< MatrixDouble > diffN_face_edge
EntPolynomialBaseCtx * cTx
MatrixDouble N_volume_bubble
MoFEMErrorCode getValueHdivAinsworthBaseImpl(MatrixDouble &pts, MatrixDouble &shape_functions, MatrixDouble &diff_shape_functions, int volume_order, std::array< int, 4 > &faces_order, std::array< int, 3 *4 > &faces_nodes, boost::function< int(int)> broken_nbfacetri_edge_hdiv, boost::function< int(int)> broken_nbfacetri_face_hdiv, boost::function< int(int)> broken_nbvolumetet_edge_hdiv, boost::function< int(int)> broken_nbvolumetet_face_hdiv, boost::function< int(int)> broken_nbvolumetet_volume_hdiv)
MoFEMErrorCode getValueH1AinsworthBase(MatrixDouble &pts)
MoFEMErrorCode getValueHdivAinsworthBrokenBase(MatrixDouble &pts)
ublas::vector< MatrixDouble > diffN_volume_face
MoFEMErrorCode getValueL2BernsteinBezierBase(MatrixDouble &pts)
MoFEMErrorCode getValueHcurlDemkowiczBase(MatrixDouble &pts)
MoFEMErrorCode getValueHcurl(MatrixDouble &pts)
Get base functions for Hcurl space.
virtual ~TetPolynomialBase()
MoFEMErrorCode getValueH1BernsteinBezierBase(MatrixDouble &pts)
static MoFEMErrorCode setDofsSideMap(const FieldSpace space, const FieldContinuity continuity, const FieldApproximationBase base, DofsSideMap &)
Set map of dof to side number.
MoFEMErrorCode getValueH1(MatrixDouble &pts)
Get base functions for H1 space.
ublas::vector< MatrixDouble > diffN_face_bubble
ublas::vector< MatrixDouble > N_volume_edge
ublas::vector< MatrixDouble > N_volume_face
MoFEMErrorCode getValueHdiv(MatrixDouble &pts)
Get base functions for Hdiv space.
MoFEMErrorCode getValueHdivAinsworthBase(MatrixDouble &pts)
MoFEMErrorCode getValueL2AinsworthBase(MatrixDouble &pts)
ublas::vector< MatrixDouble > N_face_bubble
ublas::matrix< MatrixDouble > N_face_edge
MoFEMErrorCode getValueL2(MatrixDouble &pts)
Get base functions for L2 space.
MatrixDouble diffN_volume_bubble
MoFEMErrorCode getValueHdivDemkowiczBase(MatrixDouble &pts)
TetPolynomialBase(const void *ptr=nullptr)
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
MoFEMErrorCode getValueHdivDemkowiczBrokenBase(MatrixDouble &pts)
ublas::vector< MatrixDouble > diffN_volume_edge
MoFEMErrorCode getValueHcurlAinsworthBase(MatrixDouble &pts)
static MoFEMErrorCode setDofsSideMapHdiv(const FieldContinuity continuity, const FieldApproximationBase base, DofsSideMap &dofs_side_map)
Set the Dofs Side Map Hdiv object.
base class for all interface classes
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
boost::multi_index_container< BaseCacheItem, boost::multi_index::indexed_by< boost::multi_index::hashed_unique< composite_key< BaseCacheItem, member< BaseCacheItem, int, &BaseCacheItem::order >, member< BaseCacheItem, int, &BaseCacheItem::nb_gauss_pts > > > > > BaseCacheMI
static std::array< std::map< const void *, BaseCacheMI >, LASTBASE > hdivBaseInterior
boost::multi_index_container< HDivBaseCacheItem, boost::multi_index::indexed_by< boost::multi_index::hashed_unique< composite_key< HDivBaseCacheItem, member< HDivBaseCacheItem, int, &HDivBaseCacheItem::order >, member< HDivBaseCacheItem, int, &HDivBaseCacheItem::nb_gauss_pts >, member< HDivBaseCacheItem, int, &HDivBaseCacheItem::n0 >, member< HDivBaseCacheItem, int, &HDivBaseCacheItem::n1 >, member< HDivBaseCacheItem, int, &HDivBaseCacheItem::n2 > > > > > HDivBaseFaceCacheMI
static std::array< std::map< const void *, HDivBaseFaceCacheMI >, LASTBASE > hDivBaseFace
static std::array< std::map< const void *, BaseCacheMI >, LASTBASE > l2BaseInterior
static std::array< std::map< const void *, BaseCacheMI >, LASTBASE > hdivBrokenBaseInterior