C++ Geom::free方法代码示例

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

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

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 Element::calc_area(bool precise_for_curvature)
    {
      // First some basic arithmetics.
      double ax, ay, bx, by;
      ax = vn[1]->x - vn[0]->x;
      ay = vn[1]->y - vn[0]->y;
      bx = vn[2]->x - vn[0]->x;
      by = vn[2]->y - vn[0]->y;

      this->area = 0.5*(ax*by - ay*bx);
      if (is_quad())
      {
        ax = bx; ay = by;
        bx = vn[3]->x - vn[0]->x;
        by = vn[3]->y - vn[0]->y;

        this->area = area + 0.5*(ax*by - ay*bx);
      }

      // Either the basic approximation is fine.
      if (!this->is_curved() || !precise_for_curvature)
        return;
      // Or we want to capture the curvature precisely.
      else
      {
        // Utility data.
        RefMap refmap_curv;
        RefMap refmap_straight;
        double3* tan;

        double x_center, y_center;
        this->get_center(x_center, y_center);

        for (int isurf = 0; isurf < this->nvert; isurf++)
        {
          // 0 - prepare data structures.
          int eo = g_quad_2d_std.get_edge_points(isurf, this->get_mode() == HERMES_MODE_TRIANGLE ? g_max_tri : g_max_quad, this->get_mode());
          int np = g_quad_2d_std.get_num_points(eo, this->get_mode());
          double* x_curv = new double[np];
          double* y_curv = new double[np];
          double* x_straight = new double[np];
          double* y_straight = new double[np];

          // 1 - get the x,y coordinates for the curved element.
          refmap_curv.set_active_element(this);
          Geom<double>* geometry = init_geom_surf(&refmap_curv, isurf, this->en[isurf]->marker, eo, tan);
          memcpy(x_curv, geometry->x, np*sizeof(double));
          memcpy(y_curv, geometry->y, np*sizeof(double));
          geometry->free();
          delete geometry;

          // 2. - act if there was no curvature
          CurvMap* cm_temp = this->cm;
          this->cm = nullptr;
          refmap_straight.set_active_element(this);
          geometry = init_geom_surf(&refmap_straight, isurf, this->en[isurf]->marker, eo, tan);
          memcpy(x_straight, geometry->x, np*sizeof(double));
          memcpy(y_straight, geometry->y, np*sizeof(double));
          geometry->free();
          delete geometry;

          // 3. - compare the two, get the updated area.
          double previous_distance;
          for (int i = 0; i < np; i++)
          {
            // Distance between the curved and straight edges.
            double distance_i = std::sqrt(std::pow(x_straight[i] - x_curv[i], 2.0) + std::pow(y_straight[i] - y_curv[i], 2.0));

            // Add to- or Subtract from- the area (depends on the curvature and we shall decide based on distance from the element center).
            double distance_from_center_curved = std::pow(x_center - x_curv[i], 2.0) + std::pow(y_center - y_curv[i], 2.0);
            double distance_from_center_straight = std::pow(x_center - x_straight[i], 2.0) + std::pow(y_center - y_straight[i], 2.0);
            bool add = distance_from_center_curved > distance_from_center_straight;

            // Calculate now the area delta.
            // It depends on the integration point number etc.
            double area_delta;
            if (i == 0)
            {
              double distance_along_edge = std::sqrt(std::pow(x_straight[i] - this->vn[isurf]->x, 2.0) + std::pow(y_straight[i] - this->vn[isurf]->y, 2.0));
              area_delta = distance_i * distance_along_edge * 0.5;
            }
            if (i > 0 && i < np - 1)
            {
              double distance_along_edge = std::sqrt(std::pow(x_straight[i] - x_straight[i - 1], 2.0) + std::pow(y_straight[i] - y_straight[i - 1], 2.0));
              area_delta = 0.5*(distance_i + previous_distance) * distance_along_edge;
            }
            if (i == np - 1)
            {
              double distance_along_edge = std::sqrt(std::pow(x_straight[i] - this->vn[(isurf + 1) % this->nvert]->x, 2.0) + std::pow(y_straight[i] - this->vn[(isurf + 1) % this->nvert]->y, 2.0));
              area_delta = distance_i * distance_along_edge * 0.5;
            }

            if (add)
              area += area_delta;
            else
              area -= area_delta;

            previous_distance = distance_i;
          }

//.........这里部分代码省略.........