Deal.ii-step-6

Posted on | Edited on | In | Views: 
2019.09.28 Step-6 Study

The step-6 is mainly talking about how deal.ii is doing adaptive mesh refinement.

We can run the default code and derive the results.

Default generated result 1
Default generated result 2
Default generated result 3
Default generated result 4
Default generated result 5
Default generated result 6
Default generated result 7
Default generated result 8

Now I will go through some the possible extensions thought out by myself.

  1. Computational time test
1
2
PreconditionSSOR<> preconditioner;
preconditioner.initialization(system_matrix, 1.2);

Default SSOR preconditioner-16s.

1
2
PreconditionSSOR<> preconditioner;
preconditioner.initialization(system_matrix, 1.0);

Default SSOR preconditioner with 1.0 parameter-16s.

1
2
PreconditionSSOR<> preconditioner;
preconditioner.initialization(system_matrix, 0.5);

Default SSOR preconditioner with parameter 0.12-23s.

1
2
PreconditionJacobi<> preconditioner;  // Jacobi as a preconditioner
preconditioner.initialize(system_matrix);

Jacobi preconditioner-9s

1
2
SparseILU<double> preconditioner;
preconditioner.initialize(system_matrix);

Incomplete LU decomposition-11 s.

Jacobi preconditioner is the fastest!

  1. Better mesh
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    template<int dim>
    void Step6<dim>::make_grid()
    {
    	GridGenerator::hyper_ball (triangulation);
    	// after GridGenerator::hyper_ball is called the Triangulation has
    	// a SphericalManifold with id 0. We can use it again on the interior.
    	const Point<dim> mesh_center;
    	for (const auto &cell : triangulation.active_cell_iterators())
    		if (mesh_center.distance (cell->center()) > cell->diameter()/10)
    			cell->set_all_manifold_ids (0);
    	triangulation.refine_global (1);
    	//	GridGenerator::hyper_ball(triangulation);
    	//	triangulation.refine_global(1);
    	//	//	triangulation.execute_coarsening_and_refinement(/*1*/);
    }

Mesh crack problem
In addition, I found no improvement. Moreover, I didn’t understand the meaning of the following

1
2
if (mesh_center.distance (cell->center()) > cell->diameter()/10)
			cell->set_all_manifold_ids (0);

Further improved mesh:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
template<int dim>
void Step6<dim>::make_grid()
{
	GridGenerator::hyper_ball (triangulation);
	const Point<dim> mesh_center;
	const double core_radius  = 1.0/5.0,
			inner_radius = 1.0/3.0;
	// Step 1: Shrink the inner cell
	//
	// We cannot get a circle out of the inner cell because of
	// the degeneration problem mentioned above. Rather, shrink
	// the inner cell to a core radius of 1/5 that stays
	// sufficiently far away from the place where the
	// coefficient will have a discontinuity and where we want
	// to have cell interfaces that actually lie on a circle.
	// We do this shrinking by just scaling the location of each
	// of the vertices, given that the center of the circle is
	// simply the origin of the coordinate system.
	for (const auto &cell : triangulation.active_cell_iterators())
		if (mesh_center.distance (cell->center()) < 1e-5)
		{
			for (unsigned int v=0;
					v < GeometryInfo<dim>::vertices_per_cell;
					++v)
				cell->vertex(v) *= core_radius/mesh_center.distance (cell->vertex(v));
		}
	// Step 2: Refine all cells except the central one
	for (const auto &cell : triangulation.active_cell_iterators())
		if (mesh_center.distance (cell->center()) >= 1e-5)
			cell->set_refine_flag ();
	triangulation.execute_coarsening_and_refinement ();
	// Step 3: Resize the inner children of the outer cells
	//
	// The previous step replaced each of the four outer cells
	// by its four children, but the radial distance at which we
	// have intersected is not what we want to later refinement
	// steps. Consequently, move the vertices that were just
	// created in radial direction to a place where we need
	// them.
	for (const auto &cell : triangulation.active_cell_iterators())
		for (unsigned int v=0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
		{
			const double dist = mesh_center.distance (cell->vertex(v));
			if (dist > core_radius*1.0001 && dist < 0.9999)
				cell->vertex(v) *= inner_radius/dist;
		}
	// Step 4: Apply curved manifold description
	//
	// As discussed above, we can not expect to subdivide the
	// inner four cells (or their faces) onto concentric rings,
	// but we can do so for all other cells that are located
	// outside the inner radius. To this end, we loop over all
	// cells and determine whether it is in this zone. If it
	// isn't, then we set the manifold description of the cell
	// and all of its bounding faces to the one that describes
	// the spherical manifold already introduced above and that
	// will be used for all further mesh refinement.
	for (const auto &cell : triangulation.active_cell_iterators())
	{
		bool is_in_inner_circle = false;
		for (unsigned int v=0; v < GeometryInfo<2>::vertices_per_cell; ++v)
			if (mesh_center.distance (cell->vertex(v)) < inner_radius)
			{
				is_in_inner_circle = true;
				break;
			}
		if (is_in_inner_circle == false)
			// The Triangulation already has a SphericalManifold with
			// manifold id 0 (see the documentation of
			// GridGenerator::hyper_ball) so we just attach it to the outer
			// ring here:
			cell->set_all_manifold_ids (0);
	}

}

//template<int dim>
//void Step6<dim>::make_grid()
//{
//	//	GridGenerator::hyper_ball (triangulation);
//	//	// after GridGenerator::hyper_ball is called the Triangulation has
//	//	// a SphericalManifold with id 0. We can use it again on the interior.
//	//	const Point<dim> mesh_center;
//	//	for (const auto &cell : triangulation.active_cell_iterators())
//	//		if (mesh_center.distance (cell->center()) > cell->diameter()/10)
//	//			cell->set_all_manifold_ids (0);
//	//	triangulation.refine_global (3);
//	GridGenerator::hyper_ball(triangulation);
//	triangulation.refine_global(1);
//	//	triangulation.execute_coarsening_and_refinement(/*1*/);
//}

Further improved mesh

The mesh is still not ideal in my opinion. I must miss something.

  1. Play with coefficient
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    template<int dim>
    double coefficient (const Point<dim> &p)
    {
    	if ((p[0] < 0) && (p[1] < 0))
    		return 1;
    	else if ((p[0] >= 0) && (p[1] < 0))
    		return 10;
    	else if ((p[0] < 0) && (p[1] >= 0))
    		return 100;
    	else if ((p[0] >= 0) && (p[1] >= 0))
    		return 1000;
    	else
    	{
    		Assert (false, ExcInternalError());
    		return 0;
    	}
    }

Play with coefficient 6

Future work:

  1. The FEM information flow control should recieve more attention, such as how exactly the boundary_values are applied on the final linear equations system.

  2. The “friend” and “virtual” in C++ is still not crystal clear and should be further studied.

  3. I don’t understand how the error estimation inequality is derived.
    I don't understand how the error estimation inequality is derived.

  4. Manifold unclear either.

The official tutorial document is as follows,
https://www.dealii.org/current/doxygen/deal.II/step_6.html