Testing environment on.

- Basic conjugate gradient solver for partial assembly without preconditioning.
- Conjugate gradient solver for PA with Prec.
- Diag Solver (can be used as preconditionner).
- Conjugate gradient working on element level for DG methods.

config/config.hpp.in

@@ -58,6 +58,9 @@
 // Use LAPACK routines for various dense linear algebra operations.
 // #define MFEM_USE_LAPACK
 
+// Use Eigen for math routines
+// #define MFEM_USE_EIGEN
+
 // Use thread-safe implementation. This comes at the cost of extra memory
 // allocation and de-allocation.
 // #define MFEM_THREAD_SAFE

config/config.mk.in

@@ -21,6 +21,7 @@ MFEM_USE_EXCEPTIONS  = @MFEM_USE_EXCEPTIONS@
 MFEM_USE_GZSTREAM    = @MFEM_USE_GZSTREAM@
 MFEM_USE_LIBUNWIND   = @MFEM_USE_LIBUNWIND@
 MFEM_USE_LAPACK      = @MFEM_USE_LAPACK@
+MFEM_USE_EIGEN       = @MFEM_USE_EIGEN@
 MFEM_THREAD_SAFE     = @MFEM_THREAD_SAFE@
 MFEM_USE_OPENMP      = @MFEM_USE_OPENMP@
 MFEM_USE_MEMALLOC    = @MFEM_USE_MEMALLOC@

config/defaults.mk

@@ -21,7 +21,7 @@ NOTMAC := $(subst Darwin,,$(shell uname -s))
 CXX = g++
 MPICXX = mpicxx
 
-OPTIM_FLAGS = -O3
+OPTIM_FLAGS = -O3 -pg
 DEBUG_FLAGS = -g -Wall
 
 # Destination location of make install
@@ -90,6 +90,7 @@ MFEM_USE_EXCEPTIONS  = NO
 MFEM_USE_GZSTREAM    = NO
 MFEM_USE_LIBUNWIND   = NO
 MFEM_USE_LAPACK      = NO
+MFEM_USE_EIGEN       = NO
 MFEM_THREAD_SAFE     = NO
 MFEM_USE_OPENMP      = NO
 MFEM_USE_MEMALLOC    = YES
@@ -145,6 +146,11 @@ endif
 LAPACK_OPT =
 LAPACK_LIB = $(if $(NOTMAC),-llapack -lblas,-framework Accelerate)
 
+# Eigen configuration
+EIGEN_DIR = @MFEM_DIR@/../eigen
+EIGEN_OPT = -I$(EIGEN_DIR) -std=c++11 -Wno-enum-compare
+EIGEN_LIB =
+
 # OpenMP configuration
 OPENMP_OPT = -fopenmp
 OPENMP_LIB =

examples/ex1.cpp

@@ -48,7 +48,9 @@ int main(int argc, char *argv[])
    const char *mesh_file = "../data/star.mesh";
    int order = 1;
    bool static_cond = false;
-   bool visualization = 1;
+   bool visualization = true;
+   bool use_partial_assembly = false;
+   bool use_smoother = true;
 
    OptionsParser args(argc, argv);
    args.AddOption(&mesh_file, "-m", "--mesh",
@@ -56,8 +58,12 @@ int main(int argc, char *argv[])
    args.AddOption(&order, "-o", "--order",
                   "Finite element order (polynomial degree) or -1 for"
                   " isoparametric space.");
+   args.AddOption(&use_partial_assembly, "-pa", "--partial-assembly",
+                  "-no-pa", "--no-partial-assembly", "Enable partial assembly.");
    args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
                   "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&use_smoother, "-pc", "--peconditioner", "-no-pc",
+                  "--no-preconditioner", "Use a Gauss-Seidel preconditioner.");
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
@@ -138,33 +144,56 @@ int main(int argc, char *argv[])
    // 8. Set up the bilinear form a(.,.) on the finite element space
    //    corresponding to the Laplacian operator -Delta, by adding the Diffusion
    //    domain integrator.
+   Vector B, X;
    BilinearForm *a = new BilinearForm(fespace);
    a->AddDomainIntegrator(new DiffusionIntegrator(one));
 
+   BilinearFormOperator A_pa(new PAIntegratorMap);
+   SparseMatrix A_sp;
+   if (!use_partial_assembly)
+   {
+      a->AssembleForm(A_sp);
+   }
+   else
+   {
+      // Can add a custom FESpaceIntegrator in this way:
+      // a->AddIntegrator(new PADiffusionIntegrator(new DiffusionIntegrator(one)));
+      a->AssembleForm(A_pa);
+   }
+
+   Operator *A;
+   a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
+
    // 9. Assemble the bilinear form and the corresponding linear system,
    //    applying any necessary transformations such as: eliminating boundary
    //    conditions, applying conforming constraints for non-conforming AMR,
    //    static condensation, etc.
-   if (static_cond) { a->EnableStaticCondensation(); }
-   a->Assemble();
 
-   SparseMatrix A;
-   Vector B, X;
-   a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
-
-   cout << "Size of linear system: " << A.Height() << endl;
+   cout << "Size of linear system: " << A->Height() << endl;
 
 #ifndef MFEM_USE_SUITESPARSE
    // 10. Define a simple symmetric Gauss-Seidel preconditioner and use it to
    //     solve the system A X = B with PCG.
-   GSSmoother M(A);
-   PCG(A, M, B, X, 1, 200, 1e-12, 0.0);
+
+   if (use_smoother && !use_partial_assembly)
+   {
+      GSSmoother M(A_sp);
+      PCG(*A, M, B, X, 1, 200, 1e-12, 0.0);
+   }
+   else
+   {
+      CG(*A, B, X, 1, 200, 1e-12, 0.0);
+   }
+
 #else
    // 10. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
-   UMFPackSolver umf_solver;
-   umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
-   umf_solver.SetOperator(A);
-   umf_solver.Mult(B, X);
+   if (!use_partial_assembly)
+   {
+      UMFPackSolver umf_solver;
+      umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
+      umf_solver.SetOperator(A_sp);
+      umf_solver.Mult(B, X);
+   }
 #endif
 
    // 11. Recover the solution as a finite element grid function.
@@ -190,6 +219,7 @@ int main(int argc, char *argv[])
    }
 
    // 14. Free the used memory.
+   delete A;
    delete a;
    delete b;
    delete fespace;

examples/ex16.cpp

@@ -35,6 +35,10 @@
 using namespace std;
 using namespace mfem;
 
+bool partial_assembly_mass;
+bool partial_assembly_diff;
+bool preconditioner;
+
 /** After spatial discretization, the conduction model can be written as:
  *
  *     du/dt = M^{-1}(-Ku)
@@ -54,6 +58,8 @@ protected:
    BilinearForm *M;
    BilinearForm *K;
 
+   Operator *Koper, *Moper, *Toper;
+   BilinearFormOperator Mpaop, Kpaop;
    SparseMatrix Mmat, Kmat;
    SparseMatrix *T; // T = M + dt K
    double current_dt;
@@ -83,11 +89,39 @@ public:
    virtual ~ConductionOperator();
 };
 
+class TimeDerivativeOperator : public Operator
+{
+   Operator *Moper;
+   Operator *Koper;
+   mutable Vector Kdu;
+   const double dt;
+
+public:
+   TimeDerivativeOperator(Operator *_Moper, const double _dt, Operator *_Koper)
+      : Operator(_Moper->Height(), _Moper->Width()),
+        Moper(_Moper),
+        Koper(_Koper),
+        Kdu(Height()),
+        dt(_dt) { }
+
+   virtual void Mult(const Vector &x, Vector &y) const
+   {
+      Moper->Mult(x, y);
+      Koper->Mult(x, Kdu);
+
+      Kdu *= dt;
+      y += Kdu;
+   }
+};
+
 double InitialTemperature(const Vector &x);
 
 int main(int argc, char *argv[])
 {
    // 1. Parse command-line options.
+   partial_assembly_mass = false;
+   partial_assembly_diff = false;
+   preconditioner = true;
    const char *mesh_file = "../data/star.mesh";
    int ref_levels = 2;
    int order = 2;
@@ -121,6 +155,12 @@ int main(int argc, char *argv[])
                   "Alpha coefficient.");
    args.AddOption(&kappa, "-k", "--kappa",
                   "Kappa coefficient offset.");
+   args.AddOption(&partial_assembly_mass, "-pam", "--partial-assembly-mass",
+                  "-no-pam", "--no-partial-assembly-mass", "Enable partial assembly for the mass.");
+   args.AddOption(&partial_assembly_diff, "-pad", "--partial-assembly-diff",
+                  "-no-pad", "--no-partial-assembly-diff", "Enable partial assembly for the diffusion.");
+   args.AddOption(&preconditioner, "-pc", "--peconditioner", "-no-pc",
+                  "--no-preconditioner", "Use a Gauss-Seidel preconditioner.");
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
@@ -291,22 +331,35 @@ int main(int argc, char *argv[])
 ConductionOperator::ConductionOperator(FiniteElementSpace &f, double al,
                                        double kap, const Vector &u)
    : TimeDependentOperator(f.GetTrueVSize(), 0.0), fespace(f), M(NULL), K(NULL),
+     Toper(NULL), Mpaop(new PAIntegratorMap), Kpaop(new PAIntegratorMap),
      T(NULL), current_dt(0.0), z(height)
 {
    const double rel_tol = 1e-8;
 
    M = new BilinearForm(&fespace);
-   M->AddDomainIntegrator(new MassIntegrator());
-   M->Assemble();
-   M->FormSystemMatrix(ess_tdof_list, Mmat);
+   M->AddDomainIntegrator(new MassIntegrator);
+   if (!partial_assembly_mass)
+   {
+      M->AssembleForm(Mmat);
+      M->FormSystemOperator(ess_tdof_list, Moper);
+      M_solver.SetOperator(static_cast<SparseMatrix&>(*Moper));
+   }
+   else
+   {
+      M->AssembleForm(Mpaop);
+      M->FormSystemOperator(ess_tdof_list, Moper);
+      M_solver.SetOperator(*Moper);
+   }
 
    M_solver.iterative_mode = false;
    M_solver.SetRelTol(rel_tol);
    M_solver.SetAbsTol(0.0);
    M_solver.SetMaxIter(30);
    M_solver.SetPrintLevel(0);
-   M_solver.SetPreconditioner(M_prec);
-   M_solver.SetOperator(Mmat);
+   if (preconditioner && !partial_assembly_mass)
+   {
+      M_solver.SetPreconditioner(M_prec);
+   }
 
    alpha = al;
    kappa = kap;
@@ -316,7 +369,12 @@ ConductionOperator::ConductionOperator(FiniteElementSpace &f, double al,
    T_solver.SetAbsTol(0.0);
    T_solver.SetMaxIter(100);
    T_solver.SetPrintLevel(0);
-   T_solver.SetPreconditioner(T_prec);
+   if (preconditioner &&
+       !partial_assembly_diff &&
+       !partial_assembly_mass)
+   {
+      T_solver.SetPreconditioner(T_prec);
+   }
 
    SetParameters(u);
 }
@@ -326,7 +384,7 @@ void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
    // Compute:
    //    du_dt = M^{-1}*-K(u)
    // for du_dt
-   Kmat.Mult(u, z);
+   Koper->Mult(u, z);
    z.Neg(); // z = -z
    M_solver.Mult(z, du_dt);
 }
@@ -337,14 +395,22 @@ void ConductionOperator::ImplicitSolve(const double dt,
    // Solve the equation:
    //    du_dt = M^{-1}*[-K(u + dt*du_dt)]
    // for du_dt
-   if (!T)
+   if (!T && !Toper)
    {
-      T = Add(1.0, Mmat, dt, Kmat);
       current_dt = dt;
-      T_solver.SetOperator(*T);
+      if (!partial_assembly_diff && !partial_assembly_mass)
+      {
+         T = Add(1.0, Mmat, dt, Kmat);
+         T_solver.SetOperator(*T);
+      }
+      else
+      {
+         Toper = new TimeDerivativeOperator(Moper, dt, Koper);
+         T_solver.SetOperator(*Toper);
+      }
    }
    MFEM_VERIFY(dt == current_dt, ""); // SDIRK methods use the same dt
-   Kmat.Mult(u, z);
+   Koper->Mult(u, z);
    z.Neg();
    T_solver.Mult(z, du_dt);
 }
@@ -358,21 +424,32 @@ void ConductionOperator::SetParameters(const Vector &u)
       u_alpha_gf(i) = kappa + alpha*u_alpha_gf(i);
    }
 
-   delete K;
-   K = new BilinearForm(&fespace);
-
    GridFunctionCoefficient u_coeff(&u_alpha_gf);
 
+   delete K;
+   K = new BilinearForm(&fespace);
    K->AddDomainIntegrator(new DiffusionIntegrator(u_coeff));
-   K->Assemble();
-   K->FormSystemMatrix(ess_tdof_list, Kmat);
+   if (!partial_assembly_diff)
+   {
+      K->AssembleForm(Kmat);
+      K->FormSystemOperator(ess_tdof_list, Koper);
+   }
+   else
+   {
+      K->AssembleForm(Kpaop);
+      K->FormSystemOperator(ess_tdof_list, Koper);
+   }
+
+   // re-compute on the next ImplicitSolve
    delete T;
-   T = NULL; // re-compute T on the next ImplicitSolve
+   delete Toper;
+   Toper = NULL;
+   T = NULL;
 }
 
 ConductionOperator::~ConductionOperator()
 {
-   delete T;
+   delete Toper;
    delete M;
    delete K;
 }

examples/ex1p.cpp

@@ -54,7 +54,9 @@ int main(int argc, char *argv[])
    const char *mesh_file = "../data/star.mesh";
    int order = 1;
    bool static_cond = false;
-   bool visualization = 1;
+   bool visualization = true;
+   bool use_partial_assembly = false;
+   bool use_amg = true;
 
    OptionsParser args(argc, argv);
    args.AddOption(&mesh_file, "-m", "--mesh",
@@ -62,8 +64,12 @@ int main(int argc, char *argv[])
    args.AddOption(&order, "-o", "--order",
                   "Finite element order (polynomial degree) or -1 for"
                   " isoparametric space.");
+   args.AddOption(&use_partial_assembly, "-pa", "--partial-assembly",
+                  "-no-pa", "--no-partial-assembly", "Enable partial assembly.");
    args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
                   "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&use_amg, "-pc", "--peconditioner", "-no-pc",
+                  "--no-preconditioner", "Use an algebraic multigrid preconditioner.");
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
@@ -173,30 +179,59 @@ int main(int argc, char *argv[])
    ParBilinearForm *a = new ParBilinearForm(fespace);
    a->AddDomainIntegrator(new DiffusionIntegrator(one));
 
+   // Use the global map instead
+   BilinearFormOperator A_pa(new PAIntegratorMap);
+   HypreParMatrix A_hpm;
+   Vector B, X;
+
+   if (static_cond) { a->EnableStaticCondensation(); }
    // 11. Assemble the parallel bilinear form and the corresponding linear
    //     system, applying any necessary transformations such as: parallel
    //     assembly, eliminating boundary conditions, applying conforming
    //     constraints for non-conforming AMR, static condensation, etc.
-   if (static_cond) { a->EnableStaticCondensation(); }
-   a->Assemble();
+   if (!use_partial_assembly)
+   {
+      a->AssembleForm(A_hpm);
+   }
+   else
+   {
+      a->AssembleForm(A_pa);
+   }
 
-   HypreParMatrix A;
-   Vector B, X;
+   Operator *A;
    a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
 
    if (myid == 0)
    {
-      cout << "Size of linear system: " << A.GetGlobalNumRows() << endl;
+      cout << "Size of linear system: " << A->Height() << endl;
    }
 
    // 12. Define and apply a parallel PCG solver for AX=B with the BoomerAMG
    //     preconditioner from hypre.
-   HypreSolver *amg = new HypreBoomerAMG(A);
-   HyprePCG *pcg = new HyprePCG(A);
-   pcg->SetTol(1e-12);
-   pcg->SetMaxIter(200);
-   pcg->SetPrintLevel(2);
-   pcg->SetPreconditioner(*amg);
+   Solver *pcg = NULL;
+   HypreSolver *amg = NULL;
+   if (!use_partial_assembly)
+   {
+      HyprePCG *hypre_pcg = new HyprePCG(A_hpm);
+      pcg = hypre_pcg;
+      hypre_pcg->SetTol(1e-12);
+      hypre_pcg->SetMaxIter(200);
+      hypre_pcg->SetPrintLevel(2);
+      if (use_amg)
+      {
+         amg = new HypreBoomerAMG(A_hpm);
+         hypre_pcg->SetPreconditioner(*amg);
+      }
+   }
+   else
+   {
+      CGSolver *mfem_pcg = new CGSolver(MPI_COMM_WORLD);
+      pcg = mfem_pcg;
+      mfem_pcg->SetRelTol(1e-12);
+      mfem_pcg->SetMaxIter(200);
+      mfem_pcg->SetPrintLevel(1);
+      mfem_pcg->SetOperator(*A);
+   }
    pcg->Mult(B, X);
 
    // 13. Recover the parallel grid function corresponding to X. This is the
@@ -231,6 +266,7 @@ int main(int argc, char *argv[])
    }
 
    // 16. Free the used memory.
+   delete A;
    delete pcg;
    delete amg;
    delete a;

examples/ex9.cpp

@@ -78,6 +78,7 @@ public:
 
 int main(int argc, char *argv[])
 {
+
    // 1. Parse command-line options.
    problem = 0;
    const char *mesh_file = "../data/periodic-hexagon.mesh";
@@ -178,8 +179,20 @@ int main(int argc, char *argv[])
    FunctionCoefficient inflow(inflow_function);
    FunctionCoefficient u0(u0_function);
 
+   tic_toc.Clear();
+   tic_toc.Start();
+
    BilinearForm m(&fes);
    m.AddDomainIntegrator(new MassIntegrator);
+   m.Assemble();
+   m.Finalize();
+
+   tic_toc.Stop();
+   double mass_init_time = tic_toc.RealTime(); 
+   cout << " Mass initialization time: " << mass_init_time << "s." << endl;
+   tic_toc.Clear();
+   tic_toc.Start();
+
    BilinearForm k(&fes);
    k.AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
    k.AddInteriorFaceIntegrator(
@@ -187,15 +200,20 @@ int main(int argc, char *argv[])
    k.AddBdrFaceIntegrator(
       new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
 
-   LinearForm b(&fes);
-   b.AddBdrFaceIntegrator(
-      new BoundaryFlowIntegrator(inflow, velocity, -1.0, -0.5));
 
-   m.Assemble();
-   m.Finalize();
    int skip_zeros = 0;
    k.Assemble(skip_zeros);
    k.Finalize(skip_zeros);
+
+   tic_toc.Stop();
+   double adv_init_time = tic_toc.RealTime(); 
+   cout << " Advection initialization time: " << adv_init_time << "s." << endl;
+   tic_toc.Clear();
+   tic_toc.Start();
+
+   LinearForm b(&fes);
+   b.AddBdrFaceIntegrator(
+      new BoundaryFlowIntegrator(inflow, velocity, -1.0, -0.5));
    b.Assemble();
 
    // 7. Define the initial conditions, save the corresponding grid function to
@@ -213,6 +231,10 @@ int main(int argc, char *argv[])
       u.Save(osol);
    }
 
+   tic_toc.Stop();
+   double total_init_time = mass_init_time + adv_init_time + tic_toc.RealTime(); 
+   cout << " Initialization time: " << total_init_time << "s." << endl;
+
    // Create data collection for solution output: either VisItDataCollection for
    // ascii data files, or SidreDataCollection for binary data files.
    DataCollection *dc = NULL;
@@ -270,6 +292,9 @@ int main(int argc, char *argv[])
    adv.SetTime(t);
    ode_solver->Init(adv);
 
+   tic_toc.Clear();
+   tic_toc.Start();
+
    bool done = false;
    for (int ti = 0; !done; )
    {
@@ -297,6 +322,9 @@ int main(int argc, char *argv[])
       }
    }
 
+   tic_toc.Stop();
+   cout << " Computation time: " << tic_toc.RealTime() << "s." << endl;
+
    // 9. Save the final solution. This output can be viewed later using GLVis:
    //    "glvis -m ex9.mesh -g ex9-final.gf".
    {
@@ -317,7 +345,7 @@ int main(int argc, char *argv[])
 FE_Evolution::FE_Evolution(SparseMatrix &_M, SparseMatrix &_K, const Vector &_b)
    : TimeDependentOperator(_M.Size()), M(_M), K(_K), b(_b), z(_M.Size())
 {
-   M_solver.SetPreconditioner(M_prec);
+   //M_solver.SetPreconditioner(M_prec);
    M_solver.SetOperator(M);
 
    M_solver.iterative_mode = false;
@@ -329,10 +357,31 @@ FE_Evolution::FE_Evolution(SparseMatrix &_M, SparseMatrix &_K, const Vector &_b)
 
 void FE_Evolution::Mult(const Vector &x, Vector &y) const
 {
+	/*y = 0.;
+   Vector xx(x);
+   int size = xx.Size();
+   int n = size;
+   int order = 1;
+   int dofs = (order+1)*(order+1);
+   for (int i = 0; i < n; ++i)
+   {
+      cout << "cacahuete " << i << endl;
+      xx = 0.;
+      xx(i) = 1000.;
+      // y = M^{-1} (K x + b)
+      K.Mult(xx, z);
+      for (int j = 0; j < z.Size(); ++j)
+      {
+         z(j) = abs(z(j)) < 1e-12 ? 0 : z(j);
+      }
+      z.Print(std::cout,dofs);
+      y += z;
+   }*/
    // y = M^{-1} (K x + b)
    K.Mult(x, z);
    z += b;
    M_solver.Mult(z, y);
+	// K.Mult(x, y);
 }
 
 

examples/ex9PA.cpp

@@ -0,0 +1,584 @@
+//                                MFEM Example 9
+//
+// Compile with: make ex9
+//
+// Sample runs:
+//    ex9 -m ../data/periodic-segment.mesh -p 0 -r 2 -dt 0.005
+//    ex9 -m ../data/periodic-square.mesh -p 0 -r 2 -dt 0.01 -tf 10
+//    ex9 -m ../data/periodic-hexagon.mesh -p 0 -r 2 -dt 0.01 -tf 10
+//    ex9 -m ../data/periodic-square.mesh -p 1 -r 2 -dt 0.005 -tf 9
+//    ex9 -m ../data/periodic-hexagon.mesh -p 1 -r 2 -dt 0.005 -tf 9
+//    ex9 -m ../data/amr-quad.mesh -p 1 -r 2 -dt 0.002 -tf 9
+//    ex9 -m ../data/star-q3.mesh -p 1 -r 2 -dt 0.005 -tf 9
+//    ex9 -m ../data/disc-nurbs.mesh -p 1 -r 3 -dt 0.005 -tf 9
+//    ex9 -m ../data/disc-nurbs.mesh -p 2 -r 3 -dt 0.005 -tf 9
+//    ex9 -m ../data/periodic-square.mesh -p 3 -r 4 -dt 0.0025 -tf 9 -vs 20
+//    ex9 -m ../data/periodic-cube.mesh -p 0 -r 2 -o 2 -dt 0.02 -tf 8
+//
+// Description:  This example code solves the time-dependent advection equation
+//               du/dt + v.grad(u) = 0, where v is a given fluid velocity, and
+//               u0(x)=u(0,x) is a given initial condition.
+//
+//               The example demonstrates the use of Discontinuous Galerkin (DG)
+//               bilinear forms in MFEM (face integrators), the use of explicit
+//               ODE time integrators, the definition of periodic boundary
+//               conditions through periodic meshes, as well as the use of GLVis
+//               for persistent visualization of a time-evolving solution. The
+//               saving of time-dependent data files for external visualization
+//               with VisIt (visit.llnl.gov) is also illustrated.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include <algorithm>
+
+#include "../fem/dgpabilininteg.hpp"
+#include "../fem/dgfacefunctions.hpp"
+#include "../fem/partialassemblykernel.hpp"
+
+using namespace std;
+using namespace mfem;
+
+// Choice for the problem setup. The fluid velocity, initial condition and
+// inflow boundary condition are chosen based on this parameter.
+int problem;
+
+// Velocity coefficient
+void velocity_function(const Vector &x, Vector &v);
+
+// Initial condition
+double u0_function(const Vector &x);
+
+// Inflow boundary condition
+double inflow_function(const Vector &x);
+
+// Mesh bounding box
+Vector bb_min, bb_max;
+
+
+/** A time-dependent operator for the right-hand side of the ODE. The DG weak
+    form of du/dt = -v.grad(u) is M du/dt = K u + b, where M and K are the mass
+    and advection matrices, and b describes the flow on the boundary. This can
+    be written as a general ODE, du/dt = M^{-1} (K u + b), and this class is
+    used to evaluate the right-hand side. */
+class FE_Evolution : public TimeDependentOperator
+{
+private:
+   // BilinearForm &M;
+   Operator &M;
+   Operator &K;
+   const Vector &b;
+
+   CGSolver M_solver;
+   DSmoother M_prec;
+
+   mutable Vector z;
+
+public:
+   // FE_Evolution(BilinearForm &_M, BilinearForm &_K, const Vector &_b);
+   FE_Evolution(Operator &_M, Operator &_K, const Vector &_b);
+
+   virtual void Mult(const Vector &x, Vector &y) const;
+
+   virtual ~FE_Evolution() { }
+};
+
+
+int main(int argc, char *argv[])
+{
+
+   // 1. Parse command-line options.
+   problem = 0;
+   const char *mesh_file = "../data/periodic-hexagon.mesh";
+   int ref_levels = 2;
+   int order = 3;
+   int ode_solver_type = 4;
+   double t_final = 10.0;
+   double dt = 0.01;
+   bool visualization = true;
+   bool visit = false;
+   bool binary = false;
+   int vis_steps = 5;
+
+   int precision = 8;
+   cout.precision(precision);
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&problem, "-p", "--problem",
+                  "Problem setup to use. See options in velocity_function().");
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly.");
+   args.AddOption(&order, "-o", "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&ode_solver_type, "-s", "--ode-solver",
+                  "ODE solver: 1 - Forward Euler,\n\t"
+                  "            2 - RK2 SSP, 3 - RK3 SSP, 4 - RK4, 6 - RK6.");
+   args.AddOption(&t_final, "-tf", "--t-final",
+                  "Final time; start time is 0.");
+   args.AddOption(&dt, "-dt", "--time-step",
+                  "Time step.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
+                  "--no-visit-datafiles",
+                  "Save data files for VisIt (visit.llnl.gov) visualization.");
+   args.AddOption(&binary, "-binary", "--binary-datafiles", "-ascii",
+                  "--ascii-datafiles",
+                  "Use binary (Sidre) or ascii format for VisIt data files.");
+   args.AddOption(&vis_steps, "-vs", "--visualization-steps",
+                  "Visualize every n-th timestep.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   // 2. Read the mesh from the given mesh file. We can handle geometrically
+   //    periodic meshes in this code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 3. Define the ODE solver used for time integration. Several explicit
+   //    Runge-Kutta methods are available.
+   ODESolver *ode_solver = NULL;
+   switch (ode_solver_type)
+   {
+      case 1: ode_solver = new ForwardEulerSolver; break;
+      case 2: ode_solver = new RK2Solver(1.0); break;
+      case 3: ode_solver = new RK3SSPSolver; break;
+      case 4: ode_solver = new RK4Solver; break;
+      case 6: ode_solver = new RK6Solver; break;
+      default:
+         cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
+         return 3;
+   }
+
+   // 4. Refine the mesh to increase the resolution. In this example we do
+   //    'ref_levels' of uniform refinement, where 'ref_levels' is a
+   //    command-line parameter. If the mesh is of NURBS type, we convert it to
+   //    a (piecewise-polynomial) high-order mesh.
+   for (int lev = 0; lev < ref_levels; lev++)
+   {
+      mesh->UniformRefinement();
+   }
+   if (mesh->NURBSext)
+   {
+      mesh->SetCurvature(max(order, 1));
+   }
+   if (!mesh->GetNodes())
+   {
+      mesh->SetCurvature(1);
+   }
+   mesh->GetBoundingBox(bb_min, bb_max, max(order, 1));
+
+   // 5. Define the discontinuous DG finite element space of the given
+   //    polynomial order on the refined mesh.
+   
+   DG_FECollection fec(order, dim);
+   //H1_FECollection fec(order, dim);
+   FiniteElementSpace fes(mesh, &fec);
+
+   cout << "Number of unknowns: " << fes.GetVSize() << endl;
+
+   // 6. Set up and assemble the bilinear and linear forms corresponding to the
+   //    DG discretization. The DGTraceIntegrator involves integrals over mesh
+   //    interior faces.
+   VectorFunctionCoefficient velocity(dim, velocity_function);
+   FunctionCoefficient inflow(inflow_function);
+   FunctionCoefficient u0(u0_function);
+
+   //Creating a partial assembly Kernel
+   //Maybe not the right place to initialize tensor size.
+   int ir_order = 2*order+1;
+
+   tic_toc.Clear();
+   tic_toc.Start();
+
+   //Initialization of the Mass operator
+   // BilinearFormOperator m(&fes);
+   // m.AddDomainIntegrator(new PAMassIntegrator(&fes,ir_order));
+   // m.AddDomainIntegrator(new EigenPAMassIntegrator<2>(&fes,ir_order));
+   // m.AddDomainIntegrator(new EigenPAMassIntegrator<2,EigenDomainPAK>(&fes,ir_order));
+   //BilinearForm m(&fes);
+   //m.AddIntegrator(new PADomainInt<MassEquation>(&fes,ir_order,MassEquation::ArgsEmpty{}));
+   // m.AddDomainIntegrator(new MassIntegrator());
+   // m.AddIntegrator(new PADomainInt<MassEquation,CGSolverDG>(&fes,ir_order,MassEquation::ArgsEmpty{}));
+   // m.AddIntegrator(new PADomainInt<MassEquation>(&fes,ir_order));
+   PADomainInt<MassEquation> mass(&fes,ir_order,MassEquation::ArgsEmpty{});
+   // DiagSolverDG m(fes,ir_order,mass);
+   // PACGSolver<PADomainInt<MassEquation>> m(&fes,mass);
+   DiagSolverDG prec(fes,ir_order,mass);
+   // PAPrecCGSolver<PADomainInt<MassEquation>,DiagSolverDG> m(&fes,mass,prec);
+   // CGSolverDG<PADomainInt<MassEquation>> m(fes,ir_order,mass);
+   PrecCGSolverDG<PADomainInt<MassEquation>,DiagSolverDG> m(fes,ir_order,mass,prec);
+   Operator* mo = &m;
+
+   Array<int> ess_tdof_list;
+   // SparseMatrix msp;
+   // BilinearFormOperator mbf;
+   // Operator *mo;
+   // m.AssembleForm(msp);
+   // m.AssembleForm(mbf);
+   // m.FormSystemOperator(ess_tdof_list, mo);
+
+   tic_toc.Stop();
+   double mass_init_time = tic_toc.RealTime(); 
+   cout << " Mass initialization time: " << mass_init_time << "s." << endl;
+   tic_toc.Clear();
+   tic_toc.Start();
+
+   //Initialization of the Stiffness operator
+   BilinearForm k(&fes);
+   //k.AddDomainIntegrator(new EigenPAConvectionIntegrator<2>(&fes,ir_order,velocity, -1.0));
+   // k.AddDomainIntegrator(new PAConvectionIntegrator<DummyDomainPAK>(&fes,ir_order,velocity, -1.0));
+   typename DGConvectionEquation::Args argsEq(velocity,-1.0,-0.5);
+   k.AddIntegrator(new PADomainInt<DGConvectionEquation>(&fes,ir_order,argsEq));
+   // k.AddIntegrator(new PADomainInt<DGConvectionEquation>(&fes,ir_order,velocity,-1.0));
+   // k.AddDomainIntegrator(
+   //       new PADGConvectionFaceIntegrator<DummyFacePAK>(&fes,ir_order,velocity, 1.0, -0.5));
+   // k.AddDomainIntegrator(
+   //       new PADGConvectionFaceIntegrator2<FacePAK>(&fes,ir_order,velocity, 1.0, -0.5));
+   k.AddIntegrator(new PAFaceInt<DGConvectionEquation>(&fes,ir_order,argsEq));
+   // k.AddIntegrator(new PAFaceInt<DGConvectionEquation>(&fes,ir_order,velocity, 1.0, -0.5));
+
+   BilinearFormOperator kbf;
+   Operator *ko;
+   k.AssembleForm(kbf);
+   k.FormSystemOperator(ess_tdof_list, ko);
+
+   tic_toc.Stop();
+   double adv_init_time = tic_toc.RealTime(); 
+   cout << " Advection initialization time: " << adv_init_time << "s." << endl;
+   tic_toc.Clear();
+   tic_toc.Start();
+
+   //No need to do PA
+   LinearForm b(&fes);
+   b.AddBdrFaceIntegrator(
+      new BoundaryFlowIntegrator(inflow, velocity, -1.0, -0.5));
+
+   b.Assemble();
+
+   tic_toc.Stop();
+   double total_init_time = mass_init_time + adv_init_time + tic_toc.RealTime(); 
+   cout << " Initialization time: " << total_init_time << "s." << endl;
+
+   // 7. Define the initial conditions, save the corresponding grid function to
+   //    a file and (optionally) save data in the VisIt format and initialize
+   //    GLVis visualization.
+   GridFunction u(&fes);
+   u.ProjectCoefficient(u0);
+
+   {
+      ofstream omesh("ex9.mesh");
+      omesh.precision(precision);
+      mesh->Print(omesh);
+      ofstream osol("ex9-init.gf");
+      osol.precision(precision);
+      u.Save(osol);
+   }
+
+   // Create data collection for solution output: either VisItDataCollection for
+   // ascii data files, or SidreDataCollection for binary data files.
+   DataCollection *dc = NULL;
+   if (visit)
+   {
+      if (binary)
+      {
+#ifdef MFEM_USE_SIDRE
+         dc = new SidreDataCollection("Example9", mesh);
+#else
+         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+#endif
+      }
+      else
+      {
+         dc = new VisItDataCollection("Example9", mesh);
+         dc->SetPrecision(precision);
+      }
+      dc->RegisterField("solution", &u);
+      dc->SetCycle(0);
+      dc->SetTime(0.0);
+      dc->Save();
+   }
+
+   socketstream sout;
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      sout.open(vishost, visport);
+      if (!sout)
+      {
+         cout << "Unable to connect to GLVis server at "
+              << vishost << ':' << visport << endl;
+         visualization = false;
+         cout << "GLVis visualization disabled.\n";
+      }
+      else
+      {
+         sout.precision(precision);
+         sout << "solution\n" << *mesh << u;
+         sout << "pause\n";
+         sout << flush;
+         cout << "GLVis visualization paused."
+              << " Press space (in the GLVis window) to resume it.\n";
+      }
+   }
+
+   // 8. Define the time-dependent evolution operator describing the ODE
+   //    right-hand side, and perform time-integration (looping over the time
+   //    iterations, ti, with a time-step dt).
+   FE_Evolution adv(*mo, *ko, b);
+   // FE_Evolution adv(m, k, b);
+
+   double t = 0.0;
+   adv.SetTime(t);
+   ode_solver->Init(adv);
+
+   tic_toc.Clear();
+   tic_toc.Start();
+
+   bool done = false;
+   for (int ti = 0; !done; )
+   {
+      double dt_real = min(dt, t_final - t);
+      ode_solver->Step(u, t, dt_real);
+      ti++;
+
+      //done = true;
+      done = (t >= t_final - 1e-8*dt);
+
+      if (done || ti % vis_steps == 0)
+      {
+         cout << "time step: " << ti << ", time: " << t << endl;
+
+         if (visualization)
+         {
+            sout << "solution\n" << *mesh << u << flush;
+            // sout << "screenshot\n" << "ex9-" << ti << ".png" << flush;
+         }
+
+         if (visit)
+         {
+            dc->SetCycle(ti);
+            dc->SetTime(t);
+            dc->Save();
+         }
+      }
+   }
+
+   tic_toc.Stop();
+   cout << " Computation time: " << tic_toc.RealTime() << "s." << endl;
+
+   // 9. Save the final solution. This output can be viewed later using GLVis:
+   //    "glvis -m ex9.mesh -g ex9-final.gf".
+   {
+      ofstream osol("ex9-final.gf");
+      osol.precision(precision);
+      u.Save(osol);
+   }
+
+   // 10. Free the used memory.
+   delete ode_solver;
+   delete dc;
+   //delete mo;
+   //delete ko;
+
+   return 0;
+}
+
+
+// Implementation of class FE_Evolution
+// FE_Evolution::FE_Evolution(BilinearForm &_M, BilinearForm &_K, const Vector &_b)
+//    : TimeDependentOperator(_M.Size(), 0.0), M(_M), K(_K), b(_b), z(_M.Size())
+// {
+//    //TODO have to take into account the block diagonal structure of M
+//    //M_solver.SetPreconditioner(M_prec);
+//    M_solver.SetOperator(M);
+
+//    M_solver.iterative_mode = true;
+//    M_solver.SetRelTol(1e-9);
+//    M_solver.SetAbsTol(0.0);
+//    M_solver.SetMaxIter(100);
+//    M_solver.SetPrintLevel(0);
+// }
+FE_Evolution::FE_Evolution(Operator &_M, Operator &_K, const Vector &_b)
+   : TimeDependentOperator(_M.Height(), 0.0), M(_M), K(_K), b(_b), z(_M.Height())
+{
+   //TODO have to take into account the block diagonal structure of M
+   // M_solver.SetPreconditioner(M_prec);
+   M_solver.SetOperator(M);
+
+   M_solver.iterative_mode = true;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+}
+
+void FE_Evolution::Mult(const Vector &x, Vector &y) const
+{
+   // y = 0.;
+   // Vector xx(x);
+   // int size = xx.Size();
+   // int n = size;
+   // int order = 1;
+   // int dofs = (order+1)*(order+1);
+   // for (int i = 0; i < n; ++i)
+   // {
+   //    cout << "cacahuete " << i << endl;
+   //    xx = 0.;
+   //    xx(i) = 1000.;
+   //    // y = M^{-1} (K x + b)
+   //    K.Mult(xx, z);
+   //    for (int j = 0; j < z.Size(); ++j)
+   //    {
+   //       z(j) = abs(z(j)) < 1e-12 ? 0 : z(j);
+   //    }
+   //    z.Print(std::cout,dofs);
+   //    y += z;
+   // }
+   K.Mult(x, z);
+   z += b;
+   // M_solver.Mult(z, y);
+   M.Mult(z,y);
+   // K.Mult(x, y);
+}
+
+
+// Velocity coefficient
+void velocity_function(const Vector &x, Vector &v)
+{
+   int dim = x.Size();
+
+   // map to the reference [-1,1] domain
+   Vector X(dim);
+   for (int i = 0; i < dim; i++)
+   {
+      double center = (bb_min[i] + bb_max[i]) * 0.5;
+      X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
+   }
+
+   switch (problem)
+   {
+      case 4:
+      case 0:
+      {
+         // Translations in 1D, 2D, and 3D
+         switch (dim)
+         {
+            case 1: v(0) = 1.0; break;
+            case 2: v(0) = sqrt(2./3.); v(1) = sqrt(1./3.); break;
+            // case 2: v(0) = 1+abs(X(0)); v(1) = 1+abs(X(0)); break;
+            case 3: v(0) = sqrt(3./6.); v(1) = sqrt(2./6.); v(2) = sqrt(1./6.);
+               break;
+         }
+         break;
+      }
+      case 1:
+      case 2:
+      {
+         // Clockwise rotation in 2D around the origin
+         const double w = M_PI/2;
+         switch (dim)
+         {
+            case 1: v(0) = 1.0; break;
+            case 2: v(0) = w*X(1); v(1) = -w*X(0); break;
+            case 3: v(0) = w*X(1); v(1) = -w*X(0); v(2) = 0.0; break;
+         }
+         break;
+      }
+      case 3:
+      {
+         // Clockwise twisting rotation in 2D around the origin
+         const double w = M_PI/2;
+         double d = max((X(0)+1.)*(1.-X(0)),0.) * max((X(1)+1.)*(1.-X(1)),0.);
+         d = d*d;
+         switch (dim)
+         {
+            case 1: v(0) = 1.0; break;
+            case 2: v(0) = d*w*X(1); v(1) = -d*w*X(0); break;
+            case 3: v(0) = d*w*X(1); v(1) = -d*w*X(0); v(2) = 0.0; break;
+         }
+         break;
+      }
+   }
+}
+
+// Initial condition
+double u0_function(const Vector &x)
+{
+   int dim = x.Size();
+
+   // map to the reference [-1,1] domain
+   Vector X(dim);
+   for (int i = 0; i < dim; i++)
+   {
+      double center = (bb_min[i] + bb_max[i]) * 0.5;
+      X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
+   }
+
+   switch (problem)
+   {
+      case 0:
+
+      case 1:
+      {
+         switch (dim)
+         {
+            case 1:
+               return exp(-40.*pow(X(0)-0.5,2));
+            case 2:
+            case 3:
+            {
+               double rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
+               // double rx = 0.05, ry = 0.05, cx = -0., cy = -0.5, w = 10.;
+               if (dim == 3)
+               {
+                  const double s = (1. + 0.25*cos(2*M_PI*X(2)));
+                  rx *= s;
+                  ry *= s;
+               }
+               return ( erfc(w*(X(0)-cx-rx))*erfc(-w*(X(0)-cx+rx)) *
+                        erfc(w*(X(1)-cy-ry))*erfc(-w*(X(1)-cy+ry)) )/16;
+            }
+         }
+      }
+      case 2:
+      {
+         double x_ = X(0), y_ = X(1), rho, phi;
+         rho = hypot(x_, y_);
+         phi = atan2(y_, x_);
+         return pow(sin(M_PI*rho),2)*sin(3*phi);
+      }
+      case 3:
+      {
+         const double f = M_PI;
+         return sin(f*X(0))*sin(f*X(1));
+      }
+      case 4:
+      {
+         return exp( -40*( X(0)*X(0) + X(1)*X(1) + X(2)*X(2) ) );
+      }
+   }
+   return 0.0;
+}
+
+// Inflow boundary condition (zero for the problems considered in this example)
+double inflow_function(const Vector &x)
+{
+   switch (problem)
+   {
+      case 0:
+      case 1:
+      case 2:
+      case 3: return 0.0;
+   }
+   return 0.0;
+}

examples/makefile

@@ -21,8 +21,8 @@ CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
-SEQ_EXAMPLES = ex1 ex2 ex3 ex4 ex5 ex6 ex7 ex8 ex9 ex10 ex14 ex15 ex16 ex17\
-  ex18 ex19
+SEQ_EXAMPLES = ex1 ex2 ex3 ex4 ex5 ex6 ex7 ex8 ex9 ex9PA ex10 ex14 ex15 ex16\
+ ex17 ex18 ex19
 PAR_EXAMPLES = ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex8p ex9p ex10p ex11p ex12p\
  ex13p ex14p ex15p ex16p ex17p ex18p ex19p
 

fem/bilinearform.cpp

@@ -68,6 +68,8 @@ BilinearForm::BilinearForm (FiniteElementSpace * f)
    fes = f;
    sequence = f->GetSequence();
    mat = mat_e = NULL;
+   oper = NULL;
+   oper_type = MFEM_SPARSEMAT;
    extern_bfs = 0;
    element_matrices = NULL;
    static_cond = NULL;
@@ -85,6 +87,8 @@ BilinearForm::BilinearForm (FiniteElementSpace * f, BilinearForm * bf, int ps)
    fes = f;
    sequence = f->GetSequence();
    mat_e = NULL;
+   oper = NULL;
+   oper_type = MFEM_SPARSEMAT;
    extern_bfs = 1;
    element_matrices = NULL;
    static_cond = NULL;
@@ -235,6 +239,11 @@ void BilinearForm::AddBdrFaceIntegrator(BilinearFormIntegrator *bfi,
    bfbfi_marker.Append(&bdr_marker);
 }
 
+void BilinearForm::AddIntegrator(LinearFESpaceIntegrator *bfi)
+{
+   fesi.Append(bfi);
+}
+
 void BilinearForm::ComputeElementMatrix(int i, DenseMatrix &elmat)
 {
    if (element_matrices)
@@ -320,6 +329,8 @@ void BilinearForm::Assemble (int skip_zeros)
       AllocMat();
    }
 
+   oper_type = MFEM_SPARSEMAT;
+
 #ifdef MFEM_USE_OPENMP
    int free_element_matrices = 0;
    if (!element_matrices)
@@ -505,6 +516,21 @@ void BilinearForm::Assemble (int skip_zeros)
 #endif
 }
 
+void BilinearForm::AssembleForm(SparseMatrix &A, int skip_zeros)
+{
+   Assemble(skip_zeros);
+   oper = &A;
+   oper_type = MFEM_SPARSEMAT;
+}
+
+void BilinearForm::AssembleForm(BilinearFormOperator &A)
+{
+   A.Assemble(this);
+   oper = &A;
+   oper_type = MFEM_FORMOPER;
+}
+
+
 void BilinearForm::ConformingAssemble()
 {
    // Do not remove zero entries to preserve the symmetric structure of the
@@ -548,6 +574,7 @@ void BilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list,
    const SparseMatrix *P = fes->GetConformingProlongation();
 
    FormSystemMatrix(ess_tdof_list, A);
+   oper = &A;
 
    // Transform the system and perform the elimination in B, based on the
    // essential BC values from x. Restrict the BC part of x in X, and set the
@@ -606,6 +633,29 @@ void BilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list,
    }
 }
 
+void BilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x, Vector &b,
+                                    Operator * &A, Vector &X, Vector &B,
+                                    int copy_interior)
+{
+   if (oper_type == MFEM_SPARSEMAT)
+   {
+      SparseMatrix &Amat = static_cast<SparseMatrix&>(*oper);
+      FormLinearSystem(ess_tdof_list, x, b, Amat,
+                       X, B, copy_interior);
+      SparseMatrix *M = new SparseMatrix;
+      M->MakeRef(Amat);
+      A = M;
+   }
+   else if (oper_type == MFEM_FORMOPER)
+   {
+      oper->FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
+   }
+   else
+   {
+      mfem_error("Not supported.");
+   }
+}
+
 void BilinearForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
                                     SparseMatrix &A)
 {
@@ -643,6 +693,24 @@ void BilinearForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
    }
 }
 
+void BilinearForm::FormSystemOperator(const Array<int> &ess_tdof_list,
+                                      Operator * &A)
+{
+   if (oper_type == MFEM_SPARSEMAT)
+   {
+      FormSystemMatrix(ess_tdof_list, static_cast<SparseMatrix&>(*oper));
+      A = oper;
+   }
+   else if (oper_type == MFEM_FORMOPER)
+   {
+      A = oper;
+   }
+   else
+   {
+      mfem_error("Not supported.");
+   }
+}
+
 void BilinearForm::RecoverFEMSolution(const Vector &X,
                                       const Vector &b, Vector &x)
 {
@@ -929,10 +997,10 @@ BilinearForm::~BilinearForm()
       for (k=0; k < bbfi.Size(); k++) { delete bbfi[k]; }
       for (k=0; k < fbfi.Size(); k++) { delete fbfi[k]; }
       for (k=0; k < bfbfi.Size(); k++) { delete bfbfi[k]; }
+      for (k=0; k < fesi.Size(); k++) { delete fesi[k]; }
    }
 }
 
-
 MixedBilinearForm::MixedBilinearForm (FiniteElementSpace *tr_fes,
                                       FiniteElementSpace *te_fes)
    : Matrix(te_fes->GetVSize(), tr_fes->GetVSize())
@@ -954,7 +1022,7 @@ const double & MixedBilinearForm::Elem (int i, int j) const
 
 void MixedBilinearForm::Mult (const Vector & x, Vector & y) const
 {
-   mat -> Mult (x, y);
+   oper -> Mult (x, y);
 }
 
 void MixedBilinearForm::AddMult (const Vector & x, Vector & y,
@@ -1006,6 +1074,12 @@ void MixedBilinearForm::AddTraceFaceIntegrator (BilinearFormIntegrator * bfi)
    skt.Append (bfi);
 }
 
+void MixedBilinearForm::AddIntegrator(LinearFESpaceIntegrator *integ)
+{
+   fesi.Append(integ);
+}
+
+
 void MixedBilinearForm::Assemble (int skip_zeros)
 {
    int i, k;
@@ -1089,6 +1163,21 @@ void MixedBilinearForm::Assemble (int skip_zeros)
          }
       }
    }
+
+   oper = mat;
+}
+
+void MixedBilinearForm::AssembleForm(SparseMatrix &A, int skip_zeros)
+{
+   Assemble(skip_zeros);
+   oper = mat;
+   A.MakeRef(*mat);
+}
+
+void MixedBilinearForm::AssembleForm(BilinearFormOperator &A, int skip_zeros)
+{
+   A.Assemble(this);
+   oper = &A;
 }
 
 void MixedBilinearForm::ConformingAssemble()
@@ -1113,6 +1202,8 @@ void MixedBilinearForm::ConformingAssemble()
       mat = RAP;
    }
 
+   oper = mat;
+
    height = mat->Height();
    width = mat->Width();
 }
@@ -1182,6 +1273,7 @@ MixedBilinearForm::~MixedBilinearForm()
    for (i = 0; i < dom.Size(); i++) { delete dom[i]; }
    for (i = 0; i < bdr.Size(); i++) { delete bdr[i]; }
    for (i = 0; i < skt.Size(); i++) { delete skt[i]; }
+   for (i = 0; i < fesi.Size(); i++) { delete fesi[i]; }
 }
 
 

fem/bilinearform.hpp

@@ -24,17 +24,27 @@
 namespace mfem
 {
 
+// Forward declare BilinearFormOperator
+class BilinearFormOperator;
+
 /** Class for bilinear form - "Matrix" with associated FE space and
     BLFIntegrators. */
 class BilinearForm : public Matrix
 {
 protected:
+   // TODO remove mat
    /// Sparse matrix to be associated with the form.
    SparseMatrix *mat;
 
    /// Matrix used to eliminate b.c.
    SparseMatrix *mat_e;
 
+   /// Generic operator associated with the form.
+   Operator *oper;
+
+   /// Operator type.
+   enum Type oper_type;
+
    /// FE space on which the form lives.
    FiniteElementSpace *fes;
 
@@ -58,6 +68,9 @@ protected:
    Array<BilinearFormIntegrator*> bfbfi;
    Array<Array<int>*>             bfbfi_marker;
 
+   /// Set of fespace integrators (does not matter what type)
+   Array<LinearFESpaceIntegrator*> fesi;
+
    DenseMatrix elemmat;
    Array<int>  vdofs;
 
@@ -83,7 +96,9 @@ protected:
    BilinearForm() : Matrix (0)
    {
       fes = NULL; sequence = -1;
-      mat = mat_e = NULL; extern_bfs = 0; element_matrices = NULL;
+      mat = mat_e = NULL;
+      oper = NULL; oper_type = MFEM_SPARSEMAT;
+      extern_bfs = 0; element_matrices = NULL;
       static_cond = NULL; hybridization = NULL;
       precompute_sparsity = 0;
       diag_policy = DIAG_KEEP;
@@ -151,6 +166,8 @@ public:
 
    Array<BilinearFormIntegrator*> *GetBFBFI() { return &bfbfi; }
 
+   Array<LinearFESpaceIntegrator*> *GetFESI() { return &fesi; }
+
    const double &operator()(int i, int j) { return (*mat)(i,j); }
 
    /// Returns reference to a_{ij}.
@@ -232,6 +249,9 @@ public:
    /// Adds new boundary Face Integrator.
    void AddBdrFaceIntegrator(BilinearFormIntegrator *bfi);
 
+   /// Adds a LinearFESpaceIntegrator.
+   void AddIntegrator(LinearFESpaceIntegrator *integ);
+
    /** @brief Adds new boundary Face Integrator, restricted to specific boundary
        attributes. */
    void AddBdrFaceIntegrator(BilinearFormIntegrator *bfi,
@@ -246,6 +266,9 @@ public:
    /// Assembles the form i.e. sums over all domain/bdr integrators.
    void Assemble(int skip_zeros = 1);
 
+   void AssembleForm(BilinearFormOperator &A);
+   void AssembleForm(SparseMatrix &A, int skip_zeros = 1);
+
    /// Get the finite element space prolongation matrix
    virtual const Operator *GetProlongation() const
    { return fes->GetConformingProlongation(); }
@@ -278,10 +301,17 @@ public:
 
        NOTE: If there are no transformations, @a X simply reuses the data of
              @a x. */
+   void FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x, Vector &b,
+                         Operator * &A, Vector &X, Vector &B,
+                         int copy_interior = 0);
+
    void FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x, Vector &b,
                          SparseMatrix &A, Vector &X, Vector &B,
                          int copy_interior = 0);
 
+   /// Form the linear system matrix A, see FormLinearSystem for details.
+   void FormSystemOperator(const Array<int> &ess_tdof_list, Operator * &Aoper);
+
    /// Form the linear system matrix A, see FormLinearSystem() for details.
    void FormSystemMatrix(const Array<int> &ess_tdof_list, SparseMatrix &A);
 
@@ -375,6 +405,7 @@ public:
    virtual ~BilinearForm();
 };
 
+
 /**
    Class for assembling of bilinear forms `a(u,v)` defined on different
    trial and test spaces. The assembled matrix `A` is such that
@@ -394,12 +425,14 @@ class MixedBilinearForm : public Matrix
 {
 protected:
    SparseMatrix *mat;
+   Operator *oper;
 
    FiniteElementSpace *trial_fes, *test_fes;
 
    Array<BilinearFormIntegrator*> dom;
    Array<BilinearFormIntegrator*> bdr;
    Array<BilinearFormIntegrator*> skt; // trace face integrators
+   Array<LinearFESpaceIntegrator*> fesi;
 
 public:
    MixedBilinearForm (FiniteElementSpace *tr_fes,
@@ -442,16 +475,24 @@ public:
        two adjacent volume FEs from the test space. */
    void AddTraceFaceIntegrator (BilinearFormIntegrator * bfi);
 
+   /// Add an FESpaceIntegrator
+   void AddIntegrator (LinearFESpaceIntegrator *integ);
+
    Array<BilinearFormIntegrator*> *GetDBFI() { return &dom; }
 
    Array<BilinearFormIntegrator*> *GetBBFI() { return &bdr; }
 
    Array<BilinearFormIntegrator*> *GetTFBFI() { return &skt; }
 
+   Array<LinearFESpaceIntegrator*> *GetFESI() { return &fesi; }
+
    void operator= (const double a) { *mat = a; }
 
    void Assemble (int skip_zeros = 1);
 
+   void AssembleForm(BilinearFormOperator &A, int skip_zeros = 1);
+   void AssembleForm(SparseMatrix &A, int skip_zeros = 1);
+
    /** For partially conforming trial and/or test FE spaces, complete the
        assembly process by performing A := P2^t A P1 where A is the internal
        sparse matrix; P1 and P2 are the conforming prolongation matrices of the
@@ -469,6 +510,10 @@ public:
 
    void Update();
 
+   FiniteElementSpace *TrialFESpace() const { return trial_fes; }
+
+   FiniteElementSpace *TestFESpace() const { return test_fes; }
+
    virtual ~MixedBilinearForm();
 };
 

fem/bilinearformoper.cpp

@@ -0,0 +1,316 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+// Implementation of BilinearFormOperator
+
+#include "fem.hpp"
+
+namespace mfem
+{
+
+static void BuildDofMaps(FiniteElementSpace *fespace, Array<int> *&off,
+                         Array<int> *&ind)
+{
+   // Get the total size without vdim
+   int size = 0;
+   const int vdim = fespace->GetVDim();
+   for (int e = 0; e < fespace->GetNE(); e++)
+   {
+      const FiniteElement *fe = fespace->GetFE(e);
+      size += fe->GetDof();
+   }
+   const int local_size = size * vdim;
+   const int global_size = fespace->GetVSize();
+
+   // Now we can allocate and fill the global map
+   off = new Array<int>(global_size + 1);
+   ind = new Array<int>(local_size);
+
+   Array<int> &offsets = *off;
+   Array<int> &indices = *ind;
+
+   Array<int> global_map(local_size);
+   Array<int> elem_vdof;
+
+   int offset = 0;
+   for (int e = 0; e < fespace->GetNE(); e++)
+   {
+      const FiniteElement *fe = fespace->GetFE(e);
+      const int dofs = fe->GetDof();
+      const int vdofs = dofs * vdim;
+      const TensorBasisElement *tfe = dynamic_cast<const TensorBasisElement *>(fe);
+      const Array<int> &dof_map = tfe->GetDofMap();
+
+      fespace->GetElementVDofs(e, elem_vdof);
+
+      if (dof_map.Size()==0)
+      {
+         for (int vd = 0; vd < vdim; vd++)
+            for (int i = 0; i < vdofs; i++)
+            {
+               global_map[offset + dofs*vd + i] = elem_vdof[dofs*vd + i];
+            }
+      }else{
+         for (int vd = 0; vd < vdim; vd++)
+            for (int i = 0; i < vdofs; i++)
+            {
+               global_map[offset + dofs*vd + i] = elem_vdof[dofs*vd + dof_map[i]];
+            }
+      }
+      offset += vdofs;
+   }
+
+   // Store and use a set of offsets and indices instead of this map
+
+   // Zero the offset vector
+   offsets = 0;
+
+   // Keep track of how many local dof point to its global dof
+   // Count how many times each dof gets hit
+   for (int i = 0; i < local_size; i++)
+   {
+      const int g = global_map[i];
+      ++offsets[g + 1];
+   }
+   // Aggregate the offsets
+   for (int i = 1; i <= global_size; i++)
+   {
+      offsets[i] += offsets[i - 1];
+   }
+
+   for (int i = 0; i < local_size; i++)
+   {
+      const int g = global_map[i];
+      indices[offsets[g]++] = i;
+   }
+
+   // Shift the offset vector back by one, since it was used as a
+   // counter above.
+   for (int i = global_size; i > 0; i--)
+   {
+      offsets[i] = offsets[i - 1];
+   }
+   offsets[0] = 0;
+}
+
+BilinearFormOperator::BilinearFormOperator(IntegratorMap *_map)
+   : bf(NULL), mbf(NULL),
+     trial_fes(NULL), test_fes(NULL),
+     trial_gs(false), test_gs(false),
+     trial_offsets(NULL), trial_indices(NULL),
+     test_offsets(NULL), test_indices(NULL),
+     X(NULL), Y(NULL),
+     map(_map) { }
+
+BilinearFormOperator::~BilinearFormOperator()
+{
+   delete map;
+   Clear();
+}
+
+void BilinearFormOperator::Assemble(BilinearForm *_bf)
+{
+   if (_bf != bf)
+   {
+      bf = _bf;
+      height = bf->Height();
+      width = bf->Width();
+
+      Init(bf->FESpace(), NULL);
+
+      // Delete the old integrator list -- Note that this does not
+      // delete the integrators themselves (since the original
+      // bilinear form owns these)
+      lfesi.DeleteAll();
+
+      // Add the integrators from bf->fesi
+      Array<LinearFESpaceIntegrator*> &other_fesi = *(bf->GetFESI());
+      for (int i = 0; i < other_fesi.Size(); i++)
+      {
+         lfesi.Append(other_fesi[i]);
+      }
+
+      if (map)
+      {
+         Array<BilinearFormIntegrator*> &dbfi = *(bf->GetDBFI());
+         for (int i = 0; i < dbfi.Size(); i++)
+         {
+            lfesi.Append(map->DomainIntegrator(dbfi[i]));
+         }
+      }
+   }
+
+   Assemble();
+}
+
+void BilinearFormOperator::Assemble(MixedBilinearForm *_mbf)
+{
+   if (_mbf != mbf)
+   {
+      mbf = _mbf;
+      height = bf->Height();
+      width = bf->Width();
+
+      Init(mbf->TrialFESpace(), mbf->TestFESpace());
+
+      // Delete the old integrator list -- Note that this does not
+      // delete the integrators themselves (since the original mixed
+      // bilinear form owns these)
+      lfesi.DeleteAll();
+
+      // Add the integrators from mbf->fesi
+      Array<LinearFESpaceIntegrator*> &other_fesi = *(mbf->GetFESI());
+      for (int i = 0; i < other_fesi.Size(); i++)
+      {
+         lfesi.Append(other_fesi[i]);
+      }
+
+      if (map)
+      {
+         Array<BilinearFormIntegrator*> &dbfi = *(mbf->GetDBFI());
+         for (int i = 0; i < dbfi.Size(); i++)
+         {
+            lfesi.Append(map->DomainIntegrator(dbfi[i]));
+         }
+      }
+   }
+
+   Assemble();
+}
+
+void BilinearFormOperator::Assemble()
+{
+   // Linear assembly
+   for (int i = 0; i < lfesi.Size(); i++)
+   {
+      lfesi[i]->Assemble(trial_fes, test_fes);
+   }
+}
+
+void BilinearFormOperator::Clear()
+{
+   delete trial_offsets;
+   delete trial_indices;
+
+   if (test_fes)
+   {
+      delete test_offsets;
+      delete test_indices;
+   }
+
+   if (trial_gs) delete X;
+   if (test_gs) delete Y;
+}
+
+void BilinearFormOperator::Init(FiniteElementSpace *_trial_fes,
+                           FiniteElementSpace *_test_fes)
+{
+   if ((_trial_fes != trial_fes) || (_test_fes != test_fes))
+   {
+      // Clear before recreating
+      Clear();
+
+      trial_fes = _trial_fes;
+      test_fes = _test_fes;
+      BuildDofMaps(trial_fes, trial_offsets, trial_indices);
+
+      if (test_fes != NULL)
+      {
+         BuildDofMaps(test_fes, test_offsets, test_indices);
+      }
+      else
+      {
+         // Point to the trial offsets and indices
+         test_offsets = trial_offsets;
+         test_indices = trial_indices;
+      }
+
+      X = new Vector(trial_indices->Size());
+      Y = new Vector(test_indices->Size());
+   }
+
+   const FiniteElementSpace *actual_test_fes =
+      (test_fes != NULL) ? test_fes : trial_fes;
+   trial_gs = test_gs = true;
+   if (dynamic_cast<const L2_FECollection *>(trial_fes->FEColl()))
+   {
+      trial_gs = test_gs = false;
+   }
+   else if (dynamic_cast<const L2_FECollection *>(actual_test_fes->FEColl()))
+   {
+      test_gs = false;
+   }
+}
+
+
+void BilinearFormOperator::LToEVector(const Array<int> &offsets,
+                                      const Array<int> &indices,
+                                      const Vector &v, Vector &V) const
+{
+   const int size = v.Size();
+   for (int i = 0; i < size; i++)
+   {
+      const int offset = offsets[i];
+      const int next_offset = offsets[i + 1];
+      const double dof_value = v(i);
+      for (int j = offset; j < next_offset; j++) { V(indices[j]) = dof_value; }
+   }
+}
+
+void BilinearFormOperator::EToLVector(const Array<int> &offsets,
+                                      const Array<int> &indices,
+                                      const Vector &V, Vector &v) const
+{
+   // NOTE: This method ADDS to the output v
+   const int size = v.Size();
+   for (int i = 0; i < size; i++)
+   {
+      const int offset = offsets[i];
+      const int next_offset = offsets[i + 1];
+      double dof_value = 0;
+      for (int j = offset; j < next_offset; j++) { dof_value += V(indices[j]); }
+      v(i) += dof_value;
+   }
+}
+
+void BilinearFormOperator::AddMult(const Vector &x, Vector &y) const
+{
+   if (trial_gs) { LToEVector(*trial_offsets, *trial_indices, x, *X); }
+   else { X = const_cast<Vector *>(&x); }
+
+   if (!test_gs) { Y = &y; }
+
+   *Y = 0.0;
+   for (int i = 0; i < lfesi.Size(); i++) lfesi[i]->AddMult(*X, *Y);
+   for (int i = 0; i < nlfesi.Size(); i++) nlfesi[i]->AddMult(*X, *Y);
+
+   if (test_gs) { EToLVector(*test_offsets, *test_indices, *Y, y); }
+}
+
+
+void BilinearFormOperator::AddMultTranspose(const Vector &x, Vector &y) const
+{
+   if (test_gs) { LToEVector(*test_offsets, *test_indices, x, *X); }
+   else { X = const_cast<Vector *>(&x); }
+
+   if (!trial_gs) { Y = &y; }
+
+   *Y = 0.0;
+   for (int i = 0; i < lfesi.Size(); i++) lfesi[i]->AddMultTranspose(*X, *Y);
+   for (int i = 0; i < nlfesi.Size(); i++) nlfesi[i]->AddMultTranspose(*X, *Y);
+
+   if (trial_gs) { EToLVector(*trial_offsets, *trial_indices, *Y, y); }
+}
+
+void BilinearFormOperator::Mult(const Vector &x, Vector &y) const
+{ y = 0.0; AddMult(x, y); }
+
+}

fem/bilinearformoper.hpp

@@ -0,0 +1,96 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+//
+// Defines the general object for the abstraction of bilinear and
+// nonlinear forms.
+
+#ifndef MFEM_BILINEARFORMOPER
+#define MFEM_BILINEARFORMOPER
+
+#include "../config/config.hpp"
+#include "../linalg/linalg.hpp"
+#include "fespace.hpp"
+#include "nonlininteg.hpp"
+
+namespace mfem
+{
+
+struct IntegratorMap
+{
+   virtual LinearFESpaceIntegrator *DomainIntegrator(BilinearFormIntegrator *integ) const
+   { mfem_error("Not supported."); return NULL; }
+   virtual LinearFESpaceIntegrator *InteriorFaceIntegrator(BilinearFormIntegrator *integ) const
+   { mfem_error("Not supported."); return NULL; }
+   virtual LinearFESpaceIntegrator *BdrFaceIntegrator(BilinearFormIntegrator *integ) const
+   { mfem_error("Not supported."); return NULL; }
+
+   virtual NonlinearFESpaceIntegrator *DomainIntegrator(NonlinearFormIntegrator *integ) const
+   { mfem_error("Not supported."); return NULL; }
+
+   virtual ~IntegratorMap() { }
+};
+
+class BilinearFormOperator : public Operator
+{
+protected:
+   BilinearForm *bf;  // Do not own
+   MixedBilinearForm *mbf;  // Do not own
+
+   FiniteElementSpace *trial_fes;  // Do not own
+   FiniteElementSpace *test_fes;  // Do not own
+   bool trial_gs, test_gs;
+
+   Array<int> *trial_offsets, *trial_indices;
+   Array<int> *test_offsets, *test_indices;
+   mutable Vector *X;
+   mutable Vector *Y;
+
+   Array<LinearFESpaceIntegrator*> lfesi;
+   Array<NonlinearFESpaceIntegrator*> nlfesi;
+
+   IntegratorMap *map;
+
+   // Convert between vector types before calling Mult.
+   void LToEVector(const Array<int> &offsets, const Array<int> &indices,
+                   const Vector &v, Vector &V) const;
+   void EToLVector(const Array<int> &offsets, const Array<int> &indices,
+                   const Vector &V, Vector &v) const;
+
+   void Clear();
+
+   void Init(FiniteElementSpace *_trial_fes, FiniteElementSpace *_test_fes);
+
+public:
+   // Create an empty object or assemble what is needed by the
+   // bilinear form integrators to later compute the action.
+   BilinearFormOperator(IntegratorMap *_map = NULL);
+   ~BilinearFormOperator();
+
+   void Assemble();
+   void Assemble(BilinearForm *bf);
+   void Assemble(MixedBilinearForm *bf);
+
+   /// Perform the action of the bilinear form on a vector and set y.
+   virtual void Mult(const Vector &x, Vector &y) const;
+
+   virtual const Operator *GetProlongation() const { return trial_fes->GetProlongationMatrix(); }
+   virtual const Operator *GetRestriction() const  { return trial_fes->GetRestrictionMatrix(); }
+
+   /// Perform the action of the bilinear form on a vector and add to y.
+   void AddMult(const Vector &x, Vector &y) const;
+
+   /// Perform the (transposed) action of the bilinear form on a vector and add to y.
+   void AddMultTranspose(const Vector &x, Vector &y) const;
+};
+
+}
+
+#endif

fem/bilininteg.cpp

@@ -61,6 +61,19 @@ void BilinearFormIntegrator::AssembleElementVector(
               "   is not implemented for this class.");
 }
 
+void BilinearFormIntegrator::AssembleMult(
+   const Vector &x, Vector& y)
+{
+   mfem_error("BilinearFormIntegrator::AssembleMult\n"
+              "   is not implemented fot this class.");
+}
+
+void BilinearFormIntegrator::AssembleMultTranspose(
+   const Vector &x, Vector& y)
+{
+   mfem_error("BilinearFormIntegrator::AssembleMultTranspose\n"
+              "   is not implemented fot this class.");
+}
 
 void TransposeIntegrator::AssembleElementMatrix (
    const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
@@ -697,7 +710,7 @@ double DiffusionIntegrator::ComputeFluxEnergy
       }
       else
       {
-         MQ->Eval(mq, Trans, ip);
+        MQ->Eval(mq, Trans, ip);
          energy += w * mq.InnerProduct(pointflux, pointflux);
       }
 
@@ -716,7 +729,6 @@ double DiffusionIntegrator::ComputeFluxEnergy
    return energy;
 }
 
-
 void MassIntegrator::AssembleElementMatrix
 ( const FiniteElement &el, ElementTransformation &Trans,
   DenseMatrix &elmat )
@@ -807,7 +819,6 @@ void MassIntegrator::AssembleElementMatrix2(
    }
 }
 
-
 void BoundaryMassIntegrator::AssembleFaceMatrix(
    const FiniteElement &el1, const FiniteElement &el2,
    FaceElementTransformations &Trans, DenseMatrix &elmat)
@@ -2054,7 +2065,6 @@ void DivDivIntegrator::AssembleElementMatrix(
    }
 }
 
-
 void VectorDiffusionIntegrator::AssembleElementMatrix(
    const FiniteElement &el,
    ElementTransformation &Trans,
@@ -2298,7 +2308,6 @@ void DGTraceIntegrator::AssembleFaceMatrix(const FiniteElement &el1,
       }
       ir = &IntRules.Get(Trans.FaceGeom, order);
    }
-
    for (int p = 0; p < ir->GetNPoints(); p++)
    {
       const IntegrationPoint &ip = ir->IntPoint(p);

fem/bilininteg.hpp

@@ -80,6 +80,17 @@ public:
                                     Vector &flux, Vector *d_energy = NULL)
    { return 0.0; }
 
+   /** Assemble any element or face-specific terms required for the
+       action with the bilinear form integrator. Later applied with
+       AssembleVector. */
+   virtual void AssembleOperator(const FiniteElementSpace *trial_fes,
+                                 const FiniteElementSpace *test_fes) { }
+
+   /** Compute `y = A * x` where A is the bilinear form integrator for
+       all elements/faces. */
+   virtual void AssembleMult(const Vector &fun, Vector &vect);
+   virtual void AssembleMultTranspose(const Vector &fun, Vector &vect);
+
    virtual ~BilinearFormIntegrator() { }
 };
 
@@ -1622,6 +1633,11 @@ public:
    virtual double ComputeFluxEnergy(const FiniteElement &fluxelem,
                                     ElementTransformation &Trans,
                                     Vector &flux, Vector *d_energy = NULL);
+
+   // Friend partial assembly version so it has access to the coefficients.
+   friend class PADiffusionIntegrator;
+   // TODO: Add a GetPAIntegrator method here
+   // PAIntegrator* GetPAIntegrator(type);
 };
 
 /** Class for local mass matrix assembling a(u,v) := (Q u, v) */
@@ -1649,6 +1665,9 @@ public:
                                        const FiniteElement &test_fe,
                                        ElementTransformation &Trans,
                                        DenseMatrix &elmat);
+
+   // Friend partial assembly version so it has access to the coefficients.
+   friend class PAMassIntegrator;
 };
 
 class BoundaryMassIntegrator : public MassIntegrator
@@ -1742,6 +1761,9 @@ public:
                                        const FiniteElement &test_fe,
                                        ElementTransformation &Trans,
                                        DenseMatrix &elmat);
+
+   // Friend partial assembly version so it has access to the coefficients.
+   friend class PAMassIntegrator;
 };
 
 

fem/dgfacefunctions.cpp

@@ -0,0 +1,760 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+//This file contains useful functions to compute fluxes for DG methods.
+
+#include <vector>
+#include "fem.hpp"
+#include "dalg.hpp"
+
+using std::vector;
+using std::pair;
+
+namespace mfem
+{
+
+/**
+*	Returns the canonical coordinate vectors e_1 and e_2.
+*/
+void getBaseVector2D(Vector& e1, Vector& e2)
+{
+	e1.SetSize(2);
+	e1(0) = 1;
+	e1(1) = 0;
+	e2.SetSize(2);
+	e2(0) = 0;
+	e2(1) = 1;
+}
+
+/**
+*	Returns the canonical coordinate vectors e_1, e_2 and e_3.
+*/
+void getBaseVector3D(Vector& e1, Vector& e2, Vector& e3)
+{
+	e1.SetSize(3);
+	e1(0) = 1;
+	e1(1) = 0;
+	e1(2) = 0;
+	e2.SetSize(3);
+	e2(0) = 0;
+	e2(1) = 1;
+	e2(2) = 0;
+	e3.SetSize(3);
+	e3(0) = 0;
+	e3(1) = 0;
+	e3(2) = 1;
+}
+
+/** 
+*   A function that initialize the local coordinate base for a face with
+*   indice face_ind.
+*   This returns the local face coordinate base expressed in reference
+*   element coordinate.
+*/
+// Highly dependent of the node ordering from geom.cpp
+void InitFaceCoord2D(const int face_id, IntMatrix& base)
+{
+	//Vector e1,e2;
+	//getBaseVector2D(e1,e2);
+	base.Zero();
+	switch(face_id)
+	{
+		case 0://SOUTH
+			base(0,0)= 1;//base.SetCol(0, e1);
+			base(1,1)=-1;//base.SetCol(1,-e2);
+			break;
+		case 1://EAST
+			base(1,0)= 1;//base.SetCol(0, e2);
+			base(0,1)= 1;//base.SetCol(1, e1);
+			break;
+		case 2://NORTH
+			base(0,0)=-1;//base.SetCol(0,-e1);
+			base(1,1)= 1;//base.SetCol(1, e2);
+			break;
+		case 3://WEST
+			base(1,0)=-1;//base.SetCol(0,-e2);
+			base(0,1)= 1;//base.SetCol(1, e1);
+			break;
+		default:
+			mfem_error("The face_ind exceeds the number of faces in this dimension.");
+			break;
+	}
+}
+
+// Highly dependent of the node ordering from geom.cpp
+void InitFaceCoord3D(const int face_id, IntMatrix& base)
+{
+	//Vector e1,e2,e3;
+	//getBaseVector3D(e1,e2,e3);
+	base.Zero();
+	switch(face_id)
+	{
+		case 0://BOTTOM
+			base(0,0)= 1;//base.SetCol(0, e1);
+			base(1,1)=-1;//base.SetCol(1,-e2);
+			base(2,2)=-1;//base.SetCol(2,-e3);
+			break;		
+		case 1://SOUTH
+			base(0,0)= 1;//base.SetCol(0, e1);
+			base(2,1)= 1;//base.SetCol(1, e3);
+			base(1,2)=-1;//base.SetCol(2,-e2);
+			break;
+		case 2://EAST
+			base(1,0)= 1;//base.SetCol(0, e2);
+			base(2,1)= 1;//base.SetCol(1, e3);
+			base(0,2)= 1;//base.SetCol(2, e1);
+			break;
+		case 3://NORTH
+			base(0,0)=-1;//base.SetCol(0,-e1);
+			base(2,1)= 1;//base.SetCol(1, e3);
+			base(1,2)= 1;//base.SetCol(2, e2);
+			break;
+		case 4://WEST
+			base(1,0)=-1;//base.SetCol(0,-e2);
+			base(2,1)= 1;//base.SetCol(1, e3);
+			base(0,2)=-1;//base.SetCol(2,-e1);
+			break;
+		case 5://TOP
+			base(0,0)= 1;//base.SetCol(0, e1);
+			base(1,1)= 1;//base.SetCol(1, e2);
+			base(2,2)= 1;//base.SetCol(2, e3);
+			break;
+		default:
+			mfem_error("The face_ind exceeds the number of faces in this dimension.");
+			break;
+	}
+}
+
+/** Maps the coordinate vectors of the first face to the coordinate vectors of the second face.
+*   nb_rot is the number of rotation to opperate so that the first node of each face match.
+*   The result map contains pairs of int, where the first int is the cofficient, and the 
+*   second int is the indice of the second face vector.
+*/
+// There shouldn't be any rotation in 2D.
+void GetLocalCoordMap2D(vector<pair<int,int> >& map, const int nb_rot)
+{
+	map.resize(2);
+	//First and second coordinate vectors should always be of opposite direction in 2D.
+	//TODO Maybe not
+	map[0] = pair<int,int>(-1,0);
+	map[1] = pair<int,int>(-1,1);
+}
+
+// Default parameter nb_rot=0 should be only use with a structured mesh.
+// Rotations follow the ordering of the nodes.
+/*void GetLocalCoordMap3D(vector<pair<int,int> >& map, const int nb_rot)
+{
+	map.resize(3);
+	// Normal to the face are always of opposite direction
+	map[2] = pair<int,int>(-1,2);
+	// nb_rot determines how local coordinates are oriented from one face to the other.
+	// See case 2 for an example.
+	switch(nb_rot)
+	{
+		case 0:
+			map[0] = pair<int,int>( 1,1);
+			map[1] = pair<int,int>( 1,0);
+			break;
+		case 1:
+			map[0] = pair<int,int>(-1,0);
+			map[1] = pair<int,int>( 1,1);
+			break;
+		case 2:
+			//first vector equals -1 times the second vector of the other face coordinates
+			map[0] = pair<int,int>(-1,1);
+			//second vector equals -1 times the first vector of the other face coordinates
+			map[1] = pair<int,int>(-1,0);
+			break;
+		case 3:
+			map[0] = pair<int,int>( 1,0);
+			map[1] = pair<int,int>(-1,1);
+			break;
+		default:
+			mfem_error("There shouldn't be that many rotations.");
+			break;
+	}
+}*/
+
+void GetLocalCoordMap3D(vector< pair<int,int> >& map, const int orientation)
+{
+	map.resize(3);
+	// orientation determines how local coordinates are oriented from one face to the other.
+	// See case 2 for an example.
+	switch(orientation)
+	{
+		case 0://{0, 1, 2, 3}
+			map[0] = pair<int,int>( 1,0);
+			map[1] = pair<int,int>( 1,1);
+			map[2] = pair<int,int>( 1,2);
+			break;
+		case 1://{0, 3, 2, 1}
+			map[0] = pair<int,int>( 1,1);
+			map[1] = pair<int,int>( 1,0);
+			map[2] = pair<int,int>(-1,2);
+			break;
+		case 2://{1, 2, 3, 0}
+			//first vector equals -1 times the second vector of the other face coordinates
+			map[0] = pair<int,int>(-1,1);
+			//second vector equals -1 times the first vector of the other face coordinates
+			map[1] = pair<int,int>( 1,0);
+			//third vector equals -1 times the third vector of the other face coordinates
+			map[2] = pair<int,int>( 1,2);
+			break;
+		case 3://{1, 0, 3, 2}
+			map[0] = pair<int,int>(-1,0);
+			map[1] = pair<int,int>( 1,1);
+			map[2] = pair<int,int>(-1,2);
+			break;
+		case 4://{2, 3, 0, 1}
+			map[0] = pair<int,int>(-1,0);
+			map[1] = pair<int,int>(-1,1);
+			map[2] = pair<int,int>( 1,2);
+			break;
+		case 5://{2, 1, 0, 3}
+			map[0] = pair<int,int>(-1,1);
+			map[1] = pair<int,int>(-1,0);
+			map[2] = pair<int,int>(-1,2);
+			break;
+		case 6://{3, 0, 1, 2}
+			map[0] = pair<int,int>( 1,1);
+			map[1] = pair<int,int>(-1,0);
+			map[2] = pair<int,int>( 1,2);
+			break;
+		case 7://{3, 2, 1, 0}
+			map[0] = pair<int,int>( 1,0);
+			map[1] = pair<int,int>(-1,1);
+			map[2] = pair<int,int>(-1,2);
+			break;
+		default:
+			mfem_error("There shouldn't be that many orientations.");
+			break;
+	}
+}
+
+/**
+*	Returns the change of matrix P from base_K2 to base_K1 according to the mapping map.
+*/
+void GetChangeOfBasis(const IntMatrix& base_K1, IntMatrix& base_K2,
+								const vector<pair<int,int> >& map, IntMatrix& P)
+{
+/*	int dim = map.size();
+	for (int j = 0; j < dim; j++)
+	{
+		int i = 0;
+		//we look if the vector is colinear with e_j
+		// Can be replaced by base_K2(j,i)!=0
+		while (base_K2(j,i)!=0) i++;
+		int coeff = map[i].first;
+		int ind = map[i].second;
+		for (int k = 0; k < dim; ++k)
+		{
+			P(k,j) = coeff * base_K1(k,ind);
+		}
+	}*/
+	//TODO make it valid for 3D!!!
+	int dim = base_K1.Height();
+	// for (int i = 0; i < dim; ++i)
+	// {
+	// 	int coeff = map[i].first;
+	// 	int ind   = map[i].second;
+	// 	for (int j = 0; j < dim; ++j)
+	// 	{
+	// 		int sum = 0;
+	// 		for (int k = 0; k < dim; ++k)
+	// 		{
+	// 			sum += coeff*base_K1(i,k)*base_K2(j,k);
+	// 		}
+	// 		P(ind,j) =  sum;
+	// 	}
+	// }
+	int i,j,ind;
+	double coeff;
+	for (int n = 0; n < dim; ++n)
+	{
+		i = 0;
+		while( base_K1(i,n)   == 0 ) ++i;
+		j = 0;
+		ind = map[n].second;
+		while( base_K2(j,ind) == 0 ) ++j;
+		coeff = map[n].first;
+		P(i,j) = coeff * base_K1(i,n) * base_K2(j,ind);
+	}
+}
+
+void GetChangeOfBasis2D(const int face_id1, const int face_id2, IntMatrix& P)
+{
+	// We add 8 because of C++ stupid definition of modulo
+	int nb_rot = (8 + face_id2 - face_id1 - 2)%4;
+	// if (face_id2!=-1)
+	// {
+		// cout << "face_id1=" << face_id1 << ", face_id2=" << face_id2 << ", nb_rot=" << nb_rot << endl;
+	// }
+	P.Zero();
+	switch(nb_rot)
+	{
+	case 0://Id=R^4
+		P(0,0) = 1;
+		P(1,1) = 1;
+		break;
+	case 1://R
+		P(1,0) = 1;
+		P(0,1) =-1;
+		break;
+	case 2://R²
+		P(0,0) =-1;
+		P(1,1) =-1;
+		break;
+	case 3://R³
+		P(1,0) =-1;
+		P(0,1) = 1;
+		break;
+	default:mfem_error("C++ modulo error in GetChangeOfBasis2D");
+	}
+}
+
+void GetChangeOfBasis(const int permutation, IntMatrix& P)
+{
+	int code1 = permutation/100;
+	int ind1  = code1/2;
+	int val1  = code1%2==0?-1:1;
+	int code2 = (permutation%100)/10;
+	int ind2  = code2/2;
+	int val2  = code2%2==0?-1:1;
+	int code3 = permutation%10;
+	int ind3  = code3/2;
+	int val3  = code3%2==0?-1:1;
+	P.Zero();
+	P(ind1,0) = val1;
+	P(ind2,1) = val2;
+	P(ind3,2) = val3;
+}
+
+/**
+*	Returns the face_id that identifies the face on the reference element, and nb_rot the
+*  "rotations" the face did between reference to physical spaces.
+*/
+void GetIdRotInfo(const int face_info, int& face_id, int& nb_rot){
+	int orientation = face_info % 64;
+	face_id = face_info / 64;
+	// Test if my understanding of mfem code is correct, error if not
+	//MFEM_ASSERT(orientation % 2 == 0, "Unexpected inside out face");
+	nb_rot = orientation;// / 2;
+}
+
+void GetFaceInfo(const Mesh* mesh, const int face, int& ind_elt1, int& ind_elt2, int& face_id1, int& face_id2, int& nb_rot1, int& nb_rot2)
+{
+	// We collect the indices of the two elements on the face, element1 is the master element,
+	// the one that defines the normal to the face.
+	mesh->GetFaceElements(face,&ind_elt1,&ind_elt2);
+	int info_elt1, info_elt2;
+	// We collect the informations on the face for the two elements.
+	mesh->GetFaceInfos(face,&info_elt1,&info_elt2);
+	GetIdRotInfo(info_elt1,face_id1,nb_rot1);//nb_rot1 is always 0 by convention
+	GetIdRotInfo(info_elt2,face_id2,nb_rot2);
+}
+
+/**
+*	Returns the permutation id, so that we can permute dofs to be in a structured case.
+*/
+int Permutation2D(const int face_id_trial, const int face_id_test)
+{
+	int perm = face_id_trial - face_id_test - 2;
+	perm = perm < 0 ? perm+4 : perm;
+	return perm;
+}
+
+/**
+*  Returns an integer that encrypts P.
+*/
+void Permutation3D(const int face_id1, const int face_id2, const int orientation, int& perm1, int& perm2)
+{
+	IntMatrix K1(3,3);
+	K1.Zero();
+	InitFaceCoord3D(face_id1, K1);
+	IntMatrix K2(3,3);
+	K2.Zero();
+	InitFaceCoord3D(face_id2, K2);
+	vector< pair<int,int> > map;
+	GetLocalCoordMap3D(map, orientation);
+	IntMatrix P(3,3);
+	P.Zero();
+	GetChangeOfBasis(K1, K2, map, P);
+	// cout << "orientation=" << orientation << endl;
+	// cout << P(0,0) << ", " << P(0,1) << ", " << P(0,2) << endl;
+	// cout << P(1,0) << ", " << P(1,1) << ", " << P(1,2) << endl;
+	// cout << P(2,0) << ", " << P(2,1) << ", " << P(2,2) << endl;
+	perm1  = 0;
+	// Encrypts first column
+	perm1 += 100*(0*(P(0,0)==-1) + 1*(P(0,0)==1) + 2*(P(1,0)==-1) + 3*(P(1,0)==1) + 4*(P(2,0)==-1) + 5*(P(2,0)==1));
+	// Encrypts second column
+	perm1 += 10 *(0*(P(0,1)==-1) + 1*(P(0,1)==1) + 2*(P(1,1)==-1) + 3*(P(1,1)==1) + 4*(P(2,1)==-1) + 5*(P(2,1)==1));
+	// Encrypts third column
+	perm1 +=     (0*(P(0,2)==-1) + 1*(P(0,2)==1) + 2*(P(1,2)==-1) + 3*(P(1,2)==1) + 4*(P(2,2)==-1) + 5*(P(2,2)==1));
+	// Encrypts the transposed permutation matrix in a second integer.
+	perm2  = 0;
+	perm2 += 100*(0*(P(0,0)==-1) + 1*(P(0,0)==1) + 2*(P(0,1)==-1) + 3*(P(0,1)==1) + 4*(P(0,2)==-1) + 5*(P(0,2)==1));
+	perm2 += 10 *(0*(P(1,0)==-1) + 1*(P(1,0)==1) + 2*(P(1,1)==-1) + 3*(P(1,1)==1) + 4*(P(1,2)==-1) + 5*(P(1,2)==1));
+	perm2 +=     (0*(P(2,0)==-1) + 1*(P(2,0)==1) + 2*(P(2,1)==-1) + 3*(P(2,1)==1) + 4*(P(2,2)==-1) + 5*(P(2,2)==1));
+}
+
+void GetPermutation(const int dim, const int face_id1, const int face_id2, const int orientation, int& perm1, int& perm2)
+{
+	switch(dim){
+		case 1:
+			mfem_error("Not yet implemented");
+			break;
+		case 2:
+			perm1 = Permutation2D(face_id1, face_id2);
+			perm2 = Permutation2D(face_id2, face_id1);
+			break;
+		case 3:
+			Permutation3D(face_id1, face_id2, orientation, perm1, perm2);
+			break;
+		default:
+			mfem_error("Dimension of the problem too high.");
+			break;
+	}
+}
+
+/**
+*	Hardcoded permutation due to arbitrary hardcoded orientation in geom.cpp.
+*   Will break if geom.cpp changes.
+*   This function could be improved by returning the 'permutation' parameters once,
+*   instead of recomputing them for every quadrature point...
+*/
+int GetFaceQuadIndex3D(const int face_id, const int orientation, const int qind, const int quads, Tensor<1,int>& ind_f)
+{
+	// cout << "orientation=" << orientation << endl;
+	int& k1 = ind_f(0);
+	int& k2 = ind_f(1);
+	int kf1,kf2;
+	kf1 = qind%quads;
+	kf2 = qind/quads;
+	switch(face_id)
+	{
+		case 0://BOTTOM
+			switch(orientation)
+			{
+				case 0://{0, 1, 2, 3}
+					k1 = kf1;
+					k2 = quads-1-kf2;
+					break;
+				case 1://{0, 3, 2, 1}
+					k1 = quads-1-kf2;
+					k2 = kf1;
+					break;
+				case 2://{1, 2, 3, 0}
+					k1 = quads-1-kf2;
+					k2 = quads-1-kf1;
+					break;
+				case 3://{1, 0, 3, 2}
+					k1 = quads-1-kf1;
+					k2 = quads-1-kf2;
+					break;
+				case 4://{2, 3, 0, 1}
+					k1 = quads-1-kf1;
+					k2 = kf2;
+					break;
+				case 5://{2, 1, 0, 3}
+					k1 = kf2;
+					k2 = quads-1-kf1;
+					break;
+				case 6://{3, 0, 1, 2}
+					k1 = kf2;
+					k2 = kf1;
+					break;
+				case 7://{3, 2, 1, 0}
+					k1 = kf1;
+					k2 = kf2;
+					break;
+				default:
+					mfem_error("This orientation does not exist in 3D");
+					break;
+			}
+			break;
+		case 1://SOUTH
+			switch(orientation)
+			{
+				case 0://{0, 1, 2, 3}
+					k1 = kf1;
+					k2 = kf2;
+					break;
+				case 1://{0, 3, 2, 1}
+					k1 = kf2;
+					k2 = kf1;
+					break;
+				case 2://{1, 2, 3, 0}
+					k1 = kf2;
+					k2 = quads-1-kf1;
+					break;
+				case 3://{1, 0, 3, 2}
+					k1 = quads-1-kf1;
+					k2 = kf2;
+					break;
+				case 4://{2, 3, 0, 1}
+					k1 = quads-1-kf1;
+					k2 = quads-1-kf2;
+					break;
+				case 5://{2, 1, 0, 3}
+					k1 = quads-1-kf2;
+					k2 = quads-1-kf1;
+					break;
+				case 6://{3, 0, 1, 2}
+					k1 = quads-1-kf2;
+					k2 = kf1;
+					break;
+				case 7://{3, 2, 1, 0}
+					k1 = kf1;
+					k2 = quads-1-kf2;
+					break;
+				default:
+					mfem_error("This orientation does not exist in 3D");
+					break;
+			}
+			break;
+		case 2://EAST
+			switch(orientation)
+			{
+				case 0://{0, 1, 2, 3}
+					k1 = kf1;
+					k2 = kf2;
+					break;
+				case 1://{0, 3, 2, 1}
+					k1 = kf2;
+					k2 = kf1;
+					break;
+				case 2://{1, 2, 3, 0}
+					k1 = kf2;
+					k2 = quads-1-kf1;
+					break;
+				case 3://{1, 0, 3, 2}
+					k1 = quads-1-kf1;
+					k2 = kf2;
+					break;
+				case 4://{2, 3, 0, 1}
+					k1 = quads-1-kf1;
+					k2 = quads-1-kf2;
+					break;
+				case 5://{2, 1, 0, 3}
+					k1 = quads-1-kf2;
+					k2 = quads-1-kf1;
+					break;
+				case 6://{3, 0, 1, 2}
+					k1 = quads-1-kf2;
+					k2 = kf1;
+					break;
+				case 7://{3, 2, 1, 0}
+					k1 = kf1;
+					k2 = quads-1-kf2;
+					break;
+				default:
+					mfem_error("This orientation does not exist in 3D");
+					break;
+			}
+			break;
+		case 3://NORTH
+			switch(orientation)
+			{
+				case 0://{0, 1, 2, 3}
+					k1 = quads-1-kf1;
+					k2 = kf2;
+					break;
+				case 1://{0, 3, 2, 1}
+					k1 = kf2;
+					k2 = quads-1-kf1;
+					break;
+				case 2://{1, 2, 3, 0}
+					k1 = kf2;
+					k2 = kf1;
+					break;
+				case 3://{1, 0, 3, 2}
+					k1 = kf1;
+					k2 = kf2;
+					break;
+				case 4://{2, 3, 0, 1}
+					k1 = kf1;
+					k2 = quads-1-kf2;
+					break;
+				case 5://{2, 1, 0, 3}
+					k1 = quads-1-kf2;
+					k2 = kf1;
+					break;
+				case 6://{3, 0, 1, 2}
+					k1 = quads-1-kf2;
+					k2 = quads-1-kf1;
+					break;
+				case 7://{3, 2, 1, 0}
+					k1 = quads-1-kf1;
+					k2 = quads-1-kf2;
+					break;
+				default:
+					mfem_error("This orientation does not exist in 3D");
+					break;
+			}
+			break;
+		case 4://WEST
+			switch(orientation)
+			{
+				case 0://{0, 1, 2, 3}
+					k1 = quads-1-kf1;
+					k2 = kf2;
+					break;
+				case 1://{0, 3, 2, 1}
+					k1 = kf2;
+					k2 = quads-1-kf1;
+					break;
+				case 2://{1, 2, 3, 0}
+					k1 = kf2;
+					k2 = kf1;
+					break;
+				case 3://{1, 0, 3, 2}
+					k1 = kf1;
+					k2 = kf2;
+					break;
+				case 4://{2, 3, 0, 1}
+					k1 = kf1;
+					k2 = quads-1-kf2;
+					break;
+				case 5://{2, 1, 0, 3}
+					k1 = quads-1-kf2;
+					k2 = kf1;
+					break;
+				case 6://{3, 0, 1, 2}
+					k1 = quads-1-kf2;
+					k2 = quads-1-kf1;
+					break;
+				case 7://{3, 2, 1, 0}
+					k1 = quads-1-kf1;
+					k2 = quads-1-kf2;
+					break;
+				default:
+					mfem_error("This orientation does not exist in 3D");
+					break;
+			}
+			break;
+		case 5://TOP
+			switch(orientation)
+			{
+				case 0://{0, 1, 2, 3}
+					k1 = kf1;
+					k2 = kf2;
+					break;
+				case 1://{0, 3, 2, 1}
+					k1 = kf2;
+					k2 = kf1;
+					break;
+				case 2://{1, 2, 3, 0}
+					k1 = kf2;
+					k2 = quads-1-kf1;
+					break;
+				case 3://{1, 0, 3, 2}
+					k1 = quads-1-kf1;
+					k2 = kf2;
+					break;
+				case 4://{2, 3, 0, 1}
+					k1 = quads-1-kf1;
+					k2 = quads-1-kf2;
+					break;
+				case 5://{2, 1, 0, 3}
+					k1 = quads-1-kf2;
+					k2 = quads-1-kf1;
+					break;
+				case 6://{3, 0, 1, 2}
+					k1 = quads-1-kf2;
+					k2 = kf1;
+					break;
+				case 7://{3, 2, 1, 0}
+					k1 = kf1;
+					k2 = quads-1-kf2;
+					break;
+				default:
+					mfem_error("This orientation does not exist in 3D");
+					break;
+			}
+			break;
+		default:
+			mfem_error("This face_id does not exist in 3D");
+			break;	
+	}
+	return k1 + quads*k2;
+}
+
+int GetFaceQuadIndex(const int dim, const int face_id, const int orientation, const int qind, const int quads, Tensor<1,int>& ind_f)
+{
+	int res = 0;
+	switch(dim)
+	{
+		case 1:
+			break;
+		case 2:
+            if(face_id<=1){//SOUTH or EAST (canonical ordering)
+               res = ind_f(0) = qind;
+            }else{//NORTH or WEST (counter-canonical ordering)
+               res = ind_f(0) = quads-1-qind;
+            }
+            break;
+        case 3:
+			res = GetFaceQuadIndex3D(face_id, orientation, qind, quads, ind_f);
+			break;
+		default:
+			mfem_error("Dimension too high.");
+			break;
+	}
+	return res;
+}
+
+const int GetGlobalQuadIndex(const int dim, const int face_id, const int quads, Tensor<1,int>& ind_f)
+{
+   switch(dim)
+   {
+      case 1:
+         if (face_id==0)//WEST
+         {
+            return 0;
+         }else{//EAST
+            return quads-1;
+         }
+      case 2:
+            switch(face_id)
+            {
+               case 0://SOUTH
+                  return ind_f(0);
+               case 1://EAST
+                  return quads-1 + ind_f(0)*quads;
+               case 2://NORTH
+                  return ind_f(0) + (quads-1)*quads;
+               case 3://WEST
+                  return ind_f(0)*quads;
+            }
+        case 3:
+         switch(face_id)
+         {
+            case 0://BOTTOM
+               return ind_f(0) + ind_f(1)*quads;
+            case 1://SOUTH
+               return ind_f(0) + ind_f(1)*quads*quads;
+            case 2://EAST
+               return (quads-1) + ind_f(0)*quads + ind_f(1)*quads*quads;
+            case 3://NORTH
+               return ind_f(0) + (quads-1)*quads + ind_f(1)*quads*quads;
+            case 4://WEST
+               return ind_f(0)*quads + ind_f(1)*quads*quads;
+            case 5://TOP
+               return ind_f(0) + ind_f(1)*quads + (quads-1)*quads*quads;
+         }
+      default:
+         mfem_error("Dimension too high.");
+         break;
+   }
+   return -1;
+}
+
+}

fem/dgfacefunctions.hpp

@@ -0,0 +1,119 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+//This file contains useful functions to compute fluxes for DG methods.
+
+
+#ifndef MFEM_DGFACEFUNC
+#define MFEM_DGFACEFUNC
+#include "dalg.hpp"
+
+using std::vector;
+using std::pair;
+
+namespace mfem
+{
+
+/**
+*	Returns the canonical coordinate vectors e_1 and e_2.
+*/
+void getBaseVector2D(Vector& e1, Vector& e2);
+
+/**
+*	Returns the canonical coordinate vectors e_1, e_2 and e_3.
+*/
+void getBaseVector3D(Vector& e1, Vector& e2, Vector& e3);
+
+/** A function that initialize the local coordinate base for a face with
+*   indice face_ind.
+*   This returns the local face coordinate base expressed in reference
+*   element coordinate.
+*/
+// Highly dependent of the node ordering from geom.cpp
+void InitFaceCoord2D(const int face_id, IntMatrix& base);
+
+// Highly dependent of the node ordering from geom.cpp
+void InitFaceCoord3D(const int face_id, IntMatrix& base);
+
+/** Maps the coordinate vectors of the first face to the coordinate vectors of the second face.
+*   nb_rot is the number of rotation to opperate so that the first node of each face match.
+*   The result map contains pairs of int, where the first int is a direction cofficient,
+*	 and the second int is the indice of the second face vector.
+*/
+// There shouldn't be any rotation in 2D.
+void GetLocalCoordMap2D(vector<pair<int,int> >& map, const int nb_rot = 0);
+
+// Rotations follow the ordering of the nodes.
+void GetLocalCoordMap3D(vector<pair<int,int> >& map, const int nb_rot);
+
+/**
+*	Returns the change of matrix P from base_K2 to base_K1 according to the mapping map.
+*/
+void GetChangeOfBasis(const IntMatrix& base_K1, IntMatrix& base_K2,
+								const vector<pair<int,int> >& map, IntMatrix& P);
+
+void GetChangeOfBasis(const int permutation, IntMatrix& P);
+
+/**
+*	Returns the change of coordinate from second element to first element on a 2D face.
+*/
+void GetChangeOfBasis2D(const int face_id1, const int face_id2, IntMatrix& P);
+
+/**
+*	Returns the indices, face ID, and number of rotations, of the two element sharing a face.
+*  The number of rotations is relative to the element 1, so nb_rot1 is always 0.
+*/
+void GetFaceInfo(const Mesh* mesh, const int face,
+						int& ind_elt1, int& ind_elt2,
+						int& face_id1, int& face_id2,
+						int& nb_rot1, int& nb_rot2);
+
+/**
+*	Returns the face_id that identifies the face on the reference element, and nb_rot the
+*  "rotations" the face did between reference to physical spaces.
+*/
+void GetIdRotInfo(const int face_info, int& face_id, int& nb_rot);
+
+/**
+*	Returns an integer identifying the permutation to apply to be in structured-
+*  like configuration for 2D hex meshes.
+*/
+int Permutation2D(const int face_id_trial, const int face_id_test);
+
+/**
+*	Returns an integer identifying the permutation to apply to be in structured-
+*  like configuration for 3D hex meshes.
+*/
+void Permutation3D(const int face_id1, const int face_id2, const int orientation, int& perm1, int& perm2);
+
+/**
+*	Returns an integer identifying the permutation to apply to be in structured-
+*  like configuration.
+*/
+void GetPermutation(const int dim, const int face_id1, const int face_id2, const int orientation, int& perm1, int& perm2);
+
+int GetFaceQuadIndex3D(const int face_id, const int orientation, const int qind, const int quads, Tensor<1,int>& ind_f);
+
+/**
+*	Returns the indices of a quadrature point on the face of an hex element relative to the index of the quadrature
+*  point on the reference face.
+*/
+int GetFaceQuadIndex(const int dim, const int face_id, const int orientation, const int qind, const int quads, Tensor<1,int>& ind_f);
+
+/**
+*	Returns the indices of a quadrature point on the element relative to the index of the quadrature
+*  point on the reference face.
+*/
+const int GetGlobalQuadIndex(const int dim, const int face_id, const int quads, Tensor<1,int>& ind_f);
+
+}
+
+#endif // MFEM_DGFACEFUNC

fem/dgpabilininteg.hpp

@@ -0,0 +1,174 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+
+// This file contains a prototype version for Discontinuous Galerkin Partial assembly
+
+#ifndef MFEM_DGPABILININTEG
+#define MFEM_DGPABILININTEG
+
+#include "../config/config.hpp"
+#include "bilininteg.hpp"
+#include "dalg.hpp"
+
+#include "fem.hpp"
+#include <cmath>
+#include <algorithm>
+#include "../linalg/vector.hpp"
+
+namespace mfem
+{
+
+/**
+*  The different operators available for the Kernels
+*/
+enum PAOp { BtDB, BtDG, GtDB, GtDG };
+
+
+/**
+*	A class that describes the Convection Equation using DG for Partial Assembly.
+*/
+class DGConvectionEquation
+{
+public:
+   /**
+   *  Defines the Kernel to apply to the Domain
+   */
+   static const PAOp OpName = BtDG;
+
+   /**
+   *  Defines the variables needed to build D for the Domain kernel
+   */
+   struct Args {
+      Args(VectorCoefficient& _q, double _a = 1.0, double _b = -1.0) : q(_q), a(_a), b(_b) {}
+      VectorCoefficient& q;
+      double a;
+      double b;
+   };
+
+   /**
+   *  Returns the values of the D tensor at a given integration Point.
+   */
+   void evalD(Tensor<1>& res, ElementTransformation *Tr, const IntegrationPoint& ip,
+                  const Args& args)
+   {
+      const int dim = res.size(0);
+      Vector qvec(dim);
+      const DenseMatrix& locD = Tr->AdjugateJacobian();
+      args.q.Eval(qvec, *Tr, ip);
+      for (int i = 0; i < dim; ++i)
+      {
+         double val = 0.0;
+         for (int j = 0; j < dim; ++j)
+         {
+            val += locD(i,j) * qvec(j);
+         }
+         res(i) = ip.weight * args.a * val;
+      }
+   }
+
+   /**
+   *  Returns the values of the D tensor at a given integration Point.
+   */
+   void evalD(Tensor<1>& res, ElementTransformation *Tr, const IntegrationPoint& ip,
+                  const Tensor<2>& Jac, const Args& args)
+   {
+      const int dim = res.size(0);
+      Vector qvec(dim);
+      args.q.Eval(qvec, *Tr, ip);
+      Tensor<2> Adj(dim,dim);
+      adjugate(Jac,Adj);
+      for (int i = 0; i < dim; ++i)
+      {
+         double val = 0.0;
+         for (int j = 0; j < dim; ++j)
+         {
+            val += Adj(i,j) * qvec(j);
+         }
+         res(i) = ip.weight * args.a * val;
+      }
+   }
+
+   /**
+   *  Defines the Kernel to apply to the Faces
+   */
+   static const PAOp FaceOpName = BtDB;
+
+   /**
+   *  Returns the values of the Dint and Dext tensors at a given integration Point for
+   *  each element over a face.
+   */
+   void evalFaceD(double& res11, double& res21, double& res22, double& res12,
+      const FaceElementTransformations* face_tr, const Vector& normal,
+      const IntegrationPoint& ip1, const IntegrationPoint& ip2,
+      const Args& args)
+   {
+      const int dim = normal.Size();
+      Vector qvec(dim);
+      // FIXME: qvec might be discontinuous if not constant with a periodic mesh
+      // We should then use the evaluation on Elem2 and eip2
+      args.q.Eval( qvec, *(face_tr->Elem1), ip1 );
+      const double res = qvec * normal;
+      const double a = -args.a, b = args.b;
+      res11 = ip1.weight * (   a/2 * res + b * abs(res) );
+      res21 = ip1.weight * (   a/2 * res - b * abs(res) );
+      res22 = ip1.weight * ( - a/2 * res + b * abs(res) );
+      res12 = ip1.weight * ( - a/2 * res - b * abs(res) );
+   }
+
+   void evalFaceD(double& res11, double& res21, double& res22, double& res12,
+      const FaceElementTransformations* face_tr, const Vector& normal,
+      const IntegrationPoint& ip1, const IntegrationPoint& ip2,
+      const Tensor<2>& Jac1, const Tensor<2>& Jac2,
+      const Args& args)
+   {
+      const int dim = normal.Size();
+      Vector qvec(dim);
+      // FIXME: qvec might be discontinuous if not constant with a periodic mesh
+      // We should then use the evaluation on Elem2 and eip2
+      args.q.Eval( qvec, *(face_tr->Elem1), ip1 );
+      const double res = qvec * normal;
+      const double a = -args.a, b = args.b;
+      res11 = ip1.weight * (   a/2 * res + b * abs(res) );
+      res21 = ip1.weight * (   a/2 * res - b * abs(res) );
+      res22 = ip1.weight * ( - a/2 * res + b * abs(res) );
+      res12 = ip1.weight * ( - a/2 * res - b * abs(res) );
+   }  
+};
+
+class MassEquation
+{
+public:
+   static const PAOp OpName = BtDB;
+
+   struct ArgsEmpty{};
+
+   void evalD(double& res, ElementTransformation* Tr, const IntegrationPoint& ip,
+               const Tensor<2>& Jac, ArgsEmpty args = {})
+   {
+      res = ip.weight * det(Jac);
+   }
+
+   struct ArgsCoeff
+   {
+      Coefficient& coeff;
+   };
+
+   void evalD(double& res, ElementTransformation* Tr, const IntegrationPoint& ip,
+               const Tensor<2>& Jac, ArgsCoeff& args)
+   {
+      res = args.coeff.Eval(*Tr, ip) * ip.weight * det(Jac);
+   }
+};
+
+}
+
+#endif //MFEM_DGPABILININTEG

fem/fem.hpp

@@ -26,6 +26,8 @@
 #include "linearform.hpp"
 #include "nonlinearform.hpp"
 #include "bilinearform.hpp"
+#include "bilinearformoper.hpp"
+#include "painteg.hpp"
 #include "hybridization.hpp"
 #include "datacollection.hpp"
 #include "estimators.hpp"

fem/fespace.hpp

@@ -249,6 +249,16 @@ public:
    /// Returns number of degrees of freedom.
    inline int GetNDofs() const { return ndofs; }
 
+   /// Returns number of degrees of freedom in each direction.
+   inline const int GetNDofs1d() const { return GetFE(0)->GetOrder()+1; }
+
+   /// Returns number of quadrature points in each direction.
+   inline const int GetNQuads1d(const int order) const
+   { 
+    const IntegrationRule &ir1d = IntRules.Get(Geometry::SEGMENT, order);
+    return ir1d.GetNPoints();
+  }
+
    /// Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
    inline int GetVSize() const { return vdim * ndofs; }
 
@@ -301,9 +311,9 @@ public:
    ElementTransformation *GetElementTransformation(int i) const
    { return mesh->GetElementTransformation(i); }
 
-   /** @brief Returns the transformation defining the @a i-th element in the
-       user-defined variable @a ElTr. */
-   void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
+   /** Returns the transformation defining the i-th element in the user-defined
+       variable. */
+   void GetElementTransformation(int i, IsoparametricTransformation *ElTr) const
    { mesh->GetElementTransformation(i, ElTr); }
 
    /// Returns ElementTransformation for the @a i-th boundary element.

fem/nonlininteg.hpp

@@ -13,6 +13,7 @@
 #define MFEM_NONLININTEG
 
 #include "../config/config.hpp"
+#include "fespace.hpp"
 #include "fe.hpp"
 #include "coefficient.hpp"
 
@@ -71,6 +72,67 @@ public:
    virtual ~NonlinearFormIntegrator() { }
 };
 
+
+class Integrator
+{
+protected:
+   const IntegrationRule *IntRule;
+
+public:
+   Integrator(const IntegrationRule *_IntRule = NULL) :
+      IntRule(_IntRule) { }
+
+   void SetIntegrationRule(const IntegrationRule *ir) { IntRule = ir; }
+};
+
+
+class LinearFESpaceIntegrator : public Integrator
+{
+public:
+   LinearFESpaceIntegrator(const IntegrationRule *_IntRule = NULL) :
+      Integrator(_IntRule) { }
+
+   virtual ~LinearFESpaceIntegrator() { }
+
+   /// Internally assemble the integrator for the specific trial and
+   /// test spaces (with an optional vector u for semilinear forms).
+   virtual void Assemble(FiniteElementSpace *trial_fes,
+                         FiniteElementSpace *test_fes) { }
+
+   /// Apply the action A * x = y.
+   virtual void AddMult(const Vector &x, Vector &y)
+   { mfem_error("Not supported"); }
+
+   /// Apply the transposed action A^T * x = y.
+   virtual void AddMultTranspose(const Vector &x, Vector &y)
+   { mfem_error("Not supported"); }
+};
+
+
+class NonlinearFESpaceIntegrator : public Integrator
+{
+public:
+   NonlinearFESpaceIntegrator(const IntegrationRule *_IntRule = NULL) :
+      Integrator(_IntRule) { }
+
+   virtual ~NonlinearFESpaceIntegrator() { }
+
+   /// Internally assemble the integrator for the specific trial and
+   /// test spaces (with an optional vector u for semilinear forms).
+   virtual void Assemble(FiniteElementSpace *trial_fes,
+                         FiniteElementSpace *test_fes,
+                         const Vector &u) { }
+
+   /// Apply the action A(u) * x = y.
+   virtual void AddMult(const Vector &x, Vector &y)
+   { mfem_error("Not supported"); }
+
+   /// Apply the transposed action A(u)^T * x = y.
+   virtual void AddMultTranspose(const Vector &x, Vector &y)
+   { mfem_error("Not supported"); }
+};
+
+
 /** The abstract base class BlockNonlinearFormIntegrator is
     a generalization of the NonlinearFormIntegrator class suitable
     for block state vectors. */

fem/painteg.cpp

@@ -0,0 +1,704 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+// Implementation of FESpaceIntegrators.
+
+#include "fem.hpp"
+
+namespace mfem
+{
+
+static void ComputeBasis1d(const FiniteElement *fe,
+                           const TensorBasisElement *tfe, int ir_order,
+                           DenseMatrix &shape1d)
+{
+   // Compute the 1d shape functions and gradients
+   const Poly_1D::Basis &basis1d = tfe->GetBasis1D();
+   const IntegrationRule &ir1d = IntRules.Get(Geometry::SEGMENT, ir_order);
+
+   const int quads1d = ir1d.GetNPoints();
+   const int dofs = fe->GetOrder() + 1;
+
+   shape1d.SetSize(dofs, quads1d);
+
+   Vector u(dofs);
+   for (int k = 0; k < quads1d; k++)
+   {
+      const IntegrationPoint &ip = ir1d.IntPoint(k);
+      basis1d.Eval(ip.x, u);
+      for (int i = 0; i < dofs; i++)
+      {
+         shape1d(i, k) = u(i);
+      }
+   }
+}
+
+static void ComputeBasis1d(const FiniteElement *fe,
+                           const TensorBasisElement *tfe, int ir_order,
+                           DenseMatrix &shape1d, DenseMatrix &dshape1d)
+{
+   // Compute the 1d shape functions and gradients
+   const Poly_1D::Basis &basis1d = tfe->GetBasis1D();
+   const IntegrationRule &ir1d = IntRules.Get(Geometry::SEGMENT, ir_order);
+
+   const int quads1d = ir1d.GetNPoints();
+   const int dofs = fe->GetOrder() + 1;
+
+   shape1d.SetSize(dofs, quads1d);
+   dshape1d.SetSize(dofs, quads1d);
+
+   Vector u(dofs);
+   Vector d(dofs);
+   for (int k = 0; k < quads1d; k++)
+   {
+      const IntegrationPoint &ip = ir1d.IntPoint(k);
+      basis1d.Eval(ip.x, u, d);
+      for (int i = 0; i < dofs; i++)
+      {
+         shape1d(i, k) = u(i);
+         dshape1d(i, k) = d(i);
+      }
+   }
+}
+
+void PADiffusionIntegrator::Assemble(FiniteElementSpace *_trial_fes,
+                                     FiniteElementSpace *_test_fes)
+{
+   // Assumption: trial and test fespaces are the same (no mixed forms yet)
+   fes = _trial_fes;
+
+   // Assumption: all are same finite elements
+   const FiniteElement *fe = fes->GetFE(0);
+
+   // Get the corresponding tensor basis element
+   const TensorBasisElement *tfe = dynamic_cast<const TensorBasisElement*>(fe);
+
+   // Set integration rule
+   int ir_order;
+   if (!IntRule)
+   {
+      const int dim = fe->GetDim();
+      if (fe->Space() == FunctionSpace::Pk)
+      {
+         ir_order = 2*fe->GetOrder() - 2;
+      }
+      else
+         // order = 2*fe.GetOrder() - 2;  // <-- this seems to work fine too
+      {
+         ir_order = 2*fe->GetOrder() + dim - 1;
+      }
+
+      if (fe->Space() == FunctionSpace::rQk)
+      {
+         SetIntegrationRule(&RefinedIntRules.Get(fe->GetGeomType(), ir_order));
+      }
+      else
+      {
+         SetIntegrationRule(&IntRules.Get(fe->GetGeomType(), ir_order));
+      }
+   }
+   else
+   {
+      ir_order = IntRule->GetOrder();
+   }
+
+   // Store the 1d shape functions and gradients
+   ComputeBasis1d(fes->GetFE(0), tfe, ir_order, shape1d, dshape1d);
+
+   // Create the operator
+   const int elems   = fes->GetNE();
+   const int dim     = fe->GetDim();
+   const int quads   = IntRule->GetNPoints();
+   const int entries = dim * (dim + 1) / 2;
+   Dtensor.SetSize(entries, quads, elems);
+
+   DenseMatrix invdfdx(dim, dim);
+   DenseMatrix mat(dim, dim);
+   DenseMatrix cmat(dim, dim);
+
+   Coefficient *coeff = integ->Q;
+   MatrixCoefficient *mcoeff = integ->MQ;
+   for (int e = 0; e < fes->GetNE(); e++)
+   {
+      ElementTransformation *Tr = fes->GetElementTransformation(e);
+      DenseMatrix &Dmat = Dtensor(e);
+      for (int k = 0; k < quads; k++)
+      {
+         const IntegrationPoint &ip = IntRule->IntPoint(k);
+         Tr->SetIntPoint(&ip);
+         const DenseMatrix &temp = Tr->AdjugateJacobian();
+         MultABt(temp, temp, mat);
+         mat *= ip.weight / Tr->Weight();
+
+         if (coeff != NULL)
+         {
+            const double c = coeff->Eval(*Tr, ip);
+            for (int j = 0, l = 0; j < dim; j++)
+               for (int i = j; i < dim; i++, l++)
+               {
+                  Dmat(l, k) = c * mat(i, j);
+               }
+
+         }
+         else if (mcoeff != NULL)
+         {
+            mcoeff->Eval(cmat, *Tr, ip);
+            for (int j = 0, l = 0; j < dim; j++)
+               for (int i = j; i < dim; i++, l++)
+               {
+                  Dmat(l, k) = cmat(i, j) * mat(i, j);
+               }
+
+         }
+         else
+         {
+            for (int j = 0, l = 0; j < dim; j++)
+               for (int i = j; i < dim; i++, l++)
+               {
+                  Dmat(l, k) = mat(i, j);
+               }
+         }
+
+      }
+   }
+}
+
+
+void PADiffusionIntegrator::MultSeg(const Vector &V, Vector &U)
+{
+   const int dofs1d = shape1d.Height();
+   const int quads1d = shape1d.Width();
+   const int quads = quads1d;
+   const int vdim = fes->GetVDim();
+
+   Vector Q(quads1d);
+
+   int offset = 0;
+   for (int e = 0; e < fes->GetNE(); ++e)
+   {
+      for (int vd = 0; vd < vdim; ++vd)
+      {
+         const Vector Vmat(V.GetData() + offset, dofs1d);
+         Vector Umat(U.GetData() + offset, dofs1d);
+
+         // Q_k1 = dshape_j1_k1 * V_i1
+         dshape1d.MultTranspose(Vmat, Q);
+
+         double *data_q = Q.GetData();
+         const double *data_d = Dtensor(e).GetData();
+         for (int k = 0; k < quads; ++k)
+         {
+            data_q[k] *= data_d[k];
+         }
+
+         // Q_k1 = dshape_j1_k1 * Q_k1
+         dshape1d.AddMult(Q, Umat);
+
+         // increment offset into E-vectors.
+         offset += dofs1d;
+      }
+   }
+}
+
+
+void PADiffusionIntegrator::MultQuad(const Vector &V, Vector &U)
+{
+   const int dim = 2;
+   const int terms = dim*(dim+1)/2;
+   const int vdim = fes->GetVDim();
+
+   const int dofs1d = shape1d.Height();
+   const int quads1d = shape1d.Width();
+
+   const int dofs   = dofs1d * dofs1d;
+   const int quads  = IntRule->GetNPoints();
+
+   DenseTensor QQ(quads1d, quads1d, dim);
+   DenseMatrix DQ(dofs1d, quads1d);
+
+   int offset = 0;
+   for (int e = 0; e < fes->GetNE(); ++e)
+   {
+      for (int vd = 0; vd < vdim; ++vd)
+      {
+         const DenseMatrix Vmat(V.GetData() + offset, dofs1d, dofs1d);
+         DenseMatrix Umat(U.GetData() + offset, dofs1d, dofs1d);
+
+         // DQ_j2_k1   = E_j1_j2  * dshape_j1_k1 -- contract in x direction
+         // QQ_0_k1_k2 = DQ_j2_k1 * shape_j2_k2  -- contract in y direction
+         MultAtB(Vmat, dshape1d, DQ);
+         MultAtB(DQ, shape1d, QQ(0));
+
+         // DQ_j2_k1   = E_j1_j2  * shape_j1_k1  -- contract in x direction
+         // QQ_1_k1_k2 = DQ_j2_k1 * dshape_j2_k2 -- contract in y direction
+         MultAtB(Vmat, shape1d, DQ);
+         MultAtB(DQ, dshape1d, QQ(1));
+
+         // QQ_c_k1_k2 = Dmat_c_d_k1_k2 * QQ_d_k1_k2
+         // NOTE: (k1, k2) = k -- 1d index over tensor product of quad points
+         double *data_qq = QQ(0).GetData();
+         const double *data_d = Dtensor(e).GetData();
+         for (int k = 0; k < quads; ++k)
+         {
+            const double D00 = data_d[terms*k + 0];
+            const double D01 = data_d[terms*k + 1];
+            const double D11 = data_d[terms*k + 2];
+
+            const double q0 = data_qq[0*quads + k];
+            const double q1 = data_qq[1*quads + k];
+
+            data_qq[0*quads + k] = D00 * q0 + D01 * q1;
+            data_qq[1*quads + k] = D01 * q0 + D11 * q1;
+         }
+
+         // DQ_i2_k1   = shape_i2_k2  * QQ_0_k1_k2
+         // U_i1_i2   += dshape_i1_k1 * DQ_i2_k1
+         MultABt(shape1d, QQ(0), DQ);
+         AddMultABt(dshape1d, DQ, Umat);
+
+         // DQ_i2_k1   = dshape_i2_k2 * QQ_1_k1_k2
+         // U_i1_i2   += shape_i1_k1  * DQ_i2_k1
+         MultABt(dshape1d, QQ(1), DQ);
+         AddMultABt(shape1d, DQ, Umat);
+
+         // increment offset
+         offset += dofs;
+      }
+   }
+}
+
+void PADiffusionIntegrator::MultHex(const Vector &V, Vector &U)
+{
+   const int dim = 3;
+   const int terms = dim*(dim+1)/2;
+   const int vdim = fes->GetVDim();
+
+   const int dofs1d = shape1d.Height();
+   const int quads1d = shape1d.Width();
+
+   const int dofs   = dofs1d * dofs1d * dofs1d;
+   const int quads  = IntRule->GetNPoints();
+
+   DenseMatrix Q(quads1d, dim);
+   DenseTensor QQ(quads1d, quads1d, dim);
+
+   Array<double> QQQmem(quads1d * quads1d * quads1d * dim);
+   double *data_qqq = QQQmem.GetData();
+   DenseTensor QQQ0(data_qqq + 0*quads, quads1d, quads1d, quads1d);
+   DenseTensor QQQ1(data_qqq + 1*quads, quads1d, quads1d, quads1d);
+   DenseTensor QQQ2(data_qqq + 2*quads, quads1d, quads1d, quads1d);
+
+   int offset = 0;
+   for (int e = 0; e < fes->GetNE(); ++e)
+   {
+      for (int vd = 0; vd < vdim; ++vd)
+      {
+         const DenseTensor Vmat(V.GetData() + offset, dofs1d, dofs1d, dofs1d);
+         DenseTensor Umat(U.GetData() + offset, dofs1d, dofs1d, dofs1d);
+
+         // QQQ_0_k1_k2_k3 = dshape_j1_k1 * shape_j2_k2  * shape_j3_k3  * Vmat_j1_j2_j3
+         // QQQ_1_k1_k2_k3 = shape_j1_k1  * dshape_j2_k2 * shape_j3_k3  * Vmat_j1_j2_j3
+         // QQQ_2_k1_k2_k3 = shape_j1_k1  * shape_j2_k2  * dshape_j3_k3 * Vmat_j1_j2_j3
+         QQQ0 = 0.; QQQ1 = 0.; QQQ2 = 0.;
+         for (int j3 = 0; j3 < dofs1d; ++j3)
+         {
+            QQ = 0.;
+            for (int j2 = 0; j2 < dofs1d; ++j2)
+            {
+               Q = 0.;
+               for (int j1 = 0; j1 < dofs1d; ++j1)
+               {
+                  for (int k1 = 0; k1 < quads1d; ++k1)
+                  {
+                     Q(k1, 0) += Vmat(j1, j2, j3) * dshape1d(j1, k1);
+                     Q(k1, 1) += Vmat(j1, j2, j3) * shape1d(j1, k1);
+                  }
+               }
+               for (int k2 = 0; k2 < quads1d; ++k2)
+                  for (int k1 = 0; k1 < quads1d; ++k1)
+                  {
+                     QQ(k1, k2, 0) += Q(k1, 0) * shape1d(j2, k2);
+                     QQ(k1, k2, 1) += Q(k1, 1) * dshape1d(j2, k2);
+                     QQ(k1, k2, 2) += Q(k1, 1) * shape1d(j2, k2);
+                  }
+            }
+            for (int k3 = 0; k3 < quads1d; ++k3)
+               for (int k2 = 0; k2 < quads1d; ++k2)
+                  for (int k1 = 0; k1 < quads1d; ++k1)
+                  {
+                     QQQ0(k1, k2, k3) += QQ(k1, k2, 0) * shape1d(j3, k3);
+                     QQQ1(k1, k2, k3) += QQ(k1, k2, 1) * shape1d(j3, k3);
+                     QQQ2(k1, k2, k3) += QQ(k1, k2, 2) * dshape1d(j3, k3);
+                  }
+         }
+
+         // QQQ_c_k1_k2_k3 = Dmat_c_d_k1_k2_k3 * QQQ_d_k1_k2_k3
+         // NOTE: (k1, k2, k3) = q -- 1d quad point index
+         const double *data_d = Dtensor(e).GetData();
+         for (int k = 0; k < quads; ++k)
+         {
+            const double D00 = data_d[terms*k + 0];
+            const double D01 = data_d[terms*k + 1];
+            const double D02 = data_d[terms*k + 2];
+            const double D11 = data_d[terms*k + 3];
+            const double D12 = data_d[terms*k + 4];
+            const double D22 = data_d[terms*k + 5];
+
+            const double q0 = data_qqq[0*quads + k];
+            const double q1 = data_qqq[1*quads + k];
+            const double q2 = data_qqq[2*quads + k];
+
+            data_qqq[0*quads + k] = D00 * q0 + D01 * q1 + D02 * q2;
+            data_qqq[1*quads + k] = D01 * q0 + D11 * q1 + D12 * q2;
+            data_qqq[2*quads + k] = D02 * q0 + D12 * q1 + D22 * q2;
+         }
+
+         // Apply transpose of the first operator that takes V -> QQQd -- QQQd -> U
+         for (int k3 = 0; k3 < quads1d; ++k3)
+         {
+            QQ = 0.;
+            for (int k2 = 0; k2 < quads1d; ++k2)
+            {
+               Q = 0.;
+               for (int k1 = 0; k1 < quads1d; ++k1)
+               {
+                  for (int i1 = 0; i1 < dofs1d; ++i1)
+                  {
+                     Q(i1, 0) += QQQ0(k1, k2, k3) * dshape1d(i1, k1);
+                     Q(i1, 1) += QQQ1(k1, k2, k3) * shape1d(i1, k1);
+                     Q(i1, 2) += QQQ2(k1, k2, k3) * shape1d(i1, k1);
+                  }
+               }
+               for (int i2 = 0; i2 < dofs1d; ++i2)
+                  for (int i1 = 0; i1 < dofs1d; ++i1)
+                  {
+                     QQ(i1, i2, 0) += Q(i1, 0) * shape1d(i2, k2);
+                     QQ(i1, i2, 1) += Q(i1, 1) * dshape1d(i2, k2);
+                     QQ(i1, i2, 2) += Q(i1, 2) * shape1d(i2, k2);
+                  }
+            }
+            for (int i3 = 0; i3 < dofs1d; ++i3)
+               for (int i2 = 0; i2 < dofs1d; ++i2)
+                  for (int i1 = 0; i1 < dofs1d; ++i1)
+                  {
+                     Umat(i1, i2, i3) +=
+                        QQ(i1, i2, 0) * shape1d(i3, k3) +
+                        QQ(i1, i2, 1) * shape1d(i3, k3) +
+                        QQ(i1, i2, 2) * dshape1d(i3, k3);
+                  }
+         }
+
+         // increment offset
+         offset += dofs;
+      }
+   }
+}
+
+void PADiffusionIntegrator::AddMult(const Vector &x, Vector &y)
+{
+   const int dim = fes->GetMesh()->Dimension();
+
+   switch (dim)
+   {
+   case 1: MultSeg(x, y); break;
+   case 2: MultQuad(x, y); break;
+   case 3: MultHex(x, y); break;
+   default: mfem_error("Not yet supported"); break;
+   }
+}
+
+
+void PAMassIntegrator::Assemble(FiniteElementSpace *_trial_fes,
+                                FiniteElementSpace *_test_fes)
+{
+   // Assumption: trial and test fespaces are the same (no mixed forms yet)
+   fes = _trial_fes;
+
+   // Assumption: all are same finite elements
+   const FiniteElement *fe = fes->GetFE(0);
+
+   // Get the corresponding tensor basis element
+   const TensorBasisElement *tfe = dynamic_cast<const TensorBasisElement*>(fe);
+
+   // Set integration rule
+   int ir_order;
+   if (!IntRule)
+   {
+      // int order = 2 * el.GetOrder();
+      // ir_order = 2 * fe.GetOrder() + Trans.OrderW();
+      ir_order = 2 * fe->GetOrder() + 1;
+
+      if (fe->Space() == FunctionSpace::rQk)
+      {
+         SetIntegrationRule(&RefinedIntRules.Get(fe->GetGeomType(), ir_order));
+      }
+      else
+      {
+         SetIntegrationRule(&IntRules.Get(fe->GetGeomType(), ir_order));
+      }
+   }
+   else
+   {
+      ir_order = IntRule->GetOrder();
+   }
+
+   ComputeBasis1d(fes->GetFE(0), tfe, ir_order, shape1d);
+
+   // Create the operator
+   const int nelem   = fes->GetNE();
+   const int dim     = fe->GetDim();
+   const int quads   = IntRule->GetNPoints();
+   const int vdim    = integ ? 1 : dim;
+   Dtensor.SetSize(quads, vdim, nelem);
+
+   Coefficient *coeff = NULL;
+   VectorCoefficient *vcoeff = NULL;
+   if (integ)
+   {
+      coeff = integ->Q;
+   }
+   else if (vinteg)
+   {
+      coeff = vinteg->Q;
+      vcoeff = vinteg->VQ;
+      if (vinteg->MQ != NULL) mfem_error("Not supported.");
+   }
+   DenseMatrix invdfdx(dim, dim);
+   DenseMatrix mat(dim, dim);
+   Vector cv(vdim);
+   for (int e = 0; e < fes->GetNE(); e++)
+   {
+      ElementTransformation *Tr = fes->GetElementTransformation(e);
+      DenseMatrix &Dmat = Dtensor(e);
+      for (int k = 0; k < quads; k++)
+      {
+         const IntegrationPoint &ip = IntRule->IntPoint(k);
+         Tr->SetIntPoint(&ip);
+         const double weight = ip.weight * Tr->Weight();
+         if (vcoeff != NULL)
+         {
+            vcoeff->Eval(cv, *Tr, ip);
+         }
+         for (int v = 0; v < vdim; v++)
+         {
+            Dmat(k, v) = weight;
+            if (coeff != NULL) Dmat(k, v) *= coeff->Eval(*Tr, ip);
+            else if (vcoeff != NULL)
+            {
+               Dmat(k, v) *= cv(v);
+            }
+         }
+      }
+   }
+}
+
+void PAMassIntegrator::MultSeg(const Vector &V, Vector &U)
+{
+   const int dofs1d = shape1d.Height();
+   const int quads1d = shape1d.Width();
+   const int quads = quads1d;
+   const int vdim = fes->GetVDim();
+
+   Vector Q(quads1d);
+
+   int offset = 0;
+   for (int e = 0; e < fes->GetNE(); ++e)
+   {
+      DenseMatrix &Dmat = Dtensor(e);
+      for (int vd = 0; vd < vdim; ++vd)
+      {
+         const Vector Vmat(V.GetData() + offset, dofs1d);
+         Vector Umat(U.GetData() + offset, dofs1d);
+
+         // Q_k1 = dshape_j1_k1 * V_i1
+         shape1d.MultTranspose(Vmat, Q);
+
+         double *data_q = Q.GetData();
+         const double *data_d = Dmat.GetColumn(vd);
+         for (int k = 0; k < quads; ++k) { data_q[k] *= data_d[k]; }
+
+         // Q_k1 = dshape_j1_k1 * Q_k1
+         shape1d.AddMult(Q, Umat);
+
+         // Increment offset into E-vectors.
+         offset += dofs1d;
+      }
+   }
+}
+
+void PAMassIntegrator::MultQuad(const Vector &V, Vector &U)
+{
+   const int dofs1d = shape1d.Height();
+   const int quads1d = shape1d.Width();
+
+   const int dofs   = dofs1d * dofs1d;
+   const int quads  = IntRule->GetNPoints();
+   const int vdim = fes->GetVDim();
+
+   DenseMatrix QQ(quads1d, quads1d);
+   DenseMatrix DQ(dofs1d, quads1d);
+
+   int offset = 0;
+   for (int e = 0; e < fes->GetNE(); ++e)
+   {
+      DenseMatrix &Dmat = Dtensor(e);
+      for (int vd = 0; vd < vdim; ++vd)
+      {
+         const DenseMatrix Vmat(V.GetData() + offset, dofs1d, dofs1d);
+         DenseMatrix Umat(U.GetData() + offset, dofs1d, dofs1d);
+
+         // DQ_j2_k1   = E_j1_j2  * dshape_j1_k1 -- contract in x direction
+         // QQ_0_k1_k2 = DQ_j2_k1 * shape_j2_k2  -- contract in y direction
+         MultAtB(Vmat, shape1d, DQ);
+         MultAtB(DQ, shape1d, QQ);
+
+         // QQ_c_k1_k2 = Dmat_c_d_k1_k2 * QQ_d_k1_k2
+         // NOTE: (k1, k2) = k -- 1d index over tensor product of quad points
+         double *data_qq = QQ.GetData();
+         const double *data_d = Dmat.GetColumn(vd);
+         for (int k = 0; k < quads; ++k) { data_qq[k] *= data_d[k]; }
+
+         // DQ_i2_k1   = shape_i2_k2  * QQ_0_k1_k2
+         // U_i1_i2   += dshape_i1_k1 * DQ_i2_k1
+         MultABt(shape1d, QQ, DQ);
+         AddMultABt(shape1d, DQ, Umat);
+
+         // increment offset
+         offset += dofs;
+      }
+   }
+}
+
+void PAMassIntegrator::MultHex(const Vector &V, Vector &U)
+{
+   const int dofs1d = shape1d.Height();
+   const int quads1d = shape1d.Width();
+
+   const int dofs   = dofs1d * dofs1d * dofs1d;
+   const int quads  = IntRule->GetNPoints();
+   const int vdim = fes->GetVDim();
+
+   Vector Q(quads1d);
+   DenseMatrix QQ(quads1d, quads1d);
+   DenseTensor QQQ(quads1d, quads1d, quads1d);
+
+   int offset = 0;
+   for (int e = 0; e < fes->GetNE(); ++e)
+   {
+      DenseMatrix &Dmat = Dtensor(e);
+      for (int vd = 0; vd < vdim; ++vd)
+      {
+         const DenseTensor Vmat(V.GetData() + offset, dofs1d, dofs1d, dofs1d);
+         DenseTensor Umat(U.GetData() + offset, dofs1d, dofs1d, dofs1d);
+
+         // QQQ_k1_k2_k3 = shape_j1_k1 * shape_j2_k2  * shape_j3_k3  * Vmat_j1_j2_j3
+         QQQ = 0.;
+         for (int j3 = 0; j3 < dofs1d; ++j3)
+         {
+            QQ = 0.;
+            for (int j2 = 0; j2 < dofs1d; ++j2)
+            {
+               Q = 0.;
+               for (int j1 = 0; j1 < dofs1d; ++j1)
+               {
+                  for (int k1 = 0; k1 < quads1d; ++k1)
+                  {
+                     Q(k1) += Vmat(j1, j2, j3) * shape1d(j1, k1);
+                  }
+               }
+               for (int k2 = 0; k2 < quads1d; ++k2)
+                  for (int k1 = 0; k1 < quads1d; ++k1)
+                  {
+                     QQ(k1, k2) += Q(k1) * shape1d(j2, k2);
+                  }
+            }
+            for (int k3 = 0; k3 < quads1d; ++k3)
+               for (int k2 = 0; k2 < quads1d; ++k2)
+                  for (int k1 = 0; k1 < quads1d; ++k1)
+                  {
+                     QQQ(k1, k2, k3) += QQ(k1, k2) * shape1d(j3, k3);
+                  }
+         }
+
+         // QQQ_k1_k2_k3 = Dmat_k1_k2_k3 * QQQ_k1_k2_k3
+         // NOTE: (k1, k2, k3) = q -- 1d quad point index
+         double *data_qqq = QQQ.GetData(0);
+         const double *data_d = Dmat.GetColumn(vd);
+         for (int k = 0; k < quads; ++k) { data_qqq[k] *= data_d[k]; }
+
+         // Apply transpose of the first operator that takes V -> QQQ -- QQQ -> U
+         for (int k3 = 0; k3 < quads1d; ++k3)
+         {
+            QQ = 0.;
+            for (int k2 = 0; k2 < quads1d; ++k2)
+            {
+               Q = 0.;
+               for (int k1 = 0; k1 < quads1d; ++k1)
+               {
+                  for (int i1 = 0; i1 < dofs1d; ++i1)
+                  {
+                     Q(i1) += QQQ(k1, k2, k3) * shape1d(i1, k1);
+                  }
+               }
+               for (int i2 = 0; i2 < dofs1d; ++i2)
+                  for (int i1 = 0; i1 < dofs1d; ++i1)
+                  {
+                     QQ(i1, i2) += Q(i1) * shape1d(i2, k2);
+                  }
+            }
+            for (int i3 = 0; i3 < dofs1d; ++i3)
+               for (int i2 = 0; i2 < dofs1d; ++i2)
+                  for (int i1 = 0; i1 < dofs1d; ++i1)
+                  {
+                     Umat(i1, i2, i3) += shape1d(i3, k3) * QQ(i1, i2);
+                  }
+         }
+
+         // increment offset
+         offset += dofs;
+      }
+   }
+}
+
+void PAMassIntegrator::AddMult(const Vector &x, Vector &y)
+{
+   const int dim = fes->GetMesh()->Dimension();
+
+   switch (dim)
+   {
+   case 1: MultSeg(x, y); break;
+   case 2: MultQuad(x, y); break;
+   case 3: MultHex(x, y); break;
+   default: mfem_error("Not yet supported"); break;
+   }
+}
+
+LinearFESpaceIntegrator *PAIntegratorMap::DomainIntegrator(BilinearFormIntegrator *integ) const
+{
+   {
+      DiffusionIntegrator *actual_integ = dynamic_cast<DiffusionIntegrator*>(integ);
+      if (actual_integ) { return new PADiffusionIntegrator(actual_integ); }
+   }
+   {
+      MassIntegrator *actual_integ = dynamic_cast<MassIntegrator*>(integ);
+      if (actual_integ) { return new PAMassIntegrator(actual_integ); }
+   }
+   {
+      VectorMassIntegrator *actual_integ = dynamic_cast<VectorMassIntegrator*>(integ);
+      if (actual_integ) { return new PAMassIntegrator(actual_integ); }
+   }
+
+   mfem_error("Not supported.");
+   return NULL;
+}
+
+
+}

fem/painteg.hpp

@@ -0,0 +1,90 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+// This file contains FESpaceIntegrators.
+
+
+#ifndef MFEM_PAINTEG
+#define MFEM_PAINTEG
+
+#include "../config/config.hpp"
+#include "nonlininteg.hpp"
+#include "bilinearformoper.hpp"
+
+namespace mfem
+{
+
+// These integrators use constructors based on the non-PA versions so
+// that the options are consistent. If that is not the case, the
+// friendship can be revoked and those constructors can be removed.
+
+/** Class for computing the action of (grad(u), grad(v)) from a scalar
+ * fespace using a partially assembled operator at quadrature
+ * points. */
+class PADiffusionIntegrator : public LinearFESpaceIntegrator
+{
+protected:
+   // Carry pointer in order to have access to coefficient
+   DiffusionIntegrator *integ;     // Own this
+   const FiniteElementSpace *fes;  // TODO: support mixed spaces
+   DenseTensor Dtensor;
+   DenseMatrix shape1d, dshape1d;
+
+   // Action methods
+   void MultSeg(const Vector &V, Vector &U);
+   void MultQuad(const Vector &V, Vector &U);
+   void MultHex(const Vector &V, Vector &U);
+
+public:
+   PADiffusionIntegrator(DiffusionIntegrator *_integ) : integ(_integ) {}
+   ~PADiffusionIntegrator() { delete integ; }
+
+   virtual void Assemble(FiniteElementSpace *trial_fes,
+                         FiniteElementSpace *test_fes);
+
+   virtual void AddMult(const Vector &x, Vector &y);
+};
+
+/** Class for computing the action of (u, v) from a scalar fespace
+ * using a partially assembled operator at quadrature points. */
+class PAMassIntegrator : public LinearFESpaceIntegrator
+{
+protected:
+   MassIntegrator *integ;     // Own this
+   VectorMassIntegrator *vinteg;     // Own this
+   const FiniteElementSpace *fes;  // TODO: support mixed spaces
+   DenseTensor Dtensor;
+   DenseMatrix shape1d;
+
+   // Action methods
+   void MultSeg(const Vector &V, Vector &U);
+   void MultQuad(const Vector &V, Vector &U);
+   void MultHex(const Vector &V, Vector &U);
+
+public:
+   PAMassIntegrator(MassIntegrator *_integ) : integ(_integ), vinteg(NULL) {}
+   PAMassIntegrator(VectorMassIntegrator *_vinteg) : integ(NULL), vinteg(_vinteg) {}
+   ~PAMassIntegrator() { delete integ; delete vinteg; }
+
+   virtual void Assemble(FiniteElementSpace *_trial_fes,
+                         FiniteElementSpace *_test_fes);
+
+   virtual void AddMult(const Vector &x, Vector &y);
+};
+
+struct PAIntegratorMap : public IntegratorMap
+{
+   virtual LinearFESpaceIntegrator *DomainIntegrator(BilinearFormIntegrator *integ) const;
+};
+
+}
+
+#endif

fem/partialassemblykernel.hpp

@@ -0,0 +1,266 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+// This file contains operator-based bilinear form integrators used
+// with BilinearFormOperator.
+
+#ifndef MFEM_PAK
+#define MFEM_PAK
+
+#include "fem.hpp"
+#include "../config/config.hpp"
+#include "bilininteg.hpp"
+#include "dalg.hpp"
+#include "dgfacefunctions.hpp"
+#include "domainkernels.hpp"
+#include "facekernels.hpp"
+#include "solverkernels.hpp"
+#include <iostream>
+
+namespace mfem
+{
+
+/////////////////////////////////////////////////
+//                                             //
+//                                             //
+//        PARTIAL ASSEMBLY INTEGRATORS         //
+//                                             //
+//                                             //
+/////////////////////////////////////////////////
+
+
+  /////////////////////////////
+ // Domain Kernel Interface //
+/////////////////////////////
+
+struct ElementInfo
+{
+   int dim;
+   int k;
+   int e;
+   ElementTransformation* tr;
+   IntegrationPoint ip;
+   Tensor<2> J_ek;
+};
+
+/**
+*  A partial assembly Integrator class for domain integrals.
+*  Takes an 'Equation' template parameter, that must contain 'OpName' of
+*  type 'PAOp' and a function named 'evalD', that receives a 'res' vector,
+*  the element transformation and the integration point, and then whatever
+*  is needed to compute at the point (Coefficient, VectorCoeffcient, etc...).
+*  The 'IMPL' template parameter allows to switch between different implementations
+*  of the tensor contraction kernels.
+*/
+template < typename Equation,
+            template<typename,PAOp> class IMPL = DomainMult>
+class PADomainInt
+: public LinearFESpaceIntegrator, public IMPL<Equation,Equation::OpName>, public Operator
+{
+private:
+   typedef IMPL<Equation,Equation::OpName> Op;
+
+public:
+   /**
+   *  The constructor is templated so that the argument needed for 'evalD' can be
+   *  packed arbitrarily ('evalD' with the corresponding signature must exist).
+   */
+   template <typename Args>
+   PADomainInt(FiniteElementSpace *fes, const int order, const Args& args)
+   : LinearFESpaceIntegrator(&IntRules.Get(fes->GetFE(0)->GetGeomType(), order)),
+     Op(fes,order,args),
+     Operator()
+   {
+      const int nb_elts = fes->GetNE();
+      const int quads  = IntRule->GetNPoints();
+      const FiniteElement* fe = fes->GetFE(0);
+      const int dim = fe->GetDim();
+      this->InitD(dim,quads,nb_elts);
+      Tensor<1> Jac1D(dim*dim*quads*nb_elts);
+      EvalJacobians(dim,fes,order,Jac1D);
+      Tensor<4> Jac(Jac1D.getData(),dim,dim,quads,nb_elts);
+      for (int e = 0; e < nb_elts; ++e)
+      {
+         ElementTransformation *Tr = fes->GetElementTransformation(e);
+         for (int k = 0; k < quads; ++k)
+         {
+            Tensor<2> J_ek(&Jac(0,0,k,e),dim,dim);
+            const IntegrationPoint &ip = IntRule->IntPoint(k);
+            Tr->SetIntPoint(&ip);
+            this->evalEq(dim, k, e, Tr, ip, J_ek, args);
+         }
+      }
+   }
+
+   const typename Op::DTensor& getD() const
+   {
+      return Op::getD();
+   }
+
+   /**
+   *  Applies the partial assembly operator.
+   */
+   virtual void AddMult(const Vector &fun, Vector &vect)
+   {
+      int dim = this->fes->GetFE(0)->GetDim();
+      switch(dim)
+      {
+      case 1:this->Mult1d(fun,vect); break;
+      case 2:this->Mult2d(fun,vect); break;
+      case 3:this->Mult3d(fun,vect); break;
+      default: mfem_error("More than # dimension not yet supported"); break;
+      }
+   }
+
+   virtual void Mult(const Vector &fun, Vector &vect) const{
+      int dim = this->fes->GetFE(0)->GetDim();
+      switch(dim)
+      {
+      case 1:this->Mult1d(fun,vect); break;
+      case 2:this->Mult2d(fun,vect); break;
+      case 3:this->Mult3d(fun,vect); break;
+      default: mfem_error("More than # dimension not yet supported"); break;
+      }
+   }
+};
+
+  ///////////////////////////
+ // Face Kernel Interface //
+///////////////////////////
+
+struct FaceInfo
+{
+   int dim; // The problem dimension
+   int k1, k2; // The indices of 
+   IntegrationPoint eip1, eip2; // The integration points on each element
+   Vector* normal;  // The normal to the face
+   int ind_elt1, ind_elt2; // The indices of the elements
+   int face_id1, face_id2; // The face ID for the face according to each element
+   FaceElementTransformations* face_tr; // The Face transformation
+   Tensor<2> J_e1, J_e2; // The Jacobians for each element at their respective quadrature point
+};
+
+/**
+*  A partial assembly Integrator interface class for face integrals.
+*  The template parameters have the same role as for 'PADomainInt'.
+*/
+template <typename Equation, template<typename,PAOp> class IMPL = FaceMult>
+class PAFaceInt
+: public LinearFESpaceIntegrator, public IMPL<Equation,Equation::FaceOpName>
+{
+private:
+   typedef IMPL<Equation,Equation::FaceOpName> Op;
+
+public:
+   template <typename Args>
+   PAFaceInt(FiniteElementSpace* fes, const int order, Args& args)
+   : LinearFESpaceIntegrator(&IntRules.Get(fes->GetFE(0)->GetGeomType(), order)),
+     Op(fes, order, args)
+   {
+      const int dim = fes->GetFE(0)->GetDim();
+      const int quads1d = fes->GetNQuads1d(order);
+      Mesh* mesh = fes->GetMesh();
+      const int nb_elts = fes->GetNE();
+      const int nb_faces_elt = 2*dim;
+      const int nb_faces = mesh->GetNumFaces();
+      int geom;
+      switch(dim){
+         case 1:geom = Geometry::POINT;break;
+         case 2:geom = Geometry::SEGMENT;break;
+         case 3:geom = Geometry::SQUARE;break;
+      }
+      const IntegrationRule& ir = IntRules.Get(geom, order);
+      const int quads  = ir.GetNPoints();
+      Vector qvec(dim);
+      Tensor<1> normal(dim);
+      Vector n(normal.getData(),dim);
+      // Vector n(dim);
+      this->init(dim,quads,nb_elts,nb_faces_elt);
+      // !!! Should not be recomputed... !!!
+      Tensor<1> Jac1D(dim*dim*quads*quads1d*nb_elts);
+      EvalJacobians(dim,fes,order,Jac1D);
+      Tensor<4> Jac(Jac1D.getData(),dim,dim,quads*quads1d,nb_elts);// Creating a view
+      // !!!                             !!!
+      // We have a per face approach for the fluxes
+      for (int face = 0; face < nb_faces; ++face)
+      {
+         int ind_elt1, ind_elt2;
+         int face_id1, face_id2;
+         int nb_rot1, nb_rot2;
+         GetFaceInfo(mesh, face, ind_elt1, ind_elt2, face_id1, face_id2, nb_rot1, nb_rot2);
+         FaceElementTransformations* face_tr = mesh->GetFaceElementTransformations(face);
+         int perm1, perm2;
+         // cout << "ind_elt1=" << ind_elt1 << ", face_id1=" << face_id1 << ", nb_rot1=" << nb_rot1 << ", ind_elt2=" << ind_elt2 << ", face_id2=" << face_id2 << ", nb_rot2=" << nb_rot2 << endl;
+         for (int kf = 0; kf < quads; ++kf)
+         {
+            const IntegrationPoint& ip = ir.IntPoint(kf);
+            if(ind_elt2!=-1){//Not a boundary face
+               Tensor<1,int> ind_f1(dim-1), ind_f2(dim-1);
+               // We compute the lexicographical index on each face
+               int k1 = GetFaceQuadIndex(dim,face_id1,nb_rot1,kf,quads1d,ind_f1);
+               int k2 = GetFaceQuadIndex(dim,face_id2,nb_rot2,kf,quads1d,ind_f2);
+               this->initFaceData(dim,ind_elt1,face_id1,nb_rot1,perm1,ind_elt2,face_id2,nb_rot2,perm2);
+               face_tr->Face->SetIntPoint( &ip );
+               IntegrationPoint eip1;
+               face_tr->Loc1.Transform(ip,eip1);
+               eip1.weight = ip.weight;//Sets the weight since Transform doesn't do it...
+               // face_tr->Elem1->SetIntPoint( &eip1 );
+               IntegrationPoint eip2;
+               face_tr->Loc2.Transform(ip,eip2);
+               eip2.weight = ip.weight;//Sets the weight since Transform doesn't do it...
+               // face_tr->Elem2->SetIntPoint( &eip2 );
+               int kg1 = GetGlobalQuadIndex(dim,face_id1,quads1d,ind_f1);
+               int kg2 = GetGlobalQuadIndex(dim,face_id2,quads1d,ind_f2);
+               Tensor<2> J_e1(&Jac(0,0,kg1,ind_elt1),dim,dim);
+               Tensor<2> J_e2(&Jac(0,0,kg2,ind_elt2),dim,dim);
+               Tensor<2> Adj(dim,dim);
+               adjugate(J_e1,Adj);
+               calcOrtho( Adj, face_id1, normal); // normal*determinant (risky, bug prone)
+               this->evalEq(dim,k1,k2,n,ind_elt1,face_id1,ind_elt2,face_id2,face_tr,eip1,eip2,J_e1,J_e2,args);
+               // FaceInfo face_info = {dim,k1,k2,n,ind_elt1,face_id1,ind_elt2,face_id2,face_tr,eip1,eip2,J_e1,J_e2};
+               // this->evalEq(face_info,args);
+            }else{//Boundary face
+               this->initBoundaryFaceData(ind_elt1,face_id1);
+               // TODO: Something should be done here when there is boundary conditions!
+               // D11(ind) = 0;  
+            }
+         }
+      }
+   }
+
+   // Perform the action of the BilinearFormIntegrator
+   virtual void AddMult(const Vector &fun, Vector &vect)
+   {
+      int dim = this->fes->GetFE(0)->GetDim();
+      switch(dim)
+      {
+         case 1:
+            mfem_error("Not yet implemented");
+            break;
+         case 2:
+            this->EvalInt2D(fun, vect);
+            this->EvalExt2D(fun, vect);
+            break;
+         case 3:
+            this->EvalInt3D(fun, vect);
+            this->EvalExt3D(fun, vect);
+            break;
+         default:
+            mfem_error("Face Kernel does not exist for this dimension.");
+            break;
+      }
+   }
+
+};
+
+}
+
+#endif //MFEM_PAK

fem/pbilinearform.cpp

@@ -239,6 +239,20 @@ void ParBilinearForm::Assemble(int skip_zeros)
    }
 }
 
+void ParBilinearForm::AssembleForm(BilinearFormOperator &A, int skip_zeros)
+{
+   A.Assemble(this);
+   oper = &A;
+   oper_type = MFEM_FORMOPER;
+}
+
+void ParBilinearForm::AssembleForm(HypreParMatrix &A, int skip_zeros)
+{
+   Assemble(skip_zeros);
+   oper = &A;
+   oper_type = Hypre_ParCSR;
+}
+
 void ParBilinearForm
 ::ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
                                HypreParMatrix &A, const HypreParVector &X,
@@ -324,6 +338,43 @@ void ParBilinearForm::FormLinearSystem(
    }
 }
 
+template <typename OpType>
+void ParBilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x, Vector &b,
+                                       OpType &A, Vector &X, Vector &B,
+                                       int copy_interior)
+{
+   OperatorHandle Ah;
+   FormLinearSystem(ess_tdof_list, x, b, Ah, X, B, copy_interior);
+   OpType *A_ptr = Ah.Is<OpType>();
+   MFEM_VERIFY(A_ptr, "invalid OpType used");
+   A.MakeRef(*A_ptr);
+}
+
+template <>
+void ParBilinearForm::FormLinearSystem<Operator*>(const Array<int> &ess_tdof_list, Vector &x, Vector &b,
+                                                  Operator * &A, Vector &X, Vector &B,
+                                                  int copy_interior)
+{
+   if (oper_type == Hypre_ParCSR)
+   {
+      HypreParMatrix &Amat = static_cast<HypreParMatrix&>(*oper);
+      FormLinearSystem(ess_tdof_list, x, b, Amat,
+                       X, B, copy_interior);
+      HypreParMatrix *M = new HypreParMatrix;
+      M->MakeRef(Amat);
+      A = M;
+   }
+   else if (oper_type == MFEM_FORMOPER)
+   {
+      oper->FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
+   }
+   else
+   {
+      mfem_error("Not supported.");
+   }
+}
+
+
 void ParBilinearForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
                                        OperatorHandle &A)
 {

fem/pbilinearform.hpp

@@ -70,6 +70,9 @@ public:
    /// Assemble the local matrix
    void Assemble(int skip_zeros = 1);
 
+   void AssembleForm(BilinearFormOperator &A, int skip_zeros = 1);
+   void AssembleForm(HypreParMatrix &A, int skip_zeros = 1);
+
    /// Returns the matrix assembled on the true dofs, i.e. P^t A P.
    /** The returned matrix has to be deleted by the caller. */
    HypreParMatrix *ParallelAssemble() { return ParallelAssemble(mat); }
@@ -178,14 +181,7 @@ public:
    template <typename OpType>
    void FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x, Vector &b,
                          OpType &A, Vector &X, Vector &B,
-                         int copy_interior = 0)
-   {
-      OperatorHandle Ah;
-      FormLinearSystem(ess_tdof_list, x, b, Ah, X, B, copy_interior);
-      OpType *A_ptr = Ah.Is<OpType>();
-      MFEM_VERIFY(A_ptr, "invalid OpType used");
-      A.MakeRef(*A_ptr);
-   }
+                         int copy_interior = 0);
 
    /// Form the linear system matrix @a A, see FormLinearSystem() for details.
    void FormSystemMatrix(const Array<int> &ess_tdof_list, OperatorHandle &A);

fem/solverkernels.hpp

@@ -0,0 +1,647 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+// This file contains operator-based bilinear form integrators used
+// with BilinearFormOperator.
+
+
+#ifndef MFEM_SOLVERKERNELS
+#define MFEM_SOLVERKERNELS
+
+#include "dalg.hpp"
+#include "dgpabilininteg.hpp"
+#include "tensorialfunctions.hpp"
+#include "partialassemblykernel.hpp"
+
+namespace mfem
+{
+
+template <typename Op>
+class CGSolverDG: public Operator
+{
+public:
+   typedef Tensor<2> Tensor2d;
+
+protected:
+   FiniteElementSpace& fes;
+   const Tensor2d& D;
+   Tensor2d shape1d;
+   const double treshold;
+
+public:
+   CGSolverDG(FiniteElementSpace& fes, int order, const Op& op, const double tr = 1e-10)
+   : Operator(fes.GetVSize()), fes(fes), D(op.getD()),
+     shape1d(fes.GetNDofs1d(),fes.GetNQuads1d(order)), treshold(tr)
+   {
+      ComputeBasis1d(fes.GetFE(0), order, shape1d);
+   }
+
+   virtual void Mult(const Vector &U, Vector &V) const
+   {
+      switch(fes.GetFE(0)->GetDim())
+      {
+         case 1:
+            Mult1d(U, V);
+            break;
+         case 2:
+            Mult2d(U, V);
+            break;
+         case 3:
+            Mult3d(U, V);
+            break;
+      }
+   }
+
+protected:
+   /**
+   *  The domain Kernels for BtDB in 1d,2d and 3d.
+   */
+   void Mult1d(const Vector &V, Vector &U) const;
+   void Mult2d(const Vector &V, Vector &U) const;
+   void Mult3d(const Vector &V, Vector &U) const;  
+
+};
+
+template <typename Mass, typename Prec>
+class PrecCGSolverDG: public Operator
+{
+public:
+   typedef Tensor<2> Tensor2d;
+
+protected:
+   FiniteElementSpace& fes;
+   const Tensor2d& D;
+   Prec& prec;
+   Tensor2d shape1d;
+   const double treshold;
+
+public:
+   PrecCGSolverDG(FiniteElementSpace& fes, int order, Mass& mass, Prec& prec, const double tr = 1e-10)
+   : Operator(fes.GetVSize()), fes(fes), D(mass.getD()), prec(prec),
+     shape1d(fes.GetNDofs1d(),fes.GetNQuads1d(order)), treshold(tr)
+   {
+      ComputeBasis1d(fes.GetFE(0), order, shape1d);
+   }
+
+   virtual void Mult(const Vector &U, Vector &V) const
+   {
+      switch(fes.GetFE(0)->GetDim())
+      {
+         case 1:
+            Mult1d(U, V);
+            break;
+         case 2:
+            Mult2d(U, V);
+            break;
+         case 3:
+            Mult3d(U, V);
+            break;
+      }
+   }
+
+protected:
+   void Mult1d(const Vector &V, Vector &U) const;
+   void Mult2d(const Vector &V, Vector &U) const;
+   void Mult3d(const Vector &V, Vector &U) const;  
+
+};
+
+template <typename Mass>
+class PACGSolver: public Operator
+{
+private:
+   const int dofs;
+   Mass& mass;
+   const double treshold;
+
+public:
+   PACGSolver(const FiniteElementSpace* fes, Mass& mass, const double tr = 1e-10)
+   : Operator(fes->GetVSize()), dofs(fes->GetNDofs()), mass(mass), treshold(tr)
+   {
+
+   }
+
+   virtual void Mult(const Vector &U, Vector &V) const
+   {
+      int iter = 0;
+
+      Vector b(dofs);
+      b = U;
+      Vector& x = V;
+      x = 0.0;
+      Vector r(dofs);
+      for (int i = 0; i < dofs; ++i)
+      {
+         r(i) = -b(i);
+      }
+      Vector p(dofs);
+      for (int i = 0; i < dofs; ++i)
+      {
+         p(i) = - r(i);
+      }
+      double rsold = r*r;
+      for(iter=0; iter<dofs; iter++){
+         Vector Ap(dofs);
+         Ap = 0.0;
+         // Ap = A * p
+         mass.AddMult(p,Ap);
+         const double alpha = rsold / (p*Ap);
+         for (int i = 0; i < dofs; ++i)
+         {
+            x(i) = x(i) + alpha * p(i);
+            r(i) = r(i) + alpha * Ap(i);
+         }
+         // rsnew = y*r;
+         const double rsnew = r*r;
+         if (sqrt(rsnew)<treshold) break;
+         // p = r + (rsnew/rsold) * p
+         for (int i = 0; i < dofs; ++i)
+         {
+            p(i) = -r(i) + (rsnew/rsold) * p(i);
+         }
+         rsold = rsnew;
+      }
+      // cout << "residual=" << rsold << endl;
+      // cout << "iter=" << iter << endl;
+   }
+};
+
+template <typename Mass, typename Prec>
+class PAPrecCGSolver: public Operator
+{
+private:
+   const int dofs;
+   Mass& mass;
+   Prec& prec;
+   const double treshold;
+
+public:
+   PAPrecCGSolver(const FiniteElementSpace* fes, Mass& mass, Prec& prec, const double tr = 1e-10)
+   : Operator(fes->GetVSize()), dofs(fes->GetNDofs()), mass(mass), prec(prec), treshold(tr)
+   {
+
+   }
+
+   virtual void Mult(const Vector &U, Vector &V) const
+   {
+      int iter = 0;
+
+      Vector b(dofs);
+      b = U;
+      Vector& x = V;
+      x = 0.0;
+      Vector r(dofs);
+      for (int i = 0; i < dofs; ++i)
+      {
+         r(i) = -b(i);
+      }
+      Vector y(dofs);
+      prec.Mult(r,y);
+      Vector p(dofs);
+      for (int i = 0; i < dofs; ++i)
+      {
+         p(i) = - y(i);
+      }
+      double rsold = y*r;
+      for(iter=0; iter<dofs; iter++){
+         Vector Ap(dofs);
+         Ap = 0.0;
+         mass.AddMult(p,Ap);
+         const double alpha = rsold / (p*Ap);
+         for (int i = 0; i < dofs; ++i)
+         {
+            x(i) = x(i) + alpha * p(i);
+            r(i) = r(i) + alpha * Ap(i);
+         }
+         prec.Mult(r,y);
+         const double rsnew = y*r;
+         if (sqrt(rsnew)<treshold){
+            // cout << "residual=" << rsnew << endl;
+            // cout << "iter=" << iter << endl;
+            break;  
+         }
+         const double beta = rsnew/rsold;
+         for (int i = 0; i < dofs; ++i)
+         {
+            p(i) = -y(i) + beta * p(i);
+         }
+         rsold = rsnew;
+      }
+   }
+};
+
+/**
+*  A class that implement the BtDB partial assembly Kernel
+*/
+class DiagSolverDG: public Operator
+{
+private:
+   const int nbelts;
+   const int dofs;
+   const int size;
+   Tensor<2> D;
+
+public:
+   template <typename Op>
+   DiagSolverDG(FiniteElementSpace fes, int order, Op& op, bool fast_eval = false)
+   : nbelts( fes.GetNE() ), dofs( fes.GetFE(0)->GetDof() ), size(nbelts*dofs), D(dofs,nbelts),
+     Operator(fes.GetVSize())
+   {
+      if(fast_eval)// FIXME: For some reason does not work in 2d...
+      {
+         Vector U(nbelts*dofs);
+         for (int i = 0; i < size; ++i)
+         {
+            U(i) = 1.0;
+         }
+         //If the op is diagonal, then we obtain the diagonal by multiplying by a vector of 1.
+         Vector V(D.getData(),size);
+         op.AddMult(U,V);
+      }else{
+         GetDiag(fes,order,op,D);
+      }
+   }
+
+   /**
+   *  The domain Kernels for BtDB in 1d,2d and 3d.
+   */
+   virtual void Mult(const Vector &U, Vector &V) const
+   {
+      const Tensor<2> Ut(U.GetData(), dofs, nbelts);
+      Tensor<2> Vt(V.GetData(), dofs, nbelts);
+      for (int e = 0; e < nbelts; e++)
+      {
+         for (int k = 0; k < dofs; ++k)
+         {
+            Vt(k,e) =  Ut(k,e)/D(k,e);
+         }
+      }
+   }
+
+   template <typename Vector>
+   void Mult(const int elt, const Vector &U, Vector &V) const
+   {
+      for (int k = 0; k < dofs; ++k)
+      {
+         V(k) =  U(k)/D(k,elt);
+      }
+   }
+
+};
+
+template<typename Op>
+void CGSolverDG<Op>::Mult1d(const Vector &V, Vector &U) const
+{
+   const int dofs1d = shape1d.Height();
+   const int quads1d = shape1d.Width();
+   const int quads = quads1d;
+   const int dofs = dofs1d;
+
+   for (int e = 0; e < fes.GetNE(); e++)
+   {
+      const Tensor<1> b(V.GetData() + e*dofs, dofs);
+      const Tensor<1> eD(D.getData() + e*quads, quads1d);
+      Tensor<1> x(U.GetData() + e*dofs, dofs1d);
+      x.zero();
+      Tensor<1> r(dofs);
+      r = b;
+      Tensor<1> p(dofs);
+      p = r;
+      double rsold = norm2sq(r);
+
+      for(int i=0; i<dofs; i++){
+         Tensor<1> Ap(dofs),tmp(quads1d);
+         // Ap = A * p
+         contract(shape1d,p,tmp);
+         cWiseMult(eD,tmp,tmp);
+         contract(shape1d,tmp,Ap);
+         const double alpha = rsold / dot(p,Ap);
+         // x = x + alpha * p
+         // r = r - alpha * Ap
+         for (int i = 0; i < dofs; ++i)
+         {
+            x(i) = x(i) + alpha * p(i);
+            r(i) = r(i) - alpha * Ap(i);
+         }
+         const double rsnew = norm2sq(r);
+         if (sqrt(rsnew)<treshold) break;
+         // p = r + (rsnew/rsold) * p
+         for (int i = 0; i < dofs1d; ++i)
+         {
+            p(i) = r(i) + (rsnew/rsold) * p(i);
+         }
+         rsold = rsnew;
+      }
+   }
+}
+
+template<typename Op>
+void CGSolverDG<Op>::Mult2d(const Vector &V, Vector &U) const
+{
+   const int dofs1d = shape1d.Height();
+   const int quads1d = shape1d.Width();
+   const int quads = quads1d * quads1d;
+   const int dofs = dofs1d * dofs1d;
+
+   Tensor<1> r(dofs);
+   Tensor<1> p(dofs);
+   Tensor<2> pT(p.getData(),dofs1d,dofs1d);
+   Tensor<1> Ap(dofs);
+   Tensor<2> ApT(Ap.getData(),dofs1d,dofs1d);
+   Tensor<2> tmp1(dofs1d,quads1d), tmp2(quads1d,quads1d), tmp3(quads1d,dofs1d);
+   for (int e = 0; e < fes.GetNE(); e++)
+   {
+      const Tensor<1> b(V.GetData() + e*dofs, dofs);
+      const Tensor<2> eD(D.getData() + e*quads, quads1d, quads1d);
+      Tensor<1> x(U.GetData() + e*dofs, dofs);
+      // Tensor<1> x(dofs);
+      x.zero();
+      // r = b - A x
+      for (int i = 0; i < dofs; ++i)
+      {
+         r(i) = b(i);
+      }
+      // p = r;
+      for (int i = 0; i < dofs; ++i)
+      {
+         p(i) = r(i);
+      }
+      double rsold = norm2sq(r);
+      for(int i=0; i<dofs; i++){
+         // Ap = A * p
+         contract(shape1d,pT,tmp1);
+         contract(shape1d,tmp1,tmp2);
+         cWiseMult(eD,tmp2,tmp2);
+         contractT(shape1d,tmp2,tmp3);
+         contractT(shape1d,tmp3,ApT);
+         const double alpha = rsold / dot(p,Ap);
+         // x = x + alpha * p
+         // r = r - alpha * Ap
+         for (int i = 0; i < dofs; ++i)
+         {
+            x(i) = x(i) + alpha * p(i);
+            r(i) = r(i) - alpha * Ap(i);
+         }
+         const double rsnew = norm2sq(r);
+         // cout << "elem " << e << " iter " << i << endl;
+         // cout << "residual = " << rsnew << endl;
+         if (sqrt(rsnew)<treshold) break;
+         // p = r + (rsnew/rsold) * p
+         for (int i = 0; i < dofs; ++i)
+         {
+            p(i) = r(i) + (rsnew/rsold) * p(i);
+         }
+         rsold = rsnew;
+      }
+   }
+}
+
+template<typename Op>
+void CGSolverDG<Op>::Mult3d(const Vector &V, Vector &U) const
+{
+   const int dofs1d = shape1d.Height();
+   const int quads1d = shape1d.Width();
+   const int quads = quads1d * quads1d * quads1d;
+   const int dofs = dofs1d * dofs1d * dofs1d;
+
+   for (int e = 0; e < fes.GetNE(); e++)
+   {
+      const Tensor<1> b(V.GetData() + e*dofs, dofs);
+      const Tensor<3> eD(D.getData() + e*quads, quads1d, quads1d, quads1d);
+      Tensor<1> x(U.GetData() + e*dofs, dofs);
+      x.zero();
+      Tensor<1> r(dofs);
+      r = b;
+      Tensor<1> p(dofs);
+      p = r;
+      double rsold = norm2sq(r);
+
+      for(int i=0; i<dofs; i++){
+         Tensor<3> pT(p.getData(),dofs1d,dofs1d,dofs1d);
+         Tensor<1> Ap(dofs);
+         Tensor<3> ApT(Ap.getData(),dofs1d,dofs1d,dofs1d);
+         Tensor<3> tmp1(dofs1d,dofs1d,quads1d), tmp2(dofs1d,quads1d,quads1d), tmp3(quads1d,quads1d,quads1d),
+                     tmp4(quads1d,quads1d,dofs1d), tmp5(quads1d,dofs1d,dofs1d);
+         // Ap = A * p
+         contract(shape1d,pT,tmp1);
+         contract(shape1d,tmp1,tmp2);
+         contract(shape1d,tmp2,tmp3);
+         cWiseMult(eD,tmp3,tmp3);
+         contract(shape1d,tmp3,tmp4);
+         contract(shape1d,tmp4,tmp5);
+         contract(shape1d,tmp5,ApT);
+         const double alpha = rsold / dot(p,Ap);
+         // x = x + alpha * p
+         // r = r - alpha * Ap
+         for (int i = 0; i < dofs; ++i)
+         {
+            x(i) = x(i) + alpha * p(i);
+            r(i) = r(i) - alpha * Ap(i);
+         }
+         const double rsnew = norm2sq(r);
+         if (sqrt(rsnew)<treshold) break;
+         // p = r + (rsnew/rsold) * p
+         for (int i = 0; i < dofs; ++i)
+         {
+            p(i) = r(i) + (rsnew/rsold) * p(i);
+         }
+         rsold = rsnew;
+      }
+   }
+}
+
+
+template <typename Mass, typename Prec>
+void PrecCGSolverDG<Mass,Prec>::Mult1d(const Vector &V, Vector &U) const
+{
+   const int dofs1d = shape1d.Height();
+   const int quads1d = shape1d.Width();
+   const int quads = quads1d;
+   const int dofs = dofs1d;
+
+   for (int e = 0; e < fes.GetNE(); e++)
+   {
+      const Tensor<1> b(V.GetData() + e*dofs, dofs);
+      const Tensor<1> eD(D.getData() + e*quads, quads1d);
+      Tensor<1> x(U.GetData() + e*dofs, dofs1d);
+      x.zero();
+      Tensor<1> r(dofs);
+      for (int i = 0; i < dofs; ++i)
+      {
+         r(i) = -b(i);
+      }
+      Tensor<1> y(dofs);
+      prec.Mult(e,r,y);
+      Tensor<1> p(dofs);
+      for (int i = 0; i < dofs; ++i)
+      {
+         p(i) = -y(i);
+      }
+      double rsold = dot(r,y);
+      for(int iter=0; iter<dofs; iter++){
+         Tensor<1> Ap(dofs),tmp(quads1d);
+         // Ap = A * p
+         contract(shape1d,p,tmp);
+         cWiseMult(eD,tmp,tmp);
+         contract(shape1d,tmp,Ap);
+         const double alpha = rsold / dot(p,Ap);
+         for (int i = 0; i < dofs; ++i)
+         {
+            x(i) = x(i) + alpha * p(i);
+            r(i) = r(i) + alpha * Ap(i);
+         }
+         prec.Mult(e,r,y);
+         const double rsnew = dot(r,y);
+         if (sqrt(rsnew)<treshold){
+            // cout << "residual=" << rsnew << endl;
+            // cout << "iter=" << iter << endl;            
+            break;
+         }
+         const double beta = rsnew/rsold;
+         for (int i = 0; i < dofs1d; ++i)
+         {
+            p(i) = -y(i) + beta * p(i);
+         }
+         rsold = rsnew;
+      }
+   }
+}
+
+template <typename Mass, typename Prec>
+void PrecCGSolverDG<Mass,Prec>::Mult2d(const Vector &V, Vector &U) const
+{
+   const int dofs1d = shape1d.Height();
+   const int quads1d = shape1d.Width();
+   const int quads = quads1d * quads1d;
+   const int dofs = dofs1d * dofs1d;
+
+   Tensor<1> r(dofs);
+   Tensor<1> y(dofs);
+   Tensor<1> p(dofs);
+   Tensor<2> pT(p.getData(),dofs1d,dofs1d);
+   Tensor<1> Ap(dofs);
+   Tensor<2> ApT(Ap.getData(),dofs1d,dofs1d);
+   Tensor<2> tmp1(dofs1d,quads1d), tmp2(quads1d,quads1d), tmp3(quads1d,dofs1d);
+   for (int e = 0; e < fes.GetNE(); e++)
+   {
+      const Tensor<1> b(V.GetData() + e*dofs, dofs);
+      const Tensor<2> eD(D.getData() + e*quads, quads1d, quads1d);
+      Tensor<1> x(U.GetData() + e*dofs, dofs);
+      x.zero();
+      for (int i = 0; i < dofs; ++i)
+      {
+         r(i) = -b(i);
+      }
+      prec.Mult(e,r,y);
+      for (int i = 0; i < dofs; ++i)
+      {
+         p(i) = -y(i);
+      }
+      double rsold = dot(y,r);
+      for(int iter=0; iter<dofs; iter++){
+         // Ap = A * p
+         contract(shape1d,pT,tmp1);
+         contract(shape1d,tmp1,tmp2);
+         cWiseMult(eD,tmp2,tmp2);
+         contractT(shape1d,tmp2,tmp3);
+         contractT(shape1d,tmp3,ApT);
+         const double alpha = rsold / dot(p,Ap);
+         for (int i = 0; i < dofs; ++i)
+         {
+            x(i) = x(i) + alpha * p(i);
+            r(i) = r(i) + alpha * Ap(i);
+         }
+         prec.Mult(e,r,y);
+         const double rsnew = dot(y,r);
+         if (sqrt(rsnew)<treshold){
+            // cout << "residual=" << rsnew << endl;
+            // cout << "iter=" << iter << endl;
+            break;  
+         }
+         const double beta = (rsnew/rsold);
+         for (int i = 0; i < dofs; ++i)
+         {
+            p(i) = -y(i) + beta * p(i);
+         }
+         rsold = rsnew;
+      }
+   }
+}
+
+template <typename Mass, typename Prec>
+void PrecCGSolverDG<Mass,Prec>::Mult3d(const Vector &V, Vector &U) const
+{
+   const int dofs1d = shape1d.Height();
+   const int quads1d = shape1d.Width();
+   const int quads = quads1d * quads1d * quads1d;
+   const int dofs = dofs1d * dofs1d * dofs1d;
+
+   Tensor<1> r(dofs);
+   Tensor<1> y(dofs);
+   Tensor<1> p(dofs);
+   Tensor<1> Ap(dofs);
+   Tensor<3> tmp1(dofs1d,dofs1d,quads1d), tmp2(dofs1d,quads1d,quads1d), tmp3(quads1d,quads1d,quads1d),
+               tmp4(quads1d,quads1d,dofs1d), tmp5(quads1d,dofs1d,dofs1d);
+   for (int e = 0; e < fes.GetNE(); e++)
+   {
+      const Tensor<1> b(V.GetData() + e*dofs, dofs);
+      const Tensor<3> eD(D.getData() + e*quads, quads1d, quads1d, quads1d);
+      Tensor<1> x(U.GetData() + e*dofs, dofs);
+      x.zero();
+      for (int i = 0; i < dofs; ++i)
+      {
+         r(i) = -b(i);
+      }
+      prec.Mult(e,r,y);
+      for (int i = 0; i < dofs; ++i)
+      {
+         p(i) = -y(i);
+      }
+      double rsold = dot(y,r);
+      for(int i=0; i<dofs; i++){
+         Tensor<3> pT(p.getData(),dofs1d,dofs1d,dofs1d);
+         Tensor<3> ApT(Ap.getData(),dofs1d,dofs1d,dofs1d);
+         // Ap = A * p
+         contract(shape1d,pT,tmp1);
+         contract(shape1d,tmp1,tmp2);
+         contract(shape1d,tmp2,tmp3);
+         cWiseMult(eD,tmp3,tmp3);
+         contract(shape1d,tmp3,tmp4);
+         contract(shape1d,tmp4,tmp5);
+         contract(shape1d,tmp5,ApT);
+         const double alpha = rsold / dot(p,Ap);
+         for (int i = 0; i < dofs; ++i)
+         {
+            x(i) = x(i) + alpha * p(i);
+            r(i) = r(i) + alpha * Ap(i);
+         }
+         prec.Mult(e,r,y);
+         const double rsnew = dot(y,r);
+         if (sqrt(rsnew)<treshold){
+            // cout << "residual=" << rsnew << endl;
+            // cout << "iter=" << iter << endl;
+            break;  
+         }
+         const double beta = (rsnew/rsold);
+         for (int i = 0; i < dofs; ++i)
+         {
+            p(i) = -y(i) + beta * p(i);
+         }
+         rsold = rsnew;
+      }
+   }
+}
+
+
+
+}
+
+#endif //MFEM_SOLVERKERNELS

fem/tensorialfunctions.hpp

@@ -0,0 +1,711 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+// This file contains operator-based bilinear form integrators used
+// with BilinearFormOperator.
+
+#ifndef MFEM_TENSORFUNC
+#define MFEM_TENSORFUNC
+
+#include "dalg.hpp"
+#include "fem.hpp"
+
+namespace mfem
+{
+
+/**
+* Gives the evaluation of the 1d basis functions and their derivative at one point @param x
+*/
+template <typename Tensor>
+void ComputeBasis0d(const FiniteElement *fe, double x,
+                     Tensor& shape0d, Tensor& dshape0d)
+{
+   const TensorBasisElement* tfe(dynamic_cast<const TensorBasisElement*>(fe));
+   const Poly_1D::Basis &basis0d = tfe->GetBasis1D();
+
+   const int quads0d = 1;
+   const int dofs = fe->GetOrder() + 1;
+
+   // We use Matrix and not Vector because we don't want shape0d and dshape0d to have
+   // a different treatment than shape1d and dshape1d
+   // shape0d  = Tensor(dofs, quads0d);
+   // dshape0d = Tensor(dofs, quads0d);
+
+   Vector u(dofs);
+   Vector d(dofs);
+   basis0d.Eval(x, u, d);
+   for (int i = 0; i < dofs; i++)
+   {
+      shape0d(i, 0) = u(i);
+      dshape0d(i, 0) = d(i);
+   }
+}
+
+/**
+* Gives the evaluation of the 1d basis functions and their derivative at all quadrature points
+*/
+template <typename Tensor>
+void ComputeBasis1d(const FiniteElement *fe, int order, Tensor& shape1d,
+                           Tensor& dshape1d, bool backward=false)
+{
+   const TensorBasisElement* tfe(dynamic_cast<const TensorBasisElement*>(fe));
+   const Poly_1D::Basis &basis1d = tfe->GetBasis1D();
+   const IntegrationRule &ir1d = IntRules.Get(Geometry::SEGMENT, order);
+
+   const int quads1d = ir1d.GetNPoints();
+   const int dofs = fe->GetOrder() + 1;
+
+   // shape1d  = Tensor(dofs, quads1d);
+   // dshape1d = Tensor(dofs, quads1d);
+
+   Vector u(dofs);
+   Vector d(dofs);
+   for (int k = 0; k < quads1d; k++)
+   {
+      int ind = backward ? quads1d -1 - k : k;
+      const IntegrationPoint &ip = ir1d.IntPoint(k);
+      basis1d.Eval(ip.x, u, d);
+      for (int i = 0; i < dofs; i++)
+      {
+         shape1d(i, ind) = u(i);
+         dshape1d(i, ind) = d(i);
+      }
+   }
+}
+
+/**
+* Gives the evaluation of the 1d basis functions at one point @param x
+*/
+template <typename Tensor>
+void ComputeBasis0d(const FiniteElement *fe, double x, Tensor& shape0d)
+{
+   const TensorBasisElement* tfe(dynamic_cast<const TensorBasisElement*>(fe));
+   const Poly_1D::Basis &basis0d = tfe->GetBasis1D();
+
+   const int quads0d = 1;
+   const int dofs = fe->GetOrder() + 1;
+
+   // We use Matrix and not Vector because we don't want shape0d and dshape0d to have
+   // a different treatment than shape1d and dshape1d
+   // Well... that was before, we might want to reconsider this.
+   // shape0d  = Tensor(dofs, quads0d);
+
+   Vector u(dofs);
+   Vector d(dofs);
+   basis0d.Eval(x, u, d);
+   for (int i = 0; i < dofs; i++)
+   {
+      shape0d(i, 0) = u(i);
+   }
+}
+
+/**
+* Gives the evaluation of the 1d basis functions at all quadrature points
+*/
+template <typename Tensor>
+void ComputeBasis1d(const FiniteElement *fe, int order, Tensor& shape1d, bool backward=false)
+{
+   const TensorBasisElement* tfe(dynamic_cast<const TensorBasisElement*>(fe));
+   const Poly_1D::Basis &basis1d = tfe->GetBasis1D();
+   const IntegrationRule &ir1d = IntRules.Get(Geometry::SEGMENT, order);
+
+   const int quads1d = ir1d.GetNPoints();
+   const int dofs = fe->GetOrder() + 1;
+
+   // shape1d  = Tensor(dofs, quads1d);
+
+   Vector u(dofs);
+   Vector d(dofs);
+   for (int k = 0; k < quads1d; k++)
+   {
+      int ind = backward ? quads1d -1 - k : k;
+      const IntegrationPoint &ip = ir1d.IntPoint(k);
+      basis1d.Eval(ip.x, u, d);
+      for (int i = 0; i < dofs; i++)
+      {
+         shape1d(i, ind) = u(i);
+      }
+   }
+}
+
+class TensorBasis
+{
+protected:
+   const int dim;
+
+private:
+   typedef Tensor<2> Tensor2d;
+   Tensor2d shape1d;
+
+public:
+   TensorBasis(FiniteElementSpace* fes, const int order)
+   : dim(fes->GetFE(0)->GetDim()),
+     shape1d(fes->GetNDofs1d(),fes->GetNQuads1d(order))
+   {
+      // Store the 1d shape functions and gradients
+      ComputeBasis1d(fes->GetFE(0), order, shape1d);
+   }
+
+   const Tensor2d& getB() const
+   {
+      return shape1d;
+   }
+};
+
+class FaceTensorBasis: public TensorBasis
+{
+private:
+   typedef Tensor<2> Tensor2d;
+   Tensor2d shape0d0, shape0d1;
+
+public:
+   FaceTensorBasis(FiniteElementSpace* fes, const int order)
+   : TensorBasis(fes,order),
+     shape0d0(fes->GetNDofs1d(),fes->GetNQuads1d(order)),
+     shape0d1(fes->GetNDofs1d(),fes->GetNQuads1d(order))
+   {
+      // Store the two 0d shape functions and gradients
+      // in x = 0.0
+      ComputeBasis0d(fes->GetFE(0), 0.0  , shape0d0);
+      // in x = 1.0
+      ComputeBasis0d(fes->GetFE(0), 1.0  , shape0d1);
+   }
+
+   const Tensor2d& getB0d(const int face_id) const
+   {
+      switch(dim)
+      {
+         case 1:
+         switch(face_id)
+         {
+            case 0://WEST
+               return shape0d0;
+            case 1://EAST
+               return shape0d1;
+         }
+         case 2:
+         switch(face_id)
+         {
+            case 0://SOUTH
+               return shape0d0;
+            case 1://EAST
+               return shape0d1;
+            case 2://NORTH
+               return shape0d1;
+            case 3://WEST
+               return shape0d0;
+         }
+         case 3:
+         switch(face_id)
+         {
+            case 0://BOTTOM
+               return shape0d0;
+            case 1://SOUTH
+               return shape0d0;
+            case 2://EAST
+               return shape0d1;
+            case 3://NORTH
+               return shape0d1;
+            case 4://WEST
+               return shape0d0;
+            case 5://TOP
+               return shape0d1;
+         }
+      }
+   }
+
+   const Tensor2d& getB0dTrial(const int face_id) const
+   {
+      switch(dim)
+      {
+         case 1:
+         switch(face_id)
+         {
+            case 0://WEST
+               return shape0d1;
+            case 1://EAST
+               return shape0d0;
+         }
+         case 2:
+         switch(face_id)
+         {
+            case 0://SOUTH
+               return shape0d1;
+            case 1://EAST
+               return shape0d0;
+            case 2://NORTH
+               return shape0d0;
+            case 3://WEST
+               return shape0d1;
+         }
+         case 3:
+         switch(face_id)
+         {
+            case 0://BOTTOM
+               return shape0d1;
+            case 1://SOUTH
+               return shape0d1;
+            case 2://EAST
+               return shape0d0;
+            case 3://NORTH
+               return shape0d0;
+            case 4://WEST
+               return shape0d1;
+            case 5://TOP
+               return shape0d0;
+         }
+      }
+   }
+
+   const Tensor2d& getB0dTest(const int face_id) const
+   {
+      switch(dim)
+      {
+         case 1:
+         switch(face_id)
+         {
+            case 0://WEST
+               return shape0d0;
+            case 1://EAST
+               return shape0d1;
+         }
+         case 2:
+         switch(face_id)
+         {
+            case 0://SOUTH
+               return shape0d0;
+            case 1://EAST
+               return shape0d1;
+            case 2://NORTH
+               return shape0d1;
+            case 3://WEST
+               return shape0d0;
+         }
+         case 3:
+         switch(face_id)
+         {
+            case 0://BOTTOM
+               return shape0d0;
+            case 1://SOUTH
+               return shape0d0;
+            case 2://EAST
+               return shape0d1;
+            case 3://NORTH
+               return shape0d1;
+            case 4://WEST
+               return shape0d0;
+            case 5://TOP
+               return shape0d1;
+         }
+      }
+   }
+};
+
+inline void ScatterDofs(const FiniteElementSpace* mfes, const Table& eldof, const Array<int>& dof_map,
+                  const GridFunction* nodes, const int dofs, const int dim, const int e,
+                  Tensor<2>& LexPointMat)
+{
+   if (dof_map.Size()==0)
+   {
+      if(mfes->GetOrdering()==Ordering::byVDIM)
+      {
+         for (int i = 0; i < dofs; ++i)
+         {
+            for (int j = 0; j < dim; ++j)
+            {
+               LexPointMat(j,i) = (*nodes)( ( e*dofs + i )*dim + j );
+            }
+         }
+      }else{
+         for (int i = 0; i < dofs; ++i)
+         {
+            for (int j = 0; j < dim; ++j)
+            {
+               LexPointMat(j,i) = (*nodes)( ( e*dofs + i ) + j*mfes->GetNDofs() );
+            }
+         }
+      }
+   }else{
+      if(mfes->GetOrdering()==Ordering::byVDIM)
+      {
+         for (int i = 0; i < dofs; ++i)
+         {
+            const int pivot = dof_map[i];
+            for (int j = 0; j < dim; ++j)
+            {
+               LexPointMat(j,i) = (*nodes)(eldof.GetJ()[ e*dofs + pivot ]*dim + j);
+            }
+         }
+      }else{
+         for (int i = 0; i < dofs; ++i)
+         {
+            const int pivot = dof_map[i];
+            for (int j = 0; j < dim; ++j)
+            {
+               LexPointMat(j,i) = (*nodes)(eldof.GetJ()[ e*dofs + pivot ] + j*mfes->GetNDofs());
+            }
+         }
+      }
+   }   
+}
+
+void EvalJacobians1D(const FiniteElementSpace* fes, const int order, Tensor<1>& J)
+{
+   const int dim = 1;
+
+   const Mesh* mesh = fes->GetMesh();
+   const FiniteElementSpace* mfes = mesh->GetNodalFESpace();
+   const Table& eldof = mfes->GetElementToDofTable();
+   const FiniteElement* fe = mfes->GetFE(0);
+   const TensorBasisElement* tfe = dynamic_cast<const TensorBasisElement*>(fe);
+   const Array<int>& dof_map = tfe->GetDofMap();
+   const GridFunction* nodes = mesh->GetNodes();
+   Tensor<2> shape1d(mfes->GetNDofs1d(),mfes->GetNQuads1d(order)), dshape1d(mfes->GetNDofs1d(),mfes->GetNQuads1d(order));
+   ComputeBasis1d( fe, order, shape1d, dshape1d );
+
+   const int NE = fes->GetNE();
+
+   const int quads1d = shape1d.Width();
+   const int dofs1d = shape1d.Height();
+   const int dofs = dofs1d;
+
+   Tensor<2> Jac(J.getData(),quads1d,NE);
+   Jac.zero();
+   Tensor<2> LexPointMat(dim,dofs);
+
+   Tensor<1> T0(LexPointMat.getData(),dofs1d);
+   for (int e = 0; e < NE; ++e)
+   {
+      ScatterDofs(mfes, eldof, dof_map, nodes, dofs, dim, e, LexPointMat);
+      // Computing the Jacobian with the tensor product structure
+      for (int j1 = 0; j1 < quads1d; ++j1)
+      {
+         for (int i1 = 0; i1 < dofs1d; ++i1)
+         {
+            for (int d = 0; d < dim; ++d)
+            {
+               Jac(j1,e) += T0(i1) * dshape1d(i1,j1);
+            }
+         }
+      }
+   }
+}
+
+void EvalJacobians2D(const FiniteElementSpace* fes, const int order, Tensor<1>& J)
+{
+   const int dim = 2;
+
+   const Mesh* mesh = fes->GetMesh();
+   const FiniteElementSpace* mfes = mesh->GetNodalFESpace();
+   const Table& eldof = mfes->GetElementToDofTable();
+   const FiniteElement* fe = mfes->GetFE(0);
+   const TensorBasisElement* tfe = dynamic_cast<const TensorBasisElement*>(fe);
+   const Array<int>& dof_map = tfe->GetDofMap();
+   const GridFunction* nodes = mesh->GetNodes();
+   Tensor<2> shape1d(mfes->GetNDofs1d(),mfes->GetNQuads1d(order)), dshape1d(mfes->GetNDofs1d(),mfes->GetNQuads1d(order));
+   ComputeBasis1d( fe, order, shape1d, dshape1d );
+
+   const int NE = mfes->GetNE();
+
+   const int quads1d = shape1d.Width();
+   const int dofs1d  = shape1d.Height();
+   const int dofs = dofs1d * dofs1d;
+
+
+   Tensor<5> Jac(J.getData(),dim,dim,quads1d,quads1d,NE);
+   Jac.zero();
+   Tensor<2> LexPointMat(dim,dofs);
+
+   Tensor<3> T0(LexPointMat.getData(),dim,dofs1d,dofs1d);
+   Tensor<2> T1b(dim,quads1d), T1d(dim,quads1d);
+   for (int e = 0; e < NE; ++e)
+   {
+      ScatterDofs(mfes, eldof, dof_map, nodes, dofs, dim, e, LexPointMat);
+      // Computing the Jacobian with the tensor product structure
+      for (int i2 = 0; i2 < dofs1d; ++i2)
+      {
+         T1b.zero();
+         T1d.zero();
+         for (int j1 = 0; j1 < quads1d; ++j1)
+         {
+            for (int i1 = 0; i1 < dofs1d; ++i1)
+            {
+               for (int d = 0; d < dim; ++d)
+               {
+                  T1b(d,j1) += T0(d,i1,i2) *  shape1d(i1,j1);
+                  T1d(d,j1) += T0(d,i1,i2) * dshape1d(i1,j1);
+               }
+            }
+         }
+         for (int j2 = 0; j2 < quads1d; ++j2)
+         {
+            for (int j1 = 0; j1 < quads1d; ++j1)
+            {
+               for (int d = 0; d < dim; ++d)
+               {
+                  Jac(d,0,j1,j2,e) += T1d(d,j1) *  shape1d(i2,j2);
+                  Jac(d,1,j1,j2,e) += T1b(d,j1) * dshape1d(i2,j2);
+               }
+            }
+         }
+      }
+   }
+}
+
+void EvalJacobians3D(const FiniteElementSpace* fes, const int order, Tensor<1>& J)
+{
+   const int dim = 3;
+
+   const Mesh* mesh = fes->GetMesh();
+   const FiniteElementSpace* mfes = mesh->GetNodalFESpace();
+   const Table& eldof = mfes->GetElementToDofTable();
+   const FiniteElement* fe = mfes->GetFE(0);
+   const TensorBasisElement* tfe = dynamic_cast<const TensorBasisElement*>(fe);
+   const Array<int>& dof_map = tfe->GetDofMap();
+   const GridFunction* nodes = mesh->GetNodes();
+   Tensor<2> shape1d(mfes->GetNDofs1d(),mfes->GetNQuads1d(order)), dshape1d(mfes->GetNDofs1d(),mfes->GetNQuads1d(order));
+   ComputeBasis1d( fe, order, shape1d, dshape1d );
+
+   const int NE = fes->GetNE();
+
+   const int quads1d = shape1d.Width();
+   const int dofs1d = shape1d.Height();
+   const int dofs = dofs1d * dofs1d * dofs1d;
+
+   Tensor<6> Jac(J.getData(),dim,dim,quads1d,quads1d,quads1d,NE);
+   Jac.zero();
+   Tensor<2> LexPointMat(dim,dofs);
+
+   Tensor<4> T0(LexPointMat.getData(),dim,dofs1d,dofs1d,dofs1d);
+   Tensor<2> T1b(dim,quads1d), T1d(dim,quads1d);
+   Tensor<3> T2bb(dim,quads1d,quads1d), T2db(dim,quads1d,quads1d), T2bd(dim,quads1d,quads1d);
+   for (int e = 0; e < NE; ++e)
+   {
+      ScatterDofs(mfes, eldof, dof_map, nodes, dofs, dim, e, LexPointMat);
+      // Computing the Jacobian with the tensor product structure
+      for (int i3 = 0; i3 < dofs1d; ++i3)
+      {
+         T2bb.zero();
+         T2db.zero();
+         T2bd.zero();
+         for (int i2 = 0; i2 < dofs1d; ++i2)
+         {
+            T1b.zero();
+            T1d.zero();
+            for (int j1 = 0; j1 < quads1d; ++j1)
+            {
+               for (int i1 = 0; i1 < dofs1d; ++i1)
+               {
+                  for (int d = 0; d < dim; ++d)
+                  {
+                     T1b(d,j1) += T0(d,i1,i2,i3) *  shape1d(i1,j1);
+                     T1d(d,j1) += T0(d,i1,i2,i3) * dshape1d(i1,j1);
+                  }
+               }
+            }
+            for (int j2 = 0; j2 < quads1d; ++j2)
+            {
+               for (int j1 = 0; j1 < quads1d; ++j1)
+               {
+                  for (int d = 0; d < dim; ++d)
+                  {
+                     T2bb(d,j1,j2) += T1b(d,j1) *  shape1d(i2,j2);
+                     T2bd(d,j1,j2) += T1b(d,j1) * dshape1d(i2,j2);
+                     T2db(d,j1,j2) += T1d(d,j1) *  shape1d(i2,j2);
+                  }
+               }
+            }
+         }
+         for (int j3 = 0; j3 < quads1d; ++j3)
+         {
+            for (int j2 = 0; j2 < quads1d; ++j2)
+            {
+               for (int j1 = 0; j1 < quads1d; ++j1)
+               {
+                  for (int d = 0; d < dim; ++d)
+                  {
+                     Jac(d,0,j1,j2,j3,e) += T2db(d,j1,j2) *  shape1d(i3,j3);
+                     Jac(d,1,j1,j2,j3,e) += T2bd(d,j1,j2) *  shape1d(i3,j3);
+                     Jac(d,2,j1,j2,j3,e) += T2bb(d,j1,j2) * dshape1d(i3,j3);
+                  }
+               }
+            }
+         }
+      }
+   }
+}
+
+void EvalJacobians( const int dim, const FiniteElementSpace* fes, const int order,
+                     Tensor<1>& J )
+{
+   switch(dim)
+   {
+      case 1:
+         EvalJacobians1D( fes, order, J );
+         break;
+      case 2:
+         EvalJacobians2D( fes, order, J );
+         break;
+      case 3:
+         EvalJacobians3D( fes, order, J );
+         break;
+      default:
+         mfem_error("This orientation does not exist in 3D");
+   }
+}
+
+template <int Dim, typename Op>
+void GetDiag1d(const FiniteElementSpace& fes, const int order, const Op& op, Tensor<Dim>& diag)
+{
+   const int dofs1d = fes.GetNDofs1d();
+   const int quads1d = fes.GetNQuads1d(order);
+   const int nb_elts = fes.GetNE();
+   Tensor<2> diagT(diag, dofs1d, nb_elts);
+   auto Dlin = op.getD();
+   Tensor<2> D(Dlin, quads1d, nb_elts);
+   Tensor<2> shape1d(dofs1d,quads1d);
+   ComputeBasis1d( fes.GetFE(0), order, shape1d);
+
+   diagT.zero();
+   for (int e = 0; e < nb_elts; ++e)
+   {
+      for (int j1 = 0; j1 < quads1d; ++j1)
+      {
+         for (int i1 = 0; i1 < dofs1d; ++i1)
+         {
+               diagT(i1,e) += D(j1,e) * shape1d(i1,j1) * shape1d(i1,j1);
+         }
+      }
+   }
+}
+
+template <int Dim, typename Op>
+void GetDiag2d(const FiniteElementSpace& fes, const int order, const Op& op, Tensor<Dim>& diag)
+{
+   const int dofs1d = fes.GetNDofs1d();
+   const int quads1d = fes.GetNQuads1d(order);
+   const int nb_elts = fes.GetNE();
+   Tensor<3> diagT(diag, dofs1d, dofs1d, nb_elts);
+   auto Dlin = op.getD();
+   Tensor<3> D(Dlin, quads1d, quads1d, nb_elts);
+   Tensor<2> shape1d(dofs1d,quads1d);
+   ComputeBasis1d( fes.GetFE(0), order, shape1d);
+
+   Tensor<2> T1(dofs1d,quads1d);
+
+   diagT.zero();
+   for (int e = 0; e < nb_elts; ++e)
+   {
+      T1.zero();
+      for (int j2 = 0; j2 < quads1d; ++j2)
+      {
+         for (int j1 = 0; j1 < quads1d; ++j1)
+         {
+            for (int i1 = 0; i1 < dofs1d; ++i1)
+            {
+                  T1(i1,j2) += D(j1,j2,e) * shape1d(i1,j1) * shape1d(i1,j1);
+            }
+         }
+      }
+      for (int j2 = 0; j2 < quads1d; ++j2)
+      {
+         for (int i2 = 0; i2 < dofs1d; ++i2)
+         {
+            for (int i1 = 0; i1 < dofs1d; ++i1)
+            {
+                  diagT(i1,i2,e) += T1(i1,j2) * shape1d(i2,j2) * shape1d(i2,j2);
+            }
+         }
+      }
+   }
+}
+
+template <int Dim, typename Op>
+void GetDiag3d(const FiniteElementSpace& fes, const int order, const Op& op, Tensor<Dim>& diag)
+{
+   const int dofs1d = fes.GetNDofs1d();
+   const int quads1d = fes.GetNQuads1d(order);
+   const int nb_elts = fes.GetNE();
+   Tensor<4> diagT(diag, dofs1d, dofs1d, dofs1d, nb_elts);
+   auto Dlin = op.getD();
+   Tensor<4> D(Dlin, quads1d, quads1d, quads1d, nb_elts);
+   Tensor<2> shape1d(dofs1d,quads1d);
+   ComputeBasis1d( fes.GetFE(0), order, shape1d);
+
+   Tensor<3> T1(dofs1d,quads1d,quads1d),T2(dofs1d,dofs1d,quads1d);
+
+   diagT.zero();
+   for (int e = 0; e < nb_elts; ++e)
+   {
+      T1.zero();
+      for (int j3 = 0; j3 < quads1d; ++j3)
+      {
+         for (int j2 = 0; j2 < quads1d; ++j2)
+         {
+            for (int j1 = 0; j1 < quads1d; ++j1)
+            {
+               for (int i1 = 0; i1 < dofs1d; ++i1)
+               {
+                     T1(i1,j2,j3) += D(j1,j2,j3,e) * shape1d(i1,j1) * shape1d(i1,j1);
+               }
+            }
+         }
+      }
+      T2.zero();
+      for (int j3 = 0; j3 < quads1d; ++j3)
+      {
+         for (int j2 = 0; j2 < quads1d; ++j2)
+         {
+            for (int i2 = 0; i2 < dofs1d; ++i2)
+            {
+               for (int i1 = 0; i1 < dofs1d; ++i1)
+               {
+                     T2(i1,i2,j3) += T1(i1,j2,j3) * shape1d(i2,j2) * shape1d(i2,j2);
+               }
+            }
+         }
+      }
+      for (int j3 = 0; j3 < quads1d; ++j3)
+      {
+         for (int i3 = 0; i3 < dofs1d; ++i3)
+         {
+            for (int i2 = 0; i2 < dofs1d; ++i2)
+            {
+               for (int i1 = 0; i1 < dofs1d; ++i1)
+               {
+                     diagT(i1,i2,i3,e) += T2(i1,i2,j3) * shape1d(i3,j3) * shape1d(i3,j3);
+               }
+            }
+         }
+      }
+   }
+}
+
+template <int Dim, typename Op>
+void GetDiag(const FiniteElementSpace& fes, const int order, const Op& op, Tensor<Dim>& diag)
+{
+   switch(fes.GetFE(0)->GetDim())
+   {
+      case 1:
+         GetDiag1d(fes, order, op, diag);
+         break;
+      case 2:
+         GetDiag2d(fes, order, op, diag);
+         break;
+      case 3:
+         GetDiag3d(fes, order, op, diag);
+         break;
+   }
+}
+
+}
+
+#endif

linalg/densemat.hpp

@@ -337,6 +337,19 @@ public:
 
    /// Destroys dense matrix.
    virtual ~DenseMatrix();
+
+   friend std::ostream& operator<<(std::ostream& os, const DenseMatrix& M)
+   {
+      for (int i = 0; i < M.Height(); ++i)
+      {
+        for (int j = 0; j < M.Width(); ++j)
+        {
+          os << M(i,j) << " ";   
+        }
+        os << "\n";
+      }
+      return os;
+   }
 };
 
 /// C = A + alpha*B
@@ -656,6 +669,15 @@ public:
       own_data = true;
    }
 
+   DenseTensor(double *d, int i, int j, int k)
+      : Mk(NULL, i, j)
+   {
+      nk = k;
+      tdata = d;
+      own_data = false;
+
+   }
+   
    /// Copy constructor: deep copy
    DenseTensor(const DenseTensor& other)
       : Mk(NULL, other.Mk.height, other.Mk.width), nk(other.nk), own_data(true)

linalg/operator.hpp

@@ -125,6 +125,7 @@ public:
    {
       ANY_TYPE,         ///< ID for the base class Operator, i.e. any type.
       MFEM_SPARSEMAT,   ///< ID for class SparseMatrix
+      MFEM_FORMOPER,
       Hypre_ParCSR,     ///< ID for class HypreParMatrix.
       PETSC_MATAIJ,     ///< ID for class PetscParMatrix, MATAIJ format.
       PETSC_MATIS,      ///< ID for class PetscParMatrix, MATIS format.

mesh/mesh.cpp

@@ -764,13 +764,13 @@ FaceElementTransformations *Mesh::GetBdrFaceTransformations(int BdrElemNo)
    return tr;
 }
 
-void Mesh::GetFaceElements(int Face, int *Elem1, int *Elem2)
+void Mesh::GetFaceElements(int Face, int *Elem1, int *Elem2) const
 {
    *Elem1 = faces_info[Face].Elem1No;
    *Elem2 = faces_info[Face].Elem2No;
 }
 
-void Mesh::GetFaceInfos(int Face, int *Inf1, int *Inf2)
+void Mesh::GetFaceInfos(int Face, int *Inf1, int *Inf2) const
 {
    *Inf1 = faces_info[Face].Elem1Inf;
    *Inf2 = faces_info[Face].Elem2Inf;

mesh/mesh.hpp

@@ -844,8 +844,8 @@ public:
    {
       return (faces_info[FaceNo].Elem2No >= 0);
    }
-   void GetFaceElements (int Face, int *Elem1, int *Elem2);
-   void GetFaceInfos (int Face, int *Inf1, int *Inf2);
+   void GetFaceElements (int Face, int *Elem1, int *Elem2) const;
+   void GetFaceInfos (int Face, int *Inf1, int *Inf2) const;
 
    int GetFaceGeometryType(int Face) const;
    int GetFaceElementType(int Face) const;