# FEM Mesh related errors

Bug introduced in 12.0. Fixed in 12.1.0.

I try to solve the simple heat transfer equation but there are some errors at the stage of the pre-processing.

Needs["NDSolveFEM"]
c = 1380*0.6 + 4200*0.4;
ρ = 1250.*0.6 + 1000*0.4;
λ = (0.111*0.6 + 0.56*0.4 );
Q = 3000000;
σ = 5.675*10^-8;
ϵ = 0.33;
T0 = 300.;
R0 = 0.00007;
rr = 0.0005;
h = 0.001;
MaxT=5;

Ω = ToElementMesh[Rectangle[{0., 0.}, {rr, h}],
MaxCellMeasure -> 0.5 10^-10, MeshQualityGoal -> 1];


The mesh is ElementMesh[{{0., 0.0005}, {0., 0.001}}, {QuadElement["<" 5000 ">"]}] with 5000 elements.

op = D[u[t, r, z],t] - λ/(c ρ)*(D[r^2 D[u[t, r, z], r], r]/r^2 +
D[u[t, r, z], {z, 2}]);
Γ = NeumannValue[
Piecewise[{{Q - ϵ σ (u[t, r, z]^4 - T0^4),
0 <= r <= R0}, {-ϵ σ (u[t, r, z]^4 - T0^4),
R0 < r <= rr}}], z == 0] +
NeumannValue[-ϵ σ (u[t, r, z]^4 - T0^4), z == h] +
NeumannValue[-ϵ σ (u[t, r, z]^4 - T0^4), r == rr];

fun=NDSolveValue[{op == Γ,
u[0, r, z] == T0}, u, {t, 0,
MaxT}, {r, z} ∈ Ω,
EvaluationMonitor :> (tp = t),
Method -> {"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement","IntegrationOrder" -> 4}}}, MaxSteps -> 100000,
MaxStepFraction -> 0.01]


Running the NDSolveValue leads to generation of a lot of errors:

Part::partw: Part {10154,10155,30390} of {} does not exist.

Part::partw: Part {10155,10156,30392} of {} does not exist.

Part::partw: Part {10156,10157,30394} of {} does not exist.

General::stop: Further output of Part::partw will be suppressed during this calculation.

Dot::dotsh: Tensors {NDSolveFEMFEMBoundaryConditionsDumpIntegratedShapeFunctionLookup[142]} and {{}[[{10154,10155,30390}]],{}[[{10155,10156,30392}]],{}[[{10156,10157,30394}]],{}[[{10157,10158,30396}]],{}[[{10158,10159,30398}]],{}[[{10159,10160,30400}]],{}[[{10160,10161,30402}]],{}[[{10161,10162,30404}]],{}[[{10162,10163,30406}]],{}[[{10163,10164,30408}]],{}[[{10164,10165,30410}]],{}[[{10165,10166,30412}]],{}[[{10166,10167,30414}]],<<26>>,{}[[{10193,10194,30468}]],{}[[{10194,10195,30470}]],{}[[{10195,10196,30472}]],{}[[{10196,10197,30474}]],{}[[{10197,10198,30476}]],{}[[{10198,10199,30478}]],{}[[{10199,10200,30480}]],{}[[{10200,10201,30482}]],{}[[{10201,10202,30484}]],{}[[{10202,10203,30486}]],{}[[{10203,10204,30488}]],<<92>>} have incompatible shapes.

CompiledFunction::cfta: Argument {NDSolveFEMFEMBoundaryConditionsDumpIntegratedShapeFunctionLookup[142]} at position 1 should be a rank 2 tensor of machine-size real numbers.

Transpose::nmtx: The first two levels of {{NDSolveFEMFEMBoundaryConditionsDumpIntegratedShapeFunctionLookup[142]}.{{}[[{10154,10155,30390}]],{}[[{10155,10156,30392}]],{}[[{10156,10157,30394}]],{}[[{10157,10158,30396}]],{}[[{10158,10159,30398}]],{}[[{10159,10160,30400}]],{}[[{10160,10161,30402}]],{}[[{10161,10162,30404}]],{}[[{10162,10163,30406}]],{}[[{10163,10164,30408}]],{}[[{10164,10165,30410}]],<<29>>,{}[[{10194,10195,30470}]],{}[[{10195,10196,30472}]],{}[[{10196,10197,30474}]],{}[[{10197,10198,30476}]],{}[[{10198,10199,30478}]],{}[[{10199,10200,30480}]],{}[[{10200,10201,30482}]],{}[[{10201,10202,30484}]],{}[[{10202,10203,30486}]],{}[[{10203,10204,30488}]],<<92>>}} cannot be transposed.

CompiledFunction::cfta: Argument Transpose[{{NDSolveFEMFEMBoundaryConditionsDumpIntegratedShapeFunctionLookup[142]}.{{}[[{10154,10155,30390}]],{}[[{10155,10156,30392}]],{}[[{10156,10157,30394}]],{}[[{10157,10158,30396}]],{}[[{10158,10159,30398}]],{}[[{10159,10160,30400}]],{}[[{10160,10161,30402}]],{}[[{10161,10162,30404}]],{}[[{10162,10163,30406}]],{}[[{10163,10164,30408}]],{}[[{10164,10165,30410}]],<<29>>,{}[[{10194,10195,30470}]],{}[[{10195,10196,30472}]],{}[[{10196,10197,30474}]],{}[[{10197,10198,30476}]],{}[[{10198,10199,30478}]],{}[[{10199,10200,30480}]],{}[[{10200,10201,30482}]],{}[[{10201,10202,30484}]],{}[[{10202,10203,30486}]],{}[[{10203,10204,30488}]],<<92>>}}] at position 1 should be a rank 2 tensor of machine-size real numbers.

Dot::dotsh: Tensors {NDSolveFEMFEMBoundaryConditionsDumpIntegratedShapeFunctionLookup[142]} and {{}[[{10154,10155,30390}]],{}[[{10155,10156,30392}]],{}[[{10156,10157,30394}]],{}[[{10157,10158,30396}]],{}[[{10158,10159,30398}]],{}[[{10159,10160,30400}]],{}[[{10160,10161,30402}]],{}[[{10161,10162,30404}]],{}[[{10162,10163,30406}]],{}[[{10163,10164,30408}]],{}[[{10164,10165,30410}]],{}[[{10165,10166,30412}]],{}[[{10166,10167,30414}]],<<26>>,{}[[{10193,10194,30468}]],{}[[{10194,10195,30470}]],{}[[{10195,10196,30472}]],{}[[{10196,10197,30474}]],{}[[{10197,10198,30476}]],{}[[{10198,10199,30478}]],{}[[{10199,10200,30480}]],{}[[{10200,10201,30482}]],{}[[{10201,10202,30484}]],{}[[{10202,10203,30486}]],{}[[{10203,10204,30488}]],<<92>>} have incompatible shapes.

Dot::dotsh: Tensors {NDSolveFEMFEMBoundaryConditionsDumpIntegratedShapeFunctionLookup[142]} and {{}[[{10154,10155,30390}]],{}[[{10155,10156,30392}]],{}[[{10156,10157,30394}]],{}[[{10157,10158,30396}]],{}[[{10158,10159,30398}]],{}[[{10159,10160,30400}]],{}[[{10160,10161,30402}]],{}[[{10161,10162,30404}]],{}[[{10162,10163,30406}]],{}[[{10163,10164,30408}]],{}[[{10164,10165,30410}]],{}[[{10165,10166,30412}]],{}[[{10166,10167,30414}]],<<26>>,{}[[{10193,10194,30468}]],{}[[{10194,10195,30470}]],{}[[{10195,10196,30472}]],{}[[{10196,10197,30474}]],{}[[{10197,10198,30476}]],{}[[{10198,10199,30478}]],{}[[{10199,10200,30480}]],{}[[{10200,10201,30482}]],{}[[{10201,10202,30484}]],{}[[{10202,10203,30486}]],{}[[{10203,10204,30488}]],<<92>>} have incompatible shapes.

General::stop: Further output of Dot::dotsh will be suppressed during this calculation.

CompiledFunction::cfta: Argument {{{{-1.87275*10^-8 (-8.1*10^9+Dot[<<2>>]^4)}},{{-1.87275*10^-8 (-8.1*10^9+Dot[<<2>>]^4)}},{{-1.87275*10^-8 (-8.1*10^9+Dot[<<2>>]^4)}}},{{{-1.87275*10^-8 (-8.1*10^9+Dot[<<2>>]^4)}},{{-1.87275*10^-8 (-8.1*10^9+Dot[<<2>>]^4)}},{{-1.87275*10^-8 (-8.1*10^9+Dot[<<2>>]^4)}}},{{{-1.87275*10^-8 (-8.1*10^9+Dot[<<2>>]^4)}},{{-1.87275*10^-8 (-8.1*10^9+Dot[<<2>>]^4)}},{{-1.87275*10^-8 (-8.1*10^9+Dot[<<2>>]^4)}}},<<46>>,{{{-1.87275*10^-8 (-8.1*10^9+Dot[<<2>>]^4)}},{{-1.87275*10^-8 (-8.1*10^9+Dot[<<2>>]^4)}},{{-1.87275*10^-8 (-8.1*10^9+Dot[<<2>>]^4)}}},<<92>>} at position 1 should be a rank 3 tensor of machine-size real numbers.

General::stop: Further output of CompiledFunction::cfta will be suppressed during this calculation.

AssembleMatrix::badmat: {7.04225*10^-6 MapThreadDot[{{{0.687298},{-0.0872983},{0.4}},{{-0.0872983},{0.687298},{0.4}},{{0.},{0.},{1.}}},{{{-1.87275*10^-8 Plus[<<2>>]}},{{-1.87275*10^-8 Plus[<<2>>]}},{{-1.87275*10^-8 Plus[<<2>>]}}}],7.04225*10^-6 MapThreadDot[{{{0.687298},{-0.0872983},{0.4}},{{-0.0872983},{0.687298},{0.4}},{{0.},{0.},{1.}}},{{{-1.87275*10^-8 Plus[<<2>>]}},{{-1.87275*10^-8 Plus[<<2>>]}},{{-1.87275*10^-8 Plus[<<2>>]}}}],<<48>>,<<92>>} is not a valid element matrix to be assembled.

General::stop: Further output of AssembleMatrix::badmat will be suppressed during this calculation.

Transpose::nmtx: The first two levels of {{NDSolveFEMFEMBoundaryConditionsDumpIntegratedShapeFunctionLookup[71]}.{{}[[{286,143,10722}]],{}[[{429,286,11008}]],{}[[{572,429,11293}]],{}[[{715,572,11578}]],{}[[{858,715,11863}]],{}[[{1001,858,12148}]],{}[[{1144,1001,12433}]],{}[[{1287,1144,12718}]],{}[[{1430,1287,13003}]],{}[[{1573,1430,13288}]],{}[[{1716,1573,13573}]],{}[[{1859,1716,13858}]],{}[[{2002,1859,14143}]],<<26>>,{}[[{5863,5720,21838}]],{}[[{6006,5863,22123}]],{}[[{6149,6006,22408}]],{}[[{6292,6149,22693}]],{}[[{6435,6292,22978}]],{}[[{6578,6435,23263}]],{}[[{6721,6578,23548}]],{}[[{6864,6721,23833}]],{}[[{7007,6864,24118}]],{}[[{7150,7007,24403}]],{}[[{7293,7150,24688}]],<<21>>}} cannot be transposed.

Transpose::nmtx: The first two levels of {{NDSolveFEMFEMBoundaryConditionsDumpIntegratedShapeFunctionLookup[71]}.{{}[[{1,144,10297}]],{}[[{144,287,10724}]],{}[[{287,430,11009}]],{}[[{430,573,11294}]],{}[[{573,716,11579}]],{}[[{716,859,11864}]],{}[[{859,1002,12149}]],{}[[{1002,1145,12434}]],{}[[{1145,1288,12719}]],{}[[{1288,1431,13004}]],{}[[{1431,1574,13289}]],{}[[{1574,1717,13574}]],<<27>>,{}[[{5578,5721,21554}]],{}[[{5721,5864,21839}]],{}[[{5864,6007,22124}]],{}[[{6007,6150,22409}]],{}[[{6150,6293,22694}]],{}[[{6293,6436,22979}]],{}[[{6436,6579,23264}]],{}[[{6579,6722,23549}]],{}[[{6722,6865,23834}]],{}[[{6865,7008,24119}]],{}[[{7008,7151,24404}]],<<21>>}} cannot be transposed.

General::stop: Further output of Transpose::nmtx will be suppressed during this calculation.

NDSolveValue::indexss: The DAE solver failed at t = 0.. The solver is intended for index 1 DAE systems and structural analysis indicates that the DAE is structurally singular.


As I understood, the initial errors occur because the solver tries to find the extra mesh elements above the 5000 of the basically produced. And the following errors are derived from these..

Can anybody help with a proper problem initialization?

The solver finally produces certain solution but I'm not sure that it is correct..

• Do you have a typo there? You did not normalize the heat flux by $\rho c$? Commented Mar 14, 2020 at 15:09
• @AlexTrounev, Thank you, there are missed denominators in BC. However, this does not change anything regarding the technical problems with solution. Commented Mar 16, 2020 at 4:36
• You can use linear FEM with my code. Commented Mar 16, 2020 at 11:31
• Version 12.1 has a tutorial on Heat Transfer Modeling Commented Mar 20, 2020 at 19:49

This is a bug with time dependent nonlinear NeumannValue introduced in Version 12.0. This is fixed in version 12.1 which is hopefully coming in the not to distant future at the time of this writing. There is no known top level workaround. I apologize for the inconvenience. If this is very important for you, I cloud try to implement this with the low level FEM functions - I have not tried this so I can no guarantee that this will work. Let me know what you think.

Update:

Here is a low level code.

Needs["NDSolveFEM"]
c = 1380*0.6 + 4200*0.4;
ρ = 1250.*0.6 + 1000*0.4;
λ = (0.111*0.6 + 0.56*0.4);
Q = 3000000;
σ = 5.675*10^-8;
ϵ = 0.33;
T0 = 300.;
R0 = 0.00007;
rr = 0.0005;
h = 0.001;
MaxT = 5;

mesh = ToElementMesh[Rectangle[{0., 0.}, {rr, h}],
MaxCellMeasure -> 10 10^-10, MeshQualityGoal -> 1,
"MeshElementType" -> "TriangleElement", "MeshOrder" -> 1]


Note that I have substantially coarsened the mesh. Try to fix the remaining issues first before switching back to a finer mesh.

Set up the PDE without the boundary conditions:

op = D[u[t, r, z],
t] - λ/(c ρ)*(D[r^2 D[u[t, r, z], r], r]/r^2 +
D[u[t, r, z], {z, 2}]);


The BCs have some issues: Split the the Piecewise into two BCs. For some reason the second BC introduces a real slow convergence. You'd need to experiment a bit with that.

Γ = {(*NeumannValue[
Piecewise[{{Q-ϵ σ (u[t,r,z]^4-T0^4),
0\[LessEqual]r\[LessEqual]R0},{-ϵ σ (u[t,r,z]^4-
T0^4),R0<r\[LessEqual]rr}}],
z\[Equal]0],*)(*NeumannValue[-ϵ σ (u[t,r,z]^4-
T0^4),z\[Equal]h],*)
NeumannValue[-ϵ σ (u[t, r, z]^4 - T0^4), r == rr]};


Set up FEM data:

{sdpde, sdbc, vd, sd, methodData} =
NDSolveFEMProcessPDEEquations[{op == 0, u[0, r, z] == T0},
u, {t, 0, MaxT}, {r, z} \[Element] mesh];


Now, we initialize the BCs separately:

initBCs = InitializeBoundaryConditions[vd, sd, {Γ}];
sbcs = DiscretizeBoundaryConditions[initBCs, methodData, sd];


Set up a helper function to apply the BCs during time integration:

discretizePDEResidual[t_?NumericQ, u_?VectorQ, dudt_?VectorQ] :=

Module[{l, s, d, tdpde, tbcs, nldpde, nlbcs, sdTemp},

NDSolveSetSolutionDataComponent[sd, "Time", t];
NDSolveSetSolutionDataComponent[sd, "DependentVariables", u];

s = sdpde["StiffnessMatrix"];
d = sdpde["DampingMatrix"];

tbcs = DiscretizeBoundaryConditions[initBCs, methodData, sd,
"Transient"];

nlbcs =
DiscretizeBoundaryConditions[initBCs, methodData, sd,
"Nonlinear"];

DeployBoundaryConditions[{l, s, d}, nlbcs];
DeployBoundaryConditions[{l, s, d}, tbcs];
DeployBoundaryConditions[{l, s, d}, sbcs];

d.dudt + s.u - l
]


Set up the initial conditions and the sparsity pattern:

initT0 = T0 & /@ mesh["Coordinates"];
sparsity = sdpde["DampingMatrix"]["PatternArray"];


Do the time integration. Because there were may step rejections (with the other BCs) I added the "IDA" option to reduce those. This might not be necessary once the NeumannValue issues are understood.

Monitor[tufun =
NDSolveValue[{discretizePDEResidual[t, u[t], u'[ t]] == 0,
u[0] == initT0}, u, {t, 0, MaxT}
, Method -> {
"TimeIntegration" -> {"IDA", "MaxDifferenceOrder" -> 2}
, "EquationSimplification" -> "Residual"}
, Jacobian -> {Automatic, Sparse -> sparsity}
, EvaluationMonitor :> (monitor = Row[{"t = ", CForm[t]}])
(*,MaxStepFraction\[Rule]0.01*)
], monitor]


From the result you can then construct an interpolating function with:

ufun = ElementMeshInterpolation[{tufun["Coordinates"][[1]],
methodData["ElementMesh"]}, Partition[tufun["ValuesOnGrid"], 1]]

• Thank you for your example of coding of my task. The BC with piecewise is of interest here (localized heating) but the solutions does not run with it. I've decreased the heating intensity Q for 4-5 orders of magnitude and it starts solving but I need namely such big values. Increasing the mesh density does not help.. So, something is not good with this solution.. What means these f <-> b jumps of t at the beginning of calculations? Commented Mar 14, 2020 at 7:53
• How can we define the initial step in time? It always try dance around 10^-5 making attempts with shorter steps again and again.. Commented Mar 14, 2020 at 8:03
• @user21 In what cases is this bug detected? If, for example, NeumannValue is linear and time independent? Commented Mar 14, 2020 at 12:22
• @AlexTrounev, the issue comes up if the NeumannValue is nonlinear and time depdenent. There is no is no issue for the stationary case (= the time independent case) and there is no issue for a linear time dependent NeumannValue Commented Mar 16, 2020 at 6:16
• @Rom38, try ref/StartingStepSize. You can also try to set initT0 to something else. This will set the initial condition. Another idea is to use a smooth transition for the Piecewise NeumannValue, just to see if that works. Commented Mar 16, 2020 at 6:22

We can use a linear FEM solver. It is necessary to normalize u by T0, Q and $$\sigma u^4$$ by $$\rho c T0$$. Then the code is

Needs["NDSolveFEM"]
c = 1380*0.6 + 4200*0.4;
\[Rho] = 1250.*0.6 + 1000*0.4;
\[Lambda] = (0.111*0.6 + 0.56*0.4); T0 = 300.;
Q = 3000000/(c \[Rho] T0);
qr = 0.33 T0^3 (5.675*10^-8)/(c \[Rho] );
k = \[Lambda]/(c \[Rho] );

R0 = 0.00007;
rr = 0.0005;
h = 0.001;
tm = 5 ; tau = 1/20; nmax = Round[tm/tau];
\[CapitalOmega] =
ToElementMesh[Rectangle[{0., 0.}, {rr, h}],
MeshRefinementFunction ->
Function[{vertices, area},
area > 0.0000005 (0.00001 + 10 Norm[Mean[vertices]])]];
\[CapitalOmega]["Wireframe"]
U[0][r_, z_] := 1; Do[
U[i] = NDSolveValue[{(u[r, z] - U[i - 1][r, z])/tau -
k (D[u[r, z], r, r] + 2 D[u[r, z], r]/r +
D[u[r, z], {z, 2}]) == (NeumannValue[
Q - qr  (U[i - 1][r, z]^3 u[r, z] - 1),
z == 0 && 0 < r <= R0] +
NeumannValue[-qr (U[i - 1][r, z]^3 u[r, z] - 1),
z == 0 && R0 < r <= rr] +
NeumannValue[-qr  (U[i - 1][r, z]^3 u[r, z] - 1), z == h] +
NeumannValue[-qr (U[i - 1][r, z]^3 u[r, z] - 1), r == rr])},
u, {r, z} \[Element] \[CapitalOmega]];, {i, 1, nmax}];


Temperature visualization

{DensityPlot[T0 U[nmax][r, z], {r, z} \[Element] \[CapitalOmega],
ColorFunction -> "TemperatureMap", PlotLegends -> Automatic,
PlotRange -> All, FrameLabel -> Automatic],
DensityPlot[T0 U[nmax][r, z], {r, z} \[Element] \[CapitalOmega],
ColorFunction -> "TemperatureMap", PlotLegends -> Automatic,
FrameLabel -> Automatic], ListLinePlot[Table[{i tau 1., T0 U[i][.0, .0]}, {i, 0, nmax}],
AxesLabel -> {t, T}, PlotRange -> All]}
`

• Thanks, this is good approach but I need solve the equation with strong inner heat sources. This makes the temporal dependence highly non-linear.. I'll try your workaround to see how does it work. Commented Mar 16, 2020 at 4:38
• @Rom38, try to use Alex's answer as an initial condition for the equation. Commented Mar 16, 2020 at 6:17
• @Rom38 What is strong? Now it is 3,000,000. Commented Mar 16, 2020 at 11:35