v0.16.0
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 tet_data.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4,
429 false);
430 CHKERR Tools::shapeFunMBTET(
431 &*tet_data.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin(),
432 &pts_tet(0, 0), &pts_tet(1, 0), &pts_tet(2, 0), nb_gauss_pts);
433 CHKERR TetPolynomialBase().getValue(
434 pts_tet, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
435 tet_data, "TEST_FIELD", H1, CONTINUOUS,
436 AINSWORTH_BERNSTEIN_BEZIER_BASE, NOBASE)));
437 if (tet_data.dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE).get() ==
438 tet_data.dataOnEntities[MBVERTEX][0]
439 .getNSharedPtr(AINSWORTH_BERNSTEIN_BEZIER_BASE)
440 .get())
441 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "The same pointers");
442
443 double sum = 0, diff_sum = 0;
444 std::cout << "Vertices\n";
445 std::cout << tet_data.dataOnEntities[MBVERTEX][0].getN("TEST_FIELD")
446 << std::endl;
447 std::cout << tet_data.dataOnEntities[MBVERTEX][0].getDiffN("TEST_FIELD")
448 << std::endl;
449 sum +=
450 sum_matrix(tet_data.dataOnEntities[MBVERTEX][0].getN("TEST_FIELD"));
451 diff_sum += sum_matrix(
452 tet_data.dataOnEntities[MBVERTEX][0].getDiffN("TEST_FIELD"));
453 std::cout << "Edges\n";
454 for (int ee = 0; ee < 6; ee++) {
455 std::cout << tet_data.dataOnEntities[MBEDGE][ee].getN("TEST_FIELD")
456 << std::endl;
457 std::cout << tet_data.dataOnEntities[MBEDGE][ee].getDiffN("TEST_FIELD")
458 << std::endl;
459 sum +=
460 sum_matrix(tet_data.dataOnEntities[MBEDGE][ee].getN("TEST_FIELD"));
461 diff_sum += sum_matrix(
462 tet_data.dataOnEntities[MBEDGE][ee].getDiffN("TEST_FIELD"));
463 }
464 std::cout << "Faces\n";
465 for (int ff = 0; ff < 4; ff++) {
466 std::cout << tet_data.dataOnEntities[MBTRI][ff].getN("TEST_FIELD")
467 << std::endl;
468 std::cout << tet_data.dataOnEntities[MBTRI][ff].getDiffN("TEST_FIELD")
469 << std::endl;
470 sum +=
471 sum_matrix(tet_data.dataOnEntities[MBTRI][ff].getN("TEST_FIELD"));
472 diff_sum += sum_matrix(
473 tet_data.dataOnEntities[MBTRI][ff].getDiffN("TEST_FIELD"));
474 }
475 std::cout << "Tets\n";
476 std::cout << tet_data.dataOnEntities[MBTET][0].getN("TEST_FIELD")
477 << std::endl;
478 std::cout << tet_data.dataOnEntities[MBTET][0].getDiffN("TEST_FIELD")
479 << std::endl;
480 sum += sum_matrix(tet_data.dataOnEntities[MBTET][0].getN("TEST_FIELD"));
481 diff_sum +=
482 sum_matrix(tet_data.dataOnEntities[MBTET][0].getDiffN("TEST_FIELD"));
483 std::cout << "sum " << sum << std::endl;
484 std::cout << "diff_sum " << diff_sum << std::endl;
485 if (fabs(4.38395 - sum) > eps) {
486 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
487 }
488 if (fabs(diff_sum) > eps) {
489 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
490 }
491 }
492
493 if (choice_value == HDIVTET_AINSWORTH) {
494 ierr = TetPolynomialBase().getValue(
495 pts_tet, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
496 tet_data, HDIV, CONTINUOUS, AINSWORTH_LEGENDRE_BASE)));
497 CHKERRG(ierr);
498 double sum = 0, diff_sum = 0;
499 std::cout << "Faces\n";
500 for (int ff = 0; ff < 4; ff++) {
501 std::cout << tet_data.dataOnEntities[MBTRI][ff].getN(
502 AINSWORTH_LEGENDRE_BASE)
503 << std::endl;
504 std::cout << tet_data.dataOnEntities[MBTRI][ff].getDiffN(
505 AINSWORTH_LEGENDRE_BASE)
506 << std::endl;
507 sum += sum_matrix(
508 tet_data.dataOnEntities[MBTRI][ff].getN(AINSWORTH_LEGENDRE_BASE));
509 diff_sum += sum_matrix(tet_data.dataOnEntities[MBTRI][ff].getDiffN(
510 AINSWORTH_LEGENDRE_BASE));
511 }
512 std::cout << "Tets\n";
513 std::cout << tet_data.dataOnEntities[MBTET][0].getN(
514 AINSWORTH_LEGENDRE_BASE)
515 << std::endl;
516 std::cout << tet_data.dataOnEntities[MBTET][0].getDiffN(
517 AINSWORTH_LEGENDRE_BASE)
518 << std::endl;
519 sum += sum_matrix(
520 tet_data.dataOnEntities[MBTET][0].getN(AINSWORTH_LEGENDRE_BASE));
521 diff_sum += sum_matrix(
522 tet_data.dataOnEntities[MBTET][0].getDiffN(AINSWORTH_LEGENDRE_BASE));
523 std::cout << "sum " << 1e8 * sum << std::endl;
524 std::cout << "diff_sum " << 1e8 * diff_sum << std::endl;
525 if (fabs(0.178713 - sum) > eps) {
526 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
527 }
528 if (fabs(31.6048 - diff_sum) > eps) {
529 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
530 }
531 }
532
533 if (choice_value == HDIVTET_DEMKOWICZ) {
534 ierr = TetPolynomialBase().getValue(
535 pts_tet, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
536 tet_data, HDIV, CONTINUOUS, DEMKOWICZ_JACOBI_BASE)));
537 CHKERRG(ierr);
538 double sum = 0, diff_sum = 0;
539 std::cout << "Faces\n";
540 for (int ff = 0; ff < 4; ff++) {
541 std::cout << tet_data.dataOnEntities[MBTRI][ff].getN(
542 DEMKOWICZ_JACOBI_BASE)
543 << std::endl;
544 std::cout << tet_data.dataOnEntities[MBTRI][ff].getDiffN(
545 DEMKOWICZ_JACOBI_BASE)
546 << std::endl;
547 sum += sum_matrix(
548 tet_data.dataOnEntities[MBTRI][ff].getN(DEMKOWICZ_JACOBI_BASE));
549 diff_sum += sum_matrix(
550 tet_data.dataOnEntities[MBTRI][ff].getDiffN(DEMKOWICZ_JACOBI_BASE));
551 }
552 std::cout << "Tets\n";
553 std::cout << tet_data.dataOnEntities[MBTET][0].getN(DEMKOWICZ_JACOBI_BASE)
554 << std::endl;
555 std::cout << tet_data.dataOnEntities[MBTET][0].getDiffN(
556 DEMKOWICZ_JACOBI_BASE)
557 << std::endl;
558 sum += sum_matrix(
559 tet_data.dataOnEntities[MBTET][0].getN(DEMKOWICZ_JACOBI_BASE));
560 diff_sum += sum_matrix(
561 tet_data.dataOnEntities[MBTET][0].getDiffN(DEMKOWICZ_JACOBI_BASE));
562 std::cout << "sum " << sum << std::endl;
563 std::cout << "diff_sum " << diff_sum << std::endl;
564 const double expected_result = -2.70651;
565 const double expected_diff_result = 289.421;
566 if (fabs((expected_result - sum) / expected_result) > eps) {
567 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
568 }
569 if (fabs((expected_diff_result - diff_sum) / expected_diff_result) >
570 eps) {
571 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
572 }
573 }
574
575 if (choice_value == HCURLTET_AINSWORTH) {
576 ierr = TetPolynomialBase().getValue(
577 pts_tet, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
578 tet_data, HCURL, CONTINUOUS, AINSWORTH_LEGENDRE_BASE)));
579 CHKERRG(ierr);
580 double sum = 0, diff_sum = 0;
581 std::cout << "Edges\n";
582 for (int ee = 0; ee < 6; ee++) {
583 std::cout << tet_data.dataOnEntities[MBEDGE][ee].getN(
584 AINSWORTH_LEGENDRE_BASE)
585 << std::endl;
586 std::cout << tet_data.dataOnEntities[MBEDGE][ee].getDiffN(
587 AINSWORTH_LEGENDRE_BASE)
588 << std::endl;
589 sum += sum_matrix(
590 tet_data.dataOnEntities[MBEDGE][ee].getN(AINSWORTH_LEGENDRE_BASE));
591 diff_sum += sum_matrix(tet_data.dataOnEntities[MBEDGE][ee].getDiffN(
592 AINSWORTH_LEGENDRE_BASE));
593 }
594 std::cout << "Faces\n";
595 for (int ff = 0; ff < 4; ff++) {
596 std::cout << tet_data.dataOnEntities[MBTRI][ff].getN(
597 AINSWORTH_LEGENDRE_BASE)
598 << std::endl;
599 std::cout << tet_data.dataOnEntities[MBTRI][ff].getDiffN(
600 AINSWORTH_LEGENDRE_BASE)
601 << std::endl;
602 sum += sum_matrix(
603 tet_data.dataOnEntities[MBTRI][ff].getN(AINSWORTH_LEGENDRE_BASE));
604 diff_sum += sum_matrix(tet_data.dataOnEntities[MBTRI][ff].getDiffN(
605 AINSWORTH_LEGENDRE_BASE));
606 }
607 std::cout << "Tets\n";
608 std::cout << tet_data.dataOnEntities[MBTET][0].getN(
609 AINSWORTH_LEGENDRE_BASE)
610 << std::endl;
611 std::cout << tet_data.dataOnEntities[MBTET][0].getDiffN(
612 AINSWORTH_LEGENDRE_BASE)
613 << std::endl;
614 sum += sum_matrix(
615 tet_data.dataOnEntities[MBTET][0].getN(AINSWORTH_LEGENDRE_BASE));
616 diff_sum += sum_matrix(
617 tet_data.dataOnEntities[MBTET][0].getDiffN(AINSWORTH_LEGENDRE_BASE));
618 std::cout << "sum " << sum << std::endl;
619 std::cout << "diff_sum " << diff_sum << std::endl;
620 if (fabs(-1.7798 - sum) > eps) {
621 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
622 }
623 if (fabs(-67.1793 - diff_sum) > eps) {
624 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
625 }
626 }
627
628 if (choice_value == HCURLTET_DEMKOWICZ) {
629 CHKERR TetPolynomialBase().getValue(
630 pts_tet, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
631 tet_data, HCURL, CONTINUOUS, DEMKOWICZ_JACOBI_BASE)));
632 double sum = 0, diff_sum = 0;
633 std::cout << "Edges\n";
634 for (int ee = 0; ee < 6; ee++) {
635 std::cout << tet_data.dataOnEntities[MBEDGE][ee].getN(
636 DEMKOWICZ_JACOBI_BASE)
637 << std::endl;
638 std::cout << tet_data.dataOnEntities[MBEDGE][ee].getDiffN(
639 DEMKOWICZ_JACOBI_BASE)
640 << std::endl;
641 sum += sum_matrix(
642 tet_data.dataOnEntities[MBEDGE][ee].getN(DEMKOWICZ_JACOBI_BASE));
643 diff_sum += sum_matrix(tet_data.dataOnEntities[MBEDGE][ee].getDiffN(
644 DEMKOWICZ_JACOBI_BASE));
645 }
646 std::cout << "Faces\n";
647 for (int ff = 0; ff < 4; ff++) {
648 std::cout << tet_data.dataOnEntities[MBTRI][ff].getN(
649 DEMKOWICZ_JACOBI_BASE)
650 << std::endl;
651 std::cout << tet_data.dataOnEntities[MBTRI][ff].getDiffN(
652 DEMKOWICZ_JACOBI_BASE)
653 << std::endl;
654 sum += sum_matrix(
655 tet_data.dataOnEntities[MBTRI][ff].getN(DEMKOWICZ_JACOBI_BASE));
656 diff_sum += sum_matrix(
657 tet_data.dataOnEntities[MBTRI][ff].getDiffN(DEMKOWICZ_JACOBI_BASE));
658 }
659 std::cout << "Tets\n";
660 std::cout << tet_data.dataOnEntities[MBTET][0].getN(DEMKOWICZ_JACOBI_BASE)
661 << std::endl;
662 std::cout << tet_data.dataOnEntities[MBTET][0].getDiffN(
663 DEMKOWICZ_JACOBI_BASE)
664 << std::endl;
665 sum += sum_matrix(
666 tet_data.dataOnEntities[MBTET][0].getN(DEMKOWICZ_JACOBI_BASE));
667 diff_sum += sum_matrix(
668 tet_data.dataOnEntities[MBTET][0].getDiffN(DEMKOWICZ_JACOBI_BASE));
669 std::cout << "sum " << sum << std::endl;
670 std::cout << "diff_sum " << diff_sum << std::endl;
671 if (fabs(7.35513 - sum) > eps) {
672 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
673 }
674 if (fabs(62.4549 - diff_sum) > eps) {
675 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
676 }
677 }
678
679 if (choice_value == L2TET) {
680 ierr = TetPolynomialBase().getValue(
681 pts_tet, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
682 tet_data, L2, CONTINUOUS, AINSWORTH_LEGENDRE_BASE)));
683 CHKERRG(ierr);
684 double sum = 0, diff_sum = 0;
685 std::cout << "Tets\n";
686 std::cout << tet_data.dataOnEntities[MBTET][0].getN(
687 AINSWORTH_LEGENDRE_BASE)
688 << std::endl;
689 std::cout << tet_data.dataOnEntities[MBTET][0].getDiffN(
690 AINSWORTH_LEGENDRE_BASE)
691 << std::endl;
692 sum += sum_matrix(
693 tet_data.dataOnEntities[MBTET][0].getN(AINSWORTH_LEGENDRE_BASE));
694 diff_sum += sum_matrix(
695 tet_data.dataOnEntities[MBTET][0].getDiffN(AINSWORTH_LEGENDRE_BASE));
696 std::cout << "sum " << sum << std::endl;
697 std::cout << "diff_sum " << diff_sum << std::endl;
698 if (fabs(3.60352 - sum) > eps) {
699 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
700 }
701 if (fabs(-36.9994 - diff_sum) > eps) {
702 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
703 }
704 }
705
706 EntitiesFieldData tri_data(MBTRI);
707 for (int type = MBVERTEX; type != MBMAXTYPE; type++) {
708 tri_data.spacesOnEntities[type].set(L2);
709 tri_data.spacesOnEntities[type].set(H1);
710 tri_data.spacesOnEntities[type].set(HDIV);
711 tri_data.spacesOnEntities[type].set(HCURL);
712 }
713 tri_data.dataOnEntities[MBVERTEX].resize(1);
714 tri_data.dataOnEntities[MBVERTEX][0].getOrder() = 1;
715 tri_data.dataOnEntities[MBEDGE].resize(3);
716 for (int ee = 0; ee < 3; ee++) {
717 tri_data.dataOnEntities[MBEDGE][ee].getOrder() = 3;
718 tri_data.dataOnEntities[MBEDGE][ee].getSense() = 1;
719 }
720 tri_data.dataOnEntities[MBTRI].resize(1);
721 tri_data.dataOnEntities[MBTRI][0].getOrder() = 4;
722 tri_data.dataOnEntities[MBVERTEX][0].getBBNodeOrder().resize(3, false);
723 std::fill(tri_data.dataOnEntities[MBVERTEX][0].getBBNodeOrder().begin(),
724 tri_data.dataOnEntities[MBVERTEX][0].getBBNodeOrder().end(), 4);
725
726 MatrixDouble pts_tri;
727 int tri_rule = 2;
728 nb_gauss_pts = QUAD_2D_TABLE[tri_rule]->npoints;
729 pts_tri.resize(2, nb_gauss_pts, false);
730 cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[tri_rule]->points[1], 3,
731 &pts_tri(0, 0), 1);
732 cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[tri_rule]->points[2], 3,
733 &pts_tri(1, 0), 1);
734 tri_data.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 3,
735 false);
736 {
737 double *shape_ptr =
738 &*tri_data.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
739 cblas_dcopy(3 * nb_gauss_pts, QUAD_2D_TABLE[tri_rule]->points, 1,
740 shape_ptr, 1);
741 }
742 tri_data.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).resize(3, 2, false);
743 ierr = ShapeDiffMBTRI(
744 &*tri_data.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).data().begin());
745 CHKERRG(ierr);
746
747 if (choice_value == H1TRI_AINSWORTH) {
748 CHKERR TriPolynomialBase().getValue(
749 pts_tri, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
750 tri_data, H1, CONTINUOUS, AINSWORTH_LEGENDRE_BASE, NOBASE)));
751 CHKERR check_same_matrix(
752 tri_data.dataOnEntities[MBVERTEX][0].getN(NOBASE),
753 tri_data.dataOnEntities[MBVERTEX][0].getN(AINSWORTH_LEGENDRE_BASE),
754 eps);
755
756 double sum = 0, diff_sum = 0;
757 std::cout << "Edges\n";
758 for (int ee = 0; ee < 3; ee++) {
759 std::cout << tri_data.dataOnEntities[MBEDGE][ee].getN(
760 AINSWORTH_LEGENDRE_BASE)
761 << std::endl;
762 std::cout << tri_data.dataOnEntities[MBEDGE][ee].getDiffN(
763 AINSWORTH_LEGENDRE_BASE)
764 << std::endl;
765 sum += sum_matrix(
766 tri_data.dataOnEntities[MBEDGE][ee].getN(AINSWORTH_LEGENDRE_BASE));
767 diff_sum += sum_matrix(tri_data.dataOnEntities[MBEDGE][ee].getDiffN(
768 AINSWORTH_LEGENDRE_BASE));
769 }
770 std::cout << "Face\n";
771 std::cout << tri_data.dataOnEntities[MBTRI][0].getN(
772 AINSWORTH_LEGENDRE_BASE)
773 << std::endl;
774 std::cout << tri_data.dataOnEntities[MBTRI][0].getDiffN(
775 AINSWORTH_LEGENDRE_BASE)
776 << std::endl;
777 sum += sum_matrix(
778 tri_data.dataOnEntities[MBTRI][0].getN(AINSWORTH_LEGENDRE_BASE));
779 diff_sum += sum_matrix(
780 tri_data.dataOnEntities[MBTRI][0].getDiffN(AINSWORTH_LEGENDRE_BASE));
781 std::cout << "sum " << sum << std::endl;
782 std::cout << "diff_sum " << diff_sum << std::endl;
783 if (fabs(0.805556 - sum) > eps)
784 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
785
786 if (fabs(0.0833333 - diff_sum) > eps)
787 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
788 }
789
790 if (choice_value == H1TRI_BERNSTEIN_BEZIER) {
791 CHKERR TriPolynomialBase().getValue(
792 pts_tri, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
793 tri_data, "TET_FIELD", H1, CONTINUOUS,
794 AINSWORTH_BERNSTEIN_BEZIER_BASE, NOBASE)));
795 if (tri_data.dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE).get() ==
796 tri_data.dataOnEntities[MBVERTEX][0]
797 .getNSharedPtr(AINSWORTH_BERNSTEIN_BEZIER_BASE)
798 .get())
799 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
800 "Pointers should be diffrent");
801
802 double sum = 0, diff_sum = 0;
803 std::cout << "Vertex\n";
804 std::cout << tri_data.dataOnEntities[MBVERTEX][0].getN("TET_FIELD")
805 << std::endl;
806 std::cout << tri_data.dataOnEntities[MBVERTEX][0].getDiffN("TET_FIELD")
807 << std::endl;
808 sum += sum_matrix(tri_data.dataOnEntities[MBVERTEX][0].getN("TET_FIELD"));
809 diff_sum += sum_matrix(
810 tri_data.dataOnEntities[MBVERTEX][0].getDiffN("TET_FIELD"));
811 std::cout << "Edges\n";
812 for (int ee = 0; ee < 3; ee++) {
813 std::cout << tri_data.dataOnEntities[MBEDGE][ee].getN("TET_FIELD")
814 << std::endl;
815 std::cout << tri_data.dataOnEntities[MBEDGE][ee].getDiffN("TET_FIELD")
816 << std::endl;
817 sum +=
818 sum_matrix(tri_data.dataOnEntities[MBEDGE][ee].getN("TET_FIELD"));
819 diff_sum += sum_matrix(
820 tri_data.dataOnEntities[MBEDGE][ee].getDiffN("TET_FIELD"));
821 }
822 std::cout << "Face\n";
823 std::cout << tri_data.dataOnEntities[MBTRI][0].getN("TET_FIELD")
824 << std::endl;
825 std::cout << tri_data.dataOnEntities[MBTRI][0].getDiffN("TET_FIELD")
826 << std::endl;
827 sum += sum_matrix(tri_data.dataOnEntities[MBTRI][0].getN("TET_FIELD"));
828 diff_sum +=
829 sum_matrix(tri_data.dataOnEntities[MBTRI][0].getDiffN("TET_FIELD"));
830 std::cout << "sum " << sum << std::endl;
831 std::cout << "diff_sum " << diff_sum << std::endl;
832 if (std::abs(3.01389 - sum) > eps)
833 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
834
835 if (std::abs(diff_sum) > eps)
836 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
837 "wrong result %3.4e != %3.4e", 0., diff_sum);
838 }
839
840 if (choice_value == HDIVTRI_AINSWORTH) {
841 CHKERR TriPolynomialBase().getValue(
842 pts_tri, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
843 tri_data, HDIV, CONTINUOUS, AINSWORTH_LEGENDRE_BASE, NOBASE)));
844 CHKERR check_same_matrix(
845 tri_data.dataOnEntities[MBVERTEX][0].getN(NOBASE),
846 tri_data.dataOnEntities[MBVERTEX][0].getN(AINSWORTH_LEGENDRE_BASE),
847 eps);
848
849 double sum = 0;
850 std::cout << "Face\n";
851 std::cout << tri_data.dataOnEntities[MBTRI][0].getN(
852 AINSWORTH_LEGENDRE_BASE)
853 << std::endl;
854 sum += sum_matrix(
855 tri_data.dataOnEntities[MBTRI][0].getN(AINSWORTH_LEGENDRE_BASE));
856 std::cout << "sum " << sum << std::endl;
857 if (fabs(2.21296 - sum) > eps)
858 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
859 }
860
861 if (choice_value == HDIVTRI_DEMKOWICZ) {
862 ierr = TriPolynomialBase().getValue(
863 pts_tri, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
864 tri_data, HDIV, CONTINUOUS, DEMKOWICZ_JACOBI_BASE, NOBASE)));
865 CHKERRG(ierr);
866 ierr = check_same_matrix(
867 tri_data.dataOnEntities[MBVERTEX][0].getN(NOBASE),
868 tri_data.dataOnEntities[MBVERTEX][0].getN(DEMKOWICZ_JACOBI_BASE),
869 eps);
870 CHKERRG(ierr);
871 double sum = 0;
872 std::cout << "Face\n";
873 std::cout << tri_data.dataOnEntities[MBTRI][0].getN(DEMKOWICZ_JACOBI_BASE)
874 << std::endl;
875 sum += sum_matrix(
876 tri_data.dataOnEntities[MBTRI][0].getN(DEMKOWICZ_JACOBI_BASE));
877 std::cout << "sum " << sum << std::endl;
878 const double expected_result = 28.25;
879 if (fabs((expected_result - sum) / expected_result) > eps) {
880 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
881 }
882 }
883
884 if (choice_value == HCURLTRI_AINSWORTH) {
885 ierr = TriPolynomialBase().getValue(
886 pts_tri, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
887 tri_data, HCURL, CONTINUOUS, AINSWORTH_LEGENDRE_BASE, NOBASE)));
888 CHKERRG(ierr);
889 ierr = check_same_matrix(
890 tri_data.dataOnEntities[MBVERTEX][0].getN(NOBASE),
891 tri_data.dataOnEntities[MBVERTEX][0].getN(AINSWORTH_LEGENDRE_BASE),
892 eps);
893 CHKERRG(ierr);
894 double sum = 0; //,diff_sum = 0;
895 std::cout << "Edges\n";
896 for (int ee = 0; ee < 3; ee++) {
897 std::cout << tri_data.dataOnEntities[MBEDGE][ee].getN(
898 AINSWORTH_LEGENDRE_BASE)
899 << std::endl;
900 sum += sum_matrix(
901 tri_data.dataOnEntities[MBEDGE][ee].getN(AINSWORTH_LEGENDRE_BASE));
902 }
903 std::cout << "Face\n";
904 std::cout << tri_data.dataOnEntities[MBTRI][0].getN(
905 AINSWORTH_LEGENDRE_BASE)
906 << std::endl;
907 sum += sum_matrix(
908 tri_data.dataOnEntities[MBTRI][0].getN(AINSWORTH_LEGENDRE_BASE));
909 std::cout << "sum " << sum << std::endl;
910 if (fabs(0.333333 - sum) > eps) {
911 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
912 }
913 }
914
915 if (choice_value == HCURLTRI_DEMKOWICZ) {
916 CHKERR TriPolynomialBase().getValue(
917 pts_tri, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
918 tri_data, HCURL, CONTINUOUS, DEMKOWICZ_JACOBI_BASE, NOBASE)));
919 ierr = check_same_matrix(
920 tri_data.dataOnEntities[MBVERTEX][0].getN(NOBASE),
921 tri_data.dataOnEntities[MBVERTEX][0].getN(DEMKOWICZ_JACOBI_BASE),
922 eps);
923 CHKERRG(ierr);
924 double sum = 0; //,diff_sum = 0;
925 std::cout << "Edges\n";
926 for (int ee = 0; ee < 3; ee++) {
927 std::cout << tri_data.dataOnEntities[MBEDGE][ee].getN(
928 DEMKOWICZ_JACOBI_BASE)
929 << std::endl;
930 sum += sum_matrix(
931 tri_data.dataOnEntities[MBEDGE][ee].getN(DEMKOWICZ_JACOBI_BASE));
932 }
933 std::cout << "Face\n";
934 std::cout << tri_data.dataOnEntities[MBTRI][0].getN(DEMKOWICZ_JACOBI_BASE)
935 << std::endl;
936 sum += sum_matrix(
937 tri_data.dataOnEntities[MBTRI][0].getN(DEMKOWICZ_JACOBI_BASE));
938 std::cout << "sum " << sum << std::endl;
939 if (fabs(1.04591 - sum) > eps) {
940 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
941 }
942 }
943
944 if (choice_value == L2TRI) {
945 ierr = TriPolynomialBase().getValue(
946 pts_tri, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
947 tri_data, L2, CONTINUOUS, AINSWORTH_LEGENDRE_BASE, NOBASE)));
948 CHKERRG(ierr);
949 ierr = check_same_matrix(
950 tri_data.dataOnEntities[MBVERTEX][0].getN(NOBASE),
951 tri_data.dataOnEntities[MBVERTEX][0].getN(AINSWORTH_LEGENDRE_BASE),
952 eps);
953 CHKERRG(ierr);
954 double sum = 0;
955 std::cout << "Face\n";
956 std::cout << tri_data.dataOnEntities[MBTRI][0].getN(
957 AINSWORTH_LEGENDRE_BASE)
958 << std::endl;
959 sum += sum_matrix(
960 tri_data.dataOnEntities[MBTRI][0].getN(AINSWORTH_LEGENDRE_BASE));
961 std::cout << "sum " << sum << std::endl;
962 if (fabs(0.671875 - sum) > eps) {
963 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
964 }
965 }
966
967 EntitiesFieldData edge_data(MBTRI);
968 for (int type = MBVERTEX; type != MBMAXTYPE; type++) {
969 edge_data.spacesOnEntities[type].set(L2);
970 edge_data.spacesOnEntities[type].set(H1);
971 edge_data.spacesOnEntities[type].set(HDIV);
972 edge_data.spacesOnEntities[type].set(HCURL);
973 }
974 edge_data.dataOnEntities[MBVERTEX].resize(1);
975 edge_data.dataOnEntities[MBVERTEX][0].getOrder() = 1;
976 edge_data.dataOnEntities[MBEDGE].resize(1);
977 edge_data.dataOnEntities[MBEDGE][0].getOrder() = 4;
978 edge_data.dataOnEntities[MBEDGE][0].getSense() = 1;
979 edge_data.dataOnEntities[MBVERTEX][0].getBBNodeOrder().resize(2, false);
980 std::fill(edge_data.dataOnEntities[MBVERTEX][0].getBBNodeOrder().begin(),
981 edge_data.dataOnEntities[MBVERTEX][0].getBBNodeOrder().end(), 4);
982
983 MatrixDouble pts_edge;
984 int edge_rule = 6;
985 nb_gauss_pts = QUAD_1D_TABLE[edge_rule]->npoints;
986 pts_edge.resize(2, nb_gauss_pts, false);
987 cblas_dcopy(nb_gauss_pts, &QUAD_1D_TABLE[edge_rule]->points[1], 2,
988 &pts_edge(0, 0), 1);
989 cblas_dcopy(nb_gauss_pts, QUAD_1D_TABLE[edge_rule]->weights, 1,
990 &pts_edge(1, 0), 1);
991
992 edge_data.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 2,
993 false);
994 {
995 double *shape_ptr =
996 &*edge_data.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
997 cblas_dcopy(2 * nb_gauss_pts, QUAD_1D_TABLE[edge_rule]->points, 1,
998 shape_ptr, 1);
999 }
1000
1001 if (choice_value == H1EDGE_AINSWORTH) {
1002 CHKERR EdgePolynomialBase().getValue(
1003 pts_edge, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
1004 edge_data, H1, CONTINUOUS, AINSWORTH_LEGENDRE_BASE, NOBASE)));
1005 ierr = check_same_matrix(
1006 edge_data.dataOnEntities[MBVERTEX][0].getN(NOBASE),
1007 edge_data.dataOnEntities[MBVERTEX][0].getN(AINSWORTH_LEGENDRE_BASE),
1008 eps);
1009 CHKERRG(ierr);
1010 double sum = 0, diff_sum = 0;
1011 std::cout << "Edge\n";
1012 std::cout << edge_data.dataOnEntities[MBEDGE][0].getN(
1013 AINSWORTH_LEGENDRE_BASE)
1014 << std::endl;
1015 std::cout << edge_data.dataOnEntities[MBEDGE][0].getDiffN(
1016 AINSWORTH_LEGENDRE_BASE)
1017 << std::endl;
1018 sum += sum_matrix(
1019 edge_data.dataOnEntities[MBEDGE][0].getN(AINSWORTH_LEGENDRE_BASE));
1020 diff_sum += sum_matrix(edge_data.dataOnEntities[MBEDGE][0].getDiffN(
1021 AINSWORTH_LEGENDRE_BASE));
1022 std::cout << "sum " << sum << std::endl;
1023 std::cout << "diff sum " << diff_sum << std::endl;
1024 if (std::abs(0.506122 - sum) > eps)
1025 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
1026
1027 if (std::abs(-0.571429 - diff_sum) > eps)
1028 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
1029 }
1030
1031 if (choice_value == H1EDGE_BERNSTEIN_BEZIER) {
1032 CHKERR EdgePolynomialBase().getValue(
1033 pts_edge, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
1034 edge_data, "TET_FIELD", H1, CONTINUOUS,
1035 AINSWORTH_BERNSTEIN_BEZIER_BASE, NOBASE)));
1036 if (edge_data.dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE).get() ==
1037 edge_data.dataOnEntities[MBVERTEX][0]
1038 .getNSharedPtr(AINSWORTH_BERNSTEIN_BEZIER_BASE)
1039 .get())
1040 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1041 "Should be diffrent pointers");
1042
1043 double sum = 0, diff_sum = 0;
1044 std::cout << "Vertex\n";
1045 std::cout << edge_data.dataOnEntities[MBVERTEX][0].getN("TET_FIELD")
1046 << std::endl;
1047 std::cout << edge_data.dataOnEntities[MBVERTEX][0].getDiffN("TET_FIELD")
1048 << std::endl;
1049 std::cout << "Edge\n";
1050 std::cout << edge_data.dataOnEntities[MBEDGE][0].getN("TET_FIELD")
1051 << std::endl;
1052 std::cout << edge_data.dataOnEntities[MBEDGE][0].getDiffN("TET_FIELD")
1053 << std::endl;
1054 sum +=
1055 sum_matrix(edge_data.dataOnEntities[MBVERTEX][0].getN("TET_FIELD"));
1056 diff_sum += sum_matrix(
1057 edge_data.dataOnEntities[MBVERTEX][0].getDiffN("TET_FIELD"));
1058 sum += sum_matrix(edge_data.dataOnEntities[MBEDGE][0].getN("TET_FIELD"));
1059 diff_sum +=
1060 sum_matrix(edge_data.dataOnEntities[MBEDGE][0].getDiffN("TET_FIELD"));
1061 std::cout << "sum " << sum << std::endl;
1062 std::cout << "diff sum " << diff_sum << std::endl;
1063 if (std::abs(4 - sum) > eps)
1064 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
1065
1066 if (std::abs(diff_sum) > eps)
1067 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
1068 }
1069
1070 if (choice_value == HCURLEDGE_AINSWORTH) {
1071 CHKERR EdgePolynomialBase().getValue(
1072 pts_edge,
1073 boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
1074 edge_data, HCURL, CONTINUOUS, AINSWORTH_LEGENDRE_BASE, NOBASE)));
1075 ierr = check_same_matrix(
1076 edge_data.dataOnEntities[MBVERTEX][0].getN(NOBASE),
1077 edge_data.dataOnEntities[MBVERTEX][0].getN(AINSWORTH_LEGENDRE_BASE),
1078 eps);
1079 CHKERRG(ierr);
1080 std::cout << edge_data.dataOnEntities[MBEDGE][0].getN(
1081 AINSWORTH_LEGENDRE_BASE)
1082 << std::endl;
1083 double sum = 0;
1084 sum += sum_matrix(
1085 edge_data.dataOnEntities[MBEDGE][0].getN(AINSWORTH_LEGENDRE_BASE));
1086 std::cout.precision(std::numeric_limits<double>::max_digits10 - 1);
1087 std::cout << "sum " << sum << std::endl;
1088 if (fabs(-4.174149659863944 - sum) > eps) {
1089 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
1090 }
1091 }
1092
1093 if (choice_value == HCURLEDGE_DEMKOWICZ) {
1094 CHKERR EdgePolynomialBase().getValue(
1095 pts_edge, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
1096 edge_data, HCURL, CONTINUOUS, DEMKOWICZ_JACOBI_BASE, NOBASE)));
1097 ierr = check_same_matrix(
1098 edge_data.dataOnEntities[MBVERTEX][0].getN(NOBASE),
1099 edge_data.dataOnEntities[MBVERTEX][0].getN(DEMKOWICZ_JACOBI_BASE),
1100 eps);
1101 CHKERRG(ierr);
1102 std::cout << edge_data.dataOnEntities[MBEDGE][0].getN(
1103 DEMKOWICZ_JACOBI_BASE)
1104 << std::endl;
1105 double sum = sum_matrix(
1106 edge_data.dataOnEntities[MBEDGE][0].getN(DEMKOWICZ_JACOBI_BASE));
1107 std::cout.precision(std::numeric_limits<double>::max_digits10 - 1);
1108 std::cout << "sum " << sum << std::endl;
1109 if (std::fabs(4.571428571428571 - sum) > eps) {
1110 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
1111 }
1112 }
1113
1114 if (choice_value == L2EDGE) {
1115 CHKERR EdgePolynomialBase().getValue(
1116 pts_edge,
1117 boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
1118 edge_data, L2, CONTINUOUS, DEMKOWICZ_JACOBI_BASE, NOBASE)));
1119 CHKERR EdgePolynomialBase().getValue(
1120 pts_edge,
1121 boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
1122 edge_data, L2, CONTINUOUS, AINSWORTH_LEGENDRE_BASE, NOBASE)));
1123 CHKERR EdgePolynomialBase().getValue(
1124 pts_edge,
1125 boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
1126 edge_data, L2, CONTINUOUS, AINSWORTH_LOBATTO_BASE, NOBASE)));
1127
1128 MatrixDouble diff_n =
1129 edge_data.dataOnEntities[MBEDGE][0].getN(DEMKOWICZ_JACOBI_BASE) -
1130 edge_data.dataOnEntities[MBEDGE][0].getN(AINSWORTH_LEGENDRE_BASE);
1131 double sum = sum_matrix(diff_n);
1132 std::cout << "sum " << sum << std::endl;
1133 if (std::fabs(sum) > eps) {
1134 SETERRQ(
1135 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1136 "Difference between ainsworth and demkowicz base should be zero");
1137 }
1138
1139 auto order = edge_data.dataOnEntities[MBEDGE][0].getOrder();
1140 auto &base =
1141 edge_data.dataOnEntities[MBEDGE][0].getN(AINSWORTH_LEGENDRE_BASE);
1142
1143 if (base.size1() != pts_edge.size2())
1144 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1145 "wrong number of integration points");
1146 if (base.size2() != NBEDGE_L2(order))
1147 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1148 "wrong number of base functions");
1149
1150 double ort_sum = 0;
1151 for (auto gg = 0; gg != pts_edge.size2(); ++gg) {
1152 double w = pts_edge(1, gg);
1153 for (int i = 0; i != NBEDGE_L2(order); ++i) {
1154 for (int j = 0; j != NBEDGE_L2(order); ++j) {
1155 if (i != j) {
1156 ort_sum += base(gg, i) * base(gg, j) * w;
1157 }
1158 }
1159 }
1160 }
1161 std::cout << "ort sum " << ort_sum << std::endl;
1162 if (std::fabs(ort_sum) > eps) {
1163 SETERRQ(
1164 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1165 "check orthogonality of base");
1166 }
1167
1168 auto &diff_base =
1169 edge_data.dataOnEntities[MBEDGE][0].getDiffN(AINSWORTH_LOBATTO_BASE);
1170 double ort_diff_sum = 0;
1171 for (auto gg = 0; gg != pts_edge.size2(); ++gg) {
1172 double w = pts_edge(1, gg);
1173 for (int i = 2; i != NBEDGE_L2(order); ++i) {
1174 for (int j = 2; j != NBEDGE_L2(order); ++j) {
1175 if (i != j) {
1176 ort_diff_sum += diff_base(gg, i) * diff_base(gg, j) * w;
1177 }
1178 }
1179 }
1180 }
1181 std::cout << "ort diff sum " << ort_diff_sum << std::endl;
1182 if (std::fabs(ort_diff_sum) > eps) {
1183 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1184 "check orthogonality of diff lobatto base");
1185 }
1186
1187 }
1188
1189 EntitiesFieldData quad_data(MBQUAD);
1190 for (int type = MBVERTEX; type != MBMAXTYPE; type++) {
1191 quad_data.spacesOnEntities[type].set(H1);
1192 }
1193 quad_data.dataOnEntities[MBVERTEX].resize(1);
1194 quad_data.dataOnEntities[MBVERTEX][0].getOrder() = 1;
1195 quad_data.dataOnEntities[MBEDGE].resize(4);
1196 for (int ee = 0; ee < 4; ee++) {
1197 quad_data.dataOnEntities[MBEDGE][ee].getOrder() = 4;
1198 quad_data.dataOnEntities[MBEDGE][ee].getSense() = 1;
1199 }
1200 quad_data.dataOnEntities[MBQUAD].resize(1);
1201 quad_data.dataOnEntities[MBQUAD][0].getOrder() = 6;
1202
1203 MatrixDouble pts_quad;
1204 int rule_ksi = 6;
1205 int rule_eta = 4;
1206 CHKERR Tools::outerProductOfEdgeIntegrationPtsForQuad(pts_quad, rule_ksi,
1207 rule_eta);
1208 nb_gauss_pts = pts_quad.size2();
1209 quad_data.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4,
1210 false);
1211 quad_data.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).resize(nb_gauss_pts,
1212 8, false);
1213 for (int i = 0; i < nb_gauss_pts; ++i) {
1214 quad_data.dataOnEntities[MBVERTEX][0].getN(NOBASE)(i, 0) =
1215 N_MBQUAD0(pts_quad(0, i), pts_quad(1, i));
1216 quad_data.dataOnEntities[MBVERTEX][0].getN(NOBASE)(i, 1) =
1217 N_MBQUAD1(pts_quad(0, i), pts_quad(1, i));
1218 quad_data.dataOnEntities[MBVERTEX][0].getN(NOBASE)(i, 2) =
1219 N_MBQUAD2(pts_quad(0, i), pts_quad(1, i));
1220 quad_data.dataOnEntities[MBVERTEX][0].getN(NOBASE)(i, 3) =
1221 N_MBQUAD3(pts_quad(0, i), pts_quad(1, i));
1222
1223 quad_data.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE)(i, 0) =
1224 diffN_MBQUAD0x(pts_quad(1, i));
1225 quad_data.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE)(i, 1) =
1226 diffN_MBQUAD0y(pts_quad(0, i));
1227 quad_data.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE)(i, 2) =
1228 diffN_MBQUAD1x(pts_quad(1, i));
1229 quad_data.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE)(i, 3) =
1230 diffN_MBQUAD1y(pts_quad(0, i));
1231 quad_data.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE)(i, 4) =
1232 diffN_MBQUAD2x(pts_quad(1, i));
1233 quad_data.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE)(i, 5) =
1234 diffN_MBQUAD2y(pts_quad(0, i));
1235 quad_data.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE)(i, 6) =
1236 diffN_MBQUAD3x(pts_quad(1, i));
1237 quad_data.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE)(i, 7) =
1238 diffN_MBQUAD3y(pts_quad(0, i));
1239 }
1240
1241 if (choice_value == H1QUAD) {
1242 CHKERR QuadPolynomialBase().getValue(
1243 pts_quad, boost::shared_ptr<BaseFunctionCtx>(new EntPolynomialBaseCtx(
1244 quad_data, H1, CONTINUOUS, AINSWORTH_LEGENDRE_BASE, NOBASE)));
1245 ierr = check_same_matrix(
1246 quad_data.dataOnEntities[MBVERTEX][0].getN(NOBASE),
1247 quad_data.dataOnEntities[MBVERTEX][0].getN(AINSWORTH_LEGENDRE_BASE),
1248 eps);
1249 CHKERRG(ierr);
1250 double sum = 0, diff_sum = 0;
1251 for (int ee = 0; ee < 4; ee++) {
1252 sum += sum_matrix(
1253 quad_data.dataOnEntities[MBEDGE][ee].getN(AINSWORTH_LEGENDRE_BASE));
1254 diff_sum += sum_matrix(quad_data.dataOnEntities[MBEDGE][ee].getDiffN(
1255 AINSWORTH_LEGENDRE_BASE));
1256 }
1257 sum += sum_matrix(
1258 quad_data.dataOnEntities[MBQUAD][0].getN(AINSWORTH_LEGENDRE_BASE));
1259 diff_sum += sum_matrix(quad_data.dataOnEntities[MBQUAD][0].getDiffN(
1260 AINSWORTH_LEGENDRE_BASE));
1261
1262 std::cout << sum << " " << diff_sum << endl;
1263
1264 if (std::abs(3.58249 - sum) > eps)
1265 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
1266
1267
1268 if (std::abs(-0.134694 - diff_sum) > eps)
1269 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong result");
1270
1271 }
1272 }
1273 CATCH_ERRORS;
1274
1275 MoFEM::Core::Finalize();
1276}
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
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
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
NBEDGE_L2
#define NBEDGE_L2(P)
Number of base functions on edge from L2 space.
Definition h1_hdiv_hcurl_l2.h:48
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:38
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:2064
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:759
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:1127
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