Closed gdmcbain closed 5 years ago
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
.
You are correct, only boundary tags are read:
Separate Submeshes should be created for each subdomain.
Currently there is no nice way to use them.
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?
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.
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)
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.
I've found a reference for the aforementioned insulated wire:
This might be a better example than #90 . It demonstrates both features of subdomains:
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.
I got that working in insulated.py on branch piecewise-89 #90.
I used a very quick hack in skfem.mesh.mesh.from_meshio
https://github.com/kinnala/scikit-fem/blob/a7ca4936db69ed7ce62698110c5b599a64b88306/skfem/mesh/mesh.py#L267
which will probably want some thought and then some tidying.
This is then used in the script: https://github.com/kinnala/scikit-fem/blob/a7ca4936db69ed7ce62698110c5b599a64b88306/docs/examples/insulated.py#L60 and for the subdomain-dependent coefficients as https://github.com/kinnala/scikit-fem/blob/a7ca4936db69ed7ce62698110c5b599a64b88306/docs/examples/insulated.py#L70-L71 and the subdomain-dependent linear form as https://github.com/kinnala/scikit-fem/blob/a7ca4936db69ed7ce62698110c5b599a64b88306/docs/examples/insulated.py#L77-L84 .
I haven't checked the result against the closed-form solution yet but it looks O. K.
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.
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.
Check 37be054.
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.
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.
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 asPhysical Surface
(in two dimensions). Should k be defined withElementTriP0
? And then populated using thecell_data
returned bypygmsh.generate_mesh
or the.cell_data
attribute of themeshio.Mesh
returned bymeshio.read
—I don't thinkskfem.mesh.from_meshio
currently does this, it only looks for tags on facets, e.g. for mixed boundary conditions.