from abaqus import * from abaqusConstants import * import mesh import os import uuid suffix = uuid.uuid4().hex[:8] model_name = 'Bilayer_Bending_' + suffix job_name = 'Bilayer_Job_' + suffix desktop = '' if os.name == 'nt': try: import winreg key = winreg.OpenKey( winreg.HKEY_CURRENT_USER, r'Software\Microsoft\Windows\CurrentVersion\Explorer\User Shell Folders' ) desktop = os.path.expandvars( winreg.QueryValueEx(key, 'Desktop')[0] ) winreg.CloseKey(key) except Exception: desktop = '' if not desktop: desktop = os.path.join(os.path.expanduser('~'), 'Desktop') if not os.path.isdir(desktop): os.makedirs(desktop) umat_path = os.path.join(desktop, 'Bilayer_UMAT_' + suffix + '.for') umat_source = """ SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD, 1 RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN, 2 TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME, 3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS, 4 DROT,PNEWDT,CELENT,DFGRD0,DFGRD1,NOEL, 5 NPT,LAYER,KSPT,JSTEP,KINC) INCLUDE 'ABA_PARAM.INC' CHARACTER*80 CMNAME DIMENSION STRESS(NTENS),STATEV(NSTATV), 1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS), 2 DRPLDE(NTENS),STRAN(NTENS),DSTRAN(NTENS), 3 TIME(2),PREDEF(1),DPRED(1),PROPS(NPROPS), 4 COORDS(3),DROT(3,3),DFGRD0(3,3), 5 DFGRD1(3,3),JSTEP(4) DOUBLE PRECISION EMOD,ENU,GMOD,ALAM,DSIG DOUBLE PRECISION SOLD(6) INTEGER I,J EMOD=PROPS(1) ENU=PROPS(2) GMOD=EMOD/(2.D0*(1.D0+ENU)) ALAM=EMOD*ENU/((1.D0+ENU)*(1.D0-2.D0*ENU)) DO I=1,NTENS SOLD(I)=STRESS(I) DDSDDT(I)=0.D0 DRPLDE(I)=0.D0 DO J=1,NTENS DDSDDE(I,J)=0.D0 END DO END DO DO I=1,NDI DO J=1,NDI DDSDDE(I,J)=ALAM END DO DDSDDE(I,I)=ALAM+2.D0*GMOD END DO DO I=NDI+1,NTENS DDSDDE(I,I)=GMOD END DO DO I=1,NTENS DSIG=0.D0 DO J=1,NTENS DSIG=DSIG+DDSDDE(I,J)*DSTRAN(J) END DO STRESS(I)=SOLD(I)+DSIG SSE=SSE+0.5D0*(SOLD(I)+STRESS(I))*DSTRAN(I) END DO SPD=0.D0 SCD=0.D0 RPL=0.D0 DRPLDT=0.D0 RETURN END """ with open(umat_path, 'w') as handle: handle.write(umat_source) model = mdb.Model(name=model_name) length = 80.0 width = 10.0 layer_thickness = 2.0 total_thickness = 2.0 * layer_thickness tol = 1.0e-5 sketch = model.ConstrainedSketch(name='BeamProfile', sheetSize=200.0) sketch.rectangle(point1=(0.0, 0.0), point2=(length, width)) part = model.Part( name='BilayerBeam', dimensionality=THREE_D, type=DEFORMABLE_BODY ) part.BaseSolidExtrude(sketch=sketch, depth=total_thickness) del model.sketches['BeamProfile'] plane = part.DatumPlaneByPrincipalPlane( principalPlane=XYPLANE, offset=layer_thickness ) part.PartitionCellByDatumPlane( datumPlane=part.datums[plane.id], cells=part.cells[:] ) lower_cells = part.cells.getByBoundingBox( xMin=-tol, xMax=length + tol, yMin=-tol, yMax=width + tol, zMin=-tol, zMax=layer_thickness + tol ) upper_cells = part.cells.getByBoundingBox( xMin=-tol, xMax=length + tol, yMin=-tol, yMax=width + tol, zMin=layer_thickness - tol, zMax=total_thickness + tol ) if len(lower_cells) != 1 or len(upper_cells) != 1: raise RuntimeError('Layer selection failed.') part.Set(name='LAYER_BUILTIN', cells=lower_cells) part.Set(name='LAYER_UMAT', cells=upper_cells) material = model.Material(name='BuiltinElastic') material.Elastic(table=((70000.0, 0.30),)) material = model.Material(name='UserElastic') material.UserMaterial(mechanicalConstants=(35000.0, 0.30)) for section_name, material_name, set_name in ( ('SectionBuiltin', 'BuiltinElastic', 'LAYER_BUILTIN'), ('SectionUMAT', 'UserElastic', 'LAYER_UMAT') ): model.HomogeneousSolidSection( name=section_name, material=material_name, thickness=None ) part.SectionAssignment( region=part.sets[set_name], sectionName=section_name ) part.setMeshControls( regions=part.cells[:], elemShape=HEX, technique=STRUCTURED ) part.setElementType( regions=(part.cells[:],), elemTypes=(mesh.ElemType(elemCode=C3D8, elemLibrary=STANDARD),) ) part.seedPart(size=1.0, deviationFactor=0.1, minSizeFactor=0.1) part.generateMesh() if len(part.elements) == 0: raise RuntimeError('Mesh generation failed.') assembly = model.rootAssembly assembly.DatumCsysByDefault(CARTESIAN) instance = assembly.Instance( name='Beam-1', part=part, dependent=ON ) assembly.regenerate() fixed_faces = instance.faces.getByBoundingBox( xMin=-tol, xMax=tol, yMin=-tol, yMax=width + tol, zMin=-tol, zMax=total_thickness + tol ) tip_faces = instance.faces.getByBoundingBox( xMin=length - tol, xMax=length + tol, yMin=-tol, yMax=width + tol, zMin=-tol, zMax=total_thickness + tol ) if len(fixed_faces) != 2 or len(tip_faces) != 2: raise RuntimeError('End face selection failed.') assembly.Set(name='FIXED_END', faces=fixed_faces) assembly.Set(name='LOADED_END', faces=tip_faces) model.StaticStep( name='Bending', previous='Initial', nlgeom=OFF, timePeriod=1.0, initialInc=0.05, maxInc=0.05, minInc=1.0e-8, maxNumInc=100 ) model.EncastreBC( name='Clamp', createStepName='Initial', region=assembly.sets['FIXED_END'] ) model.DisplacementBC( name='TipBending', createStepName='Bending', region=assembly.sets['LOADED_END'], u1=UNSET, u2=UNSET, u3=-1.0, ur1=UNSET, ur2=UNSET, ur3=UNSET ) for name in list(model.fieldOutputRequests.keys()): del model.fieldOutputRequests[name] model.FieldOutputRequest( name='Fields', createStepName='Bending', variables=('S', 'E', 'U', 'RF'), frequency=1 ) model.HistoryOutputRequest( name='TipHistory', createStepName='Bending', region=assembly.sets['LOADED_END'], variables=('U3', 'RF3'), frequency=1 ) job = mdb.Job( name=job_name, model=model_name, type=ANALYSIS, userSubroutine=umat_path, numCpus=4, numDomains=4, memory=80, memoryUnits=PERCENTAGE ) job.writeInput(consistencyChecking=ON) inp_path = os.path.abspath(job_name + '.inp') if not os.path.isfile(inp_path): raise RuntimeError('Input file was not written.') with open(inp_path, 'r') as handle: inp_text = handle.read().upper() if '*ELEMENT, TYPE=C3D8' not in inp_text: raise RuntimeError('Expected C3D8 elements were not found.') if '*USER MATERIAL' not in inp_text: raise RuntimeError('User material definition was not found.') if 'Viewport: 1' in session.viewports: viewport = session.viewports['Viewport: 1'] viewport.setValues(displayedObject=assembly) viewport.assemblyDisplay.setValues(mesh=ON) viewport.view.setValues(session.views['Iso']) viewport.view.fitView() print('Model: ' + model_name) print('Job: ' + job_name) print('UMAT: ' + umat_path) print('INP: ' + inp_path) print('Elements: ' + str(len(part.elements))) print('Input written. Job has not been submitted.') print('SOUTHWEITION_JOB_PREPARED') import os as _sw_os try: print('SOUTHWEITION_CWD ' + _sw_os.getcwd()) def _sw_job_model(_job): _model = getattr(_job, 'model', '') try: return str(_model.name) except Exception: return str(_model or '') _sw_models = list(mdb.models.keys()) _sw_latest_model = _sw_models[-1] if _sw_models else '' _sw_name = '' for _cand in reversed(list(mdb.jobs.keys())): if _sw_latest_model and _sw_job_model(mdb.jobs[_cand]) == _sw_latest_model: _sw_name = str(_cand) break if not _sw_name: _sw_names = list(mdb.jobs.keys()) if _sw_names: _sw_name = str(_sw_names[-1]) if _sw_name: try: mdb.jobs[_sw_name].setValues(numCpus=4, numDomains=4) except Exception: pass mdb.jobs[_sw_name].writeInput(consistencyChecking=OFF) print('SOUTHWEITION_JOB ' + str(_sw_name)) print('SOUTHWEITION_MODEL ' + str(_sw_latest_model)) print('SOUTHWEITION_INP ' + _sw_os.path.abspath(str(_sw_name) + '.inp')) print('SOUTHWEITION_JOB_PREPARED') except Exception: print('SOUTHWEITION_JOB_PREPARED')