165 * #ifndef MFMFE_DATA_H
166 * #define MFMFE_DATA_H
168 * #include <deal.II/base/
function.h>
169 * #include <deal.II/base/tensor_function.h>
174 * <a name=
"data.h-Dataandexactsolution"></a>
175 * <h3>Data and exact solution.</h3>
179 * This file declares the classes
for the given data, i.e.
180 * right-hand side, exact solution, permeability tensor and
181 * boundary conditions. For simplicity only 2
d cases are
182 * provided, but 3
d can be added straightforwardly.
193 *
class RightHandSide :
public Function<dim>
196 * RightHandSide () :
Function<dim>(1) {}
199 *
const unsigned int component = 0)
const override;
203 *
double RightHandSide<dim>::value (
const Point<dim> &p,
204 *
const unsigned int )
const
206 *
const double x =
p[0];
207 *
const double y =
p[1];
212 *
return -(x*(y*y*y*y)*6.0-(y*y)*
sin(x*y*2.0)*2.0+2.0)*(x*2.0+x*x+y*y+1.0)-
sin(x*y)*(
cos(x*y*2.0)+(x*x)*(y*y*y)*1.2E1
213 * -x*y*
sin(x*y*2.0)*2.0)*2.0-(x*2.0+2.0)*(x*2.0+(x*x)*(y*y*y*y)*3.0+y*
cos(x*y*2.0))+(x*x)*(
sin(x*y*2.0)
214 * -x*(y*y)*6.0)*
pow(x+1.0,2.0)*2.0-x*
cos(x*y)*(x*2.0+(x*x)*(y*y*y*y)*3.0+y*(
pow(
cos(x*y),2.0)*2.0-1.0))
215 * -x*y*
cos(x*y)*((x*x)*(y*y*y)*4.0+
pow(
cos(x*y),2.0)*2.0-1.0);
217 *
Assert(
false, ExcMessage(
"The RHS data for dim != 2 is not provided"));
224 *
class PressureBoundaryValues :
public Function<dim>
227 * PressureBoundaryValues () :
Function<dim>(1) {}
230 *
const unsigned int component = 0)
const override;
234 *
double PressureBoundaryValues<dim>::value (
const Point<dim> &p,
235 *
const unsigned int )
const
237 *
const double x =
p[0];
238 *
const double y =
p[1];
243 *
return (x*x*x)*(y*y*y*y)+
cos(x*y)*
sin(x*y)+x*x;
245 *
Assert(
false, ExcMessage(
"The BC data for dim != 2 is not provided"));
252 *
class ExactSolution :
public Function<dim>
255 * ExactSolution () :
Function<dim>(dim+1) {}
266 * ExactSolution<dim>::vector_value (
const Point<dim> &p,
270 * ExcDimensionMismatch (
values.size(), dim+1));
272 *
const double x =
p[0];
273 *
const double y =
p[1];
278 *
values(0) = -(x*2.0+(x*x)*(y*y*y*y)*3.0+y*
cos(x*y*2.0))*(x*2.0+x*x+y*y+1.0)-x*
sin(x*y)*(
cos(x*y*2.0)+(x*x)*(y*y*y)*4.0);
279 *
values(1) = -
sin(x*y)*(x*2.0+(x*x)*(y*y*y*y)*3.0+y*
cos(x*y*2.0))-x*(
cos(x*y*2.0)+(x*x)*(y*y*y)*4.0)*
pow(x+1.0,2.0);
283 *
Assert(
false, ExcMessage(
"The exact solution for dim != 2 is not provided"));
289 * ExactSolution<dim>::vector_gradient (
const Point<dim> &p,
292 *
const double x =
p[0];
293 *
const double y =
p[1];
298 * grads[0][0] = -(x*(y*y*y*y)*6.0-(y*y)*
sin(x*y*2.0)*2.0+2.0)*(x*2.0+x*x+y*y+1.0)-
sin(x*y)*(
cos(x*y*2.0)
299 * +(x*x)*(y*y*y)*1.2E1-x*y*
sin(x*y*2.0)*2.0)-(x*2.0+2.0)*(x*2.0+(x*x)*(y*y*y*y)*3.0
300 * +y*
cos(x*y*2.0))-x*y*
cos(x*y)*((x*x)*(y*y*y)*4.0+
pow(
cos(x*y),2.0)*2.0-1.0);
301 * grads[0][1] = -(
cos(x*y*2.0)+(x*x)*(y*y*y)*1.2E1-x*y*
sin(x*y*2.0)*2.0)*(x*2.0+x*x+y*y+1.0)
302 * -y*(x*2.0+(x*x)*(y*y*y*y)*3.0+y*
cos(x*y*2.0))*2.0-(x*x)*
cos(x*y)*((x*x)*(y*y*y)*4.0
303 * +
pow(
cos(x*y),2.0)*2.0-1.0)+(x*x)*
sin(x*y)*(
sin(x*y*2.0)-x*(y*y)*6.0)*2.0;
304 * grads[1][0] = -
sin(x*y)*(x*(y*y*y*y)*6.0-(y*y)*
sin(x*y*2.0)*2.0+2.0)-
pow(x+1.0,2.0)*(
cos(x*y*2.0)
305 * +(x*x)*(y*y*y)*1.2E1-x*y*
sin(x*y*2.0)*2.0)-x*(
cos(x*y*2.0)+(x*x)*(y*y*y)*4.0)*(x*2.0+2.0)
306 * -y*
cos(x*y)*(x*2.0+(x*x)*(y*y*y*y)*3.0+y*(
pow(
cos(x*y),2.0)*2.0-1.0));
307 * grads[1][1] = -
sin(x*y)*(
cos(x*y*2.0)+(x*x)*(y*y*y)*1.2E1-x*y*
sin(x*y*2.0)*2.0)+(x*x)*(
sin(x*y*2.0)
308 * -x*(y*y)*6.0)*
pow(x+1.0,2.0)*2.0-x*
cos(x*y)*(x*2.0+(x*x)*(y*y*y*y)*3.0
309 * +y*(
pow(
cos(x*y),2.0)*2.0-1.0));
313 *
Assert(
false, ExcMessage(
"The exact solution's gradient for dim != 2 is not provided"));
331 * KInverse<dim>::value_list (
const std::vector<
Point<dim> > &points,
335 * ExcDimensionMismatch (points.size(),
values.size()));
337 *
for (
unsigned int p=0;
p<points.size(); ++
p)
341 *
const double x = points[
p][0];
342 *
const double y = points[
p][1];
347 *
values[
p][0][0] =
pow(x+1.0,2.0)/(x*4.0+(x*x)*(y*y)-
pow(
sin(x*y),2.0)+x*(y*y)*2.0+(x*x)*6.0+(x*x*x)*4.0+x*x*x*x+y*y+1.0);
348 *
values[
p][0][1] = -
sin(x*y)/(x*4.0+(x*x)*(y*y)-
pow(
sin(x*y),2.0)+x*(y*y)*2.0+(x*x)*6.0+(x*x*x)*4.0+x*x*x*x+y*y+1.0);
349 *
values[
p][1][0] = -
sin(x*y)/(x*4.0+(x*x)*(y*y)-
pow(
sin(x*y),2.0)+x*(y*y)*2.0+(x*x)*6.0+(x*x*x)*4.0+x*x*x*x+y*y+1.0);
350 *
values[
p][1][1] = (x*2.0+x*x+y*y+1.0)/(x*4.0+(x*x)*(y*y)-
pow(
sin(x*y),2.0)+x*(y*y)*2.0+(x*x)*6.0+(x*x*x)*4.0+x*x*x*x+y*y+1.0);
353 *
Assert(
false, ExcMessage(
"The inverse of permeability tensor for dim != 2 is not provided"));
363<a name=
"ann-mfmfe.cc"></a>
364<h1>Annotated version of mfmfe.cc</h1>
387 * <a name=
"mfmfe.cc-Includefiles"></a>
388 * <h3>Include files</h3>
392 * As usual, the list of necessary header files. There is not
393 * much
new here, the files are included in order
394 * base-lac-grid-dofs-numerics followed by the
C++ headers.
397 * #include <deal.II/base/convergence_table.h>
398 * #include <deal.II/base/quadrature_lib.h>
399 * #include <deal.II/base/logstream.h>
400 * #include <deal.II/base/timer.h>
401 * #include <deal.II/base/utilities.h>
402 * #include <deal.II/base/work_stream.h>
404 * #include <deal.II/lac/full_matrix.h>
405 * #include <deal.II/lac/solver_cg.h>
406 * #include <deal.II/lac/block_sparse_matrix.h>
407 * #include <deal.II/lac/block_vector.h>
408 * #include <deal.II/lac/precondition.h>
410 * #include <deal.II/grid/grid_generator.h>
411 * #include <deal.II/grid/grid_tools.h>
412 * #include <deal.II/grid/grid_in.h>
413 * #include <deal.II/grid/
tria.h>
414 * #include <deal.II/dofs/dof_renumbering.h>
415 * #include <deal.II/dofs/dof_tools.h>
416 * #include <deal.II/fe/fe_dgq.h>
417 * #include <deal.II/fe/fe_system.h>
418 * #include <deal.II/fe/fe_tools.h>
419 * #include <deal.II/numerics/vector_tools.h>
420 * #include <deal.II/numerics/matrix_tools.h>
421 * #include <deal.II/numerics/data_out.h>
424 * #include <unordered_map>
428 * This is a header needed
for the purposes of the
429 * multipoint flux mixed method, as
it declares the
430 *
new enhanced Raviart-Thomas finite element.
433 * #include <deal.II/fe/fe_rt_bubbles.h>
437 * For the sake of readability, the classes representing
438 * data, i.e. RHS, BCs, permeability tensor and the exact
439 * solution are placed in a file data.h which is included
447 * As
always the program is in the
namespace of its own with
448 * the deal.II classes and
functions imported into
it
458 * <a name=
"mfmfe.cc-Definitionofmultipointfluxassemblydatastructures"></a>
459 * <h3>Definition of multipoint flux assembly data structures</h3>
463 * The main idea of the MFMFE method is to perform local elimination
464 * of the velocity variables in order to obtain the resulting
465 * pressure system. Since in deal.II assembly happens cell-wise,
466 * some extra work needs to be done in order to get the local
467 * mass matrices @f$A_i@f$ and the corresponding to them @f$B_i@f$.
470 *
namespace DataStructures
474 * This will be achieved by assembling cell-wise, but instead of placing
475 * the terms into a global system
matrix, they will populate node-associated
476 * full matrices. For
this, a data structure with fast lookup is crucial, hence
483 *
size_t operator()(
const Point<dim> &p)
const
486 * h1 = std::hash<double>()(
p[0]);
493 * h2 = std::hash<double>()(
p[1]);
496 * h2 = std::hash<double>()(
p[1]);
497 * h3 = std::hash<double>()(
p[2]);
498 *
return (h1 ^ (h2 << 1)) ^ h3;
500 *
Assert(
false, ExcNotImplemented());
507 * Here, the actual hash-tables are defined. We use the
C++ STL <code>unordered_map</code>,
508 * with the hash
function specified above. For convenience these are aliased as follows
512 *
using PointToMatrixMap = std::unordered_map<Point<dim>, std::map<std::pair<types::global_dof_index,types::global_dof_index>,
double>, hash_points<dim>>;
515 *
using PointToVectorMap = std::unordered_map<Point<dim>, std::map<types::global_dof_index, double>, hash_points<dim>>;
518 *
using PointToIndexMap = std::unordered_map<Point<dim>, std::set<types::global_dof_index>, hash_points<dim>>;
522 * Next, since
this particular program allows
for the use of
523 * multiple threads, the helper CopyData structures
524 * are defined. There are two kinds of these, one is
used
525 *
for the copying cell-wise contributions to the corresponging
526 * node-associated data structures...
530 *
struct NodeAssemblyCopyData
532 * PointToMatrixMap<dim> cell_mat;
533 * PointToVectorMap<dim> cell_vec;
534 * PointToIndexMap<dim> local_pres_indices;
535 * PointToIndexMap<dim> local_vel_indices;
536 * std::vector<types::global_dof_index> local_dof_indices;
541 * ... and the other one
for the actual process of
542 * local velocity elimination and assembling the global
547 *
struct NodeEliminationCopyData
560 * Similarly, two ScratchData classes are defined.
561 * One
for the assembly part, where we need
563 *
for the basis fuctions...
567 *
struct NodeAssemblyScratchData
574 * NodeAssemblyScratchData (
const NodeAssemblyScratchData &scratch_data);
578 * std::vector<unsigned int> n_faces_at_vertex;
580 *
const unsigned long num_cells;
582 * std::vector<Tensor<2,dim>> k_inverse_values;
583 * std::vector<double> rhs_values;
584 * std::vector<double> pres_bc_values;
586 * std::vector<Tensor<1,dim> > phi_u;
587 * std::vector<double> div_phi_u;
588 * std::vector<double> phi_p;
592 * NodeAssemblyScratchData<dim>::
602 * fe_face_values (fe,
607 * k_inverse_values(quad.size()),
608 * rhs_values(quad.size()),
609 * pres_bc_values(f_quad.size()),
610 * phi_u(fe.dofs_per_cell),
611 * div_phi_u(fe.dofs_per_cell),
612 * phi_p(fe.dofs_per_cell)
617 *
for (; face != endf; ++face)
618 *
for (
unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
619 * n_faces_at_vertex[face->vertex_index(v)] += 1;
623 * NodeAssemblyScratchData<dim>::
624 * NodeAssemblyScratchData (
const NodeAssemblyScratchData &scratch_data)
626 * fe_values (scratch_data.fe_values.get_fe(),
627 * scratch_data.fe_values.get_quadrature(),
630 * fe_face_values (scratch_data.fe_face_values.get_fe(),
631 * scratch_data.fe_face_values.get_quadrature(),
634 * n_faces_at_vertex(scratch_data.n_faces_at_vertex),
635 * num_cells(scratch_data.num_cells),
636 * k_inverse_values(scratch_data.k_inverse_values),
637 * rhs_values(scratch_data.rhs_values),
638 * pres_bc_values(scratch_data.pres_bc_values),
639 * phi_u(scratch_data.phi_u),
640 * div_phi_u(scratch_data.div_phi_u),
641 * phi_p(scratch_data.phi_p)
646 * ...and the other, simpler one,
for the velocity elimination and recovery
649 *
struct VertexEliminationScratchData
651 * VertexEliminationScratchData () =
default;
652 * VertexEliminationScratchData (
const VertexEliminationScratchData &scratch_data);
663 * VertexEliminationScratchData::
664 * VertexEliminationScratchData (
const VertexEliminationScratchData &scratch_data)
666 * velocity_matrix(scratch_data.velocity_matrix),
667 * pressure_rhs(scratch_data.pressure_rhs),
668 * local_pressure_solution(scratch_data.local_pressure_solution),
669 * tmp_rhs1(scratch_data.tmp_rhs1),
670 * tmp_rhs2(scratch_data.tmp_rhs2),
671 * tmp_rhs3(scratch_data.tmp_rhs3)
680 * <a name=
"mfmfe.cc-ThecodeMultipointMixedDarcyProblemcodeclasstemplate"></a>
681 * <h3>The <code>MultipointMixedDarcyProblem</code>
class template</h3>
685 * The main
class, besides the constructor and destructor, has only one
public member
686 * <code>
run()</code>, similarly to the tutorial programs. The
private members can
687 * be grouped into the ones that are
used for the cell-wise assembly, vertex elimination,
688 * pressure solve, vertex velocity recovery and postprocessing. Apart from the
689 * MFMFE-specific data structures, the rest of the members should look familiar.
693 *
class MultipointMixedDarcyProblem
696 * MultipointMixedDarcyProblem (
const unsigned int degree);
697 * ~MultipointMixedDarcyProblem ();
698 *
void run (
const unsigned int refine);
701 * DataStructures::NodeAssemblyScratchData<dim> &scratch_data,
702 * DataStructures::NodeAssemblyCopyData<dim> ©_data);
703 *
void copy_cell_to_node(
const DataStructures::NodeAssemblyCopyData<dim> ©_data);
704 *
void node_assembly();
705 *
void make_cell_centered_sp ();
706 *
void nodal_elimination(
const typename DataStructures::PointToMatrixMap<dim>::iterator &n_it,
707 * DataStructures::VertexEliminationScratchData &scratch_data,
708 * DataStructures::NodeEliminationCopyData<dim> ©_data);
709 *
void copy_node_to_system(
const DataStructures::NodeEliminationCopyData<dim> ©_data);
710 *
void pressure_assembly ();
711 *
void solve_pressure ();
712 *
void velocity_assembly (
const typename DataStructures::PointToMatrixMap<dim>::iterator &n_it,
713 * DataStructures::VertexEliminationScratchData &scratch_data,
714 * DataStructures::NodeEliminationCopyData<dim> ©_data);
715 *
void copy_node_velocity_to_global(
const DataStructures::NodeEliminationCopyData<dim> ©_data);
716 *
void velocity_recovery ();
717 *
void reset_data_structures ();
718 *
void compute_errors (
const unsigned int cycle);
719 *
void output_results (
const unsigned int cycle,
const unsigned int refine);
721 *
const unsigned int degree;
731 * std::unordered_map<Point<dim>,
FullMatrix<double>, DataStructures::hash_points<dim>> pressure_matrix;
732 * std::unordered_map<Point<dim>,
FullMatrix<double>, DataStructures::hash_points<dim>> A_inverse;
733 * std::unordered_map<Point<dim>,
Vector<double>, DataStructures::hash_points<dim>> velocity_rhs;
735 * DataStructures::PointToMatrixMap<dim> node_matrix;
736 * DataStructures::PointToVectorMap<dim> node_rhs;
738 * DataStructures::PointToIndexMap<dim> pressure_indices;
739 * DataStructures::PointToIndexMap<dim> velocity_indices;
741 *
unsigned long n_v, n_p;
753 * <a name=
"mfmfe.cc-Constructoranddestructorcodereset_data_structurescode"></a>
754 * <h4>Constructor and destructor, <code>reset_data_structures</code></h4>
758 * In the constructor of
this class, we store the
value that was
759 * passed in concerning the degree of the finite elements we shall use (a
761 * and then construct the vector valued element belonging to the space @f$V_h^k@f$ described
762 * in the introduction. The constructor also takes care of initializing the
763 * computing timer, as
it is of interest
for us how well our method performs.
767 * MultipointMixedDarcyProblem<dim>::MultipointMixedDarcyProblem (
const unsigned int degree)
780 * The destructor clears the <code>dof_handler</code> and
781 * all of the data structures we
used for the method.
785 * MultipointMixedDarcyProblem<dim>::~MultipointMixedDarcyProblem()
787 * reset_data_structures ();
788 * dof_handler.clear();
794 * This method clears all the data that was
used after one refinement
799 *
void MultipointMixedDarcyProblem<dim>::reset_data_structures ()
801 * pressure_indices.clear();
802 * velocity_indices.clear();
803 * velocity_rhs.clear();
805 * pressure_matrix.clear();
806 * node_matrix.clear();
814 * <a name=
"mfmfe.cc-Cellwiseassemblyandcreationofthelocalnodalbaseddatastructures"></a>
815 * <h4>Cell-wise assembly and creation of the local, nodal-based data structures</h4>
819 * First, the
function that copies local cell contributions to the corresponding nodal
820 * matrices and vectors is defined. It places the
values obtained from local cell integration
821 * into the correct place in a
matrix/vector corresponging to a specific node.
825 *
void MultipointMixedDarcyProblem<dim>::copy_cell_to_node(
const DataStructures::NodeAssemblyCopyData<dim> ©_data)
827 *
for (
auto m : copy_data.cell_mat)
832 *
for (
auto p : copy_data.cell_vec.at(m.
first))
835 *
for (
auto p : copy_data.local_pres_indices.at(m.
first))
838 *
for (
auto p : copy_data.local_vel_indices.at(m.
first))
847 * Second, the
function that does the cell assembly is defined. While
it is
848 * similar to the tutorial programs in a way
it uses scrath and
copy data
849 * structures, the need to localize the DOFs leads to several differences.
853 *
void MultipointMixedDarcyProblem<dim>::
855 * DataStructures::NodeAssemblyScratchData<dim> &scratch_data,
856 * DataStructures::NodeAssemblyCopyData<dim> ©_data)
858 * copy_data.cell_mat.clear();
859 * copy_data.cell_vec.clear();
860 * copy_data.local_vel_indices.clear();
861 * copy_data.local_pres_indices.clear();
863 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
864 *
const unsigned int n_q_points = scratch_data.fe_values.get_quadrature().size();
865 *
const unsigned int n_face_q_points = scratch_data.fe_face_values.get_quadrature().size();
867 * copy_data.local_dof_indices.resize(dofs_per_cell);
868 * cell->get_dof_indices (copy_data.local_dof_indices);
870 * scratch_data.fe_values.reinit (cell);
872 *
const KInverse<dim> k_inverse;
873 *
const RightHandSide<dim> rhs;
874 *
const PressureBoundaryValues<dim> pressure_bc;
876 * k_inverse.value_list (scratch_data.fe_values.get_quadrature_points(), scratch_data.k_inverse_values);
877 * rhs.value_list(scratch_data.fe_values.get_quadrature_points(), scratch_data.rhs_values);
883 * std::unordered_map<unsigned int, std::unordered_map<unsigned int, double>> div_map;
887 * One, we need to be able to
assemble the communication between velocity and
888 * pressure variables and put
it on the right place in our
final, local version
889 * of the B
matrix. This is a little messy, as such communication is not in fact
890 * local, so we
do it in two steps. First, we compute all relevant LHS and RHS
893 *
for (
unsigned int q=0; q<n_q_points; ++q)
895 *
const Point<dim> p = scratch_data.fe_values.quadrature_point(q);
897 *
for (
unsigned int k=0; k<dofs_per_cell; ++k)
899 * scratch_data.phi_u[k] = scratch_data.fe_values[velocity].value(k, q);
900 * scratch_data.div_phi_u[k] = scratch_data.fe_values[velocity].divergence (k, q);
901 * scratch_data.phi_p[k] = scratch_data.fe_values[pressure].value (k, q);
904 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
906 *
for (
unsigned int j=n_vel; j<dofs_per_cell; ++j)
908 *
double div_term = (- scratch_data.div_phi_u[i] * scratch_data.phi_p[j]
909 * - scratch_data.phi_p[i] * scratch_data.div_phi_u[j]) * scratch_data.fe_values.JxW(q);
912 * div_map[i][j] += div_term;
915 *
double source_term = -scratch_data.phi_p[i] * scratch_data.rhs_values[q] * scratch_data.fe_values.JxW(q);
917 *
if (
std::abs(scratch_data.phi_p[i]) > 1.e-12 ||
std::abs(source_term) > 1.e-12)
918 * copy_data.cell_vec[p][copy_data.local_dof_indices[i]] += source_term;
924 * Then, by making another pass, we compute the mass
matrix terms and incorporate the
925 * divergence form and RHS accordingly. This
second pass, allows us to know where
926 * the total contribution will be put in the nodal data structures, as with
this
927 * choice of quadrature rule and finite element only the basis
functions corresponding
928 * to the same quadrature points yield non-zero contribution.
931 *
for (
unsigned int q=0; q<n_q_points; ++q)
933 * std::set<types::global_dof_index> vel_indices;
934 *
const Point<dim> p = scratch_data.fe_values.quadrature_point(q);
936 *
for (
unsigned int k=0; k<dofs_per_cell; ++k)
938 * scratch_data.phi_u[k] = scratch_data.fe_values[velocity].value(k, q);
939 * scratch_data.div_phi_u[k] = scratch_data.fe_values[velocity].divergence (k, q);
940 * scratch_data.phi_p[k] = scratch_data.fe_values[pressure].value (k, q);
943 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
944 *
for (
unsigned int j=i; j<dofs_per_cell; ++j)
946 *
double mass_term = scratch_data.phi_u[i]
947 * * scratch_data.k_inverse_values[q]
948 * * scratch_data.phi_u[j]
949 * * scratch_data.fe_values.JxW(q);
953 * copy_data.cell_mat[
p][std::make_pair(copy_data.local_dof_indices[i], copy_data.local_dof_indices[j])] +=
955 * vel_indices.insert(i);
956 * copy_data.local_vel_indices[
p].insert(copy_data.local_dof_indices[j]);
960 *
for (
auto i : vel_indices)
961 *
for (auto el : div_map[i])
964 * copy_data.cell_mat[
p][
std::make_pair(copy_data.local_dof_indices[i],
965 * copy_data.local_dof_indices[el.
first])] += el.
second;
966 * copy_data.local_pres_indices[
p].insert(copy_data.local_dof_indices[el.first]);
972 * The pressure boundary conditions are computed as in @ref step_20
"step-20",
975 * std::map<types::global_dof_index,double> pres_bc;
976 *
for (
unsigned int face_no=0;
977 * face_no<GeometryInfo<dim>::faces_per_cell;
979 *
if (cell->at_boundary(face_no))
981 * scratch_data.fe_face_values.reinit (cell, face_no);
982 * pressure_bc.value_list(scratch_data.fe_face_values.get_quadrature_points(), scratch_data.pres_bc_values);
984 *
for (
unsigned int q=0; q<n_face_q_points; ++q)
985 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
987 *
double tmp = -(scratch_data.fe_face_values[velocity].value(i, q) *
988 * scratch_data.fe_face_values.normal_vector(q) *
989 * scratch_data.pres_bc_values[q] *
990 * scratch_data.fe_face_values.JxW(q));
993 * pres_bc[copy_data.local_dof_indices[i]] += tmp;
999 * ...but we distribute them to the corresponding nodal data structures
1002 *
for (
auto m : copy_data.cell_vec)
1003 *
for (unsigned
int i=0; i<dofs_per_cell; ++i)
1004 *
if (
std::abs(pres_bc[copy_data.local_dof_indices[i]]) > 1.e-12)
1005 * copy_data.cell_vec[m.first][copy_data.local_dof_indices[i]] += pres_bc[copy_data.local_dof_indices[i]];
1011 * Finally, <code>node_assembly()</code> takes care of all the
1012 * local computations via
WorkStream mechanism. Notice that the choice
1013 * of the quadrature rule here is dictated by the formulation of the
1014 * method. It has to be <code>degree+1</code> points Gauss-Lobatto
1015 *
for the
volume integrals and <code>degree</code>
for the face ones,
1016 * as mentioned in the introduction.
1019 *
template <
int dim>
1020 *
void MultipointMixedDarcyProblem<dim>::node_assembly()
1024 * dof_handler.distribute_dofs(fe);
1026 *
const std::vector<types::global_dof_index> dofs_per_component
1030 *
QGauss<dim-1> face_quad(degree);
1032 * n_v = dofs_per_component[0];
1033 * n_p = dofs_per_component[dim];
1035 * pres_rhs.reinit(n_p);
1038 * dof_handler.end(),
1040 * &MultipointMixedDarcyProblem::assemble_system_cell,
1041 * &MultipointMixedDarcyProblem::copy_cell_to_node,
1042 * DataStructures::NodeAssemblyScratchData<dim>(fe,
triangulation,quad,face_quad),
1043 * DataStructures::NodeAssemblyCopyData<dim>());
1049 * <a name=
"mfmfe.cc-Makingthesparsitypattern"></a>
1050 * <h4>Making the sparsity pattern</h4>
1054 * Having computed all the local contributions, we actually have
1055 * all the information needed to make a cell-centered sparsity
1057 * leads to a slower solution.
1060 *
template <
int dim>
1061 *
void MultipointMixedDarcyProblem<dim>::make_cell_centered_sp()
1066 * std::set<types::global_dof_index>::iterator pi_it, pj_it;
1067 *
unsigned int i, j;
1068 *
for (
auto el : node_matrix)
1070 * pi_it != pressure_indices[el.first].end();
1072 *
for (pj_it = pi_it, j = 0;
1073 * pj_it != pressure_indices[el.first].end();
1075 * dsp.add(*pi_it - n_v, *pj_it - n_v);
1079 * cell_centered_sp.copy_from(dsp);
1080 * pres_system_matrix.reinit (cell_centered_sp);
1087 * <a name=
"mfmfe.cc-Thelocaleliminationprocedure"></a>
1088 * <h4>The local elimination procedure</h4>
1092 * This
function finally performs the local elimination procedure.
1093 * Mathematically,
it follows the same idea as in computing the
1094 * Schur complement (as mentioned in the introduction) but we
do
1095 * so locally. Namely, local velocity DOFs are expressed in terms
1096 * of corresponding pressure
values, and then
used for the local
1100 *
template <
int dim>
1101 *
void MultipointMixedDarcyProblem<dim>::
1102 * nodal_elimination(
const typename DataStructures::PointToMatrixMap<dim>::iterator &n_it,
1103 * DataStructures::VertexEliminationScratchData &scratch_data,
1104 * DataStructures::NodeEliminationCopyData<dim> ©_data)
1106 *
unsigned int n_edges = velocity_indices.at((*n_it).first).size();
1107 *
unsigned int n_cells = pressure_indices.at((*n_it).first).size();
1109 * scratch_data.velocity_matrix.reinit(n_edges,n_edges);
1110 * copy_data.pressure_matrix.reinit(n_edges,n_cells);
1112 * copy_data.velocity_rhs.reinit(n_edges);
1113 * scratch_data.pressure_rhs.reinit(n_cells);
1116 * std::set<types::global_dof_index>::iterator vi_it, vj_it, p_it;
1118 *
for (vi_it = velocity_indices.at((*n_it).first).begin(), i = 0;
1119 * vi_it != velocity_indices.at((*n_it).first).end();
1123 *
for (vj_it = velocity_indices.at((*n_it).first).begin(), j = 0;
1124 * vj_it != velocity_indices.at((*n_it).first).end();
1127 * scratch_data.velocity_matrix.add(i, j, node_matrix[(*n_it).first][std::make_pair(*vi_it, *vj_it)]);
1129 * scratch_data.velocity_matrix.add(j, i, node_matrix[(*n_it).first][std::make_pair(*vi_it, *vj_it)]);
1132 *
for (p_it = pressure_indices.at((*n_it).first).begin(), j = 0;
1133 * p_it != pressure_indices.at((*n_it).first).end();
1135 * copy_data.pressure_matrix.add(i, j, node_matrix[(*n_it).first][std::make_pair(*vi_it, *p_it)]);
1137 * copy_data.velocity_rhs(i) += node_rhs.at((*n_it).first)[*vi_it];
1140 *
for (p_it = pressure_indices.at((*n_it).first).begin(), i = 0;
1141 * p_it != pressure_indices.at((*n_it).first).end();
1143 * scratch_data.pressure_rhs(i) += node_rhs.at((*n_it).first)[*p_it];
1146 * copy_data.Ainverse.reinit(n_edges,n_edges);
1148 * scratch_data.tmp_rhs1.reinit(n_edges);
1149 * scratch_data.tmp_rhs2.reinit(n_edges);
1150 * scratch_data.tmp_rhs3.reinit(n_cells);
1152 * copy_data.Ainverse.invert(scratch_data.velocity_matrix);
1153 * copy_data.node_pres_matrix.reinit(n_cells, n_cells);
1154 * copy_data.node_pres_rhs = scratch_data.pressure_rhs;
1156 * copy_data.node_pres_matrix = 0;
1157 * copy_data.node_pres_matrix.triple_product(copy_data.Ainverse,
1158 * copy_data.pressure_matrix,
1159 * copy_data.pressure_matrix,
true,
false);
1161 * copy_data.Ainverse.vmult(scratch_data.tmp_rhs1, copy_data.velocity_rhs,
false);
1162 * copy_data.pressure_matrix.Tvmult(scratch_data.tmp_rhs3, scratch_data.tmp_rhs1,
false);
1163 * copy_data.node_pres_rhs *= -1.0;
1164 * copy_data.node_pres_rhs += scratch_data.tmp_rhs3;
1166 * copy_data.p = (*n_it).first;
1172 * Each node
's pressure system is then distributed to a global pressure
1173 * system, using the indices we computed in the previous stages.
1176 * template <int dim>
1177 * void MultipointMixedDarcyProblem<dim>::
1178 * copy_node_to_system(const DataStructures::NodeEliminationCopyData<dim> ©_data)
1180 * A_inverse[copy_data.p] = copy_data.Ainverse;
1181 * pressure_matrix[copy_data.p] = copy_data.pressure_matrix;
1182 * velocity_rhs[copy_data.p] = copy_data.velocity_rhs;
1185 * std::set<types::global_dof_index>::iterator pi_it, pj_it;
1187 * for (pi_it = pressure_indices[copy_data.p].begin(), i = 0;
1188 * pi_it != pressure_indices[copy_data.p].end();
1192 * for (pj_it = pressure_indices[copy_data.p].begin(), j = 0;
1193 * pj_it != pressure_indices[copy_data.p].end();
1195 * pres_system_matrix.add(*pi_it - n_v, *pj_it - n_v, copy_data.node_pres_matrix(i, j));
1197 * pres_rhs(*pi_it - n_v) += copy_data.node_pres_rhs(i);
1205 * The @ref WorkStream mechanism is again used for the assembly
1206 * of the global system for the pressure variable, where the
1207 * previous functions are used to perform local computations.
1210 * template <int dim>
1211 * void MultipointMixedDarcyProblem<dim>::pressure_assembly()
1213 * TimerOutput::Scope t(computing_timer, "Pressure matrix assembly");
1215 * QGaussLobatto<dim> quad(degree+1);
1216 * QGauss<dim-1> face_quad(degree);
1218 * pres_rhs.reinit(n_p);
1220 * WorkStream::run(node_matrix.begin(),
1221 * node_matrix.end(),
1223 * &MultipointMixedDarcyProblem::nodal_elimination,
1224 * &MultipointMixedDarcyProblem::copy_node_to_system,
1225 * DataStructures::VertexEliminationScratchData(),
1226 * DataStructures::NodeEliminationCopyData<dim>());
1234 * <a name="mfmfe.cc-Velocitysolutionrecovery"></a>
1235 * <h4>Velocity solution recovery</h4>
1239 * After solving for the pressure variable, we want to follow
1240 * the above procedure backwards, in order to obtain the
1241 * velocity solution (again, this is similar in nature to the
1242 * Schur complement approach, see @ref step_20 "step-20", but here it is done
1243 * locally at each node). We have almost everything computed and
1244 * stored already, including inverses of local mass matrices,
1245 * so the following is a relatively straightforward implementation.
1248 * template <int dim>
1249 * void MultipointMixedDarcyProblem<dim>::
1250 * velocity_assembly (const typename DataStructures::PointToMatrixMap<dim>::iterator &n_it,
1251 * DataStructures::VertexEliminationScratchData &scratch_data,
1252 * DataStructures::NodeEliminationCopyData<dim> ©_data)
1254 * unsigned int n_edges = velocity_indices.at((*n_it).first).size();
1255 * unsigned int n_cells = pressure_indices.at((*n_it).first).size();
1257 * scratch_data.tmp_rhs1.reinit(n_edges);
1258 * scratch_data.tmp_rhs2.reinit(n_edges);
1259 * scratch_data.tmp_rhs3.reinit(n_cells);
1260 * scratch_data.local_pressure_solution.reinit(n_cells);
1262 * copy_data.vertex_vel_solution.reinit(n_edges);
1264 * std::set<types::global_dof_index>::iterator p_it;
1267 * for (p_it = pressure_indices[(*n_it).first].begin(), i = 0;
1268 * p_it != pressure_indices[(*n_it).first].end();
1270 * scratch_data.local_pressure_solution(i) = pres_solution(*p_it - n_v);
1272 * pressure_matrix[(*n_it).first].vmult(scratch_data.tmp_rhs2, scratch_data.local_pressure_solution, false);
1273 * scratch_data.tmp_rhs2 *= -1.0;
1274 * scratch_data.tmp_rhs2+=velocity_rhs[(*n_it).first];
1275 * A_inverse[(*n_it).first].vmult(copy_data.vertex_vel_solution, scratch_data.tmp_rhs2, false);
1277 * copy_data.p = (*n_it).first;
1283 * Copy nodal velocities to a global solution vector by using
1284 * local computations and indices from early stages.
1287 * template <int dim>
1288 * void MultipointMixedDarcyProblem<dim>::
1289 * copy_node_velocity_to_global(const DataStructures::NodeEliminationCopyData<dim> ©_data)
1291 * std::set<types::global_dof_index>::iterator vi_it;
1294 * for (vi_it = velocity_indices[copy_data.p].begin(), i = 0;
1295 * vi_it != velocity_indices[copy_data.p].end();
1297 * vel_solution(*vi_it) += copy_data.vertex_vel_solution(i);
1303 * Use @ref WorkStream to run everything concurrently.
1306 * template <int dim>
1307 * void MultipointMixedDarcyProblem<dim>::velocity_recovery()
1309 * TimerOutput::Scope t(computing_timer, "Velocity solution recovery");
1311 * QGaussLobatto<dim> quad(degree+1);
1312 * QGauss<dim-1> face_quad(degree);
1314 * vel_solution.reinit(n_v);
1316 * WorkStream::run(node_matrix.begin(),
1317 * node_matrix.end(),
1319 * &MultipointMixedDarcyProblem::velocity_assembly,
1320 * &MultipointMixedDarcyProblem::copy_node_velocity_to_global,
1321 * DataStructures::VertexEliminationScratchData(),
1322 * DataStructures::NodeEliminationCopyData<dim>());
1324 * solution.reinit(2);
1325 * solution.block(0) = vel_solution;
1326 * solution.block(1) = pres_solution;
1327 * solution.collect_sizes();
1335 * <a name="mfmfe.cc-Pressuresystemsolver"></a>
1336 * <h4>Pressure system solver</h4>
1340 * The solver part is trivial. We use the CG solver with no
1341 * preconditioner for simplicity.
1344 * template <int dim>
1345 * void MultipointMixedDarcyProblem<dim>::solve_pressure()
1347 * TimerOutput::Scope t(computing_timer, "Pressure CG solve");
1349 * pres_solution.reinit(n_p);
1351 * SolverControl solver_control (static_cast<int>(2.0*n_p), 1e-10);
1352 * SolverCG<> solver (solver_control);
1354 * PreconditionIdentity identity;
1355 * solver.solve(pres_system_matrix, pres_solution, pres_rhs, identity);
1363 * <a name="mfmfe.cc-Postprocessing"></a>
1364 * <h3>Postprocessing</h3>
1368 * We have two postprocessing steps here, first one computes the
1369 * errors in order to populate the convergence tables. The other
1370 * one takes care of the output of the solutions in <code>.vtk</code>
1376 * <a name="mfmfe.cc-Computeerrors"></a>
1377 * <h4>Compute errors</h4>
1381 * The implementation of this function is almost identical to @ref step_20 "step-20".
1382 * We use @ref ComponentSelectFunction as masks to use the right
1383 * solution component (velocity or pressure) and @ref integrate_difference
1384 * to compute the errors. Since we also want to compute Hdiv seminorm of the
1385 * velocity error, one must provide gradients in the <code>ExactSolution</code>
1386 * class implementation to avoid exceptions. The only noteworthy thing here
1387 * is that we again use lower order quadrature rule instead of projecting the
1388 * solution to an appropriate space in order to show superconvergence, which is
1389 * mathematically justified.
1392 * template <int dim>
1393 * void MultipointMixedDarcyProblem<dim>::compute_errors(const unsigned cycle)
1395 * TimerOutput::Scope t(computing_timer, "Compute errors");
1397 * const ComponentSelectFunction<dim> pressure_mask(dim, dim+1);
1398 * const ComponentSelectFunction<dim> velocity_mask(std::make_pair(0, dim), dim+1);
1400 * ExactSolution<dim> exact_solution;
1402 * Vector<double> cellwise_errors (triangulation.n_active_cells());
1404 * QTrapezoid<1> q_trapez;
1405 * QIterated<dim> quadrature(q_trapez,degree+2);
1406 * QGauss<dim> quadrature_super(degree);
1408 * VectorTools::integrate_difference (dof_handler, solution, exact_solution,
1409 * cellwise_errors, quadrature,
1410 * VectorTools::L2_norm,
1412 * const double p_l2_error = cellwise_errors.l2_norm();
1414 * VectorTools::integrate_difference (dof_handler, solution, exact_solution,
1415 * cellwise_errors, quadrature_super,
1416 * VectorTools::L2_norm,
1418 * const double p_l2_mid_error = cellwise_errors.l2_norm();
1420 * VectorTools::integrate_difference (dof_handler, solution, exact_solution,
1421 * cellwise_errors, quadrature,
1422 * VectorTools::L2_norm,
1424 * const double u_l2_error = cellwise_errors.l2_norm();
1426 * VectorTools::integrate_difference (dof_handler, solution, exact_solution,
1427 * cellwise_errors, quadrature,
1428 * VectorTools::Hdiv_seminorm,
1430 * const double u_hd_error = cellwise_errors.l2_norm();
1432 * const unsigned int n_active_cells=triangulation.n_active_cells();
1433 * const unsigned int n_dofs=dof_handler.n_dofs();
1435 * convergence_table.add_value("cycle", cycle);
1436 * convergence_table.add_value("cells", n_active_cells);
1437 * convergence_table.add_value("dofs", n_dofs);
1438 * convergence_table.add_value("Velocity,L2", u_l2_error);
1439 * convergence_table.add_value("Velocity,Hdiv", u_hd_error);
1440 * convergence_table.add_value("Pressure,L2", p_l2_error);
1441 * convergence_table.add_value("Pressure,L2-nodal", p_l2_mid_error);
1449 * <a name="mfmfe.cc-Outputresults"></a>
1450 * <h4>Output results</h4>
1454 * This function also follows the same idea as in @ref step_20 "step-20" tutorial
1455 * program. The only modification to it is the part involving
1456 * a convergence table.
1459 * template <int dim>
1460 * void MultipointMixedDarcyProblem<dim>::output_results(const unsigned int cycle, const unsigned int refine)
1462 * TimerOutput::Scope t(computing_timer, "Output results");
1464 * std::vector<std::string> solution_names(dim, "u");
1465 * solution_names.push_back ("p");
1466 * std::vector<DataComponentInterpretation::DataComponentInterpretation>
1467 * interpretation (dim, DataComponentInterpretation::component_is_part_of_vector);
1468 * interpretation.push_back (DataComponentInterpretation::component_is_scalar);
1470 * DataOut<dim> data_out;
1471 * data_out.add_data_vector (dof_handler, solution, solution_names, interpretation);
1472 * data_out.build_patches ();
1474 * std::ofstream output ("solution" + std::to_string(dim) + "d-" + std::to_string(cycle) + ".vtk");
1475 * data_out.write_vtk (output);
1477 * convergence_table.set_precision("Velocity,L2", 3);
1478 * convergence_table.set_precision("Velocity,Hdiv", 3);
1479 * convergence_table.set_precision("Pressure,L2", 3);
1480 * convergence_table.set_precision("Pressure,L2-nodal", 3);
1481 * convergence_table.set_scientific("Velocity,L2", true);
1482 * convergence_table.set_scientific("Velocity,Hdiv", true);
1483 * convergence_table.set_scientific("Pressure,L2", true);
1484 * convergence_table.set_scientific("Pressure,L2-nodal", true);
1485 * convergence_table.set_tex_caption("cells", "\\# cells");
1486 * convergence_table.set_tex_caption("dofs", "\\# dofs");
1487 * convergence_table.set_tex_caption("Velocity,L2", " \\|\\u - \\u_h\\|_{L^2} ");
1488 * convergence_table.set_tex_caption("Velocity,Hdiv", " \\|\\nabla\\cdot(\\u - \\u_h)\\|_{L^2} ");
1489 * convergence_table.set_tex_caption("Pressure,L2", " \\|p - p_h\\|_{L^2} ");
1490 * convergence_table.set_tex_caption("Pressure,L2-nodal", " \\|Qp - p_h\\|_{L^2} ");
1491 * convergence_table.set_tex_format("cells", "r");
1492 * convergence_table.set_tex_format("dofs", "r");
1494 * convergence_table.evaluate_convergence_rates("Velocity,L2", ConvergenceTable::reduction_rate_log2);
1495 * convergence_table.evaluate_convergence_rates("Velocity,Hdiv", ConvergenceTable::reduction_rate_log2);
1496 * convergence_table.evaluate_convergence_rates("Pressure,L2", ConvergenceTable::reduction_rate_log2);
1497 * convergence_table.evaluate_convergence_rates("Pressure,L2-nodal", ConvergenceTable::reduction_rate_log2);
1499 * std::ofstream error_table_file("error" + std::to_string(dim) + "d.tex");
1501 * if (cycle == refine-1)
1503 * convergence_table.write_text(std::cout);
1504 * convergence_table.write_tex(error_table_file);
1513 * <a name="mfmfe.cc-Runfunction"></a>
1514 * <h3>Run function</h3>
1518 * The driver method <code>run()</code>
1519 * takes care of mesh generation and arranging calls to member methods in
1520 * the right way. It also resets data structures and clear triangulation and
1521 * DOF handler as we run the method on a sequence of refinements in order
1522 * to record convergence rates.
1525 * template <int dim>
1526 * void MultipointMixedDarcyProblem<dim>::run(const unsigned int refine)
1528 * Assert(refine > 0, ExcMessage("Must at least have 1 refinement cycle!"));
1530 * dof_handler.clear();
1531 * triangulation.clear();
1532 * convergence_table.clear();
1534 * for (unsigned int cycle=0; cycle<refine; ++cycle)
1540 * We first generate the hyper cube and refine it twice
1541 * so that we could distort the grid slightly and
1542 * demonstrate the method's ability to work in such a
1554 * make_cell_centered_sp();
1555 * pressure_assembly();
1556 * solve_pressure ();
1557 * velocity_recovery ();
1558 * compute_errors (cycle);
1559 * output_results (cycle, refine);
1560 * reset_data_structures ();
1562 * computing_timer.print_summary ();
1563 * computing_timer.reset ();
1572 * <a name=
"mfmfe.cc-Thecodemaincodefunction"></a>
1573 * <h3>The <code>main</code>
function</h3>
1577 * In the main functione we pass the order of the Finite Element as an argument
1578 * to the constructor of the Multipoint Flux Mixed Darcy problem, and the number
1579 * of refinement cycles as an argument
for the
run method.
1586 *
using namespace dealii;
1587 *
using namespace MFMFE;
1591 * MultipointMixedDarcyProblem<2> mfmfe_problem(2);
1592 * mfmfe_problem.run(6);
1594 *
catch (std::exception &exc)
1596 * std::cerr << std::endl << std::endl
1597 * <<
"----------------------------------------------------"
1599 * std::cerr <<
"Exception on processing: " << std::endl
1600 * << exc.what() << std::endl
1601 * <<
"Aborting!" << std::endl
1602 * <<
"----------------------------------------------------"
1609 * std::cerr << std::endl << std::endl
1610 * <<
"----------------------------------------------------"
1612 * std::cerr <<
"Unknown exception!" << std::endl
1613 * <<
"Aborting!" << std::endl
1614 * <<
"----------------------------------------------------"
virtual void vector_gradient(const Point< dim > &p, std::vector< Tensor< 1, dim, RangeNumberType > > &gradients) const
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
virtual void vector_value(const Point< dim > &p, Vector< RangeNumberType > &values) const
static void set_thread_limit(const unsigned int max_threads=numbers::invalid_unsigned_int)
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< value_type > &values) const
face_iterator end_face() const
unsigned int n_active_cells() const
unsigned int n_vertices() const
active_face_iterator begin_active_face() const
#define Assert(cond, exc)
typename ActiveSelector::active_cell_iterator active_cell_iterator
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
void component_wise(DoFHandler< dim, spacedim > &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
@ matrix
Contents is actually a matrix.
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
constexpr T pow(const T base, const int iexp)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
void run(const std::vector< std::vector< Iterator > > &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
void copy(const T *begin, const T *end, U *dest)
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
const InputIterator OutputIterator const Function & function
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation