added a check to verify h1-l2 mixed projection should work

examples/ex1.cpp

@@ -119,8 +119,8 @@ int main(int argc, char *argv[])
    //    largest number that gives a final mesh with no more than 50,000
    //    elements.
    {
-      int ref_levels =
-         (int)floor(log(50000./mesh.GetNE())/log(2.)/dim);
+     int ref_levels = 
+       (int)floor(log(50000./mesh.GetNE())/log(2.)/dim);
       for (int l = 0; l < ref_levels; l++)
       {
          mesh.UniformRefinement();
@@ -191,7 +191,83 @@ int main(int argc, char *argv[])
    //     static condensation, etc.
    if (static_cond) { a.EnableStaticCondensation(); }
    a.Assemble();
+   
+   /// Test code
+   a.Finalize(0);
+   SparseMatrix A_mat(a.SpMat());
 
+   AssembledSparseMatrix A_assembled_mat(fespace, fespace, ElementDofOrdering::NATIVE);
+   Vector A_elem_mat_data(fespace.GetFE(0)->GetDof() * fespace.GetFE(0)->GetDim() *
+			  fespace.GetFE(0)->GetDof() * fespace.GetFE(0)->GetDim() *
+			  fespace.GetNE());
+   A_elem_mat_data = 0.;
+   DenseTensor A_elem_mat;
+   A_elem_mat.UseExternalData(A_elem_mat_data.GetData(),
+			      fespace.GetFE(0)->GetDof() * fespace.GetFE(0)->GetDim(),
+			      fespace.GetFE(0)->GetDof() * fespace.GetFE(0)->GetDim(),
+			      fespace.GetNE());
+   for (int e = 0; e < fespace.GetNE(); e++) {
+     a.ComputeElementMatrix(e, A_elem_mat(e));
+   }
+   A_assembled_mat.FillData(A_elem_mat_data);
+   A_assembled_mat.Finalize();
+   
+   double diff = 0.;
+   for (int r = 0; r < A_assembled_mat.Height(); r++) {
+     auto columns = A_assembled_mat.GetRowColumns(r);
+     for (int c = 0; c < A_assembled_mat.RowSize(r); c++) {
+       diff += (A_mat(r, columns[c]) - A_assembled_mat(r, columns[c])) *
+	 (A_mat(r, columns[c]) - A_assembled_mat(r, columns[c]));
+       if (fabs(A_mat(r, columns[c]) - A_assembled_mat(r, columns[c])) > 1.e-5) {
+	 std::cout << "(" << r << "," << c << ") " << A_mat(r, columns[c]) << " " << A_assembled_mat(r, columns[c]) << " " << (A_mat(r, columns[c]) - A_assembled_mat(r, columns[c])) << std::endl;
+       }
+     }
+   }
+   std::cout << "AssembledMatrix Diff Norm: " << sqrt(diff) << std::endl;
+
+   // test H1 to L2
+
+   {
+     L2_FECollection    fec_l2(order, dim);
+     FiniteElementSpace fespace_l2(&mesh, &fec_l2);
+     MixedBilinearForm blf(&fespace, &fespace_l2);
+     blf.AddDomainIntegrator(new MixedScalarMassIntegrator());
+     blf.Assemble(0);
+     blf.Finalize(0);
+
+     SparseMatrix B_mat(blf.SpMat());
+     AssembledSparseMatrix B_assembled_mat(fespace_l2, fespace, ElementDofOrdering::NATIVE);
+     Vector B_elem_mat_data(fespace_l2.GetFE(0)->GetDof() * fespace_l2.GetFE(0)->GetDim() *
+			    fespace.GetFE(0)->GetDof() * fespace.GetFE(0)->GetDim() *
+			    fespace.GetNE());
+     B_elem_mat_data = 0.;
+     DenseTensor B_elem_mat;
+     B_elem_mat.UseExternalData(B_elem_mat_data.GetData(),
+				fespace_l2.GetFE(0)->GetDof() * fespace_l2.GetFE(0)->GetDim(),
+				fespace.GetFE(0)->GetDof() * fespace.GetFE(0)->GetDim(),
+				fespace.GetNE());
+     for (int e = 0; e < fespace.GetNE(); e++) {
+       blf.ComputeElementMatrix(e, B_elem_mat(e));
+     }
+     B_assembled_mat.FillData(B_elem_mat_data);
+     B_assembled_mat.Finalize();
+
+     double diff = 0.;
+     for (int r = 0; r < B_assembled_mat.Height(); r++) {
+       auto columns = B_assembled_mat.GetRowColumns(r);
+       for (int c = 0; c < B_assembled_mat.RowSize(r); c++) {
+	 diff += (B_mat(r, columns[c]) - B_assembled_mat(r, columns[c])) *
+	   (B_mat(r, columns[c]) - B_assembled_mat(r, columns[c]));
+	 if (fabs(B_mat(r, columns[c]) - B_assembled_mat(r, columns[c])) > 1.e-5) {
+	   std::cout << "(" << r << "," << c << ") " << B_mat(r, columns[c]) << " " << B_assembled_mat(r, columns[c]) << " " << (B_mat(r, columns[c]) - B_assembled_mat(r, columns[c])) << std::endl;
+	 }
+       }
+     }
+     std::cout << "AssembledMatrix H1-L2 Diff Norm: " << sqrt(diff) << std::endl;   
+   }
+   
+   /// End test code
+   
    OperatorPtr A;
    Vector B, X;
    a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);

examples/ex3.cpp

@@ -106,8 +106,8 @@ int main(int argc, char *argv[])
    //    largest number that gives a final mesh with no more than 50,000
    //    elements.
    {
-      int ref_levels =
-         (int)floor(log(50000./mesh->GetNE())/log(2.)/dim);
+     int ref_levels = 
+       (int)floor(log(50000./mesh->GetNE())/log(2.)/dim);
       for (int l = 0; l < ref_levels; l++)
       {
          mesh->UniformRefinement();
@@ -169,6 +169,43 @@ int main(int argc, char *argv[])
    if (static_cond) { a->EnableStaticCondensation(); }
    a->Assemble();
 
+   /// Test code
+   a->Finalize(0);
+   SparseMatrix A_mat(a->SpMat());
+
+   AssembledSparseMatrix A_assembled_mat(*fespace, *fespace, ElementDofOrdering::NATIVE);
+   Vector A_elem_mat_data(fespace->GetFE(0)->GetDof() * fespace->GetFE(0)->GetDim() *
+			  fespace->GetFE(0)->GetDof() * fespace->GetFE(0)->GetDim() *
+			  fespace->GetNE());
+   A_elem_mat_data = 0.;
+   DenseTensor A_elem_mat;
+   A_elem_mat.UseExternalData(A_elem_mat_data.GetData(),
+			      fespace->GetFE(0)->GetDof() * fespace->GetFE(0)->GetDim(),
+			      fespace->GetFE(0)->GetDof() * fespace->GetFE(0)->GetDim(),
+			      fespace->GetNE());
+   for (int e = 0; e < fespace->GetNE(); e++) {
+     a->ComputeElementMatrix(e, A_elem_mat(e));
+   }
+   A_assembled_mat.FillData(A_elem_mat_data);
+   A_assembled_mat.Finalize();
+   
+   double diff = 0.;
+   for (int r = 0; r < A_assembled_mat.Height(); r++) {
+     auto columns = A_assembled_mat.GetRowColumns(r);
+     for (int c = 0; c < A_assembled_mat.RowSize(r); c++) {
+       const int column_ind = columns[c]  >= 0 ? columns[c] : -1-columns[c];
+       diff += (A_mat(r, column_ind) - A_assembled_mat(r, column_ind)) *
+	 (A_mat(r, column_ind) - A_assembled_mat(r, column_ind));
+       if (fabs(A_mat(r, column_ind) - A_assembled_mat(r, column_ind)) > 1.e-5) {
+	 std::cout << "(" << r << "," << column_ind << ") " << A_mat(r, column_ind) << " " << A_assembled_mat(r, column_ind) << " " << (A_mat(r, column_ind) - A_assembled_mat(r, column_ind)) << std::endl;
+       }
+     }
+   }
+   std::cout << "AssembledMatrix Diff Norm: " << sqrt(diff) << std::endl;
+	       
+   /// End test code
+
+   
    OperatorPtr A;
    Vector B, X;
    a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);

fem/restriction.cpp

@@ -1447,4 +1447,254 @@ int ToLexOrdering(const int dim, const int face_id, const int size1d,
    }
 }
 
+AssembledSparseMatrix::AssembledSparseMatrix(const mfem::FiniteElementSpace& test,   // test_elem_dofs * ne * vdim x vdim * test_ndofs
+                        const mfem::FiniteElementSpace& trial,  // trial_elem_dofs * ne * vdim x vdim * trial_ndofs
+                        mfem::ElementDofOrdering        elem_order)
+      : mfem::SparseMatrix(test.GetNDofs() * test.GetVDim(), trial.GetNDofs() * trial.GetVDim()),
+        test_fes(test),
+        trial_fes(trial),
+        test_restriction(test, elem_order),
+        trial_restriction(trial, elem_order),
+        elem_ordering(elem_order)
+  {
+    GetMemoryI().New(Height() + 1, GetMemoryI().GetMemoryType());
+
+    const int nnz = FillI();
+    GetMemoryJ().New(nnz, GetMemoryJ().GetMemoryType());
+    GetMemoryData().New(nnz, GetMemoryData().GetMemoryType());
+    FillJ();
+
+    // zero initialize the data
+    for (int i = 0; i < nnz; i++) {
+      A[i] = 0.;
+    }
+  }
+
+int AssembledSparseMatrix::FillI()
+{
+  // ElementRestriction creates a CSR matrix that maps vdof -> (dof, ne).
+  // offsets are the row offsets corresponding to a vdof
+  // indices maps a given vdof to the the assembled element matrix vector (dof * ne + d).
+  // gatherMap takes an element matrix vector offset (dof, ne) and returns the partition-local vdof (d.o.f. id).
+  [[maybe_unused]] auto& test_offsets = test_restriction.offsets;  
+  [[maybe_unused]] auto& test_indices = test_restriction.indices;    
+  [[maybe_unused]] auto& test_gatherMap  = test_restriction.gatherMap;
+  [[maybe_unused]] auto& trial_offsets   = trial_restriction.offsets;
+  [[maybe_unused]] auto& trial_indices   = trial_restriction.indices;
+  [[maybe_unused]] auto& trial_gatherMap = trial_restriction.gatherMap;
+
+  /**
+     We expect mat_ea to be of size (test_elem_dof * test_vdim, trial_elem_dof * trial_vdim, ne)
+     We assume a consistent striding from (elem_dof, vd) within each element
+
+   */
+  const int test_elem_dof  = test_fes.GetFE(0)->GetDof();
+  const int trial_elem_dof = trial_fes.GetFE(0)->GetDof();
+  const int test_vdim      = test_fes.GetVDim();
+  const int trial_vdim     = trial_fes.GetVDim();
+  const int test_ndofs     = test_fes.GetNDofs();
+
+  auto I = ReadWriteI();
+  for (int i = 0; i < test_vdim * test_ndofs; i++) {
+    I[i] = 0;
+  }
+
+  for (int test_vdof = 0; test_vdof < test_fes.GetNDofs(); test_vdof++) {
+    // Look through each element corresponding to a test_vdof
+    const int test_row_offset = test_offsets[test_vdof];
+    const int nrow_elems      = test_offsets[test_vdof + 1] - test_row_offset;
+
+    // Build temporary array to get rid of duplicates
+    mfem::Array<int> trial_vdofs(nrow_elems * trial_elem_dof);
+    trial_vdofs = -1;
+    int nnz_row = 0;
+    for (int e_index = 0; e_index < nrow_elems; e_index++) {
+      // test_indices can be negative in the case of Hcurl
+      const int                  test_index_v = test_indices[test_row_offset + e_index];
+      const int                  test_index = test_index_v >= 0 ? test_index_v : -test_index_v - 1;
+      const int                  e           = test_index / test_elem_dof;
+      [[maybe_unused]] const int test_i_elem = test_index % test_elem_dof;
+
+      // find corresponding trial_vdofs
+      mfem::Array<int> trial_elem_vdofs(trial_elem_dof);
+      for (int j = 0; j < trial_elem_dof; j++) {
+	// this might be negative
+        const auto trial_j_vdof_v = trial_gatherMap[trial_elem_dof * e + j];
+	const auto trial_j_vdof = trial_j_vdof_v >= 0 ? trial_j_vdof_v : -1 -trial_j_vdof_v;
+        trial_elem_vdofs[j]     = trial_j_vdof;
+        if (trial_vdofs.Find(trial_j_vdof) == -1) {
+          // we haven't seen this before
+          trial_vdofs[nnz_row] = trial_j_vdof;
+          nnz_row++;
+        }
+      }
+    }
+
+    // add entries to I
+    for (int vi = 0; vi < test_vdim; vi++) {
+      const auto nnz_index_v = test_fes.DofToVDof(test_vdof, vi);
+      const auto nnz_index = nnz_index_v >= 0 ? nnz_index_v : -1 -nnz_index_v;
+      I[nnz_index] = nnz_row * trial_vdim;
+    }
+  }
+
+  // Perform inclusive scan on all entries
+  int nnz = 0;
+  for (int i = 0; i < test_ndofs * trial_vdim; i++) {
+    int temp = I[i];
+    I[i]     = nnz;
+    nnz += temp;
+  }
+  I[test_ndofs * trial_vdim] = nnz;
+
+  return nnz;
+}
+
+void AssembledSparseMatrix::FillJ()
+{
+  auto I = ReadWriteI();
+  auto J = WriteJ();
+
+  [[maybe_unused]] auto& test_offsets = test_restriction.offsets;  // offsets for rows.. each row is a test_vdof
+  [[maybe_unused]] auto& test_indices = test_restriction.indices;  // returns (test_elem_dof , ne) id corresponding to a test_vdof_offset
+  [[maybe_unused]] auto& test_gatherMap  = test_restriction.gatherMap;  // returns test_vdof
+  [[maybe_unused]] auto& trial_offsets   = trial_restriction.offsets;
+  [[maybe_unused]] auto& trial_indices   = trial_restriction.indices;
+  [[maybe_unused]] auto& trial_gatherMap = trial_restriction.gatherMap;
+
+  const int                  test_elem_dof  = test_fes.GetFE(0)->GetDof();
+  const int                  trial_elem_dof = trial_fes.GetFE(0)->GetDof();
+  const int                  test_vdim      = test_fes.GetVDim();
+  const int                  trial_vdim     = trial_fes.GetVDim();
+  [[maybe_unused]] const int test_ndofs     = test_fes.GetNDofs();
+
+  const int ne = trial_fes.GetNE();
+  ea_map.SetSize(test_elem_dof * test_vdim * trial_elem_dof * trial_vdim * ne);
+  auto map_ea = Reshape(ea_map.ReadWrite(), test_elem_dof * test_vdim, trial_elem_dof * trial_vdim, ne);
+
+  // initialize J
+  for (int j = 0; j < this->J.Capacity(); j++) {
+    this->J[j] = -1;
+  }
+
+  for (int test_vdof = 0; test_vdof < test_fes.GetNDofs(); test_vdof++) {
+    // Look through each element corresponding to a test_vdof
+    const int test_row_offset = test_offsets[test_vdof];
+    const int nrow_elems      = test_offsets[test_vdof + 1] - test_row_offset;
+
+    // here we assume all the components have the same number of columns
+    const int        nnz_row = I[test_fes.DofToVDof(test_vdof, 0) + 1] - I[test_fes.DofToVDof(test_vdof, 0)];
+    mfem::Array<int> trial_vdofs(nnz_row);
+    trial_vdofs      = -1;
+    int j_vdof_index = 0;
+
+    // Build temporary array for assembled J
+    for (int e_index = 0; e_index < nrow_elems; e_index++) {
+            // test_indices can be negative in the case of Hcurl
+      const int                  test_index_v = test_indices[test_row_offset + e_index];
+      const int                  test_index = test_index_v >= 0 ? test_index_v : -test_index_v - 1;
+      const int                  e           = test_index / test_elem_dof;
+      [[maybe_unused]] const int test_i_elem = test_index % test_elem_dof;
+
+      // find corresponding trial_vdofs
+      mfem::Array<int> trial_elem_vdofs(trial_elem_dof);
+      for (int j_elem = 0; j_elem < trial_elem_dof; j_elem++) {
+	// could be negative.. but trial_elem_vdofs is a temporary array
+        const auto trial_j_vdof_v  = trial_gatherMap[trial_elem_dof * e + j_elem];
+	const auto trial_j_vdof = trial_j_vdof_v >= 0 ? trial_j_vdof_v : -1 -trial_j_vdof_v;
+        trial_elem_vdofs[j_elem] = trial_j_vdof;
+
+	// since trial_j_vdof could be negative but there are now two indices that point to the same dof (just oriented differently).. we only want to search for positive oriented indices
+        auto find_index = trial_vdofs.Find(trial_j_vdof);
+        if (find_index == -1) {
+          // we haven't seen this before
+          trial_vdofs[j_vdof_index] = trial_j_vdof;
+
+          // we can add this entry to J
+          for (int vi = 0; vi < test_vdim; vi++) {
+            const auto i_dof_offset = I[test_fes.DofToVDof(test_vdof, vi)];
+
+            // this access pattern corresnponds to j_vdof_index + vj * nnz_row
+            for (int vj = 0; vj < trial_vdim; vj++) {
+              const auto column_index = j_vdof_index + vj * nnz_row / trial_vdim;
+              const auto j_nnz_index  = i_dof_offset + column_index;
+	      // this index may be negative, but J needs to be positive
+	      const auto j_value = trial_fes.DofToVDof(trial_j_vdof_v, vj);
+              J[j_nnz_index]          = j_value >= 0 ? j_value : -1-j_value;
+            }
+          }
+
+          // write mapping from ea to csr_nnz_index (can probably optimize this)
+          for (int vi = 0; vi < test_vdim; vi++) {
+            const auto i_dof_offset = I[test_fes.DofToVDof(test_vdof, vi)];
+            for (int vj = 0; vj < trial_vdim; vj++) {
+              const auto column_index = j_vdof_index + vj * nnz_row / trial_vdim;
+	      const int index_val = i_dof_offset + column_index;
+	      const int trial_index = trial_fes.DofToVDof(trial_j_vdof_v, vj);
+	      const int orientation_factor = (test_index_v >= 0 ? 1 : -1) * (trial_index >=0 ? 1 : -1);
+              map_ea(test_i_elem + test_elem_dof * vi, j_elem + trial_elem_dof * vj, e) =
+		orientation_factor > 0 ? index_val : -1-index_val;
+            }
+          }
+
+          j_vdof_index++;
+        } else {
+          // this is a duplicate entry
+          // write mapping from ea to csr_nnz_index (can probably optimize this)
+          for (int vi = 0; vi < test_vdim; vi++) {
+            const auto i_dof_offset = I[test_fes.DofToVDof(test_vdof, vi)];
+            for (int vj = 0; vj < trial_vdim; vj++) {
+              const auto column_index = find_index + vj * nnz_row / trial_vdim;
+	      const int index_val = i_dof_offset + column_index;
+	      const int trial_index = trial_fes.DofToVDof(trial_j_vdof_v, vj);
+	      const int orientation_factor = (test_index_v >= 0 ? 1 : -1) * (trial_index >=0 ? 1 : -1);      
+              map_ea(test_i_elem + test_elem_dof * vi, j_elem + trial_elem_dof * vj, e) =
+		orientation_factor > 0 ? index_val : -1-index_val;
+            }
+          }
+        }
+      }
+    }
+  }
+}
+void AssembledSparseMatrix::FillData(const mfem::Vector& ea_data)
+{
+  auto Data = WriteData();
+
+  [[maybe_unused]] auto& test_offsets = test_restriction.offsets;  // offsets for rows.. each row is a test_vdof
+  [[maybe_unused]] auto& test_indices = test_restriction.indices;  // returns (test_elem_dof , ne) id corresponding to a test_vdof_offset
+  [[maybe_unused]] auto& test_gatherMap  = test_restriction.gatherMap;  // returns test_vdof
+  [[maybe_unused]] auto& trial_offsets   = trial_restriction.offsets;
+  [[maybe_unused]] auto& trial_indices   = trial_restriction.indices;
+  [[maybe_unused]] auto& trial_gatherMap = trial_restriction.gatherMap;
+  
+  const int                  test_elem_dof  = test_fes.GetFE(0)->GetDof();
+  const int                  trial_elem_dof = trial_fes.GetFE(0)->GetDof();
+  const int                  test_vdim      = test_fes.GetVDim();
+  const int                  trial_vdim     = trial_fes.GetVDim();
+  [[maybe_unused]] const int test_ndofs     = test_fes.GetNDofs();
+
+  const int ne = trial_fes.GetNE();
+
+  auto map_ea = Reshape(ea_map.Read(), test_elem_dof * test_vdim, trial_elem_dof * trial_vdim, ne);
+
+  auto mat_ea = Reshape(ea_data.Read(), test_elem_dof * test_vdim, trial_elem_dof * trial_vdim, ne);
+
+  // Use map_ea to take ea_data directly to CSR entry
+  for (int e = 0; e < ne; e++) {
+    for (int i_elem = 0; i_elem < test_elem_dof; i_elem++) {
+      for (int vi = 0; vi < test_vdim; vi++) {
+        for (int j_elem = 0; j_elem < trial_elem_dof; j_elem++) {
+          for (int vj = 0; vj < trial_vdim; vj++) {
+	    const auto map_ea_v = map_ea(i_elem + vi * test_elem_dof, j_elem + vj * trial_elem_dof, e);
+	    const auto map_ea_index = map_ea_v >= 0 ? map_ea_v : -1 -map_ea_v;
+            Data[map_ea_index] += (map_ea_v >= 0 ? 1 : -1) * 
+                mat_ea(i_elem + vi * test_elem_dof, j_elem + vj * trial_elem_dof, e);
+          }
+        }
+      }
+    }
+  }
+}
+
 } // namespace mfem

fem/restriction.hpp

@@ -19,6 +19,7 @@ namespace mfem
 {
 
 class FiniteElementSpace;
+  class AssembledSparseMatrix;
 enum class ElementDofOrdering;
 
 /** An enum type to specify if only e1 value is requested (SingleValued) or both
@@ -47,6 +48,8 @@ protected:
    Array<int> indices;
    Array<int> gatherMap;
 
+  friend class AssembledSparseMatrix;
+  
 public:
    ElementRestriction(const FiniteElementSpace&, ElementDofOrdering);
    void Mult(const Vector &x, Vector &y) const;
@@ -78,6 +81,7 @@ public:
    /** Fill the J and Data arrays of SparseMatrix corresponding to the sparsity
        pattern given by this ElementRestriction, and the values of ea_data. */
    void FillJAndData(const Vector &ea_data, SparseMatrix &mat) const;
+
 };
 
 /// Operator that converts L2 FiniteElementSpace L-vectors to E-vectors.
@@ -185,6 +189,44 @@ int PermuteFaceL2(const int dim, const int face_id1,
                   const int face_id2, const int orientation,
                   const int size1d, const int index);
 
+  /**
+     Creates a CSR sparse matrix from element matrices assembled usign a mfem::ElementDofOrdering
+   */
+class AssembledSparseMatrix : public mfem::SparseMatrix {
+public:
+
+  /**
+   * @brief AssembledSparseMatrix creates a SparseMatrix based on finite element spaces and ElementDofOrdering
+   *
+   * @param[in] test Test finite element space
+   * @param[in] trial Trial finite element space
+   * @param[in] elem_order ElementDofOrdering chosen for both spaces
+   */
+  AssembledSparseMatrix(const mfem::FiniteElementSpace& test,   // test_elem_dofs * ne * vdim x vdim * test_ndofs
+                        const mfem::FiniteElementSpace& trial,  // trial_elem_dofs * ne * vdim x vdim * trial_ndofs
+                        mfem::ElementDofOrdering        elem_order);
+
+  /// Updates SparseMatrix entries based on new element assembled matrices
+  virtual void FillData(const mfem::Vector& ea_data);
+
+protected:  
+  const mfem::FiniteElementSpace& test_fes;
+  const mfem::FiniteElementSpace& trial_fes;
+  // class local ElementRestriction objects
+  mfem::ElementRestriction        test_restriction;
+  mfem::ElementRestriction        trial_restriction;
+  mfem::ElementDofOrdering        elem_ordering;
+  mfem::Array<int>                ea_map;
+  
+private:
+  int  FillI();
+  void FillJ();
+
+};
+
+
 }
 
+
+
 #endif //MFEM_RESTRICTION