Fractals in matlab/octave

Codes in matlab/octave for generating some popular fractals:

Mandelbrot

mandelbrot

function mandelbrotmain()
clear variables; clc;

maxit=200;
x=[-2,1];
y=[-1 1];

xpix=601;
ypix=401;

x=linspace(x(1),x(2),xpix);
y=linspace(y(1),y(2),ypix);
[xG, yG]=meshgrid(x,y);

c=mb(maxit,xG,yG);

figure
imagesc(x,y,c);
colormap([1 1 1;0 0 0]);
axis on;
grind on;

endfunction

function count=mb(maxItr,xG,yG)
c=xG+1i*yG;
count=ones(size(c));
z=c;

for n=1:maxItr
z=z.*z+c;
inside=abs(z)<=2;
count=count+inside;
endfor
endfunction

 

Sierpinski triangle
sierpinski

clc;
clear variables;
close windows;

clf

N=500;
x=zeros(1,N);y=x;r=x;
for a=2:N
c=randi([0 2]);
r(a)=c;
switch c
case 0
x(a)=0.5*x(a-1);
y(a)=0.5*y(a-1);
case 1
x(a)=0.5*x(a-1)+.25;
y(a)=0.5*y(a-1)+sqrt(3)/4;
case 2
x(a)=0.5*x(a-1)+.5;
y(a)=0.5*y(a-1);
end
end
plot(x,y,’o’)
title(‘Sierpinski’s triangle’)
legend(sprintf(‘N=%d Iterations’,N))

 

Fractal tree
fractaltrree

function treemain
tic
clear all;
depth = 7;
figure 1;
hold on;

drawTree2(0, 0, 90, depth);
toc;
endfunction

function drawTree2(x1, y1, angle, depth)
deg_to_rad = pi / 180.0;
rot=30; #degrees :: rotation from second iteration
branchLength=100*depth;
if (depth != 0)
x2 = x1 + cos(angle * deg_to_rad) * branchLength;
y2 = y1 + sin(angle * deg_to_rad) * branchLength ;
line([x1, x2], [y1, y2],’LineWidth’,depth);
drawTree2(x2, y2, angle – rot, depth – 1);
drawTree2(x2, y2, angle + rot, depth – 1);

endif
endfunction

 

A note on fractal (from HH Hardey):
There are upper and lower limits to describe natural objects by fractals. A figure may look like a fern leaf, but continued generation of the “leaf on larger and larger scales produces only a larger and larger fern leaf. A fern “plant” will never be produced. Similarly, there is a smallest scale to which a fractal description applies. As the overall fern leaf pattern is repeated on smaller and smaller scales, eventually the scale of a single fern plant cell will be reached. A single cell does not look like the overall shape.

Using Math Kernal Library (MKL) in visual studio with Fortran

First install MKL from intel website

To use intel MKL library inside visual studio do the followings:
1. Right click the project>Properties
2. In Property page>Fortran>Libraries set (a) Runtime library : Multithreaded and (b) Use Intel MKL: Sequential
mklsetting

Once this setting is done, you can compile the file in DEBUG mode. Repeat the same process (1-2) for RELEASE mode too once more.(Note exe of DEBUG mode is far slower than exe of RELEASE mode. The settings to set are listed below.

To do the same thing MANUALLY (which gives more insight to the code) do the following in each DEBUG and RELEASE mode.
1. Property page>Fortran>General>Additional Include Directories: C:\Program Files (x86)\IntelSWTools\compilers_and_libraries_2019.4.245\windows\mkl\include
2. Property page>Fortran>Preprocessor>Preprocess source file: Yes
3. Property page>Linker>General>Additional Library Directories:C:\Program Files (x86)\IntelSWTools\compilers_and_libraries_2019.4.245\windows\mkl\lib\ia32_win
4. Property page>Linker>Input>Additional Dependencies: mkl_intel_c.lib mkl_intel_thread.lib mkl_core.lib libiomp5md.lib

The manual setting is described in the video as well:

For method to install Fortran in visual studio see this post.

Specific gravity of soil particles

Consider the following figure.
specific gravity
W1= weight of empty vessel
W2= Weight of vessel + weight of sample
W3=Weight of vessel + weight of sample+Water
W4=Weight of vessel + Water

Here,
Weight of solid = W2-W1

Similarly, from last two figures
W3-W4=Weight of solid- Weight of water
or, W3-W4= W2-W1 -Weight of water
or, Weight of water=(W2-W1)-(W3-W4)(= weight of equal volume of water as that of solid)

Now,
specific-gravity=\frac{weight of solid}{weight of equal volume of water}
i.e. \rho_s=\frac{(W2-W1)}{(W2-W1)-(W3-W4)}

Field variables in code-aster

Available variables in code aster that I tested. I will be updating this table after testing other elements.

FIELD NAME AVAILABLE VARIABLES APPLICABLE MODEL TYPES NON- APPLICABLE MODEL
EFGE_ELNO N,VY,VZ,MT, MFY,MFZ POU_D_E, 3D, 2D
SIEF_ELNO POU_D_T, POU_D_E
POU_D_TG, POU_D_T_GD
FORC_F FX FY FZ MX MY MZ
DEPL_R DX DY DZ DRX DRY MARTINI DRZ
SIEF_R SIXX SIYY SIZZ SIXY SIXZ SIYZ 3D
SIEF_R NR VY VZ MT MFY MFZ Discrete elements
SIEF_R NXX NYY NXY MXX MYY MXY QX QY HULL
SIEF_R BR Beam with wrapping
SIEF_R POU_D_T, POU_D_E
POU_D_TG, POU_D_T_GD
SIEF_ELGA 3D,2D, DKT, DST, Q4G, Q4GG, COQUE_3D, BARR, 2D_BARR POU_D_T, POU_D_E
POU_D_TG, POU_D_T_GD
SIGM_ELGA 3D,2D, DKT, DST, Q4G, Q4GG, COQUE_3D POU_D_T, POU_D_E
POU_D_TG, POU_D_T_GD
SIGM_ELNO POU_D_T, POU_D_E
POU_D_TG, POU_D_T_GD
EFGE_ELNO

Effort

3D, 2D
EFGE_ELNO 3D, 2D
FORC_NODA

nodal forces

DX DY DZ DRX DRY MARTINI DRZ
EPSI_ELGA
EPSG_ELGA 3D Hulls, plates
(except DKTG
and Q4GG), PIPES, Multifibre beams
EPME_ELGA 3D Hulls, plates
(except DKTG
and Q4GG)
EPMG_ELGA 3D Hulls, plates
(except DKTG
and Q4GG), PIPES, Multifibre beams
EPSP_ELGA 3D Hulls, plates
(except DKTG
and Q4GG)
EPSI_ELNO

 

3D
EPSG_ELNO 3D Hulls, plates
(except DKTG
and Q4GG), PIPES, Multifibre beams
EPME_ELNO 3D Hulls, plates
(except DKTG
and Q4GG)
EPMG_ELNO 3D Hulls, plates
(except DKTG
and Q4GG), PIPES, Multifibre beams
EPSP_ELNO 3D Hulls, plates
(except DKTG
and Q4GG)
DEGE_ELGA

 

Plates,
Coques1D
3D
DEGE_ELNO

 

Plates,
Coques1D
3D
VARI_ELGA BARR, 2D_BARR

Some clues:

DEGE : deformations generalized on the elements of beam or hull
ELNO: Node- calculation carried out by extrapolation with the nodes of the quantities at the points of Gauss
ELGA:Gauss point
NOEU:

EF~:Effort

SI~: Stress

Beginning with Code Aster and Salome Meca

Code aster is an opensource (free) FEM analysis software targetted to structural, thermal and aucostic analysis. It is available for windows and linux. This software can be used for non-linear analysis and thus is quite useful for engineering studies epscially for academic research. But the sofware is little difficult to use for begineers. I am sharing some introductory videos to get started with this software. You can visit the forum for discussion on your specific problem.

  • Analysis of a simple beam

  • Non-linear analysis of simple steel rod

  • Modeling hardeing of steel

  • Analysis of plate and plotting stress strain diagram