Skip to content

Commit

Permalink
meshgen: move boxtree code to separate file
Browse files Browse the repository at this point in the history
  • Loading branch information
@alexfikl
alexfikl committed Apr 9, 2024
1 parent ca924ea commit c9a0945
Show file tree
Hide file tree
Showing 3 changed files with 257 additions and 245 deletions.
2 changes: 1 addition & 1 deletion contrib/meshgen11_dealii/meshgen.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,7 +49,7 @@ constexpr double b = 1;
const d::UpdateFlags update_flags =
d::update_quadrature_points | d::update_JxW_values;

void greet() { std::cout << "Hello from meshgen11_dealii." << std::endl; }
void greet() { std::cout << "Hello from mesh generation via deal.II." << std::endl; }

// {{{ utils

Expand Down
279 changes: 35 additions & 244 deletions volumential/meshgen.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,14 +21,17 @@
"""

__doc__ = """
Mesh generation
Mesh generation.
.. autoclass:: MeshGenBase
.. autoclass:: MeshGen1D
:members:
.. autoclass:: MeshGen2D
:members:
.. autoclass:: MeshGen3D
:members:
.. autofunction:: greet
.. autofunction:: make_uniform_cubic_grid
"""

import logging
Expand All @@ -39,265 +42,53 @@
from pytools.obj_array import make_obj_array


logger = logging.getLogger(__name__)

provider = None

# {{{ meshgen Python provider


class MeshGenBase:
"""Base class for Meshgen via BoxTree.
The interface is similar to the Meshgen via Deal.II, except that
the arguments a and b can also be of higher dimensions to allow
added flexibility on choosing the bounding box.
This base class cannot be used directly since the boxtree is not
built in the constructor until dimension_specific_setup() is
provided.
In addition to the common capabilities of the deal.II implementation,
this class also implements native CL getters like get_q_points_dev()
that play well with the other libraries like boxtree.
"""

def __init__(self, degree, nlevels, a=-1, b=1, queue=None):
assert degree > 0
assert nlevels > 0
self.degree = degree
self.quadrature_formula = LegendreGaussQuadrature(degree - 1)
self.nlevels = nlevels

self.bound_a = np.array([a]).flatten()
self.bound_b = np.array([b]).flatten()
assert len(self.bound_a) == len(self.bound_b)
assert np.all(self.bound_a < self.bound_b)

self.dim = len(self.bound_a)
self.root_vertex = self.bound_a
self.root_extent = np.max(self.bound_b - self.bound_a)

if queue is None:
ctx = cl.create_some_context()
self.queue = cl.CommandQueue(ctx)
else:
self.queue = queue

# plug in dimension-specific details
self.dimension_specific_setup()
self.n_q_points = self.degree ** self.dim

self.boxtree = BoxTree()
self.boxtree.generate_uniform_boxtree(
self.queue,
nlevels=self.nlevels,
root_extent=self.root_extent,
root_vertex=self.root_vertex,
)
self.quadrature = QuadratureOnBoxTree(
self.boxtree, self.quadrature_formula
)

def dimension_specific_setup(self):
pass

def get_q_points_dev(self):
return self.quadrature.get_q_points(self.queue)

def get_q_points(self):
q_points_dev = self.get_q_points_dev()
n_all_q_points = len(q_points_dev[0])
q_points = np.zeros((n_all_q_points, self.dim))
for d in range(self.dim):
q_points[:, d] = q_points_dev[d].get(self.queue)
return q_points

def get_q_weights_dev(self):
return self.quadrature.get_q_weights(self.queue)

def get_q_weights(self):
q_weights_dev = self.get_q_weights_dev()
return q_weights_dev.get(self.queue)

def get_cell_measures_dev(self):
return self.quadrature.get_cell_measures(self.queue)

def get_cell_measures(self):
cell_measures_dev = self.get_cell_measures_dev()
return cell_measures_dev.get(self.queue)

def get_cell_centers_dev(self):
return self.quadrature.get_cell_centers(self.queue)

def get_cell_centers(self):
cell_centers_dev = self.get_cell_centers_dev()
n_active_cells = self.n_active_cells()
cell_centers = np.zeros((n_active_cells, self.dim))
for d in range(self.dim):
cell_centers[:, d] = cell_centers_dev[d].get(self.queue)
return cell_centers

def n_cells(self):
"""Note that this value can be larger than the actual number
of cells used in the boxtree. It mainly serves as the bound for
iterators on cells.
"""
return self.boxtree.nboxes

def n_active_cells(self):
return self.boxtree.n_active_boxes

def update_mesh(
self, criteria, top_fraction_of_cells, bottom_fraction_of_cells
):
# TODO
raise NotImplementedError

def print_info(self, logging_func=logger.info):
logging_func("Number of cells: " + str(self.n_cells()))
logging_func("Number of active cells: " + str(self.n_active_cells()))
logging_func("Number of quad points per cell: "
+ str(self.n_q_points))

def generate_gmsh(self, filename):
"""
# TODO
Write the active boxes as a gmsh file.
The file format specifications can be found at:
http://gmsh.info/doc/texinfo/gmsh.html#MSH-ASCII-file-format
"""
raise NotImplementedError


# }}} End meshgen Python provider

try:
logger.info("Trying to find a mesh generator..")
import volumential.meshgen_dealii # noqa: F401

provider = "meshgen_dealii"

except ImportError as e:
logger.debug(repr(e))
logger.warning("Meshgen via Deal.II is not present or unusable.")

try:
logger.info("Trying out BoxTree.TreeInteractiveBuild interface.")
import boxtree.tree_interactive_build # noqa: F401
from modepy import LegendreGaussQuadrature

provider = "meshgen_boxtree"

except ImportError as ee:
logger.debug(repr(ee))
logger.warning("Meshgen via BoxTree is not present or unusable.")
raise RuntimeError("Cannot find a usable Meshgen implementation.")

else:
# {{{ Meshgen via BoxTree
logger.info("Using Meshgen via BoxTree interface.")
from boxtree.tree_interactive_build import BoxTree, QuadratureOnBoxTree

def greet():
return "Hello from Meshgen via BoxTree!"

def make_uniform_cubic_grid(degree, nlevels=1, dim=2, queue=None, **kwargs):
"""Uniform cubic grid in [-1,1]^dim.
This function provides backward compatibility with meshgen_dealii.
"""
if queue is None:
ctx = cl.create_some_context()
queue = cl.CommandQueue(ctx)

# For meshgen_dealii compatibility
if "level" in kwargs:
nlevels = kwargs["level"]

tree = BoxTree()
tree.generate_uniform_boxtree(
queue, nlevels=nlevels, root_extent=2, root_vertex=np.zeros(dim) - 1
)
quad_rule = LegendreGaussQuadrature(degree - 1)
quad = QuadratureOnBoxTree(tree, quad_rule)
q_weights = quad.get_q_weights(queue).get(queue)
q_points_dev = quad.get_q_points(queue)
n_all_q_points = len(q_weights)
q_points = np.zeros((n_all_q_points, dim))
for d in range(dim):
q_points[:, d] = q_points_dev[d].get(queue)

# Adding a placeholder for deprecated point radii
return (q_points, q_weights, None)

class MeshGen1D(MeshGenBase):
"""Meshgen in 1D
"""

def dimension_specific_setup(self):
assert self.dim == 1

class MeshGen2D(MeshGenBase):
"""Meshgen in 2D
"""

def dimension_specific_setup(self):
if self.dim == 1:
# allow passing scalar values of a and b to the constructor
self.dim = 2
self.root_vertex = np.zeros(self.dim) + self.bound_a
else:
assert self.dim == 2

class MeshGen3D(MeshGenBase):
"""Meshgen in 3D
"""

def dimension_specific_setup(self):
if self.dim == 1:
# allow passing scalar values of a and b to the constructor
self.dim = 3
self.root_vertex = np.zeros(self.dim) + self.bound_a
else:
assert self.dim == 3

# }}} End Meshgen via BoxTree

else:
# noexcept on importing meshgen_dealii
logger.info("Using Meshgen via Deal.II interface.")
from volumential.meshgen_dealii import MeshGen2D, MeshGen3D, greet # noqa: F401
from volumential.meshgen_dealii import MeshGen1D, MeshGen2D, MeshGen3D, greet
PROVIDER = "meshgen_dealii"

def make_uniform_cubic_grid(degree, level, dim, queue=None):
from volumential.meshgen_dealii import make_uniform_cubic_grid as _mucg
return _mucg(degree, level, dim)
except ImportError:
try:
from volumential.meshgen_boxtree import (
MeshGen1D, MeshGen2D, MeshGen3D, greet)
PROVIDER = "meshgen_boxtree"
except ImportError:
raise RuntimeError(
"Cannot find a usable mesh generation implementation. "
"Install the deal.II implementation or the boxtree "
"implementation to use mesh generation.")


# FIXME: lowercase variable should be deprecated and remove
provider = PROVIDER
__all__ = (
"provider", "PROVIDER",
"MeshGen1D", "MeshGen2D", "MeshGen3D", "greet",
"make_uniform_cubic_grid",
)

logger = logging.getLogger(__name__)


# {{{ mesh utils

def build_geometry_info(ctx, queue, dim, q_order, mesh,
bbox=None, a=None, b=None):
"""Build tree, traversal and other geo info for FMM computation,
given the box mesh over/encompassing the domain.
bbox=None, a=None, b=None):
"""Build tree, traversal and other geometric information for FMM computation.
The bouding box can be specified in one of two ways:
1. via scalars a, b, dim-homogeneous ([a, b]^dim)
2. via bbox (e.g. np.array([[a1, b1], [a2, b2], [a3, b3]]))
"""

if dim == 1:
if not isinstance(mesh, MeshGen1D):
raise ValueError()

if dim == 2:
if not isinstance(mesh, MeshGen2D):
raise ValueError()

elif dim == 3:
if not isinstance(mesh, MeshGen3D):
raise ValueError()
if dim != mesh.dim:
raise ValueError(
f"Got a mesh of dimension {mesh.dim}, but expected dimension {dim}")

else:
raise ValueError("only supports 1 <= dim <= 3")
if not 1 <= dim <= 3:
raise ValueError(f"Unsupported dimension {dim}: 1 <= dim <= 3")

q_points = mesh.get_q_points()
q_weights = mesh.get_q_weights()
Expand Down
Loading

0 comments on commit c9a0945

Please sign in to comment.