dof_handler.object_dof_ptr[obj_level][d][obj_index];
const unsigned int ptr_1 =
ptr_0 +
- dof_handler.get_fe(fe_index).template n_dofs_per_object<dim>();
+ dof_handler.get_fe(fe_index).template n_dofs_per_object<dim>(0);
return {ptr_0, ptr_1};
}
{
const auto &fe = accessor.get_fe(fe_index_);
- const unsigned int //
- dofs_per_vertex = fe.n_dofs_per_vertex(), //
- dofs_per_line = fe.n_dofs_per_line(), //
- dofs_per_quad = fe.n_dofs_per_quad(), //
- dofs_per_hex = fe.n_dofs_per_hex(); //
+ const unsigned int //
+ dofs_per_vertex = fe.n_dofs_per_vertex(), //
+ dofs_per_line = fe.n_dofs_per_line(), //
+ dofs_per_quad = fe.n_dofs_per_quad(0 /*dummy*/), //
+ dofs_per_hex = fe.n_dofs_per_hex(); //
const unsigned int inner_dofs =
structdim == 1 ? dofs_per_line :
(void)fe_index;
for (unsigned int d = 0;
- d < fe.template n_dofs_per_object<structdim_>();
+ d < fe.template n_dofs_per_object<structdim_>(0);
++d, ++index)
index_value[index] = accessor.mg_dof_index(level, mapping(d));
}
(void)fe_index;
for (unsigned int d = 0;
- d < fe.template n_dofs_per_object<structdim_>();
+ d < fe.template n_dofs_per_object<structdim_>(0);
++d, ++index)
accessor.set_mg_dof_index(level, mapping(d), index_value[index]);
}
void
setup_constraint_weigths(unsigned int fe_degree)
{
+ const unsigned int face_no =
+ 0; // we assume that all faces have the same number of dofs
+
FE_Q<2> fe_q(fe_degree);
- FullMatrix<double> interpolation_matrix(fe_q.n_dofs_per_face(),
- fe_q.n_dofs_per_face());
- fe_q.get_subface_interpolation_matrix(fe_q, 0, interpolation_matrix);
+ FullMatrix<double> interpolation_matrix(fe_q.n_dofs_per_face(face_no),
+ fe_q.n_dofs_per_face(face_no));
+ fe_q.get_subface_interpolation_matrix(fe_q,
+ 0,
+ interpolation_matrix,
+ face_no);
std::vector<unsigned int> mapping =
FETools::lexicographic_to_hierarchic_numbering<1>(fe_degree);
- FullMatrix<double> mapped_matrix(fe_q.n_dofs_per_face(),
- fe_q.n_dofs_per_face());
- for (unsigned int i = 0; i < fe_q.n_dofs_per_face(); ++i)
- for (unsigned int j = 0; j < fe_q.n_dofs_per_face(); ++j)
+ FullMatrix<double> mapped_matrix(fe_q.n_dofs_per_face(face_no),
+ fe_q.n_dofs_per_face(face_no));
+ for (unsigned int i = 0; i < fe_q.n_dofs_per_face(face_no); ++i)
+ for (unsigned int j = 0; j < fe_q.n_dofs_per_face(face_no); ++j)
mapped_matrix(i, j) = interpolation_matrix(mapping[i], mapping[j]);
cudaError_t error_code =
cudaMemcpyToSymbol(internal::constraint_weights,
&mapped_matrix[0][0],
- sizeof(double) * fe_q.n_dofs_per_face() *
- fe_q.n_dofs_per_face());
+ sizeof(double) * fe_q.n_dofs_per_face(face_no) *
+ fe_q.n_dofs_per_face(face_no));
AssertCuda(error_code);
}
} // namespace internal
const FiniteElement<dim> &fe = fes[fe_index];
// cache number of finite elements and dofs_per_cell
dof_info[no].dofs_per_cell.push_back(fe.n_dofs_per_cell());
- dof_info[no].dofs_per_face.push_back(fe.n_dofs_per_face());
+ dof_info[no].dofs_per_face.push_back(fe.n_dofs_per_face(
+ 0)); // we assume that all faces have the same number of dofs
dof_info[no].dimension = dim;
dof_info[no].n_base_elements = fe.n_base_elements();
dof_info[no].n_components.resize(dof_info[no].n_base_elements);
}
const unsigned int dofs_per_face =
- cell->face(f)->get_fe(0).n_dofs_per_face();
+ dof.get_fe(0).n_dofs_per_face(f);
std::vector<types::global_dof_index> dofs_1(dofs_per_face);
std::vector<types::global_dof_index> dofs_2(dofs_per_face);
const bool fe_is_system = (n_components != 1);
const bool fe_is_primitive = fe.is_primitive();
- const unsigned int dofs_per_face = fe.n_dofs_per_face();
-
copy_data.cell = cell;
copy_data.dofs_per_cell = fe.n_dofs_per_cell();
copy_data.dofs.resize(copy_data.dofs_per_cell);
cell->get_dof_indices(copy_data.dofs);
- std::vector<types::global_dof_index> dofs_on_face_vector(dofs_per_face);
+ std::vector<types::global_dof_index> dofs_on_face_vector;
// Because CopyData objects are reused and emplace_back is
// used, dof_is_on_face, cell_matrix, and cell_vector must be
}
+ dofs_on_face_vector.resize(fe.n_dofs_per_face(face));
cell->face(face)->get_dof_indices(dofs_on_face_vector);
// for each dof on the cell, have a flag whether it is on
// the face
const FiniteElement<dim, spacedim> &fe = cell->get_fe();
const bool fe_is_system = (n_components != 1);
const bool fe_is_primitive = fe.is_primitive();
- const unsigned int dofs_per_face = fe.n_dofs_per_face();
copy_data.cell = cell;
copy_data.dofs_per_cell = fe.n_dofs_per_cell();
std::vector<number> rhs_values_scalar;
std::vector<Vector<number>> rhs_values_system;
- std::vector<types::global_dof_index> dofs_on_face_vector(dofs_per_face);
+ std::vector<types::global_dof_index> dofs_on_face_vector;
copy_data.dofs.resize(copy_data.dofs_per_cell);
cell->get_dof_indices(copy_data.dofs);
}
}
+ dofs_on_face_vector.resize(fe.n_dofs_per_face(face));
cell->face(face)->get_dof_indices(dofs_on_face_vector,
cell->active_fe_index());
// for each dof on the cell, have a
function_values);
for (unsigned int i = 0; i < fe.n_dofs_per_vertex(); ++i)
- if (component_mask[fe.face_system_to_component_index(i)
+ if (component_mask[fe.face_system_to_component_index(
+ i, direction)
.first])
boundary_values[cell->vertex_dof_index(
direction, i, cell->active_fe_index())] =
function_values(
- fe.face_system_to_component_index(i).first);
+ fe.face_system_to_component_index(i, direction)
+ .first);
}
}
else // dim > 1
dof_values_system.reserve(
dof.get_fe_collection().max_dofs_per_face());
+ // TODO: get support for each face -> PR #10764
+ const unsigned int face_no = 0;
+
// before we start with the loop over all cells create an hp::FEValues
// object that holds the interpolation points of all finite elements
// that may ever be in use
//
// to do this, we check whether the FE has support points on the
// face at all:
- if (fe.has_face_support_points())
- q_collection.push_back(
- Quadrature<dim - 1>(fe.get_unit_face_support_points()));
+ if (fe.has_face_support_points(face_no))
+ q_collection.push_back(Quadrature<dim - 1>(
+ fe.get_unit_face_support_points(face_no)));
else
{
// if not, then we should try a more clever way. the idea is
// components we are interested in are primitive (by the above
// check), we can safely put such a check in front
std::vector<Point<dim - 1>> unit_support_points(
- fe.n_dofs_per_face());
+ fe.n_dofs_per_face(face_no));
- for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i)
- if (fe.is_primitive(fe.face_to_cell_index(i, 0)))
- if (component_mask[fe.face_system_to_component_index(i)
+ for (unsigned int i = 0; i < fe.n_dofs_per_face(face_no); ++i)
+ if (fe.is_primitive(fe.face_to_cell_index(i, face_no)))
+ if (component_mask[fe.face_system_to_component_index(
+ i, face_no)
.first] == true)
- unit_support_points[i] = fe.unit_face_support_point(i);
+ unit_support_points[i] =
+ fe.unit_face_support_point(i, face_no);
q_collection.push_back(
Quadrature<dim - 1>(unit_support_points));
// element in use here has DoFs on the face at all
if ((function_map.find(boundary_component) !=
function_map.end()) &&
- (cell->get_fe().n_dofs_per_face() > 0))
+ (cell->get_fe().n_dofs_per_face(face_no) > 0))
{
// face is of the right component
x_fe_values.reinit(cell, face_no);
// get indices, physical location and boundary values of
// dofs on this face
- face_dofs.resize(fe.n_dofs_per_face());
+ face_dofs.resize(fe.n_dofs_per_face(face_no));
face->get_dof_indices(face_dofs, cell->active_fe_index());
const std::vector<Point<spacedim>> &dof_locations =
fe_values.get_quadrature_points();
{
// resize array. avoid construction of a memory
// allocating temporary if possible
- if (dof_values_system.size() < fe.n_dofs_per_face())
- dof_values_system.resize(fe.n_dofs_per_face(),
- Vector<number>(
- fe.n_components()));
+ if (dof_values_system.size() <
+ fe.n_dofs_per_face(face_no))
+ dof_values_system.resize(
+ fe.n_dofs_per_face(face_no),
+ Vector<number>(fe.n_components()));
else
- dof_values_system.resize(fe.n_dofs_per_face());
+ dof_values_system.resize(
+ fe.n_dofs_per_face(face_no));
function_map.find(boundary_component)
->second->vector_value_list(dof_locations,
unsigned int component;
if (fe.is_primitive())
component =
- fe.face_system_to_component_index(i).first;
+ fe.face_system_to_component_index(i, face_no)
+ .first;
else
{
// non-primitive case. make sure that this
// fe has only one component, so save some computations
{
// get only the one component that this function has
- dof_values_scalar.resize(fe.n_dofs_per_face());
+ dof_values_scalar.resize(fe.n_dofs_per_face(face_no));
function_map.find(boundary_component)
->second->value_list(dof_locations,
dof_values_scalar,
// Compute the degrees of
// freedom.
- for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i)
+ for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
if (((dynamic_cast<const FESystem<dim> *>(&fe) != nullptr) &&
(fe.system_to_base_index(fe.face_to_cell_index(i, face))
.first == base_indices) &&
// Compute the degrees
// of freedom.
- for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i)
+ for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
if (((dynamic_cast<const FESystem<dim> *>(&fe) !=
nullptr) &&
(fe.system_to_base_index(
for (unsigned int d = 0; d < dim; ++d)
tmp[d] = values[q_point](first_vector_component + d);
- for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i)
+ for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
if (((dynamic_cast<const FESystem<dim> *>(&fe) !=
nullptr) &&
(fe.system_to_base_index(
unsigned int index = 0;
- for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i)
+ for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
if (((dynamic_cast<const FESystem<dim> *>(&fe) !=
nullptr) &&
(fe.system_to_base_index(
{
unsigned int index = 0;
- for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i)
+ for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
if (((dynamic_cast<const FESystem<dim> *>(&fe) != nullptr) &&
(fe.system_to_base_index(fe.face_to_cell_index(i, face))
.first == base_indices) &&
for (unsigned int d = 0; d < dim; ++d)
tmp[d] = values[q_point](first_vector_component + d);
- for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i)
+ for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
if (((dynamic_cast<const FESystem<dim> *>(&fe) !=
nullptr) &&
(fe.system_to_base_index(
unsigned int index = 0;
- for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i)
+ for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
if (((dynamic_cast<const FESystem<dim> *>(&fe) !=
nullptr) &&
(fe.system_to_base_index(
unsigned int index = 0;
- for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i)
+ for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
if (((dynamic_cast<const FESystem<dim> *>(&fe) != nullptr) &&
(fe.system_to_base_index(fe.face_to_cell_index(i, face))
.first == base_indices) &&
// corresponding degrees of freedom.
const unsigned int superdegree = dof_handler.get_fe().degree;
const QGauss<dim - 1> reference_face_quadrature(2 * superdegree);
- const unsigned int dofs_per_face = dof_handler.get_fe().n_dofs_per_face();
+
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(dof_handler.get_fe().n_unique_faces(), 1);
+ const unsigned int dofs_per_face = dof_handler.get_fe().n_dofs_per_face(0);
+
const hp::FECollection<dim> &fe_collection(dof_handler.get_fe_collection());
const hp::MappingCollection<dim> mapping_collection(mapping);
hp::QCollection<dim> face_quadrature_collection;
}
const unsigned int dofs_per_face =
- cell->get_fe().n_dofs_per_face();
+ cell->get_fe().n_dofs_per_face(face);
dofs_processed.resize(dofs_per_face);
dof_values.resize(dofs_per_face);
const unsigned int superdegree = cell->get_fe().degree;
const unsigned int degree = superdegree - 1;
const unsigned int dofs_per_face =
- cell->get_fe().n_dofs_per_face();
+ cell->get_fe().n_dofs_per_face(face);
dofs_processed.resize(dofs_per_face);
dof_values.resize(dofs_per_face);
1);
unsigned int associated_edge_dof_index = 0;
- for (unsigned int face_idx = 0; face_idx < fe.n_dofs_per_face();
+ for (unsigned int face_idx = 0;
+ face_idx < fe.n_dofs_per_face(face);
++face_idx)
{
const unsigned int cell_idx =
// Loop over these quad-interior dofs.
unsigned int associated_face_dof_index = 0;
for (unsigned int quad_dof_idx = 0;
- quad_dof_idx < fe.n_dofs_per_quad();
+ quad_dof_idx < fe.n_dofs_per_quad(face);
++quad_dof_idx)
{
const unsigned int face_idx =
}
const unsigned int dofs_per_face =
- cell->get_fe().n_dofs_per_face();
+ cell->get_fe().n_dofs_per_face(face);
dofs_processed.resize(dofs_per_face);
dof_values.resize(dofs_per_face);
cell->get_fe().degree;
const unsigned int degree = superdegree - 1;
const unsigned int dofs_per_face =
- cell->get_fe().n_dofs_per_face();
+ cell->get_fe().n_dofs_per_face(face);
dofs_processed.resize(dofs_per_face);
dof_values.resize(dofs_per_face);
0};
std::vector<Vector<double>> values(fe_values.n_quadrature_points,
Vector<double>(2));
- Vector<double> dof_values(fe.n_dofs_per_face());
+ Vector<double> dof_values(fe.n_dofs_per_face(face));
// Get the values of the boundary function at the quadrature points.
{
jacobians[q_point][1][face_coordinate_direction[face]] *
jacobians[q_point][1][face_coordinate_direction[face]]);
- for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i)
+ for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
dof_values(i) +=
tmp * (normals[q_point] *
fe_values[vec].value(
}
std::vector<types::global_dof_index> face_dof_indices(
- fe.n_dofs_per_face());
+ fe.n_dofs_per_face(face));
cell->face(face)->get_dof_indices(face_dof_indices,
cell->active_fe_index());
// Copy the computed values in the AffineConstraints only, if the degree
// of freedom is not already constrained.
- for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i)
+ for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
if (!(constraints.is_constrained(face_dof_indices[i])) &&
fe.get_nonzero_components(fe.face_to_cell_index(
i,
{1, 2}, {1, 2}, {2, 0}, {2, 0}, {0, 1}, {0, 1}};
std::vector<Vector<double>> values(fe_values.n_quadrature_points,
Vector<double>(3));
- Vector<double> dof_values_local(fe.n_dofs_per_face());
+ Vector<double> dof_values_local(fe.n_dofs_per_face(face));
{
const std::vector<Point<3>> &quadrature_points =
jacobians[q_point][2][face_coordinate_directions[face][1]] *
jacobians[q_point][2][face_coordinate_directions[face][1]]));
- for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i)
+ for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
dof_values_local(i) +=
tmp * (normals[q_point] *
fe_values[vec].value(
}
std::vector<types::global_dof_index> face_dof_indices(
- fe.n_dofs_per_face());
+ fe.n_dofs_per_face(face));
cell->face(face)->get_dof_indices(face_dof_indices,
cell->active_fe_index());
- for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i)
+ for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
if (projected_dofs[face_dof_indices[i]] < fe.degree &&
fe.get_nonzero_components(fe.face_to_cell_index(
i,
for (unsigned int i = 0; i < fe_collection.size(); ++i)
mapping_collection.push_back(mapping);
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(dof_handler.get_fe().n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
// now also create a quadrature collection for the faces of a cell. fill
// it with a quadrature formula with the support points on faces for each
// FE
for (unsigned int i = 0; i < fe_collection.size(); ++i)
{
const std::vector<Point<dim - 1>> &unit_support_points =
- fe_collection[i].get_unit_face_support_points();
+ fe_collection[i].get_unit_face_support_points(face_no);
- Assert(unit_support_points.size() == fe_collection[i].n_dofs_per_face(),
+ Assert(unit_support_points.size() ==
+ fe_collection[i].n_dofs_per_face(face_no),
ExcInternalError());
face_quadrature_collection.push_back(
cell->face(face_no);
// get the indices of the dofs on this cell...
- face_dofs.resize(fe.n_dofs_per_face());
+ face_dofs.resize(fe.n_dofs_per_face(face_no));
face->get_dof_indices(face_dofs, cell->active_fe_index());
x_fe_face_values.reinit(cell, face_no);
// then identify which of them correspond to the selected set of
// vector components
for (unsigned int i = 0; i < face_dofs.size(); ++i)
- if (fe.face_system_to_component_index(i).first ==
+ if (fe.face_system_to_component_index(i, face_no).first ==
first_vector_component)
{
// find corresponding other components of vector
"Error: the finite element does not have enough components "
"to define a normal direction."));
- for (unsigned int k = 0; k < fe.n_dofs_per_face(); ++k)
+ for (unsigned int k = 0; k < fe.n_dofs_per_face(face_no);
+ ++k)
if ((k != i) &&
(face_quadrature_collection[cell->active_fe_index()]
.point(k) ==
face_quadrature_collection[cell->active_fe_index()]
.point(i)) &&
- (fe.face_system_to_component_index(k).first >=
- first_vector_component) &&
- (fe.face_system_to_component_index(k).first <
+ (fe.face_system_to_component_index(k, face_no)
+ .first >= first_vector_component) &&
+ (fe.face_system_to_component_index(k, face_no).first <
first_vector_component + dim))
vector_dofs.dof_indices
- [fe.face_system_to_component_index(k).first -
+ [fe.face_system_to_component_index(k, face_no).first -
first_vector_component] = face_dofs[k];
for (unsigned int d = 0; d < dim; ++d)
for (unsigned int i = 0; i < fe_collection.size(); ++i)
mapping_collection.push_back(mapping);
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(dof_handler.get_fe().n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
// now also create a quadrature collection for the faces of a cell. fill
// it with a quadrature formula with the support points on faces for each
// FE
for (unsigned int i = 0; i < fe_collection.size(); ++i)
{
const std::vector<Point<dim - 1>> &unit_support_points =
- fe_collection[i].get_unit_face_support_points();
+ fe_collection[i].get_unit_face_support_points(face_no);
- Assert(unit_support_points.size() == fe_collection[i].n_dofs_per_face(),
+ Assert(unit_support_points.size() ==
+ fe_collection[i].n_dofs_per_face(face_no),
ExcInternalError());
face_quadrature_collection.push_back(
cell->face(face_no);
// get the indices of the dofs on this cell...
- face_dofs.resize(fe.n_dofs_per_face());
+ face_dofs.resize(fe.n_dofs_per_face(face_no));
face->get_dof_indices(face_dofs, cell->active_fe_index());
x_fe_face_values.reinit(cell, face_no);
std::map<types::global_dof_index, double> dof_to_b_value;
unsigned int n_scalar_indices = 0;
- cell_vector_dofs.resize(fe.n_dofs_per_face());
- for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i)
+ cell_vector_dofs.resize(fe.n_dofs_per_face(face_no));
+ for (unsigned int i = 0; i < fe.n_dofs_per_face(face_no); ++i)
{
- if (fe.face_system_to_component_index(i).first >=
+ if (fe.face_system_to_component_index(i, face_no).first >=
first_vector_component &&
- fe.face_system_to_component_index(i).first <
+ fe.face_system_to_component_index(i, face_no).first <
first_vector_component + dim)
{
const unsigned int component =
- fe.face_system_to_component_index(i).first -
+ fe.face_system_to_component_index(i, face_no).first -
first_vector_component;
- n_scalar_indices =
- std::max(n_scalar_indices,
- fe.face_system_to_component_index(i).second +
- 1);
- cell_vector_dofs[fe.face_system_to_component_index(i)
+ n_scalar_indices = std::max(
+ n_scalar_indices,
+ fe.face_system_to_component_index(i, face_no).second +
+ 1);
+ cell_vector_dofs[fe.face_system_to_component_index(
+ i, face_no)
.second][component] = face_dofs[i];
const Point<dim> point = fe_values.quadrature_point(i);
for (auto f : GeometryInfo<dim>::face_indices())
if (cell->at_boundary(f))
{
- face_dof_indices.resize(cell->get_fe().n_dofs_per_face());
+ face_dof_indices.resize(cell->get_fe().n_dofs_per_face(f));
cell->face(f)->get_dof_indices(face_dof_indices,
cell->active_fe_index());
- for (unsigned int i = 0; i < cell->get_fe().n_dofs_per_face();
+ for (unsigned int i = 0; i < cell->get_fe().n_dofs_per_face(f);
++i)
// enter zero boundary values
// for all boundary nodes
{
dof_handler.object_dof_indices[i][2].resize(
dof_handler.tria->n_raw_cells(i) *
- dof_handler.get_fe().n_dofs_per_quad(),
+ dof_handler.get_fe().n_dofs_per_quad(
+ 0 /*note: in 2D there is only one quad*/),
numbers::invalid_dof_index);
dof_handler.object_dof_ptr[i][2].reserve(
for (unsigned int j = 0; j < dof_handler.tria->n_raw_cells(i) + 1;
j++)
dof_handler.object_dof_ptr[i][2].push_back(
- j * dof_handler.get_fe().n_dofs_per_quad());
+ j * dof_handler.get_fe().n_dofs_per_quad(
+ 0 /*note: in 2D there is only one quad*/));
dof_handler.cell_dof_cache_indices[i].resize(
dof_handler.tria->n_raw_cells(i) *
{
const auto &face = cell->face(face_index);
const auto n_dofs_per_quad =
- dof_handler.get_fe().n_dofs_per_quad(/*face_index*/);
+ dof_handler.get_fe().n_dofs_per_quad(face_index);
auto &n_dofs_per_quad_target =
dof_handler.object_dof_ptr[0][2][face->index() + 1];
std::make_unique<
internal::DoFHandlerImplementation::DoFLevel<2>>());
dof_handler.mg_levels.back()->dof_object.dofs =
- std::vector<types::global_dof_index>(tria.n_raw_quads(i) *
- fe.n_dofs_per_quad(),
- numbers::invalid_dof_index);
+ std::vector<types::global_dof_index>(
+ tria.n_raw_quads(i) *
+ fe.n_dofs_per_quad(0 /*note: in 2D there is only one quad*/),
+ numbers::invalid_dof_index);
}
dof_handler.mg_faces =
std::vector<types::global_dof_index>(tria.n_raw_lines() *
fe.n_dofs_per_line(),
numbers::invalid_dof_index);
- dof_handler.mg_faces->quads.dofs =
- std::vector<types::global_dof_index>(tria.n_raw_quads() *
- fe.n_dofs_per_quad(),
- numbers::invalid_dof_index);
+
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(fe.n_unique_faces(), 1);
+ dof_handler.mg_faces->quads.dofs = std::vector<types::global_dof_index>(
+ tria.n_raw_quads() * fe.n_dofs_per_quad(0 /*=face_no*/),
+ numbers::invalid_dof_index);
const unsigned int n_vertices = tria.n_vertices();
fe_slots_needed = 1;
n_face_slots +=
dof_handler.get_fe(cell->active_fe_index())
- .template n_dofs_per_object<dim - 1>();
+ .template n_dofs_per_object<dim - 1>(face);
}
else
{
fe_slots_needed = 2;
n_face_slots +=
dof_handler.get_fe(cell->active_fe_index())
- .template n_dofs_per_object<dim - 1>() +
+ .template n_dofs_per_object<dim - 1>(face) +
dof_handler
.get_fe(cell->neighbor(face)->active_fe_index())
- .template n_dofs_per_object<dim - 1>();
+ .template n_dofs_per_object<dim - 1>(
+ cell->neighbor_face_no(face));
}
// mark this face as visited
const unsigned int fe = cell->active_fe_index();
const unsigned int n_dofs =
dof_handler.get_fe(fe)
- .template n_dofs_per_object<dim - 1>();
+ .template n_dofs_per_object<dim - 1>(face);
const unsigned int offset =
dof_handler
.hp_object_fe_ptr[d][cell->face(face)->index()];
}
else
{
- unsigned int fe_1 = cell->active_fe_index();
+ unsigned int fe_1 = cell->active_fe_index();
+ unsigned int face_no_1 = face;
unsigned int fe_2 =
cell->neighbor(face)->active_fe_index();
+ unsigned int face_no_2 = cell->neighbor_face_no(face);
if (fe_2 < fe_1)
- std::swap(fe_1, fe_2);
+ {
+ std::swap(fe_1, fe_2);
+ std::swap(face_no_1, face_no_2);
+ }
const unsigned int n_dofs_1 =
dof_handler.get_fe(fe_1)
- .template n_dofs_per_object<dim - 1>();
+ .template n_dofs_per_object<dim - 1>(face_no_1);
const unsigned int n_dofs_2 =
dof_handler.get_fe(fe_2)
- .template n_dofs_per_object<dim - 1>();
+ .template n_dofs_per_object<dim - 1>(face_no_2);
const unsigned int offset =
dof_handler
if (face->at_boundary())
{
const unsigned int dofs_per_face =
- cell->get_fe().n_dofs_per_face();
+ cell->get_fe().n_dofs_per_face(iface);
dofs_on_face.resize(dofs_per_face);
face->get_dof_indices(dofs_on_face, cell->active_fe_index());
(boundary_ids.find(boundary_id) != boundary_ids.end()))
{
const unsigned int dofs_per_face =
- cell->get_fe().n_dofs_per_face();
+ cell->get_fe().n_dofs_per_face(iface);
dofs_on_face.resize(dofs_per_face);
face->get_dof_indices(dofs_on_face, cell->active_fe_index());
ensure_existence_and_return_dof_identities(
const FiniteElement<dim, spacedim> &fe1,
const FiniteElement<dim, spacedim> &fe2,
- std::unique_ptr<DoFIdentities> & identities)
+ std::unique_ptr<DoFIdentities> & identities,
+ const unsigned int face_no = numbers::invalid_unsigned_int)
{
+ Assert(structdim == 2 || face_no == numbers::invalid_unsigned_int,
+ ExcInternalError());
+
// see if we need to fill this entry, or whether it already
// exists
if (identities.get() == nullptr)
case 2:
{
identities = std::make_unique<DoFIdentities>(
- fe1.hp_quad_dof_identities(fe2));
+ fe1.hp_quad_dof_identities(fe2, face_no));
break;
}
for (unsigned int i = 0; i < identities->size(); ++i)
{
Assert((*identities)[i].first <
- fe1.template n_dofs_per_object<structdim>(),
+ fe1.template n_dofs_per_object<structdim>(face_no),
ExcInternalError());
Assert((*identities)[i].second <
- fe2.template n_dofs_per_object<structdim>(),
+ fe2.template n_dofs_per_object<structdim>(face_no),
ExcInternalError());
}
}
fe_indices,
/*codim=*/dim - 2);
+ const unsigned int most_dominating_fe_index_face_no =
+ cell->active_fe_index() == most_dominating_fe_index ?
+ q :
+ cell->neighbor_face_no(q);
+
// if we found the most dominating element, then use
// this to eliminate some of the degrees of freedom
// by identification. otherwise, the code that
quad_dof_identities
[most_dominating_fe_index][other_fe_index]
[cell->quad(q)->reference_cell_type() ==
- ReferenceCell::Type::Quad]);
+ ReferenceCell::Type::Quad],
+ most_dominating_fe_index_face_no);
for (const auto &identity : identities)
{
fe_indices,
/*codim=*/dim - 2);
+ const unsigned int most_dominating_fe_index_face_no =
+ cell->active_fe_index() == most_dominating_fe_index ?
+ q :
+ cell->neighbor_face_no(q);
+
// if we found the most dominating element, then use
// this to eliminate some of the degrees of freedom
// by identification. otherwise, the code that
quad_dof_identities
[most_dominating_fe_index][other_fe_index]
[cell->quad(q)->reference_cell_type() ==
- ReferenceCell::Type::Quad]);
+ ReferenceCell::Type::Quad],
+ most_dominating_fe_index_face_no);
for (const auto &identity : identities)
{
for (unsigned int d = 0;
d <
- dof_handler.get_fe(fe_index).n_dofs_per_quad();
+ dof_handler.get_fe(fe_index).n_dofs_per_quad(q);
++d)
{
const types::global_dof_index old_dof_index =
if (cell->quad(l)->user_flag_set())
{
for (unsigned int d = 0;
- d < dof_handler.get_fe().n_dofs_per_quad();
+ d < dof_handler.get_fe().n_dofs_per_quad(l);
++d)
{
const dealii::types::global_dof_index idx =
{
const FiniteElement<dim, spacedim> &fe = cell->get_fe();
- const unsigned int dofs_per_face = fe.n_dofs_per_face();
+ const unsigned int dofs_per_face = fe.n_dofs_per_face(face);
face_dof_indices.resize(dofs_per_face);
cell->face(face)->get_dof_indices(face_dof_indices,
cell->active_fe_index());
- for (unsigned int i = 0; i < fe.n_dofs_per_face(); ++i)
+ for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
if (!check_vector_component)
selected_dofs.add_index(face_dof_indices[i]);
else
numbers::invalid_unsigned_int)));
if (fe.is_primitive(cell_index))
{
- if (component_mask
- [fe.face_system_to_component_index(i).first] ==
- true)
+ if (component_mask[fe.face_system_to_component_index(
+ i, face)
+ .first] == true)
selected_dofs.add_index(face_dof_indices[i]);
}
else // not primitive
{
// simply take all DoFs that live on this subface
std::vector<types::global_dof_index> ldi(
- fe.n_dofs_per_face());
+ fe.n_dofs_per_face(f, child));
face->child(child)->get_dof_indices(ldi);
selected_dofs.add_indices(ldi.begin(), ldi.end());
}
for (const unsigned int f : cell->face_indices())
if (cell->at_boundary(f))
{
- const unsigned int dofs_per_face = cell->get_fe().n_dofs_per_face();
+ const unsigned int dofs_per_face =
+ cell->get_fe().n_dofs_per_face(f);
dofs_on_face.resize(dofs_per_face);
cell->face(f)->get_dof_indices(dofs_on_face,
cell->active_fe_index());
if (boundary_ids.find(cell->face(f)->boundary_id()) !=
boundary_ids.end())
{
- const unsigned int dofs_per_face = cell->get_fe().n_dofs_per_face();
+ const unsigned int dofs_per_face =
+ cell->get_fe().n_dofs_per_face(f);
dofs_on_face.resize(dofs_per_face);
cell->face(f)->get_dof_indices(dofs_on_face,
cell->active_fe_index());
// Exclude degrees of freedom on faces opposite to the vertex
exclude.resize(fe.n_dofs_per_cell());
std::fill(exclude.begin(), exclude.end(), false);
- const unsigned int dpf = fe.n_dofs_per_face();
for (const unsigned int face : cell->face_indices())
if (cell->at_boundary(face) ||
cell->neighbor(face)->level() != cell->level())
- for (unsigned int i = 0; i < dpf; ++i)
+ for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
exclude[fe.face_to_cell_index(i, face)] = true;
for (types::global_dof_index j = 0; j < indices.size(); ++j)
if (!exclude[j])
{
// Eliminate dofs on faces of the child which are on faces
// of the parent
- const unsigned int dpf = fe.n_dofs_per_face();
-
for (unsigned int d = 0; d < dim; ++d)
{
const unsigned int face =
GeometryInfo<dim>::vertex_to_face[child][d];
- for (unsigned int i = 0; i < dpf; ++i)
+ for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
exclude[fe.face_to_cell_index(i, face)] = true;
}
for (const unsigned int face :
GeometryInfo<dim>::face_indices())
if (cell->at_boundary(face))
- for (unsigned int i = 0; i < dpf; ++i)
+ for (unsigned int i = 0; i < fe.n_dofs_per_face(face);
+ ++i)
exclude[fe.face_to_cell_index(i, face)] = false;
}
// vertex
exclude.resize(fe.n_dofs_per_cell());
std::fill(exclude.begin(), exclude.end(), false);
- const unsigned int dpf = fe.n_dofs_per_face();
for (unsigned int d = 0; d < dim; ++d)
{
GeometryInfo<dim>::vertex_to_face[v][d];
const unsigned int face =
GeometryInfo<dim>::opposite_face[a_face];
- for (unsigned int i = 0; i < dpf; ++i)
+ for (unsigned int i = 0; i < fe.n_dofs_per_face(face); ++i)
{
// For each dof, get the block it is in and decide to
// exclude it or not
const FullMatrix<double> & face_interpolation_matrix,
std::vector<bool> & primary_dof_mask)
{
- Assert(fe1.n_dofs_per_face() >= fe2.n_dofs_per_face(),
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(fe1.n_unique_faces(), 1);
+ AssertDimension(fe2.n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+ (void)face_no;
+
+ Assert(fe1.n_dofs_per_face(face_no) >= fe2.n_dofs_per_face(face_no),
ExcInternalError());
- AssertDimension(primary_dof_mask.size(), fe1.n_dofs_per_face());
+ AssertDimension(primary_dof_mask.size(), fe1.n_dofs_per_face(face_no));
Assert(fe2.n_dofs_per_vertex() <= fe1.n_dofs_per_vertex(),
ExcInternalError());
Assert(fe2.n_dofs_per_line() <= fe1.n_dofs_per_line(),
ExcInternalError());
- Assert((dim < 3) || (fe2.n_dofs_per_quad() <= fe1.n_dofs_per_quad()),
+ Assert((dim < 3) ||
+ (fe2.n_dofs_per_quad(face_no) <= fe1.n_dofs_per_quad(face_no)),
ExcInternalError());
// the idea here is to designate as many DoFs in fe1 per object (vertex,
// same algorithm as above
unsigned int dofs_added = 0;
unsigned int i = 0;
- while (dofs_added < fe2.n_dofs_per_quad())
+ while (dofs_added < fe2.n_dofs_per_quad(q))
{
- Assert(i < fe1.n_dofs_per_quad(), ExcInternalError());
+ Assert(i < fe1.n_dofs_per_quad(q), ExcInternalError());
primary_dof_list.push_back(index + i);
if (check_primary_dof_list(face_interpolation_matrix,
++i;
}
- index += fe1.n_dofs_per_quad();
+ index += fe1.n_dofs_per_quad(q);
}
- AssertDimension(index, fe1.n_dofs_per_face());
- AssertDimension(primary_dof_list.size(), fe2.n_dofs_per_face());
+ AssertDimension(index, fe1.n_dofs_per_face(face_no));
+ AssertDimension(primary_dof_list.size(), fe2.n_dofs_per_face(face_no));
// finally copy the list into the mask
std::fill(primary_dof_mask.begin(), primary_dof_mask.end(), false);
const FullMatrix<double> & face_interpolation_matrix,
std::unique_ptr<std::vector<bool>> &primary_dof_mask)
{
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(fe1.n_unique_faces(), 1);
+ AssertDimension(fe2.n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
if (primary_dof_mask == nullptr)
{
primary_dof_mask =
- std::make_unique<std::vector<bool>>(fe1.n_dofs_per_face());
+ std::make_unique<std::vector<bool>>(fe1.n_dofs_per_face(face_no));
select_primary_dofs_for_face_restriction(fe1,
fe2,
face_interpolation_matrix,
const FiniteElement<dim, spacedim> & fe2,
std::unique_ptr<FullMatrix<double>> &matrix)
{
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(fe1.n_unique_faces(), 1);
+ AssertDimension(fe2.n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
if (matrix == nullptr)
{
- matrix =
- std::make_unique<FullMatrix<double>>(fe2.n_dofs_per_face(),
- fe1.n_dofs_per_face());
- fe1.get_face_interpolation_matrix(fe2, *matrix);
+ matrix = std::make_unique<FullMatrix<double>>(
+ fe2.n_dofs_per_face(face_no), fe1.n_dofs_per_face(face_no));
+ fe1.get_face_interpolation_matrix(fe2, *matrix, face_no);
}
}
const unsigned int subface,
std::unique_ptr<FullMatrix<double>> &matrix)
{
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(fe1.n_unique_faces(), 1);
+ AssertDimension(fe2.n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
if (matrix == nullptr)
{
- matrix =
- std::make_unique<FullMatrix<double>>(fe2.n_dofs_per_face(),
- fe1.n_dofs_per_face());
- fe1.get_subface_interpolation_matrix(fe2, subface, *matrix);
+ matrix = std::make_unique<FullMatrix<double>>(
+ fe2.n_dofs_per_face(face_no), fe1.n_dofs_per_face(face_no));
+ fe1.get_subface_interpolation_matrix(fe2,
+ subface,
+ *matrix,
+ face_no);
}
}
// first of all, make sure that we treat a case which is
// possible, i.e. either no dofs on the face at all or no
// anisotropic refinement
- if (cell->get_fe().n_dofs_per_face() == 0)
+ if (cell->get_fe().n_dofs_per_face(face) == 0)
continue;
Assert(cell->face(face)->refinement_case() ==
const FiniteElement<dim> &fe = cell->get_fe();
const unsigned int fe_index = cell->active_fe_index();
- const unsigned int n_dofs_on_mother = fe.n_dofs_per_face();
+ const unsigned int n_dofs_on_mother = fe.n_dofs_per_face(face);
const unsigned int n_dofs_on_children =
(5 * fe.n_dofs_per_vertex() + 12 * fe.n_dofs_per_line() +
- 4 * fe.n_dofs_per_quad());
+ 4 * fe.n_dofs_per_quad(face));
// TODO[TL]: think about this and the following in case of
// anisotropic refinement
for (unsigned int dof = 0; dof != fe.n_dofs_per_line(); ++dof)
dofs_on_mother[next_index++] =
this_face->line(line)->dof_index(dof, fe_index);
- for (unsigned int dof = 0; dof != fe.n_dofs_per_quad(); ++dof)
+ for (unsigned int dof = 0; dof != fe.n_dofs_per_quad(face);
+ ++dof)
dofs_on_mother[next_index++] =
this_face->dof_index(dof, fe_index);
AssertDimension(next_index, dofs_on_mother.size());
if (cell->neighbor_child_on_subface(face, child)
->is_artificial())
continue;
- for (unsigned int dof = 0; dof != fe.n_dofs_per_quad();
+ for (unsigned int dof = 0; dof != fe.n_dofs_per_quad(face);
++dof)
dofs_on_children.push_back(
this_face->child(child)->dof_index(dof, fe_index));
// first of all, make sure that we treat a case which is
// possible, i.e. either no dofs on the face at all or no
// anisotropic refinement
- if (cell->get_fe().n_dofs_per_face() == 0)
+ if (cell->get_fe().n_dofs_per_face(face) == 0)
continue;
Assert(cell->face(face)->refinement_case() ==
//
// so we are going to constrain the DoFs on the face
// children against the DoFs on the face itself
- primary_dofs.resize(cell->get_fe().n_dofs_per_face());
+ primary_dofs.resize(
+ cell->get_fe().n_dofs_per_face(face));
cell->face(face)->get_dof_indices(
primary_dofs, cell->active_fe_index());
// the face DoFs from the cell DoFs.
dependent_dofs.resize(
subface->get_fe(subface_fe_index)
- .n_dofs_per_face());
+ .n_dofs_per_face(face, c));
subface->get_dof_indices(dependent_dofs,
subface_fe_index);
// first get the interpolation matrix from the mother to
// the virtual dofs
- Assert(dominating_fe.n_dofs_per_face() <=
- cell->get_fe().n_dofs_per_face(),
+ Assert(dominating_fe.n_dofs_per_face(face) <=
+ cell->get_fe().n_dofs_per_face(face),
ExcInternalError());
ensure_existence_of_face_matrix(
// now compute the constraint matrix as the product
// between the inverse matrix and the dependent part
constraint_matrix.reinit(
- cell->get_fe().n_dofs_per_face() -
- dominating_fe.n_dofs_per_face(),
- dominating_fe.n_dofs_per_face());
+ cell->get_fe().n_dofs_per_face(face) -
+ dominating_fe.n_dofs_per_face(face),
+ dominating_fe.n_dofs_per_face(face));
restrict_mother_to_virtual_dependent.mmult(
constraint_matrix,
restrict_mother_to_virtual_primary_inv);
// then figure out the global numbers of primary and
// dependent dofs and apply constraints
- scratch_dofs.resize(cell->get_fe().n_dofs_per_face());
+ scratch_dofs.resize(
+ cell->get_fe().n_dofs_per_face(face));
cell->face(face)->get_dof_indices(
scratch_dofs, cell->active_fe_index());
primary_dofs.clear();
dependent_dofs.clear();
for (unsigned int i = 0;
- i < cell->get_fe().n_dofs_per_face();
+ i < cell->get_fe().n_dofs_per_face(face);
++i)
if ((*primary_dof_masks[dominating_fe_index]
[cell
dependent_dofs.push_back(scratch_dofs[i]);
AssertDimension(primary_dofs.size(),
- dominating_fe.n_dofs_per_face());
+ dominating_fe.n_dofs_per_face(face));
AssertDimension(dependent_dofs.size(),
- cell->get_fe().n_dofs_per_face() -
- dominating_fe.n_dofs_per_face());
+ cell->get_fe().n_dofs_per_face(face) -
+ dominating_fe.n_dofs_per_face(face));
filter_constraints(primary_dofs,
dependent_dofs,
// first get the interpolation matrix from the
// subface to the virtual dofs
- Assert(dominating_fe.n_dofs_per_face() <=
- subface_fe.n_dofs_per_face(),
+ Assert(dominating_fe.n_dofs_per_face(face) <=
+ subface_fe.n_dofs_per_face(face),
ExcInternalError());
ensure_existence_of_subface_matrix(
dominating_fe,
[dominating_fe_index][subface_fe_index][sf]);
constraint_matrix.reinit(
- subface_fe.n_dofs_per_face(),
- dominating_fe.n_dofs_per_face());
+ subface_fe.n_dofs_per_face(face),
+ dominating_fe.n_dofs_per_face(face));
restrict_subface_to_virtual.mmult(
constraint_matrix,
restrict_mother_to_virtual_primary_inv);
- dependent_dofs.resize(subface_fe.n_dofs_per_face());
+ dependent_dofs.resize(
+ subface_fe.n_dofs_per_face(face));
cell->face(face)->child(sf)->get_dof_indices(
dependent_dofs, subface_fe_index);
// Get DoFs on dominating and dominated side of the
// face
primary_dofs.resize(
- cell->get_fe().n_dofs_per_face());
+ cell->get_fe().n_dofs_per_face(face));
cell->face(face)->get_dof_indices(
primary_dofs, cell->active_fe_index());
break;
dependent_dofs.resize(
- neighbor->get_fe().n_dofs_per_face());
+ neighbor->get_fe().n_dofs_per_face(face));
cell->face(face)->get_dof_indices(
dependent_dofs, neighbor->active_fe_index());
// first get the interpolation matrix from main FE
// to the virtual dofs
- Assert(dominating_fe.n_dofs_per_face() <=
- cell->get_fe().n_dofs_per_face(),
+ Assert(dominating_fe.n_dofs_per_face(face) <=
+ cell->get_fe().n_dofs_per_face(face),
ExcInternalError());
ensure_existence_of_face_matrix(
// now compute the constraint matrix as the product
// between the inverse matrix and the dependent part
constraint_matrix.reinit(
- cell->get_fe().n_dofs_per_face() -
- dominating_fe.n_dofs_per_face(),
- dominating_fe.n_dofs_per_face());
+ cell->get_fe().n_dofs_per_face(face) -
+ dominating_fe.n_dofs_per_face(face),
+ dominating_fe.n_dofs_per_face(face));
restrict_mother_to_virtual_dependent.mmult(
constraint_matrix,
restrict_mother_to_virtual_primary_inv);
// then figure out the global numbers of primary and
// dependent dofs and apply constraints
scratch_dofs.resize(
- cell->get_fe().n_dofs_per_face());
+ cell->get_fe().n_dofs_per_face(face));
cell->face(face)->get_dof_indices(
scratch_dofs, cell->active_fe_index());
primary_dofs.clear();
dependent_dofs.clear();
for (unsigned int i = 0;
- i < cell->get_fe().n_dofs_per_face();
+ i < cell->get_fe().n_dofs_per_face(face);
++i)
if ((*primary_dof_masks[dominating_fe_index]
[cell->active_fe_index()])
dependent_dofs.push_back(scratch_dofs[i]);
AssertDimension(primary_dofs.size(),
- dominating_fe.n_dofs_per_face());
- AssertDimension(dependent_dofs.size(),
- cell->get_fe().n_dofs_per_face() -
- dominating_fe.n_dofs_per_face());
+ dominating_fe.n_dofs_per_face(
+ face));
+ AssertDimension(
+ dependent_dofs.size(),
+ cell->get_fe().n_dofs_per_face(face) -
+ dominating_fe.n_dofs_per_face(face));
filter_constraints(primary_dofs,
dependent_dofs,
// now do the same for another FE this is pretty
// much the same we do above to resolve h-refinement
// constraints
- Assert(dominating_fe.n_dofs_per_face() <=
- neighbor->get_fe().n_dofs_per_face(),
+ Assert(dominating_fe.n_dofs_per_face(face) <=
+ neighbor->get_fe().n_dofs_per_face(face),
ExcInternalError());
ensure_existence_of_face_matrix(
[neighbor->active_fe_index()]);
constraint_matrix.reinit(
- neighbor->get_fe().n_dofs_per_face(),
- dominating_fe.n_dofs_per_face());
+ neighbor->get_fe().n_dofs_per_face(face),
+ dominating_fe.n_dofs_per_face(face));
restrict_secondface_to_virtual.mmult(
constraint_matrix,
restrict_mother_to_virtual_primary_inv);
dependent_dofs.resize(
- neighbor->get_fe().n_dofs_per_face());
+ neighbor->get_fe().n_dofs_per_face(face));
cell->face(face)->get_dof_indices(
dependent_dofs, neighbor->active_fe_index());
Assert(face_1->n_active_fe_indices() == 1, ExcInternalError());
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(
+ face_1->get_fe(face_1->nth_active_fe_index(0)).n_unique_faces(), 1);
+ AssertDimension(
+ face_2->get_fe(face_2->nth_active_fe_index(0)).n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
// If face_2 does have children, then we need to iterate over these
// children and set periodic constraints in the inverse direction:
ExcNotImplemented());
const unsigned int dofs_per_face =
- face_1->get_fe(face_1->nth_active_fe_index(0)).n_dofs_per_face();
+ face_1->get_fe(face_1->nth_active_fe_index(0))
+ .n_dofs_per_face(face_no);
FullMatrix<double> child_transformation(dofs_per_face, dofs_per_face);
FullMatrix<double> subface_interpolation(dofs_per_face,
dofs_per_face);
// interpolated from face_1 to face_2 by multiplying from the left
// with the one that interpolates from face_2 to its child
const auto &fe = face_1->get_fe(face_1->nth_active_fe_index(0));
- fe.get_subface_interpolation_matrix(fe, c, subface_interpolation);
+ fe.get_subface_interpolation_matrix(fe,
+ c,
+ subface_interpolation,
+ face_no);
subface_interpolation.mmult(child_transformation, transformation);
set_periodicity_constraints(face_1,
"zero or equal to the number of components in the finite "
"element."));
- const unsigned int dofs_per_face = fe.n_dofs_per_face();
+ const unsigned int dofs_per_face = fe.n_dofs_per_face(face_no);
std::vector<types::global_dof_index> dofs_1(dofs_per_face);
std::vector<types::global_dof_index> dofs_2(dofs_per_face);
// Obey the component mask
if ((component_mask.n_selected_components(fe.n_components()) !=
fe.n_components()) &&
- !component_mask[fe.face_system_to_component_index(i).first])
+ !component_mask[fe.face_system_to_component_index(i, face_no)
+ .first])
continue;
// We have to be careful to treat so called "identity
const FullMatrix<double> & matrix,
const std::vector<unsigned int> & first_vector_components)
{
+ // 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;
+
Assert(matrix.m() == matrix.n(), ExcInternalError());
- const unsigned int n_dofs_per_face = fe.n_dofs_per_face();
+ const unsigned int n_dofs_per_face = fe.n_dofs_per_face(face_no);
if (matrix.m() == n_dofs_per_face)
{
Assert(matrix.m() == spacedim, ExcInternalError())
Quadrature<dim - 1>
- quadrature(fe.get_unit_face_support_points());
+ quadrature(fe.get_unit_face_support_points(face_no));
// have an array that stores the location of each vector-dof tuple we want
// to rotate.
std::vector<unsigned int>::const_iterator comp_it =
std::find(first_vector_components.begin(),
first_vector_components.end(),
- fe.face_system_to_component_index(i).first);
+ fe.face_system_to_component_index(i, face_no).first);
if (comp_it != first_vector_components.end())
{
const unsigned int first_vector_component = *comp_it;
for (unsigned int k = 0; k < n_dofs_per_face; ++k)
if ((k != i) && (quadrature.point(k) == quadrature.point(i)) &&
- (fe.face_system_to_component_index(k).first >=
+ (fe.face_system_to_component_index(k, face_no).first >=
first_vector_component) &&
- (fe.face_system_to_component_index(k).first <
+ (fe.face_system_to_component_index(k, face_no).first <
first_vector_component + spacedim))
{
- vector_dofs[fe.face_system_to_component_index(k).first -
+ vector_dofs[fe.face_system_to_component_index(k, face_no)
+ .first -
first_vector_component] = k;
n_found++;
if (n_found == dim)
const std::vector<unsigned int> & first_vector_components,
const number periodicity_factor)
{
+ // TODO: the implementation makes the assumption that all faces have the
+ // same number of dofs
+ AssertDimension(
+ face_1->get_fe(face_1->nth_active_fe_index(0)).n_unique_faces(), 1);
+ AssertDimension(
+ face_2->get_fe(face_2->nth_active_fe_index(0)).n_unique_faces(), 1);
+ const unsigned int face_no = 0;
+
static const int dim = FaceIterator::AccessorType::dimension;
static const int spacedim = FaceIterator::AccessorType::space_dimension;
{
Assert(face_1->n_active_fe_indices() == 1, ExcInternalError());
const unsigned int n_dofs_per_face =
- face_1->get_fe(face_1->nth_active_fe_index(0)).n_dofs_per_face();
+ face_1->get_fe(face_1->nth_active_fe_index(0))
+ .n_dofs_per_face(face_no);
Assert(matrix.m() == 0 ||
(first_vector_components.empty() &&
{
Assert(face_2->n_active_fe_indices() == 1, ExcInternalError());
const unsigned int n_dofs_per_face =
- face_2->get_fe(face_2->nth_active_fe_index(0)).n_dofs_per_face();
+ face_2->get_fe(face_2->nth_active_fe_index(0))
+ .n_dofs_per_face(face_no);
Assert(matrix.m() == 0 ||
(first_vector_components.empty() &&
face_2->get_fe(face_2->nth_active_fe_index(0)) :
face_1->get_fe(face_1->nth_active_fe_index(0));
- const unsigned int n_dofs_per_face = fe.n_dofs_per_face();
+ const unsigned int n_dofs_per_face = fe.n_dofs_per_face(face_no);
// Sometimes we just have nothing to do (for all finite elements, or
// systems which accidentally don't have any dofs on the boundary).
(face->boundary_id() == boundary_id)))
{
// get indices and physical location on this face
- face_dofs.resize(fe.n_dofs_per_face());
+ face_dofs.resize(fe.n_dofs_per_face(face_no));
face->get_dof_indices(face_dofs, cell->active_fe_index());
// enter those dofs into the list that match the component
for (const unsigned int f : cell->face_indices())
if (cell->at_boundary(f))
{
- const unsigned int dofs_per_face = cell->get_fe().n_dofs_per_face();
+ const unsigned int dofs_per_face =
+ cell->get_fe().n_dofs_per_face(f);
dofs_on_this_face.resize(dofs_per_face);
cell->face(f)->get_dof_indices(dofs_on_this_face,
cell->active_fe_index());
if (boundary_ids.find(cell->face(f)->boundary_id()) !=
boundary_ids.end())
{
- const unsigned int dofs_per_face = cell->get_fe().n_dofs_per_face();
+ const unsigned int dofs_per_face =
+ cell->get_fe().n_dofs_per_face(f);
dofs_on_this_face.resize(dofs_per_face);
cell->face(f)->get_dof_indices(dofs_on_this_face,
cell->active_fe_index());
ExcInternalError());
const unsigned int n_dofs_per_face =
- cell->get_fe().n_dofs_per_face();
+ cell->get_fe().n_dofs_per_face(f, c);
local_face_dof_indices.resize(n_dofs_per_face);
cell->face(f)->child(c)->get_dof_indices(
unsigned int subface = neighbor_face_no_subface_no.second;
const unsigned int n_dofs_per_face =
- cell->get_fe().n_dofs_per_face();
+ cell->get_fe().n_dofs_per_face(face_no);
local_face_dof_indices.resize(n_dofs_per_face);
cell->neighbor(f)->face(face_no)->get_dof_indices(
// original cell
{
const unsigned int n_dofs_per_face =
- cell->get_fe().n_dofs_per_face();
+ cell->get_fe().n_dofs_per_face(face_no, c);
local_face_dof_indices.resize(n_dofs_per_face);
Assert(cell->neighbor(f)
// n_dofs_per_vertex() is
// arbitrary, not the other way
// round.
- // TODO: This assumes that even in hp context, the dofs per face
- // coincide!
+ // TODO: This assumes that the dofs per face on all faces coincide!
+ const unsigned int face_no = 0;
+
+ Assert(fe.reference_cell_type() == ReferenceCell::get_hypercube(dim),
+ ExcNotImplemented());
+
unsigned int increment =
- fe.n_dofs_per_cell() - dim * fe.n_dofs_per_face();
+ fe.n_dofs_per_cell() - dim * fe.n_dofs_per_face(face_no);
while (i < fe.get_first_line_index())
row_lengths[cell_indices[i++]] += increment;
// From now on, if an object is
// In all other cases we
// subtract adjacent faces to be
// added in the loop below.
- increment = (dim > 1) ?
- fe.n_dofs_per_cell() - (dim - 1) * fe.n_dofs_per_face() :
- fe.n_dofs_per_cell() - GeometryInfo<dim>::faces_per_cell *
- fe.n_dofs_per_face();
- while (i < fe.get_first_quad_index())
+ increment =
+ (dim > 1) ?
+ fe.n_dofs_per_cell() - (dim - 1) * fe.n_dofs_per_face(face_no) :
+ fe.n_dofs_per_cell() -
+ GeometryInfo<dim>::faces_per_cell * fe.n_dofs_per_face(face_no);
+ while (i < fe.get_first_quad_index(face_no))
row_lengths[cell_indices[i++]] += increment;
// Now quads in 2D and 3D
- increment = (dim > 2) ?
- fe.n_dofs_per_cell() - (dim - 2) * fe.n_dofs_per_face() :
- fe.n_dofs_per_cell() - GeometryInfo<dim>::faces_per_cell *
- fe.n_dofs_per_face();
+ increment =
+ (dim > 2) ?
+ fe.n_dofs_per_cell() - (dim - 2) * fe.n_dofs_per_face(face_no) :
+ fe.n_dofs_per_cell() -
+ GeometryInfo<dim>::faces_per_cell * fe.n_dofs_per_face(face_no);
while (i < fe.get_first_hex_index())
row_lengths[cell_indices[i++]] += increment;
// Finally, cells in 3D
- increment = fe.n_dofs_per_cell() -
- GeometryInfo<dim>::faces_per_cell * fe.n_dofs_per_face();
+ increment = fe.n_dofs_per_cell() - GeometryInfo<dim>::faces_per_cell *
+ fe.n_dofs_per_face(face_no);
while (i < fe.n_dofs_per_cell())
row_lengths[cell_indices[i++]] += increment;
for (unsigned int local_dof = 0;
local_dof < fe.n_dofs_per_cell();
++local_dof)
- row_lengths[cell_indices[local_dof]] += fe.n_dofs_per_face();
+ row_lengths[cell_indices[local_dof]] +=
+ fe.n_dofs_per_face(face_no);
continue;
}
if (flux_coupling != DoFTools::none)
{
const unsigned int dof_increment =
- nfe.n_dofs_per_cell() - nfe.n_dofs_per_face();
+ nfe.n_dofs_per_cell() - nfe.n_dofs_per_face(face_no);
for (unsigned int local_dof = 0;
local_dof < fe.n_dofs_per_cell();
++local_dof)
neighbor->get_mg_dof_indices(neighbor_indices);
for (unsigned int local_dof = 0; local_dof < fe.n_dofs_per_cell();
++local_dof)
- row_lengths[cell_indices[local_dof]] += nfe.n_dofs_per_face();
+ row_lengths[cell_indices[local_dof]] +=
+ nfe.n_dofs_per_face(face_no);
for (unsigned int local_dof = 0; local_dof < nfe.n_dofs_per_cell();
++local_dof)
- row_lengths[neighbor_indices[local_dof]] += fe.n_dofs_per_face();
+ row_lengths[neighbor_indices[local_dof]] +=
+ fe.n_dofs_per_face(face_no);
}
}
user_flags_triangulation.load_user_flags(old_flags);
const FiniteElement<dim> &fe = cell->get_fe();
const unsigned int fe_index = cell->active_fe_index();
+
+ // TODO: This assumes that the dofs per face on all faces coincide!
+ const unsigned int face_no = 0;
+ Assert(fe.reference_cell_type() == ReferenceCell::get_hypercube(dim),
+ ExcNotImplemented());
+
Assert(couplings.n_rows() == fe.n_components(),
ExcDimensionMismatch(couplings.n_rows(), fe.n_components()));
Assert(couplings.n_cols() == fe.n_components(),
fe.first_block_of_base(base) +
mult) != DoFTools::none)
{
- increment = fe.base_element(base).n_dofs_per_cell() -
- dim * fe.base_element(base).n_dofs_per_face();
+ increment =
+ fe.base_element(base).n_dofs_per_cell() -
+ dim * fe.base_element(base).n_dofs_per_face(face_no);
row_lengths[cell_indices[i]] += increment;
}
++i;
// In all other cases we
// subtract adjacent faces to be
// added in the loop below.
- while (i < fe.get_first_quad_index())
+ while (i < fe.get_first_quad_index(face_no))
{
for (unsigned int base = 0; base < fe.n_base_elements(); ++base)
for (unsigned int mult = 0; mult < fe.element_multiplicity(base);
fe.base_element(base).n_dofs_per_cell() -
((dim > 1) ? (dim - 1) :
GeometryInfo<dim>::faces_per_cell) *
- fe.base_element(base).n_dofs_per_face();
+ fe.base_element(base).n_dofs_per_face(face_no);
row_lengths[cell_indices[i]] += increment;
}
++i;
fe.base_element(base).n_dofs_per_cell() -
((dim > 2) ? (dim - 2) :
GeometryInfo<dim>::faces_per_cell) *
- fe.base_element(base).n_dofs_per_face();
+ fe.base_element(base).n_dofs_per_face(face_no);
row_lengths[cell_indices[i]] += increment;
}
++i;
fe.first_block_of_base(base) +
mult) != DoFTools::none)
{
- increment = fe.base_element(base).n_dofs_per_cell() -
- GeometryInfo<dim>::faces_per_cell *
- fe.base_element(base).n_dofs_per_face();
+ increment =
+ fe.base_element(base).n_dofs_per_cell() -
+ GeometryInfo<dim>::faces_per_cell *
+ fe.base_element(base).n_dofs_per_face(face_no);
row_lengths[cell_indices[i]] += increment;
}
++i;
for (unsigned int local_dof = 0;
local_dof < fe.n_dofs_per_cell();
++local_dof)
- row_lengths[cell_indices[local_dof]] += fe.n_dofs_per_face();
+ row_lengths[cell_indices[local_dof]] +=
+ fe.n_dofs_per_face(face_no);
continue;
}
{
const unsigned int dof_increment =
nfe.base_element(base).n_dofs_per_cell() -
- nfe.base_element(base).n_dofs_per_face();
+ nfe.base_element(base).n_dofs_per_face(face_no);
row_lengths[cell_indices[local_dof]] += dof_increment;
}
fe.system_to_component_index(local_dof).first,
nfe.first_block_of_base(base) + mult) != DoFTools::none)
row_lengths[cell_indices[local_dof]] +=
- nfe.base_element(base).n_dofs_per_face();
+ nfe.base_element(base).n_dofs_per_face(face_no);
for (unsigned int base = 0; base < fe.n_base_elements(); ++base)
for (unsigned int mult = 0; mult < fe.element_multiplicity(base);
++mult)
nfe.system_to_component_index(local_dof).first,
fe.first_block_of_base(base) + mult) != DoFTools::none)
row_lengths[neighbor_indices[local_dof]] +=
- fe.base_element(base).n_dofs_per_face();
+ fe.base_element(base).n_dofs_per_face(face_no);
}
}
user_flags_triangulation.load_user_flags(old_flags);
continue;
const FiniteElement<dim> &fe = cell->get_fe();
const unsigned int level = cell->level();
- local_dofs.resize(fe.n_dofs_per_face());
for (const unsigned int face_no : GeometryInfo<dim>::face_indices())
if (cell->at_boundary(face_no) == true)
// Face is listed in boundary map
if (boundary_ids.find(bi) != boundary_ids.end())
{
+ local_dofs.resize(fe.n_dofs_per_face(face_no));
face->get_mg_dof_indices(level, local_dofs);
dofs_by_level[level].insert(dofs_by_level[level].end(),
local_dofs.begin(),
// get indices, physical location and boundary values of
// dofs on this face
- local_dofs.resize(fe.n_dofs_per_face());
+ local_dofs.resize(fe.n_dofs_per_face(face_no));
face->get_mg_dof_indices(level, local_dofs);
if (fe_is_system)
{
numbers::invalid_unsigned_int;
if (fe.is_primitive())
component =
- fe.face_system_to_component_index(i).first;
+ fe.face_system_to_component_index(i, face_no)
+ .first;
else
{
// Just pick the first of the components
const FiniteElement<dim, spacedim> &fe = mg_dof_handler.get_fe();
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
- const unsigned int dofs_per_face = fe.n_dofs_per_face();
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
// Do refinement face from the coarse side
if (neighbor->level() < cell->level())
{
- for (unsigned int j = 0; j < dofs_per_face; ++j)
+ for (unsigned int j = 0; j < fe.n_dofs_per_face(face_nr);
+ ++j)
cell_dofs[fe.face_to_cell_index(j, face_nr)] = true;
has_coarser_neighbor = true;
projects / dealii.git / commitdiff