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());
}


Attachment: after.eps
Description: PostScript document

Attachment: before.eps
Description: PostScript document

Reply via email to