I wanted to write some functions to convert my solution from using
FE_Q-elements to FE_Bernstein-elements (and back). For that I wrote two
functions, convert_feQ_to_feB() and convert_feB_to_feQ() (as listed in the
attachment). When converting from FE_Q-elements to FE_Bernstein-elements, I
get the same values for input and output, but on the way back I get wrong
values (currently I am just using a vector filled with the value 1), but I
can not find the reason why the conversion back is failing. Did I forget
something in the corresponding function?
Thanks!
--
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/e288a022-ff5b-4876-8ca3-9807605ef2a3%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/fe/fe.h>
#include <deal.II/fe/fe_bernstein.h>
#include <deal.II/fe/mapping_q.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/tria.h>
#include <deal.II/lac/vector.h>
#include <deal.II/numerics/fe_field_function.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/vector_tools.templates.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/numerics/solution_transfer.h>
#include <deal.II/numerics/error_estimator.h>
using namespace dealii;
template <int dim>
void convert_feB_to_feQ(const DoFHandler<dim> &dhB,
const DoFHandler<dim> &dhq,
const Vector<double> &in,
Vector<double> &out){
AssertDimension(in.size(), dhB.n_dofs());
const FiniteElement<dim> &feB = dhB.get_fe();
const FiniteElement<dim> &feq = dhq.get_fe();
const ComponentMask fe_mask(feB.get_nonzero_components(0).size(), true);
AssertDimension(fe_mask.size(), feB.get_nonzero_components(0).size());
std::vector<unsigned int> fe_to_real(fe_mask.size(),
numbers::invalid_unsigned_int);
unsigned int size = 0;
for (unsigned int i = 0; i < fe_mask.size(); ++i)
{
if (fe_mask[i])
fe_to_real[i] = size++;
}
const ComponentMask maskq(dim, true);
FullMatrix<double> transfer(feq.dofs_per_cell, feB.dofs_per_cell);
FullMatrix<double> local_transfer(feB.dofs_per_cell);
const std::vector<Point<dim>> &points = feq.get_unit_support_points();
std::vector<unsigned int> fe_to_feq(feB.dofs_per_cell,
numbers::invalid_unsigned_int);
unsigned int index = 0;
for (unsigned int i = 0; i < feB.dofs_per_cell; ++i)
if (fe_mask[feB.system_to_component_index(i).first])
fe_to_feq[i] = index++;
Assert(index == feB.dofs_per_cell, ExcNotImplemented());
for (unsigned int j = 0; j < feB.dofs_per_cell; ++j)
{
const unsigned int comp_j = feB.system_to_component_index(j).first;
if (fe_mask[comp_j])
for (unsigned int i = 0; i < points.size(); ++i)
{
if (fe_to_real[comp_j] ==
feB.system_to_component_index(i).first)
local_transfer(i, fe_to_feq[j]) =
feq.shape_value(j, points[i]);
}
}
local_transfer.invert(local_transfer);
for (unsigned int i = 0; i < feq.dofs_per_cell; ++i)
if (fe_to_feq[i] != numbers::invalid_unsigned_int)
for (unsigned int j = 0; j < feB.dofs_per_cell; ++j)
transfer(i, j) = local_transfer(fe_to_feq[i], j);
VectorTools::interpolate(dhB, dhq, transfer, in, out);
}
template <int dim>
void convert_feQ_to_feB(const DoFHandler<dim> &dhq,
const DoFHandler<dim> &dhb,
const Vector<double> &in,
Vector<double> &out){
AssertDimension(in.size(), dhq.n_dofs());
const FiniteElement<dim> &feq = dhq.get_fe();
const FiniteElement<dim> &feB = dhb.get_fe();
const ComponentMask fe_mask(feq.get_nonzero_components(0).size(), true);
AssertDimension(fe_mask.size(), feq.get_nonzero_components(0).size());
std::vector<unsigned int> fe_to_real(fe_mask.size(),
numbers::invalid_unsigned_int);
unsigned int size = 0;
for (unsigned int i = 0; i < fe_mask.size(); ++i)
{
if (fe_mask[i])
fe_to_real[i] = size++;
}
const ComponentMask maskq(dim, true);
FullMatrix<double> transfer(feq.dofs_per_cell, feB.dofs_per_cell);
FullMatrix<double> local_transfer(feB.dofs_per_cell);
const std::vector<Point<dim>> &points = feq.get_unit_support_points();
std::vector<unsigned int> fe_to_feq(feq.dofs_per_cell,
numbers::invalid_unsigned_int);
unsigned int index = 0;
for (unsigned int i = 0; i < feq.dofs_per_cell; ++i)
if (fe_mask[feq.system_to_component_index(i).first])
fe_to_feq[i] = index++;
Assert(index == feq.dofs_per_cell, ExcNotImplemented());
for (unsigned int j = 0; j < feq.dofs_per_cell; ++j)
{
const unsigned int comp_j = feq.system_to_component_index(j).first;
if (fe_mask[comp_j])
for (unsigned int i = 0; i < points.size(); ++i)
{
if (fe_to_real[comp_j] ==
feq.system_to_component_index(i).first)
local_transfer(i, fe_to_feq[j]) =
feq.shape_value(j, points[i]);
}
}
local_transfer.invert(local_transfer);
for (unsigned int i = 0; i < feq.dofs_per_cell; ++i)
if (fe_to_feq[i] != numbers::invalid_unsigned_int)
for (unsigned int j = 0; j < feB.dofs_per_cell; ++j)
transfer(i, j) = local_transfer(fe_to_feq[i], j);
VectorTools::interpolate(dhq, dhb, transfer, in, out);
}
template <int dim, typename T>
void test_interpolation(const FiniteElement<dim> &fe_1,
const FiniteElement<dim> &fe_2,
const T& f,
const size_t order_mapping,
bool distort_mesh = false,
bool print_function_values = false){
std::cout << "dim " << dim << " " << fe_1.get_name() << std::endl;
std::cout << "dim " << dim << " " << fe_2.get_name() << std::endl;
Triangulation<dim> triangulation;
GridGenerator::hyper_cube(triangulation, -0.3, 0.7);
triangulation.refine_global(dim == 2 ? 1 : 1);
MappingQ<dim> mapping(order_mapping);
DoFHandler<dim> dhq(triangulation), dhb(triangulation);
dhq.distribute_dofs(fe_1);
dhb.distribute_dofs(fe_2);
Vector<double> iq(dhq.n_dofs()), ib(dhb.n_dofs()), ib2(dhb.n_dofs());
VectorTools::interpolate(dhq, f, iq);
VectorTools::interpolate(dhb, f, ib);
convert_feQ_to_feB<dim>(dhq, dhb, iq, ib2);
convert_feB_to_feQ<dim>(dhb, dhq, ib, iq);
std::cout << ib << '\n';
std::cout << ib2 << '\n';
std::cout << iq << '\n';
std::cout << iq.l2_norm() << '\n';
}
int
main()
{
test_interpolation<2>(FE_Q<2>(2), FE_Bernstein<2>(2), ConstantFunction<2>(1.0, 1), 2);
}
