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));
30 const std::vector<util::Point> &nodes) {
31 return getShapes(mapPointToRefElem(p, nodes));
34std::vector<std::vector<double>>
36 const std::vector<util::Point> &nodes) {
38 auto ders_ref = getDerShapes(mapPointToRefElem(p, nodes));
41 std::vector<std::vector<double>> J;
42 auto detJ = getJacobian(p, nodes, &J);
47 for (
size_t i=0; i<3; i++) {
49 ders[i][0] = (ders_ref[i][0] * J[1][1] - ders_ref[i][1] * J[0][1]) / detJ;
51 ders[i][1] = (-ders_ref[i][0] * J[1][0] + ders_ref[i][1] * J[0][0]) / detJ;
57std::vector<fe::QuadData>
64 for (
auto &qd : qds) {
67 qd.d_detJ = getJacobian(qd.d_p, nodes, &(qd.d_J));
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;
79 std::vector<std::vector<double>> ders;
81 for (
size_t i=0; i<3; i++) {
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;
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;
89 ders.push_back(std::vector<double>{d1, d2});
91 qd.d_derShapes = ders;
97std::vector<fe::QuadData>
104 for (
auto &qd : qds) {
107 qd.d_w *= getJacobian(qd.d_p, nodes,
nullptr);
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;
122 return std::vector<double>{1. - p.
d_x - p.
d_y, p.
d_x, p.
d_y};
125std::vector<std::vector<double>>
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.});
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)) /
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)) /
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";
167 return {xi, eta, 0.};
171 const std::vector<util::Point> &nodes,
172 std::vector<std::vector<double>> *J) {
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};
181 return (*J)[0][0] * (*J)[1][1] - (*J)[0][1] * (*J)[1][0];
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);
195 if (!d_quads.empty())
199 if (d_quadOrder == 0) {
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.});
211 if (d_quadOrder == 1) {
223 d_quads.push_back(qd);
229 if (d_quadOrder == 2) {
239 d_quads.push_back(qd);
247 d_quads.push_back(qd);
255 d_quads.push_back(qd);
261 if (d_quadOrder == 3) {
271 d_quads.push_back(qd);
279 d_quads.push_back(qd);
287 d_quads.push_back(qd);
295 d_quads.push_back(qd);
301 if (d_quadOrder == 4) {
305 qd.
d_w = 0.5 * 0.22338158967801;
311 d_quads.push_back(qd);
313 qd.
d_w = 0.5 * 0.22338158967801;
319 d_quads.push_back(qd);
321 qd.
d_w = 0.5 * 0.22338158967801;
327 d_quads.push_back(qd);
329 qd.
d_w = 0.5 * 0.10995174365532;
335 d_quads.push_back(qd);
337 qd.
d_w = 0.5 * 0.10995174365532;
343 d_quads.push_back(qd);
345 qd.
d_w = 0.5 * 0.10995174365532;
351 d_quads.push_back(qd);
357 if (d_quadOrder == 5) {
361 qd.
d_w = 0.5 * 0.22500000000000;
367 d_quads.push_back(qd);
369 qd.
d_w = 0.5 * 0.13239415278851;
375 d_quads.push_back(qd);
377 qd.
d_w = 0.5 * 0.13239415278851;
383 d_quads.push_back(qd);
385 qd.
d_w = 0.5 * 0.13239415278851;
391 d_quads.push_back(qd);
393 qd.
d_w = 0.5 * 0.12593918054483;
399 d_quads.push_back(qd);
401 qd.
d_w = 0.5 * 0.12593918054483;
407 d_quads.push_back(qd);
409 qd.
d_w = 0.5 * 0.12593918054483;
415 d_quads.push_back(qd);
A base class which provides methods to map points to/from reference element and to compute quadrature...
double elemSize(const std::vector< util::Point > &nodes) override
Returns the area of element.
std::vector< fe::QuadData > getQuadPoints(const std::vector< util::Point > &nodes) override
Get vector of quadrature data.
double getJacobian(const util::Point &p, const std::vector< util::Point > &nodes, std::vector< std::vector< double > > *J) override
Computes the Jacobian of map .
std::vector< fe::QuadData > getQuadDatas(const std::vector< util::Point > &nodes) override
Get vector of quadrature data.
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.
std::vector< double > getShapes(const util::Point &p, const std::vector< util::Point > &nodes) override
Returns the values of shape function at point p.
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.
void init() override
Compute the quadrature points for triangle element.
TriElem(size_t order)
Constructor.
Collection of methods and data related to finite element and mesh.
Collection of methods useful in simulation.
bool isGreater(const double &a, const double &b)
Returns true if a > b.
bool isLess(const double &a, const double &b)
Returns true if a < b.
A struct to store the quadrature data. List of data are.
std::vector< double > d_shapes
Value of shape functions at quad point p.
std::vector< std::vector< double > > d_derShapes
Value of derivative of shape functions at quad point p.
double d_w
Quadrature weight.
util::Point d_p
Quadrature point in 1-d, 2-d or 3-d.
std::vector< std::vector< double > > d_J
Jacobian of the map from reference element to the element.
double d_detJ
Determinant of the Jacobian of the map from reference element to the element.
A structure to represent 3d vectors.
double d_y
the y coordinate
double d_x
the x coordinate