correct essential boundary comment

Co-authored-by: Socratis Petrides <petrides1@llnl.gov>

doc/CodeDocumentation.dox

@@ -119,6 +119,8 @@ namespace mfem {
  * - <a class="el" href="ex40p_8cpp_source.html">Example 40p</a>: parallel eikonal equation
  * - <a class="el" href="ex41_8cpp_source.html">Example 41</a>: DG/CG IMEX time-dependent advection-diffusion
  * - <a class="el" href="ex41p_8cpp_source.html">Example 41p</a>: parallel DG/CG IMEX time-dependent advection-diffusion
+ * - <a class="el" href="ex43_8cpp_source.html">Example 43</a>: sliding boundary conditions in linear elasticity
+ * - <a class="el" href="ex43p_8cpp_source.html">Example 43p</a>: parallel sliding boundary conditions in linear elasticity
  *
  * <H4>AmgX Examples</H4>
  * - Variants of Examples

examples/CMakeLists.txt

@@ -47,6 +47,7 @@ list(APPEND ALL_EXE_SRCS
   ex39.cpp
   ex40.cpp
   ex41.cpp
+  ex43.cpp
   )
 
 if (MFEM_USE_MPI)
@@ -91,6 +92,7 @@ if (MFEM_USE_MPI)
     ex39p.cpp
     ex40p.cpp
     ex41p.cpp
+    ex43p.cpp
   )
 endif()
 

examples/ex43.cpp

@@ -0,0 +1,278 @@
+//                                MFEM Example 43
+//
+// Compile with: make ex43
+//
+// Sample runs:  ex43 -m ../data/ball-nurbs.mesh -r 2
+//               ex43 -m ../data/ref-cube.mesh -r 2
+//               ex43 -m ../data/fichera.mesh
+//               ex43 -m ../data/star.mesh
+//
+// Description:  This example code solves a linear elasticity problem using
+//               Nitsche's method to enforce sliding boundary conditions. In
+//               particular, we consider a linear elastic body that is displaced
+//               in the normal direction on the entire boundary, but is free to
+//               slide in the tangential direction. This is achieved by imposing
+//               homogeneous Dirichlet boundary conditions on the normal
+//               component of the displacement, while applying homogeneous
+//               Neumann boundary conditions on the tangential components of the
+//               displacement. By enforcing a uniform, constant normal
+//               displacement on the boundary, we can simulate the effect of
+//               compressing or expanding the elastic body uniformly. These
+//               boundary conditions are applied weakly using Nitsche's method,
+//               allowing for more flexibility in handling complex geometries in
+//               either 2D or 3D.
+//
+//               The strong form is given by:
+//
+//                          −Div(σ(u)) = 0      in Ω
+//                               u ⋅ n = g      on Γ
+//                                σ(u) ⊥ n      on Γ
+//
+//               where σ(u) = λ tr(ε(u)) I + 2μ ε(u) is the stress tensor, ε(u)
+//               is the strain tensor, λ and μ are the Lamé parameters, and g is
+//               the prescribed displacement on the boundary. Here, n is the
+//               outward normal on the boundary Γ = ∂Ω.
+//
+//               The weak form using Nitsche's method is:
+//
+//               Find u ∈ V     such that   a(u,v) = b(v)   for all v ∈ V
+//
+//               where
+//
+//                          a(u,v) := ∫_Ω σ(u) : ε(v) dx
+//                                    - ∫_Γ (σ(u) n ⋅ n) (v ⋅ n) dS
+//                                    - ∫_Γ (σ(v) n ⋅ n) (u ⋅ n) dS
+//                                    + κ ∫_Γ h⁻¹ (λ + 2μ) (u ⋅ n) (v ⋅ n) dS,
+//
+//                            b(v) := - ∫_Γ σ(v) n ⋅ n g dS
+//                                    + κ ∫_Γ h⁻¹ (λ + 2μ) (v ⋅ n) g dS,
+//
+//               with κ > 0 being a penalty parameter. Here, h is a
+//               characteristic element size on the boundary. The function
+//               space V is a vector H1-conforming finite element space.
+//
+//               This example can be viewed as an alternative to Example 28.
+//               Whereas Example 28 imposes sliding boundary conditions using
+//               the general-purpose constrained system solvers found in
+//               mfem/linalg/constraints.hpp, this example employs Nitsche's
+//               method to weakly enforce the same condition by modifying the
+//               underlying variational formulation. Unlike Example 28, the
+//               approach here is specialized to isotropic linear elasticity,
+//               but it has the advantage of producing a well-conditioned SPD
+//               stiffness matrix that can be readily preconditioned with
+//               standard AMG. We recommend reviewing Example 2 before working
+//               through this example.
+
+#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";
+   real_t displ_mag = 0.1;
+   int order = 1;
+   int ref_levels = 0;
+   real_t lambda = 1.0;
+   real_t mu = 1.0;
+   real_t kappa = -1.0;
+   bool static_cond = false;
+   bool visualization = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&displ_mag, "-g", "--displ",
+                  "Magnitude of the normal displacement.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&ref_levels, "-r", "--ref_levels",
+                  "Number of uniform mesh refinements.");
+   args.AddOption(&lambda, "-l", "--lambda", "First Lamé parameter.");
+   args.AddOption(&mu, "-mu", "--mu", "Second Lamé parameter.");
+   args.AddOption(&kappa, "-k", "--kappa",
+                  "The penalty parameter, should be positive."
+                  " Negative values are replaced with (order+1)^2.");
+   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;
+   }
+   if (kappa < 0)
+   {
+      kappa = (order+1)*(order+1);
+   }
+   args.PrintOptions(cout);
+
+   // 2. Read the mesh from the given mesh file. We can handle triangular,
+   //    quadrilateral, tetrahedral or hexahedral elements with the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 3. Select the order of the finite element discretization space. For NURBS
+   //    meshes, we increase the order by degree elevation.
+   if (mesh->NURBSext)
+   {
+      mesh->DegreeElevate(order, order);
+   }
+
+   // 4. Refine the mesh to increase the resolution. In this example we do
+   //    'ref_levels' of uniform refinement.
+   for (int i = 0; i < ref_levels; i++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 5. Interpolate the geometry after refinement to control geometry error.
+   int curvature_order = max(order, 2);
+   mesh->SetCurvature(curvature_order);
+
+   // 6. Define a finite element space on the mesh. Here we use vector finite
+   //    elements, i.e. dim copies of a scalar finite element space. The vector
+   //    dimension is specified by the last argument of the FiniteElementSpace
+   //    constructor. For NURBS meshes, we use the (degree elevated) NURBS space
+   //    associated with the mesh nodes.
+   FiniteElementCollection *fec;
+   FiniteElementSpace *fespace;
+   if (mesh->NURBSext)
+   {
+      fec = NULL;
+      fespace = mesh->GetNodes()->FESpace();
+   }
+   else
+   {
+      fec = new H1_FECollection(order, dim);
+      fespace = new FiniteElementSpace(mesh, fec, dim);
+   }
+   cout << "Number of finite element unknowns: " << fespace->GetTrueVSize()
+        << endl << "Assembling: " << flush;
+
+   // 7. Mark the boundary attributes where the sliding (Nitsche) boundary
+   //    conditions are to be applied. These b.c. are imposed weakly, by adding
+   //    the appropriate boundary integrators over the marked 'ess_bdr' to the
+   //    bilinear and linear forms. Thus, no dofs are eliminated; there are no
+   //    essential boundary conditions.
+   Array<int> ess_tdof_list, ess_bdr(mesh->bdr_attributes.Max());
+   ess_bdr = 1;
+
+   // 8. 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
+   //    corresponding to the linear elasticity integrator with constant
+   //    coefficients lambda and mu.
+   ConstantCoefficient lambda_c(lambda);
+   ConstantCoefficient mu_c(mu);
+
+   BilinearForm *a = new BilinearForm(fespace);
+   a->AddDomainIntegrator(new ElasticityIntegrator(lambda_c,mu_c));
+   a->AddBdrFaceIntegrator(
+      new SlidingElasticityIntegrator(lambda_c, mu_c, kappa),
+      ess_bdr);
+
+   // 10. Set up the linear form b(.) corresponding to the Nitsche method
+   //     to impose the Dirichlet boundary conditions. Here, we set the
+   //     prescribed displacement on the Dirichlet boundary to be a constant
+   //     normal displacement of magnitude 'displ_mag'.
+   ConstantCoefficient g(displ_mag);
+
+   LinearForm *b = new LinearForm(fespace);
+   b->AddBdrFaceIntegrator(
+      new SlidingElasticityLFIntegrator(
+         g, lambda_c, mu_c, kappa), ess_bdr);
+   b->Assemble();
+
+   // 11. 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.
+   cout << "matrix ... " << flush;
+   if (static_cond) { a->EnableStaticCondensation(); }
+   a->Assemble();
+
+   SparseMatrix A;
+   Vector B, X;
+   a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
+   cout << "done." << endl;
+
+   cout << "Size of linear system: " << A.Height() << endl;
+
+#ifndef MFEM_USE_SUITESPARSE
+   // 12. Define a simple symmetric Gauss-Seidel preconditioner and use it to
+   //     solve the system Ax=b with PCG.
+   GSSmoother M(A);
+   PCG(A, M, B, X, 1, 500, 1e-12, 0.0);
+#else
+   // 12. 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
+
+   // 13. Recover the solution as a finite element grid function.
+   a->RecoverFEMSolution(X, *b, x);
+
+   // 14. For non-NURBS meshes, make the mesh curved based on the finite element
+   //     space. This means that we define the mesh elements through a fespace
+   //     based transformation of the reference element. This allows us to save
+   //     the displaced mesh as a curved mesh when using high-order finite
+   //     element displacement field. We assume that the initial mesh (read from
+   //     the file) is not higher order curved mesh compared to the chosen FE
+   //     space.
+   if (!mesh->NURBSext)
+   {
+      mesh->SetNodalFESpace(fespace);
+   }
+
+   // 15. Save the displaced mesh and the inverted solution (which gives the
+   //     backward displacements to the original grid). This output can be
+   //     viewed later using GLVis: "glvis -m displaced.mesh -g sol.gf".
+   {
+      GridFunction *nodes = mesh->GetNodes();
+      *nodes += x;
+      x *= -1;
+      ofstream mesh_ofs("displaced.mesh");
+      mesh_ofs.precision(8);
+      mesh->Print(mesh_ofs);
+      ofstream sol_ofs("sol.gf");
+      sol_ofs.precision(8);
+      x.Save(sol_ofs);
+   }
+
+   // 16. Send the above data by socket to a GLVis server. Use the "n" and "b"
+   //     keys in GLVis to visualize the displacements.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << *mesh << x << flush;
+   }
+
+   // 17. Free the used memory.
+   delete a;
+   delete b;
+   if (fec)
+   {
+      delete fespace;
+      delete fec;
+   }
+   delete mesh;
+
+   return 0;
+}

examples/ex43p.cpp

@@ -0,0 +1,332 @@
+//                                MFEM Example 43 - Parallel Version
+//
+// Compile with: make ex43p
+//
+// Sample runs:  mpirun -np 4 ex43p -m ../data/ball-nurbs.mesh -r 2
+//               mpirun -np 4 ex43p -m ../data/ref-cube.mesh -r 2
+//               mpirun -np 4 ex43p -m ../data/fichera.mesh
+//               mpirun -np 4 ex43p -m ../data/star.mesh
+//
+// Description:  This example code solves a linear elasticity problem using
+//               Nitsche's method to enforce sliding boundary conditions. In
+//               particular, we consider a linear elastic body that is displaced
+//               in the normal direction on the entire boundary, but is free to
+//               slide in the tangential direction. This is achieved by imposing
+//               homogeneous Dirichlet boundary conditions on the normal
+//               component of the displacement, while applying homogeneous
+//               Neumann boundary conditions on the tangential components of the
+//               displacement. By enforcing a uniform, constant normal
+//               displacement on the boundary, we can simulate the effect of
+//               compressing or expanding the elastic body uniformly. These
+//               boundary conditions are applied weakly using Nitsche's method,
+//               allowing for more flexibility in handling complex geometries in
+//               either 2D or 3D.
+//
+//               The strong form is given by:
+//
+//                          −Div(σ(u)) = 0      in Ω
+//                               u ⋅ n = g      on Γ
+//                                σ(u) ⊥ n      on Γ
+//
+//               where σ(u) = λ tr(ε(u)) I + 2μ ε(u) is the stress tensor, ε(u)
+//               is the strain tensor, λ and μ are the Lamé parameters, and g is
+//               the prescribed displacement on the boundary. Here, n is the
+//               outward normal on the boundary Γ = ∂Ω.
+//
+//               The weak form using Nitsche's method is:
+//
+//               Find u ∈ V     such that   a(u,v) = b(v)   for all v ∈ V
+//
+//               where
+//
+//                          a(u,v) := ∫_Ω σ(u) : ε(v) dx
+//                                    - ∫_Γ (σ(u) n ⋅ n) (v ⋅ n) dS
+//                                    - ∫_Γ (σ(v) n ⋅ n) (u ⋅ n) dS
+//                                    + κ ∫_Γ h⁻¹ (λ + 2μ) (u ⋅ n) (v ⋅ n) dS,
+//
+//                            b(v) := - ∫_Γ σ(v) n ⋅ n g dS
+//                                    + κ ∫_Γ h⁻¹ (λ + 2μ) (v ⋅ n) g dS,
+//
+//               with κ > 0 being a penalty parameter. Here, h is a
+//               characteristic element size on the boundary. The function
+//               space V is a vector H1-conforming finite element space.
+//
+//               This example can be viewed as an alternative to Example 28.
+//               Whereas Example 28 imposes sliding boundary conditions using
+//               the general-purpose constrained system solvers found in
+//               mfem/linalg/constraints.hpp, this example employs Nitsche's
+//               method to weakly enforce the same condition by modifying the
+//               underlying variational formulation. Unlike Example 28, the
+//               approach here is specialized to isotropic linear elasticity,
+//               but it has the advantage of producing a well-conditioned SPD
+//               stiffness matrix that can be readily preconditioned with
+//               standard AMG. We recommend reviewing Example 2 before working
+//               through this example.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI and HYPRE.
+   Mpi::Init(argc, argv);
+   int num_procs = Mpi::WorldSize();
+   int myid = Mpi::WorldRank();
+   Hypre::Init();
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../data/star.mesh";
+   real_t displ_mag = 0.1;
+   int order = 1;
+   int ref_levels = 0;
+   real_t lambda = 1.0;
+   real_t mu = 1.0;
+   real_t kappa = -1.0;
+   bool static_cond = false;
+   bool reorder_space = false;
+   bool visualization = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&displ_mag, "-g", "--displ",
+                  "Magnitude of the normal displacement.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&ref_levels, "-r", "--ref_levels",
+                  "Number of uniform mesh refinements.");
+   args.AddOption(&lambda, "-l", "--lambda", "First Lamé parameter.");
+   args.AddOption(&mu, "-mu", "--mu", "Second Lamé parameter.");
+   args.AddOption(&kappa, "-k", "--kappa",
+                  "The penalty parameter, should be positive."
+                  " Negative values are replaced with (order+1)^2.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&reorder_space, "-nodes", "--by-nodes", "-vdim", "--by-vdim",
+                  "Use byNODES ordering of vector space instead of byVDIM");
+   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);
+      }
+      return 1;
+   }
+   if (kappa < 0)
+   {
+      kappa = (order+1)*(order+1);
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // 3. Read the (serial) mesh from the given mesh file. We can handle triangular,
+   //    quadrilateral, tetrahedral or hexahedral elements with the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 4. Select the order of the finite element discretization space. For NURBS
+   //    meshes, we increase the order by degree elevation.
+   if (mesh->NURBSext)
+   {
+      mesh->DegreeElevate(order, order);
+   }
+
+   // 5. Refine the mesh to increase the resolution. In this example we do
+   //    'ref_levels' of uniform refinement.
+   for (int i = 0; i < ref_levels; i++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 6. Interpolate the geometry after refinement to control geometry error.
+   int curvature_order = max(order, 2);
+   mesh->SetCurvature(curvature_order);
+
+   // 7. 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;
+
+   // 8. Define a finite element space on the mesh. Here we use vector finite
+   //    elements, i.e. dim copies of a scalar finite element space. The vector
+   //    dimension is specified by the last argument of the FiniteElementSpace
+   //    constructor. For NURBS meshes, we use the (degree elevated) NURBS space
+   //    associated with the mesh nodes.
+   FiniteElementCollection *fec;
+   ParFiniteElementSpace *fespace;
+   const bool use_nodal_fespace = pmesh->NURBSext;
+   if (use_nodal_fespace)
+   {
+      fec = NULL;
+      fespace = (ParFiniteElementSpace *)pmesh->GetNodes()->FESpace();
+   }
+   else
+   {
+      fec = new H1_FECollection(order, dim);
+      if (reorder_space)
+      {
+         fespace = new ParFiniteElementSpace(pmesh, fec, dim, Ordering::byNODES);
+      }
+      else
+      {
+         fespace = new ParFiniteElementSpace(pmesh, fec, dim, Ordering::byVDIM);
+      }
+   }
+   HYPRE_BigInt size = fespace->GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of finite element unknowns: " << size << endl
+           << "Assembling: " << flush;
+   }
+
+   // 9. Mark the boundary attributes where the sliding (Nitsche) boundary
+   //    conditions are to be applied. These b.c. are imposed weakly, by adding
+   //    the appropriate boundary integrators over the marked 'ess_bdr' to the
+   //    bilinear and linear forms. Thus, no dofs are eliminated; there are no
+   //    essential boundary conditions.
+   Array<int> ess_tdof_list, ess_bdr;
+   if (pmesh->bdr_attributes.Size())
+   {
+      ess_bdr.SetSize(pmesh->bdr_attributes.Max());
+      ess_bdr = 1;
+   }
+
+   // 10. 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.
+   ParGridFunction x(fespace);
+   x = 0.0;
+
+   // 11. Set up the bilinear form a(.,.) on the finite element space
+   //     corresponding to the linear elasticity integrator with constant
+   //     coefficients lambda and mu.
+   ConstantCoefficient lambda_c(lambda);
+   ConstantCoefficient mu_c(mu);
+
+   ParBilinearForm *a = new ParBilinearForm(fespace);
+   a->AddDomainIntegrator(new ElasticityIntegrator(lambda_c,mu_c));
+   a->AddBdrFaceIntegrator(
+      new SlidingElasticityIntegrator(lambda_c, mu_c, kappa),
+      ess_bdr);
+
+   // 12. Set up the linear form b(.) corresponding to the Nitsche method
+   //     to impose the Dirichlet boundary conditions. Here, we set the
+   //     prescribed displacement on the Dirichlet boundary to be a constant
+   //     normal displacement of magnitude 'displ_mag'.
+   ConstantCoefficient g(displ_mag);
+
+   ParLinearForm *b = new ParLinearForm(fespace);
+   b->AddBdrFaceIntegrator(
+      new SlidingElasticityLFIntegrator(
+         g, lambda_c, mu_c, kappa), ess_bdr);
+   b->Assemble();
+
+   // 13. 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 (myid == 0) { cout << "matrix ... " << flush; }
+   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 << "done." << endl;
+      cout << "Size of linear system: " << A.GetGlobalNumRows() << endl;
+   }
+
+   // 14. Define and apply a parallel PCG solver for A X = B with the BoomerAMG
+   //     preconditioner from hypre.
+   HypreBoomerAMG *amg = new HypreBoomerAMG(A);
+   if (!a->StaticCondensationIsEnabled())
+   {
+      amg->SetElasticityOptions(fespace);
+   }
+   else
+   {
+      amg->SetSystemsOptions(dim, reorder_space);
+   }
+   HyprePCG *pcg = new HyprePCG(A);
+   pcg->SetTol(1e-8);
+   pcg->SetMaxIter(500);
+   pcg->SetPrintLevel(2);
+   pcg->SetPreconditioner(*amg);
+   pcg->Mult(B, X);
+
+   // 15. Recover the parallel grid function corresponding to X. This is the
+   //     local finite element solution on each processor.
+   a->RecoverFEMSolution(X, *b, x);
+
+   // 16. For non-NURBS meshes, make the mesh curved based on the finite element
+   //     space. This means that we define the mesh elements through a fespace
+   //     based transformation of the reference element.  This allows us to save
+   //     the displaced mesh as a curved mesh when using high-order finite
+   //     element displacement field. We assume that the initial mesh (read from
+   //     the file) is not higher order curved mesh compared to the chosen FE
+   //     space.
+   if (!use_nodal_fespace)
+   {
+      pmesh->SetNodalFESpace(fespace);
+   }
+
+   // 17. Save in parallel the displaced mesh and the inverted solution (which
+   //     gives the backward displacements to the original grid). This output
+   //     can be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
+   {
+      GridFunction *nodes = pmesh->GetNodes();
+      *nodes += x;
+      x *= -1;
+
+      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);
+   }
+
+   // 18. Send the above data by socket to a GLVis server.  Use the "n" and "b"
+   //     keys in GLVis to visualize the displacements.
+   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;
+   }
+
+   // 19. Free the used memory.
+   delete pcg;
+   delete amg;
+   delete a;
+   delete b;
+   if (fec)
+   {
+      delete fespace;
+      delete fec;
+   }
+   delete pmesh;
+
+   return 0;
+}

examples/makefile

@@ -22,11 +22,11 @@ MFEM_LIB_FILE = mfem_is_not_built
 
 SEQ_EXAMPLES = ex0 ex1 ex2 ex3 ex4 ex5 ex6 ex7 ex8 ex9 ex10 ex14 ex15 ex16 \
    ex17 ex18 ex19 ex20 ex21 ex22 ex23 ex24 ex25 ex26 ex27 ex28 ex29 ex30 \
-   ex31 ex33 ex34 ex36 ex37 ex38 ex39 ex40 ex41
+   ex31 ex33 ex34 ex36 ex37 ex38 ex39 ex40 ex41 ex43
 PAR_EXAMPLES = ex0p ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex8p ex9p ex10p ex11p \
    ex12p ex13p ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex21p ex22p ex24p \
    ex25p ex26p ex27p ex28p ex29p ex30p ex31p ex32p ex33p ex34p ex35p ex36p \
-   ex37p ex39p ex40p ex41p
+   ex37p ex39p ex40p ex41p ex43p
 SEQ_DEVICE_EXAMPLES = ex1 ex3 ex4 ex5 ex6 ex9 ex14 ex22 ex24 ex25 ex26 ex34
 PAR_DEVICE_EXAMPLES = ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex9p ex13p ex14p \
    ex22p ex24p ex25p ex26p ex34p ex35p

fem/bilininteg.cpp

@@ -4213,6 +4213,181 @@ void DGElasticityIntegrator::AssembleFaceMatrix(
    }
 }
 
+void SlidingElasticityIntegrator::AssembleFaceMatrix(
+   const FiniteElement &el1, const FiniteElement &el2,
+   FaceElementTransformations &Trans, DenseMatrix &elmat)
+{
+   MFEM_ASSERT(Trans.Elem2No < 0,
+               "support for interior faces is not implemented");
+
+#ifdef MFEM_THREAD_SAFE
+   // For descriptions of these variables, see the class declaration.
+   Vector shape1;
+   DenseMatrix dshape1;
+   DenseMatrix adjJ;
+   DenseMatrix dshape1_ps;
+   Vector nor;
+   Vector nL1;
+   Vector nM1;
+   Vector nt1;
+   Vector dshape1_dnM;
+   Vector dshape1_dnt;
+   DenseMatrix jmat;
+#endif
+
+   const int dim = el1.GetDim();
+   const int ndofs1 = el1.GetDof();
+   const int nvdofs = dim * ndofs1;
+
+   // Initially 'elmat' corresponds to the term:
+   //    < { sigma(u) n . ñ }, v . ñ > =
+   //    < { (lambda div(u) I + mu (grad(u) + grad(u)^T)) n . ñ }, v . ñ >
+   // But eventually, it's going to be replaced by:
+   //    elmat := -elmat + alpha*elmat^T + jmat
+   elmat.SetSize(nvdofs);
+   elmat = 0.;
+
+   const bool kappa_is_nonzero = (kappa != 0.0);
+   if (kappa_is_nonzero)
+   {
+      jmat.SetSize(nvdofs);
+      jmat = 0.;
+   }
+
+   adjJ.SetSize(dim);
+   shape1.SetSize(ndofs1);
+   dshape1.SetSize(ndofs1, dim);
+   dshape1_ps.SetSize(ndofs1, dim);
+   nor.SetSize(dim);
+   nL1.SetSize(dim);
+   nM1.SetSize(dim);
+   nt1.SetSize(dim);
+   dshape1_dnM.SetSize(ndofs1);
+   dshape1_dnt.SetSize(ndofs1);
+
+   const IntegrationRule *ir = IntRule;
+   if (ir == NULL)
+   {
+      // a simple choice for the integration order; is this OK?
+      const int order = 2 * el1.GetOrder();
+      ir = &IntRules.Get(Trans.GetGeometryType(), order);
+   }
+
+   for (int pind = 0; pind < ir->GetNPoints(); ++pind)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(pind);
+
+      // Set the integration point in the face and the neighboring elements
+      Trans.SetAllIntPoints(&ip);
+
+      // Access the neighboring element's integration point
+      const IntegrationPoint &eip1 = Trans.GetElement1IntPoint();
+
+      el1.CalcShape(eip1, shape1);
+      el1.CalcDShape(eip1, dshape1);
+
+      CalcAdjugate(Trans.Elem1->Jacobian(), adjJ);
+      Mult(dshape1, adjJ, dshape1_ps);
+
+      if (dim == 1)
+      {
+         nor(0) = 2*eip1.x - 1.0;
+      }
+      else
+      {
+         CalcOrtho(Trans.Jacobian(), nor);
+      }
+
+      if (!nt)
+      {
+         // Set ñ to the unit normal vector if not provided
+         nt1 = nor;
+         nt1 /= nt1.Norml2();
+      }
+      else
+      {
+         // Evaluate vector function ñ at integration point
+         nt->Eval(nt1, *Trans.Elem1, eip1);
+      }
+
+      const real_t W = ip.weight;
+      const real_t W1 = W / Trans.Elem1->Weight();
+      const real_t WL1 = W1 * lambda->Eval(*Trans.Elem1, eip1);
+      const real_t WM1 = W1 * mu->Eval(*Trans.Elem1, eip1);
+      nL1.Set(WL1, nor);
+      nM1.Set(WM1, nor);
+      const real_t WLM = WL1 + 2.0*WM1;
+      dshape1_ps.Mult(nM1, dshape1_dnM);
+      dshape1_ps.Mult(nt1, dshape1_dnt);
+
+      const real_t jmatcoef = kappa * (nor*nor) * WLM;
+
+      const real_t nL_dot_nt1 = nL1 * nt1;
+      for (int jm = 0, j = 0; jm < dim; ++jm)
+      {
+         for (int jdof = 0; jdof < ndofs1; ++jdof, ++j)
+         {
+            const real_t t1 = dshape1_ps(jdof, jm) * nL_dot_nt1;
+            const real_t t2 = dshape1_dnM(jdof) * nt1(jm);
+            const real_t t3 = dshape1_dnt(jdof) * nM1(jm);
+            const real_t tt = t1 + t2 + t3;
+            for (int im = 0, i = 0; im < dim; ++im)
+            {
+               for (int idof = 0; idof < ndofs1; ++idof, ++i)
+               {
+                  elmat(i, j) += tt * shape1(idof) * nt1(im);
+               }
+            }
+         }
+      }
+
+      if (kappa_is_nonzero)
+      {
+         for (int jm = 0, j = 0; jm < dim; ++jm)
+         {
+            for (int jdof = 0; jdof < ndofs1; ++jdof, ++j)
+            {
+               const real_t sj = jmatcoef * shape1(jdof) * nt1(jm);
+               for (int im = 0, i = 0; im < dim; ++im)
+               {
+                  for (int idof = 0; idof < ndofs1; ++idof, ++i)
+                  {
+                     jmat(i, j) += shape1(idof) * sj * nt1(im);
+                  }
+               }
+            }
+         }
+      }
+   }
+
+   // elmat := -elmat + alpha*elmat^t + jmat
+   if (kappa_is_nonzero)
+   {
+      for (int i = 0; i < nvdofs; ++i)
+      {
+         for (int j = 0; j < i; ++j)
+         {
+            real_t aij = elmat(i,j), aji = elmat(j,i), mij = jmat(i,j);
+            elmat(i,j) = alpha*aji - aij + mij;
+            elmat(j,i) = alpha*aij - aji + mij;
+         }
+         elmat(i,i) = (alpha - 1.)*elmat(i,i) + jmat(i,i);
+      }
+   }
+   else
+   {
+      for (int i = 0; i < nvdofs; ++i)
+      {
+         for (int j = 0; j < i; ++j)
+         {
+            real_t aij = elmat(i,j), aji = elmat(j,i);
+            elmat(i,j) = alpha*aji - aij;
+            elmat(j,i) = alpha*aij - aji;
+         }
+         elmat(i,i) *= (alpha - 1.);
+      }
+   }
+}
 
 void TraceJumpIntegrator::AssembleFaceMatrix(
    const FiniteElement &trial_face_fe, const FiniteElement &test_fe1,

fem/bilininteg.hpp

@@ -3738,6 +3738,84 @@ protected:
       DenseMatrix &elmat, DenseMatrix &jmat);
 };
 
+/** Integrator for the Nitsche elasticity form:
+    $$
+    \begin{split}
+    a(u,v)
+    &:= -\langle \sigma(u)\, \vec{n} \cdot \tilde{n},\ v \cdot \tilde{n}
+    \rangle + \alpha \langle \sigma(v)\, \vec{n} \cdot \tilde{n},\ u \cdot
+    \tilde{n} \rangle + \kappa \langle h^{-1} (\lambda + 2\mu)\, u \cdot
+    \tilde{n},\ v \cdot \tilde{n} \rangle \\
+    &= -\int_\Gamma (\sigma(u)\, n \cdot \tilde{n})(v \cdot \tilde{n})\, dS +
+    \alpha \int_\Gamma (\sigma(v)\, n \cdot \tilde{n})(u \cdot \tilde{n})\,
+    dS + \kappa \int_\Gamma h^{-1} (\lambda + 2\mu)(u \cdot \tilde{n})(v
+    \cdot \tilde{n})\, dS.
+    \end{split}
+    $$
+
+    For isotropic media,
+    $$
+    \begin{split}
+    \sigma(u) &= \lambda \nabla \cdot u I + 2 \mu \varepsilon(u) \\
+              &= \lambda \nabla \cdot u I + 2 \mu \frac{1}{2} (\nabla u + \nabla
+                 u^{\mathrm{T}}) \\
+              &= \lambda \nabla \cdot u I + \mu (\nabla u + \nabla u^{\mathrm{T}})
+    \end{split}
+    $$
+    where $I$ is the identity matrix, $\lambda$ and $\mu$ are the Lamé
+    coefficients (see ElasticityIntegrator), $\tilde{n}$ is a unit vector
+    field, $\alpha = \pm 1$ and $\kappa > 0$ are the Nitsche parameters, and
+    $u$, $v$ are the trial and test functions, respectively.
+
+    This is a '%Vector' integrator, i.e. defined for FE spaces using multiple
+    copies of a scalar FE space.
+ */
+class SlidingElasticityIntegrator : public BilinearFormIntegrator
+{
+public:
+   SlidingElasticityIntegrator(Coefficient &lambda_, Coefficient &mu_,
+                               real_t kappa_)
+      : nt(NULL), lambda(&lambda_), mu(&mu_), alpha(-1.0), kappa(kappa_) { }
+
+   SlidingElasticityIntegrator(VectorCoefficient &nt_, Coefficient &lambda_,
+                               Coefficient &mu_, real_t alpha_, real_t kappa_)
+      : nt(&nt_), lambda(&lambda_), mu(&mu_), alpha(alpha_), kappa(kappa_) { }
+
+   using BilinearFormIntegrator::AssembleFaceMatrix;
+   void AssembleFaceMatrix(const FiniteElement &el1,
+                           const FiniteElement &el2,
+                           FaceElementTransformations &Trans,
+                           DenseMatrix &elmat) override;
+
+protected:
+   VectorCoefficient *nt;
+   Coefficient *lambda, *mu;
+   real_t alpha, kappa;
+
+#ifndef MFEM_THREAD_SAFE
+   // values of all scalar basis functions for one component of u (which is a
+   // vector) at the integration point in the reference space
+   Vector shape1;
+   // values of derivatives of all scalar basis functions for one component
+   // of u (which is a vector) at the integration point in the reference space
+   DenseMatrix dshape1;
+   // Adjugate of the Jacobian of the transformation: adjJ = det(J) J^{-1}
+   DenseMatrix adjJ;
+   // gradient of shape functions in the real (physical, not reference)
+   // coordinates, scaled by det(J):
+   //    dshape_ps(jdof,jm) = sum_{t} adjJ(t,jm)*dshape(jdof,t)
+   DenseMatrix dshape1_ps;
+   Vector nor;  // nor = |weight(J_face)| n
+   Vector nL1;  // nL1 = (lambda1 * ip.weight / detJ1) nor
+   Vector nM1;  // nM1 = (mu1     * ip.weight / detJ1) nor
+   Vector nt1;  // nt1 = vector function ñ evaluated at ip1
+   Vector dshape1_dnM; // dshape1_dnM = dshape1_ps . nM1
+   Vector dshape1_dnt; // dshape1_dnt = dshape1_ps . nt1
+   // 'jmat' corresponds to the term: kappa <h⁻¹ u ⋅ ñ, v ⋅ ñ>
+   DenseMatrix jmat;
+#endif
+};
+
 /** Integrator for the DPG form:$ \langle v, [w] \rangle $ over all faces (the interface) where
     the trial variable $v$ is defined on the interface and the test variable $w$ is
     defined inside the elements, generally in a DG space. */

fem/lininteg.cpp

@@ -1054,7 +1054,154 @@ void DGElasticityDirichletLFIntegrator::AssembleRHSElementVect(
    }
 }
 
+void SlidingElasticityLFIntegrator::AssembleRHSElementVect(
+   const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
+{
+   mfem_error("SlidingElasticityLFIntegrator::AssembleRHSElementVect");
+}
 
+void SlidingElasticityLFIntegrator::AssembleRHSElementVect(
+   const FiniteElement &el, FaceElementTransformations &Tr, Vector &elvect)
+{
+   MFEM_ASSERT(Tr.Elem2No < 0, "interior boundary is not supported");
+
+#ifdef MFEM_THREAD_SAFE
+   Vector shape;
+   DenseMatrix dshape;
+   DenseMatrix adjJ;
+   DenseMatrix dshape_ps;
+   Vector nor;
+   Vector dshape_dn;
+   Vector dshape_du;
+   real_t g_val;
+   Vector nt_val;
+#endif
+
+   const int dim = el.GetDim();
+   const int ndofs = el.GetDof();
+   const int nvdofs = dim*ndofs;
+
+   elvect.SetSize(nvdofs);
+   elvect = 0.0;
+
+   adjJ.SetSize(dim);
+   shape.SetSize(ndofs);
+   dshape.SetSize(ndofs, dim);
+   dshape_ps.SetSize(ndofs, dim);
+   nor.SetSize(dim);
+   dshape_dn.SetSize(ndofs);
+   dshape_du.SetSize(ndofs);
+   nt_val.SetSize(dim);
+
+   const IntegrationRule *ir = IntRule;
+   if (ir == NULL)
+   {
+      const int order = 2*el.GetOrder(); // <-----
+      ir = &IntRules.Get(Tr.GetGeometryType(), order);
+   }
+
+   for (int pi = 0; pi < ir->GetNPoints(); ++pi)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(pi);
+
+      // Set the integration point in the face and the neighboring element
+      Tr.SetAllIntPoints(&ip);
+
+      // Access the neighboring element's integration point
+      const IntegrationPoint &eip = Tr.GetElement1IntPoint();
+
+      el.CalcShape(eip, shape);
+      el.CalcDShape(eip, dshape);
+
+      CalcAdjugate(Tr.Elem1->Jacobian(), adjJ);
+      Mult(dshape, adjJ, dshape_ps);
+
+      if (dim == 1)
+      {
+         nor(0) = 2*eip.x - 1.0;
+      }
+      else
+      {
+         CalcOrtho(Tr.Jacobian(), nor);
+      }
+
+      if (!nt)
+      {
+         // Set nt to the unit normal vector if not provided
+         nt_val = nor;
+         nt_val /= nt_val.Norml2();
+      }
+      else
+      {
+         // Evaluate the vector field using the face transformation.
+         nt->Eval(nt_val, Tr, ip);
+      }
+
+      // Evaluate the Dirichlet b.c. using the face transformation.
+      g_val = g->Eval(Tr, ip);
+
+      real_t WL, WM, jcoef;
+      {
+         const real_t W = ip.weight / Tr.Elem1->Weight();
+         WL = W * lambda->Eval(*Tr.Elem1, eip);
+         WM = W * mu->Eval(*Tr.Elem1, eip);
+         jcoef = kappa * (WL + 2.0*WM) * (nor*nor);
+         dshape_ps.Mult(nor, dshape_dn);
+         dshape_ps.Mult(nt_val, dshape_du);
+      }
+
+      // alpha < g, (lambda div(v) I + mu (grad(v) + grad(v)^T)) n . ñ > +
+      //   + kappa < h^{-1} (lambda + 2 mu) g, v . ñ >
+
+      // i = idof + ndofs * im
+      // v_phi(i,d) = delta(im,d) phi(idof)
+      // div(v_phi(i)) = dphi(idof,im)
+      // (grad(v_phi(i)))(k,l) = delta(im,k) dphi(idof,l)
+      //
+      // term 1:
+      //   alpha < g, lambda div(v_phi(i)) n . ñ > =
+      //   alpha lambda g div(v_phi(i)) (n.ñ) =
+      //   alpha lambda g dphi(idof,im) (n.ñ) --> quadrature -->
+      //   ip.weight/det(J1) alpha lambda g (nor.ñ) dshape_ps(idof,im) =
+      //   alpha * WL * g_val * (nor*nt_val) * dshape_ps(idof,im)
+      // term 2:
+      //   alpha < g, mu grad(v_phi(i)) n . ñ > =
+      //   alpha mu g ñ^T grad(v_phi(i)) n =
+      //   alpha mu g ñ(k) delta(im,k) dphi(idof,l) n(l) =
+      //   alpha mu g ñ(im) dphi(idof,l) n(l) --> quadrature -->
+      //   ip.weight/det(J1) alpha mu ñ(im) g dshape_ps(idof,l) nor(l) =
+      //   alpha * WM * g_val * nt_val(im) * dshape_dn(idof)
+      // term 3:
+      //   alpha < g, mu (grad(v_phi(i)))^T n . ñ > =
+      //   alpha mu g n^T grad(v_phi(i)) ñ =
+      //   alpha mu g n(k) delta(im,k) dphi(idof,l) ñ(l) =
+      //   alpha mu g n(im) dphi(idof,l) ñ(l) --> quadrature -->
+      //   ip.weight/det(J1) alpha mu g nor(im) dshape_ps(idof,l) ñ(l) =
+      //   alpha * WM * g_val * nor(im) * dshape_du(idof)
+      // term j:
+      //   < kappa h^{-1} (lambda + 2 mu) g, ñ . v_phi(i) > =
+      //   kappa/h (lambda + 2 mu) g ñ(k) v_phi(i,k) =
+      //   kappa/h (lambda + 2 mu) g ñ(k) delta(im,k) phi(idof) =
+      //   kappa/h (lambda + 2 mu) g ñ(im) phi(idof) --> quadrature -->
+      //      [ 1/h = |nor|/det(J1) ]
+      //   ip.weight/det(J1) |nor|^2 (lambda + 2 mu) kappa g ñ(im) phi(idof) =
+      //   jcoef * g_val * nt_val(im) * shape(idof)
+
+      WM *= alpha;
+      const real_t t1 = alpha * WL * g_val * (nor*nt_val);
+      for (int im = 0, i = 0; im < dim; ++im)
+      {
+         const real_t t2 = WM * g_val * nt_val(im);
+         const real_t t3 = WM * g_val * nor(im);
+         const real_t tj = jcoef * g_val * nt_val(im);
+         for (int idof = 0; idof < ndofs; ++idof, ++i)
+         {
+            elvect(i) += (t1*dshape_ps(idof,im) + t2*dshape_dn(idof) +
+                          t3*dshape_du(idof) + tj*shape(idof));
+         }
+      }
+   }
+}
 
 void WhiteGaussianNoiseDomainLFIntegrator::AssembleRHSElementVect
 (const FiniteElement &el,

fem/lininteg.hpp

@@ -646,6 +646,62 @@ public:
    using LinearFormIntegrator::AssembleRHSElementVect;
 };
 
+/** Boundary linear form integrator for imposing non-zero Dirichlet boundary
+    conditions, in a Nitsche elasticity formulation. Specifically, the linear
+    form is given by
+    $$
+    \begin{split}
+    b(v) &:= \alpha \int_\Gamma (\lambda\, \mathrm{div}(v)\, I + \mu (\nabla v
+    + \nabla v^{\mathrm{T}}))\, n \cdot \tilde{n}\, g\, dS + \kappa \int_\Gamma
+    h^{-1} (\lambda + 2\mu) (v \cdot \tilde{n})\, g\, dS
+    \end{split}
+    $$
+    where $g$ is the given Dirichlet data, $n$ is the unit normal, $\tilde{n}$ is
+    a unit vector field, and $\alpha = \pm 1$, $\kappa > 0$ are the Nitsche
+    parameters. The parameters $\lambda$ and $\mu$ should match the parameters
+    with the same names used in the bilinear form integrator,
+    SlidingElasticityIntegrator.
+ */
+class SlidingElasticityLFIntegrator : public LinearFormIntegrator
+{
+protected:
+   Coefficient *g;
+   VectorCoefficient *nt;
+   Coefficient *lambda, *mu;
+   real_t alpha, kappa;
+
+#ifndef MFEM_THREAD_SAFE
+   Vector shape;
+   DenseMatrix dshape;
+   DenseMatrix adjJ;
+   DenseMatrix dshape_ps;
+   Vector nor;
+   Vector dshape_dn;
+   Vector dshape_du;
+   real_t g_val;
+   Vector nt_val;
+#endif
+
+public:
+   SlidingElasticityLFIntegrator(Coefficient &g_,
+                                 Coefficient &lambda_, Coefficient &mu_,
+                                 real_t kappa_)
+      : g(&g_), nt(NULL), lambda(&lambda_), mu(&mu_), alpha(-1.0), kappa(kappa_) {}
+
+   SlidingElasticityLFIntegrator(Coefficient &g_, VectorCoefficient &nt_,
+                                 Coefficient &lambda_, Coefficient &mu_,
+                                 real_t alpha_, real_t kappa_)
+      : g(&g_), nt(&nt_), lambda(&lambda_), mu(&mu_), alpha(alpha_), kappa(kappa_) {}
+
+   void AssembleRHSElementVect(const FiniteElement &el,
+                               ElementTransformation &Tr,
+                               Vector &elvect) override;
+   void AssembleRHSElementVect(const FiniteElement &el,
+                               FaceElementTransformations &Tr,
+                               Vector &elvect) override;
+
+   using LinearFormIntegrator::AssembleRHSElementVect;
+};
 
 /** Class for spatial white Gaussian noise integration.