From ae6e392692465eb9855fcab1c3d843c97dc3dbb3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Franck=20P=C3=A9rignon?= <franck.perignon@imag.fr> Date: Wed, 17 Jul 2013 14:59:36 +0000 Subject: [PATCH] Example to test post-treatment and visu (for Caroline) --- Examples/testVisu.py | 120 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 120 insertions(+) create mode 100755 Examples/testVisu.py diff --git a/Examples/testVisu.py b/Examples/testVisu.py new file mode 100755 index 000000000..62d284b5a --- /dev/null +++ b/Examples/testVisu.py @@ -0,0 +1,120 @@ +# +# 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 -- GitLab