How to define a Function in a submesh based on function values present in an adjoining submesh?

1

Hello,

I have two submeshes - 1 and 2 as shown below. They share a common boundary between them.

Domain_Mesh - Copy

I have a function u defined on the function space, Vusb1 in the top submesh. I have obtained this function by solving a PDE in this submesh.

u = Function(Vsub1)  

# ... I have then solved for this u on Vsub1 and obtained a solution

I am now trying to define a function u2 in the bottom submesh that takes the values of u at the shared boundary and zero everywhere else. Something like this:

u2 = Function(Vsub2)
values of u2(at the shared boundary) = values of u at the shared boundary (This u is defined in the top submesh)

I know the the co-ordinates of this shared boundary and I have also marked this boundary as 1 in both the submeshes.

I’m struggling to implement this.

Thanks in advance for the help offered !

2

Hi, Please add a minimal working code example, as described in Read before posting: How do I get my question answered?
Without such an example, it requires more effort from the community to try to help you with your problem, as one has to create example code and hopefully interpret your question correctly. Especially since SubMesh is a relatively new feature, a code example will help others that might face similar problems.

3

Dear @dokken

Thanks a lot for the reply and guidance. I have built a minimal working code that describes the issue for which I have requested help. The code is attached below:

I have a function u defined on the function space, Vsub in the top submesh. I have obtained this function by solving a PDE in this submesh.

from dolfin import *
from fenics import *
from ufl import nabla_div
from ufl import sinh

r = Constant(25.2)
g = -100
C = Constant(19.34)
kappa = 0.04

# create mesh
mesh = RectangleMesh(Point(0,0), Point(100e-6, 200e-6), 100,100)

regions = MeshFunction('size_t', mesh, mesh.topology().dim())
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)

class A(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] < 100e-6            #Bottom half of domain

a = A()

class B(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] >= 100e-6           #Top half of domain
b = B()

class TopBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 200e-6)
top = TopBoundary()

class MiddleBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near (x[1], 100e-6)
middle = MiddleBoundary()

regions.set_all(0)                     
b.mark(regions, 1)                      #Top region is marked 1; bottom remains marked as 0

boundaries.set_all(0)
middle.mark(boundaries, 2)            #Middle boundary in the domain
top.mark(boundaries, 1)                 #Top most boudnary in the domain


# Working on the top submesh - To solve a PDE

# define a submesh, composed of region 1 - top region 
new_region = MeshFunction('size_t', mesh, mesh.topology().dim())
new_region.set_all(0)
new_region.array()[regions.array() == 1] = 1
submesh = SubMesh(mesh, new_region, 1)

# define new meshfunctions on this submesh 
submesh_regions = MeshFunction('size_t', submesh, submesh.topology().dim())
submesh_regions.set_all(0)
submesh_boundaries = MeshFunction('size_t', submesh, submesh.topology().dim() - 1)
submesh_boundaries.set_all(0)

middle.mark(submesh_boundaries, 2)
top.mark(submesh_boundaries, 1)

ds = Measure("ds")[submesh_boundaries]

Vsub = FunctionSpace(submesh, 'CG', 1)
u = TrialFunction(Vsub)
v = TestFunction(Vsub)
u = Function(Vsub)
dx = Measure("dx")[submesh_regions]


bcs = []
F = kappa*inner(nabla_grad(u), nabla_grad(v))*dx + g*v*ds(1) + r*sinh(C*u)*v*ds(2)

J = derivative(F, u)
problem = NonlinearVariationalProblem(F, u, bcs, J)
solver  = NonlinearVariationalSolver(problem)
solver.solve()

#Plotting solution u obtained on Vsub - defined on top submesh 
import matplotlib.pyplot as plt
from matplotlib import cm

p = plot(u, cmap = cm.inferno)
plt.colorbar(p)
plt.show()

# Working on the bottom submesh 

# define a submesh, composed of region 0 - bottom region 
new_region2 = MeshFunction('size_t', mesh, mesh.topology().dim())
new_region2.set_all(0)
new_region2.array()[regions.array() == 0] = 1
submesh2 = SubMesh(mesh, new_region2, 1)

# define new meshfunctions on this submesh 
submesh_regions2 = MeshFunction('size_t', submesh2, submesh2.topology().dim())
submesh_regions2.set_all(0)
submesh_boundaries2 = MeshFunction('size_t', submesh2, submesh2.topology().dim() - 1)
submesh_boundaries2.set_all(0)

middle.mark(submesh_boundaries2, 2)

ds2 = Measure("ds")[submesh_boundaries2]

Vsub2 = FunctionSpace(submesh2, 'CG', 1)
u2 = TrialFunction(Vsub2)
v2 = TestFunction(Vsub2)
dx2 = Measure("dx")[submesh_regions2]

u2 = Function(Vsub2)

Following this, I need to define a function u2 in the bottom submesh that takes the values of u at the shared boundary and zero everywhere else. Something like this:

values of u2(at the shared boundary) = values of u at the shared boundary (This u is defined in the top submesh)

Since, both the boundaries are marked and share the same co-ordinates, I thought of extracting the dofs and modifying the function u2 as follows:

from numpy import where
top_dofs = where(submesh_boundaries.array()==2)
bottom_dofs = where(submesh_boundaries2.array()==2)

u2.vector()[bottom_dofs[0]] = u.vector()[top_dofs[0]]

It shows the following error.

u2.vector()[bottom_dofs[0]] = u.vector()[top_dofs[0]]
IndexError: Vector index out of range

I’m doing some fundamental mistake in assigning the function values based on the dofs, and not able to update the function u2 properly.

Please help me resolving this issue. Thanks in advance !

4

Dear @Sabarish_Gopi,

As this is not my field of expertise, my solution to this problem is highly specialized to CG-1 spaces.
Maybe @cdaversin has a better solution than this.

Please note that your example is far from a minimal example,as you have included solving a non-linear PDE. Let the following code be a guide on how to produce a Minimal code example in the future.

from  dolfin  import *
import numpy as np
n = 10

mesh = UnitSquareMesh(n, n)

# Divide a mesh into two parts, x[0]<0.7 and x[0]>=0.7
marker = MeshFunction("size_t",mesh ,mesh.topology().dim(),0)
for c in  cells(mesh):
    marker[c] = 0.7 < c.midpoint().x()
marker_array = marker.array()
submesh0 = MeshView.create(marker ,0)
submesh1 = MeshView.create(marker ,1)

# Initialize  function  spaces  and  basis  functions
V0 = FunctionSpace(submesh0 , "CG", 1)
V1 = FunctionSpace(submesh1 , "CG", 1)
W = MixedFunctionSpace(V0,V1)

# Create empty function, and interpolate u=x[1] onto part 0
u = Function(W)
u.sub(0).interpolate(Expression("x[1]", degree=1))

# Create mapping from submesh to parent mesh
mv0 = submesh0.topology().mapping()[mesh.id()]
vertex_map0 = np.array(mv0.vertex_map(),dtype=np.int32)
mv1 = submesh1.topology().mapping()[mesh.id()]
vertex_map1 = np.array(mv1.vertex_map(), dtype=np.int32)
# Find which vertices are in common for the two submeshes
common_vertices = np.flatnonzero(np.isin(vertex_map0, vertex_map1))

dofs0 = vertex_to_dof_map(V0)
dofs1 = vertex_to_dof_map(V1)
# Map dofs from submesh 0 to submesh 1 through their common vertex.
# To do this with other function-spaces than CG-1, mappings
# between facets has to be added
for vertex in common_vertices:
    vertex_0 = vertex_map0[vertex]
    dof0 = dofs0[vertex]
    vertex_1 = np.flatnonzero(vertex_map1 == vertex_0)[0]
    dof1 = dofs1[vertex_1]
    u.sub(1).vector().vec().setValueLocal(dof1,  u.sub(0).vector()[dof0])

File("u0.pvd") << u.sub(0)
File("u1.pvd") << u.sub(1)
1 Like
5

Dear @dokken,

Thanks a lot for spending time to help me. I now understand how to provide a minimal working code.

In the code you have shared, it shows the following error for MeshView and MixedFunctionSpace.

NameError: name 'MeshView' is not defined

Any leads on what I’m missing ? Thanks again !

6

To use the meshview functionality was added as part of the mixed dimensional FE implementation in FEniCS (see arxiv paper).

This is added to the development version of dolfin, and can for instance be found with docker:

docker run -ti -v $(pwd):/home/fenics/shared -w /home/fenics/shared/ --rm quay.io/fenicsproject/dev:latest
1 Like
7

Dear @dokken

Is there a way around the usage of MeshView and MixedFunctionSpace provided just by mapping the dofs from submesh0 to submesh1 through the common vertices.
I have split the domain into two submeshes, also I have a function u defined on one of the submeshes.

If not, how to get started with the docker ? I’m very new to it.

Thank you very much.

8

As i have not used submeshes I do not know. There is probably possible to do something similar to what I’ve done above. If you can make your problem similar to the code I posted I might Give it a go. Your current problem is too big for me to take time in reducing it.
There is a long tutorial on docker here: https://fenics.readthedocs.io/projects/containers/en/latest/introduction.html#installing-docker

1 Like
9

Hi! I have a question about modeling: if my model consists of two regions with a common boundary and I want to use FEniCS to calculate heat transfer or mechanical problems, when drawing the mesh using Gmsh, should the common boundary be drawn as a single edge or as two overlapping edges (such as the interface between a solid and a liquid)? What are the differences between these two methods? If it should be drawn as two overlapping edges, then how can I match these two edges in FEniCS? I hope to get an answer, thank you!

10

I recently encountered a similar problem. Has the poster solved this problem? I’m looking forward to your reply.

11

This can for instance be done with a common interface in DOLFINx, with the latest developments (v0.9) on submeshes. See for instance: Discontinuity at interface using mixed domains - #2 by dokken

12

Thank you very much for your early reply! In my recent test, I found that for the two-dimensional case, points on adjacent boundaries can be included in two sub-grids, but the three-dimensional case seems to be a little different? Here is my two-dimensional test example:

import gmsh
import os
from dolfinx.io import gmshio
from mpi4py import MPI
from dolfinx.mesh import create_submesh
import numpy as np

# Create output directory
output_dir = "./data"
if not os.path.exists(output_dir):
    os.makedirs(output_dir)

gmsh.initialize()
gmsh.model.add("two_regions_3d")

# Create geometry
outer_radius = 1.0
inner_radius = 0.3

region1_full = gmsh.model.occ.addDisk(0, 0, 0, outer_radius, outer_radius)
temp_cut = gmsh.model.occ.addDisk(0, 0, 0, inner_radius, inner_radius)
region2 = gmsh.model.occ.addDisk(0, 0, 0, inner_radius, inner_radius)

region1 = gmsh.model.occ.cut([(2, region1_full)], [(2, temp_cut)])
gmsh.model.occ.synchronize()

region1_tag = region1[0][0][1]

# Set mesh parameters and create physical groups
mesh_size = 0.05
gmsh.model.mesh.setSize(gmsh.model.getEntities(0), mesh_size)

gmsh.model.addPhysicalGroup(2, [region1_tag], 11111)
gmsh.model.addPhysicalGroup(2, [region2], 22222)

# Mark boundaries
boundary1 = gmsh.model.getBoundary([(2, region1_tag)], oriented=False)
boundary2 = gmsh.model.getBoundary([(2, region2)], oriented=False)

boundary1_tags = [tag for dim, tag in boundary1]
boundary2_tags = [tag for dim, tag in boundary2]

gmsh.model.addPhysicalGroup(1, boundary1_tags, 44444)  # outer boundary
gmsh.model.addPhysicalGroup(1, boundary2_tags, 33333)  # interface

gmsh.model.mesh.generate(2)
gmsh.write(os.path.join(output_dir, "two_regions_2d.msh"))
gmsh.finalize()

# Read mesh and process
file = "./data/two_regions_2d.msh"
domain, cell_markers, facet_markers = gmshio.read_from_msh(file, MPI.COMM_WORLD, gdim=2)

# Create submeshes
tdim = domain.topology.dim
submesh1, map1, vertex_map1, geom_map1 = create_submesh(domain, tdim, cell_markers.find(11111))
submesh2, map2, vertex_map2, geom_map2 = create_submesh(domain, tdim, cell_markers.find(22222))

# Get interface vertices
fdim = tdim - 1
domain.topology.create_connectivity(fdim, 0)
face_to_vertex = domain.topology.connectivity(fdim, 0)

interface_vertices = []
interface_facets = facet_markers.find(33333)
for facet in interface_facets:
    vertices = face_to_vertex.links(facet)
    interface_vertices.extend(vertices)
interface_vertices = np.unique(interface_vertices)

# Print mesh information
print(f"Total vertices in original mesh: {domain.geometry.x.shape[0]}")
print(f"Vertices in region1: {submesh1.geometry.x.shape[0]}")
print(f"Vertices in region2: {submesh2.geometry.x.shape[0]}")
print(f"Number of interface vertices: {len(interface_vertices)}")

# Check interface vertices membership
def find_matching_vertices(coords, reference_coords, tolerance=1e-5):
    matches = 0
    for ref_coord in reference_coords:
        distances = np.linalg.norm(coords - ref_coord, axis=1)
        if np.any(distances < tolerance):
            matches += 1
    return matches

coords1 = submesh1.geometry.x
coords2 = submesh2.geometry.x
interface_coords = domain.geometry.x[interface_vertices]

matches1 = find_matching_vertices(coords1, interface_coords)
matches2 = find_matching_vertices(coords2, interface_coords)

print(f"\nInterface vertices analysis:")
print(f"- {matches1} vertices found in region1")
print(f"- {matches2} vertices found in region2")

The result is as follows:
Number of interface vertices: 38

Interface vertices analysis:

  • 38 vertices found in region1
  • 38 vertices found in region2

but for the 3D case, the situation is different:

import gmsh
import os
from dolfinx.io import gmshio
from mpi4py import MPI
from dolfinx.mesh import create_submesh
import numpy as np

# Create output directory
output_dir = "./data"
if not os.path.exists(output_dir):
    os.makedirs(output_dir)

gmsh.initialize()
gmsh.model.add("two_regions_2d")

# Create geometry
outer_radius = 1.0
inner_radius = 0.3

region1_full = gmsh.model.occ.addSphere(0, 0, 0, outer_radius, 1)
temp_cut = gmsh.model.occ.addSphere(0, 0, 0, inner_radius, 2)
region2 = gmsh.model.occ.addSphere(0, 0, 0, inner_radius, 3)

region1 = gmsh.model.occ.cut([(3, 1)], [(3, 2)], 4)
gmsh.model.occ.synchronize()

region1_tag = region1[0][0][1]

# Set mesh parameters and create physical groups
mesh_size = 0.1
gmsh.model.mesh.setSize(gmsh.model.getEntities(0), mesh_size)

gmsh.model.addPhysicalGroup(3, [region1_tag], 11111)
gmsh.model.addPhysicalGroup(3, [region2], 22222)

# Mark boundaries
boundary1 = gmsh.model.getBoundary([(3, region1_tag)], oriented=False)
boundary2 = gmsh.model.getBoundary([(3, region2)], oriented=False)

boundary1_tags = [tag for dim, tag in boundary1]
boundary2_tags = [tag for dim, tag in boundary2]

gmsh.model.addPhysicalGroup(2, boundary1_tags, 44444)  # outer boundary
gmsh.model.addPhysicalGroup(2, boundary2_tags, 33333)  # interface

gmsh.model.mesh.generate(3)
gmsh.write(os.path.join(output_dir, "two_regions_3d.msh"))
gmsh.finalize()

# Read mesh and process
file = "./data/two_regions_3d.msh"
domain, cell_markers, facet_markers = gmshio.read_from_msh(file, MPI.COMM_WORLD, gdim=3)

# Create submeshes
tdim = domain.topology.dim
submesh1, map1, vertex_map1, geom_map1 = create_submesh(domain, tdim, cell_markers.find(11111))
submesh2, map2, vertex_map2, geom_map2 = create_submesh(domain, tdim, cell_markers.find(22222))

# Get interface vertices
fdim = tdim - 1
domain.topology.create_connectivity(fdim, 0)
face_to_vertex = domain.topology.connectivity(fdim, 0)

interface_vertices = []
interface_facets = facet_markers.find(33333)
for facet in interface_facets:
    vertices = face_to_vertex.links(facet)
    interface_vertices.extend(vertices)
interface_vertices = np.unique(interface_vertices)

# Print mesh information
print(f"Total vertices in original mesh: {domain.geometry.x.shape[0]}")
print(f"Vertices in region1: {submesh1.geometry.x.shape[0]}")
print(f"Vertices in region2: {submesh2.geometry.x.shape[0]}")
print(f"Number of interface vertices: {len(interface_vertices)}")

# Check interface vertices membership
def find_matching_vertices(coords, reference_coords, tolerance=1e-5):
    matches = 0
    for ref_coord in reference_coords:
        distances = np.linalg.norm(coords - ref_coord, axis=1)
        if np.any(distances < tolerance):
            matches += 1
    return matches

coords1 = submesh1.geometry.x
coords2 = submesh2.geometry.x
interface_coords = domain.geometry.x[interface_vertices]

matches1 = find_matching_vertices(coords1, interface_coords)
matches2 = find_matching_vertices(coords2, interface_coords)

print(f"\nInterface vertices analysis:")
print(f"- {matches1} vertices found in region1")
print(f"- {matches2} vertices found in region2")

The result is:

  • 14 vertices found in region1
  • 160 vertices found in region2

How can I make the two sub-meshes have the same nodes on adjacent boundary in 3D so that I can pass the boundary conditions?

13

You should use transfer meshtags, as for instance added in scifem:
https://scientificcomputing.github.io/scifem/examples/transfer_tags_to_submesh.html
to transfer facet-tags, vertex-tag or edge tags to the submesh. If you really want to match vertex by vertex, transfering a vertex tag would be the solution.

What kind of boundary conditions do you want to pass?

14

Thank you very much for your advice. I am dealing with a Cauchy problem, so both the Dirichlet condition and the Neumann condition are needed to pass.

15

Neither would need the vertex coordinates, as long as you can transfer the tags of the relevant entity to the sub mesh (at least without any further detail about the coupling at the interface).