Support MatIS and BDDC

[pbrubeck]

Description

Support assemble(a, mat_type="is")

[connorjward]

This encodes the boundary conditions into the definition of the matrix, which is a bit bizarre. Why can't MATIS set lgmaps after the fact?

[pbrubeck]

This does not encode BCs, just the removal of repeated ghost dofs. This is the correct way of using MatIS

[pbrubeck]

The way we handle BCs on matrices requires changing the LGmap after creating the matrix. This is not possible for a MatIS, as the LGMap is what defines it. If you change the LGMap then the matrix is destroyed.

For MatIS, BCs are never imposed on the LGMap. Instead we assemble the matrix without BCs, and we zero out the BC rows and columns.

[wence-]

I think a bunch of the indexing is wrong when constructing the masked lgmaps?

[wence-]

AKA getHeightStratum(0)

[wence-]

If block is True and lgmap has block_size > 1 then this produces out of bounds indices.

[wence-]

Or conversely, if this is a mixed space where one component has block size > 1, then this masks out dofs even when running on a single process (which seems wrong).

[wence-]

I think you're trying to mask out those entries in lgmap.indices which belong to dofs that are not reachable from owned cells?

I think for a non-mixed space, this is (if block=True):

cmap = V.cell_node_map()
mask = np.setdiff1d(
    cmap.values_with_halo[cmap.iterset.size:], 
    cmap.values[:cmap.iterset.size],
)

For a mixed space, you need to index through to local_ises, I think like this?

mask_pieces = []
for i, W in enumerate(V):
    iset = V.dof_dset.local_ises[i]
    cmap = W.cell_node_map()
    to_mask = np.setdiff1d(cmap.values_with_halo[cmap.iterset.size:], cmap.values[:cmap.iterset.size])
    bs = iset.block_size
    if bs > 1:
        to_mask = np.concatenate([i + bs * to_mask for i in range(bs)])
    mask_pieces.append(iset.indices[to_mask])
mask = np.concatenate(mask_pieces)

[pbrubeck]

cell_node_map is not defined for extruded meshes, but also there is no DM for extruded. The correct thing is to execute a pyop2 kernel that visits only the dofs that are reachable from owned cells.

[pbrubeck]

Only VectorElements or TensorElements can have blocksize != 1. Any MixedElement has blocksize = 1

[pbrubeck]

@ksagiyam is it possible to get a RestrictedFunctionSpace with the unghosted FunctionSpace (the space that corresponds to DistributedMeshOverlapType.None)?

[pbrubeck]

Alternatively, if I reconstruct the FunctionSpace with the mesh without the overlap, would I get consistent dof orderings? How do I recover the mesh without overlap?

[ksagiyam]

You can make a RestrictedFunctionSpace with boundary_set=ghost entities.

DoF ordering on each entity is the same only for CG and DG spaces (assuming that cone orderings are preserved). Otherwise DoFs depend on the global vertex numbers, which are arbitrary.

I think what we want is to be able to pass ignore_halos=True to op2.Sparsity and have op2.DataSet.lgmap (and similar) return the correct lgmap in the first place depending on the value of ignore_halos?

[stefanozampini]

@pbrubeck Is this still broken? We can do mixed Poisson H(div) x L2 if we have monolithic l2gmaps