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 MoFEMErrorCode MoFEM::Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET(
14  int *faces_nodes, int *p, double *N, double *diffN, double *phi_f_e[4][3],
15  double *diff_phi_f_e[4][3], int gdim,
16  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
17  double *L, double *diffL,
18  const int dim)) {
19 
20  MoFEMFunctionBeginHot;
21  for (int ff = 0; ff < 4; ff++) {
22  if (diff_phi_f_e != NULL) {
23  ierr = Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET_ON_FACE(
24  &faces_nodes[3 * ff], p[ff], N, diffN, phi_f_e[ff], diff_phi_f_e[ff],
25  gdim, 4, base_polynomials);
26  CHKERRG(ierr);
27  } else {
28  ierr = Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET_ON_FACE(
29  &faces_nodes[3 * ff], p[ff], N, diffN, phi_f_e[ff], NULL, gdim, 4,
30  base_polynomials);
31  CHKERRG(ierr);
32  }
33  }
34  MoFEMFunctionReturnHot(0);
35 }
36 
37 MoFEMErrorCode MoFEM::Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET_ON_FACE(
38  int *faces_nodes, int p, double *N, double *diffN, double *phi_f_e[3],
39  double *diff_phi_f_e[3], int gdim, int nb,
40  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
41  double *L, double *diffL,
42  const int dim)) {
43 
44  const int face_edges_nodes[3][2] = {{0, 1}, {1, 2}, {2, 0}};
45  const int face_oposite_edges_node[] = {2, 0, 1};
46  FTensor::Index<'i', 3> i;
47  FTensor::Index<'j', 3> j;
48 
49  MoFEMFunctionBeginHot;
50  if (p < 1)
51  MoFEMFunctionReturnHot(0);
52 
53  FTensor::Tensor1<double, 3> t_edge_cross[3];
54  FTensor::Tensor1<double, 3> t_node_diff_ksi[4];
55  FTensor::Tensor1<double, 3> t_diff_ksi0i[3];
56  if (diffN != NULL) {
57  t_node_diff_ksi[0] =
58  FTensor::Tensor1<double, 3>(diffN[0], diffN[1], diffN[2]);
59  t_node_diff_ksi[1] =
60  FTensor::Tensor1<double, 3>(diffN[3], diffN[4], diffN[5]);
61  t_node_diff_ksi[2] =
62  FTensor::Tensor1<double, 3>(diffN[6], diffN[7], diffN[8]);
63  t_node_diff_ksi[3] =
64  FTensor::Tensor1<double, 3>(diffN[9], diffN[10], diffN[11]);
65  for (int ee = 0; ee < 3; ee++) {
66  const int n0 = faces_nodes[face_edges_nodes[ee][0]];
67  const int n1 = faces_nodes[face_edges_nodes[ee][1]];
68  t_diff_ksi0i[ee](i) = t_node_diff_ksi[n1](i) - t_node_diff_ksi[n0](i);
69  t_edge_cross[ee](0) = t_node_diff_ksi[n0](1) * t_node_diff_ksi[n1](2) -
70  t_node_diff_ksi[n0](2) * t_node_diff_ksi[n1](1);
71  t_edge_cross[ee](1) = t_node_diff_ksi[n0](2) * t_node_diff_ksi[n1](0) -
72  t_node_diff_ksi[n0](0) * t_node_diff_ksi[n1](2);
73  t_edge_cross[ee](2) = t_node_diff_ksi[n0](0) * t_node_diff_ksi[n1](1) -
74  t_node_diff_ksi[n0](1) * t_node_diff_ksi[n1](0);
75  }
76  } else {
77  for (int ee = 0; ee < 3; ee++) {
78  t_edge_cross[ee](0) = 1;
79  t_edge_cross[ee](1) = 0;
80  t_edge_cross[ee](2) = 0;
81  }
82  }
83  double psi_l[p + 1], diff_psi_l[3 * (p + 1)];
84  boost::shared_ptr<FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3>>
85  t_diff_phi_f_e_ptr;
86 
87  for (int ee = 0; ee != 3; ee++) {
88  const int i0 = faces_nodes[face_edges_nodes[ee][0]];
89  const int i1 = faces_nodes[face_edges_nodes[ee][1]];
90  const int iO = faces_nodes[face_oposite_edges_node[ee]];
91  FTensor::Tensor1<double *, 3> t_psi_f_e(&phi_f_e[ee][0], &phi_f_e[ee][1],
92  &phi_f_e[ee][2], 3);
93  if (diff_phi_f_e) {
94  t_diff_phi_f_e_ptr = boost::shared_ptr<
95  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3>>(
96  new FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3>(
97  &diff_phi_f_e[ee][HVEC0_0], &diff_phi_f_e[ee][HVEC0_1],
98  &diff_phi_f_e[ee][HVEC0_2], &diff_phi_f_e[ee][HVEC1_0],
99  &diff_phi_f_e[ee][HVEC1_1], &diff_phi_f_e[ee][HVEC1_2],
100  &diff_phi_f_e[ee][HVEC2_0], &diff_phi_f_e[ee][HVEC2_1],
101  &diff_phi_f_e[ee][HVEC2_2]));
102  }
103  for (int ii = 0; ii != gdim; ii++) {
104  const int node_shift = ii * nb;
105  const double n0 = N[node_shift + i0];
106  const double n1 = N[node_shift + i1];
107  const double lambda = N[node_shift + iO];
108  const double ksi0i = n1 - n0;
109  if (diff_phi_f_e) {
110  ierr = base_polynomials(p, ksi0i, &t_diff_ksi0i[ee](0), psi_l,
111  diff_psi_l, 3);
112  CHKERRG(ierr);
113  } else {
114  ierr = base_polynomials(p, ksi0i, NULL, psi_l, NULL, 3);
115  CHKERRG(ierr);
116  }
117  FTensor::Tensor0<FTensor::PackPtr<double *, 1>> t_psi_l(&psi_l[0]);
118  FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3> t_diff_psi_l(
119  &diff_psi_l[0], &diff_psi_l[p + 1], &diff_psi_l[2 * p + 2]);
120  for (int l = 0; l <= p - 1; l++) {
121  t_psi_f_e(i) = lambda * t_psi_l * t_edge_cross[ee](i);
122  if (t_diff_phi_f_e_ptr) {
123  (*t_diff_phi_f_e_ptr)(i, j) =
124  (t_node_diff_ksi[iO](j) * t_psi_l + lambda * t_diff_psi_l(j)) *
125  t_edge_cross[ee](i);
126  ++t_diff_psi_l;
127  ++(*t_diff_phi_f_e_ptr);
128  }
129  ++t_psi_f_e;
130  ++t_psi_l;
131  }
132  }
133  }
134  MoFEMFunctionReturnHot(0);
135 }
136 
137 MoFEMErrorCode MoFEM::Hdiv_Ainsworth_FaceBubbleShapeFunctions(
138  int *faces_nodes, int *p, double *N, double *diffN, double *phi_f[],
139  double *diff_phi_f[], int gdim,
140  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
141  double *L, double *diffL,
142  const int dim)) {
143 
144  MoFEMFunctionBeginHot;
145  for (int ff = 0; ff < 4; ff++) {
146  double *diff;
147  if (diff_phi_f != NULL) {
148  diff = diff_phi_f[ff];
149  } else {
150  diff = NULL;
151  }
152  ierr = Hdiv_Ainsworth_FaceBubbleShapeFunctions_ON_FACE(
153  &faces_nodes[3 * ff], p[ff], N, diffN, phi_f[ff], diff, gdim, 4,
154  base_polynomials);
155  CHKERRG(ierr);
156  }
157  MoFEMFunctionReturnHot(0);
158 }
159 
160 MoFEMErrorCode MoFEM::Hdiv_Ainsworth_FaceBubbleShapeFunctions_ON_FACE(
161  int *face_nodes, int p, double *N, double *diffN, double *phi_f,
162  double *diff_phi_f, int gdim, int nb,
163  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
164  double *L, double *diffL,
165  const int dim)) {
166  FTensor::Index<'i', 3> i;
167  FTensor::Index<'j', 3> j;
168 
169  MoFEMFunctionBeginHot;
170  if (p < 3)
171  MoFEMFunctionReturnHot(0);
172 
173  const int vert_i = face_nodes[1];
174  const int vert_j = face_nodes[2];
175  const int i0 = face_nodes[0];
176  FTensor::Tensor1<double, 3> t_cross;
177  FTensor::Tensor1<double, 3> t_node_diff_ksi[4];
178  FTensor::Tensor1<double, 3> t_diff_ksi0i;
179  FTensor::Tensor1<double, 3> t_diff_ksi0j;
180 
181  if (diffN) {
182  t_node_diff_ksi[0] =
183  FTensor::Tensor1<double, 3>(diffN[0], diffN[1], diffN[2]);
184  t_node_diff_ksi[1] =
185  FTensor::Tensor1<double, 3>(diffN[3], diffN[4], diffN[5]);
186  t_node_diff_ksi[2] =
187  FTensor::Tensor1<double, 3>(diffN[6], diffN[7], diffN[8]);
188  t_node_diff_ksi[3] =
189  FTensor::Tensor1<double, 3>(diffN[9], diffN[10], diffN[11]);
190  t_diff_ksi0i(i) = t_node_diff_ksi[vert_i](i) - t_node_diff_ksi[i0](i);
191  t_diff_ksi0j(i) = t_node_diff_ksi[vert_j](i) - t_node_diff_ksi[i0](i);
192  t_cross(0) = t_node_diff_ksi[vert_i](1) * t_node_diff_ksi[vert_j](2) -
193  t_node_diff_ksi[vert_i](2) * t_node_diff_ksi[vert_j](1);
194  t_cross(1) = t_node_diff_ksi[vert_i](2) * t_node_diff_ksi[vert_j](0) -
195  t_node_diff_ksi[vert_i](0) * t_node_diff_ksi[vert_j](2);
196  t_cross(2) = t_node_diff_ksi[vert_i](0) * t_node_diff_ksi[vert_j](1) -
197  t_node_diff_ksi[vert_i](1) * t_node_diff_ksi[vert_j](0);
198  } else {
199  t_cross(0) = 1;
200  t_cross(1) = 0;
201  t_cross(2) = 0;
202  }
203 
204  double psi_l[p + 1], diff_psi_l[3 * (p + 1)];
205  double psi_m[p + 1], diff_psi_m[3 * (p + 1)];
206  FTensor::Tensor1<double, 3> t_diff_beta_0ij(0.,0.,0.);
207 
208  FTensor::Tensor1<double *, 3> t_psi_f(&phi_f[HVEC0], &phi_f[HVEC1],
209  &phi_f[HVEC2], 3);
210 
211  boost::shared_ptr<FTensor::Tensor2<double *, 3, 3> > t_diff_phi_f_ptr;
212  if (diff_phi_f) {
213  t_diff_phi_f_ptr = boost::shared_ptr<FTensor::Tensor2<double *, 3, 3> >(
214  new FTensor::Tensor2<double *, 3, 3>(
215  &diff_phi_f[HVEC0_0], &diff_phi_f[HVEC0_1], &diff_phi_f[HVEC0_2],
216  &diff_phi_f[HVEC1_0], &diff_phi_f[HVEC1_1], &diff_phi_f[HVEC1_2],
217  &diff_phi_f[HVEC2_0], &diff_phi_f[HVEC2_1], &diff_phi_f[HVEC2_2],
218  9));
219  }
220 
221  for (int ii = 0; ii < gdim; ii++) {
222 
223  int node_shift = ii * nb;
224  const double ni = N[node_shift + vert_i];
225  const double nj = N[node_shift + vert_j];
226  const double n0 = N[node_shift + i0];
227  const double ksi0i = ni - n0;
228  const double ksi0j = nj - n0;
229  double beta_0ij = n0 * ni * nj;
230  if (diff_phi_f) {
231  t_diff_beta_0ij(i) = (ni * nj) * t_node_diff_ksi[i0](i) +
232  (n0 * nj) * t_node_diff_ksi[vert_i](i) +
233  (n0 * ni) * t_node_diff_ksi[vert_j](i);
234  ierr = base_polynomials(p, ksi0i, &t_diff_ksi0i(0), psi_l, diff_psi_l, 3);
235  CHKERRG(ierr);
236  ierr = base_polynomials(p, ksi0j, &t_diff_ksi0j(0), psi_m, diff_psi_m, 3);
237  CHKERRG(ierr);
238  } else {
239  ierr = base_polynomials(p, ksi0i, NULL, psi_l, NULL, 3);
240  CHKERRG(ierr);
241  ierr = base_polynomials(p, ksi0j, NULL, psi_m, NULL, 3);
242  CHKERRG(ierr);
243  }
244 
245  int jj = 0;
246  int oo = 0;
247  for (; oo <= p - 3; oo++) {
248  FTensor::Tensor0<double *> t_psi_l(&psi_l[0]);
249  FTensor::Tensor1<double *, 3> t_diff_psi_l(diff_psi_l, &diff_psi_l[p + 1],
250  &diff_psi_l[2 * p + 2], 1);
251  for (int l = 0; l <= oo; l++) {
252  int m = oo - l;
253  if (m >= 0) {
254  FTensor::Tensor1<double, 3> t_diff_psi_m(
255  diff_psi_m[m], diff_psi_m[p + 1 + m], diff_psi_m[2 * p + 2 + m]);
256  t_psi_f(i) = (beta_0ij * t_psi_l * psi_m[m]) * t_cross(i);
257  ++t_psi_f;
258  if (diff_phi_f) {
259  (*t_diff_phi_f_ptr)(i, j) =
260  ((t_psi_l * psi_m[m]) * t_diff_beta_0ij(j) +
261  (beta_0ij * psi_m[m]) * t_diff_psi_l(j) +
262  (beta_0ij * t_psi_l) * t_diff_psi_m(j)) *
263  t_cross(i);
264  ++(*t_diff_phi_f_ptr);
265  }
266  }
267  ++t_psi_l;
268  ++t_diff_psi_l;
269  ++jj;
270  }
271  }
272  if (jj != NBFACETRI_AINSWORTH_FACE_HDIV(p)) {
273  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
274  "wrong order %d != %d", jj, NBFACETRI_AINSWORTH_FACE_HDIV(p));
275  }
276  }
277  MoFEMFunctionReturnHot(0);
278 }
279 
280 MoFEMErrorCode MoFEM::Hdiv_Ainsworth_EdgeBasedVolumeShapeFunctions_MBTET(
281  int p, double *N, double *diffN, double *phi_v_e[6],
282  double *diff_phi_v_e[6], int gdim,
283  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
284  double *L, double *diffL,
285  const int dim)) {
286 
287  const int edges_nodes[6][2] = {{0, 1}, {1, 2}, {2, 0},
288  {0, 3}, {1, 3}, {2, 3}};
289 
290  MoFEMFunctionBeginHot;
291  if (p < 2)
292  MoFEMFunctionReturnHot(0);
293  if (diffN == NULL) {
294  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
295  }
296 
297  FTensor::Tensor1<double, 3> t_coords[4] = {
298  FTensor::Tensor1<double, 3>(0., 0., 0.),
299  FTensor::Tensor1<double, 3>(1., 0., 0.),
300  FTensor::Tensor1<double, 3>(0., 1., 0.),
301  FTensor::Tensor1<double, 3>(0., 0., 1.)};
302  FTensor::Tensor1<double *, 3> t_node_diff_ksi[4] = {
303  FTensor::Tensor1<double *, 3>(&diffN[0], &diffN[1], &diffN[2]),
304  FTensor::Tensor1<double *, 3>(&diffN[3], &diffN[4], &diffN[5]),
305  FTensor::Tensor1<double *, 3>(&diffN[6], &diffN[7], &diffN[8]),
306  FTensor::Tensor1<double *, 3>(&diffN[9], &diffN[10], &diffN[11])};
307 
308  FTensor::Index<'i', 3> i;
309  FTensor::Index<'j', 3> j;
310 
311  FTensor::Tensor1<double, 3> t_tou_e;
312  FTensor::Tensor1<double, 3> t_diff_ksi0i;
313  FTensor::Tensor1<double, 3> t_diff_beta_e;
314 
315  double psi_l[p + 1];
316  double diff_psi_l[3 * (p + 1)];
317 
318  for (int ee = 0; ee != 6; ee++) {
319  t_tou_e(i) =
320  t_coords[edges_nodes[ee][1]](i) - t_coords[edges_nodes[ee][0]](i);
321  t_diff_ksi0i(i) = t_node_diff_ksi[edges_nodes[ee][1]](i) -
322  t_node_diff_ksi[edges_nodes[ee][0]](i);
323  FTensor::Tensor1<double *, 3> t_psi_v_e(&phi_v_e[ee][0], &phi_v_e[ee][1],
324  &phi_v_e[ee][2], 3);
325  FTensor::Tensor2<double *, 3, 3> t_diff_phi_v_e(
326  &diff_phi_v_e[ee][HVEC0_0], &diff_phi_v_e[ee][HVEC0_1],
327  &diff_phi_v_e[ee][HVEC0_2], &diff_phi_v_e[ee][HVEC1_0],
328  &diff_phi_v_e[ee][HVEC1_1], &diff_phi_v_e[ee][HVEC1_2],
329  &diff_phi_v_e[ee][HVEC2_0], &diff_phi_v_e[ee][HVEC2_1],
330  &diff_phi_v_e[ee][HVEC2_2], 9);
331  for (int ii = 0; ii != gdim; ii++) {
332  const int node_shift = ii * 4;
333  const double ni = N[node_shift + edges_nodes[ee][1]];
334  const double n0 = N[node_shift + edges_nodes[ee][0]];
335  const double beta_e = ni * n0;
336  const double ksi0i = ni - n0;
337  if (diff_phi_v_e) {
338  t_diff_beta_e(i) = ni * t_node_diff_ksi[edges_nodes[ee][0]](i) +
339  t_node_diff_ksi[edges_nodes[ee][1]](i) * n0;
340  ierr =
341  base_polynomials(p, ksi0i, &t_diff_ksi0i(0), psi_l, diff_psi_l, 3);
342  CHKERRG(ierr);
343  } else {
344  ierr = base_polynomials(p, ksi0i, NULL, psi_l, NULL, 3);
345  CHKERRG(ierr);
346  }
347  FTensor::Tensor0<double *> t_psi_l(&psi_l[0]);
348  FTensor::Tensor1<double *, 3> t_diff_psi_l(
349  &diff_psi_l[0], &diff_psi_l[p + 1], &diff_psi_l[2 * p + 2], 1);
350  for (int l = 0; l <= p - 2; l++) {
351  t_psi_v_e(i) = (beta_e * t_psi_l) * t_tou_e(i);
352  ++t_psi_v_e;
353  if (diff_phi_v_e) {
354  t_diff_phi_v_e(i, j) =
355  (t_diff_beta_e(j) * t_psi_l + beta_e * t_diff_psi_l(j)) *
356  t_tou_e(i);
357  ++t_diff_phi_v_e;
358  ++t_diff_psi_l;
359  }
360  ++t_psi_l;
361  }
362  }
363  }
364 
365  MoFEMFunctionReturnHot(0);
366 }
367 
368 MoFEMErrorCode MoFEM::Hdiv_Ainsworth_FaceBasedVolumeShapeFunctions_MBTET(
369  int p, double *N, double *diffN, double *phi_v_f[], double *diff_phi_v_f[],
370  int gdim,
371  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
372  double *L, double *diffL,
373  const int dim)) {
374 
375  const int faces_nodes[4][3] = {{0, 1, 3}, {1, 2, 3}, {0, 2, 3}, {0, 1, 2}};
376 
377  MoFEMFunctionBeginHot;
378  if (p < 3)
379  MoFEMFunctionReturnHot(0);
380 
381  FTensor::Tensor1<double, 3> t_coords[4] = {
382  FTensor::Tensor1<double, 3>(0., 0., 0.),
383  FTensor::Tensor1<double, 3>(1., 0., 0.),
384  FTensor::Tensor1<double, 3>(0., 1., 0.),
385  FTensor::Tensor1<double, 3>(0., 0., 1.)};
386 
387  FTensor::Tensor1<double *, 3> t_node_diff_ksi[4] = {
388  FTensor::Tensor1<double *, 3>(&diffN[0], &diffN[1], &diffN[2]),
389  FTensor::Tensor1<double *, 3>(&diffN[3], &diffN[4], &diffN[5]),
390  FTensor::Tensor1<double *, 3>(&diffN[6], &diffN[7], &diffN[8]),
391  FTensor::Tensor1<double *, 3>(&diffN[9], &diffN[10], &diffN[11])};
392 
393  FTensor::Index<'i', 3> i;
394  FTensor::Index<'j', 3> j;
395 
396  FTensor::Tensor1<double, 3> t_tau0i[4], t_tau0j[4];
397  FTensor::Tensor1<double, 3> t_diff_ksi0i[4], t_diff_ksi0j[4];
398  for (int ff = 0; ff != 4; ff++) {
399  const int v0 = faces_nodes[ff][0];
400  const int vi = faces_nodes[ff][1];
401  const int vj = faces_nodes[ff][2];
402  t_tau0i[ff](i) = t_coords[vi](i) - t_coords[v0](i);
403  t_tau0j[ff](i) = t_coords[vj](i) - t_coords[v0](i);
404  t_diff_ksi0i[ff](i) = t_node_diff_ksi[vi](i) - t_node_diff_ksi[v0](i);
405  t_diff_ksi0j[ff](i) = t_node_diff_ksi[vj](i) - t_node_diff_ksi[v0](i);
406  }
407 
408  double psi_l[p + 1], psi_m[p + 1];
409  double diff_psi_l[3 * (p + 1)], diff_psi_m[3 * (p + 1)];
410  for (int ff = 0; ff != 4; ff++) {
411  const int v0 = faces_nodes[ff][0];
412  const int vi = faces_nodes[ff][1];
413  const int vj = faces_nodes[ff][2];
414  FTensor::Tensor1<double *, 3> t_phi_v_f(
415  &phi_v_f[ff][HVEC0], &phi_v_f[ff][HVEC1], &phi_v_f[ff][HVEC2], 3);
416  FTensor::Tensor2<double *, 3, 3> t_diff_phi_v_f(
417  &diff_phi_v_f[ff][HVEC0_0], &diff_phi_v_f[ff][HVEC0_1],
418  &diff_phi_v_f[ff][HVEC0_2], &diff_phi_v_f[ff][HVEC1_0],
419  &diff_phi_v_f[ff][HVEC1_1], &diff_phi_v_f[ff][HVEC1_2],
420  &diff_phi_v_f[ff][HVEC2_0], &diff_phi_v_f[ff][HVEC2_1],
421  &diff_phi_v_f[ff][HVEC2_2], 9);
422  for (int ii = 0; ii < gdim; ii++) {
423  const int node_shift = 4 * ii;
424  const double n0 = N[node_shift + v0];
425  const double ni = N[node_shift + vi];
426  const double nj = N[node_shift + vj];
427  const double beta_f = n0 * ni * nj;
428  FTensor::Tensor1<double, 3> t_diff_beta_f;
429  t_diff_beta_f(i) = (ni * nj) * t_node_diff_ksi[v0](i) +
430  (n0 * nj) * t_node_diff_ksi[vi](i) +
431  (n0 * ni) * t_node_diff_ksi[vj](i);
432  const double ksi0i = ni - n0;
433  const double ksi0j = nj - n0;
434  ierr = base_polynomials(p, ksi0i, &t_diff_ksi0i[ff](0), psi_l, diff_psi_l,
435  3);
436  CHKERRG(ierr);
437  ierr = base_polynomials(p, ksi0j, &t_diff_ksi0j[ff](0), psi_m, diff_psi_m,
438  3);
439  CHKERRG(ierr);
440  FTensor::Tensor1<double, 3> t_diff_a;
441  int jj = 0;
442  for (int oo = 0; oo <= p - 3; oo++) {
443  FTensor::Tensor0<double *> t_psi_l(&psi_l[0], 1);
444  FTensor::Tensor1<double *, 3> t_diff_psi_l(
445  diff_psi_l, &diff_psi_l[p + 1], &diff_psi_l[2 * p + 2], 1);
446  for (int l = 0; l <= oo; l++) {
447  int m = oo - l;
448  if (m >= 0) {
449 
450  FTensor::Tensor1<double *, 3> t_diff_psi_m(
451  &diff_psi_m[m], &diff_psi_m[p + 1 + m],
452  &diff_psi_m[2 * p + 2 + m], 1);
453  const double a = beta_f * t_psi_l * psi_m[m];
454  t_phi_v_f(i) = a * t_tau0i[ff](i);
455  ++t_phi_v_f;
456  ++jj;
457  t_phi_v_f(i) = a * t_tau0j[ff](i);
458  ++t_phi_v_f;
459  ++jj;
460 
461  t_diff_a(j) = (t_psi_l * psi_m[m]) * t_diff_beta_f(j) +
462  (beta_f * psi_m[m]) * t_diff_psi_l(j) +
463  (beta_f * t_psi_l) * t_diff_psi_m(j);
464  t_diff_phi_v_f(i, j) = t_diff_a(j) * t_tau0i[ff](i);
465  ++t_diff_phi_v_f;
466  t_diff_phi_v_f(i, j) = t_diff_a(j) * t_tau0j[ff](i);
467  ++t_diff_phi_v_f;
468 
469  ++t_psi_l;
470  ++t_diff_psi_l;
471  }
472  }
473  }
474  if (jj != NBVOLUMETET_AINSWORTH_FACE_HDIV(p)) {
475  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
476  "wrong order %d != %d", jj,
477  NBVOLUMETET_AINSWORTH_FACE_HDIV(p));
478  }
479  }
480  }
481  MoFEMFunctionReturnHot(0);
482 }
483 
484 MoFEMErrorCode MoFEM::Hdiv_Ainsworth_VolumeBubbleShapeFunctions_MBTET(
485  int p, double *N, double *diffN, double *phi_v, double *diff_phi_v,
486  int gdim,
487  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s,
488  double *L, double *diffL,
489  const int dim)) {
490 
491  MoFEMFunctionBeginHot;
492  if (p < 4)
493  MoFEMFunctionReturnHot(0);
494 
495  FTensor::Tensor1<double *, 3> t_node_diff_ksi[4] = {
496  FTensor::Tensor1<double *, 3>(&diffN[0], &diffN[1], &diffN[2]),
497  FTensor::Tensor1<double *, 3>(&diffN[3], &diffN[4], &diffN[5]),
498  FTensor::Tensor1<double *, 3>(&diffN[6], &diffN[7], &diffN[8]),
499  FTensor::Tensor1<double *, 3>(&diffN[9], &diffN[10], &diffN[11])};
500 
501  FTensor::Index<'i', 3> i;
502  FTensor::Index<'j', 3> j;
503  FTensor::Number<0> N0;
504  FTensor::Number<1> N1;
505  FTensor::Number<2> N2;
506 
507  FTensor::Tensor1<double, 3> t_diff_ksi0i;
508  FTensor::Tensor1<double, 3> t_diff_ksi0j;
509  FTensor::Tensor1<double, 3> t_diff_ksi0k;
510 
511  t_diff_ksi0i(i) = t_node_diff_ksi[1](i) - t_node_diff_ksi[0](i);
512  t_diff_ksi0j(i) = t_node_diff_ksi[2](i) - t_node_diff_ksi[0](i);
513  t_diff_ksi0k(i) = t_node_diff_ksi[3](i) - t_node_diff_ksi[0](i);
514 
515  double psi_l[p + 1];
516  double diff_psi_l[3 * (p + 1)];
517  double psi_m[p + 1];
518  double diff_psi_m[3 * (p + 1)];
519  double psi_n[p + 1];
520  double diff_psi_n[3 * (p + 1)];
521 
522  FTensor::Tensor1<double *, 3> t_phi_v(phi_v, &phi_v[HVEC1], &phi_v[HVEC2], 3);
523  FTensor::Tensor2<double *, 3, 3> t_diff_phi_v(
524  diff_phi_v, &diff_phi_v[HVEC0_1], &diff_phi_v[HVEC0_2],
525  &diff_phi_v[HVEC1_0], &diff_phi_v[HVEC1_1], &diff_phi_v[HVEC1_2],
526  &diff_phi_v[HVEC2_0], &diff_phi_v[HVEC2_1], &diff_phi_v[HVEC2_2], 9);
527 
528  FTensor::Tensor1<double, 3> t_diff_beta_v;
529  for (int ii = 0; ii < gdim; ii++) {
530  const int node_shift = ii * 4;
531  const double n0 = N[node_shift + 0];
532  const double ni = N[node_shift + 1];
533  const double nj = N[node_shift + 2];
534  const double nk = N[node_shift + 3];
535  const double ksi0i = ni - n0;
536  const double ksi0j = nj - n0;
537  const double ksi0k = nk - n0;
538  const double beta_v = n0 * ni * nj * nk;
539  t_diff_beta_v(i) = (ni * nj * nk) * t_node_diff_ksi[0](i) +
540  (n0 * nj * nk) * t_node_diff_ksi[1](i) +
541  (n0 * ni * nk) * t_node_diff_ksi[2](i) +
542  (n0 * ni * nj) * t_node_diff_ksi[3](i);
543  ierr = base_polynomials(p, ksi0i, &t_diff_ksi0i(0), psi_l, diff_psi_l, 3);
544  CHKERRG(ierr);
545  ierr = base_polynomials(p, ksi0j, &t_diff_ksi0j(0), psi_m, diff_psi_m, 3);
546  CHKERRG(ierr);
547  ierr = base_polynomials(p, ksi0k, &t_diff_ksi0k(0), psi_n, diff_psi_n, 3);
548  CHKERRG(ierr);
549 
550  FTensor::Tensor1<double, 3> t_diff_a;
551 
552  int jj = 0;
553  for (int oo = 0; oo <= p - 4; oo++) {
554  FTensor::Tensor0<double *> t_psi_l(&psi_l[0]);
555  FTensor::Tensor1<double *, 3> t_diff_psi_l(diff_psi_l, &diff_psi_l[p + 1],
556  &diff_psi_l[2 * p + 2], 1);
557  for (int l = 0; l <= oo; l++) {
558  FTensor::Tensor0<double *> t_psi_m(&psi_m[0]);
559  FTensor::Tensor1<double *, 3> t_diff_psi_m(
560  diff_psi_m, &diff_psi_m[p + 1], &diff_psi_m[2 * p + 2], 1);
561  for (int m = 0; (l + m) <= oo; m++) {
562  int n = oo - l - m;
563  if (n >= 0) {
564  FTensor::Tensor1<double, 3> t_diff_psi_n(diff_psi_n[n],
565  diff_psi_n[p + 1 + n],
566  diff_psi_n[2 * p + 2 + n]);
567  const double a = beta_v * t_psi_l * t_psi_m * psi_n[n];
568  t_phi_v(0) = a;
569  t_phi_v(1) = 0;
570  t_phi_v(2) = 0;
571  ++t_phi_v;
572  t_phi_v(0) = 0;
573  t_phi_v(1) = a;
574  t_phi_v(2) = 0;
575  ++t_phi_v;
576  t_phi_v(0) = 0;
577  t_phi_v(1) = 0;
578  t_phi_v(2) = a;
579  ++t_phi_v;
580  t_diff_a(j) = (t_psi_l * t_psi_m * psi_n[n]) * t_diff_beta_v(j) +
581  (beta_v * t_psi_m * psi_n[n]) * t_diff_psi_l(j) +
582  (beta_v * t_psi_l * psi_n[n]) * t_diff_psi_m(j) +
583  (beta_v * t_psi_l * t_psi_m) * t_diff_psi_n(j);
584  t_diff_phi_v(N0, j) = t_diff_a(j);
585  t_diff_phi_v(N1, j) = 0;
586  t_diff_phi_v(N2, j) = 0;
587  ++t_diff_phi_v;
588  t_diff_phi_v(N0, j) = 0;
589  t_diff_phi_v(N1, j) = t_diff_a(j);
590  t_diff_phi_v(N2, j) = 0;
591  ++t_diff_phi_v;
592  t_diff_phi_v(N0, j) = 0;
593  t_diff_phi_v(N1, j) = 0;
594  t_diff_phi_v(N2, j) = t_diff_a(j);
595  ++t_diff_phi_v;
596  ++jj;
597  }
598  ++t_psi_m;
599  ++t_diff_psi_m;
600  }
601  ++t_psi_l;
602  ++t_diff_psi_l;
603  }
604  }
605 
606  if (3 * jj != NBVOLUMETET_AINSWORTH_VOLUME_HDIV(p)) {
607  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
608  "wrong order %d != %d", jj,
609  NBVOLUMETET_AINSWORTH_VOLUME_HDIV(p));
610  }
611  }
612 
613  MoFEMFunctionReturnHot(0);
614 }
615 
616 MoFEMErrorCode
617 MoFEM::Hdiv_Demkowicz_Face_MBTET_ON_FACE(int *faces_nodes, int p, double *N,
618  double *diffN, double *phi_f,
619  double *diff_phi_f, int gdim, int nb) {
620  const int face_edges_nodes[3][2] = {{0, 1}, {1, 2}, {2, 0}};
621  const int face_oposite_edges_node[] = {2, 0, 1};
622 
623  MoFEMFunctionBeginHot;
624 
625  FTensor::Index<'i', 3> i;
626  FTensor::Index<'j', 3> j;
627 
628  FTensor::Tensor1<double, 3> t_cross[3];
629  FTensor::Tensor2<double, 3, 3> t_diff_cross(0.,0.,0., 0.,0.,0., 0.,0.,0.);
630  FTensor::Tensor1<double, 3> t_node_diff_ksi[4];
631  FTensor::Tensor1<double, 3> t_node_diff_sum_n0_n1;
632 
633  const int i0 = faces_nodes[0];
634  const int i1 = faces_nodes[1];
635  const int i2 = faces_nodes[2];
636  const int o[] = {faces_nodes[face_oposite_edges_node[0]],
637  faces_nodes[face_oposite_edges_node[1]],
638  faces_nodes[face_oposite_edges_node[2]]};
639 
640  FTensor::Tensor1<double, 3> t_diff_n0_p_n1;
641  FTensor::Tensor1<double, 3> t_diff_n0_p_n1_p_n2;
642 
643  if (diff_phi_f) {
644  t_node_diff_ksi[0] =
645  FTensor::Tensor1<double, 3>(diffN[0], diffN[1], diffN[2]);
646  t_node_diff_ksi[1] =
647  FTensor::Tensor1<double, 3>(diffN[3], diffN[4], diffN[5]);
648  t_node_diff_ksi[2] =
649  FTensor::Tensor1<double, 3>(diffN[6], diffN[7], diffN[8]);
650  t_node_diff_ksi[3] =
651  FTensor::Tensor1<double, 3>(diffN[9], diffN[10], diffN[11]);
652  t_diff_cross(i, j) = 0;
653  for (int ee = 0; ee != 3; ee++) {
654  int ei0 = faces_nodes[face_edges_nodes[ee][0]];
655  int ei1 = faces_nodes[face_edges_nodes[ee][1]];
656  t_cross[ee](0) = t_node_diff_ksi[ei0](1) * t_node_diff_ksi[ei1](2) -
657  t_node_diff_ksi[ei0](2) * t_node_diff_ksi[ei1](1);
658  t_cross[ee](1) = t_node_diff_ksi[ei0](2) * t_node_diff_ksi[ei1](0) -
659  t_node_diff_ksi[ei0](0) * t_node_diff_ksi[ei1](2);
660  t_cross[ee](2) = t_node_diff_ksi[ei0](0) * t_node_diff_ksi[ei1](1) -
661  t_node_diff_ksi[ei0](1) * t_node_diff_ksi[ei1](0);
662  FTensor::Tensor1<double, 3> t_diff_o(
663  diffN[3 * o[ee] + 0], diffN[3 * o[ee] + 1], diffN[3 * o[ee] + 2]);
664  t_diff_cross(i, j) += t_cross[ee](i) * t_diff_o(j);
665  // cerr << t_cross[ee](0) << " " << t_cross[ee](1) << " " <<
666  // t_cross[ee](2) << endl;
667  }
668  // cerr << endl << endl;
669  t_diff_n0_p_n1(i) = t_node_diff_ksi[i0](i) + t_node_diff_ksi[i1](i);
670  t_diff_n0_p_n1_p_n2(i) = t_diff_n0_p_n1(i) + t_node_diff_ksi[i2](i);
671  } else {
672  for (int ee = 0; ee != 3; ee++) {
673  t_cross[ee](0) = 1;
674  t_cross[ee](1) = 0;
675  t_cross[ee](2) = 0;
676  }
677  }
678 
679  FTensor::Tensor1<double *, 3> t_phi(&phi_f[HVEC0], &phi_f[HVEC1],
680  &phi_f[HVEC2], 3);
681  boost::shared_ptr<FTensor::Tensor2<double *, 3, 3> > t_diff_phi_ptr;
682  if (diff_phi_f) {
683  t_diff_phi_ptr = boost::shared_ptr<FTensor::Tensor2<double *, 3, 3> >(
684  new FTensor::Tensor2<double *, 3, 3>(
685  &diff_phi_f[HVEC0_0], &diff_phi_f[HVEC0_1], &diff_phi_f[HVEC0_2],
686  &diff_phi_f[HVEC1_0], &diff_phi_f[HVEC1_1], &diff_phi_f[HVEC1_2],
687  &diff_phi_f[HVEC2_0], &diff_phi_f[HVEC2_1], &diff_phi_f[HVEC2_2],
688  9));
689  }
690 
691  double fi[p + 1], diff_fi[3 * p + 3];
692  double fj[p + 1], diff_fj[3 * p + 3];
693  double tmp_fj[p + 1], tmp_diff_fj[3 * p + 3];
694  for (int ii = 0; ii != gdim; ii++) {
695  const int shift = ii * nb;
696  double n0 = N[shift + i0];
697  double n1 = N[shift + i1];
698  double n2 = N[shift + i2];
699  double *diff_n1 = (diff_phi_f) ? &t_node_diff_ksi[i1](0) : NULL;
700  double *diff_n0_p_n1 = (diff_phi_f) ? &t_diff_n0_p_n1(0) : NULL;
701  ierr = Jacobi_polynomials(p, 0, n1, n0 + n1, diff_n1, diff_n0_p_n1, fi,
702  diff_phi_f ? diff_fi : NULL, 3);
703  CHKERRG(ierr);
704  for (int pp = 0; pp <= p; pp++) {
705  double *diff_n2 = (diff_phi_f) ? &t_node_diff_ksi[i2](0) : NULL;
706  double *diff_n0_p_n1_p_n2 = (diff_phi_f) ? &t_diff_n0_p_n1_p_n2(0) : NULL;
707  ierr = Jacobi_polynomials(pp, 2 * pp + 1, n2, n0 + n1 + n2, diff_n2,
708  diff_n0_p_n1_p_n2, tmp_fj,
709  diff_phi_f ? tmp_diff_fj : NULL, 3);
710  CHKERRG(ierr);
711  fj[pp] = tmp_fj[pp];
712  if (diff_phi_f) {
713  diff_fj[0 * (p + 1) + pp] = tmp_diff_fj[0 * (pp + 1) + pp];
714  diff_fj[1 * (p + 1) + pp] = tmp_diff_fj[1 * (pp + 1) + pp];
715  diff_fj[2 * (p + 1) + pp] = tmp_diff_fj[2 * (pp + 1) + pp];
716  }
717  }
718  double no0 = N[shift + o[0]];
719  double no1 = N[shift + o[1]];
720  double no2 = N[shift + o[2]];
721  FTensor::Tensor1<double, 3> base0;
722  base0(i) = no0 * t_cross[0](i) + no1 * t_cross[1](i) + no2 * t_cross[2](i);
723  int jj = 0;
724  for (int oo = 0; oo < p; oo++) {
725  FTensor::Tensor0<double *> t_fi(fi);
726  FTensor::Tensor1<double *, 3> t_diff_fi(&diff_fi[0], &diff_fi[p + 1],
727  &diff_fi[2 * p + 2], 1);
728  for (int ll = 0; ll <= oo; ll++) {
729  const int mm = oo - ll;
730  if (mm >= 0) {
731  const double a = t_fi * fj[mm];
732  // cerr << ll << " " << mm << " " << a << endl;
733  t_phi(i) = a * base0(i);
734  if (diff_phi_f) {
735  FTensor::Tensor1<double *, 3> t_diff_fj(
736  &diff_fj[0 + mm], &diff_fj[p + 1 + mm],
737  &diff_fj[2 * p + 2 + mm], 1);
738  (*t_diff_phi_ptr)(i, j) =
739  a * t_diff_cross(i, j) +
740  (t_diff_fi(j) * fj[mm] + t_fi * t_diff_fj(j)) * base0(i);
741  ++(*t_diff_phi_ptr);
742  ++t_diff_fi;
743  }
744  ++t_fi;
745  ++t_phi;
746  ++jj;
747  }
748  }
749  }
750  if (jj != NBFACETRI_DEMKOWICZ_HDIV(p)) {
751  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
752  "wrong number of base functions "
753  "jj!=NBFACETRI_DEMKOWICZ_HDIV(p) "
754  "%d!=%d",
755  jj, NBFACETRI_DEMKOWICZ_HDIV(p));
756  }
757  }
758 
759  MoFEMFunctionReturnHot(0);
760 }
761 
762 MoFEMErrorCode MoFEM::Hdiv_Demkowicz_Interior_MBTET(
763  int p, double *N, double *diffN, int p_f[], double *phi_f[4],
764  double *diff_phi_f[4], double *phi_v, double *diff_phi_v, int gdim) {
765 
766  const int opposite_face_node[4] = {2, 0, 1, 3};
767  // list of zero node faces
768  const int znf[] = {0, 2, 3};
769  MoFEMFunctionBegin;
770 
771  FTensor::Index<'i', 3> i;
772  FTensor::Index<'j', 3> j;
773 
774  FTensor::Tensor1<double, 3> t_node_diff_ksi[4];
775  t_node_diff_ksi[0] =
776  FTensor::Tensor1<double, 3>(diffN[0], diffN[1], diffN[2]);
777  t_node_diff_ksi[1] =
778  FTensor::Tensor1<double, 3>(diffN[3], diffN[4], diffN[5]);
779  t_node_diff_ksi[2] =
780  FTensor::Tensor1<double, 3>(diffN[6], diffN[7], diffN[8]);
781  t_node_diff_ksi[3] =
782  FTensor::Tensor1<double, 3>(diffN[9], diffN[10], diffN[11]);
783  FTensor::Tensor1<double, 3> t_m_node_diff_ksi[4];
784  for (int ff = 0; ff != 4; ++ff) {
785  t_m_node_diff_ksi[ff](i) = -t_node_diff_ksi[ff](i);
786  }
787 
788 
789  FTensor::Tensor1<double *, 3> t_phi_v(&phi_v[HVEC0], &phi_v[HVEC1],
790  &phi_v[HVEC2], 3);
791  FTensor::Tensor2<double *, 3, 3> t_diff_phi_v(
792  &diff_phi_v[HVEC0_0], &diff_phi_v[HVEC0_1], &diff_phi_v[HVEC0_2],
793  &diff_phi_v[HVEC1_0], &diff_phi_v[HVEC1_1], &diff_phi_v[HVEC1_2],
794  &diff_phi_v[HVEC2_0], &diff_phi_v[HVEC2_1], &diff_phi_v[HVEC2_2], 9);
795 
796  MatrixDouble fk(3, p + 1), diff_fk(3, 3 * p + 3);
797 
798  for (int ii = 0; ii != gdim; ii++) {
799  const int shift = 4 * ii;
800 
801  for (int ff = 0; ff != 3; ff++) {
802  const int fff = znf[ff];
803  const int iO = opposite_face_node[fff];
804  const double nO = N[shift + iO];
805  for (int pp = 1; pp <= p; pp++) {
806  CHKERR IntegratedJacobi_polynomials(
807  pp, 2 * pp + 2, nO, 1 - nO, &t_node_diff_ksi[iO](0),
808  &t_m_node_diff_ksi[iO](0), &fk(ff, 0), &diff_fk(ff, 0), 3);
809  }
810  }
811 
812  int jj = 0;
813  for (int oo = 2; oo <= p; oo++) {
814  for (int k = 1; k != oo; k++) {
815  int OO = oo - k;
816  if (OO >= 0) {
817  int s = NBFACETRI_DEMKOWICZ_HDIV(OO - 1);
818  // Note that we do faces 0,2,3, skipping 1. All the faces which have
819  // zero node in it.
820  int nb_dofs = NBFACETRI_DEMKOWICZ_HDIV(p_f[znf[0]]);
821  int sp[] = {ii * 3 * nb_dofs + 3 * s, ii * 3 * nb_dofs + 3 * s,
822  ii * 3 * nb_dofs + 3 * s};
823  FTensor::Tensor1<double *, 3> t_phi_f[] = {
824  FTensor::Tensor1<double *, 3>(&phi_f[znf[0]][sp[0] + HVEC0],
825  &phi_f[znf[0]][sp[0] + HVEC1],
826  &phi_f[znf[0]][sp[0] + HVEC2], 3),
827  FTensor::Tensor1<double *, 3>(&phi_f[znf[1]][sp[1] + HVEC0],
828  &phi_f[znf[1]][sp[1] + HVEC1],
829  &phi_f[znf[1]][sp[1] + HVEC2], 3),
830  FTensor::Tensor1<double *, 3>(&phi_f[znf[2]][sp[2] + HVEC0],
831  &phi_f[znf[2]][sp[2] + HVEC1],
832  &phi_f[znf[2]][sp[2] + HVEC2], 3)};
833  int sdp[] = {ii * 9 * nb_dofs + 9 * s, ii * 9 * nb_dofs + 9 * s,
834  ii * 9 * nb_dofs + 9 * s};
835  FTensor::Tensor2<double *, 3, 3> t_diff_phi_f[] = {
836  FTensor::Tensor2<double *, 3, 3>(
837  &diff_phi_f[znf[0]][sdp[0] + HVEC0_0],
838  &diff_phi_f[znf[0]][sdp[0] + HVEC0_1],
839  &diff_phi_f[znf[0]][sdp[0] + HVEC0_2],
840  &diff_phi_f[znf[0]][sdp[0] + HVEC1_0],
841  &diff_phi_f[znf[0]][sdp[0] + HVEC1_1],
842  &diff_phi_f[znf[0]][sdp[0] + HVEC1_2],
843  &diff_phi_f[znf[0]][sdp[0] + HVEC2_0],
844  &diff_phi_f[znf[0]][sdp[0] + HVEC2_1],
845  &diff_phi_f[znf[0]][sdp[0] + HVEC2_2], 9),
846  FTensor::Tensor2<double *, 3, 3>(
847  &diff_phi_f[znf[1]][sdp[1] + HVEC0_0],
848  &diff_phi_f[znf[1]][sdp[1] + HVEC0_1],
849  &diff_phi_f[znf[1]][sdp[1] + HVEC0_2],
850  &diff_phi_f[znf[1]][sdp[1] + HVEC1_0],
851  &diff_phi_f[znf[1]][sdp[1] + HVEC1_1],
852  &diff_phi_f[znf[1]][sdp[1] + HVEC1_2],
853  &diff_phi_f[znf[1]][sdp[1] + HVEC2_0],
854  &diff_phi_f[znf[1]][sdp[1] + HVEC2_1],
855  &diff_phi_f[znf[1]][sdp[1] + HVEC2_2], 9),
856  FTensor::Tensor2<double *, 3, 3>(
857  &diff_phi_f[znf[2]][sdp[2] + HVEC0_0],
858  &diff_phi_f[znf[2]][sdp[2] + HVEC0_1],
859  &diff_phi_f[znf[2]][sdp[2] + HVEC0_2],
860  &diff_phi_f[znf[2]][sdp[2] + HVEC1_0],
861  &diff_phi_f[znf[2]][sdp[2] + HVEC1_1],
862  &diff_phi_f[znf[2]][sdp[2] + HVEC1_2],
863  &diff_phi_f[znf[2]][sdp[2] + HVEC2_0],
864  &diff_phi_f[znf[2]][sdp[2] + HVEC2_1],
865  &diff_phi_f[znf[2]][sdp[2] + HVEC2_2], 9)};
866  for (int ij = s; ij != NBFACETRI_DEMKOWICZ_HDIV(OO); ij++) {
867  for (int ff = 0; ff != 3; ff++) {
868  FTensor::Tensor1<double, 3> t_diff_fk(diff_fk(ff, 0 * p + k - 1),
869  diff_fk(ff, 1 * p + k - 1),
870  diff_fk(ff, 2 * p + k - 1));
871  t_phi_v(i) = fk(ff, k - 1) * t_phi_f[ff](i);
872  t_diff_phi_v(i, j) = t_diff_fk(j) * t_phi_f[ff](i) +
873  fk(ff, k - 1) * t_diff_phi_f[ff](i, j);
874  ++t_phi_v;
875  ++t_diff_phi_v;
876  ++t_phi_f[ff];
877  ++t_diff_phi_f[ff];
878  ++jj;
879  }
880  }
881  }
882  }
883  }
884  if (jj != NBVOLUMETET_DEMKOWICZ_HDIV(p)) {
885  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
886  "wrong number of base functions "
887  "jj!=NBVOLUMETET_DEMKOWICZ_HDIV(p) "
888  "%d!=%d",
889  jj, NBVOLUMETET_DEMKOWICZ_HDIV(p));
890  }
891  }
892 
893  MoFEMFunctionReturn(0);
894 }
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
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:368
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:617
FTensor::Tensor1< double, 3 >
HVEC0_2
@ HVEC0_2
Definition: definitions.h:198
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:13
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:37
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
HVEC0_1
@ HVEC0_1
Definition: definitions.h:195
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:762
HVEC0_0
@ HVEC0_0
Definition: definitions.h:192
HVEC1_1
@ HVEC1_1
Definition: definitions.h:196
HVEC1
@ HVEC1
Definition: definitions.h:186
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
HVEC2_1
@ HVEC2_1
Definition: definitions.h:197
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:151
FTensor::Number< 0 >
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
HVEC2_2
@ HVEC2_2
Definition: definitions.h:200
HVEC1_2
@ HVEC1_2
Definition: definitions.h:199
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
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:160
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:484
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:137
HVEC2_0
@ HVEC2_0
Definition: definitions.h:194
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:280
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:440
HVEC2
@ HVEC2
Definition: definitions.h:186
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
HVEC1_0
@ HVEC1_0
Definition: definitions.h:193
HVEC0
@ HVEC0
Definition: definitions.h:186
Generated by Doxygen 1.8.17 and hosted at University of Glasgow