A programming challenge with Mathematica

Question

I am interfacing an external program with Mathematica. I am creating an input file for the external program. Its about converting geometry data from a Mathematica generated graphics into a predefined format. Here is an example Geometry.

Figure 1

The geometry can be described in many ways in Mathematica. One laborious way is the following.

dat={{1.,-1.,0.},{0.,-1.,0.5},{0.,-1.,-0.5},{1.,-0.3333,0.},{0.,-0.3333,0.5},
{0.,-0.3333,-0.5},{1.,0.3333,0.},{0.,0.3333,0.5},{0.,0.3333,-0.5},{1.,1.,0.},
{0.,1.,0.5},{0.,1.,-0.5},{10.,-1.,0.},{10.,-0.3333,0.},{10.,0.3333,0.},{10.,1.,0.}};

Show[ListPointPlot3D[dat,PlotStyle->{{Red,PointSize[Large]}}],Graphics3D[{Opacity[.8],
Cyan,GraphicsComplex[dat,Polygon[{{1,2,5,4},{1,3,6,4},{2,3,6,5},{4,5,8,7},{4,6,9,7},
{5,6,9,8},{7,8,11,10},{7,9,12,10},{8,9,12,11},{1,2,3},{10,12,11},{1,4,14,13},
{4,7,15,14},{7,10,16,15}}]]}],AspectRatio->GoldenRatio]

This generates the required 3D geometry in GraphicsComplex format of MMA.

This geometry is described as the following input file for my external program.

# GEOMETRY
# x y z [m]
NODES 16
1. -1. 0.
0. -1. 0.5
0. -1. -0.5
1. -0.3333 0.
0. -0.3333 0.50. -0.3333 -0.5
1. 0.3333 0.
0. 0.3333 0.5
0. 0.3333 -0.5
1. 1. 0.
0. 1. 0.5
0. 1. -0.5
10. -1. 0.
10. -0.3333 0.
10. 0.3333 0.
10. 1. -0.
# type node_id1 node_id2 node_id3 node_id4  elem_id1 elem_id2 elem_id3 elem_id4
PANELS 14
1 1 4 5 2 4 2 10 0
1 2 5 6 3 1 5 3 10
1 3 6 4 1 2 6 10 0
1 4 7 8 5 7 5 1 0
1 5 8 9 6 4 8 6 2
1 6 9 7 4 5 9 3 0
1 7 10 11 8 8 4 11 0
1 8 11 12 9 7 9 5 11
1 9 12 10 7 8 6 11 0
2 1 2 3 1 2 3
2 10 12 11 9 8 7
10 4 1 13 14 1 3
10 7 4 14 15 4 6
10 10 7 15 16 7 9
# end of input file

Now the description I have from the documentation of this external program is pretty short. I am quoting it here.

---

1. First keyword NODES states total number of nodes. After this line there should be no comment or empty lines. Next lines consist of three values x, y and z node coordinates and number of lines must be the same as number of nodes.
2. Next keyword is PANEL and states how many panels we have. After that we have lines defining each panel. First integer defines panel type
3. ID 1 – quadrilateral panel - is defined by four nodes and four neighboring panels. Neighboring panels are panels that share same sides (pair of nodes) and is needed for velocity and pressure calculation (methods 1 and 2). Missing neighbors (for example for panels near the trailing edge) are filled with value 0 (see Figure 1).
4. ID 2 – triangular panel – is defined by three nodes and three neighboring panels.
5. ID 10 – wake panel – is quadrilateral panel defined with four nodes and with two (neighboring) panels which are located on the trailing edge (panels to which wake panel is applying Kutta condition).
6. Panel types 1 and 2 must be defined before type 10 in input file. Important to notice is the surface normal; order of nodes defining panels should be counter clockwise. By the right-hand rule if fingers are bended to follow numbering, thumb will show normal vector that should point “outwards” geometry.

---

Challenge!!

We are given with a 3D CAD model in a file called One.obj and it is exported fine in MMA.

cd = Import["One.obj"]

The output is a MMA Graphics3D object

Now I can get easily access the geometry data as MMA internally reads them.

{ver1, pol1} = cd[[1]][[2]] /. GraphicsComplex -> List;
MyPol = pol1 // First // First;
Graphics3D[GraphicsComplex[ver1,MyPol],Axes-> True]

1. How we can use the vertices and polygon information contained in ver1 and pol1 and write them in a text file as described in the input file example above. In this case we will only have ID2 type (triangular) panels.
2. Using the Mathematica triangulation how to find the surface area of this 3D object. Is there any inbuilt function that can compute surface area in MMA?
3. No need to create the wake panel or ID10 type elements right now. A input file with only triangular elements will be fine.

---

Sorry for such a long post but its a puzzle that I am trying to solve for a long time. Hope some of you expert may have the right insight to crack it.

BR

Chadchadabe

I usually try to get paid just for reading problems like this one. – Kasiekask

@belisarius Unfortunately I am from a much more underprivileged society where even an elaborate solution to this problem will surly fail to fetch any payment for the solver. But once again my apology for such a lengthy question. – Chadchadabe

Are you aware that Mathematica 7 (and maybe earlier) supports both import and export of .obj format graphics files? – Droit

@DeepYellow I was sure about it. But export import for .obj files does not solve my problem.. – Chadchadabe

Answer 1

Q1 and Q2 are easy enough that you could drop the "challenge" labels in your question. Q3 could use some clarification.

Q1

edges = cd[[1, 2, 1]];

polygons = cd[[1, 2, 2, 1, 1, 1]];

---

Update Q1

The main problem is to find the neighbor of each polygon. The following does this:

(* Split every triangle in 3 edges, with nodes in each edge sorted *)
triangleEdges = (Sort /@ Subsets[#, {2}]) & /@ polygons;

(* Generate a list of edges *)
singleEdges = Union[Flatten[triangleEdges, 1]];

(* Define a function which, given an edge (node number list), returns the bordering  *)
(* triangle numbers. It's done by working through each of the triangles' edges       *)
ClearAll[edgesNeighbors]
edgesNeighbors[_] = {};
MapIndexed[(
   edgesNeighbors[#1[[1]]] = Flatten[{edgesNeighbors[#1[[1]]], #2[[1]]}];
   edgesNeighbors[#1[[2]]] = Flatten[{edgesNeighbors[#1[[2]]], #2[[1]]}];
   edgesNeighbors[#1[[3]]] = Flatten[{edgesNeighbors[#1[[3]]], #2[[1]]}];
   ) &, triangleEdges
];

(* Build a triangle relation table. Each '1' indicates a triangle relation *)
relations = ConstantArray[0, {triangleEdges // Length, triangleEdges // Length}];
Scan[
  (n = edgesNeighbors[##]; 
     If[Length[n] == 2, 
        {n1, n2} = n; 
        relations[[n1, n2]] = 1;  relations[[n2, n1]] = 1];
   ) &, singleEdges
]

MatrixPlot[relations]

(* Build a neighborhood list *)
triangleNeigbours = 
    Table[Flatten[Position[relations[[i]], 1]], {i,triangleEdges // Length}];

(* Test: Which triangles border on triangle number 1? *)
triangleNeigbours[[1]]

(* ==> {32, 61, 83} *)

(* Check this *)
polygons[[{1, 32, 61, 83}]]

(* ==> {{1, 2, 3}, {3, 2, 52}, {1, 3, 50}, {19, 2, 1}} *)
(* Indeed, they all share an edge with #1 *)

---

You can use the low level output functions described here to output these. I'll leave the details to you (that's my challenge to you).

Q2
The area of the wing is the summed area of the individual polygons. The individual areas can be calculated as follows:

ClearAll[polygonArea];
polygonArea[pts_List] :=
 Module[{dtpts = Append[pts, pts[[1]]]},
   If[Length[pts] < 3, 
      0, 
      1/2 Sum[Det[{dtpts[[i]], dtpts[[i + 1]]}], {i, 1, Length[dtpts] - 1}]
   ]
 ]

based on this Mathworld page.

The area is signed BTW, so you may want to use Abs.

CORRECTION
The above area function is only usable for general polygons in 2D. For the area of a triangle in 3D the following can be used:

ClearAll[polygonArea];
polygonArea[pts_List?(Length[#] == 3 &)] := 
    Norm[Cross[pts[[2]] - pts[[1]], pts[[3]] - pts[[1]]]]/2

Pedagogy

Question 1 is not at all answered just by getting the edges and the polygons.It is already extracted when we form the GraphicsComplex. The main problem is that one needs to give also information how edges in each polygon is connected to other polygons. Look at the element columns of the example input file. So we need kind of marking for each polygon. Check the red marked polygons in Figure 1. – Chadchadabe

Thank you very much for the Polygon area stuff. – Chadchadabe

@Chadchadabe My apologies for reading your question too superficially. I'll dive deeper, but that'll take some time. In the meantime a question: shouldn't the 0 in the line '1 1 4 5 2 4 2 10 0' directly following PANEL be a 12 (or even a 3)? Both surfaces seem to be connected to nodes 1 and 4 and are therefore neighboring panels. – Pedagogy

Sorry to write back after such a long gap. Got a injury while playing and had a bad time really. By the way your solution is still not working for the external solver. I think the panel ordering is not going correctly. Here is an example input file that works. It is a sphere with 24 panels. mediafire.com/?0hu5ql817h2ghiy . Here the use quad panels meaning in place of triangles they have quadrilateral. Hope you help me... – Chadchadabe

@Chadchadabe Without details about how you used my functions to generate the panel list I cannot provide any help. I assume the nodes of a panel should be ordered counterclockwise and I assume the neighboring planes should be ordered in the same order as the edges. There is still the question above you haven't answered. about edges that have more than 1 connecting panels, like panels 1, 3 and 12? You really have to figure that out yourself. The details are not a Mathematica problem and I unfortunately don't have much spare time left. – Pedagogy

if you are interested check #7676532 if you can give some more insight there. – Chadchadabe