restricted dof

Asked by dbeacham on 2010-02-04

I think I'm right in thinking that dolfin can only perform n-dimensional calculations on meshes embedded in R^n, however is there an easy way to restrict the dof/manipulate the assembly such that when I do something like

a = u*v*ds
L = f*v*ds

the assembled forms have non-zero determinant? Currently the assembly gives these matrices (reordered):
      [ds stuff | 0 ]
A= [------------------]
      [ 0 | 0 ]
      [ds stuff ]
b= [----------]
      [ 0 ]

whereas I'd either like to just keep the ds stuff (although that would require another mapping between the reordered ds nodes and the original mesh) or possibly set the bottom right of A to the identity matrix (without having to manipulate the matrix myself).

Alternatively, are there any better approaches I could adopt?

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
dbeacham
Solved:
2010-02-11
Last query:
2010-02-11
Last reply:
2010-02-09
Anders Logg (logg) said : #1

On Thu, Feb 04, 2010 at 12:08:58PM -0000, dbeacham wrote:
> New question #99875 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/99875
>
> I think I'm right in thinking that dolfin can only perform n-dimensional calculations on meshes embedded in R^n, however is there an easy way to restrict the dof/manipulate the assembly such that when I do something like
>
> a = u*v*ds
> L = f*v*ds
>
> the assembled forms have non-zero determinant? Currently the assembly gives these matrices (reordered):
> [ds stuff | 0 ]
> A= [------------------]
> [ 0 | 0 ]
> [ds stuff ]
> b= [----------]
> [ 0 ]
>
> whereas I'd either like to just keep the ds stuff (although that would require another mapping between the reordered ds nodes and the original mesh) or possibly set the bottom right of A to the identity matrix (without having to manipulate the matrix myself).
>
> Alternatively, are there any better approaches I could adopt?

I don't know of a way other than to assemble the matrices and then
extract the relevant portions into a new matrix.

We have had some support earlier for restriction of the global set of
dofs to a subset but that has been broken for a while.

--
Anders

dbeacham (blackcabbage) said : #2

Thanks for this - but, could someone give me a quick pointer as to how to move around the dofmap. Using meshfn.values()/Iterators I can easily find the edge instances I wish to integrate over, but would like to know how to convert this knowledge to the indices of the dofs I am integrating over, so I can extract the relevant submatrix for computation. [Sorry if my terminology is a bit confused.]

Apart from trying to do extract the ds section of a matrix for pure boundary integrals, I also want to be able to set the pressure boundary condition at a single point in a Navier-Stokes calculation, so I would also be interested to know how the dofmapping scales with general mixed element spaces (ie Taylor-Hood (P2/P1)). [I assume there is no easier way of setting a single pressure value than this?]

Thanks, David.

Anders Logg (logg) said : #3

On Tue, Feb 09, 2010 at 09:38:37AM -0000, dbeacham wrote:
> Question #99875 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/99875
>
> Status: Answered => Open
>
> dbeacham is still having a problem:
> Thanks for this - but, could someone give me a quick pointer as to how
> to move around the dofmap. Using meshfn.values()/Iterators I can easily
> find the edge instances I wish to integrate over, but would like to know
> how to convert this knowledge to the indices of the dofs I am
> integrating over, so I can extract the relevant submatrix for
> computation. [Sorry if my terminology is a bit confused.]
>
> Apart from trying to do extract the ds section of a matrix for pure
> boundary integrals, I also want to be able to set the pressure boundary
> condition at a single point in a Navier-Stokes calculation, so I would
> also be interested to know how the dofmapping scales with general mixed
> element spaces (ie Taylor-Hood (P2/P1)). [I assume there is no easier
> way of setting a single pressure value than this?]
>
> Thanks, David.

To get the dofs on a facet, you need to call the following functions
in DofMap:

  tabulate_dofs: to get the local-to-global mapping on a cell
  num_facet_dofs: to get the number of dofs on each facet
  tabulate_facet_dofs: to get the local-to-local mapping from the facet
                       dofs to the corresponding cell dofs

For more detailed information, take a look at the documentation for
these functions in the UFC manual.

--
Anders

dbeacham (blackcabbage) said : #4

Thanks, this had put me on the right track. Just one more thing (hopefully):

In the python interface, how do I generate the 'ufc::cell const &' argument used in DofMap.local_dimension? I've tried using V.cell()/V.ufl_element() but they don't work, so I assume they're the incorrect way.

Johan Hake (johan-hake) said : #5

On Thursday 11 February 2010 06:31:18 dbeacham wrote:
> Question #99875 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/99875
>
> Status: Answered => Solved
>
> dbeacham confirmed that the question is solved:
> Thanks, this had put me on the right track. Just one more thing
> (hopefully):
>
> In the python interface, how do I generate the 'ufc::cell const &'
> argument used in DofMap.local_dimension? I've tried using
> V.cell()/V.ufl_element() but they don't work, so I assume they're the
> incorrect way.

I guess you can't as it is not possible to create an ufc::cell in the Python
interface (as far as I know...). I am also not sure we want to expose more of
the UFC interface to PyDOLFIN, and particular not ufc::cell, as it comes with
some nasty double ** attributes in its public interface.

However, in DofMap many functions come in one DOLFIN and one UFC version, but
local_dimension is not one of them. Is there a reason for this Anders? It
should be straight forward to just add it.

Johan

Anders Logg (logg) said : #6

On Thu, Feb 11, 2010 at 08:32:17AM -0800, Johan Hake wrote:
> On Thursday 11 February 2010 06:31:18 dbeacham wrote:
> > Question #99875 on DOLFIN changed:
> > https://answers.launchpad.net/dolfin/+question/99875
> >
> > Status: Answered => Solved
> >
> > dbeacham confirmed that the question is solved:
> > Thanks, this had put me on the right track. Just one more thing
> > (hopefully):
> >
> > In the python interface, how do I generate the 'ufc::cell const &'
> > argument used in DofMap.local_dimension? I've tried using
> > V.cell()/V.ufl_element() but they don't work, so I assume they're the
> > incorrect way.
>
> I guess you can't as it is not possible to create an ufc::cell in the Python
> interface (as far as I know...). I am also not sure we want to expose more of
> the UFC interface to PyDOLFIN, and particular not ufc::cell, as it comes with
> some nasty double ** attributes in its public interface.
>
> However, in DofMap many functions come in one DOLFIN and one UFC version, but
> local_dimension is not one of them. Is there a reason for this Anders? It
> should be straight forward to just add it.
>
> Johan

It is actually possible to call some of them. I have added
dolfin::Cell versions for some of UFC functions, for example
tabulate_dofs in DofMap:

  /// Tabulate the local-to-global mapping of dofs on a cell (UFC cell version)
  void tabulate_dofs(uint* dofs, const ufc::cell& ufc_cell, uint cell_index) const;

  /// Tabulate the local-to-global mapping of dofs on a cell (DOLFIN cell version)
  void tabulate_dofs(uint* dofs, const Cell& cell) const;

I don't think it would be a problem to add more of these. It's useful
for testing out algorithms in Python. See for example the attached
script which makes heavy use of low-level UFC functions from Python.
The algorithm in question has now been ported to C++ in the DOLFIN
Extrapolation class (which was something like a factor 1000 faster).

But in your case, are you sure you need to call DofMap::local_dimension? If the
dimension does not vary over the mesh, you can call
DofMap::max_local_dimension.

--
Anders

dbeacham (blackcabbage) said : #7

You're right, I can just call max_local_dimension, and I had already been using the dolfin version of tabulate_X.

Thanks for all your help - I think I now understand everything I needed to, and have been able to manipulate the assembled matrices to set point values of functions etc. My ds integrals should now just be a (relatively straightforward) variation on the theme.