@@ -299,14 +299,11 @@ def test_facet_expression(compile_args):
299299 consts = np .array ([], dtype = dtype )
300300 entity_index = np .array ([0 ], dtype = np .intc )
301301 quad_perm = np .array ([0 ], dtype = np .dtype ("uint8" ))
302- tangents = np .array ([coords [1 ] - coords [2 ], coords [2 ] - coords [0 ], coords [0 ] - coords [1 ]])
303- midpoints = np .array (
304- [
305- coords [1 ] + (coords [2 ] - coords [1 ]) / 2 ,
306- coords [0 ] + (coords [2 ] - coords [0 ]) / 2 ,
307- coords [1 ] + (coords [1 ] - coords [0 ]) / 2 ,
308- ]
309- )
302+
303+ edges = basix .topology (basix .CellType .triangle )[1 ]
304+ tangents = np .array ([coords [j ] - coords [i ] for i , j in edges ])
305+ midpoints = np .array ([(coords [i ] + coords [j ]) / 2 for i , j in edges ])
306+
310307 for i , (tangent , midpoint ) in enumerate (zip (tangents , midpoints )):
311308 # normalize tangent
312309 tangent /= np .linalg .norm (tangent )
@@ -329,7 +326,7 @@ def test_facet_expression(compile_args):
329326 assert np .isclose (np .linalg .norm (output ), 1 )
330327
331328 # Check that facet normal is pointing out of the cell
332- assert np .dot (midpoint - coords [ i ] , output ) > 0
329+ assert max ( np .dot (midpoint - c , output ) for c in coords ) > 0
333330
334331
335332def test_facet_geometry_expressions (compile_args ):
@@ -518,11 +515,12 @@ def test_mixed_mesh_expression(compile_args):
518515
519516 # Define exact normals
520517 face_normal = np .array ([0.0 , 0.0 , 1.0 ])
521- e0 = coords [2 ] - coords [1 ]
522- e1 = coords [0 ] - coords [2 ]
523- e2 = coords [1 ] - coords [0 ]
524- edges = np .array ([e0 , e1 , e2 ])
518+ topology = basix .cell .topology (basix .CellType .triangle )
519+ edges = np .array ([coords [j ] - coords [i ] for i , j in topology [1 ]])
525520 edge_normals = np .cross (edges , face_normal )[:, :2 ]
521+ for i , orient in enumerate (basix .cell .facet_orientations (basix .CellType .triangle )):
522+ if orient == 1 :
523+ edge_normals [i ] *= - 1
526524 norms = np .linalg .norm (edge_normals , axis = 1 , keepdims = True )
527525 edge_normals_normalized = edge_normals / norms
528526
@@ -554,7 +552,6 @@ def test_mixed_mesh_expression(compile_args):
554552@pytest .mark .parametrize ("facet_perm" , [0 , 1 ], ids = ["no_perm" , "perm" ])
555553@pytest .mark .parametrize ("local_index" , [0 , 1 , 2 ], ids = ["facet0" , "facet1" , "facet2" ])
556554def test_expression_facet_perm (compile_args , facet_perm , local_index ):
557-
558555 c_el = basix .ufl .element ("Lagrange" , "triangle" , 1 , shape = (2 ,))
559556 mesh = ufl .Mesh (c_el )
560557 expr = ufl .SpatialCoordinate (mesh )
@@ -586,10 +583,10 @@ def test_expression_facet_perm(compile_args, facet_perm, local_index):
586583 entity_index = np .array ([local_index ], dtype = np .intc )
587584
588585 # Perm 0, means that global facet is ordered as
589- # 1 ----2
586+ # 0 ----1
590587 # and quadrature point is not reflected
591588 # Perm 1, means that global facet is ordered as
592- # 2 ----1
589+ # 1 ----0
593590 # and quadrature point should be reflected
594591 quad_perm = np .array ([facet_perm ], dtype = np .uint8 )
595592
@@ -606,8 +603,13 @@ def test_expression_facet_perm(compile_args, facet_perm, local_index):
606603 ffi .NULL ,
607604 )
608605
606+ edges = basix .topology (basix .CellType .triangle )[1 ]
607+
609608 ordered_points = points * (facet_perm == 0 ) + (1 - points ) * (facet_perm == 1 )
610- facet = np . delete ( coords , local_index , axis = 0 )[:, : 2 ]
609+ facet = [ coords [ i ][: 2 ] for i in edges [ local_index ] ]
611610 edge = facet [1 ] - facet [0 ]
611+ print (local_index , facet )
612612 exact_value = facet [0 ] + ordered_points * edge
613- np .testing .assert_allclose (output , exact_value )
613+ print (output )
614+ print (exact_value )
615+ np .testing .assert_allclose (output , exact_value , 1e-10 )
0 commit comments