I have 3D mesh and calculate the heat equation using NDSolve. I want to access the temperature values on the individual boundary elements to solve another problem (posted here).

I wrote some code that allows me to access the temperature from the interpolatingFunction from NDSolve on the individual boundary elements using an FVM approach (triangle element and it's centroid). I tested this approach by manually integrating over a boundary, but the solution differs too much from integrating using NIntegrate.

I think the error comes from FEM calculating on the nodes, rather than using the centroid. Is there a way to access the value of the boundary nodes (from the interpolatingFunction) and some kind of area of the node as a weight for integration?

A solution would be to calculate both integrals (manually and using NIntegrate) and calculate a correction factor that I can multiply with when using the individual values.

Find my approach below.

I am generating a mesh:

xmax = 2;
ymax = 1;
zmax = 1;

cubi = Cuboid[{0, 0, 0}, {xmax, ymax, zmax}];

mesh = DiscretizeRegion[cubi, MaxCellMeasure -> 0.0001];
bmesh = ToBoundaryMesh[mesh]
Graphics3D[mesh, Axes -> True, AxesLabel -> {"x", "y", "z"}]


Gathering the required information from the boundary mesh:

triList =  bmesh["BoundaryElements"][[1, 1]];
coordsList = bmesh["Coordinates"];

Creating a list of all the triangles that I am interested in (here all the ones on the side of the cube where y==0, see this post if you want to use boundary markers to do so):

y0Tris = {}; (*List of all tris with y=0*)

(*Find all triangle elements with all point's y coordinate = 0*)
 currTri = triList[[i]];
 yp1 = coordsList[[currTri[[1]]]][[2]];
 yp2 = coordsList[[currTri[[2]]]][[2]];
 yp3 = coordsList[[currTri[[3]]]][[2]];
 If[yp1 == 0 && yp2 == 0 && yp3 == 0, AppendTo[y0Tris, currTri]]
 , {i, Length[triList]}

Calculating the areas of the relevant triangles:

y0TrisAreas = Table[
   , {currTri, y0Tris}];

Calculating the centroids of the relevant triangles:

y0TrisCentroids = Table[
   , {currTri, y0Tris}];

Now I calculate the heat equation, setting a heat flux on the left and a temperature on the right:

ks = 1;
(*solve heat equation*)
TFun = NDSolveValue[{ks*\!\(\*SubsuperscriptBox[\(\[Del]\), \({x, y, z}\), \(2\)]\(T[x, y, z]\)\) == 0
      + NeumannValue[-100, x == 0], 
      DirichletCondition[T[x, y, z] == 0, x == xmax]},
   T, {x, y, z} \[Element] mesh];


Integrating over the boundary y==0, calculating the integral of the temperature and the average temperature:

areaBoundary = NIntegrate[Piecewise[{{1, y == 0}, {0, True}}], {x, y, z} \[Element] bmesh];

timeInt = First@Timing[
    integTemp = NIntegrate[
      Piecewise[{{TFun[x, y, z], y == 0}, {0, True}}], {x, y, z} \[Element] bmesh];(*integrating over whole area T(x,y,z)dA*)
    averageTemp = integTemp/areaBoundary; (*calculating average temperature on boundary*)
Print["Integrating time: " <> ToString@timeInt];
Print["Tavg*A: " <> ToString@integTemp];
Print["Tavg: " <> ToString@averageTemp];

Integrating time: 0.640625
Tavg*A: 200.
Tavg: 100.

And finally integrating manually by accessing the value of the temperature for each triangle element at it's centroid, multiplying with it's area and summing up over all triangles:

timeIntMan = First@Timing[
    (*performing manual integration*)
    integTempManual = Plus @@ Table[
       xcoord = y0TrisCentroids[[i]][[1]];
       ycoord = y0TrisCentroids[[i]][[2]];
       zcoord = y0TrisCentroids[[i]][[3]];
       TFun[xcoord, ycoord, zcoord]*
        y0TrisAreas[[i]](*calculating for each triangle T(centroid)*area_triangle*)
       , {i,Length[y0TrisCentroids]}];(*integrating over whole are T(x,y,z)dA*)
    averageTempManual = integTempManual/areaBoundary;(*calculating average temperature on boundary*)
Print["Integrating manual time: " <> ToString@timeIntMan];
Print["Tavg*A: " <> ToString@integTempManual];
Print["Tavg: " <> ToString@averageTempManual];

Integrating manual time: 0.03125
Tavg*A: 200.
Tavg: 100.

Now this works fine for this simple case. I set up a case that is a bit more complex. Heat flux on the left like the prior example, but temperature dependent Neumann BC on the front (y==0)

ks = 1;
Tamb = 20;
alphaFreeConvection = 5;

TFun2 = NDSolveValue[{ks*\!\(\*SubsuperscriptBox[\(\[Del]\), \({x, y, z}\), \(2\)]\(T[x, y, z]\)\) == 0
      + NeumannValue[-100, x == 0]
      + NeumannValue[(T[x, y, z] - Tamb)*alphaFreeConvection, y == 0]},
   T, {x, y, z} \[Element] mesh];


Performing the integration as above gives different results. Note, that the difference between the approaches becomes far more significant for more complex cases.

Integrating time: 0.640625
Tavg*A: 60.
Tavg: 30.

Integrating manual time: 0.03125
Tavg*A: 59.9901
Tavg: 29.9951

What I want to use this for, is calculating the temperature of the boundary (y==0, x-z-plane) as a function of x. Therefore I would take all the individual temperatures of the boundary elements, associate them to their x-coordinate (here centroid of the triangle) and throw them it into an interpolation function.

EDIT - using the function @Oleksii Semenov propsed:

I tried calculating the temperature at the centroids of the triangles using barycentric coordinates to interpolate from the values on the nodes.

This approach isn't succsessful either.

Find the code of the approach below:

timeIntMan2Bary = First@Timing[
   TemprArray = TFun2["ValuesOnGrid"];(*get values of boundary nodes*)
   (*loop over each triangle and sum up*)
   integTempManual2Bary = Plus @@ Table[
      currTri = y0Tris[[i]];
      (*creating 2d coordinates for barycentric coordinates, 
      deleting y coordinate,
      would need to map 3d triangles to 2d if using more complex boundaries*)
      pLis = {
        Delete[coordsList[[currTri]][[1]], 2],
        Delete[coordsList[[currTri]][[2]], 2],
        Delete[coordsList[[currTri]][[3]], 2]
      (*calculating centroid in 2d triangle coordinates, deleting y*)
      centroid = Delete[
        , 2];
      baryWeights = 
        centroid];(*barycentric weights*)
      nodeVals =Flatten@Partition[TemprArray[[Flatten[currTri]]], 
         3];(*temperatures at nodes, triangle corners*)
      tempInterp = 
       Plus @@ (nodeVals*
          baryWeights);(*interpolate temperature from barycentric weights*)
      tempInteg = tempInterp*y0TrisAreas[[i]];
      tempInteg(*to table*)
      , {i, Length[y0Tris]}
      ];(*end table*)
   ];(*end timing*)

averageTempManual2Bary = integTempManual2Bary/
  NIntegrate[Piecewise[{{1, y == 0}, {0, True}}], {x, y, z} \[Element] bmesh];

Print["Integrating manualBary time: " <> ToString@timeIntMan2Bary];
Print["Tavg*A: " <> ToString@integTempManual2Bary];
Print["Tavg: " <> ToString@averageTempManual2Bary];

Integrating manualBary time: 0.734375
Tavg*A: 60.1539
Tavg: 30.0769

EDIT - As @user21 mentioned, it seems to be unclear, what I want to do. Now I want to try explaining with equations.

I want to turn the temperature of the boundary of the solid (y=0) into a function of the x coordinate. eq1

To do so, I was thinking to divide my x range into small segments and do the following. eq2 Thats why I want to use the node/triangle element values of the temperature to stuff them into the segments.

With this I would have T_w for discrete x coordinates. Those I could use to generate an interpolation function or else.

With this I want to solve the coupling of a 1D region with a 3D region. The problem is to match the heat going into the channel with the heat going out of the wall. eq3 Where alpha is the heat transfer coefficient and Pi*D_c is the circumference of the channel. eq4


1 Answer 1


You can get solution values in nodes by

TemprArray = TFun[[4]]

{200., 200., 200., 200., 200., 200., 200., 200., 200., 200., 200., \
, 124.242, 127.273, 127.273, 130.303, 130.303, 130.303, \
124.242, 123.794, 133.333, 136.364}

It works in the case of stationary solution. For transient solutions there is another way.

As expected the number of mesh nodes coincide with length of TemprArray



Further we can get temperature values for every triangle posed on surface y=0 by

Partition[TemprArray[[Flatten[y0Tris]]], 3]  

Calculation of $T_w(x)$

Lets introduce additional 1D mesh on $0\leq x\leq x_{max}$ and calculate average (averaging along $0\leq z\leq z_{max}$) values of surface temperature in nodes of this mesh. By using Interpolation we can obtain function $T_w(x)$

n = 10;               (*number of elements in x-direction*)
h = xmax/n;           (*spetial step size*)  
Xmesh = Table[(i - 1)*h, {i, 1, n + 1}];
Tw[x_] := NIntegrate[TFun[x, 0, z], {z, 0, zmax}]/zmax
TwP = Map[Tw, Xmesh];
TwInt = Interpolation[Transpose[{Xmesh, TwP}]]
pic = Plot[TwInt[x], {x, 0, xmax}, ImageSize -> 300, 
  AxesLabel -> {"x", "Tw"}]

enter image description here

  • 2
    $\begingroup$ It's better to use TFun["ValuesOnGrid"] - directly using Part is fragile; for example if the data structure underneath will have to be changed for some reason. $\endgroup$
    – user21
    Oct 7, 2021 at 19:12
  • $\begingroup$ Thanks! I tried using the node values to interpolate the temperature of the triangle centroids with barycentric coordinates, but that doesn't really work either. I guess I am just lacking knowloedge of FEM. I will just use my first approach and than calculate a correction factor to match with the integral from NIntegrate. $\endgroup$
    – Tobias
    Oct 8, 2021 at 11:00
  • 1
    $\begingroup$ @Tobias, it's sill not clear in your posts what you actually want to compute; unless that becomes clear people will not be able to help. For example it is well known that the method of false transients is slow, but unless you write out the equations you want to solve I for one can not help. Simplify your problem. Concerning nodal values: it's almost never a good idea in FEM to use nodal values.... $\endgroup$
    – user21
    Oct 8, 2021 at 11:35
  • $\begingroup$ @user21 thanks, I added some equations to my question. Hope that helps to understand my problem. $\endgroup$
    – Tobias
    Oct 8, 2021 at 13:05
  • 1
    $\begingroup$ @Tobias You are welcome. To my mind we should avoid of using 3D interpolation function TFun(x,y,z) for numerical integration here. It is time consuming. I will try to post today workaround. By using structural FE mesh (in z direction only) we can reach computational savings. $\endgroup$ Oct 11, 2021 at 13:07

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.