Hi all,
I have a very simple question but cannot think out:
I use step-6.cc to implement different adaptive mesh, and solve u with each
mesh. Then I want to compute :
[image: ask-7.JPG]
But there is a big gup on the results of the integration of the square of
grad u:
Cycle 0:
Number of active cells: 20
Number of degrees of freedom: 89
u square 0: 0.0305621 grad u square 0: 0.256542
Cycle 1:
Number of active cells: 44
Number of degrees of freedom: 209
u square 1: 0.0438936 grad u square 1: 1.31255
Cycle 2:
Number of active cells: 92
Number of degrees of freedom: 449
u square 2: 0.0505902 grad u square 2: 3.584
Cycle 3:
Number of active cells: 200
Number of degrees of freedom: 921
u square 3: 0.0536183 grad u square 3: 6.60495
Cycle 4:
Number of active cells: 440
Number of degrees of freedom: 2017
u square 4: 0.0555872 grad u square 4: 19.7488
Cycle 5:
Number of active cells: 956
Number of degrees of freedom: 4425
u square 5: 0.0562045 grad u square 5: 48.0924
Cycle 6:
Number of active cells: 1916
Number of degrees of freedom: 8993
u square 6: 0.0572988 grad u square 6: 117.562
Cycle 7:
Number of active cells: 3860
Number of degrees of freedom: 18353
u square 7: 0.0577537 grad u square 7: 258.867
And I don't know why the differences are so huge, and attach is my test
code.
Thanks in advance!
Best,
Chucui
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
/* ---------------------------------------------------------------------
*
* Copyright (C) 2000 - 2016 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE at
* the top level of the deal.II distribution.
*
* ---------------------------------------------------------------------
*
* Author: Wolfgang Bangerth, University of Heidelberg, 2000
*/
// @sect3{Include files}
// The first few files have already been covered in previous examples and will
// thus not be further commented on.
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/grid/tria.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <fstream>
#include <iostream>
#include <deal.II/fe/fe_q.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/numerics/error_estimator.h>
// Finally, this is as in previous programs:
using namespace dealii;
template <int dim>
class Step6
{
public:
Step6 ();
~Step6 ();
void run ();
private:
void setup_system ();
void assemble_system ();
void solve ();
void refine_grid ();
void output_results (const unsigned int cycle) const;
Triangulation<dim> triangulation;
DoFHandler<dim> dof_handler;
FE_Q<dim> fe;
ConstraintMatrix constraints;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
SparseMatrix<double> volume_matrix;
SparseMatrix<double> gradient_matrix;
Vector<double> solution;
Vector<double> system_rhs;
};
template <int dim>
double coefficient (const Point<dim> &p)
{
if (p.square() < 0.5*0.5)
return 20;
else
return 1;
}
template <int dim>
Step6<dim>::Step6 ()
:
dof_handler (triangulation),
fe (2)
{}
template <int dim>
Step6<dim>::~Step6 ()
{
dof_handler.clear ();
}
template <int dim>
void Step6<dim>::setup_system ()
{
dof_handler.distribute_dofs (fe);
solution.reinit (dof_handler.n_dofs());
system_rhs.reinit (dof_handler.n_dofs());
constraints.clear ();
DoFTools::make_hanging_node_constraints (dof_handler,
constraints);
VectorTools::interpolate_boundary_values (dof_handler,
0,
ZeroFunction<dim>(),
constraints);
constraints.close ();
DynamicSparsityPattern dsp(dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler,
dsp,
constraints,
/*keep_constrained_dofs = */ false);
sparsity_pattern.copy_from(dsp);
system_matrix.reinit (sparsity_pattern);
volume_matrix.reinit (sparsity_pattern);
gradient_matrix.reinit (sparsity_pattern);
}
template <int dim>
void Step6<dim>::assemble_system ()
{
const QGauss<dim> quadrature_formula(3);
FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
FullMatrix<double> cell_volume_matrix (dofs_per_cell, dofs_per_cell);
FullMatrix<double> cell_gradient_matrix (dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs (dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
cell_matrix = 0;
cell_volume_matrix = 0;
cell_gradient_matrix = 0;
cell_rhs = 0;
fe_values.reinit (cell);
for (unsigned int q_index=0; q_index<n_q_points; ++q_index)
{
const double current_coefficient = coefficient<dim>
(fe_values.quadrature_point (q_index));
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
cell_matrix(i,j) += (current_coefficient *
fe_values.shape_grad(i,q_index) *
fe_values.shape_grad(j,q_index) *
fe_values.JxW(q_index));
cell_volume_matrix(i,j) += (
fe_values.shape_value(i,q_index) *
fe_values.shape_value(j,q_index) *
fe_values.JxW(q_index));
cell_gradient_matrix(i,j) += (
fe_values.shape_grad(i,q_index) *
fe_values.shape_grad(j,q_index) *
fe_values.JxW(q_index));
}
cell_rhs(i) += (fe_values.shape_value(i,q_index) *
1.0 *
fe_values.JxW(q_index));
}
}
// Finally, transfer the contributions from @p cell_matrix and
// @p cell_rhs into the global objects.
cell->get_dof_indices (local_dof_indices);
constraints.distribute_local_to_global (cell_matrix,
cell_rhs,
local_dof_indices,
system_matrix,
system_rhs);
constraints.distribute_local_to_global (cell_volume_matrix,
local_dof_indices,
volume_matrix);
constraints.distribute_local_to_global (cell_gradient_matrix,
local_dof_indices,
gradient_matrix);
}
}
template <int dim>
void Step6<dim>::solve ()
{
SolverControl solver_control (1000, 1e-12);
SolverCG<> solver (solver_control);
PreconditionSSOR<> preconditioner;
preconditioner.initialize(system_matrix, 1.2);
solver.solve (system_matrix, solution, system_rhs,
preconditioner);
constraints.distribute (solution);
}
template <int dim>
void Step6<dim>::refine_grid ()
{
Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
KellyErrorEstimator<dim>::estimate (dof_handler,
QGauss<dim-1>(3),
typename FunctionMap<dim>::type(),
solution,
estimated_error_per_cell);
GridRefinement::refine_and_coarsen_fixed_number (triangulation,
estimated_error_per_cell,
0.3, 0.03);
triangulation.execute_coarsening_and_refinement ();
}
template <int dim>
void Step6<dim>::output_results (const unsigned int cycle) const
{
Assert (cycle < 10, ExcNotImplemented());
std::string filename = "grid-";
filename += ('0' + cycle);
filename += ".eps";
std::ofstream output (filename.c_str());
GridOut grid_out;
grid_out.write_eps (triangulation, output);
}
template <int dim>
void Step6<dim>::run ()
{
for (unsigned int cycle=0; cycle<8; ++cycle)
{
std::cout << "Cycle " << cycle << ':' << std::endl;
if (cycle == 0)
{
GridGenerator::hyper_ball (triangulation);
static const SphericalManifold<dim> boundary;
triangulation.set_all_manifold_ids_on_boundary(0);
triangulation.set_manifold (0, boundary);
triangulation.refine_global (1);
}
else
refine_grid ();
std::cout << " Number of active cells: "
<< triangulation.n_active_cells()
<< std::endl;
setup_system ();
std::cout << " Number of degrees of freedom: "
<< dof_handler.n_dofs()
<< std::endl;
assemble_system ();
solve ();
output_results (cycle);
double int_test_volume = volume_matrix.matrix_norm_square(solution);
double int_test_gradient = gradient_matrix.matrix_norm_square(solution);
std::cout << " u square " << cycle <<": " << int_test_volume << " grad u square " << cycle <<": " << int_test_gradient << std::endl;
}
DataOutBase::EpsFlags eps_flags;
eps_flags.z_scaling = 4;
DataOut<dim> data_out;
data_out.set_flags (eps_flags);
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (solution, "solution");
data_out.build_patches ();
std::ofstream output ("final-solution.eps");
data_out.write_eps (output);
}
int main ()
{
// The general idea behind the layout of this function is as follows: let's
// try to run the program as we did before...
try
{
Step6<2> laplace_problem_2d;
laplace_problem_2d.run ();
}
catch (std::exception &exc)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
// If the exception that was thrown somewhere was not an object of a class
// derived from the standard <code>exception</code> class, then we can't do
// anything at all. We then simply print an error message and exit.
catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
