Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
MultipointFluxMixedFiniteElementMethods.h
Go to the documentation of this file.
1 
164  *
165  * #ifndef MFMFE_DATA_H
166  * #define MFMFE_DATA_H
167  *
168  * #include <deal.II/base/function.h>
169  * #include <deal.II/base/tensor_function.h>
170  *
171  * @endcode
172  *
173  *
174  * <a name="Dataandexactsolution"></a>
175  * <h3>Data and exact solution.</h3>
176  *
177 
178  *
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 2d cases are
182  * provided, but 3d can be added straightforwardly.
183  *
184 
185  *
186  *
187  * @code
188  * namespace MFMFE
189  * {
190  * using namespace dealii;
191  *
192  * template <int dim>
193  * class RightHandSide : public Function<dim>
194  * {
195  * public:
196  * RightHandSide () : Function<dim>(1) {}
197  *
198  * virtual double value (const Point<dim> &p,
199  * const unsigned int component = 0) const;
200  * };
201  *
202  * template <int dim>
203  * double RightHandSide<dim>::value (const Point<dim> &p,
204  * const unsigned int /*component*/) const
205  * {
206  * const double x = p[0];
207  * const double y = p[1];
208  *
209  * switch (dim)
210  * {
211  * case 2:
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);
216  * default:
217  * Assert(false, ExcMessage("The RHS data for dim != 2 is not provided"));
218  * }
219  * }
220  *
221  *
222  *
223  * template <int dim>
224  * class PressureBoundaryValues : public Function<dim>
225  * {
226  * public:
227  * PressureBoundaryValues () : Function<dim>(1) {}
228  *
229  * virtual double value (const Point<dim> &p,
230  * const unsigned int component = 0) const;
231  * };
232  *
233  * template <int dim>
235  * const unsigned int /*component*/) const
236  * {
237  * const double x = p[0];
238  * const double y = p[1];
239  *
240  * switch (dim)
241  * {
242  * case 2:
243  * return (x*x*x)*(y*y*y*y)+cos(x*y)*sin(x*y)+x*x;
244  * default:
245  * Assert(false, ExcMessage("The BC data for dim != 2 is not provided"));
246  * }
247  * }
248  *
249  *
250  *
251  * template <int dim>
252  * class ExactSolution : public Function<dim>
253  * {
254  * public:
255  * ExactSolution () : Function<dim>(dim+1) {}
256  *
257  * virtual void vector_value (const Point<dim> &p,
258  * Vector<double> &value) const;
259  *
260  * virtual void vector_gradient (const Point<dim> &p,
261  * std::vector<Tensor<1,dim,double>> &grads) const;
262  * };
263  *
264  * template <int dim>
265  * void
266  * ExactSolution<dim>::vector_value (const Point<dim> &p,
267  * Vector<double> &values) const
268  * {
269  * Assert (values.size() == dim+1,
270  * ExcDimensionMismatch (values.size(), dim+1));
271  *
272  * const double x = p[0];
273  * const double y = p[1];
274  *
275  * switch (dim)
276  * {
277  * case 2:
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;
281  * break;
282  * default:
283  * Assert(false, ExcMessage("The exact solution for dim != 2 is not provided"));
284  * }
285  * }
286  *
287  * template <int dim>
288  * void
289  * ExactSolution<dim>::vector_gradient (const Point<dim> &p,
290  * std::vector<Tensor<1,dim,double>> &grads) const
291  * {
292  * const double x = p[0];
293  * const double y = p[1];
294  *
295  * switch (dim)
296  * {
297  * case 2:
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));
310  * grads[2] = 0;
311  * break;
312  * default:
313  * Assert(false, ExcMessage("The exact solution's gradient for dim != 2 is not provided"));
314  * }
315  * }
316  *
317  *
318  *
319  * template <int dim>
320  * class KInverse : public TensorFunction<2,dim>
321  * {
322  * public:
323  * KInverse () : TensorFunction<2,dim>() {}
324  *
325  * virtual void value_list (const std::vector<Point<dim> > &points,
326  * std::vector<Tensor<2,dim> > &values) const;
327  * };
328  *
329  * template <int dim>
330  * void
331  * KInverse<dim>::value_list (const std::vector<Point<dim> > &points,
332  * std::vector<Tensor<2,dim> > &values) const
333  * {
334  * Assert (points.size() == values.size(),
335  * ExcDimensionMismatch (points.size(), values.size()));
336  *
337  * for (unsigned int p=0; p<points.size(); ++p)
338  * {
339  * values[p].clear ();
340  *
341  * const double x = points[p][0];
342  * const double y = points[p][1];
343  *
344  * switch (dim)
345  * {
346  * case 2:
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);
351  * break;
352  * default:
353  * Assert(false, ExcMessage("The inverse of permeability tensor for dim != 2 is not provided"));
354  * }
355  * }
356  * }
357  * }
358  *
359  * #endif //MFMFE_DATA_H
360  * @endcode
361 
362 
363 <a name="ann-mfmfe.cc"></a>
364 <h1>Annotated version of mfmfe.cc</h1>
365  *
366  *
367  *
368  *
369  * @code
370  * /* ---------------------------------------------------------------------
371  * *
372  * * This file is part of the deal.II Code Gallery.
373  * *
374  * * ---------------------------------------------------------------------
375  * *
376  * * Author: Ilona Ambartsumyan, Eldar Khattatov, University of Pittsburgh, 2018
377  * */
378  *
379  *
380  * @endcode
381  *
382  *
383  * <a name="Includefiles"></a>
384  * <h3>Include files</h3>
385  *
386 
387  *
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.
391  *
392  * @code
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>
398  *
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>
404  *
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>
417  *
418  * #include <fstream>
419  * #include <unordered_map>
420  *
421  * @endcode
422  *
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.
426  *
427  * @code
428  * #include <deal.II/fe/fe_rt_bubbles.h>
429  *
430  * @endcode
431  *
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
435  * here
436  *
437  * @code
438  * #include "data.h"
439  *
440  * @endcode
441  *
442  * As always the program is in the namespace of its own with
443  * the deal.II classes and functions imported into it
444  *
445  * @code
446  * namespace MFMFE
447  * {
448  * using namespace dealii;
449  *
450  * @endcode
451  *
452  *
453  * <a name="Definitionofmultipointfluxassemblydatastructures"></a>
454  * <h3>Definition of multipoint flux assembly data structures</h3>
455  *
456 
457  *
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$.
463  *
464  * @code
465  * namespace DataStructures
466  * {
467  * @endcode
468  *
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
472  * the hash table, with the keys as Point<dim>
473  *
474  * @code
475  * template <int dim>
476  * struct hash_points
477  * {
478  * size_t operator()(const Point<dim> &p) const
479  * {
480  * size_t h1,h2,h3;
481  * h1 = std::hash<double>()(p[0]);
482  *
483  * switch (dim)
484  * {
485  * case 1:
486  * return h1;
487  * case 2:
488  * h2 = std::hash<double>()(p[1]);
489  * return (h1 ^ h2);
490  * case 3:
491  * h2 = std::hash<double>()(p[1]);
492  * h3 = std::hash<double>()(p[2]);
493  * return (h1 ^ (h2 << 1)) ^ h3;
494  * default:
495  * Assert(false, ExcNotImplemented());
496  * }
497  * }
498  * };
499  *
500  * @endcode
501  *
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
504  *
505  * @code
506  * template <int dim>
507  * using PointToMatrixMap = std::unordered_map<Point<dim>, std::map<std::pair<types::global_dof_index,types::global_dof_index>, double>, hash_points<dim>>;
508  *
509  * template <int dim>
510  * using PointToVectorMap = std::unordered_map<Point<dim>, std::map<types::global_dof_index, double>, hash_points<dim>>;
511  *
512  * template <int dim>
513  * using PointToIndexMap = std::unordered_map<Point<dim>, std::set<types::global_dof_index>, hash_points<dim>>;
514  *
515  * @endcode
516  *
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...
522  *
523  * @code
524  * template <int dim>
525  * struct NodeAssemblyCopyData
526  * {
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;
532  * };
533  *
534  * @endcode
535  *
536  * ... and the other one for the actual process of
537  * local velocity elimination and assembling the global
538  * pressure system:
539  *
540  * @code
541  * template <int dim>
542  * struct NodeEliminationCopyData
543  * {
544  * FullMatrix<double> node_pres_matrix;
545  * Vector<double> node_pres_rhs;
546  * FullMatrix<double> Ainverse;
547  * FullMatrix<double> pressure_matrix;
548  * Vector<double> velocity_rhs;
549  * Vector<double> vertex_vel_solution;
550  * Point<dim> p;
551  * };
552  *
553  * @endcode
554  *
555  * Similarly, two ScratchData classes are defined.
556  * One for the assembly part, where we need
557  * FEValues, FEFaceValues, Quadrature and storage
558  * for the basis fuctions...
559  *
560  * @code
561  * template <int dim>
562  * struct NodeAssemblyScratchData
563  * {
564  * NodeAssemblyScratchData (const FiniteElement<dim> &fe,
565  * const Triangulation<dim> &tria,
566  * const Quadrature<dim> &quad,
567  * const Quadrature<dim-1> &f_quad);
568  *
569  * NodeAssemblyScratchData (const NodeAssemblyScratchData &scratch_data);
570  *
571  * FEValues<dim> fe_values;
572  * FEFaceValues<dim> fe_face_values;
573  * std::vector<unsigned int> n_faces_at_vertex;
574  *
575  * const unsigned long num_cells;
576  *
577  * std::vector<Tensor<2,dim>> k_inverse_values;
578  * std::vector<double> rhs_values;
579  * std::vector<double> pres_bc_values;
580  *
581  * std::vector<Tensor<1,dim> > phi_u;
582  * std::vector<double> div_phi_u;
583  * std::vector<double> phi_p;
584  * };
585  *
586  * template <int dim>
587  * NodeAssemblyScratchData<dim>::
588  * NodeAssemblyScratchData (const FiniteElement<dim> &fe,
589  * const Triangulation<dim> &tria,
590  * const Quadrature<dim> &quad,
591  * const Quadrature<dim-1> &f_quad)
592  * :
593  * fe_values (fe,
594  * quad,
597  * fe_face_values (fe,
598  * f_quad,
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)
608  * {
609  * n_faces_at_vertex.resize(tria.n_vertices(), 0);
610  * typename Triangulation<dim>::active_face_iterator face = tria.begin_active_face(), endf = tria.end_face();
611  *
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;
615  * }
616  *
617  * template <int dim>
618  * NodeAssemblyScratchData<dim>::
619  * NodeAssemblyScratchData (const NodeAssemblyScratchData &scratch_data)
620  * :
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)
637  * {}
638  *
639  * @endcode
640  *
641  * ...and the other, simpler one, for the velocity elimination and recovery
642  *
643  * @code
644  * struct VertexEliminationScratchData
645  * {
646  * VertexEliminationScratchData () = default;
647  * VertexEliminationScratchData (const VertexEliminationScratchData &scratch_data);
648  *
649  * FullMatrix<double> velocity_matrix;
650  * Vector<double> pressure_rhs;
651  *
652  * Vector<double> local_pressure_solution;
653  * Vector<double> tmp_rhs1;
654  * Vector<double> tmp_rhs2;
655  * Vector<double> tmp_rhs3;
656  * };
657  *
658  * VertexEliminationScratchData::
659  * VertexEliminationScratchData (const VertexEliminationScratchData &scratch_data)
660  * :
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)
667  * {}
668  * }
669  *
670  *
671  *
672  * @endcode
673  *
674  *
675  * <a name="ThecodeMultipointMixedDarcyProblemcodeclasstemplate"></a>
676  * <h3>The <code>MultipointMixedDarcyProblem</code> class template</h3>
677  *
678 
679  *
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.
685  *
686  * @code
687  * template <int dim>
688  * class MultipointMixedDarcyProblem
689  * {
690  * public:
691  * MultipointMixedDarcyProblem (const unsigned int degree);
692  * ~MultipointMixedDarcyProblem ();
693  * void run (const unsigned int refine);
694  * private:
695  * void assemble_system_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
696  * DataStructures::NodeAssemblyScratchData<dim> &scratch_data,
697  * DataStructures::NodeAssemblyCopyData<dim> &copy_data);
698  * void copy_cell_to_node(const DataStructures::NodeAssemblyCopyData<dim> &copy_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> &copy_data);
704  * void copy_node_to_system(const DataStructures::NodeEliminationCopyData<dim> &copy_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> &copy_data);
710  * void copy_node_velocity_to_global(const DataStructures::NodeEliminationCopyData<dim> &copy_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);
715  *
716  * const unsigned int degree;
718  * FESystem<dim> fe;
719  * DoFHandler<dim> dof_handler;
720  * BlockVector<double> solution;
721  *
722  * SparsityPattern cell_centered_sp;
723  * SparseMatrix<double> pres_system_matrix;
724  * Vector<double> pres_rhs;
725  *
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;
729  *
730  * DataStructures::PointToMatrixMap<dim> node_matrix;
731  * DataStructures::PointToVectorMap<dim> node_rhs;
732  *
733  * DataStructures::PointToIndexMap<dim> pressure_indices;
734  * DataStructures::PointToIndexMap<dim> velocity_indices;
735  *
736  * unsigned long n_v, n_p;
737  *
738  * Vector<double> pres_solution;
739  * Vector<double> vel_solution;
740  *
741  * ConvergenceTable convergence_table;
742  * TimerOutput computing_timer;
743  * };
744  *
745  * @endcode
746  *
747  *
748  * <a name="Constructoranddestructorcodereset_data_structurescode"></a>
749  * <h4>Constructor and destructor, <code>reset_data_structures</code></h4>
750  *
751 
752  *
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
755  * degree of one would mean the use of @ref FE_RT_Bubbles(1) and @ref FE_DGQ(0)),
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.
759  *
760  * @code
761  * template <int dim>
762  * MultipointMixedDarcyProblem<dim>::MultipointMixedDarcyProblem (const unsigned int degree)
763  * :
764  * degree(degree),
765  * fe(FE_RT_Bubbles<dim>(degree), 1,
766  * FE_DGQ<dim>(degree-1), 1),
767  * dof_handler(triangulation),
768  * computing_timer(std::cout, TimerOutput::summary,
770  * {}
771  *
772  *
773  * @endcode
774  *
775  * The destructor clears the <code>dof_handler</code> and
776  * all of the data structures we used for the method.
777  *
778  * @code
779  * template <int dim>
780  * MultipointMixedDarcyProblem<dim>::~MultipointMixedDarcyProblem()
781  * {
782  * reset_data_structures ();
783  * dof_handler.clear();
784  * }
785  *
786  *
787  * @endcode
788  *
789  * This method clears all the data that was used after one refinement
790  * cycle.
791  *
792  * @code
793  * template <int dim>
794  * void MultipointMixedDarcyProblem<dim>::reset_data_structures ()
795  * {
796  * pressure_indices.clear();
797  * velocity_indices.clear();
798  * velocity_rhs.clear();
799  * A_inverse.clear();
800  * pressure_matrix.clear();
801  * node_matrix.clear();
802  * node_rhs.clear();
803  * }
804  *
805  *
806  * @endcode
807  *
808  *
809  * <a name="Cellwiseassemblyandcreationofthelocalnodalbaseddatastructures"></a>
810  * <h4>Cell-wise assembly and creation of the local, nodal-based data structures</h4>
811  *
812 
813  *
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.
817  *
818  * @code
819  * template <int dim>
820  * void MultipointMixedDarcyProblem<dim>::copy_cell_to_node(const DataStructures::NodeAssemblyCopyData<dim> &copy_data)
821  * {
822  * for (auto m : copy_data.cell_mat)
823  * {
824  * for (auto p : m.second)
825  * node_matrix[m.first][p.first] += p.second;
826  *
827  * for (auto p : copy_data.cell_vec.at(m.first))
828  * node_rhs[m.first][p.first] += p.second;
829  *
830  * for (auto p : copy_data.local_pres_indices.at(m.first))
831  * pressure_indices[m.first].insert(p);
832  *
833  * for (auto p : copy_data.local_vel_indices.at(m.first))
834  * velocity_indices[m.first].insert(p);
835  * }
836  * }
837  *
838  *
839  *
840  * @endcode
841  *
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.
845  *
846  * @code
847  * template <int dim>
848  * void MultipointMixedDarcyProblem<dim>::
849  * assemble_system_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
850  * DataStructures::NodeAssemblyScratchData<dim> &scratch_data,
851  * DataStructures::NodeAssemblyCopyData<dim> &copy_data)
852  * {
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();
857  *
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();
861  *
862  * copy_data.local_dof_indices.resize(dofs_per_cell);
863  * cell->get_dof_indices (copy_data.local_dof_indices);
864  *
865  * scratch_data.fe_values.reinit (cell);
866  *
867  * const KInverse<dim> k_inverse;
868  * const RightHandSide<dim> rhs;
869  * const PressureBoundaryValues<dim> pressure_bc;
870  *
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);
873  *
874  * const FEValuesExtractors::Vector velocity (0);
875  * const FEValuesExtractors::Scalar pressure (dim);
876  *
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;
879  *
880  * @endcode
881  *
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
886  *
887  * @code
888  * for (unsigned int q=0; q<n_q_points; ++q)
889  * {
890  * const Point<dim> p = scratch_data.fe_values.quadrature_point(q);
891  *
892  * for (unsigned int k=0; k<dofs_per_cell; ++k)
893  * {
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);
897  * }
898  *
899  * for (unsigned int i=0; i<dofs_per_cell; ++i)
900  * {
901  * for (unsigned int j=n_vel; j<dofs_per_cell; ++j)
902  * {
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);
905  *
906  * if (std::abs(div_term) > 1.e-12)
907  * div_map[i][j] += div_term;
908  * }
909  *
910  * double source_term = -scratch_data.phi_p[i] * scratch_data.rhs_values[q] * scratch_data.fe_values.JxW(q);
911  *
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;
914  * }
915  * }
916  *
917  * @endcode
918  *
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.
924  *
925  * @code
926  * for (unsigned int q=0; q<n_q_points; ++q)
927  * {
928  * std::set<types::global_dof_index> vel_indices;
929  * const Point<dim> p = scratch_data.fe_values.quadrature_point(q);
930  *
931  * for (unsigned int k=0; k<dofs_per_cell; ++k)
932  * {
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);
936  * }
937  *
938  * for (unsigned int i=0; i<dofs_per_cell; ++i)
939  * for (unsigned int j=i; j<dofs_per_cell; ++j)
940  * {
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);
945  *
946  * if (std::abs(mass_term) > 1.e-12)
947  * {
948  * copy_data.cell_mat[p][std::make_pair(copy_data.local_dof_indices[i], copy_data.local_dof_indices[j])] +=
949  * mass_term;
950  * vel_indices.insert(i);
951  * copy_data.local_vel_indices[p].insert(copy_data.local_dof_indices[j]);
952  * }
953  * }
954  *
955  * for (auto i : vel_indices)
956  * for (auto el : div_map[i])
957  * if (std::abs(el.second) > 1.e-12)
958  * {
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]);
962  * }
963  * }
964  *
965  * @endcode
966  *
967  * The pressure boundary conditions are computed as in @ref step_20 "step-20",
968  *
969  * @code
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;
973  * ++face_no)
974  * if (cell->at_boundary(face_no))
975  * {
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);
978  *
979  * for (unsigned int q=0; q<n_face_q_points; ++q)
980  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
981  * {
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));
986  *
987  * if (std::abs(tmp) > 1.e-12)
988  * pres_bc[copy_data.local_dof_indices[i]] += tmp;
989  * }
990  * }
991  *
992  * @endcode
993  *
994  * ...but we distribute them to the corresponding nodal data structures
995  *
996  * @code
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]];
1001  * }
1002  *
1003  *
1004  * @endcode
1005  *
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.
1012  *
1013  * @code
1014  * template <int dim>
1015  * void MultipointMixedDarcyProblem<dim>::node_assembly()
1016  * {
1017  * TimerOutput::Scope t(computing_timer, "Nodal assembly");
1018  *
1019  * dof_handler.distribute_dofs(fe);
1020  * DoFRenumbering::component_wise (dof_handler);
1021  * std::vector<types::global_dof_index> dofs_per_component (dim+1);
1022  * DoFTools::count_dofs_per_component (dof_handler, dofs_per_component);
1023  *
1024  * QGaussLobatto<dim> quad(degree+1);
1025  * QGauss<dim-1> face_quad(degree);
1026  *
1027  * n_v = dofs_per_component[0];
1028  * n_p = dofs_per_component[dim];
1029  *
1030  * pres_rhs.reinit(n_p);
1031  *
1032  * WorkStream::run(dof_handler.begin_active(),
1033  * dof_handler.end(),
1034  * *this,
1035  * &MultipointMixedDarcyProblem::assemble_system_cell,
1036  * &MultipointMixedDarcyProblem::copy_cell_to_node,
1037  * DataStructures::NodeAssemblyScratchData<dim>(fe, triangulation,quad,face_quad),
1038  * DataStructures::NodeAssemblyCopyData<dim>());
1039  * }
1040  *
1041  * @endcode
1042  *
1043  *
1044  * <a name="Makingthesparsitypattern"></a>
1045  * <h4>Making the sparsity pattern</h4>
1046  *
1047 
1048  *
1049  * Having computed all the local contributions, we actually have
1050  * all the information needed to make a cell-centered sparsity
1051  * pattern manually. We do this here, because @ref SparseMatrixEZ
1052  * leads to a slower solution.
1053  *
1054  * @code
1055  * template <int dim>
1056  * void MultipointMixedDarcyProblem<dim>::make_cell_centered_sp()
1057  * {
1058  * TimerOutput::Scope t(computing_timer, "Make sparsity pattern");
1059  * DynamicSparsityPattern dsp(n_p, n_p);
1060  *
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();
1066  * ++pi_it, ++i)
1067  * for (pj_it = pi_it, j = 0;
1068  * pj_it != pressure_indices[el.first].end();
1069  * ++pj_it, ++j)
1070  * dsp.add(*pi_it - n_v, *pj_it - n_v);
1071  *
1072  *
1073  * dsp.symmetrize();
1074  * cell_centered_sp.copy_from(dsp);
1075  * pres_system_matrix.reinit (cell_centered_sp);
1076  * }
1077  *
1078  *
1079  * @endcode
1080  *
1081  *
1082  * <a name="Thelocaleliminationprocedure"></a>
1083  * <h4>The local elimination procedure</h4>
1084  *
1085 
1086  *
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
1092  * pressure systems.
1093  *
1094  * @code
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> &copy_data)
1100  * {
1101  * unsigned int n_edges = velocity_indices.at((*n_it).first).size();
1102  * unsigned int n_cells = pressure_indices.at((*n_it).first).size();
1103  *
1104  * scratch_data.velocity_matrix.reinit(n_edges,n_edges);
1105  * copy_data.pressure_matrix.reinit(n_edges,n_cells);
1106  *
1107  * copy_data.velocity_rhs.reinit(n_edges);
1108  * scratch_data.pressure_rhs.reinit(n_cells);
1109  *
1110  * {
1111  * std::set<types::global_dof_index>::iterator vi_it, vj_it, p_it;
1112  * unsigned int i;
1113  * for (vi_it = velocity_indices.at((*n_it).first).begin(), i = 0;
1114  * vi_it != velocity_indices.at((*n_it).first).end();
1115  * ++vi_it, ++i)
1116  * {
1117  * unsigned int j;
1118  * for (vj_it = velocity_indices.at((*n_it).first).begin(), j = 0;
1119  * vj_it != velocity_indices.at((*n_it).first).end();
1120  * ++vj_it, ++j)
1121  * {
1122  * scratch_data.velocity_matrix.add(i, j, node_matrix[(*n_it).first][std::make_pair(*vi_it, *vj_it)]);
1123  * if (j != i)
1124  * scratch_data.velocity_matrix.add(j, i, node_matrix[(*n_it).first][std::make_pair(*vi_it, *vj_it)]);
1125  * }
1126  *
1127  * for (p_it = pressure_indices.at((*n_it).first).begin(), j = 0;
1128  * p_it != pressure_indices.at((*n_it).first).end();
1129  * ++p_it, ++j)
1130  * copy_data.pressure_matrix.add(i, j, node_matrix[(*n_it).first][std::make_pair(*vi_it, *p_it)]);
1131  *
1132  * copy_data.velocity_rhs(i) += node_rhs.at((*n_it).first)[*vi_it];
1133  * }
1134  *
1135  * for (p_it = pressure_indices.at((*n_it).first).begin(), i = 0;
1136  * p_it != pressure_indices.at((*n_it).first).end();
1137  * ++p_it, ++i)
1138  * scratch_data.pressure_rhs(i) += node_rhs.at((*n_it).first)[*p_it];
1139  * }
1140  *
1141  * copy_data.Ainverse.reinit(n_edges,n_edges);
1142  *
1143  * scratch_data.tmp_rhs1.reinit(n_edges);
1144  * scratch_data.tmp_rhs2.reinit(n_edges);
1145  * scratch_data.tmp_rhs3.reinit(n_cells);
1146  *
1147  * copy_data.Ainverse.invert(scratch_data.velocity_matrix);
1148  * copy_data.node_pres_matrix.reinit(n_cells, n_cells);
1149  * copy_data.node_pres_rhs = scratch_data.pressure_rhs;
1150  *
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);
1155  *
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;
1160  *
1161  * copy_data.p = (*n_it).first;
1162  * }
1163  *
1164  *
1165  * @endcode
1166  *
1167  * Each node's pressure system is then distributed to a global pressure
1168  * system, using the indices we computed in the previous stages.
1169  *
1170  * @code
1171  * template <int dim>
1172  * void MultipointMixedDarcyProblem<dim>::
1173  * copy_node_to_system(const DataStructures::NodeEliminationCopyData<dim> &copy_data)
1174  * {
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;
1178  *
1179  * {
1180  * std::set<types::global_dof_index>::iterator pi_it, pj_it;
1181  * unsigned int i;
1182  * for (pi_it = pressure_indices[copy_data.p].begin(), i = 0;
1183  * pi_it != pressure_indices[copy_data.p].end();
1184  * ++pi_it, ++i)
1185  * {
1186  * unsigned int j;
1187  * for (pj_it = pressure_indices[copy_data.p].begin(), j = 0;
1188  * pj_it != pressure_indices[copy_data.p].end();
1189  * ++pj_it, ++j)
1190  * pres_system_matrix.add(*pi_it - n_v, *pj_it - n_v, copy_data.node_pres_matrix(i, j));
1191  *
1192  * pres_rhs(*pi_it - n_v) += copy_data.node_pres_rhs(i);
1193  * }
1194  * }
1195  * }
1196  *
1197  *
1198  * @endcode
1199  *
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.
1203  *
1204  * @code
1205  * template <int dim>
1206  * void MultipointMixedDarcyProblem<dim>::pressure_assembly()
1207  * {
1208  * TimerOutput::Scope t(computing_timer, "Pressure matrix assembly");
1209  *
1210  * QGaussLobatto<dim> quad(degree+1);
1211  * QGauss<dim-1> face_quad(degree);
1212  *
1213  * pres_rhs.reinit(n_p);
1214  *
1215  * WorkStream::run(node_matrix.begin(),
1216  * node_matrix.end(),
1217  * *this,
1218  * &MultipointMixedDarcyProblem::nodal_elimination,
1219  * &MultipointMixedDarcyProblem::copy_node_to_system,
1220  * DataStructures::VertexEliminationScratchData(),
1221  * DataStructures::NodeEliminationCopyData<dim>());
1222  * }
1223  *
1224  *
1225  *
1226  * @endcode
1227  *
1228  *
1229  * <a name="Velocitysolutionrecovery"></a>
1230  * <h4>Velocity solution recovery</h4>
1231  *
1232 
1233  *
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.
1241  *
1242  * @code
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> &copy_data)
1248  * {
1249  * unsigned int n_edges = velocity_indices.at((*n_it).first).size();
1250  * unsigned int n_cells = pressure_indices.at((*n_it).first).size();
1251  *
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);
1256  *
1257  * copy_data.vertex_vel_solution.reinit(n_edges);
1258  *
1259  * std::set<types::global_dof_index>::iterator p_it;
1260  * unsigned int i;
1261  *
1262  * for (p_it = pressure_indices[(*n_it).first].begin(), i = 0;
1263  * p_it != pressure_indices[(*n_it).first].end();
1264  * ++p_it, ++i)
1265  * scratch_data.local_pressure_solution(i) = pres_solution(*p_it - n_v);
1266  *
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);
1271  *
1272  * copy_data.p = (*n_it).first;
1273  * }
1274  *
1275  *
1276  * @endcode
1277  *
1278  * Copy nodal velocities to a global solution vector by using
1279  * local computations and indices from early stages.
1280  *
1281  * @code
1282  * template <int dim>
1283  * void MultipointMixedDarcyProblem<dim>::
1284  * copy_node_velocity_to_global(const DataStructures::NodeEliminationCopyData<dim> &copy_data)
1285  * {
1286  * std::set<types::global_dof_index>::iterator vi_it;
1287  * unsigned int i;
1288  *
1289  * for (vi_it = velocity_indices[copy_data.p].begin(), i = 0;
1290  * vi_it != velocity_indices[copy_data.p].end();
1291  * ++vi_it, ++i)
1292  * vel_solution(*vi_it) += copy_data.vertex_vel_solution(i);
1293  * }
1294  *
1295  *
1296  * @endcode
1297  *
1298  * Use @ref WorkStream to run everything concurrently.
1299  *
1300  * @code
1301  * template <int dim>
1302  * void MultipointMixedDarcyProblem<dim>::velocity_recovery()
1303  * {
1304  * TimerOutput::Scope t(computing_timer, "Velocity solution recovery");
1305  *
1306  * QGaussLobatto<dim> quad(degree+1);
1307  * QGauss<dim-1> face_quad(degree);
1308  *
1309  * vel_solution.reinit(n_v);
1310  *
1311  * WorkStream::run(node_matrix.begin(),
1312  * node_matrix.end(),
1313  * *this,
1314  * &MultipointMixedDarcyProblem::velocity_assembly,
1315  * &MultipointMixedDarcyProblem::copy_node_velocity_to_global,
1316  * DataStructures::VertexEliminationScratchData(),
1317  * DataStructures::NodeEliminationCopyData<dim>());
1318  *
1319  * solution.reinit(2);
1320  * solution.block(0) = vel_solution;
1321  * solution.block(1) = pres_solution;
1322  * solution.collect_sizes();
1323  * }
1324  *
1325  *
1326  *
1327  * @endcode
1328  *
1329  *
1330  * <a name="Pressuresystemsolver"></a>
1331  * <h4>Pressure system solver</h4>
1332  *
1333 
1334  *
1335  * The solver part is trivial. We use the CG solver with no
1336  * preconditioner for simplicity.
1337  *
1338  * @code
1339  * template <int dim>
1340  * void MultipointMixedDarcyProblem<dim>::solve_pressure()
1341  * {
1342  * TimerOutput::Scope t(computing_timer, "Pressure CG solve");
1343  *
1344  * pres_solution.reinit(n_p);
1345  *
1346  * SolverControl solver_control (2.0*n_p, 1e-10);
1347  * SolverCG<> solver (solver_control);
1348  *
1349  * PreconditionIdentity identity;
1350  * solver.solve(pres_system_matrix, pres_solution, pres_rhs, identity);
1351  * }
1352  *
1353  *
1354  *
1355  * @endcode
1356  *
1357  *
1358  * <a name="Postprocessing"></a>
1359  * <h3>Postprocessing</h3>
1360  *
1361 
1362  *
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>
1366  * format.
1367  *
1368 
1369  *
1370  *
1371  * <a name="Computeerrors"></a>
1372  * <h4>Compute errors</h4>
1373  *
1374 
1375  *
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.
1385  *
1386  * @code
1387  * template <int dim>
1388  * void MultipointMixedDarcyProblem<dim>::compute_errors(const unsigned cycle)
1389  * {
1390  * TimerOutput::Scope t(computing_timer, "Compute errors");
1391  *
1392  * const ComponentSelectFunction<dim> pressure_mask(dim, dim+1);
1393  * const ComponentSelectFunction<dim> velocity_mask(std::make_pair(0, dim), dim+1);
1394  *
1395  * ExactSolution<dim> exact_solution;
1396  *
1397  * Vector<double> cellwise_errors (triangulation.n_active_cells());
1398  *
1399  * QTrapez<1> q_trapez;
1400  * QIterated<dim> quadrature(q_trapez,degree+2);
1401  * QGauss<dim> quadrature_super(degree);
1402  *
1403  * VectorTools::integrate_difference (dof_handler, solution, exact_solution,
1404  * cellwise_errors, quadrature,
1405  * VectorTools::L2_norm,
1406  * &pressure_mask);
1407  * const double p_l2_error = cellwise_errors.l2_norm();
1408  *
1409  * VectorTools::integrate_difference (dof_handler, solution, exact_solution,
1410  * cellwise_errors, quadrature_super,
1411  * VectorTools::L2_norm,
1412  * &pressure_mask);
1413  * const double p_l2_mid_error = cellwise_errors.l2_norm();
1414  *
1415  * VectorTools::integrate_difference (dof_handler, solution, exact_solution,
1416  * cellwise_errors, quadrature,
1417  * VectorTools::L2_norm,
1418  * &velocity_mask);
1419  * const double u_l2_error = cellwise_errors.l2_norm();
1420  *
1421  * VectorTools::integrate_difference (dof_handler, solution, exact_solution,
1422  * cellwise_errors, quadrature,
1423  * VectorTools::Hdiv_seminorm,
1424  * &velocity_mask);
1425  * const double u_hd_error = cellwise_errors.l2_norm();
1426  *
1427  * const unsigned int n_active_cells=triangulation.n_active_cells();
1428  * const unsigned int n_dofs=dof_handler.n_dofs();
1429  *
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);
1437  * }
1438  *
1439  *
1440  *
1441  * @endcode
1442  *
1443  *
1444  * <a name="Outputresults"></a>
1445  * <h4>Output results</h4>
1446  *
1447 
1448  *
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.
1452  *
1453  * @code
1454  * template <int dim>
1455  * void MultipointMixedDarcyProblem<dim>::output_results(const unsigned int cycle, const unsigned int refine)
1456  * {
1457  * TimerOutput::Scope t(computing_timer, "Output results");
1458  *
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);
1464  *
1465  * DataOut<dim> data_out;
1466  * data_out.add_data_vector (dof_handler, solution, solution_names, interpretation);
1467  * data_out.build_patches ();
1468  *
1469  * std::ofstream output ("solution" + std::to_string(dim) + "d-" + std::to_string(cycle) + ".vtk");
1470  * data_out.write_vtk (output);
1471  *
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");
1488  *
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);
1493  *
1494  * std::ofstream error_table_file("error" + std::to_string(dim) + "d.tex");
1495  *
1496  * if (cycle == refine-1)
1497  * {
1498  * convergence_table.write_text(std::cout);
1499  * convergence_table.write_tex(error_table_file);
1500  * }
1501  * }
1502  *
1503  *
1504  *
1505  * @endcode
1506  *
1507  *
1508  * <a name="Runfunction"></a>
1509  * <h3>Run function</h3>
1510  *
1511 
1512  *
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.
1518  *
1519  * @code
1520  * template <int dim>
1521  * void MultipointMixedDarcyProblem<dim>::run(const unsigned int refine)
1522  * {
1523  * Assert(refine > 0, ExcMessage("Must at least have 1 refinement cycle!"));
1524  *
1525  * dof_handler.clear();
1526  * triangulation.clear();
1527  * convergence_table.clear();
1528  *
1529  * for (unsigned int cycle=0; cycle<refine; ++cycle)
1530  * {
1531  * if (cycle == 0)
1532  * {
1533  * @endcode
1534  *
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
1538  * case.
1539  *
1540  * @code
1542  * triangulation.refine_global(2);
1544  * }
1545  * else
1546  * triangulation.refine_global(1);
1547  *
1548  * node_assembly();
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 ();
1556  *
1557  * computing_timer.print_summary ();
1558  * computing_timer.reset ();
1559  * }
1560  * }
1561  * }
1562  *
1563  *
1564  * @endcode
1565  *
1566  *
1567  * <a name="Thecodemaincodefunction"></a>
1568  * <h3>The <code>main</code> function</h3>
1569  *
1570 
1571  *
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.
1575  *
1576  * @code
1577  * int main ()
1578  * {
1579  * try
1580  * {
1581  * using namespace dealii;
1582  * using namespace MFMFE;
1583  *
1585  *
1586  * MultipointMixedDarcyProblem<2> mfmfe_problem(2);
1587  * mfmfe_problem.run(6);
1588  * }
1589  * catch (std::exception &exc)
1590  * {
1591  * std::cerr << std::endl << std::endl
1592  * << "----------------------------------------------------"
1593  * << std::endl;
1594  * std::cerr << "Exception on processing: " << std::endl
1595  * << exc.what() << std::endl
1596  * << "Aborting!" << std::endl
1597  * << "----------------------------------------------------"
1598  * << std::endl;
1599  *
1600  * return 1;
1601  * }
1602  * catch (...)
1603  * {
1604  * std::cerr << std::endl << std::endl
1605  * << "----------------------------------------------------"
1606  * << std::endl;
1607  * std::cerr << "Unknown exception!" << std::endl
1608  * << "Aborting!" << std::endl
1609  * << "----------------------------------------------------"
1610  * << std::endl;
1611  * return 1;
1612  * }
1613  *
1614  * return 0;
1615  * }
1616  * @endcode
1617 
1618 
1619 */
FE_DGQ
Definition: fe_dgq.h:112
update_quadrature_points
@ update_quadrature_points
Transformed quadrature points.
Definition: fe_update_flags.h:122
TimerOutput::Scope
Definition: timer.h:554
GridTools::volume
double volume(const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping=(StaticMappingQ1< dim, spacedim >::mapping))
Definition: grid_tools.cc:133
dealii
Definition: namespace_dealii.h:25
SparseMatrixEZ
Definition: sparse_matrix_ez.h:107
BlockVector< double >
std::sin
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5297
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
WorkStream
Definition: work_stream.h:157
Triangulation< dim >
std::pow
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
Definition: vectorization.h:5428
FEValuesExtractors::Scalar
Definition: fe_values_extractors.h:95
SparseMatrix< double >
TimerOutput::summary
@ summary
Definition: timer.h:605
GridTools::distort_random
void distort_random(const double factor, Triangulation< dim, spacedim > &triangulation, const bool keep_boundary=true)
Definition: grid_tools.cc:1013
Physics::Elasticity::Kinematics::C
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
TimerOutput::wall_times
@ wall_times
Definition: timer.h:649
DoFTools::always
@ always
Definition: dof_tools.h:236
update_values
@ update_values
Shape function values.
Definition: fe_update_flags.h:78
second
Point< 2 > second
Definition: grid_out.cc:4353
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
internal::p4est::functions
int(&) functions(const void *v1, const void *v2)
Definition: p4est_wrappers.cc:339
DoFHandler
Definition: dof_handler.h:205
FEValues
Definition: fe.h:38
WorkStream::run
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)
Definition: work_stream.h:1185
TimerOutput
Definition: timer.h:546
LAPACKSupport::one
static const types::blas_int one
Definition: lapack_support.h:183
FEValuesExtractors::Vector
Definition: fe_values_extractors.h:150
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
Tensor< 1, dim, double >
update_gradients
@ update_gradients
Shape function gradients.
Definition: fe_update_flags.h:84
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition: lapack_support.h:60
VectorTools::mean
@ mean
Definition: vector_tools_common.h:79
DynamicSparsityPattern
Definition: dynamic_sparsity_pattern.h:323
SparsityPattern
Definition: sparsity_pattern.h:865
std::abs
inline ::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &x)
Definition: vectorization.h:5450
FE_RT_Bubbles
Definition: fe_rt_bubbles.h:90
QGaussLobatto
Definition: quadrature_lib.h:76
FiniteElement< dim >
TensorFunction
Definition: tensor_function.h:57
QGauss
Definition: quadrature_lib.h:40
Triangulation::begin_active_face
active_face_iterator begin_active_face() const
Definition: tria.cc:12242
value
static const bool value
Definition: dof_tools_constraints.cc:433
update_JxW_values
@ update_JxW_values
Transformed quadrature weights.
Definition: fe_update_flags.h:129
update_normal_vectors
@ update_normal_vectors
Normal vectors.
Definition: fe_update_flags.h:136
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
internal::assemble
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
Definition: loop.h:71
DoFTools::count_dofs_per_component
void count_dofs_per_component(const DoFHandlerType &dof_handler, std::vector< types::global_dof_index > &dofs_per_component, const bool vector_valued_once=false, const std::vector< unsigned int > &target_component={})
Definition: dof_tools.cc:1892
GridGenerator::hyper_cube
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
LAPACKSupport::zero
static const types::blas_int zero
Definition: lapack_support.h:179
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Point< dim >
FullMatrix< double >
Function
Definition: function.h:151
triangulation
const typename ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Definition: p4est_wrappers.cc:69
FEFaceValues
Definition: fe.h:40
internal::TriangulationImplementation::n_cells
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:12618
Quadrature
Definition: quadrature.h:85
Vector< double >
GridRefinement::refine
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
Definition: grid_refinement.cc:41
DoFRenumbering::component_wise
void component_wise(DoFHandlerType &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
Definition: dof_renumbering.cc:633
FESystem
Definition: fe.h:44
ConvergenceTable
Definition: convergence_table.h:63
std::cos
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5324
internal::VectorOperations::copy
void copy(const T *begin, const T *end, U *dest)
Definition: vector_operations_internal.h:67
MultithreadInfo::set_thread_limit
static void set_thread_limit(const unsigned int max_threads=numbers::invalid_unsigned_int)
Definition: multithread_info.cc:116