From: hartmann Date: Tue, 10 Aug 1999 14:45:05 +0000 (+0000) Subject: use a boolean vector for the restriction_is_additive_flags to make it work also for... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b044057e050a4b0c09bf1c49c4ea1ff21e212773;p=dealii-svn.git use a boolean vector for the restriction_is_additive_flags to make it work also for mixed systems git-svn-id: https://svn.dealii.org/trunk@1659 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 c2390ec97b..7e8062fd4e 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -105,103 +105,6 @@ class FiniteElementData */ const unsigned int n_components; - /** - * This flag determines how the restriction - * of data from child cells to its mother - * is to be done. In this, it also - * determines in which way the restriction - * matrices of the derived class are to - * be used. - * - * For most elements, the mode is the - * following. Consider a 1d linear element, - * with two children and nodal values - * 1 and 2 on the first child, and 2 and 4 - * on the second child. The restriction - * to the mother child then yields the - * values 1 and four, i.e. the values on - * the mother cell can be obtained by - * pointwise interpolation, where for each - * nodal value on the mother child one - * point on exactly one child exists. - * However, already on the quadratic - * element, the midpoint on the mother - * element can be obtained from any of - * the two children, which however would - * both yield the same value due to - * continuity. What we do in practice - * is to compute them from both sides - * and set them, rather than add them up. - * This makes some things much easier. In - * practice, if a degree of freedom on - * one of the child cells yields a - * nonzero contribution to one of the - * degrees of freedom on the mother - * cell, we overwrite the value on - * the mother cell. This way, when setting - * up the restriction matrices, we do not - * have to track which child is responsible - * for setting a given value on the mother - * cell. We call this the non-additive - * mode. - * - * The other possibility would be to - * add up the contributions from the - * different children. This would mean - * that both of the inner endpoint of - * the quadratic child elements above - * would have a weight of 1/2 with - * respect to the midpoint value on - * the mother cell. However, this also - * means that we have to first compute - * the restriction to the mother cell - * by addition from the child cells, and - * afterwards set them to the global - * vector. The same process, adding - * up the local contributions to the - * global vector is not possible since - * we do not know how many coarse cells - * contribute to nodes on the boundary. - * - * In contrast to the non-additive mode - * described above, which is the simplest - * way for elements can be interpolated - * from its children, interpolation is - * not possible for piecewise constant - * elements, to name only one example. - * Here, the value on the mother cell - * has to be taken the average of the - * values on the children, i.e. all - * children contribute alike to the - * one degree of freedom. Here, we have - * to sum up the contributions of all - * child cells with the same weight, - * and the non-additive mode of above - * would only set the value on the mother - * cell to the value of one of the child - * cell, irrespective of the values on the - * other cells. - * - * Similarly, for discontinuous linear - * elements, it might be better to not - * interpolate the values at the corners - * from the child cells, but to take a - * better average, for example - * interpolating at the centers of the - * child cells; in that case, the - * contributions of the child cells - * have to be additive as well. - * - * Given these notes, the flag under - * consideration has to be set to #false# - * for the usual continuous Lagrange - * elements, and #true# for the other - * cases mentioned above. The main function - * where it is used is - * #DoFAccessor::get_interpolated_dof_values#. - */ - const bool restriction_is_additive; - /** * Default constructor. Constructs * an element @@ -218,8 +121,7 @@ class FiniteElementData FiniteElementData (const unsigned int dofs_per_vertex, const unsigned int dofs_per_line, const unsigned int n_transform_functions, - const unsigned int n_components, - const bool restriction_is_additive); + const unsigned int n_components); /** * Constructor for a 2-dimensional @@ -229,8 +131,7 @@ class FiniteElementData const unsigned int dofs_per_line, const unsigned int dofs_per_quad, const unsigned int n_transform_functions, - const unsigned int n_components, - const bool restriction_is_additive); + const unsigned int n_components); /** * Constructor for a 3-dimensional @@ -241,8 +142,7 @@ class FiniteElementData const unsigned int dofs_per_quad, const unsigned int dofs_per_hex, const unsigned int n_transform_functions, - const unsigned int n_components, - const bool restriction_is_additive); + const unsigned int n_components); /** * Declare this destructor virtual in @@ -300,7 +200,8 @@ class FiniteElementBase : public Subscriptor, * matrices explicitely after calling * this base class' constructor. */ - FiniteElementBase (const FiniteElementData &fe_data); + FiniteElementBase (const FiniteElementData &fe_data, + const vector &restriction_is_additive_flags); /** * Return a readonly reference to the @@ -392,7 +293,13 @@ class FiniteElementBase : public Subscriptor, * element. */ unsigned int component_to_base(unsigned int index) const; - + + /** + * Access the #restriction_is_additive_flag# + * field. See there for more information on + * its contents. + */ + bool restriction_is_additive (const unsigned int component) const; /** * Exception @@ -596,6 +503,104 @@ class FiniteElementBase : public Subscriptor, * element. */ vector component_to_base_table; + + /** + * This flag determines how the restriction + * of data from child cells to its mother + * is to be done. In this, it also + * determines in which way the restriction + * matrices of the derived class are to + * be used. + * + * For most elements, the mode is the + * following. Consider a 1d linear element, + * with two children and nodal values + * 1 and 2 on the first child, and 2 and 4 + * on the second child. The restriction + * to the mother child then yields the + * values 1 and four, i.e. the values on + * the mother cell can be obtained by + * pointwise interpolation, where for each + * nodal value on the mother child one + * point on exactly one child exists. + * However, already on the quadratic + * element, the midpoint on the mother + * element can be obtained from any of + * the two children, which however would + * both yield the same value due to + * continuity. What we do in practice + * is to compute them from both sides + * and set them, rather than add them up. + * This makes some things much easier. In + * practice, if a degree of freedom on + * one of the child cells yields a + * nonzero contribution to one of the + * degrees of freedom on the mother + * cell, we overwrite the value on + * the mother cell. This way, when setting + * up the restriction matrices, we do not + * have to track which child is responsible + * for setting a given value on the mother + * cell. We call this the non-additive + * mode. + * + * The other possibility would be to + * add up the contributions from the + * different children. This would mean + * that both of the inner endpoint of + * the quadratic child elements above + * would have a weight of 1/2 with + * respect to the midpoint value on + * the mother cell. However, this also + * means that we have to first compute + * the restriction to the mother cell + * by addition from the child cells, and + * afterwards set them to the global + * vector. The same process, adding + * up the local contributions to the + * global vector is not possible since + * we do not know how many coarse cells + * contribute to nodes on the boundary. + * + * In contrast to the non-additive mode + * described above, which is the simplest + * way for elements can be interpolated + * from its children, interpolation is + * not possible for piecewise constant + * elements, to name only one example. + * Here, the value on the mother cell + * has to be taken the average of the + * values on the children, i.e. all + * children contribute alike to the + * one degree of freedom. Here, we have + * to sum up the contributions of all + * child cells with the same weight, + * and the non-additive mode of above + * would only set the value on the mother + * cell to the value of one of the child + * cell, irrespective of the values on the + * other cells. + * + * Similarly, for discontinuous linear + * elements, it might be better to not + * interpolate the values at the corners + * from the child cells, but to take a + * better average, for example + * interpolating at the centers of the + * child cells; in that case, the + * contributions of the child cells + * have to be additive as well. + * + * Given these notes, the flag under + * consideration has to be set to #false# + * for the usual continuous Lagrange + * elements, and #true# for the other + * cases mentioned above. The main function + * where it is used is + * #DoFAccessor::get_interpolated_dof_values#. + * There is one flag per component. + */ + const vector restriction_is_additive_flags; }; @@ -897,7 +902,8 @@ class FiniteElement : public FiniteElementBase /** * Constructor */ - FiniteElement (const FiniteElementData &fe_data); + FiniteElement (const FiniteElementData &fe_data, + const vector restriction_is_additive_flags); /** * Destructor. Only declared to have a @@ -1619,6 +1625,17 @@ FiniteElementBase::component_to_base (unsigned int index) const return component_to_base_table[index]; } +template +inline +bool +FiniteElementBase::restriction_is_additive (const unsigned int component) const +{ + Assert(component const FiniteElementData &fe3, const unsigned int N3); + + /** + * Helper function used in the constructor: + * takes a #FiniteElement# object + * and returns an boolean vector including + * the #restriction_is_additive_flags# of + * the mixed element consisting of #N# + * elements of the sub-element #fe#. + */ + static vector + compute_restriction_is_additive_flags (const FiniteElement &fe, + const unsigned int N); + /** + * Same as above for mixed elements + * with two different sub-elements. + */ + static vector + compute_restriction_is_additive_flags (const FiniteElement &fe1, + const unsigned int N1, + const FiniteElement &fe2, + const unsigned int N2); + + /** + * Same as above for mixed elements + * with three different sub-elements. + */ + static vector + compute_restriction_is_additive_flags (const FiniteElement &fe1, + const unsigned int N1, + const FiniteElement &fe2, + const unsigned int N2, + const FiniteElement &fe3, + const unsigned int N3); + /** * This function is simply singled out of * the constructor. It sets up the @@ -494,32 +528,11 @@ FESystem::n_base_elements() const -template -inline unsigned int -FESystem::element_multiplicity(unsigned int index) const -{ - Assert (index < base_elements.size(), - ExcIndexRange(index, 0, base_elements.size())); - return base_elements[index].second; -} - - - -template -inline const FiniteElement& -FESystem::base_element(unsigned int index) const -{ - Assert (index < base_elements.size(), - ExcIndexRange(index, 0, base_elements.size())); - return *base_elements[index].first; -} - - - template template FESystem::FESystem (const FE &fe, const unsigned int n_elements) : - FiniteElement (multiply_dof_numbers(fe, n_elements)), + FiniteElement (multiply_dof_numbers(fe, n_elements), + compute_restriction_is_additive_flags (fe, n_elements)), base_elements(1) { base_elements[0] = ElementPair(new FE, n_elements); @@ -534,7 +547,9 @@ template FESystem::FESystem (const FE1 &fe1, const unsigned int n1, const FE2 &fe2, const unsigned int n2) : - FiniteElement (multiply_dof_numbers(fe1, n1, fe2, n2)), + FiniteElement (multiply_dof_numbers(fe1, n1, fe2, n2), + compute_restriction_is_additive_flags (fe1, n1, + fe2, n2)), base_elements(2) { Assert(fe1.n_transform_functions == fe2.n_transform_functions, @@ -557,7 +572,10 @@ FESystem::FESystem (const FE1 &fe1, const unsigned int n1, : FiniteElement (multiply_dof_numbers(fe1, n1, fe2, n2, - fe3, n3)), + fe3, n3), + compute_restriction_is_additive_flags (fe1, n1, + fe2, n2, + fe3, n3)), base_elements(3) { Assert(fe1.n_transform_functions == fe2.n_transform_functions, @@ -576,6 +594,31 @@ FESystem::FESystem (const FE1 &fe1, const unsigned int n1, +template +inline unsigned int +FESystem::element_multiplicity(unsigned int index) const +{ + Assert (index < base_elements.size(), + ExcIndexRange(index, 0, base_elements.size())); + return base_elements[index].second; +} + + + +template +inline const FiniteElement& +FESystem::base_element(unsigned int index) const +{ + Assert (index < base_elements.size(), + ExcIndexRange(index, 0, base_elements.size())); + return *base_elements[index].first; +} + + + + + + /*---------------------------- fe_lib.system.h ---------------------------*/ /* end of #ifndef __fe_system_H */ diff --git a/deal.II/deal.II/include/fe/q1_mapping.h b/deal.II/deal.II/include/fe/q1_mapping.h index e0638265a1..2b88ab4cb8 100644 --- a/deal.II/deal.II/include/fe/q1_mapping.h +++ b/deal.II/deal.II/include/fe/q1_mapping.h @@ -32,10 +32,10 @@ class FEQ1Mapping : public FiniteElement */ FEQ1Mapping (const unsigned int dofs_per_vertex, const unsigned int dofs_per_line, - const unsigned int dofs_per_quad =0, - const unsigned int dofs_per_hex =0, - const unsigned int n_components =1, - const bool restriction_is_additive=false); + const unsigned int dofs_per_quad, + const unsigned int dofs_per_hex, + const unsigned int n_components, + const vector restriction_is_additive_flags); /** * Return the value of the #i#th shape diff --git a/deal.II/deal.II/source/dofs/dof_accessor.cc b/deal.II/deal.II/source/dofs/dof_accessor.cc index bf8a635ce6..b39d07b7fb 100644 --- a/deal.II/deal.II/source/dofs/dof_accessor.cc +++ b/deal.II/deal.II/source/dofs/dof_accessor.cc @@ -680,9 +680,6 @@ DoFCellAccessor::get_interpolated_dof_values (const Vector &values, interpolated_values.clear (); - const bool restriction_is_additive - = dof_handler->get_fe().restriction_is_additive; - for (unsigned int child=0; child::children_per_cell; ++child) { @@ -702,18 +699,18 @@ DoFCellAccessor::get_interpolated_dof_values (const Vector &values, // end in adding up the contribution // from nodes on boundaries of // children more than once. - if (restriction_is_additive) + for (unsigned int i=0; iget_fe().system_to_component_index(i).first; + + if (dof_handler->get_fe().restriction_is_additive(component)) interpolated_values(i) += tmp2(i); - } - else - { - for (unsigned int i=0; i::FiniteElementData (const unsigned int dofs_per_vertex, const unsigned int dofs_per_quad, const unsigned int dofs_per_hex, const unsigned int n_transform_functions, - const unsigned int n_components, - const bool restriction_is_additive) : + const unsigned int n_components) : dofs_per_vertex(dofs_per_vertex), dofs_per_line(dofs_per_line), dofs_per_quad(dofs_per_quad), @@ -52,8 +51,7 @@ FiniteElementData::FiniteElementData (const unsigned int dofs_per_vertex, * dofs_per_line), total_dofs (first_hex_index+dofs_per_hex), n_transform_functions (n_transform_functions), - n_components(n_components), - restriction_is_additive(restriction_is_additive) + n_components(n_components) { Assert(dim==3, ExcDimensionMismatch(3,dim)); }; @@ -65,8 +63,7 @@ FiniteElementData::FiniteElementData (const unsigned int dofs_per_vertex, const unsigned int dofs_per_line, const unsigned int dofs_per_quad, const unsigned int n_transform_functions, - const unsigned int n_components, - const bool restriction_is_additive) : + const unsigned int n_components) : dofs_per_vertex(dofs_per_vertex), dofs_per_line(dofs_per_line), dofs_per_quad(dofs_per_quad), @@ -85,8 +82,7 @@ FiniteElementData::FiniteElementData (const unsigned int dofs_per_vertex, * dofs_per_line), total_dofs (first_quad_index+dofs_per_quad), n_transform_functions (n_transform_functions), - n_components(n_components), - restriction_is_additive(restriction_is_additive) + n_components(n_components) { Assert(dim==2, ExcDimensionMismatch(2,dim)); }; @@ -97,8 +93,7 @@ template FiniteElementData::FiniteElementData (const unsigned int dofs_per_vertex, const unsigned int dofs_per_line, const unsigned int n_transform_functions, - const unsigned int n_components, - const bool restriction_is_additive) : + const unsigned int n_components) : dofs_per_vertex(dofs_per_vertex), dofs_per_line(dofs_per_line), dofs_per_quad(0), @@ -116,8 +111,7 @@ FiniteElementData::FiniteElementData (const unsigned int dofs_per_vertex, * dofs_per_line), total_dofs (first_line_index+dofs_per_line), n_transform_functions (n_transform_functions), - n_components(n_components), - restriction_is_additive(restriction_is_additive) + n_components(n_components) { Assert(dim==1, ExcDimensionMismatch(1,dim)); }; @@ -139,8 +133,7 @@ bool FiniteElementData::operator== (const FiniteElementData &f) const (dofs_per_quad == f.dofs_per_quad) && (dofs_per_hex == f.dofs_per_hex) && (n_transform_functions == f.n_transform_functions) && - (n_components == f.n_components) && - (restriction_is_additive == f.restriction_is_additive)); + (n_components == f.n_components)); }; @@ -150,14 +143,19 @@ bool FiniteElementData::operator== (const FiniteElementData &f) const template -FiniteElementBase::FiniteElementBase (const FiniteElementData &fe_data) : +FiniteElementBase::FiniteElementBase (const FiniteElementData &fe_data, + const vector &restriction_is_additive_flags) : FiniteElementData (fe_data), system_to_component_table(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)), - component_to_base_table(n_components) + component_to_base_table(n_components), + restriction_is_additive_flags(restriction_is_additive_flags) { + Assert(restriction_is_additive_flags.size()==fe_data.n_components, + ExcWrongFieldDimension(restriction_is_additive_flags.size(),fe_data.n_components)); + for (unsigned int i=0; i::children_per_cell; ++i) { restriction[i].reinit (total_dofs, total_dofs); @@ -250,14 +248,14 @@ bool FiniteElementBase::operator == (const FiniteElementBase &f) const - - /*------------------------------- FiniteElement ----------------------*/ template -FiniteElement::FiniteElement (const FiniteElementData &fe_data) : - FiniteElementBase (fe_data) {}; +FiniteElement::FiniteElement (const FiniteElementData &fe_data, + const vector restriction_is_additive_flags) : + FiniteElementBase (fe_data, + restriction_is_additive_flags) {}; #if deal_II_dimension == 1 diff --git a/deal.II/deal.II/source/fe/fe_lib.criss_cross.cc b/deal.II/deal.II/source/fe/fe_lib.criss_cross.cc index b8fab956ff..53238f0763 100644 --- a/deal.II/deal.II/source/fe/fe_lib.criss_cross.cc +++ b/deal.II/deal.II/source/fe/fe_lib.criss_cross.cc @@ -310,7 +310,8 @@ template <> FECrissCross<1>::FECrissCross () : // set more or less invalid data - FiniteElement<1> (FiniteElementData<1> (0,0,0,0)) + FiniteElement<1> (FiniteElementData<1> (0,0,0,0), + vector()) { Assert (false, ExcNotUseful()); }; @@ -458,7 +459,8 @@ void FECrissCross<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &, template <> FECrissCross<2>::FECrissCross () : - FiniteElement<2> (FiniteElementData<2> (1,0,1,5,1)) + FiniteElement<2> (FiniteElementData<2> (1,0,1,5,1), + vector (1, false)) { interface_constraints(0,0) = 1./2.; interface_constraints(0,1) = 1./2.; diff --git a/deal.II/deal.II/source/fe/fe_lib.cubic.cc b/deal.II/deal.II/source/fe/fe_lib.cubic.cc index 9c852bc5f4..8a99f20dc9 100644 --- a/deal.II/deal.II/source/fe/fe_lib.cubic.cc +++ b/deal.II/deal.II/source/fe/fe_lib.cubic.cc @@ -21,7 +21,8 @@ template <> void FEQ3::initialize_matrices (); template <> FEQ3<1>::FEQ3 () : - FEQ1Mapping<1> (1, 2) { + FEQ1Mapping<1> (1, 2, 0, 0, 1, + vector (1, false)) { initialize_matrices (); }; @@ -29,7 +30,8 @@ FEQ3<1>::FEQ3 () : template <> FEQ3<1>::FEQ3 (const int) : - FEQ1Mapping<1> (0, 4) + FEQ1Mapping<1> (0, 4, 0, 0, 1, + vector (1, false)) { initialize_matrices (); }; @@ -240,7 +242,8 @@ void FEQ3<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterator &cell, template <> FEQ3<2>::FEQ3 () : - FEQ1Mapping<2> (1, 2, 4) + FEQ1Mapping<2> (1, 2, 4, 0, 1, + vector (1, false)) { interface_constraints(0,0) = -1.0/16.0; interface_constraints(0,1) = -1.0/16.0; @@ -264,7 +267,8 @@ FEQ3<2>::FEQ3 () : template <> FEQ3<2>::FEQ3 (const int) : - FEQ1Mapping<2> (0, 0, 16) + FEQ1Mapping<2> (0, 0, 16, 0, 1, + vector (1, false)) { initialize_matrices (); }; @@ -1603,7 +1607,8 @@ void FEQ3<2>::get_face_support_points (const DoFHandler<2>::face_iterator &face, template <> FEQ3<3>::FEQ3 () : - FEQ1Mapping<3> (1, 2, 4, 8) + FEQ1Mapping<3> (1, 2, 4, 8, 1, + vector (1, false)) { interface_constraints(0,0) = 1.0/256.0; interface_constraints(0,1) = 1.0/256.0; @@ -1865,7 +1870,8 @@ FEQ3<3>::FEQ3 () : template <> FEQ3<3>::FEQ3 (const int) : - FEQ1Mapping<3> (0, 0, 0, 64) + FEQ1Mapping<3> (0, 0, 0, 64, 1, + vector (1, false)) { initialize_matrices (); }; diff --git a/deal.II/deal.II/source/fe/fe_lib.dg.constant.cc b/deal.II/deal.II/source/fe/fe_lib.dg.constant.cc index 912999ee60..a9d8810202 100644 --- a/deal.II/deal.II/source/fe/fe_lib.dg.constant.cc +++ b/deal.II/deal.II/source/fe/fe_lib.dg.constant.cc @@ -18,7 +18,7 @@ FEDG_Q0::FEDG_Q0 () : (dim==2 ? 1 : 0), (dim==3 ? 1 : 0), 1, - true) + vector (1, true)) { for (unsigned int i=0; i::children_per_cell; ++i) { diff --git a/deal.II/deal.II/source/fe/fe_lib.linear.cc b/deal.II/deal.II/source/fe/fe_lib.linear.cc index aa22e94404..4778488d2e 100644 --- a/deal.II/deal.II/source/fe/fe_lib.linear.cc +++ b/deal.II/deal.II/source/fe/fe_lib.linear.cc @@ -22,7 +22,8 @@ template <> void FEQ1::initialize_matrices (); template <> FEQ1<1>::FEQ1 () : - FEQ1Mapping<1> (1, 0) + FEQ1Mapping<1> (1, 0, 0, 0, 1, + vector (1, false)) { initialize_matrices (); }; @@ -30,7 +31,8 @@ FEQ1<1>::FEQ1 () : template <> FEQ1<1>::FEQ1 (const int) : - FEQ1Mapping<1> (0, 2) + FEQ1Mapping<1> (0, 2, 0, 0, 1, + vector (1, false)) { initialize_matrices (); }; @@ -165,7 +167,8 @@ void FEQ1<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterator &cell, template <> FEQ1<2>::FEQ1 () : - FEQ1Mapping<2> (1, 0, 0) + FEQ1Mapping<2> (1, 0, 0, 0, 1, + vector (1, false)) { interface_constraints(0,0) = 1./2.; interface_constraints(0,1) = 1./2.; @@ -177,7 +180,8 @@ FEQ1<2>::FEQ1 () : template <> FEQ1<2>::FEQ1 (const int) : - FEQ1Mapping<2> (0, 0, 4) + FEQ1Mapping<2> (0, 0, 4, 0, 1, + vector (1, false)) { initialize_matrices (); }; @@ -433,7 +437,8 @@ void FEQ1<2>::get_unit_support_points (vector > &unit_points) const { template <> FEQ1<3>::FEQ1 () : - FEQ1Mapping<3> (1, 0, 0, 0) + FEQ1Mapping<3> (1, 0, 0, 0, 1, + vector (1, false)) { interface_constraints(0,0) = 1.0/4.0; interface_constraints(0,1) = 1.0/4.0; @@ -455,7 +460,8 @@ FEQ1<3>::FEQ1 () : template <> FEQ1<3>::FEQ1 (const int) : - FEQ1Mapping<3> (0, 0, 0, 8) + FEQ1Mapping<3> (0, 0, 0, 8, 1, + vector (1, false)) { initialize_matrices (); }; diff --git a/deal.II/deal.II/source/fe/fe_lib.quadratic.cc b/deal.II/deal.II/source/fe/fe_lib.quadratic.cc index 430eb71af2..0dd21c0f68 100644 --- a/deal.II/deal.II/source/fe/fe_lib.quadratic.cc +++ b/deal.II/deal.II/source/fe/fe_lib.quadratic.cc @@ -19,7 +19,8 @@ template <> void FEQ2::initialize_matrices (); template <> FEQ2<1>::FEQ2 () : - FEQ1Mapping<1> (1, 1) { + FEQ1Mapping<1> (1, 1, 0, 0, 1, + vector (1, false)) { initialize_matrices (); }; @@ -27,7 +28,8 @@ FEQ2<1>::FEQ2 () : template <> FEQ2<1>::FEQ2 (const int) : - FEQ1Mapping<1> (0, 3) { + FEQ1Mapping<1> (0, 3, 0, 0, 1, + vector (1, false)) { initialize_matrices (); }; @@ -121,7 +123,7 @@ void FEQ2<1>::initialize_matrices () { template <> double FEQ2<1>::shape_value(const unsigned int i, - const Point<1> &p) const + const Point<1> &p) const { Assert((i::shape_value(const unsigned int i, template <> Tensor<1,1> FEQ2<1>::shape_grad(const unsigned int i, - const Point<1> &p) const + const Point<1> &p) const { Assert((i::shape_grad(const unsigned int i, template <> Tensor<2,1> FEQ2<1>::shape_grad_grad (const unsigned int i, - const Point<1> &) const + const Point<1> &) const { Assert((i::get_face_support_points (const DoFHandler<1>::face_iterator &, template <> void FEQ2<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterator &cell, - FullMatrix &local_mass_matrix) const { + FullMatrix &local_mass_matrix) const { Assert (local_mass_matrix.n() == total_dofs, ExcWrongFieldDimension(local_mass_matrix.n(),total_dofs)); Assert (local_mass_matrix.m() == total_dofs, @@ -232,7 +234,8 @@ void FEQ2<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterator &cell, template <> FEQ2<2>::FEQ2 () : - FEQ1Mapping<2> (1, 1, 1) + FEQ1Mapping<2> (1, 1, 1, 0, 1, + vector (1, false)) { interface_constraints(0,2) = 1.0; interface_constraints(1,0) = 3./8.; @@ -249,7 +252,8 @@ FEQ2<2>::FEQ2 () : template <> FEQ2<2>::FEQ2 (const int) : - FEQ1Mapping<2> (0, 0, 9) + FEQ1Mapping<2> (0, 0, 9, 0, 1, + vector (1, false)) { initialize_matrices (); }; @@ -1121,7 +1125,8 @@ void FEQ2<2>::get_face_support_points (const DoFHandler<2>::face_iterator &face, template <> FEQ2<3>::FEQ2 () : - FEQ1Mapping<3> (1, 1, 1, 1) + FEQ1Mapping<3> (1, 1, 1, 1, 1, + vector (1, false)) { interface_constraints(0,8) = 1.0; interface_constraints(1,4) = 1.0; @@ -1208,7 +1213,8 @@ FEQ2<3>::FEQ2 () : template <> FEQ2<3>::FEQ2 (const int) : - FEQ1Mapping<3> (0, 0, 0, 27) + FEQ1Mapping<3> (0, 0, 0, 27, 1, + vector (1, false)) { initialize_matrices (); }; diff --git a/deal.II/deal.II/source/fe/fe_lib.quartic.cc b/deal.II/deal.II/source/fe/fe_lib.quartic.cc index 94e7f846aa..9ee9594842 100644 --- a/deal.II/deal.II/source/fe/fe_lib.quartic.cc +++ b/deal.II/deal.II/source/fe/fe_lib.quartic.cc @@ -19,14 +19,16 @@ template <> void FEQ4::initialize_matrices (); template <> FEQ4<1>::FEQ4 () : - FEQ1Mapping<1> (1, 3) { + FEQ1Mapping<1> (1, 3, 0, 0, 1, + vector (1, false)) { initialize_matrices (); }; template <> FEQ4<1>::FEQ4 (const int) : - FEQ1Mapping<1> (0, 5) { + FEQ1Mapping<1> (0, 5, 0, 0, 1, + vector (1, false)) { initialize_matrices (); }; @@ -218,7 +220,8 @@ void FEQ4<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterator &cell, template <> FEQ4<2>::FEQ4 () : - FEQ1Mapping<2> (1, 3, 9) + FEQ1Mapping<2> (1, 3, 9, 0, 1, + vector (1, false)) { interface_constraints(0,3) = 1.0; interface_constraints(1,0) = 35.0/128.0; @@ -251,7 +254,8 @@ FEQ4<2>::FEQ4 () : template <> FEQ4<2>::FEQ4 (const int) : - FEQ1Mapping<2> (0, 0, 25) + FEQ1Mapping<2> (0, 0, 25, 0, 1, + vector (1, false)) { initialize_matrices (); }; diff --git a/deal.II/deal.II/source/fe/fe_system.cc b/deal.II/deal.II/source/fe/fe_system.cc index 2ecc79a584..82bb713624 100644 --- a/deal.II/deal.II/source/fe/fe_system.cc +++ b/deal.II/deal.II/source/fe/fe_system.cc @@ -514,6 +514,62 @@ FESystem<2>::multiply_dof_numbers (const FiniteElementData<2> &fe1, #endif +template +vector +FESystem::compute_restriction_is_additive_flags (const FiniteElement &fe, + const unsigned int n_elements) +{ + vector tmp; + for (unsigned int i=0; i +vector +FESystem::compute_restriction_is_additive_flags (const FiniteElement &fe1, + const unsigned int N1, + const FiniteElement &fe2, + const unsigned int N2) +{ + vector tmp; + for (unsigned int i=0; i +vector +FESystem::compute_restriction_is_additive_flags (const FiniteElement &fe1, + const unsigned int N1, + const FiniteElement &fe2, + const unsigned int N2, + const FiniteElement &fe3, + const unsigned int N3) +{ + vector tmp; + for (unsigned int i=0; i double FESystem::shape_value (const unsigned int i, diff --git a/deal.II/deal.II/source/fe/q1_mapping.cc b/deal.II/deal.II/source/fe/q1_mapping.cc index 2753e1e8a7..a496850ce8 100644 --- a/deal.II/deal.II/source/fe/q1_mapping.cc +++ b/deal.II/deal.II/source/fe/q1_mapping.cc @@ -24,12 +24,12 @@ FEQ1Mapping<1>::FEQ1Mapping (const unsigned int dofs_per_vertex, const unsigned int dofs_per_quad, const unsigned int dofs_per_hex, const unsigned int n_components, - const bool restriction_is_additive) : + const vector restriction_is_additive_flags) : FiniteElement<1> (FiniteElementData<1> (dofs_per_vertex, dofs_per_line, GeometryInfo<1>::vertices_per_cell, - n_components, - restriction_is_additive)) + n_components), + restriction_is_additive_flags) { Assert (dofs_per_quad==0, ExcInvalidData()); Assert (dofs_per_hex==0, ExcInvalidData()); @@ -147,13 +147,13 @@ FEQ1Mapping<2>::FEQ1Mapping (const unsigned int dofs_per_vertex, const unsigned int dofs_per_quad, const unsigned int dofs_per_hex, const unsigned int n_components, - const bool restriction_is_additive) : + const vector restriction_is_additive_flags) : FiniteElement<2> (FiniteElementData<2> (dofs_per_vertex, dofs_per_line, dofs_per_quad, GeometryInfo<2>::vertices_per_cell, - n_components, - restriction_is_additive)) + n_components), + restriction_is_additive_flags) { Assert (dofs_per_hex == 0, ExcInvalidData()); }; @@ -321,14 +321,14 @@ FEQ1Mapping<3>::FEQ1Mapping (const unsigned int dofs_per_vertex, const unsigned int dofs_per_quad, const unsigned int dofs_per_hex, const unsigned int n_components, - const bool restriction_is_additive) : + const vector restriction_is_additive_flags) : FiniteElement<3> (FiniteElementData<3> (dofs_per_vertex, dofs_per_line, dofs_per_quad, dofs_per_hex, GeometryInfo<3>::vertices_per_cell, - n_components, - restriction_is_additive)) + n_components), + restriction_is_additive_flags) {};