From b814c0402a82162689ef055e6279778c4b96ebec Mon Sep 17 00:00:00 2001 From: ralberd Date: Wed, 15 Jan 2025 14:21:00 -0700 Subject: [PATCH 01/19] adds material_qoi output for multi blocks --- optimism/Mechanics.py | 31 ++++++++++++++++++++++++++++--- 1 file changed, 28 insertions(+), 3 deletions(-) diff --git a/optimism/Mechanics.py b/optimism/Mechanics.py index 751ee25c..619969d7 100644 --- a/optimism/Mechanics.py +++ b/optimism/Mechanics.py @@ -255,12 +255,37 @@ def compute_output_energy_densities_and_stresses(U, stateVariables, dt=0.0): stresses = stresses.at[elemIds].set(blockStresses) return energy_densities, stresses - def compute_initial_state(): return _compute_initial_state_multi_block(fs, materialModels) - - return MechanicsFunctions(compute_strain_energy, jit(compute_updated_internal_variables), jit(compute_element_stiffnesses), jit(compute_output_energy_densities_and_stresses), compute_initial_state, None, None) + def qoi_to_lagrangian_qoi(compute_material_qoi): + def L(U, gradU, Q, X, dt): + return compute_material_qoi(gradU, Q, dt) + return L + + def integrated_material_qoi(U, stateVariables, dt=0.0): + material_qoi = np.zeros((Mesh.num_elements(fs.mesh), len(fs.quadratureRule))) + for blockKey in materialModels: + elemIds = fs.mesh.blocks[blockKey] + block_qoi = FunctionSpace.integrate_over_block(fs, U, stateVariables, dt, + qoi_to_lagrangian_qoi(materialModels[blockKey].compute_material_qoi), + elemIds, + modify_element_gradient=modify_element_gradient) + material_qoi = material_qoi.at[elemIds].set(block_qoi) + return material_qoi + + def compute_output_material_qoi(U, stateVariables, dt=0.0): + material_qoi = np.zeros((Mesh.num_elements(fs.mesh), len(fs.quadratureRule))) + for blockKey in materialModels: + elemIds = fs.mesh.blocks[blockKey] + block_qoi = FunctionSpace.evaluate_on_block(fs, U, stateVariables, dt, + qoi_to_lagrangian_qoi(materialModels[blockKey].compute_material_qoi), + elemIds, + modify_element_gradient=modify_element_gradient) + material_qoi = material_qoi.at[elemIds].set(block_qoi) + return material_qoi + + return MechanicsFunctions(compute_strain_energy, jit(compute_updated_internal_variables), jit(compute_element_stiffnesses), jit(compute_output_energy_densities_and_stresses), compute_initial_state, integrated_material_qoi, jit(compute_output_material_qoi)) ###### From 09e1ab7b4210c05bd75811a34f40fb5e374b2337 Mon Sep 17 00:00:00 2001 From: ralberd Date: Wed, 15 Jan 2025 14:22:59 -0700 Subject: [PATCH 02/19] updates ReadExodusMesh to work if there is no element map --- optimism/ReadExodusMesh.py | 11 +++++++++-- optimism/test/test_ReadExodusMesh.py | 7 +++++++ 2 files changed, 16 insertions(+), 2 deletions(-) diff --git a/optimism/ReadExodusMesh.py b/optimism/ReadExodusMesh.py index facb4e40..b20d4745 100644 --- a/optimism/ReadExodusMesh.py +++ b/optimism/ReadExodusMesh.py @@ -92,11 +92,18 @@ def _read_blocks(exodusDataset): def _read_block_maps(exodusDataset, blocks): block_maps = {} - elementMap = exodusDataset.variables['elem_num_map'] + if 'elem_num_map' in exodusDataset.variables: + elementMap = exodusDataset.variables['elem_num_map'] + else: + nEle = 1 + for blockName, blockElems in blocks.items(): + nEle = nEle + len(blockElems) + elementMap = onp.arange(1, nEle) + firstElemInBlock = 0 for blockName, blockElems in blocks.items(): nElemsInBlock = len(blockElems) - block_maps[blockName] = elementMap[firstElemInBlock:firstElemInBlock+nElemsInBlock] + block_maps[blockName] = onp.array(elementMap[firstElemInBlock:firstElemInBlock + nElemsInBlock]) firstElemInBlock += nElemsInBlock return block_maps diff --git a/optimism/test/test_ReadExodusMesh.py b/optimism/test/test_ReadExodusMesh.py index 6262ef24..40ab927d 100644 --- a/optimism/test/test_ReadExodusMesh.py +++ b/optimism/test/test_ReadExodusMesh.py @@ -90,6 +90,13 @@ def test_side_set_sizes(self): self.assertEqual(self.mesh.sideSets["right"].shape[0], 2) self.assertEqual(self.mesh.sideSets["top"].shape[0], 4) + @TestFixture.unittest.skipIf(not haveNetCDF, skipMessage) + def test_block_maps_no_map(self): + mesh = ReadExodusMesh.read_exodus_mesh(pathlib.Path(__file__).parent.joinpath('square_mesh_no_map.exo')) + block1Map = mesh.block_maps['Block1'] + self.assertEqual(block1Map.size, 10) + self.assertArrayEqual(block1Map, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]) + @TestFixture.unittest.skipIf(not haveNetCDF, skipMessage) def test_block_maps(self): # ncdump output (see patch_2_blocks.txt) From 187e0377beebe28cf6887a0e58e5b149448b323b Mon Sep 17 00:00:00 2001 From: ralberd Date: Wed, 15 Jan 2025 14:23:39 -0700 Subject: [PATCH 03/19] adds material_qoi to NeoHookean material; stores det(F) --- optimism/material/Neohookean.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/optimism/material/Neohookean.py b/optimism/material/Neohookean.py index cfb5ec76..e0a52e99 100644 --- a/optimism/material/Neohookean.py +++ b/optimism/material/Neohookean.py @@ -32,9 +32,15 @@ def compute_state_new(dispGrad, internalVars, dt): density = properties.get('density') + def compute_material_qoi(dispGrad, internalVars, dt): + del internalVars + del dt + return _compute_volumetric_jacobian(dispGrad) + return MaterialModel(compute_energy_density = strain_energy, compute_initial_state = make_initial_state, compute_state_new = compute_state_new, + compute_material_qoi = compute_material_qoi, density = density) @@ -78,3 +84,7 @@ def _compute_state_new(dispGrad, internalVars, props): del dispGrad del props return internalVars + +def _compute_volumetric_jacobian(dispGrad): + F = dispGrad + np.eye(3) + return np.linalg.det(F) \ No newline at end of file From 94ac0542142d325bc33c48cd8a8c12d2e837b407 Mon Sep 17 00:00:00 2001 From: ralberd Date: Wed, 15 Jan 2025 14:24:05 -0700 Subject: [PATCH 04/19] adds ThirdMediumNeoHookean material for testing third medium formulations --- optimism/material/ThirdMediumNeoHookean.py | 65 ++++++++++++++++++++++ 1 file changed, 65 insertions(+) create mode 100644 optimism/material/ThirdMediumNeoHookean.py diff --git a/optimism/material/ThirdMediumNeoHookean.py b/optimism/material/ThirdMediumNeoHookean.py new file mode 100644 index 00000000..257b3404 --- /dev/null +++ b/optimism/material/ThirdMediumNeoHookean.py @@ -0,0 +1,65 @@ +import jax.numpy as np + +from optimism.material.MaterialModel import MaterialModel + +# props +PROPS_MU = 0 +PROPS_KAPPA = 1 +PROPS_LAMBDA = 2 + +def create_material_model_functions(properties): + + density = properties.get('density') + props = _make_properties(properties['bulk modulus'], + properties['shear modulus']) + + def strain_energy(dispGrad, internalVars, dt): + del internalVars + del dt + # return neo_hookean_energy_density(dispGrad, props) + return stabilized_neo_hookean_energy_density(dispGrad, props) + + def compute_state_new(dispGrad, internalVars, dt): + del dispGrad + del dt + return internalVars + + def compute_material_qoi(dispGrad, internalVars, dt): + del internalVars + del dt + return _compute_volumetric_jacobian(dispGrad) + + return MaterialModel(compute_energy_density = strain_energy, + compute_initial_state = make_initial_state, + compute_state_new = compute_state_new, + compute_material_qoi = compute_material_qoi, + density = density) + +def make_initial_state(): + return np.array([]) + +def _make_properties(kappa, mu): + lamda = kappa - (2.0/3.0) * mu + return np.array([mu, kappa, lamda]) + +def neo_hookean_energy_density(dispGrad, props): + F = dispGrad + np.eye(3) + J = np.linalg.det(F) + J23 = np.power(J, -2.0/3.0) + I1Bar = J23*np.tensordot(F,F) + Wvol = 0.5*props[PROPS_KAPPA]*(np.log(J)**2) + Wdev = 0.5*props[PROPS_MU]*(I1Bar - 3.0) + return Wdev + Wvol + +def stabilized_neo_hookean_energy_density(dispGrad, props): + F = dispGrad + np.eye(3) + J = np.linalg.det(F) + I1 = np.tensordot(F,F) + Wiso = props[PROPS_MU]/2.0 * (I1 - 3.0) + alpha = 1.0 + props[PROPS_MU]/props[PROPS_LAMBDA] + Wvol = props[PROPS_LAMBDA]/2.0 * (J - alpha)**2 + return Wiso + Wvol + +def _compute_volumetric_jacobian(dispGrad): + F = dispGrad + np.eye(3) + return np.linalg.det(F) \ No newline at end of file From 197fcad73e043830f5db5f091dc919de55c2d307 Mon Sep 17 00:00:00 2001 From: ralberd Date: Wed, 15 Jan 2025 14:25:03 -0700 Subject: [PATCH 05/19] adds test to try out different energy formulations for third medium --- optimism/material/test/test_ThirdMedium.py | 128 +++++++++++++++++++++ 1 file changed, 128 insertions(+) create mode 100644 optimism/material/test/test_ThirdMedium.py diff --git a/optimism/material/test/test_ThirdMedium.py b/optimism/material/test/test_ThirdMedium.py new file mode 100644 index 00000000..ec348f70 --- /dev/null +++ b/optimism/material/test/test_ThirdMedium.py @@ -0,0 +1,128 @@ +import unittest +from optimism.test.TestFixture import TestFixture +import jax +import jax.numpy as np +from optimism.material import Neohookean + +from matplotlib import pyplot as plt + +plotting = True + +def energy_neo_hookean_adagio(dispGrad, kappa, mu): + F = dispGrad + np.eye(3) + J = np.linalg.det(F) + J23 = np.power(J, -2.0/3.0) + I1Bar = J23*np.tensordot(F,F) + Wvol = 0.5*kappa*(0.5*J**2 - 0.5 - np.log(J)) + Wdev = 0.5*mu*(I1Bar - 3.0) + return Wdev + Wvol + +def energy_neo_hookean_rokos_volumetric(dispGrad, kappa, mu): + F = dispGrad + np.eye(3) + J = np.linalg.det(F) + return kappa/2.0 * np.log(J)**2 + +def energy_neo_hookean_rokos_isochoric(dispGrad, kappa, mu): + F = dispGrad + np.eye(3) + J = np.linalg.det(F) + J23 = np.power(J, -2.0/3.0) + I1Bar = J23*np.tensordot(F,F) + return mu / 2.0 * (I1Bar - 3.0) + +def energy_stable_neo_hookean(dispGrad, kappa, mu): + F = dispGrad + np.eye(3) + J = np.linalg.det(F) + I1 = np.tensordot(F,F) + Wiso = mu/2.0*(I1 - 3.0) + lam = kappa - (2.0/3.0) * mu + alpha = 1.0 + mu/lam + Wvol = lam/2.0 * (J - alpha)**2 + return Wiso + Wvol + +def energy_stable_arruda_boyce(dispGrad, kappa, mu): + F = dispGrad + np.eye(3) + J = np.linalg.det(F) + I1 = np.tensordot(F,F) + beta1 = 1.0 + beta2 = 1.0 + Wiso = mu/2.0*(I1 - 3.0) + beta1*mu/4.0*(I1**2 - 9.0) + beta2*mu/6.0*(I1**3 - 27.0) + lam = kappa - (2.0/3.0) * mu + alpha = 1.0 + mu/lam * (1.0 + 3.0*beta1 + 9.0*beta2) + Wvol = lam/2.0 * (J - alpha)**2 + return Wiso + Wvol + +class ThirdMediumModelFixture(TestFixture): + def setUp(self): + self.kappa = 100.0 + self.mu = 10.0 + + def test_energy_consistency(self): + # generate a random displacement gradient + key = jax.random.PRNGKey(1) + dispGrad = jax.random.uniform(key, (3, 3)) + + # compare to Neo-Hookean model + props = { + 'elastic modulus': 9.0*self.kappa*self.mu / (3.0*self.kappa + self.mu), + 'poisson ratio': (3.0*self.kappa - 2.0*self.mu) / 2.0 / (3.0*self.kappa + self.mu), + 'version': 'adagio' + } + ref_material = Neohookean.create_material_model_functions(props) + energy_gold = ref_material.compute_energy_density(dispGrad, ref_material.compute_initial_state(), dt=0.0) + + energy = energy_neo_hookean_adagio(dispGrad, self.kappa, self.mu) + self.assertAlmostEqual(energy, energy_gold, 12) + + def test_plot_biaxial_response(self): + n_steps = 100 + f_vals = np.linspace(0.0, 1.0, n_steps) + Fs = jax.vmap( + lambda fii: np.array( + [[fii, 0.0, 0.0], + [0.0, fii, 0.0], + [0.0, 0.0, 1.0]] + ) + )(f_vals) + + energies = np.zeros((4,n_steps)) + Jvals = np.zeros(n_steps) + for n, F in enumerate(Fs): + dispGrad = F - np.eye(3) + energies = energies.at[0,n].set(energy_neo_hookean_rokos_volumetric(dispGrad, self.kappa, self.mu)) + energies = energies.at[1,n].set(energy_neo_hookean_rokos_isochoric(dispGrad, self.kappa, self.mu)) + energies = energies.at[2,n].set(energy_stable_neo_hookean(dispGrad, self.kappa, self.mu)) + energies = energies.at[3,n].set(energy_stable_arruda_boyce(dispGrad, self.kappa, self.mu)) + Jvals = Jvals.at[n].set(np.linalg.det(F)) + + if plotting: + plt.figure(1) + fig, ax1 = plt.subplots() + + # energy + ax1.set_xlabel('F_11 = F_22 (F_33 = 1)') + ax1.set_ylabel('Energy') + legends1 = ax1.plot(f_vals, energies[0,:], linestyle='--', color='r') + legends2 = ax1.plot(f_vals, energies[1,:], linestyle='--', color='b') + legends3 = ax1.plot(f_vals, energies[2,:], linestyle='--', color='g') + legends4 = ax1.plot(f_vals, energies[3,:], linestyle='--', color='c') + + ax2 = ax1.twinx() # instantiate a second axes that shares the same x-axis + ax2.set_ylabel('J') # we already handled the x-label with ax1 + legends5 = ax2.plot(f_vals, Jvals, linestyle='--', color='k') + + ax1.legend(legends1+legends2+legends3+legends4,[ + "Rokos Neo Hookean Volumetric term", + "Rokos Neo Hookean Isochoric term", + "Stable Neohookean energy", + "Stable Arruda Boyce energy", + "Volume change J" + ], + loc=0, frameon=True) + + fig.tight_layout() # otherwise the right y-label is slightly clipped + + plt.savefig('energy_vs_compression.png') + + +if __name__ == "__main__": + unittest.main() From bcb3301f3e77da51600f0183465b3c9cd2411693 Mon Sep 17 00:00:00 2001 From: ralberd Date: Sun, 9 Feb 2025 08:01:50 -0700 Subject: [PATCH 06/19] adds quadratic lagrange triangle elements --- optimism/Interpolants.py | 141 ++++++++++++++++++-- optimism/ReadExodusMesh.py | 18 ++- optimism/test/test_Interpolants.py | 201 ++++++++++++++++++++++++----- 3 files changed, 308 insertions(+), 52 deletions(-) diff --git a/optimism/Interpolants.py b/optimism/Interpolants.py index 76f1dee9..934c63fb 100644 --- a/optimism/Interpolants.py +++ b/optimism/Interpolants.py @@ -44,9 +44,9 @@ class ShapeFunctions(eqx.Module): is the number of nodes in the element (which is equal to the number of shape functions). gradients: Values of the parametric gradients of the shape functions. - Shape is ``(nPts, nDim, nNodes)``, where ``nDim`` is the number + Shape is ``(nPts, nNodes, nDim)``, where ``nDim`` is the number of spatial dimensions. Line elements are an exception, which - have shape ``(nPts, nNdodes)``. + have shape ``(nNodes, nPts)``. """ values: Float[Array, "nq nn"] gradients: Float[Array, "nq nn nd"] @@ -59,6 +59,8 @@ def __iter__(self): LINE_ELEMENT = 0 TRIANGLE_ELEMENT = 1 TRIANGLE_ELEMENT_WITH_BUBBLE = 2 +LAGRANGE_LINE_ELEMENT = 3 +LAGRANGE_TRIANGLE_ELEMENT = 4 def make_parent_elements(degree): @@ -86,6 +88,33 @@ def get_lobatto_nodes_1d(degree): return xn +def make_lagrange_parent_element_1d(degree): + """Lagrange Interpolation points on the unit interval [0, 1]. + Only implemented for second degree + """ + if degree != 2: + raise NotImplementedError + + xn = np.array([0.0, 0.5, 1.0]) + vertexPoints = np.array([0, 2], dtype=np.int32) + interiorPoints = np.array([1], dtype=np.int32) + return ParentElement(LAGRANGE_LINE_ELEMENT, int(degree), xn, vertexPoints, None, interiorPoints) + + +def vander1d(x, degree): + x = onp.asarray(x) + A = onp.zeros((x.shape[0], degree + 1)) + dA = onp.zeros((x.shape[0], degree + 1)) + domain = [0.0, 1.0] + for i in range(degree + 1): + p = onp.polynomial.Legendre.basis(i, domain=domain) + p *= onp.sqrt(2.0*i + 1.0) # keep polynomial orthonormal + A[:, i] = p(x) + dp = p.deriv() + dA[:, i] = dp(x) + return A, dA + + def shape1d(degree, nodalPoints, evaluationPoints): """Evaluate shape functions and derivatives at points in the master element. @@ -108,18 +137,34 @@ def shape1d(degree, nodalPoints, evaluationPoints): return ShapeFunctions(shape, dshape) -def vander1d(x, degree): - x = onp.asarray(x) - A = onp.zeros((x.shape[0], degree + 1)) - dA = onp.zeros((x.shape[0], degree + 1)) - domain = [0.0, 1.0] - for i in range(degree + 1): - p = onp.polynomial.Legendre.basis(i, domain=domain) - p *= onp.sqrt(2.0*i + 1.0) # keep polynomial orthonormal - A[:, i] = p(x) - dp = p.deriv() - dA[:, i] = dp(x) - return A, dA +def shape1d_lagrange(degree, nodalPoints, evaluationPoints): + """Evaluate Lagrange shape functions and derivatives at points in the parent element. + Only implemented for second degree + + Returns: + Shape function values and shape function derivatives at ``evaluationPoints``, + in a tuple (``shape``, ``dshape``). + shapes: [nNodes, nEvalPoints] + dshapes: [nNodes, nEvalPoints] + """ + if degree != 2: + raise NotImplementedError + + denom1 = (nodalPoints[0] - nodalPoints[1]) * (nodalPoints[0] - nodalPoints[2]) + denom2 = (nodalPoints[1] - nodalPoints[0]) * (nodalPoints[1] - nodalPoints[2]) + denom3 = (nodalPoints[2] - nodalPoints[0]) * (nodalPoints[2] - nodalPoints[1]) + + shape1 = (evaluationPoints - nodalPoints[1])*(evaluationPoints - nodalPoints[2]) / denom1 + shape2 = (evaluationPoints - nodalPoints[0])*(evaluationPoints - nodalPoints[2]) / denom2 + shape3 = (evaluationPoints - nodalPoints[0])*(evaluationPoints - nodalPoints[1]) / denom3 + shape = np.stack((shape1, shape2, shape3)) + + dshape1 = (2.0*evaluationPoints - nodalPoints[2] - nodalPoints[1]) / denom1 + dshape2 = (2.0*evaluationPoints - nodalPoints[2] - nodalPoints[0]) / denom2 + dshape3 = (2.0*evaluationPoints - nodalPoints[1] - nodalPoints[0]) / denom3 + dshape = np.stack((dshape1, dshape2, dshape3)) + + return ShapeFunctions(shape, dshape) def make_parent_element_2d(degree): @@ -169,6 +214,27 @@ def make_parent_element_2d(degree): return ParentElement(TRIANGLE_ELEMENT, int(degree), points, vertexPoints, facePoints, interiorPoints) +def make_lagrange_parent_element_2d(degree): + """Lagrange interpolation points on the triangle + Only implemented for second degree triangles. + + Convention for numbering: + + 2 + o + | \ + 4 o o 1 + | \ + o--o--o + 5 3 0 + """ + if degree != 2: + raise NotImplementedError + + xn = np.array([[1.0, 0.0], [0.5, 0.5], [0.0, 1.0], [0.5, 0.0], [0.0, 0.5], [0.0, 0.0]]) + vertexPoints = np.array([0, 2, 5], dtype=np.int32) + facePoints = np.array([[0, 1, 2], [2, 4, 5], [5, 3, 0]], dtype=np.int32) + return ParentElement(LAGRANGE_TRIANGLE_ELEMENT, int(degree), xn, vertexPoints, facePoints, None) def pascal_triangle_monomials(degree): p = [] @@ -261,6 +327,49 @@ def shape2d(degree, nodalPoints, evaluationPoints): return ShapeFunctions(np.asarray(shapes), np.asarray(dshapes)) +def shape2d_lagrange(degree, nodalPoints, evaluationPoints): + """Evaluate Lagrange shape functions and derivatives at points in the parent element. + Only implemented for second degree + + Reference: + T. Hughes. "The Finite Element Method" + Appendix 3.I + """ + if degree != 2: + raise NotImplementedError + + numEvalPoints = evaluationPoints.shape[0] + r = evaluationPoints[:,0] + s = evaluationPoints[:,1] + # t = 1.0 - r - s + + shape0 = 2.0*r*r - r # r * (2.0 * r - 1.0) + shape1 = 4.0 * r * s # 4.0 * r * s + shape2 = 2.0*s*s - s # s * (2.0 * s - 1.0) + shape3 = 4.0*(r - r*s - r*r) # 4.0 * r * t + shape4 = 4.0*(s - r*s - s*s) # 4.0 * s * t + shape5 = 1.0 - 3.0*(r + s) + 4.0*r*s + 2.0*(r*r + s*s) # t * (2.0 * t - 1.0) + shape = np.stack((shape0, shape1, shape2, shape3, shape4, shape5)).T + + dshape0_dr = 4.0*r - 1.0 + dshape0_ds = np.zeros(numEvalPoints) + dshape1_dr = 4.0*s + dshape1_ds = 4.0*r + dshape2_dr = np.zeros(numEvalPoints) + dshape2_ds = 4.0*s - 1.0 + dshape3_dr = 4.0*(1.0 - s - 2.0*r) + dshape3_ds = -4.0*r + dshape4_dr = -4.0*s + dshape4_ds = 4.0*(1.0 - r - 2.0*s) + dshape5_dr = 4.0*(r + s) - 3.0 + dshape5_ds = 4.0*(r + s) - 3.0 + dshape_dr = np.stack((dshape0_dr, dshape1_dr, dshape2_dr, dshape3_dr, dshape4_dr, dshape5_dr)).T + dshape_ds = np.stack((dshape0_ds, dshape1_ds, dshape2_ds, dshape3_ds, dshape4_ds, dshape5_ds)).T + dshape = np.stack((dshape_dr, dshape_ds), axis=2) + + return ShapeFunctions(shape, dshape) + + def compute_shapes(parentElement, evaluationPoints): if parentElement.elementType == LINE_ELEMENT: return shape1d(parentElement.degree, parentElement.coordinates, evaluationPoints) @@ -268,6 +377,10 @@ def compute_shapes(parentElement, evaluationPoints): return shape2d(parentElement.degree, parentElement.coordinates, evaluationPoints) elif parentElement.elementType == TRIANGLE_ELEMENT_WITH_BUBBLE: return shape2dBubble(parentElement, evaluationPoints) + elif parentElement.elementType == LAGRANGE_LINE_ELEMENT: + return shape1d_lagrange(parentElement.degree, parentElement.coordinates, evaluationPoints) + elif parentElement.elementType == LAGRANGE_TRIANGLE_ELEMENT: + return shape2d_lagrange(parentElement.degree, parentElement.coordinates, evaluationPoints) else: raise ValueError('Unknown element type.') diff --git a/optimism/ReadExodusMesh.py b/optimism/ReadExodusMesh.py index b20d4745..bbbd77d1 100644 --- a/optimism/ReadExodusMesh.py +++ b/optimism/ReadExodusMesh.py @@ -3,11 +3,15 @@ from optimism import Interpolants import netCDF4 import numpy as onp +from enum import Enum, auto -exodusToNativeTri6NodeOrder = np.array([0, 3, 1, 5, 4, 2]) +class InterpolationType(Enum): + LOBATTO = auto() + LAGRANGE = auto() +exodusToNativeTri6NodeOrder = np.array([0, 3, 1, 5, 4, 2]) -def read_exodus_mesh(fileName): +def read_exodus_mesh(fileName, interpolationType = InterpolationType.LOBATTO): with netCDF4.Dataset(fileName) as exData: coords = _read_coordinates(exData) conns, blocks = _read_blocks(exData) @@ -17,15 +21,21 @@ def read_exodus_mesh(fileName): elementType = _read_element_type(exData).lower() if elementType == "tri3" or elementType == "tri": - basis, basis1d = Interpolants.make_parent_elements(degree = 1) + degree = 1 simplexNodesOrdinals = np.arange(coords.shape[0]) elif elementType == "tri6": - basis, basis1d = Interpolants.make_parent_elements(degree = 2) + degree = 2 simplexNodesOrdinals = _get_vertex_nodes_from_exodus_tri6_mesh(conns) conns = conns[:, exodusToNativeTri6NodeOrder] else: raise + if interpolationType == InterpolationType.LAGRANGE: + basis1d = Interpolants.make_lagrange_parent_element_1d(degree) + basis = Interpolants.make_lagrange_parent_element_2d(degree) + else: + basis, basis1d = Interpolants.make_parent_elements(degree) + return Mesh.Mesh(coords, conns, simplexNodesOrdinals, basis, basis1d, blocks, nodeSets, sideSets, blockMaps) diff --git a/optimism/test/test_Interpolants.py b/optimism/test/test_Interpolants.py index 6354b213..ff71bda9 100644 --- a/optimism/test/test_Interpolants.py +++ b/optimism/test/test_Interpolants.py @@ -7,7 +7,7 @@ from optimism import Interpolants from optimism import QuadratureRule -from . import TestFixture +import TestFixture tol = 1e-14 @@ -22,8 +22,47 @@ def generate_random_points_in_triangle(npts): points = jax.numpy.column_stack((x,y)) return onp.asarray(points) +class ElementFixture(TestFixture.TestFixture): + def assertCorrectInterpolationOn2DElement(self, element, numRandomEvalPoints, shape2dFun): + degree = element.degree + x = generate_random_points_in_triangle(numRandomEvalPoints) + polyCoeffs = onp.fliplr(onp.triu(onp.ones((degree+1,degree+1)))) -class TestInterpolants(TestFixture.TestFixture): + shape, _ = shape2dFun(degree, element.coordinates, x) + fn = onp.polynomial.polynomial.polyval2d(element.coordinates[:,0], + element.coordinates[:,1], + polyCoeffs) + fInterpolated = onp.dot(shape, fn) + + expected = onp.polynomial.polynomial.polyval2d(x[:,0], x[:,1], polyCoeffs) + self.assertArrayNear(expected, fInterpolated, 14) + + def assertCorrectGradInterpolationOn2DElement(self, element, numRandomEvalPoints, shape2dFun): + degree = element.degree + x = generate_random_points_in_triangle(numRandomEvalPoints) + poly = onp.fliplr(onp.triu(onp.ones((degree+1,degree+1)))) + + _, dShape = shape2dFun(degree, element.coordinates, x) + fn = onp.polynomial.polynomial.polyval2d(element.coordinates[:,0], + element.coordinates[:,1], + poly) + dfInterpolated = onp.einsum('qai,a->qi',dShape, fn) + + # exact x derivative + direction = 0 + DPoly = onp.polynomial.polynomial.polyder(poly, 1, scl=1, axis=direction) + expected0 = onp.polynomial.polynomial.polyval2d(x[:,0], x[:,1], DPoly) + + self.assertArrayNear(expected0, dfInterpolated[:,0], 14) + + direction = 1 + DPoly = onp.polynomial.polynomial.polyder(poly, 1, scl=1, axis=direction) + expected1 = onp.polynomial.polynomial.polyval2d(x[:,0], x[:,1], DPoly) + + self.assertArrayNear(expected1, dfInterpolated[:,1], 14) + + +class TestInterpolants(ElementFixture): def setUp(self): maxDegree = 5 @@ -117,44 +156,16 @@ def test_shape_kronecker_delta_property(self): def test_interpolation(self): - x = generate_random_points_in_triangle(1) + numRandomEvalPoints = 1 for element in self.elements: - degree = element.degree - polyCoeffs = onp.fliplr(onp.triu(onp.ones((degree+1,degree+1)))) - expected = onp.polynomial.polynomial.polyval2d(x[:,0], x[:,1], polyCoeffs) - - shape, _ = Interpolants.shape2d(degree, element.coordinates, x) - fn = onp.polynomial.polynomial.polyval2d(element.coordinates[:,0], - element.coordinates[:,1], - polyCoeffs) - fInterpolated = onp.dot(shape, fn) - self.assertArrayNear(expected, fInterpolated, 14) + self.assertCorrectInterpolationOn2DElement(element, numRandomEvalPoints, Interpolants.shape2d) def test_grad_interpolation(self): - x = generate_random_points_in_triangle(1) + numRandomEvalPoints = 1 for element in self.elements: - degree = element.degree - poly = onp.fliplr(onp.triu(onp.ones((degree+1,degree+1)))) - - _, dShape = Interpolants.shape2d(degree, element.coordinates, x) - fn = onp.polynomial.polynomial.polyval2d(element.coordinates[:,0], - element.coordinates[:,1], - poly) - dfInterpolated = onp.einsum('qai,a->qi',dShape, fn) - - # exact x derivative - direction = 0 - DPoly = onp.polynomial.polynomial.polyder(poly, 1, scl=1, axis=direction) - expected0 = onp.polynomial.polynomial.polyval2d(x[:,0], x[:,1], DPoly) - - self.assertArrayNear(expected0, dfInterpolated[:,0], 13) - - direction = 1 - DPoly = onp.polynomial.polynomial.polyder(poly, 1, scl=1, axis=direction) - expected1 = onp.polynomial.polynomial.polyval2d(x[:,0], x[:,1], DPoly) + self.assertCorrectGradInterpolationOn2DElement(element, numRandomEvalPoints, Interpolants.shape2d) - self.assertArrayNear(expected1, dfInterpolated[:,1], 13) def no_test_plot_high_order_nodes(self): degree = 10 @@ -228,5 +239,127 @@ def no_test_plot_shape_functions(self): plt.show() +class TestLagrangeInterpolants(ElementFixture): + + def setUp(self): + self.degree = 2 # only second degree is implemented + self.qr1d = QuadratureRule.create_quadrature_rule_1D(degree=2) + self.nQPts1d = len(self.qr1d) + self.qr = QuadratureRule.create_quadrature_rule_on_triangle(degree=2) + self.nQPts = len(self.qr) + + def test_degree_other_than_two_throws(self): + badDegrees = [1, 3, 4, 5] + dummyCoords = np.array([[0.0, 1.0], [1.0, 0.0]]) + for degree in badDegrees: + self.assertRaises(NotImplementedError, Interpolants.make_lagrange_parent_element_1d, degree) + self.assertRaises(NotImplementedError, Interpolants.make_lagrange_parent_element_2d, degree) + self.assertRaises(NotImplementedError, Interpolants.shape1d_lagrange, degree, dummyCoords, self.qr1d) + self.assertRaises(NotImplementedError, Interpolants.shape2d_lagrange, degree, dummyCoords, self.qr) + + def test_1D_interpolant_points(self): + element = Interpolants.make_lagrange_parent_element_1d(self.degree) + p = element.coordinates + self.assertTrue(np.all(p >= 0.0) and onp.all(p <= 1.0)) + self.assertNear(p[element.interiorNodes], 0.5, 16) + self.assertArrayNear(p[element.vertexNodes], onp.array([0.0, 1.0]), 16) + self.assertTrue(element.faceNodes is None) + + def test_tri_interpolant_points(self): + element = Interpolants.make_lagrange_parent_element_2d(self.degree) + p = element.coordinates + # x conditions + self.assertTrue(np.all(p[:,0] >= -tol)) + self.assertTrue(np.all(p[:,0] <= 1.0 + tol)) + # y conditions + self.assertTrue(np.all(p[:,1] >= -tol)) + self.assertTrue(np.all(p[:,1] <= 1. - p[:,0] + tol)) + # vertex node conditions + self.assertArrayNear(p[element.vertexNodes[0],:], onp.array([1.0, 0.0]), 16) + self.assertArrayNear(p[element.vertexNodes[1],:], onp.array([0.0, 1.0]), 16) + self.assertArrayNear(p[element.vertexNodes[2],:], onp.array([0.0, 0.0]), 16) + # face node conditions + k = element.faceNodes + self.assertTrue(np.all(p[k,0] > -tol)) + self.assertTrue(np.all(p[k,1] + p[k,0] - 1. < tol)) + self.assertArrayNear(p[k[0,:],:], onp.array([[1.0, 0.0], [0.5, 0.5], [0.0, 1.0]]), 16) + self.assertArrayNear(p[k[1,:],:], onp.array([[0.0, 1.0], [0.0, 0.5], [0.0, 0.0]]), 16) + self.assertArrayNear(p[k[2,:],:], onp.array([[0.0, 0.0], [0.5, 0.0], [1.0, 0.0]]), 16) + # interior node conditions + self.assertTrue(element.interiorNodes is None) + + def test_edge_shape_kronecker_delta_property(self): + element = Interpolants.make_lagrange_parent_element_1d(self.degree) + shapeAtNodes, _ = Interpolants.shape1d_lagrange(element.degree, element.coordinates, element.coordinates) + nNodes = element.coordinates.shape[0] + self.assertArrayNear(shapeAtNodes, np.identity(nNodes), 14) + + def test_edge_shape_partition_of_unity(self): + element = Interpolants.make_lagrange_parent_element_1d(self.degree) + shapes, _ = Interpolants.shape1d_lagrange(element.degree, element.coordinates, self.qr1d.xigauss) + self.assertArrayNear(np.sum(shapes, axis=0), np.ones(self.nQPts1d), 14) + + def test_edge_shapeGrads_partition_of_unity(self): + element = Interpolants.make_lagrange_parent_element_1d(self.degree) + _, shapeGradients = Interpolants.shape1d_lagrange(element.degree, element.coordinates, self.qr1d.xigauss) + self.assertArrayNear(np.sum(shapeGradients, axis=0), np.zeros(self.nQPts1d), 14) + + def test_edge_interpolation(self): + numRandomEvalPoints = 5 + x = onp.random.uniform(low=0.0, high=1.0, size=(numRandomEvalPoints)) # random points on unit interval + coeffs = onp.random.uniform(low=0.0, high=1.0, size=(3)) # random polynomial coefficients on unit interval + + element = Interpolants.make_lagrange_parent_element_1d(self.degree) + p = element.coordinates + shape, _ = Interpolants.shape1d_lagrange(element.degree, p, x) + fn = onp.polynomial.polynomial.polyval(p, coeffs) + fInterpolated = onp.dot(fn, shape) + + expected = onp.polynomial.polynomial.polyval(x, coeffs) + self.assertArrayNear(expected, np.array(fInterpolated), 14) + + def test_edge_grad_interpolation(self): + numRandomEvalPoints = 5 + x = onp.random.uniform(low=0.0, high=1.0, size=(numRandomEvalPoints)) # random points on unit interval + coeffs = onp.random.uniform(low=0.0, high=1.0, size=(3)) # random polynomial coefficients on unit interval + + element = Interpolants.make_lagrange_parent_element_1d(self.degree) + p = element.coordinates + _, dShape = Interpolants.shape1d_lagrange(element.degree, p, x) + dfn = onp.polynomial.polynomial.polyval(p, coeffs) + dfInterpolated = onp.dot(dfn, dShape) + + dPoly = onp.polynomial.polynomial.polyder(coeffs) + expected = onp.polynomial.polynomial.polyval(x, dPoly) + self.assertArrayNear(expected, np.array(dfInterpolated), 14) + + def test_tri_shape_kronecker_delta_property(self): + element = Interpolants.make_lagrange_parent_element_2d(self.degree) + shapeAtNodes, _ = Interpolants.shape2d_lagrange(element.degree, element.coordinates, element.coordinates) + nNodes = element.coordinates.shape[0] + self.assertArrayNear(shapeAtNodes, np.identity(nNodes), 14) + + def test_tri_shape_partition_of_unity(self): + element = Interpolants.make_lagrange_parent_element_2d(self.degree) + shapes, _ = Interpolants.shape2d_lagrange(element.degree, element.coordinates, self.qr.xigauss) + self.assertArrayNear(np.sum(shapes, axis=1), np.ones(self.nQPts), 14) + + def test_tri_shapeGrads_partition_of_unity(self): + element = Interpolants.make_lagrange_parent_element_2d(self.degree) + _, shapeGradients = Interpolants.shape2d_lagrange(element.degree, element.coordinates, self.qr.xigauss) + self.assertArrayNear(np.sum(shapeGradients, axis=1), np.zeros((self.nQPts, 2)), 14) + + def test_tri_interpolation(self): + numRandomEvalPoints = 5 + element = Interpolants.make_lagrange_parent_element_2d(self.degree) + self.assertCorrectInterpolationOn2DElement(element, numRandomEvalPoints, Interpolants.shape2d_lagrange) + + def test_tri_grad_interpolation(self): + numRandomEvalPoints = 5 + element = Interpolants.make_lagrange_parent_element_2d(self.degree) + self.assertCorrectGradInterpolationOn2DElement(element, numRandomEvalPoints, Interpolants.shape2d_lagrange) + + + if __name__ == '__main__': unittest.main() From 4ff5fb11954fa7f0f0793b02ffa4f643f4efadca Mon Sep 17 00:00:00 2001 From: ralberd Date: Mon, 10 Feb 2025 17:42:43 -0700 Subject: [PATCH 07/19] adds Lagrange elements option to create_higher_order_mesh_from_simplex_mesh --- optimism/Interpolants.py | 6 +++++- optimism/Mesh.py | 17 +++++++++++------ optimism/ReadExodusMesh.py | 9 ++------- optimism/test/test_Interpolants.py | 2 +- 4 files changed, 19 insertions(+), 15 deletions(-) diff --git a/optimism/Interpolants.py b/optimism/Interpolants.py index 934c63fb..4806bb80 100644 --- a/optimism/Interpolants.py +++ b/optimism/Interpolants.py @@ -1,9 +1,13 @@ from jaxtyping import Array, Float, Int from scipy import special +from enum import Enum, auto import numpy as onp import equinox as eqx import jax.numpy as np +class InterpolationType(Enum): + LOBATTO = auto() + LAGRANGE = auto() class ParentElement(eqx.Module): """Finite element on reference domain. @@ -234,7 +238,7 @@ def make_lagrange_parent_element_2d(degree): xn = np.array([[1.0, 0.0], [0.5, 0.5], [0.0, 1.0], [0.5, 0.0], [0.0, 0.5], [0.0, 0.0]]) vertexPoints = np.array([0, 2, 5], dtype=np.int32) facePoints = np.array([[0, 1, 2], [2, 4, 5], [5, 3, 0]], dtype=np.int32) - return ParentElement(LAGRANGE_TRIANGLE_ELEMENT, int(degree), xn, vertexPoints, facePoints, None) + return ParentElement(LAGRANGE_TRIANGLE_ELEMENT, int(degree), xn, vertexPoints, facePoints, np.array([], dtype=np.int32)) def pascal_triangle_monomials(degree): p = [] diff --git a/optimism/Mesh.py b/optimism/Mesh.py index 16061d9a..a520a83f 100644 --- a/optimism/Mesh.py +++ b/optimism/Mesh.py @@ -209,15 +209,20 @@ def create_edges(conns): return edgeConns, edges -def create_higher_order_mesh_from_simplex_mesh(mesh, order, useBubbleElement=False, copyNodeSets=False, createNodeSetsFromSideSets=False): +def create_higher_order_mesh_from_simplex_mesh(mesh, order, interpolationType = Interpolants.InterpolationType.LOBATTO, useBubbleElement=False, copyNodeSets=False, createNodeSetsFromSideSets=False): if order==1: return mesh - parentElement1d = Interpolants.make_parent_element_1d(order) - - if useBubbleElement: - basis = Interpolants.make_parent_element_2d_with_bubble(order) + if interpolationType == Interpolants.InterpolationType.LAGRANGE: + if useBubbleElement: + raise NotImplementedError + parentElement1d = Interpolants.make_lagrange_parent_element_1d(order) + basis = Interpolants.make_lagrange_parent_element_2d(order) else: - basis = Interpolants.make_parent_element_2d(order) + parentElement1d = Interpolants.make_parent_element_1d(order) + if useBubbleElement: + basis = Interpolants.make_parent_element_2d_with_bubble(order) + else: + basis = Interpolants.make_parent_element_2d(order) conns = np.zeros((num_elements(mesh), basis.coordinates.shape[0]), dtype=np.int_) diff --git a/optimism/ReadExodusMesh.py b/optimism/ReadExodusMesh.py index bbbd77d1..a20a94d7 100644 --- a/optimism/ReadExodusMesh.py +++ b/optimism/ReadExodusMesh.py @@ -3,15 +3,10 @@ from optimism import Interpolants import netCDF4 import numpy as onp -from enum import Enum, auto - -class InterpolationType(Enum): - LOBATTO = auto() - LAGRANGE = auto() exodusToNativeTri6NodeOrder = np.array([0, 3, 1, 5, 4, 2]) -def read_exodus_mesh(fileName, interpolationType = InterpolationType.LOBATTO): +def read_exodus_mesh(fileName, interpolationType = Interpolants.InterpolationType.LOBATTO): with netCDF4.Dataset(fileName) as exData: coords = _read_coordinates(exData) conns, blocks = _read_blocks(exData) @@ -30,7 +25,7 @@ def read_exodus_mesh(fileName, interpolationType = InterpolationType.LOBATTO): else: raise - if interpolationType == InterpolationType.LAGRANGE: + if interpolationType == Interpolants.InterpolationType.LAGRANGE: basis1d = Interpolants.make_lagrange_parent_element_1d(degree) basis = Interpolants.make_lagrange_parent_element_2d(degree) else: diff --git a/optimism/test/test_Interpolants.py b/optimism/test/test_Interpolants.py index ff71bda9..78943411 100644 --- a/optimism/test/test_Interpolants.py +++ b/optimism/test/test_Interpolants.py @@ -286,7 +286,7 @@ def test_tri_interpolant_points(self): self.assertArrayNear(p[k[1,:],:], onp.array([[0.0, 1.0], [0.0, 0.5], [0.0, 0.0]]), 16) self.assertArrayNear(p[k[2,:],:], onp.array([[0.0, 0.0], [0.5, 0.0], [1.0, 0.0]]), 16) # interior node conditions - self.assertTrue(element.interiorNodes is None) + self.assertTrue(element.interiorNodes.shape == 0) def test_edge_shape_kronecker_delta_property(self): element = Interpolants.make_lagrange_parent_element_1d(self.degree) From 2013f1c3a05056c59df7be4e4cf8bcef28157c9c Mon Sep 17 00:00:00 2001 From: ralberd Date: Mon, 10 Feb 2025 17:49:29 -0700 Subject: [PATCH 08/19] fix for interpolants test --- optimism/test/test_Interpolants.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/optimism/test/test_Interpolants.py b/optimism/test/test_Interpolants.py index 78943411..bc7a566f 100644 --- a/optimism/test/test_Interpolants.py +++ b/optimism/test/test_Interpolants.py @@ -286,7 +286,7 @@ def test_tri_interpolant_points(self): self.assertArrayNear(p[k[1,:],:], onp.array([[0.0, 1.0], [0.0, 0.5], [0.0, 0.0]]), 16) self.assertArrayNear(p[k[2,:],:], onp.array([[0.0, 0.0], [0.5, 0.0], [1.0, 0.0]]), 16) # interior node conditions - self.assertTrue(element.interiorNodes.shape == 0) + self.assertTrue(element.interiorNodes.shape == (0,)) def test_edge_shape_kronecker_delta_property(self): element = Interpolants.make_lagrange_parent_element_1d(self.degree) From 18b1b6d0a3519c1811a536c2eb090902f186eedf Mon Sep 17 00:00:00 2001 From: ralberd Date: Mon, 10 Feb 2025 18:29:41 -0700 Subject: [PATCH 09/19] fixes issue with TestFixture import --- optimism/test/test_Interpolants.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/optimism/test/test_Interpolants.py b/optimism/test/test_Interpolants.py index bc7a566f..57947f25 100644 --- a/optimism/test/test_Interpolants.py +++ b/optimism/test/test_Interpolants.py @@ -7,7 +7,7 @@ from optimism import Interpolants from optimism import QuadratureRule -import TestFixture +from . import TestFixture tol = 1e-14 From cf81b559347ea7a195b35910d596105778e5de96 Mon Sep 17 00:00:00 2001 From: ralberd Date: Mon, 10 Feb 2025 18:31:34 -0700 Subject: [PATCH 10/19] adds patch tests for quadratic Lagrange elements; removes duplicated code --- optimism/test/test_PatchTest.py | 244 +++++++++++++++++--------------- 1 file changed, 128 insertions(+), 116 deletions(-) diff --git a/optimism/test/test_PatchTest.py b/optimism/test/test_PatchTest.py index 9bd797d7..6c969a05 100644 --- a/optimism/test/test_PatchTest.py +++ b/optimism/test/test_PatchTest.py @@ -23,67 +23,120 @@ 'poisson ratio': nu, 'strain measure': 'linear'} +targetDispGrad = np.array([[0.1, -0.2],[0.4, -0.1]]) -class LinearPatchTestLinearElements(MeshFixture.MeshFixture): - - def setUp(self): - self.Nx = 7 - self.Ny = 7 - xRange = [0.,1.] - yRange = [0.,1.] - - self.targetDispGrad = np.array([[0.1, -0.2],[0.4, -0.1]]) - - self.mesh, self.U = self.create_mesh_and_disp(self.Nx, self.Ny, xRange, yRange, - lambda x : self.targetDispGrad.dot(x)) - quadRule = QuadratureRule.create_quadrature_rule_on_triangle(degree=1) - self.fs = FunctionSpace.construct_function_space(self.mesh, quadRule) +class PatchTestFixture(MeshFixture.MeshFixture): + def initializeLinearPatchTestFromMesh(self, mesh, quadratureDegree): + self.quadRule = QuadratureRule.create_quadrature_rule_on_triangle(degree=quadratureDegree) + self.fs = FunctionSpace.construct_function_space(mesh, self.quadRule) materialModel = MatModel.create_material_model_functions(props) - + 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(self): + def initializeQuadraticPatchTestFromMesh(self, mesh): + alpha = 2.0 + beta = 1.0 + + self.UTarget = jax.vmap( lambda x : targetDispGrad@x + + np.array([alpha * x[0]*x[0], beta * x[1]*x[1]]) )(mesh.coords) + + self.quadRule = QuadratureRule.create_quadrature_rule_on_triangle(degree=2) + self.fs = FunctionSpace.construct_function_space(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 assertLinearPatchTestPasses(self): ebcs = [FunctionSpace.EssentialBC(nodeSet='all_boundary', component=0), FunctionSpace.EssentialBC(nodeSet='all_boundary', component=1)] - dofManager = FunctionSpace.DofManager(self.fs, dim=self.U.shape[1], EssentialBCs=ebcs) - Ubc = dofManager.get_bc_values(self.U) + dofManager = FunctionSpace.DofManager(self.fs, self.UTarget.shape[1], ebcs) + Ubc = dofManager.get_bc_values(self.UTarget) - # Uu is U_unconstrained @jax.jit def objective(Uu): U = dofManager.create_field(Uu, Ubc) return self.compute_energy(U, self.internals) with Timer(name="NewtonSolve"): - Uu, solverSuccess = newton_solve(objective, dofManager.get_unknown_values(self.U)) + Uu, solverSuccess = newton_solve(objective, dofManager.get_unknown_values(self.UTarget)) self.assertTrue(solverSuccess) - self.U = dofManager.create_field(Uu, Ubc) + U = dofManager.create_field(Uu, Ubc) - dispGrads = FunctionSpace.compute_field_gradient(self.fs, self.U) + dispGrads = FunctionSpace.compute_field_gradient(self.fs, U) ne, nqpe = self.fs.vols.shape for dg in dispGrads.reshape(ne*nqpe,2,2): - self.assertArrayNear(dg, self.targetDispGrad, 14) + self.assertArrayNear(dg, targetDispGrad, 14) grad_func = jax.jit(jax.grad(objective)) - Uu = dofManager.get_unknown_values(self.U) + Uu = dofManager.get_unknown_values(U) self.assertNear(np.linalg.norm(grad_func(Uu)), 0.0, 14) + def assertQuadraticPatchTestPasses(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, solverSuccess = newton_solve(objective, 0.0*self.dofManager.get_unknown_values(self.UTarget)) + self.assertTrue(solverSuccess) + + 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) + + +class LinearPatchTestLinearElements(PatchTestFixture): + + def setUp(self): + self.Nx = 7 + self.Ny = 7 + xRange = [0.,1.] + yRange = [0.,1.] + + self.mesh, self.UTarget = self.create_mesh_and_disp(self.Nx, self.Ny, xRange, yRange, + lambda x : targetDispGrad.dot(x)) + self.initializeLinearPatchTestFromMesh(self.mesh, quadratureDegree=1) + + + def test_dirichlet_patch_test(self): + self.assertLinearPatchTestPasses() + def test_neumann_patch_test(self): ebcs = [FunctionSpace.EssentialBC(nodeSet='left', component=0), FunctionSpace.EssentialBC(nodeSet='bottom', component=1)] - self.U = np.zeros(self.U.shape) - dofManager = FunctionSpace.DofManager(self.fs, self.U.shape[1], ebcs) - Ubc = dofManager.get_bc_values(self.U) + self.UTarget = np.zeros(self.UTarget.shape) + dofManager = FunctionSpace.DofManager(self.fs, self.UTarget.shape[1], ebcs) + Ubc = dofManager.get_bc_values(self.UTarget) sigma = np.array([[1.0, 0.0], [0.0, 0.0]]) traction_func = lambda x, n: np.dot(sigma, n) @@ -98,10 +151,10 @@ def objective(Uu): return internalPotential + loadPotential with Timer(name="NewtonSolve"): - Uu, solverSuccess = newton_solve(objective, dofManager.get_unknown_values(self.U)) + Uu, solverSuccess = newton_solve(objective, dofManager.get_unknown_values(self.UTarget)) self.assertTrue(solverSuccess) - self.U = dofManager.create_field(Uu, Ubc) + self.UTarget = dofManager.create_field(Uu, Ubc) # exact solution modulus1 = (1.0 - nu**2)/E @@ -109,11 +162,11 @@ def objective(Uu): UExact = np.column_stack( ((modulus1*sigma[0, 0] + modulus2*sigma[1, 1])*self.mesh.coords[:,0], (modulus2*sigma[0, 0] + modulus1*sigma[1, 1])*self.mesh.coords[:,1]) ) - self.assertArrayNear(self.U, UExact, 14) + self.assertArrayNear(self.UTarget, UExact, 14) -class LinearPatchTestQuadraticElements(MeshFixture.MeshFixture): +class LinearPatchTestQuadraticElements(PatchTestFixture): def setUp(self): self.Nx = 3 @@ -121,53 +174,15 @@ def setUp(self): 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)) + lambda x : targetDispGrad.dot(x)) self.mesh = Mesh.create_higher_order_mesh_from_simplex_mesh(mesh, order=2, createNodeSetsFromSideSets=True) - - self.UTarget = self.mesh.coords@self.targetDispGrad.T + self.initializeLinearPatchTestFromMesh(self.mesh, quadratureDegree=2) + self.UTarget = self.mesh.coords@targetDispGrad.T - 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) - - 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)] - dofManager = FunctionSpace.DofManager(self.fs, self.UTarget.shape[1], ebcs) - Ubc = dofManager.get_bc_values(self.UTarget) - - @jax.jit - def objective(Uu): - U = dofManager.create_field(Uu, Ubc) - return self.compute_energy(U, self.internals) - - with Timer(name="NewtonSolve"): - Uu, solverSuccess = newton_solve(objective, dofManager.get_unknown_values(self.UTarget)) - self.assertTrue(solverSuccess) - - U = dofManager.create_field(Uu, Ubc) - - dispGrads = FunctionSpace.compute_field_gradient(self.fs, U) - ne, nqpe = self.fs.vols.shape - for dg in dispGrads.reshape(ne*nqpe,2,2): - self.assertArrayNear(dg, self.targetDispGrad, 14) - - grad_func = jax.jit(jax.grad(objective)) - Uu = dofManager.get_unknown_values(U) - self.assertNear(np.linalg.norm(grad_func(Uu)), 0.0, 14) + self.assertLinearPatchTestPasses() def test_dirichlet_patch_test_with_quadratic_elements_and_constant_jac_projection(self): @@ -197,14 +212,14 @@ def objective(Uu): dispGrads = FunctionSpace.compute_field_gradient(self.fs, U, modify_grad) ne, nqpe = self.fs.vols.shape for dg in dispGrads.reshape(ne*nqpe,*dispGrads.shape[2:]): - self.assertArrayNear(dg[:2,:2], self.targetDispGrad, 14) + self.assertArrayNear(dg[:2,:2], targetDispGrad, 14) grad_func = jax.jit(jax.grad(objective)) Uu = dofManager.get_unknown_values(U) self.assertNear(np.linalg.norm(grad_func(Uu)), 0.0, 14) -class QuadraticPatchTestQuadraticElements(MeshFixture.MeshFixture): +class QuadraticPatchTestQuadraticElements(PatchTestFixture): def setUp(self): self.Nx = 3 @@ -212,57 +227,54 @@ def setUp(self): 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)) + lambda x : targetDispGrad.dot(x)) self.mesh = Mesh.create_higher_order_mesh_from_simplex_mesh(mesh, order=2, createNodeSetsFromSideSets=True) + self.initializeQuadraticPatchTestFromMesh(self.mesh) - 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) + def test_dirichlet_patch_test_with_quadratic_elements(self): + self.assertQuadraticPatchTestPasses() - 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) +class LinearPatchTestQuadraticLagrangeElements(PatchTestFixture): + + def setUp(self): + self.Nx = 3 + self.Ny = 3 + xRange = [0.,1.] + yRange = [0.,1.] - 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)) + mesh, _ = self.create_mesh_and_disp(self.Nx, self.Ny, xRange, yRange, + lambda x : targetDispGrad.dot(x)) + self.mesh = Mesh.create_higher_order_mesh_from_simplex_mesh(mesh, order=2, + interpolationType=Interpolants.InterpolationType.LAGRANGE, + createNodeSetsFromSideSets=True) + self.initializeLinearPatchTestFromMesh(self.mesh, quadratureDegree=2) + self.UTarget = self.mesh.coords@targetDispGrad.T + + def test_dirichlet_patch_test_with_quadratic_elements(self): + self.assertLinearPatchTestPasses() - 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, solverSuccess = newton_solve(objective, 0.0*self.dofManager.get_unknown_values(self.UTarget)) - self.assertTrue(solverSuccess) - U = self.dofManager.create_field(Uu, Ubc) +class QuadraticPatchTestQuadraticLagrangeElements(PatchTestFixture): + + def setUp(self): + self.Nx = 3 + self.Ny = 3 + xRange = [0.,1.] + yRange = [0.,1.] + + mesh, _ = self.create_mesh_and_disp(self.Nx, self.Ny, xRange, yRange, + lambda x : targetDispGrad.dot(x)) - self.assertArrayNear(U, self.UTarget, 13) + self.mesh = Mesh.create_higher_order_mesh_from_simplex_mesh(mesh, order=2, + interpolationType=Interpolants.InterpolationType.LAGRANGE, + createNodeSetsFromSideSets=True) + self.initializeQuadraticPatchTestFromMesh(self.mesh) - 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) + def test_dirichlet_patch_test_with_quadratic_elements(self): + self.assertQuadraticPatchTestPasses() if __name__ == '__main__': From 631a0c15e2a8d9f875a25f6e92d2935a67f078ca Mon Sep 17 00:00:00 2001 From: ralberd Date: Thu, 13 Feb 2025 17:47:36 -0700 Subject: [PATCH 11/19] implements shape function second derivatives for quadratic lagrange element --- optimism/Interpolants.py | 61 ++++++++++++++++++++++++++++++ optimism/test/test_Interpolants.py | 39 +++++++++++++++++++ 2 files changed, 100 insertions(+) diff --git a/optimism/Interpolants.py b/optimism/Interpolants.py index 4806bb80..d558df74 100644 --- a/optimism/Interpolants.py +++ b/optimism/Interpolants.py @@ -374,6 +374,67 @@ def shape2d_lagrange(degree, nodalPoints, evaluationPoints): return ShapeFunctions(shape, dshape) +def shape2d_lagrange_second_derivatives(degree, nodalPoints, evaluationPoints): + """Evaluate second derivatives of Lagrange shape functions at points in the parent element. + Only implemented for second degree + + Shape of returned second derivatives is ``(nPts, nNodes, nTerms)``, where ``nTerms`` + is the number of derivative terms (3). These terms are organized as + [ (d^2 N)/(dr^2), (d^2 N)/(ds^2), (d^2 N)/(dr ds) ] + where r and s are the local triangle coordinates. + + Following shape2d_lagrange above, the shape functions are + + N0 = r * (2.0 * r - 1.0) = 2.0*r*r - r + N1 = 4.0 * r * s = 4.0 * r * s + N2 = s * (2.0 * s - 1.0) = 2.0*s*s - s + N3 = 4.0 * r * t = 4.0*(r - r*s - r*r) + N4 = 4.0 * s * t = 4.0*(s - r*s - s*s) + N5 = t * (2.0 * t - 1.0) = 1.0 - 3.0*(r + s) + 4.0*r*s + 2.0*(r*r + s*s) + + Reference: + T. Hughes. "The Finite Element Method" + Appendix 3.I + """ + if degree != 2: + raise NotImplementedError + + numEvalPoints = evaluationPoints.shape[0] + r = evaluationPoints[:,0] + s = evaluationPoints[:,1] + # t = 1.0 - r - s + + d2shape0_dr2 = 4.0 * np.ones(numEvalPoints) + d2shape0_ds2 = np.zeros(numEvalPoints) + d2shape0_drds = np.zeros(numEvalPoints) + + d2shape1_dr2 = np.zeros(numEvalPoints) + d2shape1_ds2 = np.zeros(numEvalPoints) + d2shape1_drds = 4.0 * np.ones(numEvalPoints) + + d2shape2_dr2 = np.zeros(numEvalPoints) + d2shape2_ds2 = 4.0 * np.ones(numEvalPoints) + d2shape2_drds = np.zeros(numEvalPoints) + + d2shape3_dr2 = -8.0 * np.ones(numEvalPoints) + d2shape3_ds2 = np.zeros(numEvalPoints) + d2shape3_drds = -4.0 * np.ones(numEvalPoints) + + d2shape4_dr2 = np.zeros(numEvalPoints) + d2shape4_ds2 = -8.0 * np.ones(numEvalPoints) + d2shape4_drds = -4.0 * np.ones(numEvalPoints) + + d2shape5_dr2 = 4.0 * np.ones(numEvalPoints) + d2shape5_ds2 = 4.0 * np.ones(numEvalPoints) + d2shape5_drds = 4.0 * np.ones(numEvalPoints) + + d2shape_dr2 = np.stack((d2shape0_dr2, d2shape1_dr2, d2shape2_dr2, d2shape3_dr2, d2shape4_dr2, d2shape5_dr2)).T + d2shape_ds2 = np.stack((d2shape0_ds2, d2shape1_ds2, d2shape2_ds2, d2shape3_ds2, d2shape4_ds2, d2shape5_ds2)).T + d2shape_drds = np.stack((d2shape0_drds, d2shape1_drds, d2shape2_drds, d2shape3_drds, d2shape4_drds, d2shape5_drds)).T + + return np.stack((d2shape_dr2, d2shape_ds2, d2shape_drds), axis=2) + + def compute_shapes(parentElement, evaluationPoints): if parentElement.elementType == LINE_ELEMENT: return shape1d(parentElement.degree, parentElement.coordinates, evaluationPoints) diff --git a/optimism/test/test_Interpolants.py b/optimism/test/test_Interpolants.py index 57947f25..7a36316c 100644 --- a/optimism/test/test_Interpolants.py +++ b/optimism/test/test_Interpolants.py @@ -256,6 +256,7 @@ def test_degree_other_than_two_throws(self): self.assertRaises(NotImplementedError, Interpolants.make_lagrange_parent_element_2d, degree) self.assertRaises(NotImplementedError, Interpolants.shape1d_lagrange, degree, dummyCoords, self.qr1d) self.assertRaises(NotImplementedError, Interpolants.shape2d_lagrange, degree, dummyCoords, self.qr) + self.assertRaises(NotImplementedError, Interpolants.shape2d_lagrange_second_derivatives, degree, dummyCoords, self.qr) def test_1D_interpolant_points(self): element = Interpolants.make_lagrange_parent_element_1d(self.degree) @@ -358,6 +359,44 @@ def test_tri_grad_interpolation(self): numRandomEvalPoints = 5 element = Interpolants.make_lagrange_parent_element_2d(self.degree) self.assertCorrectGradInterpolationOn2DElement(element, numRandomEvalPoints, Interpolants.shape2d_lagrange) + + def test_tri_second_grad_partition_of_unity(self): + element = Interpolants.make_lagrange_parent_element_2d(self.degree) + shapeSecondGrad = Interpolants.shape2d_lagrange_second_derivatives(element.degree, element.coordinates, self.qr.xigauss) + self.assertArrayNear(np.sum(shapeSecondGrad, axis=1), np.zeros((self.nQPts, 3)), 14) + + def test_tri_second_grad_interpolation(self): + numRandomEvalPoints = 5 + element = Interpolants.make_lagrange_parent_element_2d(self.degree) + + degree = element.degree + x = generate_random_points_in_triangle(numRandomEvalPoints) + + poly = onp.array([[1, 1, 2], [1, 1, 0], [1, 0, 0]]) + + shapeSecondGrad = Interpolants.shape2d_lagrange_second_derivatives(degree, element.coordinates, x) + fn = onp.polynomial.polynomial.polyval2d(element.coordinates[:,0], + element.coordinates[:,1], + poly) + d2fInterpolated = onp.einsum('qai,a->qi', shapeSecondGrad, fn) + + # 00 derivative + DPoly = onp.polynomial.polynomial.polyder(poly, 1, scl=1, axis=0) + D2Poly = onp.polynomial.polynomial.polyder(DPoly, 1, scl=1, axis=0) + expected00 = onp.polynomial.polynomial.polyval2d(x[:,0], x[:,1], D2Poly) + self.assertArrayNear(expected00, d2fInterpolated[:,0], 14) + + # 11 derivative + DPoly = onp.polynomial.polynomial.polyder(poly, 1, scl=1, axis=1) + D2Poly = onp.polynomial.polynomial.polyder(DPoly, 1, scl=1, axis=1) + expected11 = onp.polynomial.polynomial.polyval2d(x[:,0], x[:,1], D2Poly) + self.assertArrayNear(expected11, d2fInterpolated[:,1], 14) + + # 01 derivative + DPoly = onp.polynomial.polynomial.polyder(poly, 1, scl=1, axis=0) + D2Poly = onp.polynomial.polynomial.polyder(DPoly, 1, scl=1, axis=1) + expected01 = onp.polynomial.polynomial.polyval2d(x[:,0], x[:,1], D2Poly) + self.assertArrayNear(expected01, d2fInterpolated[:,2], 14) From 3b13bf09fdaf601674b46bcbe0c9bb9b11cf1c46 Mon Sep 17 00:00:00 2001 From: ralberd Date: Fri, 14 Feb 2025 15:38:50 -0700 Subject: [PATCH 12/19] adds functions for computing field hessians to FunctionSpace --- optimism/FunctionSpace.py | 48 +++++++++++++++++++++++++++-- optimism/test/test_FunctionSpace.py | 48 +++++++++++++++++++++++++++++ 2 files changed, 94 insertions(+), 2 deletions(-) diff --git a/optimism/FunctionSpace.py b/optimism/FunctionSpace.py index d22dff3b..113715b6 100644 --- a/optimism/FunctionSpace.py +++ b/optimism/FunctionSpace.py @@ -22,6 +22,7 @@ class FunctionSpace(eqx.Module): elements in the mesh, ``nqpe`` is the number of quadrature points per element, ``npe`` is the number of nodes per element, and ``nd`` is the spatial dimension of the domain. + For shape Hessians, nh is the number of hessian entries, set to 3 Attributes: shapes: Shape function values on each element, shape (ne, nqpe, npe) @@ -31,6 +32,9 @@ class FunctionSpace(eqx.Module): element in the domain. Shape (ne, nqpe). shapeGrads: Derivatives of the shape functions with respect to the spatial coordinates of the domain. Shape (ne, nqpe, npe, nd). + shapeHessians: Second derivatives of the shape functions with respect to the + spatial coordinates of the domain. Shape (ne, nqpe, npe, nh). + These are only implemented for quadratic Lagrange elements. mesh: The ``Mesh`` object of the domain. quadratureRule: The ``QuadratureRule`` on which to sample the shape functions. @@ -40,6 +44,7 @@ class FunctionSpace(eqx.Module): shapes: Float[Array, "ne nqpe nn"] vols: Float[Array, "ne nqpe"] shapeGrads: Float[Array, "ne nqpe nn nd"] + shapeHessians: Float[Array, "ne nqpe nn nh"] mesh: Mesh.Mesh quadratureRule: QuadratureRule.QuadratureRule isAxisymmetric: bool @@ -112,16 +117,40 @@ def construct_function_space_from_parent_element(mesh, shapeOnRef, quadratureRul isAxisymmetric = True vols = jax.vmap(el_vols, (None, 0, None, 0, None))(mesh.coords, mesh.conns, mesh.parentElement, shapes, quadratureRule.wgauss) - return FunctionSpace(shapes, vols, shapeGrads, mesh, quadratureRule, isAxisymmetric) + if mesh.parentElement.elementType == Interpolants.LAGRANGE_TRIANGLE_ELEMENT: + shapeOnRefHessians = Interpolants.shape2d_lagrange_second_derivatives(mesh.parentElement.degree, mesh.parentElement.coordinates, quadratureRule.xigauss) + shapeHessians = jax.vmap(map_element_shape_hessians, (None, 0, None, None, None))(mesh.coords, mesh.conns, mesh.parentElement, shapeOnRef.gradients, shapeOnRefHessians) + return FunctionSpace(shapes, vols, shapeGrads, shapeHessians, mesh, quadratureRule, isAxisymmetric) + else: + return FunctionSpace(shapes, vols, shapeGrads, None, mesh, quadratureRule, isAxisymmetric) def map_element_shape_grads(coordField, nodeOrdinals, parentElement, shapeGradients): Xn = coordField.take(nodeOrdinals,0) v = Xn[parentElement.vertexNodes] - J = np.column_stack((v[0] - v[2], v[1] - v[2])) + J = np.column_stack((v[0] - v[2], v[1] - v[2])) # assumes simplex element return jax.vmap(lambda dN: solve(J.T, dN.T).T)(shapeGradients) +def map_element_shape_hessians(coordField, nodeOrdinals, parentElement, shapeGradients, shapeHessians): + Xn = coordField.take(nodeOrdinals,0) + + dNdX = map_element_shape_grads(coordField, nodeOrdinals, parentElement, shapeGradients) + + def quad_point_shape_hessian(shapeGrads, dNdX, shapeHessians): + dXdXi = np.tensordot(Xn, shapeGrads, axes=[0,0]) + J = np.array([[dXdXi[0,0]**2, dXdXi[1,0]**2, 2.0*dXdXi[1,0]*dXdXi[0,0]], + [dXdXi[0,1]**2, dXdXi[1,1]**2, 2.0*dXdXi[1,1]*dXdXi[0,1]], + [dXdXi[0,0]*dXdXi[0,1], dXdXi[1,0]*dXdXi[1,1], dXdXi[1,0]*dXdXi[0,1] + dXdXi[0,0]*dXdXi[1,1]]]) + + d2X_dXi2 = np.tensordot(Xn, shapeHessians, axes=[0,0]) + b = np.array([shapeHessians[:,0].T - dNdX[:,0].T * d2X_dXi2[0,0] - dNdX[:,1].T * d2X_dXi2[1,0], + shapeHessians[:,1].T - dNdX[:,0].T * d2X_dXi2[0,1] - dNdX[:,1].T * d2X_dXi2[1,1], + shapeHessians[:,2].T - dNdX[:,0].T * d2X_dXi2[0,2] - dNdX[:,1].T * d2X_dXi2[1,2]]) + return solve(J, b).T + return jax.vmap(quad_point_shape_hessian, (0, 0, 0))(shapeGradients, dNdX, shapeHessians) # vmap over nqpe + + def compute_element_volumes(coordField, nodeOrdinals, parentElement, shapes, weights): Xn = coordField.take(nodeOrdinals,0) v = Xn[parentElement.vertexNodes] @@ -189,6 +218,12 @@ def compute_field_gradient(functionSpace, nodalField, modify_element_gradient=de return jax.vmap(compute_element_field_gradient, (None,None,0,0,0,0,None))(nodalField, functionSpace.mesh.coords, functionSpace.shapes, functionSpace.shapeGrads, functionSpace.vols, functionSpace.mesh.conns, modify_element_gradient) +def compute_field_hessian(functionSpace, nodalField): + if functionSpace.shapeHessians == None: + raise NotImplementedError + return jax.vmap(compute_element_field_hessian, (None,0,0))(nodalField, functionSpace.shapeHessians, functionSpace.mesh.conns) + + def interpolate_to_points(functionSpace, nodalField): return jax.vmap(interpolate_to_element_points, (None, 0, 0))(nodalField, functionSpace.shapes, functionSpace.mesh.conns) @@ -282,6 +317,15 @@ def compute_quadrature_point_field_gradient(u, shapeGrad): return dg +def compute_element_field_hessian(U, elemShapeHessians, elemConnectivity): + elemNodalDisps = U[elemConnectivity] + return jax.vmap(compute_quadrature_point_field_hessian, (None, 0))(elemNodalDisps, elemShapeHessians) + + +def compute_quadrature_point_field_hessian(u, shapeHess): + return np.tensordot(u, shapeHess, axes=[0,0]) + + def interpolate_to_point(elementNodalValues, shape): return np.dot(shape, elementNodalValues) diff --git a/optimism/test/test_FunctionSpace.py b/optimism/test/test_FunctionSpace.py index 7baa8017..4117c715 100644 --- a/optimism/test/test_FunctionSpace.py +++ b/optimism/test/test_FunctionSpace.py @@ -170,6 +170,10 @@ def test_linear_reproducing_multi_point_quadrature(self): self.assertTrue(np.allclose(dispGrads, exact)) + def test_field_hessian_throws(self): + self.assertRaises(NotImplementedError, FunctionSpace.compute_field_hessian, self.fs, self.U) + + def test_integrate_constant_field_multi_point_quadrature(self): integralOfOne = FunctionSpace.integrate_over_block(self.fs, self.U, @@ -284,5 +288,49 @@ def centroid(v): exact = 1.5*self.xRange[1]*self.yRange[1] self.assertAlmostEqual(f, exact, 12) + +class TestFunctionSpaceWithQuadraticLagrangeElementFixture(MeshFixture.MeshFixture): + def setUp(self): + self.quadratureRule = QuadratureRule.create_quadrature_rule_on_triangle(degree=1) + + # mesh + self.Nx = 2 + self.Ny = 2 + self.xRange = [0.,1.] + self.yRange = [0.,1.] + self.targetDispHessian = np.array([[0.1, -0.2, 0.6],[-0.4, 0.5, -0.6]]) + + # u_i = 0.5 G_ijk x_j x_k + initial_disp_func = lambda x : 0.5*self.targetDispHessian.dot(np.array([x[0]*x[0], x[1]*x[1], 2.0*x[0]*x[1]])) + mesh, _ = self.create_mesh_and_disp(self.Nx, self.Ny, self.xRange, self.yRange, + initial_disp_func) + self.mesh = Mesh.create_higher_order_mesh_from_simplex_mesh(mesh, order=2, + interpolationType=Interpolants.InterpolationType.LAGRANGE) + self.U = jax.vmap(initial_disp_func)(self.mesh.coords) + self.fs = FunctionSpace.construct_function_space(self.mesh, self.quadratureRule) + + def test_quadratic_reproducing_gradient_and_hessian(self): + # du_i / dx_p = G_ipk x_k + def compute_expected_displacement_grads(X, elemConnectivity, elemShapeFunctions): + elemNodalX = X[elemConnectivity] + def compute_qp_displacement_grad(NodalX, shapeFuncs): + x = np.dot(shapeFuncs, NodalX) + return np.array([[self.targetDispHessian[0,0]*x[0] + self.targetDispHessian[0,2]*x[1], self.targetDispHessian[0,2]*x[0] + self.targetDispHessian[0,1]*x[1]], + [self.targetDispHessian[1,0]*x[0] + self.targetDispHessian[1,2]*x[1], self.targetDispHessian[1,2]*x[0] + self.targetDispHessian[1,1]*x[1]]]) + return jax.vmap(compute_qp_displacement_grad, (None, 0))(elemNodalX, elemShapeFunctions) + + dispGrads = FunctionSpace.compute_field_gradient(self.fs, self.U) + shapeOnRef = Interpolants.compute_shapes(self.mesh.parentElement, self.quadratureRule.xigauss) + exact = jax.vmap(compute_expected_displacement_grads, (None, 0, None))(self.mesh.coords, self.mesh.conns, shapeOnRef.values) + self.assertTrue(np.allclose(dispGrads, exact)) + + dispHessians = FunctionSpace.compute_field_hessian(self.fs, self.U) + nElements = Mesh.num_elements(self.mesh) + npts = self.quadratureRule.xigauss.shape[0] + exact = np.tile(self.targetDispHessian, (nElements, npts, 1, 1)) + self.assertTrue(np.allclose(dispHessians, exact)) + + + if __name__ == '__main__': unittest.main() From b815778e5f278a6f964a62584a6c5f14e085ad1e Mon Sep 17 00:00:00 2001 From: ralberd Date: Fri, 14 Feb 2025 15:53:19 -0700 Subject: [PATCH 13/19] adds hessian argument for direct FunctionSpace construction to construct_function_space_for_adjoint --- optimism/inverse/AdjointFunctionSpace.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/optimism/inverse/AdjointFunctionSpace.py b/optimism/inverse/AdjointFunctionSpace.py index eeba2584..80713504 100644 --- a/optimism/inverse/AdjointFunctionSpace.py +++ b/optimism/inverse/AdjointFunctionSpace.py @@ -25,4 +25,4 @@ def construct_function_space_for_adjoint(coords, shapeOnRef, mesh, quadratureRul parentElement=mesh.parentElement, parentElement1d=mesh.parentElement1d, blocks=mesh.blocks, nodeSets=mesh.nodeSets, sideSets=mesh.sideSets) - return FunctionSpace.FunctionSpace(shapes, vols, shapeGrads, mesh, quadratureRule, isAxisymmetric) + return FunctionSpace.FunctionSpace(shapes, vols, shapeGrads, None, mesh, quadratureRule, isAxisymmetric) From 5f23c49300536204f83bd2f30233c727fcec70be Mon Sep 17 00:00:00 2001 From: ralberd Date: Mon, 17 Feb 2025 14:29:05 -0700 Subject: [PATCH 14/19] concatenates shape gradients and hessians into same container in order to facilitate transfer to material points. Requires ensuring material points only take correct components for displacement gradients. Requires changes to element modification functions (only updated plane strain) --- optimism/FunctionSpace.py | 24 ++---------------------- optimism/Mechanics.py | 4 +++- optimism/inverse/AdjointFunctionSpace.py | 2 +- optimism/material/LinearElastic.py | 2 +- optimism/test/test_FunctionSpace.py | 15 ++++++++------- optimism/test/test_PatchTest.py | 2 +- 6 files changed, 16 insertions(+), 33 deletions(-) diff --git a/optimism/FunctionSpace.py b/optimism/FunctionSpace.py index 113715b6..9c2752fc 100644 --- a/optimism/FunctionSpace.py +++ b/optimism/FunctionSpace.py @@ -22,7 +22,6 @@ class FunctionSpace(eqx.Module): elements in the mesh, ``nqpe`` is the number of quadrature points per element, ``npe`` is the number of nodes per element, and ``nd`` is the spatial dimension of the domain. - For shape Hessians, nh is the number of hessian entries, set to 3 Attributes: shapes: Shape function values on each element, shape (ne, nqpe, npe) @@ -32,9 +31,6 @@ class FunctionSpace(eqx.Module): element in the domain. Shape (ne, nqpe). shapeGrads: Derivatives of the shape functions with respect to the spatial coordinates of the domain. Shape (ne, nqpe, npe, nd). - shapeHessians: Second derivatives of the shape functions with respect to the - spatial coordinates of the domain. Shape (ne, nqpe, npe, nh). - These are only implemented for quadratic Lagrange elements. mesh: The ``Mesh`` object of the domain. quadratureRule: The ``QuadratureRule`` on which to sample the shape functions. @@ -44,7 +40,6 @@ class FunctionSpace(eqx.Module): shapes: Float[Array, "ne nqpe nn"] vols: Float[Array, "ne nqpe"] shapeGrads: Float[Array, "ne nqpe nn nd"] - shapeHessians: Float[Array, "ne nqpe nn nh"] mesh: Mesh.Mesh quadratureRule: QuadratureRule.QuadratureRule isAxisymmetric: bool @@ -120,9 +115,9 @@ def construct_function_space_from_parent_element(mesh, shapeOnRef, quadratureRul if mesh.parentElement.elementType == Interpolants.LAGRANGE_TRIANGLE_ELEMENT: shapeOnRefHessians = Interpolants.shape2d_lagrange_second_derivatives(mesh.parentElement.degree, mesh.parentElement.coordinates, quadratureRule.xigauss) shapeHessians = jax.vmap(map_element_shape_hessians, (None, 0, None, None, None))(mesh.coords, mesh.conns, mesh.parentElement, shapeOnRef.gradients, shapeOnRefHessians) - return FunctionSpace(shapes, vols, shapeGrads, shapeHessians, mesh, quadratureRule, isAxisymmetric) + return FunctionSpace(shapes, vols, np.concatenate((shapeGrads, shapeHessians), axis=-1), mesh, quadratureRule, isAxisymmetric) else: - return FunctionSpace(shapes, vols, shapeGrads, None, mesh, quadratureRule, isAxisymmetric) + return FunctionSpace(shapes, vols, shapeGrads, mesh, quadratureRule, isAxisymmetric) def map_element_shape_grads(coordField, nodeOrdinals, parentElement, shapeGradients): @@ -218,12 +213,6 @@ def compute_field_gradient(functionSpace, nodalField, modify_element_gradient=de return jax.vmap(compute_element_field_gradient, (None,None,0,0,0,0,None))(nodalField, functionSpace.mesh.coords, functionSpace.shapes, functionSpace.shapeGrads, functionSpace.vols, functionSpace.mesh.conns, modify_element_gradient) -def compute_field_hessian(functionSpace, nodalField): - if functionSpace.shapeHessians == None: - raise NotImplementedError - return jax.vmap(compute_element_field_hessian, (None,0,0))(nodalField, functionSpace.shapeHessians, functionSpace.mesh.conns) - - def interpolate_to_points(functionSpace, nodalField): return jax.vmap(interpolate_to_element_points, (None, 0, 0))(nodalField, functionSpace.shapes, functionSpace.mesh.conns) @@ -317,15 +306,6 @@ def compute_quadrature_point_field_gradient(u, shapeGrad): return dg -def compute_element_field_hessian(U, elemShapeHessians, elemConnectivity): - elemNodalDisps = U[elemConnectivity] - return jax.vmap(compute_quadrature_point_field_hessian, (None, 0))(elemNodalDisps, elemShapeHessians) - - -def compute_quadrature_point_field_hessian(u, shapeHess): - return np.tensordot(u, shapeHess, axes=[0,0]) - - def interpolate_to_point(elementNodalValues, shape): return np.dot(shape, elementNodalValues) diff --git a/optimism/Mechanics.py b/optimism/Mechanics.py index 619969d7..9a6cf25e 100644 --- a/optimism/Mechanics.py +++ b/optimism/Mechanics.py @@ -52,7 +52,9 @@ class NewmarkParameters(eqx.Module): def plane_strain_gradient_transformation(elemDispGrads, elemShapes, elemVols, elemNodalDisps, elemNodalCoords): - return vmap(tensor_2D_to_3D)(elemDispGrads) + def map_2D_tensor_to_3D(H): + return np.zeros((3,H.shape[1]+1)).at[ 0:H.shape[0], 0:2 ].set(H[:,0:2]).at[ 0:H.shape[0], 3: ].set(H[:,2:]) + return vmap(map_2D_tensor_to_3D)(elemDispGrads) def volume_average_J_gradient_transformation(elemDispGrads, elemVols, pShapes): diff --git a/optimism/inverse/AdjointFunctionSpace.py b/optimism/inverse/AdjointFunctionSpace.py index 80713504..eeba2584 100644 --- a/optimism/inverse/AdjointFunctionSpace.py +++ b/optimism/inverse/AdjointFunctionSpace.py @@ -25,4 +25,4 @@ def construct_function_space_for_adjoint(coords, shapeOnRef, mesh, quadratureRul parentElement=mesh.parentElement, parentElement1d=mesh.parentElement1d, blocks=mesh.blocks, nodeSets=mesh.nodeSets, sideSets=mesh.sideSets) - return FunctionSpace.FunctionSpace(shapes, vols, shapeGrads, None, mesh, quadratureRule, isAxisymmetric) + return FunctionSpace.FunctionSpace(shapes, vols, shapeGrads, mesh, quadratureRule, isAxisymmetric) diff --git a/optimism/material/LinearElastic.py b/optimism/material/LinearElastic.py index 290cbb7d..83cc1a7b 100644 --- a/optimism/material/LinearElastic.py +++ b/optimism/material/LinearElastic.py @@ -31,7 +31,7 @@ def create_material_model_functions(properties): def strain_energy(dispGrad, internalVars, dt): del internalVars del dt - strain = _strain(dispGrad) + strain = _strain(dispGrad[0:3,0:3]) return _linear_elastic_energy_density(strain, props) density = properties.get('density') diff --git a/optimism/test/test_FunctionSpace.py b/optimism/test/test_FunctionSpace.py index 4117c715..27a04125 100644 --- a/optimism/test/test_FunctionSpace.py +++ b/optimism/test/test_FunctionSpace.py @@ -170,10 +170,6 @@ def test_linear_reproducing_multi_point_quadrature(self): self.assertTrue(np.allclose(dispGrads, exact)) - def test_field_hessian_throws(self): - self.assertRaises(NotImplementedError, FunctionSpace.compute_field_hessian, self.fs, self.U) - - def test_integrate_constant_field_multi_point_quadrature(self): integralOfOne = FunctionSpace.integrate_over_block(self.fs, self.U, @@ -291,7 +287,7 @@ def centroid(v): class TestFunctionSpaceWithQuadraticLagrangeElementFixture(MeshFixture.MeshFixture): def setUp(self): - self.quadratureRule = QuadratureRule.create_quadrature_rule_on_triangle(degree=1) + self.quadratureRule = QuadratureRule.create_quadrature_rule_on_triangle(degree=2) # mesh self.Nx = 2 @@ -319,12 +315,17 @@ def compute_qp_displacement_grad(NodalX, shapeFuncs): [self.targetDispHessian[1,0]*x[0] + self.targetDispHessian[1,2]*x[1], self.targetDispHessian[1,2]*x[0] + self.targetDispHessian[1,1]*x[1]]]) return jax.vmap(compute_qp_displacement_grad, (None, 0))(elemNodalX, elemShapeFunctions) - dispGrads = FunctionSpace.compute_field_gradient(self.fs, self.U) + from optimism import Mechanics + dispGradsAndHessians = FunctionSpace.compute_field_gradient(self.fs, self.U, modify_element_gradient=Mechanics.plane_strain_gradient_transformation) + + # check displacement gradients + dispGrads = dispGradsAndHessians[:,:,0:2,0:2] shapeOnRef = Interpolants.compute_shapes(self.mesh.parentElement, self.quadratureRule.xigauss) exact = jax.vmap(compute_expected_displacement_grads, (None, 0, None))(self.mesh.coords, self.mesh.conns, shapeOnRef.values) self.assertTrue(np.allclose(dispGrads, exact)) - dispHessians = FunctionSpace.compute_field_hessian(self.fs, self.U) + # check displacement hessians + dispHessians = dispGradsAndHessians[:,:,0:2,3:] nElements = Mesh.num_elements(self.mesh) npts = self.quadratureRule.xigauss.shape[0] exact = np.tile(self.targetDispHessian, (nElements, npts, 1, 1)) diff --git a/optimism/test/test_PatchTest.py b/optimism/test/test_PatchTest.py index 6c969a05..421db4f1 100644 --- a/optimism/test/test_PatchTest.py +++ b/optimism/test/test_PatchTest.py @@ -78,7 +78,7 @@ def objective(Uu): dispGrads = FunctionSpace.compute_field_gradient(self.fs, U) ne, nqpe = self.fs.vols.shape - for dg in dispGrads.reshape(ne*nqpe,2,2): + for dg in dispGrads[:,:,0:2,0:2].reshape(ne*nqpe,2,2): self.assertArrayNear(dg, targetDispGrad, 14) grad_func = jax.jit(jax.grad(objective)) From a2a2a2ade5e9d4952b04ce2bdcf0c5ced9f2d3a4 Mon Sep 17 00:00:00 2001 From: ralberd Date: Mon, 17 Feb 2025 15:08:57 -0700 Subject: [PATCH 15/19] testing Teran invertible Neo Hookean in third medium --- optimism/material/ThirdMediumNeoHookean.py | 42 +++++++++++++++++++++- optimism/material/test/test_ThirdMedium.py | 13 +++++-- 2 files changed, 51 insertions(+), 4 deletions(-) diff --git a/optimism/material/ThirdMediumNeoHookean.py b/optimism/material/ThirdMediumNeoHookean.py index 257b3404..2d128b34 100644 --- a/optimism/material/ThirdMediumNeoHookean.py +++ b/optimism/material/ThirdMediumNeoHookean.py @@ -1,6 +1,9 @@ +import jax import jax.numpy as np from optimism.material.MaterialModel import MaterialModel +from optimism import TensorMath +from optimism import ScalarRootFind # props PROPS_MU = 0 @@ -17,7 +20,8 @@ def strain_energy(dispGrad, internalVars, dt): del internalVars del dt # return neo_hookean_energy_density(dispGrad, props) - return stabilized_neo_hookean_energy_density(dispGrad, props) + # return stabilized_neo_hookean_energy_density(dispGrad, props) + return teran_invertible_energy_density(dispGrad, props) def compute_state_new(dispGrad, internalVars, dt): del dispGrad @@ -60,6 +64,42 @@ def stabilized_neo_hookean_energy_density(dispGrad, props): Wvol = props[PROPS_LAMBDA]/2.0 * (J - alpha)**2 return Wiso + Wvol +def teran_invertible_energy_density(dispGrad, props): + J_min = 0.01 + J = TensorMath.det(dispGrad + np.identity(3)) + return jax.lax.cond(J >= J_min, _standard_energy, _stabilizing_energy_extension, dispGrad, props, J_min) + +def _standard_energy(gradDisp, props, J_min): + I1m3 = 2*np.trace(gradDisp) + np.tensordot(gradDisp, gradDisp) + Jm1 = TensorMath.detpIm1(gradDisp) + logJ = np.log1p(Jm1) + return 0.5*props[PROPS_MU]*I1m3 - props[PROPS_MU]*logJ + 0.5*props[PROPS_LAMBDA]*logJ**2 + +def _energy_from_principal_stretches(stretches, props): + J = stretches[0]*stretches[1]*stretches[2] + return 0.5*props[PROPS_MU]*(np.sum(stretches**2) - 3) - props[PROPS_MU]*np.log(J) + 0.5*props[PROPS_LAMBDA]*np.log(J)**2 + +def _stabilizing_energy_extension(dispGrad, props, J_min): + F = dispGrad + np.eye(3) + C = F.T@F + stretches_squared, _ = TensorMath.eigen_sym33_unit(C) + stretches = np.sqrt(stretches_squared) + stretches = stretches.at[0].set(np.where(np.linalg.det(F) < 0, -stretches[0], stretches[0])) + ee = stretches - 1 + I1 = ee[0] + ee[1] + ee[2] + I2 = ee[0]*ee[1] + ee[1]*ee[2] + ee[2]*ee[0] + I3 = ee[0]*ee[1]*ee[2] + solver_settings = ScalarRootFind.get_settings(x_tol=1e-8) + s, _ = ScalarRootFind.find_root(lambda x: I3*x**3 + I2*x**2 + I1*x + (1 - J_min), 0.5, np.array([0.0, 1.0]), solver_settings) + q = 1 + s*ee # series expansion point + h = np.linalg.norm(stretches - q) + v = h*ee/np.linalg.norm(ee) # h*u in the paper + W = lambda x: _energy_from_principal_stretches(x, props) + psi0, psi1 = jax.jvp(W, (q,), (v,)) + hess = jax.hessian(W)(q) + psi2 = 0.5*np.dot(v, hess.dot(v)) + return psi0 + psi1 + psi2 + def _compute_volumetric_jacobian(dispGrad): F = dispGrad + np.eye(3) return np.linalg.det(F) \ No newline at end of file diff --git a/optimism/material/test/test_ThirdMedium.py b/optimism/material/test/test_ThirdMedium.py index ec348f70..ed1bd59b 100644 --- a/optimism/material/test/test_ThirdMedium.py +++ b/optimism/material/test/test_ThirdMedium.py @@ -3,6 +3,7 @@ import jax import jax.numpy as np from optimism.material import Neohookean +from optimism.material import ThirdMediumNeoHookean from matplotlib import pyplot as plt @@ -51,6 +52,11 @@ def energy_stable_arruda_boyce(dispGrad, kappa, mu): Wvol = lam/2.0 * (J - alpha)**2 return Wiso + Wvol +def energy_teran_invertible(dispGrad, kappa, mu): + lamda = kappa - (2.0/3.0) * mu + props = np.array([mu, kappa, lamda]) + return ThirdMediumNeoHookean.teran_invertible_energy_density(dispGrad, props) + class ThirdMediumModelFixture(TestFixture): def setUp(self): self.kappa = 100.0 @@ -91,7 +97,7 @@ def test_plot_biaxial_response(self): energies = energies.at[0,n].set(energy_neo_hookean_rokos_volumetric(dispGrad, self.kappa, self.mu)) energies = energies.at[1,n].set(energy_neo_hookean_rokos_isochoric(dispGrad, self.kappa, self.mu)) energies = energies.at[2,n].set(energy_stable_neo_hookean(dispGrad, self.kappa, self.mu)) - energies = energies.at[3,n].set(energy_stable_arruda_boyce(dispGrad, self.kappa, self.mu)) + energies = energies.at[3,n].set(energy_teran_invertible(dispGrad, self.kappa, self.mu)) Jvals = Jvals.at[n].set(np.linalg.det(F)) if plotting: @@ -110,11 +116,11 @@ def test_plot_biaxial_response(self): ax2.set_ylabel('J') # we already handled the x-label with ax1 legends5 = ax2.plot(f_vals, Jvals, linestyle='--', color='k') - ax1.legend(legends1+legends2+legends3+legends4,[ + ax1.legend(legends1+legends2+legends3+legends4+legends5,[ "Rokos Neo Hookean Volumetric term", "Rokos Neo Hookean Isochoric term", "Stable Neohookean energy", - "Stable Arruda Boyce energy", + "Teran invertible Neo Hookean", "Volume change J" ], loc=0, frameon=True) @@ -124,5 +130,6 @@ def test_plot_biaxial_response(self): plt.savefig('energy_vs_compression.png') + if __name__ == "__main__": unittest.main() From 33e61ce98c5a083ce84ee09e16f69467a22c1491 Mon Sep 17 00:00:00 2001 From: ralberd Date: Wed, 19 Feb 2025 13:57:51 -0700 Subject: [PATCH 16/19] small fix for change to displacement gradient storage --- optimism/Mechanics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/optimism/Mechanics.py b/optimism/Mechanics.py index 9a6cf25e..5b487ab9 100644 --- a/optimism/Mechanics.py +++ b/optimism/Mechanics.py @@ -254,7 +254,7 @@ def compute_output_energy_densities_and_stresses(U, stateVariables, dt=0.0): elemIds = fs.mesh.blocks[blockKey] blockEnergyDensities, blockStresses = FunctionSpace.evaluate_on_block(fs, U, stateVariables, dt, output_constitutive, elemIds, modify_element_gradient=modify_element_gradient) energy_densities = energy_densities.at[elemIds].set(blockEnergyDensities) - stresses = stresses.at[elemIds].set(blockStresses) + stresses = stresses.at[elemIds].set(blockStresses[:,:,0:3,0:3]) return energy_densities, stresses def compute_initial_state(): From 56a8aab652eae459d088aa2bff9459ea66122363 Mon Sep 17 00:00:00 2001 From: ralberd Date: Wed, 19 Feb 2025 13:58:23 -0700 Subject: [PATCH 17/19] fixes adjoint function space to match changes in FunctionSpace --- optimism/inverse/AdjointFunctionSpace.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/optimism/inverse/AdjointFunctionSpace.py b/optimism/inverse/AdjointFunctionSpace.py index eeba2584..c7a43c2a 100644 --- a/optimism/inverse/AdjointFunctionSpace.py +++ b/optimism/inverse/AdjointFunctionSpace.py @@ -3,8 +3,9 @@ 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 optimism.FunctionSpace import map_element_shape_grads, map_element_shape_hessians from jax import vmap +import jax.numpy as np def construct_function_space_for_adjoint(coords, shapeOnRef, mesh, quadratureRule, mode2D='cartesian'): @@ -25,4 +26,9 @@ def construct_function_space_for_adjoint(coords, shapeOnRef, mesh, quadratureRul parentElement=mesh.parentElement, parentElement1d=mesh.parentElement1d, blocks=mesh.blocks, nodeSets=mesh.nodeSets, sideSets=mesh.sideSets) - return FunctionSpace.FunctionSpace(shapes, vols, shapeGrads, mesh, quadratureRule, isAxisymmetric) + if mesh.parentElement.elementType == Interpolants.LAGRANGE_TRIANGLE_ELEMENT: + shapeOnRefHessians = Interpolants.shape2d_lagrange_second_derivatives(mesh.parentElement.degree, mesh.parentElement.coordinates, quadratureRule.xigauss) + shapeHessians = vmap(map_element_shape_hessians, (None, 0, None, None, None))(coords, mesh.conns, mesh.parentElement, shapeOnRef.gradients, shapeOnRefHessians) + return FunctionSpace.FunctionSpace(shapes, vols, np.concatenate((shapeGrads, shapeHessians), axis=-1), mesh, quadratureRule, isAxisymmetric) + else: + return FunctionSpace.FunctionSpace(shapes, vols, shapeGrads, mesh, quadratureRule, isAxisymmetric) From 2aea3bfd252cf443081c09fffb1d016bd0cca9dc Mon Sep 17 00:00:00 2001 From: ralberd Date: Wed, 19 Feb 2025 13:58:50 -0700 Subject: [PATCH 18/19] update Neohookean to work with changes to displacement gradient storage --- optimism/material/Neohookean.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/optimism/material/Neohookean.py b/optimism/material/Neohookean.py index e0a52e99..fc9edac3 100644 --- a/optimism/material/Neohookean.py +++ b/optimism/material/Neohookean.py @@ -24,18 +24,18 @@ def create_material_model_functions(properties): def strain_energy(dispGrad, internalVars, dt): del dt - return energy_density(dispGrad, internalVars, props) + return energy_density(dispGrad[0:3,0:3], internalVars, props) def compute_state_new(dispGrad, internalVars, dt): del dt - return _compute_state_new(dispGrad, internalVars, props) + return _compute_state_new(dispGrad[0:3,0:3], internalVars, props) density = properties.get('density') def compute_material_qoi(dispGrad, internalVars, dt): del internalVars del dt - return _compute_volumetric_jacobian(dispGrad) + return _compute_volumetric_jacobian(dispGrad[0:3,0:3]) return MaterialModel(compute_energy_density = strain_energy, compute_initial_state = make_initial_state, From a350a4854939b9c1315c486681c226aaae9b4bd8 Mon Sep 17 00:00:00 2001 From: ralberd Date: Wed, 19 Feb 2025 14:01:24 -0700 Subject: [PATCH 19/19] adds regularization term to third medium material model --- optimism/material/ThirdMediumNeoHookean.py | 22 +++++++++++++------- optimism/material/test/test_ThirdMedium.py | 24 ++++++++++++++++++++++ 2 files changed, 39 insertions(+), 7 deletions(-) diff --git a/optimism/material/ThirdMediumNeoHookean.py b/optimism/material/ThirdMediumNeoHookean.py index 2d128b34..d16fbc52 100644 --- a/optimism/material/ThirdMediumNeoHookean.py +++ b/optimism/material/ThirdMediumNeoHookean.py @@ -9,19 +9,22 @@ PROPS_MU = 0 PROPS_KAPPA = 1 PROPS_LAMBDA = 2 +PROPS_REGULARIZATION_CONSTANT = 3 def create_material_model_functions(properties): density = properties.get('density') props = _make_properties(properties['bulk modulus'], - properties['shear modulus']) + properties['shear modulus'], + properties['regularization_constant']) def strain_energy(dispGrad, internalVars, dt): del internalVars del dt - # return neo_hookean_energy_density(dispGrad, props) + # return neo_hookean_energy_density(dispGrad[0:3,0:3], props) # return stabilized_neo_hookean_energy_density(dispGrad, props) - return teran_invertible_energy_density(dispGrad, props) + # return teran_invertible_energy_density(dispGrad, props) + return neo_hookean_energy_density(dispGrad[0:3,0:3], props) + _regularization_energy(dispGrad[0:2, 3:], props) def compute_state_new(dispGrad, internalVars, dt): del dispGrad @@ -31,7 +34,7 @@ def compute_state_new(dispGrad, internalVars, dt): def compute_material_qoi(dispGrad, internalVars, dt): del internalVars del dt - return _compute_volumetric_jacobian(dispGrad) + return _compute_volumetric_jacobian(dispGrad[0:3,0:3]) return MaterialModel(compute_energy_density = strain_energy, compute_initial_state = make_initial_state, @@ -42,9 +45,9 @@ def compute_material_qoi(dispGrad, internalVars, dt): def make_initial_state(): return np.array([]) -def _make_properties(kappa, mu): +def _make_properties(kappa, mu, c): lamda = kappa - (2.0/3.0) * mu - return np.array([mu, kappa, lamda]) + return np.array([mu, kappa, lamda, c]) def neo_hookean_energy_density(dispGrad, props): F = dispGrad + np.eye(3) @@ -102,4 +105,9 @@ def _stabilizing_energy_extension(dispGrad, props, J_min): def _compute_volumetric_jacobian(dispGrad): F = dispGrad + np.eye(3) - return np.linalg.det(F) \ No newline at end of file + return np.linalg.det(F) + +def _regularization_energy(dispHessian, props): + def third_order_inner_product(A, B): + return np.dot(A[:,0:2].flatten(), B[:,0:2].flatten()) + 2.0*np.dot(A[:,3], B[:,3]) + return 0.5 * props[PROPS_REGULARIZATION_CONSTANT] * third_order_inner_product(dispHessian, dispHessian) diff --git a/optimism/material/test/test_ThirdMedium.py b/optimism/material/test/test_ThirdMedium.py index ed1bd59b..4166cd23 100644 --- a/optimism/material/test/test_ThirdMedium.py +++ b/optimism/material/test/test_ThirdMedium.py @@ -78,6 +78,30 @@ def test_energy_consistency(self): energy = energy_neo_hookean_adagio(dispGrad, self.kappa, self.mu) self.assertAlmostEqual(energy, energy_gold, 12) + + def test_regularization_term(self): + props = { + 'bulk modulus': 0.0, + 'shear modulus': 0.0, + 'regularization_constant': 1.0 + } + material = ThirdMediumNeoHookean.create_material_model_functions(props) + + key = jax.random.PRNGKey(1) + dispHessian = jax.random.uniform(key, (2, 3)) + dispGradAndHessian = np.zeros((3,6)) + dispGradAndHessian.at[0:2,3:].set(dispHessian) + + internalVariables = material.compute_initial_state() + dt = 0.0 + + energy = material.compute_energy_density(dispGradAndHessian, internalVariables, dt) + + dispHessianFull = np.zeros((2,4)) + dispHessianFull.at[:,0:3].set(dispHessian) + dispHessianFull.at[:,3].set(dispHessian[:,-1]) + exact = 0.5*np.dot(dispHessianFull.flatten(), dispHessianFull.flatten()) + self.assertAlmostEqual(energy.item(), exact, 12) def test_plot_biaxial_response(self): n_steps = 100