Version 1.2.1

data/escher.mesh

@@ -115,27 +115,27 @@ vertices
 3
 -0.816497 0 -0.57735
 -0.816497 -0.816497 0
--1.70538 0 0
+-1.705385 0 0
 -0.816497 0 0.57735
 -0.816497 0.816497 0
 -0.852693 -0.852693 -1.20589
-0 0 -1.1547
+0 0 -1.154701
 0 -0.816497 -0.57735
 0 -0.816497 0.57735
-0 0 1.1547
+0 0 1.154701
 -0.852693 -0.852693 1.20589
 -0.852693 0.852693 -1.20589
 0 0.816497 -0.57735
 0 0.816497 0.57735
 -0.852693 0.852693 1.20589
 0.816497 -0.816497 0
-0 -1.70538 0
-0 1.70538 0
+0 -1.705385 0
+0 1.705385 0
 0.816497 0.816497 0
 0.852693 -0.852693 -1.20589
 0.816497 0 -0.57735
 0.816497 0 0.57735
 0.852693 -0.852693 1.20589
 0.852693 0.852693 -1.20589
 0.852693 0.852693 1.20589
-1.70538 0 0
+1.705385 0 0

examples/ex1.cpp

@@ -99,19 +99,22 @@ int main (int argc, char *argv[])
    //    using GLVis: "glvis -m refined.mesh -g sol.gf".
    {
       ofstream mesh_ofs("refined.mesh");
+      mesh_ofs.precision(8);
       mesh->Print(mesh_ofs);
       ofstream sol_ofs("sol.gf");
+      sol_ofs.precision(8);
       x.Save(sol_ofs);
    }
 
    // 9. (Optional) Send the solution by socket to a GLVis server.
    char vishost[] = "localhost";
    int  visport   = 19916;
-   osockstream sol_sock (visport, vishost);
+   osockstream sol_sock(visport, vishost);
    if (mesh->Dimension() == 2)
       sol_sock << "fem2d_gf_data\n";
    else
       sol_sock << "fem3d_gf_data\n";
+   sol_sock.precision(8);
    mesh->Print(sol_sock);
    x.Save(sol_sock);
    sol_sock.send();

examples/ex1p.cpp

@@ -145,13 +145,15 @@ int main (int argc, char *argv[])
       ofstream mesh_ofs;
       if (myid == 0)
          mesh_ofs.open("refined.mesh");
+      mesh_ofs.precision(8);
       pmesh->PrintAsOne(mesh_ofs);
       if (myid == 0)
          mesh_ofs.close();
 
       ofstream sol_ofs;
       if (myid == 0)
          sol_ofs.open("sol.gf");
+      sol_ofs.precision(8);
       x.SaveAsOne(sol_ofs);
       if (myid == 0)
          sol_ofs.close();
@@ -168,6 +170,7 @@ int main (int argc, char *argv[])
          *sol_sock << "fem2d_gf_data\n";
       else
          *sol_sock << "fem3d_gf_data\n";
+      sol_sock->precision(8);
    }
    pmesh->PrintAsOne(*sol_sock);
    x.SaveAsOne(*sol_sock);

examples/ex2.cpp

@@ -56,6 +56,14 @@ int main (int argc, char *argv[])
    mesh = new Mesh(imesh, 1, 1);
    imesh.close();
 
+   if (mesh->attributes.Max() < 2 || mesh->bdr_attributes.Max() < 2)
+   {
+      cerr << "\nInput mesh should have at least two materials and "
+           << "two boundary attributes! (See schematic in ex2.cpp)\n"
+           << endl;
+      return 3;
+   }
+
    int dim = mesh->Dimension();
 
    // 2. Refine the mesh to increase the resolution. In this example we do
@@ -170,8 +178,10 @@ int main (int argc, char *argv[])
       *nodes += x;
       x *= -1;
       ofstream mesh_ofs("displaced.mesh");
+      mesh_ofs.precision(8);
       mesh->Print(mesh_ofs);
       ofstream sol_ofs("sol.gf");
+      sol_ofs.precision(8);
       x.Save(sol_ofs);
    }
 
@@ -181,11 +191,12 @@ int main (int argc, char *argv[])
    //     the displacements.
    char vishost[] = "localhost";
    int  visport   = 19916;
-   osockstream sol_sock (visport, vishost);
+   osockstream sol_sock(visport, vishost);
    if (dim == 2)
       sol_sock << "vfem2d_gf_data\n";
    else
       sol_sock << "vfem3d_gf_data\n";
+   sol_sock.precision(8);
    mesh->Print(sol_sock);
    x.Save(sol_sock);
    sol_sock.send();

examples/ex2p.cpp

@@ -68,6 +68,15 @@ int main (int argc, char *argv[])
    mesh = new Mesh(imesh, 1, 1);
    imesh.close();
 
+   if (mesh->attributes.Max() < 2 || mesh->bdr_attributes.Max() < 2)
+   {
+      if (myid == 0)
+         cerr << "\nInput mesh should have at least two materials and "
+              << "two boundary attributes! (See schematic in ex2.cpp)\n"
+              << endl;
+      MPI_Finalize();
+      return 3;
+   }
    int dim = mesh->Dimension();
 
    // 3. Refine the serial mesh on all processors to increase the resolution. In
@@ -228,13 +237,15 @@ int main (int argc, char *argv[])
       ofstream mesh_ofs;
       if (myid == 0)
          mesh_ofs.open("displaced.mesh");
+      mesh_ofs.precision(8);
       pmesh->PrintAsOne(mesh_ofs);
       if (myid == 0)
          mesh_ofs.close();
 
       ofstream sol_ofs;
       if (myid == 0)
          sol_ofs.open("sol.gf");
+      sol_ofs.precision(8);
       x.SaveAsOne(sol_ofs);
       if (myid == 0)
          sol_ofs.close();
@@ -254,6 +265,7 @@ int main (int argc, char *argv[])
          *sol_sock << "vfem2d_gf_data\n";
       else
          *sol_sock << "vfem3d_gf_data\n";
+      sol_sock->precision(8);
    }
    pmesh->PrintAsOne(*sol_sock);
    x.SaveAsOne(*sol_sock);

examples/ex3.cpp

@@ -130,16 +130,19 @@ int main (int argc, char *argv[])
    //     later using GLVis: "glvis -m refined.mesh -g sol.gf".
    {
       ofstream mesh_ofs("refined.mesh");
+      mesh_ofs.precision(8);
       mesh->Print(mesh_ofs);
       ofstream sol_ofs("sol.gf");
+      sol_ofs.precision(8);
       dx.Save(sol_ofs);
    }
 
    // 11. (Optional) Send the solution by socket to a GLVis server.
    char vishost[] = "localhost";
    int  visport   = 19916;
-   osockstream sol_sock (visport, vishost);
+   osockstream sol_sock(visport, vishost);
    sol_sock << "vfem3d_gf_data\n";
+   sol_sock.precision(8);
    mesh->Print(sol_sock);
    dx.Save(sol_sock);
    sol_sock.send();

examples/ex3p.cpp

@@ -185,13 +185,15 @@ int main (int argc, char *argv[])
       ofstream mesh_ofs;
       if (myid == 0)
          mesh_ofs.open("refined.mesh");
+      mesh_ofs.precision(8);
       pmesh->PrintAsOne(mesh_ofs);
       if (myid == 0)
          mesh_ofs.close();
 
       ofstream sol_ofs;
       if (myid == 0)
          sol_ofs.open("sol.gf");
+      sol_ofs.precision(8);
       dx.SaveAsOne(sol_ofs);
       if (myid == 0)
          sol_ofs.close();
@@ -205,6 +207,7 @@ int main (int argc, char *argv[])
    {
       sol_sock = new osockstream(visport, vishost);
       *sol_sock << "vfem3d_gf_data\n";
+      sol_sock->precision(8);
    }
    pmesh->PrintAsOne(*sol_sock);
    dx.SaveAsOne(*sol_sock);

linalg/hypre.cpp

@@ -455,6 +455,26 @@ void HypreParMatrix::Mult(const Vector &x, Vector &y) const
    hypre_ParCSRMatrixMatvec(1.0, A, *X, 0.0, *Y);
 }
 
+void HypreParMatrix::MultTranspose(const Vector &x, Vector &y) const
+{
+   if (X == NULL)
+   {
+      X = new HypreParVector(GetGlobalNumCols(),
+                             x.GetData(),
+                             GetColStarts());
+      Y = new HypreParVector(GetGlobalNumRows(),
+                             y.GetData(),
+                             GetRowStarts());
+   }
+   else
+   {
+      X -> SetData(x.GetData());
+      Y -> SetData(y.GetData());
+   }
+
+   hypre_ParCSRMatrixMatvecT(1.0, A, *X, 0.0, *Y);
+}
+
 int HypreParMatrix::Mult(HYPRE_ParVector x, HYPRE_ParVector y,
                          double a, double b)
 {

linalg/hypre.hpp

@@ -164,6 +164,7 @@ class HypreParMatrix : public Operator
                      double alpha = 1.0, double beta = 0.0);
 
    virtual void Mult(const Vector &x, Vector &y) const;
+   virtual void MultTranspose(const Vector &x, Vector &y) const;
 
    /// Prints the locally owned rows in parallel
    void Print(const char *fname, int offi = 0, int offj = 0);
@@ -247,6 +248,7 @@ class HyprePCG : public HypreSolver
 
    /// Solve Ax=b with hypre's PCG
    virtual void Mult(const HypreParVector &b, HypreParVector &x) const;
+   using HypreSolver::Mult;
 
    virtual ~HyprePCG();
 };
@@ -285,6 +287,7 @@ class HypreGMRES : public HypreSolver
 
    /// Solve Ax=b with hypre's GMRES
    virtual void Mult (const HypreParVector &b, HypreParVector &x) const;
+   using HypreSolver::Mult;
 
    virtual ~HypreGMRES();
 };

mesh/mesh.cpp

@@ -1756,6 +1756,14 @@ void Mesh::Load(istream &input, int generate_edges, int refine)
       }
    }
 
+   if (NumOfBdrElements == 0 && Dim > 2)
+   {
+      // in 3D, generate boundary elements before we 'MarkForRefinement'
+      GetElementToFaceTable();
+      GenerateFaces();
+      GenerateBoundaryElements();
+   }
+
    if (!curved)
    {
       // check and fix element orientation
@@ -1770,8 +1778,6 @@ void Mesh::Load(istream &input, int generate_edges, int refine)
    {
       GetElementToFaceTable();
       GenerateFaces();
-      if (NumOfBdrElements == 0)
-         GenerateBoundaryElements();
       // check and fix boundary element orientation
       if ( !(curved && (meshgen & 1)) )
          CheckBdrElementOrientation();
@@ -3031,6 +3037,35 @@ void Mesh::GenerateFaces()
    }
 }
 
+STable3D *Mesh::GetFacesTable()
+{
+   STable3D *faces_tbl = new STable3D(NumOfVertices);
+   for (int i = 0; i < NumOfElements; i++)
+   {
+      const int *v = elements[i]->GetVertices();
+      switch (GetElementType(i))
+      {
+      case Element::TETRAHEDRON:
+         faces_tbl->Push(v[1], v[2], v[3]);
+         faces_tbl->Push(v[0], v[3], v[2]);
+         faces_tbl->Push(v[0], v[1], v[3]);
+         faces_tbl->Push(v[0], v[2], v[1]);
+         break;
+      case Element::HEXAHEDRON:
+         // find the face by the vertices with the smallest 3 numbers
+         // z = 0, y = 0, x = 1, y = 1, x = 0, z = 1
+         faces_tbl->Push4(v[3], v[2], v[1], v[0]);
+         faces_tbl->Push4(v[0], v[1], v[5], v[4]);
+         faces_tbl->Push4(v[1], v[2], v[6], v[5]);
+         faces_tbl->Push4(v[2], v[3], v[7], v[6]);
+         faces_tbl->Push4(v[3], v[0], v[4], v[7]);
+         faces_tbl->Push4(v[4], v[5], v[6], v[7]);
+         break;
+      }
+   }
+   return faces_tbl;
+}
+
 STable3D *Mesh::GetElementToFaceTable(int ret_ftbl)
 {
    int i, *v;

mesh/mesh.hpp

@@ -92,6 +92,7 @@ class Mesh
    void MarkTriMeshForRefinement();
    void MarkTetMeshForRefinement();
 
+   STable3D *GetFacesTable();
    STable3D *GetElementToFaceTable(int ret_ftbl = 0);
 
    /** Red refinement. Element with index i is refined. The default

mesh/pmesh.cpp

@@ -1396,6 +1396,8 @@ void ParMesh::RefineGroups(const DSTable &v_to_v, int *middle)
 
 void ParMesh::QuadUniformRefinement()
 {
+   SetState(Mesh::NORMAL);
+
    int oedge = NumOfVertices, wtls = WantTwoLevelState;
 
    if (Nodes)  // curved mesh
@@ -1484,13 +1486,15 @@ void ParMesh::QuadUniformRefinement()
 
 void ParMesh::HexUniformRefinement()
 {
+   SetState(Mesh::NORMAL);
+
    int wtls = WantTwoLevelState;
    int oedge = NumOfVertices;
    int oface = oedge + NumOfEdges;
 
    DSTable v_to_v(NumOfVertices);
    GetVertexToVertexTable(v_to_v);
-   STable3D *faces_tbl = GetElementToFaceTable(1);
+   STable3D *faces_tbl = GetFacesTable();
 
    if (Nodes)  // curved mesh
       UseTwoLevelState(1);
@@ -1611,7 +1615,7 @@ void ParMesh::HexUniformRefinement()
       }
 
       delete faces_tbl;
-      faces_tbl = GetElementToFaceTable(1);
+      faces_tbl = GetFacesTable();
       for (i = 0; i < shared_faces.Size(); i++)
       {
          shared_faces[i]->GetVertices(v);