import numpy as np
import numpy.linalg as la
import pyopencl as cl
import pyopencl.clmath
import meshmode.mesh.generation as mgen
from meshmode.discretization import Discretization
from meshmode.discretization.poly_element import PolynomialWarpAndBlendGroupFactory
import matplotlib.pyplot as plt
cl_ctx = cl.create_some_context(interactive=True)
queue = cl.CommandQueue(cl_ctx)
Choose platform: [0] <pyopencl.Platform 'Portable Computing Language' at 0x7f7d74a9f5c0> [1] <pyopencl.Platform 'Intel(R) OpenCL' at 0x195da28>
Set the environment variable PYOPENCL_CTX='0' to avoid being asked again.
order = 10
mesh = mgen.generate_icosphere(r=1, order=order, uniform_refinement_rounds=1)
discr = Discretization(cl_ctx, mesh,
PolynomialWarpAndBlendGroupFactory(order))
from meshmode.discretization.visualization import make_visualizer
vis = make_visualizer(queue, discr, order)
x, y, z = discr.nodes()
val = cl.clmath.sin(3*x.with_queue(queue))
plt.figure(figsize=(10, 7))
vis.show_scalar_in_matplotlib_3d(val, do_show=False)
grp = discr.groups[0]
val_view = grp.view(val.get(queue=queue))
nodes = grp.view(discr.nodes().get(queue=queue))
print(val_view.shape)
print(nodes.shape)
(80, 66) (3, 80, 66)
plt.plot(grp.unit_nodes[0], grp.unit_nodes[1], "o")
plt.gca().set_aspect("equal")
from modepy.matrices import vandermonde
vdm = vandermonde(grp.basis(), grp.unit_nodes)
dr_vdm, ds_vdm = vandermonde(grp.grad_basis(), grp.unit_nodes)
# Use la.solve in real code!
dr = dr_vdm @ la.inv(vdm)
ds = ds_vdm @ la.inv(vdm)
dr_nodes = np.einsum("ij,dkj->dki", dr, nodes)
ds_nodes = np.einsum("ij,dkj->dki", ds, nodes)
normals = np.cross(dr_nodes, ds_nodes, axisa=0, axisb=0, axisc=0)
area_elements = la.norm(normals, axis=0)
area = np.einsum("i, ki->", grp.weights, area_elements)
true_area = 4*np.pi
abs(area - true_area)/true_area
2.5275656016611805e-11