From e09ac7418b4b82a33d49d16beeac604edcd9cbd8 Mon Sep 17 00:00:00 2001 From: Addison Alvey-Blanco Date: Fri, 11 Oct 2024 10:07:50 -0500 Subject: [PATCH] use all 1d bases and 1d nodes to create quadrature --- meshmode/discretization/poly_element.py | 27 +++++++++++++++++++------ 1 file changed, 21 insertions(+), 6 deletions(-) diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index afc18732..3b533576 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -525,12 +525,27 @@ def basis_obj(self): @memoize_method def quadrature_rule(self): - basis = self._basis.bases[0] - nodes = self.unit_nodes_1d - mass_matrix = mp.mass_matrix(basis, nodes) - weights = np.dot(mass_matrix, - np.ones(len(basis.functions))) - quads = (mp.Quadrature(nodes, weights, exact_to=self.order),)*self.dim + from modepy.tools import reshape_array_for_tensor_product_space as fold + + quads = [] + + if self.dim != 1: + nodes_tp = fold(self.space, self._nodes) + else: + nodes_tp = self._nodes + + for idim, zipped in enumerate(zip(nodes_tp, self._basis.bases)): + nodes, basis = zipped + + # get current dimension's nodes from fastest varying axis + if self.dim != 1: + nodes = np.swapaxes(nodes, 0, idim)[:,*(0,)*(self.dim-1)] + + nodes_1d = nodes.reshape(1, -1) + mass_matrix = mp.mass_matrix(basis, nodes_1d) + weights = np.dot(mass_matrix, np.ones(len(basis.functions))) + + quads.append(mp.Quadrature(nodes_1d, weights, exact_to=self.order)) return mp.TensorProductQuadrature(quads)