Merge pull request sandialabs#67 from tupek2/main

Improved quadratic patch test

optimism/test/test_PatchTest.py

@@ -13,15 +13,18 @@
 from optimism import QuadratureRule
 from optimism.test import MeshFixture
 
+kappa = 1.0
+mu = 0.25
+
+E = 9 * kappa * mu / (3 * kappa + mu)
+nu = 0.5 * (3 * kappa - 2 * mu) / (3 * kappa + mu)
 
-E = 1.0
-nu = 0.3
 props = {'elastic modulus': E,
          'poisson ratio': nu,
          'strain measure': 'linear'}
 
 
-class PatchTest(MeshFixture.MeshFixture):
+class LinearPatchTestLinearElements(MeshFixture.MeshFixture):
 
     def setUp(self):
         self.Nx = 7
@@ -108,7 +111,7 @@ def objective(Uu):
 
 
 
-class PatchTestQuadraticElements(MeshFixture.MeshFixture):
+class LinearPatchTestQuadraticElements(MeshFixture.MeshFixture):
 
     def setUp(self):
         self.Nx = 3
@@ -196,6 +199,66 @@ def objective(Uu):
         Uu = dofManager.get_unknown_values(U)
         self.assertNear(np.linalg.norm(grad_func(Uu)), 0.0, 14)
 
+
+class QuadraticPatchTestQuadraticElements(MeshFixture.MeshFixture):
+
+    def setUp(self):
+        self.Nx = 3
+        self.Ny = 3
+        xRange = [0.,1.]
+        yRange = [0.,1.]
+
+        self.targetDispGrad = np.array([[0.1, -0.2],[0.4, -0.1]]) 
+
+        mesh, _ = self.create_mesh_and_disp(self.Nx, self.Ny, xRange, yRange,
+                                            lambda x : self.targetDispGrad.dot(x))
+
+        self.mesh = Mesh.create_higher_order_mesh_from_simplex_mesh(mesh, order=2, createNodeSetsFromSideSets=True)
+
+        alpha = 2.0
+        beta = 1.0
+
+        self.UTarget = jax.vmap( lambda x : self.targetDispGrad@x + 
+                                 np.array([alpha * x[0]*x[0], beta * x[1]*x[1]]) )(self.mesh.coords)
+
+        self.quadRule = QuadratureRule.create_quadrature_rule_on_triangle(degree=2)
+        self.fs = FunctionSpace.construct_function_space(self.mesh, self.quadRule)
+
+        materialModel = MatModel.create_material_model_functions(props)
+
+        self.b = -(2*kappa + 8*mu/3.0) * np.array([alpha, beta])
+
+        mcxFuncs = Mechanics.create_mechanics_functions(self.fs, "plane strain", materialModel)
+
+        self.compute_energy = jax.jit(mcxFuncs.compute_strain_energy)
+        self.internals = mcxFuncs.compute_initial_state()
+
+    def test_dirichlet_patch_test_with_quadratic_elements(self):
+        ebcs = [FunctionSpace.EssentialBC(nodeSet='all_boundary', component=0),
+                FunctionSpace.EssentialBC(nodeSet='all_boundary', component=1)]
+        self.dofManager = FunctionSpace.DofManager(self.fs, self.UTarget.shape[1], ebcs)
+        Ubc = self.dofManager.get_bc_values(self.UTarget)
+
+        def constant_body_force_potential(U, internals, b):
+            dtUnused = 0.0
+            f = lambda u, du, q, x, dt: -np.dot(b, u)
+            return FunctionSpace.integrate_over_block(self.fs, U, internals, dtUnused, f, slice(None))
+
+        def objective(Uu):
+            U = self.dofManager.create_field(Uu, Ubc)
+            return self.compute_energy(U, self.internals) + constant_body_force_potential(U, self.internals, self.b)
+
+        with Timer(name="NewtonSolve"):
+            Uu = newton_solve(objective, 0.0*self.dofManager.get_unknown_values(self.UTarget))
+
+        U = self.dofManager.create_field(Uu, Ubc)
+
+        self.assertArrayNear(U, self.UTarget, 13)
+
+        grad_func = jax.jit(jax.grad(objective))
+        Uu = self.dofManager.get_unknown_values(U)
+        self.assertNear(np.linalg.norm(grad_func(Uu)), 0.0, 14)
+
 
 if __name__ == '__main__':
     unittest.main()