C++ Quad2D类代码示例

void FilterFluxDensity::precalculate(int order, int mask)
{
  Quad2D* quad = quads[cur_quad];
  int np = quad->get_num_points(order, this->get_active_element()->get_mode());
  Node* node = new_node(H2D_FN_DEFAULT, np);

  sln[0]->set_quad_order(order, H2D_FN_VAL | H2D_FN_DX | H2D_FN_DY);
  sln[1]->set_quad_order(order, H2D_FN_VAL | H2D_FN_DX | H2D_FN_DY);

  double *dudx1, *dudy1, *dudx2, *dudy2;
  sln[0]->get_dx_dy_values(dudx1, dudy1);
  sln[1]->get_dx_dy_values(dudx2, dudy2);
  double *uval1 = sln[0]->get_fn_values();
  double *uval2 = sln[1]->get_fn_values();
  update_refmap();
  double *x = refmap->get_phys_x(order);

  for (int i = 0; i < np; i++)
  {
    node->values[0][0][i] = std::sqrt(sqr(dudy1[i]) + sqr(dudy2[i]) +
      sqr(dudx1[i] + ((x[i] > 1e-10) ? uval1[i] / x[i] : 0.0)) +
      sqr(dudx2[i] + ((x[i] > 1e-10) ? uval2[i] / x[i] : 0.0)));
  }

  if(nodes->present(order)) 
  {
    assert(nodes->get(order) == cur_node);
    ::free(nodes->get(order));
  }
  nodes->add(node, order);
  cur_node = node;
}

// Calculates maximum of a given function, including its coordinates.
Extremum get_peak(MeshFunction *sln)
{
  Quad2D* quad = &g_quad_2d_std;
  sln->set_quad_2d(quad);
  Element* e;
  Mesh* mesh = sln->get_mesh();
  
  scalar peak = 0.0;
  double pos_x = 0.0;
  double pos_y = 0.0;
  
  for_all_active_elements(e, mesh)
  {
    update_limit_table(e->get_mode());
    sln->set_active_element(e);
    RefMap* ru = sln->get_refmap();
    int o = sln->get_fn_order() + ru->get_inv_ref_order();
    limit_order(o);
    sln->set_quad_order(o, H2D_FN_VAL);
    scalar *uval = sln->get_fn_values();
    int np = quad->get_num_points(o);
    double* x = ru->get_phys_x(o);
    double* y = ru->get_phys_y(o);
    
    for (int i = 0; i < np; i++)
      if (uval[i] > peak) {
        peak = uval[i];
        pos_x = x[i];
        pos_y = y[i];
      }
  }

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

//------------------------------------------------------------------------------
// 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()

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

// Initialize edge marker, coordinates, tangent and normals
Geom<double>* init_geom_surf(RefMap *rm, SurfPos* surf_pos, const int order)
{
  Geom<double>* e = new Geom<double>;
  e->edge_marker = surf_pos->marker;
  e->elem_marker = rm->get_active_element()->marker;
  e->diam = rm->get_active_element()->get_diameter();
  e->id = rm->get_active_element()->en[surf_pos->surf_num]->id;
  e->x = rm->get_phys_x(order);
  e->y = rm->get_phys_y(order);
  double3 *tan;
  tan = rm->get_tangent(surf_pos->surf_num, order);

  Quad2D* quad = rm->get_quad_2d();
  int np = quad->get_num_points(order);
  e->tx = new double [np];
  e->ty = new double [np];
  e->nx = new double [np];
  e->ny = new double [np];
  for (int i = 0; i < np; i++)
  {
    e->tx[i] = tan[i][0];  e->ty[i] =   tan[i][1];
    e->nx[i] = tan[i][1];  e->ny[i] = - tan[i][0];
  }
    e->orientation = rm->get_active_element()->get_edge_orientation(surf_pos->surf_num);
	return e;
}

void SimpleFilter::precalculate(int order, int mask)
{
  if (mask & (H2D_FN_DX | H2D_FN_DY | H2D_FN_DXX | H2D_FN_DYY | H2D_FN_DXY))
    error("Filter not defined for derivatives.");

  Quad2D* quad = quads[cur_quad];
  int np = quad->get_num_points(order);
  Node* node = new_node(H2D_FN_VAL, np);

  // precalculate all solutions
  for (int i = 0; i < num; i++)
    sln[i]->set_quad_order(order, item[i]);

  for (int j = 0; j < num_components; j++)
  {
    // obtain corresponding tables
    scalar* tab[10];
    for (int i = 0; i < num; i++)
    {
      int a = 0, b = 0, mask = item[i];
      if (mask >= 0x40) { a = 1; mask >>= 6; }
      while (!(mask & 1)) { mask >>= 1; b++; }
      tab[i] = sln[i]->get_values(num_components == 1 ? a : j, b);
      if (tab[i] == NULL) error("Value of 'item%d' is incorrect in filter definition.", i+1);
    }

		Tuple<scalar*> values;
		for(int i = 0; i < this->num; i++)
			values.push_back(tab[i]);

    // apply the filter
		filter_fn(np, values, node->values[j][0]);
  }

    void SimpleFilter<Scalar>::precalculate(int order, int mask)
    {
      if(mask & (H2D_FN_DX | H2D_FN_DY | H2D_FN_DXX | H2D_FN_DYY | H2D_FN_DXY))
        throw Hermes::Exceptions::Exception("Filter not defined for derivatives.");

      Quad2D* quad = this->quads[this->cur_quad];
      int np = quad->get_num_points(order, this->element->get_mode());
      struct Function<Scalar>::Node* node = this->new_node(H2D_FN_VAL, np);

      // precalculate all solutions
      for (int i = 0; i < this->num; i++)
        this->sln[i]->set_quad_order(order, item[i]);

      for (int j = 0; j < this->num_components; j++)
      {
        // obtain corresponding tables
        Scalar* tab[H2D_MAX_COMPONENTS];
        for (int i = 0; i < this->num; i++)
        {
          int a = 0, b = 0, mask = item[i];
          if(mask >= 0x40) { a = 1; mask >>= 6; }
          while (!(mask & 1)) { mask >>= 1; b++; }
          tab[i] = this->sln[i]->get_values(this->num_components == 1 ? a : j, b);
          if(tab[i] == nullptr) throw Hermes::Exceptions::Exception("Value of 'item%d' is incorrect in filter definition.", i + 1);
        }

        Hermes::vector<Scalar*> values;
        for(int i = 0; i < this->num; i++)
          values.push_back(tab[i]);

        // apply the filter
        filter_fn(np, values, node->values[j][0]);
      }

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.
}

// Actual evaluation of volume matrix form (calculates integral)
scalar FeProblem::eval_form(WeakForm::MatrixFormVol *mfv, Tuple<Solution *> u_ext, PrecalcShapeset *fu, PrecalcShapeset *fv, RefMap *ru, RefMap *rv)
{
  // determine the integration order
  int inc = (fu->get_num_components() == 2) ? 1 : 0;
  AUTOLA_OR(Func<Ord>*, oi, wf->neq);
  for (int i = 0; i < wf->neq; i++) oi[i] = init_fn_ord(u_ext[i]->get_fn_order() + inc);
  Func<Ord>* ou = init_fn_ord(fu->get_fn_order() + inc);
  Func<Ord>* ov = init_fn_ord(fv->get_fn_order() + inc);
  ExtData<Ord>* fake_ext = init_ext_fns_ord(mfv->ext);

  double fake_wt = 1.0;
  Geom<Ord>* fake_e = init_geom_ord();
  Ord o = mfv->ord(1, &fake_wt, oi, ou, ov, fake_e, fake_ext);
  int order = ru->get_inv_ref_order();
  order += o.get_order();
  limit_order_nowarn(order);

  for (int i = 0; i < wf->neq; i++) {  oi[i]->free_ord(); delete oi[i]; }
  ou->free_ord(); delete ou;
  ov->free_ord(); delete ov;
  delete fake_e;
  fake_ext->free_ord(); delete fake_ext;

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

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

  // 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, order);
  Func<double>* u = get_fn(fu, ru, order);
  Func<double>* v = get_fn(fv, rv, order);
  ExtData<scalar>* ext = init_ext_fns(mfv->ext, rv, order);

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

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

void ViewScalarFilter::precalculate(int order, int mask)
{
    Quad2D* quad = quads[cur_quad];
    int np = quad->get_num_points(order);
    node = new_node(H2D_FN_DEFAULT, np);

    if (sln[0])
    {
        sln[0]->set_quad_order(order, H2D_FN_VAL | H2D_FN_DX | H2D_FN_DY);
        sln[0]->get_dx_dy_values(dudx1, dudy1);
        value1 = sln[0]->get_fn_values();
    }

    if (num >= 2 && sln[1])
    {
        sln[1]->set_quad_order(order, H2D_FN_VAL | H2D_FN_DX | H2D_FN_DY);
        sln[1]->get_dx_dy_values(dudx2, dudy2);
        value2 = sln[1]->get_fn_values();
    }

    if (num >= 3 && sln[2])
    {
        sln[2]->set_quad_order(order, H2D_FN_VAL | H2D_FN_DX | H2D_FN_DY);
        sln[2]->get_dx_dy_values(dudx3, dudy3);
        value3 = sln[2]->get_fn_values();
    }

    update_refmap();

    x = refmap->get_phys_x(order);
    y = refmap->get_phys_y(order);
    Element *e = refmap->get_active_element();

    material = Util::scene()->labels[atoi(mesh->get_element_markers_conversion().get_user_marker(e->marker).c_str())]->material;

    for (int i = 0; i < np; i++)
    {
        calculateVariable(i);
    }

    if (nodes->present(order))
    {
        assert(nodes->get(order) == cur_node);
        ::free(nodes->get(order));
    }
    nodes->add(node, order);
    cur_node = node;
}

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

/// Calculate number of negative solution values.
int get_num_of_neg(MeshFunction *sln)
{
	Quad2D* quad = &g_quad_2d_std;
  sln->set_quad_2d(quad);
  Element* e;
  Mesh* mesh = sln->get_mesh();

  int n = 0;

  for_all_active_elements(e, mesh)
  {
    update_limit_table(e->get_mode());
    sln->set_active_element(e);
    RefMap* ru = sln->get_refmap();
    int o = sln->get_fn_order() + ru->get_inv_ref_order();
    limit_order(o);
    sln->set_quad_order(o, H2D_FN_VAL);
    scalar *uval = sln->get_fn_values();
    int np = quad->get_num_points(o);

		for (int i = 0; i < np; i++)
			if (uval[i] < -1e-12)
				n++;
  }