I'm new to Matlab programming. I have two models created by tetrahedral meshes by using EIDORS and Netgen. Then i need to create a map based on the intersection between the tetrahedral elements by using Matlab. So, in order to find the intersection i tried using the point test method. Please refer to the link. http://steve.hollasch.net/cgindex/geometry/ptintet.html
Let's assume model1 and model2.
First, i extract the vertices from both models and created matrices.
Then, i used the point test method to calculate the intersection between the models.
For model1 i have 3000++ tetrahedral elements and model2 i have 8000++ tetrahedral elements.
In order to calculate the intersection, i looped one by one to determine which of the elements between the two models have intersection and then i created a matrix by indexing the element number of model2 to the element number of model1.
But, there appears to be zero for few elements which is not possible because all elements from model1 should atleast intersects with a few of the elements from model2.
Therefore, in the end i expect to get a matrix consists of (the element number of model1 x the element number of model2 which intersects with the corresponding elements from model1). Can help me solve this? Please refer to my code.
function [elemno,deter0,deter1,deter2,deter3,deter4] = checkp(filename1,filename2);
%/* to check whether the vertices of a layered model''s element are inside the
% tetrahedron of the generic model
%after the model is created using netgen and eidors, there will be a struct type
%named model_parameters.
%model_parameters.vtx refers to the vertices which consists of (nodes x vertices %(x,y,z))
%model_parameters.simp refers to the elements which consists (numberofelements x nodes) nodes are linked to the vertices.*/
filename1 = [filename1 '.mat'];
filename2 =[filename2 '.mat'];
first = load(filename1);
second =load(filename2);
vtx = first.model_parameters.vtx;
simp = first.model_parameters.simp;
[simpr,simpc] = size(simp);
vtx2 = second.model_parameters.vtx;
simp2= second.model_parameters.simp;
[simpr2,simpc2] = size(simp2);
%//extracting the vertices of the elements from the simplices(element)
for loop1 = 1 : simpr
elemx(loop1,1) = vtx(simp(loop1,1),1);
elemx(loop1,2) = vtx(simp(loop1,2),1);
elemx(loop1,3) = vtx(simp(loop1,3),1);
elemx(loop1,4) = vtx(simp(loop1,4),1);
elemy(loop1,1) = vtx(simp(loop1,1),2);
elemy(loop1,2) = vtx(simp(loop1,2),2);
elemy(loop1,3) = vtx(simp(loop1,3),2);
elemy(loop1,4) = vtx(simp(loop1,4),2);
elemz(loop1,1) = vtx(simp(loop1,1),3);
elemz(loop1,2) = vtx(simp(loop1,2),3);
elemz(loop1,3) = vtx(simp(loop1,3),3);
elemz(loop1,4) = vtx(simp(loop1,4),3);
end
for loop2 = 1:simpr2
elemx2(loop2,1) = vtx2(simp2(loop2,1),1);
elemx2(loop2,2) = vtx2(simp2(loop2,2),1);
elemx2(loop2,3) = vtx2(simp2(loop2,3),1);
elemx2(loop2,4) = vtx2(simp2(loop2,4),1);
elemy2(loop2,1) = vtx2(simp2(loop2,1),2);
elemy2(loop2,2) = vtx2(simp2(loop2,2),2);
elemy2(loop2,3) = vtx2(simp2(loop2,3),2);
elemy2(loop2,4) = vtx2(simp2(loop2,4),2);
elemz2(loop2,1) = vtx2(simp2(loop2,1),3);
elemz2(loop2,2) = vtx2(simp2(loop2,2),3);
elemz2(loop2,3) = vtx2(simp2(loop2,3),3);
elemz2(loop2,4) = vtx2(simp2(loop2,4),3);
end
%//point test calculation
r =[1;1;1;1];
for a = 1:simpr
m=1;
for b=1:simpr2
for n = 1:4
p = [elemx2(b,n),elemy2(b,n),elemz2(b,n)];
n1=[elemx(a,1),elemy(a,1),elemz(a,1)];
n2=[elemx(a,2),elemy(a,2),elemz(a,2)];
n3=[elemx(a,3),elemy(a,3),elemz(a,3)];
n4=[elemx(a,4),elemy(a,4),elemz(a,4)];
d0 =[n1;n2;n3;n4];
d0 =[d0 r];
d1 =[p;n2;n3;n4];
d1 =[d1 r];
d2 =[n1;p;n3;n4];
d2 =[d2 r];
d3 =[n1;n2;p;n4];
d3 =[d3 r];
d4 =[n1;n2;n3;p];
d4 =[d4 r];
deter0 = sign(det(d0));
deter1 = sign(det(d1));
deter2 = sign(det(d2));
deter3 = sign(det(d3));
deter4 = sign(det(d4));
if isequal(deter0,deter1,deter2,deter3,deter4)
elemno(a,m) = b;
m=m+1;
break;
else
continue;
end
end
end
end
