v0.13.2
Loading...
Searching...
No Matches
Functions | Variables
testing_base_functions.cpp File Reference
#include <MoFEM.hpp>
#include <quad.h>

Go to the source code of this file.

Functions

static double sum_matrix (MatrixDouble &m)
 
int main (int argc, char *argv[])
 

Variables

static char help [] = "testing interface inserting algorithm\n\n"
 

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 20 of file testing_base_functions.cpp.

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

◆ sum_matrix()

static double sum_matrix ( MatrixDouble m)
static

Definition at line 10 of file testing_base_functions.cpp.

10 {
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}
FTensor::Index< 'm', SPACE_DIM > m

Variable Documentation

◆ help

char help[] = "testing interface inserting algorithm\n\n"
static

Definition at line 8 of file testing_base_functions.cpp.