-
Notifications
You must be signed in to change notification settings - Fork 13
/
lgr_domain.cpp
201 lines (185 loc) · 7.26 KB
/
lgr_domain.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
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
#include <lgr_domain.hpp>
#include <lgr_input.hpp>
#include <lgr_state.hpp>
namespace lgr {
domain::~domain()
{
}
void
union_domain::add(std::unique_ptr<domain>&& uptr)
{
m_domains.push_back(std::move(uptr));
}
void
union_domain::mark(
hpc::device_array_vector<hpc::position<double>, node_index> const& points,
int const marker,
hpc::device_vector<int, node_index>* markers) const
{
for (auto const& uptr : m_domains) {
uptr->mark(points, marker, markers);
}
}
#ifdef HPC_ENABLE_STRONG_INDICES
void
union_domain::mark(
hpc::device_array_vector<hpc::position<double>, element_index> const& points,
material_index const marker,
hpc::device_vector<material_index, element_index>* markers) const
{
for (auto const& uptr : m_domains) {
uptr->mark(points, marker, markers);
}
}
#endif
void
union_domain::mark(
hpc::device_array_vector<hpc::position<double>, node_index> const& points,
material_index const marker,
hpc::device_vector<material_set, node_index>* markers) const
{
for (auto const& uptr : m_domains) {
uptr->mark(points, marker, markers);
}
}
template <class Index, class IsInFunctor>
HPC_NOINLINE void
collect_set(
hpc::counting_range<Index> const range,
IsInFunctor is_in_functor,
hpc::device_vector<Index, int>& set_items);
template <class Index, class IsInFunctor>
void
collect_set(
hpc::counting_range<Index> const range,
IsInFunctor is_in_functor,
hpc::device_vector<Index, int>& set_items)
{
hpc::device_vector<int, Index> offsets(range.size());
hpc::transform_inclusive_scan(hpc::device_policy(), range, offsets, hpc::plus<int>(), is_in_functor);
int const set_size = hpc::transform_reduce(hpc::device_policy(), range, int(0), hpc::plus<int>(), is_in_functor);
set_items.resize(set_size);
auto const set_items_to_items = set_items.begin();
auto const items_to_offsets = offsets.cbegin();
auto functor = [=] HPC_DEVICE(Index const item) {
if (is_in_functor(item) != 0) {
set_items_to_items[items_to_offsets[item] - 1] = item;
}
};
hpc::for_each(hpc::device_policy(), range, functor);
}
void
assign_element_materials(input const& in, state& s)
{
hpc::fill(hpc::device_policy(), s.material, material_index(0));
hpc::device_array_vector<hpc::position<double>, element_index> centroid_vector(s.elements.size());
auto const elements_to_element_nodes = s.elements * s.nodes_in_element;
auto const element_nodes_to_nodes = s.elements_to_nodes.cbegin();
auto const nodes_to_x = s.x.cbegin();
double const N = 1.0 / double(hpc::weaken(s.nodes_in_element.size()));
auto const elements_to_centroids = centroid_vector.begin();
auto centroid_functor = [=] HPC_DEVICE(element_index const element) {
auto centroid = hpc::position<double>::zero();
for (auto const element_node : elements_to_element_nodes[element]) {
node_index const node = element_nodes_to_nodes[element_node];
auto const x = nodes_to_x[node].load();
centroid = centroid + x;
}
centroid = centroid * N;
elements_to_centroids[element] = centroid;
};
hpc::for_each(hpc::device_policy(), s.elements, centroid_functor);
for (auto const material : in.materials) {
auto const& domain = in.domains[material];
if (domain) {
domain->mark(centroid_vector, material, &s.material);
}
}
}
void
compute_nodal_materials(input const& in, state& s)
{
auto const nodes_to_node_elements = s.nodes_to_node_elements.cbegin();
auto const node_elements_to_elements = s.node_elements_to_elements.cbegin();
auto const elements_to_materials = s.material.cbegin();
s.nodal_materials.resize(s.nodes.size());
auto const nodes_to_materials = s.nodal_materials.begin();
auto functor = [=] HPC_DEVICE(node_index const node) {
auto const node_elements = nodes_to_node_elements[node];
auto node_materials = material_set::none();
for (auto const node_element : node_elements) {
element_index const element = node_elements_to_elements[node_element];
material_index const element_material = elements_to_materials[element];
node_materials = node_materials | material_set(element_material);
}
nodes_to_materials[node] = node_materials;
};
hpc::for_each(hpc::device_policy(), s.nodes, functor);
for (auto const boundary : in.boundaries) {
auto const& domain = in.domains[boundary];
domain->mark(s.x, boundary, &s.nodal_materials);
}
}
void
collect_node_sets(input const& in, state& s)
{
hpc::counting_range<material_index> const all_materials(in.materials.size() + in.boundaries.size());
s.node_sets.resize(all_materials.size());
assert(s.nodal_materials.size() == s.nodes.size());
auto const nodes_to_materials = s.nodal_materials.cbegin();
for (auto const material : all_materials) {
auto is_in_functor = [=] HPC_DEVICE(node_index const node) -> int {
material_set const materials = nodes_to_materials[node];
return (materials.contains(material_set(material))) ? 1 : 0;
};
collect_set(s.nodes, is_in_functor, s.node_sets[material]);
}
}
void
collect_element_sets(input const& in, state& s)
{
s.element_sets.resize(in.materials.size());
auto const elements_to_material = s.material.cbegin();
for (auto const material : in.materials) {
auto is_in_functor = [=] HPC_DEVICE(element_index const element) -> int {
material_index const element_material = elements_to_material[element];
return (element_material == material) ? 1 : 0;
};
collect_set(s.elements, is_in_functor, s.element_sets[material]);
}
}
std::unique_ptr<domain>
epsilon_around_plane_domain(plane const& p, double eps)
{
auto out = std::make_unique<clipped_domain<all_space>>(all_space{});
out->clip({p.normal, p.origin - eps});
out->clip({-p.normal, -p.origin - eps});
return std::unique_ptr<domain>(std::move(out));
}
std::unique_ptr<domain>
sphere_domain(hpc::position<double> const origin, double const radius)
{
lgr::sphere const s{origin, radius};
auto out = std::make_unique<clipped_domain<lgr::sphere>>(s);
return std::unique_ptr<domain>(std::move(out));
}
std::unique_ptr<domain>
half_space_domain(plane const& p)
{
auto out = std::make_unique<clipped_domain<all_space>>(all_space{});
out->clip(p);
return std::unique_ptr<domain>(std::move(out));
}
std::unique_ptr<domain>
box_domain(hpc::position<double> const lower_left, hpc::position<double> const upper_right)
{
auto out = std::make_unique<clipped_domain<all_space>>(all_space{});
out->clip({hpc::vector3<double>::x_axis(), lower_left(0)});
out->clip({-hpc::vector3<double>::x_axis(), -upper_right(0)});
out->clip({hpc::vector3<double>::y_axis(), lower_left(1)});
out->clip({-hpc::vector3<double>::y_axis(), -upper_right(1)});
out->clip({hpc::vector3<double>::z_axis(), lower_left(2)});
out->clip({-hpc::vector3<double>::z_axis(), -upper_right(2)});
return std::unique_ptr<domain>(std::move(out));
}
} // namespace lgr