face_tri3.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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 std::unique_ptr<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 libmesh_make_unique<Side<Edge2,Tri3>>(this,i);
116 
117  else
118  {
119  std::unique_ptr<Elem> edge = libmesh_make_unique<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 edge;
127  }
128 }
129 
130 
131 void Tri3::connectivity(const unsigned int libmesh_dbg_var(sf),
132  const IOPackage iop,
133  std::vector<dof_id_type> & conn) const
134 {
135  libmesh_assert_less (sf, this->n_sub_elem());
136  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
137 
138  switch (iop)
139  {
140  case TECPLOT:
141  {
142  conn.resize(4);
143  conn[0] = this->node_id(0)+1;
144  conn[1] = this->node_id(1)+1;
145  conn[2] = this->node_id(2)+1;
146  conn[3] = this->node_id(2)+1;
147  return;
148  }
149 
150  case VTK:
151  {
152  conn.resize(3);
153  conn[0] = this->node_id(0);
154  conn[1] = this->node_id(1);
155  conn[2] = this->node_id(2);
156  return;
157  }
158 
159  default:
160  libmesh_error_msg("Unsupported IO package " << iop);
161  }
162 }
163 
164 
165 
166 
167 
168 
170 {
171  // 3-node triangles have the following formula for computing the area
172  return 0.5 * cross_norm(point(1) - point(0),
173  point(2) - point(0));
174 }
175 
176 
177 
178 std::pair<Real, Real> Tri3::min_and_max_angle() const
179 {
180  Point v10 ( this->point(1) - this->point(0) );
181  Point v20 ( this->point(2) - this->point(0) );
182  Point v21 ( this->point(2) - this->point(1) );
183 
184  const Real
185  len_10=v10.norm(),
186  len_20=v20.norm(),
187  len_21=v21.norm();
188 
189  const Real
190  theta0=std::acos(( v10*v20)/len_10/len_20),
191  theta1=std::acos((-v10*v21)/len_10/len_21),
192  theta2=libMesh::pi - theta0 - theta1
193  ;
194 
195  libmesh_assert_greater (theta0, 0.);
196  libmesh_assert_greater (theta1, 0.);
197  libmesh_assert_greater (theta2, 0.);
198 
199  return std::make_pair(std::min(theta0, std::min(theta1,theta2)),
200  std::max(theta0, std::max(theta1,theta2)));
201 }
202 
203 bool Tri3::contains_point (const Point & p, Real tol) const
204 {
205  // See "Barycentric Technique" section at
206  // http://www.blackpawn.com/texts/pointinpoly for details.
207 
208  // Compute vectors
209  Point v0 = this->point(1) - this->point(0);
210  Point v1 = this->point(2) - this->point(0);
211  Point v2 = p - this->point(0);
212 
213  // Compute dot products
214  Real dot00 = v0 * v0;
215  Real dot01 = v0 * v1;
216  Real dot02 = v0 * v2;
217  Real dot11 = v1 * v1;
218  Real dot12 = v1 * v2;
219 
220  // Out of plane check
221  if (std::abs(triple_product(v2, v0, v1)) / std::max(dot00, dot11) > tol)
222  return false;
223 
224  // Compute barycentric coordinates
225  Real invDenom = 1 / (dot00 * dot11 - dot01 * dot01);
226  Real u = (dot11 * dot02 - dot01 * dot12) * invDenom;
227  Real v = (dot00 * dot12 - dot01 * dot02) * invDenom;
228 
229  // Check if point is in triangle
230  return (u > -tol) && (v > -tol) && (u + v < 1 + tol);
231 }
232 
233 } // namespace libMesh
double abs(double a)
virtual Real volume() const libmesh_override
Definition: face_tri3.C:169
T cross_norm(const TypeVector< T > &b, const TypeVector< T > &c)
Definition: type_vector.h:1092
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
Real norm() const
Definition: type_vector.h:909
std::pair< Real, Real > min_and_max_angle() const
Definition: face_tri3.C:178
long double max(long double a, double b)
static const unsigned int side_nodes_map[3][2]
Definition: face_tri3.h:133
T triple_product(const TypeVector< T > &a, const TypeVector< T > &b, const TypeVector< T > &c)
Definition: type_vector.h:1051
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:1875
virtual bool contains_point(const Point &p, Real tol) const libmesh_override
Definition: face_tri3.C:203
subdomain_id_type subdomain_id() const
Definition: elem.h:1952
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:1810
virtual void connectivity(const unsigned int sf, const IOPackage iop, std::vector< dof_id_type > &conn) const libmesh_override
Definition: face_tri3.C:131
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const libmesh_override
Definition: face_tri3.C:99
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:1832
long double min(long double a, double b)
A geometric point in (x,y,z) space.
Definition: point.h:38
virtual std::unique_ptr< Elem > build_side_ptr(const unsigned int i, bool proxy) libmesh_override
Definition: face_tri3.C:109
virtual bool is_face(const unsigned int i) const libmesh_override
Definition: face_tri3.C:94
const Real pi
Definition: libmesh.h:172