04-08-2020

geometry:

    film_length : 0.5 #x length of film

    film_height : 0.5 # y height of film

    film_width : 1.0 #z width of film

    film_length : 1.5 #x length of film

    film_height : 1.5 # y height of film

    film_width : 1.5 #z width of film

    mesh: 0.5               #Mesh size in mm

    El_type: "C3D8" #Element type

UMAT: "3D_SMA_UM.f"

    init_temp: 0.0

    steps: 2  #Number of steps

    type: 1 #1- Force, 2 - Displacement, 3 - Temperature

    value: 150000

    value: 150000000 #X direction

# +####################################################################

# ############### Base Model for tensile test with cracks ###############

# #######################################################################

# importing modules

import regionToolset

import job

from sketch import *

import mesh as msh

from mesh import *

from load import *

from interaction import *

from step import *

from assembly import *

from section import *

from material import *

from part import *

from sys import path

import sys

import random

import numpy as np

from numpy import random

from job import *

# Change the path for Site packages

path.append("/home/apadhy/.virtualenvs/abaqus_py/lib/python2.7/site-packages")

import yaml

import os

import math

import copy

################# create parts ###############

def create_part(geometry, model):

    """[summary]

    Arguments:

        geometry {[dict]} -- [subdictionary geometry in Yaml file]

        model {[modelobj]} -- [Model object]

    """

    length = float(geometry["film_length"])   	# Getting length of the composite from dictionary

    height = float(geometry["film_height"])  # Getting height of film from dictionary

    width = float(geometry["film_width"])      # Getting extrusion width of all parts from dictionary

    # Defining the film

    model.ConstrainedSketch(name="film_sketch", sheetSize=200.0) ## Create a sketch named "Mid" with sheet size for display area on screen

    model.sketches["film_sketch"].rectangle(

        point1=(0.0, 0.0), point2=(length, height)) # Create a rectangle with corner points

    model.Part(dimensionality=THREE_D, name="film", type=DEFORMABLE_BODY)

    model.parts["film"].BaseSolidExtrude(

        sketch=model.sketches["film_sketch"], depth=width) # Define part as 3D deformable

	# Defining the assembly

    # Create instances of each part

    model.rootAssembly.DatumCsysByDefault(CARTESIAN)

    model.rootAssembly.Instance(

        dependent=ON, name="film-1", part=model.parts["film"]) #Needs to be a dependent object

# ############### create Materials ###############

def material(materials, model):

    """[Creates materials]

    Arguments:

        materials {[dict]} -- [Dictionary of materials in Yaml]

        model {[Model]} -- [Abaqus model object]

    """

    #Define engineering constants and variables from YAML dictionary

    Density=float(materials["Density"])

    Depvar= int(materials["Depvar"])

    Specific_h= float(materials["Specific_heat"])

    Conductivity= float(materials["Conductivity"])

    User= materials["User_material"]

    print (User)

    model.Material(name="SMA_UM")

    model.materials['SMA_UM'].Density(table=((Density, ), ))

    model.materials['SMA_UM'].Depvar(n=Depvar)

    model.materials['SMA_UM'].Conductivity(table=((Conductivity, ), ))

    model.materials['SMA_UM'].SpecificHeat(table=((Specific_h, ), ))

    #model.materials['SMA_UM'].UserMaterial(mechanicalConstants=(float(User["E"]),float(User["NU"])))

    model.materials['SMA_UM'].UserMaterial(mechanicalConstants=(float(User["IPHASE"]),float(User["MODEL"]),

    float(User["TOL"]), float(User["xi0"]),int(User["NELMTP"]),float(User["EA"]),float(User["EM"]), float(User["nu"]),float(User["alphaA"]),

    float(User["alphaM"]),float(User["Mos"]),float(User["Mof"]),float(User["Aos"]),float(User["Aof"]),float(User["H"]),float(User["rSoMH"]),

    float(User["rSoAH"]),float(User["epstr11"]),float(User["epstr22"]),float(User["epstr33"]),float(User["2epstr23"]),float(User["2epstr13"]),

    float(User["2epstr12"]), float(User["FRULE"])))

## ############### assign sections ###############

def section(geometry,model):

    """[Creates sections for materials and assigns them to element sets chosen at random]

    Arguments:

        geometry {[dict]} -- [Dictionary for geometry from Yaml file]

        materials {[dict]} -- [Dictionary for materials from Yaml file]

        model {[model]} -- [Abaqus model object]

    """

    l=geometry["film_length"]

    h=geometry["film_height"]

    w=geometry["film_width"]

    mesh=geometry["mesh"]

    l_div=int(round(l/mesh,2))

    h_div=int(round(h/mesh,2))

    w_div=int(round(w/mesh,2))

    el_total=l_div*w_div*h_div

    string=""

    for el in range(1,el_total+1):

        string+= "{},".format(el)

        if el%14==0:

            string+="\n"

    model.keywordBlock.synchVersions(storeNodesAndElements=False)

    for n, line in enumerate(model.keywordBlock.sieBlocks):

        if "*End Part" in line:

            kwds1="*Elset, elset=EALL\n"+string+"\n"

            kwds2="** Section: Section-SMA\n"

            kwds3="*Solid Section, elset=EALL, material=SMA_UM,CONTROLS=HG-1\n"

            model.keywordBlock.insert(n-1, kwds1+kwds2+kwds3)

            kwds4="*SECTION CONTROLS, NAME=HG-1, HOURGLASS=ENHANCED"

            model.keywordBlock.insert(n+1, kwds4)

def section_old(geometry, model):

    # define variables used

    ###############################################

    # Create section outer and middle from materials and width

    ###Section upper

    MM = model.parts["film"]

    c = MM.cells

    model.HomogeneousSolidSection(

            material="SMA_UM", name="SMA_Section", thickness=None)

    region = regionToolset.Region(cells=c)

    MM.SectionAssignment(region=region, sectionName="SMA_Section", offset=0.0,

                        offsetType=MIDDLE_SURFACE, offsetField="",

                        thicknessAssignment=FROM_SECTION)

########################  Create SEEDS for MESH Generation  ########################################

def edge_seeding(geometry):

    """[Creates seeds on edges in cells from mesh size]

    """

    l=geometry["film_length"]

    h=geometry["film_height"]

    w=geometry["film_width"]

    mesh=geometry["mesh"]

    l_div=int(round(l/mesh,2))

    h_div=int(round(h/mesh,2))

    w_div=int(round(w/mesh,2))

    n_total= int((l_div+1)*(w_div+1)*(h_div+1))

    el_total=l_div*w_div*h_div

    print ("Total nodes",n_total)

    print ("Total elements",el_total)

    N_coord=[]

    El_coord=dict()

    Node=dict()

    N_index=1

    El=0

    for y in range(0,h_div):

        for x in range(0,l_div):

            for z in range(0,w_div):

                coord=[[mesh*x,mesh*y,mesh*z],[mesh*(x+1),mesh*y,mesh*z],[mesh*(x+1),mesh*(y+1),mesh*z],

                             [mesh*x,mesh*(y+1),mesh*z],[mesh*x,mesh*y,mesh*(z+1)],[mesh*(x+1),mesh*y,mesh*(z+1)],

                             [mesh*(x+1),mesh*(y+1),mesh*(z+1)],[mesh*x,mesh*(y+1),mesh*(z+1)]]

                El+=1

                El_index=[]

                for n in range(len(coord)):

                    if coord[n] not in N_coord:

                        N_coord.append(coord[n])

                        Node[N_index]=coord[n]

                        El_index.append(N_index)

                        N_index+=1

                        #print (coord[n])

                        #print("###")

                    else:

                        for index,coordi in Node.items():

                            if coordi == coord[n]:

                                El_index.append(index)

                El_coord[El]=El_index

    # print(El_coord)

    # print(Node)

    Node_out= ""

    El_out=""

    for key,value in sorted(Node.items()) :

        string= "\t{},\t{}\n".format(key,str(value)[1:-1])

        #print string

        Node_out+=string

    for key,value in sorted(El_coord.items()) :

        string= "\t{},\t{}\n".format(key,str(value)[1:-1])

        #print string

        El_out+=string

    return Node_out,El_out

###########################################

#####################  Generate MESH   ###########################################

def Mesh(part,geometry,model):

    """[Mesh each part]

    Arguments:

        part {[partobj]} -- [Part object]

        elemtype {[str]} -- [Element type]

    """

    Node_out,El_out = edge_seeding(geometry)

    El_type=geometry["El_type"]

    model.keywordBlock.synchVersions(storeNodesAndElements=False)

    for n, line in enumerate(model.keywordBlock.sieBlocks):

        if "*Part" in line:

            l=line.split("\n")[0]

            kwds= l+"\n"+"*Node\n"

            model.keywordBlock.replace(n, kwds)

            kwds1="{}".format(Node_out)+"*Element, type="+El_type+"\n"+"{}\n".format(El_out)

            model.keywordBlock.insert(n, kwds1)

yaml_path = os.path.join(os.getcwd(), "3D_UM.yaml")

# load the yaml file and immediatly close it after loading

with open(yaml_path, "r") as stream:

    try:

        input=yaml.safe_load(stream)

    except yaml.YAMLError as exc:

        sys.__stderr__.writelines(exc)

inp_dir= os.path.dirname(yaml_path)

os.chdir(inp_dir)

model = mdb.models["Model-1"]                # Specify model number

geometry= input["geometry"]

materials= input["materials"]

create_part(geometry, model)

part= model.parts["film"]

material(materials, model)

section(geometry,model)

Mesh(part,geometry,model)

job_name="JOB"

model.rootAssembly.regenerate()

mdb.Job(name=job_name, model="Model-1", description="", type=ANALYSIS,

    atTime=None, waitMinutes=0, waitHours=0, queue=None, memory=90,

    memoryUnits=PERCENTAGE , getMemoryFromAnalysis=True, #

    explicitPrecision=SINGLE, nodalOutputPrecision=SINGLE, echoPrint=OFF,

    modelPrint=OFF, contactPrint=OFF, historyPrint=OFF, userSubroutine="",

    scratch="", resultsFormat=ODB, multiprocessingMode=DEFAULT, numCpus=2, numGPUs=0,numDomains=2)

try:

    mdb.jobs[job_name].writeInput(consistencyChecking=OFF)

except Exception as ex:

    sys.__stderr__.writelines("ERROR JOB IS: {} \n".format(job_name))

    sys.__stderr__.writelines(str(ex))

# +####################################################################

# ############### Base Model for tensile test with cracks ###############

# #######################################################################

# importing modules

import regionToolset

import job

from sketch import *

import mesh as msh

from mesh import *

from load import *

from interaction import *

from step import *

from assembly import *

from section import *

from material import *

from part import *

from sys import path

import sys

import random

import numpy as np

from numpy import random

from job import *

# Change the path for Site packages

path.append("/home/apadhy/.virtualenvs/abaqus_py/lib/python2.7/site-packages")

import yaml

import os

import math

import copy

################ create step ###############

def step(model):

    """Define load step with solution settings

    Arguments:

        model {[model]} -- Abaqus model object

    """

    # Default parameters for FEM simulation set here

    model.StaticStep(name='Loading', previous='Initial',

        timePeriod=100.0, maxNumInc=1000, initialInc=1.0, minInc=0.001,

        maxInc=1.0)

    model.StaticStep(name='Unload', previous='Loading',

        timePeriod=100.0, maxNumInc=1000, initialInc=1.0, minInc=0.001,

        maxInc=1.0)

    # model.StaticStep(adaptiveDampingRatio=0.05,

    #               continueDampingFactors=False, initialInc=1, maxInc=100, maxNumInc=100,

    #               minInc=0.001, name="Unload", nlgeom=OFF, previous="Loading",

    #               stabilizationMagnitude=0.0002, stabilizationMethod=DAMPING_FACTOR)

    #model.ImplicitDynamicsStep(name='Loading', previous='Initial',

    #    maxNumInc=10000, initialInc=0.1, minInc=1e-06, nlgeom=OFF)

    #model.ImplicitDynamicsStep(name='Unload', previous='Loading',

    #    maxNumInc=10000, initialInc=0.1, minInc=1e-06)

    ## AMPLITUDE

    #model.SmoothStepAmplitude(name='Amp-1', timeSpan=STEP, data=((

    #    0.0, 0.0), (1.0, 1.0)))

    #model.SmoothStepAmplitude(name='Amp-2', timeSpan=STEP, data=((

    #    0.0, 1.0), (1.0, 0.0)))

## ############### create boundary conditions ###############

def BC(loading, geometry, model):

    """Sets boundary conditions for the body

    """

    a = model.rootAssembly

    p= model.parts["FILM"]

    mesh=geometry["mesh"]/2

    model.rootAssembly.regenerate()

    value = loading["value"]    # Prescribed displacement

    l = geometry["film_length"]

    w = geometry["film_width"]

    h= geometry["film_height"]

    init_temp= loading["init_temp"]

    #Left nodes fixed in X

    #####################################

    n0= a.instances["FILM-1"].nodes.getByBoundingBox(-mesh,-mesh,-mesh,mesh,h+mesh,w+mesh)

    region = regionToolset.Region(nodes=n0)

    a.Set(nodes=n0, name="BC-fix-X")

    model.DisplacementBC(name="BC-fix-X",

                         region=region, u1=0.0, u2=UNSET,createStepName='Initial',

                         u3=UNSET, ur1=UNSET, ur2=UNSET, ur3=UNSET, amplitude='', fixed=OFF,

                         distributionType=UNIFORM, fieldName="", localCsys=None)         # Displacements set

    ####################################################

    # Bottom and top edge nodes fixed in Y and Z

    n0= a.instances["FILM-1"].nodes.getByBoundingCylinder((0,0,0),(l-mesh,0,0),mesh)

    n1= a.instances["FILM-1"].nodes.getByBoundingCylinder((0,h,w),(l-mesh,h,w),mesh)

    region = regionToolset.Region(nodes=n1+n0)

    a.Set(nodes=n1+n0, name="BC-fix-YZ")

    model.DisplacementBC(name="BC-YZ",

                         region=region, u1=UNSET, u2=0.0,createStepName='Initial',

                         u3=0.0, ur1=0.0, ur2=UNSET, ur3=UNSET, amplitude='', fixed=OFF,

                         distributionType=UNIFORM, fieldName="", localCsys=None)         # Displacements set

    # Displacements set

    ###############################################

    # Initial temperature assign to all nodes

    n2= a.instances["FILM-1"].nodes.getByBoundingBox(-mesh,-mesh,-mesh,l+mesh,h+mesh,w+mesh)

    region = regionToolset.Region(nodes=n2)

    a.Set(nodes=n2, name="BC-Temp")

    model.Temperature(name='Predefined Field-1',

        createStepName='Initial', region=region, distributionType=UNIFORM,

        crossSectionDistribution=CONSTANT_THROUGH_THICKNESS, magnitudes=(325.0,

        ))

    # right faces prescribed displacement in X

    #side1Faces = p.faces.getByBoundingBox(-0.01+l,-0.01,-0.01,l+0.01,h+0.01,w+0.01)

    s = a.instances["FILM-1"].nodes.getByBoundingBox(l-mesh,-mesh,-mesh,l+mesh,h+mesh,w+mesh)

    if loading["type"]==1:

        region = regionToolset.Region(nodes=s)

        a.Set(nodes=s, name="BC-F")

        model.ConcentratedForce(name='Load-1', createStepName='Loading',

            region=region,cf1=value,amplitude=UNSET)

        model.loads['Load-1'].setValuesInStep(stepName='Unload', cf1=0.0)

            # mdb.models['Model-1'].boundaryConditions['Bearing1'].setValuesInStep(

        #     stepName='Unload', amplitude='Amp-2')

        #model.boundaryConditions['Disp-1'].setValuesInStep(stepName='Unload', amplitude='')

    elif loading["type"]==2:

        a.Surface(side1Faces=side1Faces, name="BC-disp")

        Surf=a.surfaces["BC-disp"]

        model.DisplacementBC(name="Disp-1", createStepName='Loading',

                            region=Surf, u1=value, u2=UNSET,

                            u3=UNSET, ur1=UNSET, ur2=UNSET, ur3=UNSET, amplitude=UNSET, fixed=OFF,

                            distributionType=UNIFORM, fieldName="", localCsys=None)

        model.DisplacementBC(name="Disp-1", createStepName='Unload',

                            region=region, u1=0.0, u2=UNSET,

                            u3=UNSET, ur1=UNSET, ur2=UNSET, ur3=UNSET, amplitude=UNSET, fixed=OFF,

                            distributionType=UNIFORM, fieldName="", localCsys=None)

        #model.boundaryConditions['Disp-1'].setValuesInStep(stepName='Unload', amplitude='')

# ############### set Output requests ###############

def outputreq(model):

    """Defines output variables to be shown in Visualization after solving

    Arguments:

        model {object} -- model object

    """

    model.FieldOutputRequest(name="F-Output-1" ,createStepName='Loading',

        variables=('S', 'U','RF','SDV'))

    model.FieldOutputRequest(name="F-Output-2",createStepName='Unload',

        variables=('S', 'U','RF','SDV'))

    model.HistoryOutputRequest(name='H-Output-1',

        createStepName='Loading', variables=('IRA1', 'IRA2', 'IRA3', 'IRAR1',

        'IRAR2', 'IRAR3'))

    model.HistoryOutputRequest(name='H-Output-2',

        createStepName='Unload', variables=('IRA1', 'IRA2', 'IRA3', 'IRAR1',

        'IRAR2', 'IRAR3'))

    #'S', 'PE', 'PEEQ', 'PEMAG', 'LE', 'U', 'V', 'A', 'RF','SDV'), numIntervals=50, timeMarks=OFF)

    #####

    # SDV1 Direction of transformation indicator

    # SDV2 Martensitic volume fraction

    # SDV3 Flag for transformation 0 - No 1 - Yes

    # SDV5 Component of transformation strain in 11 direction

    # SDV6 Component of transformation strain in 22 direction

    # SDV7 Component of transformation strain in 33 direction

    # SDV8 2*Component of transformation strain in 23 direction

    # SDV9 2*Component of transformation strain in 13 direction

    # SDV10 2*Component of transformation strain in 12 direction

    # SDV14 Value of effective transformation strain

    # SDV15 Current temperature

    # SDV20 Component of the transformation direction tensor A in 11 direction

    # SDV21 Component of the transformation direction tensor A in 22 direction

    # SDV22 Component of the transformation direction tensor A in 33 direction

    # SDV23 Component of the transformation direction tensor A in 23 direction

    # SDV24 Component of the transformation direction tensor A in 13 direction

    # SDV25 Component of the transformation direction tensor A in 12 direction

    # SDV28 Value of the martensitic volume fraction at point of reversal of transformation

#####################################################

####################################################

yaml_path = os.path.join(os.getcwd(), "3D_UM.yaml")

# load the yaml file and immediatly close it after loading

with open(yaml_path, "r") as stream:

    try:

        input=yaml.safe_load(stream)

    except yaml.YAMLError as exc:

        sys.__stderr__.writelines(exc)

inp_dir= os.path.dirname(yaml_path)

INP=os.path.join(os.getcwd(),"JOB.inp")

CAE="SMA"              #os.path.join(inp_dir,(os.path.splitext(os.path.basename(yaml_path))[0]))

os.chdir(inp_dir)

Mdb()

geometry= input["geometry"]

materials= input["materials"]

loading= input["loading"]

### Function calls  #####

mdb.ModelFromInputFile(name='JOB',

    inputFileName=INP)

model = mdb.models["JOB"]   # Specify model number

step(model)

outputreq(model)

BC(loading, geometry, model)

mdb.saveAs(pathName=CAE)

job_name="SMA"

model.rootAssembly.regenerate()

sub=  os.path.join(os.getcwd(),input["UMAT"])

sys.__stderr__.writelines("JOB IS: {} \n".format(job_name))

mdb.Job(name=job_name, model="JOB", description="", type=ANALYSIS,

    atTime=None, waitMinutes=0, waitHours=0, queue=None, memory=90,

    memoryUnits=PERCENTAGE , getMemoryFromAnalysis=True, #

    explicitPrecision=SINGLE, nodalOutputPrecision=SINGLE, echoPrint=OFF,

    modelPrint=OFF, contactPrint=OFF, historyPrint=OFF, userSubroutine=sub,

    scratch="", resultsFormat=ODB, multiprocessingMode=DEFAULT, numCpus=2, numGPUs=0,numDomains=2)

# try:

#     mdb.jobs[job_name].writeInput(consistencyChecking=ON)

# except Exception as ex:

#     sys.__stderr__.writelines("ERROR JOB IS: {} \n".format(job_name))

#     sys.__stderr__.writelines(str(ex))

geometry:

    film_length : 0.5 #x length of film

    film_height : 0.5 # y height of film

    film_width : 1.0 #z width of film

    mesh: 0.25               #Mesh size in mm

UMAT: "1_Elasticity.f"

materials:

    Density: 6450

    Depvar: 100

    Conductivity: 22.0

    Specific_heat: 329.0

    User_material:

        E: 70e9

        NU: 0.33

        # IPHASE: 1   # 1: Austenite, 2: Martensite

        # MODEL: 3  # Constitutive model 1: Tanaka's exponential, 2: Bo and Lagoudas, 3: Liand and Roger's cosine model

        # TOL: 1.0e-8  # Convergence Tolerance

        # xi0: 0.0  #Initial volume of the martensitic volume fraction

        # NELMTP: 16  #Number of integration points in all SMA finite elements

        # EA: 70e9   # Young's modulus for austenite, Pa

        # EM: 30e9 # YOung's modulus for martensite, Pa

        # nu: 0.33 # Poissons ratio

        # alphaA: 22e-6    #Thermal expansion coefficient of austenite K-1

        # alphaM:  22e-6  #Thermal expansion coefficient of martensite K-1

        # Mos: 291.0   #Martensite start temperature

        # Mof: 271.0   #Martensite finish temperature

        # Aos: 295.0  #Austenite start temperature

        # Aof: 315.0   #Austenite finish temperature

        # H: 0.05    #Maximum current transformation strains

        # rSoMH: -0.35e6     #Martensite stress influence coefficient, Pa/K

        # rSoAH:  -0.35e6    #Austenite stress influence coefficient, Pa/K

        # epstr11: 0.0  #Initial value of transformation strain tensor 11 component

        # epstr22: 0.0   #Initial value of transformation strain tensor 22 component

        # epstr33: 0.0   #Initial value of transformation strain tensor 33 component

        # 2epstr23: 0.0    #2X Initial value of transformation strain tensor 23 component

        # 2epstr13: 0.0   #2X Initial value of transformation strain tensor 13 component

        # 2epstr12: 0.0   #2X Initial value of transformation strain tensor 12 component

        # FRULE: 1  #Flag for the form of the transformation tensor 1- Direction dependent 2- Direction independent

#Loading condition in mm,N

loading:

    init_temp: 0.0

    steps: 2  #Number of steps

    type: 1 #1- Force, 2 - Displacement, 3 - Temperature

    value: 150000