[FEM] 格子分割

mesh.py

# -*- coding: utf-8 -*-

## 必要なモジュールをインポート
import os
import sys
import time

## PySide系モジュール
from PySide.QtGui import *
from PySide.QtCore import *

## 演算系モジュール
import numpy as np


# =============================================================================
## GUIの構築
class GUI(QWidget):
    def __init__(self, parent=None):
        QWidget.__init__(self, parent)
        self.setObjectName('name')

        # self.setFixedSize(LENGTH, LENGTH)

        # クラスプロパティ
        # -----------------
        # 画面バッファ
        self.pixmap = QPixmap(self.size())

        # 初期画面の準備
        # --------------
        # 画面バッファの初期化
        self.refreshPixmap()
        # グリッドの表示
        painter = QPainter(self.pixmap)
        # self.drawGrid(painter)
        self.drawGeoVertics(painter)

        # メインループの準備と開始
        # -------------------------
        # self.timer = QTimer()
        # self.timer.timeout.connect(self.mainloop)
        # self.timer.start(self.interval_time)
        self.update()

    # ************************************************************* #
    # メインループ
    # ************************************************************* #
    def mainloop(self):
        """
        アニメーションのメインループ
        アルゴリズムの時間更新等はここで行う
        """
        # 画面更新
        self.update()

    def paintEvent(self, *args, **kwargs):
        """画面描画関数。self.pixmapにバッファした画像を描画する"""
        # QPainterを生成
        painter = QStylePainter(self)
        # QPainterでバッファに準備したデータを描画
        painter.drawPixmap(0, 0, self.pixmap)

    # ************************************************************* #
    # 描画系補助関数
    # ************************************************************* #
    def drawGrid(self, painter):
        """ マップを表示する関数
         (AgentやGoal,Startはメイン内で描画する"""
        w = self.size().width()
        h = self.size().height()
        dx = w / 10.
        dy = h / 10.
        for xn in range(10):
            for yn in range(10):
                painter.setBrush(QBrush(QColor(100, 100, 100), Qt.SolidPattern))
                painter.drawRect(xn * dx, yn * dy, dx, dy)

    def drawGeoVertics(self, painter):

        self.padding = self.size().width() / 10.
        xmax = self.size().width() - self.padding
        ymax = self.size().height() - self.padding
        xmin = xmax / 10.
        ymin = ymax / 10.
        p1 = np.array((xmin, ymin))
        p2 = np.array((xmax, ymin))
        p3 = np.array((xmax/3., ymax))
        p4 = np.array((xmin, ymax))
        points = [p1, p2, p3, p4]

        geo_lines = []
        geo_lines.append(GeoLine(points[0], points[1], div=25))
        geo_lines.append(GeoLine(points[1], points[2], div=25))
        geo_lines.append(GeoLine(points[2], points[3]))
        geo_lines.append(GeoLine(points[3], points[0]))


        # X座標の分割数を取得
        xi = 0
        if geo_lines[xi].div != geo_lines[xi + 2].div:
            if geo_lines[xi].div > geo_lines[xi + 2].div:
                divx = geo_lines[xi].div
            else:
                divx = geo_lines[xi + 2].div
        else:
            divx = geo_lines[xi].div

        # Y座標の分割数を取得
        yi = 2
        if geo_lines[yi].div != geo_lines[1].div:
            if geo_lines[yi].div > geo_lines[1].div:
                divy = geo_lines[yi].div
            else:
                divy = geo_lines[1].div
        else:
            divy = geo_lines[yi].div

        I = divx + 1
        J = divy + 1
        # 正規化座標を作る
        c = [[np.asarray([i * (1. / divx), j * (1. / divy)]) for j in range(J)] for i in range(I)]

        # 境界条件
        c[0][0] = p1
        c[divx][0] = p2
        c[divx][divy] = p3
        c[0][divy] = p4

        # 実座標に変換
        for xi in range(I):
            o1 = c[xi][0] = (c[divx][0] - c[0][0]) * (xi / float(divx)) + c[0][0]
            o2 = c[xi][divy] = (c[divx][divy] - c[0][divy]) * (xi / float(divx)) + c[0][divy]
            o12 = o2 - o1
            for yi in range(J):
                c[xi][yi] = o12 * (float(yi) / divy) + o1

        # エレメントを作る(divx x divy)
        elements = []
        for xi in range(divx):
            for yi in range(divy):
                d = [c[xi][yi], c[xi + 1][yi], c[xi + 1][yi + 1], c[xi][yi + 1]]
                elements.append(d)

        # エレメントがちゃんとできているか表示
        for element in elements:
            polygon = QPolygon()
            for p in element:
                polygon << QPoint(p[0], p[1])
            painter.drawPolygon(polygon)

        import csv
        with open('mesh.csv', 'w') as f:
            writer = csv.writer(f, lineterminator='\n')
            writer.writerows(elements)
            for e in elements:
                out = ['id', e]
                writer.writerow(out)
        # 画像の保存
        self.pixmap.save('pixmap.bmp')

    # ************************************************************* #
    # その他Qt関連補助関数
    # ************************************************************* #
    def refreshPixmap(self):
        """
        画面バッファの初期化関数
        """
        # 画面バッファの初期化
        self.pixmap = QPixmap(self.size())
        # 画面を塗りつぶし (おまじない)
        self.pixmap.fill(self, 0, 0)
        self.pixmap.fill(Qt.white)
        # ぺインターの生成 (おまじない)
        painter = QPainter(self.pixmap)
        # ぺインターによる初期化 (おまじない)
        painter.initFrom(self)
        pass

    def keyPressEvent(self, event):
        """子供ウィジェットの場合は、無効"""
        e = event.key()

        if e == Qt.Key_Q:
            self.close()
        else:
            pass
        self.update()

    def mousePressEvent(self, event):
        print event.pos()
        print event.posF()
        print event.button()
    pass


# =============================================================================
## GUIの起動
def main():
    app = QApplication(sys.argv)
    wnd = GUI()
    wnd.show()
    sys.exit(app.exec_())


class GeoLine(object):
    def __init__(self, p1, p2, div=3):
        """
        :param p1: ndarray(1x2)
        :param p2: ndarray(1x2)
        :param div: int, 分割数
        """
        self.p1 = p1
        self.p2 = p2
        self.div = div
        self.p12 = self.p2 - self.p1
        self.unit_p12 = self.p12 / np.linalg.norm(self.p12)

if __name__ == '__main__':
    main()


# -----------------------------------------------------------------------------
# EOF
# -----------------------------------------------------------------------------