[firedrake] Enforcing slip boundary condition strongly

[Lawrence Mitchell]

> On 10 Jul 2015, at 10:25, David Ham <David.Ham at imperial.ac.uk> wrote:
> 
> Hi Colin,
> 
> That does not appear to be what Justin is asking for. He is asking how to set a Dirichlet condition on one component of a VectorFunctionSpace. I'm not sure we can do that right now, but I think this should be fixable. I want to check what the Dolfin syntax for this is (so we maintain compatibility) and then try to fix it. (Of course someone might want to pleasantly surprise me and tell me we can do this already!)

We can't do it at the moment, the rest of this email explains why and sketches a possible solution.

OK, so there are two aspects to this, which make things somewhat tricky.

Our data model for indexing vector fields is that we store the Map to the node and stack the dofs on top.

To apply such slip BCs, where we only set one of the two dofs, we need to do the following:

- Apply the BC value to the appropriate part of the vector.  I think this is reasonably straightforward and appear to have something that works.

- Fix up the operator appropriately.  This latter option is trickier.  The way we currently do this is to drop boundary conditions on the floor when assembling the operator.  This is done by flagging the appropriate rows/columns to be dropped by putting a negative value in the Map at the appropriate place and letting PETSc do the rest.

For vector-valued fields we do block insertion, i.e. we insert a 2x2 block at each location in the matrix indexed by the base map.  So the current setup can only drop complete blocks, rather than part of one (as required in this case).

The question is, then, how to communicate the relevant bit of information into the code that inserts.

How about this:

When we splat the values in the Map we convert:

val => -(val + 1)

If we're apply the BC on part of the map we decorate it with the index offset.

We then build a local (unrolled map) and fix up the entries.  The negative offset encoding allows us to recover the correct value.

So assembly turns from:

A = build_element_tensor(...);
MatSetValuesBlockedLocal(block_rows, block_cols, A);

Into:

A = build_element_tensor(...);
{
    PetscInt rows[arity*dim];
    PetscInt cols[arity*dim];
    // Assume we're setting the bc on the second component.
    for ( i = 0; i < arity; i++ ) {
         if (row_map[i] < 0) {
             rows[i*dim + 0] = dim*(-(row_map[i] + 1));
             rows[i*dim + 1] = -1;
         } else {
              rows[i*dim + 0] = dim*row_map[i];
              rows[i*dim + 1] = dim*row_map[i] + 1;
         }
    }
    // Same for cols
    MatSetValuesLocal(rows, cols, A);
}

All the other options I came up with seem more complicated.

Lawrence
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 455 bytes
Desc: Message signed with OpenPGP using GPGMail
URL: <http://mailman.ic.ac.uk/pipermail/firedrake/attachments/20150710/a71c8d2a/attachment.sig>