Access the vertex values just below/above the internal boundary

1

Hello all,
I have a question that I wish someone could help. Truly appreciate!
I have a 2-sphere mesh that is spited into two regions, which an internal boundary as their interface. I would like to access/set the all values of a function that is below/above the internal boundary (in this case, a horizontal cut with respect to the z-axis). I ran into several difficulties:

  1. I do not know how to access them with respect to different domains. I have tried to naively access the positive/negative sides across the boundary by using ("+/-"), but it seems does not work.
  2. The ("-") seems to have no values assigned.
    Below is the code.
    Thank you!
def sphere_mesh(outer_radius,mesh_density ):
    center = Point()
    sphere_shell = Sphere(center, outer_radius)#-Sphere(center, inner_radius)
    mesh = generate_mesh(sphere_shell, mesh_density)
    
    bmesh = BoundaryMesh(mesh, "exterior")
    print(bmesh.topology().dim())
    global_normal = Expression(("x[0]", "x[1]", "x[2]"), degree=1)
    bmesh.init_cell_orientations(global_normal)
    return bmesh




def mark_region(bmesh,z_loc,bc_tol=0.1,tol=1E-10):
    class Bc_region(SubDomain):
        def inside(self,x,on_boundary):
            return   (near(x[2],-z_loc,bc_tol) or near(x[2],z_loc,bc_tol)  )

    class Bulk_region(SubDomain):
        def inside(self,x,on_boundary):
            return (x[2]<=z_loc+tol and -z_loc-tol <= x[2])# or((near(x[2],-z_loc,bc_tol) or near(x[2],z_loc,bc_tol)  ))

    
    bc_region=Bc_region()
    bulk_region=Bulk_region()
    # define the boundary 
    boundary_parts =MeshFunction('size_t', bmesh, bmesh.topology().dim()-1, 0) #FacetFunction("size_t", self.bmesh)
    #print('bc dimension:'+str(self.bmesh.topology().dim()-1))
    boundary_parts.set_all(0)
    

    bulk_parts =MeshFunction('size_t', bmesh, bmesh.topology().dim(), 0) #CellFunction("size_t", self.bmesh)
    bulk_parts.set_all(0)
    bulk_region.mark(bulk_parts,1)
    bc_region.mark(boundary_parts,2)
    


    return boundary_parts, bulk_parts


b_mesh=sphere_mesh(outer_radius=1,mesh_density=30)
mark_region(b_mesh,0.5)
P=FiniteElement('DG',b_mesh.ufl_cell(),1)
element=MixedElement([P,P])
b_space=FunctionSpace(b_mesh,element)

Init_val=Expression(('x[1]*x[0]','x[2]*x[1]'),degree=1)
projected_fuc=project(Init_val,b_space)

#get the value in region   below the boundary
fuc_below_bc=projected_fuc("+")
print(project(fuc_below_bc,b_space).compute_vertex_values()) 
fuc_above_bc=projected_fuc("-")
print(project(fuc_above_bc,b_space).compute_vertex_values())
2

Have you considered having a look at: Integrating over an interior surface - #4 by MiroK

3

Hello @dokken ,
Yes, I am aware of the post. However, I am confused on getting the values below/above the interface in each subdomain respectively, not integrating along the boundary.

4

The point is that the “+” and “-” restrictions are only to be used with dS integrals, where they usually are random. In your code you are doing several things that does not make sense:

  1. use the restrictions “+” and “-” in project. Project is a function that integrates over dx, where restrictions does not matter.
  2. You use Function.compute_vertex_values(). This function evaluates the function at the vertices of the mesh, which is not well-defined when you are using a DG space.

I’ve attached a minimal working code that extracts the values along a boundary, and gather the in two groups, those attached to the cells with the marker 0 (bottom) and marker 1 (top).

import numpy as np
from dolfin import *


mesh = UnitSquareMesh(4, 4)

V = FunctionSpace(mesh, "DG", 0)

# Make a function that is 2 if x[1] >= 0.5 else 1
expr = Expression("x[1]>=0.5 ? 2 : 1", degree=0)

u = project(expr, V)


class Top(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] >= 0.5-1e-16


top_domain = Top()
cell_marker = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
top_domain.mark(cell_marker, 1)


class MidLine(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.5)


mid = MidLine()
facet_marker = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
mid.mark(facet_marker, 2)

top_vals = []
top_coords = []

bottom_vals = []
bottom_coords = []
x = V.tabulate_dof_coordinates()
for facet in facets(mesh):
    idx = facet.index()
    if facet_marker.array()[idx] == 2:
        connected_cells = facet.entities(mesh.topology().dim())
        cell_vals = cell_marker.array()[connected_cells]
        for cell, val in zip(connected_cells, cell_vals):
            u_dofs = V.dofmap().cell_dofs(cell)
            for dof in u_dofs:
                if val == 0:
                    bottom_vals.append(u.vector().get_local()[dof])
                    bottom_coords.append(x[dof])
                else:
                    top_vals.append(u.vector().get_local()[dof])
                    top_coords.append(x[dof])
                    print(u.vector().get_local()[dof], x[dof])

bottom_coords = np.vstack(bottom_coords)
top_coords = np.vstack(bottom_coords)

print("bottom values")
for coord, val in zip(bottom_coords, bottom_vals):
    print(coord, val)

print("top values")
for coord, val in zip(top_coords, top_vals):
    print(coord, val)

File("u.pvd") << u
5

Hello @dokken ,
That seems working, and thank you for the example!