v0.14.0
Hdiv.cpp
Go to the documentation of this file.
1 /** \file Hdiv.cpp
2 
3  \brief Implementation of H-curl base
4 
5  Based on Hierarchic Finite Element Bases on Unstructured Tetrahedral
6  Meshes, by Mark Ainsworth and Joe Coyle
7  Shape functions for MBTRI/MBTET and HCurl space
8 
9 */
10 
11 using namespace MoFEM;
12 
13 boost::function<int(int)> AinsworthOrderHooks::broken_nbfacetri_edge_hdiv =
14  [](int p) { return p; };
15 boost::function<int(int)> AinsworthOrderHooks::broken_nbfacetri_face_hdiv =
16  [](int p) { return p; };
17 boost::function<int(int)> AinsworthOrderHooks::broken_nbvolumetet_edge_hdiv =
18  [](int p) { return p; };
19 boost::function<int(int)> AinsworthOrderHooks::broken_nbvolumetet_face_hdiv =
20  [](int p) { return p; };
21 boost::function<int(int)> AinsworthOrderHooks::broken_nbvolumetet_volume_hdiv =
22  [](int p) { return p; };
23 
24 MoFEMErrorCode MoFEM::Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET(
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,
29  const int dim)) {
30 
31  MoFEMFunctionBeginHot;
32 #ifdef NDEBUG
33  if (!diff_phi_f_e)
34  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
35  "expected to return derivatives");
36 #endif
37 
38  for (int ff = 0; ff < 4; ff++) {
39  CHKERR Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET_ON_FACE(
40  &faces_nodes[3 * ff], p[ff], N, diffN, phi_f_e[ff], diff_phi_f_e[ff],
41  gdim, 4, base_polynomials);
42  }
43 
44  MoFEMFunctionReturnHot(0);
45 }
46 
47 MoFEMErrorCode MoFEM::Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET_ON_FACE(
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,
52  const int dim)) {
53 
54  constexpr int face_edges_nodes[3][2] = {{0, 1}, {1, 2}, {2, 0}};
55  constexpr int face_opposite_edges_node[] = {2, 0, 1};
56  FTENSOR_INDEX(3, i);
57  FTENSOR_INDEX(3, j);
58  FTENSOR_INDEX(3, k);
59 
60  MoFEMFunctionBeginHot;
61  if (p < 1)
62  MoFEMFunctionReturnHot(0);
63 
64  FTensor::Tensor1<double, 3> t_edge_cross[3];
65  FTensor::Tensor1<double, 3> t_node_diff_ksi[4];
66  FTensor::Tensor1<double, 3> t_diff_ksi[3];
67  if (diffN) {
68  t_node_diff_ksi[0] =
69  FTensor::Tensor1<double, 3>(diffN[0], diffN[1], diffN[2]);
70  t_node_diff_ksi[1] =
71  FTensor::Tensor1<double, 3>(diffN[3], diffN[4], diffN[5]);
72  t_node_diff_ksi[2] =
73  FTensor::Tensor1<double, 3>(diffN[6], diffN[7], diffN[8]);
74  t_node_diff_ksi[3] =
75  FTensor::Tensor1<double, 3>(diffN[9], diffN[10], diffN[11]);
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);
82  }
83  } else {
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;
88  }
89  }
90  double psi_l[p + 1];
91 
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);
103  auto t_psi_l = FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(psi_l);
104  int jj = 0;
105  for (int l = 0; l <= p - 1; l++) {
106  t_psi_f_e(i) = (lambda * t_psi_l) * t_edge_cross[ee](i);
107  ++t_psi_f_e;
108  ++t_psi_l;
109  ++jj;
110  }
111 #ifdef NDEBUG
112  if (jj != NBFACETRI_AINSWORTH_EDGE_HDIV(p)) {
113  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
114  "wrong order %d != %d", jj, NBFACETRI_AINSWORTH_EDGE_HDIV(p));
115  }
116 #endif
117  }
118  }
119 
120  if (diff_phi_f_e) {
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]];
126  auto t_diff_phi_f_e = getFTensor2HVecFromPtr<3, 3>(diff_phi_f_e[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);
133  auto t_psi_l = FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(psi_l);
134  auto t_diff_psi_l = FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 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)) *
139  t_edge_cross[ee](i);
140  ++t_diff_psi_l;
141  ++t_diff_phi_f_e;
142  ++t_psi_l;
143  }
144  }
145  }
146  }
147 
148  MoFEMFunctionReturnHot(0);
149 }
150 
151 MoFEMErrorCode MoFEM::Hdiv_Ainsworth_FaceBubbleShapeFunctions(
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,
156  const int dim)) {
157 
158  MoFEMFunctionBeginHot;
159 
160 #ifdef NDEBUG
161  if (!diff_phi_f)
162  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
163  "expected to return derivatives");
164 #endif
165 
166  for (int ff = 0; ff < 4; ff++) {
167  CHKERR Hdiv_Ainsworth_FaceBubbleShapeFunctions_ON_FACE(
168  &faces_nodes[3 * ff], p[ff], N, diffN, phi_f[ff], diff_phi_f[ff], gdim,
169  4, base_polynomials);
170  }
171  MoFEMFunctionReturnHot(0);
172 }
173 
174 MoFEMErrorCode MoFEM::Hdiv_Ainsworth_FaceBubbleShapeFunctions_ON_FACE(
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,
179  const int dim)) {
180  FTENSOR_INDEX(3, i);
181  FTENSOR_INDEX(3, j);
182  FTENSOR_INDEX(3, k);
183 
184  MoFEMFunctionBeginHot;
185  if (p < 3)
186  MoFEMFunctionReturnHot(0);
187 
188  const int vert_i = face_nodes[1];
189  const int vert_j = face_nodes[2];
190  const int i0 = face_nodes[0];
191  FTensor::Tensor1<double, 3> t_cross;
192  FTensor::Tensor1<double, 3> t_node_diff_ksi[4];
193  FTensor::Tensor1<double, 3> t_diff_ksi0i;
194  FTensor::Tensor1<double, 3> t_diff_ksi0j;
195 
196  if (diffN) {
197  t_node_diff_ksi[0] =
198  FTensor::Tensor1<double, 3>(diffN[0], diffN[1], diffN[2]);
199  t_node_diff_ksi[1] =
200  FTensor::Tensor1<double, 3>(diffN[3], diffN[4], diffN[5]);
201  t_node_diff_ksi[2] =
202  FTensor::Tensor1<double, 3>(diffN[6], diffN[7], diffN[8]);
203  t_node_diff_ksi[3] =
204  FTensor::Tensor1<double, 3>(diffN[9], diffN[10], diffN[11]);
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);
209  } else {
210  t_cross(0) = 1;
211  t_cross(1) = 0;
212  t_cross(2) = 0;
213  }
214 
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);
218 
219  for (int ii = 0; ii < gdim; ii++) {
220 
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);
230 
231  int jj = 0;
232  int oo = 0;
233  for (; oo <= p - 3; oo++) {
234  auto t_psi_l = FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(psi_l);
235  for (int l = 0; l <= oo; l++) {
236  int m = oo - l;
237  if (m >= 0) {
238  t_psi_f(i) = (beta_0ij * t_psi_l * psi_m[m]) * t_cross(i);
239  ++t_psi_f;
240  }
241  ++t_psi_l;
242  ++jj;
243  }
244  }
245 #ifdef NDEBUG
246  if (jj != NBFACETRI_AINSWORTH_FACE_HDIV(p)) {
247  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
248  "wrong order %d != %d", jj, NBFACETRI_AINSWORTH_FACE_HDIV(p));
249  }
250 #endif
251  }
252 
253  if (diff_phi_f) {
254  auto t_diff_phi_f = getFTensor2HVecFromPtr<3, 3>(diff_phi_f);
255 
256  for (int ii = 0; ii < gdim; ii++) {
257 
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;
265  FTensor::Tensor1<double, 3> t_diff_beta_0ij;
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);
271 
272  int jj = 0;
273  int oo = 0;
274  for (; oo <= p - 3; oo++) {
275  auto t_psi_l = FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(psi_l);
276  auto t_diff_psi_l = FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3>{
277  &diff_psi_l[0], &diff_psi_l[p + 1], &diff_psi_l[2 * p + 2]};
278  for (int l = 0; l <= oo; l++) {
279  int m = oo - l;
280  if (m >= 0) {
281  auto t_diff_psi_m =
282  FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3>{
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)) *
288  t_cross(i);
289  ++t_diff_phi_f;
290  }
291  ++t_psi_l;
292  ++t_diff_psi_l;
293  ++jj;
294  }
295  }
296 #ifndef NDEBUG
297  if (jj != NBFACETRI_AINSWORTH_FACE_HDIV(p)) {
298  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
299  "wrong order %d != %d", jj, NBFACETRI_AINSWORTH_FACE_HDIV(p));
300  }
301 #endif
302  }
303  }
304  MoFEMFunctionReturnHot(0);
305 }
306 
307 MoFEMErrorCode MoFEM::Hdiv_Ainsworth_EdgeBasedVolumeShapeFunctions_MBTET(
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,
312  const int dim)) {
313 
314  constexpr int edges_nodes[6][2] = {{0, 1}, {1, 2}, {2, 0},
315  {0, 3}, {1, 3}, {2, 3}};
316 
317  MoFEMFunctionBeginHot;
318  if (p < 2)
319  MoFEMFunctionReturnHot(0);
320  if (diffN == NULL) {
321  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
322  }
323 
324  FTensor::Tensor1<double, 3> t_coords[4] = {
325  FTensor::Tensor1<double, 3>(0., 0., 0.),
326  FTensor::Tensor1<double, 3>(1., 0., 0.),
327  FTensor::Tensor1<double, 3>(0., 1., 0.),
328  FTensor::Tensor1<double, 3>(0., 0., 1.)};
329  FTensor::Tensor1<double *, 3> t_node_diff_ksi[4] = {
330  FTensor::Tensor1<double *, 3>(&diffN[0], &diffN[1], &diffN[2]),
331  FTensor::Tensor1<double *, 3>(&diffN[3], &diffN[4], &diffN[5]),
332  FTensor::Tensor1<double *, 3>(&diffN[6], &diffN[7], &diffN[8]),
333  FTensor::Tensor1<double *, 3>(&diffN[9], &diffN[10], &diffN[11])};
334 
335  FTENSOR_INDEX(3, i);
336  FTENSOR_INDEX(3, j);
337 
338  FTensor::Tensor1<double, 3> t_tou_e;
339  FTensor::Tensor1<double, 3> t_diff_ksi0i;
340  FTensor::Tensor1<double, 3> t_diff_beta_e;
341 
342  double psi_l[p + 1];
343  double diff_psi_l[3 * (p + 1)];
344 
345  for (int ee = 0; ee != 6; ee++) {
346  t_tou_e(i) =
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);
358  auto t_psi_l = FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(psi_l);
359  for (int l = 0; l <= p - 2; l++) {
360  t_phi_v_e(i) = (beta_e * t_psi_l) * t_tou_e(i);
361  ++t_phi_v_e;
362  ++t_psi_l;
363  }
364  }
365  }
366 
367  if (diff_phi_v_e) {
368  for (int ee = 0; ee != 6; ee++) {
369  t_tou_e(i) =
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);
373  auto t_diff_phi_v_e = getFTensor2HVecFromPtr<3, 3>(diff_phi_v_e[ee]);
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);
383  auto t_psi_l = FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(psi_l);
384  auto t_diff_psi_l = FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 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)) *
389  t_tou_e(i);
390  ++t_diff_phi_v_e;
391  ++t_diff_psi_l;
392  ++t_psi_l;
393  }
394  }
395  }
396  }
397 
398  MoFEMFunctionReturnHot(0);
399 }
400 
401 MoFEMErrorCode MoFEM::Hdiv_Ainsworth_FaceBasedVolumeShapeFunctions_MBTET(
402  int p, double *N, double *diffN, double *phi_v_f[], double *diff_phi_v_f[],
403  int gdim,
404  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
405  double *L, double *diffL,
406  const int dim)) {
407 
408  constexpr int faces_nodes[4][3] = {
409  {0, 1, 3}, {1, 2, 3}, {0, 3, 2}, {0, 2, 1}};
410 
411  MoFEMFunctionBeginHot;
412  if (p < 3)
413  MoFEMFunctionReturnHot(0);
414 
415  FTensor::Tensor1<double, 3> t_coords[4] = {
416  FTensor::Tensor1<double, 3>(0., 0., 0.),
417  FTensor::Tensor1<double, 3>(1., 0., 0.),
418  FTensor::Tensor1<double, 3>(0., 1., 0.),
419  FTensor::Tensor1<double, 3>(0., 0., 1.)};
420 
421  FTensor::Tensor1<double *, 3> t_node_diff_ksi[4] = {
422  FTensor::Tensor1<double *, 3>(&diffN[0], &diffN[1], &diffN[2]),
423  FTensor::Tensor1<double *, 3>(&diffN[3], &diffN[4], &diffN[5]),
424  FTensor::Tensor1<double *, 3>(&diffN[6], &diffN[7], &diffN[8]),
425  FTensor::Tensor1<double *, 3>(&diffN[9], &diffN[10], &diffN[11])};
426 
427  FTENSOR_INDEX(3, i);
428  FTENSOR_INDEX(3, j);
429 
430  FTensor::Tensor1<double, 3> t_tau0i[4], t_tau0j[4];
431  FTensor::Tensor1<double, 3> t_diff_ksi0i[4], t_diff_ksi0j[4];
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);
440  }
441 
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]);
449  auto t_diff_phi_v_f = getFTensor2HVecFromPtr<3, 3>(diff_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;
456  FTensor::Tensor1<double, 3> t_diff_beta_f;
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);
464  FTensor::Tensor1<double, 3> t_diff_a;
465  int jj = 0;
466  for (int oo = 0; oo <= p - 3; oo++) {
467  auto t_psi_l = FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(psi_l);
468  auto t_diff_psi_l = FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3>(
469  diff_psi_l, &diff_psi_l[p + 1], &diff_psi_l[2 * p + 2]);
470  for (int l = 0; l <= oo; l++) {
471  int m = oo - l;
472  if (m >= 0) {
473  auto t_diff_psi_m = FTensor::Tensor1<double, 3>{
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);
478  ++t_phi_v_f;
479  ++jj;
480  t_phi_v_f(i) = a * t_tau0j[ff](i);
481  ++t_phi_v_f;
482  ++jj;
483 
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);
488  ++t_diff_phi_v_f;
489  t_diff_phi_v_f(i, j) = t_diff_a(j) * t_tau0j[ff](i);
490  ++t_diff_phi_v_f;
491 
492  ++t_psi_l;
493  ++t_diff_psi_l;
494  }
495  }
496  }
497  if (jj != NBVOLUMETET_AINSWORTH_FACE_HDIV(p)) {
498  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
499  "wrong order %d != %d", jj,
500  NBVOLUMETET_AINSWORTH_FACE_HDIV(p));
501  }
502  }
503  }
504  MoFEMFunctionReturnHot(0);
505 }
506 
507 MoFEMErrorCode MoFEM::Hdiv_Ainsworth_VolumeBubbleShapeFunctions_MBTET(
508  int p, double *N, double *diffN, double *phi_v, double *diff_phi_v,
509  int gdim,
510  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
511  double *L, double *diffL,
512  const int dim)) {
513 
514  MoFEMFunctionBeginHot;
515  if (p < 4)
516  MoFEMFunctionReturnHot(0);
517 
518  FTensor::Tensor1<double *, 3> t_node_diff_ksi[4] = {
519  FTensor::Tensor1<double *, 3>(&diffN[0], &diffN[1], &diffN[2]),
520  FTensor::Tensor1<double *, 3>(&diffN[3], &diffN[4], &diffN[5]),
521  FTensor::Tensor1<double *, 3>(&diffN[6], &diffN[7], &diffN[8]),
522  FTensor::Tensor1<double *, 3>(&diffN[9], &diffN[10], &diffN[11])};
523 
524  FTensor::Index<'i', 3> i;
525  FTensor::Index<'j', 3> j;
526  FTensor::Number<0> N0;
527  FTensor::Number<1> N1;
528  FTensor::Number<2> N2;
529 
530  FTensor::Tensor1<double, 3> t_diff_ksi0i;
531  FTensor::Tensor1<double, 3> t_diff_ksi0j;
532  FTensor::Tensor1<double, 3> t_diff_ksi0k;
533 
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);
537 
538  double psi_l[p + 1];
539  double diff_psi_l[3 * (p + 1)];
540  double psi_m[p + 1];
541  double diff_psi_m[3 * (p + 1)];
542  double psi_n[p + 1];
543  double diff_psi_n[3 * (p + 1)];
544 
545  auto t_phi_v = getFTensor1FromPtr<3>(phi_v);
546  auto t_diff_phi_v = getFTensor2HVecFromPtr<3, 3>(diff_phi_v);
547 
548  FTensor::Tensor1<double, 3> t_diff_beta_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);
566 
567  FTensor::Tensor1<double, 3> t_diff_a;
568 
569  int jj = 0;
570  for (int oo = 0; oo <= p - 4; oo++) {
571  auto t_psi_l = FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(psi_l);
572  auto t_diff_psi_l = FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3>(
573  diff_psi_l, &diff_psi_l[p + 1], &diff_psi_l[2 * p + 2]);
574  for (int l = 0; l <= oo; l++) {
575  auto t_psi_m = FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(psi_m);
576  auto t_diff_psi_m = FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3>(
577  diff_psi_m, &diff_psi_m[p + 1], &diff_psi_m[2 * p + 2]);
578  for (int m = 0; (l + m) <= oo; m++) {
579  int n = oo - l - m;
580  if (n >= 0) {
581  FTensor::Tensor1<double, 3> t_diff_psi_n(diff_psi_n[n],
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];
585  t_phi_v(0) = a;
586  t_phi_v(1) = 0;
587  t_phi_v(2) = 0;
588  ++t_phi_v;
589  t_phi_v(0) = 0;
590  t_phi_v(1) = a;
591  t_phi_v(2) = 0;
592  ++t_phi_v;
593  t_phi_v(0) = 0;
594  t_phi_v(1) = 0;
595  t_phi_v(2) = a;
596  ++t_phi_v;
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;
604  ++t_diff_phi_v;
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;
608  ++t_diff_phi_v;
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);
612  ++t_diff_phi_v;
613  ++jj;
614  }
615  ++t_psi_m;
616  ++t_diff_psi_m;
617  }
618  ++t_psi_l;
619  ++t_diff_psi_l;
620  }
621  }
622 
623  if (3 * jj != NBVOLUMETET_AINSWORTH_VOLUME_HDIV(p)) {
624  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
625  "wrong order %d != %d", jj,
626  NBVOLUMETET_AINSWORTH_VOLUME_HDIV(p));
627  }
628  }
629 
630  MoFEMFunctionReturnHot(0);
631 }
632 
633 MoFEMErrorCode
634 MoFEM::Hdiv_Demkowicz_Face_MBTET_ON_FACE(int *faces_nodes, int p, double *N,
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};
639 
640  MoFEMFunctionBeginHot;
641 
642  FTensor::Index<'i', 3> i;
643  FTensor::Index<'j', 3> j;
644 
645  FTensor::Tensor1<double, 3> t_cross[3];
646  FTensor::Tensor2<double, 3, 3> t_diff_cross(0., 0., 0., 0., 0., 0., 0., 0.,
647  0.);
648  FTensor::Tensor1<double, 3> t_node_diff_ksi[4];
649  FTensor::Tensor1<double, 3> t_node_diff_sum_n0_n1;
650 
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]]};
657 
658  FTensor::Tensor1<double, 3> t_diff_n0_p_n1;
659  FTensor::Tensor1<double, 3> t_diff_n0_p_n1_p_n2;
660 
661  if (diff_phi_f) {
662  t_node_diff_ksi[0] =
663  FTensor::Tensor1<double, 3>(diffN[0], diffN[1], diffN[2]);
664  t_node_diff_ksi[1] =
665  FTensor::Tensor1<double, 3>(diffN[3], diffN[4], diffN[5]);
666  t_node_diff_ksi[2] =
667  FTensor::Tensor1<double, 3>(diffN[6], diffN[7], diffN[8]);
668  t_node_diff_ksi[3] =
669  FTensor::Tensor1<double, 3>(diffN[9], diffN[10], diffN[11]);
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);
680  FTensor::Tensor1<double, 3> t_diff_o(
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);
683  // cerr << t_cross[ee](0) << " " << t_cross[ee](1) << " " <<
684  // t_cross[ee](2) << endl;
685  }
686  // cerr << endl << endl;
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);
689  } else {
690  for (int ee = 0; ee != 3; ee++) {
691  t_cross[ee](0) = 1;
692  t_cross[ee](1) = 0;
693  t_cross[ee](2) = 0;
694  }
695  }
696 
697  FTensor::Tensor1<double *, 3> t_phi(&phi_f[HVEC0], &phi_f[HVEC1],
698  &phi_f[HVEC2], 3);
699  boost::shared_ptr<FTensor::Tensor2<double *, 3, 3>> t_diff_phi_ptr;
700  if (diff_phi_f) {
701  t_diff_phi_ptr = boost::shared_ptr<FTensor::Tensor2<double *, 3, 3>>(
702  new FTensor::Tensor2<double *, 3, 3>(
703  &diff_phi_f[HVEC0_0], &diff_phi_f[HVEC0_1], &diff_phi_f[HVEC0_2],
704  &diff_phi_f[HVEC1_0], &diff_phi_f[HVEC1_1], &diff_phi_f[HVEC1_2],
705  &diff_phi_f[HVEC2_0], &diff_phi_f[HVEC2_1], &diff_phi_f[HVEC2_2],
706  9));
707  }
708 
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;
719  ierr = Jacobi_polynomials(p, 0, n1, n0 + n1, diff_n1, diff_n0_p_n1, fi,
720  diff_phi_f ? diff_fi : NULL, 3);
721  CHKERRG(ierr);
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;
725  ierr = Jacobi_polynomials(pp, 2 * pp + 1, n2, n0 + n1 + n2, diff_n2,
726  diff_n0_p_n1_p_n2, tmp_fj,
727  diff_phi_f ? tmp_diff_fj : NULL, 3);
728  CHKERRG(ierr);
729  fj[pp] = tmp_fj[pp];
730  if (diff_phi_f) {
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];
734  }
735  }
736  double no0 = N[shift + o[0]];
737  double no1 = N[shift + o[1]];
738  double no2 = N[shift + o[2]];
739  FTensor::Tensor1<double, 3> base0;
740  base0(i) = no0 * t_cross[0](i) + no1 * t_cross[1](i) + no2 * t_cross[2](i);
741  int jj = 0;
742  for (int oo = 0; oo < p; oo++) {
743  FTensor::Tensor0<double *> t_fi(fi);
744  FTensor::Tensor1<double *, 3> t_diff_fi(&diff_fi[0], &diff_fi[p + 1],
745  &diff_fi[2 * p + 2], 1);
746  for (int ll = 0; ll <= oo; ll++) {
747  const int mm = oo - ll;
748  if (mm >= 0) {
749  const double a = t_fi * fj[mm];
750  // cerr << ll << " " << mm << " " << a << endl;
751  t_phi(i) = a * base0(i);
752  if (diff_phi_f) {
753  FTensor::Tensor1<double *, 3> t_diff_fj(
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);
759  ++(*t_diff_phi_ptr);
760  ++t_diff_fi;
761  }
762  ++t_fi;
763  ++t_phi;
764  ++jj;
765  }
766  }
767  }
768  if (jj != NBFACETRI_DEMKOWICZ_HDIV(p)) {
769  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
770  "wrong number of base functions "
771  "jj!=NBFACETRI_DEMKOWICZ_HDIV(p) "
772  "%d!=%d",
773  jj, NBFACETRI_DEMKOWICZ_HDIV(p));
774  }
775  }
776 
777  MoFEMFunctionReturnHot(0);
778 }
779 
780 MoFEMErrorCode MoFEM::Hdiv_Demkowicz_Interior_MBTET(
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) {
783 
784  const int opposite_face_node[4] = {2, 0, 1, 3};
785  // list of zero node faces
786  const int znf[] = {0, 2, 3};
787  MoFEMFunctionBegin;
788 
789  FTensor::Index<'i', 3> i;
790  FTensor::Index<'j', 3> j;
791 
792  FTensor::Tensor1<double, 3> t_node_diff_ksi[4];
793  t_node_diff_ksi[0] =
794  FTensor::Tensor1<double, 3>(diffN[0], diffN[1], diffN[2]);
795  t_node_diff_ksi[1] =
796  FTensor::Tensor1<double, 3>(diffN[3], diffN[4], diffN[5]);
797  t_node_diff_ksi[2] =
798  FTensor::Tensor1<double, 3>(diffN[6], diffN[7], diffN[8]);
799  t_node_diff_ksi[3] =
800  FTensor::Tensor1<double, 3>(diffN[9], diffN[10], diffN[11]);
801  FTensor::Tensor1<double, 3> t_m_node_diff_ksi[4];
802  for (int ff = 0; ff != 4; ++ff) {
803  t_m_node_diff_ksi[ff](i) = -t_node_diff_ksi[ff](i);
804  }
805 
806  FTensor::Tensor1<double *, 3> t_phi_v(&phi_v[HVEC0], &phi_v[HVEC1],
807  &phi_v[HVEC2], 3);
808  FTensor::Tensor2<double *, 3, 3> t_diff_phi_v(
809  &diff_phi_v[HVEC0_0], &diff_phi_v[HVEC0_1], &diff_phi_v[HVEC0_2],
810  &diff_phi_v[HVEC1_0], &diff_phi_v[HVEC1_1], &diff_phi_v[HVEC1_2],
811  &diff_phi_v[HVEC2_0], &diff_phi_v[HVEC2_1], &diff_phi_v[HVEC2_2], 9);
812 
813  MatrixDouble fk(3, p + 1), diff_fk(3, 3 * p + 3);
814 
815  for (int ii = 0; ii != gdim; ii++) {
816  const int shift = 4 * ii;
817 
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++) {
823  CHKERR IntegratedJacobi_polynomials(
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);
826  }
827  }
828 
829  int jj = 0;
830  for (int oo = 2; oo <= p; oo++) {
831  for (int k = 1; k != oo; k++) {
832  int OO = oo - k;
833  if (OO >= 0) {
834  int s = NBFACETRI_DEMKOWICZ_HDIV(OO - 1);
835  // Note that we do faces 0,2,3, skipping 1. All the faces which have
836  // zero node in it.
837  int nb_dofs = NBFACETRI_DEMKOWICZ_HDIV(p_f[znf[0]]);
838  int sp[] = {ii * 3 * nb_dofs + 3 * s, ii * 3 * nb_dofs + 3 * s,
839  ii * 3 * nb_dofs + 3 * s};
840  FTensor::Tensor1<double *, 3> t_phi_f[] = {
841  FTensor::Tensor1<double *, 3>(&phi_f[znf[0]][sp[0] + HVEC0],
842  &phi_f[znf[0]][sp[0] + HVEC1],
843  &phi_f[znf[0]][sp[0] + HVEC2], 3),
844  FTensor::Tensor1<double *, 3>(&phi_f[znf[1]][sp[1] + HVEC0],
845  &phi_f[znf[1]][sp[1] + HVEC1],
846  &phi_f[znf[1]][sp[1] + HVEC2], 3),
847  FTensor::Tensor1<double *, 3>(&phi_f[znf[2]][sp[2] + HVEC0],
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};
852  FTensor::Tensor2<double *, 3, 3> t_diff_phi_f[] = {
853  FTensor::Tensor2<double *, 3, 3>(
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),
863  FTensor::Tensor2<double *, 3, 3>(
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),
873  FTensor::Tensor2<double *, 3, 3>(
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)};
883  for (int ij = s; ij != NBFACETRI_DEMKOWICZ_HDIV(OO); ij++) {
884  for (int ff = 0; ff != 3; ff++) {
885  FTensor::Tensor1<double, 3> t_diff_fk(diff_fk(ff, 0 * p + k - 1),
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);
891  ++t_phi_v;
892  ++t_diff_phi_v;
893  ++t_phi_f[ff];
894  ++t_diff_phi_f[ff];
895  ++jj;
896  }
897  }
898  }
899  }
900  }
901  if (jj != NBVOLUMETET_DEMKOWICZ_HDIV(p)) {
902  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
903  "wrong number of base functions "
904  "jj!=NBVOLUMETET_DEMKOWICZ_HDIV(p) "
905  "%d!=%d",
906  jj, NBVOLUMETET_DEMKOWICZ_HDIV(p));
907  }
908  }
909 
910  MoFEMFunctionReturn(0);
911 }
UBlasMatrix< double >
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::Hdiv_Ainsworth_FaceBasedVolumeShapeFunctions_MBTET
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))
Definition: Hdiv.cpp:401
MoFEM::Hdiv_Demkowicz_Face_MBTET_ON_FACE
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)
Definition: Hdiv.cpp:634
FTensor::Tensor1< double, 3 >
MoFEM::AinsworthOrderHooks::broken_nbvolumetet_volume_hdiv
static boost::function< int(int)> broken_nbvolumetet_volume_hdiv
Definition: Hdiv.hpp:29
HVEC0_2
@ HVEC0_2
Definition: definitions.h:211
MoFEM::Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET
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 .
Definition: Hdiv.cpp:24
NBFACETRI_AINSWORTH_EDGE_HDIV
#define NBFACETRI_AINSWORTH_EDGE_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:130
MoFEM::Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET_ON_FACE
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 .
Definition: Hdiv.cpp:47
MoFEM::AinsworthOrderHooks::broken_nbvolumetet_face_hdiv
static boost::function< int(int)> broken_nbvolumetet_face_hdiv
Definition: Hdiv.hpp:28
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
HVEC0_1
@ HVEC0_1
Definition: definitions.h:208
NBVOLUMETET_AINSWORTH_FACE_HDIV
#define NBVOLUMETET_AINSWORTH_FACE_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:134
MoFEM::Hdiv_Demkowicz_Interior_MBTET
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)
Definition: Hdiv.cpp:780
FTENSOR_INDEX
#define FTENSOR_INDEX(DIM, I)
Definition: Templates.hpp:2011
HVEC0_0
@ HVEC0_0
Definition: definitions.h:205
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
HVEC1_1
@ HVEC1_1
Definition: definitions.h:209
HVEC1
@ HVEC1
Definition: definitions.h:199
FTensor::Tensor2< double, 3, 3 >
MoFEM::AinsworthOrderHooks::broken_nbfacetri_edge_hdiv
static boost::function< int(int)> broken_nbfacetri_edge_hdiv
Definition: Hdiv.hpp:25
HVEC2_1
@ HVEC2_1
Definition: definitions.h:210
IntegratedJacobi_polynomials
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.
Definition: base_functions.c:134
MoFEM::AinsworthOrderHooks::broken_nbvolumetet_edge_hdiv
static boost::function< int(int)> broken_nbvolumetet_edge_hdiv
Definition: Hdiv.hpp:27
FTensor::Number< 0 >
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
HVEC2_2
@ HVEC2_2
Definition: definitions.h:213
HVEC1_2
@ HVEC1_2
Definition: definitions.h:212
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
NBVOLUMETET_AINSWORTH_VOLUME_HDIV
#define NBVOLUMETET_AINSWORTH_VOLUME_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:135
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEM::getFTensor2HVecFromPtr< 3, 3 >
FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > getFTensor2HVecFromPtr< 3, 3 >(double *ptr)
Definition: Templates.hpp:968
NBFACETRI_DEMKOWICZ_HDIV
#define NBFACETRI_DEMKOWICZ_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:139
NBVOLUMETET_DEMKOWICZ_HDIV
#define NBVOLUMETET_DEMKOWICZ_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:140
MoFEM::L
VectorDouble L
Definition: Projection10NodeCoordsOnField.cpp:124
MoFEM::Hdiv_Ainsworth_FaceBubbleShapeFunctions_ON_FACE
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 .
Definition: Hdiv.cpp:174
MoFEM::AinsworthOrderHooks::broken_nbfacetri_face_hdiv
static boost::function< int(int)> broken_nbfacetri_face_hdiv
Definition: Hdiv.hpp:26
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::Tensor0
Definition: Tensor0.hpp:16
MoFEM::Hdiv_Ainsworth_VolumeBubbleShapeFunctions_MBTET
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 .
Definition: Hdiv.cpp:507
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
Jacobi_polynomials
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.
Definition: base_functions.c:67
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::Hdiv_Ainsworth_FaceBubbleShapeFunctions
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 .
Definition: Hdiv.cpp:151
HVEC2_0
@ HVEC2_0
Definition: definitions.h:207
MoFEM::Hdiv_Ainsworth_EdgeBasedVolumeShapeFunctions_MBTET
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 .
Definition: Hdiv.cpp:307
NBFACETRI_AINSWORTH_FACE_HDIV
#define NBFACETRI_AINSWORTH_FACE_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:131
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
HVEC2
@ HVEC2
Definition: definitions.h:199
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
convert.int
int
Definition: convert.py:64
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
HVEC1_0
@ HVEC1_0
Definition: definitions.h:206
HVEC0
@ HVEC0
Definition: definitions.h:199
Generated by Doxygen 1.8.17 and hosted at University of Glasgow