PeriDEM 0.3.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 - 2026 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/elemIncludes.h"
13#include "fe/bMatrix.h"
14#include "util/point.h"
15#include "util/matrix.h"
16#include "util/feElementDefs.h"
17#include "util/vecMethods.h"
18#include "nsearch/nsearch.h"
19#include <csv/csv.h>
20#include <algorithm>
21#include <cmath>
22#include <cstdlib>
23#include <fstream>
24#include <string>
25
26namespace {
27
28 int debug_id = -1;
29 const double tol = 1.0E-12;
30
31 util::SymMatrix3 strainFromB(const fe::B &B,
32 const std::vector<util::Point> &u) {
33 std::vector<double> uflat(B.nDof());
34 const int dim = B.dim();
35 for (size_t a = 0; a < u.size(); a++) {
36 uflat[dim * static_cast<int>(a)] = u[a][0];
37 if (dim > 1)
38 uflat[dim * static_cast<int>(a) + 1] = u[a][1];
39 if (dim > 2)
40 uflat[dim * static_cast<int>(a) + 2] = u[a][2];
41 }
42 std::vector<double> e(B.nStrain(), 0.);
43 for (int i = 0; i < B.nStrain(); i++)
44 for (int j = 0; j < B.nDof(); j++)
45 e[i] += B(i, j) * uflat[j];
46 util::SymMatrix3 s;
47 s(0, 0) = e[0];
48 if (dim > 1) {
49 s(1, 1) = e[1];
50 s(0, 1) = (dim == 2) ? e[2] : e[5];
51 }
52 if (dim > 2) {
53 s(2, 2) = e[2];
54 s(1, 2) = e[3];
55 s(0, 2) = e[4];
56 }
57 return s;
58 }
59
60 void readNodes(const std::string &filename,
61 std::vector<util::Point> &nodes) {
62
63 // csv reader
64 io::CSVReader<3> in(filename);
65 double x, y, z;
66 while (in.read_row(x, y, z))
67 nodes.emplace_back(x, y, z);
68 }
69
70 size_t readElements(const std::string &filename, const size_t &elem_type,
71 std::vector<size_t> &elements) {
72
73 if (elem_type == util::vtk_type_triangle) {
74 io::CSVReader<3> in(filename);
75 std::vector<size_t> ids(3, 0);
76 while (in.read_row(ids[0], ids[1], ids[2])) {
77 for (auto id: ids)
78 elements.emplace_back(id);
79 }
80
81 size_t num_vertex = util::vtk_map_element_to_num_nodes[elem_type];
82 return elements.size() / num_vertex;
83 } else if (elem_type == util::vtk_type_quad) {
84 io::CSVReader<4> in(filename);
85 std::vector<size_t> ids(4, 0);
86 while (in.read_row(ids[0], ids[1], ids[2], ids[3])) {
87 for (auto id: ids)
88 elements.emplace_back(id);
89 }
90
91 size_t num_vertex = util::vtk_map_element_to_num_nodes[elem_type];
92 return elements.size() / num_vertex;
93 } else if (elem_type == util::vtk_type_tetra) {
94 io::CSVReader<4> in(filename);
95 std::vector<size_t> ids(4, 0);
96 while (in.read_row(ids[0], ids[1], ids[2], ids[3])) {
97 for (auto id: ids)
98 elements.emplace_back(id);
99 }
100
101 size_t num_vertex = util::vtk_map_element_to_num_nodes[elem_type];
102 return elements.size() / num_vertex;
103 } else {
104 std::cerr << "Error: readElements() only supports vtk_type_triangle, vtk_type_quad, and vtk_type_tetra elem_type in testing.\n";
105 exit(1);
106 }
107 }
108
109 bool checkRefIntegration(const size_t &n, const size_t &i,
110 const size_t &j,
111 const std::vector<fe::QuadData> &qds,
112 double &I_exact) {
113
114 double I_approx = 0.;
115 for (auto qd: qds)
116 I_approx += qd.d_w * std::pow(qd.d_p.d_x, i) * std::pow(qd.d_p.d_y, j);
117
118 if (std::abs(I_exact - I_approx) > tol) {
119 std::cout << "Error in order = " << n << ". Exact integration = " << I_exact
120 << " and approximate integration = " << I_approx
121 << " of polynomial of order (i = " << i << " + j = " << j
122 << ") = " << i + j << " over reference element "
123 << "is not matching using quadrature points.\n";
124
125 return false;
126 }
127
128 return true;
129 }
130
131 bool checkRefIntegration(const size_t &n, const size_t &i,
132 const size_t &j, const size_t &k,
133 const std::vector<fe::QuadData> &qds,
134 double &I_exact) {
135
136 double I_approx = 0.;
137 for (auto qd: qds)
138 I_approx += qd.d_w * std::pow(qd.d_p.d_x, i) * std::pow(qd.d_p.d_y, j) *
139 std::pow(qd.d_p.d_z, k);
140
141 if (std::abs(I_exact - I_approx) > tol) {
142 std::cout << "Error in order = " << n << ". Exact integration = " << I_exact
143 << " and approximate integration = " << I_approx
144 << " of polynomial of order (i = " << i << " + j = " << j
145 << " + k = " << k << ") = " << i + j + k
146 << " over reference element "
147 << "is not matching using quadrature points.\n";
148
149 std::cout << "Print " << i << " " << j << " " << k
150 << " debug id = " << debug_id << "\n";
151 for (auto qd: qds)
152 std::cout << qd.printStr() << "\n";
153
154 return false;
155 }
156
157 return true;
158 }
159
160} // namespace
161
162//
163// Interface methods
164//
165
166double test::getNChooseR(size_t n, size_t r) {
167
168 if (r == 0)
169 return 1.;
170
171 double a = 1.;
172 for (size_t i = 1; i <= r; i++)
173 a *= double(n - i + 1) / double(i);
174
175 return a;
176}
177
178double test::getExactIntegrationRefTri(size_t alpha, size_t beta) {
179
180 // compute exact integration of s^\alpha t^\beta
181 double I = 0.;
182 for (size_t k = 0; k <= beta + 1; k++) {
183 if (k % 2 == 0)
184 I +=
185 test::getNChooseR(beta + 1, k) / double((alpha + 1 + k) * (beta + 1));
186 else
187 I -=
188 test::getNChooseR(beta + 1, k) / double((alpha + 1 + k) * (beta + 1));
189 }
190
191 return I;
192}
193
194double test::getExactIntegrationRefQuad(size_t alpha, size_t beta) {
195
196 // compute exact integration of s^\alpha t^\beta
197 if (alpha % 2 == 0 and beta % 2 == 0)
198 return 4. / double((alpha + 1) * (beta + 1));
199 else
200 return 0.;
201}
202
203double test::getExactIntegrationRefTet(size_t alpha, size_t beta,
204 size_t theta) {
205
206 double I = 0.;
207 for (size_t i = 0; i <= theta + 1; i++) {
208
209 double factor_i = test::getNChooseR(theta + 1, i) /
210 (double(theta + 1) * double(i + beta + 1));
211 if (i % 2 != 0)
212 factor_i = factor_i * (-1.);
213
214 for (size_t j = 0; j <= theta + beta + 2 + 1; j++) {
215
216 double factor_j =
217 test::getNChooseR(theta + beta + 2, j) / (double(j + alpha + 1));
218 if (j % 2 != 0)
219 factor_j = factor_j * (-1.);
220
221 I += factor_i * factor_j;
222 }
223 }
224
225 return I;
226}
227
228
229void test::testLineElem(size_t n, std::string filepath) { return; }
230
231void test::testTriElem(size_t n, std::string filepath) {
232
233 //
234 // Test1: We test accuracy of integrals of polynomials over reference
235 // triangle. Reference triangle {(0,0), (1,0), (0,1)}.
236 //
237 // Test2: We consider simple mesh in meshFeTest.txt over square domain
238 // [0,1]^2 and test the accuracy of polynomials over square domain.
239 //
240
241 // get Quadrature
242 auto quad = fe::TriElem(n);
243
244 //
245 // Test 1
246 //
247 size_t error_test_1 = 0;
248 {
249 // T1 (reference triangle)
250 // get quad points at reference triangle
251 std::vector<util::Point> nodes = {util::Point(), util::Point(1., 0., 0.),
252 util::Point(0., 1., 0.)};
253 std::vector<fe::QuadData> qds = quad.getQuadPoints(nodes);
254 double sum = 0.;
255 for (auto qd : qds)
256 sum += qd.d_w;
257
258 if (std::abs(sum - 0.5) > tol) {
259 std::cout << "Error in order = " << n
260 << ". Sum of quad weights is not "
261 "equal to area of reference "
262 "triangle.\n";
263 error_test_1++;
264 }
265
266 //
267 // test the exactness of integration for polynomial
268 //
269 for (size_t i = 0; i <= n; i++)
270 for (size_t j = 0; j <= n; j++) {
271
272 if (i + j > n)
273 continue;
274
275 //
276 // when {(0,0), (1,0), (0,1)}
277 //
278 nodes = {util::Point(), util::Point(1., 0., 0.),
279 util::Point(0., 1., 0.)};
280 qds = quad.getQuadPoints(nodes);
281 // test integration of polynomial f(s,t) = s^i t^j
282 // get the exact integration
283 double I_exact = test::getExactIntegrationRefTri(i, j);
284 if (!checkRefIntegration(n, i, j, qds, I_exact))
285 error_test_1++;
286
287 //
288 // when vertices are {(1,0), (0,1), (0,0)}
289 //
290 nodes = {util::Point(1., 0., 0.), util::Point(0., 1., 0.),
291 util::Point()};
292 qds = quad.getQuadPoints(nodes);
293 //
294 // After changing the order of vertices, we have got a new
295 // triangle which is in coordinate system (x,y) and we are
296 // integrating function f(x,y) = x^i y^j
297 //
298 // The quad data we have got is such that quad point is in (x,y)
299 // coordinate, weight is such that determinant of the Jacobian is
300 // included in the weight.
301 //
302 // Thus the following method for I_approx is correct.
303 if (!checkRefIntegration(n, i, j, qds, I_exact))
304 error_test_1++;
305
306 //
307 // when vertices are {(0,1), (0,0), (1,0)}
308 //
309 nodes = {util::Point(0., 1., 0.), util::Point(),
310 util::Point(1., 0., 0.)};
311 qds = quad.getQuadPoints(nodes);
312 if (!checkRefIntegration(n, i, j, qds, I_exact))
313 error_test_1++;
314 }
315 } // Test 1
316
317 //
318 // Test 2
319 //
320 size_t error_test_2 = 0;
321 {
322 static std::vector<util::Point> nodes;
323 static std::vector<size_t> elements;
324 static size_t num_vertex = 3;
325 static size_t elem_type = util::vtk_type_triangle;
326 static size_t num_elems = 0;
327 if (num_elems == 0) {
328 readNodes(filepath + "/triMesh_nodes.csv", nodes);
329 num_elems = readElements(filepath + "/triMesh_elements.csv", elem_type,
330 elements);
331 }
332
333 // loop over polynomials
334 for (size_t i = 0; i <= n; i++)
335 for (size_t j = 0; j <= n; j++) {
336
337 if (i + j > n)
338 continue;
339
340 double I_exact = 1. / (double(i + 1) * double(j + 1));
341 double I_approx = 0.;
342 // loop over elements and compute I_approx
343 for (size_t e = 0; e < num_elems; e++) {
344 std::vector<util::Point> enodes = {
345 nodes[elements[num_vertex * e + 0]],
346 nodes[elements[num_vertex * e + 1]],
347 nodes[elements[num_vertex * e + 2]]};
348 std::vector<fe::QuadData> qds = quad.getQuadPoints(enodes);
349 for (auto qd : qds)
350 I_approx +=
351 qd.d_w * std::pow(qd.d_p.d_x, i) * std::pow(qd.d_p.d_y, j);
352 }
353
354 if (std::abs(I_exact - I_approx) > tol) {
355 std::cout << "Error in order = " << n
356 << ". Exact integration = " << I_exact
357 << " and approximate integration = " << I_approx
358 << " of polynomial of order (i = " << i << " + j = " << j
359 << ") = " << i + j << " over square domain [0,1]x[0,1] "
360 << "is not matching using quadrature points.\n";
361
362 error_test_2++;
363 }
364 }
365 }
366
367 if (n == 1) {
368 std::cout << "**********************************\n";
369 std::cout << "Triangle Quadrature Test\n";
370 std::cout << "**********************************\n";
371 }
372 std::cout << "Quad order = " << n << ". ";
373 if (error_test_1 == 0)
374 std::cout << "TEST 1 : PASS. ";
375 else
376 std::cout << "TEST 1 : FAIL. ";
377 if (error_test_2 == 0)
378 std::cout << "TEST 2 : PASS. ";
379 else
380 std::cout << "TEST 2 : FAIL. ";
381 std::cout << "\n";
382}
383
384void test::testQuadElem(size_t n, std::string filepath) {
385
386 //
387 // Test1: We test accuracy of integrals of polynomials over reference
388 // quadrangle. Reference triangle {(-1,-1), (1,-1), (1,1), (-1,1)}.
389 //
390 // Test2: We consider simple mesh in meshFeTest.txt over square domain
391 // [0,1]^2 and test the accuracy of polynomials over square domain.
392 //
393
394 // get Quadrature
395 auto quad = fe::QuadElem(n);
396
397 //
398 // Test 1
399 //
400 size_t error_test_1 = 0;
401 {
402 // T1 (reference quadrangle)
403 // get quad points at reference triangle
404 std::vector<util::Point> nodes = {
405 util::Point(-1., -1., 0.), util::Point(1., -1., 0.),
406 util::Point(1., 1., 0.), util::Point(-1., 1., 0.)};
407 std::vector<fe::QuadData> qds = quad.getQuadPoints(nodes);
408 double sum = 0.;
409 for (auto qd : qds)
410 sum += qd.d_w;
411
412 if (std::abs(sum - 4.0) > tol) {
413 std::cout << "Error in order = " << n
414 << ". Sum of quad weights is not "
415 "equal to area of reference "
416 "quadrangle.\n";
417 error_test_1++;
418 }
419
420 //
421 // test the exactness of integration for polynomial
422 //
423 for (size_t i = 0; i <= 2 * n - 1; i++)
424 for (size_t j = 0; j <= 2 * n - 1; j++) {
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 // test integration of polynomial f(s,t) = s^i t^j
433 // get the exact integration
434 double I_exact = test::getExactIntegrationRefQuad(i, j);
435 if (!checkRefIntegration(n, i, j, qds, I_exact))
436 error_test_1++;
437
438 //
439 // when {(-1,1), (-1,-1), (1,-1), (1,1)}
440 //
441 nodes = {util::Point(-1., 1., 0.), util::Point(-1., -1., 0.),
442 util::Point(1., -1., 0.), util::Point(1., 1., 0.)};
443 qds = quad.getQuadPoints(nodes);
444 //
445 // After changing the order of vertices, we have got a new
446 // triangle which is in coordinate system (x,y) and we are
447 // integrating function f(x,y) = x^i y^j
448 //
449 // The quad data we have got is such that quad point is in (x,y)
450 // coordinate, weight is such that determinant of the Jacobian is
451 // included in the weight.
452 //
453 // Thus the following method for I_approx is correct.
454 if (!checkRefIntegration(n, i, j, qds, I_exact))
455 error_test_1++;
456
457 //
458 // when {(1,1), (-1,1), (-1,-1), (1,-1)}
459 //
460 nodes = {util::Point(1., 1., 0.), util::Point(-1., 1., 0.),
461 util::Point(-1., -1., 0.), util::Point(1., -1., 0.)};
462 qds = quad.getQuadPoints(nodes);
463 if (!checkRefIntegration(n, i, j, qds, I_exact))
464 error_test_1++;
465
466 //
467 // when {(1,-1), (1,1), (-1,1), (-1,-1)}
468 //
469 nodes = {util::Point(1., -1., 0.), util::Point(1., 1., 0.),
470 util::Point(-1., 1., 0.), util::Point(-1., -1., 0.)};
471 qds = quad.getQuadPoints(nodes);
472 if (!checkRefIntegration(n, i, j, qds, I_exact))
473 error_test_1++;
474 }
475 } // Test 1
476
477 //
478 // Test 2
479 //
480 size_t error_test_2 = 0;
481 {
482 static std::vector<util::Point> nodes;
483 static std::vector<size_t> elements;
484 static size_t num_vertex = 4;
485 static size_t elem_type = util::vtk_type_quad;
486 static size_t num_elems = 0;
487 if (num_elems == 0) {
488 readNodes(filepath + "/quadMesh_nodes.csv", nodes);
489 num_elems = readElements(filepath + "/quadMesh_elements.csv", elem_type,
490 elements);
491 }
492
493 // loop over polynomials
494 for (size_t i = 0; i <= 2 * n - 1; i++)
495 for (size_t j = 0; j <= 2 * n - 1; j++) {
496
497 double I_exact = 1. / (double(i + 1) * double(j + 1));
498 double I_approx = 0.;
499 // loop over elements and compute I_approx
500 for (size_t e = 0; e < num_elems; e++) {
501 std::vector<util::Point> enodes = {
502 nodes[elements[num_vertex * e + 0]],
503 nodes[elements[num_vertex * e + 1]],
504 nodes[elements[num_vertex * e + 2]],
505 nodes[elements[num_vertex * e + 3]]};
506 std::vector<fe::QuadData> qds = quad.getQuadPoints(enodes);
507 for (auto qd : qds)
508 I_approx +=
509 qd.d_w * std::pow(qd.d_p.d_x, i) * std::pow(qd.d_p.d_y, j);
510 }
511
512 if (std::abs(I_exact - I_approx) > tol) {
513 std::cout << "Error in order = " << n
514 << ". Exact integration = " << I_exact
515 << " and approximate integration = " << I_approx
516 << " of polynomial of order (i = " << i << " + j = " << j
517 << ") = " << i + j << " over square domain [0,1]x[0,1] "
518 << "is not matching using quadrature points.\n";
519
520 error_test_2++;
521 }
522 }
523 }
524
525 if (n == 1) {
526 std::cout << "**********************************\n";
527 std::cout << "Quadrangle Quadrature Test\n";
528 std::cout << "**********************************\n";
529 }
530 std::cout << "Quad order = " << n << ". ";
531 std::cout << (error_test_1 == 0 ? "TEST 1 : PASS. " : "TEST 1 : FAIL. ");
532 std::cout << (error_test_2 == 0 ? "TEST 2 : PASS. \n" : "TEST 2 : FAIL. \n");
533}
534
535void test::testTriElemTime(size_t n, size_t N) {
536
537 // get Quadrature
538 auto quad = fe::TriElem(n);
539
540 //
541 // Test 1
542 //
543 std::vector<util::Point> nodes = {util::Point(2., 2., 0.),
544 util::Point(4., 2., 0.),
545 util::Point(2., 4., 0.)};
546 size_t num_vertex = 3;
547 // std::vector<size_t> elements;
548 // for (size_t i = 0; i < 3 * N; i++)
549 // elements.emplace_back(i % 2);
550 std::vector<std::vector<size_t>> elements;
551 for (size_t i = 0; i < N; i++)
552 elements.emplace_back(std::vector<size_t>{0, 1, 2});
553
554 auto t11 = steady_clock::now();
555 // method 1: Compute quad points at each call
556 // loop over elements and compute I_approx
557 double sum = 0.;
558 for (size_t e = 0; e < N; e++) {
559 // std::vector<util::Point> enodes = {nodes[elements[num_vertex * e +
560 // 0]],
561 // nodes[elements[num_vertex * e +
562 // 1]], nodes[elements[num_vertex * e
563 // + 2]]};
564 std::vector<util::Point> enodes = {
565 nodes[elements[e][0]], nodes[elements[e][1]], nodes[elements[e][2]]};
566 std::vector<fe::QuadData> qds = quad.getQuadPoints(enodes);
567 for (auto qd : qds)
568 sum += qd.d_w * (qd.d_shapes[0] + qd.d_shapes[1] + qd.d_shapes[2]);
569 }
570 auto t12 = steady_clock::now();
571
572 // method 2: Compute quad points in the beginning and use it when needed
573 size_t num_quad_pts = 0;
574 std::vector<fe::QuadData> quad_data;
575 for (size_t e = 0; e < N; e++) {
576 std::vector<fe::QuadData> qds = quad.getQuadPoints(nodes);
577 if (e == 0)
578 num_quad_pts = qds.size();
579 for (auto qd : qds)
580 quad_data.emplace_back(qd);
581 }
582
583 auto t21 = steady_clock::now();
584 sum = 0.;
585 for (size_t e = 0; e < N; e++) {
586 for (size_t q = 0; q < num_quad_pts; q++) {
587 fe::QuadData qd = quad_data[e * num_quad_pts + q];
588 sum += qd.d_w * (qd.d_shapes[0] + qd.d_shapes[1] + qd.d_shapes[2]);
589 }
590 }
591 auto t22 = steady_clock::now();
592
593 if (n == 1 and N == 1000) {
594 std::cout << "**********************************\n";
595 std::cout << "Quadrature Time Efficiency Test\n";
596 std::cout << "**********************************\n";
597 }
598 std::cout << "Quad order = " << n << ". Num Elements = " << N << ".\n ";
599 double dt_1 = util::methods::timeDiff(t12, t11, "seconds");
600 double dt_2 = util::methods::timeDiff(t21, t22, "seconds");
601 double perc = (dt_1 - dt_2) * 100. / dt_2;
602 double qpt_mem = 13 * sizeof(double);
603 double mem2 = double(quad_data.capacity() * qpt_mem) / double(1000000);
604 std::cout << " dt1 = " << dt_1 << ", dt2 = " << dt_2 << ", perc = " << perc
605 << ". Mem saved = " << mem2 << " MB.\n";
606}
607
608void test::testTetElem(size_t n, std::string filepath) {
609
610 //
611 // Test1: We test accuracy of integrals of polynomials over reference
612 // tetrahedron. Reference element {(0,0,0), (1,0,0), (0,1,0), (0,0,1)}.
613 //
614 // Test2: We consider simple mesh in meshFeTest.txt over cubic domain
615 // [0,1]^3 and test the accuracy of polynomials over cubic domain.
616 //
617
618 // get Quadrature
619 auto quad = fe::TetElem(n);
620
621 //
622 // Test 1
623 //
624 size_t error_test_1 = 0;
625 {
626 // T1 (reference triangle)
627 // get quad points at reference triangle
628 std::vector<util::Point> nodes = {util::Point(), util::Point(1., 0., 0.),
629 util::Point(0., 1., 0.),
630 util::Point(0., 0., 1.)};
631 std::vector<fe::QuadData> qds = quad.getQuadPoints(nodes);
632
633 double sum = 0.;
634 for (auto qd : qds)
635 sum += qd.d_w;
636
637 if (std::abs(sum - 1. / 6.) > tol) {
638 std::cout << "Error in order = " << n
639 << ". Sum of quad weights is not "
640 "equal to volume of reference "
641 "tetrahedron.\n";
642 error_test_1++;
643 }
644
645 //
646 // test the exactness of integration for polynomial
647 //
648 for (size_t i = 0; i <= n; i++)
649 for (size_t j = 0; j <= n; j++)
650 for (size_t k = 0; k <= n; k++) {
651
652 if (i + j + k > n)
653 continue;
654
655 //
656 // +ve order of indices are:
657 // {0,1,2,3}; {1,2,0,3}; {2,3,0,1}; {0,3,1,2}
658
659 //
660 // when {(0,0,0), (1,0,0), (0,1,0), (0,0,1)}
661 //
662 nodes = {util::Point(), util::Point(1., 0., 0.),
663 util::Point(0., 1., 0.), util::Point(0., 0., 1.)};
664 qds = quad.getQuadPoints(nodes);
665 // test integration of polynomial f(s,t) = s^i t^j
666 // get the exact integration
667 debug_id = 0;
668 double I_exact = test::getExactIntegrationRefTet(i, j, k);
669 if (!checkRefIntegration(n, i, j, k, qds, I_exact)) {
670 error_test_1++;
671 }
672
673 //
674 // when vertices are {(1,0,0), (0,1,0), (0,0,1), (0,0,0)}
675 //
676 nodes = {util::Point(1., 0., 0.), util::Point(0., 1., 0.),
677 util::Point(0., 0., 0.), util::Point(0., 0., 1.)};
678 qds = quad.getQuadPoints(nodes);
679 //
680 // After changing the order of vertices, we have got a new
681 // triangle which is in coordinate system (x,y) and we are
682 // integrating function f(x,y) = x^i y^j
683 //
684 // The quad data we have got is such that quad point is in (x,y)
685 // coordinate, weight is such that determinant of the Jacobian is
686 // included in the weight.
687 //
688 // Thus the following method for I_approx is correct.
689 debug_id = 1;
690 if (!checkRefIntegration(n, i, j, k, qds, I_exact))
691 error_test_1++;
692
693 //
694 // when vertices are {(0,1,0), (0,0,1), (0,0,0), (1,0,0)}
695 //
696 nodes = {util::Point(0., 1., 0.), util::Point(0., 0., 1.),
697 util::Point(0., 0., 0.), util::Point(1., 0., 0.)};
698 qds = quad.getQuadPoints(nodes);
699 debug_id = 2;
700 if (!checkRefIntegration(n, i, j, k, qds, I_exact))
701 error_test_1++;
702
703 //
704 // when vertices are {(0,0,1), (0,0,0), (1,0,0), (0,1,0)}
705 //
706 nodes = {util::Point(0., 0., 0.), util::Point(0., 0., 1.),
707 util::Point(1., 0., 0.), util::Point(0., 1., 0.)};
708 qds = quad.getQuadPoints(nodes);
709 debug_id = 3;
710 if (!checkRefIntegration(n, i, j, k, qds, I_exact))
711 error_test_1++;
712 }
713 } // Test 1
714
715 //
716 // Test 2
717 //
718 size_t error_test_2 = 0;
719 if (false) {
720 static std::vector<util::Point> nodes;
721 static std::vector<size_t> elements;
722 static size_t num_vertex = 4;
723 static size_t elem_type = util::vtk_type_tetra;
724 static size_t num_elems = 0;
725 if (num_elems == 0) {
726 readNodes(filepath + "tetMesh_nodes.csv", nodes);
727 num_elems = readElements(filepath + "tetMesh_elements.csv", elem_type,
728 elements);
729 }
730
731 // loop over polynomials
732 for (size_t i = 0; i <= n; i++)
733 for (size_t j = 0; j <= n; j++)
734 for (size_t k = 0; k <= n; k++) {
735
736 if (i + j + k > n)
737 continue;
738
739 double I_exact = 1. / (double(i + 1) * double(j + 1) * double(k + 1));
740 double I_approx = 0.;
741 // loop over elements and compute I_approx
742 for (size_t e = 0; e < num_elems; e++) {
743 std::vector<util::Point> enodes = {
744 nodes[elements[num_vertex * e + 0]],
745 nodes[elements[num_vertex * e + 1]],
746 nodes[elements[num_vertex * e + 2]],
747 nodes[elements[num_vertex * e + 3]]};
748 std::vector<fe::QuadData> qds = quad.getQuadPoints(enodes);
749 for (auto qd : qds) {
750 I_approx += qd.d_w * std::pow(qd.d_p.d_x, i) *
751 std::pow(qd.d_p.d_y, j) * std::pow(qd.d_p.d_z, k);
752
753 if (false) {
754
755 std::cout << "Print " << i << " " << j << " " << k << "\n";
756 std::cout << util::io::printStr(enodes) << "\n";
757 std::vector<size_t> enode_ids = {
758 elements[num_vertex * e + 0],
759 elements[num_vertex * e + 1],
760 elements[num_vertex * e + 2],
761 elements[num_vertex * e + 3]};
762 std::cout << util::io::printStr(enode_ids) << "\n";
763 std::cout << qd.printStr() << "\n";
764 }
765
766 }
767 }
768
769 if (std::abs(I_exact - I_approx) > tol) {
770 std::cout << "Error in order = " << n
771 << ". Exact integration = " << I_exact
772 << " and approximate integration = " << I_approx
773 << " of polynomial of order (i = " << i << " + j = " << j
774 << " + k = " << k << ") = " << i + j + k
775 << " over cubic domain [0,1]x[0,1]x[0,1] "
776 << "is not matching using quadrature points.\n";
777
778 error_test_2++;
779 }
780 }
781 }
782
783 if (n == 1) {
784 std::cout << "**********************************\n";
785 std::cout << "Tetrahedron Quadrature Test\n";
786 std::cout << "**********************************\n";
787 }
788 std::cout << "Quad order = " << n << ". ";
789 if (error_test_1 == 0)
790 std::cout << "TEST 1 : PASS. ";
791 else
792 std::cout << "TEST 1 : FAIL. ";
793 // if (error_test_2 == 0)
794 // std::cout << "TEST 2 : PASS. ";
795 // else
796 // std::cout << "TEST 2 : FAIL. ";
797 std::cout << "\n";
798}
799
800void test::testPatchTri() {
801 const double ptol = 1.0e-5;
802 fe::TriElem tri(1);
803 const std::vector<util::Point> nodes = {
804 util::Point(0., 0., 0.), util::Point(1., 0., 0.),
805 util::Point(0., 1., 0.)};
806 const std::vector<util::Point> u = {
807 util::Point(0., 0., 0.), util::Point(2., 0., 0.),
808 util::Point(0., 3., 0.)};
809 auto qds = tri.getQuadDatas(nodes);
810 size_t nfail = 0;
811 for (const auto &qd : qds) {
812 auto e = strainFromB(fe::B(qd.d_derShapes, 2), u);
813 if (std::abs(e(0, 0) - 2.) > ptol || std::abs(e(1, 1) - 3.) > ptol ||
814 std::abs(e(0, 1)) > ptol)
815 nfail++;
816 }
817 std::cout << "**********************************\n";
818 std::cout << "Patch test (triangle, linear u)\n";
819 std::cout << "**********************************\n";
820 if (nfail) {
821 std::cout << "PATCH TRI : FAIL.\n";
822 exit(EXIT_FAILURE);
823 }
824 std::cout << "PATCH TRI : PASS.\n";
825}
826
827void test::testPatchQuad() {
828 const double ptol = 1.0e-5;
829 fe::QuadElem quad(1);
830 const std::vector<util::Point> nodes = {
831 util::Point(0., 0., 0.), util::Point(1., 0., 0.),
832 util::Point(1., 1., 0.), util::Point(0., 1., 0.)};
833 const std::vector<util::Point> u = {
834 util::Point(0., 0., 0.), util::Point(2., 0., 0.),
835 util::Point(2., 3., 0.), util::Point(0., 3., 0.)};
836 auto qds = quad.getQuadDatas(nodes);
837 size_t nfail = 0;
838 for (const auto &qd : qds) {
839 auto e = strainFromB(fe::B(qd.d_derShapes, 2), u);
840 if (std::abs(e(0, 0) - 2.) > ptol || std::abs(e(1, 1) - 3.) > ptol ||
841 std::abs(e(0, 1)) > ptol)
842 nfail++;
843 }
844 std::cout << "**********************************\n";
845 std::cout << "Patch test (quad, linear u)\n";
846 std::cout << "**********************************\n";
847 if (nfail) {
848 std::cout << "PATCH QUAD : FAIL.\n";
849 exit(EXIT_FAILURE);
850 }
851 std::cout << "PATCH QUAD : PASS.\n";
852}
853
854void test::testPatchTet() {
855 const double ptol = 1.0e-5;
856 fe::TetElem tet(1);
857 const std::vector<util::Point> nodes = {
858 util::Point(0., 0., 0.), util::Point(1., 0., 0.),
859 util::Point(0., 1., 0.), util::Point(0., 0., 1.)};
860 const std::vector<util::Point> u = {
861 util::Point(0., 0., 0.), util::Point(2., 0., 0.),
862 util::Point(0., 3., 0.), util::Point(0., 0., 4.)};
863 auto qds = tet.getQuadDatas(nodes);
864 size_t nfail = 0;
865 for (const auto &qd : qds) {
866 auto e = strainFromB(fe::B(qd.d_derShapes, 3), u);
867 if (std::abs(e(0, 0) - 2.) > ptol || std::abs(e(1, 1) - 3.) > ptol ||
868 std::abs(e(2, 2) - 4.) > ptol || std::abs(e(0, 1)) > ptol ||
869 std::abs(e(0, 2)) > ptol || std::abs(e(1, 2)) > ptol)
870 nfail++;
871 }
872 std::cout << "**********************************\n";
873 std::cout << "Patch test (tet, linear u)\n";
874 std::cout << "**********************************\n";
875 if (nfail) {
876 std::cout << "PATCH TET : FAIL.\n";
877 exit(EXIT_FAILURE);
878 }
879 std::cout << "PATCH TET : PASS.\n";
880}
881
882void test::testPatchTriDistorted() {
883 const double ptol = 1.0e-5;
884 fe::TriElem tri(1);
885 const std::vector<util::Point> nodes = {
886 util::Point(0., 0., 0.), util::Point(1.3, 0.2, 0.),
887 util::Point(0.15, 1.1, 0.)};
888 const std::vector<util::Point> u = {
889 util::Point(0., 0., 0.), util::Point(2.6, 0.6, 0.),
890 util::Point(0.3, 3.3, 0.)};
891 auto qds = tri.getQuadDatas(nodes);
892 size_t nfail = 0;
893 for (const auto &qd : qds) {
894 auto e = strainFromB(fe::B(qd.d_derShapes, 2), u);
895 if (std::abs(e(0, 0) - 2.) > ptol || std::abs(e(1, 1) - 3.) > ptol ||
896 std::abs(e(0, 1)) > ptol)
897 nfail++;
898 }
899 std::cout << "**********************************\n";
900 std::cout << "Patch test (distorted triangle, linear u)\n";
901 std::cout << "**********************************\n";
902 if (nfail) {
903 std::cout << "PATCH TRI DISTORTED : FAIL.\n";
904 exit(EXIT_FAILURE);
905 }
906 std::cout << "PATCH TRI DISTORTED : PASS.\n";
907}
fe::B::dim
int dim() const
Definition bMatrix.h:31
fe::B::nDof
int nDof() const
Definition bMatrix.h:30
fe::B::nStrain
int nStrain() const
Definition bMatrix.h:29
fe::BaseElem::getQuadDatas
std::vector< fe::QuadData > getQuadDatas(const std::vector< util::Point > &nodes)
Get quadrature data mapped to the physical element (N, dN/dx, x, w, J). Shared isoparametric map; typ...
Definition baseElem.cpp:61
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:29
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:70
anonymous_namespace{testFeLib.cpp}::readNodes
void readNodes(const std::string &filename, std::vector< util::Point > &nodes)
Definition testFeLib.cpp:60
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:109
anonymous_namespace{testFeLib.cpp}::strainFromB
util::SymMatrix3 strainFromB(const fe::B &B, const std::vector< util::Point > &u)
Definition testFeLib.cpp:31
anonymous_namespace{testFeLib.cpp}::debug_id
int debug_id
Definition testFeLib.cpp:28
test::testPatchTriDistorted
void testPatchTriDistorted()
Perform test on quadrature points on line elements (NOT IMPLEMENTED)
Definition testFeLib.cpp:882
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:535
test::getExactIntegrationRefTri
double getExactIntegrationRefTri(size_t alpha, size_t beta)
Computes integration of polynomial exactly over reference triangle.
Definition testFeLib.cpp:178
test::getNChooseR
double getNChooseR(size_t n, size_t r)
Computes "n choose r".
Definition testFeLib.cpp:166
test::getExactIntegrationRefQuad
double getExactIntegrationRefQuad(size_t alpha, size_t beta)
Computes integration of polynomial exactly over reference quadrangle.
Definition testFeLib.cpp:194
test::testLineElem
void testLineElem(size_t n, std::string filepath)
Perform test on quadrature points on line elements (NOT IMPLEMENTED)
Definition testFeLib.cpp:229
test::testTetElem
void testTetElem(size_t n, std::string filepath)
Perform test on quadrature points on tetrahedral elements.
Definition testFeLib.cpp:608
test::testQuadElem
void testQuadElem(size_t n, std::string filepath)
Perform test on quadrature points on quadrangle elements.
Definition testFeLib.cpp:384
test::testPatchTet
void testPatchTet()
Perform test on quadrature points on line elements (NOT IMPLEMENTED)
Definition testFeLib.cpp:854
test::testPatchQuad
void testPatchQuad()
Perform test on quadrature points on line elements (NOT IMPLEMENTED)
Definition testFeLib.cpp:827
test::testPatchTri
void testPatchTri()
Perform test on quadrature points on line elements (NOT IMPLEMENTED)
Definition testFeLib.cpp:800
test::getExactIntegrationRefTet
double getExactIntegrationRefTet(size_t alpha, size_t beta, size_t theta)
Computes integration of polynomial exactly over reference tetrahedral.
Definition testFeLib.cpp:203
test::testTriElem
void testTriElem(size_t n, std::string filepath)
Perform test on quadrature points on triangle elements.
Definition testFeLib.cpp:231
util::io::printStr
std::string printStr(const T &msg, int nt=print_default_tab)
Returns formatted string for output.
Definition io.h:96
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 vecMethods.h:309
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
util::SymMatrix3
A structure to represent 3d matrices.
Definition matrix.h:258