libMesh
inf_fe_boundary.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 
20 // C++ includes
21 
22 
23 // Local includes
24 #include "libmesh/libmesh_config.h"
25 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
26 #include "libmesh/inf_fe.h"
27 #include "libmesh/inf_fe_macro.h"
28 #include "libmesh/quadrature.h"
29 #include "libmesh/elem.h"
30 
31 namespace libMesh
32 {
33 
34 // Method for 2D, 3D -- see inf_fe_1D.C for a 1D version of this
35 template <unsigned int Dim, FEFamily T_radial, InfMapType T_base>
37  const unsigned int s,
38  const Real tolerance,
39  const std::vector<Point> * const pts,
40  const std::vector<Real> * const weights)
41 {
42  if (weights != libmesh_nullptr)
43  libmesh_not_implemented_msg("ERROR: User-specified weights for infinite elements are not implemented!");
44 
45  if (pts != libmesh_nullptr)
46  libmesh_not_implemented_msg("ERROR: User-specified points for infinite elements are not implemented!");
47 
48  // We don't do this for 1D elements!
49  libmesh_assert_not_equal_to (Dim, 1);
50 
51  libmesh_assert(inf_elem);
52  libmesh_assert(qrule);
53 
54  // Don't do this for the base
55  libmesh_assert_not_equal_to (s, 0);
56 
57  // Build the side of interest
58  const UniquePtr<const Elem> side(inf_elem->build_side_ptr(s));
59 
60  // set the element type
61  elem_type = inf_elem->type();
62 
63  // eventually initialize radial quadrature rule
64  bool radial_qrule_initialized = false;
65 
66  if (current_fe_type.radial_order != fe_type.radial_order)
67  {
68  current_fe_type.radial_order = fe_type.radial_order;
69  radial_qrule->init(EDGE2, inf_elem->p_level());
70  radial_qrule_initialized = true;
71  }
72 
73  // Initialize the face shape functions
74  if (this->get_type() != inf_elem->type() ||
75  base_fe->shapes_need_reinit() ||
76  radial_qrule_initialized)
77  this->init_face_shape_functions (qrule->get_points(), side.get());
78 
79 
80  // compute the face map
81  this->_fe_map->compute_face_map(this->dim, _total_qrule_weights, side.get());
82 
83  // make a copy of the Jacobian for integration
84  const std::vector<Real> JxW_int(this->_fe_map->get_JxW());
85 
86  // Find where the integration points are located on the
87  // full element.
88  std::vector<Point> qp; this->inverse_map (inf_elem, this->_fe_map->get_xyz(),
89  qp, tolerance);
90 
91  // compute the shape function and derivative values
92  // at the points qp
93  this->reinit (inf_elem, &qp);
94 
95  // copy back old data
96  this->_fe_map->get_JxW() = JxW_int;
97 
98 }
99 
100 
101 
102 // Method for 2D, 3D -- see inf_fe_1D.C for a 1D version of this
103 template <unsigned int Dim, FEFamily T_radial, InfMapType T_base>
105  const unsigned int,
106  const Real,
107  const std::vector<Point> * const pts,
108  const std::vector<Real> * const /*weights*/)
109 {
110  // We don't do this for 1D elements!
111  //libmesh_assert_not_equal_to (Dim, 1);
112  libmesh_not_implemented_msg("ERROR: Edge conditions for infinite elements not implemented!");
113 
114  if (pts != libmesh_nullptr)
115  libmesh_not_implemented_msg("ERROR: User-specified points for infinite elements not implemented!");
116 }
117 
118 
119 
120 
121 template <unsigned int Dim, FEFamily T_radial, InfMapType T_base>
123  const Elem * inf_side)
124 {
125  libmesh_assert(inf_side);
126 
127  // Currently, this makes only sense in 3-D!
128  libmesh_assert_equal_to (Dim, 3);
129 
130  // Initialize the radial shape functions
131  this->init_radial_shape_functions(inf_side);
132 
133  // Initialize the base shape functions
134  this->update_base_elem(inf_side);
135 
136  // Initialize the base quadrature rule
137  base_qrule->init(base_elem->type(), inf_side->p_level());
138 
139  // base_fe still corresponds to the (dim-1)-dimensional base of the InfFE object,
140  // so update the fe_base.
141  {
142  libmesh_assert_equal_to (Dim, 3);
143  base_fe.reset(FEBase::build(Dim-2, this->fe_type).release());
144  base_fe->attach_quadrature_rule(base_qrule.get());
145  }
146 
147  // initialize the shape functions on the base
148  base_fe->init_base_shape_functions(base_fe->qrule->get_points(),
149  base_elem.get());
150 
151  // the number of quadrature points
152  const unsigned int n_radial_qp =
153  cast_int<unsigned int>(som.size());
154  const unsigned int n_base_qp = base_qrule->n_points();
155  const unsigned int n_total_qp = n_radial_qp * n_base_qp;
156 
157  // the quadrature weights
158  _total_qrule_weights.resize(n_total_qp);
159 
160  // now init the shapes for boundary work
161  {
162 
163  // The element type and order to use in the base map
164  const Order base_mapping_order ( base_elem->default_order() );
165  const ElemType base_mapping_elem_type ( base_elem->type() );
166 
167  // the number of mapping shape functions
168  // (Lagrange shape functions are used for mapping in the base)
169  const unsigned int n_radial_mapping_sf =
170  cast_int<unsigned int>(radial_map.size());
171  const unsigned int n_base_mapping_shape_functions = Base::n_base_mapping_sf(base_mapping_elem_type,
172  base_mapping_order);
173 
174  const unsigned int n_total_mapping_shape_functions =
175  n_radial_mapping_sf * n_base_mapping_shape_functions;
176 
177 
178  // initialize the node and shape numbering maps
179  {
180  _radial_node_index.resize (n_total_mapping_shape_functions);
181  _base_node_index.resize (n_total_mapping_shape_functions);
182 
183  const ElemType inf_face_elem_type (inf_side->type());
184 
185  // fill the node index map
186  for (unsigned int n=0; n<n_total_mapping_shape_functions; n++)
187  {
188  compute_node_indices (inf_face_elem_type,
189  n,
190  _base_node_index[n],
191  _radial_node_index[n]);
192 
193  libmesh_assert_less (_base_node_index[n], n_base_mapping_shape_functions);
194  libmesh_assert_less (_radial_node_index[n], n_radial_mapping_sf);
195  }
196 
197  }
198 
199  // resize map data fields
200  {
201  std::vector<std::vector<Real>> & psi_map = this->_fe_map->get_psi();
202  std::vector<std::vector<Real>> & dpsidxi_map = this->_fe_map->get_dpsidxi();
203  std::vector<std::vector<Real>> & d2psidxi2_map = this->_fe_map->get_d2psidxi2();
204  psi_map.resize (n_total_mapping_shape_functions);
205  dpsidxi_map.resize (n_total_mapping_shape_functions);
206  d2psidxi2_map.resize (n_total_mapping_shape_functions);
207 
208  // if (Dim == 3)
209  {
210  std::vector<std::vector<Real>> & dpsideta_map = this->_fe_map->get_dpsideta();
211  std::vector<std::vector<Real>> & d2psidxideta_map = this->_fe_map->get_d2psidxideta();
212  std::vector<std::vector<Real>> & d2psideta2_map = this->_fe_map->get_d2psideta2();
213  dpsideta_map.resize (n_total_mapping_shape_functions);
214  d2psidxideta_map.resize (n_total_mapping_shape_functions);
215  d2psideta2_map.resize (n_total_mapping_shape_functions);
216  }
217 
218  for (unsigned int i=0; i<n_total_mapping_shape_functions; i++)
219  {
220  psi_map[i].resize (n_total_qp);
221  dpsidxi_map[i].resize (n_total_qp);
222  d2psidxi2_map[i].resize (n_total_qp);
223 
224  // if (Dim == 3)
225  {
226  std::vector<std::vector<Real>> & dpsideta_map = this->_fe_map->get_dpsideta();
227  std::vector<std::vector<Real>> & d2psidxideta_map = this->_fe_map->get_d2psidxideta();
228  std::vector<std::vector<Real>> & d2psideta2_map = this->_fe_map->get_d2psideta2();
229  dpsideta_map[i].resize (n_total_qp);
230  d2psidxideta_map[i].resize (n_total_qp);
231  d2psideta2_map[i].resize (n_total_qp);
232  }
233  }
234  }
235 
236 
237  // compute shape maps
238  {
239  const std::vector<std::vector<Real>> & S_map = (base_fe->get_fe_map()).get_phi_map();
240  const std::vector<std::vector<Real>> & Ss_map = (base_fe->get_fe_map()).get_dphidxi_map();
241 
242  std::vector<std::vector<Real>> & psi_map = this->_fe_map->get_psi();
243  std::vector<std::vector<Real>> & dpsidxi_map = this->_fe_map->get_dpsidxi();
244  std::vector<std::vector<Real>> & dpsideta_map = this->_fe_map->get_dpsideta();
245 
246  for (unsigned int rp=0; rp<n_radial_qp; rp++) // over radial qps
247  for (unsigned int bp=0; bp<n_base_qp; bp++) // over base qps
248  for (unsigned int ti=0; ti<n_total_mapping_shape_functions; ti++) // over all mapping shapes
249  {
250  // let the index vectors take care of selecting the appropriate base/radial mapping shape
251  const unsigned int bi = _base_node_index [ti];
252  const unsigned int ri = _radial_node_index[ti];
253  psi_map [ti][bp+rp*n_base_qp] = S_map [bi][bp] * radial_map [ri][rp];
254  dpsidxi_map [ti][bp+rp*n_base_qp] = Ss_map[bi][bp] * radial_map [ri][rp];
255  dpsideta_map [ti][bp+rp*n_base_qp] = S_map [bi][bp] * dradialdv_map[ri][rp];
256 
257  // second derivatives are not implemented for infinite elements
258  // d2psidxi2_map [ti][bp+rp*n_base_qp] = 0.;
259  // d2psidxideta_map [ti][bp+rp*n_base_qp] = 0.;
260  // d2psideta2_map [ti][bp+rp*n_base_qp] = 0.;
261  }
262 
263  }
264 
265  }
266 
267  // quadrature rule weights
268  {
269  const std::vector<Real> & radial_qw = radial_qrule->get_weights();
270  const std::vector<Real> & base_qw = base_qrule->get_weights();
271 
272  libmesh_assert_equal_to (radial_qw.size(), n_radial_qp);
273  libmesh_assert_equal_to (base_qw.size(), n_base_qp);
274 
275  for (unsigned int rp=0; rp<n_radial_qp; rp++)
276  for (unsigned int bp=0; bp<n_base_qp; bp++)
277  {
278  _total_qrule_weights[bp + rp*n_base_qp] = radial_qw[rp] * base_qw[bp];
279  }
280  }
281 
282 }
283 
284 
285 
286 
287 // Explicit instantiations - doesn't make sense in 1D, but as
288 // long as we only return errors, we are fine... ;-)
289 //#include "libmesh/inf_fe_instantiate_1D.h"
290 //#include "libmesh/inf_fe_instantiate_2D.h"
291 //#include "libmesh/inf_fe_instantiate_3D.h"
292 INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, reinit(const Elem *, const unsigned int, const Real, const std::vector<Point> * const, const std::vector<Real> * const));
293 INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, reinit(const Elem *, const unsigned int, const Real, const std::vector<Point> * const, const std::vector<Real> * const));
294 INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, reinit(const Elem *, const unsigned int, const Real, const std::vector<Point> * const, const std::vector<Real> * const));
295 INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, edge_reinit(const Elem *, const unsigned int, const Real, const std::vector<Point> * const, const std::vector<Real> * const));
296 INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, edge_reinit(const Elem *, const unsigned int, const Real, const std::vector<Point> * const, const std::vector<Real> * const));
297 INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, edge_reinit(const Elem *, const unsigned int, const Real, const std::vector<Point> * const, const std::vector<Real> * const));
298 INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, init_face_shape_functions(const std::vector<Point> &, const Elem *));
299 INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, init_face_shape_functions(const std::vector<Point> &, const Elem *));
300 INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, init_face_shape_functions(const std::vector<Point> &, const Elem *));
301 
302 } // namespace libMesh
303 
304 #endif //ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
void init_face_shape_functions(const std::vector< Point > &qp, const Elem *side)
Not implemented yet.
virtual UniquePtr< Elem > build_side_ptr(const unsigned int i, bool proxy=true)=0
virtual ElemType type() const =0
unsigned int p_level() const
Definition: elem.h:2422
unsigned int dim
ElemType
Defines an enum for geometric element types.
unsigned short int side
Definition: xdr_io.C:49
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, Elem *, Base::build_elem(const Elem *))
const class libmesh_nullptr_t libmesh_nullptr
The libMesh namespace provides an interface to certain functionality in the library.
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
virtual void reinit(const Elem *elem, const std::vector< Point > *const pts=libmesh_nullptr, const std::vector< Real > *const weights=libmesh_nullptr) libmesh_override
This is at the core of this class.
Definition: inf_fe.C:107
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
virtual void edge_reinit(const Elem *elem, const unsigned int edge, const Real tolerance=TOLERANCE, const std::vector< Point > *const pts=libmesh_nullptr, const std::vector< Real > *const weights=libmesh_nullptr) libmesh_override
Not implemented yet.
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.