Skip to content

Commit

Permalink
NonMatching::FEValues: accept CellAccessor
Browse files Browse the repository at this point in the history
  • Loading branch information
peterrum committed Apr 28, 2024
1 parent e9eb5ab commit 91a1336
Show file tree
Hide file tree
Showing 3 changed files with 109 additions and 12 deletions.
32 changes: 31 additions & 1 deletion include/deal.II/non_matching/fe_values.h
Original file line number Diff line number Diff line change
Expand Up @@ -223,7 +223,27 @@ namespace NonMatching
template <bool level_dof_access>
void
reinit(
const TriaIterator<DoFCellAccessor<dim, dim, level_dof_access>> &cell);
const TriaIterator<DoFCellAccessor<dim, dim, level_dof_access>> &cell,
const unsigned int q_index = numbers::invalid_unsigned_int,
const unsigned int mapping_index = numbers::invalid_unsigned_int);

/**
* Like the previous function, but for non-DoFHandler iterators. The reason
* this function exists is so that one can use hp::FEValues for
* Triangulation objects too.
*
* Since <code>cell-@>active_fe_index()</code> doesn't make sense for
* triangulation iterators, this function chooses the zero-th finite
* element, mapping, and quadrature object from the relevant constructions
* passed to the constructor of this object. The only exception is if you
* specify a value different from the default value for any of these last
* three arguments.
*/
void
reinit(const TriaIterator<CellAccessor<dim, dim>> &cell,
const unsigned int q_index = numbers::invalid_unsigned_int,
const unsigned int mapping_index = numbers::invalid_unsigned_int,
const unsigned int fe_index = numbers::invalid_unsigned_int);

/**
* Return an dealii::FEValues object reinitialized with a quadrature for the
Expand Down Expand Up @@ -258,6 +278,16 @@ namespace NonMatching
get_surface_fe_values() const;

private:
/**
* Internal function called by the reinit() functions.
*/
template <typename CellIteratorType>
void
reinit_internal(const CellIteratorType &cell,
const unsigned int q_index,
const unsigned int mapping_index,
const unsigned int fe_index);

/**
* Do work common to the constructors. The incoming QCollection should be
* quadratures integrating over $[0, 1]^{dim}$. These will be used on the
Expand Down
85 changes: 75 additions & 10 deletions source/non_matching/fe_values.cc
Original file line number Diff line number Diff line change
Expand Up @@ -148,10 +148,68 @@ namespace NonMatching
template <bool level_dof_access>
void
FEValues<dim>::reinit(
const TriaIterator<DoFCellAccessor<dim, dim, level_dof_access>> &cell)
const TriaIterator<DoFCellAccessor<dim, dim, level_dof_access>> &cell,
const unsigned int q_index,
const unsigned int mapping_index)
{
this->reinit_internal(cell,
q_index,
mapping_index,
cell->active_fe_index());
}



template <int dim>
void
FEValues<dim>::reinit(const TriaIterator<CellAccessor<dim, dim>> &cell,
const unsigned int q_index,
const unsigned int mapping_index,
const unsigned int fe_index)
{
this->reinit_internal(cell, q_index, mapping_index, fe_index);
}



template <int dim>
template <typename CellIteratorType>
void
FEValues<dim>::reinit_internal(const CellIteratorType &cell,
const unsigned int q_index_in,
const unsigned int mapping_index_in,
const unsigned int fe_index_in)
{
current_cell_location = mesh_classifier->location_to_level_set(cell);
active_fe_index = cell->active_fe_index();
this->active_fe_index = fe_index_in;

unsigned int mapping_index = mapping_index_in;
unsigned int q_index = q_index_in;
unsigned int q_index_1D = q_index_in;

if (mapping_index == numbers::invalid_unsigned_int)
{
if (mapping_collection->size() > 1)
mapping_index = active_fe_index;
else
mapping_index = 0;
}

if (q_index == numbers::invalid_unsigned_int)
{
if (fe_values_inside_full_quadrature.size() > 1)
q_index = active_fe_index;
else
q_index = 0;
}

if (q_index_1D == numbers::invalid_unsigned_int)
{
if (q_collection_1D.size() > 1)
q_index_1D = active_fe_index;
else
q_index_1D = 0;
}

// These objects were created with a quadrature based on the previous cell
// and are thus no longer valid.
Expand All @@ -163,22 +221,29 @@ namespace NonMatching
{
case LocationToLevelSet::inside:
{
fe_values_inside_full_quadrature.at(active_fe_index)->reinit(cell);
Assert((active_fe_index == mapping_index) ||
((mapping_collection->size() == 1) &&
(mapping_index == 0)),
ExcNotImplemented());
Assert(active_fe_index == q_index, ExcNotImplemented());

fe_values_inside_full_quadrature.at(q_index)->reinit(cell);
break;
}
case LocationToLevelSet::outside:
{
fe_values_outside_full_quadrature.at(active_fe_index)->reinit(cell);
Assert((active_fe_index == mapping_index) ||
((mapping_collection->size() == 1) &&
(mapping_index == 0)),
ExcNotImplemented());
Assert(active_fe_index == q_index, ExcNotImplemented());

fe_values_outside_full_quadrature.at(q_index)->reinit(cell);
break;
}
case LocationToLevelSet::intersected:
{
const unsigned int mapping_index =
mapping_collection->size() > 1 ? active_fe_index : 0;

const unsigned int q1D_index =
q_collection_1D.size() > 1 ? active_fe_index : 0;
quadrature_generator.set_1D_quadrature(q1D_index);
quadrature_generator.set_1D_quadrature(q_index_1D);
quadrature_generator.generate(cell);

const Quadrature<dim> &inside_quadrature =
Expand Down
4 changes: 3 additions & 1 deletion source/non_matching/fe_values.inst.in
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,9 @@ for (deal_II_dimension : DIMENSIONS; lda : BOOL)
{
template void FEValues<deal_II_dimension>::reinit(
const TriaIterator<
DoFCellAccessor<deal_II_dimension, deal_II_dimension, lda>> &);
DoFCellAccessor<deal_II_dimension, deal_II_dimension, lda>> &,
const unsigned int,
const unsigned int);

template void FEInterfaceValues<deal_II_dimension>::do_reinit(
const TriaIterator<
Expand Down

0 comments on commit 91a1336

Please sign in to comment.