Hi All,
I am facing the following issue: I am setting periodic boundary conditions
in all directions on a dofHandler object attached to a single element
triangulation. When I print the ConstraintMatrix I observe that a trivial
call to DoFTools::make_hanging_node_constraints(..) which I call prior to
calling DoFTools::make_periodicity_constraints gives an erroneous
ConstraintMatrix as shown below for the 8 degrees of freedom corresponding
to the 8 nodes.
0 1: 1
2 3: 1
4 5: 1
6 7: 1
1 3: 1
5 7: 1
3 7: 1
However, if the don't call ConstraintMatrix.close() (line 102 in the
minimal example) after the trivial call to
DoFTools::make_hanging_node_constraints(..), I get the correct
ConstraintMatrix. Likewise, if I don't use
make_hanging_node_constraints(..) the ConstraintMatrix is correct too
0 7: 1
1 7: 1
2 7: 1
3 7: 1
4 7: 1
5 7: 1
6 7: 1
I also checked that if I printed the ConstraintMatrix after just setting
the hanging_node_constraints, it doesn't print anything.
I wonder if I am making a mistake here. I have provided a minimal working
example file which reproduces this error.
Thanks,
Sambit
--
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].
For more options, visit https://groups.google.com/d/optout.
//Include all deal.II header file
#include <deal.II/base/quadrature.h>
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/timer.h>
#include <deal.II/base/table.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q1.h>
#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_generator.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/base/tensor_function.h>
#include <deal.II/base/point.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/utilities.h>
#include <deal.II/lac/lapack_full_matrix.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/base/parameter_handler.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/tria_boundary_lib.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/solver_bicgstab.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/lac/parallel_vector.h>
#include <deal.II/matrix_free/matrix_free.h>
#include <deal.II/matrix_free/fe_evaluation.h>
#include <deal.II/lac/slepc_solver.h>
#include <deal.II/base/config.h>
#include <deal.II/base/smartpointer.h>
//Include generic C++ headers
#include <fstream>
#include <iostream>
using namespace dealii;
int main (int argc, char *argv[])
{
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv);
const double L=20;
parallel::distributed::Triangulation<3> triangulation(MPI_COMM_WORLD);
GridGenerator::hyper_cube (triangulation, -L, L);
//mark faces
typename parallel::distributed::Triangulation<3>::active_cell_iterator cell = triangulation.begin_active(), endc = triangulation.end();
for(; cell!=endc; ++cell)
{
for(unsigned int f=0; f < GeometryInfo<3>::faces_per_cell; ++f)
{
const Point<3> face_center = cell->face(f)->center();
if(cell->face(f)->at_boundary())
{
if (std::abs(face_center[0]+(L))<1.0e-8)
cell->face(f)->set_boundary_id(1);
else if (std::abs(face_center[0]-(L))<1.0e-8)
cell->face(f)->set_boundary_id(2);
else if (std::abs(face_center[1]+(L))<1.0e-8)
cell->face(f)->set_boundary_id(3);
else if (std::abs(face_center[1]-(L))<1.0e-8)
cell->face(f)->set_boundary_id(4);
else if (std::abs(face_center[2]+(L))<1.0e-8)
cell->face(f)->set_boundary_id(5);
else if (std::abs(face_center[2]-(L))<1.0e-8)
cell->face(f)->set_boundary_id(6);
}
}
}
std::vector<GridTools::PeriodicFacePair<typename parallel::distributed::Triangulation<3>::cell_iterator> > periodicity_vector;
for (int i = 0; i < 3; ++i)
{
GridTools::collect_periodic_faces(triangulation, /*b_id1*/ 2*i+1, /*b_id2*/ 2*i+2,/*direction*/ i, periodicity_vector);
}
triangulation.add_periodicity(periodicity_vector);
FESystem<3> FE(FE_Q<3>(QGauss<1>(2)), 1); //linear shape function
DoFHandler<3> dofHandler (triangulation);
dofHandler.distribute_dofs(FE);
ConstraintMatrix perConstraints;
perConstraints.clear();
DoFTools::make_hanging_node_constraints(dofHandler, perConstraints);
perConstraints.close();
std::cout<< "Adding periodicity constraints to dofHandler "<< std::endl;
std::vector<GridTools::PeriodicFacePair<typename DoFHandler<3>::cell_iterator> > periodicity_vectorDof;
for (int i = 0; i < 3; ++i)
{
GridTools::collect_periodic_faces(dofHandler, /*b_id1*/ 2*i+1, /*b_id2*/ 2*i+2,/*direction*/ i, periodicity_vectorDof);
}
DoFTools::make_periodicity_constraints<DoFHandler<3> >(periodicity_vectorDof, perConstraints);
perConstraints.close();
std::cout<< "printing constraint matrix"<<std::endl;
perConstraints.print(std::cout);
std::cout << "number of elements: "
<< triangulation.n_global_active_cells()
<< std::endl;
}
