inf_fe.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 // Local includes
21 #include "libmesh/libmesh_config.h"
22 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
23 #include "libmesh/inf_fe.h"
25 #include "libmesh/elem.h"
27 
28 namespace libMesh
29 {
30 
31 
32 
33 // ------------------------------------------------------------
34 // InfFE class members
35 
36 
37 
38 // Constructor
39 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
41  FEBase (Dim, fet),
42 
43  _n_total_approx_sf (0),
44  _n_total_qp (0),
45 
46  base_qrule (NULL),
47  radial_qrule (NULL),
48  base_elem (NULL),
49  base_fe (NULL),
50 
51  // initialize the current_fe_type to all the same
52  // values as \p fet (since the FE families and coordinate
53  // map type should @e not change), but use an invalid order
54  // for the radial part (since this is the only order
55  // that may change!).
56  // the data structures like \p phi etc are not initialized
57  // through the constructor, but throught reinit()
58  current_fe_type ( FEType(fet.order,
59  fet.family,
61  fet.radial_family,
62  fet.inf_map) )
63 
64 {
65  // Sanity checks
66  libmesh_assert_equal_to (T_radial, fe_type.radial_family);
67  libmesh_assert_equal_to (T_map, fe_type.inf_map);
68 
69  // build the base_fe object, handle the AutoPtr
70  if (Dim != 1)
71  {
72  AutoPtr<FEBase> ap_fb(FEBase::build(Dim-1, fet));
73  base_fe = ap_fb.release();
74  }
75 }
76 
77 
78 
79 
80 // Destructor
81 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
83 {
84  // delete pointers, if necessary
85  delete base_qrule;
86  base_qrule = NULL;
87 
88  delete radial_qrule;
89  radial_qrule = NULL;
90 
91  delete base_elem;
92  base_elem = NULL;
93 
94  delete base_fe;
95  base_fe = NULL;
96 }
97 
98 
99 
100 
101 
102 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
104 {
105  libmesh_assert(q);
106  libmesh_assert(base_fe);
107 
108  const Order base_int_order = q->get_order();
109  const Order radial_int_order = static_cast<Order>(2 * (static_cast<unsigned int>(fe_type.radial_order) + 1) +2);
110  const unsigned int qrule_dim = q->get_dim();
111 
112  if (Dim != 1)
113  {
114  // build a Dim-1 quadrature rule of the type that we received
115  AutoPtr<QBase> apq( QBase::build(q->type(), qrule_dim-1, base_int_order) );
116  base_qrule = apq.release();
117  base_fe->attach_quadrature_rule(base_qrule);
118  }
119 
120  // in radial direction, always use Gauss quadrature
121  radial_qrule = new QGauss(1, radial_int_order);
122 
123  // currently not used. But maybe helpful to store the QBase*
124  // with which we initialized our own quadrature rules
125  qrule = q;
126 }
127 
128 
129 
130 
131 
132 template <unsigned int Dim, FEFamily T_radial, InfMapType T_base>
134 {
135  if (base_elem != NULL)
136  delete base_elem;
137  base_elem = Base::build_elem(inf_elem);
138 }
139 
140 
141 
142 
143 
144 
145 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
147  const std::vector<Point>* const pts,
148  const std::vector<Real>* const weights)
149 {
150  libmesh_assert(base_fe);
151  libmesh_assert(base_fe->qrule);
152  libmesh_assert_equal_to (base_fe->qrule, base_qrule);
153  libmesh_assert(radial_qrule);
154  libmesh_assert(inf_elem);
155 
156  if (pts == NULL)
157  {
158  bool init_shape_functions_required = false;
159 
160  // -----------------------------------------------------------------
161  // init the radial data fields only when the radial order changes
162  if (current_fe_type.radial_order != fe_type.radial_order)
163  {
164  current_fe_type.radial_order = fe_type.radial_order;
165 
166  // Watch out: this call to QBase->init() only works for
167  // current_fe_type = const! To allow variable Order,
168  // the init() of QBase has to be modified...
169  radial_qrule->init(EDGE2);
170 
171  // initialize the radial shape functions
172  this->init_radial_shape_functions(inf_elem);
173 
174  init_shape_functions_required=true;
175  }
176 
177 
178  bool update_base_elem_required=true;
179 
180  // -----------------------------------------------------------------
181  // update the type in accordance to the current cell
182  // and reinit if the cell type has changed or (as in
183  // the case of the hierarchics) the shape functions
184  // depend on the particular element and need a reinit
185  if ( ( Dim != 1) &&
186  ( (this->get_type() != inf_elem->type()) ||
187  (base_fe->shapes_need_reinit()) ) )
188  {
189  // store the new element type, update base_elem
190  // here. Through \p update_base_elem_required,
191  // remember whether it has to be updated (see below).
192  elem_type = inf_elem->type();
193  this->update_base_elem(inf_elem);
194  update_base_elem_required=false;
195 
196  // initialize the base quadrature rule for the new element
197  base_qrule->init(base_elem->type());
198 
199  // initialize the shape functions in the base
200  base_fe->init_base_shape_functions(base_fe->qrule->get_points(),
201  base_elem);
202 
203  init_shape_functions_required=true;
204  }
205 
206 
207  // when either the radial or base part change,
208  // we have to init the whole fields
209  if (init_shape_functions_required)
210  this->init_shape_functions (inf_elem);
211 
212  // computing the distance only works when we have the current
213  // base_elem stored. This happens when fe_type is const,
214  // the inf_elem->type remains the same. Then we have to
215  // update the base elem _here_.
216  if (update_base_elem_required)
217  this->update_base_elem(inf_elem);
218 
219  // compute dist (depends on geometry, therefore has to be updated for
220  // each and every new element), throw radial and base part together
221  this->combine_base_radial (inf_elem);
222 
223  this->_fe_map->compute_map (this->dim,_total_qrule_weights, inf_elem);
224 
225  // Compute the shape functions and the derivatives
226  // at all quadrature points.
227  this->compute_shape_functions (inf_elem,base_fe->qrule->get_points());
228  }
229 
230  else // if pts != NULL
231  {
232  // update the elem_type
233  elem_type = inf_elem->type();
234 
235  // init radial shapes
236  this->init_radial_shape_functions(inf_elem);
237 
238  // update the base
239  this->update_base_elem(inf_elem);
240 
241  // the finite element on the ifem base
242  {
243  AutoPtr<FEBase> ap_fb(FEBase::build(Dim-1, this->fe_type));
244  if (base_fe != NULL)
245  delete base_fe;
246  base_fe = ap_fb.release();
247  }
248 
249  // inite base shapes
250  base_fe->init_base_shape_functions(*pts,
251  base_elem);
252 
253  this->init_shape_functions (inf_elem);
254 
255  // combine the base and radial shapes
256  this->combine_base_radial (inf_elem);
257 
258  // weights
259  if (weights != NULL)
260  {
261  this->_fe_map->compute_map (this->dim, *weights, inf_elem);
262  }
263  else
264  {
265  std::vector<Real> dummy_weights (pts->size(), 1.);
266  this->_fe_map->compute_map (this->dim, dummy_weights, inf_elem);
267  }
268 
269  // finally compute the ifem shapes
270  this->compute_shape_functions (inf_elem,*pts);
271  }
272 
273 }
274 
275 
276 
277 
278 
279 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
280 void InfFE<Dim,T_radial,T_map>::init_radial_shape_functions(const Elem* libmesh_dbg_var(inf_elem))
281 {
282  libmesh_assert(radial_qrule);
283  libmesh_assert(inf_elem);
284 
285 
289  START_LOG("init_radial_shape_functions()", "InfFE");
290 
291 
292  // -----------------------------------------------------------------
293  // initialize most of the things related to mapping
294 
295  // The order to use in the radial map (currently independent of the element type)
296  const Order radial_mapping_order (Radial::mapping_order());
297  const unsigned int n_radial_mapping_shape_functions (Radial::n_dofs(radial_mapping_order));
298 
299 
300 
301  // -----------------------------------------------------------------
302  // initialize most of the things related to physical approximation
303 
304  const Order radial_approx_order (fe_type.radial_order);
305  const unsigned int n_radial_approx_shape_functions (Radial::n_dofs(radial_approx_order));
306 
307  const unsigned int n_radial_qp = radial_qrule->n_points();
308  const std::vector<Point>& radial_qp = radial_qrule->get_points();
309 
310 
311 
312  // -----------------------------------------------------------------
313  // resize the radial data fields
314 
315  mode.resize (n_radial_approx_shape_functions); // the radial polynomials (eval)
316  dmodedv.resize (n_radial_approx_shape_functions);
317 
318  som.resize (n_radial_qp); // the (1-v)/2 weight
319  dsomdv.resize (n_radial_qp);
320 
321  radial_map.resize (n_radial_mapping_shape_functions); // the radial map
322  dradialdv_map.resize (n_radial_mapping_shape_functions);
323 
324 
325  for (unsigned int i=0; i<n_radial_mapping_shape_functions; i++)
326  {
327  radial_map[i].resize (n_radial_qp);
328  dradialdv_map[i].resize (n_radial_qp);
329  }
330 
331 
332  for (unsigned int i=0; i<n_radial_approx_shape_functions; i++)
333  {
334  mode[i].resize (n_radial_qp);
335  dmodedv[i].resize (n_radial_qp);
336  }
337 
338 
339  // compute scalar values at radial quadrature points
340  for (unsigned int p=0; p<n_radial_qp; p++)
341  {
342  som[p] = Radial::decay (radial_qp[p](0));
343  dsomdv[p] = Radial::decay_deriv (radial_qp[p](0));
344  }
345 
346 
347  // evaluate the mode shapes in radial direction at radial quadrature points
348  for (unsigned int i=0; i<n_radial_approx_shape_functions; i++)
349  for (unsigned int p=0; p<n_radial_qp; p++)
350  {
351  mode[i][p] = InfFE<Dim,T_radial,T_map>::eval (radial_qp[p](0), radial_approx_order, i);
352  dmodedv[i][p] = InfFE<Dim,T_radial,T_map>::eval_deriv (radial_qp[p](0), radial_approx_order, i);
353  }
354 
355 
356  // evaluate the mapping functions in radial direction at radial quadrature points
357  for (unsigned int i=0; i<n_radial_mapping_shape_functions; i++)
358  for (unsigned int p=0; p<n_radial_qp; p++)
359  {
360  radial_map[i][p] = InfFE<Dim,INFINITE_MAP,T_map>::eval (radial_qp[p](0), radial_mapping_order, i);
361  dradialdv_map[i][p] = InfFE<Dim,INFINITE_MAP,T_map>::eval_deriv (radial_qp[p](0), radial_mapping_order, i);
362  }
363 
367  STOP_LOG("init_radial_shape_functions()", "InfFE");
368 
369 }
370 
371 
372 
373 
374 
375 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
377 {
378  libmesh_assert(inf_elem);
379 
380 
381  // Start logging the radial shape function initialization
382  START_LOG("init_shape_functions()", "InfFE");
383 
384 
385  // -----------------------------------------------------------------
386  // fast access to some const int's for the radial data
387  const unsigned int n_radial_mapping_sf =
388  libmesh_cast_int<unsigned int>(radial_map.size());
389  const unsigned int n_radial_approx_sf =
390  libmesh_cast_int<unsigned int>(mode.size());
391  const unsigned int n_radial_qp =
392  libmesh_cast_int<unsigned int>(som.size());
393 
394 
395  // -----------------------------------------------------------------
396  // initialize most of the things related to mapping
397 
398  // The element type and order to use in the base map
399  const Order base_mapping_order ( base_elem->default_order() );
400  const ElemType base_mapping_elem_type ( base_elem->type() );
401 
402  // the number of base shape functions used to construct the map
403  // (Lagrange shape functions are used for mapping in the base)
404  unsigned int n_base_mapping_shape_functions = Base::n_base_mapping_sf(base_mapping_elem_type,
405  base_mapping_order);
406 
407  const unsigned int n_total_mapping_shape_functions =
408  n_radial_mapping_sf * n_base_mapping_shape_functions;
409 
410 
411 
412  // -----------------------------------------------------------------
413  // initialize most of the things related to physical approximation
414 
415  unsigned int n_base_approx_shape_functions;
416  if (Dim > 1)
417  n_base_approx_shape_functions = base_fe->n_shape_functions();
418  else
419  n_base_approx_shape_functions = 1;
420 
421 
422  const unsigned int n_total_approx_shape_functions =
423  n_radial_approx_sf * n_base_approx_shape_functions;
424 
425  // update class member field
426  _n_total_approx_sf = n_total_approx_shape_functions;
427 
428 
429  // The number of the base quadrature points.
430  const unsigned int n_base_qp = base_qrule->n_points();
431 
432  // The total number of quadrature points.
433  const unsigned int n_total_qp = n_radial_qp * n_base_qp;
434 
435 
436  // update class member field
437  _n_total_qp = n_total_qp;
438 
439 
440 
441  // -----------------------------------------------------------------
442  // initialize the node and shape numbering maps
443  {
444  // these vectors work as follows: the i-th entry stores
445  // the associated base/radial node number
446  _radial_node_index.resize (n_total_mapping_shape_functions);
447  _base_node_index.resize (n_total_mapping_shape_functions);
448 
449  // similar for the shapes: the i-th entry stores
450  // the associated base/radial shape number
451  _radial_shape_index.resize (n_total_approx_shape_functions);
452  _base_shape_index.resize (n_total_approx_shape_functions);
453 
454  const ElemType inf_elem_type (inf_elem->type());
455 
456  // fill the node index map
457  for (unsigned int n=0; n<n_total_mapping_shape_functions; n++)
458  {
459  compute_node_indices (inf_elem_type,
460  n,
461  _base_node_index[n],
462  _radial_node_index[n]);
463  libmesh_assert_less (_base_node_index[n], n_base_mapping_shape_functions);
464  libmesh_assert_less (_radial_node_index[n], n_radial_mapping_sf);
465  }
466 
467  // fill the shape index map
468  for (unsigned int n=0; n<n_total_approx_shape_functions; n++)
469  {
470  compute_shape_indices (this->fe_type,
471  inf_elem_type,
472  n,
473  _base_shape_index[n],
474  _radial_shape_index[n]);
475  libmesh_assert_less (_base_shape_index[n], n_base_approx_shape_functions);
476  libmesh_assert_less (_radial_shape_index[n], n_radial_approx_sf);
477  }
478  }
479 
480 
481 
482 
483 
484  // -----------------------------------------------------------------
485  // resize the base data fields
486  dist.resize(n_base_mapping_shape_functions);
487 
488 
489 
490  // -----------------------------------------------------------------
491  // resize the total data fields
492 
493  // the phase term varies with xi, eta and zeta(v): store it for _all_ qp
494  //
495  // when computing the phase, we need the base approximations
496  // therefore, initialize the phase here, but evaluate it
497  // in combine_base_radial().
498  //
499  // the weight, though, is only needed at the radial quadrature points, n_radial_qp.
500  // but for a uniform interface to the protected data fields
501  // the weight data field (which are accessible from the outside) are expanded to n_total_qp.
502  weight.resize (n_total_qp);
503  dweightdv.resize (n_total_qp);
504  dweight.resize (n_total_qp);
505 
506  dphase.resize (n_total_qp);
507  dphasedxi.resize (n_total_qp);
508  dphasedeta.resize (n_total_qp);
509  dphasedzeta.resize (n_total_qp);
510 
511  // this vector contains the integration weights for the combined quadrature rule
512  _total_qrule_weights.resize(n_total_qp);
513 
514 
515  // -----------------------------------------------------------------
516  // InfFE's data fields phi, dphi, dphidx, phi_map etc hold the _total_
517  // shape and mapping functions, respectively
518  {
519  phi.resize (n_total_approx_shape_functions);
520  dphi.resize (n_total_approx_shape_functions);
521  dphidx.resize (n_total_approx_shape_functions);
522  dphidy.resize (n_total_approx_shape_functions);
523  dphidz.resize (n_total_approx_shape_functions);
524  dphidxi.resize (n_total_approx_shape_functions);
525 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
526  libmesh_do_once(libMesh::err << "Second derivatives for Infinite elements"
527  << " are not yet implemented!"
528  << std::endl);
529 
530  d2phi.resize (n_total_approx_shape_functions);
531  d2phidx2.resize (n_total_approx_shape_functions);
532  d2phidxdy.resize (n_total_approx_shape_functions);
533  d2phidxdz.resize (n_total_approx_shape_functions);
534  d2phidy2.resize (n_total_approx_shape_functions);
535  d2phidydz.resize (n_total_approx_shape_functions);
536  d2phidz2.resize (n_total_approx_shape_functions);
537  d2phidxi2.resize (n_total_approx_shape_functions);
538 
539  if (Dim > 1)
540  {
541  d2phidxideta.resize (n_total_approx_shape_functions);
542  d2phideta2.resize (n_total_approx_shape_functions);
543  }
544 
545  if (Dim > 2)
546  {
547  d2phidetadzeta.resize (n_total_approx_shape_functions);
548  d2phidxidzeta.resize (n_total_approx_shape_functions);
549  d2phidzeta2.resize (n_total_approx_shape_functions);
550  }
551 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
552 
553  if (Dim > 1)
554  dphideta.resize (n_total_approx_shape_functions);
555 
556  if (Dim == 3)
557  dphidzeta.resize (n_total_approx_shape_functions);
558 
559 
560 
561  std::vector<std::vector<Real> >& phi_map = this->_fe_map->get_phi_map();
562  std::vector<std::vector<Real> >& dphidxi_map = this->_fe_map->get_dphidxi_map();
563 
564  phi_map.resize (n_total_mapping_shape_functions);
565  dphidxi_map.resize (n_total_mapping_shape_functions);
566 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
567  std::vector<std::vector<Real> >& d2phidxi2_map = this->_fe_map->get_d2phidxi2_map();
568  d2phidxi2_map.resize (n_total_mapping_shape_functions);
569 
570  if (Dim > 1)
571  {
572  std::vector<std::vector<Real> >& d2phidxideta_map = this->_fe_map->get_d2phidxideta_map();
573  std::vector<std::vector<Real> >& d2phideta2_map = this->_fe_map->get_d2phideta2_map();
574  d2phidxideta_map.resize (n_total_mapping_shape_functions);
575  d2phideta2_map.resize (n_total_mapping_shape_functions);
576  }
577 
578  if (Dim == 3)
579  {
580  std::vector<std::vector<Real> >& d2phidxidzeta_map = this->_fe_map->get_d2phidxidzeta_map();
581  std::vector<std::vector<Real> >& d2phidetadzeta_map = this->_fe_map->get_d2phidetadzeta_map();
582  std::vector<std::vector<Real> >& d2phidzeta2_map = this->_fe_map->get_d2phidzeta2_map();
583  d2phidxidzeta_map.resize (n_total_mapping_shape_functions);
584  d2phidetadzeta_map.resize (n_total_mapping_shape_functions);
585  d2phidzeta2_map.resize (n_total_mapping_shape_functions);
586  }
587 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
588 
589  if (Dim > 1)
590  {
591  std::vector<std::vector<Real> >& dphideta_map = this->_fe_map->get_dphideta_map();
592  dphideta_map.resize (n_total_mapping_shape_functions);
593  }
594 
595  if (Dim == 3)
596  {
597  std::vector<std::vector<Real> >& dphidzeta_map = this->_fe_map->get_dphidzeta_map();
598  dphidzeta_map.resize (n_total_mapping_shape_functions);
599  }
600  }
601 
602 
603 
604  // -----------------------------------------------------------------
605  // collect all the for loops, where inner vectors are
606  // resized to the appropriate number of quadrature points
607  {
608  for (unsigned int i=0; i<n_total_approx_shape_functions; i++)
609  {
610  phi[i].resize (n_total_qp);
611  dphi[i].resize (n_total_qp);
612  dphidx[i].resize (n_total_qp);
613  dphidy[i].resize (n_total_qp);
614  dphidz[i].resize (n_total_qp);
615  dphidxi[i].resize (n_total_qp);
616 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
617  d2phi[i].resize (n_total_qp);
618  d2phidx2[i].resize (n_total_qp);
619  d2phidxdy[i].resize (n_total_qp);
620  d2phidxdz[i].resize (n_total_qp);
621  d2phidy2[i].resize (n_total_qp);
622  d2phidydz[i].resize (n_total_qp);
623  d2phidy2[i].resize (n_total_qp);
624  d2phidxi2[i].resize (n_total_qp);
625 
626  if (Dim > 1)
627  {
628  d2phidxideta[i].resize (n_total_qp);
629  d2phideta2[i].resize (n_total_qp);
630  }
631  if (Dim > 2)
632  {
633  d2phidxidzeta[i].resize (n_total_qp);
634  d2phidetadzeta[i].resize (n_total_qp);
635  d2phidzeta2[i].resize (n_total_qp);
636  }
637 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
638 
639  if (Dim > 1)
640  dphideta[i].resize (n_total_qp);
641 
642  if (Dim == 3)
643  dphidzeta[i].resize (n_total_qp);
644 
645  }
646 
647  for (unsigned int i=0; i<n_total_mapping_shape_functions; i++)
648  {
649  std::vector<std::vector<Real> >& phi_map = this->_fe_map->get_phi_map();
650  std::vector<std::vector<Real> >& dphidxi_map = this->_fe_map->get_dphidxi_map();
651  phi_map[i].resize (n_total_qp);
652  dphidxi_map[i].resize (n_total_qp);
653 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
654  std::vector<std::vector<Real> >& d2phidxi2_map = this->_fe_map->get_d2phidxi2_map();
655  d2phidxi2_map[i].resize (n_total_qp);
656  if (Dim > 1)
657  {
658  std::vector<std::vector<Real> >& d2phidxideta_map = this->_fe_map->get_d2phidxideta_map();
659  std::vector<std::vector<Real> >& d2phideta2_map = this->_fe_map->get_d2phideta2_map();
660  d2phidxideta_map[i].resize (n_total_qp);
661  d2phideta2_map[i].resize (n_total_qp);
662  }
663 
664  if (Dim > 2)
665  {
666  std::vector<std::vector<Real> >& d2phidxidzeta_map = this->_fe_map->get_d2phidxidzeta_map();
667  std::vector<std::vector<Real> >& d2phidetadzeta_map = this->_fe_map->get_d2phidetadzeta_map();
668  std::vector<std::vector<Real> >& d2phidzeta2_map = this->_fe_map->get_d2phidzeta2_map();
669  d2phidxidzeta_map[i].resize (n_total_qp);
670  d2phidetadzeta_map[i].resize (n_total_qp);
671  d2phidzeta2_map[i].resize (n_total_qp);
672  }
673 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
674 
675  if (Dim > 1)
676  {
677  std::vector<std::vector<Real> >& dphideta_map = this->_fe_map->get_dphideta_map();
678  dphideta_map[i].resize (n_total_qp);
679  }
680 
681  if (Dim == 3)
682  {
683  std::vector<std::vector<Real> >& dphidzeta_map = this->_fe_map->get_dphidzeta_map();
684  dphidzeta_map[i].resize (n_total_qp);
685  }
686  }
687  }
688 
689 
690 
691  {
692  // -----------------------------------------------------------------
693  // (a) compute scalar values at _all_ quadrature points -- for uniform
694  // access from the outside to these fields
695  // (b) form a std::vector<Real> which contains the appropriate weights
696  // of the combined quadrature rule!
697  const std::vector<Point>& radial_qp = radial_qrule->get_points();
698  libmesh_assert_equal_to (radial_qp.size(), n_radial_qp);
699 
700  const std::vector<Real>& radial_qw = radial_qrule->get_weights();
701  const std::vector<Real>& base_qw = base_qrule->get_weights();
702  libmesh_assert_equal_to (radial_qw.size(), n_radial_qp);
703  libmesh_assert_equal_to (base_qw.size(), n_base_qp);
704 
705  for (unsigned int rp=0; rp<n_radial_qp; rp++)
706  for (unsigned int bp=0; bp<n_base_qp; bp++)
707  {
708  weight [ bp+rp*n_base_qp ] = Radial::D (radial_qp[rp](0));
709  dweightdv[ bp+rp*n_base_qp ] = Radial::D_deriv (radial_qp[rp](0));
710 
711  _total_qrule_weights[ bp+rp*n_base_qp ] = radial_qw[rp] * base_qw[bp];
712  }
713  }
714 
715 
719  STOP_LOG("init_shape_functions()", "InfFE");
720 
721 }
722 
723 
724 
725 
726 
727 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
729 {
730  libmesh_assert(inf_elem);
731  // at least check whether the base element type is correct.
732  // otherwise this version of computing dist would give problems
733  libmesh_assert_equal_to (base_elem->type(), Base::get_elem_type(inf_elem->type()));
734 
735 
739  START_LOG("combine_base_radial()", "InfFE");
740 
741 
742  // zero the phase, since it is to be summed up
743  std::fill (dphasedxi.begin(), dphasedxi.end(), 0.);
744  std::fill (dphasedeta.begin(), dphasedeta.end(), 0.);
745  std::fill (dphasedzeta.begin(), dphasedzeta.end(), 0.);
746 
747 
748  const unsigned int n_base_mapping_sf =
749  libmesh_cast_int<unsigned int>(dist.size());
750  const Point origin = inf_elem->origin();
751 
752  // for each new infinite element, compute the radial distances
753  for (unsigned int n=0; n<n_base_mapping_sf; n++)
754  dist[n] = Point(base_elem->point(n) - origin).size();
755 
756 
757  switch (Dim)
758  {
759 
760  //------------------------------------------------------------
761  // 1D
762  case 1:
763  {
764  libmesh_not_implemented();
765  break;
766  }
767 
768 
769 
770  //------------------------------------------------------------
771  // 2D
772  case 2:
773  {
774  libmesh_not_implemented();
775  break;
776  }
777 
778 
779 
780  //------------------------------------------------------------
781  // 3D
782  case 3:
783  {
784  // fast access to the approximation and mapping shapes of base_fe
785  const std::vector<std::vector<Real> >& S = base_fe->phi;
786  const std::vector<std::vector<Real> >& Ss = base_fe->dphidxi;
787  const std::vector<std::vector<Real> >& St = base_fe->dphideta;
788  const std::vector<std::vector<Real> >& S_map = (base_fe->get_fe_map()).get_phi_map();
789  const std::vector<std::vector<Real> >& Ss_map = (base_fe->get_fe_map()).get_dphidxi_map();
790  const std::vector<std::vector<Real> >& St_map = (base_fe->get_fe_map()).get_dphideta_map();
791 
792  const unsigned int n_radial_qp = radial_qrule->n_points();
793  const unsigned int n_base_qp = base_qrule-> n_points();
794 
795  const unsigned int n_total_mapping_sf =
796  libmesh_cast_int<unsigned int>(radial_map.size()) * n_base_mapping_sf;
797 
798  const unsigned int n_total_approx_sf = Radial::n_dofs(fe_type.radial_order) * base_fe->n_shape_functions();
799 
800 
801  // compute the phase term derivatives
802  {
803  unsigned int tp=0;
804  for (unsigned int rp=0; rp<n_radial_qp; rp++) // over radial qp's
805  for (unsigned int bp=0; bp<n_base_qp; bp++) // over base qp's
806  {
807  // sum over all base shapes, to get the average distance
808  for (unsigned int i=0; i<n_base_mapping_sf; i++)
809  {
810  dphasedxi[tp] += Ss_map[i][bp] * dist[i] * radial_map [1][rp];
811  dphasedeta[tp] += St_map[i][bp] * dist[i] * radial_map [1][rp];
812  dphasedzeta[tp] += S_map [i][bp] * dist[i] * dradialdv_map[1][rp];
813  }
814 
815  tp++;
816 
817  } // loop radial and base qp's
818 
819  }
820 
821  libmesh_assert_equal_to (phi.size(), n_total_approx_sf);
822  libmesh_assert_equal_to (dphidxi.size(), n_total_approx_sf);
823  libmesh_assert_equal_to (dphideta.size(), n_total_approx_sf);
824  libmesh_assert_equal_to (dphidzeta.size(), n_total_approx_sf);
825 
826  // compute the overall approximation shape functions,
827  // pick the appropriate radial and base shapes through using
828  // _base_shape_index and _radial_shape_index
829  for (unsigned int rp=0; rp<n_radial_qp; rp++) // over radial qp's
830  for (unsigned int bp=0; bp<n_base_qp; bp++) // over base qp's
831  for (unsigned int ti=0; ti<n_total_approx_sf; ti++) // over _all_ approx_sf
832  {
833  // let the index vectors take care of selecting the appropriate base/radial shape
834  const unsigned int bi = _base_shape_index [ti];
835  const unsigned int ri = _radial_shape_index[ti];
836  phi [ti][bp+rp*n_base_qp] = S [bi][bp] * mode[ri][rp] * som[rp];
837  dphidxi [ti][bp+rp*n_base_qp] = Ss[bi][bp] * mode[ri][rp] * som[rp];
838  dphideta [ti][bp+rp*n_base_qp] = St[bi][bp] * mode[ri][rp] * som[rp];
839  dphidzeta[ti][bp+rp*n_base_qp] = S [bi][bp]
840  * (dmodedv[ri][rp] * som[rp] + mode[ri][rp] * dsomdv[rp]);
841  }
842 
843  std::vector<std::vector<Real> >& phi_map = this->_fe_map->get_phi_map();
844  std::vector<std::vector<Real> >& dphidxi_map = this->_fe_map->get_dphidxi_map();
845  std::vector<std::vector<Real> >& dphideta_map = this->_fe_map->get_dphideta_map();
846  std::vector<std::vector<Real> >& dphidzeta_map = this->_fe_map->get_dphidzeta_map();
847 
848  libmesh_assert_equal_to (phi_map.size(), n_total_mapping_sf);
849  libmesh_assert_equal_to (dphidxi_map.size(), n_total_mapping_sf);
850  libmesh_assert_equal_to (dphideta_map.size(), n_total_mapping_sf);
851  libmesh_assert_equal_to (dphidzeta_map.size(), n_total_mapping_sf);
852 
853  // compute the overall mapping functions,
854  // pick the appropriate radial and base entries through using
855  // _base_node_index and _radial_node_index
856  for (unsigned int rp=0; rp<n_radial_qp; rp++) // over radial qp's
857  for (unsigned int bp=0; bp<n_base_qp; bp++) // over base qp's
858  for (unsigned int ti=0; ti<n_total_mapping_sf; ti++) // over all mapping shapes
859  {
860  // let the index vectors take care of selecting the appropriate base/radial mapping shape
861  const unsigned int bi = _base_node_index [ti];
862  const unsigned int ri = _radial_node_index[ti];
863  phi_map [ti][bp+rp*n_base_qp] = S_map [bi][bp] * radial_map [ri][rp];
864  dphidxi_map [ti][bp+rp*n_base_qp] = Ss_map[bi][bp] * radial_map [ri][rp];
865  dphideta_map [ti][bp+rp*n_base_qp] = St_map[bi][bp] * radial_map [ri][rp];
866  dphidzeta_map[ti][bp+rp*n_base_qp] = S_map [bi][bp] * dradialdv_map[ri][rp];
867  }
868 
869 
870  break;
871  }
872 
873 
874  default:
875  libmesh_error();
876  }
877 
878 
882  STOP_LOG("combine_base_radial()", "InfFE");
883 
884 }
885 
886 
887 
888 
889 
890 
891 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
892 void InfFE<Dim,T_radial,T_map>::compute_shape_functions(const Elem*, const std::vector<Point>&)
893 {
894  libmesh_assert(radial_qrule);
895 
896 
897 
898  // Start logging the overall computation of shape functions
899  START_LOG("compute_shape_functions()", "InfFE");
900 
901 
902  const unsigned int n_total_qp = _n_total_qp;
903 
904 
905  //-------------------------------------------------------------------------
906  // Compute the shape function values (and derivatives)
907  // at the Quadrature points. Note that the actual values
908  // have already been computed via init_shape_functions
909 
910  // Compute the value of the derivative shape function i at quadrature point p
911  switch (dim)
912  {
913 
914  case 1:
915  {
916  libmesh_not_implemented();
917  break;
918  }
919 
920  case 2:
921  {
922  libmesh_not_implemented();
923  break;
924  }
925 
926  case 3:
927  {
928  const std::vector<Real>& dxidx_map = this->_fe_map->get_dxidx();
929  const std::vector<Real>& dxidy_map = this->_fe_map->get_dxidy();
930  const std::vector<Real>& dxidz_map = this->_fe_map->get_dxidz();
931 
932  const std::vector<Real>& detadx_map = this->_fe_map->get_detadx();
933  const std::vector<Real>& detady_map = this->_fe_map->get_detady();
934  const std::vector<Real>& detadz_map = this->_fe_map->get_detadz();
935 
936  const std::vector<Real>& dzetadx_map = this->_fe_map->get_dzetadx();
937  const std::vector<Real>& dzetady_map = this->_fe_map->get_dzetady();
938  const std::vector<Real>& dzetadz_map = this->_fe_map->get_dzetadz();
939 
940  // These are _all_ shape functions of this infinite element
941  for (unsigned int i=0; i<phi.size(); i++)
942  for (unsigned int p=0; p<n_total_qp; p++)
943  {
944  // dphi/dx = (dphi/dxi)*(dxi/dx) + (dphi/deta)*(deta/dx) + (dphi/dzeta)*(dzeta/dx);
945  dphi[i][p](0) =
946  dphidx[i][p] = (dphidxi[i][p]*dxidx_map[p] +
947  dphideta[i][p]*detadx_map[p] +
948  dphidzeta[i][p]*dzetadx_map[p]);
949 
950  // dphi/dy = (dphi/dxi)*(dxi/dy) + (dphi/deta)*(deta/dy) + (dphi/dzeta)*(dzeta/dy);
951  dphi[i][p](1) =
952  dphidy[i][p] = (dphidxi[i][p]*dxidy_map[p] +
953  dphideta[i][p]*detady_map[p] +
954  dphidzeta[i][p]*dzetady_map[p]);
955 
956  // dphi/dz = (dphi/dxi)*(dxi/dz) + (dphi/deta)*(deta/dz) + (dphi/dzeta)*(dzeta/dz);
957  dphi[i][p](2) =
958  dphidz[i][p] = (dphidxi[i][p]*dxidz_map[p] +
959  dphideta[i][p]*detadz_map[p] +
960  dphidzeta[i][p]*dzetadz_map[p]);
961  }
962 
963 
964  // This is the derivative of the phase term of this infinite element
965  for (unsigned int p=0; p<n_total_qp; p++)
966  {
967  // the derivative of the phase term
968  dphase[p](0) = (dphasedxi[p] * dxidx_map[p] +
969  dphasedeta[p] * detadx_map[p] +
970  dphasedzeta[p] * dzetadx_map[p]);
971 
972  dphase[p](1) = (dphasedxi[p] * dxidy_map[p] +
973  dphasedeta[p] * detady_map[p] +
974  dphasedzeta[p] * dzetady_map[p]);
975 
976  dphase[p](2) = (dphasedxi[p] * dxidz_map[p] +
977  dphasedeta[p] * detadz_map[p] +
978  dphasedzeta[p] * dzetadz_map[p]);
979 
980  // the derivative of the radial weight - varies only in radial direction,
981  // therefore dweightdxi = dweightdeta = 0.
982  dweight[p](0) = dweightdv[p] * dzetadx_map[p];
983 
984  dweight[p](1) = dweightdv[p] * dzetady_map[p];
985 
986  dweight[p](2) = dweightdv[p] * dzetadz_map[p];
987 
988  }
989 
990  break;
991  }
992 
993 
994 
995  default:
996  {
997  libmesh_error();
998  }
999  }
1000 
1001 
1002 
1003  // Stop logging the overall computation of shape functions
1004  STOP_LOG("compute_shape_functions()", "InfFE");
1005 
1006 }
1007 
1008 
1009 
1010 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
1012 {
1013  return false;
1014 }
1015 
1016 } // namespace libMesh
1017 
1018 
1019 //--------------------------------------------------------------
1020 // Explicit instantiations
1024 
1025 
1026 
1027 #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