Taken we have the following function/solution of a previous simulation:
from mpi4py import MPI
from dolfinx.mesh import CellType, create_unit_square
from dolfinx.fem import functionspace, Function
from dolfinx.io import XDMFFile
from basix.ufl import element
import numpy as np
element_type = "Lagrange"
element_dof = 2
mesh = create_unit_square(
MPI.COMM_WORLD,
points=((0.0, 0.0), (1.0, 1.0)),
n=(N, N),
CellType.triangle
)
elem = element(element_type, mesh.basix_cell(), element_dof)
V = functionspace(mesh, elem)
psi = Function(V)
# Interpolate some random function
radius = 0.3
def phi(x): # in my usecase I dont actually know phi, and only have psi
return np.sqrt((x[0] - 0.5)**2 + (x[1]-0.5)**2) - radius
psi.interpolate(phi)
How would I mark all subdomains similar to the following pseudo code?
tol = 1e-8
if psi(x) >= tol:
# Domain 1
else:
# Domain 2
The post I am referring to utilized:
locate_entities
for this, but this requires to to have phi(x), which I dont have (or do I have to read out the value of psi for every point, which seems less optimal?)
This Depends a bit on how you want to quantity psi(x) on a domain. Since a domain is a cell Where psi(x) can be both positive and negative, you need to decide on some criteria.
Suggestion:
psi(midpoint(cell_i))>=0
N/M vertices of a cell satisfies psi(vertices_i)>=tol, i=1,…,M
average value of psi in cell>=tol
Each of these approaches requires slightly different implementations.
Sorry, that I was not clear. I meant higher order function space.
But would 3 not lead to the boundary between the domains being also the boundary of the cell?
Does this not lead to require a finer mesh size and not utilize the higher order function space?
Well, you said you wanted to have a marker for the mesh that uses something similar to locate_entities. Then you would have to mark whole cells or not.
What I meant, was that I want to define some subdomains based on a the value of a function. And “normaly” this seems to be made by marking the cells. if there is a better way like ufl.conditional I am happy to try it out.
@dokken I implemented ufl.conditional and tested it with the following minimal example:
from mpi4py import MPI
from dolfinx.mesh import CellType, create_unit_square
from dolfinx.fem import functionspace, Function
from dolfinx.io import XDMFFile
from basix.ufl import element
import numpy as np
element_type = "Lagrange"
element_dof = 1 # <-- increasing this should lead in a lower error
N = 90
mesh = create_unit_square(
MPI.COMM_WORLD,
nx=N, ny = N,
cell_type = CellType.triangle
)
elem = element(element_type, mesh.basix_cell(), element_dof)
V = functionspace(mesh, elem)
psi = Function(V)
# Interpolate some random function
radius = 0.3
def phi(x): # in my usecase I dont actually know phi, and only have psi
return np.sqrt((x[0] - 0.5)**2 + (x[1]-0.5)**2) - radius
psi.interpolate(phi)
import ufl
import dolfinx.fem as fem
tol = 1e-8
condition = ufl.lt(psi, tol)
true_statement = 1
false_statement = 0
f = ufl.conditional(condition, true_statement, false_statement)
dx = ufl.Measure("dx", domain=mesh)
markers = fem.form(f * dx)
area = fem.assemble_scalar(markers)
area_true = np.pi * radius**2
err_abs = np.abs(area_true-area)
print(err_abs)
What I would expect is that element_dof i.e 5 would result in a smaller error (err_abs). But as a matter effect the err_abs does not change at all, when I change element_dof and seems only be effected by ``N`! Did I miss something?
The accuracy of the area is quite close for any degree.
You can further investigate the behavior by varying the quadrature degree as well as the element degree:
import dolfinx.fem as fem
import ufl
from mpi4py import MPI
from dolfinx.mesh import CellType, create_unit_square
from dolfinx.fem import functionspace, Function
from dolfinx.io import VTXWriter
from basix.ufl import element
import numpy as np
def compute_area(N, degree, quadrature_degree):
metadata = {"quadrature_degree": quadrature_degree}
mesh = create_unit_square(
MPI.COMM_WORLD,
nx=N, ny=N,
cell_type=CellType.triangle
)
element_type = "Lagrange"
elem = element(element_type, mesh.basix_cell(), degree)
V = functionspace(mesh, elem)
psi = Function(V)
# Interpolate some random function
radius = 0.3
def phi(x): # in my usecase I dont actually know phi, and only have psi
return np.sqrt((x[0] - 0.5)**2 + (x[1]-0.5)**2) - radius
psi.interpolate(phi)
tol = 1e-8
condition = ufl.lt(psi, tol)
true_statement = 1
false_statement = 0
f = ufl.conditional(condition, true_statement, false_statement)
dx = ufl.Measure("dx", domain=mesh)
markers = fem.form(f * dx(metadata=metadata))
area = fem.assemble_scalar(markers)
area_true = np.pi * radius**2
err_abs = np.abs(area_true-area)
return err_abs
for degree in [1, 2, 3, 4, 5]:
for quadrature_degree in [1, 2, 4, 8]:
error = compute_area(
90, quadrature_degree=quadrature_degree, degree=degree)
print(degree, quadrature_degree, error)
