dolfin team mailing list archive

[Mikael Mortensen <mikael.mortensen@xxxxxxxxx>]

Hi,

This week I have been working on migrating the periodic boundary conditions to the dofmap. I have made it to the point where I now have working code, at least for CG elements, for single and multiple periodic directions, regular and mixed functionspaces, 2D/3D and in parallel. Symmetric problems remain symmetric also for periodic domains. It is actually not all that difficult, but of coarse there are many ways this could be done. My suggested approach is implemented in the branch lp:~mikael-mortensen/dolfin/periodic_to_dofmap. The code is still a bit messy and in need of refinement and optimization. However, it works, which makes it a good starting point:-)

The implemented approach is basically done in 2 steps as follows:

Step 1) Create a facet-to-facet map linking a facet on one periodic boundary to a facet on the other. Store also the processors the facets live on.

The way it's currently implemented is as follows: take a periodic SubDomain (just like before with a map) and add this to the mesh:

periodicX = PeriodicSubDomainX()
periodicY = PeriodicSubDomainY()
mesh.add_periodic_direction(periodicX)
mesh.add_periodic_direction(periodicY)

or

mf = FacetFunction(mesh)
mf.set_all(0)
periodicX.mark(mf, 1)
periodicY.mark(mf, 2)
mesh.add_periodic_direction(mf, 1, 2)

or something like that. Only the first is currently implemented.

Adding the periodic SubDomain directly to the mesh I think has some advantages:

i) The facet-to-facet map will be common for all FunctionSpaces/dofmaps and can be computed just once. 

ii) I think it should be possible(?) do some work on mesh connectivity such that the mesh truly bites its tail. A facet on a periodic boundary would then have two adjacent cells, just like any other internal facet. The periodic boundary would as such no longer be part of the exterior domain. I may be wrong, but I think this would be important for DG methods. I have not done anything in detail here yet though.

After adding the periodicity to the mesh, the user does not have to do anything more about the periodic boundary conditions . Everything will be handled automatically by modifying the dofmaps.

Step 2) Elimination of slaves from the dofmap._dofmap. 

When a FunctionSpace now creates a dofmap on a periodic mesh everything run as normal at first. However, at the end of DofmapBuilder::build I run through the facet-to-facet map and efficiently creates a dof pair map of masters and slaves. This map is subsequently used to exchange all slaves in dofmap._dofmap with masters (and to update off_process_owners).

The major downside at this point is that we are left with a bunch of slave dofs that will never enter computations. To make it work at this point you have to ident all slave-rows from an assembled coefficient matrix (I have at the moment a call to ident_zeros in case of a periodic mesh placed at the end of Assembler::assemble). Nothing has do be done with assembled vectors or scalars. After assembly all the slave rows are empty because the slaves do not appear in the dofmap._dofmap. After identing and solving all slaves are zero, but this does not seem to affect anything like plotting. The only problem I have encountered is when using a function like normalize that performs work on the entire vector. For that purpose I have created a method "update_slaves", that copies the value of masters to slaves. This method must be called before the regular normalize.

I'm currently trying to figure out how (if possible) to get rid of these slave dofs entirely. I would very much appreciate if someone had any thoughts on this? I probably have to subclass DofMap since then at least the global_dimension will be different from the ufc_dofmap. 

I've added a new demo in my branch under demo/undocumented/periodic/python that uses this new approach if somebody cares to try it out.

Best regards

Mikael