[ANN] Gridap.jl: A feature-rich Finite Element ecosystem 100% in Julia

Package Announcements
, , ,
1

We (@santiagobadia and @fverdugo) are happy to announce Gridap.jl a new Finite Element framework written in Julia for the numerical simulation of a wide range of mathematical models governed by partial differential equations (PDEs).

Gridap.jl is designed with three main purposes:

  • to help application experts to easily simulate complex PDE-based real world problems
  • to help researchers improve productivity when developing new Finite Element techniques
  • and also for its usage in numerical PDE courses

What can Gridap.jl do ?

Gridap is able to solve different types of PDEs:

  • linear and nonlinear
  • Steady state, state-dependent, and time-dependent
  • Single field and multi-physics
  • PDEs in 1D, 2D, 3D, 4D, 5D, …

with different types of Finite Element techniques:

  • Continuous and Discontinuous Galerkin methods
  • Grad, div, and curl-conforming interpolations of arbitrary order
  • Embedded Finite Elements (AgFEM and CutFEM)

that can easily be extended.

Why Gridap.jl ?

The main motivation behind Gridap.jl is to find an improved balance between computational performance, user-experience and work-flow productivity, when working with Finite Element libraries. To this end, the library design is based on lazy data-structures that represent objects (e.g., elemental matrices and vectors) on the entire computational domain. This allows us the library developers to hide assembly loops and other core computations from the user-code leading to a very compact, user-friendly, syntax, while providing a high degree of flexibility to users to define their own FE solvers.

For instance a Poisson problem can be solved with Gridap.jl in few lines of code by using the different abstractions provided by the library:

using Gridap
# FE mesh (aka discrete model)
url = "https://raw.githubusercontent.com/gridap/Tutorials/master/models/model.json"
model = DiscreteModelFromFile(download(url,"model.json"))
# Test and trial FE spaces
order = 2; V0 = TestFESpace(
  reffe=:Lagrangian, order=order, valuetype=Float64,
  conformity=:H1, model=model, dirichlet_tags="sides")
g(x) = 2.0; Ug = TrialFESpace(V0,g)
# Volume terms
trian_Ω = Triangulation(model)
degree = 2*order; quad_Ω = CellQuadrature(trian_Ω,degree)
a_Ω(u,v) = ∇(v)⋅∇(u); b_Ω(v) = v
t_Ω = AffineFETerm(a_Ω,b_Ω,trian_Ω,quad_Ω)
# Neumann terms
neumanntags = ["circle", "triangle", "square"]
trian_Γ = BoundaryTriangulation(model,neumanntags)
quad_Γ = CellQuadrature(trian_Γ,degree)
h(x) = 3.0; b_Γ(v) = v*h
t_Γ = FESource(b_Γ,trian_Γ,quad_Γ);
# FE problem and solution
op = AffineFEOperator(Ug,V0,t_Ω,t_Γ)
solver = LinearFESolver(LUSolver())
uh = solve(solver,op)
# Output file for visualization
writevtk(trian_Ω,"results",cellfields=["uh"=>uh])

fig_uh
fig_uh957×406 51.7 KB

How to start ?

If you are further interested in the project, visit the Gridap.jl repository.

If you want to start learning how to solve PDEs with Gridap.jl, then visit our Tutorials repository.

72 Likes
2

Thanks for sharing this great package! Could you please say a word about how Gridap differs from JuliaFEM, FinEtools and JuAFEM?

1 Like
3

Each Finite Element library is unique in some way with its own features. This is like with text editors, you cannot say that X can replace Y.

In any case, Gridap is a quite large project with a decent amount of functionality ready to be used in different scenarios.

Perhaps the distinctive feature is the use of lazy arrays to work with cell-wise data and the compact user API, which is as compact as FEniCs but 100% in julia.

12 Likes
4

A diplomatic answer :smile: Anyways, always good to see more high quality FEA packages! So how easy or difficult would you say it is to implement some basic topology optimization on top of Gridap?

1 Like
5

Really impressive docs!

What about time dependent problems? Integration with DiffEq framework?

Also, if it is similar in functionality to Fenics, then that would be a useful thing to state in the Readme, as Fenics is so widely known.

1 Like
6

So how easy or difficult would you say it is to implement some basic topology optimization on top of Gridap?!

It should be possible. If you are familiar with topology optimization, we can collaborate to write a short tutorial on TO with gridap

7 Likes
7

Sounds good. I can maybe support Gridap in GitHub - JuliaTopOpt/TopOpt.jl: A package for binary and continuous, single and multi-material, truss and continuum, 2D and 3D topology optimization on unstructured meshes using automatic differentiation in Julia..

2 Likes
8

We can do time integration (method of lines)

https://github.com/gridap/GridapODEs.jl

and we are developing some tutorials.

The coupling with DiffEq is ongoing work, we have just some preliminary results.

9 Likes
9

If you specify which are the inputs and outputs you need from a FE solver, I can say how this is done in Gridap.

10

For now, TopOpt is based on JuAFEM. I want to abstract away the FEA backend though to enable supporting other packages. Once I have a minimal interface, let’s have this discussion again, hopefully soon!

3 Likes
11

I’ve actually found this package a bit earlier on this survey site, impressive work!
A few remarks:

  • I noticed Gridap.TensorValues that subtype Number to have broadcasted operations work nicely, have you considered using the julia Ref() on normal arrays instead of introducing a new abstraction? The reason being that many methods that people have implemented on AbstractArrays can then be reused in a pointwise manner.
  • Last time I tried, I remember there was no easy way to create the grid without importing a mesh-file (relying on gmsh). Do you consider domain specification and meshing out of scope for your package or have you considered adding some primitives there? I realize GridapGmsh covers this somewhat but it seems more like a low-level interface to gmsh.
  • How easy is it to hook into the assembly/solve phase? Can one for example implement dof-reorderings or other preconditioners easily?
  • Skimming the codebase, I noticed that the finite element types are hard-coded in quite a few places (namely in dof-mappings). Can you comment on how hard it would be to add a new custom finite element (nodal basis, reference-element-mapped)?
12

Last time I tried, I remember there was no easy way to create the grid without importing a mesh-file (relying on gmsh). Do you consider domain specification and meshing out of scope for your package or have you considered adding some primitives there? I realize GridapGmsh covers this somewhat but it seems more like a low-level interface to gmsh.

Mesh generation is definitively out-of-scope. We only provide a simple Cartesian mesh generator for convenience. Unstructured grids are generated elsewhere. If you have a msh file generated with GMSH, then you can read it with function GmshDiscreteModel provided by the GridapGmsh.jl package (we have tested for linear triangles and tets). The other option is to avoid generating unstructured meshes and use embedded FEM. This can be done with GridapEmbedded.jl. Here, you can specify your geometry as bolean operations between simple shapes.

2 Likes
13

I don’t know if I have understood it correctly, but I think this is not true. Finite element types are not hard-coded. The code is implemented in terms of an abstract type called ReferenceFE which can be extended with custom user-defined concrete implementations. It is true that the library provides a default local-to-global dof map, but you can implement you own if you wish. You just need to implement a new constructor FESpaceWithMyOwnDOFMap that should return some object with a well defined interface.

14

Yes, you can use the linear solver you want. For instance, for the example above you can do:

op = AffineFEOperator(Ug,V0,t_Ω,t_Γ)
A = get_matrix(op)
b = get_vector(op)
x = A\b # Solve this with your favorite linear solver
uh = FEFunction(Ug,x)
writevtk(trian_Ω,"results",cellfields=["uh"=>uh])
15

I noticed that at least one member of this forum found Gridap.jl in the PDE-software list
prior to this announcement!

That is great, and so I would like to invite all of you: include your work in the list, it seems to provide good exposure. Your software may come in handy in someone else’s work. Do let them know of its existence. File an issue, or open a PR.
Thanks,
P

13 Likes
16

Amazing work! Thank you for this great contribution Gridap.jl authors!

I’d like to jump in and share another standardization effort that is ongoing regarding meshes for FEM, spatial modeling and visualization: https://github.com/JuliaGeometry/GeometryBasics.jl/issues/15

A prototype a la Tables.jl is being sketched in Meshes.jl. Ultimately, we would like to benefit from a common interface so that we can operate on mesh objects within different frameworks.

7 Likes
17

This seems like a very nice package with good documentation and tutorials. Thank you.

I’ve run through the tutorial set without trouble. However, when I attempt create my own model.json file, as described in Tutorial 1 notebook as:

“The file "model.json" is a regular json file that includes a set of fields that describe the discrete model. It was generated by using together the GMSH mesh generator and the GridapGmsh package. First, we generate a "model.msh" file with GMSH (which contains a FE mesh and information about user-defined physical boundaries in {GMSH} format). Then, this file is converted to the Gridap-compatible "model.json" file using the conversion tools available in the GridapGmsh package. See the documentation of the GridapGmsh for more information.”

I get stuck at the GridapGmsh step. I don’t see any tools in there (documentation or source code) to perform the conversion. Can you point me in the right direction?

18

Thanks for the interest! Try this:

using Gridap
using Gridap.Io
using GridapGmsh

model = GmshDiscreteModel(“model.msh”)

fn = “model.json”
to_json_file(model,fn)

2 Likes
19

Ah yes thank you. That’s the kick I needed. I see now that it’s in the tutorial/models folder for each of the model .jl files.

20

Nice that it helps!

BTW, you can use the object “model” directly in your computation without writing it to json and reading it again. We only use the json file in the tutorials as a way of avoinding to force the users to install gmsh.