Skip to content

Commit

Permalink
Remove cell number inputs from permutation information (FEniCS#800)
Browse files Browse the repository at this point in the history 
* rename orientation -> entity_local_index in test

* Use full arrays, not rows or elements

* gitignore the tempdirs made by the tests

* update documentation

* Resize array for 1d elements

* Added wrappers to reflection and rotation information

* consts

* Use types, not autos

* More types

* Added RowMajor everywhere

* More RowMajors

* clangformat

* RowMajors in Topology.h

* Use ColumnMajor (default)

* Use col not row

Co-authored-by: Garth N. Wells <gnw20@cam.ac.uk>
  • Loading branch information
@mscroggs @garth-wells
mscroggs and garth-wells authored Feb 26, 2020
1 parent 6bc3ce2 commit 3232b87
Show file tree
Hide file tree
Showing 12 changed files with 301 additions and 201 deletions.
1 change: 1 addition & 0 deletions .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -127,3 +127,4 @@ dist/
*.egg-info
.pytest_cache
.*.swp
*_tempdir_*
96 changes: 57 additions & 39 deletions cpp/dolfinx/fem/assemble_matrix_impl.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -132,6 +132,14 @@ void fem::impl::assemble_cells(
Eigen::Matrix<PetscScalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
Ae;

const Eigen::Ref<const Eigen::Array<bool, Eigen::Dynamic, Eigen::Dynamic>>
cell_edge_reflections = mesh.topology().get_edge_reflections();
const Eigen::Ref<const Eigen::Array<bool, Eigen::Dynamic, Eigen::Dynamic>>
cell_face_reflections = mesh.topology().get_face_reflections();
const Eigen::Ref<
const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>>
cell_face_rotations = mesh.topology().get_face_rotations();

// Iterate over active cells
PetscErrorCode ierr;
for (std::int32_t cell_index : active_cells)
Expand All @@ -141,22 +149,18 @@ void fem::impl::assemble_cells(
for (int j = 0; j < gdim; ++j)
coordinate_dofs(i, j) = x_g(cell_g[pos_g[cell_index] + i], j);

const Eigen::Ref<const Eigen::Array<bool, 1, Eigen::Dynamic>>
cell_edge_reflections
= mesh.topology().get_edge_reflections(cell_index);
const Eigen::Ref<const Eigen::Array<bool, 1, Eigen::Dynamic>>
cell_face_reflections
= mesh.topology().get_face_reflections(cell_index);
const Eigen::Ref<const Eigen::Array<std::uint8_t, 1, Eigen::Dynamic>>
cell_face_rotations = mesh.topology().get_face_rotations(cell_index);

// Tabulate tensor
auto coeff_cell = coeffs.row(cell_index);
const Eigen::Ref<const Eigen::Array<bool, Eigen::Dynamic, 1>> e_ref_cell
= cell_edge_reflections.col(cell_index);
const Eigen::Ref<const Eigen::Array<bool, Eigen::Dynamic, 1>> f_ref_cell
= cell_face_reflections.col(cell_index);
const Eigen::Ref<const Eigen::Array<uint8_t, Eigen::Dynamic, 1>> f_rot_cell
= cell_face_rotations.col(cell_index);
Ae.setZero(num_dofs_per_cell0, num_dofs_per_cell1);
kernel(Ae.data(), coeff_cell.data(), constant_values.data(),
coordinate_dofs.data(), nullptr, nullptr,
cell_edge_reflections.data(), cell_face_reflections.data(),
cell_face_rotations.data());
coordinate_dofs.data(), nullptr, nullptr, e_ref_cell.data(),
f_ref_cell.data(), f_rot_cell.data());

// Zero rows/columns for essential bcs
if (!bc0.empty())
Expand Down Expand Up @@ -226,6 +230,18 @@ void fem::impl::assemble_exterior_facets(
Eigen::Matrix<PetscScalar, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
Ae;

const Eigen::Ref<
const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>>
perms = mesh.topology().get_facet_permutations();

const Eigen::Ref<const Eigen::Array<bool, Eigen::Dynamic, Eigen::Dynamic>>
cell_edge_reflections = mesh.topology().get_edge_reflections();
const Eigen::Ref<const Eigen::Array<bool, Eigen::Dynamic, Eigen::Dynamic>>
cell_face_reflections = mesh.topology().get_face_reflections();
const Eigen::Ref<
const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>>
cell_face_rotations = mesh.topology().get_face_rotations();

// Iterate over all facets
PetscErrorCode ierr;
auto f_to_c = mesh.topology().connectivity(tdim - 1, tdim);
Expand All @@ -248,28 +264,23 @@ void fem::impl::assemble_exterior_facets(
for (int j = 0; j < gdim; ++j)
coordinate_dofs(i, j) = x_g(cell_g[pos_g[cells[0]] + i], j);

// Get the permutation of the facet
const std::uint8_t perm = mesh.topology().get_facet_permutation(
cells[0], tdim - 1, local_facet);

// Get dof maps for cell
auto dmap0 = dofmap0.cell_dofs(cells[0]);
auto dmap1 = dofmap1.cell_dofs(cells[0]);

const Eigen::Ref<const Eigen::Array<bool, 1, Eigen::Dynamic>>
cell_edge_reflections = mesh.topology().get_edge_reflections(cells[0]);
const Eigen::Ref<const Eigen::Array<bool, 1, Eigen::Dynamic>>
cell_face_reflections = mesh.topology().get_face_reflections(cells[0]);
const Eigen::Ref<const Eigen::Array<std::uint8_t, 1, Eigen::Dynamic>>
cell_face_rotations = mesh.topology().get_face_rotations(cells[0]);

// Tabulate tensor
auto coeff_cell = coeffs.row(cells[0]);
const Eigen::Ref<const Eigen::Array<bool, Eigen::Dynamic, 1>> e_ref_cell
= cell_edge_reflections.col(cells[0]);
const Eigen::Ref<const Eigen::Array<bool, Eigen::Dynamic, 1>> f_ref_cell
= cell_face_reflections.col(cells[0]);
const Eigen::Ref<const Eigen::Array<uint8_t, Eigen::Dynamic, 1>> f_rot_cell
= cell_face_rotations.col(cells[0]);
const std::uint8_t perm = perms(local_facet, cells[0]);
Ae.setZero(dmap0.size(), dmap1.size());
fn(Ae.data(), coeff_cell.data(), constant_values.data(),
coordinate_dofs.data(), &local_facet, &perm,
cell_edge_reflections.data(), cell_face_reflections.data(),
cell_face_rotations.data());
coordinate_dofs.data(), &local_facet, &perm, e_ref_cell.data(),
f_ref_cell.data(), f_rot_cell.data());

// Zero rows/columns for essential bcs
if (!bc0.empty())
Expand Down Expand Up @@ -341,6 +352,18 @@ void fem::impl::assemble_interior_facets(
// Temporaries for joint dofmaps
Eigen::Array<PetscInt, Eigen::Dynamic, 1> dmapjoint0, dmapjoint1;

const Eigen::Ref<
const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>>
perms = mesh.topology().get_facet_permutations();

const Eigen::Ref<const Eigen::Array<bool, Eigen::Dynamic, Eigen::Dynamic>>
cell_edge_reflections = mesh.topology().get_edge_reflections();
const Eigen::Ref<const Eigen::Array<bool, Eigen::Dynamic, Eigen::Dynamic>>
cell_face_reflections = mesh.topology().get_face_reflections();
const Eigen::Ref<
const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>>
cell_face_rotations = mesh.topology().get_face_rotations();

// Iterate over all facets
PetscErrorCode ierr;
auto c = mesh.topology().connectivity(tdim - 1, tdim);
Expand Down Expand Up @@ -370,10 +393,7 @@ void fem::impl::assemble_interior_facets(
const std::array<int, 2> local_facet = {local_facet0, local_facet1};

const std::array<std::uint8_t, 2> perm
= {mesh.topology().get_facet_permutation(cells[0], tdim - 1,
local_facet[0]),
mesh.topology().get_facet_permutation(cells[1], tdim - 1,
local_facet[1])};
= {perms(local_facet[0], cells[0]), perms(local_facet[1], cells[1])};

for (int i = 0; i < num_dofs_g; ++i)
{
Expand Down Expand Up @@ -424,19 +444,17 @@ void fem::impl::assemble_interior_facets(
= coeff_cell1.segment(offsets[i], num_entries);
}

const Eigen::Ref<const Eigen::Array<bool, 1, Eigen::Dynamic>>
cell_edge_reflections = mesh.topology().get_edge_reflections(cells[0]);
const Eigen::Ref<const Eigen::Array<bool, 1, Eigen::Dynamic>>
cell_face_reflections = mesh.topology().get_face_reflections(cells[0]);
const Eigen::Ref<const Eigen::Array<std::uint8_t, 1, Eigen::Dynamic>>
cell_face_rotations = mesh.topology().get_face_rotations(cells[0]);

// Tabulate tensor
const Eigen::Ref<const Eigen::Array<bool, Eigen::Dynamic, 1>> e_ref_cell
= cell_edge_reflections.col(cells[0]);
const Eigen::Ref<const Eigen::Array<bool, Eigen::Dynamic, 1>> f_ref_cell
= cell_face_reflections.col(cells[0]);
const Eigen::Ref<const Eigen::Array<uint8_t, Eigen::Dynamic, 1>> f_rot_cell
= cell_face_rotations.col(cells[0]);
Ae.setZero(dmapjoint0.size(), dmapjoint1.size());
fn(Ae.data(), coeff_array.data(), constant_values.data(),
coordinate_dofs.data(), local_facet.data(), perm.data(),
cell_edge_reflections.data(), cell_face_reflections.data(),
cell_face_rotations.data());
e_ref_cell.data(), f_ref_cell.data(), f_rot_cell.data());

// Zero rows/columns for essential bcs
if (!bc0.empty())
Expand Down
92 changes: 58 additions & 34 deletions cpp/dolfinx/fem/assemble_scalar_impl.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,6 +115,14 @@ PetscScalar fem::impl::assemble_cells(
Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
coordinate_dofs(num_dofs_g, gdim);

const Eigen::Ref<const Eigen::Array<bool, Eigen::Dynamic, Eigen::Dynamic>>
cell_edge_reflections = mesh.topology().get_edge_reflections();
const Eigen::Ref<const Eigen::Array<bool, Eigen::Dynamic, Eigen::Dynamic>>
cell_face_reflections = mesh.topology().get_face_reflections();
const Eigen::Ref<
const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>>
cell_face_rotations = mesh.topology().get_face_rotations();

// Iterate over all cells
PetscScalar value(0);
for (auto& cell : active_cells)
Expand All @@ -124,17 +132,17 @@ PetscScalar fem::impl::assemble_cells(
for (int j = 0; j < gdim; ++j)
coordinate_dofs(i, j) = x_g(cell_g[pos_g[cell] + i], j);

const Eigen::Ref<const Eigen::Array<bool, 1, Eigen::Dynamic>>
cell_edge_reflections = mesh.topology().get_edge_reflections(cell);
const Eigen::Ref<const Eigen::Array<bool, 1, Eigen::Dynamic>>
cell_face_reflections = mesh.topology().get_face_reflections(cell);
const Eigen::Ref<const Eigen::Array<std::uint8_t, 1, Eigen::Dynamic>>
cell_face_rotations = mesh.topology().get_face_rotations(cell);

auto coeff_cell = coeffs.row(cell);
const Eigen::Ref<const Eigen::Array<bool, Eigen::Dynamic, 1>> e_ref_cell
= cell_edge_reflections.col(cell);
const Eigen::Ref<const Eigen::Array<bool, Eigen::Dynamic, 1>> f_ref_cell
= cell_face_reflections.col(cell);
const Eigen::Ref<const Eigen::Array<uint8_t, Eigen::Dynamic, 1>> f_rot_cell
= cell_face_rotations.col(cell);

fn(&value, coeff_cell.data(), constant_values.data(),
coordinate_dofs.data(), nullptr, nullptr, cell_edge_reflections.data(),
cell_face_reflections.data(), cell_face_rotations.data());
coordinate_dofs.data(), nullptr, nullptr, e_ref_cell.data(),
f_ref_cell.data(), f_rot_cell.data());
}

return value;
Expand Down Expand Up @@ -172,6 +180,17 @@ PetscScalar fem::impl::assemble_exterior_facets(
Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
coordinate_dofs(num_dofs_g, gdim);

const Eigen::Ref<
const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>>
perms = mesh.topology().get_facet_permutations();

const Eigen::Ref<const Eigen::Array<bool, Eigen::Dynamic, Eigen::Dynamic>>
cell_edge_reflections = mesh.topology().get_edge_reflections();
const Eigen::Ref<const Eigen::Array<bool, Eigen::Dynamic, Eigen::Dynamic>>
cell_face_reflections = mesh.topology().get_face_reflections();
const Eigen::Ref<
const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>>
cell_face_rotations = mesh.topology().get_face_rotations();
auto f_to_c = mesh.topology().connectivity(tdim - 1, tdim);
assert(f_to_c);
auto c_to_f = mesh.topology().connectivity(tdim, tdim - 1);
Expand All @@ -191,25 +210,22 @@ PetscScalar fem::impl::assemble_exterior_facets(
assert(it != (facets.data() + facets.rows()));
const int local_facet = std::distance(facets.data(), it);

const std::uint8_t perm
= mesh.topology().get_facet_permutation(cell, tdim - 1, local_facet);

for (int i = 0; i < num_dofs_g; ++i)
for (int j = 0; j < gdim; ++j)
coordinate_dofs(i, j) = x_g(cell_g[pos_g[cell] + i], j);

const Eigen::Ref<const Eigen::Array<bool, 1, Eigen::Dynamic>>
cell_edge_reflections = mesh.topology().get_edge_reflections(cell);
const Eigen::Ref<const Eigen::Array<bool, 1, Eigen::Dynamic>>
cell_face_reflections = mesh.topology().get_face_reflections(cell);
const Eigen::Ref<const Eigen::Array<std::uint8_t, 1, Eigen::Dynamic>>
cell_face_rotations = mesh.topology().get_face_rotations(cell);

auto coeff_cell = coeffs.row(cell);
const Eigen::Ref<const Eigen::Array<bool, Eigen::Dynamic, 1>> e_ref_cell
= cell_edge_reflections.col(cell);
const Eigen::Ref<const Eigen::Array<bool, Eigen::Dynamic, 1>> f_ref_cell
= cell_face_reflections.col(cell);
const Eigen::Ref<const Eigen::Array<uint8_t, Eigen::Dynamic, 1>> f_rot_cell
= cell_face_rotations.col(cell);
const std::uint8_t perm = perms(local_facet, cell);

fn(&value, coeff_cell.data(), constant_values.data(),
coordinate_dofs.data(), &local_facet, &perm,
cell_edge_reflections.data(), cell_face_reflections.data(),
cell_face_rotations.data());
coordinate_dofs.data(), &local_facet, &perm, e_ref_cell.data(),
f_ref_cell.data(), f_rot_cell.data());
}

return value;
Expand Down Expand Up @@ -250,6 +266,18 @@ PetscScalar fem::impl::assemble_interior_facets(
Eigen::Array<PetscScalar, Eigen::Dynamic, 1> coeff_array(2 * offsets.back());
assert(offsets.back() == coeffs.cols());

const Eigen::Ref<
const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>>
perms = mesh.topology().get_facet_permutations();

const Eigen::Ref<const Eigen::Array<bool, Eigen::Dynamic, Eigen::Dynamic>>
cell_edge_reflections = mesh.topology().get_edge_reflections();
const Eigen::Ref<const Eigen::Array<bool, Eigen::Dynamic, Eigen::Dynamic>>
cell_face_reflections = mesh.topology().get_face_reflections();
const Eigen::Ref<
const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>>
cell_face_rotations = mesh.topology().get_face_rotations();

auto f_to_c = mesh.topology().connectivity(tdim - 1, tdim);
assert(f_to_c);
auto c_to_f = mesh.topology().connectivity(tdim, tdim - 1);
Expand All @@ -274,16 +302,7 @@ PetscScalar fem::impl::assemble_interior_facets(
}

const std::array<std::uint8_t, 2> perm
= {mesh.topology().get_facet_permutation(cells[0], tdim - 1,
local_facet[0]),
mesh.topology().get_facet_permutation(cells[1], tdim - 1,
local_facet[1])};
const Eigen::Ref<const Eigen::Array<bool, 1, Eigen::Dynamic>>
cell_edge_reflections = mesh.topology().get_edge_reflections(cells[0]);
const Eigen::Ref<const Eigen::Array<bool, 1, Eigen::Dynamic>>
cell_face_reflections = mesh.topology().get_face_reflections(cells[0]);
const Eigen::Ref<const Eigen::Array<std::uint8_t, 1, Eigen::Dynamic>>
cell_face_rotations = mesh.topology().get_face_rotations(cells[0]);
= {perms(local_facet[0], cells[0]), perms(local_facet[1], cells[1])};

for (int i = 0; i < num_dofs_g; ++i)
{
Expand All @@ -309,6 +328,12 @@ PetscScalar fem::impl::assemble_interior_facets(
// w[coefficient][restriction][dof]
auto coeff_cell0 = coeffs.row(cells[0]);
auto coeff_cell1 = coeffs.row(cells[1]);
const Eigen::Ref<const Eigen::Array<bool, Eigen::Dynamic, 1>> e_ref_cell
= cell_edge_reflections.col(cells[0]);
const Eigen::Ref<const Eigen::Array<bool, Eigen::Dynamic, 1>> f_ref_cell
= cell_face_reflections.col(cells[0]);
const Eigen::Ref<const Eigen::Array<uint8_t, Eigen::Dynamic, 1>> f_rot_cell
= cell_face_rotations.col(cells[0]);

// Loop over coefficients
for (std::size_t i = 0; i < offsets.size() - 1; ++i)
Expand All @@ -322,8 +347,7 @@ PetscScalar fem::impl::assemble_interior_facets(
}
fn(&value, coeff_array.data(), constant_values.data(),
coordinate_dofs.data(), local_facet.data(), perm.data(),
cell_edge_reflections.data(), cell_face_reflections.data(),
cell_face_rotations.data());
e_ref_cell.data(), f_ref_cell.data(), f_rot_cell.data());
}

return value;
Expand Down
Loading

0 comments on commit 3232b87

Please sign in to comment.