*/
const unsigned int first_hex_index;
+ /**
+ * First index of dof on a line for face data.
+ */
+ const unsigned int first_face_line_index;
+
+ /**
+ * First index of dof on a quad for face data.
+ */
+ const unsigned int first_face_quad_index;
+
/**
* Total number of degrees of freedom
* on a cell. This information is
/**
* Compute system index from components.
*/
- unsigned component_to_system_index (unsigned component,
- unsigned component_index) const;
+ unsigned int component_to_system_index (unsigned int component,
+ unsigned int component_index) const;
/**
* Compute component and index from system index.
* component and second index in
* component.
*/
- pair<unsigned,unsigned> system_to_component_index (unsigned index) const;
+ pair<unsigned int,unsigned int> system_to_component_index (unsigned int index) const;
+
+ /**
+ * Compute system index from components on a face.
+ */
+ unsigned int face_component_to_system_index (unsigned int component,
+ unsigned int component_index) const;
+
+ /**
+ * Compute component and index from system
+ * index for a face.
+ *
+ * Return value contains first
+ * component and second index in
+ * component.
+ */
+ pair<unsigned int,unsigned int> face_system_to_component_index (unsigned int index) const;
/**
/**
* Map between linear dofs and component dofs.
*/
- vector< pair<unsigned, unsigned> > system_to_component_table;
+ vector< pair<unsigned int, unsigned int> > system_to_component_table;
+
+ /**
+ * Map between linear dofs and component dofs on face.
+ */
+ vector< pair<unsigned int, unsigned int> > face_system_to_component_table;
/**
* Map between component and linear dofs.
*/
- vector< vector<unsigned> > component_to_system_table;
+ vector< vector<unsigned int> > component_to_system_table;
+
+ /**
+ * Map between component and linear dofs on a face.
+ */
+ vector< vector<unsigned int> > face_component_to_system_table;
};
* components of mixed
* discretizations.
*/
- virtual const FiniteElement<dim>& base_element(unsigned index) const;
+ virtual const FiniteElement<dim>& base_element(unsigned int index) const;
/**
* Exception
};
template <int dim>
-inline unsigned
-FiniteElementBase<dim>::component_to_system_index (unsigned component,
- unsigned component_index) const
+inline unsigned int
+FiniteElementBase<dim>::component_to_system_index (unsigned int component,
+ unsigned int component_index) const
{
Assert(component<n_components, ExcInvalidIndex(component));
Assert(component_index<component_to_system_table[component].size(),
}
template <int dim>
-inline pair<unsigned,unsigned>
-FiniteElementBase<dim>::system_to_component_index (unsigned index) const
+inline pair<unsigned int,unsigned int>
+FiniteElementBase<dim>::system_to_component_index (unsigned int index) const
{
Assert(index < system_to_component_table.size(), ExcInvalidIndex(index));
return system_to_component_table[index];
}
+template <int dim>
+inline unsigned int
+FiniteElementBase<dim>::face_component_to_system_index (unsigned int component,
+ unsigned int component_index) const
+{
+ Assert(component<n_components, ExcInvalidIndex(component));
+ Assert(component_index<face_component_to_system_table[component].size(),
+ ExcInvalidIndex(component_index));
+ return face_component_to_system_table[component][component_index];
+}
+
+template <int dim>
+inline pair<unsigned int,unsigned int>
+FiniteElementBase<dim>::face_system_to_component_index (unsigned int index) const
+{
+ Assert(index < system_to_component_table.size(), ExcInvalidIndex(index));
+ return face_system_to_component_table[index];
+}
+
/*---------------------------- fe.h ---------------------------*/
/* end of #ifndef __fe_H */
#endif
* the #.h# file.
*/
void initialize();
+
+ /**
+ * Used by #initialize#.
+ */
+ void build_cell_table();
+
+ /**
+ * Used by #initialize#.
+ */
+ void build_face_table();
+
/**
*Exception.
*/
* @author Wolfgang Bangerth, 1998
*/
template <int dim>
-class VectorTools {
+class VectorTools //<dim>
+{
public:
/**
* Declare a data type which denotes a
*/
typedef map<unsigned char,const Function<dim>*> FunctionMap;
+ /**
+ * Data type for vector valued boundary function map.
+ */
+ typedef map<unsigned char,const VectorFunction<dim>*> VectorFunctionMap;
+
+// typedef map<pair<unsigned char, unsigned int> ,const Function<dim>*>
+// VectorFunctionMap;
+
/**
* Compute the interpolation of
* #function# at the support points to
Vector<double> &rhs_vector);
/**
- * Make up the list of node subject
+ * Make up the list of nodes subject
* to Dirichlet boundary conditions
- * and the values they are to be
- * assigned, by interpolation around
+ * and the values to be
+ * assigned to them, by interpolation around
* the boundary. If the
* #boundary_values# contained values
* before, the new ones are added, or
const FunctionMap &dirichlet_bc,
const Boundary<dim> &boundary,
map<int,double> &boundary_values);
-
+
+ /**
+ * Create boundary value information for vector
+ * valued functions.
+ * See the other #interpolate_boundary_values#.
+ */
+ static void interpolate_boundary_values (const DoFHandler<dim> &dof,
+ const VectorFunctionMap &dirichlet_bc,
+ const Boundary<dim> &boundary,
+ map<int,double> &boundary_values);
/**
* Project #function# to the boundary
* of the domain, using the given quadrature
dofs_per_face(GeometryInfo<dim>::vertices_per_face * dofs_per_vertex+
GeometryInfo<dim>::lines_per_face * dofs_per_line +
dofs_per_quad),
- first_line_index(GeometryInfo<dim>::vertices_per_cell * dofs_per_vertex),
+ first_line_index(GeometryInfo<dim>::vertices_per_cell
+ * dofs_per_vertex),
first_quad_index(first_line_index+
- 12*dofs_per_line),
+ GeometryInfo<dim>::lines_per_cell
+ * dofs_per_line),
first_hex_index(first_quad_index+
- GeometryInfo<dim>::faces_per_cell * dofs_per_quad),
+ GeometryInfo<dim>::faces_per_cell
+ * dofs_per_quad),
+ first_face_line_index(GeometryInfo<dim-1>::vertices_per_cell
+ * dofs_per_vertex),
+ first_face_quad_index(first_line_index+
+ GeometryInfo<dim-1>::lines_per_cell
+ * dofs_per_line),
total_dofs (first_hex_index+dofs_per_hex),
n_transform_functions (n_transform_functions),
n_components(n_components)
GeometryInfo<dim>::lines_per_cell * dofs_per_line),
first_hex_index(first_quad_index+
GeometryInfo<dim>::quads_per_cell*dofs_per_quad),
+ first_face_line_index(GeometryInfo<dim-1>::vertices_per_cell
+ * dofs_per_vertex),
+ first_face_quad_index(first_line_index+
+ GeometryInfo<dim-1>::lines_per_cell
+ * dofs_per_line),
total_dofs (first_quad_index+dofs_per_quad),
n_transform_functions (n_transform_functions),
n_components(n_components)
GeometryInfo<dim>::lines_per_cell * dofs_per_line),
first_hex_index(first_quad_index+
GeometryInfo<dim>::quads_per_cell*dofs_per_quad),
+ first_face_line_index(GeometryInfo<dim-1>::vertices_per_cell
+ * dofs_per_vertex),
+ first_face_quad_index(first_line_index+
+ GeometryInfo<dim-1>::lines_per_cell
+ * dofs_per_line),
total_dofs (first_line_index+dofs_per_line),
n_transform_functions (n_transform_functions),
n_components(n_components)
FiniteElementBase<dim>::FiniteElementBase (const FiniteElementData<dim> &fe_data) :
FiniteElementData<dim> (fe_data),
system_to_component_table(total_dofs),
- component_to_system_table(n_components, vector<unsigned>(total_dofs))
+ face_system_to_component_table(dofs_per_face),
+ component_to_system_table(n_components, vector<unsigned>(total_dofs)),
+ face_component_to_system_table(n_components, vector<unsigned>(dofs_per_face))
{
for (unsigned int i=0; i<GeometryInfo<dim>::children_per_cell; ++i)
{
system_to_component_table[j] = pair<unsigned,unsigned>(0,j);
component_to_system_table[0][j] = j;
}
+ for (unsigned int j=0 ; j<dofs_per_face ; ++j)
+ {
+ face_system_to_component_table[j] = pair<unsigned,unsigned>(0,j);
+ face_component_to_system_table[0][j] = j;
+ }
};
};
-
template <int dim>
-void FESystem<dim>::initialize ()
+void
+FESystem<dim>::build_cell_table()
{
- // Inintialize mapping from component
- // to base element
unsigned total_index = 0;
for (unsigned base=0 ; base < n_base_elements() ; ++base)
for (unsigned m = 0; m < element_multiplicity(base); ++m)
// to be thought of.
// 1. Vertices
- total_index = 0;
+ total_index = 0;
for (unsigned vertex_number= 0 ; vertex_number < GeometryInfo<dim>::vertices_per_cell ;
++vertex_number)
{
}
}
-
+ // Inintialize mapping from component
+ // to base element
// Initialize mapping from components to
// linear index. Fortunately, this is
// the inverse of what we just did.
for (unsigned sys=0 ; sys < total_dofs ; ++sys)
component_to_system_table[system_to_component_table[sys].first]
[system_to_component_table[sys].second] = sys;
+}
+
+
+template <int dim>
+void
+FESystem<dim>::build_face_table()
+{
+ unsigned total_index = 0;
+ for (unsigned base=0 ; base < n_base_elements() ; ++base)
+ for (unsigned m = 0; m < element_multiplicity(base); ++m)
+ component_to_base_table[total_index++] = base;
+
+ // Initialize index table
+ // Multi-component base elements have
+ // to be thought of.
+
+ // 1. Vertices
+ total_index = 0;
+ for (unsigned vertex_number= 0 ; vertex_number < GeometryInfo<dim>::vertices_per_face ;
+ ++vertex_number)
+ {
+ unsigned comp_start = 0;
+ for(unsigned base = 0; base < n_base_elements() ;
+ ++base)
+ {
+ for (unsigned m = 0; m < element_multiplicity(base); ++m)
+ {
+ for (unsigned local_index = 0 ;
+ local_index < base_element(base).dofs_per_vertex ;
+ ++local_index)
+ {
+ face_system_to_component_table[total_index++]
+ = pair<unsigned,unsigned>
+ (comp_start+m,
+ vertex_number*base_element(base).dofs_per_vertex+local_index);
+ }
+ }
+ comp_start += element_multiplicity(base);
+ }
+ }
+
+ // 2. Lines
+ for (unsigned line_number= 0 ; ((line_number != GeometryInfo<dim>::lines_per_face) &&
+ (GeometryInfo<dim>::lines_per_cell > 0));
+ ++line_number)
+ {
+ unsigned comp_start = 0;
+ for(unsigned base = 0; base < n_base_elements() ;
+ ++base)
+ {
+ for (unsigned m = 0; m < element_multiplicity(base); ++m)
+ {
+ for (unsigned local_index = 0 ;
+ local_index < base_element(base).dofs_per_line ;
+ ++local_index)
+ {
+ face_system_to_component_table[total_index++]
+ = pair<unsigned,unsigned>
+ (comp_start+m,
+ line_number*base_element(base).dofs_per_line
+ +local_index+base_element(base).first_face_line_index);
+ }
+ }
+ comp_start += element_multiplicity(base);
+ }
+ }
+
+ // 3. Quads
+ for (unsigned quad_number= 0 ; ((quad_number != GeometryInfo<dim>::quads_per_face) &&
+ (GeometryInfo<dim>::quads_per_cell > 0));
+ ++quad_number)
+ {
+ unsigned comp_start = 0;
+ for(unsigned base = 0; base < n_base_elements() ;
+ ++base)
+ {
+ for (unsigned m = 0; m < element_multiplicity(base); ++m)
+ {
+ for (unsigned local_index = 0 ;
+ local_index < base_element(base).dofs_per_quad ;
+ ++local_index)
+ {
+ face_system_to_component_table[total_index++]
+ = pair<unsigned,unsigned>
+ (comp_start+m,
+ quad_number*base_element(base).dofs_per_quad
+ +local_index+base_element(base).first_face_quad_index);
+ }
+ }
+ comp_start += element_multiplicity(base);
+ }
+ }
+ // Inintialize mapping from component
+ // to base element
+ // Initialize mapping from components to
+ // linear index. Fortunately, this is
+ // the inverse of what we just did.
+ for (unsigned comp=0 ; comp<n_components ; ++comp)
+ face_component_to_system_table[comp]
+ .resize(base_element(component_to_base_table[comp]).total_dofs);
+ for (unsigned sys=0 ; sys < dofs_per_face ; ++sys)
+ face_component_to_system_table[face_system_to_component_table[sys].first]
+ [face_system_to_component_table[sys].second] = sys;
+}
+
+
+
+template <int dim>
+void FESystem<dim>::initialize ()
+{
+ build_cell_table();
+ build_face_table();
// distribute the matrices of the base
// finite element to the matrices of
// this object
template <int dim>
-void FESystem<dim>::get_face_support_points (const DoFHandler<dim>::face_iterator &/*face*/,
- const Boundary<dim> &/*boundary*/,
- vector<Point<dim> > &/*support_points*/) const
+void FESystem<dim>::get_face_support_points (const DoFHandler<dim>::face_iterator & face,
+ const Boundary<dim> & boundary,
+ vector<Point<dim> > & support_points) const
{
- Assert(false, ExcNotImplemented());
-/*
Assert (support_points.size() == dofs_per_face,
ExcWrongFieldDimension (support_points.size(), dofs_per_face));
- vector<Point<dim> > base_support_points (base_element->dofs_per_face);
- base_element->get_face_support_points (face, boundary, base_support_points);
-
- for (unsigned int i=0; i<base_element->dofs_per_face; ++i)
- for (unsigned int n=0; n<n_sub_elements; ++n)
- support_points[i*n_sub_elements+n] = base_support_points[i];
-*/
-};
+ vector<Point<dim> > base_support_points (base_element(0).dofs_per_face);
+ unsigned int comp = 0;
+ for (unsigned int base=0 ; base<n_base_elements(); ++base)
+ {
+ base_support_points.resize(base_element(base).dofs_per_face);
+ base_element(base).get_face_support_points (face, boundary, base_support_points);
+ for (unsigned int inbase = 0 ; inbase < element_multiplicity(base); ++inbase)
+ {
+ for (unsigned int i=0; i<base_element(base).dofs_per_face; ++i)
+ {
+ support_points[face_component_to_system_index(comp,i)]
+ = base_support_points[i];
+ }
+
+ ++comp;
+ }
+ }
+}
// first get the local mass matrix for
// the base object
- FullMatrix<double> base_mass_matrix (base_element->total_dofs,
- base_element->total_dofs);
- base_element->get_local_mass_matrix (cell, boundary, base_mass_matrix);
+ FullMatrix<double> base_mass_matrix (base_element.total_dofs,
+ base_element.total_dofs);
+ base_element.get_local_mass_matrix (cell, boundary, base_mass_matrix);
// now distribute it to the mass matrix
// of this object
- for (unsigned int i=0; i<base_element->total_dofs; ++i)
- for (unsigned int j=0; j<base_element->total_dofs; ++j)
+ for (unsigned int i=0; i<base_element.total_dofs; ++i)
+ for (unsigned int j=0; j<base_element.total_dofs; ++j)
for (unsigned int n=0; n<n_sub_elements; ++n)
// only fill diagonals of the blocks
local_mass_matrix (i*n_sub_elements + n,
VectorTools<1>::interpolate_boundary_values (const DoFHandler<1> &,
const FunctionMap &,
const Boundary<1> &,
- map<int,double> &) {
+ map<int,double> &)
+{
+ Assert (false, ExcNotImplemented());
+};
+
+template <>
+void VectorTools<1>::interpolate_boundary_values (const DoFHandler<1> &,
+ const VectorFunctionMap&,
+ const Boundary<1>&,
+ map<int,double>&)
+{
Assert (false, ExcNotImplemented());
};
};
};
+template <int dim>
+void
+VectorTools<dim>::interpolate_boundary_values (const DoFHandler<dim> &dof,
+ const VectorFunctionMap &dirichlet_bc,
+ const Boundary<dim> &boundary,
+ map<int,double> &boundary_values)
+{
+ Assert (dirichlet_bc.find(255) == dirichlet_bc.end(),
+ ExcInvalidBoundaryIndicator());
+
+ const FiniteElement<dim> &fe = dof.get_fe();
+
+ // use two face iterators, since we need
+ // a DoF-iterator for the dof indices, but
+ // a Tria-iterator for the fe object
+ DoFHandler<dim>::active_face_iterator face = dof.begin_active_face(),
+ endf = dof.end_face();
+
+ typename VectorFunctionMap::const_iterator function_ptr;
+
+ // field to store the indices of dofs
+ vector<int> face_dofs (fe.dofs_per_face);
+ vector<Point<dim> > dof_locations (face_dofs.size(), Point<dim>());
+ vector< Vector<double> > dof_values (fe.dofs_per_face, Vector<double>(fe.n_components));
+
+ for (; face!=endf; ++face)
+ if ((function_ptr = dirichlet_bc.find(face->boundary_indicator())) !=
+ dirichlet_bc.end())
+ // face is subject to one of the
+ // bc listed in #dirichlet_bc#
+ {
+ // get indices, physical location and
+ // boundary values of dofs on this
+ // face
+ face->get_dof_indices (face_dofs);
+ fe.get_face_support_points (face, boundary, dof_locations);
+ function_ptr->second->value_list (dof_locations, dof_values);
+
+ // enter into list
+
+ for (unsigned int i=0; i<face_dofs.size(); ++i)
+ {
+ pair<unsigned int, unsigned int>
+ index = fe.face_system_to_component_index(i);
+
+ boundary_values[face_dofs[i]] = dof_values[i](index.first);
+ }
+ }
+}
+
template <int dim>