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