Restriction to a subdomain : dofs

[tsokar]

Hi,
first of all, thank you for all the work you put on this project !

I was digging in example 26 where the resolution of the problem is restricted to a subdomain through condensation of the outer part. I wanted to do something similar for an application, but defining the problem on the core only (dofs in particular).
In the example, the elements where the basis is evaluated are restricted to the core subdomain.

core_basis = Basis(mesh, basis.elem, elements=mesh.subdomains['core'])

This is fine for the assembly, but the dofs of the basis are still defined on the whole mesh : basis.dofs.N == core_basis.dofs.N.

Is there a mean to define the basis, but with both elements and dofs only in the subdomain ?

Thanks !
T.

[kinnala]

There is no such feature right now.

I have started sketching such a feature multiple times. But so far I have not been happy with the result.

Perhaps this request will finally give me enough motivation to get it actually done. My personal issue has been related to FacetBasis where I'd like to restrict to DOFs on boundary. However, that seems to be a more challenging problem than restricting DOFs on subdomain so perhaps I will solve this first.

My workaround has been to manually slice the resulting matrices.

[kinnala]

Here is a diff which takes you up to assembly:

modified   skfem/assembly/dofs.py
@@ -22,6 +22,7 @@ class DofsView:
     facet_rows: Union[ndarray, slice] = slice(None)
     edge_rows: Union[ndarray, slice] = slice(None)
     interior_rows: Union[ndarray, slice] = slice(None)
+    _cache = {}
 
     def __repr__(self):
         nnodal = len(np.unique(
@@ -194,6 +195,20 @@ class DofsView:
                              ix=self.interior_ix,
                              rows=self.interior_rows)
 
+    @property
+    def element_dofs(self):
+        if 'element_dofs' not in self._cache:
+            element_dofs = self.obj.element_dofs[:, self.interior_ix]
+            element_dofs = (np.unique(element_dofs.flatten(),
+                                      return_index=True)[1]
+                            .reshape(element_dofs.shape))
+            self._cache['element_dofs'] = element_dofs
+        return self._cache['element_dofs']
+
+    @property
+    def N(self):
+        return np.max(self.element_dofs) + 1
+
     def __getattr__(self, attr):
         return getattr(self.obj, attr)

Now you can assemble as follows:

from skfem import *
from skfem.models import laplace

m = MeshTri().refined(6).with_subdomains({
   'left': lambda x: x[0] < 0.5,
})

e = ElementTriP1()
dofs = Dofs(m, e)

In [8]: basis = Basis(m, e, dofs=dofs.get_element_dofs([0, 1]), elements=[0, 1])
Unable to calculate global DOF locations.

In [9]: laplace.assemble(basis)
Out[9]: 
<6x6 sparse matrix of type '<class 'numpy.float64'>'
   with 14 stored elements in Compressed Sparse Row format>

But now the issue is that what to do with the solution vector? Its length is only 6 so there are no values for the remaining, unsolved nodes in the mesh. All the plotting routines and export routines which expect all nodes to have something will fail unless they are modified somehow.

[tsokar]

Thanks a lot, I will try this !
Concerning the solution vector, I don't mind setting the remaining values arbitrarily to zero in order to use the export routines as is.

[kinnala]

Setting elements=[0, 1] in Basis constructor effectively disables integration over all other elements than 0 and 1. What if you then add ones in the diagonal of the disabled DOFs? That would also lead to the same result. I'll try it out.