- Creates a *.grd file
- Runs elmergrid to generate the mesh files from the *.grd file
- Creates a *.sif file (save some result to textfile, create paraview file)
- Runs elmersolver
- Reads the result from the text file and shows it in Octave
- Creates a command file for elmerpost
- Runs the command file with elmerpost
- Runs Paraview
(Does somebody know how to activate the "eye" in Paraview from the command line to display the plot?)
Octave is a free Matlab "clone" and can be found here for windows: http://octave.sourceforge.net/
Have Fun!
ElmerSimulation.m:
- Code: Select all
%this function generates the input file for a Elmer fem simulation, runs the simulation and does some postprocessing
%define paths ====================================================================================================
elmer_path = 'C:\Elmer6.2\bin';
paraview_path = '"C:\Program Files (x86)\ParaView 3.10.1\bin\paraview.exe"';
project_path = 'C:\elmerfiles';
%define voltage====================================================================================================
Ua = 1;
%create string for line seperator========================================================================================
lb = sprintf('\n');
%generate *.grd file=================================================================================================
%define path
grd_path = [project_path filesep 'geometry.grd'];
%delete old file if it exists
if isequal(exist(grd_path,'file'),2)
try
delete(grd_path);
catch
end
end
%open file
fileid = fopen(grd_path,'w');
%write file
fwrite(fileid,['##### ElmerGrid input file for structured grid generation ######' lb]);
fwrite(fileid,['Version = 210903' lb]);
fwrite(fileid,['Coordinate System = Cartesian 2D' lb]);
fwrite(fileid,['Subcell Divisions in 2D = 5 3' lb]);
fwrite(fileid,['Subcell Sizes 1 = 2000e-6 100e-6 2000e-6 100e-6 2000e-6' lb]); %sizes in x-direction
fwrite(fileid,['Subcell Sizes 2 = 100e-6 100e-6 300e-6' lb]); %sizes in y-direction
fwrite(fileid, lb);
fwrite(fileid,['Material Structure in 2D' lb]); %1: silicon, 2: metal, 0,3,4: dummy surrounding (helps with boundary definitions)
fwrite(fileid,[' 1 1 1 1 1 ' lb]);
fwrite(fileid,[' 0 2 0 3 0' lb]);
fwrite(fileid,[' 0 4 0 4 0' lb]);
fwrite(fileid,['End' lb]);
fwrite(fileid, lb);
fwrite(fileid,['Materials Interval = 1 3' lb]);
fwrite(fileid, lb);
fwrite(fileid,['Boundary Definitions' lb]);
%Each boundary requires a type that refers to the number of the boundary condition. It may then later
%be referred to when setting up the boundary conditions. The boundaries are created in such a way that
%the material int should have elements created. The out may be an empty subcell. If the value of the
%flag double is 2 or more secondary nodes on the boundary will be created. This enables the use of
%discontinuous boundary conditions.
%It is also possible to create boundary conditions on a given side of a material. The number −1 refers
%to down, −2 to right, −3 to top and −4 to left. This means that the boundary is created on the given
%side of the material that is specified by the other index.
%In order to make the creation of the boundary condition more robust there are some special numbers
%that may be used. −11 refers to smaller, −9 to bigger, and −10 to anything. For example, if the
%intmat is 5 and the outmat is −9, then the boundary is created for all materials that have a larger
%material number than 5 and have a common boundary with material 5. This makes it possible to create
%several material boundaries with one definition.
fwrite(fileid,[' #type out int double' lb]);
fwrite(fileid,[' 1 4 2 ' lb]); %left exterior metal boundary
fwrite(fileid,[' 2 2 1 ' lb]); %left inner metal boundary
fwrite(fileid,[' 3 4 3 ' lb]); %right exterior metal boundary
fwrite(fileid,[' 4 3 1 ' lb]); %right inner metal boundary
fwrite(fileid,[' 5 0 1 ' lb]); %other boundaries bulk
fwrite(fileid,[' 6 0 2 ' lb]); %other boundaries metal left
fwrite(fileid,[' 7 0 2 ' lb]); %other boundaries metal right
fwrite(fileid,['End' lb]);
fwrite(fileid, lb);
fwrite(fileid,['Numbering = Horizontal' lb]);
fwrite(fileid,['Element Degree = 2' lb]);
fwrite(fileid,['Element Innernodes = False' lb]);
fwrite(fileid,['!Triangles = True' lb]);
fwrite(fileid,['Reference Density = 10e-6' lb]);
fwrite(fileid,['Element Ratios 1 = 1 1 1 1 1' lb]);
fwrite(fileid,['Element Ratios 2 = 1 1 1' lb]);
fwrite(fileid,['Element Densities 1 = 1 1 1 1 1' lb]);
fwrite(fileid,['Element Densities 2 = 1 1 1' lb]);
%close file
fclose(fileid);
%generate *.mesh files ==============================================================================================
%delete old files if they exist
if isequal(exist([project_path filesep 'mesh.header'],'file'),2)
try
delete([project_path filesep 'mesh.header']);
delete([project_path filesep 'mesh.nodes']);
delete([project_path filesep 'mesh.elements']);
delete([project_path filesep 'mesh.boundary']);
catch
end
end
mesh_command = [ elmer_path filesep 'elmergrid.exe 1 2 ' grd_path ' -out ' project_path]; %run elmergrid.exe without arguments to get help on the command line parameters
system(mesh_command);
%generate *.sif file =================================================================================================
%define path
sif_path = [project_path filesep 'case.sif'];
%delete old file if it exists
if isequal(exist(sif_path,'file'),2)
try
delete(sif_path);
catch
end
end
%open file
fileid = fopen(sif_path,'w');
%write file
fwrite(fileid,['Header' lb]);
fwrite(fileid,[' CHECK KEYWORDS Warn' lb]);
fwrite(fileid,[' Mesh DB "' project_path '"' lb]);
fwrite(fileid,[' Include Path "' project_path '"' lb]);
fwrite(fileid,[' Results Directory "' project_path '"' lb]);
fwrite(fileid,['End' lb]);
fwrite(fileid, lb);
fwrite(fileid,['Simulation' lb]);
fwrite(fileid,[' Max Output Level = 10' lb]);
fwrite(fileid,[' Coordinate System = Cartesian' lb]);
fwrite(fileid,[' Coordinate Mapping(3) = 1 2 3' lb]);
fwrite(fileid,[' Simulation Type = Steady state' lb]);
fwrite(fileid,[' Steady State Max Iterations = 1' lb]);
fwrite(fileid,[' Output Intervals = 1' lb]);
fwrite(fileid,[' Timestepping Method = BDF' lb]);
fwrite(fileid,[' BDF Order = 1' lb]);
fwrite(fileid,[' Solver Input File = case.sif' lb]);
fwrite(fileid,[' Post File = case.ep' lb]);
fwrite(fileid,['End' lb]);
fwrite(fileid, lb);
fwrite(fileid,['Constants' lb]);
fwrite(fileid,[' Permittivity of Vacuum = 8.8542e-12' lb]);
fwrite(fileid,[' Boltzmann Constant = 1.3807e-23' lb]);
fwrite(fileid,[' Unit Charge = 1.602e-19' lb]);
fwrite(fileid,['End' lb]);
fwrite(fileid, lb);
fwrite(fileid,['Body 1' lb]);
fwrite(fileid,[' Target Bodies(1) = 1' lb]);
fwrite(fileid,[' Name = "bulk"' lb]);
fwrite(fileid,[' Equation = 1' lb]);
fwrite(fileid,[' Material = 1' lb]);
fwrite(fileid,[' Body Force = 1' lb]);
fwrite(fileid,['End' lb]);
fwrite(fileid, lb);
fwrite(fileid,['Body 2' lb]);
fwrite(fileid,[' Target Bodies(1) = 2' lb]);
fwrite(fileid,[' Name = "metal left"' lb]);
fwrite(fileid,[' Equation = 1' lb]);
fwrite(fileid,[' Material = 2' lb]);
fwrite(fileid,[' Body Force = 1' lb]);
fwrite(fileid,['End' lb]);
fwrite(fileid, lb);
fwrite(fileid,['Body 3' lb]);
fwrite(fileid,[' Target Bodies(1) = 3' lb]);
fwrite(fileid,[' Name = "metal right"' lb]);
fwrite(fileid,[' Equation = 1' lb]);
fwrite(fileid,[' Material = 2' lb]);
fwrite(fileid,[' Body Force = 1' lb]);
fwrite(fileid,['End' lb]);
fwrite(fileid, lb);
%solve for potential ---------------------------------------------------------------------
fwrite(fileid,['Solver 1' lb]);
fwrite(fileid,[' Equation = Static Current Conduction' lb]);
fwrite(fileid,['Calculate Loads = True' lb]);
fwrite(fileid,[' Variable = -dofs 1 Potential' lb]);
fwrite(fileid,[' Procedure = "StatCurrentSolve" "StatCurrentSolver"' lb]);
fwrite(fileid,[' Exec Solver = Always' lb]);
fwrite(fileid,[' Stabilize = True' lb]);
fwrite(fileid,[' Bubbles = False' lb]);
fwrite(fileid,[' Lumped Mass Matrix = False' lb]);
fwrite(fileid,[' Optimize Bandwidth = True' lb]);
fwrite(fileid,[' Steady State Convergence Tolerance = 1.0e-5' lb]);
fwrite(fileid,[' Nonlinear System Convergence Tolerance = 1.0e-8' lb]);
fwrite(fileid,[' Nonlinear System Max Iterations = 20' lb]);
fwrite(fileid,[' Nonlinear System Newton After Iterations = 3' lb]);
fwrite(fileid,[' Nonlinear System Newton After Tolerance = 1.0e-3' lb]);
fwrite(fileid,[' Nonlinear System Relaxation Factor = 1' lb]);
fwrite(fileid,[' Linear System Solver = Iterative' lb]);
fwrite(fileid,[' Linear System Iterative Method = BiCGStab' lb]);
fwrite(fileid,[' Linear System Max Iterations = 500' lb]);
fwrite(fileid,[' Linear System Convergence Tolerance = 1.0e-8' lb]);
fwrite(fileid,[' Linear System Preconditioning = ILU0' lb]);
fwrite(fileid,[' Linear System ILUT Tolerance = 1.0e-3' lb]);
fwrite(fileid,[' Linear System Abort Not Converged = False' lb]);
fwrite(fileid,[' Linear System Residual Output = 1' lb]);
fwrite(fileid,[' Linear System Precondition Recompute = 1' lb]);
fwrite(fileid,['End' lb]);
fwrite(fileid, lb);
%save some output files (for paraview ect.)------------------------------------
fwrite(fileid,['Solver 2' lb]);
fwrite(fileid,[' Exec Solver = After Simulation' lb]);
fwrite(fileid,[' Equation = "result output"' lb]);
fwrite(fileid,[' Procedure = "ResultOutputSolve" "ResultOutputSolver"' lb]);
fwrite(fileid,[' Output File Name = "paraview"' lb]);
fwrite(fileid,[' Binary Output = Logical True' lb]);
fwrite(fileid,[' Single Precision = Logical True' lb]);
%fwrite(fileid,[' Gid Format = Logical True' lb]);
%fwrite(fileid,[' Gmsh Format = Logical True' lb]);
%fwrite(fileid,[' Vtk Format = Logical True' lb]);
fwrite(fileid,[' Vtu Format = Logical True' lb]);
%fwrite(fileid,[' Dx Format = Logical True' lb]);
fwrite(fileid,[' Scalar Field 1 = String Potential' lb]);
fwrite(fileid,['End' lb]);
fwrite(fileid, lb);
%calculate flux convective flux
fwrite(fileid,['Solver 3' lb]);
fwrite(fileid,[' Exec Solver = After All' lb]);
fwrite(fileid,[' procedure = File "SaveData" "SaveScalars"' lb]);
fwrite(fileid,[' Filename = "' project_path filesep 'flux_results.dat"' lb]);
fwrite(fileid,[' variable 1 = Potential Loads' lb]);
fwrite(fileid,[' operator 1 = boundary sum' lb]);
% fwrite(fileid,['Save Coordinates(1,2) = 0 0' lb])
fwrite(fileid,['End' lb]);
fwrite(fileid, lb);
%{
%apply some operators on the data and save results---------------------
% max, min, max abs, min abs, mean, variance, deviation
% int mean, int variance
% volume
% boundary sum, boundary dofs, boundary mean, boundary max, boundary min,
% boundary max abs, boundary min abs, area, boundary int, boundary int mean.
% diffusive energy, convective energy, potential energy
% diffusive flux, convective flux, boundary int, boundary int mean, area
% ...
fwrite(fileid,['Solver 2' lb]);
fwrite(fileid,['Equation = SaveScalars' lb]);
fwrite(fileid,['Exec Solver = After All' lb]);
fwrite(fileid,['Procedure = "SaveData" "SaveScalars"' lb]);
%fwrite(fileid,['Save Flux = True' lb]);
fwrite(fileid,['Filename = "' project_path filesep 'operator_results.dat"' lb]);
fwrite(fileid,['Variable 1 = potential loads' lb]);
fwrite(fileid,['Operator 1 = boundary sum' lb]);
%fwrite(fileid,['Variable 2 =coordinate 1' lb]);
%fwrite(fileid,['Save Coordinates(1,2) = 0 0' lb])
fwrite(fileid,['End' lb]);
fwrite(fileid, lb);
%}
%{
%save some line data--------------------------------------------------------------------
fwrite(fileid,['Solver 4' lb]);
fwrite(fileid,['Exec Solver = After Simulation' lb]);
fwrite(fileid,['Equation = SaveLine' lb]);
fwrite(fileid,['Procedure = "SaveData" "SaveLine"' lb]);
fwrite(fileid,['Filename = "' project_path filesep 'line.dat"' lb]);
fwrite(fileid,['Polyline coordinates(2,2) = 0 150e-6 6200e-6 150e-6' lb]);
%fwrite(fileid,['Save Flux = True' lb]);
%fwrite(fileid,['Variable 1 = Potential' lb]);
%fwrite(fileid,['Variable 2 = coordinate 1' lb]);
%fwrite(fileid,['Variable 3 = coordinate 2' lb]);
fwrite(fileid,['End' lb]);
fwrite(fileid, lb);
%}
%equation-------------------------------------------------------------------------------------
fwrite(fileid,['Equation 1' lb]);
fwrite(fileid,[' Name = "Equation 1"' lb]);
fwrite(fileid,[' Active Solvers(1) = 1' lb]);
fwrite(fileid,['End' lb]);
fwrite(fileid, lb);
%bulk properties-----------------------------------------------------------------------
fwrite(fileid,['Material 1' lb]);
fwrite(fileid,[' Name = "silicon"' lb]);
fwrite(fileid,[' Electric Conductivity = 1' lb]);
fwrite(fileid,['End' lb]);
fwrite(fileid, lb);
%metal properties-----------------------------------------------------------------------
fwrite(fileid,['Material 2' lb]);
fwrite(fileid,[' Name = "metal"' lb]);
fwrite(fileid,[' Electric Conductivity = 2' lb]);
fwrite(fileid,['End' lb]);
fwrite(fileid, lb);
%current source
fwrite(fileid,['Body Force 1' lb]);
fwrite(fileid,[' Name = "Current Source"' lb]);
fwrite(fileid,[' Current Source = 0' lb]);
fwrite(fileid,['End' lb]);
fwrite(fileid, lb);
%boundary conditons------------------------------------------------------------------
%ground potential
fwrite(fileid,['Boundary Condition 1' lb]);
fwrite(fileid,[' Name = "ground"' lb]);
fwrite(fileid,[' Target Boundaries(1) = 1' lb]);
fwrite(fileid,[' Save Scalars = Logical True' lb]);
fwrite(fileid,[' Potential = 0' lb]);
fwrite(fileid,['End' lb]);
fwrite(fileid,lb);
%potential
fwrite(fileid,['Boundary Condition 2' lb]);
fwrite(fileid,[' Name = "external potential Ua"' lb]);
fwrite(fileid,[' Target Boundaries(1) = 3' lb]);
fwrite(fileid,[' Potential = ' num2str(Ua) lb]);
fwrite(fileid,['End' lb]);
fwrite(fileid,lb);
%insulation
fwrite(fileid,['Boundary Condition 3' lb]);
fwrite(fileid,[' Name = "insulation"' lb]);
fwrite(fileid,[' Target Boundaries(1) = 5' lb]);
fwrite(fileid,[' Current Density = 0' lb]);
fwrite(fileid,['End' lb]);
fwrite(fileid,lb);
fwrite(fileid,['Boundary Condition 4' lb]);
fwrite(fileid,[' Name = "insulation"' lb]);
fwrite(fileid,[' Target Boundaries(1) = 6' lb]);
fwrite(fileid,[' Current Density = 0' lb]);
fwrite(fileid,['End' lb]);
fwrite(fileid,lb);
fwrite(fileid,['Boundary Condition 5' lb]);
fwrite(fileid,[' Name = "insulation"' lb]);
fwrite(fileid,[' Target Boundaries(1) = 7' lb]);
fwrite(fileid,[' Current Density = 0' lb]);
fwrite(fileid,['End' lb]);
fwrite(fileid,lb);
%close file
fclose(fileid);
%generate start info for elmersolver.exe=================================================================================
%define path
solverinfo_path = [project_path filesep 'ELMERSOLVER_STARTINFO'];
%delete old file if it exists
if isequal(exist(solverinfo_path,'file'),2)
try
delete(solverinfo_path);
catch
end
end
%open file
fileid = fopen(solverinfo_path,'w');
%write file
fwrite(fileid,['case.sif' lb]);
fwrite(fileid,['1' lb])
fwrite(fileid, lb)
%close file
fclose(fileid);
%run simulation with elmersolver.exe====================================================================================
%delete old files if they exist
if isequal(exist([project_path filesep 'case.ep'],'file'),2)
try
delete([project_path filesep 'case.ep']);
catch
end
end
simulation_command = [ elmer_path filesep 'elmersolver.exe ' project_path filesep 'case.sif'];
system(simulation_command);
%read results from flux_results.dat=======================================================================================
fileid = fopen([project_path filesep 'flux_results.dat'],'r');
result_string = fgetl(fileid);
fclose(fileid);
%calculate resistance and display results==================================================================================
I = -eval(result_string);
U = Ua;
R = U / I;
%gmsgbox(['R = ' num2str(R,3) ' Ohm, I = ' num2str(I,3) ' A, Ua = ' num2str(U,3) ' V']);
['R = ' num2str(R,3) ' Ohm, I = ' num2str(I,3) ' A, Ua = ' num2str(U,3) ' V']
%generate *.txt command file for elmerpost.exe=================================================================================
%define path
post_commands_path = [project_path filesep 'post_commands.txt'];
%delete old file if it exists
if isequal(exist(post_commands_path,'file'),2)
try
delete(post_commands_path);
catch
end
end
%open file
fileid = fopen(post_commands_path,'w');
%write file
fwrite(fileid,['readfile case.ep' lb]);
fwrite(fileid, lb);
fwrite(fileid,['background 255 255 255' lb]);
fwrite(fileid, lb);
fwrite(fileid,['set DisplayStyle(ColorMesh) 1' lb]);
fwrite(fileid,['set MeshStyle 2' lb]);
fwrite(fileid,['set MeshColor "Potential"' lb]);
fwrite(fileid, lb);
fwrite(fileid,['set DisplayStyle(ColorScale) 1' lb]);
fwrite(fileid,['set ColorScaleFontColor 0' lb]);
fwrite(fileid, lb);
fwrite(fileid,['UpdateObject' lb]);
%close file
fclose(fileid);
%execute post_command.txt with elmerpost.exe===============================================================================
%get current path
original_path = pwd;
%change to project path
eval(['cd ' project_path]);
%run elmerpost.exe
system('elmerpost.exe source post_commands.txt');
%go back to original path
eval(['cd ' original_path]);
%open result file with paraview.exe (you have to click on the eye to activate the view)======================================================
paraview_command = [ paraview_path ' ' project_path filesep 'paraview0001.vtu']
system(paraview_command);
geometry.grd:
- Code: Select all
##### ElmerGrid input file for structured grid generation ######
Version = 210903
Coordinate System = Cartesian 2D
Subcell Divisions in 2D = 5 3
Subcell Sizes 1 = 2000e-6 100e-6 2000e-6 100e-6 2000e-6
Subcell Sizes 2 = 100e-6 100e-6 300e-6
Material Structure in 2D
1 1 1 1 1
0 2 0 3 0
0 4 0 4 0
End
Materials Interval = 1 3
Boundary Definitions
#type out int double
1 4 2
2 2 1
3 4 3
4 3 1
5 0 1
6 0 2
7 0 2
End
Numbering = Horizontal
Element Degree = 2
Element Innernodes = False
!Triangles = True
Reference Density = 10e-6
Element Ratios 1 = 1 1 1 1 1
Element Ratios 2 = 1 1 1
Element Densities 1 = 1 1 1 1 1
Element Densities 2 = 1 1 1
case.sif
- Code: Select all
Header
CHECK KEYWORDS Warn
Mesh DB "D:\SpiceGUI\elmerfiles"
Include Path "D:\SpiceGUI\elmerfiles"
Results Directory "D:\SpiceGUI\elmerfiles"
End
Simulation
Max Output Level = 10
Coordinate System = Cartesian
Coordinate Mapping(3) = 1 2 3
Simulation Type = Steady state
Steady State Max Iterations = 1
Output Intervals = 1
Timestepping Method = BDF
BDF Order = 1
Solver Input File = case.sif
Post File = case.ep
End
Constants
Permittivity of Vacuum = 8.8542e-12
Boltzmann Constant = 1.3807e-23
Unit Charge = 1.602e-19
End
Body 1
Target Bodies(1) = 1
Name = "bulk"
Equation = 1
Material = 1
Body Force = 1
End
Body 2
Target Bodies(1) = 2
Name = "metal left"
Equation = 1
Material = 2
Body Force = 1
End
Body 3
Target Bodies(1) = 3
Name = "metal right"
Equation = 1
Material = 2
Body Force = 1
End
Solver 1
Equation = Static Current Conduction
Calculate Loads = True
Variable = -dofs 1 Potential
Procedure = "StatCurrentSolve" "StatCurrentSolver"
Exec Solver = Always
Stabilize = True
Bubbles = False
Lumped Mass Matrix = False
Optimize Bandwidth = True
Steady State Convergence Tolerance = 1.0e-5
Nonlinear System Convergence Tolerance = 1.0e-8
Nonlinear System Max Iterations = 20
Nonlinear System Newton After Iterations = 3
Nonlinear System Newton After Tolerance = 1.0e-3
Nonlinear System Relaxation Factor = 1
Linear System Solver = Iterative
Linear System Iterative Method = BiCGStab
Linear System Max Iterations = 500
Linear System Convergence Tolerance = 1.0e-8
Linear System Preconditioning = ILU0
Linear System ILUT Tolerance = 1.0e-3
Linear System Abort Not Converged = False
Linear System Residual Output = 1
Linear System Precondition Recompute = 1
End
Solver 2
Exec Solver = After Simulation
Equation = "result output"
Procedure = "ResultOutputSolve" "ResultOutputSolver"
Output File Name = "paraview"
Binary Output = Logical True
Single Precision = Logical True
Vtu Format = Logical True
Scalar Field 1 = String Potential
End
Solver 3
Exec Solver = After All
procedure = File "SaveData" "SaveScalars"
Filename = "D:\SpiceGUI\elmerfiles\flux_results.dat"
variable 1 = Potential Loads
operator 1 = boundary sum
End
Equation 1
Name = "Equation 1"
Active Solvers(1) = 1
End
Material 1
Name = "silicon"
Electric Conductivity = 1
End
Material 2
Name = "metal"
Electric Conductivity = 2
End
Body Force 1
Name = "Current Source"
Current Source = 0
End
Boundary Condition 1
Name = "ground"
Target Boundaries(1) = 1
Save Scalars = Logical True
Potential = 0
End
Boundary Condition 2
Name = "external potential Ua"
Target Boundaries(1) = 3
Potential = 1
End
Boundary Condition 3
Name = "insulation"
Target Boundaries(1) = 5
Current Density = 0
End
Boundary Condition 4
Name = "insulation"
Target Boundaries(1) = 6
Current Density = 0
End
Boundary Condition 5
Name = "insulation"
Target Boundaries(1) = 7
Current Density = 0
End
ELMERSOLVER_STARTINFO:
- Code: Select all
case.sif
1
post_commands.txt:
- Code: Select all
readfile case.ep
background 255 255 255
set DisplayStyle(ColorMesh) 1
set MeshStyle 2
set MeshColor "Potential"
set DisplayStyle(ColorScale) 1
set ColorScaleFontColor 0
UpdateObject
