Merge branch 'master' into mixednonlin-dev

.gitignore

@@ -114,10 +114,12 @@ examples/petsc/elastic_energy.*
 
 miniapps/electromagnetics/volta
 miniapps/electromagnetics/tesla
+miniapps/electromagnetics/maxwell
 miniapps/electromagnetics/joule
 
 miniapps/electromagnetics/Volta-AMR*
 miniapps/electromagnetics/Tesla-AMR*
+miniapps/electromagnetics/Maxwell-Parallel*
 miniapps/electromagnetics/Joule_*
 
 miniapps/meshing/mobius-strip
@@ -143,3 +145,12 @@ miniapps/performance/mesh.*
 miniapps/performance/sol.*
 
 miniapps/tools/display-basis
+
+miniapps/nurbs/ex1
+miniapps/nurbs/ex1p
+miniapps/nurbs/ex11p
+miniapps/nurbs/refined.mesh
+miniapps/nurbs/mesh.*
+miniapps/nurbs/sol.*
+miniapps/nurbs/mode_*
+miniapps/nurbs/Example1*

CHANGELOG

@@ -13,6 +13,11 @@ Version 3.3.3 (development)
 
 Discretization improvements
 ---------------------------
+- Added variable order NURBS: for each space each knot vector in the mesh can
+  have a different order. The order information is now part of the finite
+  element space header in the NURBS mesh output, so NURBS meshes in the old
+  format need to be updated.
+
 - Added classes for general block nonlinear finite element operators (deriving
   from BlockNonlinearForm/ParBlockNonlinearForm). Such operators have
   assemble-based action and also support assembly of the gradient operator to
@@ -26,10 +31,22 @@ Discretization improvements
 
 New and updated examples and miniapps
 -------------------------------------
+- Added a set of miniapps specifically focused on Isogeometric Analysis (IGA) on
+  NURBS meshes in the miniapps/nurbs directory. Currently the directory contains
+  variable order NURBS versions of examples 1, 1p and 11p.
+
+- Added a new electromagnetics miniapp, Maxwell, for simulating time-domain
+  electromagnetics phenomena as a coupled first order system of equations.
+
 - Added a new serial and parallel example (ex19) that solves the quasi-static
   incompressible hyperelastic equations. The example demonstrates the use of
   block nonlinear forms as well as custom block preconditioners.
 
+New and improved solvers and preconditioners
+--------------------------------------------
+- Added symplectic integrators of orders 1-4 for systems of first order ODEs
+  derived from a Hamiltonian, see class SIASolver in linalg/ode.hpp.
+
 Miscellaneous
 -------------
 - Added a new developer tool, config/sample-runs.sh, that extracts the sample

CMakeLists.txt

@@ -405,12 +405,12 @@ endif()
 # Add 'check' target - quick test
 if (NOT MFEM_USE_MPI)
   add_custom_target(check
-    ${CMAKE_CTEST_COMMAND} -R ex1_ser -E performance -C ${CMAKE_CFG_INTDIR}
+    ${CMAKE_CTEST_COMMAND} -R '^ex1_ser' -C ${CMAKE_CFG_INTDIR}
     USES_TERMINAL)
   add_dependencies(check ex1)
 else()
   add_custom_target(check
-    ${CMAKE_CTEST_COMMAND} -R ex1p -E performance -C ${CMAKE_CFG_INTDIR}
+    ${CMAKE_CTEST_COMMAND} -R '^ex1p' -C ${CMAKE_CFG_INTDIR}
     USES_TERMINAL)
   add_dependencies(check ex1p)
 endif()

doc/CodeDocumentation.conf.in

@@ -771,6 +771,7 @@ INPUT                  = @MFEM_SOURCE_DIR@/doc/CodeDocumentation.dox \
                          @MFEM_SOURCE_DIR@/miniapps/common \
                          @MFEM_SOURCE_DIR@/miniapps/meshing \
                          @MFEM_SOURCE_DIR@/miniapps/tools \
+                         @MFEM_SOURCE_DIR@/miniapps/nurbs \
                          @MFEM_SOURCE_DIR@/miniapps/electromagnetics \
                          @MFEM_SOURCE_DIR@/miniapps/performance
 

doc/CodeDocumentation.dox

@@ -99,6 +99,7 @@ namespace mfem {
  * <H3>Miniapps</H3>
  * - <a class="el" href="volta_8cpp_source.html">Volta</a>: simple electrostatics simulation code
  * - <a class="el" href="tesla_8cpp_source.html">Tesla</a>: simple magnetostatics simulation code
+ * - <a class="el" href="maxwell_8cpp_source.html">Maxwell</a>: simple transient full-wave electromagnetics simulation code
  * - <a class="el" href="joule_8cpp_source.html">Joule</a>: transient magnetics and Joule heating miniapp
  * - <a class="el" href="mobius-strip_8cpp_source.html">Mobius Strip</a>: generate various Mobius strip-like meshes
  * - <a class="el" href="klein-bottle_8cpp_source.html">Klein Bottle</a>: generate three types of Klein bottle surfaces

examples/ex12p.cpp

@@ -109,9 +109,9 @@ int main(int argc, char *argv[])
 
    // 4. Select the order of the finite element discretization space. For NURBS
    //    meshes, we increase the order by degree elevation.
-   if (mesh->NURBSext && order > mesh->NURBSext->GetOrder())
+   if (mesh->NURBSext)
    {
-      mesh->DegreeElevate(order - mesh->NURBSext->GetOrder());
+      mesh->DegreeElevate(order, order);
    }
 
    // 5. Refine the serial mesh on all processors to increase the resolution. In

examples/ex2.cpp

@@ -84,9 +84,9 @@ int main(int argc, char *argv[])
 
    // 3. Select the order of the finite element discretization space. For NURBS
    //    meshes, we increase the order by degree elevation.
-   if (mesh->NURBSext && order > mesh->NURBSext->GetOrder())
+   if (mesh->NURBSext)
    {
-      mesh->DegreeElevate(order - mesh->NURBSext->GetOrder());
+      mesh->DegreeElevate(order, order);
    }
 
    // 4. Refine the mesh to increase the resolution. In this example we do

examples/ex2p.cpp

@@ -107,9 +107,9 @@ int main(int argc, char *argv[])
 
    // 4. Select the order of the finite element discretization space. For NURBS
    //    meshes, we increase the order by degree elevation.
-   if (mesh->NURBSext && order > mesh->NURBSext->GetOrder())
+   if (mesh->NURBSext)
    {
-      mesh->DegreeElevate(order - mesh->NURBSext->GetOrder());
+      mesh->DegreeElevate(order, order);
    }
 
    // 5. Refine the serial mesh on all processors to increase the resolution. In

examples/petsc/ex2p.cpp

@@ -121,9 +121,9 @@ int main(int argc, char *argv[])
 
    // 4. Select the order of the finite element discretization space. For NURBS
    //    meshes, we increase the order by degree elevation.
-   if (mesh->NURBSext && order > mesh->NURBSext->GetOrder())
+   if (mesh->NURBSext)
    {
-      mesh->DegreeElevate(order - mesh->NURBSext->GetOrder());
+      mesh->DegreeElevate(order, order);
    }
 
    // 5. Refine the serial mesh on all processors to increase the resolution. In

fem/bilinearform.cpp

@@ -205,6 +205,14 @@ void BilinearForm::AddDomainIntegrator (BilinearFormIntegrator * bfi)
 void BilinearForm::AddBoundaryIntegrator (BilinearFormIntegrator * bfi)
 {
    bbfi.Append (bfi);
+   bbfi_marker.Append(NULL); // NULL marker means apply everywhere
+}
+
+void BilinearForm::AddBoundaryIntegrator (BilinearFormIntegrator * bfi,
+                                          Array<int> &bdr_marker)
+{
+   bbfi.Append (bfi);
+   bbfi_marker.Append(&bdr_marker);
 }
 
 void BilinearForm::AddInteriorFaceIntegrator (BilinearFormIntegrator * bfi)
@@ -357,14 +365,41 @@ void BilinearForm::Assemble (int skip_zeros)
 
    if (bbfi.Size())
    {
+      // Which boundary attributes need to be processed?
+      Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
+                                 mesh->bdr_attributes.Max() : 0);
+      bdr_attr_marker = 0;
+      for (int k = 0; k < bbfi.Size(); k++)
+      {
+         if (bbfi_marker[k] == NULL)
+         {
+            bdr_attr_marker = 1;
+            break;
+         }
+         Array<int> &bdr_marker = *bbfi_marker[k];
+         MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
+                     "invalid boundary marker for boundary integrator #"
+                     << k << ", counting from zero");
+         for (int i = 0; i < bdr_attr_marker.Size(); i++)
+         {
+            bdr_attr_marker[i] |= bdr_marker[i];
+         }
+      }
+
       for (i = 0; i < fes -> GetNBE(); i++)
       {
+         const int bdr_attr = mesh->GetBdrAttribute(i);
+         if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
+
          const FiniteElement &be = *fes->GetBE(i);
          fes -> GetBdrElementVDofs (i, vdofs);
          eltrans = fes -> GetBdrElementTransformation (i);
          bbfi[0]->AssembleElementMatrix(be, *eltrans, elmat);
          for (int k = 1; k < bbfi.Size(); k++)
          {
+            if (bbfi_marker[k] &&
+                (*bbfi_marker[k])[bdr_attr-1] == 0) { continue; }
+
             bbfi[k]->AssembleElementMatrix(be, *eltrans, elemmat);
             elmat += elemmat;
          }

fem/bilinearform.hpp

@@ -49,6 +49,7 @@ protected:
 
    /// Set of Boundary Integrators to be applied.
    Array<BilinearFormIntegrator*> bbfi;
+   Array<Array<int>*>             bbfi_marker;
 
    /// Set of interior face Integrators to be applied.
    Array<BilinearFormIntegrator*> fbfi;
@@ -212,6 +213,11 @@ public:
    /// Adds new Boundary Integrator.
    void AddBoundaryIntegrator(BilinearFormIntegrator *bfi);
 
+   /** @brief Adds new Boundary Integrator, restricted to specific boundary
+       attributes. */
+   void AddBoundaryIntegrator(BilinearFormIntegrator * bfi,
+                              Array<int> &bdr_marker);
+
    /// Adds new interior Face Integrator.
    void AddInteriorFaceIntegrator(BilinearFormIntegrator *bfi);
 

fem/bilininteg.cpp

@@ -25,7 +25,7 @@ void BilinearFormIntegrator::AssembleElementMatrix (
    DenseMatrix &elmat )
 {
    mfem_error ("BilinearFormIntegrator::AssembleElementMatrix (...)\n"
-               "   is not implemented fot this class.");
+               "   is not implemented for this class.");
 }
 
 void BilinearFormIntegrator::AssembleElementMatrix2 (
@@ -33,7 +33,7 @@ void BilinearFormIntegrator::AssembleElementMatrix2 (
    ElementTransformation &Trans, DenseMatrix &elmat )
 {
    mfem_error ("BilinearFormIntegrator::AssembleElementMatrix2 (...)\n"
-               "   is not implemented fot this class.");
+               "   is not implemented for this class.");
 }
 
 void BilinearFormIntegrator::AssembleFaceMatrix (
@@ -41,7 +41,7 @@ void BilinearFormIntegrator::AssembleFaceMatrix (
    FaceElementTransformations &Trans, DenseMatrix &elmat)
 {
    mfem_error ("BilinearFormIntegrator::AssembleFaceMatrix (...)\n"
-               "   is not implemented fot this class.");
+               "   is not implemented for this class.");
 }
 
 void BilinearFormIntegrator::AssembleFaceMatrix(
@@ -58,7 +58,7 @@ void BilinearFormIntegrator::AssembleElementVector(
    Vector &elvect)
 {
    mfem_error("BilinearFormIntegrator::AssembleElementVector\n"
-              "   is not implemented fot this class.");
+              "   is not implemented for this class.");
 }
 
 
@@ -1212,7 +1212,7 @@ void VectorFEWeakDivergenceIntegrator::AssembleElementMatrix2(
       //   n = 2*(d-1)*k+(l-1)+(m-1)
       //
       // In the next formula we use the expressions for n with k=1, which means
-      // that the term Q/det(J) is disregard:
+      // that the term Q/det(J) is disregarded:
       int ir_order = (trial_fe.Space() == FunctionSpace::Pk) ?
                      (trial_fe.GetOrder() + test_fe.GetOrder() - 1) :
                      (trial_fe.GetOrder() + test_fe.GetOrder() + 2*(dim-2));

fem/coefficient.hpp

@@ -367,7 +367,7 @@ public:
 /// Vector coefficient defined by a vector GridFunction
 class VectorGridFunctionCoefficient : public VectorCoefficient
 {
-private:
+protected:
    GridFunction *GridFunc;
 
 public:

fem/fe.cpp

@@ -31,6 +31,7 @@ FiniteElement::FiniteElement(int D, int G, int Do, int O, int F)
    DerivType = NONE;
    DerivRangeType = SCALAR;
    DerivMapType = VALUE;
+   for (int i = 0; i < Geometry::MaxDim; i++) { Orders[i] = -1; }
 #ifndef MFEM_THREAD_SAFE
    vshape.SetSize(Dof, Dim);
 #endif
@@ -10883,6 +10884,14 @@ void ND_SegmentElement::CalcVShape(const IntegrationPoint &ip,
    obasis1d.Eval(ip.x, vshape);
 }
 
+void NURBS1DFiniteElement::SetOrder() const
+{
+   Order = kv[0]->GetOrder();
+   Dof = Order + 1;
+
+   weights.SetSize(Dof);
+   shape_x.SetSize(Dof);
+}
 
 void NURBS1DFiniteElement::CalcShape(const IntegrationPoint &ip,
                                      Vector &shape) const
@@ -10917,6 +10926,21 @@ void NURBS1DFiniteElement::CalcDShape(const IntegrationPoint &ip,
    add(sum, grad, -dsum*sum*sum, shape_x, grad);
 }
 
+void NURBS2DFiniteElement::SetOrder() const
+{
+   Orders[0] = kv[0]->GetOrder();
+   Orders[1] = kv[1]->GetOrder();
+   shape_x.SetSize(Orders[0]+1);
+   shape_y.SetSize(Orders[1]+1);
+   dshape_x.SetSize(Orders[0]+1);
+   dshape_y.SetSize(Orders[1]+1);
+
+   Order = max(Orders[0], Orders[1]);
+   Dof = (Orders[0] + 1)*(Orders[1] + 1);
+   u.SetSize(Dof);
+   weights.SetSize(Dof);
+}
+
 void NURBS2DFiniteElement::CalcShape(const IntegrationPoint &ip,
                                      Vector &shape) const
 {
@@ -10924,10 +10948,10 @@ void NURBS2DFiniteElement::CalcShape(const IntegrationPoint &ip,
    kv[1]->CalcShape(shape_y, ijk[1], ip.y);
 
    double sum = 0.0;
-   for (int o = 0, j = 0; j <= Order; j++)
+   for (int o = 0, j = 0; j <= Orders[1]; j++)
    {
       const double sy = shape_y(j);
-      for (int i = 0; i <= Order; i++, o++)
+      for (int i = 0; i <= Orders[0]; i++, o++)
       {
          sum += ( shape(o) = shape_x(i)*sy*weights(o) );
       }
@@ -10948,10 +10972,10 @@ void NURBS2DFiniteElement::CalcDShape(const IntegrationPoint &ip,
    kv[1]->CalcDShape(dshape_y, ijk[1], ip.y);
 
    sum = dsum[0] = dsum[1] = 0.0;
-   for (int o = 0, j = 0; j <= Order; j++)
+   for (int o = 0, j = 0; j <= Orders[1]; j++)
    {
       const double sy = shape_y(j), dsy = dshape_y(j);
-      for (int i = 0; i <= Order; i++, o++)
+      for (int i = 0; i <= Orders[0]; i++, o++)
       {
          sum += ( u(o) = shape_x(i)*sy*weights(o) );
 
@@ -10971,6 +10995,26 @@ void NURBS2DFiniteElement::CalcDShape(const IntegrationPoint &ip,
    }
 }
 
+//---------------------------------------------------------------------
+void NURBS3DFiniteElement::SetOrder() const
+{
+   Orders[0] = kv[0]->GetOrder();
+   Orders[1] = kv[1]->GetOrder();
+   Orders[2] = kv[2]->GetOrder();
+   shape_x.SetSize(Orders[0]+1);
+   shape_y.SetSize(Orders[1]+1);
+   shape_z.SetSize(Orders[2]+1);
+
+   dshape_x.SetSize(Orders[0]+1);
+   dshape_y.SetSize(Orders[1]+1);
+   dshape_z.SetSize(Orders[2]+1);
+
+   Order = max(max(Orders[0], Orders[1]), Orders[2]);
+   Dof = (Orders[0] + 1)*(Orders[1] + 1)*(Orders[2] + 1);
+   u.SetSize(Dof);
+   weights.SetSize(Dof);
+}
+
 void NURBS3DFiniteElement::CalcShape(const IntegrationPoint &ip,
                                      Vector &shape) const
 {
@@ -10979,13 +11023,13 @@ void NURBS3DFiniteElement::CalcShape(const IntegrationPoint &ip,
    kv[2]->CalcShape(shape_z, ijk[2], ip.z);
 
    double sum = 0.0;
-   for (int o = 0, k = 0; k <= Order; k++)
+   for (int o = 0, k = 0; k <= Orders[2]; k++)
    {
       const double sz = shape_z(k);
-      for (int j = 0; j <= Order; j++)
+      for (int j = 0; j <= Orders[1]; j++)
       {
          const double sy_sz = shape_y(j)*sz;
-         for (int i = 0; i <= Order; i++, o++)
+         for (int i = 0; i <= Orders[0]; i++, o++)
          {
             sum += ( shape(o) = shape_x(i)*sy_sz*weights(o) );
          }
@@ -11009,15 +11053,15 @@ void NURBS3DFiniteElement::CalcDShape(const IntegrationPoint &ip,
    kv[2]->CalcDShape(dshape_z, ijk[2], ip.z);
 
    sum = dsum[0] = dsum[1] = dsum[2] = 0.0;
-   for (int o = 0, k = 0; k <= Order; k++)
+   for (int o = 0, k = 0; k <= Orders[2]; k++)
    {
       const double sz = shape_z(k), dsz = dshape_z(k);
-      for (int j = 0; j <= Order; j++)
+      for (int j = 0; j <= Orders[1]; j++)
       {
          const double  sy_sz  =  shape_y(j)* sz;
          const double dsy_sz  = dshape_y(j)* sz;
          const double  sy_dsz =  shape_y(j)*dsz;
-         for (int i = 0; i <= Order; i++, o++)
+         for (int i = 0; i <= Orders[0]; i++, o++)
          {
             sum += ( u(o) = shape_x(i)*sy_sz*weights(o) );
 

fem/fe.hpp

@@ -129,10 +129,12 @@ class FiniteElement
 protected:
    int Dim,      ///< Dimension of reference space
        GeomType, ///< Geometry::Type of the reference element
-       Dof,      ///< Number of degrees of freedom
-       Order,    ///< Order/degree of the shape functions
        FuncSpace, RangeType, MapType,
        DerivType, DerivRangeType, DerivMapType;
+   mutable
+   int  Dof,      ///< Number of degrees of freedom
+        Order;    ///< Order/degree of the shape functions
+   mutable int Orders[Geometry::MaxDim]; ///< Anisotropic orders
    IntegrationRule Nodes;
 #ifndef MFEM_THREAD_SAFE
    mutable DenseMatrix vshape; // Dof x Dim
@@ -201,9 +203,17 @@ public:
    /// Returns the number of degrees of freedom in the finite element
    int GetDof() const { return Dof; }
 
-   /// Returns the order of the finite element
+   /** @brief Returns the order of the finite element. In the case of
+       anisotropic orders, returns the maximum order. */
    int GetOrder() const { return Order; }
 
+   /** @brief Returns true if the FiniteElement basis *may be using* different
+       orders/degrees in different spatial directions. */
+   bool HasAnisotropicOrders() const { return Orders[0] != -1; }
+
+   /// Returns an array containing the anisotropic orders/degrees.
+   const int *GetAnisotropicOrders() const { return Orders; }
+
    /// Returns the type of space on each element
    int Space() const { return FuncSpace; }
 
@@ -2418,8 +2428,9 @@ public:
 class NURBSFiniteElement : public ScalarFiniteElement
 {
 protected:
-   mutable Array <KnotVector*> kv;
-   mutable int *ijk, patch, elem;
+   mutable Array <const KnotVector*> kv;
+   mutable const int *ijk;
+   mutable int patch, elem;
    mutable Vector weights;
 
 public:
@@ -2434,13 +2445,15 @@ public:
    }
 
    void                 Reset      ()         const { patch = elem = -1; }
-   void                 SetIJK     (int *IJK) const { ijk = IJK; }
+   void                 SetIJK     (const int *IJK) const { ijk = IJK; }
    int                  GetPatch   ()         const { return patch; }
    void                 SetPatch   (int p)    const { patch = p; }
    int                  GetElement ()         const { return elem; }
    void                 SetElement (int e)    const { elem = e; }
-   Array <KnotVector*> &KnotVectors()         const { return kv; }
+   Array <const KnotVector*> &KnotVectors()   const { return kv; }
    Vector              &Weights    ()         const { return weights; }
+   /// Update the NURBSFiniteElement according to the currently set knot vectors
+   virtual void         SetOrder   ()         const { }
 };
 
 class NURBS1DFiniteElement : public NURBSFiniteElement
@@ -2453,6 +2466,7 @@ public:
       : NURBSFiniteElement(1, Geometry::SEGMENT, p + 1, p, FunctionSpace::Qk),
         shape_x(p + 1) { }
 
+   virtual void SetOrder() const;
    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
    virtual void CalcDShape(const IntegrationPoint &ip,
                            DenseMatrix &dshape) const;
@@ -2466,9 +2480,18 @@ protected:
 public:
    NURBS2DFiniteElement(int p)
       : NURBSFiniteElement(2, Geometry::SQUARE, (p + 1)*(p + 1), p,
-                           FunctionSpace::Qk), u(Dof),
-        shape_x(p + 1), shape_y(p + 1), dshape_x(p + 1), dshape_y(p + 1) { }
+                           FunctionSpace::Qk),
+        u(Dof), shape_x(p + 1), shape_y(p + 1), dshape_x(p + 1), dshape_y(p + 1)
+   { Orders[0] = Orders[1] = p; }
 
+   NURBS2DFiniteElement(int px, int py)
+      : NURBSFiniteElement(2, Geometry::SQUARE, (px + 1)*(py + 1),
+                           std::max(px, py), FunctionSpace::Qk),
+        u(Dof), shape_x(px + 1), shape_y(py + 1), dshape_x(px + 1),
+        dshape_y(py + 1)
+   { Orders[0] = px; Orders[1] = py; }
+
+   virtual void SetOrder() const;
    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
    virtual void CalcDShape(const IntegrationPoint &ip,
                            DenseMatrix &dshape) const;
@@ -2482,10 +2505,19 @@ protected:
 public:
    NURBS3DFiniteElement(int p)
       : NURBSFiniteElement(3, Geometry::CUBE, (p + 1)*(p + 1)*(p + 1), p,
-                           FunctionSpace::Qk), u(Dof),
-        shape_x(p + 1), shape_y(p + 1), shape_z(p + 1),
-        dshape_x(p + 1), dshape_y(p + 1), dshape_z(p + 1) { }
+                           FunctionSpace::Qk),
+        u(Dof), shape_x(p + 1), shape_y(p + 1), shape_z(p + 1),
+        dshape_x(p + 1), dshape_y(p + 1), dshape_z(p + 1)
+   { Orders[0] = Orders[1] = Orders[2] = p; }
 
+   NURBS3DFiniteElement(int px, int py, int pz)
+      : NURBSFiniteElement(3, Geometry::CUBE, (px + 1)*(py + 1)*(pz + 1),
+                           std::max(std::max(px,py),pz), FunctionSpace::Qk),
+        u(Dof), shape_x(px + 1), shape_y(py + 1), shape_z(pz + 1),
+        dshape_x(px + 1), dshape_y(py + 1), dshape_z(pz + 1)
+   { Orders[0] = px; Orders[1] = py; Orders[2] = pz; }
+
+   virtual void SetOrder() const;
    virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
    virtual void CalcDShape(const IntegrationPoint &ip,
                            DenseMatrix &dshape) const;

fem/fe_coll.cpp

@@ -260,12 +260,20 @@ FiniteElementCollection *FiniteElementCollection::New(const char *name)
    }
    else if (!strncmp(name, "NURBS", 5))
    {
-      fec = new NURBSFECollection(atoi(name + 5));
+      if (name[5] != '\0')
+      {
+         // "NURBS" + "number" --> fixed order nurbs collection
+         fec = new NURBSFECollection(atoi(name + 5));
+      }
+      else
+      {
+         // "NURBS" --> variable order nurbs collection
+         fec = new NURBSFECollection();
+      }
    }
    else
    {
-      mfem_error("FiniteElementCollection::New : "
-                 "Unknown FiniteElementCollection!");
+      MFEM_ABORT("unknown FiniteElementCollection: " << name);
    }
    MFEM_VERIFY(!strcmp(fec->Name(), name), "input name: \"" << name
                << "\" does not match the created collection name: \""
@@ -2418,16 +2426,30 @@ Local_FECollection::Local_FECollection(const char *fe_name)
 }
 
 
-void NURBSFECollection::Allocate(int Order)
+NURBSFECollection::NURBSFECollection(int Order)
 {
-   SegmentFE        = new NURBS1DFiniteElement(Order);
-   QuadrilateralFE  = new NURBS2DFiniteElement(Order);
-   ParallelepipedFE = new NURBS3DFiniteElement(Order);
+   const int order = (Order == VariableOrder) ? 1 : Order;
+   SegmentFE        = new NURBS1DFiniteElement(order);
+   QuadrilateralFE  = new NURBS2DFiniteElement(order);
+   ParallelepipedFE = new NURBS3DFiniteElement(order);
 
-   snprintf(name, 16, "NURBS%i", Order);
+   SetOrder(Order);
 }
 
-void NURBSFECollection::Deallocate()
+void NURBSFECollection::SetOrder(int Order) const
+{
+   mOrder = Order;
+   if (Order != VariableOrder)
+   {
+      snprintf(name, 16, "NURBS%i", Order);
+   }
+   else
+   {
+      snprintf(name, 16, "NURBS");
+   }
+}
+
+NURBSFECollection::~NURBSFECollection()
 {
    delete ParallelepipedFE;
    delete QuadrilateralFE;

fem/fe_coll.hpp

@@ -270,18 +270,20 @@ private:
    NURBS2DFiniteElement *QuadrilateralFE;
    NURBS3DFiniteElement *ParallelepipedFE;
 
-   char name[16];
-
-   void Allocate(int Order);
-   void Deallocate();
+   mutable int mOrder; // >= 1 or VariableOrder
+   // The 'name' can be:
+   // 1) name = "NURBS" + "number", for fixed order, or
+   // 2) name = "NURBS", for VariableOrder.
+   // The name is updated before writing it to a stream, for example, see
+   // FiniteElementSpace::Save().
+   mutable char name[16];
 
 public:
-   explicit NURBSFECollection(int Order) { Allocate(Order); }
+   enum { VariableOrder = -1 };
 
-   int GetOrder() const { return SegmentFE->GetOrder(); }
-
-   /// Change the order of the collection
-   void UpdateOrder(int Order) { Deallocate(); Allocate(Order); }
+   /** @brief The parameter @a Order must be either a positive number, for fixed
+      order, or VariableOrder (default). */
+   explicit NURBSFECollection(int Order = VariableOrder);
 
    void Reset() const
    {
@@ -290,6 +292,14 @@ public:
       ParallelepipedFE->Reset();
    }
 
+   /** @brief Get the order of the NURBS collection: either a positive number,
+       when using fixed order, or VariableOrder. */
+   int GetOrder() const { return mOrder; }
+
+   /** @brief Set the order and the name, based on the given @a Order: either a
+       positive number for fixed order, or VariableOrder. */
+   void SetOrder(int Order) const;
+
    virtual const FiniteElement *
    FiniteElementForGeometry(int GeomType) const;
 
@@ -301,7 +311,7 @@ public:
 
    FiniteElementCollection *GetTraceCollection() const;
 
-   virtual ~NURBSFECollection() { Deallocate(); }
+   virtual ~NURBSFECollection();
 };
 
 

fem/fespace.cpp

@@ -11,6 +11,7 @@
 
 // Implementation of FiniteElementSpace
 
+#include "../general/text.hpp"
 #include "../mesh/mesh_headers.hpp"
 #include "fem.hpp"
 
@@ -53,6 +54,43 @@ DofsToVDofs<Ordering::byVDIM>(int ndofs, int vdim, Array<int> &dofs)
    }
 }
 
+
+FiniteElementSpace::FiniteElementSpace()
+   : mesh(NULL), fec(NULL), vdim(0), ordering(Ordering::byNODES),
+     ndofs(0), nvdofs(0), nedofs(0), nfdofs(0), nbdofs(0),
+     fdofs(NULL), bdofs(NULL),
+     elem_dof(NULL), bdrElem_dof(NULL),
+     NURBSext(NULL), own_ext(false),
+     cP(NULL), cR(NULL), cP_is_set(false),
+     T(NULL), own_T(false),
+     sequence(0)
+{ }
+
+FiniteElementSpace::FiniteElementSpace(const FiniteElementSpace &orig,
+                                       Mesh *mesh,
+                                       const FiniteElementCollection *fec)
+{
+   mesh = mesh ? mesh : orig.mesh;
+   fec = fec ? fec : orig.fec;
+   NURBSExtension *NURBSext = NULL;
+   if (orig.NURBSext && orig.NURBSext != orig.mesh->NURBSext)
+   {
+#ifdef MFEM_USE_MPI
+      ParNURBSExtension *pNURBSext =
+         dynamic_cast<ParNURBSExtension *>(orig.NURBSext);
+      if (pNURBSext)
+      {
+         NURBSext = new ParNURBSExtension(*pNURBSext);
+      }
+      else
+#endif
+      {
+         NURBSext = new NURBSExtension(*orig.NURBSext);
+      }
+   }
+   Constructor(mesh, NURBSext, fec, orig.vdim, orig.ordering);
+}
+
 int FiniteElementSpace::GetOrder(int i) const
 {
    int GeomType = mesh->GetElementBaseGeometry(i);
@@ -861,9 +899,9 @@ void FiniteElementSpace::GetLocalDerefinementMatrices(
          {
             IntegrationPoint ip;
             ip.Set(pt, dim);
-            fe->CalcShape(ip, shape); // TODO: H(curl), etc.?
             MFEM_ASSERT(dynamic_cast<const NodalFiniteElement*>(fe),
                         "only nodal FEs are implemented");
+            fe->CalcShape(ip, shape); // TODO: H(curl), etc.?
             lR.SetRow(j, shape);
          }
       }
@@ -927,9 +965,9 @@ SparseMatrix* FiniteElementSpace::DerefinementMatrix(int old_ndofs,
    return R;
 }
 
-FiniteElementSpace::FiniteElementSpace(Mesh *mesh,
-                                       const FiniteElementCollection *fec,
-                                       int vdim, int ordering)
+void FiniteElementSpace::Constructor(Mesh *mesh, NURBSExtension *NURBSext,
+                                     const FiniteElementCollection *fec,
+                                     int vdim, int ordering)
 {
    this->mesh = mesh;
    this->fec = fec;
@@ -948,33 +986,29 @@ FiniteElementSpace::FiniteElementSpace(Mesh *mesh,
          mfem_error("FiniteElementSpace::FiniteElementSpace :\n"
                     "   NURBS FE space requires NURBS mesh.");
       }
+
+      if (NURBSext == NULL)
+      {
+         this->NURBSext = mesh->NURBSext;
+         own_ext = 0;
+      }
       else
       {
-         int Order = nurbs_fec->GetOrder();
-         if (mesh->NURBSext->GetOrder() == Order)
-         {
-            NURBSext = mesh->NURBSext;
-            own_ext = 0;
-         }
-         else
-         {
-            NURBSext = new NURBSExtension(mesh->NURBSext, Order);
-            own_ext = 1;
-         }
-         UpdateNURBS();
-         cP = cR = NULL;
-         cP_is_set = false;
-         T = NULL;
-         own_T = true;
+         this->NURBSext = NURBSext;
+         own_ext = 1;
       }
+      UpdateNURBS();
+      cP = cR = NULL;
+      cP_is_set = false;
+      T = NULL;
+      own_T = true;
    }
    else
    {
-      NURBSext = NULL;
+      this->NURBSext = NULL;
       own_ext = 0;
       Construct();
    }
-
    BuildElementToDofTable();
 }
 
@@ -1007,7 +1041,8 @@ void FiniteElementSpace::UpdateNURBS()
 
 void FiniteElementSpace::Construct()
 {
-   int i;
+   // This method should be used only for non-NURBS spaces.
+   MFEM_ASSERT(!NURBSext, "internal error");
 
    elem_dof = NULL;
    bdrElem_dof = NULL;
@@ -1046,7 +1081,7 @@ void FiniteElementSpace::Construct()
       {
          fdofs = new int[mesh->GetNFaces()+1];
          fdofs[0] = 0;
-         for (i = 0; i < mesh->GetNFaces(); i++)
+         for (int i = 0; i < mesh->GetNFaces(); i++)
          {
             nfdofs += fdof;
             // nfdofs += fec->DofForGeometry(mesh->GetFaceBaseGeometry(i));
@@ -1059,7 +1094,7 @@ void FiniteElementSpace::Construct()
    {
       bdofs = new int[mesh->GetNE()+1];
       bdofs[0] = 0;
-      for (i = 0; i < mesh->GetNE(); i++)
+      for (int i = 0; i < mesh->GetNE(); i++)
       {
          int geom = mesh->GetElementBaseGeometry(i);
          nbdofs += fec->DofForGeometry(geom);
@@ -1565,10 +1600,157 @@ void FiniteElementSpace::Update(bool want_transform)
 
 void FiniteElementSpace::Save(std::ostream &out) const
 {
-   out << "FiniteElementSpace\n"
+   int fes_format = 90; // the original format, v0.9
+   bool nurbs_unit_weights = false;
+
+   // Determine the format that should be used.
+   if (!NURBSext)
+   {
+      // TODO: if this is a variable-order FE space, use fes_format = 100.
+   }
+   else
+   {
+      const NURBSFECollection *nurbs_fec =
+         dynamic_cast<const NURBSFECollection *>(fec);
+      MFEM_VERIFY(nurbs_fec, "invalid FE collection");
+      nurbs_fec->SetOrder(NURBSext->GetOrder());
+      const double eps = 5e-14;
+      nurbs_unit_weights = (NURBSext->GetWeights().Min() >= 1.0-eps &&
+                            NURBSext->GetWeights().Max() <= 1.0+eps);
+      if (NURBSext->GetOrder() == NURBSFECollection::VariableOrder ||
+          (NURBSext != mesh->NURBSext && !nurbs_unit_weights))
+      {
+         fes_format = 100; // v1.0 format
+      }
+   }
+
+   out << (fes_format == 90 ?
+           "FiniteElementSpace\n" : "MFEM FiniteElementSpace v1.0\n")
        << "FiniteElementCollection: " << fec->Name() << '\n'
        << "VDim: " << vdim << '\n'
        << "Ordering: " << ordering << '\n';
+
+   if (fes_format == 100) // v1.0
+   {
+      if (!NURBSext)
+      {
+         // TODO: this is a variable-order FE space --> write 'element_orders'.
+      }
+      else if (NURBSext != mesh->NURBSext)
+      {
+         if (NURBSext->GetOrder() != NURBSFECollection::VariableOrder)
+         {
+            out << "NURBS_order\n" << NURBSext->GetOrder() << '\n';
+         }
+         else
+         {
+            out << "NURBS_orders\n";
+            // 1 = do not write the size, just the entries:
+            NURBSext->GetOrders().Save(out, 1);
+         }
+         // If the weights are not unit, write them to the output:
+         if (!nurbs_unit_weights)
+         {
+            out << "NURBS_weights\n";
+            NURBSext->GetWeights().Print(out, 1);
+         }
+      }
+      out << "End: MFEM FiniteElementSpace v1.0\n";
+   }
+}
+
+FiniteElementCollection *FiniteElementSpace::Load(Mesh *m, std::istream &input)
+{
+   string buff;
+   int fes_format = 0, ord;
+   FiniteElementCollection *r_fec;
+
+   Destroy();
+
+   input >> std::ws;
+   getline(input, buff);  // 'FiniteElementSpace'
+   filter_dos(buff);
+   if (buff == "FiniteElementSpace") { fes_format = 90; /* v0.9 */ }
+   else if (buff == "MFEM FiniteElementSpace v1.0") { fes_format = 100; }
+   else { MFEM_ABORT("input stream is not a FiniteElementSpace!"); }
+   getline(input, buff, ' '); // 'FiniteElementCollection:'
+   input >> std::ws;
+   getline(input, buff);
+   filter_dos(buff);
+   r_fec = FiniteElementCollection::New(buff.c_str());
+   getline(input, buff, ' '); // 'VDim:'
+   input >> vdim;
+   getline(input, buff, ' '); // 'Ordering:'
+   input >> ord;
+
+   NURBSFECollection *nurbs_fec = dynamic_cast<NURBSFECollection*>(r_fec);
+   NURBSExtension *NURBSext = NULL;
+   if (fes_format == 90) // original format, v0.9
+   {
+      if (nurbs_fec)
+      {
+         MFEM_VERIFY(m->NURBSext, "NURBS FE collection requires a NURBS mesh!");
+         const int order = nurbs_fec->GetOrder();
+         if (order != m->NURBSext->GetOrder() &&
+             order != NURBSFECollection::VariableOrder)
+         {
+            NURBSext = new NURBSExtension(m->NURBSext, order);
+         }
+      }
+   }
+   else if (fes_format == 100) // v1.0
+   {
+      while (1)
+      {
+         skip_comment_lines(input, '#');
+         MFEM_VERIFY(input.good(), "error reading FiniteElementSpace v1.0");
+         getline(input, buff);
+         filter_dos(buff);
+         if (buff == "NURBS_order" || buff == "NURBS_orders")
+         {
+            MFEM_VERIFY(nurbs_fec,
+                        buff << ": NURBS FE collection is required!");
+            MFEM_VERIFY(m->NURBSext, buff << ": NURBS mesh is required!");
+            MFEM_VERIFY(!NURBSext, buff << ": order redefinition!");
+            if (buff == "NURBS_order")
+            {
+               int order;
+               input >> order;
+               NURBSext = new NURBSExtension(m->NURBSext, order);
+            }
+            else
+            {
+               Array<int> orders;
+               orders.Load(m->NURBSext->GetNKV(), input);
+               NURBSext = new NURBSExtension(m->NURBSext, orders);
+            }
+         }
+         else if (buff == "NURBS_weights")
+         {
+            MFEM_VERIFY(NURBSext, "NURBS_weights: NURBS_orders have to be "
+                        "specified before NURBS_weights!");
+            NURBSext->GetWeights().Load(input, NURBSext->GetNDof());
+         }
+         else if (buff == "element_orders")
+         {
+            MFEM_VERIFY(!nurbs_fec, "section element_orders cannot be used "
+                        "with a NURBS FE collection");
+            MFEM_ABORT("element_orders: not implemented yet!");
+         }
+         else if (buff == "End: MFEM FiniteElementSpace v1.0")
+         {
+            break;
+         }
+         else
+         {
+            MFEM_ABORT("unknown section: " << buff);
+         }
+      }
+   }
+
+   Constructor(m, NURBSext, r_fec, vdim, ord);
+
+   return r_fec;
 }
 
 

fem/fespace.hpp

@@ -21,15 +21,21 @@
 namespace mfem
 {
 
-/** The ordering method used when the number of unknowns per mesh
-    node (vector dimension) is bigger than 1. */
+/** @brief The ordering method used when the number of unknowns per mesh node
+    (vector dimension) is bigger than 1. */
 class Ordering
 {
 public:
-   /** Ordering methods:
-       byNODES - loop first over the nodes (inner loop) then over the vector dimension (outer loop),
-       byVDIM  - loop first over the vector dimension (inner loop) then over the nodes (outer loop)  */
-   enum Type { byNODES, byVDIM };
+   /// %Ordering methods:
+   enum Type
+   {
+      byNODES, /**< loop first over the nodes (inner loop) then over the vector
+                    dimension (outer loop); symbolically it can be represented
+                    as: XXX...,YYY...,ZZZ... */
+      byVDIM   /**< loop first over the vector dimension (inner loop) then over
+                    the nodes (outer loop); symbolically it can be represented
+                    as: XYZ,XYZ,XYZ,... */
+   };
 
    template <Type Ord>
    static inline int Map(int ndofs, int vdim, int dof, int vd);
@@ -55,32 +61,33 @@ Ordering::Map<Ordering::byVDIM>(int ndofs, int vdim, int dof, int vd)
 
 class NURBSExtension;
 
-/** Class FiniteElementSpace - responsible for providing FEM view of the mesh
-    (mainly managing the set of degrees of freedom). */
+/** @brief Class FiniteElementSpace - responsible for providing FEM view of the
+    mesh, mainly managing the set of degrees of freedom. */
 class FiniteElementSpace
 {
 protected:
-   /// The mesh that FE space lives on.
+   /// The mesh that FE space lives on (not owned).
    Mesh *mesh;
 
+   /// Associated FE collection (not owned).
    const FiniteElementCollection *fec;
 
-   /// Vector dimension (number of unknowns per degree of freedom).
+   /// %Vector dimension (number of unknowns per degree of freedom).
    int vdim;
 
-   /** Type of ordering of dofs.
-       Ordering::byNODES - first nodes, then vector dimension,
-       Ordering::byVDIM  - first vector dimension, then nodes  */
+   /** Type of ordering of the vector dofs when #vdim > 1.
+       - Ordering::byNODES - first nodes, then vector dimension,
+       - Ordering::byVDIM  - first vector dimension, then nodes */
    Ordering::Type ordering;
 
-   /// Number of degrees of freedom. Number of unknowns are ndofs*vdim.
+   /// Number of degrees of freedom. Number of unknowns is #ndofs * #vdim.
    int ndofs;
 
    int nvdofs, nedofs, nfdofs, nbdofs;
    int *fdofs, *bdofs;
 
-   mutable Table *elem_dof;
-   Table *bdrElem_dof;
+   mutable Table *elem_dof; // if NURBS FE space, not owned; otherwise, owned.
+   Table *bdrElem_dof; // used only with NURBS FE spaces; not owned.
 
    Array<int> dof_elem_array, dof_ldof_array;
 
@@ -90,9 +97,9 @@ protected:
    /** Matrix representing the prolongation from the global conforming dofs to
        a set of intermediate partially conforming dofs, e.g. the dofs associated
        with a "cut" space on a non-conforming mesh. */
-   mutable SparseMatrix *cP;
+   mutable SparseMatrix *cP; // owned
    /// Conforming restriction matrix such that cR.cP=I.
-   mutable SparseMatrix *cR;
+   mutable SparseMatrix *cR; // owned
    mutable bool cP_is_set;
 
    /// Transformation to apply to GridFunctions after space Update().
@@ -127,14 +134,51 @@ protected:
    /// Calculate GridFunction restriction matrix after mesh derefinement.
    SparseMatrix* DerefinementMatrix(int old_ndofs, const Table* old_elem_dof);
 
+   /// Help function for constructors.
+   void Constructor(Mesh *mesh, NURBSExtension *ext,
+                    const FiniteElementCollection *fec,
+                    int vdim = 1, int ordering = Ordering::byNODES);
 
 public:
-   FiniteElementSpace(Mesh *mesh, const FiniteElementCollection *fec,
-                      int vdim = 1, int ordering = Ordering::byNODES);
+   /** @brief Default constructor: the object is invalid until initialized using
+       the method Load(). */
+   FiniteElementSpace();
+
+   /** @brief Copy constructor: deep copy all data from @a orig except the Mesh,
+       the FiniteElementCollection, ans some derived data. */
+   /** If the @a mesh or @a fec poiters are NULL (default), then the new
+       FiniteElementSpace will reuse the respective pointers from @a orig. If
+       any of these pointers is not NULL, the given pointer will be used instead
+       of the one used by @a orig.
+
+       @note The objects pointed to by the @a mesh and @a fec parameters must be
+       either the same objects as the ones used by @a orig, or copies of them.
+       Otherwise, the behavior is undefined.
+
+       @note Derived data objects, such as the conforming prolongation and
+       restriction matrices, and the update operator, will not be copied, even
+       if they are created in the @a orig object. */
+   FiniteElementSpace(const FiniteElementSpace &orig, Mesh *mesh = NULL,
+                      const FiniteElementCollection *fec = NULL);
+
+   FiniteElementSpace(Mesh *mesh,
+                      const FiniteElementCollection *fec,
+                      int vdim = 1, int ordering = Ordering::byNODES)
+   { Constructor(mesh, NULL, fec, vdim, ordering); }
+
+   /// Construct a NURBS FE space based on the given NURBSExtension, @a ext.
+   /** @note If the pointer @a ext is NULL, this constructor is equivalent to
+       the standard constructor with the same arguments minus the
+       NURBSExtension, @a ext. */
+   FiniteElementSpace(Mesh *mesh, NURBSExtension *ext,
+                      const FiniteElementCollection *fec,
+                      int vdim = 1, int ordering = Ordering::byNODES)
+   { Constructor(mesh, ext, fec, vdim, ordering); }
 
    /// Returns the mesh
    inline Mesh *GetMesh() const { return mesh; }
 
+   const NURBSExtension *GetNURBSext() const { return NURBSext; }
    NURBSExtension *GetNURBSext() { return NURBSext; }
    NURBSExtension *StealNURBSext();
 
@@ -160,6 +204,7 @@ public:
    /// Returns number of degrees of freedom.
    inline int GetNDofs() const { return ndofs; }
 
+   /// Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
    inline int GetVSize() const { return vdim * ndofs; }
 
    /// Return the number of vector true (conforming) dofs.
@@ -377,6 +422,10 @@ public:
 
    void Save(std::ostream &out) const;
 
+   /** @brief Read a FiniteElementSpace from a stream. The returned
+       FiniteElementCollection is owned by the caller. */
+   FiniteElementCollection *Load(Mesh *m, std::istream &input);
+
    virtual ~FiniteElementSpace();
 };
 

fem/geom.hpp

@@ -32,6 +32,7 @@ public:
    enum Type { POINT, SEGMENT, TRIANGLE, SQUARE, TETRAHEDRON, CUBE };
 
    static const int NumGeom = 6;
+   static const int MaxDim = 3;
    static const int NumBdrArray[NumGeom];
    static const char *Name[NumGeom];
    static const double Volume[NumGeom];

fem/gridfunc.cpp

@@ -30,31 +30,32 @@ using namespace std;
 GridFunction::GridFunction(Mesh *m, std::istream &input)
    : Vector()
 {
-   std::string buff;
-   int vdim;
+   fes = new FiniteElementSpace;
+   fec = fes->Load(m, input);
 
-   input >> std::ws;
-   getline(input, buff);  // 'FiniteElementSpace'
-   filter_dos(buff);
-   if (buff != "FiniteElementSpace")
+   skip_comment_lines(input, '#');
+   istream::int_type next_char = input.peek();
+   if (next_char == 'N') // First letter of "NURBS_patches"
    {
-      mfem_error("GridFunction::GridFunction():"
-                 " input stream is not a GridFunction!");
+      string buff;
+      getline(input, buff);
+      filter_dos(buff);
+      if (buff == "NURBS_patches")
+      {
+         MFEM_VERIFY(fes->GetNURBSext(),
+                     "NURBS_patches requires NURBS FE space");
+         fes->GetNURBSext()->LoadSolution(input, *this);
+      }
+      else
+      {
+         MFEM_ABORT("unknown section: " << buff);
+      }
    }
-   getline(input, buff, ' '); // 'FiniteElementCollection:'
-   input >> std::ws;
-   getline(input, buff);
-   filter_dos(buff);
-   fec = FiniteElementCollection::New(buff.c_str());
-   getline(input, buff, ' '); // 'VDim:'
-   input >> vdim;
-   getline(input, buff, ' '); // 'Ordering:'
-   int ordering;
-   input >> ordering;
-   getline(input, buff); // read the empty line
-   fes = new FiniteElementSpace(m, fec, vdim, ordering);
-   Vector::Load(input, fes->GetVSize());
-   sequence = 0;
+   else
+   {
+      Vector::Load(input, fes->GetVSize());
+   }
+   sequence = fes->GetSequence();
 }
 
 GridFunction::GridFunction(Mesh *m, GridFunction *gf_array[], int num_pieces)
@@ -2260,6 +2261,16 @@ void GridFunction::Save(std::ostream &out) const
 {
    fes->Save(out);
    out << '\n';
+#if 0
+   // Testing: write NURBS GridFunctions using "NURBS_patches" format.
+   if (fes->GetNURBSext())
+   {
+      out << "NURBS_patches\n";
+      fes->GetNURBSext()->PrintSolution(*this, out);
+      out.flush();
+      return;
+   }
+#endif
    if (fes->GetOrdering() == Ordering::byNODES)
    {
       Vector::Print(out, 1);

fem/pfespace.cpp

@@ -24,34 +24,99 @@
 namespace mfem
 {
 
+ParFiniteElementSpace::ParFiniteElementSpace(
+   const ParFiniteElementSpace &orig, ParMesh *pmesh,
+   const FiniteElementCollection *fec)
+   : FiniteElementSpace(orig, pmesh, fec)
+{
+   ParInit(pmesh ? pmesh : orig.pmesh);
+}
+
+ParFiniteElementSpace::ParFiniteElementSpace(
+   const FiniteElementSpace &orig, ParMesh &pmesh,
+   const FiniteElementCollection *fec)
+   : FiniteElementSpace(orig, &pmesh, fec)
+{
+   ParInit(&pmesh);
+}
+
+ParFiniteElementSpace::ParFiniteElementSpace(
+   ParMesh *pm, const FiniteElementSpace *global_fes, const int *partitioning,
+   const FiniteElementCollection *f)
+   : FiniteElementSpace(pm, MakeLocalNURBSext(global_fes->GetNURBSext(),
+                                              pm->NURBSext),
+                        f ? f : global_fes->FEColl(),
+                        global_fes->GetVDim(), global_fes->GetOrdering())
+{
+   ParInit(pm);
+   // For NURBS spaces, the variable-order data is contained in the
+   // NURBSExtension of 'global_fes' and inside the ParNURBSExtension of 'pm'.
+
+   // TODO: when general variable-order support is added, copy the local portion
+   // of the variable-oder data from 'global_fes' to 'this'.
+}
+
 ParFiniteElementSpace::ParFiniteElementSpace(
    ParMesh *pm, const FiniteElementCollection *f, int dim, int ordering)
    : FiniteElementSpace(pm, f, dim, ordering)
 {
-   mesh = pmesh = pm;
+   ParInit(pm);
+}
 
+ParFiniteElementSpace::ParFiniteElementSpace(
+   ParMesh *pm, NURBSExtension *ext, const FiniteElementCollection *f,
+   int dim, int ordering)
+   : FiniteElementSpace(pm, ext, f, dim, ordering)
+{
+   ParInit(pm);
+}
+
+// static method
+ParNURBSExtension *ParFiniteElementSpace::MakeLocalNURBSext(
+   const NURBSExtension *globNURBSext, const NURBSExtension *parNURBSext)
+{
+   if (globNURBSext == NULL) { return NULL; }
+   const ParNURBSExtension *pNURBSext =
+      dynamic_cast<const ParNURBSExtension*>(parNURBSext);
+   MFEM_ASSERT(pNURBSext, "need a ParNURBSExtension");
+   // make a copy of globNURBSext:
+   NURBSExtension *tmp_globNURBSext = new NURBSExtension(*globNURBSext);
+   // tmp_globNURBSext will be deleted by the following ParNURBSExtension ctor:
+   return new ParNURBSExtension(tmp_globNURBSext, pNURBSext);
+}
+
+void ParFiniteElementSpace::ParInit(ParMesh *pm)
+{
+   pmesh = pm;
    MyComm = pmesh->GetComm();
-   MPI_Comm_size(MyComm, &NRanks);
-   MPI_Comm_rank(MyComm, &MyRank);
-
-   num_face_nbr_dofs = -1;
-
+   NRanks = pmesh->GetNRanks();
+   MyRank = pmesh->GetMyRank();
+   gcomm = NULL;
    P = NULL;
    Pconf = NULL;
    R = NULL;
-   gcomm = NULL;
+   num_face_nbr_dofs = -1;
+
+   if (NURBSext && !pNURBSext())
+   {
+      // This is necessary in some cases: e.g. when the FiniteElementSpace
+      // constructor creates a serial NURBSExtension of higher order than the
+      // mesh NURBSExtension.
+      MFEM_ASSERT(own_ext, "internal error");
+
+      ParNURBSExtension *pNe = new ParNURBSExtension(
+         NURBSext, dynamic_cast<ParNURBSExtension *>(pmesh->NURBSext));
+      // serial NURBSext is destroyed by the above constructor
+      NURBSext = pNe;
+      UpdateNURBS();
+   }
 
    Construct();
 
    // Apply the ldof_signs to the elem_dof Table
    if (Conforming() && !NURBSext)
    {
-      Array<int> dofs;
-      for (int i = 0; i < elem_dof->Size(); i++)
-      {
-         dofs.MakeRef(elem_dof->GetRow(i), elem_dof->RowSize(i));
-         ApplyLDofSigns(dofs);
-      }
+      ApplyLDofSigns(*elem_dof);
    }
 }
 
@@ -59,17 +124,6 @@ void ParFiniteElementSpace::Construct()
 {
    if (NURBSext)
    {
-      if (own_ext)
-      {
-         // the FiniteElementSpace constructor created a serial
-         // NURBSExtension of higher order than the mesh NURBSExtension
-
-         ParNURBSExtension *pNe = new ParNURBSExtension(
-            NURBSext, dynamic_cast<ParNURBSExtension *>(pmesh->NURBSext));
-         // serial NURBSext is destroyed by the above constructor
-         NURBSext = pNe;
-         UpdateNURBS();
-      }
       ConstructTrueNURBSDofs();
       GenerateGlobalOffsets();
    }
@@ -98,6 +152,7 @@ void ParFiniteElementSpace::GetGroupComm(
 
    nvd = fec->DofForGeometry(Geometry::POINT);
    ned = fec->DofForGeometry(Geometry::SEGMENT);
+   // Assuming all faces are the same type:
    nfd = (fdofs) ? (fdofs[1]-fdofs[0]) : (0);
 
    if (ldof_sign)
@@ -263,6 +318,12 @@ void ParFiniteElementSpace::ApplyLDofSigns(Array<int> &dofs) const
    }
 }
 
+void ParFiniteElementSpace::ApplyLDofSigns(Table &el_dof) const
+{
+   Array<int> all_dofs(el_dof.GetJ(), el_dof.Size_of_connections());
+   ApplyLDofSigns(all_dofs);
+}
+
 void ParFiniteElementSpace::GetElementDofs(int i, Array<int> &dofs) const
 {
    if (elem_dof)

fem/pfespace.hpp

@@ -35,10 +35,10 @@ private:
    MPI_Comm MyComm;
    int NRanks, MyRank;
 
-   /// Parallel mesh.
+   /// Parallel mesh; not owned.
    ParMesh *pmesh;
 
-   /// GroupCommunicator on the local VDofs
+   /// GroupCommunicator on the local VDofs; owned.
    GroupCommunicator *gcomm;
 
    /// Number of true dofs in this processor (local true dofs).
@@ -65,19 +65,25 @@ private:
    /// The sign of the basis functions at the scalar local dofs.
    Array<int> ldof_sign;
 
-   /// The matrix P (interpolation from true dof to dof).
+   /// The matrix P (interpolation from true dof to dof); owned.
    mutable HypreParMatrix *P;
    ConformingProlongationOperator *Pconf;
 
-   /// The (block-diagonal) matrix R (restriction of dof to true dof)
+   /// The (block-diagonal) matrix R (restriction of dof to true dof); owned.
    mutable SparseMatrix *R;
 
    ParNURBSExtension *pNURBSext() const
    { return dynamic_cast<ParNURBSExtension *>(NURBSext); }
 
+   static ParNURBSExtension *MakeLocalNURBSext(
+      const NURBSExtension *globNURBSext, const NURBSExtension *parNURBSext);
+
    GroupTopology &GetGroupTopo() const
    { return (NURBSext) ? pNURBSext()->gtopo : pmesh->gtopo; }
 
+   // Auxiliary method used in constructors
+   void ParInit(ParMesh *pm);
+
    void Construct();
    void Destroy();
 
@@ -93,6 +99,7 @@ private:
    void ConstructTrueNURBSDofs();
 
    void ApplyLDofSigns(Array<int> &dofs) const;
+   void ApplyLDofSigns(Table &el_dof) const;
 
    /// Helper struct to store DOF dependencies in a parallel NC mesh.
    struct Dependency
@@ -160,9 +167,57 @@ public:
    // Local face-neighbor data: face-neighbor to ldof
    Table send_face_nbr_ldof;
 
+   /** @brief Copy constructor: deep copy all data from @a orig except the
+       ParMesh, the FiniteElementCollection, and some derived data. */
+   /** If the @a pmesh or @a fec poiters are NULL (default), then the new
+       ParFiniteElementSpace will reuse the respective pointers from @a orig. If
+       any of these pointers is not NULL, the given pointer will be used instead
+       of the one used by @a orig.
+
+       @note The objects pointed to by the @a pmesh and @a fec parameters must
+       be either the same objects as the ones used by @a orig, or copies of
+       them. Otherwise, the behavior is undefined.
+
+       @note Derived data objects, such as the parallel prolongation and
+       restriction operators, the update operator, and any of the face-neighbor
+       data, will not be copied, even if they are created in the @a orig object.
+   */
+   ParFiniteElementSpace(const ParFiniteElementSpace &orig,
+                         ParMesh *pmesh = NULL,
+                         const FiniteElementCollection *fec = NULL);
+
+   /** @brief Convert/copy the *local* (Par)FiniteElementSpace @a orig to
+       ParFiniteElementSpace: deep copy all data from @a orig except the Mesh,
+       the FiniteElementCollection, and some derived data. */
+   ParFiniteElementSpace(const FiniteElementSpace &orig, ParMesh &pmesh,
+                         const FiniteElementCollection *fec = NULL);
+
+   /** @brief Construct the *local* ParFiniteElementSpace corresponing to the
+       global FE space, @a global_fes. */
+   /** The parameter @a pm is the *local* ParMesh obtained by decomposing the
+       global Mesh used by @a global_fes. The array @a partitioning represents
+       the parallel decomposition - it maps global element ids to MPI ranks.
+       If the FiniteElementCollection, @a f, is NULL (default), the FE
+       collection used by @a global_fes will be reused. If @a f is not NULL, it
+       must be the same as, or a copy of, the FE collection used by
+       @a global_fes. */
+   ParFiniteElementSpace(ParMesh *pm, const FiniteElementSpace *global_fes,
+                         const int *partitioning,
+                         const FiniteElementCollection *f = NULL);
+
    ParFiniteElementSpace(ParMesh *pm, const FiniteElementCollection *f,
                          int dim = 1, int ordering = Ordering::byNODES);
 
+   /// Construct a NURBS FE space based on the given NURBSExtension, @a ext.
+   /** The parameter @a ext will be deleted by this constructor, replaced by a
+       ParNURBSExtension owned by the ParFiniteElementSpace.
+       @note If the pointer @a ext is NULL, this constructor is equivalent to
+       the standard constructor with the same arguments minus the
+       NURBSExtension, @a ext. */
+   ParFiniteElementSpace(ParMesh *pm, NURBSExtension *ext,
+                         const FiniteElementCollection *f,
+                         int dim = 1, int ordering = Ordering::byNODES);
+
    MPI_Comm GetComm() const { return MyComm; }
    int GetNRanks() const { return NRanks; }
    int GetMyRank() const { return MyRank; }

fem/pgridfunc.cpp

@@ -33,32 +33,37 @@ ParGridFunction::ParGridFunction(ParFiniteElementSpace *pf, HypreParVector *tv)
    Distribute(tv);
 }
 
-ParGridFunction::ParGridFunction(ParMesh *pmesh, GridFunction *gf,
-                                 int * partitioning)
+ParGridFunction::ParGridFunction(ParMesh *pmesh, const GridFunction *gf,
+                                 const int *partitioning)
 {
+   const FiniteElementSpace *glob_fes = gf->FESpace();
    // duplicate the FiniteElementCollection from 'gf'
-   fec = FiniteElementCollection::New(gf->FESpace()->FEColl()->Name());
-   fes = pfes = new ParFiniteElementSpace(pmesh, fec, gf->FESpace()->GetVDim(),
-                                          gf->FESpace()->GetOrdering());
+   fec = FiniteElementCollection::New(glob_fes->FEColl()->Name());
+   // create a local ParFiniteElementSpace from the global one:
+   fes = pfes = new ParFiniteElementSpace(pmesh, glob_fes, partitioning, fec);
    SetSize(pfes->GetVSize());
 
    if (partitioning)
    {
+      // Assumption: the map "local element id" -> "global element id" is
+      // increasing, i.e. the local numbering preserves the element order from
+      // the global numbering.
       Array<int> gvdofs, lvdofs;
       Vector lnodes;
       int element_counter = 0;
-      Mesh & mesh(*gf->FESpace()->GetMesh());
-      int MyRank;
-      MPI_Comm_rank(pfes->GetComm(), &MyRank);
-      for (int i = 0; i < mesh.GetNE(); i++)
+      const int MyRank = pfes->GetMyRank();
+      const int glob_ne = glob_fes->GetNE();
+      for (int i = 0; i < glob_ne; i++)
+      {
          if (partitioning[i] == MyRank)
          {
             pfes->GetElementVDofs(element_counter, lvdofs);
-            gf->FESpace()->GetElementVDofs(i, gvdofs);
+            glob_fes->GetElementVDofs(i, gvdofs);
             gf->GetSubVector(gvdofs, lnodes);
             SetSubVector(lvdofs, lnodes);
             element_counter++;
          }
+      }
    }
 }
 

fem/pgridfunc.hpp

@@ -57,13 +57,15 @@ public:
        processor. The ParGridFunction does not assume ownership of the data. */
    ParGridFunction(ParFiniteElementSpace *pf, GridFunction *gf);
 
-   /** Creates grid function on (all) dofs from a given vector on the true dofs,
-       i.e. P tv. */
+   /** @brief Creates grid function on (all) dofs from a given vector on the
+       true dofs, i.e. P tv. */
    ParGridFunction(ParFiniteElementSpace *pf, HypreParVector *tv);
 
-   /** Construct a ParGridFunction from the given serial GridFunction.
-       If partitioning == NULL (default), the data from 'gf' is NOT copied. */
-   ParGridFunction(ParMesh *pmesh, GridFunction *gf, int * partitioning = NULL);
+   /** @brief Construct a local ParGridFunction from the given *global*
+       GridFunction. If @a partitioning is NULL (default), the data from @a gf
+       is NOT copied. */
+   ParGridFunction(ParMesh *pmesh, const GridFunction *gf,
+                   const int *partitioning = NULL);
 
    /// Assign constant values to the ParGridFunction data.
    ParGridFunction &operator=(double value)

general/array.hpp

@@ -79,9 +79,32 @@ public:
    inline Array(T *_data, int asize, int ainc = 0)
    { data = _data; size = asize; allocsize = -asize; inc = ainc; }
 
+   /// Copy constructor: deep copy
+   Array(const Array<T> &src)
+      : BaseArray(src.size, 0, sizeof(T))
+   { std::memcpy(data, src.data, size*sizeof(T)); }
+
+   /// Copy constructor (deep copy) from an Array of convertable type
+   template <typename CT>
+   Array(const Array<CT> &src)
+      : BaseArray(src.Size(), 0, sizeof(T))
+   { for (int i = 0; i < size; i++) { (*this)[i] = T(src[i]); } }
+
    /// Destructor
    inline ~Array() { }
 
+   /// Assignment operator: deep copy
+   Array<T> &operator=(const Array<T> &src) { src.Copy(*this); return *this; }
+
+   /// Assignment operator (deep copy) from an Array of convertable type
+   template <typename CT>
+   Array<T> &operator=(const Array<CT> &src)
+   {
+      SetSize(src.Size());
+      for (int i = 0; i < size; i++) { (*this)[i] = T(src[i]); }
+      return *this;
+   }
+
    /// Return the data as 'T *'
    inline operator T *() { return (T *)data; }
 
@@ -161,7 +184,7 @@ public:
    inline void Copy(Array &copy) const
    {
       copy.SetSize(Size());
-      memcpy(copy.GetData(), data, Size()*sizeof(T));
+      std::memcpy(copy.GetData(), data, Size()*sizeof(T));
    }
 
    /// Make this Array a reference to a pointer
@@ -238,12 +261,6 @@ public:
    inline T* end() const { return (T*) data + size; }
 
    long MemoryUsage() const { return Capacity() * sizeof(T); }
-
-private:
-   /// Array copy is not supported
-   Array<T> &operator=(Array<T> &);
-   /// Array copy is not supported
-   Array(const Array<T> &);
 };
 
 template <class T>
@@ -276,15 +293,15 @@ private:
    friend void Swap<T>(Array2D<T> &, Array2D<T> &);
 
    Array<T> array1d;
-   int M, N; // number of rows and columns
+   int N; // number of columns
 
 public:
-   Array2D() { M = N = 0; }
-   Array2D(int m, int n) : array1d(m*n) { M = m; N = n; }
+   Array2D() { N = 0; }
+   Array2D(int m, int n) : array1d(m*n) { N = n; }
 
-   void SetSize(int m, int n) { array1d.SetSize(m*n); M = m; N = n; }
+   void SetSize(int m, int n) { array1d.SetSize(m*n); N = n; }
 
-   int NumRows() const { return M; }
+   int NumRows() const { return array1d.Size()/N; }
    int NumCols() const { return N; }
 
    inline const T &operator()(int i, int j) const;
@@ -326,7 +343,7 @@ public:
    */
    void Load(std::istream &in, int fmt = 0)
    {
-      if (fmt == 0) { in >> M >> N; array1d.SetSize(M*N); }
+      if (fmt == 0) { int M; in >> M >> N; array1d.SetSize(M*N); }
       array1d.Load(in, 1);
    }
 
@@ -339,14 +356,17 @@ public:
    { SetSize(new_size0,new_size1); Load(in, 1); }
 
    void Copy(Array2D &copy) const
-   { copy.M = M; copy.N = N; array1d.Copy(copy.array1d); }
+   { copy.N = N; array1d.Copy(copy.array1d); }
 
    inline void operator=(const T &a)
    { array1d = a; }
 
    /// Make this Array a reference to 'master'
    inline void MakeRef(const Array2D &master)
-   { M = master.M; N = master.N; array1d.MakeRef(master.array1d); }
+   { N = master.N; array1d.MakeRef(master.array1d); }
+
+   /// Delete all dynamically allocated memory, reseting all dimentions to zero.
+   inline void DeleteAll() { N = 0; array1d.DeleteAll(); }
 
    /// Prints array to stream with width elements per row
    void Print(std::ostream &out = mfem::out, int width = 4);

general/communication.cpp

@@ -57,6 +57,9 @@ void GroupTopology::ProcToLProc()
 
    map<int, int> proc_lproc;
 
+   // The local processor ids are assigned following the group order and within
+   // a group following their ordering in the group. In other words, the ids are
+   // assigned based on their order in the J array of group_lproc.
    int lproc_counter = 0;
    for (int i = 0; i < group_lproc.Size_of_connections(); i++)
    {
@@ -122,88 +125,142 @@ void GroupTopology::Create(ListOfIntegerSets &groups, int mpitag)
 
    ProcToLProc();
 
-   // build group_mgroup
+   // Build 'group_mgroup':
+
+   // Use aggregated neighbor communication: at most one send to and/or one
+   // receive from each neighbor.
+
    group_mgroup.SetSize(NGroups());
-   group_mgroup[0] = 0; // the local group
-
-   int send_counter = 0;
-   int recv_counter = 0;
-   for (int i = 1; i < NGroups(); i++)
-      if (groupmaster_lproc[i] != 0) // we are not the master
-      {
-         recv_counter++;
-      }
-      else
-      {
-         send_counter += group_lproc.RowSize(i)-1;
-      }
-
-   MPI_Request *requests = new MPI_Request[send_counter];
-   MPI_Status  *statuses = new MPI_Status[send_counter];
-
-   int max_recv_size = 0;
-   send_counter = 0;
-   for (int i = 1; i < NGroups(); i++)
+   MFEM_DEBUG_DO(group_mgroup = -1);
+   for (int g = 0; g < NGroups(); g++)
    {
-      if (groupmaster_lproc[i] == 0) // we are the master
-      {
-         group_mgroup[i] = i;
+      if (IAmMaster(g)) { group_mgroup[g] = g; }
+   }
 
-         for (int j = group_lproc.GetI()[i];
-              j < group_lproc.GetI()[i+1]; j++)
+   // The Table 'lproc_cgroup': for each lproc, list the groups that are owned
+   // by this rank or by that lproc.
+   Table lproc_cgroup;
+   {
+      Array<Connection> lproc_cgroup_list;
+      for (int g = 1; g < NGroups(); g++)
+      {
+         if (IAmMaster(g))
          {
-            if (group_lproc.GetJ()[j] != 0)
+            const int gs = GetGroupSize(g);
+            const int *lprocs = GetGroup(g);
+            for (int i = 0; i < gs; i++)
             {
-               MPI_Isend(group_mgroupandproc.GetRow(i),
-                         group_mgroupandproc.RowSize(i),
-                         MPI_INT,
-                         lproc_proc[group_lproc.GetJ()[j]],
-                         mpitag,
-                         MyComm,
-                         &requests[send_counter]);
-               send_counter++;
+               if (lprocs[i])
+               {
+                  lproc_cgroup_list.Append(Connection(lprocs[i],g));
+               }
             }
          }
-      }
-      else // we are not the master
-         if (max_recv_size < group_lproc.RowSize(i))
+         else
          {
-            max_recv_size = group_lproc.RowSize(i);
+            lproc_cgroup_list.Append(Connection(GetGroupMaster(g),g));
          }
+      }
+      lproc_cgroup_list.Sort();
+      lproc_cgroup_list.Unique();
+      lproc_cgroup.MakeFromList(GetNumNeighbors(), lproc_cgroup_list);
    }
-   max_recv_size++;
 
-   IntegerSet group;
-   if (recv_counter > 0)
+   // Determine size of the send-receive buffer. For each neighbor the buffer
+   // contains: <send-part><receive-part> with each part consisting of a list of
+   // groups. Each group, g, has group_lproc.RowSize(g)+2 integers: the first
+   // entry is group_lproc.RowSize(g) - the number of processors in the group,
+   // followed by the group-id in the master processor, followed by the ranks of
+   // the processors in the group.
+   Table buffer;
+   buffer.MakeI(2*lproc_cgroup.Size()-2); // excluding the "local" lproc, 0
+   for (int nbr = 1; nbr < lproc_cgroup.Size(); nbr++)
    {
-      int count;
-      MPI_Status status;
-      int *recv_buf = new int[max_recv_size];
-      for ( ; recv_counter > 0; recv_counter--)
+      const int send_row = 2*(nbr-1);
+      const int recv_row = send_row+1;
+      const int ng = lproc_cgroup.RowSize(nbr);
+      const int *g = lproc_cgroup.GetRow(nbr);
+      for (int j = 0; j < ng; j++)
       {
-         MPI_Recv(recv_buf, max_recv_size, MPI_INT,
-                  MPI_ANY_SOURCE, mpitag, MyComm, &status);
-
-         MPI_Get_count(&status, MPI_INT, &count);
-
-         group.Recreate(count-1, recv_buf+1);
-         int g = groups.Lookup(group);
-         group_mgroup[g] = recv_buf[0];
-
-         if (lproc_proc[groupmaster_lproc[g]] != status.MPI_SOURCE)
+         const int gs = group_lproc.RowSize(g[j]);
+         if (IAmMaster(g[j]))
          {
-            mfem::err << "\n\n\nGroupTopology::GroupTopology: "
-                      << MyRank() << ": ERROR\n\n\n" << endl;
-            mfem_error();
+            buffer.AddColumnsInRow(send_row, gs+2);
+         }
+         else
+         {
+            MFEM_ASSERT(GetGroupMaster(g[j]) == nbr, "internal error");
+            buffer.AddColumnsInRow(recv_row, gs+2);
+         }
+      }
+   }
+   buffer.MakeJ();
+   for (int nbr = 1; nbr < lproc_cgroup.Size(); nbr++)
+   {
+      const int send_row = 2*(nbr-1);
+      const int recv_row = send_row+1;
+      const int ng = lproc_cgroup.RowSize(nbr);
+      const int *g = lproc_cgroup.GetRow(nbr);
+      for (int j = 0; j < ng; j++)
+      {
+         const int gs = group_lproc.RowSize(g[j]);
+         if (IAmMaster(g[j]))
+         {
+            buffer.AddConnection(send_row, gs);
+            buffer.AddConnections(
+               send_row, group_mgroupandproc.GetRow(g[j]), gs+1);
+         }
+         else
+         {
+            buffer.AddColumnsInRow(recv_row, gs+2);
+         }
+      }
+   }
+   buffer.ShiftUpI();
+   Array<MPI_Request> send_requests(lproc_cgroup.Size()-1);
+   Array<MPI_Request> recv_requests(lproc_cgroup.Size()-1);
+   send_requests = MPI_REQUEST_NULL;
+   recv_requests = MPI_REQUEST_NULL;
+   for (int nbr = 1; nbr < lproc_cgroup.Size(); nbr++)
+   {
+      const int send_row = 2*(nbr-1);
+      const int recv_row = send_row+1;
+      const int send_size = buffer.RowSize(send_row);
+      const int recv_size = buffer.RowSize(recv_row);
+      if (send_size > 0)
+      {
+         MPI_Isend(buffer.GetRow(send_row), send_size, MPI_INT, lproc_proc[nbr],
+                   mpitag, MyComm, &send_requests[nbr-1]);
+      }
+      if (recv_size > 0)
+      {
+         MPI_Irecv(buffer.GetRow(recv_row), recv_size, MPI_INT, lproc_proc[nbr],
+                   mpitag, MyComm, &recv_requests[nbr-1]);
+      }
+   }
+   {
+      int idx;
+      IntegerSet group;
+      while (MPI_Waitany(recv_requests.Size(), recv_requests.GetData(), &idx,
+                         MPI_STATUS_IGNORE),
+             idx != MPI_UNDEFINED)
+      {
+         const int recv_size = buffer.RowSize(2*idx+1);
+         const int *recv_buf = buffer.GetRow(2*idx+1);
+         for (int s = 0;  s < recv_size; s += recv_buf[s]+2)
+         {
+            group.Recreate(recv_buf[s], recv_buf+s+2);
+            const int g = groups.Lookup(group);
+            MFEM_ASSERT(group_mgroup[g] == -1, "communication error");
+            group_mgroup[g] = recv_buf[s+1];
          }
       }
-      delete [] recv_buf;
    }
 
-   MPI_Waitall(send_counter, requests, statuses);
+   MPI_Waitall(send_requests.Size(), send_requests.GetData(),
+               MPI_STATUSES_IGNORE);
 
-   delete [] statuses;
-   delete [] requests;
+   // debug barrier: MPI_Barrier(MyComm);
 }
 
 void GroupTopology::Save(ostream &out) const
@@ -283,7 +340,7 @@ GroupCommunicator::GroupCommunicator(GroupTopology &gt, Mode m)
    buf_offsets = NULL;
 }
 
-void GroupCommunicator::Create(Array<int> &ldof_group)
+void GroupCommunicator::Create(const Array<int> &ldof_group)
 {
    group_ldof.MakeI(gtopo.NGroups());
    for (int i = 0; i < ldof_group.Size(); i++)
@@ -418,7 +475,7 @@ void GroupCommunicator::Finalize()
    }
 }
 
-void GroupCommunicator::SetLTDofTable(Array<int> &ldof_ltdof)
+void GroupCommunicator::SetLTDofTable(const Array<int> &ldof_ltdof)
 {
    if (group_ltdof.Size() == group_ldof.Size()) { return; }
 

general/communication.hpp

@@ -78,9 +78,9 @@ public:
    GroupTopology(const GroupTopology &gt);
    void SetComm(MPI_Comm comm) { MyComm = comm; }
 
-   MPI_Comm GetComm() { return MyComm; }
-   int MyRank() { int r; MPI_Comm_rank(MyComm, &r); return r; }
-   int NRanks() { int s; MPI_Comm_size(MyComm, &s); return s; }
+   MPI_Comm GetComm() const { return MyComm; }
+   int MyRank() const { int r; MPI_Comm_rank(MyComm, &r); return r; }
+   int NRanks() const { int s; MPI_Comm_size(MyComm, &s); return s; }
 
    void Create(ListOfIntegerSets &groups, int mpitag);
 
@@ -152,7 +152,7 @@ public:
 
    /** @brief Initialize the communicator from a local-dof to group map.
        Finalize() is called internally. */
-   void Create(Array<int> &ldof_group);
+   void Create(const Array<int> &ldof_group);
 
    /** @brief Fill-in the returned Table reference to initialize the
        GroupCommunicator then call Finalize(). */
@@ -164,7 +164,7 @@ public:
    /// Initialize the internal group_ltdof Table.
    /** This method must be called before performing operations that use local
        data layout 2, see CopyGroupToBuffer() for layout descriptions. */
-   void SetLTDofTable(Array<int> &ldof_ltdof);
+   void SetLTDofTable(const Array<int> &ldof_ltdof);
 
    /// Get a reference to the associated GroupTopology object
    GroupTopology &GetGroupTopology() { return gtopo; }

general/error.hpp

@@ -68,7 +68,7 @@ void mfem_warning(const char *msg = NULL);
    if (!(x))                                            \
    {                                                    \
       _MFEM_MESSAGE("Verification failed: ("            \
-                    << #x << ") is false: " << msg, 0); \
+                    << #x << ") is false:\n --> " << msg, 0); \
    }
 
 // Use this if the only place your variable is used is in ASSERTs
@@ -86,7 +86,7 @@ void mfem_warning(const char *msg = NULL);
    if (!(x))                                            \
    {                                                    \
       _MFEM_MESSAGE("Assertion failed: ("               \
-                    << #x << ") is false: " << msg, 0); \
+                    << #x << ") is false:\n --> " << msg, 0); \
    }
 
 // A macro that exposes its argument in debug mode only.

general/globals.cpp

@@ -13,6 +13,8 @@
 #include "../config/config.hpp"
 #include "globals.hpp"
 #include <iostream>
+#include <sstream>
+#include <iomanip>
 
 namespace mfem
 {
@@ -20,6 +22,16 @@ namespace mfem
 OutStream out(std::cout);
 OutStream err(std::cerr);
 
+
+std::string MakeParFilename(const std::string &prefix, const int myid,
+                            const std::string suffix, const int width)
+{
+   std::stringstream fname;
+   fname << prefix << std::setw(width) << std::setfill('0') << myid << suffix;
+   return fname.str();
+}
+
+
 #ifdef MFEM_USE_MPI
 
 MPI_Comm MFEM_COMM_WORLD = MPI_COMM_WORLD;

general/globals.hpp

@@ -69,6 +69,20 @@ extern OutStream out;
 extern OutStream err;
 
 
+/** @brief Construct a string of the form "<prefix><myid><suffix>" where the
+    integer @a myid is padded with leading zeros to be at least @a width digits
+    long. */
+/** This is a convenience function, e.g. to redirect mfem::out to individual
+    files for each rank, one can use:
+    \code
+       std::ofstream out_file(MakeParFilename("app_out.", myid).c_str());
+       mfem::out.SetStream(out_file);
+    \endcode
+*/
+std::string MakeParFilename(const std::string &prefix, const int myid,
+                            const std::string suffix = "", const int width = 6);
+
+
 #ifdef MFEM_USE_MPI
 
 /** @name MFEM "global" communicator functions.

general/table.hpp

@@ -64,7 +64,7 @@ public:
    /// Create a table with an upper limit for the number of connections.
    explicit Table (int dim, int connections_per_row = 3);
 
-   /** Create a table from a list of connections, see MakeFromList. */
+   /** Create a table from a list of connections, see MakeFromList(). */
    Table(int nrows, Array<Connection> &list) : size(-1), I(NULL), J(NULL)
    { MakeFromList(nrows, list); }
 

general/text.hpp

@@ -13,10 +13,11 @@
 #define MFEM_TEXT
 
 #include <istream>
+#include <iomanip>
 #include <sstream>
 #include <string>
 #include <limits>
-#include <iomanip>
+#include <algorithm>
 
 namespace mfem
 {

linalg/ode.cpp

@@ -485,4 +485,109 @@ void SDIRK33Solver::Step(Vector &x, double &t, double &dt)
    t += dt;
 }
 
+void
+SIASolver::Init(Operator &P, TimeDependentOperator & F)
+{
+   P_ = &P; F_ = &F;
+
+   dp_.SetSize(F_->Height());
+   dq_.SetSize(P_->Height());
+}
+
+void
+SIA1Solver::Step(Vector &q, Vector &p, double &t, double &dt)
+{
+   F_->SetTime(t);
+   F_->Mult(q,dp_);
+   p.Add(dt,dp_);
+
+   P_->Mult(p,dq_);
+   q.Add(dt,dq_);
+
+   t += dt;
+}
+
+void
+SIA2Solver::Step(Vector &q, Vector &p, double &t, double &dt)
+{
+   P_->Mult(p,dq_);
+   q.Add(0.5*dt,dq_);
+
+   F_->SetTime(t+0.5*dt);
+   F_->Mult(q,dp_);
+   p.Add(dt,dp_);
+
+   P_->Mult(p,dq_);
+   q.Add(0.5*dt,dq_);
+
+   t += dt;
+}
+
+SIAVSolver::SIAVSolver(int order)
+   : order_(order)
+{
+   a_.SetSize(order);
+   b_.SetSize(order);
+
+   switch (order_)
+   {
+      case 1:
+         a_[0] = 1.0;
+         b_[0] = 1.0;
+         break;
+      case 2:
+         a_[0] = 0.5;
+         a_[1] = 0.5;
+         b_[0] = 0.0;
+         b_[1] = 1.0;
+         break;
+      case 3:
+         a_[0] =  2.0/3.0;
+         a_[1] = -2.0/3.0;
+         a_[2] =  1.0;
+         b_[0] =  7.0/24.0;
+         b_[1] =  0.75;
+         b_[2] = -1.0/24.0;
+         break;
+      case 4:
+         a_[0] = (2.0+pow(2.0,1.0/3.0)+pow(2.0,-1.0/3.0))/6.0;
+         a_[1] = (1.0-pow(2.0,1.0/3.0)-pow(2.0,-1.0/3.0))/6.0;
+         a_[2] = a_[1];
+         a_[3] = a_[0];
+         b_[0] = 0.0;
+         b_[1] = 1.0/(2.0-pow(2.0,1.0/3.0));
+         b_[2] = 1.0/(1.0-pow(2.0,2.0/3.0));
+         b_[3] = b_[1];
+         break;
+      default:
+         MFEM_ASSERT(false, "Unsupported order in SIAVSolver");
+   };
+}
+
+void
+SIAVSolver::Step(Vector &q, Vector &p, double &t, double &dt)
+{
+   for (int i=0; i<order_; i++)
+   {
+      if ( b_[i] != 0.0 )
+      {
+         F_->SetTime(t);
+         if ( F_->isExplicit() )
+         {
+            F_->Mult(q, dp_);
+         }
+         else
+         {
+            F_->ImplicitSolve(b_[i] * dt, q, dp_);
+         }
+         p.Add(b_[i] * dt, dp_);
+      }
+
+      P_->Mult(p, dq_);
+      q.Add(a_[i] * dt, dq_);
+
+      t += a_[i] * dt;
+   }
+}
+
 }

linalg/ode.hpp

@@ -282,6 +282,75 @@ public:
    virtual void Step(Vector &x, double &t, double &dt);
 };
 
+/// The SIASolver class is based on the Symplectic Integration Algorithm
+/// described in "A Symplectic Integration Algorithm for Separable Hamiltonian
+/// Functions" by J. Candy and W. Rozmus, Journal of Computational Physics,
+/// Vol. 92, pages 230-256 (1991).
+
+/** The Symplectic Integration Algorithm (SIA) is designed for systems of first
+    order ODEs derived from a Hamiltonian.
+       H(q,p,t) = T(p) + V(q,t)
+    Which leads to the equations:
+       dq/dt = dT/dp
+       dp/dt = -dV/dq
+    In the integrator the operators P and F are defined to be:
+       P = dT/dp
+       F = -dV/dq
+ */
+class SIASolver
+{
+public:
+   SIASolver() : F_(NULL), P_(NULL) {}
+
+   virtual void Init(Operator &P, TimeDependentOperator & F);
+
+   virtual void Step(Vector &q, Vector &p, double &t, double &dt) = 0;
+
+   virtual void Run(Vector &q, Vector &p, double &t, double &dt, double tf)
+   {
+      while (t < tf) { Step(q, p, t, dt); }
+   }
+
+   virtual ~SIASolver() {}
+
+protected:
+   TimeDependentOperator * F_; // p_{i+1} = p_{i} + dt F(q_{i})
+   Operator              * P_; // q_{i+1} = q_{i} + dt P(p_{i+1})
+
+   mutable Vector dp_;
+   mutable Vector dq_;
+};
+
+// First order Symplectic Integration Algorithm
+class SIA1Solver : public SIASolver
+{
+public:
+   SIA1Solver() {}
+   void Step(Vector &q, Vector &p, double &t, double &dt);
+};
+
+// Second order Symplectic Integration Algorithm
+class SIA2Solver : public SIASolver
+{
+public:
+   SIA2Solver() {}
+   void Step(Vector &q, Vector &p, double &t, double &dt);
+};
+
+// Variable order Symplectic Integration Algorithm (orders 1-4)
+class SIAVSolver : public SIASolver
+{
+public:
+   SIAVSolver(int order);
+   void Step(Vector &q, Vector &p, double &t, double &dt);
+
+private:
+   int order_;
+
+   Array<double> a_;
+   Array<double> b_;
+};
+
 }
 
 #endif

makefile

@@ -93,7 +93,7 @@ EXAMPLE_SUBDIRS = sundials petsc
 EXAMPLE_DIRS := examples $(addprefix examples/,$(EXAMPLE_SUBDIRS))
 EXAMPLE_TEST_DIRS := examples
 
-MINIAPP_SUBDIRS = common electromagnetics meshing performance tools
+MINIAPP_SUBDIRS = common electromagnetics meshing performance tools nurbs
 MINIAPP_DIRS := $(addprefix miniapps/,$(MINIAPP_SUBDIRS))
 MINIAPP_TEST_DIRS := $(filter-out %/common,$(MINIAPP_DIRS))
 MINIAPP_USE_COMMON := $(addprefix miniapps/,electromagnetics tools)

mesh/mesh.cpp

@@ -622,34 +622,34 @@ void Mesh::GetLocalQuadToHexTransformation(
 }
 
 void Mesh::GetLocalFaceTransformation(
-   int face_type, int elem_type, IsoparametricTransformation &Transf, int inf)
+   int face_type, int elem_type, IsoparametricTransformation &Transf, int info)
 {
    switch (face_type)
    {
       case Element::POINT:
-         GetLocalPtToSegTransformation(Transf, inf);
+         GetLocalPtToSegTransformation(Transf, info);
          break;
 
       case Element::SEGMENT:
          if (elem_type == Element::TRIANGLE)
          {
-            GetLocalSegToTriTransformation(Transf, inf);
+            GetLocalSegToTriTransformation(Transf, info);
          }
          else
          {
             MFEM_ASSERT(elem_type == Element::QUADRILATERAL, "");
-            GetLocalSegToQuadTransformation(Transf, inf);
+            GetLocalSegToQuadTransformation(Transf, info);
          }
          break;
 
       case Element::TRIANGLE:
          MFEM_ASSERT(elem_type == Element::TETRAHEDRON, "");
-         GetLocalTriToTetTransformation(Transf, inf);
+         GetLocalTriToTetTransformation(Transf, info);
          break;
 
       case Element::QUADRILATERAL:
          MFEM_ASSERT(elem_type == Element::HEXAHEDRON, "");
-         GetLocalQuadToHexTransformation(Transf, inf);
+         GetLocalQuadToHexTransformation(Transf, info);
          break;
    }
 }
@@ -2342,7 +2342,6 @@ Mesh::Mesh(const Mesh &mesh, bool copy_nodes)
    }
 
    // Copy the vertices
-   MFEM_ASSERT(mesh.vertices.Size() == NumOfVertices, "internal MFEM error!");
    mesh.vertices.Copy(vertices);
 
    // Duplicate the boundary
@@ -2389,12 +2388,23 @@ Mesh::Mesh(const Mesh &mesh, bool copy_nodes)
    mesh.attributes.Copy(attributes);
    mesh.bdr_attributes.Copy(bdr_attributes);
 
-   // No support for NURBS meshes, yet. Need deep copy for NURBSExtension.
-   MFEM_VERIFY(mesh.NURBSext == NULL,
-               "copying NURBS meshes is not implemented");
-   NURBSext = NULL;
+   // Deep copy the NURBSExtension.
+#ifdef MFEM_USE_MPI
+   ParNURBSExtension *pNURBSext =
+      dynamic_cast<ParNURBSExtension *>(mesh.NURBSext);
+   if (pNURBSext)
+   {
+      NURBSext = new ParNURBSExtension(*pNURBSext);
+   }
+   else
+#endif
+   {
+      NURBSext = mesh.NURBSext ? new NURBSExtension(*mesh.NURBSext) : NULL;
+   }
 
    // Deep copy the NCMesh.
+   // TODO: ParNCMesh; ParMesh has a separate 'pncmesh' pointer, and 'ncmesh'
+   //       is initialized from it. Need ParNCMesh copy constructor.
    ncmesh = mesh.ncmesh ? new NCMesh(*mesh.ncmesh) : NULL;
 
    // Duplicate the Nodes, including the FiniteElementCollection and the
@@ -2406,8 +2416,7 @@ Mesh::Mesh(const Mesh &mesh, bool copy_nodes)
       FiniteElementCollection *fec_copy =
          FiniteElementCollection::New(fec->Name());
       FiniteElementSpace *fes_copy =
-         new FiniteElementSpace(this, fec_copy, fes->GetVDim(),
-                                fes->GetOrdering());
+         new FiniteElementSpace(*fes, this, fec_copy);
       Nodes = new GridFunction(fes_copy);
       Nodes->MakeOwner(fec_copy);
       *Nodes = *mesh.Nodes;
@@ -3082,6 +3091,9 @@ void Mesh::KnotInsert(Array<KnotVector *> &kv)
 
    NURBSext->KnotInsert(kv);
 
+   last_operation = Mesh::NONE; // FiniteElementSpace::Update is not supported
+   sequence++;
+
    UpdateNURBS();
 }
 
@@ -3092,13 +3104,13 @@ void Mesh::NURBSUniformRefinement()
 
    NURBSext->UniformRefinement();
 
-   last_operation = Mesh::REFINE;
+   last_operation = Mesh::NONE; // FiniteElementSpace::Update is not supported
    sequence++;
 
    UpdateNURBS();
 }
 
-void Mesh::DegreeElevate(int t)
+void Mesh::DegreeElevate(int rel_degree, int degree)
 {
    if (NURBSext == NULL)
    {
@@ -3107,15 +3119,10 @@ void Mesh::DegreeElevate(int t)
 
    NURBSext->ConvertToPatches(*Nodes);
 
-   NURBSext->DegreeElevate(t);
+   NURBSext->DegreeElevate(rel_degree, degree);
 
-   NURBSFECollection *nurbs_fec =
-      dynamic_cast<NURBSFECollection *>(Nodes->OwnFEC());
-   if (!nurbs_fec)
-   {
-      mfem_error("Mesh::DegreeElevate");
-   }
-   nurbs_fec->UpdateOrder(nurbs_fec->GetOrder() + t);
+   last_operation = Mesh::NONE; // FiniteElementSpace::Update is not supported
+   sequence++;
 
    UpdateNURBS();
 }

mesh/mesh.hpp

@@ -67,7 +67,9 @@ protected:
 
    Array<Element *> elements;
    // Vertices are only at the corners of elements, where you would expect them
-   // in the lowest-order mesh.
+   // in the lowest-order mesh. In some cases, e.g. in a Mesh that defines the
+   // patch topology for a NURBS mesh (see LoadPatchTopo()) the vertices may be
+   // empty while NumOfVertices is positive.
    Array<Vertex> vertices;
    Array<Element *> boundary;
    Array<Element *> faces;
@@ -307,10 +309,6 @@ protected:
    /// Used in GetFaceElementTransformations (...)
    void GetLocalQuadToHexTransformation (IsoparametricTransformation &loc,
                                          int i);
-   /// Used in GetFaceElementTransformations (...)
-   void GetLocalFaceTransformation(int face_type, int elem_type,
-                                   IsoparametricTransformation &Transf,
-                                   int inf);
    /** Used in GetFaceElementTransformations to account for the fact that a
        slave face occupies only a portion of its master face. */
    void ApplyLocalSlaveTransformation(IsoparametricTransformation &transf,
@@ -404,7 +402,7 @@ protected:
    void InitFromNCMesh(const NCMesh &ncmesh);
 
    /// Create from a nonconforming mesh.
-   Mesh(const NCMesh &ncmesh);
+   explicit Mesh(const NCMesh &ncmesh);
 
    /// Swaps internal data with another mesh. By default, non-geometry members
    /// like 'ncmesh' and 'NURBSExt' are only swapped when 'non_geometry' is set.
@@ -562,14 +560,14 @@ public:
    /** Creates mesh by reading a file in MFEM, netgen, or VTK format. If
        generate_edges = 0 (default) edges are not generated, if 1 edges are
        generated. */
-   Mesh(const char *filename, int generate_edges = 0, int refine = 1,
-        bool fix_orientation = true);
+   explicit Mesh(const char *filename, int generate_edges = 0, int refine = 1,
+                 bool fix_orientation = true);
 
    /** Creates mesh by reading data stream in MFEM, netgen, or VTK format. If
        generate_edges = 0 (default) edges are not generated, if 1 edges are
        generated. */
-   Mesh(std::istream &input, int generate_edges = 0, int refine = 1,
-        bool fix_orientation = true);
+   explicit Mesh(std::istream &input, int generate_edges = 0, int refine = 1,
+                 bool fix_orientation = true);
 
    /// Create a disjoint mesh from the given mesh array
    Mesh(Mesh *mesh_array[], int num_pieces);
@@ -779,10 +777,19 @@ public:
    ElementTransformation * GetBdrElementTransformation(int i);
    void GetBdrElementTransformation(int i, IsoparametricTransformation *ElTr);
 
-   /** Returns the transformation defining the given face element.
-       The transformation is stored in a user-defined variable. */
+   /** @brief Returns the transformation defining the given face element in a
+       user-defined variable. */
    void GetFaceTransformation(int i, IsoparametricTransformation *FTr);
 
+   /** @brief A helper method that constructs a transformation from the
+       reference space of a face to the reference space of an element. */
+   /** The local index of the face as a face in the element and its orientation
+       are given by the input parameter @a info, as @a info = 64*loc_face_idx +
+       loc_face_orientation. */
+   void GetLocalFaceTransformation(int face_type, int elem_type,
+                                   IsoparametricTransformation &Transf,
+                                   int info);
+
    /// Returns the transformation defining the given face element
    ElementTransformation *GetFaceTransformation(int FaceNo);
 
@@ -992,7 +999,9 @@ public:
 
    ///@{ @name NURBS mesh refinement methods
    void KnotInsert(Array<KnotVector *> &kv);
-   void DegreeElevate(int t);
+   /* For each knot vector:
+         new_degree = max(old_degree, min(old_degree + rel_degree, degree)). */
+   void DegreeElevate(int rel_degree, int degree = 16);
    ///@}
 
    /** Make sure that a quad/hex mesh is considered to be non-conforming (i.e.,

mesh/ncmesh.hpp

@@ -86,7 +86,7 @@ public:
    /** Initialize with elements from 'mesh'. If an already nonconforming mesh
        is being loaded, 'vertex_parents' must point to a stream at the appropriate
        section of the mesh file which contains the vertex hierarchy. */
-   NCMesh(const Mesh *mesh, std::istream *vertex_parents = NULL);
+   explicit NCMesh(const Mesh *mesh, std::istream *vertex_parents = NULL);
 
    NCMesh(const NCMesh &other); // deep copy
 

mesh/nurbs.hpp

@@ -61,8 +61,8 @@ public:
 
    int findKnotSpan(double u) const;
 
-   void CalcShape (Vector &shape, int i, double xi);
-   void CalcDShape(Vector &grad,  int i, double xi);
+   void CalcShape (Vector &shape, int i, double xi) const;
+   void CalcDShape(Vector &grad,  int i, double xi) const;
 
    void Difference(const KnotVector &kv, Vector &diff) const;
    void UniformRefinement(Vector &newknots) const;
@@ -86,7 +86,7 @@ class NURBSPatch
 {
 protected:
    int     ni, nj, nk, Dim;
-   double *data;
+   double *data; // the layout of data is: (Dim x ni x nj x nk)
 
    Array<KnotVector *> kv;
 
@@ -104,14 +104,16 @@ protected:
    NURBSPatch(NURBSPatch *parent, int dir, int Order, int NCP);
 
 public:
+   NURBSPatch(const NURBSPatch &orig);
    NURBSPatch(std::istream &input);
-   NURBSPatch(KnotVector *kv0, KnotVector *kv1, int dim_);
-   NURBSPatch(KnotVector *kv0, KnotVector *kv1, KnotVector *kv2, int dim_);
-   NURBSPatch(Array<KnotVector *> &kv, int dim_);
+   NURBSPatch(const KnotVector *kv0, const KnotVector *kv1, int dim_);
+   NURBSPatch(const KnotVector *kv0, const KnotVector *kv1,
+              const KnotVector *kv2, int dim_);
+   NURBSPatch(Array<const KnotVector *> &kv, int dim_);
 
    ~NURBSPatch();
 
-   void Print(std::ostream &out);
+   void Print(std::ostream &out) const;
 
    void DegreeElevate(int dir, int t);
    void KnotInsert   (int dir, const KnotVector &knot);
@@ -121,6 +123,9 @@ public:
    void DegreeElevate(int t);
    void UniformRefinement();
 
+   // Return the number of components stored in the NURBSPatch
+   int GetNC() const { return Dim; }
+   int GetNKV() const { return kv.Size(); }
    KnotVector *GetKV(int i) { return kv[i]; }
 
    // Standard B-NET access functions
@@ -135,7 +140,7 @@ public:
    void FlipDirection(int dir);
    void SwapDirections(int dir1, int dir2);
    void Rotate3D(double normal[], double angle);
-   int MakeUniformDegree();
+   int MakeUniformDegree(int degree = -1);
    friend NURBSPatch *Interpolate(NURBSPatch &p1, NURBSPatch &p2);
    friend NURBSPatch *Revolve3D(NURBSPatch &patch, double n[], double ang,
                                 int times);
@@ -157,7 +162,8 @@ class NURBSExtension
    friend class NURBSPatchMap;
 
 protected:
-   int Order;
+   int mOrder; // see GetOrder() for description
+   Array<int> mOrders;
    int NumOfKnotVectors;
    // global entity counts
    int NumOfVertices, NumOfElements, NumOfBdrElements, NumOfDofs;
@@ -197,15 +203,21 @@ protected:
 
    Array<NURBSPatch *> patches;
 
-   inline int         KnotInd(int edge);
+   inline int         KnotInd(int edge) const;
    inline KnotVector *KnotVec(int edge);
-   inline KnotVector *KnotVec(int edge, int oedge, int *okv);
+   inline const KnotVector *KnotVec(int edge) const;
+   inline const KnotVector *KnotVec(int edge, int oedge, int *okv) const;
 
    void CheckPatches();
    void CheckBdrPatches();
 
    void GetPatchKnotVectors   (int p, Array<KnotVector *> &kv);
+   void GetPatchKnotVectors   (int p, Array<const KnotVector *> &kv) const;
    void GetBdrPatchKnotVectors(int p, Array<KnotVector *> &kv);
+   void GetBdrPatchKnotVectors(int p, Array<const KnotVector *> &kv) const;
+
+   void SetOrderFromOrders();
+   void SetOrdersFromKnotVectors();
 
    // also count the global NumOfVertices and the global NumOfDofs
    void GenerateOffsets();
@@ -215,12 +227,12 @@ protected:
    void CountBdrElements();
 
    // generate the mesh elements
-   void Get2DElementTopo(Array<Element *> &elements);
-   void Get3DElementTopo(Array<Element *> &elements);
+   void Get2DElementTopo(Array<Element *> &elements) const;
+   void Get3DElementTopo(Array<Element *> &elements) const;
 
    // generate the boundary mesh elements
-   void Get2DBdrElementTopo(Array<Element *> &boundary);
-   void Get3DBdrElementTopo(Array<Element *> &boundary);
+   void Get2DBdrElementTopo(Array<Element *> &boundary) const;
+   void Get3DBdrElementTopo(Array<Element *> &boundary) const;
 
 
    // FE space generation functions
@@ -242,14 +254,16 @@ protected:
    void Generate2DBdrElementDofTable();
    void Generate3DBdrElementDofTable();
 
-   // Patch <--> FE translation functions
-   void GetPatchNets  (const Vector &Nodes);
-   void Get2DPatchNets(const Vector &Nodes);
-   void Get3DPatchNets(const Vector &Nodes);
+   // FE --> Patch translation functions
+   void GetPatchNets  (const Vector &Nodes, int vdim);
+   void Get2DPatchNets(const Vector &Nodes, int vdim);
+   void Get3DPatchNets(const Vector &Nodes, int vdim);
 
-   void SetSolutionVector  (Vector &Nodes);
-   void Set2DSolutionVector(Vector &Nodes);
-   void Set3DSolutionVector(Vector &Nodes);
+   // Patch --> FE translation functions
+   // Side effects: delete the patches, update the weights from the patches
+   void SetSolutionVector  (Vector &Nodes, int vdim);
+   void Set2DSolutionVector(Vector &Nodes, int vdim);
+   void Set3DSolutionVector(Vector &Nodes, int vdim);
 
    // determine activeVert, NumOfActiveVertices from the activeElem array
    void GenerateActiveVertices();
@@ -263,11 +277,21 @@ protected:
    NURBSExtension() { }
 
 public:
+   /// Copy constructor: deep copy
+   NURBSExtension(const NURBSExtension &orig);
    /// Read-in a NURBSExtension
    NURBSExtension(std::istream &input);
-   /** Create a NURBSExtension with elevated order by repeating the endpoints
-       of the knot vectors and using uniform weights of 1. */
-   NURBSExtension(NURBSExtension *parent, int Order);
+   /** @brief Create a NURBSExtension with elevated order by repeating the
+       endpoints of the knot vectors and using uniform weights of 1. */
+   /** If a knot vector in @a parent already has order greater than or equal to
+       @a newOrder, it will be used unmodified. */
+   NURBSExtension(NURBSExtension *parent, int newOrder);
+   /** @brief Create a NURBSExtension with elevated knot vector orders (by
+       repeating the endpoints of the knot vectors and using uniform weights of
+       1) as given by the array @a newOrders. */
+   /** If a knot vector in @a parent already has order greater than or equal to
+       the corresponding entry in @a newOrder, it will be used unmodified. */
+   NURBSExtension(NURBSExtension *parent, const Array<int> &newOrders);
    /// Construct a NURBSExtension by merging a partitioned NURBS mesh
    NURBSExtension(Mesh *mesh_array[], int num_pieces);
 
@@ -279,33 +303,39 @@ public:
 
    // Print functions
    void Print(std::ostream &out) const;
-   void PrintCharacteristics(std::ostream &out);
+   void PrintCharacteristics(std::ostream &out) const;
 
    // Meta data functions
-   int Dimension() { return patchTopo->Dimension(); }
-   int GetNP()     { return patchTopo->GetNE(); }
-   int GetNBP()    { return patchTopo->GetNBE(); }
-   int GetOrder()  { return Order; }
-   int GetNKV()    { return NumOfKnotVectors; }
+   int Dimension() const { return patchTopo->Dimension(); }
+   int GetNP()     const { return patchTopo->GetNE(); }
+   int GetNBP()    const { return patchTopo->GetNBE(); }
 
-   int GetGNV()  { return NumOfVertices; }
-   int GetNV()   { return NumOfActiveVertices; }
-   int GetGNE()  { return NumOfElements; }
-   int GetNE()   { return NumOfActiveElems; }
-   int GetGNBE() { return NumOfBdrElements; }
-   int GetNBE()  { return NumOfActiveBdrElems; }
+   /// Read-only access to the orders of all knot vectors.
+   const Array<int> &GetOrders() const { return mOrders; }
+   /** @brief If all orders are identical, return that number. Otherwise, return
+       NURBSFECollection::VariableOrder. */
+   int GetOrder() const { return mOrder; }
 
-   int GetNTotalDof() { return NumOfDofs; }
-   int GetNDof()      { return NumOfActiveDofs; }
+   int GetNKV()  const { return NumOfKnotVectors; }
 
-   // Knotvector access function
+   int GetGNV()  const { return NumOfVertices; }
+   int GetNV()   const { return NumOfActiveVertices; }
+   int GetGNE()  const { return NumOfElements; }
+   int GetNE()   const { return NumOfActiveElems; }
+   int GetGNBE() const { return NumOfBdrElements; }
+   int GetNBE()  const { return NumOfActiveBdrElems; }
+
+   int GetNTotalDof() const { return NumOfDofs; }
+   int GetNDof()      const { return NumOfActiveDofs; }
+
+   // Knotvector read-only access function
    const KnotVector *GetKnotVector(int i) const { return knotVectors[i]; }
 
    // Mesh generation functions
-   void GetElementTopo   (Array<Element *> &elements);
-   void GetBdrElementTopo(Array<Element *> &boundary);
+   void GetElementTopo   (Array<Element *> &elements) const;
+   void GetBdrElementTopo(Array<Element *> &boundary) const;
 
-   bool HavePatches() { return (patches.Size() != 0); }
+   bool HavePatches() const { return (patches.Size() != 0); }
 
    Table *GetElementDofTable() { return el_dof; }
    Table *GetBdrElementDofTable() { return bel_dof; }
@@ -314,20 +344,26 @@ public:
    void GetElementLocalToGlobal(Array<int> &lelem_elem);
 
    // Load functions
-   void LoadFE(int i, const FiniteElement *FE);
-   void LoadBE(int i, const FiniteElement *BE);
+   void LoadFE(int i, const FiniteElement *FE) const;
+   void LoadBE(int i, const FiniteElement *BE) const;
 
    const Vector &GetWeights() const { return  weights; }
    Vector       &GetWeights()       { return  weights; }
 
-   // Translation functions: from FE coordinates into to IJK patch
+   // Translation functions: from FE coordinates to IJK patch
    // format and vice versa
    void ConvertToPatches(const Vector &Nodes);
    void SetKnotsFromPatches();
    void SetCoordsFromPatches(Vector &Nodes);
 
+   // Read a GridFunction written patch-by-patch, e.g. with PrintSolution().
+   void LoadSolution(std::istream &input, GridFunction &sol) const;
+   // Write a GridFunction patch-by-patch.
+   void PrintSolution(const GridFunction &sol, std::ostream &out) const;
+
    // Refinement methods
-   void DegreeElevate(int t);
+   // new_degree = max(old_degree, min(old_degree + rel_degree, degree))
+   void DegreeElevate(int rel_degree, int degree = 16);
    void UniformRefinement();
    void KnotInsert(Array<KnotVector *> &kv);
 };
@@ -343,20 +379,24 @@ private:
    Table *Get2DGlobalElementDofTable();
    Table *Get3DGlobalElementDofTable();
 
-   void SetActive(int *partitioning, const Array<bool> &active_bel);
-   void BuildGroups(int *partitioning, const Table &elem_dof);
+   void SetActive(const int *partitioning, const Array<bool> &active_bel);
+   void BuildGroups(const int *partitioning, const Table &elem_dof);
 
 public:
    GroupTopology gtopo;
 
    Array<int> ldof_group;
 
+   ParNURBSExtension(const ParNURBSExtension &orig);
+
    ParNURBSExtension(MPI_Comm comm, NURBSExtension *parent, int *partitioning,
                      const Array<bool> &active_bel);
 
-   // create a parallel version of 'parent' with partitioning as in
-   // 'par_parent'; the 'parent' object is destroyed
-   ParNURBSExtension(NURBSExtension *parent, ParNURBSExtension *par_parent);
+   // Create a parallel version of 'parent' with partitioning as in
+   // 'par_parent'; the 'parent' object is destroyed.
+   // The 'parent' can be either a local NURBSExtension or a global one.
+   ParNURBSExtension(NURBSExtension *parent,
+                     const ParNURBSExtension *par_parent);
 
    virtual ~ParNURBSExtension() { delete [] partitioning; }
 };
@@ -366,7 +406,7 @@ public:
 class NURBSPatchMap
 {
 private:
-   NURBSExtension *Ext;
+   const NURBSExtension *Ext;
 
    int I, J, K, pOffset, opatch;
    Array<int> verts, edges, faces, oedge, oface;
@@ -381,21 +421,21 @@ private:
                           const int N1, const int N2, const int Or);
 
    // also set verts, edges, faces, orientations etc
-   void GetPatchKnotVectors   (int p, KnotVector *kv[]);
-   void GetBdrPatchKnotVectors(int p, KnotVector *kv[], int *okv);
+   void GetPatchKnotVectors   (int p, const KnotVector *kv[]);
+   void GetBdrPatchKnotVectors(int p, const KnotVector *kv[], int *okv);
 
 public:
-   NURBSPatchMap(NURBSExtension *ext) { Ext = ext; }
+   NURBSPatchMap(const NURBSExtension *ext) { Ext = ext; }
 
    int nx() { return I + 1; }
    int ny() { return J + 1; }
    int nz() { return K + 1; }
 
-   void SetPatchVertexMap(int p, KnotVector *kv[]);
-   void SetPatchDofMap   (int p, KnotVector *kv[]);
+   void SetPatchVertexMap(int p, const KnotVector *kv[]);
+   void SetPatchDofMap   (int p, const KnotVector *kv[]);
 
-   void SetBdrPatchVertexMap(int p, KnotVector *kv[], int *okv);
-   void SetBdrPatchDofMap   (int p, KnotVector *kv[], int *okv);
+   void SetBdrPatchVertexMap(int p, const KnotVector *kv[], int *okv);
+   void SetBdrPatchDofMap   (int p, const KnotVector *kv[], int *okv);
 
    inline int operator()(const int i) const;
    inline int operator[](const int i) const { return (*this)(i); }
@@ -471,7 +511,7 @@ inline const double &NURBSPatch::operator()(int i, int j, int k, int l) const
 }
 
 
-inline int NURBSExtension::KnotInd(int edge)
+inline int NURBSExtension::KnotInd(int edge) const
 {
    int kv = edge_to_knot[edge];
    return (kv >= 0) ? kv : (-1-kv);
@@ -482,7 +522,13 @@ inline KnotVector *NURBSExtension::KnotVec(int edge)
    return knotVectors[KnotInd(edge)];
 }
 
-inline KnotVector *NURBSExtension::KnotVec(int edge, int oedge, int *okv)
+inline const KnotVector *NURBSExtension::KnotVec(int edge) const
+{
+   return knotVectors[KnotInd(edge)];
+}
+
+inline const KnotVector *NURBSExtension::KnotVec(int edge, int oedge, int *okv)
+const
 {
    int kv = edge_to_knot[edge];
    if (kv >= 0)
@@ -498,6 +544,7 @@ inline KnotVector *NURBSExtension::KnotVec(int edge, int oedge, int *okv)
 }
 
 
+// static method
 inline int NURBSPatchMap::Or2D(const int n1, const int n2,
                                const int N1, const int N2, const int Or)
 {
@@ -521,7 +568,7 @@ inline int NURBSPatchMap::Or2D(const int n1, const int n2,
 
 inline int NURBSPatchMap::operator()(const int i) const
 {
-   int i1 = i - 1;
+   const int i1 = i - 1;
    switch (F(i1, I))
    {
       case 0: return verts[0];
@@ -536,7 +583,7 @@ inline int NURBSPatchMap::operator()(const int i) const
 
 inline int NURBSPatchMap::operator()(const int i, const int j) const
 {
-   int i1 = i - 1, j1 = j - 1;
+   const int i1 = i - 1, j1 = j - 1;
    switch (3*F(j1, J) + F(i1, I))
    {
       case 0: return verts[0];
@@ -559,7 +606,7 @@ inline int NURBSPatchMap::operator()(const int i, const int j, const int k)
 const
 {
    // Needs testing
-   int i1 = i - 1, j1 = j - 1, k1 = k - 1;
+   const int i1 = i - 1, j1 = j - 1, k1 = k - 1;
    switch (3*(3*F(k1, K) + F(j1, J)) + F(i1, I))
    {
       case  0: return verts[0];

mesh/pmesh.cpp

@@ -59,8 +59,11 @@ ParMesh::ParMesh(const ParMesh &pmesh, bool copy_nodes)
    // Do not copy face-neighbor data (can be generated if needed)
    have_face_nbr_data = false;
 
+   // If pmesh has a ParNURBSExtension, it was copied by the Mesh copy ctor, so
+   // there is no need to do anything here.
+
    MFEM_VERIFY(pmesh.pncmesh == NULL,
-               "copying non-conforming meshes is not implemented");
+               "copy of parallel non-conforming meshes is not implemented");
    pncmesh = NULL;
 
    // Copy the Nodes as a ParGridFunction, including the FiniteElementCollection
@@ -72,8 +75,7 @@ ParMesh::ParMesh(const ParMesh &pmesh, bool copy_nodes)
       FiniteElementCollection *fec_copy =
          FiniteElementCollection::New(fec->Name());
       ParFiniteElementSpace *pfes_copy =
-         new ParFiniteElementSpace(this, fec_copy, fes->GetVDim(),
-                                   fes->GetOrdering());
+         new ParFiniteElementSpace(*fes, *this, fec_copy);
       Nodes = new ParGridFunction(pfes_copy);
       Nodes->MakeOwner(fec_copy);
       *Nodes = *pmesh.Nodes;
@@ -173,7 +175,7 @@ ParMesh::ParMesh(MPI_Comm comm, Mesh &mesh, int *partitioning_,
    NumOfElements = element_counter;
    vertices.SetSize(NumOfVertices);
 
-   // re-enumerate the local vertices to preserve the global ordering
+   // Re-enumerate the local vertices to preserve the global ordering.
    for (i = vert_counter = 0; i < vert_global_local.Size(); i++)
       if (vert_global_local[i] >= 0)
       {
@@ -188,7 +190,9 @@ ParMesh::ParMesh(MPI_Comm comm, Mesh &mesh, int *partitioning_,
                                                   mesh.GetVertex(i));
       }
 
-   // determine elements
+   // Determine elements, enumerating the local elements to preserve the global
+   // order. This is used, e.g. by the ParGridFunction ctor that takes a global
+   // GridFunction as input parameter.
    element_counter = 0;
    elements.SetSize(NumOfElements);
    for (i = 0; i < mesh.GetNE(); i++)
@@ -639,13 +643,30 @@ ParMesh::ParMesh(MPI_Comm comm, Mesh &mesh, int *partitioning_,
 
    if (mesh.NURBSext)
    {
+      MFEM_ASSERT(mesh.GetNodes() &&
+                  mesh.GetNodes()->FESpace()->GetNURBSext() == mesh.NURBSext,
+                  "invalid NURBS mesh");
       NURBSext = new ParNURBSExtension(comm, mesh.NURBSext, partitioning,
                                        activeBdrElem);
    }
 
    if (mesh.GetNodes()) // curved mesh
    {
-      Nodes = new ParGridFunction(this, mesh.GetNodes());
+      if (!NURBSext)
+      {
+         Nodes = new ParGridFunction(this, mesh.GetNodes());
+      }
+      else
+      {
+         const FiniteElementSpace *glob_fes = mesh.GetNodes()->FESpace();
+         FiniteElementCollection *nfec =
+            FiniteElementCollection::New(glob_fes->FEColl()->Name());
+         ParFiniteElementSpace *pfes =
+            new ParFiniteElementSpace(this, nfec, glob_fes->GetVDim(),
+                                      glob_fes->GetOrdering());
+         Nodes = new ParGridFunction(pfes);
+         Nodes->MakeOwner(nfec); // Nodes will own nfec and pfes
+      }
       own_nodes = 1;
 
       Array<int> gvdofs, lvdofs;

miniapps/CMakeLists.txt

@@ -18,3 +18,4 @@ add_subdirectory(electromagnetics)
 add_subdirectory(meshing)
 add_subdirectory(performance)
 add_subdirectory(tools)
+add_subdirectory(nurbs)

miniapps/common/makefile

@@ -18,6 +18,7 @@ CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
 # MFEM_INSTALL_DIR = ../../mfem
 # CONFIG_MK = $(MFEM_INSTALL_DIR)/share/mfem/config.mk
 
+MFEM_LIB_FILE = mfem_is_not_built
 ifneq (clean,$(MAKECMDGOALS))
    -include $(CONFIG_MK)
 endif
@@ -46,6 +47,9 @@ all: $(MINIOBJS)
 $(CONFIG_MK):
 	$(error The MFEM library is not configured)
 
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
 clean:
 	rm -f *.o *~
 	rm -rf *.dSYM *.TVD.*breakpoints

miniapps/electromagnetics/CMakeLists.txt

@@ -28,6 +28,12 @@ if (MFEM_USE_MPI)
    EXTRA_HEADERS joule_solver.hpp ${MFEM_MINIAPPS_COMMON_HEADERS}
    LIBRARIES mfem mfem_miniapps_common)
 
+ add_mfem_miniapp(maxwell
+   MAIN maxwell.cpp
+   EXTRA_SOURCES maxwell_solver.cpp
+   EXTRA_HEADERS maxwell_solver.hpp ${MFEM_MINIAPPS_COMMON_HEADERS}
+   LIBRARIES mfem mfem_miniapps_common)
+
   # Add the corresponding tests to the "test" target
   add_test(NAME tesla_np=4
     COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4
@@ -48,4 +54,11 @@ if (MFEM_USE_MPI)
    -no-vis -p rod -tf 3 -m ${CMAKE_CURRENT_SOURCE_DIR}/cylinder-hex.mesh
    ${MPIEXEC_POSTFLAGS})
 
+ add_test(NAME maxwell_np=4
+   COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4
+   ${MPIEXEC_PREFLAGS}
+   $<TARGET_FILE:maxwell>
+   -no-vis -abcs "-1" -dp "-0.3 0.0 0.0 0.3 0.0 0.0 0.1 1 .5 .5"
+   ${MPIEXEC_POSTFLAGS})
+
 endif()

miniapps/electromagnetics/electromagnetics.hpp

@@ -0,0 +1,33 @@
+// 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.
+
+#ifndef MFEM_ELECTROMAGNETICS_HPP
+#define MFEM_ELECTROMAGNETICS_HPP
+
+namespace mfem
+{
+
+namespace electromagnetics
+{
+
+// Physical Constants
+
+// Permittivity of Free Space (units F/m)
+static const double epsilon0_ = 8.8541878176e-12;
+
+// Permeability of Free Space (units H/m)
+static const double mu0_ = 4.0e-7*M_PI;
+
+} // namespace electromagnetics
+
+} // namespace mfem
+
+#endif // MFEM_ELECTROMAGNETICS_HPP

miniapps/electromagnetics/makefile

@@ -22,7 +22,7 @@ MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
 SEQ_MINIAPPS =
-PAR_MINIAPPS = volta tesla joule
+PAR_MINIAPPS = volta tesla maxwell joule
 ifeq ($(MFEM_USE_MPI),NO)
    MINIAPPS = $(SEQ_MINIAPPS)
 else
@@ -62,6 +62,9 @@ volta-test-par: volta
 tesla-test-par: tesla
 	@$(call mfem-test,$<, $(RUN_MPI), Electromagnetic miniapp,\
 	-maxit 2 -cr '0 0 -0.2 0 0 0.2 0.2 0.4 1')
+maxwell-test-par: maxwell
+	@$(call mfem-test,$<, $(RUN_MPI), Electromagnetic miniapp,\
+	-abcs '-1' -dp '-0.3 0.0 0.0 0.3 0.0 0.0 0.1 1 .5 .5')
 joule-test-par: joule
 	@$(call mfem-test,$<, $(RUN_MPI), Electromagnetic miniapp,\
 	-p rod -tf 3 -m $(SRC)cylinder-hex.mesh)
@@ -79,4 +82,4 @@ clean-build:
 	rm -rf *.dSYM *.TVD.*breakpoints
 
 clean-exec:
-	@rm -rf Volta-AMR* Tesla-AMR* Joule_*
+	@rm -rf Volta-AMR* Tesla-AMR* Maxwell-Parallel* Joule_*

miniapps/electromagnetics/maxwell.cpp

@@ -0,0 +1,523 @@
+// 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.
+//
+//    ------------------------------------------------------------------
+//    Maxwell Miniapp:  Simple Full-Wave Electromagnetic Simulation Code
+//    ------------------------------------------------------------------
+//
+// This miniapp solves a simple 3D full-wave electromagnetic problem using the
+// coupled, first-order equations:
+//
+//                 epsilon dE/dt = Curl 1/mu B - sigma E - J
+//                         dB/dt = - Curl E
+//
+// The permittivity function is that of the vacuum with an optional dielectric
+// sphere. The permeability function is that of the vacuum with an optional
+// diamagnetic or paramagnetic spherical shell. The optional conductivity
+// function is also a user-defined sphere.
+//
+// The optional current density is a pulse of current in the shape of a cylinder
+// with a time dependence resembling the derivative of a Gaussian distribution.
+//
+// Boundary conditions can be 'natural' meaning zero tangential current,
+// 'Dirichlet' which sets the time-derivative of the tangential components of E,
+// or 'absorbing' (we use a simple Sommerfeld first order absorbing boundary
+// condition).
+//
+// We discretize the electric field with H(Curl) finite elements (Nedelec edge
+// elements) and the magnetic flux with H(Div) finite elements (Raviart-Thomas
+// elements).
+//
+// The symplectic time integration algorithm used below is designed to conserve
+// energy unless lossy materials or absorbing boundary conditions are used.
+// When losses are expected, the algorithm uses an implicit method which
+// includes the loss operators in the left hand side of the linear system.
+//
+// For increased accuracy the time integration order can be set to 2, 3, or 4
+// (the default is 1st order).
+//
+// Compile with: make maxwell
+//
+// Sample runs:
+//
+//   Current source in a sphere with absorbing boundary conditions:
+//     mpirun -np 4 maxwell -m ../../data/ball-nurbs.mesh -rs 2
+//                          -abcs '-1'
+//                          -dp '-0.3 0.0 0.0 0.3 0.0 0.0 0.1 1 .5 .5'
+//
+//   Current source in a metal sphere with dielectric and conducting materials:
+//     mpirun -np 4 maxwell -m ../../data/ball-nurbs.mesh -rs 2
+//                          -dbcs '-1'
+//                          -dp '-0.3 0.0 0.0 0.3 0.0 0.0 0.1 1 .5 .5'
+//                          -cs '0.0 0.0 -0.5 .2 3e6'
+//                          -ds '0.0 0.0 0.5 .2 10'
+//
+//   Current source in a metal box:
+//     mpirun -np 4 maxwell -m ../../data/fichera.mesh -rs 3
+//                          -ts 0.25 -tf 10 -dbcs '-1'
+//                          -dp '-0.5 -0.5 0.0 -0.5 -0.5 1.0 0.1 1 .5 1'
+//
+//   Current source with a mixture of absorbing and reflecting boundaries:
+//     mpirun -np 4 maxwell -m ../../data/fichera.mesh -rs 3
+//                          -ts 0.25 -tf 10
+//                          -dp '-0.5 -0.5 0.0 -0.5 -0.5 1.0 0.1 1 .5 1'
+//                          -dbcs '4 8 19 21' -abcs '5 18'
+//
+//   By default the sources and fields are all zero:
+//     mpirun -np 4 maxwell
+
+#include "maxwell_solver.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+using namespace mfem::miniapps;
+using namespace mfem::electromagnetics;
+
+// Permittivity Function
+static Vector ds_params_(0);  // Center, Radius, and Permittivity
+//                               of dielectric sphere
+double dielectric_sphere(const Vector &);
+double epsilon(const Vector &x) { return dielectric_sphere(x); }
+
+// Permeability Function
+static Vector ms_params_(0);  // Center, Inner and Outer Radii, and
+//                               Permeability of magnetic shell
+double magnetic_shell(const Vector &);
+double muInv(const Vector & x) { return 1.0/magnetic_shell(x); }
+
+// Conductivity Function
+static Vector cs_params_(0);  // Center, Radius, and Conductivity
+//                               of conductive sphere
+double conductive_sphere(const Vector &);
+double sigma(const Vector &x) { return conductive_sphere(x); }
+
+// Current Density Function
+static Vector dp_params_(0);  // Axis Start, Axis End, Rod Radius,
+//                               Total Current of Rod, and Frequency
+void dipole_pulse(const Vector &x, double t, Vector &j);
+void j_src(const Vector &x, double t, Vector &j) { dipole_pulse(x, t, j); }
+
+// dE/dt Boundary Condition: The following function returns zero but any time
+// depenent function could be used.
+void dEdtBCFunc(const Vector &x, double t, Vector &E);
+
+// The following functions return zero but they could be modified to set initial
+// conditions for the electric and magnetic fields
+void EFieldFunc(const Vector &, Vector&);
+void BFieldFunc(const Vector &, Vector&);
+
+// Scale factor between input time units and seconds
+static double tScale_ = 1e-9;  // Input time in nanosecond
+
+int SnapTimeStep(double tmax, double dtmax, double & dt);
+
+// Prints the program's logo to the given output stream
+void display_banner(ostream & os);
+
+int main(int argc, char *argv[])
+{
+   MPI_Session mpi(argc, argv);
+
+   if ( mpi.Root() ) { display_banner(cout); }
+
+   // Parse command-line options.
+   const char *mesh_file = "../../data/ball-nurbs.mesh";
+   int sOrder = 1;
+   int tOrder = 1;
+   int serial_ref_levels = 0;
+   int parallel_ref_levels = 0;
+   bool visualization = true;
+   bool visit = true;
+   double dt = 1.0e-12;
+   double dtsf = 0.95;
+   double ti = 0.0;
+   double ts = 1.0;
+   double tf = 40.0;
+
+   Array<int> abcs;
+   Array<int> dbcs;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&sOrder, "-so", "--spatial-order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&tOrder, "-to", "--temporal-order",
+                  "Time integration order.");
+   args.AddOption(&serial_ref_levels, "-rs", "--serial-ref-levels",
+                  "Number of serial refinement levels.");
+   args.AddOption(&parallel_ref_levels, "-rp", "--parallel-ref-levels",
+                  "Number of parallel refinement levels.");
+   args.AddOption(&dtsf, "-sf", "--dt-safety-factor",
+                  "Used to reduce the time step below the upper bound.");
+   args.AddOption(&ti, "-ti", "--initial-time",
+                  "Beginning of time interval to simulate (ns).");
+   args.AddOption(&tf, "-tf", "--final-time",
+                  "End of time interval to simulate (ns).");
+   args.AddOption(&ts, "-ts", "--snapshot-time",
+                  "Time between snapshots (ns).");
+   args.AddOption(&ds_params_, "-ds", "--dielectric-sphere-params",
+                  "Center, Radius, and Permittivity of Dielectric Sphere");
+   args.AddOption(&ms_params_, "-ms", "--magnetic-shell-params",
+                  "Center, Inner Radius, Outer Radius, and Permeability "
+                  "of Magnetic Shell");
+   args.AddOption(&cs_params_, "-cs", "--conductive-sphere-params",
+                  "Center, Radius, and Conductivity of Conductive Sphere");
+   args.AddOption(&dp_params_, "-dp", "--dipole-pulse-params",
+                  "Axis End Points, Radius, Amplitude, "
+                  "Pulse Center (ns), Pulse Width (ns)");
+   args.AddOption(&abcs, "-abcs", "--absorbing-bc-surf",
+                  "Absorbing Boundary Condition Surfaces");
+   args.AddOption(&dbcs, "-dbcs", "--dirichlet-bc-surf",
+                  "Dirichlet Boundary Condition Surfaces");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&visit, "-visit", "--visit", "-no-visit",
+                  "--no-visualization",
+                  "Enable or disable VisIt visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (mpi.Root())
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (mpi.Root())
+   {
+      args.PrintOptions(cout);
+   }
+
+   // Read the (serial) mesh from the given mesh file on all processors.  We can
+   // handle triangular, quadrilateral, tetrahedral, hexahedral, surface and
+   // volume meshes with the same code.
+   Mesh *mesh;
+   ifstream imesh(mesh_file);
+   if (!imesh)
+   {
+      if (mpi.Root())
+      {
+         cerr << "\nCan not open mesh file: " << mesh_file << '\n' << endl;
+      }
+      return 2;
+   }
+   mesh = new Mesh(imesh, 1, 1);
+   imesh.close();
+
+   // Project a NURBS mesh to a piecewise-quadratic curved mesh
+   if (mesh->NURBSext)
+   {
+      mesh->UniformRefinement();
+      if (serial_ref_levels > 0) { serial_ref_levels--; }
+
+      mesh->SetCurvature(2);
+   }
+
+   // Refine the serial mesh on all processors to increase the resolution. In
+   // this example we do 'ref_levels' of uniform refinement.
+   for (int l = 0; l < serial_ref_levels; l++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 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(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+
+   // Refine this mesh in parallel to increase the resolution.
+   for (int l = 0; l < parallel_ref_levels; l++)
+   {
+      pmesh.UniformRefinement();
+   }
+
+   // Create the Electromagnetic solver
+   MaxwellSolver Maxwell(pmesh, sOrder,
+                         ( ds_params_.Size() > 0 ) ? epsilon     : NULL,
+                         ( ms_params_.Size() > 0 ) ? muInv       : NULL,
+                         ( cs_params_.Size() > 0 ) ? sigma       : NULL,
+                         ( dp_params_.Size() > 0 ) ? j_src       : NULL,
+                         abcs, dbcs,
+                         (       dbcs.Size() > 0 ) ? dEdtBCFunc  : NULL
+                        );
+
+   // Display the current number of DoFs in each finite element space
+   Maxwell.PrintSizes();
+
+   // Set the initial conditions for both the electric and magnetic fields
+   VectorFunctionCoefficient EFieldCoef(3,EFieldFunc);
+   VectorFunctionCoefficient BFieldCoef(3,BFieldFunc);
+
+   Maxwell.SetInitialEField(EFieldCoef);
+   Maxwell.SetInitialBField(BFieldCoef);
+
+   // Compute the energy of the initial fields
+   double energy = Maxwell.GetEnergy();
+   if ( mpi.Root() )
+   {
+      cout << "Energy(" << ti << "ns):  " << energy << "J" << endl;
+   }
+
+   // Approximate the largest stable time step
+   double dtmax = Maxwell.GetMaximumTimeStep();
+
+   // Convert times from nanoseconds to seconds
+   ti *= tScale_;
+   tf *= tScale_;
+   ts *= tScale_;
+
+   if ( mpi.Root() )
+   {
+      cout << "Maximum Time Step:     " << dtmax / tScale_ << "ns" << endl;
+   }
+
+   // Round down the time step so that tf-ti is an integer multiple of dt
+   int nsteps = SnapTimeStep(tf-ti, dtsf * dtmax, dt);
+   if ( mpi.Root() )
+   {
+      cout << "Number of Time Steps:  " << nsteps << endl;
+      cout << "Time Step Size:        " << dt / tScale_ << "ns" << endl;
+   }
+
+   // Create the ODE solver
+   SIAVSolver siaSolver(tOrder);
+   siaSolver.Init(Maxwell.GetNegCurl(), Maxwell);
+
+
+   // Initialize GLVis visualization
+   if (visualization)
+   {
+      Maxwell.InitializeGLVis();
+   }
+
+   // Initialize VisIt visualization
+   VisItDataCollection visit_dc("Maxwell-Parallel", &pmesh);
+
+   double t = ti;
+   Maxwell.SetTime(t);
+
+   if ( visit )
+   {
+      Maxwell.RegisterVisItFields(visit_dc);
+   }
+
+   // Write initial fields to disk for VisIt
+   if ( visit )
+   {
+      Maxwell.WriteVisItFields(0);
+   }
+
+   // Send the initial condition by socket to a GLVis server.
+   if (visualization)
+   {
+      Maxwell.DisplayToGLVis();
+   }
+
+   // The main time evolution loop.
+   int it = 1;
+   while (t < tf)
+   {
+      // Run the simulation until a snapshot is needed
+      siaSolver.Run(Maxwell.GetBField(), Maxwell.GetEField(), t, dt,
+                    max(t + dt, ti + ts * it));
+
+      // Approximate the current energy if the fields
+      energy = Maxwell.GetEnergy();
+      if ( mpi.Root() )
+      {
+         cout << "Energy(" << t/tScale_ << "ns):  " << energy << "J" << endl;
+      }
+
+      // Update local DoFs with current true DoFs
+      Maxwell.SyncGridFuncs();
+
+      // Write fields to disk for VisIt
+      if ( visit )
+      {
+         Maxwell.WriteVisItFields(it);
+      }
+
+      // Send the solution by socket to a GLVis server.
+      if (visualization)
+      {
+         Maxwell.DisplayToGLVis();
+      }
+
+      it++;
+   }
+
+   return 0;
+}
+
+// Print the Maxwell ascii logo to the given ostream
+void display_banner(ostream & os)
+{
+   os << "     ___    ____                                      " << endl
+      << "    /   |  /   /                           __   __    " << endl
+      << "   /    |_/ _ /__  ___  _____  _  __ ____ |  | |  |   " << endl
+      << "  /         \\__  \\ \\  \\/  /\\ \\/ \\/ // __ \\|  | |  |   "
+      << endl
+      << " /   /|_/   // __ \\_>    <  \\     /\\  ___/|  |_|  |__ " << endl
+      << "/___/  /_  /(____  /__/\\_ \\  \\/\\_/  \\___  >____/____/ " << endl
+      << "         \\/       \\/      \\/             \\/           " << endl
+      << flush;
+}
+
+// A sphere with constant permittivity.  The sphere has a radius, center, and
+// permittivity specified on the command line and stored in ds_params_.
+double dielectric_sphere(const Vector &x)
+{
+   double r2 = 0.0;
+
+   for (int i=0; i<x.Size(); i++)
+   {
+      r2 += (x(i)-ds_params_(i))*(x(i)-ds_params_(i));
+   }
+
+   if ( sqrt(r2) <= ds_params_(x.Size()) )
+   {
+      return ds_params_(x.Size()+1) * epsilon0_;
+   }
+   return epsilon0_;
+}
+
+// A spherical shell with constant permeability.  The sphere has inner and outer
+// radii, center, and relative permeability specified on the command line and
+// stored in ms_params_.
+double magnetic_shell(const Vector &x)
+{
+   double r2 = 0.0;
+
+   for (int i=0; i<x.Size(); i++)
+   {
+      r2 += (x(i)-ms_params_(i))*(x(i)-ms_params_(i));
+   }
+
+   if ( sqrt(r2) >= ms_params_(x.Size()) &&
+        sqrt(r2) <= ms_params_(x.Size()+1) )
+   {
+      return mu0_*ms_params_(x.Size()+2);
+   }
+   return mu0_;
+}
+
+// A sphere with constant conductivity.  The sphere has a radius, center, and
+// conductivity specified on the command line and stored in ls_params_.
+double conductive_sphere(const Vector &x)
+{
+   double r2 = 0.0;
+
+   for (int i=0; i<x.Size(); i++)
+   {
+      r2 += (x(i)-cs_params_(i))*(x(i)-cs_params_(i));
+   }
+
+   if ( sqrt(r2) <= cs_params_(x.Size()) )
+   {
+      return cs_params_(x.Size()+1);
+   }
+   return 0.0;
+}
+
+// A cylindrical rod of current density.  The rod has two axis end points, a
+// radus, a current amplitude in Amperes, a center time, and a width.  All of
+// these parameters are stored in dp_params_.
+void dipole_pulse(const Vector &x, double t, Vector &j)
+{
+   MFEM_ASSERT(x.Size() == 3, "current source requires 3D space.");
+
+   j.SetSize(x.Size());
+   j = 0.0;
+
+   Vector  v(x.Size());  // Normalized Axis vector
+   Vector xu(x.Size());  // x vector relative to the axis end-point
+
+   xu = x;
+
+   for (int i=0; i<x.Size(); i++)
+   {
+      xu[i] -= dp_params_[i];
+      v[i]   = dp_params_[x.Size()+i] - dp_params_[i];
+   }
+
+   double h = v.Norml2();
+
+   if ( h == 0.0 )
+   {
+      return;
+   }
+   v /= h;
+
+   double r = dp_params_[2*x.Size()+0];
+   double a = dp_params_[2*x.Size()+1] * tScale_;
+   double b = dp_params_[2*x.Size()+2] * tScale_;
+   double c = dp_params_[2*x.Size()+3] * tScale_;
+
+   double xv = xu * v;
+
+   // Compute perpendicular vector from axis to x
+   xu.Add(-xv, v);
+
+   double xp = xu.Norml2();
+
+   if ( xv >= 0.0 && xv <= h && xp <= r )
+   {
+      j = v;
+   }
+
+   j *= a * (t - b) * exp(-0.5 * pow((t-b)/c, 2)) / (c * c);
+}
+
+void
+EFieldFunc(const Vector &x, Vector &E)
+{
+   E.SetSize(3);
+   E = 0.0;
+}
+
+void
+BFieldFunc(const Vector &x, Vector &B)
+{
+   B.SetSize(3);
+   B = 0.0;
+}
+
+void
+dEdtBCFunc(const Vector &x, double t, Vector &dE)
+{
+   dE.SetSize(3);
+   dE = 0.0;
+}
+
+int
+SnapTimeStep(double tmax, double dtmax, double & dt)
+{
+   double dsteps = tmax/dtmax;
+
+   int nsteps = pow(10,(int)ceil(log10(dsteps)));
+
+   for (int i=1; i<=5; i++)
+   {
+      int a = (int)ceil(log10(dsteps/pow(5.0,i)));
+      int nstepsi = (int)pow(5,i)*max(1,(int)pow(10,a));
+
+      nsteps = min(nsteps,nstepsi);
+   }
+
+   dt = tmax / nsteps;
+
+   return nsteps;
+}

miniapps/electromagnetics/maxwell_solver.cpp

@@ -0,0 +1,669 @@
+// 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.
+
+#include "maxwell_solver.hpp"
+
+#ifdef MFEM_USE_MPI
+
+using namespace std;
+
+namespace mfem
+{
+
+using namespace miniapps;
+
+namespace electromagnetics
+{
+
+// Used for combining scalar coefficients
+double prodFunc(double a, double b) { return a * b; }
+
+MaxwellSolver::MaxwellSolver(ParMesh & pmesh, int order,
+                             double (*eps     )(const Vector&),
+                             double (*muInv   )(const Vector&),
+                             double (*sigma   )(const Vector&),
+                             void   (*j_src   )(const Vector&, double, Vector&),
+                             Array<int> & abcs,
+                             Array<int> & dbcs,
+                             void   (*dEdt_bc )(const Vector&, double, Vector&))
+   : myid_(0),
+     num_procs_(1),
+     order_(order),
+     logging_(1),
+     dtMax_(-1.0),
+     dtScale_(1.0e6),
+     pmesh_(&pmesh),
+     HCurlFESpace_(NULL),
+     HDivFESpace_(NULL),
+     hDivMassMuInv_(NULL),
+     hCurlLosses_(NULL),
+     weakCurlMuInv_(NULL),
+     Curl_(NULL),
+     e_(NULL),
+     b_(NULL),
+     j_(NULL),
+     dedt_(NULL),
+     rhs_(NULL),
+     jd_(NULL),
+     M1Losses_(NULL),
+     M2MuInv_(NULL),
+     NegCurl_(NULL),
+     WeakCurlMuInv_(NULL),
+     E_(NULL),
+     B_(NULL),
+     HD_(NULL),
+     RHS_(NULL),
+     epsCoef_(NULL),
+     muInvCoef_(NULL),
+     sigmaCoef_(NULL),
+     etaInvCoef_(NULL),
+     eCoef_(NULL),
+     bCoef_(NULL),
+     jCoef_(NULL),
+     dEdtBCCoef_(NULL),
+     eps_(eps),
+     muInv_(muInv),
+     sigma_(sigma),
+     j_src_(j_src),
+     dEdt_bc_(dEdt_bc),
+     visit_dc_(NULL)
+{
+   // Initialize MPI variables
+   MPI_Comm_size(pmesh_->GetComm(), &num_procs_);
+   MPI_Comm_rank(pmesh_->GetComm(), &myid_);
+
+   // Define compatible parallel finite element spaces on the parallel
+   // mesh. Here we use arbitrary order H1, Nedelec, and Raviart-Thomas finite
+   // elements.
+   HCurlFESpace_ = new ND_ParFESpace(pmesh_,order,pmesh_->Dimension());
+   HDivFESpace_  = new RT_ParFESpace(pmesh_,order,pmesh_->Dimension());
+
+   this->height = HCurlFESpace_->GlobalTrueVSize();
+   this->width  = HDivFESpace_->GlobalTrueVSize();
+
+   // Check for absorbing materials or boundaries
+   lossy_ = abcs.Size() > 0 || sigma_ != NULL;
+
+   // Require implicit handling of loss terms
+   type = lossy_ ? IMPLICIT : EXPLICIT;
+
+   // Electric permittivity
+   if ( eps_ == NULL )
+   {
+      epsCoef_ = new ConstantCoefficient(epsilon0_);
+   }
+   else
+   {
+      if ( myid_ == 0 && logging_ > 0 )
+      {
+         cout << "Creating Permittivity Coefficient" << endl;
+      }
+      epsCoef_ = new FunctionCoefficient(eps_);
+   }
+
+   // Inverse of the magnetic permeability
+   if ( muInv_ == NULL )
+   {
+      muInvCoef_ = new ConstantCoefficient(1.0/mu0_);
+   }
+   else
+   {
+      if ( myid_ == 0 && logging_ > 0 )
+      {
+         cout << "Creating Permeability Coefficient" << endl;
+      }
+      muInvCoef_ = new FunctionCoefficient(muInv_);
+   }
+
+   // Electric conductivity
+   if ( sigma_ != NULL )
+   {
+      if ( myid_ == 0 && logging_ > 0 )
+      {
+         cout << "Creating Conductivity Coefficient" << endl;
+      }
+      sigmaCoef_ = new FunctionCoefficient(sigma_);
+   }
+
+   // Impedance of free space
+   if ( abcs.Size() > 0 )
+   {
+      if ( myid_ == 0 && logging_ > 0 )
+      {
+         cout << "Creating Admittance Coefficient" << endl;
+      }
+
+      abc_marker_.SetSize(pmesh.bdr_attributes.Max());
+      if ( abcs.Size() == 1 && abcs[0] < 0 )
+      {
+         // Mark all boundaries as absorbing
+         abc_marker_ = 1;
+      }
+      else
+      {
+         // Mark select boundaries as absorbing
+         abc_marker_ = 0;
+         for (int i=0; i<abcs.Size(); i++)
+         {
+            abc_marker_[abcs[i]-1] = 1;
+         }
+      }
+      etaInvCoef_ = new ConstantCoefficient(sqrt(epsilon0_/mu0_));
+   }
+
+   // Electric Field Boundary Condition
+   if ( dbcs.Size() > 0 )
+   {
+      if ( myid_ == 0 && logging_ > 0 )
+      {
+         cout << "Configuring Dirichlet BC" << endl;
+      }
+      dbc_marker_.SetSize(pmesh.bdr_attributes.Max());
+      if ( dbcs.Size() == 1 && dbcs[0] < 0 )
+      {
+         // Mark all boundaries as Dirichlet
+         dbc_marker_ = 1;
+      }
+      else
+      {
+         // Mark select boundaries as Dirichlet
+         dbc_marker_ = 0;
+         for (int i=0; i<dbcs.Size(); i++)
+         {
+            dbc_marker_[dbcs[i]-1] = 1;
+         }
+      }
+
+      HCurlFESpace_->GetEssentialTrueDofs(dbc_marker_, dbc_dofs_);
+
+      if ( dEdt_bc_ != NULL )
+      {
+         dEdtBCCoef_ = new VectorFunctionCoefficient(3,dEdt_bc_);
+      }
+      else
+      {
+         Vector ebc(3); ebc = 0.0;
+         dEdtBCCoef_ = new VectorConstantCoefficient(ebc);
+      }
+   }
+
+   // Bilinear Forms
+   if ( myid_ == 0 && logging_ > 0 )
+   {
+      cout << "Creating H(Div) Mass Operator" << endl;
+   }
+   hDivMassMuInv_ = new ParBilinearForm(HDivFESpace_);
+   hDivMassMuInv_->AddDomainIntegrator(new VectorFEMassIntegrator(*muInvCoef_));
+
+   if ( myid_ == 0 && logging_ > 0 )
+   {
+      cout << "Creating Weak Curl Operator" << endl;
+   }
+   weakCurlMuInv_ = new ParMixedBilinearForm(HDivFESpace_,HCurlFESpace_);
+   weakCurlMuInv_->AddDomainIntegrator(
+      new MixedVectorWeakCurlIntegrator(*muInvCoef_));
+
+   // Assemble Matrices
+   hDivMassMuInv_->Assemble();
+   weakCurlMuInv_->Assemble();
+
+   hDivMassMuInv_->Finalize();
+   weakCurlMuInv_->Finalize();
+
+   if ( sigmaCoef_ || etaInvCoef_ )
+   {
+      if ( myid_ == 0 && logging_ > 0 )
+      {
+         cout << "Creating H(Curl) Loss Operator" << endl;
+      }
+      hCurlLosses_  = new ParBilinearForm(HCurlFESpace_);
+      if ( sigmaCoef_ )
+      {
+         if ( myid_ == 0 && logging_ > 0 )
+         {
+            cout << "Adding domain integrator for conductive regions" << endl;
+         }
+         hCurlLosses_->AddDomainIntegrator(
+            new VectorFEMassIntegrator(*sigmaCoef_));
+      }
+      if ( etaInvCoef_ )
+      {
+         if ( myid_ == 0 && logging_ > 0 )
+         {
+            cout << "Adding boundary integrator for absorbing boundary" << endl;
+         }
+         hCurlLosses_->AddBoundaryIntegrator(
+            new VectorFEMassIntegrator(*etaInvCoef_), abc_marker_);
+      }
+      hCurlLosses_->Assemble();
+      hCurlLosses_->Finalize();
+      M1Losses_ = hCurlLosses_->ParallelAssemble();
+   }
+
+   // Create Linear Algebra Matrices
+   M2MuInv_ = hDivMassMuInv_->ParallelAssemble();
+   WeakCurlMuInv_ = weakCurlMuInv_->ParallelAssemble();
+
+   if ( myid_ == 0 && logging_ > 0 )
+   {
+      cout << "Creating discrete curl operator" << endl;
+   }
+   Curl_ = new ParDiscreteCurlOperator(HCurlFESpace_, HDivFESpace_);
+   Curl_->Assemble();
+   Curl_->Finalize();
+   NegCurl_ = Curl_->ParallelAssemble();
+   // Beware this modifies the matrix stored within the Curl_ object.
+   *NegCurl_ *= -1.0;
+
+   // Build grid functions
+   e_    = new ParGridFunction(HCurlFESpace_);
+   dedt_ = new ParGridFunction(HCurlFESpace_);
+   rhs_  = new ParGridFunction(HCurlFESpace_);
+   b_    = new ParGridFunction(HDivFESpace_);
+
+   E_ = e_->ParallelProject();
+   B_ = b_->ParallelProject();
+
+   HD_  = new HypreParVector(HDivFESpace_);
+   RHS_ = new HypreParVector(HCurlFESpace_);
+
+   // Initialize dedt to zero
+   *dedt_ = 0.0;
+
+   if ( j_src_)
+   {
+      if ( myid_ == 0 && logging_ > 0 )
+      {
+         cout << "Creating Current Source" << endl;
+      }
+      jCoef_ = new VectorFunctionCoefficient(3,j_src_);
+
+      j_  = new ParGridFunction(HCurlFESpace_);
+      j_->ProjectCoefficient(*jCoef_);
+
+      jd_ = new ParLinearForm(HCurlFESpace_);
+      jd_->AddDomainIntegrator(new VectorFEDomainLFIntegrator(*jCoef_));
+      jd_->Assemble();
+   }
+
+   dtMax_ = GetMaximumTimeStep();
+}
+
+MaxwellSolver::~MaxwellSolver()
+{
+   delete epsCoef_;
+   delete muInvCoef_;
+   delete etaInvCoef_;
+   delete jCoef_;
+   delete dEdtBCCoef_;
+
+   delete E_;
+   delete B_;
+   delete HD_;
+   delete RHS_;
+
+   delete e_;
+   delete b_;
+   delete j_;
+   delete dedt_;
+   delete rhs_;
+   delete jd_;
+
+   delete Curl_;
+
+   delete M1Losses_;
+   delete M2MuInv_;
+
+   delete hDivMassMuInv_;
+   delete hCurlLosses_;
+   delete weakCurlMuInv_;
+
+   delete HCurlFESpace_;
+   delete HDivFESpace_;
+
+   map<int, ParBilinearForm*>::iterator mit1;
+   for (mit1=a1_.begin(); mit1!=a1_.end(); mit1++)
+   {
+      int i = mit1->first;
+      delete pcg_[i];
+      delete diagScale_[i];
+      delete A1_[i];
+      delete a1_[i];
+   }
+
+   map<int, Coefficient*>::iterator mit2;
+   for (mit2=dtCoef_.begin(); mit2!=dtCoef_.end(); mit2++)
+   {
+      delete mit2->second;
+   }
+   for (mit2=dtSigmaCoef_.begin(); mit2!=dtSigmaCoef_.end(); mit2++)
+   {
+      delete mit2->second;
+   }
+   for (mit2=dtEtaInvCoef_.begin(); mit2!=dtEtaInvCoef_.end(); mit2++)
+   {
+      delete mit2->second;
+   }
+
+   map<string, socketstream*>::iterator mit3;
+   for (mit3=socks_.begin(); mit3!=socks_.end(); mit3++)
+   {
+      delete mit3->second;
+   }
+}
+
+HYPRE_Int
+MaxwellSolver::GetProblemSize()
+{
+   return HCurlFESpace_->GlobalTrueVSize();
+}
+
+void
+MaxwellSolver::PrintSizes()
+{
+   HYPRE_Int size_nd = HCurlFESpace_->GlobalTrueVSize();
+   HYPRE_Int size_rt = HDivFESpace_->GlobalTrueVSize();
+   if ( myid_ == 0 )
+   {
+      cout << "Number of H(Curl) unknowns: " << size_nd << endl;
+      cout << "Number of H(Div)  unknowns: " << size_rt << endl << flush;
+   }
+}
+
+void
+MaxwellSolver::SetInitialEField(VectorCoefficient & EFieldCoef)
+{
+   eCoef_ = &EFieldCoef;
+   e_->ProjectCoefficient(EFieldCoef);
+   e_->ParallelProject(*E_);
+}
+
+void
+MaxwellSolver::SetInitialBField(VectorCoefficient & BFieldCoef)
+{
+   bCoef_ = &BFieldCoef;
+   b_->ProjectCoefficient(BFieldCoef);
+   b_->ParallelProject(*B_);
+}
+
+void
+MaxwellSolver::Mult(const Vector &B, Vector &dEdt) const
+{
+   implicitSolve(0.0, B, dEdt);
+}
+
+void
+MaxwellSolver::ImplicitSolve(double dt, const Vector &B, Vector &dEdt)
+{
+   implicitSolve(dt, B, dEdt);
+}
+
+void
+MaxwellSolver::setupSolver(const int idt, const double dt) const
+{
+   if ( pcg_.find(idt) == pcg_.end() )
+   {
+      if ( myid_ == 0 && logging_ > 0 )
+      {
+         cout << "Creating implicit operator for dt = " << dt << endl;
+      }
+
+      a1_[idt] = new ParBilinearForm(HCurlFESpace_);
+      a1_[idt]->AddDomainIntegrator(
+         new VectorFEMassIntegrator(epsCoef_));
+
+      if ( idt != 0 )
+      {
+         dtCoef_[idt] = new ConstantCoefficient(0.5 * dt);
+         if ( sigmaCoef_ )
+         {
+            dtSigmaCoef_[idt] = new TransformedCoefficient(dtCoef_[idt],
+                                                           sigmaCoef_,
+                                                           prodFunc);
+            a1_[idt]->AddDomainIntegrator(
+               new VectorFEMassIntegrator(dtSigmaCoef_[idt]));
+         }
+         if ( etaInvCoef_ )
+         {
+            dtEtaInvCoef_[idt] = new TransformedCoefficient(dtCoef_[idt],
+                                                            etaInvCoef_,
+                                                            prodFunc);
+            a1_[idt]->AddBoundaryIntegrator(
+               new VectorFEMassIntegrator(dtEtaInvCoef_[idt]),
+               const_cast<Array<int>&>(abc_marker_));
+         }
+      }
+
+      a1_[idt]->Assemble();
+      a1_[idt]->Finalize();
+      A1_[idt] = a1_[idt]->ParallelAssemble();
+
+      diagScale_[idt] = new HypreDiagScale(*A1_[idt]);
+      pcg_[idt] = new HyprePCG(*A1_[idt]);
+      pcg_[idt]->SetTol(1.0e-12);
+      pcg_[idt]->SetMaxIter(200);
+      pcg_[idt]->SetPrintLevel(0);
+      pcg_[idt]->SetPreconditioner(*diagScale_[idt]);
+   }
+}
+
+void
+MaxwellSolver::implicitSolve(double dt, const Vector &B, Vector &dEdt) const
+{
+   int idt = hCurlLosses_ ? ((int)(dtScale_ * dt / dtMax_)) : 0;
+
+   b_->Distribute(B);
+   weakCurlMuInv_->Mult(*b_, *rhs_);
+
+   if ( hCurlLosses_ )
+   {
+      e_->Distribute(*E_);
+      hCurlLosses_->AddMult(*e_, *rhs_, -1.0);
+   }
+
+   if ( jd_ )
+   {
+      jCoef_->SetTime(t); // 't' is member data from mfem::TimeDependentOperator
+      jd_->Assemble();
+      *rhs_ -= *jd_;
+   }
+
+   if ( dEdtBCCoef_ )
+   {
+      dEdtBCCoef_->SetTime(t);
+      dedt_->ProjectBdrCoefficientTangent(*dEdtBCCoef_,
+                                          const_cast<Array<int>&>(dbc_marker_));
+   }
+
+   // Create objects and matrices for solving with the given time step
+   setupSolver(idt, dt);
+
+   // Apply essential BCs and determine true DoFs for the right hand side
+   a1_[idt]->FormLinearSystem(dbc_dofs_, *dedt_, *rhs_, *A1_[idt], dEdt, *RHS_);
+
+   // Solve for the time derivative of the electric field (true DoFs)
+   pcg_[idt]->Mult(*RHS_, dEdt);
+
+   // Distribute shared DoFs to relevant processors
+   a1_[idt]->RecoverFEMSolution(dEdt, *rhs_, *dedt_);
+}
+
+void
+MaxwellSolver::SyncGridFuncs()
+{
+   e_->Distribute(*E_);
+   b_->Distribute(*B_);
+}
+
+double
+MaxwellSolver::GetMaximumTimeStep() const
+{
+   if ( dtMax_ > 0.0 )
+   {
+      return dtMax_;
+   }
+
+   HypreParVector * v0 = new HypreParVector(HCurlFESpace_);
+   HypreParVector * v1 = new HypreParVector(HCurlFESpace_);
+   HypreParVector * u0 = new HypreParVector(HDivFESpace_);
+
+   v0->Randomize(1234);
+
+   int iter = 0, nstep = 20;
+   double dt0 = 1.0, dt1 = 1.0, change = 1.0, ptol = 0.001;
+
+   // Create Solver assuming no loss operators
+   setupSolver(0, 0.0);
+
+   // Use power method to approximate the largest eigenvalue of the update
+   // operator.
+   while ( iter < nstep && change > ptol )
+   {
+      double normV0 = InnerProduct(*v0,*v0);
+      *v0 /= sqrt(normV0);
+
+      NegCurl_->Mult(*v0,*u0);
+      M2MuInv_->Mult(*u0,*HD_);
+      NegCurl_->MultTranspose(*HD_,*RHS_);
+
+      pcg_[0]->Mult(*RHS_,*v1);
+
+      double lambda = InnerProduct(*v0,*v1);
+      dt1 = 2.0/sqrt(lambda);
+      change = fabs((dt1-dt0)/dt0);
+      dt0 = dt1;
+
+      if ( myid_ == 0 && logging_ > 1 )
+      {
+         cout << iter << ":  " << dt0 << " " << change << endl;
+      }
+
+      std::swap(v0, v1);
+
+      iter++;
+   }
+
+   delete v0;
+   delete v1;
+   delete u0;
+
+   return dt0;
+}
+
+double
+MaxwellSolver::GetEnergy() const
+{
+   double energy = 0.0;
+
+   A1_[0]->Mult(*E_,*RHS_);
+   M2MuInv_->Mult(*B_,*HD_);
+
+   energy = InnerProduct(*E_,*RHS_) + InnerProduct(*B_,*HD_);
+
+   return 0.5 * energy;
+}
+
+void
+MaxwellSolver::RegisterVisItFields(VisItDataCollection & visit_dc)
+{
+   visit_dc_ = &visit_dc;
+
+   visit_dc.RegisterField("E", e_);
+   visit_dc.RegisterField("B", b_);
+   if ( j_ )
+   {
+      visit_dc.RegisterField("J", j_);
+   }
+}
+
+void
+MaxwellSolver::WriteVisItFields(int it)
+{
+   if ( visit_dc_ )
+   {
+      if ( myid_ == 0 && logging_ > 1 )
+      { cout << "Writing VisIt files ..." << flush; }
+
+      if ( j_ )
+      {
+         jCoef_->SetTime(t);
+         j_->ProjectCoefficient(*jCoef_);
+      }
+
+      visit_dc_->SetCycle(it);
+      visit_dc_->SetTime(t);
+      visit_dc_->Save();
+
+      if ( myid_ == 0 && logging_ > 1 ) { cout << " " << endl << flush; }
+   }
+}
+
+void
+MaxwellSolver::InitializeGLVis()
+{
+   if ( myid_ == 0 && logging_ > 0 )
+   { cout << "Opening GLVis sockets." << endl << flush; }
+
+   socks_["E"] = new socketstream;
+   socks_["E"]->precision(8);
+
+   socks_["B"] = new socketstream;
+   socks_["B"]->precision(8);
+
+   if ( j_ )
+   {
+      socks_["J"] = new socketstream;
+      socks_["J"]->precision(8);
+   }
+
+   if ( myid_ == 0 && logging_ > 0 )
+   { cout << "GLVis sockets open." << endl << flush; }
+}
+
+void
+MaxwellSolver::DisplayToGLVis()
+{
+   if ( myid_ == 0 && logging_ > 1 )
+   { cout << "Sending data to GLVis ..." << flush; }
+
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+
+   int Wx = 0, Wy = 0; // window position
+   int Ww = 350, Wh = 350; // window size
+   int offx = Ww+10, offy = Wh+45; // window offsets
+
+   VisualizeField(*socks_["E"], vishost, visport,
+                  *e_, "Electric Field (E)", Wx, Wy, Ww, Wh);
+   Wx += offx;
+
+   VisualizeField(*socks_["B"], vishost, visport,
+                  *b_, "Magnetic Flux Density (B)", Wx, Wy, Ww, Wh);
+
+   if ( j_ )
+   {
+      Wx = 0;
+      Wy += offy;
+
+      jCoef_->SetTime(t); // Is member data from mfem::TimeDependentOperator
+      j_->ProjectCoefficient(*jCoef_);
+
+      VisualizeField(*socks_["J"], vishost, visport,
+                     *j_, "Current Density (J)", Wx, Wy, Ww, Wh);
+   }
+   if ( myid_ == 0 && logging_ > 1 ) { cout << " " << flush; }
+}
+
+} // namespace electromagnetics
+
+} // namespace mfem
+
+#endif // MFEM_USE_MPI

miniapps/electromagnetics/maxwell_solver.hpp

@@ -0,0 +1,176 @@
+// 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.
+
+#ifndef MFEM_MAXWELL_SOLVER
+#define MFEM_MAXWELL_SOLVER
+
+#include "../common/pfem_extras.hpp"
+#include "electromagnetics.hpp"
+
+#ifdef MFEM_USE_MPI
+
+#include <string>
+#include <map>
+
+using namespace std;
+using namespace mfem;
+
+namespace mfem
+{
+
+using miniapps::ND_ParFESpace;
+using miniapps::RT_ParFESpace;
+using miniapps::ParDiscreteCurlOperator;
+
+namespace electromagnetics
+{
+
+class MaxwellSolver : public TimeDependentOperator
+{
+public:
+   MaxwellSolver(ParMesh & pmesh, int sOrder,
+                 double (*eps     )(const Vector&),
+                 double (*muInv   )(const Vector&),
+                 double (*sigma   )(const Vector&),
+                 void   (*j_src   )(const Vector&, double, Vector&),
+                 Array<int> & abcs, Array<int> & dbcs,
+                 void   (*dEdt_bc )(const Vector&, double, Vector&));
+
+   ~MaxwellSolver();
+
+   int GetLogging() const { return logging_; }
+   void SetLogging(int logging) { logging_ = logging; }
+
+   HYPRE_Int GetProblemSize();
+
+   void PrintSizes();
+
+   void SetInitialEField(VectorCoefficient & EFieldCoef);
+   void SetInitialBField(VectorCoefficient & BFieldCoef);
+
+   void Mult(const Vector &B, Vector &dEdt) const;
+
+   void ImplicitSolve(const double dt, const Vector &x, Vector &k);
+
+   double GetMaximumTimeStep() const;
+
+   double GetEnergy() const;
+
+   Operator & GetNegCurl() { return *NegCurl_; }
+
+   Vector & GetEField() { return *E_; }
+   Vector & GetBField() { return *B_; }
+
+   void SyncGridFuncs();
+
+   void RegisterVisItFields(VisItDataCollection & visit_dc);
+
+   void WriteVisItFields(int it = 0);
+
+   void InitializeGLVis();
+
+   void DisplayToGLVis();
+
+private:
+
+   // This method alters mutable member data
+   void setupSolver(const int idt, const double dt) const;
+
+   void implicitSolve(const double dt, const Vector &x, Vector &k) const;
+
+   int myid_;
+   int num_procs_;
+   int order_;
+   int logging_;
+
+   bool lossy_;
+
+   double dtMax_;   // Maximum stable time step
+   double dtScale_; // Used to scale dt before converting to an integer
+
+   ParMesh * pmesh_;
+
+   ND_ParFESpace * HCurlFESpace_;
+   RT_ParFESpace * HDivFESpace_;
+
+   ParBilinearForm * hDivMassMuInv_;
+   ParBilinearForm * hCurlLosses_;
+   ParMixedBilinearForm * weakCurlMuInv_;
+
+   ParDiscreteCurlOperator * Curl_;
+
+   ParGridFunction * e_;    // Electric Field (HCurl)
+   ParGridFunction * b_;    // Magnetic Flux (HDiv)
+   ParGridFunction * j_;    // Volumetric Current Density (HCurl)
+   ParGridFunction * dedt_; // Time Derivative of Electric Field (HCurl)
+   ParGridFunction * rhs_;  // Dual of displacement current, rhs vector (HCurl)
+   ParLinearForm   * jd_;   // Dual of current density (HCurl)
+
+   HypreParMatrix * M1Losses_;
+   HypreParMatrix * M2MuInv_;
+   HypreParMatrix * NegCurl_;
+   HypreParMatrix * WeakCurlMuInv_;
+   HypreParVector * E_; // Current value of the electric field DoFs
+   HypreParVector * B_; // Current value of the magnetic flux DoFs
+   mutable HypreParVector * HD_; // Used in energy calculation
+   mutable HypreParVector * RHS_;
+
+   Coefficient       * epsCoef_;    // Electric Permittivity Coefficient
+   Coefficient       * muInvCoef_;  // Magnetic Permeability Coefficient
+   Coefficient       * sigmaCoef_;  // Electric Conductivity Coefficient
+   Coefficient       * etaInvCoef_; // Admittance Coefficient
+   VectorCoefficient * eCoef_;      // Initial Electric Field
+   VectorCoefficient * bCoef_;      // Initial Magnetic Flux
+   VectorCoefficient * jCoef_;      // Time dependent current density
+   VectorCoefficient * dEdtBCCoef_; // Time dependent boundary condition
+
+   double (*eps_    )(const Vector&);
+   double (*muInv_  )(const Vector&);
+   double (*sigma_  )(const Vector&);
+   void   (*j_src_  )(const Vector&, double, Vector&);
+
+   // Array of 0's and 1's marking the location of absorbing surfaces
+   Array<int> abc_marker_;
+
+   // Array of 0's and 1's marking the location of Dirichlet boundaries
+   Array<int> dbc_marker_;
+   void   (*dEdt_bc_)(const Vector&, double, Vector&);
+
+   // Dirichlet degrees of freedom
+   Array<int>   dbc_dofs_;
+
+   // High order symplectic integration requires partial time steps of differing
+   // lengths. If losses are present the system matrix includes a portion scaled
+   // by the time step. Consequently, high order time integration requires
+   // different system matrices. The following maps contain various objects that
+   // depend on the time step.
+   mutable std::map<int, ParBilinearForm *> a1_;
+   mutable std::map<int, HypreParMatrix  *> A1_;
+   mutable std::map<int, Coefficient     *> dtCoef_;
+   mutable std::map<int, Coefficient     *> dtSigmaCoef_;
+   mutable std::map<int, Coefficient     *> dtEtaInvCoef_;
+   mutable std::map<int, HypreDiagScale  *> diagScale_;
+   mutable std::map<int, HyprePCG        *> pcg_;
+
+   // Data collection used to write VisIt files
+   VisItDataCollection * visit_dc_;
+
+   // Sockets used to communicate with GLVis
+   std::map<std::string, socketstream*> socks_;
+};
+
+} // namespace electromagnetics
+
+} // namespace mfem
+
+#endif // MFEM_USE_MPI
+
+#endif // MFEM_MAXWELL_SOLVER

miniapps/electromagnetics/tesla_solver.hpp

@@ -13,6 +13,7 @@
 #define MFEM_TESLA_SOLVER
 
 #include "../common/pfem_extras.hpp"
+#include "electromagnetics.hpp"
 
 #ifdef MFEM_USE_MPI
 
@@ -32,10 +33,6 @@ using miniapps::DivergenceFreeProjector;
 namespace electromagnetics
 {
 
-// Physical Constants
-// Permeability of Free Space (units H/m)
-static double mu0_ = 4.0e-7*M_PI;
-
 class SurfaceCurrent;
 class TeslaSolver
 {

miniapps/electromagnetics/volta.cpp

@@ -65,10 +65,6 @@ using namespace std;
 using namespace mfem;
 using namespace mfem::electromagnetics;
 
-// Physical Constants
-// Permittivity of Free Space (units F/m)
-static double epsilon0_ = 8.8541878176e-12;
-
 // Permittivity Functions
 Coefficient * SetupPermittivityCoefficient();
 

miniapps/electromagnetics/volta_solver.hpp

@@ -13,6 +13,7 @@
 #define MFEM_VOLTA_SOLVER
 
 #include "../common/pfem_extras.hpp"
+#include "electromagnetics.hpp"
 
 #ifdef MFEM_USE_MPI
 

miniapps/nurbs/CMakeLists.txt

@@ -0,0 +1,37 @@
+# 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.
+
+add_mfem_miniapp(nurbs_ex1
+  MAIN ex1.cpp
+  LIBRARIES mfem)
+
+add_test(NAME nurbs_ex1_ser
+  COMMAND $<TARGET_FILE:nurbs_ex1> -no-vis)
+
+if (MFEM_USE_MPI)
+  add_mfem_miniapp(nurbs_ex1p
+    MAIN ex1p.cpp
+    LIBRARIES mfem)
+
+  add_test(NAME nurbs_ex1p_np=4
+    COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4
+    ${MPIEXEC_PREFLAGS} $<TARGET_FILE:nurbs_ex1p> -no-vis
+    ${MPIEXEC_POSTFLAGS})
+
+  add_mfem_miniapp(nurbs_ex11p
+    MAIN ex11p.cpp
+    LIBRARIES mfem)
+
+  add_test(NAME nurbs_ex11p_np=4
+    COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4
+    ${MPIEXEC_PREFLAGS} $<TARGET_FILE:nurbs_ex11p> -no-vis
+    ${MPIEXEC_POSTFLAGS})
+endif()

miniapps/nurbs/ex1.cpp

@@ -0,0 +1,232 @@
+//                          MFEM Example 1 - NURBS Version
+//
+// Compile with: make ex1
+//
+// Sample runs:  ex1 -m ../../data/square-disc.mesh
+//               ex1 -m ../../data/star.mesh
+//               ex1 -m ../../data/escher.mesh
+//               ex1 -m ../../data/fichera.mesh
+//               ex1 -m ../../data/square-disc-p2.vtk -o 2
+//               ex1 -m ../../data/square-disc-p3.mesh -o 3
+//               ex1 -m ../../data/square-disc-nurbs.mesh -o -1
+//               ex1 -m ../../data/disc-nurbs.mesh -o -1
+//               ex1 -m ../../data/pipe-nurbs.mesh -o -1
+//               ex1 -m ../../data/star-surf.mesh
+//               ex1 -m ../../data/square-disc-surf.mesh
+//               ex1 -m ../../data/inline-segment.mesh
+//               ex1 -m ../../data/amr-quad.mesh
+//               ex1 -m ../../data/amr-hex.mesh
+//               ex1 -m ../../data/fichera-amr.mesh
+//               ex1 -m ../../data/mobius-strip.mesh
+//               ex1 -m ../../data/mobius-strip.mesh -o -1 -sc
+//
+// Description:  This example code demonstrates the use of MFEM to define a
+//               simple finite element discretization of the Laplace problem
+//               -Delta u = 1 with homogeneous Dirichlet boundary conditions.
+//               Specifically, we discretize using a FE space of the specified
+//               order, or if order < 1 using an isoparametric/isogeometric
+//               space (i.e. quadratic for quadratic curvilinear mesh, NURBS for
+//               NURBS mesh, etc.)
+//
+//               The example highlights the use of mesh refinement, finite
+//               element grid functions, as well as linear and bilinear forms
+//               corresponding to the left-hand side and right-hand side of the
+//               discrete linear system. We also cover the explicit elimination
+//               of essential boundary conditions, static condensation, and the
+//               optional connection to the GLVis tool for visualization.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../../data/star.mesh";
+   bool static_cond = false;
+   bool visualization = 1;
+   Array<int> order(1);
+   order[0] = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   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 triangular,
+   //    quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
+   //    the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 3. Refine the mesh to increase the resolution. In this example we do
+   //    'ref_levels' of uniform refinement. We choose 'ref_levels' to be the
+   //    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);
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh->UniformRefinement();
+      }
+   }
+
+   // 4. Define a finite element space on the mesh. Here we use continuous
+   //    Lagrange finite elements of the specified order. If order < 1, we
+   //    instead use an isoparametric/isogeometric space.
+   FiniteElementCollection *fec;
+   NURBSExtension *NURBSext = NULL;
+   int own_fec = 0;
+
+   if (order[0] == -1) // Isoparametric
+   {
+      if (mesh->GetNodes())
+      {
+         fec = mesh->GetNodes()->OwnFEC();
+         own_fec = 0;
+         cout << "Using isoparametric FEs: " << fec->Name() << endl;
+      }
+      else
+      {
+         cout <<"Mesh does not have FEs --> Assume order 1.\n";
+         fec = new H1_FECollection(1, dim);
+         own_fec = 1;
+      }
+   }
+   else if (mesh->NURBSext && (order[0] > 0) )  // Subparametric NURBS
+   {
+      fec = new NURBSFECollection(order[0]);
+      own_fec = 1;
+      int nkv = mesh->NURBSext->GetNKV();
+
+      if (order.Size() == 1)
+      {
+         int tmp = order[0];
+         order.SetSize(nkv);
+         order = tmp;
+      }
+      if (order.Size() != nkv ) { mfem_error("Wrong number of orders set."); }
+      NURBSext = new NURBSExtension(mesh->NURBSext, order);
+   }
+   else
+   {
+      if (order.Size() > 1) { cout <<"Wrong number of orders set, needs one.\n"; }
+      fec = new H1_FECollection(abs(order[0]), dim);
+      own_fec = 1;
+   }
+   FiniteElementSpace *fespace = new FiniteElementSpace(mesh, NURBSext, fec);
+   cout << "Number of finite element unknowns: "
+        << fespace->GetTrueVSize() << endl;
+
+   // 5. Determine the list of true (i.e. conforming) essential boundary dofs.
+   //    In this example, the boundary conditions are defined by marking all
+   //    the boundary attributes from the mesh as essential (Dirichlet) and
+   //    converting them to a list of true dofs.
+   Array<int> ess_tdof_list;
+   if (mesh->bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(mesh->bdr_attributes.Max());
+      ess_bdr = 1;
+      fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   // 6. Set up the linear form b(.) which corresponds to the right-hand side of
+   //    the FEM linear system, which in this case is (1,phi_i) where phi_i are
+   //    the basis functions in the finite element fespace.
+   LinearForm *b = new LinearForm(fespace);
+   ConstantCoefficient one(1.0);
+   b->AddDomainIntegrator(new DomainLFIntegrator(one));
+   b->Assemble();
+
+   // 7. Define the solution vector x as a finite element grid function
+   //    corresponding to fespace. Initialize x with initial guess of zero,
+   //    which satisfies the boundary conditions.
+   GridFunction x(fespace);
+   x = 0.0;
+
+   // 8. Set up the bilinear form a(.,.) on the finite element space
+   //    corresponding to the Laplacian operator -Delta, by adding the Diffusion
+   //    domain integrator.
+   BilinearForm *a = new BilinearForm(fespace);
+   a->AddDomainIntegrator(new DiffusionIntegrator(one));
+
+   // 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;
+
+#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);
+#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);
+#endif
+
+   // 11. Recover the solution as a finite element grid function.
+   a->RecoverFEMSolution(X, *b, x);
+
+   // 12. Save the refined mesh and the solution. This output can be viewed 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);
+   x.Save(sol_ofs);
+
+   // 13. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << *mesh << x << flush;
+   }
+
+   // 14. Save data in the VisIt format
+   VisItDataCollection visit_dc("Example1", mesh);
+   visit_dc.RegisterField("solution", &x);
+   visit_dc.Save();
+
+   // 15. Free the used memory.
+   delete a;
+   delete b;
+   delete fespace;
+   if (own_fec) { delete fec; }
+   delete mesh;
+
+   return 0;
+}

miniapps/nurbs/ex11p.cpp

@@ -0,0 +1,397 @@
+//                    MFEM Example 11 - Parallel NURBS Version
+//
+// Compile with: make ex11p
+//
+// Sample runs:  mpirun -np 4 ex11p -m ../../data/square-disc.mesh
+//               mpirun -np 4 ex11p -m ../../data/star.mesh
+//               mpirun -np 4 ex11p -m ../../data/escher.mesh
+//               mpirun -np 4 ex11p -m ../../data/fichera.mesh
+//               mpirun -np 4 ex11p -m ../../data/square-disc-p2.vtk -o 2
+//               mpirun -np 4 ex11p -m ../../data/square-disc-p3.mesh -o 3
+//               mpirun -np 4 ex11p -m ../../data/square-disc-nurbs.mesh -o -1
+//               mpirun -np 4 ex11p -m ../../data/disc-nurbs.mesh -o -1 -n 20
+//               mpirun -np 4 ex11p -m ../../data/pipe-nurbs.mesh -o -1
+//               mpirun -np 4 ex11p -m ../../data/ball-nurbs.mesh -o 2
+//               mpirun -np 4 ex11p -m ../../data/star-surf.mesh
+//               mpirun -np 4 ex11p -m ../../data/square-disc-surf.mesh
+//               mpirun -np 4 ex11p -m ../../data/inline-segment.mesh
+//               mpirun -np 4 ex11p -m ../../data/amr-quad.mesh
+//               mpirun -np 4 ex11p -m ../../data/amr-hex.mesh
+//               mpirun -np 4 ex11p -m ../../data/mobius-strip.mesh -n 8
+//               mpirun -np 4 ex11p -m ../../data/klein-bottle.mesh -n 10
+//
+// Description:  This example code demonstrates the use of MFEM to solve the
+//               eigenvalue problem -Delta u = lambda u with homogeneous
+//               Dirichlet boundary conditions.
+//
+//               We compute a number of the lowest eigenmodes by discretizing
+//               the Laplacian and Mass operators using a FE space of the
+//               specified order, or an isoparametric/isogeometric space if
+//               order < 1 (quadratic for quadratic curvilinear mesh, NURBS for
+//               NURBS mesh, etc.)
+//
+//               The example highlights the use of the LOBPCG eigenvalue solver
+//               together with the BoomerAMG preconditioner in HYPRE, as well as
+//               optionally the SuperLU or STRUMPACK parallel direct solvers.
+//               Reusing a single GLVis visualization window for multiple
+//               eigenfunctions is also illustrated.
+//
+//               We recommend viewing Example 1 before viewing this example.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI.
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../../data/star.mesh";
+   int ser_ref_levels = 2;
+   int par_ref_levels = 1;
+   Array<int> order(1);
+   order[0] = 0;
+   int nev = 5;
+   int seed = 75;
+   bool slu_solver  = false;
+   bool sp_solver = false;
+   bool visualization = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
+                  "Number of times to refine the mesh uniformly in parallel.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&nev, "-n", "--num-eigs",
+                  "Number of desired eigenmodes.");
+   args.AddOption(&seed, "-s", "--seed",
+                  "Random seed used to initialize LOBPCG.");
+#ifdef MFEM_USE_SUPERLU
+   args.AddOption(&slu_solver, "-slu", "--superlu", "-no-slu",
+                  "--no-superlu", "Use the SuperLU Solver.");
+#endif
+#ifdef MFEM_USE_STRUMPACK
+   args.AddOption(&sp_solver, "-sp", "--strumpack", "-no-sp",
+                  "--no-strumpack", "Use the STRUMPACK Solver.");
+#endif
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (slu_solver && sp_solver)
+   {
+      if (myid == 0)
+         cout << "WARNING: Both SuperLU and STRUMPACK have been selected,"
+              << " please choose either one." << endl
+              << "         Defaulting to SuperLU." << endl;
+      sp_solver = false;
+   }
+   // The command line options are also passed to the STRUMPACK
+   // solver. So do not exit if some options are not recognized.
+   if (!sp_solver)
+   {
+      if (!args.Good())
+      {
+         if (myid == 0)
+         {
+            args.PrintUsage(cout);
+         }
+         MPI_Finalize();
+         return 1;
+      }
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // 3. Read the (serial) mesh from the given mesh file on all processors. We
+   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+   //    and volume meshes with the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 4. Refine the serial mesh on all processors to increase the resolution. In
+   //    this example we do 'ref_levels' of uniform refinement (2 by default, or
+   //    specified on the command line with -rs).
+   for (int lev = 0; lev < ser_ref_levels; lev++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution (1 time by
+   //    default, or specified on the command line with -rp). Once the parallel
+   //    mesh is defined, the serial mesh can be deleted.
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   for (int lev = 0; lev < par_ref_levels; lev++)
+   {
+      pmesh->UniformRefinement();
+   }
+
+   // 6. Define a parallel finite element space on the parallel mesh. Here we
+   //    use continuous Lagrange finite elements of the specified order. If
+   //    order < 1, we instead use an isoparametric/isogeometric space.
+   FiniteElementCollection *fec;
+   NURBSExtension *NURBSext = NULL;
+   int own_fec = 0;
+
+   if (order[0] == 0) // Isoparametric
+   {
+      if (pmesh->GetNodes())
+      {
+         fec = pmesh->GetNodes()->OwnFEC();
+         own_fec = 0;
+         cout << "Using isoparametric FEs: " << fec->Name() << endl;
+      }
+      else
+      {
+         cout <<"Mesh does not have FEs --> Assume order 1.\n";
+         fec = new H1_FECollection(1, dim);
+         own_fec = 1;
+      }
+   }
+   else if (pmesh->NURBSext && (order[0] > 0) )  // Subparametric NURBS
+   {
+      fec = new NURBSFECollection(order[0]);
+      own_fec = 1;
+      int nkv = pmesh->NURBSext->GetNKV();
+
+      if (order.Size() == 1)
+      {
+         int tmp = order[0];
+         order.SetSize(nkv);
+         order = tmp;
+      }
+      if (order.Size() != nkv ) { mfem_error("Wrong number of orders set."); }
+      NURBSext = new NURBSExtension(pmesh->NURBSext, order);
+   }
+   else
+   {
+      if (order.Size() > 1) { cout <<"Wrong number of orders set, needs one.\n"; }
+      fec = new H1_FECollection(abs(order[0]), dim);
+      own_fec = 1;
+   }
+   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh,NURBSext,fec);
+   HYPRE_Int size = fespace->GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of unknowns: " << size << endl;
+   }
+
+   // 7. Set up the parallel bilinear forms a(.,.) and m(.,.) on the finite
+   //    element space. The first corresponds to the Laplacian operator -Delta,
+   //    while the second is a simple mass matrix needed on the right hand side
+   //    of the generalized eigenvalue problem below. The boundary conditions
+   //    are implemented by elimination with special values on the diagonal to
+   //    shift the Dirichlet eigenvalues out of the computational range. After
+   //    serial and parallel assembly we extract the corresponding parallel
+   //    matrices A and M.
+   ConstantCoefficient one(1.0);
+   Array<int> ess_bdr;
+   if (pmesh->bdr_attributes.Size())
+   {
+      ess_bdr.SetSize(pmesh->bdr_attributes.Max());
+      ess_bdr = 1;
+   }
+
+   ParBilinearForm *a = new ParBilinearForm(fespace);
+   a->AddDomainIntegrator(new DiffusionIntegrator(one));
+   if (pmesh->bdr_attributes.Size() == 0)
+   {
+      // Add a mass term if the mesh has no boundary, e.g. periodic mesh or
+      // closed surface.
+      a->AddDomainIntegrator(new MassIntegrator(one));
+   }
+   a->Assemble();
+   a->EliminateEssentialBCDiag(ess_bdr, 1.0);
+   a->Finalize();
+
+   ParBilinearForm *m = new ParBilinearForm(fespace);
+   m->AddDomainIntegrator(new MassIntegrator(one));
+   m->Assemble();
+   // shift the eigenvalue corresponding to eliminated dofs to a large value
+   m->EliminateEssentialBCDiag(ess_bdr, numeric_limits<double>::min());
+   m->Finalize();
+
+   HypreParMatrix *A = a->ParallelAssemble();
+   HypreParMatrix *M = m->ParallelAssemble();
+
+#if defined(MFEM_USE_SUPERLU) || defined(MFEM_USE_STRUMPACK)
+   Operator * Arow = NULL;
+#ifdef MFEM_USE_SUPERLU
+   if (slu_solver)
+   {
+      Arow = new SuperLURowLocMatrix(*A);
+   }
+#endif
+#ifdef MFEM_USE_STRUMPACK
+   if (sp_solver)
+   {
+      Arow = new STRUMPACKRowLocMatrix(*A);
+   }
+#endif
+#endif
+
+   delete a;
+   delete m;
+
+   // 8. Define and configure the LOBPCG eigensolver and the BoomerAMG
+   //    preconditioner for A to be used within the solver. Set the matrices
+   //    which define the generalized eigenproblem A x = lambda M x.
+   Solver * precond = NULL;
+   if (!slu_solver && !sp_solver)
+   {
+      HypreBoomerAMG * amg = new HypreBoomerAMG(*A);
+      amg->SetPrintLevel(0);
+      precond = amg;
+   }
+   else
+   {
+#ifdef MFEM_USE_SUPERLU
+      if (slu_solver)
+      {
+         SuperLUSolver * superlu = new SuperLUSolver(MPI_COMM_WORLD);
+         superlu->SetPrintStatistics(false);
+         superlu->SetSymmetricPattern(true);
+         superlu->SetColumnPermutation(superlu::PARMETIS);
+         superlu->SetOperator(*Arow);
+         precond = superlu;
+      }
+#endif
+#ifdef MFEM_USE_STRUMPACK
+      if (sp_solver)
+      {
+         STRUMPACKSolver * strumpack = new STRUMPACKSolver(argc, argv, MPI_COMM_WORLD);
+         strumpack->SetPrintFactorStatistics(true);
+         strumpack->SetPrintSolveStatistics(false);
+         strumpack->SetKrylovSolver(strumpack::KrylovSolver::DIRECT);
+         strumpack->SetReorderingStrategy(strumpack::ReorderingStrategy::METIS);
+         strumpack->SetMC64Job(strumpack::MC64Job::NONE);
+         // strumpack->SetSymmetricPattern(true);
+         strumpack->SetOperator(*Arow);
+         strumpack->SetFromCommandLine();
+         precond = strumpack;
+      }
+#endif
+   }
+
+
+   HypreLOBPCG * lobpcg = new HypreLOBPCG(MPI_COMM_WORLD);
+   lobpcg->SetNumModes(nev);
+   lobpcg->SetRandomSeed(seed);
+   lobpcg->SetPreconditioner(*precond);
+   lobpcg->SetMaxIter(200);
+   lobpcg->SetTol(1e-8);
+   lobpcg->SetPrecondUsageMode(1);
+   lobpcg->SetPrintLevel(1);
+   lobpcg->SetMassMatrix(*M);
+   lobpcg->SetOperator(*A);
+
+   // 9. Compute the eigenmodes and extract the array of eigenvalues. Define a
+   //    parallel grid function to represent each of the eigenmodes returned by
+   //    the solver.
+   Array<double> eigenvalues;
+   lobpcg->Solve();
+   lobpcg->GetEigenvalues(eigenvalues);
+   ParGridFunction x(fespace);
+
+   // 10. Save the refined mesh and the modes in parallel. This output can be
+   //     viewed later using GLVis: "glvis -np <np> -m mesh -g mode".
+   {
+      ostringstream mesh_name, mode_name;
+      mesh_name << "mesh." << setfill('0') << setw(6) << myid;
+
+      ofstream mesh_ofs(mesh_name.str().c_str());
+      mesh_ofs.precision(8);
+      pmesh->Print(mesh_ofs);
+
+      for (int i=0; i<nev; i++)
+      {
+         // convert eigenvector from HypreParVector to ParGridFunction
+         x = lobpcg->GetEigenvector(i);
+
+         mode_name << "mode_" << setfill('0') << setw(2) << i << "."
+                   << setfill('0') << setw(6) << myid;
+
+         ofstream mode_ofs(mode_name.str().c_str());
+         mode_ofs.precision(8);
+         x.Save(mode_ofs);
+         mode_name.str("");
+      }
+   }
+
+   // 11. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream mode_sock(vishost, visport);
+      mode_sock.precision(8);
+
+      for (int i=0; i<nev; i++)
+      {
+         if ( myid == 0 )
+         {
+            cout << "Eigenmode " << i+1 << '/' << nev
+                 << ", Lambda = " << eigenvalues[i] << endl;
+         }
+
+         // convert eigenvector from HypreParVector to ParGridFunction
+         x = lobpcg->GetEigenvector(i);
+
+         mode_sock << "parallel " << num_procs << " " << myid << "\n"
+                   << "solution\n" << *pmesh << x << flush
+                   << "window_title 'Eigenmode " << i+1 << '/' << nev
+                   << ", Lambda = " << eigenvalues[i] << "'" << endl;
+
+         char c;
+         if (myid == 0)
+         {
+            cout << "press (q)uit or (c)ontinue --> " << flush;
+            cin >> c;
+         }
+         MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
+
+         if (c != 'c')
+         {
+            break;
+         }
+      }
+      mode_sock.close();
+   }
+
+   // 12. Free the used memory.
+   delete lobpcg;
+   delete precond;
+   delete M;
+   delete A;
+#if defined(MFEM_USE_SUPERLU) || defined(MFEM_USE_STRUMPACK)
+   delete Arow;
+#endif
+
+   delete fespace;
+   if (own_fec)
+   {
+      delete fec;
+   }
+   delete pmesh;
+
+   MPI_Finalize();
+
+   return 0;
+}

miniapps/nurbs/ex1p.cpp

@@ -0,0 +1,274 @@
+//                     MFEM Example 1 - Parallel NURBS Version
+//
+// Compile with: make ex1p
+//
+// Sample runs:  mpirun -np 4 ex1p -m ../../data/square-disc.mesh
+//               mpirun -np 4 ex1p -m ../../data/star.mesh
+//               mpirun -np 4 ex1p -m ../../data/escher.mesh
+//               mpirun -np 4 ex1p -m ../../data/fichera.mesh
+//               mpirun -np 4 ex1p -m ../../data/square-disc-p2.vtk -o 2
+//               mpirun -np 4 ex1p -m ../../data/square-disc-p3.mesh -o 3
+//               mpirun -np 4 ex1p -m ../../data/square-disc-nurbs.mesh -o -1
+//               mpirun -np 4 ex1p -m ../../data/disc-nurbs.mesh -o -1
+//               mpirun -np 4 ex1p -m ../../data/pipe-nurbs.mesh -o -1
+//               mpirun -np 4 ex1p -m ../../data/ball-nurbs.mesh -o 2
+//               mpirun -np 4 ex1p -m ../../data/star-surf.mesh
+//               mpirun -np 4 ex1p -m ../../data/square-disc-surf.mesh
+//               mpirun -np 4 ex1p -m ../../data/inline-segment.mesh
+//               mpirun -np 4 ex1p -m ../../data/amr-quad.mesh
+//               mpirun -np 4 ex1p -m ../../data/amr-hex.mesh
+//               mpirun -np 4 ex1p -m ../../data/mobius-strip.mesh
+//               mpirun -np 4 ex1p -m ../../data/mobius-strip.mesh -o -1 -sc
+//
+// Description:  This example code demonstrates the use of MFEM to define a
+//               simple finite element discretization of the Laplace problem
+//               -Delta u = 1 with homogeneous Dirichlet boundary conditions.
+//               Specifically, we discretize using a FE space of the specified
+//               order, or if order < 1 using an isoparametric/isogeometric
+//               space (i.e. quadratic for quadratic curvilinear mesh, NURBS for
+//               NURBS mesh, etc.)
+//
+//               The example highlights the use of mesh refinement, finite
+//               element grid functions, as well as linear and bilinear forms
+//               corresponding to the left-hand side and right-hand side of the
+//               discrete linear system. We also cover the explicit elimination
+//               of essential boundary conditions, static condensation, and the
+//               optional connection to the GLVis tool for visualization.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI.
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../../data/star.mesh";
+   Array<int> order(1);
+   order[0] = 1;
+   bool static_cond = false;
+   bool visualization = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // 3. Read the (serial) mesh from the given mesh file on all processors.  We
+   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+   //    and volume meshes with the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 4. Refine the serial mesh on all processors to increase the resolution. In
+   //    this example we do 'ref_levels' of uniform refinement. We choose
+   //    '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);
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh->UniformRefinement();
+      }
+   }
+
+   // 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.
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   {
+      int par_ref_levels = 2;
+      for (int l = 0; l < par_ref_levels; l++)
+      {
+         pmesh->UniformRefinement();
+      }
+   }
+
+   // 6. Define a parallel finite element space on the parallel mesh. Here we
+   //    use continuous Lagrange finite elements of the specified order. If
+   //    order < 1, we instead use an isoparametric/isogeometric space.
+   FiniteElementCollection *fec;
+   NURBSExtension *NURBSext = NULL;
+   int own_fec = 0;
+
+   if (order[0] == -1) // Isoparametric
+   {
+      if (pmesh->GetNodes())
+      {
+         fec = pmesh->GetNodes()->OwnFEC();
+         own_fec = 0;
+         cout << "Using isoparametric FEs: " << fec->Name() << endl;
+      }
+      else
+      {
+         cout <<"Mesh does not have FEs --> Assume order 1.\n";
+         fec = new H1_FECollection(1, dim);
+         own_fec = 1;
+      }
+   }
+   else if (pmesh->NURBSext && (order[0] > 0) )  // Subparametric NURBS
+   {
+      fec = new NURBSFECollection(order[0]);
+      own_fec = 1;
+      int nkv = pmesh->NURBSext->GetNKV();
+
+      if (order.Size() == 1)
+      {
+         int tmp = order[0];
+         order.SetSize(nkv);
+         order = tmp;
+      }
+      if (order.Size() != nkv ) { mfem_error("Wrong number of orders set."); }
+      NURBSext = new NURBSExtension(pmesh->NURBSext, order);
+   }
+   else
+   {
+      if (order.Size() > 1) { cout <<"Wrong number of orders set, needs one.\n"; }
+      fec = new H1_FECollection(abs(order[0]), dim);
+      own_fec = 1;
+   }
+   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh,NURBSext,fec);
+   HYPRE_Int size = fespace->GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of finite element unknowns: " << size << endl;
+   }
+
+   // 7. Determine the list of true (i.e. parallel conforming) essential
+   //    boundary dofs. In this example, the boundary conditions are defined
+   //    by marking all the boundary attributes from the mesh as essential
+   //    (Dirichlet) and converting them to a list of true dofs.
+   Array<int> ess_tdof_list;
+   if (pmesh->bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+      ess_bdr = 1;
+      fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   // 8. Set up the parallel linear form b(.) which corresponds to the
+   //    right-hand side of the FEM linear system, which in this case is
+   //    (1,phi_i) where phi_i are the basis functions in fespace.
+   ParLinearForm *b = new ParLinearForm(fespace);
+   ConstantCoefficient one(1.0);
+   b->AddDomainIntegrator(new DomainLFIntegrator(one));
+   b->Assemble();
+
+   // 9. Define the solution vector x as a parallel finite element grid function
+   //    corresponding to fespace. Initialize x with initial guess of zero,
+   //    which satisfies the boundary conditions.
+   ParGridFunction x(fespace);
+   x = 0.0;
+
+   // 10. Set up the parallel bilinear form a(.,.) on the finite element space
+   //     corresponding to the Laplacian operator -Delta, by adding the Diffusion
+   //     domain integrator.
+   ParBilinearForm *a = new ParBilinearForm(fespace);
+   a->AddDomainIntegrator(new DiffusionIntegrator(one));
+
+   // 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();
+
+   HypreParMatrix A;
+   Vector B, X;
+   a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
+
+   if (myid == 0)
+   {
+      cout << "Size of linear system: " << A.GetGlobalNumRows() << 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);
+   pcg->Mult(B, X);
+
+   // 13. Recover the parallel grid function corresponding to X. This is the
+   //     local finite element solution on each processor.
+   a->RecoverFEMSolution(X, *b, x);
+
+   // 14. Save the refined mesh and the solution in parallel. This output can
+   //     be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
+   {
+      ostringstream mesh_name, sol_name;
+      mesh_name << "mesh." << setfill('0') << setw(6) << myid;
+      sol_name << "sol." << setfill('0') << setw(6) << myid;
+
+      ofstream mesh_ofs(mesh_name.str().c_str());
+      mesh_ofs.precision(8);
+      pmesh->Print(mesh_ofs);
+
+      ofstream sol_ofs(sol_name.str().c_str());
+      sol_ofs.precision(8);
+      x.Save(sol_ofs);
+   }
+
+   // 15. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << *pmesh << x << flush;
+   }
+
+   // 16. Save data in the VisIt format
+   VisItDataCollection visit_dc("Example1-Parallel", pmesh);
+   visit_dc.RegisterField("solution", &x);
+   visit_dc.Save();
+
+   // 17. Free the used memory.
+   delete pcg;
+   delete amg;
+   delete a;
+   delete b;
+   delete fespace;
+   if (own_fec) { delete fec; }
+   delete pmesh;
+
+   MPI_Finalize();
+
+   return 0;
+}

miniapps/nurbs/makefile

@@ -0,0 +1,70 @@
+# 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.
+
+# Use the MFEM build directory
+MFEM_DIR ?= ../..
+MFEM_BUILD_DIR ?= ../..
+SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/nurbs/,)
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+# Use the MFEM install directory
+# MFEM_INSTALL_DIR = ../../mfem
+# CONFIG_MK = $(MFEM_INSTALL_DIR)/share/mfem/config.mk
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+SEQ_MINIAPPS = ex1
+PAR_MINIAPPS = ex1p ex11p
+ifeq ($(MFEM_USE_MPI),NO)
+   MINIAPPS = $(SEQ_MINIAPPS)
+else
+   MINIAPPS = $(PAR_MINIAPPS) $(SEQ_MINIAPPS)
+endif
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all clean clean-build clean-exec
+.PRECIOUS: %.o
+
+# Remove built-in rule
+%: %.cpp
+
+# Replace the default implicit rule for *.cpp files
+%: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(MFEM_LIBS)
+
+all: $(MINIAPPS)
+
+MFEM_TESTS = MINIAPPS
+include $(MFEM_TEST_MK)
+
+# Testing: Parallel vs. serial runs
+RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
+%-test-par: %
+	@$(call mfem-test,$<, $(RUN_MPI), NURBS miniapp)
+%-test-seq: %
+	@$(call mfem-test,$<,, NURBS miniapp)
+
+# Testing: "test" target and mfem-test* variables are defined in config/test.mk
+
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
+
+clean: clean-build clean-exec
+
+clean-build:
+	rm -f *.o *~ $(SEQ_MINIAPPS) $(PAR_MINIAPPS)
+	rm -rf *.dSYM *.TVD.*breakpoints
+
+clean-exec:
+	@rm -f refined.mesh mesh.* sol.* mode_*
+	@rm -rf Example1*

miniapps/nurbs/pipe-nurbs-2d.mesh

@@ -0,0 +1,70 @@
+MFEM NURBS mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# SEGMENT     = 1
+# SQUARE      = 3
+# CUBE        = 5
+#
+
+dimension
+2
+
+elements
+1
+1 3 0 1 2 3
+
+boundary
+4
+1 1 0 1
+1 1 2 3
+1 1 3 0
+1 1 1 2
+
+edges
+4
+0 0 1
+0 3 2
+1 0 3
+1 1 2
+
+vertices
+4
+
+knotvectors
+2
+2 3 0 0 0 1 1 1
+2 3 0 0 0 1 1 1
+
+weights
+1
+1
+1
+1
+
+0.7071067811865475244
+0.7071067811865475244
+
+1
+1
+
+0.7071067811865475244
+
+FiniteElementSpace
+FiniteElementCollection: NURBS2
+VDim: 2
+Ordering: 1
+
+0 0
+2 2
+1 2
+0 1
+
+2 0
+1 1
+
+0 0.5
+1.5 2
+
+1.5 0.5

miniapps/nurbs/square-nurbs.mesh

@@ -0,0 +1,57 @@
+MFEM NURBS mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# SEGMENT     = 1
+# SQUARE      = 3
+# CUBE        = 5
+#
+
+dimension
+2
+
+elements
+1
+1 3 0 1 2 3
+
+boundary
+4
+1 1 0 1
+1 1 2 3
+1 1 3 0
+1 1 1 2
+
+
+
+edges
+4
+0 0 1
+0 3 2
+1 0 3
+1 1 2
+
+vertices
+4
+
+knotvectors
+2
+1 2 0 0 1 1
+1 2 0 0 1 1
+
+weights
+1
+1
+1
+1
+
+FiniteElementSpace
+FiniteElementCollection: NURBS1
+VDim: 2
+Ordering: 1
+
+0 0
+1 0
+1 1
+0 1
+