Subject: [UFL-dev] Dirichlet boundary conditions

From: Kent-Andre Mardal <[email-masked]>
Date: Thu, 19 Jun 2008 11:24:41 +0200

Should we include Dirichlet conditions in UFL ? 

It could be done with the notion of traces: 

g = Function(...)
u = TrialFunction(...)
T = FacetTrace(...)

T*u = g 

meaning that u = g on the boundary. 
One could also have T(i) for the trace
on subdomain i. 

I think it also could be put into the UFC assembly loop in a nice way, 

assembly(...){

  ...

  assembleCells(...);
  
  assembleExteriorFacets(...);

  assembleInteriorFacets(...);

  enforceTraces(...);

}

The bad thing is that it affects both the right and left hand side so
that the matrix and the right hand side need to be assembled at the same
time (or one could make a very sparse matrix that resemble applying the
constraints). 

Kent

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Thu, 19 Jun 2008 11:38:37 +0200

>From a programming point of view, this line makes no sense to me:

> T*u = g

What it means mathematically is clear, but this can never be valid Python.

Do you mean something like:

a = u*v*dx + dot(u,n)*v*ds(0) + g*T(1)

where T(i) and ds(i) would refer to the same boundary?

> time (or one could make a very sparse matrix that resemble applying the

I don't think they'll _need to_ be assembled at the same time, this could
also be ignored during initial assembly and then done as a separate step.
In particular if it's only part of UFL and not UFC: the Assembler shouldn't
know about UFL at all, only UFC (which we may of course extend as well!).

--
Martin

---

From: Kent-Andre Mardal <[email-masked]>
Date: Thu, 19 Jun 2008 11:55:16 +0200

> What it means mathematically is clear, but this can never be valid Python.

You're right. Maybe 
T*u == g. 
Is this possible in Python ? It is in C++... 

> where T(i) and ds(i) would refer to the same boundary?

Could be but I think this is less mathematically clear 
(for one thing a would not be bilinear)

c = T*u == g. 
a = u*v*dx + dot(u,n)*v*ds(0)

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Thu, 19 Jun 2008 12:17:07 +0200

> Is this possible in Python ? It is in C++...

The assignment operator doesn't exist in Python.

Even if it did work, say in C++ or using ==, "T*u" would a a temporary object,
and "(T*u) == g" would be a call to (T*u).operator==(g) returning another
temporary object, which would then dissapear on the next line.

> a = u*v*dx + dot(u,n)*v*ds(0)

So you don't want the BC to be part of the form.
Where should it be used? Why should it be part of UFL
and not just a dolfin component?

If we are to include this in UFL, it should support arbitrary UFL expressions
in place of g, and we would also have to support component-wise conditions:

just_some_expression = (dot(g, n) + inner(grad(g), grad(g)))
c = T*(u[0]) == just_some_expression

This means we must generate code for it, and need to extend UFC to include it.

If we only support the basic "u = g" condition, I don't see the point
with adding it to UFL.

-- 
Martin

---

From: Kent-Andre Mardal <[email-masked]>
Date: Thu, 19 Jun 2008 13:43:55 +0200

> So you don't want the BC to be part of the form.

It is a form, just not a bilinear one. 

> and not just a dolfin component?

Could be that we should leave Dirichlet bc as it is. It's just
that T is in some respect very similar to ds. 

> c = T*(u[0]) == just_some_expression

This makes sense, eg in elasticity you have Neumann in one component and
Dirichlet in another on the same sub-domain. 

> If we only support the basic "u = g" condition, I don't see the point

I am not sure how much this affect UFC, but it could for instance be
nice to write

T*u*normal = g 

where g is a Function and normal is the outward normal. 

But anyway, this is just a thought that has not been thought properly
through yet... 

Kent

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Thu, 19 Jun 2008 15:39:58 +0200

> It is a form, just not a bilinear one.

UFL will have, as FFC has, operators "lhs(form)" and "rhs(form)"
which extracts the bilinear and linear integrals from a form.
It is possible that these could "do the right thing" with a form
like the one above including DBC.

We could collect them more explicitly:
  a = u*v*dx
  L = f*v*ds
  bc = T*u == g
  Axb = System(a, L, bc)

> that T is in some respect very similar to ds.

Don't misunderstand me, I think it is a good idea,
I'm just trying to filter out all problems I see.

We must keep in mind that dolfin is a C++ library,
UFL is a Python implementation, and dolfin should
only relate to UFC, not UFL directly. The exception
is perhaps a thin layer in PyDOLFIN where jit is called
from a form compiler to convert from UFL to UFC.
It is important to stick with the interfaces we agree on.

> where g is a Function and normal is the outward normal.

(As above this syntax isn't possible in the embedded language approach
we're using.
If we designed our own language from scratch (like mython?), that
could be possible.
Basically, we need to embed the information we wish to pass around into objects,
restricted by the syntax Python provides.)

That aside, we would have to define something flexible enough, which
may be difficult.
Basically, we could do something like:

  bc = T*u_components == expression

Here expression is any UFL-expression, which is compiled into some
future extension of UFC.

u_components will typically be:
  u - all values (in case of non-scalar element)
  u[0] - a single component (vector valued element, must be fixed integer)
  u[0, 1] - a single component (tensor valued element, must be fixed integers)
  dot(u, n) - normal component of u (handled differently depending on element!)

We _could_ define a small set of valid expressions for u_components
like the list above, if we are sure we can cover everything we need.

The extensions we need in UFC will be something to evaluate the
expression above given coefficients (similar to the integral classes
and tabulate_tensor), and something to tabulate which dofs to touch.

How to handle the dot(u, n), I'm not sure?

Regarding the subdomain numbering: we may want to have separate
numbering of T(i) and ds(i), to allow for BCs with some components
Dirichlet and some Neumann? Any examples of this?

Also, in many cases you'll want to just fix one point to avoid translation.
This kind of situation makes it difficult...

> But anyway, this is just a thought that has not been thought properly

Those are the best kind!

--
Martin

---

From: Kent-Andre Mardal <[email-masked]>
Date: Thu, 19 Jun 2008 16:20:29 +0200

> We could collect them more explicitly:

Yes, or even 

Axb = System ( u*v*dx == f*vdx, T*u == g) 

> (As above this syntax isn't possible in the embedded language approach

Sorry, I meant == ..

> like the list above, if we are sure we can cover everything we need.

But aren't these already in UFL. 

> Dirichlet and some Neumann? Any examples of this?

I'm not sure. 

> This kind of situation makes it difficult...

I guess one would the mark one node and have a subdomain which is really
only a node (?). 

> > But anyway, this is just a thought that has not been thought properly

Sometimes, but not allways ...

Kent

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Thu, 19 Jun 2008 17:35:26 +0200

> Axb = System ( u*v*dx == f*vdx, T*u == g)

Even better!

> But aren't these already in UFL.

The syntax for the items in the list above, yes, and we'll reuse that.
But it doesn't make sense in this context to write f.ex.
  T*exp(f)*u == g,
so the list above is a suggestion for what we define as valid in place
of u_components.

> I guess one would the mark one node and have a subdomain which is really

That should work.

I'm convinced that we can at least cover most boundary conditions
with this, but I'd like to see how it's best implemented in the
assembler before we design the UFC extension and incorporate
this into the form compilers.

--
Martin

---

From: [email-masked]
Date: Fri, 20 Jun 2008 11:25:39 +0200 (CEST)

> so the list above is a suggestion for what we define as valid in place

I see.

> this into the form compilers.

I will start playing with the Assembler shortly and will probably get back
to you
next week.

Kent

---

From: Anders Logg <[email-masked]>
Date: Fri, 20 Jun 2008 22:10:10 +0200

If we should add notation for boundary conditions to UFL, then I think
it must also be supported through the UFC interface so that the
assembler may assemble the correct matrix/vector directly (without
zeroing out entries in a separate step).

Regarding the notation, maybe it should be a constraint on the
function spaces. Take for example Poisson with u = g on the boundary:

  Find u in H^1_g : a(v, u) = L(v) for all v in H^1_0

The boundary conditions are part of the test and trial spaces
(subscripts g and 0 on H^1). We could do something like

  V = FunctionSpace("Lagrange", "triangle", 1) # test space
  W = FunctionSpace("Lagrange", "triangle", 1) # trial space
  V.restrict(g)
  W.restrict(0)

  a = dot(grad(v), grad(u))*dx
  L = v*f*dx

-- 
Anders

---

From: [email-masked]
Date: Mon, 23 Jun 2008 11:31:04 +0200 (CEST)

> zeroing out entries in a separate step).

I don't understand what you mean with this. It has to be a separate step
that is performed after the tabulate_tensor functions.

>   a = dot(grad(v), grad(u))*dx

We could do that.

Boundary conditions are usually part of the test and trial space, but this
is just short-hand
notation. Either you can say

where H^1_g are the functions in H^1 with Tu=g on \partial \Omega
(and similar with H^1_0)

or you can say

Find u in H^1 : a(v, u) = L(v) for all v in H^1 with Tu = g on \partial
\Omega
(and similar with v).

It is just that H^1_0 is so established that we hardly think of traces.

But I think

T(i) * u == g

looks good, it can also be done on the test or trial space

V = FunctionSpace("Lagrange", "triangle", 1)
T(i)* V == g

But we can also do restrict (in both cases).

Kent

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Mon, 23 Jun 2008 11:51:41 +0200

> that is performed after the tabulate_tensor functions.

The assembler cannot use UFL directly, and UFL is supposed to pass
through a form compiler. Somebody must convert the UFL to something
the C++ assembler can eat, and UFC is the interface we hide equation
details behind. Unless you have a better option?

> looks good, it can also be done on the test or trial space

I think you've still not answered the question about how this is used.
The line above won't work alone, we'll have to do at least

  bc = T(i) * u == g

which looks a bit more strange than your line alone,
and then 'bc' must be attached to an equation somehow
if we are to pass around simple objects like we do with forms.

> But we can also do restrict (in both cases).

I don't think I like the V.restrict(g) option.
One reason is purely visual, it doesn't quite
"fit" with the rest of UFL. More importantly,
I'd like to use the same element object for
coefficient functions and basis functions.
Maybe we can do:

V = FiniteElement("Lagrange", "triangle", 1)
Vg = Restricted(V, g, i) # subspace of V with u=g on dS_i
u = TrialFunction(Vg)
f = Function(V)

--
Martin

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Mon, 23 Jun 2008 11:54:01 +0200

> details behind. Unless you have a better option?

However, I agree that it should probably be a separate step in the
assembler after calling tabulate_tensor but before calling A.add.

--
Martin

---

From: [email-masked]
Date: Mon, 23 Jun 2008 13:22:13 +0200 (CEST)

> details behind. Unless you have a better option?

I agree that UFL objects cannot be used directly.

> I think you've still not answered the question about how this is used.

I'll get back to that. I have not looked at all details yet. :)

> Vg = Restricted(V, g, i) # subspace of V with u=g on dS_i

I like the abstraction/name trace better than restriction, but no hard
feelings.
It could also be 'Constraints'.

Kent

---

From: Anders Logg <[email-masked]>
Date: Mon, 23 Jun 2008 13:40:21 +0200

> assembler after calling tabulate_tensor but before calling A.add.

Yes, and the important thing is that it's communicated to the
assembler using the UFC interface (which we may need to extend).

It should not be something that is communicated from UFL --> DOLFIN
without UFC.

-- 
Anders

---

From: Anders Logg <[email-masked]>
Date: Mon, 23 Jun 2008 13:58:49 +0200

> It could also be 'Constraints'.

How about this:

A FunctionSpace is created from a FiniteElement and a set of
Constraints, so we may do things like

  element = FiniteElement("Lagrange", "triangle", 1)
  constraints = (Trace(u, s0) == g, Trace(u, s1) == h)

  V = FunctionSpace(element, constraints)

We can let __eq__ in Trace return a Constraint object.

See here also:

  http://www.geuz.org/getdp/doc/slides/getdp-10.html

-- 
Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Mon, 23 Jun 2008 14:09:57 +0200

> We can let __eq__ in Trace return a Constraint object.

Looks good to me. So should a BasisFunction always be defined on a
FunctionSpace instead of an element? Or is "BasisFunction(element)" a
shortcut for "BasisFunction(FunctionSpace(element))" then?

We still need to support components like:

constraint = Trace(u[0], s0) == g
V = FunctionSpace(element, constraint)

I don't think we have to support arbitrary expressions in place of g.

--
Martin