C++ PrecalcShapeset::get_fn_values方法代码示例

static void old_projection(Element* e, int order, double2* proj, double* old[2])
{
  int mo2 = quad2d.get_max_order();
  int np = quad2d.get_num_points(mo2);

  for (unsigned int k = 0; k < e->nvert; k++) // loop over vertices
  {
    // vertex basis functions in all integration points
    double* vd;
    int index_v = ref_map_shapeset.get_vertex_index(k);
    ref_map_pss.set_active_shape(index_v);
    ref_map_pss.set_quad_order(mo2);
    vd = ref_map_pss.get_fn_values();

    for (int m = 0; m < 2; m++)   // part 0 or 1
      for (int j = 0; j < np; j++)
        old[m][j] += proj[k][m] * vd[j];

    for (int ii = 0; ii < order - 1; ii++)
    {
      // edge basis functions in all integration points
      double* ed;
      int index_e = ref_map_shapeset.get_edge_index(k,0,ii+2);
      ref_map_pss.set_active_shape(index_e);
      ref_map_pss.set_quad_order(mo2);
      ed = ref_map_pss.get_fn_values();

      for (int m = 0; m < 2; m++)  //part 0 or 1
        for (int j = 0; j < np; j++)
          old[m][j] += proj[e->nvert + k * (order-1) + ii][m] * ed[j];
    }
  }
}

void test_edge_rotation()
{
  info("Testing edge rotation...");
  int mode = shapeset.get_mode();
  int ne = mode ? 4 : 3;
  
  for (int ori = 0; ori <= 1; ori++)
  {
    for (int order = 0; order <= shapeset.get_max_order(); order++)
    {
      double *e01, *e02, *ee1, *ee2;
      precalc.set_active_shape(shapeset.get_edge_index(0, ori, order));
      precalc.set_quad_order(quad.get_edge_points(0));
      e01 = precalc.get_fn_values(0);
      if (nc > 1) e02 = precalc.get_fn_values(1);
      
      for (int e = 1; e < ne; e++)
      {
        precalc.set_active_shape(shapeset.get_edge_index(e, ori, order));
        precalc.set_quad_order(quad.get_edge_points(e));
        ee1 = precalc.get_fn_values(0);
        if (nc > 1) ee2 = precalc.get_fn_values(1);
        
        int np = quad.get_num_points(quad.get_edge_points(0));
        if (nc == 1)
        {
          for (int i = 0; i < np; i++)
            if (!eq(e01[i], ee1[i]))
            {
              info("order=%d, ori=%d, edge=%d -- not equal to edge 0", order, ori, e);
            }  
        }
        else
        {
          for (int i = 0; i < np; i++)
          {
            double x = rot[mode][e][0][0] * ee1[i] + rot[mode][e][0][1] * ee2[i];
            double y = rot[mode][e][1][0] * ee1[i] + rot[mode][e][1][1] * ee2[i];
            if (!eq(e01[i], x) || !eq(e02[i], y))
            {
              info("order=%d, ori=%d, edge=%d -- not equal to edge 0", order, ori, e);
              printf("x comp: 0-ta %g, %d-ta %g\n", e01[i], e, x);
              printf("y comp: 0-ta %g, %d-ta %g\n\n", e02[i], e, y);
            }  
          }
        }          
      }
    }
  }
}

void RefMap::calc_phys_y(int order)
{
  // transform all y coordinates of the integration points
  int i, j, np = quad_2d->get_num_points(order);
  double* y = cur_node->phys_y[order] = new double[np];
  memset(y, 0, np * sizeof(double));
  ref_map_pss.force_transform(sub_idx, ctm);
  for (i = 0; i < nc; i++)
  {
    ref_map_pss.set_active_shape(indices[i]);
    ref_map_pss.set_quad_order(order);
    double* fn = ref_map_pss.get_fn_values();
    for (j = 0; j < np; j++)
      y[j] += coeffs[i][1] * fn[j];
  }
}

static void calc_bubble_projection(Element* e, Nurbs** nurbs, int order, double2* proj)
{
  ref_map_pss.set_active_element(e);

  int i, j, k;
  int mo2 = quad2d.get_max_order();
  int np = quad2d.get_num_points(mo2);
  int qo = e->is_quad() ? make_quad_order(order, order) : order;
  int nb = ref_map_shapeset.get_num_bubbles(qo);

  AUTOLA_OR(double2, fn, np);
  memset(fn, 0, sizeof(double2) * np);

  double* rhside[2];
  double* old[2];
  for (i = 0; i < 2; i++) {
    rhside[i] = new double[nb];
    old[i] = new double[np];
    memset(rhside[i], 0, sizeof(double) * nb);
    memset(old[i], 0, sizeof(double) * np);
  }

  // compute known part of projection (vertex and edge part)
  old_projection(e, order, proj, old);

  // fn values of both components of nonpolynomial function
  double3* pt = quad2d.get_points(mo2);
  for (j = 0; j < np; j++)  // over all integration points
  {
    double2 a;
    a[0] = ctm.m[0] * pt[j][0] + ctm.t[0];
    a[1] = ctm.m[1] * pt[j][1] + ctm.t[1];
    calc_ref_map(e, nurbs, a[0], a[1], fn[j]);
  }

  double2* result = proj + e->nvert + e->nvert * (order - 1);
  for (k = 0; k < 2; k++)
  {
    for (i = 0; i < nb; i++) // loop over bubble basis functions
    {
      // bubble basis functions in all integration points
      double *bfn;
      int index_i = ref_map_shapeset.get_bubble_indices(qo)[i];
      ref_map_pss.set_active_shape(index_i);
      ref_map_pss.set_quad_order(mo2);
      bfn = ref_map_pss.get_fn_values();

      for (j = 0; j < np; j++) // over all integration points
        rhside[k][i] += pt[j][2] * (bfn[j] * (fn[j][k] - old[k][j]));
    }

    // solve
    if (e->nvert == 3)
      cholsl(bubble_proj_matrix_tri, nb, bubble_tri_p, rhside[k], rhside[k]);
    else
      cholsl(bubble_proj_matrix_quad, nb, bubble_quad_p, rhside[k], rhside[k]);

    for (i = 0; i < nb; i++)
      result[i][k] = rhside[k][i];
  }

  for (i = 0; i < 2; i++) {
    delete [] rhside[i];
    delete [] old[i];
  }
}