@@ -17,26 +17,23 @@ def proj_bc(u, f):
1717 return proj (u , f , bcs = DirichletBC (u .function_space (), f , "on_boundary" ))
1818
1919
20- def h1_proj (u , f , bcs = None ):
21- # compute h1 projection of f into u's function
20+ def riesz_proj (u , f , bcs = None ):
21+ # compute the Riesz representative of f in u's function
2222 # space, store the result in u.
2323 v = TestFunction (u .function_space ())
24+ w = TrialFunction (u .function_space ())
2425
25- d = {H2 : grad , H1 : grad , HCurl : curl , HDiv : div , HDivDiv : div }[u .ufl_element ().sobolev_space ]
26- F = (inner (d (u - f ), d (v )) * dx
27- + inner (u - f , v ) * dx )
28- fcp = {"mode" : "vanilla" }
29- solve (F == 0 , u ,
30- bcs = bcs ,
31- solver_parameters = {"snes_type" : "ksponly" ,
32- "ksp_type" : "preonly" ,
33- "pc_type" : "cholesky" },
34- form_compiler_parameters = fcp )
35- return assemble (F (u - f ), form_compiler_parameters = fcp )** 0.5
26+ # Apply the Riesz map <==> projection in the H(d) inner-product norm
27+ d = {H2 : grad , H1 : grad , HCurl : curl , HDiv : div }[u .ufl_element ().sobolev_space ]
28+ a = inner (d (w ), d (v )) * dx + inner (w , v ) * dx
29+ L = a (v , f )
30+ solve (a == L , u , bcs = bcs )
3631
32+ return assemble (a (u - f , u - f ))** 0.5
3733
38- def h1_proj_bc (u , f ):
39- return h1_proj (u , f , bcs = DirichletBC (u .function_space (), f , "on_boundary" ))
34+
35+ def riesz_proj_bc (u , f ):
36+ return riesz_proj (u , f , bcs = DirichletBC (u .function_space (), f , "on_boundary" ))
4037
4138
4239@pytest .fixture (params = ("square" , "cube" ))
@@ -53,7 +50,7 @@ def mesh(request):
5350 ('alfeld' , 1 ),
5451 ('alfeld' , 'd' ),
5552 ('iso(2)' , 2 )])
56- @pytest .mark .parametrize ('op' , (interp , proj , proj_bc , h1_proj , h1_proj_bc ))
53+ @pytest .mark .parametrize ('op' , (interp , proj , proj_bc , riesz_proj , riesz_proj_bc ))
5754def test_projection_scalar_monomial (op , mesh , degree , variant ):
5855 if degree == 'd' :
5956 degree = mesh .geometric_dimension
@@ -112,7 +109,7 @@ def run_convergence(mh, el, deg, convrate, op):
112109 (2 , 'HCT' , 3 , 3 ),
113110 (2 , 'HCT-red' , 3 , 2 ),
114111 ])
115- @pytest .mark .parametrize ('op' , (proj , h1_proj ))
112+ @pytest .mark .parametrize ('op' , (proj , riesz_proj ))
116113def test_scalar_convergence (hierarchy , dim , el , deg , convrate , op ):
117114 if op == proj :
118115 convrate += 1
@@ -132,7 +129,7 @@ def test_scalar_convergence(hierarchy, dim, el, deg, convrate, op):
132129 (3 , 'Guzman-Neilan 1st kind H1' , 1 , 1 ),
133130 (3 , 'Guzman-Neilan H1(div)' , 3 , 2 ),
134131 ])
135- @pytest .mark .parametrize ('op' , (proj , h1_proj ))
132+ @pytest .mark .parametrize ('op' , (proj , riesz_proj ))
136133def test_piola_convergence (hierarchy , dim , el , deg , convrate , op ):
137134 if op == proj :
138135 convrate += 1
0 commit comments