face_tri.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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/face_tri.h"
22 #include "libmesh/edge_edge2.h"
23 
24 namespace libMesh
25 {
26 
27 
28 
29 
30 
31 // ------------------------------------------------------------
32 // Tri class member functions
33 dof_id_type Tri::key (const unsigned int s) const
34 {
35  libmesh_assert_less (s, this->n_sides());
36 
37  switch (s)
38  {
39  case 0:
40  return
41  this->compute_key (this->node(0),
42  this->node(1));
43 
44  case 1:
45  return
46  this->compute_key (this->node(1),
47  this->node(2));
48  case 2:
49  return
50  this->compute_key (this->node(2),
51  this->node(0));
52  }
53 
54 
55  // We will never get here... Look at the code above.
56  libmesh_error();
57  return 0;
58 }
59 
60 
61 
62 AutoPtr<Elem> Tri::side (const unsigned int i) const
63 {
64  libmesh_assert_less (i, this->n_sides());
65 
66  Elem* edge = new Edge2;
67 
68  switch (i)
69  {
70  case 0:
71  {
72  edge->set_node(0) = this->get_node(0);
73  edge->set_node(1) = this->get_node(1);
74 
75  AutoPtr<Elem> ap_edge(edge);
76  return ap_edge;
77  }
78  case 1:
79  {
80  edge->set_node(0) = this->get_node(1);
81  edge->set_node(1) = this->get_node(2);
82 
83  AutoPtr<Elem> ap_edge(edge);
84  return ap_edge;
85  }
86  case 2:
87  {
88  edge->set_node(0) = this->get_node(2);
89  edge->set_node(1) = this->get_node(0);
90 
91  AutoPtr<Elem> ap_edge(edge);
92  return ap_edge;
93  }
94  default:
95  {
96  libmesh_error();
97  }
98  }
99 
100 
101  // We will never get here... Look at the code above.
102  libmesh_error();
103  AutoPtr<Elem> ap_edge(edge);
104  return ap_edge;
105 }
106 
107 
108 
109 bool Tri::is_child_on_side(const unsigned int c,
110  const unsigned int s) const
111 {
112  libmesh_assert_less (c, this->n_children());
113  libmesh_assert_less (s, this->n_sides());
114 
115  return (c == s || c == (s+1)%3);
116 }
117 
118 
119 
121 {
122  switch (q)
123  {
124 
128  case DISTORTION:
129  case STRETCH:
130  {
131  const Node* p1 = this->get_node(0);
132  const Node* p2 = this->get_node(1);
133  const Node* p3 = this->get_node(2);
134 
135  Point v1 = (*p2) - (*p1);
136  Point v2 = (*p3) - (*p1);
137  Point v3 = (*p3) - (*p2);
138  const Real l1 = v1.size();
139  const Real l2 = v2.size();
140  const Real l3 = v3.size();
141 
142  // if one length is 0, quality is quite bad!
143  if ((l1 <=0.) || (l2 <= 0.) || (l3 <= 0.))
144  return 0.;
145 
146  const Real s1 = std::sin(std::acos(v1*v2/l1/l2)/2.);
147  v1 *= -1;
148  const Real s2 = std::sin(std::acos(v1*v3/l1/l3)/2.);
149  const Real s3 = std::sin(std::acos(v2*v3/l2/l3)/2.);
150 
151  return 8. * s1 * s2 * s3;
152 
153  }
154  default:
155  return Elem::quality(q);
156  }
157 
163  return Elem::quality(q);
164 
165 }
166 
167 
168 
169 
170 
171 
172 std::pair<Real, Real> Tri::qual_bounds (const ElemQuality q) const
173 {
174  std::pair<Real, Real> bounds;
175 
176  switch (q)
177  {
178 
179  case MAX_ANGLE:
180  bounds.first = 60.;
181  bounds.second = 90.;
182  break;
183 
184  case MIN_ANGLE:
185  bounds.first = 30.;
186  bounds.second = 60.;
187  break;
188 
189  case CONDITION:
190  bounds.first = 1.;
191  bounds.second = 1.3;
192  break;
193 
194  case JACOBIAN:
195  bounds.first = 0.5;
196  bounds.second = 1.155;
197  break;
198 
199  case SIZE:
200  case SHAPE:
201  bounds.first = 0.25;
202  bounds.second = 1.;
203  break;
204 
205  case DISTORTION:
206  bounds.first = 0.6;
207  bounds.second = 1.;
208  break;
209 
210  default:
211  libMesh::out << "Warning: Invalid quality measure chosen." << std::endl;
212  bounds.first = -1;
213  bounds.second = -1;
214  }
215 
216  return bounds;
217 }
218 
219 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:04 UTC

Hosted By:
SourceForge.net Logo