I tested the implementations with the example from step 12, but refinement
after one step ist not allowed (else I will get the error "ou are trying to
access the matrix entry with index <>", which means that either my code is
wrong, or something else is not working correctly. The question is: What is
not working correctly here...
Am Dienstag, 19. März 2019 11:11:08 UTC+1 schrieb Maxi Miller:
>
> Will take a look at it, thanks!
>
> Am Dienstag, 19. März 2019 10:59:32 UTC+1 schrieb Luca Heltai:
>>
>> > Is there a reason then that there are several examples using
>> integration_loop(), but (afaik) only one using mesh_loop?
>>
>> Yes. mesh_loop is more recent w.r.t. integration_loop.
>>
>> mesh_loop was introduced to address some of the oddities that are in
>> integration_loop, that make its use somewhat less intuitive w.r.t.
>> WorkStream::run, and was originally thought as the “kernel” of
>> integration_loop.
>>
>> The original idea was to replace the internals of integration_loop with
>> mesh_loop, but this stage of the project is not yet completed (for lack of
>> time of the original developers of mesh_loop… cough cough… :) …).
>>
>> The examples in
>>
>> https://github.com/dealii/dealii/pull/7806
>>
>> show you a few more ways to use mesh_loop, for DG methods, or using
>> Automatic Differentiation.
>>
>> Hopefully there will be more examples soon on using mesh_loop.
>>
>> L.
>>
>>
>>
--
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) 2009 - 2018 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.md at
* the top level directory of deal.II.
*
* ---------------------------------------------------------------------
*
* Author: Guido Kanschat, Texas A&M University, 2009
*/
// 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/multithread_info.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/matrix_out.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/fe/mapping_q1.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/lac/solver_richardson.h>
#include <deal.II/lac/precondition_block.h>
#include <deal.II/numerics/derivative_approximation.h>
#include <deal.II/meshworker/dof_info.h>
#include <deal.II/meshworker/integration_info.h>
#include <deal.II/meshworker/simple.h>
#include <deal.II/meshworker/loop.h>
#include <deal.II/meshworker/mesh_loop.h>
#include <copy_data.h>
#include <scratch_data.h>
// Like in all programs, we finish this section by including the needed C++
// headers and declaring we want to use objects in the dealii namespace without
// prefix.
#include <iostream>
#include <fstream>
namespace Step12
{
using namespace dealii;
// @sect3{Equation data}
//
// First, we define a class describing the inhomogeneous boundary data. Since
// only its values are used, we implement value_list(), but leave all other
// functions of Function undefined.
template <int dim>
class BoundaryValues : public Function<dim>
{
public:
BoundaryValues() = default;
virtual void value_list(const std::vector<Point<dim>> &points,
std::vector<double> & values,
const unsigned int component = 0) const override;
};
// Given the flow direction, the inflow boundary of the unit square $[0,1]^2$
// are the right and the lower boundaries. We prescribe discontinuous boundary
// values 1 and 0 on the x-axis and value 0 on the right boundary. The values
// of this function on the outflow boundaries will not be used within the DG
// scheme.
template <int dim>
void BoundaryValues<dim>::value_list(const std::vector<Point<dim>> &points,
std::vector<double> & values,
const unsigned int component) const
{
(void)component;
AssertIndexRange(component, 1);
Assert(values.size() == points.size(),
ExcDimensionMismatch(values.size(), points.size()));
for (unsigned int i = 0; i < values.size(); ++i)
{
if (points[i](0) < 0.5)
values[i] = 1.;
else
values[i] = 0.;
}
}
// Finally, a function that computes and returns the wind field
// $\beta=\beta(\mathbf x)$. As explained in the introduction, we will use a
// rotational field around the origin in 2d. In 3d, we simply leave the
// $z$-component unset (i.e., at zero), whereas the function can not be used
// in 1d in its current implementation:
template <int dim>
Tensor<1, dim> beta(const Point<dim> &p)
{
Assert(dim >= 2, ExcNotImplemented());
Point<dim> wind_field;
wind_field(0) = -p(1);
wind_field(1) = p(0);
wind_field /= wind_field.norm();
return wind_field;
}
template <int dim>
class AdvectionProblem
{
public:
AdvectionProblem();
void run();
private:
using ScratchData = MeshWorker::ScratchData<dim>;
using CopyData = MeshWorker::CopyData<1 + GeometryInfo<dim>::faces_per_cell, 1, 1 + GeometryInfo<dim>::faces_per_cell>;
void setup_system();
void assemble_system();
void solve(Vector<double> &solution);
void refine_grid();
void output_results(const unsigned int cycle) const;
void local_assemble_system(const typename DoFHandler<dim>::active_cell_iterator &cell,
ScratchData &scratch_data,
CopyData ©_data);
void local_assemble_boundary(const typename DoFHandler<dim>::active_cell_iterator &cell,
const unsigned int face_no,
ScratchData &scratch_data,
CopyData ©_data);
void local_assemble_faces(const typename DoFHandler<dim>::active_cell_iterator &interior_cell,
const unsigned int interior_face_no,
const unsigned int interior_subface_no,
const typename DoFHandler<dim>::active_cell_iterator &exterior_cell,
const unsigned int exterior_face_no,
const unsigned int exterior_subface_no,
ScratchData &scratch_data,
CopyData ©_data);
void copy_local_to_global(const CopyData ©_data);
AffineConstraints<double> hanging_node_constraints;
Triangulation<dim> triangulation;
const MappingQ1<dim> mapping;
UpdateFlags cell_flags = update_values | update_gradients | update_quadrature_points | update_JxW_values;
UpdateFlags face_flags = update_values | update_gradients | update_quadrature_points | update_normal_vectors | update_JxW_values;
FE_DGQ<dim> fe;
DoFHandler<dim> dof_handler;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> right_hand_side;
};
template <int dim>
AdvectionProblem<dim>::AdvectionProblem()
: mapping()
, fe(3)
, dof_handler(triangulation)
{}
template <int dim>
void AdvectionProblem<dim>::setup_system()
{
dof_handler.distribute_dofs(fe);
hanging_node_constraints.clear();
DoFTools::make_hanging_node_constraints(dof_handler,
hanging_node_constraints);
hanging_node_constraints.close();
DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
DoFTools::make_flux_sparsity_pattern(dof_handler, dsp);
sparsity_pattern.copy_from(dsp);
system_matrix.reinit(sparsity_pattern);
solution.reinit(dof_handler.n_dofs());
right_hand_side.reinit(dof_handler.n_dofs());
}
template <int dim>
void AdvectionProblem<dim>::assemble_system()
{
std::cout << "Using MeshWorker\n";
MeshWorker::mesh_loop(dof_handler.begin_active(),
dof_handler.end(),
std::bind(&AdvectionProblem<dim>::local_assemble_system,
this,
std::placeholders::_1,
std::placeholders::_2,
std::placeholders::_3),
std::bind(&AdvectionProblem<dim>::copy_local_to_global,
this,
std::placeholders::_1),
ScratchData(fe, QGauss<dim>(fe.degree + 1), cell_flags, QGauss<dim - 1>(fe.degree + 1), face_flags),
CopyData(fe.dofs_per_cell),
MeshWorker::AssembleFlags::assemble_own_cells | MeshWorker::AssembleFlags::assemble_boundary_faces | MeshWorker::AssembleFlags::assemble_own_interior_faces_both,
std::bind(&AdvectionProblem<dim>::local_assemble_boundary,
this,
std::placeholders::_1,
std::placeholders::_2,
std::placeholders::_3,
std::placeholders::_4),
std::bind(&AdvectionProblem<dim>::local_assemble_faces,
this,
std::placeholders::_1,
std::placeholders::_2,
std::placeholders::_3,
std::placeholders::_4,
std::placeholders::_5,
std::placeholders::_6,
std::placeholders::_7,
std::placeholders::_8)
);
}
template <int dim>
void AdvectionProblem<dim>::local_assemble_system(const typename DoFHandler<dim>::active_cell_iterator &cell,
ScratchData &scratch_data,
CopyData ©_data)
{
const auto &fev = scratch_data.reinit(cell);
const auto &JxW = scratch_data.get_JxW_values();
const auto &p = scratch_data.get_quadrature_points();
copy_data = 0;
copy_data.local_dof_indices[0] = scratch_data.get_local_dof_indices();
for(size_t q = 0; q < p.size(); ++q)
{
const Tensor<1, dim> beta_at_q_point = beta(p[q]);
for(size_t i = 0; i < fev.dofs_per_cell; ++i)
for(size_t j = 0; j < fev.dofs_per_cell; ++j)
copy_data.matrices[0](i, j) += -beta_at_q_point
* fev.shape_grad(i, q)
* fev.shape_value(j, q)
* JxW[q];
}
}
template <int dim>
void AdvectionProblem<dim>::local_assemble_boundary(const typename DoFHandler<dim>::active_cell_iterator &cell,
const unsigned int face_no,
ScratchData &scratch_data,
CopyData ©_data)
{
const auto &fev = scratch_data.reinit(cell, face_no);
const auto &JxW = scratch_data.get_JxW_values();
const auto &p = scratch_data.get_quadrature_points();
const auto &n = scratch_data.get_normal_vectors();
std::vector<double> g(p.size());
static BoundaryValues<dim> boundary_function;
boundary_function.value_list(p, g);
for (unsigned int q = 0; q < p.size(); ++q)
{
const double beta_dot_n = beta(p[q]) * n[q];
if (beta_dot_n > 0)
for (unsigned int i = 0; i < fev.dofs_per_cell; ++i)
for (unsigned int j = 0; j < fev.dofs_per_cell; ++j)
{
copy_data.matrices[0](i, j) += beta_dot_n
* fev.shape_value(j, q)
* fev.shape_value(i, q)
* JxW[q];
}
else
for (unsigned int i = 0; i < fev.dofs_per_cell; ++i)
copy_data.vectors[0](i) += -beta_dot_n
* g[q]
* fev.shape_value(i, q)
* JxW[q];
}
}
template <int dim>
void AdvectionProblem<dim>::local_assemble_faces(const typename DoFHandler<dim>::active_cell_iterator &interior_cell,
const unsigned int interior_face_no,
const unsigned int interior_subface_no,
const typename DoFHandler<dim>::active_cell_iterator &exterior_cell,
const unsigned int exterior_face_no,
const unsigned int exterior_subface_no,
ScratchData &scratch_data,
CopyData ©_data)
{
const auto &fev = scratch_data.reinit(interior_cell, interior_face_no, interior_subface_no);
const auto &JxW = scratch_data.get_JxW_values();
const auto &nfev = scratch_data.reinit_neighbour(exterior_cell, exterior_face_no, exterior_subface_no);
copy_data.local_dof_indices[interior_face_no + 1] = scratch_data.get_neighbor_dof_indices();
const auto &p = scratch_data.get_quadrature_points();
const auto &n = scratch_data.get_normal_vectors();
const auto &nn = scratch_data.get_neighbor_normal_vectors();
for(unsigned int q = 0; q < p.size(); ++q)
{
const auto beta_dot_n = beta(p[q]) * n[q];
if(beta_dot_n > 0)
{
for(unsigned int i = 0; i < fev.dofs_per_cell; ++i)
for(unsigned int j = 0; j < fev.dofs_per_cell; ++j)
{
copy_data.matrices[0](i, j) +=
beta_dot_n
* fev.shape_value(i, q)
* fev.shape_value(j, q)
* JxW[q];
}
}
else {
for(unsigned int i = 0; i < fev.dofs_per_cell; ++i)
for(unsigned int j = 0; j < nfev.dofs_per_cell; ++j)
{
copy_data.matrices[interior_face_no + 1](i, j) +=
beta_dot_n
* fev.shape_value(i, q)
* nfev.shape_value(j, q)
* JxW[q];
}
}
}
}
template <int dim>
void AdvectionProblem<dim>::copy_local_to_global(const CopyData ©_data)
{
hanging_node_constraints.distribute_local_to_global(copy_data.matrices[0],
copy_data.vectors[0],
copy_data.local_dof_indices[0],
system_matrix,
right_hand_side);
for(unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
hanging_node_constraints.distribute_local_to_global(copy_data.matrices[1 + f],
copy_data.local_dof_indices[0],
copy_data.local_dof_indices[1 + f],
system_matrix);
}
template <int dim>
void AdvectionProblem<dim>::solve(Vector<double> &solution)
{
SolverControl solver_control(1000, 1e-12);
SolverRichardson<> solver(solver_control);
// Here we create the preconditioner,
PreconditionBlockSSOR<SparseMatrix<double>> preconditioner;
// then assign the matrix to it and set the right block size:
preconditioner.initialize(system_matrix, fe.dofs_per_cell);
// After these preparations we are ready to start the linear solver.
solver.solve(system_matrix, solution, right_hand_side, preconditioner);
}
template <int dim>
void AdvectionProblem<dim>::refine_grid()
{
Vector<float> gradient_indicator(triangulation.n_active_cells());
DerivativeApproximation::approximate_gradient(mapping,
dof_handler,
solution,
gradient_indicator);
unsigned int cell_no = 0;
for (const auto &cell : dof_handler.active_cell_iterators())
gradient_indicator(cell_no++) *=
std::pow(cell->diameter(), 1 + 1.0 * dim / 2);
// Finally they serve as refinement indicator.
GridRefinement::refine_and_coarsen_fixed_number(triangulation,
gradient_indicator,
0.3,
0.1);
triangulation.execute_coarsening_and_refinement();
}
template <int dim>
void AdvectionProblem<dim>::output_results(const unsigned int cycle) const
{
const std::string filename = "sol-" + std::to_string(cycle) + ".vtu";
deallog << "Writing solution to <" << filename << ">" << std::endl;
std::ofstream vtu_output(filename);
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(solution, "u");
data_out.build_patches();
data_out.write_vtu(vtu_output);
}
template <int dim>
void AdvectionProblem<dim>::run()
{
for (unsigned int cycle = 0; cycle < 2; ++cycle)
{
deallog << "Cycle " << cycle << std::endl;
if (cycle == 0)
{
GridGenerator::hyper_cube(triangulation);
triangulation.refine_global(4);
//refine_grid();
}
else
refine_grid();
deallog << "Number of active cells: "
<< triangulation.n_active_cells() << std::endl;
setup_system();
//refine_grid();
std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl;
assemble_system();
solve(solution);
output_results(cycle);
}
}
}
int main()
{
try
{
Step12::AdvectionProblem<2> dgmethod;
dgmethod.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;
}
catch (...)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
