Piecewise function in discontinuous Galerkin space


51
views
0
24 days ago by

Hello everyone. I want to define a piecewise function in a 1st-order discontinuous Galerkin space DG1. Consider the function,
\begin{equation} f(x,y) = \begin{cases} 0 & x < 1 \\ 1 & x > 1 \end{cases} \end{equation}
Determining the nodal values in the DG1 space at the jump (i.e., `x=1`) is where I run into problems. For the elements to the left of `x=1`, I want all the nodes to have the value zero; and to the right of `x=1`, I want the element nodes to all have the value one. How do I assign the correct nodal values?

One workaround I have found is to first define an element-wise constant discontinuous space, DG0, and then to interpolate values from DG0 to DG1, as I've shown in the example code block below. Although I would like to avoid creating an extra function space if possible.

Also, is there a more concise method to assign values to a function from a numpy array than using f.vector().set_local()?

from dolfin import *
import numpy as np

# create a mesh
mesh = RectangleMesh(Point((0, 0)), Point((2, 1)), 2, 1, "left")

# function spaces
DG0 = FunctionSpace(mesh, "DG", 0)
DG1 = FunctionSpace(mesh, "DG", 1)

# coordinates for each function space
x0, y0 = DG0.tabulate_dof_coordinates().reshape((-1, 2)).T
x1, y1 = DG1.tabulate_dof_coordinates().reshape((-1, 2)).T

# array for DG0 space
a = np.zeros(len(x0))
a[x0 > 1] = 1

# functions in DG0 and DG1 spaces
f0 = Function(DG0)
f0.vector().set_local(a)
f1 = interpolate(f0, DG1)

# print
print(a)
print(f0.vector().get_local())
print(f1.vector().get_local())

Community: FEniCS Project

1 Answer


1
19 days ago by
You can avoid the DG0 space if you are willing to loop over cells

dm = V.dofmap()
arr = f.vector().get_local()

for cell in cells(mesh):
    val = 0 if cell.midpoint().x() < 1 else 1
    for dof in dm.cell_dofs(cell.index()):
        arr[dof] = val

f.vector().set_local(arr)
f.vector().apply('insert')​

Using cell.midpoint()[0] may be more future proof, but the above code probably works fine for now (not tested)
Please login to add an answer/comment or follow this question.

Similar posts:
Search »