How should the finite element method (FEM) framework in the language be extended to be more useful?

With the release of version 12.0 all fundamental FEM solvers (linear, nonlinear, stationary, transient, harmonic, parametric, eigensolver) are implemented. As many of you know I am a developer of the FEM in Mathematica. As such I do not have questions about the language or framework to ask here; my primary purpose on this site is to help you make the most of the FEM framework. However, I would like to give people on this site that are actively using the FEM framework a voice in what you think could be useful extensions/improvements for the framework.

What are suggestions for improvement or missing functionality that you think would make your work with the FEM easier?

When you write an answer, please try to be as specific as you can, possibly show code that illustrates the problem. Limit your answer to one item, multiple entries are of course OK. Try to be reasonable. Suggestions do not need to be complicated; it can be as simple as tutorial XYZ should have a sentence about ZZZ. With up votes given to various suggestions I will hopefully get an idea what is useful to most people and can prioritize accordingly. Also, please understand that I can not give a commitment that everything requested will/can be implemented and it may take some time before things requested actually see the light of day in the product.

Update 12.1:

I'd like to point point out additions to the FEM framework that fix or alleviate the requests put forward here.

• FEM Programming tutorial extensions. Here I added more examples of how to make use of the low level functions. For example there is a new section on Transient PDEs with Nonlinear Transient Coefficients with this you can model phase change for example. Another new section Transient PDEs with Integral Coefficients shows how to solve transient integral PDEs. These additions are to alleviate this request.
• There is a new tutorial NDSolve Options for Finite Elements on all possible options for the stationary finite element solver. The time dependent options will follow in a future version. This is to alleviate this and in particular this request. Where the second one is not fully fulfilled because it lacks specific application examples. This will remain the case until I get customer examples that I can share.
• OpenCascaseLink. The link provides an initial interface to OpenCascade's Computer Aided Design (CAD) engine. Among many features there is also a new boundary mesh generator called "OpenCascade" that works well for 3D symbolic boolean regions. It's not the default yet depending on how it behaves in the wild it may become the default in a future version. What also may be of interest is the capability to read and write some STEP files (AP203/AP214). This addition is to alleviate this request and partially this one.
• PDE model tutorial extensions. The PDEModels Overview shows the current PDE models available. We now have tutorials for Acoustics and HeatTransfer. Additionally, there are application examples model from Acoustics, Fluid Dynamics, Heat Transfer and Multiphysics. These are long modeling examples. Also you find links to short documentation examples on this overview page. This is certainly something we will see more of in the future. These additions are to start to address this request.
• Iterative solvers. This was not explicitly requested here, but I could imagine this is of interest to some people here too. Both the FEM Options tutorial and the FEM Usage Tips tutorial have sections on how to make use the iterative solvers.

Update 12.1.1:

• A new Mass Transport PDE model tutorial has been added. Accompanying the tutorial two application examples have been added: Microscale Simulation of Catalyst Deactivation and Catalytic Converter
• The OpenCascadeLink got a few updates and is now available from the Wolfram GitHub page

Update 12.2:

• At the Virtual Wolfram Technology conference 2020 I held a FEM Meetup where I discussed questions and suggestions collected on this page. You can see the Video of the FEM Meetup 2020.
• We started to implement a PDE modeling framework (Overview Video). The purpose of this is to make PDE setup easier. The framework consists of basic PDE terms that can be combined to more extensive 'PDE components' to make PDE models from various fields of physics. Currently implemented are Acoustics, HeatTransfer and MassTransport. What is new is that each of those are accompanied by area specific boundary conditions that evaluate to the proper NeumannValue or DirichletCondition.

This is how something then looks like:

vars = {p[t, x], t, {x}};
pars = <|"Material" -> Entity["Element", "Tungsten"]|>;
AcousticPDEComponent[vars, pars] == 
 AcousticAbsorbingValue[x == 1, vars, pars]

• All in all 32 (!) new reference pages with details about the PDE terms and components.
• There is a new mass transport model about Gas Absorption and all other PDE modeling related tutorial and monographs have been updated to make use of the PDE modeling framework and some got new sections like the Interphase Mass Transfer
• The OpenCascadeLink got a few updates, bug fixes and documentation improvements. For example, OCL can now deal with TransformedRegion and got a Torus graphics primitive.

Update 12.3:

• The "OpenCascade" boundary mesh generator is now the default for boolean regions in 3D.
• The OpenCascadeLink has been improved and extended. Example CAD models have been added. A CAD model of simple book shelf bracket and a CAD model a complicated Helical bevel gear are available.
• Working with multimaterials in PDEComponents has been made easier.

This now works:

HeatTransferPDEComponent[{T[t, x, y], t, {x, y}}, <|
  "Material" -> {{y <= 1, Entity["Element", "Tungsten"]}, {y > 1, 
     Entity["Element", "Titanium"]}}|>]

• As requested in comments under this answer, convex hull and Delaunay meshes can now also be generated ToBoundaryMesh["Coordinates" -> pts] and ToElementMesh["Coordinates" -> pts]

Update 13.0:

• The PDE modeling framework has been extended by the field of linear elastic SolidMechanics; There is a SolidMechanics monograph and a verification notebook and a set of companion functions to set up the SolidMechanicsPDEComponent, compute the SolidMechanicsStrain and SolidMechanicsStress and boundary conditions. Here is a recording of a Solid Mechanics presentation.
• There are two new electromagnetics models: Single aperture scalar diffraction and MultipleApertureVectorDiffraction
• There is the new element type PrismElement
• ElementMeshRegionProduct does the same as RegionProduct just for ElementMesh.
• ToGradedMesh generates graded (anisotropic) meshes 1D and is to be used with ElementMeshRegionProduct
• There is a speed improvement when using interpolating functions as PDE coefficients. If these interpolating functions are based on the same ElementMesh as the PDE is to be solved over, then a speed up is attained.

For example running this code:

mesh = ToElementMesh[Rectangle[], MaxCellMeasure -> 0.0005];
fun = ElementMeshInterpolation[mesh, Sqrt[Total[mesh["Coordinates"]^2, {2}]]];
RepeatedTiming[
 sol = NDSolveValue[{-Laplacian[u[x, y], {x, y}] + fun[x, y]*u[x, y] ==1, DirichletCondition[u[x, y] == 0, x == 0]}, u, 
    Element[{x, y}, mesh]];]

Takes about 0.16 seconds with version 12.3.1 and 0.08 seconds in version 13.0. (and 0.04 seconds with the coefficient Sqrt[x^2+y^2]). So using interpolating functions (with the same mesh) is much more performant then before.

Update 13.1:

The PDE modeling framework has been extended further:

• The linear elastic SolidMechanics monograph has been extended. The paragraphs added or modified are marked as such.

• The SolidMechanicsPDEComponent has been extended to allow for hyperelastic material models. Specifically, a St. Venant-Kirchhoff and a neo-Hookeam model have been implemented. There is a new monograph on Hyperelasticity that explains how these models work and shows how to add your own models.

• The basic xyzPDETerms and the MassTransportPDEComponent and HeatTransferPDEComponent now accept a parameter <|"RegionSymmetry"->"Axisymmetric"|> that generates the PDE for a truncated cylindrical coordinate system to make use of axisymmetric geometries. Other xyzPDEComponents will follow suite in the next releases.

• Besides the ref pages there are three new models that show case the use of this: Spherical Capacitor, Thermal Analysis of a Disc Brake and Radial Effects in a Tubular Reactor

• There is a new multi material model: Heat Conduction in a Multilayer Sphere

ParametricFunction can now take an option during function evaluation when the FEM is used. This is useful to give the solver needs an updated initial seed to find the solution of a highly nonlinear problem. Here is a pseudo code:

The parametric function is constructed as usual

pfun = ParametricNDSolve[FEMModel, {u[x, ...], ..}, 
   Element[{x, ..}, mesh], p];

Previously, one could not give an option when the pfun is evaluated and you can now.

pfun[pNew, "InitialSeeding" -> {u[x, ..] == oldUSolution, ..}]

This is useful for solving extremely nonlinear PDEs where the solution can not be found in one go. You can then iterate to the solution like show in the following:

{uSolution, ..} = {0 &, ..};
Do[
 pNew = step*pMax/nsteps;
 {uSolution, ..} = 
  pfun[pNew, "InitialSeeding" -> {u[x, ..] == uSolution, ..}];
 , {step, 1, nsteps, 1}]

All this is explained and showcased in a section on the Hyperleastic Material Models. One caveat: You can not use symbolic expressions for the restart of the initial seed, as that would mess up analysis that ParametricNDSolve does when called.