pointsource_pde: comparison src/dolfinx_access/minmax

-:c4f923017349
+:a0d61ac3719e
 //     // assert(adj->num_links(i) == 3); // Should not need this
 //     auto d = fp->dofmap()->cell_dofs(i);
 //     printf("%zd\n", d.size());
 // }
-// DOF coordinates (may be more than nodes, e.g., intermediate points)
+// DOF coordinates and dofmap for the function space itself (not the mesh geometry).
-auto gx = geom.x();
+// For a scalar P2 space this gives us, per cell, the 6 nodal points (3 corners + 3 edge
-// auto num_points = gx.size() / 3;
+// midpoints) directly in basix's local ordering, so we no longer need to reconstruct
-//  Map from cells to DOF indices
+// edge midpoints via f->eval.
-auto geom_dofmap = geom.dofmap();
-assert(geom_dofmap.extent(0) == (size_t)n_cells);
-assert(geom_dofmap.extent(1) == 3);
 auto fp_dofmap = fp->dofmap();
+// transpose=false gives a flat vector laid out as (num_dofs, 3) row-major, i.e.
-// auto dof_coords = fp->tabulate_dof_coordinates();
+// 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++) {
+// Coordinates of all 6 P2 nodes for this cell, in Fenics' hard-coded 3D layout.
-// Cell corners
+std::array<double, 18> x{};
-auto j0 = geom_dofmap(i, 0);
+// Nodal values of f at those same 6 points.
-auto j1 = geom_dofmap(i, 1);
+std::array<double, 6> vals{};
-auto j2 = geom_dofmap(i, 2);
-// Construct list of coordinates. Due to f->eval, we construct this is a single list
+for (std::size_t k = 0; k < 6; k++) {
-// with all the Fenics weirdness of hard-coded 3D.
+auto d = dofs[k];
-std::array<double, 9> x = {gx[3 * j0], gx[3 * j0 + 1], 0,
+vals[k] = v_array[d];
-gx[3 * j1], gx[3 * j1 + 1], 0,
+x[3 * k] = dof_coords[3 * d];
-gx[3 * j2], gx[3 * j2 + 1], 0};
+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.
 if (max) {
 res.v = -res.v;
 }
 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);
 }
 m.def("min_Function_f64_p2", &min_Function_f64_p2);
 m.def("max_Function_f64_p2", &max_Function_f64_p2);
 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")
 .def_rw("x", &CoordValuePair::x)
 .def_rw("v", &CoordValuePair::v);

comparison: src/dolfinx_access/minmax_p2.cc

src/dolfinx_access/minmax_p2.cc

Mercurial > repos > pointsource_pde / file comparison

comparison: src/dolfinx_access/minmax_p2.cc

src/dolfinx_access/minmax_p2.cc