PeriDEM 0.2.0
PeriDEM -- Peridynamics-based high-fidelity model for granular media
Loading...
Searching...
No Matches
testFeLib.cpp
Go to the documentation of this file.
1/*
2 * -------------------------------------------
3 * Copyright (c) 2021 - 2024 Prashant K. Jha
4 * -------------------------------------------
5 * PeriDEM https://github.com/prashjha/PeriDEM
6 *
7 * Distributed under the Boost Software License, Version 1.0. (See accompanying
8 * file LICENSE)
9 */
10
11#include "testFeLib.h"
12#include "fe/quadElem.h"
13#include "fe/triElem.h"
14#include "fe/tetElem.h"
15#include "util/point.h"
16#include "util/feElementDefs.h"
17#include "util/methods.h"
18#include "nsearch/nsearch.h"
19#include <csv/csv.h>
20#include <algorithm>
21#include <fstream>
22#include <string>
23
24namespace {
25
26 int debug_id = -1;
27 const double tol = 1.0E-12;
28
29 void readNodes(const std::string &filename,
30 std::vector<util::Point> &nodes) {
31
32 // csv reader
33 io::CSVReader<3> in(filename);
34 double x, y, z;
35 while (in.read_row(x, y, z))
36 nodes.emplace_back(x, y, z);
37 }
38
39 size_t readElements(const std::string &filename, const size_t &elem_type,
40 std::vector<size_t> &elements) {
41
42 if (elem_type == util::vtk_type_triangle) {
43 io::CSVReader<3> in(filename);
44 std::vector<size_t> ids(3, 0);
45 while (in.read_row(ids[0], ids[1], ids[2])) {
46 for (auto id: ids)
47 elements.emplace_back(id);
48 }
49
50 size_t num_vertex = util::vtk_map_element_to_num_nodes[elem_type];
51 return elements.size() / num_vertex;
52 } else if (elem_type == util::vtk_type_quad) {
53 io::CSVReader<4> in(filename);
54 std::vector<size_t> ids(4, 0);
55 while (in.read_row(ids[0], ids[1], ids[2], ids[3])) {
56 for (auto id: ids)
57 elements.emplace_back(id);
58 }
59
60 size_t num_vertex = util::vtk_map_element_to_num_nodes[elem_type];
61 return elements.size() / num_vertex;
62 } else if (elem_type == util::vtk_type_tetra) {
63 io::CSVReader<4> in(filename);
64 std::vector<size_t> ids(4, 0);
65 while (in.read_row(ids[0], ids[1], ids[2], ids[3])) {
66 for (auto id: ids)
67 elements.emplace_back(id);
68 }
69
70 size_t num_vertex = util::vtk_map_element_to_num_nodes[elem_type];
71 return elements.size() / num_vertex;
72 } else {
73 std::cerr << "Error: readElements() only supports vtk_type_triangle, vtk_type_quad, and vtk_type_tetra elem_type in testing.\n";
74 exit(1);
75 }
76 }
77
78 bool checkRefIntegration(const size_t &n, const size_t &i,
79 const size_t &j,
80 const std::vector<fe::QuadData> &qds,
81 double &I_exact) {
82
83 double I_approx = 0.;
84 for (auto qd: qds)
85 I_approx += qd.d_w * std::pow(qd.d_p.d_x, i) * std::pow(qd.d_p.d_y, j);
86
87 if (std::abs(I_exact - I_approx) > tol) {
88 std::cout << "Error in order = " << n << ". Exact integration = " << I_exact
89 << " and approximate integration = " << I_approx
90 << " of polynomial of order (i = " << i << " + j = " << j
91 << ") = " << i + j << " over reference element "
92 << "is not matching using quadrature points.\n";
93
94 return false;
95 }
96
97 return true;
98 }
99
100 bool checkRefIntegration(const size_t &n, const size_t &i,
101 const size_t &j, const size_t &k,
102 const std::vector<fe::QuadData> &qds,
103 double &I_exact) {
104
105 double I_approx = 0.;
106 for (auto qd: qds)
107 I_approx += qd.d_w * std::pow(qd.d_p.d_x, i) * std::pow(qd.d_p.d_y, j) *
108 std::pow(qd.d_p.d_z, k);
109
110 if (std::abs(I_exact - I_approx) > tol) {
111 std::cout << "Error in order = " << n << ". Exact integration = " << I_exact
112 << " and approximate integration = " << I_approx
113 << " of polynomial of order (i = " << i << " + j = " << j
114 << " + k = " << k << ") = " << i + j + k
115 << " over reference element "
116 << "is not matching using quadrature points.\n";
117
118 std::cout << "Print " << i << " " << j << " " << k
119 << " debug id = " << debug_id << "\n";
120 for (auto qd: qds)
121 std::cout << qd.printStr() << "\n";
122
123 return false;
124 }
125
126 return true;
127 }
128
129} // namespace
130
131//
132// Interface methods
133//
134
135double test::getNChooseR(size_t n, size_t r) {
136
137 if (r == 0)
138 return 1.;
139
140 double a = 1.;
141 for (size_t i = 1; i <= r; i++)
142 a *= double(n - i + 1) / double(i);
143
144 return a;
145}
146
147double test::getExactIntegrationRefTri(size_t alpha, size_t beta) {
148
149 // compute exact integration of s^\alpha t^\beta
150 double I = 0.;
151 for (size_t k = 0; k <= beta + 1; k++) {
152 if (k % 2 == 0)
153 I +=
154 test::getNChooseR(beta + 1, k) / double((alpha + 1 + k) * (beta + 1));
155 else
156 I -=
157 test::getNChooseR(beta + 1, k) / double((alpha + 1 + k) * (beta + 1));
158 }
159
160 return I;
161}
162
163double test::getExactIntegrationRefQuad(size_t alpha, size_t beta) {
164
165 // compute exact integration of s^\alpha t^\beta
166 if (alpha % 2 == 0 and beta % 2 == 0)
167 return 4. / double((alpha + 1) * (beta + 1));
168 else
169 return 0.;
170}
171
172double test::getExactIntegrationRefTet(size_t alpha, size_t beta,
173 size_t theta) {
174
175 double I = 0.;
176 for (size_t i = 0; i <= theta + 1; i++) {
177
178 double factor_i = test::getNChooseR(theta + 1, i) /
179 (double(theta + 1) * double(i + beta + 1));
180 if (i % 2 != 0)
181 factor_i = factor_i * (-1.);
182
183 for (size_t j = 0; j <= theta + beta + 2 + 1; j++) {
184
185 double factor_j =
186 test::getNChooseR(theta + beta + 2, j) / (double(j + alpha + 1));
187 if (j % 2 != 0)
188 factor_j = factor_j * (-1.);
189
190 I += factor_i * factor_j;
191 }
192 }
193
194 return I;
195}
196
197
198void test::testLineElem(size_t n, std::string filepath) { return; }
199
200void test::testTriElem(size_t n, std::string filepath) {
201
202 //
203 // Test1: We test accuracy of integrals of polynomials over reference
204 // triangle. Reference triangle {(0,0), (1,0), (0,1)}.
205 //
206 // Test2: We consider simple mesh in meshFeTest.txt over square domain
207 // [0,1]^2 and test the accuracy of polynomials over square domain.
208 //
209
210 // get Quadrature
211 auto quad = fe::TriElem(n);
212
213 //
214 // Test 1
215 //
216 size_t error_test_1 = 0;
217 {
218 // T1 (reference triangle)
219 // get quad points at reference triangle
220 std::vector<util::Point> nodes = {util::Point(), util::Point(1., 0., 0.),
221 util::Point(0., 1., 0.)};
222 std::vector<fe::QuadData> qds = quad.getQuadPoints(nodes);
223 double sum = 0.;
224 for (auto qd : qds)
225 sum += qd.d_w;
226
227 if (std::abs(sum - 0.5) > tol) {
228 std::cout << "Error in order = " << n
229 << ". Sum of quad weights is not "
230 "equal to area of reference "
231 "triangle.\n";
232 error_test_1++;
233 }
234
235 //
236 // test the exactness of integration for polynomial
237 //
238 for (size_t i = 0; i <= n; i++)
239 for (size_t j = 0; j <= n; j++) {
240
241 if (i + j > n)
242 continue;
243
244 //
245 // when {(0,0), (1,0), (0,1)}
246 //
247 nodes = {util::Point(), util::Point(1., 0., 0.),
248 util::Point(0., 1., 0.)};
249 qds = quad.getQuadPoints(nodes);
250 // test integration of polynomial f(s,t) = s^i t^j
251 // get the exact integration
252 double I_exact = test::getExactIntegrationRefTri(i, j);
253 if (!checkRefIntegration(n, i, j, qds, I_exact))
254 error_test_1++;
255
256 //
257 // when vertices are {(1,0), (0,1), (0,0)}
258 //
259 nodes = {util::Point(1., 0., 0.), util::Point(0., 1., 0.),
260 util::Point()};
261 qds = quad.getQuadPoints(nodes);
262 //
263 // After changing the order of vertices, we have got a new
264 // triangle which is in coordinate system (x,y) and we are
265 // integrating function f(x,y) = x^i y^j
266 //
267 // The quad data we have got is such that quad point is in (x,y)
268 // coordinate, weight is such that determinant of the Jacobian is
269 // included in the weight.
270 //
271 // Thus the following method for I_approx is correct.
272 if (!checkRefIntegration(n, i, j, qds, I_exact))
273 error_test_1++;
274
275 //
276 // when vertices are {(0,1), (0,0), (1,0)}
277 //
278 nodes = {util::Point(0., 1., 0.), util::Point(),
279 util::Point(1., 0., 0.)};
280 qds = quad.getQuadPoints(nodes);
281 if (!checkRefIntegration(n, i, j, qds, I_exact))
282 error_test_1++;
283 }
284 } // Test 1
285
286 //
287 // Test 2
288 //
289 size_t error_test_2 = 0;
290 {
291 static std::vector<util::Point> nodes;
292 static std::vector<size_t> elements;
293 static size_t num_vertex = 3;
294 static size_t elem_type = util::vtk_type_triangle;
295 static size_t num_elems = 0;
296 if (num_elems == 0) {
297 readNodes(filepath + "/triMesh_nodes.csv", nodes);
298 num_elems = readElements(filepath + "/triMesh_elements.csv", elem_type,
299 elements);
300 }
301
302 // loop over polynomials
303 for (size_t i = 0; i <= n; i++)
304 for (size_t j = 0; j <= n; j++) {
305
306 if (i + j > n)
307 continue;
308
309 double I_exact = 1. / (double(i + 1) * double(j + 1));
310 double I_approx = 0.;
311 // loop over elements and compute I_approx
312 for (size_t e = 0; e < num_elems; e++) {
313 std::vector<util::Point> enodes = {
314 nodes[elements[num_vertex * e + 0]],
315 nodes[elements[num_vertex * e + 1]],
316 nodes[elements[num_vertex * e + 2]]};
317 std::vector<fe::QuadData> qds = quad.getQuadPoints(enodes);
318 for (auto qd : qds)
319 I_approx +=
320 qd.d_w * std::pow(qd.d_p.d_x, i) * std::pow(qd.d_p.d_y, j);
321 }
322
323 if (std::abs(I_exact - I_approx) > tol) {
324 std::cout << "Error in order = " << n
325 << ". Exact integration = " << I_exact
326 << " and approximate integration = " << I_approx
327 << " of polynomial of order (i = " << i << " + j = " << j
328 << ") = " << i + j << " over square domain [0,1]x[0,1] "
329 << "is not matching using quadrature points.\n";
330
331 error_test_2++;
332 }
333 }
334 }
335
336 if (n == 1) {
337 std::cout << "**********************************\n";
338 std::cout << "Triangle Quadrature Test\n";
339 std::cout << "**********************************\n";
340 }
341 std::cout << "Quad order = " << n << ". ";
342 if (error_test_1 == 0)
343 std::cout << "TEST 1 : PASS. ";
344 else
345 std::cout << "TEST 1 : FAIL. ";
346 if (error_test_2 == 0)
347 std::cout << "TEST 2 : PASS. ";
348 else
349 std::cout << "TEST 2 : FAIL. ";
350 std::cout << "\n";
351}
352
353void test::testQuadElem(size_t n, std::string filepath) {
354
355 //
356 // Test1: We test accuracy of integrals of polynomials over reference
357 // quadrangle. Reference triangle {(-1,-1), (1,-1), (1,1), (-1,1)}.
358 //
359 // Test2: We consider simple mesh in meshFeTest.txt over square domain
360 // [0,1]^2 and test the accuracy of polynomials over square domain.
361 //
362
363 // get Quadrature
364 auto quad = fe::QuadElem(n);
365
366 //
367 // Test 1
368 //
369 size_t error_test_1 = 0;
370 {
371 // T1 (reference quadrangle)
372 // get quad points at reference triangle
373 std::vector<util::Point> nodes = {
374 util::Point(-1., -1., 0.), util::Point(1., -1., 0.),
375 util::Point(1., 1., 0.), util::Point(-1., 1., 0.)};
376 std::vector<fe::QuadData> qds = quad.getQuadPoints(nodes);
377 double sum = 0.;
378 for (auto qd : qds)
379 sum += qd.d_w;
380
381 if (std::abs(sum - 4.0) > tol) {
382 std::cout << "Error in order = " << n
383 << ". Sum of quad weights is not "
384 "equal to area of reference "
385 "quadrangle.\n";
386 error_test_1++;
387 }
388
389 //
390 // test the exactness of integration for polynomial
391 //
392 for (size_t i = 0; i <= 2 * n - 1; i++)
393 for (size_t j = 0; j <= 2 * n - 1; j++) {
394
395 //
396 // when {(-1,-1), (1,-1), (1,1), (-1,1)}
397 //
398 nodes = {util::Point(-1., -1., 0.), util::Point(1., -1., 0.),
399 util::Point(1., 1., 0.), util::Point(-1., 1., 0.)};
400 qds = quad.getQuadPoints(nodes);
401 // test integration of polynomial f(s,t) = s^i t^j
402 // get the exact integration
403 double I_exact = test::getExactIntegrationRefQuad(i, j);
404 if (!checkRefIntegration(n, i, j, qds, I_exact))
405 error_test_1++;
406
407 //
408 // when {(-1,1), (-1,-1), (1,-1), (1,1)}
409 //
410 nodes = {util::Point(-1., 1., 0.), util::Point(-1., -1., 0.),
411 util::Point(1., -1., 0.), util::Point(1., 1., 0.)};
412 qds = quad.getQuadPoints(nodes);
413 //
414 // After changing the order of vertices, we have got a new
415 // triangle which is in coordinate system (x,y) and we are
416 // integrating function f(x,y) = x^i y^j
417 //
418 // The quad data we have got is such that quad point is in (x,y)
419 // coordinate, weight is such that determinant of the Jacobian is
420 // included in the weight.
421 //
422 // Thus the following method for I_approx is correct.
423 if (!checkRefIntegration(n, i, j, qds, I_exact))
424 error_test_1++;
425
426 //
427 // when {(1,1), (-1,1), (-1,-1), (1,-1)}
428 //
429 nodes = {util::Point(1., 1., 0.), util::Point(-1., 1., 0.),
430 util::Point(-1., -1., 0.), util::Point(1., -1., 0.)};
431 qds = quad.getQuadPoints(nodes);
432 if (!checkRefIntegration(n, i, j, qds, I_exact))
433 error_test_1++;
434
435 //
436 // when {(1,-1), (1,1), (-1,1), (-1,-1)}
437 //
438 nodes = {util::Point(1., -1., 0.), util::Point(1., 1., 0.),
439 util::Point(-1., 1., 0.), util::Point(-1., -1., 0.)};
440 qds = quad.getQuadPoints(nodes);
441 if (!checkRefIntegration(n, i, j, qds, I_exact))
442 error_test_1++;
443 }
444 } // Test 1
445
446 //
447 // Test 2
448 //
449 size_t error_test_2 = 0;
450 {
451 static std::vector<util::Point> nodes;
452 static std::vector<size_t> elements;
453 static size_t num_vertex = 4;
454 static size_t elem_type = util::vtk_type_quad;
455 static size_t num_elems = 0;
456 if (num_elems == 0) {
457 readNodes(filepath + "/quadMesh_nodes.csv", nodes);
458 num_elems = readElements(filepath + "/quadMesh_elements.csv", elem_type,
459 elements);
460 }
461
462 // loop over polynomials
463 for (size_t i = 0; i <= 2 * n - 1; i++)
464 for (size_t j = 0; j <= 2 * n - 1; j++) {
465
466 double I_exact = 1. / (double(i + 1) * double(j + 1));
467 double I_approx = 0.;
468 // loop over elements and compute I_approx
469 for (size_t e = 0; e < num_elems; e++) {
470 std::vector<util::Point> enodes = {
471 nodes[elements[num_vertex * e + 0]],
472 nodes[elements[num_vertex * e + 1]],
473 nodes[elements[num_vertex * e + 2]],
474 nodes[elements[num_vertex * e + 3]]};
475 std::vector<fe::QuadData> qds = quad.getQuadPoints(enodes);
476 for (auto qd : qds)
477 I_approx +=
478 qd.d_w * std::pow(qd.d_p.d_x, i) * std::pow(qd.d_p.d_y, j);
479 }
480
481 if (std::abs(I_exact - I_approx) > tol) {
482 std::cout << "Error in order = " << n
483 << ". Exact integration = " << I_exact
484 << " and approximate integration = " << I_approx
485 << " of polynomial of order (i = " << i << " + j = " << j
486 << ") = " << i + j << " over square domain [0,1]x[0,1] "
487 << "is not matching using quadrature points.\n";
488
489 error_test_2++;
490 }
491 }
492 }
493
494 if (n == 1) {
495 std::cout << "**********************************\n";
496 std::cout << "Quadrangle Quadrature Test\n";
497 std::cout << "**********************************\n";
498 }
499 std::cout << "Quad order = " << n << ". ";
500 std::cout << (error_test_1 == 0 ? "TEST 1 : PASS. " : "TEST 1 : FAIL. ");
501 std::cout << (error_test_2 == 0 ? "TEST 2 : PASS. \n" : "TEST 2 : FAIL. \n");
502}
503
504void test::testTriElemTime(size_t n, size_t N) {
505
506 // get Quadrature
507 auto quad = fe::TriElem(n);
508
509 //
510 // Test 1
511 //
512 std::vector<util::Point> nodes = {util::Point(2., 2., 0.),
513 util::Point(4., 2., 0.),
514 util::Point(2., 4., 0.)};
515 size_t num_vertex = 3;
516 // std::vector<size_t> elements;
517 // for (size_t i = 0; i < 3 * N; i++)
518 // elements.emplace_back(i % 2);
519 std::vector<std::vector<size_t>> elements;
520 for (size_t i = 0; i < N; i++)
521 elements.emplace_back(std::vector<size_t>{0, 1, 2});
522
523 auto t11 = steady_clock::now();
524 // method 1: Compute quad points on the fly
525 // loop over elements and compute I_approx
526 double sum = 0.;
527 for (size_t e = 0; e < N; e++) {
528 // std::vector<util::Point> enodes = {nodes[elements[num_vertex * e +
529 // 0]],
530 // nodes[elements[num_vertex * e +
531 // 1]], nodes[elements[num_vertex * e
532 // + 2]]};
533 std::vector<util::Point> enodes = {
534 nodes[elements[e][0]], nodes[elements[e][1]], nodes[elements[e][2]]};
535 std::vector<fe::QuadData> qds = quad.getQuadPoints(enodes);
536 for (auto qd : qds)
537 sum += qd.d_w * (qd.d_shapes[0] + qd.d_shapes[1] + qd.d_shapes[2]);
538 }
539 auto t12 = steady_clock::now();
540
541 // method 2: Compute quad points in the beginning and use it when needed
542 size_t num_quad_pts = 0;
543 std::vector<fe::QuadData> quad_data;
544 for (size_t e = 0; e < N; e++) {
545 std::vector<fe::QuadData> qds = quad.getQuadPoints(nodes);
546 if (e == 0)
547 num_quad_pts = qds.size();
548 for (auto qd : qds)
549 quad_data.emplace_back(qd);
550 }
551
552 auto t21 = steady_clock::now();
553 sum = 0.;
554 for (size_t e = 0; e < N; e++) {
555 for (size_t q = 0; q < num_quad_pts; q++) {
556 fe::QuadData qd = quad_data[e * num_quad_pts + q];
557 sum += qd.d_w * (qd.d_shapes[0] + qd.d_shapes[1] + qd.d_shapes[2]);
558 }
559 }
560 auto t22 = steady_clock::now();
561
562 if (n == 1 and N == 1000) {
563 std::cout << "**********************************\n";
564 std::cout << "Quadrature Time Efficiency Test\n";
565 std::cout << "**********************************\n";
566 }
567 std::cout << "Quad order = " << n << ". Num Elements = " << N << ".\n ";
568 double dt_1 = util::methods::timeDiff(t12, t11, "seconds");
569 double dt_2 = util::methods::timeDiff(t21, t22, "seconds");
570 double perc = (dt_1 - dt_2) * 100. / dt_2;
571 double qpt_mem = 13 * sizeof(double);
572 double mem2 = double(quad_data.capacity() * qpt_mem) / double(1000000);
573 std::cout << " dt1 = " << dt_1 << ", dt2 = " << dt_2 << ", perc = " << perc
574 << ". Mem saved = " << mem2 << " MB.\n";
575}
576
577void test::testTetElem(size_t n, std::string filepath) {
578
579 //
580 // Test1: We test accuracy of integrals of polynomials over reference
581 // tetrahedron. Reference element {(0,0,0), (1,0,0), (0,1,0), (0,0,1)}.
582 //
583 // Test2: We consider simple mesh in meshFeTest.txt over cubic domain
584 // [0,1]^3 and test the accuracy of polynomials over cubic domain.
585 //
586
587 // get Quadrature
588 auto quad = fe::TetElem(n);
589
590 //
591 // Test 1
592 //
593 size_t error_test_1 = 0;
594 {
595 // T1 (reference triangle)
596 // get quad points at reference triangle
597 std::vector<util::Point> nodes = {util::Point(), util::Point(1., 0., 0.),
598 util::Point(0., 1., 0.),
599 util::Point(0., 0., 1.)};
600 std::vector<fe::QuadData> qds = quad.getQuadPoints(nodes);
601
602 double sum = 0.;
603 for (auto qd : qds)
604 sum += qd.d_w;
605
606 if (std::abs(sum - 1. / 6.) > tol) {
607 std::cout << "Error in order = " << n
608 << ". Sum of quad weights is not "
609 "equal to volume of reference "
610 "tetrahedron.\n";
611 error_test_1++;
612 }
613
614 //
615 // test the exactness of integration for polynomial
616 //
617 for (size_t i = 0; i <= n; i++)
618 for (size_t j = 0; j <= n; j++)
619 for (size_t k = 0; k <= n; k++) {
620
621 if (i + j + k > n)
622 continue;
623
624 //
625 // +ve order of indices are:
626 // {0,1,2,3}; {1,2,0,3}; {2,3,0,1}; {0,3,1,2}
627
628 //
629 // when {(0,0,0), (1,0,0), (0,1,0), (0,0,1)}
630 //
631 nodes = {util::Point(), util::Point(1., 0., 0.),
632 util::Point(0., 1., 0.), util::Point(0., 0., 1.)};
633 qds = quad.getQuadPoints(nodes);
634 // test integration of polynomial f(s,t) = s^i t^j
635 // get the exact integration
636 debug_id = 0;
637 double I_exact = test::getExactIntegrationRefTet(i, j, k);
638 if (!checkRefIntegration(n, i, j, k, qds, I_exact)) {
639 error_test_1++;
640 }
641
642 //
643 // when vertices are {(1,0,0), (0,1,0), (0,0,1), (0,0,0)}
644 //
645 nodes = {util::Point(1., 0., 0.), util::Point(0., 1., 0.),
646 util::Point(0., 0., 0.), util::Point(0., 0., 1.)};
647 qds = quad.getQuadPoints(nodes);
648 //
649 // After changing the order of vertices, we have got a new
650 // triangle which is in coordinate system (x,y) and we are
651 // integrating function f(x,y) = x^i y^j
652 //
653 // The quad data we have got is such that quad point is in (x,y)
654 // coordinate, weight is such that determinant of the Jacobian is
655 // included in the weight.
656 //
657 // Thus the following method for I_approx is correct.
658 debug_id = 1;
659 if (!checkRefIntegration(n, i, j, k, qds, I_exact))
660 error_test_1++;
661
662 //
663 // when vertices are {(0,1,0), (0,0,1), (0,0,0), (1,0,0)}
664 //
665 nodes = {util::Point(0., 1., 0.), util::Point(0., 0., 1.),
666 util::Point(0., 0., 0.), util::Point(1., 0., 0.)};
667 qds = quad.getQuadPoints(nodes);
668 debug_id = 2;
669 if (!checkRefIntegration(n, i, j, k, qds, I_exact))
670 error_test_1++;
671
672 //
673 // when vertices are {(0,0,1), (0,0,0), (1,0,0), (0,1,0)}
674 //
675 nodes = {util::Point(0., 0., 0.), util::Point(0., 0., 1.),
676 util::Point(1., 0., 0.), util::Point(0., 1., 0.)};
677 qds = quad.getQuadPoints(nodes);
678 debug_id = 3;
679 if (!checkRefIntegration(n, i, j, k, qds, I_exact))
680 error_test_1++;
681 }
682 } // Test 1
683
684 //
685 // Test 2
686 //
687 size_t error_test_2 = 0;
688 if (false) {
689 static std::vector<util::Point> nodes;
690 static std::vector<size_t> elements;
691 static size_t num_vertex = 4;
692 static size_t elem_type = util::vtk_type_tetra;
693 static size_t num_elems = 0;
694 if (num_elems == 0) {
695 readNodes(filepath + "tetMesh_nodes.csv", nodes);
696 num_elems = readElements(filepath + "tetMesh_elements.csv", elem_type,
697 elements);
698 }
699
700 // loop over polynomials
701 for (size_t i = 0; i <= n; i++)
702 for (size_t j = 0; j <= n; j++)
703 for (size_t k = 0; k <= n; k++) {
704
705 if (i + j + k > n)
706 continue;
707
708 double I_exact = 1. / (double(i + 1) * double(j + 1) * double(k + 1));
709 double I_approx = 0.;
710 // loop over elements and compute I_approx
711 for (size_t e = 0; e < num_elems; e++) {
712 std::vector<util::Point> enodes = {
713 nodes[elements[num_vertex * e + 0]],
714 nodes[elements[num_vertex * e + 1]],
715 nodes[elements[num_vertex * e + 2]],
716 nodes[elements[num_vertex * e + 3]]};
717 std::vector<fe::QuadData> qds = quad.getQuadPoints(enodes);
718 for (auto qd : qds) {
719 I_approx += qd.d_w * std::pow(qd.d_p.d_x, i) *
720 std::pow(qd.d_p.d_y, j) * std::pow(qd.d_p.d_z, k);
721
722 if (false) {
723
724 std::cout << "Print " << i << " " << j << " " << k << "\n";
725 std::cout << util::io::printStr(enodes) << "\n";
726 std::vector<size_t> enode_ids = {
727 elements[num_vertex * e + 0],
728 elements[num_vertex * e + 1],
729 elements[num_vertex * e + 2],
730 elements[num_vertex * e + 3]};
731 std::cout << util::io::printStr(enode_ids) << "\n";
732 std::cout << qd.printStr() << "\n";
733 }
734
735 }
736 }
737
738 if (std::abs(I_exact - I_approx) > tol) {
739 std::cout << "Error in order = " << n
740 << ". Exact integration = " << I_exact
741 << " and approximate integration = " << I_approx
742 << " of polynomial of order (i = " << i << " + j = " << j
743 << " + k = " << k << ") = " << i + j + k
744 << " over cubic domain [0,1]x[0,1]x[0,1] "
745 << "is not matching using quadrature points.\n";
746
747 error_test_2++;
748 }
749 }
750 }
751
752 if (n == 1) {
753 std::cout << "**********************************\n";
754 std::cout << "Tetrahedron Quadrature Test\n";
755 std::cout << "**********************************\n";
756 }
757 std::cout << "Quad order = " << n << ". ";
758 if (error_test_1 == 0)
759 std::cout << "TEST 1 : PASS. ";
760 else
761 std::cout << "TEST 1 : FAIL. ";
762 // if (error_test_2 == 0)
763 // std::cout << "TEST 2 : PASS. ";
764 // else
765 // std::cout << "TEST 2 : FAIL. ";
766 std::cout << "\n";
767}
fe::QuadElem
A class for mapping and quadrature related operations for bi-linear quadrangle element.
Definition quadElem.h:64
fe::TetElem
A class for mapping and quadrature related operations for linear tetrahedron element.
Definition tetElem.h:141
fe::TriElem
A class for mapping and quadrature related operations for linear triangle element.
Definition triElem.h:91
util::vtk_map_element_to_num_nodes
static int vtk_map_element_to_num_nodes[16]
Map from element type to number of nodes (for vtk)
Definition feElementDefs.h:85
util::vtk_type_triangle
static const int vtk_type_triangle
Integer flag for triangle element.
Definition feElementDefs.h:55
util::vtk_type_quad
static const int vtk_type_quad
Integer flag for quad element.
Definition feElementDefs.h:67
util::vtk_type_tetra
static const int vtk_type_tetra
Integer flag for tetrahedron element.
Definition feElementDefs.h:70
anonymous_namespace{testFeLib.cpp}::tol
const double tol
Definition testFeLib.cpp:27
anonymous_namespace{testFeLib.cpp}::readElements
size_t readElements(const std::string &filename, const size_t &elem_type, std::vector< size_t > &elements)
Definition testFeLib.cpp:39
anonymous_namespace{testFeLib.cpp}::readNodes
void readNodes(const std::string &filename, std::vector< util::Point > &nodes)
Definition testFeLib.cpp:29
anonymous_namespace{testFeLib.cpp}::checkRefIntegration
bool checkRefIntegration(const size_t &n, const size_t &i, const size_t &j, const std::vector< fe::QuadData > &qds, double &I_exact)
Definition testFeLib.cpp:78
anonymous_namespace{testFeLib.cpp}::debug_id
int debug_id
Definition testFeLib.cpp:26
test::testTriElemTime
void testTriElemTime(size_t n, size_t N)
Computes the time needed when quad data for elements are stored and when they are computed as and whe...
Definition testFeLib.cpp:504
test::getExactIntegrationRefTri
double getExactIntegrationRefTri(size_t alpha, size_t beta)
Computes integration of polynomial exactly over reference triangle.
Definition testFeLib.cpp:147
test::getNChooseR
double getNChooseR(size_t n, size_t r)
Computes "n choose r".
Definition testFeLib.cpp:135
test::getExactIntegrationRefQuad
double getExactIntegrationRefQuad(size_t alpha, size_t beta)
Computes integration of polynomial exactly over reference quadrangle.
Definition testFeLib.cpp:163
test::testLineElem
void testLineElem(size_t n, std::string filepath)
Perform test on quadrature points on line elements (NOT IMPLEMENTED)
Definition testFeLib.cpp:198
test::testTetElem
void testTetElem(size_t n, std::string filepath)
Perform test on quadrature points on tetrahedral elements.
Definition testFeLib.cpp:577
test::testQuadElem
void testQuadElem(size_t n, std::string filepath)
Perform test on quadrature points on quadrangle elements.
Definition testFeLib.cpp:353
test::getExactIntegrationRefTet
double getExactIntegrationRefTet(size_t alpha, size_t beta, size_t theta)
Computes integration of polynomial exactly over reference tetrahedral.
Definition testFeLib.cpp:172
test::testTriElem
void testTriElem(size_t n, std::string filepath)
Perform test on quadrature points on triangle elements.
Definition testFeLib.cpp:200
util::io::printStr
std::string printStr(const T &msg, int nt=print_default_tab)
Returns formatted string for output.
Definition io.h:54
util::methods::timeDiff
float timeDiff(std::chrono::steady_clock::time_point begin, std::chrono::steady_clock::time_point end, std::string unit="microseconds")
Returns difference between two times.
Definition methods.h:304
fe::QuadData
A struct to store the quadrature data. List of data are.
Definition quadData.h:23
fe::QuadData::d_shapes
std::vector< double > d_shapes
Value of shape functions at quad point p.
Definition quadData.h:37
fe::QuadData::d_w
double d_w
Quadrature weight.
Definition quadData.h:26
util::Point
A structure to represent 3d vectors.
Definition point.h:30