]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Added functionality which corrects some errors in the
authorkayser-herold <kayser-herold@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 29 Jul 2006 20:40:45 +0000 (20:40 +0000)
committerkayser-herold <kayser-herold@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 29 Jul 2006 20:40:45 +0000 (20:40 +0000)
code for the Raviart-Thomas and ABF elements. The first
correction is a sign change of the Face-DoFs, which is
required to work on more complicated domains. This fix
is currently only applied in the 2D case. The second
correction is related to the coordinate transformation
from the reference element to the global coordinate
system. This transformation now becomes the Piola
transformation, which according to Brezzi/Fortin 1992
is the correct one for subspaces of H_div.

git-svn-id: https://svn.dealii.org/trunk@13504 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/fe/fe_poly_tensor.cc
deal.II/deal.II/source/fe/fe_raviart_thomas.cc
deal.II/deal.II/source/fe/fe_system.cc

index efe91afaaf49a3d2aa4980d7555eb85dd52102a7..6df60434200fda41133d3ee795885697c4ccd516 100644 (file)
@@ -282,6 +282,39 @@ FE_PolyTensor<POLY,dim>::fill_fe_values (
         ExcDimensionMismatch(fe_data.shape_values.size(), n_dofs));
   Assert(!(flags & update_values) || fe_data.shape_values[0].size() == n_quad,
         ExcDimensionMismatch(fe_data.shape_values[0].size(), n_quad));
+
+  // Create table with sign changes, due to the special structure of the RT elements.
+  // TODO: Preliminary hack to demonstrate the overall prinicple!
+
+  // Compute eventual sign changes depending on the neighborhood
+  // between two faces.
+  std::vector<double> sign_change (n_dofs, 1.0);
+  
+#if deal_II_dimension > 1
+  const unsigned int dofs_per_face = this->dofs_per_face;
+
+  if (dim == 2)
+    {
+      for (unsigned int f = GeometryInfo<dim>::faces_per_cell / 2; 
+          f < GeometryInfo<dim>::faces_per_cell; ++f)
+       {
+         typename Triangulation<dim>::face_iterator face = cell->face (f);
+         if (!face->at_boundary ())
+           {
+             unsigned int nn = cell->neighbor_of_neighbor (f);
+             if (nn < GeometryInfo<dim>::faces_per_cell / 2)
+               {
+                 for (unsigned int j = 0; j < dofs_per_face; ++j)
+                   sign_change[f * dofs_per_face + j] = -1.0;
+               }
+           }
+       }
+    }
+  else
+    {
+      // TODO: Think about 3D!.
+    }
+#endif
   
   for (unsigned int i=0; i<n_dofs; ++i)
     {
@@ -312,8 +345,15 @@ FE_PolyTensor<POLY,dim>::fill_fe_values (
                    
                                                     // then copy over to target:
                    for (unsigned int k=0; k<n_quad; ++k)
-                     for (unsigned int d=0; d<dim; ++d)
-                       data.shape_values(first+d,k) = shape_values[k][d];
+                     {
+                       // Recompute determinant
+                       double J = 1.0;
+                       if (mapping_type == contravariant)
+                         J = data.JxW_values[k] / quadrature.weight(k);
+
+                       for (unsigned int d=0; d<dim; ++d)
+                         data.shape_values(first+d,k) = sign_change[i] * (shape_values[k][d] / J);
+                     }
                  }
                break;
              default:
@@ -351,16 +391,25 @@ FE_PolyTensor<POLY,dim>::fill_fe_values (
                
                for (unsigned int k=0; k<n_quad; ++k)
                  for (unsigned int d=0;d<dim;++d)
-                   data.shape_gradients[first+d][k] = shape_grads1[k][d];
+                   data.shape_gradients[first+d][k] = shape_grads2[k][d];
                break;
              case contravariant:
-               Assert(false, ExcNotImplemented());             
                mapping.transform_covariant(fe_data.shape_grads[i], 0,
                                            shape_grads1,
                                            mapping_data);
+
+               mapping.transform_contravariant(shape_grads1, 0,
+                                               shape_grads2,
+                                               mapping_data);
+
                for (unsigned int k=0; k<n_quad; ++k)
                  for (unsigned int d=0;d<dim;++d)
-                   data.shape_gradients[first+d][k] = shape_grads1[k][d];
+                   {
+                     // Recompute determinant
+                     double J = data.JxW_values[k] / quadrature.weight(k);
+                     data.shape_gradients[first+d][k] = sign_change[i] * 
+                       shape_grads2[k][d] / J;
+                   }
                break;
              default:
                Assert(false, ExcNotImplemented());
@@ -411,6 +460,38 @@ FE_PolyTensor<POLY,dim>::fill_fe_face_values (
 //           && dynamic_cast<const MappingCartesian<dim>*>(&mapping) != 0),
 //      ExcNotImplemented());
 //TODO: Size assertions
+
+  // Create table with sign changes, due to the special structure of the RT elements.
+  // TODO: Preliminary hack to demonstrate the overall prinicple!
+
+  // Compute eventual sign changes depending on the neighborhood
+  // between two faces.
+  std::vector<double> sign_change (n_dofs, 1.0);
+#if deal_II_dimension > 1
+  const unsigned int dofs_per_face = this->dofs_per_face;
+  
+  if (dim == 2)
+    {
+      for (unsigned int f = GeometryInfo<dim>::faces_per_cell / 2; 
+          f < GeometryInfo<dim>::faces_per_cell; ++f)
+       {
+         typename Triangulation<dim>::face_iterator face = cell->face (f);
+         if (!face->at_boundary ())
+           {
+             unsigned int nn = cell->neighbor_of_neighbor (f);
+             if (nn < GeometryInfo<dim>::faces_per_cell / 2)
+               {
+                 for (unsigned int j = 0; j < dofs_per_face; ++j)
+                   sign_change[f * dofs_per_face + j] = -1.0;
+               }
+           }
+       }
+    }
+  else
+    {
+      // TODO: Think about 3D!.
+    }
+#endif
   
   for (unsigned int i=0; i<n_dofs; ++i)
     {
@@ -441,8 +522,15 @@ FE_PolyTensor<POLY,dim>::fill_fe_face_values (
                    
                                                     // then copy over to target:
                    for (unsigned int k=0; k<n_quad; ++k)
-                     for (unsigned int d=0; d<dim; ++d)
-                       data.shape_values(first+d,k) = shape_values[k][d];
+                     {
+                       // Recompute determinant
+                       double J = 1.0;
+                       if (mapping_type == contravariant)
+                         J = data.JxW_values[k] / quadrature.weight(k);
+
+                       for (unsigned int d=0; d<dim; ++d)
+                         data.shape_values(first+d,k) = sign_change[i] * (shape_values[k][d] / J);
+                     }
                  }
                break;
              default:
@@ -481,6 +569,26 @@ FE_PolyTensor<POLY,dim>::fill_fe_face_values (
                  for (unsigned int d=0;d<dim;++d)
                    data.shape_gradients[first+d][k] = shape_grads1[k][d];
                break;
+
+             case contravariant:
+               mapping.transform_covariant(fe_data.shape_grads[i], offset,
+                                           shape_grads1,
+                                           mapping_data);
+
+               mapping.transform_contravariant(shape_grads1, 0,
+                                               shape_grads2,
+                                               mapping_data);
+
+               for (unsigned int k=0; k<n_quad; ++k)
+                 for (unsigned int d=0;d<dim;++d)
+                   {
+                     // Recompute determinant
+                     double J = data.JxW_values[k] / quadrature.weight(k);
+                     data.shape_gradients[first+d][k] = sign_change[i] * 
+                       shape_grads2[k][d] / J;
+                   }
+               break;
+               
              default:
                Assert(false, ExcNotImplemented());
            }
index bc6aaa047dbe6f529e1aafdd796de4fbee67fe37..0190e4849e401eaebe2dab7fecf53967220bc580 100644 (file)
@@ -43,7 +43,7 @@ FE_RaviartThomas<dim>::FE_RaviartThomas (const unsigned int deg)
   Assert (dim >= 2, ExcImpossibleInDim(dim));
   const unsigned int n_dofs = this->dofs_per_cell;
   
-  this->mapping_type = this->covariant;
+  this->mapping_type = this->contravariant;
                                   // First, initialize the
                                   // generalized support points and
                                   // quadrature weights, since they
@@ -449,11 +449,17 @@ FE_RaviartThomas<dim>::update_each (const UpdateFlags flags) const
   UpdateFlags out = update_default;
 
   if (flags & update_values)
-    out |= update_values             | update_covariant_transformation;
+    out |= update_values             | update_covariant_transformation
+                                     | update_contravariant_transformation 
+                                     | update_JxW_values;
   if (flags & update_gradients)
-    out |= update_gradients          | update_covariant_transformation;
+    out |= update_gradients          | update_covariant_transformation 
+                                     | update_contravariant_transformation
+                                     | update_JxW_values;
+  //TODO: Set update flags appropriately and figure out, how the second
+  // derivatives for the RT elements can be computed correctly.
   if (flags & update_second_derivatives)
-    out |= update_second_derivatives | update_covariant_transformation;
+    out |= update_second_derivatives | update_contravariant_transformation;
 
   return out;
 }
index 072150622671e2356f72cc588e64b293ae691bd1..9a85175dd1945280fcf6545668ccac23569cbda6 100644 (file)
@@ -821,6 +821,19 @@ compute_fill (const Mapping<dim>                   &mapping,
          FEValuesData<dim> &
             base_data    = fe_data.get_fe_values_data(base_no);
 
+                                          //TODO: Think about a smarter alternative
+                                          // Copy quadrature points. These
+                                          // are required for computing the
+                                          // determinant in the FEPolyTensor
+                                          // class. The determinant is
+                                          // one ingredient of the Piola
+                                          // transformation, which is applied
+                                          // to correctly map the RT space from
+                                          // the reference element to the global
+                                          // coordinate system.
+         base_data.JxW_values = data.JxW_values;
+
+
                                            // Make sure that in the
                                            // case of fill_fe_values
                                            // the data is only copied

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.