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

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 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
//.........这里部分代码省略.........

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
//.........这里部分代码省略.........

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));
//.........这里部分代码省略.........