From 0c16d874b53cd8d626b58ac1d0ae86f654e4af20 Mon Sep 17 00:00:00 2001 From: nate-sime Date: Tue, 12 Apr 2022 17:58:35 -0400 Subject: [PATCH] some cpp updates for mesh quality and bug fixes --- CMakeLists.txt | 59 ++++++++++++++++++ examples/demo_ksp_monitors.py | 9 ++- examples/demo_meshquality.py | 69 +++++++++++++++++++++ febug/cpp_/febug.cpp | 21 +++++++ febug/cpp_/meshquality.cpp | 91 ++++++++++++++++++++++++++++ febug/meshquality.py | 110 ++++++++++++++++------------------ setup.py | 81 ++++++++++++++++++++++++- 7 files changed, 376 insertions(+), 64 deletions(-) create mode 100644 CMakeLists.txt create mode 100644 examples/demo_meshquality.py create mode 100644 febug/cpp_/febug.cpp create mode 100644 febug/cpp_/meshquality.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..c707512 --- /dev/null +++ b/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() + diff --git a/examples/demo_ksp_monitors.py b/examples/demo_ksp_monitors.py index 6a88141..6f7ce7d 100644 --- a/examples/demo_ksp_monitors.py +++ b/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) \ No newline at end of file diff --git a/examples/demo_meshquality.py b/examples/demo_meshquality.py new file mode 100644 index 0000000..81faf3a --- /dev/null +++ b/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() diff --git a/febug/cpp_/febug.cpp b/febug/cpp_/febug.cpp new file mode 100644 index 0000000..ce26266 --- /dev/null +++ b/febug/cpp_/febug.cpp @@ -0,0 +1,21 @@ +#include +#include + +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); +} diff --git a/febug/cpp_/meshquality.cpp b/febug/cpp_/meshquality.cpp new file mode 100644 index 0000000..bcd4286 --- /dev/null +++ b/febug/cpp_/meshquality.cpp @@ -0,0 +1,91 @@ +#include +#include +#include +#include +#include +#include +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& x_g = mesh.geometry().x(); + std::vector 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 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 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> p0; + xt::xtensor_fixed> v1; + xt::xtensor_fixed> v2; + xt::xtensor_fixed> v3; + + auto copy_x_g = [&x_g, &gdim, &e_to_g]( + xt::xtensor_fixed>& 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>& v) + { return std::sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]); }; + + auto dot = [](const xt::xtensor_fixed>& u, + const xt::xtensor_fixed>& 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(dhangles.size(), dhangles.data()); + }); +} +} diff --git a/febug/meshquality.py b/febug/meshquality.py index 112d0d4..3ddc107 100644 --- a/febug/meshquality.py +++ b/febug/meshquality.py @@ -1,82 +1,72 @@ +import math import dolfinx.mesh -import numba import numpy as np -import tqdm +from mpi4py import MPI def dihedral_angles(mesh: dolfinx.mesh.Mesh): - mesh.topology.create_connectivity(3, 1) - mesh.topology.create_connectivity(1, 0) - c2e = mesh.topology.connectivity(3, 1) - e2v = mesh.topology.connectivity(1, 0) + import febug.cpp.meshquality + return febug.cpp.meshquality.dihedral_angle(mesh) - x = mesh.geometry.x - dhangles = [] - for c in tqdm.trange(mesh.topology.index_map(3).size_local): - edges = c2e.links(c) +def histogram_gather(data, bins, weights=None, + comm=MPI.COMM_WORLD, rank=0): + counts, edges = np.histogram( + data, bins=bins, weights=weights) - for i in range(6): - i0 = e2v.links(edges[i])[0] - i1 = e2v.links(edges[i])[1] - i2 = e2v.links(edges[5-i])[0] - i3 = e2v.links(edges[5-i])[1] + if comm.rank == rank: + counts_global = np.zeros_like(counts) + else: + counts_global = None - p0 = x[i0] - v1 = x[i1] - p0 - v2 = x[i2] - p0 - v3 = x[i3] - p0 + comm.Reduce([counts, MPI.INT], [counts_global, MPI.INT], op=MPI.SUM, + root=rank) - print(i0, i1, i2, i3) + return counts_global, edges - v1 /= np.linalg.norm(v1) - v2 /= np.linalg.norm(v2) - v3 /= np.linalg.norm(v3) - cphi = (v2.dot(v3) - v1.dot(v2)*v1.dot(v3)) / \ - (np.linalg.norm(np.cross(v1, v2)) * np.linalg.norm(np.cross(v1, v3))) +def hist_unicode(data, bins, weights=None, + cmd_width=80, title=None, comm=MPI.COMM_WORLD): + counts, edges = histogram_gather( + data, bins=bins, weights=weights, comm=comm, rank=0) - dhangle = np.arccos(cphi) - dhangles.append(dhangle) - return dhangles + if comm.rank != 0: + return + if title is not None: + print(title) -def dihedral_angles2(mesh: dolfinx.mesh.Mesh): - mesh.topology.create_connectivity(3, 1) - mesh.topology.create_connectivity(1, 0) - c2e = mesh.topology.connectivity(3, 1) - e2v = mesh.topology.connectivity(1, 0) + barv = "▏▎▍▌▋▊▉█" + barh = "█▇▆▅▄▃▂▁" - x = mesh.geometry.x + xaxis_tick_buffer = 8 - ncells = mesh.topology.index_map(3).size_local - nedges = mesh.topology.index_map(1).size_local + x_intervals = np.logspace( + 0, math.ceil(math.log10(counts.max())), cmd_width + 1) + num_full_blks = np.searchsorted(x_intervals, counts) - dhangles = np.zeros(nedges, dtype=np.double) - half_planes = [] - for c in tqdm.trange(ncells): - edges = c2e.links(c) + # Draw the main bar plots + print(f"{edges[0]:>{xaxis_tick_buffer}.2f} _", flush=True) + for j in range(1, len(edges)): + print(f"{edges[j]:>{xaxis_tick_buffer}.2f} _", end="", flush=True) + print(barv[-1] * num_full_blks[j-1]) - for i in range(6): - i0 = e2v.links(edges[i])[0] - i1 = e2v.links(edges[i])[1] - i2 = e2v.links(edges[5-i])[0] - i3 = e2v.links(edges[5-i])[1] + # Figure out x ticks and their spacing + major_ticks = np.arange(0, math.ceil(math.log10(counts.max()))+1) + tick_pos = np.searchsorted(x_intervals, 10**major_ticks) - p0 = x[i0] - v1 = x[i1] - p0 - v2 = x[i2] - p0 - v3 = x[i3] - p0 + pos = xaxis_tick_buffer + 2 + tick_str = "".join( + "┊" if idx in tick_pos else " " for idx in range(cmd_width+1)) + print(" "*pos + tick_str) - print(i0, i1, i2, i3) + num2ss = str.maketrans("-0123456789", "⁻⁰¹²³⁴⁵⁶⁷⁸⁹") + xticks_str = list(" " * (pos + cmd_width)) - v1 /= np.linalg.norm(v1) - v2 /= np.linalg.norm(v2) - v3 /= np.linalg.norm(v3) - - cphi = (v2.dot(v3) - v1.dot(v2)*v1.dot(v3)) / \ - (np.linalg.norm(np.cross(v1, v2)) * np.linalg.norm(np.cross(v1, v3))) - - dhangle = np.arccos(cphi) - dhangles.append(dhangle) - return dhangles \ No newline at end of file + for tick in major_ticks: + tick_label = f"10{str(tick).translate(num2ss)}" + offset = -1 + local_pos_left = tick_pos[tick] + pos + offset + local_pos_right = tick_pos[tick] + pos + len(tick_label) + offset + xticks_str[local_pos_left:local_pos_right] = tick_label + print("".join(xticks_str)) \ No newline at end of file diff --git a/setup.py b/setup.py index 5d4c80b..3942e58 100644 --- a/setup.py +++ b/setup.py @@ -1,9 +1,86 @@ -from setuptools import setup + +import os +import platform +import re +import subprocess +import sys +from distutils.version import LooseVersion + +from setuptools import Extension, setup +from setuptools.command.build_ext import build_ext + +if sys.version_info < (3, 7): + print("Python 3.7 or higher required, please upgrade.") + sys.exit(1) + +VERSION = "0.0.0" + +REQUIREMENTS = [] + + +class CMakeExtension(Extension): + def __init__(self, name, sourcedir=''): + Extension.__init__(self, name, sources=[]) + self.sourcedir = os.path.abspath(sourcedir) + + +class CMakeBuild(build_ext): + def run(self): + try: + out = subprocess.check_output(['cmake', '--version']) + except OSError: + raise RuntimeError("CMake must be installed to build the" + + "following extensions:" + + ", ".join(e.name for e in self.extensions)) + + if platform.system() == "Windows": + cmake_version = LooseVersion(re.search(r'version\s*([\d.]+)', + out.decode()).group(1)) + if cmake_version < '3.1.0': + raise RuntimeError("CMake >= 3.1.0 is required on Windows") + + for ext in self.extensions: + self.build_extension(ext) + + def build_extension(self, ext): + extdir = os.path.abspath(os.path.dirname( + self.get_ext_fullpath(ext.name))) + cmake_args = ['-DCMAKE_LIBRARY_OUTPUT_DIRECTORY=' + extdir, + '-DPYTHON_EXECUTABLE=' + sys.executable] + + cfg = 'Debug' if self.debug else 'Release' + build_args = ['--config', cfg] + + if platform.system() == "Windows": + cmake_args += ['-DCMAKE_LIBRARY_OUTPUT_DIRECTORY_{}={}'.format( + cfg.upper(), extdir)] + if sys.maxsize > 2** 32: + cmake_args += ['-A', 'x64'] + build_args += ['--', '/m'] + else: + cmake_args += ['-DCMAKE_BUILD_TYPE=' + cfg] + build_args += ['--', '-j3'] + + env = os.environ.copy() + import pybind11 + env['pybind11_DIR'] = pybind11.get_cmake_dir() + env['CXXFLAGS'] = '{} -DVERSION_INFO=\\"{}\\"'.format( + env.get('CXXFLAGS', ''), self.distribution.get_version()) + if not os.path.exists(self.build_temp): + os.makedirs(self.build_temp) + subprocess.check_call(['cmake', ext.sourcedir] + cmake_args, + cwd=self.build_temp, env=env) + subprocess.check_call(['cmake', '--build', '.'] + build_args, + cwd=self.build_temp, env=env) + setup( name='febug', - version='0', + version=VERSION, packages=['febug'], + ext_modules=[CMakeExtension('febug.cpp')], + cmdclass=dict(build_ext=CMakeBuild), + install_requires=REQUIREMENTS, url='', license='MIT', author='Nate Sime',