FROM SPIntegral IMPORT OnTwoTriangles;

PROCEDURE DiffTriangIntegral(
    f, g: T;
    h: Integrand; 
    order: CARDINAL;
    fZero, gZero: BOOLEAN;
  ): LONG = 
  (*
    Computes the integral of "h(ft, gt, p)" over the intersection of
    every triangle "ft" of "f.tri" and every triangle "gt" of "g.tri",
    using at least "(order+1)(order+2)/2" sample points in each
    triangle.  If "fZero" is true, assumes that "h" is zero on any
    triangle where "f" is zero; and analogously for "gZero" and "g". *)
  
  VAR integral: LONG := 0.0d0;
      fi, gi: CARDINAL;
      ft, gt: TriangleData;
  BEGIN
    IF (fZero OR gZero) AND Disjoint(f.supp, g.supp) THEN RETURN 0.0d0 END;
    WITH
      fd = f.data^, fN = NUMBER(fd),
      fT = NUMBER(f.tri.side^),
      gd = g.data^, gN = NUMBER(gd),
      gT = NUMBER(g.tri.side^)
    DO
      SortTriangles(fd); 
      SortTriangles(gd);
      FOR i:= 0 TO fT-1 DO 
        IF fi < fN AND fd[fi].face = i THEN ft := fd[fi]; INC(fi) ELSE ft := NIL END;
        IF fZero AND ft = NIL THEN
          (* "h" is zero in "ft", skip it *)
        ELSE
          WITH 
            e1 = f.tri.side[i],
            T10 = Dest(Lnext(e1)).c,
            T11 = Org(e1).c,
            T12 = Dest(e1).c
          DO   
            FOR j:= 0 TO gT-1 DO 
              IF gi < gN AND gd[gi].face = j THEN gt := gd[gi]; INC(gi) ELSE gt := NIL END;
              IF gZero AND gt = NIL THEN
                (* "h" is zero in "gt", skip it. *)
              ELSE
                WITH  
                  e2 = g.tri.side[j],
                  T20 = Dest(Lnext(e2)).c,
                  T21 = Org(e2).c,
                  T22 = Dest(e2).c
                DO 

                  PROCEDURE HT(READONLY p: Point): LONG = BEGIN RETURN h(ft, gt, p) END HT;
                  
                  BEGIN 
                    integral := integral + 
                      OnTwoTriangles(T10,T11,T12, T20,T21,T22, HT, order)
                  END
                  
                END
              END
            END
          END
        END
      END
    END;
    RETURN integral;
  END DiffTriangIntegral;