C++ Quad2D::get_points方法代码示例

// Actual evaluation of surface linear form (calculates integral)
scalar FeProblem::eval_form(WeakForm::VectorFormSurf *vfs, Tuple<Solution *> u_ext, PrecalcShapeset *fv, RefMap *rv, EdgePos* ep)
{
  // eval the form
  Quad2D* quad = fv->get_quad_2d();
  int eo = quad->get_edge_points(ep->edge);
  double3* pt = quad->get_points(eo);
  int np = quad->get_num_points(eo);

  // init geometry and jacobian*weights
  if (cache_e[eo] == NULL)
  {
    cache_e[eo] = init_geom_surf(rv, ep, eo);
    double3* tan = rv->get_tangent(ep->edge);
    cache_jwt[eo] = new double[np];
    for(int i = 0; i < np; i++)
      cache_jwt[eo][i] = pt[i][2] * tan[i][2];
  }
  Geom<double>* e = cache_e[eo];
  double* jwt = cache_jwt[eo];

  // function values and values of external functions
  AUTOLA_OR(Func<scalar>*, prev, wf->neq);
  for (int i = 0; i < wf->neq; i++) prev[i]  = init_fn(u_ext[i], rv, eo);
  Func<double>* v = get_fn(fv, rv, eo);
  ExtData<scalar>* ext = init_ext_fns(vfs->ext, rv, eo);

  scalar res = vfs->fn(np, jwt, prev, v, e, ext);

  for (int i = 0; i < wf->neq; i++) {  prev[i]->free_fn(); delete prev[i]; }
  ext->free(); delete ext;
  return 0.5 * res;
}

// Custom function to calculate drag coefficient.
double integrate_over_wall(MeshFunction* meshfn, int marker)
{
  Quad2D* quad = &g_quad_2d_std;
  meshfn->set_quad_2d(quad);

  double integral = 0.0;
  Element* e;
  Mesh* mesh = meshfn->get_mesh();

  for_all_active_elements(e, mesh)
  {
    for(int edge = 0; edge < e->nvert; edge++)
    {
      if ((e->en[edge]->bnd) && (e->en[edge]->marker == marker))
      {
        update_limit_table(e->get_mode());
        RefMap* ru = meshfn->get_refmap();

        meshfn->set_active_element(e);
        int eo = quad->get_edge_points(edge);
        meshfn->set_quad_order(eo, H2D_FN_VAL);
        scalar *uval = meshfn->get_fn_values();
        double3* pt = quad->get_points(eo);
        double3* tan = ru->get_tangent(edge);
        for (int i = 0; i < quad->get_num_points(eo); i++)
          integral += pt[i][2] * uval[i] * tan[i][2];
      }
    }
  }
  return integral * 0.5;
}

//------------------------------------------------------------------------------
// Compute marked boundary length 
//
double CalculateBoundaryLength(Mesh* mesh, int bdryMarker)
{
  // Variables declaration.
  Element* e;
  double length = 0;
  RefMap rm;
  rm.set_quad_2d(&g_quad_2d_std);
  Quad2D * quad = rm.get_quad_2d();
  int points_location;
  double3* points;
  int np;
  double3* tangents;

  // Loop through all boundary faces of all active elements.
  for_all_active_elements(e, mesh) {
    for(int edge = 0; edge < e->nvert; ++edge) {
      if ((e->en[edge]->bnd) && (e->en[edge]->marker == bdryMarker)) {
        rm.set_active_element(e);
        points_location = quad->get_edge_points(edge);
        points = quad->get_points(points_location);
        np = quad->get_num_points(points_location);
        tangents = rm.get_tangent(edge, points_location);
        for(int i = 0; i < np; i++) {
          // Weights sum up to two on every edge, therefore the division by two must be present.
          length +=  0.5 * points[i][2] * tangents[i][2];
        }
      }
    }
  }
  return length;
} // end of CalculateBoundaryLength()

scalar HcurlOrthoHP::eval_error(biform_val_t bi_fn, biform_ord_t bi_ord,
                             MeshFunction *sln1, MeshFunction *sln2, MeshFunction *rsln1, MeshFunction *rsln2,
                             RefMap *rv1,        RefMap *rv2,        RefMap *rrv1,        RefMap *rrv2)
{
  // determine the integration order
  int inc = (rsln1->get_num_components() == 2) ? 1 : 0;
  Func<Ord>* ou = init_fn_ord(rsln1->get_fn_order() + inc);
  Func<Ord>* ov = init_fn_ord(rsln2->get_fn_order() + inc);

  double fake_wt = 1.0;
  Geom<Ord>* fake_e = init_geom_ord();
  Ord o = bi_ord(1, &fake_wt, ou, ov, fake_e, NULL);
  int order = rrv1->get_inv_ref_order();
  order += o.get_order();
  limit_order(order);

  ou->free_ord(); delete ou;
  ov->free_ord(); delete ov;
  delete fake_e;

  // eval the form
  Quad2D* quad = sln1->get_quad_2d();
  double3* pt = quad->get_points(order);
  int np = quad->get_num_points(order);

  // init geometry and jacobian*weights
  Geom<double>* e = init_geom_vol(rrv1, order);
  double* jac = rrv1->get_jacobian(order);
  double* jwt = new double[np];
  for(int i = 0; i < np; i++)
    jwt[i] = pt[i][2] * jac[i];

  // function values and values of external functions
  Func<scalar>* err1 = init_fn(sln1, rv1, order);
  Func<scalar>* err2 = init_fn(sln2, rv2, order);
  Func<scalar>* v1 = init_fn(rsln1, rrv1, order);
  Func<scalar>* v2 = init_fn(rsln2, rrv2, order);

  for (int i = 0; i < np; i++)
  {
    err1->val0[i] = err1->val0[i] - v1->val0[i];
    err1->val1[i] = err1->val1[i] - v1->val1[i];
    err1->curl[i] = err1->curl[i] - v1->curl[i];
    err2->val0[i] = err2->val0[i] - v2->val0[i];
    err2->val1[i] = err2->val1[i] - v2->val1[i];
    err2->curl[i] = err2->curl[i] - v2->curl[i];
  }

  scalar res = bi_fn(np, jwt, err1, err2, e, NULL);

  e->free(); delete e;
  delete [] jwt;
  err1->free_fn(); delete err1;
  err2->free_fn(); delete err2;
  v1->free_fn(); delete v1;
  v2->free_fn(); delete v2;

  return res;
}

double KellyTypeAdapt::eval_boundary_estimator(KellyTypeAdapt::ErrorEstimatorForm* err_est_form, RefMap *rm, SurfPos* surf_pos)
{
  // determine the integration order
  int inc = (this->sln[err_est_form->i]->get_num_components() == 2) ? 1 : 0;
  Func<Ord>** oi = new Func<Ord>* [num];
  for (int i = 0; i < num; i++)
    oi[i] = init_fn_ord(this->sln[i]->get_edge_fn_order(surf_pos->surf_num) + inc);

  // Order of additional external functions.
  ExtData<Ord>* fake_ext = dp.init_ext_fns_ord(err_est_form->ext, surf_pos->surf_num);

  double fake_wt = 1.0;
  Geom<Ord>* fake_e = init_geom_ord();
  Ord o = err_est_form->ord(1, &fake_wt, oi, oi[err_est_form->i], fake_e, fake_ext);
  int order = rm->get_inv_ref_order();
  order += o.get_order();

  limit_order(order);

  // Clean up.
  for (int i = 0; i < this->num; i++)
    if (oi[i] != NULL) { oi[i]->free_ord(); delete oi[i]; }
  delete [] oi;
  delete fake_e;
  delete fake_ext;

  // eval the form
  Quad2D* quad = this->sln[err_est_form->i]->get_quad_2d();
  int eo = quad->get_edge_points(surf_pos->surf_num, order);
  double3* pt = quad->get_points(eo);
  int np = quad->get_num_points(eo);

  // init geometry and jacobian*weights
  Geom<double>* e = init_geom_surf(rm, surf_pos, eo);
  double3* tan = rm->get_tangent(surf_pos->surf_num, eo);
  double* jwt = new double[np];
  for(int i = 0; i < np; i++)
    jwt[i] = pt[i][2] * tan[i][2];

  // function values
  Func<scalar>** ui = new Func<scalar>* [num];
  for (int i = 0; i < num; i++)
    ui[i] = init_fn(this->sln[i], eo);
  ExtData<scalar>* ext = dp.init_ext_fns(err_est_form->ext, rm, eo);

  scalar res = boundary_scaling_const *
                err_est_form->value(np, jwt, ui, ui[err_est_form->i], e, ext);

  for (int i = 0; i < this->num; i++)
    if (ui[i] != NULL) { ui[i]->free_fn(); delete ui[i]; }
  delete [] ui;
  if (ext != NULL) { ext->free(); delete ext; }
  e->free(); delete e;
  delete [] jwt;

  return std::abs(0.5*res);   // Edges are parameterized from 0 to 1 while integration weights
                              // are defined in (-1, 1). Thus multiplying with 0.5 to correct
                              // the weights.
}

double KellyTypeAdapt::eval_volumetric_estimator(KellyTypeAdapt::ErrorEstimatorForm* err_est_form, RefMap *rm)
{
  // determine the integration order
  int inc = (this->sln[err_est_form->i]->get_num_components() == 2) ? 1 : 0;

  Func<Ord>** oi = new Func<Ord>* [num];
  for (int i = 0; i < num; i++)
    oi[i] = init_fn_ord(this->sln[i]->get_fn_order() + inc);

  // Order of additional external functions.
  ExtData<Ord>* fake_ext = dp.init_ext_fns_ord(err_est_form->ext);

  double fake_wt = 1.0;
  Geom<Ord>* fake_e = init_geom_ord();
  Ord o = err_est_form->ord(1, &fake_wt, oi, oi[err_est_form->i], fake_e, fake_ext);
  int order = rm->get_inv_ref_order();
  order += o.get_order();

  limit_order(order);

  // Clean up.
  for (int i = 0; i < this->num; i++)
    if (oi[i] != NULL) { oi[i]->free_ord(); delete oi[i]; }
  delete [] oi;
  delete fake_e;
  delete fake_ext;

  // eval the form
  Quad2D* quad = this->sln[err_est_form->i]->get_quad_2d();
  double3* pt = quad->get_points(order);
  int np = quad->get_num_points(order);

  // init geometry and jacobian*weights
  Geom<double>* e = init_geom_vol(rm, order);
  double* jac = rm->get_jacobian(order);
  double* jwt = new double[np];
  for(int i = 0; i < np; i++)
    jwt[i] = pt[i][2] * jac[i];

  // function values
  Func<scalar>** ui = new Func<scalar>* [num];
  
  for (int i = 0; i < num; i++)
    ui[i] = init_fn(this->sln[i], order);
  
  ExtData<scalar>* ext = dp.init_ext_fns(err_est_form->ext, rm, order);

  scalar res = volumetric_scaling_const *
                err_est_form->value(np, jwt, ui, ui[err_est_form->i], e, ext);

  for (int i = 0; i < this->num; i++)
    if (ui[i] != NULL) { ui[i]->free_fn(); delete ui[i]; }
  delete [] ui;
  if (ext != NULL) { ext->free(); delete ext; }
  e->free(); delete e;
  delete [] jwt;

  return std::abs(res);
}

int RefMap::calc_inv_ref_order()
{
  Quad2D* quad = get_quad_2d();
  int i, o, mo = quad->get_max_order();

  // check first the positivity of the jacobian
  double3* pt = quad->get_points(mo);
  double2x2* m = get_inv_ref_map(mo);
  double* jac = get_jacobian(mo);
  for (i = 0; i < quad->get_num_points(mo); i++)
    if (jac[i] <= 0.0)
      error("Element #%d is concave or badly oriented.", element->id);

  // next, estimate the "exact" value of the typical integral int_grad_u_grad_v
  // (with grad_u == grad_v == (1,1)) using the maximum integration rule
  double exact1 = 0.0;
  double exact2 = 0.0;
  for (i = 0; i < quad->get_num_points(mo); i++, m++)
  {
    exact1 += pt[i][2] * jac[i] * (sqr((*m)[0][0] + (*m)[0][1]) + sqr((*m)[1][0] + (*m)[1][1]));
    exact2 += pt[i][2] / jac[i];
  }
  // find sufficient quadrature degree
  for (o = 0; o < mo; o++)
  {
    pt = quad->get_points(o);
    m = get_inv_ref_map(o);
    jac = get_jacobian(o);
    double result1 = 0.0;
    double result2 = 0.0;
    for (i = 0; i < quad->get_num_points(o); i++, m++)
    {
      result1 += pt[i][2] * jac[i] * (sqr((*m)[0][0] + (*m)[0][1]) + sqr((*m)[1][0] + (*m)[1][1]));
      result2 += pt[i][2] / jac[i] ;
    }
    if ((fabs((exact1 - result1) / exact1) < 1e-8) &&
        (fabs((exact2 - result2) / exact2) < 1e-8)) break;
  }
  if (o >= 10) {
    warn("Element #%d is too distorted (iro ~ %d).", element->id, o);
  }
  return o;
}

double KellyTypeAdapt::eval_solution_norm(Adapt::MatrixFormVolError* form, RefMap *rm, MeshFunction* sln)
{
  // determine the integration order
  int inc = (sln->get_num_components() == 2) ? 1 : 0;
  Func<Ord>* ou = init_fn_ord(sln->get_fn_order() + inc);

  double fake_wt = 1.0;
  Geom<Ord>* fake_e = init_geom_ord();
  Ord o = form->ord(1, &fake_wt, NULL, ou, ou, fake_e, NULL);
  int order = rm->get_inv_ref_order();
  order += o.get_order();

  Solution *sol = static_cast<Solution *>(sln);
  if(sol && sol->get_type() == HERMES_EXACT) {
    limit_order_nowarn(order);
  }
  else {
    limit_order(order);
  }

  ou->free_ord(); delete ou;
  delete fake_e;

  // eval the form
  Quad2D* quad = sln->get_quad_2d();
  double3* pt = quad->get_points(order);
  int np = quad->get_num_points(order);

  // init geometry and jacobian*weights
  Geom<double>* e = init_geom_vol(rm, order);
  double* jac = rm->get_jacobian(order);
  double* jwt = new double[np];
  for(int i = 0; i < np; i++)
    jwt[i] = pt[i][2] * jac[i];

  // function values
  Func<scalar>* u = init_fn(sln, order);
  scalar res = form->value(np, jwt, NULL, u, u, e, NULL);

  e->free(); delete e;
  delete [] jwt;
  u->free_fn(); delete u;

  return std::abs(res);
}

// Preparation of mesh-functions
Func<scalar>* init_fn(MeshFunction *fu, RefMap *rm, const int order)
{
  // sanity checks
  if (fu == NULL) error("NULL MeshFunction in Func<scalar>*::init_fn().");
  if (fu->get_mesh() == NULL) error("Uninitialized MeshFunction used.");

  int nc = fu->get_num_components();
  Quad2D* quad = fu->get_quad_2d();
  fu->set_quad_order(order);
  double3* pt = quad->get_points(order);
  int np = quad->get_num_points(order);
  Func<scalar>* u = new Func<scalar>(np, nc);

  if (u->nc == 1)
  {
    u->val = new scalar [np];
    u->dx  = new scalar [np];
    u->dy  = new scalar [np];
    memcpy(u->val, fu->get_fn_values(), np * sizeof(scalar));
    memcpy(u->dx, fu->get_dx_values(), np * sizeof(scalar));
    memcpy(u->dy, fu->get_dy_values(), np * sizeof(scalar));
  }
  else if (u->nc == 2)
  {
    u->val0 = new scalar [np];
    u->val1 = new scalar [np];
    u->curl = new scalar [np];
    u->div = new scalar [np];

    memcpy(u->val0, fu->get_fn_values(0), np * sizeof(scalar));
    memcpy(u->val1, fu->get_fn_values(1), np * sizeof(scalar));

    // This works.
    scalar *dx1 = fu->get_dx_values(1);
    scalar *dy0 = fu->get_dy_values(0);
    for (int i = 0; i < np; i++) u->curl[i] = dx1[i] - dy0[i];
    // This was never tested but it should work.
    scalar *dx0 = fu->get_dx_values(0);
    scalar *dy1 = fu->get_dy_values(1);
    for (int i = 0; i < np; i++) u->div[i] = dx0[i] + dy1[i];
  }
  return u;
}

void PrecalcShapeset::precalculate(int order, int mask)
{
  int i, j, k;

  // initialization
  Quad2D* quad = get_quad_2d();
  quad->set_mode(mode);
  H2D_CHECK_ORDER(quad, order);
  int np = quad->get_num_points(order);
  double3* pt = quad->get_points(order);

  int oldmask = (cur_node != NULL) ? cur_node->mask : 0;
  int newmask = mask | oldmask;
  Node* node = new_node(newmask, np);

  // precalculate all required tables
  for (j = 0; j < num_components; j++)
  {
    for (k = 0; k < 6; k++)
    {
      if (newmask & idx2mask[k][j]) {
        if (oldmask & idx2mask[k][j])
          memcpy(node->values[j][k], cur_node->values[j][k], np * sizeof(double));
        else
          for (i = 0; i < np; i++)
            node->values[j][k][i] = shapeset->get_value(k, index, ctm->m[0] * pt[i][0] + ctm->t[0],
                                                                  ctm->m[1] * pt[i][1] + ctm->t[1], j);
      }
    }
  }
  if(nodes->present(order)) {
    assert(nodes->get(order) == cur_node);
    ::free(nodes->get(order));
  }
  nodes->add(node, order);
  cur_node = node;
}

void HcurlOrthoHP::calc_projection_errors(Element* e, int order, Solution* rsln,
                                   double herr[8][11], double perr[11])
{
  int i, j, k, son, s;
  int m = e->get_mode();
  double error;
  scalar prod;

  if (!obase_ready) calc_ortho_base();

  // select quadrature, obtain integration points and weights
  Quad2D* quad = &g_quad_2d_std;
  quad->set_mode(m);
  rsln->set_quad_2d(quad);
  double3* pt = quad->get_points(20);
  int np = quad->get_num_points(20);

  // everything is done on the reference domain
  // no reference mapping, no transformations
  rsln->enable_transform(false);

  // obtain reference solution values on all four refined sons
  scalar* rval0[4], *rval1[4], *rd1dx[4], *rd0dy[4];
  Element* base = rsln->get_mesh()->get_element(e->id);
  assert(!base->active);
  for (son = 0; son < 4; son++)
  {
    Element* e = base->sons[son];
    assert(e != NULL);
    rsln->set_active_element(e);
    rsln->set_quad_order(20);
    rval0[son] = rsln->get_fn_values(0);
    rval1[son] = rsln->get_fn_values(1);
    rd1dx[son] = rsln->get_dx_values(1);
    rd0dy[son] = rsln->get_dy_values(0);
  }

  // h-candidates: calculate products of the reference solution with orthonormal basis
  // functions on son elements, obtaining (partial) projections and their errors
  scalar proj_0[4][121];
  scalar proj_1[4][121];
  scalar proj_c[4][121];
  for (son = 0; son < 4; son++)
  {
    memset(proj_0[0], 0, sizeof(proj_0[0]));
    memset(proj_1[0], 0, sizeof(proj_1[0]));
    memset(proj_c[0], 0, sizeof(proj_c[0]));

    for (i = 0; i <= order; i++)  // h-candidates: max order equals to original element order
    {
      // update the projection to the current order
      for (j = basecnt[m][i]; j < basecnt[m][i+1]; j++)
      {
        for (k = 0, prod = 0.0; k < np; k++)
        {
          scalar rcurl = (rd1dx[son][k] - rd0dy[son][k]);
          scalar r0 = rval0[son][k];
          scalar r1 = rval1[son][k];
          prod += pt[k][2] * (( r0    * obase_0[m][8][j][k] ) +
                              ( r1    * obase_1[m][8][j][k] ) +
                              ( rcurl * obase_c[m][8][j][k] ) );
        }
        for (k = 0; k < np; k++)
        {
          proj_0[0][k] += obase_0[m][8][j][k] * prod;
          proj_1[0][k] += obase_1[m][8][j][k] * prod;
          proj_c[0][k] += obase_c[m][8][j][k] * prod;
        }
      }

      // calculate the H(curl) error of the projection
      for (k = 0, error = 0.0; k < np; k++)
      {
        scalar rcurl = (rd1dx[son][k] - rd0dy[son][k]);
        scalar r0 = rval0[son][k];
        scalar r1 = rval1[son][k];
        error += pt[k][2] *  ( sqr(r0    - proj_0[0][k]) +
                               sqr(r1    - proj_1[0][k]) +
                               sqr(rcurl - proj_c[0][k]) );
      }
      herr[son][i] = error;
    }
  }

  // aniso-candidates: calculate projections and their errors (only quadrilaterals)
  if (m) {
    const double mx[4] = { 2.0, 2.0, 1.0, 1.0};
    const double my[4] = { 1.0, 1.0, 2.0, 2.0};
    const int sons[4][2] = {{0,1},{3,2},{0,3},{1,2}};
    const int tr[4][2] = {{6,7},{6,7},{4,5},{4,5}};

    for (son = 0; son < 4; son++) // 2 sons for vertical split, 2 sons for horizontal split
    {
      memset(proj_0, 0, sizeof(proj_0));
      memset(proj_1, 0, sizeof(proj_1));
      memset(proj_c, 0, sizeof(proj_c));

      for (i = 0; i <= order+1; i++)  // h-candidates: max order equals to original element order+1
      {
        // update the projection to the current order
//.........这里部分代码省略.........

double KellyTypeAdapt::eval_interface_estimator(KellyTypeAdapt::ErrorEstimatorForm* err_est_form,
                                                RefMap *rm, SurfPos* surf_pos,
                                                LightArray<NeighborSearch*>& neighbor_searches, int neighbor_index)
{
  NeighborSearch* nbs = neighbor_searches.get(neighbor_index);
  Hermes::vector<MeshFunction*> slns;
  for (int i = 0; i < num; i++)
    slns.push_back(this->sln[i]);
  
  // Determine integration order.
  ExtData<Ord>* fake_ui = dp.init_ext_fns_ord(slns, neighbor_searches);
  
  // Order of additional external functions.
  // ExtData<Ord>* fake_ext = dp.init_ext_fns_ord(err_est_form->ext, nbs);

  // Order of geometric attributes (eg. for multiplication of a solution with coordinates, normals, etc.).
  Geom<Ord>* fake_e = new InterfaceGeom<Ord>(init_geom_ord(), nbs->neighb_el->marker, nbs->neighb_el->id, nbs->neighb_el->get_diameter());
  double fake_wt = 1.0;
  Ord o = err_est_form->ord(1, &fake_wt, fake_ui->fn, fake_ui->fn[err_est_form->i], fake_e, NULL);

  int order = rm->get_inv_ref_order();
  order += o.get_order();

  limit_order(order);

  // Clean up.
  if (fake_ui != NULL)
  {
    for (int i = 0; i < num; i++)
      delete fake_ui->fn[i];
    fake_ui->free_ord();
    delete fake_ui;
  }
  
  delete fake_e;
  
  //delete fake_ext;
  
  Quad2D* quad = this->sln[err_est_form->i]->get_quad_2d();
  int eo = quad->get_edge_points(surf_pos->surf_num, order);
  int np = quad->get_num_points(eo);
  double3* pt = quad->get_points(eo);
  
  // Init geometry and jacobian*weights (do not use the NeighborSearch caching mechanism).
  double3* tan = rm->get_tangent(surf_pos->surf_num, eo);
  double* jwt = new double[np];
  for(int i = 0; i < np; i++)
    jwt[i] = pt[i][2] * tan[i][2];
  
  Geom<double>* e = new InterfaceGeom<double>(init_geom_surf(rm, surf_pos, eo), 
                                              nbs->neighb_el->marker, 
                                              nbs->neighb_el->id, 
                                              nbs->neighb_el->get_diameter());
    
  // function values
  ExtData<scalar>* ui = dp.init_ext_fns(slns, neighbor_searches, order);
  //ExtData<scalar>* ext = dp.init_ext_fns(err_est_form->ext, nbs);

  scalar res = interface_scaling_const *
                err_est_form->value(np, jwt, ui->fn, ui->fn[err_est_form->i], e, NULL);

  if (ui != NULL) { ui->free(); delete ui; }
  //if (ext != NULL) { ext->free(); delete ext; }
  e->free(); delete e;
  delete [] jwt;

  return std::abs(0.5*res);   // Edges are parameterized from 0 to 1 while integration weights
                              // are defined in (-1, 1). Thus multiplying with 0.5 to correct
                              // the weights.
}

void L2OrthoHP::calc_projection_errors(Element* e, int order, Solution* rsln,
                                       double herr[8][11], double perr[11])
{
  int i, j, s, k, r, son;
  int m = e->get_mode();
  double error;
  Scalar prod;

  if (!obase_ready) calc_ortho_base();

  // select quadrature, obtain integration points and weights
  Quad2D* quad = &g_quad_2d_std;
  quad->set_mode(m);
  rsln->set_quad_2d(quad);
  double3* pt = quad->get_points(20);
  int np = quad->get_num_points(20);

  // everything is done on the reference domain
  // -- no reference mapping, no transformations
  rsln->enable_transform(false);

  // obtain reference solution values on all four refined sons
  Scalar* rval[4][3];
  Element* base = rsln->get_mesh()->get_element(e->id);
  assert(!base->active);
  for (son = 0; son < 4; son++)
  {
    Element* e = base->sons[son];
    assert(e != NULL);
    rsln->set_active_element(e);
    rsln->set_quad_order(20);
    rval[son][0] = rsln->get_fn_values();
    rval[son][1] = rsln->get_dx_values();
    rval[son][2] = rsln->get_dy_values();
  }

  // h-cadidates: calculate products of the reference solution with orthonormal basis
  // functions on son elements, obtaining (partial) projections and their errors
  Scalar3 proj[4][121];
  for (son = 0; son < 4; son++)
  {
    memset(proj[0], 0, sizeof(proj[0]));
    for (i = 1; i <= order; i++)
    {
      // update the projection to the current order
      for (j = basecnt[m][i-1]; j < basecnt[m][i]; j++)
      {
        for (k = 0, prod = 0.0; k < np; k++)
          prod += pt[k][2] * (rval[son][0][k] * obase[m][8][j][k][0]);

        for (k = 0; k < np; k++)
          for (r = 0; r < 3; r++)
            proj[0][k][r] += obase[m][8][j][k][r] * prod;
      }

      // calculate the error of the projection
      for (k = 0, error = 0.0; k < np; k++)
        error += pt[k][2] * (sqr(rval[son][0][k] - proj[0][k][0]));
      herr[son][i] = error;
    }
  }

  // aniso-candidates: calculate projections and their errors (only quadrilaterals)
  if (m)
  {
    const double mx[4] = { 2.0, 2.0, 1.0, 1.0};
    const double my[4] = { 1.0, 1.0, 2.0, 2.0};
    const int sons[4][2] = { {0,1}, {3,2}, {0,3}, {1,2} };
    const int tr[4][2]   = { {6,7}, {6,7}, {4,5}, {4,5} };

    for (son = 0; son < 4; son++) // 2 sons for vertical split, 2 sons for horizontal split
    {
      memset(proj, 0, sizeof(proj));
      for (i = 1; i <= order+1; i++)  // h-candidates: max order equals to original element order+1
      {
        // update the projection to the current order
        for (j = basecnt[m][i-1]; j < basecnt[m][i]; j++)
        {
          for (s = 0, prod = 0.0; s < 2; s++) // each son has 2 subsons (regular square sons)
            for (k = 0; k < np; k++)
              prod += pt[k][2] * rval[sons[son][s]][0][k] * obase[m][tr[son][s]][j][k][0];
          prod *= 0.5;

          for (s = 0; s < 2; s++)
            for (k = 0; k < np; k++)
              for (r = 0; r < 1; r++)
                proj[s][k][r] += prod * obase[m][tr[son][s]][j][k][r];
        }

        // calculate the error of the projection
        for (s = 0, error = 0.0; s < 2; s++)
          for (k = 0; k < np; k++)
            error += pt[k][2] * sqr(rval[sons[son][s]][0][k] - proj[s][k][0]);
        herr[4 + son][i] = error * 0.5;
      }
    }
  }

  // p-candidates: calculate projections and their errors
  memset(proj, 0, sizeof(proj));
//.........这里部分代码省略.........

  void ProjBasedSelector::calc_projection_errors(Element* e, const CandsInfo& info_h, const CandsInfo& info_p, const CandsInfo& info_aniso, Solution* rsln, CandElemProjError herr[4], CandElemProjError perr, CandElemProjError anisoerr[4]) {
    assert_msg(info_h.is_empty() || (H2D_GET_H_ORDER(info_h.max_quad_order) <= H2DRS_MAX_ORDER && H2D_GET_V_ORDER(info_h.max_quad_order) <= H2DRS_MAX_ORDER), "Maximum allowed order of a son of H-candidate is %d but order (H:%d,V:%d) requested.", H2DRS_MAX_ORDER, H2D_GET_H_ORDER(info_h.max_quad_order), H2D_GET_V_ORDER(info_h.max_quad_order));
    assert_msg(info_p.is_empty() || (H2D_GET_H_ORDER(info_p.max_quad_order) <= H2DRS_MAX_ORDER && H2D_GET_V_ORDER(info_p.max_quad_order) <= H2DRS_MAX_ORDER), "Maximum allowed order of a son of P-candidate is %d but order (H:%d,V:%d) requested.", H2DRS_MAX_ORDER, H2D_GET_H_ORDER(info_p.max_quad_order), H2D_GET_V_ORDER(info_p.max_quad_order));
    assert_msg(info_aniso.is_empty() || (H2D_GET_H_ORDER(info_aniso.max_quad_order) <= H2DRS_MAX_ORDER && H2D_GET_V_ORDER(info_aniso.max_quad_order) <= H2DRS_MAX_ORDER), "Maximum allowed order of a son of ANISO-candidate is %d but order (H:%d,V:%d) requested.", H2DRS_MAX_ORDER, H2D_GET_H_ORDER(info_aniso.max_quad_order), H2D_GET_V_ORDER(info_aniso.max_quad_order));

    int mode = e->get_mode();

    // select quadrature, obtain integration points and weights
    Quad2D* quad = &g_quad_2d_std;
    quad->set_mode(mode);
    rsln->set_quad_2d(quad);
    double3* gip_points = quad->get_points(H2DRS_INTR_GIP_ORDER);
    int num_gip_points = quad->get_num_points(H2DRS_INTR_GIP_ORDER);

    // everything is done on the reference domain
    rsln->enable_transform(false);

    // obtain reference solution values on all four refined sons
    scalar** rval[H2D_MAX_ELEMENT_SONS];
    Element* base_element = rsln->get_mesh()->get_element(e->id);
    assert(!base_element->active);
    for (int son = 0; son < H2D_MAX_ELEMENT_SONS; son++)
    {
      //set element
      Element* e = base_element->sons[son];
      assert(e != NULL);

      //obtain precalculated values
      rval[son] = precalc_ref_solution(son, rsln, e, H2DRS_INTR_GIP_ORDER);
    }

    //retrieve transformations
    Trf* trfs = NULL;
    int num_noni_trfs = 0;
    if (mode == HERMES_MODE_TRIANGLE) {
      trfs = tri_trf;
      num_noni_trfs = H2D_TRF_TRI_NUM;
    }
    else {
      trfs = quad_trf;
      num_noni_trfs = H2D_TRF_QUAD_NUM;
    }

    // precalculate values of shape functions
    TrfShape empty_shape_vals;
    if (!cached_shape_vals_valid[mode]) {
      precalc_ortho_shapes(gip_points, num_gip_points, trfs, num_noni_trfs, shape_indices[mode], max_shape_inx[mode], cached_shape_ortho_vals[mode]);
      precalc_shapes(gip_points, num_gip_points, trfs, num_noni_trfs, shape_indices[mode], max_shape_inx[mode], cached_shape_vals[mode]);
      cached_shape_vals_valid[mode] = true;

      //issue a warning if ortho values are defined and the selected cand_list might benefit from that but it cannot because elements do not have uniform orders
      if (!warn_uniform_orders && mode == HERMES_MODE_QUAD && !cached_shape_ortho_vals[mode][H2D_TRF_IDENTITY].empty()) {
        warn_uniform_orders = true;
        if (cand_list == H2D_H_ISO || cand_list == H2D_H_ANISO || cand_list == H2D_P_ISO || cand_list == H2D_HP_ISO || cand_list == H2D_HP_ANISO_H) {
          warn_if(!info_h.uniform_orders || !info_aniso.uniform_orders || !info_p.uniform_orders, "Possible inefficiency: %s might be more efficient if the input mesh contains elements with uniform orders strictly.", get_cand_list_str(cand_list));
        }
      }
    }
    TrfShape& svals = cached_shape_vals[mode];
    TrfShape& ortho_svals = cached_shape_ortho_vals[mode];

    //H-candidates
    if (!info_h.is_empty()) {
      Trf* p_trf_identity[1] = { &trfs[H2D_TRF_IDENTITY] };
      std::vector<TrfShapeExp>* p_trf_svals[1] = { &svals[H2D_TRF_IDENTITY] };
      std::vector<TrfShapeExp>* p_trf_ortho_svals[1] = { &ortho_svals[H2D_TRF_IDENTITY] };
      for(int son = 0; son < H2D_MAX_ELEMENT_SONS; son++) {
        scalar **sub_rval[1] = { rval[son] };
        calc_error_cand_element(mode, gip_points, num_gip_points
          , 1, &base_element->sons[son], p_trf_identity, sub_rval
          , p_trf_svals, p_trf_ortho_svals
          , info_h, herr[son]);
      }
    }

    //ANISO-candidates
    if (!info_aniso.is_empty()) {
      const int sons[4][2] = { {0,1}, {3,2}, {0,3}, {1,2} }; //indices of sons for sub-areas
      const int tr[4][2]   = { {6,7}, {6,7}, {4,5}, {4,5} }; //indices of ref. domain transformations for sub-areas
      for(int version = 0; version < 4; version++) { // 2 elements for vertical split, 2 elements for horizontal split
        Trf* sub_trfs[2] = { &trfs[tr[version][0]], &trfs[tr[version][1]] };
        Element* sub_domains[2] = { base_element->sons[sons[version][0]], base_element->sons[sons[version][1]] };
        scalar **sub_rval[2] = { rval[sons[version][0]], rval[sons[version][1]] };
        std::vector<TrfShapeExp>* sub_svals[2] = { &svals[tr[version][0]], &svals[tr[version][1]] };
        std::vector<TrfShapeExp>* sub_ortho_svals[2] = { &ortho_svals[tr[version][0]], &ortho_svals[tr[version][1]] };
        calc_error_cand_element(mode, gip_points, num_gip_points
          , 2, sub_domains, sub_trfs, sub_rval
          , sub_svals, sub_ortho_svals
          , info_aniso, anisoerr[version]);
      }
    }

    //P-candidates
    if (!info_p.is_empty()) {
      Trf* sub_trfs[4] = { &trfs[0], &trfs[1], &trfs[2], &trfs[3] };
      scalar **sub_rval[4] = { rval[0], rval[1], rval[2], rval[3] };
      std::vector<TrfShapeExp>* sub_svals[4] = { &svals[0], &svals[1], &svals[2], &svals[3] };
      std::vector<TrfShapeExp>* sub_ortho_svals[4] = { &ortho_svals[0], &ortho_svals[1], &ortho_svals[2], &ortho_svals[3] };

      calc_error_cand_element(mode, gip_points, num_gip_points
//.........这里部分代码省略.........

/// Preparation of solutions.
Func<scalar>* init_fn(Solution *fu, const int order)
{
  // Sanity checks.
  if (fu == NULL) error("NULL MeshFunction in Func<scalar>*::init_fn().");
  if (fu->get_mesh() == NULL) error("Uninitialized MeshFunction used.");

  ESpaceType space_type = fu->get_space_type();
  ESolutionType sln_type = fu->get_type();

  int nc = fu->get_num_components();
  Quad2D* quad = fu->get_quad_2d();
#ifdef H2D_SECOND_DERIVATIVES_ENABLED
  if (space_type == HERMES_H1_SPACE && sln_type != HERMES_EXACT)
    fu->set_quad_order(order, H2D_FN_ALL);
  else
#endif
    fu->set_quad_order(order);

  double3* pt = quad->get_points(order);
  int np = quad->get_num_points(order);
  Func<scalar>* u = new Func<scalar>(np, nc);

  if (u->nc == 1) {
    u->val = new scalar [np];
    u->dx  = new scalar [np];
    u->dy  = new scalar [np];
#ifdef H2D_SECOND_DERIVATIVES_ENABLED
    if (space_type == HERMES_H1_SPACE && sln_type != HERMES_EXACT)
      u->laplace = new scalar [np];
#endif
    memcpy(u->val, fu->get_fn_values(), np * sizeof(scalar));
    memcpy(u->dx, fu->get_dx_values(), np * sizeof(scalar));
    memcpy(u->dy, fu->get_dy_values(), np * sizeof(scalar));
#ifdef H2D_SECOND_DERIVATIVES_ENABLED
    if (space_type == HERMES_H1_SPACE) {
      if(sln_type == HERMES_SLN) {
        scalar *dxx = fu->get_dxx_values();
        scalar *dyy = fu->get_dyy_values();
        for (int i = 0; i < np; i++)
          u->laplace[i] = dxx[i] + dyy[i];
      }
      else if (sln_type == HERMES_CONST)
        memset(u->laplace, 0, np * sizeof(scalar));
    }
#endif
  }
  else if (u->nc == 2) {
    u->val0 = new scalar [np];
    u->val1 = new scalar [np];
    u->curl = new scalar [np];
    u->div = new scalar [np];

    memcpy(u->val0, fu->get_fn_values(0), np * sizeof(scalar));
    memcpy(u->val1, fu->get_fn_values(1), np * sizeof(scalar));

    scalar *dx1 = fu->get_dx_values(1);
    scalar *dy0 = fu->get_dy_values(0);
    for (int i = 0; i < np; i++) u->curl[i] = dx1[i] - dy0[i];
    
    scalar *dx0 = fu->get_dx_values(0);
    scalar *dy1 = fu->get_dy_values(1);
    for (int i = 0; i < np; i++) u->div[i] = dx0[i] + dy1[i];
  }
  return u;
}