From 0850769b5dc75b6e18539a10accd7d67b947f810 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Sat, 25 Jul 2015 18:00:06 +0200 Subject: [PATCH] Fixed bug in MappingFEField. Closes #1138. --- .../deal.II/numerics/vector_tools.templates.h | 82 ++++++++++++++++--- 1 file changed, 70 insertions(+), 12 deletions(-) diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index 50b9a2cfca..9a207bf594 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -7177,33 +7177,91 @@ namespace VectorTools // Once we have this, interpolate with the given finite element // to get a Mapping which is interpolatory at the support points // of FE_Q(fe.degree()) - FESystem feq(FE_Q(fe.degree), spacedim); + const FESystem *fe_system = dynamic_cast *>(&fe); + Assert(fe_system, ExcNotImplemented()); + unsigned int degree = numbers::invalid_unsigned_int; + + // Get information about the blocks + for (unsigned int i=0; icomponent_to_base_index(i).first; + Assert(degree == numbers::invalid_unsigned_int || + degree == fe_system->base_element(base_i).degree, + ExcNotImplemented()); + degree = fe_system->base_element(base_i).degree; + } + + FESystem feq(FE_Q(degree), spacedim); DH dhq(dh.get_tria()); dhq.distribute_dofs(feq); Vector eulerq(dhq.n_dofs()); const ComponentMask maskq(spacedim, true); get_position_vector(dhq, eulerq); - FullMatrix transfer(fe.dofs_per_cell); + FullMatrix transfer(fe.dofs_per_cell, feq.dofs_per_cell); + FullMatrix local_transfer(feq.dofs_per_cell); const std::vector > &points = feq.get_unit_support_points(); - // Here construct the interpolation matrix from FE_Q^spacedim to - // the FiniteElement used by euler_dof_handler. + // Here we construct the interpolation matrix from + // FE_Q^spacedim to the FiniteElement used by + // euler_dof_handler. // - // The interpolation matrix is then passed to the - // VectorTools::interpolate() function to generate + // In order to construct such interpolation matrix, we have to + // solve the following system: + // + // v_j phi_j(q_i) = w_k psi_k(q_i) = w_k delta_ki = w_i + // + // where psi_k are the basis functions for fe_q, and phi_i are + // the basis functions of the target space. + // + // Morally, we should invert the matrix T_ij = phi_i(q_j), + // however in general this matrix is not invertible, since + // there may by components which do not contribute to the + // displacement vector. Since we are not interested in those + // components, we construct a square matrix with the same + // number of components of the FE_Q system. The FE_Q system + // was constructed above in such a way that the polynomial + // degree of the FE_Q system and that of the given FE are the + // same on the cell, which should guarantee that, for the + // displacement components only, the interpolation matrix is + // invertible. We construct a mapping between indices first, + // and check that this is the case. If not, we bail out, not + // knowing what to do in this case. + + std::vector fe_to_feq(fe.dofs_per_cell, numbers::invalid_unsigned_int); + unsigned int index=0; for (unsigned int i=0; i