Skip to content

Commit

Permalink
More const correctness
Browse files Browse the repository at this point in the history
  • Loading branch information
@garth-wells
garth-wells committed Jun 8, 2024
1 parent 9235f10 commit 1301fa1
Show file tree
Hide file tree
Showing 6 changed files with 88 additions and 66 deletions.
49 changes: 26 additions & 23 deletions cpp/Contact.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -142,7 +142,7 @@ dolfinx_contact::Contact::Contact(
}
//------------------------------------------------------------------------------------------------
std::pair<std::vector<double>, std::array<std::size_t, 3>>
dolfinx_contact::Contact::qp_phys(int surface)
dolfinx_contact::Contact::qp_phys(int surface) const
{
const std::size_t num_facets = _local_facets[surface];
const std::size_t num_q_points
Expand All @@ -153,7 +153,8 @@ dolfinx_contact::Contact::qp_phys(int surface)
}
//------------------------------------------------------------------------------------------------
std::size_t dolfinx_contact::Contact::coefficients_size(
bool meshtie, std::shared_ptr<const dolfinx::fem::FunctionSpace<double>> V)
bool meshtie,
std::shared_ptr<const dolfinx::fem::FunctionSpace<double>> V) const
{
// mesh data
const std::size_t gdim = _mesh->geometry().dim(); // geometrical dimension
Expand Down Expand Up @@ -208,9 +209,9 @@ std::size_t dolfinx_contact::Contact::coefficients_size(
return std::accumulate(cstrides.cbegin(), cstrides.cend(), 0);
};
}

//------------------------------------------------------------------------------------------------
Mat dolfinx_contact::Contact::create_petsc_matrix(
const dolfinx::fem::Form<PetscScalar>& a, const std::string& type)
const dolfinx::fem::Form<PetscScalar>& a, const std::string& type) const
{

// Build standard sparsity pattern
Expand Down Expand Up @@ -311,7 +312,7 @@ void dolfinx_contact::Contact::create_distance_map(int pair)
}
//------------------------------------------------------------------------------------------------
std::pair<std::vector<PetscScalar>, int>
dolfinx_contact::Contact::pack_nx(int pair)
dolfinx_contact::Contact::pack_nx(int pair) const
{
auto [quadrature_mt, candidate_mt] = _contact_pairs[pair];
const std::shared_ptr<const dolfinx::mesh::Mesh<double>>& quadrature_mesh
Expand Down Expand Up @@ -415,7 +416,8 @@ dolfinx_contact::Contact::pack_nx(int pair)
//------------------------------------------------------------------------------------------------
dolfinx_contact::kernel_fn<PetscScalar>
dolfinx_contact::Contact::generate_kernel(
Kernel type, std::shared_ptr<const dolfinx::fem::FunctionSpace<double>> V)
Kernel type,
std::shared_ptr<const dolfinx::fem::FunctionSpace<double>> V) const
{
std::size_t max_links
= *std::max_element(_max_links.begin(), _max_links.end());
Expand Down Expand Up @@ -475,7 +477,7 @@ void dolfinx_contact::Contact::max_links(int pair)
}
//------------------------------------------------------------------------------------------------
std::pair<std::vector<PetscScalar>, int>
dolfinx_contact::Contact::pack_gap(int pair)
dolfinx_contact::Contact::pack_gap(int pair) const
{
// FIXME: This function should take in the quadrature points
// (push_forward_quadrature) of the relevant facet, and the reference
Expand Down Expand Up @@ -522,8 +524,8 @@ dolfinx_contact::Contact::pack_gap(int pair)
std::array<std::size_t, 2> shape = _reference_contact_shape[pair];

const int q_gdim = quadrature_mesh->geometry().dim();
mdspan3_t qp_span(_qp_phys[quadrature_mt].data(), num_facets, num_q_point,
q_gdim);
cmdspan3_t qp_span(_qp_phys[quadrature_mt].data(), num_facets, num_q_point,
q_gdim);

// Get information about submesh geometry and topology
std::span<const double> x_g = geometry.x();
Expand Down Expand Up @@ -598,7 +600,8 @@ dolfinx_contact::Contact::pack_gap(int pair)
//------------------------------------------------------------------------------------------------
std::pair<std::vector<PetscScalar>, int>
dolfinx_contact::Contact::pack_test_functions(
int pair, std::shared_ptr<const dolfinx::fem::FunctionSpace<double>> V)
int pair,
std::shared_ptr<const dolfinx::fem::FunctionSpace<double>> V) const
{
auto [quadrature_mt, candidate_mt] = _contact_pairs[pair];

Expand Down Expand Up @@ -776,7 +779,7 @@ void dolfinx_contact::Contact::crop_invalid_points(std::size_t pair,
//------------------------------------------------------------------------------------------------
std::pair<std::vector<PetscScalar>, int>
dolfinx_contact::Contact::pack_u_contact(
int pair, std::shared_ptr<dolfinx::fem::Function<PetscScalar>> u)
int pair, std::shared_ptr<dolfinx::fem::Function<PetscScalar>> u) const
{
dolfinx::common::Timer t("Pack contact u");
auto [quadrature_mt, candidate_mt] = _contact_pairs[pair];
Expand Down Expand Up @@ -930,8 +933,8 @@ dolfinx_contact::Contact::pack_gap_plane(int pair, double g)
// FIXME: This does not work for prism meshes
std::size_t num_q_point
= _quadrature_rule->offset()[1] - _quadrature_rule->offset()[0];
mdspan3_t qp_span(_qp_phys[quadrature_mt].data(), num_facets, num_q_point,
gdim);
cmdspan3_t 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);
const int cstride = (int)num_q_point * gdim;
for (std::size_t i = 0; i < num_facets; i++)
Expand All @@ -944,13 +947,12 @@ dolfinx_contact::Contact::pack_gap_plane(int pair, double g)
return {std::move(c), cstride};
}
//------------------------------------------------------------------------------------------------

std::pair<std::vector<PetscScalar>, int>
dolfinx_contact::Contact::pack_ny(int pair)
dolfinx_contact::Contact::pack_ny(int pair) const
{
// FIXME: This function should take in the quadrature points
// (push_forward_quadrature) of the relevant facet, and the reference points
// on the other surface (output of distance map)
// (push_forward_quadrature) of the relevant facet, and the reference
// points on the other surface (output of distance map)
auto [quadrature_mt, candidate_mt] = _contact_pairs[pair];

// Get mesh info for candidate side
Expand Down Expand Up @@ -1329,7 +1331,8 @@ void dolfinx_contact::Contact::assemble_vector(
//-----------------------------------------------------------------------------------------------
std::pair<std::vector<PetscScalar>, int>
dolfinx_contact::Contact::pack_grad_test_functions(
int pair, std::shared_ptr<const dolfinx::fem::FunctionSpace<double>> V)
int pair,
std::shared_ptr<const dolfinx::fem::FunctionSpace<double>> V) const
{
auto [quadrature_mt, candidate_mt] = _contact_pairs[pair];
// Mesh info
Expand All @@ -1352,8 +1355,8 @@ dolfinx_contact::Contact::pack_grad_test_functions(
const std::size_t num_q_points
= _quadrature_rule->offset()[1] - _quadrature_rule->offset()[0];
std::vector<double> q_points(std::size_t(num_q_points) * std::size_t(gdim));
mdspan3_t qp_span(_qp_phys[quadrature_mt].data(), num_facets, num_q_points,
gdim);
cmdspan3_t qp_span(_qp_phys[quadrature_mt].data(), num_facets, num_q_points,
gdim);

std::shared_ptr<const dolfinx::graph::AdjacencyList<int>> facet_map
= _submesh.facet_map();
Expand Down Expand Up @@ -1448,7 +1451,7 @@ dolfinx_contact::Contact::pack_grad_test_functions(
//-----------------------------------------------------------------------------------------------
std::pair<std::vector<PetscScalar>, int>
dolfinx_contact::Contact::pack_grad_u_contact(
int pair, std::shared_ptr<dolfinx::fem::Function<PetscScalar>> u)
int pair, std::shared_ptr<dolfinx::fem::Function<PetscScalar>> u) const
{
auto [quadrature_mt, candidate_mt] = _contact_pairs[pair];

Expand Down Expand Up @@ -1556,9 +1559,9 @@ void dolfinx_contact::Contact::update_submesh_geometry(
{
_submesh.update_geometry(u);
}

//-----------------------------------------------------------------------------------------------
std::size_t dolfinx_contact::Contact::num_q_points() const
{
return _quadrature_rule->offset()[1] - _quadrature_rule->offset()[0];
}
}
//-----------------------------------------------------------------------------------------------
99 changes: 59 additions & 40 deletions cpp/Contact.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@ class Contact
public:
// empty constructor
Contact() = default;

/// Constructor
/// @param[in] markers List of meshtags defining the contact surfaces
/// @param[in] surfaces Adjacency list. Links of i contains meshtag values
Expand Down Expand Up @@ -85,7 +86,7 @@ class Contact
return _quadrature_rule;
}

std::size_t max_links()
std::size_t max_links() const
{
return *std::max_element(_max_links.begin(), _max_links.end());
}
Expand All @@ -97,7 +98,7 @@ class Contact
/// if false
std::size_t coefficients_size(
bool meshtie,
std::shared_ptr<const dolfinx::fem::FunctionSpace<double>> V);
std::shared_ptr<const dolfinx::fem::FunctionSpace<double>> V) const;

/// return distance map (adjacency map mapping quadrature points on surface
/// to closest facet on other surface)
Expand All @@ -108,12 +109,13 @@ class Contact
return _facet_maps[surface];
}

/// Return the quadrature points on physical facet for each facet on surface
/// Return the quadrature points on physical facet for each facet on
/// surface.
/// @param[in] surface The index of the surface (0 or 1).
/// @returns The quadrature points and shape (num_facets, num_q_points, gdim).
/// The points are flattened row major.
std::pair<std::vector<double>, std::array<std::size_t, 3>>
qp_phys(int surface);
qp_phys(int surface) const;

/// Return the submesh of all cells containing facets of the contact surface
const SubMesh& submesh() const { return _submesh; }
Expand All @@ -135,7 +137,7 @@ class Contact
/// available types
/// @returns Mat The PETSc matrix
Mat create_petsc_matrix(const dolfinx::fem::Form<PetscScalar>& a,
const std::string& type);
const std::string& type) const;

/// Assemble matrix over exterior facets (for contact facets)
///
Expand Down Expand Up @@ -185,9 +187,9 @@ class Contact
/// packed at quadrature points. The coefficient `u` is packed at dofs.
/// @note The vector valued coefficents `gap`, `test_fn`, `u`, `u_opposite`
/// has dimension `bs == gdim`.
kernel_fn<PetscScalar>
generate_kernel(Kernel type,
std::shared_ptr<const dolfinx::fem::FunctionSpace<double>> V);
kernel_fn<PetscScalar> generate_kernel(
Kernel type,
std::shared_ptr<const dolfinx::fem::FunctionSpace<double>> V) const;

/// Compute push forward of quadrature points _qp_ref_facet to the
/// physical facet for each facet in _facet_"origin_meshtag" Creates and
Expand Down Expand Up @@ -216,61 +218,67 @@ class Contact
/// @param[out] c - gap packed on facets. c[i*cstride + gdim * k+ j]
/// contains the jth component of the Gap on the ith facet at kth
/// quadrature point
std::pair<std::vector<PetscScalar>, int> pack_gap(int pair);
std::pair<std::vector<PetscScalar>, int> pack_gap(int pair) const;

/// Compute test functions on opposite surface at quadrature points of
/// facets
/// @param[in] pair - index of contact pair
/// @param[in] gap - gap packed on facets per quadrature point
/// @todo Function signature are args docs look wrong
/// @brief Compute test functions on opposite surface at quadrature
/// points of facets
/// @param[in] pair index of contact pair
/// @param[in] gap gap packed on facets per quadrature point
/// @param[out] c - test functions packed on facets.
std::pair<std::vector<PetscScalar>, int> pack_test_functions(
int pair, std::shared_ptr<const dolfinx::fem::FunctionSpace<double>> V);
int pair,
std::shared_ptr<const dolfinx::fem::FunctionSpace<double>> V) const;

/// @todo Function signature are args docs look wrong
/// Compute gradient of test functions on opposite surface (initial
/// configuration) at quadrature points of facets
/// @param[in] pair - index of contact pair
/// @param[in] pair index of contact pair
/// @param[out] c - test functions packed on facets.
std::pair<std::vector<PetscScalar>, int> pack_grad_test_functions(
int pair, std::shared_ptr<const dolfinx::fem::FunctionSpace<double>> V);
int pair,
std::shared_ptr<const dolfinx::fem::FunctionSpace<double>> V) const;

/// Remove points from facet map with a distance larger than tol
/// in the surface or if the angle of distance vector and opposite surface is
/// too large
/// @param[in] pair - index of the contact pair
/// @param[in] gap - for computing distance
/// @param[in] n_y - normals for checking angle
/// @param[in] tol - max distance
/// @param[in] pair index of the contact pair
/// @param[in] gap for computing distance
/// @param[in] n_y normals for checking angle
/// @param[in] tol max distance
void crop_invalid_points(std::size_t pair, std::span<const double> gap,
std::span<const double> n_y, double tol);

/// @todo Function signature are args docs look wrong
/// Compute function on opposite surface at quadrature points of
/// facets
/// @param[in] pair - index of contact pair
/// @param[in] - gap packed on facets per quadrature point
/// @param[in] pair index of contact pair
/// @param[in] gap packed on facets per quadrature point
/// @param[out] c - test functions packed on facets.
std::pair<std::vector<PetscScalar>, int>
pack_u_contact(int pair,
std::shared_ptr<dolfinx::fem::Function<PetscScalar>> u);
std::shared_ptr<dolfinx::fem::Function<PetscScalar>> u) const;

/// Compute gradient of function on opposite surface at quadrature points of
/// facets
/// @param[in] pair - index of contact pair
/// @param[out] c - test functions packed on facets.
std::pair<std::vector<PetscScalar>, int>
pack_grad_u_contact(int pair,
std::shared_ptr<dolfinx::fem::Function<PetscScalar>> u);
/// @todo Function signature are args docs look wrong
/// Compute gradient of function on opposite surface at quadrature
/// points of facets
/// @param[in] pair index of contact pair
/// @param[out] u test functions packed on facets.
std::pair<std::vector<PetscScalar>, int> pack_grad_u_contact(
int pair, std::shared_ptr<dolfinx::fem::Function<PetscScalar>> u) const;

/// Compute outward surface normal at x
/// @param[in] pair - index of contact pair
/// @returns c - (normals, cstride) ny packed on facets.
std::pair<std::vector<PetscScalar>, int> pack_nx(int pair);
/// @param[in] pair index of contact pair
/// @returns (normals, cstride) ny packed on facets.
std::pair<std::vector<PetscScalar>, int> pack_nx(int pair) const;

/// Compute inward surface normal at Pi(x)
/// @param[in] pair - index of contact pair
/// @returns c - normals ny packed on facets.
std::pair<std::vector<PetscScalar>, int> pack_ny(int pair);
/// @param[in] pair index of contact pair
/// @returns c normals ny packed on facets.
std::pair<std::vector<PetscScalar>, int> pack_ny(int pair) const;

/// Pack gap with rigid surface defined by x[gdim-1] = -g.
/// @todo Make functon cont
/// @brief Pack gap with rigid surface defined by x[gdim-1] = -g.
/// g_vec = zeros(gdim), g_vec[gdim-1] = -g
/// Gap = x - g_vec
/// @param[in] pair - index of contact pair
Expand All @@ -282,7 +290,7 @@ class Contact

/// This function updates the submesh geometry for all submeshes using
/// a function given on the parent mesh
/// @param[in] u - displacement
/// @param[in] u displacement
void update_submesh_geometry(dolfinx::fem::Function<PetscScalar>& u);

/// Return number of quadrature points per facet
Expand All @@ -292,40 +300,51 @@ class Contact
private:
std::shared_ptr<QuadratureRule> _quadrature_rule; // quadrature rule
std::vector<int> _surfaces; // meshtag values for surfaces

// store index of candidate_surface for each quadrature_surface
std::vector<std::array<int, 2>> _contact_pairs;
std::shared_ptr<dolfinx::mesh::Mesh<double>> _mesh; // mesh

// _facets_maps[i] = adjacency list of closest facet on candidate surface
// for every quadrature point in _qp_phys[i] (quadrature points on every
// facet of ith surface)
std::vector<
std::shared_ptr<const dolfinx::graph::AdjacencyList<std::int32_t>>>
_facet_maps;

// reference points of the contact points on the opposite surface for each
// surface output of compute_distance_map
std::vector<std::vector<double>> _reference_contact_points;

// shape associated with _reference_contact_points
std::vector<std::array<std::size_t, 2>> _reference_contact_shape;
// _qp_phys[i] contains the quadrature points on the physical facets for
// each facet on ith surface in _surfaces

// _qp_phys[i] contains the quadrature points on the physical facets
// for each facet on ith surface in _surfaces
std::vector<std::vector<double>> _qp_phys;

// quadrature points on facets of reference cell
std::vector<double> _reference_basis;
std::array<std::size_t, 4> _reference_shape;

// maximum number of cells linked to a cell on ith surface
std::vector<std::size_t> _max_links;

// submesh containing all cells linked to facets on any of the contact
// surfaces
SubMesh _submesh;

// Adjacency list linking facets as (cell, facet) pairs to the index of the
// surface. The pairs are flattened row-major
std::shared_ptr<const dolfinx::graph::AdjacencyList<std::int32_t>>
_cell_facet_pairs;

// number of facets owned by process for each surface
std::vector<std::size_t> _local_facets;

// Contact search mode
std::vector<ContactMode> _mode;

// Search radius for ray-tracing
double _radius = -1;
};
Expand Down
1 change: 1 addition & 0 deletions cpp/KernelData.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -131,3 +131,4 @@ dolfinx_contact::cmdspan3_t dolfinx_contact::KernelData::ref_jacobians() const
{
return dolfinx_contact::cmdspan3_t(_ref_jacobians.data(), _jac_shape);
}
//-----------------------------------------------------------------------------
2 changes: 1 addition & 1 deletion cpp/SubMesh.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -154,7 +154,7 @@ dolfinx_contact::SubMesh::create_functionspace(
//-----------------------------------------------------------------------------------------------
void dolfinx_contact::SubMesh::copy_function(
dolfinx::fem::Function<PetscScalar>& u_parent,
dolfinx::fem::Function<PetscScalar>& u_sub)
dolfinx::fem::Function<PetscScalar>& u_sub) const
{
// retrieve function space on submesh
std::shared_ptr<const dolfinx::fem::FunctionSpace<double>> V_sub
Expand Down
Loading

0 comments on commit 1301fa1

Please sign in to comment.