]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add a first version of the P1 non-conforming element.
authorJaeryun Yim <jaeryun.yim@gmail.com>
Thu, 5 Nov 2015 02:46:07 +0000 (11:46 +0900)
committerJaeryun Yim <jaeryun.yim@gmail.com>
Fri, 16 Sep 2016 12:07:19 +0000 (21:07 +0900)
This version currently only covers the parameteric version of the element.

include/deal.II/fe/fe_p1nc.h [new file with mode: 0644]
source/fe/CMakeLists.txt
source/fe/fe_p1nc.cc [new file with mode: 0644]

diff --git a/include/deal.II/fe/fe_p1nc.h b/include/deal.II/fe/fe_p1nc.h
new file mode 100644 (file)
index 0000000..d4d39fe
--- /dev/null
@@ -0,0 +1,156 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii__fe_p1nc_h
+#define dealii__fe_p1nc_h
+
+#include <deal.II/base/config.h>
+#include <deal.II/fe/fe.h>
+#include <deal.II/base/quadrature.h>
+#include <deal.II/base/qprojector.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+/*!@addtogroup fe */
+/*@{*/
+
+
+// parametric version of P1 Nonconforming FE
+class FE_P1NCParametric : public FiniteElement<2,2>
+{
+private:
+
+
+  static
+  std::vector<ComponentMask>
+  get_nonzero_component();
+
+
+  static
+  std::vector<unsigned int>
+  get_dpo_vector ();
+
+
+  class InternalData : public FiniteElement<2,2>::InternalDataBase
+  {
+  public:
+
+    std::vector<std::vector<double> > shape_values;
+    std::vector<std::vector<Tensor<1,2> > > shape_gradients;
+  };
+
+
+  virtual FiniteElement<2,2>::InternalDataBase *
+  get_data (const UpdateFlags update_flags,
+            const Mapping<2,2> &mapping,
+            const Quadrature<2> &quadrature) const ;
+
+
+
+
+  virtual
+  void
+  fill_fe_values (const Mapping<2,2>                                &mapping,
+                  const Triangulation<2,2>::cell_iterator           &cell,
+                  const Quadrature<2>                               &quadrature,
+                  const Mapping<2,2>::InternalDataBase              &mapping_internal,
+                  const FiniteElement<2,2>::InternalDataBase        &fe_internal,
+                  const internal::FEValues::MappingRelatedData<2,2> &mapping_data,
+                  internal::FEValues::FiniteElementRelatedData<2,2> &output_data,
+                  const CellSimilarity::Similarity                   cell_similarity) const;
+
+
+
+
+  virtual
+  void
+  fill_fe_face_values (const Mapping<2,2>                                &mapping,
+                       const Triangulation<2,2>::cell_iterator           &cell,
+                       const unsigned int                                 face_no,
+                       const Quadrature<1>                               &quadrature,
+                       const Mapping<2,2>::InternalDataBase              &mapping_internal,
+                       const FiniteElement<2,2>::InternalDataBase        &fe_internal,
+                       const internal::FEValues::MappingRelatedData<2,2> &mapping_data,
+                       internal::FEValues::FiniteElementRelatedData<2,2> &output_data) const;
+
+
+
+
+  virtual
+  void
+  fill_fe_subface_values (const Mapping<2,2>                                &mapping,
+                          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>::InternalDataBase              &mapping_internal,
+                          const FiniteElement<2,2>::InternalDataBase        &fe_internal,
+                          const internal::FEValues::MappingRelatedData<2,2> &mapping_data,
+                          internal::FEValues::FiniteElementRelatedData<2,2> &output_data) const;
+
+
+
+public:
+  FE_P1NCParametric() ;
+
+
+
+  virtual std::string get_name () const ;
+
+
+
+  virtual UpdateFlags     update_once (const UpdateFlags flags) const ;
+
+  virtual UpdateFlags     update_each (const UpdateFlags flags) const ;
+
+  const std::vector<Point<1> > &get_unit_face_support_points () const ;
+
+  Point<1> unit_face_support_point (const unsigned index) const ;
+
+  bool has_face_support_points () const ;
+
+  virtual FiniteElement<2,2> *clone () const ;
+
+  virtual ~FE_P1NCParametric ();
+
+
+
+
+
+  virtual double shape_value (const unsigned int i,
+                              const Point<2>    &p) const ;
+
+  virtual Tensor<1,2> shape_grad (const unsigned int i,
+                                  const Point<2>    &p) const ;
+
+
+
+protected:
+
+  void initialize_constraints () ;
+
+};
+
+
+
+
+
+/** @}*/
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index 0e8f549797cc141fa10f5a1c0cf4d98d90e38372..e0542b863b783713b4bc899498a7f497c2834be4 100644 (file)
@@ -33,6 +33,7 @@ SET(_src
   fe_nothing.cc
   fe_poly.cc
   fe_poly_tensor.cc
+  fe_p1nc.cc
   fe_q_base.cc
   fe_q.cc
   fe_q_bubbles.cc
diff --git a/source/fe/fe_p1nc.cc b/source/fe/fe_p1nc.cc
new file mode 100644 (file)
index 0000000..c20bdb5
--- /dev/null
@@ -0,0 +1,477 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+#include <deal.II/fe/fe_p1nc.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+// FE_P1NCParametric
+
+std::vector<ComponentMask>
+FE_P1NCParametric::get_nonzero_component()
+{
+  const unsigned int dofs_per_cell = 4;
+  std::vector<ComponentMask> masks (dofs_per_cell);
+  for (unsigned int i=0; i<dofs_per_cell; ++i)
+    masks[i] = ComponentMask(1,true);
+  return masks;
+}
+
+
+std::vector<unsigned int>
+FE_P1NCParametric::get_dpo_vector ()
+{
+  std::vector<unsigned int> dpo(3);
+  dpo[0] = 1; //dofs per object: vertex
+  dpo[1] = 0; // line
+  dpo[2] = 0; // quad
+  return dpo;
+}
+
+
+FiniteElement<2,2>::InternalDataBase *
+FE_P1NCParametric::get_data (const UpdateFlags update_flags,
+                             const Mapping<2,2> &mapping,
+                             const Quadrature<2> &quadrature) const
+{
+
+  InternalData *data = new InternalData;
+
+
+  data->update_once = update_once(update_flags);
+  data->update_each = update_each(update_flags);
+  data->update_flags = data->update_once | data->update_each;
+
+  const UpdateFlags flags(data->update_flags);
+  const unsigned int n_q_points = quadrature.size();
+
+
+  // initialize shape value
+  if (flags & update_values)
+    data->shape_values.resize (this->dofs_per_cell,
+                               std::vector<double> (n_q_points));
+
+  // initialize shape gradient
+  if (flags & update_gradients)
+    data->shape_gradients.resize (this->dofs_per_cell,
+                                  std::vector<Tensor<1,2> > (n_q_points));
+
+
+  // update
+  if (flags & (update_values | update_gradients))
+    for (unsigned int i=0; i<n_q_points; ++i)
+      {
+        if (flags & update_values)
+          for (unsigned int k=0; k<this->dofs_per_cell; ++k)
+            {
+              data->shape_values[k][i] = shape_value(k,quadrature.point(i));
+            }
+
+
+        if (flags & update_gradients)
+          for (unsigned int k=0; k<this->dofs_per_cell; ++k)
+            {
+              data->shape_gradients[k][i] = shape_grad(k,quadrature.point(i));
+            }
+
+      }
+  return data;
+
+}
+
+
+// get fe value in real cell from fedata in reference cell
+void
+FE_P1NCParametric::fill_fe_values (const Mapping<2,2>                                &mapping,
+                                   const Triangulation<2,2>::cell_iterator           &cell,
+                                   const Quadrature<2>                               &quadrature,
+                                   const Mapping<2,2>::InternalDataBase              &mapping_internal,
+                                   const FiniteElement<2,2>::InternalDataBase        &fe_internal,
+                                   const internal::FEValues::MappingRelatedData<2,2> &mapping_data,
+                                   internal::FEValues::FiniteElementRelatedData<2,2> &output_data,
+                                   const CellSimilarity::Similarity                   cell_similarity) const
+
+{
+
+  Assert (dynamic_cast<const InternalData *> (&fe_internal) != 0, ExcInternalError());
+  const InternalData &fe_data = static_cast<const InternalData &> (fe_internal);
+
+  const UpdateFlags flags(fe_data.current_update_flags());
+
+  // gets shape function values from its input argument fedata
+  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
+    {
+      // shape value
+      if (flags & update_values)
+        for (unsigned int i=0; i<quadrature.size(); ++i)
+          output_data.shape_values(k,i) = fe_data.shape_values[k][i];
+
+
+      // shape gradient
+      if (flags & update_gradients && cell_similarity != CellSimilarity::translation)
+        mapping.transform(fe_data.shape_gradients[k],
+                          mapping_covariant,
+                          mapping_internal,
+                          output_data.shape_gradients[k]);
+
+    }
+
+}
+
+
+void
+FE_P1NCParametric::fill_fe_face_values (const Mapping<2,2>                                &mapping,
+                                        const Triangulation<2,2>::cell_iterator           &cell,
+                                        const unsigned int                                 face_no,
+                                        const Quadrature<1>                               &quadrature,
+                                        const Mapping<2,2>::InternalDataBase              &mapping_internal,
+                                        const FiniteElement<2,2>::InternalDataBase        &fe_internal,
+                                        const internal::FEValues::MappingRelatedData<2,2> &mapping_data,
+                                        internal::FEValues::FiniteElementRelatedData<2,2> &output_data) const
+
+
+{
+
+  Assert (dynamic_cast<const InternalData *> (&fe_internal) != 0, ExcInternalError());
+  const InternalData &fe_data = static_cast<const InternalData &> (fe_internal);
+
+
+  const QProjector<2>::DataSetDescriptor offset
+    = QProjector<2>::DataSetDescriptor::face (face_no,
+                                              cell->face_orientation(face_no),
+                                              cell->face_flip(face_no),
+                                              cell->face_rotation(face_no),
+                                              quadrature.size());
+
+  const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
+
+  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
+    {
+      if (flags & update_values)
+        for (unsigned int i=0; i<quadrature.size(); ++i)
+          output_data.shape_values(k,i) = fe_data.shape_values[k][i+offset];
+
+      if (flags & update_gradients)
+        mapping.transform(make_slice(fe_data.shape_gradients[k], offset, quadrature.size()),
+                          mapping_covariant,
+                          mapping_internal,
+                          output_data.shape_gradients[k]);
+    }
+
+
+}
+
+
+
+
+
+
+void
+FE_P1NCParametric::fill_fe_subface_values (const Mapping<2,2>                                &mapping,
+                                           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>::InternalDataBase              &mapping_internal,
+                                           const FiniteElement<2,2>::InternalDataBase        &fe_internal,
+                                           const internal::FEValues::MappingRelatedData<2,2> &mapping_data,
+                                           internal::FEValues::FiniteElementRelatedData<2,2> &output_data) const
+
+{
+
+  Assert (dynamic_cast<const InternalData *> (&fe_internal) != 0, ExcInternalError());
+  const InternalData &fe_data = static_cast<const InternalData &> (fe_internal);
+
+
+
+  const QProjector<2>::DataSetDescriptor offset
+    = QProjector<2>::DataSetDescriptor::subface (face_no, sub_no,
+                                                 cell->face_orientation(face_no),
+                                                 cell->face_flip(face_no),
+                                                 cell->face_rotation(face_no),
+                                                 quadrature.size(),
+                                                 cell->subface_case(face_no));
+
+  const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
+
+  for (unsigned int k=0; k<this->dofs_per_cell; ++k)
+    {
+      if (flags & update_values)
+        for (unsigned int i=0; i<quadrature.size(); ++i)
+          output_data.shape_values(k,i) = fe_data.shape_values[k][i+offset];
+
+      if (flags & update_gradients)
+        mapping.transform(make_slice(fe_data.shape_gradients[k], offset, quadrature.size()),
+                          mapping_covariant,
+                          mapping_internal,
+                          output_data.shape_gradients[k]);
+    }
+
+}
+
+
+
+FE_P1NCParametric::FE_P1NCParametric()
+  :
+  FiniteElement<2,2>(FiniteElementData<2>(get_dpo_vector(),
+                                          1,
+                                          1,
+                                          FiniteElementData<2>::L2,
+                                          numbers::invalid_unsigned_int),
+                     std::vector<bool> (1, false),
+                     get_nonzero_component())
+{
+
+  // support points: 4 vertices
+  unit_support_points.resize(4) ;
+  unit_support_points[0][0] = 0.0 ;
+  unit_support_points[0][1] = 0.0 ;
+  unit_support_points[1][0] = 1.0 ;
+  unit_support_points[1][1] = 0.0 ;
+  unit_support_points[2][0] = 0.0 ;
+  unit_support_points[2][1] = 1.0 ;
+  unit_support_points[3][0] = 1.0 ;
+  unit_support_points[3][1] = 1.0 ;
+
+  // 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 ;
+
+
+  // CONSTRAINTS MATRIX INITIALIZATION
+  initialize_constraints () ;
+}
+
+
+
+std::string FE_P1NCParametric::get_name () const
+{
+  return "FE_P1NCParametric";
+}
+
+
+
+UpdateFlags     FE_P1NCParametric::update_once (const UpdateFlags flags) const
+{
+  return (update_default | (flags & update_values));
+}
+
+UpdateFlags     FE_P1NCParametric::update_each (const UpdateFlags flags) const
+{
+  UpdateFlags out = update_default;
+
+  if (flags & update_gradients)
+    out |= update_gradients | update_covariant_transformation;
+  if (flags & update_cell_normal_vectors)
+    out |= update_cell_normal_vectors | update_JxW_values;
+
+  return out;
+}
+
+
+const std::vector<Point<1> > &FE_P1NCParametric::get_unit_face_support_points () const
+{
+
+  Assert ((unit_face_support_points.size() == 0) ||
+          (unit_face_support_points.size() == this->dofs_per_face),
+          ExcInternalError());
+  return unit_face_support_points;
+}
+
+Point<1> FE_P1NCParametric::unit_face_support_point (const unsigned index) const
+{
+  Assert (index < this->dofs_per_face,
+          ExcIndexRange (index, 0, this->dofs_per_face));
+  Assert (unit_face_support_points.size() == this->dofs_per_face,
+          ExcFEHasNoSupportPoints ());
+  return unit_face_support_points[index];
+}
+
+bool FE_P1NCParametric::has_face_support_points () const
+{
+  return (unit_face_support_points.size() != 0);
+}
+
+FiniteElement<2,2> *FE_P1NCParametric::clone () const
+{
+  return new FE_P1NCParametric(*this);
+}
+
+FE_P1NCParametric::~FE_P1NCParametric () {}
+
+// Define the exact value of the shape function for the P1NC scalar
+// function in the unit cell = (0,1)^2.
+//
+//  2-----------------|------------------3
+//  |                                    |
+//  |                                    |
+//  |                                    |
+//  |                                    |
+//  |                                    |
+//  |                                    |
+//  |                                    |
+//  |                                    |
+//  |                                    |
+//  -                                    -
+//  |                                    |
+//  |                                    |
+//  |                                    |
+//  |                                    |
+//  |                                    |
+//  |                                    |
+//  |                                    |
+//  |                                    |
+//  |                                    |
+//  0-----------------|------------------1
+//
+
+
+double FE_P1NCParametric::shape_value (const unsigned int i,
+                                       const Point<2>    &p) const
+{
+
+
+  double value = 0.0 ;
+
+  // P1NC with HALF COEFFICIENT
+  switch (i)
+    {
+    case 0:
+      value = -p[0]-p[1]+1.5 ;
+      break ;
+    case 1:
+      value =  p[0]-p[1]+0.5 ;
+      break ;
+    case 2:
+      value = -p[0]+p[1]+0.5 ;
+      break ;
+    case 3:
+      value =  p[0]+p[1]-0.5 ;
+      break ;
+    default:
+      value = 0.0 ;
+    };
+
+  return value*0.5 ;
+
+
+  // // Q1 for test
+  // switch (i)
+  //  {
+  //  case 0:
+  //    return  (1.0-p[0])*(1.0-p[1]) ;
+  //  case 1:
+  //    return  p[0]*(1.0-p[1]) ;
+  //  case 2:
+  //    return  (1.0-p[0])*p[1] ;
+  //  case 3:
+  //    return  p[0]*p[1] ;
+  //  default:
+  //    return 0 ;
+  //  };
+
+
+}
+
+
+
+
+Tensor<1,2> FE_P1NCParametric::shape_grad (const unsigned int i,
+                                           const Point<2>    &p) const
+{
+  Tensor<1,2> grad ;
+
+
+  // P1NC with HALF COEFFICIENT
+  switch (i)
+    {
+    case 0:
+      grad[0] = -1.0 ;
+      grad[1] = -1.0 ;
+      break ;
+    case 1:
+      grad[0] =  1.0 ;
+      grad[1] = -1.0 ;
+      break ;
+    case 2:
+      grad[0] = -1.0 ;
+      grad[1] =  1.0 ;
+      break ;
+    case 3:
+      grad[0] =  1.0 ;
+      grad[1] =  1.0 ;
+      break ;
+    default:
+      grad = 0.0 ;
+    };
+
+  return grad*0.5 ;
+
+
+  // // Q1 for test
+  // switch (i)
+  //  {
+  //  case 0:
+  //    grad[0] = p[1]-1.0 ;
+  //    grad[1] = p[0]-1.0 ;
+  //    return grad ;
+  //  case 1:
+  //    grad[0] = 1.0-p[1] ;
+  //    grad[1] = -p[0] ;
+  //    return grad ;
+  //  case 2:
+  //    grad[0] = -p[1] ;
+  //    grad[1] = 1.0-p[0] ;
+  //    return grad ;
+  //  case 3:
+  //    grad[0] =  p[1] ;
+  //    grad[1] =  p[0] ;
+  //    return grad ;
+  //  default:
+  //    return grad ;
+  //  };
+
+
+}
+
+
+// 2015 05 13 CONSTRAINTS MATRIX FOR HANGING NODE
+// BASED ON 'fe_q_base.cc'
+
+void FE_P1NCParametric::initialize_constraints ()
+{
+
+  std::vector<Point<1> > constraint_points;
+  // Add midpoint
+  constraint_points.push_back (Point<1> (0.5));
+
+  // Now construct relation between destination (child) and source (mother)
+  // dofs.
+
+  interface_constraints
+  .TableBase<2,double>::reinit (interface_constraints_size());
+
+
+  interface_constraints(0,0) = 0.5 ;
+  interface_constraints(0,1) = 0.5 ;
+
+}
+
+
+
+DEAL_II_NAMESPACE_CLOSE

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.