Mercurial > repos > pointsource_pde / file comparison
comparison: src/dolfinx_access/function.cc
src/dolfinx_access/function.cc
- changeset 3
- c3a4f4bb87f7
- parent 1
- a4137aedcb3a
equal
deleted
inserted
replaced
| 1:a4137aedcb3a | 3:c3a4f4bb87f7 |
|---|---|
| 15 using namespace dolfinx::fem; | 15 using namespace dolfinx::fem; |
| 16 using namespace dolfinx::geometry; | 16 using namespace dolfinx::geometry; |
| 17 using namespace basix::element; | 17 using namespace basix::element; |
| 18 using namespace dolfinx::graph; | 18 using namespace dolfinx::graph; |
| 19 using namespace dolfinx::la; | 19 using namespace dolfinx::la; |
| 20 namespace cell = basix::cell; | |
| 21 | 20 |
| 22 namespace dolfinx_access { | 21 namespace dolfinx_access { |
| 23 extern void drop_Function_f64(Function_f64* f) { | 22 extern void drop_Function_f64(Function_f64* f) { |
| 24 delete f; | 23 delete f; |
| 25 } | 24 } |
| 27 FunctionInfo info_Function_f64(Function_f64 const* f) { | 26 FunctionInfo info_Function_f64(Function_f64 const* f) { |
| 28 auto fp = f->function_space(); | 27 auto fp = f->function_space(); |
| 29 auto el = fp->element()->basix_element(); | 28 auto el = fp->element()->basix_element(); |
| 30 return FunctionInfo{ | 29 return FunctionInfo{ |
| 31 .domain_dim = (uint32_t)fp->mesh()->geometry().dim(), | 30 .domain_dim = (uint32_t)fp->mesh()->geometry().dim(), |
| 32 .codomain_dim = (uint32_t)fp->value_size(), | 31 .codomain_dim = (uint32_t)fp->element()->value_size(), |
| 33 .order = (uint32_t)fp->element()->basix_element().degree(), | 32 .order = (uint32_t)fp->element()->basix_element().degree(), |
| 34 .triangular_mesh = el.cell_type() == cell::type::triangle, | 33 .triangular_mesh = el.cell_type() == basix::cell::type::triangle, |
| 35 }; | 34 }; |
| 36 } | 35 } |
| 37 | 36 |
| 38 void check_compat_Function_f64(Function_f64 const* f, Function_f64 const* g) { | 37 void check_compat_Function_f64(Function_f64 const* f, Function_f64 const* g) { |
| 39 if (f->function_space() != g->function_space()) { | 38 if (f->function_space() != g->function_space()) { |
| 40 throw "Incompatible function spaces"; | 39 throw "Incompatible function spaces"; |
| 41 } | 40 } |
| 42 } | 41 } |
| 43 | 42 |
| 44 std::map<std::weak_ptr<const mesh::Mesh<double>>, std::shared_ptr<BoundingBoxTree<double>>, | 43 std::map<std::weak_ptr<const dolfinx::mesh::Mesh<double>>, |
| 45 std::owner_less<std::weak_ptr<const mesh::Mesh<double>>>> | 44 std::shared_ptr<BoundingBoxTree<double>>, |
| 45 std::owner_less<std::weak_ptr<const dolfinx::mesh::Mesh<double>>>> | |
| 46 bb_tree_cache; | 46 bb_tree_cache; |
| 47 | 47 |
| 48 /// Returns the cell containg x | 48 /// Returns the cell containg x |
| 49 int32_t cell_Mesh_f64(std::shared_ptr<const mesh::Mesh<double>> mesh, | 49 int32_t cell_Mesh_f64(std::shared_ptr<const dolfinx::mesh::Mesh<double>> mesh, |
| 50 const std::array<double, 3>& x) { | 50 const std::array<double, 3>& x) { |
| 51 | 51 |
| 52 std::weak_ptr<const mesh::Mesh<double>> mesh_weak = mesh; | 52 std::weak_ptr<const dolfinx::mesh::Mesh<double>> mesh_weak = mesh; |
| 53 std::shared_ptr<BoundingBoxTree<double>> bb_tree = bb_tree_cache[mesh_weak]; | 53 std::shared_ptr<BoundingBoxTree<double>> bb_tree = bb_tree_cache[mesh_weak]; |
| 54 if (bb_tree == nullptr) { | 54 if (bb_tree == nullptr) { |
| 55 auto topo = mesh->topology(); | 55 auto topo = mesh->topology(); |
| 56 // auto dim = mesh->geometry().dim(); | 56 // auto dim = mesh->geometry().dim(); |
| 57 auto dim = topo->dim(); | 57 auto dim = topo->dim(); |
| 66 int local_cells = index_map->size_local(); | 66 int local_cells = index_map->size_local(); |
| 67 // We don't really do MPI; pointsource_algs is not designed for it. | 67 // We don't really do MPI; pointsource_algs is not designed for it. |
| 68 assert(local_cells == index_map->size_global()); | 68 assert(local_cells == index_map->size_global()); |
| 69 std::vector<std::int32_t> r(local_cells); | 69 std::vector<std::int32_t> r(local_cells); |
| 70 std::iota(r.begin(), r.end(), 0); | 70 std::iota(r.begin(), r.end(), 0); |
| 71 bb_tree = std::make_shared<BoundingBoxTree<double>>(BoundingBoxTree(*mesh, dim, r)); | 71 bb_tree = |
| 72 std::make_shared<BoundingBoxTree<double>>(BoundingBoxTree(*mesh, dim, 0.0, r)); | |
| 72 bb_tree_cache[mesh_weak] = bb_tree; | 73 bb_tree_cache[mesh_weak] = bb_tree; |
| 73 } | 74 } |
| 74 // Find colliding bounding boxes. This is an AdjancencyTree. | 75 // Find colliding bounding boxes. This is an AdjancencyTree. |
| 75 // | 76 // |
| 76 // C++ and Dolfinx are just vomit-inducing 🤮. We need some weird const double span here, | 77 // C++ and Dolfinx are just vomit-inducing 🤮. We need some weird const double span here, |
| 141 f->eval(x, {6, 3}, {{i, i, i, i, i, i}}, std::span(res), {6, 1}); | 142 f->eval(x, {6, 3}, {{i, i, i, i, i, i}}, std::span(res), {6, 1}); |
| 142 return res; | 143 return res; |
| 143 } | 144 } |
| 144 | 145 |
| 145 rust::Slice<const double> data_Function_f64(Function_f64 const* f) { | 146 rust::Slice<const double> data_Function_f64(Function_f64 const* f) { |
| 146 auto span = f->x()->array(); | 147 const std::vector<double>& v = f->x()->array(); |
| 147 return rust::Slice(span.data(), span.size()); | 148 return rust::Slice<const double>(v.data(), v.size()); |
| 148 } | 149 } |
| 149 | 150 |
| 150 rust::Slice<double> data_mut_Function_f64(Function_f64* f) { | 151 rust::Slice<double> data_mut_Function_f64(Function_f64* f) { |
| 151 auto span = f->x()->mutable_array(); | 152 std::vector<double>& v = f->x()->array(); |
| 152 return rust::Slice(span.data(), span.size()); | 153 return rust::Slice<double>(v.data(), v.size()); |
| 153 } | 154 } |
| 154 | 155 |
| 155 Function_f64* similar_Function_f64(Function_f64 const* f) { | 156 Function_f64* similar_Function_f64(Function_f64 const* f) { |
| 156 return new Function_f64(f->function_space()); | 157 return new Function_f64(f->function_space()); |
| 157 } | 158 } |
| 158 | 159 |
| 159 Function_f64* clone_Function_f64(Function_f64 const* f) { | 160 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 // return new Function_f64(f->function_space(), |
| 162 // std::make_shared<Vector<double>>(*f->x())); | |
| 163 auto f_copy = new Function_f64(f->function_space()); | |
| 164 auto fx = f->x(); | |
| 165 auto x_new = f_copy->x(); | |
| 166 const std::vector<double>& a_old = fx->array(); | |
| 167 std::vector<double>& a_new = x_new->array(); | |
| 168 std::copy(a_old.begin(), a_old.end(), a_new.begin()); | |
| 169 x_new->scatter_fwd(); | |
| 170 return f_copy; | |
| 171 } | |
| 172 | |
| 173 void scatter_fwd_Function_f64(Function_f64* f) { | |
| 174 f->x()->scatter_fwd(); | |
| 161 } | 175 } |
| 162 | 176 |
| 163 // PyObject* py_similar_Function_f64(Function_f64 const* f) { | 177 // PyObject* py_similar_Function_f64(Function_f64 const* f) { |
| 164 // return nanobind::cast(Function_f64(f->function_space())).release().ptr(); | 178 // return nanobind::cast(Function_f64(f->function_space())).release().ptr(); |
| 165 // } | 179 // } |
