]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
add restriction is additive flag and change the constructors
authorhartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 9 Aug 1999 15:45:18 +0000 (15:45 +0000)
committerhartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 9 Aug 1999 15:45:18 +0000 (15:45 +0000)
git-svn-id: https://svn.dealii.org/trunk@1647 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe.h
deal.II/deal.II/include/fe/q1_mapping.h
deal.II/deal.II/source/fe/fe.cc
deal.II/deal.II/source/fe/fe_lib.dg.constant.cc
deal.II/deal.II/source/fe/q1_mapping.cc

index 39618580d9665af694140748c163254c3c491cc6..c2390ec97bb063397fdb2b04939731cb269a481d 100644 (file)
@@ -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
index 60b9ff60c26a61b076cbe2cf22af03838e99bc93..e0638265a13dd0a4a5bcd1cc87886b258e97cc31 100644 (file)
@@ -19,9 +19,7 @@
  * are implemented here and do not have to be taken care of later.
  */
 template <int dim>
-class FEQ1Mapping
-  :
-  public FiniteElement<dim>
+class FEQ1Mapping : public FiniteElement<dim>
 {
   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
index ded20956be2333b030f79a231164858c89e472b6..a84a261f7ec26a8c8970fbc6afc55a41e5bc8091 100644 (file)
@@ -28,7 +28,8 @@ FiniteElementData<dim>::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<dim>::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<dim>::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<dim>::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 <int dim>
 FiniteElementData<dim>::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<dim>::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<dim>::operator== (const FiniteElementData<dim> &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));
 };
 
 
index 10aef4f63ab1915ac50e3b215db439c26b3e56df..5e9e2b368d8af6b734643ab72a7ec7259e5d3a5f 100644 (file)
 
 
 
-#if deal_II_dimension == 1
 
-template <>
-FEDG_Q0<1>::FEDG_Q0 () :
-               FEQ1Mapping<1> (0, 1)
+template <int dim>
+FEDG_Q0<dim>::FEDG_Q0 () :
+               FEQ1Mapping<dim> (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<GeometryInfo<dim>::children_per_cell; ++i)
+    { 
+      restriction[i](0,0) = 1./GeometryInfo<dim>::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 <int dim>
 inline
 double
index 8c1d44d1f5aa8b1cc7b90b096b41dd393e1de4c7..2753e1e8a753d9ee705ed456e70c1cc4a0a8256b 100644 (file)
 
 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))
 {};
 
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.