# Future enhancements for the finite element method

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:

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

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:

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.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. Cudos go to the FrontEnd team for implementing this for me.

The PDE modeling framework has been extended further:

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

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.

• Something I always want is direct access to the internals of a package. Sometimes it simply doesn’t make sense for you or WRI to try to implement everything. Instead I think it’d be really cool if all the work you’ve done here could be easily reused by someone like Henrik to implement their own types of solvers that simply aren’t general enough to qualify for being included in the primary FEM package. Jun 11, 2019 at 7:27
• @b3m2a1, you comment make is sound like the FEM package internals are not documented; however, they are fully documented here and specific function have there ref pages. Also, Henrik has made use of the low level FEM functions, e.g. here. So, generally, I'd say low level package things are documented. If you think something is missing, let me know. Jun 13, 2019 at 11:10
• I use Ansys Mechanical APDL heavily every single day at work but unfortunately I haven't updated my mma since V10. How exactly would you compete with a software like Ansys?
– Öskå
Jun 23, 2019 at 18:04
• @Öskå, your comment does not really relate to my question. For the highlights of the language have a look at the new in 11 and the new in 12. I can not speak for Ansys as eons have past since I last used it. Jun 24, 2019 at 4:08
• @user21, in future can we use isogeometric analysis method in Mathematica (build-in) for high order problems, e.g. cahn-hilliard? Jul 20, 2019 at 23:06

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.

• Hm. Sounds as if this would amount to generating a MeshRefinementFunction from a posteriori error analysis of an already computed solution. Should be doable for linear elliptic problems but might turn out to be a horrible job for nonlinear or transient PDE. May 27, 2019 at 10:03
• @Henrik Schumacher You are right, it should be a difficult task. In COMSOL, however, I use it for a nonlinear time-dependent problem. It contains first time derivative and the Laplacian linear terms plus nonlinear ones. It works. May 27, 2019 at 11:02
• Don't get me wrong: I think adaptive refinement would be great! I just meant to damp expectations a bit, in particular towards convection-dominated parabolic PDE (like Navier-Stokes) and even more towards hyperbolic PDE. For those one might be forced to use different spatial discretizations for different time slices. May 27, 2019 at 11:42
• @AlexeiBoulbitch Adaptive refinement strategy should be very useful for large scale problems! May 28, 2019 at 19:53
• @ABCDEMMM Yes, exactly. Especially in multiscale problems where one simulates an object with a characteristic size, L, bur this object has important features with the characteristic size, say, L/100or even L/1000. It is exactly the case I face with. May 29, 2019 at 8:36

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

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.

• Thanks for your comments and kind words. Concerning your first request, have a look at the PDEModels/tutorial/AcousticsTimeDomain and PDEModels/tutorial/AcousticsFrequencyDomain (the in product version looks better) and let me know if you think this is the right direction. Thinks like AcousticModel could be packaged up for usage. Now, this is not structural mechanics yet, but hopefully we will get there. May 31, 2019 at 4:48
• I like the first idea of predefined common physical models, like AcousticModel from tutorial pages. And @user21, thank you for your effort on extensive documentation, these acoustic tutorials are amazing. May 31, 2019 at 7:33
• @user21 will we have an isogeometric analysis solver in future Mathematica Version? Mar 23, 2021 at 11:28
• @ABCDEMMM, not in the forseeable future. But you can roll your own. Everything is there. Mar 23, 2021 at 12:20

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["NDSolveFEM"]

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,
groups = bmesh["BoundaryElementMarkerUnion"];
temp = Most[Range[0, 1, 1/(Length[groups])]];
colors = ColorData["BrightBands"][#] & /@ temp;
bmesh["Wireframe"["MeshElementStyle" -> FaceForm /@ colors]]


mesh = ToElementMesh[region,
Through[{Min, Mean}[Join @@ mesh["Quality"]]]

{0.0458246, 0.695077}

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


Update 12.2 (user21):

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:

• Could you elaborate a bit on what you mean. Is it more that Boolean representations of geometries do not work well or more that the quality of a mesh (e.g. a simple geometry) is not satisfactory. Perhaps you have an example to add to your suggestion? May 29, 2019 at 10:41

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:

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.

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

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.

# 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…

• Thanks for the suggestion. This would mean that I'd need to write a code that transforms higher order derivatives into systems of equations, like done here May 30, 2019 at 6:29
• Forgot to mention that I also added the example to the help system. You can find it by clicking on the message NDSolve::femcmsd and following the link or by going to FEMDocumentation/ref/message/InitializePDECoefficients/femcmsd May 30, 2019 at 6:48
• as per question mathematica.stackexchange.com/q/91150/1089 :-) Jun 5, 2019 at 17:46
• @user21 does V13 support for PDE Whose Spatial Derivative Order Exceeds 2? Dec 9, 2021 at 21:26
• @ABCDEMMM, no, not for the FEM. The problem is , to be honest, I have not come up with a way to do this. If I do not have a method/idea how this could achieved I am afraid I will not be able to solve this issue. I am sure there is a way, but I have not come across it yet. Dec 10, 2021 at 7:57

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

• Have you seen the Finite Element Programming tutorial that explains how the solvers work? In case you have, could you comment a bit more specifically what you think should be improved for this tutorial that would be helpful for me. May 27, 2019 at 12:25
• Thanks for the link. I did not see it, may be it will be more clear than previously present tutorial. However, the example of fully solved complicated enough task (instead of the primitive examples in help) is desired in such tutorial to see the whole approach in complex. May 28, 2019 at 5:14
• OK, have a look and let me know what could be improved if there is anything. May 28, 2019 at 5:18
• Did you already have a chance to glance over the tutorial? Is this what you are looking for or should something else be added? May 30, 2019 at 5:58
• @user21, yes, I've looked this. It is much better than the previous one but there are just few examples of additional options and even if they are present the effect of most options is not clear. The above answer of Alexei Boulbitch nicely explains this point May 30, 2019 at 6:31

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

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

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.

• @user21 Nice addition to ParametricFunction. Where can I find this capability in the Mathematica documentation? In particular, is this new option consistent with differentiating ParametricFunction with respect to p? Sep 6, 2022 at 2:07
• @bbgodfrey, the only documentation there is is in the section on hyper elastic materials in the solid mechanics documentation; the reason is this feature only tested with FEM. That values for the reinitilization can only be numeric values, precisely because allowing symbolic values will fiddle with the differentiation of the object. Right now (V13.1) I doubt that I will develop this any further as the time stepping approach gives you an automatic adaptive step length. The step size in the initial seeding approach is fixed, or manual at best. Sep 6, 2022 at 8:11
• So, while this approach is more convenient to set up it's not as general as the time stepping. Sep 6, 2022 at 8:11

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.

• You mean things like this? May 30, 2019 at 5:05
• @user21 I also want this feature. If can use the operator such as $I[f]=\int f(x+z) dP(dz)$ or $max_{e} [f(x+e)]$(jump) in NDSolve,most of partial integro-differential equations can be solved easily... May 31, 2019 at 9:08
• @user21 Thanks for the link, yes that's exactly what I had in mind. The link you provide seems to require an advanced understanding of FEM in Mathematica. I wish we could run it as easily as solving the regular Dirchlet problem. Next, I also had in mind the possibility of building non-linear terms with it. May 31, 2019 at 16:57
• The problem is that for nonlocal operators, quite different techniques have to be applied. Since the stiffness matrix will be a dense matrix, one has to invest a lot of effort in implementing (i) matrix-free matrix-vector multiplication in an efficient way and (ii) a good preconditioning technique. For (i), one has to apply methods such as Barnes-Hut tree codes, fast multipole method, or hierarchical matrices. That is definitely a non-trivial task and essentially non of the classical tools for FEM can be reused. For (ii), see here. Jun 2, 2019 at 11:44
• Btw.: I know this so well that this is hard work because I am actually working towards that direction for a project of mine. So maybe, I will post something somewhere here someday. But I doubt that nonlocal operators have a ''market'' sufficiently large so that developping a general-purpose solver for Mathematica is justified. Jun 2, 2019 at 11:49

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.

• Paul, I very much sympathize with your request and I forwarded this to the relevant developer. However, this is not something I can decide or implement as much as I'd like to. If you feel strongly about this then you could perhaps contact support AT wolfram.com and put forward your request. This way there will be a record of demand for such functionality. Feb 20, 2020 at 5:40
• @user21. Thanks for your kind reply. I will go ahead and make such a request, and we'll see what happens. Feb 20, 2020 at 14:22
• Have made a request as suggested and await a reply. In the meantime, I can report that I was able to code my own StreamPlot3D function, (it turned out to be easier than expected). A quick sketch: start at a seed point and step along field lines in small steps of a distance variable "t". At each step, store the new position as a 3D space point. Also store the stepping variable, t. Then create 3 separate 2D lists, one each for (x,t), (y,t) and (z,t) which are each made into Interpolating functions. Then a simple ParametricPlot3D works a treat. Can superpose on graphic of region as required. Feb 23, 2020 at 12:32

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

• I have trouble understanding your suggestion. Could you elaborate a bit on what you mean? Are you suggesting that given a PDE and a mesh the solver should automatically find symmetries and exploit those? May 28, 2019 at 4:23
• @user21 : Well, that would be great, but, except for contact symmetries and very specific discrete symmetries, this seems out of reach. Rather, given a user-supplied list of symmetries, minimize deviation from them. Simplest example: isotropic 2-d wave equation on square grids. Rotational symmetry may or may not be important in an application but (except for stupidly short time steps) is not a feature of the discrete solutions. If a user specifies it, perturb the difference equations to minimize violation of that symmetry. May 28, 2019 at 6:25
• @EricTowers You are basically asking for trouble. My user icon should remind you of that... ;) Mar 23, 2020 at 8:23

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

• As I understood, all the data are stored in kernel. Notebook can have it just if you want to see it as numbers of as figure. Which coupling do you mean? You can run separate kernels for each solution with corresponding mesh.. May 30, 2019 at 4:58
• Thanks for you suggestion Tm, are perhaps missing an example on how to export and import meshes and interpolating functions? I could add this to the section dealing with large scale FEM in the Finite Element usage tips tutorial. Would that help? May 30, 2019 at 5:55
• @user21 Yes, I think examples of import/export of meshing and interpolations functions would be useful and also how to set up a parallel example. It would help me understand the process better. My broader long term question is what If I want to run a long transient simulation on a large mesh based off of 3D CAD where I anticipate the final solution size to be 10x my RAM, then how would MMA respond to that situation? Most of the solvers I use just pump new time step files into a solution directory freeing up memory for the next time step. MMA, I suspect, would be different. May 31, 2019 at 5:23

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.

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.

• While I very much sympathize with your request, this is nothing in my control. This needs to be done by the developer of Eigensystem which I am not. But there is more complication, the FEAST eigensolver used is the one from MKL (intel) and it depends on them when they make these new features available. But I'll pass the request on. Jun 22, 2020 at 18:47
• @user21 can OPEN-MPI be supported in future version? Jun 3, 2021 at 1:06
• @ABCDEMMM, this is not my choice, the developer of Eigensystem will have to make that choice. Also, as far as I know the OpenMPI version is not available in the MKL-FEAST version we use. What you could do is write an email to the support and suggest it's addition. Perhaps you can include an example to support your case. Jun 3, 2021 at 3:50

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

• You mean like this: ElementMeshInterpolation[mesh,#]& - could you add an example that explains in what cases this would be useful. Jan 11, 2022 at 9:50
• Thanks for the quick reply. Because ElementMeshInterpolation[mesh,par] depends linearly on par sometimes it would be useful to get a form ElementMeshInterpolation[mesh].par. For the simpe case of a 1D-mesh  ElementMeshInterpolation[mesh] would be a set of local triangle-functions, which could be evaluated once without knowing par. Jan 11, 2022 at 10:40
• hm, I am not sure I understand what you are looking for. Including a small example would go a long way - perhaps what you want is already there. Jan 11, 2022 at 11:34
• @user21: Setting the parameters to all permutations of {1,0,...,0 } would give the list of basisfunctions I'm interested . Later I'll add a small example Jan 11, 2022 at 11:36

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,

<< NDSolveFEM

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.

• (-1) I don't think this is a proper answer for this post, it should be posted as a separate question. And, the underlying request is essentially a duplicate of first request in Hugh's answer, isn't it? Then, related: mathematica.stackexchange.com/q/178847/1871 Jun 1, 2019 at 15:58
• @xzczd. I'm sorry that I cannot agree with you. The essence here is how to set boundary conditions in this case. And I cannot find any slight hint in the document. I mean it will be beneficial for the users to provide more illustrating examples, like in COMSOL, and would be further enhancements. Jun 2, 2019 at 1:03
• @xzczd, Thank you very much to point this post to me: mathematica.stackexchange.com/q/178847/1871. It is very similar to my questions. As you can see in my post that I have used the same tricks to mesh the region. However, the result is still not right . I cannot figure out the reason. Jun 2, 2019 at 13:33