Hello dealii-team,
I have created a failing test case in serial when hanging nodes+periodic
constraints are mixed.
*The algorithm of the test is follows: (attached *
*minimalExampleHangingNodesPeriodicConstraintsBug.cc)*
1) Create a hypercube (-20,20) with origin at the center
2) Set periodic boundary conditions on the hypercube
3) Perform two levels of mesh refinement:- a) 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 two constraintMatrices: *constraints* and
*onlyHangingNodeConstraints.*
* constraints:* both hanging node and periodic constraints
*onlyHangingNodeConstraints:* just hanging node constraints
5) Create two nodal vectors vec1 and vec2 with nodal value =
nodal_coordinate.norm(). Set and distribute vec1 using *constraints. *Set
and distribute vec2 using *onlyHangingNodeConstraints.*
6) Print l2 norm of vec1 and vec2.
*The issue:*
Since the nodal value = nodal_coordinate.norm() is intrinsically periodic
(the domain being a hypercube with origin at the center), I would expect
vec1 and vec2 to have the same l2 norm. However I get different l2 norm
values:
L2 Norm vec1 (distributed using combined constraints): 178.494
L2 Norm vec2 (distributed using only hanging node constraints): 176.271
If I modify step 3) to: a) two step uniform refinement hypercube to get 64
cells ,and b) then pick one of the interior cells and refine only that cell
which creates hanging nodes only in the interior. Total (71 cells). The l2
norm of vec1 and vec2 matches in this case:
L2 Norm vec1 (distributed using combined constraints): 278.505
L2 Norm vec2 (distributed using only hanging node constraints): 278.505
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 dealii+unsubscr...@googlegroups.com.
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);
//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]-L)<1.0e-5)
cell->face(f)->set_boundary_id(1);
else if (std::abs(face_center[0]+L)<1.0e-5)
cell->face(f)->set_boundary_id(2);
else if (std::abs(face_center[1]-L)<1.0e-5)
cell->face(f)->set_boundary_id(3);
else if (std::abs(face_center[1]+L)<1.0e-5)
cell->face(f)->set_boundary_id(4);
else if (std::abs(face_center[2]-L)<1.0e-5)
cell->face(f)->set_boundary_id(5);
else if (std::abs(face_center[2]+L)<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);
std::cout << "Periodic Facepairs size: " << periodicity_vector.size() << std::endl;
//FIXME: two cases of mesh refinement, both having two levels of refinement to create hanging nodes
//Case 0 which has hanging node constraint on periodic face is the buggy case
//Case 1 which has hanging node constraints in the interior works fine
const unsigned int Case=0;
if (Case==0)
{
triangulation.refine_global(1);
//pick the first cell and refine
cell = triangulation.begin_active();
cell->set_refine_flag();
triangulation.execute_coarsening_and_refinement();
}
else if (Case==1)
{
triangulation.refine_global(2);
//pick an interior cell and refine again
cell = triangulation.begin_active();
endc = triangulation.end();
for (; cell!=endc; ++cell)
{
Point<3> cellCentroid=cell->center();
if (std::fabs(cellCentroid[0]-L/4.0)<1e-5
&& std::fabs(cellCentroid[1]-L/4.0)<1e-5
&& std::fabs(cellCentroid[2]-L/4.0)<1e-5)
cell->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 combined constraint matrix
ConstraintMatrix constraints;
constraints.clear();
DoFTools::make_hanging_node_constraints(dofHandler, constraints);
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);
///create only hanging node constraint matrix
ConstraintMatrix onlyHangingNodeConstraints;onlyHangingNodeConstraints.clear();
DoFTools::make_hanging_node_constraints(dofHandler, onlyHangingNodeConstraints);
onlyHangingNodeConstraints.close();
Vector<double> vec1;
vec1.reinit(dofHandler.n_dofs());
vec1=0;
Vector<double> vec2;
vec2.reinit(dofHandler.n_dofs());
vec2=0;
const unsigned int vertices_per_cell=GeometryInfo<3>::vertices_per_cell;
std::vector<bool> vertex_touched(dofHandler.get_triangulation().n_vertices(),
false);
DoFHandler<3>::active_cell_iterator
cellDof = dofHandler.begin_active(),
endcDof = dofHandler.end();
for (; cellDof!=endcDof; ++cellDof)
for (unsigned int i=0; i<vertices_per_cell; ++i)
{
const unsigned global_vertex_no = cellDof->vertex_index(i);
if (vertex_touched[global_vertex_no])
continue;
vertex_touched[global_vertex_no]=true;
Point<3> nodalCoor = cellDof->vertex(i);
const unsigned int globalDofIndex=cellDof->vertex_dof_index(i,0);
if(!constraints.is_constrained(globalDofIndex))
vec1[globalDofIndex]=nodalCoor.norm();
if(!onlyHangingNodeConstraints.is_constrained(globalDofIndex))
vec2[globalDofIndex]=nodalCoor.norm();
}
constraints.distribute(vec1);
std::cout<<"L2 Norm vec1 (distributed using combined constraints): "<<vec1.l2_norm()<<std::endl;
onlyHangingNodeConstraints.distribute(vec2);
std::cout<<"L2 Norm vec2 (distributed using only hanging node constraints): "<<vec2.l2_norm()<<std::endl;
}
