[P]arallel [Hi]gh-order [Li]brary for [P]DEs  Latest
Parallel High-Order Library for PDEs through hp-adaptive Discontinuous Galerkin methods
operators.cpp
1 #include <deal.II/base/conditional_ostream.h>
2 #include <deal.II/base/parameter_handler.h>
3 
4 #include <deal.II/base/qprojector.h>
5 #include <deal.II/base/geometry_info.h>
6 
7 #include <deal.II/grid/tria.h>
8 
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>
14 
15 
16 #include <deal.II/dofs/dof_handler.h>
17 
18 #include <deal.II/hp/q_collection.h>
19 #include <deal.II/hp/mapping_collection.h>
20 #include <deal.II/hp/fe_values.h>
21 
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>
27 
28 #include <Epetra_RowMatrixTransposer.h>
29 #include <AztecOO.h>
30 
31 #include "ADTypes.hpp"
32 #include <Sacado.hpp>
33 #include <CoDiPack/include/codi.hpp>
34 
35 #include "operators.h"
36 
37 namespace PHiLiP {
38 namespace OPERATOR {
39 
40 //Constructor
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)
52 {}
53 
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)
59 {
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();
66 
67  if constexpr (dim==1)
68  return basis_x;
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];
78  }
79  }
80  }
81  }
82  return tens_prod;
83  }
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];
95  }
96  }
97  }
98  }
99  }
100  }
101  return tens_prod;
102  }
103 }
104 
105 template <int dim, int n_faces, typename real>
107  const int nstate,
108  const dealii::FullMatrix<double> &basis_x,
109  const dealii::FullMatrix<double> &basis_y,
110  const dealii::FullMatrix<double> &basis_z)
111 {
112  //assert that each basis matrix is of size (rows x columns)
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();
119 
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;
126 
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);
134 
135 
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];
143  }
144  }
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];
149  }
150  }
151  }
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];
156  }
157  }
158  }
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];
166  }
167  }
168  }
169  return tens_prod;
170 }
171 
172 template <int dim, int n_faces, typename real>
174 {
175  if ((n==0)||(n==1))
176  return 1;
177  else
178  return n*compute_factorial(n-1);
179 }
180 
181 /**********************************
182 *
183 * Sum Factorization class
184 *
185 **********************************/
186 //Constructor
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)
192  : OperatorsBase<dim,n_faces,real>::OperatorsBase(nstate_input, max_degree_input, grid_degree_input)
193 {}
194 
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,
202  const bool adding,
203  const double factor)
204 {
205  //assert that each basis matrix is of size (rows x columns)
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());
215  }
216  if constexpr (dim == 2){
217  assert(rows_x * rows_y == output_vect.size());
218  assert(columns_x * columns_y == input_vect.size());
219  }
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());
223  }
224 
225  if constexpr (dim==1){
226  for(unsigned int iquad=0; iquad<rows_x; iquad++){
227  if(!adding)
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];
231  }
232  }
233  }
234  if constexpr (dim==2){
235  //convert the input vector to matrix
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];//jdof runs fastest (x) idof slowest (y)
240  }
241  }
242  dealii::FullMatrix<double> temp(rows_x, columns_y);
243  basis_x.mmult(temp, input_mat);//apply x tensor product
244  dealii::FullMatrix<double> output_mat(rows_y, rows_x);
245  basis_y.mTmult(output_mat, temp);//apply y tensor product
246  //convert mat back to vect
247  for(unsigned int iquad=0; iquad<rows_y; iquad++){
248  for(unsigned int jquad=0; jquad<rows_x; jquad++){
249  if(adding)
250  output_vect[iquad * rows_x + jquad] += factor * output_mat[iquad][jquad];
251  else
252  output_vect[iquad * rows_x + jquad] = factor * output_mat[iquad][jquad];
253  }
254  }
255 
256  }
257  if constexpr (dim==3){
258  //convert vect to mat first
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];//kdof runs fastest (x) idof slowest (z)
265  }
266  }
267  }
268  dealii::FullMatrix<double> temp(rows_x, columns_y * columns_z);
269  basis_x.mmult(temp, input_mat);//apply x tensor product
270  //convert to have y dofs ie/ change the stride
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];//extract y runs second fastest
276  }
277  }
278  }
279  dealii::FullMatrix<double> temp3(rows_y, rows_x * columns_z);
280  basis_y.mmult(temp3, temp2);//apply y tensor product
281  dealii::FullMatrix<double> temp4(columns_z, rows_x * rows_y);
282  //convert to have z dofs ie/ change the stride
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];//extract z runs slowest
287  }
288  }
289  }
290  dealii::FullMatrix<double> output_mat(rows_z, rows_x * rows_y);
291  basis_z.mmult(output_mat, temp4);
292  //convert mat to vect
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;
297  if(adding)
298  output_vect[quad_index] += factor * output_mat[iquad][kquad * rows_y + jquad];
299  else
300  output_vect[quad_index] = factor * output_mat[iquad][kquad * rows_y + jquad];
301  }
302  }
303  }
304  }
305 }
306 
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,
312  const bool adding,
313  const double factor)
314 {
315  this->matrix_vector_mult(input_vect, output_vect, basis_x, basis_x, basis_x, adding, factor);
316 }
317 
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,
325  const bool adding,
326  const double factor)
327 {
328  if(face_number == 0)
329  this->matrix_vector_mult(input_vect, output_vect, basis_surf[0], basis_vol, basis_vol, adding, factor);
330  if(face_number == 1)
331  this->matrix_vector_mult(input_vect, output_vect, basis_surf[1], basis_vol, basis_vol, adding, factor);
332  if(face_number == 2)
333  this->matrix_vector_mult(input_vect, output_vect, basis_vol, basis_surf[0], basis_vol, adding, factor);
334  if(face_number == 3)
335  this->matrix_vector_mult(input_vect, output_vect, basis_vol, basis_surf[1], basis_vol, adding, factor);
336  if(face_number == 4)
337  this->matrix_vector_mult(input_vect, output_vect, basis_vol, basis_vol, basis_surf[0], adding, factor);
338  if(face_number == 5)
339  this->matrix_vector_mult(input_vect, output_vect, basis_vol, basis_vol, basis_surf[1], adding, factor);
340 }
341 
342 
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,
351  const bool adding,
352  const double factor)
353 {
354  if(face_number == 0)
355  this->inner_product(input_vect, weight_vect, output_vect, basis_surf[0], basis_vol, basis_vol, adding, factor);
356  if(face_number == 1)
357  this->inner_product(input_vect, weight_vect, output_vect, basis_surf[1], basis_vol, basis_vol, adding, factor);
358  if(face_number == 2)
359  this->inner_product(input_vect, weight_vect, output_vect, basis_vol, basis_surf[0], basis_vol, adding, factor);
360  if(face_number == 3)
361  this->inner_product(input_vect, weight_vect, output_vect, basis_vol, basis_surf[1], basis_vol, adding, factor);
362  if(face_number == 4)
363  this->inner_product(input_vect, weight_vect, output_vect, basis_vol, basis_vol, basis_surf[0], adding, factor);
364  if(face_number == 5)
365  this->inner_product(input_vect, weight_vect, output_vect, basis_vol, basis_vol, basis_surf[1], adding, factor);
366 }
367 
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)
374 {
375  divergence_matrix_vector_mult(input_vect, output_vect,
376  basis, basis, basis,
377  gradient_basis, gradient_basis, gradient_basis);
378 }
379 
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)
390 {
391  for(int idim=0; idim<dim;idim++){
392  if(idim==0)
393  this->matrix_vector_mult(input_vect[idim], output_vect,
394  gradient_basis_x,
395  basis_y,
396  basis_z,
397  false);//first one doesn't add in the divergence
398  if(idim==1)
399  this->matrix_vector_mult(input_vect[idim], output_vect,
400  basis_x,
401  gradient_basis_y,
402  basis_z,
403  true);
404  if(idim==2)
405  this->matrix_vector_mult(input_vect[idim], output_vect,
406  basis_x,
407  basis_y,
408  gradient_basis_z,
409  true);
410  }
411 }
412 
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)
419 {
420  gradient_matrix_vector_mult(input_vect, output_vect,
421  basis, basis, basis,
422  gradient_basis, gradient_basis, gradient_basis);
423 }
424 
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)
435 {
436  for(int idim=0; idim<dim;idim++){
437 // output_vect[idim].resize(input_vect.size());
438  if(idim==0)
439  this->matrix_vector_mult(input_vect, output_vect[idim],
440  gradient_basis_x,
441  basis_y,
442  basis_z,
443  false);//first one doesn't add in the divergence
444  if(idim==1)
445  this->matrix_vector_mult(input_vect, output_vect[idim],
446  basis_x,
447  gradient_basis_y,
448  basis_z,
449  false);
450  if(idim==2)
451  this->matrix_vector_mult(input_vect, output_vect[idim],
452  basis_x,
453  basis_y,
454  gradient_basis_z,
455  false);
456  }
457 }
458 
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,
467  const bool adding,
468  const double factor)
469 {
470  //assert that each basis matrix is of size (rows x columns)
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();
477  //Note the assertion has columns to output and rows to input
478  //bc we transpose the basis inputted for the inner product
479  if constexpr (dim == 1){
480  assert(rows_x == input_vect.size());
481  assert(columns_x == output_vect.size());
482  }
483  if constexpr (dim == 2){
484  assert(rows_x * rows_y == input_vect.size());
485  assert(columns_x * columns_y == output_vect.size());
486  }
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());
490  }
491  assert(weight_vect.size() == input_vect.size());
492 
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);
496 
497  //set as the transpose as inputed basis
498  //found an issue with Tadd for arbitrary size so I manually do it here.
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];
502  }
503  }
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];
507  }
508  }
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];
512  }
513  }
514 
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];
518  }
519 
520  this->matrix_vector_mult(new_input_vect, output_vect, basis_x_trans, basis_y_trans, basis_z_trans, adding, factor);
521 }
522 
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,
529  const bool adding,
530  const double factor)
531 {
532  this->inner_product(input_vect, weight_vect, output_vect, basis_x, basis_x, basis_x, adding, factor);
533 }
534 
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)
542 {
543  assert(input_mat[0].m() == output_vect.size());
544 
545  dealii::FullMatrix<double> output_mat(input_mat[0].m(), input_mat[0].n());
546  for(int idim=0; idim<dim; idim++){
547  two_pt_flux_Hadamard_product(input_mat[idim], output_mat, basis, weights, idim);
548  if constexpr(dim==1){
549  for(unsigned int row=0; row<input_mat[0].m(); row++){//n^d rows
550  for(unsigned int col=0; col<basis.m(); col++){//only need to sum n columns
551  const unsigned int col_index = col;
552  output_vect[row] += scaling * output_mat[row][col_index];//scaled by 2.0 for 2pt flux
553  }
554  }
555  }
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++){
562  if(idim==0){
563  const unsigned int col_index = col + irow * size_1D;
564  output_vect[row_index] += scaling * output_mat[row_index][col_index];//scaled by 2.0 for 2pt flux
565  }
566  if(idim==1){
567  const unsigned int col_index = col * size_1D + jrow;
568  output_vect[row_index] += scaling * output_mat[row_index][col_index];//scaled by 2.0 for 2pt flux
569  }
570  }
571  }
572  }
573  }
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++){
581  if(idim==0){
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];//scaled by 2.0 for 2pt flux
584  }
585  if(idim==1){
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];//scaled by 2.0 for 2pt flux
588  }
589  if(idim==2){
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];//scaled by 2.0 for 2pt flux
592  }
593  }
594  }
595  }
596  }
597  }
598  }
599 }
600 
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)//scaling is unit_ref_normal[dim_not_zero]
611 {
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;
615 
616  dealii::FullMatrix<double> output_mat(input_mat.m(), input_mat.n());
617  two_pt_flux_Hadamard_product(input_mat, output_mat, surf_basis[iface_1D], weights, dim_not_zero);
618  if constexpr(dim==1){
619  for(unsigned int row=0; row<surf_basis[iface_1D].m(); row++){//n rows
620  for(unsigned int col=0; col<surf_basis[iface_1D].n(); col++){//only need to sum n columns
621  output_vect_vol[col] += scaling
622  * output_mat[row][col];//scaled by 2.0 for 2pt flux
623  output_vect_surf[row] -= scaling //minus because skew-symm form
624  * output_mat[row][col];//scaled by 2.0 for 2pt flux
625  }
626  }
627  }
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++){
635  if(dim_not_zero==0){
636  const unsigned int col_index = col + irow * size_1D_col;
637  output_vect_vol[col_index] += scaling * output_mat[row_index][col_index];//scaled by 2.0 for 2pt flux
638  output_vect_surf[row_index] -= scaling * output_mat[row_index][col_index];//scaled by 2.0 for 2pt flux
639  }
640  if(dim_not_zero==1){
641  const unsigned int col_index = col * size_1D_col + irow;
642  output_vect_vol[col_index] += scaling * output_mat[row_index][col_index];//scaled by 2.0 for 2pt flux
643  output_vect_surf[row_index] -= scaling * output_mat[row_index][col_index];//scaled by 2.0 for 2pt flux
644  }
645  }
646  }
647  }
648  }
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++){
657  if(dim_not_zero==0){
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];//scaled by 2.0 for 2pt flux
660  output_vect_surf[row_index] -= scaling * output_mat[row_index][col_index];//scaled by 2.0 for 2pt flux
661  }
662  if(dim_not_zero==1){
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];//scaled by 2.0 for 2pt flux
665  output_vect_surf[row_index] -= scaling * output_mat[row_index][col_index];//scaled by 2.0 for 2pt flux
666  }
667  if(dim_not_zero==2){
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];//scaled by 2.0 for 2pt flux
670  output_vect_surf[row_index] -= scaling * output_mat[row_index][col_index];//scaled by 2.0 for 2pt flux
671  }
672  }
673  }
674  }
675  }
676  }
677 }
678 
679 
680 
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,
687  const int direction)
688 {
689  assert(input_mat.size() == output_mat.size());
690  const unsigned int size = basis.n();
691  assert(size == weights.size());
692 
693  if constexpr(dim == 1){
694  Hadamard_product(input_mat, basis, output_mat);
695  }
696  if constexpr(dim == 2){
697  //In the general case, the basis is non-square (think surface lifting functions).
698  //We assume the other directions are square but variable in the basis of interest.
699  const unsigned int rows = basis.m();
700  assert(rows == input_mat.m());
701  if(direction == 0){
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);
706  //fill index range for diagonal blocks of rize rows_x x columns_x
707  std::iota(row_index.begin(), row_index.end(), idiag*rows);
708  std::iota(col_index.begin(), col_index.end(), idiag*size);
709  //extract diagonal block from input matrix
710  local_block.extract_submatrix_from(input_mat, row_index, col_index);
711  dealii::FullMatrix<double> local_Hadamard(rows, size);
712  Hadamard_product(local_block, basis, local_Hadamard);
713  //scale by the diagonal weight from tensor product
714  local_Hadamard *= weights[idiag];
715  //write the values into diagonal block of output
716  local_Hadamard.scatter_matrix_to(row_index, col_index, output_mat);
717  }
718  }
719  if(direction == 1){
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]
727  * weights[kdiag];
728  }
729  }
730  }
731 
732  }
733  }
734  if constexpr(dim == 3){
735  const unsigned int rows = basis.m();
736  if(direction == 0){
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);
743  //fill index range for diagonal blocks of rize rows_x x columns_x
744  std::iota(row_index.begin(), row_index.end(), idiag*rows);
745  std::iota(col_index.begin(), col_index.end(), idiag*size);
746  //extract diagonal block from input matrix
747  local_block.extract_submatrix_from(input_mat, row_index, col_index);
748  dealii::FullMatrix<double> local_Hadamard(rows, size);
749  Hadamard_product(local_block, basis, local_Hadamard);
750  //scale by the diagonal weight from tensor product
751  local_Hadamard *= weights[kdiag];
752  kdiag++;
753  const unsigned int jdiag = idiag / size;
754  local_Hadamard *= weights[jdiag];
755  //write the values into diagonal block of output
756  local_Hadamard.scatter_matrix_to(row_index, col_index, output_mat);
757  }
758  }
759  if(direction == 1){
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]
768  * weights[zdiag]
769  * weights[kdiag];
770  }
771  }
772  }
773  }
774  }
775  if(direction == 2){
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]
784  * weights[kdiag]
785  * weights[jdiag];
786  }
787  }
788  }
789  }
790  }
791  }
792 }
793 
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)
799 {
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());
804 
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];
809  }
810  }
811 }
812 
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)
819 {
820  //Note that for all directions, the rows vector should always be the same.
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;
827  }
828  }
829  }
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;
838  //direction 0
839  const unsigned int col_index_0 = idiag * columns_size;
840  columns[array_index][0] = col_index_0 + kdiag;
841  //direction 1
842  const unsigned int col_index_1 = kdiag * columns_size;
843  columns[array_index][1] = col_index_1 + jdiag;
844  }
845  }
846  }
847  }
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
856  + ldiag;
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;
862  //direction 0
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;
866  //direction 1
867  const unsigned int col_index_1 = ldiag * columns_size;
868  columns[array_index][1] = col_index_1 + kdiag + idiag * columns_size * columns_size;
869  //direction 2
870  const unsigned int col_index_2 = ldiag * columns_size * columns_size;
871  columns[array_index][2] = col_index_2 + kdiag + jdiag * columns_size;
872  }
873  }
874  }
875  }
876  }
877 }
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)
887 {
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];
892  }
893  }
894  }
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){
899  counter = 0;
900  }
901  //direction 0
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];
904  //direction 1
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];
907  }
908  }
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){
913  counter = 0;
914  }
915  //direction 0
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];
919  //direction 1
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];
923  //direction 2
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];
927  }
928  }
929 
930 }
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)
938 {
939  //Note that for all directions, the rows vector should always be the same.
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;
946  }
947  }
948  }
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 ;
953 
954  const unsigned int row_index = idiag;
955  rows[array_index] = row_index;
956  //direction 0
957  if(dim_not_zero == 0){
958  const unsigned int col_index_0 = idiag * columns_size;
959  columns[array_index] = col_index_0 + jdiag;
960  }
961  //direction 1
962  if(dim_not_zero == 1){
963  const unsigned int col_index_1 = jdiag * columns_size;
964  columns[array_index] = col_index_1 + idiag;
965  }
966  }
967  }
968  }
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
975  + kdiag;
976  const unsigned int row_index = idiag * columns_size + jdiag;
977  rows[array_index] = row_index;
978  //direction 0
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;
983  }
984  //direction 1
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;
988  }
989  //direction 2
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;
993  }
994  }
995  }
996  }
997  }
998 }
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)
1009 {
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];
1014  }
1015  }
1016  }
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){
1021  counter = 0;
1022  }
1023  //direction 0
1024  if(dim_not_zero == 0){
1025  basis_sparse[rows[index]][counter] = basis[0][columns[index]%columns_size_1D]//oneD surf basis only 1 point
1026  * weights[rows[index]%columns_size_1D];
1027  }
1028  //direction 1
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];
1032  }
1033  }
1034  }
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){
1039  counter = 0;
1040  }
1041  //direction 0
1042  if(dim_not_zero == 0){
1043  basis_sparse[rows[index]][counter] = basis[0][columns[index]%columns_size_1D]//oneD surf basis only 1 point
1044  * weights[rows[index]%columns_size_1D]
1045  * weights[rows[index]/columns_size_1D];
1046  }
1047  //direction 1
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];
1052  }
1053  //direction 2
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];
1058  }
1059  }
1060  }
1061 }
1062 
1063 /*******************************************
1064  *
1065  * VOLUME OPERATORS FUNCTIONS
1066  *
1067  *
1068  ******************************************/
1069 
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)
1075  : SumFactorizedOperators<dim,n_faces,real>::SumFactorizedOperators(nstate_input, max_degree_input, grid_degree_input)
1076 {
1077  //Initialize to the max degrees
1078  current_degree = max_degree_input;
1079 }
1080 
1081 template <int dim, int n_faces, typename real>
1083  const dealii::FESystem<1,1> &finite_element,
1084  const dealii::Quadrature<1> &quadrature)
1085 {
1086  const unsigned int n_quad_pts = quadrature.size();
1087  const unsigned int n_dofs = finite_element.dofs_per_cell;
1088  //allocate the basis at volume cubature
1089  this->oneD_vol_operator.reinit(n_quad_pts, n_dofs);
1090  //loop and store
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;
1095  //Basis function idof of poly degree idegree evaluated at cubature node qpoint.
1096  this->oneD_vol_operator[iquad][idof] = finite_element.shape_value_component(idof,qpoint,istate);
1097  }
1098  }
1099 }
1100 
1101 template <int dim, int n_faces, typename real>
1103  const dealii::FESystem<1,1> &finite_element,
1104  const dealii::Quadrature<1> &quadrature)
1105 {
1106  const unsigned int n_quad_pts = quadrature.size();
1107  const unsigned int n_dofs = finite_element.dofs_per_cell;
1108  //allocate the basis at volume cubature
1109  this->oneD_grad_operator.reinit(n_quad_pts, n_dofs);
1110  //loop and store
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;
1115  //Basis function idof of poly degree idegree evaluated at cubature node qpoint.
1116  this->oneD_grad_operator[iquad][idof] = finite_element.shape_grad_component(idof,qpoint,istate)[0];
1117  }
1118  }
1119 }
1120 
1121 template <int dim, int n_faces, typename real>
1123  const dealii::FESystem<1,1> &finite_element,
1124  const dealii::Quadrature<0> &face_quadrature)
1125 {
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;
1129  //loop and store
1130  for(unsigned int iface=0; iface<n_faces_1D; iface++){
1131  //allocate the facet operator
1132  this->oneD_surf_operator[iface].reinit(n_face_quad_pts, n_dofs);
1133  //sum factorized operators use a 1D element.
1134  const dealii::Quadrature<1> quadrature = dealii::QProjector<1>::project_to_face(dealii::ReferenceCell::get_hypercube(1),
1135  face_quadrature,
1136  iface);
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;
1141  //Basis function idof of poly degree idegree evaluated at cubature node qpoint.
1142  this->oneD_surf_operator[iface][iquad][idof] = finite_element.shape_value_component(idof,qpoint,istate);
1143  }
1144  }
1145  }
1146 }
1147 
1148 template <int dim, int n_faces, typename real>
1150  const dealii::FESystem<1,1> &finite_element,
1151  const dealii::Quadrature<0> &face_quadrature)
1152 {
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;
1156  //loop and store
1157  for(unsigned int iface=0; iface<n_faces_1D; iface++){
1158  //allocate the facet operator
1159  this->oneD_surf_grad_operator[iface].reinit(n_face_quad_pts, n_dofs);
1160  //sum factorized operators use a 1D element.
1161  const dealii::Quadrature<1> quadrature = dealii::QProjector<1>::project_to_face(dealii::ReferenceCell::get_hypercube(1),
1162  face_quadrature,
1163  iface);
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;
1168  //Basis function idof of poly degree idegree evaluated at cubature node qpoint.
1169  this->oneD_surf_grad_operator[iface][iquad][idof] = finite_element.shape_grad_component(idof,qpoint,istate)[0];
1170  }
1171  }
1172  }
1173 }
1174 
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)
1180  : SumFactorizedOperators<dim,n_faces,real>::SumFactorizedOperators(nstate_input, max_degree_input, grid_degree_input)
1181 {
1182  //Initialize to the max degrees
1183  current_degree = max_degree_input;
1184 }
1185 
1186 template <int dim, int n_faces, typename real>
1188  const dealii::FESystem<1,1> &finite_element,
1189  const dealii::Quadrature<1> &quadrature)
1190 {
1191  const unsigned int n_quad_pts = quadrature.size();
1192  const unsigned int n_dofs = finite_element.dofs_per_cell;
1193  //allocate
1194  this->oneD_vol_operator.reinit(n_quad_pts, n_dofs);
1195  //loop and store
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;
1201  //Basis function idof of poly degree idegree evaluated at cubature node qpoint multiplied by quad weight.
1202  this->oneD_vol_operator[iquad][idof] = quad_weights[iquad] * finite_element.shape_value_component(idof,qpoint,istate);
1203  }
1204  }
1205 }
1206 
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)
1212  : SumFactorizedOperators<dim,n_faces,real>::SumFactorizedOperators(nstate_input, max_degree_input, grid_degree_input)
1213 {
1214  //Initialize to the max degrees
1215  current_degree = max_degree_input;
1216 }
1217 
1218 template <int dim, int n_faces, typename real>
1220  const dealii::FESystem<1,1> &finite_element,
1221  const dealii::Quadrature<1> &quadrature)
1222 {
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 ();
1226  //allocate
1227  this->oneD_vol_operator.reinit(n_dofs,n_dofs);
1228  //loop and store
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;
1233  double value = 0.0;
1234  for (unsigned int iquad=0; iquad<n_quad_pts; ++iquad) {
1235  const dealii::Point<1> qpoint = quadrature.point(iquad);
1236  value +=
1237  finite_element.shape_value_component(itest,qpoint,istate_test)
1238  * finite_element.shape_value_component(itrial,qpoint,istate_trial)
1239  * quad_weights[iquad];
1240  }
1241 
1242  this->oneD_vol_operator[itrial][itest] = 0.0;
1243  this->oneD_vol_operator[itest][itrial] = 0.0;
1244  if(istate_test==istate_trial) {
1245  this->oneD_vol_operator[itrial][itest] = value;
1246  this->oneD_vol_operator[itest][itrial] = value;
1247  }
1248  }
1249  }
1250 }
1251 
1252 template <int dim, int n_faces, typename real>
1254  const int nstate,
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)
1259 
1260 {
1261  const unsigned int n_shape_fns = n_dofs / nstate;
1262  assert(nstate*pow(basis.oneD_vol_operator.m() / nstate, dim) == n_dofs);
1263  dealii::FullMatrix<double> mass_matrix_dim(n_dofs);
1264  dealii::FullMatrix<double> basis_dim(n_dofs);
1265  basis_dim = basis.tensor_product_state(
1266  nstate,
1267  basis.oneD_vol_operator,
1268  basis.oneD_vol_operator,
1269  basis.oneD_vol_operator);
1270  //loop and store
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) {
1274  double value = 0.0;
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]
1280  * det_Jac[iquad]
1281  * quad_weights[iquad];
1282  }
1283  mass_matrix_dim[trial_index][test_index] = value;
1284  mass_matrix_dim[test_index][trial_index] = value;
1285  }
1286  }
1287  }
1288  return mass_matrix_dim;
1289 }
1290 
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)
1297  : SumFactorizedOperators<dim,n_faces,real>::SumFactorizedOperators(nstate_input, max_degree_input, grid_degree_input)
1298  , store_skew_symmetric_form(store_skew_symmetric_form_input)
1299 {
1300  //Initialize to the max degrees
1301  current_degree = max_degree_input;
1302 }
1303 
1304 template <int dim, int n_faces, typename real>
1306  const dealii::FESystem<1,1> &finite_element,
1307  const dealii::Quadrature<1> &quadrature)
1308 {
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 ();
1312  //allocate
1313  this->oneD_vol_operator.reinit(n_dofs,n_dofs);
1314  //loop and store
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;
1319  double value = 0.0;
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]//since it's a 1D operator
1324  * quad_weights[iquad];
1325  }
1326  if(istate == istate_test){
1327  this->oneD_vol_operator[itest][idof] = value;
1328  }
1329  }
1330  }
1332  //allocate
1333  oneD_skew_symm_vol_oper.reinit(n_dofs,n_dofs);
1334  //solve
1335  for(unsigned int idof=0; idof<n_dofs; idof++){
1336  for(unsigned int jdof=0; jdof<n_dofs; jdof++){
1337  oneD_skew_symm_vol_oper[idof][jdof] = this->oneD_vol_operator[idof][jdof]
1338  - this->oneD_vol_operator[jdof][idof];
1339  }
1340  }
1341  }
1342 }
1343 
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)
1349  : SumFactorizedOperators<dim,n_faces,real>::SumFactorizedOperators(nstate_input, max_degree_input, grid_degree_input)
1350 {
1351  //Initialize to the max degrees
1352  current_degree = max_degree_input;
1353 }
1354 
1355 template <int dim, int n_faces, typename real>
1357  const dealii::FESystem<1,1> &finite_element,
1358  const dealii::Quadrature<1> &quadrature)
1359 {
1360  const unsigned int n_dofs = finite_element.dofs_per_cell;
1361  local_mass<dim,n_faces,real> mass_matrix(this->nstate, this->max_degree, this->max_grid_degree);
1362  mass_matrix.build_1D_volume_operator(finite_element, quadrature);
1364  stiffness.build_1D_volume_operator(finite_element, quadrature);
1365  //allocate
1366  this->oneD_vol_operator.reinit(n_dofs,n_dofs);
1367  dealii::FullMatrix<double> inv_mass(n_dofs);
1368  inv_mass.invert(mass_matrix.oneD_vol_operator);
1369  //solves
1370  inv_mass.mmult(this->oneD_vol_operator, stiffness.oneD_vol_operator);
1371 }
1372 
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)
1378  : SumFactorizedOperators<dim,n_faces,real>::SumFactorizedOperators(nstate_input, max_degree_input, grid_degree_input)
1379 {
1380  //Initialize to the max degrees
1381  current_degree = max_degree_input;
1382 }
1383 
1384 template <int dim, int n_faces, typename real>
1386  const dealii::FESystem<1,1> &finite_element,
1387  const dealii::Quadrature<1> &quadrature)
1388 {
1389  const unsigned int n_dofs = finite_element.dofs_per_cell;
1390  //allocate
1391  this->oneD_vol_operator.reinit(n_dofs,n_dofs);
1392  //set as identity
1393  for(unsigned int idof=0; idof<n_dofs; idof++){
1394  this->oneD_vol_operator[idof][idof] = 1.0;//set it equal to identity
1395  }
1396  //get modal basis differential operator
1398  diff_oper.build_1D_volume_operator(finite_element, quadrature);
1399  //loop and solve
1400  for(unsigned int idegree=0; idegree< this->max_degree; idegree++){
1401  dealii::FullMatrix<double> derivative_p_temp(n_dofs, n_dofs);
1402  derivative_p_temp.add(1.0, this->oneD_vol_operator);
1403  diff_oper.oneD_vol_operator.mmult(this->oneD_vol_operator, derivative_p_temp);
1404  }
1405 }
1406 
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,
1413  : SumFactorizedOperators<dim,n_faces,real>::SumFactorizedOperators(nstate_input, max_degree_input, grid_degree_input)
1414  , FR_param_type(FR_param_input)
1415 {
1416  //Initialize to the max degrees
1417  current_degree = max_degree_input;
1418  //get the FR corrcetion parameter value
1420 }
1421 
1422 template <int dim, int n_faces, typename real>
1424  const unsigned int curr_cell_degree,
1425  double &c)
1426 {
1427  const double pfact = this->compute_factorial(curr_cell_degree);
1428  const double pfact2 = this->compute_factorial(2.0 * curr_cell_degree);
1429  double cp = pfact2/(pow(pfact,2));//since ref element [0,1]
1430  c = 2.0 * (curr_cell_degree+1)/( curr_cell_degree*((2.0*curr_cell_degree+1.0)*(pow(pfact*cp,2))));
1431  c/=2.0;//since orthonormal
1432 }
1433 template <int dim, int n_faces, typename real>
1435  const unsigned int curr_cell_degree,
1436  double &c)
1437 {
1438  const double pfact = this->compute_factorial(curr_cell_degree);
1439  const double pfact2 = this->compute_factorial(2.0 * 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))));
1442  c/=2.0;//since orthonormal
1443 }
1444 template <int dim, int n_faces, typename real>
1446  const unsigned int curr_cell_degree,
1447  double &c)
1448 {
1449  const double pfact = this->compute_factorial(curr_cell_degree);
1450  const double pfact2 = this->compute_factorial(2.0 * 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));
1453  c/=2.0;//since orthonormal
1454 }
1455 template <int dim, int n_faces, typename real>
1457  const unsigned int curr_cell_degree,
1458  double &c)
1459 {
1460  get_c_negative_FR_parameter(curr_cell_degree, c);
1461  c/=2.0;
1462 }
1463 template <int dim, int n_faces, typename real>
1465  const unsigned int curr_cell_degree,
1466  double &c)
1467 {
1468  if(curr_cell_degree == 2){
1469  c = 0.186;
1470 // c = 0.173;//RK33
1471  }
1472  else if(curr_cell_degree == 3){
1473  c = 3.67e-3;
1474  }
1475  else if(curr_cell_degree == 4){
1476  c = 4.79e-5;
1477 // c = 4.92e-5;//RK33
1478  }
1479  else if(curr_cell_degree == 5){
1480  c = 4.24e-7;
1481  }
1482  else{
1483  this->pcout << "ERROR: cPlus values are only defined for p=2 through p=5. Aborting..." << std::endl;
1484  std::abort();
1485  }
1486 
1487  c/=2.0;//since orthonormal
1488  c/=pow(pow(2.0,curr_cell_degree),2);//since ref elem [0,1]
1489 }
1490 
1491 template <int dim, int n_faces, typename real>
1493  const unsigned int curr_cell_degree,
1494  double &c)
1495 {
1497  if(FR_param_type == FR_enum::cHU || FR_param_type == FR_enum::cHULumped){
1498  get_Huynh_g2_parameter(curr_cell_degree, FR_param);
1499  }
1500  else if(FR_param_type == FR_enum::cSD){
1501  get_spectral_difference_parameter(curr_cell_degree, c);
1502  }
1503  else if(FR_param_type == FR_enum::cNegative){
1504  get_c_negative_FR_parameter(curr_cell_degree, c);
1505  }
1506  else if(FR_param_type == FR_enum::cNegative2){
1507  get_c_negative_divided_by_two_FR_parameter(curr_cell_degree, c);
1508  }
1509  else if(FR_param_type == FR_enum::cDG){
1510  //DG case is the 0.0 case.
1511  c = 0.0;
1512  }
1513  else if(FR_param_type == FR_enum::c10Thousand){
1514  //Set the value to 10000 for arbitrary high-numbers.
1515  c = 10000.0;
1516  }
1517  else if(FR_param_type == FR_enum::cPlus){
1518  get_c_plus_parameter(curr_cell_degree, c);
1519  }
1520 }
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,
1526  const double c,
1527  dealii::FullMatrix<double> &Flux_Reconstruction_operator)
1528 {
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);
1534 }
1535 
1536 
1537 template <int dim, int n_faces, typename real>
1539  const dealii::FESystem<1,1> &finite_element,
1540  const dealii::Quadrature<1> &quadrature)
1541 {
1542  const unsigned int n_dofs = finite_element.dofs_per_cell;
1543  //allocate the volume operator
1544  this->oneD_vol_operator.reinit(n_dofs, n_dofs);
1545  //build the FR correction operator
1546  derivative_p<dim,n_faces,real> pth_derivative(this->nstate, this->max_degree, this->max_grid_degree);
1547  pth_derivative.build_1D_volume_operator(finite_element, quadrature);
1548 
1549  local_mass<dim,n_faces,real> local_Mass_Matrix(this->nstate, this->max_degree, this->max_grid_degree);
1550  local_Mass_Matrix.build_1D_volume_operator(finite_element, quadrature);
1551  //solves
1552  build_local_Flux_Reconstruction_operator(local_Mass_Matrix.oneD_vol_operator, pth_derivative.oneD_vol_operator, n_dofs, FR_param, this->oneD_vol_operator);
1553 }
1554 
1555 template <int dim, int n_faces, typename real>
1557  const int nstate,
1558  const unsigned int n_dofs,
1559  dealii::FullMatrix<double> &pth_deriv,
1560  dealii::FullMatrix<double> &mass_matrix)
1561 {
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);
1566  if(idim==0){
1567  pth_deriv_dim = this->tensor_product_state(nstate, pth_deriv, identity, identity);
1568  }
1569  if(idim==1){
1570  pth_deriv_dim = this->tensor_product_state(nstate, identity, pth_deriv, identity);
1571  }
1572  if(idim==2){
1573  pth_deriv_dim = this->tensor_product_state(nstate, identity, identity, pth_deriv);
1574  }
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);
1580  }
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);
1586  if(idim==0){
1587  pth_deriv_dim = this->tensor_product_state(nstate, pth_deriv, pth_deriv, identity);
1588  }
1589  if(idim==1){
1590  pth_deriv_dim = this->tensor_product_state(nstate, identity, pth_deriv, pth_deriv);
1591  }
1592  if(idim==2){
1593  pth_deriv_dim = this->tensor_product_state(nstate, pth_deriv, identity, pth_deriv);
1594  }
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);
1600  }
1601  }
1602  if constexpr (dim == 3){
1603  double FR_param_cubed = pow(FR_param,3.0);
1604  dealii::FullMatrix<double> pth_deriv_dim(n_dofs);
1605  pth_deriv_dim = this->tensor_product_state(nstate, pth_deriv, pth_deriv, pth_deriv);
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);
1611  }
1612 
1613  return Flux_Reconstruction_operator;
1614 }
1615 
1616 template <int dim, int n_faces, typename real>
1618  const dealii::FullMatrix<double> &local_Mass_Matrix,
1619  const int nstate,
1620  const unsigned int n_dofs)
1621 {
1622  dealii::FullMatrix<double> dim_FR_operator(n_dofs);
1623  if constexpr (dim == 1){
1624  dim_FR_operator = this->oneD_vol_operator;
1625  }
1626  if (dim >= 2){
1627  dealii::FullMatrix<double> FR1(n_dofs);
1628  FR1 = this->tensor_product_state(nstate, this->oneD_vol_operator, local_Mass_Matrix, local_Mass_Matrix);
1629  dealii::FullMatrix<double> FR2(n_dofs);
1630  FR2 = this->tensor_product_state(nstate, local_Mass_Matrix, this->oneD_vol_operator, local_Mass_Matrix);
1631  dealii::FullMatrix<double> FR_cross1(n_dofs);
1632  FR_cross1 = this->tensor_product_state(nstate, this->oneD_vol_operator, this->oneD_vol_operator, local_Mass_Matrix);
1633  dim_FR_operator.add(1.0, FR1, 1.0, FR2, 1.0, FR_cross1);
1634  }
1635  if constexpr (dim == 3){
1636  dealii::FullMatrix<double> FR3(n_dofs);
1637  FR3 = this->tensor_product_state(nstate, local_Mass_Matrix, local_Mass_Matrix, this->oneD_vol_operator);
1638  dealii::FullMatrix<double> FR_cross2(n_dofs);
1639  FR_cross2 = this->tensor_product_state(nstate, this->oneD_vol_operator, local_Mass_Matrix, this->oneD_vol_operator);
1640  dealii::FullMatrix<double> FR_cross3(n_dofs);
1641  FR_cross3 = this->tensor_product_state(nstate, local_Mass_Matrix, this->oneD_vol_operator, this->oneD_vol_operator);
1642  dealii::FullMatrix<double> FR_triple(n_dofs);
1643  FR_triple = this->tensor_product_state(nstate, this->oneD_vol_operator, this->oneD_vol_operator, this->oneD_vol_operator);
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);
1646  }
1647  return dim_FR_operator;
1648 
1649 }
1650 
1651 
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,
1657  const Parameters::AllParameters::Flux_Reconstruction_Aux FR_param_aux_input)
1658  : local_Flux_Reconstruction_operator<dim,n_faces,real>::local_Flux_Reconstruction_operator(nstate_input, max_degree_input, grid_degree_input, Parameters::AllParameters::Flux_Reconstruction::cDG)
1659  , FR_param_aux_type(FR_param_aux_input)
1660 {
1661  //Initialize to the max degrees
1662  current_degree = max_degree_input;
1663  //get the FR corrcetion parameter value
1665 }
1666 
1667 template <int dim, int n_faces, typename real>
1669  const unsigned int curr_cell_degree,
1670  double &k)
1671 {
1673  if(FR_param_aux_type == FR_Aux_enum::kHU){
1674  this->get_Huynh_g2_parameter(curr_cell_degree, k);
1675  }
1676  else if(FR_param_aux_type == FR_Aux_enum::kSD){
1677  this->get_spectral_difference_parameter(curr_cell_degree, k);
1678  }
1679  else if(FR_param_aux_type == FR_Aux_enum::kNegative){
1680  this->get_c_negative_FR_parameter(curr_cell_degree, k);
1681  }
1682  else if(FR_param_aux_type == FR_Aux_enum::kNegative2){//knegative divided by 2
1683  this->get_c_negative_divided_by_two_FR_parameter(curr_cell_degree, k);
1684  }
1685  else if(FR_param_aux_type == FR_Aux_enum::kDG){
1686  k = 0.0;
1687  }
1688  else if(FR_param_aux_type == FR_Aux_enum::k10Thousand){
1689  k = 10000.0;
1690  }
1691  else if(FR_param_aux_type == FR_Aux_enum::kPlus){
1692  this->get_c_plus_parameter(curr_cell_degree, k);
1693  }
1694 }
1695 template <int dim, int n_faces, typename real>
1697  const dealii::FESystem<1,1> &finite_element,
1698  const dealii::Quadrature<1> &quadrature)
1699 {
1700  const unsigned int n_dofs = finite_element.dofs_per_cell;
1701  //build the FR correction operator
1702  derivative_p<dim,n_faces,real> pth_derivative(this->nstate, this->max_degree, this->max_grid_degree);
1703  pth_derivative.build_1D_volume_operator(finite_element, quadrature);
1704  local_mass<dim,n_faces,real> local_Mass_Matrix(this->nstate, this->max_degree, this->max_grid_degree);
1705  local_Mass_Matrix.build_1D_volume_operator(finite_element, quadrature);
1706  //allocate the volume operator
1707  this->oneD_vol_operator.reinit(n_dofs, n_dofs);
1708  //solves
1709  this->build_local_Flux_Reconstruction_operator(local_Mass_Matrix.oneD_vol_operator, pth_derivative.oneD_vol_operator, n_dofs, FR_param_aux, this->oneD_vol_operator);
1710 }
1711 
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)
1717  : SumFactorizedOperators<dim,n_faces,real>::SumFactorizedOperators(nstate_input, max_degree_input, grid_degree_input)
1718 {
1719  //Initialize to the max degrees
1720  current_degree = max_degree_input;
1721 }
1722 
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)
1728 {
1729  norm_matrix_inverse.mTmult(volume_projection, integral_vol_basis);
1730 }
1731 template <int dim, int n_faces, typename real>
1733  const dealii::FESystem<1,1> &finite_element,
1734  const dealii::Quadrature<1> &quadrature)
1735 {
1736  const unsigned int n_dofs = finite_element.dofs_per_cell;
1737  const unsigned int n_quad_pts = quadrature.size();
1738  vol_integral_basis<dim,n_faces,real> integral_vol_basis(this->nstate, this->max_degree, this->max_grid_degree);
1739  integral_vol_basis.build_1D_volume_operator(finite_element, quadrature);
1740  local_mass<dim,n_faces,real> local_Mass_Matrix(this->nstate, this->max_degree, this->max_grid_degree);
1741  local_Mass_Matrix.build_1D_volume_operator(finite_element, quadrature);
1742  dealii::FullMatrix<double> mass_inv(n_dofs);
1743  mass_inv.invert(local_Mass_Matrix.oneD_vol_operator);
1744  //allocate the volume operator
1745  this->oneD_vol_operator.reinit(n_dofs, n_quad_pts);
1746  //solves
1747  compute_local_vol_projection_operator(mass_inv, integral_vol_basis.oneD_vol_operator, this->oneD_vol_operator);
1748 }
1749 
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)
1757  : vol_projection_operator<dim,n_faces,real>::vol_projection_operator(nstate_input, max_degree_input, grid_degree_input)
1758  , store_transpose(store_transpose_input)
1759  , FR_param_type(FR_param_input)
1760 {
1761  //Initialize to the max degrees
1762  current_degree = max_degree_input;
1763 }
1764 
1765 template <int dim, int n_faces, typename real>
1767  const dealii::FESystem<1,1> &finite_element,
1768  const dealii::Quadrature<1> &quadrature)
1769 {
1770  const unsigned int n_dofs = finite_element.dofs_per_cell;
1771  const unsigned int n_quad_pts = quadrature.size();
1772  vol_integral_basis<dim,n_faces,real> integral_vol_basis(this->nstate, this->max_degree, this->max_grid_degree);
1773  integral_vol_basis.build_1D_volume_operator(finite_element, quadrature);
1774  FR_mass_inv<dim,n_faces,real> local_FR_Mass_Matrix_inv(this->nstate, this->max_degree, this->max_grid_degree, FR_param_type);
1775  local_FR_Mass_Matrix_inv.build_1D_volume_operator(finite_element, quadrature);
1776  //allocate the volume operator
1777  this->oneD_vol_operator.reinit(n_dofs, n_quad_pts);
1778  //solves
1779  this->compute_local_vol_projection_operator(local_FR_Mass_Matrix_inv.oneD_vol_operator, integral_vol_basis.oneD_vol_operator, this->oneD_vol_operator);
1780 
1781  if(store_transpose){
1782  oneD_transpose_vol_operator.reinit(n_quad_pts, n_dofs);
1783  for(unsigned int idof=0; idof<n_dofs; idof++){
1784  for(unsigned int iquad=0; iquad<n_quad_pts; iquad++){
1785  oneD_transpose_vol_operator[iquad][idof] = this->oneD_vol_operator[idof][iquad];
1786  }
1787  }
1788  }
1789 }
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)
1797  : vol_projection_operator<dim,n_faces,real>::vol_projection_operator(nstate_input, max_degree_input, grid_degree_input)
1798  , store_transpose(store_transpose_input)
1799  , FR_param_type(FR_param_input)
1800 {
1801  //Initialize to the max degrees
1802  current_degree = max_degree_input;
1803 }
1804 
1805 template <int dim, int n_faces, typename real>
1807  const dealii::FESystem<1,1> &finite_element,
1808  const dealii::Quadrature<1> &quadrature)
1809 {
1810  const unsigned int n_dofs = finite_element.dofs_per_cell;
1811  const unsigned int n_quad_pts = quadrature.size();
1812  vol_integral_basis<dim,n_faces,real> integral_vol_basis(this->nstate, this->max_degree, this->max_grid_degree);
1813  integral_vol_basis.build_1D_volume_operator(finite_element, quadrature);
1814  FR_mass_inv_aux<dim,n_faces,real> local_FR_Mass_Matrix_inv(this->nstate, this->max_degree, this->max_grid_degree, FR_param_type);
1815  local_FR_Mass_Matrix_inv.build_1D_volume_operator(finite_element, quadrature);
1816  //allocate the volume operator
1817  this->oneD_vol_operator.reinit(n_dofs, n_quad_pts);
1818  //solves
1819  this->compute_local_vol_projection_operator(local_FR_Mass_Matrix_inv.oneD_vol_operator, integral_vol_basis.oneD_vol_operator, this->oneD_vol_operator);
1820 
1821  if(store_transpose){
1822  oneD_transpose_vol_operator.reinit(n_quad_pts, n_dofs);
1823  for(unsigned int idof=0; idof<n_dofs; idof++){
1824  for(unsigned int iquad=0; iquad<n_quad_pts; iquad++){
1825  oneD_transpose_vol_operator[iquad][idof] = this->oneD_vol_operator[idof][iquad];
1826  }
1827  }
1828  }
1829 }
1830 
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,
1837  : SumFactorizedOperators<dim,n_faces,real>::SumFactorizedOperators(nstate_input, max_degree_input, grid_degree_input)
1838  , FR_param_type(FR_param_input)
1839 {
1840  //Initialize to the max degrees
1841  current_degree = max_degree_input;
1842 }
1843 
1844 template <int dim, int n_faces, typename real>
1846  const dealii::FESystem<1,1> &finite_element,
1847  const dealii::Quadrature<1> &quadrature)
1848 {
1849  const unsigned int n_dofs = finite_element.dofs_per_cell;
1850  local_mass<dim,n_faces,real> local_Mass_Matrix(this->nstate, this->max_degree, this->max_grid_degree);
1851  local_Mass_Matrix.build_1D_volume_operator(finite_element, quadrature);
1853  local_FR_oper.build_1D_volume_operator(finite_element, quadrature);
1854  dealii::FullMatrix<double> FR_mass_matrix(n_dofs);
1855  FR_mass_matrix.add(1.0, local_Mass_Matrix.oneD_vol_operator, 1.0, local_FR_oper.oneD_vol_operator);
1856  //allocate the volume operator
1857  this->oneD_vol_operator.reinit(n_dofs, n_dofs);
1858  //solves
1859  this->oneD_vol_operator.invert(FR_mass_matrix);
1860 }
1861 
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,
1868  : SumFactorizedOperators<dim,n_faces,real>::SumFactorizedOperators(nstate_input, max_degree_input, grid_degree_input)
1869  , FR_param_type(FR_param_input)
1870 {
1871  //Initialize to the max degrees
1872  current_degree = max_degree_input;
1873 }
1874 
1875 template <int dim, int n_faces, typename real>
1877  const dealii::FESystem<1,1> &finite_element,
1878  const dealii::Quadrature<1> &quadrature)
1879 {
1880  const unsigned int n_dofs = finite_element.dofs_per_cell;
1881  local_mass<dim,n_faces,real> local_Mass_Matrix(this->nstate, this->max_degree, this->max_grid_degree);
1882  local_Mass_Matrix.build_1D_volume_operator(finite_element, quadrature);
1884  local_FR_oper.build_1D_volume_operator(finite_element, quadrature);
1885  dealii::FullMatrix<double> FR_mass_matrix(n_dofs);
1886  FR_mass_matrix.add(1.0, local_Mass_Matrix.oneD_vol_operator, 1.0, local_FR_oper.oneD_vol_operator);
1887  //allocate the volume operator
1888  this->oneD_vol_operator.reinit(n_dofs, n_dofs);
1889  //solves
1890  this->oneD_vol_operator.invert(FR_mass_matrix);
1891 }
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,
1898  : SumFactorizedOperators<dim,n_faces,real>::SumFactorizedOperators(nstate_input, max_degree_input, grid_degree_input)
1899  , FR_param_type(FR_param_input)
1900 {
1901  //Initialize to the max degrees
1902  current_degree = max_degree_input;
1903 }
1904 
1905 template <int dim, int n_faces, typename real>
1907  const dealii::FESystem<1,1> &finite_element,
1908  const dealii::Quadrature<1> &quadrature)
1909 {
1910  const unsigned int n_dofs = finite_element.dofs_per_cell;
1911  local_mass<dim,n_faces,real> local_Mass_Matrix(this->nstate, this->max_degree, this->max_grid_degree);
1912  local_Mass_Matrix.build_1D_volume_operator(finite_element, quadrature);
1914  local_FR_oper.build_1D_volume_operator(finite_element, quadrature);
1915  dealii::FullMatrix<double> FR_mass_matrix(n_dofs);
1916  FR_mass_matrix.add(1.0, local_Mass_Matrix.oneD_vol_operator, 1.0, local_FR_oper.oneD_vol_operator);
1917  //allocate the volume operator
1918  this->oneD_vol_operator.reinit(n_dofs, n_dofs);
1919  //solves
1920  this->oneD_vol_operator.add(1.0, FR_mass_matrix);
1921 }
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,
1928  : SumFactorizedOperators<dim,n_faces,real>::SumFactorizedOperators(nstate_input, max_degree_input, grid_degree_input)
1929  , FR_param_type(FR_param_input)
1930 {
1931  //Initialize to the max degrees
1932  current_degree = max_degree_input;
1933 }
1934 
1935 template <int dim, int n_faces, typename real>
1937  const dealii::FESystem<1,1> &finite_element,
1938  const dealii::Quadrature<1> &quadrature)
1939 {
1940  const unsigned int n_dofs = finite_element.dofs_per_cell;
1941  local_mass<dim,n_faces,real> local_Mass_Matrix(this->nstate, this->max_degree, this->max_grid_degree);
1942  local_Mass_Matrix.build_1D_volume_operator(finite_element, quadrature);
1944  local_FR_oper.build_1D_volume_operator(finite_element, quadrature);
1945  dealii::FullMatrix<double> FR_mass_matrix(n_dofs);
1946  FR_mass_matrix.add(1.0, local_Mass_Matrix.oneD_vol_operator, 1.0, local_FR_oper.oneD_vol_operator);
1947  //allocate the volume operator
1948  this->oneD_vol_operator.reinit(n_dofs, n_dofs);
1949  //solves
1950  this->oneD_vol_operator.add(1.0, FR_mass_matrix);
1951 }
1952 
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)
1958  : SumFactorizedOperators<dim,n_faces,real>::SumFactorizedOperators(nstate_input, max_degree_input, grid_degree_input)
1959 {
1960  //Initialize to the max degrees
1961  current_degree = max_degree_input;
1962 }
1963 
1964 template <int dim, int n_faces, typename real>
1966  const dealii::FESystem<1,1> &finite_element,
1967  const dealii::Quadrature<1> &quadrature)
1968 {
1969  const unsigned int n_quad_pts = quadrature.size();
1970  const unsigned int n_dofs = finite_element.dofs_per_cell;
1971  //allocate
1972  this->oneD_grad_operator.reinit(n_quad_pts, n_dofs);
1973  //loop and store
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];
1981  }
1982  }
1983 }
1984 
1985 /*************************************
1986 *
1987 * SURFACE OPERATORS
1988 *
1989 *************************************/
1990 
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)
1996  : SumFactorizedOperators<dim,n_faces,real>::SumFactorizedOperators(nstate_input, max_degree_input, grid_degree_input)
1997 {
1998  //Initialize to the max degrees
1999  current_degree = max_degree_input;
2000 }
2001 
2002 template <int dim, int n_faces, typename real>
2004  const dealii::FESystem<1,1> &finite_element,
2005  const dealii::Quadrature<0> &face_quadrature)
2006 {
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 ();
2011  //loop and store
2012  for(unsigned int iface=0; iface<n_faces_1D; iface++){
2013  //allocate the facet operator
2014  this->oneD_surf_operator[iface].reinit(n_face_quad_pts, n_dofs);
2015  //sum factorized operators use a 1D element.
2016  const dealii::Quadrature<1> quadrature = dealii::QProjector<1>::project_to_face(dealii::ReferenceCell::get_hypercube(1),
2017  face_quadrature,
2018  iface);
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;
2023  //Basis function idof of poly degree idegree evaluated at cubature node qpoint.
2024  this->oneD_surf_operator[iface][iquad][idof] = finite_element.shape_value_component(idof,qpoint,istate)
2025  * quad_weights[iquad];
2026  }
2027  }
2028  }
2029 }
2030 
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)
2036  : SumFactorizedOperators<dim,n_faces,real>::SumFactorizedOperators(nstate_input, max_degree_input, grid_degree_input)
2037 {
2038  //Initialize to the max degrees
2039  current_degree = max_degree_input;
2040 }
2041 
2042 template <int dim, int n_faces, typename real>
2044  const dealii::FESystem<1,1> &finite_element,
2045  const dealii::Quadrature<1> &quadrature)
2046 {
2047  const unsigned int n_dofs = finite_element.dofs_per_cell;
2048  local_mass<dim,n_faces,real> local_Mass_Matrix(this->nstate, this->max_degree, this->max_grid_degree);
2049  local_Mass_Matrix.build_1D_volume_operator(finite_element, quadrature);
2050  //allocate the volume operator
2051  this->oneD_vol_operator.reinit(n_dofs, n_dofs);
2052  //solves
2053  this->oneD_vol_operator.add(1.0, local_Mass_Matrix.oneD_vol_operator);
2054 }
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)
2061 {
2062  dealii::FullMatrix<double> norm_inv(n_dofs);
2063  norm_inv.invert(norm_matrix);
2064  norm_inv.mTmult(lifting, face_integral);
2065 }
2066 template <int dim, int n_faces, typename real>
2068  const dealii::FESystem<1,1> &finite_element,
2069  const dealii::Quadrature<0> &face_quadrature)
2070 {
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;
2074  //create surface integral of basis functions
2075  face_integral_basis<dim,n_faces,real> basis_int_facet(this->nstate, this->max_degree, this->max_grid_degree);
2076  basis_int_facet.build_1D_surface_operator(finite_element, face_quadrature);
2077  //loop and store
2078  for(unsigned int iface=0; iface<n_faces_1D; iface++){
2079  //allocate the facet operator
2080  this->oneD_surf_operator[iface].reinit(n_dofs, n_face_quad_pts);
2081  build_local_surface_lifting_operator(n_dofs, this->oneD_vol_operator, basis_int_facet.oneD_surf_operator[iface], this->oneD_surf_operator[iface]);
2082  }
2083 }
2084 
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,
2091  : lifting_operator<dim,n_faces,real>::lifting_operator(nstate_input, max_degree_input, grid_degree_input)
2092  , FR_param_type(FR_param_input)
2093 {
2094  //Initialize to the max degrees
2095  current_degree = max_degree_input;
2096 }
2097 
2098 template <int dim, int n_faces, typename real>
2100  const dealii::FESystem<1,1> &finite_element,
2101  const dealii::Quadrature<1> &quadrature)
2102 {
2103  const unsigned int n_dofs = finite_element.dofs_per_cell;
2104  local_mass<dim,n_faces,real> local_Mass_Matrix(this->nstate, this->max_degree, this->max_grid_degree);
2105  local_Mass_Matrix.build_1D_volume_operator(finite_element, quadrature);
2107  local_FR.build_1D_volume_operator(finite_element, quadrature);
2108  //allocate the volume operator
2109  this->oneD_vol_operator.reinit(n_dofs, n_dofs);
2110  //solves
2111  this->oneD_vol_operator.add(1.0, local_Mass_Matrix.oneD_vol_operator);
2112  this->oneD_vol_operator.add(1.0, local_FR.oneD_vol_operator);
2113 }
2114 template <int dim, int n_faces, typename real>
2116  const dealii::FESystem<1,1> &finite_element,
2117  const dealii::Quadrature<0> &face_quadrature)
2118 {
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;
2122  //create surface integral of basis functions
2123  face_integral_basis<dim,n_faces,real> basis_int_facet(this->nstate, this->max_degree, this->max_grid_degree);
2124  basis_int_facet.build_1D_surface_operator(finite_element, face_quadrature);
2125  //loop and store
2126  for(unsigned int iface=0; iface<n_faces_1D; iface++){
2127  //allocate the facet operator
2128  this->oneD_surf_operator[iface].reinit(n_dofs, n_face_quad_pts);
2129  this->build_local_surface_lifting_operator(n_dofs, this->oneD_vol_operator, basis_int_facet.oneD_surf_operator[iface], this->oneD_surf_operator[iface]);
2130  }
2131 }
2132 
2133 
2134 /******************************************************************************
2135 *
2136 * METRIC MAPPING OPERATORS
2137 *
2138 ******************************************************************************/
2139 
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)
2145  : SumFactorizedOperators<dim,n_faces,real>::SumFactorizedOperators(nstate_input, max_degree_input, 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)
2148 {
2149  //Initialize to the max degrees
2150  current_degree = max_degree_input;
2151  current_grid_degree = grid_degree_input;
2152 }
2153 
2154 template <int dim, int n_faces, typename real>
2156  const dealii::FESystem<1,1> &finite_element,
2157  const dealii::Quadrature<1> &quadrature)
2158 {
2159  assert(finite_element.dofs_per_cell == quadrature.size());//checks collocation
2160  mapping_shape_functions_grid_nodes.build_1D_volume_operator(finite_element, quadrature);
2161  mapping_shape_functions_grid_nodes.build_1D_gradient_operator(finite_element, quadrature);
2162 }
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)
2168 {
2169  mapping_shape_functions_flux_nodes.build_1D_volume_operator(finite_element, quadrature);
2170  mapping_shape_functions_flux_nodes.build_1D_gradient_operator(finite_element, quadrature);
2171  mapping_shape_functions_flux_nodes.build_1D_surface_operator(finite_element, face_quadrature);
2172  mapping_shape_functions_flux_nodes.build_1D_surface_gradient_operator(finite_element, face_quadrature);
2173 
2174 }
2175 template <int dim, int n_faces, typename real>
2177  const dealii::FESystem<1,1> &finite_element,
2178  const dealii::Quadrature<1> &quadrature)
2179 {
2180  mapping_shape_functions_flux_nodes.build_1D_volume_operator(finite_element, quadrature);
2181  mapping_shape_functions_flux_nodes.build_1D_gradient_operator(finite_element, quadrature);
2182 }
2183 
2184 /***********************************************
2185 *
2186 * METRIC DET JAC AND COFACTOR
2187 *
2188 ************************************************/
2189 //Constructor
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)
2198  : SumFactorizedOperators<dim,n_faces,real>::SumFactorizedOperators(nstate_input, max_degree_input, grid_degree_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)
2202 {}
2203 
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)
2209 {
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];
2213  }
2214  }
2215 
2216 }
2217 
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)
2223 {
2224  for(int idim=0; idim<dim; idim++){
2225  phys[idim] = metric_cofactor[idim] * ref;
2226 // phys[idim] = 0.0;
2227 // for(int idim2=0; idim2<dim; idim2++){
2228 // phys[idim] += metric_cofactor[idim][idim2]
2229 // * ref[idim2];
2230 // }
2231  }
2232 }
2233 
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)
2239 {
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];
2247  }
2248  }
2249  }
2250 
2251 }
2252 
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)
2259 {
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]
2265  * ref[idim2];
2266  }
2267  }
2268  phys[iquad] /= phys[iquad].norm();
2269  }
2270 }
2271 
2272 template <typename real, int dim, int n_faces>
2274  const unsigned int n_quad_pts,
2275  const unsigned int /*n_metric_dofs*/,//dofs of metric basis. NOTE: this is the number of mapping support points
2276  const std::array<std::vector<real>,dim> &mapping_support_points,
2278 {
2279  det_Jac_vol.resize(n_quad_pts);
2280  //compute determinant of metric Jacobian
2282  n_quad_pts,
2283  mapping_support_points,
2284  mapping_basis.mapping_shape_functions_flux_nodes.oneD_vol_operator,
2285  mapping_basis.mapping_shape_functions_flux_nodes.oneD_vol_operator,
2286  mapping_basis.mapping_shape_functions_flux_nodes.oneD_vol_operator,
2287  mapping_basis.mapping_shape_functions_flux_nodes.oneD_grad_operator,
2288  mapping_basis.mapping_shape_functions_flux_nodes.oneD_grad_operator,
2289  mapping_basis.mapping_shape_functions_flux_nodes.oneD_grad_operator,
2290  det_Jac_vol);
2291 }
2292 
2293 template <typename real, int dim, int n_faces>
2295  const unsigned int n_quad_pts,
2296  const unsigned int n_metric_dofs,//dofs of metric basis. NOTE: this is the number of mapping support points
2297  const std::array<std::vector<real>,dim> &mapping_support_points,
2299  const bool use_invariant_curl_form)
2300 {
2301  det_Jac_vol.resize(n_quad_pts);
2302  for(int idim=0; idim<dim; idim++){
2303  for(int jdim=0; jdim<dim; jdim++){
2304  metric_cofactor_vol[idim][jdim].resize(n_quad_pts);
2305  }
2306  }
2307  //compute determinant of metric Jacobian
2309  n_quad_pts,
2310  mapping_support_points,
2311  mapping_basis.mapping_shape_functions_flux_nodes.oneD_vol_operator,
2312  mapping_basis.mapping_shape_functions_flux_nodes.oneD_vol_operator,
2313  mapping_basis.mapping_shape_functions_flux_nodes.oneD_vol_operator,
2314  mapping_basis.mapping_shape_functions_flux_nodes.oneD_grad_operator,
2315  mapping_basis.mapping_shape_functions_flux_nodes.oneD_grad_operator,
2316  mapping_basis.mapping_shape_functions_flux_nodes.oneD_grad_operator,
2317  det_Jac_vol);
2318  //compute the metric cofactor
2320  n_quad_pts,
2321  n_metric_dofs,
2322  mapping_support_points,
2323  mapping_basis.mapping_shape_functions_grid_nodes.oneD_vol_operator,
2324  mapping_basis.mapping_shape_functions_grid_nodes.oneD_vol_operator,
2325  mapping_basis.mapping_shape_functions_grid_nodes.oneD_vol_operator,
2326  mapping_basis.mapping_shape_functions_flux_nodes.oneD_vol_operator,
2327  mapping_basis.mapping_shape_functions_flux_nodes.oneD_vol_operator,
2328  mapping_basis.mapping_shape_functions_flux_nodes.oneD_vol_operator,
2329  mapping_basis.mapping_shape_functions_grid_nodes.oneD_grad_operator,
2330  mapping_basis.mapping_shape_functions_grid_nodes.oneD_grad_operator,
2331  mapping_basis.mapping_shape_functions_grid_nodes.oneD_grad_operator,
2332  mapping_basis.mapping_shape_functions_flux_nodes.oneD_grad_operator,
2333  mapping_basis.mapping_shape_functions_flux_nodes.oneD_grad_operator,
2334  mapping_basis.mapping_shape_functions_flux_nodes.oneD_grad_operator,
2336  use_invariant_curl_form);
2337 
2339  for(int idim=0; idim<dim; idim++){
2340  flux_nodes_vol[idim].resize(n_quad_pts);
2341  this->matrix_vector_mult(mapping_support_points[idim],
2342  flux_nodes_vol[idim],
2343  mapping_basis.mapping_shape_functions_flux_nodes.oneD_vol_operator,
2344  mapping_basis.mapping_shape_functions_flux_nodes.oneD_vol_operator,
2345  mapping_basis.mapping_shape_functions_flux_nodes.oneD_vol_operator);
2346  }
2347  }
2348 }
2349 
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,//dofs of metric basis. NOTE: this is the number of mapping support points
2355  const std::array<std::vector<real>,dim> &mapping_support_points,
2357  const bool use_invariant_curl_form)
2358 {
2359  det_Jac_surf.resize(n_quad_pts);
2360  for(int idim=0; idim<dim; idim++){
2361  for(int jdim=0; jdim<dim; jdim++){
2362  metric_cofactor_surf[idim][jdim].resize(n_quad_pts);
2363  }
2364  }
2365  //compute determinant of metric Jacobian
2367  n_quad_pts,
2368  mapping_support_points,
2369  (iface == 0) ? mapping_basis.mapping_shape_functions_flux_nodes.oneD_surf_operator[0] :
2370  ((iface == 1) ? mapping_basis.mapping_shape_functions_flux_nodes.oneD_surf_operator[1] :
2371  mapping_basis.mapping_shape_functions_flux_nodes.oneD_vol_operator),
2372  (iface == 2) ? mapping_basis.mapping_shape_functions_flux_nodes.oneD_surf_operator[0] :
2373  ((iface == 3) ? mapping_basis.mapping_shape_functions_flux_nodes.oneD_surf_operator[1] :
2374  mapping_basis.mapping_shape_functions_flux_nodes.oneD_vol_operator),
2375  (iface == 4) ? mapping_basis.mapping_shape_functions_flux_nodes.oneD_surf_operator[0] :
2376  ((iface == 5) ? mapping_basis.mapping_shape_functions_flux_nodes.oneD_surf_operator[1] :
2377  mapping_basis.mapping_shape_functions_flux_nodes.oneD_vol_operator),
2378  (iface == 0) ? mapping_basis.mapping_shape_functions_flux_nodes.oneD_surf_grad_operator[0] :
2379  ((iface == 1) ? mapping_basis.mapping_shape_functions_flux_nodes.oneD_surf_grad_operator[1] :
2380  mapping_basis.mapping_shape_functions_flux_nodes.oneD_grad_operator),
2381  (iface == 2) ? mapping_basis.mapping_shape_functions_flux_nodes.oneD_surf_grad_operator[0] :
2382  ((iface == 3) ? mapping_basis.mapping_shape_functions_flux_nodes.oneD_surf_grad_operator[1] :
2383  mapping_basis.mapping_shape_functions_flux_nodes.oneD_grad_operator),
2384  (iface == 4) ? mapping_basis.mapping_shape_functions_flux_nodes.oneD_surf_grad_operator[0] :
2385  ((iface == 5) ? mapping_basis.mapping_shape_functions_flux_nodes.oneD_surf_grad_operator[1] :
2386  mapping_basis.mapping_shape_functions_flux_nodes.oneD_grad_operator),
2387  det_Jac_surf);
2388  //compute the metric cofactor
2390  n_quad_pts,
2391  n_metric_dofs,
2392  mapping_support_points,
2393  mapping_basis.mapping_shape_functions_grid_nodes.oneD_vol_operator,
2394  mapping_basis.mapping_shape_functions_grid_nodes.oneD_vol_operator,
2395  mapping_basis.mapping_shape_functions_grid_nodes.oneD_vol_operator,
2396  (iface == 0) ? mapping_basis.mapping_shape_functions_flux_nodes.oneD_surf_operator[0] :
2397  ((iface == 1) ? mapping_basis.mapping_shape_functions_flux_nodes.oneD_surf_operator[1] :
2398  mapping_basis.mapping_shape_functions_flux_nodes.oneD_vol_operator),
2399  (iface == 2) ? mapping_basis.mapping_shape_functions_flux_nodes.oneD_surf_operator[0] :
2400  ((iface == 3) ? mapping_basis.mapping_shape_functions_flux_nodes.oneD_surf_operator[1] :
2401  mapping_basis.mapping_shape_functions_flux_nodes.oneD_vol_operator),
2402  (iface == 4) ? mapping_basis.mapping_shape_functions_flux_nodes.oneD_surf_operator[0] :
2403  ((iface == 5) ? mapping_basis.mapping_shape_functions_flux_nodes.oneD_surf_operator[1] :
2404  mapping_basis.mapping_shape_functions_flux_nodes.oneD_vol_operator),
2405  mapping_basis.mapping_shape_functions_grid_nodes.oneD_grad_operator,
2406  mapping_basis.mapping_shape_functions_grid_nodes.oneD_grad_operator,
2407  mapping_basis.mapping_shape_functions_grid_nodes.oneD_grad_operator,
2408  (iface == 0) ? mapping_basis.mapping_shape_functions_flux_nodes.oneD_surf_grad_operator[0] :
2409  ((iface == 1) ? mapping_basis.mapping_shape_functions_flux_nodes.oneD_surf_grad_operator[1] :
2410  mapping_basis.mapping_shape_functions_flux_nodes.oneD_grad_operator),
2411  (iface == 2) ? mapping_basis.mapping_shape_functions_flux_nodes.oneD_surf_grad_operator[0] :
2412  ((iface == 3) ? mapping_basis.mapping_shape_functions_flux_nodes.oneD_surf_grad_operator[1] :
2413  mapping_basis.mapping_shape_functions_flux_nodes.oneD_grad_operator),
2414  (iface == 4) ? mapping_basis.mapping_shape_functions_flux_nodes.oneD_surf_grad_operator[0] :
2415  ((iface == 5) ? mapping_basis.mapping_shape_functions_flux_nodes.oneD_surf_grad_operator[1] :
2416  mapping_basis.mapping_shape_functions_flux_nodes.oneD_grad_operator),
2418  use_invariant_curl_form);
2419 
2421  for(int iface=0; iface<n_faces; iface++){
2422  for(int idim=0; idim<dim; idim++){
2423  flux_nodes_surf[iface][idim].resize(n_quad_pts);
2424  this->matrix_vector_mult(mapping_support_points[idim],
2425  flux_nodes_surf[iface][idim],
2426  (iface == 0) ? mapping_basis.mapping_shape_functions_flux_nodes.oneD_surf_operator[0] :
2427  ((iface == 1) ? mapping_basis.mapping_shape_functions_flux_nodes.oneD_surf_operator[1] :
2428  mapping_basis.mapping_shape_functions_flux_nodes.oneD_vol_operator),
2429  (iface == 2) ? mapping_basis.mapping_shape_functions_flux_nodes.oneD_surf_operator[0] :
2430  ((iface == 3) ? mapping_basis.mapping_shape_functions_flux_nodes.oneD_surf_operator[1] :
2431  mapping_basis.mapping_shape_functions_flux_nodes.oneD_vol_operator),
2432  (iface == 4) ? mapping_basis.mapping_shape_functions_flux_nodes.oneD_surf_operator[0] :
2433  ((iface == 5) ? mapping_basis.mapping_shape_functions_flux_nodes.oneD_surf_operator[1] :
2434  mapping_basis.mapping_shape_functions_flux_nodes.oneD_vol_operator));
2435  }
2436  }
2437  }
2438 }
2439 
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)
2451 {
2452  for(int idim=0; idim<dim; idim++){
2453  for(int jdim=0; jdim<dim; jdim++){
2454  // std::vector<real> output_vect(pow(n_quad_pts_1D,dim));
2455  std::vector<real> output_vect(n_quad_pts);
2456  if(jdim == 0)
2457  this->matrix_vector_mult(mapping_support_points[idim], output_vect,
2458  grad_basis_x_flux_nodes,
2459  basis_y_flux_nodes,
2460  basis_z_flux_nodes);
2461  if(jdim == 1)
2462  this->matrix_vector_mult(mapping_support_points[idim], output_vect,
2463  basis_x_flux_nodes,
2464  grad_basis_y_flux_nodes,
2465  basis_z_flux_nodes);
2466  if(jdim == 2)
2467  this->matrix_vector_mult(mapping_support_points[idim], output_vect,
2468  basis_x_flux_nodes,
2469  basis_y_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];
2473  }
2474  }
2475  }
2476 }
2477 
2478 template <typename real, int dim, int n_faces>
2480  const unsigned int n_quad_pts,//number volume 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)
2489 {
2490  //mapping support points must be passed as a vector[dim][n_metric_dofs]
2491  assert(pow(this->max_grid_degree+1,dim) == mapping_support_points[0].size());
2492 
2493  std::vector<dealii::Tensor<2,dim,double>> Jacobian_flux_nodes(n_quad_pts);
2494  this->build_metric_Jacobian(n_quad_pts,
2495  mapping_support_points,
2496  basis_x_flux_nodes,
2497  basis_y_flux_nodes,
2498  basis_z_flux_nodes,
2499  grad_basis_x_flux_nodes,
2500  grad_basis_y_flux_nodes,
2501  grad_basis_z_flux_nodes,
2502  Jacobian_flux_nodes);
2503  if(store_Jacobian){
2504  for(int idim=0; idim<dim; idim++){
2505  for(int jdim=0; jdim<dim; jdim++){
2506  metric_Jacobian_vol_cubature[idim][jdim].resize(n_quad_pts);
2507  for(unsigned int iquad=0; iquad<n_quad_pts; iquad++){
2508  metric_Jacobian_vol_cubature[idim][jdim][iquad] = Jacobian_flux_nodes[iquad][idim][jdim];
2509  }
2510  }
2511  }
2512  }
2513 
2514  for(unsigned int iquad=0; iquad<n_quad_pts; iquad++){
2515  det_metric_Jac[iquad] = dealii::determinant(Jacobian_flux_nodes[iquad]);
2516  //check for valid cell
2517  if(det_metric_Jac[iquad] <= 1e-14){
2518  std::cout<<"The determinant of the Jacobian is negative. Aborting..."<<std::endl;
2519  std::abort();
2520  }
2521  }
2522 }
2523 template <typename real, int dim, int n_faces>
2525  const unsigned int n_quad_pts,//number flux pts
2526  const unsigned int n_metric_dofs,//dofs of metric basis. NOTE: this is the number of mapping support points
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)
2542 {
2543  //mapping support points must be passed as a vector[dim][n_metric_dofs]
2544  //Solve for Cofactor
2545  if (dim == 1){//constant for 1D
2546  std::fill(metric_cofactor[0][0].begin(), metric_cofactor[0][0].end(), 1.0);
2547  }
2548  if (dim == 2){
2549  std::vector<dealii::Tensor<2,dim,double>> Jacobian_flux_nodes(n_quad_pts);
2550  this->build_metric_Jacobian(n_quad_pts,
2551  mapping_support_points,
2552  basis_x_flux_nodes,
2553  basis_y_flux_nodes,
2554  basis_z_flux_nodes,
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];
2564  }
2565  }
2566  if (dim == 3){
2567  compute_local_3D_cofactor(n_metric_dofs,
2568  n_quad_pts,
2569  mapping_support_points,
2570  basis_x_grid_nodes,
2571  basis_y_grid_nodes,
2572  basis_z_grid_nodes,
2573  basis_x_flux_nodes,
2574  basis_y_flux_nodes,
2575  basis_z_flux_nodes,
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,
2582  metric_cofactor,
2583  use_invariant_curl_form);
2584  }
2585 }
2586 
2587 template <typename real, int dim, int n_faces>
2589  const unsigned int n_metric_dofs,
2590  const unsigned int /*n_quad_pts*/,
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)
2606 {
2607  //Conservative Curl Form
2608 
2609  //compute grad Xm, NOTE: it is the same as the Jacobian at Grid Nodes
2610  std::vector<dealii::Tensor<2,dim,real>> grad_Xm(n_metric_dofs);//gradient of mapping support points at Grid Nodes
2611  this->build_metric_Jacobian(n_metric_dofs,
2612  mapping_support_points,
2613  basis_x_grid_nodes,
2614  basis_y_grid_nodes,
2615  basis_z_grid_nodes,
2616  grad_basis_x_grid_nodes,
2617  grad_basis_y_grid_nodes,
2618  grad_basis_z_grid_nodes,
2619  grad_Xm);
2620 
2621  //Store Xl * grad Xm using fact the mapping shape functions collocated on Grid Nodes.
2622  //Also, we only store the ones needed, not all to reduce computational cost.
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);
2626 
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);
2630 
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);
2634 
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];
2640  }
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];
2645  }
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];
2650  }
2651 
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];
2656  }
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];
2661  }
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];
2666  }
2667 
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];
2672  }
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];
2677  }
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];
2682  }
2683  }
2684 
2685  //Compute metric Cofactor via conservative curl form at flux nodes.
2686  //C11
2687  this->matrix_vector_mult(z_dy_dzeta, metric_cofactor[0][0],
2688  basis_x_flux_nodes, grad_basis_y_flux_nodes, basis_z_flux_nodes, false, -1.0);
2689  this->matrix_vector_mult(z_dy_deta, metric_cofactor[0][0],
2690  basis_x_flux_nodes, basis_y_flux_nodes, grad_basis_z_flux_nodes, true);
2691  //C12
2692  this->matrix_vector_mult(z_dy_dzeta, metric_cofactor[0][1],
2693  grad_basis_x_flux_nodes, basis_y_flux_nodes, basis_z_flux_nodes);
2694  this->matrix_vector_mult(z_dy_dxi, metric_cofactor[0][1],
2695  basis_x_flux_nodes, basis_y_flux_nodes, grad_basis_z_flux_nodes, true, -1.0);
2696  //C13
2697  this->matrix_vector_mult(z_dy_deta, metric_cofactor[0][2],
2698  grad_basis_x_flux_nodes, basis_y_flux_nodes, basis_z_flux_nodes, false, -1.0);
2699  this->matrix_vector_mult(z_dy_dxi, metric_cofactor[0][2],
2700  basis_x_flux_nodes, grad_basis_y_flux_nodes, basis_z_flux_nodes, true);
2701 
2702  //C21
2703  this->matrix_vector_mult(x_dz_dzeta, metric_cofactor[1][0],
2704  basis_x_flux_nodes, grad_basis_y_flux_nodes, basis_z_flux_nodes, false, -1.0);
2705  this->matrix_vector_mult(x_dz_deta, metric_cofactor[1][0],
2706  basis_x_flux_nodes, basis_y_flux_nodes, grad_basis_z_flux_nodes, true);
2707  //C22
2708  this->matrix_vector_mult(x_dz_dzeta, metric_cofactor[1][1],
2709  grad_basis_x_flux_nodes, basis_y_flux_nodes, basis_z_flux_nodes);
2710  this->matrix_vector_mult(x_dz_dxi, metric_cofactor[1][1],
2711  basis_x_flux_nodes, basis_y_flux_nodes, grad_basis_z_flux_nodes, true, -1.0);
2712  //C23
2713  this->matrix_vector_mult(x_dz_deta, metric_cofactor[1][2],
2714  grad_basis_x_flux_nodes, basis_y_flux_nodes, basis_z_flux_nodes, false, -1.0);
2715  this->matrix_vector_mult(x_dz_dxi, metric_cofactor[1][2],
2716  basis_x_flux_nodes, grad_basis_y_flux_nodes, basis_z_flux_nodes, true);
2717 
2718  //C31
2719  this->matrix_vector_mult(y_dx_dzeta, metric_cofactor[2][0],
2720  basis_x_flux_nodes, grad_basis_y_flux_nodes, basis_z_flux_nodes, false, -1.0);
2721  this->matrix_vector_mult(y_dx_deta, metric_cofactor[2][0],
2722  basis_x_flux_nodes, basis_y_flux_nodes, grad_basis_z_flux_nodes, true);
2723  //C32
2724  this->matrix_vector_mult(y_dx_dzeta, metric_cofactor[2][1],
2725  grad_basis_x_flux_nodes, basis_y_flux_nodes, basis_z_flux_nodes);
2726  this->matrix_vector_mult(y_dx_dxi, metric_cofactor[2][1],
2727  basis_x_flux_nodes, basis_y_flux_nodes, grad_basis_z_flux_nodes, true, -1.0);
2728  //C33
2729  this->matrix_vector_mult(y_dx_deta, metric_cofactor[2][2],
2730  grad_basis_x_flux_nodes, basis_y_flux_nodes, basis_z_flux_nodes, false, -1.0);
2731  this->matrix_vector_mult(y_dx_dxi, metric_cofactor[2][2],
2732  basis_x_flux_nodes, grad_basis_y_flux_nodes, basis_z_flux_nodes, true);
2733 
2734 }
2735 
2736 /**********************************
2737 *
2738 * Sum Factorization STATE class
2739 *
2740 **********************************/
2741 //Constructor
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)
2746  : SumFactorizedOperators<dim,n_faces,real>::SumFactorizedOperators(nstate, max_degree_input, grid_degree_input)
2747 {}
2748 
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)
2753  : SumFactorizedOperatorsState<dim,nstate,n_faces,real>::SumFactorizedOperatorsState(max_degree_input, grid_degree_input)
2754 {
2755  //Initialize to the max degrees
2756  current_degree = max_degree_input;
2757 }
2758 
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)
2763 {
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;
2767  //Note thate the flux basis should only have one state in the finite element.
2768  //loop and store
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;
2774  if(ishape == 0)
2775  this->oneD_vol_state_operator[istate].reinit(n_quad_pts, n_shape_fns);
2776 
2777  //Basis function idof of poly degree idegree evaluated at cubature node qpoint.
2778  this->oneD_vol_state_operator[istate][iquad][ishape] = finite_element.shape_value_component(idof,qpoint,istate);
2779  }
2780  }
2781 }
2782 
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)
2787 {
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;
2791  //loop and store
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;
2797  if(ishape == 0)
2798  this->oneD_grad_state_operator[istate].reinit(n_quad_pts, n_shape_fns);
2799 
2800  this->oneD_grad_state_operator[istate][iquad][ishape] = finite_element.shape_grad_component(idof, qpoint, istate)[0];
2801  }
2802  }
2803 }
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)
2808 {
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;
2813  //loop and store
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),
2816  face_quadrature,
2817  iface);
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;
2822  if(ishape == 0)
2823  this->oneD_surf_state_operator[istate][iface].reinit(n_face_quad_pts, n_shape_fns);
2824 
2825  this->oneD_surf_state_operator[istate][iface][iquad][ishape] = finite_element.shape_value_component(idof,quadrature.point(iquad),istate);
2826  }
2827  }
2828  }
2829 }
2830 
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)
2835  : SumFactorizedOperatorsState<dim,nstate,n_faces,real>::SumFactorizedOperatorsState(max_degree_input, grid_degree_input)
2836 {
2837  //Initialize to the max degrees
2838  current_degree = max_degree_input;
2839 }
2840 
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)
2845 {
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);
2849  //Note thate the flux basis should only have one state in the finite element.
2850  //loop and store
2851  for(int istate=0; istate<nstate; istate++){
2852  //allocate the basis at volume cubature
2853  this->oneD_vol_state_operator[istate].reinit(n_quad_pts, n_dofs);
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++){
2857  //Basis function idof of poly degree idegree evaluated at cubature node qpoint.
2858  this->oneD_vol_state_operator[istate][iquad][idof] = finite_element.shape_value_component(idof,qpoint,0);
2859  }
2860  }
2861  }
2862 }
2863 
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)
2868 {
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);
2872  //loop and store
2873  for(int istate=0; istate<nstate; istate++){
2874  //allocate
2875  this->oneD_grad_state_operator[istate].reinit(n_quad_pts, n_dofs);
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++){
2879  this->oneD_grad_state_operator[istate][iquad][idof] = finite_element.shape_grad_component(idof, qpoint, 0)[0];
2880  }
2881  }
2882  }
2883 }
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)
2888 {
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;
2892  //loop and store
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),
2895  face_quadrature,
2896  iface);
2897  for(int istate=0; istate<nstate; istate++){
2898  this->oneD_surf_state_operator[istate][iface].reinit(n_face_quad_pts, n_dofs);
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);
2902  }
2903  }
2904  }
2905  }
2906 }
2907 
2908 
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)
2913  : flux_basis_functions_state<dim,nstate,n_faces,real>::flux_basis_functions_state(max_degree_input, grid_degree_input)
2914 {
2915  //Initialize to the max degrees
2916  current_degree = max_degree_input;
2917 }
2918 
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)
2923 {
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;
2927  //loop and store
2928  const std::vector<double> &quad_weights = quadrature.get_weights ();
2929  for(int istate_flux=0; istate_flux<nstate; istate_flux++){
2930  //allocate
2931  this->oneD_vol_state_operator[istate_flux].reinit(n_dofs, n_quad_pts);
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++){
2935  double value = 0.0;
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]
2940  * this->oneD_grad_state_operator[istate_flux][iquad][idof];
2941  }
2942  this->oneD_vol_state_operator[istate_flux][itest][idof] = value;
2943  }
2944  }
2945  }
2946 }
2947 
2948 
2950 
2952 
2958 
2975 
2976 //template class basis_at_facet_cubature <PHILIP_DIM, 2*PHILIP_DIM>;
2980 
2982 
2984 //template class vol_metric_operators <double,PHILIP_DIM, 2*PHILIP_DIM>;
2985 //template class vol_determinant_metric_Jacobian<PHILIP_DIM, 2*PHILIP_DIM>;
2986 //template class vol_metric_cofactor<PHILIP_DIM, 2*PHILIP_DIM>;
2987 //template class surface_metric_cofactor<PHILIP_DIM, 2*PHILIP_DIM>;
2988 //
2994 
3005 
3006 } // OPERATOR namespace
3007 } // PHiLiP namespace
3008 
The FLUX basis functions separated by nstate with n shape functions.
Definition: operators.h:1348
unsigned int current_degree
Stores the degree of the current poly degree.
Definition: operators.h:732
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...
Definition: operators.cpp:602
unsigned int current_degree
Stores the degree of the current poly degree.
Definition: operators.h:1325
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.
Definition: operators.cpp:1253
basis_functions< dim, n_faces, real > mapping_shape_functions_grid_nodes
Object of mapping shape functions evaluated at grid nodes.
Definition: operators.h:1042
dealii::Tensor< 2, dim, std::vector< real > > metric_cofactor_vol
The volume metric cofactor matrix.
Definition: operators.h:1165
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.
Definition: operators.cpp:2273
void build_1D_volume_state_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
Definition: operators.cpp:2760
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.
Definition: operators.cpp:2056
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
Definition: operators.cpp:1082
const Parameters::AllParameters::Flux_Reconstruction FR_param_type
Flux reconstruction parameter type.
Definition: operators.h:738
std::array< dealii::FullMatrix< double >, 2 > oneD_surf_operator
Stores the one dimensional surface operator.
Definition: operators.h:360
The metric independent inverse of the FR mass matrix .
Definition: operators.h:782
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.
Definition: operators.cpp:536
bool store_transpose
Flag is store transpose operator.
Definition: operators.h:735
unsigned int current_degree
Stores the degree of the current poly degree.
Definition: operators.h:426
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.
Definition: operators.cpp:2191
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
Definition: operators.cpp:1766
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.
Definition: operators.cpp:1668
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...
Definition: operators.cpp:932
The integration of gradient of solution basis.
Definition: operators.h:883
unsigned int current_degree
Stores the degree of the current poly degree.
Definition: operators.h:816
unsigned int current_degree
Stores the degree of the current poly degree.
Definition: operators.h:550
void build_1D_gradient_state_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
Definition: operators.cpp:2865
Sum Factorization derived class.
Definition: operators.h:131
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.
Definition: operators.cpp:2086
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...
Definition: operators.cpp:196
unsigned int current_grid_degree
Stores the degree of the current grid degree.
Definition: operators.h:1039
const unsigned int max_grid_degree
Max grid degree.
Definition: operators.h:68
unsigned int current_degree
Stores the degree of the current poly degree.
Definition: operators.h:863
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
Definition: operators.cpp:1356
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.
Definition: operators.cpp:344
const Parameters::AllParameters::Flux_Reconstruction_Aux FR_param_type
Flux reconstruction parameter type.
Definition: operators.h:769
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
Definition: operators.cpp:1845
dealii::FullMatrix< double > oneD_grad_operator
Stores the one dimensional gradient operator.
Definition: operators.h:363
dealii::ConditionalOStream pcout
Parallel std::cout that only outputs on mpi_rank==0.
Definition: operators.h:122
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.
Definition: operators.cpp:1408
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.
Definition: operators.cpp:55
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
Definition: operators.cpp:1732
The metric independent FR mass matrix for auxiliary equation .
Definition: operators.h:852
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.
Definition: operators.cpp:1923
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.
Definition: operators.cpp:106
bool store_transpose
Flag is store transpose operator.
Definition: operators.h:766
const Parameters::AllParameters::Flux_Reconstruction_Aux FR_param_type
Flux reconstruction parameter type.
Definition: operators.h:819
const int nstate
Number of states.
Definition: operators.h:70
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.
Definition: operators.cpp:2254
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.
Definition: operators.cpp:795
double compute_factorial(double n)
Standard function to compute factorial of a number.
Definition: operators.cpp:173
Projection operator corresponding to basis functions onto -norm for auxiliary equation.
Definition: operators.h:751
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.
Definition: operators.cpp:524
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
Definition: operators.cpp:1538
unsigned int current_degree
Stores the degree of the current poly degree.
Definition: operators.h:893
Files for the baseline physics.
Definition: ADTypes.hpp:10
void build_1D_surface_state_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 0 > &face_quadrature)
Assembles the one dimensional operator.
Definition: operators.cpp:2885
double FR_param_aux
Flux reconstruction paramater value.
Definition: operators.h:672
lifting_operator(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input)
Constructor.
Definition: operators.cpp:2032
unsigned int current_degree
Stores the degree of the current poly degree.
Definition: operators.h:1036
void build_1D_surface_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 0 > &face_quadrature)
Assembles the one dimensional operator.
Definition: operators.cpp:2067
const unsigned int max_degree
Max polynomial degree.
Definition: operators.h:66
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.
Definition: operators.cpp:2588
SumFactorizedOperatorsState(const unsigned int max_degree_input, const unsigned int grid_degree_input)
Constructor.
Definition: operators.cpp:2743
unsigned int current_degree
Stores the degree of the current poly degree.
Definition: operators.h:926
std::array< dealii::FullMatrix< double >, 2 > oneD_surf_grad_operator
Stores the one dimensional surface gradient operator.
Definition: operators.h:366
dealii::Tensor< 2, dim, std::vector< real > > metric_cofactor_surf
The facet metric cofactor matrix, for ONE face.
Definition: operators.h:1168
dealii::FullMatrix< double > oneD_transpose_vol_operator
Stores the transpose of the operator for fast weight-adjusted solves.
Definition: operators.h:777
unsigned int current_degree
Stores the degree of the current poly degree.
Definition: operators.h:993
In order to have all state operators be arrays of array, we template by dim, type, nstate, and number of faces.
Definition: operators.h:1292
unsigned int current_degree
Stores the degree of the current poly degree.
Definition: operators.h:951
basis_functions_state(const unsigned int max_degree_input, const unsigned int grid_degree_input)
Constructor.
Definition: operators.cpp:2750
Flux_Reconstruction
Type of correction in Flux Reconstruction.
dealii::FullMatrix< double > oneD_skew_symm_vol_oper
Skew-symmetric volume operator .
Definition: operators.h:494
const bool store_Jacobian
Flag if store metric Jacobian at flux nodes.
Definition: operators.h:1099
void get_Huynh_g2_parameter(const unsigned int curr_cell_degree, double &c)
Evaluates Huynh&#39;s g2 parameter for flux reconstruction.
Definition: operators.cpp:1423
-th order modal derivative of basis fuctions, ie/
Definition: operators.h:519
face_integral_basis(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input)
Constructor.
Definition: operators.cpp:1992
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.
Definition: operators.cpp:2155
That is Quadrature Weights multiplies with basis_at_vol_cubature.
Definition: operators.h:416
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.
Definition: operators.cpp:1653
const bool store_surf_flux_nodes
Flag if store metric Jacobian at flux nodes.
Definition: operators.h:1105
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.
Definition: operators.cpp:2524
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.
Definition: operators.cpp:2219
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.
Definition: operators.cpp:426
This is the solution basis , the modal differential opertaor commonly seen in DG defined as ...
Definition: operators.h:499
ESFR correction matrix without jac dependence.
Definition: operators.h:539
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.
Definition: operators.cpp:1863
std::vector< real > det_Jac_vol
The determinant of the metric Jacobian at volume cubature nodes.
Definition: operators.h:1171
Local mass matrix without jacobian dependence.
Definition: operators.h:436
vol_projection_operator(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input)
Constructor.
Definition: operators.cpp:1713
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.
Definition: operators.cpp:1893
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.
Definition: operators.cpp:1751
The DG lifting operator is defined as the operator that lifts inner products of polynomials of some o...
Definition: operators.h:941
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
Definition: operators.cpp:1696
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...
Definition: operators.cpp:1456
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.
Definition: operators.cpp:2043
unsigned int current_degree
Stores the degree of the current poly degree.
Definition: operators.h:391
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.
Definition: operators.cpp:2441
unsigned int current_degree
Stores the degree of the current poly degree.
Definition: operators.h:763
unsigned int current_degree
Stores the degree of the current poly degree.
Definition: operators.h:1392
virtual void build_1D_volume_state_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
Definition: operators.cpp:2842
dealii::FullMatrix< double > oneD_transpose_vol_operator
Stores the transpose of the operator for fast weight-adjusted solves.
Definition: operators.h:746
dealii::Tensor< 2, dim, std::vector< real > > metric_Jacobian_vol_cubature
Stores the metric Jacobian at flux nodes.
Definition: operators.h:1177
void get_spectral_difference_parameter(const unsigned int curr_cell_degree, double &c)
Evaluates the spectral difference parameter for flux reconstruction.
Definition: operators.cpp:1434
void build_1D_gradient_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
Definition: operators.cpp:1102
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
Definition: operators.cpp:1187
The ESFR lifting operator.
Definition: operators.h:982
const Parameters::AllParameters::Flux_Reconstruction_Aux FR_param_type
Flux reconstruction parameter type.
Definition: operators.h:866
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...
Definition: operators.cpp:1522
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...
Definition: operators.cpp:1000
double FR_param
Flux reconstruction paramater value.
Definition: operators.h:556
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.
Definition: operators.cpp:682
unsigned int current_degree
Stores the degree of the current poly degree.
Definition: operators.h:666
void build_1D_gradient_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
Definition: operators.cpp:1965
void build_1D_surface_state_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 0 > &face_quadrature)
Assembles the one dimensional operator.
Definition: operators.cpp:2805
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
Definition: operators.cpp:1219
Base metric operators class that stores functions used in both the volume and on surface.
Definition: operators.h:1086
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.
Definition: operators.cpp:319
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.
Definition: operators.cpp:2479
const Parameters::AllParameters::Flux_Reconstruction FR_param_type
Flux reconstruction parameter type.
Definition: operators.h:842
ESFR correction matrix for AUX EQUATION without jac dependence.
Definition: operators.h:655
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...
Definition: operators.cpp:369
void build_1D_gradient_state_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
Definition: operators.cpp:2784
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.
Definition: operators.cpp:1292
The mapping shape functions evaluated at the desired nodes (facet set included in volume grid nodes f...
Definition: operators.h:1026
SumFactorizedOperators(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input)
Precompute 1D operator in constructor.
Definition: operators.cpp:188
modal_basis_differential_operator(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input)
Constructor.
Definition: operators.cpp:1345
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.
Definition: operators.cpp:2099
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. ...
Definition: operators.cpp:381
void get_FR_correction_parameter(const unsigned int curr_cell_degree, double &c)
Gets the FR correction parameter for the primary equation and stores.
Definition: operators.cpp:1492
const bool store_vol_flux_nodes
Flag if store metric Jacobian at flux nodes.
Definition: operators.h:1102
OperatorsBase(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input)
Constructor.
Definition: operators.cpp:42
The basis functions separated by nstate with n shape functions.
Definition: operators.h:1316
basis_functions(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input)
Constructor.
Definition: operators.cpp:1071
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...
Definition: operators.cpp:414
std::array< dealii::FullMatrix< double >, nstate > oneD_vol_state_operator
Stores the one dimensional volume operator.
Definition: operators.h:1301
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.
Definition: operators.cpp:2176
The metric independent FR mass matrix .
Definition: operators.h:828
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...
Definition: operators.cpp:460
const Parameters::AllParameters::Flux_Reconstruction FR_param_type
Flux reconstruction parameter type.
Definition: operators.h:553
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.
Definition: operators.cpp:2205
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
Definition: operators.cpp:1385
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.
Definition: operators.cpp:1791
dealii::FullMatrix< double > oneD_vol_operator
Stores the one dimensional volume operator.
Definition: operators.h:355
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...
Definition: operators.cpp:879
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.
Definition: operators.cpp:2351
unsigned int current_degree
Stores the degree of the current poly degree.
Definition: operators.h:793
void build_1D_surface_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 0 > &quadrature)
Assembles the one dimensional operator.
Definition: operators.cpp:1122
derivative_p(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input)
Constructor.
Definition: operators.cpp:1374
std::vector< real > det_Jac_surf
The determinant of the metric Jacobian at facet cubature nodes.
Definition: operators.h:1174
basis_functions< dim, n_faces, real > mapping_shape_functions_flux_nodes
Object of mapping shape functions evaluated at flux nodes.
Definition: operators.h:1045
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...
Definition: operators.cpp:1445
unsigned int current_degree
Stores the degree of the current poly degree.
Definition: operators.h:704
const Parameters::AllParameters::Flux_Reconstruction_Aux FR_param_aux_type
Flux reconstruction parameter type.
Definition: operators.h:669
The metric independent inverse of the FR mass matrix for auxiliary equation .
Definition: operators.h:805
void get_c_plus_parameter(const unsigned int curr_cell_degree, double &c)
Gets the FR correction parameter corresponding to the maximum timestep.
Definition: operators.cpp:1464
std::array< std::array< dealii::FullMatrix< double >, 2 >, nstate > oneD_surf_state_operator
Stores the one dimensional surface operator.
Definition: operators.h:1304
Projection operator corresponding to basis functions onto -norm.
Definition: operators.h:720
void build_1D_surface_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 0 > &face_quadrature)
Assembles the one dimensional operator.
Definition: operators.cpp:2115
unsigned int current_degree
Stores the degree of the current poly degree.
Definition: operators.h:446
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)...
Definition: operators.cpp:1556
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.
Definition: operators.cpp:2294
void build_1D_surface_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 0 > &face_quadrature)
Assembles the one dimensional operator.
Definition: operators.cpp:2003
const bool store_skew_symmetric_form
Flag to store the skew symmetric form .
Definition: operators.h:486
std::array< dealii::Tensor< 1, dim, std::vector< real > >, n_faces > flux_nodes_surf
Stores the physical facet flux nodes.
Definition: operators.h:1183
"Stiffness" operator used in DG Strong form.
Definition: operators.h:1383
const Parameters::AllParameters::Flux_Reconstruction FR_param_type
Flux reconstruction parameter type.
Definition: operators.h:996
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).
Definition: operators.cpp:1724
unsigned int current_degree
Stores the degree of the current poly degree.
Definition: operators.h:839
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
Definition: operators.cpp:1305
std::array< dealii::FullMatrix< double >, nstate > oneD_grad_state_operator
Stores the one dimensional gradient operator.
Definition: operators.h:1307
mapping_shape_functions(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input)
Constructor.
Definition: operators.cpp:2141
local_flux_basis_stiffness(const unsigned int max_degree_input, const unsigned int grid_degree_input)
Constructor.
Definition: operators.cpp:2910
Operator base class.
Definition: operators.h:53
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...
Definition: operators.cpp:1617
unsigned int current_degree
Stores the degree of the current poly degree.
Definition: operators.h:509
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
Definition: operators.cpp:1806
vol_integral_basis(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input)
Constructor.
Definition: operators.cpp:1176
const Parameters::AllParameters::Flux_Reconstruction FR_param_type
Flux reconstruction parameter type.
Definition: operators.h:796
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
Definition: operators.cpp:1906
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.
Definition: operators.cpp:2164
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.
Definition: operators.cpp:2235
vol_integral_gradient_basis(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input)
Constructor.
Definition: operators.cpp:1954
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.
Definition: operators.cpp:308
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.
Definition: operators.cpp:1832
void build_1D_surface_gradient_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 0 > &quadrature)
Assembles the one dimensional operator.
Definition: operators.cpp:1149
unsigned int current_degree
Stores the degree of the current poly degree.
Definition: operators.h:1357
dealii::Tensor< 1, dim, std::vector< real > > flux_nodes_vol
Stores the physical volume flux nodes.
Definition: operators.h:1180
The surface integral of test functions.
Definition: operators.h:916
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
Definition: operators.cpp:1936
unsigned int current_degree
Stores the degree of the current poly degree.
Definition: operators.h:529
Projection operator corresponding to basis functions onto M-norm (L2).
Definition: operators.h:694
Local stiffness matrix without jacobian dependence.
Definition: operators.h:472
flux_basis_functions_state(const unsigned int max_degree_input, const unsigned int grid_degree_input)
Constructor.
Definition: operators.cpp:2832
void build_1D_volume_state_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
Definition: operators.cpp:2920
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...
Definition: operators.cpp:814
unsigned int current_degree
Stores the degree of the current poly degree.
Definition: operators.h:483
local_mass(const int nstate_input, const unsigned int max_degree_input, const unsigned int grid_degree_input)
Constructor.
Definition: operators.cpp:1208
void build_1D_volume_operator(const dealii::FESystem< 1, 1 > &finite_element, const dealii::Quadrature< 1 > &quadrature)
Assembles the one dimensional operator.
Definition: operators.cpp:1876