Dear Prof. Bangerth,
Thank you for your reply. My apologies for deleting my initial post as I
realized the mistake I was making right after posting it here. Basically I
was doing something like this
ConstraintMatrix perConstraints;
perConstraints.clear();
perConstraints.close();
...........then add periodic constraints
which according the description in the manual for close() is probably
wrong as that should be called at the end after setting all constraints. I
have attached the minimal working example of the same which now doesn't
include the DoFTools::make_hanging_node_constraints because that was not
causing the above issue.
Best,
Sambit
On Sunday, December 10, 2017 at 4:49:11 PM UTC-5, Wolfgang Bangerth wrote:
>
> On 12/10/2017 02:42 PM, Sambit Das wrote:
> >
> > I wanted to update DoFTools::make_hanging_node_constraints is not the
> issue
> > above. Even I comment out that line and call ConstraintMatrix.close()
> prior to
> > calling make_periodicity_constraints, I still get the same error. This
> means
> > calling ConstraintMatrix.close() on a constraintMatrix with no
> constraints is
> > leading to the error.
>
> That should be easy to test. Can you create a minimal working example that
> demonstrates this and that we can turn into a testcase after the bug is
> fixed?
>
> You've got a good start with your program, but keep removing things that
> are
> not necessary while making sure that the bug still appears. *Everything*
> that's not necessary should go to expose the problem as clearly as
> possible!
>
> Best
> W.
>
>
> --
> ------------------------------------------------------------------------
> Wolfgang Bangerth email: [email protected]
> <javascript:>
> www: http://www.math.colostate.edu/~bangerth/
>
>
--
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();
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;
}
