Hello Marc, I am trying to get the test working and unfortunately it is going a lot less smoothly than expected. Please find attached the test as solution_transfer_05.cc. Only h-refinement occurs. I have attached the before and after grids, which look OK. It is still that error
8134: -------------------------------------------------------- 8134: An error occurred in line <405> of file </home/ddong/Libraries/dealii/source/distributed/solution_transfer.cc> in function 8134: 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 = 3; VectorType = dealii::LinearAlgebra::distributed::Vector<double>; DoFHandlerType = dealii::hp::DoFHandler<3, 3>; typename dealii::parallel::distributed::Triangulation<dim, DoFHandlerType:: space_dimension>::cell_iterator = dealii::TriaIterator<dealii::CellAccessor<3, 3> >; typename dealii::parallel::distributed::Triangulation<dim, DoFHandlerType:: space_dimension>::CellStatus = dealii::Triangulation<3, 3>::CellStatus] 8134: The violated condition was: 8134: dof_handler->get_fe(fe_index).dofs_per_cell == it_dofvalues->size() 8134: Additional information: 8134: The transferred data was packed with a different number of dofs than the currently registered FE object assigned to the DoFHandler has. 8134: -------------------------------------------------------- Note that if I decide to refine the cells (uncomment line 150) instead of coarsening them (line 149), the solution transfer is successful. Playing around with various refinement/coarsening, it seems a bit arbitrary when it does work and when it won't. Am I calling anything out of order? Doug On Sunday, August 18, 2019 at 9:44:50 PM UTC-4, Marc Fehling wrote: > > Hi Doug, > > when dealing with distributed meshes, ownership of cells change and we may > not know which finite element lives on cells that the process got recently > assigned to. Thus, we need to transfer each cell's `active_fe_index`, which > we do automatically during coarsening and refinement. However, you set > `active_fe_indices` after refinement happened, which works in the serial > case, but no longer in the parallel one. Before executing refinement, you > need to set `future_fe_indices` that describe to which finite element your > cell will be assigned to, and you need to do that before refinement > happened! This should resolve both issues. > > Further, you initialize `LinearAlgebra::distrtibuted::Vector` objects > without any parallel distribution by using this constructor. > <https://www.dealii.org/current/doxygen/deal.II/classLinearAlgebra_1_1distributed_1_1Vector.html#a3be6c4ce529bb9b6c13eb831d0a86f55> > > Try using one a different one. > > Please see `tests/mpi/solution_transfer_04.cc > <https://github.com/dealii/dealii/blob/master/tests/mpi/solution_transfer_04.cc>` > > and `tests/mpi/p_coarsening_and_refinement.cc > <https://github.com/dealii/dealii/blob/master/tests/mpi/p_refinement_and_coarsening.cc>` > > for working examples (I guess we should provide one using an actual > `SolutionTransfer` object as well), that should hopefully be applicable to > your problem. This is a recently added feature: If you have any suggestions > for improvement or encounter more problems, feel free to message us! > > Marc > -- 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/8601d9b7-76f3-42c6-8972-d91ff3b640bc%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/mapping_q1.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_refinement.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/la_parallel_vector.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/distributed/solution_transfer.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/grid/grid_out.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 print_mesh_info(const dealii::parallel::distributed::Triangulation<dim> &triangulation, const std::string &filename)
{
if (dim == 2) {
std::ofstream out (filename);
dealii::GridOut grid_out;
grid_out.write_eps (triangulation, out);
std::cout << " written to " << filename << std::endl << std::endl;
}
}
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_dgq;
for (unsigned int deg = 1; deg <= max_degree; ++deg) {
fe_dgq.push_back(FE_DGQ<dim>(deg));
}
hp::DoFHandler<dim> dgq_dof_handler(tria);
MappingQGeneric<dim> mapping(1);
// refine a few cells
//++cell;
//++cell;
for (auto cell = tria.begin_active(); cell != tria.end(); ++cell) {
if (!(cell->is_locally_owned())) continue;
cell->set_refine_flag();
}
tria.prepare_coarsening_and_refinement();
tria.execute_coarsening_and_refinement();
print_mesh_info(tria, "before.eps");
// randomly assign FE orders
unsigned int counter = 0;
for (auto cell = dgq_dof_handler.begin_active(); cell != dgq_dof_handler.end(); ++cell, ++counter) {
if (!(cell->is_locally_owned())) continue;
if (counter < 15)
cell->set_active_fe_index(1);
else
cell->set_active_fe_index(Testing::rand() % max_degree);
}
dgq_dof_handler.distribute_dofs(fe_dgq);
IndexSet dgq_locally_owned_dofs = dgq_dof_handler.locally_owned_dofs();
IndexSet dgq_locally_relevant_dofs;
dealii::DoFTools::extract_locally_relevant_dofs(dgq_dof_handler, dgq_locally_relevant_dofs);
IndexSet dgq_ghost_dofs = dgq_locally_relevant_dofs;
dgq_ghost_dofs.subtract_set(dgq_locally_owned_dofs);
LinearAlgebra::distributed::Vector<double> dgq_solution;
dgq_solution.reinit(dgq_locally_owned_dofs, dgq_ghost_dofs, MPI_COMM_WORLD);
VectorTools::interpolate(mapping, dgq_dof_handler, function, dgq_solution);
parallel::distributed::SolutionTransfer<dim, LinearAlgebra::distributed::Vector<double>, hp::DoFHandler<dim>> dgq_soltrans(dgq_dof_handler);
// test b): do some coarsening and
// refinement
LinearAlgebra::distributed::Vector<double> dgq_old_solution = dgq_solution;
// h-refinement
counter = 0;
for (auto cell = tria.begin_active(); cell != tria.end(); ++cell, ++counter) {
if (!(cell->is_locally_owned())) continue;
//if (counter % 2 == 0) cell->set_coarsen_flag();
//else cell->set_coarsen_flag();
////else cell->set_refine_flag();
//continue;
if (counter > 120) {
cell->set_coarsen_flag();
//cell->set_refine_flag();
} else if (Testing::rand() % 3 == 0) {
cell->set_refine_flag();
} else if (Testing::rand() % 3 == 1) {
cell->set_coarsen_flag();
}
}
// p-refinement
counter = 0;
//for (auto cell = dgq_dof_handler.begin_active(); cell != dgq_dof_handler.end(); ++cell, ++counter) {
// if (!(cell->is_locally_owned())) continue;
// if (counter > 20 && counter < 90)
// cell->set_future_fe_index(0);
// else
// cell->set_future_fe_index(Testing::rand() % max_degree);
//}
tria.prepare_coarsening_and_refinement();
dgq_soltrans.prepare_for_coarsening_and_refinement(dgq_old_solution);
tria.execute_coarsening_and_refinement();
print_mesh_info(tria, "after.eps");
dgq_dof_handler.distribute_dofs(fe_dgq);
dgq_locally_owned_dofs = dgq_dof_handler.locally_owned_dofs();
dealii::DoFTools::extract_locally_relevant_dofs(dgq_dof_handler, dgq_locally_relevant_dofs);
dgq_ghost_dofs = dgq_locally_relevant_dofs;
dgq_ghost_dofs.subtract_set(dgq_locally_owned_dofs);
dgq_solution.reinit(dgq_locally_owned_dofs, dgq_ghost_dofs, MPI_COMM_WORLD);
dgq_soltrans.interpolate(dgq_solution);
// check correctness by comparing the values
// on points of QGauss of order 2.
{
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
cell = dgq_dof_handler.begin_active(),
endc = dgq_dof_handler.end();
for (; cell != endc; ++cell)
{
if (!(cell->is_locally_owned())) continue;
hp_fe_val.reinit(cell, 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(function.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();
// No 1D since parallel::distributed::Triangulation has not be implemented for 1D
// 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());
}
after.eps
Description: PostScript document
before.eps
Description: PostScript document
