The attach is the code for refine global.
在 2019年2月11日星期一 UTC+8上午11:21:00，[email protected]写道：
>
> Dear Prof. Bangerth,
>>
>
> Thanks for your quick reply!
>
> > You mean you are wondering why the "grad u square" term grows with the
> size of the problem but "u square" does not?
>
> Yes, that's what I am confused. Our u_exact have exact "grad u_exact
> square" and "u_exact square". And our numerical solution "u" and its
> "grad u square" and "u square" should be approximation of the exact one
> though with adaptive mesh. So I think "grad u square" should not have so
> huge jump with different adaptive meshes, because they are all the
> approximation of "grad u_exact square".
>
> > I don't know either (and the one place I looked at in your code seems
> correct to me), but what happens if you do global refinement? Does it
> increase with a fixed ratio from one step to the next?
>
> I change the code in attaching, and the output is:
> 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: 80
> Number of degrees of freedom: 337
> u square 1: 0.0489011 grad u square 1: 0.341164
> Cycle 2:
> Number of active cells: 320
> Number of degrees of freedom: 1313
> u square 2: 0.0540201 grad u square 2: 0.361384
> Cycle 3:
> Number of active cells: 1280
> Number of degrees of freedom: 5185
> u square 3: 0.0549874 grad u square 3: 0.365557
> Cycle 4:
> Number of active cells: 5120
> Number of degrees of freedom: 20609
> u square 4: 0.0553762 grad u square 4: 0.367069
> Cycle 5:
> Number of active cells: 20480
> Number of degrees of freedom: 82177
> u square 5: 0.0555421 grad u square 5: 0.367712
> Cycle 6:
> Number of active cells: 81920
> Number of degrees of freedom: 328193
> u square 6: 0.0556186 grad u square 6: 0.367986
> Cycle 7:
> Number of active cells: 327680
> Number of degrees of freedom: 1311745
> u square 7: 0.0556576 grad u square 7: 0.368119
>
>
> and the ratio
>
> [image: ask-8-1.JPG]
> are:
> u_square: 1.8410 2.4038 1.3149 1.2287 1.1168
> 0.9720
> grad_u_square: 2.0652 2.2766 1.4646 1.2336 1.2843 0.9891
>
> And the results without adaptive mesh seem like approximate the exact
> results correctly. So I am confused why the results of "grad u square" with
> adaptive mesh are so strange? Maybe the reason is that we refinre mesh
> where the gradient of u is huge, so only "grad_u_square" can be very
> different with changing adaptive meshes? If so, how to deal with this
> problem? How to approximate "grad u_ecaxt square" correctly with changing
> adaptive mesh? Maybe change the value smaller of scale of refine cells and
> coarsen cells, to make the change smaller?
>
> Thank you very much!
>
> 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 [email protected].
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 (1000000, 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 ();
triangulation.refine_global (1);
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;
}
