diff --git a/optimism/Interpolants.py b/optimism/Interpolants.py index 395f289..75ed693 100644 --- a/optimism/Interpolants.py +++ b/optimism/Interpolants.py @@ -69,7 +69,7 @@ def get_lobatto_nodes_1d(degree): p = onp.polynomial.Legendre.basis(degree, domain=[0.0, 1.0]) dp = p.deriv() xInterior = dp.roots() - xn = np.hstack(([0.0], xInterior, [1.0])) + xn = np.hstack((np.array([0.0]), xInterior, np.array([1.0]))) return xn diff --git a/optimism/inverse/AdjointFunctionSpace.py b/optimism/inverse/AdjointFunctionSpace.py new file mode 100644 index 0000000..8c6ac0c --- /dev/null +++ b/optimism/inverse/AdjointFunctionSpace.py @@ -0,0 +1,30 @@ +from optimism import FunctionSpace +from optimism import Interpolants +from optimism import Mesh +from optimism.FunctionSpace import compute_element_volumes +from optimism.FunctionSpace import compute_element_volumes_axisymmetric +from optimism.FunctionSpace import map_element_shape_grads +from jax import vmap + +def construct_function_space_for_adjoint(coords, mesh, quadratureRule, mode2D='cartesian'): + + shapeOnRef = Interpolants.compute_shapes(mesh.parentElement, quadratureRule.xigauss) + + shapes = vmap(lambda elConns, elShape: elShape, (0, None))(mesh.conns, shapeOnRef.values) + + shapeGrads = vmap(map_element_shape_grads, (None, 0, None, None))(coords, mesh.conns, mesh.parentElement, shapeOnRef.gradients) + + if mode2D == 'cartesian': + el_vols = compute_element_volumes + isAxisymmetric = False + elif mode2D == 'axisymmetric': + el_vols = compute_element_volumes_axisymmetric + isAxisymmetric = True + vols = vmap(el_vols, (None, 0, None, 0, None))(coords, mesh.conns, mesh.parentElement, shapes, quadratureRule.wgauss) + + # unpack mesh and remake a mesh to make sure we get all the AD + mesh = Mesh.Mesh(coords=coords, conns=mesh.conns, simplexNodesOrdinals=mesh.simplexNodesOrdinals, + parentElement=mesh.parentElement, parentElement1d=mesh.parentElement1d, blocks=mesh.blocks, + nodeSets=mesh.nodeSets, sideSets=mesh.sideSets) + + return FunctionSpace.FunctionSpace(shapes, vols, shapeGrads, mesh, quadratureRule, isAxisymmetric) diff --git a/optimism/inverse/MechanicsInverse.py b/optimism/inverse/MechanicsInverse.py new file mode 100644 index 0000000..9b85f1c --- /dev/null +++ b/optimism/inverse/MechanicsInverse.py @@ -0,0 +1,99 @@ +from collections import namedtuple + +from optimism.JaxConfig import * +from optimism import FunctionSpace +from optimism import Interpolants +from optimism.TensorMath import tensor_2D_to_3D + +IvsUpdateInverseFunctions = namedtuple('IvsUpdateInverseFunctions', + ['ivs_update_jac_ivs_prev', + 'ivs_update_jac_disp_vjp', + 'ivs_update_jac_coords_vjp']) + +PathDependentResidualInverseFunctions = namedtuple('PathDependentResidualInverseFunctions', + ['residual_jac_ivs_prev_vjp', + 'residual_jac_coords_vjp']) + +ResidualInverseFunctions = namedtuple('ResidualInverseFunctions', + ['residual_jac_coords_vjp']) + +def _compute_element_field_gradient(U, elemShapeGrads, elemConnectivity, modify_element_gradient): + elemNodalDisps = U[elemConnectivity] + elemGrads = vmap(FunctionSpace.compute_quadrature_point_field_gradient, (None, 0))(elemNodalDisps, elemShapeGrads) + elemGrads = modify_element_gradient(elemGrads) + return elemGrads + +def _compute_field_gradient(shapeGrads, conns, nodalField, modify_element_gradient): + return vmap(_compute_element_field_gradient, (None,0,0,None))(nodalField, shapeGrads, conns, modify_element_gradient) + +def _compute_updated_internal_variables_gradient(dispGrads, states, dt, compute_state_new, output_shape): + dgQuadPointRavel = dispGrads.reshape(dispGrads.shape[0]*dispGrads.shape[1],*dispGrads.shape[2:]) + stQuadPointRavel = states.reshape(states.shape[0]*states.shape[1],*states.shape[2:]) + statesNew = vmap(compute_state_new, (0, 0, None))(dgQuadPointRavel, stQuadPointRavel, dt) + return statesNew.reshape(output_shape) + + +def create_ivs_update_inverse_functions(functionSpace, mode2D, materialModel, pressureProjectionDegree=None): + fs = functionSpace + shapeOnRef = Interpolants.compute_shapes(fs.mesh.parentElement, fs.quadratureRule.xigauss) + + if mode2D == 'plane strain': + grad_2D_to_3D = vmap(tensor_2D_to_3D) + elif mode2D == 'axisymmetric': + raise NotImplementedError + + modify_element_gradient = grad_2D_to_3D + if pressureProjectionDegree is not None: + raise NotImplementedError + + def compute_partial_ivs_update_partial_ivs_prev(U, stateVariables, dt=0.0): + dispGrads = _compute_field_gradient(fs.shapeGrads, fs.mesh.conns, U, modify_element_gradient) + update_gradient = jacfwd(materialModel.compute_state_new, argnums=1) + grad_shape = stateVariables.shape + (stateVariables.shape[2],) + return _compute_updated_internal_variables_gradient(dispGrads, stateVariables, dt,\ + update_gradient, grad_shape) + + def compute_ivs_update_parameterized(U, stateVariables, coords, dt=0.0): + shapeGrads = vmap(FunctionSpace.map_element_shape_grads, (None, 0, None, None))(coords, fs.mesh.conns, fs.mesh.parentElement, shapeOnRef.gradients) + dispGrads = _compute_field_gradient(shapeGrads, fs.mesh.conns, U, modify_element_gradient) + update_func = materialModel.compute_state_new + output_shape = stateVariables.shape + return _compute_updated_internal_variables_gradient(dispGrads, stateVariables, dt,\ + update_func, output_shape) + + compute_partial_ivs_update_partial_coords = jit(lambda u, ivs, x, av, dt=0.0: + vjp(lambda z: compute_ivs_update_parameterized(u, ivs, z, dt), x)[1](av)[0]) + + def compute_ivs_update(U, stateVariables, dt=0.0): + dispGrads = _compute_field_gradient(fs.shapeGrads, fs.mesh.conns, U, modify_element_gradient) + update_func = materialModel.compute_state_new + output_shape = stateVariables.shape + return _compute_updated_internal_variables_gradient(dispGrads, stateVariables, dt,\ + update_func, output_shape) + + compute_partial_ivs_update_partial_disp = jit(lambda x, ivs, av, dt=0.0: + vjp(lambda z: compute_ivs_update(z, ivs, dt), x)[1](av)[0]) + + return IvsUpdateInverseFunctions(jit(compute_partial_ivs_update_partial_ivs_prev), + compute_partial_ivs_update_partial_disp, + compute_partial_ivs_update_partial_coords + ) + +def create_path_dependent_residual_inverse_functions(energyFunction): + + compute_partial_residual_partial_ivs_prev = jit(lambda u, q, iv, x, vx: + vjp(lambda z: grad(energyFunction, 0)(u, q, z, x), iv)[1](vx)[0]) + + compute_partial_residual_partial_coords = jit(lambda u, q, iv, x, vx: + vjp(lambda z: grad(energyFunction, 0)(u, q, iv, z), x)[1](vx)[0]) + + return PathDependentResidualInverseFunctions(compute_partial_residual_partial_ivs_prev, + compute_partial_residual_partial_coords + ) + +def create_residual_inverse_functions(energyFunction): + + compute_partial_residual_partial_coords = jit(lambda u, q, x, vx: + vjp(lambda z: grad(energyFunction, 0)(u, q, z), x)[1](vx)[0]) + + return ResidualInverseFunctions(compute_partial_residual_partial_coords) \ No newline at end of file diff --git a/optimism/inverse/test/FiniteDifferenceFixture.py b/optimism/inverse/test/FiniteDifferenceFixture.py new file mode 100644 index 0000000..c2dece5 --- /dev/null +++ b/optimism/inverse/test/FiniteDifferenceFixture.py @@ -0,0 +1,67 @@ +from optimism.test.MeshFixture import MeshFixture +from collections import namedtuple +import numpy as onp + +class FiniteDifferenceFixture(MeshFixture): + def assertFiniteDifferenceCheckHasVShape(self, errors, tolerance=1e-6): + minError = min(errors) + self.assertLess(minError, tolerance, "Smallest finite difference error not less than tolerance.") + self.assertLess(minError, errors[0], "Finite difference error does not decrease from initial step size.") + self.assertLess(minError, errors[-1], "Finite difference error does not increase after reaching minimum. Try more finite difference steps.") + + def build_direction_vector(self, numDesignVars, seed=123): + + onp.random.seed(seed) + directionVector = onp.random.uniform(-1.0, 1.0, numDesignVars) + normVector = directionVector / onp.linalg.norm(directionVector) + + return onp.array(normVector) + + def compute_finite_difference_error(self, stepSize, initialParameters): + storedState = self.forward_solve(initialParameters) + originalObjective = self.compute_objective_function(storedState, initialParameters) + gradient = self.compute_gradient(storedState, initialParameters) + + directionVector = self.build_direction_vector(initialParameters.shape[0]) + directionalDerivative = onp.tensordot(directionVector, gradient, axes=1) + + perturbedParameters = initialParameters + stepSize * directionVector + storedState = self.forward_solve(perturbedParameters) + perturbedObjective = self.compute_objective_function(storedState, perturbedParameters) + + fd_value = (perturbedObjective - originalObjective) / stepSize + error = abs(directionalDerivative - fd_value) + + return error + + def compute_finite_difference_errors(self, stepSize, steps, initialParameters, printOutput=True): + storedState = self.forward_solve(initialParameters) + originalObjective = self.compute_objective_function(storedState, initialParameters) + gradient = self.compute_gradient(storedState, initialParameters) + + directionVector = self.build_direction_vector(initialParameters.shape[0]) + directionalDerivative = onp.tensordot(directionVector, gradient, axes=1) + + fd_values = [] + errors = [] + for i in range(0, steps): + perturbedParameters = initialParameters + stepSize * directionVector + storedState = self.forward_solve(perturbedParameters) + perturbedObjective = self.compute_objective_function(storedState, perturbedParameters) + + fd_value = (perturbedObjective - originalObjective) / stepSize + fd_values.append(fd_value) + + error = abs(directionalDerivative - fd_value) + errors.append(error) + + stepSize *= 1e-1 + + if printOutput: + print("\n grad'*dir | FD approx | abs error") + print("--------------------------------------------------------------------------------") + for i in range(0, steps): + print(f" {directionalDerivative} | {fd_values[i]} | {errors[i]}") + + return errors + diff --git a/optimism/inverse/test/test_Hyperelastic_gradient_checks.py b/optimism/inverse/test/test_Hyperelastic_gradient_checks.py new file mode 100644 index 0000000..593051d --- /dev/null +++ b/optimism/inverse/test/test_Hyperelastic_gradient_checks.py @@ -0,0 +1,216 @@ +import unittest + +import jax +import jax.numpy as np +import numpy as onp +from scipy.sparse import linalg + +from optimism import EquationSolver as EqSolver +from optimism import QuadratureRule +from optimism import FunctionSpace +from optimism import Mechanics +from optimism import Objective +from optimism import Mesh +from optimism.material import Neohookean + +from optimism.inverse.test.FiniteDifferenceFixture import FiniteDifferenceFixture + +from optimism.inverse import MechanicsInverse +from optimism.inverse import AdjointFunctionSpace + +class NeoHookeanGlobalMeshAdjointSolveFixture(FiniteDifferenceFixture): + def setUp(self): + dispGrad0 = np.array([[0.4, -0.2], + [-0.04, 0.68]]) + self.initialMesh, self.U = self.create_mesh_and_disp(4,4,[0.,1.],[0.,1.], + lambda x: dispGrad0@x) + + shearModulus = 0.855 # MPa + bulkModulus = 100*shearModulus # MPa + youngModulus = 9.0*bulkModulus*shearModulus / (3.0*bulkModulus + shearModulus) + poissonRatio = (3.0*bulkModulus - 2.0*shearModulus) / 2.0 / (3.0*bulkModulus + shearModulus) + props = { + 'elastic modulus': youngModulus, + 'poisson ratio': poissonRatio, + 'version': 'coupled' + } + self.materialModel = Neohookean.create_material_model_functions(props) + + self.quadRule = QuadratureRule.create_quadrature_rule_on_triangle(degree=1) + + self.EBCs = [FunctionSpace.EssentialBC(nodeSet='all_boundary', component=0), + FunctionSpace.EssentialBC(nodeSet='all_boundary', component=1)] + + self.steps = 2 + + def forward_solve(self, parameters): + + self.mesh = Mesh.construct_mesh_from_basic_data(parameters.reshape(self.initialMesh.coords.shape),\ + self.initialMesh.conns, self.initialMesh.blocks,\ + self.initialMesh.nodeSets, self.initialMesh.sideSets) + + functionSpace = FunctionSpace.construct_function_space(self.mesh, self.quadRule) + mechFuncs = Mechanics.create_mechanics_functions(functionSpace, + "plane strain", + self.materialModel) + self.dofManager = FunctionSpace.DofManager(functionSpace, 2, self.EBCs) + Ubc = self.dofManager.get_bc_values(self.U) + + Ubc_inc = Ubc / self.steps + ivs = mechFuncs.compute_initial_state() + p = Objective.Params(bc_data=np.zeros(Ubc.shape), state_data=ivs) + + def compute_energy(Uu, p): + U = self.dofManager.create_field(Uu, p.bc_data) + internalVariables = p.state_data + return mechFuncs.compute_strain_energy(U, internalVariables) + + Uu = 0.0*self.dofManager.get_unknown_values(self.U) + self.objective = Objective.Objective(compute_energy, Uu, p) + + storedState = [] + storedState.append((Uu, p)) + + for step in range(1, self.steps+1): + p = Objective.param_index_update(p, 0, step*Ubc_inc) + Uu = EqSolver.nonlinear_equation_solve(self.objective, Uu, p, EqSolver.get_settings(), useWarmStart=False) + storedState.append((Uu, p)) + + return storedState + + def strain_energy_objective(self, storedState, parameters): + coords = parameters.reshape(self.mesh.coords.shape) + adjoint_func_space = AdjointFunctionSpace.construct_function_space_for_adjoint(coords, self.mesh, self.quadRule) + mech_funcs = Mechanics.create_mechanics_functions(adjoint_func_space, mode2D='plane strain', materialModel=self.materialModel) + + def energy_function_coords(Uu, p, coords): + U = self.dofManager.create_field(Uu, p.bc_data) + internal_vars = p[1] + return mech_funcs.compute_strain_energy(U, internal_vars) + + return energy_function_coords(storedState[-1][0], storedState[-1][1], parameters) + + def strain_energy_gradient(self, storedState, parameters): + return jax.grad(self.strain_energy_objective, argnums=1)(storedState, parameters) + + def total_work_increment(self, Uu, Uu_prev, p, p_prev, coords, nodal_forces): + index = (self.mesh.nodeSets['left'], 1) # arbitrarily choosing left side nodeset for reaction force + + U = self.dofManager.create_field(Uu, p.bc_data) + force = np.array(nodal_forces(U, p, coords).at[index].get()) + disp = U.at[index].get() + + U_prev = self.dofManager.create_field(Uu_prev, p_prev.bc_data) + force_prev = np.array(nodal_forces(U_prev, p_prev, coords).at[index].get()) + disp_prev = U_prev.at[index].get() + + return 0.5*np.tensordot((force + force_prev),(disp - disp_prev), axes=1) + + def total_work_objective(self, storedState, parameters): + parameters = parameters.reshape(self.mesh.coords.shape) + + def energy_function_all_dofs(U, p, coords): + adjoint_func_space = AdjointFunctionSpace.construct_function_space_for_adjoint(coords, self.mesh, self.quadRule) + mech_funcs = Mechanics.create_mechanics_functions(adjoint_func_space, mode2D='plane strain',materialModel=self.materialModel) + internal_vars = p[1] + return mech_funcs.compute_strain_energy(U, internal_vars) + + nodal_forces = jax.jit(jax.grad(energy_function_all_dofs, argnums=0)) + + val = 0.0 + for step in range(1, self.steps+1): + Uu = storedState[step][0] + Uu_prev = storedState[step-1][0] + p = storedState[step][1] + p_prev = storedState[step-1][1] + val += self.total_work_increment(Uu, Uu_prev, p, p_prev, parameters, nodal_forces) + + return val + + def total_work_gradient_just_jax(self, storedState, parameters): + return jax.grad(self.total_work_objective, argnums=1)(storedState, parameters) + + def total_work_gradient_with_adjoint(self, storedState, parameters): + def energy_function_coords(Uu, p, coords): + adjoint_func_space = AdjointFunctionSpace.construct_function_space_for_adjoint(coords, self.mesh, self.quadRule) + mech_funcs = Mechanics.create_mechanics_functions(adjoint_func_space, mode2D='plane strain', materialModel=self.materialModel) + U = self.dofManager.create_field(Uu, p.bc_data) + internal_vars = p[1] + return mech_funcs.compute_strain_energy(U, internal_vars) + + def energy_function_all_dofs(U, p, coords): + adjoint_func_space = AdjointFunctionSpace.construct_function_space_for_adjoint(coords, self.mesh, self.quadRule) + mech_funcs = Mechanics.create_mechanics_functions(adjoint_func_space, mode2D='plane strain',materialModel=self.materialModel) + internal_vars = p[1] + return mech_funcs.compute_strain_energy(U, internal_vars) + + nodal_forces = jax.jit(jax.grad(energy_function_all_dofs, argnums=0)) + + residualInverseFuncs = MechanicsInverse.create_residual_inverse_functions(energy_function_coords) + + compute_df = jax.grad(self.total_work_increment, (0, 1, 4)) + + parameters = parameters.reshape(self.mesh.coords.shape) + gradient = np.zeros(parameters.shape) + adjointLoad = np.zeros(storedState[0][0].shape) + + for step in reversed(range(1, self.steps+1)): + Uu = storedState[step][0] + p = storedState[step][1] + + Uu_prev = storedState[step-1][0] + p_prev = storedState[step-1][1] + + df_du, df_dun, df_dx = compute_df(Uu, Uu_prev, p, p_prev, parameters, nodal_forces) + + adjointLoad -= df_du + + n = self.dofManager.get_unknown_size() + self.objective.p = p # have to update parameters to get precond to work + self.objective.update_precond(Uu) # update preconditioner for use in cg (will converge in 1 iteration as long as the preconditioner is not approximate) + dRdu = linalg.LinearOperator((n, n), lambda V: onp.asarray(self.objective.hessian_vec(Uu, V))) + dRdu_decomp = linalg.LinearOperator((n, n), lambda V: onp.asarray(self.objective.apply_precond(V))) + adjointVector = linalg.cg(dRdu, onp.array(adjointLoad, copy=False), tol=1e-10, atol=0.0, M=dRdu_decomp)[0] + + gradient += df_dx + gradient += residualInverseFuncs.residual_jac_coords_vjp(Uu, p, parameters, adjointVector) + + adjointLoad = -df_dun + # The action dRdU * lambda (the same as lambda^T * dRdU) - For Hyperelastic the residual is not dependent on Uu_prev, so don't actually need this term + + return gradient.ravel() + + + + def test_self_adjoint_gradient(self): + self.compute_objective_function = self.strain_energy_objective + self.compute_gradient = self.strain_energy_gradient + + stepSize = 1e-7 + + error = self.compute_finite_difference_error(stepSize, self.initialMesh.coords.ravel()) + self.assertLessEqual(error, 1e-7) + + @unittest.expectedFailure + def test_non_self_adjoint_gradient_without_adjoint_solve(self): + self.compute_objective_function = self.total_work_objective + self.compute_gradient = self.total_work_gradient_just_jax + + stepSize = 1e-7 + + error = self.compute_finite_difference_error(stepSize, self.initialMesh.coords.ravel()) + self.assertLessEqual(error, 1e-7) + + def test_non_self_adjoint_gradient_with_adjoint_solve(self): + self.compute_objective_function = self.total_work_objective + self.compute_gradient = self.total_work_gradient_with_adjoint + + stepSize = 1e-7 + + error = self.compute_finite_difference_error(stepSize, self.initialMesh.coords.ravel()) + self.assertLessEqual(error, 1e-7) + + + +if __name__ == '__main__': + unittest.main() \ No newline at end of file diff --git a/optimism/inverse/test/test_J2Plastic_gradient_checks.py b/optimism/inverse/test/test_J2Plastic_gradient_checks.py new file mode 100644 index 0000000..63f116d --- /dev/null +++ b/optimism/inverse/test/test_J2Plastic_gradient_checks.py @@ -0,0 +1,310 @@ +import unittest + +import jax +import jax.numpy as np +import numpy as onp +from scipy.sparse import linalg + +from optimism import EquationSolver as EqSolver +from optimism import QuadratureRule +from optimism import FunctionSpace +from optimism import Mechanics +from optimism import Objective +from optimism import Mesh +from optimism.material import J2Plastic as J2 + +from optimism.inverse.test.FiniteDifferenceFixture import FiniteDifferenceFixture + +from optimism.inverse import MechanicsInverse +from optimism.inverse import AdjointFunctionSpace + +class J2GlobalMeshAdjointSolveFixture(FiniteDifferenceFixture): + def setUp(self): + dispGrad0 = np.array([[0.4, -0.2], + [-0.04, 0.68]]) + self.initialMesh, self.U = self.create_mesh_and_disp(4,4,[0.,1.],[0.,1.], + lambda x: dispGrad0@x) + + E = 100.0 + poisson = 0.321 + H = 1e-2 * E + Y0 = 0.3 * E + props = { + 'elastic modulus': E, + 'poisson ratio': poisson, + 'yield strength': Y0, + 'kinematics': 'small deformations', + 'hardening model': 'linear', + 'hardening modulus': H + } + self.materialModel = J2.create_material_model_functions(props) + + self.quadRule = QuadratureRule.create_quadrature_rule_on_triangle(degree=1) + + self.EBCs = [FunctionSpace.EssentialBC(nodeSet='all_boundary', component=0), + FunctionSpace.EssentialBC(nodeSet='all_boundary', component=1)] + + self.steps = 2 + + def forward_solve(self, parameters): + self.mesh = Mesh.construct_mesh_from_basic_data(parameters.reshape(self.initialMesh.coords.shape),\ + self.initialMesh.conns, self.initialMesh.blocks,\ + self.initialMesh.nodeSets, self.initialMesh.sideSets) + + functionSpace = FunctionSpace.construct_function_space(self.mesh, self.quadRule) + mechFuncs = Mechanics.create_mechanics_functions(functionSpace, + "plane strain", + self.materialModel) + self.dofManager = FunctionSpace.DofManager(functionSpace, 2, self.EBCs) + Ubc = self.dofManager.get_bc_values(self.U) + + Ubc_inc = Ubc / self.steps + ivs = mechFuncs.compute_initial_state() + p = Objective.Params(bc_data=np.zeros(Ubc.shape), state_data=ivs) + + def compute_energy(Uu, p): + U = self.dofManager.create_field(Uu, p.bc_data) + internalVariables = p.state_data + return mechFuncs.compute_strain_energy(U, internalVariables) + + Uu = 0.0*self.dofManager.get_unknown_values(self.U) + self.objective = Objective.Objective(compute_energy, Uu, p) + + storedState = [] + storedState.append((Uu, p)) + + for step in range(1, self.steps+1): + p = Objective.param_index_update(p, 0, step*Ubc_inc) + Uu = EqSolver.nonlinear_equation_solve(self.objective, Uu, p, EqSolver.get_settings(), useWarmStart=False) + U = self.dofManager.create_field(Uu, p.bc_data) + ivs = mechFuncs.compute_updated_internal_variables(U, p.state_data) + p = Objective.param_index_update(p, 1, ivs) + storedState.append((Uu, p)) + + return storedState + + def total_work_increment(self, Uu, Uu_prev, ivs, ivs_prev, p, p_prev, coordinates, nodal_forces): + index = (self.mesh.nodeSets['left'], 1) # arbitrarily choosing left side nodeset for reaction force + + U = self.dofManager.create_field(Uu, p.bc_data) + force = np.array(nodal_forces(U, ivs, coordinates).at[index].get()) + disp = U.at[index].get() + + U_prev = self.dofManager.create_field(Uu_prev, p_prev.bc_data) + force_prev = np.array(nodal_forces(U_prev, ivs_prev, coordinates).at[index].get()) + disp_prev = U_prev.at[index].get() + + return 0.5*np.tensordot((force + force_prev),(disp - disp_prev), axes=1) + + def total_work_objective(self, storedState, parameters): + parameters = parameters.reshape(self.mesh.coords.shape) + + def energy_function_all_dofs(U, ivs, coords): + adjoint_func_space = AdjointFunctionSpace.construct_function_space_for_adjoint(coords, self.mesh, self.quadRule) + mech_funcs = Mechanics.create_mechanics_functions(adjoint_func_space, mode2D='plane strain',materialModel=self.materialModel) + return mech_funcs.compute_strain_energy(U, ivs) + + nodal_forces = jax.jit(jax.grad(energy_function_all_dofs, argnums=0)) + + val = 0.0 + for step in range(1, self.steps+1): + Uu = storedState[step][0] + Uu_prev = storedState[step-1][0] + p = storedState[step][1] + p_prev = storedState[step-1][1] + ivs = p.state_data + ivs_prev = p_prev.state_data + + val += self.total_work_increment(Uu, Uu_prev, ivs, ivs_prev, p, p_prev, parameters, nodal_forces) + + return val + + def total_work_gradient(self, storedState, parameters): + + def energy_function_coords(Uu, p, ivs_prev, coords): + adjoint_func_space = AdjointFunctionSpace.construct_function_space_for_adjoint(coords, self.mesh, self.quadRule) + mech_funcs = Mechanics.create_mechanics_functions(adjoint_func_space, mode2D='plane strain', materialModel=self.materialModel) + U = self.dofManager.create_field(Uu, p.bc_data) + return mech_funcs.compute_strain_energy(U, ivs_prev) + + def energy_function_all_dofs(U, ivs, coords): + adjoint_func_space = AdjointFunctionSpace.construct_function_space_for_adjoint(coords, self.mesh, self.quadRule) + mech_funcs = Mechanics.create_mechanics_functions(adjoint_func_space, mode2D='plane strain',materialModel=self.materialModel) + return mech_funcs.compute_strain_energy(U, ivs) + + nodal_forces = jax.jit(jax.grad(energy_function_all_dofs, argnums=0)) + + functionSpace = FunctionSpace.construct_function_space(self.mesh, self.quadRule) + ivsUpdateInverseFuncs = MechanicsInverse.create_ivs_update_inverse_functions(functionSpace, + "plane strain", + self.materialModel) + + residualInverseFuncs = MechanicsInverse.create_path_dependent_residual_inverse_functions(energy_function_coords) + + compute_df = jax.grad(self.total_work_increment, (0, 1, 2, 3, 6)) + + parameters = parameters.reshape(self.mesh.coords.shape) + gradient = np.zeros(parameters.shape) + mu = np.zeros(storedState[0][1].state_data.shape) + adjointLoad = np.zeros(storedState[0][0].shape) + + for step in reversed(range(1, self.steps+1)): + Uu = storedState[step][0] + Uu_prev = storedState[step-1][0] + p = storedState[step][1] + p_prev = storedState[step-1][1] + ivs = p.state_data + ivs_prev = p_prev.state_data + + dc_dcn = ivsUpdateInverseFuncs.ivs_update_jac_ivs_prev(self.dofManager.create_field(Uu, p.bc_data), ivs_prev) + + df_du, df_dun, df_dc, df_dcn, df_dx = compute_df(Uu, Uu_prev, ivs, ivs_prev, p, p_prev, parameters, nodal_forces) + + mu += df_dc + adjointLoad -= df_du + adjointLoad -= self.dofManager.get_unknown_values(ivsUpdateInverseFuncs.ivs_update_jac_disp_vjp(self.dofManager.create_field(Uu, p.bc_data), ivs_prev, mu)) + + n = self.dofManager.get_unknown_size() + p_objective = Objective.Params(bc_data=p.bc_data, state_data=p_prev.state_data) # remember R is a function of ivs_prev + self.objective.p = p_objective + self.objective.update_precond(Uu) # update preconditioner for use in cg (will converge in 1 iteration as long as the preconditioner is not approximate) + dRdu = linalg.LinearOperator((n, n), lambda V: onp.asarray(self.objective.hessian_vec(Uu, V))) + dRdu_decomp = linalg.LinearOperator((n, n), lambda V: onp.asarray(self.objective.apply_precond(V))) + adjointVector = linalg.cg(dRdu, onp.array(adjointLoad, copy=False), tol=1e-10, atol=0.0, M=dRdu_decomp)[0] + + gradient += df_dx + gradient += residualInverseFuncs.residual_jac_coords_vjp(Uu, p, ivs_prev, parameters, adjointVector) + gradient += ivsUpdateInverseFuncs.ivs_update_jac_coords_vjp(self.dofManager.create_field(Uu, p.bc_data), ivs_prev, parameters, mu) + + mu = np.einsum('ijk,ijkn->ijn', mu, dc_dcn) + mu += residualInverseFuncs.residual_jac_ivs_prev_vjp(Uu, p, ivs_prev, parameters, adjointVector) + mu += df_dcn + + adjointLoad = -df_dun + + return gradient.ravel() + + def compute_L2_norm_difference(self, uSteps, ivsSteps, bcsSteps, coordinates, nodal_forces): + index = (self.mesh.nodeSets['left'], 1) # arbitrarily choosing left side nodeset for reaction force + + numerator = 0.0 + denominator= 0.0 + for i in range(0, len(self.targetSteps)): + step = self.targetSteps[i] + Uu = uSteps[step] + bc_data = bcsSteps[step] + ivs = ivsSteps[step] + + U = self.dofManager.create_field(Uu, bc_data) + force = np.sum(np.array(nodal_forces(U, ivs, coordinates).at[index].get())) + + diff = force - self.targetForces[i] + numerator += diff*diff + denominator += self.targetForces[i]*self.targetForces[i] + + return np.sqrt(numerator/denominator) + + def target_curve_objective(self, storedState, parameters): + parameters = parameters.reshape(self.mesh.coords.shape) + + def energy_function_all_dofs(U, ivs, coords): + adjoint_func_space = AdjointFunctionSpace.construct_function_space_for_adjoint(coords, self.mesh, self.quadRule) + mech_funcs = Mechanics.create_mechanics_functions(adjoint_func_space, mode2D='plane strain',materialModel=self.materialModel) + return mech_funcs.compute_strain_energy(U, ivs) + + nodal_forces = jax.jit(jax.grad(energy_function_all_dofs, argnums=0)) + + uSteps = np.stack([storedState[i][0] for i in range(0, self.steps+1)], axis=0) + ivsSteps = np.stack([storedState[i][1].state_data for i in range(0, self.steps+1)], axis=0) + bcsSteps = np.stack([storedState[i][1].bc_data for i in range(0, self.steps+1)], axis=0) + + return self.compute_L2_norm_difference(uSteps, ivsSteps, bcsSteps, parameters, nodal_forces) + + def target_curve_gradient(self, storedState, parameters): + + def energy_function_coords(Uu, p, ivs_prev, coords): + adjoint_func_space = AdjointFunctionSpace.construct_function_space_for_adjoint(coords, self.mesh, self.quadRule) + mech_funcs = Mechanics.create_mechanics_functions(adjoint_func_space, mode2D='plane strain', materialModel=self.materialModel) + U = self.dofManager.create_field(Uu, p.bc_data) + return mech_funcs.compute_strain_energy(U, ivs_prev) + + def energy_function_all_dofs(U, ivs, coords): + adjoint_func_space = AdjointFunctionSpace.construct_function_space_for_adjoint(coords, self.mesh, self.quadRule) + mech_funcs = Mechanics.create_mechanics_functions(adjoint_func_space, mode2D='plane strain',materialModel=self.materialModel) + return mech_funcs.compute_strain_energy(U, ivs) + + nodal_forces = jax.jit(jax.grad(energy_function_all_dofs, argnums=0)) + + functionSpace = FunctionSpace.construct_function_space(self.mesh, self.quadRule) + ivsUpdateInverseFuncs = MechanicsInverse.create_ivs_update_inverse_functions(functionSpace, + "plane strain", + self.materialModel) + + residualInverseFuncs = MechanicsInverse.create_path_dependent_residual_inverse_functions(energy_function_coords) + + parameters = parameters.reshape(self.mesh.coords.shape) + + # derivatives of F + uSteps = np.stack([storedState[i][0] for i in range(0, self.steps+1)], axis=0) + ivsSteps = np.stack([storedState[i][1].state_data for i in range(0, self.steps+1)], axis=0) + bcsSteps = np.stack([storedState[i][1].bc_data for i in range(0, self.steps+1)], axis=0) + df_du, df_dc, gradient = jax.grad(self.compute_L2_norm_difference, (0, 1, 3))(uSteps, ivsSteps, bcsSteps, parameters, nodal_forces) + + mu = np.zeros(ivsSteps[0].shape) + adjointLoad = np.zeros(uSteps[0].shape) + + for step in reversed(range(1, self.steps+1)): + Uu = uSteps[step] + p = storedState[step][1] + p_prev = storedState[step-1][1] + ivs_prev = ivsSteps[step-1] + + dc_dcn = ivsUpdateInverseFuncs.ivs_update_jac_ivs_prev(self.dofManager.create_field(Uu, p.bc_data), ivs_prev) + + mu += df_dc[step] + adjointLoad -= df_du[step] + adjointLoad -= self.dofManager.get_unknown_values(ivsUpdateInverseFuncs.ivs_update_jac_disp_vjp(self.dofManager.create_field(Uu, p.bc_data), ivs_prev, mu)) + + n = self.dofManager.get_unknown_size() + p_objective = Objective.Params(bc_data=p.bc_data, state_data=p_prev.state_data) # remember R is a function of ivs_prev + self.objective.p = p_objective + self.objective.update_precond(Uu) # update preconditioner for use in cg (will converge in 1 iteration as long as the preconditioner is not approximate) + dRdu = linalg.LinearOperator((n, n), lambda V: onp.asarray(self.objective.hessian_vec(Uu, V))) + dRdu_decomp = linalg.LinearOperator((n, n), lambda V: onp.asarray(self.objective.apply_precond(V))) + adjointVector = linalg.cg(dRdu, onp.array(adjointLoad, copy=False), tol=1e-10, atol=0.0, M=dRdu_decomp)[0] + + gradient += residualInverseFuncs.residual_jac_coords_vjp(Uu, p, ivs_prev, parameters, adjointVector) + gradient += ivsUpdateInverseFuncs.ivs_update_jac_coords_vjp(self.dofManager.create_field(Uu, p.bc_data), ivs_prev, parameters, mu) + + mu = np.einsum('ijk,ijkn->ijn', mu, dc_dcn) + mu += residualInverseFuncs.residual_jac_ivs_prev_vjp(Uu, p, ivs_prev, parameters, adjointVector) + + adjointLoad = np.zeros(storedState[0][0].shape) + + return gradient.ravel() + + def test_total_work_gradient_with_adjoint_solve(self): + self.compute_objective_function = self.total_work_objective + self.compute_gradient = self.total_work_gradient + + stepSize = 1e-7 + + error = self.compute_finite_difference_error(stepSize, self.initialMesh.coords.ravel()) + self.assertLessEqual(error, 1e-7) + + def test_target_curve_gradient_with_adjoint_solve(self): + self.compute_objective_function = self.target_curve_objective + self.compute_gradient = self.target_curve_gradient + + self.targetSteps = [1, 2] + self.targetForces = [4.5, 5.5] # [4.542013626078756, 5.7673988583067555] actual forces + + stepSize = 1e-8 + + error = self.compute_finite_difference_error(stepSize, self.initialMesh.coords.ravel()) + self.assertLessEqual(error, 1e-6) + + + +if __name__ == '__main__': + unittest.main() diff --git a/optimism/inverse/test/test_J2Plastic_inverse.py b/optimism/inverse/test/test_J2Plastic_inverse.py new file mode 100644 index 0000000..ee733d7 --- /dev/null +++ b/optimism/inverse/test/test_J2Plastic_inverse.py @@ -0,0 +1,277 @@ +import unittest + +import jax +import jax.numpy as np + +from optimism import EquationSolver as EqSolver +from optimism import FunctionSpace +from optimism.material import J2Plastic as J2 +from optimism import Mechanics +from optimism import Mesh +from optimism import Objective +from optimism import QuadratureRule +from optimism import TensorMath + +from optimism.test.TestFixture import TestFixture +from optimism.test.MeshFixture import MeshFixture + +from optimism.inverse import MechanicsInverse + +def make_disp_grad_from_strain(strain): + return linalg.expm(strain) - np.identity(3) + +def small_strain_linear_hardening_analytic_gradients(disp_grad, state, state_previous, props): + sqrt32 = np.sqrt(3./2.) + mu = 0.5*props['elastic modulus']/(1.0 + props['poisson ratio']) + c0 = 3.*mu + props['hardening modulus'] + + vol_projection = np.zeros((9,9)).at[(0,0,0,4,4,4,8,8,8),(0,4,8,0,4,8,0,4,8)].set(1.) + dev_projection = np.eye(9,9) - vol_projection/3. + sym_projection = np.zeros((9,9)).at[(1,1,2,2,3,3,5,5,6,6,7,7),(1,3,2,6,1,3,5,7,2,6,5,7)].set(0.5).at[(0,4,8),(0,4,8)].set(1.) + + eqps = state[0] + deqps = eqps - state_previous[0] + + eps_pn = state_previous[1:] + eps = TensorMath.sym(disp_grad) + eps_e_tr = eps - eps_pn.reshape((3,3)) + + dev_eps_e_tr = TensorMath.dev(eps_e_tr) + dev_eps_e_tr_squared = np.tensordot(dev_eps_e_tr, dev_eps_e_tr) + n_tr = 1./np.sqrt(dev_eps_e_tr_squared) * dev_eps_e_tr + + deqps_deps = (2.*mu * sqrt32/ c0)*n_tr + deps_p_deps = (sqrt32 * deqps / np.sqrt(dev_eps_e_tr_squared)) * (np.tensordot(dev_projection,sym_projection,axes=1) - np.kron(n_tr.ravel(),n_tr.ravel()).reshape(9,9)) + deps_p_deps += sqrt32 * np.kron((2.*mu * sqrt32/ c0) * n_tr.ravel(), n_tr.ravel()).reshape(9,9) + + deqps_deqps_n = 3.*mu / c0 + deqps_deps_p_n = -(2.*mu * sqrt32/ c0)*n_tr + deqps_dc_n = np.append(deqps_deqps_n, deqps_deps_p_n.ravel()) + + deps_p_deqps_n = (3.*mu / c0 - 1.) * sqrt32 * n_tr + deps_p_deps_p_n = np.eye(9,9) - (sqrt32 * deqps / np.sqrt(dev_eps_e_tr_squared)) * (dev_projection - np.kron(n_tr.ravel(),n_tr.ravel()).reshape(9,9)) + deps_p_deps_p_n -= sqrt32 * np.kron((2.*mu * sqrt32/ c0) * n_tr.ravel(), n_tr.ravel()).reshape(9,9) + deps_p_dc_n = np.hstack((deps_p_deqps_n.ravel()[:,None],deps_p_deps_p_n)) + + return np.vstack((deqps_deps.ravel(), deps_p_deps)), np.vstack((deqps_dc_n, deps_p_dc_n)) + + +class J2MaterialPointUpdateGradsFixture(TestFixture): + def setUp(self): + E = 100.0 + poisson = 0.321 + Y0 = 0.2*E + H = 1.0e-2*E + + self.props = {'elastic modulus': E, + 'poisson ratio': poisson, + 'yield strength': Y0, + 'kinematics': 'small deformations', + 'hardening model': 'linear', + 'hardening modulus': H} + + materialModel = J2.create_material_model_functions(self.props) + + self.compute_state_new = jax.jit(materialModel.compute_state_new) + self.compute_initial_state = materialModel.compute_initial_state + + self.compute_state_new_derivs = jax.jit(jax.jacfwd(self.compute_state_new, (0, 1))) + + def test_jax_computation_of_state_derivs_at_elastic_step(self): + dt = 1.0 + + initial_state = self.compute_initial_state() + dc_dugrad, dc_dc_n = self.compute_state_new_derivs(np.zeros((3,3)), initial_state, dt) + + # Check derivative sizes + self.assertEqual(dc_dugrad.shape, (10,3,3)) + self.assertEqual(dc_dc_n.shape, (10,10)) + + # check derivatives w.r.t. displacement grad + self.assertArrayNear(dc_dugrad.ravel(), np.zeros((10,3,3)).ravel(), 12) + + # check derivatives w.r.t. previous step internal variables + self.assertArrayNear(dc_dc_n.ravel(), np.eye(10).ravel(), 12) + + def test_jax_computation_of_state_derivs_at_plastic_step(self): + dt = 1.0 + # get random displacement gradient + key = jax.random.PRNGKey(0) + dispGrad = jax.random.uniform(key, (3, 3)) + initial_state = self.compute_initial_state() + + state = self.compute_state_new(dispGrad, initial_state, dt) + dc_dugrad, dc_dc_n = self.compute_state_new_derivs(dispGrad, initial_state, dt) + + # Compute data for gold values (assuming small strain and linear kinematics) + dc_dugrad_gold, dc_dc_n_gold = small_strain_linear_hardening_analytic_gradients(dispGrad, state, initial_state, self.props) + + # Check derivative sizes + self.assertEqual(dc_dugrad.shape, (10,3,3)) + self.assertEqual(dc_dc_n.shape, (10,10)) + + # check derivatives w.r.t. displacement grad + self.assertArrayNear(dc_dugrad[0,:,:].ravel(), dc_dugrad_gold[0,:].ravel(), 12) + self.assertArrayNear(dc_dugrad[1:,:,:].ravel(), dc_dugrad_gold[1:,:].ravel() , 12) + + # check derivatives w.r.t. previous step internal variables + self.assertNear(dc_dc_n[0,0], dc_dc_n_gold[0,0], 12) + self.assertArrayNear(dc_dc_n[0,1:], dc_dc_n_gold[0,1:], 12) + self.assertArrayNear(dc_dc_n[1:,0], dc_dc_n_gold[1:,0], 12) + self.assertArrayNear(dc_dc_n[1:,1:].ravel(), dc_dc_n_gold[1:,1:].ravel() , 12) + + +class J2GlobalMeshUpdateGradsFixture(MeshFixture): + @classmethod + def setUpClass(cls): + dispGrad0 = np.array([[0.4, -0.2], + [-0.04, 0.68]]) + mf = MeshFixture() + cls.mesh, cls.U = mf.create_mesh_and_disp(4,4,[0.,1.],[0.,1.], + lambda x: dispGrad0@x) + + E = 100.0 + poisson = 0.321 + H = 1e-2 * E + Y0 = 0.3 * E + + cls.props = {'elastic modulus': E, + 'poisson ratio': poisson, + 'yield strength': Y0, + 'kinematics': 'small deformations', + 'hardening model': 'linear', + 'hardening modulus': H} + + cls.materialModel = J2.create_material_model_functions(cls.props) + cls.quadRule = QuadratureRule.create_quadrature_rule_on_triangle(degree=1) + cls.fs = FunctionSpace.construct_function_space(cls.mesh, cls.quadRule) + cls.mechFuncs = Mechanics.create_mechanics_functions(cls.fs, + "plane strain", + cls.materialModel) + cls.ivs_prev = cls.mechFuncs.compute_initial_state() + + EBCs = [FunctionSpace.EssentialBC(nodeSet='all_boundary', component=0), + FunctionSpace.EssentialBC(nodeSet='all_boundary', component=1)] + cls.dofManager = FunctionSpace.DofManager(cls.fs, 2, EBCs) + cls.Ubc = cls.dofManager.get_bc_values(cls.U) + + p = Objective.Params(None, cls.ivs_prev, None, None, None) + UuGuess = 0.0*cls.dofManager.get_unknown_values(cls.U) + + def compute_energy(Uu, p): + U = cls.dofManager.create_field(Uu, cls.Ubc) + internalVariables = p[1] + return cls.mechFuncs.compute_strain_energy(U, internalVariables) + + objective = Objective.Objective(compute_energy, UuGuess, p) + cls.Uu = EqSolver.nonlinear_equation_solve(objective, UuGuess, p, EqSolver.get_settings(), useWarmStart=False) + U = cls.dofManager.create_field(cls.Uu, cls.Ubc) + cls.ivs = cls.mechFuncs.compute_updated_internal_variables(U, cls.ivs_prev) + + def test_state_derivs_at_elastic_step(self): + + def update_internal_vars_test(Uu, ivs_prev): + U = self.dofManager.create_field(Uu) + ivs = self.mechFuncs.compute_updated_internal_variables(U, ivs_prev) + return ivs + + Uu = 0.0*self.dofManager.get_unknown_values(self.U) + + update_internal_variables_derivs = jax.jacfwd(update_internal_vars_test, (0,1)) + dc_du, dc_dc_n = update_internal_variables_derivs(Uu, self.ivs_prev) + + nElems = Mesh.num_elements(self.mesh) + nQpsPerElem = QuadratureRule.len(self.quadRule) + nIntVars = 10 + nFreeDofs = Uu.shape[0] + + self.assertEqual(dc_du.shape, (nElems,nQpsPerElem,nIntVars,nFreeDofs)) + self.assertEqual(dc_dc_n.shape, (nElems,nQpsPerElem,nIntVars,nElems,nQpsPerElem,nIntVars)) + + for i in range(0,nElems): + for j in range(0,nQpsPerElem): + self.assertArrayNear(dc_du[i,j,:,:].ravel(), np.zeros((nIntVars,nFreeDofs)).ravel(), 12) + self.assertArrayNear(dc_dc_n[i,j,:,i,j,:].ravel(), np.eye(nIntVars).ravel(), 12) + + def test_state_derivs_at_plastic_step(self): + + def update_internal_vars_test(U, ivs_prev): + ivs = self.mechFuncs.compute_updated_internal_variables(U, ivs_prev) + return ivs + + update_internal_variables_derivs = jax.jacfwd(update_internal_vars_test, (0,1)) + + U = self.dofManager.create_field(self.Uu, self.Ubc) + dc_du, dc_dc_n = update_internal_variables_derivs(U, self.ivs_prev) + + nElems = Mesh.num_elements(self.mesh) + nQpsPerElem = QuadratureRule.len(self.quadRule) + nIntVars = 10 + nDims = 2 + nNodes = Mesh.num_nodes(self.mesh) + + self.assertEqual(dc_du.shape, (nElems,nQpsPerElem,nIntVars,nNodes,nDims)) + self.assertEqual(dc_dc_n.shape, (nElems,nQpsPerElem,nIntVars,nElems,nQpsPerElem,nIntVars)) + + for i in range(0,nElems): + for j in range(0,nQpsPerElem): + state = self.ivs[i,j,:] + initial_state = self.ivs_prev[i,j,:] + + conn = self.mesh.conns[i] + Uele = U[conn] + shapeGrads = self.fs.shapeGrads[i,j,:,:] + dispGrad = TensorMath.tensor_2D_to_3D(np.tensordot(Uele,shapeGrads,axes=[0,0])) + Be_mat = np.zeros((9,6)).\ + at[(0,3),(0,1)].set(shapeGrads[0,0]).at[(1,4),(0,1)].set(shapeGrads[0,1]).\ + at[(0,3),(2,3)].set(shapeGrads[1,0]).at[(1,4),(2,3)].set(shapeGrads[1,1]).\ + at[(0,3),(4,5)].set(shapeGrads[2,0]).at[(1,4),(4,5)].set(shapeGrads[2,1]) + + dc_dugrad_gold, dc_dc_n_gold = small_strain_linear_hardening_analytic_gradients(dispGrad, state, initial_state, self.props) + + dc_duele_gold = np.tensordot(dc_dugrad_gold, Be_mat, axes=1) + dc_er_du = dc_du[i,j,:,:,:] + + self.assertArrayNear(dc_er_du[:,conn,:].ravel(), dc_duele_gold.reshape(10,3,2).ravel(), 10) + self.assertArrayNear(np.delete(dc_er_du, conn, axis=1).ravel(), np.zeros((10,nNodes-conn.shape[0],2)).ravel(), 10) + + for p in range(0,nElems): + for q in range(0,nQpsPerElem): + if(i == p and j == q): + self.assertArrayNear(dc_dc_n[i,j,:,i,j,:].ravel(), dc_dc_n_gold.ravel(), 10) + else: + self.assertArrayNear(dc_dc_n[i,j,:,p,q,:].ravel(), np.zeros((nIntVars,nIntVars)).ravel(), 10) + + def test_state_derivs_computed_locally_at_plastic_step(self): + + ivsUpdateInverseFuncs = MechanicsInverse.create_ivs_update_inverse_functions(self.fs, + "plane strain", + self.materialModel) + + U = self.dofManager.create_field(self.Uu, self.Ubc) + dc_dc_n = ivsUpdateInverseFuncs.ivs_update_jac_ivs_prev(U, self.ivs_prev) + + nElems = Mesh.num_elements(self.mesh) + nQpsPerElem = QuadratureRule.len(self.quadRule) + nIntVars = 10 + + self.assertEqual(dc_dc_n.shape, (nElems,nQpsPerElem,nIntVars,nIntVars)) + + for i in range(0,nElems): + for j in range(0,nQpsPerElem): + state = self.ivs[i,j,:] + initial_state = self.ivs_prev[i,j,:] + + conn = self.mesh.conns[i] + Uele = U[conn] + shapeGrads = self.fs.shapeGrads[i,j,:,:] + dispGrad = TensorMath.tensor_2D_to_3D(np.tensordot(Uele,shapeGrads,axes=[0,0])) + + _, dc_dc_n_gold = small_strain_linear_hardening_analytic_gradients(dispGrad, state, initial_state, self.props) + + self.assertArrayNear(dc_dc_n[i,j,:,:].ravel(), dc_dc_n_gold.ravel(), 10) + + + +if __name__ == '__main__': + unittest.main()