From: guido Date: Tue, 16 Feb 1999 13:40:14 +0000 (+0000) Subject: Most functions in FESystem implemented, the rest may be shouldnt X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c8d631f42563a2dffa07a824e81ea3497dc838f3;p=dealii-svn.git Most functions in FESystem implemented, the rest may be shouldnt git-svn-id: https://svn.dealii.org/trunk@813 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 a7ac0c2930..22ae825564 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -435,6 +435,11 @@ class FiniteElementBase : * Map between linear dofs and component dofs. */ vector< pair > system_to_component_table; + + /** + * Map between component and linear dofs. + */ + vector< vector > component_to_system_table; }; @@ -1130,11 +1135,13 @@ class FiniteElement : public FiniteElementBase { template inline unsigned -FiniteElementBase::component_to_system_index (unsigned /*component*/, +FiniteElementBase::component_to_system_index (unsigned component, unsigned component_index) const { - Assert(false, ExcNotImplemented()); - return component_index; + Assert(component diff --git a/deal.II/deal.II/include/fe/fe_system.h b/deal.II/deal.II/include/fe/fe_system.h index 5a5beaabf8..d702635683 100644 --- a/deal.II/deal.II/include/fe/fe_system.h +++ b/deal.II/deal.II/include/fe/fe_system.h @@ -101,7 +101,7 @@ class FESystem : public FiniteElement /** * Constructor for mixed * discretizations with two - * components. + * base elements. * * See the other constructor. */ @@ -317,7 +317,7 @@ class FESystem : public FiniteElement vector > &normal_vectors) const; /** - *Number of different composing + * Number of different base * elements of this object. * * Since these objects can have @@ -327,7 +327,7 @@ class FESystem : public FiniteElement * of finite elements composed * into this structure. */ - unsigned n_component_elements() const; + unsigned n_base_elements() const; /** * How often is a composing element used. @@ -372,6 +372,21 @@ class FESystem : public FiniteElement */ vector< ElementPair > base_elements; + /** + * The base element establishing a + * component. + * + * This table converts a + * component number to the + * #base_element# number. While + * component information contains + * multiplicity of base elements, + * the result allows access to + * shape functions of the base + * element. + * + */ + vector component_to_base_table; /** * Helper function used in the constructor: @@ -382,8 +397,7 @@ class FESystem : public FiniteElement * multiplied by #n#. Don't touch the * number of functions for the * transformation from unit to real - * cell. - */ + * cell. */ static FiniteElementData multiply_dof_numbers (const FiniteElementData &fe_data, const unsigned int N); @@ -421,7 +435,7 @@ class FESystem : public FiniteElement template inline unsigned -FESystem::n_component_elements() const +FESystem::n_base_elements() const { return base_elements.size(); } @@ -447,7 +461,8 @@ template template FESystem::FESystem (const FE &fe, const unsigned int n_elements) : FiniteElement (multiply_dof_numbers(fe, n_elements)), - base_elements(1) + base_elements(1), + component_to_base_table(n_components) { base_elements[0] = ElementPair(new FE, n_elements); base_elements[0].first -> subscribe (); @@ -460,7 +475,8 @@ FESystem::FESystem (const FE1 &fe1, const unsigned int n1, const FE2 &fe2, const unsigned int n2) : FiniteElement (multiply_dof_numbers(fe1, n1, fe2, n2)), - base_elements(2) + base_elements(2), + component_to_base_table(n_components) { base_elements[0] = ElementPair(new FE1, n1); base_elements[0].first -> subscribe (); diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index f8e488b8fd..ce0a9290e5 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -125,7 +125,8 @@ bool FiniteElementData::operator== (const FiniteElementData<2> &f) const template <> FiniteElementBase<1>::FiniteElementBase (const FiniteElementData<1> &fe_data) : FiniteElementData<1> (fe_data), - system_to_component_table(total_dofs) + system_to_component_table(total_dofs), + component_to_system_table(n_components, vector(total_dofs)) { const unsigned int dim=1; for (unsigned int i=0; i::children_per_cell; ++i) @@ -137,7 +138,10 @@ FiniteElementBase<1>::FiniteElementBase (const FiniteElementData<1> &fe_data) : interface_constraints(0,0)=1.; for (unsigned int j=0 ; j(0,j); + { + system_to_component_table[j] = pair(0,j); + component_to_system_table[0][j] = j; + } }; #endif @@ -148,7 +152,8 @@ FiniteElementBase<1>::FiniteElementBase (const FiniteElementData<1> &fe_data) : template <> FiniteElementBase<2>::FiniteElementBase (const FiniteElementData<2> &fe_data) : FiniteElementData<2> (fe_data), - system_to_component_table(total_dofs) + system_to_component_table(total_dofs), + component_to_system_table(n_components, vector(total_dofs)) { const unsigned int dim=2; for (unsigned int i=0; i::children_per_cell; ++i) @@ -159,7 +164,10 @@ FiniteElementBase<2>::FiniteElementBase (const FiniteElementData<2> &fe_data) : interface_constraints.reinit (dofs_per_vertex+2*dofs_per_line, 2*dofs_per_vertex+dofs_per_line); for (unsigned int j=0 ; j(0,j); + { + system_to_component_table[j] = pair(0,j); + component_to_system_table[0][j] = j; + } }; #endif diff --git a/deal.II/deal.II/source/fe/fe_system.cc b/deal.II/deal.II/source/fe/fe_system.cc index fc4d87be23..94991e57d4 100644 --- a/deal.II/deal.II/source/fe/fe_system.cc +++ b/deal.II/deal.II/source/fe/fe_system.cc @@ -23,30 +23,39 @@ FESystem::~FESystem () template void FESystem::initialize () { + // 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) + component_to_base_table[total_index++] = base; + // Initialize index table - // Multi-component base elements have to be thought of. + // Multi-component base elements have + // to be thought of. + // 1. Vertices - unsigned total_index = 0; + total_index = 0; for (unsigned vertex_number= 0 ; vertex_number < GeometryInfo<2>::vertices_per_cell ; ++vertex_number) { unsigned comp_start = 0; - for(unsigned comp = 0; comp < n_component_elements() ; - ++comp) + for(unsigned base = 0; base < n_base_elements() ; + ++base) { - for (unsigned m = 0; m < element_multiplicity(comp); ++m) + for (unsigned m = 0; m < element_multiplicity(base); ++m) { for (unsigned local_index = 0 ; - local_index < base_element(comp).dofs_per_vertex ; + local_index < base_element(base).dofs_per_vertex ; ++local_index) { system_to_component_table[total_index++] = pair (comp_start+m, - vertex_number*base_element(comp).dofs_per_vertex+local_index); + vertex_number*base_element(base).dofs_per_vertex+local_index); } } - comp_start += element_multiplicity(comp); + comp_start += element_multiplicity(base); } } @@ -55,23 +64,23 @@ void FESystem::initialize () ++line_number) { unsigned comp_start = 0; - for(unsigned comp = 0; comp < n_component_elements() ; - ++comp) + for(unsigned base = 0; base < n_base_elements() ; + ++base) { - for (unsigned m = 0; m < element_multiplicity(comp); ++m) + for (unsigned m = 0; m < element_multiplicity(base); ++m) { for (unsigned local_index = 0 ; - local_index < base_element(comp).dofs_per_line ; + local_index < base_element(base).dofs_per_line ; ++local_index) { system_to_component_table[total_index++] = pair (comp_start+m, - line_number*base_element(comp).dofs_per_line - +local_index+base_element(comp).first_line_index); + line_number*base_element(base).dofs_per_line + +local_index+base_element(base).first_line_index); } } - comp_start += element_multiplicity(comp); + comp_start += element_multiplicity(base); } } @@ -80,23 +89,23 @@ void FESystem::initialize () ++quad_number) { unsigned comp_start = 0; - for(unsigned comp = 0; comp < n_component_elements() ; - ++comp) + for(unsigned base = 0; base < n_base_elements() ; + ++base) { - for (unsigned m = 0; m < element_multiplicity(comp); ++m) + for (unsigned m = 0; m < element_multiplicity(base); ++m) { for (unsigned local_index = 0 ; - local_index < base_element(comp).dofs_per_quad ; + local_index < base_element(base).dofs_per_quad ; ++local_index) { system_to_component_table[total_index++] = pair (comp_start+m, - quad_number*base_element(comp).dofs_per_quad - +local_index+base_element(comp).first_quad_index); + quad_number*base_element(base).dofs_per_quad + +local_index+base_element(base).first_quad_index); } } - comp_start += element_multiplicity(comp); + comp_start += element_multiplicity(base); } } @@ -105,25 +114,37 @@ void FESystem::initialize () ++hex_number) { unsigned comp_start = 0; - for(unsigned comp = 0; comp < n_component_elements() ; - ++comp, comp_start += element_multiplicity(comp)) + for(unsigned base = 0; base < n_base_elements() ; + ++base) { - for (unsigned m = 0; m < element_multiplicity(comp); ++m) + for (unsigned m = 0; m < element_multiplicity(base); ++m) { for (unsigned local_index = 0 ; - local_index < base_element(comp).dofs_per_hex ; + local_index < base_element(base).dofs_per_hex ; ++local_index) { system_to_component_table[total_index++] = pair (comp_start+m, - hex_number*base_element(comp).dofs_per_hex - +local_index+base_element(comp).first_hex_index); + hex_number*base_element(base).dofs_per_hex + +local_index+base_element(base).first_hex_index); } } + comp_start += element_multiplicity(base); + } } + // Initialize mapping from components to + // linear index. Fortunately, this is + // the inverse of what we just did. + for (unsigned comp=0 ; comp::multiply_dof_numbers (const FiniteElementData<2> &fe1, template -double FESystem::shape_value (const unsigned int i, - const Point &p) const +double +FESystem::shape_value (const unsigned int i, + const Point &p) const { - Assert(false, ExcNotImplemented()); - return 0.; - Assert((ishape_value (i / n_sub_elements, p); + pair comp = system_to_component_index(i); + + return base_element(component_to_base_table[comp.first]) + .shape_value(comp.second, p); }; @@ -233,12 +254,12 @@ Tensor<1,dim> FESystem::shape_grad (const unsigned int i, const Point &p) const { - Assert(false, ExcNotImplemented()); - return Tensor<1,dim>(); - Assert((ishape_grad (i / n_sub_elements, p); + pair comp = system_to_component_index(i); + + return base_element(component_to_base_table[comp.first]) + .shape_grad(comp.second, p); }; @@ -248,20 +269,23 @@ Tensor<2,dim> FESystem::shape_grad_grad (const unsigned int i, const Point &p) const { - Assert(false, ExcNotImplemented()); - return Tensor<2,dim>(); - Assert((ishape_grad_grad (i / n_sub_elements, p); + + pair comp = system_to_component_index(i); + + return base_element(component_to_base_table[comp.first]) + .shape_grad_grad(comp.second, p); }; template -void FESystem::get_unit_support_points (vector > &support_points) const +void FESystem::get_unit_support_points (vector > &/*support_points*/) const { -/* Assert (support_points.size() == total_dofs, + Assert(false, ExcNotImplemented()); +/* + Assert (support_points.size() == total_dofs, ExcWrongFieldDimension (support_points.size(), total_dofs)); vector > base_support_points (base_element->total_dofs); @@ -276,9 +300,9 @@ void FESystem::get_unit_support_points (vector > &support_points template -void FESystem::get_support_points (const DoFHandler::cell_iterator &cell, - const Boundary &boundary, - vector > &support_points) const +void FESystem::get_support_points (const DoFHandler::cell_iterator &/*cell*/, + const Boundary &/*boundary*/, + vector > &/*support_points*/) const { Assert(false, ExcNotImplemented()); /* @@ -297,9 +321,9 @@ void FESystem::get_support_points (const DoFHandler::cell_iterator &ce 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()); /* @@ -318,9 +342,9 @@ void FESystem::get_face_support_points (const DoFHandler::face_iterato template -void FESystem::get_local_mass_matrix (const DoFHandler::cell_iterator &cell, - const Boundary &boundary, - dFMatrix &local_mass_matrix) const +void FESystem::get_local_mass_matrix (const DoFHandler::cell_iterator &/*cell*/, + const Boundary &/*boundary*/, + dFMatrix &/*local_mass_matrix*/) const { Assert(false, ExcNotImplemented()); /* @@ -350,19 +374,17 @@ void FESystem::get_local_mass_matrix (const DoFHandler::cell_iterator template -double FESystem::shape_value_transform (const unsigned int i, - const Point &p) const +double FESystem::shape_value_transform (const unsigned int /*i*/, + const Point &/*p*/) const { Assert(false, ExcNotImplemented()); return 0.; - -// return base_element->shape_value_transform (i, p); }; template -Tensor<1,dim> FESystem::shape_grad_transform (const unsigned int i, - const Point &p) const +Tensor<1,dim> FESystem::shape_grad_transform (const unsigned int /*i*/, + const Point &/*p*/) const { Assert(false, ExcNotImplemented()); return Tensor<1,dim>(); @@ -373,10 +395,10 @@ Tensor<1,dim> FESystem::shape_grad_transform (const unsigned int i, template -void FESystem::get_face_jacobians (const DoFHandler::face_iterator &face, - const Boundary &boundary, - const vector > &unit_points, - vector &face_jacobi_determinants) const +void FESystem::get_face_jacobians (const DoFHandler::face_iterator &/*face*/, + const Boundary &/*boundary*/, + const vector > &/*unit_points*/, + vector &/*face_jacobi_determinants*/) const { Assert(false, ExcNotImplemented()); @@ -386,10 +408,10 @@ void FESystem::get_face_jacobians (const DoFHandler::face_iterator &fa template -void FESystem::get_subface_jacobians (const DoFHandler::face_iterator &face, - const unsigned int subface_no, - const vector > &unit_points, - vector &face_jacobi_determinants) const +void FESystem::get_subface_jacobians (const DoFHandler::face_iterator &/*face*/, + const unsigned int /*subface_no*/, + const vector > &/*unit_points*/, + vector &/*face_jacobi_determinants*/) const { Assert(false, ExcNotImplemented()); @@ -400,11 +422,11 @@ void FESystem::get_subface_jacobians (const DoFHandler::face_iterator template -void FESystem::get_normal_vectors (const DoFHandler::cell_iterator &cell, - const unsigned int face_no, - const Boundary &boundary, - const vector > &unit_points, - vector > &normal_vectors) const +void FESystem::get_normal_vectors (const DoFHandler::cell_iterator &/*cell*/, + const unsigned int /*face_no*/, + const Boundary &/*boundary*/, + const vector > &/*unit_points*/, + vector > &/*normal_vectors*/) const { Assert(false, ExcNotImplemented()); @@ -414,11 +436,11 @@ void FESystem::get_normal_vectors (const DoFHandler::cell_iterator &ce template -void FESystem::get_normal_vectors (const DoFHandler::cell_iterator &cell, - const unsigned int face_no, - const unsigned int subface_no, - const vector > &unit_points, - vector > &normal_vectors) const +void FESystem::get_normal_vectors (const DoFHandler::cell_iterator &/*cell*/, + const unsigned int /*face_no*/, + const unsigned int /*subface_no*/, + const vector > &/*unit_points*/, + vector > &/*normal_vectors*/) const { Assert(false, ExcNotImplemented()); diff --git a/deal.II/lac/include/lac/dsmatrix.h b/deal.II/lac/include/lac/dsmatrix.h index e08e21ccd4..30484ce44a 100644 --- a/deal.II/lac/include/lac/dsmatrix.h +++ b/deal.II/lac/include/lac/dsmatrix.h @@ -475,7 +475,8 @@ class dSMatrix /** * Set the element #(i,j)# to #value#. * Throws an error if the entry does - * not exist. + * not exist. Still, it is allowed to store + * zero values in non-existent fields. */ void set (const unsigned int i, const unsigned int j, const double value); @@ -483,7 +484,8 @@ class dSMatrix /** * Add #value# to the element #(i,j)#. * Throws an error if the entry does - * not exist. + * not exist. Still, it is allowed to store + * zero values in non-existent fields. */ void add (const unsigned int i, const unsigned int j, const double value); @@ -784,9 +786,12 @@ inline void dSMatrix::set (const unsigned int i, const unsigned int j, const double value) { Assert (cols != 0, ExcMatrixNotInitialized()); - Assert (cols->operator()(i,j) != -1, + Assert ((cols->operator()(i,j) != -1) || (value == 0.), ExcInvalidIndex(i,j)); - val[cols->operator()(i,j)] = value; + + const int index = cols->operator()(i,j); + + if (index >= 0) val[index] = value; }; @@ -795,9 +800,12 @@ inline void dSMatrix::add (const unsigned int i, const unsigned int j, const double value) { Assert (cols != 0, ExcMatrixNotInitialized()); - Assert (cols->operator()(i,j) != -1, + Assert ((cols->operator()(i,j) != -1) || (value == 0.), ExcInvalidIndex(i,j)); - val[cols->operator()(i,j)] += value; + + const int index = cols->operator()(i,j); + + if (index >= 0) val[index] += value; };