drwells
commented
1 year ago overview
@boyceg and I are going to use this issue to develop a roadmap for this project. If we want to integrate this with IBAMR then the best path forward is to use the same mechanisms we use for fdl::IFEDMethod to associate FE information (DoFHandlers, vectors, etc) with SAMRAI patches. All that association code is already written and reduces this project's scope considerably: once you have an fdl::PatchMap then you can use plain serial data structures. This makes development and testing pretty easy.
roadmap
- find all intersections between each finite difference stencil and the mesh a. in IBAMR, this is done element-by-element, with some bookkeeping to look only over a small subset of the mesh elements b. note that it is critical that intersections not be double-counted. in an element-by-element approach, we need to make sure we don't double-count intersections that occur on element edges or nodes
- evaluate jump conditions at those intersections and evaluate the related correction terms, which appear as RHS terms in structured Cartesian grid equations a. lowest order jump conditions are just the normal and tangential components of surface forces b. IBAMR currently uses the element geometry to describe the surface normals, but it might be better to use a globally continuous description of normal vectors c. it is not clear to me that we need to project the jump conditions onto basis functions that live on the finite element surface mesh; this is what was in Amin's initial implementation, but the DG scheme that Amin and Michael are using now is closer to just evaluating the jump conditions pointwise from the surface force F and the surface normal/tangential directions
- tests a. we can do basic nontrivial tests just by applying forces on a stationary boundary and solving the N-S equations b. also we want to have unit tests for the basic intersection algorithms c. we can always start an IIM implementation using IB for velocity interpolation d. but there is no reason to do a fancy corrected velocity interpolation scheme until jump conditions are imposed in the N-S solver
Implementation
High-level view
Like IFEDMethod (which has an additional layer for
NodalInteraction/ElementalInteraction), we should write these classes with
three distinct layers:
- A class which inherits from
IBStrategywhich is solely responsible for book-keeping and hooking into the rest of IBAMR. - Data structures which efficiently store and grant access via iteration to intersections.
- Functions for computing intersections or using intersections to compute forces on the Cartesian grid.
Importantly (for the purposes of testing and clean code) a layer will only call
the layers under it: i.e., layer 3 will never call layer 1 (nor should it know
anything about fdl::IIMethod). This architecture is summarized in fiddle's readme.
Layer 1: fdl::IIMethod
Most of the book-keeping work is already done by fiddle: see https://github.com/drwells/fiddle/blob/master/include/fiddle/interaction/ifed_method_base.h. That class can set things up and determine which elements intersect which patches.
We should do this last since it will be relatively easy to bolt the pieces together.
Layer 2: data structures (#184)
One such data structure is PatchIntersectionMap, which builds on PatchMap to
store all the intersections between the FE mesh and a given Patch. The
iterators are also responsible doing the easy to get wrong index calculations
with SideIndex.
We might need more but fundamentally we need something to let us easily loop over intersections so that we can implement lots of IIM algorithms.
next steps
- [ ] Use the new intersection algorithms to populate
PatchIntersectionMap. To start, we should add a constructor similar toNodalPatchMaps. It might be easier to set this up with some auxiliary functions first before we write the class constructor.
Layer 3: intersection algorithms (#179 and more) and IIM algorithms
geometric primitives
We should implement basic intersection algorithms. For starters we should write code which intersects a ray with a triangle or line.
In principle, we might want to cover stencils which intersect cells an arbitrary number of times, but for now lets just stick with single intersections.
There's a couple of different ways to implement this:
- We could set up an
rtreeand then intersect each stencil with thertree: that would quickly tell us which elements it might intersect - We could loop over all elements and for each find which stencils intersect: this would require doing a second step to ensure that we do not double-count any stencils.
- We could provide a line/triangle to a function as well as some description of
the stencil and return
std::optional<double>(i.e., return nothing if no intersection exists and otherwise the intersection's convex combination coefficient)
I think some combination of 1 and 2 is the right way forward:
- get the physical-space bounding box for an element
- convert it into an index-space bounding box
- check just the stencils (defined as pairs of cell centers) which start or end in the index-space bounding box
If we use std::optional (available as std_cxx17::optional in deal.II) then that function can find either (like Haskell's Maybe: I don't know the jargon for this in C++ yet) Nothing or Just c in which c is the convex combination coefficient.
IIM algorithms
These should be similar to the functions in interaction/interaction_utilities.h.
normal vectors
There are a couple of different ways to do this: e.g., Phong shading. The best
way to implement all of these in a general way is with a finite element field.
For example, with a tet mesh we can do FESystem<2, 3>(FE_SimplexDGP<2, 3>(0), 3)
to implement standard normal vectors (3 values per element, cell-centered)
or FESystem<2, 3>(FE_SimplexDGP<2, 3>(p), 3) to do various generalizations of
Phong shading (with p = 1 being the standard linear interpolation). Possible
function signature:
template <int dim, int spacedim = dim>
LinearAlgebra::distributed::Vector<double>
compute_phong_normal_field(const DoFHandler<dim, spacedim> &dof_handler,
const Mapping<dim, spacedim> &mapping);Like everything else in fiddle we should write these functions first, debug them, and then hook them into a bigger solver.
The major advantage to using FE fields for normal vectors is that we can utilize the preexisting data transfer mechanisms in fiddle to get data to the IB data partitioning easily as computing Phong normal vectors requires additional ghost data only present in the normal FEM data partitioning.
next steps
- [ ] surface meshes in deal.II are well-oriented, so we are guaranteed that the normal vectors on each cell are all either pointing in or out. These vectors are available through
FEValues. Since we might want to use multiple representations of the normal vectors we should start by expressing the default piecewise-constant normal vectors as anFESystem<dim, spacedim>(FE_DGQ<dim, spacedim>(0), spacedim)orFESystem<dim, spacedim>(FE_SimplexDGP<dim, spacedim>(0), spacedim)FE field. - [ ] We can generalize the previous step to general discontinuous or continuous vector fields.
level sets
From a conversation with Qi last week: we need a way, when doing velocity interpolation on the interface, to distinguish between grid points outside or inside the interface. Our conversation convinced me that the only reasonable way to do this is with a level set method in which we tag cells outside the interface with 1 and those inside with -1. In particular, we'll need to tag side-centered data and cell-centered data independently.
literature
- https://research.nvidia.com/sites/default/files/pubs/2019-03_Cool-Patches%3A-A/Chapter_08.pdf
- TODO: can we easily combine quaternions with FEM to simplify calculation of normal vectors?
- https://mathweb.ucsd.edu/~sbuss/ResearchWeb/spheremean/paper.pdf
- moonolith https://bitbucket.org/zulianp/moonolith/src/master/ implements some intersection algorithms we might want
ryanleaf1000
From a discussion today. This is not something I have time to complete myself but I still want to write down some notes.
The fundamental geometric problem we need to solve is ray intersection: we need to find 1) the points where a finite difference stencil intersects the mesh and 2) the corresponding unique element. If we get those then we can start building IIM algorithms.
There's a lot of interesting things to investigate:
@ryanleaf1000 A good place to start would be to look at how to intersect rays and elements in serial: i.e., given a
fdl::PatchMap, compute a new structure describing where the intersections are with SAMRAI indices, deal.II elements, points on the structure, and points in the reference configuration. We should get this working before doing any physics. This is similar to https://github.com/drwells/fiddle/blob/057b76976b93b3f5c140a1a63b905576de41d26f/source/interaction/interaction_utilities.cc#L202-L269