hp::FEValues<3>& hp_fe_values,
const Function<3>& boundary_function,
const unsigned int first_vector_component,
- std::vector<double>& dof_values)
+ std::vector<double>& dof_values,
+ std::vector<bool>& dofs_processed)
{
const double tol = 0.5 * cell->get_fe ().degree * 1e-13 / cell->face (face)->line (line)->diameter ();
const unsigned int dim = 3;
// objects.
const FEValues<dim>&
fe_values = hp_fe_values.get_present_fe_values ();
+ const FiniteElement<dim>& fe = cell->get_fe ();
const std::vector< DerivativeForm<1,dim,spacedim> > &
jacobians = fe_values.get_jacobians ();
const std::vector<Point<dim> >&
std::vector<Point<dim> > tangentials (fe_values.n_quadrature_points);
std::vector<Vector<double> > values (fe_values.n_quadrature_points,
- Vector<double> (cell->get_fe ().n_components ()));
+ Vector<double> (fe.n_components ()));
// Get boundary function values
// at quadrature points.
const std::vector<Point<dim> >&
reference_quadrature_points = fe_values.get_quadrature ().get_points ();
- const unsigned int superdegree = cell->get_fe ().degree;
- const unsigned int degree = superdegree - 1;
+ const unsigned int degree = fe.degree - 1;
+ std::pair<unsigned int, unsigned int> base_indices (0, 0);
+
+ if (dynamic_cast<const FESystem<dim>*> (&cell->get_fe ()) != 0)
+ {
+ unsigned int fe_index = 0;
+ unsigned int fe_index_old = 0;
+ unsigned int i = 0;
+
+ for (; i < fe.n_base_elements (); ++i)
+ {
+ fe_index_old = fe_index;
+ fe_index += fe.element_multiplicity (i) * fe.base_element (i).n_components ();
+
+ if (fe_index >= first_vector_component)
+ break;
+ }
+
+ base_indices.first = i;
+ base_indices.second = (first_vector_component - fe_index_old) / fe.base_element (i).n_components ();
+ }
// coordinate directions of
// the edges of the face.
// Compute the degrees of
// freedom.
- for (unsigned int i = 0; i <= degree; ++i)
- dof_values[i + line * superdegree]
- += (fe_values.JxW (q_point)
- * (values[q_point] (first_vector_component) * tangentials[q_point] (0)
- + values[q_point] (first_vector_component + 1) * tangentials[q_point] (1)
- + values[q_point] (first_vector_component + 2) * tangentials[q_point] (2))
- * (fe_values[vec].value (cell->get_fe ().face_to_cell_index (i + line * superdegree,
- face),
- q_point)
- * tangentials[q_point])
- / std::sqrt (jacobians[q_point][0][edge_coordinate_direction[face][line]]
- * jacobians[q_point][0][edge_coordinate_direction[face][line]]
- + jacobians[q_point][1][edge_coordinate_direction[face][line]]
- * jacobians[q_point][1][edge_coordinate_direction[face][line]]
- + jacobians[q_point][2][edge_coordinate_direction[face][line]]
- * jacobians[q_point][2][edge_coordinate_direction[face][line]]));
+ for (unsigned int i = 0; i < fe.dofs_per_face; ++i)
+ if (((dynamic_cast<const FESystem<dim>*> (&fe) != 0)
+ && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).first == base_indices)
+ && (fe.base_element (base_indices.first).face_to_cell_index (line * fe.degree, face)
+ <= fe.system_to_base_index (fe.face_to_cell_index (i, face)).second)
+ && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).second
+ <= fe.base_element (base_indices.first).face_to_cell_index
+ ((line + 1) * fe.degree - 1, face)))
+ || ((dynamic_cast<const FE_Nedelec<dim>*> (&fe) != 0) && (line * fe.degree <= i)
+ && (i < (line + 1) * fe.degree)))
+ {
+ dof_values[i]
+ += (fe_values.JxW (q_point)
+ * (values[q_point] (first_vector_component) * tangentials[q_point] (0)
+ + values[q_point] (first_vector_component + 1) * tangentials[q_point] (1)
+ + values[q_point] (first_vector_component + 2) * tangentials[q_point] (2))
+ * (fe_values[vec].value (fe.face_to_cell_index (i, face), q_point) * tangentials[q_point])
+ / std::sqrt (jacobians[q_point][0][edge_coordinate_direction[face][line]]
+ * jacobians[q_point][0][edge_coordinate_direction[face][line]]
+ + jacobians[q_point][1][edge_coordinate_direction[face][line]]
+ * jacobians[q_point][1][edge_coordinate_direction[face][line]]
+ + jacobians[q_point][2][edge_coordinate_direction[face][line]]
+ * jacobians[q_point][2][edge_coordinate_direction[face][line]]));
+
+ if (q_point == 0)
+ dofs_processed[i] = true;
+ }
}
}
hp::FEValues<dim>&,
const Function<dim>&,
const unsigned int,
- std::vector<double>&)
+ std::vector<double>&,
+ std::vector<bool>&)
{
Assert (false, ExcInternalError ());
}
hp::FEValues<dim>& hp_fe_values,
const Function<dim>& boundary_function,
const unsigned int first_vector_component,
- std::vector<double>& dof_values)
+ std::vector<double>& dof_values,
+ std::vector<bool>& dofs_processed)
{
const unsigned int spacedim = dim;
hp_fe_values.reinit (cell, cell->active_fe_index ()
// objects.
const FEValues<dim>&
fe_values = hp_fe_values.get_present_fe_values ();
+ const FiniteElement<dim>& fe = cell->get_fe ();
const std::vector< DerivativeForm<1,dim,spacedim> > &
jacobians = fe_values.get_jacobians ();
+ const std::vector<Point<dim> >&
+ quadrature_points = fe_values.get_quadrature_points ();
+ const unsigned int degree = fe.degree - 1;
+ std::pair<unsigned int, unsigned int> base_indices (0, 0);
+
+ if (dynamic_cast<const FESystem<dim>*> (&cell->get_fe ()) != 0)
+ {
+ unsigned int fe_index = 0;
+ unsigned int fe_index_old = 0;
+ unsigned int i = 0;
+
+ for (; i < fe.n_base_elements (); ++i)
+ {
+ fe_index_old = fe_index;
+ fe_index += fe.element_multiplicity (i) * fe.base_element (i).n_components ();
+
+ if (fe_index >= first_vector_component)
+ break;
+ }
+
+ base_indices.first = i;
+ base_indices.second = (first_vector_component - fe_index_old) / fe.base_element (i).n_components ();
+ }
std::vector<Vector<double> >
- values (fe_values.n_quadrature_points, Vector<double> (cell->get_fe ().n_components ()));
+ values (fe_values.n_quadrature_points, Vector<double> (fe.n_components ()));
+
+ // Get boundary function
+ // values at quadrature
+ // points.
+ boundary_function.vector_value_list (quadrature_points, values);
switch (dim)
{
case 2:
{
- const std::vector<Point<dim> >&
- quadrature_points = fe_values.get_quadrature_points ();
std::vector<Point<dim> >
tangentials (fe_values.n_quadrature_points);
- // Get boundary function
- // values at quadrature
- // points.
- boundary_function.vector_value_list (quadrature_points, values);
-
const std::vector<Point<dim> >&
reference_quadrature_points = fe_values.get_quadrature ().get_points ();
- const unsigned int degree = cell->get_fe ().degree - 1;
// coordinate directions
// of the face.
/= std::sqrt (tangentials[q_point].square ());
// Compute the degrees
// of freedom.
- for (unsigned int i = 0; i <= degree; ++i)
- dof_values[i]
- += fe_values.JxW (q_point)
- * (values[q_point] (first_vector_component)
- * tangentials[q_point] (0)
- + values[q_point] (first_vector_component + 1)
- * tangentials[q_point] (1))
- * (fe_values[vec].value (cell->get_fe ().face_to_cell_index (i, face),
- q_point) * tangentials[q_point])
- / std::sqrt (jacobians[q_point][0][face_coordinate_direction[face]]
- * jacobians[q_point][0][face_coordinate_direction[face]]
- + jacobians[q_point][1][face_coordinate_direction[face]]
- * jacobians[q_point][1][face_coordinate_direction[face]]);
+ for (unsigned int i = 0; i < fe.dofs_per_face; ++i)
+ if (((dynamic_cast<const FESystem<dim>*> (&fe) != 0)
+ && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).first == base_indices))
+ || (dynamic_cast<const FE_Nedelec<dim>*> (&fe) != 0))
+ {
+ dof_values[i]
+ += fe_values.JxW (q_point)
+ * (values[q_point] (first_vector_component)
+ * tangentials[q_point] (0)
+ + values[q_point] (first_vector_component + 1)
+ * tangentials[q_point] (1))
+ * (fe_values[vec].value (fe.face_to_cell_index (i, face), q_point)
+ * tangentials[q_point])
+ / std::sqrt (jacobians[q_point][0][face_coordinate_direction[face]]
+ * jacobians[q_point][0][face_coordinate_direction[face]]
+ + jacobians[q_point][1][face_coordinate_direction[face]]
+ * jacobians[q_point][1][face_coordinate_direction[face]]);
+
+ if (q_point == 0)
+ dofs_processed[i] = true;
+ }
}
break;
case 3:
{
- const std::vector<Point<dim> >&
- quadrature_points = fe_values.get_quadrature_points ();
-
- // Get boundary function
- // values at quadrature
- // points.
- boundary_function.vector_value_list (quadrature_points, values);
-
const FEValuesExtractors::Vector vec (first_vector_component);
- const unsigned int superdegree = cell->get_fe ().degree;
- const unsigned int degree = superdegree - 1;
FullMatrix<double>
- assembling_matrix (degree * superdegree,
+ assembling_matrix (degree * fe.degree,
dim * fe_values.n_quadrature_points);
Vector<double> assembling_vector (assembling_matrix.n ());
Vector<double> cell_rhs (assembling_matrix.m ());
for (unsigned int d = 0; d < dim; ++d)
tmp[d] = values[q_point] (first_vector_component + d);
-
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j <= degree; ++j)
- tmp -= dof_values[(i + 2) * superdegree + j]
- * fe_values[vec].value (cell->get_fe ().face_to_cell_index
- ((i + 2) * superdegree + j,
- face), q_point);
+
+ for (unsigned int i = 0; i < fe.dofs_per_face; ++i)
+ if (((dynamic_cast<const FESystem<dim>*> (&fe) != 0)
+ && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).first == base_indices)
+ && (fe.base_element (base_indices.first).face_to_cell_index (2 * fe.degree, face)
+ <= fe.system_to_base_index (fe.face_to_cell_index (i, face)).second)
+ && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).second
+ < fe.base_element (base_indices.first).face_to_cell_index (4 * fe.degree, face)))
+ || ((dynamic_cast<const FE_Nedelec<dim>*> (&fe) != 0) && (2 * fe.degree <= i)
+ && (i < 4 * fe.degree)))
+ tmp -= dof_values[i] * fe_values[vec].value (fe.face_to_cell_index (i, face), q_point);
const double JxW
= std::sqrt (fe_values.JxW (q_point)
// the face.
for (unsigned int d = 0; d < dim; ++d)
assembling_vector (dim * q_point + d) = JxW * tmp[d];
-
- for (unsigned int i = 0; i <= degree; ++i)
- for (unsigned int j = 0; j < degree; ++j)
+
+ unsigned int index = 0;
+
+ for (unsigned int i = 0; i < fe.dofs_per_face; ++i)
+ if (((dynamic_cast<const FESystem<dim>*> (&fe) != 0)
+ && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).first == base_indices)
+ && (fe.base_element (base_indices.first).face_to_cell_index
+ (GeometryInfo<dim>::lines_per_face * fe.degree, face)
+ <= fe.system_to_base_index (fe.face_to_cell_index (i, face)).second)
+ && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).second
+ < fe.base_element (base_indices.first).face_to_cell_index
+ ((degree + GeometryInfo<dim>::lines_per_face) * fe.degree, face)))
+ || ((dynamic_cast<const FE_Nedelec<dim>*> (&fe) != 0)
+ && (GeometryInfo<dim>::lines_per_face * fe.degree <= i)
+ && (i < (degree + GeometryInfo<dim>::lines_per_face) * fe.degree)))
{
const Tensor<1, dim> shape_value
- = (JxW
- * fe_values[vec].value (cell->get_fe ()
- .face_to_cell_index
- ((i + GeometryInfo<dim>::lines_per_face)
- * degree
- + j
- + GeometryInfo<dim>::lines_per_face,
- face),
- q_point));
+ = (JxW * fe_values[vec].value (fe.face_to_cell_index (i, face),
+ q_point));
for (unsigned int d = 0; d < dim; ++d)
- assembling_matrix (i * degree + j,
- dim * q_point + d)
- = shape_value[d];
+ assembling_matrix (index, dim * q_point + d) = shape_value[d];
+
+ ++index;
}
}
// Store the computed
// values.
- for (unsigned int i = 0; i <= degree; ++i)
- for (unsigned int j = 0; j < degree; ++j)
- dof_values[(i + GeometryInfo<dim>::lines_per_face)
- * degree + j + GeometryInfo<dim>::lines_per_face]
- = solution (i * degree + j);
+ {
+ unsigned int index = 0;
+
+ for (unsigned int i = 0; i < fe.dofs_per_face; ++i)
+ if (((dynamic_cast<const FESystem<dim>*> (&fe) != 0)
+ && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).first == base_indices)
+ && (fe.base_element (base_indices.first).face_to_cell_index
+ (GeometryInfo<dim>::lines_per_face * fe.degree, face)
+ <= fe.system_to_base_index (fe.face_to_cell_index (i, face)).second)
+ && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).second
+ < fe.base_element (base_indices.first).face_to_cell_index
+ ((degree + GeometryInfo<dim>::lines_per_face) * fe.degree, face)))
+ || ((dynamic_cast<const FE_Nedelec<dim>*> (&fe) != 0)
+ && (GeometryInfo<dim>::lines_per_face * fe.degree <= i)
+ && (i < (degree + GeometryInfo<dim>::lines_per_face) * fe.degree)))
+ {
+ dof_values[i] = solution (index);
+ dofs_processed[i] = true;
+ ++index;
+ }
+ }
// Now we do the
// same as above
for (unsigned int d = 0; d < dim; ++d)
tmp[d] = values[q_point] (first_vector_component + d);
- for (unsigned int i = 0; i < 2; ++i)
- for (unsigned int j = 0; j <= degree; ++j)
- tmp
- -= dof_values[i * superdegree + j]
- * fe_values[vec].value (cell->get_fe ().face_to_cell_index
- (i * superdegree + j, face), q_point);
+ for (unsigned int i = 0; i < fe.dofs_per_face; ++i)
+ if (((dynamic_cast<const FESystem<dim>*> (&fe) != 0)
+ && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).first == base_indices)
+ && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).second
+ < fe.base_element (base_indices.first).face_to_cell_index (2 * fe.degree, face)))
+ || ((dynamic_cast<const FE_Nedelec<dim>*> (&fe) != 0) && (i < 2 * fe.degree)))
+ tmp -= dof_values[i] * fe_values[vec].value (fe.face_to_cell_index (i, face), q_point);
const double JxW
= std::sqrt (fe_values.JxW (q_point)
for (unsigned int d = 0; d < dim; ++d)
assembling_vector (dim * q_point + d) = JxW * tmp[d];
- for (unsigned int i = 0; i < degree; ++i)
- for (unsigned int j = 0; j <= degree; ++j)
+ unsigned int index = 0;
+
+ for (unsigned int i = 0; i < fe.dofs_per_face; ++i)
+ if (((dynamic_cast<const FESystem<dim>*> (&fe) != 0)
+ && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).first == base_indices)
+ && (fe.base_element (base_indices.first).face_to_cell_index
+ ((degree + GeometryInfo<dim>::lines_per_face) * fe.degree, face)
+ <= fe.system_to_base_index (fe.face_to_cell_index (i, face)).second))
+ || ((dynamic_cast<const FE_Nedelec<dim>*> (&fe) != 0)
+ && ((degree + GeometryInfo<dim>::lines_per_face) * fe.degree <= i)))
{
const Tensor<1, dim> shape_value
- = (JxW
- * fe_values[vec].value (cell->get_fe ().face_to_cell_index
- ((i + degree + GeometryInfo<dim>::lines_per_face)
- * superdegree + j, face), q_point));
+ = JxW * fe_values[vec].value (fe.face_to_cell_index (i, face),
+ q_point);
for (unsigned int d = 0; d < dim; ++d)
- assembling_matrix (i * superdegree + j, dim * q_point + d)
- = shape_value[d];
+ assembling_matrix (index, dim * q_point + d) = shape_value[d];
+
+ ++index;
}
}
assembling_matrix.vmult (cell_rhs, assembling_vector);
cell_matrix_inv.vmult (solution, cell_rhs);
- for (unsigned int i = 0; i < degree; ++i)
- for (unsigned int j = 0; j <= degree; ++j)
- dof_values[(i + degree + GeometryInfo<dim>::lines_per_face) * superdegree + j]
- = solution (i * superdegree + j);
+ unsigned int index = 0;
+
+ for (unsigned int i = 0; i < fe.dofs_per_face; ++i)
+ if (((dynamic_cast<const FESystem<dim>*> (&fe) != 0)
+ && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).first == base_indices)
+ && (fe.base_element (base_indices.first).face_to_cell_index
+ ((degree + GeometryInfo<dim>::lines_per_face) * fe.degree, face)
+ <= fe.system_to_base_index (fe.face_to_cell_index (i, face)).second))
+ || ((dynamic_cast<const FE_Nedelec<dim>*> (&fe) != 0)
+ && ((degree + GeometryInfo<dim>::lines_per_face) * fe.degree <= i)))
+ {
+ dof_values[i] = solution (index);
+ dofs_processed[i] = true;
+ ++index;
+ }
break;
}
update_quadrature_points |
update_values);
+ std::vector<bool> dofs_processed (dofs_per_face);
std::vector<double> dof_values (dofs_per_face);
std::vector<unsigned int> face_dof_indices (dofs_per_face);
typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active ();
// element. If the FE
// is a FESystem, we
// cannot check this.
- if (dynamic_cast<const FESystem<dim>*> (&cell->get_fe ()) == 0) {
+ if (dynamic_cast<const FESystem<dim>*> (&cell->get_fe ()) == 0)
+ {
typedef FiniteElement<dim> FEL;
AssertThrow (dynamic_cast<const FE_Nedelec<dim>*> (&cell->get_fe ()) != 0,
}
for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
+ {
dof_values[dof] = 0.0;
+ dofs_processed[dof] = false;
+ }
// Compute the
// projection of the
::compute_face_projection_curl_conforming (cell, face, fe_face_values,
boundary_function,
first_vector_component,
- dof_values);
+ dof_values, dofs_processed);
cell->face (face)->get_dof_indices (face_dof_indices,
cell->active_fe_index ());
// if the degree of
// freedom is not
// already constrained.
- const double tol = 1e-13;
+ const double tol = 0.5 * cell->get_fe ().degree * 1e-13 / cell->face (face)->diameter ();
for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
- if (constraints.can_store_line (face_dof_indices[dof])
- && !(constraints.is_constrained (face_dof_indices[dof])))
+ if (dofs_processed[dof] && constraints.can_store_line (face_dof_indices[dof])
+ && !(constraints.is_constrained (face_dof_indices[dof])))
{
constraints.add_line (face_dof_indices[dof]);
{
const QGauss<dim - 2> reference_edge_quadrature (2 * superdegree);
const unsigned int degree = superdegree - 1;
- const unsigned int n_dofs = dof_handler.n_dofs ();
hp::QCollection<dim> edge_quadrature_collection;
for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
update_JxW_values |
update_quadrature_points |
update_values);
- std::vector<double> computed_constraints (n_dofs);
- std::vector<int> projected_dofs (n_dofs);
-
- for (unsigned int dof = 0; dof < n_dofs; ++dof) {
- computed_constraints[dof] = 0.0;
- projected_dofs[dof] = -1;
- }
for (; cell != dof_handler.end (); ++cell)
if (cell->at_boundary () && cell->is_locally_owned ())
}
for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
+ {
dof_values[dof] = 0.0;
-
- cell->face (face)->get_dof_indices (face_dof_indices,
- cell->active_fe_index ());
+ dofs_processed[dof] = false;
+ }
// First we compute the
// projection on the
// edges.
for (unsigned int line = 0;
line < GeometryInfo<3>::lines_per_face; ++line)
- {
- // If we have reached
- // this edge through
- // another cell
- // before, we do not
- // do here anything
- // unless we have a
- // good reason, i.e.
- // a higher
- // polynomial degree.
- if (projected_dofs[face_dof_indices[line * superdegree]]
- <
- (int) degree)
- {
- // Compute the
- // projection of
- // the boundary
- // function on the
- // edge.
- internals
- ::compute_edge_projection (cell, face, line,
- fe_edge_values,
- boundary_function,
- first_vector_component,
- dof_values);
- // Mark the
- // projected
- // degrees of
- // freedom.
- for (unsigned int dof = line * superdegree;
- dof < (line + 1) * superdegree; ++dof)
- projected_dofs[face_dof_indices[dof]] = degree;
- }
-
- // If we have
- // computed the
- // values in a
- // previous step of
- // the loop, we just
- // copy the values in
- // the local vector.
- else
- for (unsigned int dof = line * superdegree;
- dof < (line + 1) * superdegree; ++dof)
- dof_values[dof] = computed_constraints[face_dof_indices[dof]];
- }
+ internals
+ ::compute_edge_projection (cell, face, line, fe_edge_values,
+ boundary_function,
+ first_vector_component,
+ dof_values, dofs_processed);
// If there are higher
// order shape
// still some work
// left.
if (degree > 0)
- {
- // Compute the
- // projection of the
- // boundary function
- // on the interior of
- // the face.
- internals
- ::compute_face_projection_curl_conforming (cell, face, fe_face_values,
- boundary_function,
- first_vector_component,
- dof_values);
-
- // Mark the projected
- // degrees of
- // freedom.
- for (unsigned int dof = GeometryInfo<dim>::lines_per_face * superdegree;
- dof < dofs_per_face; ++dof)
- projected_dofs[face_dof_indices[dof]] = degree;
- }
+ internals
+ ::compute_face_projection_curl_conforming (cell, face, fe_face_values,
+ boundary_function,
+ first_vector_component,
+ dof_values,
+ dofs_processed);
// Store the computed
// values in the global
// vector.
+
+ cell->face (face)->get_dof_indices (face_dof_indices,
+ cell->active_fe_index ());
+
const double tol = 0.5 * superdegree * 1e-13 / cell->face (face)->diameter ();
for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
- if (std::abs (computed_constraints[face_dof_indices[dof]] - dof_values[dof]) > tol)
- computed_constraints[face_dof_indices[dof]] = dof_values[dof];
- }
+ if (dofs_processed[dof] && constraints.can_store_line (face_dof_indices[dof])
+ && !(constraints.is_constrained (face_dof_indices[dof])))
+ {
+ constraints.add_line (face_dof_indices[dof]);
- // Add the computed constraints
- // to the constraint matrix, if
- // the degree of freedom is not
- // already constrained.
- for (unsigned int dof = 0; dof < n_dofs; ++dof)
- if ((projected_dofs[dof] != -1) && constraints.can_store_line (dof)
- && !(constraints.is_constrained (dof)))
- {
- constraints.add_line (dof);
- constraints.set_inhomogeneity (dof, computed_constraints[dof]);
- }
+ if (std::abs (dof_values[dof]) > tol)
+ constraints.set_inhomogeneity (face_dof_indices[dof], dof_values[dof]);
+ }
+ }
break;
}
update_JxW_values |
update_quadrature_points |
update_values);
+ std::vector<bool> dofs_processed;
std::vector<double> dof_values;
std::vector<unsigned int> face_dof_indices;
typename hp::DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active ();
}
const unsigned int dofs_per_face = cell->get_fe ().dofs_per_face;
-
+
+ dofs_processed.resize (dofs_per_face);
dof_values.resize (dofs_per_face);
for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
+ {
dof_values[dof] = 0.0;
+ dofs_processed[dof] = false;
+ }
internals
::compute_face_projection_curl_conforming (cell, face, fe_face_values,
boundary_function,
first_vector_component,
- dof_values);
+ dof_values, dofs_processed);
face_dof_indices.resize (dofs_per_face);
cell->face (face)->get_dof_indices (face_dof_indices,
cell->active_fe_index ());
const double tol = 0.5 * cell->get_fe ().degree * 1e-13 / cell->face (face)->diameter ();
for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
- if (constraints.can_store_line (face_dof_indices[dof])
- && !(constraints.is_constrained (face_dof_indices[dof])))
+ if (dofs_processed[dof] && constraints.can_store_line (face_dof_indices[dof])
+ && !(constraints.is_constrained (face_dof_indices[dof])))
{
constraints.add_line (face_dof_indices[dof]);
case 3:
{
- const unsigned int n_dofs = dof_handler.n_dofs ();
hp::QCollection<dim> edge_quadrature_collection;
for (unsigned int i = 0; i < fe_collection.size (); ++i)
update_JxW_values |
update_quadrature_points |
update_values);
- std::vector<double> computed_constraints (n_dofs);
- std::vector<int> projected_dofs (n_dofs);
-
- for (unsigned int dof = 0; dof < n_dofs; ++dof) {
- computed_constraints[dof] = 0.0;
- projected_dofs[dof] = -1;
- }
for (; cell != dof_handler.end (); ++cell)
if (cell->at_boundary () && cell->is_locally_owned ())
const unsigned int superdegree = cell->get_fe ().degree;
const unsigned int degree = superdegree - 1;
const unsigned int dofs_per_face = cell->get_fe ().dofs_per_face;
-
+
+ dofs_processed.resize (dofs_per_face);
dof_values.resize (dofs_per_face);
for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
+ {
dof_values[dof] = 0.0;
-
- face_dof_indices.resize (dofs_per_face);
- cell->face (face)->get_dof_indices (face_dof_indices,
- cell->active_fe_index ());
+ dofs_processed[dof] = false;
+ }
for (unsigned int line = 0;
line < GeometryInfo<dim>::lines_per_face; ++line)
- {
- if (projected_dofs[face_dof_indices[line * superdegree]]
- <
- (int) degree)
- {
- internals
- ::compute_edge_projection (cell, face, line,
- fe_edge_values,
- boundary_function,
- first_vector_component,
- dof_values);
-
- for (unsigned int dof = line * superdegree;
- dof < (line + 1) * superdegree; ++dof)
- projected_dofs[face_dof_indices[dof]] = degree;
- }
-
- else
- for (unsigned int dof = line * superdegree;
- dof < (line + 1) * superdegree; ++dof)
- dof_values[dof] = computed_constraints[face_dof_indices[dof]];
- }
+ internals
+ ::compute_edge_projection (cell, face, line, fe_edge_values,
+ boundary_function,
+ first_vector_component,
+ dof_values, dofs_processed);
+ // If there are higher
+ // order shape
+ // functions, there is
+ // still some work
+ // left.
if (degree > 0)
- {
- internals
- ::compute_face_projection_curl_conforming (cell, face, fe_face_values,
- boundary_function,
- first_vector_component,
- dof_values);
-
- for (unsigned int dof = GeometryInfo<dim>::lines_per_face * superdegree;
- dof < dofs_per_face; ++dof)
- projected_dofs[face_dof_indices[dof]] = degree;
- }
+ internals
+ ::compute_face_projection_curl_conforming (cell, face, fe_face_values,
+ boundary_function,
+ first_vector_component,
+ dof_values, dofs_processed);
+
+ face_dof_indices.resize (dofs_per_face);
+ cell->face (face)->get_dof_indices (face_dof_indices,
+ cell->active_fe_index ());
+
const double tol = 0.5 * superdegree * 1e-13 / cell->face (face)->diameter ();
for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
- if (std::abs (computed_constraints[face_dof_indices[dof]] - dof_values[dof]) > tol)
- computed_constraints[face_dof_indices[dof]] = dof_values[dof];
- }
+ if (dofs_processed[dof] && constraints.can_store_line (face_dof_indices[dof])
+ && !(constraints.is_constrained (face_dof_indices[dof])))
+ {
+ constraints.add_line (face_dof_indices[dof]);
- for (unsigned int dof = 0; dof < n_dofs; ++dof)
- if ((projected_dofs[dof] != -1) && constraints.can_store_line (dof)
- && !(constraints.is_constrained (dof)))
- {
- constraints.add_line (dof);
- constraints.set_inhomogeneity (dof, computed_constraints[dof]);
- }
+ if (std::abs (dof_values[dof]) > tol)
+ constraints.set_inhomogeneity (face_dof_indices[dof], dof_values[dof]);
+ }
+ }
break;
}