Dear deal.ii community,
I am trying to interpolate a distributed solution vector onto a globally
refined distributed mesh. Actually I need to globally refine more than once
but I believe I can iterate this procedure (is there a more efficient way?
The finite element can be vector or scalar valued, and the solution vector
a block vector)
However, despite looking at tutorials and the documentation I get an error
(segmentation fault) without further explanation. I did not manage to track
down the error - probably I am using the interface not correctly. Anyone
has experience with this?
A minimal example is attached. (I used deal.ii 9.2.0 but it should work
with 9.1.1, too.)
Best regards,
Konrad
--
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/9978ebc5-4ff9-411f-b892-d8220a2146ec%40googlegroups.com.
##
# CMake script for the step-40 tutorial program:
##
# Set the name of the project and target:
SET(TARGET "interpolate")
# 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.2.0 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 AND NOT DEAL_II_PETSC_WITH_COMPLEX) OR
DEAL_II_WITH_TRILINOS) OR NOT DEAL_II_WITH_P4EST) # 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/base/function.h>
#include <deal.II/base/timer.h>
#include <deal.II/lac/generic_linear_algebra.h>
namespace LA {
using namespace dealii::LinearAlgebraTrilinos;
} // namespace LA
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/vector.h>
#include <deal.II/dofs/dof_accessor.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_values.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/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/index_set.h>
#include <deal.II/base/utilities.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/distributed/grid_refinement.h>
#include <deal.II/distributed/solution_transfer.h>
#include <deal.II/distributed/tria.h>
#include <fstream>
#include <iostream>
using namespace dealii;
template <int dim> class MyScalarFunction : public Function<dim> {
public:
MyScalarFunction() : Function<dim>() {}
virtual double value(const Point<dim> &p,
const unsigned int component = 0) const override;
virtual void value_list(const std::vector<Point<dim>> &points,
std::vector<double> &values,
const unsigned int component = 0) const override;
};
template <int dim>
double MyScalarFunction<dim>::value(const Point<dim> &p,
const unsigned int /*component*/) const {
double return_value = 0;
for (unsigned int d = 0; d < dim; ++d)
return_value += (p(d) - 0.5) * (p(d) - 0.5);
return return_value;
}
template <int dim>
void MyScalarFunction<dim>::value_list(
const std::vector<Point<dim>> &points, std::vector<double> &values,
const unsigned int /*component = 0*/) const {
Assert(points.size() == values.size(),
ExcDimensionMismatch(points.size(), values.size()));
for (unsigned int p = 0; p < points.size(); ++p) {
values[p] = value(points[p]);
} // end ++p
}
template <int dim> class RefineInterpolate {
public:
RefineInterpolate();
void run();
private:
void setup_system();
void project_on_distributed_mesh();
void refine_and_interpolate_on_distributed_mesh();
void output(const std::string &) const;
MPI_Comm mpi_communicator;
parallel::distributed::Triangulation<dim> triangulation;
FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
IndexSet locally_owned_dofs;
IndexSet locally_relevant_dofs;
AffineConstraints<double> constraints;
LA::MPI::Vector locally_relevant_solution;
ConditionalOStream pcout;
TimerOutput computing_timer;
};
template <int dim>
RefineInterpolate<dim>::RefineInterpolate()
: mpi_communicator(MPI_COMM_WORLD),
triangulation(mpi_communicator,
typename Triangulation<dim>::MeshSmoothing(
Triangulation<dim>::smoothing_on_refinement |
Triangulation<dim>::smoothing_on_coarsening)),
fe(1), dof_handler(triangulation),
pcout(std::cout,
(Utilities::MPI::this_mpi_process(mpi_communicator) == 0)),
computing_timer(mpi_communicator, pcout, TimerOutput::summary,
TimerOutput::wall_times) {}
template <int dim> void RefineInterpolate<dim>::setup_system() {
TimerOutput::Scope t(computing_timer, "setup");
dof_handler.distribute_dofs(fe);
locally_owned_dofs = dof_handler.locally_owned_dofs();
DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
locally_relevant_solution.reinit(locally_owned_dofs,
locally_relevant_dofs,
mpi_communicator);
constraints.clear();
constraints.reinit(locally_relevant_dofs);
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
constraints.close();
}
template <int dim> void RefineInterpolate<dim>::project_on_distributed_mesh() {
TimerOutput::Scope t(computing_timer, "project on distributed mesh");
MyScalarFunction<dim> my_function;
QGauss<dim> quad_rule(3);
LA::MPI::Vector locally_owned_solution;
locally_owned_solution.reinit(locally_owned_dofs,
mpi_communicator);
VectorTools::project(dof_handler, constraints, quad_rule, my_function,
locally_owned_solution);
locally_relevant_solution = locally_owned_solution;
}
template <int dim>
void RefineInterpolate<dim>::refine_and_interpolate_on_distributed_mesh() {
TimerOutput::Scope t(computing_timer, "refine and interpolate");
parallel::distributed::SolutionTransfer<dim, LA::MPI::Vector>
solution_transfer(dof_handler);
/*
* Mark everything for refinement.
*/
{
for (typename Triangulation<dim>::active_cell_iterator cell =
triangulation.begin_active();
cell != triangulation.end(); ++cell)
if (cell->is_locally_owned())
cell->set_refine_flag();
}
/*
* Prepare the triangulation for refinement.
*/
triangulation.prepare_coarsening_and_refinement();
/*
* Transfer locally_relevant_solution to vector that
* provides access to locally relevant DoFs.
*/
LA::MPI::Vector locally_relevant_solution_copy;
locally_relevant_solution_copy.reinit(
locally_owned_dofs, locally_relevant_dofs, mpi_communicator);
locally_relevant_solution_copy = locally_relevant_solution;
/*
* Prepare the refinement in the transfer object,
* locally_relevant_old_solution is the source.
*/
solution_transfer.prepare_for_coarsening_and_refinement(
locally_relevant_solution_copy);
/*
* Now actually refine the mesh
*/
triangulation.execute_coarsening_and_refinement();
/*
* Overwrite the information about local dofs after refinement
*/
locally_owned_dofs = dof_handler.locally_owned_dofs();
DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
/*
* New locally_owned_solution from new dofs.
*/
LA::MPI::Vector locally_owned_solution;
locally_owned_solution.reinit(locally_owned_dofs, mpi_communicator);
/*
* Now interpolate to new mesh.
*/
solution_transfer.interpolate(locally_owned_solution);
/*
* Take care of constraints.
*/
constraints.distribute(locally_owned_solution);
/*
* The new solution vector with locally relevant entries
* must be reinitialized.
*/
locally_relevant_solution.reinit(locally_owned_dofs, locally_relevant_dofs,
mpi_communicator);
locally_relevant_solution = locally_owned_solution;
}
template <int dim>
void RefineInterpolate<dim>::output(const std::string &add_str) const {
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(locally_relevant_solution, "u");
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();
std::string filename("solution_");
filename +=
add_str + "." +
Utilities::int_to_string(triangulation.locally_owned_subdomain(), 4);
std::ofstream output(filename + ".vtu");
data_out.write_vtu(output);
// pvtu-record for all local outputs
if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0) {
std::vector<std::string> local_filenames(
Utilities::MPI::n_mpi_processes(mpi_communicator), "solution_");
for (unsigned int i = 0;
i < Utilities::MPI::n_mpi_processes(mpi_communicator); ++i) {
local_filenames[i] +=
add_str + "." + Utilities::int_to_string(i, 4) + ".vtu";
}
std::string master_file = filename + ".pvtu";
std::ofstream master_output(master_file);
data_out.write_pvtu_record(master_output, local_filenames);
}
}
template <int dim> void RefineInterpolate<dim>::run() {
pcout << "Running with "
<< "Trilinos"
<< " on " << Utilities::MPI::n_mpi_processes(mpi_communicator)
<< " MPI rank(s)..." << std::endl;
GridGenerator::hyper_cube(triangulation);
triangulation.refine_global(3);
setup_system();
pcout << " Number of active cells: "
<< triangulation.n_global_active_cells() << std::endl
<< " Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl;
project_on_distributed_mesh();
{
TimerOutput::Scope t(computing_timer, "output");
output("before");
}
refine_and_interpolate_on_distributed_mesh();
{
TimerOutput::Scope t(computing_timer, "output");
output("after");
}
computing_timer.print_summary();
computing_timer.reset();
pcout << std::endl;
}
int main(int argc, char *argv[]) {
try {
using namespace dealii;
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
RefineInterpolate<2> refine_interpolate;
refine_interpolate.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;
}
