]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
the new trace element
authorangela.klewinghaus <angela.klewinghaus@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 25 Feb 2013 16:08:49 +0000 (16:08 +0000)
committerangela.klewinghaus <angela.klewinghaus@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 25 Feb 2013 16:08:49 +0000 (16:08 +0000)
git-svn-id: https://svn.dealii.org/branches/branch_face_elements@28556 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/source/fe/fe_trace.cc [new file with mode: 0644]
deal.II/source/fe/fe_trace.inst.in [new file with mode: 0644]

diff --git a/deal.II/source/fe/fe_trace.cc b/deal.II/source/fe/fe_trace.cc
new file mode 100644 (file)
index 0000000..390e0e0
--- /dev/null
@@ -0,0 +1,247 @@
+#include <deal.II/grid/tria.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/dofs/dof_accessor.h>
+
+#include <deal.II/base/config.h>
+#include <deal.II/base/tensor_product_polynomials.h>
+#include <deal.II/fe/fe_poly_face.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/data_out_faces.h>
+#include <fstream>
+#include <iostream>
+
+//using namespace dealii;
+
+#include <deal.II/fe/fe_face.h>
+#include <deal.II/fe/fe_poly_face.templates.h>
+#include <sstream>
+
+DEAL_II_NAMESPACE_OPEN
+
+/** FE-Class for continuous face functions
+* adaption of FE_Q and FE_FaceQ
+*/
+
+template <int dim, int spacedim=dim>
+class FE_TraceQ : public FE_PolyFace<TensorProductPolynomials<dim-1>, dim, spacedim>
+{
+  public:
+                                     /**
+                                      * Constructor for tensor product
+                                      * polynomials of degree
+                                      * <tt>p</tt>. The shape
+                                      * functions created using this
+                                      * constructor correspond to
+                                      * Legendre polynomials in each
+                                      * coordinate direction.
+                                      */
+    FE_TraceQ(unsigned int p);
+
+    virtual FiniteElement<dim,spacedim>* clone() const;
+
+                                     /**
+                                      * Return a string that uniquely
+                                      * identifies a finite
+                                      * element. This class returns
+                                      * <tt>FE_DGQ<dim>(degree)</tt> , with
+                                      * <tt>dim</tt> and <tt>degree</tt>
+                                      * replaced by appropriate
+                                      * values.
+                                      */
+    virtual std::string get_name () const;
+
+                                     /**
+                                      * Check for non-zero values on a face.
+                                      *
+                                      * This function returns
+                                      * @p true, if the shape
+                                      * function @p shape_index has
+                                      * non-zero values on the face
+                                      * @p face_index.
+                                      *
+                                      * Implementation of the
+                                      * interface in
+                                      * FiniteElement
+                                      */
+    virtual bool has_support_on_face (const unsigned int shape_index,
+                                      const unsigned int face_index) const;
+
+  private:
+                                     /**
+                                      * Return vector with dofs per
+                                      * vertex, line, quad, hex.
+                                      */
+    static std::vector<unsigned int> get_dpo_vector (const unsigned int deg);
+};
+
+
+template <int dim, int spacedim>
+FE_TraceQ<dim,spacedim>::FE_TraceQ (const unsigned int degree)
+                :
+                FE_PolyFace<TensorProductPolynomials<dim-1>, dim, spacedim> (
+                  TensorProductPolynomials<dim-1>(Polynomials::LagrangeEquidistant::generate_complete_basis(degree)),
+                  FiniteElementData<dim>(get_dpo_vector(degree), 1, degree, FiniteElementData<dim>::L2),
+                  std::vector<bool>(1,true))
+{  
+  std::vector<unsigned int> renumber (this->dofs_per_face);
+  FETools::hierarchic_to_lexicographic_numbering<dim-1> (degree, renumber);
+  this->poly_space.set_numbering(renumber);
+}
+
+
+template <int dim, int spacedim>
+FiniteElement<dim,spacedim>*
+FE_TraceQ<dim,spacedim>::clone() const
+{
+  return new FE_TraceQ<dim,spacedim>(this->degree);
+}
+
+
+template <int dim, int spacedim>
+std::string
+FE_TraceQ<dim,spacedim>::get_name () const
+{
+                                   // note that the
+                                   // FETools::get_fe_from_name
+                                   // function depends on the
+                                   // particular format of the string
+                                   // this function returns, so they
+                                   // have to be kept in synch
+
+  std::ostringstream namebuf;
+  namebuf << "FE_TraceQ<" << dim << ">(" << this->degree << ")";
+
+  return namebuf.str();
+}
+
+template <int dim, int spacedim>
+bool
+FE_TraceQ<dim,spacedim>::has_support_on_face (const unsigned int shape_index,
+                                const unsigned int face_index) const
+{
+  Assert (shape_index < this->dofs_per_cell,
+          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));
+
+                                   // let's see whether this is a
+                                   // vertex
+  if (shape_index < this->first_line_index)
+    {
+                                       // for Q elements, there is
+                                       // one dof per vertex, so
+                                       // shape_index==vertex_number. check
+                                       // whether this vertex is
+                                       // on the given face. thus,
+                                       // for each face, give a
+                                       // list of vertices
+      const unsigned int vertex_no = shape_index;
+      Assert (vertex_no < GeometryInfo<dim>::vertices_per_cell,
+              ExcInternalError());
+
+      for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
+        if (GeometryInfo<dim>::face_to_cell_vertices(face_index, v) == vertex_no)
+          return true;
+
+      return false;
+    }
+  else if (shape_index < this->first_quad_index)
+                                     // ok, dof is on a line
+    {
+      const unsigned int line_index
+        = (shape_index - this->first_line_index) / this->dofs_per_line;
+      Assert (line_index < GeometryInfo<dim>::lines_per_cell,
+              ExcInternalError());
+
+                                       // in 2d, the line is the
+                                       // face, so get the line
+                                       // index
+      if (dim == 2)
+        return (line_index == face_index);
+      else if (dim == 3)
+        {
+                                           // see whether the
+                                           // given line is on the
+                                           // given face.
+          for (unsigned int l=0; l<GeometryInfo<dim>::lines_per_face; ++l)
+            if (GeometryInfo<3>::face_to_cell_lines(face_index, l) == line_index)
+              return true;
+
+          return false;
+        }
+      else
+        Assert (false, ExcNotImplemented());
+    }
+  else if (shape_index < this->first_hex_index)
+                                     // dof is on a quad
+    {
+      const unsigned int quad_index
+        = (shape_index - this->first_quad_index) / this->dofs_per_quad;
+      Assert (static_cast<signed int>(quad_index) <
+              static_cast<signed int>(GeometryInfo<dim>::quads_per_cell),
+              ExcInternalError());
+
+                                       // in 2d, cell bubble are
+                                       // zero on all faces. but
+                                       // we have treated this
+                                       // case above already
+      Assert (dim != 2, ExcInternalError());
+
+                                       // in 3d,
+                                       // quad_index=face_index
+      if (dim == 3)
+        return (quad_index == face_index);
+      else
+        Assert (false, ExcNotImplemented());
+    }
+  else
+                                     // dof on hex
+    {
+                                       // can only happen in 3d,
+                                       // but this case has
+                                       // already been covered
+                                       // above
+      Assert (false, ExcNotImplemented());
+      return false;
+    }
+
+                                   // we should not have gotten here
+  Assert (false, ExcInternalError());
+  return false;
+}
+
+template <int dim, int spacedim>
+std::vector<unsigned int>
+FE_TraceQ<dim,spacedim>::get_dpo_vector (const unsigned int deg)
+{
+  std::vector<unsigned int> dpo(dim+1, 1U);
+  dpo[dim]=0;
+  dpo[0]=1;
+  for (unsigned int i=1;i<dim;++i)
+    dpo[i] = dpo[i-1]*(deg-1);
+  return dpo;
+}
+
+// explicit instantiations
+#include "fe_trace.inst"
+
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/deal.II/source/fe/fe_trace.inst.in b/deal.II/source/fe/fe_trace.inst.in
new file mode 100644 (file)
index 0000000..22e13a6
--- /dev/null
@@ -0,0 +1,21 @@
+//---------------------------------------------------------------------------
+//    $Id: fe_face.inst.in 22676 2010-11-10 15:21:40Z kronbichler $
+//    Version: $Name$
+//
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+
+for (deal_II_dimension : DIMENSIONS)
+  {
+#if deal_II_dimension > 1
+    template class FE_TraceQ<deal_II_dimension>;
+#endif
+  }
+

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.