Future enhancements for the finite element method

Question

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 14.0:

The PDE modeling framework has been extended further:

• The SolidMechanicsPDEComponent now supports transversely isotropic material models. This is useful for modeling fiber reinforced material or material properties that are not axis aligned, such as wood.

• The SolidMechanicsPDEComponent now supports more hyperelastic material models. Compressible and nearly incompressible models are available. Multiplicative decompositions are now supported, to allow coupling of several fields of physics. We have given special attention to make model calibration dead easy.

• New application models have been added to showcase the functionality, for example a Hyperelastic Model Comparison or the Hygroscopic Swelling that shows hyperelasticity with moisture and thermal effects.

• We started to explore electromagnetism. As a first step we support electrostatics and added an electrostatics monograph and functions such as ElectrostaticPDEComponent. We will see more of this field in the future.

• We added support for fluid dynamics. A monograph on Laminar Flow is available. This shows how to set up the Navier-Stokes equations, couple them with other fields of physics. Non-Newtonian flow is supported with build in models, that can be extended. The main function is FluidFowPDEComponent

• The physicists among you will hopefully enjoy the new SchrodingerPDEComponent

• Many other updates went into the PDEModels.

FEM related updates:

• There is the new DiscontinuousInterpolatingFunction, that serves two purposes: To handle discontinuous functions and when the derivative of a function is discontinuous. This is best explained with an example.

We have a two material region:

Needs["NDSolve`FEM`"]
mesh = ToElementMesh[
   ToBoundaryMesh[
    "Coordinates" -> {{0, 0}, {1, 0}, {1, 1/2}, {1, 1}, {0, 1}, {0, 
       1/2}}, "BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 
         4}, {4, 5}, {5, 6}, {6, 1}, {3, 6}}]}], MaxCellMeasure -> 1, 
   "RegionMarker" -> {{{1/2, 1/3}, 1}, {{1/2, 2/3}, 2}}];
mesh["Wireframe"[
  "MeshElementStyle" -> {FaceForm[LightGreen], FaceForm[LightRed]}]]

Now, we evaluate a discontinuous function over the mesh:

function = If[ElementMarker == 1, 0, 1];
(* old behaviour *)
ifunOLD = 
  EvaluateOnElementMesh[{x, y}, function, mesh, 
   "DiscontinuousInterpolation" -> False];
ifunNEW = EvaluateOnElementMesh[{x, y}, function, mesh];

Plot3D[{ifunOLD[x, y], ifunNEW[x, y]}, {x, 0, 1}, {y, 0, 1}, 
 PlotRange -> All]

You see the new function sharply represents the discontinuity.

You can precisely control which parts of the boundary belong to which sub-region. Here is a more elaborate example.

Update 13.3:

The PDE modeling framework has been extended further:

• The SolidMechanicsPDEComponent now supports axisymmetric models.

• The SolidMechanicsPDEComponent now supports the nearly incompressible Yeoh material model. The Hyperelasticity monograph has been extended to explain the Yeoh model and build the infrastructure for more hyperelastic materials to be added.

The following PDE model application examples have been added to the PDEModels Overview:

• A Passive Dew Condenser (Fluid Dynamics + Heat Transfer)

• A Contactless Anemometer (Heat Transfer)

• A Quantum Ring makes the first example in the new 'Physics' category. This section also showcases two more elaborate documentation examples: One on Avoided Crossings and one on Eigenfunction clustering

• In the solid mechanics domain, there is the Biaxial Tensile Test of Hyperelastic Tissue and the Layered Vascular Vessel with Yeoh Constitutive Model both featuring hyperelastic material models. The Disc Brake has been updated and now also features an axisymmetric solid mechanics model. Note, that the parameter ModelFunction has been renamed to SolidMechanicsModelFunction. The same rename holds true for the parameters ModelForm which was renamed to SolidMechanicsModelForm Old code will continue to work, but you should be aware of this change. Further more functionality to compute ElasticEnergy has been added.

Updates for the FEM code include the following:

It has be come simpler to specify mesh elements with markers. Previously, each element needed it's own marker. There is now a simplification if each element uses the same marker the the marker needs to be specified only once. This now works:

Needs["NDSolve`FEM`"]
bmesh = ToBoundaryMesh[
"Coordinates" -> {{0, 0}, {1, 0}, {1, 1}, {0, 1}},
"BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 4}}, 1],
LineElement[{{4, 1}}, 2]}]; 

BoundaryUnitNormal is now documented.

Update 13.2:

One of the most important updates in 13.2 is actually a system wide improvement. Consider this code:

NDSolveValue[{Laplacian[u[x, y], {x, y}] == 1}, 
 u[x, y], {x, y} \[Element] Disk[]]

Note the circle-i

at the end of two of the three messages. This information link is only present iff the error message given has a message reference page. If you click on the circle-i, it will take you to the respective message reference page. There you will find hints for isolating and solving the issue at hand. Again, this circle-i is only issued if there is such a message ref page. This should help to better understand why NDSolve issued a message and how to fix it. Kudos go to the FrontEnd team for implementing this for me.

The PDE modeling framework has been extended further:

• Almost all PDEComponents and PDETerms now support "RegionSymmetry"->"Axisymmetric". The only exception is the SolidMechanicsPDEComponent which will provide "RegionSymmetry" support in the next version (13.3)

• The SolidMechanics monograph now has a section explaining the minimal number of constraints needed for solid mechanics modeling.

• The SolidMechanicsPDEComponent now supports Kirchhoff stresses.

• A new acoustics model, the Helmholtz Resonator has been added.

Here is a Video on the usage of axisymmetry and hyperelasticity

A new workflow for creating finite element meshes from elevation data (contours) has been added.

OpenCascadeLink has been updated to use OpenCascade 7.6.0

• now supports for closed B-spline conversions as long as the knots and weights are set to Automatic.

• now allows for surface mesh simplifications (elimination of unnecessary edges and faces) after boolean options.

• now allows shape simplifications

• now provides shape healing capabilities

Overview presentation: If you'd like to get a very rough overview of the FEM functionality in version 13.2 have a look at the PDE Modeling overview talk I gave.

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-Hookean 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 SolidMechanicsPDEComponent has been extended to handle Rayleigh Damping.

• 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.

When you use ClearSystemCache[] then the FEM case will also be cleared in some cases this can reduce the memory footprint.

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 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 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.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.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.

Answer 1

A useful feature that I regularly use in COMSOL and would like to be able to use in Mma is the "AdaptiveMeshRefinement" (as it is called in COMSOL).

This means that COMSOL makes a mesh. With this mesh, it solves the problem. Then it evaluates a function that characterizes the steepness of the solution. Typically, it is the gradient of the solution squared, but it can also be a user-defined one. Then COMSOL transforms the previous mesh such that it becomes denser in the place, where this function has a higher value, and which may grow coarser in regions where this function is smaller. Then it solves the problem with a new mesh. It repeats such a refinement several times.

The number of mesh refinements during one run can be adjusted. One controls the refinement by specific parameters. One of them, for example, can define how many times the mesh size decreases (or increases). Another may determine the way of the division of the mesh cell.

Let us note that in COMSOL one does not really allow for varying all such parameters, and some tuning settings do not function, but some of their combination do work, and I use them. Yet, I did not see anything like this in MMA. However, I feel it be advantageous.

Answer 2

I think user21 needs to be congratulated for developing the finite element method and for asking this question. My thoughts are as follows:

1. The purpose of finite elements is to solve differential equations on complex geometries.

2. The goal of the Wolfram Language is simple, if ambitious: have everything be right there, in the language, and be as automatic as possible. Quote from blog by Stephen Wolfram May 21, 2019 here.

3. There is a large industrial usage of finite elements for engineering. Stress and dynamics being possibly the big users.

There are three stages in a finite element calculation. Preprocessing, Solving and Postprocessing.

The Wolfram Language ought to be good at preprocessing and sorting out the differential equations. However, this is difficult and does not correspond to Wolfram's point in 2 above. To solve stress problems you have to coerce textbook equations into this form

where the $ c_{i j}$ are 3 by 3 matrices. I have tried but failed to do this although user21 has provided a working version here. First request: can we make formulating equations and coercing them into the correct form straightforward. Examples would be helpful. I will perhaps post elsewhere where I have got stuck in this process. Also, there are variants of the stress equations and nonlinear stress problems that need to be formulated.

The other issue with preprocessing is making a good mesh. This means building a good solid model and meshing. At the moment this means discretizing early using BoundaryDiscretizeRegion which does not lead to a good mesh. Further we only have second order meshes and calculating stress requires the derivatives of the displacements. Thus the stresses have only first order interpolation. Either we need higher order mesh interpolation or the ability to use very fine meshes. This is along the lines of the h -p question Second request: more solid modelling and meshing capability.

The solving stage is up to the Wolfram language numerics. Will they be capable of solving industrial engineering solutions mentioned in point 3 above? This is very much a policy question for Wolfram. Big engineering problems or only toy problems by comparison.

Finally a comment on post processing. This is where the Wolfram Language is good. You don't have to learn a new language. This is a strong point for developing finite elements in the Wolfram Language.

Finally a comment on solving fluid problems. As I understand it these are the really big problems for which no mesh is adequate. Solving fluid flow at large Reynolds numbers is not usually done in finite elements but in a finite difference formulation. A vast range of turbulence models are used the simplest being $k-\epsilon$ used with wall functions. Is this outside the scope of what is being considered?

Update 12.1 (user21):

Please see:

• the update PDEModels Overview page
• the Using OpenCascadeLink tutorial

Update 12.2 (user21):

A framework for PDE terms and components has been added. With more fields to come in the future.

Update 13.0 (user21):

A linear elastic SolidMechanics PDE modeling frame work has been added.

Update 13.1 (user21):

• Hyperelasticity has been added to the solid mechanics.

• Rayleigh Damping has been added.

• 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.

Answer 3

In my opinion one thing that is still missing for really useful FEM framework is better quality of meshing (of Boolean representations of geometries) in 3D (ToElementMesh). I know this is not an easy task, but I would still like to include it on wishlist.

For example:

Get["NDSolve`FEM`"]

box = Cuboid[{0, 0, 0}, {1, 1, 1}];
holes = Thread@Ball[{{1., 0.5, 0.5}, {1., 1., 0.5}, {1., 1., 1.}}, 0.2];
reg = Fold[RegionDifference, box, holes];
bounds = RegionBounds[reg];

mesh = ToElementMesh[
  reg,
  bounds,
  MaxCellMeasure -> 0.05
]

Through[{Min, Mean}[Join @@ mesh["Quality"]]]
(* {0.000165709, 0.319868} *)

mesh["Wireframe"[
  "MeshElement" -> "MeshElements",
  "MeshElementStyle" -> FaceForm@LightBlue
]]

The resulting mesh has quite poor quality.

Update 12.1 (user21):

In version 12.1 you can use:

bmesh = ToBoundaryMesh[region, 
   "BoundaryMeshGenerator" -> {"OpenCascade"}];
groups = bmesh["BoundaryElementMarkerUnion"];
temp = Most[Range[0, 1, 1/(Length[groups])]];
colors = ColorData["BrightBands"][#] & /@ temp;
bmesh["Wireframe"["MeshElementStyle" -> FaceForm /@ colors]]

mesh = ToElementMesh[region, 
   "BoundaryMeshGenerator" -> {"OpenCascade"}];
Through[{Min, Mean}[Join @@ mesh["Quality"]]]

{0.0458246, 0.695077}

mesh["Wireframe"["MeshElement" -> "MeshElements", 
  "MeshElementStyle" -> FaceForm@LightBlue]]

Update 12.2 (user21):

More graphics primitives and operations have been added. Please see the updated Using OpenCascadeLink tutorial.

Update 12.3 (user21):

The "OpenCascade" boundary mesh generator is now the default for boolean regions in 3D. So this is now sufficient:

ToElementMesh[region];

Besides that the OpenCascadeLink has been improved and extended. Among other things it now has two application examples showing the creation of CAD models: a simple book shelf bracket and a complicated Helical bevel gear:

Answer 4

It is obligatory that I make a wish for finite elements on immersed curves and surfaces. This has a plethora of applications in geometry processing, but also in physics, chemistry and microbiology. Here is a short, incomplete list of posts that could have been solved easier with surface FEM:

1. How to estimate geodesics on discrete surfaces?

2. Smoothing 3D contours as post processing

3. Can Mathematica solve Plateau's problem (finding a minimal surface with specified boundary)?

4. How to apply different equations to different parts of a geometry in PDE?

Surface FEM can be added with reasonable effort because first order elements can be implemented straightforwardly with essentially the same techniques as for full-dimensional domains. Also the data types for the meshes are already out there.

Answer 5

I think that it might be beneficial to write down the tutorial describing the ways to choose and to fine-tune the solvers used. This proposal is close to that of @Rom38, but slightly differs from his one.

The point is that different equations require different fine-tuning methods. Technically, I can imagine that one can demonstrate a few methods on one equation, other few ones on another and so on. Like this, one will be able to show all the main techniques.

It will be ideal if one gives these techniques with some comments explaining why he has applied this or that method. However, I guess that sometimes one knows why the way is suitable, but in some instances, one needs simply to try. The fact that there is no clear indication of what to apply in this case is also advantageous to write directly as the explanation.

Anyway, it would be of great advantage for the users to have various examples of such fine-tuning approaches before the eyes.

One problem here is that the developer (user21) has in mind particular examples of equations, and actually, we see these examples in the existing tutorials. We, however, deal with other examples of equations challenging to solve. And it is for these equations that we need some specific fine-tuning.

I propose that we can post examples of nonlinear equations that we can imagine to be of general interest, or mail them to the user21 as examples. This will enable user21 to collect a pool of equations to take examples.

Writing such a tutorial is in no case simple. I guess that it is a task for a considerable time. After all, one has to (1) collect many examples and (2) solve them all. However, I believe that such a tutorial will has a potential to make FEM in MMA to be a real working instrument.

Update 12.1 (user21):

Please see:

• the new Finite Element Options tutorial

While this tutorial does not address all issues mentioned here it forms a basis by collecting all options for (stationary) FEM in one place and explaining what they are for and where to find more information. This is at least an overview of what one can try to do to solve stubborn PDEs.

Answer 6

Support for PDE Whose Spatial Derivative Order Exceeds 2

I've been stopped in v9 for a long time and don't consider myself as somebody actively using the FEM framework, but since nobody has mentioned this for so long, I'd like to add. According to the FEM-related question coming out here, this seems to be the most needed missing functionality. Just search femcmsd in this site, you'll see… only 9 related posts? Well, perhaps the keyword is not always included…

Answer 7

I guess, that one of the best improvements will be the detailed guide "how it works". I mean, for example the step-by-step solution of let's say transient 2D (or even 3D) heat transfer equation with heat sources (or anything else) with application of the main perfomance tweaks (mesh configuration, submethods with comments about effects, etc).

The primitive examples that present now are not clear about details of configuration..

Update 12.1 (user21):

Please see:

• the (updated) Finite Element Programming tutorial
• the new Finite Element Options tutorial

Answer 8

Anisotropic Meshing

The following is a feature request for anisotropic meshing. A proper mesh is as or more important than just having the proper equations to obtain accurate simulation results.

The problem arises when one needs to capture either very close or very far away from some feature in the main model. Naïvely meshing the model with the uniform mesh size will cause the model size to blow up. The problem is exacerbated by going to higher model dimensions in both space and time.

Fortunately, the FEM is often quite robust to element aspect ratios for many types of physics. This allows one to use very flat or very stretched-out elements while simultaneously reducing model size while maintaining accuracy. Boundary layer meshing and infinite domain elements are types of anisotropic meshing commonly found in FEM packages.

Without boundary layer meshing, one will often over-predict shear stresses, heat, and mass transfer rates at the wall in fluid flow problems. Without infinite domain elements, the slopes of dependent variables will be in error due to truncation.

I have used anisotropic meshing to solve various problems on Mathematica Stackexchange, as shown in the following list.

1. 1D Meshes
   • Transient
     • 1D mesh generation for PDE solution
     • Wrong solution from multi-materials FEM NDSolve
     • Neumann boundary condition ignored
     • Dirichlet Condition at Infinity
     • Defining mesh size for NDSolve
2. 2D Meshes
   • Steady-State
     • Mathematica vs. MATLAB: why am I getting different results for PDE with non-constant boundary condition?
     • Improving mesh and NDSolve solution convergence
     • PDE system. convection dominated, method AffineCovariantNewton failed, etc
     • Laplace's equation in spherical coordinates
   • Transient
     • Controlling dynamic time step size in NDSolveValue
     • How to model diffusion through a membrane?
     • Mass Transport FEM Using Quad Mesh
     • NDSolve with equation system with unknown functions defined on different domains
3. 3D-Meshes
   • Create graded mesh
   • Stationary
     • How to Improve FEM Solution with NDSolve?
     • 3D FEM Vector Potential

Anisotropic meshing of complex geometries in 3D

Of the examples I showed in the above bullet list, the geometries were of a simple variety. The most straightforward way to create a boundary layer mesh from a tetrahedral mesh would be to extrude a prism layer. The current FEM Solver does not accept prism elements. So, either the Solver needs to be extended to accept prism elements or a procedure to split prisms into tetrahedra will be required. I am not sure that either option is simple. I have used commercial software to split prisms into tetrahedra, but I did not see too many options available in my cursory search online. Perhaps LaGrit could be used to perform the splitting operation (see Python documentation for grid2grid_prismtotet3).

Update 13.0 (user21):

This has been implemented with ToGradedMesh and ElementMeshRegionProduct. This is probably a particularly fine example of something that has demonstrated it's usefulness here on SE and made it into the product - Thank you Tim for your continuous support of the FEM in the Wolfram language.

Answer 9

I see one more expansion of MMA tools in the FEM for nonlinear PDEs. This is a "Parametric Continuation."

The point is that provided equation has a parameter, say, eps varying from 0 to 1 one starts its solution with eps=0 and MMA solves the equation while gradually increasing the parameter in steps until eps=1. Each next solution takes the result of the previous one as the initial seed.

The main idea is that one can have a nonlinear equation that is much too complex to be solved directly. However, by introducing the parameter eps one can sometimes transform it into a solvable one. Then gradually increasing eps sometimes it is possible to slowly "pull" the solution to eps=1, which is the initial objective.

Update 13.1 (user21)

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.

Another option is to use a pseudo time integration. This has the advantage that the adaptive time stepping of the time integration can be used, but the disadvantage is the need to set up the problem as a time dependent problem and also to think about initial conditions. (This is also shown in the Hyperelasticity monograph mentioned above) My experience is that the time dependent works more generally, but the parametric restart is quicker to set up.

Answer 10

I would greatly appreciate some support for non-local operators. What I have in mind are the fractional powers of the Laplace operator that now appear quite frequently in modeling non-standard diffusions.

Answer 11

This is really about the visualisation of FEM results in 3D, but I post it here since it is related. We have StreamPlot which plots 2D streamlines, which I have used for steady-state results in 2D. And we have VectorPlot (for 2D) and VectorPlot3D (for 3D). Something which would be very useful (and rather natural) would be a StreamPlot3D function. While the differences between VectorPlot and StreamPlot are quite subtle, I have found StreamPlot to be more helpful in my recent applications to 2D (it can be harnessed to make quite sparsely populated plots with seeded streamlines). It would be great to have the analogue for 3D (since VectorPlot3D can make very busy plots which are hard to interpret). Thanks.

Update V13.0 (user21)

There are now VectorDisplacementPlot and VectorDisplacementPlot3D that should be useful for you.

Answer 12

Decouple the Notebook from the Mesh and Solution by Creating Separate Directories

If the vision is to have Mathematica ultimately solve industrial scale problems, then the meshes and solutions will become huge especially when dealing with 3D transients or Lagrangian particle tracing data. I believe real value of the notebook is to document and capture the simulation workflow and not as a storage mechanism for the mesh and solution. Indeed, one small notebook could drive many meshes and solutions by simply pointing to another directory.

Answer 13

I've long wanted to specify problem symmetries and have the mesh and equations modified to support those symmetries. I.e., modified to minimize solution deviation from the given symmetries. (There's probably a "Galerkin with symmetry-preserving basis" hiding in here somewhere...)

Answer 14

I would like a named migration term be added to the MassTransportPDEComponent so all terms of the Smoluchowski or Fokker-Planck equations can be modeled. It appears that currently Mathematica (v 12.2) is missing the ability to add a named migration term except by hand, and such solution strategy is not explicitly described in the documentation. One would expect an example in the documentation showing how to solve particles diffusing under a force or potential using such terms.

Further details and examples are described in this question. Posting this here was suggested in this answer.

Motivation and explicit problems of diffusion under potentials one want to solve with this are described in chapter 4 here.

Update 12.3 (user21)

Examples for both the Fokker–Planck and the Smoluchowski have been added.

Answer 15

It would be nice to update the FEAST solver to the latest (4.0 as of 2020) version to allow for non-Hermitian problems and to benefit from the performance enhancements.

Answer 16

Using Finite Element Method more and more often I would wish an extended version of ElementMeshInterpolation[mesh,par]

To my understanding the "object" ElementMeshInterpolation[mesh] should provide the set of basic elementfunctions. For practical use a splitted definition in the form

`ElementMeshInterpolation[mesh][par]`

would be very useful!

Thanks

Answer 17

I think the setting of boundary conditions should be enhanced for some typical cases, and models should be given in the documents like in COMSOL， although there is a monongraph on acoustics. But I think it is far from sufficient.

I have been working on waveguide mode analysis using FEM in Mathematica for a week, but I haven't succeeded until now.

The optical fiber-like waveguide is featured with different refractive index in core and in clad, and the interface between the core and the clad should have the boundary condition of Dz (the normal component of D) and En (the tangential component of E) are continuous. But I don't know how to express this kind of boundary condition in Mma. I think this is of course different in Neumann, Dirichlet and Robin conditions.

The physical model is described below.

For Helmholtz equation for optical waveguide: $\nabla _{\{x,y,z\}}^2\text{E[x,y,z]+$\epsilon $(}\frac{2 \pi }{\lambda })^2\text{ E[x, y, z]=0}$

Assuming that $\text{E[x,y,z]=E[x,y]*Exp(i*$\beta $*z)}$,

$\nabla _{\{x,y,z\}}^2\text{E[x,y,z]+$\epsilon $(}\frac{2\pi }{\lambda })^2\text{*E[x,y,z]=0}$ becomes to
$\nabla _{\{x,y,\}}^2\text{E[x,y]+($\epsilon $(}\frac{2\pi }{\lambda })^2-\beta ^2\text{)*E[x,y]=0}$

$\epsilon$ is different for core and cladding, i.e. , $\epsilon$core and $\epsilon$clad, respectively.

The boundary conditions at the interface should be : (1) the tangential of the E, i.e. Et, is continuous. (2) the normal component of the D, i.e. Dn, is continous, in which D=$\epsilon$*E. In the cylindrical coordinates of {r, $\theta$, z}, the boundary condition should be Ez and E$\theta$ should be the same, and Dr should be the same at both side of the interface. These contidions are my main concern when using FEM for the analysis of the eigenmode. Although they can be formulated easily in some special cases such as in rectangular or circular waveguide, but I'd like to try a more general form.

Here is my unsuccessful try.( Mma 12.0, Win 10)

To make the mesh points on the boundary, it can be used like this,

<< NDSolve`FEM`

r = 0.8;
outerCirclePoints = 
    With[{r = 2.}, 
      Table[{r Cos[θ], r Sin[θ]}, {θ, 
          Range[0, 2 π, 0.05 π] // 
  Most}]];       (* the outer circle  *)
innerCirclePoints = 
    With[{r = r}, 
      Table[{r Cos[θ], r Sin[θ]}, {θ, 
          Range[0, 2 π, 0.08 π] // 
  Most}]];                             (* the inner circle *)

bmesh = ToBoundaryMesh[
      "Coordinates" -> Join[outerCirclePoints, innerCirclePoints], 
      "BoundaryElements" -> {LineElement[
            Riffle[Range[Length@outerCirclePoints], 
                RotateLeft[Range[Length@outerCirclePoints], 1]] // 
              Partition[#, 2] &], 
          LineElement[
            Riffle[Range[Length@outerCirclePoints + 1, 
                  Length@Join[outerCirclePoints, innerCirclePoints]], 
                RotateLeft[
                  Range[Length@outerCirclePoints + 1, 
                   Length@Join[outerCirclePoints,innerCirclePoints]],1]] //Partition[#,2] &]}];                                                     
    mesh = ToElementMesh[bmesh];
{bmesh["Wireframe"], mesh["Wireframe"]}
 (* generate the boundary and element mesh, to make the mesh points \
on the outer and inner circles   *)

glass = 1.45^2; air = 1.; k0 = (2 π)/1.55;
ϵ[x_, y_] := 
     If[x^2 + y^2 <= r^2, glass, air]

   

helm = \!\(\*SubsuperscriptBox[\(∇\), \({x, y}\), \(2\)]\(u[x,y]\)\) + ϵ[x, y]*k0^2*u[x, y];
boundary = DirichletCondition[u[x, y] == 0., True];

(*region=ImplicitRegion[x^2+y^2≤2.^2,{x,y}];*)

{vals, funs} = NDEigensystem[{helm, boundary}, u[x, y], {x, y} ∈ mesh, 1,Method -> {"Eigensystem" -> {"FEAST","Interval" -> {k0^2, glass* k0^2}}}];
vals

 Table[Plot3D[funs[[i]], {x, y} ∈ mesh, PlotRange -> All, 
    PlotLabel -> vals[[i]]], {i, Length[vals]}]

Although the profile in the figure seems right, but the eigenvalue is not right, since I can check it using analytical solutions here.

---

Edit 1

I notice their is a very closely related post here, where PML is employed. However, there are some bug there, and it couldnot properly run.

Are there some more examples? Thank you in advance.