some cpp updates for mesh quality and bug fixes

CMakeLists.txt

@@ -0,0 +1,59 @@
+cmake_minimum_required(VERSION 3.10)
+
+PROJECT(febug_pybind11)
+
+find_package(pybind11 REQUIRED CONFIG HINTS ${PYBIND11_DIR} ${PYBIND11_ROOT}
+  $ENV{PYBIND11_DIR} $ENV{PYBIND11_ROOT})
+
+# Set C++ standard before finding pybind11
+set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+set(CMAKE_CXX_EXTENSIONS OFF)
+
+# Create the binding library
+pybind11_add_module(cpp SHARED
+   febug/cpp_/meshquality.cpp
+   febug/cpp_/febug.cpp)
+
+target_link_libraries(cpp PUBLIC dolfinx)
+target_include_directories(cpp PUBLIC dolfinx)
+
+# Get python include-dirs
+execute_process(
+  COMMAND ${Python3_EXECUTABLE} -c "import dolfinx.wrappers, sys; sys.stdout.write(str(dolfinx.wrappers.get_include_path()))"
+  OUTPUT_VARIABLE DOLFINX_PY_DIR
+  RESULT_VARIABLE DOLFINX_PY_COMMAND_RESULT OUTPUT_STRIP_TRAILING_WHITESPACE)
+
+if (DOLFINX_PY_DIR)
+  message(STATUS "Adding ${DOLFINX_PY_DIR} to include directories")
+  target_include_directories(cpp PRIVATE ${DOLFINX_PY_DIR})
+endif()
+
+# Add strict compiler flags
+include(CheckCXXCompilerFlag)
+CHECK_CXX_COMPILER_FLAG("-Wall -Werror -pedantic" HAVE_PEDANTIC)
+if (HAVE_PEDANTIC)
+  target_compile_options(cpp PRIVATE -Wall;-Werror;-pedantic)
+endif()
+
+
+# Add to CMake search path
+set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_CURRENT_SOURCE_DIR})
+
+# Add DOLFINX libraries and other config
+find_package(DOLFINX REQUIRED)
+if (DOLFINX_FOUND)
+   target_include_directories(cpp PRIVATE ${DOLFINX_INCLUDE_DIR})
+
+# Find petsc4py through python
+execute_process(
+  COMMAND ${PYTHON_EXECUTABLE} -c "import petsc4py; print(petsc4py.get_include())"
+  OUTPUT_VARIABLE PETSC4PY_INCLUDE_DIR
+  RESULT_VARIABLE PETSC4PY_NOT_FOUND
+  ERROR_QUIET
+  OUTPUT_STRIP_TRAILING_WHITESPACE
+  )
+target_include_directories(cpp PRIVATE ${PETSC4PY_INCLUDE_DIR})
+
+endif()
+

examples/demo_ksp_monitors.py

@@ -13,6 +13,11 @@
     [16, 16, 16], dolfinx.mesh.CellType.tetrahedron,
     dolfinx.mesh.GhostMode.shared_facet)
 
+import febug.meshquality
+febug.meshquality.hist_unicode(
+    np.degrees(febug.meshquality.dihedral_angles(mesh)),
+    bins=np.linspace(0, 180, 10), cmd_width=80, title="Dihedral angles")
+
 x = ufl.SpatialCoordinate(mesh)
 f = ufl.sin(ufl.pi*x[0])*ufl.sin(ufl.pi*x[1])*ufl.sin(ufl.pi*x[2])
 
@@ -49,6 +54,6 @@
 
 import febug.monitors
 solver.setMonitor(febug.monitors.monitor_mpl())
-# solver.setMonitor(febug.monitors.monitor_unicode_graph())
-solver.setMonitor(febug.monitors.monitor_text_petsc())
+solver.setMonitor(febug.monitors.monitor_unicode_graph())
+# solver.setMonitor(febug.monitors.monitor_text_petsc())
 solver.solve(b, uh.vector)

examples/demo_meshquality.py

@@ -0,0 +1,69 @@
+import numpy as np
+import dolfinx
+from mpi4py import MPI
+import febug.meshquality
+import matplotlib.pyplot as plt
+
+
+def generate_initial_mesh(algorithm):
+    import gmsh
+    comm = MPI.COMM_WORLD
+    if comm.rank == 0:
+        gmsh.initialize()
+        gmsh.option.setNumber("General.Terminal", 0)
+        gmsh.option.setNumber("Mesh.Algorithm3D", algorithm)
+        model = gmsh.model()
+        model_name = "mesh"
+        model.add(model_name)
+        model.setCurrent(model_name)
+
+        model.occ.addSphere(0, 0, 0, 1)
+
+        model.occ.synchronize()
+        model.mesh.generate(3)
+
+        x = dolfinx.io.extract_gmsh_geometry(model, model_name=model_name)
+        element_types, element_tags, node_tags = model.mesh.getElements(dim=3)
+        assert len(element_types) == 1
+        etype = comm.bcast(element_types[0], root=0)
+        name, dim, order, num_nodes, local_coords, num_first_order_nodes = \
+            model.mesh.getElementProperties(etype)
+        cells = node_tags[0].reshape(-1, num_nodes) - 1
+        num_nodes = comm.bcast(cells.shape[1], root=0)
+    else:
+        num_nodes = comm.bcast(None, root=0)
+        etype = comm.bcast(None, root=0)
+        cells, x = np.empty([0, num_nodes]), np.empty([0, 3])
+
+    mesh = dolfinx.mesh.create_mesh(
+        comm, cells, x,
+        dolfinx.io.ufl_mesh_from_gmsh(etype, x.shape[1]))
+    return mesh
+
+
+for algorithm in (1, 3, 4, 7, 10):
+    counts = []
+    edges = np.linspace(0, 180, 50)
+    for j in range(3):
+        if j == 0:
+            mesh = generate_initial_mesh(algorithm)
+        else:
+            mesh.topology.create_connectivity(1, 0)
+            mesh = dolfinx.mesh.refine(mesh, redistribute=True)
+
+        dangles = np.degrees(febug.meshquality.dihedral_angles(mesh))
+        count, _ = febug.meshquality.histogram_gather(
+            dangles, bins=edges)
+
+        if mesh.comm.rank == 0:
+            plt.bar(edges[:-1], count, width=edges[1:] - edges[:-1], zorder=-j)
+
+    if mesh.comm.rank == 0:
+        plt.gca().set_yscale("log")
+        plt.grid()
+        plt.xlim((edges.min(), edges.max()))
+        plt.xlabel("dihedral angle")
+        plt.ylabel("frequency")
+        plt.title(f"Algorithm {algorithm}")
+        plt.savefig(f"algorithm {algorithm}.png", bbox_inches="tight")
+        plt.clf()

febug/cpp_/febug.cpp

@@ -0,0 +1,21 @@
+#include <iostream>
+#include <pybind11/pybind11.h>
+
+namespace py = pybind11;
+
+namespace febug_wrappers
+{
+void meshquality(py::module& m);
+} // namespace febug_wrappers
+
+PYBIND11_MODULE(cpp, m)
+{
+  // Create module for C++ wrappers
+  m.doc() = "Febug Python interface";
+//  m.attr("__version__") = "0.0.0";
+
+  // Create mpc submodule [mpc]
+  py::module meshquality = m.def_submodule(
+    "meshquality", "mesh quality module");
+  febug_wrappers::meshquality(meshquality);
+}

febug/cpp_/meshquality.cpp

@@ -0,0 +1,91 @@
+#include <dolfinx/common/IndexMap.h>
+#include <dolfinx/common/math.h>
+#include <dolfinx/mesh/Mesh.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/numpy.h>
+#include <pybind11/stl.h>
+namespace py = pybind11;
+
+namespace febug_wrappers
+{
+void meshquality(py::module& m)
+{
+
+  m.def("dihedral_angle", [](dolfinx::mesh::Mesh& mesh)
+  {
+    mesh.topology().create_connectivity(3, 1);
+    mesh.topology().create_connectivity(1, 0);
+    const auto& c2e = mesh.topology().connectivity(3, 1);
+    const auto& e2v = mesh.topology().connectivity(1, 0);
+
+    const xtl::span<const double>& x_g = mesh.geometry().x();
+    std::vector<std::int32_t> entity_list(
+        mesh.topology().index_map(0)->size_local()
+        + mesh.topology().index_map(0)->num_ghosts());
+    std::iota(std::begin(entity_list), std::end(entity_list), 0);
+    const xt::xtensor<std::int32_t, 2> e_to_g = entities_to_geometry(
+        mesh, 0, entity_list, false);
+
+    const std::int32_t num_cells = mesh.topology().index_map(3)->size_local();
+    std::vector<double> dhangles(6*num_cells, 0.0);
+
+    if (mesh.geometry().dim() != 3)
+      throw std::runtime_error("Only supported for 3d meshes of tetrahedra");
+    const std::size_t gdim = 3;
+
+    xt::xtensor_fixed<double, xt::xshape<3>> p0;
+    xt::xtensor_fixed<double, xt::xshape<3>> v1;
+    xt::xtensor_fixed<double, xt::xshape<3>> v2;
+    xt::xtensor_fixed<double, xt::xshape<3>> v3;
+
+    auto copy_x_g = [&x_g, &gdim, &e_to_g](
+        xt::xtensor_fixed<double, xt::xshape<3>>& v,
+        const std::size_t idx)
+    {
+      const std::size_t pos = gdim * e_to_g(idx, 0);
+      for (std::size_t i=0; i < gdim; ++i)
+        v[i] = x_g[pos + i];
+    };
+
+    auto norm2 = [](const xt::xtensor_fixed<double, xt::xshape<3>>& v)
+    { return std::sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]); };
+
+    auto dot = [](const xt::xtensor_fixed<double, xt::xshape<3>>& u,
+                  const xt::xtensor_fixed<double, xt::xshape<3>>& v)
+    { return u[0]*v[0] + u[1]*v[1] + u[2]*v[2]; };
+
+    for (std::int32_t c = 0; c < num_cells; ++c)
+    {
+      const auto& edges = c2e->links(c);
+
+      for (std::int32_t i = 0; i < 6; ++i)
+      {
+        const auto i0 = e2v->links(edges[i])[0];
+        const auto i1 = e2v->links(edges[i])[1];
+        const auto i2 = e2v->links(edges[5-i])[0];
+        const auto i3 = e2v->links(edges[5-i])[1];
+
+        copy_x_g(p0, i0);
+        copy_x_g(v1, i1);
+        copy_x_g(v2, i2);
+        copy_x_g(v3, i3);
+
+        v1 -= p0;
+        v2 -= p0;
+        v3 -= p0;
+
+        v1 /= norm2(v1);
+        v2 /= norm2(v2);
+        v3 /= norm2(v3);
+
+        double cphi = (dot(v2, v3) - dot(v1, v2)*dot(v1, v3));
+        cphi /= (norm2(dolfinx::math::cross(v1, v2)) * norm2(dolfinx::math::cross(v1, v3)));
+
+        double dhangle = std::acos(cphi);
+        dhangles[c*6 + i] = dhangle;
+      }
+    }
+    return py::array_t<double>(dhangles.size(), dhangles.data());
+  });
+}
+}