Skip to content

Commit

Permalink
Replace typedefs
Browse files Browse the repository at this point in the history
  • Loading branch information
garth-wells committed Jun 9, 2024
1 parent ac69084 commit 1793cc3
Show file tree
Hide file tree
Showing 18 changed files with 148 additions and 117 deletions.
13 changes: 8 additions & 5 deletions cpp/Contact.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -361,11 +361,12 @@ dolfinx_contact::Contact::pack_nx(int pair) const
= dolfinx::mesh::cell_type_to_basix_type(topology->cell_type());
auto [facet_normalsb, n_shape]
= basix::cell::facet_outward_normals<double>(cell_type);
cmdspan2_t facet_normals(facet_normalsb.data(), n_shape);
mdspan_t<const double, 2> facet_normals(facet_normalsb.data(), n_shape);

// Working memory for loop
std::vector<double> coordinate_dofsb(num_dofs_g * gdim);
cmdspan2_t coordinate_dofs(coordinate_dofsb.data(), num_dofs_g, gdim);
mdspan_t<const double, 2> coordinate_dofs(coordinate_dofsb.data(), num_dofs_g,
gdim);
std::array<double, 9> Jb;
std::array<double, 9> Kb;
mdspan2_t J(Jb.data(), gdim, tdim);
Expand Down Expand Up @@ -532,7 +533,8 @@ dolfinx_contact::Contact::pack_gap(int pair) const
const int tdim = topology->dim();

std::vector<double> coordinate_dofsb(num_dofs_g * gdim);
cmdspan2_t coordinate_dofs(coordinate_dofsb.data(), num_dofs_g, gdim);
mdspan_t<const double, 2> coordinate_dofs(coordinate_dofsb.data(), num_dofs_g,
gdim);
std::array<double, 3> coordb;
mdspan2_t coord(coordb.data(), 1, gdim);

Expand Down Expand Up @@ -1021,11 +1023,12 @@ dolfinx_contact::Contact::pack_ny(int pair) const
candidate_mesh->topology()->cell_type());
auto [facet_normalsb, n_shape]
= basix::cell::facet_outward_normals<double>(cell_type);
cmdspan2_t facet_normals(facet_normalsb.data(), n_shape);
mdspan_t<const double, 2> facet_normals(facet_normalsb.data(), n_shape);

// Working memory for loop
std::vector<double> coordinate_dofsb(num_dofs_g * gdim);
cmdspan2_t coordinate_dofs(coordinate_dofsb.data(), num_dofs_g, gdim);
mdspan_t<const double, 2> coordinate_dofs(coordinate_dofsb.data(), num_dofs_g,
gdim);
std::array<double, 9> Jb;
std::array<double, 9> Kb;
mdspan2_t J(Jb.data(), gdim, tdim);
Expand Down
12 changes: 7 additions & 5 deletions cpp/KernelData.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -92,7 +92,8 @@ dolfinx_contact::KernelData::KernelData(
double dolfinx_contact::KernelData::compute_first_facet_jacobian(
const std::size_t facet_index, dolfinx_contact::mdspan2_t J,
dolfinx_contact::mdspan2_t K, dolfinx_contact::mdspan2_t J_tot,
std::span<double> detJ_scratch, dolfinx_contact::cmdspan2_t coords) const
std::span<double> detJ_scratch,
dolfinx_contact::mdspan_t<const double, 2> coords) const
{
dolfinx_contact::cmdspan4_t full_basis(_c_basis_values.data(),
_c_basis_shape);
Expand All @@ -112,12 +113,13 @@ double dolfinx_contact::KernelData::compute_first_facet_jacobian(
}
//-----------------------------------------------------------------------------
void dolfinx_contact::KernelData::update_normal(
std::span<double> n, dolfinx_contact::cmdspan2_t K,
std::span<double> n, dolfinx_contact::mdspan_t<const double, 2> K,
const std::size_t local_index) const
{
return _update_normal(
n, K, dolfinx_contact::cmdspan2_t(_facet_normals.data(), _normals_shape),
local_index);
return _update_normal(n, K,
dolfinx_contact::mdspan_t<const double, 2>(
_facet_normals.data(), _normals_shape),
local_index);
}
//-----------------------------------------------------------------------------
std::span<const double>
Expand Down
21 changes: 11 additions & 10 deletions cpp/KernelData.h
Original file line number Diff line number Diff line change
Expand Up @@ -20,12 +20,13 @@

namespace dolfinx_contact
{
using jac_fn = std::function<double(double, mdspan2_t, mdspan2_t, mdspan2_t,
std::span<double>, cmdspan2_t, s_cmdspan2_t,
cmdspan2_t)>;
using jac_fn = std::function<double(
double, mdspan2_t, mdspan2_t, mdspan2_t, std::span<double>,
mdspan_t<const double, 2>, s_cmdspan2_t, mdspan_t<const double, 2>)>;

using normal_fn = std::function<void(std::span<double>, cmdspan2_t, cmdspan2_t,
const std::size_t)>;
using normal_fn
= std::function<void(std::span<double>, mdspan_t<const double, 2>,
mdspan_t<const double, 2>, const std::size_t)>;

class KernelData
{
Expand Down Expand Up @@ -100,9 +101,9 @@ class KernelData
const std::vector<std::size_t>& offsets_array() const { return _offsets; }

// Return reference facet normals
cmdspan2_t facet_normals() const
mdspan_t<const double, 2> facet_normals() const
{
return cmdspan2_t(_facet_normals.data(), _normals_shape);
return mdspan_t<const double, 2>(_facet_normals.data(), _normals_shape);
}

/// Compute the following jacobians on a given facet:
Expand All @@ -119,7 +120,7 @@ class KernelData
double update_jacobian(std::size_t q, std::size_t facet_index, double detJ,
mdspan2_t J, mdspan2_t K, mdspan2_t J_tot,
std::span<double> detJ_scratch,
cmdspan2_t coords) const
mdspan_t<const double, 2> coords) const
{
cmdspan4_t full_basis(_c_basis_values.data(), _c_basis_shape);
const std::size_t q_pos = _qp_offsets[facet_index] + q;
Expand Down Expand Up @@ -148,13 +149,13 @@ class KernelData
double compute_first_facet_jacobian(std::size_t facet_index, mdspan2_t J,
mdspan2_t K, mdspan2_t J_tot,
std::span<double> detJ_scratch,
cmdspan2_t coords) const;
mdspan_t<const double, 2> coords) const;

/// update normal
/// @param[in, out] n The facet normal
/// @param[in] K The inverse Jacobian
/// @param[in] local_index The facet index local to the cell
void update_normal(std::span<double> n, cmdspan2_t K,
void update_normal(std::span<double> n, mdspan_t<const double, 2> K,
std::size_t local_index) const;

/// Return quadrature weights for the i-th facet
Expand Down
9 changes: 5 additions & 4 deletions cpp/QuadratureRule.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,8 @@ dolfinx_contact::QuadratureRule::QuadratureRule(dolfinx::mesh::CellType ct,
std::vector<double>& q_weights = quadrature.back();
std::size_t num_points = q_weights.size();
std::size_t pt_shape = quadrature.front().size() / num_points;
cmdspan2_t qp(quadrature.front().data(), num_points, pt_shape);
mdspan_t<const double, 2> qp(quadrature.front().data(), num_points,
pt_shape);

_points = std::vector<double>(num_points * _num_sub_entities * _tdim);
_entity_offset = std::vector<std::size_t>(_num_sub_entities + 1, 0);
Expand Down Expand Up @@ -91,7 +92,7 @@ dolfinx_contact::QuadratureRule::QuadratureRule(dolfinx::mesh::CellType ct,

auto [sub_geomb, sub_geom_shape]
= basix::cell::sub_entity_geometry<double>(b_ct, dim, i);
cmdspan2_t coords(sub_geomb.data(), sub_geom_shape);
mdspan_t<const double, 2> coords(sub_geomb.data(), sub_geom_shape);

// Push forward quadrature point from reference entity to reference
// entity on cell
Expand Down Expand Up @@ -132,10 +133,10 @@ std::size_t QuadratureRule::num_points(int i) const
return _entity_offset[i + 1] - _entity_offset[i];
}
//-----------------------------------------------------------------------------------------------
cmdspan2_t QuadratureRule::points(int i) const
mdspan_t<const double, 2> QuadratureRule::points(int i) const
{
assert(i < _num_sub_entities);
cmdspan2_t all_points(_points.data(), _weights.size(), _tdim);
mdspan_t<const double, 2> all_points(_points.data(), _weights.size(), _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 Down
7 changes: 1 addition & 6 deletions cpp/QuadratureRule.h
Original file line number Diff line number Diff line change
Expand Up @@ -21,8 +21,6 @@ using mdspan_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<

using mdspan2_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
double, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
using cmdspan2_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
const double, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 2>>;
using mdspan3_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
double, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 3>>;
using mdspan4_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
Expand All @@ -35,9 +33,6 @@ using s_cmdspan2_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
using s_cmdspan3_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
const double, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 3>,
stdex::layout_stride>;
using s_cmdspan4_t = MDSPAN_IMPL_STANDARD_NAMESPACE::mdspan<
const double, MDSPAN_IMPL_STANDARD_NAMESPACE::dextents<std::size_t, 4>,
stdex::layout_stride>;

class QuadratureRule
{
Expand Down Expand Up @@ -85,7 +80,7 @@ class QuadratureRule

/// Return the quadrature points for the ith entity
/// @param[in] i The local entity index
cmdspan2_t points(int i) const;
mdspan_t<const double, 2> points(int i) const;

/// Return the quadrature weights for the ith entity
/// @param[in] i The local entity index
Expand Down
4 changes: 2 additions & 2 deletions cpp/RayTracing.h
Original file line number Diff line number Diff line change
Expand Up @@ -273,7 +273,7 @@ int raytracing_cell(
auto dphi = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
basis, std::pair{1, tdim + 1}, 0,
MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent, 0);
cmdspan2_t coords(coordinate_dofs.data(), cmap.dim(), gdim);
mdspan_t<const double, 2> coords(coordinate_dofs.data(), cmap.dim(), gdim);
mdspan2_t _xk(x_k.data(), 1, gdim);
for (int k = 0; k < max_iter; ++k)
{
Expand Down Expand Up @@ -500,7 +500,7 @@ compute_ray(const dolfinx::mesh::Mesh<double>& mesh,
std::span<const double, tdim - 1> xi, std::span<double, tdim> X)
{
const std::vector<int>& facet = facets[facet_index];
dolfinx_contact::cmdspan2_t x(xb.data(), x_shape);
dolfinx_contact::mdspan_t<const double, 2> x(xb.data(), x_shape);
for (std::size_t i = 0; i < tdim; ++i)
{
X[i] = x(facet.front(), i);
Expand Down
12 changes: 7 additions & 5 deletions cpp/coefficients.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -15,9 +15,9 @@ using namespace dolfinx_contact;
void dolfinx_contact::transformed_push_forward(
const dolfinx::fem::FiniteElement<double>* element,
cmdspan4_t reference_basis, std::vector<double>& element_basisb,
mdspan3_t basis_values, cmdspan2_t J, cmdspan2_t K, double detJ,
std::size_t basis_offset, std::size_t q, std::int32_t cell,
std::span<const std::uint32_t> cell_info)
mdspan3_t basis_values, mdspan_t<const double, 2> J,
mdspan_t<const double, 2> K, double detJ, std::size_t basis_offset,
std::size_t q, std::int32_t cell, std::span<const std::uint32_t> cell_info)
{
const std::function<void(const std::span<PetscScalar>&,
const std::span<const std::uint32_t>&, std::int32_t,
Expand All @@ -28,8 +28,10 @@ void dolfinx_contact::transformed_push_forward(
// Get push forward function
auto push_forward_fn
= element->basix_element()
.map_fn<dolfinx_contact::mdspan2_t, dolfinx_contact::cmdspan2_t,
dolfinx_contact::cmdspan2_t, dolfinx_contact::cmdspan2_t>();
.map_fn<dolfinx_contact::mdspan2_t,
dolfinx_contact::mdspan_t<const double, 2>,
dolfinx_contact::mdspan_t<const double, 2>,
dolfinx_contact::mdspan_t<const double, 2>>();
mdspan2_t element_basis(element_basisb.data(), basis_values.extent(1),
basis_values.extent(2));

Expand Down
6 changes: 3 additions & 3 deletions cpp/coefficients.h
Original file line number Diff line number Diff line change
Expand Up @@ -30,9 +30,9 @@ namespace dolfinx_contact
void transformed_push_forward(
const dolfinx::fem::FiniteElement<double>* element,
cmdspan4_t reference_basis, std::vector<double>& element_basisb,
mdspan3_t basis_values, cmdspan2_t J, cmdspan2_t K, double detJ,
std::size_t basis_offset, std::size_t q, std::int32_t cell,
std::span<const std::uint32_t> cell_info);
mdspan3_t basis_values, mdspan_t<const double, 2> J,
mdspan_t<const double, 2> K, double detJ, std::size_t basis_offset,
std::size_t q, std::int32_t cell, std::span<const std::uint32_t> cell_info);

/// @brief Pack a coefficient at quadrature points.
///
Expand Down
18 changes: 12 additions & 6 deletions cpp/contact_kernels.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -72,7 +72,8 @@ dolfinx_contact::generate_contact_kernel(
const std::uint32_t tdim = kd.tdim();

// NOTE: DOLFINx has 3D input coordinate dofs
cmdspan2_t coord(coordinate_dofs, kd.num_coordinate_dofs(), 3);
mdspan_t<const double, 2> coord(coordinate_dofs, kd.num_coordinate_dofs(),
3);

// Create data structures for jacobians
// We allocate more memory than required, but its better for the compiler
Expand Down Expand Up @@ -221,7 +222,8 @@ dolfinx_contact::generate_contact_kernel(
const std::uint32_t tdim = kd.tdim();

// NOTE: DOLFINx has 3D input coordinate dofs
cmdspan2_t coord(coordinate_dofs, kd.num_coordinate_dofs(), 3);
mdspan_t<const double, 2> coord(coordinate_dofs, kd.num_coordinate_dofs(),
3);

// Create data structures for jacobians
// We allocate more memory than required, but its better for the compiler
Expand Down Expand Up @@ -387,7 +389,8 @@ dolfinx_contact::generate_contact_kernel(
const std::uint32_t tdim = kd.tdim();

// NOTE: DOLFINx has 3D input coordinate dofs
cmdspan2_t coord(coordinate_dofs, kd.num_coordinate_dofs(), 3);
mdspan_t<const double, 2> coord(coordinate_dofs, kd.num_coordinate_dofs(),
3);

// Create data structures for jacobians
// We allocate more memory than required, but its better for the compiler
Expand Down Expand Up @@ -571,7 +574,8 @@ dolfinx_contact::generate_contact_kernel(
const std::uint32_t tdim = kd.tdim();

// NOTE: DOLFINx has 3D input coordinate dofs
cmdspan2_t coord(coordinate_dofs, kd.num_coordinate_dofs(), 3);
mdspan_t<const double, 2> coord(coordinate_dofs, kd.num_coordinate_dofs(),
3);

// Create data structures for jacobians
// We allocate more memory than required, but its better for the compiler
Expand Down Expand Up @@ -802,7 +806,8 @@ dolfinx_contact::generate_contact_kernel(
const std::uint32_t tdim = kd.tdim();

// NOTE: DOLFINx has 3D input coordinate dofs
cmdspan2_t coord(coordinate_dofs, kd.num_coordinate_dofs(), 3);
mdspan_t<const double, 2> coord(coordinate_dofs, kd.num_coordinate_dofs(),
3);

// Create data structures for jacobians
// We allocate more memory than required, but its better for the compiler
Expand Down Expand Up @@ -999,7 +1004,8 @@ dolfinx_contact::generate_contact_kernel(
const std::uint32_t tdim = kd.tdim();

// NOTE: DOLFINx has 3D input coordinate dofs
cmdspan2_t coord(coordinate_dofs, kd.num_coordinate_dofs(), 3);
mdspan_t<const double, 2> coord(coordinate_dofs, kd.num_coordinate_dofs(),
3);

// Create data structures for jacobians
// We allocate more memory than required, but its better for the compiler
Expand Down
10 changes: 4 additions & 6 deletions cpp/elasticity.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@

//-----------------------------------------------------------------------------
void dolfinx_contact::compute_normal_strain_basis(
mdspan2_t epsn, mdspan2_t tr, dolfinx_contact::cmdspan2_t K,
mdspan2_t epsn, mdspan2_t tr, dolfinx_contact::mdspan_t<const double, 2> K,
dolfinx_contact::s_cmdspan3_t dphi, std::array<double, 3> n_surf,
std::span<const double> n_phys, std::size_t q_pos)
{
Expand Down Expand Up @@ -42,11 +42,9 @@ void dolfinx_contact::compute_normal_strain_basis(
}
}
//-----------------------------------------------------------------------------
void dolfinx_contact::compute_sigma_n_basis(mdspan3_t sig_n, cmdspan2_t K,
s_cmdspan3_t dphi,
std::span<const double> n,
double mu, double lmbda,
std::size_t q_pos)
void dolfinx_contact::compute_sigma_n_basis(
mdspan3_t sig_n, mdspan_t<const double, 2> K, s_cmdspan3_t dphi,
std::span<const double> n, double mu, double lmbda, std::size_t q_pos)
{
const std::size_t ndofs_cell = sig_n.extent(0);
const std::size_t bs = K.extent(1);
Expand Down
9 changes: 5 additions & 4 deletions cpp/elasticity.h
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,8 @@ namespace dolfinx_contact
/// @param[in] n_1 1st normal vector, typically n_surf
/// @param[in] n_2 2nd normal vector, typically n_phys
/// @param[in] q_pos offset of quadrature point for accessing dphi
void compute_normal_strain_basis(mdspan2_t epsn, mdspan2_t tr, cmdspan2_t K,
void compute_normal_strain_basis(mdspan2_t epsn, mdspan2_t tr,
mdspan_t<const double, 2> K,
s_cmdspan3_t s_dphi, std::array<double, 3> n_1,
std::span<const double> n_2,
std::size_t q_pos);
Expand All @@ -46,9 +47,9 @@ void compute_normal_strain_basis(mdspan2_t epsn, mdspan2_t tr, cmdspan2_t K,
/// @param[in] mu The poisson ratio
/// @param[in] lmbda The 1st Lame parameter
/// @param[in] q_pos The offset of quadrature point for accessing dphi
void compute_sigma_n_basis(mdspan3_t sig_n, cmdspan2_t K, s_cmdspan3_t dphi,
std::span<const double> n, double mu, double lmbda,
std::size_t q_pos);
void compute_sigma_n_basis(mdspan3_t sig_n, mdspan_t<const double, 2> K,
s_cmdspan3_t dphi, std::span<const double> n,
double mu, double lmbda, std::size_t q_pos);

/// @brief Compute sigma(u)*n from grad(u)
///
Expand Down
10 changes: 6 additions & 4 deletions cpp/geometric_quantities.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ void dolfinx_contact::pull_back_nonaffine(
std::span<double> X, std::span<double> work_array,
std::span<const double> x,
const dolfinx::fem::CoordinateElement<double>& cmap,
cmdspan2_t cell_geometry, double tol, const int max_it)
mdspan_t<const double, 2> cell_geometry, double tol, const int max_it)
{
assert((std::size_t)cmap.dim() == cell_geometry.extent(0));
// Temporary data structures for Newton iteration
Expand Down Expand Up @@ -95,9 +95,10 @@ void dolfinx_contact::pull_back_nonaffine(

std::array<double, 3> dolfinx_contact::push_forward_facet_normal(
std::span<double> work_array, std::span<const double> x, std::size_t gdim,
std::size_t tdim, cmdspan2_t coordinate_dofs, const std::size_t facet_index,
std::size_t tdim, mdspan_t<const double, 2> coordinate_dofs,
const std::size_t facet_index,
const dolfinx::fem::CoordinateElement<double>& cmap,
cmdspan2_t reference_normals)
mdspan_t<const double, 2> reference_normals)
{
const std::array<std::size_t, 4> c_shape = cmap.tabulate_shape(1, 1);
const std::size_t basis_size
Expand Down Expand Up @@ -156,7 +157,8 @@ std::array<double, 3> dolfinx_contact::push_forward_facet_normal(
//-----------------------------------------------------------------------------
double
dolfinx_contact::compute_circumradius(const dolfinx::mesh::Mesh<double>& mesh,
double detJ, cmdspan2_t coordinate_dofs)
double detJ,
mdspan_t<const double, 2> coordinate_dofs)
{
const dolfinx::mesh::CellType cell_type = mesh.topology()->cell_type();
const int gdim = mesh.geometry().dim();
Expand Down
Loading

0 comments on commit 1793cc3

Please sign in to comment.