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

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

Hosted By:
SourceForge.net Logo