v0.14.0
Loading...
Searching...
No Matches
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
11using namespace MoFEM;
12
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
21 for (int ff = 0; ff < 4; ff++) {
22 if (diff_phi_f_e != NULL) {
24 &faces_nodes[3 * ff], p[ff], N, diffN, phi_f_e[ff], diff_phi_f_e[ff],
25 gdim, 4, base_polynomials);
27 } else {
29 &faces_nodes[3 * ff], p[ff], N, diffN, phi_f_e[ff], NULL, gdim, 4,
30 base_polynomials);
32 }
33 }
35}
36
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};
48
50 if (p < 1)
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<
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 }
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 }
135}
136
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
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 }
153 &faces_nodes[3 * ff], p[ff], N, diffN, phi_f[ff], diff, gdim, 4,
154 base_polynomials);
155 CHKERRG(ierr);
156 }
158}
159
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)) {
168
170 if (p < 3)
172
173 const int vert_i = face_nodes[1];
174 const int vert_j = face_nodes[2];
175 const int i0 = face_nodes[0];
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> >(
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 }
278}
279
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
291 if (p < 2)
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
310
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);
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]);
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
366}
367
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
378 if (p < 3)
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
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];
415 &phi_v_f[ff][HVEC0], &phi_v_f[ff][HVEC1], &phi_v_f[ff][HVEC2], 3);
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);
441 int jj = 0;
442 for (int oo = 0; oo <= p - 3; oo++) {
443 FTensor::Tensor0<double *> t_psi_l(&psi_l[0], 1);
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
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 }
475 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
476 "wrong order %d != %d", jj,
478 }
479 }
480 }
482}
483
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
492 if (p < 4)
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
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);
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
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]);
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,
610 }
611 }
612
614}
615
617MoFEM::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
624
627
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);
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> >(
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]];
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++) {
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) {
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",
756 }
757 }
758
760}
761
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};
770
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);
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++) {
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};
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[] = {
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),
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),
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",
890 }
891 }
892
894}
static Index< 'o', 3 > o
static Number< 2 > N2
static Index< 'p', 3 > p
static Number< 1 > N1
static Number< 0 > N0
constexpr double a
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.
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.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
@ HVEC0
Definition: definitions.h:186
@ HVEC1
Definition: definitions.h:186
@ HVEC2
Definition: definitions.h:186
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
@ HVEC1_1
Definition: definitions.h:196
@ HVEC0_1
Definition: definitions.h:195
@ HVEC1_0
Definition: definitions.h:193
@ HVEC2_1
Definition: definitions.h:197
@ HVEC1_2
Definition: definitions.h:199
@ HVEC2_2
Definition: definitions.h:200
@ HVEC2_0
Definition: definitions.h:194
@ HVEC0_2
Definition: definitions.h:198
@ HVEC0_0
Definition: definitions.h:192
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
const int dim
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
#define NBVOLUMETET_DEMKOWICZ_HDIV(P)
#define NBVOLUMETET_AINSWORTH_FACE_HDIV(P)
#define NBFACETRI_DEMKOWICZ_HDIV(P)
#define NBFACETRI_AINSWORTH_FACE_HDIV(P)
#define NBVOLUMETET_AINSWORTH_VOLUME_HDIV(P)
FTensor::Index< 'i', SPACE_DIM > i
static double lambda
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
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
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
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
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
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
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
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
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
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
const int N
Definition: speed_test.cpp:3