Dear dealii community, I am working on a gradient enhanced crystal plasticity code where the displacement field is continuous but the slip fields can be discontinuous. To this end, I am using a hp::FECollection<dim>. For example, in the case of two crystals with one slip system I have a hp::FECollection<2> with
fe_collection[0] = [FE_Q<2> FE_Q<2> FE_Q<2> FE_Nothing<2>] fe_collection[1] = [FE_Q<2> FE_Q<2> FE_Nothing<2> FE_Q<2> ] where the first two FE_Q<2> correspond to the displacement field in 2D and the latter to the slip field. My reference problem consist of the unit square under simple shear with periodic boundary conditions on the x-direction. When I call the DoFTools::make_periodicity_constraints<dim, dim> method, I am getting the following error when I have subdomains with different values of fe_index An error occurred in line <1690> of file </calculate/temp/iwtm009/spack-stage-dealii-9.3.3-h7vrbckbqhxjllhyiwpmxjdp3242wjvv/spack-src/include/deal.II/dofs/dof_accessor.templates.h> in function const dealii::FiniteElement<dim, spacedim>& dealii::DoFAccessor<structdim, dim, spacedim, lda>::get_fe(unsigned int) const [with int structdim = 1; int dim = 2; int spacedim = 2; bool level_dof_access = false] The violated condition was: fe_index_is_active(fe_index) == true Additional information: This function can only be called for active FE indices I wrote a minimum working example based on my code to showcase the error, which I have attached to this post. The executable takes a true or false argument to assign different material identifiers to different subdomains or not. If there are no subdomains and therefore, and more importantly, the size of the hp::FECollection<dim> is 1, the code runs with no problem. My question is if my algorithm is wrong or does the method DoFTools::make_periodicity_constraints<dim, dim> does not currently support hp::FECollection<dim> with sizes bigger than one. Cheers, Jose -- 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 on the web visit https://groups.google.com/d/msgid/dealii/b4ea5417-bc21-4541-ac57-e5d71383a7b9n%40googlegroups.com.
#include <deal.II/base/index_set.h> #include <deal.II/distributed/tria.h> #include <deal.II/dofs/dof_handler.h> #include <deal.II/dofs/dof_renumbering.h> #include <deal.II/dofs/dof_tools.h> #include <deal.II/fe/fe_q.h> #include <deal.II/fe/fe_nothing.h> #include <deal.II/fe/fe_system.h> #include <deal.II/fe/fe_values_extractors.h> #include <deal.II/grid/tria.h> #include <deal.II/grid/grid_generator.h> #include <deal.II/grid/grid_out.h> #include <deal.II/grid/grid_tools.h> #include <deal.II/hp/fe_collection.h> #include <string> namespace Tests { template<int dim> class MakePeriodicityConstraints { public: MakePeriodicityConstraints(const bool flag_set_material_ids); void run(); private: dealii::parallel::distributed::Triangulation<dim> triangulation; dealii::DoFHandler<dim> dof_handler; dealii::hp::FECollection<dim> fe_collection; const double width; const double height; const double n_fields_per_domain; double n_domains; const bool flag_set_material_ids; void make_grid(); void setup_dofs(); }; template<int dim> MakePeriodicityConstraints<dim>::MakePeriodicityConstraints( const bool flag_set_material_ids) : triangulation(MPI_COMM_WORLD, typename dealii::Triangulation<dim>::MeshSmoothing( dealii::Triangulation<dim>::smoothing_on_refinement | dealii::Triangulation<dim>::smoothing_on_coarsening)), dof_handler(triangulation), width(1.0), height(1.0), n_fields_per_domain(2), flag_set_material_ids(flag_set_material_ids) {} template<int dim> void MakePeriodicityConstraints<dim>::run() { make_grid(); setup_dofs(); } template<int dim> void MakePeriodicityConstraints<dim>::make_grid() { std::vector<unsigned int> repetitions(2, 10); dealii::GridGenerator::subdivided_hyper_rectangle( triangulation, repetitions, dealii::Point<dim>(0,0), dealii::Point<dim>(width, height), true); std::vector<dealii::GridTools::PeriodicFacePair< typename dealii::parallel::distributed::Triangulation<dim>::cell_iterator>> periodicity_vector; dealii::GridTools::collect_periodic_faces(triangulation, 0, 1, 0, periodicity_vector); triangulation.add_periodicity(periodicity_vector); triangulation.refine_global(1); // Set material ids for (const auto &cell : triangulation.active_cell_iterators()) if (cell->is_locally_owned()) { if (flag_set_material_ids) { if (std::fabs(cell->center()[1]) < height/3.0) cell->set_material_id(0); else if (std::fabs(cell->center()[1]) < 2.0*height/3.0) cell->set_material_id(1); else cell->set_material_id(2); } else cell->set_material_id(0); } // Set the active finite elemente index of each cell for (const auto &cell : dof_handler.active_cell_iterators()) if (cell->is_locally_owned()) cell->set_active_fe_index(cell->material_id()); // Get the number of domains std::set<dealii::types::material_id> domain_set; for (const auto &cell : triangulation.active_cell_iterators()) if (cell->is_locally_owned()) if (!domain_set.count(cell->material_id())) domain_set.emplace(cell->material_id()); domain_set = dealii::Utilities::MPI::compute_set_union(domain_set, MPI_COMM_WORLD); n_domains = domain_set.size(); } template<int dim> void MakePeriodicityConstraints<dim>::setup_dofs() { for (dealii::types::material_id i = 0; i < n_domains; ++i) { std::vector<const dealii::FiniteElement<dim>*> finite_elements; if (false) for (dealii::types::material_id j = 0; j < n_domains; ++j) for (unsigned int k = 0; k < dim; ++k) if (i == j) finite_elements.push_back( new dealii::FE_Q<dim>(2)); else finite_elements.push_back( new dealii::FE_Nothing<dim>()); else for (unsigned int j = 0; j < dim; ++j) finite_elements.push_back( new dealii::FE_Q<dim>(2)); for (dealii::types::material_id j = 0; j < n_domains; ++j) for (unsigned int k = 0; k < n_fields_per_domain; ++k) if (i == j) finite_elements.push_back( new dealii::FE_Q<dim>(1)); else finite_elements.push_back( new dealii::FE_Nothing<dim>()); fe_collection.push_back( dealii::FESystem<dim>( finite_elements, std::vector<unsigned int>(finite_elements.size(), 1))); for (auto finite_element: finite_elements) delete finite_element; finite_elements.clear(); } dof_handler.distribute_dofs(fe_collection); dealii::DoFRenumbering::Cuthill_McKee(dof_handler); dealii::IndexSet locally_owned_dofs = dof_handler.locally_owned_dofs(); dealii::IndexSet locally_relevant_dofs; dealii::DoFTools::extract_locally_relevant_dofs( dof_handler, locally_relevant_dofs); dealii::AffineConstraints<double> hanging_node_constraints; dealii::AffineConstraints<double> affine_constraints; hanging_node_constraints.clear(); { hanging_node_constraints.reinit(locally_relevant_dofs); dealii::DoFTools::make_hanging_node_constraints( dof_handler, hanging_node_constraints); } hanging_node_constraints.close(); affine_constraints.clear(); { affine_constraints.reinit(locally_relevant_dofs); affine_constraints.merge(hanging_node_constraints); } affine_constraints.close(); std::vector< dealii::GridTools:: PeriodicFacePair<typename dealii::DoFHandler<dim>::cell_iterator>> periodicity_vector; dealii::GridTools::collect_periodic_faces( dof_handler, 0, 1, 0, periodicity_vector); dealii::AffineConstraints<double> tmp_affine_constraints; tmp_affine_constraints.clear(); { tmp_affine_constraints.reinit(locally_relevant_dofs); tmp_affine_constraints.merge(hanging_node_constraints); dealii::DoFTools::make_periodicity_constraints<dim, dim>( periodicity_vector, tmp_affine_constraints); } tmp_affine_constraints.close(); affine_constraints.merge( tmp_affine_constraints, dealii::AffineConstraints<double>::MergeConflictBehavior::right_object_wins); } } // namespace Test int main(int argc, char *argv[]) { try { dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization( argc, argv, dealii::numbers::invalid_unsigned_int); bool flag_set_material_ids; switch (argc) { case 1: flag_set_material_ids = true; break; case 2: { const std::string arg(argv[1]); if (arg == "true") flag_set_material_ids = true; else if (arg == "false") flag_set_material_ids = false; else AssertThrow(false, dealii::ExcNotImplemented()); } break; default: AssertThrow(false, dealii::ExcNotImplemented()); } Tests::MakePeriodicityConstraints<2> test(flag_set_material_ids); test.run(); } 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; }