Merge branch 'master' into lp-error-est-dev

Author: tzanio

CHANGELOG

@@ -35,6 +35,12 @@ Discretization improvements
 
 - Added support for simplices in GSLIB-FindPoints.
 
+- Added support for H1 and L2 element matrix assembly in the mass, convection,
+  diffusion, transpose, and the face DG trace integrators. This is compatible
+  with GPU device execution and is illustrated in Example 9/9p, see the option
+  '-ea'. When enabled, this level of assembly stores independent dense matrices
+  for the elements, and independent dense matrices for the faces in the DG case.
+
 - Added new partial assembly kernels for H(div) bilinear forms, as well as
   VectorFEDivergenceIntegrator.
 

examples/ex9.cpp

@@ -19,6 +19,7 @@
 //
 // Device sample runs:
 //    ex9 -pa
+//    ex9 -ea
 //    ex9 -pa -m ../data/periodic-cube.mesh
 //    ex9 -pa -m ../data/periodic-cube.mesh -d cuda
 //
@@ -142,6 +143,7 @@ int main(int argc, char *argv[])
    int ref_levels = 2;
    int order = 3;
    bool pa = false;
+   bool ea = false;
    const char *device_config = "cpu";
    int ode_solver_type = 4;
    double t_final = 10.0;
@@ -166,6 +168,8 @@ int main(int argc, char *argv[])
                   "Order (degree) of the finite elements.");
    args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
                   "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&ea, "-ea", "--element-assembly", "-no-ea",
+                  "--no-element-assembly", "Enable Element Assembly.");
    args.AddOption(&device_config, "-d", "--device",
                   "Device configuration string, see Device::Configure().");
    args.AddOption(&ode_solver_type, "-s", "--ode-solver",
@@ -269,6 +273,11 @@ int main(int argc, char *argv[])
       m.SetAssemblyLevel(AssemblyLevel::PARTIAL);
       k.SetAssemblyLevel(AssemblyLevel::PARTIAL);
    }
+   else if (ea)
+   {
+      m.SetAssemblyLevel(AssemblyLevel::ELEMENT);
+      k.SetAssemblyLevel(AssemblyLevel::ELEMENT);
+   }
    m.AddDomainIntegrator(new MassIntegrator);
    k.AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
    k.AddInteriorFaceIntegrator(
@@ -429,8 +438,9 @@ FE_Evolution::FE_Evolution(BilinearForm &_M, BilinearForm &_K, const Vector &_b)
    : TimeDependentOperator(_M.Height()), M(_M), K(_K), b(_b), z(_M.Height())
 {
    bool pa = M.GetAssemblyLevel() == AssemblyLevel::PARTIAL;
+   bool ea = M.GetAssemblyLevel() == AssemblyLevel::ELEMENT;
    Array<int> ess_tdof_list;
-   if (pa)
+   if (pa || ea)
    {
       M_prec = new OperatorJacobiSmoother(M, ess_tdof_list);
       M_solver.SetOperator(M);

examples/ex9p.cpp

@@ -19,6 +19,7 @@
 //
 // Device sample runs:
 //    mpirun -np 4 ex9p -pa
+//    mpirun -np 4 ex9p -ea
 //    mpirun -np 4 ex9p -pa -m ../data/periodic-cube.mesh
 //    mpirun -np 4 ex9p -pa -m ../data/periodic-cube.mesh -d cuda
 //
@@ -161,6 +162,7 @@ int main(int argc, char *argv[])
    int par_ref_levels = 0;
    int order = 3;
    bool pa = false;
+   bool ea = false;
    const char *device_config = "cpu";
    int ode_solver_type = 4;
    double t_final = 10.0;
@@ -188,6 +190,8 @@ int main(int argc, char *argv[])
                   "Order (degree) of the finite elements.");
    args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
                   "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&ea, "-ea", "--element-assembly", "-no-ea",
+                  "--no-element-assembly", "Enable Element Assembly.");
    args.AddOption(&device_config, "-d", "--device",
                   "Device configuration string, see Device::Configure().");
    args.AddOption(&ode_solver_type, "-s", "--ode-solver",
@@ -319,6 +323,11 @@ int main(int argc, char *argv[])
       m->SetAssemblyLevel(AssemblyLevel::PARTIAL);
       k->SetAssemblyLevel(AssemblyLevel::PARTIAL);
    }
+   else if (ea)
+   {
+      m->SetAssemblyLevel(AssemblyLevel::ELEMENT);
+      k->SetAssemblyLevel(AssemblyLevel::ELEMENT);
+   }
    m->AddDomainIntegrator(new MassIntegrator);
    k->AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
    k->AddInteriorFaceIntegrator(
@@ -556,8 +565,9 @@ FE_Evolution::FE_Evolution(ParBilinearForm &_M, ParBilinearForm &_K,
      z(_M.Height())
 {
    bool pa = _M.GetAssemblyLevel()==AssemblyLevel::PARTIAL;
+   bool ea = _M.GetAssemblyLevel()==AssemblyLevel::ELEMENT;
 
-   if (pa)
+   if (pa || ea)
    {
       M.Reset(&_M, false);
       K.Reset(&_K, false);
@@ -571,7 +581,7 @@ FE_Evolution::FE_Evolution(ParBilinearForm &_M, ParBilinearForm &_K,
    M_solver.SetOperator(*M);
 
    Array<int> ess_tdof_list;
-   if (pa)
+   if (pa || ea)
    {
       M_prec = new OperatorJacobiSmoother(_M, ess_tdof_list);
       dg_solver = NULL;

fem/CMakeLists.txt

@@ -13,15 +13,20 @@ set(SRCS
   bilinearform.cpp
   bilinearform_ext.cpp
   bilininteg.cpp
-  bilininteg_convection.cpp
-  bilininteg_dgtrace.cpp
-  bilininteg_diffusion.cpp
+  bilininteg_convection_pa.cpp
+  bilininteg_convection_ea.cpp
+  bilininteg_dgtrace_pa.cpp
+  bilininteg_dgtrace_ea.cpp
+  bilininteg_diffusion_pa.cpp
+  bilininteg_diffusion_ea.cpp
   bilininteg_divergence.cpp
   bilininteg_hcurl.cpp
   bilininteg_hdiv.cpp
   bilininteg_vectorfe.cpp
   bilininteg_gradient.cpp
-  bilininteg_mass.cpp
+  bilininteg_mass_pa.cpp
+  bilininteg_mass_ea.cpp
+  bilininteg_transpose_ea.cpp
   bilininteg_vecdiffusion.cpp
   bilininteg_vecmass.cpp
   coefficient.cpp

fem/bilinearform.cpp

@@ -126,8 +126,7 @@ void BilinearForm::SetAssemblyLevel(AssemblyLevel assembly_level)
          // Use the original BilinearForm implementation for now
          break;
       case AssemblyLevel::ELEMENT:
-         mfem_error("Element assembly not supported yet... stay tuned!");
-         // ext = new EABilinearFormExtension(this);
+         ext = new EABilinearFormExtension(this);
          break;
       case AssemblyLevel::PARTIAL:
          ext = new PABilinearFormExtension(this);