Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
fe_field_function.templates.h
1 // ---------------------------------------------------------------------
2 // @f$Id: fe_field_function.templates.h 30040 2013-07-18 17:06:48Z maier @f$
3 //
4 // Copyright (C) 2007 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 
18 #include <deal.II/base/utilities.h>
19 #include <deal.II/base/logstream.h>
20 #include <deal.II/grid/grid_tools.h>
21 #include <deal.II/hp/fe_collection.h>
22 #include <deal.II/hp/fe_values.h>
23 #include <deal.II/hp/mapping_collection.h>
24 #include <deal.II/hp/q_collection.h>
25 #include <deal.II/fe/fe_values.h>
26 #include <deal.II/numerics/fe_field_function.h>
27 
28 
30 
31 namespace Functions
32 {
33 
34  template <int dim, typename DH, typename VECTOR>
36  const VECTOR &myv,
37  const Mapping<dim> &mymapping)
38  :
39  Function<dim>(mydh.get_fe().n_components()),
40  dh(&mydh, "FEFieldFunction"),
41  data_vector(myv),
42  mapping(mymapping),
43  cell_hint(dh->end()),
44  n_components(mydh.get_fe().n_components())
45  {
46  }
47 
48 
49 
50  template <int dim, typename DH, typename VECTOR>
51  void
53  set_active_cell(const typename DH::active_cell_iterator &newcell)
54  {
55  cell_hint.get() = newcell;
56  }
57 
58 
59 
60  template <int dim, typename DH, typename VECTOR>
62  Vector<double> &values) const
63  {
64  Assert (values.size() == n_components,
65  ExcDimensionMismatch(values.size(), n_components));
66  typename DH::active_cell_iterator cell = cell_hint.get();
67  if (cell == dh->end())
68  cell = dh->begin_active();
69 
70  boost::optional<Point<dim> >
71  qp = get_reference_coordinates (cell, p);
72  if (!qp)
73  {
74  const std::pair<typename DH::active_cell_iterator, Point<dim> > my_pair
75  = GridTools::find_active_cell_around_point (mapping, *dh, p);
76  AssertThrow (my_pair.first->is_locally_owned(),
77  ExcPointNotAvailableHere());
78 
79  cell = my_pair.first;
80  qp = my_pair.second;
81  }
82 
83  cell_hint.get() = cell;
84 
85  // Now we can find out about the point
86  Quadrature<dim> quad(qp.get());
87  FEValues<dim> fe_v(mapping, cell->get_fe(), quad,
89  fe_v.reinit(cell);
90  std::vector< Vector<double> > vvalues (1, values);
91  fe_v.get_function_values(data_vector, vvalues);
92  values = vvalues[0];
93  }
94 
95 
96 
97  template <int dim, typename DH, typename VECTOR>
98  double
100  const unsigned int comp) const
101  {
102  Vector<double> values(n_components);
103  vector_value(p, values);
104  return values(comp);
105  }
106 
107 
108 
109  template <int dim, typename DH, typename VECTOR>
110  void
113  std::vector<Tensor<1,dim> > &gradients) const
114  {
115  Assert (gradients.size() == n_components,
116  ExcDimensionMismatch(gradients.size(), n_components));
117  typename DH::active_cell_iterator cell = cell_hint.get();
118  if (cell == dh->end())
119  cell = dh->begin_active();
120 
121  boost::optional<Point<dim> >
122  qp = get_reference_coordinates (cell, p);
123  if (!qp)
124  {
125  const std::pair<typename DH::active_cell_iterator, Point<dim> > my_pair
126  = GridTools::find_active_cell_around_point (mapping, *dh, p);
127  AssertThrow (my_pair.first->is_locally_owned(),
128  ExcPointNotAvailableHere());
129 
130  cell = my_pair.first;
131  qp = my_pair.second;
132  }
133 
134  cell_hint.get() = cell;
135 
136  // Now we can find out about the point
137  Quadrature<dim> quad(qp.get());
138  FEValues<dim> fe_v(mapping, cell->get_fe(), quad,
140  fe_v.reinit(cell);
141  std::vector< std::vector<Tensor<1,dim> > > vgrads
142  (1, std::vector<Tensor<1,dim> >(n_components) );
143  fe_v.get_function_grads(data_vector, vgrads);
144  gradients = vgrads[0];
145  }
146 
147 
148 
149  template <int dim, typename DH, typename VECTOR>
153  const unsigned int comp) const
154  {
155  std::vector<Tensor<1,dim> > grads(n_components);
156  vector_gradient(p, grads);
157  return grads[comp];
158  }
159 
160 
161 
162  template <int dim, typename DH, typename VECTOR>
163  void
166  Vector<double> &values) const
167  {
168  Assert (values.size() == n_components,
169  ExcDimensionMismatch(values.size(), n_components));
170  typename DH::active_cell_iterator cell = cell_hint.get();
171  if (cell == dh->end())
172  cell = dh->begin_active();
173 
174  boost::optional<Point<dim> >
175  qp = get_reference_coordinates (cell, p);
176  if (!qp)
177  {
178  const std::pair<typename DH::active_cell_iterator, Point<dim> > my_pair
179  = GridTools::find_active_cell_around_point (mapping, *dh, p);
180  AssertThrow (my_pair.first->is_locally_owned(),
181  ExcPointNotAvailableHere());
182 
183  cell = my_pair.first;
184  qp = my_pair.second;
185  }
186 
187  cell_hint.get() = cell;
188 
189  // Now we can find out about the point
190  Quadrature<dim> quad(qp.get());
191  FEValues<dim> fe_v(mapping, cell->get_fe(), quad,
193  fe_v.reinit(cell);
194  std::vector< Vector<double> > vvalues (1, values);
195  fe_v.get_function_laplacians(data_vector, vvalues);
196  values = vvalues[0];
197  }
198 
199 
200 
201  template <int dim, typename DH, typename VECTOR>
203  (const Point<dim> &p, const unsigned int comp) const
204  {
205  Vector<double> lap(n_components);
206  vector_laplacian(p, lap);
207  return lap[comp];
208  }
209 
210 
211  // Now the list versions
212  // ==============================
213 
214  template <int dim, typename DH, typename VECTOR>
215  void
217  vector_value_list (const std::vector<Point< dim > > &points,
218  std::vector< Vector<double> > &values) const
219  {
220  Assert(points.size() == values.size(),
221  ExcDimensionMismatch(points.size(), values.size()));
222 
223  std::vector<typename DH::active_cell_iterator > cells;
224  std::vector<std::vector<Point<dim> > > qpoints;
225  std::vector<std::vector<unsigned int> > maps;
226 
227  unsigned int ncells = compute_point_locations(points, cells, qpoints, maps);
228  hp::MappingCollection<dim> mapping_collection (mapping);
229  hp::FECollection<dim> fe_collection (dh->get_fe ());
230  hp::QCollection<dim> quadrature_collection;
231  // Create quadrature collection
232  for (unsigned int i=0; i<ncells; ++i)
233  {
234  // Number of quadrature points on this cell
235  unsigned int nq = qpoints[i].size();
236  // Construct a quadrature formula
237  std::vector< double > ww(nq, 1./((double) nq));
238 
239  quadrature_collection.push_back (Quadrature<dim> (qpoints[i], ww));
240  }
241  // Get a function value object
242  hp::FEValues<dim> fe_v(mapping_collection, fe_collection, quadrature_collection,
243  update_values);
244  // Now gather all the informations we need
245  for (unsigned int i=0; i<ncells; ++i)
246  {
247  fe_v.reinit(cells[i], i, 0);
248  const unsigned int nq = qpoints[i].size();
249  std::vector< Vector<double> > vvalues (nq, Vector<double>(n_components));
250  fe_v.get_present_fe_values ().get_function_values(data_vector, vvalues);
251  for (unsigned int q=0; q<nq; ++q)
252  values[maps[i][q]] = vvalues[q];
253  }
254  }
255 
256 
257 
258  template <int dim, typename DH, typename VECTOR>
259  void
261  value_list (const std::vector<Point< dim > > &points,
262  std::vector< double > &values,
263  const unsigned int component) const
264  {
265  Assert(points.size() == values.size(),
266  ExcDimensionMismatch(points.size(), values.size()));
267  std::vector< Vector<double> > vvalues(points.size(), Vector<double>(n_components));
268  vector_value_list(points, vvalues);
269  for (unsigned int q=0; q<points.size(); ++q)
270  values[q] = vvalues[q](component);
271  }
272 
273 
274 
275  template <int dim, typename DH, typename VECTOR>
276  void
278  vector_gradient_list (const std::vector<Point< dim > > &points,
279  std::vector<
280  std::vector< Tensor<1,dim> > > &values) const
281  {
282  Assert(points.size() == values.size(),
283  ExcDimensionMismatch(points.size(), values.size()));
284 
285  std::vector<typename DH::active_cell_iterator > cells;
286  std::vector<std::vector<Point<dim> > > qpoints;
287  std::vector<std::vector<unsigned int> > maps;
288 
289  unsigned int ncells = compute_point_locations(points, cells, qpoints, maps);
290  hp::MappingCollection<dim> mapping_collection (mapping);
291  hp::FECollection<dim> fe_collection (dh->get_fe ());
292  hp::QCollection<dim> quadrature_collection;
293  // Create quadrature collection
294  for (unsigned int i=0; i<ncells; ++i)
295  {
296  // Number of quadrature points on this cell
297  unsigned int nq = qpoints[i].size();
298  // Construct a quadrature formula
299  std::vector< double > ww(nq, 1./((double) nq));
300 
301  quadrature_collection.push_back (Quadrature<dim> (qpoints[i], ww));
302  }
303  // Get a function value object
304  hp::FEValues<dim> fe_v(mapping_collection, fe_collection, quadrature_collection,
306  // Now gather all the informations we need
307  for (unsigned int i=0; i<ncells; ++i)
308  {
309  fe_v.reinit(cells[i], i, 0);
310  const unsigned int nq = qpoints[i].size();
311  std::vector< std::vector<Tensor<1,dim> > >
312  vgrads (nq, std::vector<Tensor<1,dim> >(n_components));
313  fe_v.get_present_fe_values ().get_function_grads(data_vector, vgrads);
314  for (unsigned int q=0; q<nq; ++q)
315  values[maps[i][q]] = vgrads[q];
316  }
317  }
318 
319  template <int dim, typename DH, typename VECTOR>
320  void
322  gradient_list (const std::vector<Point< dim > > &points,
323  std::vector< Tensor<1,dim> > &values,
324  const unsigned int component) const
325  {
326  Assert(points.size() == values.size(),
327  ExcDimensionMismatch(points.size(), values.size()));
328  std::vector< std::vector<Tensor<1,dim> > >
329  vvalues(points.size(), std::vector<Tensor<1,dim> >(n_components));
330  vector_gradient_list(points, vvalues);
331  for (unsigned int q=0; q<points.size(); ++q)
332  values[q] = vvalues[q][component];
333  }
334 
335 
336  template <int dim, typename DH, typename VECTOR>
337  void
339  vector_laplacian_list (const std::vector<Point< dim > > &points,
340  std::vector< Vector<double> > &values) const
341  {
342  Assert(points.size() == values.size(),
343  ExcDimensionMismatch(points.size(), values.size()));
344 
345  std::vector<typename DH::active_cell_iterator > cells;
346  std::vector<std::vector<Point<dim> > > qpoints;
347  std::vector<std::vector<unsigned int> > maps;
348 
349  unsigned int ncells = compute_point_locations(points, cells, qpoints, maps);
350  hp::MappingCollection<dim> mapping_collection (mapping);
351  hp::FECollection<dim> fe_collection (dh->get_fe ());
352  hp::QCollection<dim> quadrature_collection;
353  // Create quadrature collection
354  for (unsigned int i=0; i<ncells; ++i)
355  {
356  // Number of quadrature points on this cell
357  unsigned int nq = qpoints[i].size();
358  // Construct a quadrature formula
359  std::vector< double > ww(nq, 1./((double) nq));
360 
361  quadrature_collection.push_back (Quadrature<dim> (qpoints[i], ww));
362  }
363  // Get a function value object
364  hp::FEValues<dim> fe_v(mapping_collection, fe_collection, quadrature_collection,
366  // Now gather all the informations we need
367  for (unsigned int i=0; i<ncells; ++i)
368  {
369  fe_v.reinit(cells[i], i, 0);
370  const unsigned int nq = qpoints[i].size();
371  std::vector< Vector<double> > vvalues (nq, Vector<double>(n_components));
372  fe_v.get_present_fe_values ().get_function_laplacians(data_vector, vvalues);
373  for (unsigned int q=0; q<nq; ++q)
374  values[maps[i][q]] = vvalues[q];
375  }
376  }
377 
378  template <int dim, typename DH, typename VECTOR>
379  void
381  laplacian_list (const std::vector<Point< dim > > &points,
382  std::vector< double > &values,
383  const unsigned int component) const
384  {
385  Assert(points.size() == values.size(),
386  ExcDimensionMismatch(points.size(), values.size()));
387  std::vector< Vector<double> > vvalues(points.size(), Vector<double>(n_components));
388  vector_laplacian_list(points, vvalues);
389  for (unsigned int q=0; q<points.size(); ++q)
390  values[q] = vvalues[q](component);
391  }
392 
393 
394 
395  template <int dim, typename DH, typename VECTOR>
397  compute_point_locations(const std::vector<Point<dim> > &points,
398  std::vector<typename DH::active_cell_iterator > &cells,
399  std::vector<std::vector<Point<dim> > > &qpoints,
400  std::vector<std::vector<unsigned int> > &maps) const
401  {
402  // How many points are here?
403  const unsigned int np = points.size();
404 
405  // Reset output maps.
406  cells.clear();
407  qpoints.clear();
408  maps.clear();
409 
410  // Now the easy case.
411  if (np==0) return 0;
412 
413  // Keep track of the points we
414  // found
415  std::vector<bool> point_flags(np, false);
416 
417  // Set this to true until all
418  // points have been classified
419  bool left_over = true;
420 
421  // Current quadrature point
422  typename DH::active_cell_iterator cell = cell_hint.get();
423  if (cell == dh->end())
424  cell = dh->begin_active();
425 
426  {
427  // see if the point is
428  // inside the
429  // cell. there are two
430  // ways that
431  // transform_real_to_unit_cell
432  // can indicate that a
433  // point is outside: by
434  // returning
435  // coordinates that lie
436  // outside the
437  // reference cell, or
438  // by throwing an
439  // exception. handle
440  // both
441  boost::optional<Point<dim> >
442  qp = get_reference_coordinates (cell, points[0]);
443  if (!qp)
444  {
445  const std::pair<typename DH::active_cell_iterator, Point<dim> >
447  (mapping, *dh, points[0]);
448  AssertThrow (my_pair.first->is_locally_owned(),
449  ExcPointNotAvailableHere());
450 
451  cell = my_pair.first;
452  qp = my_pair.second;
453  point_flags[0] = true;
454  }
455 
456  // Put in the first point.
457  cells.push_back(cell);
458  qpoints.push_back(std::vector<Point<dim> >(1, qp.get()));
459  maps.push_back(std::vector<unsigned int> (1, 0));
460  }
461 
462 
463  // Check if we need to do anything else
464  if (points.size() > 1)
465  left_over = true;
466  else
467  left_over = false;
468 
469 
470  // This is the first index of a non processed point
471  unsigned int first_outside = 1;
472 
473  // And this is the index of the current cell
474  unsigned int c = 0;
475 
476  while (left_over == true)
477  {
478  // Assume this is the last one
479  left_over = false;
480  Assert(first_outside < np,
481  ExcIndexRange(first_outside, 0, np));
482 
483  // If we found one in this cell, keep looking in the same cell
484  for (unsigned int p=first_outside; p<np; ++p)
485  if (point_flags[p] == false)
486  {
487  // same logic as above
488  const boost::optional<Point<dim> >
489  qp = get_reference_coordinates (cells[c], points[p]);
490  if (qp)
491  {
492  point_flags[p] = true;
493  qpoints[c].push_back(qp.get());
494  maps[c].push_back(p);
495  }
496  else
497  {
498  // Set things up for next round
499  if (left_over == false)
500  first_outside = p;
501  left_over = true;
502  }
503  }
504  // If we got here and there is
505  // no left over, we are
506  // done. Else we need to find
507  // the next cell
508  if (left_over == true)
509  {
510  const std::pair<typename DH::active_cell_iterator, Point<dim> > my_pair
511  = GridTools::find_active_cell_around_point (mapping, *dh, points[first_outside]);
512  AssertThrow (my_pair.first->is_locally_owned(),
513  ExcPointNotAvailableHere());
514 
515  cells.push_back(my_pair.first);
516  qpoints.push_back(std::vector<Point<dim> >(1, my_pair.second));
517  maps.push_back(std::vector<unsigned int>(1, first_outside));
518  c++;
519  point_flags[first_outside] = true;
520  // And check if we can exit the loop now
521  if (first_outside == np-1)
522  left_over = false;
523  }
524  }
525 
526  // Augment of one the number of cells
527  ++c;
528  // Debug Checking
529  Assert(c == cells.size(), ExcInternalError());
530 
531  Assert(c == maps.size(),
532  ExcDimensionMismatch(c, maps.size()));
533 
534  Assert(c == qpoints.size(),
535  ExcDimensionMismatch(c, qpoints.size()));
536 
537 #ifdef DEBUG
538  unsigned int qps = 0;
539  // The number of points in all
540  // the cells must be the same as
541  // the number of points we
542  // started off from.
543  for (unsigned int n=0; n<c; ++n)
544  {
545  Assert(qpoints[n].size() == maps[n].size(),
546  ExcDimensionMismatch(qpoints[n].size(), maps[n].size()));
547  qps += qpoints[n].size();
548  }
549  Assert(qps == np,
550  ExcDimensionMismatch(qps, np));
551 #endif
552 
553  return c;
554  }
555 
556 
557  template <int dim, typename DH, typename VECTOR>
558  boost::optional<Point<dim> >
560  get_reference_coordinates (const typename DH::active_cell_iterator &cell,
561  const Point<dim> &point) const
562  {
563  try
564  {
565  Point<dim> qp = mapping.transform_real_to_unit_cell(cell, point);
567  return qp;
568  else
569  return boost::optional<Point<dim> >();
570  }
571  catch (const typename Mapping<dim>::ExcTransformationFailed &)
572  {
573  // transformation failed, so
574  // assume the point is
575  // outside
576  return boost::optional<Point<dim> >();
577  }
578  }
579 
580 }
581 
582 DEAL_II_NAMESPACE_CLOSE
boost::optional< Point< dim > > get_reference_coordinates(const typename DH::active_cell_iterator &cell, const Point< dim > &point) const
Shape function values.
virtual void vector_laplacian_list(const std::vector< Point< dim > > &points, std::vector< Vector< double > > &values) const
void reinit(const TriaIterator< DoFCellAccessor< DH, lda > > cell, const unsigned int q_index=numbers::invalid_unsigned_int, const unsigned int mapping_index=numbers::invalid_unsigned_int, const unsigned int fe_index=numbers::invalid_unsigned_int)
unsigned int compute_point_locations(const std::vector< Point< dim > > &points, std::vector< typename DH::active_cell_iterator > &cells, std::vector< std::vector< Point< dim > > > &qpoints, std::vector< std::vector< unsigned int > > &maps) const
FEFieldFunction(const DH &dh, const VECTOR &data_vector, const Mapping< dim > &mapping=StaticMappingQ1< dim >::mapping)
#define AssertThrow(cond, exc)
Definition: exceptions.h:362
virtual void vector_value_list(const std::vector< Point< dim > > &points, std::vector< Vector< double > > &values) const
virtual double value(const Point< dim > &p, const unsigned int component=0) const
#define Assert(cond, exc)
Definition: exceptions.h:299
const ::FEValues< dim, spacedim > & get_present_fe_values() const
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const
void reinit(const TriaIterator< DoFCellAccessor< DH, level_dof_access > > cell)
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
virtual void vector_gradient(const Point< dim > &p, std::vector< Tensor< 1, dim > > &gradients) const
virtual void vector_value(const Point< dim > &p, Vector< double > &values) const
Second derivatives of shape functions.
void push_back(const FiniteElement< dim, spacedim > &new_fe)
virtual void laplacian_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const
void set_active_cell(const typename DH::active_cell_iterator &newcell)
virtual void vector_gradient_list(const std::vector< Point< dim > > &p, std::vector< std::vector< Tensor< 1, dim > > > &gradients) const
Shape function gradients.
virtual void vector_laplacian(const Point< dim > &p, Vector< double > &values) const
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const
Container< dim, spacedim >::active_cell_iterator find_active_cell_around_point(const Container< dim, spacedim > &container, const Point< spacedim > &p)
::ExceptionBase & ExcInternalError()
std::size_t size() const
virtual void gradient_list(const std::vector< Point< dim > > &p, std::vector< Tensor< 1, dim > > &gradients, const unsigned int component=0) const