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;
}

Reply via email to