--- /dev/null
+New: FE_SimplexP and FE_SimplexDGP now support cubic elements.
+<br>
+(David Wells and Peter Munch, 2023/10/11)
namespace
{
+ // TODO: replace this with QProjector once QProjector learns how to project
+ // quadrature points onto separate faces of simplices
+ template <int dim>
+ Point<dim>
+ face_midpoint(const unsigned int face_no)
+ {
+ const auto reference_cell = ReferenceCells::get_simplex<dim>();
+ const auto face_reference_cell =
+ reference_cell.face_reference_cell(face_no);
+
+ Point<dim> midpoint;
+ for (const auto face_vertex_no : face_reference_cell.vertex_indices())
+ {
+ const auto vertex_no = reference_cell.face_to_cell_vertices(
+ face_no,
+ face_vertex_no,
+ ReferenceCell::default_combined_face_orientation());
+
+ midpoint += reference_cell.template vertex<dim>(vertex_no);
+ }
+
+ midpoint /= reference_cell.face_reference_cell(0).n_vertices();
+ return midpoint;
+ }
+
/**
* Helper function to set up the dpo vector of FE_SimplexP for a given @p dim and
* @p degree.
return {1, 0};
case 2:
return {1, 1};
+ case 3:
+ return {1, 2};
default:
Assert(false, ExcNotImplemented());
}
return {1, 0, 0};
case 2:
return {1, 1, 0};
+ case 3:
+ return {1, 2, 1};
default:
Assert(false, ExcNotImplemented());
}
return {1, 0, 0, 0};
case 2:
return {1, 1, 0, 0};
+ case 3:
+ return {1, 2, 1, 0};
default:
Assert(false, ExcNotImplemented());
}
std::vector<Point<dim>>
unit_support_points_fe_p(const unsigned int degree)
{
+ Assert(dim != 0, ExcInternalError());
+ Assert(degree <= 3, ExcNotImplemented());
std::vector<Point<dim>> unit_points;
- // If we do dim - 1 we can get here in dim = 0:
- if (dim == 0)
- {
- unit_points.emplace_back();
- return unit_points;
- }
+ const auto reference_cell = ReferenceCells::get_simplex<dim>();
- const auto reference_cell = ReferenceCells::get_simplex<dim>();
// Piecewise constants are a special case: use a support point at the
// centroid and only the centroid
if (degree == 0)
return unit_points;
}
- Assert(degree <= 2, ExcNotImplemented());
- for (const auto &vertex_no : reference_cell.vertex_indices())
- unit_points.emplace_back(reference_cell.template vertex<dim>(vertex_no));
-
- if (degree == 2)
- for (const auto &line_no : reference_cell.line_indices())
- {
- const auto v0 = reference_cell.template vertex<dim>(
- reference_cell.line_to_cell_vertices(line_no, 0));
- const auto v1 = reference_cell.template vertex<dim>(
- reference_cell.line_to_cell_vertices(line_no, 1));
- unit_points.emplace_back((v0 + v1) / 2.0);
- }
+ // otherwise write everything as linear combinations of vertices
+ const auto dpo = get_dpo_vector_fe_p(dim, degree);
+ Assert(dpo.size() == dim + 1, ExcInternalError());
+ Assert(dpo[0] == 1, ExcNotImplemented());
+ // vertices:
+ for (const unsigned int d : reference_cell.vertex_indices())
+ unit_points.push_back(reference_cell.template vertex<dim>(d));
+ // lines:
+ for (const unsigned int l : reference_cell.line_indices())
+ {
+ const Point<dim> p0 =
+ unit_points[reference_cell.line_to_cell_vertices(l, 0)];
+ const Point<dim> p1 =
+ unit_points[reference_cell.line_to_cell_vertices(l, 1)];
+ for (unsigned int p = 0; p < dpo[1]; ++p)
+ unit_points.push_back((double(dpo[1] - p) / (dpo[1] + 1)) * p0 +
+ (double(p + 1) / (dpo[1] + 1)) * p1);
+ }
+ // quads:
+ if (dim > 1 && dpo[2] > 0)
+ {
+ Assert(dpo[2] == 1, ExcNotImplemented());
+ if (dim == 2)
+ unit_points.push_back(reference_cell.template barycenter<dim>());
+ if (dim == 3)
+ for (const unsigned int f : reference_cell.face_indices())
+ unit_points.push_back(face_midpoint<dim>(f));
+ }
return unit_points;
}
+ template <>
+ std::vector<Point<0>>
+ unit_support_points_fe_p(const unsigned int /*degree*/)
+ {
+ return {Point<0>()};
+ }
+
/**
* Set up a vector that contains the unit (reference) cell's faces support
* points for FE_SimplexP and sufficiently similar elements.
FullMatrix<double>
constraints_fe_p<2>(const unsigned int degree)
{
- const unsigned int dim = 2;
-
- Assert(degree <= 2, ExcNotImplemented());
+ constexpr int dim = 2;
+ Assert(degree <= 3, ExcNotImplemented());
// the following implements the 2d case
// (the 3d case is not implemented yet)
unit_support_points_fe_p<dim>(degree),
unit_face_support_points_fe_p<dim>(degree, FiniteElementData<dim>::H1),
constraints_fe_p<dim>(degree))
-{}
+{
+ if (degree > 2)
+ for (unsigned int i = 0; i < this->n_dofs_per_line(); ++i)
+ this->adjust_line_dof_index_for_line_orientation_table[i] =
+ this->n_dofs_per_line() - 1 - i - i;
+ // We do not support multiple DoFs per quad yet
+ Assert(degree <= 3, ExcNotImplemented());
+}
{
#ifdef SIMPLEX
FE_SimplexP<dim> fe(degree);
- QGaussSimplex<dim> quadrature(degree + 2);
+ QGaussSimplex<dim> quadrature(degree + 1);
#else
FE_Q<dim> fe(degree);
- QGauss<dim> quadrature(degree + 2);
+ QGauss<dim> quadrature(degree + 1);
#endif
deallog << "FE = " << fe.get_name() << std::endl;
+ deallog << std::setprecision(4);
double previous_error = 1.0;
reference_cell.template get_default_linear_mapping<dim>();
dummy.close();
- const bool l2_projection = false;
+ SparsityPattern sparsity_pattern;
+ {
+ DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
+ DoFTools::make_sparsity_pattern(dof_handler, dsp);
+ sparsity_pattern.copy_from(dsp);
+ }
- if (l2_projection)
- {
- VectorTools::project(
- mapping, dof_handler, dummy, quadrature, function, solution);
- }
- else
- {
- SparsityPattern sparsity_pattern;
- {
- DynamicSparsityPattern dsp(dof_handler.n_dofs(),
- dof_handler.n_dofs());
- DoFTools::make_sparsity_pattern(dof_handler, dsp);
- sparsity_pattern.copy_from(dsp);
- }
-
- SparseMatrix<double> h1_matrix(sparsity_pattern);
- SparseMatrix<double> laplace_matrix(sparsity_pattern);
-
- MatrixCreator::create_mass_matrix(mapping,
- dof_handler,
- quadrature,
- h1_matrix);
- MatrixCreator::create_laplace_matrix(mapping,
- dof_handler,
- quadrature,
- laplace_matrix);
-
- h1_matrix.add(0.0, laplace_matrix);
-
- Vector<double> rhs(solution.size());
- VectorTools::create_right_hand_side(
- mapping, dof_handler, quadrature, ForcingH1<dim>(), rhs);
-
- SolverControl solver_control(1000,
- 1e-12 * rhs.l2_norm(),
- false,
- false);
- SolverCG<Vector<double>> cg(solver_control);
-
- cg.solve(h1_matrix, solution, rhs, PreconditionIdentity());
- }
+ SparseMatrix<double> h1_matrix(sparsity_pattern);
+ SparseMatrix<double> laplace_matrix(sparsity_pattern);
+
+ MatrixCreator::create_mass_matrix(mapping,
+ dof_handler,
+ quadrature,
+ h1_matrix);
+ MatrixCreator::create_laplace_matrix(mapping,
+ dof_handler,
+ quadrature,
+ laplace_matrix);
+
+ h1_matrix.add(0.0, laplace_matrix);
+
+ Vector<double> rhs(solution.size());
+ VectorTools::create_right_hand_side(
+ mapping, dof_handler, quadrature, ForcingH1<dim>(), rhs);
+
+ SolverControl solver_control(1000, 1e-12 * rhs.l2_norm(), false, false);
+ SolverCG<Vector<double>> cg(solver_control);
+
+ cg.solve(h1_matrix, solution, rhs, PreconditionIdentity());
VectorTools::integrate_difference(mapping,
dof_handler,
solution,
function,
cell_errors,
- Quadrature<dim>(
- fe.get_unit_support_points()),
- VectorTools::Linfty_norm);
-
- const double max_error =
- *std::max_element(cell_errors.begin(), cell_errors.end());
- deallog << "max error = " << max_error << std::endl;
- if (max_error != 0.0)
- deallog << "ratio = " << previous_error / max_error << std::endl;
- previous_error = max_error;
+ quadrature,
+ VectorTools::L2_norm);
+ const double L2_error =
+ VectorTools::compute_global_error(tria,
+ cell_errors,
+ VectorTools::L2_norm);
+
+ deallog << "L2 error = " << L2_error << std::endl;
+ if (L2_error != 0.0)
+ deallog << "ratio = " << previous_error / L2_error << std::endl;
+ previous_error = L2_error;
#if 0
if (dim == 2)
test<2>(1);
test<2>(2);
+ test<2>(3);
test<3>(1);
test<3>(2);
+ test<3>(3);
}
DEAL::FE = FE_SimplexP<2>(1)
-DEAL::max error = 0.332362
-DEAL::ratio = 3.00877
-DEAL::max error = 0.100426
-DEAL::ratio = 3.30953
-DEAL::max error = 0.0271946
-DEAL::ratio = 3.69285
-DEAL::max error = 0.00747939
-DEAL::ratio = 3.63594
+DEAL::L2 error = 0.05224
+DEAL::ratio = 19.14
+DEAL::L2 error = 0.02949
+DEAL::ratio = 1.772
+DEAL::L2 error = 0.007768
+DEAL::ratio = 3.796
+DEAL::L2 error = 0.001951
+DEAL::ratio = 3.981
DEAL::FE = FE_SimplexP<2>(2)
-DEAL::max error = 0.0430400
-DEAL::ratio = 23.2342
-DEAL::max error = 0.0116055
-DEAL::ratio = 3.70859
-DEAL::max error = 0.00232162
-DEAL::ratio = 4.99887
-DEAL::max error = 0.000336599
-DEAL::ratio = 6.89729
+DEAL::L2 error = 0.004078
+DEAL::ratio = 245.2
+DEAL::L2 error = 0.002064
+DEAL::ratio = 1.975
+DEAL::L2 error = 0.0003512
+DEAL::ratio = 5.879
+DEAL::L2 error = 5.040e-05
+DEAL::ratio = 6.967
+DEAL::FE = FE_SimplexP<2>(3)
+DEAL::L2 error = 0.004061
+DEAL::ratio = 246.3
+DEAL::L2 error = 0.0002464
+DEAL::ratio = 16.48
+DEAL::L2 error = 1.603e-05
+DEAL::ratio = 15.38
+DEAL::L2 error = 1.009e-06
+DEAL::ratio = 15.88
DEAL::FE = FE_SimplexP<3>(1)
-DEAL::max error = 0.620494
-DEAL::ratio = 1.61162
-DEAL::max error = 0.259138
-DEAL::ratio = 2.39445
-DEAL::max error = 0.0619310
-DEAL::ratio = 4.18431
-DEAL::max error = 0.0151022
-DEAL::ratio = 4.10080
+DEAL::L2 error = 0.09826
+DEAL::ratio = 10.18
+DEAL::L2 error = 0.02938
+DEAL::ratio = 3.344
+DEAL::L2 error = 0.007387
+DEAL::ratio = 3.977
+DEAL::L2 error = 0.001763
+DEAL::ratio = 4.191
DEAL::FE = FE_SimplexP<3>(2)
-DEAL::max error = 0.0550999
-DEAL::ratio = 18.1489
-DEAL::max error = 0.0466852
-DEAL::ratio = 1.18024
-DEAL::max error = 0.00502563
-DEAL::ratio = 9.28942
-DEAL::max error = 0.000810078
-DEAL::ratio = 6.20388
+DEAL::L2 error = 0.01087
+DEAL::ratio = 91.95
+DEAL::L2 error = 0.004092
+DEAL::ratio = 2.657
+DEAL::L2 error = 0.0007005
+DEAL::ratio = 5.842
+DEAL::L2 error = 0.0001018
+DEAL::ratio = 6.880
+DEAL::FE = FE_SimplexP<3>(3)
+DEAL::L2 error = 0.006416
+DEAL::ratio = 155.9
+DEAL::L2 error = 0.0003557
+DEAL::ratio = 18.03
+DEAL::L2 error = 2.324e-05
+DEAL::ratio = 15.31
+DEAL::L2 error = 1.447e-06
+DEAL::ratio = 16.06
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_values.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria.h>
#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_control.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparsity_pattern.h>
#include <deal.II/lac/vector.h>
#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/matrix_creator.h>
#include <deal.II/numerics/vector_tools_integrate_difference.h>
#include <deal.II/numerics/vector_tools_project.h>
+#include <deal.II/numerics/vector_tools_rhs.h>
#include "../tests.h"
{
FE_SimplexDGP<dim> fe(degree);
deallog << "FE = " << fe.get_name() << std::endl;
+ deallog << std::setprecision(4);
QGaussSimplex<dim> quadrature(degree + 1);
double previous_error = 1.0;
const auto &mapping =
reference_cell.template get_default_linear_mapping<dim>();
dummy.close();
- VectorTools::project(
- mapping, dof_handler, dummy, quadrature, function, solution);
+
+ SparsityPattern sparsity_pattern;
+ {
+ DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
+ DoFTools::make_sparsity_pattern(dof_handler, dsp);
+ sparsity_pattern.copy_from(dsp);
+ }
+
+ SparseMatrix<double> mass_matrix(sparsity_pattern);
+ MatrixCreator::create_mass_matrix(mapping,
+ dof_handler,
+ quadrature,
+ mass_matrix);
+ Vector<double> rhs(solution.size());
+ VectorTools::create_right_hand_side(
+ mapping, dof_handler, quadrature, function, rhs);
+
+ SolverControl solver_control(1000, 1e-12 * rhs.l2_norm(), false, false);
+ SolverCG<Vector<double>> cg(solver_control);
+
+ cg.solve(mass_matrix, solution, rhs, PreconditionIdentity());
VectorTools::integrate_difference(mapping,
dof_handler,
solution,
function,
cell_errors,
- Quadrature<dim>(
- fe.get_unit_support_points()),
- VectorTools::Linfty_norm);
-
- const double max_error =
- *std::max_element(cell_errors.begin(), cell_errors.end());
- deallog << "max error = " << max_error << std::endl;
- if (max_error != 0.0)
- deallog << "ratio = " << previous_error / max_error << std::endl;
- previous_error = max_error;
+ quadrature,
+ VectorTools::L2_norm);
+ const double L2_error =
+ VectorTools::compute_global_error(tria,
+ cell_errors,
+ VectorTools::L2_norm);
+
+ deallog << "L2 error = " << L2_error << std::endl;
+ if (L2_error != 0.0)
+ deallog << "ratio = " << previous_error / L2_error << std::endl;
+ previous_error = L2_error;
#if 0
if (dim == 2)
test<2>(1);
test<2>(2);
+ test<2>(3);
test<3>(1);
test<3>(2);
+ test<3>(3);
}
DEAL::FE = FE_SimplexDGP<2>(1)
-DEAL::max error = 0.113889
-DEAL::ratio = 8.78048
-DEAL::max error = 0.0302392
-DEAL::ratio = 3.76627
-DEAL::max error = 0.00767275
-DEAL::ratio = 3.94112
-DEAL::max error = 0.00192529
-DEAL::ratio = 3.98525
+DEAL::L2 error = 0.009794
+DEAL::ratio = 102.1
+DEAL::L2 error = 0.002479
+DEAL::ratio = 3.952
+DEAL::L2 error = 0.0006215
+DEAL::ratio = 3.988
+DEAL::L2 error = 0.0001555
+DEAL::ratio = 3.997
DEAL::FE = FE_SimplexDGP<2>(2)
-DEAL::max error = 0.0174164
-DEAL::ratio = 57.4171
-DEAL::max error = 0.00227404
-DEAL::ratio = 7.65880
-DEAL::max error = 0.000287343
-DEAL::ratio = 7.91404
-DEAL::max error = 3.60148e-05
-DEAL::ratio = 7.97847
+DEAL::L2 error = 0.0002017
+DEAL::ratio = 4957.
+DEAL::L2 error = 2.334e-05
+DEAL::ratio = 8.644
+DEAL::L2 error = 2.854e-06
+DEAL::ratio = 8.178
+DEAL::L2 error = 3.547e-07
+DEAL::ratio = 8.046
+DEAL::FE = FE_SimplexDGP<2>(3)
+DEAL::L2 error = 0.0001932
+DEAL::ratio = 5175.
+DEAL::L2 error = 1.224e-05
+DEAL::ratio = 15.79
+DEAL::L2 error = 7.677e-07
+DEAL::ratio = 15.95
+DEAL::L2 error = 4.802e-08
+DEAL::ratio = 15.99
DEAL::FE = FE_SimplexDGP<3>(1)
-DEAL::max error = 0.208217
-DEAL::ratio = 4.80269
-DEAL::max error = 0.0629409
-DEAL::ratio = 3.30813
-DEAL::max error = 0.0164605
-DEAL::ratio = 3.82376
-DEAL::max error = 0.00416117
-DEAL::ratio = 3.95573
+DEAL::L2 error = 0.01241
+DEAL::ratio = 80.59
+DEAL::L2 error = 0.003150
+DEAL::ratio = 3.939
+DEAL::L2 error = 0.0007879
+DEAL::ratio = 3.998
+DEAL::L2 error = 0.0001968
+DEAL::ratio = 4.004
DEAL::FE = FE_SimplexDGP<3>(2)
-DEAL::max error = 0.0415553
-DEAL::ratio = 24.0643
-DEAL::max error = 0.00641441
-DEAL::ratio = 6.47844
-DEAL::max error = 0.000842371
-DEAL::ratio = 7.61471
-DEAL::max error = 0.000106584
-DEAL::ratio = 7.90336
+DEAL::L2 error = 0.001633
+DEAL::ratio = 612.3
+DEAL::L2 error = 0.0002072
+DEAL::ratio = 7.883
+DEAL::L2 error = 2.599e-05
+DEAL::ratio = 7.971
+DEAL::L2 error = 3.252e-06
+DEAL::ratio = 7.993
+DEAL::FE = FE_SimplexDGP<3>(3)
+DEAL::L2 error = 0.0002064
+DEAL::ratio = 4845.
+DEAL::L2 error = 1.309e-05
+DEAL::ratio = 15.77
+DEAL::L2 error = 8.214e-07
+DEAL::ratio = 15.94
+DEAL::L2 error = 5.138e-08
+DEAL::ratio = 15.99
projects / dealii.git / commitdiff