Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

MGTwoLevelTransfer (p): add fallback FE case #12924

Draft
wants to merge 2 commits into
base: master
Choose a base branch
from

Conversation

peterrum
Copy link
Member

@peterrum peterrum commented Nov 7, 2021

... to allow more type of FEs (that are currently not supported by MatrixFree, e.g. RT or Nedelec).

Most changed lines are related to changed indentation (due to extra if-statements). These are the extra operations:

{
lexicographic_numbering_fine[fe_index_pair.second].resize(
dof_handler_fine.get_fe(fe_index_pair.first.second)
.n_dofs_per_cell());
lexicographic_numbering_coarse[fe_index_pair.second].resize(
dof_handler_coarse.get_fe(fe_index_pair.first.first)
.n_dofs_per_cell());
for (unsigned int i = 0;
i <
lexicographic_numbering_fine[fe_index_pair.second].size();
++i)
lexicographic_numbering_fine[fe_index_pair.second][i] = i;
for (unsigned int i = 0;
i < lexicographic_numbering_coarse[fe_index_pair.second]
.size();
++i)
lexicographic_numbering_coarse[fe_index_pair.second][i] = i;
}
}

{
transfer.n_components = 1;
const auto &fe_fine =
dof_handler_fine.get_fe(fe_index_pair.second);
const auto &fe_coarse =
dof_handler_coarse.get_fe(fe_index_pair.first);
{
FullMatrix<double> matrix(fe_fine.n_dofs_per_cell(),
fe_coarse.n_dofs_per_cell());
FETools::get_projection_matrix(fe_coarse, fe_fine, matrix);
transfer.schemes[fe_index_no].prolongation_matrix.resize(
fe_fine.n_dofs_per_cell() * fe_coarse.n_dofs_per_cell());
for (unsigned int i = 0, k = 0;
i < fe_coarse.n_dofs_per_cell();
++i)
for (unsigned int j = 0; j < fe_fine.n_dofs_per_cell();
++j, ++k)
transfer.schemes[fe_index_no].prolongation_matrix[k] =
matrix(j, i);
}
{
FullMatrix<double> matrix(fe_coarse.n_dofs_per_cell(),
fe_fine.n_dofs_per_cell());
FETools::get_projection_matrix(fe_fine, fe_coarse, matrix);
transfer.schemes[fe_index_no].restriction_matrix.resize(
fe_fine.n_dofs_per_cell() * fe_coarse.n_dofs_per_cell());
for (unsigned int i = 0, k = 0;
i < fe_coarse.n_dofs_per_cell();
++i)
for (unsigned int j = 0; j < fe_fine.n_dofs_per_cell();
++j, ++k)
transfer.schemes[fe_index_no].restriction_matrix[k] =
matrix(i, j);
}
}

depends on #12923

Copy link
Member

@kronbichler kronbichler left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks like a good extension to me, but there are some questions I have regarding vector-valued cases (did we support them before)? See below.

Comment on lines -1334 to +1346
const auto dummy_quadrature =
reference_cell.template get_gauss_type_quadrature<dim>(1);
internal::MatrixFreeFunctions::ShapeInfo<Number> shape_info;
shape_info.reinit(dummy_quadrature, fe_fine, 0);
lexicographic_numbering = shape_info.lexicographic_numbering;
AssertThrow(false, ExcNotImplemented());
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there an assumption about a single component somewhere else or are you now restricting to a smaller class of finite elements in the hypercube case? The original code of ShapeInfo allowed multiple components, but you ask for n_components()==1 now, which is more restrictive than fe.n_base_elements()==1 I found as assertions elsewhere.

Comment on lines +1452 to +1465
std::vector<unsigned int> renumbering(fe->n_dofs_per_cell());
{
AssertIndexRange(fe->n_dofs_per_vertex(), 2);
renumbering[0] = 0;
for (unsigned int i = 0; i < fe->dofs_per_line; ++i)
renumbering[i + fe->n_dofs_per_vertex()] =
GeometryInfo<1>::vertices_per_cell *
fe->n_dofs_per_vertex() +
i;
if (fe->n_dofs_per_vertex() > 0)
renumbering[fe->n_dofs_per_cell() -
fe->n_dofs_per_vertex()] =
fe->n_dofs_per_vertex();
}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I know this is pre-existing, but isn't this the regular hierarchic-to-lexicographic renumbering you would get from ShapeInfo::lexicographic_numbering or maybe the inverse? I know it is a bit of a hammer to call that function, but this is not performance-critical code and 1D anyway, and I find it more documenting, and is what we do below.

Comment on lines +1824 to +1843
else
{
const auto dummy_quadrature =
reference_cell.template get_gauss_type_quadrature<dim>(1);

internal::MatrixFreeFunctions::ShapeInfo<Number> shape_info;
shape_info.reinit(dummy_quadrature,
dof_handler_fine.get_fe(
fe_index_pair.first.second),
0);
lexicographic_numbering_fine[fe_index_pair.second] =
shape_info.lexicographic_numbering;

shape_info.reinit(dummy_quadrature,
dof_handler_coarse.get_fe(
fe_index_pair.first.first),
0);
lexicographic_numbering_coarse[fe_index_pair.second] =
shape_info.lexicographic_numbering;
}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't see how this else clause differs from the other places, except that it uses a different dummy quadrature. But we never use anything from the quadrature (as it is a dummy), so I don't see why we can't use the one from the if clause also for non-hypercubes. Again, not the change in question here, but maybe we can fix it. If not, I would prefer to only set the quadrature formula in an if/else setup.

@peterrum peterrum added this to the Release 9.5 milestone Apr 22, 2022
@peterrum peterrum changed the title MGTwoLevelTransfer (p): add fallback FE case [WIP] MGTwoLevelTransfer (p): add fallback FE case Dec 10, 2022
@kronbichler
Copy link
Member

@peterrum what should we do about this PR? I suggest we move the target milestone, it seems like too much work the release now.

@peterrum peterrum removed this from the Release 9.5 milestone Jun 19, 2023
@peterrum
Copy link
Member Author

@kronbichler I have removed the label.

@tamiko tamiko changed the title [WIP] MGTwoLevelTransfer (p): add fallback FE case MGTwoLevelTransfer (p): add fallback FE case Jul 4, 2023
@tamiko tamiko marked this pull request as draft July 4, 2023 18:02
@adamqc
Copy link
Contributor

adamqc commented Nov 18, 2023

Hi @peterrum, I have reviewed this PR, however, given my limited familiarity with MGTransfers, I'm uncertain about the next steps. Could you please provide some guidance on how to proceed and complete it?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

3 participants