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.

`ClearAll["Global`*"] SetDirectory[NotebookDirectory[]]; Needs["NDSolve`FEM`"]; 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}]}; pde // MatrixForm // TraditionalForm \[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"]; md["DegreesOfFreedom"] Show[mesh2["Wireframe"]] ClearAll[rhs] 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}}(*-{{Subscript[b, 1]},{Subscript[b, 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, 0}}](*{{\[PartialD]Subscript[b, 1]/\[PartialD]u,\[PartialD]Subscript[b, 1]/\[PartialD]v},{\[PartialD]Subscript[b, 2]/\[PartialD]u,\[PartialD]Subscript[b, 2]/\[PartialD]v}}*); 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] 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}}(*-{{Subscript[b, 1]},{Subscript[b, 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, y]}}](*{{\[PartialD]Subscript[b, \ 1]/\[PartialD]u,\[PartialD]Subscript[b, \ 1]/\[PartialD]v},{\[PartialD]Subscript[b, \ 2]/\[PartialD]u,\[PartialD]Subscript[b, 2]/\[PartialD]v}}*); `

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?