Find boundary elements

You should be able to iterate over the relevant cells and then add the
dof associated with that cell:

   dm = V.dofmap()
   bdr_elements = []
   for cell in cells(mesh):
      if not some_criteria:
         continue
      bdr_elements.extend(dm.cell_dofs(cell.index()))
   bdry_elements np.array(bdr_elements)

The iterator can be over a SubsetIterator too. Thene you can give a
MeshFunction with a certain label, and you will get an iterator over
that subset.

Johan

On 05/21/2012 12:15 PM, K. Hoffmann wrote:
> New question #197941 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/197941
>
> Hi
>
> I have defined a piecewise constant function (in 2D)
>
> ------------------
> f = Function(FunctionSpace(mesh, "DG", 0))
> ------------------
>
> I need this function to attain a certain value on the boundary elements. How do I find the element numbers of the boundary elements, such that I can define the function values on the boundary like this?
>
> ------------------
> f.vector()[bdry_elements] = something
> -------------------
>
> Thanks in advance
>
>
>

On 05/21/2012 01:40 PM, Johan Hake wrote:
> Question #197941 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/197941
>
> Status: Open => Answered
>
> Johan Hake proposed the following answer:
> You should be able to iterate over the relevant cells and then add the
> dof associated with that cell:
>
> dm = V.dofmap()
> bdr_elements = []
> for cell in cells(mesh):
> if not some_criteria:
> continue
> bdr_elements.extend(dm.cell_dofs(cell.index()))
> bdry_elements np.array(bdr_elements)
>
> The iterator can be over a SubsetIterator too. Thene you can give a
> MeshFunction with a certain label, and you will get an iterator over
> that subset.

If you also make:

   bdr_elements

a set, this algorithm should be able to be generalized to other elements
too, where some dofs are shared among cells. The syntax to extend is
different but that is trivial.

Johan

> Johan
>
> On 05/21/2012 12:15 PM, K. Hoffmann wrote:
>> New question #197941 on DOLFIN:
>> https://answers.launchpad.net/dolfin/+question/197941
>>
>> Hi
>>
>> I have defined a piecewise constant function (in 2D)
>>
>> ------------------
>> f = Function(FunctionSpace(mesh, "DG", 0))
>> ------------------
>>
>> I need this function to attain a certain value on the boundary elements. How do I find the element numbers of the boundary elements, such that I can define the function values on the boundary like this?
>>
>> ------------------
>> f.vector()[bdry_elements] = something
>> -------------------
>>
>> Thanks in advance
>>
>>
>>
>

Maybe I should have noted that I'm quite new to dolfin, so your code seems a bit complex to me. Anyway, I think I get the main idea of your suggestions Johan.

I'm still puzzled about which criteria to use. The criteria could use the dof of each element, but how can I find this? It seems to me that cell_dofs(N) returns the nodes of element N. Is there a way to get the dof of these nodes then?

On 05/21/2012 04:25 PM, K. Hoffmann wrote:
> Question #197941 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/197941
>
> Status: Answered => Open
>
> K. Hoffmann is still having a problem: Maybe I should have noted that
> I'm quite new to dolfin, so your code seems a bit complex to me.
> Anyway, I think I get the main idea of your suggestions Johan.

The main idea is that the information that maps geometry to doegrees of
freedoms, the indices of the vector, is kept within the FunctionSpace.
More specific within the DofMap object. You use the dof map object to go
from cell to global degree of freedom.

> I'm still puzzled about which criteria to use. The criteria could
> use the dof of each element, but how can I find this?

Not sure what you mean with criteria.

> It seems to me
> that cell_dofs(N) returns the nodes of element N. Is there a way to
> get the dof of these nodes then?

In a DG0 element there are only 1 dof per element, which for serial runs
happens to be the same as the cell index.

Johan

I never got it to work using dof. Instead I now use the coordinates of the vertices of each element (which are found using tabulate_coordinates) to do the determination. This is not very elegant, but fairly easy to implement since I use a circular geometry. Thanks for the help Johan.

On 05/24/2012 09:50 AM, K. Hoffmann wrote:
> Question #197941 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/197941
>
> Status: Answered => Solved
>
> K. Hoffmann confirmed that the question is solved:
> I never got it to work using dof. Instead I now use the coordinates of
> the vertices of each element (which are found using
> tabulate_coordinates) to do the determination. This is not very elegant,
> but fairly easy to implement since I use a circular geometry. Thanks for
> the help Johan.
>

Ok, good you found a way. Your problem might relate to not being able to
find the cells on the boundary. Also note that DG0 elements does not
have any dofs at the boundary.

This snippet tabulate the dofs at a certain boundary, given by a
FacetFunction. It is a bit tricky as one need to go from face to cell to
face again, to find the dofs.

Johan

   import numpy as np
   from dolfin import *
   # Assume the existence of face_subdomains

   dofs = set()
   # Tabulate the local map
   local_dofs = np.zeros(3, dtype=np.uintc)
   facet_to_cell_dofs = []
   for local_facet in range(4):
     dm.tabulate_facet_dofs(local_dofs, local_facet)
     facet_to_cell_dofs.append(list(local_dofs))
   facet_to_cell_dofs = np.array(facet_to_cell_dofs)

   # Iterate over relevant faces
   for face in SubsetIterator(face_subdomains, 2):
     # Get cell only one as the face is a boundary face
     cell = Cell(mesh, face.entities(3)[0])
     # Get local face index
     local_face_index = (cell.entities(2) == \
                         face.index()).nonzero()[0][0]
     # Get cell dofs and pick the correct slice of these
     dofs.update(dm.cell_dofs(cell.index())\
                 [facet_to_cell_dofs[local_face_index]])

   dofs = np.array(dofs)