libMesh
fe_xyz_shape_1D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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  libmesh_error_msg("We'll never get here!");
90  return 0.;
91 
92 }
93 
94 
95 
96 template <>
98  const Order,
99  const unsigned int,
100  const unsigned int,
101  const Point &)
102 {
103  libmesh_error_msg("XYZ polynomials require the element \nbecause the centroid is needed.");
104  return 0.;
105 }
106 
107 
108 
109 template <>
111  const Order libmesh_dbg_var(order),
112  const unsigned int i,
113  const unsigned int libmesh_dbg_var(j),
114  const Point & point_in)
115 {
116  libmesh_assert(elem);
117  libmesh_assert_less_equal (i, order + elem->p_level());
118 
119  // only d()/dxi in 1D!
120 
121  libmesh_assert_equal_to (j, 0);
122 
123  Point centroid = elem->centroid();
124  Real max_distance = 0.;
125  for (unsigned int p = 0; p < elem->n_nodes(); p++)
126  {
127  const Real distance = std::abs(centroid(0) - elem->point(p)(0));
128  max_distance = std::max(distance, max_distance);
129  }
130 
131  const Real x = point_in(0);
132  const Real xc = centroid(0);
133  const Real dx = (x - xc)/max_distance;
134 
135  // monomials. since they are hierarchic we only need one case block.
136  switch (i)
137  {
138  case 0:
139  return 0.;
140 
141  case 1:
142  return 1.;
143 
144  case 2:
145  return 2.*dx/max_distance;
146 
147  case 3:
148  return 3.*dx*dx/max_distance;
149 
150  case 4:
151  return 4.*dx*dx*dx/max_distance;
152 
153  default:
154  Real val = i;
155  for (unsigned int index = 1; index != i; ++index)
156  val *= dx;
157  return val/max_distance;
158  }
159 
160  libmesh_error_msg("We'll never get here!");
161  return 0.;
162 }
163 
164 
165 
166 template <>
168  const Order,
169  const unsigned int,
170  const unsigned int,
171  const Point &)
172 {
173  libmesh_error_msg("XYZ polynomials require the element \nbecause the centroid is needed.");
174  return 0.;
175 }
176 
177 
178 
179 template <>
181  const Order libmesh_dbg_var(order),
182  const unsigned int i,
183  const unsigned int libmesh_dbg_var(j),
184  const Point & point_in)
185 {
186  libmesh_assert(elem);
187  libmesh_assert_less_equal (i, order + elem->p_level());
188 
189  // only d2()/dxi2 in 1D!
190 
191  libmesh_assert_equal_to (j, 0);
192 
193  Point centroid = elem->centroid();
194  Real max_distance = 0.;
195  for (unsigned int p = 0; p < elem->n_nodes(); p++)
196  {
197  const Real distance = std::abs(centroid(0) - elem->point(p)(0));
198  max_distance = std::max(distance, max_distance);
199  }
200 
201  const Real x = point_in(0);
202  const Real xc = centroid(0);
203  const Real dx = (x - xc)/max_distance;
204  const Real dist2 = pow(max_distance,2.);
205 
206  // monomials. since they are hierarchic we only need one case block.
207  switch (i)
208  {
209  case 0:
210  case 1:
211  return 0.;
212 
213  case 2:
214  return 2./dist2;
215 
216  case 3:
217  return 6.*dx/dist2;
218 
219  case 4:
220  return 12.*dx*dx/dist2;
221 
222  default:
223  Real val = 2.;
224  for (unsigned int index = 2; index != i; ++index)
225  val *= (index+1) * dx;
226  return val/dist2;
227  }
228 
229  libmesh_error_msg("We'll never get here!");
230  return 0.;
231 }
232 
233 } // namespace libMesh
double abs(double a)
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
unsigned int p_level() const
Definition: elem.h:2422
ElemType
Defines an enum for geometric element types.
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
The libMesh namespace provides an interface to certain functionality in the library.
long double max(long double a, double b)
Real distance(const Point &p)
libmesh_assert(j)
virtual unsigned int n_nodes() const =0
PetscErrorCode Vec x
PetscErrorCode Vec Mat libmesh_dbg_var(j)
double pow(double a, int b)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Point & point(const unsigned int i) const
Definition: elem.h:1809
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
virtual Point centroid() const
Definition: elem.C:446
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)