1 #include <deal.II/base/conditional_ostream.h> 2 #include <deal.II/base/parameter_handler.h> 4 #include <deal.II/base/qprojector.h> 5 #include <deal.II/base/geometry_info.h> 7 #include <deal.II/grid/tria.h> 9 #include <deal.II/fe/fe_dgq.h> 10 #include <deal.II/fe/fe_dgp.h> 11 #include <deal.II/fe/fe_system.h> 12 #include <deal.II/fe/mapping_fe_field.h> 13 #include <deal.II/fe/mapping_q1_eulerian.h> 16 #include <deal.II/dofs/dof_handler.h> 18 #include <deal.II/hp/q_collection.h> 19 #include <deal.II/hp/mapping_collection.h> 20 #include <deal.II/hp/fe_values.h> 22 #include <deal.II/lac/vector.h> 23 #include <deal.II/lac/sparsity_pattern.h> 24 #include <deal.II/lac/trilinos_sparse_matrix.h> 25 #include <deal.II/lac/trilinos_vector.h> 26 #include <deal.II/lac/identity_matrix.h> 28 #include <Epetra_RowMatrixTransposer.h> 31 #include "ADTypes.hpp" 33 #include <CoDiPack/include/codi.hpp> 35 #include "operators.h" 41 template <
int dim,
int n_faces,
typename real>
43 const int nstate_input,
44 const unsigned int max_degree_input,
45 const unsigned int grid_degree_input)
46 : max_degree(max_degree_input)
47 , max_grid_degree(grid_degree_input)
48 , nstate(nstate_input)
49 , max_grid_degree_check(grid_degree_input)
50 , mpi_communicator(MPI_COMM_WORLD)
51 , pcout(
std::cout, dealii::Utilities::MPI::this_mpi_process(mpi_communicator)==0)
54 template <
int dim,
int n_faces,
typename real>
56 const dealii::FullMatrix<double> &basis_x,
57 const dealii::FullMatrix<double> &basis_y,
58 const dealii::FullMatrix<double> &basis_z)
60 const unsigned int rows_x = basis_x.m();
61 const unsigned int columns_x = basis_x.n();
62 const unsigned int rows_y = basis_y.m();
63 const unsigned int columns_y = basis_y.n();
64 const unsigned int rows_z = basis_z.m();
65 const unsigned int columns_z = basis_z.n();
69 if constexpr (dim==2){
70 dealii::FullMatrix<double> tens_prod(rows_x * rows_y, columns_x * columns_y);
71 for(
unsigned int jdof=0; jdof<rows_y; jdof++){
72 for(
unsigned int kdof=0; kdof<rows_x; kdof++){
73 for(
unsigned int ndof=0; ndof<columns_y; ndof++){
74 for(
unsigned int odof=0; odof<columns_x; odof++){
75 const unsigned int index_row = rows_x * jdof + kdof;
76 const unsigned int index_col = columns_x * ndof + odof;
77 tens_prod[index_row][index_col] = basis_x[kdof][odof] * basis_y[jdof][ndof];
84 if constexpr (dim==3){
85 dealii::FullMatrix<double> tens_prod(rows_x * rows_y * rows_z, columns_x * columns_y * columns_z);
86 for(
unsigned int idof=0; idof<rows_z; idof++){
87 for(
unsigned int jdof=0; jdof<rows_y; jdof++){
88 for(
unsigned int kdof=0; kdof<rows_x; kdof++){
89 for(
unsigned int mdof=0; mdof<columns_z; mdof++){
90 for(
unsigned int ndof=0; ndof<columns_y; ndof++){
91 for(
unsigned int odof=0; odof<columns_x; odof++){
92 const unsigned int index_row = rows_x * rows_y * idof + rows_x * jdof + kdof;
93 const unsigned int index_col = columns_x * columns_y * mdof + columns_x * ndof + odof;
94 tens_prod[index_row][index_col] = basis_x[kdof][odof] * basis_y[jdof][ndof] * basis_z[idof][mdof];
105 template <
int dim,
int n_faces,
typename real>
108 const dealii::FullMatrix<double> &basis_x,
109 const dealii::FullMatrix<double> &basis_y,
110 const dealii::FullMatrix<double> &basis_z)
113 const unsigned int rows_x = basis_x.m();
114 const unsigned int columns_x = basis_x.n();
115 const unsigned int rows_y = basis_y.m();
116 const unsigned int columns_y = basis_y.n();
117 const unsigned int rows_z = basis_z.m();
118 const unsigned int columns_z = basis_z.n();
120 const unsigned int rows_1state_x = rows_x /
nstate;
121 const unsigned int columns_1state_x = columns_x /
nstate;
122 const unsigned int rows_1state_y = rows_y /
nstate;
123 const unsigned int columns_1state_y = columns_y /
nstate;
124 const unsigned int rows_1state_z = rows_z /
nstate;
125 const unsigned int columns_1state_z = columns_z /
nstate;
127 const unsigned int rows_all_states = (dim == 1) ? rows_1state_x * nstate :
128 ( (dim == 2) ? rows_1state_x * rows_1state_y *
nstate :
129 rows_1state_x * rows_1state_y * rows_1state_z *
nstate);
130 const unsigned int columns_all_states = (dim == 1) ? columns_1state_x * nstate :
131 ( (dim == 2) ? columns_1state_x * columns_1state_y *
nstate :
132 columns_1state_x * columns_1state_y * columns_1state_z *
nstate);
133 dealii::FullMatrix<double> tens_prod(rows_all_states, columns_all_states);
136 for(
int istate=0; istate<
nstate; istate++){
137 dealii::FullMatrix<double> basis_x_1state(rows_1state_x, columns_1state_x);
138 dealii::FullMatrix<double> basis_y_1state(rows_1state_y, columns_1state_y);
139 dealii::FullMatrix<double> basis_z_1state(rows_1state_z, columns_1state_z);
140 for(
unsigned int r=0; r<rows_1state_x; r++){
141 for(
unsigned int c=0; c<columns_1state_x; c++){
142 basis_x_1state[r][c] = basis_x[istate*rows_1state_x + r][istate*columns_1state_x + c];
145 if constexpr(dim>=2){
146 for(
unsigned int r=0; r<rows_1state_y; r++){
147 for(
unsigned int c=0; c<columns_1state_y; c++){
148 basis_y_1state[r][c] = basis_y[istate*rows_1state_y + r][istate*columns_1state_y + c];
152 if constexpr(dim>=3){
153 for(
unsigned int r=0; r<rows_1state_z; r++){
154 for(
unsigned int c=0; c<columns_1state_z; c++){
155 basis_z_1state[r][c] = basis_z[istate*rows_1state_z + r][istate*columns_1state_z + c];
159 const unsigned int r1state = (dim == 1) ? rows_1state_x : ( (dim==2) ? rows_1state_x * rows_1state_y : rows_1state_x * rows_1state_y * rows_1state_z);
160 const unsigned int c1state = (dim == 1) ? columns_1state_x : ( (dim==2) ? columns_1state_x * columns_1state_y : columns_1state_x * columns_1state_y * columns_1state_z);
161 dealii::FullMatrix<double> tens_prod_1state(r1state, c1state);
162 tens_prod_1state =
tensor_product(basis_x_1state, basis_y_1state, basis_z_1state);
163 for(
unsigned int r=0; r<r1state; r++){
164 for(
unsigned int c=0; c<c1state; c++){
165 tens_prod[istate*r1state + r][istate*c1state + c] = tens_prod_1state[r][c];
172 template <
int dim,
int n_faces,
typename real>
187 template <
int dim,
int n_faces,
typename real>
189 const int nstate_input,
190 const unsigned int max_degree_input,
191 const unsigned int grid_degree_input)
195 template <
int dim,
int n_faces,
typename real>
197 const std::vector<real> &input_vect,
198 std::vector<real> &output_vect,
199 const dealii::FullMatrix<double> &basis_x,
200 const dealii::FullMatrix<double> &basis_y,
201 const dealii::FullMatrix<double> &basis_z,
206 const unsigned int rows_x = basis_x.m();
207 const unsigned int rows_y = basis_y.m();
208 const unsigned int rows_z = basis_z.m();
209 const unsigned int columns_x = basis_x.n();
210 const unsigned int columns_y = basis_y.n();
211 const unsigned int columns_z = basis_z.n();
212 if constexpr (dim == 1){
213 assert(rows_x == output_vect.size());
214 assert(columns_x == input_vect.size());
216 if constexpr (dim == 2){
217 assert(rows_x * rows_y == output_vect.size());
218 assert(columns_x * columns_y == input_vect.size());
220 if constexpr (dim == 3){
221 assert(rows_x * rows_y * rows_z == output_vect.size());
222 assert(columns_x * columns_y * columns_z == input_vect.size());
225 if constexpr (dim==1){
226 for(
unsigned int iquad=0; iquad<rows_x; iquad++){
228 output_vect[iquad] = 0.0;
229 for(
unsigned int jquad=0; jquad<columns_x; jquad++){
230 output_vect[iquad] += factor * basis_x[iquad][jquad] * input_vect[jquad];
234 if constexpr (dim==2){
236 dealii::FullMatrix<double> input_mat(columns_x, columns_y);
237 for(
unsigned int idof=0; idof<columns_y; idof++){
238 for(
unsigned int jdof=0; jdof<columns_x; jdof++){
239 input_mat[jdof][idof] = input_vect[idof * columns_x + jdof];
242 dealii::FullMatrix<double> temp(rows_x, columns_y);
243 basis_x.mmult(temp, input_mat);
244 dealii::FullMatrix<double> output_mat(rows_y, rows_x);
245 basis_y.mTmult(output_mat, temp);
247 for(
unsigned int iquad=0; iquad<rows_y; iquad++){
248 for(
unsigned int jquad=0; jquad<rows_x; jquad++){
250 output_vect[iquad * rows_x + jquad] += factor * output_mat[iquad][jquad];
252 output_vect[iquad * rows_x + jquad] = factor * output_mat[iquad][jquad];
257 if constexpr (dim==3){
259 dealii::FullMatrix<double> input_mat(columns_x, columns_y * columns_z);
260 for(
unsigned int idof=0; idof<columns_z; idof++){
261 for(
unsigned int jdof=0; jdof<columns_y; jdof++){
262 for(
unsigned int kdof=0; kdof<columns_x; kdof++){
263 const unsigned int dof_index = idof * columns_x * columns_y + jdof * columns_x + kdof;
264 input_mat[kdof][idof * columns_y + jdof] = input_vect[dof_index];
268 dealii::FullMatrix<double> temp(rows_x, columns_y * columns_z);
269 basis_x.mmult(temp, input_mat);
271 dealii::FullMatrix<double> temp2(columns_y, rows_x * columns_z);
272 for(
unsigned int iquad=0; iquad<rows_x; iquad++){
273 for(
unsigned int idof=0; idof<columns_z; idof++){
274 for(
unsigned int jdof=0; jdof<columns_y; jdof++){
275 temp2[jdof][iquad * columns_z + idof] = temp[iquad][idof * columns_y + jdof];
279 dealii::FullMatrix<double> temp3(rows_y, rows_x * columns_z);
280 basis_y.mmult(temp3, temp2);
281 dealii::FullMatrix<double> temp4(columns_z, rows_x * rows_y);
283 for(
unsigned int iquad=0; iquad<rows_x; iquad++){
284 for(
unsigned int idof=0; idof<columns_z; idof++){
285 for(
unsigned int jquad=0; jquad<rows_y; jquad++){
286 temp4[idof][iquad * rows_y + jquad] = temp3[jquad][iquad * columns_z + idof];
290 dealii::FullMatrix<double> output_mat(rows_z, rows_x * rows_y);
291 basis_z.mmult(output_mat, temp4);
293 for(
unsigned int iquad=0; iquad<rows_z; iquad++){
294 for(
unsigned int jquad=0; jquad<rows_y; jquad++){
295 for(
unsigned int kquad=0; kquad<rows_x; kquad++){
296 const unsigned int quad_index = iquad * rows_x * rows_y + jquad * rows_x + kquad;
298 output_vect[quad_index] += factor * output_mat[iquad][kquad * rows_y + jquad];
300 output_vect[quad_index] = factor * output_mat[iquad][kquad * rows_y + jquad];
307 template <
int dim,
int n_faces,
typename real>
309 const std::vector<real> &input_vect,
310 std::vector<real> &output_vect,
311 const dealii::FullMatrix<double> &basis_x,
315 this->
matrix_vector_mult(input_vect, output_vect, basis_x, basis_x, basis_x, adding, factor);
318 template <
int dim,
int n_faces,
typename real>
320 const unsigned int face_number,
321 const std::vector<real> &input_vect,
322 std::vector<real> &output_vect,
323 const std::array<dealii::FullMatrix<double>,2> &basis_surf,
324 const dealii::FullMatrix<double> &basis_vol,
329 this->
matrix_vector_mult(input_vect, output_vect, basis_surf[0], basis_vol, basis_vol, adding, factor);
331 this->
matrix_vector_mult(input_vect, output_vect, basis_surf[1], basis_vol, basis_vol, adding, factor);
333 this->
matrix_vector_mult(input_vect, output_vect, basis_vol, basis_surf[0], basis_vol, adding, factor);
335 this->
matrix_vector_mult(input_vect, output_vect, basis_vol, basis_surf[1], basis_vol, adding, factor);
337 this->
matrix_vector_mult(input_vect, output_vect, basis_vol, basis_vol, basis_surf[0], adding, factor);
339 this->
matrix_vector_mult(input_vect, output_vect, basis_vol, basis_vol, basis_surf[1], adding, factor);
343 template <
int dim,
int n_faces,
typename real>
345 const unsigned int face_number,
346 const std::vector<real> &input_vect,
347 const std::vector<real> &weight_vect,
348 std::vector<real> &output_vect,
349 const std::array<dealii::FullMatrix<double>,2> &basis_surf,
350 const dealii::FullMatrix<double> &basis_vol,
355 this->
inner_product(input_vect, weight_vect, output_vect, basis_surf[0], basis_vol, basis_vol, adding, factor);
357 this->
inner_product(input_vect, weight_vect, output_vect, basis_surf[1], basis_vol, basis_vol, adding, factor);
359 this->
inner_product(input_vect, weight_vect, output_vect, basis_vol, basis_surf[0], basis_vol, adding, factor);
361 this->
inner_product(input_vect, weight_vect, output_vect, basis_vol, basis_surf[1], basis_vol, adding, factor);
363 this->
inner_product(input_vect, weight_vect, output_vect, basis_vol, basis_vol, basis_surf[0], adding, factor);
365 this->
inner_product(input_vect, weight_vect, output_vect, basis_vol, basis_vol, basis_surf[1], adding, factor);
368 template <
int dim,
int n_faces,
typename real>
370 const dealii::Tensor<1,dim,std::vector<real>> &input_vect,
371 std::vector<real> &output_vect,
372 const dealii::FullMatrix<double> &basis,
373 const dealii::FullMatrix<double> &gradient_basis)
377 gradient_basis, gradient_basis, gradient_basis);
380 template <
int dim,
int n_faces,
typename real>
382 const dealii::Tensor<1,dim,std::vector<real>> &input_vect,
383 std::vector<real> &output_vect,
384 const dealii::FullMatrix<double> &basis_x,
385 const dealii::FullMatrix<double> &basis_y,
386 const dealii::FullMatrix<double> &basis_z,
387 const dealii::FullMatrix<double> &gradient_basis_x,
388 const dealii::FullMatrix<double> &gradient_basis_y,
389 const dealii::FullMatrix<double> &gradient_basis_z)
391 for(
int idim=0; idim<dim;idim++){
413 template <
int dim,
int n_faces,
typename real>
415 const std::vector<real> &input_vect,
416 dealii::Tensor<1,dim,std::vector<real>> &output_vect,
417 const dealii::FullMatrix<double> &basis,
418 const dealii::FullMatrix<double> &gradient_basis)
422 gradient_basis, gradient_basis, gradient_basis);
425 template <
int dim,
int n_faces,
typename real>
427 const std::vector<real> &input_vect,
428 dealii::Tensor<1,dim,std::vector<real>> &output_vect,
429 const dealii::FullMatrix<double> &basis_x,
430 const dealii::FullMatrix<double> &basis_y,
431 const dealii::FullMatrix<double> &basis_z,
432 const dealii::FullMatrix<double> &gradient_basis_x,
433 const dealii::FullMatrix<double> &gradient_basis_y,
434 const dealii::FullMatrix<double> &gradient_basis_z)
436 for(
int idim=0; idim<dim;idim++){
459 template <
int dim,
int n_faces,
typename real>
461 const std::vector<real> &input_vect,
462 const std::vector<real> &weight_vect,
463 std::vector<real> &output_vect,
464 const dealii::FullMatrix<double> &basis_x,
465 const dealii::FullMatrix<double> &basis_y,
466 const dealii::FullMatrix<double> &basis_z,
471 const unsigned int rows_x = basis_x.m();
472 const unsigned int rows_y = basis_y.m();
473 const unsigned int rows_z = basis_z.m();
474 const unsigned int columns_x = basis_x.n();
475 const unsigned int columns_y = basis_y.n();
476 const unsigned int columns_z = basis_z.n();
479 if constexpr (dim == 1){
480 assert(rows_x == input_vect.size());
481 assert(columns_x == output_vect.size());
483 if constexpr (dim == 2){
484 assert(rows_x * rows_y == input_vect.size());
485 assert(columns_x * columns_y == output_vect.size());
487 if constexpr (dim == 3){
488 assert(rows_x * rows_y * rows_z == input_vect.size());
489 assert(columns_x * columns_y * columns_z == output_vect.size());
491 assert(weight_vect.size() == input_vect.size());
493 dealii::FullMatrix<double> basis_x_trans(columns_x, rows_x);
494 dealii::FullMatrix<double> basis_y_trans(columns_y, rows_y);
495 dealii::FullMatrix<double> basis_z_trans(columns_z, rows_z);
499 for(
unsigned int row=0; row<rows_x; row++){
500 for(
unsigned int col=0; col<columns_x; col++){
501 basis_x_trans[col][row] = basis_x[row][col];
504 for(
unsigned int row=0; row<rows_y; row++){
505 for(
unsigned int col=0; col<columns_y; col++){
506 basis_y_trans[col][row] = basis_y[row][col];
509 for(
unsigned int row=0; row<rows_z; row++){
510 for(
unsigned int col=0; col<columns_z; col++){
511 basis_z_trans[col][row] = basis_z[row][col];
515 std::vector<double> new_input_vect(input_vect.size());
516 for(
unsigned int iquad=0; iquad<input_vect.size(); iquad++){
517 new_input_vect[iquad] = input_vect[iquad] * weight_vect[iquad];
520 this->
matrix_vector_mult(new_input_vect, output_vect, basis_x_trans, basis_y_trans, basis_z_trans, adding, factor);
523 template <
int dim,
int n_faces,
typename real>
525 const std::vector<real> &input_vect,
526 const std::vector<real> &weight_vect,
527 std::vector<real> &output_vect,
528 const dealii::FullMatrix<double> &basis_x,
532 this->
inner_product(input_vect, weight_vect, output_vect, basis_x, basis_x, basis_x, adding, factor);
535 template <
int dim,
int n_faces,
typename real>
537 const dealii::Tensor<1,dim,dealii::FullMatrix<real>> &input_mat,
538 std::vector<real> &output_vect,
539 const std::vector<real> &weights,
540 const dealii::FullMatrix<double> &basis,
541 const double scaling)
543 assert(input_mat[0].m() == output_vect.size());
545 dealii::FullMatrix<double> output_mat(input_mat[0].m(), input_mat[0].n());
546 for(
int idim=0; idim<dim; idim++){
548 if constexpr(dim==1){
549 for(
unsigned int row=0; row<input_mat[0].m(); row++){
550 for(
unsigned int col=0; col<basis.m(); col++){
551 const unsigned int col_index = col;
552 output_vect[row] += scaling * output_mat[row][col_index];
556 if constexpr(dim==2){
557 const unsigned int size_1D = basis.m();
558 for(
unsigned int irow=0; irow<size_1D; irow++){
559 for(
unsigned int jrow=0; jrow<size_1D; jrow++){
560 const unsigned int row_index = irow * size_1D + jrow;
561 for(
unsigned int col=0; col<size_1D; col++){
563 const unsigned int col_index = col + irow * size_1D;
564 output_vect[row_index] += scaling * output_mat[row_index][col_index];
567 const unsigned int col_index = col * size_1D + jrow;
568 output_vect[row_index] += scaling * output_mat[row_index][col_index];
574 if constexpr(dim==3){
575 const unsigned int size_1D = basis.m();
576 for(
unsigned int irow=0; irow<size_1D; irow++){
577 for(
unsigned int jrow=0; jrow<size_1D; jrow++){
578 for(
unsigned int krow=0; krow<size_1D; krow++){
579 const unsigned int row_index = irow * size_1D * size_1D + jrow * size_1D + krow;
580 for(
unsigned int col=0; col<size_1D; col++){
582 const unsigned int col_index = col + irow * size_1D * size_1D + jrow * size_1D;
583 output_vect[row_index] += scaling * output_mat[row_index][col_index];
586 const unsigned int col_index = col * size_1D + krow + irow * size_1D * size_1D;
587 output_vect[row_index] += scaling * output_mat[row_index][col_index];
590 const unsigned int col_index = col * size_1D * size_1D + krow + jrow * size_1D;
591 output_vect[row_index] += scaling * output_mat[row_index][col_index];
601 template <
int dim,
int n_faces,
typename real>
603 const dealii::FullMatrix<real> &input_mat,
604 std::vector<real> &output_vect_vol,
605 std::vector<real> &output_vect_surf,
606 const std::vector<real> &weights,
607 const std::array<dealii::FullMatrix<double>,2> &surf_basis,
608 const unsigned int iface,
609 const unsigned int dim_not_zero,
610 const double scaling)
612 assert(input_mat.n() == output_vect_vol.size());
613 assert(input_mat.m() == output_vect_surf.size());
614 const unsigned int iface_1D = iface % 2;
616 dealii::FullMatrix<double> output_mat(input_mat.m(), input_mat.n());
618 if constexpr(dim==1){
619 for(
unsigned int row=0; row<surf_basis[iface_1D].m(); row++){
620 for(
unsigned int col=0; col<surf_basis[iface_1D].n(); col++){
621 output_vect_vol[col] += scaling
622 * output_mat[row][col];
623 output_vect_surf[row] -= scaling
624 * output_mat[row][col];
628 if constexpr(dim==2){
629 const unsigned int size_1D_row = surf_basis[iface_1D].m();
630 const unsigned int size_1D_col = surf_basis[iface_1D].n();
631 for(
unsigned int irow=0; irow<size_1D_col; irow++){
632 for(
unsigned int jrow=0; jrow<size_1D_row; jrow++){
633 const unsigned int row_index = irow * size_1D_row + jrow;
634 for(
unsigned int col=0; col<size_1D_col; col++){
636 const unsigned int col_index = col + irow * size_1D_col;
637 output_vect_vol[col_index] += scaling * output_mat[row_index][col_index];
638 output_vect_surf[row_index] -= scaling * output_mat[row_index][col_index];
641 const unsigned int col_index = col * size_1D_col + irow;
642 output_vect_vol[col_index] += scaling * output_mat[row_index][col_index];
643 output_vect_surf[row_index] -= scaling * output_mat[row_index][col_index];
649 if constexpr(dim==3){
650 const unsigned int size_1D_row = surf_basis[iface_1D].m();
651 const unsigned int size_1D_col = surf_basis[iface_1D].n();
652 for(
unsigned int irow=0; irow<size_1D_col; irow++){
653 for(
unsigned int jrow=0; jrow<size_1D_col; jrow++){
654 for(
unsigned int krow=0; krow<size_1D_row; krow++){
655 const unsigned int row_index = irow * size_1D_row * size_1D_col + jrow * size_1D_row + krow;
656 for(
unsigned int col=0; col<size_1D_col; col++){
658 const unsigned int col_index = col + irow * size_1D_col * size_1D_col + jrow * size_1D_col;
659 output_vect_vol[col_index] += scaling * output_mat[row_index][col_index];
660 output_vect_surf[row_index] -= scaling * output_mat[row_index][col_index];
663 const unsigned int col_index = col * size_1D_col + jrow + irow * size_1D_col * size_1D_col;
664 output_vect_vol[col_index] += scaling * output_mat[row_index][col_index];
665 output_vect_surf[row_index] -= scaling * output_mat[row_index][col_index];
668 const unsigned int col_index = col * size_1D_col * size_1D_col + jrow + irow * size_1D_col;
669 output_vect_vol[col_index] += scaling * output_mat[row_index][col_index];
670 output_vect_surf[row_index] -= scaling * output_mat[row_index][col_index];
681 template <
int dim,
int n_faces,
typename real>
683 const dealii::FullMatrix<real> &input_mat,
684 dealii::FullMatrix<real> &output_mat,
685 const dealii::FullMatrix<double> &basis,
686 const std::vector<real> &weights,
689 assert(input_mat.size() == output_mat.size());
690 const unsigned int size = basis.n();
691 assert(size == weights.size());
693 if constexpr(dim == 1){
696 if constexpr(dim == 2){
699 const unsigned int rows = basis.m();
700 assert(rows == input_mat.m());
702 for(
unsigned int idiag=0; idiag<size; idiag++){
703 dealii::FullMatrix<double> local_block(rows, size);
704 std::vector<unsigned int> row_index(rows);
705 std::vector<unsigned int> col_index(size);
707 std::iota(row_index.begin(), row_index.end(), idiag*rows);
708 std::iota(col_index.begin(), col_index.end(), idiag*size);
710 local_block.extract_submatrix_from(input_mat, row_index, col_index);
711 dealii::FullMatrix<double> local_Hadamard(rows, size);
714 local_Hadamard *= weights[idiag];
716 local_Hadamard.scatter_matrix_to(row_index, col_index, output_mat);
720 for(
unsigned int idiag=0; idiag<rows; idiag++){
721 const unsigned int row_index = idiag * size;
722 for(
unsigned int jdiag=0; jdiag<size; jdiag++){
723 const unsigned int col_index = jdiag * size;
724 for(
unsigned int kdiag=0; kdiag<size; kdiag++){
725 output_mat[row_index + kdiag][col_index + kdiag] = basis[idiag][jdiag]
726 * input_mat[row_index + kdiag][col_index + kdiag]
734 if constexpr(dim == 3){
735 const unsigned int rows = basis.m();
737 unsigned int kdiag=0;
738 for(
unsigned int idiag=0; idiag< size * size; idiag++){
739 if(kdiag==size) kdiag = 0;
740 dealii::FullMatrix<double> local_block(rows, size);
741 std::vector<unsigned int> row_index(rows);
742 std::vector<unsigned int> col_index(size);
744 std::iota(row_index.begin(), row_index.end(), idiag*rows);
745 std::iota(col_index.begin(), col_index.end(), idiag*size);
747 local_block.extract_submatrix_from(input_mat, row_index, col_index);
748 dealii::FullMatrix<double> local_Hadamard(rows, size);
751 local_Hadamard *= weights[kdiag];
753 const unsigned int jdiag = idiag / size;
754 local_Hadamard *= weights[jdiag];
756 local_Hadamard.scatter_matrix_to(row_index, col_index, output_mat);
760 for(
unsigned int zdiag=0; zdiag<size; zdiag++){
761 for(
unsigned int idiag=0; idiag<rows; idiag++){
762 const unsigned int row_index = zdiag * size * rows + idiag * size;
763 for(
unsigned int jdiag=0; jdiag<size; jdiag++){
764 const unsigned int col_index = zdiag * size * size + jdiag * size;
765 for(
unsigned int kdiag=0; kdiag<size; kdiag++){
766 output_mat[row_index + kdiag][col_index + kdiag] = basis[idiag][jdiag]
767 * input_mat[row_index + kdiag][col_index + kdiag]
776 for(
unsigned int row_block=0; row_block<rows; row_block++){
777 for(
unsigned int col_block=0; col_block<size; col_block++){
778 const unsigned int row_index = row_block * size * size;
779 const unsigned int col_index = col_block * size * size;
780 for(
unsigned int jdiag=0; jdiag<size; jdiag++){
781 for(
unsigned int kdiag=0; kdiag<size; kdiag++){
782 output_mat[row_index + jdiag * size + kdiag][col_index + jdiag * size + kdiag] = basis[row_block][col_block]
783 * input_mat[row_index + jdiag * size + kdiag][col_index + jdiag * size + kdiag]
794 template <
int dim,
int n_faces,
typename real>
796 const dealii::FullMatrix<real> &input_mat1,
797 const dealii::FullMatrix<real> &input_mat2,
798 dealii::FullMatrix<real> &output_mat)
800 const unsigned int rows = input_mat1.m();
801 const unsigned int columns = input_mat1.n();
802 assert(rows == input_mat2.m());
803 assert(columns == input_mat2.n());
805 for(
unsigned int irow=0; irow<rows; irow++){
806 for(
unsigned int icol=0; icol<columns; icol++){
807 output_mat[irow][icol] = input_mat1[irow][icol]
808 * input_mat2[irow][icol];
813 template <
int dim,
int n_faces,
typename real>
815 const unsigned int rows_size,
816 const unsigned int columns_size,
817 std::vector<std::array<unsigned int,dim>> &rows,
818 std::vector<std::array<unsigned int,dim>> &columns)
821 if constexpr(dim == 1){
822 for(
unsigned int irow=0; irow<rows_size; irow++){
823 for(
unsigned int icol=0; icol<columns_size; icol++){
824 const unsigned int array_index = irow * rows_size + icol;
825 rows[array_index][0] = irow;
826 columns[array_index][0] = icol;
830 if constexpr(dim == 2){
831 for(
unsigned int idiag=0; idiag<rows_size; idiag++){
832 for(
unsigned int jdiag=0; jdiag<columns_size; jdiag++){
833 for(
unsigned int kdiag=0; kdiag<columns_size; kdiag++){
834 const unsigned int array_index = idiag * rows_size * columns_size + jdiag * columns_size + kdiag;
835 const unsigned int row_index = idiag * rows_size;
836 rows[array_index][0] = row_index + jdiag;
837 rows[array_index][1] = row_index + jdiag;
839 const unsigned int col_index_0 = idiag * columns_size;
840 columns[array_index][0] = col_index_0 + kdiag;
842 const unsigned int col_index_1 = kdiag * columns_size;
843 columns[array_index][1] = col_index_1 + jdiag;
848 if constexpr(dim == 3){
849 for(
unsigned int idiag=0; idiag<rows_size; idiag++){
850 for(
unsigned int jdiag=0; jdiag<columns_size; jdiag++){
851 for(
unsigned int kdiag=0; kdiag<columns_size; kdiag++){
852 for(
unsigned int ldiag=0; ldiag<columns_size; ldiag++){
853 const unsigned int array_index = idiag * rows_size * columns_size * columns_size
854 + jdiag * columns_size * columns_size
855 + kdiag * columns_size
857 const unsigned int row_index = idiag * rows_size * columns_size
858 + jdiag * columns_size;
859 rows[array_index][0] = row_index + kdiag;
860 rows[array_index][1] = row_index + kdiag;
861 rows[array_index][2] = row_index + kdiag;
863 const unsigned int col_index_0 = idiag * columns_size * columns_size
864 + jdiag * columns_size;
865 columns[array_index][0] = col_index_0 + ldiag;
867 const unsigned int col_index_1 = ldiag * columns_size;
868 columns[array_index][1] = col_index_1 + kdiag + idiag * columns_size * columns_size;
870 const unsigned int col_index_2 = ldiag * columns_size * columns_size;
871 columns[array_index][2] = col_index_2 + kdiag + jdiag * columns_size;
878 template <
int dim,
int n_faces,
typename real>
880 const unsigned int rows_size_1D,
881 const unsigned int columns_size_1D,
882 const std::vector<std::array<unsigned int,dim>> &rows,
883 const std::vector<std::array<unsigned int,dim>> &columns,
884 const dealii::FullMatrix<double> &basis,
885 const std::vector<double> &weights,
886 std::array<dealii::FullMatrix<double>,dim> &basis_sparse)
888 if constexpr(dim == 1){
889 for(
unsigned int irow=0; irow<rows_size_1D; irow++){
890 for(
unsigned int icol=0; icol<columns_size_1D; icol++){
891 basis_sparse[0][irow][icol] = basis[irow][icol];
895 if constexpr(dim == 2){
896 const unsigned int total_size = rows.size();
897 for(
unsigned int index=0, counter=0; index<total_size; index++, counter++){
898 if(counter == columns_size_1D){
902 basis_sparse[0][rows[index][0]][counter] = basis[rows[index][0]%rows_size_1D][columns[index][0]%columns_size_1D]
903 * weights[rows[index][1]/columns_size_1D];
905 basis_sparse[1][rows[index][1]][counter] = basis[rows[index][1]/rows_size_1D][columns[index][1]/columns_size_1D]
906 * weights[rows[index][0]%columns_size_1D];
909 if constexpr(dim == 3){
910 const unsigned int total_size = rows.size();
911 for(
unsigned int index=0, counter=0; index<total_size; index++, counter++){
912 if(counter == columns_size_1D){
916 basis_sparse[0][rows[index][0]][counter] = basis[rows[index][0]%rows_size_1D][columns[index][0]%columns_size_1D]
917 * weights[(rows[index][1]/columns_size_1D)%columns_size_1D]
918 * weights[rows[index][2]/columns_size_1D/columns_size_1D];
920 basis_sparse[1][rows[index][1]][counter] = basis[(rows[index][1]/rows_size_1D)%rows_size_1D][(columns[index][1]/columns_size_1D)%columns_size_1D]
921 * weights[rows[index][0]%columns_size_1D]
922 * weights[rows[index][2]/columns_size_1D/columns_size_1D];
924 basis_sparse[2][rows[index][2]][counter] = basis[rows[index][2]/rows_size_1D/rows_size_1D][columns[index][2]/columns_size_1D/columns_size_1D]
925 * weights[rows[index][0]%columns_size_1D]
926 * weights[(rows[index][1]/columns_size_1D)%columns_size_1D];
931 template <
int dim,
int n_faces,
typename real>
933 const unsigned int rows_size,
934 const unsigned int columns_size,
935 std::vector<unsigned int> &rows,
936 std::vector<unsigned int> &columns,
937 const int dim_not_zero)
940 if constexpr(dim == 1){
941 for(
unsigned int irow=0; irow<rows_size; irow++){
942 for(
unsigned int icol=0; icol<columns_size; icol++){
943 const unsigned int array_index = irow * columns_size + icol;
944 rows[array_index] = irow;
945 columns[array_index] = icol;
949 if constexpr(dim == 2){
950 for(
unsigned int idiag=0; idiag<rows_size; idiag++){
951 for(
unsigned int jdiag=0; jdiag<columns_size; jdiag++){
952 const unsigned int array_index = idiag * columns_size + jdiag ;
954 const unsigned int row_index = idiag;
955 rows[array_index] = row_index;
957 if(dim_not_zero == 0){
958 const unsigned int col_index_0 = idiag * columns_size;
959 columns[array_index] = col_index_0 + jdiag;
962 if(dim_not_zero == 1){
963 const unsigned int col_index_1 = jdiag * columns_size;
964 columns[array_index] = col_index_1 + idiag;
969 if constexpr(dim == 3){
970 for(
unsigned int idiag=0; idiag<columns_size; idiag++){
971 for(
unsigned int jdiag=0; jdiag<columns_size; jdiag++){
972 for(
unsigned int kdiag=0; kdiag<columns_size; kdiag++){
973 const unsigned int array_index = idiag * columns_size * columns_size
974 + jdiag * columns_size
976 const unsigned int row_index = idiag * columns_size + jdiag;
977 rows[array_index] = row_index;
979 if(dim_not_zero == 0){
980 const unsigned int col_index_0 = idiag * columns_size * columns_size
981 + jdiag * columns_size;
982 columns[array_index] = col_index_0 + kdiag;
985 if(dim_not_zero == 1){
986 const unsigned int col_index_1 = kdiag * columns_size;
987 columns[array_index] = col_index_1 + jdiag + idiag * columns_size * columns_size;
990 if(dim_not_zero == 2){
991 const unsigned int col_index_2 = kdiag * columns_size * columns_size;
992 columns[array_index] = col_index_2 + jdiag + idiag * columns_size;
999 template <
int dim,
int n_faces,
typename real>
1001 const unsigned int rows_size,
1002 const unsigned int columns_size_1D,
1003 const std::vector<unsigned int> &rows,
1004 const std::vector<unsigned int> &columns,
1005 const dealii::FullMatrix<double> &basis,
1006 const std::vector<double> &weights,
1007 dealii::FullMatrix<double> &basis_sparse,
1008 const int dim_not_zero)
1010 if constexpr(dim == 1){
1011 for(
unsigned int irow=0; irow<rows_size; irow++){
1012 for(
unsigned int icol=0; icol<columns_size_1D; icol++){
1013 basis_sparse[irow][icol] = basis[irow][icol];
1017 if constexpr(dim == 2){
1018 const unsigned int total_size = rows.size();
1019 for(
unsigned int index=0, counter=0; index<total_size; index++, counter++){
1020 if(counter == columns_size_1D){
1024 if(dim_not_zero == 0){
1025 basis_sparse[rows[index]][counter] = basis[0][columns[index]%columns_size_1D]
1026 * weights[rows[index]%columns_size_1D];
1029 if(dim_not_zero == 1){
1030 basis_sparse[rows[index]][counter] = basis[0][columns[index]/columns_size_1D]
1031 * weights[rows[index]%columns_size_1D];
1035 if constexpr(dim == 3){
1036 const unsigned int total_size = rows.size();
1037 for(
unsigned int index=0, counter=0; index<total_size; index++, counter++){
1038 if(counter == columns_size_1D){
1042 if(dim_not_zero == 0){
1043 basis_sparse[rows[index]][counter] = basis[0][columns[index]%columns_size_1D]
1044 * weights[rows[index]%columns_size_1D]
1045 * weights[rows[index]/columns_size_1D];
1048 if(dim_not_zero == 1){
1049 basis_sparse[rows[index]][counter] = basis[0][(columns[index]/columns_size_1D)%columns_size_1D]
1050 * weights[rows[index]%columns_size_1D]
1051 * weights[rows[index]/columns_size_1D];
1054 if(dim_not_zero == 2){
1055 basis_sparse[rows[index]][counter] = basis[0][columns[index]/columns_size_1D/columns_size_1D]
1056 * weights[rows[index]%columns_size_1D]
1057 * weights[rows[index]/columns_size_1D];
1070 template <
int dim,
int n_faces,
typename real>
1072 const int nstate_input,
1073 const unsigned int max_degree_input,
1074 const unsigned int grid_degree_input)
1081 template <
int dim,
int n_faces,
typename real>
1083 const dealii::FESystem<1,1> &finite_element,
1084 const dealii::Quadrature<1> &quadrature)
1086 const unsigned int n_quad_pts = quadrature.size();
1087 const unsigned int n_dofs = finite_element.dofs_per_cell;
1091 for(
unsigned int iquad=0; iquad<n_quad_pts; iquad++){
1092 const dealii::Point<1> qpoint = quadrature.point(iquad);
1093 for(
unsigned int idof=0; idof<n_dofs; idof++){
1094 const int istate = finite_element.system_to_component_index(idof).first;
1096 this->
oneD_vol_operator[iquad][idof] = finite_element.shape_value_component(idof,qpoint,istate);
1101 template <
int dim,
int n_faces,
typename real>
1103 const dealii::FESystem<1,1> &finite_element,
1104 const dealii::Quadrature<1> &quadrature)
1106 const unsigned int n_quad_pts = quadrature.size();
1107 const unsigned int n_dofs = finite_element.dofs_per_cell;
1111 for(
unsigned int iquad=0; iquad<n_quad_pts; iquad++){
1112 const dealii::Point<1> qpoint = quadrature.point(iquad);
1113 for(
unsigned int idof=0; idof<n_dofs; idof++){
1114 const int istate = finite_element.system_to_component_index(idof).first;
1116 this->
oneD_grad_operator[iquad][idof] = finite_element.shape_grad_component(idof,qpoint,istate)[0];
1121 template <
int dim,
int n_faces,
typename real>
1123 const dealii::FESystem<1,1> &finite_element,
1124 const dealii::Quadrature<0> &face_quadrature)
1126 const unsigned int n_face_quad_pts = face_quadrature.size();
1127 const unsigned int n_dofs = finite_element.dofs_per_cell;
1128 const unsigned int n_faces_1D = n_faces / dim;
1130 for(
unsigned int iface=0; iface<n_faces_1D; iface++){
1134 const dealii::Quadrature<1> quadrature = dealii::QProjector<1>::project_to_face(dealii::ReferenceCell::get_hypercube(1),
1137 for(
unsigned int iquad=0; iquad<n_face_quad_pts; iquad++){
1138 const dealii::Point<1> qpoint = quadrature.point(iquad);
1139 for(
unsigned int idof=0; idof<n_dofs; idof++){
1140 const int istate = finite_element.system_to_component_index(idof).first;
1142 this->
oneD_surf_operator[iface][iquad][idof] = finite_element.shape_value_component(idof,qpoint,istate);
1148 template <
int dim,
int n_faces,
typename real>
1150 const dealii::FESystem<1,1> &finite_element,
1151 const dealii::Quadrature<0> &face_quadrature)
1153 const unsigned int n_face_quad_pts = face_quadrature.size();
1154 const unsigned int n_dofs = finite_element.dofs_per_cell;
1155 const unsigned int n_faces_1D = n_faces / dim;
1157 for(
unsigned int iface=0; iface<n_faces_1D; iface++){
1161 const dealii::Quadrature<1> quadrature = dealii::QProjector<1>::project_to_face(dealii::ReferenceCell::get_hypercube(1),
1164 for(
unsigned int iquad=0; iquad<n_face_quad_pts; iquad++){
1165 const dealii::Point<1> qpoint = quadrature.point(iquad);
1166 for(
unsigned int idof=0; idof<n_dofs; idof++){
1167 const int istate = finite_element.system_to_component_index(idof).first;
1169 this->
oneD_surf_grad_operator[iface][iquad][idof] = finite_element.shape_grad_component(idof,qpoint,istate)[0];
1175 template <
int dim,
int n_faces,
typename real>
1177 const int nstate_input,
1178 const unsigned int max_degree_input,
1179 const unsigned int grid_degree_input)
1186 template <
int dim,
int n_faces,
typename real>
1188 const dealii::FESystem<1,1> &finite_element,
1189 const dealii::Quadrature<1> &quadrature)
1191 const unsigned int n_quad_pts = quadrature.size();
1192 const unsigned int n_dofs = finite_element.dofs_per_cell;
1196 const std::vector<double> &quad_weights = quadrature.get_weights ();
1197 for(
unsigned int iquad=0; iquad<n_quad_pts; iquad++){
1198 const dealii::Point<1> qpoint = quadrature.point(iquad);
1199 for(
unsigned int idof=0; idof<n_dofs; idof++){
1200 const int istate = finite_element.system_to_component_index(idof).first;
1202 this->
oneD_vol_operator[iquad][idof] = quad_weights[iquad] * finite_element.shape_value_component(idof,qpoint,istate);
1207 template <
int dim,
int n_faces,
typename real>
1209 const int nstate_input,
1210 const unsigned int max_degree_input,
1211 const unsigned int grid_degree_input)
1218 template <
int dim,
int n_faces,
typename real>
1220 const dealii::FESystem<1,1> &finite_element,
1221 const dealii::Quadrature<1> &quadrature)
1223 const unsigned int n_quad_pts = quadrature.size();
1224 const unsigned int n_dofs = finite_element.dofs_per_cell;
1225 const std::vector<double> &quad_weights = quadrature.get_weights ();
1229 for (
unsigned int itest=0; itest<n_dofs; ++itest) {
1230 const int istate_test = finite_element.system_to_component_index(itest).first;
1231 for (
unsigned int itrial=itest; itrial<n_dofs; ++itrial) {
1232 const int istate_trial = finite_element.system_to_component_index(itrial).first;
1234 for (
unsigned int iquad=0; iquad<n_quad_pts; ++iquad) {
1235 const dealii::Point<1> qpoint = quadrature.point(iquad);
1237 finite_element.shape_value_component(itest,qpoint,istate_test)
1238 * finite_element.shape_value_component(itrial,qpoint,istate_trial)
1239 * quad_weights[iquad];
1244 if(istate_test==istate_trial) {
1252 template <
int dim,
int n_faces,
typename real>
1255 const unsigned int n_dofs,
const unsigned int n_quad_pts,
1257 const std::vector<double> &det_Jac,
1258 const std::vector<double> &quad_weights)
1261 const unsigned int n_shape_fns = n_dofs /
nstate;
1263 dealii::FullMatrix<double> mass_matrix_dim(n_dofs);
1264 dealii::FullMatrix<double> basis_dim(n_dofs);
1271 for(
int istate=0; istate<
nstate; istate++){
1272 for(
unsigned int itest=0; itest<n_shape_fns; ++itest){
1273 for (
unsigned int itrial=itest; itrial<n_shape_fns; ++itrial) {
1275 const unsigned int trial_index = istate*n_shape_fns + itrial;
1276 const unsigned int test_index = istate*n_shape_fns + itest;
1277 for (
unsigned int iquad=0; iquad<n_quad_pts; ++iquad) {
1278 value += basis_dim[iquad][test_index]
1279 * basis_dim[iquad][trial_index]
1281 * quad_weights[iquad];
1283 mass_matrix_dim[trial_index][test_index] = value;
1284 mass_matrix_dim[test_index][trial_index] = value;
1288 return mass_matrix_dim;
1291 template <
int dim,
int n_faces,
typename real>
1293 const int nstate_input,
1294 const unsigned int max_degree_input,
1295 const unsigned int grid_degree_input,
1296 const bool store_skew_symmetric_form_input)
1298 , store_skew_symmetric_form(store_skew_symmetric_form_input)
1304 template <
int dim,
int n_faces,
typename real>
1306 const dealii::FESystem<1,1> &finite_element,
1307 const dealii::Quadrature<1> &quadrature)
1309 const unsigned int n_quad_pts = quadrature.size();
1310 const unsigned int n_dofs = finite_element.dofs_per_cell;
1311 const std::vector<double> &quad_weights = quadrature.get_weights ();
1315 for(
unsigned int itest=0; itest<n_dofs; itest++){
1316 const int istate_test = finite_element.system_to_component_index(itest).first;
1317 for(
unsigned int idof=0; idof<n_dofs; idof++){
1318 const int istate = finite_element.system_to_component_index(idof).first;
1320 for(
unsigned int iquad=0; iquad<n_quad_pts; iquad++){
1321 const dealii::Point<1> qpoint = quadrature.point(iquad);
1322 value += finite_element.shape_value_component(itest,qpoint,istate_test)
1323 * finite_element.shape_grad_component(idof, qpoint, istate)[0]
1324 * quad_weights[iquad];
1326 if(istate == istate_test){
1335 for(
unsigned int idof=0; idof<n_dofs; idof++){
1336 for(
unsigned int jdof=0; jdof<n_dofs; jdof++){
1344 template <
int dim,
int n_faces,
typename real>
1346 const int nstate_input,
1347 const unsigned int max_degree_input,
1348 const unsigned int grid_degree_input)
1355 template <
int dim,
int n_faces,
typename real>
1357 const dealii::FESystem<1,1> &finite_element,
1358 const dealii::Quadrature<1> &quadrature)
1360 const unsigned int n_dofs = finite_element.dofs_per_cell;
1367 dealii::FullMatrix<double> inv_mass(n_dofs);
1373 template <
int dim,
int n_faces,
typename real>
1375 const int nstate_input,
1376 const unsigned int max_degree_input,
1377 const unsigned int grid_degree_input)
1384 template <
int dim,
int n_faces,
typename real>
1386 const dealii::FESystem<1,1> &finite_element,
1387 const dealii::Quadrature<1> &quadrature)
1389 const unsigned int n_dofs = finite_element.dofs_per_cell;
1393 for(
unsigned int idof=0; idof<n_dofs; idof++){
1400 for(
unsigned int idegree=0; idegree< this->
max_degree; idegree++){
1401 dealii::FullMatrix<double> derivative_p_temp(n_dofs, n_dofs);
1407 template <
int dim,
int n_faces,
typename real>
1409 const int nstate_input,
1410 const unsigned int max_degree_input,
1411 const unsigned int grid_degree_input,
1414 , FR_param_type(FR_param_input)
1422 template <
int dim,
int n_faces,
typename real>
1424 const unsigned int curr_cell_degree,
1429 double cp = pfact2/(pow(pfact,2));
1430 c = 2.0 * (curr_cell_degree+1)/( curr_cell_degree*((2.0*curr_cell_degree+1.0)*(pow(pfact*cp,2))));
1433 template <
int dim,
int n_faces,
typename real>
1435 const unsigned int curr_cell_degree,
1440 double cp = pfact2/(pow(pfact,2));
1441 c = 2.0 * (curr_cell_degree)/( (curr_cell_degree+1.0)*((2.0*curr_cell_degree+1.0)*(pow(pfact*cp,2))));
1444 template <
int dim,
int n_faces,
typename real>
1446 const unsigned int curr_cell_degree,
1451 double cp = pfact2/(pow(pfact,2));
1452 c = - 2.0 / ( pow((2.0*curr_cell_degree+1.0)*(pow(pfact*cp,2)),1.0));
1455 template <
int dim,
int n_faces,
typename real>
1457 const unsigned int curr_cell_degree,
1463 template <
int dim,
int n_faces,
typename real>
1465 const unsigned int curr_cell_degree,
1468 if(curr_cell_degree == 2){
1472 else if(curr_cell_degree == 3){
1475 else if(curr_cell_degree == 4){
1479 else if(curr_cell_degree == 5){
1483 this->
pcout <<
"ERROR: cPlus values are only defined for p=2 through p=5. Aborting..." << std::endl;
1488 c/=pow(pow(2.0,curr_cell_degree),2);
1491 template <
int dim,
int n_faces,
typename real>
1493 const unsigned int curr_cell_degree,
1521 template <
int dim,
int n_faces,
typename real>
1523 const dealii::FullMatrix<double> &local_Mass_Matrix,
1524 const dealii::FullMatrix<double> &pth_derivative,
1525 const unsigned int n_dofs,
1527 dealii::FullMatrix<double> &Flux_Reconstruction_operator)
1529 dealii::FullMatrix<double> derivative_p_temp(n_dofs);
1530 derivative_p_temp.add(c, pth_derivative);
1531 dealii::FullMatrix<double> Flux_Reconstruction_operator_temp(n_dofs);
1532 derivative_p_temp.Tmmult(Flux_Reconstruction_operator_temp, local_Mass_Matrix);
1533 Flux_Reconstruction_operator_temp.mmult(Flux_Reconstruction_operator, pth_derivative);
1537 template <
int dim,
int n_faces,
typename real>
1539 const dealii::FESystem<1,1> &finite_element,
1540 const dealii::Quadrature<1> &quadrature)
1542 const unsigned int n_dofs = finite_element.dofs_per_cell;
1555 template <
int dim,
int n_faces,
typename real>
1558 const unsigned int n_dofs,
1559 dealii::FullMatrix<double> &pth_deriv,
1560 dealii::FullMatrix<double> &mass_matrix)
1562 dealii::FullMatrix<double> Flux_Reconstruction_operator(n_dofs);
1563 dealii::FullMatrix<double> identity (dealii::IdentityMatrix(pth_deriv.m()));
1564 for(
int idim=0; idim<dim; idim++){
1565 dealii::FullMatrix<double> pth_deriv_dim(n_dofs);
1575 dealii::FullMatrix<double> derivative_p_temp(n_dofs);
1576 derivative_p_temp.add(
FR_param, pth_deriv_dim);
1577 dealii::FullMatrix<double> Flux_Reconstruction_operator_temp(n_dofs);
1578 derivative_p_temp.Tmmult(Flux_Reconstruction_operator_temp, mass_matrix);
1579 Flux_Reconstruction_operator_temp.mmult(Flux_Reconstruction_operator, pth_deriv_dim,
true);
1581 if constexpr (dim>=2){
1582 const int deriv_2p_loop = (dim==2) ? 1 : dim;
1583 double FR_param_sqrd = pow(
FR_param,2.0);
1584 for(
int idim=0; idim<deriv_2p_loop; idim++){
1585 dealii::FullMatrix<double> pth_deriv_dim(n_dofs);
1595 dealii::FullMatrix<double> derivative_p_temp(n_dofs);
1596 derivative_p_temp.add(FR_param_sqrd, pth_deriv_dim);
1597 dealii::FullMatrix<double> Flux_Reconstruction_operator_temp(n_dofs);
1598 derivative_p_temp.Tmmult(Flux_Reconstruction_operator_temp, mass_matrix);
1599 Flux_Reconstruction_operator_temp.mmult(Flux_Reconstruction_operator, pth_deriv_dim,
true);
1602 if constexpr (dim == 3){
1603 double FR_param_cubed = pow(
FR_param,3.0);
1604 dealii::FullMatrix<double> pth_deriv_dim(n_dofs);
1606 dealii::FullMatrix<double> derivative_p_temp(n_dofs);
1607 derivative_p_temp.add(FR_param_cubed, pth_deriv_dim);
1608 dealii::FullMatrix<double> Flux_Reconstruction_operator_temp(n_dofs);
1609 derivative_p_temp.Tmmult(Flux_Reconstruction_operator_temp, mass_matrix);
1610 Flux_Reconstruction_operator_temp.mmult(Flux_Reconstruction_operator, pth_deriv_dim,
true);
1613 return Flux_Reconstruction_operator;
1616 template <
int dim,
int n_faces,
typename real>
1618 const dealii::FullMatrix<double> &local_Mass_Matrix,
1620 const unsigned int n_dofs)
1622 dealii::FullMatrix<double> dim_FR_operator(n_dofs);
1623 if constexpr (dim == 1){
1627 dealii::FullMatrix<double> FR1(n_dofs);
1629 dealii::FullMatrix<double> FR2(n_dofs);
1631 dealii::FullMatrix<double> FR_cross1(n_dofs);
1633 dim_FR_operator.add(1.0, FR1, 1.0, FR2, 1.0, FR_cross1);
1635 if constexpr (dim == 3){
1636 dealii::FullMatrix<double> FR3(n_dofs);
1638 dealii::FullMatrix<double> FR_cross2(n_dofs);
1640 dealii::FullMatrix<double> FR_cross3(n_dofs);
1642 dealii::FullMatrix<double> FR_triple(n_dofs);
1644 dim_FR_operator.add(1.0, FR3, 1.0, FR_cross2, 1.0, FR_cross3);
1645 dim_FR_operator.add(1.0, FR_triple);
1647 return dim_FR_operator;
1652 template <
int dim,
int n_faces,
typename real>
1654 const int nstate_input,
1655 const unsigned int max_degree_input,
1656 const unsigned int grid_degree_input,
1659 , FR_param_aux_type(FR_param_aux_input)
1667 template <
int dim,
int n_faces,
typename real>
1669 const unsigned int curr_cell_degree,
1695 template <
int dim,
int n_faces,
typename real>
1697 const dealii::FESystem<1,1> &finite_element,
1698 const dealii::Quadrature<1> &quadrature)
1700 const unsigned int n_dofs = finite_element.dofs_per_cell;
1712 template <
int dim,
int n_faces,
typename real>
1714 const int nstate_input,
1715 const unsigned int max_degree_input,
1716 const unsigned int grid_degree_input)
1723 template <
int dim,
int n_faces,
typename real>
1725 const dealii::FullMatrix<double> &norm_matrix_inverse,
1726 const dealii::FullMatrix<double> &integral_vol_basis,
1727 dealii::FullMatrix<double> &volume_projection)
1729 norm_matrix_inverse.mTmult(volume_projection, integral_vol_basis);
1731 template <
int dim,
int n_faces,
typename real>
1733 const dealii::FESystem<1,1> &finite_element,
1734 const dealii::Quadrature<1> &quadrature)
1736 const unsigned int n_dofs = finite_element.dofs_per_cell;
1737 const unsigned int n_quad_pts = quadrature.size();
1742 dealii::FullMatrix<double> mass_inv(n_dofs);
1750 template <
int dim,
int n_faces,
typename real>
1752 const int nstate_input,
1753 const unsigned int max_degree_input,
1754 const unsigned int grid_degree_input,
1756 const bool store_transpose_input)
1758 , store_transpose(store_transpose_input)
1759 , FR_param_type(FR_param_input)
1765 template <
int dim,
int n_faces,
typename real>
1767 const dealii::FESystem<1,1> &finite_element,
1768 const dealii::Quadrature<1> &quadrature)
1770 const unsigned int n_dofs = finite_element.dofs_per_cell;
1771 const unsigned int n_quad_pts = quadrature.size();
1783 for(
unsigned int idof=0; idof<n_dofs; idof++){
1784 for(
unsigned int iquad=0; iquad<n_quad_pts; iquad++){
1790 template <
int dim,
int n_faces,
typename real>
1792 const int nstate_input,
1793 const unsigned int max_degree_input,
1794 const unsigned int grid_degree_input,
1796 const bool store_transpose_input)
1805 template <
int dim,
int n_faces,
typename real>
1807 const dealii::FESystem<1,1> &finite_element,
1808 const dealii::Quadrature<1> &quadrature)
1810 const unsigned int n_dofs = finite_element.dofs_per_cell;
1811 const unsigned int n_quad_pts = quadrature.size();
1823 for(
unsigned int idof=0; idof<n_dofs; idof++){
1824 for(
unsigned int iquad=0; iquad<n_quad_pts; iquad++){
1831 template <
int dim,
int n_faces,
typename real>
1833 const int nstate_input,
1834 const unsigned int max_degree_input,
1835 const unsigned int grid_degree_input,
1844 template <
int dim,
int n_faces,
typename real>
1846 const dealii::FESystem<1,1> &finite_element,
1847 const dealii::Quadrature<1> &quadrature)
1849 const unsigned int n_dofs = finite_element.dofs_per_cell;
1854 dealii::FullMatrix<double> FR_mass_matrix(n_dofs);
1862 template <
int dim,
int n_faces,
typename real>
1864 const int nstate_input,
1865 const unsigned int max_degree_input,
1866 const unsigned int grid_degree_input,
1875 template <
int dim,
int n_faces,
typename real>
1877 const dealii::FESystem<1,1> &finite_element,
1878 const dealii::Quadrature<1> &quadrature)
1880 const unsigned int n_dofs = finite_element.dofs_per_cell;
1885 dealii::FullMatrix<double> FR_mass_matrix(n_dofs);
1892 template <
int dim,
int n_faces,
typename real>
1894 const int nstate_input,
1895 const unsigned int max_degree_input,
1896 const unsigned int grid_degree_input,
1905 template <
int dim,
int n_faces,
typename real>
1907 const dealii::FESystem<1,1> &finite_element,
1908 const dealii::Quadrature<1> &quadrature)
1910 const unsigned int n_dofs = finite_element.dofs_per_cell;
1915 dealii::FullMatrix<double> FR_mass_matrix(n_dofs);
1922 template <
int dim,
int n_faces,
typename real>
1924 const int nstate_input,
1925 const unsigned int max_degree_input,
1926 const unsigned int grid_degree_input,
1935 template <
int dim,
int n_faces,
typename real>
1937 const dealii::FESystem<1,1> &finite_element,
1938 const dealii::Quadrature<1> &quadrature)
1940 const unsigned int n_dofs = finite_element.dofs_per_cell;
1945 dealii::FullMatrix<double> FR_mass_matrix(n_dofs);
1953 template <
int dim,
int n_faces,
typename real>
1955 const int nstate_input,
1956 const unsigned int max_degree_input,
1957 const unsigned int grid_degree_input)
1964 template <
int dim,
int n_faces,
typename real>
1966 const dealii::FESystem<1,1> &finite_element,
1967 const dealii::Quadrature<1> &quadrature)
1969 const unsigned int n_quad_pts = quadrature.size();
1970 const unsigned int n_dofs = finite_element.dofs_per_cell;
1974 const std::vector<double> &quad_weights = quadrature.get_weights ();
1975 for(
unsigned int itest=0; itest<n_dofs; itest++){
1976 const int istate_test = finite_element.system_to_component_index(itest).first;
1977 for(
unsigned int iquad=0; iquad<n_quad_pts; iquad++){
1978 const dealii::Point<1> qpoint = quadrature.point(iquad);
1979 this->
oneD_grad_operator[iquad][itest] = finite_element.shape_grad_component(itest, qpoint, istate_test)[0]
1980 * quad_weights[iquad];
1991 template <
int dim,
int n_faces,
typename real>
1993 const int nstate_input,
1994 const unsigned int max_degree_input,
1995 const unsigned int grid_degree_input)
2002 template <
int dim,
int n_faces,
typename real>
2004 const dealii::FESystem<1,1> &finite_element,
2005 const dealii::Quadrature<0> &face_quadrature)
2007 const unsigned int n_face_quad_pts = face_quadrature.size();
2008 const unsigned int n_dofs = finite_element.dofs_per_cell;
2009 const unsigned int n_faces_1D = n_faces / dim;
2010 const std::vector<double> &quad_weights = face_quadrature.get_weights ();
2012 for(
unsigned int iface=0; iface<n_faces_1D; iface++){
2016 const dealii::Quadrature<1> quadrature = dealii::QProjector<1>::project_to_face(dealii::ReferenceCell::get_hypercube(1),
2019 for(
unsigned int iquad=0; iquad<n_face_quad_pts; iquad++){
2020 const dealii::Point<1> qpoint = quadrature.point(iquad);
2021 for(
unsigned int idof=0; idof<n_dofs; idof++){
2022 const int istate = finite_element.system_to_component_index(idof).first;
2024 this->
oneD_surf_operator[iface][iquad][idof] = finite_element.shape_value_component(idof,qpoint,istate)
2025 * quad_weights[iquad];
2031 template <
int dim,
int n_faces,
typename real>
2033 const int nstate_input,
2034 const unsigned int max_degree_input,
2035 const unsigned int grid_degree_input)
2042 template <
int dim,
int n_faces,
typename real>
2044 const dealii::FESystem<1,1> &finite_element,
2045 const dealii::Quadrature<1> &quadrature)
2047 const unsigned int n_dofs = finite_element.dofs_per_cell;
2055 template <
int dim,
int n_faces,
typename real>
2057 const unsigned int n_dofs,
2058 const dealii::FullMatrix<double> &norm_matrix,
2059 const dealii::FullMatrix<double> &face_integral,
2060 dealii::FullMatrix<double> &lifting)
2062 dealii::FullMatrix<double> norm_inv(n_dofs);
2063 norm_inv.invert(norm_matrix);
2064 norm_inv.mTmult(lifting, face_integral);
2066 template <
int dim,
int n_faces,
typename real>
2068 const dealii::FESystem<1,1> &finite_element,
2069 const dealii::Quadrature<0> &face_quadrature)
2071 const unsigned int n_face_quad_pts = face_quadrature.size();
2072 const unsigned int n_dofs = finite_element.dofs_per_cell;
2073 const unsigned int n_faces_1D = n_faces / dim;
2078 for(
unsigned int iface=0; iface<n_faces_1D; iface++){
2085 template <
int dim,
int n_faces,
typename real>
2087 const int nstate_input,
2088 const unsigned int max_degree_input,
2089 const unsigned int grid_degree_input,
2092 , FR_param_type(FR_param_input)
2098 template <
int dim,
int n_faces,
typename real>
2100 const dealii::FESystem<1,1> &finite_element,
2101 const dealii::Quadrature<1> &quadrature)
2103 const unsigned int n_dofs = finite_element.dofs_per_cell;
2114 template <
int dim,
int n_faces,
typename real>
2116 const dealii::FESystem<1,1> &finite_element,
2117 const dealii::Quadrature<0> &face_quadrature)
2119 const unsigned int n_face_quad_pts = face_quadrature.size();
2120 const unsigned int n_dofs = finite_element.dofs_per_cell;
2121 const unsigned int n_faces_1D = n_faces / dim;
2126 for(
unsigned int iface=0; iface<n_faces_1D; iface++){
2140 template <
int dim,
int n_faces,
typename real>
2142 const int nstate_input,
2143 const unsigned int max_degree_input,
2144 const unsigned int grid_degree_input)
2146 , mapping_shape_functions_grid_nodes(nstate_input, max_degree_input, grid_degree_input)
2147 , mapping_shape_functions_flux_nodes(nstate_input, max_degree_input, grid_degree_input)
2154 template <
int dim,
int n_faces,
typename real>
2156 const dealii::FESystem<1,1> &finite_element,
2157 const dealii::Quadrature<1> &quadrature)
2159 assert(finite_element.dofs_per_cell == quadrature.size());
2163 template <
int dim,
int n_faces,
typename real>
2165 const dealii::FESystem<1,1> &finite_element,
2166 const dealii::Quadrature<1> &quadrature,
2167 const dealii::Quadrature<0> &face_quadrature)
2175 template <
int dim,
int n_faces,
typename real>
2177 const dealii::FESystem<1,1> &finite_element,
2178 const dealii::Quadrature<1> &quadrature)
2190 template <
typename real,
int dim,
int n_faces>
2192 const int nstate_input,
2193 const unsigned int max_degree_input,
2194 const unsigned int grid_degree_input,
2195 const bool store_vol_flux_nodes_input,
2196 const bool store_surf_flux_nodes_input,
2197 const bool store_Jacobian_input)
2199 , store_Jacobian(store_Jacobian_input)
2200 , store_vol_flux_nodes(store_vol_flux_nodes_input)
2201 , store_surf_flux_nodes(store_surf_flux_nodes_input)
2204 template <
typename real,
int dim,
int n_faces>
2206 const dealii::Tensor<1,dim,real> &phys,
2207 const dealii::Tensor<2,dim,real> &metric_cofactor,
2208 dealii::Tensor<1,dim,real> &ref)
2210 for(
int idim=0; idim<dim; idim++){
2211 for(
int idim2=0; idim2<dim; idim2++){
2212 ref[idim] += metric_cofactor[idim2][idim] * phys[idim2];
2218 template <
typename real,
int dim,
int n_faces>
2220 const dealii::Tensor<1,dim,real> &ref,
2221 const dealii::Tensor<2,dim,real> &metric_cofactor,
2222 dealii::Tensor<1,dim,real> &phys)
2224 for(
int idim=0; idim<dim; idim++){
2225 phys[idim] = metric_cofactor[idim] * ref;
2234 template <
typename real,
int dim,
int n_faces>
2236 const dealii::Tensor<1,dim,std::vector<real>> &phys,
2237 const dealii::Tensor<2,dim,std::vector<real>> &metric_cofactor,
2238 dealii::Tensor<1,dim,std::vector<real>> &ref)
2240 assert(phys[0].size() == metric_cofactor[0][0].size());
2241 const unsigned int n_pts = phys[0].size();
2242 for(
int idim=0; idim<dim; idim++){
2243 ref[idim].resize(n_pts);
2244 for(
int idim2=0; idim2<dim; idim2++){
2245 for(
unsigned int ipt=0; ipt<n_pts; ipt++){
2246 ref[idim][ipt] += metric_cofactor[idim2][idim][ipt] * phys[idim2][ipt];
2253 template <
typename real,
int dim,
int n_faces>
2255 const unsigned int n_quad_pts,
2256 const dealii::Tensor<1,dim,real> &ref,
2257 const dealii::Tensor<2,dim,std::vector<real>> &metric_cofactor,
2258 std::vector<dealii::Tensor<1,dim,real>> &phys)
2260 for(
unsigned int iquad=0; iquad<n_quad_pts; iquad++){
2261 for(
int idim=0; idim<dim; idim++){
2262 phys[iquad][idim] = 0.0;
2263 for(
int idim2=0; idim2<dim; idim2++){
2264 phys[iquad][idim] += metric_cofactor[idim][idim2][iquad]
2268 phys[iquad] /= phys[iquad].norm();
2272 template <
typename real,
int dim,
int n_faces>
2274 const unsigned int n_quad_pts,
2275 const unsigned int ,
2276 const std::array<std::vector<real>,dim> &mapping_support_points,
2283 mapping_support_points,
2293 template <
typename real,
int dim,
int n_faces>
2295 const unsigned int n_quad_pts,
2296 const unsigned int n_metric_dofs,
2297 const std::array<std::vector<real>,dim> &mapping_support_points,
2299 const bool use_invariant_curl_form)
2302 for(
int idim=0; idim<dim; idim++){
2303 for(
int jdim=0; jdim<dim; jdim++){
2310 mapping_support_points,
2322 mapping_support_points,
2336 use_invariant_curl_form);
2339 for(
int idim=0; idim<dim; idim++){
2350 template <
typename real,
int dim,
int n_faces>
2352 const unsigned int iface,
2353 const unsigned int n_quad_pts,
2354 const unsigned int n_metric_dofs,
2355 const std::array<std::vector<real>,dim> &mapping_support_points,
2357 const bool use_invariant_curl_form)
2360 for(
int idim=0; idim<dim; idim++){
2361 for(
int jdim=0; jdim<dim; jdim++){
2368 mapping_support_points,
2392 mapping_support_points,
2418 use_invariant_curl_form);
2421 for(
int iface=0; iface<n_faces; iface++){
2422 for(
int idim=0; idim<dim; idim++){
2440 template <
typename real,
int dim,
int n_faces>
2442 const unsigned int n_quad_pts,
2443 const std::array<std::vector<real>,dim> &mapping_support_points,
2444 const dealii::FullMatrix<double> &basis_x_flux_nodes,
2445 const dealii::FullMatrix<double> &basis_y_flux_nodes,
2446 const dealii::FullMatrix<double> &basis_z_flux_nodes,
2447 const dealii::FullMatrix<double> &grad_basis_x_flux_nodes,
2448 const dealii::FullMatrix<double> &grad_basis_y_flux_nodes,
2449 const dealii::FullMatrix<double> &grad_basis_z_flux_nodes,
2450 std::vector<dealii::Tensor<2,dim,real>> &local_Jac)
2452 for(
int idim=0; idim<dim; idim++){
2453 for(
int jdim=0; jdim<dim; jdim++){
2455 std::vector<real> output_vect(n_quad_pts);
2458 grad_basis_x_flux_nodes,
2460 basis_z_flux_nodes);
2464 grad_basis_y_flux_nodes,
2465 basis_z_flux_nodes);
2470 grad_basis_z_flux_nodes);
2471 for(
unsigned int iquad=0; iquad<n_quad_pts; iquad++){
2472 local_Jac[iquad][idim][jdim] = output_vect[iquad];
2478 template <
typename real,
int dim,
int n_faces>
2480 const unsigned int n_quad_pts,
2481 const std::array<std::vector<real>,dim> &mapping_support_points,
2482 const dealii::FullMatrix<double> &basis_x_flux_nodes,
2483 const dealii::FullMatrix<double> &basis_y_flux_nodes,
2484 const dealii::FullMatrix<double> &basis_z_flux_nodes,
2485 const dealii::FullMatrix<double> &grad_basis_x_flux_nodes,
2486 const dealii::FullMatrix<double> &grad_basis_y_flux_nodes,
2487 const dealii::FullMatrix<double> &grad_basis_z_flux_nodes,
2488 std::vector<real> &det_metric_Jac)
2491 assert(pow(this->
max_grid_degree+1,dim) == mapping_support_points[0].size());
2493 std::vector<dealii::Tensor<2,dim,double>> Jacobian_flux_nodes(n_quad_pts);
2495 mapping_support_points,
2499 grad_basis_x_flux_nodes,
2500 grad_basis_y_flux_nodes,
2501 grad_basis_z_flux_nodes,
2502 Jacobian_flux_nodes);
2504 for(
int idim=0; idim<dim; idim++){
2505 for(
int jdim=0; jdim<dim; jdim++){
2507 for(
unsigned int iquad=0; iquad<n_quad_pts; iquad++){
2514 for(
unsigned int iquad=0; iquad<n_quad_pts; iquad++){
2515 det_metric_Jac[iquad] = dealii::determinant(Jacobian_flux_nodes[iquad]);
2517 if(det_metric_Jac[iquad] <= 1e-14){
2518 std::cout<<
"The determinant of the Jacobian is negative. Aborting..."<<std::endl;
2523 template <
typename real,
int dim,
int n_faces>
2525 const unsigned int n_quad_pts,
2526 const unsigned int n_metric_dofs,
2527 const std::array<std::vector<real>,dim> &mapping_support_points,
2528 const dealii::FullMatrix<double> &basis_x_grid_nodes,
2529 const dealii::FullMatrix<double> &basis_y_grid_nodes,
2530 const dealii::FullMatrix<double> &basis_z_grid_nodes,
2531 const dealii::FullMatrix<double> &basis_x_flux_nodes,
2532 const dealii::FullMatrix<double> &basis_y_flux_nodes,
2533 const dealii::FullMatrix<double> &basis_z_flux_nodes,
2534 const dealii::FullMatrix<double> &grad_basis_x_grid_nodes,
2535 const dealii::FullMatrix<double> &grad_basis_y_grid_nodes,
2536 const dealii::FullMatrix<double> &grad_basis_z_grid_nodes,
2537 const dealii::FullMatrix<double> &grad_basis_x_flux_nodes,
2538 const dealii::FullMatrix<double> &grad_basis_y_flux_nodes,
2539 const dealii::FullMatrix<double> &grad_basis_z_flux_nodes,
2540 dealii::Tensor<2,dim,std::vector<real>> &metric_cofactor,
2541 const bool use_invariant_curl_form)
2546 std::fill(metric_cofactor[0][0].begin(), metric_cofactor[0][0].end(), 1.0);
2549 std::vector<dealii::Tensor<2,dim,double>> Jacobian_flux_nodes(n_quad_pts);
2551 mapping_support_points,
2555 grad_basis_x_flux_nodes,
2556 grad_basis_y_flux_nodes,
2557 grad_basis_z_flux_nodes,
2558 Jacobian_flux_nodes);
2559 for(
unsigned int iquad=0; iquad<n_quad_pts; iquad++){
2560 metric_cofactor[0][0][iquad] = Jacobian_flux_nodes[iquad][1][1];
2561 metric_cofactor[1][0][iquad] = - Jacobian_flux_nodes[iquad][0][1];
2562 metric_cofactor[0][1][iquad] = - Jacobian_flux_nodes[iquad][1][0];
2563 metric_cofactor[1][1][iquad] = Jacobian_flux_nodes[iquad][0][0];
2569 mapping_support_points,
2576 grad_basis_x_grid_nodes,
2577 grad_basis_y_grid_nodes,
2578 grad_basis_z_grid_nodes,
2579 grad_basis_x_flux_nodes,
2580 grad_basis_y_flux_nodes,
2581 grad_basis_z_flux_nodes,
2583 use_invariant_curl_form);
2587 template <
typename real,
int dim,
int n_faces>
2589 const unsigned int n_metric_dofs,
2590 const unsigned int ,
2591 const std::array<std::vector<real>,dim> &mapping_support_points,
2592 const dealii::FullMatrix<double> &basis_x_grid_nodes,
2593 const dealii::FullMatrix<double> &basis_y_grid_nodes,
2594 const dealii::FullMatrix<double> &basis_z_grid_nodes,
2595 const dealii::FullMatrix<double> &basis_x_flux_nodes,
2596 const dealii::FullMatrix<double> &basis_y_flux_nodes,
2597 const dealii::FullMatrix<double> &basis_z_flux_nodes,
2598 const dealii::FullMatrix<double> &grad_basis_x_grid_nodes,
2599 const dealii::FullMatrix<double> &grad_basis_y_grid_nodes,
2600 const dealii::FullMatrix<double> &grad_basis_z_grid_nodes,
2601 const dealii::FullMatrix<double> &grad_basis_x_flux_nodes,
2602 const dealii::FullMatrix<double> &grad_basis_y_flux_nodes,
2603 const dealii::FullMatrix<double> &grad_basis_z_flux_nodes,
2604 dealii::Tensor<2,dim,std::vector<real>> &metric_cofactor,
2605 const bool use_invariant_curl_form)
2610 std::vector<dealii::Tensor<2,dim,real>> grad_Xm(n_metric_dofs);
2612 mapping_support_points,
2616 grad_basis_x_grid_nodes,
2617 grad_basis_y_grid_nodes,
2618 grad_basis_z_grid_nodes,
2623 std::vector<real> z_dy_dxi(n_metric_dofs);
2624 std::vector<real> z_dy_deta(n_metric_dofs);
2625 std::vector<real> z_dy_dzeta(n_metric_dofs);
2627 std::vector<real> x_dz_dxi(n_metric_dofs);
2628 std::vector<real> x_dz_deta(n_metric_dofs);
2629 std::vector<real> x_dz_dzeta(n_metric_dofs);
2631 std::vector<real> y_dx_dxi(n_metric_dofs);
2632 std::vector<real> y_dx_deta(n_metric_dofs);
2633 std::vector<real> y_dx_dzeta(n_metric_dofs);
2635 for(
unsigned int grid_node=0; grid_node<n_metric_dofs; grid_node++){
2636 z_dy_dxi[grid_node] = mapping_support_points[2][grid_node] * grad_Xm[grid_node][1][0];
2637 if(use_invariant_curl_form){
2638 z_dy_dxi[grid_node] = 0.5 * z_dy_dxi[grid_node]
2639 - 0.5 * mapping_support_points[1][grid_node] * grad_Xm[grid_node][2][0];
2641 z_dy_deta[grid_node] = mapping_support_points[2][grid_node] * grad_Xm[grid_node][1][1];
2642 if(use_invariant_curl_form){
2643 z_dy_deta[grid_node] = 0.5 * z_dy_deta[grid_node]
2644 - 0.5 * mapping_support_points[1][grid_node] * grad_Xm[grid_node][2][1];
2646 z_dy_dzeta[grid_node] = mapping_support_points[2][grid_node] * grad_Xm[grid_node][1][2];
2647 if(use_invariant_curl_form){
2648 z_dy_dzeta[grid_node] = 0.5 * z_dy_dzeta[grid_node]
2649 - 0.5 * mapping_support_points[1][grid_node] * grad_Xm[grid_node][2][2];
2652 x_dz_dxi[grid_node] = mapping_support_points[0][grid_node] * grad_Xm[grid_node][2][0];
2653 if(use_invariant_curl_form){
2654 x_dz_dxi[grid_node] = 0.5 * x_dz_dxi[grid_node]
2655 - 0.5 * mapping_support_points[2][grid_node] * grad_Xm[grid_node][1][0];
2657 x_dz_deta[grid_node] = mapping_support_points[0][grid_node] * grad_Xm[grid_node][2][1];
2658 if(use_invariant_curl_form){
2659 x_dz_deta[grid_node] = 0.5 * x_dz_deta[grid_node]
2660 - 0.5 * mapping_support_points[2][grid_node] * grad_Xm[grid_node][1][1];
2662 x_dz_dzeta[grid_node] = mapping_support_points[0][grid_node] * grad_Xm[grid_node][2][2];
2663 if(use_invariant_curl_form){
2664 x_dz_dzeta[grid_node] = 0.5 * x_dz_dzeta[grid_node]
2665 - 0.5 * mapping_support_points[2][grid_node] * grad_Xm[grid_node][1][2];
2668 y_dx_dxi[grid_node] = mapping_support_points[1][grid_node] * grad_Xm[grid_node][0][0];
2669 if(use_invariant_curl_form){
2670 y_dx_dxi[grid_node] = 0.5 * y_dx_dxi[grid_node]
2671 - 0.5 * mapping_support_points[0][grid_node] * grad_Xm[grid_node][1][0];
2673 y_dx_deta[grid_node] = mapping_support_points[1][grid_node] * grad_Xm[grid_node][0][1];
2674 if(use_invariant_curl_form){
2675 y_dx_deta[grid_node] = 0.5 * y_dx_deta[grid_node]
2676 - 0.5 * mapping_support_points[0][grid_node] * grad_Xm[grid_node][1][1];
2678 y_dx_dzeta[grid_node] = mapping_support_points[1][grid_node] * grad_Xm[grid_node][0][2];
2679 if(use_invariant_curl_form){
2680 y_dx_dzeta[grid_node] = 0.5 * y_dx_dzeta[grid_node]
2681 - 0.5 * mapping_support_points[0][grid_node] * grad_Xm[grid_node][1][2];
2688 basis_x_flux_nodes, grad_basis_y_flux_nodes, basis_z_flux_nodes,
false, -1.0);
2690 basis_x_flux_nodes, basis_y_flux_nodes, grad_basis_z_flux_nodes,
true);
2693 grad_basis_x_flux_nodes, basis_y_flux_nodes, basis_z_flux_nodes);
2695 basis_x_flux_nodes, basis_y_flux_nodes, grad_basis_z_flux_nodes,
true, -1.0);
2698 grad_basis_x_flux_nodes, basis_y_flux_nodes, basis_z_flux_nodes,
false, -1.0);
2700 basis_x_flux_nodes, grad_basis_y_flux_nodes, basis_z_flux_nodes,
true);
2704 basis_x_flux_nodes, grad_basis_y_flux_nodes, basis_z_flux_nodes,
false, -1.0);
2706 basis_x_flux_nodes, basis_y_flux_nodes, grad_basis_z_flux_nodes,
true);
2709 grad_basis_x_flux_nodes, basis_y_flux_nodes, basis_z_flux_nodes);
2711 basis_x_flux_nodes, basis_y_flux_nodes, grad_basis_z_flux_nodes,
true, -1.0);
2714 grad_basis_x_flux_nodes, basis_y_flux_nodes, basis_z_flux_nodes,
false, -1.0);
2716 basis_x_flux_nodes, grad_basis_y_flux_nodes, basis_z_flux_nodes,
true);
2720 basis_x_flux_nodes, grad_basis_y_flux_nodes, basis_z_flux_nodes,
false, -1.0);
2722 basis_x_flux_nodes, basis_y_flux_nodes, grad_basis_z_flux_nodes,
true);
2725 grad_basis_x_flux_nodes, basis_y_flux_nodes, basis_z_flux_nodes);
2727 basis_x_flux_nodes, basis_y_flux_nodes, grad_basis_z_flux_nodes,
true, -1.0);
2730 grad_basis_x_flux_nodes, basis_y_flux_nodes, basis_z_flux_nodes,
false, -1.0);
2732 basis_x_flux_nodes, grad_basis_y_flux_nodes, basis_z_flux_nodes,
true);
2742 template <
int dim,
int nstate,
int n_faces,
typename real>
2744 const unsigned int max_degree_input,
2745 const unsigned int grid_degree_input)
2749 template <
int dim,
int nstate,
int n_faces,
typename real>
2751 const unsigned int max_degree_input,
2752 const unsigned int grid_degree_input)
2759 template <
int dim,
int nstate,
int n_faces,
typename real>
2761 const dealii::FESystem<1,1> &finite_element,
2762 const dealii::Quadrature<1> &quadrature)
2764 const unsigned int n_quad_pts = quadrature.size();
2765 const unsigned int n_dofs = finite_element.dofs_per_cell;
2766 const unsigned int n_shape_fns = n_dofs /
nstate;
2769 for(
unsigned int iquad=0; iquad<n_quad_pts; iquad++){
2770 const dealii::Point<1> qpoint = quadrature.point(iquad);
2771 for(
unsigned int idof=0; idof<n_dofs; idof++){
2772 const unsigned int istate = finite_element.system_to_component_index(idof).first;
2773 const unsigned int ishape = finite_element.system_to_component_index(idof).second;
2778 this->
oneD_vol_state_operator[istate][iquad][ishape] = finite_element.shape_value_component(idof,qpoint,istate);
2783 template <
int dim,
int nstate,
int n_faces,
typename real>
2785 const dealii::FESystem<1,1> &finite_element,
2786 const dealii::Quadrature<1> &quadrature)
2788 const unsigned int n_quad_pts = quadrature.size();
2789 const unsigned int n_dofs = finite_element.dofs_per_cell;
2790 const unsigned int n_shape_fns = n_dofs /
nstate;
2792 for(
unsigned int iquad=0; iquad<n_quad_pts; iquad++){
2793 const dealii::Point<1> qpoint = quadrature.point(iquad);
2794 for(
unsigned int idof=0; idof<n_dofs; idof++){
2795 const unsigned int istate = finite_element.system_to_component_index(idof).first;
2796 const unsigned int ishape = finite_element.system_to_component_index(idof).second;
2800 this->
oneD_grad_state_operator[istate][iquad][ishape] = finite_element.shape_grad_component(idof, qpoint, istate)[0];
2804 template <
int dim,
int nstate,
int n_faces,
typename real>
2806 const dealii::FESystem<1,1> &finite_element,
2807 const dealii::Quadrature<0> &face_quadrature)
2809 const unsigned int n_face_quad_pts = face_quadrature.size();
2810 const unsigned int n_dofs = finite_element.dofs_per_cell;
2811 const unsigned int n_faces_1D = n_faces / dim;
2812 const unsigned int n_shape_fns = n_dofs /
nstate;
2814 for(
unsigned int iface=0; iface<n_faces_1D; iface++){
2815 const dealii::Quadrature<1> quadrature = dealii::QProjector<1>::project_to_face(dealii::ReferenceCell::get_hypercube(1),
2818 for(
unsigned int iquad=0; iquad<n_face_quad_pts; iquad++){
2819 for(
unsigned int idof=0; idof<n_dofs; idof++){
2820 const unsigned int istate = finite_element.system_to_component_index(idof).first;
2821 const unsigned int ishape = finite_element.system_to_component_index(idof).second;
2825 this->
oneD_surf_state_operator[istate][iface][iquad][ishape] = finite_element.shape_value_component(idof,quadrature.point(iquad),istate);
2831 template <
int dim,
int nstate,
int n_faces,
typename real>
2833 const unsigned int max_degree_input,
2834 const unsigned int grid_degree_input)
2841 template <
int dim,
int nstate,
int n_faces,
typename real>
2843 const dealii::FESystem<1,1> &finite_element,
2844 const dealii::Quadrature<1> &quadrature)
2846 const unsigned int n_quad_pts = quadrature.size();
2847 const unsigned int n_dofs = finite_element.dofs_per_cell;
2848 assert(n_dofs == n_quad_pts);
2851 for(
int istate=0; istate<
nstate; istate++){
2854 for(
unsigned int iquad=0; iquad<n_quad_pts; iquad++){
2855 const dealii::Point<1> qpoint = quadrature.point(iquad);
2856 for(
unsigned int idof=0; idof<n_dofs; idof++){
2864 template <
int dim,
int nstate,
int n_faces,
typename real>
2866 const dealii::FESystem<1,1> &finite_element,
2867 const dealii::Quadrature<1> &quadrature)
2869 const unsigned int n_quad_pts = quadrature.size();
2870 const unsigned int n_dofs = finite_element.dofs_per_cell;
2871 assert(n_dofs == n_quad_pts);
2873 for(
int istate=0; istate<
nstate; istate++){
2876 for(
unsigned int iquad=0; iquad<n_quad_pts; iquad++){
2877 const dealii::Point<1> qpoint = quadrature.point(iquad);
2878 for(
unsigned int idof=0; idof<n_dofs; idof++){
2884 template <
int dim,
int nstate,
int n_faces,
typename real>
2886 const dealii::FESystem<1,1> &finite_element,
2887 const dealii::Quadrature<0> &face_quadrature)
2889 const unsigned int n_face_quad_pts = face_quadrature.size();
2890 const unsigned int n_dofs = finite_element.dofs_per_cell;
2891 const unsigned int n_faces_1D = n_faces / dim;
2893 for(
unsigned int iface=0; iface<n_faces_1D; iface++){
2894 const dealii::Quadrature<1> quadrature = dealii::QProjector<1>::project_to_face(dealii::ReferenceCell::get_hypercube(1),
2897 for(
int istate=0; istate<
nstate; istate++){
2899 for(
unsigned int iquad=0; iquad<n_face_quad_pts; iquad++){
2900 for(
unsigned int idof=0; idof<n_dofs; idof++){
2901 this->
oneD_surf_state_operator[istate][iface][iquad][idof] = finite_element.shape_value_component(idof,quadrature.point(iquad),0);
2909 template <
int dim,
int nstate,
int n_faces,
typename real>
2911 const unsigned int max_degree_input,
2912 const unsigned int grid_degree_input)
2919 template <
int dim,
int nstate,
int n_faces,
typename real>
2921 const dealii::FESystem<1,1> &finite_element,
2922 const dealii::Quadrature<1> &quadrature)
2924 const unsigned int n_quad_pts = quadrature.size();
2925 const unsigned int n_dofs_flux = quadrature.size();
2926 const unsigned int n_dofs = finite_element.dofs_per_cell;
2928 const std::vector<double> &quad_weights = quadrature.get_weights ();
2929 for(
int istate_flux=0; istate_flux<
nstate; istate_flux++){
2932 for(
unsigned int itest=0; itest<n_dofs; itest++){
2933 const int istate_test = finite_element.system_to_component_index(itest).first;
2934 for(
unsigned int idof=0; idof<n_dofs_flux; idof++){
2936 for(
unsigned int iquad=0; iquad<n_quad_pts; iquad++){
2937 const dealii::Point<1> qpoint = quadrature.point(iquad);
2938 value += finite_element.shape_value_component(itest, qpoint, istate_test)
2939 * quad_weights[iquad]
The FLUX basis functions separated by nstate with n shape functions.
unsigned int current_degree
Stores the degree of the current poly degree.
void surface_two_pt_flux_Hadamard_product(const dealii::FullMatrix< real > &input_mat, std::vector< real > &output_vect_vol, std::vector< real > &output_vect_surf, const std::vector< real > &weights, const std::array< dealii::FullMatrix< double >, 2 > &surf_basis, const unsigned int iface, const unsigned int dim_not_zero, const double scaling=2.0)
Computes the surface cross Hadamard products for skew-symmetric form from Eq. (15) in Chan...
unsigned int current_degree
Stores the degree of the current poly degree.
dealii::FullMatrix< double > build_dim_mass_matrix(const int nstate, const unsigned int n_dofs, const unsigned int n_quad_pts, basis_functions< dim, n_faces, real > &basis, const std::vector< double > &det_Jac, const std::vector< double > &quad_weights)
Assemble the dim mass matrix on the fly with metric Jacobian dependence.
basis_functions< dim, n_faces, real > mapping_shape_functions_grid_nodes
Object of mapping shape functions evaluated at grid nodes.
dealii::Tensor< 2, dim, std::vector< real > > metric_cofactor_vol
The volume metric cofactor matrix.
void build_determinant_volume_metric_Jacobian(const unsigned int n_quad_pts, const unsigned int n_metric_dofs, const std::array< std::vector< real >, dim > &mapping_support_points, mapping_shape_functions< dim, n_faces, real > &mapping_basis)
Builds just the determinant of the volume metric determinant.
void build_1D_volume_state_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
void build_local_surface_lifting_operator(const unsigned int n_dofs, const dealii::FullMatrix< double > &norm_matrix, const dealii::FullMatrix< double > &face_integral, dealii::FullMatrix< double > &lifting)
Builds the local lifting operator.
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
const Parameters::AllParameters::Flux_Reconstruction FR_param_type
Flux reconstruction parameter type.
std::array< dealii::FullMatrix< double >, 2 > oneD_surf_operator
Stores the one dimensional surface operator.
The metric independent inverse of the FR mass matrix .
void divergence_two_pt_flux_Hadamard_product(const dealii::Tensor< 1, dim, dealii::FullMatrix< real >> &input_mat, std::vector< real > &output_vect, const std::vector< real > &weights, const dealii::FullMatrix< double > &basis, const double scaling=2.0)
Computes the divergence of the 2pt flux Hadamard products, then sums the rows.
bool store_transpose
Flag is store transpose operator.
unsigned int current_degree
Stores the degree of the current poly degree.
metric_operators(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input, const bool store_vol_flux_nodes_input=false, const bool store_surf_flux_nodes_input=false, const bool store_Jacobian_input=false)
Constructor.
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
void get_FR_aux_correction_parameter(const unsigned int curr_cell_degree, double &k)
Gets the FR correction parameter for the auxiliary equations and stores.
void sum_factorized_Hadamard_surface_sparsity_pattern(const unsigned int rows_size, const unsigned int columns_size, std::vector< unsigned int > &rows, std::vector< unsigned int > &columns, const int dim_not_zero)
Computes the rows and columns vectors with non-zero indices for surface sum-factorized Hadamard produ...
The integration of gradient of solution basis.
unsigned int current_degree
Stores the degree of the current poly degree.
unsigned int current_degree
Stores the degree of the current poly degree.
void build_1D_gradient_state_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
Sum Factorization derived class.
lifting_operator_FR(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input, const Parameters::AllParameters::Flux_Reconstruction FR_param_input)
Constructor.
void matrix_vector_mult(const std::vector< real > &input_vect, std::vector< real > &output_vect, const dealii::FullMatrix< double > &basis_x, const dealii::FullMatrix< double > &basis_y, const dealii::FullMatrix< double > &basis_z, const bool adding=false, const double factor=1.0) override
Computes a matrix-vector product using sum-factorization. Pass the one-dimensional basis...
unsigned int current_grid_degree
Stores the degree of the current grid degree.
const unsigned int max_grid_degree
Max grid degree.
unsigned int current_degree
Stores the degree of the current poly degree.
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
void inner_product_surface_1D(const unsigned int face_number, const std::vector< real > &input_vect, const std::vector< real > &weight_vect, std::vector< real > &output_vect, const std::array< dealii::FullMatrix< double >, 2 > &basis_surf, const dealii::FullMatrix< double > &basis_vol, const bool adding=false, const double factor=1.0)
Apply sum-factorization inner product on a surface.
const Parameters::AllParameters::Flux_Reconstruction_Aux FR_param_type
Flux reconstruction parameter type.
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
dealii::FullMatrix< double > oneD_grad_operator
Stores the one dimensional gradient operator.
dealii::ConditionalOStream pcout
Parallel std::cout that only outputs on mpi_rank==0.
local_Flux_Reconstruction_operator(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input, const Parameters::AllParameters::Flux_Reconstruction FR_param_input)
Constructor.
dealii::FullMatrix< double > tensor_product(const dealii::FullMatrix< double > &basis_x, const dealii::FullMatrix< double > &basis_y, const dealii::FullMatrix< double > &basis_z)
Returns the tensor product of matrices passed.
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
The metric independent FR mass matrix for auxiliary equation .
FR_mass_aux(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input, const Parameters::AllParameters::Flux_Reconstruction_Aux FR_param_input)
Constructor.
dealii::FullMatrix< double > tensor_product_state(const int nstate, const dealii::FullMatrix< double > &basis_x, const dealii::FullMatrix< double > &basis_y, const dealii::FullMatrix< double > &basis_z)
Returns the tensor product of matrices passed, but makes it sparse diagonal by state.
bool store_transpose
Flag is store transpose operator.
const Parameters::AllParameters::Flux_Reconstruction_Aux FR_param_type
Flux reconstruction parameter type.
const int nstate
Number of states.
void transform_reference_unit_normal_to_physical_unit_normal(const unsigned int n_quad_pts, const dealii::Tensor< 1, dim, real > &ref, const dealii::Tensor< 2, dim, std::vector< real >> &metric_cofactor, std::vector< dealii::Tensor< 1, dim, real >> &phys)
Given a reference tensor, return the physical tensor.
void Hadamard_product(const dealii::FullMatrix< real > &input_mat1, const dealii::FullMatrix< real > &input_mat2, dealii::FullMatrix< real > &output_mat)
Computes a single Hadamard product.
double compute_factorial(double n)
Standard function to compute factorial of a number.
Projection operator corresponding to basis functions onto -norm for auxiliary equation.
void inner_product_1D(const std::vector< real > &input_vect, const std::vector< real > &weight_vect, std::vector< real > &output_vect, const dealii::FullMatrix< double > &basis_x, const bool adding=false, const double factor=1.0)
Apply the inner product operation using the 1D operator in each direction.
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
unsigned int current_degree
Stores the degree of the current poly degree.
Files for the baseline physics.
void build_1D_surface_state_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 0 > &face_quadrature)
Assembles the one dimensional operator.
double FR_param_aux
Flux reconstruction paramater value.
lifting_operator(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input)
Constructor.
unsigned int current_degree
Stores the degree of the current poly degree.
void build_1D_surface_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 0 > &face_quadrature)
Assembles the one dimensional operator.
const unsigned int max_degree
Max polynomial degree.
void compute_local_3D_cofactor(const unsigned int n_metric_dofs, const unsigned int n_quad_pts, const std::array< std::vector< real >, dim > &mapping_support_points, const dealii::FullMatrix< double > &basis_x_grid_nodes, const dealii::FullMatrix< double > &basis_y_grid_nodes, const dealii::FullMatrix< double > &basis_z_grid_nodes, const dealii::FullMatrix< double > &basis_x_flux_nodes, const dealii::FullMatrix< double > &basis_y_flux_nodes, const dealii::FullMatrix< double > &basis_z_flux_nodes, const dealii::FullMatrix< double > &grad_basis_x_grid_nodes, const dealii::FullMatrix< double > &grad_basis_y_grid_nodes, const dealii::FullMatrix< double > &grad_basis_z_grid_nodes, const dealii::FullMatrix< double > &grad_basis_x_flux_nodes, const dealii::FullMatrix< double > &grad_basis_y_flux_nodes, const dealii::FullMatrix< double > &grad_basis_z_flux_nodes, dealii::Tensor< 2, dim, std::vector< real >> &metric_cofactor, const bool use_invariant_curl_form=false)
Computes local 3D cofactor matrix.
unsigned int current_degree
Stores the degree of the current poly degree.
std::array< dealii::FullMatrix< double >, 2 > oneD_surf_grad_operator
Stores the one dimensional surface gradient operator.
dealii::Tensor< 2, dim, std::vector< real > > metric_cofactor_surf
The facet metric cofactor matrix, for ONE face.
dealii::FullMatrix< double > oneD_transpose_vol_operator
Stores the transpose of the operator for fast weight-adjusted solves.
unsigned int current_degree
Stores the degree of the current poly degree.
unsigned int current_degree
Stores the degree of the current poly degree.
basis_functions_state(const unsigned int max_degree_input, const unsigned int grid_degree_input)
Constructor.
Flux_Reconstruction
Type of correction in Flux Reconstruction.
dealii::FullMatrix< double > oneD_skew_symm_vol_oper
Skew-symmetric volume operator .
const bool store_Jacobian
Flag if store metric Jacobian at flux nodes.
void get_Huynh_g2_parameter(const unsigned int curr_cell_degree, double &c)
Evaluates Huynh's g2 parameter for flux reconstruction.
-th order modal derivative of basis fuctions, ie/
face_integral_basis(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input)
Constructor.
void build_1D_shape_functions_at_grid_nodes(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Constructs the volume operator and gradient operator.
That is Quadrature Weights multiplies with basis_at_vol_cubature.
local_Flux_Reconstruction_operator_aux(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input, const Parameters::AllParameters::Flux_Reconstruction_Aux FR_param_aux_input)
Constructor.
const bool store_surf_flux_nodes
Flag if store metric Jacobian at flux nodes.
void build_local_metric_cofactor_matrix(const unsigned int n_quad_pts, const unsigned int n_metric_dofs, const std::array< std::vector< real >, dim > &mapping_support_points, const dealii::FullMatrix< double > &basis_x_grid_nodes, const dealii::FullMatrix< double > &basis_y_grid_nodes, const dealii::FullMatrix< double > &basis_z_grid_nodes, const dealii::FullMatrix< double > &basis_x_flux_nodes, const dealii::FullMatrix< double > &basis_y_flux_nodes, const dealii::FullMatrix< double > &basis_z_flux_nodes, const dealii::FullMatrix< double > &grad_basis_x_grid_nodes, const dealii::FullMatrix< double > &grad_basis_y_grid_nodes, const dealii::FullMatrix< double > &grad_basis_z_grid_nodes, const dealii::FullMatrix< double > &grad_basis_x_flux_nodes, const dealii::FullMatrix< double > &grad_basis_y_flux_nodes, const dealii::FullMatrix< double > &grad_basis_z_flux_nodes, dealii::Tensor< 2, dim, std::vector< real >> &metric_cofactor, const bool use_invariant_curl_form=false)
Called on the fly and returns the metric cofactor at cubature nodes.
void transform_reference_to_physical(const dealii::Tensor< 1, dim, real > &ref, const dealii::Tensor< 2, dim, real > &metric_cofactor, dealii::Tensor< 1, dim, real > &phys)
Given a reference tensor, return the physical tensor.
void gradient_matrix_vector_mult(const std::vector< real > &input_vect, dealii::Tensor< 1, dim, std::vector< real >> &output_vect, const dealii::FullMatrix< double > &basis_x, const dealii::FullMatrix< double > &basis_y, const dealii::FullMatrix< double > &basis_z, const dealii::FullMatrix< double > &gradient_basis_x, const dealii::FullMatrix< double > &gradient_basis_y, const dealii::FullMatrix< double > &gradient_basis_z)
Computes the gradient of a scalar using sum-factorization.
This is the solution basis , the modal differential opertaor commonly seen in DG defined as ...
ESFR correction matrix without jac dependence.
FR_mass_inv_aux(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input, const Parameters::AllParameters::Flux_Reconstruction_Aux FR_param_input)
Constructor.
std::vector< real > det_Jac_vol
The determinant of the metric Jacobian at volume cubature nodes.
Local mass matrix without jacobian dependence.
vol_projection_operator(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input)
Constructor.
FR_mass(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input, const Parameters::AllParameters::Flux_Reconstruction FR_param_input)
Constructor.
vol_projection_operator_FR(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input, const Parameters::AllParameters::Flux_Reconstruction FR_param_input, const bool store_transpose_input=false)
Constructor.
The DG lifting operator is defined as the operator that lifts inner products of polynomials of some o...
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
void get_c_negative_divided_by_two_FR_parameter(const unsigned int curr_cell_degree, double &c)
Evaluates the flux reconstruction parameter at the bottom limit where the scheme is unstable...
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &face_quadrature)
Assembles the one dimensional norm operator that it is lifted onto.
unsigned int current_degree
Stores the degree of the current poly degree.
void build_metric_Jacobian(const unsigned int n_quad_pts, const std::array< std::vector< real >, dim > &mapping_support_points, const dealii::FullMatrix< double > &basis_x_flux_nodes, const dealii::FullMatrix< double > &basis_y_flux_nodes, const dealii::FullMatrix< double > &basis_z_flux_nodes, const dealii::FullMatrix< double > &grad_basis_x_flux_nodes, const dealii::FullMatrix< double > &grad_basis_y_flux_nodes, const dealii::FullMatrix< double > &grad_basis_z_flux_nodes, std::vector< dealii::Tensor< 2, dim, real >> &local_Jac)
Builds the metric Jacobian evaluated at a vector of points.
unsigned int current_degree
Stores the degree of the current poly degree.
unsigned int current_degree
Stores the degree of the current poly degree.
virtual void build_1D_volume_state_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
dealii::FullMatrix< double > oneD_transpose_vol_operator
Stores the transpose of the operator for fast weight-adjusted solves.
dealii::Tensor< 2, dim, std::vector< real > > metric_Jacobian_vol_cubature
Stores the metric Jacobian at flux nodes.
void get_spectral_difference_parameter(const unsigned int curr_cell_degree, double &c)
Evaluates the spectral difference parameter for flux reconstruction.
void build_1D_gradient_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
The ESFR lifting operator.
const Parameters::AllParameters::Flux_Reconstruction_Aux FR_param_type
Flux reconstruction parameter type.
void build_local_Flux_Reconstruction_operator(const dealii::FullMatrix< double > &local_Mass_Matrix, const dealii::FullMatrix< double > &pth_derivative, const unsigned int n_dofs, const double c, dealii::FullMatrix< double > &Flux_Reconstruction_operator)
Computes a single local Flux_Reconstruction operator (ESFR correction operator) on the fly for a loca...
void sum_factorized_Hadamard_surface_basis_assembly(const unsigned int rows_size, const unsigned int columns_size_1D, const std::vector< unsigned int > &rows, const std::vector< unsigned int > &columns, const dealii::FullMatrix< double > &basis, const std::vector< double > &weights, dealii::FullMatrix< double > &basis_sparse, const int dim_not_zero)
Constructs the basis operator storing all non-zero entries for a "sum-factorized" surface Hadamard p...
double FR_param
Flux reconstruction paramater value.
void two_pt_flux_Hadamard_product(const dealii::FullMatrix< real > &input_mat, dealii::FullMatrix< real > &output_mat, const dealii::FullMatrix< double > &basis, const std::vector< real > &weights, const int direction)
Computes the Hadamard product ONLY for 2pt flux calculations.
unsigned int current_degree
Stores the degree of the current poly degree.
void build_1D_gradient_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
void build_1D_surface_state_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 0 > &face_quadrature)
Assembles the one dimensional operator.
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
Base metric operators class that stores functions used in both the volume and on surface.
void matrix_vector_mult_surface_1D(const unsigned int face_number, const std::vector< real > &input_vect, std::vector< real > &output_vect, const std::array< dealii::FullMatrix< double >, 2 > &basis_surf, const dealii::FullMatrix< double > &basis_vol, const bool adding=false, const double factor=1.0)
Apply sum-factorization matrix vector multiplication on a surface.
void build_determinant_metric_Jacobian(const unsigned int n_quad_pts, const std::array< std::vector< real >, dim > &mapping_support_points, const dealii::FullMatrix< double > &basis_x_flux_nodes, const dealii::FullMatrix< double > &basis_y_flux_nodes, const dealii::FullMatrix< double > &basis_z_flux_nodes, const dealii::FullMatrix< double > &grad_basis_x_flux_nodes, const dealii::FullMatrix< double > &grad_basis_y_flux_nodes, const dealii::FullMatrix< double > &grad_basis_z_flux_nodes, std::vector< real > &det_metric_Jac)
Assembles the determinant of metric Jacobian.
const Parameters::AllParameters::Flux_Reconstruction FR_param_type
Flux reconstruction parameter type.
ESFR correction matrix for AUX EQUATION without jac dependence.
void divergence_matrix_vector_mult_1D(const dealii::Tensor< 1, dim, std::vector< real >> &input_vect, std::vector< real > &output_vect, const dealii::FullMatrix< double > &basis, const dealii::FullMatrix< double > &gradient_basis)
Computes the divergence using sum-factorization where the basis are the same in each direction...
void build_1D_gradient_state_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
local_basis_stiffness(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input, const bool store_skew_symmetric_form_input=false)
Constructor.
The mapping shape functions evaluated at the desired nodes (facet set included in volume grid nodes f...
SumFactorizedOperators(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input)
Precompute 1D operator in constructor.
modal_basis_differential_operator(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input)
Constructor.
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &face_quadrature)
Assembles the one dimensional norm operator that it is lifted onto.
void divergence_matrix_vector_mult(const dealii::Tensor< 1, dim, std::vector< real >> &input_vect, std::vector< real > &output_vect, const dealii::FullMatrix< double > &basis_x, const dealii::FullMatrix< double > &basis_y, const dealii::FullMatrix< double > &basis_z, const dealii::FullMatrix< double > &gradient_basis_x, const dealii::FullMatrix< double > &gradient_basis_y, const dealii::FullMatrix< double > &gradient_basis_z)
Computes the divergence using the sum factorization matrix-vector multiplication. ...
void get_FR_correction_parameter(const unsigned int curr_cell_degree, double &c)
Gets the FR correction parameter for the primary equation and stores.
const bool store_vol_flux_nodes
Flag if store metric Jacobian at flux nodes.
OperatorsBase(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input)
Constructor.
The basis functions separated by nstate with n shape functions.
basis_functions(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input)
Constructor.
void gradient_matrix_vector_mult_1D(const std::vector< real > &input_vect, dealii::Tensor< 1, dim, std::vector< real >> &output_vect, const dealii::FullMatrix< double > &basis, const dealii::FullMatrix< double > &gradient_basis)
Computes the gradient of a scalar using sum-factorization where the basis are the same in each direct...
void build_1D_shape_functions_at_volume_flux_nodes(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Constructs the volume and volume gradient operator.
The metric independent FR mass matrix .
void inner_product(const std::vector< real > &input_vect, const std::vector< real > &weight_vect, std::vector< real > &output_vect, const dealii::FullMatrix< double > &basis_x, const dealii::FullMatrix< double > &basis_y, const dealii::FullMatrix< double > &basis_z, const bool adding=false, const double factor=1.0) override
Computes the inner product between a matrix and a vector multiplied by some weight function...
const Parameters::AllParameters::Flux_Reconstruction FR_param_type
Flux reconstruction parameter type.
void transform_physical_to_reference(const dealii::Tensor< 1, dim, real > &phys, const dealii::Tensor< 2, dim, real > &metric_cofactor, dealii::Tensor< 1, dim, real > &ref)
Given a physical tensor, return the reference tensor.
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
vol_projection_operator_FR_aux(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input, const Parameters::AllParameters::Flux_Reconstruction_Aux FR_param_input, const bool store_transpose_input=false)
Constructor.
dealii::FullMatrix< double > oneD_vol_operator
Stores the one dimensional volume operator.
void sum_factorized_Hadamard_basis_assembly(const unsigned int rows_size_1D, const unsigned int columns_size_1D, const std::vector< std::array< unsigned int, dim >> &rows, const std::vector< std::array< unsigned int, dim >> &columns, const dealii::FullMatrix< double > &basis, const std::vector< double > &weights, std::array< dealii::FullMatrix< double >, dim > &basis_sparse)
Constructs the basis operator storing all non-zero entries for a "sum-factorized" Hadamard product...
Flux_Reconstruction_Aux
Type of correction in Flux Reconstruction for the auxiliary variables.
void build_facet_metric_operators(const unsigned int iface, const unsigned int n_quad_pts, const unsigned int n_metric_dofs, const std::array< std::vector< real >, dim > &mapping_support_points, mapping_shape_functions< dim, n_faces, real > &mapping_basis, const bool use_invariant_curl_form=false)
Builds the facet metric operators.
unsigned int current_degree
Stores the degree of the current poly degree.
void build_1D_surface_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 0 > &quadrature)
Assembles the one dimensional operator.
derivative_p(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input)
Constructor.
std::vector< real > det_Jac_surf
The determinant of the metric Jacobian at facet cubature nodes.
basis_functions< dim, n_faces, real > mapping_shape_functions_flux_nodes
Object of mapping shape functions evaluated at flux nodes.
void get_c_negative_FR_parameter(const unsigned int curr_cell_degree, double &c)
Evaluates the flux reconstruction parameter at the bottom limit where the scheme is unstable...
unsigned int current_degree
Stores the degree of the current poly degree.
const Parameters::AllParameters::Flux_Reconstruction_Aux FR_param_aux_type
Flux reconstruction parameter type.
The metric independent inverse of the FR mass matrix for auxiliary equation .
void get_c_plus_parameter(const unsigned int curr_cell_degree, double &c)
Gets the FR correction parameter corresponding to the maximum timestep.
Projection operator corresponding to basis functions onto -norm.
void build_1D_surface_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 0 > &face_quadrature)
Assembles the one dimensional operator.
unsigned int current_degree
Stores the degree of the current poly degree.
dealii::FullMatrix< double > build_dim_Flux_Reconstruction_operator_directly(const int nstate, const unsigned int n_dofs, dealii::FullMatrix< double > &pth_deriv, dealii::FullMatrix< double > &mass_matrix)
Computes the dim sized flux reconstruction operator for general Mass Matrix (needed for curvilinear)...
void build_volume_metric_operators(const unsigned int n_quad_pts, const unsigned int n_metric_dofs, const std::array< std::vector< real >, dim > &mapping_support_points, mapping_shape_functions< dim, n_faces, real > &mapping_basis, const bool use_invariant_curl_form=false)
Builds the volume metric operators.
void build_1D_surface_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 0 > &face_quadrature)
Assembles the one dimensional operator.
const bool store_skew_symmetric_form
Flag to store the skew symmetric form .
std::array< dealii::Tensor< 1, dim, std::vector< real > >, n_faces > flux_nodes_surf
Stores the physical facet flux nodes.
"Stiffness" operator used in DG Strong form.
const Parameters::AllParameters::Flux_Reconstruction FR_param_type
Flux reconstruction parameter type.
void compute_local_vol_projection_operator(const dealii::FullMatrix< double > &norm_matrix_inverse, const dealii::FullMatrix< double > &integral_vol_basis, dealii::FullMatrix< double > &volume_projection)
Computes a single local projection operator on some space (norm).
unsigned int current_degree
Stores the degree of the current poly degree.
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
mapping_shape_functions(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input)
Constructor.
local_flux_basis_stiffness(const unsigned int max_degree_input, const unsigned int grid_degree_input)
Constructor.
dealii::FullMatrix< double > build_dim_Flux_Reconstruction_operator(const dealii::FullMatrix< double > &local_Mass_Matrix, const int nstate, const unsigned int n_dofs)
Computes the dim sized flux reconstruction operator with simplified tensor product form...
unsigned int current_degree
Stores the degree of the current poly degree.
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
vol_integral_basis(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input)
Constructor.
const Parameters::AllParameters::Flux_Reconstruction FR_param_type
Flux reconstruction parameter type.
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
void build_1D_shape_functions_at_flux_nodes(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature, const dealii::Quadrature< 0 > &face_quadrature)
Constructs the volume, gradient, surface, and surface gradient operator.
void transform_physical_to_reference_vector(const dealii::Tensor< 1, dim, std::vector< real >> &phys, const dealii::Tensor< 2, dim, std::vector< real >> &metric_cofactor, dealii::Tensor< 1, dim, std::vector< real >> &ref)
Given a physical tensor of vector of points, return the reference tensor of vector.
vol_integral_gradient_basis(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input)
Constructor.
void matrix_vector_mult_1D(const std::vector< real > &input_vect, std::vector< real > &output_vect, const dealii::FullMatrix< double > &basis_x, const bool adding=false, const double factor=1.0)
Apply the matrix vector operation using the 1D operator in each direction.
FR_mass_inv(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input, const Parameters::AllParameters::Flux_Reconstruction FR_param_input)
Constructor.
void build_1D_surface_gradient_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 0 > &quadrature)
Assembles the one dimensional operator.
unsigned int current_degree
Stores the degree of the current poly degree.
dealii::Tensor< 1, dim, std::vector< real > > flux_nodes_vol
Stores the physical volume flux nodes.
The surface integral of test functions.
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
unsigned int current_degree
Stores the degree of the current poly degree.
Projection operator corresponding to basis functions onto M-norm (L2).
Local stiffness matrix without jacobian dependence.
flux_basis_functions_state(const unsigned int max_degree_input, const unsigned int grid_degree_input)
Constructor.
void build_1D_volume_state_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
void sum_factorized_Hadamard_sparsity_pattern(const unsigned int rows_size, const unsigned int columns_size, std::vector< std::array< unsigned int, dim >> &rows, std::vector< std::array< unsigned int, dim >> &columns)
Computes the rows and columns vectors with non-zero indices for sum-factorized Hadamard products...
unsigned int current_degree
Stores the degree of the current poly degree.
local_mass(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input)
Constructor.
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.