Merge pull request #25 from jpdean/refactor_fluid_solvers

Refactoring

cg_dg_advec_diffusion.py

@@ -150,19 +150,13 @@ def create_mesh(comm, h):
 msh_to_sm_1[sm_1_to_msh] = np.arange(len(sm_1_to_msh))
 
 # Create integration entities for the interface integral
-fdim = tdim - 1
-msh.topology.create_connectivity(tdim, fdim)
-msh.topology.create_connectivity(fdim, tdim)
-facet_imap = msh.topology.index_map(fdim)
-c_to_f = msh.topology.connectivity(tdim, fdim)
-f_to_c = msh.topology.connectivity(fdim, tdim)
+interface_facets = ft.indices[ft.values == bound_ids["interface"]]
 domain_0_cells = ct.indices[ct.values == vol_ids["omega_0"]]
 domain_1_cells = ct.indices[ct.values == vol_ids["omega_1"]]
-interface_facets = ft.indices[ft.values == bound_ids["interface"]]
 interface_entities, msh_to_sm_0, msh_to_sm_1 = \
     compute_interface_integration_entities(
-        interface_facets, domain_0_cells, domain_1_cells, c_to_f,
-        f_to_c, facet_imap, msh_to_sm_0, msh_to_sm_1)
+        msh, interface_facets, domain_0_cells, domain_1_cells,
+        msh_to_sm_0, msh_to_sm_1)
 
 # Compute integration entities for boundary terms
 boundary_entites = compute_integration_domains(
@@ -284,7 +278,7 @@ def create_mesh(comm, h):
 # is enforced weakly by the DG scheme.
 ft_sm_1 = convert_facet_tags(msh, submesh_1, sm_1_to_msh, ft)
 bound_facets_sm_1 = ft_sm_1.indices[ft_sm_1.values == bound_ids["boundary_1"]]
-bound_dofs = fem.locate_dofs_topological(V_1, fdim, bound_facets_sm_1)
+bound_dofs = fem.locate_dofs_topological(V_1, tdim - 1, bound_facets_sm_1)
 u_bc_1 = fem.Function(V_1)
 u_bc_1.interpolate(u_e)
 bc_1 = fem.dirichletbc(u_bc_1, bound_dofs)

hdg_conv_diff.py

@@ -9,12 +9,12 @@
 import numpy as np
 from petsc4py import PETSc
 from dolfinx.cpp.mesh import cell_num_entities
-from utils import norm_L2
+from utils import norm_L2, compute_cell_boundary_integration_entities
 
 
 def u_e(x):
     "Function to represent the exact solution"
-    if type(x) == ufl.SpatialCoordinate:
+    if isinstance(x, ufl.SpatialCoordinate):
         module = ufl
     else:
         module = np
@@ -58,17 +58,13 @@ def boundary(x):
 ubar, vbar = ufl.TrialFunction(Vbar), ufl.TestFunction(Vbar)
 
 # Create integration entities and define integration measures. We want
-# to integrate around each element boundary so we loop over cells and
-# add each facet of the cells as (cell, local facet) pairs
-# TODO Create without Python loop using numpy to improve performance
-all_facets = 0  # Tag
-facet_integration_entities = []
-for cell in range(msh.topology.index_map(tdim).size_local):
-    for local_facet in range(num_cell_facets):
-        facet_integration_entities.extend([cell, local_facet])
+# to integrate around each element boundary, so we call the following
+# convenience function:
+cell_boundary_facets = compute_cell_boundary_integration_entities(msh)
 dx_c = ufl.Measure("dx", domain=msh)
+all_facets = 0  # Tag
 ds_c = ufl.Measure("ds",
-                   subdomain_data=[(all_facets, facet_integration_entities)],
+                   subdomain_data=[(all_facets, cell_boundary_facets)],
                    domain=msh)
 dx_f = ufl.Measure("dx", domain=facet_mesh)