Skip to content

Commit

Permalink
Simplify
Browse files Browse the repository at this point in the history
  • Loading branch information
@garth-wells
garth-wells committed Jun 18, 2024
1 parent a5fa3d6 commit f66c2dd
Show file tree
Hide file tree
Showing 4 changed files with 60 additions and 69 deletions.
85 changes: 34 additions & 51 deletions cpp/Contact.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,11 +36,8 @@ dolfinx::graph::AdjacencyList<std::int32_t> create_cell_facet_pairs(
// Get the facets ids with tag links[i]
std::vector<std::int32_t> facets = marker->find(tag);

// Surface index
// int index = surfaces.offsets()[s] + i;

// Compute the (cell, local index) pairs for each facet
auto [cell_facet_pairs, num_local] = compute_active_entities(
std::vector<std::int32_t> cell_facet_pairs = compute_active_entities(
mesh, facets, dolfinx::fem::IntegralType::exterior_facet);

// Add to list of (cell, local index) pairs
Expand All @@ -50,9 +47,6 @@ dolfinx::graph::AdjacencyList<std::int32_t> create_cell_facet_pairs(

// Add to offset
offsets.push_back(offsets.back() + cell_facet_pairs.size());

// Store number of owned facets in this surface
// _local_facets[index] = num_local;
}
}

Expand Down Expand Up @@ -122,13 +116,10 @@ void compute_linked_cells(
} // namespace

//----------------------------------------------------------------------------
Contact::Contact(
const std::vector<
std::shared_ptr<const dolfinx::mesh::MeshTags<std::int32_t>>>& markers,
const dolfinx::graph::AdjacencyList<std::int32_t>& surfaces,
const std::vector<std::array<int, 2>>& contact_pairs,
std::shared_ptr<const dolfinx::mesh::Mesh<double>> mesh,
const std::vector<ContactMode>& mode, int q_deg)
Contact::Contact(const dolfinx::graph::AdjacencyList<std::int32_t>& surfaces,
const std::vector<std::array<int, 2>>& contact_pairs,
std::shared_ptr<const dolfinx::mesh::Mesh<double>> mesh,
const std::vector<ContactMode>& mode, int q_deg)
: _quadrature_rule(mesh->topology()->cell_type(), q_deg,
mesh->topology()->dim() - 1,
basix::quadrature::type::Default),
Expand All @@ -137,40 +128,10 @@ Contact::Contact(
_reference_contact_points(contact_pairs.size()),
_reference_contact_shape(contact_pairs.size()),
_qp_phys(surfaces.array().size()), _max_links(contact_pairs.size()),
_cell_facet_pairs(create_cell_facet_pairs(*mesh, markers, surfaces)),
_submesh(*mesh, _cell_facet_pairs.array()),
_local_facets(surfaces.array().size()), _mode(mode)
_cell_facet_pairs(surfaces), _submesh(*mesh, surfaces.array()),
_mode(mode)
{
assert(_mesh);

if (markers.size() != (std::size_t)surfaces.num_nodes())
throw std::runtime_error("maker and surfaces have different sizes.");

// for (int i = 0; i < surfaces.num_nodes(); ++i)
// {
// if (surfaces.num_links(i) != 1)
// throw std::runtime_error("Num surfces links !=0 : "
// + std::to_string(surfaces.num_links(i)));
// }

// TODO: remove the below code, compute num_local directly from
// _cell_facet_pairs
for (std::size_t s = 0; s < markers.size(); ++s)
{
std::shared_ptr<const dolfinx::mesh::MeshTags<int>> marker = markers[s];
std::span<const int> links = surfaces.links(int(s));
for (std::size_t i = 0; i < links.size(); ++i)
{
std::vector<std::int32_t> facets = marker->find(links[i]);
auto [cell_facet_pairs, num_local] = compute_active_entities(
*mesh, facets, dolfinx::fem::IntegralType::exterior_facet);

// store how many facets are owned by the process
int index = surfaces.offsets()[s] + i;
_local_facets[index] = num_local;
}
}

auto topology = mesh->topology();
std::int32_t num_cells = topology->index_map(topology->dim())->size_local();
std::vector<std::size_t> _num_local_facets;
Expand All @@ -181,19 +142,41 @@ Contact::Contact(
// Count the number of owned (cells, local facet index pairs). List
// is sorted by cell index, so we could backwards from the end of
// the list for efficiency.
_num_local_facets.push_back(0);
_local_facets.push_back(0);
for (std::int32_t i = (std::int32_t)facets.size() - 2; i >= 0; i -= 2)
{
if (facets[i] < num_cells)
{
_num_local_facets.back() = i / 2 + 1;
_local_facets.back() = i / 2 + 1;
break;
}
}
}

if (_local_facets != _num_local_facets)
throw std::runtime_error("data mis-match");
}
//----------------------------------------------------------------------------
Contact::Contact(
const std::vector<
std::shared_ptr<const dolfinx::mesh::MeshTags<std::int32_t>>>& markers,
const dolfinx::graph::AdjacencyList<std::int32_t>& surfaces,
const std::vector<std::array<int, 2>>& contact_pairs,
std::shared_ptr<const dolfinx::mesh::Mesh<double>> mesh,
const std::vector<ContactMode>& mode, int q_deg)
: Contact(create_cell_facet_pairs(*mesh, markers, surfaces), contact_pairs,
mesh, mode, q_deg)
// _quadrature_rule(mesh->topology()->cell_type(), q_deg,
// mesh->topology()->dim() - 1,
// basix::quadrature::type::Default),
// _contact_pairs(contact_pairs), _mesh(mesh),
// _facet_maps(contact_pairs.size()),
// _reference_contact_points(contact_pairs.size()),
// _reference_contact_shape(contact_pairs.size()),
// _qp_phys(surfaces.array().size()), _max_links(contact_pairs.size()),
// _cell_facet_pairs(create_cell_facet_pairs(*mesh, markers, surfaces)),
// _submesh(*mesh, _cell_facet_pairs.array()),
// _local_facets(surfaces.array().size()), _mode(mode)
{
if (markers.size() != (std::size_t)surfaces.num_nodes())
throw std::runtime_error("maker and surfaces have different sizes.");
}
//----------------------------------------------------------------------------
std::pair<std::vector<double>, std::array<std::size_t, 3>>
Expand Down
20 changes: 18 additions & 2 deletions cpp/Contact.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,23 @@ class Contact
public:
/// Constructor
///
/// @todo Simplify the constructor arguments
/// @param[in] surfaces The facets for each surfaces, stored as
/// flattened (cell index, local facet index) pairs for each surface.
/// The facet pairs for surface `i` are returned by
/// `surface.links(i)`.
/// @param[in] contact_pairs list of pairs `(i, j)` marking the `i`th
/// and `j`th surface in surfaces->array() as a contact pair
/// @param[in] mesh TODO
/// @param[in] mode Contact detection algorithm for each pair
/// @param[in] q_deg The quadrature degree.
Contact(const dolfinx::graph::AdjacencyList<std::int32_t>& surfaces,
const std::vector<std::array<int, 2>>& contact_pairs,
std::shared_ptr<const dolfinx::mesh::Mesh<double>> mesh,
const std::vector<ContactMode>& mode, int q_deg = 3);

/// Constructor
///
/// @todo Remove this constructor
///
/// @param[in] markers Tags `markers[i]` for surface `i`.
/// @param[in] surfaces For surface `i`, `surfaces.links(i)` gives the
Expand Down Expand Up @@ -305,7 +321,7 @@ class Contact
// quadrature rule
QuadratureRule _quadrature_rule;

// store index of candidate_surface for each quadrature_surface
// store index of candidate_surface for each quadrature_surface
std::vector<std::array<int, 2>> _contact_pairs;

std::shared_ptr<const dolfinx::mesh::Mesh<double>> _mesh;
Expand Down
22 changes: 7 additions & 15 deletions cpp/utils.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -655,28 +655,24 @@ dolfinx_contact::get_update_normal(
}
}
//----------------------------------------------------------------------------
std::pair<std::vector<std::int32_t>, std::size_t>
dolfinx_contact::compute_active_entities(
std::vector<std::int32_t> dolfinx_contact::compute_active_entities(
const dolfinx::mesh::Mesh<double>& mesh,
std::span<const std::int32_t> entities, dolfinx::fem::IntegralType integral)
{

int tdim = mesh.topology()->dim();
const int size_local = mesh.topology()->index_map(tdim)->size_local();
switch (integral)
{
case dolfinx::fem::IntegralType::cell:
{
std::vector<std::int32_t> active_entities(entities.begin(), entities.end());
dolfinx::radix_sort(std::span<std::int32_t>(active_entities));
auto cell_it = std::upper_bound(active_entities.begin(),
active_entities.end(), size_local);
std::size_t num_local = std::distance(active_entities.begin(), cell_it);
return std::make_pair(active_entities, num_local);
return active_entities;
}
case dolfinx::fem::IntegralType::exterior_facet:
{
auto topology = mesh.topology();
assert(topology);
int tdim = mesh.topology()->dim();

auto f_to_c = topology->connectivity(tdim - 1, tdim);
assert(f_to_c);
auto c_to_f = topology->connectivity(tdim, tdim - 1);
Expand All @@ -701,28 +697,24 @@ dolfinx_contact::compute_active_entities(
dolfinx::argsort_radix<std::int32_t>(cells, perm);

// Sort cells in ascending order
std::size_t num_local = 0;
std::vector<std::int32_t> active_entities;
active_entities.reserve(2 * entities.size());
for (auto p : perm)
{
active_entities.push_back(cells[p]);
active_entities.push_back(facets[p]);
if (cells[p] < size_local)
num_local += 1;
}

return {std::move(active_entities), num_local};
return active_entities;
}
default:
throw std::invalid_argument(
"Integral type not supported. Note that this function "
"has not been implemented for interior facets.");
}

return {};
// return {};
}

//-------------------------------------------------------------------------------------
dolfinx::graph::AdjacencyList<std::int32_t>
dolfinx_contact::entities_to_geometry_dofs(
Expand Down
2 changes: 1 addition & 1 deletion cpp/utils.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -270,7 +270,7 @@ get_update_normal(const dolfinx::fem::CoordinateElement<double>& cmap);
/// each facet in `entities`, the corresponding (cell, local facet) pair
/// is computed, where "local facet" is the index if the facet relative
/// to the cell. Returned list if sorted by cell index.
std::pair<std::vector<std::int32_t>, std::size_t>
std::vector<std::int32_t>
compute_active_entities(const dolfinx::mesh::Mesh<double>& mesh,
std::span<const std::int32_t> entities,
dolfinx::fem::IntegralType integral);
Expand Down

0 comments on commit f66c2dd

Please sign in to comment.