diff --git a/optimism/FunctionSpace.py b/optimism/FunctionSpace.py index d22dff3b..9c2752fc 100644 --- a/optimism/FunctionSpace.py +++ b/optimism/FunctionSpace.py @@ -112,16 +112,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, np.concatenate((shapeGrads, shapeHessians), axis=-1), mesh, quadratureRule, isAxisymmetric) + else: + return FunctionSpace(shapes, vols, shapeGrads, 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] diff --git a/optimism/Interpolants.py b/optimism/Interpolants.py index 76f1dee9..d558df74 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. @@ -44,9 +48,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 +63,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 +92,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 +141,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 +218,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, np.array([], dtype=np.int32)) def pascal_triangle_monomials(degree): p = [] @@ -261,6 +331,110 @@ 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 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) @@ -268,6 +442,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/Mechanics.py b/optimism/Mechanics.py index 751ee25c..5b487ab9 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): @@ -252,15 +254,40 @@ 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(): 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)) ###### 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 b20d4745..a20a94d7 100644 --- a/optimism/ReadExodusMesh.py +++ b/optimism/ReadExodusMesh.py @@ -6,8 +6,7 @@ exodusToNativeTri6NodeOrder = np.array([0, 3, 1, 5, 4, 2]) - -def read_exodus_mesh(fileName): +def read_exodus_mesh(fileName, interpolationType = Interpolants.InterpolationType.LOBATTO): with netCDF4.Dataset(fileName) as exData: coords = _read_coordinates(exData) conns, blocks = _read_blocks(exData) @@ -17,15 +16,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 == Interpolants.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/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) 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/material/Neohookean.py b/optimism/material/Neohookean.py index cfb5ec76..fc9edac3 100644 --- a/optimism/material/Neohookean.py +++ b/optimism/material/Neohookean.py @@ -24,17 +24,23 @@ 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[0:3,0:3]) + 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 diff --git a/optimism/material/ThirdMediumNeoHookean.py b/optimism/material/ThirdMediumNeoHookean.py new file mode 100644 index 00000000..d16fbc52 --- /dev/null +++ b/optimism/material/ThirdMediumNeoHookean.py @@ -0,0 +1,113 @@ +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 +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['regularization_constant']) + + def strain_energy(dispGrad, internalVars, dt): + del internalVars + del dt + # 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 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 + del dt + return internalVars + + def compute_material_qoi(dispGrad, internalVars, dt): + del internalVars + del dt + return _compute_volumetric_jacobian(dispGrad[0:3,0:3]) + + 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, c): + lamda = kappa - (2.0/3.0) * mu + return np.array([mu, kappa, lamda, c]) + +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 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) + +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 new file mode 100644 index 00000000..4166cd23 --- /dev/null +++ b/optimism/material/test/test_ThirdMedium.py @@ -0,0 +1,159 @@ +import unittest +from optimism.test.TestFixture import TestFixture +import jax +import jax.numpy as np +from optimism.material import Neohookean +from optimism.material import ThirdMediumNeoHookean + +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 + +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 + 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_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 + 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_teran_invertible(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+legends5,[ + "Rokos Neo Hookean Volumetric term", + "Rokos Neo Hookean Isochoric term", + "Stable Neohookean energy", + "Teran invertible Neo Hookean", + "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() diff --git a/optimism/test/test_FunctionSpace.py b/optimism/test/test_FunctionSpace.py index 7baa8017..27a04125 100644 --- a/optimism/test/test_FunctionSpace.py +++ b/optimism/test/test_FunctionSpace.py @@ -284,5 +284,54 @@ 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=2) + + # 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) + + 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)) + + # 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)) + self.assertTrue(np.allclose(dispHessians, exact)) + + + if __name__ == '__main__': unittest.main() diff --git a/optimism/test/test_Interpolants.py b/optimism/test/test_Interpolants.py index 6354b213..7a36316c 100644 --- a/optimism/test/test_Interpolants.py +++ b/optimism/test/test_Interpolants.py @@ -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.assertCorrectGradInterpolationOn2DElement(element, numRandomEvalPoints, Interpolants.shape2d) - 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.assertArrayNear(expected1, dfInterpolated[:,1], 13) def no_test_plot_high_order_nodes(self): degree = 10 @@ -228,5 +239,166 @@ 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) + 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) + 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.shape == (0,)) + + 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) + + 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) + + + if __name__ == '__main__': unittest.main() diff --git a/optimism/test/test_PatchTest.py b/optimism/test/test_PatchTest.py index 9bd797d7..421db4f1 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) + 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)) - 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) - - self.b = -(2*kappa + 8*mu/3.0) * np.array([alpha, beta]) - - mcxFuncs = Mechanics.create_mechanics_functions(self.fs, "plane strain", materialModel) + def test_dirichlet_patch_test_with_quadratic_elements(self): + self.assertQuadraticPatchTestPasses() - 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__':