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

◆ 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}
m
FTensor::Index< 'm', 3 > m
Definition shallow_wave.cpp:80

Variable Documentation

◆ help

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

Definition at line 8 of file testing_base_functions.cpp.

Generated by Doxygen 1.12.0 and hosted at University of Glasgow