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