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