search

dealii / dealii

The development repository for the deal.II finite element library
https://www.dealii.org
Other
1.34k stars 741 forks source link

TransfiniteInterpolationManifold::initialize with <2,3> triangulations #12029

Open dianeguignard opened 3 years ago

dianeguignard commented 3 years ago

Dear all,

I encounter a dimension issue when using the initialize function of TransfiniteInterpolationManifold for <2,3> triangulations.

The error can be reproduced using the code below (in practice, I am using this for the mesh of a circular domain, but since hyper_ball is not implemented for <2,3>, it requires a few more lines of code). Note that I did not have this issue with a previous version of deal.ii (I tested on version 9.1.1), which does not have the function GridTools::affine_cell_approximation() producing the assertion error.

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

using namespace dealii;

int main()
{

    Triangulation <2,3> tria;
    GridGenerator::hyper_cube(tria);

    TransfiniteInterpolationManifold<2,3> inner_manifold;
    inner_manifold.initialize(tria);

    return 0;
}
drwells commented 3 years ago

Definitely a bug - thanks for reporting this.

The problem stems from https://github.com/dealii/dealii/blob/master/include/deal.II/fe/mapping_q_internal.h#L899: it looks like we set up an affine approximation but with dim = 2 and spacedim = 3 we supply all 9 points on the quadrilateral (instead of just 4), so a bounds check fails.

bangerth commented 3 years ago

For easier bug hunting, can one of you attach the actual error message?

dianeguignard commented 3 years ago

Please find the error message below. As David said, 9 points are passed to the function affine_cell_approximation (variable vertices) while GeometryInfo<2>::vertices_per_cell = 4. However, it seems that vertices[v] is only used for 0 <= v GeometryInfo< dim>::vertices_per_cell in affine_cell_approximation.

An error occurred in line <289> of file </home/dguignar/Desktop/Prog/dealii/dealii-master/source/grid/grid_tools.cc> in function
    std::pair<dealii::DerivativeForm<1, dim, spacedim>, dealii::Tensor<1, spacedim> > dealii::GridTools::affine_cell_approximation(const dealii::ArrayView<const dealii::Point<spacedim> >&) [with int dim = 2; int spacedim = 3]
The violated condition was: 
    static_cast<typename ::dealii::internal::argument_type<void( typename std::common_type<decltype(vertices.size()), decltype(GeometryInfo<dim>::vertices_per_cell)>::type)>::type>(vertices.size()) == static_cast<typename ::dealii::internal::argument_type<void( typename std::common_type<decltype(vertices.size()), decltype(GeometryInfo<dim>::vertices_per_cell)>::type)>::type>(GeometryInfo<dim>::vertices_per_cell)
Additional information: 
    Dimension 9 not equal to 4.

Stacktrace:
-----------
#0  /home/dguignar/Desktop/Prog/projects/local/lib/libdeal_II.g.so.9.3.0-pre: std::pair<dealii::DerivativeForm<1, 2, 3, double>, dealii::Tensor<1, 3, double> > dealii::GridTools::affine_cell_approximation<2, 3>(dealii::ArrayView<dealii::Point<3, double> const, dealii::MemorySpace::Host> const&)
#1  /home/dguignar/Desktop/Prog/projects/local/lib/libdeal_II.g.so.9.3.0-pre: dealii::internal::MappingQGenericImplementation::InverseQuadraticApproximation<2, 3>::InverseQuadraticApproximation(std::vector<dealii::Point<3, double>, std::allocator<dealii::Point<3, double> > > const&, std::vector<dealii::Point<2, double>, std::allocator<dealii::Point<2, double> > > const&)
#2  /home/dguignar/Desktop/Prog/projects/local/lib/libdeal_II.g.so.9.3.0-pre: void __gnu_cxx::new_allocator<dealii::internal::MappingQGenericImplementation::InverseQuadraticApproximation<2, 3> >::construct<dealii::internal::MappingQGenericImplementation::InverseQuadraticApproximation<2, 3>, std::vector<dealii::Point<3, double>, std::allocator<dealii::Point<3, double> > >&, std::vector<dealii::Point<2, double>, std::allocator<dealii::Point<2, double> > >&>(dealii::internal::MappingQGenericImplementation::InverseQuadraticApproximation<2, 3>*, std::vector<dealii::Point<3, double>, std::allocator<dealii::Point<3, double> > >&, std::vector<dealii::Point<2, double>, std::allocator<dealii::Point<2, double> > >&)
#3  /home/dguignar/Desktop/Prog/projects/local/lib/libdeal_II.g.so.9.3.0-pre: void std::allocator_traits<std::allocator<dealii::internal::MappingQGenericImplementation::InverseQuadraticApproximation<2, 3> > >::construct<dealii::internal::MappingQGenericImplementation::InverseQuadraticApproximation<2, 3>, std::vector<dealii::Point<3, double>, std::allocator<dealii::Point<3, double> > >&, std::vector<dealii::Point<2, double>, std::allocator<dealii::Point<2, double> > >&>(std::allocator<dealii::internal::MappingQGenericImplementation::InverseQuadraticApproximation<2, 3> >&, dealii::internal::MappingQGenericImplementation::InverseQuadraticApproximation<2, 3>*, std::vector<dealii::Point<3, double>, std::allocator<dealii::Point<3, double> > >&, std::vector<dealii::Point<2, double>, std::allocator<dealii::Point<2, double> > >&)
#4  /home/dguignar/Desktop/Prog/projects/local/lib/libdeal_II.g.so.9.3.0-pre: void std::vector<dealii::internal::MappingQGenericImplementation::InverseQuadraticApproximation<2, 3>, std::allocator<dealii::internal::MappingQGenericImplementation::InverseQuadraticApproximation<2, 3> > >::_M_realloc_insert<std::vector<dealii::Point<3, double>, std::allocator<dealii::Point<3, double> > >&, std::vector<dealii::Point<2, double>, std::allocator<dealii::Point<2, double> > >&>(__gnu_cxx::__normal_iterator<dealii::internal::MappingQGenericImplementation::InverseQuadraticApproximation<2, 3>*, std::vector<dealii::internal::MappingQGenericImplementation::InverseQuadraticApproximation<2, 3>, std::allocator<dealii::internal::MappingQGenericImplementation::InverseQuadraticApproximation<2, 3> > > >, std::vector<dealii::Point<3, double>, std::allocator<dealii::Point<3, double> > >&, std::vector<dealii::Point<2, double>, std::allocator<dealii::Point<2, double> > >&)
#5  /home/dguignar/Desktop/Prog/projects/local/lib/libdeal_II.g.so.9.3.0-pre: void std::vector<dealii::internal::MappingQGenericImplementation::InverseQuadraticApproximation<2, 3>, std::allocator<dealii::internal::MappingQGenericImplementation::InverseQuadraticApproximation<2, 3> > >::emplace_back<std::vector<dealii::Point<3, double>, std::allocator<dealii::Point<3, double> > >&, std::vector<dealii::Point<2, double>, std::allocator<dealii::Point<2, double> > >&>(std::vector<dealii::Point<3, double>, std::allocator<dealii::Point<3, double> > >&, std::vector<dealii::Point<2, double>, std::allocator<dealii::Point<2, double> > >&)
#6  /home/dguignar/Desktop/Prog/projects/local/lib/libdeal_II.g.so.9.3.0-pre: dealii::TransfiniteInterpolationManifold<2, 3>::initialize(dealii::Triangulation<2, 3> const&)
#7  ./Test: main