CGGTonsorialBubbleBase.cpp

Go to the documentation of this file.

/** \file CGGTonsorialBubbleBase.hpp
 
  \brief Implementation of tonsorial bubble base div(v) = 0.
 
  Implementation is based and motiveted by \cite cockburn2010new. This base
  is used to approximate stresses using Hdiv base with weakly enforced
  symmetry.
 
*/
 
#include <MoFEM.hpp>
#include <h1_hdiv_hcurl_l2.h>
using namespace MoFEM;
#include <CGGTonsorialBubbleBase.hpp>
 
using namespace FTensor;
 
namespace EshelbianPlasticity {
 
template <typename T>
MoFEMErrorCode CGG_BubbleBase_MBTET_Impl(const int p, const T *N,
                                         const T *diffN,
                                         Tensor2<PackPtr<T *, 9>, 3, 3> &t_phi,
                                         const int gdim) {
  MoFEMFunctionBeginHot;
 
  // FIXME: In this implementation symmetry of mix derivatives is not exploited
 
  if (p < 1)
    MoFEMFunctionReturnHot(0);
 
  Index<'i', 3> i;
  Index<'j', 3> j;
  Index<'k', 3> k;
  Index<'m', 3> m;
  Index<'f', 3> f;
 
  Tensor1<T, 3> t_diff_n[4];
  {
    Tensor1<PackPtr<const T *, 3>, 3> t_diff_n_tmp(&diffN[0], &diffN[1],
                                                   &diffN[2]);
    for (int ii = 0; ii != 4; ++ii) {
      t_diff_n[ii](i) = t_diff_n_tmp(i);
      ++t_diff_n_tmp;
    }
  }
 
  constexpr AinsworthHdivOptions::EdgeCoordinate ccg_ksi_type =
      AinsworthHdivOptions::edge_coordinate;
 
  auto get_ksi = [](const T n0, const T n1) {
    if constexpr (ccg_ksi_type == AinsworthHdivOptions::EdgeCoordinate::PAPER) {
      return n1 - n0;
    } else if constexpr (ccg_ksi_type ==
                         AinsworthHdivOptions::EdgeCoordinate::ENDPOINT_O) {
      return T(1) - T(2) * n0;
 
    } else if constexpr (ccg_ksi_type ==
                         AinsworthHdivOptions::EdgeCoordinate::ENDPOINT_I) {
      return T(2) * n1 - T(1);
    }
    return T(0);
  };
 
  auto get_diff_ksi = [&](const int n0, const int n1) {
    Tensor1<T, 3> t_diff_ksi;
    for (int dd = 0; dd != 3; ++dd) {
      if constexpr (ccg_ksi_type ==
                    AinsworthHdivOptions::EdgeCoordinate::PAPER) {
        t_diff_ksi(dd) = t_diff_n[n1](dd) - t_diff_n[n0](dd);
      } else if constexpr (ccg_ksi_type ==
                           AinsworthHdivOptions::EdgeCoordinate::ENDPOINT_O) {
        t_diff_ksi(dd) = -T(2) * t_diff_n[n0](dd);
      } else if constexpr (ccg_ksi_type ==
                           AinsworthHdivOptions::EdgeCoordinate::ENDPOINT_I) {
        t_diff_ksi(dd) = T(2) * t_diff_n[n1](dd);
      }
    }
    return t_diff_ksi;
  };
 
  Tensor1<T, 3> t_diff_ksi[3];
  for (int ii = 0; ii != 3; ++ii)
    t_diff_ksi[ii] = get_diff_ksi(0, ii + 1);
 
  int lp = p >= 2 ? p - 2 + 1 : 0;
  UBlasVector<T> l[3] = {VectorDouble(lp + 1), VectorDouble(lp + 1),
                         VectorDouble(lp + 1)};
  UBlasMatrix<T> diff_l[3] = {MatrixDouble(3, lp + 1), MatrixDouble(3, lp + 1),
                              MatrixDouble(3, lp + 1)};
  UBlasMatrix<T> diff2_l[3] = {MatrixDouble(9, lp + 1), MatrixDouble(9, lp + 1),
                               MatrixDouble(9, lp + 1)};
 
  for (int ii = 0; ii != 3; ++ii)
    diff2_l[ii].clear();
 
  for (int gg = 0; gg != gdim; ++gg) {
 
    const int node_shift = gg * 4;
 
    for (int ii = 0; ii != 3; ++ii) {
 
      auto &t_diff_ksi_ii = t_diff_ksi[ii];
      auto &l_ii = l[ii];
      auto &diff_l_ii = diff_l[ii];
      auto &diff2_l_ii = diff2_l[ii];
      const T ksi_ii = get_ksi(N[node_shift], N[node_shift + ii + 1]);
      CHKERR Legendre_polynomials(lp, ksi_ii, &t_diff_ksi_ii(0),
                                  &*l_ii.data().begin(),
                                  &*diff_l_ii.data().begin(), 3);
 
      for (int l = 1; l < lp; ++l) {
        const T a = ((2 * (T)l + 1) / ((T)l + 1));
        const T b = ((T)l / ((T)l + 1));
        for (int d0 = 0; d0 != 3; ++d0)
          for (int d1 = 0; d1 != 3; ++d1) {
            const int r = 3 * d0 + d1;
            diff2_l_ii(r, l + 1) = a * (t_diff_ksi_ii(d0) * diff_l_ii(d1, l) +
                                        t_diff_ksi_ii(d1) * diff_l_ii(d0, l) +
                                        ksi_ii * diff2_l_ii(r, l)) -
                                   b * diff2_l_ii(r, l - 1);
          }
      }
    }
 
    const T n[] = {N[node_shift + 0], N[node_shift + 1], N[node_shift + 2],
                   N[node_shift + 3]};
 
    Tensor2<T, 3, 3> t_bk;
    Tensor3<T, 3, 3, 3> t_bk_diff;
    const int tab[4][4] = {
        {1, 2, 3, 0}, {2, 3, 0, 1}, {3, 0, 1, 2}, {0, 1, 2, 3}};
    t_bk(i, j) = 0;
    t_bk_diff(i, j, k) = 0;
    for (int ii = 0; ii != 4; ++ii) {
      const int i0 = tab[ii][0];
      const int i1 = tab[ii][1];
      const int i2 = tab[ii][2];
      const int i3 = tab[ii][3];
      auto &t_diff_n_i0 = t_diff_n[i0];
      auto &t_diff_n_i1 = t_diff_n[i1];
      auto &t_diff_n_i2 = t_diff_n[i2];
      auto &t_diff_n_i3 = t_diff_n[i3];
      Tensor2<T, 3, 3> t_k;
      t_k(i, j) = t_diff_n_i3(i) * t_diff_n_i3(j);
      const T b = n[i0] * n[i1] * n[i2];
      t_bk(i, j) += b * t_k(i, j);
      Tensor1<T, 3> t_diff_b;
      t_diff_b(i) = t_diff_n_i0(i) * n[i1] * n[i2] +
                    t_diff_n_i1(i) * n[i0] * n[i2] +
                    t_diff_n_i2(i) * n[i0] * n[i1];
      t_bk_diff(i, j, k) += t_k(i, j) * t_diff_b(k);
    }
 
    int zz = 0;
    for (int o = p - 2 + 1; o <= p - 2 + 1; ++o) {
 
      for (int ii = 0; ii <= o; ++ii)
        for (int jj = 0; (ii + jj) <= o; ++jj) {
 
          const int kk = o - ii - jj;
 
          auto get_diff_l = [&](const int y, const int i) {
            return Tensor1<T, 3>(diff_l[y](0, i), diff_l[y](1, i),
                                 diff_l[y](2, i));
          };
          auto get_diff2_l = [&](const int y, const int i) {
            return Tensor2<T, 3, 3>(
                diff2_l[y](0, i), diff2_l[y](1, i), diff2_l[y](2, i),
                diff2_l[y](3, i), diff2_l[y](4, i), diff2_l[y](5, i),
                diff2_l[y](6, i), diff2_l[y](7, i), diff2_l[y](8, i));
          };
 
          auto l_i = l[0][ii];
          auto t_diff_i = get_diff_l(0, ii);
          auto t_diff2_i = get_diff2_l(0, ii);
          auto l_j = l[1][jj];
          auto t_diff_j = get_diff_l(1, jj);
          auto t_diff2_j = get_diff2_l(1, jj);
          auto l_k = l[2][kk];
          auto t_diff_k = get_diff_l(2, kk);
          auto t_diff2_k = get_diff2_l(2, kk);
 
          Tensor1<T, 3> t_diff_l2;
          t_diff_l2(i) = t_diff_i(i) * l_j * l_k + t_diff_j(i) * l_i * l_k +
                         t_diff_k(i) * l_i * l_j;
          Tensor2<T, 3, 3> t_diff2_l2;
          t_diff2_l2(i, j) =
              t_diff2_i(i, j) * l_j * l_k + t_diff_i(i) * t_diff_j(j) * l_k +
              t_diff_i(i) * l_j * t_diff_k(j) +
 
              t_diff2_j(i, j) * l_i * l_k + t_diff_j(i) * t_diff_i(j) * l_k +
              t_diff_j(i) * l_i * t_diff_k(j) +
 
              t_diff2_k(i, j) * l_i * l_j + t_diff_k(i) * t_diff_i(j) * l_j +
              t_diff_k(i) * l_i * t_diff_j(j);
 
          for (int dd = 0; dd != 3; ++dd) {
 
            Tensor2<T, 3, 3> t_axial_diff;
            t_axial_diff(i, j) = 0;
            for (int mm = 0; mm != 3; ++mm)
              t_axial_diff(dd, mm) = t_diff_l2(mm);
 
            FTensor::Tensor3<T, 3, 3, 3> t_levi_civita;
            t_levi_civita(i, j, k) = FTensor::levi_civita(i, j, k);
 
            Tensor3<T, 3, 3, 3> t_A_diff;
            t_A_diff(i, j, k) =
                levi_civita<double>(i, j, m) * t_axial_diff(m, k);
            Tensor2<T, 3, 3> t_curl_A;
            t_curl_A(i, j) = levi_civita<double>(j, m, f) * t_A_diff(i, f, m);
            Tensor3<T, 3, 3, 3> t_curl_A_bK_diff;
            t_curl_A_bK_diff(i, j, k) = t_curl_A(i, m) * t_bk_diff(m, j, k);
 
            Tensor3<T, 3, 3, 3> t_axial_diff2;
            t_axial_diff2(i, j, k) = 0;
            for (int mm = 0; mm != 3; ++mm)
              for (int nn = 0; nn != 3; ++nn)
                t_axial_diff2(dd, mm, nn) = t_diff2_l2(mm, nn);
            Tensor4<T, 3, 3, 3, 3> t_A_diff2;
            t_A_diff2(i, j, k, f) =
                levi_civita<double>(i, j, m) * t_axial_diff2(m, k, f);
            Tensor3<T, 3, 3, 3> t_curl_A_diff2;
            t_curl_A_diff2(i, j, k) =
                levi_civita<double>(j, m, f) * t_A_diff2(i, f, m, k);
            Tensor3<T, 3, 3, 3> t_curl_A_diff2_bK;
            t_curl_A_diff2_bK(i, j, k) = t_curl_A_diff2(i, m, k) * t_bk(m, j);
 
            t_phi(i, j) =
                levi_civita<double>(j, m, f) *
                (t_curl_A_bK_diff(i, f, m) + t_curl_A_diff2_bK(i, f, m));
 
            ++t_phi;
            ++zz;
          }
        }
    }
#ifndef NDEBUG
    if (zz != NBVOLUMETET_CCG_BUBBLE(p))
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Wrong number of base functions %d != %d", zz,
              NBVOLUMETET_CCG_BUBBLE(p));
#endif
  }
 
  MoFEMFunctionReturnHot(0);
}
 
MoFEMErrorCode CGG_BubbleBase_MBTET(const int p, const double *N,
                                    const double *diffN,
                                    Tensor2<PackPtr<double *, 9>, 3, 3> &t_phi,
                                    const int gdim) {
  return CGG_BubbleBase_MBTET_Impl<double>(p, N, diffN, t_phi, gdim);
}
 
MoFEMErrorCode CGG_BubbleBase_MBTET(
    const int p, const std::complex<double> *N,
    const std::complex<double> *diffN,
    FTensor::Tensor2<FTensor::PackPtr<std::complex<double> *, 9>, 3, 3> &t_phi,
    const int gdim) {
  return CGG_BubbleBase_MBTET_Impl<std::complex<double>>(p, N, diffN, t_phi,
                                                         gdim);
}
 
} // namespace EshelbianPlasticity