Example to test post-treatment and visu (for Caroline)

Examples/testVisu.py

+#
+# Example of field creation/discretisation and output
+# to test post-treatment (vtk, hdf5 ...)
+#
+# Usage :
+# python testVisu.py
+# or
+# mpirun -np NB python testVisu.py
+# NB = number of mpi processes
+
+
+import parmepy as pp
+import math
+from parmepy.fields.continuous import Field
+from parmepy.domain.obstacle.sphere import Sphere
+from parmepy.operator.analytic import Analytic
+from parmepy.domain.obstacle.plates import Plates
+pi = math.pi
+from parmepy.operator.penalization import Penalization
+from parmepy.operator.monitors.printer import Printer
+from parmepy.problem.simulation import Simulation
+cos = math.cos
+sin = math.sin
+nb = 33
+Lx = Ly = Lz = 2
+
+
+## ======== 3D ======== ##
+
+# Domain (box-shaped)
+dom = pp.Box(dimension=3, length=[Lx, Ly, Lz], origin=[-1., -1., -1.])
+
+# global resolution for the grid
+resol3D = [nb, nb, nb]
+
+# Continuous fields (a scalar and a vector)
+
+
+# Example of function to init the scalar field
+def func3D(x, y, z, t):
+    return cos(x) + sin(y) + z
+
+
+# Example of function to init the velocity field
+# 1 everywhere ...
+def func3D_v(x, y, z, t):
+    return 1, 1, 1
+
+scal3D = Field(domain=dom, name='Toto')
+velo3D = Field(domain=dom, name='Tutu', isVector=True)
+
+# --  Penalisation operator, with ghost points --
+# First obstacle : a sphere
+sphere = Sphere(dom, position=[0., 0., 0.], radius=0.5, porousLayers=[0.13])
+# Second obstacle : planes normal to direction 1
+plates = Plates(dom, normal_dir=0, epsilon=0.1)
+# the operator
+penal3D = Penalization([velo3D, scal3D], [sphere, plates], coeff=[1e6, 10],
+                       resolutions={velo3D: resol3D, scal3D: resol3D},
+                       ghosts=[2, 2, 2])
+
+
+# setup :
+# ---> create a topology (including local meshes)
+# ---> discretize velo3D and scal3D on this topology
+penal3D.setUp()
+
+# Monitor : print scal3D to vtk file.
+printer = Printer(fields=[velo3D], frequency=1)
+printer.setUp()
+
+# ---  Get some useful parameters/variables ---
+# - the topology created during penalisation setup -
+topo3D = scal3D.discreteFields.keys()[0]
+# equivalent to :
+# topo3D = dom2.topologies.values()[0]
+# topo3D = velo3D.discreteFields.keys()[0]
+print topo3D
+# Local mesh coordinates:
+mesh = topo3D.mesh.coords
+# indices of points where computation is performed :
+ind = topo3D.mesh.iCompute
+#
+# - discrete fields corresponding to this topology -
+velo_d = velo3D.discreteFields[topo3D]
+scal_d = scal3D.discreteFields[topo3D]
+# numpy arrays :
+tab_velo_x = velo_d.data[0]
+# equivalent to velo_d[0]
+tab_velo_y = velo_d.data[1]
+
+print velo_d
+print scal_d
+
+# Print the whole field, including ghosts:
+print scal_d[0]
+# field values, excluding ghosts:
+print scal_d[0][ind]
+
+# Operators to initialize velocity and scalar fields
+opS = Analytic(scal3D, formula=func3D, topo=topo3D, resolutions=None)
+opV = Analytic(velo3D, formula=func3D_v, topo=topo3D, resolutions=None)
+opS.setUp()
+opV.setUp()
+
+# Simulation parameters:
+simu = Simulation(tinit=0.0, tend=1.0, timeStep=1e-3)
+
+# initialize scal3D and velo3D
+opV.apply(simu)
+opS.apply(simu)
+# Penalize ...
+penal3D.apply(simu)
+printer.apply(simu)
+
+
+# result :
+# for each mpi process rk
+# out_rk_iter_000... initial state
+# out_rk_iter_001... current state