Function CalcDifference

comp/postproc.cpp:1078–1178

  template <class SCAL>
  void CalcDifference (const S_GridFunction<SCAL> & u1,
		       const S_GridFunction<SCAL> & u2,
		       shared_ptr<BilinearFormIntegrator> bli1,
		       shared_ptr<BilinearFormIntegrator> bli2,
		       FlatVector<double> & diff,
		       int domain, LocalHeap & lh)
  {
    shared_ptr<MeshAccess> ma = u1.GetFESpace()->GetMeshAccess();
    BaseStatusHandler::PushStatus ("Calc Difference");

    const FESpace & fes1 = *u1.GetFESpace();
    const FESpace & fes2 = *u2.GetFESpace();

    bool bound1 = bli1->BoundaryForm();
    bool bound2 = bli2->BoundaryForm();


    if(bound1!=bound2) 
      {
	cout << " ERROR: CalcDifference :: bli1->BoundaryForm != bl2.BoundaryForm there is something wrong?" << endl; 
	diff = 0; 
	return; 
      } 

    int ne      = bound1 ? ma->GetNSE() : ma->GetNE();
    int dim1    = fes1.GetDimension();
    int dim2    = fes2.GetDimension();
    int dimflux1 = bli1->DimFlux();
    int dimflux2 = bli2->DimFlux();

    if(dimflux1 != dimflux2) 
      { 
	cout << " ERROR: CalcDifference :: dimflux1 != dimflux2 !!!!! -> set diff = 0" << endl; 
	diff = 0; 
	return; 	
      } 

    bool applyd1 = 0;
    bool applyd2 = 0;

    Array<int> dnums1;
    Array<int> dnums2;

    // double sum = 0;
    for (int i = 0; i < ne; i++)
      {
        ElementId ei(bound1 ? BND : VOL, i);
	HeapReset hr (lh);
        BaseStatusHandler::SetThreadPercentage ( 100.0*i / ne );

	int eldom = ma->GetElIndex(ei);
        // bound1 ? ma->GetSElIndex(i) : ma->GetElIndex(i);
	
	if ((domain != -1) && (domain != eldom))
	  continue;

	const FiniteElement & fel1 = fes1.GetFE (ei, lh);
	const FiniteElement & fel2 = fes2.GetFE (ei, lh);