Hello,
I am trying to use the parallel::distributed::SolutionTransfer<dim,
LinearAlgebra::distributed::Vector<double>, hp::DoFHandler<dim>> class and
I can't seem to use it similarly to the serial version.
I looked through the tests/distributed_grids/solution_transfer_0*.cc tests
and none of them seem to be testing for the hp refinement. The code I'm
working on might be a bit large. Instead, I recreated the errors by copying
the tests/hp/solution_transfer.cc, which contains SolutionTransfer for hp,
and modified it to use the parallel distributed version of it.
Please find attached a proposed test. Line 160-175 are commented out and
does the p-refinement. Therefore, the error below is for h-refinement only,
but with an hp object. Oddly enough, my code works fine for h-refinement,
but I get the same error below when I do p-refinement.
110: --------------------------------------------------------
110: An error occurred in line <401> of file
</home/ddong/Libraries/dealii/source/distributed/solution_transfer.cc> in
function
110: void dealii::parallel::distributed::SolutionTransfer<dim,
VectorType, DoFHandlerType>::unpack_callback(const typename
dealii::parallel::distributed::Triangulation<dim, DoFHandlerType::
space_dimension>::cell_iterator&, typename
dealii::parallel::distributed::Triangulation<dim, DoFHandlerType::
space_dimension>::CellStatus, const
boost::iterator_range<__gnu_cxx::__normal_iterator<const char*,
std::vector<char, std::allocator<char> > > >&, std::vector<VectorType*>&)
[with int dim = 2; VectorType =
dealii::LinearAlgebra::distributed::Vector<double>; DoFHandlerType =
dealii::hp::DoFHandler<2, 2>; typename
dealii::parallel::distributed::Triangulation<dim, DoFHandlerType::
space_dimension>::cell_iterator =
dealii::TriaIterator<dealii::CellAccessor<2, 2> >; typename
dealii::parallel::distributed::Triangulation<dim, DoFHandlerType::
space_dimension>::CellStatus = dealii::Triangulation<2, 2>::CellStatus]
110: The violated condition was:
110: dof_handler->get_fe(fe_index).dofs_per_cell == it_dofvalues->size()
110: Additional information:
110: The transferred data was packed with a different number of dofs
than thecurrently registered FE object assigned to the DoFHandler has.
110: --------------------------------------------------------
Now, if I uncomment those lines, and therefore do p-refinement in the
attached test, I get
110: --------------------------------------------------------
110: An error occurred in line <382> of file
</home/ddong/Libraries/dealii/source/distributed/solution_transfer.cc> in
function
110: void dealii::parallel::distributed::SolutionTransfer<dim,
VectorType, DoFHandlerType>::unpack_callback(const typename
dealii::parallel::distributed::Triangulation<dim, DoFHandlerType::
space_dimension>::cell_iterator&, typename
dealii::parallel::distributed::Triangulation<dim, DoFHandlerType::
space_dimension>::CellStatus, const
boost::iterator_range<__gnu_cxx::__normal_iterator<const char*,
std::vector<char, std::allocator<char> > > >&, std::vector<VectorType*>&)
[with int dim = 2; VectorType =
dealii::LinearAlgebra::distributed::Vector<double>; DoFHandlerType =
dealii::hp::DoFHandler<2, 2>; typename
dealii::parallel::distributed::Triangulation<dim, DoFHandlerType::
space_dimension>::cell_iterator =
dealii::TriaIterator<dealii::CellAccessor<2, 2> >; typename
dealii::parallel::distributed::Triangulation<dim, DoFHandlerType::
space_dimension>::CellStatus = dealii::Triangulation<2, 2>::CellStatus]
110: The violated condition was:
110: cell->child(child_index)->active_fe_index() == fe_index
110: Additional information:
110: This exception -- which is used in many places in the library --
usually indicates that some condition which the author of the code thought
must be satisfied at a certain point in an algorithm, is not fulfilled. An
example would be that the first part of an algorithm sorts elements of an
array in ascending order, and a second part of the algorithm later
encounters an element that is not larger than the previous one.
110:
110: There is usually not very much you can do if you encounter such an
exception since it indicates an error in deal.II, not in your own program.
Try to come up with the smallest possible program that still demonstrates
the error and contact the deal.II mailing lists with it to obtain help.
110: --------------------------------------------------------
I went through some of the deal.II source code and it's a bit out of my
depth. Please let me know if I am somehow not p-refining correctly, or if I
can be of any help to fix this.
Best regards,
Doug
--
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.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/1836ac96-4f6c-4d90-ad2f-5e68cbabe623%40googlegroups.com.
// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2018 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.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------
#include <deal.II/base/function.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/mapping_q1.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/hp/dof_handler.h>
#include <deal.II/hp/fe_collection.h>
#include <deal.II/hp/fe_values.h>
#include <deal.II/lac/vector.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/distributed/solution_transfer.h>
#include <deal.II/distributed/tria.h>
#include <iostream>
#include <vector>
#include "../tests.h"
// a linear function that should be transferred exactly with Q1 and Q2
// elements
template <int dim>
class MyFunction : public Function<dim>
{
public:
MyFunction()
: Function<dim>(){};
virtual double
value(const Point<dim> &p, const unsigned int) const
{
double f = 0.25 + 2 * p[0];
if (dim > 1)
f += 0.772 * p[1];
if (dim > 2)
f -= 3.112 * p[2];
return f;
};
};
template <int dim>
void
transfer(std::ostream &out)
{
MyFunction<dim> function;
parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
GridGenerator::hyper_cube(tria);
tria.refine_global(5 - dim);
const unsigned int max_degree = 6 - dim;
hp::FECollection<dim> fe_q;
hp::FECollection<dim> fe_dgq;
for (unsigned int deg = 1; deg <= max_degree; ++deg)
{
fe_q.push_back(FE_Q<dim>(deg));
fe_dgq.push_back(FE_DGQ<dim>(deg));
}
hp::DoFHandler<dim> q_dof_handler(tria);
hp::DoFHandler<dim> dgq_dof_handler(tria);
LinearAlgebra::distributed::Vector<double> (q_solution);
LinearAlgebra::distributed::Vector<double> (dgq_solution);
MappingQGeneric<dim> mapping(1);
// refine a few cells
typename Triangulation<dim>::active_cell_iterator cell = tria.begin_active(),
endc = tria.end();
++cell;
++cell;
for (; cell != endc; ++cell)
cell->set_refine_flag();
tria.prepare_coarsening_and_refinement();
tria.execute_coarsening_and_refinement();
// randomly assign FE orders
unsigned int counter = 0;
{
typename hp::DoFHandler<dim>::active_cell_iterator
cell = q_dof_handler.begin_active(),
celldg = dgq_dof_handler.begin_active(), endc = q_dof_handler.end();
for (; cell != endc; ++cell, ++celldg, ++counter)
{
if (counter < 15)
cell->set_active_fe_index(1);
else
cell->set_active_fe_index(Testing::rand() % max_degree);
if (counter < 15)
celldg->set_active_fe_index(1);
else
celldg->set_active_fe_index(Testing::rand() % max_degree);
}
}
q_dof_handler.distribute_dofs(fe_q);
q_solution.reinit(q_dof_handler.n_dofs());
dgq_dof_handler.distribute_dofs(fe_dgq);
dgq_solution.reinit(dgq_dof_handler.n_dofs());
AffineConstraints<double> cm;
cm.close();
VectorTools::interpolate(mapping, q_dof_handler, function, q_solution);
VectorTools::interpolate(mapping, dgq_dof_handler, function, dgq_solution);
parallel::distributed::SolutionTransfer<dim, LinearAlgebra::distributed::Vector<double>, hp::DoFHandler<dim>> q_soltrans(q_dof_handler);
parallel::distributed::SolutionTransfer<dim, LinearAlgebra::distributed::Vector<double>, hp::DoFHandler<dim>> dgq_soltrans(dgq_dof_handler);
// test b): do some coarsening and
// refinement
counter = 0;
cell = tria.begin_active();
endc = tria.end();
for (; cell != endc; ++cell, ++counter)
{
if (counter > 120)
cell->set_coarsen_flag();
else if (Testing::rand() % 3 == 0)
cell->set_refine_flag();
else if (Testing::rand() % 3 == 3)
cell->set_coarsen_flag();
}
LinearAlgebra::distributed::Vector<double> q_old_solution = q_solution, dgq_old_solution = dgq_solution;
tria.prepare_coarsening_and_refinement();
q_soltrans.prepare_for_coarsening_and_refinement(q_old_solution);
dgq_soltrans.prepare_for_coarsening_and_refinement(dgq_old_solution);
tria.execute_coarsening_and_refinement();
counter = 0;
//{
// typename hp::DoFHandler<dim>::active_cell_iterator
// cell = q_dof_handler.begin_active(),
// celldg = dgq_dof_handler.begin_active(), endc = q_dof_handler.end();
// for (; cell != endc; ++cell, ++celldg, ++counter)
// {
// if (counter > 20 && counter < 90)
// cell->set_active_fe_index(0);
// else
// cell->set_active_fe_index(Testing::rand() % max_degree);
// if (counter > 20 && counter < 90)
// celldg->set_active_fe_index(0);
// else
// celldg->set_active_fe_index(Testing::rand() % max_degree);
// }
//}
q_dof_handler.distribute_dofs(fe_q);
dgq_dof_handler.distribute_dofs(fe_dgq);
q_solution.reinit(q_dof_handler.n_dofs());
dgq_solution.reinit(dgq_dof_handler.n_dofs());
q_soltrans.interpolate(q_solution);
dgq_soltrans.interpolate(dgq_solution);
// check correctness by comparing the values
// on points of QGauss of order 2.
MyFunction<dim> func;
{
double error = 0;
const hp::QCollection<dim> quad(QGauss<dim>(2));
hp::FEValues<dim> hp_fe_val(fe_q,
quad,
update_values | update_quadrature_points);
std::vector<double> vals(quad[0].size());
typename hp::DoFHandler<dim>::active_cell_iterator cell = q_dof_handler
.begin_active(),
endc =
q_dof_handler.end();
for (; cell != endc; ++cell)
{
hp_fe_val.reinit(cell, 0);
const FEValues<dim> &fe_val = hp_fe_val.get_present_fe_values();
fe_val.get_function_values(q_solution, vals);
for (unsigned int q = 0; q < fe_val.n_quadrature_points; ++q)
{
error +=
std::fabs(func.value(fe_val.quadrature_point(q), 0) - vals[q]);
}
}
deallog << "Error in interpolating hp FE_Q: " << error << std::endl;
}
{
double error = 0;
const hp::QCollection<dim> quad(QGauss<dim>(2));
hp::FEValues<dim> hp_fe_val(fe_dgq,
quad,
update_values | update_quadrature_points);
std::vector<double> vals(quad[0].size());
typename hp::DoFHandler<dim>::active_cell_iterator
celldg = dgq_dof_handler.begin_active(),
endc = dgq_dof_handler.end();
for (; celldg != endc; ++celldg)
{
hp_fe_val.reinit(celldg, 0);
const FEValues<dim> &fe_val = hp_fe_val.get_present_fe_values();
fe_val.get_function_values(dgq_solution, vals);
for (unsigned int q = 0; q < fe_val.n_quadrature_points; ++q)
{
error +=
std::fabs(func.value(fe_val.quadrature_point(q), 0) - vals[q]);
}
}
deallog << "Error in interpolating hp FE_DGQ: " << error << std::endl;
}
}
int
main(int argc, char *argv[])
{
#ifdef DEAL_II_WITH_MPI
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
#else
(void)argc;
(void)argv;
#endif
initlog();
// deallog << " 1D solution transfer" << std::endl;
// transfer<1>(deallog.get_file_stream());
deallog << " 2D solution transfer" << std::endl;
transfer<2>(deallog.get_file_stream());
deallog << " 3D solution transfer" << std::endl;
transfer<3>(deallog.get_file_stream());
}
