Hello,
I am facing an issue when I use hanging nodes with periodic boundary
conditions. I have reproduced the error in the attached minimal example
where I do the following steps:
1) Create a hypercube
2) Set appropriate boundary_ids on the faces of the hypercube for periodic
boundary conditions, call collect_periodic_faces and
triangulation.add_periodicity.
3) Perform two levels of mesh refinement:- a) first one step uniform
refinement hypercube to get 8 cells ,and b) then pick one of the 8 cells
and refine only that cell which creates hanging nodes on the faces. Finally
I get 15 cells (see attached image).
4) Create constraintMatrix which includes both hanging node and periodic
constraints.
*The issue:*
If I flip the boundary ids of the face pairs while marking the faces of
the hypercube in step-2 (change Ltemp=L in line 62 to Ltemp=-L), the size
of the constraint matrix created in step-4 changes. Further investigation
of the constraint matrix entries in the two cases revealed
that the number of identity_constraints are different- in one case it gives
19 identity_constraints and 17 in the other case. Manually counting the
numbering of identity_constraints using the attached image gives 19.
Is it correct to expect the size of the constraint matrix to remain
unchanged in the above case?
Additional note: The minimal example doesn't throw any error in the debug
mode. I am using deal ii version 8.5.1 and serial triangulation.
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;
Triangulation<3,3> triangulation;
GridGenerator::hyper_cube (triangulation, -L, L);
const double Ltemp=L;//FIXME: Size of constraint matrix changes between Ltemp=L, and Ltemp=-L;
//mark faces
typename Triangulation<3,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]-Ltemp)<1.0e-5)
cell->face(f)->set_boundary_id(1);
else if (std::abs(face_center[0]+Ltemp)<1.0e-5)
cell->face(f)->set_boundary_id(2);
else if (std::abs(face_center[1]-Ltemp)<1.0e-5)
cell->face(f)->set_boundary_id(3);
else if (std::abs(face_center[1]+Ltemp)<1.0e-5)
cell->face(f)->set_boundary_id(4);
else if (std::abs(face_center[2]-Ltemp)<1.0e-5)
cell->face(f)->set_boundary_id(5);
else if (std::abs(face_center[2]+Ltemp)<1.0e-5)
cell->face(f)->set_boundary_id(6);
}
}
}
std::vector<GridTools::PeriodicFacePair<typename Triangulation<3,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);
//two level mesh refinement
triangulation.refine_global(1);
typename Triangulation<3,3>::active_cell_iterator cellBegin = triangulation.begin_active();
cellBegin->set_refine_flag();
triangulation.execute_coarsening_and_refinement();
std::cout << "number of elements: "
<< triangulation.n_global_active_cells()
<< std::endl;
/////
FESystem<3> FE(FE_Q<3>(QGaussLobatto<1>(2)), 1); //linear shape function
DoFHandler<3> dofHandler (triangulation);
dofHandler.distribute_dofs(FE);
///creating constraint matrix
ConstraintMatrix constraints;
constraints.clear();
std::cout<< "Adding hanging node constraints... "<< std::endl;
DoFTools::make_hanging_node_constraints(dofHandler, constraints);
std::cout<< "Adding periodicity constraints... "<< 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, constraints);
constraints.close();
std::cout<< "size of constraint matrix: "<< constraints.n_constraints()<<std::endl;
std::cout<< "printing constraint matrix ...."<<std::endl;
constraints.print(std::cout);
/////write mesh
std::cout<<"writing mesh .." <<std::endl;
DataOut<3> data_out;
data_out.attach_dof_handler(dofHandler);
data_out.build_patches ();
std::ofstream output("mesh.vtu");
data_out.write_vtu (output);
}
