v0.14.0
testing_base_functions.cpp
Go to the documentation of this file.
1 
2 
3 #include <MoFEM.hpp>
4 #include <quad.h>
5 
6 using namespace MoFEM;
7 
8 static char help[] = "testing interface inserting algorithm\n\n";
9 
10 static double sum_matrix(MatrixDouble &m) {
11  double s = 0;
12  for (unsigned int ii = 0; ii < m.size1(); ii++) {
13  for (unsigned int jj = 0; jj < m.size2(); jj++) {
14  s += m(ii, jj);
15  }
16  }
17  return s;
18 }
19 
20 int main(int argc, char *argv[]) {
21 
22  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
23 
24  try {
25 
26  moab::Core mb_instance;
27  moab::Interface &moab = mb_instance;
28  MoFEM::Core core(moab);
29  // MoFEM::Interface& m_field = core;
30 
31  enum bases {
32  LEGENDREPOLYNOMIAL,
33  LOBATTOPOLYNOMIAL,
34  JACOBIPOLYNOMIAL,
35  INTEGRATEDJACOBIPOLYNOMIAL,
36  H1TET_AINSWORTH,
37  H1TET_BERNSTEIN_BEZIER,
38  HDIVTET_AINSWORTH,
39  HDIVTET_DEMKOWICZ,
40  HCURLTET_AINSWORTH,
41  HCURLTET_DEMKOWICZ,
42  L2TET,
43  H1TRI_AINSWORTH,
44  H1TRI_BERNSTEIN_BEZIER,
45  H1QUAD,
46  HDIVTRI_AINSWORTH,
47  HDIVTRI_DEMKOWICZ,
48  HCURLTRI_AINSWORTH,
49  HCURLTRI_DEMKOWICZ,
50  L2TRI,
51  H1EDGE_AINSWORTH,
52  H1EDGE_BERNSTEIN_BEZIER,
53  HCURLEDGE_AINSWORTH,
54  HCURLEDGE_DEMKOWICZ,
55  L2EDGE,
56  H1FLATPRIS,
57  H1FATPRISM,
58  LASTOP
59  };
60 
61  const char *list[] = {"legendre",
62  "lobatto",
63  "jacobi",
64  "integrated_jacobi",
65  "h1tet_ainsworth",
66  "h1tet_bernstein_bezier",
67  "hdivtet_ainsworth",
68  "hdivtet_demkowicz",
69  "hcurltet_ainsworth",
70  "hcurltet_demkowicz",
71  "l2tet",
72  "h1tri_ainsworth",
73  "h1tri_bernstein_bezier",
74  "h1quad",
75  "hdivtri_ainsworth",
76  "hdivtri_demkowicz",
77  "hcurltri_ainsworth",
78  "hcurltri_demkowicz",
79  "l2tri",
80  "h1edge_ainsworth",
81  "h1edge_bernstein_bezier",
82  "hcurledge_ainsworth",
83  "hcurledge_demkowicz",
84  "l2edge",
85  "h1flatprism",
86  "h1fatprism"};
87 
88  PetscBool flg;
89  PetscInt choice_value = LEGENDREPOLYNOMIAL;
90  ierr = PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list, LASTOP,
91  &choice_value, &flg);
92  CHKERRG(ierr);
93  if (flg != PETSC_TRUE) {
94  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "base not set");
95  }
96 
97  MatrixDouble pts_1d(1, 3);
98  pts_1d(0, 0) = -0.5;
99  pts_1d(0, 1) = 0.;
100  pts_1d(0, 2) = +0.5;
101 
102  boost::shared_ptr<MatrixDouble> base_ptr(new MatrixDouble);
103  boost::shared_ptr<MatrixDouble> diff_base_ptr(new MatrixDouble);
104 
105  const double eps = 1e-3;
106 
107  if (choice_value == LEGENDREPOLYNOMIAL) {
108  double diff_s = 1;
109  ierr = LegendrePolynomial().getValue(
110  pts_1d, boost::shared_ptr<BaseFunctionCtx>(new LegendrePolynomialCtx(
111  4, &diff_s, 1, base_ptr, diff_base_ptr)));
112  CHKERRG(ierr);
113 
114  std::cout << "LegendrePolynomial\n";
115  std::cout << pts_1d << std::endl;
116  std::cout << *base_ptr << std::endl;
117  std::cout << *diff_base_ptr << std::endl;
118  double sum = sum_matrix(*base_ptr);
119  double diff_sum = sum_matrix(*diff_base_ptr);
120  std::cout << sum << std::endl;
121  std::cout << diff_sum << std::endl;
122  if (fabs(2.04688 - sum) > eps) {
123  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
124  }
125  if (fabs(2.25 - diff_sum) > eps) {
126  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
127  }
128  }
129 
130  pts_1d.resize(1, 11, false);
131  for (int ii = 0; ii != 11; ii++) {
132  pts_1d(0, ii) = 2 * ((double)ii / 10) - 1;
133  }
134 
135  boost::shared_ptr<MatrixDouble> kernel_base_ptr(new MatrixDouble);
136  boost::shared_ptr<MatrixDouble> diff_kernel_base_ptr(new MatrixDouble);
137 
138  if (choice_value == LOBATTOPOLYNOMIAL) {
139  double diff_s = 1;
140  ierr = LobattoPolynomial().getValue(
141  pts_1d, boost::shared_ptr<BaseFunctionCtx>(new LobattoPolynomialCtx(
142  7 + 2, &diff_s, 1, base_ptr, diff_base_ptr)));
143  CHKERRG(ierr);
144  ierr = KernelLobattoPolynomial().getValue(
145  pts_1d,
146  boost::shared_ptr<BaseFunctionCtx>(new KernelLobattoPolynomialCtx(
147  7, &diff_s, 1, kernel_base_ptr, diff_kernel_base_ptr)));
148  CHKERRG(ierr);
149  for (int ii = 0; ii != 11; ii++) {
150  double s = pts_1d(0, ii);
151  std::cerr << "lobatto_plot " << s << " " << (*base_ptr)(ii, 1) << " "
152  << (*diff_base_ptr)(ii, 1) << " " << (*kernel_base_ptr)(ii, 1)
153  << " " << (*diff_kernel_base_ptr)(ii, 1) << " "
154  << (*kernel_base_ptr)(ii, 1) * (1 - s * s) << " "
155  << (*kernel_base_ptr)(ii, 1) * (-2 * s) +
156  (*diff_kernel_base_ptr)(ii, 1) * (1 - s * s)
157  << " " << std::endl;
158  }
159  std::cout << "LobattoPolynomial\n";
160  std::cout << pts_1d << std::endl;
161  std::cout << *base_ptr << std::endl;
162  std::cout << *diff_base_ptr << std::endl;
163  {
164  double sum = sum_matrix(*base_ptr);
165  double diff_sum = sum_matrix(*diff_base_ptr);
166  std::cout << sum << std::endl;
167  std::cout << diff_sum << std::endl;
168  if (fabs(9.39466 - sum) > eps) {
169  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
170  }
171  if (fabs(14.0774 - diff_sum) > eps) {
172  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
173  }
174  }
175  {
176  double sum = sum_matrix(*kernel_base_ptr);
177  double diff_sum = sum_matrix(*diff_kernel_base_ptr);
178  std::cout << sum << std::endl;
179  std::cout << diff_sum << std::endl;
180  if (fabs(-13.9906 * 4 - sum) > eps) {
181  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
182  }
183  if (fabs(-101.678 * 4 - diff_sum) > eps) {
184  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
185  }
186  }
187  }
188 
189  if (choice_value == JACOBIPOLYNOMIAL) {
190 
191  int n = 21;
192  MatrixDouble pts_1d(1, n);
193  pts_1d.resize(1, n, false);
194  MatrixDouble pts_1d_t(1, n);
195  for (int ii = 0; ii != n; ii++) {
196  pts_1d(0, ii) = (double)ii / (n - 1);
197  pts_1d_t(0, ii) = 1;
198  }
199 
200  base_ptr->clear();
201  diff_base_ptr->clear();
202 
203  double diff_x = 1;
204  double diff_t = 0;
205  ierr = JacobiPolynomial().getValue(
206  pts_1d, pts_1d_t,
207  boost::shared_ptr<BaseFunctionCtx>(new JacobiPolynomialCtx(
208  5, &diff_x, &diff_t, 1, 0, base_ptr, diff_base_ptr)));
209  CHKERRG(ierr);
210  for (int ii = 0; ii != n; ii++) {
211  double s = pts_1d(0, ii);
212  std::cerr << "jacobi_plot " << s << " " << (*base_ptr)(ii, 4) << " "
213  << (*diff_base_ptr)(ii, 4) << std::endl;
214  }
215  for (int ii = 0; ii != n - 1; ii++) {
216  double s = (pts_1d(0, ii) + pts_1d(0, ii + 1)) / 2;
217  std::cerr << "jacobi_diff_plot_fd " << s << " "
218  << ((*base_ptr)(ii + 1, 4) - (*base_ptr)(ii, 4)) /
219  (1. / (n - 1))
220  << std::endl;
221  }
222  std::cout << "JacobiPolynomial\n";
223  std::cout << pts_1d << std::endl;
224  std::cout << *base_ptr << std::endl;
225  std::cout << *diff_base_ptr << std::endl;
226  {
227  double sum = sum_matrix(*base_ptr);
228  double diff_sum = sum_matrix(*diff_base_ptr);
229  std::cout << sum << std::endl;
230  std::cout << diff_sum << std::endl;
231  if (fabs(23.2164 - sum) > eps) {
232  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
233  }
234  if (fabs(167.995 - diff_sum) > eps) {
235  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
236  }
237  }
238  }
239 
240  if (choice_value == INTEGRATEDJACOBIPOLYNOMIAL) {
241 
242  int n = 21;
243  MatrixDouble pts_1d(1, n);
244  pts_1d.resize(1, n, false);
245  MatrixDouble pts_1d_t(1, n);
246  for (int ii = 0; ii != n; ii++) {
247  pts_1d(0, ii) = (double)ii / 20.;
248  pts_1d_t(0, ii) = 1 - pts_1d(0, ii);
249  }
250 
251  base_ptr->clear();
252  diff_base_ptr->clear();
253 
254  double diff_x = 1;
255  double diff_t = 0;
256  ierr = IntegratedJacobiPolynomial().getValue(
257  pts_1d, pts_1d_t,
258  boost::shared_ptr<BaseFunctionCtx>(new IntegratedJacobiPolynomialCtx(
259  6, &diff_x, &diff_t, 1, 1, base_ptr, diff_base_ptr)));
260  CHKERRG(ierr);
261  for (int ii = 0; ii != n; ii++) {
262  double s = pts_1d(0, ii);
263  std::cerr << "integrated_jacobi_plot " << s << " " << (double)ii / 20.
264  << " " << (*base_ptr)(ii, 1) << " " << (*base_ptr)(ii, 2)
265  << " " << (*base_ptr)(ii, 3) << " " << (*base_ptr)(ii, 4)
266  << " " << (*base_ptr)(ii, 5) << endl;
267  ;
268  // << " " << (*diff_base_ptr)(ii, 4) << std::endl;
269  }
270  std::cout << "IntegratedJacobiPolynomial\n";
271  std::cout << pts_1d << std::endl;
272  std::cout << *base_ptr << std::endl;
273  std::cout << *diff_base_ptr << std::endl;
274  {
275  double sum = sum_matrix(*base_ptr);
276  double diff_sum = sum_matrix(*diff_base_ptr);
277  std::cout << sum << std::endl;
278  std::cout << diff_sum << std::endl;
279  // if(fabs(7.1915-sum)>eps) {
280  // SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"wrong result");
281  // }
282  // if(fabs(23.2164-diff_sum)>eps) {
283  // SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"wrong result");
284  // }
285  }
286  }
287 
288  EntitiesFieldData tet_data(MBTET);
289  for (int type = MBVERTEX; type != MBMAXTYPE; type++) {
290  tet_data.spacesOnEntities[type].set(L2);
291  tet_data.spacesOnEntities[type].set(H1);
292  tet_data.spacesOnEntities[type].set(HDIV);
293  tet_data.spacesOnEntities[type].set(HCURL);
294  }
295  tet_data.dataOnEntities[MBVERTEX].resize(1);
296  tet_data.dataOnEntities[MBVERTEX][0].getOrder() = 1;
297  tet_data.dataOnEntities[MBEDGE].resize(6);
298  for (int ee = 0; ee < 6; ee++) {
299  tet_data.dataOnEntities[MBEDGE][ee].getOrder() = 3;
300  tet_data.dataOnEntities[MBEDGE][ee].getSense() = 1;
301  }
302  tet_data.dataOnEntities[MBTRI].resize(4);
303  for (int ff = 0; ff < 4; ff++) {
304  tet_data.dataOnEntities[MBTRI][ff].getOrder() = 4;
305  tet_data.dataOnEntities[MBTRI][ff].getSense() = 1;
306  }
307  tet_data.dataOnEntities[MBTET].resize(1);
308  tet_data.dataOnEntities[MBTET][0].getOrder() = 5;
309  tet_data.dataOnEntities[MBVERTEX][0].getBBNodeOrder().resize(4, false);
310  std::fill(tet_data.dataOnEntities[MBVERTEX][0].getBBNodeOrder().begin(),
311  tet_data.dataOnEntities[MBVERTEX][0].getBBNodeOrder().end(), 5);
312 
313  MatrixDouble pts_tet;
314  int tet_rule = 2;
315  int nb_gauss_pts = QUAD_3D_TABLE[tet_rule]->npoints;
316  pts_tet.resize(3, nb_gauss_pts, false);
317  cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[tet_rule]->points[1], 4,
318  &pts_tet(0, 0), 1);
319  cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[tet_rule]->points[2], 4,
320  &pts_tet(1, 0), 1);
321  cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[tet_rule]->points[3], 4,
322  &pts_tet(2, 0), 1);
323  tet_data.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4,
324  false);
325  {
326  double *shape_ptr =
327  &*tet_data.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
328  cblas_dcopy(4 * nb_gauss_pts, QUAD_3D_TABLE[tet_rule]->points, 1,
329  shape_ptr, 1);
330  }
331  tet_data.facesNodes.resize(4, 3);
332  const int cannonical_tet_face[4][3] = {
333  {0, 1, 3}, {1, 2, 3}, {0, 3, 2}, {0, 2, 1}};
334  for (int ff = 0; ff < 4; ff++) {
335  for (int nn = 0; nn < 3; nn++) {
336  tet_data.facesNodes(ff, nn) = cannonical_tet_face[ff][nn];
337  }
338  }
339 
340  if (choice_value == H1TET_AINSWORTH) {
341 
342  tet_data.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4,
343  false);
344  ierr = ShapeMBTET(
345  &*tet_data.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin(),
346  &pts_tet(0, 0), &pts_tet(1, 0), &pts_tet(2, 0), nb_gauss_pts);
347  CHKERRG(ierr);
348 
349  ierr = TetPolynomialBase().getValue(
350  pts_tet,
351  boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
352  tet_data, H1, CONTINUOUS, AINSWORTH_LEGENDRE_BASE, NOBASE)));
353  CHKERRG(ierr);
354  if (tet_data.dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE).get() !=
355  tet_data.dataOnEntities[MBVERTEX][0]
356  .getNSharedPtr(AINSWORTH_LEGENDRE_BASE)
357  .get()) {
358  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
359  "Different pointers");
360  }
361 
362  double sum = 0, diff_sum = 0;
363  std::cout << "Edges\n";
364  for (int ee = 0; ee < 6; ee++) {
365  std::cout << tet_data.dataOnEntities[MBEDGE][ee].getN(
366  AINSWORTH_LEGENDRE_BASE)
367  << std::endl;
368  std::cout << tet_data.dataOnEntities[MBEDGE][ee].getDiffN(
369  AINSWORTH_LEGENDRE_BASE)
370  << std::endl;
371  sum += sum_matrix(
372  tet_data.dataOnEntities[MBEDGE][ee].getN(AINSWORTH_LEGENDRE_BASE));
373  diff_sum += sum_matrix(tet_data.dataOnEntities[MBEDGE][ee].getDiffN(
374  AINSWORTH_LEGENDRE_BASE));
375  }
376  std::cout << "Faces\n";
377  for (int ff = 0; ff < 4; ff++) {
378  std::cout << tet_data.dataOnEntities[MBTRI][ff].getN(
379  AINSWORTH_LEGENDRE_BASE)
380  << std::endl;
381  std::cout << tet_data.dataOnEntities[MBTRI][ff].getDiffN(
382  AINSWORTH_LEGENDRE_BASE)
383  << std::endl;
384  sum += sum_matrix(
385  tet_data.dataOnEntities[MBTRI][ff].getN(AINSWORTH_LEGENDRE_BASE));
386  diff_sum += sum_matrix(tet_data.dataOnEntities[MBTRI][ff].getDiffN(
387  AINSWORTH_LEGENDRE_BASE));
388  }
389  std::cout << "Tets\n";
390  std::cout << tet_data.dataOnEntities[MBTET][0].getN(
391  AINSWORTH_LEGENDRE_BASE)
392  << std::endl;
393  std::cout << tet_data.dataOnEntities[MBTET][0].getDiffN(
394  AINSWORTH_LEGENDRE_BASE)
395  << std::endl;
396  sum += sum_matrix(
397  tet_data.dataOnEntities[MBTET][0].getN(AINSWORTH_LEGENDRE_BASE));
398  diff_sum += sum_matrix(
399  tet_data.dataOnEntities[MBTET][0].getDiffN(AINSWORTH_LEGENDRE_BASE));
400  std::cout << "sum " << sum << std::endl;
401  std::cout << "diff_sum " << diff_sum << std::endl;
402  if (fabs(1.3509 - sum) > eps) {
403  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
404  }
405  if (fabs(0.233313 - diff_sum) > eps) {
406  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
407  }
408  }
409 
410  if (choice_value == H1TET_BERNSTEIN_BEZIER) {
411 
412  tet_data.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4,
413  false);
414  CHKERR Tools::shapeFunMBTET(
415  &*tet_data.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin(),
416  &pts_tet(0, 0), &pts_tet(1, 0), &pts_tet(2, 0), nb_gauss_pts);
417  CHKERR TetPolynomialBase().getValue(
418  pts_tet, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
419  tet_data, "TEST_FIELD", H1, CONTINUOUS,
420  AINSWORTH_BERNSTEIN_BEZIER_BASE, NOBASE)));
421  if (tet_data.dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE).get() ==
422  tet_data.dataOnEntities[MBVERTEX][0]
423  .getNSharedPtr(AINSWORTH_BERNSTEIN_BEZIER_BASE)
424  .get())
425  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "The same pointers");
426 
427  double sum = 0, diff_sum = 0;
428  std::cout << "Vertices\n";
429  std::cout << tet_data.dataOnEntities[MBVERTEX][0].getN("TEST_FIELD")
430  << std::endl;
431  std::cout << tet_data.dataOnEntities[MBVERTEX][0].getDiffN("TEST_FIELD")
432  << std::endl;
433  sum +=
434  sum_matrix(tet_data.dataOnEntities[MBVERTEX][0].getN("TEST_FIELD"));
435  diff_sum += sum_matrix(
436  tet_data.dataOnEntities[MBVERTEX][0].getDiffN("TEST_FIELD"));
437  std::cout << "Edges\n";
438  for (int ee = 0; ee < 6; ee++) {
439  std::cout << tet_data.dataOnEntities[MBEDGE][ee].getN("TEST_FIELD")
440  << std::endl;
441  std::cout << tet_data.dataOnEntities[MBEDGE][ee].getDiffN("TEST_FIELD")
442  << std::endl;
443  sum +=
444  sum_matrix(tet_data.dataOnEntities[MBEDGE][ee].getN("TEST_FIELD"));
445  diff_sum += sum_matrix(
446  tet_data.dataOnEntities[MBEDGE][ee].getDiffN("TEST_FIELD"));
447  }
448  std::cout << "Faces\n";
449  for (int ff = 0; ff < 4; ff++) {
450  std::cout << tet_data.dataOnEntities[MBTRI][ff].getN("TEST_FIELD")
451  << std::endl;
452  std::cout << tet_data.dataOnEntities[MBTRI][ff].getDiffN("TEST_FIELD")
453  << std::endl;
454  sum +=
455  sum_matrix(tet_data.dataOnEntities[MBTRI][ff].getN("TEST_FIELD"));
456  diff_sum += sum_matrix(
457  tet_data.dataOnEntities[MBTRI][ff].getDiffN("TEST_FIELD"));
458  }
459  std::cout << "Tets\n";
460  std::cout << tet_data.dataOnEntities[MBTET][0].getN("TEST_FIELD")
461  << std::endl;
462  std::cout << tet_data.dataOnEntities[MBTET][0].getDiffN("TEST_FIELD")
463  << std::endl;
464  sum += sum_matrix(tet_data.dataOnEntities[MBTET][0].getN("TEST_FIELD"));
465  diff_sum +=
466  sum_matrix(tet_data.dataOnEntities[MBTET][0].getDiffN("TEST_FIELD"));
467  std::cout << "sum " << sum << std::endl;
468  std::cout << "diff_sum " << diff_sum << std::endl;
469  if (fabs(4.38395 - sum) > eps) {
470  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
471  }
472  if (fabs(diff_sum) > eps) {
473  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
474  }
475  }
476 
477  if (choice_value == HDIVTET_AINSWORTH) {
478  ierr = TetPolynomialBase().getValue(
479  pts_tet, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
480  tet_data, HDIV, CONTINUOUS, AINSWORTH_LEGENDRE_BASE)));
481  CHKERRG(ierr);
482  double sum = 0, diff_sum = 0;
483  std::cout << "Faces\n";
484  for (int ff = 0; ff < 4; ff++) {
485  std::cout << tet_data.dataOnEntities[MBTRI][ff].getN(
486  AINSWORTH_LEGENDRE_BASE)
487  << std::endl;
488  std::cout << tet_data.dataOnEntities[MBTRI][ff].getDiffN(
489  AINSWORTH_LEGENDRE_BASE)
490  << std::endl;
491  sum += sum_matrix(
492  tet_data.dataOnEntities[MBTRI][ff].getN(AINSWORTH_LEGENDRE_BASE));
493  diff_sum += sum_matrix(tet_data.dataOnEntities[MBTRI][ff].getDiffN(
494  AINSWORTH_LEGENDRE_BASE));
495  }
496  std::cout << "Tets\n";
497  std::cout << tet_data.dataOnEntities[MBTET][0].getN(
498  AINSWORTH_LEGENDRE_BASE)
499  << std::endl;
500  std::cout << tet_data.dataOnEntities[MBTET][0].getDiffN(
501  AINSWORTH_LEGENDRE_BASE)
502  << std::endl;
503  sum += sum_matrix(
504  tet_data.dataOnEntities[MBTET][0].getN(AINSWORTH_LEGENDRE_BASE));
505  diff_sum += sum_matrix(
506  tet_data.dataOnEntities[MBTET][0].getDiffN(AINSWORTH_LEGENDRE_BASE));
507  std::cout << "sum " << 1e8 * sum << std::endl;
508  std::cout << "diff_sum " << 1e8 * diff_sum << std::endl;
509  if (fabs(0.188636 - sum) > eps) {
510  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
511  }
512  if (fabs(32.9562 - diff_sum) > eps) {
513  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
514  }
515  }
516 
517  if (choice_value == HDIVTET_DEMKOWICZ) {
518  ierr = TetPolynomialBase().getValue(
519  pts_tet, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
520  tet_data, HDIV, CONTINUOUS, DEMKOWICZ_JACOBI_BASE)));
521  CHKERRG(ierr);
522  double sum = 0, diff_sum = 0;
523  std::cout << "Faces\n";
524  for (int ff = 0; ff < 4; ff++) {
525  std::cout << tet_data.dataOnEntities[MBTRI][ff].getN(
526  DEMKOWICZ_JACOBI_BASE)
527  << std::endl;
528  std::cout << tet_data.dataOnEntities[MBTRI][ff].getDiffN(
529  DEMKOWICZ_JACOBI_BASE)
530  << std::endl;
531  sum += sum_matrix(
532  tet_data.dataOnEntities[MBTRI][ff].getN(DEMKOWICZ_JACOBI_BASE));
533  diff_sum += sum_matrix(
534  tet_data.dataOnEntities[MBTRI][ff].getDiffN(DEMKOWICZ_JACOBI_BASE));
535  }
536  std::cout << "Tets\n";
537  std::cout << tet_data.dataOnEntities[MBTET][0].getN(DEMKOWICZ_JACOBI_BASE)
538  << std::endl;
539  std::cout << tet_data.dataOnEntities[MBTET][0].getDiffN(
540  DEMKOWICZ_JACOBI_BASE)
541  << std::endl;
542  sum += sum_matrix(
543  tet_data.dataOnEntities[MBTET][0].getN(DEMKOWICZ_JACOBI_BASE));
544  diff_sum += sum_matrix(
545  tet_data.dataOnEntities[MBTET][0].getDiffN(DEMKOWICZ_JACOBI_BASE));
546  std::cout << "sum " << sum << std::endl;
547  std::cout << "diff_sum " << diff_sum << std::endl;
548  const double expected_result = -2.70651;
549  const double expected_diff_result = 289.421;
550  if (fabs((expected_result - sum) / expected_result) > eps) {
551  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
552  }
553  if (fabs((expected_diff_result - diff_sum) / expected_diff_result) >
554  eps) {
555  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
556  }
557  }
558 
559  if (choice_value == HCURLTET_AINSWORTH) {
560  ierr = TetPolynomialBase().getValue(
561  pts_tet, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
562  tet_data, HCURL, CONTINUOUS, AINSWORTH_LEGENDRE_BASE)));
563  CHKERRG(ierr);
564  double sum = 0, diff_sum = 0;
565  std::cout << "Edges\n";
566  for (int ee = 0; ee < 6; ee++) {
567  std::cout << tet_data.dataOnEntities[MBEDGE][ee].getN(
568  AINSWORTH_LEGENDRE_BASE)
569  << std::endl;
570  std::cout << tet_data.dataOnEntities[MBEDGE][ee].getDiffN(
571  AINSWORTH_LEGENDRE_BASE)
572  << std::endl;
573  sum += sum_matrix(
574  tet_data.dataOnEntities[MBEDGE][ee].getN(AINSWORTH_LEGENDRE_BASE));
575  diff_sum += sum_matrix(tet_data.dataOnEntities[MBEDGE][ee].getDiffN(
576  AINSWORTH_LEGENDRE_BASE));
577  }
578  std::cout << "Faces\n";
579  for (int ff = 0; ff < 4; ff++) {
580  std::cout << tet_data.dataOnEntities[MBTRI][ff].getN(
581  AINSWORTH_LEGENDRE_BASE)
582  << std::endl;
583  std::cout << tet_data.dataOnEntities[MBTRI][ff].getDiffN(
584  AINSWORTH_LEGENDRE_BASE)
585  << std::endl;
586  sum += sum_matrix(
587  tet_data.dataOnEntities[MBTRI][ff].getN(AINSWORTH_LEGENDRE_BASE));
588  diff_sum += sum_matrix(tet_data.dataOnEntities[MBTRI][ff].getDiffN(
589  AINSWORTH_LEGENDRE_BASE));
590  }
591  std::cout << "Tets\n";
592  std::cout << tet_data.dataOnEntities[MBTET][0].getN(
593  AINSWORTH_LEGENDRE_BASE)
594  << std::endl;
595  std::cout << tet_data.dataOnEntities[MBTET][0].getDiffN(
596  AINSWORTH_LEGENDRE_BASE)
597  << std::endl;
598  sum += sum_matrix(
599  tet_data.dataOnEntities[MBTET][0].getN(AINSWORTH_LEGENDRE_BASE));
600  diff_sum += sum_matrix(
601  tet_data.dataOnEntities[MBTET][0].getDiffN(AINSWORTH_LEGENDRE_BASE));
602  std::cout << "sum " << sum << std::endl;
603  std::cout << "diff_sum " << diff_sum << std::endl;
604  if (fabs(-1.7798 - sum) > eps) {
605  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
606  }
607  if (fabs(-67.1793 - diff_sum) > eps) {
608  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
609  }
610  }
611 
612  if (choice_value == HCURLTET_DEMKOWICZ) {
613  CHKERR TetPolynomialBase().getValue(
614  pts_tet, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
615  tet_data, HCURL, CONTINUOUS, DEMKOWICZ_JACOBI_BASE)));
616  double sum = 0, diff_sum = 0;
617  std::cout << "Edges\n";
618  for (int ee = 0; ee < 6; ee++) {
619  std::cout << tet_data.dataOnEntities[MBEDGE][ee].getN(
620  DEMKOWICZ_JACOBI_BASE)
621  << std::endl;
622  std::cout << tet_data.dataOnEntities[MBEDGE][ee].getDiffN(
623  DEMKOWICZ_JACOBI_BASE)
624  << std::endl;
625  sum += sum_matrix(
626  tet_data.dataOnEntities[MBEDGE][ee].getN(DEMKOWICZ_JACOBI_BASE));
627  diff_sum += sum_matrix(tet_data.dataOnEntities[MBEDGE][ee].getDiffN(
628  DEMKOWICZ_JACOBI_BASE));
629  }
630  std::cout << "Faces\n";
631  for (int ff = 0; ff < 4; ff++) {
632  std::cout << tet_data.dataOnEntities[MBTRI][ff].getN(
633  DEMKOWICZ_JACOBI_BASE)
634  << std::endl;
635  std::cout << tet_data.dataOnEntities[MBTRI][ff].getDiffN(
636  DEMKOWICZ_JACOBI_BASE)
637  << std::endl;
638  sum += sum_matrix(
639  tet_data.dataOnEntities[MBTRI][ff].getN(DEMKOWICZ_JACOBI_BASE));
640  diff_sum += sum_matrix(
641  tet_data.dataOnEntities[MBTRI][ff].getDiffN(DEMKOWICZ_JACOBI_BASE));
642  }
643  std::cout << "Tets\n";
644  std::cout << tet_data.dataOnEntities[MBTET][0].getN(DEMKOWICZ_JACOBI_BASE)
645  << std::endl;
646  std::cout << tet_data.dataOnEntities[MBTET][0].getDiffN(
647  DEMKOWICZ_JACOBI_BASE)
648  << std::endl;
649  sum += sum_matrix(
650  tet_data.dataOnEntities[MBTET][0].getN(DEMKOWICZ_JACOBI_BASE));
651  diff_sum += sum_matrix(
652  tet_data.dataOnEntities[MBTET][0].getDiffN(DEMKOWICZ_JACOBI_BASE));
653  std::cout << "sum " << sum << std::endl;
654  std::cout << "diff_sum " << diff_sum << std::endl;
655  if (fabs(7.35513 - sum) > eps) {
656  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
657  }
658  if (fabs(62.4549 - diff_sum) > eps) {
659  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
660  }
661  }
662 
663  if (choice_value == L2TET) {
664  ierr = TetPolynomialBase().getValue(
665  pts_tet, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
666  tet_data, L2, CONTINUOUS, AINSWORTH_LEGENDRE_BASE)));
667  CHKERRG(ierr);
668  double sum = 0, diff_sum = 0;
669  std::cout << "Tets\n";
670  std::cout << tet_data.dataOnEntities[MBTET][0].getN(
671  AINSWORTH_LEGENDRE_BASE)
672  << std::endl;
673  std::cout << tet_data.dataOnEntities[MBTET][0].getDiffN(
674  AINSWORTH_LEGENDRE_BASE)
675  << std::endl;
676  sum += sum_matrix(
677  tet_data.dataOnEntities[MBTET][0].getN(AINSWORTH_LEGENDRE_BASE));
678  diff_sum += sum_matrix(
679  tet_data.dataOnEntities[MBTET][0].getDiffN(AINSWORTH_LEGENDRE_BASE));
680  std::cout << "sum " << sum << std::endl;
681  std::cout << "diff_sum " << diff_sum << std::endl;
682  if (fabs(3.60352 - sum) > eps) {
683  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
684  }
685  if (fabs(-36.9994 - diff_sum) > eps) {
686  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
687  }
688  }
689 
690  EntitiesFieldData tri_data(MBTRI);
691  for (int type = MBVERTEX; type != MBMAXTYPE; type++) {
692  tri_data.spacesOnEntities[type].set(L2);
693  tri_data.spacesOnEntities[type].set(H1);
694  tri_data.spacesOnEntities[type].set(HDIV);
695  tri_data.spacesOnEntities[type].set(HCURL);
696  }
697  tri_data.dataOnEntities[MBVERTEX].resize(1);
698  tri_data.dataOnEntities[MBVERTEX][0].getOrder() = 1;
699  tri_data.dataOnEntities[MBEDGE].resize(3);
700  for (int ee = 0; ee < 3; ee++) {
701  tri_data.dataOnEntities[MBEDGE][ee].getOrder() = 3;
702  tri_data.dataOnEntities[MBEDGE][ee].getSense() = 1;
703  }
704  tri_data.dataOnEntities[MBTRI].resize(1);
705  tri_data.dataOnEntities[MBTRI][0].getOrder() = 4;
706  tri_data.dataOnEntities[MBVERTEX][0].getBBNodeOrder().resize(3, false);
707  std::fill(tri_data.dataOnEntities[MBVERTEX][0].getBBNodeOrder().begin(),
708  tri_data.dataOnEntities[MBVERTEX][0].getBBNodeOrder().end(), 4);
709 
710  MatrixDouble pts_tri;
711  int tri_rule = 2;
712  nb_gauss_pts = QUAD_2D_TABLE[tri_rule]->npoints;
713  pts_tri.resize(2, nb_gauss_pts, false);
714  cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[tri_rule]->points[1], 3,
715  &pts_tri(0, 0), 1);
716  cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[tri_rule]->points[2], 3,
717  &pts_tri(1, 0), 1);
718  tri_data.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 3,
719  false);
720  {
721  double *shape_ptr =
722  &*tri_data.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
723  cblas_dcopy(3 * nb_gauss_pts, QUAD_2D_TABLE[tri_rule]->points, 1,
724  shape_ptr, 1);
725  }
726  tri_data.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).resize(3, 2, false);
727  ierr = ShapeDiffMBTRI(
728  &*tri_data.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).data().begin());
729  CHKERRG(ierr);
730 
731  if (choice_value == H1TRI_AINSWORTH) {
732  CHKERR TriPolynomialBase().getValue(
733  pts_tri, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
734  tri_data, H1, CONTINUOUS, AINSWORTH_LEGENDRE_BASE, NOBASE)));
735  if (tri_data.dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE).get() !=
736  tri_data.dataOnEntities[MBVERTEX][0]
737  .getNSharedPtr(AINSWORTH_LEGENDRE_BASE)
738  .get())
739  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
740  "Different pointers");
741 
742  double sum = 0, diff_sum = 0;
743  std::cout << "Edges\n";
744  for (int ee = 0; ee < 3; ee++) {
745  std::cout << tri_data.dataOnEntities[MBEDGE][ee].getN(
746  AINSWORTH_LEGENDRE_BASE)
747  << std::endl;
748  std::cout << tri_data.dataOnEntities[MBEDGE][ee].getDiffN(
749  AINSWORTH_LEGENDRE_BASE)
750  << std::endl;
751  sum += sum_matrix(
752  tri_data.dataOnEntities[MBEDGE][ee].getN(AINSWORTH_LEGENDRE_BASE));
753  diff_sum += sum_matrix(tri_data.dataOnEntities[MBEDGE][ee].getDiffN(
754  AINSWORTH_LEGENDRE_BASE));
755  }
756  std::cout << "Face\n";
757  std::cout << tri_data.dataOnEntities[MBTRI][0].getN(
758  AINSWORTH_LEGENDRE_BASE)
759  << std::endl;
760  std::cout << tri_data.dataOnEntities[MBTRI][0].getDiffN(
761  AINSWORTH_LEGENDRE_BASE)
762  << std::endl;
763  sum += sum_matrix(
764  tri_data.dataOnEntities[MBTRI][0].getN(AINSWORTH_LEGENDRE_BASE));
765  diff_sum += sum_matrix(
766  tri_data.dataOnEntities[MBTRI][0].getDiffN(AINSWORTH_LEGENDRE_BASE));
767  std::cout << "sum " << sum << std::endl;
768  std::cout << "diff_sum " << diff_sum << std::endl;
769  if (fabs(0.805556 - sum) > eps)
770  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
771 
772  if (fabs(0.0833333 - diff_sum) > eps)
773  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
774  }
775 
776  if (choice_value == H1TRI_BERNSTEIN_BEZIER) {
777  CHKERR TriPolynomialBase().getValue(
778  pts_tri, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
779  tri_data, "TET_FIELD", H1, CONTINUOUS,
780  AINSWORTH_BERNSTEIN_BEZIER_BASE, NOBASE)));
781  if (tri_data.dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE).get() ==
782  tri_data.dataOnEntities[MBVERTEX][0]
783  .getNSharedPtr(AINSWORTH_BERNSTEIN_BEZIER_BASE)
784  .get())
785  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
786  "Pointers should be diffrent");
787 
788  double sum = 0, diff_sum = 0;
789  std::cout << "Vertex\n";
790  std::cout << tri_data.dataOnEntities[MBVERTEX][0].getN("TET_FIELD")
791  << std::endl;
792  std::cout << tri_data.dataOnEntities[MBVERTEX][0].getDiffN("TET_FIELD")
793  << std::endl;
794  sum += sum_matrix(tri_data.dataOnEntities[MBVERTEX][0].getN("TET_FIELD"));
795  diff_sum += sum_matrix(
796  tri_data.dataOnEntities[MBVERTEX][0].getDiffN("TET_FIELD"));
797  std::cout << "Edges\n";
798  for (int ee = 0; ee < 3; ee++) {
799  std::cout << tri_data.dataOnEntities[MBEDGE][ee].getN("TET_FIELD")
800  << std::endl;
801  std::cout << tri_data.dataOnEntities[MBEDGE][ee].getDiffN("TET_FIELD")
802  << std::endl;
803  sum +=
804  sum_matrix(tri_data.dataOnEntities[MBEDGE][ee].getN("TET_FIELD"));
805  diff_sum += sum_matrix(
806  tri_data.dataOnEntities[MBEDGE][ee].getDiffN("TET_FIELD"));
807  }
808  std::cout << "Face\n";
809  std::cout << tri_data.dataOnEntities[MBTRI][0].getN("TET_FIELD")
810  << std::endl;
811  std::cout << tri_data.dataOnEntities[MBTRI][0].getDiffN("TET_FIELD")
812  << std::endl;
813  sum += sum_matrix(tri_data.dataOnEntities[MBTRI][0].getN("TET_FIELD"));
814  diff_sum +=
815  sum_matrix(tri_data.dataOnEntities[MBTRI][0].getDiffN("TET_FIELD"));
816  std::cout << "sum " << sum << std::endl;
817  std::cout << "diff_sum " << diff_sum << std::endl;
818  if (std::abs(3.01389 - sum) > eps)
819  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
820 
821  if (std::abs(diff_sum) > eps)
822  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
823  "wrong result %3.4e != $3.4e", 0, diff_sum);
824  }
825 
826  if (choice_value == HDIVTRI_AINSWORTH) {
827  CHKERR TriPolynomialBase().getValue(
828  pts_tri, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
829  tri_data, HDIV, CONTINUOUS, AINSWORTH_LEGENDRE_BASE, NOBASE)));
830  if (tri_data.dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE).get() !=
831  tri_data.dataOnEntities[MBVERTEX][0]
832  .getNSharedPtr(AINSWORTH_LEGENDRE_BASE)
833  .get())
834  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
835  "Different pointers");
836 
837  double sum = 0;
838  std::cout << "Face\n";
839  std::cout << tri_data.dataOnEntities[MBTRI][0].getN(
840  AINSWORTH_LEGENDRE_BASE)
841  << std::endl;
842  sum += sum_matrix(
843  tri_data.dataOnEntities[MBTRI][0].getN(AINSWORTH_LEGENDRE_BASE));
844  std::cout << "sum " << sum << std::endl;
845  if (fabs(1.93056 - sum) > eps)
846  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
847  }
848 
849  if (choice_value == HDIVTRI_DEMKOWICZ) {
850  ierr = TriPolynomialBase().getValue(
851  pts_tri, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
852  tri_data, HDIV, CONTINUOUS, DEMKOWICZ_JACOBI_BASE, NOBASE)));
853  CHKERRG(ierr);
854  if (tri_data.dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE).get() !=
855  tri_data.dataOnEntities[MBVERTEX][0]
856  .getNSharedPtr(DEMKOWICZ_JACOBI_BASE)
857  .get()) {
858  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
859  "Different pointers");
860  }
861  double sum = 0;
862  std::cout << "Face\n";
863  std::cout << tri_data.dataOnEntities[MBTRI][0].getN(DEMKOWICZ_JACOBI_BASE)
864  << std::endl;
865  sum += sum_matrix(
866  tri_data.dataOnEntities[MBTRI][0].getN(DEMKOWICZ_JACOBI_BASE));
867  std::cout << "sum " << sum << std::endl;
868  const double expected_result = 28.25;
869  if (fabs((expected_result - sum) / expected_result) > eps) {
870  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
871  }
872  }
873 
874  if (choice_value == HCURLTRI_AINSWORTH) {
875  ierr = TriPolynomialBase().getValue(
876  pts_tri, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
877  tri_data, HCURL, CONTINUOUS, AINSWORTH_LEGENDRE_BASE, NOBASE)));
878  CHKERRG(ierr);
879  if (tri_data.dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE).get() !=
880  tri_data.dataOnEntities[MBVERTEX][0]
881  .getNSharedPtr(AINSWORTH_LEGENDRE_BASE)
882  .get()) {
883  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
884  "Different pointers");
885  }
886  double sum = 0; //,diff_sum = 0;
887  std::cout << "Edges\n";
888  for (int ee = 0; ee < 3; ee++) {
889  std::cout << tri_data.dataOnEntities[MBEDGE][ee].getN(
890  AINSWORTH_LEGENDRE_BASE)
891  << std::endl;
892  sum += sum_matrix(
893  tri_data.dataOnEntities[MBEDGE][ee].getN(AINSWORTH_LEGENDRE_BASE));
894  }
895  std::cout << "Face\n";
896  std::cout << tri_data.dataOnEntities[MBTRI][0].getN(
897  AINSWORTH_LEGENDRE_BASE)
898  << std::endl;
899  sum += sum_matrix(
900  tri_data.dataOnEntities[MBTRI][0].getN(AINSWORTH_LEGENDRE_BASE));
901  std::cout << "sum " << sum << std::endl;
902  if (fabs(0.333333 - sum) > eps) {
903  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
904  }
905  }
906 
907  if (choice_value == HCURLTRI_DEMKOWICZ) {
908  CHKERR TriPolynomialBase().getValue(
909  pts_tri, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
910  tri_data, HCURL, CONTINUOUS, DEMKOWICZ_JACOBI_BASE, NOBASE)));
911  if (tri_data.dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE).get() !=
912  tri_data.dataOnEntities[MBVERTEX][0]
913  .getNSharedPtr(DEMKOWICZ_JACOBI_BASE)
914  .get()) {
915  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
916  "Different pointers");
917  }
918  double sum = 0; //,diff_sum = 0;
919  std::cout << "Edges\n";
920  for (int ee = 0; ee < 3; ee++) {
921  std::cout << tri_data.dataOnEntities[MBEDGE][ee].getN(
922  DEMKOWICZ_JACOBI_BASE)
923  << std::endl;
924  sum += sum_matrix(
925  tri_data.dataOnEntities[MBEDGE][ee].getN(DEMKOWICZ_JACOBI_BASE));
926  }
927  std::cout << "Face\n";
928  std::cout << tri_data.dataOnEntities[MBTRI][0].getN(DEMKOWICZ_JACOBI_BASE)
929  << std::endl;
930  sum += sum_matrix(
931  tri_data.dataOnEntities[MBTRI][0].getN(DEMKOWICZ_JACOBI_BASE));
932  std::cout << "sum " << sum << std::endl;
933  if (fabs(1.04591 - sum) > eps) {
934  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
935  }
936  }
937 
938  if (choice_value == L2TRI) {
939  ierr = TriPolynomialBase().getValue(
940  pts_tri, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
941  tri_data, L2, CONTINUOUS, AINSWORTH_LEGENDRE_BASE, NOBASE)));
942  CHKERRG(ierr);
943  if (tri_data.dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE).get() !=
944  tri_data.dataOnEntities[MBVERTEX][0]
945  .getNSharedPtr(AINSWORTH_LEGENDRE_BASE)
946  .get()) {
947  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
948  "Different pointers");
949  }
950  double sum = 0;
951  std::cout << "Face\n";
952  std::cout << tri_data.dataOnEntities[MBTRI][0].getN(
953  AINSWORTH_LEGENDRE_BASE)
954  << std::endl;
955  sum += sum_matrix(
956  tri_data.dataOnEntities[MBTRI][0].getN(AINSWORTH_LEGENDRE_BASE));
957  std::cout << "sum " << sum << std::endl;
958  if (fabs(0.671875 - sum) > eps) {
959  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
960  }
961  }
962 
963  EntitiesFieldData edge_data(MBTRI);
964  for (int type = MBVERTEX; type != MBMAXTYPE; type++) {
965  edge_data.spacesOnEntities[type].set(L2);
966  edge_data.spacesOnEntities[type].set(H1);
967  edge_data.spacesOnEntities[type].set(HDIV);
968  edge_data.spacesOnEntities[type].set(HCURL);
969  }
970  edge_data.dataOnEntities[MBVERTEX].resize(1);
971  edge_data.dataOnEntities[MBVERTEX][0].getOrder() = 1;
972  edge_data.dataOnEntities[MBEDGE].resize(1);
973  edge_data.dataOnEntities[MBEDGE][0].getOrder() = 4;
974  edge_data.dataOnEntities[MBEDGE][0].getSense() = 1;
975  edge_data.dataOnEntities[MBVERTEX][0].getBBNodeOrder().resize(2, false);
976  std::fill(edge_data.dataOnEntities[MBVERTEX][0].getBBNodeOrder().begin(),
977  edge_data.dataOnEntities[MBVERTEX][0].getBBNodeOrder().end(), 4);
978 
979  MatrixDouble pts_edge;
980  int edge_rule = 6;
981  nb_gauss_pts = QUAD_1D_TABLE[edge_rule]->npoints;
982  pts_edge.resize(2, nb_gauss_pts, false);
983  cblas_dcopy(nb_gauss_pts, &QUAD_1D_TABLE[edge_rule]->points[1], 2,
984  &pts_edge(0, 0), 1);
985  cblas_dcopy(nb_gauss_pts, QUAD_1D_TABLE[edge_rule]->weights, 1,
986  &pts_edge(1, 0), 1);
987 
988  edge_data.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 2,
989  false);
990  {
991  double *shape_ptr =
992  &*edge_data.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
993  cblas_dcopy(2 * nb_gauss_pts, QUAD_1D_TABLE[edge_rule]->points, 1,
994  shape_ptr, 1);
995  }
996 
997  if (choice_value == H1EDGE_AINSWORTH) {
998  CHKERR EdgePolynomialBase().getValue(
999  pts_edge, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
1000  edge_data, H1, CONTINUOUS, AINSWORTH_LEGENDRE_BASE, NOBASE)));
1001  if (edge_data.dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE).get() !=
1002  edge_data.dataOnEntities[MBVERTEX][0]
1003  .getNSharedPtr(AINSWORTH_LEGENDRE_BASE)
1004  .get()) {
1005  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1006  "Different pointers");
1007  }
1008  double sum = 0, diff_sum = 0;
1009  std::cout << "Edge\n";
1010  std::cout << edge_data.dataOnEntities[MBEDGE][0].getN(
1011  AINSWORTH_LEGENDRE_BASE)
1012  << std::endl;
1013  std::cout << edge_data.dataOnEntities[MBEDGE][0].getDiffN(
1014  AINSWORTH_LEGENDRE_BASE)
1015  << std::endl;
1016  sum += sum_matrix(
1017  edge_data.dataOnEntities[MBEDGE][0].getN(AINSWORTH_LEGENDRE_BASE));
1018  diff_sum += sum_matrix(edge_data.dataOnEntities[MBEDGE][0].getDiffN(
1019  AINSWORTH_LEGENDRE_BASE));
1020  std::cout << "sum " << sum << std::endl;
1021  std::cout << "diff sum " << diff_sum << std::endl;
1022  if (std::abs(0.506122 - sum) > eps)
1023  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
1024 
1025  if (std::abs(2.85714 - diff_sum) > eps)
1026  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
1027  }
1028 
1029  if (choice_value == H1EDGE_BERNSTEIN_BEZIER) {
1030  CHKERR EdgePolynomialBase().getValue(
1031  pts_edge, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
1032  edge_data, "TET_FIELD", H1, CONTINUOUS,
1033  AINSWORTH_BERNSTEIN_BEZIER_BASE, NOBASE)));
1034  if (edge_data.dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE).get() ==
1035  edge_data.dataOnEntities[MBVERTEX][0]
1036  .getNSharedPtr(AINSWORTH_BERNSTEIN_BEZIER_BASE)
1037  .get())
1038  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1039  "Should be diffrent pointers");
1040 
1041  double sum = 0, diff_sum = 0;
1042  std::cout << "Vertex\n";
1043  std::cout << edge_data.dataOnEntities[MBVERTEX][0].getN("TET_FIELD")
1044  << std::endl;
1045  std::cout << edge_data.dataOnEntities[MBVERTEX][0].getDiffN("TET_FIELD")
1046  << std::endl;
1047  std::cout << "Edge\n";
1048  std::cout << edge_data.dataOnEntities[MBEDGE][0].getN("TET_FIELD")
1049  << std::endl;
1050  std::cout << edge_data.dataOnEntities[MBEDGE][0].getDiffN("TET_FIELD")
1051  << std::endl;
1052  sum +=
1053  sum_matrix(edge_data.dataOnEntities[MBVERTEX][0].getN("TET_FIELD"));
1054  diff_sum += sum_matrix(
1055  edge_data.dataOnEntities[MBVERTEX][0].getDiffN("TET_FIELD"));
1056  sum += sum_matrix(edge_data.dataOnEntities[MBEDGE][0].getN("TET_FIELD"));
1057  diff_sum +=
1058  sum_matrix(edge_data.dataOnEntities[MBEDGE][0].getDiffN("TET_FIELD"));
1059  std::cout << "sum " << sum << std::endl;
1060  std::cout << "diff sum " << diff_sum << std::endl;
1061  if (std::abs(4 - sum) > eps)
1062  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
1063 
1064  if (std::abs(diff_sum) > eps)
1065  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
1066  }
1067 
1068  if (choice_value == HCURLEDGE_AINSWORTH) {
1069  CHKERR EdgePolynomialBase().getValue(
1070  pts_edge,
1071  boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
1072  edge_data, HCURL, CONTINUOUS, AINSWORTH_LEGENDRE_BASE, NOBASE)));
1073  if (edge_data.dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE).get() !=
1074  edge_data.dataOnEntities[MBVERTEX][0]
1075  .getNSharedPtr(AINSWORTH_LEGENDRE_BASE)
1076  .get()) {
1077  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1078  "Different pointers");
1079  }
1080  std::cout << edge_data.dataOnEntities[MBEDGE][0].getN(
1081  AINSWORTH_LEGENDRE_BASE)
1082  << std::endl;
1083  double sum = 0;
1084  sum += sum_matrix(
1085  edge_data.dataOnEntities[MBEDGE][0].getN(AINSWORTH_LEGENDRE_BASE));
1086  std::cout.precision(std::numeric_limits<double>::max_digits10 - 1);
1087  std::cout << "sum " << sum << std::endl;
1088  if (fabs(-4.174149659863944 - sum) > eps) {
1089  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
1090  }
1091  }
1092 
1093  if (choice_value == HCURLEDGE_DEMKOWICZ) {
1094  CHKERR EdgePolynomialBase().getValue(
1095  pts_edge, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
1096  edge_data, HCURL, CONTINUOUS, DEMKOWICZ_JACOBI_BASE, NOBASE)));
1097  if (edge_data.dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE).get() !=
1098  edge_data.dataOnEntities[MBVERTEX][0]
1099  .getNSharedPtr(DEMKOWICZ_JACOBI_BASE)
1100  .get()) {
1101  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1102  "Different pointers");
1103  }
1104  std::cout << edge_data.dataOnEntities[MBEDGE][0].getN(
1105  DEMKOWICZ_JACOBI_BASE)
1106  << std::endl;
1107  double sum = sum_matrix(
1108  edge_data.dataOnEntities[MBEDGE][0].getN(DEMKOWICZ_JACOBI_BASE));
1109  std::cout.precision(std::numeric_limits<double>::max_digits10 - 1);
1110  std::cout << "sum " << sum << std::endl;
1111  if (std::fabs(4.571428571428571 - sum) > eps) {
1112  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
1113  }
1114  }
1115 
1116  if (choice_value == L2EDGE) {
1117  CHKERR EdgePolynomialBase().getValue(
1118  pts_edge,
1119  boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
1120  edge_data, L2, CONTINUOUS, DEMKOWICZ_JACOBI_BASE, NOBASE)));
1121  CHKERR EdgePolynomialBase().getValue(
1122  pts_edge,
1123  boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
1124  edge_data, L2, CONTINUOUS, AINSWORTH_LEGENDRE_BASE, NOBASE)));
1125  CHKERR EdgePolynomialBase().getValue(
1126  pts_edge,
1127  boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
1128  edge_data, L2, CONTINUOUS, AINSWORTH_LOBATTO_BASE, NOBASE)));
1129 
1130  MatrixDouble diff_n =
1131  edge_data.dataOnEntities[MBEDGE][0].getN(DEMKOWICZ_JACOBI_BASE) -
1132  edge_data.dataOnEntities[MBEDGE][0].getN(AINSWORTH_LEGENDRE_BASE);
1133  double sum = sum_matrix(diff_n);
1134  std::cout << "sum " << sum << std::endl;
1135  if (std::fabs(sum) > eps) {
1136  SETERRQ(
1137  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1138  "Difference between ainsworth and demkowicz base should be zero");
1139  }
1140 
1141  auto order = edge_data.dataOnEntities[MBEDGE][0].getOrder();
1142  auto &base =
1143  edge_data.dataOnEntities[MBEDGE][0].getN(AINSWORTH_LEGENDRE_BASE);
1144 
1145  if (base.size1() != pts_edge.size2())
1146  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1147  "wrong number of integration points");
1148  if (base.size2() != NBEDGE_L2(order))
1149  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1150  "wrong number of base functions");
1151 
1152  double ort_sum = 0;
1153  for (auto gg = 0; gg != pts_edge.size2(); ++gg) {
1154  double w = pts_edge(1, gg);
1155  for (int i = 0; i != NBEDGE_L2(order); ++i) {
1156  for (int j = 0; j != NBEDGE_L2(order); ++j) {
1157  if (i != j) {
1158  ort_sum += base(gg, i) * base(gg, j) * w;
1159  }
1160  }
1161  }
1162  }
1163  std::cout << "ort sum " << ort_sum << std::endl;
1164  if (std::fabs(ort_sum) > eps) {
1165  SETERRQ(
1166  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1167  "check orthogonality of base");
1168  }
1169 
1170  auto &diff_base =
1171  edge_data.dataOnEntities[MBEDGE][0].getDiffN(AINSWORTH_LOBATTO_BASE);
1172  double ort_diff_sum = 0;
1173  for (auto gg = 0; gg != pts_edge.size2(); ++gg) {
1174  double w = pts_edge(1, gg);
1175  for (int i = 2; i != NBEDGE_L2(order); ++i) {
1176  for (int j = 2; j != NBEDGE_L2(order); ++j) {
1177  if (i != j) {
1178  ort_diff_sum += diff_base(gg, i) * diff_base(gg, j) * w;
1179  }
1180  }
1181  }
1182  }
1183  std::cout << "ort diff sum " << ort_diff_sum << std::endl;
1184  if (std::fabs(ort_diff_sum) > eps) {
1185  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1186  "check orthogonality of diff lobatto base");
1187  }
1188 
1189  }
1190 
1191  EntitiesFieldData quad_data(MBQUAD);
1192  for (int type = MBVERTEX; type != MBMAXTYPE; type++) {
1193  quad_data.spacesOnEntities[type].set(H1);
1194  }
1195  quad_data.dataOnEntities[MBVERTEX].resize(1);
1196  quad_data.dataOnEntities[MBVERTEX][0].getOrder() = 1;
1197  quad_data.dataOnEntities[MBEDGE].resize(4);
1198  for (int ee = 0; ee < 4; ee++) {
1199  quad_data.dataOnEntities[MBEDGE][ee].getOrder() = 4;
1200  quad_data.dataOnEntities[MBEDGE][ee].getSense() = 1;
1201  }
1202  quad_data.dataOnEntities[MBQUAD].resize(1);
1203  quad_data.dataOnEntities[MBQUAD][0].getOrder() = 6;
1204 
1205  MatrixDouble pts_quad;
1206  int rule_ksi = 6;
1207  int rule_eta = 4;
1208  CHKERR Tools::outerProductOfEdgeIntegrationPtsForQuad(pts_quad, rule_ksi,
1209  rule_eta);
1210  nb_gauss_pts = pts_quad.size2();
1211  quad_data.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4,
1212  false);
1213  quad_data.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).resize(nb_gauss_pts,
1214  8, false);
1215  for (int i = 0; i < nb_gauss_pts; ++i) {
1216  quad_data.dataOnEntities[MBVERTEX][0].getN(NOBASE)(i, 0) =
1217  N_MBQUAD0(pts_quad(0, i), pts_quad(1, i));
1218  quad_data.dataOnEntities[MBVERTEX][0].getN(NOBASE)(i, 1) =
1219  N_MBQUAD1(pts_quad(0, i), pts_quad(1, i));
1220  quad_data.dataOnEntities[MBVERTEX][0].getN(NOBASE)(i, 2) =
1221  N_MBQUAD2(pts_quad(0, i), pts_quad(1, i));
1222  quad_data.dataOnEntities[MBVERTEX][0].getN(NOBASE)(i, 3) =
1223  N_MBQUAD3(pts_quad(0, i), pts_quad(1, i));
1224 
1225  quad_data.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE)(i, 0) =
1226  diffN_MBQUAD0x(pts_quad(1, i));
1227  quad_data.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE)(i, 1) =
1228  diffN_MBQUAD0y(pts_quad(0, i));
1229  quad_data.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE)(i, 2) =
1230  diffN_MBQUAD1x(pts_quad(1, i));
1231  quad_data.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE)(i, 3) =
1232  diffN_MBQUAD1y(pts_quad(0, i));
1233  quad_data.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE)(i, 4) =
1234  diffN_MBQUAD2x(pts_quad(1, i));
1235  quad_data.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE)(i, 5) =
1236  diffN_MBQUAD2y(pts_quad(0, i));
1237  quad_data.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE)(i, 6) =
1238  diffN_MBQUAD3x(pts_quad(1, i));
1239  quad_data.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE)(i, 7) =
1240  diffN_MBQUAD3y(pts_quad(0, i));
1241  }
1242 
1243  if (choice_value == H1QUAD) {
1244  CHKERR QuadPolynomialBase().getValue(
1245  pts_quad, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
1246  quad_data, H1, CONTINUOUS, AINSWORTH_LEGENDRE_BASE, NOBASE)));
1247  if (quad_data.dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE).get() !=
1248  quad_data.dataOnEntities[MBVERTEX][0]
1249  .getNSharedPtr(AINSWORTH_LEGENDRE_BASE)
1250  .get()) {
1251  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1252  "Different pointers");
1253  }
1254  double sum = 0, diff_sum = 0;
1255  for (int ee = 0; ee < 4; ee++) {
1256  sum += sum_matrix(
1257  quad_data.dataOnEntities[MBEDGE][ee].getN(AINSWORTH_LEGENDRE_BASE));
1258  diff_sum += sum_matrix(quad_data.dataOnEntities[MBEDGE][ee].getDiffN(
1259  AINSWORTH_LEGENDRE_BASE));
1260  }
1261  sum += sum_matrix(
1262  quad_data.dataOnEntities[MBQUAD][0].getN(AINSWORTH_LEGENDRE_BASE));
1263  diff_sum += sum_matrix(quad_data.dataOnEntities[MBQUAD][0].getDiffN(
1264  AINSWORTH_LEGENDRE_BASE));
1265 
1266  std::cout << sum << " " << diff_sum << endl;
1267 
1268  if (std::abs(3.58249 - sum) > eps)
1269  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
1270 
1271 
1272  if (std::abs(-0.134694 - diff_sum) > eps)
1273  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
1274 
1275  }
1276  }
1277  CATCH_ERRORS;
1278 
1279  MoFEM::Core::Finalize();
1280 }
MoFEM::LegendrePolynomial
Calculating Legendre base functions.
Definition: LegendrePolynomial.hpp:45
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
diffN_MBQUAD1x
#define diffN_MBQUAD1x(y)
Definition: fem_tools.h:63
MoFEM::LegendrePolynomialCtx
Class used to give arguments to Legendre base functions.
Definition: LegendrePolynomial.hpp:17
MoFEM::EntPolynomialBaseCtx
Class used to pass element data to calculate base functions on tet,triangle,edge.
Definition: EntPolynomialBaseCtx.hpp:22
NOBASE
@ NOBASE
Definition: definitions.h:59
diffN_MBQUAD1y
#define diffN_MBQUAD1y(x)
Definition: fem_tools.h:64
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
QUAD_2D_TABLE
static QUAD *const QUAD_2D_TABLE[]
Definition: quad.h:175
MoFEM::JacobiPolynomial
Calculating Legendre base functions.
Definition: JacobiPolynomial.hpp:51
MoFEM::LobattoPolynomial
Calculating Lobatto base functions.
Definition: LobattoPolynomial.hpp:35
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
diffN_MBQUAD2x
#define diffN_MBQUAD2x(y)
Definition: fem_tools.h:65
MoFEM.hpp
MoFEM::KernelLobattoPolynomial::getValue
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Definition: LobattoPolynomial.cpp:59
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
N_MBQUAD2
#define N_MBQUAD2(x, y)
quad shape function
Definition: fem_tools.h:59
MoFEM::EntitiesFieldData::spacesOnEntities
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
Definition: EntitiesFieldData.hpp:49
MOFEM_IMPOSSIBLE_CASE
@ MOFEM_IMPOSSIBLE_CASE
Definition: definitions.h:35
MoFEM::EntitiesFieldData::facesNodes
MatrixInt facesNodes
nodes on finite element faces
Definition: EntitiesFieldData.hpp:45
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
ShapeDiffMBTRI
PetscErrorCode ShapeDiffMBTRI(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:194
QUAD_1D_TABLE
static QUAD *const QUAD_1D_TABLE[]
Definition: quad.h:164
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::LegendrePolynomial::getValue
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Definition: LegendrePolynomial.cpp:22
diffN_MBQUAD3x
#define diffN_MBQUAD3x(y)
Definition: fem_tools.h:67
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
help
static char help[]
Definition: testing_base_functions.cpp:8
MoFEM::LobattoPolynomial::getValue
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Definition: LobattoPolynomial.cpp:24
MoFEM::JacobiPolynomialCtx
Class used to give arguments to Legendre base functions.
Definition: JacobiPolynomial.hpp:17
MoFEM::TriPolynomialBase
Calculate base functions on triangle.
Definition: TriPolynomialBase.hpp:16
NBEDGE_L2
#define NBEDGE_L2(P)
Number of base functions on edge from L2 space.
Definition: h1_hdiv_hcurl_l2.h:48
double
MoFEM::KernelLobattoPolynomialCtx
Class used to give arguments to Kernel Lobatto base functions.
Definition: LobattoPolynomial.hpp:51
convert.type
type
Definition: convert.py:64
sum_matrix
static double sum_matrix(MatrixDouble &m)
Definition: testing_base_functions.cpp:10
MoFEM::TriPolynomialBase::getValue
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Definition: TriPolynomialBase.cpp:1130
QUAD_::npoints
int npoints
Definition: quad.h:29
N_MBQUAD1
#define N_MBQUAD1(x, y)
quad shape function
Definition: fem_tools.h:58
diffN_MBQUAD2y
#define diffN_MBQUAD2y(x)
Definition: fem_tools.h:66
MoFEM::QuadPolynomialBase::getValue
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Definition: QuadPolynomialBase.cpp:430
AINSWORTH_LOBATTO_BASE
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
MoFEM::LobattoPolynomialCtx
Class used to give arguments to Lobatto base functions.
Definition: LobattoPolynomial.hpp:17
MoFEM::IntegratedJacobiPolynomialCtx
Definition: JacobiPolynomial.hpp:63
diffN_MBQUAD0x
#define diffN_MBQUAD0x(y)
Definition: fem_tools.h:61
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
QUAD_3D_TABLE
static QUAD *const QUAD_3D_TABLE[]
Definition: quad.h:187
MoFEM::QuadPolynomialBase
Calculate base functions on triangle.
Definition: QuadPolynomialBase.hpp:21
diffN_MBQUAD0y
#define diffN_MBQUAD0y(x)
Definition: fem_tools.h:62
MoFEM::EdgePolynomialBase::getValue
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Definition: EdgePolynomialBase.cpp:15
diffN_MBQUAD3y
#define diffN_MBQUAD3y(x)
Definition: fem_tools.h:68
AINSWORTH_BERNSTEIN_BEZIER_BASE
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
convert.n
n
Definition: convert.py:82
MoFEM::Tools::shapeFunMBTET
static MoFEMErrorCode shapeFunMBTET(double *shape, const double *ksi, const double *eta, const double *zeta, const double nb)
Calculate shape functions on tetrahedron.
Definition: Tools.hpp:738
MoFEM::EdgePolynomialBase
Calculate base functions on tetrahedral.
Definition: EdgePolynomialBase.hpp:19
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
MoFEM::JacobiPolynomial::getValue
MoFEMErrorCode getValue(MatrixDouble &pts_x, MatrixDouble &pts_t, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Definition: JacobiPolynomial.cpp:50
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
main
int main(int argc, char *argv[])
Definition: testing_base_functions.cpp:20
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
N_MBQUAD0
#define N_MBQUAD0(x, y)
quad shape function
Definition: fem_tools.h:57
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::PetscOptionsGetEList
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
Definition: DeprecatedPetsc.hpp:203
MoFEM::TetPolynomialBase::getValue
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Definition: TetPolynomialBase.cpp:2021
MoFEM::EntitiesFieldData::dataOnEntities
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
Definition: EntitiesFieldData.hpp:57
sdf_wavy_2d.w
int w
Definition: sdf_wavy_2d.py:6
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
N_MBQUAD3
#define N_MBQUAD3(x, y)
quad shape function
Definition: fem_tools.h:60
MoFEM::Tools::outerProductOfEdgeIntegrationPtsForQuad
static MoFEMErrorCode outerProductOfEdgeIntegrationPtsForQuad(MatrixDouble &pts, const int edge0, const int edge1)
Definition: Tools.cpp:610
MoFEM::EntitiesFieldData
data structure for finite element entity
Definition: EntitiesFieldData.hpp:40
MoFEM::KernelLobattoPolynomial
Calculating Lobatto base functions.
Definition: LobattoPolynomial.hpp:69
CONTINUOUS
@ CONTINUOUS
Regular field.
Definition: definitions.h:100
FractureMechanics::LASTOP
@ LASTOP
Definition: CrackFrontElement.hpp:23
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
MoFEM::IntegratedJacobiPolynomial::getValue
MoFEMErrorCode getValue(MatrixDouble &pts_x, MatrixDouble &pts_t, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Definition: JacobiPolynomial.cpp:70
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
MoFEM::IntegratedJacobiPolynomial
Definition: JacobiPolynomial.hpp:79
MoFEM::TetPolynomialBase
Calculate base functions on tetrahedral.
Definition: TetPolynomialBase.hpp:17
ShapeMBTET
PetscErrorCode ShapeMBTET(double *N, const double *G_X, const double *G_Y, const double *G_Z, int DIM)
calculate shape functions
Definition: fem_tools.c:306
quad.h
Generated by Doxygen 1.8.17 and hosted at University of Glasgow