Fixing copyright year

examples/test_surf_diff.cpp

@@ -23,8 +23,6 @@ int main(int argc, char *argv[])
    int mesh_type = 4; // Default to Quadrilateral mesh
    int ref_levels = 0;
    bool static_cond = false;
-   bool pa = false;
-   const char *device_config = "cpu";
    bool visualization = true;
 
    OptionsParser args(argc, argv);
@@ -37,10 +35,6 @@ int main(int argc, char *argv[])
                   " isoparametric space.");
    args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
                   "--no-static-condensation", "Enable static condensation.");
-   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
-                  "--no-partial-assembly", "Enable Partial Assembly.");
-   args.AddOption(&device_config, "-d", "--device",
-                  "Device configuration string, see Device::Configure().");
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
@@ -52,18 +46,13 @@ int main(int argc, char *argv[])
    }
    args.PrintOptions(cout);
 
-   // 2. Enable hardware devices such as GPUs, and programming models such as
-   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
-   Device device(device_config);
-   device.Print();
-
-   // 3. Read the mesh from the given mesh file. We can handle triangular,
+   // 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();
 
-   // 4. Refine the mesh to increase the resolution. In this example we do
+   // 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++)
    {
@@ -73,7 +62,7 @@ int main(int argc, char *argv[])
    mesh->SetCurvature(3);
    mesh->Transform(trans);
 
-   // 5. Define a finite element space on the mesh. Here we use continuous
+   // 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;
@@ -98,7 +87,7 @@ int main(int argc, char *argv[])
    cout << "Number of finite element unknowns: "
         << fespace.GetTrueVSize() << endl;
 
-   // 6. Determine the list of true (i.e. conforming) essential boundary dofs.
+   // 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.
@@ -110,7 +99,7 @@ int main(int argc, char *argv[])
       fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
    }
 
-   // 7. Set up the linear form b(.) which corresponds to the right-hand side of
+   // 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);
@@ -118,24 +107,23 @@ int main(int argc, char *argv[])
    b.AddDomainIntegrator(new DomainLFIntegrator(one));
    b.Assemble();
 
-   // 8. Define the solution vector x as a finite element grid function
+   // 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;
 
-   // 9. Set up the bilinear form a(.,.) on the finite element space
+   // 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);
-   if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
    a.AddDomainIntegrator(new DiffusionIntegrator(sigma));
 
-   // 10. 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.
+   // 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();
 
@@ -145,7 +133,7 @@ int main(int argc, char *argv[])
 
    cout << "Size of linear system: " << A->Height() << endl;
 
-   // 11. Solve the linear system A X = B.
+   // 10. Solve the linear system A X = B.
    if (!pa)
    {
       // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
@@ -165,7 +153,7 @@ int main(int argc, char *argv[])
       }
    }
 
-   // 12. Recover the solution as a finite element grid function.
+   // 11. Recover the solution as a finite element grid function.
    a.RecoverFEMSolution(X, b, x);
 
    FunctionCoefficient uCoef(uExact);
@@ -173,7 +161,7 @@ int main(int argc, char *argv[])
 
    mfem::out << "|u - u_h|_2 = " << err << endl;
 
-   // 13. Save the refined mesh and the solution. This output can be viewed later
+   // 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);
@@ -182,7 +170,7 @@ int main(int argc, char *argv[])
    sol_ofs.precision(8);
    x.Save(sol_ofs);
 
-   // 14. Send the solution by socket to a GLVis server.
+   // 13. Send the solution by socket to a GLVis server.
    if (visualization)
    {
       char vishost[] = "localhost";
@@ -192,7 +180,7 @@ int main(int argc, char *argv[])
       sol_sock << "solution\n" << *mesh << x << flush;
    }
 
-   // 15. Free the used memory.
+   // 14. Free the used memory.
    if (delete_fec)
    {
       delete fec;

fem/bilininteg.cpp

@@ -552,20 +552,33 @@ 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);
-   DenseMatrix invdfdx(spaceDim, 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(spaceDim, spaceDim);
+   dshapedxt_m.SetSize(nd, MQ ? spaceDim : 0);
+   M.SetSize(MQ ? spaceDim : 0);
    D.SetSize(VQ ? VQ->GetVDim() : 0);
 #endif
    elmat.SetSize(nd);
 
    const IntegrationRule *ir = IntRule ? IntRule : &GetRule(el, el);
-   DenseMatrix dshapedxt_m(nd, spaceDim);
 
    elmat = 0.0;
    for (int i = 0; i < ir->GetNPoints(); i++)
@@ -581,9 +594,9 @@ void DiffusionIntegrator::AssembleElementMatrix
       Mult(dshape, Trans.AdjugateJacobian(), dshapedxt);
       if (MQ)
       {
-         MQ->Eval(invdfdx, Trans, ip);
-         invdfdx *= w;
-         Mult(dshapedxt, invdfdx, dshapedxt_m);
+         MQ->Eval(M, Trans, ip);
+         M *= w;
+         Mult(dshapedxt, M, dshapedxt_m);
          AddMultABt(dshapedxt_m, dshapedxt, elmat);
       }
       else if (VQ)
@@ -614,23 +627,39 @@ 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(spaceDim, 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);
    dshapedxt.SetSize(tr_nd, spaceDim);
    te_dshape.SetSize(te_nd, dim);
    te_dshapedxt.SetSize(te_nd, spaceDim);
-   invdfdx.SetSize(spaceDim, 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);
 
    const IntegrationRule *ir = IntRule ? IntRule : &GetRule(trial_fe, test_fe);
-   DenseMatrix te_dshapedxt_m(te_nd, spaceDim);
 
    elmat = 0.0;
    for (int i = 0; i < ir->GetNPoints(); i++)
@@ -648,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_dshapedxt_m);
-         AddMultABt(te_dshapedxt_m, dshapedxt, elmat);
+         MQ->Eval(M, Trans, ip);
+         M *= w;
+         Mult(te_dshapedxt, M, dshapedxt_m);
+         AddMultABt(dshapedxt_m, dshapedxt, elmat);
       }
       else if (VQ)
       {
@@ -677,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);
 
@@ -721,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];
             }
          }
       }
@@ -753,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();
@@ -774,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);
          }
       }
    }
@@ -810,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);
@@ -849,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

@@ -1,4 +1,4 @@
-// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// 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.
 //