20 {
21
23
24 try {
25
26 moab::Core mb_instance;
27 moab::Interface &moab = mb_instance;
29
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 H1FLATPRIS,
56 H1FATPRISM,
58 };
59
60 const char *list[] = {"legendre",
61 "lobatto",
62 "jacobi",
63 "integrated_jacobi",
64 "h1tet_ainsworth",
65 "h1tet_bernstein_bezier",
66 "hdivtet_ainsworth",
67 "hdivtet_demkowicz",
68 "hcurltet_ainsworth",
69 "hcurltet_demkowicz",
70 "l2tet",
71 "h1tri_ainsworth",
72 "h1tri_bernstein_bezier",
73 "h1quad",
74 "hdivtri_ainsworth",
75 "hdivtri_demkowicz",
76 "hcurltri_ainsworth",
77 "hcurltri_demkowicz",
78 "l2tri",
79 "h1edge_ainsworth",
80 "h1edge_bernstein_bezier",
81 "hcurledge_ainsworth",
82 "hcurledge_demkowicz",
83 "h1flatprism",
84 "h1fatprism"};
85
86 PetscBool flg;
87 PetscInt choice_value = LEGENDREPOLYNOMIAL;
89 &choice_value, &flg);
91 if (flg != PETSC_TRUE) {
93 }
94
96 pts_1d(0, 0) = -0.5;
97 pts_1d(0, 1) = 0.;
98 pts_1d(0, 2) = +0.5;
99
100 boost::shared_ptr<MatrixDouble> base_ptr(
new MatrixDouble);
101 boost::shared_ptr<MatrixDouble> diff_base_ptr(
new MatrixDouble);
102
103 const double eps = 1e-3;
104
105 if (choice_value == LEGENDREPOLYNOMIAL) {
106 double diff_s = 1;
109 4, &diff_s, 1, base_ptr, diff_base_ptr)));
111
112 std::cout << "LegendrePolynomial\n";
113 std::cout << pts_1d << std::endl;
114 std::cout << *base_ptr << std::endl;
115 std::cout << *diff_base_ptr << std::endl;
118 std::cout << sum << std::endl;
119 std::cout << diff_sum << std::endl;
120 if (fabs(2.04688 - sum) >
eps) {
122 }
123 if (fabs(2.25 - diff_sum) >
eps) {
125 }
126 }
127
128 pts_1d.resize(1, 11, false);
129 for (int ii = 0; ii != 11; ii++) {
130 pts_1d(0, ii) = 2 * ((
double)ii / 10) - 1;
131 }
132
133 boost::shared_ptr<MatrixDouble> kernel_base_ptr(
new MatrixDouble);
134 boost::shared_ptr<MatrixDouble> diff_kernel_base_ptr(
new MatrixDouble);
135
136 if (choice_value == LOBATTOPOLYNOMIAL) {
137 double diff_s = 1;
140 7 + 2, &diff_s, 1, base_ptr, diff_base_ptr)));
143 pts_1d,
145 7, &diff_s, 1, kernel_base_ptr, diff_kernel_base_ptr)));
147 for (int ii = 0; ii != 11; ii++) {
148 double s = pts_1d(0, ii);
149 std::cerr << "lobatto_plot " << s << " " << (*base_ptr)(ii, 1) << " "
150 << (*diff_base_ptr)(ii, 1) << " " << (*kernel_base_ptr)(ii, 1)
151 << " " << (*diff_kernel_base_ptr)(ii, 1) << " "
152 << (*kernel_base_ptr)(ii, 1) * (1 - s * s) << " "
153 << (*kernel_base_ptr)(ii, 1) * (-2 * s) +
154 (*diff_kernel_base_ptr)(ii, 1) * (1 - s * s)
155 << " " << std::endl;
156 }
157 std::cout << "LobattoPolynomial\n";
158 std::cout << pts_1d << std::endl;
159 std::cout << *base_ptr << std::endl;
160 std::cout << *diff_base_ptr << std::endl;
161 {
164 std::cout << sum << std::endl;
165 std::cout << diff_sum << std::endl;
166 if (fabs(9.39466 - sum) >
eps) {
168 }
169 if (fabs(14.0774 - diff_sum) >
eps) {
171 }
172 }
173 {
175 double diff_sum =
sum_matrix(*diff_kernel_base_ptr);
176 std::cout << sum << std::endl;
177 std::cout << diff_sum << std::endl;
178 if (fabs(-13.9906 * 4 - sum) >
eps) {
180 }
181 if (fabs(-101.678 * 4 - diff_sum) >
eps) {
183 }
184 }
185 }
186
187 if (choice_value == JACOBIPOLYNOMIAL) {
188
191 pts_1d.resize(1,
n,
false);
193 for (
int ii = 0; ii !=
n; ii++) {
194 pts_1d(0, ii) = (
double)ii / (
n - 1);
195 pts_1d_t(0, ii) = 1;
196 }
197
198 base_ptr->clear();
199 diff_base_ptr->clear();
200
201 double diff_x = 1;
202 double diff_t = 0;
204 pts_1d, pts_1d_t,
206 5, &diff_x, &diff_t, 1, 0, base_ptr, diff_base_ptr)));
208 for (
int ii = 0; ii !=
n; ii++) {
209 double s = pts_1d(0, ii);
210 std::cerr << "jacobi_plot " << s << " " << (*base_ptr)(ii, 4) << " "
211 << (*diff_base_ptr)(ii, 4) << std::endl;
212 }
213 for (
int ii = 0; ii !=
n - 1; ii++) {
214 double s = (pts_1d(0, ii) + pts_1d(0, ii + 1)) / 2;
215 std::cerr << "jacobi_diff_plot_fd " << s << " "
216 << ((*base_ptr)(ii + 1, 4) - (*base_ptr)(ii, 4)) /
218 << std::endl;
219 }
220 std::cout << "JacobiPolynomial\n";
221 std::cout << pts_1d << std::endl;
222 std::cout << *base_ptr << std::endl;
223 std::cout << *diff_base_ptr << std::endl;
224 {
227 std::cout << sum << std::endl;
228 std::cout << diff_sum << std::endl;
229 if (fabs(23.2164 - sum) >
eps) {
231 }
232 if (fabs(167.995 - diff_sum) >
eps) {
234 }
235 }
236 }
237
238 if (choice_value == INTEGRATEDJACOBIPOLYNOMIAL) {
239
242 pts_1d.resize(1,
n,
false);
244 for (
int ii = 0; ii !=
n; ii++) {
245 pts_1d(0, ii) = (
double)ii / 20.;
246 pts_1d_t(0, ii) = 1 - pts_1d(0, ii);
247 }
248
249 base_ptr->clear();
250 diff_base_ptr->clear();
251
252 double diff_x = 1;
253 double diff_t = 0;
255 pts_1d, pts_1d_t,
257 6, &diff_x, &diff_t, 1, 1, base_ptr, diff_base_ptr)));
259 for (
int ii = 0; ii !=
n; ii++) {
260 double s = pts_1d(0, ii);
261 std::cerr <<
"integrated_jacobi_plot " << s <<
" " << (
double)ii / 20.
262 << " " << (*base_ptr)(ii, 1) << " " << (*base_ptr)(ii, 2)
263 << " " << (*base_ptr)(ii, 3) << " " << (*base_ptr)(ii, 4)
264 << " " << (*base_ptr)(ii, 5) << endl;
265 ;
266
267 }
268 std::cout << "IntegratedJacobiPolynomial\n";
269 std::cout << pts_1d << std::endl;
270 std::cout << *base_ptr << std::endl;
271 std::cout << *diff_base_ptr << std::endl;
272 {
275 std::cout << sum << std::endl;
276 std::cout << diff_sum << std::endl;
277
278
279
280
281
282
283 }
284 }
285
287 for (
int type = MBVERTEX;
type != MBMAXTYPE;
type++) {
288 tet_data.spacesOnEntities[
type].set(
L2);
289 tet_data.spacesOnEntities[
type].set(
H1);
290 tet_data.spacesOnEntities[
type].set(
HDIV);
291 tet_data.spacesOnEntities[
type].set(
HCURL);
292 }
293 tet_data.dataOnEntities[MBVERTEX].resize(1);
294 tet_data.dataOnEntities[MBVERTEX][0].getOrder() = 1;
295 tet_data.dataOnEntities[MBEDGE].resize(6);
296 for (int ee = 0; ee < 6; ee++) {
297 tet_data.dataOnEntities[MBEDGE][ee].getOrder() = 3;
298 tet_data.dataOnEntities[MBEDGE][ee].getSense() = 1;
299 }
300 tet_data.dataOnEntities[MBTRI].resize(4);
301 for (int ff = 0; ff < 4; ff++) {
302 tet_data.dataOnEntities[MBTRI][ff].getOrder() = 4;
303 tet_data.dataOnEntities[MBTRI][ff].getSense() = 1;
304 }
305 tet_data.dataOnEntities[MBTET].resize(1);
306 tet_data.dataOnEntities[MBTET][0].getOrder() = 5;
307 tet_data.dataOnEntities[MBVERTEX][0].getBBNodeOrder().resize(4, false);
308 std::fill(tet_data.dataOnEntities[MBVERTEX][0].getBBNodeOrder().begin(),
309 tet_data.dataOnEntities[MBVERTEX][0].getBBNodeOrder().end(), 5);
310
312 int tet_rule = 2;
314 pts_tet.resize(3, nb_gauss_pts, false);
315 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[tet_rule]->points[1], 4,
316 &pts_tet(0, 0), 1);
317 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[tet_rule]->points[2], 4,
318 &pts_tet(1, 0), 1);
319 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[tet_rule]->points[3], 4,
320 &pts_tet(2, 0), 1);
321 tet_data.dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 4,
322 false);
323 {
324 double *shape_ptr =
325 &*tet_data.dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
326 cblas_dcopy(4 * nb_gauss_pts,
QUAD_3D_TABLE[tet_rule]->points, 1,
327 shape_ptr, 1);
328 }
329 tet_data.facesNodes.resize(4, 3);
330 const int cannonical_tet_face[4][3] = {
331 {0, 1, 3}, {1, 2, 3}, {0, 3, 2}, {0, 2, 1}};
332 for (int ff = 0; ff < 4; ff++) {
333 for (int nn = 0; nn < 3; nn++) {
334 tet_data.facesNodes(ff, nn) = cannonical_tet_face[ff][nn];
335 }
336 }
337
338 if (choice_value == H1TET_AINSWORTH) {
339
340 tet_data.dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 4,
341 false);
343 &*tet_data.dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin(),
344 &pts_tet(0, 0), &pts_tet(1, 0), &pts_tet(2, 0), nb_gauss_pts);
346
351 if (tet_data.dataOnEntities[MBVERTEX][0].getNSharedPtr(
NOBASE).get() !=
352 tet_data.dataOnEntities[MBVERTEX][0]
354 .get()) {
356 "Different pointers");
357 }
358
359 double sum = 0, diff_sum = 0;
360 std::cout << "Edges\n";
361 for (int ee = 0; ee < 6; ee++) {
362 std::cout << tet_data.dataOnEntities[MBEDGE][ee].getN(
364 << std::endl;
365 std::cout << tet_data.dataOnEntities[MBEDGE][ee].getDiffN(
367 << std::endl;
370 diff_sum +=
sum_matrix(tet_data.dataOnEntities[MBEDGE][ee].getDiffN(
372 }
373 std::cout << "Faces\n";
374 for (int ff = 0; ff < 4; ff++) {
375 std::cout << tet_data.dataOnEntities[MBTRI][ff].getN(
377 << std::endl;
378 std::cout << tet_data.dataOnEntities[MBTRI][ff].getDiffN(
380 << std::endl;
383 diff_sum +=
sum_matrix(tet_data.dataOnEntities[MBTRI][ff].getDiffN(
385 }
386 std::cout << "Tets\n";
387 std::cout << tet_data.dataOnEntities[MBTET][0].getN(
389 << std::endl;
390 std::cout << tet_data.dataOnEntities[MBTET][0].getDiffN(
392 << std::endl;
397 std::cout << "sum " << sum << std::endl;
398 std::cout << "diff_sum " << diff_sum << std::endl;
399 if (fabs(1.3509 - sum) >
eps) {
401 }
402 if (fabs(0.233313 - diff_sum) >
eps) {
404 }
405 }
406
407 if (choice_value == H1TET_BERNSTEIN_BEZIER) {
408
409 tet_data.dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 4,
410 false);
412 &*tet_data.dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin(),
413 &pts_tet(0, 0), &pts_tet(1, 0), &pts_tet(2, 0), nb_gauss_pts);
416 tet_data,
"TEST_FIELD",
H1,
418 if (tet_data.dataOnEntities[MBVERTEX][0].getNSharedPtr(
NOBASE).get() ==
419 tet_data.dataOnEntities[MBVERTEX][0]
421 .get())
423
424 double sum = 0, diff_sum = 0;
425 std::cout << "Vertices\n";
426 std::cout << tet_data.dataOnEntities[MBVERTEX][0].getN("TEST_FIELD")
427 << std::endl;
428 std::cout << tet_data.dataOnEntities[MBVERTEX][0].getDiffN("TEST_FIELD")
429 << std::endl;
430 sum +=
431 sum_matrix(tet_data.dataOnEntities[MBVERTEX][0].getN(
"TEST_FIELD"));
433 tet_data.dataOnEntities[MBVERTEX][0].getDiffN("TEST_FIELD"));
434 std::cout << "Edges\n";
435 for (int ee = 0; ee < 6; ee++) {
436 std::cout << tet_data.dataOnEntities[MBEDGE][ee].getN("TEST_FIELD")
437 << std::endl;
438 std::cout << tet_data.dataOnEntities[MBEDGE][ee].getDiffN("TEST_FIELD")
439 << std::endl;
440 sum +=
441 sum_matrix(tet_data.dataOnEntities[MBEDGE][ee].getN(
"TEST_FIELD"));
443 tet_data.dataOnEntities[MBEDGE][ee].getDiffN("TEST_FIELD"));
444 }
445 std::cout << "Faces\n";
446 for (int ff = 0; ff < 4; ff++) {
447 std::cout << tet_data.dataOnEntities[MBTRI][ff].getN("TEST_FIELD")
448 << std::endl;
449 std::cout << tet_data.dataOnEntities[MBTRI][ff].getDiffN("TEST_FIELD")
450 << std::endl;
451 sum +=
452 sum_matrix(tet_data.dataOnEntities[MBTRI][ff].getN(
"TEST_FIELD"));
454 tet_data.dataOnEntities[MBTRI][ff].getDiffN("TEST_FIELD"));
455 }
456 std::cout << "Tets\n";
457 std::cout << tet_data.dataOnEntities[MBTET][0].getN("TEST_FIELD")
458 << std::endl;
459 std::cout << tet_data.dataOnEntities[MBTET][0].getDiffN("TEST_FIELD")
460 << std::endl;
461 sum +=
sum_matrix(tet_data.dataOnEntities[MBTET][0].getN(
"TEST_FIELD"));
462 diff_sum +=
463 sum_matrix(tet_data.dataOnEntities[MBTET][0].getDiffN(
"TEST_FIELD"));
464 std::cout << "sum " << sum << std::endl;
465 std::cout << "diff_sum " << diff_sum << std::endl;
466 if (fabs(4.38395 - sum) >
eps) {
468 }
469 if (fabs(diff_sum) >
eps) {
471 }
472 }
473
474 if (choice_value == HDIVTET_AINSWORTH) {
479 double sum = 0, diff_sum = 0;
480 std::cout << "Faces\n";
481 for (int ff = 0; ff < 4; ff++) {
482 std::cout << tet_data.dataOnEntities[MBTRI][ff].getN(
484 << std::endl;
485 std::cout << tet_data.dataOnEntities[MBTRI][ff].getDiffN(
487 << std::endl;
490 diff_sum +=
sum_matrix(tet_data.dataOnEntities[MBTRI][ff].getDiffN(
492 }
493 std::cout << "Tets\n";
494 std::cout << tet_data.dataOnEntities[MBTET][0].getN(
496 << std::endl;
497 std::cout << tet_data.dataOnEntities[MBTET][0].getDiffN(
499 << std::endl;
504 std::cout << "sum " << 1e8 * sum << std::endl;
505 std::cout << "diff_sum " << 1e8 * diff_sum << std::endl;
506 if (fabs(0.188636 - sum) >
eps) {
508 }
509 if (fabs(32.9562 - diff_sum) >
eps) {
511 }
512 }
513
514 if (choice_value == HDIVTET_DEMKOWICZ) {
519 double sum = 0, diff_sum = 0;
520 std::cout << "Faces\n";
521 for (int ff = 0; ff < 4; ff++) {
522 std::cout << tet_data.dataOnEntities[MBTRI][ff].getN(
524 << std::endl;
525 std::cout << tet_data.dataOnEntities[MBTRI][ff].getDiffN(
527 << std::endl;
532 }
533 std::cout << "Tets\n";
535 << std::endl;
536 std::cout << tet_data.dataOnEntities[MBTET][0].getDiffN(
538 << std::endl;
543 std::cout << "sum " << sum << std::endl;
544 std::cout << "diff_sum " << diff_sum << std::endl;
545 const double expected_result = -2.70651;
546 const double expected_diff_result = 289.421;
547 if (fabs((expected_result - sum) / expected_result) >
eps) {
549 }
550 if (fabs((expected_diff_result - diff_sum) / expected_diff_result) >
553 }
554 }
555
556 if (choice_value == HCURLTET_AINSWORTH) {
561 double sum = 0, diff_sum = 0;
562 std::cout << "Edges\n";
563 for (int ee = 0; ee < 6; ee++) {
564 std::cout << tet_data.dataOnEntities[MBEDGE][ee].getN(
566 << std::endl;
567 std::cout << tet_data.dataOnEntities[MBEDGE][ee].getDiffN(
569 << std::endl;
572 diff_sum +=
sum_matrix(tet_data.dataOnEntities[MBEDGE][ee].getDiffN(
574 }
575 std::cout << "Faces\n";
576 for (int ff = 0; ff < 4; ff++) {
577 std::cout << tet_data.dataOnEntities[MBTRI][ff].getN(
579 << std::endl;
580 std::cout << tet_data.dataOnEntities[MBTRI][ff].getDiffN(
582 << std::endl;
585 diff_sum +=
sum_matrix(tet_data.dataOnEntities[MBTRI][ff].getDiffN(
587 }
588 std::cout << "Tets\n";
589 std::cout << tet_data.dataOnEntities[MBTET][0].getN(
591 << std::endl;
592 std::cout << tet_data.dataOnEntities[MBTET][0].getDiffN(
594 << std::endl;
599 std::cout << "sum " << sum << std::endl;
600 std::cout << "diff_sum " << diff_sum << std::endl;
601 if (fabs(-1.7798 - sum) >
eps) {
603 }
604 if (fabs(-67.1793 - diff_sum) >
eps) {
606 }
607 }
608
609 if (choice_value == HCURLTET_DEMKOWICZ) {
613 double sum = 0, diff_sum = 0;
614 std::cout << "Edges\n";
615 for (int ee = 0; ee < 6; ee++) {
616 std::cout << tet_data.dataOnEntities[MBEDGE][ee].getN(
618 << std::endl;
619 std::cout << tet_data.dataOnEntities[MBEDGE][ee].getDiffN(
621 << std::endl;
624 diff_sum +=
sum_matrix(tet_data.dataOnEntities[MBEDGE][ee].getDiffN(
626 }
627 std::cout << "Faces\n";
628 for (int ff = 0; ff < 4; ff++) {
629 std::cout << tet_data.dataOnEntities[MBTRI][ff].getN(
631 << std::endl;
632 std::cout << tet_data.dataOnEntities[MBTRI][ff].getDiffN(
634 << std::endl;
639 }
640 std::cout << "Tets\n";
642 << std::endl;
643 std::cout << tet_data.dataOnEntities[MBTET][0].getDiffN(
645 << std::endl;
650 std::cout << "sum " << sum << std::endl;
651 std::cout << "diff_sum " << diff_sum << std::endl;
652 if (fabs(7.35513 - sum) >
eps) {
654 }
655 if (fabs(62.4549 - diff_sum) >
eps) {
657 }
658 }
659
660 if (choice_value == L2TET) {
665 double sum = 0, diff_sum = 0;
666 std::cout << "Tets\n";
667 std::cout << tet_data.dataOnEntities[MBTET][0].getN(
669 << std::endl;
670 std::cout << tet_data.dataOnEntities[MBTET][0].getDiffN(
672 << std::endl;
677 std::cout << "sum " << sum << std::endl;
678 std::cout << "diff_sum " << diff_sum << std::endl;
679 if (fabs(3.60352 - sum) >
eps) {
681 }
682 if (fabs(-36.9994 - diff_sum) >
eps) {
684 }
685 }
686
688 for (
int type = MBVERTEX;
type != MBMAXTYPE;
type++) {
689 tri_data.spacesOnEntities[
type].set(
L2);
690 tri_data.spacesOnEntities[
type].set(
H1);
691 tri_data.spacesOnEntities[
type].set(
HDIV);
692 tri_data.spacesOnEntities[
type].set(
HCURL);
693 }
694 tri_data.dataOnEntities[MBVERTEX].resize(1);
695 tri_data.dataOnEntities[MBVERTEX][0].getOrder() = 1;
696 tri_data.dataOnEntities[MBEDGE].resize(3);
697 for (int ee = 0; ee < 3; ee++) {
698 tri_data.dataOnEntities[MBEDGE][ee].getOrder() = 3;
699 tri_data.dataOnEntities[MBEDGE][ee].getSense() = 1;
700 }
701 tri_data.dataOnEntities[MBTRI].resize(1);
702 tri_data.dataOnEntities[MBTRI][0].getOrder() = 4;
703 tri_data.dataOnEntities[MBVERTEX][0].getBBNodeOrder().resize(3, false);
704 std::fill(tri_data.dataOnEntities[MBVERTEX][0].getBBNodeOrder().begin(),
705 tri_data.dataOnEntities[MBVERTEX][0].getBBNodeOrder().end(), 4);
706
708 int tri_rule = 2;
710 pts_tri.resize(2, nb_gauss_pts, false);
711 cblas_dcopy(nb_gauss_pts, &
QUAD_2D_TABLE[tri_rule]->points[1], 3,
712 &pts_tri(0, 0), 1);
713 cblas_dcopy(nb_gauss_pts, &
QUAD_2D_TABLE[tri_rule]->points[2], 3,
714 &pts_tri(1, 0), 1);
715 tri_data.dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 3,
716 false);
717 {
718 double *shape_ptr =
719 &*tri_data.dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
720 cblas_dcopy(3 * nb_gauss_pts,
QUAD_2D_TABLE[tri_rule]->points, 1,
721 shape_ptr, 1);
722 }
723 tri_data.dataOnEntities[MBVERTEX][0].getDiffN(
NOBASE).resize(3, 2,
false);
725 &*tri_data.dataOnEntities[MBVERTEX][0].getDiffN(
NOBASE).data().begin());
727
728 if (choice_value == H1TRI_AINSWORTH) {
732 if (tri_data.dataOnEntities[MBVERTEX][0].getNSharedPtr(
NOBASE).get() !=
733 tri_data.dataOnEntities[MBVERTEX][0]
735 .get())
737 "Different pointers");
738
739 double sum = 0, diff_sum = 0;
740 std::cout << "Edges\n";
741 for (int ee = 0; ee < 3; ee++) {
742 std::cout << tri_data.dataOnEntities[MBEDGE][ee].getN(
744 << std::endl;
745 std::cout << tri_data.dataOnEntities[MBEDGE][ee].getDiffN(
747 << std::endl;
750 diff_sum +=
sum_matrix(tri_data.dataOnEntities[MBEDGE][ee].getDiffN(
752 }
753 std::cout << "Face\n";
754 std::cout << tri_data.dataOnEntities[MBTRI][0].getN(
756 << std::endl;
757 std::cout << tri_data.dataOnEntities[MBTRI][0].getDiffN(
759 << std::endl;
764 std::cout << "sum " << sum << std::endl;
765 std::cout << "diff_sum " << diff_sum << std::endl;
766 if (fabs(0.805556 - sum) >
eps)
768
769 if (fabs(0.0833333 - diff_sum) >
eps)
771 }
772
773 if (choice_value == H1TRI_BERNSTEIN_BEZIER) {
776 tri_data,
"TET_FIELD",
H1,
778 if (tri_data.dataOnEntities[MBVERTEX][0].getNSharedPtr(
NOBASE).get() ==
779 tri_data.dataOnEntities[MBVERTEX][0]
781 .get())
783 "Pointers should be diffrent");
784
785 double sum = 0, diff_sum = 0;
786 std::cout << "Vertex\n";
787 std::cout << tri_data.dataOnEntities[MBVERTEX][0].getN("TET_FIELD")
788 << std::endl;
789 std::cout << tri_data.dataOnEntities[MBVERTEX][0].getDiffN("TET_FIELD")
790 << std::endl;
791 sum +=
sum_matrix(tri_data.dataOnEntities[MBVERTEX][0].getN(
"TET_FIELD"));
793 tri_data.dataOnEntities[MBVERTEX][0].getDiffN("TET_FIELD"));
794 std::cout << "Edges\n";
795 for (int ee = 0; ee < 3; ee++) {
796 std::cout << tri_data.dataOnEntities[MBEDGE][ee].getN("TET_FIELD")
797 << std::endl;
798 std::cout << tri_data.dataOnEntities[MBEDGE][ee].getDiffN("TET_FIELD")
799 << std::endl;
800 sum +=
801 sum_matrix(tri_data.dataOnEntities[MBEDGE][ee].getN(
"TET_FIELD"));
803 tri_data.dataOnEntities[MBEDGE][ee].getDiffN("TET_FIELD"));
804 }
805 std::cout << "Face\n";
806 std::cout << tri_data.dataOnEntities[MBTRI][0].getN("TET_FIELD")
807 << std::endl;
808 std::cout << tri_data.dataOnEntities[MBTRI][0].getDiffN("TET_FIELD")
809 << std::endl;
810 sum +=
sum_matrix(tri_data.dataOnEntities[MBTRI][0].getN(
"TET_FIELD"));
811 diff_sum +=
812 sum_matrix(tri_data.dataOnEntities[MBTRI][0].getDiffN(
"TET_FIELD"));
813 std::cout << "sum " << sum << std::endl;
814 std::cout << "diff_sum " << diff_sum << std::endl;
815 if (std::abs(3.01389 - sum) >
eps)
817
818 if (std::abs(diff_sum) >
eps)
820 "wrong result %3.4e != $3.4e", 0, diff_sum);
821 }
822
823 if (choice_value == HDIVTRI_AINSWORTH) {
827 if (tri_data.dataOnEntities[MBVERTEX][0].getNSharedPtr(
NOBASE).get() !=
828 tri_data.dataOnEntities[MBVERTEX][0]
830 .get())
832 "Different pointers");
833
834 double sum = 0;
835 std::cout << "Face\n";
836 std::cout << tri_data.dataOnEntities[MBTRI][0].getN(
838 << std::endl;
841 std::cout << "sum " << sum << std::endl;
842 if (fabs(1.93056 - sum) >
eps)
844 }
845
846 if (choice_value == HDIVTRI_DEMKOWICZ) {
851 if (tri_data.dataOnEntities[MBVERTEX][0].getNSharedPtr(
NOBASE).get() !=
852 tri_data.dataOnEntities[MBVERTEX][0]
854 .get()) {
856 "Different pointers");
857 }
858 double sum = 0;
859 std::cout << "Face\n";
861 << std::endl;
864 std::cout << "sum " << sum << std::endl;
865 const double expected_result = 28.25;
866 if (fabs((expected_result - sum) / expected_result) >
eps) {
868 }
869 }
870
871 if (choice_value == HCURLTRI_AINSWORTH) {
876 if (tri_data.dataOnEntities[MBVERTEX][0].getNSharedPtr(
NOBASE).get() !=
877 tri_data.dataOnEntities[MBVERTEX][0]
879 .get()) {
881 "Different pointers");
882 }
883 double sum = 0;
884 std::cout << "Edges\n";
885 for (int ee = 0; ee < 3; ee++) {
886 std::cout << tri_data.dataOnEntities[MBEDGE][ee].getN(
888 << std::endl;
891 }
892 std::cout << "Face\n";
893 std::cout << tri_data.dataOnEntities[MBTRI][0].getN(
895 << std::endl;
898 std::cout << "sum " << sum << std::endl;
899 if (fabs(0.333333 - sum) >
eps) {
901 }
902 }
903
904 if (choice_value == HCURLTRI_DEMKOWICZ) {
908 if (tri_data.dataOnEntities[MBVERTEX][0].getNSharedPtr(
NOBASE).get() !=
909 tri_data.dataOnEntities[MBVERTEX][0]
911 .get()) {
913 "Different pointers");
914 }
915 double sum = 0;
916 std::cout << "Edges\n";
917 for (int ee = 0; ee < 3; ee++) {
918 std::cout << tri_data.dataOnEntities[MBEDGE][ee].getN(
920 << std::endl;
923 }
924 std::cout << "Face\n";
926 << std::endl;
929 std::cout << "sum " << sum << std::endl;
930 if (fabs(1.04591 - sum) >
eps) {
932 }
933 }
934
935 if (choice_value == L2TRI) {
940 if (tri_data.dataOnEntities[MBVERTEX][0].getNSharedPtr(
NOBASE).get() !=
941 tri_data.dataOnEntities[MBVERTEX][0]
943 .get()) {
945 "Different pointers");
946 }
947 double sum = 0;
948 std::cout << "Face\n";
949 std::cout << tri_data.dataOnEntities[MBTRI][0].getN(
951 << std::endl;
954 std::cout << "sum " << sum << std::endl;
955 if (fabs(0.671875 - sum) >
eps) {
957 }
958 }
959
961 for (
int type = MBVERTEX;
type != MBMAXTYPE;
type++) {
962 edge_data.spacesOnEntities[
type].set(
L2);
963 edge_data.spacesOnEntities[
type].set(
H1);
964 edge_data.spacesOnEntities[
type].set(
HDIV);
965 edge_data.spacesOnEntities[
type].set(
HCURL);
966 }
967 edge_data.dataOnEntities[MBVERTEX].resize(1);
968 edge_data.dataOnEntities[MBVERTEX][0].getOrder() = 1;
969 edge_data.dataOnEntities[MBEDGE].resize(1);
970 edge_data.dataOnEntities[MBEDGE][0].getOrder() = 4;
971 edge_data.dataOnEntities[MBEDGE][0].getSense() = 1;
972 edge_data.dataOnEntities[MBVERTEX][0].getBBNodeOrder().resize(2, false);
973 std::fill(edge_data.dataOnEntities[MBVERTEX][0].getBBNodeOrder().begin(),
974 edge_data.dataOnEntities[MBVERTEX][0].getBBNodeOrder().end(), 4);
975
977 int edge_rule = 6;
979 pts_edge.resize(1, nb_gauss_pts, false);
980 cblas_dcopy(nb_gauss_pts, &
QUAD_1D_TABLE[edge_rule]->points[1], 2,
981 &pts_edge(0, 0), 1);
982 edge_data.dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 2,
983 false);
984 {
985 double *shape_ptr =
986 &*edge_data.dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
987 cblas_dcopy(2 * nb_gauss_pts,
QUAD_1D_TABLE[edge_rule]->points, 1,
988 shape_ptr, 1);
989 }
990
991 if (choice_value == H1EDGE_AINSWORTH) {
995 if (edge_data.dataOnEntities[MBVERTEX][0].getNSharedPtr(
NOBASE).get() !=
996 edge_data.dataOnEntities[MBVERTEX][0]
998 .get()) {
1000 "Different pointers");
1001 }
1002 double sum = 0, diff_sum = 0;
1003 std::cout << "Edge\n";
1004 std::cout << edge_data.dataOnEntities[MBEDGE][0].getN(
1006 << std::endl;
1007 std::cout << edge_data.dataOnEntities[MBEDGE][0].getDiffN(
1009 << std::endl;
1012 diff_sum +=
sum_matrix(edge_data.dataOnEntities[MBEDGE][0].getDiffN(
1014 std::cout << "sum " << sum << std::endl;
1015 std::cout << "diff sum " << diff_sum << std::endl;
1016 if (std::abs(0.506122 - sum) >
eps)
1018
1019 if (std::abs(2.85714 - diff_sum) >
eps)
1021 }
1022
1023 if (choice_value == H1EDGE_BERNSTEIN_BEZIER) {
1026 edge_data,
"TET_FIELD",
H1,
1028 if (edge_data.dataOnEntities[MBVERTEX][0].getNSharedPtr(
NOBASE).get() ==
1029 edge_data.dataOnEntities[MBVERTEX][0]
1031 .get())
1033 "Should be diffrent pointers");
1034
1035 double sum = 0, diff_sum = 0;
1036 std::cout << "Vertex\n";
1037 std::cout << edge_data.dataOnEntities[MBVERTEX][0].getN("TET_FIELD")
1038 << std::endl;
1039 std::cout << edge_data.dataOnEntities[MBVERTEX][0].getDiffN("TET_FIELD")
1040 << std::endl;
1041 std::cout << "Edge\n";
1042 std::cout << edge_data.dataOnEntities[MBEDGE][0].getN("TET_FIELD")
1043 << std::endl;
1044 std::cout << edge_data.dataOnEntities[MBEDGE][0].getDiffN("TET_FIELD")
1045 << std::endl;
1046 sum +=
1047 sum_matrix(edge_data.dataOnEntities[MBVERTEX][0].getN(
"TET_FIELD"));
1049 edge_data.dataOnEntities[MBVERTEX][0].getDiffN("TET_FIELD"));
1050 sum +=
sum_matrix(edge_data.dataOnEntities[MBEDGE][0].getN(
"TET_FIELD"));
1051 diff_sum +=
1052 sum_matrix(edge_data.dataOnEntities[MBEDGE][0].getDiffN(
"TET_FIELD"));
1053 std::cout << "sum " << sum << std::endl;
1054 std::cout << "diff sum " << diff_sum << std::endl;
1055 if (std::abs(4 - sum) >
eps)
1057
1058 if (std::abs(diff_sum) >
eps)
1060 }
1061
1062 if (choice_value == HCURLEDGE_AINSWORTH) {
1067 if (edge_data.dataOnEntities[MBVERTEX][0].getNSharedPtr(
NOBASE).get() !=
1068 edge_data.dataOnEntities[MBVERTEX][0]
1070 .get()) {
1072 "Different pointers");
1073 }
1074 std::cout << edge_data.dataOnEntities[MBEDGE][0].getN(
1076 << std::endl;
1077 int sum = 0;
1080 std::cout << "sum " << sum << std::endl;
1081 if (fabs(-4 - sum) >
eps) {
1083 }
1084 }
1085
1086 if (choice_value == HCURLEDGE_DEMKOWICZ) {
1090 if (edge_data.dataOnEntities[MBVERTEX][0].getNSharedPtr(
NOBASE).get() !=
1091 edge_data.dataOnEntities[MBVERTEX][0]
1093 .get()) {
1095 "Different pointers");
1096 }
1097 std::cout << edge_data.dataOnEntities[MBEDGE][0].getN(
1099 << std::endl;
1100 int sum = 0;
1103 std::cout << "sum " << sum << std::endl;
1104 if (fabs(4 - sum) >
eps) {
1106 }
1107 }
1108
1110 for (
int type = MBVERTEX;
type != MBMAXTYPE;
type++) {
1111 quad_data.spacesOnEntities[
type].set(
H1);
1112 }
1113 quad_data.dataOnEntities[MBVERTEX].resize(1);
1114 quad_data.dataOnEntities[MBVERTEX][0].getOrder() = 1;
1115 quad_data.dataOnEntities[MBEDGE].resize(4);
1116 for (int ee = 0; ee < 4; ee++) {
1117 quad_data.dataOnEntities[MBEDGE][ee].getOrder() = 4;
1118 quad_data.dataOnEntities[MBEDGE][ee].getSense() = 1;
1119 }
1120 quad_data.dataOnEntities[MBQUAD].resize(1);
1121 quad_data.dataOnEntities[MBQUAD][0].getOrder() = 6;
1122
1124 int rule_ksi = 6;
1125 int rule_eta = 4;
1127 rule_eta);
1128 nb_gauss_pts = pts_quad.size2();
1129 quad_data.dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 4,
1130 false);
1131 quad_data.dataOnEntities[MBVERTEX][0].getDiffN(
NOBASE).resize(nb_gauss_pts,
1132 8, false);
1133 for (
int i = 0;
i < nb_gauss_pts; ++
i) {
1134 quad_data.dataOnEntities[MBVERTEX][0].getN(
NOBASE)(
i, 0) =
1136 quad_data.dataOnEntities[MBVERTEX][0].getN(
NOBASE)(
i, 1) =
1138 quad_data.dataOnEntities[MBVERTEX][0].getN(
NOBASE)(
i, 2) =
1140 quad_data.dataOnEntities[MBVERTEX][0].getN(
NOBASE)(
i, 3) =
1142
1143 quad_data.dataOnEntities[MBVERTEX][0].getDiffN(
NOBASE)(
i, 0) =
1145 quad_data.dataOnEntities[MBVERTEX][0].getDiffN(
NOBASE)(
i, 1) =
1147 quad_data.dataOnEntities[MBVERTEX][0].getDiffN(
NOBASE)(
i, 2) =
1149 quad_data.dataOnEntities[MBVERTEX][0].getDiffN(
NOBASE)(
i, 3) =
1151 quad_data.dataOnEntities[MBVERTEX][0].getDiffN(
NOBASE)(
i, 4) =
1153 quad_data.dataOnEntities[MBVERTEX][0].getDiffN(
NOBASE)(
i, 5) =
1155 quad_data.dataOnEntities[MBVERTEX][0].getDiffN(
NOBASE)(
i, 6) =
1157 quad_data.dataOnEntities[MBVERTEX][0].getDiffN(
NOBASE)(
i, 7) =
1159 }
1160
1161 if (choice_value == H1QUAD) {
1165 if (quad_data.dataOnEntities[MBVERTEX][0].getNSharedPtr(
NOBASE).get() !=
1166 quad_data.dataOnEntities[MBVERTEX][0]
1168 .get()) {
1170 "Different pointers");
1171 }
1172 double sum = 0, diff_sum = 0;
1173 for (int ee = 0; ee < 4; ee++) {
1176 diff_sum +=
sum_matrix(quad_data.dataOnEntities[MBEDGE][ee].getDiffN(
1178 }
1181 diff_sum +=
sum_matrix(quad_data.dataOnEntities[MBQUAD][0].getDiffN(
1183
1184 std::cout << sum << " " << diff_sum << endl;
1185
1186 if (std::abs(3.58249 - sum) >
eps)
1188
1189
1190 if (std::abs(-0.134694 - diff_sum) >
eps)
1192
1193 }
1194 }
1196
1198}
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
@ L2
field with C-1 continuity
@ HCURL
field with continuous tangents
@ HDIV
field with continuous normal traction
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
@ MOFEM_DATA_INCONSISTENCY
#define CHKERR
Inline error check.
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'i', SPACE_DIM > i
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
static QUAD *const QUAD_2D_TABLE[]
static QUAD *const QUAD_1D_TABLE[]
static QUAD *const QUAD_3D_TABLE[]
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Calculate base functions on tetrahedral.
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Class used to pass element data to calculate base functions on tet,triangle,edge.
data structure for finite element entity
MoFEMErrorCode getValue(MatrixDouble &pts_x, MatrixDouble &pts_t, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Class used to give arguments to Legendre base functions.
Calculating Legendre base functions.
MoFEMErrorCode getValue(MatrixDouble &pts_x, MatrixDouble &pts_t, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Class used to give arguments to Kernel Lobatto base functions.
Calculating Lobatto base functions.
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Class used to give arguments to Legendre base functions.
Calculating Legendre base functions.
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Class used to give arguments to Lobatto base functions.
Calculating Lobatto base functions.
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Calculate base functions on triangle.
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Calculate base functions on tetrahedral.
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Calculate base functions on triangle.
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
static double sum_matrix(MatrixDouble &m)