uniform_refinement_estimator.c

一个用来实现偏微分方程中网格的计算库

      projected_solutions[i] = NumericVector<Number>::build().release();
      projected_solutions[i]->init(system.solution->size(), system.solution->size());
      system.solution->localize(*projected_solutions[i],
                                system.get_dof_map().get_send_list());
    }

  // Get the uniformly refined solution.

  if (_es)
    es.solve();
  else
    system_list[0]->solve();
  
  // Get the error in the uniformly refined solution(s).

  for (unsigned int i=0; i != system_list.size(); ++i)
    {
      System &system = *system_list[i];

      unsigned int n_vars = system.n_vars();

      DofMap &dof_map = system.get_dof_map();

      std::vector<float> &_component_scale =
        (*component_scales)[&system];

      NumericVector<Number> *projected_solution = projected_solutions[i];

      // Loop over all the variables in the system
      for (unsigned int var=0; var<n_vars; var++)
        {
          Real var_scale = 1.0;

          // Possibly skip this variable
          if (!_component_scale.empty())
            {
	      if (_component_scale[var] == 0.0)
                continue;
              else
                var_scale = _component_scale[var];
            }

          // Get the error vector to fill for this system and variable
          ErrorVector *err_vec = error_per_cell;
          if (!err_vec)
            {
              libmesh_assert(errors_per_cell);
	      err_vec =
		(*errors_per_cell)[std::make_pair(&system,var)];
            }

          // The type of finite element to use for this variable
          const FEType& fe_type = dof_map.variable_type (var);
      
          // Finite element object for each fine element
          AutoPtr<FEBase> fe (FEBase::build (dim, fe_type));

          // Build and attach an appropriate quadrature rule
          AutoPtr<QBase> qrule = fe_type.default_quadrature_rule(dim);
          fe->attach_quadrature_rule (qrule.get());
      
          const std::vector<Real>&  JxW = fe->get_JxW();
          const std::vector<std::vector<Real> >& phi = fe->get_phi();
          const std::vector<std::vector<RealGradient> >& dphi =
            fe->get_dphi();
#ifdef ENABLE_SECOND_DERIVATIVES
          const std::vector<std::vector<RealTensor> >& d2phi =
            fe->get_d2phi();
#endif

          // The global DOF indices for the fine element
          std::vector<unsigned int> dof_indices;
      
          // Iterate over all the active elements in the fine mesh
          // that live on this processor.
          MeshBase::const_element_iterator       elem_it  = mesh.active_local_elements_begin();
          const MeshBase::const_element_iterator elem_end = mesh.active_local_elements_end(); 

          for (; elem_it != elem_end; ++elem_it)
	    {
	      // e is necessarily an active element on the local processor
	      const Elem* elem = *elem_it;

              // Find the element id for the corresponding coarse grid element
              const Elem* coarse = elem;
              unsigned int e_id = coarse->id();
              while (e_id >= max_coarse_elem_id)
                {
                  libmesh_assert (coarse->parent());
                  coarse = coarse->parent();
                  e_id = coarse->id();
                }
          
              double L2normsq = 0., H1seminormsq = 0., H2seminormsq = 0.;

              // reinitialize the element-specific data
              // for the current element
              fe->reinit (elem);

              // Get the local to global degree of freedom maps
              dof_map.dof_indices (elem, dof_indices, var);

              // The number of quadrature points
              const unsigned int n_qp = qrule->n_points();

              // The number of shape functions
              const unsigned int n_sf = dof_indices.size();

              //
              // Begin the loop over the Quadrature points.
              //
              for (unsigned int qp=0; qp<n_qp; qp++)
                {
                  Number u_fine = 0., u_coarse = 0.;

                  Gradient grad_u_fine, grad_u_coarse;
#ifdef ENABLE_SECOND_DERIVATIVES
                  Tensor grad2_u_fine, grad2_u_coarse;
#endif

                  // Compute solution values at the current
                  // quadrature point.  This reqiures a sum
                  // over all the shape functions evaluated
                  // at the quadrature point.
                  for (unsigned int i=0; i<n_sf; i++)
                    {
                      u_fine            += phi[i][qp]*system.current_solution (dof_indices[i]);
                      u_coarse          += phi[i][qp]*(*projected_solution) (dof_indices[i]);
                      grad_u_fine       += dphi[i][qp]*system.current_solution (dof_indices[i]);
                      grad_u_coarse     += dphi[i][qp]*(*projected_solution) (dof_indices[i]);
#ifdef ENABLE_SECOND_DERIVATIVES
                      grad2_u_fine      += d2phi[i][qp]*system.current_solution (dof_indices[i]);
                      grad2_u_coarse    += d2phi[i][qp]*(*projected_solution) (dof_indices[i]);
#endif
                    }

                  // Compute the value of the error at this quadrature point
                  const Number val_error = u_fine - u_coarse;

                  // Add the squares of the error to each contribution
		  L2normsq += JxW[qp] * var_scale *
                    libmesh_norm(val_error);
                  libmesh_assert (L2normsq     >= 0.);


                  // Compute the value of the error in the gradient at this
                  // quadrature point
                  if (_sobolev_order > 0)
                    {
                      Gradient grad_error = grad_u_fine - grad_u_coarse;

                      H1seminormsq += JxW[qp] * var_scale *
                        grad_error.size_sq();
                      libmesh_assert (H1seminormsq >= 0.);
                    }

#ifdef ENABLE_SECOND_DERIVATIVES
                  // Compute the value of the error in the hessian at this
                  // quadrature point
                  if (_sobolev_order > 1)
                    {
                      Tensor grad2_error = grad2_u_fine - grad2_u_coarse;

		      H2seminormsq += JxW[qp] * var_scale *
                        grad2_error.size_sq();
                      libmesh_assert (H2seminormsq >= 0.);
                    }
#endif
                } // end qp loop

              (*err_vec)[e_id] += L2normsq;
              if (_sobolev_order > 0)
                (*err_vec)[e_id] += H1seminormsq;
              if (_sobolev_order > 1)
                (*err_vec)[e_id] += H2seminormsq;
            } // End loop over active local elements
        } // End loop over variables

      // Don't bother projecting the solution; we'll restore from backup
      // after coarsening
      system.project_solution_on_reinit() = false;
    }


  // Uniformly coarsen the mesh, without projecting the solution
  libmesh_assert (number_h_refinements > 0 || number_p_refinements > 0);

  for (unsigned int i = 0; i != number_h_refinements; ++i)
    {
      mesh_refinement.uniformly_coarsen(1);
      // FIXME - should the reinits here be necessary? - RHS
      es.reinit();
    }
      
  for (unsigned int i = 0; i != number_p_refinements; ++i)
    {
      mesh_refinement.uniformly_p_coarsen(1);
      es.reinit();
    }

  // We should be back where we started
  libmesh_assert(n_coarse_elem == mesh.n_elem());

  // Each processor has now computed the error contribuions
  // for its local elements.  We need to sum the vector
  // and then take the square-root of each component.  Note
  // that we only need to sum if we are running on multiple
  // processors, and we only need to take the square-root
  // if the value is nonzero.  There will in general be many
  // zeros for the inactive elements.

  if (error_per_cell)
    {
      // First sum the vector of estimated error values
      this->reduce_error(*error_per_cell);

      // Compute the square-root of each component.
      START_LOG("std::sqrt()", "UniformRefinementEstimator");
      for (unsigned int i=0; i<error_per_cell->size(); i++)
        if ((*error_per_cell)[i] != 0.)
          (*error_per_cell)[i] = std::sqrt((*error_per_cell)[i]);
      STOP_LOG("std::sqrt()", "UniformRefinementEstimator");
    }
  else
    {
      for (ErrorMap::iterator i = errors_per_cell->begin();
           i != errors_per_cell->end(); ++i)
        {
          ErrorVector *e = i->second;
          // First sum the vector of estimated error values
          this->reduce_error(*e);

          // Compute the square-root of each component.
          START_LOG("std::sqrt()", "UniformRefinementEstimator");
          for (unsigned int i=0; i<e->size(); i++)
            if ((*e)[i] != 0.)
              (*e)[i] = std::sqrt((*e)[i]);
          STOP_LOG("std::sqrt()", "UniformRefinementEstimator");
        }
    }

  // Restore old solutions and clean up the heap
  for (unsigned int i=0; i != system_list.size(); ++i)
    {
      System &system = *system_list[i];

      system.project_solution_on_reinit() = old_projection_settings[i];
  
      // Restore the coarse solution vectors and delete their copies
      *system.solution = *coarse_solutions[i];
      delete coarse_solutions[i];
      *system.current_local_solution = *coarse_local_solutions[i];
      delete coarse_local_solutions[i];
      delete projected_solutions[i];

      for (System::vectors_iterator vec = system.vectors_begin(); vec !=
           system.vectors_end(); ++vec)
        {
          // The (string) name of this vector
          const std::string& var_name = vec->first;

          system.get_vector(var_name) = *coarse_vectors[i][var_name];

          coarse_vectors[i][var_name]->clear();
          delete coarse_vectors[i][var_name];
        }
    }

  if (!_component_scales)
    delete component_scales;
}

#endif // #ifdef ENABLE_AMR