Mathematica面临的编程挑战 [英] A programming challenge with Mathematica

查看：334
本文介绍了Mathematica面临的编程挑战的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！

问题描述

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

图1 b
$ b



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}}; 
 
显示[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中生成所需的3D几何图形 MMA格式。



几何被描述为下面的输入文件为我的外部程序。

 ＃GEOMETRY 
 ＃xyz [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.3333 0.5 
 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 
面板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 
＃输入结束文件

现在我从这个外部程序的文档描述很短。






  1. 第一个关键字 NODES 个国家/地区总数
    个节点。在此行之后，不应有任何评论或空行。下一行由
    三个值x，y和z节点坐标组成，行数必须与节点数
    相同。
  2. 下一个关键字 PANEL ，并说明我们有多少面板。之后，我们有定义每个面板的
    行。第一个整数定义面板类型
    3. - 由四个节点和四个相邻的面板定义。
    相邻面板是共享相同边（节点对）的面板，需要
    速度和压力计算（方法1和2）。丢失的邻居（例如对于后边附近的
    面板）填充值为0（请参阅图1 ）。 - 三角面板 - 由三个节点和三个相邻面板定义。 面板 - 是四边形面板，由四个节点和两个位于后边缘的
    （相邻）面板定义（面板的唤醒面板为应用Kutta条件的
    ）。 >
  3. 面板类型1和2必须在输入文件的类型10之前定义。
    需要注意的是表面法线;节点定义面板的顺序应该是逆时针
    。按照右手法则，如果手指弯曲按照编号进行操作，
    拇指将显示应该指向向外几何体的法线向量。





挑战!!



我们在一个叫做输出这些。我会留下细节给你（这是我的挑战）。



Q2

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

pre code ClearAll [polygonArea];
polygonArea [pts_List]：=
Module [{dtpts = Append [pts，pts [[1]]]}，
如果[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 [[3]]，pts [[3]] -  pts [[1]]] / 2


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

解决方案

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

这篇关于Mathematica面临的编程挑战的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！

查看全文
相关文章
其他开发最新文章
热门教程
热门工具
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆