AMS 562 Final Project—Spherical Triangular

Posted 2021-09-10 位错梦的寐

AMS 562 Final Project—Spherical Triangular
Meshes
Due Dec. 20, 2018 5.00 pm
1 Introduction
Mesh is an important concept in computational science; it is used in many
different applications, for instance, most PDE based engineering problems. A
mesh is a tessellation of a domain with simple geometry objects—referred as
elements or cells.
In general, mesh can be grouped into the following two families:

1. structured grids, and
2. unstructured meshes.
   The first one is topologically equivalent to unit square (in 2D) with uniformly
   distributed (in some weighted sense) grid lines along x and y directions. Therefore,
   each of the grid points can be accessed with a pair of indices, i.e. x-/y- index,
   see below.
   Structured grids are hard to extended to general domains that are needed in
   typical engineering problems. This is particularly true for structural mechanics,
   whereas the structural components can be very complicated.
   1
   AMS 562 Final Project—Spherical Triangular Meshes 2 Data Structure
   With this consideration, people have come up with the unstructured meshes,
   which can easily represent complex geometries; see the picture below for a surface
   mesh on a dragon shape object.
   For unstructured meshes, one of the popular choices of element types for surface
   domain tessellation is triangles, because they are simple and robust (tessellating
   a surface with triangles is much easier than with other cells, e.g. quadrilateral
   cells).
   However, with unstructured grids, it’s impossible to access the grid points with
   x-, y-, (z-) index. Therefore, special data structures are needed (later).
   For this project, we will only deal with unstructured triangular
   meshes on spherical surfaces.
3. Data Structure
   2.1 Core Concept for Storing Points and Triangles
   For unstructured meshes, one cannot store, like in structured grids, xyz coordinates
   separately. Typically, a n by 3 storage is utilized, where each entry of
   the array stands for a physical point in the Cartesian coordinate system. For
   instance, you can store five arbitrary points in the following “matrix” format.
   x0 y0 z0
   x1 y1 z1
   x2 y2 z2
   x3 y3 z3
   x4 y4 z4
   In order to represent the mesh, we use a connectivity table that forms the
   connection relation between triangles and coordinates IDs. A connectivity table
   2
   AMS 562 Final Project—Spherical Triangular Meshes 2 Data Structure
   (for triangles) is an m by 3 integer storage, where m is the number of triangles,
   and the three integer IDs are the node IDs in the coordinates that form the
   triangle.
   2.2 An Example
   Let’s consider a unit square plane, we then can split the plane into two triangles
   along one of the diagonal (say, [0 0] and [1 1]). Therefore, we have 4 points and
4. triangles and the mesh can be represented by (assume xy plane)
   points =
   0.0 0.0
   1.0 0.0
   1.0 1.0
   0.0 1.0
   triangles =
   
5. 1 2
6. 3 0
   In this example, we have 4 points with the coordinates listed in points, then in
   the connectivity table, i.e. triangles, the first triangle contains nodes 0, 1 and
   2, while 2, 3 and 1 for the second one.
   Hint draw this with pencil and paper.
   2.3 Store The Mesh with C++ STL
   For both coordinates and connectivity table, we will use the std::vector as
   the container, and for each single point and triangle, we will use std::array.
   Therefore, the overall data structures are:
   // spherical coordinates
   class SphCo {
   public:
   ...
   private:
   std::vector<std::array<double, 3>> _pts;
   };
   // connectivity
   class Triangles {
   public:
   ...
   private:
   std::vector<std::array<int, 3>> _conn;
   };
   3
   AMS 562 Final Project—Spherical Triangular Meshes 3 Your Tasks
   With this implementation, we know that when you loop through the coordinates
   and connectivity table, the entries are of types std::array<double, 3> and
   std::array<int, 3>, resp.
   2.4 Remark Regarding Triangle Node Ordering
   For the spherical mesh, you can assume that the node ordering of triangles
   following the right-hand rule results vectors that always point to outward normal
   directions with respect to the spherical triangular faces.
7. Your Tasks
   For this project, you do NOT need to implement any of the data structures
   listed above. They are already implemented and ready to use. What you need
   to do is to implement several algorithms in order to perform the following tasks.
   3.1 Determine Node to Triangle Adjacency
   Essentially, a connectivity table defines the adjacency information between
   triangles and nodes, i.e. a row of the table is the adjacent nodes of that triangle.
   Starting from this, we can build its dual graph, i.e. node to cell adjacency, so
   that, given a node, we can find all triangles contain it.
   Unlike the connectivity table, for unstructured meshes, the node to cell adjacency
   information, in general, cannot be stored in a matrix-like storage, we, then, can
   use
   // vector of vector
   // where the inner vector is the list of element IDs that
   // are shared by a specific node
   std::vector<std::vector<int>> n2e_adj;
   In order to build the adjacency information, you need to implement the following
   algorithm.
   Algorithm I: Determine the node to cell adj
   Inputs: Connectivity table conn and number of points np
   Output: Node to triangle adjacency mapping, adj
   Initialize adj with size np
   for each triangle Ti
   in conn, do
   for each node vj in Ti
   , do
   4
   AMS 562 Final Project—Spherical Triangular Meshes 3 Your Tasks
   push back triangle ID i to list adjj
   end for
   end for
   3.2 Compute the Averaged Outward Normal Vectors
   Given a triangle, we can compute its normal vector by using the cross product,
   i.e. for a triangle ?ABC, the normal vector n~ is computed by
   n~ = AB~ × AC~
   where × denotes cross product.
   Since we already know that the triangle is oriented in the way that the right-hand
   rule gives the outward normal, directly using the formula above will give you
   the correct answer.
   We want first compute and store all outward normals of triangles, then use the
   adjacent mapping to average them on nodes.
   Algorithm II.1 : Compute the outward normals of spherical triangles
   Inputs: Connectivity table conn and coordinates co
   Output: Outward normals of triangles, f_nrms
   Initialize f_nrms with size ne = conn.size()
   for each triangle Ti
   in conn, do
   Get coordinates A, B, C = co(Ti,0), co(Ti,1), co(Ti,2)
   Assign f_nrmsi with the outward normal vector
   end for
   Normalize f_nrms so that each of the vector has Euclidean length 1
   In order to easily normalized the vectors, you should consider use SphCo to store
   f_nrms.
   Algorithm II.2 : Compute the averaged normal on nodes
   Inputs: Node to cell mapping adj, np, and f_nrms
   Output: Outward normals of nodes, n_nrms
   Initialize n_nrms with size np
   for each node vj , do
   Get the local adjacency information adjj
   Assign n_nrmsj to be the average of f_nrms(adjj )
   end for
   Normalize n_nrms so that each of the vector has Euclidean length 1
   5
   AMS 562 Final Project—Spherical Triangular Meshes 4 Requirements
   3.3 Determine the Computation Errors
   For points located on the spherical surface, their outward normal directions are
   uniquely defined. In this case, we want to see how much error we have introduced
   from our computation. The metric we use is
   ~a ·
   ~b = k~akk~bk cos(θ)
   where θ is the angle between ~a and ~b. We, then, have
   θ = cos?1

~a ·
~b
k~akk~bk
!
Since we know that both ~a and ~b are normalized, we can simplify the computation
by
θ = cos?1

~a ·
~b

Algorithm III: Compute the angle errors
Inputs: Analytic normals (coordinates) exact and numerical normals
num
Output: Error angles arccos_theta
Initialize arccos_theta with size exact.size()
for each node vj , do
Assign arccos_thetaj with the error angle given exactj and numj
normal vectors
end for
4 Requirements
You need to implement the algorithms list above in manip.cpp file and make
sure you pass all tests. There are three test cases, which are list below.
6
AMS 562 Final Project—Spherical Triangular Meshes 5 Hints
Note that there are 4 for loops, you must implement at least two of them with
std::for_each with lambda or traditional functors.
You will gain extra credits by using std::for_each for all iteration loops.
5 Hints
There are some sample usages of std::for_each and lambda in co.hpp and
conn.hpp. Also, type make debug to use Valgrind during development. Please
look at the comments in manip.hpp carefully.
If you get error values of nan, then it means you probably have some bug in
adj array (you have at least a node with no adjacent triangles). If you get error
values of ~ 180 degree, then it means you mess up with the node ordering so
that the right-hand rule gives you inner normal directions.

WX：codehelp