Mercurial > repos > pointsource_pde / file comparison

comparison: src/dolfinx_access/function.cc

src/dolfinx_access/function.cc

changeset 1
a4137aedcb3a
child 3
c3a4f4bb87f7
equal deleted inserted replaced
0:7ec1cfe19a24 1:a4137aedcb3a
1 #include <array>
2 #include <cmath>
3 #include <span>
4
5 #include <basix/finite-element.h>
6 #include <dolfinx/fem/Function.h>
7 #include <nanobind/nanobind.h>
8
9 #include "dolfinx_access/function.h"
10 // #include "mpi_proto.h"
11 #include "dolfinx/geometry/BoundingBoxTree.h"
12 #include "pointsource_pde/src/dolfinx_access.rs.h"
13 #include "rust/cxx.h"
14
15 using namespace dolfinx::fem;
16 using namespace dolfinx::geometry;
17 using namespace basix::element;
18 using namespace dolfinx::graph;
19 using namespace dolfinx::la;
20 namespace cell = basix::cell;
21
22 namespace dolfinx_access {
23 extern void drop_Function_f64(Function_f64* f) {
24 delete f;
25 }
26
27 FunctionInfo info_Function_f64(Function_f64 const* f) {
28 auto fp = f->function_space();
29 auto el = fp->element()->basix_element();
30 return FunctionInfo{
31 .domain_dim = (uint32_t)fp->mesh()->geometry().dim(),
32 .codomain_dim = (uint32_t)fp->value_size(),
33 .order = (uint32_t)fp->element()->basix_element().degree(),
34 .triangular_mesh = el.cell_type() == cell::type::triangle,
35 };
36 }
37
38 void check_compat_Function_f64(Function_f64 const* f, Function_f64 const* g) {
39 if (f->function_space() != g->function_space()) {
40 throw "Incompatible function spaces";
41 }
42 }
43
44 std::map<std::weak_ptr<const mesh::Mesh<double>>, std::shared_ptr<BoundingBoxTree<double>>,
45 std::owner_less<std::weak_ptr<const mesh::Mesh<double>>>>
46 bb_tree_cache;
47
48 /// Returns the cell containg x
49 int32_t cell_Mesh_f64(std::shared_ptr<const mesh::Mesh<double>> mesh,
50 const std::array<double, 3>& x) {
51
52 std::weak_ptr<const mesh::Mesh<double>> mesh_weak = mesh;
53 std::shared_ptr<BoundingBoxTree<double>> bb_tree = bb_tree_cache[mesh_weak];
54 if (bb_tree == nullptr) {
55 auto topo = mesh->topology();
56 // auto dim = mesh->geometry().dim();
57 auto dim = topo->dim();
58
59 // This fails as it inorrectly calls mesh.topology_mutable(), so need to
60 // copy-paste the implementation below, with the appropriate fix.
61 //*bb_tree = new BoundingBoxTree(*mesh, dim);
62 // And this also doesn't work due to the privacy of `range`.
63 //*bb_tree = new BoundingBoxTree(*mesh, dim, BoundingBoxTree::range(mesh->topology(),
64 // dim
65 auto index_map = topo->index_map(dim);
66 int local_cells = index_map->size_local();
67 // We don't really do MPI; pointsource_algs is not designed for it.
68 assert(local_cells == index_map->size_global());
69 std::vector<std::int32_t> r(local_cells);
70 std::iota(r.begin(), r.end(), 0);
71 bb_tree = std::make_shared<BoundingBoxTree<double>>(BoundingBoxTree(*mesh, dim, r));
72 bb_tree_cache[mesh_weak] = bb_tree;
73 }
74 // Find colliding bounding boxes. This is an AdjancencyTree.
75 //
76 // C++ and Dolfinx are just vomit-inducing 🤮. We need some weird const double span here,
77 // but then just standard x in compute_first_collding_cell below.
78 auto colliding_bbs_adj = compute_collisions(*bb_tree.get(), std::span<const double>(x));
79 // Since we compute the collisions for just one node, the `array` of all links
80 // of the adjancency list, will corresponding exactly the ccells whose bounding boxes
81 // collide with the point.
82 auto colliding_bbs = colliding_bbs_adj.array();
83 // Within the colliding bounding boxes, find the first cell that actually collides.
84 int32_t local_cell = compute_first_colliding_cell(*mesh, colliding_bbs, x, 1e-9);
85 // local_cell may be -1 if it was not found by the local process, so combine
86 // results with MPI_MAX.
87 int32_t global_cell = local_cell;
88 MPI_Allreduce(&local_cell, &global_cell, 1, MPI_INT32_T, MPI_MAX, mesh->comm());
89
90 return global_cell;
91 }
92
93 int32_t cell_FunctionSpace_f64(FunctionSpace<double> const* v,
94 const std::array<double, 3>& x) {
95 return cell_Mesh_f64(v->mesh(), x);
96 }
97
98 int32_t cell_Function_f64(Function_f64 const* f, const std::array<double, 3>& x) {
99 return cell_Mesh_f64(f->function_space()->mesh(), x);
100 }
101
102 std::array<double, 3 * 3> cell_coords_Function_f64_triangle(Function_f64 const* f,
103 int32_t cell) {
104 auto geom = f->function_space()->mesh()->geometry();
105 auto gx = geom.x();
106 auto geom_dofmap = geom.dofmap();
107 assert(geom_dofmap.extent(1) == 3);
108
109 auto j0 = geom_dofmap(cell, 0);
110 auto j1 = geom_dofmap(cell, 1);
111 auto j2 = geom_dofmap(cell, 2);
112
113 // Construct list of coordinates. Due to f->eval, we construct this as a single list
114 // with all the Fenics weirdness of hard-coded 3D.
115 return {gx[3 * j0], gx[3 * j0 + 1], 0, /* */
116 gx[3 * j1], gx[3 * j1 + 1], 0, /* */
117 gx[3 * j2], gx[3 * j2 + 1], 0};
118 }
119
120 // We assume `f` to be real-valued.
121 double eval_Function_f64_cell(Function_f64 const* f, const std::array<double, 3>& x,
122 int32_t cell) {
123 if (cell < 0) {
124 return 0.0;
125 } else {
126 std::array<double, 1> res;
127 f->eval(x, {1, 3}, {{cell}}, std::span(res), {1, 1});
128 return res[0];
129 }
130 }
131
132 double eval_Function_f64(Function_f64 const* f, const std::array<double, 3>& x) {
133 return eval_Function_f64_cell(f, x, cell_Function_f64(f, x));
134 }
135
136 std::array<double, 6> eval_Function_f64_cell_6(Function_f64 const* f,
137 const std::array<double, 3 * 6>& x,
138 int32_t cell) {
139 std::array<double, 6> res;
140 auto i = cell;
141 f->eval(x, {6, 3}, {{i, i, i, i, i, i}}, std::span(res), {6, 1});
142 return res;
143 }
144
145 rust::Slice<const double> data_Function_f64(Function_f64 const* f) {
146 auto span = f->x()->array();
147 return rust::Slice(span.data(), span.size());
148 }
149
150 rust::Slice<double> data_mut_Function_f64(Function_f64* f) {
151 auto span = f->x()->mutable_array();
152 return rust::Slice(span.data(), span.size());
153 }
154
155 Function_f64* similar_Function_f64(Function_f64 const* f) {
156 return new Function_f64(f->function_space());
157 }
158
159 Function_f64* clone_Function_f64(Function_f64 const* f) {
160 return new Function_f64(f->function_space(), std::make_shared<Vector<double>>(*f->x()));
161 }
162
163 // PyObject* py_similar_Function_f64(Function_f64 const* f) {
164 // return nanobind::cast(Function_f64(f->function_space())).release().ptr();
165 // }
166
167 } // namespace dolfinx_access

Mercurial Repository: pointsource_pde

mercurial