tether force evaluation in IIMethod

[cpuelz]

IBAMR/examples/IIM/ex9/example.cpp

Lines 212 to 225 in 79502c7

	for (unsigned int d = 0; d < NDIM; ++d)
	{
	const auto el_begin = mesh_bndry.active_local_elements_begin();
	const auto el_end = mesh_bndry.active_local_elements_end();
	for (auto el_it = el_begin; el_it != el_end; ++el_it)
	{
	const Elem* elem_bndry = *el_it;
	if ((elem_bndry->contains_point(X)) && !((vol_beam_bndry_info->has_boundary_id(elem, side, 0))))
	{
	F(d) = Tau_new_beam_surface_system->point_value(d, X, elem_bndry); //&side_elem);
	}
	}
	}

Any thoughts about how to make this evaluation faster? This seems to be a bottleneck for large surface meshes.

[cpuelz]

not sure if there is a way to avoid looping over local surface mesh elements during each evaluation.

[boyceg]

We must avoid looping over all boundary elements just to evaluate a single point force. Somehow we need either to rely on existing libMesh data structures to map between surface and volume elements, or we need to cache that stuff ourselves.

[drwells]

An easy partial fix is to exchange the loop order: we can make the loop over dimensions the innermost one to cut the amount of work by 2/3rds.

[AminKolah]

We must avoid looping over all boundary elements just to evaluate a single point force. Somehow we need either to rely on existing libMesh data structures to map between surface and volume elements, or we need to cache that stuff ourselves.

We must avoid looping over all boundary elements just to evaluate a single point force. Somehow we need either to rely on existing libMesh data structures to map between surface and volume elements, or we need to cache that stuff ourselves.

But this is not a single point force. In fact, the traction is imposed on all side-sets of the flexible flapping wing in this example except for a single smaller side that has a zero displacement (clamped to the housing).

[AminKolah]

An easy partial fix is to exchange the loop order: we can make the loop over dimensions the innermost one to cut the amount of work by 2/3rds.

Yes, this would be helpful, especially if you a dense TET surface mesh with a lot of elements.

[AminKolah]

Sorry by "TET", I meant surface triangulation based on a volumetric tetrahedral mesh

[cpuelz]

I guess we can store the info we get from the interior_parent() function beforehand so we have a one-to-one correspondence between volumetric elements on the surface and their corresponding surface elements.

[AminKolah]

I guess we can store the info we get from the interior_parent() function beforehand so we have a one-to-one correspondence between volumetric elements on the surface and their corresponding surface elements.

Okay I guess I'm now getting what @boyceg was also referring to. Such function would be somehow the inverse of interior_parent(). I guess it's possible to create a mapping between the two representations and reuse it throughout. I found the contains_point to be the convenient choice at the time.

[boyceg]

I guess we can store the info we get from the interior_parent() function beforehand so we have a one-to-one correspondence between volumetric elements on the surface and their corresponding surface elements.

It seems like you want two mappings: surface element to {volumetric element + side} and {volumetric element + side} to surface element. Should be easy to track if libMesh does not already provide it.

[boyceg]

I guess we can store the info we get from the interior_parent() function beforehand so we have a one-to-one correspondence between volumetric elements on the surface and their corresponding surface elements.

Okay I guess I'm now getting what @boyceg was also referring to. Such function would be somehow the inverse of interior_parent(). I guess it's possible to create a mapping between the two representations and reuse it throughout. I found the contains_point to be the convenient choice at the time.

The implementation is OK for 2D models, and using contains_point is simpler than setting up the mappings. The problem is that in 3D, the runtime complexity of looping over all surface elements dominates all other operators. Briefly, if you are using an N-by-N-by-N Cartesian grid, then there are O(N^2) surface elements, and this function needs to be evaluated at least once per surface element in each time step. Consequently, the runtime complexity is O(N^4) --- which is asymptotically worse than solving the 3D incompressible Navier-Stokes equations.

[boyceg]

(In contrast, in 2D, there are O(N) surface elements, and so the runtime complexity is O(N^2) --- the same as solving the fluid equations --- and if there are not too many surface elements in a model, you might not even notice this function.)

[AminKolah]

I guess we can store the info we get from the interior_parent() function beforehand so we have a one-to-one correspondence between volumetric elements on the surface and their corresponding surface elements.

Okay I guess I'm now getting what @boyceg was also referring to. Such function would be somehow the inverse of interior_parent(). I guess it's possible to create a mapping between the two representations and reuse it throughout. I found the contains_point to be the convenient choice at the time.

The implementation is OK for 2D models, and using contains_point is simpler than setting up the mappings. The problem is that in 3D, the runtime complexity of looping over all surface elements dominates all other operators. Briefly, if you are using an N-by-N-by-N Cartesian grid, then there are O(N^2) surface elements, and this function needs to be evaluated at least once per surface element in each time step. Consequently, the runtime complexity is O(N^4) --- which is asymptotically worse than solving the 3D incompressible Navier-Stokes equations.

This makes perfect sense, particularly for situations in which the structure expands through the entire Cartesian domain (such as a flexible tube). In example 9, however, the flexible wing where we need to apply the fluid traction is only a small portion of the computational domain, so maybe that's why it hasn't been too bad to loop over the entire surface elements for this particular example.

[boyceg]

Yes, exactly --- for smaller 2D models this is not a big deal. But I think that @cpuelz is trying to use this same code in a larger 3D model and is running into performance problems.

[cpuelz]

so I built an unordered map like this:

        const auto el_begin = bdry_mesh.active_local_elements_begin();
        const auto el_end = bdry_mesh.active_local_elements_end();
        for (auto el_it = el_begin; el_it != el_end; ++el_it)
        {
            const Elem* surface_elem = *el_it;
            const Elem* volume_elem = surface_elem->interior_parent();
            const libMesh::dof_id_type surface_elem_id = surface_elem->id();
            const libMesh::dof_id_type volume_elem_id = volume_elem->id();
            volume_id_to_surface_id[volume_elem_id] = surface_elem_id;
        }

and I evaluate the tethering force function on the volumetric mesh like this:

    {        
        MeshBase& mesh_bndry = bndry_G_systems->get_mesh();
        
        const libMesh::dof_id_signed_type volume_elem_id = elem->id();
        const libMesh::dof_id_signed_type surface_elem_id = volume_id_to_surface_id[volume_elem_id];
        const Elem* elem_bndry = mesh_bndry.elem(surface_elem_id);
        TBOX_ASSERT(elem_bndry->contains_point(X));
        
        for (unsigned int d = 0; d < NDIM; ++d)
        {
            F(d) = Tau_new_surface_system->point_value(d, X, elem_bndry);
        }
            
        return;
    } // tether_force_function

but the assertion I have there is thrown. I suppose either I am not building the map correctly or not using the correct element labeling.

[boyceg]

There could be multiple side elements associated with the same volume element.

[boyceg]

I think you might need to make the key {volume element ID, side index}

[cpuelz]

I call the interior_parent() function for each surface element to get the volume element ID key. I would think that key would be unique for each surface element.

[boyceg]

Yes but you could have multiple surface elements associated with each interior parent -- e.g. the map could be pointing back to the "wong" surface element.

[cpuelz]

disregard that comment. I can visualize a case where you have multiple sides associated with a volume element. I was thinking that would not be the case because my surface mesh is just the vessel lumen, but maybe something funny is happening.

[abarret]

If you have a boundary mesh that is synced with the volumetric mesh via the BoundaryInfo object, it has functions to map between nodes in the interior and nodes on the boundary.
https://mooseframework.inl.gov/docs/doxygen/libmesh/classlibMesh_1_1BoundaryInfo.html#a060ab76ef9acf346dc59ef93ac71411d

[boyceg]

If you want to maintain your own mapping, multimap should work --- you would need to use contains_point only if there is more than one side element associated with the volume element.

[drwells]

Do one of you want to try to turn this into a patch?

[cpuelz]

I have some code I could put in to this example that I think works.