[MATLAB] 自作カラーマップと数字状に穴の開いた有限要素法解析モデルを利用して2022年の始まりを祝う

# 挨拶 & 雑談
2022年, あけましておめでとうございます (記事を書き始めた時は正月休みでした). 本年もMATLAB Student Ambassadorをよろしくお願いいたします. さて, [Twitter](https://twitter.com/hibs_MATLAB_Amb)における, 2022年最初の投稿はこんな感じでした.

<blockquote class="twitter-tweet"><p lang="ja" dir="ltr">Happy New Year!!<br>新年, 明けましておめでとうございます.<br>本年もよろしくお願いいたします.<br><br>2022年の干支にちなんで, 年末に行った有限要素法によるトラ柄の研究の成果です.<a href="https://twitter.com/hashtag/MATLAB?src=hash&amp;ref_src=twsrc%5Etfw">#MATLAB</a> <a href="https://t.co/d18CJMWBlq">pic.twitter.com/d18CJMWBlq</a></p>&mdash; ヒビヤ@MATLAB Ambassador (@hibs_MATLAB_Amb) <a href="https://twitter.com/hibs_MATLAB_Amb/status/1476931411990900750?ref_src=twsrc%5Etfw">December 31, 2021</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>

これは私のスウェーデン留学時代からの超個人的な習慣です. スウェーデンに留学していたときに構造力学の授業があり, 友達と「授業で学んだことを生かして何か作りたいね～」なんて話をしてました. その結果, こんなものが完成しました. 当時は「日本では正月になると新春とか言われるんだよー」と話しながらcolormapをspringにしてましたが, MATLAB Student Ambassadorになった2022年には, もう少し頑張ろうという事で, 2022年の干支であるトラを頑張ってカラーマップで表現してみました!!

この記事では, 使用したカラーマップに加えて, 数字を模した穴のある有限要素法解析モデルをMATALBで作るコードを紹介したいと思います!!

# 利用したバージョン及び製品
本記事に掲載されているコードは**R2021b**において検証を行っています.
また, 本記事に掲載されているコードをMATLABで実行するためには, [**Partial Differential Equation Toolbox**](https://jp.mathworks.com/products/pde.html)が必要になります.

# MATLAB における colormap の作成
MATLABでは自分で[colormap](https://jp.mathworks.com/help/matlab/ref/colormap.html)を作成することができます. 方法は簡単で, 0.0~1.0の範囲の値から成る3列の行列を作成し, 1列目に赤, 2列目に緑, 3列目に緑の強度を入力します ([MATLAB Documentation](https://jp.mathworks.com/help/matlab/ref/colormap.html)より). 最大値が255ではないところは要注意です!! 今回は, トラを表現するために, 以下のcolormapを作成しました.

```math
\begin{align}
red(x) &= \left\{
\begin{array}{ll}
0 & (x \lt 0.2) \\
(5x-1)^4 & (0.2 \leq x < 0.4) \\
1 & (0.4 \leq x < 0.6) \\
(4-5x)^4 & (0.6 \leq x < 0.8) \\
0 & (x > 0.8)
\end{array}\right. \\
green(x) &= \left\{
\begin{array}{ll}
0 & (x \lt 0.2) \\
0.8431(5x-1)^4 & (0.2 \leq x < 0.4) \\
0.8431 & (0.4 \leq x < 0.6) \\
0.8431(4-5x)^4 & (0.6 \leq x < 0.8) \\
0 & (x > 0.8)
\end{array}\right. \\
blue(x) &= 0
\end{align}
```

```code:matlab 
n = 1000;
cmap = zeros(n,3);
color1 = [0,0,0];
color2 = [255/255,215/255,0/255];

grad_start1 = 0.2;
grad_finish1 = 0.4;
grad_start2 = 0.6;
grad_finish2 = 0.8;
order = 4;

for ii=1:n
    if (ii/n < grad_start1)
        cmap(ii,:) = color1;
    elseif (ii/n < grad_finish1)
        cmap(ii,:) = color2*(1-((grad_finish1*n-ii)/n/(grad_finish1-grad_start1))^order);
    elseif (ii/n < grad_start2)
        cmap(ii,:) = color2;
    elseif (ii/n < grad_finish2)
        cmap(ii,:) = color2*((grad_finish2*n-ii)/n/(grad_finish2-grad_start2))^order;
    else
        cmap(ii,:) = color1;
    end
end
cmap(cmap > 1) = 1;
```

![tiger_color_map.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1635383/832f9170-14f0-0657-18bf-2aacb2e5d778.png)

# 穴あき有限要素法モデルの作成
次に, 数字を模した穴をあけた解析モデルを作成しました. これには、本記事の最後に載せた自作関数, crate_NumModel, を利用しています. この関数では, 穴に表示する1つの整数Numと穴の大きさを指定するパラメータradius_rateを受け取ってから, [2-D Geometry Creation at Command Line](https://jp.mathworks.com/help/pde/ug/create-geometry-at-the-command-line.html)のルールに基づいた配列を作成し, [decsg関数](https://jp.mathworks.com/help/pde/ug/decsg.html)を使ってモデルを作成しています.

今回は, 2022年を祝うための解析モデルが必要なため, 以下を作成しました.

![HappyNewYear2022.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1635383/fa8ae161-c1c2-db6c-dcb9-8155577e9eb9.png)

今のところ桁数に上限は設けいていないため, より大きな数字も出力が可能です.

![matlab_unicode.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1635383/d951d663-53e6-5af6-8c67-c5a30387acde.png)

# Partial Differential Equation Toolboxを利用した有限要素法解析

この部分については, 別の記事([[Blender+MATLAB] 有限要素法解析を行い, グラデーションのついた文字を作る!!](https://qiita.com/hibs_MATLAB_Amb/items/04009f25b98a18017bba))に書かれておりますので, ぜひご覧ください!!

基本的に解析において必要なのは,

1. 材料定数(ヤング率, ポアソン比)の決定 ([structuralProperties](https://jp.mathworks.com/help/pde/ug/pde.structuralmodel.structuralproperties.html))
1. 負荷する荷重の大きさと方向の決定 ([structuralBoundaryLoad](https://jp.mathworks.com/help/pde/ug/pde.structuralmodel.structuralboundaryload.html))
1. 固定する部分の決定 ([structuralBC](https://jp.mathworks.com/help/pde/ug/pde.structuralmodel.structuralbc.html))
1. メッシュの作成 ([generateMesh](https://jp.mathworks.com/help/pde/ug/generate-mesh.html))

となります. 今回は以下のように設定の上で解析を行いました.

![HappyNewYear2022_Model.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1635383/7bcefbc2-75fe-d617-ddec-abed77371c93.png)

![HappyNewYear2022_Mesh.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1635383/8799f4bc-041b-b253-598a-a78e89695bea.png)


```code:matlab
% Material
E = 197*10^9; % ヤング率 [Pa]
nu = 0.3; % ポアソン比
structuralProperties(model,'YoungsModulus',E,'PoissonsRatio',nu);

% Boundary Condition
structuralBoundaryLoad(model,'Edge',3,'SurfaceTraction',[100;0]);
structuralBoundaryLoad(model,'Edge',4,'SurfaceTraction',[0;100]);
structuralBC(model,'Edge',1,'XDisplacement',0);
structuralBC(model,'Edge',2,'YDisplacement',0);

% Generate Mesh
mesh = generateMesh(model,'Hmax',1/70,'Hmin',1/70,'Hgrad',2.0,'GeometricOrder','linear');

% Analysis
R = solve(model);
```

# 結果
このように, カラーマップ及び解析モデルを作成して有限要素法解析を行った結果, 見事, トラのような柄(???)の解析結果を得られました!!

```code:matlab
figure
pdeplot(model,'XYData',R.Stress.xx,'ColorMap',cmap,'ColorBar','off')
c = colorbar('southoutside');
c.Label.String = 'Stress [MPa]';
axis equal
title 'Normal Stress Along x-Direction';
set(gca,'FontSize',14);

figure
pdeplot(model,'XYData',R.Strain.yy,'ColorMap',cmap,'ColorBar','off')
colorbar('southoutside')
axis equal
title 'Normal Strain Along y-Direction';
set(gca,'FontSize',14);
```

![HappyNewYear2022_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1635383/9062f93f-99b9-b16a-3eb3-a649650a09e2.png)

![HappyNewYear2022_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1635383/b22d7a68-518f-6030-8355-ab3e4e1eff65.png)


# 穴あきの解析モデルを作成する自作関数
**R2021b**でのみ検証を行っております.

```code:matlab
function model = create_NumModel(Num,radius_rate)

%
% create_NumModel - Create FEM Model including Number 
% 
%   model = create_NumModel(num,Radius_Rate)
%   
%   -Inputs-
%   Num         - An Integer including in the model
%   Radius_Rate - Define the Size of Circle (Default : 1)
%   
%   -Outputs-
%   model   - FEM Model
%
%   --- Usage Example ---
%   
%   model = create_NumModel(2022,1);
%
%   ---------------------
%
%

switch nargin
    case 1
        radius_rate = 1;
end



if (Num ~= round(Num))
    % Judge the Num is Integer or not
    error('Error!! Number should be an Integer.');
elseif (numel(Num) > 1)
    % Judge the Num is Integer or not
    error('Error!! Number should be One Integer.');
else
    % Sort the Input Number
    numchar = char(num2str(round(Num)));
    n = length(numchar);
    num = zeros(1,n);

    for ii =1:n
        if (numchar(ii) == '0')
            num(ii) = 10;
        else
            num(ii) = str2double(numchar(ii));
        end
    end
    
    % Data of relationship between Number and Whole Location
    ND = zeros(10,13);
    ND( 1,:) = [0 0 1 0 1 0 0 1 0 1 0 0 1];
    ND( 2,:) = [1 1 1 0 1 1 1 1 1 0 1 1 1];
    ND( 3,:) = [1 1 1 0 1 1 1 1 0 1 1 1 1];
    ND( 4,:) = [1 0 1 1 1 1 1 1 0 1 0 0 1];
    ND( 5,:) = [1 1 1 1 0 1 1 1 0 1 1 1 1];
    ND( 6,:) = [1 1 1 1 0 1 1 1 1 1 1 1 1];
    ND( 7,:) = [1 1 1 1 1 0 0 1 0 1 0 0 1];
    ND( 8,:) = [1 1 1 1 1 1 1 1 1 1 1 1 1];
    ND( 9,:) = [1 1 1 1 1 1 1 1 0 1 1 1 1];
    ND(10,:) = [1 1 1 1 1 1 0 1 1 1 1 1 1];
    
    % Size Data
    l = 1; % Plane Size
    r = radius_rate*(l/4)/10; % Radius
    w = 15; w2 = w/2; % Define Margin Rate

    % Hole Location
    hole_data = zeros(13,2);
    hole_data(1,:) = [-l/w l/w2];
    hole_data(2,:) = [0 l/w2];
    hole_data(3,:) = [l/w l/w2];
    hole_data(4,:) = [-l/w l/w];
    hole_data(5,:) = [l/w l/w];
    hole_data(6,:) = [-l/w 0];
    hole_data(7,:) = [0 0];
    hole_data(8,:) = [l/w,0];
    hole_data(9,:) = [-l/w -l/w];
    hole_data(10,:) = [l/w -l/w];
    hole_data(11,:) = [-l/w -l/w2];
    hole_data(12,:) = [0 -l/w2];
    hole_data(13,:) = [l/w -l/w2];
    
    % Create Model
    model = createpde('structural','static-planestress');
    R11 = [3,4,0,l*n/4,l*n/4,0,l/4,l/4,-l/4,-l/4]';
    gm = R11;
    sf = 'R11';
    ns = "R11";
    holenum = 0;
    for ii=1:n
        for jj=1:13
            if (ND(num(ii),jj) == 1)
                C1 = zeros(length(R11),1);
                C1(1:4,1) = [1,(2*ii-1)*l/8+hole_data(jj,1),hole_data(jj,2),r]';
                gm = [gm C1];
                holenum = holenum+1;
                str = sprintf('C%d',holenum);
                sf = sprintf('%s-%s',sf,str);
                str = sprintf("C%d",holenum);
                ns = [ns;str];
            end
        end
    end
    %sf = 'R1-C1-C2-C3-C4-C5-C6-C7-C8-C9-C10-C11';
    ns = char(ns);
    ns = ns';
    g = decsg(gm,sf,ns);
    geometryFromEdges(model,g);
end
```

# 最後に
なぜ, 年始を過ぎてまで一生懸命書いたかというと, 来年の干支である卯(ウサギ)は, 白一色というイメージが強いので, カラーマップの作成の必要がなくなってしまうためでした. 最後まで読んでいただきありがとうございました.