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