diff --git a/examples/multi_material_tutorial/uniaxial_two_material.py b/examples/multi_material_tutorial/uniaxial_two_material.py new file mode 100644 index 00000000..0304be64 --- /dev/null +++ b/examples/multi_material_tutorial/uniaxial_two_material.py @@ -0,0 +1,167 @@ +import numpy as onp # as is "old numpy" + +from optimism.JaxConfig import * +from optimism import EquationSolver +from optimism import FunctionSpace +from optimism import Interpolants +from optimism.material import J2Plastic +from optimism.material import Neohookean +from optimism import Mechanics +from optimism import Mesh +from optimism import Objective +from optimism import QuadratureRule +from optimism import SparseMatrixAssembler +from optimism import VTKWriter + + +# set up the mesh +w = 0.1 +L = 1.0 +nodesInX = 5 # must be odd +nodesInY = 15 +xRange = [0.0, w] +yRange = [0.0, L] +mesh = Mesh.construct_structured_mesh(nodesInX, nodesInY, xRange, yRange) + +def triangle_centroid(vertices): + return np.average(vertices, axis=0) + +blockIds = -1*onp.ones(Mesh.num_elements(mesh), dtype=onp.int64) +for e, t in enumerate(mesh.conns): + vertices = mesh.coords[t,:] + if triangle_centroid(vertices)[0] < w/2: + blockIds[e] = 0 + else: + blockIds[e] = 1 + +blocks = {'soft': np.flatnonzero(np.array(blockIds) == 0), + 'hard': np.flatnonzero(np.array(blockIds) == 1)} + +mesh = Mesh.mesh_with_blocks(mesh, blocks) + +nodeTol = 1e-8 +nodeSets = {'left': np.flatnonzero(mesh.coords[:,0] < xRange[0] + nodeTol), + 'right': np.flatnonzero(mesh.coords[:,0] > xRange[1] - nodeTol), + 'bottom': np.flatnonzero(mesh.coords[:,1] < yRange[0] + nodeTol), + 'top': np.flatnonzero(mesh.coords[:,1] > yRange[1] - nodeTol)} + +mesh = Mesh.mesh_with_nodesets(mesh, nodeSets) + +# set the essential boundary conditions and create the a DofManager to +# handle the indexing between unknowns and degrees of freedom. +EBCs = [Mesh.EssentialBC(nodeSet='left', field=0), + Mesh.EssentialBC(nodeSet='bottom', field=1), + Mesh.EssentialBC(nodeSet='top', field = 1)] + +fieldShape = mesh.coords.shape +dofManager = Mesh.DofManager(mesh, fieldShape, EBCs) + +# write blocks and bcs to paraview output to check things are correct +writer = VTKWriter.VTKWriter(mesh, baseFileName='check_problem_setup') + +writer.add_cell_field(name='block_id', cellData=blockIds, + fieldType=VTKWriter.VTKFieldType.SCALARS, + dataType=VTKWriter.VTKDataType.INT) + +bcs = np.array(dofManager.isBc, dtype=np.int64) +writer.add_nodal_field(name='bcs', nodalData=bcs, fieldType=VTKWriter.VTKFieldType.VECTORS, + dataType=VTKWriter.VTKDataType.INT) + +writer.write() + +# create the function space +order = 2*mesh.masterElement.degree +quadRule = QuadratureRule.create_quadrature_rule_on_triangle(degree=2*(order - 1)) +fs = FunctionSpace.construct_function_space(mesh, quadRule) + +# create the material models +ESoft = 5.0 +nuSoft = 0.25 +props = {'elastic modulus': ESoft, + 'poisson ratio': nuSoft} +softMaterialModel = Neohookean.create_material_model_functions(props) + +E = 10.0 +nu = 0.25 +Y0 = 0.01*E +H = E/100 +props = {'elastic modulus': E, + 'poisson ratio': nu, + 'yield strength': Y0, + 'hardening model': 'linear', + 'hardening modulus': H} +hardMaterialModel = J2Plastic.create_material_model_functions(props) + +materialModels = {'soft': softMaterialModel, 'hard': hardMaterialModel} + +# mechanics functions +mechanicsFunctions = Mechanics.create_multi_block_mechanics_functions(fs, mode2D="plane strain", materialModels=materialModels) + +# helper function to fill in nodal values of essential boundary conditions +def get_ubcs(p): + appliedDisp = p[0] + EbcIndex = (mesh.nodeSets['top'], 1) + V = np.zeros(fieldShape).at[EbcIndex].set(appliedDisp) + return dofManager.get_bc_values(V) + +# helper function to go from unknowns to full DoF array +def create_field(Uu, p): + return dofManager.create_field(Uu, get_ubcs(p)) + +# write the energy to minimize +def energy_function(Uu, p): + U = create_field(Uu, p) + internalVariables = p.state_data + return mechanicsFunctions.compute_strain_energy(U, internalVariables) + +# Tell objective how to assemble preconditioner matrix +def assemble_sparse_preconditioner_matrix(Uu, p): + U = create_field(Uu, p) + internalVariables = p.state_data + elementStiffnesses = mechanicsFunctions.compute_element_stiffnesses(U, internalVariables) + return SparseMatrixAssembler.assemble_sparse_stiffness_matrix( + elementStiffnesses, mesh.conns, dofManager) + +# solver settings +solverSettings = EquationSolver.get_settings(max_cumulative_cg_iters=100, + max_trust_iters=1000, + use_preconditioned_inner_product_for_cg=True) + +precondStrategy = Objective.PrecondStrategy(assemble_sparse_preconditioner_matrix) + +# initialize unknown displacements to zero +Uu = dofManager.get_unknown_values(np.zeros(fieldShape)) + +# set initial values of parameters +appliedDisp = 0.0 +state = mechanicsFunctions.compute_initial_state() +p = Objective.Params(appliedDisp, state) + +# Construct an objective object for the equation solver to work on +objective = Objective.ScaledObjective(energy_function, Uu, p, precondStrategy=precondStrategy) + +# increment the applied displacement +appliedDisp = L*0.003 +p = Objective.param_index_update(p, 0, appliedDisp) + +# Find unknown displacements by minimizing the objective +Uu = EquationSolver.nonlinear_equation_solve(objective, Uu, p, solverSettings) + +# update the state variables in the new equilibrium configuration +U = create_field(Uu, p) +state = mechanicsFunctions.compute_updated_internal_variables(U, p.state_data) +p = Objective.param_index_update(p, 1, state) + +# write solution data to VTK file for post-processing +writer = VTKWriter.VTKWriter(mesh, baseFileName='uniaxial_two_material') + +U = create_field(Uu, p) +writer.add_nodal_field(name='displacement', nodalData=U, fieldType=VTKWriter.VTKFieldType.VECTORS) + +_, stresses = mechanicsFunctions.compute_output_energy_densities_and_stresses( + U, state) +cellStresses = FunctionSpace.project_quadrature_field_to_element_field(fs, stresses) +writer.add_cell_field(name='stress', cellData=cellStresses, + fieldType=VTKWriter.VTKFieldType.TENSORS) + +writer.write()