Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Interpolation to and from sub-meshes for co-dimension 0 #3114

Merged
merged 41 commits into from
Apr 23, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
41 commits
Select commit Hold shift + click to select a range
a51a9fb
First attempt on submesh interpolation
jorgensd Mar 22, 2024
2368272
ADd some more conversion
jorgensd Mar 24, 2024
dca89f4
UPdate more interfaces
jorgensd Mar 24, 2024
4e59abb
Expose python arrays of mappings to C++.
jorgensd Mar 24, 2024
62fd08f
Add more asserts to please mypy
jorgensd Mar 24, 2024
8232486
Fix C++ interface
jorgensd Mar 24, 2024
62a49fa
Set some default zero-spans in C++ + fix hyperelasticity demo
jorgensd Mar 24, 2024
0739027
Set atol depending on default floating type
jorgensd Mar 24, 2024
7111ed1
Use exact solution on submesh when interpolating back to avoid accumu…
jorgensd Mar 24, 2024
71f9adf
Use exact solution in mapping back due to same issue
jorgensd Mar 24, 2024
af38559
Less strict tolerance
jorgensd Mar 24, 2024
5d48051
Merge branch 'main' into dokken/submesh_codim0_interpolation
jorgensd Apr 8, 2024
2fb6ee8
Apply suggestions from code review
jorgensd Apr 8, 2024
eaa9e43
Add default arguments in C++ and add documentation to tests
jorgensd Apr 8, 2024
46f0118
Fix submesh C++ wrapping
jorgensd Apr 8, 2024
72579a9
Renaming + updating docstring
jorgensd Apr 8, 2024
45cf1ed
Fix conversion of u
jorgensd Apr 8, 2024
d176828
More renaming
jorgensd Apr 8, 2024
460acb0
Fix documentation
jorgensd Apr 8, 2024
2f06795
Apply suggestions from code review
jorgensd Apr 10, 2024
ab0641e
Merge branch 'main' into dokken/submesh_codim0_interpolation
jorgensd Apr 12, 2024
a37cb6c
Merge branch 'main' into dokken/submesh_codim0_interpolation
jorgensd Apr 13, 2024
61a4e28
Merge branch 'main' into dokken/submesh_codim0_interpolation
jorgensd Apr 14, 2024
f8052ed
Improve documentation
jorgensd Apr 14, 2024
c5df0d6
Ruff format docstring
jorgensd Apr 14, 2024
72cbdad
Merge branch 'main' into dokken/submesh_codim0_interpolation
jorgensd Apr 16, 2024
0744255
Merge branch 'main' into dokken/submesh_codim0_interpolation
jorgensd Apr 19, 2024
ae5d27f
Merge branch 'main' into dokken/submesh_codim0_interpolation
jhale Apr 19, 2024
abcbdae
Merge branch 'main' into dokken/submesh_codim0_interpolation
jorgensd Apr 19, 2024
ed8aba3
Merge branch 'main' into dokken/submesh_codim0_interpolation
jhale Apr 22, 2024
c3fb557
Merge branch 'main' into dokken/submesh_codim0_interpolation
jorgensd Apr 22, 2024
cce45bd
Try less strict atol
jorgensd Apr 22, 2024
a1a57ac
Modify test to avoid squaring functions
jorgensd Apr 22, 2024
0d390f2
Ruff formatting
jorgensd Apr 22, 2024
4d3f535
Merge branch 'main' into dokken/submesh_codim0_interpolation
jorgensd Apr 22, 2024
2a04e4d
Scatter forward data
jorgensd Apr 22, 2024
a908db1
Modify test again
jorgensd Apr 22, 2024
4cd78fc
Try scaling down gradient
jorgensd Apr 22, 2024
0d81933
Fix test
jorgensd Apr 22, 2024
ae3156b
Merge branch 'main' into dokken/submesh_codim0_interpolation
jhale Apr 23, 2024
4953137
clang-format
jhale Apr 23, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion cpp/demo/hyperelasticity/main.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -264,7 +264,7 @@ int main(int argc, char* argv[])

auto sigma = fem::Function<T>(S);
sigma.name = "cauchy_stress";
sigma.interpolate(sigma_expression);
sigma.interpolate(sigma_expression, *mesh);

// Save solution in VTK format
io::VTKFile file_u(mesh->comm(), "u.pvd", "w");
Expand Down
3 changes: 2 additions & 1 deletion cpp/demo/interpolation_different_meshes/main.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -74,7 +74,8 @@ int main(int argc, char* argv[])
*u_hex->function_space()->mesh(),
*u_hex->function_space()->element(),
*u_tet->function_space()->mesh(), padding);
u_hex->interpolate(*u_tet, nmm_interpolation_data);
constexpr std::span<const std::int32_t> cell_map;
u_hex->interpolate(*u_tet, cell_map, nmm_interpolation_data);

#ifdef HAS_ADIOS2
io::VTXWriter<double> write_tet(mesh_tet->comm(), "u_tet.bp", {u_tet});
Expand Down
64 changes: 51 additions & 13 deletions cpp/dolfinx/fem/Function.h
Original file line number Diff line number Diff line change
Expand Up @@ -153,28 +153,34 @@ class Function
/// @brief Interpolate a provided Function.
/// @param[in] v The function to be interpolated
/// @param[in] cells The cells to interpolate on
/// @param[in] cell_map A map from cells in the mesh associated with `this`
/// function to cells in mesh associated with `v`
/// @param[in] nmm_interpolation_data Auxiliary data to interpolate on
/// nonmatching meshes. This data can be generated with
/// generate_nonmatching_meshes_interpolation_data (optional).
void interpolate(
const Function<value_type, geometry_type>& v,
std::span<const std::int32_t> cells,
std::span<const std::int32_t> cell_map,
const std::tuple<std::span<const std::int32_t>,
std::span<const std::int32_t>,
std::span<const geometry_type>,
std::span<const std::int32_t>>& nmm_interpolation_data
= {})
{
fem::interpolate(*this, v, cells, nmm_interpolation_data);
fem::interpolate(*this, v, cells, cell_map, nmm_interpolation_data);
}

/// @brief Interpolate a provided Function.
/// @param[in] v The function to be interpolated
/// @param[in] cell_map Map from cells in self to cell indices in `v`
/// @param[in] nmm_interpolation_data Auxiliary data to interpolate on
/// nonmatching meshes. This data can be generated with
/// generate_nonmatching_meshes_interpolation_data (optional).
void interpolate(
const Function<value_type, geometry_type>& v,
std::span<const std::int32_t> cell_map = std::span<const std::int32_t>()
= {},
const std::tuple<std::span<const std::int32_t>,
std::span<const std::int32_t>,
std::span<const geometry_type>,
Expand All @@ -184,12 +190,12 @@ class Function
assert(_function_space);
assert(_function_space->mesh());
int tdim = _function_space->mesh()->topology()->dim();
auto cell_map = _function_space->mesh()->topology()->index_map(tdim);
assert(cell_map);
std::int32_t num_cells = cell_map->size_local() + cell_map->num_ghosts();
auto cell_imap = _function_space->mesh()->topology()->index_map(tdim);
assert(cell_imap);
std::int32_t num_cells = cell_imap->size_local() + cell_imap->num_ghosts();
std::vector<std::int32_t> cells(num_cells, 0);
std::iota(cells.begin(), cells.end(), 0);
interpolate(v, cells, nmm_interpolation_data);
interpolate(v, cells, cell_map, nmm_interpolation_data);
}

/// Interpolate an expression function on a list of cells
Expand Down Expand Up @@ -283,8 +289,13 @@ class Function
/// `FiniteElement::interpolation_points()` for the element associated
/// with `u`.
/// @param[in] cells The cells to interpolate on
/// @param[in] expr_mesh The mesh to evaluate the expression on
/// @param[in] cell_map Map from `cells` to cells in expression
void interpolate(const Expression<value_type, geometry_type>& e,
std::span<const std::int32_t> cells)
std::span<const std::int32_t> cells,
const dolfinx::mesh::Mesh<geometry_type>& expr_mesh,
std::span<const std::int32_t> cell_map
= std::span<const std::int32_t>())
{
// Check that spaces are compatible
assert(_function_space);
Expand Down Expand Up @@ -331,8 +342,25 @@ class Function

// Evaluate Expression at points
assert(_function_space->mesh());
e.eval(*_function_space->mesh(), cells, fdata,
{num_cells, num_points * value_size});

std::vector<std::int32_t> cells_expr;
cells_expr.reserve(num_cells);
// Get mesh and check if mesh is shared
if (auto mesh_v = _function_space->mesh();
expr_mesh.topology() == mesh_v->topology())
{
cells_expr.insert(cells_expr.end(), cells.begin(), cells.end());
}
// If meshes are different and input mapping is given
else if (!cell_map.empty())
{
std::transform(cells.begin(), cells.end(), std::back_inserter(cells_expr),
[&cell_map](std::int32_t c) { return cell_map[c]; });
}
else
std::runtime_error("Meshes are different and no cell map is provided");

e.eval(expr_mesh, cells_expr, fdata, {num_cells, num_points * value_size});

// Reshape evaluated data to fit interpolate
// Expression returns matrix of shape (num_cells, num_points *
Expand All @@ -356,17 +384,27 @@ class Function

/// Interpolate an Expression (based on UFL) on all cells
/// @param[in] e The function to be interpolated
void interpolate(const Expression<value_type, geometry_type>& e)
/// @param[in] expr_mesh Mesh the Expression `e` is defined on.
/// @param[in] cell_map Map from `cells` to cells in expression if
/// receiving function is defined on a different mesh than the expression
void interpolate(const Expression<value_type, geometry_type>& e,
const dolfinx::mesh::Mesh<geometry_type>& expr_mesh,
std::span<const std::int32_t> cell_map
jorgensd marked this conversation as resolved.
Show resolved Hide resolved
= std::span<const std::int32_t>() = {})
{
assert(_function_space);
assert(_function_space->mesh());
const int tdim = _function_space->mesh()->topology()->dim();
auto cell_map = _function_space->mesh()->topology()->index_map(tdim);
assert(cell_map);
std::int32_t num_cells = cell_map->size_local() + cell_map->num_ghosts();
auto cell_imap = _function_space->mesh()->topology()->index_map(tdim);
assert(cell_imap);
std::int32_t num_cells = cell_imap->size_local() + cell_imap->num_ghosts();
std::vector<std::int32_t> cells(num_cells, 0);
std::iota(cells.begin(), cells.end(), 0);
interpolate(e, cells);

if (cell_map.size() == 0)
interpolate(e, cells, expr_mesh, cell_map);
else
interpolate(e, cells, expr_mesh, std::span(cells.data(), cells.size()));
}

/// @brief Evaluate the Function at points.
Expand Down
Loading
Loading