@vv  Thank you for the solution!

@Rouben Rostamian  
Thank you for your interest and original solution!

Here is my short manual solution without trigonometry for an arbitrary right triangle:

We use the well-known formula for the radius of a circle inscribed in a right triangle  (a + b - c) / 2 . Further (see the figure) triangles  AEF = FKC , since they are similar, and the inscribed circles in them are the same. Therefore,  F  divides the segment  AC  in the ratio a : b . We get  AF = a*c / (a + b) . The similarity coefficient of triangle  AEF  with respect to triangle  ABC  is  k = c / (a+b) .
So we get the final answer  r = (a+b-c)*c / (2*(a+b))   (for the right triangle with sides 3, 4, 5 radius r = 5 / 7 )

Of course my and your answers for right triangles are the same:

is((a+b-c)*c/(a+b)/2=c/(a/(c-b)+2+b/(c-a))) assuming c^2=a^2+b^2;

                                       true
  

@minhthien2016  It is not difficult to do. I already wrote in one of the comments that each element of the list  L , i.e.  [x, y, z, u, v, w]  means a tetrahedron, in which  x, y, z  are the lengths of 3 edges coming out of one vertex,  u is opposite z ,  v is opposite y , w is opposite x  . 

@aaadoors  The condition you are writing about is necessary but not sufficient for the existence of such a tetrahedron. Read about the Cayley-Menger determinant in Wikipedia  https://en.wikipedia.org/wiki/Cayley%E2%80%93Menger_determinant  

Counterexample to your statement:
Bottom base -  (15,16,17).  Side faces -  (15,8,9), (17,9,10), (16,10,8)

@minhthien  You asked about the number of all solutions. The following code finds this number. This number is 4 times smaller than the length of the list  L , because each solution is represented in this list by 4 copies, depending on the choice of the tetrahedron vertex (the first 3 numbers in each sublist of  L  are the lengths of the 3 edges coming out of one vertex).

restart;
with(combinat):
primes := [seq(ithprime(i), i = 1 .. 20)];
S:=choose(primes, 6):
F:=(x,y,z,u,v,w)->`if`(LinearAlgebra:-Determinant(<0,1,1,1,1; 1,0,x^2,y^2,z^2; 1,x^2,0,u^2,v^2; 1,y^2,u^2,0,w^2; 1,z^2,v^2,w^2,0>)>0,true,false): k:=0:
for s in S do
for t in [seq(seq([s1[],s2[]], s2=permute(convert(s,set) minus convert(s1,set))), s1=choose(s, 3))] do
x,y,z,u,v,w:=t[];
if F(t[]) and min(x+y+z,x+v+u,w+y+u,w+v+z)>max(x+w,y+v,z+u) then k:=k+1; L[k]:=t;  fi;
od: od:
L:=convert(L, list): 
k/4;

@vv  Now I understand, thank you.

@vv  You still haven't explained why your code is skipping solutions (see my comment above "Missed solutions")

@vv  My solution  [3, 5, 13, 7, 11, 17]  is correct. Below is a detailed check. Tetrahedron ABCS with vertex S and base ABC. AS=3, BS=13, CS=5, AB=11, AC=7, BC=17. There is no triangle (5,13,7) here.

restart; # Check solution [3, 5, 13, 7, 11, 17]
f:=(X,Y)->sqrt((X[1]-Y[1])^2+(X[2]-Y[2])^2+(X[3]-Y[3])^2): # Distance beetween points X and Y
A:=[0,0,0]: B:=[x,y,0]: C:=[7,0,0]: 
sol1:=solve({f(A,B)=11,f(B,C)=17,y>0},explicit);
B:=eval([x,y,0],sol1);
S:=[x,y,z]:
sol2:=solve({f(A,S)=3,f(B,S)=13,f(C,S)=5,z>0},explicit);
S:=eval([x,y,z],sol2);
with(plottools): with(plots):
T:=textplot3d([[A[],"A"],[B[],"B"],[C[],"C"],[S[],"S"]], color=black, font=[times,16], align={left,above}):
display(curve([A,B,C,A,S,C]),line(S,B),T, color=red, thickness=2, scaling=constrained, axes=none);
f(A,B), f(A,C), f(B,C), f(A,S), f(B,S), f(C,S); # Check of lengths of sides

@vv   It seems your code is missing some solutions. In my answer for the first 10 primes, the first solution found is  [ 3, 5, 13, 7, 11, 17]   (22 solutions in total). But when I ran your code (slightly modified) to solve the same problem, I only got 5 solutions, and the solution above is missing:

restart;
primes := [seq(ithprime(i), i = 1 .. 10)]:
with(combinat):
comb:= choose(primes, 6):
IsTetra:=proc(L::list(posint))
  local d:=Matrix(4, {(1,2)=L[1], (1,3)=L[2], (2,3)=L[3], (1,4)=L[4], (2,4)=L[5], (3,4)=L[6]}, shape='symmetric');
  Matrix(3, (i,j)-> d[4,i]^2+d[4,j]^2-d[i,j]^2);  # Schoenberg
  LinearAlgebra:-IsDefinite(%, query = 'positive_definite')
end proc:
n:=0:
for s in comb do
  if IsTetra(s) then 
     n:=n+1; L[n]:=s;
  fi:
od:
L:=convert(L,list);
nops(L);

@nm  You have not responded to my answer in any way, probably you are not interested in this approach. I still think that the image of several stream lines together with several arrows on them very clearly characterizes the vector field for the corresponding system of differential equations. I have significantly edited my answer, adding several arrows.

@vv  Edges  x, y, z  come out from the same vertex, u is opposite z ,  v is opposite y , w is opposite x  . 

@Rouben Rostamian  Thanks for finding this example!  I tested the  Into_4_Equal_Areas  procedure on this triangle and it found both solutions. To get symbolic solutions I set in the code Digits=15 and then set in identify command to evalf[8] , because unfortunately fsolve doesn't work very well, when calculating with 15  digits it returns  0.166666631418626  instead  1/6=0.166666666666667 .

L:=[[[0,0],[1,0]],[[1,0],[1/2,8/9]],[[1/2,8/9],[0,0]]];
Sol:=Into_4_Equal_Areas(L);
identify(evalf[8](Sol));
plots:-display(<Pic(L,Sol[1]) | Pic(L,Sol[2])>);

@Alfred_F  Thank you for your interest! This is a very important question about the number of essential solutions for a given flat region (solutions are called essential if they do not transform into one another during movements). It seems that for regular polygons, whose number of sides is divisible by 4, there are infinitely many essential solutions. I do not know of any other example where the number of essential solutions would be greater than 1.

@Rouben Rostamian  Thank you for your interest and high rating. The procedure will work correctly if there is a solution in which the cutting lines intersect the boundary of the domain at exactly 4 points. If the domain is "strongly non-convex", then this condition is not satisfied. For example, in the example below, a solution exists (it is easy to find without using the procedure), but the procedure does not work here, because the cutting lines intersect the boundary at more than 4 points.

The latest result for Maple <=2021 is also not quite correct, part of the parabola lies above the straight line (plotted in Maple 2018.2):

plots:-display(plots:-implicitplot(y=x^2, x=-1..1, y=0..(x+1)/2),plot((x+1)/2, x=-1..1), color=[red,green]);

Edit. Here is another example of incorrect plotting (Maple 2018.2):

restart;
plots:-display(plots:-implicitplot(y=-4*x^2+1, x=-1..1, y=0..(x+1)/2),plot((x+1)/2, x=-1..1), color=[red,green], scaling=constrained, size=[500,300]);