- edited description
Dirichlet boundary conditions of the form dot(b, n): Integral of type cell cannot contain a ReferenceNormal and segfault
The following code produces a segfault:
from dolfin import *
b = Constant([1.0, 2.0])
mesh = UnitSquareMesh(20, 20)
V = FunctionSpace(mesh, 'CG', 1)
n = FacetNormal(mesh)
bcs = DirichletBC(V, dot(b, n), 'on_boundary')
Right before that:
Integral of type cell cannot contain a ReferenceNormal.
Comments (17)
-
reporter -
This is a compound issue. It is the projection of a FacetFunction that fails. The boundary form dot(b,n) gets projected using dot(b,n)vdx. Changing dx to ds will pass the compilation and assemble. So we could expect an error message in UFL. Then it would be nifty to be able to project a boundary form to a Function.
The former is a bug either in UFL or in DOLFIN (Martin?), the latter is a feature request.
-
The former is this one:
https://bitbucket.org/fenics-project/ufl/issue/6/add-error-message-if-facet-normal-is-used
The latter is a feature request for generating expressions from ufl code that can be evaluated in boundary dofs. Something like that is on my horizon but not my highest priority.
-
- marked as minor
-
Note that normal can't be evaluated also in cell integrals on boundary mesh (where it could work)
from dolfin import * mesh = UnitSquareMesh(4, 4) bmesh = BoundaryMesh(mesh, 'exterior') n = FacetNormal(mesh) assemble(inner(n, n)*dx, mesh=bmesh)
error: ‘n1’ was not declared in this scope error: ‘n0’ was not declared in this scope
-
In the meantime the enclosed hack can be used for a calculation of RT1 outter normal which can be used in DirichletBC
from dolfin import * #mesh = UnitSquareMesh(4, 4) # Produces funny values in corners mesh = EllipseMesh(Point(1.0, 2.0), [1.0, 2.2], 0.2) RT1 = FunctionSpace(mesh, 'RT', 1) n = FacetNormal(mesh) # Not working #n = project(n, RT1) #n = interpolate(n, RT1) n_trial, n_test = TrialFunction(RT1), TestFunction(RT1) a = inner(n_trial, n_test)*ds L = inner(n, n_test)*ds A = assemble(a, keep_diagonal=True) A.ident_zeros() # Regularize the matrix b = assemble(L) #parameters['allow_extrapolation'] = True # May be needed for evaluation tests below n = Function(RT1) #parameters['allow_extrapolation'] = False solve(A, n.vector(), b, 'cg', 'amg') # Check values on facet midpoints #for f in facets(mesh): # if f.exterior(): # p = f.midpoint() # print f.index(), [v.index() for v in vertices(f)], n(p) # Check values at vertices #for f in facets(mesh): # if f.exterior(): # for v in vertices(f): # p = f.midpoint() # print v.index(), n(p) # Check correctness of mean value area = assemble(Constant(1.)*ds, mesh=mesh) nds = assemble(inner(n, n)*ds) print "Average value of RT1 normal on boundary:", nds / area print "Plotting. Note that RT1 field is first interpolated onto CG1." plot(n) interactive() print "Now you have approximate RT1 outter normal field which can be used in DirichletBC."
-
Thanks for the workaround Jan! Are you sure that this is working correctly in 3-D as well? When I replace
mesh = EllipseMesh(Point(1.0, 2.0), [1.0, 2.2], 0.2)
withmesh = SphereMesh(Point(1.0, 2.0, 0.0), 1, 0.1)
the normal vectors in the plot look odd... -
@taerath
-
RT1 normal is definitely only an approximation to proper normal field which is piece-wise constant field on facets.
-
Plots of RT fields may be funny as interpolation to plottable CG1 done behind the scene is often oscillatory near boundary to my opinion. See also this question by @dnarnold.
-
-
I believe there are at least two quite different problems here. One is that it is not possible to compile expressions like
dot(v, n)*ds
where v is an RT on a manifold and n is the facet normal. This means that you cannot solve a mixed Poisson problem with inhomogeneous Dirichlet data on manifold. This is the issue discussed here.
The second issue is not a problem of compilation. The code runs, but gives incorrect results when solving a mixed Poisson problem with inhomogeneous Neumann on a manifold. This issue is discussed here.
-
@dnarnold But it is not obvious that any of your issues make @taerath's normals looking funny. Remember, the example does not concern manifold.
-
@blechta I need the facet normals to calculate the flux of a vector field through a surface. When I'm using the above workaround, the results look very noisy (compared to the results from a different code), so I don't think the problem is just with the plotting. Is there a particular reason to approximate the normals with RT? For testing, I replaced
FunctionSpace(mesh, 'RT', 1)
withVectorFunctionSpace(mesh, 'CG', 2)
. Both the plot and the results of my flux calculation look a lot better with that (though still not perfect). -
@taerath No good reason generally. Just it may be in a good duality to what you're multiplying and integrating.
Yes, I'm not fully convinced that these are just post-processing oscillations. Maybe the procedure for calculation of RT1 normal introduces spurious oscillations or there is some bug connected with Doug's examples. But these seem to be manifold exclusive.
-
For what it's worth, here is a different workaround that I came up with. At least for my specific case, it gives much better results. It probably makes all sorts of implicit assumptions about my mesh and problem though:
from dolfin import * import numpy as np def get_facet_normal(bmesh): '''Manually calculate FacetNormal function''' if not bmesh.type().dim() == 2: raise ValueError('Only works for 2-D mesh') vertices = bmesh.coordinates() cells = bmesh.cells() vec1 = vertices[cells[:, 1]] - vertices[cells[:, 0]] vec2 = vertices[cells[:, 2]] - vertices[cells[:, 0]] normals = np.cross(vec1, vec2) normals /= np.sqrt((normals**2).sum(axis=1))[:, np.newaxis] # Ensure outward pointing normal bmesh.init_cell_orientations(Expression(('x[0]', 'x[1]', 'x[2]'))) normals[bmesh.cell_orientations() == 1] *= -1 V = VectorFunctionSpace(bmesh, 'DG', 0) norm = Function(V) nv = norm.vector() for n in (0,1,2): dofmap = V.sub(n).dofmap() for i in xrange(dofmap.global_dimension()): dof_indices = dofmap.cell_dofs(i) assert len(dof_indices) == 1 nv[dof_indices[0]] = normals[i, n] return norm mesh = SphereMesh(Point(1.0, 2.0, 0.0), 1, 0.1) bmesh = BoundaryMesh(mesh, 'exterior') n = get_facet_normal(bmesh) # Check values on facet midpoints for f in facets(bmesh): p = f.midpoint() print f.index(), [v.index() for v in vertices(f)], n(p) # Check correctness of mean value area = assemble(Constant(1.)*dx, mesh=bmesh) nds = assemble(inner(n, n)*dx) print "Average value of RT1 normal on boundary:", nds / area plot(n) interactive()
-
- changed milestone to 1.6
-
-
- removed milestone
Removing milestone: 1.7 (automated comment)
-
reporter - edited description
- changed title to Dirichlet boundary conditions of the form dot(b, n): Integral of type cell cannot contain a ReferenceNormal and segfault
Update MWE
- Log in to comment