From: hartmann Date: Mon, 9 Aug 1999 15:45:18 +0000 (+0000) Subject: add restriction is additive flag and change the constructors X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=85b71e55a222a3fd6ffde6159156638a32b326fb;p=dealii-svn.git add restriction is additive flag and change the constructors git-svn-id: https://svn.dealii.org/trunk@1647 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 39618580d9..c2390ec97b 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -100,15 +100,115 @@ class FiniteElementData /** - * Number of components and dimension of the image space. + * Number of components and dimension of + * the image space. */ 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 + * Default constructor. Constructs + * an element * which is not so useful. Checking * #total_dofs# is therefore a good way to - * check if something went wrong. */ + * check if something went wrong. + */ FiniteElementData (); /** @@ -118,7 +218,8 @@ 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 unsigned int n_components, + const bool restriction_is_additive); /** * Constructor for a 2-dimensional @@ -128,7 +229,8 @@ 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 unsigned int n_components, + const bool restriction_is_additive); /** * Constructor for a 3-dimensional @@ -139,7 +241,8 @@ 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 unsigned int n_components, + const bool restriction_is_additive); /** * Declare this destructor virtual in @@ -251,7 +354,8 @@ class FiniteElementBase : public Subscriptor, unsigned int component_index) const; /** - * Compute component and index from system index. + * Compute component and index from + * system index. * * Return value contains first * component and second index in diff --git a/deal.II/deal.II/include/fe/q1_mapping.h b/deal.II/deal.II/include/fe/q1_mapping.h index 60b9ff60c2..e0638265a1 100644 --- a/deal.II/deal.II/include/fe/q1_mapping.h +++ b/deal.II/deal.II/include/fe/q1_mapping.h @@ -19,9 +19,7 @@ * are implemented here and do not have to be taken care of later. */ template -class FEQ1Mapping - : - public FiniteElement +class FEQ1Mapping : public FiniteElement { public: /** @@ -33,10 +31,11 @@ class FEQ1Mapping * shall be zero. */ 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 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); /** * Return the value of the #i#th shape diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index ded20956be..a84a261f7e 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -28,7 +28,8 @@ FiniteElementData::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 unsigned int n_components, + const bool restriction_is_additive) : dofs_per_vertex(dofs_per_vertex), dofs_per_line(dofs_per_line), dofs_per_quad(dofs_per_quad), @@ -51,7 +52,8 @@ 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) + n_components(n_components), + restriction_is_additive(restriction_is_additive) { Assert(dim==3, ExcDimensionMismatch(3,dim)); }; @@ -63,7 +65,8 @@ 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 unsigned int n_components, + const bool restriction_is_additive) : dofs_per_vertex(dofs_per_vertex), dofs_per_line(dofs_per_line), dofs_per_quad(dofs_per_quad), @@ -82,7 +85,8 @@ 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) + n_components(n_components), + restriction_is_additive(restriction_is_additive) { Assert(dim==2, ExcDimensionMismatch(2,dim)); }; @@ -93,7 +97,8 @@ 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 unsigned int n_components, + const bool restriction_is_additive) : dofs_per_vertex(dofs_per_vertex), dofs_per_line(dofs_per_line), dofs_per_quad(0), @@ -111,7 +116,8 @@ 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) + n_components(n_components), + restriction_is_additive(restriction_is_additive) { Assert(dim==1, ExcDimensionMismatch(1,dim)); }; @@ -133,7 +139,8 @@ 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)); + (n_components == f.n_components) && + (restriction_is_additive == f.restriction_is_additive)); }; 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 10aef4f63a..5e9e2b368d 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 @@ -10,39 +10,28 @@ -#if deal_II_dimension == 1 -template <> -FEDG_Q0<1>::FEDG_Q0 () : - FEQ1Mapping<1> (0, 1) +template +FEDG_Q0::FEDG_Q0 () : + FEQ1Mapping (0, + (dim==1 ? 1 : 0), + (dim==2 ? 1 : 0), + (dim==3 ? 1 : 0), + 1, + true) { - // for restriction and prolongation matrices: - // note that we do not add up all the - // contributions since then we would get - // two summands per vertex in 1d (four - // in 2d, etc), but only one per line dof. - // We could accomplish for that by dividing - // the vertex dof values by 2 (4, etc), but - // would get into trouble at the boundary - // of the domain since there only one - // cell contributes to a vertex. Rather, - // we do not add up the contributions but - // set them right into the matrices! - - // The restriction matrices got crazy values - // as it is yet not clear how they should work - // in the DG(0) case. In general - // the use of the restriction matrices - // is not yet finally decided about, too. - restriction[0](0,0) = 1e8; - restriction[1](0,0) = 1e8; - - prolongation[0](0,0) = 1.0; - prolongation[1](0,0) = 1.0; + for (unsigned int i=0; i::children_per_cell; ++i) + { + restriction[i](0,0) = 1./GeometryInfo::children_per_cell; + prolongation[i](0,0) = 1.0; + } }; +#if deal_II_dimension == 1 + + template <> void FEDG_Q0<1>::get_face_support_points (const DoFHandler<1>::face_iterator &, @@ -53,40 +42,6 @@ FEDG_Q0<1>::get_face_support_points (const DoFHandler<1>::face_iterator &, #endif - - -#if deal_II_dimension == 2 - -template <> -FEDG_Q0<2>::FEDG_Q0 () : - FEQ1Mapping<2> (0, 0, 1) -{ - // The restriction matrices got crazy values - // as it is yet not clear how they should work - // in the DG(0) case. In general - // the use of the restriction matrices - // is not yet finally decided about, too. - restriction[0](0,0) = 1e8; - restriction[1](0,0) = 1e8; - restriction[2](0,0) = 1e8; - restriction[3](0,0) = 1e8; - - prolongation[0](0,0) = 1.0; - - prolongation[1](0,0) = 1.0; - - prolongation[2](0,0) = 1.0; - - prolongation[3](0,0) = 1.0; -}; - - - -#endif - - - - template inline double diff --git a/deal.II/deal.II/source/fe/q1_mapping.cc b/deal.II/deal.II/source/fe/q1_mapping.cc index 8c1d44d1f5..2753e1e8a7 100644 --- a/deal.II/deal.II/source/fe/q1_mapping.cc +++ b/deal.II/deal.II/source/fe/q1_mapping.cc @@ -20,14 +20,16 @@ template <> FEQ1Mapping<1>::FEQ1Mapping (const unsigned int dofs_per_vertex, - const unsigned int dofs_per_line, - const unsigned int dofs_per_quad, - const unsigned int dofs_per_hex, - const unsigned int n_components) : + const unsigned int dofs_per_line, + const unsigned int dofs_per_quad, + const unsigned int dofs_per_hex, + const unsigned int n_components, + const bool restriction_is_additive) : FiniteElement<1> (FiniteElementData<1> (dofs_per_vertex, dofs_per_line, GeometryInfo<1>::vertices_per_cell, - n_components)) + n_components, + restriction_is_additive)) { Assert (dofs_per_quad==0, ExcInvalidData()); Assert (dofs_per_hex==0, ExcInvalidData()); @@ -141,15 +143,17 @@ void FEQ1Mapping<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell, template <> FEQ1Mapping<2>::FEQ1Mapping (const unsigned int dofs_per_vertex, - const unsigned int dofs_per_line, - const unsigned int dofs_per_quad, - const unsigned int dofs_per_hex, - const unsigned int n_components) : + const unsigned int dofs_per_line, + const unsigned int dofs_per_quad, + const unsigned int dofs_per_hex, + const unsigned int n_components, + const bool restriction_is_additive) : FiniteElement<2> (FiniteElementData<2> (dofs_per_vertex, dofs_per_line, dofs_per_quad, GeometryInfo<2>::vertices_per_cell, - n_components)) + n_components, + restriction_is_additive)) { Assert (dofs_per_hex == 0, ExcInvalidData()); }; @@ -313,16 +317,18 @@ void FEQ1Mapping<2>::get_normal_vectors (const DoFHandler<2>::cell_iterator &cel template <> FEQ1Mapping<3>::FEQ1Mapping (const unsigned int dofs_per_vertex, - const unsigned int dofs_per_line, - const unsigned int dofs_per_quad, - const unsigned int dofs_per_hex, - const unsigned int n_components) : + const unsigned int dofs_per_line, + const unsigned int dofs_per_quad, + const unsigned int dofs_per_hex, + const unsigned int n_components, + const bool restriction_is_additive) : FiniteElement<3> (FiniteElementData<3> (dofs_per_vertex, dofs_per_line, dofs_per_quad, dofs_per_hex, GeometryInfo<3>::vertices_per_cell, - n_components)) + n_components, + restriction_is_additive)) {};