From b419b3155b4f3f122a74f88eb048eb165b0f764c Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Fri, 20 Jun 2014 22:55:32 -0500
Subject: [PATCH] Refactor towards more knowledge in element groups
---
doc/mesh.rst | 11 +-
meshmode/discretization/__init__.py | 248 ++++++++++++++++++-----
meshmode/discretization/poly_element.py | 240 ++++------------------
meshmode/discretization/visualization.py | 8 +-
meshmode/mesh/__init__.py | 60 ++++--
meshmode/mesh/generation.py | 4 +-
meshmode/mesh/io.py | 4 +-
7 files changed, 295 insertions(+), 280 deletions(-)
diff --git a/doc/mesh.rst b/doc/mesh.rst
index a6b9a5b..04f06b4 100644
--- a/doc/mesh.rst
+++ b/doc/mesh.rst
@@ -1,7 +1,7 @@
Mesh management
===============
-.. automodule:: pytential.mesh
+.. automodule:: meshmode.mesh
.. autoclass:: MeshElementGroup
:members:
@@ -14,16 +14,11 @@ Mesh management
Mesh generation
---------------
-.. automodule:: pytential.mesh.generation
+.. automodule:: meshmode.mesh.generation
Mesh input/output
-----------------
-.. automodule:: pytential.mesh.io
-
-Mesh connections (for interpolation)
-------------------------------------
-
-.. automodule:: pytential.mesh.connection
+.. automodule:: meshmode.mesh.io
.. vim: sw=4
diff --git a/meshmode/discretization/__init__.py b/meshmode/discretization/__init__.py
index f2a5ff9..e7a09c2 100644
--- a/meshmode/discretization/__init__.py
+++ b/meshmode/discretization/__init__.py
@@ -22,19 +22,96 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
-
import numpy as np
-from pytools import memoize_method
-
+from pytools import memoize_method, memoize_method_nested
+import loopy as lp
+import pyopencl as cl
__doc__ = """
-
+.. autoclass:: ElementGroupBase
+.. autoclass:: ElementGroupFactory
.. autoclass:: Discretization
"""
+# {{{ element group base
+
+class ElementGroupBase(object):
+ """Container for the :class:`QBXGrid` data corresponding to
+ one :class:`meshmode.mesh.MeshElementGroup`.
+
+ .. attribute :: mesh_el_group
+ .. attribute :: node_nr_base
+ """
+
+ def __init__(self, mesh_el_group, order, node_nr_base):
+ """
+ :arg mesh_el_group: an instance of
+ :class:`meshmode.mesh.MeshElementGroup`
+ """
+ self.mesh_el_group = mesh_el_group
+ self.order = order
+ self.node_nr_base = node_nr_base
+
+ @property
+ def nelements(self):
+ return self.mesh_el_group.nelements
+
+ @property
+ def nunit_nodes(self):
+ return self.unit_nodes.shape[-1]
+
+ @property
+ def nnodes(self):
+ return self.nunit_nodes * self.nelements
+
+ @property
+ def dim(self):
+ return self.mesh_el_group.dim
+
+ def _nodes(self):
+ # Not cached, because the global nodes array is what counts.
+ # This is just used to build that.
+
+ return np.tensordot(
+ self.mesh_el_group.nodes,
+ self._from_mesh_interp_matrix(),
+ (-1, -1))
+
+ def view(self, global_array):
+ return global_array[
+ ..., self.node_nr_base:self.node_nr_base + self.nnodes] \
+ .reshape(
+ global_array.shape[:-1]
+ + (self.nelements, self.nunit_nodes))
+
+# }}}
+
+
+# {{{ group factories
+
+class ElementGroupFactory(object):
+ """
+ .. function:: __call__(mesh_ele_group, node_nr_base)
+ """
+
+
+class OrderBasedGroupFactory(ElementGroupFactory):
+ def __init__(self, order):
+ self.order = order
+
+ def __call__(self, mesh_el_group, node_nr_base):
+ if not isinstance(mesh_el_group, self.mesh_group_class):
+ raise TypeError("only mesh element groups of type '%s' "
+ "are supported" % self.mesh_group_class.__name__)
+
+ return self.group_class(mesh_el_group, self.order, node_nr_base)
+
+# }}}
+
+
class Discretization(object):
- """Abstract interface for discretizations.
+ """An unstructured composite discretization.
.. attribute:: real_dtype
@@ -46,7 +123,11 @@ class Discretization(object):
.. attribute:: ambient_dim
- .. method:: empty(dtype, queue=None, extra_dims=None)
+ .. attribute :: nnodes
+
+ .. attribute :: groups
+
+ .. method:: empty(dtype, queue=None, extra_dims=None)
.. method:: nodes()
@@ -57,51 +138,124 @@ class Discretization(object):
.. method:: quad_weights(queue)
shape: ``(nnodes)``
-
- .. rubric:: Layer potential source discretizations only
-
- .. method:: preprocess_optemplate(name, expr)
-
- .. method:: op_group_features(expr)
-
- Return a characteristic tuple by which operators that can be
- executed together can be grouped.
-
- *expr* is a subclass of
- :class:`pymbolic.primitives.LayerPotentialOperatorBase`.
"""
- @memoize_method
- def _integral_op(self):
- from pytential import sym, bind
- return bind(self, sym.integral(sym.var("integrand")))
-
- def integral(self, queue, x):
- return self._integral_op()(queue, integrand=x)
-
- @memoize_method
- def _norm_op(self, num_components):
- from pytential import sym, bind
- if num_components is not None:
- from pymbolic.primitives import make_sym_vector
- v = make_sym_vector("integrand", num_components)
- integrand = sym.real(np.dot(sym.conj(v), v))
+ def __init__(self, cl_ctx, mesh, group_factory, real_dtype=np.float64):
+ """
+ :arg order: A polynomial-order-like parameter passed unmodified to
+ :attr:`group_class`. See subclasses for more precise definition.
+ """
+
+ self.cl_context = cl_ctx
+
+ self.mesh = mesh
+ self.nnodes = 0
+ self.groups = []
+ for mg in mesh.groups:
+ ng = group_factory(mg, self.nnodes)
+ self.groups.append(ng)
+ self.nnodes += ng.nnodes
+
+ self.real_dtype = np.dtype(real_dtype)
+ self.complex_dtype = {
+ np.float32: np.complex64,
+ np.float64: np.complex128
+ }[self.real_dtype.type]
+
+ @property
+ def dim(self):
+ return self.mesh.dim
+
+ @property
+ def ambient_dim(self):
+ return self.mesh.ambient_dim
+
+ def empty(self, dtype, queue=None, extra_dims=None):
+ if queue is None:
+ first_arg = self.cl_context
else:
- integrand = sym.abs(sym.var("integrand"))**2
-
- return bind(self, sym.integral(integrand))
+ first_arg = queue
+
+ shape = (self.nnodes,)
+ if extra_dims is not None:
+ shape = extra_dims + shape
+
+ return cl.array.empty(first_arg, shape, dtype=dtype)
+
+ def num_reference_derivative(
+ self, queue, ref_axes, vec):
+ @memoize_method_nested
+ def knl():
+ knl = lp.make_kernel(
+ """{[k,i,j]:
+ 0<=k<nelements and
+ 0<=i,j<ndiscr_nodes}""",
+ "result[k,i] = sum(j, diff_mat[i, j] * vec[k, j])",
+ default_offset=lp.auto, name="diff")
+
+ knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
+ return lp.tag_inames(knl, dict(k="g.0"))
+
+ result = self.empty(vec.dtype)
+
+ for grp in self.groups:
+ mat = None
+ for ref_axis in ref_axes:
+ next_mat = grp.diff_matrices()[ref_axis]
+ if mat is None:
+ mat = next_mat
+ else:
+ mat = np.dot(next_mat, mat)
+
+ knl()(queue, diff_mat=mat, result=grp.view(result), vec=grp.view(vec))
+
+ return result
+
+ def quad_weights(self, queue):
+ @memoize_method_nested
+ def knl():
+ knl = lp.make_kernel(
+ "{[k,i]: 0<=k<nelements and 0<=i<ndiscr_nodes}",
+ "result[k,i] = weights[i]",
+ name="quad_weights")
+
+ knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
+ return lp.tag_inames(knl, dict(k="g.0"))
+
+ result = self.empty(self.real_dtype)
+ for grp in self.groups:
+ knl()(queue, result=grp.view(result), weights=grp.weights)
+ return result
- def norm(self, queue, x):
- from pymbolic.geometric_algebra import MultiVector
- if isinstance(x, MultiVector):
- x = x.as_vector(np.object)
-
- num_components = None
- if isinstance(x, np.ndarray):
- num_components, = x.shape
+ @memoize_method
+ def nodes(self):
+ @memoize_method_nested
+ def knl():
+ knl = lp.make_kernel(
+ """{[d,k,i,j]:
+ 0<=d<dims and
+ 0<=k<nelements and
+ 0<=i<ndiscr_nodes and
+ 0<=j<nmesh_nodes}""",
+ """
+ result[d, k, i] = \
+ sum(j, resampling_mat[i, j] * nodes[d, k, j])
+ """,
+ name="nodes")
+
+ knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
+ return lp.tag_inames(knl, dict(k="g.0"))
+
+ result = self.empty(self.real_dtype, extra_dims=(self.ambient_dim,))
+
+ with cl.CommandQueue(self.cl_context) as queue:
+ for grp in self.groups:
+ meg = grp.mesh_el_group
+ knl()(queue,
+ resampling_mat=grp.resampling_matrix(),
+ result=grp.view(result), nodes=meg.nodes)
+
+ return result
- norm_op = self._norm_op(num_components)
- from math import sqrt
- return sqrt(norm_op(queue, integrand=x))
# vim: fdm=marker
diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py
index 114d6ac..43445c6 100644
--- a/meshmode/discretization/poly_element.py
+++ b/meshmode/discretization/poly_element.py
@@ -27,6 +27,7 @@ import numpy as np
#import numpy.linalg as la
from pytools import memoize_method, memoize_method_nested
from meshmode.discretization import Discretization
+from meshmode.mesh import SimplexElementGroup as _MeshSimplexElementGroup
import pyopencl as cl
@@ -46,70 +47,44 @@ __doc__ = """
# loopy sees strides that it doesn't expect and complains.
-# {{{ element group base
+from meshmode.discretization import ElementGroupBase
-class PolynomialElementGroupBase(object):
- """Container for the :class:`QBXGrid` data corresponding to
- one :class:`pytential.mesh.MeshElementGroup`.
- .. attribute :: grid
- .. attribute :: mesh_el_group
- .. attribute :: node_nr_base
- """
-
- def __init__(self, grid, mesh_el_group, order, node_nr_base):
- """
- :arg grid: an instance of :class:`QBXGridBase`
- :arg mesh_el_group: an instance of
- :class:`pytential.mesh.MeshElementGroup`
- """
- self.grid = grid
- self.mesh_el_group = mesh_el_group
- self.order = order
- self.node_nr_base = node_nr_base
-
- @property
- def nelements(self):
- return self.mesh_el_group.nelements
+# {{{ concrete element groups
- @property
- def nunit_nodes(self):
- return self.unit_nodes.shape[-1]
-
- @property
- def nnodes(self):
- return self.nunit_nodes * self.nelements
+class PolynomialSimplexElementGroupBase(ElementGroupBase):
+ def basis(self):
+ return mp.simplex_onb(self.dim, self.order)
@memoize_method
- def _from_mesh_interp_matrix(self):
+ def from_mesh_interp_matrix(self):
meg = self.mesh_el_group
return mp.resampling_matrix(
mp.simplex_onb(meg.dim, meg.order),
self.unit_nodes,
meg.unit_nodes)
- def _nodes(self):
- # Not cached, because the global nodes array is what counts.
- # This is just used to build that.
-
- return np.tensordot(
- self.mesh_el_group.nodes,
- self._from_mesh_interp_matrix(),
- (-1, -1))
-
- def view(self, global_array):
- return global_array[
- ..., self.node_nr_base:self.node_nr_base + self.nnodes] \
- .reshape(
- global_array.shape[:-1]
- + (self.nelements, self.nunit_nodes))
+ @memoize_method
+ def diff_matrices(self):
+ result = mp.differentiation_matrices(
+ self.basis(),
+ mp.grad_simplex_onb(self.dim, self.order),
+ self.unit_nodes)
-# }}}
+ if not isinstance(result, tuple):
+ return (result,)
+ else:
+ return result
+ @memoize_method
+ def resampling_matrix(self):
+ meg = self.mesh_el_group
+ return mp.resampling_matrix(
+ mp.simplex_onb(self.dim, meg.order),
+ self.unit_nodes, meg.unit_nodes)
-# {{{ concrete element groups
-class PolynomialQuadratureElementGroup(PolynomialElementGroupBase):
+class QuadratureSimplexElementGroup(PolynomialSimplexElementGroupBase):
@memoize_method
def _quadrature_rule(self):
dims = self.mesh_el_group.dim
@@ -129,7 +104,7 @@ class PolynomialQuadratureElementGroup(PolynomialElementGroupBase):
return self._quadrature_rule().weights
-class PolynomialWarpAndBlendElementGroup(PolynomialElementGroupBase):
+class PolynomialWarpAndBlendElementGroup(PolynomialSimplexElementGroupBase):
@property
@memoize_method
def unit_nodes(self):
@@ -144,165 +119,34 @@ class PolynomialWarpAndBlendElementGroup(PolynomialElementGroupBase):
# }}}
-# {{{ discretization
-
-class PolynomialElementDiscretizationBase(Discretization):
- """An (unstructured) composite polynomial discretization without
- any specific opinion on how to evaluate layer potentials.
-
- .. attribute :: mesh
- .. attribute :: groups
- .. attribute :: nnodes
-
- .. autoattribute :: nodes
- """
+# {{{ group factories
- def __init__(self, cl_ctx, mesh, order, real_dtype=np.float64):
- """
- :arg order: A polynomial-order-like parameter passed unmodified to
- :attr:`group_class`. See subclasses for more precise definition.
- """
+class ElementGroupFactory(object):
+ pass
- self.cl_context = cl_ctx
- self.mesh = mesh
- self.nnodes = 0
- self.groups = []
- for mg in mesh.groups:
- ng = self.group_class(self, mg, order, self.nnodes)
- self.groups.append(ng)
- self.nnodes += ng.nnodes
-
- self.real_dtype = np.dtype(real_dtype)
- self.complex_dtype = {
- np.float32: np.complex64,
- np.float64: np.complex128
- }[self.real_dtype.type]
-
- @property
- def dim(self):
- return self.mesh.dim
-
- @property
- def ambient_dim(self):
- return self.mesh.ambient_dim
-
- def empty(self, dtype, queue=None, extra_dims=None):
- if queue is None:
- first_arg = self.cl_context
- else:
- first_arg = queue
-
- shape = (self.nnodes,)
- if extra_dims is not None:
- shape = extra_dims + shape
+class OrderBasedGroupFactory(ElementGroupFactory):
+ def __init__(self, order):
+ self.order = order
- return cl.array.empty(first_arg, shape, dtype=dtype)
+ def __call__(self, mesh_el_group, node_nr_base):
+ if not isinstance(mesh_el_group, self.mesh_group_class):
+ raise TypeError("only mesh element groups of type '%s' "
+ "are supported" % self.mesh_group_class.__name__)
- @memoize_method
- def _diff_matrices(self, grp):
- result = mp.differentiation_matrices(
- mp.simplex_onb(self.dim, grp.order),
- mp.grad_simplex_onb(self.dim, grp.order),
- grp.unit_nodes)
+ return self.group_class(mesh_el_group, self.order, node_nr_base)
- if not isinstance(result, tuple):
- return (result,)
- else:
- return result
- def num_reference_derivative(
- self, queue, ref_axes, vec):
- @memoize_method_nested
- def knl():
- knl = lp.make_kernel(
- """{[k,i,j]:
- 0<=k<nelements and
- 0<=i,j<ndiscr_nodes}""",
- "result[k,i] = sum(j, diff_mat[i, j] * vec[k, j])",
- default_offset=lp.auto, name="diff")
-
- knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
- return lp.tag_inames(knl, dict(k="g.0"))
-
- result = self.empty(vec.dtype)
-
- for grp in self.groups:
- mat = None
- for ref_axis in ref_axes:
- next_mat = self._diff_matrices(grp)[ref_axis]
- if mat is None:
- mat = next_mat
- else:
- mat = np.dot(next_mat, mat)
-
- knl()(queue, diff_mat=mat, result=grp.view(result), vec=grp.view(vec))
-
- return result
-
- def quad_weights(self, queue):
- @memoize_method_nested
- def knl():
- knl = lp.make_kernel(
- "{[k,i]: 0<=k<nelements and 0<=i<ndiscr_nodes}",
- "result[k,i] = weights[i]",
- name="quad_weights")
-
- knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
- return lp.tag_inames(knl, dict(k="g.0"))
-
- result = self.empty(self.real_dtype)
- for grp in self.groups:
- knl()(queue, result=grp.view(result), weights=grp.weights)
- return result
+class QuadratureSimplexGroupFactory(OrderBasedGroupFactory):
+ mesh_group_class = _MeshSimplexElementGroup
+ group_class = QuadratureSimplexElementGroup
- @memoize_method
- def _resampling_matrix(self, grp):
- meg = grp.mesh_el_group
- return mp.resampling_matrix(
- mp.simplex_onb(self.dim, meg.order),
- grp.unit_nodes, meg.unit_nodes)
- @memoize_method
- def nodes(self):
- @memoize_method_nested
- def knl():
- knl = lp.make_kernel(
- """{[d,k,i,j]:
- 0<=d<dims and
- 0<=k<nelements and
- 0<=i<ndiscr_nodes and
- 0<=j<nmesh_nodes}""",
- """
- result[d, k, i] = \
- sum(j, resampling_mat[i, j] * nodes[d, k, j])
- """,
- name="nodes")
-
- knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
- return lp.tag_inames(knl, dict(k="g.0"))
-
- result = self.empty(self.real_dtype, extra_dims=(self.ambient_dim,))
-
- with cl.CommandQueue(self.cl_context) as queue:
- for grp in self.groups:
- meg = grp.mesh_el_group
- knl()(queue,
- resampling_mat=self._resampling_matrix(grp),
- result=grp.view(result), nodes=meg.nodes)
-
- return result
+class PolynomialWarpAndBlendGroupFactory(OrderBasedGroupFactory):
+ mesh_group_class = _MeshSimplexElementGroup
+ group_class = PolynomialWarpAndBlendElementGroup
# }}}
-class PolynomialQuadratureElementDiscretization(
- PolynomialElementDiscretizationBase):
- group_class = PolynomialQuadratureElementGroup
-
-
-class PolynomialWarpAndBlendElementDiscretization(
- PolynomialElementDiscretizationBase):
- group_class = PolynomialWarpAndBlendElementGroup
-
# vim: fdm=marker
diff --git a/meshmode/discretization/visualization.py b/meshmode/discretization/visualization.py
index 359acab..9bdf0e0 100644
--- a/meshmode/discretization/visualization.py
+++ b/meshmode/discretization/visualization.py
@@ -205,10 +205,12 @@ class Visualizer(object):
def make_visualizer(queue, discr, vis_order):
+ from meshmode.discretization import Discretization
from meshmode.discretization.poly_element import \
- PolynomialWarpAndBlendElementDiscretization
- vis_discr = PolynomialWarpAndBlendElementDiscretization(
- discr.cl_context, discr.mesh, order=vis_order,
+ PolynomialWarpAndBlendGroupFactory
+ vis_discr = Discretization(
+ discr.cl_context, discr.mesh,
+ PolynomialWarpAndBlendGroupFactory(vis_order),
real_dtype=discr.real_dtype)
from meshmode.discretization.connection import \
make_same_mesh_connection
diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py
index 2ca445b..3ff935a 100644
--- a/meshmode/mesh/__init__.py
+++ b/meshmode/mesh/__init__.py
@@ -70,31 +70,12 @@ class MeshElementGroup(Record):
The nodes are assumed to be mapped versions of *unit_nodes*.
:arg unit_nodes: ``[dim, nunit_nodes]``
The unit nodes of which *nodes* is a mapped
- version. If unspecified, the nodes from
- :func:`modepy.warp_and_blend_nodes` for *dim*
- are assumed. These must be in unit coordinates
- as defined in :mod:`modepy.nodes`.
- :arg dim: only used if *unit_nodes* is None, to get
- the default unit nodes.
+ version.
Do not supply *element_nr_base* and *node_nr_base*, they will be
automatically assigned.
"""
- if unit_nodes is None:
- if dim is None:
- raise TypeError("'dim' must be passed "
- "if 'unit_nodes' is not passed")
-
- unit_nodes = mp.warp_and_blend_nodes(dim, order)
-
- dims = unit_nodes.shape[0]
-
- if vertex_indices.shape[-1] != dims+1:
- raise ValueError("vertex_indices has wrong number of vertices per "
- "element. expected: %d, got: %d" % (dims+1,
- vertex_indices.shape[-1]))
-
Record.__init__(self,
order=order,
vertex_indices=vertex_indices,
@@ -127,6 +108,45 @@ class MeshElementGroup(Record):
def nunit_nodes(self):
return self.unit_nodes.shape[-1]
+
+class SimplexElementGroup(MeshElementGroup):
+ def __init__(self, order, vertex_indices, nodes,
+ element_nr_base=None, node_nr_base=None,
+ unit_nodes=None, dim=None):
+ """
+ :arg order: the mamximum total degree used for interpolation.
+ :arg nodes: ``[ambient_dim, nelements, nunit_nodes]``
+ The nodes are assumed to be mapped versions of *unit_nodes*.
+ :arg unit_nodes: ``[dim, nunit_nodes]``
+ The unit nodes of which *nodes* is a mapped
+ version. If unspecified, the nodes from
+ :func:`modepy.warp_and_blend_nodes` for *dim*
+ are assumed. These must be in unit coordinates
+ as defined in :mod:`modepy.nodes`.
+ :arg dim: only used if *unit_nodes* is None, to get
+ the default unit nodes.
+
+ Do not supply *element_nr_base* and *node_nr_base*, they will be
+ automatically assigned.
+ """
+
+ if unit_nodes is None:
+ if dim is None:
+ raise TypeError("'dim' must be passed "
+ "if 'unit_nodes' is not passed")
+
+ unit_nodes = mp.warp_and_blend_nodes(dim, order)
+
+ dims = unit_nodes.shape[0]
+
+ if vertex_indices.shape[-1] != dims+1:
+ raise ValueError("vertex_indices has wrong number of vertices per "
+ "element. expected: %d, got: %d" % (dims+1,
+ vertex_indices.shape[-1]))
+
+ MeshElementGroup.__init__(self, order, vertex_indices, nodes,
+ element_nr_base, node_nr_base, unit_nodes, dim)
+
# }}}
diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py
index e163f11..cbd4363 100644
--- a/meshmode/mesh/generation.py
+++ b/meshmode/mesh/generation.py
@@ -165,8 +165,8 @@ def make_curve_mesh(curve_f, element_boundaries, order):
t = el_starts[:, np.newaxis] + el_lengths[:, np.newaxis]*nodes_01
nodes = curve_f(t.ravel()).reshape(vertices.shape[0], nelements, -1)
- from meshmode.mesh import Mesh, MeshElementGroup
- egroup = MeshElementGroup(
+ from meshmode.mesh import Mesh, SimplexElementGroup
+ egroup = SimplexElementGroup(
order,
vertex_indices=np.vstack([
np.arange(nelements),
diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py
index c7758b6..160153c 100644
--- a/meshmode/mesh/io.py
+++ b/meshmode/mesh/io.py
@@ -122,7 +122,7 @@ class GmshMeshReceiver(GmshMeshReceiverBase):
# }}}
- from meshmode.mesh import MeshElementGroup, Mesh
+ from meshmode.mesh import SimplexElementGroup, Mesh
for group_el_type, ngroup_elements in el_type_hist.iteritems():
if group_el_type.dimensions != mesh_bulk_dim:
@@ -150,7 +150,7 @@ class GmshMeshReceiver(GmshMeshReceiverBase):
unit_nodes = (np.array(group_el_type.lexicographic_node_tuples(),
dtype=np.float64).T/group_el_type.order)*2 - 1
- groups.append(MeshElementGroup(
+ groups.append(SimplexElementGroup(
group_el_type.order,
vertex_indices,
nodes,
--
GitLab