gridap / gridap.jl

Grid-based approximation of partial differential equations in Julia

Documentation
Build Status
Community
Citation

What

Gridap provides a set of tools for the grid-based approximation of partial differential equations (PDEs) written in the Julia programming language. The library currently supports linear and nonlinear PDE systems for scalar and vector fields, single and multi-field problems, conforming and nonconforming finite element (FE) discretizations, on structured and unstructured meshes of simplices and n-cubes. It also provides methods for time integration. Gridap is extensible and modular. One can implement new FE spaces, new reference elements, use external mesh generators, linear solvers, post-processing tools, etc. See, e.g., the list of available Gridap plugins.

Gridap has a very expressive API allowing one to solve complex PDEs with very few lines of code. The user can write the underlying weak form with a syntax almost 1:1 to the mathematical notation, and Gridap generates an efficient FE assembly loop automatically by leveraging the Julia JIT compiler. For instance, the weak form for an interior penalty DG method for the Poisson equation can be simply specified as:

a(u,v) =
  ∫( ∇(v)⋅∇(u) )*dΩ +
  ∫( (γ/h)*v*u - v*(n_Γ⋅∇(u)) - (n_Γ⋅∇(v))*u )*dΓ +
  ∫(
    (γ/h)*jump(v*n_Λ)⋅jump(u*n_Λ) -
    jump(v*n_Λ)⋅mean(∇(u)) -
    mean(∇(v))⋅jump(u*n_Λ)
    )*dΛ

l(v) =
  ∫( v*f )*dΩ +
  ∫( (γ/h)*v*u - (n_Γ⋅∇(v))*u )*dΓ

See the complete code here. As an example for multi-field PDEs, this is how the weak form for the Stokes equation with Neumann boundary conditions can be specified:

a((u,p),(v,q)) =
  ∫( ∇(v)⊙∇(u) - (∇⋅v)*p + q*(∇⋅u) )*dΩ

l((v,q)) =
  ∫( v⋅f + q*g )*dΩ +
  ∫( v⋅(n_Γ⋅∇u) - (n_Γ⋅v)*p )*dΓ

See the complete code here.

Documentation

• STABLE — Documentation for the most recently tagged version of Gridap.jl.
• DEVEL — Documentation for the in-development version of Gridap.

Tutorials

A hands-on user-guide to the library is available as a set of tutorials. They are available as Jupyter notebooks and html pages.

Installation

Gridap is a registered package in the official Julia package registry. Thus, the installation of Gridap is straight forward using the Julia's package manager. Open the Julia REPL, type ] to enter package mode, and install as follows

pkg> add Gridap

Plugins

• GridapDistributed Distributed-memory extension of Gridap.
• GridapEmbedded Embedded finite elements in Julia.
• GridapGmsh Generate a FE mesh with GMSH and use it in Gridap.
• GridapMakie Makie plotting recipes for Gridap.
• GridapPardiso Use the Intel Pardiso MKL direct sparse solver in Gridap.
• GridapPETSc Use PETSc linear and nonlinear solvers in Gridap.

Examples

These are some popular PDEs solved with the Gridap library. Examples taken from the Gridap Tutorials.

Poisson equation	Linear elasticity	Hyper-elasticity	p-Laplacian
Poisson eq. with DG	Darcy eq. with RT	Incompressible Navier-Stokes	Isotropic damage

Known issues

Since Julia 1.6 onwards we have noticed large first call latencies of Gridap.jl codes with the default compiler optimization level (i.e., -O2). In general, while developing code, but specially if you are noting high first call latencies, we recommend to run julia with the -O1 flag. For production runs use -O2 or -O3.

Gridap community

You can ask questions and interact with the Gridap community on the Julia Slack channel #gridap (see here how to join). or our gitter.

Contributing to Gridap

Gridap is a collaborative project open to contributions. If you want to contribute, please take into account:

• Before opening a PR with a significant contribution, contact the project administrators, e.g., by writing a message in our gitter chat or by opening an issue describing what you are willing to implement. Wait for feed-back.
• Carefully read and follow the instructions in the CONTRIBUTING.md file.
• Carefully read and follow the instructions in the CODE_OF_CONDUCT.md file.
• Open a PR with your contribution.

Want to help? We have a number of issues waiting for help. You can start contributing to the Gridap project by solving some of those issues.

How to cite Gridap

In order to give credit to the Gridap contributors, we simply ask you to cite the references below in any publication in which you have made use of the Gridap project. If you are using other Gridap sub-packages, please cite them as indicated in their repositories.

@article{Badia2020,
  doi = {10.21105/joss.02520},
  url = {https://doi.org/10.21105/joss.02520},
  year = {2020},
  publisher = {The Open Journal},
  volume = {5},
  number = {52},
  pages = {2520},
  author = {Santiago Badia and Francesc Verdugo},
  title = {Gridap: An extensible Finite Element toolbox in Julia},
  journal = {Journal of Open Source Software}
}

@article{Verdugo2022,
  doi = {10.1016/j.cpc.2022.108341},
  url = {https://doi.org/10.1016/j.cpc.2022.108341},
  year = {2022},
  month = jul,
  publisher = {Elsevier {BV}},
  volume = {276},
  pages = {108341},
  author = {Francesc Verdugo and Santiago Badia},
  title = {The software design of Gridap: A Finite Element package based on the Julia {JIT} compiler},
  journal = {Computer Physics Communications}
}

Contact

Please, contact the project administrators, Santiago Badia, Francesc Verdugo, and Alberto F. Martin for further questions about licenses and terms of use.

gridap.jl's Issues

Renaming Numa

Issue to find the final name for the project

Dof ordering for vector/tensor-valued fields

This issue is only relevant for vector/tensor-valued fields

Currently, we are using different orderings for the dofs in different places of the code:

1. LagrangianDOFBasis uses "node major" ordering (i.e., first component for all nodes, second component for all nodes, etc.)
2. CLagrangianFESpace and DLagrangianFESpace use "component major" ordering. (i.e., all components for first node, then all components for second node, etc).
3. ConformingFESpace: I am not sure what ordering is being used. To ckeck.

Perhaps it is a good idea to unify things and use the same criterion everywhere.

I would choose "component major" since it is quite standard, and some AMG solvers require that the DOFs are ordered in this way.

I am wandering, why LagrangianDOFBasis uses "node major" ordering. Is there any reason?

What to do with `cellsize`?

We have three options:

1. Remove function cellsize from CellValue and CellArray since now the CachedArray is re-sizable, and therefore cellsize it is not longer needed (i.e. we don't need to preallocate to a maximum size since we can resize on the fly)

2. Keep it and cellsize means the true maximum size. For some concrete implementations it would be difficult / expensive to answer this query. (now we follow this approach in the code)

3. Keep it and cellsize means an upper bound of the maxmum size. This is less restrictive, but one has to be careful when computing cellsize of CellArrayFromBinaryBroadCastOperation

I am for option 1, but I am not sure if cellsize can be useful in some context in the future

Geometry clean up

Place to keep comments that will arise during this work

Add support for 1.0

Currently, Gridap requires Julia 1.1. However, only one line (see below) really neends Julia 1.1. We could workaround this to enable Julia 1.0.

1.0 is nicer than 1.1 since the former is a LTS version.

Gridap.jl/src/Geometry/CartesianGeometry.jl

Line 289 in edf618f

cellpointgids = cellpointlids .+ cellgid

Polytope for wedge

Is polytope giving the right nfaces for wedges?

Refactoring of CellValues and Fields

Place to keep miscellaneous comments about the ongoing re-factoring of CellValues and Co.

The re-factoring is taking place in the branch refactoring_cell_values in folders src/REFACTORED and test/REFACTORED

Some notes by @santiagobadia to be discussed

CODE COMMENTS:

• Can we do a trick similar to what is being done for gradient, for the derivative of an operator with respect to the variable, i.e.,

import Gridap: D
a(u,v) = ...
da(u,v,du) = ...
D(::typeof(a)) = da

Here we are assuming that D is in this case d()/du, where u is always the unknown of our FE problem. We must probably create a function that only depends on one unknown.

a(u,v,w,...) = ...

Given v*,w*,...

aa(u) = a(u,v*,w*,...)

and know

D(aa) is the well-defined total derivative.

Does it have sense in our framework? I have to think about it.

• We could now replace the hand-written definition of D for a given a with its automatic differentiation.

• Which are the pre-requisites (rigidity) in the nonlinear FE expressions? What should be respected when defining the nonlinear terms, etc?

• We should keep a consistent way to define derivatives. E.g.,

D(a(u)) = Da(u,du)

i.e., the direction always the last one, whereas the point in which evaluated the first one.

• There are many ways to define FE terms. E.g., AffineFETerm. Is that really needed? Can we define a one-way implementation, less noisy. And probably unify linear and nonlinear terms, using the previous trick. Is there a performance motivation?

• Don't you think that the multivalue wrapper complicates things? E.g., when creating a cell array with an integer per cell, when we index the array what we get is not the integer but a multivalue. Is that correct? Then, the user must extract the integer from the multivalue wrapper.

• What is the motivation for C and DFiniteElementSpace? What is the difference with respect to ConformingFESpace?

FURTHER DEVELOPMENTS:

• Multi-field solvers, e.g., Navier-Stokes with Qp / Q-p-1
• Interface-coupling, e.g., FSI interaction
• Time-depedent problem with time-stepping (all in driver)
• High order tets
• Non-nodal FEs (RT, Nedelec)
• Hybridized FEs (HDG, ...)

DOCUMENTATION/TUTORIALS:

• Create tutorials w/ Jupyter notebooks
• First, high level FE solvers, using the ones in tests
• Second, low level tutorials about rationale of types in Gridap, e.g., cell-wise structs.
• We can also create folders with markdown text that it is automatically uploaded at Github
• Last step, generate documentation with built-in Julia tools, check it compiles, etc

PRESENTATIONS AND PAPERS:

Proceeding about FEMPAR and Gridap, focused on Gridap for Canberra meeting. Proceedings must be sent by mid-August.

CLEANING:

Which parts of the code are dirty and require some maintenance?

• Polytope and node array structures. Eliminate old node array and use always the new one, now missing the new one for heterogenous order in space dims. Clarify the APIs of polytope, etc.

Change Signature of LagrangianRefFE

Adopt following signature

LagrangianRefFE(::Type{T},p::Polytope{D},orders) where {D,T}

Future developments

Interpolator (Combine field, dofbasis, geomap)
FE Space (global numbering, Dirichlet boundary conditions, multiphysics)
Mesh (Structured mesh, Dirichlet boundary, queries, ... GMesh.jl)
FE Operator
Nonlinear operators with Automatic Differentiation
Assembler
VTK IO
Linear algebra (direct solver, CG, PETSc.jl)

Code organization

Discussion to find a proper way of organizing code

Avoid low granularity array views

Justification of why we have to work with views carefully

Add support for Neuman and Robin boundary conditions

This basic functionality is still to be implemented

Work being done in branch surface_integration

Implement second order derivatives (Implement Laplacian)

Required for VMS, Explicit error estimators etc

Problem with fun that provides Int used to create a Map

We have a problem if we define a function that returns an Int and then we use it to create a map MapFromBroadcastUnaryOp.

fun(x::Point{2}) = 1
reffe = fesp.reffe
dofb = reffe.dofbasis
trian = fesp.trian
phi = geomap(trian)
uphys = fun ∘ phi

We have to fix it.

Homogenize node_and_comp_to_dof

The types LagrangianDOFBasis and CLagrangianFESpace uses different strategies to represent the relation between nodes, components, and dofs (see the field node_and_comp_to_dof in both structs).

It would be nice to homogenize things. I would adopt the strategy of CLagrangianFESpace

Dupplicated Polytopes.jl and PolytopesTests.jl files

@santiagobadia I have just realized that we have duplicated files:
src/RefFEs/Polytopes.jl
src/Geometry/Polytopes.jl
on the one hand, and
test/RefFEsTests/PolytopesTests.jl
test/PolytopesTests.jl
on the other hand.

The ones that we want to keep are:
src/RefFEs/Polytopes.jl
test/RefFEsTests/PolytopesTests.jl
They are the ones being used and tested in the library.

The other two:
src/Geometry/Polytopes.jl
test/PolytopesTests.jl
seem to be there by accident, they are not being used and need to be deleted.

Please, make sure that you are working on the good ones and that you have not relevant changes in the wrong ones.

Thanks!

Find transofrmation between permuted and non-permuted nodes

For a given order, a given Polytope, and a given permutation (as defined in Polytope), provide the transofrmation between permuted and non-permuted interior nodes on top of the Polytope.

Speed up compilation time in complex lazy expression trees on CellMaps

Issue motivated by discussion with @santiagobadia on 2019-06-27

The definition of bilinear and linear forms in the drivers leads to objects that represent complex operation trees. At the moment, in order to allow efficient iteration / indexing of these objects we use a lot of type parameters to achieve type stability. In practical cases, the name of such types can expand the order of 1000 characters, which makes difficult to understand the code and can lead to a considerable compilation overhead (to be better quantified).

Luckily, iteration in the resulting CellMap object is not needed. Only a global call to evaluate with a given CellPoints object is needed. Thus we don't really need type stability. Allowing type instabilty in the concerned types will simplify things a lot.

We have two options:

1. Remove type parameters from ALL the objects representing operations between CellMap objects.

2. Remove type parameters only in the results from operations between FEBasis and FEFunctions since they are the involved in the evaluation of linear and bilinear forms.

To be discussed.

Extend FE assemble for multiple terms

Needed to solve problems with weak forms that include terms that are integrated over different domains.

This is a subtask needed in issue #29

Work being done in branch extend_assembler_to_multiple_fields

Develop a new NodesArray for general polytopes

@santiagobadia

It seems that there is a problem in NodesArray for simplices polytopes:

using Gridap
const t = TET_AXIS
polytope = Polytope((t,t)) # This is a triangle
orders=[1,1]
na = NodesArray(polytope,orders)
length(na.coordinates) == 4 # Wrong? Should be 3, right?

The length of na.coordinates seems to be wrong. Could you confirm that this is a bug?

CellField structs and FE interpolation

I am thinking about the implementation of FE interpolation of functions.

We already have a Field type, which is somehow a field in the physical space.

Thus, we need to create a CellFieldFromAnalyticalFieldWithGeomap that, given an AnalyticalField and a Geomap, is able to evaluate the field in a set of points in the reference cell.

@fverdugo Do we agree on this?

Next, we should use it in evaluate for DofBasis in RefFE.jl. I am not sure the current version is OK, since the restriction of the previous object to a cell should be a Field, otherwise it would not work. I guess it is acceptable.

As far as I can see, there is not any struct now that I could use for these purposes.

Extension to non-oriented facets

For the moment, we have assumed oriented facets

To extend CellFieldFromComposeExtended to CellBases

The functionality provided by CellFieldFromComposeExtended has to be extended to support CellBasis objects.

For instance, this would be required to remove the commented part in line:

Gridap.jl/test/FESpacesTests/NonLinearFEOperatorsTests.jl

Line 53 in c077720

jac(u,v,du) = a(u,v,du) # + varinner(v,ν(du)*grad(u))

Essentially, one needs to extend CellFieldFromComposeExtended to CellMapFromComposeExtended, which in turn will require to extend FieldFromComposeExtended to MapFromComposeExtended.

EDIT

As a part of this re factoring, it would be useful that we can create CellFields giving more than one auxiliary cell field. (useful for Multi-field problems). E.g.,

νfun(x,u1,u2) = (u2*u1+1.0)*x[1]
ν(u1,u2) = CellField(trian,νfun,u1,u2)

Problems with coords

I am having trouble with the coordinates extracted from a trinagulation.

@fverdugo can you take a look at FESpaceTests l.42

D=2
nparts1d = 2
nparts = nparts1d*ones(Int64,D)
nparts_t = tuple(nparts...)
grid = CartesianGrid(partition=nparts_t,domain=(0,1,0,1),order=1) # domain, and order are optional
trian = triangulation(grid)
coords = cellcoordinates(trian)
size(coords)
IndexStyle(coords)
coords.lid_to_gid
coords.gid_to_val
coords.cv
length(coords)
isa(coords, AbstractArray)

provide reasonable values. On the other hand, coords is wrong, it only iterates till cell 2 for a 2x2 trian. Weird behaviour, it is a sub-type of AbstractArray, I can get all indices OK, but when I iterate, it does not work... The iterate for a AbstractArray is in Julia...

FullGridGraph not strictly needed in some places

E.g., for constructing a BoundaryGrid, now we use a FullGridGraph but only a portion is needed. The same for constructing a conforming fe space.

Numerical integration on the grid skeleton

Work done in branch integration_on_skeleton

UnstructuredGrids is an unregistered package

As discussed in JuliaLang/Pkg.jl#681, a package cannot depend on a package that depends on an unregistered package.

Now Numa depends on the unregistered package UnstructuredGrids.

Thus, we need to register UnstructuredGrids (either in the General register or in our own register, first option preferred)

Once this is done, delete the hack included in .travis.yml in order to allow the usage of the unregistered package UnstructuredGrids

Geometry

New API for geometry-related objects available

Tests with julia nighlty builds are failing

To be explored why...

Implement CurlGradMonomialBasis

Needed for RT-FEs. The signature should be like this:

b = CurlGradMonomialBasis(T,3) # 3,2,2 x 2,3,2 x 2,2,3

Work on branch cell_values_and_cell_maps_refactoring

Misc. comments on the developments done at branch cell_values_and_cell_maps_refactoring

Fixing small issues in CellValues.Append

@santiagobadia I have added some comments in the source code of CellValues.Append with possible improvements.

Iteration over CellMaps still leads to allocation of objects at each step

This should be fixed. The solution is to make mutable some of the types in Maps.

Non-continuous FEFunctions for high order FESpaces

@santiagobadia I have realized that high order fe spaces are not working properly. E.g.:

model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=(2,2))
order = 2
tags = [1,2,3,4,6,5]
fespace = ConformingFESpace(Float64,model,order,tags)
fun(x) = x[1] + x[2]
uh = interpolate(fespace,fun)
writevtk(trian,"trian", nref=4, cellfields=["uh"=>uh])

leads to a non-continuous FEFunction (here the model has 2x2 cells but I am visualizing with a finer mesh)

Everything OK for the same example with order = 1

perhaps the problem is in the generation of the cell_eqclass.

Comments after 25th April meeting

Merge CellValue and CellArray, since CellArray = CellValue{T<:AbstractArray}, and only implement for CellArray those methods in which efficiency is important, and replace CellData by CellValue

Replace range_sizeby a return_size that comes from the Point{D} array in the evaluation method

Create CellFieldFromExtended to represent functions in the reference space

Create CellInterpolator <: CellVector that composes a CellFieldFromExtended and CellValue{D<:DOFBasis} (constant at the beginning)

Performance improvements in generation of FESpaces

Some improvements in the performance of FESpaces are needed, in particular:

• Improve efficiency of low level routines that generate dofs, do the interpolation etc (i.e. introduce function barriers and avoid fine grained allocations). Fixed since PR #54
• Improve efficiency of the interpolation on Dirichlet boundary. (i.e., only loop over boundary cells, some times several tags can have the same function associated to them)

Add way to compute reaction forces

Issue motivated by discussion with @santiagobadia on 2019-06-27

In some applications (structural mechanics, surface-coupled problems) it will needed to compute reaction forces. To this end we can introduce a new factory DirichletFESpace function that returns a fe space where the Dirichlet dofs are free and the rest are constrained to zero (i.e., to opposite than in TestFESpace). When using this fe space as the test space in a FEOperator the resulting residual vector will be related with the reaction forces.

Further developments in Polytope and RefFE for high order FEM

1. For the moment I do not have a NodesArray. It is quite simple to create one, but I am not sure it is needed. In fact, it was not in RefFE. Adding it should be straighforward.

2. For the moment, I keep the old one because for anisotropic order the new version is not working yet, we consider same order in all dimensions. On the other hand, the generation of monomials is only working for n-cubes and n-tets. We should probably think about changes in the polynomial machinery to include other cases.

As a result, there are two missing parts:

• Create heterogeneous order nodes on n-cubes using the new machinery and eliminate the old NodesArray

• Create a new NodesArray if needed. The only thing I have not implemented is the set of nodes in the closure of a n-face. Is that needed? It seems that it is not used in the code but I guess it will be needed in facet integration for non-conforming DG methods

• Think about more general monomial generation for high order FEM

I still need to define:

• Provide a general concept of orientation that will be general for whatever dimension or polytope Since this is a topic by its own, I (@fverdugo ) have created a specific issue for it (see issue #51 )

ConstantCellArray out of bounds

In the current implementation, ConstantCellArrays provides the same array even when exceding its length. Is this the behavior we want?

Style guidelines

Discussion to set a consistent coding style

Polytope to Grid conversion working only for particular cases

The conversion currently assumes that all faces on the newly created grid have the same extrusion code. That is, it only works for hexs and tets at the moment, not for wedges, pyramids, etc. If one tries to perform the conversion for a case that violates this assumption, an error is issued by a @notimplemented

How to compute the "change of surface" for facet integration in arbitrary dimensions?

Given the D-1 tangent vectors (i.e. the Jacobian of a geomap from D-1 to D dimensions) how to compute the scalar value representing the corresponding "change of surface" ?

In D=3, I am taking the 2-norm of the cross product of the D-1 tangent vectors

In D=2, I am taking the 2 norm of the tangent vector

I even have a formula for D=4....

But, how to do it in general? is there any way?

Wrong numerical integration

It seems that there is a problem with the numerical integration

D=2
nparts1d = 2
nparts = nparts1d*ones(Int64,D)
nparts_t = tuple(nparts...)
grid = CartesianGrid(partition=nparts_t,domain=(0,1,0,1),order=1) # domain, and order are optional
trian = triangulation(grid) # Generates the Triangulation associated with this grid
quad = quadrature(trian,order=2)
fun(x::Point{2}) = 1.0
ksca = integrate(fun, trian, quad)
# ksca = 0.25

where it should provide ksca = 1.0

Module MultiNonLinearFEOperatorsTests crashes for order > 1

Module MultiNonLinearFEOperatorsTests crashes for order > 1.

Extension to multiple fields

Place to store misc comments that will arise in the extension to multiple fields.

Field values should not be AbstractArrays

Some hacks where introduced in CellValues to circumvent the fact the VectorValue and TensorValue are in fact arrays. This would be simplified if those types are not arrays.

• provide a new implementation of the types in FieldValues so that they do not extend Abstract Array. This will be done by composing an static array instead of defining constants from static arrays. This work can be isolated in a package that defines the types and the associated operations.

Misc. improvements for FESpaces, Assemblers, and FEOperators

Design improvements:

• Improve the expressivity of the abstract interface of FEOperator so that an efficient non-linear FE solver can be implemented only in terms of this abstract interface.
• Allow Neuman boundary conditions, and terms that are integrated over surfaces
• The abstract interface of Assembler assumes 1 field, 1 term (this is ok for the moment)
• Extension to multi-fields

Efficiency improvements:

• Improve efficiency of low level routines that generate dofs, do the interpolation etc (i.e. introduce function barriers and avoid fine grained allocations).
• Improve efficiency of the interpolation on Dirichlet boundary. (i.e., only loop over boundary cells, some times several tags can have the same function associated to them)

Apply is implemented now upto 6 CellValue arguments

The apply function can take up to 6 arguments for the moment. This can be easily generalized by using generated functions.

Provide a general concept of orientation that will be general for whatever dimension or polytope

Place to discuss on this topic