PeriDEM 0.2.0
PeriDEM -- Peridynamics-based high-fidelity model for granular media
Loading...
Searching...
No Matches
tetElem.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 "tetElem.h"
12#include "util/function.h"
13#include "util/matrix.h"
14#include <iostream>
15#include <util/io.h>
16
17#include "util/feElementDefs.h" // global definition of elements
18
19namespace {
20
21void checkPoint(const std::vector<double> &p, const std::vector<util::Point> &nodes) {
22
23 // check to see if p is in reference tet element
24 bool check = false;
25 if (util::isLess(p[0], -1.0E-5) || util::isLess(p[1], -1.0E-5) || util::isLess(p[2], -1.0E-5) ||
26 util::isGreater(p[0], 1. + 1.0E-5) ||
27 util::isGreater(p[1], 1. + 1.0E-5) ||
28 util::isGreater(p[2], 1. + 1.0E-5)) {
29
30 check = true;
31 }
32
33 if (!check) {
34
35 // check if projection of point in x, y, z plane is within the limit
36 if (util::isGreater(p[0], 1. + 1.0E-5 - p[1]))
37 check = true;
38
39 if (util::isGreater(p[1], 1. + 1.0E-5 - p[2]))
40 check = true;
41
42 if (util::isGreater(p[2], 1. + 1.0E-5 - p[0]))
43 check = true;
44 }
45
46 if (check) {
47 std::cerr << "Error: Point p = ("
48 << p[0] << ", " << p[1] << ", " << p[2]
49 << ") does not belong to reference tet element = {("
50 << nodes[0].d_x << ", " << nodes[0].d_y << ", " << nodes[0].d_z
51 << "), ("
52 << nodes[1].d_x << "," << nodes[1].d_y << ", " << nodes[1].d_z
53 << "), ("
54 << nodes[2].d_x << "," << nodes[2].d_y << ", " << nodes[2].d_z
55 << "), ("
56 << nodes[3].d_x << "," << nodes[3].d_y << ", " << nodes[3].d_z
57 << ")}.\n"
58 << "Coordinates in reference element are: "
59 << "xi = " << p[0]
60 << ", eta = " << p[1]
61 << ", zeta = " << p[2] << "\n";
62 exit(1);
63 }
64}
65
66} // anonymous namespace
67
68fe::TetElem::TetElem(size_t order)
69 : fe::BaseElem(order, util::vtk_type_triangle) {
70
71 if (d_quadOrder > 3) {
72 std::cout << "Error: For linear tet element, we only support upto 3 quad "
73 "order approximation.\n";
74 exit(1);
75 }
76
77 // compute quad data
78 this->init();
79}
80
81double fe::TetElem::elemSize(const std::vector<util::Point> &nodes) {
82 // volume of tet element is (1/6) a * (b x c),
83 // where a = v2 - v1, b = v3 - v1, c = v4 - v1
84 auto a = nodes[1] - nodes[0];
85 auto b = nodes[2] - nodes[0];
86 auto c = nodes[3] - nodes[0];
87 return (1. / 6.) * a * (b.cross(c));
88}
89
90std::vector<double> fe::TetElem::getShapes(
91 const util::Point &p, const std::vector<util::Point> &nodes) {
92 return getShapes(mapPointToRefElem(p, nodes));
93}
94
95std::vector<std::vector<double>> fe::TetElem::getDerShapes(
96 const util::Point &p, const std::vector<util::Point> &nodes) {
97
98 // get derivatives of shape function in reference tet element
99 auto ders_ref = getDerShapes(mapPointToRefElem(p, nodes));
100
101 // get Jacobian and its determinant
102 std::vector<std::vector<double>> J;
103 auto detJ = getJacobian(p, nodes, &J);
104
105 auto J_inv = util::inv(J);
106
107 // to hold derivatives
108 std::vector<std::vector<double>> ders(ders_ref.size(),
109 std::vector<double>(3, 0.));
110
111 // grad N_i = J_inv * grad N_i^ref
112 for (size_t i = 0; i < 4; i++)
113 ders[i] = util::dot(J_inv, ders_ref[i]);
114
115 return ders;
116}
117
118std::vector<fe::QuadData> fe::TetElem::getQuadDatas(
119 const std::vector<util::Point> &nodes) {
120
121 // copy quad data associated to reference element
122 auto qds = d_quads;
123
124 // modify data
125 for (auto &qd : qds) {
126
127 // get Jacobian and determinant
128 qd.d_detJ = getJacobian(qd.d_p, nodes, &(qd.d_J));
129
130 // transform quad weight
131 qd.d_w *= qd.d_detJ;
132
133 // map point to triangle
134 qd.d_p = util::Point();
135 for (size_t i=0; i<4; i++) {
136 qd.d_p.d_x += qd.d_shapes[i] * nodes[i].d_x;
137 qd.d_p.d_y += qd.d_shapes[i] * nodes[i].d_y;
138 qd.d_p.d_z += qd.d_shapes[i] * nodes[i].d_z;
139 }
140
141 // get inverse of Jacobian
142 auto J_inv = util::inv(qd.d_J);
143
144 // to hold derivatives
145 std::vector<std::vector<double>> ders(qd.d_derShapes.size(),
146 std::vector<double>(3, 0.));
147
148 // grad N_i = J_inv * grad N_i^ref
149 for (size_t i = 0; i < 4; i++)
150 ders[i] = util::dot(J_inv, qd.d_derShapes[i]);
151
152 qd.d_derShapes = ders;
153 }
154
155 return qds;
156}
157
158std::vector<fe::QuadData> fe::TetElem::getQuadPoints(
159 const std::vector<util::Point> &nodes) {
160 // copy quad data associated to reference element
161 auto qds = d_quads;
162
163 // modify data
164 for (auto &qd : qds) {
165
166 // transform quad weight
167 qd.d_w *= getJacobian(qd.d_p, nodes, nullptr);
168
169 // map point to triangle
170 qd.d_p = util::Point();
171 for (size_t i=0; i<4; i++) {
172 qd.d_p.d_x += qd.d_shapes[i] * nodes[i].d_x;
173 qd.d_p.d_y += qd.d_shapes[i] * nodes[i].d_y;
174 qd.d_p.d_z += qd.d_shapes[i] * nodes[i].d_z;
175 }
176 }
177
178 return qds;
179}
180
181std::vector<double> fe::TetElem::getShapes(const util::Point &p) {
182 // N1 = 1 - xi - eta - zeta, N2 = xi, N3 = eta, N4 = zeta
183 return std::vector<double>{1. - p.d_x - p.d_y - p.d_z, p.d_x, p.d_y, p.d_z};
184}
185
186std::vector<std::vector<double>> fe::TetElem::getDerShapes(
187 const util::Point &p) {
188
189 // d N1/d xi = -1, d N1/d eta = -1, d N1/d zeta = -1,
190 // d N2/ d xi = 1, d N2/d eta = 0, d N2/d zeta = 0,
191 // d N3/ d xi = 0, d N3/d eta = 1, d N3/d zeta = 0,
192 // d N4/ d xi = 0, d N4/d eta = 0, d N4/d zeta = 1,
193 std::vector<std::vector<double>> r;
194 r.push_back(std::vector<double>{-1., -1., -1.});
195 r.push_back(std::vector<double>{1., 0., 0.});
196 r.push_back(std::vector<double>{0., 1., 0.});
197 r.push_back(std::vector<double>{0., 0., 1.});
198
199 return r;
200}
201
202util::Point
203fe::TetElem::mapPointToRefElem(
204 const util::Point &p, const std::vector<util::Point> &nodes) {
205
206 // get Jacobian matrix and compute its transpose
207 std::vector<std::vector<double>> J(3, std::vector<double>(3, 0.));
208 auto detJ = getJacobian(p, nodes, &J);
209
210 // get transpose of Jacobian
211 auto B = util::transpose(J);
212 auto detB = detJ;
213
214 // get inverse of B
215 auto B_inv = util::inv(B);
216
217 // get vector from first vertex to point p
218 std::vector<double> vec_p = {p.d_x - nodes[0].d_x,
219 p.d_y - nodes[0].d_y,
220 p.d_z - nodes[0].d_z};
221 // multiply B_inv to vector to transform point
222 auto p_ref = util::dot(B_inv, vec_p);
223
224 // check point
225 checkPoint(p_ref, nodes);
226
227 if (util::isLess(p_ref[0], 0.)) p_ref[0] = 0.;
228 if (util::isLess(p_ref[1], 0.)) p_ref[1] = 0.;
229 if (util::isLess(p_ref[2], 0.)) p_ref[2] = 0.;
230 if (util::isGreater(p_ref[0], 1.)) p_ref[0] = 1.;
231 if (util::isGreater(p_ref[1], 1.)) p_ref[1] = 1.;
232 if (util::isGreater(p_ref[2], 1.)) p_ref[2] = 1.;
233
234 return util::Point(p_ref);
235}
236
237double fe::TetElem::getJacobian(const util::Point &p,
238 const std::vector<util::Point> &nodes,
239 std::vector<std::vector<double>> *J) {
240 if (J != nullptr) {
241 J->resize(3);
242 (*J)[0] = std::vector<double>{nodes[1].d_x - nodes[0].d_x,
243 nodes[1].d_y - nodes[0].d_y,
244 nodes[1].d_z - nodes[0].d_z};
245 (*J)[1] = std::vector<double>{nodes[2].d_x - nodes[0].d_x,
246 nodes[2].d_y - nodes[0].d_y,
247 nodes[2].d_z - nodes[0].d_z};
248 (*J)[2] = std::vector<double>{nodes[3].d_x - nodes[0].d_x,
249 nodes[3].d_y - nodes[0].d_y,
250 nodes[3].d_z - nodes[0].d_z};
251
252 return util::det(*J);
253 } else {
254 std::vector<std::vector<double>> J_local;
255 J_local.resize(3);
256 J_local[0] = std::vector<double>{nodes[1].d_x - nodes[0].d_x,
257 nodes[1].d_y - nodes[0].d_y,
258 nodes[1].d_z - nodes[0].d_z};
259 J_local[1] = std::vector<double>{nodes[2].d_x - nodes[0].d_x,
260 nodes[2].d_y - nodes[0].d_y,
261 nodes[2].d_z - nodes[0].d_z};
262 J_local[2] = std::vector<double>{nodes[3].d_x - nodes[0].d_x,
263 nodes[3].d_y - nodes[0].d_y,
264 nodes[3].d_z - nodes[0].d_z};
265
266 return util::det(J_local);
267 }
268}
269
270void fe::TetElem::init() {
271 //
272 // compute quad data for reference triangle with vertex at
273 // (0,0), (1,0), (0,1)
274 //
275
276 if (!d_quads.empty()) return;
277
278 // no point in zeroth order
279 if (d_quadOrder == 0) {
280 d_quads.resize(0);
281 }
282
283 // 3x3 identity matrix
284 std::vector<std::vector<double>> ident_mat;
285 ident_mat.push_back(std::vector<double>{1., 0., 0.});
286 ident_mat.push_back(std::vector<double>{0., 1., 0.});
287 ident_mat.push_back(std::vector<double>{0., 0., 1.});
288
289 //
290 // These datas are from LibMesh code
291 // See: https://libmesh.github.io/doxygen/quadrature__gauss__3D_8C_source.html
292
293 //
294 // first order quad points for triangle
295 //
296 if (d_quadOrder == 1) {
297 d_quads.clear();
298 fe::QuadData qd;
299 qd.d_w = 1. / 6.;
300 qd.d_p = util::Point(1. / 4., 1. / 4., 1. / 4.);
301 // N1 = 1 - xi - eta, N2 = xi, N3 = eta
302 qd.d_shapes = getShapes(qd.d_p);
303 // d N1/d xi = -1, d N1/d eta = -1, d N2/ d xi = 1, d N2/d eta = 0,
304 // d N3/ d xi = 0, d N3/d eta = 1
305 qd.d_derShapes = getDerShapes(qd.d_p);
306 qd.d_J = ident_mat;
307 qd.d_detJ = 1.;
308 d_quads.push_back(qd);
309 }
310
311 //
312 // second order quad points for triangle
313 //
314 if (d_quadOrder == 2) {
315 d_quads.clear();
316 fe::QuadData qd;
317
318 double w = 1. / 24.;
319 double a = 0.585410196624969;
320 double b = 0.138196601125011;
321 // point 1
322 qd.d_w = w;
323 qd.d_p = util::Point(a, b, b);
324 qd.d_shapes = getShapes(qd.d_p);
325 qd.d_derShapes = getDerShapes(qd.d_p);
326 qd.d_J = ident_mat;
327 qd.d_detJ = 1.;
328 d_quads.push_back(qd);
329 // point 2
330 qd.d_w = w;
331 qd.d_p = util::Point(b, a, b);
332 qd.d_shapes = getShapes(qd.d_p);
333 qd.d_derShapes = getDerShapes(qd.d_p);
334 qd.d_J = ident_mat;
335 qd.d_detJ = 1.;
336 d_quads.push_back(qd);
337 // point 3
338 qd.d_w = w;
339 qd.d_p = util::Point(b, b, a);
340 qd.d_shapes = getShapes(qd.d_p);
341 qd.d_derShapes = getDerShapes(qd.d_p);
342 qd.d_J = ident_mat;
343 qd.d_detJ = 1.;
344 d_quads.push_back(qd);
345 // point 4
346 qd.d_w = w;
347 qd.d_p = util::Point(b, b, b);
348 qd.d_shapes = getShapes(qd.d_p);
349 qd.d_derShapes = getDerShapes(qd.d_p);
350 qd.d_J = ident_mat;
351 qd.d_detJ = 1.;
352 d_quads.push_back(qd);
353 }
354
355 //
356 // third order quad points for triangle
357 //
358 if (d_quadOrder == 3) {
359 d_quads.clear();
360 fe::QuadData qd;
361
362 double w1 = -2. / 15.;
363 double w2 = 0.075;
364
365 double a = 0.25;
366 double b = 0.5;
367 double c = 1. / 6.;
368
369 // point 1
370 qd.d_w = w1;
371 qd.d_p = util::Point(a, a, a);
372 qd.d_shapes = getShapes(qd.d_p);
373 qd.d_derShapes = getDerShapes(qd.d_p);
374 qd.d_J = ident_mat;
375 qd.d_detJ = 1.;
376 d_quads.push_back(qd);
377 // point 2
378 qd.d_w = w2;
379 qd.d_p = util::Point(b, c, c);
380 qd.d_shapes = getShapes(qd.d_p);
381 qd.d_derShapes = getDerShapes(qd.d_p);
382 qd.d_J = ident_mat;
383 qd.d_detJ = 1.;
384 d_quads.push_back(qd);
385 // point 3
386 qd.d_w = w2;
387 qd.d_p = util::Point(c, b, c);
388 qd.d_shapes = getShapes(qd.d_p);
389 qd.d_derShapes = getDerShapes(qd.d_p);
390 qd.d_J = ident_mat;
391 qd.d_detJ = 1.;
392 d_quads.push_back(qd);
393 // point 4
394 qd.d_w = w2;
395 qd.d_p = util::Point(c, c, b);
396 qd.d_shapes = getShapes(qd.d_p);
397 qd.d_derShapes = getDerShapes(qd.d_p);
398 qd.d_J = ident_mat;
399 qd.d_detJ = 1.;
400 d_quads.push_back(qd);
401 // point 5
402 qd.d_w = w2;
403 qd.d_p = util::Point(c, c, c);
404 qd.d_shapes = getShapes(qd.d_p);
405 qd.d_derShapes = getDerShapes(qd.d_p);
406 qd.d_J = ident_mat;
407 qd.d_detJ = 1.;
408 d_quads.push_back(qd);
409 }
410}
fe::BaseElem
A base class which provides methods to map points to/from reference element and to compute quadrature...
Definition baseElem.h:84
fe::BaseElem::d_quadOrder
size_t d_quadOrder
Order of quadrature point integration approximation.
Definition baseElem.h:243
fe::TetElem::init
void init() override
Compute the quadrature points for triangle element.
Definition tetElem.cpp:270
fe::TetElem::getDerShapes
std::vector< std::vector< double > > getDerShapes(const util::Point &p, const std::vector< util::Point > &nodes) override
Returns the values of derivative of shape function at point p.
Definition tetElem.cpp:95
fe::TetElem::TetElem
TetElem(size_t order)
Constructor.
Definition tetElem.cpp:68
fe::TetElem::getQuadDatas
std::vector< fe::QuadData > getQuadDatas(const std::vector< util::Point > &nodes) override
Get vector of quadrature data.
Definition tetElem.cpp:118
fe::TetElem::getQuadPoints
std::vector< fe::QuadData > getQuadPoints(const std::vector< util::Point > &nodes) override
Get vector of quadrature data.
Definition tetElem.cpp:158
fe::TetElem::elemSize
double elemSize(const std::vector< util::Point > &nodes) override
Returns the volume of element.
Definition tetElem.cpp:81
fe::TetElem::getJacobian
double getJacobian(const util::Point &p, const std::vector< util::Point > &nodes, std::vector< std::vector< double > > *J) override
Computes the Jacobian of map .
Definition tetElem.cpp:237
fe::TetElem::mapPointToRefElem
util::Point mapPointToRefElem(const util::Point &p, const std::vector< util::Point > &nodes) override
Maps point p in a given element to the reference element.
Definition tetElem.cpp:203
fe::TetElem::getShapes
std::vector< double > getShapes(const util::Point &p, const std::vector< util::Point > &nodes) override
Returns the values of shape function at point p.
Definition tetElem.cpp:90
anonymous_namespace{tetElem.cpp}::checkPoint
void checkPoint(const std::vector< double > &p, const std::vector< util::Point > &nodes)
Definition tetElem.cpp:21
fe
Collection of methods and data related to finite element and mesh.
Definition baseElem.h:17
util
Collection of methods useful in simulation.
Definition feElementDefs.h:14
util::isGreater
bool isGreater(const double &a, const double &b)
Returns true if a > b.
Definition function.cpp:15
util::det
double det(const std::vector< std::vector< double > > &m)
Computes the determinant of matrix.
Definition matrix.cpp:75
util::inv
std::vector< std::vector< double > > inv(const std::vector< std::vector< double > > &m)
Computes the determinant of matrix.
Definition matrix.cpp:93
util::isLess
bool isLess(const double &a, const double &b)
Returns true if a < b.
Definition function.cpp:20
util::dot
std::vector< double > dot(const std::vector< std::vector< double > > &m, const std::vector< double > &v)
Computes the dot product between matrix and vector.
Definition matrix.cpp:38
util::transpose
std::vector< std::vector< double > > transpose(const std::vector< std::vector< double > > &m)
Computes the tranpose of matrix.
Definition matrix.cpp:56
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_derShapes
std::vector< std::vector< double > > d_derShapes
Value of derivative of shape functions at quad point p.
Definition quadData.h:50
fe::QuadData::d_w
double d_w
Quadrature weight.
Definition quadData.h:26
fe::QuadData::d_p
util::Point d_p
Quadrature point in 1-d, 2-d or 3-d.
Definition quadData.h:29
fe::QuadData::d_J
std::vector< std::vector< double > > d_J
Jacobian of the map from reference element to the element.
Definition quadData.h:58
fe::QuadData::d_detJ
double d_detJ
Determinant of the Jacobian of the map from reference element to the element.
Definition quadData.h:64
util::Point
A structure to represent 3d vectors.
Definition point.h:30
util::Point::d_y
double d_y
the y coordinate
Definition point.h:36
util::Point::d_z
double d_z
the z coordinate
Definition point.h:39
util::Point::d_x
double d_x
the x coordinate
Definition point.h:33