More debug outputs.

data/inline-hex.mesh

@@ -1,9 +1,9 @@
 MFEM INLINE mesh v1.0
 
 type = hex
-nx = 4
-ny = 4
-nz = 4
-sx = 1.0
+nx = 2
+ny = 1
+nz = 1
+sx = 2.0
 sy = 1.0
 sz = 1.0

data/inline-quad.mesh

@@ -1,7 +1,7 @@
 MFEM INLINE mesh v1.0
 
 type = quad
-nx = 4
-ny = 4
-sx = 1.0
+nx = 2
+ny = 1
+sx = 2.0
 sy = 1.0

examples/ex1p.cpp

@@ -64,6 +64,7 @@ int main(int argc, char *argv[])
    const char *mesh_file = "../data/star.mesh";
    int order = 1;
    bool static_cond = false;
+   bool hybridization = false;
    bool pa = false;
    const char *device_config = "cpu";
    bool visualization = true;
@@ -76,6 +77,8 @@ int main(int argc, char *argv[])
                   " isoparametric space.");
    args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
                   "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&hybridization, "-hb", "--hybridization", "-no-hb",
+                  "--no-hybridization", "Enable hybridization.");
    args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
                   "--no-partial-assembly", "Enable Partial Assembly.");
    args.AddOption(&device_config, "-d", "--device",
@@ -114,21 +117,25 @@ int main(int argc, char *argv[])
    //    'ref_levels' to be the largest number that gives a final mesh with no
    //    more than 10,000 elements.
    {
-      int ref_levels =
-         (int)floor(log(10000./mesh->GetNE())/log(2.)/dim);
+      int ref_levels = 0;
+         //(int)floor(log(10000./mesh->GetNE())/log(2.)/dim);
       for (int l = 0; l < ref_levels; l++)
       {
          mesh->UniformRefinement();
       }
    }
 
+   Array<int> refs;
+   refs.Append(0);
+   mesh->GeneralRefinement(refs);
+
    // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
    //    this mesh further in parallel to increase the resolution. Once the
    //    parallel mesh is defined, the serial mesh can be deleted.
    ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
    delete mesh;
    {
-      int par_ref_levels = 2;
+      int par_ref_levels = 0;
       for (int l = 0; l < par_ref_levels; l++)
       {
          pmesh->UniformRefinement();
@@ -199,13 +206,31 @@ int main(int argc, char *argv[])
    //     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(); }
+   FiniteElementCollection *hfec = NULL;
+   ParFiniteElementSpace *hfes = NULL;
+   if (static_cond)
+   {
+      a->EnableStaticCondensation();
+   }
+   else if (hybridization)
+   {
+      hfec = new DG_Interface_FECollection(order, dim);
+      hfes = new ParFiniteElementSpace(pmesh, hfec);
+      a->EnableHybridization(hfes, new TraceJumpIntegrator(),
+                             ess_tdof_list);
+   }
    a->Assemble();
 
    OperatorPtr A;
    Vector B, X;
    a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
 
+   HYPRE_Int glob_size = A->Height();
+   if (myid == 0)
+   {
+      cout << "Size of linear system: " << glob_size << endl;
+   }
+
    // 13. Solve the linear system A X = B.
    //     * With full assembly, use the BoomerAMG preconditioner from hypre.
    //     * With partial assembly, use no preconditioner, for now.
@@ -252,6 +277,8 @@ int main(int argc, char *argv[])
    }
 
    // 17. Free the used memory.
+   delete hfes;
+   delete hfec;
    delete a;
    delete b;
    delete fespace;

examples/ex4p.cpp

@@ -106,14 +106,18 @@ int main(int argc, char *argv[])
    //    'ref_levels' to be the largest number that gives a final mesh with no
    //    more than 1,000 elements.
    {
-      int ref_levels =
-         (int)floor(log(1000./mesh->GetNE())/log(2.)/dim);
+      int ref_levels = 0;
+         //(int)floor(log(1000./mesh->GetNE())/log(2.)/dim);
       for (int l = 0; l < ref_levels; l++)
       {
          mesh->UniformRefinement();
       }
    }
 
+   Array<int> refs;
+   refs.Append(0);
+   mesh->GeneralRefinement(refs);
+
    // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
    //    this mesh further in parallel to increase the resolution. Once the
    //    parallel mesh is defined, the serial mesh can be deleted. Tetrahedral
@@ -122,7 +126,7 @@ int main(int argc, char *argv[])
    ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
    delete mesh;
    {
-      int par_ref_levels = 2;
+      int par_ref_levels = 0;
       for (int l = 0; l < par_ref_levels; l++)
       {
          pmesh->UniformRefinement();

fem/eltrans.hpp

@@ -102,6 +102,8 @@ public:
    virtual int TransformBack(const Vector &pt, IntegrationPoint &ip) = 0;
 
    virtual ~ElementTransformation() { }
+
+   virtual void Print(std::ostream &out = mfem::out) const {}
 };
 
 
@@ -331,6 +333,9 @@ public:
    }
 
    virtual ~IsoparametricTransformation() { }
+
+   virtual void Print(std::ostream &out = mfem::out) const
+   { PointMat.PrintMatlab(out); }
 };
 
 class IntegrationPointTransformation

fem/hybridization.cpp

@@ -9,6 +9,8 @@
 // terms of the GNU Lesser General Public License (as published by the Free
 // Software Foundation) version 2.1 dated February 1999.
 
+#include <fstream>
+
 #include "hybridization.hpp"
 #include "gridfunc.hpp"
 
@@ -78,7 +80,7 @@ void Hybridization::ConstructC()
             c_pfes->ExchangeFaceNbrData();
             c_num_face_nbr_dofs = c_pfes->GetFaceNbrVSize();
          }
-#ifdef MFEM_DEBUG_HERE
+#if 1//def MFEM_DEBUG_HERE
          MFEM_WARNING('[' << c_pfes->GetMyRank() <<
                       "] num_shared_slave_faces = " << num_shared_slave_faces
                       << ", glob_num_shared_slave_faces = "
@@ -94,6 +96,10 @@ void Hybridization::ConstructC()
    const int c_vsize = c_fes->GetVSize();
    Ct = new SparseMatrix(num_hat_dofs, c_vsize + c_num_face_nbr_dofs);
 
+   char fname[100];
+   sprintf(fname, "hybrid-%03d.txt", pmesh->GetMyRank());
+   std::ofstream dbg(fname);
+
    if (c_bfi)
    {
       const int skip_zeros = 1;
@@ -106,6 +112,14 @@ void Hybridization::ConstructC()
          FTr = mesh->GetInteriorFaceTransformations(i);
          if (!FTr) { continue; }
 
+         dbg << "Face:\n";
+         FTr->Face->Print(dbg);
+         dbg << "Loc1:\n";
+         FTr->Loc1.Transf.Print(dbg);
+         dbg << "Loc2:\n";
+         FTr->Loc2.Transf.Print(dbg);
+         dbg << "\n";
+
          int o1 = hat_offsets[FTr->Elem1No];
          int s1 = hat_offsets[FTr->Elem1No+1] - o1;
          int o2 = hat_offsets[FTr->Elem2No];
@@ -127,6 +141,10 @@ void Hybridization::ConstructC()
          // zero-out small elements in elmat
          elmat.Threshold(1e-12 * elmat.MaxMaxNorm());
          Ct->AddSubMatrix(vdofs, c_vdofs, elmat, skip_zeros);
+
+         dbg << "Ct submatrix:\n";
+         elmat.PrintMatlab(dbg);
+         dbg << std::endl;
       }
 #ifdef MFEM_USE_MPI
       if (pmesh)
@@ -145,8 +163,9 @@ void Hybridization::ConstructC()
                face_fe = c_fes->GetFaceElement(face_no);
                c_fes->GetFaceVDofs(face_no, c_vdofs);
             }
-            else
+            else // ghost face
             {
+               dbg << "Ghost face!\n";
                const int fill2 = false; // only need side "1" data
                FTr = pmesh->GetSharedFaceTransformations(i, fill2);
                face_fe = c_pfes->GetFaceNbrFaceFE(face_no);
@@ -157,6 +176,13 @@ void Hybridization::ConstructC()
                   c_vdofs[j] += c_vsize;
                }
             }
+
+            dbg << "Face:\n";
+            FTr->Face->Print(dbg);
+            dbg << "Loc1:\n";
+            FTr->Loc1.Transf.Print(dbg);
+            dbg << "\n";
+
             int o1 = hat_offsets[FTr->Elem1No];
             int s1 = hat_offsets[FTr->Elem1No+1] - o1;
             vdofs.SetSize(s1);
@@ -168,7 +194,19 @@ void Hybridization::ConstructC()
             c_bfi->AssembleFaceMatrix(*face_fe, *fe, *fe, *FTr, elmat);
             // zero-out small elements in elmat
             elmat.Threshold(1e-12 * elmat.MaxMaxNorm());
+            /*if (pmesh->GetMyRank() == 1 && ghost_sface)
+            {
+               mfem::out << "TEST\n";
+               for (int i = 0; i < elmat.Height(); i++)
+               {
+                  std::swap(elmat(i,0), elmat(i,1));
+               }
+            }*/
             Ct->AddSubMatrix(vdofs, c_vdofs, elmat, skip_zeros);
+
+            dbg << "Shared face " << i << " Ct submatrix:\n";
+            elmat.PrintMatlab(dbg);
+            dbg << std::endl;
          }
          if (glob_num_shared_slave_faces)
          {

fem/pfespace.cpp

@@ -1132,13 +1132,16 @@ void ParFiniteElementSpace::GetFaceNbrFaceVDofs(int i, Array<int> &vdofs) const
    // ParMesh::GetSharedFace() such that i >= Mesh::GetNumFaces(), i.e. 'i' is
    // the index of a ghost.
    MFEM_ASSERT(Nonconforming() && i >= pmesh->GetNumFaces(), "");
+
    int el1, el2, inf1, inf2;
    pmesh->GetFaceElements(i, &el1, &el2);
    el2 = -1 - el2;
    pmesh->GetFaceInfos(i, &inf1, &inf2);
-   MFEM_ASSERT(0 <= el2 && el2 < face_nbr_element_dof.Size(), "");
+
+   MFEM_ASSERT(el2 >= 0 && el2 < face_nbr_element_dof.Size(), "");
    const int nd = face_nbr_element_dof.RowSize(el2);
    const int *vol_vdofs = face_nbr_element_dof.GetRow(el2);
+
    const Element *face_nbr_el = pmesh->face_nbr_elements[el2];
    Geometry::Type geom = face_nbr_el->GetGeometryType();
    const int face_dim = Geometry::Dimension[geom]-1;
@@ -1174,6 +1177,7 @@ const FiniteElement *ParFiniteElementSpace::GetFaceNbrFaceFE(int i) const
    MFEM_ASSERT(Nonconforming() && !NURBSext, "");
    Geometry::Type geom = (pmesh->Dimension() == 2) ?
                          Geometry::SEGMENT : Geometry::SQUARE;
+   // TODO: fix the above in nc-prism-dev    ^
    return fec->FiniteElementForGeometry(geom);
 }
 

mesh/pmesh.cpp

@@ -2341,10 +2341,10 @@ GetSharedFaceTransformations(int sf, bool fill2)
 
       if (face_type == Element::SEGMENT && fill2)
       {
-         // fix slave orientation in 2D: flip Loc2 to match Loc1 and Face
+/*         // fix slave orientation in 2D: flip Loc2 to match Loc1 and Face
          DenseMatrix &pm = FaceElemTr.Loc2.Transf.GetPointMat();
          std::swap(pm(0,0), pm(0,1));
-         std::swap(pm(1,0), pm(1,1));
+         std::swap(pm(1,0), pm(1,1));*/
       }
    }
 

mesh/pncmesh.cpp

@@ -1024,7 +1024,7 @@ void ParNCMesh::GetFaceNeighbors(ParMesh &pmesh)
    // same on different processors, this is important for ExchangeFaceNbrData)
    fnbr.Sort();
    fnbr.Unique();
-   fnbr.Sort(compare_ranks_indices);
+   fnbr.Sort(compare_ranks_indices); // TODO: lambda?
 
    // put the ranks into 'face_nbr_group'
    for (int i = 0; i < fnbr.Size(); i++)
@@ -1192,16 +1192,20 @@ void ParNCMesh::GetFaceNeighbors(ParMesh &pmesh)
             fi.Elem2No = -1 - fnbr_index[fi.Elem2No - NElements];
 
             const DenseMatrix* pm = &sf.point_matrix;
-            if (!sloc && Dim == 3)
+            if (!sloc)
             {
-               // ghost slave in 3D needs flipping orientation
-               DenseMatrix* pm2 = new DenseMatrix(*pm);
-               std::swap((*pm2)(0,1), (*pm2)(0,3));
-               std::swap((*pm2)(1,1), (*pm2)(1,3));
-               aux_pm_store.Append(pm2);
+               if (Dim == 3)
+               {
+                  // ghost slave in 3D needs flipping orientation
+                  DenseMatrix* pm2 = new DenseMatrix(*pm);
+                  std::swap((*pm2)(0,1), (*pm2)(0,3));
+                  std::swap((*pm2)(1,1), (*pm2)(1,3));
+                  aux_pm_store.Append(pm2);
 
-               fi.Elem2Inf ^= 1;
-               pm = pm2;
+                  pm = pm2;
+               }
+
+               //fi.Elem2Inf ^= 1;
 
                // The problem is that sf.point_matrix is designed for P matrix
                // construction and always has orientation relative to the slave