yade-users team mailing list archive

[Robert Caulk <question626277@xxxxxxxxxxxxxxxxxxxxx>]

Question #626277 on Yade changed:
https://answers.launchpad.net/yade/+question/626277

Description changed to:
Hello,

I am trying to visualize jcfpm bond break information in paraview using
the cracks vtkrecorder. Yade generates a cracks file filled with the
correct information as well as cracks.vtu files that appear to be empty.
My vtkrecorder is active with the argument recorders=('spheres', 'intr',
'stress', 'jcfpm', 'cracks'), and the recordCracks argument is set to
true within Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM().  I can
visualize the spheres and interaction files in paraview, but the
cracks.vtu files are empty. Am I missing something?

Distribution: Ubuntu xenial 16.04 LTS
Yade versions tried: yadedaily and master branch of yade/trunk 

Here is an MWE, assuming there is an 'out' folder in the working
directory for the vtk files. Any help is appreciated.

#!/usr/bin/python
# -*- coding: utf-8 -*-
from __future__ import division

from yade import plot,pack,timing
import time, sys, os, copy
import numpy as np
import matplotlib.pyplot as plt


young=65e9
poisson=0.3

sigmaT=6e6
frictionAngle=atan(0.8)
epsCrackOnset=1e-4
relDuctility=30
finalFricDegree = 19
density=5000
doModes=2
output='./out/010517'+str(doModes)
intRadius=1.329
dtSafety=.8
damping=0.4
strainRateTension=.05
strainRateCompression=.5
setSpeeds=True
# 1=tension, 2=compression (ANDed; 3=both)

iterper=100
specimenLength=.077
sphereRadius=0.0015
R = specimenLength / 4
numSpheres=10000
# isotropic confinement (should be negative)
isoPrestress=0
cohesion=50e6

jCFmat = O.materials.append(JCFpmMat(young=young, cohesion=cohesion,
density=density, frictionAngle=radians(finalFricDegree),
tensileStrength=sigmaT, poisson=poisson, label='JCFmat',
jointNormalStiffness=2.5e6,jointShearStiffness=1e6,jointCohesion=1e6))

sps=SpherePack()
sp = pack.randomDensePack(pack.inCylinder((0,0,-.5*specimenLength),(0,0,.5*specimenLength), R), radius=sphereRadius, spheresInCell=numSpheres, rRelFuzz=0.25, returnSpherePack=True)
sp.toSimulation(material=jCFmat)

bb=uniaxialTestFeatures()
negIds,posIds,axis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']
O.dt=dtSafety*PWaveTimeStep()
print 'Timestep',O.dt

mm,mx=[pt[axis] for pt in aabbExtrema()]
coord_25,coord_50,coord_75=mm+.25*(mx-mm),mm+.5*(mx-mm),mm+.75*(mx-mm)
area_25,area_50,area_75=approxSectionArea(coord_25,axis),approxSectionArea(coord_50,axis),approxSectionArea(coord_75,axis)

O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius,label='is2aabb'),],verletDist=.05*sphereRadius),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intRadius,label='ss2sc')],
        [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1, label='jcf')],
        [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(Key="Test", smoothJoint=True,label='interactionLaw', recordCracks=True ),Law2_ScGeom_FrictPhys_CundallStrack()],
    ),
    NewtonIntegrator(damping=damping,label='damper'),
    #CpmStateUpdater(realPeriod=.5),
    UniaxialStrainer(strainRate=strainRateTension,axis=axis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=False,blockRotations=False,setSpeeds=setSpeeds,label='strainer'),
    PyRunner(virtPeriod=1e-6/strainRateTension,realPeriod=1,command='addPlotData()',label='plotDataCollector',initRun=True),
    PyRunner(realPeriod=4,command='stopIfDamaged()',label='damageChecker'),
    VTKRecorder(iterPeriod=iterper,initRun=True,fileName=(output+'-'),recorders=['spheres','intr','stress', 'bstresses','jcfpm', 'cracks'])
]

plot.plots={'eps':('sigma',)}

O.saveTmp('initial');

O.timingEnabled=False

global mode
mode='tension' if doModes & 1 else 'compression'

def initTest():
    global mode
    print "init"
    if O.iter>0:
        O.wait();
        O.loadTmp('initial')
        print "Reversing plot data"; plot.reverseData()
    else: plot.plot(subPlots=False)
    strainer.strainRate=abs(strainRateTension) if mode=='tension' else -abs(strainRateCompression)
    try:
        from yade import qt
        renderer=qt.Renderer()
        renderer.dispScale=(1000,1000,1000) if mode=='tension' else (100,100,100)
    except ImportError: pass
    print "init done, will now run."
    O.step(); # to create initial contacts
    # now reset the interaction radius and go ahead
    ss2sc.interactionDetectionFactor=1.
    is2aabb.aabbEnlargeFactor=1.

    O.run()

def stopIfDamaged():
    global mode
    if O.iter<2 or not plot.data.has_key('sigma'): return # do nothing at the very beginning
    sigma,eps=plot.data['sigma'],plot.data['eps']
    extremum=max(sigma) if (strainer.strainRate>0) else min(sigma)
    minMaxRatio=0.5 if mode=='tension' else 0.5
    if extremum==0: return
    # uncomment to get graph for the very first time stopIfDamaged() is called
    #eudoxos.estimatePoissonYoung(principalAxis=axis,stress=strainer.avgStress,plot=True,cutoff=0.3)
    import sys;	sys.stdout.flush()
    if abs(sigma[-1]/extremum)<minMaxRatio or abs(strainer.strain)>(5e-3 if isoPrestress==0 else 5e-2):
        if mode=='tension' and doModes & 2: # only if compression is enabled
            mode='compression'
            O.save('/tmp/uniax-tension.yade.gz')
            print "Saved /tmp/uniax-tension.yade.gz (for use with interaction-histogram.py and uniax-post.py)"
            print "Damaged, switching to compression... "; O.pause()
            # important! initTest must be launched in a separate thread;
            # otherwise O.load would wait for the iteration to finish,
            # but it would wait for initTest to return and deadlock would result
            import thread; thread.start_new_thread(initTest,())
            return
        else:
            print "Damaged, stopping."
            ft,fc=max(sigma),min(sigma)
            #pltfcft[z] = abs(fc/ft)
            #print 'Strengths fc=%g, ft=%g, |fc/ft|=%g'%(fc,ft,abs(fc/ft))
            #title=O.tags['description'] if 'description' in O.tags.keys() # else O.tags['params']
            #print 'gnuplot',plot.saveGnuplot(O.tags['id'],title=title)
            print 'Bye.'
            O.pause()
            #sys.exit(0) # results in some threading exception

def addPlotData():
    yade.plot.addData({'t':O.time,'i':O.iter,'eps':strainer.strain,'sigma':strainer.avgStress+isoPrestress,
                       'sigma.25':forcesOnCoordPlane(coord_25,axis)[axis]/area_25+isoPrestress,
                       'sigma.50':forcesOnCoordPlane(coord_50,axis)[axis]/area_50+isoPrestress,
                       'sigma.75':forcesOnCoordPlane(coord_75,axis)[axis]/area_75+isoPrestress,
                       })
plot.plot(subPlots=False)
#O.run()
initTest()

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.