1

我正在用 Mathematica 连接一个外部程序。我正在为外部程序创建一个输入文件。它关于将 Mathematica 生成的图形中的几何数据转换为预定义的格式。这是一个示例几何。

图1

图1

几何可以在 Mathematica 中以多种方式描述。一种费力的方法如下。

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]

GraphicsComplex这会以 MMA 格式 生成所需的 3D 几何图形。在此处输入图像描述

该几何图形被描述为我的外部程序的以下输入文件。

# 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

现在我从这个外部程序的文档中得到的描述很短。我在这里引用它。

  1. 第一个关键字NODES表示节点总数。此行之后不应有注释或空行。下一行由三个值 x、y 和 z 节点坐标组成，行数必须与节点数相同。
  2. 下一个关键字是PANEL并说明我们有多少个面板。之后，我们有定义每个面板的行。第一个整数定义面板类型
  3. ID 1 -四边形面板- 由四个节点和四个相邻面板定义。相邻面板是共享相同侧面（节点对）的面板，并且需要用于速度和压力计算（方法 1 和 2）。缺失的邻居（例如后缘附近的面板）用值 0 填充（参见图 1）。
  4. ID 2 –三角形面板– 由三个节点和三个相邻面板定义。
  5. ID 10 –尾流面板– 是由四个节点和两个（相邻）面板定义的四边形面板，它们位于后缘（尾流面板应用 Kutta 条件的面板）。
  6. 面板类型 1 和 2 必须在输入文件中的类型 10 之前定义。需要注意的是表面法线；定义面板的节点顺序应该是逆时针的。根据右手法则，如果手指弯曲以跟随编号，拇指将显示应该指向“向外”几何图形的法线向量。

挑战！！

我们在一个名为One.obj的文件中给出了一个 3D CAD 模型，它可以在 MMA 中很好地导出。

cd = Import["One.obj"]

输出是一个 MMAGraphics3D对象 在此处输入图像描述

现在，当 MMA 在内部读取几​​何数据时，我可以轻松访问它们。

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

在此处输入图像描述

  1. 我们如何使用包含在文本文件中的顶点和多边形信息，ver1并将pol1它们写入文本文件，如上面输入文件示例中所述。在这种情况下，我们将只有ID2类型（三角形）面板。
  2. 使用 Mathematica 三角测量如何找到这个 3D 对象的表面积。是否有任何内置函数可以计算 MMA 中的表面积？
  3. 现在不需要创建唤醒面板或ID10类型的元素。只有三角形元素的输入文件就可以了。

很抱歉这么长的帖子，但这是一个我试图解决很长时间的难题。希望你们中的一些专家可能有正确的见解来破解它。

BR

4

1 回答 1

5

Q1 和 Q2 很简单，您可以在问题中去掉“挑战”标签。Q3 可以做一些澄清。

第一季度 

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

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

更新第一季度

主要问题是找到每个多边形的邻居。以下是这样做的：

(* 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 *)

您可以使用此处描述的低电平输出功能来输出这些。我会把细节留给你（这是我对你的挑战）。

Q2
机翼的面积是各个多边形的总面积。各个区域可以计算如下：

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}]
   ]
 ]

基于这个 Mathworld 页面。

该区域已签署 BTW，因此您可能需要使用Abs.

更正
上述面积函数仅适用于2D中的一般多边形。对于3D三角形的面积，可以使用以下公式：

ClearAll[polygonArea];
polygonArea[pts_List?(Length[#] == 3 &)] := 
    Norm[Cross[pts[[2]] - pts[[1]], pts[[3]] - pts[[1]]]]/2
于 2011-09-14T20:51:14.680 回答