I'm trying to model a 2D nonlinear anisotropic heat transfer using FEM with element markers. I got stuck with errors and have no idea what to do to make it work...

I have simple toy-model: two connected squares with bottom side of the right heated through a constant heat flux and a left side of a left square is kept at constant temperature. Each of the squares has different material properties. From the FEM perspective they are distinguished through ElementMarkers. I've built this model using procedure described in the Heat Transfer tutorial (PDEModels/tutorial/HeatTransfer/HeatTransfer in Mathematica help) in "Modeling Heat Sources using Element Markers" section. Heat transfer FEM solving equations, as well as all the methods described below, come from the same Heat Transfer tutorial. I'm using Mathematica 12.1.

According to the Heat Transfer tutorial simulating anisotropic heat transfer is as easy as inputing k as a matrix rather than a single value. It works well with constants (you just need to remember to define k separately and not in "parameters" replace rules' set).

Again, according to Heat Transfer tutorial simulating nonlinear heat transfer (changing properties of a material with respect to its temperature) is as easy as defining the property as a function depending on temperature. And indeed it works well.

Finally, to give different properties to different parts of the model you just need a piecewise function with ElementMarkers. Again, on its own it works great.

The problem comes when combining all three methods. I tried combinations of two of the three methods and they all work - you can do nonlinear anisotropic without element markers (same properties for the whole simulation), nonlinear with markers, and anisotropic with markers. But when you try all three it gives errors.

Here is the code for the whole simulation and errors (errors as a picture as well as it is nicer to look at image than the code). Can you give me a clue what to do?


HeatTransferModel[T_, X_List, k_, \[Rho]_, Cp_, Velocity_, Source_] :=
  Module[{V, Q, a = k},
  V = If[Velocity === "NoFlow", 
    0, \[Rho]*Cp*Velocity.Inactive[Grad][T, X]];
  Q = If[Source === "NoSource", 0, Source];
  If[ FreeQ[a, _?VectorQ], a = a*IdentityMatrix[Length[X]]];
  If[ VectorQ[a], a = DiagonalMatrix[a]];
  (* Note the - sign in the operator *)
  a = PiecewiseExpand[Piecewise[{{-a, True}}]];
  Inactive[Div][a.Inactive[Grad][T, X], X] + V - Q]

coordinates = {
   {0, 0}, {0, 1}, {1, 1}, {1, 0},
   {1, 0}, {1, 1}, {2, 1}, {2, 0}
el = LineElement[{{1, 2}, {2, 3}, {3, 4}, {4, 1}, {5, 6}, {6, 7}, {7, 
     8}, {8, 5}}];
bmesh = ToBoundaryMesh["Coordinates" -> coordinates, 
   "BoundaryElements" -> {el}];

mesh = ToElementMesh[bmesh, 
   "RegionMarker" -> {{{0.5, 0.5}, 1}, {{1.5, 0.5}, 2}}];
  "MeshElementStyle" -> {Directive[FaceForm[Red]], 

\[CapitalGamma]temp = DirichletCondition[T[x, y] == 300, x == 0];
\[CapitalGamma]heater = NeumannValue[1000, y == 0 && 1 <= x <= 2];

k2[T_] := 600 - T;

parameters = {
   \[Rho] -> 
    Piecewise[{{1, ElementMarker == 1}, {2, ElementMarker == 2}}],
   Cp -> Piecewise[{{3, ElementMarker == 1}, {4, ElementMarker == 2}}]
kSet = Piecewise[{{{{5, 0}, {0, 50}}, 
     ElementMarker == 1}, {{{k2[T[x, y]], 0}, {0, 6}}, 
     ElementMarker == 2}}];

pde = {HeatTransferModel[T[x, y], {x, y}, kSet, \[Rho], Cp, "NoFlow", 
      "NoSource"] == \[CapitalGamma]heater, \[CapitalGamma]temp} /. 
Tfun = NDSolveValue[pde, T, {x, y} \[Element] mesh]


InitializePDECoefficients::femcnmd: The PDE coefficient {-(-(\[Piecewise]   {<<2>>} Equal[<<2>>]
{<<2>>} Equal[<<2>>]
0   True

)).{(T^(1,0))[x,y],(T^(0,1))[x,y]},-(-(\[Piecewise] {<<2>>} Equal[<<2>>]
{<<2>>} Equal[<<2>>]
0   True

)).{(T^(1,0))[x,y],(T^(0,1))[x,y]}} does not evaluate to a numeric matrix of dimensions {2,1} at the coordinate {0.295534,0.46115}; it evaluated to {{-0.5,-5.},{-0.5,-5.}} instead.

NDSolveValue::fempslf: The linearization process in PDESolve failed.

1 Answer 1


There seems to be a problem with the linearization of the PDE when it contains a Piecewise.

Try to use

HeatTransferModel[T_, X_List, k_, \[Rho]_, Cp_, Velocity_, Source_] :=
  Module[{V, Q, a = k}, 
  V = If[Velocity === "NoFlow", 
    0, \[Rho]*Cp*Velocity . Inactive[Grad][T, X]];
  Q = If[Source === "NoSource", 0, Source];
  If[FreeQ[a, _?VectorQ], a = a*IdentityMatrix[Length[X]]];
  If[VectorQ[a], a = DiagonalMatrix[a]];
  (*Note the-sign in the operator*)(*a=PiecewiseExpand[Piecewise[{{-a,
  Inactive[Div][-a . Inactive[Grad][T, X], X] + V - Q]

with this kSet

kSet = {{If[ElementMarker == 1, 5, k2[T[x, y]] // Evaluate], 0}, {0, 
   If[ElementMarker == 1, 50, 6]}}

This works in V12.3 and V12.2 - I no longer have 12.1 installed.

  • $\begingroup$ Unfortunately neither If nor Piecewise work. It still gives the same error. Only the form of PDE coefficient in the error changed, everything else is the same. But I must say I didn't think of this on my own ;) $\endgroup$
    – mgolunski
    Commented Sep 17, 2021 at 11:22
  • $\begingroup$ @m_i-k_i, see update. I think you need to comment out the PiecewiseEpand too. $\endgroup$
    – user21
    Commented Sep 17, 2021 at 12:30
  • $\begingroup$ It worked. Thank you. For anyone looking at this - do not forget about the minus sign before a in last Div! (or replace PiecewiseExpand line with a = -a). The solution works also with Piecewise instead of If. $\endgroup$
    – mgolunski
    Commented Sep 17, 2021 at 13:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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