Dear deal-ii users and developers,
I am trying to use the FETools::extrapolate method on a distributed
triangulation that is adaptively refined. However, I get an exception
inside this function every time my triangulation contains hanging nodes
at ghost interfaces. I have attached a minimal working example that
illustrates this on deal.II version 9.6.1. It should be compiled in
debug mode as for usual dealii example codes, and then run with "mpirun
-np 2 main" with two mpi ranks to generate the exception. The macro
NO_GHOST_HANGING_NODES can be commented out to check that
FETools::extrapolate works correctly when there are no hanging nodes at
ghost interfaces.
Before opening a bug, I would like to check with you that I am not doing
something wrong in this test program. I was careful to set up my
distributed vectors as large as possible (i.e. with all relevant dofs)
to see if the problem was coming from there, and of course I updated the
ghost values of the "coarse" vector before calling FETools::extrapolate.
If I am not mistaken, the exception is thrown when an internal vector
with all relevant dofs is compressed line 1455 of
fe_tools_extrapolate.templates.h, but of course this internal vector
depends on the input vector I give to extrapolate so the problem could
very much come from my code. I am ready to provide any further
information that may be useful! Thanks in advance.
Best regards,
Guilhem
--
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 visit
https://groups.google.com/d/msgid/dealii/31128be5-20ae-4207-bb06-a58c2832d8fa%40gmail.com.
##
# CMake script for the step-1 tutorial program:
##
# Set the name of the project and target:
set(TARGET "main")
# Declare all source files the target consists of. Here, this is only
# the one step-X.cc file, but as you expand your project you may wish
# to add other source files as well. If your project becomes much larger,
# you may want to either replace the following statement by something like
# file(GLOB_RECURSE TARGET_SRC "source/*.cc")
# file(GLOB_RECURSE TARGET_INC "include/*.h")
# set(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
# or switch altogether to the large project CMakeLists.txt file discussed
# in the "CMake in user projects" page accessible from the "User info"
# page of the documentation.
set(TARGET_SRC
${TARGET}.cc
)
# Usually, you will not need to modify anything beyond this point...
cmake_minimum_required(VERSION 3.13.4)
find_package(deal.II 9.6.1
HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
)
if(NOT ${deal.II_FOUND})
message(FATAL_ERROR "\n"
"*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
"You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to
cmake\n"
"or set an environment variable \"DEAL_II_DIR\" that contains this path."
)
endif()
deal_ii_initialize_cached_variables()
project(${TARGET})
deal_ii_invoke_autopilot()
#include <deal.II/base/function_lib.h>
#include <deal.II/base/multithread_info.h>
#include <deal.II/base/point.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_tools.h>
#include <deal.II/grid/filtered_iterator.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/lac/la_parallel_vector.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <filesystem>
// #define NO_GHOST_HANGING_NODES
using namespace dealii;
int main(int argc, char *argv[]) {
try {
Utilities::MPI::MPI_InitFinalize mpi_obj(argc, argv, 1);
/**
* Generate a triangulation with ghost hanging nodes when #(MPI ranks)=2
*/
parallel::distributed::Triangulation<2> tria(
MPI_COMM_WORLD, Triangulation<2>::limit_level_difference_at_vertices,
parallel::distributed::Triangulation<2>::construct_multigrid_hierarchy);
GridGenerator::hyper_rectangle(tria, Point<2>(-1.,-1.), Point<2>(1.,1.));
tria.refine_global(2);
for(auto &cell : tria.active_cell_iterators() | IteratorFilters::LocallyOwnedCell()) {
auto center = cell->barycenter();
#ifdef NO_GHOST_HANGING_NODES
if((center-Point<2>(0.,0.)).norm()<0.5)
#else
if((center-Point<2>(0.,0.)).norm()<0.5 && center[0]*center[1]>0)
#endif
cell->set_refine_flag();
}
tria.prepare_coarsening_and_refinement();
tria.execute_coarsening_and_refinement();
/**
* DoFs, constraints and vecs on two FE_Q spaces with different degrees
*/
FE_Q<2> fe[2] = {FE_Q<2>(2), FE_Q<2>(4)};
DoFHandler<2> dof_handlers[2];
AffineConstraints<double> constraints[2];
LinearAlgebra::distributed::Vector<double> vecs[2];
DataOut<2> data_out[2];
for(unsigned int i=0; i<2; ++i) {
dof_handlers[i].reinit(tria);
dof_handlers[i].distribute_dofs(fe[i]);
data_out[i].attach_dof_handler(dof_handlers[i]);
auto locally_owned_dofs = dof_handlers[i].locally_owned_dofs();
auto locally_relevant_dofs = DoFTools::extract_locally_relevant_dofs(dof_handlers[i]);
constraints[i].clear();
constraints[i].reinit(locally_owned_dofs, locally_relevant_dofs);
DoFTools::make_hanging_node_constraints(dof_handlers[i], constraints[i]);
constraints[i].close();
vecs[i].reinit(locally_owned_dofs, locally_relevant_dofs, MPI_COMM_WORLD);
}
/**
* Try extrapolating a simple cosine function from the first DoFHandler to the second
*/
std::filesystem::create_directory("results");
DataOutBase::VtkFlags flags;
flags.write_higher_order_cells = true;
for(unsigned int i=0; i<2; ++i) {
if(i==0)
VectorTools::interpolate(dof_handlers[0], Functions::CosineFunction<2>(), vecs[0]);
else {
vecs[0].update_ghost_values();
FETools::extrapolate(dof_handlers[0], vecs[0], dof_handlers[1], constraints[1], vecs[1]);
vecs[1].update_ghost_values();
}
data_out[i].add_data_vector(vecs[i], "u");
data_out[i].build_patches(fe[i].tensor_degree());
data_out[i].set_flags(flags);
data_out[i].write_vtu_with_pvtu_record(
"results/", "sol", i, MPI_COMM_WORLD, 1);
}
}
catch(std::exception &exc) {
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch(...) {
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
