Verified Commit fa92caae authored by Lars Bilke's avatar Lars Bilke

Format python: Tests/.

parent 3a6977e0
from paraview.simple import *
from paraview import coprocessing
#--------------------------------------------------------------
# --------------------------------------------------------------
# Code generated from cpstate.py to create the CoProcessor.
# ParaView 5.2.0 64 bits
# ----------------------- CoProcessor definition -----------------------
def CreateCoProcessor():
def _CreatePipeline(coprocessor, datadescription):
class Pipeline:
# state file generated using paraview version 5.2.0
# ----------------------------------------------------------------
# setup the data processing pipelines
# ----------------------------------------------------------------
#### disable automatic camera reset on 'Show'
paraview.simple._DisableFirstRenderCameraReset()
# create a new 'XML Unstructured Grid Reader'
# create a producer from a simulation input
cube_1e0_pcs_0_ts_0_t_0000000vtu = coprocessor.CreateProducer(datadescription, 'input')
# create a new 'Contour'
contour1 = Contour(Input=cube_1e0_pcs_0_ts_0_t_0000000vtu)
contour1.ContourBy = ['POINTS', 'D1_left_front_N1_right']
contour1.ComputeScalars = 1
contour1.Isosurfaces = [1.0, 1.0750344444444444, 1.1500688888888888, 1.2251033333333334, 1.3001377777777778, 1.3751722222222222, 1.4502066666666669, 1.5252411111111113, 1.6002755555555557, 1.67531]
contour1.PointMergeMethod = 'Uniform Binning'
# create a new 'Parallel PolyData Writer'
parallelPolyDataWriter1 = servermanager.writers.XMLPPolyDataWriter(Input=contour1)
# register the writer with coprocessor
# and provide it with information such as the filename to use,
# how frequently to write the data, etc.
coprocessor.RegisterWriter(parallelPolyDataWriter1, filename='cube_1e1_%t.pvtp', freq=1)
# ----------------------------------------------------------------
# finally, restore active source
SetActiveSource(contour1)
# ----------------------------------------------------------------
return Pipeline()
class CoProcessor(coprocessing.CoProcessor):
def CreatePipeline(self, datadescription):
self.Pipeline = _CreatePipeline(self, datadescription)
coprocessor = CoProcessor()
# these are the frequencies at which the coprocessor updates.
freqs = {'input': [1]}
coprocessor.SetUpdateFrequencies(freqs)
return coprocessor
#--------------------------------------------------------------
def _CreatePipeline(coprocessor, datadescription):
class Pipeline:
# state file generated using paraview version 5.2.0
# ----------------------------------------------------------------
# setup the data processing pipelines
# ----------------------------------------------------------------
#### disable automatic camera reset on 'Show'
paraview.simple._DisableFirstRenderCameraReset()
# create a new 'XML Unstructured Grid Reader'
# create a producer from a simulation input
cube_1e0_pcs_0_ts_0_t_0000000vtu = coprocessor.CreateProducer(
datadescription, "input"
)
# create a new 'Contour'
contour1 = Contour(Input=cube_1e0_pcs_0_ts_0_t_0000000vtu)
contour1.ContourBy = ["POINTS", "D1_left_front_N1_right"]
contour1.ComputeScalars = 1
contour1.Isosurfaces = [
1.0,
1.0750344444444444,
1.1500688888888888,
1.2251033333333334,
1.3001377777777778,
1.3751722222222222,
1.4502066666666669,
1.5252411111111113,
1.6002755555555557,
1.67531,
]
contour1.PointMergeMethod = "Uniform Binning"
# create a new 'Parallel PolyData Writer'
parallelPolyDataWriter1 = servermanager.writers.XMLPPolyDataWriter(
Input=contour1
)
# register the writer with coprocessor
# and provide it with information such as the filename to use,
# how frequently to write the data, etc.
coprocessor.RegisterWriter(
parallelPolyDataWriter1, filename="cube_1e1_%t.pvtp", freq=1
)
# ----------------------------------------------------------------
# finally, restore active source
SetActiveSource(contour1)
# ----------------------------------------------------------------
return Pipeline()
class CoProcessor(coprocessing.CoProcessor):
def CreatePipeline(self, datadescription):
self.Pipeline = _CreatePipeline(self, datadescription)
coprocessor = CoProcessor()
# these are the frequencies at which the coprocessor updates.
freqs = {"input": [1]}
coprocessor.SetUpdateFrequencies(freqs)
return coprocessor
# --------------------------------------------------------------
# Global variables that will hold the pipeline for each timestep
# Creating the CoProcessor object, doesn't actually create the ParaView pipeline.
# It will be automatically setup when coprocessor.UpdateProducers() is called the
# first time.
coprocessor = CreateCoProcessor()
#--------------------------------------------------------------
# --------------------------------------------------------------
# Enable Live-Visualizaton with ParaView
coprocessor.EnableLiveVisualization(False, 1)
# ---------------------- Data Selection method ----------------------
def RequestDataDescription(datadescription):
"Callback to populate the request for current timestep"
global coprocessor
......@@ -86,8 +106,10 @@ def RequestDataDescription(datadescription):
# pipeline.
coprocessor.LoadRequestedData(datadescription)
# ------------------------ Processing method ------------------------
def DoCoProcessing(datadescription):
"Callback to do co-processing for current timestep"
global coprocessor
......
import OpenGeoSys
from math import pi, sin, cos, sinh, cosh
a = 2.0*pi/3.0
a = 2.0 * pi / 3.0
# analytical solution used to set the Dirichlet BCs
def solution(x, y):
return sin(a*x) * sinh(a*y)
return sin(a * x) * sinh(a * y)
# gradient of the analytical solution used to set the Neumann BCs
def grad_solution(x, y):
return a * cos(a*x) * sinh(a*y), \
a * sin(a*x) * cosh(a*y)
return a * cos(a * x) * sinh(a * y), a * sin(a * x) * cosh(a * y)
# Dirichlet BCs
class BCTop(OpenGeoSys.BoundaryCondition):
......@@ -21,6 +21,7 @@ class BCTop(OpenGeoSys.BoundaryCondition):
value = solution(x, y)
return (True, value)
class BCLeft(OpenGeoSys.BoundaryCondition):
def getDirichletBCValue(self, t, coords, node_id, primary_vars):
x, y, z = coords
......@@ -28,6 +29,7 @@ class BCLeft(OpenGeoSys.BoundaryCondition):
value = solution(x, y)
return (True, value)
class BCBottom(OpenGeoSys.BoundaryCondition):
def getDirichletBCValue(self, t, coords, node_id, primary_vars):
x, y, z = coords
......@@ -35,13 +37,14 @@ class BCBottom(OpenGeoSys.BoundaryCondition):
value = solution(x, y)
return (True, value)
# Neumann BC
class BCRight(OpenGeoSys.BoundaryCondition):
def getFlux(self, t, coords, primary_vars):
x, y, z = coords
assert x == 1.0 and z == 0.0
value = grad_solution(x, y)[0]
Jac = [ 0.0 ] # value does not depend on primary variable
Jac = [0.0] # value does not depend on primary variable
return (True, value, Jac)
......
import OpenGeoSys
from math import pi, sin
a = 2.0*pi
b = 2.0*pi
a = 2.0 * pi
b = 2.0 * pi
def solution(x, y):
return sin(a*x-pi/2.0) * sin(b*y-pi/2.0)
return sin(a * x - pi / 2.0) * sin(b * y - pi / 2.0)
# - laplace(solution) = source term
def laplace_solution(x, y):
return a*a * sin(a*x-pi/2.0) * sin(b*y-pi/2.0) + b*b * sin(a*x-pi/2.0) * sin(b*y-pi/2.0)
return a * a * sin(a * x - pi / 2.0) * sin(b * y - pi / 2.0) + b * b * sin(
a * x - pi / 2.0
) * sin(b * y - pi / 2.0)
# source term for the benchmark
class SinXSinYSourceTerm(OpenGeoSys.SourceTerm):
def getFlux(self, t, coords, primary_vars):
x, y, z = coords
value = laplace_solution(x,y)
Jac = [ 0.0 ]
value = laplace_solution(x, y)
Jac = [0.0]
return (value, Jac)
# instantiate source term object referenced in OpenGeoSys' prj file
sinx_siny_source_term = SinXSinYSourceTerm()
#--------------------------------------------------------------
# --------------------------------------------------------------
# Global timestep output options
timeStepToStartOutputAt = 0
......@@ -12,15 +12,15 @@ rescale_lookuptable = False
requestSpecificArrays = False
# a root directory under which all Catalyst output goes
rootDirectory = ''
rootDirectory = ""
# makes a cinema D index table
make_cinema_table = False
#--------------------------------------------------------------
# --------------------------------------------------------------
# Code generated from cpstate.py to create the CoProcessor.
# paraview version 5.8.0
#--------------------------------------------------------------
# --------------------------------------------------------------
from paraview.simple import *
from paraview import coprocessing
......@@ -46,13 +46,13 @@ def CreateCoProcessor():
paraview.simple._DisableFirstRenderCameraReset()
# create a new 'Superquadric'
ogs_output = coprocessor.CreateProducer(datadescription, 'input')
ogs_output = coprocessor.CreateProducer(datadescription, "input")
# create a new 'Clip'
clip1 = Clip(Input=ogs_output)
clip1.ClipType = 'Plane'
clip1.HyperTreeGridClipper = 'Plane'
clip1.Scalars = ['POINTS', 'D1_left_bottom_N1_right']
clip1.ClipType = "Plane"
clip1.HyperTreeGridClipper = "Plane"
clip1.Scalars = ["POINTS", "D1_left_bottom_N1_right"]
clip1.Value = 1.3376572416618493
# init the 'Plane' selected for 'ClipType'
......@@ -63,13 +63,20 @@ def CreateCoProcessor():
# create a new 'Contour'
contour1 = Contour(Input=ogs_output)
contour1.ContourBy = ['POINTS', 'D1_left_bottom_N1_right']
contour1.ContourBy = ["POINTS", "D1_left_bottom_N1_right"]
contour1.Isosurfaces = [
1.0, 1.075034942591522, 1.1500698851830442, 1.2251048277745662,
1.3001397703660882, 1.3751747129576102, 1.4502096555491324,
1.5252445981406544, 1.6002795407321764, 1.6753144833236986
1.0,
1.075034942591522,
1.1500698851830442,
1.2251048277745662,
1.3001397703660882,
1.3751747129576102,
1.4502096555491324,
1.5252445981406544,
1.6002795407321764,
1.6753144833236986,
]
contour1.PointMergeMethod = 'Uniform Binning'
contour1.PointMergeMethod = "Uniform Binning"
# ----------------------------------------------------------------
# finally, restore active source
......@@ -77,31 +84,35 @@ def CreateCoProcessor():
# ----------------------------------------------------------------
# Now any catalyst writers
xMLPUnstructuredGridWriter1 = servermanager.writers.XMLPUnstructuredGridWriter(
Input=clip1)
xMLPUnstructuredGridWriter1 = (
servermanager.writers.XMLPUnstructuredGridWriter(Input=clip1)
)
coprocessor.RegisterWriter(
xMLPUnstructuredGridWriter1,
filename='square_1e1_neumann_clip1_%t.pvtu',
filename="square_1e1_neumann_clip1_%t.pvtu",
freq=1,
paddingamount=0,
DataMode='Appended',
HeaderType='UInt64',
DataMode="Appended",
HeaderType="UInt64",
EncodeAppendedData=False,
CompressorType='None',
CompressionLevel='6')
CompressorType="None",
CompressionLevel="6",
)
xMLPPolyDataWriter1 = servermanager.writers.XMLPPolyDataWriter(
Input=contour1)
Input=contour1
)
coprocessor.RegisterWriter(
xMLPPolyDataWriter1,
filename='square_1e1_neumann_contour1_%t.pvtp',
filename="square_1e1_neumann_contour1_%t.pvtp",
freq=1,
paddingamount=0,
DataMode='Appended',
HeaderType='UInt64',
DataMode="Appended",
HeaderType="UInt64",
EncodeAppendedData=False,
CompressorType='None',
CompressionLevel='6')
CompressorType="None",
CompressionLevel="6",
)
return Pipeline()
......@@ -113,8 +124,7 @@ def CreateCoProcessor():
# these are the frequencies at which the coprocessor updates.
freqs = {}
coprocessor.SetUpdateFrequencies(freqs)
coprocessor.SetInitialOutputOptions(timeStepToStartOutputAt,
forceOutputAtFirstCall)
coprocessor.SetInitialOutputOptions(timeStepToStartOutputAt, forceOutputAtFirstCall)
if rootDirectory:
coprocessor.SetRootDirectory(rootDirectory)
......@@ -125,14 +135,14 @@ def CreateCoProcessor():
return coprocessor
#--------------------------------------------------------------
# --------------------------------------------------------------
# Global variable that will hold the pipeline for each timestep
# Creating the CoProcessor object, doesn't actually create the ParaView pipeline.
# It will be automatically setup when coprocessor.UpdateProducers() is called the
# first time.
coprocessor = CreateCoProcessor()
#--------------------------------------------------------------
# --------------------------------------------------------------
# Enable Live-Visualizaton with ParaView and the update frequency
coprocessor.EnableLiveVisualization(False, 1)
......@@ -163,10 +173,12 @@ def DoCoProcessing(datadescription):
coprocessor.WriteData(datadescription)
# Write image capture (Last arg: rescale lookup table), if appropriate.
coprocessor.WriteImages(datadescription,
rescale_lookuptable=rescale_lookuptable,
image_quality=0,
padding_amount=imageFileNamePadding)
coprocessor.WriteImages(
datadescription,
rescale_lookuptable=rescale_lookuptable,
image_quality=0,
padding_amount=imageFileNamePadding,
)
# Live Visualization, if enabled.
coprocessor.DoLiveVisualization(datadescription, "localhost", 22222)
......@@ -8,10 +8,11 @@ def createStressArray(points):
sigma.SetNumberOfTuples(points.GetNumberOfPoints())
for i in range(points.GetNumberOfPoints()):
y = points.GetPoint(i)[1]
sigma.SetTuple4(i, 2200 * 9.81 * 0.8 * (y - 100),
-0.5 * 1e6 + 2200 * 9.81 * (y - 100), 0, 0)
#sigma.SetTuple4(i, 4*i, 4*i + 1, 4*i + 2, 4*i + 3) # for debugging
sigma.SetName('nonequilibrium_stress')
sigma.SetTuple4(
i, 2200 * 9.81 * 0.8 * (y - 100), -0.5 * 1e6 + 2200 * 9.81 * (y - 100), 0, 0
)
# sigma.SetTuple4(i, 4*i, 4*i + 1, 4*i + 2, 4*i + 3) # for debugging
sigma.SetName("nonequilibrium_stress")
return sigma
......@@ -32,14 +33,15 @@ def writeDataToFile(mesh, filename):
w.SetInputData(mesh)
w.Write()
def addStressArraysToMesh(mesh):
mesh.GetPointData().AddArray(createStressArrayForNodes(mesh))
mesh.GetCellData().AddArray(createStressArrayForCells(mesh))
r = vtkXMLUnstructuredGridReader()
r.SetFileName('soil_column.vtu')
r.SetFileName("soil_column.vtu")
r.Update()
m = r.GetOutput()
addStressArraysToMesh(m)
writeDataToFile(m, 'soil_column_nonequilibrium_sigma.vtu')
writeDataToFile(m, "soil_column_nonequilibrium_sigma.vtu")
......@@ -12,10 +12,10 @@ def distribute_points_evenly(c2):
assert np.sqrt(c2.shape[0]).is_integer()
CELLS_PER_DIRECTION = int(np.sqrt(c2.shape[0])) - 1
r2 = np.sqrt(c2[:,0]**2 + c2[:,1]**2)
alpha2 = np.arctan2(c2[:,1], c2[:,0])
r2 = np.sqrt(c2[:, 0] ** 2 + c2[:, 1] ** 2)
alpha2 = np.arctan2(c2[:, 1], c2[:, 0])
bins = [ [] for _ in range(CELLS_PER_DIRECTION+1) ]
bins = [[] for _ in range(CELLS_PER_DIRECTION + 1)]
nbin = np.round(r2 * CELLS_PER_DIRECTION).astype(int)
for node, b in enumerate(nbin):
bins[b].append(node)
......@@ -53,10 +53,10 @@ a = np.empty(coords.shape[0])
for i, (x, y, z) in enumerate(coords):
if x > y:
R = np.sqrt(1 + (y/x)**2)
R = np.sqrt(1 + (y / x) ** 2)
a[i] = 1.0 / R
elif x < y:
R = np.sqrt(1 + (x/y)**2)
R = np.sqrt(1 + (x / y) ** 2)
a[i] = 1.0 / R
else:
a[i] = np.sqrt(0.5)
......@@ -81,28 +81,32 @@ writer.SetInputData(grid)
writer.Write()
# extract boundary of unit circle
R_squared = new_coords[:,0]**2 + new_coords[:,1]**2
phi = np.arctan2(new_coords[:,1], new_coords[:,0])
R_squared = new_coords[:, 0] ** 2 + new_coords[:, 1] ** 2
phi = np.arctan2(new_coords[:, 1], new_coords[:, 0])
idcs = np.where(abs(R_squared - 1) < 1e-8)[0]
idcs = idcs[np.argsort(phi[idcs])] # sorted with ascending polar angle
with open(out_geom, "w") as fh:
fh.write("""<?xml version="1.0" encoding="ISO-8859-1"?>
fh.write(
"""<?xml version="1.0" encoding="ISO-8859-1"?>
<OpenGeoSysGLI>
<name>geom</name>
<points>
<point id="0" x="0" y="0" z="0" name="center"/>
<point id="1" x="0" y="1" z="0" name="top"/>
<point id="2" x="1" y="0" z="0"/>
""")
"""
)
for i, (x, y, z) in enumerate(new_coords[idcs]):
fh.write(' <point id="{}" x="{}" y="{}" z="{}" />\n'.format(
i+3, x, y, z))
fh.write(
' <point id="{}" x="{}" y="{}" z="{}" />\n'.format(i + 3, x, y, z)
)
fh.write("""
fh.write(
"""
</points>
<polylines>
......@@ -114,11 +118,14 @@ with open(out_geom, "w") as fh:
<pnt>0</pnt>
<pnt>2</pnt>
</polyline>
<polyline id="2" name="outer">\n""")
<polyline id="2" name="outer">\n"""
)
for i in range(len(idcs)):
fh.write(" <pnt>{}</pnt>\n".format(i+3))
fh.write(" <pnt>{}</pnt>\n".format(i + 3))
fh.write("""</polyline>
fh.write(
"""</polyline>
</polylines>
</OpenGeoSysGLI>""")
</OpenGeoSysGLI>"""
)
......@@ -11,26 +11,26 @@ class HertzContactBC(OpenGeoSys.BoundaryCondition):
def __init__(self):
super(HertzContactBC, self).__init__()
self._first_node = None # ID of the first node of this BC's geometry
self._first_node = None # ID of the first node of this BC's geometry
self._t_old = START_TIME - 1.0 # time of previous invocation of this BC
self._boundary_x_coords = [] # the x coordinates of all boundary nodes
self._boundary_x_coords = [] # the x coordinates of all boundary nodes
def _init_timestep(self):
"""Initializes the internal state at the beginning of a new timestep."""
self._a_range = [ 0.0, SPHERE_RADIUS ] # range of possible contact radii
self._a_est = SPHERE_RADIUS # estimated contact radius
self._max_x_with_y_excess = -1.0 # maximum x value where a node is above the contact line
self._a_range = [0.0, SPHERE_RADIUS] # range of possible contact radii
self._a_est = SPHERE_RADIUS # estimated contact radius
self._tendency_is_up = True # whether the contact area is supposed to grow in the current iteration
self._a_curr = 0.0 # the radius of the contact area in this iteration (continuously updated)
self._a_prev = 0.0 # the radius of the contact area that was prescripted in the previous iteration
self._max_x_with_y_excess = (
-1.0
) # maximum x value where a node is above the contact line
self._iteration = 0 # the nonlinear solver iteration number
self._tendency_is_up = True # whether the contact area is supposed to grow in the current iteration
self._a_curr = 0.0 # the radius of the contact area in this iteration (continuously updated)
self._a_prev = 0.0 # the radius of the contact area that was prescripted in the previous iteration
self._iteration = 0 # the nonlinear solver iteration number
def _init_iteration(self):
"""Initializes the internal state at the beginning of a new nonlinear solver iteration."""
......@@ -66,9 +66,11 @@ class HertzContactBC(OpenGeoSys.BoundaryCondition):
self._iteration += 1
print("BC: a_est={:.4f}, a_prev={:.4f} ({:.4f}, {:.4f})".format(
self._a_est, self._a_prev, self._a_range[0], self._a_range[1]))
print(
"BC: a_est={:.4f}, a_prev={:.4f} ({:.4f}, {:.4f})".format(
self._a_est, self._a_prev, self._a_range[0], self._a_range[1]
)
)
def getDirichletBCValue(self, t, coords, node_id, primary_vars):
if self._t_old < t:
......@@ -93,10 +95,15 @@ class HertzContactBC(OpenGeoSys.BoundaryCondition):
try:
# check that we are at the outer boundary
assert abs(x**2 + y**2 + z**2 - SPHERE_RADIUS**2) < 1e-15
assert abs(x ** 2 + y ** 2 + z ** 2 - SPHERE_RADIUS ** 2) < 1e-15
except:
print("assert abs(x**2 + y**2 + z**2 - 1.0) < 1e-15",
x, y, z, x**2 + y**2 + z**2 - SPHERE_RADIUS**2)
print(
"assert abs(x**2 + y**2 + z**2 - 1.0) < 1e-15",
x,
y,
z,
x ** 2 + y ** 2 + z ** 2 - SPHERE_RADIUS ** 2,
)
raise
y_deformed = y + uy
......@@ -128,20 +135,17 @@ class HertzContactBC(OpenGeoSys.BoundaryCondition):
res = (True, y_top - y)
elif self._boundary_x_coords_are_initialized(t):
idx = self._boundary_x_coords.index(x)
if idx != 0 and self._boundary_x_coords[idx-1] == self._a_prev:
if idx != 0 and self._boundary_x_coords[idx - 1] == self._a_prev: