C++ Traverse::finish方法代码示例

/// Calculates the absolute error between sln1 and sln2 using function fn
double calc_abs_error(double (*fn)(MeshFunction*, MeshFunction*, RefMap*, RefMap*), MeshFunction* sln1,
                      MeshFunction* sln2)
{
  // sanity checks
  if (fn == NULL) error("error norm function is NULL in calc_abs_error().");
  if (sln1 == NULL) error("sln1 is NULL in calc_abs_error().");
  if (sln2 == NULL) error("sln2 is NULL in calc_abs_error().");

  Quad2D* quad = &g_quad_2d_std;
  sln1->set_quad_2d(quad);
  sln2->set_quad_2d(quad);

  Mesh* meshes[2] = { sln1->get_mesh(), sln2->get_mesh() };
  Transformable* tr[2] = { sln1, sln2 };
  Traverse trav;
  trav.begin(2, meshes, tr);

  double error = 0.0;
  Element** ee;
  while ((ee = trav.get_next_state(NULL, NULL)) != NULL)
  {
    update_limit_table(ee[0]->get_mode());

    RefMap* ru = sln1->get_refmap();
    RefMap* rv = sln2->get_refmap();

    error += fn(sln1, sln2, ru, rv);
  }
  trav.finish();
  return sqrt(error);
}

/// Calculates the absolute error between sln1 and sln2 using function fn
double calc_error(double (*fn)(MeshFunction*, MeshFunction*, int, QuadPt3D*), MeshFunction *sln1, MeshFunction *sln2) {
	_F_
	Mesh *meshes[2] = { sln1->get_mesh(), sln2->get_mesh() };
	Transformable *tr[2] = { sln1, sln2 };
	Traverse trav;
	trav.begin(2, meshes, tr);

	double error = 0.0;
	Element **ee;
	while ((ee = trav.get_next_state(NULL, NULL)) != NULL) {
		ElementMode3D mode = ee[0]->get_mode();

		RefMap *ru = sln1->get_refmap();
		Ord3 order = max(sln1->get_fn_order(), sln2->get_fn_order()) + ru->get_inv_ref_order();
		order.limit();

		Quad3D *quad = get_quadrature(mode);
		int np = quad->get_num_points(order);
		QuadPt3D *pt = quad->get_points(order);

		error += fn(sln1, sln2, np, pt);
	}
	trav.finish();

	return error > H3D_TINY ? sqrt(error) : error;		// do not ruin the precision by taking the sqrt
}

void Filter::init() {
	_F_
	// construct the union mesh, if necessary
	Mesh *meshes[4] = {
		sln[0]->get_mesh(),
		(num >= 2) ? sln[1]->get_mesh() : NULL,
		(num >= 3) ? sln[2]->get_mesh() : NULL,
		(num >= 4) ? sln[3]->get_mesh() : NULL
	};

	mesh = meshes[0];
	unimesh = false;

	// FIXME: better detection of same meshes
	for (int i = 1; i < num; i++)
		if (meshes[0] != meshes[1]) unimesh = true;

	if (unimesh) {
		Traverse trav;
		trav.begin(num, meshes);
		mesh = new Mesh;
		MEM_CHECK(mesh);
		unidata = trav.construct_union_mesh(mesh);
		trav.finish();
	}

	refmap->set_mesh(mesh);

	// misc init
	num_components = 1;
	memset(tables, 0, sizeof(tables));
	memset(sln_sub, 0, sizeof(sln_sub));
}

void DiscreteProblem::precalculate_sparse_structure(Page** pages)
{
  int i, j, m, n;
  AsmList* al = new AsmList[neq];

  // init multi-mesh traversal
  Mesh** meshes = new Mesh*[neq];
  for (i = 0; i < neq; i++)
    meshes[i] = spaces[i]->get_mesh();
  Traverse trav;
  trav.begin(neq, meshes);

  // loop through all triangles
  Element** e;
  while ((e = trav.get_next_state(NULL, NULL)) != NULL)
  {
    // obtain assembly lists for the element at all spaces
    for (i = 0; i < neq; i++)
      spaces[i]->get_element_assembly_list(e[i], al + i);
      // todo: neziskavat znova, pokud se element nezmenil

    // go through all equation-blocks of the local stiffness matrix
    for (m = 0; m < neq; m++)
    {
      for (n = 0; n < neq; n++)
      {
        BiForm* bf = biform[m] + n;
        if (bf->sym == NULL && bf->unsym == NULL && bf->surf == NULL) continue;

        // pretend assembling of the element stiffness matrix
        for (j = 0; j < al[n].cnt; j++)
        {
          // skip dirichlet dofs in 'j'
          if (al[n].dof[j] < 0) continue;

          for (i = 0; i < al[m].cnt; i++)
          {
            // skip dirichlet dofs in 'i'
            if (al[m].dof[i] < 0) continue;

            // register the corresponding nonzero matrix element
            page_add_ij(pages, al[n].dof[j], al[m].dof[i]); // column-oriented (UMFPACK)
            //page_add_ij(pages, al[m].dof[i], al[n].dof[j]); // row-oriented (PETSc)
          }
        }
      }
    }
  }

  trav.finish();
  delete [] meshes;
  delete [] al;
}

    void Filter<Scalar>::init()
    {
      // construct the union mesh, if necessary
      const Mesh* meshes[10];
      for(int i = 0; i < this->num; i++)
        meshes[i] = this->sln[i]->get_mesh();
      this->mesh = meshes[0];
      unimesh = false;

      for (int i = 1; i < num; i++)

      {
        if(meshes[i] == NULL)
        {
          this->warn("You may be initializing a Filter with Solution that is missing a Mesh.");
          throw Hermes::Exceptions::Exception("this->meshes[%d] is NULL in Filter<Scalar>::init().", i);
        }
        if(meshes[i]->get_seq() != this->mesh->get_seq())
        {
          unimesh = true;
          break;
        }
      }

      if(unimesh)
      {
        Traverse trav;
        trav.begin(num, meshes);
        this->mesh = new Mesh;
        unidata = trav.construct_union_mesh(const_cast<Hermes2D::Mesh*>(this->mesh));
        trav.finish();
      }

      // misc init
      this->num_components = 1;
      this->order = 0;

      memset(sln_sub, 0, sizeof(sln_sub));
      set_quad_2d(&g_quad_2d_std);

      this->deleteSolutions = false;
    }

bool adapt_velocity(Mesh* mesh,  Mesh* rmesh, MeshFunction* sln, MeshFunction* rsln, Space* space)
{
  int i, j;
  
  int ne = mesh->get_max_element_id() + 1;
  double *elem_error = new double[ne];
  memset(elem_error, 0, sizeof(double) * ne);
  
  double error;

  int num = mesh->get_num_active_elements();
  ElemList* elist = new ElemList[num]; 

  Quad2D* quad = &g_quad_2d_std;
  sln->set_quad_2d(quad);
  rsln->set_quad_2d(quad);
    
  Mesh* meshes[2] = { mesh, rmesh };
  Transformable* tr[2] = { sln, rsln };
  Traverse trav;
  trav.begin(2, meshes, tr);

  Element** ee;
  while ((ee = trav.get_next_state(NULL, NULL)) != NULL)
  {
    RefMap* ru = sln->get_refmap();
    RefMap* rr = rsln->get_refmap();
    elem_error[ee[0]->id] += int_l2_error(sln, rsln, ru, rr);
    
  }
  trav.finish();

  Element* e;     
  double total_err = 0.0;
  int n = 0;
  for_all_active_elements(e, mesh)
  {
    elist[n].id = e->id;
    elist[n].error = elem_error[e->id];
    total_err += elem_error[e->id];
    n++;
  }

void Filter::init()
{
  // construct the union mesh, if necessary
  Mesh* meshes[10];
	for(int i = 0; i < this->num; i++)
		meshes[i] = this->sln[i]->get_mesh();
  mesh = meshes[0];
  unimesh = false;

  for (int i = 1; i < num; i++)
	{
    if (meshes[i] == NULL) {
      warn("You may be initializing a Filter with Solution that is missing a Mesh.");
      error("this->meshes[%d] is NULL in Filter::init().", i);
    }
    if (meshes[i]->get_seq() != mesh->get_seq())
    {
			unimesh = true;
			break;
		}
  }

  if (unimesh)
  {
    Traverse trav;
    trav.begin(num, meshes);
    mesh = new Mesh;
    unidata = trav.construct_union_mesh(mesh);
    trav.finish();
  }

  // misc init
  num_components = 1;
  order = 0;

  for(int i = 0; i < 10; i++)
      tables[i] = new LightArray<LightArray<Node*>*>;

  memset(sln_sub, 0, sizeof(sln_sub));
  set_quad_2d(&g_quad_2d_std);
}

double HcurlOrthoHP::calc_error_n(int n, ...)
{
  int i, j, k;

  if (n != num) error("Wrong number of solutions.");

  // obtain solutions and bilinear forms
  va_list ap;
  va_start(ap, n);
  for (i = 0; i < n; i++) {
    sln[i] = va_arg(ap, Solution*);
    sln[i]->set_quad_2d(&g_quad_2d_std);
  }
  for (i = 0; i < n; i++) {
    rsln[i] = va_arg(ap, Solution*);
    rsln[i]->set_quad_2d(&g_quad_2d_std);
  }
  va_end(ap);

  // prepare multi-mesh traversal and error arrays
  AUTOLA_OR(Mesh*, meshes, 2*num);
  AUTOLA_OR(Transformable*, tr, 2*num);
  Traverse trav;
  nact = 0;
  for (i = 0; i < num; i++)
  {
    meshes[i] = sln[i]->get_mesh();
    meshes[i+num] = rsln[i]->get_mesh();
    tr[i] = sln[i];
    tr[i+num] = rsln[i];

    nact += sln[i]->get_mesh()->get_num_active_elements();

    int max = meshes[i]->get_max_element_id();
    if (errors[i] != NULL) delete [] errors[i];
    errors[i] = new double[max];
    memset(errors[i], 0, sizeof(double) * max);
  }

  double total_norm = 0.0;
  AUTOLA_OR(double, norms, num);
  memset(norms, 0, num*sizeof(double));
  double total_error = 0.0;
  if (esort != NULL) delete [] esort;
  esort = new int2[nact];

  Element** ee;
  trav.begin(2*num, meshes, tr);
  while ((ee = trav.get_next_state(NULL, NULL)) != NULL)
  {
    for (i = 0; i < num; i++)
    {
      RefMap* rmi = sln[i]->get_refmap();
      RefMap* rrmi = rsln[i]->get_refmap();
      for (j = 0; j < num; j++)
      {
        RefMap* rmj = sln[j]->get_refmap();
        RefMap* rrmj = rsln[j]->get_refmap();
        double e, t;
        if (form[i][j] != NULL)
        {
          #ifndef COMPLEX
          e = fabs(eval_error(form[i][j], ord[i][j], sln[i], sln[j], rsln[i], rsln[j], rmi, rmj, rrmi, rrmj));
          t = fabs(eval_norm(form[i][j], ord[i][j], rsln[i], rsln[j], rrmi, rrmj));
          #else
          e = std::abs(eval_error(form[i][j], ord[i][j], sln[i], sln[j], rsln[i], rsln[j], rmi, rmj, rrmi, rrmj));
          t = std::abs(eval_norm(form[i][j], ord[i][j], rsln[i], rsln[j], rrmi, rrmj));
          #endif

          norms[i] += t;
          total_norm  += t;
          total_error += e;
          errors[i][ee[i]->id] += e;
        }
      }
    }
  }
  trav.finish();

  Element* e;
  k = 0;
  for (i = 0; i < num; i++)
    for_all_active_elements(e, meshes[i]) {
      esort[k][0] = e->id;
      esort[k++][1] = i;
      errors[i][e->id] /= norms[i];
    }

//.........这里部分代码省略.........
                stage.fns[i]->set_transform(neighbor_searches.get(ns_index)->original_central_el_transform);

              rm->set_transform(neighbor_searches.get(ns_index)->original_central_el_transform);

              
              /* END COPY FROM DISCRETE_PROBLEM.CPP */
            }
            
            /* BEGIN COPY FROM DISCRETE_PROBLEM.CPP */
            
            if (multimesh)
              // Delete the multimesh tree;
              delete root;
            
            // Delete the neighbor_searches array.
            for(unsigned int j = 0; j < neighbor_searches.get_size(); j++) 
              if(neighbor_searches.present(j))
                delete neighbor_searches.get(j);
              
            /* END COPY FROM DISCRETE_PROBLEM.CPP */
            
          }
        }
      }

      if (calc_norm)
      {
        double nrm = eval_solution_norm(error_form[i][i], rm, sln[i]);
        norms[i] += nrm;
        total_norm += nrm;
      }

      errors_components[i] += err;
      total_error += err;
      errors[i][ee[i]->id] += err;

      ee[i]->visited = true;
    }
  }
  trav.finish();

  // Store the calculation for each solution component separately.
  if(component_errors != NULL)
  {
    component_errors->clear();
    for (int i = 0; i < num; i++)
    {
      if((error_flags & HERMES_TOTAL_ERROR_MASK) == HERMES_TOTAL_ERROR_ABS)
        component_errors->push_back(sqrt(errors_components[i]));
      else if ((error_flags & HERMES_TOTAL_ERROR_MASK) == HERMES_TOTAL_ERROR_REL)
        component_errors->push_back(sqrt(errors_components[i]/norms[i]));
      else
      {
        error("Unknown total error type (0x%x).", error_flags & HERMES_TOTAL_ERROR_MASK);
        return -1.0;
      }
    }
  }

  tmr.tick();
  error_time = tmr.accumulated();

  // Make the error relative if needed.
  if ((error_flags & HERMES_ELEMENT_ERROR_MASK) == HERMES_ELEMENT_ERROR_REL)
  {
    for (int i = 0; i < num; i++)
    {
      Element* e;
      for_all_active_elements(e, stage.meshes[i])
        errors[i][e->id] /= norms[i];
    }
  }

  this->errors_squared_sum = total_error;

  // Element error mask is used here, because this variable is used in the adapt()
  // function, where the processed error (sum of errors of processed element errors)
  // is matched to this variable.
  if ((error_flags & HERMES_TOTAL_ERROR_MASK) == HERMES_ELEMENT_ERROR_REL)
    errors_squared_sum /= total_norm;

  // Prepare an ordered list of elements according to an error.
  fill_regular_queue(&(stage.meshes.front()));
  have_errors = true;

  if (calc_norm)
    delete [] norms;
  delete [] errors_components;

  // Return error value.
  if ((error_flags & HERMES_TOTAL_ERROR_MASK) == HERMES_TOTAL_ERROR_ABS)
    return sqrt(total_error);
  else if ((error_flags & HERMES_TOTAL_ERROR_MASK) == HERMES_TOTAL_ERROR_REL)
    return sqrt(total_error / total_norm);
  else
  {
    error("Unknown total error type (0x%x).", error_flags & HERMES_TOTAL_ERROR_MASK);
    return -1.0;
  }
}

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

      //// assemble volume linear forms ////////////////////////////////////////
      if (rhs != NULL)
      {
        for (unsigned int ww = 0; ww < s->vfvol.size(); ww++)
        {
          WeakForm::VectorFormVol* vfv = s->vfvol[ww];
          if (isempty[vfv->i]) continue;
          if (vfv->area != H2D_ANY && !wf->is_in_area(marker, vfv->area)) continue;
          m = vfv->i;  fv = spss[m];  am = &al[m];

          for (int i = 0; i < am->cnt; i++)
          {
            if (am->dof[i] < 0) continue;
            fv->set_active_shape(am->idx[i]);
            rhs->add(am->dof[i], eval_form(vfv, u_ext, fv, refmap + m) * am->coef[i]);
          }
        }
      }

      // assemble surface integrals now: loop through boundary edges of the element
      for (unsigned int edge = 0; edge < e0->nvert; edge++)
      {
        if (!bnd[edge]) continue;
        marker = ep[edge].marker;

        // obtain the list of shape functions which are nonzero on this edge
        for (unsigned int i = 0; i < s->idx.size(); i++) {
          if (e[i] == NULL) continue;
          int j = s->idx[i];
          if ((nat[j] = (spaces[j]->bc_type_callback(marker) == BC_NATURAL)))
            spaces[j]->get_edge_assembly_list(e[i], edge, al + j);
        }

        // assemble surface matrix forms ///////////////////////////////////
        if (jac != NULL)
        {
          for (unsigned int ww = 0; ww < s->mfsurf.size(); ww++)
          {
            WeakForm::MatrixFormSurf* mfs = s->mfsurf[ww];
            if (isempty[mfs->i] || isempty[mfs->j]) continue;
            if (mfs->area != H2D_ANY && !wf->is_in_area(marker, mfs->area)) continue;
            m = mfs->i;  fv = spss[m];  am = &al[m];
            n = mfs->j;  fu = pss[n];   an = &al[n];

            if (!nat[m] || !nat[n]) continue;
            ep[edge].base = trav.get_base();
            ep[edge].space_v = spaces[m];
            ep[edge].space_u = spaces[n];

            scalar bi, **mat = get_matrix_buffer(std::max(am->cnt, an->cnt));
            for (int i = 0; i < am->cnt; i++)
            {
              if ((k = am->dof[i]) < 0) continue;
              fv->set_active_shape(am->idx[i]);
              for (int j = 0; j < an->cnt; j++)
              {
                fu->set_active_shape(an->idx[j]);
                bi = eval_form(mfs, u_ext, fu, fv, refmap+n, refmap+m, ep+edge) * an->coef[j] * am->coef[i];
                if (an->dof[j] >= 0) mat[i][j] = bi;
              }
            }
            jac->add(am->cnt, an->cnt, mat, am->dof, an->dof);
          }
        }
        // assemble surface linear forms /////////////////////////////////////
        if (rhs != NULL)
        {
          for (unsigned ww = 0; ww < s->vfsurf.size(); ww++)
          {
            WeakForm::VectorFormSurf* vfs = s->vfsurf[ww];
            if (isempty[vfs->i]) continue;
            if (vfs->area != H2D_ANY && !wf->is_in_area(marker, vfs->area)) continue;
            m = vfs->i;  fv = spss[m];  am = &al[m];

            if (!nat[m]) continue;
            ep[edge].base = trav.get_base();
            ep[edge].space_v = spaces[m];

            for (int i = 0; i < am->cnt; i++)
            {
              if (am->dof[i] < 0) continue;
              fv->set_active_shape(am->idx[i]);
              rhs->add(am->dof[i], eval_form(vfs, u_ext, fv, refmap+m, ep+edge) * am->coef[i]);
            }
          }
        }
      }
      delete_cache();
    }
    trav.finish();
  }

  for (int i = 0; i < wf->neq; i++) delete spss[i];
  delete [] buffer;
  buffer = NULL;
  mat_size = 0;

}

void FeProblem::create(SparseMatrix* mat)
{
  assert(mat != NULL);

  if (is_up_to_date())
  {
    verbose("Reusing matrix sparse structure.");
    mat->zero();
    return;
  }

  // spaces have changed: create the matrix from scratch
  mat->free();

  int ndof = get_num_dofs();
  mat->prealloc(ndof);

  AUTOLA_CL(AsmList, al, wf->neq);
  AUTOLA_OR(Mesh*, meshes, wf->neq);
  bool **blocks = wf->get_blocks();

  // init multi-mesh traversal
  for (int i = 0; i < wf->neq; i++)
    meshes[i] = spaces[i]->get_mesh();

  Traverse trav;
  trav.begin(wf->neq, meshes);

  // loop through all elements
  Element **e;
  while ((e = trav.get_next_state(NULL, NULL)) != NULL)
  {
    // obtain assembly lists for the element at all spaces
    for (int i = 0; i < wf->neq; i++)
      // TODO: do not get the assembly list again if the element was not changed
      if (e[i] != NULL)
        spaces[i]->get_element_assembly_list(e[i], al + i);

    // go through all equation-blocks of the local stiffness matrix
    for (int m = 0; m < wf->neq; m++)
      for (int n = 0; n < wf->neq; n++)
        if (blocks[m][n] && e[m] != NULL && e[n] != NULL)
        {
          AsmList *am = al + m;
          AsmList *an = al + n;

          // pretend assembling of the element stiffness matrix
          // register nonzero elements
          for (int i = 0; i < am->cnt; i++)
            if (am->dof[i] >= 0)
              for (int j = 0; j < an->cnt; j++)
                if (an->dof[j] >= 0)
                  mat->pre_add_ij(am->dof[i], an->dof[j]);
        }
  }

  trav.finish();
  delete [] blocks;

  mat->alloc();

  // save space seq numbers and weakform seq number, so we can detect their changes
  for (int i = 0; i < wf->neq; i++)
    sp_seq[i] = spaces[i]->get_seq();
  wf_seq = wf->get_seq();

  struct_changed = true;
  have_matrix = true;
}

//.........这里部分代码省略.........
                  fu->set_active_shape(an->idx[j]);
                  mat[i][j] += bf->unsym(fu, fv, refmap+n, refmap+m) * coef;
                }
              }
            }
          }

          // insert the local stiffness matrix into the global one
          insert_matrix(mat, am->dof, an->dof, am->cnt, an->cnt);

          // insert also the off-diagonal (anti-)symmetric block, if required
          if (tra)
          {
            if (biform[n][m].unsym == BF_ANTISYM) chsgn(mat, am->cnt, an->cnt);
            transpose(mat, am->cnt, an->cnt);
            insert_matrix(mat, an->dof, am->dof, an->cnt, am->cnt);

            // we also need to take care of the RHS...
            for (j = 0; j < am->cnt; j++)
              if (am->dof[j] < 0)
                for (i = 0; i < an->cnt; i++)
                  if (an->dof[i] >= 0)
                    Dir[an->dof[i]] -= mat[i][j];
          }
        }
      }

      // assemble rhs (linear form)
      if (!liform[m].lf) continue;
      for (i = 0; i < am->cnt; i++)
      {
        if (am->dof[i] < 0) continue;
        fv->set_active_shape(am->idx[i]);
        RHS[am->dof[i]] += liform[m].lf(fv, refmap+m) * am->coef[i];
      }
    }

    // assemble surface integrals now: loop through boundary edges of the element
    if (rhsonly) continue; // fixme
    for (int edge = 0; edge < e[0]->nvert; edge++)
    {
      if (!bnd[edge]) continue;

      // obtain the list of shape functions which are nonzero on this edge
      for (i = 0; i < neq; i++)
        if ((nat[i] = (spaces[i]->bc_type_callback(ep[edge].marker) == BC_NATURAL)))
          spaces[i]->get_edge_assembly_list(e[i], edge, al + i);

      // loop through the equation-blocks
      for (m = 0, am = al; m < neq; m++, am++)
      {
        if (!nat[m]) continue;
        fv = spss[m];
        ep[edge].base = trav.get_base();
        ep[edge].space_v = spaces[m];
        for (n = 0, an = al; n < neq; n++, an++)
        {
          if (!nat[n]) continue;
          BiForm* bf = biform[m] + n;
          if (!bf->surf) continue;
          fu = pss[n];
          ep[edge].space_u = spaces[n];

          // assemble the surface part of the bilinear form
          scalar bi, **mat = get_matrix_buffer(std::max(am->cnt, an->cnt));
          for (i = 0; i < am->cnt; i++)
          {
            if ((k = am->dof[i]) < 0) continue;
            fv->set_active_shape(am->idx[i]);
            for (j = 0; j < an->cnt; j++)
            {
              fu->set_active_shape(an->idx[j]);
              bi = bf->surf(fu, fv, refmap+n, refmap+m, ep+edge) * an->coef[j] * am->coef[i];
              if (an->dof[j] >= 0) mat[i][j] = bi; else Dir[k] -= bi;
            }
          }
          insert_matrix(mat, am->dof, an->dof, am->cnt, an->cnt);
        }

        // assemble the surface part of the linear form
        if (!liform[m].surf) continue;
        for (i = 0; i < am->cnt; i++)
        {
          if (am->dof[i] < 0) continue;
          fv->set_active_shape(am->idx[i]);
          RHS[am->dof[i]] += liform[m].surf(fv, refmap+m, ep+edge) * am->coef[i];
        }
      }
    }
  }

  trav.finish();
  for (i = 0; i < ndofs; i++)
    RHS[i] += Dir[i];

  if (!quiet) verbose("  (time: %g sec)", end_time());
  for (i = 0; i < neq; i++) delete spss[i];
  delete [] buffer;
  delete [] refmap;
}

void Vectorizer::process_solution(MeshFunction* xsln, int xitem, MeshFunction* ysln, int yitem, double eps)
{
  // sanity check
  if (xsln == NULL || ysln == NULL) error("One of the solutions is NULL in Vectorizer:process_solution().");


  lock_data();
  TimePeriod cpu_time;

  // initialization
  this->xsln = xsln;
  this->ysln = ysln;
  this->xitem = xitem;
  this->yitem = yitem;
  this->eps = eps;
  nv = nt = ne = nd = 0;
  del_slot = -1;

  Mesh* meshes[2] = { xsln->get_mesh(), ysln->get_mesh() };
  if (meshes[0] == NULL || meshes[1] == NULL) {
    error("One of the meshes is NULL in Vectorizer:process_solution().");
  }

  Transformable* fns[2] = { xsln, ysln };
  Traverse trav;

  // estimate the required number of vertices and triangles
  // (based on the assumption that the linear mesh will be
  // about four-times finer than the original mesh).
  int nn = meshes[0]->get_num_elements() + meshes[1]->get_num_elements();
  int ev = std::max(32 * nn, 10000);
  int et = std::max(64 * nn, 20000);
  int ee = std::max(24 * nn, 7500);
  int ed = ee;

  lin_init_array(verts, double4, cv, ev);
  lin_init_array(tris, int3, ct, et);
  lin_init_array(edges, int3, ce, ee);
  lin_init_array(dashes, int2, cd, ed);

  info = (int4*) malloc(sizeof(int4) * cv);

  // initialize the hash table
  int size = 0x1000;
  while (size*2 < cv) size *= 2;
  hash_table = (int*) malloc(sizeof(int) * size);
  memset(hash_table, 0xff, sizeof(int) * size);
  mask = size-1;


  // select the linearization quadrature
  Quad2D *old_quad_x, *old_quad_y;
  old_quad_x = xsln->get_quad_2d();
  old_quad_y = ysln->get_quad_2d();

  xsln->set_quad_2d((Quad2D*) &quad_lin);
  ysln->set_quad_2d((Quad2D*) &quad_lin);

  if (!xitem) error("Parameter 'xitem' cannot be zero.");
  if (!yitem) error("Parameter 'yitem' cannot be zero.");
  get_gv_a_b(xitem, xia, xib);
  get_gv_a_b(yitem, yia, yib);
  if (xib >= 6) error("Invalid value of paremeter 'xitem'.");
  if (yib >= 6) error("Invalid value of paremeter 'yitem'.");

  max = 1e-10;
  trav.begin(2, meshes, fns);
  Element** e;
  while ((e = trav.get_next_state(NULL, NULL)) != NULL)
  {
    xsln->set_quad_order(0, xitem);
    ysln->set_quad_order(0, yitem);
    scalar* xval = xsln->get_values(xia, xib);
    scalar* yval = ysln->get_values(yia, yib);

    for (unsigned int i = 0; i < e[0]->nvert; i++)
    {
      double fx = getvalx(i);
      double fy = getvaly(i);
      if (fabs(sqrt(fx*fx + fy*fy)) > max) max = fabs(sqrt(fx*fx + fy*fy));
    }
  }
  trav.finish();

  trav.begin(2, meshes, fns);
  // process all elements of the mesh
  while ((e = trav.get_next_state(NULL, NULL)) != NULL)
  {
    xsln->set_quad_order(0, xitem);
    ysln->set_quad_order(0, yitem);
    scalar* xval = xsln->get_values(xia, xib);
    scalar* yval = ysln->get_values(yia, yib);

    double* x = xsln->get_refmap()->get_phys_x(0);
    double* y = ysln->get_refmap()->get_phys_y(0);

    int iv[4];
    for (unsigned int i = 0; i < e[0]->nvert; i++)
    {
      double fx = getvalx(i);
//.........这里部分代码省略.........