ufl
ufl copied to clipboard
FEniCS
Transformed form arguments
In this PR we propose a new class, ufl.Transformed, and several supporting classes.
Transformed will let finite element analysis programs such as FEniCS and Firedrake monolithicaly treat, e.g.:
- DirichletBC,
- EquationBC (domain dofs and boundary dofs satisfy distinct PDEs),
- Strong zero-normal velocity condition on tilted plane boundaries (Stokes problem).
Mathematically, Transformed represents a linear transformation, {v_} = {L}{v}, where {v} is a discretised (test/trial) function that is represented by a FromArgumant v defined on V,
and {L} is a linear transformation operator; for now we only consider cases where
{L} is diagonal, having the same discrete dimensionality as {v}, and in the following we denote
the diagonal entries of {L} by {l}.
In terms of the dimensionallity, {l} could be represented in UFL by a Coefficient (defined
on V), but the values in {l} are just scalar values attached to each basis, and
nothing more; values between two nodes, for instance, is of no interest, nor the gradient of this symbolic object.
This motivated us to introduce TopologicalCoefficient class that is like Coefficient, but has no knowledge of coordinates, and TopologicalCoefficient naturally necessitates TopologicalFunctionSpace and TopologicalMesh.
If we define V as:
V = FunctionSpace(msh, element),
we define:
tV = TopologicalFunctionSpace(tmsh, element),
where tmsh is a TopologicalMesh corresponding to msh. We can then have {l} represented by l: TopologicalCoeffient defined on tV.
The proposed usage of this new class in FEniCS/Firedrake is:
v_ = Transformed(v, l),
and a form compiler is to multiply {l}_i when i^th basis appear in the context of v; e.g. if v is a test function, {l}_i is to be multiplied when assembling the i^th row of a residual vector/jacobian matrix.
If values in {l} are 0 or 1, the above represents a projection of v onto a smaller finite element space spanned by basis functions whose associated values in {l} are set 1.
This is best illustrated in the following example:
Ex: Pseudo FEniCS/Firedrake code to solve poisson's equation with Dirichlet BC, u = g.
v = TestFunction(V)
u = TrialFunction(V)
l0 = TopologicalCoefficient( 1 for interior basis and 0 for boundary basis )
l1 = TopologicalCoefficient( 0 for interior basis and 1 for boundary basis )
v0 = Transformed(v, l0)
v1 = Transformed(v, l1)
a = inner(grad(u), grad(v0)) * dx + inner((u - g), v1) * ds
L = ...
solve(a == L, ...)
(By appropriately using Transformed(u, l0) and Transformed(u, l1), one can also symmetrise the system.)
The above example shows that Transformed lets users have a fine control of the basis set they assemble equations for. As such, Transformed, to some extent, can be regarded as a generalisation of Indexed, and indeed the implementation of Transformed mimics that of Indexed and Transformed is defined as a terminal_modifier.
Though one can solve the above example using DirichletBC class, Transformed will allow for higher flexibility in forming equations. For instance, defining:
a = B0(u, v0) * dx + B1(u, v1) * dx
in the above example lets one solve B0 in the domain and B1 on the boundary. Of note here is that the second term is only assembled into the rows of residual/Jacobian associated with boundary dofs, but it is a cell integration; this is useful when applying zero-normal velocity boundary condition in Stokes problem, as the PDE tested against the tangential component of test function has a cell integral, but it has to be solved on the boundary.
What is the mathematical meaning of this? You talk about a linear operator, ok, but then you do a magic trick and jump from an operator to FE function. Any concept in UFL should have a good mathematical meaning irrespective of the notion of degrees-of-freedom.
You say this is used to treat DirichletBC, and you provide an example. How is that a strongly enforced boundary condition? I see weakly imposed BC because your measure is still ds, you'd need to integrate wrt. a point measure to have a strong BC. Even if l1 has only one expansion coefficient 1 and all other zeros, the integral will have a contribution from the support of that one basis function.
I think this is best explained using function space decomposition. In the above example we decompose V as V = V0 + V1, where V0 is a space of functions that are zero on the boundary and V1 is the rest. We then decompose v accordingly as v = v0 + v1 where v0 = Transformed(v, l0) and v1 = Transformed(v, l1), l0 and l1 encoding this decomposition.
This again compares to Indexed, though there are some differences, too. A vector function, for instance, can be decomposed as v = (vx, 0) + (0, vy) and Indexed takes out vx or vy (unlike Transformed, Indexed also changes shape). In this sense I think l0 and l1 behaves like MultiIndex.
In the above Dirichlet BC is strongly applied via Mu = Mg instead of u = g on boundary, where M is some mass matrix composed by integration over the boundary. In my understanding a strong BC is applied when we have a perfect representation of g in V on the boundary. It often reduces to setting u = g pointwise, but this is not a requirement.
If values of l1 are all zero except one, we can probably use inner(u1, v1) * ds where u1 = Transformed(u, l1), but I would need to know the original mathematical expression of the boundary condition which would make us use such l1.
If you have V = V0 + V1 and the meaning of + there is a direct sum of vector spaces, then to decompose arbitrary v \in V what you need to do is a projection in V-norm, i.e. v0 = project(v, V0) and v1 = project(v, V1). In general, you can imagine a linear operator P: V -> V0 which represents the projection. Please note this was all in infinite-dimensional setting, no discretization chosen.
What you are trying to say is that you can discretize the projection operator P resulting in a diagonal matrix? Or that for each operator P there exists exactly one vector l0 such that the action of P can be represented with some action of l0?
Indexing in UFL and Indexed incl. is always about finite-dimensional spaces, as is R^2 in your "vector example".
Strictly speaking, V by itself must be zero on the whole Dirichlet BC boundary in order to be a valid function space (recall closedness to addition). Mathematically, we always formulate PDE on spaces with zero boundary conditions. Non-zero boundary value must be first lifted to the whole domain and then appears on the RHS (and it appears also with integrals over dx also, depends on the bilinear form).
Yes, I think your mathematical interpretation is right: the mathematical meaning of Transformed is a projection, and we have l carry data that encodes P, which becomes significant only when we discretise. So the idea is as the following. l carries tV just as v carries V, and l is a symbol that encodes basis selection just as v is a symbol that encodes test/trial function coeff. associated with each basis given V (even though we don't explicitly mention basis in UFL). When assembling, we finally make explicit mention to the basis, and l is realised as an array {l}, while v is realised as an array of basis coeffs. (of the same length). Then v_tran = Transformed(v, l) represents in discrete setting:
v_tran = \sum_i ( v_i * \phi_i * {l}_i ), where \phi_i is the i^th basis.
If values in {l} are 0 or 1, the above is a projection of v.
What I tried to mean was that Indexed still takes out infinite dimensional functions vx and vy, which compares to that Transformed takes out infinite dimensional (test/trial) function, say, v0, that is zero on boundary.
I'm sorry that I was careless when writing up the poisson's example. I had g==0 in my mind so a should be given as:
a = inner(grad(u), grad(v0)) * dx + inner(u , v1) * ds,
though one can also deal with nonzero g naturally.
I have updated the documentation to describe what each new class does and added an example to illustrate the usage of Transformed. Happy to discuss this further.
@ksagiyam Can you please explain again what is this about, but in terms of a pure mathematical language? I am still having difficulties to understand the mathematical concept behind.
If this is about projection, why is not project() enough? (and project() is not a part of UFL DSL).
@ksagiyam Can you please explain again what is this about, but in terms of a pure mathematical language? I am still having difficulties to understand the mathematical concept behind.
If this is about projection, why is not
project()enough? (andproject()is not a part of UFL DSL).
There's now several pages of maths in the documentation changes which are part of this PR. Perhaps you could explain what is unclear about that explanation so it can be improved.
Let us say u \in H^1(Omega) and let us have a compactly supported function which lives on a part of boundary of Omega, b \in L^2(Gamma).
What I can understand is a product tr(u)*b evaluated on the boundary, where tr: H^1(Omega) -> L^2(Gamma) is the trace operator. This product has the meaning of filtrating "u" on the boundary into the supp(b).
Now if we accept ()*ds implicitly means an application of a trace operator, how would Transformed be different from inner(u*b, v) * ds? Or maybe inner(u*b, v)*dP? where dP would be a sum of point measures which give meaning to how strong you would like the filtration be.
Is this product what you'd like to form?
And besides, projection into a subset of a mesh is not a dof-wise multiplication 0s and 1s that live on the boundary, I do not follow the argumentation that Transformed is a projection (which should involve solve against mass matrix, right?)
I, personally, do not think that there should be in UFL a concept that can be explained only using the notion of degree-of-freedom.
Thanks. I've kept talking about a standard Dirichlet problem as I thought it would best illustrate the feature, but one application that we are actually interested in is a Stokes problem with zero-normal velocity boundary condition on tilted plane boundary. Here, on the boundary, we have zero normal velocity, but the Stokes equation (domain equation) must still be tested against the tangential component of the test function, dot(v, t) t where v is a test function active on the boundary and t is the boundary tangent. We must use dx as the Stokes equations is a domain equation, but we cannot just say ... dot(v, t) t * dx as this cell integration would involve test functions that are not active on the boundary. Thus we need to do ... dot(Transformed(v, boundary), t) t * dx to only involve test functions active on the boundary, which, I believe, is something we cannot currently do equivalent thing of.
Dof-wise multiplication by 0s and 1s is indeed a projection; if the problem is of dimension N, and if there are M boundary dofs, we have:
v = \sum_{i=0}^{N-M} v_i \phi_i + \sum_{i=N-M}^N v_i \phi_i
where \phi_i represents a basis. Multiplying 0 to boundary dofs, we get:
v_ = \sum_{i=0}^{N-M} v_i \phi_i.
Transformed(v, topo_coeff) can be viewed as v transformed (or projected) to live on a modified function space, and topo_coeff encodes this function space modification, and this is all we see in a symbolic UFL form expression, and we do not mention degrees of freedom here.
A problem solving environment (and form compilers) realises the topo_coeff (defined on a TopologicalFunctionSpace) as a finite dimensional thing, just as it realises a Coefficient (defined on a FunctionSpace) as a finite dimensional thing.
So I think this new class fits well in UFL.
Thanks for the contribution @ksagiyam! I'm generally positive to exciting new features, but as @michalhabera , I still have a difficulty in understanding the concept that Transformed represents. Let's try to sort this out.
Above you write:
"We must use dx as the Stokes equations is a domain equation, but we cannot just say ... dot(v, t) t * dx as this cell integration would involve test functions that are not active on the boundary."
I do not understand the reasoning here. In particular, if t is the boundary normal, then integrating t over cells are not well defined in the first place.
Are topological mesh, topological function space and topological coefficients standard math terms? I'm not familiar with these - could you provide some references? Is a topological mesh a graph or a mesh topology, or something different?
Thank you, and please let me clarify this as I was careless.
Any vector field t that coincides with the boundary tangent on the boundary and makes dot(v, t) t live on a right space should be allowed at this point, as dot(v, t) t is meant to be a test function that has no degree of freedom in the normal direction on the boundary, and, even if we have zero-normal-velocity condition, the Stokes equation must still be tested against this test function.
A TopologicalMesh is meant to be a Mesh with fixed coordinate removed. I think we can view this as a graph with topological_dimension and geometric_dimension attached. Those terms, I believe, are not standard math words, but I attempted to explain them in the doc. I'm happy to discuss name changes, too.
Taking into account the comments from @michalhabera, we are aiming to come up with a better interpretation of Transformed class, and I am hoping that we can give an update on this sometime next week. Thank you.
