]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add get_linear_shape().
authorJaeryun Yim <jaeryun.yim@gmail.com>
Thu, 14 Jul 2016 12:04:15 +0000 (21:04 +0900)
committerJaeryun Yim <jaeryun.yim@gmail.com>
Fri, 16 Sep 2016 12:07:22 +0000 (21:07 +0900)
+ Edit private/protected.

include/deal.II/fe/fe_p1nc.h
source/fe/fe_p1nc.cc

index ff32ebcff6edfe44d2d013e55010dc1cc9b3eb67..3a95887c77f0f54ba552060aea6ff69cf2d5c260 100644 (file)
@@ -31,8 +31,36 @@ DEAL_II_NAMESPACE_OPEN
 
 class FE_P1NC : public FiniteElement<2,2>
 {
-private:
 
+public:
+  /**
+   * Constructor for nonparametric version of P1 nonconforming element.
+   */
+  FE_P1NC() ;
+
+  /**
+   * Return the name of the class for the element.
+   */
+  virtual std::string get_name () const ;
+
+  /**
+   * Return the update flags which are needed.
+   */
+  virtual UpdateFlags     requires_update_flags (const UpdateFlags flags) const ;
+
+  /**
+   * Copy constructor.
+   */
+  virtual FiniteElement<2,2> *clone () const ;
+
+  /**
+   * Destructor.
+   */
+  virtual ~FE_P1NC ();
+
+
+
+private:
 
   static
   std::vector<ComponentMask>
@@ -47,6 +75,19 @@ private:
   get_dpo_vector ();
 
 
+  /**
+   * Compute the linear shape functions phi(x,y) = ax + by + c
+   * such that each midpoint value on two connecting edges is a half,
+   * and two other midpoint values are all zero.
+   */
+  static
+  void
+  get_linear_shape (const Triangulation<2,2>::cell_iterator &cell,
+                    std::vector<double> &a,
+                    std::vector<double> &b,
+                    std::vector<double> &c);
+
+
 
   /**
    * Do the work which is needed before cellwise data computation.
@@ -65,13 +106,13 @@ private:
   /**
    * Compute the data on the current cell.
    */
-   virtual
+  virtual
   void
   fill_fe_values (const Triangulation<2,2>::cell_iterator           &cell,
                   const CellSimilarity::Similarity                   ,
                   const Quadrature<2>                               &quadrature,
                   const Mapping<2,2>                                &mapping,
-                  const Mapping<2,2>::InternalDataBase              &,
+                  const Mapping<2,2>::InternalDataBase &,
                   const internal::FEValues::MappingRelatedData<2,2> &mapping_data,
                   const FiniteElement<2,2>::InternalDataBase        &fe_internal,
                   internal::FEValues::FiniteElementRelatedData<2,2> &output_data) const;
@@ -112,35 +153,6 @@ private:
 
 
 
-public:
-  /**
-   * Constructor for nonparametric version of P1 nonconforming element.
-   */
-  FE_P1NC() ;
-
-  /**
-   * Return the name of the class for the element.
-   */
-  virtual std::string get_name () const ;
-
-  /**
-   * Return the update flags which are needed.
-   */
-  virtual UpdateFlags     requires_update_flags (const UpdateFlags flags) const ;
-
-  /**
-   * Copy constructor.
-   */
-  virtual FiniteElement<2,2> *clone () const ;
-
-  /**
-   * Destructor.
-   */
-  virtual ~FE_P1NC ();
-
-
-protected:
-
   /**
    * Create the constraints matrix for hanging edges.
    */
index 161b461344a2c1073991e99bd5f60a18f54bbc5c..e618997cdfc160e73d8d8dc15eb66c67dbeec0ff 100644 (file)
 DEAL_II_NAMESPACE_OPEN
 
 
+FE_P1NC::FE_P1NC()
+  :
+  FiniteElement<2,2>(FiniteElementData<2>(get_dpo_vector(),
+                                          1,
+                                          1,
+                                          FiniteElementData<2>::L2),
+                     std::vector<bool> (1, false),
+                     get_nonzero_component())
+{
+
+  // face support points: 2 end vertices
+  unit_face_support_points.resize(2);
+  unit_face_support_points[0][0] = 0.0 ;
+  unit_face_support_points[1][0] = 1.0 ;
+
+
+  // initialize constraints matrix
+  initialize_constraints () ;
+}
+
+
+
+std::string FE_P1NC::get_name () const
+{
+  return "FE_P1NC";
+}
+
+
+UpdateFlags     FE_P1NC::requires_update_flags (const UpdateFlags flags) const
+{
+  UpdateFlags out = update_default;
+
+  if (flags & update_values)
+    out |= update_values | update_quadrature_points ;
+  if (flags & update_gradients)
+    out |= update_gradients ;
+  if (flags & update_cell_normal_vectors)
+    out |= update_cell_normal_vectors | update_JxW_values;
+  if (flags & update_hessians)
+    out |= update_hessians;
+  if (flags & update_hessians)
+    out |= update_hessians;
+
+  return out;
+}
+
+
+FiniteElement<2,2> *FE_P1NC::clone () const
+{
+  return new FE_P1NC(*this);
+}
+
+FE_P1NC::~FE_P1NC () {}
+
+
 std::vector<ComponentMask>
 FE_P1NC::get_nonzero_component()
 {
@@ -41,44 +96,12 @@ FE_P1NC::get_dpo_vector ()
 }
 
 
-FiniteElement<2,2>::InternalDataBase *
-FE_P1NC::get_data (const UpdateFlags update_flags,
-                                const Mapping<2,2> &,
-                                const Quadrature<2> &,
-                               dealii::internal::FEValues::FiniteElementRelatedData<2,2> &) const
-{
-  FiniteElement<2,2>::InternalDataBase *data = new FiniteElement<2,2>::InternalDataBase;
-
-  data->update_each = requires_update_flags(update_flags);
-
-  return data;
-
-
-
-
-}
-
-
-
 void
-FE_P1NC::fill_fe_values (const Triangulation<2,2>::cell_iterator           &cell,
-                                      const CellSimilarity::Similarity                   ,
-                                      const Quadrature<2>                               &quadrature,
-                                      const Mapping<2,2>                                &mapping,
-                                      const Mapping<2,2>::InternalDataBase              &,
-                                      const internal::FEValues::MappingRelatedData<2,2> &mapping_data,
-                                      const FiniteElement<2,2>::InternalDataBase        &fe_internal,
-                                      internal::FEValues::FiniteElementRelatedData<2,2> &output_data) const
+FE_P1NC::get_linear_shape (const Triangulation<2,2>::cell_iterator &cell,
+                           std::vector<double> &a,
+                           std::vector<double> &b,
+                           std::vector<double> &c)
 {
-  const UpdateFlags flags(fe_internal.update_each) ;
-
-  const unsigned int n_q_points = mapping_data.quadrature_points.size();
-
-  std::vector<double> values(flags & update_values ? this->dofs_per_cell : 0);
-  std::vector<Tensor<1,2> > grads(flags & update_gradients ? this->dofs_per_cell : 0);
-
-
-
   // edge midpoints
   std::vector<Point<2> > mpt(4) ;
 
@@ -100,10 +123,8 @@ FE_P1NC::fill_fe_values (const Triangulation<2,2>::cell_iterator           &cell
   cpt(1) = (mpt[0](1) + mpt[1](1) + mpt[2](1) + mpt[3](1))/4.0 ;
 
 
-  // basis functions with a half value: phi(x,y) = ax + by + c
-  std::vector<double> a(4), b(4), c(4) ;
-  double det ;
 
+  double det ;
   det = (mpt[0](0)-mpt[1](0))*(mpt[2](1)-mpt[3](1)) - (mpt[2](0)-mpt[3](0))*(mpt[0](1)-mpt[1](1)) ;
 
   a[0] = ((mpt[2](1)-mpt[3](1))*(0.5) -(mpt[0](1)-mpt[1](1))*(0.5))/det ;
@@ -121,6 +142,51 @@ FE_P1NC::fill_fe_values (const Triangulation<2,2>::cell_iterator           &cell
   c[2] = 0.25 - cpt(0)*a[2] - cpt(1)*b[2] ;
   c[3] = 0.25 - cpt(0)*a[3] - cpt(1)*b[3] ;
 
+}
+
+
+
+
+FiniteElement<2,2>::InternalDataBase *
+FE_P1NC::get_data (const UpdateFlags update_flags,
+                   const Mapping<2,2> &,
+                   const Quadrature<2> &,
+                   dealii::internal::FEValues::FiniteElementRelatedData<2,2> &) const
+{
+  FiniteElement<2,2>::InternalDataBase *data = new FiniteElement<2,2>::InternalDataBase;
+
+  data->update_each = requires_update_flags(update_flags);
+
+  return data;
+
+
+
+
+}
+
+
+
+void
+FE_P1NC::fill_fe_values (const Triangulation<2,2>::cell_iterator           &cell,
+                         const CellSimilarity::Similarity                   ,
+                         const Quadrature<2>                               &quadrature,
+                         const Mapping<2,2>                                &mapping,
+                         const Mapping<2,2>::InternalDataBase &,
+                         const internal::FEValues::MappingRelatedData<2,2> &mapping_data,
+                         const FiniteElement<2,2>::InternalDataBase        &fe_internal,
+                         internal::FEValues::FiniteElementRelatedData<2,2> &output_data) const
+{
+  const UpdateFlags flags(fe_internal.update_each) ;
+
+  const unsigned int n_q_points = mapping_data.quadrature_points.size();
+
+  std::vector<double> values(flags & update_values ? this->dofs_per_cell : 0);
+  std::vector<Tensor<1,2> > grads(flags & update_gradients ? this->dofs_per_cell : 0);
+
+
+  // linear shape with a half value: phi(x,y) = ax + by + c
+  std::vector<double> a(4), b(4), c(4);
+  get_linear_shape (cell, a, b, c);
 
 
   // compute basis functions
@@ -145,22 +211,22 @@ FE_P1NC::fill_fe_values (const Triangulation<2,2>::cell_iterator           &cell
 
       }
 
-      // hessian
-      std::vector<Tensor<2,2> > hessians(flags & update_hessians ? this->dofs_per_cell : 0);
-      if (flags & update_hessians)
-      {
-        for (unsigned int i=0; i<n_q_points; ++i)
+  // hessian
+  std::vector<Tensor<2,2> > hessians(flags & update_hessians ? this->dofs_per_cell : 0);
+  if (flags & update_hessians)
+    {
+      for (unsigned int i=0; i<n_q_points; ++i)
         {
-            for (unsigned int k=0; k<this->dofs_per_cell; ++k)
+          for (unsigned int k=0; k<this->dofs_per_cell; ++k)
             {
-                hessians[k][0][0] = 0.0 ;
-                hessians[k][0][1] = 0.0 ;
-                hessians[k][1][0] = 0.0 ;
-                hessians[k][1][1] = 0.0 ;
-                output_data.shape_hessians[k][i] = hessians[k];
+              hessians[k][0][0] = 0.0 ;
+              hessians[k][0][1] = 0.0 ;
+              hessians[k][1][0] = 0.0 ;
+              hessians[k][1][1] = 0.0 ;
+              output_data.shape_hessians[k][i] = hessians[k];
             }
         }
-      }
+    }
 
   // When this function is called for the graphic out,
   // MappingRelatedData does not work properly.
@@ -185,13 +251,13 @@ FE_P1NC::fill_fe_values (const Triangulation<2,2>::cell_iterator           &cell
 
 void
 FE_P1NC::fill_fe_face_values (const Triangulation<2,2>::cell_iterator           &cell,
-                                           const unsigned int                                 face_no,
-                                           const Quadrature<1>                                             &quadrature,
-                                           const Mapping<2,2>                                         &mapping,
-                                           const Mapping<2,2>::InternalDataBase              &mapping_internal,
-                                           const dealii::internal::FEValues::MappingRelatedData<2,2> &mapping_data,
-                                           const InternalDataBase                                              &fe_internal,
-                                           dealii::internal::FEValues::FiniteElementRelatedData<2,2> &output_data) const
+                              const unsigned int                                 face_no,
+                              const Quadrature<1>                                             &quadrature,
+                              const Mapping<2,2>                                         &mapping,
+                              const Mapping<2,2>::InternalDataBase              &mapping_internal,
+                              const dealii::internal::FEValues::MappingRelatedData<2,2> &mapping_data,
+                              const InternalDataBase                                              &fe_internal,
+                              dealii::internal::FEValues::FiniteElementRelatedData<2,2> &output_data) const
 
 {
   const UpdateFlags flags(fe_internal.update_each) ;
@@ -202,49 +268,9 @@ FE_P1NC::fill_fe_face_values (const Triangulation<2,2>::cell_iterator
   std::vector<Tensor<1,2> > grads(flags & update_gradients ? this->dofs_per_cell : 0);
 
 
-
-  // edge midpoints
-  std::vector<Point<2> > mpt(4) ;
-
-  mpt[0](0) = (cell->vertex(0)(0) + cell->vertex(2)(0))/2.0 ;
-  mpt[0](1) = (cell->vertex(0)(1) + cell->vertex(2)(1))/2.0 ;
-
-  mpt[1](0) = (cell->vertex(1)(0) + cell->vertex(3)(0))/2.0 ;
-  mpt[1](1) = (cell->vertex(1)(1) + cell->vertex(3)(1))/2.0 ;
-
-  mpt[2](0) = (cell->vertex(0)(0) + cell->vertex(1)(0))/2.0 ;
-  mpt[2](1) = (cell->vertex(0)(1) + cell->vertex(1)(1))/2.0 ;
-
-  mpt[3](0) = (cell->vertex(2)(0) + cell->vertex(3)(0))/2.0 ;
-  mpt[3](1) = (cell->vertex(2)(1) + cell->vertex(3)(1))/2.0 ;
-
-  // center point
-  Point<2> cpt ;
-  cpt(0) = (mpt[0](0) + mpt[1](0) + mpt[2](0) + mpt[3](0))/4.0 ;
-  cpt(1) = (mpt[0](1) + mpt[1](1) + mpt[2](1) + mpt[3](1))/4.0 ;
-
-
-  // basis functions with a half value: phi(x,y) = ax + by + c
-  std::vector<double> a(4), b(4), c(4) ;
-  double det ;
-
-  det = (mpt[0](0)-mpt[1](0))*(mpt[2](1)-mpt[3](1)) - (mpt[2](0)-mpt[3](0))*(mpt[0](1)-mpt[1](1)) ;
-
-  a[0] = ((mpt[2](1)-mpt[3](1))*(0.5) -(mpt[0](1)-mpt[1](1))*(0.5))/det ;
-  a[1] = ((mpt[2](1)-mpt[3](1))*(-0.5) -(mpt[0](1)-mpt[1](1))*(0.5))/det ;
-  a[2] = ((mpt[2](1)-mpt[3](1))*(0.5) -(mpt[0](1)-mpt[1](1))*(-0.5))/det ;
-  a[3] = ((mpt[2](1)-mpt[3](1))*(-0.5) -(mpt[0](1)-mpt[1](1))*(-0.5))/det ;
-
-  b[0] = (-(mpt[2](0)-mpt[3](0))*(0.5) +(mpt[0](0)-mpt[1](0))*(0.5))/det ;
-  b[1] = (-(mpt[2](0)-mpt[3](0))*(-0.5) +(mpt[0](0)-mpt[1](0))*(0.5))/det ;
-  b[2] = (-(mpt[2](0)-mpt[3](0))*(0.5) +(mpt[0](0)-mpt[1](0))*(-0.5))/det ;
-  b[3] = (-(mpt[2](0)-mpt[3](0))*(-0.5) +(mpt[0](0)-mpt[1](0))*(-0.5))/det ;
-
-  c[0] = 0.25 - cpt(0)*a[0] - cpt(1)*b[0] ;
-  c[1] = 0.25 - cpt(0)*a[1] - cpt(1)*b[1] ;
-  c[2] = 0.25 - cpt(0)*a[2] - cpt(1)*b[2] ;
-  c[3] = 0.25 - cpt(0)*a[3] - cpt(1)*b[3] ;
-
+  // linear shape with a half value: phi(x,y) = ax + by + c
+  std::vector<double> a(4), b(4), c(4);
+  get_linear_shape (cell, a, b, c);
 
 
   // compute basis functions
@@ -276,26 +302,26 @@ FE_P1NC::fill_fe_face_values (const Triangulation<2,2>::cell_iterator
   if (n_q_points==0)
     {
       Quadrature<2> cellquadrature = QProjector<2>::project_to_face(quadrature, face_no) ;
-      for(unsigned int i=0;i<cellquadrature.size();++i)
-               for (unsigned int k=0; k<this->dofs_per_cell; ++k)
-               {
-                   if (flags & update_values)
-                   {
-                   Point<2> realquadrature ;
-
-                   realquadrature = mapping.transform_unit_to_real_cell(cell, cellquadrature.point(i)) ;
-                   values[k] = a[k]*realquadrature(0) + b[k]*realquadrature(1) + c[k] ;
-                   output_data.shape_values[k][i] = values[k];
-                   }
-
-                   if (flags & update_gradients)
-                   {
-                   grads[k][0] = a[k] ;
-                   grads[k][1] = b[k] ;
-                   output_data.shape_gradients[k][i] = grads[k];
-                   }
-
-               }
+      for (unsigned int i=0; i<cellquadrature.size(); ++i)
+        for (unsigned int k=0; k<this->dofs_per_cell; ++k)
+          {
+            if (flags & update_values)
+              {
+                Point<2> realquadrature ;
+
+                realquadrature = mapping.transform_unit_to_real_cell(cell, cellquadrature.point(i)) ;
+                values[k] = a[k]*realquadrature(0) + b[k]*realquadrature(1) + c[k] ;
+                output_data.shape_values[k][i] = values[k];
+              }
+
+            if (flags & update_gradients)
+              {
+                grads[k][0] = a[k] ;
+                grads[k][1] = b[k] ;
+                output_data.shape_gradients[k][i] = grads[k];
+              }
+
+          }
 
     }
 
@@ -308,14 +334,14 @@ FE_P1NC::fill_fe_face_values (const Triangulation<2,2>::cell_iterator
 
 void
 FE_P1NC::fill_fe_subface_values (const Triangulation<2,2>::cell_iterator           &cell,
-                                              const unsigned int                                                  face_no,
-                                              const unsigned int                                                   sub_no,
-                                              const Quadrature<1>                                             &quadrature,
-                                              const Mapping<2,2>                                         &mapping,
-                                              const Mapping<2,2>::InternalDataBase              &mapping_internal,
-                                              const dealii::internal::FEValues::MappingRelatedData<2,2> &mapping_data,
-                                              const InternalDataBase                                              &fe_internal,
-                                              dealii::internal::FEValues::FiniteElementRelatedData<2,2> &output_data) const
+                                 const unsigned int                                                  face_no,
+                                 const unsigned int                                                   sub_no,
+                                 const Quadrature<1>                                             &quadrature,
+                                 const Mapping<2,2>                                         &mapping,
+                                 const Mapping<2,2>::InternalDataBase              &mapping_internal,
+                                 const dealii::internal::FEValues::MappingRelatedData<2,2> &mapping_data,
+                                 const InternalDataBase                                              &fe_internal,
+                                 dealii::internal::FEValues::FiniteElementRelatedData<2,2> &output_data) const
 
 {
   const UpdateFlags flags(fe_internal.update_each) ;
@@ -326,49 +352,9 @@ FE_P1NC::fill_fe_subface_values (const Triangulation<2,2>::cell_iterator
   std::vector<Tensor<1,2> > grads(flags & update_gradients ? this->dofs_per_cell : 0);
 
 
-
-  // edge midpoints
-  std::vector<Point<2> > mpt(4) ;
-
-  mpt[0](0) = (cell->vertex(0)(0) + cell->vertex(2)(0))/2.0 ;
-  mpt[0](1) = (cell->vertex(0)(1) + cell->vertex(2)(1))/2.0 ;
-
-  mpt[1](0) = (cell->vertex(1)(0) + cell->vertex(3)(0))/2.0 ;
-  mpt[1](1) = (cell->vertex(1)(1) + cell->vertex(3)(1))/2.0 ;
-
-  mpt[2](0) = (cell->vertex(0)(0) + cell->vertex(1)(0))/2.0 ;
-  mpt[2](1) = (cell->vertex(0)(1) + cell->vertex(1)(1))/2.0 ;
-
-  mpt[3](0) = (cell->vertex(2)(0) + cell->vertex(3)(0))/2.0 ;
-  mpt[3](1) = (cell->vertex(2)(1) + cell->vertex(3)(1))/2.0 ;
-
-  // center point
-  Point<2> cpt ;
-  cpt(0) = (mpt[0](0) + mpt[1](0) + mpt[2](0) + mpt[3](0))/4.0 ;
-  cpt(1) = (mpt[0](1) + mpt[1](1) + mpt[2](1) + mpt[3](1))/4.0 ;
-
-
-  // basis functions with a half value: phi(x,y) = ax + by + c
-  std::vector<double> a(4), b(4), c(4) ;
-  double det ;
-
-  det = (mpt[0](0)-mpt[1](0))*(mpt[2](1)-mpt[3](1)) - (mpt[2](0)-mpt[3](0))*(mpt[0](1)-mpt[1](1)) ;
-
-  a[0] = ((mpt[2](1)-mpt[3](1))*(0.5) -(mpt[0](1)-mpt[1](1))*(0.5))/det ;
-  a[1] = ((mpt[2](1)-mpt[3](1))*(-0.5) -(mpt[0](1)-mpt[1](1))*(0.5))/det ;
-  a[2] = ((mpt[2](1)-mpt[3](1))*(0.5) -(mpt[0](1)-mpt[1](1))*(-0.5))/det ;
-  a[3] = ((mpt[2](1)-mpt[3](1))*(-0.5) -(mpt[0](1)-mpt[1](1))*(-0.5))/det ;
-
-  b[0] = (-(mpt[2](0)-mpt[3](0))*(0.5) +(mpt[0](0)-mpt[1](0))*(0.5))/det ;
-  b[1] = (-(mpt[2](0)-mpt[3](0))*(-0.5) +(mpt[0](0)-mpt[1](0))*(0.5))/det ;
-  b[2] = (-(mpt[2](0)-mpt[3](0))*(0.5) +(mpt[0](0)-mpt[1](0))*(-0.5))/det ;
-  b[3] = (-(mpt[2](0)-mpt[3](0))*(-0.5) +(mpt[0](0)-mpt[1](0))*(-0.5))/det ;
-
-  c[0] = 0.25 - cpt(0)*a[0] - cpt(1)*b[0] ;
-  c[1] = 0.25 - cpt(0)*a[1] - cpt(1)*b[1] ;
-  c[2] = 0.25 - cpt(0)*a[2] - cpt(1)*b[2] ;
-  c[3] = 0.25 - cpt(0)*a[3] - cpt(1)*b[3] ;
-
+  // linear shape with a half value: phi(x,y) = ax + by + c
+  std::vector<double> a(4), b(4), c(4);
+  get_linear_shape (cell, a, b, c);
 
 
   // compute basis functions
@@ -400,88 +386,30 @@ FE_P1NC::fill_fe_subface_values (const Triangulation<2,2>::cell_iterator
   if (n_q_points==0)
     {
       Quadrature<2> cellquadrature = QProjector<2>::project_to_subface(quadrature, face_no, sub_no) ;
-      for(unsigned int i=0;i<cellquadrature.size();++i)
-               for (unsigned int k=0; k<this->dofs_per_cell; ++k)
-               {
-                   if (flags & update_values)
-                   {
-                   Point<2> realquadrature ;
-
-                   realquadrature = mapping.transform_unit_to_real_cell(cell, cellquadrature.point(i)) ;
-                   values[k] = a[k]*realquadrature(0) + b[k]*realquadrature(1) + c[k] ;
-                   output_data.shape_values[k][i] = values[k];
-                   }
-
-                   if (flags & update_gradients)
-                   {
-                   grads[k][0] = a[k] ;
-                   grads[k][1] = b[k] ;
-                   output_data.shape_gradients[k][i] = grads[k];
-                   }
-
-               }
-
-    }
-}
-
-
-
-FE_P1NC::FE_P1NC()
-  :
-  FiniteElement<2,2>(FiniteElementData<2>(get_dpo_vector(),
-                                          1,
-                                          1,
-                                          FiniteElementData<2>::L2),
-                     std::vector<bool> (1, false),
-                     get_nonzero_component())
-{
-
-  // face support points: 2 end vertices
-  unit_face_support_points.resize(2);
-  unit_face_support_points[0][0] = 0.0 ;
-  unit_face_support_points[1][0] = 1.0 ;
-
-
-  // initialize constraints matrix
-  initialize_constraints () ;
-}
-
-
-
-std::string FE_P1NC::get_name () const
-{
-  return "FE_P1NC";
-}
+      for (unsigned int i=0; i<cellquadrature.size(); ++i)
+        for (unsigned int k=0; k<this->dofs_per_cell; ++k)
+          {
+            if (flags & update_values)
+              {
+                Point<2> realquadrature ;
 
+                realquadrature = mapping.transform_unit_to_real_cell(cell, cellquadrature.point(i)) ;
+                values[k] = a[k]*realquadrature(0) + b[k]*realquadrature(1) + c[k] ;
+                output_data.shape_values[k][i] = values[k];
+              }
 
-UpdateFlags     FE_P1NC::requires_update_flags (const UpdateFlags flags) const
-{
-    UpdateFlags out = update_default;
-
-    if (flags & update_values)
-      out |= update_values | update_quadrature_points ;
-    if (flags & update_gradients)
-      out |= update_gradients ;
-    if (flags & update_cell_normal_vectors)
-      out |= update_cell_normal_vectors | update_JxW_values;
-    if (flags & update_hessians)
-      out |= update_hessians;
-    if (flags & update_hessians)
-      out |= update_hessians;
-  return out;
-}
+            if (flags & update_gradients)
+              {
+                grads[k][0] = a[k] ;
+                grads[k][1] = b[k] ;
+                output_data.shape_gradients[k][i] = grads[k];
+              }
 
+          }
 
-FiniteElement<2,2> *FE_P1NC::clone () const
-{
-  return new FE_P1NC(*this);
+    }
 }
 
-FE_P1NC::~FE_P1NC () {}
-
-
-
 
 void FE_P1NC::initialize_constraints ()
 {

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.