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>
29 DEAL_II_NAMESPACE_OPEN
34 template <
int dim,
typename DH,
typename VECTOR>
39 Function<dim>(mydh.get_fe().n_components()),
40 dh(&mydh,
"FEFieldFunction"),
44 n_components(mydh.get_fe().n_components())
50 template <
int dim,
typename DH,
typename VECTOR>
55 cell_hint.get() = newcell;
60 template <
int dim,
typename DH,
typename VECTOR>
66 typename DH::active_cell_iterator cell = cell_hint.get();
67 if (cell == dh->end())
68 cell = dh->begin_active();
70 boost::optional<Point<dim> >
71 qp = get_reference_coordinates (cell, p);
74 const std::pair<typename DH::active_cell_iterator, Point<dim> > my_pair
77 ExcPointNotAvailableHere());
83 cell_hint.get() = cell;
90 std::vector< Vector<double> > vvalues (1, values);
91 fe_v.get_function_values(data_vector, vvalues);
97 template <
int dim,
typename DH,
typename VECTOR>
100 const unsigned int comp)
const
103 vector_value(p, values);
109 template <
int dim,
typename DH,
typename VECTOR>
115 Assert (gradients.size() == n_components,
117 typename DH::active_cell_iterator cell = cell_hint.get();
118 if (cell == dh->end())
119 cell = dh->begin_active();
121 boost::optional<Point<dim> >
122 qp = get_reference_coordinates (cell, p);
125 const std::pair<typename DH::active_cell_iterator, Point<dim> > my_pair
128 ExcPointNotAvailableHere());
130 cell = my_pair.first;
134 cell_hint.get() = cell;
141 std::vector< std::vector<Tensor<1,dim> > > vgrads
143 fe_v.get_function_grads(data_vector, vgrads);
144 gradients = vgrads[0];
149 template <
int dim,
typename DH,
typename VECTOR>
153 const unsigned int comp)
const
155 std::vector<Tensor<1,dim> > grads(n_components);
156 vector_gradient(p, grads);
162 template <
int dim,
typename DH,
typename VECTOR>
170 typename DH::active_cell_iterator cell = cell_hint.get();
171 if (cell == dh->end())
172 cell = dh->begin_active();
174 boost::optional<Point<dim> >
175 qp = get_reference_coordinates (cell, p);
178 const std::pair<typename DH::active_cell_iterator, Point<dim> > my_pair
181 ExcPointNotAvailableHere());
183 cell = my_pair.first;
187 cell_hint.get() = cell;
194 std::vector< Vector<double> > vvalues (1, values);
195 fe_v.get_function_laplacians(data_vector, vvalues);
201 template <
int dim,
typename DH,
typename VECTOR>
206 vector_laplacian(p, lap);
214 template <
int dim,
typename DH,
typename VECTOR>
220 Assert(points.size() == values.size(),
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;
227 unsigned int ncells = compute_point_locations(points, cells, qpoints, maps);
232 for (
unsigned int i=0; i<ncells; ++i)
235 unsigned int nq = qpoints[i].size();
237 std::vector< double > ww(nq, 1./((
double) nq));
245 for (
unsigned int i=0; i<ncells; ++i)
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));
251 for (
unsigned int q=0; q<nq; ++q)
252 values[maps[i][q]] = vvalues[q];
258 template <
int dim,
typename DH,
typename VECTOR>
262 std::vector< double > &values,
263 const unsigned int component)
const
265 Assert(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);
275 template <
int dim,
typename DH,
typename VECTOR>
282 Assert(points.size() == values.size(),
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;
289 unsigned int ncells = compute_point_locations(points, cells, qpoints, maps);
294 for (
unsigned int i=0; i<ncells; ++i)
297 unsigned int nq = qpoints[i].size();
299 std::vector< double > ww(nq, 1./((
double) nq));
307 for (
unsigned int i=0; i<ncells; ++i)
309 fe_v.
reinit(cells[i], i, 0);
310 const unsigned int nq = qpoints[i].size();
311 std::vector< std::vector<Tensor<1,dim> > >
314 for (
unsigned int q=0; q<nq; ++q)
315 values[maps[i][q]] = vgrads[q];
319 template <
int dim,
typename DH,
typename VECTOR>
324 const unsigned int component)
const
326 Assert(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];
336 template <
int dim,
typename DH,
typename VECTOR>
342 Assert(points.size() == values.size(),
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;
349 unsigned int ncells = compute_point_locations(points, cells, qpoints, maps);
354 for (
unsigned int i=0; i<ncells; ++i)
357 unsigned int nq = qpoints[i].size();
359 std::vector< double > ww(nq, 1./((
double) nq));
367 for (
unsigned int i=0; i<ncells; ++i)
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));
373 for (
unsigned int q=0; q<nq; ++q)
374 values[maps[i][q]] = vvalues[q];
378 template <
int dim,
typename DH,
typename VECTOR>
382 std::vector< double > &values,
383 const unsigned int component)
const
385 Assert(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);
395 template <
int dim,
typename DH,
typename VECTOR>
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
403 const unsigned int np = points.size();
415 std::vector<bool> point_flags(np,
false);
419 bool left_over =
true;
422 typename DH::active_cell_iterator cell = cell_hint.get();
423 if (cell == dh->end())
424 cell = dh->begin_active();
441 boost::optional<Point<dim> >
442 qp = get_reference_coordinates (cell, points[0]);
445 const std::pair<typename DH::active_cell_iterator, Point<dim> >
447 (mapping, *dh, points[0]);
449 ExcPointNotAvailableHere());
451 cell = my_pair.first;
453 point_flags[0] =
true;
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));
464 if (points.size() > 1)
471 unsigned int first_outside = 1;
476 while (left_over ==
true)
480 Assert(first_outside < np,
484 for (
unsigned int p=first_outside; p<np; ++p)
485 if (point_flags[p] ==
false)
488 const boost::optional<Point<dim> >
489 qp = get_reference_coordinates (cells[c], points[p]);
492 point_flags[p] =
true;
493 qpoints[c].push_back(qp.get());
494 maps[c].push_back(p);
499 if (left_over ==
false)
508 if (left_over ==
true)
510 const std::pair<typename DH::active_cell_iterator, Point<dim> > my_pair
513 ExcPointNotAvailableHere());
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));
519 point_flags[first_outside] =
true;
521 if (first_outside == np-1)
534 Assert(c == qpoints.size(),
538 unsigned int qps = 0;
543 for (
unsigned int n=0; n<c; ++n)
545 Assert(qpoints[n].size() == maps[n].size(),
547 qps += qpoints[n].size();
557 template <
int dim,
typename DH,
typename VECTOR>
558 boost::optional<Point<dim> >
565 Point<dim> qp = mapping.transform_real_to_unit_cell(cell, point);
569 return boost::optional<Point<dim> >();
576 return boost::optional<Point<dim> >();
582 DEAL_II_NAMESPACE_CLOSE
boost::optional< Point< dim > > get_reference_coordinates(const typename DH::active_cell_iterator &cell, const Point< dim > &point) const
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)
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
void push_back(const FiniteElement< dim, spacedim > &new_fe)
#define Assert(cond, exc)
const ::FEValues< dim, spacedim > & get_present_fe_values() const
void reinit(const TriaIterator< DoFCellAccessor< DH, level_dof_access > > cell)
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const
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.
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
::ExceptionBase & ExcInternalError()
virtual void gradient_list(const std::vector< Point< dim > > &p, std::vector< Tensor< 1, dim > > &gradients, const unsigned int component=0) const