I’m trying to iterate the facets connected to each vertex. My approach seems to work in serial execution but gets stuck (runs forever without any error messages) when run in parallel.
for i in range(len(x.vector())):
v = MeshEntity(mesh, 0, dof_to_vertex_map(V)[i])
for f in facets(v):
if f.exterior():
n = f.normal()
x.vector()[i] += n[0]
I also tried the other way round - going from facets to vertices - with the same result:
for f in facets(mesh):
if f.exterior():
n = f.normal()
for v in vertices(f):
i = vertex_to_dof_map(V)[v.index()]
x.vector()[i] += n[0]
Am I doing something wrong or does anybody know a solution?
Thanks and BR!
thanks for the quick reply! I adjusted the code (and will post a MWE at the end) but it doesn’t really achieve what I’m trying to do. The sum of the function values seems to only consider the facet attached to the vertex in the local scope - I hope the following picture does a better job explaining what I mean
I got it by running the code below with 2 processes and marked the assumed domain distribution with the yellow line.
I tried to just sum the values (of “glob_z” in my case) with MPI.sum but found out, that the indexes are different on each process.
So first of all: do you know a smart way to solve the problem?
… and secondly: why are the indexes different? Does “vertex_to_dof_map” not convert to global dof indexes? And if so, why does it work the other way round with “local_to_global_index(dof)”?
… and thirdly: can you recommend a documentation that explains the differences in scope (local/global), relation to indexes or just generally working with fenics in parallel?
Thanks again and BR!
My code example:
import dolfin as df
import numpy as np
comm = df.MPI.comm_world
rank = comm.rank
mesh = df.UnitCubeMesh(5, 5, 5)
mesh.init()
top = mesh.topology()
c02 = top(0,2)
V = df.FunctionSpace(mesh, "Lagrange", 1)
v2d = df.vertex_to_dof_map(V)
glob_z = np.zeros(V.dim())
entity_it = df.cpp.mesh.MeshEntityIterator(mesh, 0)
for vertex in entity_it:
faces = c02(vertex.index())
for face in faces:
facet = df.Facet(mesh, face)
if facet.exterior():
n = facet.normal()
glob_z[v2d[vertex.index()]] += n[2]
dm = V.dofmap()
my_first, my_last = dm.ownership_range()
unowned = dm.local_to_global_unowned()
dofs = list(filter(lambda dof: dm.local_to_global_index(dof) not in unowned, range(my_last-my_first)))
z = df.Function(V, name="z")
loc_z = z.vector().get_local()
loc_z[:] = glob_z[dofs]
z.vector().set_local(loc_z)
z.vector().apply("insert")
outfile = df.XDMFFile(comm, "z.xdmf").write(z)
Hi again @Maeggis.
The reason for it not setting the ghost values is because they are not inserted due to the unowned filter. I know how to access the ghosts of a FunctionSpace, and alter those dofs, but Im not sure how one can access the ghost vertices. I’ll think about if there is a clever way to do this.
On a slightly different subject. What is the usecase for your code? This Seems very much like a projection of the FacetNormal(mesh)[2] to a CG 1 space.
thanks again for your reply and effort!
I’m trying to pre-calculate a factor that depends on the incident angle (i.e. the surface normal and a direction defined as vector) so that I can use it in an iteration.
I tried an Expression approach like:
class BoundaryNormals(UserExpression):
def __init__(self, mesh, **kwargs):
super().__init__(kwargs)
self.mesh = mesh
def eval_cell(self, values, x, ufc_cell):
values[0] = 0.0
c = Cell(mesh, ufc_cell.index)
for facet in facets(c):
if facet.exterior():
n = facet.normal()
values[0] += n[2]
def value_shape(self):
return ()
But never got it to work correctly. So I went on to do it “manually” in order to directly apply the angular calculations - and, well, it worked but unfortunately not with multiple processors.
Following up on your suggestion with FacetNormal I tried the following MWE. It seems to work in serial execution (even though I expected a value range [-1,1] instead of [-2,2] but I can deal with that ) but running it in parallel gives “nan” values at some DOFs. I tried both interpolation functions. Is there something missing or wrong?
Thanks and BR!
import dolfin as df
mesh = df.UnitCubeMesh(5, 5, 5)
mesh.init()
V = df.FunctionSpace(mesh, "Lagrange", 1)
V3 = df.VectorFunctionSpace(mesh, "Lagrange", 1)
z = df.Function(V, name="z")
q = df.TestFunction(V)
h = df.FacetArea(mesh)
n = df.FacetNormal(mesh)
v = df.interpolate(df.Expression(('0', '0', '1'), degree=1), V3)
i = df.assemble(df.dot(v, n)*q/h*df.ds)
i = df.Function(V, i)
z = df.interpolate(i, V)
# df.LagrangeInterpolator.interpolate(z, i)
df.XDMFFile(df.MPI.comm_world, "z.xdmf").write(z)
import dolfin as df
class top(df.SubDomain):
def inside(self, x, on_boundary):
return x[2] > 1-1e-6 and on_boundary
mesh = df.UnitCubeMesh(10,10,10)
mesh.init()
mf = df.MeshFunction("size_t", mesh, mesh.topology().dim()-1)
mf.set_all(0)
top().mark(mf, 1)
# Custom measure for top plane, as the normal at sharp corners are really well defined.
dTop = df.ds(domain=mesh, subdomain_data=mf, subdomain_id=1)
V = df.FunctionSpace(mesh, "Lagrange", 2)
z = df.Function(V, name="z")
p,q = df.TrialFunction(V), df.TestFunction(V)
h = df.FacetArea(mesh)
n = df.FacetNormal(mesh)
L = df.assemble(n[2]*q*dTop)
A = df.assemble(p*q*dTop,keep_diagonal=True)
A.ident_zeros()
df.solve(A, z.vector(), L)
df.XDMFFile(df.MPI.comm_world, "z.xdmf").write(z)
I need it to work on various geometries without distinct surfaces but the results I got with your approach and “ds” instaed of “dTop” look really good.
The only issue is really sharp edges, where the facet normal is really gonna be weird to project onto a CG1 space. For slightly curved edges, you can just use boundary markers for sets of boundaries to get this.
