How to represent piecewise functions

General
1

Hello everyone, I was calculating a 1D problem and applying different objective functions in different domains, I used the sympy.Piecewise function, but the oscillation value appeared in the calculation result:

from fenics import *
import sympy as sym
x = sym.symbols('x[0]')

mesh = Mesh()
with XDMFFile("mesh/mesh.xdmf") as infile:
    infile.read(mesh)
    
Q = FunctionSpace(mesh, 'P', 1)

M = sym.Piecewise((100, x <= 2e-6), (0, True))
M = sym.printing.ccode(M)
M = Expression(M, degree=2)

t_total = 200   
num_steps = 400   
dt = t_total / num_steps  

u = TrialFunction(Q)
v = TestFunction(Q)

ini_u = Expression('0', degree=1)
u_n = interpolate(ini_u, Q)

F = ((u - u_n)/dt)*v*dx - ((2-u/1000)*M)*v*dx
a, L = lhs(F), rhs(F)

u = Function(Q)
t = 0
for n in range(num_steps):
    
    t += dt

    solve(a == L, u, [])
        
    u_n.assign(u)
    
    set_log_level(LogLevel.PROGRESS)
    progress += 1

It can be seen that there are two “anomalous points” on both sides of the segmented function, a maximum value and a minimum value, which is obviously incorrect, I think it may be caused by the segmented function is not smooth, how should I solve this problem?
image

2

See Gibbs phenomenon - Wikipedia.

Standard finite element funtion spaces cannot exactly represent intra-cell discontinuities. Consider aligning your mesh with the material coefficient discontinuity. Or interpolate your material coefficient into a DG space which aligns with the mesh.

3

Thank you, I know the problem, if I still want to use CG cell type, how do I write the code?
Or it’s acceptable to have a smooth transition

4

Ok, the problem has been solved, thank you!