Skip to content

Commit

Permalink
Remove quadrature setter
Browse files Browse the repository at this point in the history
  • Loading branch information
@garth-wells
garth-wells committed Jun 16, 2024
1 parent 754b197 commit 61c0b86
Show file tree
Hide file tree
Showing 3 changed files with 46 additions and 38 deletions.
62 changes: 37 additions & 25 deletions cpp/Contact.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,10 @@ dolfinx_contact::Contact::Contact(
const std::vector<std::array<int, 2>>& contact_pairs,
std::shared_ptr<const dolfinx::mesh::Mesh<double>> mesh,
const std::vector<ContactMode>& mode, const int q_deg)
: _surfaces(surfaces.array()), _contact_pairs(contact_pairs), _mesh(mesh),
: _quadrature_rule(mesh->topology()->cell_type(), q_deg,
mesh->topology()->dim() - 1,
basix::quadrature::type::Default),
_surfaces(surfaces.array()), _contact_pairs(contact_pairs), _mesh(mesh),
_mode(mode)
{
std::size_t num_surfaces = surfaces.array().size();
Expand All @@ -96,21 +99,28 @@ dolfinx_contact::Contact::Contact(
std::shared_ptr<const dolfinx::graph::AdjacencyList<int>> c_to_f
= topology->connectivity(tdim, fdim);
assert(c_to_f);

// storing number of facets owned locally for each surface
_local_facets.resize(num_surfaces);

// used to store map from quadrature to candidate surface for each contact
// pair
_facet_maps.resize(contact_pairs.size());

// store physical quadrature points for each surface
_qp_phys.resize(num_surfaces);

// reference points on opposite surface
_reference_contact_points.resize(contact_pairs.size());

// shape of reference points on opposite surface
_reference_contact_shape.resize(contact_pairs.size());

// store max number of links for each quadrature surface
_max_links.resize(contact_pairs.size());
// Create adjacency list linking facets as (cell, facet) pairs to the index of
// the surface. The pairs are flattened row-major

// Create adjacency list linking facets as (cell, facet) pairs to the
// index of the surface. The pairs are flattened row-major
std::vector<std::int32_t> all_facet_pairs;
std::vector<std::int32_t> offsets;
offsets.push_back(0);
Expand All @@ -133,19 +143,21 @@ dolfinx_contact::Contact::Contact(
}
}
_submesh = SubMesh(*mesh, all_facet_pairs);

_cell_facet_pairs
= std::make_shared<dolfinx::graph::AdjacencyList<std::int32_t>>(
std::move(all_facet_pairs), std::move(offsets));
_quadrature_rule = std::make_shared<QuadratureRule>(
topology->cell_type(), q_deg, fdim, basix::quadrature::type::Default);

// _quadrature_rule = std::make_shared<QuadratureRule>(
// topology->cell_type(), q_deg, fdim, basix::quadrature::type::Default);
}
//----------------------------------------------------------------------------
std::pair<std::vector<double>, std::array<std::size_t, 3>>
dolfinx_contact::Contact::qp_phys(int surface) const
{
const std::size_t num_facets = _local_facets[surface];
const std::size_t num_q_points
= _quadrature_rule->offset()[1] - _quadrature_rule->offset()[0];
= _quadrature_rule.offset()[1] - _quadrature_rule.offset()[0];
const std::size_t gdim = _mesh->geometry().dim();
std::array<std::size_t, 3> shape = {num_facets, num_q_points, gdim};
return {_qp_phys[surface], shape};
Expand All @@ -166,7 +178,7 @@ std::size_t dolfinx_contact::Contact::coefficients_size(
error::check_cell_type(_mesh->topology()->cell_type());

const std::size_t num_q_points
= _quadrature_rule->offset()[1] - _quadrature_rule->offset()[0];
= _quadrature_rule.offset()[1] - _quadrature_rule.offset()[0];
const std::size_t max_links
= *std::max_element(_max_links.begin(), _max_links.end());

Expand Down Expand Up @@ -285,7 +297,7 @@ void dolfinx_contact::Contact::create_distance_map(int pair)
// Compute facet map
[[maybe_unused]] auto [adj, reference_x, shape] = compute_distance_map(
*quadrature_mesh, quadrature_facets, *candidate_mesh, submesh_facets,
*_quadrature_rule, _mode[pair], _radius);
_quadrature_rule, _mode[pair], _radius);

_facet_maps[pair]
= std::make_shared<dolfinx::graph::AdjacencyList<std::int32_t>>(adj);
Expand All @@ -297,15 +309,15 @@ void dolfinx_contact::Contact::create_distance_map(int pair)
const dolfinx::fem::CoordinateElement<double>& cmap
= candidate_mesh->geometry().cmap();
std::tie(_reference_basis, _reference_shape)
= tabulate(cmap, *_quadrature_rule);
= tabulate(cmap, _quadrature_rule);

// NOTE: This function should be moved somwhere else, or return the actual
// points such that we compuld send them in to compute_distance_map.
// Compute quadrature points on physical facet _qp_phys_"origin_meshtag"
create_q_phys(quadrature_mt);

// Update maximum number of connected cells
_max_links[pair] = _quadrature_rule->num_points(0);
_max_links[pair] = _quadrature_rule.num_points(0);
}
//------------------------------------------------------------------------------------------------
std::pair<std::vector<PetscScalar>, int>
Expand Down Expand Up @@ -333,7 +345,7 @@ dolfinx_contact::Contact::pack_nx(int pair) const

// num quadrature pints
error::check_cell_type(quadrature_mesh->topology()->cell_type());
const std::size_t num_q_points = _quadrature_rule->num_points(0);
const std::size_t num_q_points = _quadrature_rule.num_points(0);
std::vector<PetscScalar> normals(num_facets * num_q_points * gdim, 0.0);
const int cstride = (int)num_q_points * gdim;

Expand All @@ -342,8 +354,8 @@ dolfinx_contact::Contact::pack_nx(int pair) const
return {std::move(normals), cstride};

// Get all quadrature points
const std::vector<double>& q_points = _quadrature_rule->points();
assert(_quadrature_rule->tdim() == (std::size_t)tdim);
const std::vector<double>& q_points = _quadrature_rule.points();
assert(_quadrature_rule.tdim() == (std::size_t)tdim);
const std::array<std::size_t, 2> shape
= {q_points.size() / tdim, (std::size_t)tdim};

Expand Down Expand Up @@ -418,7 +430,7 @@ dolfinx_contact::Contact::generate_kernel(
{
std::size_t max_links
= *std::max_element(_max_links.begin(), _max_links.end());
return generate_contact_kernel(type, V, *_quadrature_rule, max_links);
return generate_contact_kernel(type, V, _quadrature_rule, max_links);
}
//------------------------------------------------------------------------------------------------
void dolfinx_contact::Contact::create_q_phys(int origin_meshtag)
Expand All @@ -429,7 +441,7 @@ void dolfinx_contact::Contact::create_q_phys(int origin_meshtag)
= _submesh.get_submesh_tuples(_cell_facet_pairs->links(origin_meshtag));
auto mesh_sub = _submesh.mesh();
const std::size_t gdim = mesh_sub->geometry().dim();
const std::vector<size_t>& qp_offsets = _quadrature_rule->offset();
const std::vector<size_t>& qp_offsets = _quadrature_rule.offset();
_qp_phys[origin_meshtag].resize((qp_offsets[1] - qp_offsets[0])
* (submesh_facets.size() / 2) * gdim);
compute_physical_points(
Expand Down Expand Up @@ -494,7 +506,7 @@ dolfinx_contact::Contact::pack_gap(int pair) const
// NOTE: Assumes same number of quadrature points on all facets
error::check_cell_type(candidate_mesh->topology()->cell_type());
const std::size_t num_q_point
= _quadrature_rule->offset()[1] - _quadrature_rule->offset()[0];
= _quadrature_rule.offset()[1] - _quadrature_rule.offset()[0];
const dolfinx::mesh::Geometry<double>& geometry = candidate_mesh->geometry();
const int gdim = geometry.dim();

Expand Down Expand Up @@ -646,7 +658,7 @@ dolfinx_contact::Contact::pack_test_functions(
const std::size_t max_links
= *std::max_element(_max_links.begin(), _max_links.end());
const std::size_t num_q_points
= _quadrature_rule->offset()[1] - _quadrature_rule->offset()[0];
= _quadrature_rule.offset()[1] - _quadrature_rule.offset()[0];
const std::size_t bs = element->block_size();
const auto cstride = int(num_q_points * max_links * b_shape[2] * bs);
std::vector<PetscScalar> cb(
Expand Down Expand Up @@ -737,7 +749,7 @@ void dolfinx_contact::Contact::crop_invalid_points(std::size_t pair,
auto [quadrature_mt, candidate_mt] = _contact_pairs[pair];
const std::size_t num_facets = _local_facets[quadrature_mt];
const std::size_t num_q_points
= _quadrature_rule->offset()[1] - _quadrature_rule->offset()[0];
= _quadrature_rule.offset()[1] - _quadrature_rule.offset()[0];
const std::size_t gdim = _mesh->geometry().dim();
std::shared_ptr<const dolfinx::graph::AdjacencyList<std::int32_t>>
candidate_map = _facet_maps[pair];
Expand Down Expand Up @@ -806,7 +818,7 @@ dolfinx_contact::Contact::pack_u_contact(

// Output vector
const std::size_t num_q_points
= _quadrature_rule->offset()[1] - _quadrature_rule->offset()[0];
= _quadrature_rule.offset()[1] - _quadrature_rule.offset()[0];
std::vector<PetscScalar> c(num_facets * num_q_points * bs_element, 0.0);
const auto cstride = int(num_q_points * bs_element);

Expand Down Expand Up @@ -921,14 +933,14 @@ dolfinx_contact::Contact::pack_gap_plane(int pair, double g)
const dolfinx::fem::CoordinateElement<double>& cmap
= _mesh->geometry().cmap();
std::tie(_reference_basis, _reference_shape)
= tabulate(cmap, *_quadrature_rule);
= tabulate(cmap, _quadrature_rule);

// Compute quadrature points on physical facet _qp_phys_"quadrature_mt"
create_q_phys(quadrature_mt);
const std::size_t num_facets = _local_facets[quadrature_mt];
// FIXME: This does not work for prism meshes
std::size_t num_q_point
= _quadrature_rule->offset()[1] - _quadrature_rule->offset()[0];
= _quadrature_rule.offset()[1] - _quadrature_rule.offset()[0];
mdspan_t<const double, 3> qp_span(_qp_phys[quadrature_mt].data(), num_facets,
num_q_point, gdim);
std::vector<PetscScalar> c(num_facets * num_q_point * gdim, 0.0);
Expand Down Expand Up @@ -964,7 +976,7 @@ dolfinx_contact::Contact::pack_ny(int pair) const
// create coefficient vector
const dolfinx::mesh::Geometry<double>& geometry = candidate_mesh->geometry();
const int gdim = geometry.dim();
const std::size_t num_q_points = _quadrature_rule->num_points(0);
const std::size_t num_q_points = _quadrature_rule.num_points(0);
std::vector<PetscScalar> normals(num_facets * num_q_points * gdim, 0.0);
const int cstride = (int)num_q_points * gdim;

Expand Down Expand Up @@ -1347,7 +1359,7 @@ dolfinx_contact::Contact::pack_grad_test_functions(
= _facet_maps[pair];
const std::size_t num_facets = _local_facets[quadrature_mt];
const std::size_t num_q_points
= _quadrature_rule->offset()[1] - _quadrature_rule->offset()[0];
= _quadrature_rule.offset()[1] - _quadrature_rule.offset()[0];
std::vector<double> q_points(std::size_t(num_q_points) * std::size_t(gdim));
mdspan_t<const double, 3> qp_span(_qp_phys[quadrature_mt].data(), num_facets,
num_q_points, gdim);
Expand Down Expand Up @@ -1465,7 +1477,7 @@ dolfinx_contact::Contact::pack_grad_u_contact(
= _facet_maps[pair];
const std::size_t num_facets = _local_facets[quadrature_mt];
const std::size_t num_q_points
= _quadrature_rule->offset()[1] - _quadrature_rule->offset()[0];
= _quadrature_rule.offset()[1] - _quadrature_rule.offset()[0];
// Output vector
const auto cstride = int(num_q_points * bs_element * gdim);
std::vector<PetscScalar> c(num_facets * cstride, 0.0);
Expand Down Expand Up @@ -1556,6 +1568,6 @@ void dolfinx_contact::Contact::update_submesh_geometry(
//-----------------------------------------------------------------------------------------------
std::size_t dolfinx_contact::Contact::num_q_points() const
{
return _quadrature_rule->offset()[1] - _quadrature_rule->offset()[0];
return _quadrature_rule.offset()[1] - _quadrature_rule.offset()[0];
}
//-----------------------------------------------------------------------------------------------
20 changes: 9 additions & 11 deletions cpp/Contact.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,8 +35,8 @@ namespace dolfinx_contact
class Contact
{
public:
/// empty constructor
Contact() = default;
// /// empty constructor
// Contact() = default;

/// Constructor
/// @param[in] markers List of meshtags defining the contact surfaces
Expand All @@ -52,7 +52,9 @@ class 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, const int q_deg = 3);
const std::vector<ContactMode>& mode,
// QuadratureRule q_rule,
const int q_deg = 3);

/// Return meshtag value for surface with index surface
/// @param[in] surface - the index of the surface
Expand All @@ -78,14 +80,8 @@ class Contact
/// @return TODO
std::size_t local_facets(int s) const { return _local_facets[s]; }

/// set quadrature rule
void set_quadrature_rule(QuadratureRule q_rule)
{
_quadrature_rule = std::make_shared<QuadratureRule>(q_rule);
}

/// return quadrature rule
const QuadratureRule& quadrature_rule() const { return *_quadrature_rule; }
const QuadratureRule& quadrature_rule() const { return _quadrature_rule; }

/// @brief TODO
/// @return TODO
Expand Down Expand Up @@ -311,7 +307,9 @@ class Contact
std::size_t num_q_points() const;

private:
std::shared_ptr<QuadratureRule> _quadrature_rule; // quadrature rule
// quadrature rule
QuadratureRule _quadrature_rule;

std::vector<int> _surfaces; // meshtag values for surfaces

// store index of candidate_surface for each quadrature_surface
Expand Down
2 changes: 0 additions & 2 deletions python/dolfinx_contact/wrappers.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -240,8 +240,6 @@ NB_MODULE(cpp, m)
})
.def("coefficients_size", &dolfinx_contact::Contact::coefficients_size,
nb::arg("meshtie"), nb::arg("V"))
.def("set_quadrature_rule",
&dolfinx_contact::Contact::set_quadrature_rule)
.def("set_search_radius", &dolfinx_contact::Contact::set_search_radius)
.def("generate_kernel",
[](dolfinx_contact::Contact& self, dolfinx_contact::Kernel type,
Expand Down

0 comments on commit 61c0b86

Please sign in to comment.