Finding dof locations

[kinnala]

One should be able to easily find global DOF locations using facilities such as InteriorBasis.get_dofs.

• Make local DOF locations to each element
• Map these to global coordinates
• Make nice interface for fetching them

[gdmcbain]

In #113, I noticed another case where this hurts. Besides coordinate-dependent Dirichlet data #112, in one of the Laplace examples, the exact solution was constructed using Mesh.p which only gave nodal values.

scikit-fem/docs/examples/ex13.py

Line 63 in 5f1687d

u_exact = 2 * np.arctan2(mesh.p[1, :], mesh.p[0, :]) / np.pi

This is gotten around in #113 by L²-projection, which might be a better idea anyway. That's easy as we do have access to the coordinates of the quadrature points when writing a linear form.

[kinnala]

I guess one could also implement Dirichlet conditions via an explicit L2 projection. This might introduce consistency error but I'm not sure.

[gdmcbain]

I hadn't thought of that Dirichlet via L2. Should work? Why do you think that it might introduce consistency error? I suppose I would have usually thought that L²-projection was "truer" (in some functional-analytic sense which I probably can't precise) that mere nodal interpolation. Is it that to avoid increasing the consistency error a different kind of projection is called for, suited to the correct trace space?

Would this end up giving identical answers (assuming all linear algebraic systems solved exactly) to the Lagrange-multiplier method, except that it's done in two stages (boundary then interior)? My understanding is that that effectively finds the flux such that the boundary values weakly satisfy the Dirichlet condition, i.e. w.r.t. the L² inner product, and I think that in a finite element space that set of boundary values is unique.

One to add to the sequence of examples proposed in #112 anyway, perhaps.

[gdmcbain]

It looks like the coordinates of the DOFs in ElementTriP2 are given by those of the nodes of the mesh followed by the midpoints of the facets:

np.hstack([mesh.p, mesh.p[:, mesh.facets].mean(axis=1)])

This matches, for example, what's obtained by L² projection using the w.x available inside a linear_form.

[kinnala]

I wonder if this should be closed now that we found a way to do it through L^2 projection and create an interface for performing the projections instead?

[kinnala]

Closing this since it's less relevant after adding utils.L2_projection.

[gdmcbain]

So if dof locations were to be made available through the API, where would that be? I figure through the InteriorBasis (haven't thought about FacetBasis), since they depend on both the Mesh and the Element. What about a property ib.x? Raising NotImplementedError by default, returning self.mesh.p for P1 and Q1 and as above for P2 (and maybe Q2, need to check)? L2 projection is good, but I feel that it loses some of the simplicity of the early examples. Plus serialization #158.

[kinnala]

I can try to sketch how this would work but I first need to understand what's actually useful. User probably wants to have global DOF locations corresponding to some subset of facets (e.g. some entry in Mesh.boundaries) and then the corresponding indexing to fill that information to the solution vector?

[gdmcbain]

User probably wants to have global DOF locations corresponding to some subset of facets (e.g. some entry in Mesh.boundaries)

Yes, that would be the main application.

An occasional one would be comparison with exact solutions as in

scikit-fem/docs/examples/ex13.py

Lines 42 to 44 in 0eb2dce

	u_exact = L2_projection(lambda x, y: 2 * np.arctan2(y, x) / np.pi, basis)
	u_error = u - u_exact
	error_L2 = np.sqrt(u_error @ M @ u_error)

Then that first line could be replaced with

u_exact = 2 * np.arctan(mesh.p[1], mesh.p[0]) / np.pi

[kinnala]

Status update. Now in this commit I added ElementTriP2.doflocs attribute: 778bb42.

Then you can do stuff like this:

In [1]: from skfem import *                                                                   
In [2]: m = MeshTri()                                                                                                                   
In [3]: e = ElementTriP2()                                                                    
In [4]: ib = InteriorBasis(m, e) 
In [5]: ib.mapping.F(e.doflocs.T)                                                            
Out[5]: 
array([[[0. , 1. , 0. , 0.5, 0.5, 0. ],
        [1. , 0. , 1. , 0.5, 0.5, 1. ]],

       [[0. , 0. , 1. , 0. , 0.5, 0.5],
        [0. , 1. , 1. , 0.5, 1. , 0.5]]])
In [6]: ib.element_dofs.T                                                                    
Out[6]: 
array([[0, 1, 2, 4, 6, 5],
       [1, 2, 3, 6, 8, 7]])

Now there is a correspondence between the locations in Out[5] and global DOF indices in Out[6].

In order to use this easily, we need to implement a (preferably vectorized) mapping from Out[6] to Out[5]. This can be done (conceptually) in a similar way as in MeshTri._build_mappings we build f2t.

[gdmcbain]

Oh, that's interesting. That's completely different to what I was going to try, which is just what was sketched above on 2019-02-12.

Note that mappings like that required from Out[6] to Out[5] can take a while, as profiled for #214.

Is

np.hstack([mesh.p, mesh.p[:, mesh.facets].mean(axis=1)])

cheating? (For ElementTriP2, and analogous for others.) Are there circumstances under which it could break? I envisaged adding a test that

This matches, for example, what's obtained by L² projection using the w.x available inside a linear_form.

holds.

[kinnala]

Yes that'd work in this case but I'm afraid those expressions can become tedious for complex elements. It also requires relying on some particular ordering of e.g. mesh.facets when going for #23 (for example ElementTriP3) which I'm planning to do soon.

I think we can reuse the information in GlobalBasis._build_dofnum to do the matching. There we are also creating the global DOF indexing for arbitrary elements.

[gdmcbain]

Ah, yes, that's obvious now, thanks.