Loading [MathJax]/extensions/TeX/newcommand.js
 Reference documentation for deal.II version 9.4.1
\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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
MultipointFluxMixedFiniteElementMethods.h
Go to the documentation of this file.
1
164 *
165 * #ifndef MFMFE_DATA_H
166 * #define MFMFE_DATA_H
167 *
168 * #include <deal.II/base/function.h>
169 * #include <deal.II/base/tensor_function.h>
170 *
171 * @endcode
172 *
173 *
174 * <a name="Dataandexactsolution"></a>
175 * <h3>Data and exact solution.</h3>
176 *
177
178 *
179 * This file declares the classes for the given data, i.e.
180 * right-hand side, exact solution, permeability tensor and
181 * boundary conditions. For simplicity only 2d cases are
182 * provided, but 3d can be added straightforwardly.
183 *
184
185 *
186 *
187 * @code
188 * namespace MFMFE
189 * {
190 * using namespace dealii;
191 *
192 * template <int dim>
193 * class RightHandSide : public Function<dim>
194 * {
195 * public:
196 * RightHandSide () : Function<dim>(1) {}
197 *
198 * virtual double value (const Point<dim> &p,
199 * const unsigned int component = 0) const;
200 * };
201 *
202 * template <int dim>
203 * double RightHandSide<dim>::value (const Point<dim> &p,
204 * const unsigned int /*component*/) const
205 * {
206 * const double x = p[0];
207 * const double y = p[1];
208 *
209 * switch (dim)
210 * {
211 * case 2:
212 * return -(x*(y*y*y*y)*6.0-(y*y)*sin(x*y*2.0)*2.0+2.0)*(x*2.0+x*x+y*y+1.0)-sin(x*y)*(cos(x*y*2.0)+(x*x)*(y*y*y)*1.2E1
213 * -x*y*sin(x*y*2.0)*2.0)*2.0-(x*2.0+2.0)*(x*2.0+(x*x)*(y*y*y*y)*3.0+y*cos(x*y*2.0))+(x*x)*(sin(x*y*2.0)
214 * -x*(y*y)*6.0)*pow(x+1.0,2.0)*2.0-x*cos(x*y)*(x*2.0+(x*x)*(y*y*y*y)*3.0+y*(pow(cos(x*y),2.0)*2.0-1.0))
215 * -x*y*cos(x*y)*((x*x)*(y*y*y)*4.0+pow(cos(x*y),2.0)*2.0-1.0);
216 * default:
217 * Assert(false, ExcMessage("The RHS data for dim != 2 is not provided"));
218 * }
219 * }
220 *
221 *
222 *
223 * template <int dim>
224 * class PressureBoundaryValues : public Function<dim>
225 * {
226 * public:
227 * PressureBoundaryValues () : Function<dim>(1) {}
228 *
229 * virtual double value (const Point<dim> &p,
230 * const unsigned int component = 0) const;
231 * };
232 *
233 * template <int dim>
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;
259 *
260 * virtual void vector_gradient (const Point<dim> &p,
261 * std::vector<Tensor<1,dim,double>> &grads) const;
262 * };
263 *
264 * template <int dim>
265 * void
266 * ExactSolution<dim>::vector_value (const Point<dim> &p,
267 * Vector<double> &values) const
268 * {
269 * Assert (values.size() == dim+1,
270 * ExcDimensionMismatch (values.size(), dim+1));
271 *
272 * const double x = p[0];
273 * const double y = p[1];
274 *
275 * switch (dim)
276 * {
277 * case 2:
278 * values(0) = -(x*2.0+(x*x)*(y*y*y*y)*3.0+y*cos(x*y*2.0))*(x*2.0+x*x+y*y+1.0)-x*sin(x*y)*(cos(x*y*2.0)+(x*x)*(y*y*y)*4.0);
279 * values(1) = -sin(x*y)*(x*2.0+(x*x)*(y*y*y*y)*3.0+y*cos(x*y*2.0))-x*(cos(x*y*2.0)+(x*x)*(y*y*y)*4.0)*pow(x+1.0,2.0);
280 * values(2) = (x*x*x)*(y*y*y*y)+cos(x*y)*sin(x*y)+x*x;
281 * break;
282 * default:
283 * Assert(false, ExcMessage("The exact solution for dim != 2 is not provided"));
284 * }
285 * }
286 *
287 * template <int dim>
288 * void
289 * ExactSolution<dim>::vector_gradient (const Point<dim> &p,
290 * std::vector<Tensor<1,dim,double>> &grads) const
291 * {
292 * const double x = p[0];
293 * const double y = p[1];
294 *
295 * switch (dim)
296 * {
297 * case 2:
298 * grads[0][0] = -(x*(y*y*y*y)*6.0-(y*y)*sin(x*y*2.0)*2.0+2.0)*(x*2.0+x*x+y*y+1.0)-sin(x*y)*(cos(x*y*2.0)
299 * +(x*x)*(y*y*y)*1.2E1-x*y*sin(x*y*2.0)*2.0)-(x*2.0+2.0)*(x*2.0+(x*x)*(y*y*y*y)*3.0
300 * +y*cos(x*y*2.0))-x*y*cos(x*y)*((x*x)*(y*y*y)*4.0+pow(cos(x*y),2.0)*2.0-1.0);
301 * grads[0][1] = -(cos(x*y*2.0)+(x*x)*(y*y*y)*1.2E1-x*y*sin(x*y*2.0)*2.0)*(x*2.0+x*x+y*y+1.0)
302 * -y*(x*2.0+(x*x)*(y*y*y*y)*3.0+y*cos(x*y*2.0))*2.0-(x*x)*cos(x*y)*((x*x)*(y*y*y)*4.0
303 * +pow(cos(x*y),2.0)*2.0-1.0)+(x*x)*sin(x*y)*(sin(x*y*2.0)-x*(y*y)*6.0)*2.0;
304 * grads[1][0] = -sin(x*y)*(x*(y*y*y*y)*6.0-(y*y)*sin(x*y*2.0)*2.0+2.0)-pow(x+1.0,2.0)*(cos(x*y*2.0)
305 * +(x*x)*(y*y*y)*1.2E1-x*y*sin(x*y*2.0)*2.0)-x*(cos(x*y*2.0)+(x*x)*(y*y*y)*4.0)*(x*2.0+2.0)
306 * -y*cos(x*y)*(x*2.0+(x*x)*(y*y*y*y)*3.0+y*(pow(cos(x*y),2.0)*2.0-1.0));
307 * grads[1][1] = -sin(x*y)*(cos(x*y*2.0)+(x*x)*(y*y*y)*1.2E1-x*y*sin(x*y*2.0)*2.0)+(x*x)*(sin(x*y*2.0)
308 * -x*(y*y)*6.0)*pow(x+1.0,2.0)*2.0-x*cos(x*y)*(x*2.0+(x*x)*(y*y*y*y)*3.0
309 * +y*(pow(cos(x*y),2.0)*2.0-1.0));
310 * grads[2] = 0;
311 * break;
312 * default:
313 * Assert(false, ExcMessage("The exact solution's gradient for dim != 2 is not provided"));
314 * }
315 * }
316 *
317 *
318 *
319 * template <int dim>
320 * class KInverse : public TensorFunction<2,dim>
321 * {
322 * public:
323 * KInverse () : TensorFunction<2,dim>() {}
324 *
325 * virtual void value_list (const std::vector<Point<dim> > &points,
326 * std::vector<Tensor<2,dim> > &values) const;
327 * };
328 *
329 * template <int dim>
330 * void
331 * KInverse<dim>::value_list (const std::vector<Point<dim> > &points,
332 * std::vector<Tensor<2,dim> > &values) const
333 * {
334 * Assert (points.size() == values.size(),
335 * ExcDimensionMismatch (points.size(), values.size()));
336 *
337 * for (unsigned int p=0; p<points.size(); ++p)
338 * {
339 * values[p].clear ();
340 *
341 * const double x = points[p][0];
342 * const double y = points[p][1];
343 *
344 * switch (dim)
345 * {
346 * case 2:
347 * values[p][0][0] = pow(x+1.0,2.0)/(x*4.0+(x*x)*(y*y)-pow(sin(x*y),2.0)+x*(y*y)*2.0+(x*x)*6.0+(x*x*x)*4.0+x*x*x*x+y*y+1.0);
348 * values[p][0][1] = -sin(x*y)/(x*4.0+(x*x)*(y*y)-pow(sin(x*y),2.0)+x*(y*y)*2.0+(x*x)*6.0+(x*x*x)*4.0+x*x*x*x+y*y+1.0);
349 * values[p][1][0] = -sin(x*y)/(x*4.0+(x*x)*(y*y)-pow(sin(x*y),2.0)+x*(y*y)*2.0+(x*x)*6.0+(x*x*x)*4.0+x*x*x*x+y*y+1.0);
350 * values[p][1][1] = (x*2.0+x*x+y*y+1.0)/(x*4.0+(x*x)*(y*y)-pow(sin(x*y),2.0)+x*(y*y)*2.0+(x*x)*6.0+(x*x*x)*4.0+x*x*x*x+y*y+1.0);
351 * break;
352 * default:
353 * Assert(false, ExcMessage("The inverse of permeability tensor for dim != 2 is not provided"));
354 * }
355 * }
356 * }
357 * }
358 *
359 * #endif //MFMFE_DATA_H
360 * @endcode
361
362
363<a name="ann-mfmfe.cc"></a>
364<h1>Annotated version of mfmfe.cc</h1>
365 *
366 *
367 *
368 *
369 * @code
370 * /* ---------------------------------------------------------------------
371 * *
372 * * This file is part of the deal.II Code Gallery.
373 * *
374 * * ---------------------------------------------------------------------
375 * *
376 * * Author: Ilona Ambartsumyan, Eldar Khattatov, University of Pittsburgh, 2018
377 * */
378 *
379 *
380 * @endcode
381 *
382 *
383 * <a name="Includefiles"></a>
384 * <h3>Include files</h3>
385 *
386
387 *
388 * As usual, the list of necessary header files. There is not
389 * much new here, the files are included in order
390 * base-lac-grid-dofs-numerics followed by the C++ headers.
391 *
392 * @code
393 * #include <deal.II/base/convergence_table.h>
394 * #include <deal.II/base/quadrature_lib.h>
395 * #include <deal.II/base/logstream.h>
396 * #include <deal.II/base/timer.h>
397 * #include <deal.II/base/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);
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*/
Definition: fe_dgq.h:111
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
Definition: tensor.h:503
@ wall_times
Definition: timer.h:653
@ summary
Definition: timer.h:609
face_iterator end_face() const
unsigned int n_active_cells() const
unsigned int n_vertices() const
active_face_iterator begin_active_face() const
Definition: vector.h:109
@ 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.
Point< 2 > second
Definition: grid_out.cc:4604
#define Assert(cond, exc)
Definition: exceptions.h:1473
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:438
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={})
Definition: dof_tools.cc:1888
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, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
Definition: grid_tools.cc:139
void distort_random(const double factor, Triangulation< dim, spacedim > &triangulation, const bool keep_boundary=true, const unsigned int seed=boost::random::mt19937::default_seed)
Definition: grid_tools.cc:2259
@ 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:462
void run(const Iterator &begin, const typename identity< Iterator >::type &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)
Definition: work_stream.h:474
void run(const std::vector< std::vector< Iterator > > &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition: work_stream.h:1275
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:13734
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
::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