PolyDEAL
 
Loading...
Searching...
No Matches
assembly_check.cc
Go to the documentation of this file.
1#include <deal.II/grid/grid_generator.h>
2#include <deal.II/grid/grid_out.h>
3
4#include <deal.II/lac/sparse_direct.h>
5#include <deal.II/lac/sparse_matrix.h>
6
7#include <deal.II/numerics/data_out.h>
8
9#include <agglomeration_handler.h>
10
11#include <algorithm>
12
13
14template <int dim>
15class RightHandSide : public Function<dim>
16{
17public:
18 RightHandSide()
19 : Function<dim>()
20 {}
21
22 virtual void
23 value_list(const std::vector<Point<dim>> &points,
24 std::vector<double> &values,
25 const unsigned int /*component*/) const override;
26};
27
28
29template <int dim>
30void
31RightHandSide<dim>::value_list(const std::vector<Point<dim>> &points,
32 std::vector<double> &values,
33 const unsigned int /*component*/) const
34{
35 for (unsigned int i = 0; i < values.size(); ++i)
36 values[i] = 8. * numbers::PI * numbers::PI *
37 std::sin(2. * numbers::PI * points[i][0]) *
38 std::sin(2. * numbers::PI * points[i][1]);
39}
40
41
42template <int dim>
43class Poisson
44{
45private:
46 void
47 make_grid();
48 void
49 setup_agglomeration();
50 void
51 distribute_jumps_and_averages(
52 FEFaceValues<dim> &fe_face,
53 FEFaceValues<dim> &fe_face1,
54 const typename DoFHandler<dim>::active_cell_iterator &cell,
55 const unsigned int f);
56 void
57 assemble_system();
58 void
59 solve();
60 void
61 output_results();
62
63
64 Triangulation<dim> tria;
65 MappingQ<dim> mapping;
66 FE_DGQ<dim> dg_fe;
67 DoFHandler<dim> classical_dh;
68 std::unique_ptr<AgglomerationHandler<dim>> ah;
69 SparsityPattern sparsity;
70 SparseMatrix<double> system_matrix;
71 Vector<double> solution;
72 Vector<double> system_rhs;
73 std::unique_ptr<GridTools::Cache<dim>> cached_tria;
74 std::unique_ptr<const Function<dim>> rhs_function;
75
76public:
77 Poisson();
78 void
79 run();
80
81 double penalty = 3.;
82};
83
84
85
86template <int dim>
87Poisson<dim>::Poisson()
88 : mapping(1)
89 , dg_fe(1)
90 , classical_dh(tria)
91{}
92
93template <int dim>
94void
95Poisson<dim>::make_grid()
96{
97 GridGenerator::hyper_cube(tria, -1, 1);
98 // tria.refine_global(2);
99 tria.refine_global(1);
100 classical_dh.distribute_dofs(dg_fe);
101 cached_tria = std::make_unique<GridTools::Cache<dim>>(tria, mapping);
102 rhs_function = std::make_unique<const RightHandSide<dim>>();
103}
104
105
106
107template <int dim>
108void
109Poisson<dim>::setup_agglomeration()
110{
111 std::vector<types::global_cell_index> idxs_to_be_agglomerated = {0, 1, 2, 3};
112
113 std::vector<typename Triangulation<dim>::active_cell_iterator>
114 cells_to_be_agglomerated;
115
116 PolyUtils::collect_cells_for_agglomeration(tria,
117 idxs_to_be_agglomerated,
118 cells_to_be_agglomerated);
119
120 std::vector<types::global_cell_index> idxs_to_be_agglomerated2 = {4, 5, 6, 7};
121
122 std::vector<typename Triangulation<dim>::active_cell_iterator>
123 cells_to_be_agglomerated2;
124 PolyUtils::collect_cells_for_agglomeration(tria,
125 idxs_to_be_agglomerated2,
126 cells_to_be_agglomerated2);
127
128 std::vector<types::global_cell_index> idxs_to_be_agglomerated3 = {8,
129 9,
130 10,
131 11};
132
133 std::vector<typename Triangulation<dim>::active_cell_iterator>
134 cells_to_be_agglomerated3;
135 PolyUtils::collect_cells_for_agglomeration(tria,
136 idxs_to_be_agglomerated3,
137 cells_to_be_agglomerated3);
138
139 std::vector<types::global_cell_index> idxs_to_be_agglomerated4 = {12,
140 13,
141 14,
142 15};
143
144 std::vector<typename Triangulation<dim>::active_cell_iterator>
145 cells_to_be_agglomerated4;
146 PolyUtils::collect_cells_for_agglomeration(tria,
147 idxs_to_be_agglomerated4,
148 cells_to_be_agglomerated4);
149 // Agglomerate the cells just stored
150 ah = std::make_unique<AgglomerationHandler<dim>>(*cached_tria);
151 ah->define_agglomerate(cells_to_be_agglomerated);
152 ah->define_agglomerate(cells_to_be_agglomerated2);
153 ah->define_agglomerate(cells_to_be_agglomerated3);
154 ah->define_agglomerate(cells_to_be_agglomerated4);
155 ah->distribute_agglomerated_dofs(dg_fe);
156 ah->create_agglomeration_sparsity_pattern(sparsity);
157 // // sanity check
158 // for (const auto &cell : ah->agglo_dh.active_cell_iterators())
159 // {
160 // const unsigned int agglo_faces_per_cell =
161 // ah->n_agglomerated_faces_per_cell(cell);
162 // std::cout << "Number of agglomerated faces for cell "
163 // << cell->active_cell_index() << " is " <<
164 // agglo_faces_per_cell
165 // << std::endl;
166
167 std::ofstream out("grid_coarsen.vtk");
168 GridOut grid_out;
169 grid_out.write_vtk(tria, out);
170 std::cout << "Grid written " << std::endl;
171}
172
173template <int dim>
174void
175Poisson<dim>::distribute_jumps_and_averages(
176 FEFaceValues<dim> &fe_face,
177 FEFaceValues<dim> &fe_face1,
178 const typename DoFHandler<dim>::active_cell_iterator &cell,
179 const unsigned int f)
180{
181 if (cell->face(f)->at_boundary())
182 {
183 std::cout << "at boundary" << std::endl;
184 fe_face.reinit(cell, f);
185 const unsigned int dofs_per_cell = fe_face.dofs_per_cell;
186 FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
187 std::vector<types::global_dof_index> local_dof_indices_bdary_cell(
188 dofs_per_cell);
189 const double hf = cell->face(f)->measure();
190
191 // Get normal vectors seen from each agglomeration.
192 const auto &normals = fe_face.get_normal_vectors();
193 cell_matrix = 0.;
194 for (unsigned int q_index : fe_face.quadrature_point_indices())
195 {
196 for (unsigned int i = 0; i < dofs_per_cell; ++i)
197 {
198 for (unsigned int j = 0; j < dofs_per_cell; ++j)
199 {
200 cell_matrix(i, j) +=
201 (-fe_face.shape_value(i, q_index) *
202 fe_face.shape_grad(j, q_index) * normals[q_index] -
203 fe_face.shape_grad(i, q_index) * normals[q_index] *
204 fe_face.shape_value(j, q_index) +
205 (penalty / hf) * fe_face.shape_value(i, q_index) *
206 fe_face.shape_value(j, q_index)) *
207 fe_face.JxW(q_index);
208 }
209 }
210 }
211
212 // distribute DoFs
213 cell->get_dof_indices(local_dof_indices_bdary_cell);
214 system_matrix.add(local_dof_indices_bdary_cell, cell_matrix);
215 }
216 else
217 {
218 fe_face.reinit(cell, f);
219 fe_face1.reinit(cell->neighbor(f), cell->neighbor_of_neighbor(f));
220
221 const unsigned int dofs_per_cell = fe_face.dofs_per_cell;
222 std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
223 std::vector<types::global_dof_index> local_dof_indices1(dofs_per_cell);
224
225 const double hf = cell->face(f)->measure();
226 FullMatrix<double> A_00(dofs_per_cell, dofs_per_cell);
227 FullMatrix<double> A_01(dofs_per_cell, dofs_per_cell);
228 FullMatrix<double> A_10(dofs_per_cell, dofs_per_cell);
229 FullMatrix<double> A_11(dofs_per_cell, dofs_per_cell);
230
231
232 // Get normal vectors seen from each agglomeration.
233 const auto &normals = fe_face1.get_normal_vectors();
234 // A00
235 for (unsigned int q_index : fe_face.quadrature_point_indices())
236 {
237 std::cout << normals[q_index] << std::endl;
238 for (unsigned int i = 0; i < dofs_per_cell; ++i)
239 {
240 for (unsigned int j = 0; j < dofs_per_cell; ++j)
241 {
242 A_00(i, j) +=
243 (-0.5 * fe_face.shape_grad(i, q_index) * normals[q_index] *
244 fe_face.shape_value(j, q_index) -
245 0.5 * fe_face.shape_grad(j, q_index) * normals[q_index] *
246 fe_face.shape_value(i, q_index) +
247 (penalty / hf) * fe_face.shape_value(i, q_index) *
248 fe_face.shape_value(j, q_index)) *
249 fe_face.JxW(q_index);
250 }
251 }
252 }
253 // A10
254 for (unsigned int q_index : fe_face.quadrature_point_indices())
255 {
256 for (unsigned int i = 0; i < dofs_per_cell; ++i)
257 {
258 for (unsigned int j = 0; j < dofs_per_cell; ++j)
259 {
260 A_01(i, j) +=
261 (+0.5 * fe_face.shape_grad(i, q_index) * normals[q_index] *
262 fe_face1.shape_value(j, q_index) -
263 0.5 * fe_face1.shape_grad(j, q_index) * normals[q_index] *
264 fe_face.shape_value(i, q_index) -
265 (penalty / hf) * fe_face.shape_value(i, q_index) *
266 fe_face1.shape_value(j, q_index)) *
267 fe_face1.JxW(q_index);
268 }
269 }
270 }
271 // A01
272 for (unsigned int q_index : fe_face.quadrature_point_indices())
273 {
274 for (unsigned int i = 0; i < dofs_per_cell; ++i)
275 {
276 for (unsigned int j = 0; j < dofs_per_cell; ++j)
277 {
278 A_10(i, j) +=
279 (+0.5 * fe_face1.shape_grad(i, q_index) * normals[q_index] *
280 fe_face.shape_value(j, q_index) -
281 0.5 * fe_face.shape_grad(j, q_index) * normals[q_index] *
282 fe_face1.shape_value(i, q_index) -
283 (penalty / hf) * fe_face1.shape_value(i, q_index) *
284 fe_face.shape_value(j, q_index)) *
285 fe_face1.JxW(q_index);
286 }
287 }
288 }
289 // A11
290 for (unsigned int q_index : fe_face.quadrature_point_indices())
291 {
292 for (unsigned int i = 0; i < dofs_per_cell; ++i)
293 {
294 for (unsigned int j = 0; j < dofs_per_cell; ++j)
295 {
296 A_11(i, j) +=
297 (-0.5 * fe_face1.shape_grad(i, q_index) * normals[q_index] *
298 fe_face1.shape_value(j, q_index) -
299 0.5 * fe_face1.shape_grad(j, q_index) * normals[q_index] *
300 fe_face1.shape_value(i, q_index) +
301 (penalty / hf) * fe_face1.shape_value(i, q_index) *
302 fe_face1.shape_value(j, q_index)) *
303 fe_face1.JxW(q_index);
304 }
305 }
306 }
307
308 // distribute DoFs accordingly
309 cell->get_dof_indices(local_dof_indices);
310 // Need to convert the neighbor to a valid "dh" iterator
311 typename DoFHandler<dim>::cell_iterator neigh_dh(*(cell->neighbor(f)),
312 &classical_dh);
313 // typename DoFHandler<dim>::cell_iterator neigh_dh(*(cell->neighbor(f)),
314 // &(ah->agglo_dh));
315 std::cout << "Neighbor is " << neigh_dh->active_cell_index() << std::endl;
316 neigh_dh->get_dof_indices(local_dof_indices1);
317
318 system_matrix.add(local_dof_indices, A_00);
319 system_matrix.add(local_dof_indices, local_dof_indices1, A_01);
320 system_matrix.add(local_dof_indices1, local_dof_indices, A_10);
321 system_matrix.add(local_dof_indices1, A_11);
322 }
323}
324
325template <int dim>
326void
327Poisson<dim>::assemble_system()
328{
329 DynamicSparsityPattern dsp(classical_dh.n_dofs());
330 DoFTools::make_flux_sparsity_pattern(classical_dh, dsp);
331 sparsity.copy_from(dsp);
332
333 system_matrix.reinit(sparsity);
334 solution.reinit(classical_dh.n_dofs());
335 system_rhs.reinit(classical_dh.n_dofs());
336
337 const unsigned int quadrature_degree = 3;
338 // ah->set_quadrature_degree(quadrature_degree);
339 // ah->set_agglomeration_flags(update_values | update_JxW_values |
340 // update_gradients | update_quadrature_points);
341 FEFaceValues<dim> fe_face(mapping,
342 dg_fe,
343 QGauss<dim - 1>(quadrature_degree),
344 update_values | update_JxW_values |
345 update_gradients | update_quadrature_points |
346 update_normal_vectors);
347
348
349 FEValues<dim> fe_values(mapping,
350 dg_fe,
351 QGauss<dim>(quadrature_degree),
352 update_values | update_JxW_values | update_gradients |
353 update_quadrature_points);
354
355 FEFaceValues<dim> fe_face1(mapping,
356 dg_fe,
357 QGauss<dim - 1>(quadrature_degree),
358 update_values | update_JxW_values |
359 update_gradients | update_quadrature_points |
360 update_normal_vectors);
361 const unsigned int dofs_per_cell = dg_fe.n_dofs_per_cell();
362
363 FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
364 Vector<double> cell_rhs(dofs_per_cell);
365
366 std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
367
368
369 // Master cells of the mesh
370 // for (const auto &cell :
371 // ah->agglo_dh.active_cell_iterators() |
372 // IteratorFilters::ActiveFEIndexEqualTo(ah->AggloIndex::master))
373 // {
374 // std::cout << "Cell with idx: " << cell->active_cell_index() <<
375 // std::endl; cell_matrix = 0; cell_rhs = 0;
376 // const auto &agglo_values = ah->reinit(cell);
377
378 // const auto &q_points = agglo_values.get_quadrature_points();
379 // const unsigned int n_qpoints = q_points.size();
380 // std::vector<double> rhs(n_qpoints);
381 // rhs_function->value_list(q_points, rhs);
382
383 // for (unsigned int q_index : agglo_values.quadrature_point_indices())
384 // {
385 // for (unsigned int i = 0; i < dofs_per_cell; ++i)
386 // {
387 // for (unsigned int j = 0; j < dofs_per_cell; ++j)
388 // {
389 // cell_matrix(i, j) += agglo_values.shape_grad(i, q_index) *
390 // agglo_values.shape_grad(j, q_index) *
391 // agglo_values.JxW(q_index);
392 // }
393 // cell_rhs(i) += agglo_values.shape_value(i, q_index) *
394 // rhs[q_index] * agglo_values.JxW(q_index);
395 // }
396 // }
397 // cell->get_dof_indices(local_dof_indices);
398
399 // // distribute DoFs
400 // for (unsigned int i = 0; i < dofs_per_cell; ++i)
401 // {
402 // for (unsigned int j = 0; j < dofs_per_cell; ++j)
403 // {
404 // system_matrix.add(local_dof_indices[i],
405 // local_dof_indices[j],
406 // cell_matrix(i, j));
407 // }
408 // system_rhs(local_dof_indices[i]) += cell_rhs(i);
409 // }
410
411
412
413 // // Face terms for master cells, i.e.: terms that are agglomerated
414 // const unsigned int n_agglomerated_faces =
415 // ah->n_faces(cell);
416 // std::cout << "Number of agglomerated faces per agglomeration: "
417 // << n_agglomerated_faces << std::endl;
418 // for (unsigned int f = 0; f < n_agglomerated_faces; ++f)
419 // {
420 // // We need to retrieve information about the deal.II face shared
421 // with
422 // // this agglomerated face. This can be asked to master_neighbors,
423 // that
424 // // takes precisely the master cell of the agglomeration and the
425 // index
426 // // of the **agglomerated** face.
427 // const auto &[local_face_idx,
428 // neighboring_cell,
429 // local_face_idx_outside] =
430 // ah->master_neighbors[{cell, f}];
431
432 // const auto &fe_isv = ah->reinit(cell, f);
433 // fe_face.reinit(neighboring_cell, local_face_idx_outside);
434
435 // std::vector<types::global_dof_index> local_dof_indices1(
436 // dofs_per_cell);
437
438 // const double hf =
439 // neighboring_cell->face(local_face_idx_outside)->measure();
440 // FullMatrix<double> A_00(dofs_per_cell, dofs_per_cell);
441 // FullMatrix<double> A_01(dofs_per_cell, dofs_per_cell);
442 // FullMatrix<double> A_10(dofs_per_cell, dofs_per_cell);
443 // FullMatrix<double> A_11(dofs_per_cell, dofs_per_cell);
444
445
446
447 // // Get normal vectors seen from each agglomeration.
448 // const auto &normals = fe_isv.get_normal_vectors();
449 // // A00
450 // for (unsigned int q_index : fe_isv.quadrature_point_indices())
451 // {
452 // std::cout << normals[q_index] << std::endl;
453 // for (unsigned int i = 0; i < dofs_per_cell; ++i)
454 // {
455 // for (unsigned int j = 0; j < dofs_per_cell; ++j)
456 // {
457 // A_00(i, j) +=
458 // (-0.5 * fe_isv.shape_grad(i, q_index) *
459 // normals[q_index] * fe_isv.shape_value(j, q_index)
460 // -
461 // 0.5 * fe_isv.shape_grad(j, q_index) *
462 // normals[q_index] * fe_isv.shape_value(i, q_index)
463 // +
464 // (penalty / hf) * fe_isv.shape_value(i, q_index) *
465 // fe_isv.shape_value(j, q_index)) *
466 // fe_isv.JxW(q_index);
467 // }
468 // }
469 // }
470 // // A10
471 // for (unsigned int q_index : fe_isv.quadrature_point_indices())
472 // {
473 // for (unsigned int i = 0; i < dofs_per_cell; ++i)
474 // {
475 // for (unsigned int j = 0; j < dofs_per_cell; ++j)
476 // {
477 // A_01(i, j) +=
478 // (-0.5 * fe_face.shape_grad(i, q_index) *
479 // normals[q_index] * fe_isv.shape_value(j, q_index)
480 // +
481 // 0.5 * fe_isv.shape_grad(j, q_index) *
482 // normals[q_index] * fe_face.shape_value(i, q_index)
483 // -
484 // (penalty / hf) * fe_face.shape_value(i, q_index) *
485 // fe_isv.shape_value(j, q_index)) *
486 // fe_isv.JxW(q_index);
487 // }
488 // }
489 // }
490 // // A01
491 // for (unsigned int q_index : fe_isv.quadrature_point_indices())
492 // {
493 // for (unsigned int i = 0; i < dofs_per_cell; ++i)
494 // {
495 // for (unsigned int j = 0; j < dofs_per_cell; ++j)
496 // {
497 // A_10(i, j) +=
498 // (-0.5 * fe_isv.shape_grad(i, q_index) *
499 // normals[q_index] * fe_face.shape_value(j, q_index)
500 // -
501 // 0.5 * fe_face.shape_grad(j, q_index) *
502 // normals[q_index] * fe_isv.shape_value(i, q_index)
503 // -
504 // (penalty / hf) * fe_isv.shape_value(i, q_index) *
505 // fe_face.shape_value(j, q_index)) *
506 // fe_isv.JxW(q_index);
507 // }
508 // }
509 // }
510 // // A11
511 // for (unsigned int q_index : fe_isv.quadrature_point_indices())
512 // {
513 // for (unsigned int i = 0; i < dofs_per_cell; ++i)
514 // {
515 // for (unsigned int j = 0; j < dofs_per_cell; ++j)
516 // {
517 // A_11(i, j) +=
518 // (-0.5 * fe_face.shape_grad(i, q_index) *
519 // normals[q_index] * fe_face.shape_value(j, q_index)
520 // -
521 // 0.5 * fe_face.shape_grad(j, q_index) *
522 // normals[q_index] * fe_face.shape_value(i, q_index)
523 // +
524 // (penalty / hf) * fe_face.shape_value(i, q_index) *
525 // fe_face.shape_value(j, q_index)) *
526 // fe_isv.JxW(q_index);
527 // }
528 // }
529 // }
530
531 // // distribute DoFs accordingly
532 // cell->get_dof_indices(local_dof_indices);
533 // // Need to convert the neighbor to a valid "dh" iterator
534 // typename DoFHandler<dim>::cell_iterator neigh_dh(*neighboring_cell,
535 // &(ah->agglo_dh));
536 // std::cout << "Neighbor is " << neigh_dh->active_cell_index()
537 // << std::endl;
538 // neigh_dh->get_dof_indices(local_dof_indices1);
539
540 // system_matrix.add(local_dof_indices, A_00);
541 // system_matrix.add(local_dof_indices, local_dof_indices1, A_01);
542 // system_matrix.add(local_dof_indices1, local_dof_indices, A_10);
543 // system_matrix.add(local_dof_indices1, A_11);
544 // }
545 // }
546
547
548
549 // Loop over standard deal.II cells
550 for (const auto &cell : classical_dh.active_cell_iterators())
551 {
552 std::cout << "Standard cell with idx: " << cell->active_cell_index()
553 << std::endl;
554 cell_matrix = 0;
555 cell_rhs = 0;
556 // const auto &agglo_values = ah->reinit(cell);
557 fe_values.reinit(cell);
558
559 const auto &q_points = fe_values.get_quadrature_points();
560 const unsigned int n_qpoints = q_points.size();
561 std::vector<double> rhs(n_qpoints);
562 rhs_function->value_list(q_points, rhs);
563
564 for (unsigned int q_index : fe_values.quadrature_point_indices())
565 {
566 for (unsigned int i = 0; i < dofs_per_cell; ++i)
567 {
568 for (unsigned int j = 0; j < dofs_per_cell; ++j)
569 {
570 cell_matrix(i, j) += fe_values.shape_grad(i, q_index) *
571 fe_values.shape_grad(j, q_index) *
572 fe_values.JxW(q_index);
573 }
574 cell_rhs(i) += fe_values.shape_value(i, q_index) * rhs[q_index] *
575 fe_values.JxW(q_index);
576 }
577 }
578
579 // distribute DoFs
580 cell->get_dof_indices(local_dof_indices);
581 for (unsigned int i = 0; i < dofs_per_cell; ++i)
582 {
583 for (unsigned int j = 0; j < dofs_per_cell; ++j)
584 {
585 system_matrix.add(local_dof_indices[i],
586 local_dof_indices[j],
587 cell_matrix(i, j));
588 }
589 system_rhs(local_dof_indices[i]) += cell_rhs(i);
590 }
591
592
593 // // Add the contributions from the faces of the standard deal.II mesh.
594 // const auto &inadmissible_faces_indices = ah->inadmissible_faces[cell];
595 // if (inadmissible_faces_indices.empty())
596 // {
597 // std::cout << "All faces are admissible for this cell" << std::endl;
598 // // all faces are admissible, you can add contributions.
599 // for (const auto f : cell->face_indices())
600 // {
601 // distribute_jumps_and_averages(fe_face, fe_face1, cell, f);
602 // }
603 // }
604 // else
605 // {
606 // Some faces are constrained, you need to loop only over the ones
607 // that are unconstrained.
608 for (const auto f : cell->face_indices())
609 {
610 // if (std::find(inadmissible_faces_indices.begin(),
611 // inadmissible_faces_indices.end(),
612 // f) == inadmissible_faces_indices.end())
613 // {
614 std::cout << "Face index " << f << " is admissible!" << std::endl;
615
616 distribute_jumps_and_averages(fe_face, fe_face1, cell, f);
617 // }
618 }
619 // }
620 }
621}
622
623
624template <int dim>
625void
626Poisson<dim>::solve()
627{
628 SparseDirectUMFPACK A_direct;
629 A_direct.initialize(system_matrix);
630 A_direct.vmult(solution, system_rhs);
631}
632
633
634
635template <int dim>
636void
637Poisson<dim>::output_results()
638{
639 const std::string filename = "agglomerated_Poisson.vtu";
640 std::ofstream output(filename);
641
642 DataOut<dim> data_out;
643 // data_out.attach_dof_handler(ah->agglo_dh);
644 data_out.attach_dof_handler(classical_dh);
645 data_out.add_data_vector(solution, "u", DataOut<dim>::type_dof_data);
646 data_out.build_patches(mapping);
647 data_out.write_vtu(output);
648}
649
650template <int dim>
651void
652Poisson<dim>::run()
653{
654 make_grid();
655 // setup_agglomeration();
656 assemble_system();
657 {
658 const std::string filename = "standar_matrix.txt";
659 std::ofstream output(filename);
660 system_matrix.print(output);
661 }
662 solve();
663 output_results();
664}
665
666int
667main()
668{
669 Poisson<2> poisson_problem;
670 poisson_problem.run();
671
672 return 0;
673}
agglomeration_handler.h
main
int main()
Definition assembly_check.cc:667
DataOut::write_vtu
void write_vtu(std::ostream &out) const
DataOut::attach_dof_handler
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
DataOut::add_data_vector
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation={})
DataOut::build_patches
virtual void build_patches(const unsigned int n_subdivisions=0)
FEFaceValues::reinit
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell, const unsigned int face_no)
FEValues::get_quadrature_points
const std::vector< Point< spacedim > > & get_quadrature_points() const
FEFaceValues::dofs_per_cell
const unsigned int dofs_per_cell
FEFaceValues::get_normal_vectors
const std::vector< Tensor< 1, spacedim > > & get_normal_vectors() const
FEFaceValues::quadrature_point_indices
std_cxx20::ranges::iota_view< unsigned int, unsigned int > quadrature_point_indices() const
FEFaceValues::shape_grad
const Tensor< 1, spacedim > & shape_grad(const unsigned int i, const unsigned int q_point) const
FEFaceValues::JxW
double JxW(const unsigned int q_point) const
FEFaceValues::shape_value
const double & shape_value(const unsigned int i, const unsigned int q_point) const
FEValues::reinit
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
FE_DGQ
GridOut::write_vtk
void write_vtk(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
SparseDirectUMFPACK::initialize
void initialize(const SparsityPattern &sparsity_pattern)
SparseDirectUMFPACK::vmult
void vmult(Vector< double > &dst, const Vector< double > &src) const
Poisson::distribute_jumps_and_averages
void distribute_jumps_and_averages(FEFaceValues< dim > &fe_face, FEFaceValues< dim > &fe_face1, const typename DoFHandler< dim >::active_cell_iterator &cell, const unsigned int f)
Definition assembly_check.cc:175
Poisson::mapping
MappingQ< dim > mapping
Definition assembly_check.cc:65
Poisson::sparsity
SparsityPattern sparsity
Definition assembly_check.cc:69
Poisson::system_matrix
SparseMatrix< double > system_matrix
Definition assembly_check.cc:70
Poisson::system_rhs
Vector< double > system_rhs
Definition assembly_check.cc:72
Poisson::output_results
void output_results()
Definition assembly_check.cc:637
Poisson::ah
std::unique_ptr< AgglomerationHandler< dim > > ah
Definition assembly_check.cc:68
Poisson::dg_fe
FE_DGQ< dim > dg_fe
Definition assembly_check.cc:66
Poisson::classical_dh
DoFHandler< dim > classical_dh
Definition assembly_check.cc:67
Poisson::setup_agglomeration
void setup_agglomeration()
Definition assembly_check.cc:109
Poisson::cached_tria
std::unique_ptr< GridTools::Cache< dim > > cached_tria
Definition assembly_check.cc:73
Poisson::tria
Triangulation< dim > tria
Definition assembly_check.cc:64
Poisson::assemble_system
void assemble_system()
Definition assembly_check.cc:327
Poisson::make_grid
void make_grid()
Definition assembly_check.cc:95
Poisson::rhs_function
std::unique_ptr< const Function< dim > > rhs_function
Definition assembly_check.cc:74
Poisson::solution
Vector< double > solution
Definition assembly_check.cc:71
Poisson::penalty
double penalty
Definition assembly_check.cc:81
RightHandSide::value_list
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int) const override
Definition assembly_check.cc:31
grid_generator.h
grid_out.h
DoFTools::make_flux_sparsity_pattern
void make_flux_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern)
update_values
update_values
update_normal_vectors
update_normal_vectors
update_JxW_values
update_JxW_values
update_gradients
update_gradients
update_quadrature_points
update_quadrature_points
GridGenerator::hyper_cube
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
dealii::PolyUtils::collect_cells_for_agglomeration
void collect_cells_for_agglomeration(const Triangulation< dim, spacedim > &tria, const std::vector< types::global_cell_index > &cell_idxs, std::vector< typename Triangulation< dim, spacedim >::active_cell_iterator > &cells_to_be_agglomerated)
Definition poly_utils.h:523
numbers::PI
static constexpr double PI
std::sin
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
data_out.h
sparse_direct.h
sparse_matrix.h