Deal.ii-step-3

Posted on | Edited on | In | Views: 
2019.09.23 Step-3 Study

The step-3 is the first for really calculating something in Deal.II and mainly talking about the solving process of the Poisson’s equation as an example.

We can run the default code and derive the results.

Default generated results

This example only considers the displacement in the x direction. So the DoF is the number of the nodes when using Q1 elements.

Now I will go through all the possible extensions suggested on the Deal.II online tutorial.

  1. Changing the mesh to hyper shell used previously and L-shape.
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
// hyper shell mesh, added by zzd
void Step3::make_grid()
{
	const Point<2> center(1, 0);
	const double   inner_radius = 0.5, outer_radius = 1.0;
	GridGenerator::hyper_shell(
			triangulation, center, inner_radius, outer_radius, 5);

	for (unsigned int step = 0; step < 3; ++step)
	{
		for (auto &cell : triangulation.active_cell_iterators())
			for (unsigned int v = 0; v < GeometryInfo<2>::vertices_per_cell; ++v)
			{
				const double distance_from_center =
						center.distance(cell->vertex(v));

				if (std::fabs(distance_from_center - inner_radius) < 1e-10)
				{
					cell->set_refine_flag();
					break;
				}
			}
		triangulation.execute_coarsening_and_refinement();
	}
}

Hyper_shell_results

We can see from the result that the problem we raised in the end of step-2 has some negative influence here. See the points of the mesh has been detached from its neighbours (see the sharp point at the top of the “cone” shape). Therefore, I need to further investigate this problem.

1
2
3
4
5
6
7
8
9
// hyper L shape, added by zzd
void Step3::make_grid()
{
    const unsigned int initial_global_refinement       = 4;
	const double   left = -1.0, right = 1.0;
	GridGenerator::hyper_L(
			triangulation, left, right, false);
	triangulation.refine_global(initial_global_refinement);
}

L_shape_results

The L-shape looks good here.

  1. Changing the boundary condition.
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    std::map<types::global_dof_index, double> boundary_values;
    	VectorTools::interpolate_boundary_values(dof_handler,
    			0,
    			ConstantFunction<2> (100),    // added by zzd
    			//Functions::ZeroFunction<2>(),
    			boundary_values);
    	MatrixTools::apply_boundary_values(boundary_values,
    			system_matrix,
    			solution,
    			system_rhs);

Changing the boundary conditions

  1. Modifying the type of the boundary conditions 
    1
    2
    3
    4
    5
    6
    7
    8
    void Step3::make_grid(int refinements)
    {
    	GridGenerator::hyper_cube(triangulation, -1, 1);
    	triangulation.begin_active()->face(0)->set_boundary_id(1); // added by zzd
    	triangulation.refine_global(refinements);
    	std::cout << "Number of active cells: " << triangulation.n_active_cells()
                		<< std::endl;
    }

Modifying the boundary conditions

4.Observe convergence

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
void Step3::output_results() const
{
	DataOut<2> data_out;
	data_out.attach_dof_handler(dof_handler);
	data_out.add_data_vector(solution, "solution");
	data_out.build_patches();
	std::ofstream output("solution.gpl");
	data_out.write_gnuplot(output);

	// observe the convergence, added by zzd
	std::cout << "Solution at (1/3,1/3): " \
			<< VectorTools::point_value (dof_handler, solution, \
					Point<2>(1./3, 1./3)) \
					<< std::endl;
}

Observing convergence

This is the convergence info. with 10 refinements.

We can see it is same as the 9th refinements, which imply the simulation has been converged successfully.

However, this mesh (1024*1024=104,8576, over 1 million) is too big to be shown by gnuplot on my laptop.

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