Skip to content

Commit

Permalink
Code simplifications and tidying (#182)
Browse files Browse the repository at this point in the history 
* Formatting udpates

* Improve readability

* Simplifying Quadrature class.

* Remove quadrature setter

* Simplify constructor

* Improve docs

* Work-around

* Simplify

* Simplify

* Simplify

* Updates

* Simplifications

* Simplify

* Remove commented constructor

* More simplifications

* Improve docs

* Attempt at fixing docs

* Doc update

* Add debug check

* Remove unused code

* Remove check

* Add check

* Division fix

* Reverse order

* Test search

* Type fix

* Test

* Clean up

* Simplify

* Update cpp/Contact.cpp

Co-authored-by: Jørgen Schartum Dokken <[email protected]>

* Simplify

* Update wrapper

---------

Co-authored-by: Jørgen Schartum Dokken <[email protected]>
  • Loading branch information
@garth-wells @jorgensd
garth-wells and jorgensd authored Jun 19, 2024
1 parent 3e5ed8f commit e642d82
Show file tree
Hide file tree
Showing 66 changed files with 1,912 additions and 1,288 deletions.
17 changes: 9 additions & 8 deletions .github/workflows/test-app.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,15 +52,16 @@ jobs:
needs: lint
runs-on: ubuntu-latest
container: ghcr.io/fenics/test-env:current-mpich
strategy:
matrix:
# Complex mode not working
# run_mode: ["Release", "Debug"]
petsc_mode: ["real"]
petsc_precision: ["64"]
petsc_int: ["32"]
# strategy:
# matrix:
# # Complex mode not working
# # run_mode: ["Release", "Debug"]
# petsc_mode: ["real"]
# petsc_precision: ["64"]
# petsc_int: ["32"]
env:
PETSC_ARCH: linux-gnu-${{ matrix.petsc_mode }}${{ matrix.petsc_precision}}-${{ matrix.petsc_int }}
# PETSC_ARCH: linux-gnu-${{ matrix.petsc_mode }}${{ matrix.petsc_precision}}-${{ matrix.petsc_int }}
PETSC_ARCH: linux-gnu-real64-32
# OMPI_ALLOW_RUN_AS_ROOT: 1
# OMPI_ALLOW_RUN_AS_ROOT_CONFIRM: 1
# OMPI_MCA_rmaps_base_oversubscribe: 1
Expand Down
589 changes: 318 additions & 271 deletions cpp/Contact.cpp
View file

Large diffs are not rendered by default.

76 changes: 41 additions & 35 deletions cpp/Contact.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,15 +35,31 @@ namespace dolfinx_contact
class Contact
{
public:
/// empty constructor
Contact() = default;
/// Constructor
///
/// @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
/// @param[in] markers List of meshtags defining the contact surfaces
/// @param[in] surfaces Adjacency list. Links of i contains meshtag values
/// associated with ith meshtag in markers
/// @param[in] contact_pairs list of pairs (i, j) marking the ith and jth
/// surface in surfaces->array() as a contact pair
///
/// @todo Remove this constructor
///
/// @param[in] markers Tags `markers[i]` for surface `i`.
/// @param[in] surfaces For surface `i`, `surfaces.links(i)` gives the
/// tag values of the facets in `markers[i]` that make up surface `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.
Expand All @@ -52,15 +68,11 @@ 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);

/// Return meshtag value for surface with index surface
/// @param[in] surface - the index of the surface
int surface_mt(int surface) const { return _surfaces[surface]; }
const std::vector<ContactMode>& mode, int q_deg = 3);

/// Return contact pair
/// @param[in] pair - the index of the contact pair
const std::array<int, 2>& contact_pair(int pair) const
std::array<int, 2> contact_pair(int pair) const
{
return _contact_pairs[pair];
}
Expand All @@ -70,22 +82,16 @@ class Contact
/// @return TODO
std::span<const std::int32_t> active_entities(int s) const
{
return _cell_facet_pairs->links(s);
return _cell_facet_pairs.links(s);
}

/// Return number of facets in surface s owned by the process
/// @param[in] s TODO
/// @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 @@ -144,7 +150,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 @@ -219,8 +225,8 @@ class Contact
/// _qp_phys, _phi_ref_facets
void create_distance_map(int pair);

/// Compute and pack the gap function for each quadrature point the set of
/// facets.
/// Compute and pack the gap function for each quadrature point the
/// set of facets.
///
/// For a set of facets; go through the quadrature points on each
/// facet find the closest facet on the other surface and compute the
Expand Down Expand Up @@ -252,9 +258,9 @@ class Contact
pack_grad_test_functions(int pair,
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
/// 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
Expand All @@ -266,6 +272,7 @@ class Contact
///
/// Compute function on opposite surface at quadrature points of
/// facets
///
/// @param[in] pair TODO
/// @param[in] u TODO
/// @return TODO
Expand Down Expand Up @@ -311,8 +318,8 @@ class Contact
std::size_t num_q_points() const;

private:
std::shared_ptr<QuadratureRule> _quadrature_rule; // quadrature rule
std::vector<int> _surfaces; // meshtag values for surfaces
// quadrature rule
QuadratureRule _quadrature_rule;

// store index of candidate_surface for each quadrature_surface
std::vector<std::array<int, 2>> _contact_pairs;
Expand All @@ -330,7 +337,7 @@ class Contact
// surface output of compute_distance_map
std::vector<std::vector<double>> _reference_contact_points;

// shape associated with _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
Expand All @@ -344,15 +351,14 @@ class Contact
// maximum number of cells linked to a cell on ith surface
std::vector<std::size_t> _max_links;

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

// 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;

Expand Down
4 changes: 2 additions & 2 deletions cpp/MeshTie.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,8 +24,8 @@ class MeshTie : public Contact
///
/// @param[in] markers List of meshtags defining the connected
/// surfaces
/// @param[in] surfaces Adjacency list. Links of i contains meshtag
/// values associated with ith meshtag in markers
/// @param[in] surfaces Links of i contains meshtag values associated
/// with ith meshtag in markers
/// @param[in] connected_pairs list of pairs (i, j) marking the ith
/// and jth surface in surfaces->array() as a pair of connected
/// surfaces
Expand Down
76 changes: 47 additions & 29 deletions cpp/QuadratureRule.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,17 +16,18 @@
using namespace dolfinx_contact;

//----------------------------------------------------------------------------
QuadratureRule::QuadratureRule(dolfinx::mesh::CellType ct, int degree, int dim,
basix::quadrature::type type)
: _cell_type(ct), _degree(degree), _type(type), _dim(dim)
QuadratureRule::QuadratureRule(dolfinx::mesh::CellType cell, int degree,
int dim, basix::quadrature::type type)
: _cell_type(cell), _degree(degree), _type(type), _dim(dim),
_num_sub_entities(basix::cell::num_sub_entities(
dolfinx::mesh::cell_type_to_basix_type(cell), dim))
{

basix::cell::type b_ct = dolfinx::mesh::cell_type_to_basix_type(ct);
_num_sub_entities = basix::cell::num_sub_entities(b_ct, dim);
_tdim = basix::cell::topological_dimension(b_ct);
assert(dim <= 3);
// If cell dimension no pushing forward
if (_tdim == std::size_t(dim))
basix::cell::type b_ct = dolfinx::mesh::cell_type_to_basix_type(cell);
int tdim = basix::cell::topological_dimension(b_ct);

// If cell dimension no push-forward
if (tdim == dim)
{
std::array<std::vector<double>, 2> quadrature
= basix::quadrature::make_quadrature<double>(
Expand All @@ -37,7 +38,7 @@ QuadratureRule::QuadratureRule(dolfinx::mesh::CellType ct, int degree, int dim,
mdspan_t<const double, 2> qp(quadrature.front().data(), num_points,
pt_shape);

_points = std::vector<double>(num_points * _num_sub_entities * _tdim);
_points = std::vector<double>(num_points * _num_sub_entities * tdim);
_entity_offset = std::vector<std::size_t>(_num_sub_entities + 1, 0);
_weights = std::vector<double>(num_points * _num_sub_entities);
for (std::int32_t i = 0; i < _num_sub_entities; i++)
Expand All @@ -46,9 +47,8 @@ QuadratureRule::QuadratureRule(dolfinx::mesh::CellType ct, int degree, int dim,
for (std::size_t j = 0; j < num_points; ++j)
{
_weights[i * num_points * _num_sub_entities + j] = q_weights[j];
for (std::size_t k = 0; k < _tdim; ++k)
_points[i * num_points * _num_sub_entities + j * _tdim + k]
= qp(j, k);
for (int k = 0; k < tdim; ++k)
_points[i * num_points * _num_sub_entities + j * tdim + k] = qp(j, k);
}
}
}
Expand All @@ -73,17 +73,17 @@ QuadratureRule::QuadratureRule(dolfinx::mesh::CellType ct, int degree, int dim,
type, et, basix::polyset::type::standard, degree);
const std::vector<double>& q_weights = quadrature.back();
const std::vector<double>& q_points = quadrature.front();
const std::size_t num_points = q_weights.size();
const std::size_t tdim = q_points.size() / q_weights.size();
num_points_per_entity[i] = num_points;

const std::array<std::size_t, 4> e_tab_shape
= entity_element.tabulate_shape(0, num_points);
num_points_per_entity[i] = q_weights.size();

std::array<std::size_t, 4> e_tab_shape
= entity_element.tabulate_shape(0, q_weights.size());
std::vector<double> reference_entity_b(std::reduce(
e_tab_shape.cbegin(), e_tab_shape.cend(), 1, std::multiplies{}));

entity_element.tabulate(0, q_points, {num_points, tdim},
reference_entity_b);
entity_element.tabulate(
0, q_points, {q_weights.size(), q_points.size() / q_weights.size()},
reference_entity_b);

mdspan_t<const double, 4> basis_full(reference_entity_b.data(),
e_tab_shape);
Expand All @@ -95,30 +95,37 @@ QuadratureRule::QuadratureRule(dolfinx::mesh::CellType ct, int degree, int dim,
= basix::cell::sub_entity_geometry<double>(b_ct, dim, i);
mdspan_t<const double, 2> coords(sub_geomb.data(), sub_geom_shape);

// Push forward quadrature point from reference entity to reference
// entity on cell
const std::size_t offset = _points.size();
_points.resize(_points.size() + num_points * coords.extent(1));
mdspan_t<double, 2> entity_qp(_points.data() + offset, num_points,
// Push forward quadrature point from reference entity to
// reference entity on cell
std::size_t offset = _points.size();
_points.resize(_points.size() + q_weights.size() * coords.extent(1));
mdspan_t<double, 2> entity_qp(_points.data() + offset, q_weights.size(),
coords.extent(1));
assert(coords.extent(1) == _tdim);

assert(coords.extent(1) == (std::size_t)tdim);
dolfinx::math::dot(phi, coords, entity_qp);
const std::size_t weights_offset = _weights.size();
std::size_t weights_offset = _weights.size();
_weights.resize(_weights.size() + q_weights.size());
std::copy(q_weights.cbegin(), q_weights.cend(),
std::next(_weights.begin(), weights_offset));
}

_entity_offset = std::vector<std::size_t>(_num_sub_entities + 1, 0);
std::partial_sum(num_points_per_entity.begin(), num_points_per_entity.end(),
std::next(_entity_offset.begin()));
}
}
//----------------------------------------------------------------------------
const std::vector<double>& QuadratureRule::points() const { return _points; }
//----------------------------------------------------------------------------
const std::vector<double>& QuadratureRule::weights() const { return _weights; }
//----------------------------------------------------------------------------
int QuadratureRule::dim() const { return _dim; }
//----------------------------------------------------------------------------
dolfinx::mesh::CellType QuadratureRule::cell_type(int i) const
{
basix::cell::type b_ct = dolfinx::mesh::cell_type_to_basix_type(_cell_type);
assert(i < _num_sub_entities);

basix::cell::type et = basix::cell::sub_entity_type(b_ct, _dim, i);
return dolfinx::mesh::cell_type_from_basix_type(et);
}
Expand All @@ -133,10 +140,16 @@ std::size_t QuadratureRule::num_points(int i) const
return _entity_offset[i + 1] - _entity_offset[i];
}
//----------------------------------------------------------------------------
std::size_t QuadratureRule::tdim() const
{
return dolfinx::mesh::cell_dim(_cell_type);
}
//----------------------------------------------------------------------------
mdspan_t<const double, 2> QuadratureRule::points(int i) const
{
assert(i < _num_sub_entities);
mdspan_t<const double, 2> all_points(_points.data(), _weights.size(), _tdim);
mdspan_t<const double, 2> all_points(_points.data(), _weights.size(),
this->tdim());
return MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
all_points, std::pair(_entity_offset[i], _entity_offset[i + 1]),
MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
Expand All @@ -149,3 +162,8 @@ std::span<const double> QuadratureRule::weights(int i) const
_entity_offset[i + 1] - _entity_offset[i]);
}
//----------------------------------------------------------------------------
const std::vector<std::size_t>& QuadratureRule::offset() const
{
return _entity_offset;
}
//----------------------------------------------------------------------------
Loading

0 comments on commit e642d82

Please sign in to comment.