Evaluate the jump discontinuity along the interface of two elements #4859
Replies: 3 comments 8 replies
-
|
I'm not sure that the term "evaluate" is well-defined here. |
Beta Was this translation helpful? Give feedback.
-
|
Isnt |
Beta Was this translation helpful? Give feedback.
-
|
I work together with @apalha and we made a simple example to check the definitions. We started with a simple function (on quads, since we are more familiar with those), using the following script. from firedrake import *
import numpy as np
np.set_printoptions(precision=2, suppress=True) # To reduce visual noise
# Create a simple unit square mesh with 2x1 elements, using quadrilaterals.
mesh = UnitSquareMesh(2, 1, quadrilateral=True)
# Create a function in the discontinuous space.
Wh = FunctionSpace(mesh, "DQ", 1)
vals = np.zeros(Wh.dim())
vals[0] = 1.0
vals[1] = 1.0
vals[6] = 1.0
vals[7] = 1.0
u_h = Function(Wh, val=vals, name="u_h")This creates an element-wise linear (in x) function which matches at the interface. Its gradient (in the x-direction) is constant on each element, where the value of the gradient on the left element is about 3.5 and on the left -3.5. This is also visible in the two figures (first the function, then the gradient).
If we now compute the 'jump' of the gradient like this # Compute (half) the jump of the gradient across facets.
n = FacetNormal(mesh)
Jump_gradu = dot(grad(u_h), n)
Vjump = FunctionSpace(mesh, "HDiv Trace", 0)
jump_vals = project(Jump_gradu, Vjump)
print(jump_vals.dat.data_ro[0])
print(jump_vals.dat.data_ro)the output is 3.4641016151377553
[ 3.46 0. -3.46 0. -3.46 0. 0. ]Note that the actual jump is about 7.0, while the average is 0. The other two non-zero entries are related to the boundary facets, which we ignore here (as @dham mentioned, these are probably best set to 0, though that would be too much for this example. In our original case, this is already the case.). So the numbering of facets seems to first number internal facets, then boundary facets. If we instead define the function on the same mesh as Wh = FunctionSpace(mesh, "DQ", 1)
vals = np.zeros(Wh.dim())
vals[0] = 1.0
vals[1] = 1.0
vals[6] = -1.0
vals[7] = -1.0
u_h = Function(Wh, val=vals, name="u_h")the function and its gradient look like If we then compute the 'jump' of the gradient like this again # Compute (half) the jump of the gradient across facets.
n = FacetNormal(mesh)
Jump_gradu = dot(grad(u_h), n)
Vjump = FunctionSpace(mesh, "HDiv Trace", 0)
jump_vals = project(Jump_gradu, Vjump)
print(jump_vals.dat.data_ro[0])
print(jump_vals.dat.data_ro)we instead get -3.3053706985194236e-17
[-0. 0. -3.46 0. 3.46 -0. 0. ]Note that the actual jump is 0, while the average is -3.5. The other two non-zero entries are again related to the boundary facets, one of which has now changed sign (as expected). In conclusion:
Thank you both for your help in figuring this out. |
Beta Was this translation helpful? Give feedback.




Uh oh!
There was an error while loading. Please reload this page.
-
Dear all,
We need to evaluate the jump discontinuity along the interface of two elements, i.e., evaluate the following expression at points
We are stuck at evaluating the last expression. We can create it, but it is only usable in an assembly, not evaluation.
Is there any way to evaluate an expression on the + and - side of an element?
Thank you for your help.
Kind regards,
-artur palha
Beta Was this translation helpful? Give feedback.
All reactions