Problems in command (updateMaterialStage)

Forum for asking and answering questions related to use of the OpenSeesPy module

Moderators: silvia, selimgunay, Moderators

Post Reply
francojj
Posts: 7
Joined: Sat May 23, 2020 9:22 am

Problems in command (updateMaterialStage)

Post by francojj » Mon Aug 03, 2020 9:15 pm

Dear Michael H. I copied line by line the code of a three-dimensional soil column from the TCL source code to PY, but I have problems with the results during the material update.

The results in the design spectrum without executing the material update command exactly match me from the TCL source code with PY.

Code: Select all

ops.updateMaterialStage('-material', k, '-stage', 0)
Image


While when I activate the material update the results vary in the design spectrum

Code: Select all

ops.updateMaterialStage('-material', k, '-stage', 1)
Image


Code: Select all


# =============================================================================
#         ######### Crear nodos del suelo ##########
# =============================================================================
        # creación de nodos de presión de poros
        ops.model('basic', '-ndm', 3, '-ndf', 4)
                
        with open(path + '/Post-proceso/' + perfil + '/ppNodesInfo.dat', 'w') as f:
            count = 0.0
            yCoord=0.0
            nodos = []
            dryNode = []
            altura_nf = espesor_total - nf
            
            
            for k in range(capas):
                for j in range(0, int(nNodeY[k]), 4):
                    ops.node(j+count+1,0.0, yCoord, 0.0)
                    ops.node(j+count+2, 0.0, yCoord, sElemZ)
                    ops.node(j+count+3, sElemX, yCoord, sElemZ)                    
                    ops.node(j+count+4, sElemX, yCoord, 0.0) 
                               
                    nodos.append(str(j+count+1))
                    nodos.append(str(j+count+2))
                    nodos.append(str(j+count+3))
                    nodos.append(str(j+count+4))
                      
                    #designate node sobre la superficie de agua  
                    if yCoord>=altura_nf:
                        dryNode.append(j+count+1)
                        dryNode.append(j+count+2)
                        dryNode.append(j+count+3)
                        dryNode.append(j+count+4)
 
                    yCoord=(yCoord+sElemY[k])
                   

                count=(count+nNodeY[k])
     
        print("Finished creating all soil nodes...")

# =============================================================================
#         ####### Condiciones de contorno en la base de la columna #########
# =============================================================================      
        ops.fix(1,*[0,1,1,0])
        ops.fix(2,*[0,1,1,0])
        ops.fix(3,*[0,1,1,0])
        ops.fix(4,*[0,1,1,0])
        ops.equalDOF(1,2,1)
        ops.equalDOF(1,3,1)
        ops.equalDOF(1,4,1)
        
        print('Fin de creación de nodos de la base de la columna\n\n') 
        
# =============================================================================
#         ####### Condiciones de contorno en los nudos restantes #########
# =============================================================================            
        
        count=0
        for k in range(5,int(nNodeT+1),4):
            ops.equalDOF(k, k + 1, *[1, 2, 3])
            ops.equalDOF(k, k + 2, *[1, 2, 3])
            ops.equalDOF(k, k + 3, *[1, 2, 3])
            
            
        print('Fin de creación equalDOF para nodos de presión de poros\n\n')
        
        for j in range(len(dryNode)):
            ops.fix(dryNode[j], *[0, 0, 0, 1])

        print("Finished creating all soil boundary conditions...")
        
# =============================================================================
#         ####### crear elemento y material de suelo #########
# =============================================================================
        
        cargas = []
        for j in range(capas):
            pendiente = permutaciones[i][9][j] 
            slope = math.atan(pendiente / 100)

            if tipo_suelo == 'No cohesivo':
                if float(surf) > 0:
                    ops.nDMaterial('PressureDependMultiYield02', j + 1, 3.0, rho, Gr, Br, fric, gmax, refpress, presscoef, ptang, 
                                   cc1, cc3, cd1, cd3, float(surf), 5.0, 3.0, *[1.0,0.0], ev, *[0.9, 0.02, 0.7, 101.0])
                else:
                    ops.nDMaterial('PressureDependMultiYield02', j + 1,3.0, rho, Gr, Br, fric, gmax, refpress, presscoef, ptang, 
                                   cc1, cc3, cd1, cd3,float(surf), *permutaciones[i][29][j], 5.0, 3.0, *[1.0,0.0], ev, *[0.9, 0.02, 0.7, 101.0])
            elif tipo_suelo == 'Cohesivo':
                
                if float(surf) > 0:
                     ops.nDMaterial('PressureIndependMultiYield', j + 1, 3.0, rho, Gr, Br, 
                                coh, gmax, fric, refpress, presscoef, float(surf))
                else:
                     ops.nDMaterial('PressureIndependMultiYield', j + 1, 3.0, rho, Gr, Br, 
                                coh, gmax, fric, refpress, presscoef, float(surf), *permutaciones[i][29][j])
           
        
            cargas.append([0.0, -9.81 * math.cos(slope), -9.81 * math.sin(slope)])
            
            
            
        print('Fin de la creación de material de suelo\n\n')
       
        
#-----------------------------------------------------------------------------------------
#  5. CREATE SOIL ELEMENTS
#-----------------------------------------------------------------------------------------
               
        count=0
        alpha=1.5e-6
        
        with open(path + '/Post-proceso/' + perfil + '/ppElemInfo.dat', 'w') as f:
            # crear elemento de suelo
            for k in range(capas):
                for j in range(int(nElemY[k])):
                    nI = 4 * (j + count + 1) - 3
                    nJ = nI + 1
                    nK = nI + 2
                    nL = nI + 3
                    nM = nI + 4
                    nN = nI + 5
                    nO = nI + 6
                    nP = nI + 7
                    f.write(str(j+count+1) +'\t'+ str(nI) + '\t'+ str(nJ) + '\t'+ str(nK) + '\t'+ str(nL) + '\t'+ 
                            str(nM) + '\t'+ str(nN) + '\t'+ str(nO) + '\t'+ str(nP) + '\n')
                    
                    Bc = permutaciones[i][14][k]
                    
                    ops.element('SSPbrickUP', (j + count + 1), *[nI, nJ, nK, nL, nM, nN, nO, nP], int(k+1), float(Bc), 
                                1.0, 1.0, 1.0, 1.0, float(ev), alpha,cargas[k][0], cargas[k][1], cargas[k][2])
                    
                    
                count=(count+int(nElemY[k]))
        print('Fin de la creación del elemento del suelo\n\n')
        
        #win.ui.progressBar.setValue(25)
        
        
# =============================================================================
#         ######### Amortiguamiento de Lysmer ##########
# =============================================================================       
        
        ops.model('basic', '-ndm', 3, '-ndf', 3)
        
        # definir nodos y coordenadas del amortiguamiento
        dashF = nNodeT + 1
        dashX = nNodeT + 2
        dashZ = nNodeT + 3
        
        ops.node(dashF, 0.0, 0.0, 0.0)
        ops.node(dashX, 0.0, 0.0, 0.0)
        ops.node(dashZ, 0.0, 0.0, 0.0)
        
        # definir restricciones para los nodos de amortiguamiento
        ops.fix(dashF, 1, 1, 1)
        ops.fix(dashX, 0, 1, 1)
        ops.fix(dashZ, 1, 1, 0)
        
        # definir equalDOF para el amortiguamiento en la base del suelo
        ops.equalDOF(1, dashX, 1)
        ops.equalDOF(1, dashZ, 3)
        
        print('Fin de la creación de condiciones de contorno de los nodos de amortiguamiento\n\n')
        
        # definir el material de amortiguamiento
        colArea = sElemX * sElemZ
        dashpotCoeff = vel * den * colArea
        ops.uniaxialMaterial('Viscous', capas + 1, dashpotCoeff, 1)
        
        # definir el elemento 
        ops. element('zeroLength', nElemT + 1, *[dashF, dashX], '-mat', capas + 1, '-dir', *[1])
        ops. element('zeroLength', nElemT + 2, *[dashF, dashZ], '-mat', capas + 1, '-dir', *[3])
        
        print('Fin de la creación del elemento de amortiguamiento\n\n')
        
        
#-----------------------------------------------------------------------------------------
#  9. DEFINE ANALYSIS PARAMETERS
#-----------------------------------------------------------------------------------------    
          
        # amortiguamiento de Rayleigh
        # frecuencia menor
        omega1 = 2 * math.pi * 0.2 
        # frecuencia mayor
        omega2 = 2 * math.pi * 20
        
        
        a0 = 2.0 * (amort / 100) * omega1 * omega2 / (omega1 + omega2)
        a1 = 2.0 * (amort / 100) / (omega1 + omega2)
        print('Coeficientes de amortiguamiento' +'\n'+ 'a0: ' + format(a0, '.6f') + '\n' + 'a1: ' + format(a1, '.6f') + '\n\n')
        
        #win.ui.progressBar.setValue(35)
        
              
# =============================================================================
#         ######## Determinación de análisis estático #########
# =============================================================================
        #---DETERMINE STABLE ANALYSIS TIME STEP USING CFL CONDITION
        # se determina a partir de un análisis transitorio de largo tiempo
        duration = nstep * dt
        
        # tamaño mínimo del elemento y velocidad máxima
        minSize = sElemY[0]
        vsMax = permutaciones[i][11][0]
        for j in range(1,capas):
            if sElemY[j] < minSize:
                minSize = sElemY[j]
            if permutaciones[i][11][j] > vsMax:
                vsMax = permutaciones[i][11][j]
        
        # trial analysis time step        
        kTrial = minSize / (vsMax**0.5)
       
        # tiempo de análisis y pasos de tiempo
        if dt <= kTrial:
            nStep = nstep
            dT = dt
        else:
            nStep = int(math.floor(duration / kTrial) + 1)
            dT = duration / nStep
            
        print('Número de pasos en el análisis: ' + str(nStep) + '\n')
        print('Incremento de tiempo: ' + str(dT) + '\n\n')

        
#----------------------------------------------------------------------------------------
#  7. GRAVITY ANALYSIS
#-----------------------------------------------------------------------------------------
        ops.model('basic', '-ndm', 3, '-ndf', 4)
           
        for k in range(capas):
            ops.updateMaterialStage('-material', int(k+1), '-stage', 0)
        
        # algoritmo de análisis estático
        ops.constraints('Penalty', 1e14, 1e14)
        ops.test('NormDispIncr', 1.0e-5, 30, 1)
        ops.algorithm('Newton')
        ops.numberer('Plain')
        ops.system('SparseGeneral')
        ops.integrator('Newmark', 0.5, 0.25)
        ops.analysis('Transient')
        
        print('Inicio de análisis estático elástico\n\n')
        
        ops.start()
        ops.analyze(20, 5.0e2)
      
        print('Fin de análisis estático elástico\n\n')

      
# =============================================================================
          
        for k in range(capas):
            ops.updateMaterialStage('-material', int(k+1), '-stage', 1)

# =============================================================================
        [size=150][size=85][size=50][/size][/size][/size]
        
        # plastic gravity loading
        print('Inicio de análisis estático plástico\n\n')
        
        ok=ops.analyze(40, 5.0e-2)
            
        if ok != 0:
            error = 'Error de convergencia en análisis estático de modelo' + str(perfil) + '\n\n'
            print(error)
     
            break

        print('Fin de análisis estático plástico\n\n')

#-----------------------------------------------------------------------------------------
#  11. UPDATE ELEMENT PERMEABILITY VALUES FOR POST-GRAVITY ANALYSIS
#-----------------------------------------------------------------------------------------

        ini= 1
        aum = 0
        sum = 0
        for j in range(capas):
            #Layer 3
            ops.setParameter('-val', permutaciones[i][16][j], ['-eleRange', int(ini+aum),int(nElemY[j]+sum)], 'xPerm')
            ops.setParameter('-val', permutaciones[i][17][j], ['-eleRange', int(ini+aum),int(nElemY[j]+sum)], 'yPerm')
            ops.setParameter('-val', permutaciones[i][16][j], ['-eleRange', int(ini+aum),int(nElemY[j]+sum)], 'zPerm')
            
            ini = nElemY[j] + sum
            sum += nElemY[j]
            aum = 1
            
        print( "Finished updating permeabilities for dynamic analysis...")
        

# =============================================================================
#         ########### Grabadores dinámicos ##########
# =============================================================================

        ops.setTime(0.0)
        ops.wipeAnalysis()
        ops.remove('recorders')
        
        # tiempo de la grabadora
        recDT = 10 * dt

        # crear grabadores de análisis dinámico para nodos
        ops.recorder('Node', '-file', path_disp + 'displacementxyz.out', '-time', '-dT', recDT, '-node', *nodos, '-dof', 1,2,3, 'disp')
        ops.recorder('Node', '-file', path_acel + 'accelerationxyz.out', '-time', '-dT', recDT, '-node', *nodos, '-dof', 1,2,3, 'accel')

        
# =============================================================================
#         ######### Determinación de análisis dinámico ##########
# =============================================================================
        
        # objeto de serie temporal para el historial de fuerza
        
        path_vel = path + '/Pre-proceso/' + perfil + '/VelHistoria.txt'
        
 
        ops.timeSeries('Path', 1, '-dt', dt, '-filePath', path_vel, '-factor', dashpotCoeff)
        
   
        ops.pattern('Plain', 10, 1)
        ops.load(1, 1.0, 0.0, 0.0, 0.0)
        
        print('Fin de creación de carga dinámica\n\n')
       

        # algoritmo de análisis dinámico
         ops.constraints('Penalty', 1e14, 1e14)
         ops.test('NormDispIncr', 1.0e-3, 55, 1)
         ops.algorithm('KrylovNewton')
         ops.numberer('Plain')
         ops.system('SparseGeneral')
         ops.integrator('Newmark', 0.5, 0.25)
         ops.analysis('Transient')=========================================
#       ops.rayleigh(a0, a1, 0.0, 0.0)
# =============================================================================
        
        print('Inicio de análisis dinámico\n\n')
        #win.ui.progressBar.setValue(85)
        ok = ops.analyze(nStep, dT)


I also upload the code that I am transforming into TCL format:

https://colab.research.google.com/drive ... sp=sharing


REGARDS
Top
mhscott
Posts: 874
Joined: Tue Jul 06, 2004 3:38 pm
Location: Corvallis, Oregon USA
Contact:
Contact mhscott

Re: Problems in command (updateMaterialStage)

Post by mhscott » Fri Aug 14, 2020 6:18 pm

Lo siento. Voy a hacerlo pronto.
Top
mhscott
Posts: 874
Joined: Tue Jul 06, 2004 3:38 pm
Location: Corvallis, Oregon USA
Contact:
Contact mhscott

Re: Problems in command (updateMaterialStage)

Post by mhscott » Fri Aug 14, 2020 6:27 pm

I think it's related to this: https://wp.me/pbejzW-Si

Inside setParameter for the Pressure...02.cpp file, there's a call to atoi(matTag). In your python input, use str(1) instead of 1 in the updateMaterialStage command.
Top
francojj
Posts: 7
Joined: Sat May 23, 2020 9:22 am

Re: Problems in command (updateMaterialStage)

Post by francojj » Sun Aug 16, 2020 5:34 am

Thanks for your reply Ph.D. Michael,
I tried to emit the arguments using the (str) function in the UpdateMaterial command, but I got a syntax error when running the code.

Code: Select all

WARNING updateMaterialStage: value must be integer
Traceback (most recent call last):
  File "VArena1.py", line 682, in <module>
    analisis_opensees(path, permuta)
  File "VArena1.py", line 462, in analisis_opensees
    ops.updateMaterialStage('-material', int(k+1), '-stage', str(1))
opensees.OpenSeesError: See stderr outputt

Thanks for your great contribution
Top
mhscott
Posts: 874
Joined: Tue Jul 06, 2004 3:38 pm
Location: Corvallis, Oregon USA
Contact:
Contact mhscott

Re: Problems in command (updateMaterialStage)

Post by mhscott » Sun Aug 16, 2020 10:34 am

OK, thanks for checking. The str() thing is not the problem. :(

Please send me input files (one Tcl, one Python) of updateMaterialStage called on a single element.
Top
mhscott
Posts: 874
Joined: Tue Jul 06, 2004 3:38 pm
Location: Corvallis, Oregon USA
Contact:
Contact mhscott

Re: Problems in command (updateMaterialStage)

Post by mhscott » Sun Aug 16, 2020 12:18 pm

I'm not sure what's up, but can you try changing "-stage" to "updateMaterialStage" in your Python command?
Top
francojj
Posts: 7
Joined: Sat May 23, 2020 9:22 am

Re: Problems in command (updateMaterialStage)

Post by francojj » Mon Aug 17, 2020 10:08 pm

Doctor, I checked the code again at the TCL and PY of a single type of soil and there is no variation, I do not know because in the elastic state if the results coincide while the plastic state does not.

https://github.com/fjjimenez1/Effect-site-column-3D.git



Make the mentioned changes to the updateMaterialStage command and I keep getting errors:

Code: Select all

WARNING updateMaterialStage: Only accept parameter '-stage' for now
Traceback (most recent call last):

  File "C:\Users\Franco\Desktop\V_TFT\Geometry\Field_PY.py", line 552, in <module>
    analisis_opensees(path, permuta)
  File "C:\Users\Franco\Desktop\V_TFT\Geometry\Field_PY.py", line 361, in analisis_opensees
    ops.updateMaterialStage('-material', 1, 'updateMaterialStage',1)

OpenSeesError: See stderr output
Regards
Top
Post Reply

Return to “OpenSeesPy”