Thank you for responding Wolfgang Bangerth.
yes, I have compiled the code on the linux platform with deal.ii version
9.2.0; but when I run the program, there is a issure:
--------------------------------------------------------
An error occurred in line <1138> of file
</public3/home/sc52217/deal.ii_work/include/deal.II/dofs/dof_accessor.templates.h>
in function
static dealii::types::global_dof_index
dealii::internal::DoFAccessorImplementation::Implementation::get_vertex_dof_index(const
dealii::DoFHandler<dim, spacedim>&, unsigned int, unsigned int, unsigned
int) [with int dim = 2; int spacedim = 2; dealii::types::global_dof_index =
unsigned int]
The violated condition was:
static_cast<typename ::dealii::internal::argument_type<void( typename
std::common_type<decltype(local_index),
decltype(dof_handler.get_fe().dofs_per_vertex)>::type)>::
type>(local_index) < static_cast<typename
::dealii::internal::argument_type<void( typename
std::common_type<decltype(local_index),
decltype(dof_handler.get_fe().dofs_per_vertex)>::type)>::
type>(dof_handler.get_fe().dofs_per_vertex)
Additional information:
Index 1 is not in the half-open range [0,1).
Stacktrace:
-----------
#0 ./test: unsigned int
dealii::internal::DoFAccessorImplementation::Implementation::get_vertex_dof_index<2,
2>(dealii::DoFHandler<2, 2> const&, unsigned int, unsigned int, unsigned
int)
#1 ./test: dealii::DoFAccessor<2, dealii::DoFHandler<2, 2>,
false>::vertex_dof_index(unsigned int, unsigned int, unsigned int) const
#2 ./test: Step23::ModificationLaplaceEquation<2>::move_mesh()
#3 ./test: Step23::ModificationLaplaceEquation<2>::run()
#4 ./test: main
--------------------------------------------------------
I don't understand which index in the sentence: 'Index 1 is not in the
half-open range [0,1).' refer to.
there newest source is included at appendix.
Thank you,
Chen
在2020年7月8日星期三 UTC+8 上午12:38:40<Wolfgang Bangerth> 写道：
> On 7/6/20 11:25 PM, [email protected] wrote:
> > my platform is macOS Mojave 10.14.6 Beta, the deal.ii vesion is 9.1.0
>
> The functions you're trying to call (FEValues::quadrature_point_indices()
> and
> FEValues::dof_indices()) were only introduced in deal.II 9.2. You either
> need
> to upgrade to that version, or rely on functions that were already present
> in 9.1.
>
> Best
> W.
>
>
> --
> ------------------------------------------------------------------------
> Wolfgang Bangerth email: [email protected]
> www: http://www.math.colostate.edu/~bangerth/
>
>
--
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].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/42c19e9f-01c3-4c13-8484-8062b07908ffn%40googlegroups.com.
/* ---------------------------------------------------------------------
*
* Copyright (C) 2013 - 2019 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: Timo Heister, Texas A&M University, 2013
*/
// This tutorial program is odd in the sense that, unlike for most other
// steps, the introduction already provides most of the information on how to
// use the various strategies to generate meshes. Consequently, there is
// little that remains to be commented on here, and we intersperse the code
// with relatively little text. In essence, the code here simply provides a
// reference implementation of what has already been described in the
// introduction.
// @sect3{Include files}
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/manifold_lib.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/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/utilities.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.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 <iostream>
#include <fstream>
#include <map>
#include <vector>
#include <algorithm>
namespace Step23
{
using namespace dealii;
template <int dim>
class ModificationLaplaceEquation
{
public:
ModificationLaplaceEquation();
void run();
private:
void setup_system();
void solve_deltay();
void assemble_system();
void move_mesh();
void output_mesh();
Triangulation<dim> triangulation;
FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
AffineConstraints<double> constraints;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution_deltay, old_solution_deltay;
Vector<double> system_rhs;
double time_step;
double time;
unsigned int timestep_number;
};
template <int dim>
class InitialValuesU : public Function<dim>
{
public:
virtual double value(const Point<dim> & /*p*/,
const unsigned int component = 0) const override
{
(void)component;
Assert(component == 0, ExcIndexRange(component, 0, 1));
return 0;
}
};
template <int dim>
class InitialValuesDeltay : public Function<dim>
{
public:
virtual double value(const Point<dim> & /*p*/,
const unsigned int component = 0) const override
{
(void)component;
Assert(component == 0, ExcIndexRange(component, 0, 1));
return 0;
}
};
template <int dim>
class RightHandSide : public Function<dim>
{
public:
virtual double value(const Point<dim> & /*p*/,
const unsigned int component = 0) const override
{
(void)component;
Assert(component == 0, ExcIndexRange(component, 0, 1));
return 0;
}
};
template <int dim>
class BoundaryValuesDeltay : public Function<dim>
{
public:
virtual double value(const Point<dim> &/*p*/,
const unsigned int component = 0) const override
{
(void)component;
Assert(component == 0, ExcIndexRange(component, 0, 1));
return 0.5 * std::sin(this->get_time() * 2 * numbers::PI);
}
};
template <int dim>
class BoundaryValuesC : public Function<dim>
{
public:
virtual double value(const Point<dim> &/*p*/,
const unsigned int component = 0) const override
{
(void)component;
Assert(component == 0, ExcIndexRange(component, 0, 1));
return (std::cos(this->get_time() * 2 * numbers::PI) * numbers::PI);
}
};
template <int dim>
ModificationLaplaceEquation<dim>::ModificationLaplaceEquation()
: fe(1), dof_handler(triangulation), time_step(0.005), time(time_step), timestep_number(1)
{
}
template <int dim>
void ModificationLaplaceEquation<dim>::setup_system()
{
GridIn<dim> gridin;
gridin.attach_triangulation(triangulation);
std::ifstream f("test1.msh");
gridin.read_msh(f);
std::cout << "Number of active cells: " << triangulation.n_active_cells()
<< std::endl;
dof_handler.distribute_dofs(fe);
std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl
<< std::endl;
DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler, dsp);
sparsity_pattern.copy_from(dsp);
system_matrix.reinit(sparsity_pattern);
solution_deltay.reinit(dof_handler.n_dofs());
old_solution_deltay.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
assemble_system();
constraints.close();
}
template <int dim>
void ModificationLaplaceEquation<dim>::assemble_system()
{
QGauss<dim> quadrature_formula(fe.degree + 1); // 高斯插值阶，比有限元函数高一阶
FEValues<dim> fe_values(fe,
quadrature_formula,
update_values | update_gradients | update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
std::vector<double> Area;
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
for (const auto &cell : dof_handler.active_cell_iterators())
{
Area.emplace_back(cell->measure());
}
auto maxArea = *max_element(Area.begin(), Area.end());
auto minArea = *min_element(Area.begin(), Area.end());
for (const auto &cell : dof_handler.active_cell_iterators())
{
fe_values.reinit(cell);
cell_matrix = 0;
for (const unsigned int q_index : fe_values.quadrature_point_indices())
{
for (const unsigned int i : fe_values.dof_indices())
for (const unsigned int j : fe_values.dof_indices())
cell_matrix(i, j) += (1 + (maxArea - minArea)/cell->measure() ) *
(fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
fe_values.JxW(q_index)); // dx
}
cell->get_dof_indices(local_dof_indices);
for (const unsigned int i : fe_values.dof_indices())
for (const unsigned int j : fe_values.dof_indices())
system_matrix.add(local_dof_indices[i],
local_dof_indices[j],
cell_matrix(i, j));
}
std::map<types::global_dof_index, double> boundary_values;
BoundaryValuesDeltay<dim> boundary_values_deltay_function;
boundary_values_deltay_function.set_time(time);
VectorTools::interpolate_boundary_values(dof_handler,
0,
Functions::ZeroFunction<dim>(),
boundary_values);
VectorTools::interpolate_boundary_values(dof_handler,
1,
Functions::ZeroFunction<dim>(),
boundary_values);
VectorTools::interpolate_boundary_values(dof_handler,
2,
Functions::ZeroFunction<dim>(),
boundary_values);
VectorTools::interpolate_boundary_values(dof_handler,
3,
boundary_values_deltay_function,
boundary_values);
MatrixTools::apply_boundary_values(boundary_values,
system_matrix,
solution_deltay,
system_rhs);
}
template <int dim>
void ModificationLaplaceEquation<dim>::move_mesh()
{
std::cout << " Moving mesh..." << std::endl;
std::vector<bool> vertex_touched(triangulation.n_vertices(), false);
for (auto &cell : dof_handler.active_cell_iterators())
for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
if (vertex_touched[cell->vertex_index(v)] == false)
{
vertex_touched[cell->vertex_index(v)] = true;
Point<dim> vertex_displacement;
// for (unsigned int d = 0; d < dim; ++d)
vertex_displacement[0] = 0.;
vertex_displacement[1] = solution_deltay(cell->vertex_dof_index(v, 1));
cell->vertex(v) += vertex_displacement;
}
}
template <int dim>
void ModificationLaplaceEquation<dim>::output_mesh()
{
const std::string filename =
"mesh-" + Utilities::int_to_string(timestep_number, 6) + ".vtu";
std::ofstream out(filename);
GridOut grid_out;
grid_out.write_vtk(triangulation, out);
std::cout << " written to " << filename << std::endl
<< std::endl;
}
template <int dim>
void ModificationLaplaceEquation<dim>::solve_deltay()
{
SolverControl solver_control(1000, 1e-8 * system_rhs.l2_norm());
SolverCG<Vector<double>> cg(solver_control);
cg.solve(system_matrix, solution_deltay, system_rhs, PreconditionIdentity());
std::cout << " deltay-equation: " << solver_control.last_step()
<< " CG iterations." << std::endl;
}
template <int dim>
void ModificationLaplaceEquation<dim>::run()
{
setup_system();
VectorTools::project(dof_handler,
constraints,
QGauss<dim>(fe.degree + 1),
InitialValuesDeltay<dim>(),
old_solution_deltay);
VectorTools::project(dof_handler,
constraints,
QGauss<dim>(fe.degree + 1),
InitialValuesDeltay<dim>(),
solution_deltay);
Vector<double> tmp(solution_deltay.size());
Vector<double> forcing_terms(solution_deltay.size());
// for (; time <= 1; time += time_step, ++timestep_number)
// {
time = 0.25;
timestep_number = 50;
std::cout << "Time step " << timestep_number << " at t=" << time
<< std::endl;
assemble_system();
solve_deltay();
move_mesh();
output_mesh();
old_solution_deltay = solution_deltay;
// }
}
} // namespace Step23
int main()
{
try
{
using namespace Step23;
ModificationLaplaceEquation<2> ModificationLaplace_equation_solver;
ModificationLaplace_equation_solver.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;
}
