--- a/src/dolfinx_access/minmax_p2.cc Mon Jun 29 00:20:42 2026 +0300
+++ b/src/dolfinx_access/minmax_p2.cc Sun Oct 04 13:02:26 2026 -0500
@@ -83,53 +83,41 @@
// printf("%zd\n", d.size());
// }
- // DOF coordinates (may be more than nodes, e.g., intermediate points)
- auto gx = geom.x();
- // auto num_points = gx.size() / 3;
- // Map from cells to DOF indices
- auto geom_dofmap = geom.dofmap();
- assert(geom_dofmap.extent(0) == (size_t)n_cells);
- assert(geom_dofmap.extent(1) == 3);
+ // DOF coordinates and dofmap for the function space itself (not the mesh geometry).
+ // For a scalar P2 space this gives us, per cell, the 6 nodal points (3 corners + 3 edge
+ // midpoints) directly in basix's local ordering, so we no longer need to reconstruct
+ // edge midpoints via f->eval.
auto fp_dofmap = fp->dofmap();
-
- // auto dof_coords = fp->tabulate_dof_coordinates();
+ // transpose=false gives a flat vector laid out as (num_dofs, 3) row-major, i.e.
+ // dof_coords[3*d + component].
+ auto dof_coords = fp->tabulate_dof_coordinates(false);
+ auto v_array = v->array(); // local coefficient vector
CoordValuePair res{{{0.0, 0.0}}, INFINITY};
- // printf("%zd %zd\n", v->array().size(), gx.size());
+ for (auto i = 0; i < n_cells; i++) {
+ auto dofs = fp_dofmap->cell_dofs(i);
+ assert(dofs.size() == 6);
- for (auto i = 0; i < n_cells; i++) {
- // Cell corners
- auto j0 = geom_dofmap(i, 0);
- auto j1 = geom_dofmap(i, 1);
- auto j2 = geom_dofmap(i, 2);
+ // Coordinates of all 6 P2 nodes for this cell, in Fenics' hard-coded 3D layout.
+ std::array<double, 18> x{};
+ // Nodal values of f at those same 6 points.
+ std::array<double, 6> vals{};
- // Construct list of coordinates. Due to f->eval, we construct this is a single list
- // with all the Fenics weirdness of hard-coded 3D.
- std::array<double, 9> x = {gx[3 * j0], gx[3 * j0 + 1], 0,
- gx[3 * j1], gx[3 * j1 + 1], 0,
- gx[3 * j2], gx[3 * j2 + 1], 0};
+ for (std::size_t k = 0; k < 6; k++) {
+ auto d = dofs[k];
+ vals[k] = v_array[d];
+ x[3 * k] = dof_coords[3 * d];
+ x[3 * k + 1] = dof_coords[3 * d + 1];
+ x[3 * k + 2] = 0.0;
+ }
// Find the solution on the simplex.
- auto newres = minmax_dolfinx_p2_cell(f, x, i, max);
+ auto newres = minmax_dolfinx_p2(x, vals, max);
- // printf("%f [%f, %f] [%f %f %f %f %f %f] <%f, %f> <%f, %f> <‸%f, %f>\n", newres.v,
- // newres.x[0], newres.x[1], values[0], values[1], values[2], values[3],
- // values[4], values[5], x[0], x[1], x[3], x[4], x[6], x[7]);
if (newres.v < res.v) {
res = newres;
}
-
- // // Full DOF data
- // auto d = fp_dofmap->cell_dofs(i);
- // // TODO: array() gives just local part of vector
- // assert(d.size() == 6);
- // for (auto k = 0; k < d.size(); k++) {
- // auto j = d[k];
- // // printf("%d %d\n", 3 * j + 1, nodes.size());
- // assert(2 * j + 1 <= nodes.size());
- // auto x = {nodes[3 * j], nodes[3 * j + 1]};
- // }
}
// Negate result if maximising.
@@ -139,6 +127,25 @@
return res;
}
+ void cell_dofs_Function_f64_p2(Function<double> const* f, int32_t cell,
+ std::array<double, 18>& x, std::array<double, 6>& vals) {
+ auto fp = f->function_space();
+ auto fp_dofmap = fp->dofmap();
+ auto dofs = fp_dofmap->cell_dofs(cell);
+ assert(dofs.size() == 6);
+ auto dof_coords = fp->tabulate_dof_coordinates(false);
+ auto v = f->x();
+ auto v_array = v->array(); // local coefficient vector
+
+ for (std::size_t k = 0; k < 6; k++) {
+ auto d = dofs[k];
+ vals[k] = v_array[d];
+ x[3 * k] = dof_coords[3 * d];
+ x[3 * k + 1] = dof_coords[3 * d + 1];
+ x[3 * k + 2] = 0.0;
+ }
+ }
+
CoordValuePair min_Function_f64_p2(Function<double> const* f) {
return minmax_Function_f64_p2(f, false);
}
@@ -154,6 +161,7 @@
m.def("minmax_Function_f64_p2", &minmax_Function_f64_p2);
m.def("eval_Function_f64", &eval_Function_f64);
m.def("cell_Function_f64", &cell_Function_f64);
+ m.def("cell_dofs_Function_f64_p2", &cell_dofs_Function_f64_p2);
// m.def("cell_Mesh_f64", &cell_Mesh_f64);
m.def("cell_FunctionSpace_f64", &cell_FunctionSpace_f64);
nanobind::class_<CoordValuePair>(m, "CoordValuePair")