face_tri3.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 // C++ includes
19 
20 // Local includes
21 #include "libmesh/side.h"
22 #include "libmesh/edge_edge2.h"
23 #include "libmesh/face_tri3.h"
24 
25 namespace libMesh
26 {
27 
28 
29 
30 // ------------------------------------------------------------
31 // Tri3 class static member initializations
32 const unsigned int Tri3::side_nodes_map[3][2] =
33  {
34  {0, 1}, // Side 0
35  {1, 2}, // Side 1
36  {2, 0} // Side 2
37  };
38 
39 
40 #ifdef LIBMESH_ENABLE_AMR
41 
42 const float Tri3::_embedding_matrix[4][3][3] =
43  {
44  // embedding matrix for child 0
45  {
46  // 0 1 2
47  {1.0, 0.0, 0.0}, // 0
48  {0.5, 0.5, 0.0}, // 1
49  {0.5, 0.0, 0.5} // 2
50  },
51 
52  // embedding matrix for child 1
53  {
54  // 0 1 2
55  {0.5, 0.5, 0.0}, // 0
56  {0.0, 1.0, 0.0}, // 1
57  {0.0, 0.5, 0.5} // 2
58  },
59 
60  // embedding matrix for child 2
61  {
62  // 0 1 2
63  {0.5, 0.0, 0.5}, // 0
64  {0.0, 0.5, 0.5}, // 1
65  {0.0, 0.0, 1.0} // 2
66  },
67 
68  // embedding matrix for child 3
69  {
70  // 0 1 2
71  {0.5, 0.5, 0.0}, // 0
72  {0.0, 0.5, 0.5}, // 1
73  {0.5, 0.0, 0.5} // 2
74  }
75  };
76 
77 #endif
78 
79 
80 
81 // ------------------------------------------------------------
82 // Tri3 class member functions
83 
84 bool Tri3::is_vertex(const unsigned int) const
85 {
86  return true;
87 }
88 
89 bool Tri3::is_edge(const unsigned int) const
90 {
91  return false;
92 }
93 
94 bool Tri3::is_face(const unsigned int) const
95 {
96  return false;
97 }
98 
99 bool Tri3::is_node_on_side(const unsigned int n,
100  const unsigned int s) const
101 {
102  libmesh_assert_less (s, n_sides());
103  for (unsigned int i = 0; i != 2; ++i)
104  if (side_nodes_map[s][i] == n)
105  return true;
106  return false;
107 }
108 
109 UniquePtr<Elem> Tri3::build_side_ptr (const unsigned int i,
110  bool proxy)
111 {
112  libmesh_assert_less (i, this->n_sides());
113 
114  if (proxy)
115  return UniquePtr<Elem>(new Side<Edge2,Tri3>(this,i));
116 
117  else
118  {
119  Elem * edge = new Edge2;
120  edge->subdomain_id() = this->subdomain_id();
121 
122  // Set the nodes
123  for (unsigned n=0; n<edge->n_nodes(); ++n)
124  edge->set_node(n) = this->node_ptr(Tri3::side_nodes_map[i][n]);
125 
126  return UniquePtr<Elem>(edge);
127  }
128 
129  libmesh_error_msg("We'll never get here!");
130  return UniquePtr<Elem>();
131 }
132 
133 
134 void Tri3::connectivity(const unsigned int libmesh_dbg_var(sf),
135  const IOPackage iop,
136  std::vector<dof_id_type> & conn) const
137 {
138  libmesh_assert_less (sf, this->n_sub_elem());
139  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
140 
141  switch (iop)
142  {
143  case TECPLOT:
144  {
145  conn.resize(4);
146  conn[0] = this->node_id(0)+1;
147  conn[1] = this->node_id(1)+1;
148  conn[2] = this->node_id(2)+1;
149  conn[3] = this->node_id(2)+1;
150  return;
151  }
152 
153  case VTK:
154  {
155  conn.resize(3);
156  conn[0] = this->node_id(0);
157  conn[1] = this->node_id(1);
158  conn[2] = this->node_id(2);
159  return;
160  }
161 
162  default:
163  libmesh_error_msg("Unsupported IO package " << iop);
164  }
165 }
166 
167 
168 
169 
170 
171 
173 {
174  // 3-node triangles have the following formula for computing the area
175  return 0.5 * cross_norm(point(1) - point(0),
176  point(2) - point(0));
177 }
178 
179 
180 
181 std::pair<Real, Real> Tri3::min_and_max_angle() const
182 {
183  Point v10 ( this->point(1) - this->point(0) );
184  Point v20 ( this->point(2) - this->point(0) );
185  Point v21 ( this->point(2) - this->point(1) );
186 
187  const Real
188  len_10=v10.norm(),
189  len_20=v20.norm(),
190  len_21=v21.norm();
191 
192  const Real
193  theta0=std::acos(( v10*v20)/len_10/len_20),
194  theta1=std::acos((-v10*v21)/len_10/len_21),
195  theta2=libMesh::pi - theta0 - theta1
196  ;
197 
198  libmesh_assert_greater (theta0, 0.);
199  libmesh_assert_greater (theta1, 0.);
200  libmesh_assert_greater (theta2, 0.);
201 
202  return std::make_pair(std::min(theta0, std::min(theta1,theta2)),
203  std::max(theta0, std::max(theta1,theta2)));
204 }
205 
206 bool Tri3::contains_point (const Point & p, Real tol) const
207 {
208  // See "Barycentric Technique" section at
209  // http://www.blackpawn.com/texts/pointinpoly for details.
210 
211  // Compute vectors
212  Point v0 = this->point(1) - this->point(0);
213  Point v1 = this->point(2) - this->point(0);
214  Point v2 = p - this->point(0);
215 
216  // Compute dot products
217  Real dot00 = v0 * v0;
218  Real dot01 = v0 * v1;
219  Real dot02 = v0 * v2;
220  Real dot11 = v1 * v1;
221  Real dot12 = v1 * v2;
222 
223  // Out of plane check
224  if (std::abs(triple_product(v2, v0, v1)) / std::max(dot00, dot11) > tol)
225  return false;
226 
227  // Compute barycentric coordinates
228  Real invDenom = 1 / (dot00 * dot11 - dot01 * dot01);
229  Real u = (dot11 * dot02 - dot01 * dot12) * invDenom;
230  Real v = (dot00 * dot12 - dot01 * dot02) * invDenom;
231 
232  // Check if point is in triangle
233  return (u > -tol) && (v > -tol) && (u + v < 1 + tol);
234 }
235 
236 } // namespace libMesh
double abs(double a)
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:1795
virtual Real volume() const libmesh_override
Definition: face_tri3.C:172
T cross_norm(const TypeVector< T > &b, const TypeVector< T > &c)
Definition: type_vector.h:1084
virtual bool is_edge(const unsigned int i) const libmesh_override
Definition: face_tri3.C:89
virtual unsigned int n_sub_elem() const libmesh_override
Definition: face_tri3.h:73
The base class for all geometric element types.
Definition: elem.h:86
Real norm() const
Definition: type_vector.h:903
std::pair< Real, Real > min_and_max_angle() const
Definition: face_tri3.C:181
long double max(long double a, double b)
static const unsigned int side_nodes_map[3][2]
Definition: face_tri3.h:133
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
virtual unsigned int n_nodes() const =0
T triple_product(const TypeVector< T > &a, const TypeVector< T > &b, const TypeVector< T > &c)
Definition: type_vector.h:1043
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:1730
virtual bool contains_point(const Point &p, Real tol) const libmesh_override
Definition: face_tri3.C:206
subdomain_id_type subdomain_id() const
Definition: elem.h:1805
PetscErrorCode Vec Mat libmesh_dbg_var(j)
Proxy class for efficiently representing an Elem&#39;s side.
Definition: side.h:48
virtual UniquePtr< Elem > build_side_ptr(const unsigned int i, bool proxy) libmesh_override
Definition: face_tri3.C:109
static const float _embedding_matrix[4][3][3]
Definition: face_tri3.h:177
virtual unsigned int n_sides() const libmesh_override
Definition: face_tri.h:86
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Point & point(const unsigned int i) const
Definition: elem.h:1667
virtual void connectivity(const unsigned int sf, const IOPackage iop, std::vector< dof_id_type > &conn) const libmesh_override
Definition: face_tri3.C:134
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const libmesh_override
Definition: face_tri3.C:99
A 1D geometric element with 2 nodes.
Definition: edge_edge2.h:43
virtual bool is_vertex(const unsigned int i) const libmesh_override
Definition: face_tri3.C:84
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:1689
long double min(long double a, double b)
A geometric point in (x,y,z) space.
Definition: point.h:38
virtual bool is_face(const unsigned int i) const libmesh_override
Definition: face_tri3.C:94
const Real pi
Definition: libmesh.h:172