boundaries of subdomains

[gatling-nrl]

I propose extending the class Mesh to provide some ability to return the boundary of subdomains, even for those cases where the subdomain boundary is not used as a boundary condition.

scikit-fem/docs/examples/ex13.py

Lines 62 to 63 in 71eb03e

	for port, boundary in mesh.boundaries.items():
	fbasis = FacetBasis(mesh, elements, facets=boundary)

In Ex13 you use the boundary of the domain. But suppose the domain had an internal boundary between two subdomains. Would the following changes be a good idea?

Mesh.boundary_cells(self, subdomain=None, external=False)  # doesn't exist at present
Mesh.boundary_facets(self, subdomain=None)
Mesh.boundary_edges(self, subdomain=None)  # aren't facets and edges the same thing?
Mesh.boundary_nodes(self, subdomain=None)

Where external=False specifies you want the cells with edges on the boundary and non-boundary nodes inside the subdomain and external=True means cells with edges on the boundary and non-boundary nodes outside the subdomain.

Would these additions make it trivial to extend ex13 to measurements of flux across internal subdomain boundaries?

Is it worth extending the Mesh.*_satisfying() functions to have subdomain awareness?

If these are good ideas, I'll start a branch for them.

[gdmcbain]

This isn't exactly the same as 'internal boundaries in with_boundaries' #683 but some of that might help.

[kinnala]

Mesh.boundary_edges(self, subdomain=None)  # aren't facets and edges the same thing?

Nomenclature: In scikit-fem, edges exist only for three-dimensional meshes so that facets are something always shared between two elements of the mesh. In particular, we refer to the edges of triangular and quadrilateral meshes as facets.

[kinnala]

Mesh.boundary_facets(self, subdomain=None)

This idea is fine. However, implementing this may not be completely straightforward.

The reason why we are able to find those boundary faces in a reasonable time is that the second row of f2t contains -1 if the facet is on the boundary. This is done in _build_entities and done only once for the whole mesh.

I think the easiest way to hack this together would be to create another Mesh object consisting of only the subdomain elements, using Mesh._reix (or similar), and then call Mesh.boundary_facets() on the resulting object.

In particular, implementing boundary_edges(subdomain='blah') for MeshTet1 and MeshHex1 will be otherwise pretty painful. However, I think boundary_edges(subdomain='blah') may not be needed if the use case is to calculate the flux?

[kinnala]

Here is one idea for the implementation:

m = MeshTri().refined(4).with_subdomains({'mid': lambda x: x[0] <= 0.5})
v, c = np.unique(m.t2f[:, m.subdomains['mid']].flatten(), return_counts=True)
facets = v[np.argwhere(c == 1)].flatten()

I tried drawing facets:

In [29]: basis = Basis(m, ElementTriP1())                                     
In [30]: dofs = basis.get_dofs(facets)                                        
In [31]: u = basis.zeros()                                                    
In [32]: u[dofs] = 1                                                          
In [33]: m.draw(point_data={'u': u}).show() 

[kinnala]

Maybe use nonzero instead of argwhere.

[kinnala]

Mesh.boundary_cells(self, subdomain=None, external=False)  # doesn't exist at present

Will this include only elements having a facet on the boundary or also elements that have only one node/edge at the boundary?

The former can be done by a simple m.f2t[0, m.boundary_facets()]. The latter is somewhat nontrivial to do efficiently.

Edit: We should call it boundary_elements so that it is in line with elements_satisfying etc.

[gdmcbain]

There doesn't seem to be. I'm dealing with an "electrode" where the conductivity of the electrode might be >1e6 and of the surrounding medium 1. This makes the gradients inside the electrode very small (zero for an ideal conductor).

The domain integral then would be over the "surrounding medium". This is much more accurate than a boundary integral. See ex13 and the ANZIAM paper.

[gdmcbain]

My fall back plan is essentially that exactly. Since generalized ohms law is governing equation, for any fixed conductivity the current will be a linear function of voltage. So the worst case is to solve the problem 2 times to get that linear function, and then a 3rd time to get the specified current. I have the luxury of relatively simple domain shapes and meshes, but I'm running many thousands of variations so adding any kind of non-linear search to the solution is probably to expensive.

Not necessarily. This kind of technique (‘Schur complement’) can be used on a large scale; e.g. one approach to incompressibility in fluid mechanics is to construct a function that for a given pressure field computes the divergence field, then one ‘solves’ that function to find the pressure that gives zero divergence everywhere. That function is affine (linear plus constant) so can be converted to linear and solved with conjugate gradient, MINRES, &c.. This is what's done in ex30.

scikit-fem/docs/examples/ex30.py

Lines 96 to 109 in be9642f

	def dilatation(pressure: np.ndarray) -> np.ndarray:
	"""compute the dilatation corresponding to a guessed pressure"""
	return -B @ flow(pressure)
	dilatation0 = dilatation(basis['p'].zeros())
	K = LinearOperator((basis['p'].N,) * 2,
	lambda p: dilatation(p) - dilatation0,
	dtype=dilatation0.dtype)
	pressure = solve(K, -dilatation0,
	solver=solver_iter_krylov(minres),
	M=build_pc_diag(asm(mass, basis['p'])))

[kinnala]

I'm trying to get a handle on how skfem.FacetBasis works today.

Facet integration is implemented differently than in almost any FEM codes I have read. I learned the technique from a colleague of mine several years ago (credits in the README). Consider a triangle mesh. The idea is to first take the reference element of the boundary cell (line) and create quadrature points. Then we map those points to the facets of the triangle mesh using MappingAffine.G and finally back to the reference triangle using MappingAffine.invF. Now we have a set of quadrature points on the correct facet of the reference triangle and can evaluate the local basis.

The benefit is that we automatically get gradient information which is normally a challenge (based on many questions in scicomp StackExchange). I guess the downside is that by mapping back and forth you do a bit of extra work.

However, you can think of it as follows: by changing Basis to FacetBasis you are changing only the integration domain to the boundary. Most of the forms are compatible with both with the exception that FacetBasis includes w.n for outward normal.

[gatling-nrl]

Like that?

[kinnala]

Yes.