Reference documentation for deal.II version 9.6.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\}}\)
Loading...
Searching...
No Matches
MultipointFluxMixedFiniteElementMethods.h
Go to the documentation of this file.
1
168 *  
169 *   #ifndef MFMFE_DATA_H
170 *   #define MFMFE_DATA_H
171 *  
172 *   #include <deal.II/base/function.h>
173 *   #include <deal.II/base/tensor_function.h>
174 *  
175 * @endcode
176 *
177 *
178 * <a name="data.h-Dataandexactsolution"></a>
179 * <h3>Data and exact solution.</h3>
180 *
181
182 *
183 * This file declares the classes for the given data, i.e.
184 * right-hand side, exact solution, permeability tensor and
185 * boundary conditions. For simplicity only 2d cases are
186 * provided, but 3d can be added straightforwardly.
187 *
188
189 *
190 *
191 * @code
192 *   namespace MFMFE
193 *   {
194 *   using namespace dealii;
195 *  
196 *   template <int dim>
197 *   class RightHandSide : public Function<dim>
198 *   {
199 *   public:
200 *   RightHandSide () : Function<dim>(1) {}
201 *  
202 *   virtual double value (const Point<dim> &p,
203 *   const unsigned int component = 0) const override;
204 *   };
205 *  
206 *   template <int dim>
207 *   double RightHandSide<dim>::value (const Point<dim> &p,
208 *   const unsigned int /*component*/) const
209 *   {
210 *   const double x = p[0];
211 *   const double y = p[1];
212 *  
213 *   switch (dim)
214 *   {
215 *   case 2:
216 *   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
217 *   -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)
218 *   -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))
219 *   -x*y*cos(x*y)*((x*x)*(y*y*y)*4.0+pow(cos(x*y),2.0)*2.0-1.0);
220 *   default:
221 *   Assert(false, ExcMessage("The RHS data for dim != 2 is not provided"));
222 *   }
223 *   }
224 *  
225 *  
226 *  
227 *   template <int dim>
228 *   class PressureBoundaryValues : public Function<dim>
229 *   {
230 *   public:
231 *   PressureBoundaryValues () : Function<dim>(1) {}
232 *  
233 *   virtual double value (const Point<dim> &p,
234 *   const unsigned int component = 0) const override;
235 *   };
236 *  
237 *   template <int dim>
238 *   double PressureBoundaryValues<dim>::value (const Point<dim> &p,
239 *   const unsigned int /*component*/) const
240 *   {
241 *   const double x = p[0];
242 *   const double y = p[1];
243 *  
244 *   switch (dim)
245 *   {
246 *   case 2:
247 *   return (x*x*x)*(y*y*y*y)+cos(x*y)*sin(x*y)+x*x;
248 *   default:
249 *   Assert(false, ExcMessage("The BC data for dim != 2 is not provided"));
250 *   }
251 *   }
252 *  
253 *  
254 *  
255 *   template <int dim>
256 *   class ExactSolution : public Function<dim>
257 *   {
258 *   public:
259 *   ExactSolution () : Function<dim>(dim+1) {}
260 *  
261 *   virtual void vector_value (const Point<dim> &p,
262 *   Vector<double> &value) const override;
263 *  
264 *   virtual void vector_gradient (const Point<dim> &p,
265 *   std::vector<Tensor<1,dim,double>> &grads) const override;
266 *   };
267 *  
268 *   template <int dim>
269 *   void
270 *   ExactSolution<dim>::vector_value (const Point<dim> &p,
271 *   Vector<double> &values) const
272 *   {
273 *   Assert (values.size() == dim+1,
274 *   ExcDimensionMismatch (values.size(), dim+1));
275 *  
276 *   const double x = p[0];
277 *   const double y = p[1];
278 *  
279 *   switch (dim)
280 *   {
281 *   case 2:
282 *   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);
283 *   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);
284 *   values(2) = (x*x*x)*(y*y*y*y)+cos(x*y)*sin(x*y)+x*x;
285 *   break;
286 *   default:
287 *   Assert(false, ExcMessage("The exact solution for dim != 2 is not provided"));
288 *   }
289 *   }
290 *  
291 *   template <int dim>
292 *   void
293 *   ExactSolution<dim>::vector_gradient (const Point<dim> &p,
294 *   std::vector<Tensor<1,dim,double>> &grads) const
295 *   {
296 *   const double x = p[0];
297 *   const double y = p[1];
298 *  
299 *   switch (dim)
300 *   {
301 *   case 2:
302 *   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)
303 *   +(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
304 *   +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);
305 *   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)
306 *   -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
307 *   +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;
308 *   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)
309 *   +(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)
310 *   -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));
311 *   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)
312 *   -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
313 *   +y*(pow(cos(x*y),2.0)*2.0-1.0));
314 *   grads[2] = 0;
315 *   break;
316 *   default:
317 *   Assert(false, ExcMessage("The exact solution's gradient for dim != 2 is not provided"));
318 *   }
319 *   }
320 *  
321 *  
322 *  
323 *   template <int dim>
324 *   class KInverse : public TensorFunction<2,dim>
325 *   {
326 *   public:
327 *   KInverse () : TensorFunction<2,dim>() {}
328 *  
329 *   virtual void value_list (const std::vector<Point<dim> > &points,
330 *   std::vector<Tensor<2,dim> > &values) const override;
331 *   };
332 *  
333 *   template <int dim>
334 *   void
335 *   KInverse<dim>::value_list (const std::vector<Point<dim> > &points,
336 *   std::vector<Tensor<2,dim> > &values) const
337 *   {
338 *   Assert (points.size() == values.size(),
339 *   ExcDimensionMismatch (points.size(), values.size()));
340 *  
341 *   for (unsigned int p=0; p<points.size(); ++p)
342 *   {
343 *   values[p].clear ();
344 *  
345 *   const double x = points[p][0];
346 *   const double y = points[p][1];
347 *  
348 *   switch (dim)
349 *   {
350 *   case 2:
351 *   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);
352 *   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);
353 *   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);
354 *   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);
355 *   break;
356 *   default:
357 *   Assert(false, ExcMessage("The inverse of permeability tensor for dim != 2 is not provided"));
358 *   }
359 *   }
360 *   }
361 *   }
362 *  
363 *   #endif //MFMFE_DATA_H
364 * @endcode
365
366
367<a name="ann-mfmfe.cc"></a>
368<h1>Annotated version of mfmfe.cc</h1>
369 *
370 *
371 *
372 *
373 * @code
374 *   /* -----------------------------------------------------------------------------
375 *   *
376 *   * SPDX-License-Identifier: LGPL-2.1-or-later
377 *   * Copyright (C) 2018 Ilona Ambartsumyan
378 *   * Copyright (C) 2018 Eldar Khattatov
379 *   *
380 *   * This file is part of the deal.II code gallery.
381 *   *
382 *   * -----------------------------------------------------------------------------
383 *   *
384 *   * Author: Ilona Ambartsumyan, Eldar Khattatov, University of Pittsburgh, 2018
385 *   */
386 *  
387 *  
388 * @endcode
389 *
390 *
391 * <a name="mfmfe.cc-Includefiles"></a>
392 * <h3>Include files</h3>
393 *
394
395 *
396 * As usual, the list of necessary header files. There is not
397 * much new here, the files are included in order
398 * base-lac-grid-dofs-numerics followed by the C++ headers.
399 *
400 * @code
401 *   #include <deal.II/base/convergence_table.h>
402 *   #include <deal.II/base/quadrature_lib.h>
403 *   #include <deal.II/base/logstream.h>
404 *   #include <deal.II/base/timer.h>
405 *   #include <deal.II/base/utilities.h>
406 *   #include <deal.II/base/work_stream.h>
407 *  
408 *   #include <deal.II/lac/full_matrix.h>
409 *   #include <deal.II/lac/solver_cg.h>
410 *   #include <deal.II/lac/block_sparse_matrix.h>
411 *   #include <deal.II/lac/block_vector.h>
412 *   #include <deal.II/lac/precondition.h>
413 *  
414 *   #include <deal.II/grid/grid_generator.h>
415 *   #include <deal.II/grid/grid_tools.h>
416 *   #include <deal.II/grid/grid_in.h>
417 *   #include <deal.II/grid/tria.h>
418 *   #include <deal.II/dofs/dof_renumbering.h>
419 *   #include <deal.II/dofs/dof_tools.h>
420 *   #include <deal.II/fe/fe_dgq.h>
421 *   #include <deal.II/fe/fe_system.h>
422 *   #include <deal.II/fe/fe_tools.h>
423 *   #include <deal.II/numerics/vector_tools.h>
424 *   #include <deal.II/numerics/matrix_tools.h>
425 *   #include <deal.II/numerics/data_out.h>
426 *  
427 *   #include <fstream>
428 *   #include <unordered_map>
429 *  
430 * @endcode
431 *
432 * This is a header needed for the purposes of the
433 * multipoint flux mixed method, as it declares the
434 * new enhanced Raviart-Thomas finite element.
435 *
436 * @code
437 *   #include <deal.II/fe/fe_rt_bubbles.h>
438 *  
439 * @endcode
440 *
441 * For the sake of readability, the classes representing
442 * data, i.e. RHS, BCs, permeability tensor and the exact
443 * solution are placed in a file data.h which is included
444 * here
445 *
446 * @code
447 *   #include "data.h"
448 *  
449 * @endcode
450 *
451 * As always the program is in the namespace of its own with
452 * the deal.II classes and functions imported into it
453 *
454 * @code
455 *   namespace MFMFE
456 *   {
457 *   using namespace dealii;
458 *  
459 * @endcode
460 *
461 *
462 * <a name="mfmfe.cc-Definitionofmultipointfluxassemblydatastructures"></a>
463 * <h3>Definition of multipoint flux assembly data structures</h3>
464 *
465
466 *
467 * The main idea of the MFMFE method is to perform local elimination
468 * of the velocity variables in order to obtain the resulting
469 * pressure system. Since in deal.II assembly happens cell-wise,
470 * some extra work needs to be done in order to get the local
471 * mass matrices @f$A_i@f$ and the corresponding to them @f$B_i@f$.
472 *
473 * @code
474 *   namespace DataStructures
475 *   {
476 * @endcode
477 *
478 * This will be achieved by assembling cell-wise, but instead of placing
479 * the terms into a global system matrix, they will populate node-associated
480 * full matrices. For this, a data structure with fast lookup is crucial, hence
481 * the hash table, with the keys as Point<dim>
482 *
483 * @code
484 *   template <int dim>
485 *   struct hash_points
486 *   {
487 *   size_t operator()(const Point<dim> &p) const
488 *   {
489 *   size_t h1,h2,h3;
490 *   h1 = std::hash<double>()(p[0]);
491 *  
492 *   switch (dim)
493 *   {
494 *   case 1:
495 *   return h1;
496 *   case 2:
497 *   h2 = std::hash<double>()(p[1]);
498 *   return (h1 ^ h2);
499 *   case 3:
500 *   h2 = std::hash<double>()(p[1]);
501 *   h3 = std::hash<double>()(p[2]);
502 *   return (h1 ^ (h2 << 1)) ^ h3;
503 *   default:
504 *   Assert(false, ExcNotImplemented());
505 *   }
506 *   }
507 *   };
508 *  
509 * @endcode
510 *
511 * Here, the actual hash-tables are defined. We use the C++ STL <code>unordered_map</code>,
512 * with the hash function specified above. For convenience these are aliased as follows
513 *
514 * @code
515 *   template <int dim>
516 *   using PointToMatrixMap = std::unordered_map<Point<dim>, std::map<std::pair<types::global_dof_index,types::global_dof_index>, double>, hash_points<dim>>;
517 *  
518 *   template <int dim>
519 *   using PointToVectorMap = std::unordered_map<Point<dim>, std::map<types::global_dof_index, double>, hash_points<dim>>;
520 *  
521 *   template <int dim>
522 *   using PointToIndexMap = std::unordered_map<Point<dim>, std::set<types::global_dof_index>, hash_points<dim>>;
523 *  
524 * @endcode
525 *
526 * Next, since this particular program allows for the use of
527 * multiple threads, the helper CopyData structures
528 * are defined. There are two kinds of these, one is used
529 * for the copying cell-wise contributions to the corresponding
530 * node-associated data structures...
531 *
532 * @code
533 *   template <int dim>
534 *   struct NodeAssemblyCopyData
535 *   {
536 *   PointToMatrixMap<dim> cell_mat;
537 *   PointToVectorMap<dim> cell_vec;
538 *   PointToIndexMap<dim> local_pres_indices;
539 *   PointToIndexMap<dim> local_vel_indices;
540 *   std::vector<types::global_dof_index> local_dof_indices;
541 *   };
542 *  
543 * @endcode
544 *
545 * ... and the other one for the actual process of
546 * local velocity elimination and assembling the global
547 * pressure system:
548 *
549 * @code
550 *   template <int dim>
551 *   struct NodeEliminationCopyData
552 *   {
553 *   FullMatrix<double> node_pres_matrix;
554 *   Vector<double> node_pres_rhs;
555 *   FullMatrix<double> Ainverse;
556 *   FullMatrix<double> pressure_matrix;
557 *   Vector<double> velocity_rhs;
558 *   Vector<double> vertex_vel_solution;
559 *   Point<dim> p;
560 *   };
561 *  
562 * @endcode
563 *
564 * Similarly, two ScratchData classes are defined.
565 * One for the assembly part, where we need
566 * FEValues, FEFaceValues, Quadrature and storage
567 * for the basis functions...
568 *
569 * @code
570 *   template <int dim>
571 *   struct NodeAssemblyScratchData
572 *   {
573 *   NodeAssemblyScratchData (const FiniteElement<dim> &fe,
574 *   const Triangulation<dim> &tria,
575 *   const Quadrature<dim> &quad,
576 *   const Quadrature<dim-1> &f_quad);
577 *  
578 *   NodeAssemblyScratchData (const NodeAssemblyScratchData &scratch_data);
579 *  
580 *   FEValues<dim> fe_values;
581 *   FEFaceValues<dim> fe_face_values;
582 *   std::vector<unsigned int> n_faces_at_vertex;
583 *  
584 *   const unsigned long num_cells;
585 *  
586 *   std::vector<Tensor<2,dim>> k_inverse_values;
587 *   std::vector<double> rhs_values;
588 *   std::vector<double> pres_bc_values;
589 *  
590 *   std::vector<Tensor<1,dim> > phi_u;
591 *   std::vector<double> div_phi_u;
592 *   std::vector<double> phi_p;
593 *   };
594 *  
595 *   template <int dim>
596 *   NodeAssemblyScratchData<dim>::
597 *   NodeAssemblyScratchData (const FiniteElement<dim> &fe,
598 *   const Triangulation<dim> &tria,
599 *   const Quadrature<dim> &quad,
600 *   const Quadrature<dim-1> &f_quad)
601 *   :
602 *   fe_values (fe,
603 *   quad,
606 *   fe_face_values (fe,
607 *   f_quad,
610 *   num_cells(tria.n_active_cells()),
611 *   k_inverse_values(quad.size()),
612 *   rhs_values(quad.size()),
613 *   pres_bc_values(f_quad.size()),
614 *   phi_u(fe.dofs_per_cell),
615 *   div_phi_u(fe.dofs_per_cell),
616 *   phi_p(fe.dofs_per_cell)
617 *   {
618 *   n_faces_at_vertex.resize(tria.n_vertices(), 0);
619 *   typename Triangulation<dim>::active_face_iterator face = tria.begin_active_face(), endf = tria.end_face();
620 *  
621 *   for (; face != endf; ++face)
622 *   for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
623 *   n_faces_at_vertex[face->vertex_index(v)] += 1;
624 *   }
625 *  
626 *   template <int dim>
627 *   NodeAssemblyScratchData<dim>::
628 *   NodeAssemblyScratchData (const NodeAssemblyScratchData &scratch_data)
629 *   :
630 *   fe_values (scratch_data.fe_values.get_fe(),
631 *   scratch_data.fe_values.get_quadrature(),
634 *   fe_face_values (scratch_data.fe_face_values.get_fe(),
635 *   scratch_data.fe_face_values.get_quadrature(),
638 *   n_faces_at_vertex(scratch_data.n_faces_at_vertex),
639 *   num_cells(scratch_data.num_cells),
640 *   k_inverse_values(scratch_data.k_inverse_values),
641 *   rhs_values(scratch_data.rhs_values),
642 *   pres_bc_values(scratch_data.pres_bc_values),
643 *   phi_u(scratch_data.phi_u),
644 *   div_phi_u(scratch_data.div_phi_u),
645 *   phi_p(scratch_data.phi_p)
646 *   {}
647 *  
648 * @endcode
649 *
650 * ...and the other, simpler one, for the velocity elimination and recovery
651 *
652 * @code
653 *   struct VertexEliminationScratchData
654 *   {
655 *   VertexEliminationScratchData () = default;
656 *   VertexEliminationScratchData (const VertexEliminationScratchData &scratch_data);
657 *  
658 *   FullMatrix<double> velocity_matrix;
659 *   Vector<double> pressure_rhs;
660 *  
661 *   Vector<double> local_pressure_solution;
662 *   Vector<double> tmp_rhs1;
663 *   Vector<double> tmp_rhs2;
664 *   Vector<double> tmp_rhs3;
665 *   };
666 *  
667 *   VertexEliminationScratchData::
668 *   VertexEliminationScratchData (const VertexEliminationScratchData &scratch_data)
669 *   :
670 *   velocity_matrix(scratch_data.velocity_matrix),
671 *   pressure_rhs(scratch_data.pressure_rhs),
672 *   local_pressure_solution(scratch_data.local_pressure_solution),
673 *   tmp_rhs1(scratch_data.tmp_rhs1),
674 *   tmp_rhs2(scratch_data.tmp_rhs2),
675 *   tmp_rhs3(scratch_data.tmp_rhs3)
676 *   {}
677 *   }
678 *  
679 *  
680 *  
681 * @endcode
682 *
683 *
684 * <a name="mfmfe.cc-ThecodeMultipointMixedDarcyProblemcodeclasstemplate"></a>
685 * <h3>The <code>MultipointMixedDarcyProblem</code> class template</h3>
686 *
687
688 *
689 * The main class, besides the constructor and destructor, has only one public member
690 * <code>run()</code>, similarly to the tutorial programs. The private members can
691 * be grouped into the ones that are used for the cell-wise assembly, vertex elimination,
692 * pressure solve, vertex velocity recovery and postprocessing. Apart from the
693 * MFMFE-specific data structures, the rest of the members should look familiar.
694 *
695 * @code
696 *   template <int dim>
697 *   class MultipointMixedDarcyProblem
698 *   {
699 *   public:
700 *   MultipointMixedDarcyProblem (const unsigned int degree);
701 *   ~MultipointMixedDarcyProblem ();
702 *   void run (const unsigned int refine);
703 *   private:
704 *   void assemble_system_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
705 *   DataStructures::NodeAssemblyScratchData<dim> &scratch_data,
706 *   DataStructures::NodeAssemblyCopyData<dim> &copy_data);
707 *   void copy_cell_to_node(const DataStructures::NodeAssemblyCopyData<dim> &copy_data);
708 *   void node_assembly();
709 *   void make_cell_centered_sp ();
710 *   void nodal_elimination(const typename DataStructures::PointToMatrixMap<dim>::iterator &n_it,
711 *   DataStructures::VertexEliminationScratchData &scratch_data,
712 *   DataStructures::NodeEliminationCopyData<dim> &copy_data);
713 *   void copy_node_to_system(const DataStructures::NodeEliminationCopyData<dim> &copy_data);
714 *   void pressure_assembly ();
715 *   void solve_pressure ();
716 *   void velocity_assembly (const typename DataStructures::PointToMatrixMap<dim>::iterator &n_it,
717 *   DataStructures::VertexEliminationScratchData &scratch_data,
718 *   DataStructures::NodeEliminationCopyData<dim> &copy_data);
719 *   void copy_node_velocity_to_global(const DataStructures::NodeEliminationCopyData<dim> &copy_data);
720 *   void velocity_recovery ();
721 *   void reset_data_structures ();
722 *   void compute_errors (const unsigned int cycle);
723 *   void output_results (const unsigned int cycle, const unsigned int refine);
724 *  
725 *   const unsigned int degree;
727 *   FESystem<dim> fe;
728 *   DoFHandler<dim> dof_handler;
729 *   BlockVector<double> solution;
730 *  
731 *   SparsityPattern cell_centered_sp;
732 *   SparseMatrix<double> pres_system_matrix;
733 *   Vector<double> pres_rhs;
734 *  
735 *   std::unordered_map<Point<dim>, FullMatrix<double>, DataStructures::hash_points<dim>> pressure_matrix;
736 *   std::unordered_map<Point<dim>, FullMatrix<double>, DataStructures::hash_points<dim>> A_inverse;
737 *   std::unordered_map<Point<dim>, Vector<double>, DataStructures::hash_points<dim>> velocity_rhs;
738 *  
739 *   DataStructures::PointToMatrixMap<dim> node_matrix;
740 *   DataStructures::PointToVectorMap<dim> node_rhs;
741 *  
742 *   DataStructures::PointToIndexMap<dim> pressure_indices;
743 *   DataStructures::PointToIndexMap<dim> velocity_indices;
744 *  
745 *   unsigned long n_v, n_p;
746 *  
747 *   Vector<double> pres_solution;
748 *   Vector<double> vel_solution;
749 *  
750 *   ConvergenceTable convergence_table;
751 *   TimerOutput computing_timer;
752 *   };
753 *  
754 * @endcode
755 *
756 *
757 * <a name="mfmfe.cc-Constructoranddestructorcodereset_data_structurescode"></a>
758 * <h4>Constructor and destructor, <code>reset_data_structures</code></h4>
759 *
760
761 *
762 * In the constructor of this class, we store the value that was
763 * passed in concerning the degree of the finite elements we shall use (a
764 * degree of one would mean the use of @ref FE_RT_Bubbles(1) and @ref FE_DGQ(0)),
765 * and then construct the vector valued element belonging to the space @f$V_h^k@f$ described
766 * in the introduction. The constructor also takes care of initializing the
767 * computing timer, as it is of interest for us how well our method performs.
768 *
769 * @code
770 *   template <int dim>
771 *   MultipointMixedDarcyProblem<dim>::MultipointMixedDarcyProblem (const unsigned int degree)
772 *   :
773 *   degree(degree),
774 *   fe(FE_RT_Bubbles<dim>(degree), 1,
775 *   FE_DGQ<dim>(degree-1), 1),
776 *   dof_handler(triangulation),
777 *   computing_timer(std::cout, TimerOutput::summary,
779 *   {}
780 *  
781 *  
782 * @endcode
783 *
784 * The destructor clears the <code>dof_handler</code> and
785 * all of the data structures we used for the method.
786 *
787 * @code
788 *   template <int dim>
789 *   MultipointMixedDarcyProblem<dim>::~MultipointMixedDarcyProblem()
790 *   {
791 *   reset_data_structures ();
792 *   dof_handler.clear();
793 *   }
794 *  
795 *  
796 * @endcode
797 *
798 * This method clears all the data that was used after one refinement
799 * cycle.
800 *
801 * @code
802 *   template <int dim>
803 *   void MultipointMixedDarcyProblem<dim>::reset_data_structures ()
804 *   {
805 *   pressure_indices.clear();
806 *   velocity_indices.clear();
807 *   velocity_rhs.clear();
808 *   A_inverse.clear();
809 *   pressure_matrix.clear();
810 *   node_matrix.clear();
811 *   node_rhs.clear();
812 *   }
813 *  
814 *  
815 * @endcode
816 *
817 *
818 * <a name="mfmfe.cc-Cellwiseassemblyandcreationofthelocalnodalbaseddatastructures"></a>
819 * <h4>Cell-wise assembly and creation of the local, nodal-based data structures</h4>
820 *
821
822 *
823 * First, the function that copies local cell contributions to the corresponding nodal
824 * matrices and vectors is defined. It places the values obtained from local cell integration
825 * into the correct place in a matrix/vector corresponding to a specific node.
826 *
827 * @code
828 *   template <int dim>
829 *   void MultipointMixedDarcyProblem<dim>::copy_cell_to_node(const DataStructures::NodeAssemblyCopyData<dim> &copy_data)
830 *   {
831 *   for (auto m : copy_data.cell_mat)
832 *   {
833 *   for (auto p : m.second)
834 *   node_matrix[m.first][p.first] += p.second;
835 *  
836 *   for (auto p : copy_data.cell_vec.at(m.first))
837 *   node_rhs[m.first][p.first] += p.second;
838 *  
839 *   for (auto p : copy_data.local_pres_indices.at(m.first))
840 *   pressure_indices[m.first].insert(p);
841 *  
842 *   for (auto p : copy_data.local_vel_indices.at(m.first))
843 *   velocity_indices[m.first].insert(p);
844 *   }
845 *   }
846 *  
847 *  
848 *  
849 * @endcode
850 *
851 * Second, the function that does the cell assembly is defined. While it is
852 * similar to the tutorial programs in a way it uses scrath and copy data
853 * structures, the need to localize the DOFs leads to several differences.
854 *
855 * @code
856 *   template <int dim>
857 *   void MultipointMixedDarcyProblem<dim>::
858 *   assemble_system_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
859 *   DataStructures::NodeAssemblyScratchData<dim> &scratch_data,
860 *   DataStructures::NodeAssemblyCopyData<dim> &copy_data)
861 *   {
862 *   copy_data.cell_mat.clear();
863 *   copy_data.cell_vec.clear();
864 *   copy_data.local_vel_indices.clear();
865 *   copy_data.local_pres_indices.clear();
866 *  
867 *   const unsigned int dofs_per_cell = fe.dofs_per_cell;
868 *   const unsigned int n_q_points = scratch_data.fe_values.get_quadrature().size();
869 *   const unsigned int n_face_q_points = scratch_data.fe_face_values.get_quadrature().size();
870 *  
871 *   copy_data.local_dof_indices.resize(dofs_per_cell);
872 *   cell->get_dof_indices (copy_data.local_dof_indices);
873 *  
874 *   scratch_data.fe_values.reinit (cell);
875 *  
876 *   const KInverse<dim> k_inverse;
877 *   const RightHandSide<dim> rhs;
878 *   const PressureBoundaryValues<dim> pressure_bc;
879 *  
880 *   k_inverse.value_list (scratch_data.fe_values.get_quadrature_points(), scratch_data.k_inverse_values);
881 *   rhs.value_list(scratch_data.fe_values.get_quadrature_points(), scratch_data.rhs_values);
882 *  
883 *   const FEValuesExtractors::Vector velocity (0);
884 *   const FEValuesExtractors::Scalar pressure (dim);
885 *  
886 *   const unsigned int n_vel = dim*Utilities::pow(degree+1,dim);
887 *   std::unordered_map<unsigned int, std::unordered_map<unsigned int, double>> div_map;
888 *  
889 * @endcode
890 *
891 * One, we need to be able to assemble the communication between velocity and
892 * pressure variables and put it on the right place in our final, local version
893 * of the B matrix. This is a little messy, as such communication is not in fact
894 * local, so we do it in two steps. First, we compute all relevant LHS and RHS
895 *
896 * @code
897 *   for (unsigned int q=0; q<n_q_points; ++q)
898 *   {
899 *   const Point<dim> p = scratch_data.fe_values.quadrature_point(q);
900 *  
901 *   for (unsigned int k=0; k<dofs_per_cell; ++k)
902 *   {
903 *   scratch_data.phi_u[k] = scratch_data.fe_values[velocity].value(k, q);
904 *   scratch_data.div_phi_u[k] = scratch_data.fe_values[velocity].divergence (k, q);
905 *   scratch_data.phi_p[k] = scratch_data.fe_values[pressure].value (k, q);
906 *   }
907 *  
908 *   for (unsigned int i=0; i<dofs_per_cell; ++i)
909 *   {
910 *   for (unsigned int j=n_vel; j<dofs_per_cell; ++j)
911 *   {
912 *   double div_term = (- scratch_data.div_phi_u[i] * scratch_data.phi_p[j]
913 *   - scratch_data.phi_p[i] * scratch_data.div_phi_u[j]) * scratch_data.fe_values.JxW(q);
914 *  
915 *   if (std::abs(div_term) > 1.e-12)
916 *   div_map[i][j] += div_term;
917 *   }
918 *  
919 *   double source_term = -scratch_data.phi_p[i] * scratch_data.rhs_values[q] * scratch_data.fe_values.JxW(q);
920 *  
921 *   if (std::abs(scratch_data.phi_p[i]) > 1.e-12 || std::abs(source_term) > 1.e-12)
922 *   copy_data.cell_vec[p][copy_data.local_dof_indices[i]] += source_term;
923 *   }
924 *   }
925 *  
926 * @endcode
927 *
928 * Then, by making another pass, we compute the mass matrix terms and incorporate the
929 * divergence form and RHS accordingly. This second pass, allows us to know where
930 * the total contribution will be put in the nodal data structures, as with this
931 * choice of quadrature rule and finite element only the basis functions corresponding
932 * to the same quadrature points yield non-zero contribution.
933 *
934 * @code
935 *   for (unsigned int q=0; q<n_q_points; ++q)
936 *   {
937 *   std::set<types::global_dof_index> vel_indices;
938 *   const Point<dim> p = scratch_data.fe_values.quadrature_point(q);
939 *  
940 *   for (unsigned int k=0; k<dofs_per_cell; ++k)
941 *   {
942 *   scratch_data.phi_u[k] = scratch_data.fe_values[velocity].value(k, q);
943 *   scratch_data.div_phi_u[k] = scratch_data.fe_values[velocity].divergence (k, q);
944 *   scratch_data.phi_p[k] = scratch_data.fe_values[pressure].value (k, q);
945 *   }
946 *  
947 *   for (unsigned int i=0; i<dofs_per_cell; ++i)
948 *   for (unsigned int j=i; j<dofs_per_cell; ++j)
949 *   {
950 *   double mass_term = scratch_data.phi_u[i]
951 *   * scratch_data.k_inverse_values[q]
952 *   * scratch_data.phi_u[j]
953 *   * scratch_data.fe_values.JxW(q);
954 *  
955 *   if (std::abs(mass_term) > 1.e-12)
956 *   {
957 *   copy_data.cell_mat[p][std::make_pair(copy_data.local_dof_indices[i], copy_data.local_dof_indices[j])] +=
958 *   mass_term;
959 *   vel_indices.insert(i);
960 *   copy_data.local_vel_indices[p].insert(copy_data.local_dof_indices[j]);
961 *   }
962 *   }
963 *  
964 *   for (auto i : vel_indices)
965 *   for (auto el : div_map[i])
966 *   if (std::abs(el.second) > 1.e-12)
967 *   {
968 *   copy_data.cell_mat[p][std::make_pair(copy_data.local_dof_indices[i],
969 *   copy_data.local_dof_indices[el.first])] += el.second;
970 *   copy_data.local_pres_indices[p].insert(copy_data.local_dof_indices[el.first]);
971 *   }
972 *   }
973 *  
974 * @endcode
975 *
976 * The pressure boundary conditions are computed as in @ref step_20 "step-20",
977 *
978 * @code
979 *   std::map<types::global_dof_index,double> pres_bc;
980 *   for (unsigned int face_no=0;
981 *   face_no<GeometryInfo<dim>::faces_per_cell;
982 *   ++face_no)
983 *   if (cell->at_boundary(face_no))
984 *   {
985 *   scratch_data.fe_face_values.reinit (cell, face_no);
986 *   pressure_bc.value_list(scratch_data.fe_face_values.get_quadrature_points(), scratch_data.pres_bc_values);
987 *  
988 *   for (unsigned int q=0; q<n_face_q_points; ++q)
989 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
990 *   {
991 *   double tmp = -(scratch_data.fe_face_values[velocity].value(i, q) *
992 *   scratch_data.fe_face_values.normal_vector(q) *
993 *   scratch_data.pres_bc_values[q] *
994 *   scratch_data.fe_face_values.JxW(q));
995 *  
996 *   if (std::abs(tmp) > 1.e-12)
997 *   pres_bc[copy_data.local_dof_indices[i]] += tmp;
998 *   }
999 *   }
1000 *  
1001 * @endcode
1002 *
1003 * ...but we distribute them to the corresponding nodal data structures
1004 *
1005 * @code
1006 *   for (auto m : copy_data.cell_vec)
1007 *   for (unsigned int i=0; i<dofs_per_cell; ++i)
1008 *   if (std::abs(pres_bc[copy_data.local_dof_indices[i]]) > 1.e-12)
1009 *   copy_data.cell_vec[m.first][copy_data.local_dof_indices[i]] += pres_bc[copy_data.local_dof_indices[i]];
1010 *   }
1011 *  
1012 *  
1013 * @endcode
1014 *
1015 * Finally, <code>node_assembly()</code> takes care of all the
1016 * local computations via WorkStream mechanism. Notice that the choice
1017 * of the quadrature rule here is dictated by the formulation of the
1018 * method. It has to be <code>degree+1</code> points Gauss-Lobatto
1019 * for the volume integrals and <code>degree</code> for the face ones,
1020 * as mentioned in the introduction.
1021 *
1022 * @code
1023 *   template <int dim>
1024 *   void MultipointMixedDarcyProblem<dim>::node_assembly()
1025 *   {
1026 *   TimerOutput::Scope t(computing_timer, "Nodal assembly");
1027 *  
1028 *   dof_handler.distribute_dofs(fe);
1029 *   DoFRenumbering::component_wise (dof_handler);
1030 *   const std::vector<types::global_dof_index> dofs_per_component
1031 *   = DoFTools::count_dofs_per_fe_component (dof_handler);
1032 *  
1033 *   QGaussLobatto<dim> quad(degree+1);
1034 *   QGauss<dim-1> face_quad(degree);
1035 *  
1036 *   n_v = dofs_per_component[0];
1037 *   n_p = dofs_per_component[dim];
1038 *  
1039 *   pres_rhs.reinit(n_p);
1040 *  
1041 *   WorkStream::run(dof_handler.begin_active(),
1042 *   dof_handler.end(),
1043 *   *this,
1044 *   &MultipointMixedDarcyProblem::assemble_system_cell,
1045 *   &MultipointMixedDarcyProblem::copy_cell_to_node,
1046 *   DataStructures::NodeAssemblyScratchData<dim>(fe, triangulation,quad,face_quad),
1047 *   DataStructures::NodeAssemblyCopyData<dim>());
1048 *   }
1049 *  
1050 * @endcode
1051 *
1052 *
1053 * <a name="mfmfe.cc-Makingthesparsitypattern"></a>
1054 * <h4>Making the sparsity pattern</h4>
1055 *
1056
1057 *
1058 * Having computed all the local contributions, we actually have
1059 * all the information needed to make a cell-centered sparsity
1060 * pattern manually. We do this here, because @ref SparseMatrixEZ
1061 * leads to a slower solution.
1062 *
1063 * @code
1064 *   template <int dim>
1065 *   void MultipointMixedDarcyProblem<dim>::make_cell_centered_sp()
1066 *   {
1067 *   TimerOutput::Scope t(computing_timer, "Make sparsity pattern");
1068 *   DynamicSparsityPattern dsp(n_p, n_p);
1069 *  
1070 *   std::set<types::global_dof_index>::iterator pi_it, pj_it;
1071 *   unsigned int i, j;
1072 *   for (auto el : node_matrix)
1073 *   for (pi_it = pressure_indices[el.first].begin(), i = 0;
1074 *   pi_it != pressure_indices[el.first].end();
1075 *   ++pi_it, ++i)
1076 *   for (pj_it = pi_it, j = 0;
1077 *   pj_it != pressure_indices[el.first].end();
1078 *   ++pj_it, ++j)
1079 *   dsp.add(*pi_it - n_v, *pj_it - n_v);
1080 *  
1081 *  
1082 *   dsp.symmetrize();
1083 *   cell_centered_sp.copy_from(dsp);
1084 *   pres_system_matrix.reinit (cell_centered_sp);
1085 *   }
1086 *  
1087 *  
1088 * @endcode
1089 *
1090 *
1091 * <a name="mfmfe.cc-Thelocaleliminationprocedure"></a>
1092 * <h4>The local elimination procedure</h4>
1093 *
1094
1095 *
1096 * This function finally performs the local elimination procedure.
1097 * Mathematically, it follows the same idea as in computing the
1098 * Schur complement (as mentioned in the introduction) but we do
1099 * so locally. Namely, local velocity DOFs are expressed in terms
1100 * of corresponding pressure values, and then used for the local
1101 * pressure systems.
1102 *
1103 * @code
1104 *   template <int dim>
1105 *   void MultipointMixedDarcyProblem<dim>::
1106 *   nodal_elimination(const typename DataStructures::PointToMatrixMap<dim>::iterator &n_it,
1107 *   DataStructures::VertexEliminationScratchData &scratch_data,
1108 *   DataStructures::NodeEliminationCopyData<dim> &copy_data)
1109 *   {
1110 *   unsigned int n_edges = velocity_indices.at((*n_it).first).size();
1111 *   unsigned int n_cells = pressure_indices.at((*n_it).first).size();
1112 *  
1113 *   scratch_data.velocity_matrix.reinit(n_edges,n_edges);
1114 *   copy_data.pressure_matrix.reinit(n_edges,n_cells);
1115 *  
1116 *   copy_data.velocity_rhs.reinit(n_edges);
1117 *   scratch_data.pressure_rhs.reinit(n_cells);
1118 *  
1119 *   {
1120 *   std::set<types::global_dof_index>::iterator vi_it, vj_it, p_it;
1121 *   unsigned int i;
1122 *   for (vi_it = velocity_indices.at((*n_it).first).begin(), i = 0;
1123 *   vi_it != velocity_indices.at((*n_it).first).end();
1124 *   ++vi_it, ++i)
1125 *   {
1126 *   unsigned int j;
1127 *   for (vj_it = velocity_indices.at((*n_it).first).begin(), j = 0;
1128 *   vj_it != velocity_indices.at((*n_it).first).end();
1129 *   ++vj_it, ++j)
1130 *   {
1131 *   scratch_data.velocity_matrix.add(i, j, node_matrix[(*n_it).first][std::make_pair(*vi_it, *vj_it)]);
1132 *   if (j != i)
1133 *   scratch_data.velocity_matrix.add(j, i, node_matrix[(*n_it).first][std::make_pair(*vi_it, *vj_it)]);
1134 *   }
1135 *  
1136 *   for (p_it = pressure_indices.at((*n_it).first).begin(), j = 0;
1137 *   p_it != pressure_indices.at((*n_it).first).end();
1138 *   ++p_it, ++j)
1139 *   copy_data.pressure_matrix.add(i, j, node_matrix[(*n_it).first][std::make_pair(*vi_it, *p_it)]);
1140 *  
1141 *   copy_data.velocity_rhs(i) += node_rhs.at((*n_it).first)[*vi_it];
1142 *   }
1143 *  
1144 *   for (p_it = pressure_indices.at((*n_it).first).begin(), i = 0;
1145 *   p_it != pressure_indices.at((*n_it).first).end();
1146 *   ++p_it, ++i)
1147 *   scratch_data.pressure_rhs(i) += node_rhs.at((*n_it).first)[*p_it];
1148 *   }
1149 *  
1150 *   copy_data.Ainverse.reinit(n_edges,n_edges);
1151 *  
1152 *   scratch_data.tmp_rhs1.reinit(n_edges);
1153 *   scratch_data.tmp_rhs2.reinit(n_edges);
1154 *   scratch_data.tmp_rhs3.reinit(n_cells);
1155 *  
1156 *   copy_data.Ainverse.invert(scratch_data.velocity_matrix);
1157 *   copy_data.node_pres_matrix.reinit(n_cells, n_cells);
1158 *   copy_data.node_pres_rhs = scratch_data.pressure_rhs;
1159 *  
1160 *   copy_data.node_pres_matrix = 0;
1161 *   copy_data.node_pres_matrix.triple_product(copy_data.Ainverse,
1162 *   copy_data.pressure_matrix,
1163 *   copy_data.pressure_matrix, true, false);
1164 *  
1165 *   copy_data.Ainverse.vmult(scratch_data.tmp_rhs1, copy_data.velocity_rhs, false);
1166 *   copy_data.pressure_matrix.Tvmult(scratch_data.tmp_rhs3, scratch_data.tmp_rhs1, false);
1167 *   copy_data.node_pres_rhs *= -1.0;
1168 *   copy_data.node_pres_rhs += scratch_data.tmp_rhs3;
1169 *  
1170 *   copy_data.p = (*n_it).first;
1171 *   }
1172 *  
1173 *  
1174 * @endcode
1175 *
1176 * Each node's pressure system is then distributed to a global pressure
1177 * system, using the indices we computed in the previous stages.
1178 *
1179 * @code
1180 *   template <int dim>
1181 *   void MultipointMixedDarcyProblem<dim>::
1182 *   copy_node_to_system(const DataStructures::NodeEliminationCopyData<dim> &copy_data)
1183 *   {
1184 *   A_inverse[copy_data.p] = copy_data.Ainverse;
1185 *   pressure_matrix[copy_data.p] = copy_data.pressure_matrix;
1186 *   velocity_rhs[copy_data.p] = copy_data.velocity_rhs;
1187 *  
1188 *   {
1189 *   std::set<types::global_dof_index>::iterator pi_it, pj_it;
1190 *   unsigned int i;
1191 *   for (pi_it = pressure_indices[copy_data.p].begin(), i = 0;
1192 *   pi_it != pressure_indices[copy_data.p].end();
1193 *   ++pi_it, ++i)
1194 *   {
1195 *   unsigned int j;
1196 *   for (pj_it = pressure_indices[copy_data.p].begin(), j = 0;
1197 *   pj_it != pressure_indices[copy_data.p].end();
1198 *   ++pj_it, ++j)
1199 *   pres_system_matrix.add(*pi_it - n_v, *pj_it - n_v, copy_data.node_pres_matrix(i, j));
1200 *  
1201 *   pres_rhs(*pi_it - n_v) += copy_data.node_pres_rhs(i);
1202 *   }
1203 *   }
1204 *   }
1205 *  
1206 *  
1207 * @endcode
1208 *
1209 * The @ref WorkStream mechanism is again used for the assembly
1210 * of the global system for the pressure variable, where the
1211 * previous functions are used to perform local computations.
1212 *
1213 * @code
1214 *   template <int dim>
1215 *   void MultipointMixedDarcyProblem<dim>::pressure_assembly()
1216 *   {
1217 *   TimerOutput::Scope t(computing_timer, "Pressure matrix assembly");
1218 *  
1219 *   QGaussLobatto<dim> quad(degree+1);
1220 *   QGauss<dim-1> face_quad(degree);
1221 *  
1222 *   pres_rhs.reinit(n_p);
1223 *  
1224 *   WorkStream::run(node_matrix.begin(),
1225 *   node_matrix.end(),
1226 *   *this,
1227 *   &MultipointMixedDarcyProblem::nodal_elimination,
1228 *   &MultipointMixedDarcyProblem::copy_node_to_system,
1229 *   DataStructures::VertexEliminationScratchData(),
1230 *   DataStructures::NodeEliminationCopyData<dim>());
1231 *   }
1232 *  
1233 *  
1234 *  
1235 * @endcode
1236 *
1237 *
1238 * <a name="mfmfe.cc-Velocitysolutionrecovery"></a>
1239 * <h4>Velocity solution recovery</h4>
1240 *
1241
1242 *
1243 * After solving for the pressure variable, we want to follow
1244 * the above procedure backwards, in order to obtain the
1245 * velocity solution (again, this is similar in nature to the
1246 * Schur complement approach, see @ref step_20 "step-20", but here it is done
1247 * locally at each node). We have almost everything computed and
1248 * stored already, including inverses of local mass matrices,
1249 * so the following is a relatively straightforward implementation.
1250 *
1251 * @code
1252 *   template <int dim>
1253 *   void MultipointMixedDarcyProblem<dim>::
1254 *   velocity_assembly (const typename DataStructures::PointToMatrixMap<dim>::iterator &n_it,
1255 *   DataStructures::VertexEliminationScratchData &scratch_data,
1256 *   DataStructures::NodeEliminationCopyData<dim> &copy_data)
1257 *   {
1258 *   unsigned int n_edges = velocity_indices.at((*n_it).first).size();
1259 *   unsigned int n_cells = pressure_indices.at((*n_it).first).size();
1260 *  
1261 *   scratch_data.tmp_rhs1.reinit(n_edges);
1262 *   scratch_data.tmp_rhs2.reinit(n_edges);
1263 *   scratch_data.tmp_rhs3.reinit(n_cells);
1264 *   scratch_data.local_pressure_solution.reinit(n_cells);
1265 *  
1266 *   copy_data.vertex_vel_solution.reinit(n_edges);
1267 *  
1268 *   std::set<types::global_dof_index>::iterator p_it;
1269 *   unsigned int i;
1270 *  
1271 *   for (p_it = pressure_indices[(*n_it).first].begin(), i = 0;
1272 *   p_it != pressure_indices[(*n_it).first].end();
1273 *   ++p_it, ++i)
1274 *   scratch_data.local_pressure_solution(i) = pres_solution(*p_it - n_v);
1275 *  
1276 *   pressure_matrix[(*n_it).first].vmult(scratch_data.tmp_rhs2, scratch_data.local_pressure_solution, false);
1277 *   scratch_data.tmp_rhs2 *= -1.0;
1278 *   scratch_data.tmp_rhs2+=velocity_rhs[(*n_it).first];
1279 *   A_inverse[(*n_it).first].vmult(copy_data.vertex_vel_solution, scratch_data.tmp_rhs2, false);
1280 *  
1281 *   copy_data.p = (*n_it).first;
1282 *   }
1283 *  
1284 *  
1285 * @endcode
1286 *
1287 * Copy nodal velocities to a global solution vector by using
1288 * local computations and indices from early stages.
1289 *
1290 * @code
1291 *   template <int dim>
1292 *   void MultipointMixedDarcyProblem<dim>::
1293 *   copy_node_velocity_to_global(const DataStructures::NodeEliminationCopyData<dim> &copy_data)
1294 *   {
1295 *   std::set<types::global_dof_index>::iterator vi_it;
1296 *   unsigned int i;
1297 *  
1298 *   for (vi_it = velocity_indices[copy_data.p].begin(), i = 0;
1299 *   vi_it != velocity_indices[copy_data.p].end();
1300 *   ++vi_it, ++i)
1301 *   vel_solution(*vi_it) += copy_data.vertex_vel_solution(i);
1302 *   }
1303 *  
1304 *  
1305 * @endcode
1306 *
1307 * Use @ref WorkStream to run everything concurrently.
1308 *
1309 * @code
1310 *   template <int dim>
1311 *   void MultipointMixedDarcyProblem<dim>::velocity_recovery()
1312 *   {
1313 *   TimerOutput::Scope t(computing_timer, "Velocity solution recovery");
1314 *  
1315 *   QGaussLobatto<dim> quad(degree+1);
1316 *   QGauss<dim-1> face_quad(degree);
1317 *  
1318 *   vel_solution.reinit(n_v);
1319 *  
1320 *   WorkStream::run(node_matrix.begin(),
1321 *   node_matrix.end(),
1322 *   *this,
1323 *   &MultipointMixedDarcyProblem::velocity_assembly,
1324 *   &MultipointMixedDarcyProblem::copy_node_velocity_to_global,
1325 *   DataStructures::VertexEliminationScratchData(),
1326 *   DataStructures::NodeEliminationCopyData<dim>());
1327 *  
1328 *   solution.reinit(2);
1329 *   solution.block(0) = vel_solution;
1330 *   solution.block(1) = pres_solution;
1331 *   solution.collect_sizes();
1332 *   }
1333 *  
1334 *  
1335 *  
1336 * @endcode
1337 *
1338 *
1339 * <a name="mfmfe.cc-Pressuresystemsolver"></a>
1340 * <h4>Pressure system solver</h4>
1341 *
1342
1343 *
1344 * The solver part is trivial. We use the CG solver with no
1345 * preconditioner for simplicity.
1346 *
1347 * @code
1348 *   template <int dim>
1349 *   void MultipointMixedDarcyProblem<dim>::solve_pressure()
1350 *   {
1351 *   TimerOutput::Scope t(computing_timer, "Pressure CG solve");
1352 *  
1353 *   pres_solution.reinit(n_p);
1354 *  
1355 *   SolverControl solver_control (static_cast<int>(2.0*n_p), 1e-10);
1356 *   SolverCG<> solver (solver_control);
1357 *  
1358 *   PreconditionIdentity identity;
1359 *   solver.solve(pres_system_matrix, pres_solution, pres_rhs, identity);
1360 *   }
1361 *  
1362 *  
1363 *  
1364 * @endcode
1365 *
1366 *
1367 * <a name="mfmfe.cc-Postprocessing"></a>
1368 * <h3>Postprocessing</h3>
1369 *
1370
1371 *
1372 * We have two postprocessing steps here, first one computes the
1373 * errors in order to populate the convergence tables. The other
1374 * one takes care of the output of the solutions in <code>.vtk</code>
1375 * format.
1376 *
1377
1378 *
1379 *
1380 * <a name="mfmfe.cc-Computeerrors"></a>
1381 * <h4>Compute errors</h4>
1382 *
1383
1384 *
1385 * The implementation of this function is almost identical to @ref step_20 "step-20".
1386 * We use @ref ComponentSelectFunction as masks to use the right
1387 * solution component (velocity or pressure) and @ref integrate_difference
1388 * to compute the errors. Since we also want to compute Hdiv seminorm of the
1389 * velocity error, one must provide gradients in the <code>ExactSolution</code>
1390 * class implementation to avoid exceptions. The only noteworthy thing here
1391 * is that we again use lower order quadrature rule instead of projecting the
1392 * solution to an appropriate space in order to show superconvergence, which is
1393 * mathematically justified.
1394 *
1395 * @code
1396 *   template <int dim>
1397 *   void MultipointMixedDarcyProblem<dim>::compute_errors(const unsigned cycle)
1398 *   {
1399 *   TimerOutput::Scope t(computing_timer, "Compute errors");
1400 *  
1401 *   const ComponentSelectFunction<dim> pressure_mask(dim, dim+1);
1402 *   const ComponentSelectFunction<dim> velocity_mask(std::make_pair(0, dim), dim+1);
1403 *  
1404 *   ExactSolution<dim> exact_solution;
1405 *  
1406 *   Vector<double> cellwise_errors (triangulation.n_active_cells());
1407 *  
1408 *   QTrapezoid<1> q_trapez;
1409 *   QIterated<dim> quadrature(q_trapez,degree+2);
1410 *   QGauss<dim> quadrature_super(degree);
1411 *  
1412 *   VectorTools::integrate_difference (dof_handler, solution, exact_solution,
1413 *   cellwise_errors, quadrature,
1414 *   VectorTools::L2_norm,
1415 *   &pressure_mask);
1416 *   const double p_l2_error = cellwise_errors.l2_norm();
1417 *  
1418 *   VectorTools::integrate_difference (dof_handler, solution, exact_solution,
1419 *   cellwise_errors, quadrature_super,
1420 *   VectorTools::L2_norm,
1421 *   &pressure_mask);
1422 *   const double p_l2_mid_error = cellwise_errors.l2_norm();
1423 *  
1424 *   VectorTools::integrate_difference (dof_handler, solution, exact_solution,
1425 *   cellwise_errors, quadrature,
1426 *   VectorTools::L2_norm,
1427 *   &velocity_mask);
1428 *   const double u_l2_error = cellwise_errors.l2_norm();
1429 *  
1430 *   VectorTools::integrate_difference (dof_handler, solution, exact_solution,
1431 *   cellwise_errors, quadrature,
1432 *   VectorTools::Hdiv_seminorm,
1433 *   &velocity_mask);
1434 *   const double u_hd_error = cellwise_errors.l2_norm();
1435 *  
1436 *   const unsigned int n_active_cells=triangulation.n_active_cells();
1437 *   const unsigned int n_dofs=dof_handler.n_dofs();
1438 *  
1439 *   convergence_table.add_value("cycle", cycle);
1440 *   convergence_table.add_value("cells", n_active_cells);
1441 *   convergence_table.add_value("dofs", n_dofs);
1442 *   convergence_table.add_value("Velocity,L2", u_l2_error);
1443 *   convergence_table.add_value("Velocity,Hdiv", u_hd_error);
1444 *   convergence_table.add_value("Pressure,L2", p_l2_error);
1445 *   convergence_table.add_value("Pressure,L2-nodal", p_l2_mid_error);
1446 *   }
1447 *  
1448 *  
1449 *  
1450 * @endcode
1451 *
1452 *
1453 * <a name="mfmfe.cc-Outputresults"></a>
1454 * <h4>Output results</h4>
1455 *
1456
1457 *
1458 * This function also follows the same idea as in @ref step_20 "step-20" tutorial
1459 * program. The only modification to it is the part involving
1460 * a convergence table.
1461 *
1462 * @code
1463 *   template <int dim>
1464 *   void MultipointMixedDarcyProblem<dim>::output_results(const unsigned int cycle, const unsigned int refine)
1465 *   {
1466 *   TimerOutput::Scope t(computing_timer, "Output results");
1467 *  
1468 *   std::vector<std::string> solution_names(dim, "u");
1469 *   solution_names.push_back ("p");
1470 *   std::vector<DataComponentInterpretation::DataComponentInterpretation>
1471 *   interpretation (dim, DataComponentInterpretation::component_is_part_of_vector);
1472 *   interpretation.push_back (DataComponentInterpretation::component_is_scalar);
1473 *  
1474 *   DataOut<dim> data_out;
1475 *   data_out.add_data_vector (dof_handler, solution, solution_names, interpretation);
1476 *   data_out.build_patches ();
1477 *  
1478 *   std::ofstream output ("solution" + std::to_string(dim) + "d-" + std::to_string(cycle) + ".vtk");
1479 *   data_out.write_vtk (output);
1480 *  
1481 *   convergence_table.set_precision("Velocity,L2", 3);
1482 *   convergence_table.set_precision("Velocity,Hdiv", 3);
1483 *   convergence_table.set_precision("Pressure,L2", 3);
1484 *   convergence_table.set_precision("Pressure,L2-nodal", 3);
1485 *   convergence_table.set_scientific("Velocity,L2", true);
1486 *   convergence_table.set_scientific("Velocity,Hdiv", true);
1487 *   convergence_table.set_scientific("Pressure,L2", true);
1488 *   convergence_table.set_scientific("Pressure,L2-nodal", true);
1489 *   convergence_table.set_tex_caption("cells", "\\# cells");
1490 *   convergence_table.set_tex_caption("dofs", "\\# dofs");
1491 *   convergence_table.set_tex_caption("Velocity,L2", " \\|\\u - \\u_h\\|_{L^2} ");
1492 *   convergence_table.set_tex_caption("Velocity,Hdiv", " \\|\\nabla\\cdot(\\u - \\u_h)\\|_{L^2} ");
1493 *   convergence_table.set_tex_caption("Pressure,L2", " \\|p - p_h\\|_{L^2} ");
1494 *   convergence_table.set_tex_caption("Pressure,L2-nodal", " \\|Qp - p_h\\|_{L^2} ");
1495 *   convergence_table.set_tex_format("cells", "r");
1496 *   convergence_table.set_tex_format("dofs", "r");
1497 *  
1498 *   convergence_table.evaluate_convergence_rates("Velocity,L2", ConvergenceTable::reduction_rate_log2);
1499 *   convergence_table.evaluate_convergence_rates("Velocity,Hdiv", ConvergenceTable::reduction_rate_log2);
1500 *   convergence_table.evaluate_convergence_rates("Pressure,L2", ConvergenceTable::reduction_rate_log2);
1501 *   convergence_table.evaluate_convergence_rates("Pressure,L2-nodal", ConvergenceTable::reduction_rate_log2);
1502 *  
1503 *   std::ofstream error_table_file("error" + std::to_string(dim) + "d.tex");
1504 *  
1505 *   if (cycle == refine-1)
1506 *   {
1507 *   convergence_table.write_text(std::cout);
1508 *   convergence_table.write_tex(error_table_file);
1509 *   }
1510 *   }
1511 *  
1512 *  
1513 *  
1514 * @endcode
1515 *
1516 *
1517 * <a name="mfmfe.cc-Runfunction"></a>
1518 * <h3>Run function</h3>
1519 *
1520
1521 *
1522 * The driver method <code>run()</code>
1523 * takes care of mesh generation and arranging calls to member methods in
1524 * the right way. It also resets data structures and clear triangulation and
1525 * DOF handler as we run the method on a sequence of refinements in order
1526 * to record convergence rates.
1527 *
1528 * @code
1529 *   template <int dim>
1530 *   void MultipointMixedDarcyProblem<dim>::run(const unsigned int refine)
1531 *   {
1532 *   Assert(refine > 0, ExcMessage("Must at least have 1 refinement cycle!"));
1533 *  
1534 *   dof_handler.clear();
1535 *   triangulation.clear();
1536 *   convergence_table.clear();
1537 *  
1538 *   for (unsigned int cycle=0; cycle<refine; ++cycle)
1539 *   {
1540 *   if (cycle == 0)
1541 *   {
1542 * @endcode
1543 *
1544 * We first generate the hyper cube and refine it twice
1545 * so that we could distort the grid slightly and
1546 * demonstrate the method's ability to work in such a
1547 * case.
1548 *
1549 * @code
1552 *   GridTools::distort_random (0.3, triangulation, true);
1553 *   }
1554 *   else
1556 *  
1557 *   node_assembly();
1558 *   make_cell_centered_sp();
1559 *   pressure_assembly();
1560 *   solve_pressure ();
1561 *   velocity_recovery ();
1562 *   compute_errors (cycle);
1563 *   output_results (cycle, refine);
1564 *   reset_data_structures ();
1565 *  
1566 *   computing_timer.print_summary ();
1567 *   computing_timer.reset ();
1568 *   }
1569 *   }
1570 *   }
1571 *  
1572 *  
1573 * @endcode
1574 *
1575 *
1576 * <a name="mfmfe.cc-Thecodemaincodefunction"></a>
1577 * <h3>The <code>main</code> function</h3>
1578 *
1579
1580 *
1581 * In the main functione we pass the order of the Finite Element as an argument
1582 * to the constructor of the Multipoint Flux Mixed Darcy problem, and the number
1583 * of refinement cycles as an argument for the run method.
1584 *
1585 * @code
1586 *   int main ()
1587 *   {
1588 *   try
1589 *   {
1590 *   using namespace dealii;
1591 *   using namespace MFMFE;
1592 *  
1594 *  
1595 *   MultipointMixedDarcyProblem<2> mfmfe_problem(2);
1596 *   mfmfe_problem.run(6);
1597 *   }
1598 *   catch (std::exception &exc)
1599 *   {
1600 *   std::cerr << std::endl << std::endl
1601 *   << "----------------------------------------------------"
1602 *   << std::endl;
1603 *   std::cerr << "Exception on processing: " << std::endl
1604 *   << exc.what() << std::endl
1605 *   << "Aborting!" << std::endl
1606 *   << "----------------------------------------------------"
1607 *   << std::endl;
1608 *  
1609 *   return 1;
1610 *   }
1611 *   catch (...)
1612 *   {
1613 *   std::cerr << std::endl << std::endl
1614 *   << "----------------------------------------------------"
1615 *   << std::endl;
1616 *   std::cerr << "Unknown exception!" << std::endl
1617 *   << "Aborting!" << std::endl
1618 *   << "----------------------------------------------------"
1619 *   << std::endl;
1620 *   return 1;
1621 *   }
1622 *  
1623 *   return 0;
1624 *   }
1625 * @endcode
1626
1627
1628*/
virtual void vector_gradient(const Point< dim > &p, std::vector< Tensor< 1, dim, RangeNumberType > > &gradients) const
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
virtual void vector_value(const Point< dim > &p, Vector< RangeNumberType > &values) const
static void set_thread_limit(const unsigned int max_threads=numbers::invalid_unsigned_int)
Definition point.h:111
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< value_type > &values) const
@ wall_times
Definition timer.h:651
void refine_global(const unsigned int times=1)
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
#define Assert(cond, exc)
typename ActiveSelector::active_cell_iterator active_cell_iterator
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
void component_wise(DoFHandler< dim, spacedim > &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
std::vector< types::global_dof_index > count_dofs_per_fe_component(const DoFHandler< dim, spacedim > &dof_handler, const bool vector_valued_once=false, const std::vector< unsigned int > &target_component={})
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
double volume(const Triangulation< dim, spacedim > &tria)
void distort_random(const double factor, Triangulation< dim, spacedim > &triangulation, const bool keep_boundary=true, const unsigned int seed=boost::random::mt19937::default_seed)
@ matrix
Contents is actually a matrix.
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * begin(VectorType &V)
constexpr T pow(const T base, const int iexp)
Definition utilities.h:966
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
void run(const std::vector< std::vector< Iterator > > &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition tria.cc:14883
void copy(const T *begin, const T *end, U *dest)
int(& functions)(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
Definition loop.h:70
STL namespace.
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation