From 7c2c6d6413bccc2a9a411b7e9b242971435f9185 Mon Sep 17 00:00:00 2001
From: Bogdan Enache <enache2@illinois.edu>
Date: Wed, 6 Dec 2017 19:42:30 -0600
Subject: [PATCH] Cleaned things, placated flake
---
examples/advection/var-velocity.py | 20 +++---
examples/quad/flux_quad.py | 91 ---------------------------
examples/quad/mass_quad.py | 90 ---------------------------
examples/quad/quad.py | 80 ------------------------
examples/quad/stiffness_quad.py | 99 ------------------------------
grudge/discretization.py | 2 +-
grudge/execution.py | 20 +++---
grudge/models/advection.py | 15 +++--
grudge/symbolic/operators.py | 46 ++++----------
grudge/symbolic/primitives.py | 6 +-
10 files changed, 40 insertions(+), 429 deletions(-)
delete mode 100644 examples/quad/flux_quad.py
delete mode 100644 examples/quad/mass_quad.py
delete mode 100644 examples/quad/quad.py
delete mode 100644 examples/quad/stiffness_quad.py
diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py
index 97e793fd..80d656c4 100644
--- a/examples/advection/var-velocity.py
+++ b/examples/advection/var-velocity.py
@@ -30,8 +30,6 @@ logger = logging.getLogger(__name__)
from grudge import sym, bind, DGDiscretizationWithBoundaries
from pytools.obj_array import join_fields
-import meshmode.discretization.poly_element as poly_element
-
def main(write_output=True, order=4):
cl_ctx = cl.create_some_context()
@@ -45,15 +43,12 @@ def main(write_output=True, order=4):
dt_factor = 5
h = 1/10
- final_time = 5
-
-
sym_x = sym.nodes(2)
advec_v = join_fields(-1*sym_x[1], sym_x[0]) / 2
- flux_type = "central"
+ flux_type = "upwind"
source_center = np.array([0.1, 0.1])
source_width = 0.05
@@ -71,11 +66,12 @@ def main(write_output=True, order=4):
from grudge.models.advection import VariableCoefficientAdvectionOperator
- discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, quad_min_degrees={"product": 2*order})
+ discr = DGDiscretizationWithBoundaries(cl_ctx, mesh,
+ order=order, quad_min_degrees={"product": 2*order})
+
op = VariableCoefficientAdvectionOperator(2, advec_v,
- u_analytic(sym.nodes(dim, sym.BTAG_ALL)),quad_tag="product",
+ u_analytic(sym.nodes(dim, sym.BTAG_ALL)), quad_tag="product",
flux_type=flux_type)
- #op = sym.interp("vol", sym.DOFDesc("vol", "product"))(sym.var("u"))
bound_op = bind(discr, op.sym_operator())
@@ -84,7 +80,7 @@ def main(write_output=True, order=4):
def rhs(t, u):
return bound_op(queue, t=t, u=u)
-
+ final_time = 5
dt = dt_factor * h/order**2
nsteps = (final_time // dt) + 1
dt = final_time/nsteps + 1e-15
@@ -101,9 +97,9 @@ def main(write_output=True, order=4):
if isinstance(event, dt_stepper.StateComputed):
step += 1
- print(step)
+ if step % 10 == 0:
+ print(step)
- #print(step, event.t, norm(queue, u=event.state_component[0]))
vis.write_vtk_file("fld-%04d.vtu" % step,
[("u", event.state_component)])
diff --git a/examples/quad/flux_quad.py b/examples/quad/flux_quad.py
deleted file mode 100644
index 3d7f9a32..00000000
--- a/examples/quad/flux_quad.py
+++ /dev/null
@@ -1,91 +0,0 @@
-# Copyright (C) 2007 Andreas Kloeckner
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-
-import numpy as np
-import pyopencl as cl # noqa
-import pyopencl.array # noqa
-import pyopencl.clmath # noqa
-
-import pytest # noqa
-
-from pyopencl.tools import ( # noqa
- pytest_generate_tests_for_pyopencl as pytest_generate_tests)
-
-import logging
-logger = logging.getLogger(__name__)
-
-from grudge import sym, bind, DGDiscretizationWithBoundaries
-
-import numpy.linalg as la
-
-
-
-
-def main(order=1):
- cl_ctx = cl.create_some_context()
- queue = cl.CommandQueue(cl_ctx)
-
- dim = 2
-
-
- from meshmode.mesh.generation import n_gon, generate_regular_rect_mesh
- #mesh = n_gon(3,np.array([0.1,0.2,0.4,0.5,0.6,0.7,0.9]))
- mesh = generate_regular_rect_mesh(a=(-0.5, -0.5), b=(0.5, 0.5),
- n=(3, 3), order=7)
-
-
-
-
- pquad_dd = sym.DOFDesc("vol", "product")
-
- discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, quad_min_degrees={"product": 2 * order} )
-
-
- magic_thing = sym.DOFDesc(sym.FACE_RESTR_ALL, "product")
- #magic_thing = sym.DOFDesc(sym.FACE_RESTR_ALL)
- def fluxy_stuff(u):
- return sym.interp(u.dd, magic_thing)(u.int)
-
-
- ones = bind(discr, sym.nodes(dim)[0] / sym.nodes(dim)[0])(queue, t=0)
- nds = sym.nodes(dim)
- op = sym.FaceMassOperator(magic_thing, "vol")(fluxy_stuff(sym.int_tpair(nds, "product")))
- #op = (sym.int_tpair(nds, "product").int)
- #op = sym.FaceMassOperator()(fluxy_stuff(sym.int_tpair(nds)))
- print(op)
-
- print(sym.pretty(op))
- u = bind(discr, op)(queue, t=0)
-
-
- print(u)
-
-
-
-
- #vis.write_vtk_file("fld-000o.vtu", [ ("u", o) ])
- #vis.write_vtk_file("fld-000u.vtu", [ ("u", u) ])
-
-
-
-
-
-
-
-if __name__ == "__main__":
- main(4)
-
-
diff --git a/examples/quad/mass_quad.py b/examples/quad/mass_quad.py
deleted file mode 100644
index 05c306bd..00000000
--- a/examples/quad/mass_quad.py
+++ /dev/null
@@ -1,90 +0,0 @@
-# Copyright (C) 2007 Andreas Kloeckner
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-
-import numpy as np
-import pyopencl as cl # noqa
-import pyopencl.array # noqa
-import pyopencl.clmath # noqa
-
-import pytest # noqa
-
-from pyopencl.tools import ( # noqa
- pytest_generate_tests_for_pyopencl as pytest_generate_tests)
-
-import logging
-logger = logging.getLogger(__name__)
-
-from grudge import sym, bind, DGDiscretizationWithBoundaries
-
-import numpy.linalg as la
-
-
-
-
-def main(order=1):
- cl_ctx = cl.create_some_context()
- queue = cl.CommandQueue(cl_ctx)
-
- dim = 2
-
-
- from meshmode.mesh.generation import n_gon, generate_regular_rect_mesh
- #mesh = n_gon(3,np.array([0.1,0.2,0.4,0.5,0.6,0.7,0.9]))
- mesh = generate_regular_rect_mesh(a=(-0.5, -0.5), b=(0.5, 0.5),
- n=(3, 3), order=7)
-
-
-
-
- pquad_dd = sym.DOFDesc("vol", "product")
- to_pquad = sym.interp("vol", pquad_dd)
-
- def f(x):
- x = x[0] / x[0]
- #return sym.MassOperator()(x)
- return sym.MassOperator(pquad_dd, "vol")(to_pquad(x))
- #return np.dot(from_pquad_stiffness_t, to_pquad(x))
- #return sym.interp("vol", pquad_dd)(sym.sin(3*x[0]))
-
- discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, quad_min_degrees={"product": 2 * order} )
-
-
-
-
-
-
-
- ones = bind(discr, sym.nodes(dim)[0] / sym.nodes(dim)[0])(queue, t=0)
- u = bind(discr, f(sym.nodes(dim)))(queue, t=0)
-
- print(np.dot(ones.T, u))
-
-
-
-
- #vis.write_vtk_file("fld-000o.vtu", [ ("u", o) ])
- #vis.write_vtk_file("fld-000u.vtu", [ ("u", u) ])
-
-
-
-
-
-
-
-if __name__ == "__main__":
- main(4)
-
-
diff --git a/examples/quad/quad.py b/examples/quad/quad.py
deleted file mode 100644
index 67c41fa2..00000000
--- a/examples/quad/quad.py
+++ /dev/null
@@ -1,80 +0,0 @@
-# Copyright (C) 2007 Andreas Kloeckner
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-
-import numpy as np
-import pyopencl as cl # noqa
-import pyopencl.array # noqa
-import pyopencl.clmath # noqa
-
-import pytest # noqa
-
-from pyopencl.tools import ( # noqa
- pytest_generate_tests_for_pyopencl as pytest_generate_tests)
-
-import logging
-logger = logging.getLogger(__name__)
-
-from grudge import sym, bind, DGDiscretizationWithBoundaries
-
-import numpy.linalg as la
-
-
-
-
-def main(order=1):
- cl_ctx = cl.create_some_context()
- queue = cl.CommandQueue(cl_ctx)
-
- dim = 2
-
-
- from meshmode.mesh.generation import n_gon, generate_regular_rect_mesh
- #mesh = n_gon(3,np.array([0.1,0.2,0.4,0.5,0.6,0.7,0.9]))
- mesh = generate_regular_rect_mesh(a=(-0.5, -0.5), b=(0.5, 0.5),
- n=(3, 3), order=7)
-
-
-
-
- pquad_dd = sym.DOFDesc("vol", "product")
-
- def f(x):
- #return sym.sin(3*x[0])
- return sym.interp("vol", pquad_dd)(sym.sin(3*x[0]))
-
- discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, quad_min_degrees={"product": 2 * order} )
-
- u = bind(discr, f(sym.nodes(dim)))(queue, t=0)
-
-
-
- from meshmode.discretization.visualization import make_visualizer
- #vis = make_visualizer(queue, discr.quad_volume_discr['product'], 2*order)
- #vis = make_visualizer(queue, discr.volume_discr, order)
-
- #vis.write_vtk_file("fld-000.vtu", [ ("u", u) ])
- print(u)
-
-
-
-
-
-
-
-if __name__ == "__main__":
- main(4)
-
-
diff --git a/examples/quad/stiffness_quad.py b/examples/quad/stiffness_quad.py
deleted file mode 100644
index ff334c4e..00000000
--- a/examples/quad/stiffness_quad.py
+++ /dev/null
@@ -1,99 +0,0 @@
-# Copyright (C) 2007 Andreas Kloeckner
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-
-import numpy as np
-import pyopencl as cl # noqa
-import pyopencl.array # noqa
-import pyopencl.clmath # noqa
-
-import pytest # noqa
-
-from pyopencl.tools import ( # noqa
- pytest_generate_tests_for_pyopencl as pytest_generate_tests)
-
-import logging
-logger = logging.getLogger(__name__)
-
-from grudge import sym, bind, DGDiscretizationWithBoundaries
-
-import numpy.linalg as la
-
-
-
-
-def main(order=1):
- cl_ctx = cl.create_some_context()
- queue = cl.CommandQueue(cl_ctx)
-
- dim = 2
-
-
- from meshmode.mesh.generation import n_gon, generate_regular_rect_mesh
- #mesh = n_gon(3,np.array([0.1,0.2,0.4,0.5,0.6,0.7,0.9]))
- mesh = generate_regular_rect_mesh(a=(-0.5, -0.5), b=(0.5, 0.5),
- n=(3, 3), order=7)
-
-
-
-
- pquad_dd = sym.DOFDesc("vol", "product")
- to_pquad = sym.interp("vol", pquad_dd)
- from_pquad_stiffness_t = sym.stiffness_t(dim, pquad_dd, "vol")
- from_pquad_mass = sym.MassOperator()
-
- def lin(x):
- return x[0]
-
- def non_lin(x):
- return x[1] * x[0] * x[0] + 3* x[0] * x[0] - 4
-
- def f(x, fun):
- #return sym.MassOperator()(x)
- #return np.dot(sym.stiffness_t(2), to_pquad(fun(x)))
- return np.dot(from_pquad_stiffness_t, fun(to_pquad(x)))
-
- def reg(x, fun):
- #return sym.MassOperator()(x)
- #return np.dot(sym.stiffness_t(2), to_pquad(fun(x)))
- return np.dot(sym.stiffness_t(2), fun(x))
-
- discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order, quad_min_degrees={"product": 2 * order} )
-
-
- #ones = bind(discr, sym.nodes(dim)[0] / sym.nodes(dim)[0])(queue, t=0)
- nds = sym.nodes(dim)
- #u = bind(discr, to_pquad(nds))(queue, t=0)
- #u = bind(discr, reg(nds, lin))(queue, t=0)
- u = bind(discr, f(nds, lin))(queue, t=0)
-
- print(u)
-
-
-
-
- #vis.write_vtk_file("fld-000o.vtu", [ ("u", o) ])
- #vis.write_vtk_file("fld-000u.vtu", [ ("u", u) ])
-
-
-
-
-
-
-
-if __name__ == "__main__":
- main(4)
-
-
diff --git a/grudge/discretization.py b/grudge/discretization.py
index 2babc3a4..6f5cae13 100644
--- a/grudge/discretization.py
+++ b/grudge/discretization.py
@@ -1,6 +1,6 @@
from __future__ import division, absolute_import
-__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
+__copyright__ = "Copyright (C) 2015-2017 Andreas Kloeckner, Bogdan Enache"
__license__ = """
Permission is hereby granted, free of charge, to any person obtaining a copy
diff --git a/grudge/execution.py b/grudge/execution.py
index b950a961..3fd47c3c 100644
--- a/grudge/execution.py
+++ b/grudge/execution.py
@@ -1,6 +1,6 @@
from __future__ import division, absolute_import
-__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
+__copyright__ = "Copyright (C) 2015-2017 Andreas Kloeckner, Bogdan Enache"
__license__ = """
Permission is hereby granted, free of charge, to any person obtaining a copy
@@ -206,7 +206,6 @@ class ExecutionMapper(mappers.Evaluator,
knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
return lp.tag_inames(knl, dict(k="g.0"))
-
in_discr = self.discrwb.discr_from_dd(op.dd_in)
out_discr = self.discrwb.discr_from_dd(op.dd_out)
@@ -217,24 +216,22 @@ class ExecutionMapper(mappers.Evaluator,
for i in range(len(out_discr.groups)):
in_grp = in_discr.groups[i]
out_grp = out_discr.groups[i]
- # FIXME: This cache thing
+ # FIXME: This cache thing
#cache_key = "elwise_linear", grp, op, field.dtype
#try:
- #matrix = self.bound_op.operator_data_cache[cache_key]
+ #matrix = self.bound_op.operator_data_cache[cache_key]
#except KeyError:
matrix = (
- cl.array.to_device(
- self.queue,
- np.asarray(op.matrix(out_grp, in_grp), dtype=field.dtype))
- .with_queue(None))
+ cl.array.to_device(
+ self.queue,
+ np.asarray(op.matrix(out_grp, in_grp), dtype=field.dtype))
+ .with_queue(None))
- #self.bound_op.operator_data_cache[cache_key] = matrix
+ #self.bound_op.operator_data_cache[cache_key] = matrix
knl()(self.queue, mat=matrix, result=out_grp.view(result),
vec=in_grp.view(field))
-
-
return result
def map_elementwise_max(self, op, field_expr):
@@ -366,7 +363,6 @@ class ExecutionMapper(mappers.Evaluator,
# This should be unified with map_elementwise_linear, which should
# be extended to support batching.
-
# FIXME: Enable
# assert repr_op.dd_in == repr_op.dd_out
assert repr_op.dd_in.domain_tag == repr_op.dd_out.domain_tag
diff --git a/grudge/models/advection.py b/grudge/models/advection.py
index 4d678832..8d18b032 100644
--- a/grudge/models/advection.py
+++ b/grudge/models/advection.py
@@ -3,7 +3,7 @@
from __future__ import division, absolute_import
-__copyright__ = "Copyright (C) 2009 Andreas Kloeckner"
+__copyright__ = "Copyright (C) 2009-2017 Andreas Kloeckner, Bogdan Enache"
__license__ = """
Permission is hereby granted, free of charge, to any person obtaining a copy
@@ -137,13 +137,14 @@ class VariableCoefficientAdvectionOperator(HyperbolicOperator):
self.v = v
self.inflow_u = inflow_u
self.flux_type = flux_type
+
if quad_tag is None:
self.quad_tag = sym.QTAG_NONE
else:
self.quad_tag = quad_tag
def flux(self, u):
-
+
normal = sym.normal(u. dd, self.ambient_dim)
surf_v = sym.interp("vol", u.dd)(self.v)
@@ -173,6 +174,8 @@ class VariableCoefficientAdvectionOperator(HyperbolicOperator):
bc_in = self.inflow_u
all_faces_dd = sym.DOFDesc(sym.FACE_RESTR_ALL, self.quad_tag)
+ boundary_dd = sym.DOFDesc(sym.BTAG_ALL, self.quad_tag)
+
def flux(pair):
return sym.interp(pair.dd, all_faces_dd)(
self.flux(pair))
@@ -180,17 +183,19 @@ class VariableCoefficientAdvectionOperator(HyperbolicOperator):
quad_dd = sym.DOFDesc("vol", self.quad_tag)
if self.quad_tag == sym.QTAG_NONE:
- to_quad = lambda x : x
+ def to_quad(expr):
+ return expr
else:
to_quad = sym.interp("vol", quad_dd)
from_quad_stiffness_t = sym.stiffness_t(self.ambient_dim, quad_dd, "vol")
+ quad_prod = (to_quad(self.v)*to_quad(u))
return sym.InverseMassOperator()(
- np.dot(from_quad_stiffness_t, to_quad(self.v)*to_quad(u))
+ np.dot(from_quad_stiffness_t, quad_prod)
- sym.FaceMassOperator(all_faces_dd, "vol")(
flux(sym.int_tpair(u, self.quad_tag))
- + flux(sym.bv_tpair(sym.DOFDesc(sym.BTAG_ALL, self.quad_tag), u, bc_in))
+ + flux(sym.bv_tpair(boundary_dd, u, bc_in))
# FIXME: Add back support for inflow/outflow tags
#+ flux(sym.bv_tpair(self.inflow_tag, u, bc_in))
diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py
index ff909f7f..0ef8d722 100644
--- a/grudge/symbolic/operators.py
+++ b/grudge/symbolic/operators.py
@@ -2,7 +2,7 @@
from __future__ import division, absolute_import
-__copyright__ = "Copyright (C) 2008 Andreas Kloeckner"
+__copyright__ = "Copyright (C) 2008-2017 Andreas Kloeckner, Bogdan Enache"
__license__ = """
Permission is hereby granted, free of charge, to any person obtaining a copy
@@ -101,7 +101,6 @@ class InterpolationOperator(Operator):
official_dd_in = _sym().as_dofdesc(dd_in)
official_dd_out = _sym().as_dofdesc(dd_out)
if official_dd_in == official_dd_out:
- from warnings import warn
raise ValueError("Interpolating from {} to {}"
" does not do anything.".format(official_dd_in, official_dd_out))
@@ -152,7 +151,7 @@ class DiffOperatorBase(Operator):
dd_in = _sym().DD_VOLUME
if dd_out is None:
dd_out = dd_in.with_qtag(_sym().QTAG_NONE)
- if dd_out == "vol": #FIXME I'm pretty sure doing it this way is wrong
+ if dd_out == "vol": # FIXME I'm pretty sure doing it this way is wrong
dd_out = _sym().DD_VOLUME
if dd_out.uses_quadrature():
raise ValueError("differentiation outputs are not on "
@@ -235,31 +234,21 @@ class RefStiffnessTOperator(RefDiffOperatorBase):
mapper_method = intern("map_ref_stiffness_t")
@staticmethod
- def matrices(out_element_group, in_element_group):
+ def matrices(out_elem_grp, in_elem_grp):
#FIXME: Is this right? Do we expose weird bugs if we delete this?
- if in_element_group == out_element_group:
- assert in_element_group.is_orthogonal_basis()
- mmat = in_element_group.mass_matrix()
- return [dmat.T.dot(mmat.T) for dmat in in_element_group.diff_matrices()]
-
- n_quad_nodes = in_element_group.nunit_nodes
- n_reg_nodes = out_element_group.nunit_nodes
+ if in_elem_grp == out_elem_grp:
+ assert in_elem_grp.is_orthogonal_basis()
+ mmat = in_elem_grp.mass_matrix()
+ return [dmat.T.dot(mmat.T) for dmat in in_elem_grp.diff_matrices()]
from modepy import vandermonde
- vand = vandermonde(out_element_group.basis(), out_element_group.unit_nodes)
- grad_vand = vandermonde(out_element_group.grad_basis(), in_element_group.unit_nodes)
+ vand = vandermonde(out_elem_grp.basis(), out_elem_grp.unit_nodes)
+ grad_vand = vandermonde(out_elem_grp.grad_basis(), in_elem_grp.unit_nodes)
vand_inv_T = np.linalg.inv(vand).T
- grad_basis = out_element_group.grad_basis()
- weights = in_element_group.weights
- nodes = in_element_group.unit_nodes
+ weights = in_elem_grp.weights
- f_res = [np.zeros((n_reg_nodes, n_quad_nodes)) for d in range(out_element_group.dim)]
- for d in range(out_element_group.dim):
- for i in range(n_reg_nodes):
- for j in range(n_quad_nodes):
- f_res[d][i,j] = weights[j] * np.dot(vand_inv_T[i,:], grad_vand[d][j,:])
- return f_res
+ return np.einsum('c,bz,acz->abc', weights, vand_inv_T, grad_vand)
# }}}
@@ -394,25 +383,14 @@ class RefMassOperator(RefMassOperatorBase):
if out_element_group == in_element_group:
return in_element_group.mass_matrix()
- n_quad_nodes = in_element_group.nunit_nodes
- n_reg_nodes = out_element_group.nunit_nodes
-
from modepy import vandermonde
vand = vandermonde(out_element_group.basis(), out_element_group.unit_nodes)
o_vand = vandermonde(out_element_group.basis(), in_element_group.unit_nodes)
vand_inv_T = np.linalg.inv(vand).T
- grad_basis = out_element_group.grad_basis()
weights = in_element_group.weights
- nodes = in_element_group.unit_nodes
-
- f_res = np.zeros((n_reg_nodes, n_quad_nodes))
- for i in range(n_reg_nodes):
- for j in range(n_quad_nodes):
- f_res[i,j] = weights[j] * np.dot(vand_inv_T[i,:], o_vand[j,:])
-
- return f_res
+ return np.einsum('c,bz,acz->abc', weights, vand_inv_T, o_vand)
mapper_method = intern("map_ref_mass")
diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py
index 9234c225..00f7b8d3 100644
--- a/grudge/symbolic/primitives.py
+++ b/grudge/symbolic/primitives.py
@@ -2,7 +2,7 @@
from __future__ import division, absolute_import
-__copyright__ = "Copyright (C) 2008 Andreas Kloeckner"
+__copyright__ = "Copyright (C) 2008-2017 Andreas Kloeckner, Bogdan Enache"
__license__ = """
Permission is hereby granted, free of charge, to any person obtaining a copy
@@ -666,10 +666,6 @@ def int_tpair(expression, qtag=None):
return TracePair(q_dd, i, e)
- #i = cse(_sym().interp("vol", "int_faces")(expression))
- #e = cse(_sym().OppositeInteriorFaceSwap()(i))
- #return TracePair("int_faces", i, e)
-
def bdry_tpair(dd, interior, exterior):
"""
--
GitLab