jump_error_estimator.c

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

                      parent->neighbor(n_p)->
		        active_family_tree_by_neighbor(active_neighbors,
                                                       parent);
                      // Compute the flux to each active neighbor
                      for (unsigned int a=0; 
                           a != active_neighbors.size(); ++a)
                        {
                          const Elem *f = active_neighbors[a];
                      // FIXME - what about when f->level <
                      // parent->level()??
                          if (f->level() >= parent->level())
                            {
                              fine_elem = f;
                              coarse_elem = parent;
                              Ucoarse = Uparent;

		              dof_map.dof_indices (fine_elem, dof_indices_fine, var);
		              const unsigned int n_dofs_fine = dof_indices_fine.size();
                              Ufine.resize(n_dofs_fine);
			
			      for (unsigned int i=0; i<n_dofs_fine; i++)
			        Ufine(i) = system.current_solution(dof_indices_fine[i]);
                              this->reinit_sides();
                              this->internal_side_integration();

                              error_per_cell[fine_elem->id()] += fine_error;
                              error_per_cell[coarse_elem->id()] += coarse_error;

                              // Keep track of the number of internal flux
                              // sides found on each element
                              n_flux_faces[fine_elem->id()]++;
                              n_flux_faces[coarse_elem->id()] += this->coarse_n_flux_faces_increment();
                            }
                        }
		    }
		  else if (integrate_boundary_sides)
		    {
                      fine_elem = parent;
                      Ufine = Uparent;

                      // Reinitialize shape functions on the fine element side
                      fe_fine->reinit (fine_elem, fine_side);

                      if (this->boundary_side_integration())
                        {
                          error_per_cell[fine_elem->id()] += fine_error;
                          n_flux_faces[fine_elem->id()]++;
                        }
                    } 
		}
	    }
#endif // #ifdef ENABLE_AMR

          // If we do any more flux integration, e will be the fine element
          fine_elem = e;

	  // Loop over the neighbors of element e
	  for (unsigned int n_e=0; n_e<e->n_neighbors(); n_e++)
	    {
              fine_side = n_e;

	      if (e->neighbor(n_e) != NULL) // e is not on the boundary
		{
		  const Elem* f           = e->neighbor(n_e);
		  const unsigned int f_id = f->id();

		  // Compute flux jumps if we are in case 1 or case 2.
		  if ((f->active() && (f->level() == e->level()) && (e_id < f_id))
		      || (f->level() < e->level()))
		    {		    
                      // f is now the coarse element
                      coarse_elem = f;

		      // Get the DOF indices for the two elements
		      dof_map.dof_indices (fine_elem, dof_indices_fine, var);
		      dof_map.dof_indices (coarse_elem, dof_indices_coarse, var);

		      // The number of DOFS on each element
		      const unsigned int n_dofs_fine = dof_indices_fine.size();
		      const unsigned int n_dofs_coarse = dof_indices_coarse.size();
                      Ufine.resize(n_dofs_fine);
                      Ucoarse.resize(n_dofs_coarse);

		      // The local solutions on each element
		      for (unsigned int i=0; i<n_dofs_fine; i++)
			Ufine(i) = system.current_solution(dof_indices_fine[i]);
		      for (unsigned int i=0; i<n_dofs_coarse; i++)
			Ucoarse(i) = system.current_solution(dof_indices_coarse[i]);
			
                      this->reinit_sides();
                      this->internal_side_integration();

                      error_per_cell[fine_elem->id()] += fine_error;
                      error_per_cell[coarse_elem->id()] += coarse_error;

                      // Keep track of the number of internal flux
                      // sides found on each element
                      n_flux_faces[fine_elem->id()]++;
                      n_flux_faces[coarse_elem->id()] += this->coarse_n_flux_faces_increment();
		    } // end if (case1 || case2)
		} // if (e->neigbor(n_e) != NULL)

	      // Otherwise, e is on the boundary.  If it happens to
	      // be on a Dirichlet boundary, we need not do anything.
	      // On the other hand, if e is on a Neumann (flux) boundary
	      // with grad(u).n = g, we need to compute the additional residual
	      // (h * \int |g - grad(u_h).n|^2 dS)^(1/2).
	      // We can only do this with some knowledge of the boundary
	      // conditions, i.e. the user must have attached an appropriate
	      // BC function.
	      else
		{
		  if (integrate_boundary_sides)
		    {
                      // Reinitialize shape functions on the fine element side
                      fe_fine->reinit (fine_elem, fine_side);

		      // Get the DOF indices 
		      dof_map.dof_indices (fine_elem, dof_indices_fine, var);

		      // The number of DOFS on each element
		      const unsigned int n_dofs_fine = dof_indices_fine.size();
                      Ufine.resize(n_dofs_fine);

                      for (unsigned int i=0; i<n_dofs_fine; i++)
                        Ufine(i) = system.current_solution(dof_indices_fine[i]);

                      if (this->boundary_side_integration())
                        {
                          error_per_cell[fine_elem->id()] += fine_error;
                          n_flux_faces[fine_elem->id()]++;
                        }
                    } // end if _bc_function != NULL
		} // end if (e->neighbor(n_e) == NULL)
	    } // end loop over neighbors
	} // End loop over active local elements
    } // End loop over variables


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

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

  // Compute the square-root of each component.
  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]);


  if (this->scale_by_n_flux_faces)
    {
      // Sum the vector of flux face counts
      this->reduce_error(n_flux_faces);

      // Sanity check: Make sure the number of flux faces is
      // always an integer value
#ifdef DEBUG
      for (unsigned int i=0; i<n_flux_faces.size(); ++i)
	libmesh_assert (n_flux_faces[i] == static_cast<float>(static_cast<unsigned int>(n_flux_faces[i])) );
#endif
  
      // Scale the error by the number of flux faces for each element
      for (unsigned int i=0; i<n_flux_faces.size(); ++i)
	{
	  if (n_flux_faces[i] == 0.0) // inactive or non-local element
	    continue;
      
	  //std::cout << "Element " << i << " has " << n_flux_faces[i] << " flux faces." << std::endl;
	  error_per_cell[i] /= static_cast<Real>(n_flux_faces[i]); 
	}
    }
  
  STOP_LOG("estimate_error()", "JumpErrorEstimator");
}



void
JumpErrorEstimator::reinit_sides ()
{
  // The master quadrature point locations on the coarse element
  std::vector<Point> qp_coarse;

  // Reinitialize shape functions on the fine element side
  fe_fine->reinit (fine_elem, fine_side);

  // Get the physical locations of the fine element quadrature points
  std::vector<Point> qface_point = fe_fine->get_xyz();

  // Find their locations on the coarse element
  FEInterface::inverse_map (coarse_elem->dim(), fe_coarse->get_fe_type(),
                            coarse_elem, qface_point, qp_coarse);

  // Calculate the coarse element shape functions at those locations
  fe_coarse->reinit (coarse_elem, &qp_coarse);
}



float JumpErrorEstimator::coarse_n_flux_faces_increment ()
{
  // Keep track of the number of internal flux sides found on each
  // element
  unsigned int dim = coarse_elem->dim();

  const unsigned int divisor =
    1 << (dim-1)*(fine_elem->level() - coarse_elem->level());

  // With a difference of n levels between fine and coarse elements,
  // we compute a fractional flux face for the coarse element by adding:
  // 1/2^n in 2D
  // 1/4^n in 3D
  // each time.  This code will get hit 2^n times in 2D and 4^n
  // times in 3D so that the final flux face count for the coarse
  // element will be an integer value.

  return 1.0 / static_cast<Real>(divisor);
}