Dear dealii community,
I am still trying to decipher what causes the error. It occurs in the first
assert of make_periodicity_constraints
<https://www.dealii.org/9.3.3/doxygen/deal.II/dof__tools__constraints_8cc_source.html#l02288>
:
AssertDimension(
face_1->get_fe(face_1->nth_active_fe_index(0)).n_unique_faces(), 1);
So I created vectors to store the active_fe_index and the
nth_active_fe_index(0), replicated the fe_index_is_active assert and the
get_fe() call:
active_fe_index.reinit(triangulation.n_active_cells());
nth_active_fe_index.reinit(triangulation.n_active_cells());
active_fe_index = -1.0;
nth_active_fe_index = -1.0;
for (const auto &cell :
dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
{
active_fe_index(cell->active_cell_index()) = cell->active_fe_index();
if (cell->at_boundary())
for (const auto &face_id : cell->face_indices())
if (cell->face(face_id)->at_boundary())
{
AssertThrow(cell->face(face_id)->get_active_fe_indices().size() == 1, dealii
::ExcMessage("Error!"));
nth_active_fe_index(cell->active_cell_index()) = cell->face(face_id)->
nth_active_fe_index(0);
AssertThrow(cell->face(face_id)->fe_index_is_active(cell->face(face_id)->
nth_active_fe_index(0)), dealii::ExcMessage("Error!"));
cell->face(face_id)->get_fe(cell->face(face_id)->nth_active_fe_index(0));
}
}
I don't get any error from this loop and the from the visual inspection of
.vtu file containing both vectors one can see that they indeed coincide.
Maybe the problem is coming from GridTools::collect_periodic_faces as the
face iterators are created by it?
Furthermore I also noticed that if no global refinement is performed the
executable runs without error in serial but still gets the error in
parallel.
Attached is the updated MWE. A true/false argument can be passed to the
executable to perform one global refinement or not.
Cheers,
Jose
On Wednesday, June 29, 2022 at 4:45:04 PM UTC+2 jose.a...@gmail.com wrote:
> 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/41e28efd-4d16-4774-8aaa-b7eb5946b232n%40googlegroups.com.
#include <deal.II/base/conditional_ostream.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/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 <deal.II/numerics/data_out.h>
#include <string>
namespace Tests
{
template<int dim>
class MWE
{
public:
MWE(const bool flag_refine_global);
void run();
private:
dealii::ConditionalOStream pcout;
dealii::parallel::distributed::Triangulation<dim>
triangulation;
dealii::DoFHandler<dim> dof_handler;
dealii::hp::FECollection<dim> fe_collection;
dealii::Vector<float> active_fe_index;
dealii::Vector<float> nth_active_fe_index;
const double height;
const double width;
unsigned int n_crystals;
const unsigned int n_slips;
const bool flag_refine_global;
void make_grid();
void setup_dofs();
void output();
};
template<int dim>
MWE<dim>::MWE(const bool flag_refine_global)
:
pcout(std::cout,
dealii::Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0),
triangulation(MPI_COMM_WORLD,
typename dealii::Triangulation<dim>::MeshSmoothing(
dealii::Triangulation<dim>::smoothing_on_refinement |
dealii::Triangulation<dim>::smoothing_on_coarsening)),
dof_handler(triangulation),
height(1.0),
width(1.0),
n_slips(2),
flag_refine_global(flag_refine_global)
{}
template<int dim>
void MWE<dim>::run()
{
make_grid();
setup_dofs();
}
template<int dim>
void MWE<dim>::make_grid()
{
std::vector<unsigned int> repetitions(2, 20);
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(flag_refine_global);
// Set material ids
for (const auto &cell : triangulation.active_cell_iterators())
if (cell->is_locally_owned())
{
if (std::fabs(cell->center()[1]) < height/2.0)
cell->set_material_id(0);
else
cell->set_material_id(1);
}
this->pcout << "Triangulation:"
<< std::endl
<< " Number of active cells = "
<< triangulation.n_global_active_cells()
<< std::endl << std::endl;
}
template<int dim>
void MWE<dim>::setup_dofs()
{
std::set<dealii::types::material_id> crystal_id_set;
for (const auto &cell : triangulation.active_cell_iterators())
if (cell->is_locally_owned())
if (!crystal_id_set.count(cell->material_id()))
crystal_id_set.emplace(cell->material_id());
crystal_id_set =
dealii::Utilities::MPI::compute_set_union(crystal_id_set,
MPI_COMM_WORLD);
n_crystals = crystal_id_set.size();
// 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());
const unsigned int displacement_finite_element_degree = 2;
const unsigned int slip_finite_element_degree = 1;
for (dealii::types::material_id i = 0; i < n_crystals; ++i)
{
std::vector<const dealii::FiniteElement<dim>*> finite_elements;
// A
if (false)
for (dealii::types::material_id j = 0; j < n_crystals; ++j)
for (unsigned int k = 0; k < dim; ++k)
if (i == j)
finite_elements.push_back(
new dealii::FE_Q<dim>(displacement_finite_element_degree));
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>(displacement_finite_element_degree));
// B
for (dealii::types::material_id j = 0; j < n_crystals; ++j)
for (unsigned int k = 0; k < n_slips; ++k)
if (i == j)
finite_elements.push_back(
new dealii::FE_Q<dim>(slip_finite_element_degree));
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);
output();
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);
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> hanging_node_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();
dealii::AffineConstraints<double> affine_constraints;
affine_constraints.clear();
{
affine_constraints.reinit(locally_relevant_dofs);
affine_constraints.merge(hanging_node_constraints);
dealii::DoFTools::make_periodicity_constraints<dim, dim>(
periodicity_vector,
affine_constraints);
}
affine_constraints.close();
}
template<int dim>
void MWE<dim>::output()
{
active_fe_index.reinit(triangulation.n_active_cells());
nth_active_fe_index.reinit(triangulation.n_active_cells());
active_fe_index = -1.0;
nth_active_fe_index = -1.0;
for (const auto &cell :
dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
{
active_fe_index(cell->active_cell_index()) =
cell->active_fe_index();
if (cell->at_boundary())
for (const auto &face_id : cell->face_indices())
if (cell->face(face_id)->at_boundary())
{
AssertThrow(
cell->face(face_id)->get_active_fe_indices().size() == 1,
dealii::ExcMessage("Error!"));
nth_active_fe_index(cell->active_cell_index()) =
cell->face(face_id)->nth_active_fe_index(0);
AssertThrow(cell->face(face_id)->fe_index_is_active(
cell->face(face_id)->nth_active_fe_index(0)),
dealii::ExcMessage("Error!"));
cell->face(face_id)->get_fe(
cell->face(face_id)->nth_active_fe_index(0));
}
}
dealii::DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(active_fe_index,
"active_fe_index");
data_out.add_data_vector(nth_active_fe_index,
"nth_active_fe_index");
data_out.build_patches();
data_out.write_vtu_in_parallel(
"triangulation.vtu",
MPI_COMM_WORLD);
}
} // 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_refine_global;
switch (argc)
{
case 1:
flag_refine_global = true;
break;
case 2:
{
const std::string arg(argv[1]);
if (arg == "true")
flag_refine_global = true;
else if (arg == "false")
flag_refine_global = false;
else
AssertThrow(false, dealii::ExcNotImplemented());
}
break;
default:
AssertThrow(false, dealii::ExcNotImplemented());
}
Tests::MWE<2> test(flag_refine_global);
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;
}
