Skip to content

Adjust code to build with ReferenceCell being a template <dim>. #6904

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
bangerth wants to merge 1 commit into
base: main
Choose a base branch
from
+80 −31
Open

Adjust code to build with ReferenceCell being a template <dim>. #6904

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 9 additions & 0 deletions include/aspect/compat.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -444,6 +444,15 @@ namespace aspect
const double outer_radius);
#endif

// deal.II 9.8 made ReferenceCell a template class, whereas older versions
// had it as a non-template class. This is a problem.
// Rather than litter our own code base with #ifdefs, we can just define the
// templated class variant here for older deal.II versions, and then we can
// use the same code in all versions.
#if !DEAL_II_VERSION_GTE(9,8,0)
template <int dim> using ReferenceCell = dealii::ReferenceCell;
#endif

}

#endif
10 changes: 6 additions & 4 deletions include/aspect/simulator/assemblers/interface.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -487,7 +487,6 @@ namespace aspect
*/
namespace Assemblers
{

/**
* For a reference cell (which is typically obtained by asking the finite
* element to be used), determine how many interface matrices are needed.
Expand All @@ -497,25 +496,28 @@ namespace aspect
* accommodates the fact that the neighbors of a cell can all be refined,
* though they can only be refined once.
*/
template <int dim>
unsigned int
n_interface_matrices (const ReferenceCell &reference_cell);
n_interface_matrices (const ReferenceCell<dim> &reference_cell);

/**
* For a given reference cell, and a given face we are currently
* assembling on, return which element of an array of size
* `n_interface_matrices(reference_cell)` to use.
*/
template <int dim>
unsigned int
nth_interface_matrix (const ReferenceCell &reference_cell,
nth_interface_matrix (const ReferenceCell<dim> &reference_cell,
const unsigned int face);

/**
* For a given reference cell, and a given face and sub-face we are
* currently assembling on, return which element of an array of size
* `n_interface_matrices(reference_cell)` to use.
*/
template <int dim>
unsigned int
nth_interface_matrix (const ReferenceCell &reference_cell,
nth_interface_matrix (const ReferenceCell<dim> &reference_cell,
const unsigned int face,
const unsigned int sub_face);

Expand Down
25 changes: 22 additions & 3 deletions source/simulator/assemblers/interface.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -501,8 +501,9 @@ namespace aspect

namespace Assemblers
{
template <int dim>
unsigned int
n_interface_matrices (const ReferenceCell &reference_cell)
n_interface_matrices (const ReferenceCell<dim> &reference_cell)
{
// The current implementation assumes that all faces are
// the same; so no wedges or pyramids please.
Expand All @@ -520,8 +521,9 @@ namespace aspect



template <int dim>
unsigned int
nth_interface_matrix (const ReferenceCell &reference_cell,
nth_interface_matrix (const ReferenceCell<dim> &reference_cell,
const unsigned int face)
{
AssertIndexRange (face, reference_cell.n_faces());
Expand All @@ -531,8 +533,9 @@ namespace aspect



template <int dim>
unsigned int
nth_interface_matrix (const ReferenceCell &reference_cell,
nth_interface_matrix (const ReferenceCell<dim> &reference_cell,
const unsigned int face,
const unsigned int sub_face)
{
Expand Down Expand Up @@ -649,7 +652,23 @@ namespace aspect
template class Interface<dim>; \
template class AdvectionStabilizationInterface<dim>; \
template class Manager<dim>; \
\
template \
unsigned int \
n_interface_matrices (const ReferenceCell<dim> &reference_cell); \
\
template \
unsigned int \
nth_interface_matrix (const ReferenceCell<dim> &reference_cell, \
const unsigned int face); \
\
template \
unsigned int \
nth_interface_matrix (const ReferenceCell<dim> &reference_cell, \
const unsigned int face, \
const unsigned int sub_face); \
}

ASPECT_INSTANTIATE(INSTANTIATE)

#undef INSTANTIATE
Expand Down
67 changes: 43 additions & 24 deletions source/simulator/introspection.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -140,24 +140,33 @@ namespace aspect
template <int dim>
typename Introspection<dim>::Quadratures
setup_quadratures (const Parameters<dim> &parameters,
const ReferenceCell reference_cell)
const ReferenceCell<dim> reference_cell)
{
const auto get_gauss_type_quadrature = [&reference_cell](const unsigned int degree)
{
#if DEAL_II_VERSION_GTE(9,8,0)
return reference_cell.get_gauss_type_quadrature(degree);
#else
return reference_cell.get_gauss_type_quadrature<dim>(degree);
#endif
};

typename Introspection<dim>::Quadratures quadratures;

quadratures.velocities = reference_cell.get_gauss_type_quadrature<dim>(parameters.stokes_velocity_degree+1);
quadratures.pressure = reference_cell.get_gauss_type_quadrature<dim>(parameters.use_equal_order_interpolation_for_stokes
?
parameters.stokes_velocity_degree+1
:
parameters.stokes_velocity_degree);
quadratures.temperature = reference_cell.get_gauss_type_quadrature<dim>(parameters.temperature_degree+1);
quadratures.compositional_field_max = reference_cell.get_gauss_type_quadrature<dim>(parameters.max_composition_degree+1);
quadratures.velocities = get_gauss_type_quadrature(parameters.stokes_velocity_degree+1);
quadratures.pressure = get_gauss_type_quadrature(parameters.use_equal_order_interpolation_for_stokes
?
parameters.stokes_velocity_degree+1
:
parameters.stokes_velocity_degree);
quadratures.temperature = get_gauss_type_quadrature(parameters.temperature_degree+1);
quadratures.compositional_field_max = get_gauss_type_quadrature(parameters.max_composition_degree+1);
for (const auto &degree: parameters.composition_degrees)
quadratures.compositional_fields.emplace_back(reference_cell.get_gauss_type_quadrature<dim>(degree+1));
quadratures.system = reference_cell.get_gauss_type_quadrature<dim>(std::max({parameters.stokes_velocity_degree,
parameters.temperature_degree,
parameters.max_composition_degree
}) + 1);
quadratures.compositional_fields.emplace_back(get_gauss_type_quadrature(degree+1));
quadratures.system = get_gauss_type_quadrature(std::max({parameters.stokes_velocity_degree,
parameters.temperature_degree,
parameters.max_composition_degree
}) + 1);

return quadratures;
}
Expand All @@ -167,19 +176,29 @@ namespace aspect
template <int dim>
typename Introspection<dim>::FaceQuadratures
setup_face_quadratures (const Parameters<dim> &parameters,
const ReferenceCell reference_cell)
const ReferenceCell<dim> reference_cell)
{
const auto get_gauss_type_face_quadrature = [&reference_cell](const unsigned int degree)
{
#if DEAL_II_VERSION_GTE(9,8,0)
return reference_cell.face_reference_cell(0).get_gauss_type_quadrature(degree);
#else
return reference_cell.face_reference_cell(0).get_gauss_type_quadrature<dim-1>(degree);
#endif
};


typename Introspection<dim>::FaceQuadratures quadratures;

quadratures.velocities = reference_cell.face_reference_cell(0).get_gauss_type_quadrature<dim-1>(parameters.stokes_velocity_degree+1);
quadratures.pressure = reference_cell.face_reference_cell(0).get_gauss_type_quadrature<dim-1>(parameters.use_equal_order_interpolation_for_stokes
?
parameters.stokes_velocity_degree+1
:
parameters.stokes_velocity_degree);
quadratures.temperature = reference_cell.face_reference_cell(0).get_gauss_type_quadrature<dim-1>(parameters.temperature_degree+1);
quadratures.compositional_fields = reference_cell.face_reference_cell(0).get_gauss_type_quadrature<dim-1>(parameters.max_composition_degree+1);
quadratures.system = reference_cell.face_reference_cell(0).get_gauss_type_quadrature<dim-1>(
quadratures.velocities = get_gauss_type_face_quadrature(parameters.stokes_velocity_degree+1);
quadratures.pressure = get_gauss_type_face_quadrature(parameters.use_equal_order_interpolation_for_stokes
?
parameters.stokes_velocity_degree+1
:
parameters.stokes_velocity_degree);
quadratures.temperature = get_gauss_type_face_quadrature(parameters.temperature_degree+1);
quadratures.compositional_fields = get_gauss_type_face_quadrature(parameters.max_composition_degree+1);
quadratures.system = get_gauss_type_face_quadrature(
std::max({parameters.stokes_velocity_degree, parameters.temperature_degree, parameters.max_composition_degree}) + 1);

return quadratures;
Expand Down
Loading