Dear Deal.II Community,
I am trying to use parallel::distributed::Triangulation with
hp::DoFHandler. I run into the below error in
DoFTools::make_hanging_node_constraints..
I work on shape optimization using embedding domain discretization
techniques. In my problems, I have to repeatedly refine the mesh and
approximate the geometry. Then I assign different FESystem for the cells
within and outside the geometry.
I have implemented the code in serial with serial triangulation and
hp::DoFHandler and the program works.
Now I wish to solve larger problems so I have decided to parallelize the
code and hence I use parallel::distributed::Triangulation.
Here I run into a problem in system_setup, specifically in the
function DoFTools::make_hanging_node_constraints.
Hear is the error message.
--------------------------------------------------------
An error occurred in line <1119> of file
</calculate/temp/iwtm009/spack-stage-dealii-9.1.1-holxtax5bwynpjaccyvwkli2ybveslrl/spack-src/source/dofs/dof_tools_constraints.cc>
in function
void dealii::DoFTools::internal::make_hp_hanging_node_constraints(const
DoFHandlerType&, dealii::AffineConstraints<number>&) [with DoFHandlerType =
dealii::hp::DoFHandler<2, 2>; number = double]
The violated condition was:
cell->face(face)->child(c)->n_active_fe_indices() == 1
Additional information:
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.
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.
Stacktrace:
-----------
#0
/opt/dealii-9.1.1/opt/spack/linux-ubuntu18.04-haswell/gcc-7.4.0/dealii-9.1.1-holxtax5bwynpjaccyvwkli2ybveslrl/lib/libdeal_II.g.so.9.1.1:
void
dealii::DoFTools::internal::make_hp_hanging_node_constraints<dealii::hp::DoFHandler<2,
2>, double>(dealii::hp::DoFHandler<2, 2> const&,
dealii::AffineConstraints<double>&)
#1
/opt/dealii-9.1.1/opt/spack/linux-ubuntu18.04-haswell/gcc-7.4.0/dealii-9.1.1-holxtax5bwynpjaccyvwkli2ybveslrl/lib/libdeal_II.g.so.9.1.1:
void
dealii::DoFTools::make_hanging_node_constraints<dealii::hp::DoFHandler<2,
2>, double>(dealii::hp::DoFHandler<2, 2> const&,
dealii::AffineConstraints<double>&)
#2 hp-dist-test: TestProblem<2>::setup_system()
#3 hp-dist-test: TestProblem<2>::run()
#4 hp-dist-test: main
--------------------------------------------------------
This is my setup_system function.
*template<int dim>*
*void TestProblem<dim>::setup_system() {*
* dof_handler.initialize(triangulation, fe_collection);*
* for (auto cell : dof_handler.active_cell_iterators()) {*
* if (!cell->is_locally_owned())*
* continue;*
* if (cell->user_flag_set())*
* cell->set_active_fe_index(0);*
* else*
* cell->set_active_fe_index(1);*
* }*
* dof_handler.distribute_dofs(fe_collection);*
* locally_owned_dofs = dof_handler.locally_owned_dofs();*
* DoFTools::extract_locally_relevant_dofs(dof_handler,
locally_relevant_dofs);*
* hanging_node_constraints.clear();*
* hanging_node_constraints.reinit(locally_relevant_dofs);*
* DoFTools::make_hanging_node_constraints(dof_handler,*
* hanging_node_constraints);*
* hanging_node_constraints.close();*
*}*
I did some debugging and found that the Assert fail
cell->face(face)->child(c)->n_active_fe_indices()
== 1 is failing on the *ghost cells*. To avoid this error I tried to
set_active_fe_index for the ghost cells, but this is not allowed.
Please let me know where I am going wrong.
I have attached a MWE. You can reproduce the error with mpirun for np >= 2
cores (the code will work for np=1, because no ghost cells).
Thank you,
Chaitanya Dev.
--
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].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/626cb4ca-8971-4511-a4d1-3ed255557388%40googlegroups.com.
##
# CMake script
##
# Set the name of the project and target:
SET(TARGET "hp-dist-test")
# 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 2.8.12)
FIND_PACKAGE(deal.II 9.1.1 QUIET
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()
#
# Are all dependencies fulfilled?
#
IF(NOT (DEAL_II_WITH_PETSC OR DEAL_II_WITH_TRILINOS) OR NOT DEAL_II_WITH_P4EST
OR DEAL_II_PETSC_WITH_COMPLEX) # keep in one line
MESSAGE(FATAL_ERROR "
Error! This tutorial requires a deal.II library that was configured with the
following options:
DEAL_II_WITH_PETSC = ON
DEAL_II_PETSC_WITH_COMPLEX = OFF
DEAL_II_WITH_P4EST = ON
or
DEAL_II_WITH_TRILINOS = ON
DEAL_II_WITH_P4EST = ON
However, the deal.II library found at ${DEAL_II_PATH} was configured with these
options
DEAL_II_WITH_PETSC = ${DEAL_II_WITH_PETSC}
DEAL_II_PETSC_WITH_COMPLEX = ${DEAL_II_PETSC_WITH_COMPLEX}
DEAL_II_WITH_P4EST = ${DEAL_II_WITH_P4EST}
DEAL_II_WITH_TRILINOS = ${DEAL_II_WITH_TRILINOS}
which conflict with the requirements.
One or both of the aforementioned combinations of prerequisites are not met by
your installation, but at least one is required for this tutorial step."
)
ENDIF()
DEAL_II_INITIALIZE_CACHED_VARIABLES()
SET(CLEAN_UP_FILES *.log *.gmv *.gnuplot *.gpl *.eps *.pov *.vtk *.ucd *.d2
*.vtu *.pvtu)
PROJECT(${TARGET})
DEAL_II_INVOKE_AUTOPILOT()
#include <deal.II/lac/generic_linear_algebra.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/distributed/grid_refinement.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.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/fe/fe_system.h>
#include <deal.II/fe/fe_nothing.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/index_set.h>
using namespace dealii;
template<int dim>
class TestProblem {
public:
TestProblem();
~TestProblem();
void run();
private:
void make_grid();
void setup_system();
void output_results() const;
MPI_Comm mpi_communicator;
const unsigned int n_mpi_processes;
const unsigned int this_mpi_process;
ConditionalOStream pcout;
parallel::distributed::Triangulation<dim> triangulation;
hp::DoFHandler<dim> dof_handler;
FESystem<dim> void_fe;
FESystem<dim> solid_fe;
IndexSet locally_owned_dofs;
IndexSet locally_relevant_dofs;
hp::FECollection<dim> fe_collection;
AffineConstraints<double> hanging_node_constraints;
};
template<int dim>
TestProblem<dim>::TestProblem() :
mpi_communicator(MPI_COMM_WORLD), n_mpi_processes(
Utilities::MPI::n_mpi_processes(mpi_communicator)), this_mpi_process(
Utilities::MPI::this_mpi_process(mpi_communicator)), pcout(
std::cout, (this_mpi_process == 0)), triangulation(
mpi_communicator,
typename Triangulation<dim>::MeshSmoothing(
Triangulation<dim>::smoothing_on_refinement
| Triangulation<dim>::smoothing_on_coarsening)), dof_handler(
triangulation), void_fe(FE_Nothing<dim>(), dim), solid_fe(
FE_Q<dim>(1), dim) {
fe_collection.push_back(void_fe);
fe_collection.push_back(solid_fe);
}
template<int dim>
TestProblem<dim>::~TestProblem() {
dof_handler.clear();
}
template<int dim>
void TestProblem<dim>::make_grid() {
// Generate triangulation
GridGenerator::hyper_rectangle(triangulation,
(dim == 3 ? Point<dim>(0.0, 0.0, 0.0) : Point<dim>(0.0, 0.0)),
(dim == 3 ? Point<dim>(1.0, 1.0, 1.0) : Point<dim>(1.0, 1.0)),
true);
triangulation.refine_global(3);
// Refine all the cells with y > 0.5 twice
unsigned int n_ref = 0;
while (n_ref < 2) {
for (auto cell : triangulation.active_cell_iterators()) {
if (!cell->is_locally_owned())
continue;
if (cell->center()(1) < 0.5)
cell->set_refine_flag();
}
triangulation.prepare_coarsening_and_refinement();
triangulation.execute_coarsening_and_refinement();
++n_ref;
}
// Set user flag to assign differnt fe_index.
for (auto cell : triangulation.active_cell_iterators()) {
if (!cell->is_locally_owned())
continue;
if (cell->center()(0) < 0.5)
cell->set_user_flag();
}
}
template<int dim>
void TestProblem<dim>::setup_system() {
dof_handler.initialize(triangulation, fe_collection);
// set_fe_index
for (auto cell : dof_handler.active_cell_iterators()) {
// Do nothing is not locally owned.
if (!cell->is_locally_owned())
continue;
if (cell->user_flag_set())
cell->set_active_fe_index(0);
else
cell->set_active_fe_index(1);
}
// setup_dof
dof_handler.distribute_dofs(fe_collection);
locally_owned_dofs = dof_handler.locally_owned_dofs();
DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
// Set constraints
hanging_node_constraints.clear();
hanging_node_constraints.reinit(locally_relevant_dofs);
// This below block of code is to show the error.
// The Assert fails in make_hp_hanging_node_constraints() in dof_tools_constraints.cc
// Here I want to show this assert fails only for he ghost cells.
{
for (auto cell : dof_handler.active_cell_iterators()) {
if (cell->is_artificial())
continue;
for (unsigned int face = 0;
face < GeometryInfo<dim>::faces_per_cell; ++face)
if (cell->face(face)->has_children()) {
if (cell->get_fe().dofs_per_face == 0)
continue;
for (unsigned int c = 0; c < cell->face(face)->n_children();
++c) {
if (cell->face(face)->child(c)->n_active_fe_indices()
!= 1)
std::cout << __FILE__ << " " << __LINE__ << " P:"
<< this_mpi_process
<< " active fe index of face: "
<< cell->face(face)->child(c)->n_active_fe_indices()
<< " is ghost: " << cell->is_ghost()
<< std::endl;
}
}
}
}
DoFTools::make_hanging_node_constraints(dof_handler,
hanging_node_constraints);
hanging_node_constraints.close();
}
template<int dim>
void TestProblem<dim>::output_results() const {
DataOut<dim, hp::DoFHandler<dim> > data_out;
data_out.attach_dof_handler(dof_handler);
Vector<float> subdomain(triangulation.n_active_cells());
for (unsigned int i = 0; i < subdomain.size(); ++i)
subdomain(i) = triangulation.locally_owned_subdomain();
data_out.add_data_vector(subdomain, "subdomain");
data_out.build_patches();
Vector<float> flags(triangulation.n_active_cells());
unsigned int idx = 0;
for (auto cell : triangulation.active_cell_iterators()) {
if (cell->user_flag_set())
flags[idx] = 1.0;
++idx;
}
data_out.add_data_vector(flags, "flags");
data_out.build_patches();
const std::string filename = ("solution."
+ Utilities::int_to_string(triangulation.locally_owned_subdomain(),
4));
std::ofstream output((filename + ".vtu"));
data_out.write_vtu(output);
if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0) {
std::vector<std::string> filenames;
for (unsigned int i = 0;
i < Utilities::MPI::n_mpi_processes(mpi_communicator); ++i)
filenames.push_back(
"solution." + Utilities::int_to_string(i, 4) + ".vtu");
std::ofstream master_output("solution.pvtu");
data_out.write_pvtu_record(master_output, filenames);
}
}
template<int dim>
void TestProblem<dim>::run() {
make_grid();
// Output only to visualize the locally_owned_subdomain and the user_set_flags
output_results();
setup_system();
}
int main(int argc, char *argv[]) {
try {
using namespace dealii;
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
TestProblem<2> test_problem;
test_problem.run();
} 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;
}
