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