About the application of boundary conditions at overlapping degrees of freedom

General
,
1

Following from this tutorial, it looks like the way to assign boundary conditions in scenarios where the boundary conditions overlap at the corners is to let the degrees of freedom overlap, and that this is also true for different boundary conditions. I wanted to make sure of whether or not that is always the case, or if I have misinterpreted the code.

Right now, I am currently struggling with a Helmholtz decomposition problem, and I implemented a conditional for whether or not to have the degrees of freedom to interpolate our boundary condition to should overlap or not. To help explain my code, it solves for \psi_I and \chi_I to obtain \textbf{V}_h = \textbf{V} - \hat{\textbf{k}}\times \nabla \psi_I + \nabla \chi_I (in case there is confusion on \psi_I, I believe \hat{\textbf{k}}\times \nabla \psi_I = \nabla \times \textbf{A}). Here is my code (the boundary conditions are handled in the class):

from mpi4py import MPI
import dolfinx as df
import ufl
from petsc4py import PETSc
import numpy as np
from dolfinx.fem.petsc import LinearProblem

class bc_setter:
    def __init__(self):
        self.x = ufl.SpatialCoordinate(mesh)
        self.axis = 0 # initial axis of integration
        #self.t = ufl.perp(ufl.FacetNormal(mesh))
        # initial x_c
        self.x_c = df.fem.Constant(mesh, PETSc.ScalarType(0.0))
        
        # define all boundary facets for integration
        self.boundaries = [
            (1, lambda x: np.isclose(x[1], -Ly)), # Southern boundary
            (2, lambda x: np.isclose(x[0], Lx)),  # Eastern boundary
            (3, lambda x: np.isclose(x[1], Ly)),  # Northern boundary
            (4, lambda x: np.isclose(x[0], -Lx)), # Western boundary
                ]

        facet_indices, facet_markers = [], []
        for marker, locator in self.boundaries:
            facets = df.mesh.locate_entities(mesh, fdim, locator)
            facet_indices.append(facets)
            facet_markers.append(np.full_like(facets, marker))

        facet_indices = np.hstack(facet_indices).astype(np.int32)
        facet_markers = np.hstack(facet_markers).astype(np.int32)
        sorted_facets = np.argsort(facet_indices)
        self.facet_tag = df.mesh.meshtags(
            mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets]
            )

        self.ds = ufl.Measure("ds", domain=mesh, subdomain_data=self.facet_tag) # boundary integration variable
        
        # initialize integral as self variable
        self.compiled_form = df.fem.form(-ufl.dot(n,V_h)
                                            * ufl.conditional(self.x[0]<self.x_c,1,0) * ufl.ds)
   
    # Function called by interpolate() to evaluate integral
    def cumul_int_expr(self,x):
        out_arr = np.zeros(x.shape[1])
        for i in range(x.shape[1]):
            self.x_c.value = x[self.axis][i]
            cumul_int_x_eval = df.fem.assemble_scalar(self.compiled_form)
            cixe_refined = mesh.comm.allreduce(cumul_int_x_eval, op = MPI.SUM)
            out_arr[i] = cixe_refined
        return out_arr
        
    def set_psi_bcs(self):
        psi_h_bcs = []
        bf = []
        for boundary in self.boundaries:
            boundary_marker = boundary[0]
             # create constant from previous psih_s endpoint value
            if boundary_marker == 1:
                psi_h_s0 = 0
            else:
                psi_h_s0 = psi_h_s.x.array[df.fem.locate_dofs_topological(S,fdim,self.facet_tag.find(boundary_marker-1))][np.argsort(S.tabulate_dof_coordinates()[df.fem.locate_dofs_topological(S,fdim,self.facet_tag.find(boundary_marker-1))][:,self.axis])][endpoint_index]
        
            if boundary_marker%2 == 0:
                self.axis = 1 # y-axis
                s = self.x[1] 
            else:
                self.axis = 0 # x-axis
                s = self.x[0]

            if boundary_marker <= 2: # positive integration direction
                endpoint_index = -1 
                self.compiled_form = df.fem.form(-ufl.dot(n,V_h)
                                            * ufl.conditional(s<self.x_c,1,0) * self.ds(boundary_marker))         
            else: # negative integration direction
                endpoint_index = 0
                self.compiled_form = df.fem.form(-ufl.dot(n,V_h)
                                            * ufl.conditional(s<self.x_c,0,1) * self.ds(boundary_marker))         

            # ---- psi_h(s) = psi_h(s_0) - int_{s_0}^s dot(V_h,n) ds ----
            psi_h_s = df.fem.Function(S)
            psi_h_s.interpolate(self.cumul_int_expr)
            psi_h_s.x.array[:] = psi_h_s0 + psi_h_s.x.array[:] 

            compiled_form_n = df.fem.form(ufl.dot(ufl.FacetNormal(mesh),V_h)* self.ds(boundary_marker))
            bf.append(mesh.comm.allreduce(df.fem.assemble_scalar(compiled_form_n), op = MPI.SUM)) # append boundary flux
            # interpolate psi_h(s) values onto boundary dofs 
            if overlapping_facets:
                interpolation_facets = self.facet_tag.find(boundary_marker)
                interpolation_dofs = df.fem.locate_dofs_topological(S,fdim,interpolation_facets)
            else:
                integration_facets = self.facet_tag.find(boundary_marker)
                integration_dofs = df.fem.locate_dofs_topological(S,fdim,integration_facets)
                sorted_indices = np.argsort(S.tabulate_dof_coordinates()[integration_dofs][:,self.axis])
                dof_to_remove = integration_dofs[sorted_indices][endpoint_index]
                A = integration_dofs.tolist()
                A.remove(dof_to_remove)
                interpolation_dofs = np.array(A,dtype = np.int32)
            psi_h_bcs.append( df.fem.dirichletbc(psi_h_s, interpolation_dofs) )

        # Check if V_h used is nondivergent : FAILS
        print("fluxes of V_h along 4 boundaries = ", bf)
        print("net V_h flux (supposed to be ~ 0) = ", sum(bf))
        
        return psi_h_bcs

Nx,Ny = 60,60 # Number of points in x and y
family = "Lagrange"
order = 1
overlapping_facets = False
Lx, Ly = 4,4

mesh = df.mesh.create_rectangle(MPI.COMM_WORLD, [[-Lx,-Ly],[Lx,Ly]] ,[Nx, Ny])
tdim = mesh.topology.dim
fdim = tdim-1
mesh.topology.create_connectivity(tdim-1, tdim)

# Define normal and tangential unit vectors
n = ufl.FacetNormal(mesh)
t = ufl.as_vector([-n[1], n[0]])

V_space = df.fem.functionspace(mesh, (family, order, (2,))) # Vector space
wind = df.fem.Function(V_space)
wind.interpolate(lambda x: (3*(x[0]-x[1]),3*x[0]-x[1]))

S = df.fem.functionspace(mesh, (family, order)) # Scalar space

# ---- Poisson Equation: psi_I ----
# del^2 psi_I = -curl(wind)
# psi_I = 0 on boundary
Z_expr = df.fem.Expression(wind[1].dx(0) - wind[0].dx(1),S.element.interpolation_points)
Z = df.fem.Function(S)
Z.interpolate(Z_expr) # Z = curl(wind)
boundary_facets = df.mesh.exterior_facet_indices(mesh.topology)
boundary_dofs = df.fem.locate_dofs_topological(S, tdim-1, boundary_facets)
bc = df.fem.dirichletbc(df.default_scalar_type(0),boundary_dofs,S)

u = ufl.TrialFunction(S)
v = ufl.TestFunction(S)
a = ufl.inner(ufl.grad(u),ufl.grad(v)) * ufl.dx
L = -Z*v*ufl.dx
problem = LinearProblem(a=a,L=L,bcs = [bc], petsc_options_prefix="Streamfunction_I",petsc_options={"ksp_type": "preonly", "pc_type": "lu"}) 
psi_I = problem.solve()

# ---- Poisson Equation: chi_I ----
# del^2 chi_I = -divergence(wind)
# chi_I = 0 on boundary
d_expr = df.fem.Expression(ufl.div(wind),S.element.interpolation_points)
d = df.fem.Function(S)
d.interpolate(d_expr) # Evaluate divergence(wind)
L = -d*v*ufl.dx
problem = LinearProblem(a=a,L=L,bcs = [bc], petsc_options_prefix="Velocity_potential_I",petsc_options={"ksp_type": "preonly", "pc_type": "lu"}) 
chi_I = problem.solve()

# ---- V_h = wind - V_I ----
V_I = df.fem.Function(V_space)
V_I_expr = df.fem.Expression(ufl.perp(ufl.grad(psi_I)) + ufl.grad(chi_I), V_space.element.interpolation_points)
V_I.interpolate(V_I_expr)
V_h_true = df.fem.Function(V_space)
V_h_expr = df.fem.Expression(wind - V_I, V_space.element.interpolation_points)
V_h_true.interpolate(V_h_expr) # This variable will store the true V_h to compare to later

V_h = V_h_true.copy() # This variable will change during runtime (right now it is set to V_h_true)

# ---- Laplace Equation: psi_h ----
# del^2 psi_I = 0
# psi_h(s) = psi_h(s_0) - int_{s_0}^s dot(V_h, n) ds on boundary
bc_obj = bc_setter()
psi_h_s = bc_obj.set_psi_bcs()
f = df.fem.Constant(mesh, df.default_scalar_type(0)); L = -f*v*ufl.dx
problem = LinearProblem(a = a, L = L, bcs = psi_h_s, petsc_options_prefix="Streamfunction_h",petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
psi_h_0 = problem.solve()

Which is the correct approach? Do I let the boundary dofs interpolated to for each boundary condition overlap, or not?

2

You’re talking about setting Dirichlet BC’s on Lagrange basis functions, right? This specification is important as:

  • For many other spaces, the dofs do not lay on the mesh-nodes, so corners are not problematic.
  • For other kind of BC’s (Robin, Neumann), you integrate over facets, so a Dirichlet condition set at a corner does not conflict.

But then I’m afraid the answer is somewhat unsatisfactorily: it depends on what you want…

You can only set a nodal value once. If multiple Dirichlet BC’s target the same node, FEniCS can only enforce satisfaction of one of those conditions. I suppose it will enfore the last one in the list. Optimally, you’d solve problems where the boundary conditions from either side of the shared node are not in conflict. But this is not a mathematical requirement (Dicherlet conditions are allowed to jump). Then there is no general answer to which side should take precedence.

For indicating the boundaries, it is generally safest to tag the whole boundary (resulting in those ‘overlapping’ dofs), as this will be robust irrespective of mesh/order refinement, and number of decimal points in the meshing software. Maybe that’s what triggered your question?

3

That is an important delineation I was not aware of. I know very little about functionspaces in general. Is there a list showing which ones do not require me being conscious of the corners?

Oh, if that is how it works, then this will be much more simple to work around. Only my lower left corner (for my case) will need to be taken into account for my specific case.

I know how to tag the whole boundary all at once, but I do not know how to define my boundary condition over all the dofs all at once for my particular case. So, I stuck to doing it in pieces where I used geometric arguments to tag each of the 4 boundaries separately, and cumulatively integrated over them one at a time as you see in my code. Opting to do things that way triggered the question of “What happens at the points where the dofs overlap?”

For clarity, my boundary condition is defined as
\psi_h(s) = \psi_h(s_0) - \int_{s_0}^s \textbf{V}_h \cdot \hat{\textbf{n}} ds
where these are cumulative open-path boundary integrals over each of the 4 boundary paths, s. There might be a way to perform the cumulative integration over the full loop
\psi_h(s) = -\oint_{s_0}^s \textbf{V}_h \cdot \hat{\textbf{n}}ds
in Fenics (or Fenicsx), I just do not know how. If I could do it that way, my code would, yes, be more robust as it would cut down on a lot more uncertainty, and my code would also become much cleaner.

4

I would suggest having a look at:

and

which contains most of the spaces that have integral moments

1 Like
5

Ah, I didn’t mean the whole boundary, I indeed meant the whole boundary segment (as opposed to tagging the boundary segment minus a corner node). So precisely as you are doing.

1 Like