1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2020 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
19
20 // Local includes
21 #include "libmesh/edge_edge3.h"
22 #include "libmesh/enum_io_package.h"
23 #include "libmesh/enum_order.h"
24
25 namespace libMesh
26 {
27
28 // Edge3 class static member initializations
29 const int Edge3::num_nodes;
30 const int Edge3::num_sides;
31 const int Edge3::num_edges;
32 const int Edge3::num_children;
33 const int Edge3::nodes_per_side;
34 const int Edge3::nodes_per_edge;
35
36 #ifdef LIBMESH_ENABLE_AMR
37
38 const float Edge3::_embedding_matrix[2][3][3] =
39 {
40 // embedding matrix for child 0
41 {
42 // 0 1 2
43 {1.0, 0.0, 0.0}, // left
44 {0.0, 0.0, 1.0}, // right
45 {0.375,-0.125,0.75} // middle
46 },
47
48 // embedding matrix for child 1
49 {
50 // 0 1 2
51 {0.0, 0.0, 1.0}, // left
52 {0.0, 1.0, 0.0}, // right
53 {-0.125,0.375,0.75} // middle
54 }
55 };
56
57 #endif
58
is_vertex(const unsigned int i)59 bool Edge3::is_vertex(const unsigned int i) const
60 {
61 if (i < 2)
62 return true;
63 return false;
64 }
65
is_edge(const unsigned int i)66 bool Edge3::is_edge(const unsigned int i) const
67 {
68 if (i < 2)
69 return false;
70 return true;
71 }
72
is_face(const unsigned int)73 bool Edge3::is_face(const unsigned int ) const
74 {
75 return false;
76 }
77
is_node_on_side(const unsigned int n,const unsigned int s)78 bool Edge3::is_node_on_side(const unsigned int n,
79 const unsigned int s) const
80 {
81 libmesh_assert_less (s, 2);
82 libmesh_assert_less (n, Edge3::num_nodes);
83 return (s == n);
84 }
85
is_node_on_edge(const unsigned int,const unsigned int libmesh_dbg_var (e))86 bool Edge3::is_node_on_edge(const unsigned int,
87 const unsigned int libmesh_dbg_var(e)) const
88 {
89 libmesh_assert_equal_to (e, 0);
90 return true;
91 }
92
93
94
has_affine_map()95 bool Edge3::has_affine_map() const
96 {
97 return (this->point(2).relative_fuzzy_equals
98 ((this->point(0) + this->point(1))/2));
99 }
100
101
102
default_order()103 Order Edge3::default_order() const
104 {
105 return SECOND;
106 }
107
108
109
connectivity(const unsigned int sc,const IOPackage iop,std::vector<dof_id_type> & conn)110 void Edge3::connectivity(const unsigned int sc,
111 const IOPackage iop,
112 std::vector<dof_id_type> & conn) const
113 {
114 libmesh_assert_less_equal (sc, 1);
115 libmesh_assert_less (sc, this->n_sub_elem());
116 libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
117
118 // Create storage
119 conn.resize(2);
120
121 switch (iop)
122 {
123 case TECPLOT:
124 {
125 switch (sc)
126 {
127 case 0:
128 conn[0] = this->node_id(0)+1;
129 conn[1] = this->node_id(2)+1;
130 return;
131
132 case 1:
133 conn[0] = this->node_id(2)+1;
134 conn[1] = this->node_id(1)+1;
135 return;
136
137 default:
138 libmesh_error_msg("Invalid sc = " << sc);
139 }
140 }
141
142
143 case VTK:
144 {
145 conn.resize(3);
146 conn[0] = this->node_id(0);
147 conn[1] = this->node_id(1);
148 conn[2] = this->node_id(2);
149 return;
150
151 /*
152 switch (sc)
153 {
154 case 0:
155 conn[0] = this->node_id(0);
156 conn[1] = this->node_id(2);
157
158 return;
159
160 case 1:
161 conn[0] = this->node_id(2);
162 conn[1] = this->node_id(1);
163
164 return;
165
166 default:
167 libmesh_error_msg("Invalid sc = " << sc);
168 }
169 */
170 }
171
172 default:
173 libmesh_error_msg("Unsupported IO package " << iop);
174 }
175 }
176
177
178
179 std::pair<unsigned short int, unsigned short int>
second_order_child_vertex(const unsigned int)180 Edge3::second_order_child_vertex (const unsigned int) const
181 {
182 return std::pair<unsigned short int, unsigned short int>(0,0);
183 }
184
185
186
volume()187 Real Edge3::volume () const
188 {
189 // Finding the (exact) length of a general quadratic element
190 // is a surprisingly complicated formula.
191 Point A = this->point(0) + this->point(1) - 2*this->point(2);
192 Point B = (this->point(1) - this->point(0))/2;
193
194 const Real a = A.norm_sq();
195 const Real c = B.norm_sq();
196
197 // Degenerate straight line case
198 if (a == 0.)
199 return 2. * std::sqrt(c);
200
201 const Real b = -2.*std::abs(A*B);
202
203 const Real sqrt_term1 = a - b + c;
204 const Real sqrt_term2 = a + b + c;
205
206 // Fall back on straight line case instead of computing nan
207 // Note: b can be positive or negative so we have to check both cases.
208 if (sqrt_term1 < 0. || sqrt_term2 < 0.)
209 return 2. * std::sqrt(c);
210
211 const Real r1 = std::sqrt(sqrt_term1);
212 const Real r2 = std::sqrt(sqrt_term2);
213 const Real rsum = r1 + r2;
214 const Real root_a = std::sqrt(a);
215 const Real b_over_root_a = b / root_a;
216
217 // Pre-compute the denominator of the log term. If it's zero, fall
218 // back on the straight line case.
219 const Real log_term_denom = -root_a - 0.5*b_over_root_a + r2;
220 if (log_term_denom == 0. || rsum == 0.)
221 return 2. * std::sqrt(c);
222
223 Real log1p_arg = 2*(root_a - b/rsum) / log_term_denom;
224
225 return
226 0.5*(rsum + b_over_root_a*b_over_root_a/rsum +
227 (c-0.25*b_over_root_a*b_over_root_a)*std::log1p(log1p_arg)/root_a);
228 }
229
230
231
loose_bounding_box()232 BoundingBox Edge3::loose_bounding_box () const
233 {
234 // This might be a curved line through 2-space or 3-space, in which
235 // case the full bounding box can be larger than the bounding box of
236 // just the nodes.
237 Point pmin, pmax;
238
239 for (unsigned d=0; d<LIBMESH_DIM; ++d)
240 {
241 Real center = this->point(2)(d);
242 Real hd = std::max(std::abs(center - this->point(0)(d)),
243 std::abs(center - this->point(1)(d)));
244
245 pmin(d) = center - hd;
246 pmax(d) = center + hd;
247 }
248
249 return BoundingBox(pmin, pmax);
250 }
251
252
253
key()254 dof_id_type Edge3::key () const
255 {
256 return this->compute_key(this->node_id(2));
257 }
258
259
260 } // namespace libMesh
261