Skip to content

Instantly share code, notes, and snippets.

@AbdAlazezAhmed
Last active December 7, 2023 10:12
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save AbdAlazezAhmed/2a75a6f993ffacadeb4b277d7254d222 to your computer and use it in GitHub Desktop.
Save AbdAlazezAhmed/2a75a6f993ffacadeb4b277d7254d222 to your computer and use it in GitHub Desktop.
GSoC23: Discontinuous Galerkin Infrastructure For the finite element toolbox Ferrite.jl

About the project

Project: Discontinuous Galerkin Infrastructure For the finite element toolbox Ferrite.jl

Student: Abdulaziz Hamid

Mentor: Dennis Ogiermann & Fredrik Ekre

Organization: The Julia Language

Description: This project implements the infrastructure required to use discountinuous Galerkin method in the finite element toolbox Ferrite.jl.

Pull requests and their corresponding state

# Title State
729 ConstraintHandler for discontinuous interpolations merged
710 Cross-elements coupling for DiscontinuousLagrange sparsity patterns. merged
747 InterfaceCache & InterfaceIterator merged
751 Make benchmarks work with master branch merged
743 InterfaceValues for DG interface integration merged
779 Implement FaceQuadratureRule for RefPrism and RefPyramid merged
790 Arbitrary order Lagrange interpolations for hypercubes (with arbitrary basis) and triangle (equidistant basis only) pending
787 DG heat equation tutorial pending

Implementation details

Working with constraints

Ferrite uses ConstraintHandler to strongly enforce Dirichlet, affine and periodic dirichlet boundary conditions. Dirichlet boundary condition sets the nodal values at the specified nodes. This can be enforced strongly by modifying the stiffness matrix or weakly by using penalty terms.

ConstraintHandler didn't have proper support for DiscontinuousLagrange interpolations as the facedof_indices,edgedof_indices and vertexdof_indices aren't dispatched for discontinuous interpolation as all the degrees of freedom (dof) are moved to the cell interior.

The solution was creating new methods dirichlet_*dof_indices which default to *dof_indices for continuous interpolation and using it instead in ConstraintHandler.

PR #729: ConstraintHandler for discontinuous interpolations (merged)

Sparsity patterns

create_sparsity_pattern creates a SparseMatrixCSC with stored values according to the dof numbering provided. This takes into account the coupling between fields inside the element but not cross-element coupling. However, when using DG methods the interface integrals result in cross-element coupling.

This was solved by modifying the functions to accept optional topology and cross_coupling arguments to add the cross-element coupling entries to the sparse matrix.

PR #710: Cross-elements coupling for DiscontinuousLagrange sparsity patterns. (merged)

Iterators

Ferrite provides CellIterator and CellCache which are used to conveniently iterate over all, or a subset, of the cells in a grid. However, iterating over interfaces had to be done using two cell iterators and a neighborhood check, which isn't convinient nor effecient.

This was solved by introducting InterfaceIterator which iterates over all the interfaces of the grid effeciently by utilizing the topology.

PR #747: InterfaceCache & InterfaceIterator (merged)

Benchmarks

Ferrite has benchmarking scripts to analyze its performance and help tracking down performance regression issues. These benchmarks were developed for continuous interpolations so adding support for DG methods require adding their benchmarks. However, these scripts were not updated to work with the master branch.

The scripts were updated to work with master branch in PR #751: Make benchmarks work with master branch (merged)

The benchmarks for discontinuous interpolations were added in PR #743: InterfaceValues for DG interface integration (pending)

Integration

Ferrite chaches shape function values and gradients at the quadrature points to optimize assembly performance using FEValues structs. Doing interface integrals involves using values from two different cells and has methods such as jumps and averages which are not present in FaceValues. Thus, implementing InterfaceValues was needed.

Also, when doing interface integrals the quadrature points must be synced to the same spatial coordinates on both sides of the interface, which is not the case by defauls. Thus, reinitializing InterfaceValues should transform the quadrature points from one side of the interface to another either using permutations or a transformation matrix. Using permutations should have better performance but only the hypercubes quadrature points are guaranteed to have a predetermined permutation, as simplices use tables not quadrature point generators.

Quadrature points transformations were implemented in PR #779: Implement FaceQuadratureRule for RefPrism and RefPyramid (merged)

InterfaceValues was implemented with quadrature points transformations in PR #743: InterfaceValues for DG interface integration (pending)

Heat equation example

paraview screenshot of the solution of heat equation over a cube with Dirichlet boundary condition

Part of the project was adding an example of using the interior penalty method in the heat equation example. The reason for choosing this example is its simplicity, as it's considered the hello-world of FEM.

The weak form solved in the example is

$$\int_\Omega \nabla u \cdot \nabla \delta u dΩ - \int_\Gamma [[u]] \cdot \{\nabla \delta u\} + [[\delta u]] \cdot \{\nabla u\} dΓ + \int_\Gamma \mu [[u]] ⋅ [[\delta u]] dΓ = \int_\Omega \delta u dΩ,\\$$

where $\delta u$ is the test function.

The example is in PR #787: DG heat equation tutorial (pending)

Arbitrary order interpolations

Another part of the project is utilizing sum-factorization to optimize assembly for higher order elements. This requires implementing arbitrary order interpolation as a first step.

Arbitrary order Lagrange and discontinuous Lagrange interpolations were implemented in PR #790: Arbitrary order Lagrange interpolations for hypercubes (with arbitrary basis) and triangle (equidistant basis only) (pending)

Face skeleton construction

Most of the DG related methods utilize topology. Thus, optimizing topology construction can affect the performance of the overall routine.

PR #759: FaceSkeleton Vroom 1 (pending) speeds up the construction of FaceSkeleton by utilizing neighborhood information.

Future work

  • Implement sum factorization for higher order integration as in Issue #389

Acknowledgement

I gratefully acknowledge the guidance of my mentors Fredrik Ekre and Dennis Ogiermann, and Ferrite members for always helping me throughout the project. Thanks ❤️

Links

Ferrite source code: https://github.com/Ferrite-FEM/Ferrite.jl

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment