Fixing copyright year

examples/test_surf_diff.cpp

@@ -0,0 +1,322 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+Mesh * GetMesh(int type);
+
+void trans(const Vector &x, Vector &r);
+
+void sigmaFunc(const Vector &x, DenseMatrix &s);
+
+double uExact(const Vector &x)
+{
+   return (0.25 * (2.0 + x[0]) - x[2]) * (x[2] + 0.25 * (2.0 + x[0]));
+}
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   int order = 3;
+   int mesh_type = 4; // Default to Quadrilateral mesh
+   int ref_levels = 0;
+   bool static_cond = false;
+   bool visualization = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_type, "-mt", "--mesh-type",
+                  "Mesh type: 3 - Triangular, 4 - Quadrilateral.");
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly in serial.");
+   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 = GetMesh(mesh_type);
+   int dim = mesh->Dimension();
+
+   // 3. Refine the mesh to increase the resolution. In this example we do
+   //    'ref_levels' of uniform refinement.
+   for (int l = 0; l < ref_levels; l++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   mesh->SetCurvature(3);
+   mesh->Transform(trans);
+
+   // 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;
+   bool delete_fec;
+   if (order > 0)
+   {
+      fec = new H1_FECollection(order, dim);
+      delete_fec = true;
+   }
+   else if (mesh->GetNodes())
+   {
+      fec = mesh->GetNodes()->OwnFEC();
+      delete_fec = false;
+      cout << "Using isoparametric FEs: " << fec->Name() << endl;
+   }
+   else
+   {
+      fec = new H1_FECollection(order = 1, dim);
+      delete_fec = true;
+   }
+   FiniteElementSpace fespace(mesh, 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(&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(&fespace);
+   MatrixFunctionCoefficient sigma(3, sigmaFunc);
+   a.AddDomainIntegrator(new DiffusionIntegrator(sigma));
+
+   // 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();
+
+   OperatorPtr A;
+   Vector B, X;
+   a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
+
+   cout << "Size of linear system: " << A->Height() << endl;
+
+   // 10. Solve the linear system A X = B.
+   if (!pa)
+   {
+      // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
+      GSSmoother M((SparseMatrix&)(*A));
+      PCG(*A, M, B, X, 1, 200, 1e-12, 0.0);
+   }
+   else // Jacobi preconditioning in partial assembly mode
+   {
+      if (UsesTensorBasis(fespace))
+      {
+         OperatorJacobiSmoother M(a, ess_tdof_list);
+         PCG(*A, M, B, X, 1, 400, 1e-12, 0.0);
+      }
+      else
+      {
+         CG(*A, B, X, 1, 400, 1e-12, 0.0);
+      }
+   }
+
+   // 11. Recover the solution as a finite element grid function.
+   a.RecoverFEMSolution(X, b, x);
+
+   FunctionCoefficient uCoef(uExact);
+   double err = x.ComputeL2Error(uCoef);
+
+   mfem::out << "|u - u_h|_2 = " << err << endl;
+
+   // 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. Free the used memory.
+   if (delete_fec)
+   {
+      delete fec;
+   }
+   delete mesh;
+
+   return 0;
+}
+
+Mesh * GetMesh(int type)
+{
+   Mesh * mesh = NULL;
+
+   if (type == 3)
+   {
+      mesh = new Mesh(2, 12, 16, 8, 3);
+
+      mesh->AddVertex(-1.0, -1.0, 0.0);
+      mesh->AddVertex( 1.0, -1.0, 0.0);
+      mesh->AddVertex( 1.0,  1.0, 0.0);
+      mesh->AddVertex(-1.0,  1.0, 0.0);
+      mesh->AddVertex(-1.0, -1.0, 1.0);
+      mesh->AddVertex( 1.0, -1.0, 1.0);
+      mesh->AddVertex( 1.0,  1.0, 1.0);
+      mesh->AddVertex(-1.0,  1.0, 1.0);
+      mesh->AddVertex( 0.0, -1.0, 0.5);
+      mesh->AddVertex( 1.0,  0.0, 0.5);
+      mesh->AddVertex( 0.0,  1.0, 0.5);
+      mesh->AddVertex(-1.0,  0.0, 0.5);
+
+      mesh->AddTriangle(0, 1, 8);
+      mesh->AddTriangle(1, 5, 8);
+      mesh->AddTriangle(5, 4, 8);
+      mesh->AddTriangle(4, 0, 8);
+      mesh->AddTriangle(1, 2, 9);
+      mesh->AddTriangle(2, 6, 9);
+      mesh->AddTriangle(6, 5, 9);
+      mesh->AddTriangle(5, 1, 9);
+      mesh->AddTriangle(2, 3, 10);
+      mesh->AddTriangle(3, 7, 10);
+      mesh->AddTriangle(7, 6, 10);
+      mesh->AddTriangle(6, 2, 10);
+      mesh->AddTriangle(3, 0, 11);
+      mesh->AddTriangle(0, 4, 11);
+      mesh->AddTriangle(4, 7, 11);
+      mesh->AddTriangle(7, 3, 11);
+
+      mesh->AddBdrSegment(0, 1, 1);
+      mesh->AddBdrSegment(1, 2, 1);
+      mesh->AddBdrSegment(2, 3, 1);
+      mesh->AddBdrSegment(3, 0, 1);
+      mesh->AddBdrSegment(5, 4, 2);
+      mesh->AddBdrSegment(6, 5, 2);
+      mesh->AddBdrSegment(7, 6, 2);
+      mesh->AddBdrSegment(4, 7, 2);
+   }
+   else if (type == 4)
+   {
+      mesh = new Mesh(2, 8, 4, 8, 3);
+
+      mesh->AddVertex(-1.0, -1.0, 0.0);
+      mesh->AddVertex( 1.0, -1.0, 0.0);
+      mesh->AddVertex( 1.0,  1.0, 0.0);
+      mesh->AddVertex(-1.0,  1.0, 0.0);
+      mesh->AddVertex(-1.0, -1.0, 1.0);
+      mesh->AddVertex( 1.0, -1.0, 1.0);
+      mesh->AddVertex( 1.0,  1.0, 1.0);
+      mesh->AddVertex(-1.0,  1.0, 1.0);
+
+      mesh->AddQuad(0, 1, 5, 4);
+      mesh->AddQuad(1, 2, 6, 5);
+      mesh->AddQuad(2, 3, 7, 6);
+      mesh->AddQuad(3, 0, 4, 7);
+
+      mesh->AddBdrSegment(0, 1, 1);
+      mesh->AddBdrSegment(1, 2, 1);
+      mesh->AddBdrSegment(2, 3, 1);
+      mesh->AddBdrSegment(3, 0, 1);
+      mesh->AddBdrSegment(5, 4, 2);
+      mesh->AddBdrSegment(6, 5, 2);
+      mesh->AddBdrSegment(7, 6, 2);
+      mesh->AddBdrSegment(4, 7, 2);
+   }
+   else
+   {
+      MFEM_ABORT("Unrecognized mesh type " << type << "!");
+   }
+   mesh->FinalizeTopology();
+
+   return mesh;
+}
+
+void trans(const Vector &x, Vector &r)
+{
+   r.SetSize(3);
+
+   double tol = 1e-6;
+   double theta = 0.0;
+   if (fabs(x[1] + 1.0) < tol)
+   {
+      theta = 0.25 * M_PI * (x[0] - 2.0);
+   }
+   else if (fabs(x[0] - 1.0) < tol)
+   {
+      theta = 0.25 * M_PI * x[1];
+   }
+   else if (fabs(x[1] - 1.0) < tol)
+   {
+      theta = 0.25 * M_PI * (2.0 - x[0]);
+   }
+   else if (fabs(x[0] + 1.0) < tol)
+   {
+      theta = 0.25 * M_PI * (4.0 - x[1]);
+   }
+   else
+   {
+      cout << "side not recognized "
+           << x[0] << " " << x[1] << " " << x[2] << endl;
+   }
+
+   r[0] = cos(theta);
+   r[1] = sin(theta);
+   r[2] = 0.25 * (2.0 * x[2] - 1.0) * (r[0] + 2.0);
+}
+
+void sigmaFunc(const Vector &x, DenseMatrix &s)
+{
+   s.SetSize(3);
+   double a = 17.0 - 2.0 * x[0] * (1.0 + x[0]);
+   s(0,0) = 0.5 + x[0] * x[0] * (8.0 / a - 0.5);
+   s(0,1) = x[0] * x[1] * (8.0 / a - 0.5);
+   s(0,2) = 0.0;
+   s(1,0) = s(0,1);
+   s(1,1) = 0.5 * x[0] * x[0] + 8.0 * x[1] * x[1] / a;
+   s(1,2) = 0.0;
+   s(2,0) = 0.0;
+   s(2,1) = 0.0;
+   s(2,2) = a / 32.0;
+}

fem/bilininteg.cpp

@@ -552,13 +552,28 @@ void DiffusionIntegrator::AssembleElementMatrix
    bool square = (dim == spaceDim);
    double w;
 
+   if (VQ)
+   {
+      MFEM_VERIFY(VQ->GetVDim() == spaceDim,
+                  "Unexpected dimension for VectorCoefficient");
+   }
+   if (MQ)
+   {
+      MFEM_VERIFY(MQ->GetWidth() == spaceDim,
+                  "Unexpected width for MatrixCoefficient");
+      MFEM_VERIFY(MQ->GetHeight() == spaceDim,
+                  "Unexpected height for MatrixCoefficient");
+   }
+
 #ifdef MFEM_THREAD_SAFE
-   DenseMatrix dshape(nd,dim), dshapedxt(nd,spaceDim), invdfdx(dim,spaceDim);
+   DenseMatrix dshape(nd, dim), dshapedxt(nd, spaceDim);
+   DenseMatrix dshapedxt_m(nd, MQ ? spaceDim : 0);
    Vector D(VQ ? VQ->GetVDim() : 0);
 #else
-   dshape.SetSize(nd,dim);
-   dshapedxt.SetSize(nd,spaceDim);
-   invdfdx.SetSize(dim,spaceDim);
+   dshape.SetSize(nd, dim);
+   dshapedxt.SetSize(nd, spaceDim);
+   dshapedxt_m.SetSize(nd, MQ ? spaceDim : 0);
+   M.SetSize(MQ ? spaceDim : 0);
    D.SetSize(VQ ? VQ->GetVDim() : 0);
 #endif
    elmat.SetSize(nd);
@@ -579,10 +594,10 @@ void DiffusionIntegrator::AssembleElementMatrix
       Mult(dshape, Trans.AdjugateJacobian(), dshapedxt);
       if (MQ)
       {
-         MQ->Eval(invdfdx, Trans, ip);
-         invdfdx *= w;
-         Mult(dshapedxt, invdfdx, dshape);
-         AddMultABt(dshape, dshapedxt, elmat);
+         MQ->Eval(M, Trans, ip);
+         M *= w;
+         Mult(dshapedxt, M, dshapedxt_m);
+         AddMultABt(dshapedxt_m, dshapedxt, elmat);
       }
       else if (VQ)
       {
@@ -612,10 +627,25 @@ void DiffusionIntegrator::AssembleElementMatrix2(
    bool square = (dim == spaceDim);
    double w;
 
+   if (VQ)
+   {
+      MFEM_VERIFY(VQ->GetVDim() == spaceDim,
+                  "Unexpected dimension for VectorCoefficient");
+   }
+   if (MQ)
+   {
+      MFEM_VERIFY(MQ->GetWidth() == spaceDim,
+                  "Unexpected width for MatrixCoefficient");
+      MFEM_VERIFY(MQ->GetHeight() == spaceDim,
+                  "Unexpected height for MatrixCoefficient");
+   }
+
 #ifdef MFEM_THREAD_SAFE
    DenseMatrix dshape(tr_nd, dim), dshapedxt(tr_nd, spaceDim);
    DenseMatrix te_dshape(te_nd, dim), te_dshapedxt(te_nd, spaceDim);
    DenseMatrix invdfdx(dim, spaceDim);
+   DenseMatrix dshapedxt_m(te_nd, MQ ? spaceDim : 0);
+   DenseMatrix M(MQ ? spaceDim : 0);
    Vector D(VQ ? VQ->GetVDim() : 0);
 #else
    dshape.SetSize(tr_nd, dim);
@@ -623,6 +653,8 @@ void DiffusionIntegrator::AssembleElementMatrix2(
    te_dshape.SetSize(te_nd, dim);
    te_dshapedxt.SetSize(te_nd, spaceDim);
    invdfdx.SetSize(dim, spaceDim);
+   dshapedxt_m.SetSize(te_nd, MQ ? spaceDim : 0);
+   M.SetSize(MQ ? spaceDim : 0);
    D.SetSize(VQ ? VQ->GetVDim() : 0);
 #endif
    elmat.SetSize(te_nd, tr_nd);
@@ -645,10 +677,10 @@ void DiffusionIntegrator::AssembleElementMatrix2(
       // invdfdx, dshape, and te_dshape no longer needed
       if (MQ)
       {
-         MQ->Eval(invdfdx, Trans, ip);
-         invdfdx *= w;
-         Mult(te_dshapedxt, invdfdx, te_dshape);
-         AddMultABt(te_dshape, dshapedxt, elmat);
+         MQ->Eval(M, Trans, ip);
+         M *= w;
+         Mult(te_dshapedxt, M, dshapedxt_m);
+         AddMultABt(dshapedxt_m, dshapedxt, elmat);
       }
       else if (VQ)
       {
@@ -674,24 +706,34 @@ void DiffusionIntegrator::AssembleElementVector(
 {
    int nd = el.GetDof();
    int dim = el.GetDim();
+   int spaceDim = Tr.GetSpaceDim();
    double w;
 
    if (VQ)
    {
-      MFEM_VERIFY(VQ->GetVDim() == dim, "Unexpected dimension for VectorCoefficient");
+      MFEM_VERIFY(VQ->GetVDim() == spaceDim,
+                  "Unexpected dimension for VectorCoefficient");
+   }
+   if (MQ)
+   {
+      MFEM_VERIFY(MQ->GetWidth() == spaceDim,
+                  "Unexpected width for MatrixCoefficient");
+      MFEM_VERIFY(MQ->GetHeight() == spaceDim,
+                  "Unexpected height for MatrixCoefficient");
    }
 
 #ifdef MFEM_THREAD_SAFE
-   DenseMatrix dshape(nd,dim), invdfdx(dim), mq(dim);
+   DenseMatrix dshape(nd,dim), invdfdx(dim, spaceDim), M(MQ ? spaceDim : 0);
    Vector D(VQ ? VQ->GetVDim() : 0);
 #else
    dshape.SetSize(nd,dim);
-   invdfdx.SetSize(dim);
-   mq.SetSize(dim);
+   invdfdx.SetSize(dim, spaceDim);
+   M.SetSize(MQ ? spaceDim : 0);
    D.SetSize(VQ ? VQ->GetVDim() : 0);
 #endif
    vec.SetSize(dim);
-   pointflux.SetSize(dim);
+   vecdxt.SetSize((VQ || MQ) ? spaceDim : 0);
+   pointflux.SetSize(spaceDim);
 
    elvect.SetSize(nd);
 
@@ -718,19 +760,19 @@ void DiffusionIntegrator::AssembleElementVector(
       }
       else
       {
-         dshape.MultTranspose(elfun, pointflux);
-         invdfdx.MultTranspose(pointflux, vec);
+         dshape.MultTranspose(elfun, vec);
+         invdfdx.MultTranspose(vec, vecdxt);
          if (MQ)
          {
-            MQ->Eval(mq, Tr, ip);
-            mq.Mult(vec, pointflux);
+            MQ->Eval(M, Tr, ip);
+            M.Mult(vecdxt, pointflux);
          }
          else
          {
             VQ->Eval(D, Tr, ip);
-            for (int j=0; j<dim; ++j)
+            for (int j=0; j<spaceDim; ++j)
             {
-               pointflux[j] *= D[j];
+               pointflux[j] = D[j] * vecdxt[j];
             }
          }
       }
@@ -750,14 +792,32 @@ void DiffusionIntegrator::ComputeElementFlux
    dim = el.GetDim();
    spaceDim = Trans.GetSpaceDim();
 
+   if (VQ)
+   {
+      MFEM_VERIFY(VQ->GetVDim() == spaceDim,
+                  "Unexpected dimension for VectorCoefficient");
+   }
+   if (MQ)
+   {
+      MFEM_VERIFY(MQ->GetWidth() == spaceDim,
+                  "Unexpected width for MatrixCoefficient");
+      MFEM_VERIFY(MQ->GetHeight() == spaceDim,
+                  "Unexpected height for MatrixCoefficient");
+   }
+
 #ifdef MFEM_THREAD_SAFE
    DenseMatrix dshape(nd,dim), invdfdx(dim, spaceDim);
+   DenseMatrix M(MQ ? spaceDim : 0);
+   Vector D(VQ ? VQ->GetVDim() : 0);
 #else
    dshape.SetSize(nd,dim);
    invdfdx.SetSize(dim, spaceDim);
+   M.SetSize(MQ ? spaceDim : 0);
+   D.SetSize(VQ ? VQ->GetVDim() : 0);
 #endif
    vec.SetSize(dim);
-   pointflux.SetSize(spaceDim);
+   vecdxt.SetSize(spaceDim);
+   pointflux.SetSize(MQ ? spaceDim : 0);
 
    const IntegrationRule &ir = fluxelem.GetNodes();
    fnd = ir.GetNPoints();
@@ -771,28 +831,38 @@ void DiffusionIntegrator::ComputeElementFlux
 
       Trans.SetIntPoint (&ip);
       CalcInverse(Trans.Jacobian(), invdfdx);
-      invdfdx.MultTranspose(vec, pointflux);
+      invdfdx.MultTranspose(vec, vecdxt);
 
-      if (!MQ)
+      if (!MQ && !VQ)
       {
          if (Q && with_coef)
          {
-            pointflux *= Q->Eval(Trans,ip);
+            vecdxt *= Q->Eval(Trans,ip);
          }
          for (j = 0; j < spaceDim; j++)
          {
-            flux(fnd*j+i) = pointflux(j);
+            flux(fnd*j+i) = vecdxt(j);
          }
       }
       else
       {
-         // assuming dim == spaceDim
-         MFEM_ASSERT(dim == spaceDim, "TODO");
-         MQ->Eval(invdfdx, Trans, ip);
-         invdfdx.Mult(pointflux, vec);
+         if (MQ)
+         {
+            MQ->Eval(M, Trans, ip);
+            M.Mult(vecdxt, pointflux);
+         }
+         else
+         {
+            VQ->Eval(D, Trans, ip);
+            for (int j=0; j<spaceDim; ++j)
+            {
+               pointflux[j] = D[j] * vecdxt[j];
+            }
+
+         }
          for (j = 0; j < dim; j++)
          {
-            flux(fnd*j+i) = vec(j);
+            flux(fnd*j+i) = pointflux(j);
          }
       }
    }
@@ -807,13 +877,13 @@ double DiffusionIntegrator::ComputeFluxEnergy
    int spaceDim = Trans.GetSpaceDim();
 
 #ifdef MFEM_THREAD_SAFE
-   DenseMatrix mq;
+   DenseMatrix M;
 #endif
 
    shape.SetSize(nd);
    pointflux.SetSize(spaceDim);
-   if (d_energy) { vec.SetSize(dim); }
-   if (MQ) { mq.SetSize(dim); }
+   if (d_energy) { vec.SetSize(spaceDim); }
+   if (MQ) { M.SetSize(spaceDim); }
 
    int order = 2 * fluxelem.GetOrder(); // <--
    const IntegrationRule *ir = &IntRules.Get(fluxelem.GetGeomType(), order);
@@ -846,8 +916,9 @@ double DiffusionIntegrator::ComputeFluxEnergy
       }
       else
       {
-         MQ->Eval(mq, Trans, ip);
-         energy += w * mq.InnerProduct(pointflux, pointflux);
+         MFEM_ASSERT(dim == spaceDim, "TODO");
+         MQ->Eval(M, Trans, ip);
+         energy += w * M.InnerProduct(pointflux, pointflux);
       }
 
       if (d_energy)

fem/bilininteg.hpp

@@ -1902,9 +1902,9 @@ protected:
    MatrixCoefficient *MQ;
 
 private:
-   Vector vec, pointflux, shape;
+   Vector vec, vecdxt, pointflux, shape;
 #ifndef MFEM_THREAD_SAFE
-   DenseMatrix dshape, dshapedxt, invdfdx, mq;
+   DenseMatrix dshape, dshapedxt, invdfdx, M, dshapedxt_m;
    DenseMatrix te_dshape, te_dshapedxt;
    Vector D;
 #endif

fem/tmop.cpp

@@ -656,6 +656,7 @@ double TMOP_Metric_315::EvalW(const DenseMatrix &Jpt) const
    ie.SetJacobian(Jpt.GetData());
    const double c1 = ie.Get_I3b() - 1.0;
    return c1*c1;
+   //return c1*c1*c1*c1*c1*c1;
 }
 
 void TMOP_Metric_315::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
@@ -664,6 +665,7 @@ void TMOP_Metric_315::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
    // P = 2*(I3b - 1)*dI3b
    ie.SetJacobian(Jpt.GetData());
    P.Set(2*(ie.Get_I3b() - 1.0), ie.Get_dI3b());
+   //P.Set(6*(ie.Get_I3b() - 1.0), ie.Get_dI3b());
 }
 
 void TMOP_Metric_315::AssembleH(const DenseMatrix &Jpt,
@@ -677,6 +679,9 @@ void TMOP_Metric_315::AssembleH(const DenseMatrix &Jpt,
    ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
    ie.Assemble_TProd(2*weight, ie.Get_dI3b(), A.GetData());
    ie.Assemble_ddI3b(2*weight*(ie.Get_I3b() - 1.0), A.GetData());
+   //ie.Assemble_TProd(6*weight, ie.Get_dI3b(), A.GetData());
+   //ie.Assemble_ddI3b(6*weight*(ie.Get_I3b() - 1.0), A.GetData());
+   
 }
 
 double TMOP_Metric_316::EvalW(const DenseMatrix &Jpt) const

fem/tmop_tools.cpp

@@ -197,10 +197,19 @@ void SerialAdvectorCGOper::Mult(const Vector &ind, Vector &di_dt) const
    di_dt = 0.0;
    CGSolver lin_solver;
    DSmoother prec;
+   /*
+   FGMRESSolver lin_solver;
+   GMRESSolver prec;
+   prec.SetMaxIter(50);
+   prec.SetRelTol(0.0);
+   prec.SetAbsTol(0.0);    
+   prec.SetOperator(M.SpMat());
+   */
+   
    lin_solver.SetPreconditioner(prec);
    lin_solver.SetOperator(M.SpMat());
    lin_solver.SetRelTol(1e-12); lin_solver.SetAbsTol(0.0);
-   lin_solver.SetMaxIter(100);
+   lin_solver.SetMaxIter(200);
    lin_solver.SetPrintLevel(0);
    lin_solver.Mult(rhs, di_dt);
 }

linalg/densemat.cpp

@@ -2497,13 +2497,13 @@ void MultADBt(const DenseMatrix &A, const Vector &D,
 
 void AddMultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
 {
-#ifdef MFEM_DEBUG
+  #ifdef MFEM_DEBUG
    if (A.Height() != ABt.Height() || B.Height() != ABt.Width() ||
        A.Width() != B.Width())
    {
       mfem_error("AddMultABt(...): dimension mismatch");
    }
-#endif
+  #endif
 
 #ifdef MFEM_USE_LAPACK
    static char transa = 'N', transb = 'T';

miniapps/meshing/mesh-optimizer.cpp

@@ -311,7 +311,7 @@ int main(int argc, char *argv[])
    x0 = x;
 
    // 11. Form the integrator that uses the chosen metric and target.
-   double tauval = -0.1;
+   double tauval = -0.002;
    TMOP_QualityMetric *metric = NULL;
    switch (metric_id)
    {

tests/unit/fem/test_surf_blf.cpp

@@ -0,0 +1,219 @@
+// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "mfem.hpp"
+#include "catch.hpp"
+
+using namespace mfem;
+
+namespace surf_blf
+{
+
+Mesh * GetCylMesh(int type);
+
+void trans_skew_cyl(const Vector &x, Vector &r);
+
+void sigmaFunc(const Vector &x, DenseMatrix &s);
+
+double uExact(const Vector &x)
+{
+   return (0.25 * (2.0 + x[0]) - x[2]) * (x[2] + 0.25 * (2.0 + x[0]));
+}
+
+TEST_CASE("Embedded Surface Diffusion",
+          "[DiffusionIntegrator]")
+{
+   for (int type = (int) Element::TRIANGLE;
+        type <= (int) Element::QUADRILATERAL;
+        type++)
+   {
+      int order = 3;
+
+      Mesh *mesh = GetCylMesh(type);
+      int dim = mesh->Dimension();
+
+      mesh->SetCurvature(3);
+      mesh->Transform(trans_skew_cyl);
+
+      H1_FECollection fec(order, dim);
+      FiniteElementSpace fespace(mesh, &fec);
+
+      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);
+      }
+
+      LinearForm b(&fespace);
+      ConstantCoefficient one(1.0);
+      b.AddDomainIntegrator(new DomainLFIntegrator(one));
+      b.Assemble();
+
+      GridFunction x(&fespace);
+      x = 0.0;
+
+      BilinearForm a(&fespace);
+      MatrixFunctionCoefficient sigma(3, sigmaFunc);
+      a.AddDomainIntegrator(new DiffusionIntegrator(sigma));
+      a.Assemble();
+
+      OperatorPtr A;
+      Vector B, X;
+      a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
+
+      GSSmoother M((SparseMatrix&)(*A));
+      PCG(*A, M, B, X, 0, 200, 1e-12, 0.0);
+
+      a.RecoverFEMSolution(X, b, x);
+
+      FunctionCoefficient uCoef(uExact);
+      double err = x.ComputeL2Error(uCoef);
+
+      REQUIRE(err < 1.5e-3);
+
+      delete mesh;
+   }
+}
+
+Mesh * GetCylMesh(int type)
+{
+   Mesh * mesh = NULL;
+
+   if (type == (int) Element::TRIANGLE)
+   {
+      mesh = new Mesh(2, 12, 16, 8, 3);
+
+      mesh->AddVertex(-1.0, -1.0, 0.0);
+      mesh->AddVertex( 1.0, -1.0, 0.0);
+      mesh->AddVertex( 1.0,  1.0, 0.0);
+      mesh->AddVertex(-1.0,  1.0, 0.0);
+      mesh->AddVertex(-1.0, -1.0, 1.0);
+      mesh->AddVertex( 1.0, -1.0, 1.0);
+      mesh->AddVertex( 1.0,  1.0, 1.0);
+      mesh->AddVertex(-1.0,  1.0, 1.0);
+      mesh->AddVertex( 0.0, -1.0, 0.5);
+      mesh->AddVertex( 1.0,  0.0, 0.5);
+      mesh->AddVertex( 0.0,  1.0, 0.5);
+      mesh->AddVertex(-1.0,  0.0, 0.5);
+
+      mesh->AddTriangle(0, 1, 8);
+      mesh->AddTriangle(1, 5, 8);
+      mesh->AddTriangle(5, 4, 8);
+      mesh->AddTriangle(4, 0, 8);
+      mesh->AddTriangle(1, 2, 9);
+      mesh->AddTriangle(2, 6, 9);
+      mesh->AddTriangle(6, 5, 9);
+      mesh->AddTriangle(5, 1, 9);
+      mesh->AddTriangle(2, 3, 10);
+      mesh->AddTriangle(3, 7, 10);
+      mesh->AddTriangle(7, 6, 10);
+      mesh->AddTriangle(6, 2, 10);
+      mesh->AddTriangle(3, 0, 11);
+      mesh->AddTriangle(0, 4, 11);
+      mesh->AddTriangle(4, 7, 11);
+      mesh->AddTriangle(7, 3, 11);
+
+      mesh->AddBdrSegment(0, 1, 1);
+      mesh->AddBdrSegment(1, 2, 1);
+      mesh->AddBdrSegment(2, 3, 1);
+      mesh->AddBdrSegment(3, 0, 1);
+      mesh->AddBdrSegment(5, 4, 2);
+      mesh->AddBdrSegment(6, 5, 2);
+      mesh->AddBdrSegment(7, 6, 2);
+      mesh->AddBdrSegment(4, 7, 2);
+   }
+   else if (type == (int) Element::QUADRILATERAL)
+   {
+      mesh = new Mesh(2, 8, 4, 8, 3);
+
+      mesh->AddVertex(-1.0, -1.0, 0.0);
+      mesh->AddVertex( 1.0, -1.0, 0.0);
+      mesh->AddVertex( 1.0,  1.0, 0.0);
+      mesh->AddVertex(-1.0,  1.0, 0.0);
+      mesh->AddVertex(-1.0, -1.0, 1.0);
+      mesh->AddVertex( 1.0, -1.0, 1.0);
+      mesh->AddVertex( 1.0,  1.0, 1.0);
+      mesh->AddVertex(-1.0,  1.0, 1.0);
+
+      mesh->AddQuad(0, 1, 5, 4);
+      mesh->AddQuad(1, 2, 6, 5);
+      mesh->AddQuad(2, 3, 7, 6);
+      mesh->AddQuad(3, 0, 4, 7);
+
+      mesh->AddBdrSegment(0, 1, 1);
+      mesh->AddBdrSegment(1, 2, 1);
+      mesh->AddBdrSegment(2, 3, 1);
+      mesh->AddBdrSegment(3, 0, 1);
+      mesh->AddBdrSegment(5, 4, 2);
+      mesh->AddBdrSegment(6, 5, 2);
+      mesh->AddBdrSegment(7, 6, 2);
+      mesh->AddBdrSegment(4, 7, 2);
+   }
+   else
+   {
+      MFEM_ABORT("Unrecognized mesh type " << type << "!");
+   }
+   mesh->FinalizeTopology();
+
+   return mesh;
+}
+
+void trans_skew_cyl(const Vector &x, Vector &r)
+{
+   r.SetSize(3);
+
+   double tol = 1e-6;
+   double theta = 0.0;
+   if (fabs(x[1] + 1.0) < tol)
+   {
+      theta = 0.25 * M_PI * (x[0] - 2.0);
+   }
+   else if (fabs(x[0] - 1.0) < tol)
+   {
+      theta = 0.25 * M_PI * x[1];
+   }
+   else if (fabs(x[1] - 1.0) < tol)
+   {
+      theta = 0.25 * M_PI * (2.0 - x[0]);
+   }
+   else if (fabs(x[0] + 1.0) < tol)
+   {
+      theta = 0.25 * M_PI * (4.0 - x[1]);
+   }
+   else
+   {
+      mfem::out << "side not recognized "
+                << x[0] << " " << x[1] << " " << x[2] << std::endl;
+   }
+
+   r[0] = cos(theta);
+   r[1] = sin(theta);
+   r[2] = 0.25 * (2.0 * x[2] - 1.0) * (r[0] + 2.0);
+}
+
+void sigmaFunc(const Vector &x, DenseMatrix &s)
+{
+   s.SetSize(3);
+   double a = 17.0 - 2.0 * x[0] * (1.0 + x[0]);
+   s(0,0) = 0.5 + x[0] * x[0] * (8.0 / a - 0.5);
+   s(0,1) = x[0] * x[1] * (8.0 / a - 0.5);
+   s(0,2) = 0.0;
+   s(1,0) = s(0,1);
+   s(1,1) = 0.5 * x[0] * x[0] + 8.0 * x[1] * x[1] / a;
+   s(1,2) = 0.0;
+   s(2,0) = 0.0;
+   s(2,1) = 0.0;
+   s(2,2) = a / 32.0;
+}
+
+} // namespace surf_blf