C++ PrecalcShapeset类代码示例

void DiscreteProblem::precalc_equi_coefs()
{
  int i, m;
  memset(equi, 0, sizeof(double) * ndofs);
  verbose("Precalculating equilibration coefficients...");

  RefMap refmap;
  AsmList al;
  Element* e;

  for (m = 0; m < neq; m++)
  {
    PrecalcShapeset* fu = pss[m];
    BiForm* bf = biform[m] + m;
    Mesh* mesh = spaces[m]->get_mesh();

    for_all_active_elements(e, mesh)
    {
      update_limit_table(e->get_mode());
      fu->set_active_element(e);
      refmap.set_active_element(e);

      spaces[m]->get_element_assembly_list(e, &al);
      for (i = 0; i < al.cnt; i++)
      {
        if (al.dof[i] < 0) continue;
        fu->set_active_shape(al.idx[i]);
        scalar sy = 0.0, un = 0.0;
        if (bf->unsym) un = bf->unsym(fu, fu, &refmap, &refmap);
        if (bf->sym)   sy = bf->sym  (fu, fu, &refmap, &refmap);
        #ifndef COMPLEX
        equi[al.dof[i]] += (sy + un) * sqr(al.coef[i]);
        #else
        equi[al.dof[i]] += 0;//std::norm(sy + un) * sqr(al.coef[i]);
        #endif
      }
    }
  }

int main(int argc, char* argv[])
{
  hermes2d_initialize(&argc, argv);

  info("SHAPESET TESTER");
  info("num_components = %d", shapeset.get_num_components());
  info("max_order = %d", shapeset.get_max_order());
  
  precalc.set_quad_2d(&quad);
  for (int mode = 0; mode <= 1; mode++)
  {
    shapeset.set_mode(mode);
    quad.set_mode(mode);
    precalc.set_mode(mode);
    info(mode ? "\nTESTING QUADS\n" : "\nTESTING TRIANGLES\n");
    
    //test_orders(&shapeset);
    test_edge_rotation();
    //test_edge_orientation(&shapeset);
    //test_num_bubbles(&shapeset);
  }
  
  //info("\nALL OK!\n");

  printf("\n");
  hermes2d_finalize();
  return 0;
}

static void precalculate_cholesky_projection_matrices_bubble()
{
  // *** triangles ***
  ref_map_pss.set_mode(MODE_TRIANGLE);
  int order = ref_map_shapeset.get_max_order();

  // calculate projection matrix of maximum order
  int nb = ref_map_shapeset.get_num_bubbles(order);
  int* indices = ref_map_shapeset.get_bubble_indices(order);
  bubble_proj_matrix_tri = calculate_bubble_projection_matrix(nb, indices);

  // cholesky factorization of the matrix
  bubble_tri_p = new double[nb];
  choldc(bubble_proj_matrix_tri, nb, bubble_tri_p);

  // *** quads ***
  ref_map_pss.set_mode(MODE_QUAD);
  order = ref_map_shapeset.get_max_order();
  order = make_quad_order(order, order);

  // calculate projection matrix of maximum order
  nb = ref_map_shapeset.get_num_bubbles(order);
  indices = ref_map_shapeset.get_bubble_indices(order);
  bubble_proj_matrix_quad = calculate_bubble_projection_matrix(nb, indices);

  // cholesky factorization of the matrix
  bubble_quad_p = new double[nb];
  choldc(bubble_proj_matrix_quad, nb, bubble_quad_p);
}

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];
  }
}

void RefMap::set_active_element(Element* e)
{
  if (e != element) free();

  ref_map_pss.set_active_element(e);
  quad_2d->set_mode(e->get_mode());
  num_tables = quad_2d->get_num_tables();
  assert(num_tables <= H2D_MAX_TABLES);

  if (e == element) return;
  element = e;

  reset_transform();
  update_cur_node();

  is_const = !element->is_curved() &&
             (element->is_triangle() || is_parallelogram());

  // prepare the shapes and coefficients of the reference map
  int j, k = 0;
  for (unsigned int i = 0; i < e->nvert; i++)
    indices[k++] = ref_map_shapeset.get_vertex_index(i);

  // straight-edged element
  if (e->cm == NULL)
  {
    for (unsigned int i = 0; i < e->nvert; i++)
    {
      lin_coeffs[i][0] = e->vn[i]->x;
      lin_coeffs[i][1] = e->vn[i]->y;
    }
    coeffs = lin_coeffs;
    nc = e->nvert;
  }
  else // curvilinear element - edge and bubble shapes
  {
    int o = e->cm->order;
    for (unsigned int i = 0; i < e->nvert; i++)
      for (j = 2; j <= o; j++)
        indices[k++] = ref_map_shapeset.get_edge_index(i, 0, j);

    if (e->is_quad()) o = H2D_MAKE_QUAD_ORDER(o, o);
    memcpy(indices + k, ref_map_shapeset.get_bubble_indices(o),
           ref_map_shapeset.get_num_bubbles(o) * sizeof(int));

    coeffs = e->cm->coeffs;
    nc = e->cm->nc;
  }

  // calculate the order of the inverse reference map
  if (element->iro_cache == -1 && quad_2d->get_max_order() > 1)
  {
    element->iro_cache = is_const ? 0 : calc_inv_ref_order();
  }
  inv_ref_order = element->iro_cache;

  // constant inverse reference map
  if (is_const) calc_const_inv_ref_map(); else const_jacobian = 0.0;
}

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

  int i, j, k;
  int mo1 = quad1d.get_max_order();
  int np = quad1d.get_num_points(mo1);
  int ne = order - 1;
  int mode = e->get_mode();

  assert(np <= 15 && ne <= 10);
  double2 fn[15];
  double rhside[2][10];
  memset(fn, 0, sizeof(double2) * np);
  memset(rhside[0], 0, sizeof(double) * ne);
  memset(rhside[1], 0, sizeof(double) * ne);

  double a_1, a_2, b_1, b_2;
  a_1 = ctm.m[0] * ref_vert[mode][edge][0] + ctm.t[0];
  a_2 = ctm.m[1] * ref_vert[mode][edge][1] + ctm.t[1];
  b_1 = ctm.m[0] * ref_vert[mode][e->next_vert(edge)][0] + ctm.t[0];
  b_2 = ctm.m[1] * ref_vert[mode][e->next_vert(edge)][1] + ctm.t[1];

  // values of nonpolynomial function in two vertices
  double2 fa, fb;
  calc_ref_map(e, nurbs, a_1, a_2, fa);
  calc_ref_map(e, nurbs, b_1, b_2, fb);

  double2* pt = quad1d.get_points(mo1);
  for (j = 0; j < np; j++) // over all integration points
  {
    double2 x, v;
    double t = pt[j][0];
    edge_coord(e, edge, t, x, v);
    calc_ref_map(e, nurbs, x[0], x[1], fn[j]);

    for (k = 0; k < 2; k++)
      fn[j][k] = fn[j][k] - (fa[k] + (t+1)/2.0 * (fb[k] - fa[k]));
  }

  double2* result = proj + e->nvert + edge * (order - 1);
  for (k = 0; k < 2; k++)
  {
    for (i = 0; i < ne; i++)
    {
      for (j = 0; j < np; j++)
      {
        double t = pt[j][0];
        double fi = lob[i+2](t);
        rhside[k][i] += pt[j][1] * (fi * fn[j][k]);
      }
    }
    // solve
    cholsl(edge_proj_matrix, ne, edge_p, rhside[k], rhside[k]);
    for (i = 0; i < ne; i++)
      result[i][k] = rhside[k][i];
  }
}

void RefMap::calc_inv_ref_map(int order)
{
  assert(quad_2d != NULL);
  int i, j, np = quad_2d->get_num_points(order);

  // construct jacobi matrices of the direct reference map for all integration points

  AUTOLA_OR(double2x2, m, np);
  memset(m, 0, m.size);
  ref_map_pss.force_transform(sub_idx, ctm);
  for (i = 0; i < nc; i++)
  {
    double *dx, *dy;
    ref_map_pss.set_active_shape(indices[i]);
    ref_map_pss.set_quad_order(order);
    ref_map_pss.get_dx_dy_values(dx, dy);
    for (j = 0; j < np; j++)
    {
      m[j][0][0] += coeffs[i][0] * dx[j];
      m[j][0][1] += coeffs[i][0] * dy[j];
      m[j][1][0] += coeffs[i][1] * dx[j];
      m[j][1][1] += coeffs[i][1] * dy[j];
    }
  }

  // calculate the jacobian and inverted matrix
  double trj = get_transform_jacobian();
  double2x2* irm = cur_node->inv_ref_map[order] = new double2x2[np];
  double* jac = cur_node->jacobian[order] = new double[np];
  for (i = 0; i < np; i++)
  {
    jac[i] = (m[i][0][0] * m[i][1][1] - m[i][0][1] * m[i][1][0]);
    double ij = 1.0 / jac[i];
    error_if(!finite(ij), "1/jac[%d] is infinity when calculating inv. ref. map for order %d (jac=%g)", i, order);
    assert_msg(ij == ij, "1/jac[%d] is NaN when calculating inv. ref. map for order %d (jac=%g)", i, order);

    // invert and transpose the matrix
    irm[i][0][0] =  m[i][1][1] * ij;
    irm[i][0][1] = -m[i][1][0] * ij;
    irm[i][1][0] = -m[i][0][1] * ij;
    irm[i][1][1] =  m[i][0][0] * ij;

    jac[i] *= trj;
  }
}

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);
            }  
          }
        }          
      }
    }
  }
}

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 RefMap::calc_second_ref_map(int order)
{
  assert(quad_2d != NULL);
  int i, j, np = quad_2d->get_num_points(order);

  AUTOLA_OR(double3x2, k, np);
  memset(k, 0, k.size);
  ref_map_pss.force_transform(sub_idx, ctm);
  for (i = 0; i < nc; i++)
  {
    double *dxy, *dxx, *dyy;
    ref_map_pss.set_active_shape(indices[i]);
    ref_map_pss.set_quad_order(order, H2D_FN_ALL);
    dxx = ref_map_pss.get_dxx_values();
    dyy = ref_map_pss.get_dyy_values();
    dxy = ref_map_pss.get_dxy_values();
    for (j = 0; j < np; j++)
    {
      k[j][0][0] += coeffs[i][0] * dxx[j];
      k[j][0][1] += coeffs[i][1] * dxx[j];
      k[j][1][0] += coeffs[i][0] * dxy[j];
      k[j][1][1] += coeffs[i][1] * dxy[j];
      k[j][2][0] += coeffs[i][0] * dyy[j];
      k[j][2][1] += coeffs[i][1] * dyy[j];
    }
  }

  double3x2* mm = cur_node->second_ref_map[order] = new double3x2[np];
  double2x2* m = get_inv_ref_map(order);
  for (j = 0; j < np; j++)
  {
    double a, b;
    // coefficients in second derivative with respect to xx
    a = sqr(m[j][0][0])*k[j][0][0] + 2*m[j][0][0]*m[j][0][1]*k[j][1][0] + sqr(m[j][0][1])*k[j][2][0];
    b = sqr(m[j][0][0])*k[j][0][1] + 2*m[j][0][0]*m[j][0][1]*k[j][1][1] + sqr(m[j][0][1])*k[j][2][1];
    mm[j][0][0] = -(a * m[j][0][0] + b * m[j][1][0]); // du/dx
    mm[j][0][1] = -(a * m[j][0][1] + b * m[j][1][1]); // du/dy

    // coefficients in second derivative with respect to xy
    a = m[j][0][0]*m[j][1][0]*k[j][0][0] + (m[j][0][1]*m[j][1][0] + m[j][0][0]*m[j][1][1])*k[j][1][0] + m[j][0][1]*m[j][1][1]*k[j][2][0];
    b = m[j][0][0]*m[j][1][0]*k[j][0][1] + (m[j][0][1]*m[j][1][0] + m[j][0][0]*m[j][1][1])*k[j][1][1] + m[j][0][1]*m[j][1][1]*k[j][2][1];
    mm[j][1][0] = -(a * m[j][0][0] + b * m[j][1][0]); // du/dx
    mm[j][1][1] = -(a * m[j][0][1] + b * m[j][1][1]); // du/dy

    // coefficients in second derivative with respect to yy
    a = sqr(m[j][1][0])*k[j][0][0] + 2*m[j][1][0]*m[j][1][1]*k[j][1][0] + sqr(m[j][1][1])*k[j][2][0];
    b = sqr(m[j][1][0])*k[j][0][1] + 2*m[j][1][0]*m[j][1][1]*k[j][1][1] + sqr(m[j][1][1])*k[j][2][1];
    mm[j][2][0] = -(a * m[j][0][0] + b * m[j][1][0]); // du/dx
    mm[j][2][1] = -(a * m[j][0][1] + b * m[j][1][1]); // du/dy
  }
}

void CurvMap::update_refmap_coefs(Element* e)
{
  ref_map_pss.set_quad_2d(&quad2d);
  //ref_map_pss.set_active_element(e);

  // calculation of projection matrices
  if (edge_proj_matrix == NULL) precalculate_cholesky_projection_matrix_edge();
  if (bubble_proj_matrix_tri == NULL) precalculate_cholesky_projection_matrices_bubble();

  ref_map_pss.set_mode(e->get_mode());
  ref_map_shapeset.set_mode(e->get_mode());

  // allocate projection coefficients
  int nv = e->nvert;
  int ne = order - 1;
  int qo = e->is_quad() ? make_quad_order(order, order) : order;
  int nb = ref_map_shapeset.get_num_bubbles(qo);
  nc = nv + nv*ne + nb;
  if (coefs != NULL) delete [] coefs;
  coefs = new double2[nc];

  // WARNING: do not change the format of the array 'coefs'. If it changes,
  // RefMap::set_active_element() has to be changed too.

  Nurbs** nurbs;
  if (toplevel == false)
  {
    ref_map_pss.set_active_element(e);
    ref_map_pss.set_transform(part);
    nurbs = parent->cm->nurbs;
  }
  else
  {
    ref_map_pss.reset_transform();
    nurbs = e->cm->nurbs;
  }
  ctm = *(ref_map_pss.get_ctm());
  ref_map_pss.reset_transform(); // fixme - do we need this?

  // calculation of new projection coefficients
  ref_map_projection(e, nurbs, order, coefs);
}

void RefMap::set_quad_2d(Quad2D* quad_2d)
{
  free();
  this->quad_2d = quad_2d;
  ref_map_pss.set_quad_2d(quad_2d);
}

void RefMap::calc_tangent(int edge, int eo)
{
  int i, j;
  int np = quad_2d->get_num_points(eo);
  double3* tan = cur_node->tan[edge] = new double3[np];
  int a = edge, b = element->next_vert(edge);

  if (!element->is_curved())
  {
    // straight edges: the tangent at each point is just the edge length
    tan[0][0] = element->vn[b]->x - element->vn[a]->x;
    tan[0][1] = element->vn[b]->y - element->vn[a]->y;
    tan[0][2] = sqrt(sqr(tan[0][0]) + sqr(tan[0][1]));
    double inorm = 1.0 / tan[0][2];
    tan[0][0] *= inorm;
    tan[0][1] *= inorm;
    tan[0][2] *= (edge == 0 || edge == 2) ? ctm->m[0] : ctm->m[1];

    for (i = 1; i < np; i++)
      memcpy(tan+i, tan, sizeof(double3));
  }
  else
  {
    // construct jacobi matrices of the direct reference map at integration points along the edge
    static double2x2 m[15];
    assert(np <= 15);
    memset(m, 0, np*sizeof(double2x2));
    ref_map_pss.force_transform(sub_idx, ctm);
    for (i = 0; i < nc; i++)
    {
      double *dx, *dy;
      ref_map_pss.set_active_shape(indices[i]);
      ref_map_pss.set_quad_order(eo);
      ref_map_pss.get_dx_dy_values(dx, dy);
      for (j = 0; j < np; j++)
      {
        m[j][0][0] += coeffs[i][0] * dx[j];
        m[j][0][1] += coeffs[i][0] * dy[j];
        m[j][1][0] += coeffs[i][1] * dx[j];
        m[j][1][1] += coeffs[i][1] * dy[j];
      }
    }

    // multiply them by the vector of the reference edge
    double2* v1 = ref_map_shapeset.get_ref_vertex(a);
    double2* v2 = ref_map_shapeset.get_ref_vertex(b);
    double ex = (*v2)[0] - (*v1)[0];
    double ey = (*v2)[1] - (*v1)[1];
    for (i = 0; i < np; i++)
    {
      double3& t = tan[i];
      t[0] = m[i][0][0]*ex + m[i][0][1]*ey;
      t[1] = m[i][1][0]*ex + m[i][1][1]*ey;
      t[2] = sqrt(sqr(t[0]) + sqr(t[1]));
      double inorm = 1.0 / t[2];
      t[0] *= inorm;
      t[1] *= inorm;
      t[2] *= (edge == 0 || edge == 2) ? ctm->m[0] : ctm->m[1];
    }
  }
}

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];
  }
}