Dear Daniel,
Here is the small example you asked for. If you run it as it is, it will
work. If you can 1 to 2 in line 107, it crashes. Im run dealii-8.2.1, if
that is needed.
Best,
Joel
On Tuesday, November 1, 2016 at 11:44:48 AM UTC+1, Daniel Arndt wrote:
>
> Joel,
>
> I have code that runs fine with periodic boundary condition for
>> a GridGenerator::parallelepiped. Also when I change the mesh to
>> a GridGenerator::subdivided_parallelepiped and has only 1 as
>> n_subdivisions, it also works. But when I change it to 2, I get the
>> following error:
>> [...]
>>
> GridTools::collect_periodic_faces[1] tries to match the faces on the
> boundaries specified by the two boundary_ids on the coarsest level and here
> finds a face that can't be matched.
> Can you give us a minimal example just creating that mesh and trying to
> compute periodic boundary conditions?
> From just your description, it is hard to say what the problem is, but I
> would not expect that subdividing should be an issue.
>
> Best,
> Daniel
>
> [1]
> https://dealii.org/8.4.1/doxygen/deal.II/namespaceGridTools.html#aaeadfc0053429f542fbfd48d192b94f0
>
--
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) 1999 - 2013 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/grid/tria.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_dgq.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/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.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/compressed_sparsity_pattern.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/numerics/data_out.h>
#include <fstream>
#include <iostream>
#include <math.h>
#include <deal.II/base/logstream.h>
using namespace dealii;
template <int dim>
class Step4
{
public:
Step4 ();
void run ();
private:
void make_grid ();
void setup_system();
void assemble_system ();
void solve ();
//void output_results () const;
Triangulation<dim> triangulation;
FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
ConstraintMatrix constraints;
Vector<double> solution;
Vector<double> system_rhs;
};
// basis
Point<3> basis1(4.0,0.0,0.0);
Point<3> basis2(0.0,4.0,0.0);
Point<3> basis3(0.0,0.0,4.0);
template <int dim>
Step4<dim>::Step4 ()
:
fe (1),
dof_handler (triangulation)
{}
template <int dim>
void Step4<dim>::make_grid ()
{
Point<dim> basis [dim] = {basis1,basis2,basis3};
// change 1 to 2 -> fail
GridGenerator::subdivided_parallelepiped(triangulation,1,basis,true);
//GridGenerator::parallelepiped(triangulation,basis,true);
// Translate the main grid so that midpoint is the center
const Point<dim> translation = 0.5*(basis1+basis2+basis3);
GridTools::shift(-translation,triangulation);
// refinement
//triangulation.refine_global(2);
std::ofstream out ("./grid.vtk");
GridOut grid_out;
grid_out.write_vtk (triangulation, out);
std::cout << "Grid has been written to grid.vtk" << std::endl;
}
template <int dim>
void Step4<dim>::setup_system ()
{
dof_handler.distribute_dofs (fe);
std::cout << " Number of degrees of freedom: "
<< dof_handler.n_dofs()
<< std::endl;
constraints.clear ();
bool periodic = true;
bool x_dir = true;
bool y_dir = true;
bool z_dir = true;
if (periodic == true)
{
if (x_dir == true)
{
Tensor<1,dim> offset;
offset[0]=basis1[0];
offset[1]=basis1[1];
offset[2]=basis1[2];
DoFTools::make_periodicity_constraints(dof_handler,0,1,0,offset,constraints);
}
else if (x_dir == false)
{
VectorTools::interpolate_boundary_values (dof_handler,0,ZeroFunction<dim>(),constraints);
VectorTools::interpolate_boundary_values (dof_handler,1,ZeroFunction<dim>(),constraints);
}
if (y_dir == true)
{
Tensor<1,dim> offset;
offset[0]=basis2[0];
offset[1]=basis2[1];
offset[2]=basis2[2];
DoFTools::make_periodicity_constraints(dof_handler,2,3,1,offset,constraints);
}
else if (y_dir == false)
{
VectorTools::interpolate_boundary_values (dof_handler,2,ZeroFunction<dim>(),constraints);
VectorTools::interpolate_boundary_values (dof_handler,3,ZeroFunction<dim>(),constraints);
}
if (z_dir == true)
{
Tensor<1,dim> offset;
offset[0]=basis3[0];
offset[1]=basis3[1];
offset[2]=basis3[2];
DoFTools::make_periodicity_constraints(dof_handler,4,5,2,offset,constraints);
}
else if (z_dir == false)
{
VectorTools::interpolate_boundary_values (dof_handler,4,ZeroFunction<dim>(),constraints);
VectorTools::interpolate_boundary_values (dof_handler,5,ZeroFunction<dim>(),constraints);
}
}
else
{
VectorTools::interpolate_boundary_values (dof_handler,0,ZeroFunction<dim>(),constraints);
}
DoFTools::make_hanging_node_constraints (dof_handler,constraints);
constraints.close ();
solution.reinit (dof_handler.n_dofs());
system_rhs.reinit (dof_handler.n_dofs());
CompressedSparsityPattern c_sparsity(dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler,
c_sparsity,
constraints,
/*keep_constrained_dofs = */ false);
sparsity_pattern.copy_from(c_sparsity);
system_matrix.reinit (sparsity_pattern);
}
template <int dim>
void Step4<dim>::assemble_system ()
{
const 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);
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);
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) += (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);
constraints.distribute_local_to_global (cell_matrix,
cell_rhs,
local_dof_indices,
system_matrix,
system_rhs);
}
}
template <int dim>
void Step4<dim>::solve ()
{
SolverControl solver_control (1000, 1e-6);
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;
constraints.distribute (solution);
}
template <int dim>
void Step4<dim>::run ()
{
std::cout << "Solving problem in " << dim << " space dimensions." << std::endl;
make_grid();
setup_system ();
assemble_system ();
solve ();
}
int main ()
{
deallog.depth_console (0);
Step4<3> test;
test.run();
return 0;
}
