piecewise-definition of coefficients on subdomains

[gdmcbain]

Consider a Laplace equation with variable coefficient, ∇⋅ (k ∇ u) = 0, with k piecewise-constant on ‘regions’ or ‘subdomains’; e.g. steady conduction of heat between two blocks of different thermal conductivity (with no thermal contact resistance).

I think that such problems could also be done with mortar methods (which I see from ex04 are implemented, at least for two-dimensional problems, via InterfaceMesh1D); however, monolithic approaches can be convenient too. The issue here is to assign different values to k depending on which submesh an element belongs to, e.g. as defined in Gmsh as Physical Surface (in two dimensions). Should k be defined with ElementTriP0? And then populated using the cell_data returned by pygmsh.generate_mesh or the .cell_data attribute of the meshio.Mesh returned by meshio.read—I don't think skfem.mesh.from_meshio currently does this, it only looks for tags on facets, e.g. for mixed boundary conditions.

[kinnala]

The additional parameter w given to asm, i.e.

asm(bilinf, basis, w=w)

has shape number-of-elements x number-of-quadrature-points.

Each row corresponds to a single element. If the material field is constant within each element, you can simply create 1D array of length mesh.t.shape[1] (number-of-elements) and replicate this by the number of quadrature points. The result can be passed to asm.

[kinnala]

You are correct, only boundary tags are read:

scikit-fem/skfem/mesh/mesh.py

Line 213 in 0195b4e

def from_meshio(cls: Type[MeshType], meshdata) -> MeshType:

Separate Submeshes should be created for each subdomain.

Currently there is no nice way to use them.

[kinnala]

Ok, now that I meditated this for a while: I think we should create a method to GlobalBasis (say GlobalBasis.indicator) which takes in Submesh and outputs 2-d array which can be passed to asm via the keyword argument w.

This 2-d array represents the indicator function, i.e. it is 1 if the respective quadrature point belongs to the Submesh and 0 otherwise.

These indicator function arrays can be multiplied with specified values and summed together as necessary.

What do you think?

[gdmcbain]

That would be very useful for switching on a term in a subdomain, yes. Sticking with the example of the heat equation for concreteness, it's very common for tbe source term to be confined to a region. The Joule-heated wire in an insulator is a textbook case.

For the piecewise-constant conductivity though, this style would involve one asm for each region and there could be a few of those, so the P0 approach, simplified as above ('simply create') would be more convenient.

I think both idioms should be available and both work well side by side. Both are available in the FreeFem++ language, for example.

[kinnala]

Here is a method for InteriorBasis:

def indicator(self, submesh):
    ind = np.zeros((self.nelems, len(self.W)))
    ind[submesh.t, :] = 1.0
    return ind

A similar one could be created for FacetBasis to indicate facet subsets.

Then you could create Submeshes and input them to InteriorBasis.indicator:

left_ind = basis.indicator(m.submesh(lambda x, y: x<0.0))
right_ind = basis.indicator(m.submesh(lambda x, y: x>0.0))

During assembly you could define piecewise constant fields like this:

asm(bilinf, basis, w = 3.0*left_ind +2.0*right_ind)

[gdmcbain]

Sorry, I've been distracted away from this by two other aspects of the heat equation: trying to find a platform-independent implementation of a sparse Cholesky decomposition #86 (sksparse.cholmod.cholesky is easy in Linux, but it relies on SuiteSparse which I haven't be able to get installed under Microsoft Windows) and also a deadline for multicolumn right-hand sides #13.

Ideally the implementation of piecewise-definition of coefficients on subdomains will permit user-level code that is independent of the number of 'materials' (i.e. subdomains). If the mesh is generated in pygmsh or otherwise read in with meshio from a MSH2 .msh or MEDIT .mesh file, I think there will be a one-dimensional int array mapping elements to subdomains (cell_data in meshio).

This might require an extension of Mesh.from_meshio if Mesh.load is to be used (rather than direct recourse to pygmsh.generate_mesh or meshio.read). Would it be reasonable to simply attach (without any processing) new attributes to the returned Mesh corresponding to whatever the argument meshdata has in its point_data, cell_data, and field_data so that the user can access them? There might be lots of reasons to want to do this as a way to get data from up the pipeline into scikit-fem; e.g. loading in a finite element solution from a previous run or even another code. I might launch a new issue for this.

I'm envisaging the user supplying another mapping from subdomain to value of k. I think this could be turned into an elemental array for the coefficient by indexing. Say (for five elements divided into two subdomains):

region = np.array([0, 0, 0, 1, 1])  # elements -> subdomain
k = np.array([3., 2.])  # subdomain -> thermal conductivity

then indexing composes these two mappings; i.e.

 k[region]  # elements -> thermal conductivity

is

array([3., 3., 3., 2., 2.])

which can be expanded to the quadrature points for assembly as above.

[gdmcbain]

I've found a reference for the aforementioned insulated wire:

• Carslaw, H. S., & J. C. Jaeger (1959). Conduction of Heat in Solids (2nd ed.). Oxford University Press. §7.2.V, pp 191–192

This might be a better example than #90 . It demonstrates both features of subdomains:

• subdomain-dependent coefficient (the thermal conductivity)
• variational forms only active in a subdomain (the Joule heating)

In variational form, with a slight change of notation (so that h becomes a heat transfer coefficient), the problem is (grad u, k grad T) + h [u, T] = (u, A) where the square brackets denote the inner product over the boundary and k = k0 for 0 < r < a and k1 for a < r < b and A = A0 in 0 < r < a and 0 in a < r < b.

[gdmcbain]

I got that working in insulated.py on branch piecewise-89 #90.

I used a very quick hack in skfem.mesh.mesh.from_meshio

scikit-fem/skfem/mesh/mesh.py

Line 267 in a7ca493

mesh.cell_data = meshdata.cell_data

which will probably want some thought and then some tidying.

This is then used in the script:

scikit-fem/docs/examples/insulated.py

Line 60 in a7ca493

regions = mesh.cell_data['triangle']['gmsh:physical'] - 1

and for the subdomain-dependent coefficients as

scikit-fem/docs/examples/insulated.py

Lines 70 to 71 in a7ca493

	replicate = partial(np.tile, reps=(len(basis.W), 1))
	L = asm(conduction, basis, w=replicate(thermal_conductivity[regions]).T)

and the subdomain-dependent linear form as

scikit-fem/docs/examples/insulated.py

Lines 77 to 84 in a7ca493

	in_wire = (regions == 0).astype(float)
	@linear_form
	def generation(v, dv, w):
	return w.w * v
	f = joule_heating * asm(generation, basis, w=replicate(in_wire).T)

.

I haven't checked the result against the closed-form solution yet but it looks O. K.

[kinnala]

Yes, I can see why one would like to have all data outputted by meshio accessible.

I propose that we create a single attribute Mesh.external or Mesh.ext which contains all the field data, cell data, etc. of the respective meshio object (e.g. Mesh.ext.cell_data) and where the word external or ext emphasizes that we aim not to document its contents within scikit-fem since it might depend on the format.

Still, I will probably at some point improve the wrapper from_meshio to create Submesh-objects from subdomains since I find them to be convenient regarding my future plans. I'll probably move the gmsh-specific parts of from_meshio into a separate method at some point.

[gdmcbain]

Yes, I was starting to think of something like a .external attribute too. It could very well be that what's returned by meshio.read changes; the resolution of nschloe/meshio#175 doesn't seem like the last word on how to serialize subdomains.

[kinnala]

Check 37be054.

[gdmcbain]

Thanks for the review here. This is a very useful feature.

I'll probably continue to tweak the user-level code in the insulated-wire example #90 as I begin to make real use of the feature myself, but I think this issue can be closed for now.

Later I might also try comparing this monolithic approach with one based MeshInterface1d for the same boundary value problem.

[kinnala]

Sure, I've been actually planning to write a short article on the adaptive finite element solution of problems where we have multiple subdomains with different diffusion coefficients. It might be interesting to compare the standard adaptive approach and adaptive Nitsche coupling: Nitsche might end up in better accuracy with less DOF's in some specific cases.

My PhD thesis (https://aaltodoc.aalto.fi/handle/123456789/31486) was heavily based on the use of MeshInterface1d but I'm still not happy with the interface, especially when having more than two subdomains. I've been planning to address this problem at some point.