18 #ifndef _deal2__vector_tools_templates_h
19 #define _deal2__vector_tools_templates_h
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>
63 #include <deal.II/base/std_cxx1x/array.h>
77 template <
class VECTOR,
int dim,
int spacedim,
template <
int,
int>
class DH>
79 const DH<dim,spacedim> &dof,
83 Assert (dof.get_fe().n_components() ==
function.n_components,
85 function.n_components));
89 const bool fe_is_system = (n_components != 1);
91 typename DH<dim,spacedim>::active_cell_iterator cell = dof.begin_active(),
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)
112 unit_support_points[fe_index] = fe[fe_index].get_unit_support_points();
113 Assert (unit_support_points[fe_index].size() != 0,
114 ExcNonInterpolatingFE());
135 std::vector<std::vector<types::global_dof_index> > dofs_of_rep_points(fe.size());
138 std::vector<std::vector<types::global_dof_index> > dof_to_rep_index_table(fe.size());
140 std::vector<unsigned int> n_rep_points (fe.size(), 0);
142 for (
unsigned int fe_index=0; fe_index<fe.size(); ++fe_index)
144 for (
unsigned int i=0; i<fe[fe_index].dofs_per_cell; ++i)
146 bool representative=
true;
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]])
157 dof_to_rep_index_table[fe_index].push_back(j-1);
158 representative=
false;
165 dof_to_rep_index_table[fe_index].push_back(dofs_of_rep_points[fe_index].size());
167 dofs_of_rep_points[fe_index].push_back(i);
168 ++n_rep_points[fe_index];
172 Assert(dofs_of_rep_points[fe_index].size()==n_rep_points[fe_index],
174 Assert(dof_to_rep_index_table[fe_index].size()==fe[fe_index].dofs_per_cell,
178 const unsigned int max_rep_points = *std::max_element (n_rep_points.begin(),
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);
189 std::vector<std::vector<double> > function_values_scalar(fe.size());
190 std::vector<std::vector<Vector<double> > > function_values_system(fe.size());
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]));
205 for (; cell!=endc; ++cell)
206 if (cell->is_locally_owned())
208 const unsigned int fe_index = cell->active_fe_index();
214 const std::vector<Point<spacedim> > &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]];
224 dofs_on_cell.resize (fe[fe_index].dofs_per_cell);
225 cell->get_dof_indices (dofs_on_cell);
233 function_values_system[fe_index]
234 .resize (n_rep_points[fe_index],
236 function.vector_value_list (rep_points,
237 function_values_system[fe_index]);
241 for (
unsigned int i=0; i<fe[fe_index].dofs_per_cell; ++i)
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];
247 = function_values_system[fe_index][rep_dof](component);
255 function_values_scalar[fe_index].resize (n_rep_points[fe_index]);
256 function.value_list (rep_points,
257 function_values_scalar[fe_index],
262 for (
unsigned int i=0; i<fe[fe_index].dofs_per_cell; ++i)
264 = function_values_scalar[fe_index][dof_to_rep_index_table[fe_index][i]];
267 vec.compress(VectorOperation::insert);
271 template <
class VECTOR,
class DH>
283 template <
int dim,
class InVector,
class OutVector,
int spacedim>
288 const InVector &data_1,
294 std::vector<short unsigned int> touch_count (dof_2.
n_dofs(), 0);
295 std::vector<types::global_dof_index> local_dof_indices (dof_2.
get_fe().dofs_per_cell);
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();
301 for (; h != endh; ++h, ++l)
303 h->get_dof_values(data_1, cell_data_1);
304 transfer.
vmult(cell_data_2, cell_data_1);
306 l->get_dof_indices (local_dof_indices);
309 for (
unsigned int j=0; j<dof_2.
get_fe().dofs_per_cell; ++j)
311 data_2(local_dof_indices[j]) += cell_data_2(j);
317 ++touch_count[local_dof_indices[j]];
324 for (
unsigned int i=0; i<dof_2.
n_dofs(); ++i)
326 Assert (touch_count[i] != 0,
329 data_2(i) /= touch_count[i];
343 interpolate_zero_boundary_values (
const DH &dof_handler,
344 std::map<types::global_dof_index,double> &boundary_values)
346 const unsigned int dim = DH::dimension;
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))
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)
387 boundary_values[face_dof_indices[i]] = 0.;
394 template <
int dim,
int spacedim,
395 template <
int,
int>
class DH,
398 interpolate_to_different_mesh (
const DH<dim, spacedim> &dof1,
400 const DH<dim, spacedim> &dof2,
404 ExcMessage (
"The two containers must represent triangulations that "
405 "have the same coarse meshes"));
413 interpolate_to_different_mesh(intergridmap, u1, dummy, u2);
418 template <
int dim,
int spacedim,
419 template <
int,
int>
class DH,
422 interpolate_to_different_mesh (
const DH<dim, spacedim> &dof1,
424 const DH<dim, spacedim> &dof2,
429 ExcMessage (
"The two containers must represent triangulations that "
430 "have the same coarse meshes"));
435 interpolate_to_different_mesh(intergridmap, u1, constraints, u2);
440 template <
int dim,
int spacedim,
441 template <
int,
int>
class DH,
444 interpolate_to_different_mesh (
const InterGridMap<DH<dim, spacedim> > &intergridmap,
468 typename DH<dim,spacedim>::cell_iterator cell1 = dof1.
begin();
469 const typename DH<dim,spacedim>::cell_iterator endc1 = dof1.end();
471 for (; cell1 != endc1; ++cell1)
473 const typename DH<dim,spacedim>::cell_iterator cell2 = intergridmap[cell1];
477 if (cell1->level() != cell2->level())
480 if (!cell1->active() && !cell2->active())
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"));
487 cache.
reinit(cell1->get_fe().dofs_per_cell);
491 cell1->get_interpolated_dof_values(u1, cache);
492 cell2->set_dof_values_by_interpolation(cache, u2);
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,
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)
517 if (enforce_zero_boundary ==
true)
523 interpolate_zero_boundary_values (dof, boundary_values);
527 if (project_to_boundary_first ==
true)
536 const std::vector<types::boundary_id>
537 used_boundary_indicators = dof.get_tria().get_boundary_indicators();
540 for (
unsigned int i=0; i<used_boundary_indicators.size(); ++i)
541 boundary_functions[used_boundary_indicators[i]] = &
function;
552 bool constraints_and_b_v_are_compatible (
const ConstraintMatrix &constraints,
553 std::map<types::global_dof_index,double> &boundary_values)
555 for (std::map<types::global_dof_index,double>::iterator it=boundary_values.begin();
556 it != boundary_values.end(); ++it)
571 template <
int dim,
int spacedim,
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,
579 const Q_or_QC<dim> &quadrature,
582 const bool enforce_zero_boundary,
583 const Q_or_QC<dim-1> &q_boundary,
584 const bool project_to_boundary_first)
586 Assert (dof.get_fe().n_components() ==
function.n_components,
588 function.n_components));
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);
598 const bool constraints_are_compatible =
599 constraints_and_b_v_are_compatible(constraints, boundary_values);
607 !constraints_are_compatible);
609 sparsity.copy_from (csp);
619 if (constraints_are_compatible)
623 mass_matrix,
function, tmp,
625 if (boundary_values.size() > 0)
627 mass_matrix, vec, tmp,
634 mass_matrix,
function, tmp);
636 mass_matrix, vec, tmp,
638 constraints.
condense(mass_matrix, tmp);
651 cg.solve (mass_matrix, vec, tmp, prec);
658 for (
unsigned int i=0; i<vec.size(); ++i)
659 vec_result(i) = vec(i);
664 template <
int dim,
class VECTOR,
int spacedim>
671 const bool enforce_zero_boundary,
673 const bool project_to_boundary_first)
675 do_project (mapping, dof, constraints, quadrature,
676 function, vec_result,
677 enforce_zero_boundary, q_boundary,
678 project_to_boundary_first);
682 template <
int dim,
class VECTOR,
int spacedim>
688 const bool enforce_zero_boundary,
690 const bool project_to_boundary_first)
693 enforce_zero_boundary, q_boundary, project_to_boundary_first);
698 template <
int dim,
class VECTOR,
int spacedim>
705 const bool enforce_zero_boundary,
707 const bool project_to_boundary_first)
709 do_project (mapping, dof, constraints, quadrature,
710 function, vec_result,
711 enforce_zero_boundary, q_boundary,
712 project_to_boundary_first);
716 template <
int dim,
class VECTOR,
int spacedim>
722 const bool enforce_zero_boundary,
724 const bool project_to_boundary_first)
727 dof, constraints, quadrature,
function, vec,
728 enforce_zero_boundary, q_boundary, project_to_boundary_first);
732 template <
int dim,
int spacedim>
755 std::vector<types::global_dof_index> dofs (dofs_per_cell);
758 typename DoFHandler<dim,spacedim>::active_cell_iterator
760 endc = dof_handler.
end();
764 std::vector<double> rhs_values(n_q_points);
766 for (; cell!=endc; ++cell)
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] *
781 cell->get_dof_indices (dofs);
783 for (
unsigned int i=0; i<dofs_per_cell; ++i)
784 rhs_vector(dofs[i]) += cell_vector(i);
790 std::vector<Vector<double> > rhs_values(n_q_points,
793 for (; cell!=endc; ++cell)
806 for (
unsigned int point=0; point<n_q_points; ++point)
807 for (
unsigned int i=0; i<dofs_per_cell; ++i)
809 const unsigned int component
812 cell_vector(i) += rhs_values[point](component) *
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)
827 cell_vector(i) += rhs_values[point](comp_i) *
833 cell->get_dof_indices (dofs);
835 for (
unsigned int i=0; i<dofs_per_cell; ++i)
836 rhs_vector(dofs[i]) += cell_vector(i);
843 template <
int dim,
int spacedim>
850 rhs_function, rhs_vector);
856 template <
int dim,
int spacedim>
880 typename hp::DoFHandler<dim,spacedim>::active_cell_iterator
882 endc = dof_handler.
end();
886 std::vector<double> rhs_values;
888 for (; cell!=endc; ++cell)
896 rhs_values.resize (n_q_points);
897 dofs.resize (dofs_per_cell);
898 cell_vector.reinit (dofs_per_cell);
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] *
911 cell->get_dof_indices (dofs);
913 for (
unsigned int i=0; i<dofs_per_cell; ++i)
914 rhs_vector(dofs[i]) += cell_vector(i);
920 std::vector<Vector<double> > rhs_values;
922 for (; cell!=endc; ++cell)
930 rhs_values.resize (n_q_points,
932 dofs.resize (dofs_per_cell);
933 cell_vector.reinit (dofs_per_cell);
943 if (cell->get_fe().is_primitive ())
945 for (
unsigned int point=0; point<n_q_points; ++point)
946 for (
unsigned int i=0; i<dofs_per_cell; ++i)
948 const unsigned int component
949 = cell->get_fe().system_to_component_index(i).first;
951 cell_vector(i) += rhs_values[point](component) *
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])
965 cell_vector(i) += rhs_values[point](comp_i) *
971 cell->get_dof_indices (dofs);
973 for (
unsigned int i=0; i<dofs_per_cell; ++i)
974 rhs_vector(dofs[i]) += cell_vector(i);
981 template <
int dim,
int spacedim>
988 dof_handler, quadrature,
989 rhs_function, rhs_vector);
995 template <
int dim,
int spacedim>
1004 ExcMessage (
"This function only works for scalar finite elements"));
1008 std::pair<typename DoFHandler<dim,spacedim>::active_cell_iterator,
Point<spacedim> >
1016 fe_values.
reinit(cell_point.first);
1018 const unsigned int dofs_per_cell = dof_handler.
get_fe().dofs_per_cell;
1020 std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1021 cell_point.first->get_dof_indices (local_dof_indices);
1023 for (
unsigned int i=0; i<dofs_per_cell; i++)
1024 rhs_vector(local_dof_indices[i]) = fe_values.shape_value(i,0);
1029 template <
int dim,
int spacedim>
1039 template <
int dim,
int spacedim>
1048 ExcMessage (
"This function only works for scalar finite elements"));
1052 std::pair<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator,
Point<spacedim> >
1058 FEValues<dim> fe_values(mapping[cell_point.first->active_fe_index()],
1060 fe_values.
reinit(cell_point.first);
1062 const unsigned int dofs_per_cell = cell_point.first->get_fe().dofs_per_cell;
1064 std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1065 cell_point.first->get_dof_indices (local_dof_indices);
1067 for (
unsigned int i=0; i<dofs_per_cell; i++)
1068 rhs_vector(local_dof_indices[i]) = fe_values.shape_value(i,0);
1073 template <
int dim,
int spacedim>
1086 template <
int dim,
int spacedim>
1096 ExcMessage (
"This function only works for vector-valued finite elements."));
1100 std::pair<typename DoFHandler<dim,spacedim>::active_cell_iterator,
Point<spacedim> >
1109 fe_values.
reinit(cell_point.first);
1111 const unsigned int dofs_per_cell = dof_handler.
get_fe().dofs_per_cell;
1113 std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1114 cell_point.first->get_dof_indices (local_dof_indices);
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);
1122 template <
int dim,
int spacedim>
1129 p, orientation, rhs_vector);
1133 template <
int dim,
int spacedim>
1143 ExcMessage (
"This function only works for vector-valued finite elements."));
1147 std::pair<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator,
Point<spacedim> >
1154 FEValues<dim> fe_values(mapping[cell_point.first->active_fe_index()],
1156 fe_values.
reinit(cell_point.first);
1158 const unsigned int dofs_per_cell = cell_point.first->get_fe().dofs_per_cell;
1160 std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
1161 cell_point.first->get_dof_indices (local_dof_indices);
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);
1169 template <
int dim,
int spacedim>
1177 p, orientation, rhs_vector);
1182 template <
int dim,
int spacedim>
1189 const std::set<types::boundary_id> &boundary_indicators)
1208 std::vector<types::global_dof_index> dofs (dofs_per_cell);
1211 typename DoFHandler<dim,spacedim>::active_cell_iterator cell = dof_handler.
begin_active(),
1212 endc = dof_handler.
end();
1214 if (n_components==1)
1216 std::vector<double> rhs_values(n_q_points);
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())
1224 boundary_indicators.end())))
1226 fe_values.
reinit(cell, face);
1228 const std::vector<double> &weights = fe_values.
get_JxW_values ();
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] *
1238 cell->get_dof_indices (dofs);
1240 for (
unsigned int i=0; i<dofs_per_cell; ++i)
1241 rhs_vector(dofs[i]) += cell_vector(i);
1246 std::vector<Vector<double> > rhs_values(n_q_points,
Vector<double>(n_components));
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())
1254 boundary_indicators.end())))
1256 fe_values.
reinit(cell, face);
1258 const std::vector<double> &weights = fe_values.
get_JxW_values ();
1267 for (
unsigned int point=0; point<n_q_points; ++point)
1268 for (
unsigned int i=0; i<dofs_per_cell; ++i)
1270 const unsigned int component
1273 cell_vector(i) += rhs_values[point](component) *
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)
1289 += rhs_values[point](comp_i) *
1295 cell->get_dof_indices (dofs);
1297 for (
unsigned int i=0; i<dofs_per_cell; ++i)
1298 rhs_vector(dofs[i]) += cell_vector(i);
1305 template <
int dim,
int spacedim>
1311 const std::set<types::boundary_id> &boundary_indicators)
1315 rhs_function, rhs_vector,
1316 boundary_indicators);
1321 template <
int dim,
int spacedim>
1328 const std::set<types::boundary_id> &boundary_indicators)
1348 typename hp::DoFHandler<dim,spacedim>::active_cell_iterator
1350 endc = dof_handler.
end();
1352 if (n_components==1)
1354 std::vector<double> rhs_values;
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())
1362 boundary_indicators.end())))
1364 x_fe_values.
reinit(cell, face);
1370 rhs_values.resize (n_q_points);
1372 const std::vector<double> &weights = fe_values.
get_JxW_values ();
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] *
1382 dofs.resize(dofs_per_cell);
1383 cell->get_dof_indices (dofs);
1385 for (
unsigned int i=0; i<dofs_per_cell; ++i)
1386 rhs_vector(dofs[i]) += cell_vector(i);
1391 std::vector<Vector<double> > rhs_values;
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())
1399 boundary_indicators.end())))
1401 x_fe_values.
reinit(cell, face);
1409 const std::vector<double> &weights = fe_values.
get_JxW_values ();
1416 if (cell->get_fe().is_primitive ())
1418 for (
unsigned int point=0; point<n_q_points; ++point)
1419 for (
unsigned int i=0; i<dofs_per_cell; ++i)
1421 const unsigned int component
1422 = cell->get_fe().system_to_component_index(i).first;
1424 cell_vector(i) += rhs_values[point](component) *
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])
1440 += rhs_values[point](comp_i) *
1445 dofs.resize(dofs_per_cell);
1446 cell->get_dof_indices (dofs);
1448 for (
unsigned int i=0; i<dofs_per_cell; ++i)
1449 rhs_vector(dofs[i]) += cell_vector(i);
1456 template <
int dim,
int spacedim>
1462 const std::set<types::boundary_id> &boundary_indicators)
1465 dof_handler, quadrature,
1466 rhs_function, rhs_vector,
1467 boundary_indicators);
1484 template <
int,
int>
class M_or_MC>
1486 void do_interpolate_boundary_values (
const M_or_MC<DH::dimension, DH::space_dimension> &,
1489 std::map<types::global_dof_index,double> &boundary_values,
1491 const ::internal::int2type<1>)
1493 const unsigned int dim = DH::dimension;
1494 const unsigned int spacedim=DH::space_dimension;
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 "
1504 if (function_map.size() == 0)
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)
1513 (function_map.find(cell->face(direction)->boundary_indicator()) != function_map.end()))
1516 = *function_map.find(cell->face(direction)->boundary_indicator())->second;
1526 ComponentMask::ExcNoComponentSelected());
1546 = boundary_function.
value (cell->vertex(direction));
1548 boundary_function.
vector_value (cell->vertex(direction),
1553 boundary_values[cell->
1554 vertex_dof_index(direction,i,
1555 cell->active_fe_index())]
1563 template <
int,
int>
class M_or_MC>
1566 do_interpolate_boundary_values (
const M_or_MC<DH::dimension, DH::space_dimension> &mapping,
1569 std::map<types::global_dof_index,double> &boundary_values,
1571 const ::internal::int2type<DH::dimension>)
1573 const unsigned int dim = DH::dimension;
1574 const unsigned int spacedim=DH::space_dimension;
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 "
1584 if (function_map.size() == 0)
1588 ExcInvalidBoundaryIndicator());
1590 const unsigned int n_components = DoFTools::n_components(dof);
1591 const bool fe_is_system = (n_components != 1);
1594 i!=function_map.end(); ++i)
1595 Assert (n_components == i->second->n_components,
1599 std::vector<types::global_dof_index> face_dofs;
1600 face_dofs.reserve (DoFTools::max_dofs_per_face(dof));
1602 std::vector<Point<spacedim> > dof_locations;
1603 dof_locations.reserve (DoFTools::max_dofs_per_face(dof));
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));
1618 for (
unsigned int f=0; f<finite_elements.size(); ++f)
1665 typename DH::active_cell_iterator cell = dof.begin_active(),
1667 for (; cell!=endc; ++cell)
1668 if (!cell->is_artificial())
1669 for (
unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
1680 for (
unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
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)
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 "
1694 const typename DH::face_iterator face = cell->face(face_no);
1700 if ((function_map.find(boundary_component) != function_map.end())
1702 (cell->get_fe().dofs_per_face > 0))
1705 x_fe_values.reinit(cell, face_no);
1706 const ::FEFaceValues<dim,spacedim> &fe_values =
1707 x_fe_values.get_present_fe_values();
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 ();
1726 function_map.find(boundary_component)->second
1727 ->vector_value_list (dof_locations, dof_values_system);
1733 for (
unsigned int i=0; i<face_dofs.size(); ++i)
1735 unsigned int component;
1744 const unsigned int cell_i
1767 for (
unsigned int c=0; c<n_components; ++c)
1769 Assert (component_mask[c] ==
false,
1770 FETools::ExcFENotPrimitive());
1780 if (component_mask[component] ==
true)
1781 boundary_values[face_dofs[i]] = dof_values_system[i](component);
1789 function_map.find(boundary_component)->second
1790 ->value_list (dof_locations, dof_values_scalar, 0);
1794 for (
unsigned int i=0; i<face_dofs.size(); ++i)
1795 boundary_values[face_dofs[i]] = dof_values_scalar[i];
1810 std::map<types::global_dof_index,double> &boundary_values,
1813 do_interpolate_boundary_values (mapping, dof, function_map, boundary_values,
1826 std::map<types::global_dof_index,double> &boundary_values,
1830 function_map[boundary_component] = &boundary_function;
1836 template <
int dim,
int spacedim>
1841 std::map<types::global_dof_index,double> &boundary_values,
1844 do_interpolate_boundary_values (mapping, dof, function_map, boundary_values,
1856 std::map<types::global_dof_index,double> &boundary_values,
1860 dof, boundary_component,
1861 boundary_function, boundary_values, component_mask);
1870 std::map<types::global_dof_index,double> &boundary_values,
1875 boundary_values, component_mask);
1894 std::map<types::global_dof_index,double> boundary_values;
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)
1905 constraints.
add_line (boundary_value->first);
1907 boundary_value->second);
1925 function_map[boundary_component] = &boundary_function;
1942 dof, boundary_component,
1943 boundary_function, constraints, component_mask);
1958 constraints, component_mask);
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>
1974 do_project_boundary_values (
const M_or_MC<dim, spacedim> &mapping,
1975 const DH<dim, spacedim> &dof,
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)
1996 if (component_mapping.size() == 0)
1998 AssertDimension (dof.get_fe().n_components(), boundary_functions.begin()->second->n_components);
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;
2007 AssertDimension (dof.get_fe().n_components(), component_mapping.size());
2009 std::vector<types::global_dof_index> dof_to_boundary_mapping;
2010 std::set<types::boundary_id> selected_boundary_components;
2012 i!=boundary_functions.end(); ++i)
2013 selected_boundary_components.insert (i->first);
2016 dof_to_boundary_mapping);
2019 if (dof.n_boundary_dofs (boundary_functions) == 0)
2024 dof.n_boundary_dofs (boundary_functions));
2027 dof_to_boundary_mapping,
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)
2057 if (cell->at_boundary(f))
2060 level = cell->level();
2078 mass_matrix, boundary_functions,
2096 std::vector<bool> excluded_dofs(
mass_matrix.m(),
false);
2098 double max_element = 0.;
2104 if (
mass_matrix.diag_element(i) < 1.e-8 * max_element)
2106 filtered_mass_matrix.add_constraint(i, 0.);
2109 excluded_dofs[i] =
true;
2126 filtered_precondition.
initialize(prec,
true);
2128 cg.solve (filtered_mass_matrix, boundary_projection, rhs, filtered_precondition);
2130 filtered_precondition.
clear();
2132 for (
unsigned int i=0; i<dof_to_boundary_mapping.size(); ++i)
2134 && ! excluded_dofs[dof_to_boundary_mapping[i]])
2145 boundary_values[i] = boundary_projection(dof_to_boundary_mapping[i]);
2150 template <
int dim,
int spacedim>
2156 std::map<types::global_dof_index,double> &boundary_values,
2157 std::vector<unsigned int> component_mapping)
2159 do_project_boundary_values(mapping, dof, boundary_functions, q,
2160 boundary_values, component_mapping);
2165 template <
int dim,
int spacedim>
2170 std::map<types::global_dof_index,double> &boundary_values,
2171 std::vector<unsigned int> component_mapping)
2174 boundary_values, component_mapping);
2179 template <
int dim,
int spacedim>
2184 std::map<types::global_dof_index,double> &boundary_values,
2185 std::vector<unsigned int> component_mapping)
2187 do_project_boundary_values (mapping, dof,
2195 template <
int dim,
int spacedim>
2199 std::map<types::global_dof_index,double> &boundary_values,
2200 std::vector<unsigned int> component_mapping)
2213 template <
int dim,
int spacedim>
2220 std::vector<unsigned int> component_mapping)
2222 std::map<types::global_dof_index,double> boundary_values;
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)
2231 constraints.
add_line (boundary_value->first);
2233 boundary_value->second);
2240 template <
int dim,
int spacedim>
2246 std::vector<unsigned int> component_mapping)
2249 constraints, component_mapping);
2270 for (
unsigned int i=0; i<dim; ++i)
2275 bool operator < (const VectorDoFTuple<dim> &other)
const
2277 for (
unsigned int i=0; i<dim; ++i)
2278 if (dof_indices[i] < other.dof_indices[i])
2280 else if (dof_indices[i] > other.dof_indices[i])
2287 for (
unsigned int i=0; i<dim; ++i)
2288 if (dof_indices[i] != other.dof_indices[i])
2296 return ! (*
this == other);
2305 for (
unsigned int d=0; d<dim; ++d)
2306 out << vdt.dof_indices[d] << (d < dim-1 ?
" " :
"");
2329 add_constraint (
const VectorDoFTuple<dim> &dof_indices,
2399 if (std::fabs(constraining_vector[0]) > std::fabs(constraining_vector[1]) + 1e-10)
2405 constraints.
add_line (dof_indices.dof_indices[0]);
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]);
2420 constraints.
add_line (dof_indices.dof_indices[1]);
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]);
2434 if ((std::fabs(constraining_vector[0]) >= std::fabs(constraining_vector[1])+1e-10)
2436 (std::fabs(constraining_vector[0]) >= std::fabs(constraining_vector[2])+2e-10))
2442 constraints.
add_line (dof_indices.dof_indices[0]);
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]);
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]);
2457 else if ((std::fabs(constraining_vector[1])+1e-10 >= std::fabs(constraining_vector[0]))
2459 (std::fabs(constraining_vector[1]) >= std::fabs(constraining_vector[2])+1e-10))
2465 constraints.
add_line (dof_indices.dof_indices[1]);
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]);
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]);
2486 constraints.
add_line (dof_indices.dof_indices[2]);
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]);
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]);
2530 add_tangentiality_constraints (
const VectorDoFTuple<dim> &dof_indices,
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;
2553 for (
unsigned int d=0; d<dim; ++d)
2554 if (d != largest_component)
2559 constraints.
add_line (dof_indices.dof_indices[d]);
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]);
2615 if ((std::fabs(tmp[0]) > std::fabs(tmp[1]))
2617 (std::fabs(tmp[0]) > std::fabs(tmp[2])))
2622 if ((std::fabs(tmp[1]) > std::fabs(tmp[2])))
2623 std::swap (tmp[0], tmp[1]);
2625 std::swap (tmp[0], tmp[2]);
2629 else if ((std::fabs(tmp[1]) > std::fabs(tmp[0]))
2631 (std::fabs(tmp[1]) > std::fabs(tmp[2])))
2636 if ((std::fabs(tmp[0]) > std::fabs(tmp[2])))
2637 std::swap (tmp[1], tmp[0]);
2639 std::swap (tmp[1], tmp[2]);
2648 if ((std::fabs(tmp[0]) > std::fabs(tmp[1])))
2649 std::swap (tmp[2], tmp[0]);
2651 std::swap (tmp[2], tmp[1]);
2665 cross_product (orthonormals[0], vector, tmp);
2666 cross_product (orthonormals[1], vector, orthonormals[0]);
2683 template<
typename cell_iterator>
2685 compute_edge_projection (
const cell_iterator &cell,
2686 const unsigned int face,
2687 const unsigned int line,
2690 const unsigned int first_vector_component,
2691 std::vector<double> &dof_values,
2692 std::vector<bool> &dofs_processed)
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;
2708 const std::vector< DerivativeForm<1,dim,spacedim> > &
2710 const std::vector<Point<dim> > &
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);
2725 if (
dynamic_cast<const FESystem<dim>*
> (&cell->get_fe ()) != 0)
2727 unsigned int fe_index = 0;
2728 unsigned int fe_index_old = 0;
2733 fe_index_old = fe_index;
2736 if (fe_index >= first_vector_component)
2740 base_indices.first = i;
2747 edge_coordinate_direction
2771 Point<dim> shifted_reference_point_1 = reference_quadrature_points[q_point];
2772 Point<dim> shifted_reference_point_2 = reference_quadrature_points[q_point];
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]
2779 .transform_unit_to_real_cell (cell,
2780 shifted_reference_point_1)
2783 .transform_unit_to_real_cell (cell,
2784 shifted_reference_point_2))
2786 tangentials[q_point]
2787 /= std::sqrt (tangentials[q_point].square ());
2798 ((line + 1) * fe.
degree - 1, face)))
2800 && (i < (line + 1) * fe.
degree)))
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]]));
2816 dofs_processed[i] =
true;
2824 template<
int dim,
typename cell_iterator>
2826 compute_edge_projection (
const cell_iterator &,
2832 std::vector<double> &,
2833 std::vector<bool> &)
2842 template<
int dim,
typename cell_iterator>
2844 compute_face_projection_curl_conforming (
const cell_iterator &cell,
2845 const unsigned int face,
2848 const unsigned int first_vector_component,
2849 std::vector<double> &dof_values,
2850 std::vector<bool> &dofs_processed)
2852 const unsigned int spacedim = dim;
2853 hp_fe_values.
reinit (cell, cell->active_fe_index ()
2860 const std::vector< DerivativeForm<1,dim,spacedim> > &
2862 const std::vector<Point<dim> > &
2864 const unsigned int degree = fe.
degree - 1;
2865 std::pair<unsigned int, unsigned int> base_indices (0, 0);
2867 if (
dynamic_cast<const FESystem<dim>*
> (&cell->get_fe ()) != 0)
2869 unsigned int fe_index = 0;
2870 unsigned int fe_index_old = 0;
2875 fe_index_old = fe_index;
2878 if (fe_index >= first_vector_component)
2882 base_indices.first = i;
2886 std::vector<Vector<double> >
2898 const double tol = 0.5 * cell->face (face)->diameter () / cell->get_fe ().degree;
2899 std::vector<Point<dim> >
2902 const std::vector<Point<dim> > &
2903 reference_quadrature_points = fe_values.
get_quadrature ().get_points ();
2919 for (
unsigned int q_point = 0;
2927 = reference_quadrature_points[q_point];
2929 = reference_quadrature_points[q_point];
2931 shifted_reference_point_1 (face_coordinate_direction[face])
2933 shifted_reference_point_2 (face_coordinate_direction[face])
2935 tangentials[q_point]
2937 .transform_unit_to_real_cell (cell,
2938 shifted_reference_point_1)
2941 .transform_unit_to_real_cell (cell,
2942 shifted_reference_point_2))
2944 tangentials[q_point]
2945 /= std::sqrt (tangentials[q_point].square ());
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))
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]]);
2967 dofs_processed[i] =
true;
2978 assembling_matrix (degree * fe.
degree,
2983 assembling_matrix.m ());
2985 assembling_matrix.m ());
3006 for (
unsigned int q_point = 0;
3020 for (
unsigned int d = 0; d < dim; ++d)
3021 tmp[d] = values[q_point] (first_vector_component + d);
3032 tmp -= dof_values[i] * fe_values[vec].value (fe.
face_to_cell_index (i, face), q_point);
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]]
3039 jacobians[q_point][1][global_face_coordinate_directions[face][0]]
3040 * jacobians[q_point][1][global_face_coordinate_directions[face][0]]
3042 jacobians[q_point][2][global_face_coordinate_directions[face][0]]
3043 * jacobians[q_point][2][global_face_coordinate_directions[face][0]])
3045 (jacobians[q_point][0][global_face_coordinate_directions[face][1]]
3046 * jacobians[q_point][0][global_face_coordinate_directions[face][1]]
3048 jacobians[q_point][1][global_face_coordinate_directions[face][1]]
3049 * jacobians[q_point][1][global_face_coordinate_directions[face][1]]
3051 jacobians[q_point][2][global_face_coordinate_directions[face][1]]
3052 * jacobians[q_point][2][global_face_coordinate_directions[face][1]])));
3063 for (
unsigned int d = 0; d < dim; ++d)
3064 assembling_vector (dim * q_point + d) = JxW * tmp[d];
3066 unsigned int index = 0;
3085 for (
unsigned int d = 0; d < dim; ++d)
3086 assembling_matrix (index, dim * q_point + d) = shape_value[d];
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);
3104 unsigned int index = 0;
3119 dof_values[i] = solution (index);
3120 dofs_processed[i] =
true;
3127 for (
unsigned int q_point = 0;
3132 for (
unsigned int d = 0; d < dim; ++d)
3133 tmp[d] = values[q_point] (first_vector_component + d);
3141 tmp -= dof_values[i] * fe_values[vec].value (fe.
face_to_cell_index (i, face), q_point);
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]]
3148 jacobians[q_point][1][global_face_coordinate_directions[face][0]]
3149 * jacobians[q_point][1][global_face_coordinate_directions[face][0]]
3151 jacobians[q_point][2][global_face_coordinate_directions[face][0]]
3152 * jacobians[q_point][2][global_face_coordinate_directions[face][0]])
3154 (jacobians[q_point][0][global_face_coordinate_directions[face][1]]
3155 * jacobians[q_point][0][global_face_coordinate_directions[face][1]]
3157 jacobians[q_point][1][global_face_coordinate_directions[face][1]]
3158 * jacobians[q_point][1][global_face_coordinate_directions[face][1]]
3160 jacobians[q_point][2][global_face_coordinate_directions[face][1]]
3161 * jacobians[q_point][2][global_face_coordinate_directions[face][1]])));
3163 for (
unsigned int d = 0; d < dim; ++d)
3164 assembling_vector (dim * q_point + d) = JxW * tmp[d];
3166 unsigned int index = 0;
3181 for (
unsigned int d = 0; d < dim; ++d)
3182 assembling_matrix (index, dim * q_point + d) = shape_value[d];
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);
3193 unsigned int index = 0;
3204 dof_values[i] = solution (index);
3205 dofs_processed[i] =
true;
3225 const unsigned int first_vector_component,
3243 const unsigned int superdegree = dof_handler.
get_fe ().
degree;
3244 const QGauss<dim - 1> reference_face_quadrature (2 * superdegree);
3250 for (
unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
3255 face_quadrature_collection,
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);
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)
3293 typename FEL::ExcInterpolationNotImplemented ());
3296 for (
unsigned int dof = 0; dof < dofs_per_face; ++dof)
3298 dof_values[dof] = 0.0;
3299 dofs_processed[dof] =
false;
3307 ::compute_face_projection_curl_conforming (cell, face, fe_face_values,
3309 first_vector_component,
3310 dof_values, dofs_processed);
3311 cell->face (face)->get_dof_indices (face_dof_indices,
3312 cell->active_fe_index ());
3320 for (
unsigned int dof = 0; dof < dofs_per_face; ++dof)
3321 if (dofs_processed[dof] && constraints.
can_store_line (face_dof_indices[dof])
3324 constraints.
add_line (face_dof_indices[dof]);
3326 if (std::abs (dof_values[dof]) > 1e-13)
3336 const QGauss<dim - 2> reference_edge_quadrature (2 * superdegree);
3337 const unsigned int degree = superdegree - 1;
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)
3345 (reference_edge_quadrature, line), face));
3348 edge_quadrature_collection,
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)
3377 typename FEL::ExcInterpolationNotImplemented ());
3380 for (
unsigned int dof = 0; dof < dofs_per_face; ++dof)
3382 dof_values[dof] = 0.0;
3383 dofs_processed[dof] =
false;
3389 for (
unsigned int line = 0;
3390 line < GeometryInfo<3>::lines_per_face; ++line)
3392 ::compute_edge_projection (cell, face, line, fe_edge_values,
3394 first_vector_component,
3395 dof_values, dofs_processed);
3404 ::compute_face_projection_curl_conforming (cell, face, fe_face_values,
3406 first_vector_component,
3413 cell->face (face)->get_dof_indices (face_dof_indices,
3414 cell->active_fe_index ());
3416 for (
unsigned int dof = 0; dof < dofs_per_face; ++dof)
3417 if (dofs_processed[dof] && constraints.
can_store_line (face_dof_indices[dof])
3420 constraints.
add_line (face_dof_indices[dof]);
3422 if (std::abs (dof_values[dof]) > 1e-13)
3441 const unsigned int first_vector_component,
3450 for (
unsigned int i = 0; i < fe_collection.size (); ++i)
3453 reference_face_quadrature (2 * fe_collection[i].degree);
3455 for (
unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
3461 face_quadrature_collection,
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 ();
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)
3491 typename FEL::ExcInterpolationNotImplemented ());
3494 const unsigned int dofs_per_face = cell->
get_fe ().dofs_per_face;
3496 dofs_processed.resize (dofs_per_face);
3497 dof_values.resize (dofs_per_face);
3499 for (
unsigned int dof = 0; dof < dofs_per_face; ++dof)
3501 dof_values[dof] = 0.0;
3502 dofs_processed[dof] =
false;
3506 ::compute_face_projection_curl_conforming (cell, face, fe_face_values,
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 ());
3514 for (
unsigned int dof = 0; dof < dofs_per_face; ++dof)
3515 if (dofs_processed[dof] && constraints.
can_store_line (face_dof_indices[dof])
3518 constraints.
add_line (face_dof_indices[dof]);
3520 if (std::abs (dof_values[dof]) > 1e-13)
3532 for (
unsigned int i = 0; i < fe_collection.size (); ++i)
3535 reference_edge_quadrature (2 * fe_collection[i].degree);
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)
3546 edge_quadrature_collection,
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)
3568 typename FEL::ExcInterpolationNotImplemented ());
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;
3575 dofs_processed.resize (dofs_per_face);
3576 dof_values.resize (dofs_per_face);
3578 for (
unsigned int dof = 0; dof < dofs_per_face; ++dof)
3580 dof_values[dof] = 0.0;
3581 dofs_processed[dof] =
false;
3584 for (
unsigned int line = 0;
3585 line < GeometryInfo<dim>::lines_per_face; ++line)
3587 ::compute_edge_projection (cell, face, line, fe_edge_values,
3589 first_vector_component,
3590 dof_values, dofs_processed);
3596 ::compute_face_projection_curl_conforming (cell, face, fe_face_values,
3598 first_vector_component,
3599 dof_values, dofs_processed);
3602 face_dof_indices.resize (dofs_per_face);
3603 cell->face (face)->get_dof_indices (face_dof_indices,
3604 cell->active_fe_index ());
3606 for (
unsigned int dof = 0; dof < dofs_per_face; ++dof)
3607 if (dofs_processed[dof] && constraints.
can_store_line (face_dof_indices[dof])
3610 constraints.
add_line (face_dof_indices[dof]);
3612 if (std::abs (dof_values[dof]) > 1e-13)
3630 template <
typename cell_iterator>
3632 compute_face_projection_div_conforming (
const cell_iterator &cell,
3633 const unsigned int face,
3635 const unsigned int first_vector_component,
3648 std::vector<Vector<double> >
3654 const std::vector<Point<2> > &
3664 for (
unsigned int d = 0; d < 2; ++d)
3665 tmp += normals[q_point][d] * values[q_point] (d);
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]]);
3674 dof_values (i) += tmp * (normals[q_point]
3678 std::vector<types::global_dof_index> face_dof_indices (fe.
dofs_per_face);
3680 cell->face (face)->get_dof_indices (face_dof_indices, cell->active_fe_index ());
3687 constraints.
add_line (face_dof_indices[i]);
3689 if (std::abs (dof_values (i)) > 1e-14)
3695 template<
int dim,
typename cell_iterator>
3697 compute_face_projection_div_conforming (
const cell_iterator &,
3710 template<
typename cell_iterator>
3712 compute_face_projection_div_conforming (
const cell_iterator &cell,
3713 const unsigned int face,
3715 const unsigned int first_vector_component,
3718 std::vector<double> &dof_values,
3719 std::vector<types::global_dof_index> &projected_dofs)
3735 std::vector<Vector<double> >
3740 const std::vector<Point<3> > &
3750 for (
unsigned int d = 0; d < 3; ++d)
3751 tmp += normals[q_point][d] * values[q_point] (d);
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]]));
3768 dof_values_local (i) += tmp * (normals[q_point]
3772 std::vector<types::global_dof_index> face_dof_indices (fe.
dofs_per_face);
3774 cell->face (face)->get_dof_indices (face_dof_indices, cell->active_fe_index ());
3777 if (projected_dofs[face_dof_indices[i]] < fe.
degree)
3779 dof_values[face_dof_indices[i]] = dof_values_local (i);
3780 projected_dofs[face_dof_indices[i]] = fe.
degree;
3787 template<
int dim,
typename cell_iterator>
3789 compute_face_projection_div_conforming (
const cell_iterator &,
3795 std::vector<double> &,
3796 std::vector<types::global_dof_index> &)
3806 const unsigned int first_vector_component,
3812 const unsigned int spacedim = dim;
3834 for (
unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
3838 hp::FEValues<dim> fe_values (mapping_collection, fe_collection, quadrature_collection,
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)
3864 if (
dynamic_cast<const FESystem<dim>*
> (&cell->get_fe ()) == 0)
3869 typename FEL::ExcInterpolationNotImplemented ());
3872 fe_values.
reinit (cell, face + cell->active_fe_index ()
3875 const std::vector<DerivativeForm<1,dim,spacedim> > &
3878 fe_face_values.
reinit (cell, face);
3879 internals::compute_face_projection_div_conforming (cell, face,
3881 first_vector_component,
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);
3913 for (
unsigned int dof = 0; dof < n_dofs; ++dof)
3914 projected_dofs[dof] = 0;
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)
3928 if (
dynamic_cast<const FESystem<dim>*
> (&cell->get_fe ()) == 0)
3933 typename FEL::ExcInterpolationNotImplemented ());
3936 fe_values.
reinit (cell, face + cell->active_fe_index ()
3939 const std::vector<DerivativeForm<1,dim ,spacedim> > &
3942 fe_face_values.
reinit (cell, face);
3943 internals::compute_face_projection_div_conforming (cell, face,
3945 first_vector_component,
3947 jacobians, dof_values,
3951 for (
unsigned int dof = 0; dof < n_dofs; ++dof)
3952 if ((projected_dofs[dof] != 0) && !(constraints.
is_constrained (dof)))
3956 if (std::abs (dof_values[dof]) > 1e-14)
3972 const unsigned int first_vector_component,
3978 const unsigned int spacedim = dim;
3983 for (
unsigned int i = 0; i < fe_collection.
size (); ++i)
3985 const QGauss<dim - 1> quadrature (fe_collection[i].degree + 1);
3987 face_quadrature_collection.push_back (quadrature);
3989 for (
unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
3999 hp::FEValues<dim> fe_values (mapping_collection, fe_collection, quadrature_collection,
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)
4018 if (
dynamic_cast<const FESystem<dim>*
> (&cell->get_fe ()) == 0)
4023 typename FEL::ExcInterpolationNotImplemented ());
4026 fe_values.
reinit (cell, face + cell->active_fe_index ()
4029 const std::vector<DerivativeForm<1,dim,spacedim> > &
4032 fe_face_values.
reinit (cell, face);
4033 internals::compute_face_projection_div_conforming (cell, face,
4035 first_vector_component,
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);
4050 for (
unsigned int dof = 0; dof < n_dofs; ++dof)
4051 projected_dofs[dof] = 0;
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)
4065 if (
dynamic_cast<const FESystem<dim>*
> (&cell->get_fe ()) == 0)
4070 typename FEL::ExcInterpolationNotImplemented ());
4073 fe_values.
reinit (cell, face + cell->active_fe_index ()
4076 const std::vector<DerivativeForm<1,dim,spacedim> > &
4079 fe_face_values.
reinit (cell, face);
4080 internals::compute_face_projection_div_conforming (cell, face,
4082 first_vector_component,
4084 jacobians, dof_values,
4088 for (
unsigned int dof = 0; dof < n_dofs; ++dof)
4089 if ((projected_dofs[dof] != 0) && !(constraints.
is_constrained (dof)))
4093 if (std::abs (dof_values[dof]) > 1e-14)
4107 template <
int dim,
template <
int,
int>
class DH,
int spacedim>
4110 const unsigned int first_vector_component,
4111 const std::set<types::boundary_id> &boundary_ids,
4116 ExcMessage (
"This function is not useful in 1d because it amounts "
4117 "to imposing Dirichlet values on the vector-valued "
4120 std::vector<types::global_dof_index> face_dofs;
4126 for (
unsigned int i=0; i<fe_collection.size(); ++i)
4133 for (
unsigned int i=0; i<fe_collection.size(); ++i)
4135 const std::vector<
Point<dim-1> > &
4136 unit_support_points = fe_collection[i].get_unit_face_support_points();
4138 Assert (unit_support_points.size() == fe_collection[i].dofs_per_face,
4141 face_quadrature_collection
4148 face_quadrature_collection,
4164 std::multimap<internal::VectorDoFTuple<dim>,
4165 std::pair<Tensor<1,dim>,
typename DH<dim,spacedim>::active_cell_iterator> >
4168 DoFToNormalsMap dof_to_normals_map;
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;
4178 if (boundary_ids.find(cell->face(face_no)->boundary_indicator())
4179 != boundary_ids.end())
4182 typename DH<dim,spacedim>::face_iterator face = cell->face(face_no);
4186 face->get_dof_indices (face_dofs, cell->active_fe_index());
4188 x_fe_face_values.
reinit (cell, face_no);
4193 for (
unsigned int i=0; i<face_dofs.size(); ++i)
4195 first_vector_component)
4199 vector_dofs.dof_indices[0] = face_dofs[i];
4202 ExcMessage(
"Error: the finite element does not have enough components "
4203 "to define a normal direction."));
4208 (face_quadrature_collection[cell->active_fe_index()].point(k) ==
4209 face_quadrature_collection[cell->active_fe_index()].point(i))
4212 first_vector_component)
4215 first_vector_component + dim))
4217 first_vector_component]
4220 for (
unsigned int d=0; d<dim; ++d)
4221 Assert (vector_dofs.dof_indices[d] < dof_handler.n_dofs(),
4271 = (cell->face(face_no)->get_boundary()
4272 .normal_vector (cell->face(face_no),
4275 normal_vector *= -1;
4276 Assert (std::fabs(normal_vector.
norm() - 1) < 1e-14,
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();
4286 .insert (std::make_pair (vector_dofs,
4287 std::make_pair (normal_vector,
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;
4303 typename DoFToNormalsMap::const_iterator
4304 p = dof_to_normals_map.begin();
4306 while (p != dof_to_normals_map.end())
4311 typename DoFToNormalsMap::const_iterator same_dof_range[2]
4313 for (++p; p != dof_to_normals_map.end(); ++p)
4314 if (p->first != same_dof_range[0]->first)
4316 same_dof_range[1] = p;
4319 if (p == dof_to_normals_map.end())
4320 same_dof_range[1] = dof_to_normals_map.end();
4322 #ifdef DEBUG_NO_NORMAL_FLUX
4323 std::cout <<
"For dof indices <" << p->first <<
">, found the following normals"
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
4340 <
typename DH<dim,spacedim>::active_cell_iterator,
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);
4355 = cell_to_normals_map[q->second.second].first;
4356 const unsigned int old_count
4357 = cell_to_normals_map[q->second.second].second;
4363 cell_to_normals_map[q->second.second]
4364 = std::make_pair ((old_normal * old_count + q->second.first) / (old_count + 1),
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 <<
')'
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,
4393 switch (max_n_contributions_per_cell)
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();
4421 for (
unsigned int d=0; d<dim; ++d)
4422 if (std::fabs(normal[d]) < 1e-13)
4424 normal /= normal.
norm();
4428 dof_indices = same_dof_range[0]->first;
4429 internal::add_constraint (dof_indices, normal,
4448 Assert (cell_to_normals_map.size() == 1,
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];
4465 Assert (std::fabs(determinant (t)) > 1e-3,
4466 ExcMessage(
"Found a set of normal vectors that are nearly collinear."));
4473 for (
unsigned int i=0; i<dim; ++i)
4475 ->first.dof_indices[i])
4478 same_dof_range[0]->first.dof_indices[i]))
4480 constraints.
add_line (same_dof_range[0]->first.dof_indices[i]);
4501 std::map<typename DH<dim,spacedim>::active_cell_iterator, std::list<Tensor<1,dim> > >
4503 CellContributions cell_contributions;
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);
4534 std::list<Tensor<1,dim> > tangential_vectors;
4535 for (
typename CellContributions::const_iterator
4536 contribution = cell_contributions.begin();
4537 contribution != cell_contributions.end();
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:"
4546 t = contribution->second.begin();
4547 t != contribution->second.end();
4549 std::cout <<
" " << *t << std::endl;
4554 if (contribution->second.size() < dim-1)
4559 unsigned int index=0;
4561 t = contribution->second.begin();
4562 t != contribution->second.end();
4564 normals[index] = *t;
4586 cross_product (tangent, normals[0], normals[dim-2]);
4593 ExcMessage(
"Two normal vectors from adjacent faces are almost "
4595 tangent /= tangent.
norm();
4597 tangential_vectors.push_back (tangent);
4605 const Tensor<1,dim> first_tangent = tangential_vectors.front();
4606 typename std::list<Tensor<1,dim> >::iterator
4607 t = tangential_vectors.begin();
4609 for (; t != tangential_vectors.end(); ++t)
4610 if (*t * first_tangent < 0)
4617 t = tangential_vectors.begin();
4618 t != tangential_vectors.end();
4620 average_tangent += *t;
4621 average_tangent /= average_tangent.
norm();
4626 dof_indices = same_dof_range[0]->first;
4627 internal::add_tangentiality_constraints (dof_indices,
4640 struct PointComparator
4642 bool operator ()(
const std_cxx1x::array<types::global_dof_index,dim> &p1,
4643 const std_cxx1x::array<types::global_dof_index,dim> &p2)
4645 for (
unsigned int d=0; d<dim; ++d)
4655 template <
int dim,
template <
int,
int>
class DH,
int spacedim>
4658 const unsigned int first_vector_component,
4659 const std::set<types::boundary_id> &boundary_ids,
4665 first_vector_component,
4667 no_normal_flux_constraints,
4673 std::set<std_cxx1x::array<types::global_dof_index,dim>, PointComparator<dim> > vector_dofs;
4674 std::vector<types::global_dof_index> face_dofs;
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;
4682 if (boundary_ids.find(cell->face(face_no)->boundary_indicator())
4683 != boundary_ids.end())
4686 typename DH<dim,spacedim>::face_iterator face=cell->face(face_no);
4690 face->get_dof_indices (face_dofs, cell->active_fe_index());
4692 unsigned int n_scalar_indices = 0;
4699 std::max(n_scalar_indices,
4709 for (
unsigned int i=0; i<n_scalar_indices; ++i)
4710 vector_dofs.insert(cell_vector_dofs[i]);
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)
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]))
4725 is_constrained[d] =
true;
4727 ++n_total_constraints_found;
4730 is_constrained[d] =
false;
4731 if (n_constraints > 0)
4737 if (n_constraints > 1)
4739 for (
unsigned int d=0; d<dim; ++d)
4749 unsigned constrained_index = -1;
4750 for (
unsigned int d=0; d<dim; ++d)
4751 if (is_constrained[d])
4753 constrained_index = d;
4757 const std::vector<std::pair<types::global_dof_index, double> > *constrained
4758 = no_normal_flux_constraints.get_constraint_entries((*it)[constrained_index]);
4762 for (
unsigned int c=0; c<constrained->size(); ++c)
4765 for (
unsigned int d=0; d<dim; ++d)
4766 if ((*constrained)[c].first == (*it)[d])
4769 normal[index] = (*constrained)[c].second;
4771 for (
unsigned int d=0; d<dim; ++d)
4773 if (is_constrained[d])
4775 const unsigned int new_index = (*it)[d];
4779 if (std::abs(normal[d]) > 1e-13)
4780 constraints.
add_entry(new_index, (*it)[constrained_index],
4787 no_normal_flux_constraints.n_constraints());
4794 template <
int dim,
int spacedim>
4797 IDScratchData (const ::hp::MappingCollection<dim,spacedim> &mapping,
4798 const ::hp::FECollection<dim,spacedim> &fe,
4799 const ::hp::QCollection<dim> &q,
4804 void resize_vectors (
const unsigned int n_q_points,
4805 const unsigned int n_components);
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;
4812 std::vector<::Vector<double> > psi_values;
4813 std::vector<std::vector<Tensor<1,spacedim> > > psi_grads;
4814 std::vector<double> psi_scalar;
4816 std::vector<double> tmp_values;
4817 std::vector<Tensor<1,spacedim> > tmp_gradients;
4823 template <
int dim,
int spacedim>
4826 const ::hp::FECollection<dim,spacedim> &fe,
4827 const ::hp::QCollection<dim> &q,
4830 x_fe_values(mapping, fe, q, update_flags)
4833 template <
int dim,
int spacedim>
4834 IDScratchData<dim,spacedim>::IDScratchData (
const IDScratchData &data)
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())
4842 template <
int dim,
int spacedim>
4844 IDScratchData<dim,spacedim>::resize_vectors (
const unsigned int n_q_points,
4845 const unsigned int n_components)
4847 function_values.resize (n_q_points,
4849 function_grads.resize (n_q_points,
4852 weight_values.resize (n_q_points);
4853 weight_vectors.resize (n_q_points,
4856 psi_values.resize (n_q_points,
4858 psi_grads.resize (n_q_points,
4860 psi_scalar.resize (n_q_points);
4862 tmp_values.resize (n_q_points);
4863 tmp_gradients.resize (n_q_points);
4870 template <
int dim,
int spacedim>
4876 const double exponent,
4877 const unsigned int n_components,
4878 IDScratchData<dim,spacedim> &data)
4880 const bool fe_is_system = (n_components != 1);
4882 const unsigned int n_q_points = fe_values.n_quadrature_points;
4888 data.weight_vectors);
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];
4899 for (
unsigned int k=0; k<n_q_points; ++k)
4900 data.weight_vectors[k] = 1.;
4915 exact_solution.
value_list (fe_values.get_quadrature_points(),
4917 for (
unsigned int i=0; i<n_q_points; ++i)
4918 data.psi_values[i](0) = data.tmp_values[i];
4922 for (
unsigned int q=0; q<n_q_points; ++q)
4923 data.psi_values[q] -= data.function_values[q];
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];
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]*
4953 fe_values.normal_vector(q))*
4954 fe_values.normal_vector(q));
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];
4966 for (
unsigned int q=0; q<n_q_points; ++q)
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);
4979 for (
unsigned int q=0; q<n_q_points; ++q)
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);
4989 if (!(update_flags & update_gradients))
4990 diff = std::pow(diff, 1./exponent);
4996 for (
unsigned int q=0; q<n_q_points; ++q)
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);
5006 diff=std::sqrt(diff);
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)));
5031 for (
unsigned int q=0; q<n_q_points; ++q)
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);
5039 diff = std::pow(diff, 1./exponent);
5044 for (
unsigned int q=0; q<n_q_points; ++q)
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);
5052 diff = std::sqrt(diff);
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));
5080 template <
int dim,
class InVector,
class OutVector,
class DH,
int spacedim>
5083 do_integrate_difference (const ::hp::MappingCollection<dim,spacedim> &mapping,
5085 const InVector &fe_function,
5087 OutVector &difference,
5088 const ::hp::QCollection<dim> &q,
5091 const double exponent_1)
5096 double exponent = exponent_1;
5098 const unsigned int n_components = dof.get_fe().n_components();
5106 difference.reinit (dof.get_tria().n_active_cells());
5132 if (spacedim == dim+1)
5133 update_flags |=
UpdateFlags (update_normal_vectors);
5141 if (spacedim == dim+1)
5142 update_flags |=
UpdateFlags (update_normal_vectors);
5151 IDScratchData<dim,spacedim> data(mapping, fe_collection, q, update_flags);
5154 typename DH::active_cell_iterator cell = dof.begin_active(),
5156 for (
unsigned int index=0; cell != endc; ++cell, ++index)
5157 if (cell->is_locally_owned())
5160 data.x_fe_values.reinit (cell);
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);
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);
5172 integrate_difference_inner (exact_solution, norm, weight,
5173 update_flags, exponent,
5174 n_components, data);
5179 difference(index) = 0;
5187 template <
int dim,
class InVector,
class OutVector,
int spacedim>
5191 const InVector &fe_function,
5193 OutVector &difference,
5197 const double exponent)
5201 dof, fe_function, exact_solution,
5203 norm, weight, exponent);
5207 template <
int dim,
class InVector,
class OutVector,
int spacedim>
5210 const InVector &fe_function,
5212 OutVector &difference,
5216 const double exponent)
5220 dof, fe_function, exact_solution,
5222 norm, weight, exponent);
5227 template <
int dim,
class InVector,
class OutVector,
int spacedim>
5230 const ::hp::DoFHandler<dim,spacedim> &dof,
5231 const InVector &fe_function,
5233 OutVector &difference,
5234 const ::hp::QCollection<dim> &q,
5237 const double exponent)
5241 dof, fe_function, exact_solution,
5243 norm, weight, exponent);
5247 template <
int dim,
class InVector,
class OutVector,
int spacedim>
5250 const InVector &fe_function,
5252 OutVector &difference,
5253 const ::hp::QCollection<dim> &q,
5256 const double exponent)
5260 dof, fe_function, exact_solution,
5262 norm, weight, exponent);
5267 template <
int dim,
class InVector,
int spacedim>
5270 const InVector &fe_function,
5284 template <
int dim,
class InVector,
int spacedim>
5288 const InVector &fe_function,
5301 const std::pair<typename DoFHandler<dim,spacedim>::active_cell_iterator,
Point<spacedim> >
5305 ExcPointNotAvailableHere());
5311 FEValues<dim> fe_values(mapping, fe, quadrature, update_values);
5312 fe_values.
reinit(cell_point.first);
5320 difference(0) = exact_function.
value(point);
5324 for (
unsigned int i=0; i<difference.
size(); ++i)
5325 difference(i) -= u_value[0](i);
5329 template <
int dim,
class InVector,
int spacedim>
5332 const InVector &fe_function,
5345 template <
int dim,
class InVector,
int spacedim>
5348 const InVector &fe_function,
5360 template <
int dim,
class InVector,
int spacedim>
5363 const InVector &fe_function,
5373 template <
int dim,
class InVector,
int spacedim>
5376 const InVector &fe_function,
5386 template <
int dim,
class InVector,
int spacedim>
5390 const InVector &fe_function,
5402 const std::pair<typename DoFHandler<dim,spacedim>::active_cell_iterator,
Point<spacedim> >
5407 ExcPointNotAvailableHere());
5414 FEValues<dim> fe_values(mapping, fe, quadrature, update_values);
5415 fe_values.
reinit(cell_point.first);
5426 template <
int dim,
class InVector,
int spacedim>
5430 const InVector &fe_function,
5442 const std::pair<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator,
Point<spacedim> >
5446 ExcPointNotAvailableHere());
5453 hp_fe_values.
reinit(cell_point.first);
5465 template <
int dim,
class InVector,
int spacedim>
5469 const InVector &fe_function,
5475 ExcMessage (
"Finite element is not scalar as is necessary for this function"));
5480 const std::pair<typename DoFHandler<dim,spacedim>::active_cell_iterator,
Point<spacedim> >
5484 ExcPointNotAvailableHere());
5490 FEValues<dim> fe_values(mapping, fe, quadrature, update_values);
5491 fe_values.
reinit(cell_point.first);
5495 std::vector<double> u_value(1);
5502 template <
int dim,
class InVector,
int spacedim>
5506 const InVector &fe_function,
5512 ExcMessage (
"Finite element is not scalar as is necessary for this function"));
5517 const std::pair<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator,
Point<spacedim> >
5521 ExcPointNotAvailableHere());
5528 hp_fe_values.
reinit(cell_point.first);
5533 std::vector<double> u_value(1);
5541 template <
class VECTOR>
5544 const std::vector<bool> &p_select)
5546 if (p_select.size() == 0)
5549 v.add( - v.mean_value() );
5563 const unsigned int n = v.size();
5565 Assert(p_select.size() == n,
5568 typename VECTOR::value_type s = 0.;
5569 unsigned int counter = 0;
5570 for (
unsigned int i=0; i<n; ++i)
5578 Assert (n == 0 || counter > 0, ComponentMask::ExcNoComponentSelected());
5582 for (
unsigned int i=0; i<n; ++i)
5590 template <
int dim,
class InVector,
int spacedim>
5596 const unsigned int component)
5607 typename DoFHandler<dim,spacedim>::active_cell_iterator cell;
5608 std::vector<Vector<double> > values(quadrature.
size(),
5615 if (cell->is_locally_owned())
5618 fe.get_function_values(v, values);
5619 for (
unsigned int k=0; k< quadrature.
size(); ++k)
5621 mean += fe.JxW(k) * values[k](component);
5626 #ifdef DEAL_II_WITH_P4EST
5633 double my_values[2] = {
mean, area };
5634 double global_values[2];
5636 MPI_Allreduce (&my_values, &global_values, 2, MPI_DOUBLE,
5638 p_d_triangulation->get_communicator());
5640 mean = global_values[0];
5641 area = global_values[1];
5649 template <
int dim,
class InVector,
int spacedim>
5654 const unsigned int component)
5660 DEAL_II_NAMESPACE_CLOSE
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
double get_inhomogeneity(const size_type line) const
static const unsigned int invalid_unsigned_int
#define AssertDimension(dim1, dim2)
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 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)
unsigned int n_selected_components(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
const unsigned int dofs_per_cell
const Point< spacedim > & normal_vector(const unsigned int i) const
std::map< types::boundary_id, const Function< dim > * > type
const std::vector< Point< dim-1 > > & get_unit_face_support_points() const
const std::vector< std::pair< size_type, double > > * get_constraint_entries(const size_type line) const
::ExceptionBase & ExcMessage(std::string arg1)
#define AssertIndexRange(index, range)
cell_iterator end() const
const unsigned int degree
void initialize(const MATRIX &m, bool expect_constrained_source=false)
const FEFaceValues< dim, spacedim > & get_present_fe_values() const
const IndexSet & get_local_lines() const
unsigned int size() const
Transformed quadrature points.
#define AssertThrow(cond, exc)
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
bool is_primitive(const unsigned int i) const
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
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
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)
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
T sum(const T &t, const MPI_Comm &mpi_communicator)
unsigned int global_dof_index
virtual void reinit(const size_type N, const bool fast=false)
#define Assert(cond, exc)
types::global_dof_index n_dofs() const
void condense(const SparsityPattern &uncondensed, SparsityPattern &condensed) const
const ::FEValues< dim, spacedim > & get_present_fe_values() const
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
const ComponentMask & get_nonzero_components(const unsigned int i) const
::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 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
unsigned int n_components() const
std::pair< unsigned int, unsigned int > face_system_to_component_index(const unsigned int index) const
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
const GridClass & get_destination_grid() const
const unsigned int n_components
const Mapping< dim, spacedim > & get_mapping() const
unsigned int first_selected_component(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
const unsigned int dofs_per_cell
const unsigned int n_quadrature_points
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 >())
std::ostream & operator<<(std::ostream &os, const Vector< number > &v)
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
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
Shape function gradients.
const unsigned int dofs_per_face
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 ¶meters=typename BaseClass::AdditionalData())
::ExceptionBase & ExcNotImplemented()
void add_constraint(const size_type i, const double v)
const GridClass & get_source_grid() const
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
::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
const types::boundary_id internal_face_boundary_id
void set_inhomogeneity(const size_type line, const double value)
const types::global_dof_index invalid_dof_index
const Triangulation< dim, spacedim > & get_tria() const
::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 push_back(const Quadrature< dim > &new_quadrature)
double JxW(const unsigned int quadrature_point) const
const std::vector< Point< spacedim > > & get_normal_vectors() const
cell_iterator end() const
const FiniteElement< dim, spacedim > & get_fe() const
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)