Mercurial > repos > pointsource_pde / file comparison
comparison: src/dolfinx_access/minmax_p2.cc
src/dolfinx_access/minmax_p2.cc
- changeset 1
- a4137aedcb3a
- child 3
- c3a4f4bb87f7
equal
deleted
inserted
replaced
| 0:7ec1cfe19a24 | 1:a4137aedcb3a |
|---|---|
| 1 /// Minimisation code for second-order polynomial models from Fenics. | |
| 2 #include <array> | |
| 3 #include <cmath> | |
| 4 | |
| 5 #include <basix/finite-element.h> | |
| 6 #include <dolfinx/fem/Function.h> | |
| 7 #include <dolfinx/mesh/cell_types.h> | |
| 8 #include <nanobind/nanobind.h> | |
| 9 #include <nanobind/stl/array.h> // needed for array conversions in nanobind bindings. | |
| 10 | |
| 11 #include "dolfinx_access/minmax_p2.h" | |
| 12 #include "pointsource_pde/src/dolfinx_access.rs.h" | |
| 13 | |
| 14 using namespace dolfinx::fem; | |
| 15 namespace cell = basix::cell; | |
| 16 using namespace basix::element; | |
| 17 using namespace dolfinx_access; | |
| 18 | |
| 19 namespace dolfinx_access { | |
| 20 | |
| 21 // Since Fenics documentation is poor, and it does not seem to provide an easy way to access | |
| 22 // an internal P2 model (that would not involve vomotting your guts out 🤮 dealing with the | |
| 23 // putrid horror that C++ has become (and its documentation written by lawyers and bureacrats, | |
| 24 // which is even worse than that of Fenics), our idea here is to simply assume that f is a | |
| 25 // second-order polynomial on each cell, and then evaluate it on that cell at 6 points, to | |
| 26 // construct a new model based on our existing P2 model code that is also used for minimising | |
| 27 // arbitrary sums. | |
| 28 inline CoordValuePair minmax_Function_f64_p2(Function<double> const* f, bool max) { | |
| 29 auto v = f->x(); | |
| 30 // Mesh | |
| 31 auto fp = f->function_space(); | |
| 32 auto mesh = fp->mesh(); | |
| 33 auto geom = mesh->geometry(); | |
| 34 auto topo = mesh->topology(); | |
| 35 | |
| 36 // Check that we're working with P2 elements | |
| 37 // | |
| 38 // Fenics defaults to GLL warped midpoint placement; see | |
| 39 // https://docs.fenicsproject.org/dolfinx/v0.7.0.post0/python/demos/demo_lagrange_variants.html | |
| 40 // This is not a problem with our current implementation, which works with an | |
| 41 // second-order | |
| 42 // polynomials due to duplication of effort (see above, why), but if we would directly | |
| 43 // take the cell quadrature matrix from Fenics, we should restrict el.lagrange_variant() = | |
| 44 // lagrange_variant::equispaced. | |
| 45 auto el = fp->element()->basix_element(); | |
| 46 // printf("%d %d %d %d\n", el.degree(), el.cell_type(), el.lagrange_variant(), | |
| 47 // el.family()); | |
| 48 if (el.degree() != 2 || el.cell_type() != cell::type::triangle || | |
| 49 /*el.lagrange_variant() != lagrange_variant::equispaced ||*/ el.family() != | |
| 50 family::P) { | |
| 51 throw "Only equispaced Lagrange second-order polynomial elements are supported"; | |
| 52 } | |
| 53 if (cell_num_entities(mesh::CellType::triangle, 0) != 3) { | |
| 54 throw "A triangle should have three vertices"; | |
| 55 } | |
| 56 | |
| 57 // Check that we are dealing with a scalar function | |
| 58 if (fp->value_shape().size() != 0) { | |
| 59 throw "Only scalar functions are supported, obviously."; | |
| 60 } | |
| 61 | |
| 62 // Check that we're in two dimensions | |
| 63 if (topo->dim() != 2 || geom.dim() != 2) { | |
| 64 throw "Only two-dimensional meshes are supported"; | |
| 65 } | |
| 66 | |
| 67 // Check that there are no other types of eleemnts | |
| 68 auto entity_types = topo->entity_types(2); | |
| 69 for (auto& t : entity_types) { | |
| 70 if (t != mesh::CellType::triangle) { | |
| 71 throw "Only triangular meshes are supported"; | |
| 72 } | |
| 73 } | |
| 74 | |
| 75 auto adj = topo->connectivity(2, 0); | |
| 76 auto n_cells = adj->num_nodes(); | |
| 77 | |
| 78 // assert(n_cells == n_triangles); | |
| 79 | |
| 80 // printf("cell count? %d\n", n_cells); | |
| 81 // for (auto i = 0; i < n_cells; i++) { | |
| 82 // // assert(adj->num_links(i) == 3); // Should not need this | |
| 83 // auto d = fp->dofmap()->cell_dofs(i); | |
| 84 // printf("%zd\n", d.size()); | |
| 85 // } | |
| 86 | |
| 87 // DOF coordinates (may be more than nodes, e.g., intermediate points) | |
| 88 auto gx = geom.x(); | |
| 89 // auto num_points = gx.size() / 3; | |
| 90 // Map from cells to DOF indices | |
| 91 auto geom_dofmap = geom.dofmap(); | |
| 92 assert(geom_dofmap.extent(0) == (size_t)n_cells); | |
| 93 assert(geom_dofmap.extent(1) == 3); | |
| 94 auto fp_dofmap = fp->dofmap(); | |
| 95 | |
| 96 // auto dof_coords = fp->tabulate_dof_coordinates(); | |
| 97 | |
| 98 CoordValuePair res{{{0.0, 0.0}}, INFINITY}; | |
| 99 | |
| 100 // printf("%zd %zd\n", v->array().size(), gx.size()); | |
| 101 | |
| 102 for (auto i = 0; i < n_cells; i++) { | |
| 103 // Cell corners | |
| 104 auto j0 = geom_dofmap(i, 0); | |
| 105 auto j1 = geom_dofmap(i, 1); | |
| 106 auto j2 = geom_dofmap(i, 2); | |
| 107 | |
| 108 // Construct list of coordinates. Due to f->eval, we construct this is a single list | |
| 109 // with all the Fenics weirdness of hard-coded 3D. | |
| 110 std::array<double, 9> x = {gx[3 * j0], gx[3 * j0 + 1], 0, | |
| 111 gx[3 * j1], gx[3 * j1 + 1], 0, | |
| 112 gx[3 * j2], gx[3 * j2 + 1], 0}; | |
| 113 | |
| 114 // Find the solution on the simplex. | |
| 115 auto newres = minmax_dolfinx_p2_cell(f, x, i, max); | |
| 116 | |
| 117 // printf("%f [%f, %f] [%f %f %f %f %f %f] <%f, %f> <%f, %f> <‸%f, %f>\n", newres.v, | |
| 118 // newres.x[0], newres.x[1], values[0], values[1], values[2], values[3], | |
| 119 // values[4], values[5], x[0], x[1], x[3], x[4], x[6], x[7]); | |
| 120 if (newres.v < res.v) { | |
| 121 res = newres; | |
| 122 } | |
| 123 | |
| 124 // // Full DOF data | |
| 125 // auto d = fp_dofmap->cell_dofs(i); | |
| 126 // // TODO: array() gives just local part of vector | |
| 127 // assert(d.size() == 6); | |
| 128 // for (auto k = 0; k < d.size(); k++) { | |
| 129 // auto j = d[k]; | |
| 130 // // printf("%d %d\n", 3 * j + 1, nodes.size()); | |
| 131 // assert(2 * j + 1 <= nodes.size()); | |
| 132 // auto x = {nodes[3 * j], nodes[3 * j + 1]}; | |
| 133 // } | |
| 134 } | |
| 135 | |
| 136 // Negate result if maximising. | |
| 137 if (max) { | |
| 138 res.v = -res.v; | |
| 139 } | |
| 140 return res; | |
| 141 } | |
| 142 | |
| 143 CoordValuePair min_Function_f64_p2(Function<double> const* f) { | |
| 144 return minmax_Function_f64_p2(f, false); | |
| 145 } | |
| 146 | |
| 147 CoordValuePair max_Function_f64_p2(Function<double> const* f) { | |
| 148 return minmax_Function_f64_p2(f, true); | |
| 149 } | |
| 150 | |
| 151 // Create a Python module as well | |
| 152 NB_MODULE(_dolfinx_access, m) { | |
| 153 m.def("min_Function_f64_p2", &min_Function_f64_p2); | |
| 154 m.def("max_Function_f64_p2", &max_Function_f64_p2); | |
| 155 m.def("minmax_Function_f64_p2", &minmax_Function_f64_p2); | |
| 156 m.def("eval_Function_f64", &eval_Function_f64); | |
| 157 m.def("cell_Function_f64", &cell_Function_f64); | |
| 158 // m.def("cell_Mesh_f64", &cell_Mesh_f64); | |
| 159 m.def("cell_FunctionSpace_f64", &cell_FunctionSpace_f64); | |
| 160 nanobind::class_<CoordValuePair>(m, "CoordValuePair") | |
| 161 .def_rw("x", &CoordValuePair::x) | |
| 162 .def_rw("v", &CoordValuePair::v); | |
| 163 } | |
| 164 | |
| 165 } // namespace dolfinx_access | |
| 166 | |
| 167 // Forward declaration for manual registration | |
| 168 extern "C" PyObject* PyInit__dolfinx_access(); |
