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 // Local includes
19 #include "libmesh/side.h"
20 #include "libmesh/edge_edge2.h"
21 #include "libmesh/face_tri3.h"
23 #include "libmesh/enum_order.h"
24 
25 namespace libMesh
26 {
27 
28 
29 
30 // ------------------------------------------------------------
31 // Tri3 class static member initializations
32 const int Tri3::num_nodes;
33 const int Tri3::num_sides;
34 const int Tri3::num_children;
35 const int Tri3::nodes_per_side;
36 
38  {
39  {0, 1}, // Side 0
40  {1, 2}, // Side 1
41  {2, 0} // Side 2
42  };
43 
44 
45 #ifdef LIBMESH_ENABLE_AMR
46 
48  {
49  // embedding matrix for child 0
50  {
51  // 0 1 2
52  {1.0, 0.0, 0.0}, // 0
53  {0.5, 0.5, 0.0}, // 1
54  {0.5, 0.0, 0.5} // 2
55  },
56 
57  // embedding matrix for child 1
58  {
59  // 0 1 2
60  {0.5, 0.5, 0.0}, // 0
61  {0.0, 1.0, 0.0}, // 1
62  {0.0, 0.5, 0.5} // 2
63  },
64 
65  // embedding matrix for child 2
66  {
67  // 0 1 2
68  {0.5, 0.0, 0.5}, // 0
69  {0.0, 0.5, 0.5}, // 1
70  {0.0, 0.0, 1.0} // 2
71  },
72 
73  // embedding matrix for child 3
74  {
75  // 0 1 2
76  {0.5, 0.5, 0.0}, // 0
77  {0.0, 0.5, 0.5}, // 1
78  {0.5, 0.0, 0.5} // 2
79  }
80  };
81 
82 #endif
83 
84 
85 
86 // ------------------------------------------------------------
87 // Tri3 class member functions
88 
89 bool Tri3::is_vertex(const unsigned int) const
90 {
91  return true;
92 }
93 
94 bool Tri3::is_edge(const unsigned int) const
95 {
96  return false;
97 }
98 
99 bool Tri3::is_face(const unsigned int) const
100 {
101  return false;
102 }
103 
104 bool Tri3::is_node_on_side(const unsigned int n,
105  const unsigned int s) const
106 {
107  libmesh_assert_less (s, n_sides());
108  return std::find(std::begin(side_nodes_map[s]),
110  n) != std::end(side_nodes_map[s]);
111 }
112 
113 std::vector<unsigned>
114 Tri3::nodes_on_side(const unsigned int s) const
115 {
116  libmesh_assert_less(s, n_sides());
117  return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s])};
118 }
119 
121 {
122  return FIRST;
123 }
124 
125 std::unique_ptr<Elem> Tri3::build_side_ptr (const unsigned int i,
126  bool proxy)
127 {
128  libmesh_assert_less (i, this->n_sides());
129 
130  if (proxy)
131  return libmesh_make_unique<Side<Edge2,Tri3>>(this,i);
132 
133  else
134  {
135  std::unique_ptr<Elem> edge = libmesh_make_unique<Edge2>();
136  edge->subdomain_id() = this->subdomain_id();
137 
138  // Set the nodes
139  for (unsigned n=0; n<edge->n_nodes(); ++n)
140  edge->set_node(n) = this->node_ptr(Tri3::side_nodes_map[i][n]);
141 
142  return edge;
143  }
144 }
145 
146 
147 
148 void Tri3::build_side_ptr (std::unique_ptr<Elem> & side,
149  const unsigned int i)
150 {
151  this->simple_build_side_ptr<Tri3>(side, i, EDGE2);
152 }
153 
154 
155 
156 void Tri3::connectivity(const unsigned int libmesh_dbg_var(sf),
157  const IOPackage iop,
158  std::vector<dof_id_type> & conn) const
159 {
160  libmesh_assert_less (sf, this->n_sub_elem());
161  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
162 
163  switch (iop)
164  {
165  case TECPLOT:
166  {
167  conn.resize(4);
168  conn[0] = this->node_id(0)+1;
169  conn[1] = this->node_id(1)+1;
170  conn[2] = this->node_id(2)+1;
171  conn[3] = this->node_id(2)+1;
172  return;
173  }
174 
175  case VTK:
176  {
177  conn.resize(3);
178  conn[0] = this->node_id(0);
179  conn[1] = this->node_id(1);
180  conn[2] = this->node_id(2);
181  return;
182  }
183 
184  default:
185  libmesh_error_msg("Unsupported IO package " << iop);
186  }
187 }
188 
189 
190 
191 
192 
193 
195 {
196  // 3-node triangles have the following formula for computing the area
197  return 0.5 * cross_norm(point(1) - point(0),
198  point(2) - point(0));
199 }
200 
201 
202 
203 std::pair<Real, Real> Tri3::min_and_max_angle() const
204 {
205  Point v10 ( this->point(1) - this->point(0) );
206  Point v20 ( this->point(2) - this->point(0) );
207  Point v21 ( this->point(2) - this->point(1) );
208 
209  const Real
210  len_10=v10.norm(),
211  len_20=v20.norm(),
212  len_21=v21.norm();
213 
214  const Real
215  theta0=std::acos(( v10*v20)/len_10/len_20),
216  theta1=std::acos((-v10*v21)/len_10/len_21),
217  theta2=libMesh::pi - theta0 - theta1
218  ;
219 
220  libmesh_assert_greater (theta0, 0.);
221  libmesh_assert_greater (theta1, 0.);
222  libmesh_assert_greater (theta2, 0.);
223 
224  return std::make_pair(std::min(theta0, std::min(theta1,theta2)),
225  std::max(theta0, std::max(theta1,theta2)));
226 }
227 
228 bool Tri3::contains_point (const Point & p, Real tol) const
229 {
230  // See "Barycentric Technique" section at
231  // http://www.blackpawn.com/texts/pointinpoly for details.
232 
233  // Compute vectors
234  Point v0 = this->point(1) - this->point(0);
235  Point v1 = this->point(2) - this->point(0);
236  Point v2 = p - this->point(0);
237 
238  // Compute dot products
239  Real dot00 = v0 * v0;
240  Real dot01 = v0 * v1;
241  Real dot02 = v0 * v2;
242  Real dot11 = v1 * v1;
243  Real dot12 = v1 * v2;
244 
245  // Out of plane check
246  if (std::abs(triple_product(v2, v0, v1)) / std::max(dot00, dot11) > tol)
247  return false;
248 
249  // Compute barycentric coordinates
250  Real invDenom = 1 / (dot00 * dot11 - dot01 * dot01);
251  Real u = (dot11 * dot02 - dot01 * dot12) * invDenom;
252  Real v = (dot00 * dot12 - dot01 * dot02) * invDenom;
253 
254  // Check if point is in triangle
255  return (u > -tol) && (v > -tol) && (u + v < 1 + tol);
256 }
257 
260 {
261  return Elem::loose_bounding_box();
262 }
263 
264 
265 } // namespace libMesh
double abs(double a)
static const int num_nodes
Definition: face_tri3.h:148
virtual unsigned int n_sub_elem() const override
Definition: face_tri3.h:81
Real norm() const
Definition: type_vector.h:912
static const unsigned int side_nodes_map[num_sides][nodes_per_side]
Definition: face_tri3.h:157
static const int num_sides
Definition: face_tri3.h:149
T cross_norm(const TypeVector< T > &b, const TypeVector< T > &c)
Definition: type_vector.h:1097
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const override
Definition: face_tri3.C:104
unsigned short int side
Definition: xdr_io.C:50
virtual BoundingBox loose_bounding_box() const
Definition: elem.C:2857
std::pair< Real, Real > min_and_max_angle() const
Definition: face_tri3.C:203
virtual std::vector< unsigned int > nodes_on_side(const unsigned int s) const override
Definition: face_tri3.C:114
virtual BoundingBox loose_bounding_box() const override
Definition: face_tri3.C:259
IterBase * end
virtual void connectivity(const unsigned int sf, const IOPackage iop, std::vector< dof_id_type > &conn) const override
Definition: face_tri3.C:156
long double max(long double a, double b)
virtual bool is_edge(const unsigned int i) const override
Definition: face_tri3.C:94
virtual unsigned int n_sides() const override final
Definition: face_tri.h:92
T triple_product(const TypeVector< T > &a, const TypeVector< T > &b, const TypeVector< T > &c)
Definition: type_vector.h:1054
virtual bool contains_point(const Point &p, Real tol) const override
Definition: face_tri3.C:228
static const int nodes_per_side
Definition: face_tri3.h:151
virtual std::unique_ptr< Elem > build_side_ptr(const unsigned int i, bool proxy) override
Definition: face_tri3.C:125
virtual bool is_vertex(const unsigned int i) const override
Definition: face_tri3.C:89
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
subdomain_id_type subdomain_id() const
Definition: elem.h:2034
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:1957
virtual Order default_order() const override
Definition: face_tri3.C:120
static const int num_children
Definition: face_tri3.h:150
virtual bool is_face(const unsigned int i) const override
Definition: face_tri3.C:99
virtual Real volume() const override
Definition: face_tri3.C:194
long double min(long double a, double b)
A geometric point in (x,y,z) space.
Definition: point.h:38
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:1914
const Point & point(const unsigned int i) const
Definition: elem.h:1892
static const float _embedding_matrix[num_children][num_nodes][num_nodes]
Definition: face_tri3.h:205
std::unique_ptr< Elem > side(const unsigned int i) const
Definition: elem.h:2202
const Real pi
Definition: libmesh.h:233