Hey!
I'm trying to get deal.II and mhdeal (an MHD implementation with deal.II on
GitHub) to work with adaptive mesh refinement. The problem I'm running into
is adaptive mesh refinement and NAN errors after a couple of iterations.
Below I will copy/paste the relevant (I think) parts of the code that I'm
trying to get to work. I'm new to deal.II so a lot of my development is
looking at what others have done but I'm not sure what is relevant and what
isn't. Hopefully, you can help me diagnose where the problem lies.
Step 1: Set up the triangulation
*GridGenerator::subdivided_hyper_rectangle(triangulation,
parameters.refinements, parameters.corner_a, parameters.corner_b, true);
std::vector<dealii::GridTools::PeriodicFacePair<
dealii::TriaIterator<dealii::CellAccessor<DIMENSION> > > > matched_pairs;
for (std::vector<std::array<int, 3> >::const_iterator it =
parameters.periodic_boundaries.begin(); it !=
parameters.periodic_boundaries.end(); it++)
dealii::GridTools::collect_periodic_faces(triangulation, (*it)[0],
(*it)[1], (*it)[2], matched_pairs);
triangulation.add_periodicity(matched_pairs);*
Step 2: Calculate error and flag cells for refinement/coarsening
*Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
KellyErrorEstimator<dim>::estimate(dof_handler,
QGauss<dim - 1>(this->parameters.quadrature_order + 1),
{}, solution,
estimated_error_per_cell);
GridRefinement::refine_and_coarsen_fixed_fraction(triangulation,
estimated_error_per_cell, this->parameters.refine_threshold,
this->parameters.coarsen_threshold);*
Step 3: Enforce the same level of refinement on the boundary
*for (typename Triangulation<dim>::active_cell_iterator cell =
triangulation.begin_active(); cell != triangulation.end(); ++cell) {
for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
++face) { if (cell->has_periodic_neighbor(face)) { if
(cell->refine_flag_set()) {
cell->periodic_neighbor(face)->clear_coarsen_flag();
cell->periodic_neighbor(face)->set_refine_flag(); } else if
(cell->periodic_neighbor(face)->refine_flag_set()) {
cell->clear_coarsen_flag(); cell->set_refine_flag(); }
else if (cell->coarsen_flag_set()) {
cell->periodic_neighbor(face)->clear_refine_flag();
cell->periodic_neighbor(face)->set_coarsen_flag(); } else if
(cell->periodic_neighbor(face)->coarsen_flag_set()) {
cell->clear_refine_flag(); cell->set_coarsen_flag(); }
} } } triangulation.prepare_coarsening_and_refinement();*
Step 4: Set constraints
*dof_handler.clear(); dof_handler.distribute_dofs(fe);
DoFTools::extract_locally_relevant_dofs(dof_handler,
locally_relevant_dofs); locally_owned_dofs =
dof_handler.locally_owned_dofs(); constraints.clear();
constraints.reinit(locally_relevant_dofs);
DoFTools::make_hanging_node_constraints(dof_handler, constraints); if
(this->parameters.periodic_boundaries.size() > 0) {
std::vector<GridTools::PeriodicFacePair<typename
DoFHandler<dim>::cell_iterator> > matched_pairs; for
(std::vector<std::array<int, 3> >::const_iterator it =
parameters.periodic_boundaries.begin(); it !=
parameters.periodic_boundaries.end(); it++)
GridTools::collect_periodic_faces(dof_handler, (*it)[0], (*it)[1],
(*it)[2], matched_pairs);
DoFTools::make_periodicity_constraints<DoFHandler<dim>>(matched_pairs,
constraints); } DynamicSparsityPattern dsp(locally_relevant_dofs);
DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false);
constraints.close();*
--
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/91061305-373f-403c-8e54-76cd9f76f9a0n%40googlegroups.com.
