I work on FEM in Mathematica and want to speed up the computation when dealing with very fined mesh (i.e., a huge number of DOF).

I found that a very much longer time is consumed when the Coefficients of PDE contains interpolation functions, compared to the case that the Coefficients are constants. To be more specific, let consider the following exemplary code.



width = 1.0; heigth = 2.0; ydatum = 0.0;
\[ScriptCapitalR] = 
  ImplicitRegion[0 <= x <= width && ydatum <= y <= heigth + ydatum , {x, y}];
RegionPlot[\[ScriptCapitalR], AspectRatio -> Automatic, ImageSize -> Small]

pde = {Div[{{\[Lambda] + 2 \[Mu], 0}, {0, \[Mu]}}.Grad[
       u[x, y], {x, y}] + {{0, \[Lambda]}, {\[Mu], 0}}.Grad[
       v[x, y], {x, y}], {x, y}], 
   Div[{{0, \[Mu]}, {\[Lambda], 0}}.Grad[
       u[x, y], {x, y}] + {{\[Mu], 0}, {0, \[Lambda] + 2 \[Mu]}}.Grad[
       v[x, y], {x, y}], {x, y}]};

\[CapitalGamma]g = {DirichletCondition[u[x, y] == 0, x == 0], 
       DirichletCondition[v[x, y] == 0, y == ydatum], 
       DirichletCondition[v[x, y] == 0.2, y == heigth + ydatum]};

\[Lambda] = 1.0; \[Mu] = 1.0;

{dPDE, dBC, vd, sd, md} = 
  ProcessPDEEquations[{pde == {0, 0}, \[CapitalGamma]g}, {u, 
    v}, {x, y} \[Element] \[ScriptCapitalR], 
   Method -> {"FiniteElement", 
     "MeshOptions" -> {"ImproveBoundaryPosition" -> False, 
       "MaxCellMeasure" -> 0.0001, "MeshOrder" -> 2}}];
mesh2 = md["ElementMesh"];
linearLoad = dPDE["LoadVector"];
linearStiffness = dPDE["StiffnessMatrix"];
diriPos = dBC["DirichletRows"];
offsets = md["IncidentOffsets"];


tol = 10^-8;
err = 1.0;
iter = 0;
\[Alpha] = 1.0;
rhs[uIn_] := Module[{uOld}, uOld = uIn;
  While[err >= tol && iter <= 20,
   u0 = ElementMeshInterpolation[{mesh2}, 
     uOld[[offsets[[1]] + 1 ;; offsets[[2]]]]];
   v0 = ElementMeshInterpolation[{mesh2}, 
     uOld[[offsets[[2]] + 1 ;; offsets[[3]]]]];
   nlPdeCoeff = InitializePDECoefficients[vd, sd,
     "LoadCoefficients" -> -{{0}, {-1}},
     "LoadDerivativeCoefficients" -> {{{(\[Lambda] + 2 \[Mu]) D[
            u0[x, y], x] + \[Lambda] D[v0[x, y], 
            y], \[Mu] (D[u0[x, y], y] + 
            D[v0[x, y], x])}}, {{\[Mu] (D[u0[x, y], y] + 
            D[v0[x, y], x]), \[Lambda] D[u0[x, y], 
            x] + (\[Lambda] + 2 \[Mu]) D[v0[x, y], y]}}},
     "ReactionCoefficients" -> {{0, 0}, {0, 

   t1 = AbsoluteTime[];
   nlsys = DiscretizePDE[nlPdeCoeff, md, sd];(*time ~s*)
   time1 = AbsoluteTime[] - t1;

   nlLoad = nlsys["LoadVector"];
   nlStiffness = nlsys["StiffnessMatrix"];
   ns = nlStiffness + linearStiffness;
   nl = nlLoad + linearLoad;
   nl[[diriPos]] = {0.};
   ns[[diriPos, All]] = 0.;
   ns[[All, diriPos]] = 0.;
   (ns[[#, #]] = 1.) & /@ diriPos;
   dU = LinearSolve[ns, nl];
   Print[iter, " Residual: ", Norm[nl, Infinity], "  Correction: ", 
    Norm[dU, Infinity], "  Time: ", time1];
   err = Norm[nl, Infinity]; iter++;
   uOld = uOld + \[Alpha]*dU;];

uOld = ConstantArray[{0.}, md["DegreesOfFreedom"]];
uOld[[diriPos]] = dBC["DirichletValues"];

uNew = rhs[uOld];

The "LoadDerivativeCoefficients" already contain the interpolated function u0 and v0, the other Coefficients are constants. After timing each command during running, I found the DiscretizePDE is most time consuming. For example, in this case, the DiscretizePDE command will cost ~6 seconds during each iteration.

On the other hand, if we set InitializePDECoefficients as:

nlPdeCoeff = InitializePDECoefficients[vd, sd,
   "LoadCoefficients" -> -{{0}, {-v0[x, y]^2}},
   "LoadDerivativeCoefficients" -> {{{(\[Lambda] + 2 \[Mu]) D[
          u0[x, y], x] + \[Lambda] D[v0[x, y], 
          y], \[Mu] (D[u0[x, y], y] + 
          D[v0[x, y], x])}}, {{\[Mu] (D[u0[x, y], y] + 
          D[v0[x, y], x]), \[Lambda] D[u0[x, y], 
          x] + (\[Lambda] + 2 \[Mu]) D[v0[x, y], y]}}},
   "ReactionCoefficients" -> {{0, 
      0}, {0, -2 v0[x, 

Now the time of running DiscretizePDE is ~10 seconds. If we refine the mesh, the consumed time will increase. (The above code is running with 16 cores.)

Since I need to run a lot of such code. Is there a way to cut of the time that the DiscretizedPDE needs when InitializePDECoefficients contains Interpolation functions?

  • $\begingroup$ There are several issues. The compiled FEM code needs to call the Mathematica Kernel to evaluate the interpolation function. This is slow. And also the interpolation code must find the element which you are integrating over. There is not much that can be done until a proper nonlinear FEM solver is implemented. Then one could also think a bit about how the evaluation over an interpolation function could be made faster if the mes in the interpolation function is the same as in the DiscretizePDE $\endgroup$ – user21 Nov 14 '17 at 3:33

DiscretizePDE should be the part where all the matrices get assembled. That is not entirely for free, in particular for finite elements of order two. Apart from the linear solves, this is where the numbers have to be crunched. If the coeffients are not constant, quadrature rules have to be used in order to compute the matrix entries accurately.

Still, here are some general suggestions:

  • If you do not need that much accuracy, you might consider reducing the order to 1.

  • Is it even acceptable to coarsen the mesh? Maybe you can create a mesh with fewer elements by better adapting it to the geometry and to properties of the solution? (There are some examples in the documentation that show you how that is possible.)

  • Since you already dug quite deeply into the internal structure of Mathematica's FEM and since you are going to run many solves: Can you reuse some of the generated data such as the mesh? Mass, stiffness, or reaction matrices?

I ran your code with mesh order 1 and I have to agree that this is ridiculously slow (~2.5 sec on a quad core; seemingly this is also not well parallelized) for a mesh with less than 32000 triangles. The "LoadDerivativeCoefficients" is the bottleneck here. In the worst case, Mathematica runs an instance of NIntegrate over each triangle.

Note that an InterpolatingFunction needs a certain overhead to determine in which element its argument lies. If you know what the degrees of freedom mean, you can try to assemble the "LoadDerivativeCoefficients" by hand with some 3- or 4-point Gauss quadrature (which you can find here). Doing that elementwise with the coefficients from uOld circumvents the lookup that InterpolatingFunctionss have to perform. Really, ElementMeshInterpolation is only useful for scarce evaluation, e.g., for getting a plot in the end.

  • $\begingroup$ Please see comment above. $\endgroup$ – user21 Nov 14 '17 at 3:33
  • $\begingroup$ Thank you. I try to write a FEM code in Mathematica from scratch and see if this helps in saving computational cost. $\endgroup$ – Wilhelm Jan 16 '18 at 21:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.