AEi(解析連続実行)

AEi(Analysis Execution interpreter )は if 文制御、for 文ループを備え、簡単な記述でCSVファイルからデータを読み込み、変数に代入することができる独自言語のインタープリタです。

この変数を、解析データの元となるテンプレートデータに展開することで解析データを作成し、解析を実行します。

簡単な処理の流れを以下に記載します。

1)テンプレートデータ
解析データの元となるテンプレートデータに「$◯◯◯」という$変数名で変化させていくデータを記述する。

2)CSVファイルの数値を$変数に代入する
この変化させていく$変数にCSVファイルから読み込んだ数値や文字列などを代入し、それらの値でテンプレートデータを書き換え、解析用のデータを作成する。

3)解析を実行する
作成されたデータを解析コードに渡して解析を実行する。

進化計算ライブラリDEAPでレイリー減衰パラメータ探査

前回の「遺伝的アルゴリズムによるレイリー減衰のパラメータ探査」で試した、一定減衰に相当するレイリー減衰のパラメータ探査を、今回はPythonの進化計算ライブラリのDEAPを使ってやってみます。

プログラムは、データサイエンス教育総合研究所のサイトを参考にさせていただきました。参考にしたというより、目的関数のところを一定減衰Shakeとレイリー減衰応答解析の差分2乗和とする関数に置き換えただけです。

前回は、遺伝的アルゴリズムをプログラミングしてパラメータを探しましたが、DEAPを使うとアルゴリズムの詳細をプログラミングしなくてもパラメータ探査できます。

実行結果は以下のとおり。

最も良い個体は [2.572249424150136, 6.2771604118935755]で、そのときの目的関数の値は (7.337222985431525,)

上記の個体の数値がレイリー減衰パラメータである二つの振動数になります。

このレイリー減衰のパラメータで応答計算し、評価位置で応答スペクトルを比較したものが下図です。

DEAPによるレイリー減衰パラメータ探査

一定減衰とレイリー減衰の応答性状がよく合致しています。

参考までに、全く的外れなレーリー減衰のパラメータとして、0.02Hzと50Hzを指定した場合の比較は下図になります。

不適切なレイリー減衰パラメータの応答スペクトル比較

以下のプログラムでは、この0.02Hzから50Hzを遺伝子の範囲としています。

#個体の各遺伝子を決めるために使用
import random
import numpy as np
import pandas as pd

#DEAPの中にある必要なモジュールをインポート
from deap import base
from deap import creator
from deap import tools
from deap import algorithms

from FreeField import FreeField

#最小化問題として設定(-1.0で最小化、1.0で最大化問題)
#面食らうかもしれませんが、目的が最小化の場合は以下のように設定します。
creator.create("FitnessMin", base.Fitness, weights=(-1.0,))

#個体の定義(list型と指定しただけで、中身の遺伝子は後で入れる)
#これも面食らうかもしれませんが、個体を準備したと思ってください。
creator.create("Individual", list, fitness=creator.FitnessMin)

#目的関数の定義。必ずreturnの後に,をつける
#レイリー減衰の応答解析とターゲットとなる一定減衰Shakeの差分2乗和を計算する
def resp(individual):
    #一定減衰Shakeのターゲットスペクトルの読み込み
    df_target_spec = pd.read_csv("shake.csv")
    #レイリー減衰応答解析のスペクトル計算
    FreeField_spec = np.array([FreeField(individual).resp()])
    object = np.sum((df_target_spec["shake"] - FreeField_spec[0]) ** 2)
    return object,

#各種関数の設定を行います
#交叉、選択、突然変異などには、DEAPのToolbox内にある関数を利用
toolbox = base.Toolbox()
#random.uniformの別名をattribute関数として設定。各個体の遺伝子の中身を決める関数(各遺伝子は0.02Hz~50Hzのランダムな値)
toolbox.register("attribute", random.uniform, 0.02,50)
#individualという関数を設定。それぞれの個体に含まれる2個の遺伝子をattributeにより決めるよ、ということ。
toolbox.register("individual", tools.initRepeat, creator.Individual, toolbox.attribute, 2)
#集団の個体数を設定するための関数を準備
toolbox.register("population", tools.initRepeat, list, toolbox.individual)
#トーナメント方式で次世代に子を残す親を選択(tornsizeは各トーナメントに参加する個体の数)
toolbox.register("select", tools.selTournament, tournsize=5)
#交叉関数の設定。ブレンド交叉という手法を採用
toolbox.register("mate", tools.cxBlend,alpha=0.2)
#突然変異関数の設定。indpbは各遺伝子が突然変異を起こす確率。muとsigmaは変異の平均と標準偏差
toolbox.register("mutate", tools.mutGaussian, mu=[0.0, 0.0], sigma=[20.0, 20.0], indpb=0.2)
#評価したい関数の設定(目的関数のこと)
toolbox.register("evaluate", resp)

#以下でパラメータの設定
#今回は最も単純な遺伝的アルゴリズムの手法を採用

#乱数を固定
random.seed(64)
#何世代まで行うか
NGEN = 30
#集団の個体数
POP = 10
#交叉確率
CXPB = 0.9
#個体が突然変異を起こす確率
MUTPB = 0.1

#集団は80個体という情報の設定
pop = toolbox.population(n=POP)
#集団内の個体それぞれの適応度(目的関数の値)を計算
for individual in pop:
    individual.fitness.values = toolbox.evaluate(individual)
#パレート曲線上の個体(つまり、良い結果の個体)をhofという変数に格納
hof = tools.ParetoFront()

#今回は最も単純なSimple GAという進化戦略を採用
algorithms.eaSimple(pop, toolbox, cxpb=CXPB, mutpb=MUTPB, ngen=NGEN, halloffame=hof)

#最終的な集団(pop)からベストな個体を1体選出する関数
best_ind = tools.selBest(pop, 1)[0]
#結果表示
print("最も良い個体は %sで、そのときの目的関数の値は %s" % (best_ind, best_ind.fitness.values))

以下、関連リンク
DEAPの本家サイト

教師なし機械学習で地震波をクラスタリングする

Pythonの機械学習ライブラリ、scikit-learnの教師なし機械学習で、地震波をクラスタリングしてみました。

地震波は、K-NETから2011年3月11日の地震波をダウンロードして、予めランダムに選択した地震波の応答スペクトルを計算しておき、これをデータにして5つのクラスタに分けます。

以下がプログラムです。

import glob
import numpy as np
import pandas as pd
from sklearn.cluster import KMeans
import matplotlib.pyplot as plt
import geopandas as gpd


if __name__ == '__main__':

    # 予め計算しておいた各記録波の応答スペクトルを読込む
    files = glob.glob('*.spec.csv')
    feature = np.array([])
    for file in files:
        df = pd.read_csv(file, index_col=None)
        feature = np.append(feature, df['acc'])
        # feature = np.append(feature, df['acc']/df['acc'][0])

    # 1次元numpy配列を2次元に変換する(応答スペクトルの計算ポイント数300)
    feature = feature.reshape(-1, 300)  

    # 5種類のグループにクラスタリングする
    model = KMeans(n_clusters=5, init='random').fit(feature)

    # 学習結果のラベルを取得する
    labels = model.labels_
    
    # 応答スペクトルをクラスタ毎に色分けしてプロットする
    colors = ["red", "green", "blue", "magenta", "cyan"]
    plt.grid()
    plt.xscale('log')
    plt.xlabel('Period(sec)')
    plt.ylabel('Response Acceleration Spectrum')
    for label, v, file in zip(labels, feature, files):
        plt.plot(df['period(sec)'], v, color=colors[label])
    plt.show()


    # 各観測点の情報を読込む
    geo = pd.read_csv('geo.csv', index_col='name')

    # クラスタ毎に色分けして地図上にプロットする
    jpnShp = gpd.read_file('./gm-jpn-all_u_2_2/coastl_jpn.shp')
    ax = jpnShp.plot(figsize=(10, 20), color='black', linewidth=.5)
    for label, file in zip(labels, files):
        plt.scatter(geo.at[file[:6],'Long'], geo.at[file[:6],'Lat'], marker='o', s=10, color=colors[label], alpha=1.0, linewidths = 1)
    plt.show()

データを準備したり、結果を描画するための処理が多いのですが、機械学習自体は以下のたった1行で完結します。

model = KMeans(n_clusters=5, init='random').fit(feature)

このプログラムでグループ分けされた応答スペクトルを描画すると以下のようになります。

応答スペクトルのクラスタリング(1)

少し分かりづらいですが、greenのグループが応答スペクトルが大きいグループ、blueのグループが応答スペクトルが小さいグループで、blue < cyan < red < magenta < green の順で大きくなっています。

これらの地震波を観測点ごとに、上図と同じ色分けで日本地図に描画すると以下のようになります。
(地図データは国土地理院のデータを使用しています。)

観測点プロット(1)

容易に想像できますが、震源に近い観測点で大きな応答スペクトルになっていて、遠いところで小さな応答スペクトルになっています。

さて、地震波はその特徴として周期特性がありますが、この周期特性を元にしてクラスタリングする為に、応答スペクトルを地震波の最大値で基準化してグループ分けしてみます。

地震波の時間刻みが0.01秒であれば、解像度は50Hzで、応答スペクトルの最短周期の0.02秒が地震波の最大値になるので、この値で応答スペクトル全体を徐算します。

プログラムは以下の1行のコメントを外して入れ替えます。

# feature = np.append(feature, df['acc'])
feature = np.append(feature, df['acc']/df['acc'][0])

この結果を描画したものが以下です。

応答スペクトルのクラスタリング(2)

redが長周期成分が卓越しているグループ、greenが短周期成分が卓越しているグループで、green < cyan < magenta < blue < red の順で卓越している周期成分が長周期寄りになります。

最初のグループ分けと同様にこれらを日本地図に描画したものが以下です。

観測点プロット(2)

最初のケースと比較すると、各グループの際立った地域特性は発見できません。これは、観測波の周期特性はその点の地盤の特性が大きく影響するためだと思います。

北海道や関西に長周期寄りの観測波が多いのは何か理由があるのかもしれません。

scikit-learnの公式サイトには豊富なサンプルがあるので、直ぐに機械学習を試すことができます。

※2021/10/30追記
長周期成分が卓越している記録波が遠方に多いのは、長周期成分は遠方まで到達しやすいからではないかとのコメントをいただきました。

ラージマス法の周波数応答解析

伝達関数

complex_K2はVoigtモデルを想定した複素剛性

$$k^* =(1+2ih)k$$

complex_Kは位相特性を改良した複素剛性です。

$$k^* =k(1−2h^2 +2ih\sqrt{1−h^2})$$

式から分かるとおり、h<<1.0なら同等です。

import numpy as np
import math
from matplotlib import pyplot as plt


def complex_K(k, h):
    return complex(k*(1.0-2.0*h**2.0), k*(2.0*h*math.sqrt(1.0-h**2.0)))


def complex_K2(k, h):
    return complex(k, 2.0*k*h)


if __name__ =="__main__":

    m1 = 100.0
    m2 = 100.0
    m3 = 100.0E6
    k1 = 20000
    k2 = 20000
    k3 = 0.00
    h1 = 0.05
    h2 = 0.05
    h3 = 0.0

    M = np.array([[m1, 0, 0],
                [0, m2, 0],
                [0, 0, m3]])

    cmp_k1 = complex_K(k1, h1)
    cmp_k2 = complex_K(k2, h2)
    cmp_k3 = complex_K(k3, h3)

    K = np.array([[cmp_k1, -cmp_k1, 0],
                [-cmp_k1, cmp_k1+cmp_k2, -cmp_k2],
                [0,-cmp_k2,cmp_k2+cmp_k3]])

    freq = np.arange(0.01, 5, 0.01)

    X = []
    for f in freq:
        omega = 2.0 * np.pi * f
        omega2 = omega**2.0
        cmp_K = K - omega2 * M
        inv_cmp_K = np.linalg.inv(cmp_K)
        cmp_dsp = np.dot(inv_cmp_K, np.array([0, 0, m3]))
        cmp_acc = -omega2 * cmp_dsp
        X.append(list(abs(cmp_acc)))

    X1 = [r[0] for r in X]
    X2 = [r[1] for r in X]
    X3 = [r[2] for r in X]
    plt.plot(freq, X1, label="Node 1")
    plt.plot(freq, X2, label="Node 2")
    plt.plot(freq, X3, label="Node 3")
    plt.title("Transfer function")
    plt.xlabel("Frequency(Hz)")
    plt.ylabel("Amplitude")
    plt.legend()
    plt.grid()
    plt.show()
    plt.close()    

遺伝的アルゴリズムによるレイリー減衰のパラメータ探査

遺伝的アルゴリズム(genetic algorithm、略称:GA)は、古くからある近似解探査法の一つで、パラメータを生物の遺伝子に見立て、良い遺伝子を選択して交叉と突然変異を行いながら、世代交代を繰り返して最適なパラメータを探しあてる手法です。

有名なところでは、新幹線のN700系のフロントノーズの設計に使われました。

以下は一定減衰(一次元波動論)の応答スペクトルをターゲットにして、レイリー減衰(質点系)の振動数パラメータを遺伝子としてGAで最適な解を探査したものです。

図はGAで計算した各世代の最優良パラメータにおける応答スペクトルの進化をアニメgifで表しています。

各世代の応答スペクトルの進化

モデルの概要は以下です。

解析モデル概要

このアプローチでは遺伝子が二つで、あっという間に優良解が見つかってしまいます。

しかし、各部別レイリー減衰で物性ごとにレイリー減衰を設定する場合は遺伝子の数を増やすことになるので、こういった探査法が有益だと思います。

また、水平動だけでなく、上下動も同時に評価して最適解を求めたいときは、優良遺伝子の判定を工夫することで同様に探査が可能です。

ところで、質点系レイリー減衰のモデルは、外部減衰の接続先が半無限性を模擬するダンパー端部の固定点です。この固定点と、質点系モデルとの相対速度は工学的には無意味なので、減衰が適切に評価されませんが、このアプローチは、そうしたモデルの矛盾も吸収してレイリー減衰のパラメータを決定していることになります。

以下、関連リンク

# coding:utf-8

import numpy as np
from numpy import random as rand
import pandas as pd
from FreeField import FreeField
import matplotlib.pyplot as plt

from PIL import Image
import glob

class GeneticAlgorithm(object):
    """遺伝的アルゴリズム"""

    def __init__(self):
        """初期化"""

        """定数"""
        # 染色体数
        self.N = 2
        # 個体数
        self.M = 10
        # エリート選択個数
        self.n_elite = 3
        self.p_mutation = 0.10

        """変数"""
        # 第g世代評価降順インデックス
        self.val_sorted_index = []
        # 選択された遺伝子のインデックス
        self.selected_index = []
        # 次世代遺伝子
        self.next_generation = []
        # 各世代の最良応答スペクトル
        self.best_spec = []

        """初期値"""
        # 第g世代遺伝子群
        self.X_g = np.array([[rand.uniform(0.2,50.0),       # f1
                             rand.uniform(0.2,50.0)]        # f2
                             for i in range(self.M)])

        """ターゲットスペクトル"""
        self.df_target_spec = pd.read_csv("shake.csv")


    def assess(self):
        """評価"""
        # ターゲットスペクトルとの差分2乗和を評価値とする
        FreeField_spec = np.array([FreeField(x).resp() for x in self.X_g])
        R = []
        for spec in FreeField_spec:
            R.append(np.sum((self.df_target_spec["shake"] - spec) ** 2))
        R = np.array(R)
        self.val_sorted_index = np.argsort(R)[::1]
        self.best_spec.append(FreeField_spec[self.val_sorted_index[0]])


    def select(self):
        """選択"""
        # エリート選択で N 個
        elite_index = self.val_sorted_index[0:self.n_elite]
        self.selected_index = np.r_[elite_index]
        print("selected_index",self.selected_index)


    def crossover1p(self):
        """交叉"""
        ng = []
        for i in range(self.M - self.n_elite):
            p_index = rand.choice(self.selected_index, 2, replace=False)
            new_g = np.r_[self.X_g[p_index[0]][0], self.X_g[p_index[1]][1]]
            ng.append(new_g)

        ng = np.array(ng)
        self.next_generation = np.r_[ng, self.X_g[self.selected_index[0:self.n_elite]]]


    def mutation(self):
        """突然変異"""
        # 確率的に選ばれた染色体をランダム値に変更
        # 突然変異を起こす染色体の数
        n_mutation = rand.binomial(n=self.M * self.N, p=self.p_mutation)

        # 突然変異
        for i in range(n_mutation):
            m = rand.randint(self.M)
            n = rand.randint(self.N)
            self.next_generation[m][n] = rand.uniform(0.2,50.0)


    def alternation(self):
        """世代交代"""
        self.X_g = self.next_generation


if __name__ == '__main__':

    """設定"""
    # 計算世代
    G = 50

    """計算"""
    ga = GeneticAlgorithm()
    for g in range(G):
        # 評価
        ga.assess()
        # 選択
        ga.select()
        # 交叉
        ga.crossover1p()
        # 突然変異
        ga.mutation()
        # 世代交代
        ga.alternation()

    for i, spec in enumerate(ga.best_spec):
        # グラフ描画
        plt.plot(ga.df_target_spec["t"], ga.df_target_spec["shake"], label="Target")
        plt.plot(ga.df_target_spec["t"], spec, label="GA")
        plt.title("GA = "+str(i))
        plt.xlabel("Period(sec)")
        plt.ylabel("Acceleration Response Spectrum $(m/s^{2})$")
        plt.xscale("log")
        plt.ylim(0.0, 20.0)
        plt.legend()
        plt.grid()
        png = "spec_{:0>2d}.png".format(i)
        plt.savefig(png)
        plt.clf()

    # GIF 作成
    files = sorted(glob.glob('*.png'))  
    images = list(map(lambda file : Image.open(file) , files))
    images[0].save('GA.gif' , save_all = True , append_images = images[1:] , duration = 1000)
# coding:utf-8
import numpy as np
from numpy.lib.function_base import append
import pandas as pd
import math

class FreeField(object):

    def __init__(self, x):
        self.f1 = x[0]
        self.f2 = x[1]


    # 時刻歴応答解析
    def resp(self):
        f1 = self.f1
        f2 = self.f2

        model = pd.read_csv("model.csv")
        h1 = model["h"][0]
        h2 = model["h"][0]

        size = len(model["M"])
        # 質量マトリクス
        m = np.diag(list(map(lambda x: x / 9.80665, model["M"])))
        # 剛性マトリクス
        k = np.zeros((size, size))
        kdata = size -1
        for i in range(kdata):
            i1 = i
            i2 = i + 1
            k[i1][i1] = k[i1][i1] + model["K"][i]
            k[i2][i2] = k[i2][i2] + model["K"][i]
            k[i1][i2] = k[i1][i2] - model["K"][i]
            k[i2][i1] = k[i2][i1] - model["K"][i]
        # 減衰マトリクス(レイリー減衰)
        w1 = 2.0 * np.pi * f1
        w2 = 2.0 * np.pi * f2
        alpha = (2.0 * w1 * w2 * (h2 * w1 - h1 * w2)) / (w1 ** 2 - w2 ** 2)
        beta = (2.0 * (h1 * w1 - h2 * w2)) / (w1 ** 2 - w2 ** 2)
        c = np.array(alpha * m + beta * k)

        # 底面ダンパー(2E入力用)
        c[kdata][kdata] = c[kdata][kdata] + model["C"][0]

        # 入力波読み込み
        df = pd.read_csv(model["accfile"][0])
        acc0 = df["acc"]
        dt = 0.01
        beta = 0.25
        unit_vector = np.ones(size)
        pre_acc0 = 0.0
        time = 0.0
        dis = np.zeros(size)
        vel = np.zeros(size)
        acc = np.zeros(size)
        ddis = np.zeros(size)
        dvel = np.zeros(size)
        dacc = np.zeros(size)
        acc_history = {}
        for i in range(0, size):
            acc_history[i] = []
        time_history = []
        # Newmarkβ法による数値解析(増分変位による表現)
        for i in range(0, len(acc0)):
            kbar = k + (1.0/(2.0*beta*dt)) * c + (1.0/(beta*dt**2.0)) * m
            dp1 = -1.0 * m.dot(unit_vector) * (acc0[i] - pre_acc0)
            dp2 = m.dot((1.0/(beta*dt))*vel + (1.0/(2.0*beta))*acc)
            dp3 = c.dot((1.0/(2.0*beta))*vel + (1.0/(4.0*beta)-1.0)*acc*dt)
            dp = dp1 + dp2 + dp3
            ddis = np.linalg.inv(kbar).dot(dp)
            dvel = (1.0/(2.0*beta*dt))*ddis - (1.0/(2.0*beta))*vel - ((1.0/(4.0*beta)-1.0))*acc*dt
            dacc = (1.0/(beta*dt**2.0))*ddis - (1.0/(beta*dt))*vel - (1.0/(2.0*beta))*acc
            dis += ddis
            vel += dvel
            acc += dacc
            acc_abs = acc + [acc0[i] for n in range(0,size)]
            [acc_history[i].append(x) for i, x in enumerate(acc_abs)]
            time_history.append(time)
            time += dt
            pre_acc0 = acc0[i]

        # 評価対象として最上部の応答スペクトルを計算する
        spec = self.duhm_spec(acc_history[0])
        return(spec)


    def duhm_spec(self, acc):
        spec = np.empty(0)
        dlta = (math.log(5.0) - math.log(0.02)) / float(300 - 1)
        for i in range(1, 300 + 1):
            tmp = math.log(0.02) + dlta * float(i - 1)
            period = math.exp(tmp)
            zmax = self.duhm(0.01, 0.05, period, acc)
            spec = np.append(spec, zmax)
        return(spec)


    def duhm(self, dt, h, T, a):
        """duhamel integral"""
        w = 2.0 * math.pi / T
        w2 = w * w
        hw = h * w
        wd = w * math.sqrt( 1.0 - h * h )
        wdt = wd * dt
        e = math.exp( -hw * dt )
        cwdt = math.cos( wdt )
        swdt = math.sin( wdt )
        e11 = e * ( cwdt - hw * swdt / wd )
        e12 = -e * ( wd * wd + hw * hw ) * swdt / wd
        e21 = e * swdt / wd
        e22 = e * ( cwdt + hw * swdt / wd )
        ss = - hw * swdt - wd * cwdt
        cc = - hw * cwdt + wd * swdt
        s1 = ( e * ss + wd ) / w2
        c1 = ( e * cc + hw ) / w2
        s2 = ( e * dt * ss + hw * s1 + wd * c1 ) / w2
        c2 = ( e * dt * cc + hw * c1 - wd * s1 ) /w2
        s3 = dt * s1 - s2
        c3 = dt * c1 - c2
        g11 = ( -hw * s3 + wd * c3 ) / wdt
        g12 = ( -hw * s2 + wd * c2 ) / wdt
        g21 = s3 / wdt
        g22 = s2 / wdt
        dx = 0.0
        x = 0.0
        accmax = 0.0
        for m in range(1, len(a)):
            dxf = dx
            xf = x
            ddym = a[ m ]
            ddyf = a[ m - 1 ]
            dx = e11 * dxf + e12 * xf + g11 * ddym + g12 * ddyf
            x = e21 * dxf + e22 * xf + g21 * ddym + g22 * ddyf
            ddx = 2.0 * hw * dx + w2 * x
            if abs(ddx) > abs(accmax):
                accmax = abs(ddx)

        return(accmax)

使用したデータも含め、プログラムはこちらにも置いています。

Femap API Results Browsing Object の SetColumns メソッドの不適切な挙動

FEMAP 2020.1 のAPIでは、以下のプログラムでFEMAPがフリーズします。

# coding: utf-8

import pythoncom
import Pyfemap

existObj = pythoncom.connect(Pyfemap.model.CLSID)
app = Pyfemap.model(existObj)

elLIST = [1,  2,  3,  4,  5,  6,  7]
pLIST0 = [1.0,2.0,3.0,4.0,5.0,6.0,7.0]
pLIST1 = [1.0,2.0,3.0,4.0,5.0,6.0,7.0]
pLIST2 = [1.0,2.0,3.0,4.0,5.0,6.0,7.0]
pLIST3 = [1.0,2.0,3.0,4.0,5.0,6.0,7.0]
pLIST4 = [1.0,2.0,3.0,4.0,5.0,6.0,7.0]
pLIST5 = [1.0,2.0,3.0,4.0,5.0,6.0,7.0]
pLIST6 = [1.0,2.0,3.0,4.0,5.0,6.0,7.0]
pLIST7 = [1.0,2.0,3.0,4.0,5.0,6.0,7.0]
pLIST8 = [1.0,2.0,3.0,4.0,5.0,6.0,7.0]

setID = 10
outSet = app.feOutputSet
outSet.title = "test"
outSet.value = 0
outSet.analysis = 0
outSet.Put (setID)
RBO = app.feResults
pLIST_ALL = [pLIST0, pLIST1, pLIST2, pLIST3, pLIST4, pLIST5, pLIST6, pLIST7, pLIST8]
ttl = "Corner Pressure Face 1 Set " + str(setID)
_, nCol = RBO.AddElemWithCornerColumnsV2(setID, 24000000, 24000001, 24000002, 24000003, 24000004, 24000005, 24000006, 24000007, 24000008, ttl, 4, True)
_ = RBO.SetColumns(len(nCol), nCol, len(elLIST), elLIST, pLIST_ALL)

RBO.Save()
_ = outSet.Put(-1)

しかし、以下のプログラムなら正常に動作します。

# coding: utf-8

import pythoncom
import Pyfemap

existObj = pythoncom.connect(Pyfemap.model.CLSID)
app = Pyfemap.model(existObj)

elLIST = [1,  2,  3,  4,  5,  6,  7,  8]
pLIST0 = [1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0]
pLIST1 = [1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0]
pLIST2 = [1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0]
pLIST3 = [1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0]
pLIST4 = [1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0]
pLIST5 = [1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0]
pLIST6 = [1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0]
pLIST7 = [1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0]
pLIST8 = [1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0]

setID = 10
outSet = app.feOutputSet
outSet.title = "test"
outSet.value = 0
outSet.analysis = 0
outSet.Put (setID)
RBO = app.feResults
pLIST_ALL = [pLIST0, pLIST1, pLIST2, pLIST3, pLIST4, pLIST5, pLIST6, pLIST7, pLIST8]
ttl = "Corner Pressure Face 1 Set " + str(setID)
_, nCol = RBO.AddElemWithCornerColumnsV2(setID, 24000000, 24000001, 24000002, 24000003, 24000004, 24000005, 24000006, 24000007, 24000008, ttl, 4, True)
_ = RBO.SetColumns(len(nCol), nCol, len(elLIST), elLIST, pLIST_ALL)

RBO.Save()
_ = outSet.Put(-1)

つまり、SetColumns メソッドは、値を設定する要素数が7以下の場合はFEMAPがフリーズし、8以上の場合は正常に動作する不適切な挙動があります。

これを回避するためには、以下のプログラムのように、SetColumn メソッドを使って、各列ごとに値を設定する必要があります。

# coding: utf-8

import pythoncom
import Pyfemap

existObj = pythoncom.connect(Pyfemap.model.CLSID)
app = Pyfemap.model(existObj)

elLIST = [1,  2,  3,  4,  5,  6,  7]
pLIST0 = [1.0,2.0,3.0,4.0,5.0,6.0,7.0]
pLIST1 = [1.0,2.0,3.0,4.0,5.0,6.0,7.0]
pLIST2 = [1.0,2.0,3.0,4.0,5.0,6.0,7.0]
pLIST3 = [1.0,2.0,3.0,4.0,5.0,6.0,7.0]
pLIST4 = [1.0,2.0,3.0,4.0,5.0,6.0,7.0]
pLIST5 = [1.0,2.0,3.0,4.0,5.0,6.0,7.0]
pLIST6 = [1.0,2.0,3.0,4.0,5.0,6.0,7.0]
pLIST7 = [1.0,2.0,3.0,4.0,5.0,6.0,7.0]
pLIST8 = [1.0,2.0,3.0,4.0,5.0,6.0,7.0]

setID = 10
outSet = app.feOutputSet
outSet.title = "test"
outSet.value = 0
outSet.analysis = 0
outSet.Put (setID)
RBO = app.feResults
pLIST_ALL = [pLIST0, pLIST1, pLIST2, pLIST3, pLIST4, pLIST5, pLIST6, pLIST7, pLIST8]
ttl = "Corner Pressure Face 1 Set " + str(setID)
_, nCol = RBO.AddElemWithCornerColumnsV2(setID, 24000000, 24000001, 24000002, 24000003, 24000004, 24000005, 24000006, 24000007, 24000008, ttl, 4, True)
_ = RBO.SetColumn(nCol[0], len(elLIST), elLIST, pLIST0)
_ = RBO.SetColumn(nCol[1], len(elLIST), elLIST, pLIST1)
_ = RBO.SetColumn(nCol[2], len(elLIST), elLIST, pLIST2)
_ = RBO.SetColumn(nCol[3], len(elLIST), elLIST, pLIST3)
_ = RBO.SetColumn(nCol[4], len(elLIST), elLIST, pLIST4)
_ = RBO.SetColumn(nCol[5], len(elLIST), elLIST, pLIST5)
_ = RBO.SetColumn(nCol[6], len(elLIST), elLIST, pLIST6)
_ = RBO.SetColumn(nCol[7], len(elLIST), elLIST, pLIST7)
_ = RBO.SetColumn(nCol[8], len(elLIST), elLIST, pLIST8)

RBO.Save()
_ = outSet.Put(-1)

Fortranで作成したdylibをPythonで呼び出す

既存のFortran資産をPythonから利用する手段として有効です。
以下はFortranで作成したヘロンの公式のサブルーチンをPythonから使用する例です。

環境は以下のとおりです。

OS
 macOS 10.13.6
fortran コンパイラ
 GNU Fortran (GCC) 4.9.2 20141029 (prerelease)
python
 Python 3.7.3

Fortranサブルーチンは以下です。

subroutine heron(a,b,c,s)
    implicit none
    real(8),intent(in) :: a,b,c
    real(8),intent(out):: s
    real(8) ss
    ss = (a+b+c)/2.0
    s = sqrt(ss*(ss-a)*(ss-b)*(ss-c))
end subroutine

gfortranによるコンパイルは以下のとおりです。

gfortran -shared -fPIC -o heron.dylib heron.f90

実行するPythonは以下です。

from ctypes import * 

heron = cdll.LoadLibrary("heron.dylib")

heron.heron_.argtypes = [POINTER(c_double),POINTER(c_double),POINTER(c_double),POINTER(c_double)]
heron.heron_.restype = c_void_p  

a = 3.0
b = 4.0
c = 5.0
s = 0.0
a = c_double(a)
b = c_double(b)
c = c_double(c)
s = c_double(s)
 
heron.heron_(byref(a),byref(b),byref(c),byref(s)) 
print(s.value)

戻り値の引数、s は Python 側で予め宣言しておかないと「not defined」になるので、とりあえず 0.0 を代入しています。(Fortran側でFunctionにすれば良いのですが、gfortanでは上手くいかず)

上記はアーキテクチャやFortranコンパイラによって書き方が異なるのでご注意ください。

python-docxのtable内にtableを配置する

Microsoft Wordの文章で、表に中のセルに更に表を埋め込む場合があると思います。

from docx import Document
from docx.shared import Cm
from docx.enum.text import WD_ALIGN_PARAGRAPH
from docx.enum.table import WD_TABLE_ALIGNMENT
from docx.enum.section import WD_ORIENT

doc = Document()
section = doc.sections[-1]
new_width, new_height = section.page_height, section.page_width
section.orientation = WD_ORIENT.LANDSCAPE
section.page_width = new_width
section.page_height = new_height

table = doc.add_table(3, 3, style='Table Grid')
j = 0
for row in table.rows:
    for cell in row.cells:
        if j == 4:
            tbl_in_tbl = cell.add_table(3, 3)
            k = 0
            for mini_row in tbl_in_tbl.rows:
                for mini_cell in mini_row.cells:
                    p = mini_cell.add_paragraph()
                    p.alignment=WD_ALIGN_PARAGRAPH.CENTER
                    r = p.add_run()
                    s = str(k+11) + ".jpg"
                    r.add_picture(s, width=Cm(1.5), height=Cm(1))
                    k += 1
        else:
            p = cell.add_paragraph()
            p.alignment=WD_ALIGN_PARAGRAPH.CENTER
            r = p.add_run()
            s = str(j+1) + ".jpg"
            r.add_picture(s, width=Cm(5), height=Cm(4))
        j += 1

doc.save("test.docx")

手元の環境では、セル内に表示したtableに対するstyleは設定できませんでした。

FEMAP API で行う処理の性能を向上する

brown hourglass on brown wooden table
Photo by Mike on Pexels.com

FEMAPのAPIを使用して処理を行う際、対象とする要素や節点の数が多いとかなり時間がかかることがあります。

FEMAPはユーザインタフェースが表示されている間は、モデルの状態を反映するために絶えずアップデートされています。

処理実行中にアップデートの必要がなければ feAppVisible(False) でユーザインタフェース自体を非表示にしてしまうか、 feAppLock でロックしてしまうかのどちらかで処理速度を改善できます。

処理の内容により改善の度合いは異なりますが、多くの場合10倍〜20倍程度の改善が見込めます。

処理完了後は feAppVisible(True) や feAppUnlock の呼び出しをお忘れなく。

pythonで Microsoft Word 書類を grep する

docx2txtはwordファイルからテキストを抽出するライブラリです。

以下はwordファイルから特定の文字を抽出するサンプルです。

import docx2txt

strings = docx2txt.process("000497699.docx")
lines = strings.split("\n")

moji = "厚生"

for line in lines:
    if line.find(moji) >= 0:
        print(line)

対象のwordファイルは厚生労働省のサイトにある、総合研究報告書というファイルを使っています。

以下、出力結果です。

厚生労働科学研究費
 厚生労働行政推進調査事業費
 厚生労働科学研究費
 厚生労働行政推進調査事業費様式A(10)
    厚生労働大臣
 厚生労働科学研究費
 厚生労働行政推進調査事業費
 厚生労働科学研究費
 厚生労働行政推進調査事業費
   上記補助事業について、厚生労働科学研究費補助金等取扱規程(平成10年4月9日厚生省告示第130号)第16条第3項の規定に基づき下記のとおり研究成果を報告します。
  ※規程19条第2項及び第3項に従い、事業完了後5年以内に、その結果又は経過の全部若しくは一部を刊行し、又は書籍、雑誌、新聞等に掲載した場合には、その刊行物又はその別刷一部を添えて厚生労働大臣等に届け出ること。
 研究代表者 厚生 太郎
       厚生太郎
        研究代表者 厚生 太郎 ○○○○○病院長
     厚生労働行政の課題との関連性を含めて記入すること。                                   
           なお、ヒトゲノム・遺伝子解析研究に関する倫理指針(平成25年文部科学省・厚生労働省・経済産業省告示第1号)、人を対象とする医学系研究に関する倫理指針(平成26年文部科学省・厚生労働省告示第3号)、遺伝子治療等臨床研究に関する指針(平成31年厚生労働省告示第48号)、厚生労働省の所管する実施機関における動物実験等の実施に関する基本指針(平成18年6月1日付厚生労働省大臣官房厚生科学課長通知)及び申請者が所属する研究機関で定めた倫理規定等を遵守するとともに、あらかじめ当該研究機関の長等の承認、届出、確認等が必要な研究については、研究開始前に所定の手続を行うこと。

Femap API Results Browsing Object のサンプル

簡単なサンプルです。

◆特定のアウトプットベクトルの値を取得する
アウトプットセット1、アウトプットベクトル8006の値を取得する

import pythoncom
import Pyfemap as Pyfemap
from Pyfemap import constants

existObj = pythoncom.connect(Pyfemap.model.CLSID)
app = Pyfemap.model(existObj)

ID_list = []
entity_set = app.feSet
RBO = app.feResults
RBO.VectorEntitiesV2(1, 8006, entity_set.ID, True)

entity_count = entity_set.Count()
_ = entity_set.First()
count = 0
while count < entity_count:
    ID_list.append(entity_set.CurrentID)
    _ = entity_set.Next()
    count += 1

for ID in ID_list:
    _, f = RBO.EntityValueV2(1, 8006, ID)
    print(ID, f)

◆変位ベクトルを格納する
アウトプットセット1、節点1と2に変位ベクトル(9,000,001、9,000,002、9,000,003、9,000,004)を格納する

import pythoncom
import Pyfemap as Pyfemap
from Pyfemap import constants

existObj = pythoncom.connect(Pyfemap.model.CLSID)
app = Pyfemap.model(existObj)

node = []
disp_x = []
disp_y = []
disp_z = []

node.append(1)
node.append(2)

disp_x.append(1.0)
disp_x.append(0.5)

disp_y.append(1.0)
disp_y.append(2.0)

disp_z.append(0.5)
disp_z.append(3.0)

RBO = app.feResults

nCol = RBO.AddVectorAtNodeColumnsV2(1, 9000001, 9000002, 9000003, 9000004, "USER Displacement", 1, True)
_ = RBO.SetVectorAtNodeColumnsV2(nCol[1], len(node), node, disp_x, disp_y, disp_z)
RBO.Save()

◆BEAM要素のA端、B端のアウトプットベクトルを格納する
要素番号1,2のBEAM要素に端部の応力を格納する

import pythoncom
import Pyfemap as Pyfemap
from Pyfemap import constants

existObj = pythoncom.connect(Pyfemap.model.CLSID)
app = Pyfemap.model(existObj)

elem = []
stressA = []
stressB = []

elem.append(1)
elem.append(2)

stressA.append(-10.0)
stressB.append(10.0)
stressA.append(-20.0)
stressB.append(20.0)

RBO = app.feResults

_, nCol = RBO.AddScalarAtBeamColumnsV2(1, 9000101, 9000102, "USER Stress", 2, 4, True, False)
print(nCol)
rc = RBO.SetColumn(nCol[0], len(elem), elem, stressA)
rc = RBO.SetColumn(nCol[1], len(elem), elem, stressB)
RBO.Save()