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

Behavior of FEEvaluation::reinit() #13687

Closed
peterrum opened this issue May 7, 2022 · 0 comments · Fixed by #14565
Closed

Behavior of FEEvaluation::reinit() #13687

peterrum opened this issue May 7, 2022 · 0 comments · Fixed by #14565

Comments

@peterrum
Copy link
Member

peterrum commented May 7, 2022

The function FEEvaluation::reinit() can be called either with a cell-batch index or with an array of indices of cells that should make up a new macro cell. Calling the first function and after that the second one works. But it does not work in the opposite way around, due to the assert:

--------------------------------------------------------
An error occurred in line <6765> of file </home/munch/sw-index/dealii2/include/deal.II/matrix_free/fe_evaluation.h> in function
    void dealii::FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number, VectorizedArrayType>::reinit(unsigned int) [with int dim = 2; int fe_degree = -1; int n_q_points_1d = 0; int n_components_ = 1; Number = double; VectorizedArrayType = dealii::VectorizedArray<double, 8>]
The violated condition was: 
    this->mapped_geometry == nullptr
Additional information: 
    FEEvaluation was initialized without a matrix-free object. Integer
    indexing is not possible

Stacktrace:
-----------
#0  ./main: dealii::FEEvaluation<2, -1, 0, 1, double, dealii::VectorizedArray<double, 8ul> >::reinit(unsigned int)
#1  ./main: main

What do we want? Should we allow switching between modi? Or should we add an additional assert in the second function?

Minimal test:

#include <deal.II/fe/fe_q.h>

#include <deal.II/grid/grid_generator.h>

#include <deal.II/lac/affine_constraints.h>

#include <deal.II/matrix_free/fe_evaluation.h>
#include <deal.II/matrix_free/matrix_free.h>

using namespace dealii;

int
main()
{
  const unsigned int dim    = 2;
  const unsigned int degree = 2;

  FE_Q<dim> fe(degree);

  std::cout << fe.n_dofs_per_cell() << std::endl;

  Triangulation<dim> tria;
  GridGenerator::hyper_cube(tria);

  DoFHandler<dim> dof_handler(tria);
  dof_handler.distribute_dofs(fe);

  AffineConstraints<double> constraints;

  MatrixFree<dim, double> mf;
  mf.reinit(dof_handler, constraints, QGauss<dim>(degree + 1));

  FEEvaluation<dim, -1, 0, 1, double> phi(mf);

  phi.reinit(0);
  phi.reinit(std::array<unsigned int, VectorizedArray<double>::size()>{{0}});
  phi.reinit(0);
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging a pull request may close this issue.

1 participant