How to Apply Piecewise Function In Weak Form

variational formulation
1

Hi, I have a linearly interpolated set of points that I would like to apply to a mesh. Ideally, I could use an array to produce hard boundaries between areas with and without volumetric heat flux. The effect would be completely identical to defining boundaries with functions such as:

class Top(fn.SubDomain):
    def inside(self, x, on_boundary):
        tol = 1e-10
        return on_boundary and fn.near(x[2], h, tol)

etc.

In my current iteration, I have:

import fenics as fn
import numpy as np
import scipy.constants as sc
from scipy.interpolate import RegularGridInterpolator as rgi
import matplotlib.pyplot as plt 


# Parameters
# in SI
# Material and environmental constants for Quartz and the laser setup
B = 0.02  # x-dimension of the box in meters
l = 0.02 # y-dimension of the box in meters
h = 0.01  # z-dimension of the box in meters (height of the material)
w = h/10 # half-thickness of inner heated layer
nB = nl = 20 # number of facets per meter in x and z dimensions
nh = 10 # number of facets per meter in y dimension

# Material and surrounding properties
k = 1 # Thermal conductivity in W/m·K

# Mesh Creation
mesh = fn.BoxMesh(fn.Point(0, 0, 0), fn.Point(B, l, h), int(nB), int(nl), int(nh))

# Function Space Initialization
V = fn.FunctionSpace(mesh, "CG", 2)

# Heat Structure
p_density = 1000 # watts per cubic meter
p = p_density / (B*l*2*w)

qq2 = fn.Constant(p)

pattx = patty = pattz = 200
patt3D = np.zeros((pattx, patty, pattz))
heatrange_min, heatrange_max = int(pattz*2/5), int(pattz*3/5)
patt3D[:,:,heatrange_min:heatrange_max] = p

x = np.linspace(0, B, pattx)
y = np.linspace(0, l, patty)
z = np.linspace(0, h, pattz)
interp = rgi((x, y, z), patt3D, "nearest")

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

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

# Visualize Vertical Flux Distribution
n = 200
x_vals = np.linspace(0, h, n)
x_vals2 = np.linspace(0, h, pattz)
flux_vals = [qq((0, 0, xi)) for xi in x_vals]

fig = plt.figure()
ax1 = fig.add_subplot(111)

ax1.plot(x_vals, flux_vals, color = "red", linestyle="dotted")
ax1.plot(x_vals2, patt3D[0, 0, :], color = "blue", alpha=0.5)
ax1.vlines([h/2 - w, h/2 + w], ymin=min(flux_vals)*1.2, ymax=max(flux_vals)*1.2, color = "green")

plt.show()

Where qq is my final (unfortunately continuous) function covering the mesh. Is there a way to have the data applied to the mesh without it being continuous, but rather remaining accurate as possible (to the mesh’s level of resolution) in the weak form?

Thanks in advance :slight_smile:

2

Use a DG-0 function to represent the data (piece wise constant per cell)

1 Like