diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py index 77f57a95619d3879433096b304b351f4670c8160..04b4c635ef61a5a1ad96d36b5684812e5b34d4f9 100644 --- a/grudge/symbolic/mappers/__init__.py +++ b/grudge/symbolic/mappers/__init__.py @@ -825,6 +825,9 @@ class StringifyMapper(pymbolic.mapper.stringifier.StringifyMapper): def map_ones(self, expr, enclosing_prec): return "Ones:" + self._format_dd(expr.dd) + def map_signed_face_ones(self, expr, enclosing_prec): + return "SignedOnes:" + self._format_dd(expr.dd) + # {{{ geometry data def map_node_coordinate_component(self, expr, enclosing_prec): diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py index 2318b70fbca1e4d4d0fa3428d491d53e10baae7e..8948c31718126687cf1b6584eba0ab2a8976a01e 100644 --- a/grudge/symbolic/primitives.py +++ b/grudge/symbolic/primitives.py @@ -533,7 +533,9 @@ def parametrization_derivative(ambient_dim, dim=None, dd=None): dim = ambient_dim if dim == 0: - return MultiVector(np.array([_SignedFaceOnes(dd)])) + from pymbolic.geometric_algebra import get_euclidean_space + return MultiVector(_SignedFaceOnes(dd), + space=get_euclidean_space(ambient_dim)) from pytools import product return product( @@ -621,14 +623,8 @@ def area_element(ambient_dim, dim=None, dd=None): "area_element", cse_scope.DISCRETIZATION) -def mv_normal(dd, ambient_dim, dim=None): - """Exterior unit normal as a :class:`~pymbolic.geometric_algebra.MultiVector`.""" - +def surface_normal(ambient_dim, dim=None, dd=None): dd = as_dofdesc(dd) - - if not dd.is_trace(): - raise ValueError("may only request normals on boundaries") - if dim is None: dim = ambient_dim - 1 @@ -639,23 +635,43 @@ def mv_normal(dd, ambient_dim, dim=None): / area_element(ambient_dim, dim, dd=dd) # Dorst Section 3.7.2 - mv = pder << pder.I.inv() + return cse(pder << pder.I.inv(), + "surface_normal", + cse_scope.DISCRETIZATION) - if dim == 0: - # NOTE: when the mesh is 0D, we do not have a clear notion of - # `exterior normal`, so we just take the tangent of the 1D element, - # interpolate it to the element faces and make it signed - tangent = parametrization_derivative( - ambient_dim, dim=1, dd=DD_VOLUME) - tangent = tangent / sqrt(tangent.norm_squared()) - from grudge.symbolic.operators import interp - interp = interp(DD_VOLUME, dd) - mv = MultiVector(np.array([ - mv.as_scalar() * interp(t) for t in tangent.as_vector() - ])) +def mv_normal(dd, ambient_dim, dim=None): + """Exterior unit normal as a :class:`~pymbolic.geometric_algebra.MultiVector`.""" + dd = as_dofdesc(dd) + if not dd.is_trace(): + raise ValueError("may only request normals on boundaries") + + if dim is None: + dim = ambient_dim - 1 + + if dim == ambient_dim - 1: + return surface_normal(ambient_dim, dim=dim, dd=dd) - return cse(mv, "normal", cse_scope.DISCRETIZATION) + # NOTE: In the case of (d - 2)-dimensional curves, we don't really have + # enough information on the face to decide what an "exterior face normal" + # is (e.g the "normal" to a 1D curve in 3D space is actually a + # "normal plane") + # + # The trick done here is that we take the pseudoscalar in the + # volume, move it to the curve and then contract the face tangent with + # it to get the correct face normal vector. + assert dim == ambient_dim - 2 + + from grudge.symbolic.operators import interp + pder_volm = ( + pseudoscalar(ambient_dim, dim + 1, dd=DD_VOLUME) + .map(interp(DD_VOLUME, dd))) + pder_face = pseudoscalar(ambient_dim, dim, dd=dd) + mv = cse(pder_face << pder_volm, "face_normal", cse_scope.DISCRETIZATION) + + return cse(mv / sqrt(mv.norm_squared()), + "unit_face_normal", + cse_scope.DISCRETIZATION) def normal(dd, ambient_dim, dim=None, quadrature_tag=None): diff --git a/test/test_grudge.py b/test/test_grudge.py index 81c55f041b6bf1c5ca1c065f6071124fbc1e132b..c44909803d15f3990993ab1837e6b260a119af1e 100644 --- a/test/test_grudge.py +++ b/test/test_grudge.py @@ -267,7 +267,112 @@ def test_mass_surface(ctx_factory, name): or eoc_inv.order_estimate() >= (builder.order - 1.0) -@pytest.mark.parametrize("dim", [1, 2, 3]) +@pytest.mark.parametrize("name", ["2-1-ellipse", "spheroid"]) +def test_face_normal_surface(ctx_factory, name): + """Check that face normals are orthogonal to the surface normal""" + + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + + # {{{ cases + + import mesh_data as data + if name == "2-1-ellipse": + builder = data.EllipseMeshBuilder(radius=3.1, aspect_ratio=2.0) + elif name == "spheroid": + builder = data.SpheroidMeshBuilder() + else: + raise ValueError("unknown geometry name: %s" % name) + + # }}} + + # {{{ + + mesh = builder.get_mesh(builder.resolutions[-1], builder.mesh_order) + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=builder.order) + + dv = sym.DD_VOLUME + df = sym.as_dofdesc(sym.FACE_RESTR_ALL) + + surface_normal_sym = sym.surface_normal( + mesh.ambient_dim, dim=mesh.dim, dd=dv).map( + sym.interp(dv, df)).as_vector() + surface_normal = bind(discr, surface_normal_sym)(queue) + + face_normal_sym = sym.normal(df, mesh.ambient_dim, dim=mesh.dim - 1) + face_normal = bind(discr, face_normal_sym)(queue) + + n_dot_n = cl.array.sum( + sum(nv * nf for nv, nf in zip(surface_normal, face_normal)) + ).get(queue) + + t_dot_n = 0 + if mesh.ambient_dim == 3: + # NOTE: there's only one face tangent in 3d + face_tangent_sym = sym.forward_metric_derivative_vector( + mesh.ambient_dim, 0, dd=df) + face_tangent = bind(discr, face_tangent_sym)(queue) + + t_dot_n = cl.array.sum( + sum(tf * nf for tf, nf in zip(face_tangent, face_normal)) + ).get(queue) + + # }}} + + logger.info("dot: %.5e %.5e", n_dot_n, t_dot_n) + assert abs(n_dot_n) < 1.0e-12 + assert abs(t_dot_n) < 1.0e-12 + + +@pytest.mark.parametrize("name", ["2-1-ellipse", "spheroid"]) +def test_face_mass_surface(ctx_factory, name): + """Check that integrating the exterior normals over all faces gives 0.""" + + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + + # {{{ cases + + import mesh_data as data + if name == "2-1-ellipse": + builder = data.EllipseMeshBuilder(radius=3.1, aspect_ratio=2.0) + elif name == "spheroid": + builder = data.SpheroidMeshBuilder() + else: + raise ValueError("unknown geometry name: %s" % name) + + # }}} + + # {{{ convergence + + from pytools.convergence import EOCRecorder + eoc = EOCRecorder() + + for resolution in builder.resolutions: + mesh = builder.get_mesh(resolution, builder.mesh_order) + discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=builder.order) + + # compute surface area + dd = sym.DD_VOLUME + df = sym.as_dofdesc(sym.FACE_RESTR_ALL) + sym_op = sym.NodalSum(dd)(sym.FaceMassOperator(df, dd)( + sym.normal(df, mesh.ambient_dim, dim=mesh.dim - 1))) + + value = la.norm(bind(discr, sym_op)(queue)) + logger.info("integral %.5e", value) + + # compute max element size + h_max = bind(discr, sym.h_max_from_volume(discr.ambient_dim, dd=dd))(queue) + eoc.add_data_point(h_max, value) + + # }}} + + logger.info("\n%s", str(eoc)) + assert eoc.max_error() < 1.0e-12 \ + or eoc.order_estimate() >= (builder.order - 1.0) + + +@pytest.mark.parametrize("dim", [1, 2]) def test_tri_diff_mat(ctx_factory, dim, order=4): """Check differentiation matrix along the coordinate axes on a disk