Hello Davit,
This does not work because FE_Nothing does not have a support points. I
needed something similar so I implemented it in my code, you can take a
look at
https://github.com/adamantine-sim/adamantine/blob/efde0a1818ae4f2c476aee56c55a187c7b3ece08/source/experimental_data_utils.cc#L21-L65
You may need to slightly adapt it for your use case.
Best,
Bruno
On Thursday, March 20, 2025 at 11:53:40 AM UTC-4 dav.gyu...@gmail.com wrote:
Hi all,
I have encountered a problem when trying to obtain the support points from
a dof_handler.
The Finite Element used is an hp::FECollection<dim> with two elements: 1)
FE_System<dim> fe(FE_Q<dim>(1), dim, FE_Q<dim>(1), 1), and 2)
fe_nothing(FE_Nothing<dim>(), dim, FE_Nothing<dim>(), 1). The domain is
split into two regions - active and inactive - hence the use of FE_Nothing.
Then, calling DoFTools::map_dofs_to_support_points(MappingQ1<dim>(),
dof_handler) gives an error in *debug mode*, but *no *errors in *release
mode*!
The error message is as follows:
--------------------------------------------------------
Calling MPI_Abort now.
To break execution in a GDB session, execute 'break MPI_Abort' before
running. You can also put the following into your ~/.gdbinit:
set breakpoint pending on
break MPI_Abort
set breakpoint pending auto
--------------------------------------------------------
An error occurred in line <2791> of file
</proj/petsc/dealii-9.5.2/source/fe/fe_values.cc> in function
dealii::FEValuesBase<dim, spacedim>::FEValuesBase(unsigned int,
unsigned int, dealii::UpdateFlags, const dealii::Mapping<dim, spacedim>&,
const dealii::FiniteElement<dim, spacedim>&) [with int dim = 3; int
spacedim = 3]
The violated condition was:
n_q_points > 0
Additional information:
There is nothing useful you can do with an FEValues object when using
a quadrature formula with zero quadrature points!
Stacktrace:
-----------
#0 /proj/petsc/dealii_debug/lib/libdeal_II.g.so.9.5.2:
dealii::FEValuesBase<3, 3>::FEValuesBase(unsigned int, unsigned int,
dealii::UpdateFlags, dealii::Mapping<3, 3> const&, dealii::FiniteElement<3,
3> const&)
#1 /proj/petsc/dealii_debug/lib/libdeal_II.g.so.9.5.2: dealii::FEValues<3,
3>::FEValues(dealii::Mapping<3, 3> const&, dealii::FiniteElement<3, 3>
const&, dealii::Quadrature<3> const&, dealii::UpdateFlags)
#2 /proj/petsc/dealii_debug/lib/libdeal_II.g.so.9.5.2: dealii::FEValues<3,
3>::FEValues(dealii::Mapping<3, 3> const&, dealii::FiniteElement<3, 3>
const&, dealii::hp::QCollection<3> const&, dealii::UpdateFlags)
#3 /proj/petsc/dealii_debug/lib/libdeal_II.g.so.9.5.2:
std::_MakeUniq<dealii::FEValues<3, 3> >::__single_object
std::make_unique<dealii::FEValues<3, 3>, dealii::Mapping<3, 3> const&,
dealii::FiniteElement<3, 3> const&, dealii::hp::QCollection<3> const&,
dealii::UpdateFlags const&>(dealii::Mapping<3, 3> const&,
dealii::FiniteElement<3, 3> const&, dealii::hp::QCollection<3> const&,
dealii::UpdateFlags const&)
#4 /proj/petsc/dealii_debug/lib/libdeal_II.g.so.9.5.2:
dealii::hp::FEValuesBase<3, 3, dealii::FEValues<3, 3>
>::select_fe_values(unsigned int, unsigned int, unsigned int)
#5 /proj/petsc/dealii_debug/lib/libdeal_II.g.so.9.5.2: void
dealii::hp::FEValues<3,
3>::reinit<false>(dealii::TriaIterator<dealii::DoFCellAccessor<3, 3, false>
> const&, unsigned int, unsigned int, unsigned int)
#6 /proj/petsc/dealii_debug/lib/libdeal_II.g.so.9.5.2:
#7 /proj/petsc/dealii_debug/lib/libdeal_II.g.so.9.5.2: std::map<unsigned
int, dealii::Point<3, double>, std::less<unsigned int>,
std::allocator<std::pair<unsigned int const, dealii::Point<3, double> > > >
dealii::DoFTools::map_dofs_to_support_points<3, 3>(dealii::Mapping<3, 3>
const&, dealii::DoFHandler<3, 3> const&, dealii::ComponentMask const&)
#8 test: main
--------------------------------------------------------
Which I assume has to do with the presence of FE_Nothing.
Minimum code to reproduce the issue:
int main(int argc, char *argv[])
{
try
{
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
MPI_Comm mpi_communicator(MPI_COMM_WORLD);
const unsigned int dim = 3;
parallel::distributed::Triangulation<dim> triangulation
(mpi_communicator);
GridGenerator::hyper_cube(triangulation, 0, 1);
triangulation.refine_global(4);
FESystem<dim> fe(FE_Q<dim>(1), dim, FE_Q<dim>(1), 1);
FESystem<dim> fe_nothing(FE_Nothing<dim>(), dim, FE_Nothing<
dim>(), 1);
hp::FECollection<dim> fe_collection;
fe_collection.push_back(fe);
fe_collection.push_back(fe_nothing);
DoFHandler<dim> dof_handler(triangulation);
for (const auto &cell : dof_handler.active_cell_iterators())
{
if (cell->is_locally_owned())
{
if (cell->center()[0] < 0.5)
cell->set_active_fe_index(0);
else
cell->set_active_fe_index(1);
}
}
dof_handler.distribute_dofs(fe_collection);
MappingQ1<dim> mapping;
std::map<types::global_dof_index, Point<dim> > support_points;
support_points = DoFTools::map_dofs_to_support_points(mapping,
dof_handler);
}
catch (std::exception &exc)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
I would appreciate any help! Let me know if there is any other information
required about the case.
Best wishes,
Davit Gyulamiryan
--
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 visit
https://groups.google.com/d/msgid/dealii/acef7e0f-9dbb-4b9a-bdd8-a956e9519fa5n%40googlegroups.com.
