Merge branch 'master' into najlkin/misc-impro

fem/bilininteg.hpp

@@ -3362,7 +3362,7 @@ using ConservativeDGTraceIntegrator = DGTraceIntegrator;
 /** Integrator that represents the face terms used for the non-conservative
     DG discretization of the convection equation:
     $$
-      -\alpha \langle \rho_u (u \cdot n) \{v\},[w] \rangle + \beta \langle \rho_u |u \cdot n| [v],[w] \rangle.
+      -\alpha \langle \rho_u [v], (u \cdot n) \{w\} \rangle + \beta \langle \rho_u [v], |u \cdot n| [w] \rangle.
     $$
     This integrator can be used with together with ConvectionIntegrator to
     implement an upwind DG discretization in non-conservative form, see ex9 and

fem/coefficient.hpp

@@ -261,6 +261,9 @@ public:
    /// Evaluate the coefficient at @a ip.
    real_t Eval(ElementTransformation &T,
                const IntegrationPoint &ip) override;
+
+   const std::function<real_t(const Vector &)>& GetFunction() const { return Function; }
+   const std::function<real_t(const Vector &,real_t)>& GetTDFunction() const { return TDFunction; }
 };
 
 /// A common base class for returning individual components of the domain's
@@ -773,6 +776,10 @@ public:
    void Eval(Vector &V, ElementTransformation &T,
              const IntegrationPoint &ip) override;
 
+   const std::function<void(const Vector &, Vector &)>& GetFunction() const { return Function; }
+   const std::function<void(const Vector &, real_t, Vector &)>& GetTDFunction()
+   const { return TDFunction; }
+
    virtual ~VectorFunctionCoefficient() { }
 };
 
@@ -995,7 +1002,10 @@ public:
    DeltaCoefficient& GetDeltaCoefficient() { return d; }
 
    void SetScale(real_t s) { d.SetScale(s); }
+   void SetTol(real_t tol) { d.SetTol(tol); }
+
    void SetDirection(const Vector& d_);
+   void GetDirection(Vector &d_) { d_ = dir; }
 
    void SetDeltaCenter(const Vector& center) { d.SetDeltaCenter(center); }
    void GetDeltaCenter(Vector& center) { d.GetDeltaCenter(center); }

fem/gridfunc.cpp

@@ -2358,6 +2358,87 @@ void GridFunction::ProjectDeltaCoefficient(DeltaCoefficient &delta_coeff,
    }
 }
 
+void GridFunction::ProjectVectorDeltaCoefficient(VectorDeltaCoefficient
+                                                 &vdelta_coeff,
+                                                 real_t &integral)
+{
+   DeltaCoefficient &delta_coeff = vdelta_coeff.GetDeltaCoefficient();
+
+   if (!fes->GetNE())
+   {
+      integral = 0.0;
+      return;
+   }
+
+   Mesh *mesh = fes->GetMesh();
+   const int dim = mesh->Dimension();
+   const real_t *center = delta_coeff.Center();
+   const real_t *vert = mesh->GetVertex(0);
+   real_t min_dist, dist;
+   int v_idx = 0;
+
+   // find the vertex closest to the center of the delta function
+   min_dist = Distance(center, vert, dim);
+   for (int i = 0; i < mesh->GetNV(); i++)
+   {
+      vert = mesh->GetVertex(i);
+      dist = Distance(center, vert, dim);
+      if (dist < min_dist)
+      {
+         min_dist = dist;
+         v_idx = i;
+      }
+   }
+
+   (*this) = 0.0;
+   integral = 0.0;
+
+   if (min_dist >= delta_coeff.Tol())
+   {
+      return;
+   }
+
+   // find the elements that have 'v_idx' as a vertex
+   Vector dir;
+   vdelta_coeff.GetDirection(dir);
+   MFEM_ASSERT(fes->GetVDim() == dir.Size(),
+               "Vector dimension of the grid function does not match the coefficient.");
+
+   MassIntegrator Mi(*delta_coeff.Weight());
+   DenseMatrix loc_mass;
+   Array<int> dofs, vdofs, vertices;
+   Vector vals, vvals, loc_mass_vals;
+   for (int i = 0; i < mesh->GetNE(); i++)
+   {
+      mesh->GetElementVertices(i, vertices);
+      for (int j = 0; j < vertices.Size(); j++)
+         if (vertices[j] == v_idx)
+         {
+            const FiniteElement *fe = fes->GetFE(i);
+            MFEM_ASSERT(fe->GetRangeType() == FiniteElement::SCALAR,
+                        "Implemented only for scalar finite elements.");
+
+            Mi.AssembleElementMatrix(*fe, *fes->GetElementTransformation(i),
+                                     loc_mass);
+            vals.SetSize(fe->GetDof());
+            fe->ProjectDelta(j, vals);
+            fes->GetElementDofs(i, dofs);
+            vvals.SetSize(vals.Size());
+            for (int d = 0; d < dir.Size(); d++)
+            {
+               vdofs = dofs;
+               fes->DofsToVDofs(d, vdofs);
+               vvals.Set(dir[d], vals);
+               SetSubVector(vdofs, vvals);
+            }
+            loc_mass_vals.SetSize(vals.Size());
+            loc_mass.Mult(vals, loc_mass_vals);
+            integral += loc_mass_vals.Sum(); // partition of unity basis
+            break;
+         }
+   }
+}
+
 void GridFunction::ProjectCoefficient(Coefficient &coeff)
 {
    DeltaCoefficient *delta_c = dynamic_cast<DeltaCoefficient *>(&coeff);
@@ -2415,7 +2496,10 @@ void GridFunction::ProjectCoefficient(Coefficient &coeff)
 
       ProjectDeltaCoefficient(*delta_c, integral);
 
-      (*this) *= (delta_c->Scale() / integral);
+      if (integral > 0.)
+      {
+         (*this) *= (delta_c->Scale() / integral);
+      }
    }
 }
 
@@ -2445,51 +2529,66 @@ void GridFunction::ProjectCoefficient(
 
 void GridFunction::ProjectCoefficient(VectorCoefficient &vcoeff)
 {
-   if (fes->GetNURBSext() == NULL)
+   VectorDeltaCoefficient *vdelta_c =
+      dynamic_cast<VectorDeltaCoefficient *>(&vcoeff);
+
+   if (vdelta_c == NULL)
    {
-      int i;
-      Array<int> vdofs;
-      Vector vals;
-
-      DofTransformation * doftrans = NULL;
-
-      for (i = 0; i < fes->GetNE(); i++)
+      if (fes->GetNURBSext() == NULL)
       {
-         doftrans = fes->GetElementVDofs(i, vdofs);
-         vals.SetSize(vdofs.Size());
-         fes->GetFE(i)->Project(vcoeff, *fes->GetElementTransformation(i), vals);
-         if (doftrans)
+         Array<int> vdofs;
+         Vector vals;
+         DofTransformation * doftrans = NULL;
+
+         for (int i = 0; i < fes->GetNE(); i++)
          {
-            doftrans->TransformPrimal(vals);
+            doftrans = fes->GetElementVDofs(i, vdofs);
+            vals.SetSize(vdofs.Size());
+            fes->GetFE(i)->Project(vcoeff, *fes->GetElementTransformation(i), vals);
+            if (doftrans)
+            {
+               doftrans->TransformPrimal(vals);
+            }
+            SetSubVector(vdofs, vals);
          }
-         SetSubVector(vdofs, vals);
+      }
+      else
+      {
+         // Define and assemble linear form
+         LinearForm b(fes);
+         b.AddDomainIntegrator(new VectorFEDomainLFIntegrator(vcoeff));
+         b.Assemble();
+
+         // Define and assemble bilinear form
+         BilinearForm a(fes);
+         a.AddDomainIntegrator(new VectorFEMassIntegrator());
+         a.Assemble();
+
+         // Set solver and preconditioner
+         SparseMatrix A(a.SpMat());
+         GSSmoother  prec(A);
+         CGSolver cg;
+         cg.SetOperator(A);
+         cg.SetPreconditioner(prec);
+         cg.SetRelTol(1e-12);
+         cg.SetMaxIter(1000);
+         cg.SetPrintLevel(0);
+
+         // Solve and get solution
+         *this = 0.0;
+         cg.Mult(b,*this);
       }
    }
    else
    {
-      // Define and assemble linear form
-      LinearForm b(fes);
-      b.AddDomainIntegrator(new VectorFEDomainLFIntegrator(vcoeff));
-      b.Assemble();
+      real_t integral;
 
-      // Define and assemble bilinear form
-      BilinearForm a(fes);
-      a.AddDomainIntegrator(new VectorFEMassIntegrator());
-      a.Assemble();
+      ProjectVectorDeltaCoefficient(*vdelta_c, integral);
 
-      // Set solver and preconditioner
-      SparseMatrix A(a.SpMat());
-      GSSmoother  prec(A);
-      CGSolver cg;
-      cg.SetOperator(A);
-      cg.SetPreconditioner(prec);
-      cg.SetRelTol(1e-12);
-      cg.SetMaxIter(1000);
-      cg.SetPrintLevel(0);
-
-      // Solve and get solution
-      *this = 0.0;
-      cg.Mult(b,*this);
+      if (integral > 0.)
+      {
+         (*this) *= (vdelta_c->GetDeltaCoefficient().Scale() / integral);
+      }
    }
 }
 

fem/gridfunc.hpp

@@ -53,6 +53,11 @@ protected:
    void ProjectDeltaCoefficient(DeltaCoefficient &delta_coeff,
                                 real_t &integral);
 
+   // Project the vector delta coefficient without scaling and return the
+   // (local) integral of the projection.
+   void ProjectVectorDeltaCoefficient(VectorDeltaCoefficient &vdelta_coeff,
+                                      real_t &integral);
+
    // Sum fluxes to vertices and count element contributions
    void SumFluxAndCount(BilinearFormIntegrator &blfi,
                         GridFunction &flux,

fem/linearform.hpp

@@ -179,16 +179,29 @@ public:
    /// Access all integrators added with AddBoundaryIntegrator().
    Array<LinearFormIntegrator*> *GetBLFI() { return &boundary_integs; }
 
+   /// Access all boundary markers added with AddBoundaryIntegrator().
+   /** If no marker was specified when the integrator was added, the
+       corresponding pointer (to Array<int>) will be NULL. */
+   Array<Array<int>*> *GetBLFI_Marker() { return &boundary_integs_marker; }
+
    /// Access all integrators added with AddBdrFaceIntegrator().
-   Array<LinearFormIntegrator*> *GetFLFI() { return &boundary_face_integs; }
+   Array<LinearFormIntegrator*> *GetBFLFI() { return &boundary_face_integs; }
+
+   /// (DEPRECATED) Access all integrators added with AddBdrFaceIntegrator().
+   /** @deprecated Use GetBFLFI() instead. */
+   MFEM_DEPRECATED Array<LinearFormIntegrator*> *GetFLFI() { return &boundary_face_integs; }
 
    /// Access all integrators added with AddInteriorFaceIntegrator().
    Array<LinearFormIntegrator*> *GetIFLFI() { return &interior_face_integs; }
 
-   /** @brief Access all boundary markers added with AddBdrFaceIntegrator().
-       If no marker was specified when the integrator was added, the
+   /// Access all boundary markers added with AddBdrFaceIntegrator().
+   /** If no marker was specified when the integrator was added, the
        corresponding pointer (to Array<int>) will be NULL. */
-   Array<Array<int>*> *GetFLFI_Marker() { return &boundary_face_integs_marker; }
+   Array<Array<int>*> *GetBFLFI_Marker() { return &boundary_face_integs_marker; }
+
+   //// (DEPRECATED) Access all boundary markers added with AddBdrFaceIntegrator().
+   /** @deprecated Use GetBFLFI_Marker() instead. */
+   MFEM_DEPRECATED Array<Array<int>*> *GetFLFI_Marker() { return &boundary_face_integs_marker; }
 
    /// @brief Which assembly algorithm to use: the new device-compatible fast
    /// assembly (true), or the legacy CPU-only algorithm (false).

fem/lor/lor.cpp

@@ -87,7 +87,7 @@ int LORBase::GetLOROrder() const
    return (type == L2 || type == RT) ? 0 : 1;
 }
 
-void LORBase::ConstructLocalDofPermutation(Array<int> &perm_) const
+void LORBase::ConstructLocalDofPermutation() const
 {
    FESpaceType type = GetFESpaceType();
    MFEM_VERIFY(type != H1 && type != L2, "");
@@ -107,7 +107,7 @@ void LORBase::ConstructLocalDofPermutation(Array<int> &perm_) const
 
    using GeomRef = std::pair<Geometry::Type, int>;
    std::map<GeomRef, int> point_matrices_offsets;
-   perm_.SetSize(fes_lor.GetVSize());
+   l_perm.SetSize(fes_lor.GetVSize());
 
    Array<int> vdof_ho, vdof_lor;
    for (int ilor=0; ilor<mesh_lor.GetNE(); ++ilor)
@@ -134,7 +134,7 @@ void LORBase::ConstructLocalDofPermutation(Array<int> &perm_) const
 
       if (type == L2)
       {
-         perm_[vdof_lor[0]] = vdof_ho[lor_index];
+         l_perm[vdof_lor[0]] = vdof_ho[lor_index];
          continue;
       }
 
@@ -166,7 +166,7 @@ void LORBase::ConstructLocalDofPermutation(Array<int> &perm_) const
                int s4 = idof_ho < 0 ? -1 : 1;
                int s = s1*s2*s3*s4;
                i = absdof(idof_ho);
-               perm_[absdof(idof_lor)] = s < 0 ? -1-absdof(i) : absdof(i);
+               l_perm[absdof(idof_lor)] = s < 0 ? -1-absdof(i) : absdof(i);
             }
          }
       };
@@ -217,6 +217,8 @@ void LORBase::ConstructDofPermutation() const
       return;
    }
 
+   ConstructLocalDofPermutation();
+
 #ifdef MFEM_USE_MPI
    ParFiniteElementSpace *pfes_ho
       = dynamic_cast<ParFiniteElementSpace*>(&fes_ho);
@@ -224,8 +226,6 @@ void LORBase::ConstructDofPermutation() const
       = dynamic_cast<ParFiniteElementSpace*>(&GetFESpace());
    if (pfes_ho && pfes_lor)
    {
-      Array<int> l_perm;
-      ConstructLocalDofPermutation(l_perm);
       perm.SetSize(pfes_lor->GetTrueVSize());
       for (int i=0; i<l_perm.Size(); ++i)
       {
@@ -245,7 +245,7 @@ void LORBase::ConstructDofPermutation() const
    else
 #endif
    {
-      ConstructLocalDofPermutation(perm);
+      perm = l_perm;
    }
 }
 
@@ -255,6 +255,12 @@ const Array<int> &LORBase::GetDofPermutation() const
    return perm;
 }
 
+const Array<int> &LORBase::GetLocalDofPermutation() const
+{
+   if (l_perm.Size() == 0) { ConstructLocalDofPermutation(); }
+   return l_perm;
+}
+
 bool LORBase::HasSameDofNumbering() const
 {
    FESpaceType type = GetFESpaceType();
@@ -277,8 +283,7 @@ void LORBase::SetupProlongationAndRestriction()
 {
    if (!HasSameDofNumbering())
    {
-      Array<int> p;
-      ConstructLocalDofPermutation(p);
+      const Array<int> &p = GetLocalDofPermutation();
       fes->CopyProlongationAndRestriction(fes_ho, &p);
    }
    else

fem/lor/lor.hpp

@@ -70,12 +70,12 @@ protected:
    BilinearForm *a = nullptr;
    class BatchedLORAssembly *batched_lor = nullptr;
    OperatorHandle A;
-   mutable Array<int> perm;
+   mutable Array<int> perm, l_perm;
 
    /// Constructs the local DOF (ldof) permutation. In parallel this is used as
    /// an intermediate step in computing the DOF permutation (see
    /// ConstructDofPermutation and GetDofPermutation).
-   void ConstructLocalDofPermutation(Array<int> &perm_) const;
+   void ConstructLocalDofPermutation() const;
 
    /// Construct the permutation that maps LOR DOFs to high-order DOFs. See
    /// GetDofPermutation.
@@ -127,6 +127,16 @@ public:
    /// LOR dof, @a perm[i] is the index of the corresponding HO dof.
    const Array<int> &GetDofPermutation() const;
 
+   /// @brief Returns the local permutation that maps LOR DOFs to high-order DOFs.
+   ///
+   /// This permutation is constructed the first time it is requested, and then
+   /// is cached. When running in parallel, only permutation of the local DOFs
+   /// is returned. These can be inconsistent with the orientation of the global
+   /// DOFs and GetDofPermutation() should be used instead to get a consistent
+   /// global permutation.
+   /// @see GetDofPermutation()
+   const Array<int> &GetLocalDofPermutation() const;
+
    /// Returns the low-order refined finite element space.
    FiniteElementSpace &GetFESpace() const;
 

general/optparser.cpp

@@ -48,6 +48,35 @@ int isValidAsInt(char * s)
    return 1;
 }
 
+int isValidAsUnsigned(char * s)
+{
+   if ( s == NULL || *s == '\0' )
+   {
+      return 0;   //Empty string
+   }
+
+   if ( *s == '+' )
+   {
+      ++s;
+   }
+
+   if ( *s == '\0')
+   {
+      return 0;   //sign character only
+   }
+
+   while (*s)
+   {
+      if ( !isdigit(*s) )
+      {
+         return 0;
+      }
+      ++s;
+   }
+
+   return 1;
+}
+
 int isValidAsDouble(char * s)
 {
    // A valid floating point number for atof using the "C" locale is formed by
@@ -200,6 +229,10 @@ void OptionsParser::Parse()
                   isValid = isValidAsInt(argv[i]);
                   *(int *)(options[j].var_ptr) = atoi(argv[i++]);
                   break;
+               case UNSIGNED:
+                  isValid = isValidAsUnsigned(argv[i]);
+                  *(unsigned *)(options[j].var_ptr) = strtoul(argv[i++], NULL, 0);
+                  break;
                case DOUBLE:
                   isValid = isValidAsDouble(argv[i]);
                   *(real_t *)(options[j].var_ptr) = atof(argv[i++]);
@@ -279,6 +312,10 @@ void OptionsParser::WriteValue(const Option &opt, std::ostream &os)
          os << *(int *)(opt.var_ptr);
          break;
 
+      case UNSIGNED:
+         os << *(unsigned *)(opt.var_ptr);
+         break;
+
       case DOUBLE:
          os << *(real_t *)(opt.var_ptr);
          break;
@@ -288,7 +325,7 @@ void OptionsParser::WriteValue(const Option &opt, std::ostream &os)
          break;
 
       case STD_STRING:
-         out << *(std::string *)(opt.var_ptr);
+         os << *(std::string *)(opt.var_ptr);
          break;
 
       case ARRAY:
@@ -408,9 +445,9 @@ void OptionsParser::PrintHelp(ostream &os) const
    static const char *seprtr = ", ";
    static const char *descr_sep = "\n\t";
    static const char *line_sep = "";
-   static const char *types[] = { " <int>", " <double>", " <string>",
-                                  " <string>", "", "", " '<int>...'",
-                                  " '<double>...'"
+   static const char *types[] = { " <int>", " <unsigned>", " <double>",
+                                  " <string>", " <string>", "", "",
+                                  " '<int>...'", " '<double>...'"
                                 };
 
    os << indent << "-h" << seprtr << "--help" << descr_sep

general/optparser.hpp

@@ -31,7 +31,7 @@ class Vector;
 class OptionsParser
 {
 public:
-   enum OptionType { INT, DOUBLE, STRING, STD_STRING, ENABLE, DISABLE, ARRAY, VECTOR };
+   enum OptionType { INT, UNSIGNED, DOUBLE, STRING, STD_STRING, ENABLE, DISABLE, ARRAY, VECTOR };
 
 private:
    struct Option
@@ -98,6 +98,14 @@ public:
                             required));
    }
 
+   /// Add an unsigned integer option and set 'var' to receive the value.
+   void AddOption(unsigned *var, const char *short_name, const char *long_name,
+                  const char *description, bool required = false)
+   {
+      options.Append(Option(UNSIGNED, var, short_name, long_name, description,
+                            required));
+   }
+
    /// Add a double option and set 'var' to receive the value.
    void AddOption(real_t *var, const char *short_name, const char *long_name,
                   const char *description, bool required = false)

linalg/densemat.cpp

@@ -1432,7 +1432,7 @@ void DenseMatrix::Diag(real_t c, int n)
    }
 }
 
-void DenseMatrix::Diag(real_t *diag, int n)
+void DenseMatrix::Diag(const real_t *diag, int n)
 {
    SetSize(n);
 
@@ -4398,6 +4398,57 @@ DenseTensor &DenseTensor::operator=(const DenseTensor &other)
    return *this;
 }
 
+DenseTensor &DenseTensor::operator+=(const real_t *m)
+{
+   int s = SizeI() * SizeJ() * SizeK();
+   for (int i = 0; i < s; i++)
+   {
+      tdata[i] += m[i];
+   }
+   return *this;
+}
+
+DenseTensor &DenseTensor::operator+=(const DenseTensor &m)
+{
+   MFEM_ASSERT(SizeI() == m.SizeI() && SizeJ() == m.SizeJ() &&
+               SizeK() == m.SizeK(),
+               "incompatible tensor sizes.");
+   return *this += m.Data();
+}
+
+DenseTensor &DenseTensor::operator-=(const DenseTensor &m)
+{
+   MFEM_ASSERT(SizeI() == m.SizeI() && SizeJ() == m.SizeJ() &&
+               SizeK() == m.SizeK(),
+               "incompatible tensor sizes.");
+
+   int s = SizeI() * SizeJ() * SizeK();
+   for (int i = 0; i < s; i++)
+   {
+      tdata[i] -= m.tdata[i];
+   }
+   return *this;
+}
+
+DenseTensor &DenseTensor::operator*=(real_t c)
+{
+   int s = SizeI() * SizeJ() * SizeK();
+   for (int i = 0; i < s; i++)
+   {
+      tdata[i] *= c;
+   }
+   return *this;
+}
+
+void DenseTensor::Neg()
+{
+   int s = SizeI() * SizeJ() * SizeK();
+   for (int i = 0; i < s; i++)
+   {
+      tdata[i] = -tdata[i];
+   }
+}
+
 void BatchLUFactor(DenseTensor &Mlu, Array<int> &P, const real_t TOL)
 {
    BatchedLinAlg::LUFactor(Mlu, P);

linalg/densemat.hpp

@@ -341,7 +341,9 @@ public:
    /// Creates n x n diagonal matrix with diagonal elements c
    void Diag(real_t c, int n);
    /// Creates n x n diagonal matrix with diagonal given by diag
-   void Diag(real_t *diag, int n);
+   void Diag(const real_t *diag, int n);
+   /// Creates a diagonal matrix with diagonal given by diag
+   void Diag(const Vector &diag) { Diag(diag.GetData(), diag.Size()); }
 
    /// (*this) = (*this)^t
    void Transpose();
@@ -1208,6 +1210,16 @@ public:
    /// Copy assignment operator (performs a deep copy)
    DenseTensor &operator=(const DenseTensor &other);
 
+   DenseTensor &operator+=(const real_t *m);
+   DenseTensor &operator+=(const DenseTensor &m);
+
+   DenseTensor &operator-=(const DenseTensor &m);
+
+   DenseTensor &operator*=(real_t c);
+
+   /// (*this) = -(*this)
+   void Neg();
+
    DenseMatrix &operator()(int k)
    {
       MFEM_ASSERT_INDEX_IN_RANGE(k, 0, SizeK());

linalg/sparsemat.cpp

@@ -3234,6 +3234,36 @@ SparseMatrix &SparseMatrix::operator+=(const SparseMatrix &B)
    return (*this);
 }
 
+void SparseMatrix::AddSubSparseMatrix(const real_t a, const SparseMatrix &B,
+                                      int ioff, int joff)
+{
+   MFEM_ASSERT(height >= ioff + B.height && width >= joff + B.width,
+               "Mismatch of this matrix size and the sub-matrix.  This height = "
+               << height << ", width = " << width << ", B.height = "
+               << B.height << ", B.width = " << B.width
+               << ", row offset = " << ioff << ", column offset = " << joff);
+
+   for (int i = 0; i < B.height; i++)
+   {
+      SetColPtr(i+ioff);
+      if (B.Rows)
+      {
+         for (RowNode *aux = B.Rows[i]; aux != NULL; aux = aux->Prev)
+         {
+            _Add_(aux->Column+joff, a*aux->Value);
+         }
+      }
+      else
+      {
+         for (int j = B.I[i]; j < B.I[i+1]; j++)
+         {
+            _Add_(B.J[j]+joff, a*B.A[j]);
+         }
+      }
+      ClearColPtr();
+   }
+}
+
 void SparseMatrix::Add(const real_t a, const SparseMatrix &B)
 {
    for (int i = 0; i < height; i++)

linalg/sparsemat.hpp

@@ -641,6 +641,12 @@ public:
        error if '*this' is finalized and 'B' has larger sparsity pattern. */
    SparseMatrix &operator+=(const SparseMatrix &B);
 
+   /** @brief Add the sparse sub-matrix 'B' to '*this' with the row offset 'ioff'
+       and column offset 'joff'. This operation will cause an error
+       if '*this' is finalized and 'B' has larger sparsity pattern. */
+   void AddSubSparseMatrix(const real_t a, const SparseMatrix &B, int ioff=0,
+                           int joff=0);
+
    /** @brief Add the sparse matrix 'B' scaled by the scalar 'a' into '*this'.
        Only entries in the sparsity pattern of '*this' are added. */
    void Add(const real_t a, const SparseMatrix &B);

linalg/vector.hpp

@@ -698,23 +698,27 @@ inline real_t Distance(const Vector &x, const Vector &y)
 
 inline real_t Vector::DistanceSquaredTo(const real_t *p) const
 {
+   MFEM_ASSERT(data.HostIsValid(), "Data not on host.");
    return DistanceSquared(data, p, size);
 }
 
 inline real_t Vector::DistanceSquaredTo(const Vector &p) const
 {
    MFEM_ASSERT(p.Size() == Size(), "Incompatible vector sizes.");
+   MFEM_ASSERT(data.HostIsValid() && p.data.HostIsValid(), "Data not on host.");
    return DistanceSquared(data, p.data, size);
 }
 
 inline real_t Vector::DistanceTo(const real_t *p) const
 {
+   MFEM_ASSERT(data.HostIsValid(), "Data not on host.");
    return Distance(data, p, size);
 }
 
 inline real_t Vector::DistanceTo(const Vector &p) const
 {
    MFEM_ASSERT(p.Size() == Size(), "Incompatible vector sizes.");
+   MFEM_ASSERT(data.HostIsValid() && p.data.HostIsValid(), "Data not on host.");
    return Distance(data, p.data, size);
 }
 

mesh/mesh.cpp

@@ -1101,37 +1101,41 @@ void Mesh::GetFaceElementTransformations(int FaceNo,
 #endif
 }
 
-FaceElementTransformations *Mesh::GetInteriorFaceTransformations(int FaceNo)
+FaceElementTransformations *Mesh::GetInteriorFaceTransformations(int FaceNo,
+                                                                 int mask)
 {
    GetInteriorFaceTransformations(FaceNo, FaceElemTr, Transformation,
-                                  Transformation2);
+                                  Transformation2, mask);
    return (FaceElemTr.geom == Geometry::INVALID) ? nullptr : &FaceElemTr;
 }
 
 void Mesh::GetInteriorFaceTransformations(int FaceNo,
                                           FaceElementTransformations &FElTr,
                                           IsoparametricTransformation &ElTr1,
-                                          IsoparametricTransformation &ElTr2) const
+                                          IsoparametricTransformation &ElTr2,
+                                          int mask) const
 {
    if (faces_info[FaceNo].Elem2No < 0)
    {
       FElTr.SetGeometryType(Geometry::INVALID);
       return;
    }
-   GetFaceElementTransformations(FaceNo, FElTr, ElTr1, ElTr2);
+   GetFaceElementTransformations(FaceNo, FElTr, ElTr1, ElTr2, mask);
 }
 
-FaceElementTransformations *Mesh::GetBdrFaceTransformations(int BdrElemNo)
+FaceElementTransformations *Mesh::GetBdrFaceTransformations(int BdrElemNo,
+                                                            int mask)
 {
    GetBdrFaceTransformations(BdrElemNo, FaceElemTr, Transformation,
-                             Transformation2);
+                             Transformation2, mask);
    return (FaceElemTr.geom == Geometry::INVALID) ? nullptr : &FaceElemTr;
 }
 
 void Mesh::GetBdrFaceTransformations(int BdrElemNo,
                                      FaceElementTransformations &FElTr,
                                      IsoparametricTransformation &ElTr1,
-                                     IsoparametricTransformation &ElTr2) const
+                                     IsoparametricTransformation &ElTr2,
+                                     int mask) const
 {
    // Check if the face is interior, shared, or nonconforming.
    int fn = GetBdrElementFaceIndex(BdrElemNo);
@@ -1140,7 +1144,7 @@ void Mesh::GetBdrFaceTransformations(int BdrElemNo,
       FElTr.SetGeometryType(Geometry::INVALID);
       return;
    }
-   GetFaceElementTransformations(fn, FElTr, ElTr1, ElTr2, 21);
+   GetFaceElementTransformations(fn, FElTr, ElTr1, ElTr2, mask);
    FElTr.Attribute = boundary[BdrElemNo]->GetAttribute();
    FElTr.ElementNo = BdrElemNo;
    FElTr.ElementType = ElementTransformation::BDR_FACE;

mesh/mesh.hpp

@@ -1816,28 +1816,32 @@ public:
    /// @note The returned object is owned by the class and is shared, i.e.,
    /// calling this function resets pointers obtained from previous calls.
    /// Also, this pointer should NOT be deleted by the caller.
-   FaceElementTransformations *GetInteriorFaceTransformations(int FaceNo);
+   FaceElementTransformations *GetInteriorFaceTransformations(int FaceNo,
+                                                              int mask = 31);
 
    /// @brief Variant of GetInteriorFaceTransformations using a user allocated
    /// FaceElementTransformations object.
    void GetInteriorFaceTransformations(int FaceNo,
                                        FaceElementTransformations &FElTr,
                                        IsoparametricTransformation &ElTr1,
-                                       IsoparametricTransformation &ElTr2) const;
+                                       IsoparametricTransformation &ElTr2,
+                                       int mask = 31) const;
 
    /// @brief Builds the transformation defining the given boundary face.
    ///
    /// @note The returned object is owned by the class and is shared, i.e.,
    /// calling this function resets pointers obtained from previous calls.
    /// Also, this pointer should NOT be deleted by the caller.
-   FaceElementTransformations *GetBdrFaceTransformations(int BdrElemNo);
+   FaceElementTransformations *GetBdrFaceTransformations(int BdrElemNo,
+                                                         int mask = 21);
 
    /// @brief Variant of GetBdrFaceTransformations using a user allocated
    /// FaceElementTransformations object.
    void GetBdrFaceTransformations(int BdrElemNo,
                                   FaceElementTransformations &FElTr,
                                   IsoparametricTransformation &ElTr1,
-                                  IsoparametricTransformation &ElTr2) const;
+                                  IsoparametricTransformation &ElTr2,
+                                  int mask = 21) const;
 
    /// @}