fe_xyz_boundary.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 
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::err << "ERROR: This method only makes sense for 2D/3D elements!"
48  << std::endl;
49  libmesh_error();
50 }
51 
52 
53 //-------------------------------------------------------
54 // Full specialization for 1D when this is a useless method
55 template <>
56 void FEXYZ<1>::reinit(const Elem*,
57  const unsigned int,
58  const Real,
59  const std::vector<Point>* const,
60  const std::vector<Real>* const)
61 {
62  libMesh::err << "ERROR: This method only makes sense for 2D/3D elements!"
63  << std::endl;
64  libmesh_error();
65 }
66 
67 
68 //-------------------------------------------------------
69 // Method for 2D, 3D
70 template <unsigned int Dim>
71 void FEXYZ<Dim>::reinit(const Elem* elem,
72  const unsigned int s,
73  const Real,
74  const std::vector<Point>* const pts,
75  const std::vector<Real>* const weights)
76 {
77  libmesh_assert(elem);
78  libmesh_assert (this->qrule != NULL || pts != NULL);
79  // We don't do this for 1D elements!
80  libmesh_assert_not_equal_to (Dim, 1);
81 
82  // Build the side of interest
83  const AutoPtr<Elem> side(elem->build_side(s));
84 
85  // Initialize the shape functions at the user-specified
86  // points
87  if (pts != NULL)
88  {
89  // We can't get away without recomputing shape functions next
90  // time
91  this->shapes_on_quadrature = false;
92 
93  // Set the element type
94  this->elem_type = elem->type();
95 
96  // Initialize the face shape functions
97  this->_fe_map->template init_face_shape_functions<Dim>(*pts, side.get());
98  if (weights != NULL)
99  {
100  this->compute_face_values (elem, side.get(), *weights);
101  }
102  else
103  {
104  std::vector<Real> dummy_weights (pts->size(), 1.);
105  // Compute data on the face for integration
106  this->compute_face_values (elem, side.get(), dummy_weights);
107  }
108  }
109  else
110  {
111  // initialize quadrature rule
112  this->qrule->init(side->type(), elem->p_level());
113 
114  {
115  // Set the element type
116  this->elem_type = elem->type();
117 
118  // Initialize the face shape functions
119  this->_fe_map->template init_face_shape_functions<Dim>(this->qrule->get_points(), side.get());
120  }
121  // We can't get away without recomputing shape functions next
122  // time
123  this->shapes_on_quadrature = false;
124  // Compute data on the face for integration
125  this->compute_face_values (elem, side.get(), this->qrule->get_weights());
126  }
127 }
128 
129 
130 
131 template <unsigned int Dim>
133  const Elem* side,
134  const std::vector<Real>& qw)
135 {
136  libmesh_assert(elem);
137  libmesh_assert(side);
138 
139  START_LOG("compute_face_values()", "FEXYZ");
140 
141  // The number of quadrature points.
142  const std::size_t n_qp = qw.size();
143 
144  // Number of shape functions in the finite element approximation
145  // space.
146  const unsigned int n_approx_shape_functions =
147  this->n_shape_functions(this->get_type(),
148  this->get_order());
149 
150  // Resize the shape functions and their gradients
151  this->phi.resize (n_approx_shape_functions);
152  this->dphi.resize (n_approx_shape_functions);
153  this->dphidx.resize (n_approx_shape_functions);
154  this->dphidy.resize (n_approx_shape_functions);
155  this->dphidz.resize (n_approx_shape_functions);
156 
157  for (unsigned int i=0; i<n_approx_shape_functions; i++)
158  {
159  this->phi[i].resize (n_qp);
160  this->dphi[i].resize (n_qp);
161  this->dphidx[i].resize (n_qp);
162  this->dphidy[i].resize (n_qp);
163  this->dphidz[i].resize (n_qp);
164  }
165 
166  this->_fe_map->compute_face_map(this->dim, qw, side);
167 
168  const std::vector<libMesh::Point>& xyz = this->_fe_map->get_xyz();
169 
170  switch (this->dim)
171  {
172  // A 2D finite element living in either 2D or 3D space.
173  // This means the boundary is a 1D finite element, i.e.
174  // and EDGE2 or EDGE3.
175  case 2:
176  {
177  // compute the shape function values & gradients
178  for (unsigned int i=0; i<n_approx_shape_functions; i++)
179  for (std::size_t p=0; p<n_qp; p++)
180  {
181  this->phi[i][p] = FE<Dim,XYZ>::shape (elem, this->fe_type.order, i, xyz[p]);
182 
183  this->dphi[i][p](0) =
184  this->dphidx[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 0, xyz[p]);
185 
186  this->dphi[i][p](1) =
187  this->dphidy[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 1, xyz[p]);
188 
189 #if LIBMESH_DIM == 3
190  this->dphi[i][p](2) = // can only assign to the Z component if LIBMESH_DIM==3
191 #endif
192  this->dphidz[i][p] = 0.;
193  }
194 
195  // done computing face values
196  break;
197  }
198 
199  // A 3D finite element living in 3D space.
200  case 3:
201  {
202  // compute the shape function values & gradients
203  for (unsigned int i=0; i<n_approx_shape_functions; i++)
204  for (std::size_t p=0; p<n_qp; p++)
205  {
206  this->phi[i][p] = FE<Dim,XYZ>::shape (elem, this->fe_type.order, i, xyz[p]);
207 
208  this->dphi[i][p](0) =
209  this->dphidx[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 0, xyz[p]);
210 
211  this->dphi[i][p](1) =
212  this->dphidy[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 1, xyz[p]);
213 
214  this->dphi[i][p](2) =
215  this->dphidz[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 2, xyz[p]);
216  }
217 
218  // done computing face values
219  break;
220  }
221 
222 
223  default:
224  libmesh_error();
225 
226  }
227 
228  STOP_LOG("compute_face_values()", "FEXYZ");
229 }
230 
231 
232 
233 
234 //--------------------------------------------------------------
235 // Explicit instantiations (doesn't make sense in 1D!) using fe_macro.h's macro
236 
237 template class FEXYZ<2>;
238 template class FEXYZ<3>;
239 
240 
241 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo