fe_xyz_shape_1D.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 
19 // C++ inlcludes
20 
21 // Local includes
22 #include "libmesh/fe.h"
23 #include "libmesh/elem.h"
24 
25 
26 
27 
28 // Anonymous namespace for persistant variables.
29 // This allows us to determine when the centroid needs
30 // to be recalculated.
31 namespace
32 {
33  using namespace libMesh;
34 
35  static dof_id_type old_elem_id = DofObject::invalid_id;
36  static libMesh::Point centroid;
37  static Real max_distance;
38 }
39 
40 
41 namespace libMesh
42 {
43 
44 
45 template <>
47  const Order,
48  const unsigned int,
49  const Point&)
50 {
51  libMesh::err << "XYZ polynomials require the element\n"
52  << "because the centroid is needed."
53  << std::endl;
54 
55  libmesh_error();
56  return 0.;
57 }
58 
59 
60 
61 template <>
63  const Order libmesh_dbg_var(order),
64  const unsigned int i,
65  const Point& point_in)
66 {
67  libmesh_assert(elem);
68  libmesh_assert_less_equal (i, order + elem->p_level());
69 
70  // Only recompute the centroid if the element
71  // has changed from the last one we computed.
72  // This avoids repeated centroid calculations
73  // when called in succession with the same element.
74  if (elem->id() != old_elem_id)
75  {
76  centroid = elem->centroid();
77  old_elem_id = elem->id();
78  max_distance = 0.;
79  for (unsigned int p = 0; p < elem->n_nodes(); p++)
80  {
81  const Real distance = std::abs(centroid(0) - elem->point(p)(0));
82  max_distance = std::max(distance, max_distance);
83  }
84  }
85 
86  // Using static globals for old_elem_id, etc. will fail
87  // horribly with more than one thread.
88  libmesh_assert_equal_to (libMesh::n_threads(), 1);
89 
90  const Real x = point_in(0);
91  const Real xc = centroid(0);
92  const Real dx = (x - xc)/max_distance;
93 
94  // monomials. since they are hierarchic we only need one case block.
95  switch (i)
96  {
97  case 0:
98  return 1.;
99 
100  case 1:
101  return dx;
102 
103  case 2:
104  return dx*dx;
105 
106  case 3:
107  return dx*dx*dx;
108 
109  case 4:
110  return dx*dx*dx*dx;
111 
112  default:
113  Real val = 1.;
114  for (unsigned int index = 0; index != i; ++index)
115  val *= dx;
116  return val;
117  }
118 
119  libmesh_error();
120  return 0.;
121 
122 }
123 
124 
125 
126 template <>
128  const Order,
129  const unsigned int,
130  const unsigned int,
131  const Point&)
132 {
133  libMesh::err << "XYZ polynomials require the element\n"
134  << "because the centroid is needed."
135  << std::endl;
136 
137  libmesh_error();
138  return 0.;
139 }
140 
141 
142 
143 template <>
145  const Order libmesh_dbg_var(order),
146  const unsigned int i,
147  const unsigned int libmesh_dbg_var(j),
148  const Point& point_in)
149 {
150  libmesh_assert(elem);
151  libmesh_assert_less_equal (i, order + elem->p_level());
152 
153  // only d()/dxi in 1D!
154 
155  libmesh_assert_equal_to (j, 0);
156 
157  // Only recompute the centroid if the element
158  // has changed from the last one we computed.
159  // This avoids repeated centroid calculations
160  // when called in succession with the same element.
161  if (elem->id() != old_elem_id)
162  {
163  centroid = elem->centroid();
164  old_elem_id = elem->id();
165  max_distance = 0.;
166  for (unsigned int p = 0; p < elem->n_nodes(); p++)
167  {
168  const Real distance = std::abs(centroid(0) - elem->point(p)(0));
169  max_distance = std::max(distance, max_distance);
170  }
171  }
172 
173  // Using static globals for old_elem_id, etc. will fail
174  // horribly with more than one thread.
175  libmesh_assert_equal_to (libMesh::n_threads(), 1);
176 
177  const Real x = point_in(0);
178  const Real xc = centroid(0);
179  const Real dx = (x - xc)/max_distance;
180 
181  // monomials. since they are hierarchic we only need one case block.
182  switch (i)
183  {
184  case 0:
185  return 0.;
186 
187  case 1:
188  return 1.;
189 
190  case 2:
191  return 2.*dx/max_distance;
192 
193  case 3:
194  return 3.*dx*dx/max_distance;
195 
196  case 4:
197  return 4.*dx*dx*dx/max_distance;
198 
199  default:
200  Real val = i;
201  for (unsigned int index = 1; index != i; ++index)
202  val *= dx;
203  return val/max_distance;
204  }
205 
206  libmesh_error();
207  return 0.;
208 }
209 
210 
211 
212 template <>
214  const Order,
215  const unsigned int,
216  const unsigned int,
217  const Point&)
218 {
219  libMesh::err << "XYZ polynomials require the element\n"
220  << "because the centroid is needed."
221  << std::endl;
222 
223  libmesh_error();
224  return 0.;
225 }
226 
227 
228 
229 template <>
231  const Order libmesh_dbg_var(order),
232  const unsigned int i,
233  const unsigned int libmesh_dbg_var(j),
234  const Point& point_in)
235 {
236  libmesh_assert(elem);
237  libmesh_assert_less_equal (i, order + elem->p_level());
238 
239  // only d2()/dxi2 in 1D!
240 
241  libmesh_assert_equal_to (j, 0);
242 
243  // Only recompute the centroid if the element
244  // has changed from the last one we computed.
245  // This avoids repeated centroid calculations
246  // when called in succession with the same element.
247  if (elem->id() != old_elem_id)
248  {
249  centroid = elem->centroid();
250  old_elem_id = elem->id();
251  max_distance = 0.;
252  for (unsigned int p = 0; p < elem->n_nodes(); p++)
253  {
254  const Real distance = std::abs(centroid(0) - elem->point(p)(0));
255  max_distance = std::max(distance, max_distance);
256  }
257  }
258 
259  // Using static globals for old_elem_id, etc. will fail
260  // horribly with more than one thread.
261  libmesh_assert_equal_to (libMesh::n_threads(), 1);
262 
263  const Real x = point_in(0);
264  const Real xc = centroid(0);
265  const Real dx = (x - xc)/max_distance;
266  const Real dist2 = pow(max_distance,2.);
267 
268  // monomials. since they are hierarchic we only need one case block.
269  switch (i)
270  {
271  case 0:
272  case 1:
273  return 0.;
274 
275  case 2:
276  return 2./dist2;
277 
278  case 3:
279  return 6.*dx/dist2;
280 
281  case 4:
282  return 12.*dx*dx/dist2;
283 
284  default:
285  Real val = 2.;
286  for (unsigned int index = 2; index != i; ++index)
287  val *= (index+1) * dx;
288  return val/dist2;
289  }
290 
291  libmesh_error();
292  return 0.;
293 }
294 
295 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo