clc;
clear all;
close all;

%import data
volu=importdata('C_objem.txt');
elemnum=volu.data(:,1);        %element number
volu=volu.data(:,2);           %element volume [mm3]

s1=importdata('C_s1.txt');   
s1=s1.data(:,2);                   %maximal principal stress [MPa]
s2=importdata('C_s2.txt');
s2=s2.data(:,2);                   %middle principal stress [MPa]
s3=importdata('C_s3.txt');
s3=s3.data(:,2);                   %minimal principal stress [MPa]

elem_node=importdata('C_elem_uzl_uprava.txt'); %elements and its assigned nodes table

%Weibull parameters
m=7.19;
sigma0=473.8;

%calculation of probability of failure
for i=1:length(elemnum)
      V=nnz(elem_node(i,2:end));
    for j=1:8                       
        n=elem_node(i,j+1);
        if n == 0
            BB(i,j)=0;
            
        else
        sigman=@(fi,theta) (s1(n).*(cos(fi)).^2+s2(n).*(sin(fi)).^2).*(sin(theta)).^2+s3(n).*(cos(theta)).^2;
        sigmae=@(fi,theta) (sigman(fi,theta)+abs(sigman(fi,theta)))./2;
        
        bb=@(theta,fi)(2.*m+1).*1./(2.*pi.*sigma0.^m).*sin(theta).*sigmae(fi,theta).^m.*volu(i)./V;
        BB(i,j)=quad2d(bb,0,pi/2,0,2*pi,'AbsTol',1e-10);
        end;
    end;
end;

B=4*sum(sum(BB));
Pf=(1-exp(-B))*100;

disp('Riziko porušenia súčasti B je')
B
disp('Pravdepodobnosť porušenia súčasti Pf v percentách je')
Pf
