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