// $Id$
// Version: $Name: $
//
-// Copyright (C) 2005, 2006, 2007, 2008, 2009 by the deal.II authors
+// Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
Assert (dof.get_fe().n_components() == function.n_components,
ExcDimensionMismatch(dof.get_fe().n_components(),
function.n_components));
-
+
const hp::FECollection<DH::dimension,DH::space_dimension> fe (dof.get_fe());
const unsigned int n_components = fe.n_components();
const bool fe_is_system = (n_components != 1);
-
+
typename DH::active_cell_iterator cell = dof.begin_active(),
endc = dof.end();
Assert (unit_support_points[fe_index].size() != 0,
ExcNonInterpolatingFE());
}
-
- // Find the support points
+
+ // Find the support points
// on a cell that
- // are multiply mentioned in
+ // are multiply mentioned in
// unit_support_points.
// Mark the first representative
// of each multiply mentioned
// the following vector collects all dofs i,
// 0<=i<fe.dofs_per_cell, for that
- // unit_support_points[i]
+ // unit_support_points[i]
// is a representative one. i.e.
// the following vector collects all rep dofs.
// the position of a rep dof within this vector
std::vector<std::vector<unsigned int> > dof_to_rep_index_table(fe.size());
std::vector<unsigned int> n_rep_points (fe.size(), 0);
-
+
for (unsigned int fe_index=0; fe_index<fe.size(); ++fe_index)
{
for (unsigned int i=0; i<fe[fe_index].dofs_per_cell; ++i)
// support points that are placed
// one after the other.
for (unsigned int j=dofs_of_rep_points[fe_index].size(); j>0; --j)
- if (unit_support_points[fe_index][i]
+ if (unit_support_points[fe_index][i]
== unit_support_points[fe_index][dofs_of_rep_points[fe_index][j-1]])
{
dof_to_rep_index_table[fe_index].push_back(j-1);
representative=false;
break;
}
-
+
if (representative)
{
// rep_index=dofs_of_rep_points.size()
Assert(dof_to_rep_index_table[fe_index].size()==fe[fe_index].dofs_per_cell,
ExcInternalError());
}
-
+
const unsigned int max_rep_points = *std::max_element (n_rep_points.begin(),
n_rep_points.end());
std::vector<unsigned int> dofs_on_cell (fe.max_dofs_per_cell());
std::vector<Point<DH::space_dimension> > rep_points (max_rep_points);
-
+
// get space for the values of the
// function at the rep support points.
//
// Transformed support points are computed by
// FEValues
hp::MappingCollection<dim,DH::space_dimension> mapping_collection (mapping);
-
+
hp::FEValues<dim, DH::space_dimension> fe_values (mapping_collection,
fe, support_quadrature, update_quadrature_points);
-
+
for (; cell!=endc; ++cell)
{
const unsigned int fe_index = cell->active_fe_index();
-
+
// for each cell:
// get location of finite element
// support_points
fe_values.reinit(cell);
const std::vector<Point<DH::space_dimension> >& support_points =
fe_values.get_present_fe_values().get_quadrature_points();
-
+
// pick out the representative
// support points
rep_points.resize (dofs_of_rep_points[fe_index].size());
rep_points[j] = support_points[dofs_of_rep_points[fe_index][j]];
// get indices of the dofs on this cell
- dofs_on_cell.resize (fe[fe_index].dofs_per_cell);
+ dofs_on_cell.resize (fe[fe_index].dofs_per_cell);
cell->get_dof_indices (dofs_on_cell);
// values to the global
// vector
for (unsigned int i=0; i<fe[fe_index].dofs_per_cell; ++i)
- vec(dofs_on_cell[i])
+ vec(dofs_on_cell[i])
= function_values_scalar[fe_index][dof_to_rep_index_table[fe_index][i]];
}
}
VECTOR &vec)
{
Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
- interpolate(StaticMappingQ1<DH::dimension, DH::space_dimension>::mapping,
+ interpolate(StaticMappingQ1<DH::dimension, DH::space_dimension>::mapping,
dof, function, vec);
}
std::vector<short unsigned int> touch_count (dof_2.n_dofs(), 0);
std::vector<unsigned int> local_dof_indices (dof_2.get_fe().dofs_per_cell);
-
+
typename DoFHandler<dim,spacedim>::active_cell_iterator h = dof_1.begin_active();
typename DoFHandler<dim,spacedim>::active_cell_iterator l = dof_2.begin_active();
const typename DoFHandler<dim,spacedim>::cell_iterator endh = dof_1.end();
-
+
for(; h != endh; ++h, ++l)
{
h->get_dof_values(data_1, cell_data_1);
transfer.vmult(cell_data_2, cell_data_1);
l->get_dof_indices (local_dof_indices);
-
+
// distribute cell vector
- for (unsigned int j=0; j<dof_2.get_fe().dofs_per_cell; ++j)
+ for (unsigned int j=0; j<dof_2.get_fe().dofs_per_cell; ++j)
{
data_2(local_dof_indices[j]) += cell_data_2(j);
// count, how often we have
// added to this dof
Assert (touch_count[local_dof_indices[j]] < 255,
- ExcInternalError());
+ ExcInternalError());
++touch_count[local_dof_indices[j]];
};
};
{
Assert (touch_count[i] != 0,
ExcInternalError());
-
+
data_2(i) /= touch_count[i];
};
}
#else
-
+
template <int dim, int spacedim>
void
interpolate_zero_boundary_values (const dealii::DoFHandler<dim,spacedim> &dof_handler,
boundary_values[face_dof_indices[i]] = 0.;
}
}
-
+
#endif
}
}
Assert (vec_result.size() == dof.n_dofs(),
ExcDimensionMismatch (vec_result.size(), dof.n_dofs()));
-
+
// make up boundary values
std::map<unsigned int,double> boundary_values;
-
- if (enforce_zero_boundary == true)
+
+ if (enforce_zero_boundary == true)
// no need to project boundary
// values, but enforce
// homogeneous boundary values
// anyway
internal::VectorTools::
interpolate_zero_boundary_values (dof, boundary_values);
-
+
else
// no homogeneous boundary values
if (project_to_boundary_first == true)
prec.initialize(mass_matrix, 1.2);
// solve
cg.solve (mass_matrix, vec, tmp, prec);
-
+
// distribute solution
constraints.distribute (vec);
Assert (rhs_vector.size() == dof_handler.n_dofs(),
ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs()));
rhs_vector = 0;
-
+
UpdateFlags update_flags = UpdateFlags(update_values |
update_quadrature_points |
update_JxW_values);
const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
n_q_points = fe_values.n_quadrature_points,
n_components = fe.n_components();
-
+
std::vector<unsigned int> dofs (dofs_per_cell);
Vector<double> cell_vector (dofs_per_cell);
if (n_components==1)
{
std::vector<double> rhs_values(n_q_points);
-
- for (; cell!=endc; ++cell)
+
+ for (; cell!=endc; ++cell)
{
fe_values.reinit(cell);
-
+
const std::vector<double> &weights = fe_values.get_JxW_values ();
rhs_function.value_list (fe_values.get_quadrature_points(),
rhs_values);
-
+
cell_vector = 0;
for (unsigned int point=0; point<n_q_points; ++point)
- for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
cell_vector(i) += rhs_values[point] *
fe_values.shape_value(i,point) *
weights[point];
-
+
cell->get_dof_indices (dofs);
-
+
for (unsigned int i=0; i<dofs_per_cell; ++i)
rhs_vector(dofs[i]) += cell_vector(i);
}
-
+
}
else
{
std::vector<Vector<double> > rhs_values(n_q_points,
Vector<double>(n_components));
-
- for (; cell!=endc; ++cell)
+
+ for (; cell!=endc; ++cell)
{
fe_values.reinit(cell);
-
+
const std::vector<double> &weights = fe_values.get_JxW_values ();
rhs_function.vector_value_list (fe_values.get_quadrature_points(),
rhs_values);
-
+
cell_vector = 0;
// Use the faster code if the
// FiniteElement is primitive
{
const unsigned int component
= fe.system_to_component_index(i).first;
-
+
cell_vector(i) += rhs_values[point](component) *
fe_values.shape_value(i,point) *
weights[point];
weights[point];
}
}
-
+
cell->get_dof_indices (dofs);
-
+
for (unsigned int i=0; i<dofs_per_cell; ++i)
rhs_vector(dofs[i]) += cell_vector(i);
}
Assert (rhs_vector.size() == dof_handler.n_dofs(),
ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs()));
rhs_vector = 0;
-
+
UpdateFlags update_flags = UpdateFlags(update_values |
update_quadrature_points |
update_JxW_values);
hp::FEValues<dim,spacedim> x_fe_values (mapping, fe, quadrature, update_flags);
const unsigned int n_components = fe.n_components();
-
+
std::vector<unsigned int> dofs (fe.max_dofs_per_cell());
Vector<double> cell_vector (fe.max_dofs_per_cell());
if (n_components==1)
{
std::vector<double> rhs_values;
-
- for (; cell!=endc; ++cell)
+
+ for (; cell!=endc; ++cell)
{
x_fe_values.reinit(cell);
const FEValues<dim,spacedim> &fe_values = x_fe_values.get_present_fe_values();
-
+
const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
n_q_points = fe_values.n_quadrature_points;
rhs_values.resize (n_q_points);
dofs.resize (dofs_per_cell);
cell_vector.reinit (dofs_per_cell);
-
+
const std::vector<double> &weights = fe_values.get_JxW_values ();
rhs_function.value_list (fe_values.get_quadrature_points(),
rhs_values);
-
+
cell_vector = 0;
for (unsigned int point=0; point<n_q_points; ++point)
- for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
cell_vector(i) += rhs_values[point] *
fe_values.shape_value(i,point) *
weights[point];
-
+
cell->get_dof_indices (dofs);
-
+
for (unsigned int i=0; i<dofs_per_cell; ++i)
rhs_vector(dofs[i]) += cell_vector(i);
}
-
+
}
else
{
std::vector<Vector<double> > rhs_values;
-
+
for (; cell!=endc; ++cell)
{
x_fe_values.reinit(cell);
const FEValues<dim,spacedim> &fe_values = x_fe_values.get_present_fe_values();
-
+
const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
n_q_points = fe_values.n_quadrature_points;
rhs_values.resize (n_q_points,
Vector<double>(n_components));
dofs.resize (dofs_per_cell);
cell_vector.reinit (dofs_per_cell);
-
+
const std::vector<double> &weights = fe_values.get_JxW_values ();
rhs_function.vector_value_list (fe_values.get_quadrature_points(),
rhs_values);
-
+
cell_vector = 0;
-
+
// Use the faster code if the
// FiniteElement is primitive
if (cell->get_fe().is_primitive ())
{
const unsigned int component
= cell->get_fe().system_to_component_index(i).first;
-
+
cell_vector(i) += rhs_values[point](component) *
fe_values.shape_value(i,point) *
weights[point];
weights[point];
}
}
-
+
cell->get_dof_indices (dofs);
-
+
for (unsigned int i=0; i<dofs_per_cell; ++i)
rhs_vector(dofs[i]) += cell_vector(i);
}
Vector<double> &rhs_vector)
{
Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
- create_right_hand_side(hp::StaticMappingQ1<dim,spacedim>::mapping_collection,
+ create_right_hand_side(hp::StaticMappingQ1<dim,spacedim>::mapping_collection,
dof_handler, quadrature,
rhs_function, rhs_vector);
}
ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs()));
Assert (dof_handler.get_fe().n_components() == 1,
ExcMessage ("This function only works for scalar finite elements"));
-
+
rhs_vector = 0;
std::pair<typename DoFHandler<dim,spacedim>::active_cell_iterator, Point<spacedim> >
ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs()));
Assert (dof_handler.get_fe().n_components() == 1,
ExcMessage ("This function only works for scalar finite elements"));
-
+
rhs_vector = 0;
std::pair<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator, Point<spacedim> >
ExcDimensionMismatch(fe.n_components(), rhs_function.n_components));
Assert (rhs_vector.size() == dof_handler.n_dofs(),
ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs()));
-
+
rhs_vector = 0;
-
+
UpdateFlags update_flags = UpdateFlags(update_values |
update_quadrature_points |
update_JxW_values);
const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
n_q_points = fe_values.n_quadrature_points,
n_components = fe.n_components();
-
+
std::vector<unsigned int> dofs (dofs_per_cell);
Vector<double> cell_vector (dofs_per_cell);
if (n_components==1)
{
std::vector<double> rhs_values(n_q_points);
-
+
for (; cell!=endc; ++cell)
for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
if (cell->face(face)->at_boundary () &&
boundary_indicators.end()))
{
fe_values.reinit(cell, face);
-
+
const std::vector<double> &weights = fe_values.get_JxW_values ();
rhs_function.value_list (fe_values.get_quadrature_points(), rhs_values);
-
+
cell_vector = 0;
for (unsigned int point=0; point<n_q_points; ++point)
- for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
cell_vector(i) += rhs_values[point] *
fe_values.shape_value(i,point) *
weights[point];
-
+
cell->get_dof_indices (dofs);
-
+
for (unsigned int i=0; i<dofs_per_cell; ++i)
rhs_vector(dofs[i]) += cell_vector(i);
}
else
{
std::vector<Vector<double> > rhs_values(n_q_points, Vector<double>(n_components));
-
- for (; cell!=endc; ++cell)
+
+ for (; cell!=endc; ++cell)
for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
if (cell->face(face)->at_boundary () &&
(boundary_indicators.find (cell->face(face)->boundary_indicator())
boundary_indicators.end()))
{
fe_values.reinit(cell, face);
-
+
const std::vector<double> &weights = fe_values.get_JxW_values ();
rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values);
-
+
cell_vector = 0;
-
+
// Use the faster code if the
// FiniteElement is primitive
if (fe.is_primitive ())
- {
+ {
for (unsigned int point=0; point<n_q_points; ++point)
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const unsigned int component
= fe.system_to_component_index(i).first;
-
+
cell_vector(i) += rhs_values[point](component) *
fe_values.shape_value(i,point) *
weights[point];
weights[point];
}
}
-
+
cell->get_dof_indices (dofs);
-
+
for (unsigned int i=0; i<dofs_per_cell; ++i)
rhs_vector(dofs[i]) += cell_vector(i);
}
ExcDimensionMismatch(fe.n_components(), rhs_function.n_components));
Assert (rhs_vector.size() == dof_handler.n_dofs(),
ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs()));
-
+
rhs_vector = 0;
-
+
UpdateFlags update_flags = UpdateFlags(update_values |
update_quadrature_points |
update_JxW_values);
hp::FEFaceValues<dim> x_fe_values (mapping, fe, quadrature, update_flags);
const unsigned int n_components = fe.n_components();
-
+
std::vector<unsigned int> dofs (fe.max_dofs_per_cell());
Vector<double> cell_vector (fe.max_dofs_per_cell());
if (n_components==1)
{
std::vector<double> rhs_values;
-
+
for (; cell!=endc; ++cell)
for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
if (cell->face(face)->at_boundary () &&
const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
n_q_points = fe_values.n_quadrature_points;
rhs_values.resize (n_q_points);
-
+
const std::vector<double> &weights = fe_values.get_JxW_values ();
rhs_function.value_list (fe_values.get_quadrature_points(), rhs_values);
-
+
cell_vector = 0;
for (unsigned int point=0; point<n_q_points; ++point)
- for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
cell_vector(i) += rhs_values[point] *
fe_values.shape_value(i,point) *
weights[point];
-
+
cell->get_dof_indices (dofs);
-
+
for (unsigned int i=0; i<dofs_per_cell; ++i)
rhs_vector(dofs[i]) += cell_vector(i);
}
else
{
std::vector<Vector<double> > rhs_values;
-
- for (; cell!=endc; ++cell)
+
+ for (; cell!=endc; ++cell)
for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
if (cell->face(face)->at_boundary () &&
(boundary_indicators.find (cell->face(face)->boundary_indicator())
const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
n_q_points = fe_values.n_quadrature_points;
rhs_values.resize (n_q_points, Vector<double>(n_components));
-
+
const std::vector<double> &weights = fe_values.get_JxW_values ();
rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values);
-
+
cell_vector = 0;
-
+
// Use the faster code if the
// FiniteElement is primitive
if (cell->get_fe().is_primitive ())
- {
+ {
for (unsigned int point=0; point<n_q_points; ++point)
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const unsigned int component
= cell->get_fe().system_to_component_index(i).first;
-
+
cell_vector(i) += rhs_values[point](component) *
fe_values.shape_value(i,point) *
weights[point];
weights[point];
}
}
-
+
cell->get_dof_indices (dofs);
-
+
for (unsigned int i=0; i<dofs_per_cell; ++i)
rhs_vector(dofs[i]) += cell_vector(i);
}
ExcMessage ("The number of components in the mask has to be either "
"zero or equal to the number of components in the finite "
"element."));
-
+
// check whether boundary values at
// the left or right boundary of
// the line are
// right.
const unsigned int direction = boundary_component;
Assert (direction < 2, ExcInvalidBoundaryIndicator());
-
+
// first find the outermost active
// cell by first traversing the coarse
// grid to its end and then going
typename DH::cell_iterator outermost_cell = dof.begin(0);
while (outermost_cell->neighbor(direction).state() == IteratorState::valid)
outermost_cell = outermost_cell->neighbor(direction);
-
+
while (outermost_cell->has_children())
outermost_cell = outermost_cell->child(direction);
else
boundary_function.vector_value (outermost_cell->vertex(direction),
function_values);
-
+
for (unsigned int i=0; i<fe.dofs_per_vertex; ++i)
if (component_mask[fe.face_system_to_component_index(i).first])
boundary_values[outermost_cell->vertex_dof_index(direction,i)]
// immediately
if (function_map.size() == 0)
return;
-
+
Assert (function_map.find(255) == function_map.end(),
ExcInvalidBoundaryIndicator());
for (unsigned int f=0; f<finite_elements.size(); ++f)
{
const FiniteElement<dim> &fe = finite_elements[f];
-
+
// generate a quadrature rule
// on the face from the unit
// support points. this will be
if (component_mask[fe.face_system_to_component_index(i).first]
== true)
unit_support_points[i] = fe.unit_face_support_point(i);
-
+
q_collection.push_back (Quadrature<dim-1>(unit_support_points));
- }
+ }
}
// now that we have a q_collection
// object with all the right
hp::MappingCollection<dim> mapping_collection (mapping);
hp::FEFaceValues<dim> x_fe_values (mapping_collection, finite_elements, q_collection,
update_quadrature_points);
-
+
typename DH::active_cell_iterator cell = dof.begin_active(),
endc = dof.end();
for (; cell!=endc; ++cell)
"values that correspond to primitive (scalar) base "
"elements"));
}
-
+
typename DH::face_iterator face = cell->face(face_no);
const unsigned char boundary_component = face->boundary_indicator();
- if (function_map.find(boundary_component) != function_map.end())
+ if (function_map.find(boundary_component) != function_map.end())
{
// face is of the right component
x_fe_values.reinit(cell, face_no);
face->get_dof_indices (face_dofs, cell->active_fe_index());
const std::vector<Point<DH::space_dimension> > &dof_locations
= fe_values.get_quadrature_points ();
-
+
if (fe_is_system)
{
// resize
Vector<double>(fe.n_components()));
else
dof_values_system.resize (fe.dofs_per_face);
-
+
function_map.find(boundary_component)->second
->vector_value_list (dof_locations, dof_values_system);
-
+
// enter those dofs
// into the list that
// match the
-
fe.get_nonzero_components(cell_i).begin());
}
-
+
if (component_mask[component] == true)
boundary_values[face_dofs[i]] = dof_values_system[i](component);
- }
+ }
}
else
// fe has only one component,
dof_values_scalar.resize (fe.dofs_per_face);
function_map.find(boundary_component)->second
->value_list (dof_locations, dof_values_scalar, 0);
-
+
// enter into list
-
+
for (unsigned int i=0; i<face_dofs.size(); ++i)
boundary_values[face_dofs[i]] = dof_values_scalar[i];
}
}
#endif
-
+
template <class DH>
void
const std::vector<bool> &component_mask)
{
Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
- interpolate_boundary_values(StaticMappingQ1<DH::dimension,DH::space_dimension>::mapping,
+ interpolate_boundary_values(StaticMappingQ1<DH::dimension,DH::space_dimension>::mapping,
dof, boundary_component,
boundary_function, boundary_values, component_mask);
}
const std::vector<bool> &component_mask)
{
Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
- interpolate_boundary_values(StaticMappingQ1<DH::dimension,DH::space_dimension>::mapping,
+ interpolate_boundary_values(StaticMappingQ1<DH::dimension,DH::space_dimension>::mapping,
dof, function_map,
boundary_values, component_mask);
}
ExcMessage ("The number of components in the mask has to be either "
"zero or equal to the number of components in the finite "
"element."));
-
+
// check whether boundary values at
// the left or right boundary of
// the line are
// right.
const unsigned int direction = boundary_component;
Assert (direction < 2, ExcInvalidBoundaryIndicator());
-
+
// first find the outermost active
// cell by first traversing the coarse
// grid to its end and then going
typename DH::cell_iterator outermost_cell = dof.begin(0);
while (outermost_cell->neighbor(direction).state() == IteratorState::valid)
outermost_cell = outermost_cell->neighbor(direction);
-
+
while (outermost_cell->has_children())
outermost_cell = outermost_cell->child(direction);
else
boundary_function.vector_value (outermost_cell->vertex(direction),
function_values);
-
+
for (unsigned int i=0; i<fe.dofs_per_vertex; ++i)
if (component_mask[fe.face_system_to_component_index(i).first])
{
// entries in the line here?
const unsigned int row = outermost_cell->vertex_dof_index(direction,i);
constraints.add_line (row);
- constraints.set_inhomogeneity (row,
+ constraints.set_inhomogeneity (row,
function_values(fe.face_system_to_component_index(i).first));
}
}
// Implementation for 1D
template <class DH>
void
-VectorTools::interpolate_boundary_values
+VectorTools::interpolate_boundary_values
(const Mapping<DH::dimension, DH::space_dimension> &mapping,
const DH &dof,
const typename FunctionMap<DH::space_dimension>::type &function_map,
template <class DH>
void
-VectorTools::interpolate_boundary_values
+VectorTools::interpolate_boundary_values
(const Mapping<DH::dimension, DH::space_dimension> &mapping,
const DH &dof,
const typename FunctionMap<DH::space_dimension>::type &function_map,
// immediately
if (function_map.size() == 0)
return;
-
+
Assert (function_map.find(255) == function_map.end(),
ExcInvalidBoundaryIndicator());
for (unsigned int f=0; f<finite_elements.size(); ++f)
{
const FiniteElement<dim> &fe = finite_elements[f];
-
+
// generate a quadrature rule on the
// face from the unit support
// points. this will be used to
if (component_mask[fe.face_system_to_component_index(i).first]
== true)
unit_support_points[i] = fe.unit_face_support_point(i);
-
+
q_collection.push_back (Quadrature<dim-1>(unit_support_points));
- }
+ }
}
// now that we have a q_collection object
// with all the right quadrature points,
hp::MappingCollection<dim> mapping_collection (mapping);
hp::FEFaceValues<dim> x_fe_values (mapping_collection, finite_elements, q_collection,
update_quadrature_points);
-
+
typename DH::active_cell_iterator cell = dof.begin_active(),
endc = dof.end();
for (; cell!=endc; ++cell)
"values that correspond to primitive (scalar) base "
"elements"));
}
-
+
typename DH::face_iterator face = cell->face(face_no);
const unsigned char boundary_component = face->boundary_indicator();
- if (function_map.find(boundary_component) != function_map.end())
+ if (function_map.find(boundary_component) != function_map.end())
{
// face is of the right
// component
face->get_dof_indices (face_dofs, cell->active_fe_index());
const std::vector<Point<DH::space_dimension> > &dof_locations
= fe_values.get_quadrature_points ();
-
+
if (fe_is_system)
{
// resize array. avoid
Vector<double>(fe.n_components()));
else
dof_values_system.resize (fe.dofs_per_face);
-
+
function_map.find(boundary_component)->second
->vector_value_list (dof_locations, dof_values_system);
-
+
// enter those dofs into
// the list that match the
// component
-
fe.get_nonzero_components(cell_i).begin());
}
-
+
if (component_mask[component] == true)
{
// TODO: check whether we should clear
constraints.set_inhomogeneity (face_dofs[i],
dof_values_system[i](component));
}
- }
+ }
}
else
// fe has only one component,
dof_values_scalar.resize (fe.dofs_per_face);
function_map.find(boundary_component)->second
->value_list (dof_locations, dof_values_scalar, 0);
-
+
// enter into list
-
+
for (unsigned int i=0; i<face_dofs.size(); ++i)
// TODO: check whether we should clear
// the current line first...
template <class DH>
void
-VectorTools::interpolate_boundary_values
+VectorTools::interpolate_boundary_values
(const Mapping<DH::dimension, DH::space_dimension> &mapping,
const DH &dof,
const unsigned char boundary_component,
}
#endif
-
+
template <class DH>
void
-VectorTools::interpolate_boundary_values
+VectorTools::interpolate_boundary_values
(const DH &dof,
const unsigned char boundary_component,
const Function<DH::space_dimension> &boundary_function,
const std::vector<bool> &component_mask)
{
Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
- interpolate_boundary_values(StaticMappingQ1<DH::dimension,DH::space_dimension>::mapping,
+ interpolate_boundary_values(StaticMappingQ1<DH::dimension,DH::space_dimension>::mapping,
dof, boundary_component,
boundary_function, constraints, component_mask);
}
template <class DH>
void
-VectorTools::interpolate_boundary_values
+VectorTools::interpolate_boundary_values
(const DH &dof,
const typename FunctionMap<DH::space_dimension>::type &function_map,
ConstraintMatrix &constraints,
const std::vector<bool> &component_mask)
{
Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
- interpolate_boundary_values(StaticMappingQ1<DH::dimension,DH::space_dimension>::mapping,
+ interpolate_boundary_values(StaticMappingQ1<DH::dimension,DH::space_dimension>::mapping,
dof, function_map,
constraints, component_mask);
}
}
else
AssertDimension (dof.get_fe().n_components(), component_mapping.size());
-
+
std::vector<unsigned int> dof_to_boundary_mapping;
std::set<unsigned char> selected_boundary_components;
for (typename FunctionMap<spacedim>::type::const_iterator i=boundary_functions.begin();
i!=boundary_functions.end(); ++i)
selected_boundary_components.insert (i->first);
-
+
DoFTools::map_dof_to_boundary_indices (dof, selected_boundary_components,
dof_to_boundary_mapping);
-
+
// Done if no degrees of freedom on
// the boundary
if (dof.n_boundary_dofs (boundary_functions) == 0)
if (dim>=3)
{
#ifdef DEBUG
-// Assert that there are no hanging nodes at the boundary
+// Assert that there are no hanging nodes at the boundary
int level = -1;
for (typename DoFHandler<dim,spacedim>::active_cell_iterator cell = dof.begin_active();
cell != dof.end(); ++cell)
if (level == -1)
level = cell->level();
else
- {
+ {
Assert (level == cell->level(), ExcNotImplemented());
}
}
#endif
}
sparsity.compress();
-
+
// make mass matrix and right hand side
SparseMatrix<double> mass_matrix(sparsity);
Vector<double> rhs(sparsity.n_rows());
- MatrixCreator::create_boundary_mass_matrix (mapping, dof, q,
+ MatrixCreator::create_boundary_mass_matrix (mapping, dof, q,
mass_matrix, boundary_functions,
rhs, dof_to_boundary_mapping, (const Function<spacedim>*) 0,
component_mapping);
// eliminate them here.
//TODO: Maybe we should figure out if the element really needs this
-
+
FilteredMatrix<Vector<double> > filtered_mass_matrix(mass_matrix, true);
FilteredMatrix<Vector<double> > filtered_precondition;
std::vector<bool> excluded_dofs(mass_matrix.m(), false);
for (unsigned int i=0;i<mass_matrix.m();++i)
if (mass_matrix.diag_element(i) > max_element)
max_element = mass_matrix.diag_element(i);
-
+
for (unsigned int i=0;i<mass_matrix.m();++i)
if (mass_matrix.diag_element(i) < 1.e-8 * max_element)
{
mass_matrix.diag_element(i) = 1.;
excluded_dofs[i] = true;
}
-
+
Vector<double> boundary_projection (rhs.size());
// Allow for a maximum of 5*n
// solve
cg.solve (filtered_mass_matrix, boundary_projection, rhs, filtered_precondition);
filtered_precondition.apply_constraints(boundary_projection, true);
- filtered_precondition.clear();
+ filtered_precondition.clear();
// fill in boundary values
for (unsigned int i=0; i<dof_to_boundary_mapping.size(); ++i)
if (dof_to_boundary_mapping[i] != DoFHandler<dim,spacedim>::invalid_dof_index
}
else
AssertDimension (dof.get_fe().n_components(), component_mapping.size());
-
+
std::vector<unsigned int> dof_to_boundary_mapping;
std::set<unsigned char> selected_boundary_components;
for (typename FunctionMap<spacedim>::type::const_iterator i=boundary_functions.begin();
i!=boundary_functions.end(); ++i)
selected_boundary_components.insert (i->first);
-
+
DoFTools::map_dof_to_boundary_indices (dof, selected_boundary_components,
dof_to_boundary_mapping);
-
+
// Done if no degrees of freedom on
// the boundary
if (dof.n_boundary_dofs (boundary_functions) == 0)
if (dim>=3)
{
#ifdef DEBUG
-// Assert that there are no hanging nodes at the boundary
+// Assert that there are no hanging nodes at the boundary
int level = -1;
for (typename DoFHandler<dim,spacedim>::active_cell_iterator cell = dof.begin_active();
cell != dof.end(); ++cell)
if (level == -1)
level = cell->level();
else
- {
+ {
Assert (level == cell->level(), ExcNotImplemented());
}
}
#endif
}
sparsity.compress();
-
+
// make mass matrix and right hand side
SparseMatrix<double> mass_matrix(sparsity);
Vector<double> rhs(sparsity.n_rows());
- MatrixCreator::create_boundary_mass_matrix (mapping, dof, q,
+ MatrixCreator::create_boundary_mass_matrix (mapping, dof, q,
mass_matrix, boundary_functions,
rhs, dof_to_boundary_mapping, (const Function<spacedim>*) 0,
component_mapping);
// eliminate them here.
//TODO: Maybe we should figure out if the element really needs this
-
+
FilteredMatrix<Vector<double> > filtered_mass_matrix(mass_matrix, true);
FilteredMatrix<Vector<double> > filtered_precondition;
std::vector<bool> excluded_dofs(mass_matrix.m(), false);
for (unsigned int i=0;i<mass_matrix.m();++i)
if (mass_matrix.diag_element(i) > max_element)
max_element = mass_matrix.diag_element(i);
-
+
for (unsigned int i=0;i<mass_matrix.m();++i)
if (mass_matrix.diag_element(i) < 1.e-8 * max_element)
{
mass_matrix.diag_element(i) = 1.;
excluded_dofs[i] = true;
}
-
+
Vector<double> boundary_projection (rhs.size());
// Allow for a maximum of 5*n
// solve
cg.solve (filtered_mass_matrix, boundary_projection, rhs, filtered_precondition);
filtered_precondition.apply_constraints(boundary_projection, true);
- filtered_precondition.clear();
+ filtered_precondition.clear();
// fill in boundary values
for (unsigned int i=0; i<dof_to_boundary_mapping.size(); ++i)
if (dof_to_boundary_mapping[i] != DoFHandler<dim,spacedim>::invalid_dof_index
for (unsigned int i=0; i<dim; ++i)
dof_indices[i] = numbers::invalid_unsigned_int;
}
-
-
+
+
bool operator < (const VectorDoFTuple<dim> &other) const
{
for (unsigned int i=0; i<dim; ++i)
add_constraint (const VectorDoFTuple<dim> &dof_indices,
const Tensor<1,dim> &constraining_vector,
ConstraintMatrix &constraints)
- {
+ {
// choose the DoF that has the
// largest component in the
// constraining_vector as the
(std::fabs(constraining_vector[0]) >= std::fabs(constraining_vector[2])))
{
constraints.add_line (dof_indices.dof_indices[0]);
-
+
if (std::fabs (constraining_vector[1]/constraining_vector[0])
> std::numeric_limits<double>::epsilon())
constraints.add_entry (dof_indices.dof_indices[0],
dof_indices.dof_indices[1],
-constraining_vector[1]/constraining_vector[2]);
}
-
+
break;
}
default:
Assert (false, ExcNotImplemented());
}
-
+
}
void
compute_orthonormal_vectors (const Tensor<1,dim> &vector,
Tensor<1,dim> (&orthonormals)[dim-1])
- {
+ {
switch (dim)
{
case 3:
default:
Assert (false, ExcNotImplemented());
- }
+ }
}
-
+
}
}
ExcMessage ("This function is not useful in 1d because it amounts "
"to imposing Dirichlet values on the vector-valued "
"quantity."));
-
+
const FiniteElement<dim> &fe = dof_handler.get_fe();
std::vector<unsigned int> face_dofs (fe.dofs_per_face);
std::vector<Point<spacedim> > dof_locations (fe.dofs_per_face);
- // have a map that stores normal vectors
- // for each vector-dof tuple we want to
- // constrain. since we can get at the same
- // vector dof tuple more than once (for
- // example if it is located at a vertex
- // that we visit from all adjacent cells),
- // we will want to average later on the
- // normal vectors computed on different
- // cells as described in the documentation
- // of this function. however, we can only
- // average if the contributions came from
- // different cells, whereas we want to
- // constrain twice or more in case the
- // contributions came from different faces
- // of the same cell. consequently, we also
- // have to store which cell a normal vector
- // was computed on
- typedef
+ // have a map that stores normal
+ // vectors for each vector-dof
+ // tuple we want to
+ // constrain. since we can get at
+ // the same vector dof tuple more
+ // than once (for example if it is
+ // located at a vertex that we
+ // visit from all adjacent cells),
+ // we will want to average later on
+ // the normal vectors computed on
+ // different cells as described in
+ // the documentation of this
+ // function. however, we can only
+ // average if the contributions
+ // came from different cells,
+ // whereas we want to constrain
+ // twice or more in case the
+ // contributions came from
+ // different faces of the same cell
+ // (i.e. constrain not just the
+ // *average normal direction* but
+ // *all normal directions* we
+ // find). consequently, we also
+ // have to store which cell a
+ // normal vector was computed on
+ typedef
std::multimap<internal::VectorTools::VectorDoFTuple<dim>,
std::pair<Tensor<1,dim>, typename DH<dim,spacedim>::active_cell_iterator> >
DoFToNormalsMap;
// vectors at the locations
// where they are defined:
const std::vector<Point<dim-1> > &
- unit_support_points = fe.get_unit_face_support_points();
+ unit_support_points = fe.get_unit_face_support_points();
Quadrature<dim-1> aux_quad (unit_support_points);
FEFaceValues<dim> fe_values (mapping, fe, aux_quad,
update_normal_vectors);
// components of vector
internal::VectorTools::VectorDoFTuple<dim> vector_dofs;
vector_dofs.dof_indices[0] = face_dofs[i];
-
+
for (unsigned int k=0; k<fe.dofs_per_face; ++k)
if ((k != i)
&&
vector_dofs.dof_indices[fe.face_system_to_component_index(k).first -
first_vector_component]
= face_dofs[k];
-
+
for (unsigned int d=0; d<dim; ++d)
Assert (vector_dofs.dof_indices[d] < dof_handler.n_dofs(),
ExcInternalError());
// dofs
typename DoFToNormalsMap::const_iterator
p = dof_to_normals_map.begin();
-
+
while (p != dof_to_normals_map.end())
{
// first find the range of entries in
// average over all normal vectors
case 1:
{
-
+
// compute the average
// normal vector from all
// the ones that have the
x = cell_to_normals_map.begin();
x != cell_to_normals_map.end(); ++x)
normal += x->second.first;
- normal /= normal.norm();
-
+ normal /= normal.norm();
+
// then construct constraints
// from this:
const internal::VectorTools::VectorDoFTuple<dim> &
- dof_indices = same_dof_range[0]->first;
+ dof_indices = same_dof_range[0]->first;
internal::VectorTools::add_constraint (dof_indices, normal,
constraints);
for (unsigned int i=0; i<dim; ++i, ++x)
for (unsigned int j=0; j<dim; ++j)
t[i][j] = x->second.first[j];
-
+
Assert (std::fabs(determinant (t)) > 1e-3,
ExcMessage("Found a set of normal vectors that are nearly collinear."));
}
-
+
// so all components of
// this vector dof are
// constrained. enter
break;
}
-
+
// this is the case of an
// edge contribution in 3d,
q != same_dof_range[1]; ++q)
cell_contributions[q->second.second].push_back (q->second.first);
Assert (cell_contributions.size() >= 1, ExcInternalError());
-
+
// now for each cell that
// has contributed
// determine the number
Assert (std::fabs (tangent.norm()-1) < 1e-12,
ExcInternalError());
-
+
tangential_vectors.push_back (tangent);
}
if (*t * first_tangent < 0)
*t *= -1;
}
-
+
// now compute the
// average tangent and
// normalize it
// necessary, so have a read-write
// version of it:
double exponent = exponent_1;
-
+
const unsigned int n_components = dof.get_fe().n_components();
const bool fe_is_system = (n_components != 1);
}
difference.reinit (dof.get_tria().n_active_cells());
-
+
switch (norm)
{
case dealii::VectorTools::L2_norm:
default:
break;
}
-
+
UpdateFlags update_flags = UpdateFlags (update_quadrature_points |
update_JxW_values);
switch (norm)
default:
update_flags |= UpdateFlags (update_values);
break;
- }
+ }
dealii::hp::FECollection<dim,spacedim> fe_collection (dof.get_fe());
dealii::hp::FEValues<dim,spacedim> x_fe_values(mapping, fe_collection, q, update_flags);
const unsigned int max_n_q_points = q.max_n_quadrature_points ();
-
+
std::vector< dealii::Vector<double> >
function_values (max_n_q_points, dealii::Vector<double>(n_components));
std::vector<std::vector<Tensor<1,spacedim> > >
// scalar functions
std::vector<double> tmp_values (max_n_q_points);
std::vector<Tensor<1,spacedim> > tmp_gradients (max_n_q_points);
-
+
// loop over all cells
typename DH::active_cell_iterator cell = dof.begin_active(),
endc = dof.end();
const dealii::FEValues<dim, spacedim> &fe_values = x_fe_values.get_present_fe_values ();
const unsigned int n_q_points = fe_values.n_quadrature_points;
-
+
// resize all out scratch
// arrays to the number of
// quadrature points we use
for (unsigned int k=0;k<n_q_points;++k)
weight_vectors[k] = 1.;
}
-
-
+
+
if (update_flags & update_values)
{
// first compute the exact solution
for (unsigned int i=0; i<n_q_points; ++i)
psi_values[i](0) = tmp_values[i];
}
-
+
// then subtract finite element
// fe_function
fe_values.get_function_values (fe_function, function_values);
tmp_gradients);
for (unsigned int i=0; i<n_q_points; ++i)
psi_grads[i][0] = tmp_gradients[i];
- }
-
+ }
+
// then subtract finite element
// function_grads
fe_values.get_function_grads (fe_function, function_grads);
for (unsigned int q=0; q<n_q_points; ++q)
psi_grads[q][k] -= function_grads[q][k];
}
-
+
switch (norm)
{
case dealii::VectorTools::mean:
psi_scalar[q] += std::pow(psi_values[q](k)*psi_values[q](k),
exponent/2.)
* weight_vectors[q](k);
-
+
// Integrate
diff = std::inner_product (psi_scalar.begin(), psi_scalar.end(),
fe_values.get_JxW_values().begin(),
psi_scalar[q] += std::pow(psi_grads[q][k] * psi_grads[q][k],
exponent/2.)
* weight_vectors[q](k);
-
+
diff += std::inner_product (psi_scalar.begin(), psi_scalar.end(),
fe_values.get_JxW_values().begin(),
0.0);
for (unsigned int d=0;d<dim;++d)
t = std::max(t,std::fabs(psi_grads[q][k][d])
* weight_vectors[q](k));
-
+
psi_scalar[q] = std::max(psi_scalar[q],t);
}
Assert(GeometryInfo<dim>::distance_to_unit_cell(cell_point.second) < 1e-10,
ExcInternalError());
-
+
const Quadrature<dim>
quadrature (GeometryInfo<dim>::project_to_unit_cell(cell_point.second));
FEValues<dim> fe_values(mapping, fe, quadrature, update_values);
difference(0) = exact_function.value(point);
else
exact_function.vector_value(point, difference);
-
+
for (unsigned int i=0; i<difference.size(); ++i)
difference(i) -= u_value[0](i);
}
Assert(GeometryInfo<dim>::distance_to_unit_cell(cell_point.second) < 1e-10,
ExcInternalError());
-
+
const Quadrature<dim>
quadrature (GeometryInfo<dim>::project_to_unit_cell(cell_point.second));
FEValues<dim> fe_values(mapping, fe, quadrature, update_values);
ExcDimensionMismatch (v.size(), dof.n_dofs()));
Assert (component < dof.get_fe().n_components(),
ExcIndexRange(component, 0, dof.get_fe().n_components()));
-
+
FEValues<dim> fe(mapping, dof.get_fe(), quadrature,
UpdateFlags(update_JxW_values
| update_values));
typename DoFHandler<dim,spacedim>::active_cell_iterator c;
std::vector<Vector<double> > values(quadrature.size(),
Vector<double> (dof.get_fe().n_components()));
-
+
double mean = 0.;
double area = 0.;
// Compute mean value
area += fe.JxW(k);
};
};
-
+
return (mean/area);
}
projects / dealii-svn.git / commitdiff