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

double KellyTypeAdapt::calc_err_internal(Hermes::vector<Solution *> slns,
                                         Hermes::vector<double>* component_errors,
                                         unsigned int error_flags)
{    
  int n = slns.size();
  error_if (n != this->num,
            "Wrong number of solutions.");

  TimePeriod tmr;

  for (int i = 0; i < n; i++)
  {
    this->sln[i] = slns[i];
    sln[i]->set_quad_2d(&g_quad_2d_std);
  }

  have_coarse_solutions = true;
  
  WeakForm::Stage stage;

  num_act_elems = 0;
  for (int i = 0; i < num; i++)
  {
    stage.meshes.push_back(sln[i]->get_mesh());
    stage.fns.push_back(sln[i]);

    num_act_elems += stage.meshes[i]->get_num_active_elements();
    int max = stage.meshes[i]->get_max_element_id();

    if (errors[i] != NULL) delete [] errors[i];
    errors[i] = new double[max];
    memset(errors[i], 0.0, sizeof(double) * max);
  }
/*
  for (unsigned int i = 0; i < error_estimators_vol.size(); i++)
    trset.insert(error_estimators_vol[i].ext.begin(), error_estimators_vol[i].ext.end());
  for (unsigned int i = 0; i < error_estimators_surf.size(); i++)
    trset.insert(error_estimators_surf[i].ext.begin(), error_estimators_surf[i].ext.end());
*/

  double total_norm = 0.0;

  bool calc_norm = false;
  if ((error_flags & HERMES_ELEMENT_ERROR_MASK) == HERMES_ELEMENT_ERROR_REL ||
      (error_flags & HERMES_TOTAL_ERROR_MASK) == HERMES_TOTAL_ERROR_REL) calc_norm = true;

  double *norms = NULL;
  if (calc_norm)
  {
    norms = new double[num];
    memset(norms, 0.0, num * sizeof(double));
  }

  double *errors_components = new double[num];
  memset(errors_components, 0.0, num * sizeof(double));
  this->errors_squared_sum = 0.0;
  double total_error = 0.0;

  bool bnd[4];          // FIXME: magic number - maximal possible number of element surfaces
  SurfPos surf_pos[4];
  Element **ee;
  Traverse trav;

  // Reset the e->visited status of each element of each mesh (most likely it will be set to true from
  // the latest assembling procedure).
  if (ignore_visited_segments)
  {
    for (int i = 0; i < num; i++)
    {
      Element* e;
      for_all_active_elements(e, stage.meshes[i])
        e->visited = false;
    }
  }
  
  //WARNING: AD HOC debugging parameter.
  bool multimesh = false;

  // Begin the multimesh traversal.
  trav.begin(num, &(stage.meshes.front()), &(stage.fns.front()));
  while ((ee = trav.get_next_state(bnd, surf_pos)) != NULL)
  {   
    // Go through all solution components.
    for (int i = 0; i < num; i++)
    {
      if (ee[i] == NULL)
        continue;
      
      // Set maximum integration order for use in integrals, see limit_order()
      update_limit_table(ee[i]->get_mode());

      RefMap *rm = sln[i]->get_refmap();

      double err = 0.0;

      // Go through all volumetric error estimators.
      for (unsigned int iest = 0; iest < error_estimators_vol.size(); iest++)
      {
        // Skip current error estimator if it is assigned to a different component or geometric area
        // different from that of the current active element.
//.........这里部分代码省略.........

void FeProblem::assemble(_Vector *rhs, _Matrix *jac, Tuple<Solution*> u_ext)
{
  // Sanity checks.
  if (!have_spaces) error("You have to call set_spaces() before calling assemble().");
  for (int i=0; i<this->wf->neq; i++) {
    if (this->spaces[i] == NULL) error("A space is NULL in assemble().");
  }

  int k, m, n, marker;
  AUTOLA_CL(AsmList, al, wf->neq);
  AsmList *am, *an;
  bool bnd[4];
  AUTOLA_OR(bool, nat, wf->neq);
  AUTOLA_OR(bool, isempty, wf->neq);
  EdgePos ep[4];
  reset_warn_order();

  // create slave pss's for test functions, init quadrature points
  AUTOLA_OR(PrecalcShapeset*, spss, wf->neq);
  PrecalcShapeset *fu, *fv;
  AUTOLA_CL(RefMap, refmap, wf->neq);
  for (int i = 0; i < wf->neq; i++)
  {
    spss[i] = new PrecalcShapeset(pss[i]);
    pss [i]->set_quad_2d(&g_quad_2d_std);
    spss[i]->set_quad_2d(&g_quad_2d_std);
    refmap[i].set_quad_2d(&g_quad_2d_std);
  }

  // initialize buffer
  buffer = NULL;
  mat_size = 0;
  get_matrix_buffer(9);

  // obtain a list of assembling stages
  std::vector<WeakForm::Stage> stages;
  wf->get_stages(spaces, NULL, stages, jac == NULL);

  // Loop through all assembling stages -- the purpose of this is increased performance
  // in multi-mesh calculations, where, e.g., only the right hand side uses two meshes.
  // In such a case, the matrix forms are assembled over one mesh, and only the rhs
  // traverses through the union mesh. On the other hand, if you don't use multi-mesh
  // at all, there will always be only one stage in which all forms are assembled as usual.
  Traverse trav;
  for (unsigned ss = 0; ss < stages.size(); ss++)
  {
    WeakForm::Stage* s = &stages[ss];
    for (unsigned i = 0; i < s->idx.size(); i++)
      s->fns[i] = pss[s->idx[i]];
    for (unsigned i = 0; i < s->ext.size(); i++)
      s->ext[i]->set_quad_2d(&g_quad_2d_std);
    trav.begin(s->meshes.size(), &(s->meshes.front()), &(s->fns.front()));

    // assemble one stage
    Element** e;
    while ((e = trav.get_next_state(bnd, ep)) != NULL)
    {
      // find a non-NULL e[i]
      Element* e0;
      for (unsigned int i = 0; i < s->idx.size(); i++)
        if ((e0 = e[i]) != NULL) break;
      if (e0 == NULL) continue;

      // set maximum integration order for use in integrals, see limit_order()
      update_limit_table(e0->get_mode());

      // obtain assembly lists for the element at all spaces, set appropriate mode for each pss
      memset(isempty, 0, sizeof(bool) * wf->neq);
      for (unsigned int i = 0; i < s->idx.size(); i++)
      {
        int j = s->idx[i];
        if (e[i] == NULL) { isempty[j] = true; continue; }
        spaces[j]->get_element_assembly_list(e[i], al+j);

        spss[j]->set_active_element(e[i]);
        spss[j]->set_master_transform();
        refmap[j].set_active_element(e[i]);
        refmap[j].force_transform(pss[j]->get_transform(), pss[j]->get_ctm());

        u_ext[j]->set_active_element(e[i]);
        u_ext[j]->force_transform(pss[j]->get_transform(), pss[j]->get_ctm());
      }
      marker = e0->marker;

      init_cache();
      //// assemble volume matrix forms //////////////////////////////////////
      if (jac != NULL)
      {
        for (unsigned ww = 0; ww < s->mfvol.size(); ww++)
        {
          WeakForm::MatrixFormVol* mfv = s->mfvol[ww];
          if (isempty[mfv->i] || isempty[mfv->j]) continue;
          if (mfv->area != H2D_ANY && !wf->is_in_area(marker, mfv->area)) continue;
          m = mfv->i;  fv = spss[m];  am = &al[m];
          n = mfv->j;  fu = pss[n];   an = &al[n];
          bool tra = (m != n) && (mfv->sym != 0);
          bool sym = (m == n) && (mfv->sym == 1);

          // assemble the local stiffness matrix for the form mfv
          scalar bi, **mat = get_matrix_buffer(std::max(am->cnt, an->cnt));
//.........这里部分代码省略.........

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

void DiscreteProblem::assemble_matrix_and_rhs(bool rhsonly)
{
  int i, j, k, l, m, n;
  bool bnd[4], nat[neq];
  EdgePos ep[4];

  if (!ndofs) return;
  warned_order = false;

  if (!rhsonly)
  {
    alloc_matrix_values();
    if (!quiet) { verbose("Assembling stiffness matrix..."); begin_time(); }
  }
  else
  {
    memset(RHS, 0, sizeof(scalar) * ndofs);
    if (!quiet) { verbose("Assembling RHS..."); begin_time(); }
  }

  // create slave pss's for test functions, init quadrature points
  PrecalcShapeset* spss[neq];
  PrecalcShapeset *fu, *fv;
  for (i = 0; i < neq; i++)
  {
    spss[i] = new PrecalcShapeset(pss[i]);
    pss [i]->set_quad_2d(&g_quad_2d_std);
    spss[i]->set_quad_2d(&g_quad_2d_std);
  }

  // initialize buffer
  buffer = NULL;
  mat_size = 0;
  get_matrix_buffer(9);

  // initialize assembly lists, refmap
  AsmList al[neq], *am, *an;

  RefMap* refmap = new RefMap[neq];
  for (i = 0; i < neq; i++)
    refmap[i].set_quad_2d(&g_quad_2d_std);
  for (i = 0; i < num_extern; i++)
    extern_fns[i]->set_quad_2d(&g_quad_2d_std);

  // init multi-mesh traversal
  int nm = neq + num_extern;
  Mesh* meshes[nm];
  Transformable* fn[nm];
  for (i = 0; i < neq; i++)
    meshes[i] = spaces[i]->get_mesh();
  memcpy(fn, pss, neq * sizeof(Transformable*));
  for (i = 0; i < num_extern; i++) {
    meshes[neq+i] = extern_fns[i]->get_mesh();
    fn[neq+i] = extern_fns[i];
  }
  // todo: kdyz maji nektere slozky stejnou sit, at sdili i refmapy
  //  - ale to bysme potrebovali slave RefMap

  // loop through all elements
  Element** e;
  Traverse trav;
  trav.begin(nm, meshes, fn);
  while ((e = trav.get_next_state(bnd, ep)) != NULL)
  {
    // set maximum integration order for use in integrals, see limit_order()
    update_limit_table(e[0]->get_mode());

    // obtain assembly lists for the element at all spaces, set appropriate mode for each pss
    for (i = 0; i < neq; i++)
    {
      spaces[i]->get_element_assembly_list(e[i], al + i);
      // todo: neziskavat znova, pokud se element nezmenil
      if (is_equi)
        for (j = 0; j < al[i].cnt; j++)
          if (al[i].dof[j] >= 0)
            al[i].coef[j] /= equi[al[i].dof[j]];

      spss[i]->set_active_element(e[i]);
      spss[i]->set_master_transform();
      refmap[i].set_active_element(e[i]);
      refmap[i].force_transform(pss[i]->get_transform(), pss[i]->get_ctm());
    }

    // go through all equation-blocks of the element stiffness matrix, assemble volume integrals
    for (m = 0, am = al; m < neq; m++, am++)
    {
      fv = spss[m];
      if (!rhsonly)
      {
        for (n = 0, an = al; n < neq; n++, an++)
        {
          fu = pss[n];
          BiForm* bf = biform[m] + n;
          if (!bf->sym && !bf->unsym) continue;
          if (bf->unsym == BF_SYM || bf->unsym == BF_ANTISYM) continue;
          bool tra = (biform[n][m].unsym == BF_SYM || biform[n][m].unsym == BF_ANTISYM);

          // assemble the (m,n)-block of the stiffness matrix
          scalar sy, un, **mat = get_matrix_buffer(std::max(am->cnt, an->cnt));
          for (i = 0; i < am->cnt; i++)
//.........这里部分代码省略.........

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);
//.........这里部分代码省略.........