fe_xyz_boundary.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 
20 // C++ includes
21 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
22 #include <cmath> // for std::sqrt
23 
24 
25 // Local includes
26 #include "libmesh/libmesh_common.h"
27 #include "libmesh/fe.h"
28 #include "libmesh/quadrature.h"
29 #include "libmesh/elem.h"
31 
32 namespace libMesh
33 {
34 
35 
36 
37 
38 //-------------------------------------------------------
39 // Full specialization for 0D when this is a useless method
40 template <>
41 void FEXYZ<0>::reinit(const Elem *,
42  const unsigned int,
43  const Real,
44  const std::vector<Point> * const,
45  const std::vector<Real> * const)
46 {
47  libmesh_error_msg("ERROR: This method only makes sense for 2D/3D elements!");
48 }
49 
50 
51 //-------------------------------------------------------
52 // Full specialization for 1D when this is a useless method
53 template <>
54 void FEXYZ<1>::reinit(const Elem *,
55  const unsigned int,
56  const Real,
57  const std::vector<Point> * const,
58  const std::vector<Real> * const)
59 {
60  libmesh_error_msg("ERROR: This method only makes sense for 2D/3D elements!");
61 }
62 
63 
64 //-------------------------------------------------------
65 // Method for 2D, 3D
66 template <unsigned int Dim>
67 void FEXYZ<Dim>::reinit(const Elem * elem,
68  const unsigned int s,
69  const Real,
70  const std::vector<Point> * const pts,
71  const std::vector<Real> * const weights)
72 {
73  libmesh_assert(elem);
74  libmesh_assert (this->qrule != nullptr || pts != nullptr);
75  // We don't do this for 1D elements!
76  libmesh_assert_not_equal_to (Dim, 1);
77 
78  // Build the side of interest
79  const std::unique_ptr<const Elem> side(elem->build_side_ptr(s));
80 
81  // Initialize the shape functions at the user-specified
82  // points
83  if (pts != nullptr)
84  {
85  // We can't get away without recomputing shape functions next
86  // time
87  this->shapes_on_quadrature = false;
88 
89  // Set the element type
90  this->elem_type = elem->type();
91 
92  // Initialize the face shape functions
93  this->_fe_map->template init_face_shape_functions<Dim>(*pts, side.get());
94  if (weights != nullptr)
95  {
96  this->compute_face_values (elem, side.get(), *weights);
97  }
98  else
99  {
100  std::vector<Real> dummy_weights (pts->size(), 1.);
101  // Compute data on the face for integration
102  this->compute_face_values (elem, side.get(), dummy_weights);
103  }
104  }
105  else
106  {
107  // initialize quadrature rule
108  this->qrule->init(side->type(), elem->p_level());
109 
110  {
111  // Set the element type
112  this->elem_type = elem->type();
113 
114  // Initialize the face shape functions
115  this->_fe_map->template init_face_shape_functions<Dim>(this->qrule->get_points(), side.get());
116  }
117  // We can't get away without recomputing shape functions next
118  // time
119  this->shapes_on_quadrature = false;
120  // Compute data on the face for integration
121  this->compute_face_values (elem, side.get(), this->qrule->get_weights());
122  }
123 }
124 
125 
126 
127 template <unsigned int Dim>
129  const Elem * side,
130  const std::vector<Real> & qw)
131 {
132  libmesh_assert(elem);
133  libmesh_assert(side);
134 
135  LOG_SCOPE("compute_face_values()", "FEXYZ");
136 
137  // The number of quadrature points.
138  const std::size_t n_qp = qw.size();
139 
140  // Number of shape functions in the finite element approximation
141  // space.
142  const unsigned int n_approx_shape_functions =
143  this->n_shape_functions(this->get_type(),
144  this->get_order());
145 
146  // Resize the shape functions and their gradients
147  this->phi.resize (n_approx_shape_functions);
148  this->dphi.resize (n_approx_shape_functions);
149  this->dphidx.resize (n_approx_shape_functions);
150  this->dphidy.resize (n_approx_shape_functions);
151  this->dphidz.resize (n_approx_shape_functions);
152 
153  for (unsigned int i=0; i<n_approx_shape_functions; i++)
154  {
155  this->phi[i].resize (n_qp);
156  this->dphi[i].resize (n_qp);
157  this->dphidx[i].resize (n_qp);
158  this->dphidy[i].resize (n_qp);
159  this->dphidz[i].resize (n_qp);
160  }
161 
162  this->_fe_map->compute_face_map(this->dim, qw, side);
163 
164  const std::vector<libMesh::Point> & xyz = this->_fe_map->get_xyz();
165 
166  switch (this->dim)
167  {
168  // A 2D finite element living in either 2D or 3D space.
169  // This means the boundary is a 1D finite element, i.e.
170  // and EDGE2 or EDGE3.
171  case 2:
172  {
173  // compute the shape function values & gradients
174  for (unsigned int i=0; i<n_approx_shape_functions; i++)
175  for (std::size_t p=0; p<n_qp; p++)
176  {
177  this->phi[i][p] = FE<Dim,XYZ>::shape (elem, this->fe_type.order, i, xyz[p]);
178 
179  this->dphi[i][p](0) =
180  this->dphidx[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 0, xyz[p]);
181 
182  this->dphi[i][p](1) =
183  this->dphidy[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 1, xyz[p]);
184 
185 #if LIBMESH_DIM == 3
186  this->dphi[i][p](2) = // can only assign to the Z component if LIBMESH_DIM==3
187 #endif
188  this->dphidz[i][p] = 0.;
189  }
190 
191  // done computing face values
192  break;
193  }
194 
195  // A 3D finite element living in 3D space.
196  case 3:
197  {
198  // compute the shape function values & gradients
199  for (unsigned int i=0; i<n_approx_shape_functions; i++)
200  for (std::size_t p=0; p<n_qp; p++)
201  {
202  this->phi[i][p] = FE<Dim,XYZ>::shape (elem, this->fe_type.order, i, xyz[p]);
203 
204  this->dphi[i][p](0) =
205  this->dphidx[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 0, xyz[p]);
206 
207  this->dphi[i][p](1) =
208  this->dphidy[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 1, xyz[p]);
209 
210  this->dphi[i][p](2) =
211  this->dphidz[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 2, xyz[p]);
212  }
213 
214  // done computing face values
215  break;
216  }
217 
218  default:
219  libmesh_error_msg("Invalid dim " << this->dim);
220  }
221 }
222 
223 
224 
225 
226 //--------------------------------------------------------------
227 // Explicit instantiations (doesn't make sense in 1D!) using fe_macro.h's macro
228 
229 template class FEXYZ<2>;
230 template class FEXYZ<3>;
231 
232 
233 } // namespace libMesh
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
void compute_face_values(const Elem *elem, const Elem *side, const std::vector< Real > &weights)
unsigned short int side
Definition: xdr_io.C:50
The base class for all geometric element types.
Definition: elem.h:100
unsigned int p_level() const
Definition: elem.h:2555
virtual std::unique_ptr< Elem > build_side_ptr(const unsigned int i, bool proxy=true)=0
Template class which generates the different FE families and orders.
Definition: fe.h:89
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual ElemType type() const =0
virtual void reinit(const Elem *elem, const std::vector< Point > *const pts=nullptr, const std::vector< Real > *const weights=nullptr) override
Definition: fe.h:842