Finite element calculations on subdomains

[jamiebramwell]

@samuelpmishLLNL and I are interested in doing finite element calculations on subdomains (either volumetric or surface elements) of the mesh. This could be useful for:

• Different physics in different domains
• Boundary integral contributions
• Non-virtual function material evaluations

I was pointed to #2005 where these ideas are explored. It seems like there are many ways to approach this including:

• Defining finite element spaces directly on a subdomain (approach taken in FESpaces on Subdomains #2005)
• Extracting a sub-mesh from an existing mesh (an extension of something like skin_mesh )
  • It would be best if this new mesh maintained connection maps to the parent mesh.
• A generated "mesh view"-type class that defines iterators over a subset of the mesh.

Any thoughts on these ideas? I think I am partial to the sub-mesh concept since it avoids touching all of the current mesh logic. However, creating the connections back to the parent mesh could be messy. Unfortunately something like that is required if a user wants to pursue monolithic-type solves. I look forward to hearing everyone's thoughts!

[tzanio]

ping: @mlstowell and @jakubcerveny

[sshiraiwa]

I would like to share the path we tool in Petra-M (https://github.com/piScope/PetraM_Base/tree/master/petram, https://piscope.psfc.mit.edu/index.php/Petra-M), which is essentially GUI interface to use MFEM in this context, and what I noticed. We took the sub mesh approach as you described. To do this, we assign different attribute and bdr_attribute numbers to all domains and boundaries. This part is done when generating a mother mesh, so that when MFEM load a mesh file, it is all set in such a way. Then, we diagnose the connectivity of domains and boundaries to find low order element numbering. For example, we look for edges which is shared between surface 1 and 2, collect all of them, and assign a unique number to each of such a group. We did the same thing for vertices. We call it mesh extension (https://github.com/piScope/PetraM_Base/blob/master/petram/mesh/mesh_extension.py). This numbering allows a user to choose the sub-domains to in our GUI. A typical work-flow is that a user can create a 3D geometry and mesh, then, choose where on the mesh you want to define your equation. Then we generate a partial mesh for the geometry a user selected to build PDE system. This was done in here ( https://github.com/piScope/PetraM_Base/blob/master/petram/mesh/partial_mesh.py), and so far it works fine for conforming mesh. Not sure how much extra work I need to use non-conforming mesh. We also wanted to couple simulaitons in different domains. So, we generated a degree-of-freedom mapper. (https://github.com/piScope/PetraM_Base/blob/master/petram/helper/dof_map.py) This routine compares DoFs defined on two meshes, and find which DoFs in on FES is the same as a DoFs on the other FES. Of course,. FES is supposed to be the same type and order. The vector elements (ND and RT) sometimes has multiple DoFs at the same place, and we needed to solve the small linear system. All written in Python using PyMFEM. So it maybe rather straight forward to write the same logic in C++, although I used a convenient Scipy and Numpy routines such as KD-Tree. There are two things i didn’t do 1) When generating the partial mesh (or better call it sub-mesh?), we did not record which mesh element is taken from the mother mesh. This book-keeping would be useful to generate the mapping of DoFs between meshes. However. we decided to use ParMesh without doing anything special, the partial mesh is partitioned differently compared to the mother mesh, so this book keeping might have tuned out to be rather complicated. 2) We haven’t support an adaptive mesh refinement so far. I naively think that we need to do the refinement on the mother mesh level, and propagate it to sub-meshes. Hope this helps to figure out the path forward. Best,

[v-dobrev]

Here is another related discussion: #1729.

[termi-official]

Hi,

first I want to throw in (geometric) Schwarz methods as another potential application (although we can already implement them via sub-meshes).

Do you really think we can dodge touching or duplicating the existing mesh logic with sub-meshes? As soon as we have subdomains and adaptivity (and therefore load-balancing) things seem to be a bit tricky with the sub-mesh extraction and transferring the fields, as we artificially additional bookkeeping about the elements involved in the load balancing. If we ignore adaptivity for one second, then I agree that the sub-mesh may be simpler than other approaches. But with adaptivity involved I guess the mesh logic has to be touched or duplicated with any approach taken.
As I just saw this discussion right now, please note that I have some more details post here #2005 (comment) .

@jakubcerveny should we shift the discussion from #2005 to this thread to keep things together and stall the PR for now?

Best regards,
Dennis

[jakubcerveny]

Hi Dennis, OK lets move the discussion here.

One note about adaptivity: already now, with NCMesh, the Mesh object serves as a kind of proxy that is regenerated every time the mesh is refined or rebalanced. I'm not saying this is optimal (a lot of time is spent generating many of the tables in the Mesh class, this could probably be improved by filling them directly by NCMesh), but there is a precedent. As any real AMR with subdomains will somehow need to work with the NCMesh object, I'd say a simple approach is to extract multiple Mesh representations from the NCMesh at every update.

I also think the (Par)FiniteElementSpace should not know anything about subdomains, it should continue to work the same over a Mesh or some compatible object.

[termi-official]

I think we are on the same page Jakub. :) The (Par)FiniteElementSpace should not know anything explicit about subdomains, hence the code I started in the branch is probably a really bad idea, as it violates the law of Demeter.

My new idea was not about extending, but refactoring the (Par)FESpace with some new lightweight class which handles the communication between the (Par)Mesh and (Par)FESpace to avoid the copies. I.e. the (Par)FESpace will at no point directly call or see the (Par)Mesh. Still, I am not 100% confident if this will turn out to be faster in the end. I sadly could not find much time to think this through the end yet. My hope here is that we can economize on memory allocations.

[v-dobrev]

In the interest of moving the discussion forward: is anyone interested in writing down a high-level design proposal for this? For example, a list of new classes and methods with short descriptions of intended behavior.

[termi-official]

I can put together a proposal for the mesh view based approach, based on what I learned from exploring subdomains in #2005 . However, going down this road we will probably end up again with the issue why I stopped development on this branch for now: We have to somehow interface with the communication infrastructure in ParMesh. I have not found a reasonable good way to accomplish this, yet.
Anyway, if there is interest I can put a proposal together.

[termi-official]

I just gave it a shot and sketched a first software model.

Try 1:

UML Source

Note that in addition to this we may have to touch the IO logic to preserve subdomain information.

Edit 1: Probably I should elaborate a bit. I basically introduced a (Par)MeshView as a drop-in replacement for the (Par)Mesh inside the (Par)FiniteElementSpace. The (Par)FullMeshView basically acts as a direct proxy to the (Par)Mesh and should not come with any overhead to the current implementation. Obviously subdomain will not be free in terms of performance, as we need to pay for the overhead of tracking subdomain information. To accomplish this two additional (Par)MeshViews are introduced. One to handle "normal" subdomains in a strict sense (like having some solid surrounded by a fluid) and one fancy where the subdomain has codimension > 0, which may come handy for problems involving traces or other surfaces.

For the serial case I can easily refactor the existing branch #2005 to make it work. However, in parallel things get a little bit more tricky, since we have to handle cases where non-local information (through the ParMesh's MPI graph) is required. Think about the prolongation/restriction matrix construction.