[Libmesh-users] old_solution projection in equation_systems.reinit

[Junchao Z. <jun...@gm...>]

Hello,
  I slightly modified adaptivity_ex2.C ( see the outline below). I find
after the mesh is refined, the old solution is not correctly projected to
the new mesh.
  I print out old solution in each refinement step. At the beginning, time
= 0.025, r_step = 0, old_local_solution->print() prints the solution at
time = 0, which is correct. However, in next refinement (i.e, r_step =
1), old_local_solution->print()'s output has boundary values at time=0.025,
which instead should be values at time=0, I think. The wrong
old_local_solution will affect matrix assembly, where we call
system.old_solution()
to compute u_old and grad_u_old. I tried to debug the code, but it is hard
to get how equation_systems.reinit () works in libmesh.
  Is there something missing in my code? Thank you!

// Outline of main()
system.time = 0;
Create a DirichletBoundary object and add it to equation_system;
equation_systems.init ();

for each time step {
  system.time += dt;
  equation_systems.parameters.set<Real> ("time") = system.time;
  equation_systems.parameters.set<Real> ("dt")   = dt;
  system.reinit_constraints();
  *system.old_local_solution = *system.current_local_solution;
  max_r_steps = 2;
  for (unsigned int r_step=0; r_step<max_r_steps; r_step++) {
     system.old_local_solution->print();
     system.solve();
     system.solution->print();

     if (r_step+1 != max_r_steps) {
       mesh_refinement.uniformly_refine(1);
       // This call reinitializes the \p EquationSystems object for
       // the newly refined mesh.  One of the steps in the
       // reinitialization is projecting the \p solution,
       // \p old_solution, etc... vectors from the old mesh to
       // the current one.
       equation_systems.reinit ();
     }
  }
}

// system.old_solution is used in assemble_cd()
for (unsigned int l=0; l<phi.size(); l++) {
  u_old += phi[l][qp]*system.old_solution (dof_indices[l]);
  grad_u_old.add_scaled (dphi[l][qp], system.old_solution (dof_indices[l]));
}

Fe(i) += ...u_old... + grad_u_old...

--Junchao Zhang