known, but variable, coefficients

[gatling-nrl]

When working with variable coefficients, such as σ in

it is confusing how to specify the coefficient in the skfem syntax. Suppose I have a unit square mesh with a circular subdomain at the center (0.5,0.5) with radius of 0.1.

I could

...
basis = skfem.CellBasis(mesh, skfem.ElementTriP1())
...
def sigma_func(x):
    sigma_background = 1  # arbitrarily
    sigma_subdomain = 2  # arbitrarily
    x0 = (.5, .5)
    r = .1
    subdomain_region = np.sqrt((x[0]-x0[0])**2 + (x[1]-x0[1])**2) < r
    return np.where(subdomain_region , sigma_subdomain, sigma_background)
sigma_projected = skfem.utils.projection(sigma_func, basis)
sigma_interpolated = basis.interpolate(sigma_projected)  # can be used in skfem.asm

Or I could

basis_p0 = basis.with_element(skfem.ElementTriP0())
sigma_projected = basis_p0.zeros()
sigma_projected[:] = 1
sigma_projected[basis_p0.get_dofs(elements='subdomain')] = 2  # assuming my mesh has a region labelled 'subdomain'
sigma_interpolated = basis_p0.interpolate(sigma_projected)  # can be used in skfem.asm

(Or at least I thought I could. Above code doesn't work or report errors for me.)

But can I combine these? The subdomain label and skfem.utils.projection so as to be more generic and applicable to other than P1?

[kinnala]

It comes down to what your assumptions on sigma are and how much error you are willing to stomach. I think it doesn't make sense to approximate a discontinuous function using ElementTriP1 because it is continuous.

[kinnala]

It is also possible to project to piecewise constants. Then it looks like this:

[kinnala]

My suggestion is:

• if you have mesh conforming to the boundary of the subdomain, then you can do option 2.
• if you do not have a mesh conforming to the boundary of the domain, then do projection to ElementTriP0

Projection to ElementTriP1 doesn't make sense in most cases.

[gatling-nrl]

Not only that, but now I realize if you are trying to find sigma, e.g. by inverse methods, then using P1 isn't free. For the same mesh there are more unknowns in P1 than in P0. Of course, then if smoothness of sigma is used for regularization perhaps P1 is more desirable.

[gatling-nrl]

It is also possible to project to piecewise constants. Then it looks like this:

Can you show the code for that plot?

[kinnala]

I took your first snippet and replaced P1 by P0. Plot is done using plot(basis_p0, sigma_projected) from skfem.visuals.matplotlib.

[gdmcbain]

This problem is solved in ex17 'insulated wire'.