inline std::pair<unsigned int, unsigned int>
FiniteElement<dim, spacedim>::face_system_to_component_index(
const unsigned int index,
- const unsigned int) const
+ const unsigned int face_no) const
{
- AssertIndexRange(index, face_system_to_component_table[0].size());
+ AssertIndexRange(
+ index,
+ face_system_to_component_table[this->n_unique_faces() == 1 ? 0 : face_no]
+ .size());
// in debug mode, check whether the
// function is primitive, since
//
// in 1d, the face index is equal
// to the cell index
- Assert(is_primitive(this->face_to_cell_index(index, 0)),
+ Assert(is_primitive(this->face_to_cell_index(index, face_no)),
(typename FiniteElement<dim, spacedim>::ExcShapeFunctionNotPrimitive(
index)));
- return face_system_to_component_table[0][index];
+ return face_system_to_component_table[this->n_unique_faces() == 1 ?
+ 0 :
+ face_no][index];
}
inline std::pair<std::pair<unsigned int, unsigned int>, unsigned int>
FiniteElement<dim, spacedim>::face_system_to_base_index(
const unsigned int index,
- const unsigned int) const
+ const unsigned int face_no) const
{
- AssertIndexRange(index, face_system_to_base_table[0].size());
- return face_system_to_base_table[0][index];
+ AssertIndexRange(
+ index,
+ face_system_to_base_table[this->n_unique_faces() == 1 ? 0 : face_no]
+ .size());
+ return face_system_to_base_table[this->n_unique_faces() == 1 ? 0 : face_no]
+ [index];
}
// are enumerated on the reference cell
if (cell_dof_index < this->get_first_line_index())
return GeometryPrimitive::vertex;
- else if (cell_dof_index < this->get_first_quad_index())
+ else if (cell_dof_index < this->get_first_quad_index(0))
return GeometryPrimitive::line;
else if (cell_dof_index < this->get_first_hex_index())
return GeometryPrimitive::quad;
case 3:
{
// Fill data for quad shape functions
- if (this->n_dofs_per_quad() != 0)
+ if (this->n_dofs_per_quad(face_no) != 0)
{
const unsigned int foffset =
- this->get_first_quad_index() +
- this->n_dofs_per_quad() * face_no;
- for (unsigned int k = 0; k < this->n_dofs_per_quad(); ++k)
+ this->get_first_quad_index(face_no);
+ for (unsigned int k = 0; k < this->n_dofs_per_quad(face_no);
+ ++k)
output_data.shape_values(foffset + k, i) =
fe_data
- .shape_values[k + this->get_first_face_quad_index()]
- [i];
+ .shape_values[k + this->get_first_face_quad_index(
+ face_no)][i];
}
}
DEAL_II_FALLTHROUGH;
line) *
this->n_dofs_per_line() +
k,
- i) = fe_data.shape_values
- [k + (line * this->n_dofs_per_line()) +
- this->get_first_face_line_index()][i];
+ i) =
+ fe_data.shape_values
+ [k + (line * this->n_dofs_per_line()) +
+ this->get_first_face_line_index(face_no)][i];
}
}
}
*/
virtual std::vector<std::pair<unsigned int, unsigned int>>
hp_quad_dof_identities(const FiniteElement<dim, spacedim> &fe_other,
- const unsigned int = 0) const override;
+ const unsigned int face_no = 0) const override;
/**
* @copydoc FiniteElement::compare_for_domination()
*/
template <int structdim>
std::vector<std::pair<unsigned int, unsigned int>>
- hp_object_dof_identities(const FiniteElement<dim, spacedim> &fe_other) const;
+ hp_object_dof_identities(const FiniteElement<dim, spacedim> &fe_other,
+ const unsigned int face_no = 0) const;
/**
* Usually: Fields of cell-independent data.
std::vector<std::pair<unsigned int, unsigned int>>
& face_system_to_component_table,
const FiniteElement<dim, spacedim> &finite_element,
- const bool do_tensor_product = true);
+ const bool do_tensor_product = true,
+ const unsigned int face_no = 0 /*TODO*/);
} // namespace Compositing
for (unsigned int i = 0; i < fes.size(); i++)
if (multiplicities[i] > 0)
{
+ // TODO: the implementation makes the assumption that all faces have
+ // the same number of dofs -> don't construct DPO but
+ // PrecomputedFiniteElementData
+ AssertDimension(fes[i]->n_unique_quads(), 1);
+
multiplied_dofs_per_vertex +=
fes[i]->n_dofs_per_vertex() * multiplicities[i];
multiplied_dofs_per_line +=
fes[i]->n_dofs_per_line() * multiplicities[i];
multiplied_dofs_per_quad +=
- fes[i]->n_dofs_per_quad() * multiplicities[i];
+ fes[i]->n_dofs_per_quad(0) * multiplicities[i];
multiplied_dofs_per_hex +=
fes[i]->n_dofs_per_hex() * multiplicities[i];
for (unsigned int base = 0; base < fes.size(); ++base)
for (unsigned int m = 0; m < multiplicities[base]; ++m)
for (unsigned int local_index = 0;
- local_index < fes[base]->n_dofs_per_quad();
+ local_index < fes[base]->n_dofs_per_quad(quad_number);
++local_index, ++total_index)
{
const unsigned int index_in_base =
- (fes[base]->n_dofs_per_quad() * quad_number + local_index +
- fes[base]->get_first_quad_index());
+ local_index + fes[base]->get_first_quad_index(quad_number);
Assert(index_in_base < fes[base]->n_dofs_per_cell(),
ExcInternalError());
comp_start +=
fes[base]->n_components() * do_tensor_product)
for (unsigned int local_index = 0;
- local_index < fes[base]->n_dofs_per_quad();
+ local_index < fes[base]->n_dofs_per_quad(quad_number);
++local_index, ++total_index)
{
const unsigned int index_in_base =
- (fes[base]->n_dofs_per_quad() * quad_number + local_index +
- fes[base]->get_first_quad_index());
+ local_index + fes[base]->get_first_quad_index(quad_number);
Assert(comp_start + fes[base]->n_components() <=
retval[total_index].size(),
fe.base_element(base).n_components() *
do_tensor_product)
for (unsigned int local_index = 0;
- local_index < fe.base_element(base).n_dofs_per_quad();
+ local_index <
+ fe.base_element(base).n_dofs_per_quad(quad_number);
++local_index, ++total_index)
{
const unsigned int index_in_base =
- (fe.base_element(base).n_dofs_per_quad() * quad_number +
- local_index +
- fe.base_element(base).get_first_quad_index());
+ local_index +
+ fe.base_element(base).get_first_quad_index(quad_number);
system_to_base_table[total_index] =
std::make_pair(std::make_pair(base, m), index_in_base);
std::vector<std::pair<unsigned int, unsigned int>>
& face_system_to_component_table,
const FiniteElement<dim, spacedim> &fe,
- const bool do_tensor_product)
+ const bool do_tensor_product,
+ const unsigned int face_no)
{
// Initialize index tables. do this in the same way as done for the cell
// tables, except that we now loop over the objects of faces
unsigned int total_index = 0;
for (unsigned int vertex_number = 0;
vertex_number <
- ReferenceCell::internal::Info::get_face(fe.reference_cell_type(), 0)
+ ReferenceCell::internal::Info::get_face(fe.reference_cell_type(),
+ face_no)
.n_vertices();
++vertex_number)
{
{
const unsigned int comp_in_base =
fe.base_element(base)
- .face_system_to_component_index(face_index_in_base)
+ .face_system_to_component_index(face_index_in_base,
+ face_no)
.first;
const unsigned int comp = comp_start + comp_in_base;
const unsigned int face_index_in_comp =
fe.base_element(base)
- .face_system_to_component_index(face_index_in_base)
+ .face_system_to_component_index(face_index_in_base,
+ face_no)
.second;
face_system_to_component_table[total_index] =
std::make_pair(comp, face_index_in_comp);
// 2. Lines
for (unsigned int line_number = 0;
line_number <
- ReferenceCell::internal::Info::get_face(fe.reference_cell_type(), 0)
+ ReferenceCell::internal::Info::get_face(fe.reference_cell_type(),
+ face_no)
.n_lines();
++line_number)
{
fe.base_element(base).get_first_line_index());
const unsigned int face_index_in_base =
- (fe.base_element(base).get_first_face_line_index() +
+ (fe.base_element(base).get_first_face_line_index(face_no) +
fe.base_element(base).n_dofs_per_line() * line_number +
local_index);
{
const unsigned int comp_in_base =
fe.base_element(base)
- .face_system_to_component_index(face_index_in_base)
+ .face_system_to_component_index(face_index_in_base,
+ face_no)
.first;
const unsigned int comp = comp_start + comp_in_base;
const unsigned int face_index_in_comp =
fe.base_element(base)
- .face_system_to_component_index(face_index_in_base)
+ .face_system_to_component_index(face_index_in_base,
+ face_no)
.second;
face_system_to_component_table[total_index] =
std::make_pair(comp, face_index_in_comp);
fe.base_element(base).n_components() *
do_tensor_product)
for (unsigned int local_index = 0;
- local_index < fe.base_element(base).n_dofs_per_quad();
+ local_index < fe.base_element(base).n_dofs_per_quad(face_no);
++local_index, ++total_index)
{
// do everything alike for this type of object
const unsigned int index_in_base =
(local_index +
- fe.base_element(base).get_first_quad_index());
+ fe.base_element(base).get_first_quad_index(face_no));
const unsigned int face_index_in_base =
- (fe.base_element(base).get_first_face_quad_index() +
+ (fe.base_element(base).get_first_face_quad_index(face_no) +
local_index);
face_system_to_base_table[total_index] =
{
const unsigned int comp_in_base =
fe.base_element(base)
- .face_system_to_component_index(face_index_in_base)
+ .face_system_to_component_index(face_index_in_base,
+ face_no)
.first;
const unsigned int comp = comp_start + comp_in_base;
const unsigned int face_index_in_comp =
fe.base_element(base)
- .face_system_to_component_index(face_index_in_base)
+ .face_system_to_component_index(face_index_in_base,
+ face_no)
.second;
face_system_to_component_table[total_index] =
std::make_pair(comp, face_index_in_comp);
non_primitive_index;
}
}
- Assert(total_index == fe.n_dofs_per_face(), ExcInternalError());
+ Assert(total_index == fe.n_dofs_per_face(face_no), ExcInternalError());
Assert(total_index == face_system_to_component_table.size(),
ExcInternalError());
Assert(total_index == face_system_to_base_table.size(),
const unsigned int face_fine,
const double threshold)
{
+ const unsigned int face_no = face_coarse;
+
Assert(face_coarse == 0, ExcNotImplemented());
Assert(face_fine == 0, ExcNotImplemented());
const unsigned int nc = GeometryInfo<dim>::max_children_per_face;
- const unsigned int n = fe.n_dofs_per_face();
+ const unsigned int n = fe.n_dofs_per_face(face_no);
const unsigned int nd = fe.n_components();
const unsigned int degree = fe.degree;
{
unsigned int face_dof = 0;
for (unsigned int i = 0;
- i <
- ReferenceCell::internal::Info::get_face(fe.reference_cell_type(), 0)
- .n_vertices();
+ i < ReferenceCell::internal::Info::get_face(fe.reference_cell_type(),
+ face_no)
+ .n_vertices();
++i)
{
const unsigned int offset_c =
}
}
- for (unsigned int i = 1;
- i <=
- ReferenceCell::internal::Info::get_face(fe.reference_cell_type(), 0)
- .n_lines();
+ for (unsigned int i = 1; i <= ReferenceCell::internal::Info::get_face(
+ fe.reference_cell_type(), face_no)
+ .n_lines();
++i)
{
const unsigned int offset_c =
if (dim == 3)
{
- const unsigned int offset_c =
- fe.get_first_quad_index() + face_coarse * fe.n_dofs_per_quad();
- const unsigned int offset_f =
- fe.get_first_quad_index() + face_fine * fe.n_dofs_per_quad();
- for (unsigned int j = 0; j < fe.n_dofs_per_quad(); ++j)
+ const unsigned int offset_c = fe.get_first_quad_index(face_coarse);
+ const unsigned int offset_f = fe.get_first_quad_index(face_fine);
+ for (unsigned int j = 0; j < fe.n_dofs_per_quad(face_no); ++j)
{
face_c_dofs[face_dof] = offset_c + j;
face_f_dofs[face_dof] = offset_f + j;
++face_dof;
}
}
- Assert(face_dof == fe.n_dofs_per_face(), ExcInternalError());
+ Assert(face_dof == fe.n_dofs_per_face(face_no), ExcInternalError());
}
// Set up meshes, one with a single
compute_all_non_local_data(dof2, u2_relevant);
// exclude dofs on more refined ghosted cells
- const FiniteElement<dim, spacedim> &fe = dof2.get_fe();
- const unsigned int dofs_per_face = fe.n_dofs_per_face();
- if (dofs_per_face > 0)
+ const FiniteElement<dim, spacedim> &fe = dof2.get_fe();
+ if (fe.max_dofs_per_face() > 0)
{
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
std::vector<types::global_dof_index> indices(dofs_per_cell);
const typename DoFHandler<dim, spacedim>::cell_iterator
neighbor = cell->neighbor(face);
if (neighbor->level() != cell->level())
- for (unsigned int i = 0; i < dofs_per_face; ++i)
+ for (unsigned int i = 0; i < fe.n_dofs_per_face(face);
+ ++i)
{
const types::global_dof_index index =
indices[fe.face_to_cell_index(i, face)];
for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j)
system_to_component_table[j] = std::pair<unsigned, unsigned>(0, j);
- face_system_to_component_table.resize(1);
- face_system_to_component_table[0].resize(this->n_dofs_per_face());
- for (unsigned int j = 0; j < this->n_dofs_per_face(); ++j)
- face_system_to_component_table[0][j] =
- std::pair<unsigned, unsigned>(0, j);
+ face_system_to_component_table.resize(this->n_unique_faces());
+ for (unsigned int f = 0; f < this->n_unique_faces(); ++f)
+ {
+ face_system_to_component_table[f].resize(this->n_dofs_per_face(f));
+ for (unsigned int j = 0; j < this->n_dofs_per_face(f); ++j)
+ face_system_to_component_table[f][j] =
+ std::pair<unsigned, unsigned>(0, j);
+ }
}
for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j)
system_to_base_table[j] = std::make_pair(std::make_pair(0U, 0U), j);
- face_system_to_base_table.resize(1);
- face_system_to_base_table[0].resize(this->n_dofs_per_face());
- for (unsigned int j = 0; j < this->n_dofs_per_face(); ++j)
- face_system_to_base_table[0][j] = std::make_pair(std::make_pair(0U, 0U), j);
+ face_system_to_base_table.resize(this->n_unique_faces());
+ for (unsigned int f = 0; f < this->n_unique_faces(); ++f)
+ {
+ face_system_to_base_table[f].resize(this->n_dofs_per_face(f));
+ for (unsigned int j = 0; j < this->n_dofs_per_face(f); ++j)
+ face_system_to_base_table[f][j] =
+ std::make_pair(std::make_pair(0U, 0U), j);
+ }
// Fill with default value; may be changed by constructor of derived class.
base_to_block_indices.reinit(1, 1);
FullMatrix<double>());
}
+
if (dim == 3)
{
- adjust_quad_dof_index_for_face_orientation_table.resize(1);
+ adjust_quad_dof_index_for_face_orientation_table.resize(
+ this->n_unique_quads());
- adjust_quad_dof_index_for_face_orientation_table[0] =
- Table<2, int>(this->n_dofs_per_quad(), 8);
- adjust_quad_dof_index_for_face_orientation_table[0].fill(0);
+ for (unsigned int f = 0; f < this->n_unique_quads(); ++f)
+ {
+ adjust_quad_dof_index_for_face_orientation_table[f] = Table<2, int>(
+ this->n_dofs_per_quad(f),
+ ReferenceCell::internal::Info::get_cell(this->reference_cell_type())
+ .face_reference_cell_type(f) == ReferenceCell::Type::Quad ?
+ 8 :
+ 6);
+ adjust_quad_dof_index_for_face_orientation_table[f].fill(0);
+ }
}
- unit_face_support_points.resize(1);
- generalized_face_support_points.resize(1);
+ unit_face_support_points.resize(this->n_unique_faces());
+ generalized_face_support_points.resize(this->n_unique_faces());
}
const auto &refence_cell =
ReferenceCell::internal::Info::get_cell(this->reference_cell_type());
- AssertIndexRange(face_index, this->n_dofs_per_face());
+ AssertIndexRange(face_index, this->n_dofs_per_face(face));
AssertIndexRange(face, refence_cell.n_faces());
// TODO: we could presumably solve the 3d case below using the
// other than standard orientation
if ((face_orientation != true) || (face_flip != false) ||
(face_rotation != false))
- Assert((this->n_dofs_per_line() <= 1) && (this->n_dofs_per_quad() <= 1),
+ Assert((this->n_dofs_per_line() <= 1) && (this->n_dofs_per_quad(face) <= 1),
ExcMessage(
"The function in this base class can not handle this case. "
"Rather, the derived class you are using must provide "
// we need to distinguish between DoFs on vertices, lines and in 3d quads.
// do so in a sequence of if-else statements
- if (face_index < this->get_first_face_line_index())
+ if (face_index < this->get_first_face_line_index(face))
// DoF is on a vertex
{
// get the number of the vertex on the face that corresponds to this DoF,
this->n_dofs_per_vertex() +
dof_index_on_vertex);
}
- else if (face_index < this->get_first_face_quad_index())
+ else if (face_index < this->get_first_face_quad_index(face))
// DoF is on a face
{
// do the same kind of translation as before. we need to only consider
// DoFs on the lines, i.e., ignoring those on the vertices
- const unsigned int index = face_index - this->get_first_face_line_index();
+ const unsigned int index =
+ face_index - this->get_first_face_line_index(face);
const unsigned int face_line = index / this->n_dofs_per_line();
const unsigned int dof_index_on_line = index % this->n_dofs_per_line();
Assert(dim >= 3, ExcInternalError());
// ignore vertex and line dofs
- const unsigned int index = face_index - this->get_first_face_quad_index();
+ const unsigned int index =
+ face_index - this->get_first_face_quad_index(face);
- return (this->get_first_quad_index() + face * this->n_dofs_per_quad() +
- index);
+ return (this->get_first_quad_index(face) + index);
}
}
unsigned int
FiniteElement<dim, spacedim>::adjust_quad_dof_index_for_face_orientation(
const unsigned int index,
- const unsigned int face_no,
+ const unsigned int face,
const bool face_orientation,
const bool face_flip,
const bool face_rotation) const
{
- (void)face_no;
-
// general template for 1D and 2D: not
// implemented. in fact, the function
// shouldn't even be called unless we are
// in 3d), so we don't need the table, but
// the function should also not have been
// called
- AssertIndexRange(index, this->n_dofs_per_quad());
- Assert(adjust_quad_dof_index_for_face_orientation_table[0].n_elements() ==
- 8 * this->n_dofs_per_quad(),
+ AssertIndexRange(index, this->n_dofs_per_quad(face));
+ Assert(adjust_quad_dof_index_for_face_orientation_table
+ [this->n_unique_quads() == 1 ? 0 : face]
+ .n_elements() ==
+ (ReferenceCell::internal::Info::get_cell(this->reference_cell_type())
+ .face_reference_cell_type(face) == ReferenceCell::Type::Quad ?
+ 8 :
+ 6) *
+ this->n_dofs_per_quad(face),
ExcInternalError());
- return index + adjust_quad_dof_index_for_face_orientation_table[0](
- index, 4 * face_orientation + 2 * face_flip + face_rotation);
+ return index +
+ adjust_quad_dof_index_for_face_orientation_table
+ [this->n_unique_quads() == 1 ? 0 : face](
+ index, 4 * face_orientation + 2 * face_flip + face_rotation);
}
FiniteElement<dim, spacedim>::constraints_are_implemented(
const internal::SubfaceCase<dim> &subface_case) const
{
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(this->n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
if (subface_case == internal::SubfaceCase<dim>::case_isotropic)
- return (this->n_dofs_per_face() == 0) || (interface_constraints.m() != 0);
+ return (this->n_dofs_per_face(face_no) == 0) ||
+ (interface_constraints.m() != 0);
else
return false;
}
FiniteElement<dim, spacedim>::constraints(
const internal::SubfaceCase<dim> &subface_case) const
{
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(this->n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+ (void)face_no;
+
(void)subface_case;
Assert(subface_case == internal::SubfaceCase<dim>::case_isotropic,
ExcMessage("Constraints for this element are only implemented "
"(which is always the case in 2d, and in 3d requires "
"that the neighboring cell of a coarse cell presents "
"exactly four children on the common face)."));
- Assert((this->n_dofs_per_face() == 0) || (interface_constraints.m() != 0),
+ Assert((this->n_dofs_per_face(face_no) == 0) ||
+ (interface_constraints.m() != 0),
ExcMessage("The finite element for which you try to obtain "
"hanging node constraints does not appear to "
"implement them."));
TableIndices<2>
FiniteElement<dim, spacedim>::interface_constraints_size() const
{
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(this->n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
switch (dim)
{
case 1:
return {0U, 0U};
case 2:
return {this->n_dofs_per_vertex() + 2 * this->n_dofs_per_line(),
- this->n_dofs_per_face()};
+ this->n_dofs_per_face(face_no)};
case 3:
return {5 * this->n_dofs_per_vertex() + 12 * this->n_dofs_per_line() +
- 4 * this->n_dofs_per_quad(),
- this->n_dofs_per_face()};
+ 4 * this->n_dofs_per_quad(face_no),
+ this->n_dofs_per_face(face_no)};
default:
Assert(false, ExcNotImplemented());
}
template <int dim, int spacedim>
const std::vector<Point<dim - 1>> &
FiniteElement<dim, spacedim>::get_unit_face_support_points(
- const unsigned int) const
+ const unsigned int face_no) const
{
// a finite element may define
// support points, but only if
// there are as many as there are
// degrees of freedom on a face
- Assert((unit_face_support_points[0].size() == 0) ||
- (unit_face_support_points[0].size() == this->n_dofs_per_face()),
+ Assert((unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
+ .size() == 0) ||
+ (unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
+ .size() == this->n_dofs_per_face(face_no)),
ExcInternalError());
- return unit_face_support_points[0];
+ return unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no];
}
template <int dim, int spacedim>
bool
-FiniteElement<dim, spacedim>::has_face_support_points(const unsigned int) const
+FiniteElement<dim, spacedim>::has_face_support_points(
+ const unsigned int face_no) const
{
- return (unit_face_support_points[0].size() != 0);
+ return (unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
+ .size() != 0);
}
template <int dim, int spacedim>
Point<dim - 1>
-FiniteElement<dim, spacedim>::unit_face_support_point(const unsigned int index,
- const unsigned int) const
+FiniteElement<dim, spacedim>::unit_face_support_point(
+ const unsigned int index,
+ const unsigned int face_no) const
{
- AssertIndexRange(index, this->n_dofs_per_face());
- Assert(unit_face_support_points[0].size() == this->n_dofs_per_face(),
+ AssertIndexRange(index, this->n_dofs_per_face(face_no));
+ Assert(unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
+ .size() == this->n_dofs_per_face(face_no),
ExcFEHasNoSupportPoints());
- return unit_face_support_points[0][index];
+ return unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
+ [index];
}
initialize_restriction();
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(this->n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
// TODO[TL]: for anisotropic refinement we will probably need a table of
// submatrices with an array for each refine case
std::vector<FullMatrix<double>> face_embeddings(
1 << (dim - 1),
- FullMatrix<double>(this->n_dofs_per_face(), this->n_dofs_per_face()));
+ FullMatrix<double>(this->n_dofs_per_face(face_no),
+ this->n_dofs_per_face(face_no)));
// TODO: Something goes wrong there. The error of the least squares fit
// is to large ...
// FETools::compute_face_embedding_matrices(*this, face_embeddings.data(), 0,
// 0);
- this->interface_constraints.reinit((1 << (dim - 1)) * this->n_dofs_per_face(),
- this->n_dofs_per_face());
+ this->interface_constraints.reinit((1 << (dim - 1)) *
+ this->n_dofs_per_face(face_no),
+ this->n_dofs_per_face(face_no));
unsigned int target_row = 0;
for (const auto &face_embedding : face_embeddings)
for (unsigned int i = 0; i < face_embedding.m(); ++i)
QGauss<dim> cell_quadrature(deg + 2);
const unsigned int n_interior_points = cell_quadrature.size();
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(this->n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
unsigned int n_face_points = (dim > 1) ? 1 : 0;
// compute (deg+1)^(dim-1)
for (unsigned int d = 1; d < dim; ++d)
this->generalized_support_points.resize(
GeometryInfo<dim>::faces_per_cell * n_face_points + n_interior_points);
- this->generalized_face_support_points[0].resize(n_face_points);
+ this->generalized_face_support_points[face_no].resize(n_face_points);
// These might be required when the faces contribution is computed
for (unsigned int k = 0; k < n_face_points; ++k)
{
- this->generalized_face_support_points[0][k] = face_points.point(k);
+ this->generalized_face_support_points[face_no][k] =
+ face_points.point(k);
// Compute its quadrature
// contribution for each
// moment.
for (unsigned int k = 0; k < n_face_points; ++k)
for (unsigned int i_child = 0; i_child < this->n_dofs_per_cell();
++i_child)
- for (unsigned int i_face = 0; i_face < this->n_dofs_per_face();
+ for (unsigned int i_face = 0;
+ i_face < this->n_dofs_per_face(face);
++i_face)
{
// The quadrature
// subcell are NOT
// transformed, so we
// have to do it here.
- this->restriction[iso][child](face * this->n_dofs_per_face() +
- i_face,
- i_child) +=
+ this->restriction[iso][child](
+ face * this->n_dofs_per_face(face) + i_face, i_child) +=
Utilities::fixed_power<dim - 1>(.5) * q_sub.weight(k) *
cached_values_face(i_child, k) *
this->shape_value_component(
- face * this->n_dofs_per_face() + i_face,
+ face * this->n_dofs_per_face(face) + i_face,
q_sub.point(k),
GeometryInfo<dim>::unit_normal_direction[face]);
}
polynomials[dd] = new AnisotropicPolynomials<dim>(poly);
}
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(this->n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
QGauss<dim> q_cell(rt_order + 1);
const unsigned int start_cell_dofs =
- GeometryInfo<dim>::faces_per_cell * this->n_dofs_per_face();
+ GeometryInfo<dim>::faces_per_cell * this->n_dofs_per_face(face_no);
// Store shape values, since the
// evaluation suffers if not
for (unsigned int k = 0; k < n_face_points; ++k)
for (unsigned int i = 0; i < boundary_weights.size(1); ++i)
{
- nodal_values[i + face * this->n_dofs_per_face()] +=
+ nodal_values[i + face * this->n_dofs_per_face(face)] +=
boundary_weights(k, i) *
support_point_values[face * n_face_points + k][GeometryInfo<
dim>::unit_normal_direction[face]];
}
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(this->n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
const unsigned int start_cell_dofs =
- GeometryInfo<dim>::faces_per_cell * this->n_dofs_per_face();
+ GeometryInfo<dim>::faces_per_cell * this->n_dofs_per_face(face_no);
const unsigned int start_cell_points =
GeometryInfo<dim>::faces_per_cell * n_face_points;
this->reinit_restriction_and_prolongation_matrices(true, true);
FETools::compute_embedding_matrices(*this, this->prolongation, true, 1.);
+ AssertDimension(this->n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
FullMatrix<double> face_embeddings[GeometryInfo<dim>::max_children_per_face];
for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
- face_embeddings[i].reinit(this->n_dofs_per_face(), this->n_dofs_per_face());
+ face_embeddings[i].reinit(this->n_dofs_per_face(face_no),
+ this->n_dofs_per_face(face_no));
FETools::compute_face_embedding_matrices(*this, face_embeddings, 0, 0, 1.);
- this->interface_constraints.reinit((1 << (dim - 1)) * this->n_dofs_per_face(),
- this->n_dofs_per_face());
+ this->interface_constraints.reinit((1 << (dim - 1)) *
+ this->n_dofs_per_face(face_no),
+ this->n_dofs_per_face(face_no));
unsigned int target_row = 0;
for (unsigned int d = 0; d < GeometryInfo<dim>::max_children_per_face; ++d)
for (unsigned int i = 0; i < face_embeddings[d].m(); ++i)
// initialize_support_points()
if (test_values_face.size() == 0)
{
- for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i)
+ for (unsigned int i = 0; i < this->n_dofs_per_face(f); ++i)
nodal_values[dbase + i] =
support_point_values[pbase + i]
[GeometryInfo<dim>::unit_normal_direction[f]];
- pbase += this->n_dofs_per_face();
+ pbase += this->n_dofs_per_face(f);
}
else
{
- for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i)
+ for (unsigned int i = 0; i < this->n_dofs_per_face(f); ++i)
{
double s = 0.;
for (unsigned int k = 0; k < test_values_face.size(); ++k)
}
pbase += test_values_face.size();
}
- dbase += this->n_dofs_per_face();
+ dbase += this->n_dofs_per_face(f);
}
+ AssertDimension(this->n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+ (void)face_no;
+
AssertDimension(dbase,
- this->n_dofs_per_face() * GeometryInfo<dim>::faces_per_cell);
+ this->n_dofs_per_face(face_no) *
+ GeometryInfo<dim>::faces_per_cell);
AssertDimension(pbase,
this->generalized_support_points.size() -
test_values_cell.size());
// considered later. In 2D, we can use point values.
QGauss<dim - 1> face_points(deg + 1);
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(this->n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
// Copy the quadrature formula to the face points.
- this->generalized_face_support_points[0].resize(face_points.size());
+ this->generalized_face_support_points[face_no].resize(face_points.size());
for (unsigned int k = 0; k < face_points.size(); ++k)
- this->generalized_face_support_points[0][k] = face_points.point(k);
+ this->generalized_face_support_points[face_no][k] = face_points.point(k);
// In the interior, we only test with polynomials of degree up to
// deg-2, thus we use deg points. Note that deg>=1 and the lowest
k < face_points.size() * GeometryInfo<dim>::faces_per_cell;
++k)
this->generalized_support_points[k] = faces.point(
- k + QProjector<dim>::DataSetDescriptor::face(this->reference_cell_type(),
- 0,
- true,
- false,
- false,
- this->n_dofs_per_face()));
+ k +
+ QProjector<dim>::DataSetDescriptor::face(this->reference_cell_type(),
+ 0,
+ true,
+ false,
+ false,
+ this->n_dofs_per_face(face_no)));
// Currently, for backward compatibility, we do not use moments, but
// point values on faces in 2D. In 3D, this is impossible, since the
const FiniteElement<dim, spacedim> &x_source_fe,
const unsigned int subface,
FullMatrix<double> & interpolation_matrix,
- const unsigned int) const
+ const unsigned int face_no) const
{
- Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(),
+ Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
ExcDimensionMismatch(interpolation_matrix.m(),
- x_source_fe.n_dofs_per_face()));
+ x_source_fe.n_dofs_per_face(face_no)));
// see if source is a Bernstein element
if (const FE_Bernstein<dim, spacedim> *source_fe =
{
// have this test in here since a table of size 2x0 reports its size as
// 0x0
- Assert(interpolation_matrix.n() == this->n_dofs_per_face(),
+ Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
ExcDimensionMismatch(interpolation_matrix.n(),
- this->n_dofs_per_face()));
+ this->n_dofs_per_face(face_no)));
// Make sure that the element for which the DoFs should be constrained
// is the one with the higher polynomial degree. Actually the procedure
// produced in that case might lead to problems in the hp procedures,
// which use this method.
Assert(
- this->n_dofs_per_face() <= source_fe->n_dofs_per_face(),
+ this->n_dofs_per_face(face_no) <= source_fe->n_dofs_per_face(face_no),
(typename FiniteElement<dim,
spacedim>::ExcInterpolationNotImplemented()));
const Quadrature<dim - 1> quad_face_support(
FE_Q<dim, spacedim>(QIterated<1>(QTrapez<1>(), source_fe->degree))
- .get_unit_face_support_points());
+ .get_unit_face_support_points(face_no));
// Rule of thumb for FP accuracy, that can be expected for a given
// polynomial degree. This value is used to cut off values close to
0,
subface);
- for (unsigned int i = 0; i < source_fe->n_dofs_per_face(); ++i)
+ for (unsigned int i = 0; i < source_fe->n_dofs_per_face(face_no); ++i)
{
const Point<dim> &p = subface_quadrature.point(i);
- for (unsigned int j = 0; j < this->n_dofs_per_face(); ++j)
+ for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
{
double matrix_entry =
this->shape_value(this->face_to_cell_index(j, 0), p);
// make sure that the row sum of each of the matrices is 1 at this
// point. this must be so since the shape functions sum up to 1
- for (unsigned int j = 0; j < source_fe->n_dofs_per_face(); ++j)
+ for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j)
{
double sum = 0.;
- for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i)
+ for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
sum += interpolation_matrix(j, i);
Assert(std::fabs(sum - 1) < eps, ExcInternalError());
// equivalencies to be recorded
return std::vector<std::pair<unsigned int, unsigned int>>();
}
- else if (fe_other.n_dofs_per_face() == 0)
+ else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
{
// if the other element has no elements on faces at all,
// then it would be impossible to enforce any kind of
// If the system is not primitive, these have not been initialized by
// FiniteElement
this->system_to_component_table.resize(this->n_dofs_per_cell());
- this->face_system_to_component_table[0].resize(this->n_dofs_per_face());
FETools::Compositing::build_cell_tables(this->system_to_base_table,
this->system_to_component_table,
*this,
false);
- FETools::Compositing::build_face_tables(
- this->face_system_to_base_table[0],
- this->face_system_to_component_table[0],
- *this,
- false);
+ this->face_system_to_component_table.resize(this->n_unique_faces());
+
+ for (unsigned int face_no = 0; face_no < this->n_unique_faces(); ++face_no)
+ {
+ this->face_system_to_component_table[0].resize(
+ this->n_dofs_per_face(face_no));
+
+
+ FETools::Compositing::build_face_tables(
+ this->face_system_to_base_table[face_no],
+ this->face_system_to_component_table[face_no],
+ *this,
+ false,
+ face_no);
+ }
}
// restriction and prolongation matrices are built on demand
// this FE sits on the boundary or not. Thus for moment just copy support
// points from fe system:
{
- this->unit_support_points = fe_system->unit_support_points;
- this->unit_face_support_points[0] = fe_system->unit_face_support_points[0];
+ this->unit_support_points = fe_system->unit_support_points;
+ this->unit_face_support_points = fe_system->unit_face_support_points;
}
// take adjust_quad_dof_index_for_face_orientation_table from FESystem:
const FiniteElement<dim, spacedim> &x_source_fe,
const unsigned int subface,
FullMatrix<double> & interpolation_matrix,
- const unsigned int) const
+ const unsigned int face_no) const
{
// this function is similar to the respective method in FE_Q
- Assert(interpolation_matrix.n() == this->n_dofs_per_face(),
+ Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
ExcDimensionMismatch(interpolation_matrix.n(),
- this->n_dofs_per_face()));
- Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(),
+ this->n_dofs_per_face(face_no)));
+ Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
ExcDimensionMismatch(interpolation_matrix.m(),
- x_source_fe.n_dofs_per_face()));
+ x_source_fe.n_dofs_per_face(face_no)));
// see if source is a FaceQ element
if (const FE_FaceQ<dim, spacedim> *source_fe =
// produced in that case might lead to problems in the hp procedures,
// which use this method.
Assert(
- this->n_dofs_per_face() <= source_fe->n_dofs_per_face(),
+ this->n_dofs_per_face(face_no) <= source_fe->n_dofs_per_face(face_no),
(typename FiniteElement<dim,
spacedim>::ExcInterpolationNotImplemented()));
// generate a quadrature with the unit face support points.
const Quadrature<dim - 1> face_quadrature(
- source_fe->get_unit_face_support_points());
+ source_fe->get_unit_face_support_points(face_no));
// Rule of thumb for FP accuracy, that can be expected for a given
// polynomial degree. This value is used to cut off values close to
// compute the interpolation matrix by simply taking the value at the
// support points.
- for (unsigned int i = 0; i < source_fe->n_dofs_per_face(); ++i)
+ for (unsigned int i = 0; i < source_fe->n_dofs_per_face(face_no); ++i)
{
const Point<dim - 1> p =
subface == numbers::invalid_unsigned_int ?
GeometryInfo<dim - 1>::child_to_cell_coordinates(
face_quadrature.point(i), subface);
- for (unsigned int j = 0; j < this->n_dofs_per_face(); ++j)
+ for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
{
double matrix_entry = this->poly_space.compute_value(j, p);
// make sure that the row sum of each of the matrices is 1 at this
// point. this must be so since the shape functions sum up to 1
- for (unsigned int j = 0; j < source_fe->n_dofs_per_face(); ++j)
+ for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j)
{
double sum = 0.;
- for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i)
+ for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
sum += interpolation_matrix(j, i);
Assert(std::fabs(sum - 1) < eps, ExcInternalError());
const unsigned int shape_index,
const unsigned int face_index) const
{
- return (face_index == (shape_index / this->n_dofs_per_face()));
+ return (face_index == (shape_index / this->n_dofs_per_face(face_index)));
}
// equivalencies to be recorded
return std::vector<std::pair<unsigned int, unsigned int>>();
}
- else if (fe_other.n_dofs_per_face() == 0)
+ else if (fe_other.n_unique_faces() == 1 &&
+ fe_other.n_dofs_per_face(0) == 0)
{
// if the other element has no elements on faces at all,
// then it would be impossible to enforce any kind of
// equivalencies to be recorded
return std::vector<std::pair<unsigned int, unsigned int>>();
}
- else if (fe_other.n_dofs_per_face() == 0)
+ else if (fe_other.n_unique_faces() == 1 &&
+ fe_other.n_dofs_per_face(0) == 0)
{
// if the other element has no elements on faces at all,
// then it would be impossible to enforce any kind of
const FiniteElement<1, spacedim> &x_source_fe,
const unsigned int /*subface*/,
FullMatrix<double> &interpolation_matrix,
- const unsigned int) const
+ const unsigned int face_no) const
{
(void)x_source_fe;
- Assert(interpolation_matrix.n() == this->n_dofs_per_face(),
+ (void)face_no;
+
+ Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
ExcDimensionMismatch(interpolation_matrix.n(),
- this->n_dofs_per_face()));
- Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(),
+ this->n_dofs_per_face(face_no)));
+ Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
ExcDimensionMismatch(interpolation_matrix.m(),
- x_source_fe.n_dofs_per_face()));
+ x_source_fe.n_dofs_per_face(face_no)));
interpolation_matrix(0, 0) = 1.;
}
const unsigned int shape_index,
const unsigned int face_index) const
{
- return (face_index == (shape_index / this->n_dofs_per_face()));
+ return (face_index == (shape_index / this->n_dofs_per_face(face_index)));
}
const FiniteElement<dim, spacedim> &x_source_fe,
const unsigned int subface,
FullMatrix<double> & interpolation_matrix,
- const unsigned int) const
+ const unsigned int face_no) const
{
// this function is similar to the respective method in FE_Q
- Assert(interpolation_matrix.n() == this->n_dofs_per_face(),
+ Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
ExcDimensionMismatch(interpolation_matrix.n(),
- this->n_dofs_per_face()));
- Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(),
+ this->n_dofs_per_face(face_no)));
+ Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
ExcDimensionMismatch(interpolation_matrix.m(),
- x_source_fe.n_dofs_per_face()));
+ x_source_fe.n_dofs_per_face(face_no)));
// see if source is a FaceP element
if (const FE_FaceP<dim, spacedim> *source_fe =
// produced in that case might lead to problems in the hp procedures,
// which use this method.
Assert(
- this->n_dofs_per_face() <= source_fe->n_dofs_per_face(),
+ this->n_dofs_per_face(face_no) <= source_fe->n_dofs_per_face(face_no),
(typename FiniteElement<dim,
spacedim>::ExcInterpolationNotImplemented()));
const double eps = 2e-13 * (this->degree + 1) * (dim - 1);
FullMatrix<double> mass(face_quadrature.size(),
- source_fe->n_dofs_per_face());
+ source_fe->n_dofs_per_face(face_no));
for (unsigned int k = 0; k < face_quadrature.size(); ++k)
{
GeometryInfo<dim - 1>::child_to_cell_coordinates(
face_quadrature.point(k), subface);
- for (unsigned int j = 0; j < source_fe->n_dofs_per_face(); ++j)
+ for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j)
mass(k, j) = source_fe->poly_space.compute_value(j, p);
}
Householder<double> H(mass);
Vector<double> v_in(face_quadrature.size());
- Vector<double> v_out(source_fe->n_dofs_per_face());
+ Vector<double> v_out(source_fe->n_dofs_per_face(face_no));
// compute the interpolation matrix by evaluating on the fine side and
// then solving the least squares problem
- for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i)
+ for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
{
for (unsigned int k = 0; k < face_quadrature.size(); ++k)
{
(void)result;
Assert(result < 1e-12, FETools::ExcLeastSquaresError(result));
- for (unsigned int j = 0; j < source_fe->n_dofs_per_face(); ++j)
+ for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j)
{
double matrix_entry = v_out(j);
{
Table<2, bool> constant_modes(1, this->n_dofs_per_cell());
for (unsigned int face : GeometryInfo<dim>::face_indices())
- constant_modes(0, face * this->n_dofs_per_face()) = true;
+ constant_modes(0, face * this->n_dofs_per_face(face)) = true;
return std::pair<Table<2, bool>, std::vector<unsigned int>>(
constant_modes, std::vector<unsigned int>(1, 0));
}
#endif
FullMatrix<double> face_embeddings[GeometryInfo<dim>::max_children_per_face];
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(this->n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
- face_embeddings[i].reinit(this->n_dofs_per_face(), this->n_dofs_per_face());
+ face_embeddings[i].reinit(this->n_dofs_per_face(face_no),
+ this->n_dofs_per_face(face_no));
FETools::compute_face_embedding_matrices<dim, double>(
*this,
case 2:
{
- this->interface_constraints.reinit(2 * this->n_dofs_per_face(),
- this->n_dofs_per_face());
+ this->interface_constraints.reinit(2 * this->n_dofs_per_face(face_no),
+ this->n_dofs_per_face(face_no));
for (unsigned int i = 0; i < GeometryInfo<2>::max_children_per_face;
++i)
- for (unsigned int j = 0; j < this->n_dofs_per_face(); ++j)
- for (unsigned int k = 0; k < this->n_dofs_per_face(); ++k)
- this->interface_constraints(i * this->n_dofs_per_face() + j,
+ for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
+ for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k)
+ this->interface_constraints(i * this->n_dofs_per_face(face_no) +
+ j,
k) = face_embeddings[i](j, k);
break;
case 3:
{
- this->interface_constraints.reinit(4 * (this->n_dofs_per_face() -
- this->degree),
- this->n_dofs_per_face());
+ this->interface_constraints.reinit(
+ 4 * (this->n_dofs_per_face(face_no) - this->degree),
+ this->n_dofs_per_face(face_no));
unsigned int target_row = 0;
for (unsigned int i = 0; i < 2; ++i)
for (unsigned int j = this->degree; j < 2 * this->degree;
++j, ++target_row)
- for (unsigned int k = 0; k < this->n_dofs_per_face(); ++k)
+ for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k)
this->interface_constraints(target_row, k) =
face_embeddings[2 * i](j, k);
for (unsigned int j = 3 * this->degree;
j < GeometryInfo<3>::lines_per_face * this->degree;
++j, ++target_row)
- for (unsigned int k = 0; k < this->n_dofs_per_face(); ++k)
+ for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k)
this->interface_constraints(target_row, k) =
face_embeddings[i](j, k);
for (unsigned int k = i * this->degree;
k < (i + 1) * this->degree;
++k, ++target_row)
- for (unsigned int l = 0; l < this->n_dofs_per_face(); ++l)
+ for (unsigned int l = 0; l < this->n_dofs_per_face(face_no);
+ ++l)
this->interface_constraints(target_row, l) =
face_embeddings[i + 2 * j](k, l);
for (unsigned int k = (i + 2) * this->degree;
k < (i + 3) * this->degree;
++k, ++target_row)
- for (unsigned int l = 0; l < this->n_dofs_per_face(); ++l)
+ for (unsigned int l = 0; l < this->n_dofs_per_face(face_no);
+ ++l)
this->interface_constraints(target_row, l) =
face_embeddings[2 * i + j](k, l);
++i)
for (unsigned int j =
GeometryInfo<3>::lines_per_face * this->degree;
- j < this->n_dofs_per_face();
+ j < this->n_dofs_per_face(face_no);
++j, ++target_row)
- for (unsigned int k = 0; k < this->n_dofs_per_face(); ++k)
+ for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k)
this->interface_constraints(target_row, k) =
face_embeddings[i](j, k);
{
const int dim = 2;
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(this->n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
// Create polynomial basis.
const std::vector<Polynomials::Polynomial<double>> &lobatto_polynomials =
Polynomials::Lobatto::generate_complete_basis(order + 1);
QProjector<dim>::project_to_all_faces(this->reference_cell_type(),
reference_edge_quadrature);
- this->generalized_face_support_points[0].resize(n_edge_points);
+ this->generalized_face_support_points[face_no].resize(n_edge_points);
// Create face support points.
for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
- this->generalized_face_support_points[0][q_point] =
+ this->generalized_face_support_points[face_no][q_point] =
reference_edge_quadrature.point(q_point);
if (order > 0)
boundary_weights(q_point, i) =
reference_edge_quadrature.weight(q_point) *
lobatto_polynomials_grad[i + 1].value(
- this->generalized_face_support_points[0][q_point](0));
+ this->generalized_face_support_points[face_no][q_point](0));
}
for (unsigned int q_point = 0; q_point < n_interior_points; ++q_point)
{
const int dim = 3;
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(this->n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
// Create polynomial basis.
const std::vector<Polynomials::Polynomial<double>> &lobatto_polynomials =
Polynomials::Lobatto::generate_complete_basis(order + 1);
boundary_weights.reinit(n_edge_points + n_face_points,
2 * (order + 1) * order);
- this->generalized_face_support_points[0].resize(4 * n_edge_points +
- n_face_points);
+ this->generalized_face_support_points[face_no].resize(4 * n_edge_points +
+ n_face_points);
this->generalized_support_points.resize(n_boundary_points +
n_interior_points);
for (unsigned int line = 0;
line < GeometryInfo<dim - 1>::lines_per_cell;
++line)
- this->generalized_face_support_points[0][line * n_edge_points +
- q_point] =
+ this
+ ->generalized_face_support_points[face_no][line * n_edge_points +
+ q_point] =
edge_quadrature.point(
QProjector<dim - 1>::DataSetDescriptor::face(
ReferenceCell::get_hypercube(dim - 1),
boundary_weights(q_point, i) =
reference_edge_quadrature.weight(q_point) *
lobatto_polynomials_grad[i + 1].value(
- this->generalized_face_support_points[0][q_point](1));
+ this->generalized_face_support_points[face_no][q_point](1));
}
// Create support points on faces.
for (unsigned int q_point = 0; q_point < n_face_points; ++q_point)
{
- this
- ->generalized_face_support_points[0][q_point + 4 * n_edge_points] =
+ this->generalized_face_support_points[face_no]
+ [q_point + 4 * n_edge_points] =
reference_face_quadrature.point(q_point);
for (unsigned int i = 0; i <= order; ++i)
boundary_weights(q_point + n_edge_points, 2 * (i * order + j)) =
reference_face_quadrature.weight(q_point) *
lobatto_polynomials_grad[i].value(
- this->generalized_face_support_points[0][q_point +
- 4 * n_edge_points](
- 0)) *
+ this->generalized_face_support_points
+ [face_no][q_point + 4 * n_edge_points](0)) *
lobatto_polynomials[j + 2].value(
- this->generalized_face_support_points[0][q_point +
- 4 * n_edge_points](
- 1));
+ this->generalized_face_support_points
+ [face_no][q_point + 4 * n_edge_points](1));
boundary_weights(q_point + n_edge_points,
2 * (i * order + j) + 1) =
reference_face_quadrature.weight(q_point) *
lobatto_polynomials_grad[i].value(
- this->generalized_face_support_points[0][q_point +
- 4 * n_edge_points](
- 1)) *
+ this->generalized_face_support_points
+ [face_no][q_point + 4 * n_edge_points](1)) *
lobatto_polynomials[j + 2].value(
- this->generalized_face_support_points[0][q_point +
- 4 * n_edge_points](
- 0));
+ this->generalized_face_support_points
+ [face_no][q_point + 4 * n_edge_points](0));
}
}
else
{
- this->generalized_face_support_points[0].resize(4 * n_edge_points);
+ this->generalized_face_support_points[face_no].resize(4 * n_edge_points);
this->generalized_support_points.resize(
GeometryInfo<dim>::lines_per_cell * n_edge_points);
for (unsigned int line = 0;
line < GeometryInfo<dim - 1>::lines_per_cell;
++line)
- this->generalized_face_support_points[0][line * n_edge_points +
- q_point] =
+ this
+ ->generalized_face_support_points[face_no][line * n_edge_points +
+ q_point] =
edge_quadrature.point(
QProjector<dim - 1>::DataSetDescriptor::face(
ReferenceCell::get_hypercube(dim - 1),
FE_Nedelec<dim>::get_face_interpolation_matrix(
const FiniteElement<dim> &source,
FullMatrix<double> & interpolation_matrix,
- const unsigned int) const
+ const unsigned int face_no) const
{
+ (void)face_no;
// this is only implemented, if the
// source FE is also a
// Nedelec element
AssertThrow((source.get_name().find("FE_Nedelec<") == 0) ||
(dynamic_cast<const FE_Nedelec<dim> *>(&source) != nullptr),
(typename FiniteElement<dim>::ExcInterpolationNotImplemented()));
- Assert(interpolation_matrix.m() == source.n_dofs_per_face(),
+ Assert(interpolation_matrix.m() == source.n_dofs_per_face(face_no),
ExcDimensionMismatch(interpolation_matrix.m(),
- source.n_dofs_per_face()));
- Assert(interpolation_matrix.n() == this->n_dofs_per_face(),
+ source.n_dofs_per_face(face_no)));
+ Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
ExcDimensionMismatch(interpolation_matrix.n(),
- this->n_dofs_per_face()));
+ this->n_dofs_per_face(face_no)));
// ok, source is a Nedelec element, so
// we will be able to do the work
// lead to problems in the
// hp procedures, which use this
// method.
- Assert(this->n_dofs_per_face() <= source_fe.n_dofs_per_face(),
+ Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
(typename FiniteElement<dim>::ExcInterpolationNotImplemented()));
interpolation_matrix = 0;
const FiniteElement<dim> &source,
const unsigned int subface,
FullMatrix<double> & interpolation_matrix,
- const unsigned int) const
+ const unsigned int face_no) const
{
// this is only implemented, if the
// source FE is also a
AssertThrow((source.get_name().find("FE_Nedelec<") == 0) ||
(dynamic_cast<const FE_Nedelec<dim> *>(&source) != nullptr),
typename FiniteElement<dim>::ExcInterpolationNotImplemented());
- Assert(interpolation_matrix.m() == source.n_dofs_per_face(),
+ Assert(interpolation_matrix.m() == source.n_dofs_per_face(face_no),
ExcDimensionMismatch(interpolation_matrix.m(),
- source.n_dofs_per_face()));
- Assert(interpolation_matrix.n() == this->n_dofs_per_face(),
+ source.n_dofs_per_face(face_no)));
+ Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
ExcDimensionMismatch(interpolation_matrix.n(),
- this->n_dofs_per_face()));
+ this->n_dofs_per_face(face_no)));
// ok, source is a Nedelec element, so
// we will be able to do the work
// lead to problems in the
// hp procedures, which use this
// method.
- Assert(this->n_dofs_per_face() <= source_fe.n_dofs_per_face(),
+ Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
(typename FiniteElement<dim>::ExcInterpolationNotImplemented()));
interpolation_matrix = 0.0;
// Perform projection-based interpolation
{
case 2:
{
- for (unsigned int dof = 0; dof < this->n_dofs_per_face(); ++dof)
+ for (unsigned int dof = 0; dof < this->n_dofs_per_face(face_no);
+ ++dof)
for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
++q_point)
{
Vector<double> solution(source_fe.degree - 1);
Vector<double> system_rhs(source_fe.degree - 1);
- for (unsigned int dof = 0; dof < this->n_dofs_per_face(); ++dof)
+ for (unsigned int dof = 0; dof < this->n_dofs_per_face(face_no);
+ ++dof)
{
system_rhs = 0.0;
{0.0, 1.0},
{1.0, 1.0}};
- for (unsigned int dof = 0; dof < this->n_dofs_per_face(); ++dof)
+ for (unsigned int dof = 0; dof < this->n_dofs_per_face(face_no);
+ ++dof)
for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
++q_point)
{
GeometryInfo<dim>::lines_per_face);
Vector<double> tmp(GeometryInfo<dim>::lines_per_face);
- for (unsigned int dof = 0; dof < this->n_dofs_per_face(); ++dof)
+ for (unsigned int dof = 0; dof < this->n_dofs_per_face(face_no);
+ ++dof)
{
system_rhs = 0.0;
system_rhs.reinit(system_matrix_inv.m(), 2);
tmp.reinit(2);
- for (unsigned int dof = 0; dof < this->n_dofs_per_face(); ++dof)
+ for (unsigned int dof = 0; dof < this->n_dofs_per_face(face_no);
+ ++dof)
{
system_rhs = 0.0;
const std::vector<Vector<double>> &support_point_values,
std::vector<double> & nodal_values) const
{
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(this->n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
const unsigned int deg = this->degree - 1;
Assert(support_point_values.size() == this->generalized_support_points.size(),
ExcDimensionMismatch(support_point_values.size(),
system_matrix(i, j) +=
boundary_weights(q_point, j) *
lobatto_polynomials_grad[i + 1].value(
- this->generalized_face_support_points[0][q_point](0));
+ this->generalized_face_support_points[face_no][q_point](
+ 0));
FullMatrix<double> system_matrix_inv(this->degree - 1,
this->degree - 1);
system_matrix(i, j) +=
boundary_weights(q_point, j) *
lobatto_polynomials_grad[i + 1].value(
- this->generalized_face_support_points[0][q_point](1));
+ this->generalized_face_support_points[face_no][q_point](
+ 1));
FullMatrix<double> system_matrix_inv(this->degree - 1,
this->degree - 1);
2 * (k * (this->degree - 1) + l)) *
legendre_polynomials[i].value(
this->generalized_face_support_points
- [0][q_point + 4 * n_edge_points](0)) *
+ [face_no][q_point + 4 * n_edge_points](0)) *
lobatto_polynomials[j + 2].value(
this->generalized_face_support_points
- [0][q_point + 4 * n_edge_points](1));
+ [face_no][q_point + 4 * n_edge_points](1));
system_matrix_inv.reinit(system_matrix.m(), system_matrix.m());
system_matrix_inv.invert(system_matrix);
const unsigned int faces_per_cell = GeometryInfo<dim>::faces_per_cell;
const unsigned int n_line_dofs = this->n_dofs_per_line() * lines_per_cell;
- const unsigned int n_face_dofs = this->n_dofs_per_quad() * faces_per_cell;
+
+ // we assume that all quads have the same numer of dofs
+ const unsigned int n_face_dofs = this->n_dofs_per_quad(0) * faces_per_cell;
const UpdateFlags flags(data.update_each);
const unsigned int n_q_points = quadrature.size();
// Loop through quad points:
for (unsigned int m = 0; m < faces_per_cell; ++m)
{
- const unsigned int shift_m(m * this->n_dofs_per_quad());
+ // we assume that all quads have the same numer of dofs
+ const unsigned int shift_m(m * this->n_dofs_per_quad(0));
// Calculate the offsets for each face-based shape function:
//
// Type-1 (gradients)
const unsigned int nn = cell->neighbor_face_no(f);
if (nn < GeometryInfo<dim>::faces_per_cell / 2)
- for (unsigned int j = 0; j < fe.n_dofs_per_face(); ++j)
+ for (unsigned int j = 0; j < fe.n_dofs_per_face(f); ++j)
{
const unsigned int cell_j = fe.face_to_cell_index(j, f);
- Assert(f * fe.n_dofs_per_face() + j < face_sign.size(),
+ Assert(f * fe.n_dofs_per_face(f) + j < face_sign.size(),
ExcInternalError());
Assert(mapping_kind.size() == 1 ||
cell_j < mapping_kind.size(),
if ((mapping_kind.size() > 1 ?
mapping_kind[cell_j] :
mapping_kind[0]) == mapping_raviart_thomas)
- face_sign[f * fe.n_dofs_per_face() + j] = -1.0;
+ face_sign[f * fe.n_dofs_per_face(f) + j] = -1.0;
}
}
}
const FiniteElement<dim, spacedim> &x_source_fe,
const unsigned int subface,
FullMatrix<double> & interpolation_matrix,
- const unsigned int) const
+ const unsigned int face_no) const
{
- Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(),
+ Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
ExcDimensionMismatch(interpolation_matrix.m(),
- x_source_fe.n_dofs_per_face()));
+ x_source_fe.n_dofs_per_face(face_no)));
// see if source is a Q element
if (const FE_Q_Base<PolynomialType, dim, spacedim> *source_fe =
{
// have this test in here since a table of size 2x0 reports its size as
// 0x0
- Assert(interpolation_matrix.n() == this->n_dofs_per_face(),
+ Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
ExcDimensionMismatch(interpolation_matrix.n(),
- this->n_dofs_per_face()));
+ this->n_dofs_per_face(face_no)));
// Make sure that the element for which the DoFs should be constrained
// is the one with the higher polynomial degree. Actually the procedure
// produced in that case might lead to problems in the hp procedures,
// which use this method.
Assert(
- this->n_dofs_per_face() <= source_fe->n_dofs_per_face(),
+ this->n_dofs_per_face(face_no) <= source_fe->n_dofs_per_face(face_no),
(typename FiniteElement<dim,
spacedim>::ExcInterpolationNotImplemented()));
// generate a point on this cell and evaluate the shape functions there
const Quadrature<dim - 1> quad_face_support(
- source_fe->get_unit_face_support_points());
+ source_fe->get_unit_face_support_points(face_no));
// Rule of thumb for FP accuracy, that can be expected for a given
// polynomial degree. This value is used to cut off values close to
quad_face_support,
0,
subface);
- for (unsigned int i = 0; i < source_fe->n_dofs_per_face(); ++i)
+ for (unsigned int i = 0; i < source_fe->n_dofs_per_face(face_no); ++i)
{
const Point<dim> &p = subface_quadrature.point(i);
- for (unsigned int j = 0; j < this->n_dofs_per_face(); ++j)
+ for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
{
double matrix_entry =
this->shape_value(this->face_to_cell_index(j, 0), p);
// make sure that the row sum of each of the matrices is 1 at this
// point. this must be so since the shape functions sum up to 1
- for (unsigned int j = 0; j < source_fe->n_dofs_per_face(); ++j)
+ for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j)
{
double sum = 0.;
- for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i)
+ for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
sum += interpolation_matrix(j, i);
Assert(std::fabs(sum - 1) < eps, ExcInternalError());
// equivalencies to be recorded
return std::vector<std::pair<unsigned int, unsigned int>>();
}
- else if (fe_other.n_dofs_per_face() == 0)
+ else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
{
// if the other element has no elements on faces at all,
// then it would be impossible to enforce any kind of
// equivalencies to be recorded
return std::vector<std::pair<unsigned int, unsigned int>>();
}
- else if (fe_other.n_dofs_per_face() == 0)
+ else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
{
// if the other element has no elements on faces at all,
// then it would be impossible to enforce any kind of
// equivalencies to be recorded
return std::vector<std::pair<unsigned int, unsigned int>>();
}
- else if (fe_other.n_dofs_per_face() == 0)
+ else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
{
// if the other element has no elements on faces at all,
// then it would be impossible to enforce any kind of
if (dim == 1)
return;
- this->unit_face_support_points[0].resize(
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(this->n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
+ this->unit_face_support_points[face_no].resize(
Utilities::fixed_power<dim - 1>(q_degree + 1));
// find renumbering of faces and assign from values of quadrature
// The only thing we have to do is reorder the points from tensor
// product order to the order in which we enumerate DoFs on cells
- this->unit_face_support_points[0].resize(support_quadrature.size());
+ this->unit_face_support_points[face_no].resize(support_quadrature.size());
for (unsigned int k = 0; k < support_quadrature.size(); ++k)
- this->unit_face_support_points[0][face_index_map[k]] =
+ this->unit_face_support_points[face_no][face_index_map[k]] =
support_quadrature.point(k);
}
if (dim < 3)
return;
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(this->n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
Assert(this->adjust_quad_dof_index_for_face_orientation_table[0]
- .n_elements() == 8 * this->n_dofs_per_quad(),
+ .n_elements() == 8 * this->n_dofs_per_quad(face_no),
ExcInternalError());
const unsigned int n = q_degree - 1;
- Assert(n * n == this->n_dofs_per_quad(), ExcInternalError());
+ Assert(n * n == this->n_dofs_per_quad(face_no), ExcInternalError());
// the dofs on a face are connected to a n x n matrix. for example, for
// degree==4 we have the following dofs on a quad
// rotated and mirrored numbers.
- for (unsigned int local = 0; local < this->n_dofs_per_quad(); ++local)
+ for (unsigned int local = 0; local < this->n_dofs_per_quad(face_no); ++local)
// face support points are in lexicographic ordering with x running
// fastest. invert that (y running fastest)
{
unsigned int i = local % n, j = local / n;
// face_orientation=false, face_flip=false, face_rotation=false
- this->adjust_quad_dof_index_for_face_orientation_table[0](local, 0) =
+ this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
+ 0) =
j + i * n - local;
// face_orientation=false, face_flip=false, face_rotation=true
- this->adjust_quad_dof_index_for_face_orientation_table[0](local, 1) =
+ this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
+ 1) =
i + (n - 1 - j) * n - local;
// face_orientation=false, face_flip=true, face_rotation=false
- this->adjust_quad_dof_index_for_face_orientation_table[0](local, 2) =
+ this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
+ 2) =
(n - 1 - j) + (n - 1 - i) * n - local;
// face_orientation=false, face_flip=true, face_rotation=true
- this->adjust_quad_dof_index_for_face_orientation_table[0](local, 3) =
+ this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
+ 3) =
(n - 1 - i) + j * n - local;
// face_orientation=true, face_flip=false, face_rotation=false
- this->adjust_quad_dof_index_for_face_orientation_table[0](local, 4) = 0;
+ this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
+ 4) = 0;
// face_orientation=true, face_flip=false, face_rotation=true
- this->adjust_quad_dof_index_for_face_orientation_table[0](local, 5) =
+ this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
+ 5) =
j + (n - 1 - i) * n - local;
// face_orientation=true, face_flip=true, face_rotation=false
- this->adjust_quad_dof_index_for_face_orientation_table[0](local, 6) =
+ this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
+ 6) =
(n - 1 - i) + (n - 1 - j) * n - local;
// face_orientation=true, face_flip=true, face_rotation=true
- this->adjust_quad_dof_index_for_face_orientation_table[0](local, 7) =
+ this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
+ 7) =
(n - 1 - j) + i * n - local;
}
const bool face_flip,
const bool face_rotation) const
{
- AssertIndexRange(face_index, this->n_dofs_per_face());
+ AssertIndexRange(face_index, this->n_dofs_per_face(face));
AssertIndexRange(face, GeometryInfo<dim>::faces_per_cell);
// TODO: we could presumably solve the 3d case below using the
// we need to distinguish between DoFs on vertices, lines and in 3d quads.
// do so in a sequence of if-else statements
- if (face_index < this->get_first_face_line_index())
+ if (face_index < this->get_first_face_line_index(face))
// DoF is on a vertex
{
// get the number of the vertex on the face that corresponds to this DoF,
this->n_dofs_per_vertex() +
dof_index_on_vertex);
}
- else if (face_index < this->get_first_face_quad_index())
+ else if (face_index < this->get_first_face_quad_index(face))
// DoF is on a face
{
// do the same kind of translation as before. we need to only consider
// DoFs on the lines, i.e., ignoring those on the vertices
- const unsigned int index = face_index - this->get_first_face_line_index();
+ const unsigned int index =
+ face_index - this->get_first_face_line_index(face);
const unsigned int face_line = index / this->n_dofs_per_line();
const unsigned int dof_index_on_line = index % this->n_dofs_per_line();
Assert(dim >= 3, ExcInternalError());
// ignore vertex and line dofs
- const unsigned int index = face_index - this->get_first_face_quad_index();
+ const unsigned int index =
+ face_index - this->get_first_face_quad_index(face);
// the same is true here as above for the 3d case -- someone will
// just have to draw a bunch of pictures. in the meantime,
// we can implement the Q2 case in which it is simple
- Assert((this->n_dofs_per_quad() <= 1) ||
+ Assert((this->n_dofs_per_quad(face) <= 1) ||
((face_orientation == true) && (face_flip == false) &&
(face_rotation == false)),
ExcNotImplemented());
- return (this->get_first_quad_index() + face * this->n_dofs_per_quad() +
- index);
+ return (this->get_first_quad_index(face) + index);
}
}
// first, special-case interior shape functions, since they have no support
// no-where on the boundary
- if (((dim == 2) && (shape_index >= this->get_first_quad_index())) ||
+ if (((dim == 2) &&
+ (shape_index >= this->get_first_quad_index(0 /*first quad*/))) ||
((dim == 3) && (shape_index >= this->get_first_hex_index())))
return false;
return false;
}
- else if (shape_index < this->get_first_quad_index())
+ else if (shape_index < this->get_first_quad_index(0 /*first quad*/))
// ok, dof is on a line
{
const unsigned int line_index =
// dof is on a quad
{
const unsigned int quad_index =
- (shape_index - this->get_first_quad_index()) / this->n_dofs_per_quad();
+ (shape_index - this->get_first_quad_index(0)) /
+ this->n_dofs_per_quad(face_index); // this won't work
Assert(static_cast<signed int>(quad_index) <
static_cast<signed int>(GeometryInfo<dim>::quads_per_cell),
ExcInternalError());
// one is an FE_Nothing.
if (dynamic_cast<const FE_Q_Hierarchical<dim> *>(&fe_other) != nullptr)
{
- const unsigned int this_dpq = this->n_dofs_per_quad();
- const unsigned int other_dpq = fe_other.n_dofs_per_quad();
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(this->n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
+ const unsigned int this_dpq = this->n_dofs_per_quad(face_no);
+ const unsigned int other_dpq = fe_other.n_dofs_per_quad(face_no);
// we deal with hierarchical 1d polynomials where dofs are enumerated
// increasingly. Thus we return a vector of pairs for the first N-1, where
FE_Q_Hierarchical<dim>::get_face_interpolation_matrix(
const FiniteElement<dim> &x_source_fe,
FullMatrix<double> & interpolation_matrix,
- const unsigned int) const
+ const unsigned int face_no) const
{
// this is only implemented, if the
// source FE is also a
nullptr),
(typename FiniteElement<dim>::ExcInterpolationNotImplemented()));
- Assert(interpolation_matrix.n() == this->n_dofs_per_face(),
+ Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
ExcDimensionMismatch(interpolation_matrix.n(),
- this->n_dofs_per_face()));
- Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(),
+ this->n_dofs_per_face(face_no)));
+ Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
ExcDimensionMismatch(interpolation_matrix.m(),
- x_source_fe.n_dofs_per_face()));
+ x_source_fe.n_dofs_per_face(face_no)));
// ok, source is a Q_Hierarchical element, so
// we will be able to do the work
// lead to problems in the
// hp procedures, which use this
// method.
- Assert(this->n_dofs_per_face() <= source_fe.n_dofs_per_face(),
+ Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
(typename FiniteElement<dim>::ExcInterpolationNotImplemented()));
interpolation_matrix = 0;
// element, which corresponds to 1 on diagonal of the matrix.
// DoFs which correspond to higher polynomials
// are zeroed (zero rows in the matrix).
- for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i)
+ for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
interpolation_matrix(i, i) = 1;
break;
const FiniteElement<dim> &x_source_fe,
const unsigned int subface,
FullMatrix<double> & interpolation_matrix,
- const unsigned int) const
+ const unsigned int face_no) const
{
// this is only implemented, if the
// source FE is also a
nullptr),
(typename FiniteElement<dim>::ExcInterpolationNotImplemented()));
- Assert(interpolation_matrix.n() == this->n_dofs_per_face(),
+ Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
ExcDimensionMismatch(interpolation_matrix.n(),
- this->n_dofs_per_face()));
- Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(),
+ this->n_dofs_per_face(face_no)));
+ Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
ExcDimensionMismatch(interpolation_matrix.m(),
- x_source_fe.n_dofs_per_face()));
+ x_source_fe.n_dofs_per_face(face_no)));
// ok, source is a Q_Hierarchical element, so
// we will be able to do the work
// lead to problems in the
// hp procedures, which use this
// method.
- Assert(this->n_dofs_per_face() <= source_fe.n_dofs_per_face(),
+ Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
(typename FiniteElement<dim>::ExcInterpolationNotImplemented()));
switch (dim)
interpolation_matrix(1, 0) = 0.5;
interpolation_matrix(1, 1) = 0.5;
- for (unsigned int dof = 2; dof < this->n_dofs_per_face();)
+ for (unsigned int dof = 2;
+ dof < this->n_dofs_per_face(face_no);)
{
interpolation_matrix(1, dof) = -1.0;
dof = dof + 2;
int factorial_i = 1;
- for (unsigned int i = 2; i < this->n_dofs_per_face(); ++i)
+ for (unsigned int i = 2; i < this->n_dofs_per_face(face_no);
+ ++i)
{
interpolation_matrix(i, i) = std::pow(0.5, i);
factorial_i *= i;
int factorial_j = factorial_i;
int factorial_ij = 1;
- for (unsigned int j = i + 1; j < this->n_dofs_per_face();
+ for (unsigned int j = i + 1;
+ j < this->n_dofs_per_face(face_no);
++j)
{
factorial_ij *= j - i;
interpolation_matrix(0, 0) = 0.5;
interpolation_matrix(0, 1) = 0.5;
- for (unsigned int dof = 2; dof < this->n_dofs_per_face();)
+ for (unsigned int dof = 2;
+ dof < this->n_dofs_per_face(face_no);)
{
interpolation_matrix(0, dof) = -1.0;
dof = dof + 2;
int factorial_i = 1;
- for (unsigned int i = 2; i < this->n_dofs_per_face(); ++i)
+ for (unsigned int i = 2; i < this->n_dofs_per_face(face_no);
+ ++i)
{
interpolation_matrix(i, i) = std::pow(0.5, i);
factorial_i *= i;
int factorial_j = factorial_i;
int factorial_ij = 1;
- for (unsigned int j = i + 1; j < this->n_dofs_per_face();
+ for (unsigned int j = i + 1;
+ j < this->n_dofs_per_face(face_no);
++j)
{
factorial_ij *= j - i;
{
const unsigned int codim = dim - 1;
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(this->n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
// number of points: (degree+1)^codim
unsigned int n = this->degree + 1;
for (unsigned int i = 1; i < codim; ++i)
n *= this->degree + 1;
- this->generalized_face_support_points[0].resize(n);
+ this->generalized_face_support_points[face_no].resize(n);
Point<codim> p;
else
p(2) = .5;
}
- this->generalized_face_support_points[0][face_renumber[k++]] = p;
+ this->generalized_face_support_points[face_no][face_renumber[k++]] =
+ p;
}
}
for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
h2l[next_index++] = n + 2 + i;
// inside quad
- Assert(fe.n_dofs_per_quad() ==
+ Assert(fe.n_dofs_per_quad(0 /*only one quad in 2D*/) ==
fe.n_dofs_per_line() * fe.n_dofs_per_line(),
ExcInternalError());
for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
h2l[next_index++] = (2 + i) * n2 + n + 1;
+ // TODO: the implementation makes the assumption that all faces have
+ // the same number of dofs
+ AssertDimension(fe.n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+ (void)face_no;
+
// inside quads
- Assert(fe.n_dofs_per_quad() ==
+ Assert(fe.n_dofs_per_quad(face_no) ==
fe.n_dofs_per_line() * fe.n_dofs_per_line(),
ExcInternalError());
// left face
// inside hex
Assert(fe.n_dofs_per_hex() ==
- fe.n_dofs_per_quad() * fe.n_dofs_per_line(),
+ fe.n_dofs_per_quad(face_no) * fe.n_dofs_per_line(),
ExcInternalError());
for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
for (unsigned int j = 0; j < fe.n_dofs_per_line(); ++j)
// shape functions, since they
// have no support no-where on
// the boundary
- if (((dim == 2) && (shape_index >= this->get_first_quad_index())) ||
+ if (((dim == 2) && (shape_index >=
+ this->get_first_quad_index(0 /*only one quad in 2D*/))) ||
((dim == 3) && (shape_index >= this->get_first_hex_index())))
return false;
return true;
return false;
}
- else if (shape_index < this->get_first_quad_index())
+ else if (shape_index < this->get_first_quad_index(0))
// ok, dof is on a line
{
const unsigned int line_index =
// dof is on a quad
{
const unsigned int quad_index =
- (shape_index - this->get_first_quad_index()) / this->n_dofs_per_quad();
+ (shape_index - this->get_first_quad_index(0 /*first quad*/)) /
+ this->n_dofs_per_quad(face_index);
Assert(static_cast<signed int>(quad_index) <
static_cast<signed int>(GeometryInfo<dim>::quads_per_cell),
ExcInternalError());
FETools::compute_embedding_matrices(*this, this->prolongation);
initialize_restriction();
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(this->n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
// TODO[TL]: for anisotropic refinement we will probably need a table of
// submatrices with an array for each refine case
FullMatrix<double> face_embeddings[GeometryInfo<dim>::max_children_per_face];
for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
- face_embeddings[i].reinit(this->n_dofs_per_face(), this->n_dofs_per_face());
+ face_embeddings[i].reinit(this->n_dofs_per_face(face_no),
+ this->n_dofs_per_face(face_no));
FETools::compute_face_embedding_matrices<dim, double>(*this,
face_embeddings,
0,
0);
- this->interface_constraints.reinit((1 << (dim - 1)) * this->n_dofs_per_face(),
- this->n_dofs_per_face());
+ this->interface_constraints.reinit((1 << (dim - 1)) *
+ this->n_dofs_per_face(face_no),
+ this->n_dofs_per_face(face_no));
unsigned int target_row = 0;
for (unsigned int d = 0; d < GeometryInfo<dim>::max_children_per_face; ++d)
for (unsigned int i = 0; i < face_embeddings[d].m(); ++i)
QGauss<dim> cell_quadrature(deg + 1);
const unsigned int n_interior_points = (deg > 0) ? cell_quadrature.size() : 0;
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(this->n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
unsigned int n_face_points = (dim > 1) ? 1 : 0;
// compute (deg+1)^(dim-1)
for (unsigned int d = 1; d < dim; ++d)
this->generalized_support_points.resize(
GeometryInfo<dim>::faces_per_cell * n_face_points + n_interior_points);
- this->generalized_face_support_points[0].resize(n_face_points);
+ this->generalized_face_support_points[face_no].resize(n_face_points);
// Number of the point being entered
unsigned int current = 0;
boundary_weights.reinit(n_face_points, legendre.n());
- // Assert (face_points.size() == this->n_dofs_per_face(),
- // ExcInternalError());
-
for (unsigned int k = 0; k < n_face_points; ++k)
{
- this->generalized_face_support_points[0][k] = face_points.point(k);
+ this->generalized_face_support_points[face_no][k] =
+ face_points.point(k);
// Compute its quadrature
// contribution for each
// moment.
for (unsigned int k = 0; k < n_face_points; ++k)
for (unsigned int i_child = 0; i_child < this->n_dofs_per_cell();
++i_child)
- for (unsigned int i_face = 0; i_face < this->n_dofs_per_face();
+ for (unsigned int i_face = 0;
+ i_face < this->n_dofs_per_face(face);
++i_face)
{
// The quadrature
// subcell are NOT
// transformed, so we
// have to do it here.
- this->restriction[iso][child](face * this->n_dofs_per_face() +
- i_face,
- i_child) +=
+ this->restriction[iso][child](
+ face * this->n_dofs_per_face(face) + i_face, i_child) +=
Utilities::fixed_power<dim - 1>(.5) * q_sub.weight(k) *
cached_values_on_face(i_child, k) *
this->shape_value_component(
- face * this->n_dofs_per_face() + i_face,
+ face * this->n_dofs_per_face(face) + i_face,
q_sub.point(k),
GeometryInfo<dim>::unit_normal_direction[face]);
}
polynomials[dd] = std::make_unique<AnisotropicPolynomials<dim>>(poly);
}
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(this->n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
QGauss<dim> q_cell(this->degree);
const unsigned int start_cell_dofs =
- GeometryInfo<dim>::faces_per_cell * this->n_dofs_per_face();
+ GeometryInfo<dim>::faces_per_cell * this->n_dofs_per_face(face_no);
// Store shape values, since the
// evaluation suffers if not
for (unsigned int k = 0; k < n_face_points; ++k)
for (unsigned int i = 0; i < boundary_weights.size(1); ++i)
{
- nodal_values[i + face * this->n_dofs_per_face()] +=
+ nodal_values[i + face * this->n_dofs_per_face(face)] +=
boundary_weights(k, i) *
support_point_values[face * n_face_points + k](
GeometryInfo<dim>::unit_normal_direction[face]);
}
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(this->n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
const unsigned int start_cell_dofs =
- GeometryInfo<dim>::faces_per_cell * this->n_dofs_per_face();
+ GeometryInfo<dim>::faces_per_cell * this->n_dofs_per_face(face_no);
const unsigned int start_cell_points =
GeometryInfo<dim>::faces_per_cell * n_face_points;
for (unsigned int i = 0; i < nc; ++i)
this->prolongation[ref_case - 1][i].reinit(n_dofs, n_dofs);
}
+
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(this->n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
// Fill prolongation matrices with embedding operators
FETools::compute_embedding_matrices(*this, this->prolongation);
// TODO[TL]: for anisotropic refinement we will probably need a table of
// submatrices with an array for each refine case
FullMatrix<double> face_embeddings[GeometryInfo<dim>::max_children_per_face];
for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
- face_embeddings[i].reinit(this->n_dofs_per_face(), this->n_dofs_per_face());
+ face_embeddings[i].reinit(this->n_dofs_per_face(face_no),
+ this->n_dofs_per_face(face_no));
FETools::compute_face_embedding_matrices<dim, double>(*this,
face_embeddings,
0,
0);
- this->interface_constraints.reinit((1 << (dim - 1)) * this->n_dofs_per_face(),
- this->n_dofs_per_face());
+ this->interface_constraints.reinit((1 << (dim - 1)) *
+ this->n_dofs_per_face(face_no),
+ this->n_dofs_per_face(face_no));
unsigned int target_row = 0;
for (unsigned int d = 0; d < GeometryInfo<dim>::max_children_per_face; ++d)
for (unsigned int i = 0; i < face_embeddings[d].m(); ++i)
void
FE_RaviartThomasNodal<dim>::initialize_support_points(const unsigned int deg)
{
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(this->n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
this->generalized_support_points.resize(this->n_dofs_per_cell());
- this->generalized_face_support_points[0].resize(this->n_dofs_per_face());
+ this->generalized_face_support_points[face_no].resize(
+ this->n_dofs_per_face(face_no));
// Number of the point being entered
unsigned int current = 0;
if (dim > 1)
{
QGauss<dim - 1> face_points(deg + 1);
- Assert(face_points.size() == this->n_dofs_per_face(), ExcInternalError());
- for (unsigned int k = 0; k < this->n_dofs_per_face(); ++k)
- this->generalized_face_support_points[0][k] = face_points.point(k);
+ Assert(face_points.size() == this->n_dofs_per_face(face_no),
+ ExcInternalError());
+ for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k)
+ this->generalized_face_support_points[face_no][k] =
+ face_points.point(k);
Quadrature<dim> faces =
QProjector<dim>::project_to_all_faces(this->reference_cell_type(),
face_points);
- for (unsigned int k = 0;
- k < this->n_dofs_per_face() * GeometryInfo<dim>::faces_per_cell;
+ for (unsigned int k = 0; k < this->n_dofs_per_face(face_no) *
+ GeometryInfo<dim>::faces_per_cell;
++k)
- this->generalized_support_points[k] = faces.point(
- k +
- QProjector<dim>::DataSetDescriptor::face(this->reference_cell_type(),
- 0,
- true,
- false,
- false,
- this->n_dofs_per_face()));
-
- current = this->n_dofs_per_face() * GeometryInfo<dim>::faces_per_cell;
+ this->generalized_support_points[k] =
+ faces.point(k + QProjector<dim>::DataSetDescriptor::face(
+ this->reference_cell_type(),
+ 0,
+ true,
+ false,
+ false,
+ this->n_dofs_per_face(face_no)));
+
+ current =
+ this->n_dofs_per_face(face_no) * GeometryInfo<dim>::faces_per_cell;
}
if (deg == 0)
unsigned int fbase = 0;
unsigned int f = 0;
for (; f < GeometryInfo<dim>::faces_per_cell;
- ++f, fbase += this->n_dofs_per_face())
+ ++f, fbase += this->n_dofs_per_face(f))
{
- for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i)
+ for (unsigned int i = 0; i < this->n_dofs_per_face(f); ++i)
{
nodal_values[fbase + i] = support_point_values[fbase + i](
GeometryInfo<dim>::unit_normal_direction[f]);
std::vector<std::pair<unsigned int, unsigned int>>
FE_RaviartThomasNodal<dim>::hp_quad_dof_identities(
const FiniteElement<dim> &fe_other,
- const unsigned int) const
+ const unsigned int face_no) const
{
// we can presently only compute
// these identities if both FEs are
// this works exactly like the line
// case above
- const unsigned int p = this->n_dofs_per_quad();
- const unsigned int q = fe_q_other->n_dofs_per_quad();
+ const unsigned int p = this->n_dofs_per_quad(face_no);
+
+ AssertDimension(fe_q_other->n_unique_faces(), 1);
+ const unsigned int q = fe_q_other->n_dofs_per_quad(0);
std::vector<std::pair<unsigned int, unsigned int>> identities;
FE_RaviartThomasNodal<dim>::get_face_interpolation_matrix(
const FiniteElement<dim> &x_source_fe,
FullMatrix<double> & interpolation_matrix,
- const unsigned int) const
+ const unsigned int face_no) const
{
// this is only implemented, if the
// source FE is also a
&x_source_fe) != nullptr),
typename FiniteElement<dim>::ExcInterpolationNotImplemented());
- Assert(interpolation_matrix.n() == this->n_dofs_per_face(),
+ Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
ExcDimensionMismatch(interpolation_matrix.n(),
- this->n_dofs_per_face()));
- Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(),
+ this->n_dofs_per_face(face_no)));
+ Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
ExcDimensionMismatch(interpolation_matrix.m(),
- x_source_fe.n_dofs_per_face()));
+ x_source_fe.n_dofs_per_face(face_no)));
// ok, source is a RaviartThomasNodal element, so
// we will be able to do the work
// lead to problems in the
// hp procedures, which use this
// method.
- Assert(this->n_dofs_per_face() <= source_fe.n_dofs_per_face(),
+ Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
typename FiniteElement<dim>::ExcInterpolationNotImplemented());
// generate a quadrature
// which returns the support
// points on the face.
Quadrature<dim - 1> quad_face_support(
- source_fe.generalized_face_support_points[0]);
+ source_fe.generalized_face_support_points[face_no]);
// Rule of thumb for FP accuracy,
// that can be expected for a
quad_face_support,
0);
- for (unsigned int i = 0; i < source_fe.n_dofs_per_face(); ++i)
+ for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
{
const Point<dim> &p = face_projection.point(i);
- for (unsigned int j = 0; j < this->n_dofs_per_face(); ++j)
+ for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
{
double matrix_entry =
this->shape_value_component(this->face_to_cell_index(j, 0), p, 0);
// this point. this must be so
// since the shape functions sum up
// to 1
- for (unsigned int j = 0; j < source_fe.n_dofs_per_face(); ++j)
+ for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
{
double sum = 0.;
- for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i)
+ for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
sum += interpolation_matrix(j, i);
Assert(std::fabs(sum - 1) < 2e-13 * this->degree * (dim - 1),
const FiniteElement<dim> &x_source_fe,
const unsigned int subface,
FullMatrix<double> & interpolation_matrix,
- const unsigned int) const
+ const unsigned int face_no) const
{
// this is only implemented, if the
// source FE is also a
&x_source_fe) != nullptr),
typename FiniteElement<dim>::ExcInterpolationNotImplemented());
- Assert(interpolation_matrix.n() == this->n_dofs_per_face(),
+ Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
ExcDimensionMismatch(interpolation_matrix.n(),
- this->n_dofs_per_face()));
- Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(),
+ this->n_dofs_per_face(face_no)));
+ Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
ExcDimensionMismatch(interpolation_matrix.m(),
- x_source_fe.n_dofs_per_face()));
+ x_source_fe.n_dofs_per_face(face_no)));
// ok, source is a RaviartThomasNodal element, so
// we will be able to do the work
// lead to problems in the
// hp procedures, which use this
// method.
- Assert(this->n_dofs_per_face() <= source_fe.n_dofs_per_face(),
+ Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
typename FiniteElement<dim>::ExcInterpolationNotImplemented());
// generate a quadrature
// which returns the support
// points on the face.
Quadrature<dim - 1> quad_face_support(
- source_fe.generalized_face_support_points[0]);
+ source_fe.generalized_face_support_points[face_no]);
// Rule of thumb for FP accuracy,
// that can be expected for a
0,
subface);
- for (unsigned int i = 0; i < source_fe.n_dofs_per_face(); ++i)
+ for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
{
const Point<dim> &p = subface_projection.point(i);
- for (unsigned int j = 0; j < this->n_dofs_per_face(); ++j)
+ for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
{
double matrix_entry =
this->shape_value_component(this->face_to_cell_index(j, 0), p, 0);
// this point. this must be so
// since the shape functions sum up
// to 1
- for (unsigned int j = 0; j < source_fe.n_dofs_per_face(); ++j)
+ for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
{
double sum = 0.;
- for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i)
+ for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
sum += interpolation_matrix(j, i);
Assert(std::fabs(sum - 1) < 2e-13 * this->degree * (dim - 1),
for (unsigned int i = 0; i < nc; ++i)
this->prolongation[ref_case - 1][i].reinit(n_dofs, n_dofs);
}
+
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(this->n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
// Fill prolongation matrices with embedding operators
// set tolerance to 1, as embedding error accumulate quickly
FETools::compute_embedding_matrices(*this, this->prolongation, true, 1.0);
FullMatrix<double> face_embeddings[GeometryInfo<dim>::max_children_per_face];
for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
- face_embeddings[i].reinit(this->n_dofs_per_face(), this->n_dofs_per_face());
+ face_embeddings[i].reinit(this->n_dofs_per_face(face_no),
+ this->n_dofs_per_face(face_no));
FETools::compute_face_embedding_matrices<dim, double>(*this,
face_embeddings,
0,
0);
- this->interface_constraints.reinit((1 << (dim - 1)) * this->n_dofs_per_face(),
- this->n_dofs_per_face());
+ this->interface_constraints.reinit((1 << (dim - 1)) *
+ this->n_dofs_per_face(face_no),
+ this->n_dofs_per_face(face_no));
unsigned int target_row = 0;
for (unsigned int d = 0; d < GeometryInfo<dim>::max_children_per_face; ++d)
for (unsigned int i = 0; i < face_embeddings[d].m(); ++i)
void
FE_RT_Bubbles<dim>::initialize_support_points(const unsigned int deg)
{
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(this->n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
this->generalized_support_points.resize(this->n_dofs_per_cell());
- this->generalized_face_support_points[0].resize(this->n_dofs_per_face());
+ this->generalized_face_support_points[face_no].resize(
+ this->n_dofs_per_face(face_no));
// Index of the point being entered
unsigned int current = 0;
if (dim > 1)
{
QGaussLobatto<dim - 1> face_points(deg + 1);
- Assert(face_points.size() == this->n_dofs_per_face(), ExcInternalError());
- for (unsigned int k = 0; k < this->n_dofs_per_face(); ++k)
- this->generalized_face_support_points[0][k] = face_points.point(k);
+ Assert(face_points.size() == this->n_dofs_per_face(face_no),
+ ExcInternalError());
+ for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k)
+ this->generalized_face_support_points[face_no][k] =
+ face_points.point(k);
Quadrature<dim> faces =
QProjector<dim>::project_to_all_faces(this->reference_cell_type(),
face_points);
- for (unsigned int k = 0;
- k < this->n_dofs_per_face() * GeometryInfo<dim>::faces_per_cell;
+ for (unsigned int k = 0; k < this->n_dofs_per_face(face_no) *
+ GeometryInfo<dim>::faces_per_cell;
++k)
- this->generalized_support_points[k] = faces.point(
- k +
- QProjector<dim>::DataSetDescriptor::face(this->reference_cell_type(),
- 0,
- true,
- false,
- false,
- this->n_dofs_per_face()));
-
- current = this->n_dofs_per_face() * GeometryInfo<dim>::faces_per_cell;
+ this->generalized_support_points[k] =
+ faces.point(k + QProjector<dim>::DataSetDescriptor::face(
+ this->reference_cell_type(),
+ 0,
+ true,
+ false,
+ false,
+ this->n_dofs_per_face(face_no)));
+
+ current =
+ this->n_dofs_per_face(face_no) * GeometryInfo<dim>::faces_per_cell;
}
if (deg == 1)
unsigned int fbase = 0;
unsigned int f = 0;
for (; f < GeometryInfo<dim>::faces_per_cell;
- ++f, fbase += this->n_dofs_per_face())
+ ++f, fbase += this->n_dofs_per_face(f))
{
- for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i)
+ for (unsigned int i = 0; i < this->n_dofs_per_face(f); ++i)
{
nodal_values[fbase + i] = support_point_values[fbase + i](
GeometryInfo<dim>::unit_normal_direction[f]);
// the DoFs within their own numbering. thus, translate to
// the base element numbering and then back
const std::pair<std::pair<unsigned int, unsigned int>, unsigned int>
- face_base_index = this->face_system_to_base_index(face_dof_index);
+ face_base_index = this->face_system_to_base_index(face_dof_index, face);
const unsigned int base_face_to_cell_index =
this->base_element(face_base_index.first.first)
void
FESystem<dim, spacedim>::build_interface_constraints()
{
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(this->n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
// check whether all base elements implement their interface constraint
// matrices. if this is not the case, then leave the interface costraints of
// this composed element empty as well; however, the rest of the element is
// data type, first value in pair is (base element,instance of base
// element), second is index within this instance
const std::pair<std::pair<unsigned int, unsigned int>, unsigned int>
- n_index = this->face_system_to_base_table[0][n];
+ n_index = this->face_system_to_base_table[face_no][n];
// likewise for the m index. this is more complicated due to the
// strange ordering we have for the dofs on the refined faces.
const unsigned int tmp1 =
2 * this->n_dofs_per_vertex() + index_in_line;
m_index.first =
- this->face_system_to_base_table[0][tmp1].first;
+ this->face_system_to_base_table[face_no][tmp1].first;
// what we are still missing is the index of m within the
// base elements interface_constraints table
// dof, we can construct the rest: tmp2 will denote the
// index of this shape function among the line shape
// functions:
- Assert(this->face_system_to_base_table[0][tmp1].second >=
- 2 * base_element(m_index.first.first)
- .n_dofs_per_vertex(),
- ExcInternalError());
+ Assert(
+ this->face_system_to_base_table[face_no][tmp1].second >=
+ 2 *
+ base_element(m_index.first.first).n_dofs_per_vertex(),
+ ExcInternalError());
const unsigned int tmp2 =
- this->face_system_to_base_table[0][tmp1].second -
+ this->face_system_to_base_table[face_no][tmp1].second -
2 * base_element(m_index.first.first).n_dofs_per_vertex();
Assert(tmp2 < base_element(m_index.first.first)
.n_dofs_per_line(),
const unsigned int tmp1 =
4 * this->n_dofs_per_vertex() + index_in_line;
m_index.first =
- this->face_system_to_base_table[0][tmp1].first;
+ this->face_system_to_base_table[face_no][tmp1].first;
- Assert(this->face_system_to_base_table[0][tmp1].second >=
- 4 * base_element(m_index.first.first)
- .n_dofs_per_vertex(),
- ExcInternalError());
+ Assert(
+ this->face_system_to_base_table[face_no][tmp1].second >=
+ 4 *
+ base_element(m_index.first.first).n_dofs_per_vertex(),
+ ExcInternalError());
const unsigned int tmp2 =
- this->face_system_to_base_table[0][tmp1].second -
+ this->face_system_to_base_table[face_no][tmp1].second -
4 * base_element(m_index.first.first).n_dofs_per_vertex();
Assert(tmp2 < base_element(m_index.first.first)
.n_dofs_per_line(),
const unsigned int index_in_quad =
(m - 5 * this->n_dofs_per_vertex() -
12 * this->n_dofs_per_line()) %
- this->n_dofs_per_quad();
- Assert(index_in_quad < this->n_dofs_per_quad(),
+ this->n_dofs_per_quad(face_no);
+ Assert(index_in_quad < this->n_dofs_per_quad(face_no),
ExcInternalError());
const unsigned int sub_quad =
((m - 5 * this->n_dofs_per_vertex() -
12 * this->n_dofs_per_line()) /
- this->n_dofs_per_quad());
+ this->n_dofs_per_quad(face_no));
Assert(sub_quad < 4, ExcInternalError());
const unsigned int tmp1 = 4 * this->n_dofs_per_vertex() +
4 * this->n_dofs_per_line() +
index_in_quad;
- Assert(tmp1 < this->face_system_to_base_table[0].size(),
+ Assert(tmp1 <
+ this->face_system_to_base_table[face_no].size(),
ExcInternalError());
m_index.first =
- this->face_system_to_base_table[0][tmp1].first;
-
- Assert(this->face_system_to_base_table[0][tmp1].second >=
- 4 * base_element(m_index.first.first)
- .n_dofs_per_vertex() +
- 4 * base_element(m_index.first.first)
- .n_dofs_per_line(),
- ExcInternalError());
+ this->face_system_to_base_table[face_no][tmp1].first;
+
+ Assert(
+ this->face_system_to_base_table[face_no][tmp1].second >=
+ 4 * base_element(m_index.first.first)
+ .n_dofs_per_vertex() +
+ 4 *
+ base_element(m_index.first.first).n_dofs_per_line(),
+ ExcInternalError());
const unsigned int tmp2 =
- this->face_system_to_base_table[0][tmp1].second -
+ this->face_system_to_base_table[face_no][tmp1].second -
4 *
base_element(m_index.first.first).n_dofs_per_vertex() -
4 * base_element(m_index.first.first).n_dofs_per_line();
Assert(tmp2 < base_element(m_index.first.first)
- .n_dofs_per_quad(),
+ .n_dofs_per_quad(face_no),
ExcInternalError());
m_index.second =
5 *
base_element(m_index.first.first).n_dofs_per_vertex() +
12 * base_element(m_index.first.first).n_dofs_per_line() +
- base_element(m_index.first.first).n_dofs_per_quad() *
+ base_element(m_index.first.first)
+ .n_dofs_per_quad(face_no) *
sub_quad +
tmp2;
}
// FiniteElement
this->system_to_component_table.resize(this->n_dofs_per_cell());
- this->face_system_to_component_table.resize(1);
- this->face_system_to_component_table[0].resize(this->n_dofs_per_face());
-
FETools::Compositing::build_cell_tables(this->system_to_base_table,
this->system_to_component_table,
this->component_to_base_table,
*this);
- FETools::Compositing::build_face_tables(
- this->face_system_to_base_table[0],
- this->face_system_to_component_table[0],
- *this);
+ this->face_system_to_component_table.resize(this->n_unique_faces());
+
+ for (unsigned int face_no = 0; face_no < this->n_unique_faces(); ++face_no)
+ {
+ this->face_system_to_component_table[0].resize(
+ this->n_dofs_per_face(face_no));
+
+ FETools::Compositing::build_face_tables(
+ this->face_system_to_base_table[face_no],
+ this->face_system_to_component_table[face_no],
+ *this,
+ true,
+ face_no);
+ }
}
// now initialize interface constraints, support points, and other tables.
// initialize face support points (for dim==2,3). same procedure as above
if (dim > 1)
init_tasks += Threads::new_task([&]() {
- // if one of the base elements has no support points, then it makes no
- // sense to define support points for the composed element. In that case,
- // return an empty array to demonstrate that fact (note that we ask
- // whether the base element has no support points at all, not only none on
- // the face!)
- //
- // on the other hand, if there is an element that simply has no degrees of
- // freedom on the face at all, then we don't care whether it has support
- // points or not. this is, for example, the case for the stable Stokes
- // element Q(p)^dim \times DGP(p-1).
- for (unsigned int base_el = 0; base_el < this->n_base_elements();
- ++base_el)
- if (!base_element(base_el).has_support_points() &&
- (base_element(base_el).n_dofs_per_face() > 0))
- {
- this->unit_face_support_points[0].resize(0);
- return;
- }
-
+ for (unsigned int face_no = 0; face_no < this->n_unique_faces();
+ ++face_no)
+ {
+ // if one of the base elements has no support points, then it makes
+ // no sense to define support points for the composed element. In
+ // that case, return an empty array to demonstrate that fact (note
+ // that we ask whether the base element has no support points at
+ // all, not only none on the face!)
+ //
+ // on the other hand, if there is an element that simply has no
+ // degrees of freedom on the face at all, then we don't care whether
+ // it has support points or not. this is, for example, the case for
+ // the stable Stokes element Q(p)^dim \times DGP(p-1).
+ bool flag_has_no_support_points = false;
+
+ for (unsigned int base_el = 0; base_el < this->n_base_elements();
+ ++base_el)
+ if (!base_element(base_el).has_support_points() &&
+ (base_element(base_el).n_dofs_per_face(face_no) > 0))
+ {
+ this->unit_face_support_points[face_no].resize(0);
+ flag_has_no_support_points = true;
+ break;
+ }
- // generate unit face support points from unit support points of sub
- // elements
- this->unit_face_support_points[0].resize(this->n_dofs_per_face());
- for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i)
- {
- const unsigned int base_i =
- this->face_system_to_base_table[0][i].first.first;
- const unsigned int index_in_base =
- this->face_system_to_base_table[0][i].second;
+ if (flag_has_no_support_points)
+ continue;
- Assert(index_in_base <
- base_element(base_i).unit_face_support_points[0].size(),
- ExcInternalError());
+ // generate unit face support points from unit support points of sub
+ // elements
+ this->unit_face_support_points[face_no].resize(
+ this->n_dofs_per_face(face_no));
- this->unit_face_support_points[0][i] =
- base_element(base_i).unit_face_support_points[0][index_in_base];
+ for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
+ {
+ const unsigned int base_i =
+ this->face_system_to_base_table[face_no][i].first.first;
+ const unsigned int index_in_base =
+ this->face_system_to_base_table[face_no][i].second;
+
+ Assert(
+ index_in_base <
+ base_element(base_i).unit_face_support_points[face_no].size(),
+ ExcInternalError());
+
+ this->unit_face_support_points[face_no][i] =
+ base_element(base_i)
+ .unit_face_support_points[face_no][index_in_base];
+ }
}
});
// initialize quad dof index permutation in 3d and higher
if (dim >= 3)
init_tasks += Threads::new_task([&]() {
- // the array into which we want to write should have the correct size
- // already.
- Assert(this->adjust_quad_dof_index_for_face_orientation_table[0]
- .n_elements() == 8 * this->n_dofs_per_quad(),
- ExcInternalError());
-
- // to obtain the shifts for this composed element, copy the shift
- // information of the base elements
- unsigned int index = 0;
- for (unsigned int b = 0; b < this->n_base_elements(); ++b)
+ for (unsigned int face_no = 0; face_no < this->n_unique_faces();
+ ++face_no)
{
- const Table<2, int> &temp =
- this->base_element(b)
- .adjust_quad_dof_index_for_face_orientation_table[0];
- for (unsigned int c = 0; c < this->element_multiplicity(b); ++c)
+ // the array into which we want to write should have the correct size
+ // already.
+ Assert(this->adjust_quad_dof_index_for_face_orientation_table[face_no]
+ .n_elements() == 8 * this->n_dofs_per_quad(face_no),
+ ExcInternalError());
+
+ // to obtain the shifts for this composed element, copy the shift
+ // information of the base elements
+ unsigned int index = 0;
+ for (unsigned int b = 0; b < this->n_base_elements(); ++b)
{
- for (unsigned int i = 0; i < temp.size(0); ++i)
- for (unsigned int j = 0; j < 8; ++j)
- this->adjust_quad_dof_index_for_face_orientation_table[0](
- index + i, j) = temp(i, j);
- index += temp.size(0);
+ const Table<2, int> &temp =
+ this->base_element(b)
+ .adjust_quad_dof_index_for_face_orientation_table[face_no];
+ for (unsigned int c = 0; c < this->element_multiplicity(b); ++c)
+ {
+ for (unsigned int i = 0; i < temp.size(0); ++i)
+ for (unsigned int j = 0; j < 8; ++j)
+ this->adjust_quad_dof_index_for_face_orientation_table
+ [face_no](index + i, j) = temp(i, j);
+ index += temp.size(0);
+ }
}
+ Assert(index == this->n_dofs_per_quad(face_no), ExcInternalError());
}
- Assert(index == this->n_dofs_per_quad(), ExcInternalError());
// additionally compose the permutation information for lines
Assert(this->adjust_line_dof_index_for_line_orientation_table.size() ==
this->n_dofs_per_line(),
ExcInternalError());
- index = 0;
+ unsigned int index = 0;
for (unsigned int b = 0; b < this->n_base_elements(); ++b)
{
const std::vector<int> &temp2 =
FullMatrix<double> & interpolation_matrix,
const unsigned int face_no) const
{
- Assert(interpolation_matrix.n() == this->n_dofs_per_face(),
+ Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
ExcDimensionMismatch(interpolation_matrix.n(),
- this->n_dofs_per_face()));
- Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(),
+ this->n_dofs_per_face(face_no)));
+ Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
ExcDimensionMismatch(interpolation_matrix.m(),
- x_source_fe.n_dofs_per_face()));
+ x_source_fe.n_dofs_per_face(face_no)));
// since dofs for each base are independent, we only have to stack things up
// from base element to base element
ExcNotImplemented());
// get the interpolation from the bases
- base_to_base_interpolation.reinit(base_other.n_dofs_per_face(),
- base.n_dofs_per_face());
+ base_to_base_interpolation.reinit(base_other.n_dofs_per_face(face_no),
+ base.n_dofs_per_face(face_no));
base.get_face_interpolation_matrix(base_other,
base_to_base_interpolation,
face_no);
// now translate entries. we'd like to have something like
// face_base_to_system_index, but that doesn't exist. rather, all we
// have is the reverse. well, use that then
- for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i)
- if (this->face_system_to_base_index(i).first ==
+ for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
+ if (this->face_system_to_base_index(i, face_no).first ==
std::make_pair(base_index, multiplicity))
- for (unsigned int j = 0; j < fe_other_system->n_dofs_per_face();
+ for (unsigned int j = 0;
+ j < fe_other_system->n_dofs_per_face(face_no);
++j)
- if (fe_other_system->face_system_to_base_index(j).first ==
+ if (fe_other_system->face_system_to_base_index(j, face_no)
+ .first ==
std::make_pair(base_index_other, multiplicity_other))
interpolation_matrix(j, i) = base_to_base_interpolation(
- fe_other_system->face_system_to_base_index(j).second,
- this->face_system_to_base_index(i).second);
+ fe_other_system->face_system_to_base_index(j, face_no)
+ .second,
+ this->face_system_to_base_index(i, face_no).second);
// advance to the next base element for this and the other fe_system;
// see if we can simply advance the multiplicity by one, or if have to
(dynamic_cast<const FESystem<dim, spacedim> *>(&x_source_fe) != nullptr),
(typename FiniteElement<dim, spacedim>::ExcInterpolationNotImplemented()));
- Assert(interpolation_matrix.n() == this->n_dofs_per_face(),
+ Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
ExcDimensionMismatch(interpolation_matrix.n(),
- this->n_dofs_per_face()));
- Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(),
+ this->n_dofs_per_face(face_no)));
+ Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
ExcDimensionMismatch(interpolation_matrix.m(),
- x_source_fe.n_dofs_per_face()));
+ x_source_fe.n_dofs_per_face(face_no)));
// since dofs for each base are independent, we only have to stack things up
// from base element to base element
ExcNotImplemented());
// get the interpolation from the bases
- base_to_base_interpolation.reinit(base_other.n_dofs_per_face(),
- base.n_dofs_per_face());
+ base_to_base_interpolation.reinit(base_other.n_dofs_per_face(face_no),
+ base.n_dofs_per_face(face_no));
base.get_subface_interpolation_matrix(base_other,
subface,
base_to_base_interpolation,
// now translate entries. we'd like to have something like
// face_base_to_system_index, but that doesn't exist. rather, all we
// have is the reverse. well, use that then
- for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i)
- if (this->face_system_to_base_index(i).first ==
+ for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
+ if (this->face_system_to_base_index(i, face_no).first ==
std::make_pair(base_index, multiplicity))
- for (unsigned int j = 0; j < fe_other_system->n_dofs_per_face();
+ for (unsigned int j = 0;
+ j < fe_other_system->n_dofs_per_face(face_no);
++j)
- if (fe_other_system->face_system_to_base_index(j).first ==
+ if (fe_other_system->face_system_to_base_index(j, face_no)
+ .first ==
std::make_pair(base_index_other, multiplicity_other))
interpolation_matrix(j, i) = base_to_base_interpolation(
- fe_other_system->face_system_to_base_index(j).second,
- this->face_system_to_base_index(i).second);
+ fe_other_system->face_system_to_base_index(j, face_no)
+ .second,
+ this->face_system_to_base_index(i, face_no).second);
// advance to the next base element for this and the other fe_system;
// see if we can simply advance the multiplicity by one, or if have to
template <int structdim>
std::vector<std::pair<unsigned int, unsigned int>>
FESystem<dim, spacedim>::hp_object_dof_identities(
- const FiniteElement<dim, spacedim> &fe_other) const
+ const FiniteElement<dim, spacedim> &fe_other,
+ const unsigned int face_no) const
{
// since dofs on each subobject (vertex, line, ...) are ordered such that
// first come all from the first base element all multiplicities, then
base_identities = base.hp_line_dof_identities(base_other);
break;
case 2:
- base_identities = base.hp_quad_dof_identities(base_other);
+ base_identities =
+ base.hp_quad_dof_identities(base_other, face_no);
break;
default:
Assert(false, ExcNotImplemented());
std::vector<std::pair<unsigned int, unsigned int>>
FESystem<dim, spacedim>::hp_quad_dof_identities(
const FiniteElement<dim, spacedim> &fe_other,
- const unsigned int) const
+ const unsigned int face_no) const
{
- return hp_object_dof_identities<2>(fe_other);
+ return hp_object_dof_identities<2>(fe_other, face_no);
}
const unsigned int index,
const unsigned int face_no) const
{
- AssertIndexRange(index, this->n_dofs_per_face());
- Assert((this->unit_face_support_points[0].size() ==
- this->n_dofs_per_face()) ||
- (this->unit_face_support_points[0].size() == 0),
- (typename FiniteElement<dim, spacedim>::ExcFEHasNoSupportPoints()));
+ AssertIndexRange(index, this->n_dofs_per_face(face_no));
+ Assert(
+ (this->unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
+ .size() == this->n_dofs_per_face(face_no)) ||
+ (this->unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
+ .size() == 0),
+ (typename FiniteElement<dim, spacedim>::ExcFEHasNoSupportPoints()));
// let's see whether we have the information pre-computed
- if (this->unit_face_support_points[0].size() != 0)
- return this->unit_face_support_points[0][index];
+ if (this->unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
+ .size() != 0)
+ return this
+ ->unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
+ [index];
else
// no. ask the base element whether it would like to provide this
// information
return (
- base_element(this->face_system_to_base_index(index).first.first)
- .unit_face_support_point(this->face_system_to_base_index(index).second,
- face_no));
+ base_element(this->face_system_to_base_index(index, face_no).first.first)
+ .unit_face_support_point(
+ this->face_system_to_base_index(index, face_no).second, face_no));
}
std::vector<std::pair<unsigned int, unsigned int>>
&face_system_to_component_table,
const FiniteElement<deal_II_dimension, deal_II_space_dimension> &fe,
- const bool do_tensor_product);
+ const bool do_tensor_product,
+ const unsigned int face_no);
\}
FETools::hierarchic_to_lexicographic_numbering<dim - 1>(degree));
// Initialize face support points
- this->unit_face_support_points[0] = fe_q.get_unit_face_support_points();
+ AssertDimension(this->n_unique_faces(), fe_q.n_unique_faces());
+ for (unsigned int face_no = 0; face_no < this->n_unique_faces(); ++face_no)
+ this->unit_face_support_points[face_no] =
+ fe_q.get_unit_face_support_points(face_no);
// initialize unit support points (this makes it possible to assign initial
// values to FE_TraceQ). Note that we simply take the points of fe_q but
const FiniteElement<dim, spacedim> &x_source_fe,
const unsigned int subface,
FullMatrix<double> & interpolation_matrix,
- const unsigned int) const
+ const unsigned int face_no) const
{
// this is the code from FE_FaceQ
- Assert(interpolation_matrix.n() == this->n_dofs_per_face(),
+ Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
ExcDimensionMismatch(interpolation_matrix.n(),
- this->n_dofs_per_face()));
- Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(),
+ this->n_dofs_per_face(face_no)));
+ Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
ExcDimensionMismatch(interpolation_matrix.m(),
- x_source_fe.n_dofs_per_face()));
+ x_source_fe.n_dofs_per_face(face_no)));
// see if source is a FaceQ element
if (const FE_TraceQ<dim, spacedim> *source_fe =
{
fe_q.get_subface_interpolation_matrix(source_fe->fe_q,
subface,
- interpolation_matrix);
+ interpolation_matrix,
+ face_no);
}
else if (dynamic_cast<const FE_Nothing<dim> *>(&x_source_fe) != nullptr)
{
projects / dealii.git / commitdiff