v0.14.0
CGGTonsorialBubbleBase.cpp
Go to the documentation of this file.
1 /** \file CGGTonsorialBubbleBase.hpp
2 
3  \brief Implementation of tonsorial bubble base div(v) = 0.
4 
5  Implementation is based and motiveted by \cite cockburn2010new. This base
6  is used to approximate stresses using Hdiv base with weakly enforced
7  symmetry.
8 
9 */
10 
11 #include <MoFEM.hpp>
12 #include <h1_hdiv_hcurl_l2.h>
13 using namespace MoFEM;
15 
16 using namespace FTensor;
17 
18 namespace EshelbianPlasticity {
19 
20 MoFEMErrorCode CGG_BubbleBase_MBTET(const int p, const double *N,
21  const double *diffN,
22  Tensor2<PackPtr<double *, 9>, 3, 3> &t_phi,
23  const int gdim) {
25 
26  // FIXME: In this implementation symmetry of mix derivatives is not exploited
27 
28  if (p < 1)
30 
36 
37  Tensor1<double, 3> t_diff_n[4];
38  {
39  Tensor1<PackPtr<const double *, 3>, 3> t_diff_n_tmp(&diffN[0], &diffN[1],
40  &diffN[2]);
41  for (int ii = 0; ii != 4; ++ii) {
42  t_diff_n[ii](i) = t_diff_n_tmp(i);
43  ++t_diff_n_tmp;
44  }
45  }
46 
47  Tensor1<double, 3> t_diff_ksi[3];
48  for (int ii = 0; ii != 3; ++ii)
49  t_diff_ksi[ii](i) = t_diff_n[ii + 1](i) - t_diff_n[0](i);
50 
51  int lp = p >= 2 ? p - 2 + 1 : 0;
52  VectorDouble l[3] = {VectorDouble(lp + 1), VectorDouble(lp + 1),
53  VectorDouble(lp + 1)};
54  MatrixDouble diff_l[3] = {MatrixDouble(3, lp + 1), MatrixDouble(3, lp + 1),
55  MatrixDouble(3, lp + 1)};
56  MatrixDouble diff2_l[3] = {MatrixDouble(9, lp + 1), MatrixDouble(9, lp + 1),
57  MatrixDouble(9, lp + 1)};
58 
59  for (int ii = 0; ii != 3; ++ii)
60  diff2_l[ii].clear();
61 
62  for (int gg = 0; gg != gdim; ++gg) {
63 
64  const int node_shift = gg * 4;
65 
66  for (int ii = 0; ii != 3; ++ii) {
67 
68  auto &t_diff_ksi_ii = t_diff_ksi[ii];
69  auto &l_ii = l[ii];
70  auto &diff_l_ii = diff_l[ii];
71  auto &diff2_l_ii = diff2_l[ii];
72 
73  double ksi_ii = N[node_shift + ii + 1] - N[node_shift + 0];
74 
75  CHKERR Legendre_polynomials(lp, ksi_ii, &t_diff_ksi_ii(0),
76  &*l_ii.data().begin(),
77  &*diff_l_ii.data().begin(), 3);
78 
79  for (int l = 1; l < lp; ++l) {
80  const double a = ((2 * (double)l + 1) / ((double)l + 1));
81  const double b = ((double)l / ((double)l + 1));
82  for (int d0 = 0; d0 != 3; ++d0)
83  for (int d1 = 0; d1 != 3; ++d1) {
84  const int r = 3 * d0 + d1;
85  diff2_l_ii(r, l + 1) = a * (t_diff_ksi_ii(d0) * diff_l_ii(d1, l) +
86  t_diff_ksi_ii(d1) * diff_l_ii(d0, l) +
87  ksi_ii * diff2_l_ii(r, l)) -
88  b * diff2_l_ii(r, l - 1);
89  }
90  }
91  }
92 
93  const double n[] = {N[node_shift + 0], N[node_shift + 1], N[node_shift + 2],
94  N[node_shift + 3]};
95 
97  Tensor3<double, 3, 3, 3> t_bk_diff;
98  const int tab[4][4] = {
99  {1, 2, 3, 0}, {2, 3, 0, 1}, {3, 0, 1, 2}, {0, 1, 2, 3}};
100  t_bk(i, j) = 0;
101  t_bk_diff(i, j, k) = 0;
102  for (int ii = 0; ii != 3; ++ii) {
103  const int i0 = tab[ii][0];
104  const int i1 = tab[ii][1];
105  const int i2 = tab[ii][2];
106  const int i3 = tab[ii][3];
107  auto &t_diff_n_i0 = t_diff_n[i0];
108  auto &t_diff_n_i1 = t_diff_n[i1];
109  auto &t_diff_n_i2 = t_diff_n[i2];
110  auto &t_diff_n_i3 = t_diff_n[i3];
112  t_k(i, j) = t_diff_n_i3(i) * t_diff_n_i3(j);
113  const double b = n[i0] * n[i1] * n[i2];
114  t_bk(i, j) += b * t_k(i, j);
115  Tensor1<double, 3> t_diff_b;
116  t_diff_b(i) = t_diff_n_i0(i) * n[i1] * n[i2] +
117  t_diff_n_i1(i) * n[i0] * n[i2] +
118  t_diff_n_i2(i) * n[i0] * n[i1];
119  t_bk_diff(i, j, k) += t_k(i, j) * t_diff_b(k);
120  }
121 
122  int zz = 0;
123  for (int o = p - 2 + 1; o <= p - 2 + 1; ++o) {
124 
125  for (int ii = 0; ii <= o; ++ii)
126  for (int jj = 0; (ii + jj) <= o; ++jj) {
127 
128  const int kk = o - ii - jj;
129 
130  auto get_diff_l = [&](const int y, const int i) {
131  return Tensor1<double, 3>(diff_l[y](0, i), diff_l[y](1, i),
132  diff_l[y](2, i));
133  };
134  auto get_diff2_l = [&](const int y, const int i) {
135  return Tensor2<double, 3, 3>(
136  diff2_l[y](0, i), diff2_l[y](1, i), diff2_l[y](2, i),
137  diff2_l[y](3, i), diff2_l[y](4, i), diff2_l[y](5, i),
138  diff2_l[y](6, i), diff2_l[y](7, i), diff2_l[y](8, i));
139  };
140 
141  auto l_i = l[0][ii];
142  auto t_diff_i = get_diff_l(0, ii);
143  auto t_diff2_i = get_diff2_l(0, ii);
144  auto l_j = l[1][jj];
145  auto t_diff_j = get_diff_l(1, jj);
146  auto t_diff2_j = get_diff2_l(1, jj);
147  auto l_k = l[2][kk];
148  auto t_diff_k = get_diff_l(2, kk);
149  auto t_diff2_k = get_diff2_l(2, kk);
150 
151  Tensor1<double, 3> t_diff_l2;
152  t_diff_l2(i) = t_diff_i(i) * l_j * l_k + t_diff_j(i) * l_i * l_k +
153  t_diff_k(i) * l_i * l_j;
154  Tensor2<double, 3, 3> t_diff2_l2;
155  t_diff2_l2(i, j) =
156  t_diff2_i(i, j) * l_j * l_k + t_diff_i(i) * t_diff_j(j) * l_k +
157  t_diff_i(i) * l_j * t_diff_k(j) +
158 
159  t_diff2_j(i, j) * l_i * l_k + t_diff_j(i) * t_diff_i(j) * l_k +
160  t_diff_j(i) * l_i * t_diff_k(j) +
161 
162  t_diff2_k(i, j) * l_i * l_j + t_diff_k(i) * t_diff_i(j) * l_j +
163  t_diff_k(i) * l_i * t_diff_j(j);
164 
165  for (int dd = 0; dd != 3; ++dd) {
166 
167  Tensor2<double, 3, 3> t_axial_diff;
168  t_axial_diff(i, j) = 0;
169  for (int mm = 0; mm != 3; ++mm)
170  t_axial_diff(dd, mm) = t_diff_l2(mm);
171 
172  Tensor3<double, 3, 3, 3> t_A_diff;
173  t_A_diff(i, j, k) = levi_civita(i, j, m) * t_axial_diff(m, k);
174  Tensor2<double, 3, 3> t_curl_A;
175  t_curl_A(i, j) = levi_civita(j, m, f) * t_A_diff(i, f, m);
176  Tensor3<double, 3, 3, 3> t_curl_A_bK_diff;
177  t_curl_A_bK_diff(i, j, k) = t_curl_A(i, m) * t_bk_diff(m, j, k);
178 
179  Tensor3<double, 3, 3, 3> t_axial_diff2;
180  t_axial_diff2(i, j, k) = 0;
181  for (int mm = 0; mm != 3; ++mm)
182  for (int nn = 0; nn != 3; ++nn)
183  t_axial_diff2(dd, mm, nn) = t_diff2_l2(mm, nn);
184  Tensor4<double, 3, 3, 3, 3> t_A_diff2;
185  t_A_diff2(i, j, k, f) =
186  levi_civita(i, j, m) * t_axial_diff2(m, k, f);
187  Tensor3<double, 3, 3, 3> t_curl_A_diff2;
188  t_curl_A_diff2(i, j, k) =
189  levi_civita(j, m, f) * t_A_diff2(i, f, m, k);
190  Tensor3<double, 3, 3, 3> t_curl_A_diff2_bK;
191  t_curl_A_diff2_bK(i, j, k) = t_curl_A_diff2(i, m, k) * t_bk(m, j);
192 
193  t_phi(i, j) = levi_civita(j, m, f) * (t_curl_A_bK_diff(i, f, m) +
194  t_curl_A_diff2_bK(i, f, m));
195 
196  ++t_phi;
197  ++zz;
198  }
199  }
200  }
201  if (zz != NBVOLUMETET_CCG_BUBBLE(p))
202  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
203  "Wrong number of base functions %d != %d", zz,
205  }
206 
208 }
209 
210 } // namespace EshelbianPlasticity
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
FTensor
JSON compatible output.
Definition: Christof_constructor.hpp:6
FTensor::Tensor1< double, 3 >
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
EshelbianPlasticity
Definition: CGGTonsorialBubbleBase.hpp:11
MoFEM.hpp
FTensor::levi_civita
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
Definition: Levi_Civita.hpp:617
sdf.r
int r
Definition: sdf.py:8
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
FTensor::Tensor3
Definition: Tensor3_value.hpp:12
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
a
constexpr double a
Definition: approx_sphere.cpp:30
double
FTensor::PackPtr
Definition: FTensor.hpp:54
EshelbianPlasticity::CGG_BubbleBase_MBTET
MoFEMErrorCode CGG_BubbleBase_MBTET(const int p, const double *N, const double *diffN, FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > &phi, const int gdim)
Calculate CGGT tonsorial bubble base.
Definition: CGGTonsorialBubbleBase.cpp:20
Legendre_polynomials
PetscErrorCode Legendre_polynomials(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Calculate Legendre approximation basis.
Definition: base_functions.c:15
FTensor::Tensor4
Definition: Tensor4_value.hpp:18
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
FTensor::Index< 'i', 3 >
convert.n
n
Definition: convert.py:82
N
const int N
Definition: speed_test.cpp:3
FTensor::dd
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
CGGTonsorialBubbleBase.hpp
Implementation of tonsorial bubble base div(v) = 0.
NBVOLUMETET_CCG_BUBBLE
#define NBVOLUMETET_CCG_BUBBLE(P)
Bubble function for CGG H div space.
Definition: CGGTonsorialBubbleBase.hpp:19
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
h1_hdiv_hcurl_l2.h
Functions to approximate hierarchical spaces.
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21