Here, I attached simple case that I encounter same error message. Most of code is just same at that of step-5, but I am using mesh-file that is made from other software. (I attach my mesh file too)
At constructor, there is a issue.
template <int dim>
Step5<dim>::Step5 () :
triangulation (Triangulation<dim>::maximum_smoothing),* //If I don't have
this line, the code works well with anisotropy refinement. *
*However, if I use Maximum_smoothing for my triangulation, I run into
same error message.*
fe (1),
dof_handler (triangulation)
{}
2017년 3월 31일 금요일 오후 6시 8분 16초 UTC-5, Wolfgang Bangerth 님의 말:
>
> On 03/31/2017 01:30 PM, Jaekwang Kim wrote:
> >
> > I found that what was the problem.
> >
> > I was using 'maximum' Meshsmoothing which does not allow
> >
> > So, It was not a bug of deal.ii but I was using wrong smoothing method
> > that does not go along with anisotropy,
>
> Well, it is still a bug that the library did not handle this in a more
> gracefully way, i.e., that you did not get a polite error that suggests
> what is happening here and why.
>
> So I would still be interested in having a small testcase that shows the
> problem and that would allow us to look into this some more...
>
> Best
> W.
>
> --
> ------------------------------------------------------------------------
> Wolfgang Bangerth email: bang...@colostate.edu
> <javascript:>
> 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 dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
##
# CMake script for the step-5 tutorial program:
##
# Set the name of the project and target:
SET(TARGET "step-5")
# Declare all source files the target consists of. Here, this is only
# the one step-X.cc file, but as you expand your project you may wish
# to add other source files as well. If your project becomes much larger,
# you may want to either replace the following statement by something like
# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc")
# FILE(GLOB_RECURSE TARGET_INC "include/*.h")
# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
# or switch altogether to the large project CMakeLists.txt file discussed
# in the "CMake in user projects" page accessible from the "User info"
# page of the documentation.
SET(TARGET_SRC
${TARGET}.cc
)
# Usually, you will not need to modify anything beyond this point...
CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)
FIND_PACKAGE(deal.II 8.4 QUIET
HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
)
IF(NOT ${deal.II_FOUND})
MESSAGE(FATAL_ERROR "\n"
"*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
"You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to
cmake\n"
"or set an environment variable \"DEAL_II_DIR\" that contains this path."
)
ENDIF()
DEAL_II_INITIALIZE_CACHED_VARIABLES()
PROJECT(${TARGET})
DEAL_II_INVOKE_AUTOPILOT()
mesh_stretch.inp
Description: chemical/gamess-input
/* ---------------------------------------------------------------------
*
* Copyright (C) 1999 - 2015 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, 1999
*/
#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/grid/grid_refinement.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.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/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/fe/mapping_q1.h>
#include <deal.II/numerics/derivative_approximation.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/manifold_lib.h>
#include <fstream>
#include <iostream>
#include <sstream>
using namespace dealii;
template <int dim>
class Step5
{
public:
Step5 ();
void run ();
private:
void setup_system ();
void assemble_system ();
void solve ();
void output_results (const unsigned int cycle) const;
void refine_grid (const unsigned int cycle);
void set_anisotropic_flags ();
Triangulation<dim> triangulation;
const MappingQ1<dim> mapping;
FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> system_rhs;
};
template <int dim>
class Coefficient : public Function<dim>
{
public:
Coefficient () : Function<dim>() {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
virtual void value_list (const std::vector<Point<dim> > &points,
std::vector<double> &values,
const unsigned int component = 0) const;
};
template <int dim>
double Coefficient<dim>::value (const Point<dim> &p,
const unsigned int /*component*/) const
{
if (p.square() < 0.5*0.5)
return 20;
else
return 1;
}
template <int dim>
void Coefficient<dim>::value_list (const std::vector<Point<dim> > &points,
std::vector<double> &values,
const unsigned int component) const
{
Assert (values.size() == points.size(),
ExcDimensionMismatch (values.size(), points.size()));
Assert (component == 0,
ExcIndexRange (component, 0, 1));
const unsigned int n_points = points.size();
for (unsigned int i=0; i<n_points; ++i)
{
if (points[i].square() < 0.5*0.5)
values[i] = 20;
else
values[i] = 1;
}
}
template <int dim>
Step5<dim>::Step5 () :
triangulation (Triangulation<dim>::maximum_smoothing),
fe (1),
dof_handler (triangulation)
{}
template <int dim>
void Step5<dim>::setup_system ()
{
dof_handler.distribute_dofs (fe);
std::cout << " Number of degrees of freedom: "
<< dof_handler.n_dofs()
<< std::endl;
DynamicSparsityPattern dsp(dof_handler.n_dofs());
DoFTools::make_sparsity_pattern (dof_handler, dsp);
sparsity_pattern.copy_from(dsp);
system_matrix.reinit (sparsity_pattern);
solution.reinit (dof_handler.n_dofs());
system_rhs.reinit (dof_handler.n_dofs());
}
template <int dim>
void Step5<dim>::assemble_system ()
{
QGauss<dim> quadrature_formula(2);
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);
Vector<double> cell_rhs (dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
const Coefficient<dim> coefficient;
std::vector<double> coefficient_values (n_q_points);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
cell_matrix = 0;
cell_rhs = 0;
fe_values.reinit (cell);
coefficient.value_list (fe_values.get_quadrature_points(),
coefficient_values);
for (unsigned int q_index=0; q_index<n_q_points; ++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) += (coefficient_values[q_index] *
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));
}
cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
for (unsigned int j=0; j<dofs_per_cell; ++j)
system_matrix.add (local_dof_indices[i],
local_dof_indices[j],
cell_matrix(i,j));
system_rhs(local_dof_indices[i]) += cell_rhs(i);
}
}
std::map<types::global_dof_index,double> boundary_values;
VectorTools::interpolate_boundary_values (dof_handler,
1,
ConstantFunction<dim>(1),
boundary_values);
VectorTools::interpolate_boundary_values (dof_handler,
2,
ConstantFunction<dim>(1),
boundary_values);
VectorTools::interpolate_boundary_values (dof_handler,
3,
ZeroFunction<dim>(),
boundary_values);
VectorTools::interpolate_boundary_values (dof_handler,
4,
ConstantFunction<dim>(1),
boundary_values);
VectorTools::interpolate_boundary_values (dof_handler,
5,
ZeroFunction<dim>(),
boundary_values);
MatrixTools::apply_boundary_values (boundary_values,
system_matrix,
solution,
system_rhs);
}
template <int dim>
void Step5<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);
std::cout << " " << solver_control.last_step()
<< " CG iterations needed to obtain convergence."
<< std::endl;
}
template <int dim>
void Step5<dim>::refine_grid (const unsigned int cycle)
{
//triangulation.refine_global (1);
if (cycle <1 )
{
triangulation.refine_global (1);
}
else
{
Vector<float> gradient_indicator (triangulation.n_active_cells());
DerivativeApproximation::approximate_gradient (mapping,
dof_handler,
solution,
gradient_indicator);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no)
gradient_indicator(cell_no)*=std::pow(cell->diameter(), 1+1.0*dim/2);
GridRefinement::refine_and_coarsen_fixed_number (triangulation,
gradient_indicator,
0.3, 0.0);
set_anisotropic_flags();
triangulation.execute_coarsening_and_refinement ();
}
}
template <int dim>
void Step5<dim>::set_anisotropic_flags ()
{
const double anisotropic_threshold_ratio=1.0;
std::cout << " anisotropic refinement is being conducted : " << std::endl;
QGauss<dim> quadrature (2);
QGauss<dim-1> face_quadrature (2);
UpdateFlags face_update_flags
= UpdateFlags(update_values | update_JxW_values);
FEFaceValues<dim> fe_v_face (mapping, fe, face_quadrature, face_update_flags);
FESubfaceValues<dim> fe_v_subface (mapping, fe, face_quadrature, face_update_flags);
FEFaceValues<dim> fe_v_face_neighbor (mapping, fe, face_quadrature, update_values);
typename DoFHandler<dim>::active_cell_iterator cell=dof_handler.begin_active(),
endc=dof_handler.end();
for (; cell!=endc; ++cell)
if (cell->refine_flag_set()) // Refine Flagged 셀에만 각 면별로 점프를 계산해서, 어떤 면을 Refine할 것인지를 결정한다.
{
Point<dim> jump;
Point<dim> area;
for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
{
typename DoFHandler<dim>::face_iterator face = cell->face(face_no);
if (!face->at_boundary())
{
Assert (cell->neighbor(face_no).state() == IteratorState::valid, ExcInternalError());
typename DoFHandler<dim>::cell_iterator neighbor = cell->neighbor(face_no);
std::vector<double> u (fe_v_face.n_quadrature_points);
std::vector<double> u_neighbor (fe_v_face.n_quadrature_points);
if (face->has_children())
{
unsigned int neighbor2=cell->neighbor_face_no(face_no);
for (unsigned int subface_no=0; subface_no<face->number_of_children(); ++subface_no)
{
typename DoFHandler<dim>::cell_iterator neighbor_child = cell->neighbor_child_on_subface(face_no,subface_no);
Assert (!neighbor_child->has_children(), ExcInternalError());
fe_v_subface.reinit (cell, face_no, subface_no);
fe_v_face_neighbor.reinit (neighbor_child, neighbor2);
fe_v_subface.get_function_values(solution, u);
fe_v_face_neighbor.get_function_values(solution, u_neighbor);
const std::vector<double> &JxW = fe_v_subface.get_JxW_values ();
for (unsigned int x=0; x<fe_v_subface.n_quadrature_points; ++x)
{
jump[face_no/2]+=std::fabs(u[x]-u_neighbor[x])*JxW[x];
area[face_no/2]+=JxW[x];
}
}
} // 단순히 해당 면에서의 점프를 계산하는 과정이다
else
{
if (!cell->neighbor_is_coarser(face_no))
{
unsigned int neighbor2=cell->neighbor_of_neighbor(face_no);
fe_v_face.reinit (cell, face_no);
fe_v_face_neighbor.reinit (neighbor, neighbor2);
fe_v_face.get_function_values(solution, u);
fe_v_face_neighbor.get_function_values(solution, u_neighbor);
const std::vector<double> &JxW = fe_v_face.get_JxW_values ();
for (unsigned int x=0; x<fe_v_face.n_quadrature_points; ++x)
{
jump[face_no/2]+=std::fabs(u[x]-u_neighbor[x])*JxW[x];
area[face_no/2]+=JxW[x];
}
}
else //i.e. neighbor is coarser than cell
{
std::pair<unsigned int,unsigned int> neighbor_face_subface
= cell->neighbor_of_coarser_neighbor(face_no);
Assert (neighbor_face_subface.first<GeometryInfo<dim>::faces_per_cell, ExcInternalError());
Assert (neighbor_face_subface.second<neighbor->face(neighbor_face_subface.first)->number_of_children(),
ExcInternalError());
Assert (neighbor->neighbor_child_on_subface(neighbor_face_subface.first, neighbor_face_subface.second)
== cell, ExcInternalError());
fe_v_face.reinit (cell, face_no);
fe_v_subface.reinit (neighbor, neighbor_face_subface.first,
neighbor_face_subface.second);
fe_v_face.get_function_values(solution, u);
fe_v_subface.get_function_values(solution, u_neighbor);
const std::vector<double> &JxW = fe_v_face.get_JxW_values ();
for (unsigned int x=0; x<fe_v_face.n_quadrature_points; ++x)
{
jump[face_no/2]+=std::fabs(u[x]-u_neighbor[x])*JxW[x];
area[face_no/2]+=JxW[x];
}
}
}
}
}
double average_jumps[dim];
double sum_of_average_jumps=0.;
for (unsigned int i=0; i<dim; ++i)
{
average_jumps[i] = jump(i)/area(i);
sum_of_average_jumps += average_jumps[i];
}
for (unsigned int i=0; i<dim; ++i)
{
if (average_jumps[i] > anisotropic_threshold_ratio*(sum_of_average_jumps-average_jumps[i]))
cell->set_refine_flag(RefinementCase<dim>::cut_axis(i));
}
}
}
template <int dim>
void Step5<dim>::output_results (const unsigned int cycle) const
{
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (solution, "solution");
data_out.build_patches ();
DataOutBase::EpsFlags eps_flags;
eps_flags.z_scaling = 4;
eps_flags.azimut_angle = 40;
eps_flags.turn_angle = 10;
data_out.set_flags (eps_flags);
std::ostringstream filename;
filename << "solution-"
<< cycle
<< ".vtk";
std::ofstream output (filename.str().c_str());
data_out.write_vtk (output);
}
template <int dim>
void Step5<dim>::run ()
{
{
std::vector<unsigned int> subdivisions (dim, 1);
subdivisions[0] = 4;
GridIn<dim> grid_in;
grid_in.attach_triangulation (triangulation);
std::ifstream input_file("mesh_stretch.inp");
Assert (dim==2, ExcInternalError());
grid_in.read_ucd (input_file);
const Point<2> center (0,0);
static const SphericalManifold<dim> manifold_description;
triangulation.set_manifold (1, manifold_description);
}
for (unsigned int cycle=0; cycle<8; ++cycle)
{
std::cout << "Cycle " << cycle << ':' << std::endl;
if (cycle != 0)
refine_grid (cycle);
std::cout << " Number of active cells: "
<< triangulation.n_active_cells()
<< std::endl
<< " Total number of cells: "
<< triangulation.n_cells()
<< std::endl;
setup_system ();
assemble_system ();
solve ();
output_results (cycle);
}
}
int main ()
{
Step5<2> laplace_problem_2d;
laplace_problem_2d.run ();
/*
Coefficient<2> coefficient;
std::vector<Point<2> > points (2);
std::vector<double> coefficient_values (1);
coefficient.value_list (points, coefficient_values);
*/
return 0;
}
