What is the meaning of FacetNormal for a mesh of triangles in 3D? Does it produce the normals to the facets of the triangles, i.e. the edges, or does it produce the normals to the triangles themselves?
I would suggest doing an experiment...
The reason I am asking is because I can't get my forms to compile, but I think it may be for another reason.
e.g.
RT = FunctionSpace(mesh, "RT", 1) V = TestFunction(RT) v0 = Function(RT) n = FacetNormal(mesh) print assemble(dot(V,n)*ds).array()
will work in 2D, but not for a 2D mesh in 3D. I get this:
f3663d84167207464.h: In member function ‘virtual void ffc_form_cf3a99dda1744 51d2a168d3f3663d84167207464_exterior_facet_integral_0_otherwise::tabulate_te nsor(double*, const double* const*, const double*, std::size_t) const’: /tmp/tmp7NpdaA2013-6-28-14-13_instant_ffeb0925148b9c5f074d5b4a5822570ca3eb7d fa/cf3a99dda174451d2a168d3f3663d84167207464/ffc_form_cf3a99dda174451d2a168d3 f3663d84167207464.h:789:9: error: ‘cell_orientation’ was not declared in thi s scope
Now, I think this must be specific to RT elements, as it seems OK with CG1 Vectors.
You could ask Martin Alnaes on mailing-list.
Has a resolution been found for this? The situation seems to be that you cannot solve even the simplest Dirichlet problem for the Poisson equation in mixed form on a manifold with boundary if the given Dirichlet data (meaning the values of the scalar variable) are not zero. For example, suppose we want to solve Laplacian p = 0 on a hemisphere, with p = 1 on the equator (so the exact solution is the constant 1. The mixed method solution is:
from dolfin import * mesh = Mesh("hemisphere_mesh.xml.gz") global_normal = Expression(("x[0]", "x[1]", "x[2]")) mesh.init_cell_orientations(global_normal) n = FacetNormal(mesh) V0 = FunctionSpace(mesh, "RT", 1) V1 = FunctionSpace(mesh, "DG", 0) V = MixedFunctionSpace([V0, V1]) (u, p) = TrialFunctions(V) (v, q) = TestFunctions(V) a = (dot(u, v) - p*div(v) - div(u)*q) * dx L = dot(v, n) *ds up = Function(V) solve(a == L, up)
(The needed mesh file can be downloaded from here.) But this fails in the instant compilation with the message Chris quoted: error: ‘cell_orientation’ was not declared in this scope.
What's the status with this? Is there a work-around?
Doug, I believe that issue with your code is distinct from both FFC issue 10 and 14 and should be reported.
Reported as issue 22.
FacetNormal represents the normal to the facet(s). In particular for a triangle (in 2D or 3D or nD), this would be the normal to the edges.
That makes sense - I think the problem I was having was more due to an ffc bug/issue at the time.
Nevertheless, there are an infinite number of vectors which are normal to the triangle edge in 3D. I suppose FacetNormal is therefore also in-plane with the triangle?
What are the implications for curved manifolds where the FacetNormal of two neighbouring cells on the same edge may be different? How does this work for R-T elements?
Yes, the normal in-plane with the triangle.
The implications are probably many. For RT for instance, I imagine that you'll only enforce normal continuity approximately.
