I tried to implement this idea with the following code:
TrilinosWrappers::MPI::Vector distributed_q_solution(locally_owned_dofs,
mpi_communicator);
TrilinosWrappers::MPI::Vector distributed_old_q_solution(
locally_owned_dofs, mpi_communicator);
TrilinosWrappers::MPI::Vector distributed_really_old_q_solution(
locally_owned_dofs, mpi_communicator);
TrilinosWrappers::MPI::Vector distributed_local_present_solution(
locally_owned_dofs, mpi_communicator);
TrilinosWrappers::MPI::Vector distributed_local_old_solution(
locally_owned_dofs, mpi_communicator);
TrilinosWrappers::MPI::Vector distributed_local_really_old_solution(
locally_owned_dofs, mpi_communicator);
distributed_local_present_solution = present_solution;
distributed_local_old_solution = old_solution;
distributed_local_really_old_solution = really_old_solution;
//Transfer data from FE_Bernstein to FE_Q
const ComponentMask fe_mask(fe.get_nonzero_components(0).size(), true);
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 FESystem<dim, dim> *fe_system =
dynamic_cast<const FESystem<dim, dim> *>(&fe);
Assert(fe_system, ExcNotImplemented());
unsigned int degree = numbers::invalid_unsigned_int;
// Get information about the blocks
for (unsigned int i = 0; i < fe_mask.size(); ++i)
if (fe_mask[i])
{
const unsigned int base_i =
fe_system->component_to_base_index(i).first;
Assert(degree == numbers::invalid_unsigned_int ||
degree == fe_system->base_element(base_i).degree,
ExcNotImplemented());
Assert(fe_system->base_element(base_i).is_primitive(),
ExcNotImplemented());
degree = fe_system->base_element(base_i).degree;
}
FESystem<dim> feq(FE_Q<dim>(fe_degree_to_use), n_components);
DoFHandler<dim> dhq(dof_handler.get_triangulation());
dhq.distribute_dofs(feq);
FullMatrix<double> transfer_forward(feq.dofs_per_cell, fe.dofs_per_cell
);
FullMatrix<double> local_transfer_forward(fe.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++;
// If index is not the same as feq.dofs_per_cell, we won't
// know how to invert the resulting matrix. Bail out.
Assert(index == fe.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] ==
fe.system_to_component_index(i).first)
local_transfer_forward(i, fe_to_feq[j]) =
feq.shape_value(j, points[i]);
}
}
// Now we construct the rectangular interpolation matrix. This
// one is filled only with the information from the components
// of the displacement. The rest is set to zero.
local_transfer_forward.invert(local_transfer_forward);
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 < fe.dofs_per_cell; ++j)
transfer_forward(i, j) = local_transfer_forward(fe_to_feq[i
], j);
// The interpolation matrix is then passed to the
// VectorTools::interpolate() function to generate the correct
// interpolation.
pcout << "Data transferred to feq, interpolating\n";
VectorTools::interpolate(dof_handler, dhq, transfer_forward,
distributed_local_present_solution, distributed_q_solution);
VectorTools::interpolate(dof_handler, dhq, transfer_forward,
distributed_local_old_solution, distributed_old_q_solution);
VectorTools::interpolate(dof_handler, dhq, transfer_forward,
distributed_local_really_old_solution, distributed_really_old_q_solution);
pcout << "Vectortools::interpolate done\n";
parallel::distributed::SolutionTransfer<dim, TrilinosWrappers::MPI::
Vector, DoFHandler<dim>> solution_transfer(dhq);
parallel::distributed::SolutionTransfer<dim, TrilinosWrappers::MPI::
Vector, DoFHandler<dim>> old_solution_transfer(dhq);
parallel::distributed::SolutionTransfer<dim, TrilinosWrappers::MPI::
Vector, DoFHandler<dim>> really_old_solution_transfer(dhq);
pcout << "SolutionTransfer created\n";
solution_transfer.prepare_for_coarsening_and_refinement(
distributed_q_solution);
old_solution_transfer.prepare_for_coarsening_and_refinement(
distributed_old_q_solution);
really_old_solution_transfer.prepare_for_coarsening_and_refinement(
distributed_really_old_q_solution);
triangulation.execute_coarsening_and_refinement();
setup_system(true);
print_status_update(pcout, std::string("Interpolation\n"), true);
TrilinosWrappers::MPI::Vector distributed_solution(locally_owned_dofs,
mpi_communicator);
TrilinosWrappers::MPI::Vector distributed_old_solution(
locally_owned_dofs, mpi_communicator);
TrilinosWrappers::MPI::Vector distributed_really_old_solution(
locally_owned_dofs, mpi_communicator);
print_status_update(pcout, std::string("Doing interpolation\n"), true);
solution_transfer.interpolate(distributed_q_solution);
old_solution_transfer.interpolate(distributed_old_q_solution);
really_old_solution_transfer.interpolate(
distributed_really_old_q_solution);
pcout << "Data interpolated, transferring back\n";
//Transferring back
size = 0;
for (unsigned int i = 0; i < fe_mask.size(); ++i)
{
if (fe_mask[i])
fe_to_real[i] = size++;
}
// Get information about the blocks
for (unsigned int i = 0; i < fe_mask.size(); ++i)
if (fe_mask[i])
{
const unsigned int base_i =
fe_system->component_to_base_index(i).first;
Assert(degree == numbers::invalid_unsigned_int ||
degree == fe_system->base_element(base_i).degree,
ExcNotImplemented());
Assert(fe_system->base_element(base_i).is_primitive(),
ExcNotImplemented());
degree = fe_system->base_element(base_i).degree;
}
FullMatrix<double> transfer_backward(fe.dofs_per_cell, feq.dofs_per_cell
);
FullMatrix<double> local_transfer_backward(feq.dofs_per_cell);
index = 0;
for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
if (fe_mask[fe.system_to_component_index(i).first])
fe_to_feq[i] = index++;
// If index is not the same as feq.dofs_per_cell, we won't
// know how to invert the resulting matrix. Bail out.
Assert(index == feq.dofs_per_cell, ExcNotImplemented());
for (unsigned int j = 0; j < fe.dofs_per_cell; ++j)
{
const unsigned int comp_j = fe.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_backward(i, fe_to_feq[j]) =
fe.shape_value(j, points[i]);
}
}
// Now we construct the rectangular interpolation matrix. This
// one is filled only with the information from the components
// of the displacement. The rest is set to zero.
local_transfer_backward.invert(local_transfer_backward);
for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
if (fe_to_feq[i] != numbers::invalid_unsigned_int)
for (unsigned int j = 0; j < feq.dofs_per_cell; ++j)
transfer_backward(i, j) = local_transfer_backward(fe_to_feq[
i], j);
// The interpolation matrix is then passed to the
// VectorTools::interpolate() function to generate the correct
// interpolation.
VectorTools::interpolate(dhq, dof_handler, transfer_backward,
distributed_q_solution, distributed_solution);
VectorTools::interpolate(dhq, dof_handler, transfer_backward,
distributed_old_q_solution, distributed_old_solution);
VectorTools::interpolate(dhq, dof_handler, transfer_backward,
distributed_really_old_q_solution, distributed_really_old_solution);
but it fails during solution_transfer.interpolate() with the error:
An error occurred in line <121> of file </opt/dealii/include/deal.II/dofs/
dof_levels.h> in function
const global_dof_index* dealii::internal::DoFHandlerImplementation::
DoFLevel<<anonymous> >::get_cell_cache_start(unsigned int, unsigned int)
const [with int dim = 2; dealii::types::global_dof_index = unsigned int]
The violated condition was:
obj_index * dofs_per_cell + dofs_per_cell <= cell_dof_indices_cache.size
()
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.
I can try to create an MWE, if necessary, but are there any obvious bugs in
the code?
Thanks!
Am Dienstag, 21. Mai 2019 14:08:32 UTC+2 schrieb Maxi Miller:
>
> According to
> https://www.dealii.org/developer/doxygen/deal.II/classFE__Bernstein.html#ad881d4e04f5e699e4a8bda5b4e9f5265
>
> that is impossible. Would an alternative solution be to interpolate
> everything onto a grid using FE_Q-elements, transfer the result, and
> interpolate the results back?
>
> Am Dienstag, 21. Mai 2019 03:41:58 UTC+2 schrieb Wolfgang Bangerth:
>>
>> On 5/20/19 5:18 AM, 'Maxi Miller' via deal.II User Group wrote:
>> > Trying to transfer solutions from an old grid to a new grid with
>> > parallel::distributed::SolutionTransfer<dim,
>> TrilinosWrappers::MPI::Vector,
>> > DoFHandler<dim>> solution_transfer() and
>> > solution_transfer.prepare_for_coarsening_and_refinement() fails with
>> > |
>> > Anerror occurred inline <71>of file
>> > <~/Downloads/git-files/dealii/source/fe/fe_bernstein.cc>infunction
>> >
>> constdealii::FullMatrix<double>&dealii::FE_Bernstein<dim,spacedim>::get_restriction_matrix(unsignedint,constdealii::RefinementCase<dim>&)const[withintdim
>>
>>
>> > =2;intspacedim =2]
>> > Theviolated condition was:
>> > false
>> > Additionalinformation:
>> > Youare trying to access the matrices that describe how to restrict a
>> finite
>> > element functionfromthe children of one cell to the finite element
>> space
>> > definedon their parent (i.e.,the
>> > 'restriction'or'projection'matrices).However,the current finite element
>> can
>> > either notdefine thissort of operation,orit has notyet been
>> implemented.
>> > |
>> >
>> >
>> > Which function do I have to modify to make that work?
>>
>> You will have to implement the restriction matrices for the Bernstein
>> element.
>> These are usually set in the constructor of the finite element classes.
>> Take a
>> look at the other FE classes for examples.
>>
>> The error you see is simply because you are accessing a field that has
>> never
>> been initialized.
>>
>> Best
>> W.
>>
>> --
>> ------------------------------------------------------------------------
>> Wolfgang Bangerth email: bang...@colostate.edu
>> www: http://www.math.colostate.edu/~bangerth/
>>
>>
--
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/fa00b08a-8f8f-44d4-bfba-16d0c3735cf4%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
