What is the best method for interpolating a 3d array as an equation for heat generation source over a mesh?

variational formulation
1

I’m not sure what the best method for completing this process, so I apologize if my attempt is wildly inaccurate.

The issue is as follows:

  • A cubic mesh receives a heat load in the form of a volume flux
  • This flux is derived from a dataset (simply a 3D array, where at every point in the mesh, a discrete heat load is determined
  • The number of cells in the 3D array is not equal to those in the mesh, but the relative dimension are all equal (they are both cubes)

How might one go about creating either an expression or simply interpolating the array and putting it in such a form such that the weak form solver will recognize the appropriate heat load at each point in space?

Here are code samples to give a better idea of what the question is:

# Heat Structure
patt2D = [[0, 0, 1, 2, 3, 3],
          [0, 0, 1, 2, 2, 3],
          [0, 0, 0, 1, 2, 2],
          [0, 0, 0, 0, 1, 1],
          [0, 0, 0, 0, 0, 0],
          [0, 0, 0, 0, 0, 0]]

empty = [[0 for _ in range(6)] for _ in range(6)]
thickness = 2
patt3D = np.array([patt2D if i < thickness else empty for i in range(6)])

from scipy.interpolate import RegularGridInterpolator as rgi

x = np.arange(6)
y = np.arange(6)
z = np.arange(6)

interp = rgi((x, y, z), patt3D)
Vq = ???

# Weak form 
F = (
    fn.inner(k * fn.grad(u), fn.grad(v)) * (dx)  # Base FEM terms over the entire volume
    + (sc.sigma * epsilon * (q(u_0) - T_s4) * v) * (ds) # Radiation term for all surfaces
    + r * (u - T_s) * v * (ds)  # Convection term for all surfaces
    - Vq * v * (dx) # Heat generation in volume
    )

a, L = fn.lhs(F), fn.rhs(F)

# Write solution to file
file = fn.File("Solutions/output_heat_gen_cube.pvd")

u = fn.Function(V)     # new unknown function

last_eps = 0
while eps > tol and iter < maxiter:
    iter += itermod
    fn.solve(a == L, u)
    diff = np.array(u.vector()) - np.array(u_0.vector())
    eps = np.linalg.norm(diff, ord=np.Inf)
    print('iter=%d: norm=%g' % (iter, eps))

    u_0.assign(u)   # update for next iteration

    last_eps = eps
    file << u, iter

This allowed some level of interpolation of the points, however, it wasn’t completely successful in covering the entire mesh in some value (or approximation of the value) that can be recognized as a heat load.

Is there a way of doing an expression fit of sorts instead, such that it can approximate that heat load with an expression?

Thank you in advance for any help.

2

See for instance: How to impose a pressure distribution from a data file as a boundary condition (loading)? - #5 by dokken

You can tabulate the dof coordinates of your space, and figure out what position in the 3D data-set each entry has.

3

Thanks for your reply, Dokken. Would the data on the mesh be continuous wit this technique? I haven’t been able to convert that 2D example to 3D yet, but I’m beginning to wonder if it’s actually what I want/need.

Is there no way to generate a continuous expression estimate from the data in to “lay over” the mesh afterwards?

4

If you use a continuous space, like first order Lagrange, you will get the dof coordinates of the degrees of freedom (at mesh vertices) and assign data to those, then get a linear variation between mesh nodes.

5

I have attempted to implement the 3D interpretation of what you sent:


from scipy.interpolate import RegularGridInterpolator as rgi
import matplotlib.pyplot as plt

data = np.array(patt3D.flatten())
coords = np.array([[x, y, z] for x in range(6) for y in range(6) for z in range(6)])
ind2 = np.lexsort((coords[:, 2], coords[:, 1], coords[:, 0]))

p = fn.Function(V)
dof_coords = V.tabulate_dof_coordinates()
ind1 = np.lexsort((dof_coords[:, 0], dof_coords[:, 1], dof_coords[:, 2]))

for (i2, i1) in zip(ind2, ind1):
    p.vector()[i1] = data[i2]

This now follows the creation of patt3D and precedes the weak form.

It seems to not correctly find the corresponding index location, and it doesn’t map over the whole mesh. Do you have a solution to this issue? I also don’t see any sort of interpolation occurring with this method, as in, there will be “unfluxed” facets between patches of what should be a continuous application of flux using this method, it seems? I could be wrong, however, it clearly isn’t working as intended right now.

EDIT: I have used RegularGridInterpolator instead (as I had intended from the beginning), however it seems to be somewhat effective now.


from scipy.interpolate import RegularGridInterpolator as rgi
import matplotlib.pyplot as plt 

x, y, z = np.arange(0, 6), np.arange(0, 6), np.arange(0, 6)
interp = rgi((x, y, z), patt3D.T, "linear")

qq = fn.Function(V)
q_dof_coords = V.tabulate_dof_coordinates()

for i, coord in enumerate(q_dof_coords):
    qq.vector()[i] = interp(coord)

This – under the conditions that the mesh and the patt3D have the same dimensions (6x6x6pts in this case) and I don’t iterate it past the first solution (my radiative boundary uses Picard iteration for the T^4 nonlinearity) – has garnered some success. These limitations are pretty problematic still as iteration never stabilizes, but rather blows up.

Do you have any insights on what could be causing this instability?