subject:"\[deal.II\] Mesh refinement with periodic boundary conditions"

Re: [deal.II] Mesh refinement with periodic boundary conditions

2020-02-28 Thread Andrew Davis

The commit that I am using is: 489c9d198d4ff0654f7d1f69f09404a2613b7954

I will work on a minimal example, but in the mean time here are the 
functions I use to create the sparsity pattern and assemble the system. The 
code works perfectly if either I turn off mesh refinement but keep periodic 
BCs OR keep mesh refinement but instead prescribe in/outflow boundary 
conditions.  I'm solving the advection equation c_{t} + \nabla \cdot (u c) 
= f for c with prescribed velocity u and forcing f using Crank-Nicolson as 
my time integrator.

I'm using MeshWorker to assemble the system. 

The code is working in serial for global refinement.

*Set up sparsity pattern:*

template
void SeaIceModel::SetupDoFs() {
  // set up the DoFs
  dofHandler.distribute_dofs(feSystem);

  // reset the constraints
  constraints.clear();
  
  DoFTools::make_hanging_node_constraints(dofHandler, constraints);

  std::vector::cell_iterator> 
> periodicityVector;
  GridTools::collect_periodic_faces(dofHandler, 0, 1, 0, periodicityVector);
  GridTools::collect_periodic_faces(dofHandler, 2, 3, 1, periodicityVector);
  DoFTools::make_periodicity_constraints 
>(periodicityVector, constraints);

  constraints.close();

  // get the number of degrees of freedom for each block
  std::vector ndofsPerBlock(nblocks);
  DoFTools::count_dofs_per_block(dofHandler, ndofsPerBlock);

  BlockDynamicSparsityPattern dsp(nblocks, nblocks);

  // initalize the sparsity pattern
  for( unsigned int i=0; i
void SeaIceModel::AssembleSystem(double const prevTime, double const 
currTime) {
  // reset the system matrix and rhs to zero
  systemMatrix = 0.0;
  systemRHS = 0.0;

  MeshWorker::IntegrationInfoBox infoBox;

  const unsigned int nGaussPoints = feSystem.degree + 1;
  infoBox.initialize_gauss_quadrature(nGaussPoints, nGaussPoints, 
nGaussPoints);

  infoBox.initialize_update_flags();
  infoBox.add_update_flags_cell(update_quadrature_points | update_values | 
update_gradients | update_JxW_values);
  infoBox.add_update_flags(update_normal_vectors | update_quadrature_points 
| update_values | update_JxW_values, false, true, true, true);

  infoBox.initialize(feSystem, mapping);

  MeshWorker::DoFInfo dofInfo(dofHandler);

  MeshWorker::Assembler::SystemSimple, 
BlockVector >
  assembler;
  assembler.initialize(systemMatrix, systemRHS);

  // objects to integrate over cells and boundaries---these objects have 
overloaded operator() functions for the MeshWorker
  SystemCellIntegration cellIntegration(eqnBlocks, extractors, 
solutionIndexMap, oldSolution, solution, std::tuple(prevTime, currTime, timeStep), nfields);
  SystemBoundaryIntegration boundaryIntegration(eqnBlocks, extractors, 
solutionIndexMap, oldSolution, solution, std::tuple(prevTime, currTime, timeStep), nfields);
  SystemFaceIntegration faceIntegration(eqnBlocks, extractors, 
solutionIndexMap, oldSolution, solution, std::tuple(prevTime, currTime, timeStep), nfields);

  MeshWorker::loop, 
MeshWorker::IntegrationInfoBox 
> (dofHandler.begin_active(), dofHandler.end(), dofInfo, infoBox,
cellIntegration, boundaryIntegration, faceIntegration, assembler);
}


On Friday, February 28, 2020 at 11:19:57 AM UTC-5, Daniel Arndt wrote:
>
> Andrew,
>
> Still hard to say what the problem is. Which commit are you using?
> We have quite a number of tests that make sure that adaptive grid 
> refinement also works with periodic boundary conditions.
> Can you provide a full minimal example that shows the problem you need 
> your workaround for?
>
> How are you assembling the linear system? Are you using MeshWorker, 
> MatrixFree or do you assemble matrices on your own using FEFaceValues?
> Is your code working in serial and for global refinement?
>
> Best,
> Daniel
>
>
> Am Fr., 28. Feb. 2020 um 08:29 Uhr schrieb Andrew Davis  >:
>
>> I just say a slight problem and changed this line the declaration of 
>> periodicityVector to:
>>
>> std::vector> Triangulation::cell_iterator> > periodicityVector;
>>
>> But I'm getting the same behavior.
>>
>> On Friday, February 28, 2020 at 8:08:19 AM UTC-5, Andrew Davis wrote:
>>>
>>> Yes, here is the code I used to tell the triangulation about the 
>>> periodic boundaries: 
>>>
>>>   // create a mesh on [0,1] with colorized boundaries
>>>   GridGenerator::hyper_cube(triangulation, 0.0, 1.0, true);
>>>
>>>   // collect the periodic boundary pairs
>>>   std::vector>> parallel::distributed::Triangulation::cell_iterator> > 
>>> periodicityVector;
>>>   GridTools::collect_periodic_faces(triangulation, 0, 1, 0, 
>>> periodicityVector);
>>>   GridTools::collect_periodic_faces(triangulation, 2, 3, 1, 
>>> periodicityVector);
>>>
>>>   triangulation.add_periodicity(periodicityVector);
>>>
>>>   // refine one less then maximum refinement level globally
>>>   triangulation.refine_global(maxGridLevel-1);
>>>
>>> Perhaps I'm doing it wrong?
>>>
>>> On Thursday, February 27, 2020 at 10:47:09 PM UTC-5, Daniel Arndt wrote:

 Andrew,

 Did you tell the triangulation that it

Re: [deal.II] Mesh refinement with periodic boundary conditions

2020-02-28 Thread Daniel Arndt

Andrew,

Still hard to say what the problem is. Which commit are you using?
We have quite a number of tests that make sure that adaptive grid
refinement also works with periodic boundary conditions.
Can you provide a full minimal example that shows the problem you need your
workaround for?

How are you assembling the linear system? Are you using MeshWorker,
MatrixFree or do you assemble matrices on your own using FEFaceValues?
Is your code working in serial and for global refinement?

Best,
Daniel


Am Fr., 28. Feb. 2020 um 08:29 Uhr schrieb Andrew Davis <
andrew@gmail.com>:

> I just say a slight problem and changed this line the declaration of
> periodicityVector to:
>
> std::vector Triangulation::cell_iterator>
> > periodicityVector;
>
> But I'm getting the same behavior.
>
> On Friday, February 28, 2020 at 8:08:19 AM UTC-5, Andrew Davis wrote:
>>
>> Yes, here is the code I used to tell the triangulation about the periodic
>> boundaries:
>>
>>   // create a mesh on [0,1] with colorized boundaries
>>   GridGenerator::hyper_cube(triangulation, 0.0, 1.0, true);
>>
>>   // collect the periodic boundary pairs
>>   std::vector> parallel::distributed::Triangulation::cell_iterator> >
>> periodicityVector;
>>   GridTools::collect_periodic_faces(triangulation, 0, 1, 0,
>> periodicityVector);
>>   GridTools::collect_periodic_faces(triangulation, 2, 3, 1,
>> periodicityVector);
>>
>>   triangulation.add_periodicity(periodicityVector);
>>
>>   // refine one less then maximum refinement level globally
>>   triangulation.refine_global(maxGridLevel-1);
>>
>> Perhaps I'm doing it wrong?
>>
>> On Thursday, February 27, 2020 at 10:47:09 PM UTC-5, Daniel Arndt wrote:
>>>
>>> Andrew,
>>>
>>> Did you tell the triangulation that it should take periodic boundaries
>>> into account? My best bet (with the information given) would be that you
>>> didn't call parallel::distributed::Triangulation::add_periodicity()
>>> 
>>> as described in
>>> https://www.dealii.org/current/doxygen/deal.II/step_45.html.
>>>
>>> Best,
>>> Daniel
>>>
>>> Am Do., 27. Feb. 2020 um 17:45 Uhr schrieb Andrew Davis <
>>> andr...@gmail.com>:
>>>
 HI all,

 I'm trying to implement a time-dependent convection equation using DG
 elements with periodic boundary conditions on and adaptive mesh.

 However, when I try to adapt the mesh using periodic using this code:

   // estimate the error in each cell
   Vector estimatedErrorPerCell(triangulation.n_active_cells());
   //KellyErrorEstimator::estimate(dofHandler,
 QGauss(feSystem.degree + 1), {}, solution.block(0),
 estimatedErrorPerCell);
   KellyErrorEstimator::estimate(dofHandler,
 QGauss(feSystem.degree + 1), {}, solution, estimatedErrorPerCell);

   // refine cells with the largest estimated error (80 per cent of the
 error) and coarsens those cells with the smallest error (10 per cent of the
 error)
   GridRefinement::refine_and_coarsen_fixed_fraction(triangulation,
 estimatedErrorPerCell, 0.8, 0.1);

   // to prevent the decrease in time step size limit the maximal
 refinement depth of the meshlimit the maximal refinement depth of the mesh
   if( triangulation.n_levels()>maxGridLevel ) {
 for( auto& cell :
 triangulation.active_cell_iterators_on_level(maxGridLevel) ) {
 cell->clear_refine_flag(); }
   }

   // store previous and current solution on the coarser mesh
   std::vector > tempSolution(2);
   tempSolution[0] = oldSolution; tempSolution[1] = solution;

   // transfer the solution vectors from the old mesh to the new one
   SolutionTransfer > transfer(dofHandler);

   // prepare for the refinement---we have to prepare the transfer as
 well
   *triangulation.prepare_**coarsening_and_refinement();* // for some
 reason, this doesn't tell the cells across a periodic boundary that they
 need to be refined
   transfer.prepare_for_coarsening_and_refinement(tempSolution);

   // do the refinement and reset the DoFs
   triangulation.execute_coarsening_and_refinement();
   SetupSystem();

   std::vector > temp(2);
   ReinitializeVector(temp[0]); ReinitializeVector(temp[1]);
   transfer.interpolate(tempSolution, temp);

   oldSolution = temp[0]; solution = temp[1];

   constraints.distribute(oldSolution);
   constraints.distribute(solution);

 *I get this error:*

 
 An error occurred in line <992> of file
  in function
 void
 dealii::{anonymous}::update_periodic_face_map_recursively(const typename
 dealii::Triangulation::cell_iterator&, const typename
 dealii::Triangulation::cell_iterator&, unsigned int,
 unsigned int, const std::bitset<3>&,

Re: [deal.II] Mesh refinement with periodic boundary conditions

2020-02-28 Thread Andrew Davis

I just say a slight problem and changed this line the declaration of 
periodicityVector to:

std::vector::cell_iterator> 
> periodicityVector;

But I'm getting the same behavior.

On Friday, February 28, 2020 at 8:08:19 AM UTC-5, Andrew Davis wrote:
>
> Yes, here is the code I used to tell the triangulation about the periodic 
> boundaries: 
>
>   // create a mesh on [0,1] with colorized boundaries
>   GridGenerator::hyper_cube(triangulation, 0.0, 1.0, true);
>
>   // collect the periodic boundary pairs
>   std::vector parallel::distributed::Triangulation::cell_iterator> > 
> periodicityVector;
>   GridTools::collect_periodic_faces(triangulation, 0, 1, 0, 
> periodicityVector);
>   GridTools::collect_periodic_faces(triangulation, 2, 3, 1, 
> periodicityVector);
>
>   triangulation.add_periodicity(periodicityVector);
>
>   // refine one less then maximum refinement level globally
>   triangulation.refine_global(maxGridLevel-1);
>
> Perhaps I'm doing it wrong?
>
> On Thursday, February 27, 2020 at 10:47:09 PM UTC-5, Daniel Arndt wrote:
>>
>> Andrew,
>>
>> Did you tell the triangulation that it should take periodic boundaries 
>> into account? My best bet (with the information given) would be that you 
>> didn't call parallel::distributed::Triangulation::add_periodicity() 
>> 
>> as described in 
>> https://www.dealii.org/current/doxygen/deal.II/step_45.html.
>>
>> Best,
>> Daniel
>>
>> Am Do., 27. Feb. 2020 um 17:45 Uhr schrieb Andrew Davis <
>> andr...@gmail.com>:
>>
>>> HI all,
>>>
>>> I'm trying to implement a time-dependent convection equation using DG 
>>> elements with periodic boundary conditions on and adaptive mesh.
>>>
>>> However, when I try to adapt the mesh using periodic using this code:
>>>
>>>   // estimate the error in each cell
>>>   Vector estimatedErrorPerCell(triangulation.n_active_cells());
>>>   //KellyErrorEstimator::estimate(dofHandler, 
>>> QGauss(feSystem.degree + 1), {}, solution.block(0), 
>>> estimatedErrorPerCell);
>>>   KellyErrorEstimator::estimate(dofHandler, 
>>> QGauss(feSystem.degree + 1), {}, solution, estimatedErrorPerCell);
>>>
>>>   // refine cells with the largest estimated error (80 per cent of the 
>>> error) and coarsens those cells with the smallest error (10 per cent of the 
>>> error)
>>>   GridRefinement::refine_and_coarsen_fixed_fraction(triangulation, 
>>> estimatedErrorPerCell, 0.8, 0.1);
>>>
>>>   // to prevent the decrease in time step size limit the maximal 
>>> refinement depth of the meshlimit the maximal refinement depth of the mesh
>>>   if( triangulation.n_levels()>maxGridLevel ) {
>>> for( auto& cell : 
>>> triangulation.active_cell_iterators_on_level(maxGridLevel) ) { 
>>> cell->clear_refine_flag(); }
>>>   }
>>>
>>>   // store previous and current solution on the coarser mesh
>>>   std::vector > tempSolution(2);
>>>   tempSolution[0] = oldSolution; tempSolution[1] = solution;
>>>
>>>   // transfer the solution vectors from the old mesh to the new one
>>>   SolutionTransfer > transfer(dofHandler);
>>>
>>>   // prepare for the refinement---we have to prepare the transfer as well
>>>   *triangulation.prepare_**coarsening_and_refinement();* // for some 
>>> reason, this doesn't tell the cells across a periodic boundary that they 
>>> need to be refined
>>>   transfer.prepare_for_coarsening_and_refinement(tempSolution);
>>>
>>>   // do the refinement and reset the DoFs
>>>   triangulation.execute_coarsening_and_refinement();
>>>   SetupSystem();
>>>
>>>   std::vector > temp(2);
>>>   ReinitializeVector(temp[0]); ReinitializeVector(temp[1]);
>>>   transfer.interpolate(tempSolution, temp);
>>>
>>>   oldSolution = temp[0]; solution = temp[1];
>>>
>>>   constraints.distribute(oldSolution);
>>>   constraints.distribute(solution);
>>>
>>> *I get this error:*
>>>
>>> 
>>> An error occurred in line <992> of file 
>>>  in function
>>> void dealii::{anonymous}::update_periodic_face_map_recursively(const 
>>> typename dealii::Triangulation::cell_iterator&, const 
>>> typename dealii::Triangulation::cell_iterator&, unsigned 
>>> int, unsigned int, const std::bitset<3>&, std::map>> dealii::Triangulation::cell_iterator, unsigned int>, 
>>> std::pair>> spacedim>::cell_iterator, unsigned int>, std::bitset<3> > >&) [with int dim 
>>> = 2; int spacedim = 2; typename dealii::Triangulation>> spacedim>::cell_iterator = dealii::TriaIterator 
>>> >]
>>> The violated condition was:
>>> dim == 1 || std::abs(cell_1->level() - cell_2->level()) < 2
>>> Additional information:
>>> This exception -- which is used in many places in the library -- 
>>> usually indicates that some condition which the author of the code thought 
>>> must be satisfied at a certain point in an algorithm, is not fulfilled. An 
>>> example would be that the first part of an algorithm sorts elements of

Re: [deal.II] Mesh refinement with periodic boundary conditions

2020-02-28 Thread Andrew Davis

Yes, here is the code I used to tell the triangulation about the periodic 
boundaries: 

  // create a mesh on [0,1] with colorized boundaries
  GridGenerator::hyper_cube(triangulation, 0.0, 1.0, true);

  // collect the periodic boundary pairs
  std::vector::cell_iterator> > 
periodicityVector;
  GridTools::collect_periodic_faces(triangulation, 0, 1, 0, 
periodicityVector);
  GridTools::collect_periodic_faces(triangulation, 2, 3, 1, 
periodicityVector);

  triangulation.add_periodicity(periodicityVector);

  // refine one less then maximum refinement level globally
  triangulation.refine_global(maxGridLevel-1);

Perhaps I'm doing it wrong?

On Thursday, February 27, 2020 at 10:47:09 PM UTC-5, Daniel Arndt wrote:
>
> Andrew,
>
> Did you tell the triangulation that it should take periodic boundaries 
> into account? My best bet (with the information given) would be that you 
> didn't call parallel::distributed::Triangulation::add_periodicity() 
> 
> as described in 
> https://www.dealii.org/current/doxygen/deal.II/step_45.html.
>
> Best,
> Daniel
>
> Am Do., 27. Feb. 2020 um 17:45 Uhr schrieb Andrew Davis  >:
>
>> HI all,
>>
>> I'm trying to implement a time-dependent convection equation using DG 
>> elements with periodic boundary conditions on and adaptive mesh.
>>
>> However, when I try to adapt the mesh using periodic using this code:
>>
>>   // estimate the error in each cell
>>   Vector estimatedErrorPerCell(triangulation.n_active_cells());
>>   //KellyErrorEstimator::estimate(dofHandler, 
>> QGauss(feSystem.degree + 1), {}, solution.block(0), 
>> estimatedErrorPerCell);
>>   KellyErrorEstimator::estimate(dofHandler, 
>> QGauss(feSystem.degree + 1), {}, solution, estimatedErrorPerCell);
>>
>>   // refine cells with the largest estimated error (80 per cent of the 
>> error) and coarsens those cells with the smallest error (10 per cent of the 
>> error)
>>   GridRefinement::refine_and_coarsen_fixed_fraction(triangulation, 
>> estimatedErrorPerCell, 0.8, 0.1);
>>
>>   // to prevent the decrease in time step size limit the maximal 
>> refinement depth of the meshlimit the maximal refinement depth of the mesh
>>   if( triangulation.n_levels()>maxGridLevel ) {
>> for( auto& cell : 
>> triangulation.active_cell_iterators_on_level(maxGridLevel) ) { 
>> cell->clear_refine_flag(); }
>>   }
>>
>>   // store previous and current solution on the coarser mesh
>>   std::vector > tempSolution(2);
>>   tempSolution[0] = oldSolution; tempSolution[1] = solution;
>>
>>   // transfer the solution vectors from the old mesh to the new one
>>   SolutionTransfer > transfer(dofHandler);
>>
>>   // prepare for the refinement---we have to prepare the transfer as well
>>   *triangulation.prepare_**coarsening_and_refinement();* // for some 
>> reason, this doesn't tell the cells across a periodic boundary that they 
>> need to be refined
>>   transfer.prepare_for_coarsening_and_refinement(tempSolution);
>>
>>   // do the refinement and reset the DoFs
>>   triangulation.execute_coarsening_and_refinement();
>>   SetupSystem();
>>
>>   std::vector > temp(2);
>>   ReinitializeVector(temp[0]); ReinitializeVector(temp[1]);
>>   transfer.interpolate(tempSolution, temp);
>>
>>   oldSolution = temp[0]; solution = temp[1];
>>
>>   constraints.distribute(oldSolution);
>>   constraints.distribute(solution);
>>
>> *I get this error:*
>>
>> 
>> An error occurred in line <992> of file 
>>  in function
>> void dealii::{anonymous}::update_periodic_face_map_recursively(const 
>> typename dealii::Triangulation::cell_iterator&, const 
>> typename dealii::Triangulation::cell_iterator&, unsigned 
>> int, unsigned int, const std::bitset<3>&, std::map> dealii::Triangulation::cell_iterator, unsigned int>, 
>> std::pair> spacedim>::cell_iterator, unsigned int>, std::bitset<3> > >&) [with int dim 
>> = 2; int spacedim = 2; typename dealii::Triangulation> spacedim>::cell_iterator = dealii::TriaIterator 
>> >]
>> The violated condition was:
>> dim == 1 || std::abs(cell_1->level() - cell_2->level()) < 2
>> Additional information:
>> This exception -- which is used in many places in the library -- 
>> usually indicates that some condition which the author of the code thought 
>> must be satisfied at a certain point in an algorithm, is not fulfilled. An 
>> example would be that the first part of an algorithm sorts elements of an 
>> array in ascending order, and a second part of the algorithm later 
>> encounters an element that is not larger than the previous one.
>>
>> There is usually not very much you can do if you encounter such an 
>> exception since it indicates an error in deal.II, not in your own program. 
>> Try to come up with the smallest possible program that still demonstrates 
>> the error and contact the deal.II mailing lists with it to obtain

Re: [deal.II] Mesh refinement with periodic boundary conditions

2020-02-27 Thread Daniel Arndt

Andrew,

Did you tell the triangulation that it should take periodic boundaries into
account? My best bet (with the information given) would be that you didn't
call parallel::distributed::Triangulation::add_periodicity()

as described in https://www.dealii.org/current/doxygen/deal.II/step_45.html.

Best,
Daniel

Am Do., 27. Feb. 2020 um 17:45 Uhr schrieb Andrew Davis <
andrew@gmail.com>:

> HI all,
>
> I'm trying to implement a time-dependent convection equation using DG
> elements with periodic boundary conditions on and adaptive mesh.
>
> However, when I try to adapt the mesh using periodic using this code:
>
>   // estimate the error in each cell
>   Vector estimatedErrorPerCell(triangulation.n_active_cells());
>   //KellyErrorEstimator::estimate(dofHandler,
> QGauss(feSystem.degree + 1), {}, solution.block(0),
> estimatedErrorPerCell);
>   KellyErrorEstimator::estimate(dofHandler,
> QGauss(feSystem.degree + 1), {}, solution, estimatedErrorPerCell);
>
>   // refine cells with the largest estimated error (80 per cent of the
> error) and coarsens those cells with the smallest error (10 per cent of the
> error)
>   GridRefinement::refine_and_coarsen_fixed_fraction(triangulation,
> estimatedErrorPerCell, 0.8, 0.1);
>
>   // to prevent the decrease in time step size limit the maximal
> refinement depth of the meshlimit the maximal refinement depth of the mesh
>   if( triangulation.n_levels()>maxGridLevel ) {
> for( auto& cell :
> triangulation.active_cell_iterators_on_level(maxGridLevel) ) {
> cell->clear_refine_flag(); }
>   }
>
>   // store previous and current solution on the coarser mesh
>   std::vector > tempSolution(2);
>   tempSolution[0] = oldSolution; tempSolution[1] = solution;
>
>   // transfer the solution vectors from the old mesh to the new one
>   SolutionTransfer > transfer(dofHandler);
>
>   // prepare for the refinement---we have to prepare the transfer as well
>   *triangulation.prepare_**coarsening_and_refinement();* // for some
> reason, this doesn't tell the cells across a periodic boundary that they
> need to be refined
>   transfer.prepare_for_coarsening_and_refinement(tempSolution);
>
>   // do the refinement and reset the DoFs
>   triangulation.execute_coarsening_and_refinement();
>   SetupSystem();
>
>   std::vector > temp(2);
>   ReinitializeVector(temp[0]); ReinitializeVector(temp[1]);
>   transfer.interpolate(tempSolution, temp);
>
>   oldSolution = temp[0]; solution = temp[1];
>
>   constraints.distribute(oldSolution);
>   constraints.distribute(solution);
>
> *I get this error:*
>
> 
> An error occurred in line <992> of file
>  in function
> void dealii::{anonymous}::update_periodic_face_map_recursively(const
> typename dealii::Triangulation::cell_iterator&, const
> typename dealii::Triangulation::cell_iterator&, unsigned
> int, unsigned int, const std::bitset<3>&, std::map dealii::Triangulation::cell_iterator, unsigned int>,
> std::pair spacedim>::cell_iterator, unsigned int>, std::bitset<3> > >&) [with int dim
> = 2; int spacedim = 2; typename dealii::Triangulation spacedim>::cell_iterator = dealii::TriaIterator
> >]
> The violated condition was:
> dim == 1 || std::abs(cell_1->level() - cell_2->level()) < 2
> Additional information:
> This exception -- which is used in many places in the library --
> usually indicates that some condition which the author of the code thought
> must be satisfied at a certain point in an algorithm, is not fulfilled. An
> example would be that the first part of an algorithm sorts elements of an
> array in ascending order, and a second part of the algorithm later
> encounters an element that is not larger than the previous one.
>
> There is usually not very much you can do if you encounter such an
> exception since it indicates an error in deal.II, not in your own program.
> Try to come up with the smallest possible program that still demonstrates
> the error and contact the deal.II mailing lists with it to obtain help.
>
> For some reason the refinement does not seem to know that the boundaries
> are periodic and that it should refine the cells on the other side of the
> domain. I tried changing the refinement step so that instead of just
> calling triangulation.prepare_coarsening_and_refinement();, I created my
> own function that mimics the distributed triangulation (calling the
> distributed version of enforce_mesh_balance_over_periodic_boundaries)
> that does the following:
>
> while( true ) {
> triangulation.prepare_coarsening_and_refinement();
> triangulation.update_periodic_face_map();
> bool changed =
> enforce_mesh_balance_over_periodic_boundaries(triangulation);
> if( !changed ) { break; }
>   }
>
> This fixed the error but I'm getting the strange behavior---it looks like
> the cell face fluxes are no longer being

[deal.II] Mesh refinement with periodic boundary conditions

2020-02-27 Thread Andrew Davis


HI all,

I'm trying to implement a time-dependent convection equation using DG 
elements with periodic boundary conditions on and adaptive mesh.

However, when I try to adapt the mesh using periodic using this code:

  // estimate the error in each cell
  Vector estimatedErrorPerCell(triangulation.n_active_cells());
  //KellyErrorEstimator::estimate(dofHandler, 
QGauss(feSystem.degree + 1), {}, solution.block(0), 
estimatedErrorPerCell);
  KellyErrorEstimator::estimate(dofHandler, 
QGauss(feSystem.degree + 1), {}, solution, estimatedErrorPerCell);

  // refine cells with the largest estimated error (80 per cent of the 
error) and coarsens those cells with the smallest error (10 per cent of the 
error)
  GridRefinement::refine_and_coarsen_fixed_fraction(triangulation, 
estimatedErrorPerCell, 0.8, 0.1);

  // to prevent the decrease in time step size limit the maximal refinement 
depth of the meshlimit the maximal refinement depth of the mesh
  if( triangulation.n_levels()>maxGridLevel ) {
for( auto& cell : 
triangulation.active_cell_iterators_on_level(maxGridLevel) ) { 
cell->clear_refine_flag(); }
  }

  // store previous and current solution on the coarser mesh
  std::vector > tempSolution(2);
  tempSolution[0] = oldSolution; tempSolution[1] = solution;

  // transfer the solution vectors from the old mesh to the new one
  SolutionTransfer > transfer(dofHandler);

  // prepare for the refinement---we have to prepare the transfer as well
  *triangulation.prepare_**coarsening_and_refinement();* // for some 
reason, this doesn't tell the cells across a periodic boundary that they 
need to be refined
  transfer.prepare_for_coarsening_and_refinement(tempSolution);

  // do the refinement and reset the DoFs
  triangulation.execute_coarsening_and_refinement();
  SetupSystem();

  std::vector > temp(2);
  ReinitializeVector(temp[0]); ReinitializeVector(temp[1]);
  transfer.interpolate(tempSolution, temp);

  oldSolution = temp[0]; solution = temp[1];

  constraints.distribute(oldSolution);
  constraints.distribute(solution);

*I get this error:*


An error occurred in line <992> of file 
 in function
void dealii::{anonymous}::update_periodic_face_map_recursively(const 
typename dealii::Triangulation::cell_iterator&, const 
typename dealii::Triangulation::cell_iterator&, unsigned 
int, unsigned int, const std::bitset<3>&, std::map::cell_iterator, unsigned int>, 
std::pair::cell_iterator, unsigned int>, std::bitset<3> > >&) [with int dim 
= 2; int spacedim = 2; typename dealii::Triangulation::cell_iterator = dealii::TriaIterator 
>]
The violated condition was:
dim == 1 || std::abs(cell_1->level() - cell_2->level()) < 2
Additional information:
This exception -- which is used in many places in the library -- 
usually indicates that some condition which the author of the code thought 
must be satisfied at a certain point in an algorithm, is not fulfilled. An 
example would be that the first part of an algorithm sorts elements of an 
array in ascending order, and a second part of the algorithm later 
encounters an element that is not larger than the previous one.

There is usually not very much you can do if you encounter such an 
exception since it indicates an error in deal.II, not in your own program. 
Try to come up with the smallest possible program that still demonstrates 
the error and contact the deal.II mailing lists with it to obtain help.

For some reason the refinement does not seem to know that the boundaries 
are periodic and that it should refine the cells on the other side of the 
domain. I tried changing the refinement step so that instead of just 
calling triangulation.prepare_coarsening_and_refinement();, I created my 
own function that mimics the distributed triangulation (calling the 
distributed version of enforce_mesh_balance_over_periodic_boundaries) that 
does the following:

while( true ) {
triangulation.prepare_coarsening_and_refinement();
triangulation.update_periodic_face_map();
bool changed = 
enforce_mesh_balance_over_periodic_boundaries(triangulation);
if( !changed ) { break; }
  }

This fixed the error but I'm getting the strange behavior---it looks like 
the cell face fluxes are no longer being computed correctly

Does anyone have any idea about what the issue could be? I know I left out 
a lot of code, but hopefully I included the important bits.

Thanks,
Andy
ReplyForward

 
 


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

Re: [deal.II] Mesh refinement with periodic boundary conditions

Re: [deal.II] Mesh refinement with periodic boundary conditions

Re: [deal.II] Mesh refinement with periodic boundary conditions

Re: [deal.II] Mesh refinement with periodic boundary conditions

Re: [deal.II] Mesh refinement with periodic boundary conditions

[deal.II] Mesh refinement with periodic boundary conditions

6 matches

Site Navigation

Mail list logo

Footer information