I am trying to solve a simple diffusion problem on a 1D IntervalMesh with different diffusion coefficients in two different regions. For this, I am following the tutorial Defining subdomains for different materials using locate_entities(). However, for some reason, one of my cells is not marked by the following code:
# Define tolerance
tol = 1E-14
# Define mesh quantities
mesh_n = 101
mesh_min = -2E-4
mesh_max = 1E-3
# Create mesh and define function space
mesh = create_interval(MPI.COMM_WORLD, mesh_n, (mesh_min, mesh_max))
# Define subdomains for c-region and a-region
def subdomain_c(x):
return x[0] <= tol
def subdomain_a(x):
return x[0] >= -tol
# Mark cells based on their correspondence to a region
cells_subdomain_c = locate_entities(mesh, mesh.topology.dim, subdomain_c)
cells_subdomain_a = locate_entities(mesh, mesh.topology.dim, subdomain_a)
print(cells_subdomain_c)
print(cells_subdomain_a)
So obviously, the dof #16 does not correspond to any of the two regions, even though its x-coordinate should return true for one of the two functions subdomain_c/a(). What did I miss here and how do I code this in a more robust way?
which shows you that the 16th cell consist of the vertices with coordinate -9.90099010e-06 and 1.98019802e-06
which means that the cell only satisfies subdomain_a for the 1.98019802e-06, while subdomain_c only holds for -9.90099010e-06.
To be marked with locate_entities, the condition has to hold for all vertices.
Thank you for clearing this up @dokken! I’ve finally had some time to work on this a bit, however, I am still not 100% happy with how I am handling this. Basically, I see two solutions for my problem.
A) I can assign a magic number to mesh_n, e.g. 102, to make sure that there are vertices at exactly x = 0, in which case my code works as expected. However, magic numbers are usually to be avoided and the code shouldn’t just work for certain combinations of mesh_min, mesh_max and mesh_n.
B) I can set my tolerance to half the cell size, i.e. tol = (mesh_max - mesh_min) / (2 * mesh_n) to capture the vertices properly. However, there is an edge case for a cell with vertices at +tol and -tol, which would return true for both subdomains. Also, I do not necessarily have a mesh with constant cell size at all times, in which case this approach would fail.
I feel like there should be a more elegant way to do what I want, but I am not experienced enough to see it. My question, essentially, is how to implement a robust solution to my problem that doesn’t depend on magic numbers and is free of edge cases.