basis.get_elements()

[gatling-nrl]

Now that basis.get_dofs() is coming along.... can we add basis.get_elements()

The use case is

scikit-fem/docs/examples/ex17.py

Lines 56 to 58 in 08a695d

	conductivity = basis.zero_w()
	for subdomain, elements in mesh.subdomains.items():
	conductivity[elements] = thermal_conductivity[subdomain]

But what if I don't have labeled subdomains?

I propose that we could

sigma = basis.zeros_like('elements')
copper_area = basis.get_elements(lambda x: x[0] < 0.5)
sigma[copper_area] = 5.96e7 # S/m

with a function very similar to get_dofs.

This sigma is easily plotted with matplotlib.pyplot.tripcolor(..., facecolors=sigma), but can it be passed to the Forms? Is it easy to use as the coefficient for Laplace equation?

[kinnala]

I feel zero_w is not really needed anymore. Same applies to the proposed zeros_like.

You can initialize another Basis with ElementTriP0 (piecewise constant function) and then zeros will give you the correct shape. This is how we should think about piecewise constant functions: FE functions defined by piecewise constant basis. There is Mesh.elements_satisfying so the second step also exists.

If you want to use sigma inside a form, then you make sure both Basis objects have the same number of integration points so they are compatible. Easiest way is to initialize the P0 basis like basis0=basis.with_element(ElementTriP0()). Then you can call basis0.interpolate(sigma) and pass the resulting object to assemble as a kwarg.

[kinnala]

After I get DofsView fixed the following is also an option:

basis0 = basis.with_element(ElementQuad0())
conductivity = basis0.zeros()
conductivity[basis0.get_dofs(elements='core')] = 101.
conductivity[basis0.get_dofs(elements='annulus')] = 11.

[gatling-nrl]

Is this how it should all come together? (with the flatten() being dropped soon)

    # Laplace's equation with variable coef
    @skfem.BilinearForm
    def laplace(u, v, w):
        return dot(w['sigma'] * grad(u), grad(v))
    mesh = skfem.io.from_meshio('mesh.msh')
    basis1 = skfem.Basis(mesh, skfem.ElementTriP1()) # piecewise linear
    b_dofs = basis1.get_dofs() # all boundary dofs by default
    b_values = basis1.zeros() # will have shape of full solution, but gets masked by b_dofs later
    b_values[basis1.get_dofs(elements='left').flatten()] = 1
    basis0 = basis1.with_element(skfem.ElementTriP0()) # piecewise constant
    sigma = basis0.zeros()
    sigma[basis0.get_dofs(elements=lambda x: x[0]<0.5).flatten()] = 1
    sigma[basis0.get_dofs(elements=lambda x: x[0]>=0.5).flatten()] = 2
    A = laplace.assemble(basis1, sigma=basis0.interpolate(sigma))
    u = skfem.solve(*skfem.condense(A, D=b_dofs, x=b_values))

[gatling-nrl]

In this line

L = asm(conduction, basis, conductivity=basis0.interpolate(conductivity))

is the basis0 supposed to be just basis? I'm not sure I understand why conductivity needs to be interpolated on the basis that created it. Or maybe I'm not understanding what interpolate returns.

[kinnala]

My example is correct. interpolate will take an array of DOF values (basis function coefficients) and interpolate the discrete function at the quadrature points. conductivity needs to be interpolated using P0 (piecewise constant) FEM basis functions in basis0. There simply is no way basis can figure out in what finite element basis conductivity is defined. interpolate will return DiscreteField which is essentially conductivity evaluated at quadrature points. After a field is interpolated at the quadrature points none of this matters, all discrete functions can be added or multiplied together if the quadrature points match.

[kinnala]

There are generally two ways one might represent functions in FEM. One is values at nodes (conductivity) and one is values at quadrature points (basis0.interpolate(conductivity)). You can go back and forth using L2 projection and interpolation.