On Sat, Nov 3, 2012 at 8:28 PM, Mikael Mortensen <mikael.morten...@gmail.com> wrote: > 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:-) >
Great. > 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 would eliminate the dofs in the DofMap. We can make the global_dimension member data rather than going through the ufc dofmap as it does at present. That way the global_dimension can be set by the DofMap class. Garth > 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 > > > _______________________________________________ > Mailing list: https://launchpad.net/~dolfin > Post to : dolfin@lists.launchpad.net > Unsubscribe : https://launchpad.net/~dolfin > More help : https://help.launchpad.net/ListHelp _______________________________________________ Mailing list: https://launchpad.net/~dolfin Post to : dolfin@lists.launchpad.net Unsubscribe : https://launchpad.net/~dolfin More help : https://help.launchpad.net/ListHelp
