Re-activate commented lines in hyperelastic demo (#3146)

* Re-activate commented lines in hyperelastic demo

Resolves #2584

* Updates

* Fix constants

cpp/demo/hyperelasticity/hyperelasticity.py

@@ -7,6 +7,7 @@
 from basix.ufl import element
 from ufl import (
     Coefficient,
+    Constant,
     FunctionSpace,
     Identity,
     Mesh,
@@ -15,8 +16,10 @@
     derivative,
     det,
     diff,
+    ds,
     dx,
     grad,
+    inner,
     ln,
     tr,
     variable,
@@ -39,8 +42,8 @@
 
 # Functions
 u = Coefficient(V)  # Displacement from previous iteration
-# B = Coefficient(element)        # Body force per unit volume
-# T = Coefficient(element)        # Traction force on the boundary
+B = Constant(mesh, shape=(3,))  # Body force per unit volume
+T = Constant(mesh, shape=(3,))  # Traction force on the boundary
 
 # Now, we can define the kinematic quantities involved in the model:
 
@@ -72,7 +75,7 @@
 psi = (mu / 2) * (Ic - 3) - mu * ln(J) + (lmbda / 2) * (ln(J)) ** 2
 
 # Total potential energy
-Pi = psi * dx  # - inner(B, u) * dx - inner(T, u) * ds
+Pi = psi * dx - inner(B, u) * dx - inner(T, u) * ds
 
 # First variation of Pi (directional derivative about u in the direction
 # of v)

cpp/demo/hyperelasticity/main.cpp

@@ -156,12 +156,17 @@ int main(int argc, char* argv[])
     auto V = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace(
         functionspace_form_hyperelasticity_F_form, "u", mesh));
 
+    auto B = std::make_shared<fem::Constant<T>>(std::vector<T>{0, 0, 0});
+    auto traction = std::make_shared<fem::Constant<T>>(std::vector<T>{0, 0, 0});
+
     // Define solution function
     auto u = std::make_shared<fem::Function<T>>(V);
-    auto a = std::make_shared<fem::Form<T>>(fem::create_form<T>(
-        *form_hyperelasticity_J_form, {V, V}, {{"u", u}}, {}, {}));
-    auto L = std::make_shared<fem::Form<T>>(fem::create_form<T>(
-        *form_hyperelasticity_F_form, {V}, {{"u", u}}, {}, {}));
+    auto a = std::make_shared<fem::Form<T>>(
+        fem::create_form<T>(*form_hyperelasticity_J_form, {V, V}, {{"u", u}},
+                            {{"B", B}, {"T", traction}}, {}));
+    auto L = std::make_shared<fem::Form<T>>(
+        fem::create_form<T>(*form_hyperelasticity_F_form, {V}, {{"u", u}},
+                            {{"B", B}, {"T", traction}}, {}));
 
     auto u_rotation = std::make_shared<fem::Function<T>>(V);
     u_rotation->interpolate(
@@ -239,7 +244,8 @@ int main(int argc, char* argv[])
     newton_solver.atol = 10 * std::numeric_limits<T>::epsilon();
 
     la::petsc::Vector _u(la::petsc::create_vector_wrap(*u->x()), false);
-    newton_solver.solve(_u.vec());
+    auto [niter, success] = newton_solver.solve(_u.vec());
+    std::cout << "Number of Newton iterations: " << niter << std::endl;
 
     // Compute Cauchy stress. Construct appropriate Basix element for
     // stress.