IIM roadmap

[drwells]

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:

1. how do we handle small elements which do not have any corresponding finite difference stencils?
2. how can we do higher-order elements where the intersection code is really hard?
3. what happens if a stencil intersects multiple elements?

@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

fiddle/source/interaction/interaction_utilities.cc

Lines 202 to 269 in 057b769

	template <int dim, int spacedim, typename Scalar>
	void
	count_quadrature_points_internal(
	const int qp_data_idx,
	PatchMap<dim, spacedim> &patch_map,
	const Mapping<dim, spacedim> &position_mapping,
	const std::vector<unsigned char> &quadrature_indices,
	const std::vector<Quadrature<dim>> &quadratures)
	{
	check_quadratures(quadrature_indices,
	quadratures,
	patch_map.get_triangulation());
	// We probably don't need more than 16 quadrature rules
	boost::container::small_vector<std::unique_ptr<FEValues<dim, spacedim>>, 16>
	all_position_fe_values;
	// PatchMap only supports looping over DoFHandler iterators, so we need to
	// make one and never use it explicitly
	const Triangulation<dim, spacedim> &tria = patch_map.get_triangulation();
	// No mixed meshes yet
	Assert(tria.get_reference_cells().size() == 1, ExcNotImplemented());
	const ReferenceCell reference_cell = tria.get_reference_cells().front();
	FE_Nothing<dim, spacedim> fe_nothing(reference_cell);
	DoFHandler<dim, spacedim> dof_handler(tria);
	dof_handler.distribute_dofs(fe_nothing);
	for (const Quadrature<dim> &quad : quadratures)
	{
	all_position_fe_values.emplace_back(
	std::make_unique<FEValues<dim, spacedim>>(
	position_mapping, fe_nothing, quad, update_quadrature_points));
	}
	for (unsigned int patch_n = 0; patch_n < patch_map.size(); ++patch_n)
	{
	auto patch = patch_map.get_patch(patch_n);
	tbox::Pointer<pdat::CellData<spacedim, Scalar>> qp_data =
	patch->getPatchData(qp_data_idx);
	Assert(qp_data, ExcMessage("Type mismatch"));
	Assert(qp_data->getDepth() == 1, ExcMessage("depth should be 1"));
	const hier::Box<spacedim> &patch_box = patch->getBox();
	tbox::Pointer<geom::CartesianPatchGeometry<spacedim>> patch_geom =
	patch->getPatchGeometry();
	Assert(patch_geom, ExcMessage("Type mismatch"));
	auto iter = patch_map.begin(patch_n, dof_handler);
	const auto end = patch_map.end(patch_n, dof_handler);
	for (; iter != end; ++iter)
	{
	const auto cell = *iter;
	const auto quad_index =
	quadrature_indices[cell->active_cell_index()];
	FEValues<dim, spacedim> &position_fe_values =
	*all_position_fe_values[quad_index];
	position_fe_values.reinit(cell);
	for (const Point<spacedim> &q_point :
	position_fe_values.get_quadrature_points())
	{
	const hier::Index<spacedim> i =
	IBTK::IndexUtilities::getCellIndex(q_point,
	patch_geom,
	patch_box);
	if (patch_box.contains(i))
	(*qp_data)(i) += Scalar(1);
	}
	}
	}
	}

[drwells]

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

1. 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
2. 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
3. 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:

1. A class which inherits from IBStrategy which is solely responsible for book-keeping and hooking into the rest of IBAMR.
2. Data structures which efficiently store and grant access via iteration to intersections.
3. 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 to NodalPatchMaps. 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:

1. We could set up an rtree and then intersect each stencil with the rtree:
   that would quickly tell us which elements it might intersect
2. 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.
3. 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:

1. get the physical-space bounding box for an element
2. convert it into an index-space bounding box
3. 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 an FESystem<dim, spacedim>(FE_DGQ<dim, spacedim>(0), spacedim) or FESystem<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]

Hello David,

I pretty much finished testing intersection code. Here are the two functions that I want to add to fiddle:
bool
intersect_line_with_edge(std::vector<std::pair<double, Point<1>> >& t_vals,
const DoFHandler<1, 2>::active_cell_iterator elem,
const dealii::MappingQ<1, 2> &mapping,
dealii::Point<2> r,
dealii::Tensor<1,2> q,
const double tol = 0.0);
bool
intersect_line_with_face(std::vector<std::pair<double, Point<2>> >& t_vals,
const typename DoFHandler<2, 3>::active_cell_iterator elem,
const dealii::MappingQ<2, 3> &mapping,
dealii::Point<3> r,
dealii::Tensor<1,3> q,
const double tol=0);
Where would you prefer me to add these two functions to? grid_utilities.cc? In the two functions, I do loop over interface elements.

[drwells]

Thanks. How about we put these in interaction/interaction_utilities.h and open a pull request?

[ryanleaf1000]

Hello David,
I opened a pull request yesterday, although I don't think it's ready to be merged. I believe that there are likely still modifications needed for it to be completely ready for force spreading. But please feel free to share your comments, thoughts, or any preferences you would prefer me to modify?

[drwells]

Lets move that discussion over there - this issue is just the roadmap.