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=
"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;
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);
224 *
class PressureBoundaryValues :
public Function<dim>
227 * PressureBoundaryValues () :
Function<dim>(1) {}
230 *
const unsigned int component = 0)
const;
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;
252 *
class ExactSolution :
public Function<dim>
255 * ExactSolution () :
Function<dim>(dim+1) {}
257 *
virtual void vector_value (
const Point<dim> &p,
260 *
virtual void vector_gradient (
const Point<dim> &p,
266 * ExactSolution<dim>::vector_value (
const Point<dim> &p,
269 *
Assert (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);
280 * values(2) = (x*x*x)*(y*y*y*y)+
cos(x*y)*
sin(x*y)+x*x;
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"));
325 *
virtual void value_list (
const std::vector<
Point<dim> > &points,
331 * KInverse<dim>::value_list (
const std::vector<
Point<dim> > &points,
334 *
Assert (points.size() == values.size(),
337 *
for (
unsigned int p=0; p<points.size(); ++p)
339 * values[p].clear ();
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>
383 * <a name=
"Includefiles"></a>
384 * <h3>Include files</h3>
388 * As usual, the list of necessary header files. There is not
389 * much
new here, the files are included in order
390 * base-lac-grid-dofs-numerics followed by the
C++ headers.
393 * #include <deal.II/base/convergence_table.h>
394 * #include <deal.II/base/quadrature_lib.h>
395 * #include <deal.II/base/logstream.h>
396 * #include <deal.II/base/timer.h>
397 * #include <deal.II/base/work_stream.h>
399 * #include <deal.II/lac/full_matrix.h>
400 * #include <deal.II/lac/solver_cg.h>
401 * #include <deal.II/lac/block_sparse_matrix.h>
402 * #include <deal.II/lac/block_vector.h>
403 * #include <deal.II/lac/precondition.h>
405 * #include <deal.II/grid/grid_generator.h>
406 * #include <deal.II/grid/grid_tools.h>
407 * #include <deal.II/grid/grid_in.h>
408 * #include <deal.II/grid/tria.h>
409 * #include <deal.II/dofs/dof_renumbering.h>
410 * #include <deal.II/dofs/dof_tools.h>
411 * #include <deal.II/fe/fe_dgq.h>
412 * #include <deal.II/fe/fe_system.h>
413 * #include <deal.II/fe/fe_tools.h>
414 * #include <deal.II/numerics/vector_tools.h>
415 * #include <deal.II/numerics/matrix_tools.h>
416 * #include <deal.II/numerics/data_out.h>
419 * #include <unordered_map>
423 * This is a header needed
for the purposes of the
424 * multipoint flux mixed method, as it declares the
425 *
new enhanced Raviart-Thomas finite element.
428 * #include <deal.II/fe/fe_rt_bubbles.h>
432 * For the sake of readability, the classes representing
433 * data, i.e. RHS, BCs, permeability tensor and the exact
434 * solution are placed in a file data.h which is included
442 * As
always the program is in the
namespace of its own with
443 * the deal.II classes and
functions imported into it
453 * <a name=
"Definitionofmultipointfluxassemblydatastructures"></a>
454 * <h3>Definition of multipoint flux assembly data structures</h3>
458 * The main idea of the MFMFE method is to perform local elimination
459 * of the velocity variables in order to obtain the resulting
460 * pressure system. Since in deal.II assembly happens cell-wise,
461 * some extra work needs to be done in order to get the local
462 * mass matrices @f$A_i@f$ and the corresponding to them @f$B_i@f$.
465 *
namespace DataStructures
469 * This will be achieved by assembling cell-wise, but instead of placing
470 * the terms into a global system
matrix, they will populate node-associated
471 * full matrices. For
this, a data structure with fast lookup is crucial, hence
478 *
size_t operator()(
const Point<dim> &p)
const
481 * h1 = std::hash<double>()(p[0]);
488 * h2 = std::hash<double>()(p[1]);
491 * h2 = std::hash<double>()(p[1]);
492 * h3 = std::hash<double>()(p[2]);
493 *
return (h1 ^ (h2 << 1)) ^ h3;
502 * Here, the actual hash-tables are defined. We use the
C++ STL <code>unordered_map</code>,
503 * with the hash
function specified above. For convenience these are aliased as follows
507 *
using PointToMatrixMap = std::unordered_map<Point<dim>, std::map<std::pair<types::global_dof_index,types::global_dof_index>,
double>, hash_points<dim>>;
510 *
using PointToVectorMap = std::unordered_map<Point<dim>, std::map<types::global_dof_index, double>, hash_points<dim>>;
513 *
using PointToIndexMap = std::unordered_map<Point<dim>, std::set<types::global_dof_index>, hash_points<dim>>;
517 * Next, since
this particular program allows
for the use of
518 * multiple threads, the helper CopyData structures
519 * are defined. There are two kinds of these,
one is used
520 *
for the copying cell-wise contributions to the corresponging
521 * node-associated data structures...
525 *
struct NodeAssemblyCopyData
527 * PointToMatrixMap<dim> cell_mat;
528 * PointToVectorMap<dim> cell_vec;
529 * PointToIndexMap<dim> local_pres_indices;
530 * PointToIndexMap<dim> local_vel_indices;
531 * std::vector<types::global_dof_index> local_dof_indices;
536 * ... and the other
one for the actual process of
537 * local velocity elimination and assembling the global
542 *
struct NodeEliminationCopyData
555 * Similarly, two ScratchData classes are defined.
556 * One
for the assembly part, where we need
558 *
for the basis fuctions...
562 *
struct NodeAssemblyScratchData
569 * NodeAssemblyScratchData (
const NodeAssemblyScratchData &scratch_data);
573 * std::vector<unsigned int> n_faces_at_vertex;
575 *
const unsigned long num_cells;
577 * std::vector<Tensor<2,dim>> k_inverse_values;
578 * std::vector<double> rhs_values;
579 * std::vector<double> pres_bc_values;
581 * std::vector<Tensor<1,dim> > phi_u;
582 * std::vector<double> div_phi_u;
583 * std::vector<double> phi_p;
587 * NodeAssemblyScratchData<dim>::
597 * fe_face_values (fe,
601 * num_cells(tria.n_active_cells()),
602 * k_inverse_values(quad.size()),
603 * rhs_values(quad.size()),
604 * pres_bc_values(f_quad.size()),
605 * phi_u(fe.dofs_per_cell),
606 * div_phi_u(fe.dofs_per_cell),
607 * phi_p(fe.dofs_per_cell)
609 * n_faces_at_vertex.resize(tria.n_vertices(), 0);
612 *
for (; face != endf; ++face)
613 *
for (
unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
614 * n_faces_at_vertex[face->vertex_index(v)] += 1;
618 * NodeAssemblyScratchData<dim>::
619 * NodeAssemblyScratchData (
const NodeAssemblyScratchData &scratch_data)
621 * fe_values (scratch_data.fe_values.get_fe(),
622 * scratch_data.fe_values.get_quadrature(),
625 * fe_face_values (scratch_data.fe_face_values.get_fe(),
626 * scratch_data.fe_face_values.get_quadrature(),
629 * n_faces_at_vertex(scratch_data.n_faces_at_vertex),
630 * num_cells(scratch_data.num_cells),
631 * k_inverse_values(scratch_data.k_inverse_values),
632 * rhs_values(scratch_data.rhs_values),
633 * pres_bc_values(scratch_data.pres_bc_values),
634 * phi_u(scratch_data.phi_u),
635 * div_phi_u(scratch_data.div_phi_u),
636 * phi_p(scratch_data.phi_p)
641 * ...and the other, simpler
one,
for the velocity elimination and recovery
644 *
struct VertexEliminationScratchData
646 * VertexEliminationScratchData () =
default;
647 * VertexEliminationScratchData (
const VertexEliminationScratchData &scratch_data);
658 * VertexEliminationScratchData::
659 * VertexEliminationScratchData (
const VertexEliminationScratchData &scratch_data)
661 * velocity_matrix(scratch_data.velocity_matrix),
662 * pressure_rhs(scratch_data.pressure_rhs),
663 * local_pressure_solution(scratch_data.local_pressure_solution),
664 * tmp_rhs1(scratch_data.tmp_rhs1),
665 * tmp_rhs2(scratch_data.tmp_rhs2),
666 * tmp_rhs3(scratch_data.tmp_rhs3)
675 * <a name=
"ThecodeMultipointMixedDarcyProblemcodeclasstemplate"></a>
676 * <h3>The <code>MultipointMixedDarcyProblem</code>
class template</h3>
680 * The main
class, besides the constructor and destructor, has only
one public member
681 * <code>
run()</code>, similarly to the tutorial programs. The
private members can
682 * be grouped into the ones that are used
for the cell-wise assembly, vertex elimination,
683 * pressure solve, vertex velocity recovery and postprocessing. Apart from the
684 * MFMFE-specific data structures, the rest of the members should look familiar.
688 *
class MultipointMixedDarcyProblem
691 * MultipointMixedDarcyProblem (
const unsigned int degree);
692 * ~MultipointMixedDarcyProblem ();
696 * DataStructures::NodeAssemblyScratchData<dim> &scratch_data,
697 * DataStructures::NodeAssemblyCopyData<dim> ©_data);
698 *
void copy_cell_to_node(
const DataStructures::NodeAssemblyCopyData<dim> ©_data);
699 *
void node_assembly();
700 *
void make_cell_centered_sp ();
701 *
void nodal_elimination(
const typename DataStructures::PointToMatrixMap<dim>::iterator &n_it,
702 * DataStructures::VertexEliminationScratchData &scratch_data,
703 * DataStructures::NodeEliminationCopyData<dim> ©_data);
704 *
void copy_node_to_system(
const DataStructures::NodeEliminationCopyData<dim> ©_data);
705 *
void pressure_assembly ();
706 *
void solve_pressure ();
707 *
void velocity_assembly (
const typename DataStructures::PointToMatrixMap<dim>::iterator &n_it,
708 * DataStructures::VertexEliminationScratchData &scratch_data,
709 * DataStructures::NodeEliminationCopyData<dim> ©_data);
710 *
void copy_node_velocity_to_global(
const DataStructures::NodeEliminationCopyData<dim> ©_data);
711 *
void velocity_recovery ();
712 *
void reset_data_structures ();
713 *
void compute_errors (
const unsigned int cycle);
714 *
void output_results (
const unsigned int cycle,
const unsigned int refine);
716 *
const unsigned int degree;
726 * std::unordered_map<Point<dim>,
FullMatrix<double>, DataStructures::hash_points<dim>> pressure_matrix;
727 * std::unordered_map<Point<dim>,
FullMatrix<double>, DataStructures::hash_points<dim>> A_inverse;
728 * std::unordered_map<Point<dim>,
Vector<double>, DataStructures::hash_points<dim>> velocity_rhs;
730 * DataStructures::PointToMatrixMap<dim> node_matrix;
731 * DataStructures::PointToVectorMap<dim> node_rhs;
733 * DataStructures::PointToIndexMap<dim> pressure_indices;
734 * DataStructures::PointToIndexMap<dim> velocity_indices;
736 *
unsigned long n_v, n_p;
748 * <a name=
"Constructoranddestructorcodereset_data_structurescode"></a>
749 * <h4>Constructor and destructor, <code>reset_data_structures</code></h4>
753 * In the constructor of
this class, we store the
value that was
754 * passed in concerning the degree of the finite elements we shall use (a
756 * and then construct the vector valued element belonging to the space @f$V_h^k@f$ described
757 * in the introduction. The constructor also takes care of initializing the
758 * computing timer, as it is of interest
for us how well our method performs.
762 * MultipointMixedDarcyProblem<dim>::MultipointMixedDarcyProblem (
const unsigned int degree)
775 * The destructor clears the <code>dof_handler</code> and
776 * all of the data structures we used
for the method.
780 * MultipointMixedDarcyProblem<dim>::~MultipointMixedDarcyProblem()
782 * reset_data_structures ();
783 * dof_handler.clear();
789 * This method clears all the data that was used after
one refinement
794 *
void MultipointMixedDarcyProblem<dim>::reset_data_structures ()
796 * pressure_indices.clear();
797 * velocity_indices.clear();
798 * velocity_rhs.clear();
800 * pressure_matrix.clear();
801 * node_matrix.clear();
809 * <a name=
"Cellwiseassemblyandcreationofthelocalnodalbaseddatastructures"></a>
810 * <h4>Cell-wise assembly and creation of the local, nodal-based data structures</h4>
814 * First, the
function that copies local cell contributions to the corresponding nodal
815 * matrices and vectors is defined. It places the values obtained from local cell integration
816 * into the correct place in a
matrix/vector corresponging to a specific node.
820 *
void MultipointMixedDarcyProblem<dim>::copy_cell_to_node(
const DataStructures::NodeAssemblyCopyData<dim> ©_data)
822 *
for (
auto m : copy_data.cell_mat)
824 *
for (
auto p : m.second)
825 * node_matrix[m.first][p.first] += p.second;
827 *
for (
auto p : copy_data.cell_vec.at(m.first))
828 * node_rhs[m.first][p.first] += p.second;
830 *
for (
auto p : copy_data.local_pres_indices.at(m.first))
831 * pressure_indices[m.first].insert(p);
833 *
for (
auto p : copy_data.local_vel_indices.at(m.first))
834 * velocity_indices[m.first].insert(p);
842 * Second, the
function that does the cell assembly is defined. While it is
843 * similar to the tutorial programs in a way it uses scrath and
copy data
844 * structures, the need to localize the DOFs leads to several differences.
848 *
void MultipointMixedDarcyProblem<dim>::
850 * DataStructures::NodeAssemblyScratchData<dim> &scratch_data,
851 * DataStructures::NodeAssemblyCopyData<dim> ©_data)
853 * copy_data.cell_mat.clear();
854 * copy_data.cell_vec.clear();
855 * copy_data.local_vel_indices.clear();
856 * copy_data.local_pres_indices.clear();
858 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
859 *
const unsigned int n_q_points = scratch_data.fe_values.get_quadrature().size();
860 *
const unsigned int n_face_q_points = scratch_data.fe_face_values.get_quadrature().size();
862 * copy_data.local_dof_indices.resize(dofs_per_cell);
863 * cell->get_dof_indices (copy_data.local_dof_indices);
865 * scratch_data.fe_values.reinit (cell);
867 *
const KInverse<dim> k_inverse;
868 *
const RightHandSide<dim> rhs;
869 *
const PressureBoundaryValues<dim> pressure_bc;
871 * k_inverse.value_list (scratch_data.fe_values.get_quadrature_points(), scratch_data.k_inverse_values);
872 * rhs.value_list(scratch_data.fe_values.get_quadrature_points(), scratch_data.rhs_values);
877 *
const unsigned int n_vel = dim*
pow(degree+1,dim);
878 * std::unordered_map<unsigned int, std::unordered_map<unsigned int, double>> div_map;
882 * One, we need to be able to
assemble the communication between velocity and
883 * pressure variables and put it on the right place in our
final, local version
884 * of the B
matrix. This is a little messy, as such communication is not in fact
885 * local, so we
do it in two steps. First, we compute all relevant LHS and RHS
888 *
for (
unsigned int q=0; q<n_q_points; ++q)
890 *
const Point<dim> p = scratch_data.fe_values.quadrature_point(q);
892 *
for (
unsigned int k=0; k<dofs_per_cell; ++k)
894 * scratch_data.phi_u[k] = scratch_data.fe_values[velocity].value(k, q);
895 * scratch_data.div_phi_u[k] = scratch_data.fe_values[velocity].divergence (k, q);
896 * scratch_data.phi_p[k] = scratch_data.fe_values[pressure].value (k, q);
899 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
901 *
for (
unsigned int j=n_vel; j<dofs_per_cell; ++j)
903 *
double div_term = (- scratch_data.div_phi_u[i] * scratch_data.phi_p[j]
904 * - scratch_data.phi_p[i] * scratch_data.div_phi_u[j]) * scratch_data.fe_values.JxW(q);
906 *
if (std::abs(div_term) > 1.e-12)
907 * div_map[i][j] += div_term;
910 *
double source_term = -scratch_data.phi_p[i] * scratch_data.rhs_values[q] * scratch_data.fe_values.JxW(q);
912 *
if (std::abs(scratch_data.phi_p[i]) > 1.e-12 ||
std::abs(source_term) > 1.e-12)
913 * copy_data.cell_vec[p][copy_data.local_dof_indices[i]] += source_term;
919 * Then, by making another pass, we compute the mass
matrix terms and incorporate the
920 * divergence form and RHS accordingly. This
second pass, allows us to know where
921 * the total contribution will be put in the nodal data structures, as with
this
922 * choice of quadrature rule and finite element only the basis
functions corresponding
923 * to the same quadrature points yield non-
zero contribution.
926 *
for (
unsigned int q=0; q<n_q_points; ++q)
928 * std::set<types::global_dof_index> vel_indices;
929 *
const Point<dim> p = scratch_data.fe_values.quadrature_point(q);
931 *
for (
unsigned int k=0; k<dofs_per_cell; ++k)
933 * scratch_data.phi_u[k] = scratch_data.fe_values[velocity].value(k, q);
934 * scratch_data.div_phi_u[k] = scratch_data.fe_values[velocity].divergence (k, q);
935 * scratch_data.phi_p[k] = scratch_data.fe_values[pressure].value (k, q);
938 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
939 *
for (
unsigned int j=i; j<dofs_per_cell; ++j)
941 *
double mass_term = scratch_data.phi_u[i]
942 * * scratch_data.k_inverse_values[q]
943 * * scratch_data.phi_u[j]
944 * * scratch_data.fe_values.JxW(q);
946 *
if (std::abs(mass_term) > 1.e-12)
948 * copy_data.cell_mat[p][std::make_pair(copy_data.local_dof_indices[i], copy_data.local_dof_indices[j])] +=
950 * vel_indices.insert(i);
951 * copy_data.local_vel_indices[p].insert(copy_data.local_dof_indices[j]);
955 *
for (
auto i : vel_indices)
956 *
for (
auto el : div_map[i])
957 *
if (std::abs(el.second) > 1.e-12)
959 * copy_data.cell_mat[p][std::make_pair(copy_data.local_dof_indices[i],
960 * copy_data.local_dof_indices[el.first])] += el.second;
961 * copy_data.local_pres_indices[p].insert(copy_data.local_dof_indices[el.first]);
967 * The pressure boundary conditions are computed as in @ref step_20
"step-20",
970 * std::map<types::global_dof_index,double> pres_bc;
971 *
for (
unsigned int face_no=0;
972 * face_no<GeometryInfo<dim>::faces_per_cell;
974 *
if (cell->at_boundary(face_no))
976 * scratch_data.fe_face_values.reinit (cell, face_no);
977 * pressure_bc.value_list(scratch_data.fe_face_values.get_quadrature_points(), scratch_data.pres_bc_values);
979 *
for (
unsigned int q=0; q<n_face_q_points; ++q)
980 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
982 *
double tmp = -(scratch_data.fe_face_values[velocity].value(i, q) *
983 * scratch_data.fe_face_values.normal_vector(q) *
984 * scratch_data.pres_bc_values[q] *
985 * scratch_data.fe_face_values.JxW(q));
987 *
if (std::abs(tmp) > 1.e-12)
988 * pres_bc[copy_data.local_dof_indices[i]] += tmp;
994 * ...but we distribute them to the corresponding nodal data structures
997 *
for (
auto m : copy_data.cell_vec)
998 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
999 *
if (std::abs(pres_bc[copy_data.local_dof_indices[i]]) > 1.e-12)
1000 * copy_data.cell_vec[m.first][copy_data.local_dof_indices[i]] += pres_bc[copy_data.local_dof_indices[i]];
1006 * Finally, <code>node_assembly()</code> takes care of all the
1007 * local computations via
WorkStream mechanism. Notice that the choice
1008 * of the quadrature rule here is dictated by the formulation of the
1009 * method. It has to be <code>degree+1</code> points Gauss-Lobatto
1010 *
for the
volume integrals and <code>degree</code>
for the face ones,
1011 * as mentioned in the introduction.
1014 *
template <
int dim>
1015 *
void MultipointMixedDarcyProblem<dim>::node_assembly()
1019 * dof_handler.distribute_dofs(fe);
1021 * std::vector<types::global_dof_index> dofs_per_component (dim+1);
1025 *
QGauss<dim-1> face_quad(degree);
1027 * n_v = dofs_per_component[0];
1028 * n_p = dofs_per_component[dim];
1030 * pres_rhs.reinit(n_p);
1033 * dof_handler.end(),
1035 * &MultipointMixedDarcyProblem::assemble_system_cell,
1036 * &MultipointMixedDarcyProblem::copy_cell_to_node,
1037 * DataStructures::NodeAssemblyScratchData<dim>(fe,
triangulation,quad,face_quad),
1038 * DataStructures::NodeAssemblyCopyData<dim>());
1044 * <a name=
"Makingthesparsitypattern"></a>
1045 * <h4>Making the sparsity pattern</h4>
1049 * Having computed all the local contributions, we actually have
1050 * all the information needed to make a cell-centered sparsity
1052 * leads to a slower solution.
1055 *
template <
int dim>
1056 *
void MultipointMixedDarcyProblem<dim>::make_cell_centered_sp()
1061 * std::set<types::global_dof_index>::iterator pi_it, pj_it;
1062 *
unsigned int i, j;
1063 *
for (
auto el : node_matrix)
1064 *
for (pi_it = pressure_indices[el.first].begin(), i = 0;
1065 * pi_it != pressure_indices[el.first].end();
1067 *
for (pj_it = pi_it, j = 0;
1068 * pj_it != pressure_indices[el.first].end();
1070 * dsp.add(*pi_it - n_v, *pj_it - n_v);
1074 * cell_centered_sp.copy_from(dsp);
1075 * pres_system_matrix.reinit (cell_centered_sp);
1082 * <a name=
"Thelocaleliminationprocedure"></a>
1083 * <h4>The local elimination procedure</h4>
1087 * This
function finally performs the local elimination procedure.
1088 * Mathematically, it follows the same idea as in computing the
1089 * Schur complement (as mentioned in the introduction) but we
do
1090 * so locally. Namely, local velocity DOFs are expressed in terms
1091 * of corresponding pressure values, and then used
for the local
1095 *
template <
int dim>
1096 *
void MultipointMixedDarcyProblem<dim>::
1097 * nodal_elimination(
const typename DataStructures::PointToMatrixMap<dim>::iterator &n_it,
1098 * DataStructures::VertexEliminationScratchData &scratch_data,
1099 * DataStructures::NodeEliminationCopyData<dim> ©_data)
1101 *
unsigned int n_edges = velocity_indices.at((*n_it).first).size();
1102 *
unsigned int n_cells = pressure_indices.at((*n_it).first).size();
1104 * scratch_data.velocity_matrix.reinit(n_edges,n_edges);
1105 * copy_data.pressure_matrix.reinit(n_edges,
n_cells);
1107 * copy_data.velocity_rhs.reinit(n_edges);
1108 * scratch_data.pressure_rhs.reinit(
n_cells);
1111 * std::set<types::global_dof_index>::iterator vi_it, vj_it, p_it;
1113 *
for (vi_it = velocity_indices.at((*n_it).first).begin(), i = 0;
1114 * vi_it != velocity_indices.at((*n_it).first).end();
1118 *
for (vj_it = velocity_indices.at((*n_it).first).begin(), j = 0;
1119 * vj_it != velocity_indices.at((*n_it).first).end();
1122 * scratch_data.velocity_matrix.add(i, j, node_matrix[(*n_it).first][std::make_pair(*vi_it, *vj_it)]);
1124 * scratch_data.velocity_matrix.add(j, i, node_matrix[(*n_it).first][std::make_pair(*vi_it, *vj_it)]);
1127 *
for (p_it = pressure_indices.at((*n_it).first).begin(), j = 0;
1128 * p_it != pressure_indices.at((*n_it).first).end();
1130 * copy_data.pressure_matrix.add(i, j, node_matrix[(*n_it).first][std::make_pair(*vi_it, *p_it)]);
1132 * copy_data.velocity_rhs(i) += node_rhs.at((*n_it).first)[*vi_it];
1135 *
for (p_it = pressure_indices.at((*n_it).first).begin(), i = 0;
1136 * p_it != pressure_indices.at((*n_it).first).end();
1138 * scratch_data.pressure_rhs(i) += node_rhs.at((*n_it).first)[*p_it];
1141 * copy_data.Ainverse.reinit(n_edges,n_edges);
1143 * scratch_data.tmp_rhs1.reinit(n_edges);
1144 * scratch_data.tmp_rhs2.reinit(n_edges);
1145 * scratch_data.tmp_rhs3.reinit(
n_cells);
1147 * copy_data.Ainverse.invert(scratch_data.velocity_matrix);
1149 * copy_data.node_pres_rhs = scratch_data.pressure_rhs;
1151 * copy_data.node_pres_matrix = 0;
1152 * copy_data.node_pres_matrix.triple_product(copy_data.Ainverse,
1153 * copy_data.pressure_matrix,
1154 * copy_data.pressure_matrix,
true,
false);
1156 * copy_data.Ainverse.vmult(scratch_data.tmp_rhs1, copy_data.velocity_rhs,
false);
1157 * copy_data.pressure_matrix.Tvmult(scratch_data.tmp_rhs3, scratch_data.tmp_rhs1,
false);
1158 * copy_data.node_pres_rhs *= -1.0;
1159 * copy_data.node_pres_rhs += scratch_data.tmp_rhs3;
1161 * copy_data.p = (*n_it).first;
1167 * Each node
's pressure system is then distributed to a global pressure
1168 * system, using the indices we computed in the previous stages.
1171 * template <int dim>
1172 * void MultipointMixedDarcyProblem<dim>::
1173 * copy_node_to_system(const DataStructures::NodeEliminationCopyData<dim> ©_data)
1175 * A_inverse[copy_data.p] = copy_data.Ainverse;
1176 * pressure_matrix[copy_data.p] = copy_data.pressure_matrix;
1177 * velocity_rhs[copy_data.p] = copy_data.velocity_rhs;
1180 * std::set<types::global_dof_index>::iterator pi_it, pj_it;
1182 * for (pi_it = pressure_indices[copy_data.p].begin(), i = 0;
1183 * pi_it != pressure_indices[copy_data.p].end();
1187 * for (pj_it = pressure_indices[copy_data.p].begin(), j = 0;
1188 * pj_it != pressure_indices[copy_data.p].end();
1190 * pres_system_matrix.add(*pi_it - n_v, *pj_it - n_v, copy_data.node_pres_matrix(i, j));
1192 * pres_rhs(*pi_it - n_v) += copy_data.node_pres_rhs(i);
1200 * The @ref WorkStream mechanism is again used for the assembly
1201 * of the global system for the pressure variable, where the
1202 * previous functions are used to perform local computations.
1205 * template <int dim>
1206 * void MultipointMixedDarcyProblem<dim>::pressure_assembly()
1208 * TimerOutput::Scope t(computing_timer, "Pressure matrix assembly");
1210 * QGaussLobatto<dim> quad(degree+1);
1211 * QGauss<dim-1> face_quad(degree);
1213 * pres_rhs.reinit(n_p);
1215 * WorkStream::run(node_matrix.begin(),
1216 * node_matrix.end(),
1218 * &MultipointMixedDarcyProblem::nodal_elimination,
1219 * &MultipointMixedDarcyProblem::copy_node_to_system,
1220 * DataStructures::VertexEliminationScratchData(),
1221 * DataStructures::NodeEliminationCopyData<dim>());
1229 * <a name="Velocitysolutionrecovery"></a>
1230 * <h4>Velocity solution recovery</h4>
1234 * After solving for the pressure variable, we want to follow
1235 * the above procedure backwards, in order to obtain the
1236 * velocity solution (again, this is similar in nature to the
1237 * Schur complement approach, see @ref step_20 "step-20", but here it is done
1238 * locally at each node). We have almost everything computed and
1239 * stored already, including inverses of local mass matrices,
1240 * so the following is a relatively straightforward implementation.
1243 * template <int dim>
1244 * void MultipointMixedDarcyProblem<dim>::
1245 * velocity_assembly (const typename DataStructures::PointToMatrixMap<dim>::iterator &n_it,
1246 * DataStructures::VertexEliminationScratchData &scratch_data,
1247 * DataStructures::NodeEliminationCopyData<dim> ©_data)
1249 * unsigned int n_edges = velocity_indices.at((*n_it).first).size();
1250 * unsigned int n_cells = pressure_indices.at((*n_it).first).size();
1252 * scratch_data.tmp_rhs1.reinit(n_edges);
1253 * scratch_data.tmp_rhs2.reinit(n_edges);
1254 * scratch_data.tmp_rhs3.reinit(n_cells);
1255 * scratch_data.local_pressure_solution.reinit(n_cells);
1257 * copy_data.vertex_vel_solution.reinit(n_edges);
1259 * std::set<types::global_dof_index>::iterator p_it;
1262 * for (p_it = pressure_indices[(*n_it).first].begin(), i = 0;
1263 * p_it != pressure_indices[(*n_it).first].end();
1265 * scratch_data.local_pressure_solution(i) = pres_solution(*p_it - n_v);
1267 * pressure_matrix[(*n_it).first].vmult(scratch_data.tmp_rhs2, scratch_data.local_pressure_solution, false);
1268 * scratch_data.tmp_rhs2 *= -1.0;
1269 * scratch_data.tmp_rhs2+=velocity_rhs[(*n_it).first];
1270 * A_inverse[(*n_it).first].vmult(copy_data.vertex_vel_solution, scratch_data.tmp_rhs2, false);
1272 * copy_data.p = (*n_it).first;
1278 * Copy nodal velocities to a global solution vector by using
1279 * local computations and indices from early stages.
1282 * template <int dim>
1283 * void MultipointMixedDarcyProblem<dim>::
1284 * copy_node_velocity_to_global(const DataStructures::NodeEliminationCopyData<dim> ©_data)
1286 * std::set<types::global_dof_index>::iterator vi_it;
1289 * for (vi_it = velocity_indices[copy_data.p].begin(), i = 0;
1290 * vi_it != velocity_indices[copy_data.p].end();
1292 * vel_solution(*vi_it) += copy_data.vertex_vel_solution(i);
1298 * Use @ref WorkStream to run everything concurrently.
1301 * template <int dim>
1302 * void MultipointMixedDarcyProblem<dim>::velocity_recovery()
1304 * TimerOutput::Scope t(computing_timer, "Velocity solution recovery");
1306 * QGaussLobatto<dim> quad(degree+1);
1307 * QGauss<dim-1> face_quad(degree);
1309 * vel_solution.reinit(n_v);
1311 * WorkStream::run(node_matrix.begin(),
1312 * node_matrix.end(),
1314 * &MultipointMixedDarcyProblem::velocity_assembly,
1315 * &MultipointMixedDarcyProblem::copy_node_velocity_to_global,
1316 * DataStructures::VertexEliminationScratchData(),
1317 * DataStructures::NodeEliminationCopyData<dim>());
1319 * solution.reinit(2);
1320 * solution.block(0) = vel_solution;
1321 * solution.block(1) = pres_solution;
1322 * solution.collect_sizes();
1330 * <a name="Pressuresystemsolver"></a>
1331 * <h4>Pressure system solver</h4>
1335 * The solver part is trivial. We use the CG solver with no
1336 * preconditioner for simplicity.
1339 * template <int dim>
1340 * void MultipointMixedDarcyProblem<dim>::solve_pressure()
1342 * TimerOutput::Scope t(computing_timer, "Pressure CG solve");
1344 * pres_solution.reinit(n_p);
1346 * SolverControl solver_control (2.0*n_p, 1e-10);
1347 * SolverCG<> solver (solver_control);
1349 * PreconditionIdentity identity;
1350 * solver.solve(pres_system_matrix, pres_solution, pres_rhs, identity);
1358 * <a name="Postprocessing"></a>
1359 * <h3>Postprocessing</h3>
1363 * We have two postprocessing steps here, first one computes the
1364 * errors in order to populate the convergence tables. The other
1365 * one takes care of the output of the solutions in <code>.vtk</code>
1371 * <a name="Computeerrors"></a>
1372 * <h4>Compute errors</h4>
1376 * The implementation of this function is almost identical to @ref step_20 "step-20".
1377 * We use @ref ComponentSelectFunction as masks to use the right
1378 * solution component (velocity or pressure) and @ref integrate_difference
1379 * to compute the errors. Since we also want to compute Hdiv seminorm of the
1380 * velocity error, one must provide gradients in the <code>ExactSolution</code>
1381 * class implementation to avoid exceptions. The only noteworthy thing here
1382 * is that we again use lower order quadrature rule instead of projecting the
1383 * solution to an appropriate space in order to show superconvergence, which is
1384 * mathematically justified.
1387 * template <int dim>
1388 * void MultipointMixedDarcyProblem<dim>::compute_errors(const unsigned cycle)
1390 * TimerOutput::Scope t(computing_timer, "Compute errors");
1392 * const ComponentSelectFunction<dim> pressure_mask(dim, dim+1);
1393 * const ComponentSelectFunction<dim> velocity_mask(std::make_pair(0, dim), dim+1);
1395 * ExactSolution<dim> exact_solution;
1397 * Vector<double> cellwise_errors (triangulation.n_active_cells());
1399 * QTrapez<1> q_trapez;
1400 * QIterated<dim> quadrature(q_trapez,degree+2);
1401 * QGauss<dim> quadrature_super(degree);
1403 * VectorTools::integrate_difference (dof_handler, solution, exact_solution,
1404 * cellwise_errors, quadrature,
1405 * VectorTools::L2_norm,
1407 * const double p_l2_error = cellwise_errors.l2_norm();
1409 * VectorTools::integrate_difference (dof_handler, solution, exact_solution,
1410 * cellwise_errors, quadrature_super,
1411 * VectorTools::L2_norm,
1413 * const double p_l2_mid_error = cellwise_errors.l2_norm();
1415 * VectorTools::integrate_difference (dof_handler, solution, exact_solution,
1416 * cellwise_errors, quadrature,
1417 * VectorTools::L2_norm,
1419 * const double u_l2_error = cellwise_errors.l2_norm();
1421 * VectorTools::integrate_difference (dof_handler, solution, exact_solution,
1422 * cellwise_errors, quadrature,
1423 * VectorTools::Hdiv_seminorm,
1425 * const double u_hd_error = cellwise_errors.l2_norm();
1427 * const unsigned int n_active_cells=triangulation.n_active_cells();
1428 * const unsigned int n_dofs=dof_handler.n_dofs();
1430 * convergence_table.add_value("cycle", cycle);
1431 * convergence_table.add_value("cells", n_active_cells);
1432 * convergence_table.add_value("dofs", n_dofs);
1433 * convergence_table.add_value("Velocity,L2", u_l2_error);
1434 * convergence_table.add_value("Velocity,Hdiv", u_hd_error);
1435 * convergence_table.add_value("Pressure,L2", p_l2_error);
1436 * convergence_table.add_value("Pressure,L2-nodal", p_l2_mid_error);
1444 * <a name="Outputresults"></a>
1445 * <h4>Output results</h4>
1449 * This function also follows the same idea as in @ref step_20 "step-20" tutorial
1450 * program. The only modification to it is the part involving
1451 * a convergence table.
1454 * template <int dim>
1455 * void MultipointMixedDarcyProblem<dim>::output_results(const unsigned int cycle, const unsigned int refine)
1457 * TimerOutput::Scope t(computing_timer, "Output results");
1459 * std::vector<std::string> solution_names(dim, "u");
1460 * solution_names.push_back ("p");
1461 * std::vector<DataComponentInterpretation::DataComponentInterpretation>
1462 * interpretation (dim, DataComponentInterpretation::component_is_part_of_vector);
1463 * interpretation.push_back (DataComponentInterpretation::component_is_scalar);
1465 * DataOut<dim> data_out;
1466 * data_out.add_data_vector (dof_handler, solution, solution_names, interpretation);
1467 * data_out.build_patches ();
1469 * std::ofstream output ("solution" + std::to_string(dim) + "d-" + std::to_string(cycle) + ".vtk");
1470 * data_out.write_vtk (output);
1472 * convergence_table.set_precision("Velocity,L2", 3);
1473 * convergence_table.set_precision("Velocity,Hdiv", 3);
1474 * convergence_table.set_precision("Pressure,L2", 3);
1475 * convergence_table.set_precision("Pressure,L2-nodal", 3);
1476 * convergence_table.set_scientific("Velocity,L2", true);
1477 * convergence_table.set_scientific("Velocity,Hdiv", true);
1478 * convergence_table.set_scientific("Pressure,L2", true);
1479 * convergence_table.set_scientific("Pressure,L2-nodal", true);
1480 * convergence_table.set_tex_caption("cells", "\\# cells");
1481 * convergence_table.set_tex_caption("dofs", "\\# dofs");
1482 * convergence_table.set_tex_caption("Velocity,L2", " \\|\\u - \\u_h\\|_{L^2} ");
1483 * convergence_table.set_tex_caption("Velocity,Hdiv", " \\|\\nabla\\cdot(\\u - \\u_h)\\|_{L^2} ");
1484 * convergence_table.set_tex_caption("Pressure,L2", " \\|p - p_h\\|_{L^2} ");
1485 * convergence_table.set_tex_caption("Pressure,L2-nodal", " \\|Qp - p_h\\|_{L^2} ");
1486 * convergence_table.set_tex_format("cells", "r");
1487 * convergence_table.set_tex_format("dofs", "r");
1489 * convergence_table.evaluate_convergence_rates("Velocity,L2", ConvergenceTable::reduction_rate_log2);
1490 * convergence_table.evaluate_convergence_rates("Velocity,Hdiv", ConvergenceTable::reduction_rate_log2);
1491 * convergence_table.evaluate_convergence_rates("Pressure,L2", ConvergenceTable::reduction_rate_log2);
1492 * convergence_table.evaluate_convergence_rates("Pressure,L2-nodal", ConvergenceTable::reduction_rate_log2);
1494 * std::ofstream error_table_file("error" + std::to_string(dim) + "d.tex");
1496 * if (cycle == refine-1)
1498 * convergence_table.write_text(std::cout);
1499 * convergence_table.write_tex(error_table_file);
1508 * <a name="Runfunction"></a>
1509 * <h3>Run function</h3>
1513 * The driver method <code>run()</code>
1514 * takes care of mesh generation and arranging calls to member methods in
1515 * the right way. It also resets data structures and clear triangulation and
1516 * DOF handler as we run the method on a sequence of refinements in order
1517 * to record convergence rates.
1520 * template <int dim>
1521 * void MultipointMixedDarcyProblem<dim>::run(const unsigned int refine)
1523 * Assert(refine > 0, ExcMessage("Must at least have 1 refinement cycle!"));
1525 * dof_handler.clear();
1526 * triangulation.clear();
1527 * convergence_table.clear();
1529 * for (unsigned int cycle=0; cycle<refine; ++cycle)
1535 * We first generate the hyper cube and refine it twice
1536 * so that we could distort the grid slightly and
1537 * demonstrate the method's ability to work in such a
1549 * make_cell_centered_sp();
1550 * pressure_assembly();
1551 * solve_pressure ();
1552 * velocity_recovery ();
1553 * compute_errors (cycle);
1554 * output_results (cycle,
refine);
1555 * reset_data_structures ();
1557 * computing_timer.print_summary ();
1558 * computing_timer.reset ();
1567 * <a name=
"Thecodemaincodefunction"></a>
1568 * <h3>The <code>main</code>
function</h3>
1572 * In the main functione we pass the order of the Finite Element as an argument
1573 * to the constructor of the Multipoint Flux Mixed Darcy problem, and the number
1574 * of refinement cycles as an argument
for the
run method.
1581 *
using namespace dealii;
1582 *
using namespace MFMFE;
1586 * MultipointMixedDarcyProblem<2> mfmfe_problem(2);
1587 * mfmfe_problem.run(6);
1589 *
catch (std::exception &exc)
1591 * std::cerr << std::endl << std::endl
1592 * <<
"----------------------------------------------------"
1594 * std::cerr <<
"Exception on processing: " << std::endl
1595 * << exc.what() << std::endl
1596 * <<
"Aborting!" << std::endl
1597 * <<
"----------------------------------------------------"
1604 * std::cerr << std::endl << std::endl
1605 * <<
"----------------------------------------------------"
1607 * std::cerr <<
"Unknown exception!" << std::endl
1608 * <<
"Aborting!" << std::endl
1609 * <<
"----------------------------------------------------"