CDMATH/tests/ressources/ConversionScripts/weird_mesh.py

   1 from __future__ import division\r
   2 import MEDLoader as ml\r
   3 from math import pi, cos, sin\r
   4 \r
   5 \r
   6 def getWeirdMesh(xmin, xmax, ymin, ymax, lx, ly, str, func=None, extr=-1, unpo=False):\r
   7     polys = []; ic = 0\r
   8     for iy, y in enumerate(ly[:-1]):\r
   9         for ix, x in enumerate(lx[:-1]):\r
  10             ic += 1; dx = lx[ix + 1] - x; dy = ly[iy + 1] - y\r
  11             c = str[ic % len(str)] if func is None else str[int(min(max(func(xmin + (xmax - xmin) * (x + dx / 2 - lx[0]) / (lx[-1] - lx[0]),\r
  12                                                                              ymin + (ymax - ymin) * (y + dy / 2 - ly[0]) / (ly[-1] - ly[0])), 0), len(str) - 1))]\r
  13             if c == "-":  # deux rectangles (horizontale)\r
  14                 polys.extend([[[x + dx * i, y + dy * (j + k) / 2] for (i, j) in [(0, 0), (0, 1), (1, 1), (1, 0)]] for k in range(2)])\r
  15             elif c == "|":  # deux rectangles (verticale)\r
  16                 polys.extend([[[x + dx * (i + k) / 2, y + j] for (i, j) in [(0, 0), (0, 1), (1, 1), (1, 0)]] for k in range(2)])\r
  17             elif c == "+" or c == "#" or c.isdigit():  # petits carres\r
  18                 n = 2 if c == "+" else 3 if c == "#" else int(c)\r
  19                 polys.extend([[[x + dx * (i + k) / n, y + dy * (j + l) / n] for (i, j) in [(0, 0), (0, 1), (1, 1), (1, 0)]] for k in range(int(n)) for l in range(int(n))])\r
  20             elif c == "x" or c == "*":  # triangles en croix\r
  21                 f = 2 * pi / 8\r
  22                 l = [[x + dx / 2 * (1 + cos(f * k) / max(abs(cos(f * k)), abs(sin(f * k)))), y + dy / 2 * (1 + sin(f * k) / max(abs(cos(f * k)), abs(sin(f * k))))] for k in range(-1, 8, 1 + (c == "x"))]\r
  23                 polys.extend([[l[i], [x + dx / 2, y + dy / 2], l[i + 1]] for i in range(len(l) - 1)])\r
  24             elif c == "/" or c == "\\":\r
  25                 polys.extend([[[x + dx * i, y + dy * (1 - j if c == "/" else j)] for (i, j) in [(0, 1), (1, 0), p]] for p in [(0, 0), (1, 1)]])\r
  26             else:  # par defaut -> carre vide\r
  27                 polys.append([[x + dx * i, y + dy * j] for (i, j) in [(0, 0), (0, 1), (1, 1), (1, 0)]])\r
  28 \r
  29     polys = [[[xmin + (x - lx[0]) * (xmax - xmin) / (lx[-1] - lx[0]), ymin + (y - ly[0]) * (ymax - ymin) / (ly[-1] - ly[0])] for [x, y] in p] for p in polys]\r
  30 \r
  31     mesh = ml.MEDCouplingUMesh("mesh", 2)\r
  32     mesh.allocateCells(len(polys))\r
  33     off = 0\r
  34     for p in polys:\r
  35         mesh.insertNextCell(ml.NORM_POLYGON, len(p), [off + i for i in range(len(p))])\r
  36         off += len(p)\r
  37     mesh.finishInsertingCells()\r
  38 \r
  39     pts = [p for i in range(len(polys)) for p in polys[i]]\r
  40     co = ml.DataArrayDouble(pts, len(pts), 2)\r
  41     mesh.setCoords(co)\r
  42 \r
  43     mesh.changeSpaceDimension(3)\r
  44     mesh.orientCorrectly2DCells([0., 0., -1.], False)\r
  45     mesh.changeSpaceDimension(2)\r
  46     mesh.mergeNodes(1e-8)\r
  47     mesh.conformize2D(1e-8)\r
  48     if unpo: mesh.unPolyze()\r
  49 \r
  50     if extr > 0:\r
  51         mesh.changeSpaceDimension(3)\r
  52         mesh1d = ml.MEDCouplingUMesh("ex", 1)\r
  53         mesh1d.allocateCells(extr)\r
  54         coo = [0., 0., 0.]\r
  55         for i in range(extr):\r
  56             mesh1d.insertNextCell(ml.NORM_SEG2, 2, [i, i + 1])\r
  57             coo.extend([0., 0., float(i + 1) / float(extr)])\r
  58         coo_arr = ml.DataArrayDouble(coo, extr + 1, 3)\r
  59         mesh1d.setCoords(coo_arr)\r
  60         mesh = mesh.buildExtrudedMesh(mesh1d, 0)\r
  61         mesh.setName("mesh")\r
  62 \r
  63     mf, desc, descIndx, revDesc, revDescIndx = mesh.buildDescendingConnectivity()\r
  64     mf.setName("mesh")\r
  65     mm = ml.MEDFileUMesh.New()\r
  66     mm.setMeshAtLevel(0, mesh)\r
  67     mm.setMeshAtLevel(-1, mf)\r
  68 \r
  69     bf = mf.computeCellCenterOfMass()\r
  70     noms_cl = [["left", "right"]]\r
  71     if extr > 0: noms_cl.append(["front", "back"])\r
  72     noms_cl.append(["down", "up"])\r
  73     for i in range(2 + (extr > 0)):  # frontieres en x puis y (puis en z si extr)\r
  74         for j in range(2):  # min puis max\r
  75             g = []\r
  76             for idx, b in enumerate(bf):\r
  77                 if abs(b[i] - [[xmin, xmax], [ymin, ymax], [0., 1.]][i][j]) < 1e-5:\r
  78                     g.append(idx)\r
  79             grp = ml.DataArrayInt64.New(g)\r
  80             grp.setName(noms_cl[i][j])\r
  81             mm.addGroup(-1, grp)\r
  82 \r
  83     return mm\r
  84 \r
  85 if __name__ == "__main__":\r
  86 \r
  87     xmin, xmax, ymin, ymax = 0., 1., 0., 1.\r
  88     n = 4\r
  89 \r
  90     # maillage non conforme\r
  91     getWeirdMesh(xmin, xmax, ymin, ymax, range(n + 1), range(n + 1), "o7/43ox5*2").write("mesh1.med", 2)\r
  92     # maillage strech en x\r
  93     getWeirdMesh(xmin, xmax, ymin, ymax, [i ** 3 for i in range(n + 1)], range(n + 1), "o").write("mesh2.med", 2)\r
  94     # maillage random\r
  95     import random\r
  96     paterns = "o123456789+x*|/\\"\r
  97     f = lambda x, y: random.randint(0, len(paterns))\r
  98     getWeirdMesh(xmin, xmax, ymin, ymax, range(n + 1), range(n + 1), paterns, f).write("mesh3.med", 2)\r
  99     # maillage raffine en fonction d'un champ spatial\r
 100     paterns = "o123456789"\r
 101     f = lambda x, y: 9 * (x ** 2 + y ** 2)\r
 102     getWeirdMesh(xmin, xmax, ymin, ymax, range(n + 1), range(n + 1), paterns, f).write("mesh4.med", 2)\r
 103     # maillage 3D\r
 104     getWeirdMesh(xmin, xmax, ymin, ymax, range(n + 1), range(n + 1), "o7/43ox5*2", extr=4).write("mesh5.med", 2)\r
 105     # maillage raffinement local type fvca\r
 106     f = lambda x, y: 2 * (int(x > 0.5 and y < 0.5) + int(x > 0.75 and y < 0.25))\r
 107     getWeirdMesh(xmin, xmax, ymin, ymax, range(n + 1), range(n + 1), paterns, f).write("mesh6.med", 2)\r
 108     # maillage checkerboard\r
 109     m = n - (n % 2) + 1\r
 110     getWeirdMesh(xmin, xmax, ymin, ymax, range(m + 1), range(m + 1), "o+").write("mesh7.med", 2)\r