tether force evaluation in ```IIMethod``` about ibamr HOT 23 OPEN

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.

from ibamr.

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.

from ibamr.

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.

from ibamr.

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

from ibamr.

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

from ibamr.

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.

from ibamr.

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

from ibamr.

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.

from ibamr.

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).

from ibamr.

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.

from ibamr.

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

from ibamr.

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.

from ibamr.

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.

from ibamr.

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.

from ibamr.

(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.)

from ibamr.

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.

from ibamr.

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.

from ibamr.

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.

from ibamr.

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.

from ibamr.

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

from ibamr.

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.

from ibamr.

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

from ibamr.

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

from ibamr.