* @return A FullMatrix<double> representing the assembled mass matrix.
*
*/
- FullMatrix<double>
+ template <typename Number = double>
+ FullMatrix<Number>
create_1d_cell_mass_matrix(
const FiniteElement<1> &fe,
- const double &h,
+ const Number &h,
const std::pair<bool, bool> include_endpoints = {true, true},
std::vector<unsigned int> numbering = std::vector<unsigned int>());
* numbering of the degrees of freedom. If empty, the
* numbering of the finite element is used.
*/
- FullMatrix<double>
+ template <typename Number = double>
+ FullMatrix<Number>
create_1d_cell_laplace_matrix(
const FiniteElement<1> &fe,
- const double &h,
+ const Number &h,
const std::pair<bool, bool> include_endpoints = {true, true},
std::vector<unsigned int> numbering = std::vector<unsigned int>());
* if the number of cells is above 10 an exception will be thrown.
*
*/
- FullMatrix<double>
+ template <typename Number = double>
+ FullMatrix<Number>
create_1D_discretization_matrix(
- FullMatrix<double> &cell_matrix,
+ FullMatrix<Number> &cell_matrix,
const unsigned int &n_cells,
const unsigned int &overlap,
const std::pair<bool, bool> include_endpoints = {true, true});
*
* @return A full matrix representing the 1D ghost penalty matrix.
*/
- FullMatrix<double>
+ template <typename Number = double>
+ FullMatrix<Number>
create_1d_ghost_penalty_matrix(
const FiniteElement<1> &fe,
- const double h,
- std::vector<double> coefficients = std::vector<double>());
+ const Number h,
+ std::vector<Number> coefficients = std::vector<Number>());
*
* @return A full matrix representing the 1D ghost penalty.
*/
- FullMatrix<double>
+ template <typename Number = double>
+ FullMatrix<Number>
create_1d_ghost_penalty_matrix(
const std::vector<Polynomials::Polynomial<double>>
&polynomial_basis_derivative,
fe, quadrature, boundary_ids, cell_extent, n_overlap);
}
- FullMatrix<double>
+ template <typename Number>
+ FullMatrix<Number>
+ create_1d_cell_mass_matrix(const FiniteElement<1> &fe,
+ const Number &h,
+ const std::pair<bool, bool> include_endpoints,
+ std::vector<unsigned int> numbering)
+ {
+ if (dynamic_cast<const FE_DGQ<1> *>(&fe) == nullptr &&
+ numbering.size() == 0)
+ {
+ Assert(
+ include_endpoints.first == true && include_endpoints.second == true,
+ ExcMessage(
+ "You tried to generate a 1D mass matrix with excluding boundary "
+ "dofs for a non-DGQ element without providing a numbering."));
+ }
+
+ if (numbering.size() == 0)
+ {
+ numbering.resize(fe.dofs_per_cell);
+ std::iota(numbering.begin(), numbering.end(), 0);
+ }
+
+ const unsigned int degree = fe.degree;
+ const unsigned int n_dofs_per_cell = fe.dofs_per_cell;
+ const Number &JxW = h;
+ QGauss<1> quadrature(degree + 1);
+
+ FullMatrix<Number> cell_matrix(n_dofs_per_cell, n_dofs_per_cell);
+ cell_matrix = 0;
+
+ unsigned int start_dof = include_endpoints.first ? 0 : 1;
+ unsigned int end_dof =
+ include_endpoints.second ? n_dofs_per_cell : n_dofs_per_cell - 1;
+ const unsigned int shift = include_endpoints.first ? 0 : 1;
+
+ for (unsigned int i = start_dof; i < end_dof; ++i)
+ for (unsigned int j = start_dof; j < end_dof; ++j)
+ for (unsigned int q = 0; q < quadrature.size(); ++q)
+ cell_matrix(i - shift, j - shift) +=
+ (fe.shape_value(numbering[i], quadrature.point(q)) *
+ fe.shape_value(numbering[j], quadrature.point(q))) *
+ JxW * quadrature.weight(q);
+
+ return cell_matrix;
+ }
+
+ template <typename Number>
+ FullMatrix<Number>
create_1d_cell_laplace_matrix(const FiniteElement<1> &fe,
- const double &h,
+ const Number &h,
const std::pair<bool, bool> include_endpoints,
std::vector<unsigned int> numbering)
{
const unsigned int degree = fe.degree;
const unsigned int n_dofs_per_cell = fe.dofs_per_cell;
- const double &JxW = h;
+ const Number &JxW = h;
QGauss<1> quadrature(degree + 1);
- FullMatrix<double> cell_matrix(n_dofs_per_cell, n_dofs_per_cell);
+ FullMatrix<Number> cell_matrix(n_dofs_per_cell, n_dofs_per_cell);
cell_matrix = 0;
unsigned int start_dof = include_endpoints.first ? 0 : 1;
return cell_matrix;
}
- FullMatrix<double>
- create_1D_discretization_matrix(FullMatrix<double> &cell_matrix,
+ template <typename Number>
+ FullMatrix<Number>
+ create_1D_discretization_matrix(FullMatrix<Number> &cell_matrix,
const unsigned int &n_cells,
const unsigned int &overlap,
const std::pair<bool, bool> include_endpoints)
if (!include_endpoints.second)
n_total_dofs -= 1;
- FullMatrix<double> result_matrix(n_total_dofs, n_total_dofs);
+ FullMatrix<Number> result_matrix(n_total_dofs, n_total_dofs);
result_matrix = 0;
const unsigned int left_shift = include_endpoints.first ? 0 : 1;
- FullMatrix<double>
+ template <typename Number>
+ FullMatrix<Number>
create_1d_ghost_penalty_matrix(const FiniteElement<1> &fe,
- const double h,
- std::vector<double> coefficients)
+ const Number h,
+ std::vector<Number> coefficients)
{
Assert(dynamic_cast<const FE_Q<1> *>(&fe) != nullptr, ExcNotImplemented());
Assert(h > 0, ExcMessage("Provided element size h is negative"));
}
- FullMatrix<double> penalty_matrix =
+ FullMatrix<Number> penalty_matrix =
create_1d_ghost_penalty_matrix(polynomial_basis[1]);
penalty_matrix *= coefficients[0];
for (unsigned int k = 2; k < degree + 1; ++k)
{
- FullMatrix<double> kth_matrix =
+ FullMatrix<Number> kth_matrix =
create_1d_ghost_penalty_matrix(polynomial_basis[k]);
penalty_matrix.add(coefficients[k - 1], kth_matrix);
}
}
-
- FullMatrix<double>
+ template <typename Number>
+ FullMatrix<Number>
create_1d_ghost_penalty_matrix(
const std::vector<Polynomials::Polynomial<double>>
&polynomial_basis_derivative,
const unsigned int n_total_dofs = 2 * n_dofs_per_cell - overlap;
const unsigned int shift = n_dofs_per_cell - overlap;
- FullMatrix<double> penalty_matrix(n_total_dofs, n_total_dofs);
+ FullMatrix<Number> penalty_matrix(n_total_dofs, n_total_dofs);
std::vector<double> values_left(n_dofs_per_cell);
std::vector<double> values_right(n_dofs_per_cell);
void
-print_matrix(FullMatrix<double> &matrix)
+print_matrix(const FullMatrix<double> &matrix, const std::string &label)
{
- for (unsigned int i = 0; i < matrix.m(); ++i)
- {
- for (unsigned int j = 0; j < matrix.n(); ++j)
- deallog << std::setw(10) << std::fixed << std::setprecision(4)
- << matrix(i, j) << " ";
-
- deallog << std::endl;
- }
- deallog << std::endl;
+ std::stringstream sstring;
+ sstring << label << "\n";
+ matrix.print_formatted(sstring, 4, true, 10, "0.");
+ deallog << sstring.str() << std::endl;
}
template <int fe_degree>
void
test()
{
- FE_Q<1> fe_q(fe_degree);
- FE_Q<1> fe_dgq(fe_degree);
- auto mass_matrix =
+ FE_Q<1> fe_q(fe_degree);
+ FE_DGQ<1> fe_dgq(fe_degree);
+ auto mass_matrix =
TensorProductMatrixCreator::create_1d_cell_mass_matrix(fe_q,
- 1,
+ 1.,
{true, true});
- mass_matrix.print_formatted(std::cout, 4, true, 10, "0.");
+ print_matrix(mass_matrix,
+ "Mass matrix for FE_Q<1>(" + std::to_string(fe_degree) + "):");
- std::cout << std::endl;
FEEvaluation<1, fe_degree, fe_degree + 1, 1, double> fe_eval(
auto mass_matrix_renumbered =
TensorProductMatrixCreator::create_1d_cell_mass_matrix(
fe_q, 1., {true, true}, fe_eval.get_internal_dof_numbering());
- mass_matrix_renumbered.print_formatted(std::cout, 3, true, 10, "0.");
- std::cout << std::endl;
+ print_matrix(mass_matrix_renumbered,
+ "Mass matrix for FE_Q<1>(" + std::to_string(fe_degree) +
+ ") renumbered with FEEvaluation:");
auto mass_matrix_dgq =
TensorProductMatrixCreator::create_1d_cell_mass_matrix(fe_dgq,
1.,
{true, true});
+ print_matrix(mass_matrix_dgq,
+ "Mass matrix for FE_DGQ<1>(" + std::to_string(fe_degree) + "):");
+
for (unsigned int i = 0; i < mass_matrix_dgq.m(); ++i)
for (unsigned int j = 0; j < mass_matrix_dgq.n(); ++j)
if (std::abs(mass_matrix_dgq(i, j) - mass_matrix_renumbered(i, j)) >
1e-12)
- std::cout << "Different at " << i << " " << j << std::endl;
+ deallog << "Different at " << i << " " << j << std::endl;
- std::cout << std::endl;
auto mass_matrix_patch =
TensorProductMatrixCreator::create_1D_discretization_matrix(
- mass_matrix_renumbered, 4, 1, {true, true});
+ mass_matrix_renumbered, 2, 1, {true, true});
- mass_matrix_patch.print_formatted(std::cout, 4, true, 10, "0.");
- std::cout << std::endl;
+ print_matrix(mass_matrix_patch, "Mass matrix patch:");
auto laplace_matrix =
TensorProductMatrixCreator::create_1d_cell_laplace_matrix(fe_q,
1.,
{true, true});
- laplace_matrix.print_formatted(std::cout, 4, true, 10, "0.");
- std::cout << std::endl;
+ print_matrix(laplace_matrix,
+ "Laplace matrix for FE_Q<1>(" + std::to_string(fe_degree) +
+ "):");
auto laplace_matrix_renumbered =
TensorProductMatrixCreator::create_1d_cell_laplace_matrix(
fe_q, 1., {true, true}, fe_eval.get_internal_dof_numbering());
- laplace_matrix_renumbered.print_formatted(std::cout, 4, true, 10, "0.");
- std::cout << std::endl;
+ print_matrix(laplace_matrix_renumbered,
+ "Laplace matrix for FE_Q<1>(" + std::to_string(fe_degree) +
+ ") renumbered with FEEvaluation:");
}
deallog << std::fixed;
test<1>();
+ test<2>();
test<3>();
}
\ No newline at end of file
-DEAL::1.0000 -2.0000 1.0000
-DEAL::-2.0000 4.0000 -2.0000
-DEAL::1.0000 -2.0000 1.0000
-DEAL::
-DEAL::51.0000 -130.3444 160.3444 -62.0000 -32.3607 12.3607 1.0000
-DEAL::-130.3444 338.1966 -425.0000 148.8854 125.0000 -69.0983 12.3607
-DEAL::160.3444 -425.0000 561.8034 -208.8854 -180.9017 125.0000 -32.3607
-DEAL::-62.0000 148.8854 -208.8854 244.0000 -208.8854 148.8854 -62.0000
-DEAL::-32.3607 125.0000 -180.9017 -208.8854 561.8034 -425.0000 160.3444
-DEAL::12.3607 -69.0983 125.0000 148.8854 -425.0000 338.1966 -130.3444
-DEAL::1.0000 12.3607 -32.3607 -62.0000 160.3444 -130.3444 51.0000
-DEAL::
+DEAL::Mass matrix for FE_Q<1>(1):
+3.3333e-01 1.6667e-01
+1.6667e-01 3.3333e-01
+
+DEAL::Mass matrix for FE_Q<1>(1) renumbered with FEEvaluation:
+3.3333e-01 1.6667e-01
+1.6667e-01 3.3333e-01
+
+DEAL::Mass matrix for FE_DGQ<1>(1):
+3.3333e-01 1.6667e-01
+1.6667e-01 3.3333e-01
+
+DEAL::Mass matrix patch:
+3.3333e-01 1.6667e-01 0.
+1.6667e-01 6.6667e-01 1.6667e-01
+ 0. 1.6667e-01 3.3333e-01
+
+DEAL::Laplace matrix for FE_Q<1>(1):
+1.0000e+00 -1.0000e+00
+-1.0000e+00 1.0000e+00
+
+DEAL::Laplace matrix for FE_Q<1>(1) renumbered with FEEvaluation:
+1.0000e+00 -1.0000e+00
+-1.0000e+00 1.0000e+00
+
+DEAL::Mass matrix for FE_Q<1>(2):
+1.3333e-01 -3.3333e-02 6.6667e-02
+-3.3333e-02 1.3333e-01 6.6667e-02
+6.6667e-02 6.6667e-02 5.3333e-01
+
+DEAL::Mass matrix for FE_Q<1>(2) renumbered with FEEvaluation:
+1.3333e-01 6.6667e-02 -3.3333e-02
+6.6667e-02 5.3333e-01 6.6667e-02
+-3.3333e-02 6.6667e-02 1.3333e-01
+
+DEAL::Mass matrix for FE_DGQ<1>(2):
+1.3333e-01 6.6667e-02 -3.3333e-02
+6.6667e-02 5.3333e-01 6.6667e-02
+-3.3333e-02 6.6667e-02 1.3333e-01
+
+DEAL::Mass matrix patch:
+1.3333e-01 6.6667e-02 -3.3333e-02 0. 0.
+6.6667e-02 5.3333e-01 6.6667e-02 0. 0.
+-3.3333e-02 6.6667e-02 2.6667e-01 6.6667e-02 -3.3333e-02
+ 0. 0. 6.6667e-02 5.3333e-01 6.6667e-02
+ 0. 0. -3.3333e-02 6.6667e-02 1.3333e-01
+
+DEAL::Laplace matrix for FE_Q<1>(2):
+2.3333e+00 3.3333e-01 -2.6667e+00
+3.3333e-01 2.3333e+00 -2.6667e+00
+-2.6667e+00 -2.6667e+00 5.3333e+00
+
+DEAL::Laplace matrix for FE_Q<1>(2) renumbered with FEEvaluation:
+2.3333e+00 -2.6667e+00 3.3333e-01
+-2.6667e+00 5.3333e+00 -2.6667e+00
+3.3333e-01 -2.6667e+00 2.3333e+00
+
+DEAL::Mass matrix for FE_Q<1>(3):
+7.1429e-02 1.1905e-02 2.6620e-02 -2.6620e-02
+1.1905e-02 7.1429e-02 -2.6620e-02 2.6620e-02
+2.6620e-02 -2.6620e-02 3.5714e-01 5.9524e-02
+-2.6620e-02 2.6620e-02 5.9524e-02 3.5714e-01
+
+DEAL::Mass matrix for FE_Q<1>(3) renumbered with FEEvaluation:
+7.1429e-02 2.6620e-02 -2.6620e-02 1.1905e-02
+2.6620e-02 3.5714e-01 5.9524e-02 -2.6620e-02
+-2.6620e-02 5.9524e-02 3.5714e-01 2.6620e-02
+1.1905e-02 -2.6620e-02 2.6620e-02 7.1429e-02
+
+DEAL::Mass matrix for FE_DGQ<1>(3):
+7.1429e-02 2.6620e-02 -2.6620e-02 1.1905e-02
+2.6620e-02 3.5714e-01 5.9524e-02 -2.6620e-02
+-2.6620e-02 5.9524e-02 3.5714e-01 2.6620e-02
+1.1905e-02 -2.6620e-02 2.6620e-02 7.1429e-02
+
+DEAL::Mass matrix patch:
+7.1429e-02 2.6620e-02 -2.6620e-02 1.1905e-02 0. 0. 0.
+2.6620e-02 3.5714e-01 5.9524e-02 -2.6620e-02 0. 0. 0.
+-2.6620e-02 5.9524e-02 3.5714e-01 2.6620e-02 0. 0. 0.
+1.1905e-02 -2.6620e-02 2.6620e-02 1.4286e-01 2.6620e-02 -2.6620e-02 1.1905e-02
+ 0. 0. 0. 2.6620e-02 3.5714e-01 5.9524e-02 -2.6620e-02
+ 0. 0. 0. -2.6620e-02 5.9524e-02 3.5714e-01 2.6620e-02
+ 0. 0. 0. 1.1905e-02 -2.6620e-02 2.6620e-02 7.1429e-02
+
+DEAL::Laplace matrix for FE_Q<1>(3):
+4.3333e+00 -1.6667e-01 -4.8784e+00 7.1175e-01
+-1.6667e-01 4.3333e+00 7.1175e-01 -4.8784e+00
+-4.8784e+00 7.1175e-01 8.3333e+00 -4.1667e+00
+7.1175e-01 -4.8784e+00 -4.1667e+00 8.3333e+00
+
+DEAL::Laplace matrix for FE_Q<1>(3) renumbered with FEEvaluation:
+4.3333e+00 -4.8784e+00 7.1175e-01 -1.6667e-01
+-4.8784e+00 8.3333e+00 -4.1667e+00 7.1175e-01
+7.1175e-01 -4.1667e+00 8.3333e+00 -4.8784e+00
+-1.6667e-01 7.1175e-01 -4.8784e+00 4.3333e+00
+