@@ -45,21 +45,27 @@ def test_submesh_assemble_mixed_scalar_facet_integral():
4545 mesh .mark_entities (indicator_function , 999 )
4646 subm = Submesh (mesh , dmcommon .CELL_SETS_LABEL , 999 , mesh .topological_dimension ())
4747 subm .init ()
48- V0 = FunctionSpace (mesh , "CG " , 1 )
49- V1 = FunctionSpace (subm , "CG " , 1 )
48+ V0 = FunctionSpace (mesh , "DQ " , 1 , variant = "equispaced" )
49+ V1 = FunctionSpace (subm , "DQ " , 1 , variant = "equispaced" )
5050 V = V0 * V1
5151 u = TrialFunction (V )
5252 v = TestFunction (V )
5353 u0 , u1 = split (u )
5454 v0 , v1 = split (v )
5555 dS0 = Measure ("dS" , domain = mesh )
5656 ds1 = Measure ("ds" , domain = subm )
57- a_p = inner (u1 ('|' ), v0 ('+' )) * dS0
57+ a_p = inner (u1 ('|' ), v0 ('+' )) * dS0 + inner ( u0 ( '+' ), v1 ( '|' )) * ds1 ( 5 )
5858 A_p = assemble (a_p , mat_type = "nest" )
59- print (A_p .M [0 ][0 ].values )
60- print (A_p .M [0 ][1 ].values )
61- print (A_p .M [1 ][0 ].values )
62- print (A_p .M [1 ][1 ].values )
59+ assert np .allclose (A_p .M .sparsity [0 ][0 ].nnz , [1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 ]) # bc nodes
60+ assert np .allclose (A_p .M .sparsity [0 ][1 ].nnz , [4 , 4 , 4 , 4 , 4 , 4 , 4 , 4 ])
61+ assert np .allclose (A_p .M .sparsity [1 ][0 ].nnz , [8 , 8 , 8 , 8 ])
62+ assert np .allclose (A_p .M .sparsity [1 ][1 ].nnz , [1 , 1 , 1 , 1 ]) # bc nodes
63+ M10 = np .array ([[0. , 0. , 0. , 0. , 0. , 0. , 1. / 3. , 1. / 6. ],
64+ [0. , 0. , 0. , 0. , 0. , 0. , 1. / 6. , 1. / 3. ],
65+ [0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ],
66+ [0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ]])
67+ assert np .allclose (A_p .M [0 ][1 ].values , np .transpose (M10 ))
68+ assert np .allclose (A_p .M [1 ][0 ].values , M10 )
6369 b_p = inner (u1 ('|' ), v0 ('+' )) * ds1 (5 )
6470 B_p = assemble (b_p , mat_type = "nest" )
6571 print (B_p .M [0 ][0 ].values )
@@ -72,3 +78,71 @@ def test_submesh_assemble_mixed_scalar_facet_integral():
7278 print (B_m .M [0 ][1 ].values )
7379 print (B_m .M [1 ][0 ].values )
7480 print (B_m .M [1 ][1 ].values )
81+
82+
83+ def test_submesh_assemble_cell_cell_cell_facet_integral ():
84+ # +-------+-------+-------+-------+
85+ # | | | | |
86+ # | | 555 | | mesh
87+ # | | | | |
88+ # +-------+-------+-------+-------+
89+ # +-------+-------+
90+ # | | |
91+ # | | 555 mesh_l
92+ # | | |
93+ # +-------+-------+
94+ # +-------+-------+
95+ # | | |
96+ # 555 | | mesh_r
97+ # | | |
98+ # +-------+-------+
99+ # +-------+
100+ # | |
101+ # 555 | mesh_rl
102+ # | |
103+ # +-------+
104+ mesh = RectangleMesh (4 , 1 , 4. , 1. , quadrilateral = True )
105+ x , y = SpatialCoordinate (mesh )
106+ label_int = 555
107+ label_l = 81100
108+ label_r = 80011
109+ label_rl = 80010
110+ HDivTrace0 = FunctionSpace (mesh , "HDiv Trace" , 0 )
111+ DG0 = FunctionSpace (mesh , "DG" , 0 )
112+ f_int = Function (HDivTrace0 ).interpolate (conditional (And (x > 1.9 , x < 2.1 ), 1 , 0 ))
113+ f_l = Function (DG0 ).interpolate (conditional (x < 2. , 1 , 0 ))
114+ f_r = Function (DG0 ).interpolate (conditional (x > 2. , 1 , 0 ))
115+ f_rl = Function (DG0 ).interpolate (conditional (And (x > 2. , x < 3. ), 1 , 0 ))
116+ mesh = RelabeledMesh (mesh , [f_int , f_l , f_r , f_rl ], [label_int , label_l , label_r , label_rl ])
117+ x , y = SpatialCoordinate (mesh )
118+ mesh_l = Submesh (mesh , dmcommon .CELL_SETS_LABEL , label_l , mesh .topological_dimension ())
119+ mesh_r = Submesh (mesh , dmcommon .CELL_SETS_LABEL , label_r , mesh .topological_dimension ())
120+ mesh_rl = Submesh (mesh_r , dmcommon .CELL_SETS_LABEL , label_rl , mesh .topological_dimension ())
121+ dS = Measure ("dS" , domain = mesh )
122+ ds_l = Measure ("ds" , domain = mesh_l )
123+ ds_r = Measure ("ds" , domain = mesh_r )
124+ ds_rl = Measure ("ds" , domain = mesh_rl )
125+ n_l = FacetNormal (mesh_l )
126+ n_r = FacetNormal (mesh_r )
127+ n_rl = FacetNormal (mesh_rl )
128+ assert assemble (dot (n_rl + n_l , n_rl + n_l ) * ds_rl (label_int )) < 1.e-32
129+ assert assemble (dot (n_rl + n_l , n_rl + n_l ) * ds_r (label_int )) < 1.e-32
130+ assert assemble (dot (n_rl + n_l , n_rl + n_l ) * ds_l (label_int )) < 1.e-32
131+ assert assemble (dot (n_rl + n_l , n_rl + n_l ) * dS (label_int )) < 1.e-32
132+ V_l = FunctionSpace (mesh_l , "DQ" , 1 , variant = 'equispaced' )
133+ V_rl = FunctionSpace (mesh_rl , "DQ" , 1 , variant = 'equispaced' )
134+ V = V_l * V_rl
135+ u_l , u_rl = TrialFunctions (V )
136+ v_l , v_rl = TestFunctions (V )
137+ a = inner (u_rl ('|' ), v_l ('|' )) * ds_l (label_int ) + inner (u_l ('|' ), v_rl ('|' )) * ds_rl (label_int )
138+ A = assemble (a )
139+ print (A .M [0 ][0 ].values )
140+ print (A .M [0 ][1 ].values )
141+ print (A .M [1 ][0 ].values )
142+ print (A .M [1 ][1 ].values )
143+ b = inner (u_rl ('|' ), v_l ('|' )) * dS (label_int ) + inner (u_l ('|' ), v_rl ('|' )) * dS (label_int )
144+ B = assemble (b )
145+ print (B .M [0 ][0 ].values )
146+ print (B .M [0 ][1 ].values )
147+ print (B .M [1 ][0 ].values )
148+ print (B .M [1 ][1 ].values )
0 commit comments