Dear all,
I encountered the following issue while trying to solve a Maxwell
eigenvalue problem in dealii v 9.1.1.
When using an FESystem<3> containing one FENedelec and other elements and
trying to set boundary values
via
VectorTools::project_boundary_values_curl_conforming_l2
where the argument of first_vector_component is set to the starting
component of the FENedelec.
The problem has two parts, the first is solved the second is not and input
is appreciated.
Part 1 (solved):
When the FESystem contains two FENedelec of different degree.
Running the attached code (main.cc) ends in debug mode with the
error
The violated condition was:
associated_edge_dofs == degree + 1
Additional information:
Error: Unexpected number of 3D edge DoFs
triggered by line 31 (main.cc).
The error is caused by the fact, that in vector_tools.templates.h the
degree is calculated using
fe.degree - 1
which for FESystem is the wrong value as it is the maximum degree and not
the degree of the selected element.
It seems to me, that a fix for this is to use
fe.base_element(base_indices.first).degree - 1;
a few lines later once the user indicated base_element is known.
Once this is done, in debug mode in the attached example line 31 still
gives an error since
in VectorTools::compute_face_projection_curl_conforming_l2 a 0-by-0 matrix
for the face dofs is to be inverted.
I am not sure why this pops up as it did not appear with the original
version using only FENedelec of lowest order.
Anyways this seems to be easy to fix. A patch for vector_tools.templates.h
(in dealii v9.1.1) is attached
- it is not tested that the computed result is correct, only that this
fixes the wrongly chosen degree and produces no errors.
Part 2 (input needed):
Unfortunately this did not solve the original problem in full. If the
FESystem does not only contain FENedelec, but another FE_Q
and after applying the patch the original error "Unexpected number of 3D
edge DoFs"
still pops up in line 42 of the attached example, this time not because the
wrong degree is selected, but because the
if(((dynamic_cast<const FESystem<dim> *>(&fe) !=
nullptr) ...
statement in line 5440 of the patched vector_tools.templates.h is never
true in this case.
Any input on why this might happen is appreciated - in particular whether
the assumptions made on the ordering of the
dofs is correct for such FESystems.
Thanks
Winni
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/1eea5fd0-6acb-40e1-b55d-43d94d8fe22d%40googlegroups.com.
--- deal9.1.1/vector_tools.templates.h 2019-05-27 04:07:17.000000000 +0200
+++ vector_tools.templates.h 2019-10-29 15:59:30.547095179 +0100
@@ -5345,10 +5345,7 @@
// Initialize the required objects.
const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values();
- // Store degree as fe.degree-1
- // For nedelec elements FE_Nedelec<dim> (0) returns fe.degree = 1.
- const unsigned int degree = fe.degree - 1;
-
+
const std::vector<Point<dim>> &quadrature_points =
fe_values.get_quadrature_points();
std::vector<Vector<double>> values(fe_values.n_quadrature_points,
@@ -5391,6 +5388,11 @@
base_indices.second = (first_vector_component - fe_index_old) /
fe.base_element(i).n_components();
}
+ // Store degree as fe.degree-1
+ // For nedelec elements FE_Nedelec<dim> (0) returns fe.degree = 1.
+ // For FESystem get the degree from the base_element
+ // indicated by the first_vector_component
+ const unsigned int degree = fe.base_element(base_indices.first).degree - 1;
// Find DoFs we want to constrain:
// There are fe.dofs_per_line DoFs associated with the
@@ -5606,7 +5608,6 @@
const FiniteElement<dim> & fe = cell->get_fe();
const std::vector<Point<dim>> &quadrature_points =
fe_face_values.get_quadrature_points();
- const unsigned int degree = fe.degree - 1;
std::vector<Vector<double>> values(fe_face_values.n_quadrature_points,
Vector<double>(fe.n_components()));
@@ -5646,7 +5647,8 @@
base_indices.second = (first_vector_component - fe_index_old) /
fe.base_element(i).n_components();
}
-
+ const unsigned int degree = fe.base_element(base_indices.first).degree - 1;
+
switch (dim)
{
case 2:
@@ -5994,9 +5996,11 @@
}
// Solve lienar system:
- face_matrix_inv.invert(face_matrix);
- face_matrix_inv.vmult(face_solution, face_rhs);
-
+ if(associated_face_dofs > 0)
+ {
+ face_matrix_inv.invert(face_matrix);
+ face_matrix_inv.vmult(face_solution, face_rhs);
+ }
// Store computed DoFs:
for (unsigned int associated_face_dof = 0;
#include <deal.II/grid/grid_generator.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_nedelec.h>
#include <deal.II/fe/mapping_q.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/numerics/vector_tools.h>
using namespace dealii;
int main(int argc, char **argv) {
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
const unsigned int dim=3;
Triangulation<dim> triangulation;
GridGenerator::hyper_cube(triangulation, 0, 1, true);
MappingQ<dim> mapping(1);
ConstraintMatrix constraints;
//Cases with only Nedelec Elements
//Fails
FESystem<dim> fe_system(FE_Nedelec<dim>(1), 1, FE_Nedelec<dim>(0), 1);
//Runs
//FESystem<dim> fe_system(FE_Nedelec<dim>(0), 1, FE_Nedelec<dim>(0), 1);
DoFHandler<dim> dof_handler(triangulation);
dof_handler.distribute_dofs(fe_system);
constraints.clear();
VectorTools::project_boundary_values_curl_conforming_l2(dof_handler, dim,
Functions::ZeroFunction<dim>(dim + dim), 0, constraints, mapping);
constraints.close();
//Other potentially useful cases of FESystem including Nedelec Elements
//Fails
FESystem<dim> fe_system2(FE_Q<dim>(2), 1,FE_Nedelec<dim>(0), 1);
DoFHandler<dim> dof_handler2(triangulation);
dof_handler2.distribute_dofs(fe_system2);
constraints.clear();
VectorTools::project_boundary_values_curl_conforming_l2(dof_handler2, 1,
Functions::ZeroFunction<dim>(dim + 1), 0, constraints, mapping);
constraints.close();
return 0;
}
