So this turned out to be quite an odyssey, but I was finally able to solve the problem for different cases and I want to share my solution for future reference. The packages needed in addition to JuAFEM.jl are CoherentStructures.jl and Distances.jl. My solution fixes all corner nodes completely, so all degrees of freedom (dofs) at corners are zero.
For problems where there is only 1 dof per node, the functions of CoherentStructures can be used very well. A ctx::gridcontext object has to be created either directly via the grid constructors of the package or indirectly by transferring all necessary objects. Note that gridcontext objects are hardcoded only accepting a DoFHandler with a scalar field named “T”. Creating it indirectly requires an object “PointLocator”. As far as I know, the PointLocator does not work for all element types yet, but for the implementation of only periodic boundary conditions it is not necessary, so you can create a helper object which does nothing and is only required for the gridcontext constructor. Once this is finished you need a predicate function to specify how far away opposite nodes are from each other and then the BoundaryData function will find corresponding degrees of freedom. Using applyBCS will give you the reduced stiffness matrix (PBC are applied to the equation system by adding the coupled rows and columns together), the reduced equation system can be solved and the reduced solution vector can be retransformed by undoBCS.
For problems with more than 1 dof per node but only 1 field, for example solving elastic displacement problems in 3d (ux, uy, uz for each node has to be calculated), some changes have to be made. You still have to create a gridcontext object with a helper DoFHandler T to find opposite DoFs. To calculate the reduced stiffness matrix and also maybe the reduced force vector or residual, I wrote a custom function, which does the additions now for each block instead of single rows and columns.
For problems with several degrees of freedom and also several fields in the dofhandler, it gets very complicated. Let us assume we have a 2d field u and a 1d field p in the dofhandler. First we have to realise, that the close!(dh) function of JuAFEM will perturb our dofs: In the equation system we solve, row 1 does not necessarily correspond to node 1, dof 1. Instead, a new order is chosen and the WriteVTK automatically returns this order to the original, so it is easy to miss. The reason for this step is that we usually want to create a sparsity pattern for our stiffness matrix: It is better to solve a band matrix than any sparse matrix, so from a numerical point of view, this step can be very useful (However this applies only to direct solvers like the backslash-operator “\”, iterative solvers should not be affected). Even worse, not all dofs of a single node are saved consecutively, only those dofs, which are part of the same field. So now our solution vector might look like this: [u5x u5y u8x u8y p5 p8 u3x u3y p3 and so on…]. I will call these “global” gdofs, so gdof 1 is here the first u dof of node 5.
To apply PBCs, I did the following steps:
-
Create a helper array nodedoflist, where the (i,j)-entry will save the gdof number of node i and dof j.
-
Calculate the reduced stiffness matrix and residual by using CoherentStructures.BCTable and nodedoflist. The reduced system has at least ordered all dofs of a single node consecutively.
-
Solve the reduced system with an iterative solver.
-
Split the reduced solution vector into the single dofs: for n dofs, every n-th component is the same of dof of a different node.
-
Use undoBCS of CoherentStructures, to return the reduced sub-solutions to the unreduced sub-solutions.
-
Construct the full solution vector by using the nodedoflist table in reverse: find the corresponding node and dof for all “complete” dofs in the solution vector.
So yeah overall, it took quite a lot time, but I am satisfied with the results. If you have any questions regarding this issue or want to look at my custom written functions, feel free to contact me. I currently dont want to make them completly public, as my solution feels more like a workaround than good programming pratice. I also want to thank the creators of CoherentStructures, as the contact with them was very helpful for me!