-
Notifications
You must be signed in to change notification settings - Fork 13
/
lgr_triangle.cpp
119 lines (113 loc) · 5.28 KB
/
lgr_triangle.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
#include <hpc_array.hpp>
#include <lgr_element_specific_inline.hpp>
#include <lgr_state.hpp>
#include <lgr_triangle.hpp>
namespace lgr {
void
initialize_triangle_V(state& s)
{
auto const element_nodes_to_nodes = s.elements_to_nodes.cbegin();
auto const nodes_to_x = s.x.cbegin();
auto const points_to_V = s.V.begin();
auto const elements_to_element_nodes = s.elements * s.nodes_in_element;
auto const elements_to_points = s.elements * s.points_in_element;
auto functor = [=] HPC_DEVICE(element_index const element) {
constexpr point_in_element_index fp(0);
auto const element_nodes = elements_to_element_nodes[element];
using l_t = node_in_element_index;
hpc::array<node_index, 3> nodes;
hpc::array<hpc::position<double>, 3> x;
for (int i = 0; i < 3; ++i) {
auto const node = element_nodes_to_nodes[element_nodes[l_t(i)]];
nodes[i] = node;
x[i] = nodes_to_x[node].load();
}
auto const area = triangle_area(x);
assert(area > 0.0);
auto const volume = area * hpc::length<double>(1.0);
points_to_V[elements_to_points[element][fp]] = volume;
};
hpc::for_each(hpc::device_policy(), s.elements, functor);
}
void
initialize_triangle_grad_N(state& s)
{
auto const element_nodes_to_nodes = s.elements_to_nodes.cbegin();
auto const nodes_to_x = s.x.cbegin();
auto const points_to_V = s.V.cbegin();
auto const point_nodes_to_grad_N = s.grad_N.begin();
auto const elements_to_element_nodes = s.elements * s.nodes_in_element;
auto const points_to_point_nodes = s.points * s.nodes_in_element;
auto const elements_to_points = s.elements * s.points_in_element;
auto functor = [=] HPC_DEVICE(element_index const element) {
constexpr point_in_element_index fp(0);
auto const element_nodes = elements_to_element_nodes[element];
auto const point = elements_to_points[element][fp];
auto const point_nodes = points_to_point_nodes[point];
using l_t = node_in_element_index;
hpc::array<hpc::position<double>, 3> x;
for (int i = 0; i < 3; ++i) {
auto const node = element_nodes_to_nodes[element_nodes[l_t(i)]];
x[i] = nodes_to_x[node].load();
}
auto const volume = points_to_V[point];
auto const area = volume / hpc::length<double>(1.0);
auto const grad_N = triangle_basis_gradients(x, area);
for (int i = 0; i < 3; ++i) {
point_nodes_to_grad_N[point_nodes[l_t(i)]] = grad_N[i];
}
};
hpc::for_each(hpc::device_policy(), s.elements, functor);
}
void
update_triangle_h_min_inball(input const&, state& s)
{
auto const point_nodes_to_grad_N = s.grad_N.cbegin();
auto const elements_to_h_min = s.h_min.begin();
auto const points_to_point_nodes = s.points * s.nodes_in_element;
auto const nodes_in_element = s.nodes_in_element;
auto const elements_to_points = s.elements * s.points_in_element;
auto functor = [=] HPC_DEVICE(element_index const element) {
/* find the radius of the inscribed circle.
first fun fact: the area of a triangle equals one half
times the radius of the inscribed circle times the perimeter
of the triangle, where the perimeter is the sum of its
edge lengths.
second fun fact: the magnitude of the gradient of the basis function
of a triangle's node is equal to the length of the opposite edge
divided by twice the triangle area
third fun fact: when solving for the radius, area cancels out
of the top and bottom of the division.
*/
constexpr point_in_element_index fp(0);
auto const point = elements_to_points[element][fp];
auto const point_nodes = points_to_point_nodes[point];
decltype(hpc::length<double>() / hpc::area<double>()) perimeter_over_twice_area = 0.0;
for (auto const i : nodes_in_element) {
auto const grad_N = point_nodes_to_grad_N[point_nodes[i]].load();
auto const edge_length_over_twice_area = norm(grad_N);
perimeter_over_twice_area += edge_length_over_twice_area;
}
auto const radius = 1.0 / perimeter_over_twice_area;
elements_to_h_min[element] = 2.0 * radius;
};
hpc::for_each(hpc::device_policy(), s.elements, functor);
}
void
update_triangle_h_art(state& s)
{
double const C_geom = std::sqrt(4.0 / std::sqrt(3.0));
auto const points_to_V = s.V.cbegin();
auto const elements_to_h_art = s.h_art.begin();
auto const elements_to_points = s.elements * s.points_in_element;
auto functor = [=] HPC_DEVICE(element_index const element) {
hpc::area<double> area = 0.0;
for (auto const point : elements_to_points[element]) {
area += (points_to_V[point] / hpc::length<double>(1.0));
}
auto const h_art = C_geom * sqrt(area);
elements_to_h_art[element] = h_art;
};
hpc::for_each(hpc::device_policy(), s.elements, functor);
}
} // namespace lgr