Hi there,
I've been playing with mesh refinement and solution transfer with the
Raviart-Thomas field, but found a very strange thing: in areas where the
mesh is coarsened, the solution becomes wrong.
In a very simple 2D test, I first created a uniform RT field, i.e. (1, 0),
with the help of `VectorTools::interpolate`, and then tried to refine the
left half of a square domain, while coarsening the rest. Of course, I also
tried to transfer the solution onto the new mesh. This function
`refine_mesh` is as follows.
```
template <int dim>
void
TestProblem<dim>::refine_mesh()
{
unsigned int cell_no = 0;
for (const auto &cell : dof_handler.active_cell_iterators())
{
Point<dim> center = cell->center();
if (center[0]<0.5)
{
cell->set_refine_flag();
}
else
{
cell->set_coarsen_flag();
}
}
SolutionTransfer<dim> solution_trans(dof_handler);
Vector<double> previous_solution;
previous_solution = solution;
triangulation.prepare_coarsening_and_refinement();
solution_trans.prepare_for_coarsening_and_refinement(previous_solution);
triangulation.execute_coarsening_and_refinement();
setup_dofs();
solution_trans.interpolate(previous_solution, solution);
constraints.distribute(solution);
}
```
As you can see in the attached plot, however, the transferred solution in
the coarsened area is wrong but OK in the refined part. I'm not sure what
goes wrong here in my code and would very much appreciate it if you could
give me a hint. A minimal but complete code to reproduce this plot is also
attached.
[image: visit0000.png]
Thanks.
Charlie
--
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/be6e661d-a3a9-4d74-b1a2-015d168ff9f1n%40googlegroups.com.
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/function.h>
#include <deal.II/base/utilities.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/block_sparsity_pattern.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/filtered_iterator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_raviart_thomas.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/numerics/derivative_approximation.h>
#include <deal.II/numerics/solution_transfer.h>
#include <fstream>
#include <iostream>
#include <string>
namespace Test
{
using namespace dealii;
template <int dim>
class InitialValues : public Function<dim>
{
public:
InitialValues()
: Function<dim>(dim)
{}
virtual void vector_value(const Point<dim> &p,
Vector<double> & value) const override;
};
template <int dim>
void InitialValues<dim>::vector_value(const Point<dim> &p,
Vector<double> & values) const
{
values(0) = 1.0;
values(1) = 0.0;
}
template <int dim>
class TestProblem
{
public:
TestProblem();
void run();
private:
void setup_dofs();
void output_results(int);
void refine_mesh();
Triangulation<dim> triangulation;
const FE_RaviartThomas<dim> fe;
DoFHandler<dim> dof_handler;
AffineConstraints<double> constraints;
Vector<double> solution;
Vector<double> rhs;
};
template <int dim>
TestProblem<dim>::TestProblem()
:
fe(0)
,
dof_handler(triangulation)
{}
template <int dim>
void TestProblem<dim>::setup_dofs()
{
dof_handler.distribute_dofs (fe);
std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl;
solution.reinit(dof_handler.n_dofs());
rhs.reinit(dof_handler.n_dofs());
constraints.clear();
DoFTools::make_hanging_node_constraints (dof_handler, constraints);
constraints.close();
}
template <int dim>
void TestProblem<dim>::output_results(int out_index)
{
std::vector<std::string> solution_names(dim, "u");
std::vector<DataComponentInterpretation::DataComponentInterpretation>
interpretation(dim,
DataComponentInterpretation::component_is_part_of_vector);
DataOut<dim> data_out;
data_out.add_data_vector(dof_handler,
solution,
solution_names,
interpretation);
data_out.build_patches();
const std::string filename = "sol-" + Utilities::int_to_string(out_index, 2) + ".vtu";
std::ofstream output(filename);
data_out.write_vtu(output);
}
template <int dim>
void
TestProblem<dim>::refine_mesh()
{
unsigned int cell_no = 0;
for (const auto &cell : dof_handler.active_cell_iterators())
{
Point<dim> center = cell->center();
if (center[0]<0.5)
{
cell->set_refine_flag();
}
else
{
cell->set_coarsen_flag();
}
}
SolutionTransfer<dim> solution_trans(dof_handler);
Vector<double> previous_solution;
previous_solution = solution;
triangulation.prepare_coarsening_and_refinement();
solution_trans.prepare_for_coarsening_and_refinement(previous_solution);
triangulation.execute_coarsening_and_refinement();
setup_dofs();
solution_trans.interpolate(previous_solution, solution);
constraints.distribute(solution);
}
template <int dim>
void TestProblem<dim>::run()
{
GridGenerator::hyper_cube(triangulation);
triangulation.refine_global(4);
setup_dofs();
{
const FEValuesExtractors::Vector velocities(0);
const FEValuesExtractors::Scalar pressure(dim);
VectorTools::interpolate(dof_handler,
InitialValues<dim>(),
solution,
ComponentMask());
constraints.distribute(solution);
}
int refine_step = 0;
while (refine_step <3)
{
refine_mesh();
output_results(refine_step);
refine_step++;
}
}
}
int main(int argc, char *argv[])
{
using namespace Test;
using namespace dealii;
const int dim = 2;
TestProblem<dim> test_problem;
test_problem.run();
return 0;
}
##
# CMake script for the step-32 tutorial program:
##
# Set the name of the project and target:
SET(TARGET "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
)
# Define the output that should be cleaned:
SET(CLEAN_UP_FILES *.vtu *.pvtu *.visit)
# 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_MPI OR NOT DEAL_II_WITH_P4EST OR NOT DEAL_II_WITH_TRILINOS)
# 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_MPI = ON
DEAL_II_WITH_P4EST = ON
DEAL_II_WITH_TRILINOS = ON
However, the deal.II library found at ${DEAL_II_PATH} was configured with these
options
DEAL_II_WITH_MPI = ${DEAL_II_WITH_MPI}
DEAL_II_WITH_P4EST = ${DEAL_II_WITH_P4EST}
DEAL_II_WITH_TRILINOS = ${DEAL_II_WITH_TRILINOS}
which conflict with the requirements."
)
ENDIF()
DEAL_II_INITIALIZE_CACHED_VARIABLES()
PROJECT(${TARGET})
DEAL_II_INVOKE_AUTOPILOT()
