]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add support for face and subface values. Implement a few other missing things. Fix...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 18 Sep 2002 15:00:11 +0000 (15:00 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 18 Sep 2002 15:00:11 +0000 (15:00 +0000)
git-svn-id: https://svn.dealii.org/trunk@6452 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/fe/fe_nedelec.cc
deal.II/deal.II/source/fe/fe_nedelec_1d.cc
deal.II/deal.II/source/fe/fe_nedelec_3d.cc

index 2f283743cf17eee068cf1405091ce4ecc39984f0..2db8b0705f9037c9ae568d11b53b1d29318c04ca 100644 (file)
@@ -874,9 +874,9 @@ FE_Nedelec<dim>::fill_fe_values (const Mapping<dim>                   &mapping,
       std::vector<Tensor<2,dim> > shape_grads1 (n_q_points);
       std::vector<Tensor<2,dim> > shape_grads2 (n_q_points);
 
-      Assert (data.shape_values.n_rows() == dofs_per_cell * dim,
+      Assert (data.shape_gradients.n_rows() == dofs_per_cell * dim,
              ExcInternalError());
-      Assert (data.shape_values.n_cols() == n_q_points,
+      Assert (data.shape_gradients.n_cols() == n_q_points,
              ExcInternalError());
 
                                        // loop over all shape
@@ -934,107 +934,295 @@ FE_Nedelec<dim>::fill_fe_values (const Mapping<dim>                   &mapping,
     compute_2nd (mapping, cell, 0, mapping_data, fe_data, data);
   
   fe_data.first_cell = false;
-}
+};
 
 
 
 template <int dim>
 void
-FE_Nedelec<dim>::fill_fe_face_values (const Mapping<dim>                   &/*mapping*/,
-                                     const typename DoFHandler<dim>::cell_iterator &/*cell*/,
-                                     const unsigned int                    /*face*/,
-                                     const Quadrature<dim-1>              &/*quadrature*/,
-                                     typename Mapping<dim>::InternalDataBase       &/*mapping_data*/,
-                                     typename Mapping<dim>::InternalDataBase       &/*fedata*/,
-                                     FEValuesData<dim>                    &/*data*/) const
+FE_Nedelec<dim>::fill_fe_face_values (const Mapping<dim>                   &mapping,
+                                     const typename DoFHandler<dim>::cell_iterator &cell,
+                                     const unsigned int                    face,
+                                     const Quadrature<dim-1>              &quadrature,
+                                     typename Mapping<dim>::InternalDataBase       &mapping_data,
+                                     typename Mapping<dim>::InternalDataBase       &fedata,
+                                     FEValuesData<dim>                    &data) const
 {
-//TODO!!  
-//                                // convert data object to internal
-//                                // data for this class. fails with
-//                                // an exception if that is not
-//                                // possible
-//    InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
-
-//                                // offset determines which data set
-//                                // to take (all data sets for all
-//                                // faces are stored contiguously)
-//    const unsigned int offset = face * quadrature.n_quadrature_points;
-  
-//    const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
+                                  // convert data object to internal
+                                  // data for this class. fails with
+                                  // an exception if that is not
+                                  // possible
+  InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
+
+                                   // offset determines which data set
+                                  // to take (all data sets for all
+                                  // faces are stored contiguously)
+  const unsigned int offset = face * quadrature.n_quadrature_points;
+
+                                  // get the flags indicating the
+                                  // fields that have to be filled
+  const UpdateFlags flags(fe_data.current_update_flags());
+
+  const unsigned int n_q_points = quadrature.n_quadrature_points;
+                                 
+                                  // fill shape function
+                                  // values. these are vector-valued,
+                                  // so we have to transform
+                                  // them. since the output format
+                                  // (in data.shape_values) is a
+                                  // sequence of doubles (one for
+                                  // each non-zero shape function
+                                  // value, and for each quadrature
+                                  // point, rather than a sequence of
+                                  // small vectors, we have to use a
+                                  // number of conversions
+  if (flags & update_values)
+    {
+      Assert (fe_data.shape_values.n_cols() ==
+              GeometryInfo<dim>::faces_per_cell * n_q_points,
+              ExcInternalError());
+      
+      std::vector<Tensor<1,dim> > shape_values (n_q_points);
 
-//    for (unsigned int k=0; k<dofs_per_cell; ++k)
-//      {
-//        for (unsigned int i=0;i<quadrature.n_quadrature_points;++i)
-//     if (flags & update_values)
-//       data.shape_values(k,i) = fe_data.shape_values[k][i+offset];
+      Assert (data.shape_values.n_rows() == dofs_per_cell * dim,
+             ExcInternalError());
+      Assert (data.shape_values.n_cols() == n_q_points,
+             ExcInternalError());
       
-//        if (flags & update_gradients)
-//    {
-//       Assert (data.shape_gradients[k].size() + offset <=
-//               fe_data.shape_gradients[k].size(),
-//               ExcInternalError());    
-//     mapping.transform_covariant(data.shape_gradients[k],
-//                                 fe_data.shape_gradients[k],
-//                                 mapping_data, offset);
-// }  
-//      }
-
-//    if (flags & update_second_derivatives)
-//      compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
+      for (unsigned int k=0; k<dofs_per_cell; ++k)
+       {
+                                          // first transform shape
+                                          // values...
+         Assert (fe_data.shape_values[k].size() == n_q_points,
+                 ExcInternalError());
+         mapping.transform_covariant(&*shape_values.begin(),
+                                      &*shape_values.end(),
+                                      fe_data.shape_values[k].begin()+offset,
+                                      mapping_data);
+
+                                          // then copy over to target:
+         for (unsigned int q=0; q<n_q_points; ++q)
+           for (unsigned int d=0; d<dim; ++d)
+             data.shape_values[k*dim+d][q] = shape_values[q][d];
+       };
+    };
   
-//    fe_data.first_cell = false;
+      
+  if (flags & update_gradients)
+    {
+      Assert (fe_data.shape_gradients.n_cols() ==
+              GeometryInfo<dim>::faces_per_cell * n_q_points,
+              ExcInternalError());
+
+      std::vector<Tensor<2,dim> > shape_grads1 (n_q_points);
+      std::vector<Tensor<2,dim> > shape_grads2 (n_q_points);
+
+      Assert (data.shape_gradients.n_rows() == dofs_per_cell * dim,
+             ExcInternalError());
+      Assert (data.shape_gradients.n_cols() == n_q_points,
+             ExcInternalError());
+
+                                       // loop over all shape
+                                       // functions, and treat the
+                                       // gradients of each shape
+                                       // function at all quadrature
+                                       // points
+      for (unsigned int k=0; k<dofs_per_cell; ++k)
+       {
+                                           // treat the gradients of
+                                           // this particular shape
+                                           // function at all
+                                           // q-points. if Dv is the
+                                           // gradient of the shape
+                                           // function on the unit
+                                           // cell, then
+                                           // (J^-T)Dv(J^-1) is the
+                                           // value we want to have on
+                                           // the real cell. so, we
+                                           // will have to apply a
+                                           // covariant transformation
+                                           // to Dv twice. since the
+                                           // interface only allows
+                                           // multiplication with
+                                           // (J^-1) from the right,
+                                           // we have to trick a
+                                           // little in between
+         Assert (fe_data.shape_gradients[k].size() == n_q_points,
+                 ExcInternalError());
+                                           // do first transformation
+         mapping.transform_covariant(&*shape_grads1.begin(),
+                                      &*shape_grads1.end(),
+                                      fe_data.shape_gradients[k].begin()+offset,
+                                      mapping_data);
+                                           // transpose matrix
+          for (unsigned int q=0; q<n_q_points; ++q)
+            shape_grads2[q] = transpose(shape_grads1[q]);
+                                           // do second transformation
+         mapping.transform_covariant(&*shape_grads1.begin(),
+                                      &*shape_grads1.end(),
+                                      &*shape_grads2.begin(),
+                                      mapping_data);
+                                           // transpose back
+          for (unsigned int q=0; q<n_q_points; ++q)
+            shape_grads2[q] = transpose(shape_grads1[q]);
+          
+                                          // then copy over to target:
+         for (unsigned int q=0; q<n_q_points; ++q)
+           for (unsigned int d=0; d<dim; ++d)
+             data.shape_gradients[k*dim+d][q] = shape_grads2[q][d];
+       };
+    }
+
+  if (flags & update_second_derivatives)
+    compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
+  
+  fe_data.first_cell = false;
 }
 
 
 
 template <int dim>
 void
-FE_Nedelec<dim>::fill_fe_subface_values (const Mapping<dim>                   &/*mapping*/,
-                                        const typename DoFHandler<dim>::cell_iterator &/*cell*/,
-                                        const unsigned int                    /*face*/,
-                                        const unsigned int                    /*subface*/,
-                                        const Quadrature<dim-1>              &/*quadrature*/,
-                                        typename Mapping<dim>::InternalDataBase       &/*mapping_data*/,
-                                        typename Mapping<dim>::InternalDataBase       &/*fedata*/,
-                                        FEValuesData<dim>                    &/*data*/) const
+FE_Nedelec<dim>::fill_fe_subface_values (const Mapping<dim>                   &mapping,
+                                        const typename DoFHandler<dim>::cell_iterator &cell,
+                                        const unsigned int                    face,
+                                        const unsigned int                    subface,
+                                        const Quadrature<dim-1>              &quadrature,
+                                        typename Mapping<dim>::InternalDataBase       &mapping_data,
+                                        typename Mapping<dim>::InternalDataBase       &fedata,
+                                        FEValuesData<dim>                    &data) const
 {
-//TODO!!  
-//                                // convert data object to internal
-//                                // data for this class. fails with
-//                                // an exception if that is not
-//                                // possible
-//    InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
-
-//                                // offset determines which data set
-//                                // to take (all data sets for all
-//                                // sub-faces are stored contiguously)
-//    const unsigned int offset = (face * GeometryInfo<dim>::subfaces_per_face + subface)
-//                           * quadrature.n_quadrature_points;
-
-//    const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
-
-//    for (unsigned int k=0; k<dofs_per_cell; ++k)
-//      {
-//        for (unsigned int i=0;i<quadrature.n_quadrature_points;++i)
-//     if (flags & update_values)
-//       data.shape_values(k,i) = fe_data.shape_values[k][i+offset];
+                                  // convert data object to internal
+                                  // data for this class. fails with
+                                  // an exception if that is not
+                                  // possible
+  InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
+
+                                   // offset determines which data set
+                                  // to take (all data sets for all
+                                  // faces are stored contiguously)
+  const unsigned int offset = ((face * GeometryInfo<dim>::subfaces_per_face + subface)
+                               * quadrature.n_quadrature_points);
+
+                                  // get the flags indicating the
+                                  // fields that have to be filled
+  const UpdateFlags flags(fe_data.current_update_flags());
+
+  const unsigned int n_q_points = quadrature.n_quadrature_points;
+                                 
+                                  // fill shape function
+                                  // values. these are vector-valued,
+                                  // so we have to transform
+                                  // them. since the output format
+                                  // (in data.shape_values) is a
+                                  // sequence of doubles (one for
+                                  // each non-zero shape function
+                                  // value, and for each quadrature
+                                  // point, rather than a sequence of
+                                  // small vectors, we have to use a
+                                  // number of conversions
+  if (flags & update_values)
+    {
+      Assert (fe_data.shape_values.n_cols() ==
+              GeometryInfo<dim>::faces_per_cell * n_q_points,
+              ExcInternalError());
       
-//        if (flags & update_gradients)
-//  {
-//    Assert (data.shape_gradients[k].size() + offset<=
-//               fe_data.shape_gradients[k].size(),
-//               ExcInternalError());
-//     mapping.transform_covariant(data.shape_gradients[k],
-//                                 fe_data.shape_gradients[k],
-//                                 mapping_data, offset);
-// };  
-//      }
+      std::vector<Tensor<1,dim> > shape_values (n_q_points);
+
+      Assert (data.shape_values.n_rows() == dofs_per_cell * dim,
+             ExcInternalError());
+      Assert (data.shape_values.n_cols() == n_q_points,
+             ExcInternalError());
+      
+      for (unsigned int k=0; k<dofs_per_cell; ++k)
+       {
+                                          // first transform shape
+                                          // values...
+         Assert (fe_data.shape_values[k].size() == n_q_points,
+                 ExcInternalError());
+         mapping.transform_covariant(&*shape_values.begin(),
+                                      &*shape_values.end(),
+                                      fe_data.shape_values[k].begin()+offset,
+                                      mapping_data);
+
+                                          // then copy over to target:
+         for (unsigned int q=0; q<n_q_points; ++q)
+           for (unsigned int d=0; d<dim; ++d)
+             data.shape_values[k*dim+d][q] = shape_values[q][d];
+       };
+    };
   
-//    if (flags & update_second_derivatives)
-//      compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
+      
+  if (flags & update_gradients)
+    {
+      Assert (fe_data.shape_gradients.n_cols() ==
+              GeometryInfo<dim>::faces_per_cell * n_q_points,
+              ExcInternalError());
+
+      std::vector<Tensor<2,dim> > shape_grads1 (n_q_points);
+      std::vector<Tensor<2,dim> > shape_grads2 (n_q_points);
+
+      Assert (data.shape_gradients.n_rows() == dofs_per_cell * dim,
+             ExcInternalError());
+      Assert (data.shape_gradients.n_cols() == n_q_points,
+             ExcInternalError());
+
+                                       // loop over all shape
+                                       // functions, and treat the
+                                       // gradients of each shape
+                                       // function at all quadrature
+                                       // points
+      for (unsigned int k=0; k<dofs_per_cell; ++k)
+       {
+                                           // treat the gradients of
+                                           // this particular shape
+                                           // function at all
+                                           // q-points. if Dv is the
+                                           // gradient of the shape
+                                           // function on the unit
+                                           // cell, then
+                                           // (J^-T)Dv(J^-1) is the
+                                           // value we want to have on
+                                           // the real cell. so, we
+                                           // will have to apply a
+                                           // covariant transformation
+                                           // to Dv twice. since the
+                                           // interface only allows
+                                           // multiplication with
+                                           // (J^-1) from the right,
+                                           // we have to trick a
+                                           // little in between
+         Assert (fe_data.shape_gradients[k].size() == n_q_points,
+                 ExcInternalError());
+                                           // do first transformation
+         mapping.transform_covariant(&*shape_grads1.begin(),
+                                      &*shape_grads1.end(),
+                                      fe_data.shape_gradients[k].begin()+offset,
+                                      mapping_data);
+                                           // transpose matrix
+          for (unsigned int q=0; q<n_q_points; ++q)
+            shape_grads2[q] = transpose(shape_grads1[q]);
+                                           // do second transformation
+         mapping.transform_covariant(&*shape_grads1.begin(),
+                                      &*shape_grads1.end(),
+                                      &*shape_grads2.begin(),
+                                      mapping_data);
+                                           // transpose back
+          for (unsigned int q=0; q<n_q_points; ++q)
+            shape_grads2[q] = transpose(shape_grads1[q]);
+          
+                                          // then copy over to target:
+         for (unsigned int q=0; q<n_q_points; ++q)
+           for (unsigned int d=0; d<dim; ++d)
+             data.shape_gradients[k*dim+d][q] = shape_grads2[q][d];
+       };
+    }
+
+  if (flags & update_second_derivatives)
+    compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
   
-//    fe_data.first_cell = false;
-}
+  fe_data.first_cell = false;
+};
 
 
 
@@ -1066,16 +1254,31 @@ FE_Nedelec<dim>::has_support_on_face (const unsigned int shape_index,
          ExcIndexRange (shape_index, 0, this->dofs_per_cell));
   Assert (face_index < GeometryInfo<dim>::faces_per_cell,
          ExcIndexRange (face_index, 0, GeometryInfo<dim>::faces_per_cell));
-         
-//TODO: fix for higher orders. correct now for lowest order, all dimensions  
-//TODO!!
-// can this be done in a way that is dimension and degree independent?
+
+  switch (degree)
+    {
+      case 1:
+      {
+                                         // only on non-adjacent faces
+                                         // are the values actually
+                                         // zero. list these in a
+                                         // table
+        const unsigned int opposite_faces_2d[GeometryInfo<2>::faces_per_cell]
+          = { 2, 3, 0, 1};
+        const unsigned int opposite_faces_3d[GeometryInfo<3>::faces_per_cell]
+          = { 1, 0, 4, 5, 2, 3};
+        switch (dim)
+          {
+            case 2:  return (face_index != opposite_faces_2d[shape_index]);
+            case 3:  return (face_index != opposite_faces_3d[shape_index]);
+            default: Assert (false, ExcNotImplemented());
+          };
+      };
+      
+      default:  // other degree
+            Assert (false, ExcNotImplemented());
+    };
   
-                                  // all degrees of freedom are on
-                                  // lines, so also on a face. the
-                                  // question is whether it has
-                                  // support on this particular face
-  Assert (false, ExcNotImplemented());
   return true;
 }
 
index 10cf287383a3a260e88c7ce4eb59b916c2094602..1babaddbd8c097161d25cc88e1c9803173cf701e 100644 (file)
@@ -13,8 +13,8 @@
 
 
 
-// only compile this file if in 1d. note that Nedelec elemets do not
-// make much sense in 1d, so this file only contains dummy
+// only compile this file if in 1d. note that Nedelec elements do not
+// make much sense in 1d anyway, so this file only contains dummy
 // implementations to avoid linker errors due to missing symbols
 #if deal_II_dimension == 1
 
index 38638df49fa7e8aea9ab4c5571a264f4de2858e1..6afa6fbe65168f5eadd56678d23a718f5cea0aed 100644 (file)
@@ -24,8 +24,6 @@
 // the child cell are obtained from the degrees of freedom on the
 // mother cell
 //
-// TODO: [Anna] check whether the following paragraph is correct. if so, then please multiply the values in the eight following matrices by two
-
 // note the following: since the shape functions themselves and not
 // only the gradients are transformed using the mapping object from
 // the unit cell to the real cell, the actual values of the function

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.