patch_recovery_error_estimator.c

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

	  // here unless asserts are active.
	  const std::vector<std::vector<Real> >&         phi     = fe->get_phi();
#endif
	  const std::vector<std::vector<RealGradient> > *dphi = NULL;
          if (error_estimator._sobolev_order == 1)
            dphi = &(fe->get_dphi());
#ifdef ENABLE_SECOND_DERIVATIVES
	  const std::vector<std::vector<RealTensor> >  *d2phi = NULL;
          if (error_estimator._sobolev_order == 2)
            d2phi = &(fe->get_d2phi());
#endif
      
	  // global DOF indices
	  std::vector<unsigned int> dof_indices;

	  // Compute the approprite size for the patch projection matrices
	  // and vectors; 
	  unsigned int matsize = element_order + 1;
	  if (dim > 1)
	    {
	      matsize *= (element_order + 2);
	      matsize /= 2;
	    }
	  if (dim > 2)
	    {
	      matsize *= (element_order + 3);
	      matsize /= 3;
	    }
	  	  
	  DenseMatrix<Number> Kp(matsize,matsize);
	  DenseVector<Number>
	    Fx(matsize), Pu_x_h(matsize), // Also xx
	    Fy(matsize), Pu_y_h(matsize), // Also yy
	    Fz(matsize), Pu_z_h(matsize); // Also zz
          DenseVector<Number> Fxy, Pu_xy_h, Fxz, Pu_xz_h, Fyz, Pu_yz_h;
          if (error_estimator._sobolev_order == 2)
            {
              Fxy.resize(matsize); Pu_xy_h.resize(matsize);
              Fxz.resize(matsize); Pu_xz_h.resize(matsize);
              Fyz.resize(matsize); Pu_yz_h.resize(matsize);
            }

	  //------------------------------------------------------
	  // Loop over each element in the patch and compute their
	  // contribution to the patch gradient projection.
	  Patch::const_iterator        patch_it  = patch.begin();
	  const Patch::const_iterator  patch_end = patch.end();

	  for (; patch_it != patch_end; ++patch_it)
	    {
	      // The pth element in the patch
	      const Elem* e_p = *patch_it;

	      // Reinitialize the finite element data for this element
	      fe->reinit (e_p);

	      // Get the global DOF indices for the current variable
	      // in the current element
	      dof_map.dof_indices (e_p, dof_indices, var);
	      libmesh_assert (dof_indices.size() == phi.size());

	      const unsigned int n_dofs = dof_indices.size();
	      const unsigned int n_qp   = qrule->n_points();

	      // Compute the projection components from this cell.
	      // \int_{Omega_e} \psi_i \psi_j = \int_{Omega_e} du_h/dx_k \psi_i
	      for (unsigned int qp=0; qp<n_qp; qp++)
		{
		  // Construct the shape function values for the patch projection
		  std::vector<Real> psi(specpoly(dim, element_order, q_point[qp], matsize));
		  
		  // Patch matrix contribution
		  for (unsigned int i=0; i<Kp.m(); i++)
		    for (unsigned int j=0; j<Kp.n(); j++)
		      Kp(i,j) += JxW[qp]*psi[i]*psi[j];

                  if (error_estimator._sobolev_order == 1)
                    {
		      // Compute the gradient on the current patch element
		      // at the quadrature point
		      Gradient grad_u_h;

		      for (unsigned int i=0; i<n_dofs; i++)
		        grad_u_h.add_scaled ((*dphi)[i][qp],
					     system.current_solution(dof_indices[i]));
		  
		      // Patch RHS contributions
		      for (unsigned int i=0; i<psi.size(); i++)
		        {
		          Fx(i) += JxW[qp]*grad_u_h(0)*psi[i];
		          Fy(i) += JxW[qp]*grad_u_h(1)*psi[i];
		          Fz(i) += JxW[qp]*grad_u_h(2)*psi[i];
		        }
                    }
                  else if (error_estimator._sobolev_order == 2)
                    {
#ifdef ENABLE_SECOND_DERIVATIVES
		      // Compute the hessian on the current patch element
		      // at the quadrature point
                      Tensor hess_u_h;

		      for (unsigned int i=0; i<n_dofs; i++)
		        hess_u_h.add_scaled ((*d2phi)[i][qp],
					     system.current_solution(dof_indices[i]));

		      // Patch RHS contributions
		      for (unsigned int i=0; i<psi.size(); i++)
		        {
		          Fx(i)  += JxW[qp]*hess_u_h(0,0)*psi[i];
		          Fy(i)  += JxW[qp]*hess_u_h(1,1)*psi[i];
		          Fz(i)  += JxW[qp]*hess_u_h(2,2)*psi[i];
		          Fxy(i) += JxW[qp]*hess_u_h(0,1)*psi[i];
		          Fxz(i) += JxW[qp]*hess_u_h(0,2)*psi[i];
		          Fyz(i) += JxW[qp]*hess_u_h(1,2)*psi[i];
		        }
#else
		      std::cerr << "ERROR:  --enable-second-derivatives is required\n"
				<< "        for _sobolev_order == 2!\n";
		      libmesh_error();
#endif
                    }
		} // end quadrature loop
	    } // end patch loop
	  

	  
	  //--------------------------------------------------
	  // Now we have fully assembled the projection system
	  // for this patch.  Project the gradient components.
	  // MAY NEED TO USE PARTIAL PIVOTING!
	  Kp.lu_solve (Fx, Pu_x_h);
	  Kp.lu_solve (Fy, Pu_y_h);
	  Kp.lu_solve (Fz, Pu_z_h);
          if (error_estimator._sobolev_order == 2)
            {
              Kp.lu_solve(Fxy, Pu_xy_h);
              Kp.lu_solve(Fxz, Pu_xz_h);
              Kp.lu_solve(Fyz, Pu_yz_h);
            }
	  
	  //--------------------------------------------------
	  // Finally, estimate the error in the current variable
	  // for the current element by computing ||P grad_u_h - grad_u_h||
          // or ||P hess_u_h - hess_u_h|| in the infinity (max) norm

	  fe->reinit(elem);
	  //reinitialize element
	  
	  dof_map.dof_indices (elem, dof_indices, var);
	  const unsigned int n_dofs = dof_indices.size();
	  
	  // For linear elments, grad is a constant, so we need to compute
	  // grad u_h once on the element.  Also as G_H u_h - gradu_h is linear
	  // on an element, it assumes its maximum at a vertex of the element
	  Real element_error = 0;
	  // we approximate the max norm by sampling over a set of points
	  // in future we may add specialized routines for specific cases
	  // or use some optimization package
	  const Order qorder = element_order;

	  // build a "fake" quadrature rule for the element
	  QGrid samprule (dim, qorder);
	  fe->attach_quadrature_rule (&samprule);
	  fe->reinit(elem);
	      
	  const unsigned int n_sp = samprule.n_points();
	  for (unsigned int sp=0; sp< n_sp; sp++)
	    {
	      std::vector<Number> temperr(6,0.0); // x,y,z or xx,yy,zz,xy,xz,yz
	  
              if (error_estimator._sobolev_order == 1)
                {
	          // Compute the gradient at the current sample point
	          Gradient grad_u_h;
	  
	          for (unsigned int i=0; i<n_dofs; i++)
	            grad_u_h.add_scaled ((*dphi)[i][sp],
			                 system.current_solution(dof_indices[i]));
	          // Compute the phi values at the current sample point
	          std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
	          for (unsigned int i=0; i<matsize; i++)
	            {
	              temperr[0] += psi[i]*Pu_x_h(i);
	              temperr[1] += psi[i]*Pu_y_h(i);
	              temperr[2] += psi[i]*Pu_z_h(i);
	            }
	          temperr[0] -= grad_u_h(0);
	          temperr[1] -= grad_u_h(1);
	          temperr[2] -= grad_u_h(2);
                }
              else if (error_estimator._sobolev_order == 2)
                {
#ifdef ENABLE_SECOND_DERIVATIVES
	          // Compute the Hessian at the current sample point
	          Tensor hess_u_h;
	  
	          for (unsigned int i=0; i<n_dofs; i++)
	            hess_u_h.add_scaled ((*d2phi)[i][sp],
			                 system.current_solution(dof_indices[i]));
	          // Compute the phi values at the current sample point
	          std::vector<Real> psi(specpoly(dim, element_order, q_point[sp], matsize));
	          for (unsigned int i=0; i<matsize; i++)
	            {
	              temperr[0] += psi[i]*Pu_x_h(i);
	              temperr[1] += psi[i]*Pu_y_h(i);
	              temperr[2] += psi[i]*Pu_z_h(i);
	              temperr[3] += psi[i]*Pu_xy_h(i);
	              temperr[4] += psi[i]*Pu_xz_h(i);
	              temperr[5] += psi[i]*Pu_yz_h(i);
	            }
	          temperr[0] -= hess_u_h(0,0);
	          temperr[1] -= hess_u_h(1,1);
	          temperr[2] -= hess_u_h(2,2);
	          temperr[3] -= hess_u_h(0,1);
	          temperr[4] -= hess_u_h(0,2);
	          temperr[5] -= hess_u_h(1,2);
#else
		      std::cerr << "ERROR:  --enable-second-derivatives is required\n"
				<< "        for _sobolev_order == 2!\n";
		      libmesh_error();
#endif
                }
              for (unsigned int i=0; i != 6; ++i)
                element_error = std::max(element_error, std::abs(temperr[i]));

	    } // end sample_point_loop

	  // The patch error estimator works element-by-element --
	  // there is no need to get a mutex on error_per_cell!
	  const int e_id=elem->id();
	  error_per_cell[e_id] += element_error;	  
	} // end variable loop  
    } // end element loop
}