Hi, all.
I'd like ask a question implementing periodic boundary condition.
Modifying step-22 tutorial code, I want to solve a stokes flow at Couette
configuration
I plan to solve it at (x,y) in [0,5] x [0,1], rectangular domain
with boundary condition
\vec{u} (x,1)= (1,0) (moving dragging plate)
\vec{u} (x,0)= (0,0) (stationary bottom plate)
and periodic boundary condition at both inlet and outlet
\vec{u}(0,y)=\vec{u}(5,y)
I think in deal.ii, this can be implemented as
{
constraints.clear ();
FEValuesExtractors::Vector velocities(0);
DoFTools::make_hanging_node_constraints (dof_handler,
constraints);
// PERIODIC boundary condition
std::vector<GridTools::PeriodicFacePair
<typename DoFHandler<dim>::cell_iterator> > make_pair;
GridTools::collect_periodic_faces(dof_handler,1,3
/*direction*/,0
,make_pair);
VectorTools::interpolate_boundary_values (dof_handler,
4,
BoundaryValues<dim>(),
constraints,
fe.component_mask(velocities));
VectorTools::interpolate_boundary_values (dof_handler,
2,
ZeroFunction<dim>(dim+1),
constraints,
fe.component_mask(velocities));
}
However, as I run code, I get error message that
An error occurred in line <3751> of file
</Users/jaekwangkim/Program/dealii-8.5.0/source/grid/grid_tools.cc> in
function
void dealii::GridTools::collect_periodic_faces(const MeshType &, const
types::boundary_id, const types::boundary_id, const int,
std::vector<PeriodicFacePair<typename MeshType::cell_iterator> > &, const
Tensor<1, MeshType::space_dimension> &, const FullMatrix<double> &)
[MeshType = dealii::DoFHandler<2, 2>]
The violated condition was:
pairs1.size() == pairs2.size()
Additional information:
Unmatched faces on periodic boundaries
Stacktrace:
-----------
when it tried to compile following part
* GridTools::collect_periodic_faces(dof_handler,1,3*
* /*direction*/,0*
* ,make_pair);*
I would appreciate any comments or thoughts for possible reason of the
problem.
I attach minimalized problem that has same problem.
Thanks
--
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.
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/function.h>
#include <deal.II/base/utilities.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/block_sparse_matrix.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/grid/tria.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/tria_boundary_lib.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/grid_tools.h> //For the periodic boundary condition
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_renumbering.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_system.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/numerics/error_estimator.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/lac/sparse_ilu.h>
#include <iostream>
#include <fstream>
#include <sstream>
namespace Step22
{
using namespace dealii;
template <int dim>
struct InnerPreconditioner;
template <>
struct InnerPreconditioner<2>
{
typedef SparseDirectUMFPACK type;
};
template <>
struct InnerPreconditioner<3>
{
typedef SparseILU<double> type;
};
template <int dim>
class StokesProblem
{
public:
StokesProblem (const unsigned int degree);
void run ();
private:
void define_mesh ();
void setup_dofs ();
void assemble_system ();
void solve ();
void output_results ();
const unsigned int degree;
Triangulation<dim> triangulation;
FESystem<dim> fe;
DoFHandler<dim> dof_handler;
ConstraintMatrix constraints;
BlockSparsityPattern sparsity_pattern;
BlockSparseMatrix<double> system_matrix;
BlockVector<double> solution;
BlockVector<double> system_rhs;
std_cxx11::shared_ptr<typename InnerPreconditioner<dim>::type> A_preconditioner;
};
template <int dim>
class BoundaryValues : public Function<dim>
{
public:
BoundaryValues () : Function<dim>(dim+1) {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
virtual void vector_value (const Point<dim> &p,
Vector<double> &value) const;
};
template <int dim>
double
BoundaryValues<dim>::value (const Point<dim> &p,
const unsigned int component) const
{
Assert (component < this->n_components,
ExcIndexRange (component, 0, this->n_components));
if (component == 0)
return 1;
return 0;
}
template <int dim>
void
BoundaryValues<dim>::vector_value (const Point<dim> &p,
Vector<double> &values) const
{
for (unsigned int c=0; c<this->n_components; ++c)
values(c) = BoundaryValues<dim>::value (p, c);
}
template <int dim>
class RightHandSide : public Function<dim>
{
public:
RightHandSide () : Function<dim>(dim+1) {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
virtual void vector_value (const Point<dim> &p,
Vector<double> &value) const;
};
template <int dim>
double
RightHandSide<dim>::value (const Point<dim> &/*p*/,
const unsigned int /*component*/) const
{
return 0;
}
template <int dim>
void
RightHandSide<dim>::vector_value (const Point<dim> &p,
Vector<double> &values) const
{
for (unsigned int c=0; c<this->n_components; ++c)
values(c) = RightHandSide<dim>::value (p, c);
}
template <class MatrixType, class PreconditionerType>
class InverseMatrix : public Subscriptor
{
public:
InverseMatrix (const MatrixType &m,
const PreconditionerType &preconditioner);
void vmult (Vector<double> &dst,
const Vector<double> &src) const;
private:
const SmartPointer<const MatrixType> matrix;
const SmartPointer<const PreconditionerType> preconditioner;
};
template <class MatrixType, class PreconditionerType>
InverseMatrix<MatrixType,PreconditionerType>::InverseMatrix
(const MatrixType &m,
const PreconditionerType &preconditioner)
:
matrix (&m),
preconditioner (&preconditioner)
{}
template <class MatrixType, class PreconditionerType>
void InverseMatrix<MatrixType,PreconditionerType>::vmult
(Vector<double> &dst,
const Vector<double> &src) const
{
SolverControl solver_control (src.size(), 1e-6*src.l2_norm());
SolverCG<> cg (solver_control);
dst = 0;
cg.solve (*matrix, dst, src, *preconditioner);
}
template <class PreconditionerType>
class SchurComplement : public Subscriptor
{
public:
SchurComplement (const BlockSparseMatrix<double> &system_matrix,
const InverseMatrix<SparseMatrix<double>, PreconditionerType> &A_inverse);
void vmult (Vector<double> &dst,
const Vector<double> &src) const;
private:
const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
const SmartPointer<const InverseMatrix<SparseMatrix<double>, PreconditionerType> > A_inverse;
mutable Vector<double> tmp1, tmp2;
};
template <class PreconditionerType>
SchurComplement<PreconditionerType>::SchurComplement
(const BlockSparseMatrix<double> &system_matrix,
const InverseMatrix<SparseMatrix<double>,PreconditionerType> &A_inverse)
:
system_matrix (&system_matrix),
A_inverse (&A_inverse),
tmp1 (system_matrix.block(0,0).m()),
tmp2 (system_matrix.block(0,0).m())
{}
template <class PreconditionerType>
void SchurComplement<PreconditionerType>::vmult (Vector<double> &dst,
const Vector<double> &src) const
{
system_matrix->block(0,1).vmult (tmp1, src);
A_inverse->vmult (tmp2, tmp1);
system_matrix->block(1,0).vmult (dst, tmp2);
}
template <int dim>
StokesProblem<dim>::StokesProblem (const unsigned int degree)
:
degree (degree),
triangulation (Triangulation<dim>::maximum_smoothing),
fe (FE_Q<dim>(degree+1), dim,
FE_Q<dim>(degree), 1),
dof_handler (triangulation)
{}
template <int dim>
void StokesProblem<dim>::setup_dofs ()
{
A_preconditioner.reset ();
system_matrix.clear ();
dof_handler.distribute_dofs (fe);
DoFRenumbering::Cuthill_McKee (dof_handler);
std::vector<unsigned int> block_component (dim+1,0);
block_component[dim] = 1;
DoFRenumbering::component_wise (dof_handler, block_component);
{
constraints.clear ();
FEValuesExtractors::Vector velocities(0);
DoFTools::make_hanging_node_constraints (dof_handler,
constraints);
// PERIODIC boundary condition
std::vector<GridTools::PeriodicFacePair
<typename DoFHandler<dim>::cell_iterator> > make_pair;
GridTools::collect_periodic_faces(dof_handler,1,3
/*direction*/,0
,make_pair);
VectorTools::interpolate_boundary_values (dof_handler,
4,
BoundaryValues<dim>(),
constraints,
fe.component_mask(velocities));
VectorTools::interpolate_boundary_values (dof_handler,
2,
ZeroFunction<dim>(dim+1),
constraints,
fe.component_mask(velocities));
}
constraints.close ();
std::vector<types::global_dof_index> dofs_per_block (2);
DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, block_component);
const unsigned int n_u = dofs_per_block[0], n_p = dofs_per_block[1];
{
BlockDynamicSparsityPattern dsp (2,2);
dsp.block(0,0).reinit (n_u, n_u);
dsp.block(1,0).reinit (n_p, n_u);
dsp.block(0,1).reinit (n_u, n_p);
dsp.block(1,1).reinit (n_p, n_p);
dsp.collect_sizes();
DoFTools::make_sparsity_pattern (dof_handler, dsp, constraints, false);
sparsity_pattern.copy_from (dsp);
}
system_matrix.reinit (sparsity_pattern);
solution.reinit (2);
solution.block(0).reinit (n_u);
solution.block(1).reinit (n_p);
solution.collect_sizes ();
system_rhs.reinit (2);
system_rhs.block(0).reinit (n_u);
system_rhs.block(1).reinit (n_p);
system_rhs.collect_sizes ();
}
template <int dim>
void StokesProblem<dim>::assemble_system ()
{
system_matrix=0;
system_rhs=0;
QGauss<dim> quadrature_formula(degree+2);
FEValues<dim> fe_values (fe, quadrature_formula,
update_values |
update_quadrature_points |
update_JxW_values |
update_gradients);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
FullMatrix<double> local_matrix (dofs_per_cell, dofs_per_cell);
Vector<double> local_rhs (dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
const RightHandSide<dim> right_hand_side;
std::vector<Vector<double> > rhs_values (n_q_points,
Vector<double>(dim+1));
const FEValuesExtractors::Vector velocities (0);
const FEValuesExtractors::Scalar pressure (dim);
std::vector<SymmetricTensor<2,dim> > symgrad_phi_u (dofs_per_cell);
std::vector<double> div_phi_u (dofs_per_cell);
std::vector<double> phi_p (dofs_per_cell);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
fe_values.reinit (cell);
local_matrix = 0;
local_rhs = 0;
right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
rhs_values);
for (unsigned int q=0; q<n_q_points; ++q)
{
for (unsigned int k=0; k<dofs_per_cell; ++k)
{
symgrad_phi_u[k] = fe_values[velocities].symmetric_gradient (k, q);
div_phi_u[k] = fe_values[velocities].divergence (k, q);
phi_p[k] = fe_values[pressure].value (k, q);
}
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
for (unsigned int j=0; j<=i; ++j)
{
local_matrix(i,j) += (2 * (symgrad_phi_u[i] * symgrad_phi_u[j])
- div_phi_u[i] * phi_p[j]
- phi_p[i] * div_phi_u[j]
+ phi_p[i] * phi_p[j])
* fe_values.JxW(q);
}
const unsigned int component_i =
fe.system_to_component_index(i).first;
local_rhs(i) += fe_values.shape_value(i,q) *
rhs_values[q](component_i) *
fe_values.JxW(q);
}
}
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=i+1; j<dofs_per_cell; ++j)
local_matrix(i,j) = local_matrix(j,i);
cell->get_dof_indices (local_dof_indices);
constraints.distribute_local_to_global (local_matrix, local_rhs,
local_dof_indices,
system_matrix, system_rhs);
}
A_preconditioner
= std_cxx11::shared_ptr<typename InnerPreconditioner<dim>::type>(new typename InnerPreconditioner<dim>::type());
A_preconditioner->initialize (system_matrix.block(0,0),
typename InnerPreconditioner<dim>::type::AdditionalData());
}
template <int dim>
void StokesProblem<dim>::solve ()
{
const InverseMatrix<SparseMatrix<double>,
typename InnerPreconditioner<dim>::type>
A_inverse (system_matrix.block(0,0), *A_preconditioner);
Vector<double> tmp (solution.block(0).size());
{
Vector<double> schur_rhs (solution.block(1).size());
A_inverse.vmult (tmp, system_rhs.block(0));
system_matrix.block(1,0).vmult (schur_rhs, tmp);
schur_rhs -= system_rhs.block(1);
SchurComplement<typename InnerPreconditioner<dim>::type>
schur_complement (system_matrix, A_inverse);
SolverControl solver_control (solution.block(1).size(),
1e-6*schur_rhs.l2_norm());
SolverCG<> cg (solver_control);
SparseILU<double> preconditioner;
preconditioner.initialize (system_matrix.block(1,1),
SparseILU<double>::AdditionalData());
InverseMatrix<SparseMatrix<double>,SparseILU<double> >
m_inverse (system_matrix.block(1,1), preconditioner);
cg.solve (schur_complement, solution.block(1), schur_rhs,
m_inverse);
constraints.distribute (solution);
}
{
system_matrix.block(0,1).vmult (tmp, solution.block(1));
tmp *= -1;
tmp += system_rhs.block(0);
A_inverse.vmult (solution.block(0), tmp);
constraints.distribute (solution);
}
}
template <int dim>
void
StokesProblem<dim>::output_results ()
{
std::vector<std::string> solution_names (dim, "velocity");
solution_names.push_back ("pressure");
std::vector<DataComponentInterpretation::DataComponentInterpretation>
data_component_interpretation
(dim, DataComponentInterpretation::component_is_part_of_vector);
data_component_interpretation
.push_back (DataComponentInterpretation::component_is_scalar);
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (solution, solution_names,
DataOut<dim>::type_dof_data,
data_component_interpretation);
data_out.build_patches ();
std::ostringstream filename;
filename << "solution.vtk";
std::ofstream output (filename.str().c_str());
data_out.write_vtk (output);
}
template <int dim>
void StokesProblem<dim>::define_mesh ()
{
const Point<2> end (5.0,1.0);
const Point<2> start (0.0,0.0);
GridGenerator::hyper_rectangle (triangulation, start, end, false);
for (typename Triangulation<dim>::active_cell_iterator
cell=triangulation.begin_active();
cell!=triangulation.end(); ++cell)
{
for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
{
// You need this line to ensure that this doesn't affect anything other than bondary faces!
if (cell->face(f)->at_boundary() == true)
{
if ( std::abs(cell->face(f)->center()[0]-0.5)< 1e-10 )
{
cell->face(f)->set_boundary_id(10);
}
// Boundary faces
static const double tol = 1e-12;
if ( std::abs(cell->face(f)->center()[0]-0.0)< tol )
{
cell->face(f)->set_boundary_id(1); // Will be periodic to bd_id=3
}
if ( std::abs(cell->face(f)->center()[0]-1.0)< tol )
{
cell->face(f)->set_boundary_id(3); // Will be periodic with bd_id=1
}
if ( std::abs(cell->face(f)->center()[1]-0.0)< tol )
{
cell->face(f)->set_boundary_id(2); // Dirichlet BC
}
if ( std::abs(cell->face(f)->center()[1]-1.0)< tol )
{
cell->face(f)->set_boundary_id(4); // Dirichlet BC
}
}
}
} //end of cycle over cells
triangulation.refine_global (2);
}
template <int dim>
void StokesProblem<dim>::run ()
{
define_mesh();
setup_dofs ();
assemble_system ();
solve ();
output_results ();
}
}
int main ()
{
try
{
using namespace dealii;
using namespace Step22;
StokesProblem<2> flow_problem(1);
flow_problem.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;
}
