Create SubDomain based on Function values

mesh
1

Hello everyone,

My objective is to export a part of the mesh where a condition on the values of a function is valid.
To do so, I tried to create a SubMesh based on a SubDomain. However I don’t know how to create a SubDomain based on a condition depending on the values of a Function. I know how to create a SubDomain based on explicit spatial coordinates, but how should I proceed if I don’t know explicitly the expression defining my SubDomain of interest?

Here is a minimal working example of what I am trying to achieve:

from dolfin import *

import numpy as np
import matplotlib.pyplot as plt


nx = ny = 16
mesh = RectangleMesh(Point(-2, -2), Point(2,2), nx, ny)
V = FunctionSpace(mesh, 'CG', 2)
R = 1

# Expression chosen arbitrarily for the example (in practice I don't know explicitly this expression)
u_expr = Expression('(pow(x[0],2) + pow(x[1],2) <= pow(R,2))  ? 1.0 : 0.0', degree=2,R=R)
u = interpolate(u_expr, V)

plot(u)
plot(mesh)

class RegionOfInterest(SubDomain):
    def inside(self,x,on_boundary):
        # How to return the points where u == 1.0 ? if we don't know explicitly the expression defining u
        # return x where u==1.0
        return (pow(x[0],2) + pow(x[1],2) <= pow(R,2)) # explicit expression
    
marker = MeshFunction("size_t", mesh, mesh.topology().dim(), 0) # cells marker
region = RegionOfInterest()
region.mark(marker, 1) # marking left region cells with 1

submesh_region = SubMesh(mesh, region)

pl, ax = plt.subplots(); fig = plt.gcf(); fig.set_size_inches(16, 4)
plt.subplot(1, 2, 1); p = plot(mesh)
plt.subplot(1, 2, 2); p = plot(submesh_region)

xdmf_filename = XDMFFile(MPI.comm_world, "RegionOfInterest.xdmf")
xdmf_filename.write(submesh_region)

Don’t hesitate to propose another approach if you have better ideas on how to export parts of the mesh where a condition on a Function is valid.

Thank you for your time

2

Do the points corresponding to x in the inside method for a SubDomain, correspond to the coordinates of the vertices/Nodes of the mesh or some other dof?

3

In Paraview, there is a threshold filter, that will only show the domain in a specified range. I would strongly suggest using this capability, by saving your solution as a pvd or xdmf file.

1 Like
4

Thanks a lot, I’ll have a look at it.
It makes sense to do it while post-processing