From: guido Date: Fri, 12 Mar 1999 09:16:20 +0000 (+0000) Subject: mixed discretization with Dirichlet boundary works X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=406a368d06bf450a1588bd026ed408941e0a5be0;p=dealii-svn.git mixed discretization with Dirichlet boundary works git-svn-id: https://svn.dealii.org/trunk@997 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index bc8c05745a..459711ffe9 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -70,6 +70,16 @@ class FiniteElementData */ 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 @@ -239,8 +249,8 @@ class FiniteElementBase : /** * 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. @@ -249,7 +259,23 @@ class FiniteElementBase : * component and second index in * component. */ - pair system_to_component_index (unsigned index) const; + pair 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 face_system_to_component_index (unsigned int index) const; /** @@ -430,12 +456,22 @@ class FiniteElementBase : /** * Map between linear dofs and component dofs. */ - vector< pair > system_to_component_table; + vector< pair > system_to_component_table; + + /** + * Map between linear dofs and component dofs on face. + */ + vector< pair > face_system_to_component_table; /** * Map between component and linear dofs. */ - vector< vector > component_to_system_table; + vector< vector > component_to_system_table; + + /** + * Map between component and linear dofs on a face. + */ + vector< vector > face_component_to_system_table; }; @@ -1333,7 +1369,7 @@ class FiniteElement : public FiniteElementBase * components of mixed * discretizations. */ - virtual const FiniteElement& base_element(unsigned index) const; + virtual const FiniteElement& base_element(unsigned int index) const; /** * Exception @@ -1362,9 +1398,9 @@ class FiniteElement : public FiniteElementBase }; template -inline unsigned -FiniteElementBase::component_to_system_index (unsigned component, - unsigned component_index) const +inline unsigned int +FiniteElementBase::component_to_system_index (unsigned int component, + unsigned int component_index) const { Assert(component::component_to_system_index (unsigned component, } template -inline pair -FiniteElementBase::system_to_component_index (unsigned index) const +inline pair +FiniteElementBase::system_to_component_index (unsigned int index) const { Assert(index < system_to_component_table.size(), ExcInvalidIndex(index)); return system_to_component_table[index]; } +template +inline unsigned int +FiniteElementBase::face_component_to_system_index (unsigned int component, + unsigned int component_index) const +{ + Assert(component +inline pair +FiniteElementBase::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 diff --git a/deal.II/deal.II/include/fe/fe_system.h b/deal.II/deal.II/include/fe/fe_system.h index 19e51b3aaa..8d6d5d3859 100644 --- a/deal.II/deal.II/include/fe/fe_system.h +++ b/deal.II/deal.II/include/fe/fe_system.h @@ -445,6 +445,17 @@ class FESystem // * the #.h# file. */ void initialize(); + + /** + * Used by #initialize#. + */ + void build_cell_table(); + + /** + * Used by #initialize#. + */ + void build_face_table(); + /** *Exception. */ diff --git a/deal.II/deal.II/include/numerics/vectors.h b/deal.II/deal.II/include/numerics/vectors.h index 21e09386ba..2e1d5e432a 100644 --- a/deal.II/deal.II/include/numerics/vectors.h +++ b/deal.II/deal.II/include/numerics/vectors.h @@ -243,7 +243,8 @@ enum NormType { * @author Wolfgang Bangerth, 1998 */ template -class VectorTools { +class VectorTools // +{ public: /** * Declare a data type which denotes a @@ -259,6 +260,14 @@ class VectorTools { */ typedef map*> FunctionMap; + /** + * Data type for vector valued boundary function map. + */ + typedef map*> VectorFunctionMap; + +// typedef map ,const Function*> +// VectorFunctionMap; + /** * Compute the interpolation of * #function# at the support points to @@ -329,10 +338,10 @@ class VectorTools { Vector &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 @@ -347,7 +356,16 @@ class VectorTools { const FunctionMap &dirichlet_bc, const Boundary &boundary, map &boundary_values); - + + /** + * Create boundary value information for vector + * valued functions. + * See the other #interpolate_boundary_values#. + */ + static void interpolate_boundary_values (const DoFHandler &dof, + const VectorFunctionMap &dirichlet_bc, + const Boundary &boundary, + map &boundary_values); /** * Project #function# to the boundary * of the domain, using the given quadrature diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index f28668aeb9..f4dbb5fdde 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -36,11 +36,19 @@ FiniteElementData::FiniteElementData (const unsigned int dofs_per_vertex, dofs_per_face(GeometryInfo::vertices_per_face * dofs_per_vertex+ GeometryInfo::lines_per_face * dofs_per_line + dofs_per_quad), - first_line_index(GeometryInfo::vertices_per_cell * dofs_per_vertex), + first_line_index(GeometryInfo::vertices_per_cell + * dofs_per_vertex), first_quad_index(first_line_index+ - 12*dofs_per_line), + GeometryInfo::lines_per_cell + * dofs_per_line), first_hex_index(first_quad_index+ - GeometryInfo::faces_per_cell * dofs_per_quad), + GeometryInfo::faces_per_cell + * dofs_per_quad), + first_face_line_index(GeometryInfo::vertices_per_cell + * dofs_per_vertex), + first_face_quad_index(first_line_index+ + GeometryInfo::lines_per_cell + * dofs_per_line), total_dofs (first_hex_index+dofs_per_hex), n_transform_functions (n_transform_functions), n_components(n_components) @@ -67,6 +75,11 @@ FiniteElementData::FiniteElementData (const unsigned int dofs_per_vertex, GeometryInfo::lines_per_cell * dofs_per_line), first_hex_index(first_quad_index+ GeometryInfo::quads_per_cell*dofs_per_quad), + first_face_line_index(GeometryInfo::vertices_per_cell + * dofs_per_vertex), + first_face_quad_index(first_line_index+ + GeometryInfo::lines_per_cell + * dofs_per_line), total_dofs (first_quad_index+dofs_per_quad), n_transform_functions (n_transform_functions), n_components(n_components) @@ -91,6 +104,11 @@ FiniteElementData::FiniteElementData (const unsigned int dofs_per_vertex, GeometryInfo::lines_per_cell * dofs_per_line), first_hex_index(first_quad_index+ GeometryInfo::quads_per_cell*dofs_per_quad), + first_face_line_index(GeometryInfo::vertices_per_cell + * dofs_per_vertex), + first_face_quad_index(first_line_index+ + GeometryInfo::lines_per_cell + * dofs_per_line), total_dofs (first_line_index+dofs_per_line), n_transform_functions (n_transform_functions), n_components(n_components) @@ -120,7 +138,9 @@ template FiniteElementBase::FiniteElementBase (const FiniteElementData &fe_data) : FiniteElementData (fe_data), system_to_component_table(total_dofs), - component_to_system_table(n_components, vector(total_dofs)) + face_system_to_component_table(dofs_per_face), + component_to_system_table(n_components, vector(total_dofs)), + face_component_to_system_table(n_components, vector(dofs_per_face)) { for (unsigned int i=0; i::children_per_cell; ++i) { @@ -160,6 +180,11 @@ FiniteElementBase::FiniteElementBase (const FiniteElementData &fe_data system_to_component_table[j] = pair(0,j); component_to_system_table[0][j] = j; } + for (unsigned int j=0 ; j(0,j); + face_component_to_system_table[0][j] = j; + } }; diff --git a/deal.II/deal.II/source/fe/fe_system.cc b/deal.II/deal.II/source/fe/fe_system.cc index 3bd2efb204..c7cbb58a48 100644 --- a/deal.II/deal.II/source/fe/fe_system.cc +++ b/deal.II/deal.II/source/fe/fe_system.cc @@ -19,12 +19,10 @@ FESystem::~FESystem () }; - template -void FESystem::initialize () +void +FESystem::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) @@ -35,7 +33,7 @@ void FESystem::initialize () // to be thought of. // 1. Vertices - total_index = 0; + total_index = 0; for (unsigned vertex_number= 0 ; vertex_number < GeometryInfo::vertices_per_cell ; ++vertex_number) { @@ -137,7 +135,8 @@ void FESystem::initialize () } } - + // Inintialize mapping from component + // to base element // Initialize mapping from components to // linear index. Fortunately, this is // the inverse of what we just did. @@ -148,8 +147,120 @@ void FESystem::initialize () 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 +void +FESystem::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::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 + (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::lines_per_face) && + (GeometryInfo::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 + (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::quads_per_face) && + (GeometryInfo::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 + (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 +void FESystem::initialize () +{ + build_cell_table(); + build_face_table(); // distribute the matrices of the base // finite element to the matrices of // this object @@ -324,23 +435,31 @@ void FESystem::get_support_points (const DoFHandler::cell_iterator &/* template -void FESystem::get_face_support_points (const DoFHandler::face_iterator &/*face*/, - const Boundary &/*boundary*/, - vector > &/*support_points*/) const +void FESystem::get_face_support_points (const DoFHandler::face_iterator & face, + const Boundary & boundary, + vector > & support_points) const { - Assert(false, ExcNotImplemented()); -/* Assert (support_points.size() == dofs_per_face, ExcWrongFieldDimension (support_points.size(), dofs_per_face)); - vector > base_support_points (base_element->dofs_per_face); - base_element->get_face_support_points (face, boundary, base_support_points); - - for (unsigned int i=0; idofs_per_face; ++i) - for (unsigned int n=0; n > base_support_points (base_element(0).dofs_per_face); + unsigned int comp = 0; + for (unsigned int base=0 ; base::get_local_mass_matrix (const DoFHandler::cell_iterator // first get the local mass matrix for // the base object - FullMatrix base_mass_matrix (base_element->total_dofs, - base_element->total_dofs); - base_element->get_local_mass_matrix (cell, boundary, base_mass_matrix); + FullMatrix 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; itotal_dofs; ++i) - for (unsigned int j=0; jtotal_dofs; ++j) + for (unsigned int i=0; i::interpolate_boundary_values (const DoFHandler<1> &, const FunctionMap &, const Boundary<1> &, - map &) { + map &) +{ + Assert (false, ExcNotImplemented()); +}; + +template <> +void VectorTools<1>::interpolate_boundary_values (const DoFHandler<1> &, + const VectorFunctionMap&, + const Boundary<1>&, + map&) +{ Assert (false, ExcNotImplemented()); }; @@ -292,6 +302,56 @@ VectorTools::interpolate_boundary_values (const DoFHandler &dof, }; }; +template +void +VectorTools::interpolate_boundary_values (const DoFHandler &dof, + const VectorFunctionMap &dirichlet_bc, + const Boundary &boundary, + map &boundary_values) +{ + Assert (dirichlet_bc.find(255) == dirichlet_bc.end(), + ExcInvalidBoundaryIndicator()); + + const FiniteElement &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::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 face_dofs (fe.dofs_per_face); + vector > dof_locations (face_dofs.size(), Point()); + vector< Vector > dof_values (fe.dofs_per_face, Vector(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 + index = fe.face_system_to_component_index(i); + + boundary_values[face_dofs[i]] = dof_values[i](index.first); + } + } +} + template