fe_xyz_shape_1D.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 
19 // C++ includes
20 
21 // Local includes
22 #include "libmesh/fe.h"
23 #include "libmesh/elem.h"
24 
25 
26 
27 namespace libMesh
28 {
29 
30 
31 template <>
33  const Order,
34  const unsigned int,
35  const Point &)
36 {
37  libmesh_error_msg("XYZ polynomials require the element \n because the centroid is needed.");
38  return 0.;
39 }
40 
41 
42 
43 template <>
44 Real FE<1,XYZ>::shape(const Elem * elem,
45  const Order libmesh_dbg_var(order),
46  const unsigned int i,
47  const Point & point_in)
48 {
49  libmesh_assert(elem);
50  libmesh_assert_less_equal (i, order + elem->p_level());
51 
52  Point centroid = elem->centroid();
53  Real max_distance = 0.;
54  for (unsigned int p = 0; p < elem->n_nodes(); p++)
55  {
56  const Real distance = std::abs(centroid(0) - elem->point(p)(0));
57  max_distance = std::max(distance, max_distance);
58  }
59 
60  const Real x = point_in(0);
61  const Real xc = centroid(0);
62  const Real dx = (x - xc)/max_distance;
63 
64  // monomials. since they are hierarchic we only need one case block.
65  switch (i)
66  {
67  case 0:
68  return 1.;
69 
70  case 1:
71  return dx;
72 
73  case 2:
74  return dx*dx;
75 
76  case 3:
77  return dx*dx*dx;
78 
79  case 4:
80  return dx*dx*dx*dx;
81 
82  default:
83  Real val = 1.;
84  for (unsigned int index = 0; index != i; ++index)
85  val *= dx;
86  return val;
87  }
88 }
89 
90 
91 
92 template <>
94  const Order,
95  const unsigned int,
96  const unsigned int,
97  const Point &)
98 {
99  libmesh_error_msg("XYZ polynomials require the element \nbecause the centroid is needed.");
100  return 0.;
101 }
102 
103 
104 
105 template <>
107  const Order libmesh_dbg_var(order),
108  const unsigned int i,
109  const unsigned int libmesh_dbg_var(j),
110  const Point & point_in)
111 {
112  libmesh_assert(elem);
113  libmesh_assert_less_equal (i, order + elem->p_level());
114 
115  // only d()/dxi in 1D!
116 
117  libmesh_assert_equal_to (j, 0);
118 
119  Point centroid = elem->centroid();
120  Real max_distance = 0.;
121  for (unsigned int p = 0; p < elem->n_nodes(); p++)
122  {
123  const Real distance = std::abs(centroid(0) - elem->point(p)(0));
124  max_distance = std::max(distance, max_distance);
125  }
126 
127  const Real x = point_in(0);
128  const Real xc = centroid(0);
129  const Real dx = (x - xc)/max_distance;
130 
131  // monomials. since they are hierarchic we only need one case block.
132  switch (i)
133  {
134  case 0:
135  return 0.;
136 
137  case 1:
138  return 1.;
139 
140  case 2:
141  return 2.*dx/max_distance;
142 
143  case 3:
144  return 3.*dx*dx/max_distance;
145 
146  case 4:
147  return 4.*dx*dx*dx/max_distance;
148 
149  default:
150  Real val = i;
151  for (unsigned int index = 1; index != i; ++index)
152  val *= dx;
153  return val/max_distance;
154  }
155 }
156 
157 
158 
159 template <>
161  const Order,
162  const unsigned int,
163  const unsigned int,
164  const Point &)
165 {
166  libmesh_error_msg("XYZ polynomials require the element \nbecause the centroid is needed.");
167  return 0.;
168 }
169 
170 
171 
172 template <>
174  const Order libmesh_dbg_var(order),
175  const unsigned int i,
176  const unsigned int libmesh_dbg_var(j),
177  const Point & point_in)
178 {
179  libmesh_assert(elem);
180  libmesh_assert_less_equal (i, order + elem->p_level());
181 
182  // only d2()/dxi2 in 1D!
183 
184  libmesh_assert_equal_to (j, 0);
185 
186  Point centroid = elem->centroid();
187  Real max_distance = 0.;
188  for (unsigned int p = 0; p < elem->n_nodes(); p++)
189  {
190  const Real distance = std::abs(centroid(0) - elem->point(p)(0));
191  max_distance = std::max(distance, max_distance);
192  }
193 
194  const Real x = point_in(0);
195  const Real xc = centroid(0);
196  const Real dx = (x - xc)/max_distance;
197  const Real dist2 = pow(max_distance,2.);
198 
199  // monomials. since they are hierarchic we only need one case block.
200  switch (i)
201  {
202  case 0:
203  case 1:
204  return 0.;
205 
206  case 2:
207  return 2./dist2;
208 
209  case 3:
210  return 6.*dx/dist2;
211 
212  case 4:
213  return 12.*dx*dx/dist2;
214 
215  default:
216  Real val = 2.;
217  for (unsigned int index = 2; index != i; ++index)
218  val *= (index+1) * dx;
219  return val/dist2;
220  }
221 }
222 
223 } // namespace libMesh
double abs(double a)
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
The base class for all geometric element types.
Definition: elem.h:100
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
unsigned int p_level() const
Definition: elem.h:2555
long double max(long double a, double b)
virtual unsigned int n_nodes() const =0
double pow(double a, int b)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
A geometric point in (x,y,z) space.
Definition: point.h:38
const Point & point(const unsigned int i) const
Definition: elem.h:1892
virtual Point centroid() const
Definition: elem.C:344
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)