3D FEM Mesh – Accessing values of individual boundary elements

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:

Needs["NDSolve`FEM`"] 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*) Do[  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[    Area[Triangle[      {coordsList[[currTri[[1]]]],       coordsList[[currTri[[2]]]],       coordsList[[currTri[[3]]]]}      ]]    , {currTri, y0Tris}]; 

Calculating the centroids of the relevant triangles:

y0TrisCentroids = Table[    RegionCentroid[Triangle[      {coordsList[[currTri[[1]]]],       coordsList[[currTri[[2]]]],       coordsList[[currTri[[3]]]]}      ]]    , {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.