Reference documentation for deal.II version 8.1.0
vector_tools.templates.h
1 // ---------------------------------------------------------------------
2 // @f$Id: vector_tools.templates.h 31932 2013-12-08 02:15:54Z heister @f$
3 //
4 // Copyright (C) 2005 - 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 #ifndef _deal2__vector_tools_templates_h
19 #define _deal2__vector_tools_templates_h
20 
21 #include <deal.II/base/derivative_form.h>
22 #include <deal.II/base/function.h>
23 #include <deal.II/base/quadrature.h>
24 #include <deal.II/lac/vector.h>
25 #include <deal.II/lac/block_vector.h>
26 #include <deal.II/lac/parallel_vector.h>
27 #include <deal.II/lac/parallel_block_vector.h>
28 #include <deal.II/lac/petsc_vector.h>
29 #include <deal.II/lac/petsc_block_vector.h>
30 #include <deal.II/lac/trilinos_vector.h>
31 #include <deal.II/lac/trilinos_block_vector.h>
32 #include <deal.II/lac/sparse_matrix.h>
33 #include <deal.II/lac/compressed_sparsity_pattern.h>
34 #include <deal.II/lac/compressed_simple_sparsity_pattern.h>
35 #include <deal.II/lac/precondition.h>
36 #include <deal.II/lac/solver_cg.h>
37 #include <deal.II/lac/vector_memory.h>
38 #include <deal.II/lac/filtered_matrix.h>
39 #include <deal.II/lac/constraint_matrix.h>
40 #include <deal.II/grid/tria_iterator.h>
41 #include <deal.II/grid/tria_accessor.h>
42 #include <deal.II/grid/tria_boundary.h>
43 #include <deal.II/grid/grid_tools.h>
44 #include <deal.II/grid/intergrid_map.h>
45 #include <deal.II/dofs/dof_accessor.h>
46 #include <deal.II/dofs/dof_handler.h>
47 #include <deal.II/dofs/dof_tools.h>
48 #include <deal.II/fe/fe.h>
49 #include <deal.II/fe/fe_nedelec.h>
50 #include <deal.II/fe/fe_raviart_thomas.h>
51 #include <deal.II/fe/fe_system.h>
52 #include <deal.II/fe/fe_tools.h>
53 #include <deal.II/fe/fe_values.h>
54 #include <deal.II/fe/fe_nothing.h>
55 #include <deal.II/fe/mapping_q1.h>
56 #include <deal.II/hp/dof_handler.h>
57 #include <deal.II/hp/fe_values.h>
58 #include <deal.II/hp/mapping_collection.h>
59 #include <deal.II/hp/q_collection.h>
60 #include <deal.II/numerics/vector_tools.h>
61 #include <deal.II/numerics/matrix_tools.h>
62 
63 #include <deal.II/base/std_cxx1x/array.h>
64 #include <numeric>
65 #include <algorithm>
66 #include <vector>
67 #include <cmath>
68 #include <limits>
69 #include <set>
70 
72 
73 
74 namespace VectorTools
75 {
76 
77  template <class VECTOR, int dim, int spacedim, template <int,int> class DH>
78  void interpolate (const Mapping<dim,spacedim> &mapping,
79  const DH<dim,spacedim> &dof,
80  const Function<spacedim> &function,
81  VECTOR &vec)
82  {
83  Assert (dof.get_fe().n_components() == function.n_components,
84  ExcDimensionMismatch(dof.get_fe().n_components(),
85  function.n_components));
86 
87  const hp::FECollection<dim,spacedim> fe (dof.get_fe());
88  const unsigned int n_components = fe.n_components();
89  const bool fe_is_system = (n_components != 1);
90 
91  typename DH<dim,spacedim>::active_cell_iterator cell = dof.begin_active(),
92  endc = dof.end();
93 
94  // For FESystems many of the
95  // unit_support_points will appear
96  // multiple times, as a point may be
97  // unit_support_point for several of the
98  // components of the system. The following
99  // is rather complicated, but at least
100  // attempts to avoid evaluating the
101  // vectorfunction multiple times at the
102  // same point on a cell.
103  //
104  // note that we have to set up all of the
105  // following arrays for each of the
106  // elements in the FECollection (which
107  // means only once if this is for a regular
108  // DoFHandler)
109  std::vector<std::vector<Point<dim> > > unit_support_points (fe.size());
110  for (unsigned int fe_index=0; fe_index<fe.size(); ++fe_index)
111  {
112  unit_support_points[fe_index] = fe[fe_index].get_unit_support_points();
113  Assert (unit_support_points[fe_index].size() != 0,
114  ExcNonInterpolatingFE());
115  }
116 
117 
118  // Find the support points on a cell that
119  // are mentioned multiple times in
120  // unit_support_points. Mark the first
121  // representative of each support point
122  // mentioned multiple times by appending
123  // its dof index to dofs_of_rep_points.
124  // Each multiple point gets to know the dof
125  // index of its representative point by the
126  // dof_to_rep_dof_table.
127 
128  // the following vector collects all dofs i,
129  // 0<=i<fe.dofs_per_cell, for that
130  // unit_support_points[i]
131  // is a representative one. i.e.
132  // the following vector collects all rep dofs.
133  // the position of a rep dof within this vector
134  // is called rep index.
135  std::vector<std::vector<types::global_dof_index> > dofs_of_rep_points(fe.size());
136  // the following table converts a dof i
137  // to the rep index.
138  std::vector<std::vector<types::global_dof_index> > dof_to_rep_index_table(fe.size());
139 
140  std::vector<unsigned int> n_rep_points (fe.size(), 0);
141 
142  for (unsigned int fe_index=0; fe_index<fe.size(); ++fe_index)
143  {
144  for (unsigned int i=0; i<fe[fe_index].dofs_per_cell; ++i)
145  {
146  bool representative=true;
147  // the following loop is looped
148  // the other way round to get
149  // the minimal effort of
150  // O(fe.dofs_per_cell) for multiple
151  // support points that are placed
152  // one after the other.
153  for (unsigned int j=dofs_of_rep_points[fe_index].size(); j>0; --j)
154  if (unit_support_points[fe_index][i]
155  == unit_support_points[fe_index][dofs_of_rep_points[fe_index][j-1]])
156  {
157  dof_to_rep_index_table[fe_index].push_back(j-1);
158  representative=false;
159  break;
160  }
161 
162  if (representative)
163  {
164  // rep_index=dofs_of_rep_points.size()
165  dof_to_rep_index_table[fe_index].push_back(dofs_of_rep_points[fe_index].size());
166  // dofs_of_rep_points[rep_index]=i
167  dofs_of_rep_points[fe_index].push_back(i);
168  ++n_rep_points[fe_index];
169  }
170  }
171 
172  Assert(dofs_of_rep_points[fe_index].size()==n_rep_points[fe_index],
173  ExcInternalError());
174  Assert(dof_to_rep_index_table[fe_index].size()==fe[fe_index].dofs_per_cell,
175  ExcInternalError());
176  }
177 
178  const unsigned int max_rep_points = *std::max_element (n_rep_points.begin(),
179  n_rep_points.end());
180  std::vector<types::global_dof_index> dofs_on_cell (fe.max_dofs_per_cell());
181  std::vector<Point<spacedim> > rep_points (max_rep_points);
182 
183  // get space for the values of the
184  // function at the rep support points.
185  //
186  // have two versions, one for system fe
187  // and one for scalar ones, to take the
188  // more efficient one respectively
189  std::vector<std::vector<double> > function_values_scalar(fe.size());
190  std::vector<std::vector<Vector<double> > > function_values_system(fe.size());
191 
192  // Make a quadrature rule from support points
193  // to feed it into FEValues
194  hp::QCollection<dim> support_quadrature;
195  for (unsigned int fe_index=0; fe_index<fe.size(); ++fe_index)
196  support_quadrature.push_back (Quadrature<dim>(unit_support_points[fe_index]));
197 
198  // Transformed support points are computed by
199  // FEValues
200  hp::MappingCollection<dim,spacedim> mapping_collection (mapping);
201 
202  hp::FEValues<dim,spacedim> fe_values (mapping_collection,
203  fe, support_quadrature, update_quadrature_points);
204 
205  for (; cell!=endc; ++cell)
206  if (cell->is_locally_owned())
207  {
208  const unsigned int fe_index = cell->active_fe_index();
209 
210  // for each cell:
211  // get location of finite element
212  // support_points
213  fe_values.reinit(cell);
214  const std::vector<Point<spacedim> > &support_points =
215  fe_values.get_present_fe_values().get_quadrature_points();
216 
217  // pick out the representative
218  // support points
219  rep_points.resize (dofs_of_rep_points[fe_index].size());
220  for (unsigned int j=0; j<dofs_of_rep_points[fe_index].size(); ++j)
221  rep_points[j] = support_points[dofs_of_rep_points[fe_index][j]];
222 
223  // get indices of the dofs on this cell
224  dofs_on_cell.resize (fe[fe_index].dofs_per_cell);
225  cell->get_dof_indices (dofs_on_cell);
226 
227 
228  if (fe_is_system)
229  {
230  // get function values at
231  // these points. Here: get
232  // all components
233  function_values_system[fe_index]
234  .resize (n_rep_points[fe_index],
235  Vector<double> (fe[fe_index].n_components()));
236  function.vector_value_list (rep_points,
237  function_values_system[fe_index]);
238  // distribute the function
239  // values to the global
240  // vector
241  for (unsigned int i=0; i<fe[fe_index].dofs_per_cell; ++i)
242  {
243  const unsigned int component
244  = fe[fe_index].system_to_component_index(i).first;
245  const unsigned int rep_dof=dof_to_rep_index_table[fe_index][i];
246  vec(dofs_on_cell[i])
247  = function_values_system[fe_index][rep_dof](component);
248  }
249  }
250  else
251  {
252  // get first component only,
253  // which is the only component
254  // in the function anyway
255  function_values_scalar[fe_index].resize (n_rep_points[fe_index]);
256  function.value_list (rep_points,
257  function_values_scalar[fe_index],
258  0);
259  // distribute the function
260  // values to the global
261  // vector
262  for (unsigned int i=0; i<fe[fe_index].dofs_per_cell; ++i)
263  vec(dofs_on_cell[i])
264  = function_values_scalar[fe_index][dof_to_rep_index_table[fe_index][i]];
265  }
266  }
267  vec.compress(VectorOperation::insert);
268  }
269 
270 
271  template <class VECTOR, class DH>
272  void interpolate (const DH &dof,
273  const Function<DH::space_dimension> &function,
274  VECTOR &vec)
275  {
277  dof, function, vec);
278  }
279 
280 
281 
282 
283  template <int dim, class InVector, class OutVector, int spacedim>
284  void
286  const DoFHandler<dim,spacedim> &dof_2,
287  const FullMatrix<double> &transfer,
288  const InVector &data_1,
289  OutVector &data_2)
290  {
291  Vector<double> cell_data_1(dof_1.get_fe().dofs_per_cell);
292  Vector<double> cell_data_2(dof_2.get_fe().dofs_per_cell);
293 
294  std::vector<short unsigned int> touch_count (dof_2.n_dofs(), 0); //TODO: check on datatype... kinda strange (UK)
295  std::vector<types::global_dof_index> local_dof_indices (dof_2.get_fe().dofs_per_cell);
296 
297  typename DoFHandler<dim,spacedim>::active_cell_iterator h = dof_1.begin_active();
298  typename DoFHandler<dim,spacedim>::active_cell_iterator l = dof_2.begin_active();
299  const typename DoFHandler<dim,spacedim>::cell_iterator endh = dof_1.end();
300 
301  for (; h != endh; ++h, ++l)
302  {
303  h->get_dof_values(data_1, cell_data_1);
304  transfer.vmult(cell_data_2, cell_data_1);
305 
306  l->get_dof_indices (local_dof_indices);
307 
308  // distribute cell vector
309  for (unsigned int j=0; j<dof_2.get_fe().dofs_per_cell; ++j)
310  {
311  data_2(local_dof_indices[j]) += cell_data_2(j);
312 
313  // count, how often we have
314  // added to this dof
315  Assert (touch_count[local_dof_indices[j]] < numbers::internal_face_boundary_id,
316  ExcInternalError());
317  ++touch_count[local_dof_indices[j]];
318  }
319  }
320 
321  // compute the mean value of the
322  // sum which we have placed in each
323  // entry of the output vector
324  for (unsigned int i=0; i<dof_2.n_dofs(); ++i)
325  {
326  Assert (touch_count[i] != 0,
327  ExcInternalError());
328 
329  data_2(i) /= touch_count[i];
330  }
331  }
332 
333 
334  namespace internal
335  {
341  template <class DH>
342  void
343  interpolate_zero_boundary_values (const DH &dof_handler,
344  std::map<types::global_dof_index,double> &boundary_values)
345  {
346  const unsigned int dim = DH::dimension;
347 
348  // loop over all boundary faces
349  // to get all dof indices of
350  // dofs on the boundary. note
351  // that in 3d there are cases
352  // where a face is not at the
353  // boundary, yet one of its
354  // lines is, and we should
355  // consider the degrees of
356  // freedom on it as boundary
357  // nodes. likewise, in 2d and
358  // 3d there are cases where a
359  // cell is only at the boundary
360  // by one vertex. nevertheless,
361  // since we do not support
362  // boundaries with dimension
363  // less or equal to dim-2, each
364  // such boundary dof is also
365  // found from some other face
366  // that is actually wholly on
367  // the boundary, not only by
368  // one line or one vertex
369  typename DH::active_cell_iterator
370  cell = dof_handler.begin_active(),
371  endc = dof_handler.end();
372  std::vector<types::global_dof_index> face_dof_indices;
373  for (; cell!=endc; ++cell)
374  for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
375  if (cell->at_boundary(f))
376  {
377  face_dof_indices.resize (cell->get_fe().dofs_per_face);
378  cell->face(f)->get_dof_indices (face_dof_indices,
379  cell->active_fe_index());
380  for (unsigned int i=0; i<cell->get_fe().dofs_per_face; ++i)
381  // enter zero boundary values
382  // for all boundary nodes
383  //
384  // we need not care about
385  // vector valued elements here,
386  // since we set all components
387  boundary_values[face_dof_indices[i]] = 0.;
388  }
389  }
390  }
391 
392 
393 
394  template <int dim, int spacedim,
395  template <int,int> class DH,
396  class Vector>
397  void
398  interpolate_to_different_mesh (const DH<dim, spacedim> &dof1,
399  const Vector &u1,
400  const DH<dim, spacedim> &dof2,
401  Vector &u2)
402  {
404  ExcMessage ("The two containers must represent triangulations that "
405  "have the same coarse meshes"));
406 
407  InterGridMap<DH<dim, spacedim> > intergridmap;
408  intergridmap.make_mapping(dof1, dof2);
409 
410  ConstraintMatrix dummy;
411  dummy.close();
412 
413  interpolate_to_different_mesh(intergridmap, u1, dummy, u2);
414  }
415 
416 
417 
418  template <int dim, int spacedim,
419  template <int,int> class DH,
420  class Vector>
421  void
422  interpolate_to_different_mesh (const DH<dim, spacedim> &dof1,
423  const Vector &u1,
424  const DH<dim, spacedim> &dof2,
425  const ConstraintMatrix &constraints,
426  Vector &u2)
427  {
429  ExcMessage ("The two containers must represent triangulations that "
430  "have the same coarse meshes"));
431 
432  InterGridMap<DH<dim, spacedim> > intergridmap;
433  intergridmap.make_mapping(dof1, dof2);
434 
435  interpolate_to_different_mesh(intergridmap, u1, constraints, u2);
436  }
437 
438 
439 
440  template <int dim, int spacedim,
441  template <int,int> class DH,
442  class Vector>
443  void
444  interpolate_to_different_mesh (const InterGridMap<DH<dim, spacedim> > &intergridmap,
445  const Vector &u1,
446  const ConstraintMatrix &constraints,
447  Vector &u2)
448  {
449  const DH<dim, spacedim> &dof1 = intergridmap.get_source_grid();
450  const DH<dim, spacedim> &dof2 = intergridmap.get_destination_grid();
451 
452  Assert(u1.size()==dof1.n_dofs(),
453  ExcDimensionMismatch(u1.size(), dof1.n_dofs()));
454  Assert(u2.size()==dof2.n_dofs(),
455  ExcDimensionMismatch(u2.size(),dof2.n_dofs()));
456 
457  Vector cache;
458 
459  // Looping over the finest common
460  // mesh, this means that source and
461  // destination cells have to be on the
462  // same level and at least one has to
463  // be active.
464  //
465  // Therefor, loop over all cells
466  // (active and inactive) of the source
467  // grid ..
468  typename DH<dim,spacedim>::cell_iterator cell1 = dof1.begin();
469  const typename DH<dim,spacedim>::cell_iterator endc1 = dof1.end();
470 
471  for (; cell1 != endc1; ++cell1)
472  {
473  const typename DH<dim,spacedim>::cell_iterator cell2 = intergridmap[cell1];
474 
475  // .. and skip if source and destination
476  // cells are not on the same level ..
477  if (cell1->level() != cell2->level())
478  continue;
479  // .. or none of them is active.
480  if (!cell1->active() && !cell2->active())
481  continue;
482 
483  Assert(cell1->get_fe().get_name() ==
484  cell2->get_fe().get_name(),
485  ExcMessage ("Source and destination cells need to use the same finite element"));
486 
487  cache.reinit(cell1->get_fe().dofs_per_cell);
488 
489  // Get and set the corresponding
490  // dof_values by interpolation.
491  cell1->get_interpolated_dof_values(u1, cache);
492  cell2->set_dof_values_by_interpolation(cache, u2);
493  }
494 
495  // Apply hanging node constraints.
496  constraints.distribute(u2);
497  }
498 
499 
500  namespace
501  {
505  template <int dim, int spacedim,
506  template <int,int> class DH,
507  template <int,int> class M_or_MC,
508  template <int> class Q_or_QC>
509  void project_compute_b_v (const M_or_MC<dim, spacedim> &mapping,
510  const DH<dim,spacedim> &dof,
511  const Function<spacedim> &function,
512  const bool enforce_zero_boundary,
513  const Q_or_QC<dim-1> &q_boundary,
514  const bool project_to_boundary_first,
515  std::map<types::global_dof_index,double> &boundary_values)
516  {
517  if (enforce_zero_boundary == true)
518  // no need to project boundary
519  // values, but enforce
520  // homogeneous boundary values
521  // anyway
522  internal::
523  interpolate_zero_boundary_values (dof, boundary_values);
524 
525  else
526  // no homogeneous boundary values
527  if (project_to_boundary_first == true)
528  // boundary projection required
529  {
530  // set up a list of boundary
531  // functions for the
532  // different boundary
533  // parts. We want the
534  // function to hold on
535  // all parts of the boundary
536  const std::vector<types::boundary_id>
537  used_boundary_indicators = dof.get_tria().get_boundary_indicators();
538 
539  typename FunctionMap<spacedim>::type boundary_functions;
540  for (unsigned int i=0; i<used_boundary_indicators.size(); ++i)
541  boundary_functions[used_boundary_indicators[i]] = &function;
542  project_boundary_values (mapping, dof, boundary_functions, q_boundary,
543  boundary_values);
544  }
545  }
546 
547 
552  bool constraints_and_b_v_are_compatible (const ConstraintMatrix &constraints,
553  std::map<types::global_dof_index,double> &boundary_values)
554  {
555  for (std::map<types::global_dof_index,double>::iterator it=boundary_values.begin();
556  it != boundary_values.end(); ++it)
557  if (constraints.is_constrained(it->first))
558 //TODO: This looks wrong -- shouldn't it be ==0 in the first condition and && ?
559  if (!(constraints.get_constraint_entries(it->first)->size() > 0
560  ||
561  (constraints.get_inhomogeneity(it->first) == it->second)))
562  return false;
563 
564  return true;
565  }
566 
567 
571  template <int dim, int spacedim,
572  class Vector,
573  template <int,int> class DH,
574  template <int,int> class M_or_MC,
575  template <int> class Q_or_QC>
576  void do_project (const M_or_MC<dim, spacedim> &mapping,
577  const DH<dim,spacedim> &dof,
578  const ConstraintMatrix &constraints,
579  const Q_or_QC<dim> &quadrature,
580  const Function<spacedim> &function,
581  Vector &vec_result,
582  const bool enforce_zero_boundary,
583  const Q_or_QC<dim-1> &q_boundary,
584  const bool project_to_boundary_first)
585  {
586  Assert (dof.get_fe().n_components() == function.n_components,
587  ExcDimensionMismatch(dof.get_fe().n_components(),
588  function.n_components));
589  Assert (vec_result.size() == dof.n_dofs(),
590  ExcDimensionMismatch (vec_result.size(), dof.n_dofs()));
591 
592  // make up boundary values
593  std::map<types::global_dof_index,double> boundary_values;
594  project_compute_b_v(mapping, dof, function, enforce_zero_boundary,
595  q_boundary, project_to_boundary_first, boundary_values);
596 
597  // check if constraints are compatible (see below)
598  const bool constraints_are_compatible =
599  constraints_and_b_v_are_compatible(constraints, boundary_values);
600 
601  // set up mass matrix and right hand side
602  ::Vector<double> vec (dof.n_dofs());
603  SparsityPattern sparsity;
604  {
605  CompressedSimpleSparsityPattern csp (dof.n_dofs(), dof.n_dofs());
606  DoFTools::make_sparsity_pattern (dof, csp, constraints,
607  !constraints_are_compatible);
608 
609  sparsity.copy_from (csp);
610  }
611  SparseMatrix<double> mass_matrix (sparsity);
612  ::Vector<double> tmp (mass_matrix.n());
613 
614  // If the constraint matrix does not conflict with the given boundary
615  // values (i.e., it either does not contain boundary values or it contains
616  // the same as boundary_values), we can let it call
617  // distribute_local_to_global straight away, otherwise we need to first
618  // interpolate the boundary values and then condense the matrix and vector
619  if (constraints_are_compatible)
620  {
621  const Function<spacedim> *dummy = 0;
622  MatrixCreator::create_mass_matrix (mapping, dof, quadrature,
623  mass_matrix, function, tmp,
624  dummy, constraints);
625  if (boundary_values.size() > 0)
626  MatrixTools::apply_boundary_values (boundary_values,
627  mass_matrix, vec, tmp,
628  true);
629  }
630  else
631  {
632  // create mass matrix and rhs at once, which is faster.
633  MatrixCreator::create_mass_matrix (mapping, dof, quadrature,
634  mass_matrix, function, tmp);
635  MatrixTools::apply_boundary_values (boundary_values,
636  mass_matrix, vec, tmp,
637  true);
638  constraints.condense(mass_matrix, tmp);
639  }
640 
641  // Allow for a maximum of 5*n steps to reduce the residual by 10^-12. n
642  // steps may not be sufficient, since roundoff errors may accumulate for
643  // badly conditioned matrices
644  ReductionControl control(5*tmp.size(), 0., 1e-12, false, false);
645  GrowingVectorMemory<> memory;
646  SolverCG<> cg(control,memory);
647 
648  PreconditionSSOR<> prec;
649  prec.initialize(mass_matrix, 1.2);
650 
651  cg.solve (mass_matrix, vec, tmp, prec);
652 
653  constraints.distribute (vec);
654 
655  // copy vec into vec_result. we can't use vec_result itself above, since
656  // it may be of another type than Vector<double> and that wouldn't
657  // necessarily go together with the matrix and other functions
658  for (unsigned int i=0; i<vec.size(); ++i)
659  vec_result(i) = vec(i);
660  }
661  }
662 
663 
664  template <int dim, class VECTOR, int spacedim>
665  void project (const Mapping<dim, spacedim> &mapping,
666  const DoFHandler<dim,spacedim> &dof,
667  const ConstraintMatrix &constraints,
668  const Quadrature<dim> &quadrature,
669  const Function<spacedim> &function,
670  VECTOR &vec_result,
671  const bool enforce_zero_boundary,
672  const Quadrature<dim-1> &q_boundary,
673  const bool project_to_boundary_first)
674  {
675  do_project (mapping, dof, constraints, quadrature,
676  function, vec_result,
677  enforce_zero_boundary, q_boundary,
678  project_to_boundary_first);
679  }
680 
681 
682  template <int dim, class VECTOR, int spacedim>
684  const ConstraintMatrix &constraints,
685  const Quadrature<dim> &quadrature,
686  const Function<spacedim> &function,
687  VECTOR &vec,
688  const bool enforce_zero_boundary,
689  const Quadrature<dim-1> &q_boundary,
690  const bool project_to_boundary_first)
691  {
692  project(StaticMappingQ1<dim,spacedim>::mapping, dof, constraints, quadrature, function, vec,
693  enforce_zero_boundary, q_boundary, project_to_boundary_first);
694  }
695 
696 
697 
698  template <int dim, class VECTOR, int spacedim>
700  const hp::DoFHandler<dim,spacedim> &dof,
701  const ConstraintMatrix &constraints,
702  const hp::QCollection<dim> &quadrature,
703  const Function<spacedim> &function,
704  VECTOR &vec_result,
705  const bool enforce_zero_boundary,
706  const hp::QCollection<dim-1> &q_boundary,
707  const bool project_to_boundary_first)
708  {
709  do_project (mapping, dof, constraints, quadrature,
710  function, vec_result,
711  enforce_zero_boundary, q_boundary,
712  project_to_boundary_first);
713  }
714 
715 
716  template <int dim, class VECTOR, int spacedim>
718  const ConstraintMatrix &constraints,
719  const hp::QCollection<dim> &quadrature,
720  const Function<spacedim> &function,
721  VECTOR &vec,
722  const bool enforce_zero_boundary,
723  const hp::QCollection<dim-1> &q_boundary,
724  const bool project_to_boundary_first)
725  {
727  dof, constraints, quadrature, function, vec,
728  enforce_zero_boundary, q_boundary, project_to_boundary_first);
729  }
730 
731 
732  template <int dim, int spacedim>
734  const DoFHandler<dim,spacedim> &dof_handler,
735  const Quadrature<dim> &quadrature,
736  const Function<spacedim> &rhs_function,
737  Vector<double> &rhs_vector)
738  {
739  const FiniteElement<dim,spacedim> &fe = dof_handler.get_fe();
740  Assert (fe.n_components() == rhs_function.n_components,
741  ExcDimensionMismatch(fe.n_components(), rhs_function.n_components));
742  Assert (rhs_vector.size() == dof_handler.n_dofs(),
743  ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs()));
744  rhs_vector = 0;
745 
746  UpdateFlags update_flags = UpdateFlags(update_values |
749  FEValues<dim,spacedim> fe_values (mapping, fe, quadrature, update_flags);
750 
751  const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
752  n_q_points = fe_values.n_quadrature_points,
753  n_components = fe.n_components();
754 
755  std::vector<types::global_dof_index> dofs (dofs_per_cell);
756  Vector<double> cell_vector (dofs_per_cell);
757 
758  typename DoFHandler<dim,spacedim>::active_cell_iterator
759  cell = dof_handler.begin_active(),
760  endc = dof_handler.end();
761 
762  if (n_components==1)
763  {
764  std::vector<double> rhs_values(n_q_points);
765 
766  for (; cell!=endc; ++cell)
767  {
768  fe_values.reinit(cell);
769 
770  const std::vector<double> &weights = fe_values.get_JxW_values ();
771  rhs_function.value_list (fe_values.get_quadrature_points(),
772  rhs_values);
773 
774  cell_vector = 0;
775  for (unsigned int point=0; point<n_q_points; ++point)
776  for (unsigned int i=0; i<dofs_per_cell; ++i)
777  cell_vector(i) += rhs_values[point] *
778  fe_values.shape_value(i,point) *
779  weights[point];
780 
781  cell->get_dof_indices (dofs);
782 
783  for (unsigned int i=0; i<dofs_per_cell; ++i)
784  rhs_vector(dofs[i]) += cell_vector(i);
785  }
786 
787  }
788  else
789  {
790  std::vector<Vector<double> > rhs_values(n_q_points,
791  Vector<double>(n_components));
792 
793  for (; cell!=endc; ++cell)
794  {
795  fe_values.reinit(cell);
796 
797  const std::vector<double> &weights = fe_values.get_JxW_values ();
798  rhs_function.vector_value_list (fe_values.get_quadrature_points(),
799  rhs_values);
800 
801  cell_vector = 0;
802  // Use the faster code if the
803  // FiniteElement is primitive
804  if (fe.is_primitive ())
805  {
806  for (unsigned int point=0; point<n_q_points; ++point)
807  for (unsigned int i=0; i<dofs_per_cell; ++i)
808  {
809  const unsigned int component
810  = fe.system_to_component_index(i).first;
811 
812  cell_vector(i) += rhs_values[point](component) *
813  fe_values.shape_value(i,point) *
814  weights[point];
815  }
816  }
817  else
818  {
819  // Otherwise do it the way
820  // proposed for vector valued
821  // elements
822  for (unsigned int point=0; point<n_q_points; ++point)
823  for (unsigned int i=0; i<dofs_per_cell; ++i)
824  for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
825  if (fe.get_nonzero_components(i)[comp_i])
826  {
827  cell_vector(i) += rhs_values[point](comp_i) *
828  fe_values.shape_value_component(i,point,comp_i) *
829  weights[point];
830  }
831  }
832 
833  cell->get_dof_indices (dofs);
834 
835  for (unsigned int i=0; i<dofs_per_cell; ++i)
836  rhs_vector(dofs[i]) += cell_vector(i);
837  }
838  }
839  }
840 
841 
842 
843  template <int dim, int spacedim>
845  const Quadrature<dim> &quadrature,
846  const Function<spacedim> &rhs_function,
847  Vector<double> &rhs_vector)
848  {
850  rhs_function, rhs_vector);
851  }
852 
853 
854 
855 
856  template <int dim, int spacedim>
858  const hp::DoFHandler<dim,spacedim> &dof_handler,
859  const hp::QCollection<dim> &quadrature,
860  const Function<spacedim> &rhs_function,
861  Vector<double> &rhs_vector)
862  {
863  const hp::FECollection<dim,spacedim> &fe = dof_handler.get_fe();
864  Assert (fe.n_components() == rhs_function.n_components,
865  ExcDimensionMismatch(fe.n_components(), rhs_function.n_components));
866  Assert (rhs_vector.size() == dof_handler.n_dofs(),
867  ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs()));
868  rhs_vector = 0;
869 
870  UpdateFlags update_flags = UpdateFlags(update_values |
873  hp::FEValues<dim,spacedim> x_fe_values (mapping, fe, quadrature, update_flags);
874 
875  const unsigned int n_components = fe.n_components();
876 
877  std::vector<types::global_dof_index> dofs (fe.max_dofs_per_cell());
878  Vector<double> cell_vector (fe.max_dofs_per_cell());
879 
880  typename hp::DoFHandler<dim,spacedim>::active_cell_iterator
881  cell = dof_handler.begin_active(),
882  endc = dof_handler.end();
883 
884  if (n_components==1)
885  {
886  std::vector<double> rhs_values;
887 
888  for (; cell!=endc; ++cell)
889  {
890  x_fe_values.reinit(cell);
891 
892  const FEValues<dim,spacedim> &fe_values = x_fe_values.get_present_fe_values();
893 
894  const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
895  n_q_points = fe_values.n_quadrature_points;
896  rhs_values.resize (n_q_points);
897  dofs.resize (dofs_per_cell);
898  cell_vector.reinit (dofs_per_cell);
899 
900  const std::vector<double> &weights = fe_values.get_JxW_values ();
901  rhs_function.value_list (fe_values.get_quadrature_points(),
902  rhs_values);
903 
904  cell_vector = 0;
905  for (unsigned int point=0; point<n_q_points; ++point)
906  for (unsigned int i=0; i<dofs_per_cell; ++i)
907  cell_vector(i) += rhs_values[point] *
908  fe_values.shape_value(i,point) *
909  weights[point];
910 
911  cell->get_dof_indices (dofs);
912 
913  for (unsigned int i=0; i<dofs_per_cell; ++i)
914  rhs_vector(dofs[i]) += cell_vector(i);
915  }
916 
917  }
918  else
919  {
920  std::vector<Vector<double> > rhs_values;
921 
922  for (; cell!=endc; ++cell)
923  {
924  x_fe_values.reinit(cell);
925 
926  const FEValues<dim,spacedim> &fe_values = x_fe_values.get_present_fe_values();
927 
928  const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
929  n_q_points = fe_values.n_quadrature_points;
930  rhs_values.resize (n_q_points,
931  Vector<double>(n_components));
932  dofs.resize (dofs_per_cell);
933  cell_vector.reinit (dofs_per_cell);
934 
935  const std::vector<double> &weights = fe_values.get_JxW_values ();
936  rhs_function.vector_value_list (fe_values.get_quadrature_points(),
937  rhs_values);
938 
939  cell_vector = 0;
940 
941  // Use the faster code if the
942  // FiniteElement is primitive
943  if (cell->get_fe().is_primitive ())
944  {
945  for (unsigned int point=0; point<n_q_points; ++point)
946  for (unsigned int i=0; i<dofs_per_cell; ++i)
947  {
948  const unsigned int component
949  = cell->get_fe().system_to_component_index(i).first;
950 
951  cell_vector(i) += rhs_values[point](component) *
952  fe_values.shape_value(i,point) *
953  weights[point];
954  }
955  }
956  else
957  {
958  // Otherwise do it the way proposed
959  // for vector valued elements
960  for (unsigned int point=0; point<n_q_points; ++point)
961  for (unsigned int i=0; i<dofs_per_cell; ++i)
962  for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
963  if (cell->get_fe().get_nonzero_components(i)[comp_i])
964  {
965  cell_vector(i) += rhs_values[point](comp_i) *
966  fe_values.shape_value_component(i,point,comp_i) *
967  weights[point];
968  }
969  }
970 
971  cell->get_dof_indices (dofs);
972 
973  for (unsigned int i=0; i<dofs_per_cell; ++i)
974  rhs_vector(dofs[i]) += cell_vector(i);
975  }
976  }
977  }
978 
979 
980 
981  template <int dim, int spacedim>
983  const hp::QCollection<dim> &quadrature,
984  const Function<spacedim> &rhs_function,
985  Vector<double> &rhs_vector)
986  {
988  dof_handler, quadrature,
989  rhs_function, rhs_vector);
990  }
991 
992 
993 
994 
995  template <int dim, int spacedim>
997  const DoFHandler<dim,spacedim> &dof_handler,
998  const Point<spacedim> &p,
999  Vector<double> &rhs_vector)
1000  {
1001  Assert (rhs_vector.size() == dof_handler.n_dofs(),
1002  ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs()));
1003  Assert (dof_handler.get_fe().n_components() == 1,
1004  ExcMessage ("This function only works for scalar finite elements"));
1005 
1006  rhs_vector = 0;
1007 
1008  std::pair<typename DoFHandler<dim,spacedim>::active_cell_iterator, Point<spacedim> >
1009  cell_point =
1010  GridTools::find_active_cell_around_point (mapping, dof_handler, p);
1011 
1013 
1014  FEValues<dim,spacedim> fe_values(mapping, dof_handler.get_fe(),
1016  fe_values.reinit(cell_point.first);
1017 
1018  const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
1019 
1020  std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1021  cell_point.first->get_dof_indices (local_dof_indices);
1022 
1023  for (unsigned int i=0; i<dofs_per_cell; i++)
1024  rhs_vector(local_dof_indices[i]) = fe_values.shape_value(i,0);
1025  }
1026 
1027 
1028 
1029  template <int dim, int spacedim>
1031  const Point<spacedim> &p,
1032  Vector<double> &rhs_vector)
1033  {
1035  p, rhs_vector);
1036  }
1037 
1038 
1039  template <int dim, int spacedim>
1041  const hp::DoFHandler<dim,spacedim> &dof_handler,
1042  const Point<spacedim> &p,
1043  Vector<double> &rhs_vector)
1044  {
1045  Assert (rhs_vector.size() == dof_handler.n_dofs(),
1046  ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs()));
1047  Assert (dof_handler.get_fe().n_components() == 1,
1048  ExcMessage ("This function only works for scalar finite elements"));
1049 
1050  rhs_vector = 0;
1051 
1052  std::pair<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator, Point<spacedim> >
1053  cell_point =
1054  GridTools::find_active_cell_around_point (mapping, dof_handler, p);
1055 
1057 
1058  FEValues<dim> fe_values(mapping[cell_point.first->active_fe_index()],
1059  cell_point.first->get_fe(), q, UpdateFlags(update_values));
1060  fe_values.reinit(cell_point.first);
1061 
1062  const unsigned int dofs_per_cell = cell_point.first->get_fe().dofs_per_cell;
1063 
1064  std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1065  cell_point.first->get_dof_indices (local_dof_indices);
1066 
1067  for (unsigned int i=0; i<dofs_per_cell; i++)
1068  rhs_vector(local_dof_indices[i]) = fe_values.shape_value(i,0);
1069  }
1070 
1071 
1072 
1073  template <int dim, int spacedim>
1075  const Point<spacedim> &p,
1076  Vector<double> &rhs_vector)
1077  {
1079  dof_handler,
1080  p, rhs_vector);
1081  }
1082 
1083 
1084 
1085 
1086  template <int dim, int spacedim>
1088  const DoFHandler<dim,spacedim> &dof_handler,
1089  const Point<spacedim> &p,
1090  const Point<dim> &orientation,
1091  Vector<double> &rhs_vector)
1092  {
1093  Assert (rhs_vector.size() == dof_handler.n_dofs(),
1094  ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs()));
1095  Assert (dof_handler.get_fe().n_components() == dim,
1096  ExcMessage ("This function only works for vector-valued finite elements."));
1097 
1098  rhs_vector = 0;
1099 
1100  std::pair<typename DoFHandler<dim,spacedim>::active_cell_iterator, Point<spacedim> >
1101  cell_point =
1102  GridTools::find_active_cell_around_point (mapping, dof_handler, p);
1103 
1105 
1106  const FEValuesExtractors::Vector vec (0);
1107  FEValues<dim,spacedim> fe_values(mapping, dof_handler.get_fe(),
1109  fe_values.reinit(cell_point.first);
1110 
1111  const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
1112 
1113  std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1114  cell_point.first->get_dof_indices (local_dof_indices);
1115 
1116  for (unsigned int i=0; i<dofs_per_cell; i++)
1117  rhs_vector(local_dof_indices[i]) = orientation * fe_values[vec].value(i,0);
1118  }
1119 
1120 
1121 
1122  template <int dim, int spacedim>
1124  const Point<spacedim> &p,
1125  const Point<dim> &orientation,
1126  Vector<double> &rhs_vector)
1127  {
1129  p, orientation, rhs_vector);
1130  }
1131 
1132 
1133  template <int dim, int spacedim>
1135  const hp::DoFHandler<dim,spacedim> &dof_handler,
1136  const Point<spacedim> &p,
1137  const Point<dim> &orientation,
1138  Vector<double> &rhs_vector)
1139  {
1140  Assert (rhs_vector.size() == dof_handler.n_dofs(),
1141  ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs()));
1142  Assert (dof_handler.get_fe().n_components() == dim,
1143  ExcMessage ("This function only works for vector-valued finite elements."));
1144 
1145  rhs_vector = 0;
1146 
1147  std::pair<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator, Point<spacedim> >
1148  cell_point =
1149  GridTools::find_active_cell_around_point (mapping, dof_handler, p);
1150 
1152 
1153  const FEValuesExtractors::Vector vec (0);
1154  FEValues<dim> fe_values(mapping[cell_point.first->active_fe_index()],
1155  cell_point.first->get_fe(), q, UpdateFlags(update_values));
1156  fe_values.reinit(cell_point.first);
1157 
1158  const unsigned int dofs_per_cell = cell_point.first->get_fe().dofs_per_cell;
1159 
1160  std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1161  cell_point.first->get_dof_indices (local_dof_indices);
1162 
1163  for (unsigned int i=0; i<dofs_per_cell; i++)
1164  rhs_vector(local_dof_indices[i]) = orientation * fe_values[vec].value(i,0);
1165  }
1166 
1167 
1168 
1169  template <int dim, int spacedim>
1171  const Point<spacedim> &p,
1172  const Point<dim> &orientation,
1173  Vector<double> &rhs_vector)
1174  {
1176  dof_handler,
1177  p, orientation, rhs_vector);
1178  }
1179 
1180 
1181 
1182  template <int dim, int spacedim>
1183  void
1185  const DoFHandler<dim,spacedim> &dof_handler,
1186  const Quadrature<dim-1> &quadrature,
1187  const Function<spacedim> &rhs_function,
1188  Vector<double> &rhs_vector,
1189  const std::set<types::boundary_id> &boundary_indicators)
1190  {
1191  const FiniteElement<dim> &fe = dof_handler.get_fe();
1192  Assert (fe.n_components() == rhs_function.n_components,
1193  ExcDimensionMismatch(fe.n_components(), rhs_function.n_components));
1194  Assert (rhs_vector.size() == dof_handler.n_dofs(),
1195  ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs()));
1196 
1197  rhs_vector = 0;
1198 
1199  UpdateFlags update_flags = UpdateFlags(update_values |
1202  FEFaceValues<dim> fe_values (mapping, fe, quadrature, update_flags);
1203 
1204  const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
1205  n_q_points = fe_values.n_quadrature_points,
1206  n_components = fe.n_components();
1207 
1208  std::vector<types::global_dof_index> dofs (dofs_per_cell);
1209  Vector<double> cell_vector (dofs_per_cell);
1210 
1211  typename DoFHandler<dim,spacedim>::active_cell_iterator cell = dof_handler.begin_active(),
1212  endc = dof_handler.end();
1213 
1214  if (n_components==1)
1215  {
1216  std::vector<double> rhs_values(n_q_points);
1217 
1218  for (; cell!=endc; ++cell)
1219  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
1220  if (cell->face(face)->at_boundary () &&
1221  (boundary_indicators.empty() ||
1222  (boundary_indicators.find (cell->face(face)->boundary_indicator())
1223  !=
1224  boundary_indicators.end())))
1225  {
1226  fe_values.reinit(cell, face);
1227 
1228  const std::vector<double> &weights = fe_values.get_JxW_values ();
1229  rhs_function.value_list (fe_values.get_quadrature_points(), rhs_values);
1230 
1231  cell_vector = 0;
1232  for (unsigned int point=0; point<n_q_points; ++point)
1233  for (unsigned int i=0; i<dofs_per_cell; ++i)
1234  cell_vector(i) += rhs_values[point] *
1235  fe_values.shape_value(i,point) *
1236  weights[point];
1237 
1238  cell->get_dof_indices (dofs);
1239 
1240  for (unsigned int i=0; i<dofs_per_cell; ++i)
1241  rhs_vector(dofs[i]) += cell_vector(i);
1242  }
1243  }
1244  else
1245  {
1246  std::vector<Vector<double> > rhs_values(n_q_points, Vector<double>(n_components));
1247 
1248  for (; cell!=endc; ++cell)
1249  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
1250  if (cell->face(face)->at_boundary () &&
1251  (boundary_indicators.empty() ||
1252  (boundary_indicators.find (cell->face(face)->boundary_indicator())
1253  !=
1254  boundary_indicators.end())))
1255  {
1256  fe_values.reinit(cell, face);
1257 
1258  const std::vector<double> &weights = fe_values.get_JxW_values ();
1259  rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values);
1260 
1261  cell_vector = 0;
1262 
1263  // Use the faster code if the
1264  // FiniteElement is primitive
1265  if (fe.is_primitive ())
1266  {
1267  for (unsigned int point=0; point<n_q_points; ++point)
1268  for (unsigned int i=0; i<dofs_per_cell; ++i)
1269  {
1270  const unsigned int component
1271  = fe.system_to_component_index(i).first;
1272 
1273  cell_vector(i) += rhs_values[point](component) *
1274  fe_values.shape_value(i,point) *
1275  weights[point];
1276  }
1277  }
1278  else
1279  {
1280  // And the full featured
1281  // code, if vector valued
1282  // FEs are used
1283  for (unsigned int point=0; point<n_q_points; ++point)
1284  for (unsigned int i=0; i<dofs_per_cell; ++i)
1285  for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
1286  if (fe.get_nonzero_components(i)[comp_i])
1287  {
1288  cell_vector(i)
1289  += rhs_values[point](comp_i) *
1290  fe_values.shape_value_component(i,point,comp_i) *
1291  weights[point];
1292  }
1293  }
1294 
1295  cell->get_dof_indices (dofs);
1296 
1297  for (unsigned int i=0; i<dofs_per_cell; ++i)
1298  rhs_vector(dofs[i]) += cell_vector(i);
1299  }
1300  }
1301  }
1302 
1303 
1304 
1305  template <int dim, int spacedim>
1306  void
1308  const Quadrature<dim-1> &quadrature,
1309  const Function<spacedim> &rhs_function,
1310  Vector<double> &rhs_vector,
1311  const std::set<types::boundary_id> &boundary_indicators)
1312  {
1314  quadrature,
1315  rhs_function, rhs_vector,
1316  boundary_indicators);
1317  }
1318 
1319 
1320 
1321  template <int dim, int spacedim>
1322  void
1324  const hp::DoFHandler<dim,spacedim> &dof_handler,
1325  const hp::QCollection<dim-1> &quadrature,
1326  const Function<spacedim> &rhs_function,
1327  Vector<double> &rhs_vector,
1328  const std::set<types::boundary_id> &boundary_indicators)
1329  {
1330  const hp::FECollection<dim> &fe = dof_handler.get_fe();
1331  Assert (fe.n_components() == rhs_function.n_components,
1332  ExcDimensionMismatch(fe.n_components(), rhs_function.n_components));
1333  Assert (rhs_vector.size() == dof_handler.n_dofs(),
1334  ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs()));
1335 
1336  rhs_vector = 0;
1337 
1338  UpdateFlags update_flags = UpdateFlags(update_values |
1341  hp::FEFaceValues<dim> x_fe_values (mapping, fe, quadrature, update_flags);
1342 
1343  const unsigned int n_components = fe.n_components();
1344 
1345  std::vector<types::global_dof_index> dofs (fe.max_dofs_per_cell());
1346  Vector<double> cell_vector (fe.max_dofs_per_cell());
1347 
1348  typename hp::DoFHandler<dim,spacedim>::active_cell_iterator
1349  cell = dof_handler.begin_active(),
1350  endc = dof_handler.end();
1351 
1352  if (n_components==1)
1353  {
1354  std::vector<double> rhs_values;
1355 
1356  for (; cell!=endc; ++cell)
1357  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
1358  if (cell->face(face)->at_boundary () &&
1359  (boundary_indicators.empty() ||
1360  (boundary_indicators.find (cell->face(face)->boundary_indicator())
1361  !=
1362  boundary_indicators.end())))
1363  {
1364  x_fe_values.reinit(cell, face);
1365 
1366  const FEFaceValues<dim> &fe_values = x_fe_values.get_present_fe_values();
1367 
1368  const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
1369  n_q_points = fe_values.n_quadrature_points;
1370  rhs_values.resize (n_q_points);
1371 
1372  const std::vector<double> &weights = fe_values.get_JxW_values ();
1373  rhs_function.value_list (fe_values.get_quadrature_points(), rhs_values);
1374 
1375  cell_vector = 0;
1376  for (unsigned int point=0; point<n_q_points; ++point)
1377  for (unsigned int i=0; i<dofs_per_cell; ++i)
1378  cell_vector(i) += rhs_values[point] *
1379  fe_values.shape_value(i,point) *
1380  weights[point];
1381 
1382  dofs.resize(dofs_per_cell);
1383  cell->get_dof_indices (dofs);
1384 
1385  for (unsigned int i=0; i<dofs_per_cell; ++i)
1386  rhs_vector(dofs[i]) += cell_vector(i);
1387  }
1388  }
1389  else
1390  {
1391  std::vector<Vector<double> > rhs_values;
1392 
1393  for (; cell!=endc; ++cell)
1394  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
1395  if (cell->face(face)->at_boundary () &&
1396  (boundary_indicators.empty() ||
1397  (boundary_indicators.find (cell->face(face)->boundary_indicator())
1398  !=
1399  boundary_indicators.end())))
1400  {
1401  x_fe_values.reinit(cell, face);
1402 
1403  const FEFaceValues<dim> &fe_values = x_fe_values.get_present_fe_values();
1404 
1405  const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
1406  n_q_points = fe_values.n_quadrature_points;
1407  rhs_values.resize (n_q_points, Vector<double>(n_components));
1408 
1409  const std::vector<double> &weights = fe_values.get_JxW_values ();
1410  rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values);
1411 
1412  cell_vector = 0;
1413 
1414  // Use the faster code if the
1415  // FiniteElement is primitive
1416  if (cell->get_fe().is_primitive ())
1417  {
1418  for (unsigned int point=0; point<n_q_points; ++point)
1419  for (unsigned int i=0; i<dofs_per_cell; ++i)
1420  {
1421  const unsigned int component
1422  = cell->get_fe().system_to_component_index(i).first;
1423 
1424  cell_vector(i) += rhs_values[point](component) *
1425  fe_values.shape_value(i,point) *
1426  weights[point];
1427  }
1428  }
1429  else
1430  {
1431  // And the full featured
1432  // code, if vector valued
1433  // FEs are used
1434  for (unsigned int point=0; point<n_q_points; ++point)
1435  for (unsigned int i=0; i<dofs_per_cell; ++i)
1436  for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
1437  if (cell->get_fe().get_nonzero_components(i)[comp_i])
1438  {
1439  cell_vector(i)
1440  += rhs_values[point](comp_i) *
1441  fe_values.shape_value_component(i,point,comp_i) *
1442  weights[point];
1443  }
1444  }
1445  dofs.resize(dofs_per_cell);
1446  cell->get_dof_indices (dofs);
1447 
1448  for (unsigned int i=0; i<dofs_per_cell; ++i)
1449  rhs_vector(dofs[i]) += cell_vector(i);
1450  }
1451  }
1452  }
1453 
1454 
1455 
1456  template <int dim, int spacedim>
1457  void
1459  const hp::QCollection<dim-1> &quadrature,
1460  const Function<spacedim> &rhs_function,
1461  Vector<double> &rhs_vector,
1462  const std::set<types::boundary_id> &boundary_indicators)
1463  {
1465  dof_handler, quadrature,
1466  rhs_function, rhs_vector,
1467  boundary_indicators);
1468  }
1469 
1470 
1471 
1472 // ----------- interpolate_boundary_values for std::map --------------------
1473 
1474  namespace
1475  {
1476  // interpolate boundary values in
1477  // 1d. in higher dimensions, we
1478  // use FEValues to figure out
1479  // what to do on faces, but in 1d
1480  // faces are points and it is far
1481  // easier to simply work on
1482  // individual vertices
1483  template <class DH,
1484  template <int,int> class M_or_MC>
1485  static inline
1486  void do_interpolate_boundary_values (const M_or_MC<DH::dimension, DH::space_dimension> &,
1487  const DH &dof,
1488  const typename FunctionMap<DH::space_dimension>::type &function_map,
1489  std::map<types::global_dof_index,double> &boundary_values,
1490  const ComponentMask &component_mask,
1491  const ::internal::int2type<1>)
1492  {
1493  const unsigned int dim = DH::dimension;
1494  const unsigned int spacedim=DH::space_dimension;
1495 
1496  Assert (component_mask.represents_n_components(dof.get_fe().n_components()),
1497  ExcMessage ("The number of components in the mask has to be either "
1498  "zero or equal to the number of components in the finite "
1499  "element."));
1500 
1501  // if for whatever reason we were
1502  // passed an empty map, return
1503  // immediately
1504  if (function_map.size() == 0)
1505  return;
1506 
1507  for (typename DH::active_cell_iterator cell = dof.begin_active();
1508  cell != dof.end(); ++cell)
1509  for (unsigned int direction=0;
1510  direction<GeometryInfo<dim>::faces_per_cell; ++direction)
1511  if (cell->at_boundary(direction)
1512  &&
1513  (function_map.find(cell->face(direction)->boundary_indicator()) != function_map.end()))
1514  {
1515  const Function<DH::space_dimension> &boundary_function
1516  = *function_map.find(cell->face(direction)->boundary_indicator())->second;
1517 
1518  // get the FE corresponding to this
1519  // cell
1520  const FiniteElement<dim,spacedim> &fe = cell->get_fe();
1521  Assert (fe.n_components() == boundary_function.n_components,
1523  boundary_function.n_components));
1524 
1525  Assert (component_mask.n_selected_components(fe.n_components()) > 0,
1526  ComponentMask::ExcNoComponentSelected());
1527 
1528  // now set the value of
1529  // the vertex degree of
1530  // freedom. setting
1531  // also creates the
1532  // entry in the map if
1533  // it did not exist
1534  // beforehand
1535  //
1536  // save some time by
1537  // requesting values
1538  // only once for each
1539  // point, irrespective
1540  // of the number of
1541  // components of the
1542  // function
1543  ::Vector<double> function_values (fe.n_components());
1544  if (fe.n_components() == 1)
1545  function_values(0)
1546  = boundary_function.value (cell->vertex(direction));
1547  else
1548  boundary_function.vector_value (cell->vertex(direction),
1549  function_values);
1550 
1551  for (unsigned int i=0; i<fe.dofs_per_vertex; ++i)
1552  if (component_mask[fe.face_system_to_component_index(i).first])
1553  boundary_values[cell->
1554  vertex_dof_index(direction,i,
1555  cell->active_fe_index())]
1556  = function_values(fe.face_system_to_component_index(i).first);
1557  }
1558  }
1559 
1560 
1561 
1562  template <class DH,
1563  template <int,int> class M_or_MC>
1564  static inline
1565  void
1566  do_interpolate_boundary_values (const M_or_MC<DH::dimension, DH::space_dimension> &mapping,
1567  const DH &dof,
1568  const typename FunctionMap<DH::space_dimension>::type &function_map,
1569  std::map<types::global_dof_index,double> &boundary_values,
1570  const ComponentMask &component_mask,
1571  const ::internal::int2type<DH::dimension>)
1572  {
1573  const unsigned int dim = DH::dimension;
1574  const unsigned int spacedim=DH::space_dimension;
1575 
1576  Assert (component_mask.represents_n_components(dof.get_fe().n_components()),
1577  ExcMessage ("The number of components in the mask has to be either "
1578  "zero or equal to the number of components in the finite "
1579  "element."));
1580 
1581 
1582  // if for whatever reason we were passed an empty map, return
1583  // immediately
1584  if (function_map.size() == 0)
1585  return;
1586 
1587  Assert (function_map.find(numbers::internal_face_boundary_id) == function_map.end(),
1588  ExcInvalidBoundaryIndicator());
1589 
1590  const unsigned int n_components = DoFTools::n_components(dof);
1591  const bool fe_is_system = (n_components != 1);
1592 
1593  for (typename FunctionMap<spacedim>::type::const_iterator i=function_map.begin();
1594  i!=function_map.end(); ++i)
1595  Assert (n_components == i->second->n_components,
1596  ExcDimensionMismatch(n_components, i->second->n_components));
1597 
1598  // field to store the indices
1599  std::vector<types::global_dof_index> face_dofs;
1600  face_dofs.reserve (DoFTools::max_dofs_per_face(dof));
1601 
1602  std::vector<Point<spacedim> > dof_locations;
1603  dof_locations.reserve (DoFTools::max_dofs_per_face(dof));
1604 
1605  // array to store the values of the boundary function at the boundary
1606  // points. have two arrays for scalar and vector functions to use the
1607  // more efficient one respectively
1608  std::vector<double> dof_values_scalar;
1609  std::vector<::Vector<double> > dof_values_system;
1610  dof_values_scalar.reserve (DoFTools::max_dofs_per_face (dof));
1611  dof_values_system.reserve (DoFTools::max_dofs_per_face (dof));
1612 
1613  // before we start with the loop over all cells create an hp::FEValues
1614  // object that holds the interpolation points of all finite elements
1615  // that may ever be in use
1616  ::hp::FECollection<dim,spacedim> finite_elements (dof.get_fe());
1617  ::hp::QCollection<dim-1> q_collection;
1618  for (unsigned int f=0; f<finite_elements.size(); ++f)
1619  {
1620  const FiniteElement<dim,spacedim> &fe = finite_elements[f];
1621 
1622  // generate a quadrature rule on the face from the unit support
1623  // points. this will be used to obtain the quadrature points on the
1624  // real cell's face
1625  //
1626  // to do this, we check whether the FE has support points on the
1627  // face at all:
1628  if (fe.has_face_support_points())
1630  else
1631  {
1632  // if not, then we should try a more clever way. the idea is
1633  // that a finite element may not offer support points for all
1634  // its shape functions, but maybe only some. if it offers
1635  // support points for the components we are interested in in
1636  // this function, then that's fine. if not, the function we call
1637  // in the finite element will raise an exception. the support
1638  // points for the other shape functions are left uninitialized
1639  // (well, initialized by the default constructor), since we
1640  // don't need them anyway.
1641  //
1642  // As a detour, we must make sure we only query
1643  // face_system_to_component_index if the index corresponds to a
1644  // primitive shape function. since we know that all the
1645  // components we are interested in are primitive (by the above
1646  // check), we can safely put such a check in front
1647  std::vector<Point<dim-1> > unit_support_points (fe.dofs_per_face);
1648 
1649  for (unsigned int i=0; i<fe.dofs_per_face; ++i)
1650  if (fe.is_primitive (fe.face_to_cell_index(i,0)))
1651  if (component_mask[fe.face_system_to_component_index(i).first]
1652  == true)
1653  unit_support_points[i] = fe.unit_face_support_point(i);
1654 
1655  q_collection.push_back (Quadrature<dim-1>(unit_support_points));
1656  }
1657  }
1658  // now that we have a q_collection object with all the right quadrature
1659  // points, create an hp::FEFaceValues object that we can use to evaluate
1660  // the boundary values at
1661  ::hp::MappingCollection<dim,spacedim> mapping_collection (mapping);
1662  ::hp::FEFaceValues<dim,spacedim> x_fe_values (mapping_collection, finite_elements, q_collection,
1664 
1665  typename DH::active_cell_iterator cell = dof.begin_active(),
1666  endc = dof.end();
1667  for (; cell!=endc; ++cell)
1668  if (!cell->is_artificial())
1669  for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
1670  ++face_no)
1671  {
1672  const FiniteElement<dim,spacedim> &fe = cell->get_fe();
1673 
1674  // we can presently deal only with primitive elements for
1675  // boundary values. this does not preclude us using
1676  // non-primitive elements in components that we aren't
1677  // interested in, however. make sure that all shape functions
1678  // that are non-zero for the components we are interested in,
1679  // are in fact primitive
1680  for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
1681  {
1682  const ComponentMask &nonzero_component_array
1683  = cell->get_fe().get_nonzero_components (i);
1684  for (unsigned int c=0; c<n_components; ++c)
1685  if ((nonzero_component_array[c] == true)
1686  &&
1687  (component_mask[c] == true))
1688  Assert (cell->get_fe().is_primitive (i),
1689  ExcMessage ("This function can only deal with requested boundary "
1690  "values that correspond to primitive (scalar) base "
1691  "elements"));
1692  }
1693 
1694  const typename DH::face_iterator face = cell->face(face_no);
1695  const types::boundary_id boundary_component = face->boundary_indicator();
1696 
1697  // see if this face is part of the boundaries for which we are
1698  // supposed to do something, and also see if the finite element
1699  // in use here has DoFs on the face at all
1700  if ((function_map.find(boundary_component) != function_map.end())
1701  &&
1702  (cell->get_fe().dofs_per_face > 0))
1703  {
1704  // face is of the right component
1705  x_fe_values.reinit(cell, face_no);
1706  const ::FEFaceValues<dim,spacedim> &fe_values =
1707  x_fe_values.get_present_fe_values();
1708 
1709  // get indices, physical location and boundary values of
1710  // dofs on this face
1711  face_dofs.resize (fe.dofs_per_face);
1712  face->get_dof_indices (face_dofs, cell->active_fe_index());
1713  const std::vector<Point<spacedim> > &dof_locations
1714  = fe_values.get_quadrature_points ();
1715 
1716  if (fe_is_system)
1717  {
1718  // resize array. avoid construction of a memory
1719  // allocating temporary if possible
1720  if (dof_values_system.size() < fe.dofs_per_face)
1721  dof_values_system.resize (fe.dofs_per_face,
1722  ::Vector<double>(fe.n_components()));
1723  else
1724  dof_values_system.resize (fe.dofs_per_face);
1725 
1726  function_map.find(boundary_component)->second
1727  ->vector_value_list (dof_locations, dof_values_system);
1728 
1729  // enter those dofs into the list that match the
1730  // component signature. avoid the usual complication
1731  // that we can't just use *_system_to_component_index
1732  // for non-primitive FEs
1733  for (unsigned int i=0; i<face_dofs.size(); ++i)
1734  {
1735  unsigned int component;
1736  if (fe.is_primitive())
1737  component = fe.face_system_to_component_index(i).first;
1738  else
1739  {
1740  // non-primitive case. make sure that this
1741  // particular shape function _is_ primitive, and
1742  // get at it's component. use usual trick to
1743  // transfer face dof index to cell dof index
1744  const unsigned int cell_i
1745  = (dim == 1 ?
1746  i
1747  :
1748  (dim == 2 ?
1749  (i<2*fe.dofs_per_vertex ? i : i+2*fe.dofs_per_vertex)
1750  :
1751  (dim == 3 ?
1752  (i<4*fe.dofs_per_vertex ?
1753  i
1754  :
1755  (i<4*fe.dofs_per_vertex+4*fe.dofs_per_line ?
1756  i+4*fe.dofs_per_vertex
1757  :
1758  i+4*fe.dofs_per_vertex+8*fe.dofs_per_line))
1759  :
1761  Assert (cell_i < fe.dofs_per_cell, ExcInternalError());
1762 
1763  // make sure that if this is not a primitive
1764  // shape function, then all the corresponding
1765  // components in the mask are not set
1766  if (!fe.is_primitive(cell_i))
1767  for (unsigned int c=0; c<n_components; ++c)
1768  if (fe.get_nonzero_components(cell_i)[c])
1769  Assert (component_mask[c] == false,
1770  FETools::ExcFENotPrimitive());
1771 
1772  // let's pick the first of possibly more than
1773  // one non-zero components. if shape function is
1774  // non-primitive, then we will ignore the result
1775  // in the following anyway, otherwise there's
1776  // only one non-zero component which we will use
1777  component = fe.get_nonzero_components(cell_i).first_selected_component();
1778  }
1779 
1780  if (component_mask[component] == true)
1781  boundary_values[face_dofs[i]] = dof_values_system[i](component);
1782  }
1783  }
1784  else
1785  // fe has only one component, so save some computations
1786  {
1787  // get only the one component that this function has
1788  dof_values_scalar.resize (fe.dofs_per_face);
1789  function_map.find(boundary_component)->second
1790  ->value_list (dof_locations, dof_values_scalar, 0);
1791 
1792  // enter into list
1793 
1794  for (unsigned int i=0; i<face_dofs.size(); ++i)
1795  boundary_values[face_dofs[i]] = dof_values_scalar[i];
1796  }
1797  }
1798  }
1799  } // end of interpolate_boundary_values
1800  } // end of namespace internal
1801 
1802 
1803 
1804  template <class DH>
1805  void
1806 
1808  const DH &dof,
1809  const typename FunctionMap<DH::space_dimension>::type &function_map,
1810  std::map<types::global_dof_index,double> &boundary_values,
1811  const ComponentMask &component_mask_)
1812  {
1813  do_interpolate_boundary_values (mapping, dof, function_map, boundary_values,
1814  component_mask_,
1816  }
1817 
1818 
1819 
1820  template <class DH>
1821  void
1823  const DH &dof,
1824  const types::boundary_id boundary_component,
1825  const Function<DH::space_dimension> &boundary_function,
1826  std::map<types::global_dof_index,double> &boundary_values,
1827  const ComponentMask &component_mask)
1828  {
1829  typename FunctionMap<DH::space_dimension>::type function_map;
1830  function_map[boundary_component] = &boundary_function;
1831  interpolate_boundary_values (mapping, dof, function_map, boundary_values,
1832  component_mask);
1833  }
1834 
1835 
1836  template <int dim, int spacedim>
1837  void
1839  const hp::DoFHandler<dim,spacedim> &dof,
1840  const typename FunctionMap<spacedim>::type &function_map,
1841  std::map<types::global_dof_index,double> &boundary_values,
1842  const ComponentMask &component_mask_)
1843  {
1844  do_interpolate_boundary_values (mapping, dof, function_map, boundary_values,
1845  component_mask_,
1847  }
1848 
1849 
1850 
1851  template <class DH>
1852  void
1854  const types::boundary_id boundary_component,
1855  const Function<DH::space_dimension> &boundary_function,
1856  std::map<types::global_dof_index,double> &boundary_values,
1857  const ComponentMask &component_mask)
1858  {
1860  dof, boundary_component,
1861  boundary_function, boundary_values, component_mask);
1862  }
1863 
1864 
1865 
1866  template <class DH>
1867  void
1869  const typename FunctionMap<DH::space_dimension>::type &function_map,
1870  std::map<types::global_dof_index,double> &boundary_values,
1871  const ComponentMask &component_mask)
1872  {
1874  dof, function_map,
1875  boundary_values, component_mask);
1876  }
1877 
1878 
1879 
1880 
1881 // ----------- interpolate_boundary_values for ConstraintMatrix --------------
1882 
1883 
1884 
1885  template <class DH>
1886  void
1889  const DH &dof,
1890  const typename FunctionMap<DH::space_dimension>::type &function_map,
1891  ConstraintMatrix &constraints,
1892  const ComponentMask &component_mask_)
1893  {
1894  std::map<types::global_dof_index,double> boundary_values;
1895  interpolate_boundary_values (mapping, dof, function_map,
1896  boundary_values, component_mask_);
1897  std::map<types::global_dof_index,double>::const_iterator boundary_value =
1898  boundary_values.begin();
1899  for ( ; boundary_value !=boundary_values.end(); ++boundary_value)
1900  {
1901  if (constraints.can_store_line (boundary_value->first)
1902  &&
1903  !constraints.is_constrained(boundary_value->first))
1904  {
1905  constraints.add_line (boundary_value->first);
1906  constraints.set_inhomogeneity (boundary_value->first,
1907  boundary_value->second);
1908  }
1909  }
1910  }
1911 
1912 
1913 
1914  template <class DH>
1915  void
1918  const DH &dof,
1919  const types::boundary_id boundary_component,
1920  const Function<DH::space_dimension> &boundary_function,
1921  ConstraintMatrix &constraints,
1922  const ComponentMask &component_mask)
1923  {
1924  typename FunctionMap<DH::space_dimension>::type function_map;
1925  function_map[boundary_component] = &boundary_function;
1926  interpolate_boundary_values (mapping, dof, function_map, constraints,
1927  component_mask);
1928  }
1929 
1930 
1931 
1932  template <class DH>
1933  void
1935  (const DH &dof,
1936  const types::boundary_id boundary_component,
1937  const Function<DH::space_dimension> &boundary_function,
1938  ConstraintMatrix &constraints,
1939  const ComponentMask &component_mask)
1940  {
1942  dof, boundary_component,
1943  boundary_function, constraints, component_mask);
1944  }
1945 
1946 
1947 
1948  template <class DH>
1949  void
1951  (const DH &dof,
1952  const typename FunctionMap<DH::space_dimension>::type &function_map,
1953  ConstraintMatrix &constraints,
1954  const ComponentMask &component_mask)
1955  {
1957  dof, function_map,
1958  constraints, component_mask);
1959  }
1960 
1961 
1962 
1963 
1964 // -------- implementation for project_boundary_values with std::map --------
1965 
1966 
1967  namespace
1968  {
1969  template <int dim, int spacedim,
1970  template <int,int> class DH,
1971  template <int,int> class M_or_MC,
1972  template <int> class Q_or_QC>
1973  void
1974  do_project_boundary_values (const M_or_MC<dim, spacedim> &mapping,
1975  const DH<dim, spacedim> &dof,
1976  const typename FunctionMap<spacedim>::type &boundary_functions,
1977  const Q_or_QC<dim-1> &q,
1978  std::map<types::global_dof_index,double> &boundary_values,
1979  std::vector<unsigned int> component_mapping)
1980  {
1981  // in 1d, projection onto the 0d end points == interpolation
1982  if (dim == 1)
1983  {
1984  Assert (component_mapping.size() == 0, ExcNotImplemented());
1985  interpolate_boundary_values (mapping, dof, boundary_functions,
1986  boundary_values, ComponentMask());
1987  return;
1988  }
1989 
1990  //TODO:[?] In project_boundary_values, no condensation of sparsity
1991  // structures, matrices and right hand sides or distribution of
1992  // solution vectors is performed. This is ok for dim<3 because then
1993  // there are no constrained nodes on the boundary, but is not
1994  // acceptable for higher dimensions. Fix this.
1995 
1996  if (component_mapping.size() == 0)
1997  {
1998  AssertDimension (dof.get_fe().n_components(), boundary_functions.begin()->second->n_components);
1999  // I still do not see why i
2000  // should create another copy
2001  // here
2002  component_mapping.resize(dof.get_fe().n_components());
2003  for (unsigned int i= 0 ; i < component_mapping.size() ; ++i)
2004  component_mapping[i] = i;
2005  }
2006  else
2007  AssertDimension (dof.get_fe().n_components(), component_mapping.size());
2008 
2009  std::vector<types::global_dof_index> dof_to_boundary_mapping;
2010  std::set<types::boundary_id> selected_boundary_components;
2011  for (typename FunctionMap<spacedim>::type::const_iterator i=boundary_functions.begin();
2012  i!=boundary_functions.end(); ++i)
2013  selected_boundary_components.insert (i->first);
2014 
2015  DoFTools::map_dof_to_boundary_indices (dof, selected_boundary_components,
2016  dof_to_boundary_mapping);
2017 
2018  // Done if no degrees of freedom on the boundary
2019  if (dof.n_boundary_dofs (boundary_functions) == 0)
2020  return;
2021 
2022  // set up sparsity structure
2023  CompressedSparsityPattern c_sparsity(dof.n_boundary_dofs (boundary_functions),
2024  dof.n_boundary_dofs (boundary_functions));
2026  boundary_functions,
2027  dof_to_boundary_mapping,
2028  c_sparsity);
2029  SparsityPattern sparsity;
2030  sparsity.copy_from(c_sparsity);
2031 
2032 
2033 
2034  // note: for three or more dimensions, there
2035  // may be constrained nodes on the boundary
2036  // in this case the boundary mass matrix has
2037  // to be condensed and the solution is to
2038  // be distributed afterwards, which is not
2039  // yet implemented. The reason for this is
2040  // that we cannot simply use the condense
2041  // family of functions, since the matrices
2042  // and vectors do not use the global
2043  // numbering but rather the boundary
2044  // numbering, i.e. the condense function
2045  // needs to use another indirection. There
2046  // should be not many technical problems,
2047  // but it needs to be implemented
2048  if (dim>=3)
2049  {
2050 #ifdef DEBUG
2051  // Assert that there are no hanging nodes at the boundary
2052  int level = -1;
2053  for (typename DH<dim,spacedim>::active_cell_iterator cell = dof.begin_active();
2054  cell != dof.end(); ++cell)
2055  for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
2056  {
2057  if (cell->at_boundary(f))
2058  {
2059  if (level == -1)
2060  level = cell->level();
2061  else
2062  {
2063  Assert (level == cell->level(), ExcNotImplemented());
2064  }
2065  }
2066  }
2067 #endif
2068  }
2069  sparsity.compress();
2070 
2071 
2072  // make mass matrix and right hand side
2073  SparseMatrix<double> mass_matrix(sparsity);
2074  Vector<double> rhs(sparsity.n_rows());
2075 
2076 
2078  mass_matrix, boundary_functions,
2079  rhs, dof_to_boundary_mapping, (const Function<spacedim> *) 0,
2080  component_mapping);
2081 
2082  // For certain weird elements,
2083  // there might be degrees of
2084  // freedom on the boundary, but
2085  // their shape functions do not
2086  // have support there. Let's
2087  // eliminate them here.
2088 
2089  // The Bogner-Fox-Schmidt element
2090  // is an example for those.
2091 
2092 //TODO: Maybe we should figure out if the element really needs this
2093 
2094  FilteredMatrix<Vector<double> > filtered_mass_matrix(mass_matrix, true);
2095  FilteredMatrix<Vector<double> > filtered_precondition;
2096  std::vector<bool> excluded_dofs(mass_matrix.m(), false);
2097 
2098  double max_element = 0.;
2099  for (unsigned int i=0; i<mass_matrix.m(); ++i)
2100  if (mass_matrix.diag_element(i) > max_element)
2101  max_element = mass_matrix.diag_element(i);
2102 
2103  for (unsigned int i=0; i<mass_matrix.m(); ++i)
2104  if (mass_matrix.diag_element(i) < 1.e-8 * max_element)
2105  {
2106  filtered_mass_matrix.add_constraint(i, 0.);
2107  filtered_precondition.add_constraint(i, 0.);
2108  mass_matrix.diag_element(i) = 1.;
2109  excluded_dofs[i] = true;
2110  }
2111 
2112  Vector<double> boundary_projection (rhs.size());
2113 
2114  // Allow for a maximum of 5*n
2115  // steps to reduce the residual by
2116  // 10^-12. n steps may not be
2117  // sufficient, since roundoff
2118  // errors may accumulate for badly
2119  // conditioned matrices
2120  ReductionControl control(5*rhs.size(), 0., 1.e-12, false, false);
2121  GrowingVectorMemory<> memory;
2122  SolverCG<> cg(control,memory);
2123 
2124  PreconditionSSOR<> prec;
2125  prec.initialize(mass_matrix, 1.2);
2126  filtered_precondition.initialize(prec, true);
2127  // solve
2128  cg.solve (filtered_mass_matrix, boundary_projection, rhs, filtered_precondition);
2129  filtered_precondition.apply_constraints(boundary_projection, true);
2130  filtered_precondition.clear();
2131  // fill in boundary values
2132  for (unsigned int i=0; i<dof_to_boundary_mapping.size(); ++i)
2133  if (dof_to_boundary_mapping[i] != DoFHandler<dim,spacedim>::invalid_dof_index
2134  && ! excluded_dofs[dof_to_boundary_mapping[i]])
2135  {
2136  Assert(numbers::is_finite(boundary_projection(dof_to_boundary_mapping[i])), ExcNumberNotFinite());
2137 
2138  // this dof is on one of the
2139  // interesting boundary parts
2140  //
2141  // remember: i is the global dof
2142  // number, dof_to_boundary_mapping[i]
2143  // is the number on the boundary and
2144  // thus in the solution vector
2145  boundary_values[i] = boundary_projection(dof_to_boundary_mapping[i]);
2146  }
2147  }
2148  }
2149 
2150  template <int dim, int spacedim>
2151  void
2153  const DoFHandler<dim, spacedim> &dof,
2154  const typename FunctionMap<spacedim>::type &boundary_functions,
2155  const Quadrature<dim-1> &q,
2156  std::map<types::global_dof_index,double> &boundary_values,
2157  std::vector<unsigned int> component_mapping)
2158  {
2159  do_project_boundary_values(mapping, dof, boundary_functions, q,
2160  boundary_values, component_mapping);
2161  }
2162 
2163 
2164 
2165  template <int dim, int spacedim>
2166  void
2168  const typename FunctionMap<spacedim>::type &boundary_functions,
2169  const Quadrature<dim-1> &q,
2170  std::map<types::global_dof_index,double> &boundary_values,
2171  std::vector<unsigned int> component_mapping)
2172  {
2174  boundary_values, component_mapping);
2175  }
2176 
2177 
2178 
2179  template <int dim, int spacedim>
2181  const hp::DoFHandler<dim,spacedim> &dof,
2182  const typename FunctionMap<spacedim>::type &boundary_functions,
2183  const hp::QCollection<dim-1> &q,
2184  std::map<types::global_dof_index,double> &boundary_values,
2185  std::vector<unsigned int> component_mapping)
2186  {
2187  do_project_boundary_values (mapping, dof,
2188  boundary_functions,
2189  q, boundary_values,
2190  component_mapping);
2191  }
2192 
2193 
2194 
2195  template <int dim, int spacedim>
2197  const typename FunctionMap<spacedim>::type &boundary_function,
2198  const hp::QCollection<dim-1> &q,
2199  std::map<types::global_dof_index,double> &boundary_values,
2200  std::vector<unsigned int> component_mapping)
2201  {
2203  boundary_function,
2204  q, boundary_values,
2205  component_mapping);
2206  }
2207 
2208 
2209  // ----- implementation for project_boundary_values with ConstraintMatrix -----
2210 
2211 
2212 
2213  template <int dim, int spacedim>
2214  void
2216  const DoFHandler<dim,spacedim> &dof,
2217  const typename FunctionMap<spacedim>::type &boundary_functions,
2218  const Quadrature<dim-1> &q,
2219  ConstraintMatrix &constraints,
2220  std::vector<unsigned int> component_mapping)
2221  {
2222  std::map<types::global_dof_index,double> boundary_values;
2223  project_boundary_values (mapping, dof, boundary_functions, q,
2224  boundary_values, component_mapping);
2225  std::map<types::global_dof_index,double>::const_iterator boundary_value =
2226  boundary_values.begin();
2227  for ( ; boundary_value !=boundary_values.end(); ++boundary_value)
2228  {
2229  if (!constraints.is_constrained(boundary_value->first))
2230  {
2231  constraints.add_line (boundary_value->first);
2232  constraints.set_inhomogeneity (boundary_value->first,
2233  boundary_value->second);
2234  }
2235  }
2236  }
2237 
2238 
2239 
2240  template <int dim, int spacedim>
2241  void
2243  const typename FunctionMap<spacedim>::type &boundary_functions,
2244  const Quadrature<dim-1> &q,
2245  ConstraintMatrix &constraints,
2246  std::vector<unsigned int> component_mapping)
2247  {
2249  constraints, component_mapping);
2250  }
2251 
2252 
2253 
2254 
2255  namespace internal
2256  {
2263  template <int dim>
2265  {
2266  types::global_dof_index dof_indices[dim];
2267 
2268  VectorDoFTuple ()
2269  {
2270  for (unsigned int i=0; i<dim; ++i)
2271  dof_indices[i] = numbers::invalid_dof_index;
2272  }
2273 
2274 
2275  bool operator < (const VectorDoFTuple<dim> &other) const
2276  {
2277  for (unsigned int i=0; i<dim; ++i)
2278  if (dof_indices[i] < other.dof_indices[i])
2279  return true;
2280  else if (dof_indices[i] > other.dof_indices[i])
2281  return false;
2282  return false;
2283  }
2284 
2285  bool operator == (const VectorDoFTuple<dim> &other) const
2286  {
2287  for (unsigned int i=0; i<dim; ++i)
2288  if (dof_indices[i] != other.dof_indices[i])
2289  return false;
2290 
2291  return true;
2292  }
2293 
2294  bool operator != (const VectorDoFTuple<dim> &other) const
2295  {
2296  return ! (*this == other);
2297  }
2298  };
2299 
2300 
2301  template <int dim>
2302  std::ostream &operator << (std::ostream &out,
2303  const VectorDoFTuple<dim> &vdt)
2304  {
2305  for (unsigned int d=0; d<dim; ++d)
2306  out << vdt.dof_indices[d] << (d < dim-1 ? " " : "");
2307  return out;
2308  }
2309 
2310 
2311 
2327  template <int dim>
2328  void
2329  add_constraint (const VectorDoFTuple<dim> &dof_indices,
2330  const Tensor<1,dim> &constraining_vector,
2331  ConstraintMatrix &constraints)
2332  {
2333 
2334  // choose the DoF that has the
2335  // largest component in the
2336  // constraining_vector as the
2337  // one to be constrained as
2338  // this makes the process
2339  // stable in cases where the
2340  // constraining_vector has the
2341  // form n=(1,0) or n=(0,1)
2342  //
2343  // we get constraints of the form
2344  // x0 = a_1*x1 + a2*x2 + ...
2345  // if one of the weights is
2346  // essentially zero then skip
2347  // this part. the ConstraintMatrix
2348  // can also deal with cases like
2349  // x0 = 0
2350  // if necessary
2351  //
2352  // there is a problem if we have a
2353  // normal vector of the form
2354  // (a,a,small) or (a,a,a). Depending on
2355  // round-off we may choose the first or
2356  // second component (or third, in the
2357  // latter case) as the largest one, and
2358  // depending on our choice one or
2359  // another degree of freedom will be
2360  // constrained. On a single processor
2361  // this is not much of a problem, but
2362  // it's a nightmare when we run in
2363  // parallel and two processors disagree
2364  // on which DoF should be
2365  // constrained. This led to an
2366  // incredibly difficult to find bug in
2367  // @ref step_32 "step-32" when running in parallel
2368  // with 9 or more processors.
2369  //
2370  // in practice, such normal vectors of
2371  // the form (a,a,small) or (a,a,a)
2372  // happen not infrequently since they
2373  // lie on the diagonals where vertices
2374  // frequently happen to land upon mesh
2375  // refinement if one starts from a
2376  // symmetric and regular body. we work
2377  // around this problem in the following
2378  // way: if we have a normal vector of
2379  // the form (a,b) (similarly algorithm
2380  // in 3d), we choose 'a' as the largest
2381  // coefficient not if a>b but if
2382  // a>b+1e-10. this shifts the problem
2383  // away from the frequently visited
2384  // diagonal to a line that is off the
2385  // diagonal. there will of course be
2386  // problems where the exact values of a
2387  // and b differ by exactly 1e-10 and we
2388  // get into the same instability, but
2389  // from a practical viewpoint such
2390  // problems should be much rarer. in
2391  // particular, meshes have to be very
2392  // very fine for a vertex to land on
2393  // this line if the original body had a
2394  // vertex on the diagonal as well
2395  switch (dim)
2396  {
2397  case 2:
2398  {
2399  if (std::fabs(constraining_vector[0]) > std::fabs(constraining_vector[1]) + 1e-10)
2400  {
2401  if (!constraints.is_constrained(dof_indices.dof_indices[0])
2402  &&
2403  constraints.can_store_line(dof_indices.dof_indices[0]))
2404  {
2405  constraints.add_line (dof_indices.dof_indices[0]);
2406 
2407  if (std::fabs (constraining_vector[1]/constraining_vector[0])
2408  > std::numeric_limits<double>::epsilon())
2409  constraints.add_entry (dof_indices.dof_indices[0],
2410  dof_indices.dof_indices[1],
2411  -constraining_vector[1]/constraining_vector[0]);
2412  }
2413  }
2414  else
2415  {
2416  if (!constraints.is_constrained(dof_indices.dof_indices[1])
2417  &&
2418  constraints.can_store_line(dof_indices.dof_indices[1]))
2419  {
2420  constraints.add_line (dof_indices.dof_indices[1]);
2421 
2422  if (std::fabs (constraining_vector[0]/constraining_vector[1])
2423  > std::numeric_limits<double>::epsilon())
2424  constraints.add_entry (dof_indices.dof_indices[1],
2425  dof_indices.dof_indices[0],
2426  -constraining_vector[0]/constraining_vector[1]);
2427  }
2428  }
2429  break;
2430  }
2431 
2432  case 3:
2433  {
2434  if ((std::fabs(constraining_vector[0]) >= std::fabs(constraining_vector[1])+1e-10)
2435  &&
2436  (std::fabs(constraining_vector[0]) >= std::fabs(constraining_vector[2])+2e-10))
2437  {
2438  if (!constraints.is_constrained(dof_indices.dof_indices[0])
2439  &&
2440  constraints.can_store_line(dof_indices.dof_indices[0]))
2441  {
2442  constraints.add_line (dof_indices.dof_indices[0]);
2443 
2444  if (std::fabs (constraining_vector[1]/constraining_vector[0])
2445  > std::numeric_limits<double>::epsilon())
2446  constraints.add_entry (dof_indices.dof_indices[0],
2447  dof_indices.dof_indices[1],
2448  -constraining_vector[1]/constraining_vector[0]);
2449 
2450  if (std::fabs (constraining_vector[2]/constraining_vector[0])
2451  > std::numeric_limits<double>::epsilon())
2452  constraints.add_entry (dof_indices.dof_indices[0],
2453  dof_indices.dof_indices[2],
2454  -constraining_vector[2]/constraining_vector[0]);
2455  }
2456  }
2457  else if ((std::fabs(constraining_vector[1])+1e-10 >= std::fabs(constraining_vector[0]))
2458  &&
2459  (std::fabs(constraining_vector[1]) >= std::fabs(constraining_vector[2])+1e-10))
2460  {
2461  if (!constraints.is_constrained(dof_indices.dof_indices[1])
2462  &&
2463  constraints.can_store_line(dof_indices.dof_indices[1]))
2464  {
2465  constraints.add_line (dof_indices.dof_indices[1]);
2466 
2467  if (std::fabs (constraining_vector[0]/constraining_vector[1])
2468  > std::numeric_limits<double>::epsilon())
2469  constraints.add_entry (dof_indices.dof_indices[1],
2470  dof_indices.dof_indices[0],
2471  -constraining_vector[0]/constraining_vector[1]);
2472 
2473  if (std::fabs (constraining_vector[2]/constraining_vector[1])
2474  > std::numeric_limits<double>::epsilon())
2475  constraints.add_entry (dof_indices.dof_indices[1],
2476  dof_indices.dof_indices[2],
2477  -constraining_vector[2]/constraining_vector[1]);
2478  }
2479  }
2480  else
2481  {
2482  if (!constraints.is_constrained(dof_indices.dof_indices[2])
2483  &&
2484  constraints.can_store_line(dof_indices.dof_indices[2]))
2485  {
2486  constraints.add_line (dof_indices.dof_indices[2]);
2487 
2488  if (std::fabs (constraining_vector[0]/constraining_vector[2])
2489  > std::numeric_limits<double>::epsilon())
2490  constraints.add_entry (dof_indices.dof_indices[2],
2491  dof_indices.dof_indices[0],
2492  -constraining_vector[0]/constraining_vector[2]);
2493 
2494  if (std::fabs (constraining_vector[1]/constraining_vector[2])
2495  > std::numeric_limits<double>::epsilon())
2496  constraints.add_entry (dof_indices.dof_indices[2],
2497  dof_indices.dof_indices[1],
2498  -constraining_vector[1]/constraining_vector[2]);
2499  }
2500  }
2501 
2502  break;
2503  }
2504 
2505  default:
2506  Assert (false, ExcNotImplemented());
2507  }
2508  }
2509 
2510 
2528  template <int dim>
2529  void
2530  add_tangentiality_constraints (const VectorDoFTuple<dim> &dof_indices,
2531  const Tensor<1,dim> &tangent_vector,
2532  ConstraintMatrix &constraints)
2533  {
2534 
2535  // choose the DoF that has the
2536  // largest component in the
2537  // tangent_vector as the
2538  // independent component, and
2539  // then constrain the others to
2540  // it. specifically, if, say,
2541  // component 0 of the tangent
2542  // vector t is largest by
2543  // magnitude, then
2544  // x1=t[1]/t[0]*x_0, etc.
2545  unsigned int largest_component = 0;
2546  for (unsigned int d=1; d<dim; ++d)
2547  if (std::fabs(tangent_vector[d]) > std::fabs(tangent_vector[largest_component]) + 1e-10)
2548  largest_component = d;
2549 
2550  // then constrain all of the
2551  // other degrees of freedom in
2552  // terms of the one just found
2553  for (unsigned int d=0; d<dim; ++d)
2554  if (d != largest_component)
2555  if (!constraints.is_constrained(dof_indices.dof_indices[d])
2556  &&
2557  constraints.can_store_line(dof_indices.dof_indices[d]))
2558  {
2559  constraints.add_line (dof_indices.dof_indices[d]);
2560 
2561  if (std::fabs (tangent_vector[d]/tangent_vector[largest_component])
2562  > std::numeric_limits<double>::epsilon())
2563  constraints.add_entry (dof_indices.dof_indices[d],
2564  dof_indices.dof_indices[largest_component],
2565  tangent_vector[d]/tangent_vector[largest_component]);
2566  }
2567  }
2568 
2569 
2570 
2578  template <int dim>
2579  void
2580  compute_orthonormal_vectors (const Tensor<1,dim> &vector,
2581  Tensor<1,dim> (&orthonormals)[dim-1])
2582  {
2583  switch (dim)
2584  {
2585  case 3:
2586  {
2587  // to do this in 3d, take
2588  // one vector that is
2589  // guaranteed to be not
2590  // aligned with the
2591  // average tangent and
2592  // form the cross
2593  // product. this yields
2594  // one vector that is
2595  // certainly
2596  // perpendicular to the
2597  // tangent; then take the
2598  // cross product between
2599  // this vector and the
2600  // tangent and get one
2601  // vector that is
2602  // perpendicular to both
2603 
2604  // construct a
2605  // temporary vector
2606  // by swapping the
2607  // larger two
2608  // components and
2609  // flipping one
2610  // sign; this can
2611  // not be collinear
2612  // with the average
2613  // tangent
2614  Tensor<1,dim> tmp = vector;
2615  if ((std::fabs(tmp[0]) > std::fabs(tmp[1]))
2616  &&
2617  (std::fabs(tmp[0]) > std::fabs(tmp[2])))
2618  {
2619  // entry zero
2620  // is the
2621  // largest
2622  if ((std::fabs(tmp[1]) > std::fabs(tmp[2])))
2623  std::swap (tmp[0], tmp[1]);
2624  else
2625  std::swap (tmp[0], tmp[2]);
2626 
2627  tmp[0] *= -1;
2628  }
2629  else if ((std::fabs(tmp[1]) > std::fabs(tmp[0]))
2630  &&
2631  (std::fabs(tmp[1]) > std::fabs(tmp[2])))
2632  {
2633  // entry one
2634  // is the
2635  // largest
2636  if ((std::fabs(tmp[0]) > std::fabs(tmp[2])))
2637  std::swap (tmp[1], tmp[0]);
2638  else
2639  std::swap (tmp[1], tmp[2]);
2640 
2641  tmp[1] *= -1;
2642  }
2643  else
2644  {
2645  // entry two
2646  // is the
2647  // largest
2648  if ((std::fabs(tmp[0]) > std::fabs(tmp[1])))
2649  std::swap (tmp[2], tmp[0]);
2650  else
2651  std::swap (tmp[2], tmp[1]);
2652 
2653  tmp[2] *= -1;
2654  }
2655 
2656  // make sure the two vectors
2657  // are indeed not collinear
2658  Assert (std::fabs(vector * tmp / vector.norm() / tmp.norm())
2659  <
2660  (1-1e-12),
2661  ExcInternalError());
2662 
2663  // now compute the
2664  // two normals
2665  cross_product (orthonormals[0], vector, tmp);
2666  cross_product (orthonormals[1], vector, orthonormals[0]);
2667 
2668  break;
2669  }
2670 
2671  default:
2672  Assert (false, ExcNotImplemented());
2673  }
2674  }
2675  }
2676 
2677 
2678  namespace internals
2679  {
2680  // This function computes the
2681  // projection of the boundary
2682  // function on edges for 3D.
2683  template<typename cell_iterator>
2684  void
2685  compute_edge_projection (const cell_iterator &cell,
2686  const unsigned int face,
2687  const unsigned int line,
2688  hp::FEValues<3> &hp_fe_values,
2689  const Function<3> &boundary_function,
2690  const unsigned int first_vector_component,
2691  std::vector<double> &dof_values,
2692  std::vector<bool> &dofs_processed)
2693  {
2694  const double tol = 0.5 * cell->face (face)->line (line)->diameter () / cell->get_fe ().degree;
2695  const unsigned int dim = 3;
2696  const unsigned int spacedim = 3;
2697 
2698  hp_fe_values.reinit
2699  (cell,
2700  (cell->active_fe_index () * GeometryInfo<dim>::faces_per_cell + face)
2702 
2703  // Initialize the required
2704  // objects.
2705  const FEValues<dim> &
2706  fe_values = hp_fe_values.get_present_fe_values ();
2707  const FiniteElement<dim> &fe = cell->get_fe ();
2708  const std::vector< DerivativeForm<1,dim,spacedim> > &
2709  jacobians = fe_values.get_jacobians ();
2710  const std::vector<Point<dim> > &
2711  quadrature_points = fe_values.get_quadrature_points ();
2712 
2713  std::vector<Point<dim> > tangentials (fe_values.n_quadrature_points);
2714  std::vector<Vector<double> > values (fe_values.n_quadrature_points,
2715  Vector<double> (fe.n_components ()));
2716 
2717  // Get boundary function values
2718  // at quadrature points.
2719  boundary_function.vector_value_list (quadrature_points, values);
2720 
2721  const std::vector<Point<dim> > &
2722  reference_quadrature_points = fe_values.get_quadrature ().get_points ();
2723  std::pair<unsigned int, unsigned int> base_indices (0, 0);
2724 
2725  if (dynamic_cast<const FESystem<dim>*> (&cell->get_fe ()) != 0)
2726  {
2727  unsigned int fe_index = 0;
2728  unsigned int fe_index_old = 0;
2729  unsigned int i = 0;
2730 
2731  for (; i < fe.n_base_elements (); ++i)
2732  {
2733  fe_index_old = fe_index;
2734  fe_index += fe.element_multiplicity (i) * fe.base_element (i).n_components ();
2735 
2736  if (fe_index >= first_vector_component)
2737  break;
2738  }
2739 
2740  base_indices.first = i;
2741  base_indices.second = (first_vector_component - fe_index_old) / fe.base_element (i).n_components ();
2742  }
2743 
2744  // coordinate directions of
2745  // the edges of the face.
2746  const unsigned int
2747  edge_coordinate_direction
2750  = { { 2, 2, 1, 1 },
2751  { 2, 2, 1, 1 },
2752  { 0, 0, 2, 2 },
2753  { 0, 0, 2, 2 },
2754  { 1, 1, 0, 0 },
2755  { 1, 1, 0, 0 }
2756  };
2757  const FEValuesExtractors::Vector vec (first_vector_component);
2758 
2759  // The interpolation for the
2760  // lowest order edge shape
2761  // functions is just the mean
2762  // value of the tangential
2763  // components of the boundary
2764  // function on the edge.
2765  for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points;
2766  ++q_point)
2767  {
2768  // Therefore compute the
2769  // tangential of the edge at
2770  // the quadrature point.
2771  Point<dim> shifted_reference_point_1 = reference_quadrature_points[q_point];
2772  Point<dim> shifted_reference_point_2 = reference_quadrature_points[q_point];
2773 
2774  shifted_reference_point_1 (edge_coordinate_direction[face][line]) += tol;
2775  shifted_reference_point_2 (edge_coordinate_direction[face][line]) -= tol;
2776  tangentials[q_point]
2777  = (0.5 *
2778  (fe_values.get_mapping ()
2779  .transform_unit_to_real_cell (cell,
2780  shifted_reference_point_1)
2781  -
2782  fe_values.get_mapping ()
2783  .transform_unit_to_real_cell (cell,
2784  shifted_reference_point_2))
2785  / tol);
2786  tangentials[q_point]
2787  /= std::sqrt (tangentials[q_point].square ());
2788 
2789  // Compute the degrees of
2790  // freedom.
2791  for (unsigned int i = 0; i < fe.dofs_per_face; ++i)
2792  if (((dynamic_cast<const FESystem<dim>*> (&fe) != 0)
2793  && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).first == base_indices)
2794  && (fe.base_element (base_indices.first).face_to_cell_index (line * fe.degree, face)
2795  <= fe.system_to_base_index (fe.face_to_cell_index (i, face)).second)
2796  && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).second
2797  <= fe.base_element (base_indices.first).face_to_cell_index
2798  ((line + 1) * fe.degree - 1, face)))
2799  || ((dynamic_cast<const FE_Nedelec<dim>*> (&fe) != 0) && (line * fe.degree <= i)
2800  && (i < (line + 1) * fe.degree)))
2801  {
2802  dof_values[i]
2803  += (fe_values.JxW (q_point)
2804  * (values[q_point] (first_vector_component) * tangentials[q_point] (0)
2805  + values[q_point] (first_vector_component + 1) * tangentials[q_point] (1)
2806  + values[q_point] (first_vector_component + 2) * tangentials[q_point] (2))
2807  * (fe_values[vec].value (fe.face_to_cell_index (i, face), q_point) * tangentials[q_point])
2808  / std::sqrt (jacobians[q_point][0][edge_coordinate_direction[face][line]]
2809  * jacobians[q_point][0][edge_coordinate_direction[face][line]]
2810  + jacobians[q_point][1][edge_coordinate_direction[face][line]]
2811  * jacobians[q_point][1][edge_coordinate_direction[face][line]]
2812  + jacobians[q_point][2][edge_coordinate_direction[face][line]]
2813  * jacobians[q_point][2][edge_coordinate_direction[face][line]]));
2814 
2815  if (q_point == 0)
2816  dofs_processed[i] = true;
2817  }
2818  }
2819  }
2820 
2821  // dummy implementation of above
2822  // function for all other
2823  // dimensions
2824  template<int dim, typename cell_iterator>
2825  void
2826  compute_edge_projection (const cell_iterator &,
2827  const unsigned int,
2828  const unsigned int,
2830  const Function<dim> &,
2831  const unsigned int,
2832  std::vector<double> &,
2833  std::vector<bool> &)
2834  {
2835  Assert (false, ExcInternalError ());
2836  }
2837 
2838  // This function computes the
2839  // projection of the boundary
2840  // function on the interior of
2841  // faces.
2842  template<int dim, typename cell_iterator>
2843  void
2844  compute_face_projection_curl_conforming (const cell_iterator &cell,
2845  const unsigned int face,
2846  hp::FEValues<dim> &hp_fe_values,
2847  const Function<dim> &boundary_function,
2848  const unsigned int first_vector_component,
2849  std::vector<double> &dof_values,
2850  std::vector<bool> &dofs_processed)
2851  {
2852  const unsigned int spacedim = dim;
2853  hp_fe_values.reinit (cell, cell->active_fe_index ()
2855  // Initialize the required
2856  // objects.
2857  const FEValues<dim> &
2858  fe_values = hp_fe_values.get_present_fe_values ();
2859  const FiniteElement<dim> &fe = cell->get_fe ();
2860  const std::vector< DerivativeForm<1,dim,spacedim> > &
2861  jacobians = fe_values.get_jacobians ();
2862  const std::vector<Point<dim> > &
2863  quadrature_points = fe_values.get_quadrature_points ();
2864  const unsigned int degree = fe.degree - 1;
2865  std::pair<unsigned int, unsigned int> base_indices (0, 0);
2866 
2867  if (dynamic_cast<const FESystem<dim>*> (&cell->get_fe ()) != 0)
2868  {
2869  unsigned int fe_index = 0;
2870  unsigned int fe_index_old = 0;
2871  unsigned int i = 0;
2872 
2873  for (; i < fe.n_base_elements (); ++i)
2874  {
2875  fe_index_old = fe_index;
2876  fe_index += fe.element_multiplicity (i) * fe.base_element (i).n_components ();
2877 
2878  if (fe_index >= first_vector_component)
2879  break;
2880  }
2881 
2882  base_indices.first = i;
2883  base_indices.second = (first_vector_component - fe_index_old) / fe.base_element (i).n_components ();
2884  }
2885 
2886  std::vector<Vector<double> >
2887  values (fe_values.n_quadrature_points, Vector<double> (fe.n_components ()));
2888 
2889  // Get boundary function
2890  // values at quadrature
2891  // points.
2892  boundary_function.vector_value_list (quadrature_points, values);
2893 
2894  switch (dim)
2895  {
2896  case 2:
2897  {
2898  const double tol = 0.5 * cell->face (face)->diameter () / cell->get_fe ().degree;
2899  std::vector<Point<dim> >
2900  tangentials (fe_values.n_quadrature_points);
2901 
2902  const std::vector<Point<dim> > &
2903  reference_quadrature_points = fe_values.get_quadrature ().get_points ();
2904 
2905  // coordinate directions
2906  // of the face.
2907  const unsigned int
2908  face_coordinate_direction[GeometryInfo<dim>::faces_per_cell]
2909  = { 1, 1, 0, 0 };
2910  const FEValuesExtractors::Vector vec (first_vector_component);
2911 
2912  // The interpolation for
2913  // the lowest order face
2914  // shape functions is just
2915  // the mean value of the
2916  // tangential components
2917  // of the boundary function
2918  // on the edge.
2919  for (unsigned int q_point = 0;
2920  q_point < fe_values.n_quadrature_points; ++q_point)
2921  {
2922  // Therefore compute the
2923  // tangential of the
2924  // face at the quadrature
2925  // point.
2926  Point<dim> shifted_reference_point_1
2927  = reference_quadrature_points[q_point];
2928  Point<dim> shifted_reference_point_2
2929  = reference_quadrature_points[q_point];
2930 
2931  shifted_reference_point_1 (face_coordinate_direction[face])
2932  += tol;
2933  shifted_reference_point_2 (face_coordinate_direction[face])
2934  -= tol;
2935  tangentials[q_point]
2936  = (fe_values.get_mapping ()
2937  .transform_unit_to_real_cell (cell,
2938  shifted_reference_point_1)
2939  -
2940  fe_values.get_mapping ()
2941  .transform_unit_to_real_cell (cell,
2942  shifted_reference_point_2))
2943  / tol;
2944  tangentials[q_point]
2945  /= std::sqrt (tangentials[q_point].square ());
2946  // Compute the degrees
2947  // of freedom.
2948  for (unsigned int i = 0; i < fe.dofs_per_face; ++i)
2949  if (((dynamic_cast<const FESystem<dim>*> (&fe) != 0)
2950  && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).first == base_indices))
2951  || (dynamic_cast<const FE_Nedelec<dim>*> (&fe) != 0))
2952  {
2953  dof_values[i]
2954  += fe_values.JxW (q_point)
2955  * (values[q_point] (first_vector_component)
2956  * tangentials[q_point] (0)
2957  + values[q_point] (first_vector_component + 1)
2958  * tangentials[q_point] (1))
2959  * (fe_values[vec].value (fe.face_to_cell_index (i, face), q_point)
2960  * tangentials[q_point])
2961  / std::sqrt (jacobians[q_point][0][face_coordinate_direction[face]]
2962  * jacobians[q_point][0][face_coordinate_direction[face]]
2963  + jacobians[q_point][1][face_coordinate_direction[face]]
2964  * jacobians[q_point][1][face_coordinate_direction[face]]);
2965 
2966  if (q_point == 0)
2967  dofs_processed[i] = true;
2968  }
2969  }
2970 
2971  break;
2972  }
2973 
2974  case 3:
2975  {
2976  const FEValuesExtractors::Vector vec (first_vector_component);
2978  assembling_matrix (degree * fe.degree,
2979  dim * fe_values.n_quadrature_points);
2980  Vector<double> assembling_vector (assembling_matrix.n ());
2981  Vector<double> cell_rhs (assembling_matrix.m ());
2982  FullMatrix<double> cell_matrix (assembling_matrix.m (),
2983  assembling_matrix.m ());
2984  FullMatrix<double> cell_matrix_inv (assembling_matrix.m (),
2985  assembling_matrix.m ());
2986  Vector<double> solution (cell_matrix.m ());
2987 
2988  // Get coordinate directions
2989  // of the face.
2990  const unsigned int
2991  global_face_coordinate_directions[GeometryInfo<3>::faces_per_cell][2]
2992  = { { 1, 2 },
2993  { 1, 2 },
2994  { 2, 0 },
2995  { 2, 0 },
2996  { 0, 1 },
2997  { 0, 1 }
2998  };
2999 
3000  // The projection is divided into two steps. In the first step we
3001  // project the boundary function on the horizontal shape functions.
3002  // Then the boundary function is projected on the vertical shape
3003  // functions. We begin with the horizontal shape functions and set
3004  // up a linear system of equations to get the values for degrees of
3005  // freedom associated with the interior of the face.
3006  for (unsigned int q_point = 0;
3007  q_point < fe_values.n_quadrature_points; ++q_point)
3008  {
3009  // The right hand
3010  // side of the
3011  // corresponding problem
3012  // is the residual
3013  // of the boundary
3014  // function and
3015  // the already
3016  // interpolated part
3017  // on the edges.
3018  Tensor<1, dim> tmp;
3019 
3020  for (unsigned int d = 0; d < dim; ++d)
3021  tmp[d] = values[q_point] (first_vector_component + d);
3022 
3023  for (unsigned int i = 0; i < fe.dofs_per_face; ++i)
3024  if (((dynamic_cast<const FESystem<dim>*> (&fe) != 0)
3025  && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).first == base_indices)
3026  && (fe.base_element (base_indices.first).face_to_cell_index (2 * fe.degree, face)
3027  <= fe.system_to_base_index (fe.face_to_cell_index (i, face)).second)
3028  && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).second
3029  < fe.base_element (base_indices.first).face_to_cell_index (4 * fe.degree, face)))
3030  || ((dynamic_cast<const FE_Nedelec<dim>*> (&fe) != 0) && (2 * fe.degree <= i)
3031  && (i < 4 * fe.degree)))
3032  tmp -= dof_values[i] * fe_values[vec].value (fe.face_to_cell_index (i, face), q_point);
3033 
3034  const double JxW
3035  = std::sqrt (fe_values.JxW (q_point)
3036  / ((jacobians[q_point][0][global_face_coordinate_directions[face][0]]
3037  * jacobians[q_point][0][global_face_coordinate_directions[face][0]]
3038  +
3039  jacobians[q_point][1][global_face_coordinate_directions[face][0]]
3040  * jacobians[q_point][1][global_face_coordinate_directions[face][0]]
3041  +
3042  jacobians[q_point][2][global_face_coordinate_directions[face][0]]
3043  * jacobians[q_point][2][global_face_coordinate_directions[face][0]])
3044  *
3045  (jacobians[q_point][0][global_face_coordinate_directions[face][1]]
3046  * jacobians[q_point][0][global_face_coordinate_directions[face][1]]
3047  +
3048  jacobians[q_point][1][global_face_coordinate_directions[face][1]]
3049  * jacobians[q_point][1][global_face_coordinate_directions[face][1]]
3050  +
3051  jacobians[q_point][2][global_face_coordinate_directions[face][1]]
3052  * jacobians[q_point][2][global_face_coordinate_directions[face][1]])));
3053 
3054  // In the weak form
3055  // the right hand
3056  // side function
3057  // is multiplicated
3058  // by the horizontal
3059  // shape functions
3060  // defined in the
3061  // interior of
3062  // the face.
3063  for (unsigned int d = 0; d < dim; ++d)
3064  assembling_vector (dim * q_point + d) = JxW * tmp[d];
3065 
3066  unsigned int index = 0;
3067 
3068  for (unsigned int i = 0; i < fe.dofs_per_face; ++i)
3069  if (((dynamic_cast<const FESystem<dim>*> (&fe) != 0)
3070  && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).first == base_indices)
3071  && (fe.base_element (base_indices.first).face_to_cell_index
3073  <= fe.system_to_base_index (fe.face_to_cell_index (i, face)).second)
3074  && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).second
3075  < fe.base_element (base_indices.first).face_to_cell_index
3076  ((degree + GeometryInfo<dim>::lines_per_face) * fe.degree, face)))
3077  || ((dynamic_cast<const FE_Nedelec<dim>*> (&fe) != 0)
3079  && (i < (degree + GeometryInfo<dim>::lines_per_face) * fe.degree)))
3080  {
3081  const Tensor<1, dim> shape_value
3082  = (JxW * fe_values[vec].value (fe.face_to_cell_index (i, face),
3083  q_point));
3084 
3085  for (unsigned int d = 0; d < dim; ++d)
3086  assembling_matrix (index, dim * q_point + d) = shape_value[d];
3087 
3088  ++index;
3089  }
3090  }
3091 
3092  // Create the system matrix by multiplying the assembling matrix
3093  // with its transposed and the right hand side vector by mutliplying
3094  // the assembling matrix with the assembling vector. Invert the
3095  // system matrix.
3096  assembling_matrix.mTmult (cell_matrix, assembling_matrix);
3097  cell_matrix_inv.invert (cell_matrix);
3098  assembling_matrix.vmult (cell_rhs, assembling_vector);
3099  cell_matrix_inv.vmult (solution, cell_rhs);
3100 
3101  // Store the computed
3102  // values.
3103  {
3104  unsigned int index = 0;
3105 
3106  for (unsigned int i = 0; i < fe.dofs_per_face; ++i)
3107  if (((dynamic_cast<const FESystem<dim>*> (&fe) != 0)
3108  && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).first == base_indices)
3109  && (fe.base_element (base_indices.first).face_to_cell_index
3111  <= fe.system_to_base_index (fe.face_to_cell_index (i, face)).second)
3112  && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).second
3113  < fe.base_element (base_indices.first).face_to_cell_index
3114  ((degree + GeometryInfo<dim>::lines_per_face) * fe.degree, face)))
3115  || ((dynamic_cast<const FE_Nedelec<dim>*> (&fe) != 0)
3117  && (i < (degree + GeometryInfo<dim>::lines_per_face) * fe.degree)))
3118  {
3119  dof_values[i] = solution (index);
3120  dofs_processed[i] = true;
3121  ++index;
3122  }
3123  }
3124 
3125  // Now we do the same as above with the vertical shape functions
3126  // instead of the horizontal ones.
3127  for (unsigned int q_point = 0;
3128  q_point < fe_values.n_quadrature_points; ++q_point)
3129  {
3130  Tensor<1, dim> tmp;
3131 
3132  for (unsigned int d = 0; d < dim; ++d)
3133  tmp[d] = values[q_point] (first_vector_component + d);
3134 
3135  for (unsigned int i = 0; i < fe.dofs_per_face; ++i)
3136  if (((dynamic_cast<const FESystem<dim>*> (&fe) != 0)
3137  && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).first == base_indices)
3138  && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).second
3139  < fe.base_element (base_indices.first).face_to_cell_index (2 * fe.degree, face)))
3140  || ((dynamic_cast<const FE_Nedelec<dim>*> (&fe) != 0) && (i < 2 * fe.degree)))
3141  tmp -= dof_values[i] * fe_values[vec].value (fe.face_to_cell_index (i, face), q_point);
3142 
3143  const double JxW
3144  = std::sqrt (fe_values.JxW (q_point)
3145  / ((jacobians[q_point][0][global_face_coordinate_directions[face][0]]
3146  * jacobians[q_point][0][global_face_coordinate_directions[face][0]]
3147  +
3148  jacobians[q_point][1][global_face_coordinate_directions[face][0]]
3149  * jacobians[q_point][1][global_face_coordinate_directions[face][0]]
3150  +
3151  jacobians[q_point][2][global_face_coordinate_directions[face][0]]
3152  * jacobians[q_point][2][global_face_coordinate_directions[face][0]])
3153  *
3154  (jacobians[q_point][0][global_face_coordinate_directions[face][1]]
3155  * jacobians[q_point][0][global_face_coordinate_directions[face][1]]
3156  +
3157  jacobians[q_point][1][global_face_coordinate_directions[face][1]]
3158  * jacobians[q_point][1][global_face_coordinate_directions[face][1]]
3159  +
3160  jacobians[q_point][2][global_face_coordinate_directions[face][1]]
3161  * jacobians[q_point][2][global_face_coordinate_directions[face][1]])));
3162 
3163  for (unsigned int d = 0; d < dim; ++d)
3164  assembling_vector (dim * q_point + d) = JxW * tmp[d];
3165 
3166  unsigned int index = 0;
3167 
3168  for (unsigned int i = 0; i < fe.dofs_per_face; ++i)
3169  if (((dynamic_cast<const FESystem<dim>*> (&fe) != 0)
3170  && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).first == base_indices)
3171  && (fe.base_element (base_indices.first).face_to_cell_index
3172  ((degree + GeometryInfo<dim>::lines_per_face) * fe.degree, face)
3173  <= fe.system_to_base_index (fe.face_to_cell_index (i, face)).second))
3174  || ((dynamic_cast<const FE_Nedelec<dim>*> (&fe) != 0)
3175  && ((degree + GeometryInfo<dim>::lines_per_face) * fe.degree <= i)))
3176  {
3177  const Tensor<1, dim> shape_value
3178  = JxW * fe_values[vec].value (fe.face_to_cell_index (i, face),
3179  q_point);
3180 
3181  for (unsigned int d = 0; d < dim; ++d)
3182  assembling_matrix (index, dim * q_point + d) = shape_value[d];
3183 
3184  ++index;
3185  }
3186  }
3187 
3188  assembling_matrix.mTmult (cell_matrix, assembling_matrix);
3189  cell_matrix_inv.invert (cell_matrix);
3190  assembling_matrix.vmult (cell_rhs, assembling_vector);
3191  cell_matrix_inv.vmult (solution, cell_rhs);
3192 
3193  unsigned int index = 0;
3194 
3195  for (unsigned int i = 0; i < fe.dofs_per_face; ++i)
3196  if (((dynamic_cast<const FESystem<dim>*> (&fe) != 0)
3197  && (fe.system_to_base_index (fe.face_to_cell_index (i, face)).first == base_indices)
3198  && (fe.base_element (base_indices.first).face_to_cell_index
3199  ((degree + GeometryInfo<dim>::lines_per_face) * fe.degree, face)
3200  <= fe.system_to_base_index (fe.face_to_cell_index (i, face)).second))
3201  || ((dynamic_cast<const FE_Nedelec<dim>*> (&fe) != 0)
3202  && ((degree + GeometryInfo<dim>::lines_per_face) * fe.degree <= i)))
3203  {
3204  dof_values[i] = solution (index);
3205  dofs_processed[i] = true;
3206  ++index;
3207  }
3208 
3209  break;
3210  }
3211 
3212  default:
3213  Assert (false, ExcNotImplemented ());
3214  }
3215  }
3216  }
3217 
3218 
3219 
3220 
3221  template <int dim>
3222  void
3223 
3225  const unsigned int first_vector_component,
3226  const Function<dim> &boundary_function,
3227  const types::boundary_id boundary_component,
3228  ConstraintMatrix &constraints,
3229  const Mapping<dim> &mapping)
3230  {
3231  // Projection-based interpolation is performed in two (in 2D) respectively
3232  // three (in 3D) steps. First the tangential component of the function is
3233  // interpolated on each edge. This gives the values for the degrees of
3234  // freedom corresponding to the edge shape functions. Now we are done for
3235  // 2D, but in 3D we possibly have also degrees of freedom, which are
3236  // located in the interior of the faces. Therefore we compute the residual
3237  // of the function describing the boundary values and the interpolated
3238  // part, which we have computed in the last step. On the faces there are
3239  // two kinds of shape functions, the horizontal and the vertical
3240  // ones. Thus we have to solve two linear systems of equations of size
3241  // <tt>degree * (degree + 1)<tt> to obtain the values for the
3242  // corresponding degrees of freedom.
3243  const unsigned int superdegree = dof_handler.get_fe ().degree;
3244  const QGauss<dim - 1> reference_face_quadrature (2 * superdegree);
3245  const unsigned int dofs_per_face = dof_handler.get_fe ().dofs_per_face;
3246  hp::FECollection<dim> fe_collection (dof_handler.get_fe ());
3247  hp::MappingCollection<dim> mapping_collection (mapping);
3248  hp::QCollection<dim> face_quadrature_collection;
3249 
3250  for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
3251  face_quadrature_collection.push_back
3252  (QProjector<dim>::project_to_face (reference_face_quadrature, face));
3253 
3254  hp::FEValues<dim> fe_face_values (mapping_collection, fe_collection,
3255  face_quadrature_collection,
3259  update_values);
3260 
3261  std::vector<bool> dofs_processed (dofs_per_face);
3262  std::vector<double> dof_values (dofs_per_face);
3263  std::vector<types::global_dof_index> face_dof_indices (dofs_per_face);
3264  typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active ();
3265 
3266  switch (dim)
3267  {
3268  case 2:
3269  {
3270  for (; cell != dof_handler.end (); ++cell)
3271  if (cell->at_boundary () && cell->is_locally_owned ())
3272  for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
3273  if (cell->face (face)->boundary_indicator () == boundary_component)
3274  {
3275  // if the FE is a
3276  // FE_Nothing object
3277  // there is no work to
3278  // do
3279  if (dynamic_cast<const FE_Nothing<dim>*> (&cell->get_fe ()) != 0)
3280  return;
3281 
3282  // This is only
3283  // implemented, if the
3284  // FE is a Nedelec
3285  // element. If the FE
3286  // is a FESystem, we
3287  // cannot check this.
3288  if (dynamic_cast<const FESystem<dim>*> (&cell->get_fe ()) == 0)
3289  {
3290  typedef FiniteElement<dim> FEL;
3291  AssertThrow (dynamic_cast<const FE_Nedelec<dim>*> (&cell->get_fe ()) != 0,
3292 
3293  typename FEL::ExcInterpolationNotImplemented ());
3294  }
3295 
3296  for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
3297  {
3298  dof_values[dof] = 0.0;
3299  dofs_processed[dof] = false;
3300  }
3301 
3302  // Compute the
3303  // projection of the
3304  // boundary function on
3305  // the edge.
3306  internals
3307  ::compute_face_projection_curl_conforming (cell, face, fe_face_values,
3308  boundary_function,
3309  first_vector_component,
3310  dof_values, dofs_processed);
3311  cell->face (face)->get_dof_indices (face_dof_indices,
3312  cell->active_fe_index ());
3313 
3314  // Add the computed
3315  // constraints to the
3316  // constraint matrix,
3317  // if the degree of
3318  // freedom is not
3319  // already constrained.
3320  for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
3321  if (dofs_processed[dof] && constraints.can_store_line (face_dof_indices[dof])
3322  && !(constraints.is_constrained (face_dof_indices[dof])))
3323  {
3324  constraints.add_line (face_dof_indices[dof]);
3325 
3326  if (std::abs (dof_values[dof]) > 1e-13)
3327  constraints.set_inhomogeneity (face_dof_indices[dof], dof_values[dof]);
3328  }
3329  }
3330 
3331  break;
3332  }
3333 
3334  case 3:
3335  {
3336  const QGauss<dim - 2> reference_edge_quadrature (2 * superdegree);
3337  const unsigned int degree = superdegree - 1;
3338  hp::QCollection<dim> edge_quadrature_collection;
3339 
3340  for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
3341  for (unsigned int line = 0; line < GeometryInfo<dim>::lines_per_face; ++line)
3342  edge_quadrature_collection.push_back
3345  (reference_edge_quadrature, line), face));
3346 
3347  hp::FEValues<dim> fe_edge_values (mapping_collection, fe_collection,
3348  edge_quadrature_collection,
3352  update_values);
3353 
3354  for (; cell != dof_handler.end (); ++cell)
3355  if (cell->at_boundary () && cell->is_locally_owned ())
3356  for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
3357  if (cell->face (face)->boundary_indicator () == boundary_component)
3358  {
3359  // if the FE is a
3360  // FE_Nothing object
3361  // there is no work to
3362  // do
3363  if (dynamic_cast<const FE_Nothing<dim>*> (&cell->get_fe ()) != 0)
3364  return;
3365 
3366  // This is only
3367  // implemented, if the
3368  // FE is a Nedelec
3369  // element. If the FE is
3370  // a FESystem we cannot
3371  // check this.
3372  if (dynamic_cast<const FESystem<dim>*> (&cell->get_fe ()) == 0)
3373  {
3374  typedef FiniteElement<dim> FEL;
3375 
3376  AssertThrow (dynamic_cast<const FE_Nedelec<dim>*> (&cell->get_fe ()) != 0,
3377  typename FEL::ExcInterpolationNotImplemented ());
3378  }
3379 
3380  for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
3381  {
3382  dof_values[dof] = 0.0;
3383  dofs_processed[dof] = false;
3384  }
3385 
3386  // First we compute the
3387  // projection on the
3388  // edges.
3389  for (unsigned int line = 0;
3390  line < GeometryInfo<3>::lines_per_face; ++line)
3391  internals
3392  ::compute_edge_projection (cell, face, line, fe_edge_values,
3393  boundary_function,
3394  first_vector_component,
3395  dof_values, dofs_processed);
3396 
3397  // If there are higher
3398  // order shape
3399  // functions, there is
3400  // still some work
3401  // left.
3402  if (degree > 0)
3403  internals
3404  ::compute_face_projection_curl_conforming (cell, face, fe_face_values,
3405  boundary_function,
3406  first_vector_component,
3407  dof_values,
3408  dofs_processed);
3409 
3410  // Store the computed
3411  // values in the global
3412  // vector.
3413  cell->face (face)->get_dof_indices (face_dof_indices,
3414  cell->active_fe_index ());
3415 
3416  for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
3417  if (dofs_processed[dof] && constraints.can_store_line (face_dof_indices[dof])
3418  && !(constraints.is_constrained (face_dof_indices[dof])))
3419  {
3420  constraints.add_line (face_dof_indices[dof]);
3421 
3422  if (std::abs (dof_values[dof]) > 1e-13)
3423  constraints.set_inhomogeneity (face_dof_indices[dof], dof_values[dof]);
3424  }
3425  }
3426 
3427  break;
3428  }
3429 
3430  default:
3431  Assert (false, ExcNotImplemented ());
3432  }
3433  }
3434 
3435 
3436 
3437  template <int dim>
3438  void
3439 
3441  const unsigned int first_vector_component,
3442  const Function<dim> &boundary_function,
3443  const types::boundary_id boundary_component,
3444  ConstraintMatrix &constraints,
3445  const hp::MappingCollection<dim> &mapping_collection)
3446  {
3447  hp::FECollection<dim> fe_collection (dof_handler.get_fe ());
3448  hp::QCollection<dim> face_quadrature_collection;
3449 
3450  for (unsigned int i = 0; i < fe_collection.size (); ++i)
3451  {
3452  const QGauss<dim - 1>
3453  reference_face_quadrature (2 * fe_collection[i].degree);
3454 
3455  for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
3456  face_quadrature_collection.push_back
3457  (QProjector<dim>::project_to_face (reference_face_quadrature, face));
3458  }
3459 
3460  hp::FEValues<dim> fe_face_values (mapping_collection, fe_collection,
3461  face_quadrature_collection,
3465  update_values);
3466  std::vector<bool> dofs_processed;
3467  std::vector<double> dof_values;
3468  std::vector<types::global_dof_index> face_dof_indices;
3469  typename hp::DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active ();
3470 
3471  switch (dim)
3472  {
3473  case 2:
3474  {
3475  for (; cell != dof_handler.end (); ++cell)
3476  if (cell->at_boundary () && cell->is_locally_owned ())
3477  for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
3478  if (cell->face (face)->boundary_indicator () == boundary_component)
3479  {
3480  // if the FE is a FE_Nothing object there is no work to do
3481  if (dynamic_cast<const FE_Nothing<dim>*> (&cell->get_fe ()) != 0)
3482  return;
3483 
3484  // This is only implemented, if the FE is a Nedelec
3485  // element. If the FE is a FESystem we cannot check this.
3486  if (dynamic_cast<const FESystem<dim>*> (&cell->get_fe ()) == 0)
3487  {
3488  typedef FiniteElement<dim> FEL;
3489 
3490  AssertThrow (dynamic_cast<const FE_Nedelec<dim>*> (&cell->get_fe ()) != 0,
3491  typename FEL::ExcInterpolationNotImplemented ());
3492  }
3493 
3494  const unsigned int dofs_per_face = cell->get_fe ().dofs_per_face;
3495 
3496  dofs_processed.resize (dofs_per_face);
3497  dof_values.resize (dofs_per_face);
3498 
3499  for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
3500  {
3501  dof_values[dof] = 0.0;
3502  dofs_processed[dof] = false;
3503  }
3504 
3505  internals
3506  ::compute_face_projection_curl_conforming (cell, face, fe_face_values,
3507  boundary_function,
3508  first_vector_component,
3509  dof_values, dofs_processed);
3510  face_dof_indices.resize (dofs_per_face);
3511  cell->face (face)->get_dof_indices (face_dof_indices,
3512  cell->active_fe_index ());
3513 
3514  for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
3515  if (dofs_processed[dof] && constraints.can_store_line (face_dof_indices[dof])
3516  && !(constraints.is_constrained (face_dof_indices[dof])))
3517  {
3518  constraints.add_line (face_dof_indices[dof]);
3519 
3520  if (std::abs (dof_values[dof]) > 1e-13)
3521  constraints.set_inhomogeneity (face_dof_indices[dof], dof_values[dof]);
3522  }
3523  }
3524 
3525  break;
3526  }
3527 
3528  case 3:
3529  {
3530  hp::QCollection<dim> edge_quadrature_collection;
3531 
3532  for (unsigned int i = 0; i < fe_collection.size (); ++i)
3533  {
3534  const QGauss<dim - 2>
3535  reference_edge_quadrature (2 * fe_collection[i].degree);
3536 
3537  for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
3538  for (unsigned int line = 0; line < GeometryInfo<dim>::lines_per_face; ++line)
3539  edge_quadrature_collection.push_back
3541  (QProjector<dim - 1>::project_to_face (reference_edge_quadrature, line),
3542  face));
3543  }
3544 
3545  hp::FEValues<dim> fe_edge_values (mapping_collection, fe_collection,
3546  edge_quadrature_collection,
3550  update_values);
3551 
3552  for (; cell != dof_handler.end (); ++cell)
3553  if (cell->at_boundary () && cell->is_locally_owned ())
3554  for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
3555  if (cell->face (face)->boundary_indicator () == boundary_component)
3556  {
3557  // if the FE is a FE_Nothing object there is no work to do
3558  if (dynamic_cast<const FE_Nothing<dim>*> (&cell->get_fe ()) != 0)
3559  return;
3560 
3561  // This is only implemented, if the FE is a Nedelec
3562  // element. If the FE is a FESystem we cannot check this.
3563  if (dynamic_cast<const FESystem<dim>*> (&cell->get_fe ()) == 0)
3564  {
3565  typedef FiniteElement<dim> FEL;
3566 
3567  AssertThrow (dynamic_cast<const FE_Nedelec<dim>*> (&cell->get_fe ()) != 0,
3568  typename FEL::ExcInterpolationNotImplemented ());
3569  }
3570 
3571  const unsigned int superdegree = cell->get_fe ().degree;
3572  const unsigned int degree = superdegree - 1;
3573  const unsigned int dofs_per_face = cell->get_fe ().dofs_per_face;
3574 
3575  dofs_processed.resize (dofs_per_face);
3576  dof_values.resize (dofs_per_face);
3577 
3578  for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
3579  {
3580  dof_values[dof] = 0.0;
3581  dofs_processed[dof] = false;
3582  }
3583 
3584  for (unsigned int line = 0;
3585  line < GeometryInfo<dim>::lines_per_face; ++line)
3586  internals
3587  ::compute_edge_projection (cell, face, line, fe_edge_values,
3588  boundary_function,
3589  first_vector_component,
3590  dof_values, dofs_processed);
3591 
3592  // If there are higher order shape functions, there is still
3593  // some work left.
3594  if (degree > 0)
3595  internals
3596  ::compute_face_projection_curl_conforming (cell, face, fe_face_values,
3597  boundary_function,
3598  first_vector_component,
3599  dof_values, dofs_processed);
3600 
3601 
3602  face_dof_indices.resize (dofs_per_face);
3603  cell->face (face)->get_dof_indices (face_dof_indices,
3604  cell->active_fe_index ());
3605 
3606  for (unsigned int dof = 0; dof < dofs_per_face; ++dof)
3607  if (dofs_processed[dof] && constraints.can_store_line (face_dof_indices[dof])
3608  && !(constraints.is_constrained (face_dof_indices[dof])))
3609  {
3610  constraints.add_line (face_dof_indices[dof]);
3611 
3612  if (std::abs (dof_values[dof]) > 1e-13)
3613  constraints.set_inhomogeneity (face_dof_indices[dof], dof_values[dof]);
3614  }
3615  }
3616 
3617  break;
3618  }
3619 
3620  default:
3621  Assert (false, ExcNotImplemented ());
3622  }
3623  }
3624 
3625 
3626  namespace internals
3627  {
3628  // This function computes the projection of the boundary function on the
3629  // boundary in 2d.
3630  template <typename cell_iterator>
3631  void
3632  compute_face_projection_div_conforming (const cell_iterator &cell,
3633  const unsigned int face,
3634  const FEFaceValues<2> &fe_values,
3635  const unsigned int first_vector_component,
3636  const Function<2> &boundary_function,
3637  const std::vector<DerivativeForm<1,2,2> > &jacobians,
3638  ConstraintMatrix &constraints)
3639  {
3640  // Compute the intergral over the product of the normal components of
3641  // the boundary function times the normal components of the shape
3642  // functions supported on the boundary.
3643  const FEValuesExtractors::Vector vec (first_vector_component);
3644  const FiniteElement<2> &fe = cell->get_fe ();
3645  const std::vector<Point<2> > &normals = fe_values.get_normal_vectors ();
3646  const unsigned int
3647  face_coordinate_direction[GeometryInfo<2>::faces_per_cell] = {1, 1, 0, 0};
3648  std::vector<Vector<double> >
3649  values (fe_values.n_quadrature_points, Vector<double> (2));
3650  Vector<double> dof_values (fe.dofs_per_face);
3651 
3652  // Get the values of the boundary function at the quadrature points.
3653  {
3654  const std::vector<Point<2> > &
3655  quadrature_points = fe_values.get_quadrature_points ();
3656 
3657  boundary_function.vector_value_list (quadrature_points, values);
3658  }
3659 
3660  for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point)
3661  {
3662  double tmp = 0.0;
3663 
3664  for (unsigned int d = 0; d < 2; ++d)
3665  tmp += normals[q_point][d] * values[q_point] (d);
3666 
3667  tmp *= fe_values.JxW (q_point)
3668  * std::sqrt (jacobians[q_point][0][face_coordinate_direction[face]]
3669  * jacobians[q_point][0][face_coordinate_direction[face]]
3670  + jacobians[q_point][1][face_coordinate_direction[face]]
3671  * jacobians[q_point][1][face_coordinate_direction[face]]);
3672 
3673  for (unsigned int i = 0; i < fe.dofs_per_face; ++i)
3674  dof_values (i) += tmp * (normals[q_point]
3675  * fe_values[vec].value (fe.face_to_cell_index (i, face), q_point));
3676  }
3677 
3678  std::vector<types::global_dof_index> face_dof_indices (fe.dofs_per_face);
3679 
3680  cell->face (face)->get_dof_indices (face_dof_indices, cell->active_fe_index ());
3681 
3682  // Copy the computed values in the ConstraintMatrix only, if the degree
3683  // of freedom is not already constrained.
3684  for (unsigned int i = 0; i < fe.dofs_per_face; ++i)
3685  if (!(constraints.is_constrained (face_dof_indices[i])))
3686  {
3687  constraints.add_line (face_dof_indices[i]);
3688 
3689  if (std::abs (dof_values (i)) > 1e-14)
3690  constraints.set_inhomogeneity (face_dof_indices[i], dof_values (i));
3691  }
3692  }
3693 
3694  // dummy implementation of above function for all other dimensions
3695  template<int dim, typename cell_iterator>
3696  void
3697  compute_face_projection_div_conforming (const cell_iterator &,
3698  const unsigned int,
3699  const FEFaceValues<dim> &,
3700  const unsigned int,
3701  const Function<dim> &,
3702  const std::vector<DerivativeForm<1,dim,dim> > &,
3703  ConstraintMatrix &)
3704  {
3705  Assert (false, ExcNotImplemented ());
3706  }
3707 
3708  // This function computes the projection of the boundary function on the
3709  // boundary in 3d.
3710  template<typename cell_iterator>
3711  void
3712  compute_face_projection_div_conforming (const cell_iterator &cell,
3713  const unsigned int face,
3714  const FEFaceValues<3> &fe_values,
3715  const unsigned int first_vector_component,
3716  const Function<3> &boundary_function,
3717  const std::vector<DerivativeForm<1,3,3> > &jacobians,
3718  std::vector<double> &dof_values,
3719  std::vector<types::global_dof_index> &projected_dofs)
3720  {
3721  // Compute the intergral over the product of the normal components of
3722  // the boundary function times the normal components of the shape
3723  // functions supported on the boundary.
3724  const FEValuesExtractors::Vector vec (first_vector_component);
3725  const FiniteElement<3> &fe = cell->get_fe ();
3726  const std::vector<Point<3> > &normals = fe_values.get_normal_vectors ();
3727  const unsigned int
3728  face_coordinate_directions[GeometryInfo<3>::faces_per_cell][2] = {{1, 2},
3729  {1, 2},
3730  {2, 0},
3731  {2, 0},
3732  {0, 1},
3733  {0, 1}
3734  };
3735  std::vector<Vector<double> >
3736  values (fe_values.n_quadrature_points, Vector<double> (3));
3737  Vector<double> dof_values_local (fe.dofs_per_face);
3738 
3739  {
3740  const std::vector<Point<3> > &
3741  quadrature_points = fe_values.get_quadrature_points ();
3742 
3743  boundary_function.vector_value_list (quadrature_points, values);
3744  }
3745 
3746  for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; ++q_point)
3747  {
3748  double tmp = 0.0;
3749 
3750  for (unsigned int d = 0; d < 3; ++d)
3751  tmp += normals[q_point][d] * values[q_point] (d);
3752 
3753  tmp *= fe_values.JxW (q_point)
3754  * std::sqrt ((jacobians[q_point][0][face_coordinate_directions[face][0]]
3755  * jacobians[q_point][0][face_coordinate_directions[face][0]]
3756  + jacobians[q_point][1][face_coordinate_directions[face][0]]
3757  * jacobians[q_point][1][face_coordinate_directions[face][0]]
3758  + jacobians[q_point][2][face_coordinate_directions[face][0]]
3759  * jacobians[q_point][2][face_coordinate_directions[face][0]])
3760  * (jacobians[q_point][0][face_coordinate_directions[face][1]]
3761  * jacobians[q_point][0][face_coordinate_directions[face][1]]
3762  + jacobians[q_point][1][face_coordinate_directions[face][1]]
3763  * jacobians[q_point][1][face_coordinate_directions[face][1]]
3764  + jacobians[q_point][2][face_coordinate_directions[face][1]]
3765  * jacobians[q_point][2][face_coordinate_directions[face][1]]));
3766 
3767  for (unsigned int i = 0; i < fe.dofs_per_face; ++i)
3768  dof_values_local (i) += tmp * (normals[q_point]
3769  * fe_values[vec].value (fe.face_to_cell_index (i, face), q_point));
3770  }
3771 
3772  std::vector<types::global_dof_index> face_dof_indices (fe.dofs_per_face);
3773 
3774  cell->face (face)->get_dof_indices (face_dof_indices, cell->active_fe_index ());
3775 
3776  for (unsigned int i = 0; i < fe.dofs_per_face; ++i)
3777  if (projected_dofs[face_dof_indices[i]] < fe.degree)
3778  {
3779  dof_values[face_dof_indices[i]] = dof_values_local (i);
3780  projected_dofs[face_dof_indices[i]] = fe.degree;
3781  }
3782  }
3783 
3784  // dummy implementation of above
3785  // function for all other
3786  // dimensions
3787  template<int dim, typename cell_iterator>
3788  void
3789  compute_face_projection_div_conforming (const cell_iterator &,
3790  const unsigned int,
3791  const FEFaceValues<dim> &,
3792  const unsigned int,
3793  const Function<dim> &,
3794  const std::vector<DerivativeForm<1,dim,dim> > &,
3795  std::vector<double> &,
3796  std::vector<types::global_dof_index> &)
3797  {
3798  Assert (false, ExcNotImplemented ());
3799  }
3800  }
3801 
3802 
3803  template <int dim>
3804  void
3806  const unsigned int first_vector_component,
3807  const Function<dim> &boundary_function,
3808  const types::boundary_id boundary_component,
3809  ConstraintMatrix &constraints,
3810  const Mapping<dim> &mapping)
3811  {
3812  const unsigned int spacedim = dim;
3813  // Interpolate the normal components
3814  // of the boundary functions. Since
3815  // the Raviart-Thomas elements are
3816  // constructed from a Lagrangian
3817  // basis, it suffices to compute
3818  // the integral over the product
3819  // of the normal components of the
3820  // boundary function times the
3821  // normal components of the shape
3822  // functions supported on the
3823  // boundary.
3824  const FiniteElement<dim> &fe = dof_handler.get_fe ();
3825  QGauss<dim - 1> face_quadrature (fe.degree + 1);
3826  FEFaceValues<dim> fe_face_values (mapping, fe, face_quadrature, update_JxW_values |
3829  update_values);
3830  hp::FECollection<dim> fe_collection (fe);
3831  hp::MappingCollection<dim> mapping_collection (mapping);
3832  hp::QCollection<dim> quadrature_collection;
3833 
3834  for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
3835  quadrature_collection.push_back (QProjector<dim>::project_to_face (face_quadrature,
3836  face));
3837 
3838  hp::FEValues<dim> fe_values (mapping_collection, fe_collection, quadrature_collection,
3840 
3841  switch (dim)
3842  {
3843  case 2:
3844  {
3845  for (typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active ();
3846  cell != dof_handler.end (); ++cell)
3847  if (cell->at_boundary ())
3848  for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
3849  if (cell->face (face)->boundary_indicator () == boundary_component)
3850  {
3851  // if the FE is a
3852  // FE_Nothing object
3853  // there is no work to
3854  // do
3855  if (dynamic_cast<const FE_Nothing<dim>*> (&cell->get_fe ()) != 0)
3856  return;
3857 
3858  // This is only
3859  // implemented, if the
3860  // FE is a Raviart-Thomas
3861  // element. If the FE is
3862  // a FESystem we cannot
3863  // check this.
3864  if (dynamic_cast<const FESystem<dim>*> (&cell->get_fe ()) == 0)
3865  {
3866  typedef FiniteElement<dim> FEL;
3867 
3868  AssertThrow (dynamic_cast<const FE_RaviartThomas<dim>*> (&cell->get_fe ()) != 0,
3869  typename FEL::ExcInterpolationNotImplemented ());
3870  }
3871 
3872  fe_values.reinit (cell, face + cell->active_fe_index ()
3874 
3875  const std::vector<DerivativeForm<1,dim,spacedim> > &
3876  jacobians = fe_values.get_present_fe_values ().get_jacobians ();
3877 
3878  fe_face_values.reinit (cell, face);
3879  internals::compute_face_projection_div_conforming (cell, face,
3880  fe_face_values,
3881  first_vector_component,
3882  boundary_function,
3883  jacobians,
3884  constraints);
3885  }
3886 
3887  break;
3888  }
3889 
3890  case 3:
3891  {
3892  // In three dimensions the
3893  // edges between two faces
3894  // are treated twice.
3895  // Therefore we store the
3896  // computed values in a
3897  // vector and copy them over
3898  // in the ConstraintMatrix
3899  // after all values have been
3900  // computed.
3901  // If we have two values for
3902  // one edge, we choose the one,
3903  // which was computed with the
3904  // higher order element.
3905  // If both elements are of the
3906  // same order, we just keep the
3907  // first value and do not
3908  // compute a second one.
3909  const unsigned int &n_dofs = dof_handler.n_dofs ();
3910  std::vector<double> dof_values (n_dofs);
3911  std::vector<types::global_dof_index> projected_dofs (n_dofs);
3912 
3913  for (unsigned int dof = 0; dof < n_dofs; ++dof)
3914  projected_dofs[dof] = 0;
3915 
3916  for (typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active ();
3917  cell != dof_handler.end (); ++cell)
3918  if (cell->at_boundary ())
3919  for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
3920  if (cell->face (face)->boundary_indicator () == boundary_component)
3921  {
3922  // This is only
3923  // implemented, if the
3924  // FE is a Raviart-Thomas
3925  // element. If the FE is
3926  // a FESystem we cannot
3927  // check this.
3928  if (dynamic_cast<const FESystem<dim>*> (&cell->get_fe ()) == 0)
3929  {
3930  typedef FiniteElement<dim> FEL;
3931 
3932  AssertThrow (dynamic_cast<const FE_RaviartThomas<dim>*> (&cell->get_fe ()) != 0,
3933  typename FEL::ExcInterpolationNotImplemented ());
3934  }
3935 
3936  fe_values.reinit (cell, face + cell->active_fe_index ()
3938 
3939  const std::vector<DerivativeForm<1,dim ,spacedim> > &
3940  jacobians = fe_values.get_present_fe_values ().get_jacobians ();
3941 
3942  fe_face_values.reinit (cell, face);
3943  internals::compute_face_projection_div_conforming (cell, face,
3944  fe_face_values,
3945  first_vector_component,
3946  boundary_function,
3947  jacobians, dof_values,
3948  projected_dofs);
3949  }
3950 
3951  for (unsigned int dof = 0; dof < n_dofs; ++dof)
3952  if ((projected_dofs[dof] != 0) && !(constraints.is_constrained (dof)))
3953  {
3954  constraints.add_line (dof);
3955 
3956  if (std::abs (dof_values[dof]) > 1e-14)
3957  constraints.set_inhomogeneity (dof, dof_values[dof]);
3958  }
3959 
3960  break;
3961  }
3962 
3963  default:
3964  Assert (false, ExcNotImplemented ());
3965  }
3966  }
3967 
3968 
3969  template <int dim>
3970  void
3972  const unsigned int first_vector_component,
3973  const Function<dim> &boundary_function,
3974  const types::boundary_id boundary_component,
3975  ConstraintMatrix &constraints,
3976  const hp::MappingCollection<dim, dim> &mapping_collection)
3977  {
3978  const unsigned int spacedim = dim;
3979  const hp::FECollection<dim> &fe_collection = dof_handler.get_fe ();
3980  hp::QCollection<dim - 1> face_quadrature_collection;
3981  hp::QCollection<dim> quadrature_collection;
3982 
3983  for (unsigned int i = 0; i < fe_collection.size (); ++i)
3984  {
3985  const QGauss<dim - 1> quadrature (fe_collection[i].degree + 1);
3986 
3987  face_quadrature_collection.push_back (quadrature);
3988 
3989  for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
3990  quadrature_collection.push_back (QProjector<dim>::project_to_face (quadrature,
3991  face));
3992  }
3993 
3994  hp::FEFaceValues<dim> fe_face_values (mapping_collection, fe_collection,
3995  face_quadrature_collection, update_JxW_values |
3998  update_values);
3999  hp::FEValues<dim> fe_values (mapping_collection, fe_collection, quadrature_collection,
4001 
4002  switch (dim)
4003  {
4004  case 2:
4005  {
4006  for (typename hp::DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active ();
4007  cell != dof_handler.end (); ++cell)
4008  if (cell->at_boundary ())
4009  for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
4010  if (cell->face (face)->boundary_indicator () == boundary_component)
4011  {
4012  // This is only
4013  // implemented, if the
4014  // FE is a Raviart-Thomas
4015  // element. If the FE is
4016  // a FESystem we cannot
4017  // check this.
4018  if (dynamic_cast<const FESystem<dim>*> (&cell->get_fe ()) == 0)
4019  {
4020  typedef FiniteElement<dim> FEL;
4021 
4022  AssertThrow (dynamic_cast<const FE_RaviartThomas<dim>*> (&cell->get_fe ()) != 0,
4023  typename FEL::ExcInterpolationNotImplemented ());
4024  }
4025 
4026  fe_values.reinit (cell, face + cell->active_fe_index ()
4028 
4029  const std::vector<DerivativeForm<1,dim,spacedim> > &
4030  jacobians = fe_values.get_present_fe_values ().get_jacobians ();
4031 
4032  fe_face_values.reinit (cell, face);
4033  internals::compute_face_projection_div_conforming (cell, face,
4034  fe_face_values.get_present_fe_values (),
4035  first_vector_component,
4036  boundary_function,
4037  jacobians,
4038  constraints);
4039  }
4040 
4041  break;
4042  }
4043 
4044  case 3:
4045  {
4046  const unsigned int &n_dofs = dof_handler.n_dofs ();
4047  std::vector<double> dof_values (n_dofs);
4048  std::vector<types::global_dof_index> projected_dofs (n_dofs);
4049 
4050  for (unsigned int dof = 0; dof < n_dofs; ++dof)
4051  projected_dofs[dof] = 0;
4052 
4053  for (typename hp::DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active ();
4054  cell != dof_handler.end (); ++cell)
4055  if (cell->at_boundary ())
4056  for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
4057  if (cell->face (face)->boundary_indicator () == boundary_component)
4058  {
4059  // This is only
4060  // implemented, if the
4061  // FE is a Raviart-Thomas
4062  // element. If the FE is
4063  // a FESystem we cannot
4064  // check this.
4065  if (dynamic_cast<const FESystem<dim>*> (&cell->get_fe ()) == 0)
4066  {
4067  typedef FiniteElement<dim> FEL;
4068 
4069  AssertThrow (dynamic_cast<const FE_RaviartThomas<dim>*> (&cell->get_fe ()) != 0,
4070  typename FEL::ExcInterpolationNotImplemented ());
4071  }
4072 
4073  fe_values.reinit (cell, face + cell->active_fe_index ()
4075 
4076  const std::vector<DerivativeForm<1,dim,spacedim> > &
4077  jacobians = fe_values.get_present_fe_values ().get_jacobians ();
4078 
4079  fe_face_values.reinit (cell, face);
4080  internals::compute_face_projection_div_conforming (cell, face,
4081  fe_face_values.get_present_fe_values (),
4082  first_vector_component,
4083  boundary_function,
4084  jacobians, dof_values,
4085  projected_dofs);
4086  }
4087 
4088  for (unsigned int dof = 0; dof < n_dofs; ++dof)
4089  if ((projected_dofs[dof] != 0) && !(constraints.is_constrained (dof)))
4090  {
4091  constraints.add_line (dof);
4092 
4093  if (std::abs (dof_values[dof]) > 1e-14)
4094  constraints.set_inhomogeneity (dof, dof_values[dof]);
4095  }
4096 
4097  break;
4098  }
4099 
4100  default:
4101  Assert (false, ExcNotImplemented ());
4102  }
4103  }
4104 
4105 
4106 
4107  template <int dim, template <int, int> class DH, int spacedim>
4108  void
4109  compute_no_normal_flux_constraints (const DH<dim,spacedim> &dof_handler,
4110  const unsigned int first_vector_component,
4111  const std::set<types::boundary_id> &boundary_ids,
4112  ConstraintMatrix &constraints,
4113  const Mapping<dim, spacedim> &mapping)
4114  {
4115  Assert (dim > 1,
4116  ExcMessage ("This function is not useful in 1d because it amounts "
4117  "to imposing Dirichlet values on the vector-valued "
4118  "quantity."));
4119 
4120  std::vector<types::global_dof_index> face_dofs;
4121 
4122  // create FE and mapping collections for all elements in use by this
4123  // DoFHandler
4124  hp::FECollection<dim,spacedim> fe_collection (dof_handler.get_fe());
4125  hp::MappingCollection<dim,spacedim> mapping_collection;
4126  for (unsigned int i=0; i<fe_collection.size(); ++i)
4127  mapping_collection.push_back (mapping);
4128 
4129  // now also create a quadrature collection for the faces of a cell. fill
4130  // it with a quadrature formula with the support points on faces for each
4131  // FE
4132  hp::QCollection<dim-1> face_quadrature_collection;
4133  for (unsigned int i=0; i<fe_collection.size(); ++i)
4134  {
4135  const std::vector<Point<dim-1> > &
4136  unit_support_points = fe_collection[i].get_unit_face_support_points();
4137 
4138  Assert (unit_support_points.size() == fe_collection[i].dofs_per_face,
4139  ExcInternalError());
4140 
4141  face_quadrature_collection
4142  .push_back (Quadrature<dim-1> (unit_support_points));
4143  }
4144 
4145  // now create the object with which we will generate the normal vectors
4146  hp::FEFaceValues<dim,spacedim> x_fe_face_values (mapping_collection,
4147  fe_collection,
4148  face_quadrature_collection,
4149  update_q_points |
4151 
4152  // have a map that stores normal vectors for each vector-dof tuple we want
4153  // to constrain. since we can get at the same vector dof tuple more than
4154  // once (for example if it is located at a vertex that we visit from all
4155  // adjacent cells), we will want to average later on the normal vectors
4156  // computed on different cells as described in the documentation of this
4157  // function. however, we can only average if the contributions came from
4158  // different cells, whereas we want to constrain twice or more in case the
4159  // contributions came from different faces of the same cell
4160  // (i.e. constrain not just the *average normal direction* but *all normal
4161  // directions* we find). consequently, we also have to store which cell a
4162  // normal vector was computed on
4163  typedef
4164  std::multimap<internal::VectorDoFTuple<dim>,
4165  std::pair<Tensor<1,dim>, typename DH<dim,spacedim>::active_cell_iterator> >
4166  DoFToNormalsMap;
4167 
4168  DoFToNormalsMap dof_to_normals_map;
4169 
4170  // now loop over all cells and all faces
4171  typename DH<dim,spacedim>::active_cell_iterator
4172  cell = dof_handler.begin_active(),
4173  endc = dof_handler.end();
4174  for (; cell!=endc; ++cell)
4175  if (!cell->is_artificial())
4176  for (unsigned int face_no=0; face_no < GeometryInfo<dim>::faces_per_cell;
4177  ++face_no)
4178  if (boundary_ids.find(cell->face(face_no)->boundary_indicator())
4179  != boundary_ids.end())
4180  {
4181  const FiniteElement<dim> &fe = cell->get_fe ();
4182  typename DH<dim,spacedim>::face_iterator face = cell->face(face_no);
4183 
4184  // get the indices of the dofs on this cell...
4185  face_dofs.resize (fe.dofs_per_face);
4186  face->get_dof_indices (face_dofs, cell->active_fe_index());
4187 
4188  x_fe_face_values.reinit (cell, face_no);
4189  const FEFaceValues<dim> &fe_values = x_fe_face_values.get_present_fe_values();
4190 
4191  // then identify which of them correspond to the selected set of
4192  // vector components
4193  for (unsigned int i=0; i<face_dofs.size(); ++i)
4194  if (fe.face_system_to_component_index(i).first ==
4195  first_vector_component)
4196  {
4197  // find corresponding other components of vector
4198  internal::VectorDoFTuple<dim> vector_dofs;
4199  vector_dofs.dof_indices[0] = face_dofs[i];
4200 
4201  Assert(first_vector_component+dim<=fe.n_components(),
4202  ExcMessage("Error: the finite element does not have enough components "
4203  "to define a normal direction."));
4204 
4205  for (unsigned int k=0; k<fe.dofs_per_face; ++k)
4206  if ((k != i)
4207  &&
4208  (face_quadrature_collection[cell->active_fe_index()].point(k) ==
4209  face_quadrature_collection[cell->active_fe_index()].point(i))
4210  &&
4211  (fe.face_system_to_component_index(k).first >=
4212  first_vector_component)
4213  &&
4214  (fe.face_system_to_component_index(k).first <
4215  first_vector_component + dim))
4216  vector_dofs.dof_indices[fe.face_system_to_component_index(k).first -
4217  first_vector_component]
4218  = face_dofs[k];
4219 
4220  for (unsigned int d=0; d<dim; ++d)
4221  Assert (vector_dofs.dof_indices[d] < dof_handler.n_dofs(),
4222  ExcInternalError());
4223 
4224  // we need the normal vector on this face. we know that it
4225  // is a vector of length 1 but at least with higher order
4226  // mappings it isn't always possible to guarantee that
4227  // each component is exact up to zero tolerance. in
4228  // particular, as shown in the deal.II/no_flux_06 test, if
4229  // we just take the normal vector as given by the
4230  // fe_values object, we can get entries in the normal
4231  // vectors of the unit cube that have entries up to
4232  // several times 1e-14.
4233  //
4234  // the problem with this is that this later yields
4235  // constraints that are circular (e.g., in the testcase,
4236  // we get constraints of the form
4237  //
4238  // x22 = 2.93099e-14*x21 + 2.93099e-14*x23
4239  // x21 = -2.93099e-14*x22 + 2.93099e-14*x21
4240  //
4241  // in both of these constraints, the small numbers should
4242  // be zero and the constraints should simply be
4243  // x22 = x21 = 0
4244  //
4245  // to achieve this, we utilize that we know that the
4246  // normal vector has (or should have) length 1 and that we
4247  // can simply set small elements to zero (without having
4248  // to check that they are small *relative to something
4249  // else*). we do this and then normalize the length of the
4250  // vector back to one, just to be on the safe side
4251  //
4252  // one more point: we would like to use the "real" normal
4253  // vector here, as provided by the boundary description
4254  // and as opposed to what we get from the FEValues object.
4255  // we do this in the immediately next line, but as is
4256  // obvious, the boundary only has a vague idea which side
4257  // of a cell it is on -- indicated by the face number. in
4258  // other words, it may provide the inner or outer normal.
4259  // by and large, there is no harm from this, since the
4260  // tangential vector we compute is still the same. however,
4261  // we do average over normal vectors from adjacent cells
4262  // and if they have recorded normal vectors from the inside
4263  // once and from the outside the other time, then this
4264  // averaging is going to run into trouble. as a consequence
4265  // we ask the mapping after all for its normal vector,
4266  // but we only ask it so that we can possibly correct the
4267  // sign of the normal vector provided by the boundary
4268  // if they should point in different directions. this is the
4269  // case in tests/deal.II/no_flux_11.
4270  Point<dim> normal_vector
4271  = (cell->face(face_no)->get_boundary()
4272  .normal_vector (cell->face(face_no),
4273  fe_values.quadrature_point(i)));
4274  if (normal_vector * fe_values.normal_vector(i) < 0)
4275  normal_vector *= -1;
4276  Assert (std::fabs(normal_vector.norm() - 1) < 1e-14,
4277  ExcInternalError());
4278  for (unsigned int d=0; d<dim; ++d)
4279  if (std::fabs(normal_vector[d]) < 1e-13)
4280  normal_vector[d] = 0;
4281  normal_vector /= normal_vector.norm();
4282 
4283  // now enter the (dofs,(normal_vector,cell)) entry into
4284  // the map
4285  dof_to_normals_map
4286  .insert (std::make_pair (vector_dofs,
4287  std::make_pair (normal_vector,
4288  cell)));
4289 #ifdef DEBUG_NO_NORMAL_FLUX
4290  std::cout << "Adding normal vector:" << std::endl
4291  << " dofs=" << vector_dofs << std::endl
4292  << " cell=" << cell << " at " << cell->center() << std::endl
4293  << " normal=" << normal_vector << std::endl;
4294 #endif
4295  }
4296  }
4297 
4298  // Now do something with the collected information. To this end, loop
4299  // through all sets of pairs (dofs,normal_vector) and identify which
4300  // entries belong to the same set of dofs and then do as described in the
4301  // documentation, i.e. either average the normal vector or don't for this
4302  // particular set of dofs
4303  typename DoFToNormalsMap::const_iterator
4304  p = dof_to_normals_map.begin();
4305 
4306  while (p != dof_to_normals_map.end())
4307  {
4308  // first find the range of entries in the multimap that corresponds to
4309  // the same vector-dof tuple. as usual, we define the range
4310  // half-open. the first entry of course is 'p'
4311  typename DoFToNormalsMap::const_iterator same_dof_range[2]
4312  = { p };
4313  for (++p; p != dof_to_normals_map.end(); ++p)
4314  if (p->first != same_dof_range[0]->first)
4315  {
4316  same_dof_range[1] = p;
4317  break;
4318  }
4319  if (p == dof_to_normals_map.end())
4320  same_dof_range[1] = dof_to_normals_map.end();
4321 
4322 #ifdef DEBUG_NO_NORMAL_FLUX
4323  std::cout << "For dof indices <" << p->first << ">, found the following normals"
4324  << std::endl;
4325  for (typename DoFToNormalsMap::const_iterator
4326  q = same_dof_range[0];
4327  q != same_dof_range[1]; ++q)
4328  std::cout << " " << q->second.first
4329  << " from cell " << q->second.second
4330  << std::endl;
4331 #endif
4332 
4333 
4334  // now compute the reverse mapping: for each of the cells that
4335  // contributed to the current set of vector dofs, add up the normal
4336  // vectors. the values of the map are pairs of normal vectors and
4337  // number of cells that have contributed
4338  typedef
4339  std::map
4340  <typename DH<dim,spacedim>::active_cell_iterator,
4341  std::pair<Tensor<1,dim>, unsigned int> >
4342  CellToNormalsMap;
4343 
4344  CellToNormalsMap cell_to_normals_map;
4345  for (typename DoFToNormalsMap::const_iterator
4346  q = same_dof_range[0];
4347  q != same_dof_range[1]; ++q)
4348  if (cell_to_normals_map.find (q->second.second)
4349  == cell_to_normals_map.end())
4350  cell_to_normals_map[q->second.second]
4351  = std::make_pair (q->second.first, 1U);
4352  else
4353  {
4354  const Tensor<1,dim> old_normal
4355  = cell_to_normals_map[q->second.second].first;
4356  const unsigned int old_count
4357  = cell_to_normals_map[q->second.second].second;
4358 
4359  Assert (old_count > 0, ExcInternalError());
4360 
4361  // in the same entry, store again the now averaged normal vector
4362  // and the new count
4363  cell_to_normals_map[q->second.second]
4364  = std::make_pair ((old_normal * old_count + q->second.first) / (old_count + 1),
4365  old_count + 1);
4366  }
4367  Assert (cell_to_normals_map.size() >= 1, ExcInternalError());
4368 
4369 #ifdef DEBUG_NO_NORMAL_FLUX
4370  std::cout << " cell_to_normals_map:" << std::endl;
4371  for (typename CellToNormalsMap::const_iterator
4372  x = cell_to_normals_map.begin();
4373  x != cell_to_normals_map.end(); ++x)
4374  std::cout << " " << x->first << " -> ("
4375  << x->second.first << ',' << x->second.second << ')'
4376  << std::endl;
4377 #endif
4378 
4379  // count the maximum number of contributions from each cell
4380  unsigned int max_n_contributions_per_cell = 1;
4381  for (typename CellToNormalsMap::const_iterator
4382  x = cell_to_normals_map.begin();
4383  x != cell_to_normals_map.end(); ++x)
4384  max_n_contributions_per_cell
4385  = std::max (max_n_contributions_per_cell,
4386  x->second.second);
4387 
4388  // verify that each cell can have only contributed at most dim times,
4389  // since that is the maximum number of faces that come together at a
4390  // single place
4391  Assert (max_n_contributions_per_cell <= dim, ExcInternalError());
4392 
4393  switch (max_n_contributions_per_cell)
4394  {
4395  // first deal with the case that a number of cells all have
4396  // registered that they have a normal vector defined at the
4397  // location of a given vector dof, and that each of them have
4398  // encountered this vector dof exactly once while looping over all
4399  // their faces. as stated in the documentation, this is the case
4400  // where we want to simply average over all normal vectors
4401  //
4402  // the typical case is in 2d where multiple cells meet at one
4403  // vertex sitting on the boundary. same in 3d for a vertex that
4404  // is associated with only one of the boundary indicators passed
4405  // to this function
4406  case 1:
4407  {
4408 
4409  // compute the average normal vector from all the ones that have
4410  // the same set of dofs. we could add them up and divide them by
4411  // the number of additions, or simply normalize them right away
4412  // since we want them to have unit length anyway
4413  Tensor<1,dim> normal;
4414  for (typename CellToNormalsMap::const_iterator
4415  x = cell_to_normals_map.begin();
4416  x != cell_to_normals_map.end(); ++x)
4417  normal += x->second.first;
4418  normal /= normal.norm();
4419 
4420  // normalize again
4421  for (unsigned int d=0; d<dim; ++d)
4422  if (std::fabs(normal[d]) < 1e-13)
4423  normal[d] = 0;
4424  normal /= normal.norm();
4425 
4426  // then construct constraints from this:
4428  dof_indices = same_dof_range[0]->first;
4429  internal::add_constraint (dof_indices, normal,
4430  constraints);
4431 
4432  break;
4433  }
4434 
4435 
4436  // this is the slightly more complicated case that a single cell has
4437  // contributed with exactly DIM normal vectors to the same set of
4438  // vector dofs. this is what happens in a corner in 2d and 3d (but
4439  // not on an edge in 3d, where we have only 2, i.e. <DIM,
4440  // contributions. Here we do not want to average the normal
4441  // vectors. Since we have DIM contributions, let's assume (and
4442  // verify) that they are in fact all linearly independent; in that
4443  // case, all vector components are constrained and we need to set
4444  // them to zero
4445  case dim:
4446  {
4447  // assert that indeed only a single cell has contributed
4448  Assert (cell_to_normals_map.size() == 1,
4449  ExcInternalError());
4450 
4451  // check linear independence by computing the determinant of the
4452  // matrix created from all the normal vectors. if they are
4453  // linearly independent, then the determinant is nonzero. if they
4454  // are orthogonal, then the matrix is in fact equal to 1 (since
4455  // they are all unit vectors); make sure the determinant is larger
4456  // than 1e-3 to avoid cases where cells are degenerate
4457  {
4458  Tensor<2,dim> t;
4459 
4460  typename DoFToNormalsMap::const_iterator x = same_dof_range[0];
4461  for (unsigned int i=0; i<dim; ++i, ++x)
4462  for (unsigned int j=0; j<dim; ++j)
4463  t[i][j] = x->second.first[j];
4464 
4465  Assert (std::fabs(determinant (t)) > 1e-3,
4466  ExcMessage("Found a set of normal vectors that are nearly collinear."));
4467  }
4468 
4469  // so all components of this vector dof are constrained. enter
4470  // this into the constraint matrix
4471  //
4472  // ignore dofs already constrained
4473  for (unsigned int i=0; i<dim; ++i)
4474  if (!constraints.is_constrained (same_dof_range[0]
4475  ->first.dof_indices[i])
4476  &&
4477  constraints.can_store_line(
4478  same_dof_range[0]->first.dof_indices[i]))
4479  {
4480  constraints.add_line (same_dof_range[0]->first.dof_indices[i]);
4481  // no add_entries here
4482  }
4483 
4484  break;
4485  }
4486 
4487 
4488  // this is the case of an edge contribution in 3d, i.e. the vector
4489  // is constrained in two directions but not the third.
4490  default:
4491  {
4492  Assert (dim >= 3, ExcNotImplemented());
4493  Assert (max_n_contributions_per_cell == 2, ExcInternalError());
4494 
4495  // as described in the documentation, let us first collect what
4496  // each of the cells contributed at the current point. we use a
4497  // std::list instead of a std::set (which would be more natural)
4498  // because std::set requires that the stored elements are
4499  // comparable with operator<
4500  typedef
4501  std::map<typename DH<dim,spacedim>::active_cell_iterator, std::list<Tensor<1,dim> > >
4502  CellContributions;
4503  CellContributions cell_contributions;
4504 
4505  for (typename DoFToNormalsMap::const_iterator
4506  q = same_dof_range[0];
4507  q != same_dof_range[1]; ++q)
4508  cell_contributions[q->second.second].push_back (q->second.first);
4509  Assert (cell_contributions.size() >= 1, ExcInternalError());
4510 
4511  // now for each cell that has contributed determine the number of
4512  // normal vectors it has contributed. we currently only implement
4513  // if this is dim-1 for all cells (if a single cell has
4514  // contributed dim, or if all adjacent cells have contributed 1
4515  // normal vector, this is already handled above).
4516  //
4517  // we only implement the case that all cells contribute
4518  // dim-1 because we assume that we are following an edge
4519  // of the domain (think: we are looking at a vertex
4520  // located on one of the edges of a refined cube where the
4521  // boundary indicators of the two adjacent faces of the
4522  // cube are both listed in the set of boundary indicators
4523  // passed to this function). in that case, all cells along
4524  // that edge of the domain are assumed to have contributed
4525  // dim-1 normal vectors. however, there are cases where
4526  // this assumption is not justified (see the lengthy
4527  // explanation in test no_flux_12.cc) and in those cases
4528  // we simply ignore the cell that contributes only
4529  // once. this is also discussed at length in the
4530  // documentation of this function.
4531  //
4532  // for each contributing cell compute the tangential vector that
4533  // remains unconstrained
4534  std::list<Tensor<1,dim> > tangential_vectors;
4535  for (typename CellContributions::const_iterator
4536  contribution = cell_contributions.begin();
4537  contribution != cell_contributions.end();
4538  ++contribution)
4539  {
4540 #ifdef DEBUG_NO_NORMAL_FLUX
4541  std::cout << " Treating edge case with dim-1 contributions." << std::endl
4542  << " Looking at cell " << contribution->first
4543  << " which has contributed these normal vectors:"
4544  << std::endl;
4545  for (typename std::list<Tensor<1,dim> >::const_iterator
4546  t = contribution->second.begin();
4547  t != contribution->second.end();
4548  ++t)
4549  std::cout << " " << *t << std::endl;
4550 #endif
4551 
4552  // as mentioned above, simply ignore cells that only
4553  // contribute once
4554  if (contribution->second.size() < dim-1)
4555  continue;
4556 
4557  Tensor<1,dim> normals[dim-1];
4558  {
4559  unsigned int index=0;
4560  for (typename std::list<Tensor<1,dim> >::const_iterator
4561  t = contribution->second.begin();
4562  t != contribution->second.end();
4563  ++t, ++index)
4564  normals[index] = *t;
4565  Assert (index == dim-1, ExcInternalError());
4566  }
4567 
4568  // calculate the tangent as the outer product of the normal
4569  // vectors. since these vectors do not need to be orthogonal
4570  // (think, for example, the case of the deal.II/no_flux_07
4571  // test: a sheared cube in 3d, with Q2 elements, where we have
4572  // constraints from the two normal vectors of two faces of the
4573  // sheared cube that are not perpendicular to each other), we
4574  // have to normalize the outer product
4575  Tensor<1,dim> tangent;
4576  switch (dim)
4577  {
4578  case 3:
4579  // take cross product between normals[0] and
4580  // normals[1]. write it in the current form (with [dim-2])
4581  // to make sure that compilers don't warn about
4582  // out-of-bounds accesses -- the warnings are bogus since
4583  // we get here only for dim==3, but at least one isn't
4584  // quite smart enough to notice this and warns when
4585  // compiling the function in 2d
4586  cross_product (tangent, normals[0], normals[dim-2]);
4587  break;
4588  default:
4589  Assert (false, ExcNotImplemented());
4590  }
4591 
4592  Assert (std::fabs (tangent.norm()) > 1e-12,
4593  ExcMessage("Two normal vectors from adjacent faces are almost "
4594  "parallel."));
4595  tangent /= tangent.norm();
4596 
4597  tangential_vectors.push_back (tangent);
4598  }
4599 
4600  // go through the list of tangents and make sure that they all
4601  // roughly point in the same direction as the first one (i.e. have
4602  // an angle less than 90 degrees); if they don't then flip their
4603  // sign
4604  {
4605  const Tensor<1,dim> first_tangent = tangential_vectors.front();
4606  typename std::list<Tensor<1,dim> >::iterator
4607  t = tangential_vectors.begin();
4608  ++t;
4609  for (; t != tangential_vectors.end(); ++t)
4610  if (*t * first_tangent < 0)
4611  *t *= -1;
4612  }
4613 
4614  // now compute the average tangent and normalize it
4615  Tensor<1,dim> average_tangent;
4616  for (typename std::list<Tensor<1,dim> >::const_iterator
4617  t = tangential_vectors.begin();
4618  t != tangential_vectors.end();
4619  ++t)
4620  average_tangent += *t;
4621  average_tangent /= average_tangent.norm();
4622 
4623  // now all that is left is that we add the constraints that the
4624  // vector is parallel to the tangent
4626  dof_indices = same_dof_range[0]->first;
4627  internal::add_tangentiality_constraints (dof_indices,
4628  average_tangent,
4629  constraints);
4630  }
4631  }
4632  }
4633  }
4634 
4635 
4636 
4637  namespace
4638  {
4639  template <int dim>
4640  struct PointComparator
4641  {
4642  bool operator ()(const std_cxx1x::array<types::global_dof_index,dim> &p1,
4643  const std_cxx1x::array<types::global_dof_index,dim> &p2)
4644  {
4645  for (unsigned int d=0; d<dim; ++d)
4646  if (p1[d] < p2[d])
4647  return true;
4648  return false;
4649  }
4650  };
4651  }
4652 
4653 
4654 
4655  template <int dim, template <int, int> class DH, int spacedim>
4656  void
4657  compute_normal_flux_constraints (const DH<dim,spacedim> &dof_handler,
4658  const unsigned int first_vector_component,
4659  const std::set<types::boundary_id> &boundary_ids,
4660  ConstraintMatrix &constraints,
4661  const Mapping<dim, spacedim> &mapping)
4662  {
4663  ConstraintMatrix no_normal_flux_constraints(constraints.get_local_lines());
4665  first_vector_component,
4666  boundary_ids,
4667  no_normal_flux_constraints,
4668  mapping);
4669 
4670  // Extract a list that collects all vector components that belong to the
4671  // same node (scalar basis function). When creating that list, we use an
4672  // array of dim components that stores the global degree of freedom.
4673  std::set<std_cxx1x::array<types::global_dof_index,dim>, PointComparator<dim> > vector_dofs;
4674  std::vector<types::global_dof_index> face_dofs;
4675 
4676  std::vector<std_cxx1x::array<types::global_dof_index,dim> > cell_vector_dofs;
4677  for (typename DH<dim,spacedim>::active_cell_iterator cell =
4678  dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
4679  if (!cell->is_artificial())
4680  for (unsigned int face_no=0; face_no < GeometryInfo<dim>::faces_per_cell;
4681  ++face_no)
4682  if (boundary_ids.find(cell->face(face_no)->boundary_indicator())
4683  != boundary_ids.end())
4684  {
4685  const FiniteElement<dim> &fe = cell->get_fe();
4686  typename DH<dim,spacedim>::face_iterator face=cell->face(face_no);
4687 
4688  // get the indices of the dofs on this cell...
4689  face_dofs.resize (fe.dofs_per_face);
4690  face->get_dof_indices (face_dofs, cell->active_fe_index());
4691 
4692  unsigned int n_scalar_indices = 0;
4693  cell_vector_dofs.resize(fe.dofs_per_face);
4694  for (unsigned int i=0; i<fe.dofs_per_face; ++i)
4695  if (fe.face_system_to_component_index(i).first >= first_vector_component &&
4696  fe.face_system_to_component_index(i).first < first_vector_component + dim)
4697  {
4698  n_scalar_indices =
4699  std::max(n_scalar_indices,
4700  fe.face_system_to_component_index(i).second+1);
4701  cell_vector_dofs[fe.face_system_to_component_index(i).second]
4702  [fe.face_system_to_component_index(i).first-first_vector_component]
4703  = face_dofs[i];
4704  }
4705 
4706  // now we identified the vector indices on the cell, so next
4707  // insert them into the set (it would be expensive to directly
4708  // insert incomplete points into the set)
4709  for (unsigned int i=0; i<n_scalar_indices; ++i)
4710  vector_dofs.insert(cell_vector_dofs[i]);
4711  }
4712 
4713  // iterate over the list of all vector components we found and see if we
4714  // can find constrained ones
4715  unsigned int n_total_constraints_found = 0;
4716  for (typename std::set<std_cxx1x::array<types::global_dof_index,dim>,
4717  PointComparator<dim> >::const_iterator it=vector_dofs.begin();
4718  it!=vector_dofs.end(); ++it)
4719  {
4720  unsigned int n_constraints = 0;
4721  bool is_constrained[dim];
4722  for (unsigned int d=0; d<dim; ++d)
4723  if (no_normal_flux_constraints.is_constrained((*it)[d]))
4724  {
4725  is_constrained[d] = true;
4726  ++n_constraints;
4727  ++n_total_constraints_found;
4728  }
4729  else
4730  is_constrained[d] = false;
4731  if (n_constraints > 0)
4732  {
4733  // if more than one no-flux constraint is present, we need to
4734  // constrain all vector degrees of freedom (we are in a corner
4735  // where several faces meet and to get a continuous FE solution we
4736  // need to set all conditions to zero).
4737  if (n_constraints > 1)
4738  {
4739  for (unsigned int d=0; d<dim; ++d)
4740  constraints.add_line((*it)[d]);
4741  continue;
4742  }
4743 
4744  // ok, this is a no-flux constraint, so get the index of the dof
4745  // that is currently constrained and make it unconstrained. The
4746  // constraint indices will get the normal that contain the other
4747  // indices.
4748  Tensor<1,dim> normal;
4749  unsigned constrained_index = -1;
4750  for (unsigned int d=0; d<dim; ++d)
4751  if (is_constrained[d])
4752  {
4753  constrained_index = d;
4754  normal[d] = 1.;
4755  }
4756  AssertIndexRange(constrained_index, dim);
4757  const std::vector<std::pair<types::global_dof_index, double> > *constrained
4758  = no_normal_flux_constraints.get_constraint_entries((*it)[constrained_index]);
4759  // find components to which this index is constrained to
4760  Assert(constrained != 0, ExcInternalError());
4761  Assert(constrained->size() < dim, ExcInternalError());
4762  for (unsigned int c=0; c<constrained->size(); ++c)
4763  {
4764  int index = -1;
4765  for (unsigned int d=0; d<dim; ++d)
4766  if ((*constrained)[c].first == (*it)[d])
4767  index = d;
4768  Assert (index != -1, ExcInternalError());
4769  normal[index] = (*constrained)[c].second;
4770  }
4771  for (unsigned int d=0; d<dim; ++d)
4772  {
4773  if (is_constrained[d])
4774  continue;
4775  const unsigned int new_index = (*it)[d];
4776  if (!constraints.is_constrained(new_index))
4777  {
4778  constraints.add_line(new_index);
4779  if (std::abs(normal[d]) > 1e-13)
4780  constraints.add_entry(new_index, (*it)[constrained_index],
4781  -normal[d]);
4782  }
4783  }
4784  }
4785  }
4786  AssertDimension(n_total_constraints_found,
4787  no_normal_flux_constraints.n_constraints());
4788  }
4789 
4790 
4791 
4792  namespace internal
4793  {
4794  template <int dim, int spacedim>
4796  {
4797  IDScratchData (const ::hp::MappingCollection<dim,spacedim> &mapping,
4798  const ::hp::FECollection<dim,spacedim> &fe,
4799  const ::hp::QCollection<dim> &q,
4800  const UpdateFlags update_flags);
4801 
4802  IDScratchData (const IDScratchData &data);
4803 
4804  void resize_vectors (const unsigned int n_q_points,
4805  const unsigned int n_components);
4806 
4807  std::vector<::Vector<double> > function_values;
4808  std::vector<std::vector<Tensor<1,spacedim> > > function_grads;
4809  std::vector<double> weight_values;
4810  std::vector<::Vector<double> > weight_vectors;
4811 
4812  std::vector<::Vector<double> > psi_values;
4813  std::vector<std::vector<Tensor<1,spacedim> > > psi_grads;
4814  std::vector<double> psi_scalar;
4815 
4816  std::vector<double> tmp_values;
4817  std::vector<Tensor<1,spacedim> > tmp_gradients;
4818 
4819  ::hp::FEValues<dim,spacedim> x_fe_values;
4820  };
4821 
4822 
4823  template <int dim, int spacedim>
4825  ::IDScratchData(const ::hp::MappingCollection<dim,spacedim> &mapping,
4826  const ::hp::FECollection<dim,spacedim> &fe,
4827  const ::hp::QCollection<dim> &q,
4828  const UpdateFlags update_flags)
4829  :
4830  x_fe_values(mapping, fe, q, update_flags)
4831  {}
4832 
4833  template <int dim, int spacedim>
4835  :
4836  x_fe_values(data.x_fe_values.get_mapping_collection(),
4837  data.x_fe_values.get_fe_collection(),
4838  data.x_fe_values.get_quadrature_collection(),
4839  data.x_fe_values.get_update_flags())
4840  {}
4841 
4842  template <int dim, int spacedim>
4843  void
4844  IDScratchData<dim,spacedim>::resize_vectors (const unsigned int n_q_points,
4845  const unsigned int n_components)
4846  {
4847  function_values.resize (n_q_points,
4848  ::Vector<double>(n_components));
4849  function_grads.resize (n_q_points,
4850  std::vector<Tensor<1,spacedim> >(n_components));
4851 
4852  weight_values.resize (n_q_points);
4853  weight_vectors.resize (n_q_points,
4854  ::Vector<double>(n_components));
4855 
4856  psi_values.resize (n_q_points,
4857  ::Vector<double>(n_components));
4858  psi_grads.resize (n_q_points,
4859  std::vector<Tensor<1,spacedim> >(n_components));
4860  psi_scalar.resize (n_q_points);
4861 
4862  tmp_values.resize (n_q_points);
4863  tmp_gradients.resize (n_q_points);
4864  }
4865 
4866 
4867  // avoid compiling inner function for many vector types when we always
4868  // really do the same thing by putting the main work into this helper
4869  // function
4870  template <int dim, int spacedim>
4871  double
4872  integrate_difference_inner (const Function<spacedim> &exact_solution,
4873  const NormType &norm,
4874  const Function<spacedim> *weight,
4875  const UpdateFlags update_flags,
4876  const double exponent,
4877  const unsigned int n_components,
4878  IDScratchData<dim,spacedim> &data)
4879  {
4880  const bool fe_is_system = (n_components != 1);
4881  const ::FEValues<dim, spacedim> &fe_values = data.x_fe_values.get_present_fe_values ();
4882  const unsigned int n_q_points = fe_values.n_quadrature_points;
4883 
4884  if (weight!=0)
4885  {
4886  if (weight->n_components>1)
4887  weight->vector_value_list (fe_values.get_quadrature_points(),
4888  data.weight_vectors);
4889  else
4890  {
4891  weight->value_list (fe_values.get_quadrature_points(),
4892  data.weight_values);
4893  for (unsigned int k=0; k<n_q_points; ++k)
4894  data.weight_vectors[k] = data.weight_values[k];
4895  }
4896  }
4897  else
4898  {
4899  for (unsigned int k=0; k<n_q_points; ++k)
4900  data.weight_vectors[k] = 1.;
4901  }
4902 
4903 
4904  if (update_flags & update_values)
4905  {
4906  // first compute the exact solution (vectors) at the quadrature
4907  // points try to do this as efficient as possible by avoiding a
4908  // second virtual function call in case the function really has only
4909  // one component
4910  if (fe_is_system)
4911  exact_solution.vector_value_list (fe_values.get_quadrature_points(),
4912  data.psi_values);
4913  else
4914  {
4915  exact_solution.value_list (fe_values.get_quadrature_points(),
4916  data.tmp_values);
4917  for (unsigned int i=0; i<n_q_points; ++i)
4918  data.psi_values[i](0) = data.tmp_values[i];
4919  }
4920 
4921  // then subtract finite element fe_function
4922  for (unsigned int q=0; q<n_q_points; ++q)
4923  data.psi_values[q] -= data.function_values[q];
4924  }
4925 
4926  // Do the same for gradients, if required
4927  if (update_flags & update_gradients)
4928  {
4929  // try to be a little clever to avoid recursive virtual function
4930  // calls when calling gradient_list for functions that are really
4931  // scalar functions
4932  if (fe_is_system)
4933  exact_solution.vector_gradient_list (fe_values.get_quadrature_points(),
4934  data.psi_grads);
4935  else
4936  {
4937  exact_solution.gradient_list (fe_values.get_quadrature_points(),
4938  data.tmp_gradients);
4939  for (unsigned int i=0; i<n_q_points; ++i)
4940  data.psi_grads[i][0] = data.tmp_gradients[i];
4941  }
4942 
4943  // then subtract finite element function_grads. We need to be
4944  // careful in the codimension one case, since there we only have
4945  // tangential gradients in the finite element function, not the full
4946  // gradient. This is taken care of, by subtracting the normal
4947  // component of the gradient from the exact function.
4948  if (update_flags & update_normal_vectors)
4949  for (unsigned int k=0; k<n_components; ++k)
4950  for (unsigned int q=0; q<n_q_points; ++q)
4951  data.psi_grads[q][k] -= (data.function_grads[q][k] +
4952  (data.psi_grads[q][k]* // (f.n) n
4953  fe_values.normal_vector(q))*
4954  fe_values.normal_vector(q));
4955  else
4956  for (unsigned int k=0; k<n_components; ++k)
4957  for (unsigned int q=0; q<n_q_points; ++q)
4958  data.psi_grads[q][k] -= data.function_grads[q][k];
4959  }
4960 
4961  double diff = 0;
4962  switch (norm)
4963  {
4964  case mean:
4965  // Compute values in quadrature points and integrate
4966  for (unsigned int q=0; q<n_q_points; ++q)
4967  {
4968  double sum = 0;
4969  for (unsigned int k=0; k<n_components; ++k)
4970  sum += data.psi_values[q](k) * data.weight_vectors[q](k);
4971  diff += sum * fe_values.JxW(q);
4972  }
4973  break;
4974 
4975  case Lp_norm:
4976  case L1_norm:
4977  case W1p_norm:
4978  // Compute values in quadrature points and integrate
4979  for (unsigned int q=0; q<n_q_points; ++q)
4980  {
4981  double sum = 0;
4982  for (unsigned int k=0; k<n_components; ++k)
4983  sum += std::pow(data.psi_values[q](k)*data.psi_values[q](k),
4984  exponent/2.) * data.weight_vectors[q](k);
4985  diff += sum * fe_values.JxW(q);
4986  }
4987 
4988  // Compute the root only, if no derivative values are added later
4989  if (!(update_flags & update_gradients))
4990  diff = std::pow(diff, 1./exponent);
4991  break;
4992 
4993  case L2_norm:
4994  case H1_norm:
4995  // Compute values in quadrature points and integrate
4996  for (unsigned int q=0; q<n_q_points; ++q)
4997  {
4998  double sum = 0;
4999  for (unsigned int k=0; k<n_components; ++k)
5000  sum += data.psi_values[q](k) * data.psi_values[q](k) *
5001  data.weight_vectors[q](k);
5002  diff += sum * fe_values.JxW(q);
5003  }
5004  // Compute the root only, if no derivative values are added later
5005  if (norm == L2_norm)
5006  diff=std::sqrt(diff);
5007  break;
5008 
5009  case Linfty_norm:
5010  case W1infty_norm:
5011  for (unsigned int q=0; q<n_q_points; ++q)
5012  for (unsigned int k=0; k<n_components; ++k)
5013  diff = std::max (diff, std::abs(data.psi_values[q](k)*
5014  data.weight_vectors[q](k)));
5015  break;
5016 
5017  case H1_seminorm:
5018  case W1p_seminorm:
5019  case W1infty_seminorm:
5020  break;
5021 
5022  default:
5023  Assert (false, ExcNotImplemented());
5024  break;
5025  }
5026 
5027  switch (norm)
5028  {
5029  case W1p_seminorm:
5030  case W1p_norm:
5031  for (unsigned int q=0; q<n_q_points; ++q)
5032  {
5033  double sum = 0;
5034  for (unsigned int k=0; k<n_components; ++k)
5035  sum += std::pow(data.psi_grads[q][k]*data.psi_grads[q][k],
5036  exponent/2.) * data.weight_vectors[q](k);
5037  diff += sum * fe_values.JxW(q);
5038  }
5039  diff = std::pow(diff, 1./exponent);
5040  break;
5041 
5042  case H1_seminorm:
5043  case H1_norm:
5044  for (unsigned int q=0; q<n_q_points; ++q)
5045  {
5046  double sum = 0;
5047  for (unsigned int k=0; k<n_components; ++k)
5048  sum += (data.psi_grads[q][k] * data.psi_grads[q][k]) *
5049  data.weight_vectors[q](k);
5050  diff += sum * fe_values.JxW(q);
5051  }
5052  diff = std::sqrt(diff);
5053  break;
5054 
5055  case W1infty_seminorm:
5056  case W1infty_norm:
5057  {
5058  double t = 0;
5059  for (unsigned int q=0; q<n_q_points; ++q)
5060  for (unsigned int k=0; k<n_components; ++k)
5061  for (unsigned int d=0; d<dim; ++d)
5062  t = std::max(t, std::abs(data.psi_grads[q][k][d]) *
5063  data.weight_vectors[q](k));
5064 
5065  // then add seminorm to norm if that had previously been computed
5066  diff += t;
5067  }
5068  break;
5069  default:
5070  break;
5071  }
5072 
5073  // append result of this cell to the end of the vector
5075  return diff;
5076  }
5077 
5078 
5079 
5080  template <int dim, class InVector, class OutVector, class DH, int spacedim>
5081  static
5082  void
5083  do_integrate_difference (const ::hp::MappingCollection<dim,spacedim> &mapping,
5084  const DH &dof,
5085  const InVector &fe_function,
5086  const Function<spacedim> &exact_solution,
5087  OutVector &difference,
5088  const ::hp::QCollection<dim> &q,
5089  const NormType &norm,
5090  const Function<spacedim> *weight,
5091  const double exponent_1)
5092  {
5093  // we mark the "exponent" parameter to this function "const" since it is
5094  // strictly incoming, but we need to set it to something different later
5095  // on, if necessary, so have a read-write version of it:
5096  double exponent = exponent_1;
5097 
5098  const unsigned int n_components = dof.get_fe().n_components();
5099 
5100  if (weight!=0)
5101  {
5102  Assert ((weight->n_components==1) || (weight->n_components==n_components),
5103  ExcDimensionMismatch(weight->n_components, n_components));
5104  }
5105 
5106  difference.reinit (dof.get_tria().n_active_cells());
5107 
5108  switch (norm)
5109  {
5110  case L2_norm:
5111  case H1_seminorm:
5112  case H1_norm:
5113  exponent = 2.;
5114  break;
5115 
5116  case L1_norm:
5117  exponent = 1.;
5118  break;
5119 
5120  default:
5121  break;
5122  }
5123 
5126  switch (norm)
5127  {
5128  case H1_seminorm:
5129  case W1p_seminorm:
5130  case W1infty_seminorm:
5131  update_flags |= UpdateFlags (update_gradients);
5132  if (spacedim == dim+1)
5133  update_flags |= UpdateFlags (update_normal_vectors);
5134 
5135  break;
5136 
5137  case H1_norm:
5138  case W1p_norm:
5139  case W1infty_norm:
5140  update_flags |= UpdateFlags (update_gradients);
5141  if (spacedim == dim+1)
5142  update_flags |= UpdateFlags (update_normal_vectors);
5143  // no break!
5144 
5145  default:
5146  update_flags |= UpdateFlags (update_values);
5147  break;
5148  }
5149 
5150  ::hp::FECollection<dim,spacedim> fe_collection (dof.get_fe());
5151  IDScratchData<dim,spacedim> data(mapping, fe_collection, q, update_flags);
5152 
5153  // loop over all cells
5154  typename DH::active_cell_iterator cell = dof.begin_active(),
5155  endc = dof.end();
5156  for (unsigned int index=0; cell != endc; ++cell, ++index)
5157  if (cell->is_locally_owned())
5158  {
5159  // initialize for this cell
5160  data.x_fe_values.reinit (cell);
5161 
5162  const ::FEValues<dim, spacedim> &fe_values = data.x_fe_values.get_present_fe_values ();
5163  const unsigned int n_q_points = fe_values.n_quadrature_points;
5164  data.resize_vectors (n_q_points, n_components);
5165 
5166  if (update_flags & update_values)
5167  fe_values.get_function_values (fe_function, data.function_values);
5168  if (update_flags & update_gradients)
5169  fe_values.get_function_grads (fe_function, data.function_grads);
5170 
5171  difference(index) =
5172  integrate_difference_inner (exact_solution, norm, weight,
5173  update_flags, exponent,
5174  n_components, data);
5175  }
5176  else
5177  // the cell is a ghost cell or is artificial. write a zero into the
5178  // corresponding value of the returned vector
5179  difference(index) = 0;
5180  }
5181 
5182  } // namespace internal
5183 
5184 
5185 
5186 
5187  template <int dim, class InVector, class OutVector, int spacedim>
5188  void
5190  const DoFHandler<dim,spacedim> &dof,
5191  const InVector &fe_function,
5192  const Function<spacedim> &exact_solution,
5193  OutVector &difference,
5194  const Quadrature<dim> &q,
5195  const NormType &norm,
5196  const Function<spacedim> *weight,
5197  const double exponent)
5198  {
5199  internal
5200  ::do_integrate_difference (hp::MappingCollection<dim,spacedim>(mapping),
5201  dof, fe_function, exact_solution,
5202  difference, hp::QCollection<dim>(q),
5203  norm, weight, exponent);
5204  }
5205 
5206 
5207  template <int dim, class InVector, class OutVector, int spacedim>
5208  void
5210  const InVector &fe_function,
5211  const Function<spacedim> &exact_solution,
5212  OutVector &difference,
5213  const Quadrature<dim> &q,
5214  const NormType &norm,
5215  const Function<spacedim> *weight,
5216  const double exponent)
5217  {
5218  internal
5220  dof, fe_function, exact_solution,
5221  difference, hp::QCollection<dim>(q),
5222  norm, weight, exponent);
5223  }
5224 
5225 
5226 
5227  template <int dim, class InVector, class OutVector, int spacedim>
5228  void
5229  integrate_difference (const ::hp::MappingCollection<dim,spacedim> &mapping,
5230  const ::hp::DoFHandler<dim,spacedim> &dof,
5231  const InVector &fe_function,
5232  const Function<spacedim> &exact_solution,
5233  OutVector &difference,
5234  const ::hp::QCollection<dim> &q,
5235  const NormType &norm,
5236  const Function<spacedim> *weight,
5237  const double exponent)
5238  {
5239  internal
5240  ::do_integrate_difference (hp::MappingCollection<dim,spacedim>(mapping),
5241  dof, fe_function, exact_solution,
5242  difference, q,
5243  norm, weight, exponent);
5244  }
5245 
5246 
5247  template <int dim, class InVector, class OutVector, int spacedim>
5248  void
5249  integrate_difference (const ::hp::DoFHandler<dim,spacedim> &dof,
5250  const InVector &fe_function,
5251  const Function<spacedim> &exact_solution,
5252  OutVector &difference,
5253  const ::hp::QCollection<dim> &q,
5254  const NormType &norm,
5255  const Function<spacedim> *weight,
5256  const double exponent)
5257  {
5258  internal
5260  dof, fe_function, exact_solution,
5261  difference, q,
5262  norm, weight, exponent);
5263  }
5264 
5265 
5266 
5267  template <int dim, class InVector, int spacedim>
5268  void
5270  const InVector &fe_function,
5271  const Function<spacedim> &exact_function,
5272  Vector<double> &difference,
5273  const Point<spacedim> &point)
5274  {
5276  dof,
5277  fe_function,
5278  exact_function,
5279  difference,
5280  point);
5281  }
5282 
5283 
5284  template <int dim, class InVector, int spacedim>
5285  void
5287  const DoFHandler<dim,spacedim> &dof,
5288  const InVector &fe_function,
5289  const Function<spacedim> &exact_function,
5290  Vector<double> &difference,
5291  const Point<spacedim> &point)
5292  {
5293  const FiniteElement<dim> &fe = dof.get_fe();
5294 
5295  Assert(difference.size() == fe.n_components(),
5296  ExcDimensionMismatch(difference.size(), fe.n_components()));
5297 
5298  // first find the cell in which this point
5299  // is, initialize a quadrature rule with
5300  // it, and then a FEValues object
5301  const std::pair<typename DoFHandler<dim,spacedim>::active_cell_iterator, Point<spacedim> >
5302  cell_point = GridTools::find_active_cell_around_point (mapping, dof, point);
5303 
5304  AssertThrow(cell_point.first->is_locally_owned(),
5305  ExcPointNotAvailableHere());
5306  Assert(GeometryInfo<dim>::distance_to_unit_cell(cell_point.second) < 1e-10,
5307  ExcInternalError());
5308 
5309  const Quadrature<dim>
5310  quadrature (GeometryInfo<dim>::project_to_unit_cell(cell_point.second));
5311  FEValues<dim> fe_values(mapping, fe, quadrature, update_values);
5312  fe_values.reinit(cell_point.first);
5313 
5314  // then use this to get at the values of
5315  // the given fe_function at this point
5316  std::vector<Vector<double> > u_value(1, Vector<double> (fe.n_components()));
5317  fe_values.get_function_values(fe_function, u_value);
5318 
5319  if (fe.n_components() == 1)
5320  difference(0) = exact_function.value(point);
5321  else
5322  exact_function.vector_value(point, difference);
5323 
5324  for (unsigned int i=0; i<difference.size(); ++i)
5325  difference(i) -= u_value[0](i);
5326  }
5327 
5328 
5329  template <int dim, class InVector, int spacedim>
5330  void
5332  const InVector &fe_function,
5333  const Point<spacedim> &point,
5334  Vector<double> &value)
5335  {
5336 
5338  dof,
5339  fe_function,
5340  point,
5341  value);
5342  }
5343 
5344 
5345  template <int dim, class InVector, int spacedim>
5346  void
5348  const InVector &fe_function,
5349  const Point<spacedim> &point,
5350  Vector<double> &value)
5351  {
5353  dof,
5354  fe_function,
5355  point,
5356  value);
5357  }
5358 
5359 
5360  template <int dim, class InVector, int spacedim>
5361  double
5363  const InVector &fe_function,
5364  const Point<spacedim> &point)
5365  {
5367  dof,
5368  fe_function,
5369  point);
5370  }
5371 
5372 
5373  template <int dim, class InVector, int spacedim>
5374  double
5376  const InVector &fe_function,
5377  const Point<spacedim> &point)
5378  {
5380  dof,
5381  fe_function,
5382  point);
5383  }
5384 
5385 
5386  template <int dim, class InVector, int spacedim>
5387  void
5389  const DoFHandler<dim,spacedim> &dof,
5390  const InVector &fe_function,
5391  const Point<spacedim> &point,
5392  Vector<double> &value)
5393  {
5394  const FiniteElement<dim> &fe = dof.get_fe();
5395 
5396  Assert(value.size() == fe.n_components(),
5397  ExcDimensionMismatch(value.size(), fe.n_components()));
5398 
5399  // first find the cell in which this point
5400  // is, initialize a quadrature rule with
5401  // it, and then a FEValues object
5402  const std::pair<typename DoFHandler<dim,spacedim>::active_cell_iterator, Point<spacedim> >
5403  cell_point
5404  = GridTools::find_active_cell_around_point (mapping, dof, point);
5405 
5406  AssertThrow(cell_point.first->is_locally_owned(),
5407  ExcPointNotAvailableHere());
5408  Assert(GeometryInfo<dim>::distance_to_unit_cell(cell_point.second) < 1e-10,
5409  ExcInternalError());
5410 
5411  const Quadrature<dim>
5412  quadrature (GeometryInfo<dim>::project_to_unit_cell(cell_point.second));
5413 
5414  FEValues<dim> fe_values(mapping, fe, quadrature, update_values);
5415  fe_values.reinit(cell_point.first);
5416 
5417  // then use this to get at the values of
5418  // the given fe_function at this point
5419  std::vector<Vector<double> > u_value(1, Vector<double> (fe.n_components()));
5420  fe_values.get_function_values(fe_function, u_value);
5421 
5422  value = u_value[0];
5423  }
5424 
5425 
5426  template <int dim, class InVector, int spacedim>
5427  void
5429  const hp::DoFHandler<dim,spacedim> &dof,
5430  const InVector &fe_function,
5431  const Point<spacedim> &point,
5432  Vector<double> &value)
5433  {
5434  const hp::FECollection<dim, spacedim> &fe = dof.get_fe();
5435 
5436  Assert(value.size() == fe.n_components(),
5437  ExcDimensionMismatch(value.size(), fe.n_components()));
5438 
5439  // first find the cell in which this point
5440  // is, initialize a quadrature rule with
5441  // it, and then a FEValues object
5442  const std::pair<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator, Point<spacedim> >
5443  cell_point = GridTools::find_active_cell_around_point (mapping, dof, point);
5444 
5445  AssertThrow(cell_point.first->is_locally_owned(),
5446  ExcPointNotAvailableHere());
5447  Assert(GeometryInfo<dim>::distance_to_unit_cell(cell_point.second) < 1e-10,
5448  ExcInternalError());
5449 
5450  const Quadrature<dim>
5451  quadrature (GeometryInfo<dim>::project_to_unit_cell(cell_point.second));
5452  hp::FEValues<dim, spacedim> hp_fe_values(mapping, fe, hp::QCollection<dim>(quadrature), update_values);
5453  hp_fe_values.reinit(cell_point.first);
5454  const FEValues<dim, spacedim> &fe_values = hp_fe_values.get_present_fe_values();
5455 
5456  // then use this to get at the values of
5457  // the given fe_function at this point
5458  std::vector<Vector<double> > u_value(1, Vector<double> (fe.n_components()));
5459  fe_values.get_function_values(fe_function, u_value);
5460 
5461  value = u_value[0];
5462  }
5463 
5464 
5465  template <int dim, class InVector, int spacedim>
5466  double
5468  const DoFHandler<dim,spacedim> &dof,
5469  const InVector &fe_function,
5470  const Point<spacedim> &point)
5471  {
5472  const FiniteElement<dim> &fe = dof.get_fe();
5473 
5474  Assert(fe.n_components() == 1,
5475  ExcMessage ("Finite element is not scalar as is necessary for this function"));
5476 
5477  // first find the cell in which this point
5478  // is, initialize a quadrature rule with
5479  // it, and then a FEValues object
5480  const std::pair<typename DoFHandler<dim,spacedim>::active_cell_iterator, Point<spacedim> >
5481  cell_point = GridTools::find_active_cell_around_point (mapping, dof, point);
5482 
5483  AssertThrow(cell_point.first->is_locally_owned(),
5484  ExcPointNotAvailableHere());
5485  Assert(GeometryInfo<dim>::distance_to_unit_cell(cell_point.second) < 1e-10,
5486  ExcInternalError());
5487 
5488  const Quadrature<dim>
5489  quadrature (GeometryInfo<dim>::project_to_unit_cell(cell_point.second));
5490  FEValues<dim> fe_values(mapping, fe, quadrature, update_values);
5491  fe_values.reinit(cell_point.first);
5492 
5493  // then use this to get at the values of
5494  // the given fe_function at this point
5495  std::vector<double> u_value(1);
5496  fe_values.get_function_values(fe_function, u_value);
5497 
5498  return u_value[0];
5499  }
5500 
5501 
5502  template <int dim, class InVector, int spacedim>
5503  double
5505  const hp::DoFHandler<dim,spacedim> &dof,
5506  const InVector &fe_function,
5507  const Point<spacedim> &point)
5508  {
5509  const hp::FECollection<dim, spacedim> &fe = dof.get_fe();
5510 
5511  Assert(fe.n_components() == 1,
5512  ExcMessage ("Finite element is not scalar as is necessary for this function"));
5513 
5514  // first find the cell in which this point
5515  // is, initialize a quadrature rule with
5516  // it, and then a FEValues object
5517  const std::pair<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator, Point<spacedim> >
5518  cell_point = GridTools::find_active_cell_around_point (mapping, dof, point);
5519 
5520  AssertThrow(cell_point.first->is_locally_owned(),
5521  ExcPointNotAvailableHere());
5522  Assert(GeometryInfo<dim>::distance_to_unit_cell(cell_point.second) < 1e-10,
5523  ExcInternalError());
5524 
5525  const Quadrature<dim>
5526  quadrature (GeometryInfo<dim>::project_to_unit_cell(cell_point.second));
5527  hp::FEValues<dim, spacedim> hp_fe_values(mapping, fe, hp::QCollection<dim>(quadrature), update_values);
5528  hp_fe_values.reinit(cell_point.first);
5529  const FEValues<dim, spacedim> &fe_values = hp_fe_values.get_present_fe_values();
5530 
5531  // then use this to get at the values of
5532  // the given fe_function at this point
5533  std::vector<double> u_value(1);
5534  fe_values.get_function_values(fe_function, u_value);
5535 
5536  return u_value[0];
5537  }
5538 
5539 
5540 
5541  template <class VECTOR>
5542  void
5544  const std::vector<bool> &p_select)
5545  {
5546  if (p_select.size() == 0)
5547  {
5548  // In case of an empty boolean mask operate on the whole vector:
5549  v.add( - v.mean_value() );
5550  }
5551  else
5552  {
5553  // This function is not implemented for distributed vectors, so
5554  // if v is not a boring Vector or BlockVector:
5555  Assert( dynamic_cast<Vector<double> *>(& v)
5556  || dynamic_cast<Vector<float> *>(& v)
5557  || dynamic_cast<Vector<long double> *>(& v)
5558  || dynamic_cast<BlockVector<double> *>(& v)
5559  || dynamic_cast<BlockVector<float> *>(& v)
5560  || dynamic_cast<BlockVector<long double> *>(& v),
5561  ExcNotImplemented());
5562 
5563  const unsigned int n = v.size();
5564 
5565  Assert(p_select.size() == n,
5566  ExcDimensionMismatch(p_select.size(), n));
5567 
5568  typename VECTOR::value_type s = 0.;
5569  unsigned int counter = 0;
5570  for (unsigned int i=0; i<n; ++i)
5571  if (p_select[i])
5572  {
5573  s += v(i);
5574  ++counter;
5575  }
5576  // Error out if we have not constrained anything. Note that in this
5577  // case the vector v is always nonempty.
5578  Assert (n == 0 || counter > 0, ComponentMask::ExcNoComponentSelected());
5579 
5580  s /= counter;
5581 
5582  for (unsigned int i=0; i<n; ++i)
5583  if (p_select[i])
5584  v(i) -= s;
5585  }
5586  }
5587 
5588 
5589 
5590  template <int dim, class InVector, int spacedim>
5591  double
5593  const DoFHandler<dim,spacedim> &dof,
5594  const Quadrature<dim> &quadrature,
5595  const InVector &v,
5596  const unsigned int component)
5597  {
5598  Assert (v.size() == dof.n_dofs(),
5599  ExcDimensionMismatch (v.size(), dof.n_dofs()));
5600  Assert (component < dof.get_fe().n_components(),
5601  ExcIndexRange(component, 0, dof.get_fe().n_components()));
5602 
5603  FEValues<dim,spacedim> fe(mapping, dof.get_fe(), quadrature,
5605  | update_values));
5606 
5607  typename DoFHandler<dim,spacedim>::active_cell_iterator cell;
5608  std::vector<Vector<double> > values(quadrature.size(),
5609  Vector<double> (dof.get_fe().n_components()));
5610 
5611  double mean = 0.;
5612  double area = 0.;
5613  // Compute mean value
5614  for (cell = dof.begin_active(); cell != dof.end(); ++cell)
5615  if (cell->is_locally_owned())
5616  {
5617  fe.reinit (cell);
5618  fe.get_function_values(v, values);
5619  for (unsigned int k=0; k< quadrature.size(); ++k)
5620  {
5621  mean += fe.JxW(k) * values[k](component);
5622  area += fe.JxW(k);
5623  }
5624  }
5625 
5626 #ifdef DEAL_II_WITH_P4EST
5627  // if this was a distributed DoFHandler, we need to do the reduction
5628  // over the entire domain
5630  p_d_triangulation
5631  = dynamic_cast<const parallel::distributed::Triangulation<dim,spacedim> *>(&dof.get_tria()))
5632  {
5633  double my_values[2] = { mean, area };
5634  double global_values[2];
5635 
5636  MPI_Allreduce (&my_values, &global_values, 2, MPI_DOUBLE,
5637  MPI_SUM,
5638  p_d_triangulation->get_communicator());
5639 
5640  mean = global_values[0];
5641  area = global_values[1];
5642  }
5643 #endif
5644 
5645  return (mean/area);
5646  }
5647 
5648 
5649  template <int dim, class InVector, int spacedim>
5650  double
5652  const Quadrature<dim> &quadrature,
5653  const InVector &v,
5654  const unsigned int component)
5655  {
5656  return compute_mean_value(StaticMappingQ1<dim,spacedim>::mapping, dof, quadrature, v, component);
5657  }
5658 }
5659 
5660 DEAL_II_NAMESPACE_CLOSE
5661 
5662 #endif
Transformed quadrature weights.
void reinit(const TriaIterator< DoFCellAccessor< DH, level_dof_access > > cell, const unsigned int face_no)
double shape_value_component(const unsigned int function_no, const unsigned int point_no, const unsigned int component) const
Shape function values.
double get_inhomogeneity(const size_type line) const
static const unsigned int invalid_unsigned_int
Definition: types.h:191
void create_boundary_right_hand_side(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, const Quadrature< dim-1 > &quadrature, const Function< spacedim > &rhs_function, Vector< double > &rhs_vector, const std::set< types::boundary_id > &boundary_indicators)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:858
types::global_dof_index n_dofs() const
active_cell_iterator begin_active(const unsigned int level=0) const
unsigned int max_dofs_per_cell() const
bool have_same_coarse_mesh(const Triangulation< dim, spacedim > &mesh_1, const Triangulation< dim, spacedim > &mesh_2)
bool can_store_line(const size_type line_index) const
bool is_constrained(const size_type index) 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)
void map_dof_to_boundary_indices(const DH &dof_handler, std::vector< types::global_dof_index > &mapping)
void interpolate(const Mapping< dim, spacedim > &mapping, const DH< dim, spacedim > &dof, const Function< spacedim > &function, VECTOR &vec)
unsigned int n_selected_components(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
size_type m() const
const unsigned int dofs_per_cell
Definition: fe_values.h:1418
const Point< spacedim > & normal_vector(const unsigned int i) const
std::map< types::boundary_id, const Function< dim > * > type
Definition: function_map.h:56
const std::vector< Point< dim-1 > > & get_unit_face_support_points() const
Volume element.
const std::vector< std::pair< size_type, double > > * get_constraint_entries(const size_type line) const
::ExceptionBase & ExcMessage(std::string arg1)
size_type n() const
#define AssertIndexRange(index, range)
Definition: exceptions.h:888
cell_iterator end() const
const unsigned int degree
Definition: fe_base.h:287
void initialize(const MATRIX &m, bool expect_constrained_source=false)
const IndexSet & get_local_lines() const
double compute_mean_value(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim > &quadrature, const InVector &v, const unsigned int component)
void create_right_hand_side(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, const Quadrature< dim > &quadrature, const Function< spacedim > &rhs_function, Vector< double > &rhs_vector)
unsigned int size() const
Transformed quadrature points.
#define AssertThrow(cond, exc)
Definition: exceptions.h:362
void integrate_difference(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const InVector &fe_function, const Function< spacedim > &exact_solution, OutVector &difference, const Quadrature< dim > &q, const NormType &norm, const Function< spacedim > *weight, const double exponent)
void vmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
virtual const FiniteElement< dim, spacedim > & base_element(const unsigned int index) const
const unsigned int dofs_per_line
Definition: fe_base.h:220
bool is_primitive(const unsigned int i) const
Definition: fe.h:2237
const std::vector< double > & get_JxW_values() const
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const
const Point< spacedim > & quadrature_point(const unsigned int i) const
virtual unsigned int face_to_cell_index(const unsigned int face_dof_index, const unsigned int face, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false) const
bool is_finite(const double x)
virtual void vector_value(const Point< dim > &p, Vector< double > &values) const
const Quadrature< dim > & get_quadrature() const
bool represents_n_components(const unsigned int n) const
const hp::FECollection< dim, spacedim > & get_fe() const
unsigned int element_multiplicity(const unsigned int index) const
Definition: fe.h:2072
void add_line(const size_type line)
unsigned int size() const
void get_function_values(const InputVector &fe_function, std::vector< number > &values) const
const std::vector< DerivativeForm< 1, dim, spacedim > > & get_jacobians() const
const std::vector< Point< spacedim > > & get_quadrature_points() const
size_type n() const
void compute_no_normal_flux_constraints(const DH< dim, spacedim > &dof_handler, const unsigned int first_vector_component, const std::set< types::boundary_id > &boundary_ids, ConstraintMatrix &constraints, const Mapping< dim, spacedim > &mapping)
virtual Point< dim-1 > unit_face_support_point(const unsigned int index) const
void push_back(const FiniteElement< dim, spacedim > &new_fe)
void add_entry(const size_type line, const size_type column, const double value)
number diag_element(const size_type i) const
void point_value(const DoFHandler< dim, spacedim > &dof, const InVector &fe_function, const Point< spacedim > &point, Vector< double > &value)
unsigned int global_dof_index
Definition: types.h:100
virtual void reinit(const size_type N, const bool fast=false)
#define Assert(cond, exc)
Definition: exceptions.h:299
types::global_dof_index n_dofs() const
UpdateFlags
void condense(const SparsityPattern &uncondensed, SparsityPattern &condensed) const
const ::FEValues< dim, spacedim > & get_present_fe_values() const
void project_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const typename FunctionMap< spacedim >::type &boundary_functions, const Quadrature< dim-1 > &q, std::map< types::global_dof_index, double > &boundary_values, std::vector< unsigned int > component_mapping)
void reinit(const TriaIterator< DoFCellAccessor< DH, level_dof_access > > cell)
virtual void gradient_list(const std::vector< Point< dim > > &points, std::vector< Tensor< 1, dim > > &gradients, const unsigned int component=0) const
std::size_t size() const
const ComponentMask & get_nonzero_components(const unsigned int i) const
Definition: fe.h:2215
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
virtual void vector_value_list(const std::vector< Point< dim > > &points, std::vector< Vector< double > > &values) const
void project_boundary_values_curl_conforming(const DoFHandler< dim > &dof_handler, const unsigned int first_vector_component, const Function< dim > &boundary_function, const types::boundary_id boundary_component, ConstraintMatrix &constraints, const Mapping< dim > &mapping)
void interpolate_boundary_values(const Mapping< DH::dimension, DH::space_dimension > &mapping, const DH &dof, const typename FunctionMap< DH::space_dimension >::type &function_map, std::map< types::global_dof_index, double > &boundary_values, const ComponentMask &component_mask_)
void make_boundary_sparsity_pattern(const DH &dof, const std::vector< types::global_dof_index > &dof_to_boundary_mapping, SparsityPattern &sparsity_pattern)
void project_boundary_values_div_conforming(const DoFHandler< dim > &dof_handler, const unsigned int first_vector_component, const Function< dim > &boundary_function, const types::boundary_id boundary_component, ConstraintMatrix &constraints, const Mapping< dim > &mapping)
virtual void vector_gradient_list(const std::vector< Point< dim > > &points, std::vector< std::vector< Tensor< 1, dim > > > &gradients) const
unsigned int n_base_elements() const
Definition: fe.h:2062
void create_point_source_vector(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, const Point< spacedim > &p, Vector< double > &rhs_vector)
unsigned int n_components() const
size_type m() const
std::pair< unsigned int, unsigned int > face_system_to_component_index(const unsigned int index) const
Definition: fe.h:2098
void apply_constraints(VECTOR &v, const bool matrix_is_symmetric) const
std::pair< std::pair< unsigned int, unsigned int >, unsigned int > system_to_base_index(const unsigned int index) const
Definition: fe.h:2127
iterator begin()
const unsigned int n_components
Definition: function.h:130
const Mapping< dim, spacedim > & get_mapping() const
void project(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const ConstraintMatrix &constraints, const Quadrature< dim > &quadrature, const Function< spacedim > &function, VECTOR &vec_result, const bool enforce_zero_boundary, const Quadrature< dim-1 > &q_boundary, const bool project_to_boundary_first)
unsigned int first_selected_component(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
const unsigned int dofs_per_cell
Definition: fe_base.h:271
const unsigned int n_quadrature_points
Definition: fe_values.h:1411
void point_difference(const DoFHandler< dim, spacedim > &dof, const InVector &fe_function, const Function< spacedim > &exact_function, Vector< double > &difference, const Point< spacedim > &point)
const double & shape_value(const unsigned int function_no, const unsigned int point_no) const
void create_boundary_mass_matrix(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim-1 > &q, SparseMatrix< double > &matrix, const typename FunctionMap< spacedim >::type &boundary_functions, Vector< double > &rhs_vector, std::vector< types::global_dof_index > &dof_to_boundary_mapping, const Function< spacedim > *const weight=0, std::vector< unsigned int > component_mapping=std::vector< unsigned int >())
void mTmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
std::ostream & operator<<(std::ostream &os, const Vector< number > &v)
Definition: vector.h:1546
void copy_from(const size_type n_rows, const size_type n_cols, const ForwardIterator begin, const ForwardIterator end, const bool optimize_diagonal) DEAL_II_DEPRECATED
::ExceptionBase & ExcNumberNotFinite()
void make_mapping(const GridClass &source_grid, const GridClass &destination_grid)
unsigned int n_components() const
void subtract_mean_value(VECTOR &v, const std::vector< bool > &p_select)
active_cell_iterator begin_active(const unsigned int level=0) const
std::pair< unsigned int, unsigned int > system_to_component_index(const unsigned int index) const
Definition: fe.h:2048
Shape function gradients.
Normal vectors.
size_type n_rows() const
void solve(const MATRIX &A, VECTOR &x, const VECTOR &b, const PRECONDITIONER &precondition)
const unsigned int dofs_per_face
Definition: fe_base.h:264
void distribute(const VectorType &condensed, VectorType &uncondensed) const
void make_sparsity_pattern(const DH &dof, SparsityPattern &sparsity_pattern, const ConstraintMatrix &constraints=ConstraintMatrix(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
void initialize(const MATRIX &A, const typename BaseClass::AdditionalData &parameters=typename BaseClass::AdditionalData())
::ExceptionBase & ExcNotImplemented()
void add_constraint(const size_type i, const double v)
void create_mass_matrix(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim > &q, SparseMatrix< number > &matrix, const Function< spacedim > *const a=0, const ConstraintMatrix &constraints=ConstraintMatrix())
const unsigned int dofs_per_vertex
Definition: fe_base.h:214
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
virtual double value(const Point< dim > &p, const unsigned int component=0) const
unsigned char boundary_id
Definition: types.h:128
const types::boundary_id internal_face_boundary_id
Definition: types.h:253
void set_inhomogeneity(const size_type line, const double value)
const types::global_dof_index invalid_dof_index
Definition: types.h:217
const Triangulation< dim, spacedim > & get_tria() const
Container< dim, spacedim >::active_cell_iterator find_active_cell_around_point(const Container< dim, spacedim > &container, const Point< spacedim > &p)
::ExceptionBase & ExcInternalError()
void reinit(const TriaIterator< DoFCellAccessor< DH, lda > > cell, const unsigned int face_no, 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)
bool has_face_support_points() const
void compute_normal_flux_constraints(const DH< dim, spacedim > &dof_handler, const unsigned int first_vector_component, const std::set< types::boundary_id > &boundary_ids, ConstraintMatrix &constraints, const Mapping< dim, spacedim > &mapping)
void push_back(const Quadrature< dim > &new_quadrature)
Definition: q_collection.h:251
void apply_boundary_values(const std::map< types::global_dof_index, double > &boundary_values, SparseMatrix< number > &matrix, Vector< number > &solution, Vector< number > &right_hand_side, const bool eliminate_columns=true)
double JxW(const unsigned int quadrature_point) const
const std::vector< Point< spacedim > > & get_normal_vectors() const
cell_iterator end() const
real_type norm() const
real_type norm() const
const FiniteElement< dim, spacedim > & get_fe() const