Re: [Libmesh-users] problems with mesh refinement

[Michael P. <pov...@in...>]

This is a function that does the calculation.
Would you like to get also the data?

void QuantumDensity::calculate_convergent_density()
{


  build_k_grid();

  eq = new EquationSystems(*kmesh);

  eq->add_system<LinearImplicitSystem> ("k-integration");
 
  system = &(eq->get_system<LinearImplicitSystem> ("k-integration"));
 
  system->add_variable("u", SECOND);

 
  //!system vector that contains charge density in k space from previous 
iteration
  NumericVector<Number>& old_density = system->add_vector("old density");

 
  eq->init();

 
 
 //----this is to calculate my function and put it into the system solution
  calculate_at_each_k_point();
 
 
 
  calculate_density();

 

  prepare_system_solution();
//---untill here
 
 

  if (opt.k_domain_refinement)
    {
      //----------------------------------------
      //refinement block
      //---------------------------------------
    
    

      old_density = * (system->solution);
     

      MeshRefinement mesh_refinement(*kmesh);

     
      double norm_of_error = opt.relative_accuracy;

    



      for ( ; (norm_of_error >=  opt.relative_accuracy) ;  )
    {//for

   

      if (opt.uniform_refinement)
        mesh_refinement.uniformly_refine(1);
      else
        {
         
          ErrorVector error;
         
         
          KellyErrorEstimator error_estimator;
         
          Tensor2Gen RotM_inv =  transform_matrix.transpose() ;
         
          rotate_mesh(kmesh,  RotM_inv );
         
          error_estimator.estimate_error (*system,error);
         
          rotate_mesh(kmesh, transform_matrix);
         
         
          mesh_refinement.flag_elements_by_error_fraction 
(error,opt.refine_fraction,0.0, 10);
         
         
          mesh_refinement.refine_and_coarsen_elements();
         
        
         
        
        
        

          eq->reinit();

          calculate_at_each_k_point();

          calculate_density();

          prepare_system_solution();


          old_density.add(-1.0, *(system->solution));
         
          old_density.close();

          double x1 = old_density.linfty_norm();

          system->solution->close();

          double x2 = ( system->solution->linfty_norm() );
         
          norm_of_error = x1/x2;


          old_density = *(system->solution);

          std::cout << "k space grid has " << kmesh->n_nodes() << " 
nodes " << flush;
          std::cout <<  "quantum density error " << norm_of_error << 
endl << flush;

         
        }

    }

   
     



    }//end of refinement block


}



John Peterson wrote:
> Michael Povolotskyi writes:
>  > Dear Libmesh developers,
>  > I use Libmesh in order to integrate a function over a 2D mesh.
>  > The idea is to use the adaptive mesh refinement of Libmesh.
>  > 
>  > 
>  > In order to do this, I define a ficticious equation system and assign my 
>  > function as a solution of it.
>  > Then I integrate the function.
>  > After that, the mesh is adaptively refined, the function is recalculated 
>  > on a new mesh and the integration is performed again.
>  > This is done untill the convergence is reached.
>  > 
>  > For the sake of accuracy I use a second order mesh.
>  > 
>  > If I compile a code in the optimized mode, everything works.
>  > However, in the debug mode I get a run-time error:
>  > 
>  > Bad element ID = 49, Bad neighbor ID = 13ERROR: Ancestor Element at 
>  > level 2 found ancestor neighbor at level 1
>  > [0] src/mesh/mesh.C, line 471, compiled Jan  8 2007 at 15:43:37
>
> We would need to look at your code to be sure, but it sounds like
> you are maybe not calling equation_systems.reinit() after each adaptive
> step?  This is causing find_neighbors() to get confused.
>
> -John
>
>  > What can it be?
>  > If you need a more detailed information, e.g. mesh files, I should be 
>  > happy to send it to you.
>  > 
>  > Thank you,
>  > Michael.
>
>   


-- 
------------------------------------------------------------
Michael Povolotskyi, Ph.D.
University of Rome "Tor Vergata"
Department of Electronic Engineering
Viale Politecnico, 1 - 00133 Rome - Italy 

Phone + 39 06 72597781
Fax   + 39 06 2020519
http://www.optolab.uniroma2.it/pages/moshe/moshe.html
-------------------------------------------------------------