]> https://gitweb.dealii.org/ - dealii.git/commitdiff
mapping merge
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Tue, 6 Mar 2001 12:20:49 +0000 (12:20 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Tue, 6 Mar 2001 12:20:49 +0000 (12:20 +0000)
git-svn-id: https://svn.dealii.org/trunk@4119 0785d39b-7218-0410-832d-ea1e28bc413d

tests/fe/Makefile.in
tests/fe/fe_data_test.cc [new file with mode: 0644]
tests/fe/mapping.cc [new file with mode: 0644]
tests/fe/shapes.cc [new file with mode: 0644]

index b44d0abb69d8663101e9e9caece37518afa5048f..7f60aa37bf12dc163b4a500ee08d286f3e6c676d 100644 (file)
@@ -89,6 +89,7 @@ fe_data_test-o-files = $(fe_data_test-cc-files:.cc=.o)
 endif
 
 fe_data_test: $(fe_data_test-o-files) $(libraries)
+       @echo =====linking======= $<
        @$(CXX) $(flags) -o $@ $^
 
 ############################################################
@@ -102,6 +103,7 @@ shapes-o-files = $(shapes-cc-files:.cc=.o)
 endif
 
 shapes: $(shapes-o-files) $(libraries)
+       @echo =====linking======= $<
        @$(CXX) $(flags) -o $@ $^
 
 ############################################################
@@ -115,6 +117,7 @@ derivatives-o-files = $(derivatives-cc-files:.cc=.o)
 endif
 
 derivatives: $(derivatives-o-files) $(libraries)
+       @echo =====linking======= $<
        @$(CXX) $(flags) -o $@ $^
 
 ############################################################
@@ -128,6 +131,7 @@ embedding-o-files = $(embedding-cc-files:.cc=.o)
 endif
 
 embedding: $(embedding-o-files) $(libraries)
+       @echo =====linking======= $<
        @$(CXX) $(flags) -o $@ $^
 
 ############################################################
@@ -154,6 +158,7 @@ mapping-o-files = $(mapping-cc-files:.cc=.o)
 endif
 
 mapping: $(mapping-o-files) $(libraries)
+       @echo =====linking======= $<
        @$(CXX) $(flags) -o $@ $^
 
 ############################################################
@@ -167,6 +172,7 @@ internals-o-files = $(internals-cc-files:.cc=.o)
 endif
 
 internals: $(internals-o-files) $(libraries)
+       @echo =====linking======= $<
        @$(CXX) $(flags) -o $@ $^
 
 ############################################################
@@ -180,6 +186,7 @@ performance-o-files = $(performance-cc-files:.cc=.o)
 endif
 
 performance: $(performance-o-files) $(libraries)
+       @echo =====linking======= $<
        @$(CXX) $(flags) -o $@ $^
 
 ############################################################
diff --git a/tests/fe/fe_data_test.cc b/tests/fe/fe_data_test.cc
new file mode 100644 (file)
index 0000000..cd7c381
--- /dev/null
@@ -0,0 +1,78 @@
+//----------------------------  fe_data_test.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1998, 1999, 2000, 2001 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.
+//
+//----------------------------  fe_data_test.cc  ---------------------------
+
+
+#include <iostream>
+#include <fstream>
+
+#include <base/logstream.h>
+#include <grid/tria_iterator.h>
+#include <dofs/dof_accessor.h>
+#include <fe/fe_q.h>
+#include <fe/fe_system.h>
+
+
+template<int dim>
+void test_fe_datas()
+{
+  vector<FiniteElement<dim> *> fe_datas;
+  fe_datas.push_back(new FE_Q<dim> (1));
+  fe_datas.push_back(new FE_Q<dim> (4));
+  fe_datas.push_back(new FESystem<dim>(FE_Q<dim> (2), 2));
+  fe_datas.push_back(new FESystem<dim>(FE_Q<dim> (1), 2,
+                                      FE_Q<dim> (2), 1));
+                                  // for this an assertion thrown in
+                                  // FESystem::build_interface_constraints
+                                  // for the following calls for dim=3
+//  fe_datas.push_back(new FESystem<dim>(FE_Q<dim> (3), 2));
+//  fe_datas.push_back(new FESystem<dim>(FE_Q<dim> (1), 2,
+//                                    FE_Q<dim> (3), 1));
+//  fe_datas.push_back(new FESystem<dim>(FE_Q<dim> (4), 2));
+  
+  deallog << endl << "dim=" << dim << endl;
+  for (unsigned int n=0; n<fe_datas.size(); ++n)
+    {
+      FiniteElementData<dim> *fe_data=fe_datas[n];
+      deallog << "fe_data[" << n <<"]:" << endl;
+      deallog << "dofs_per_vertex=" << fe_data->dofs_per_vertex << endl;
+      deallog << "dofs_per_line=" << fe_data->dofs_per_line << endl;
+      deallog << "dofs_per_quad=" << fe_data->dofs_per_quad << endl;
+      deallog << "dofs_per_hex=" << fe_data->dofs_per_hex << endl;
+      deallog << "first_line_index=" << fe_data->first_line_index << endl;
+      deallog << "first_quad_index=" << fe_data->first_quad_index << endl;
+      deallog << "first_hex_index=" << fe_data->first_hex_index << endl;
+      deallog << "first_face_line_index=" << fe_data->first_face_line_index << endl;
+      deallog << "first_face_quad_index=" << fe_data->first_face_quad_index << endl;
+      deallog << "dofs_per_face=" << fe_data->dofs_per_face << endl;
+      deallog << "dofs_per_cell=" << fe_data->dofs_per_cell << endl;
+      deallog << "components=" << fe_data->components << endl;
+    }
+
+                                  // delete all FiniteElementDatas
+  for (unsigned int i=0; i<fe_datas.size(); ++i)
+    delete fe_datas[i];
+}
+
+int main(int,char)
+{
+  ofstream logfile("fe_data_test.dat");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  test_fe_datas<1>();
+  test_fe_datas<2>();
+  test_fe_datas<3>();
+}
+
+
+
diff --git a/tests/fe/mapping.cc b/tests/fe/mapping.cc
new file mode 100644 (file)
index 0000000..ff693b2
--- /dev/null
@@ -0,0 +1,502 @@
+// $Id$
+// Copyright (C) 2001 Ralf Hartmann
+//
+// Show the shape functions implemented and computes the area of cells.
+
+#include <base/quadrature_lib.h>
+#include <base/logstream.h>
+#include <lac/vector.h>
+#include <grid/tria.h>
+#include <grid/tria_iterator.h>
+#include <dofs/dof_accessor.h>
+#include <grid/grid_generator.h>
+#include <grid/grid_out.h>
+#include <grid/tria_boundary_lib.h>
+#include <fe/mapping_q1.h>
+#include <fe/mapping_q.h>
+#include <fe/fe_q.h>
+#include <fe/fe_values.h>
+#include <fe/fe.h>
+#include <vector>
+#include <fstream>
+#include <string>
+#include <strstream>
+
+
+
+char fname[50];
+
+template<int dim>
+inline void
+plot_transformation(Mapping<dim> &mapping,
+                   FiniteElement<dim> &fe,
+                   DoFHandler<dim>::cell_iterator &cell,
+                   const char* name)
+{
+  const unsigned int div = 7;
+
+  QTrapez<1> q_trapez;
+  QIterated<dim> q(q_trapez, div);
+  FEValues<dim> fe_values(mapping, fe, q,
+                         UpdateFlags(update_q_points
+                                     | update_JxW_values));
+
+  fe_values.reinit(cell);
+  const vector<double> &JxW=fe_values.get_JxW_values();
+  
+  ofstream gnuplot(name);
+
+  unsigned int k=0;
+  for (unsigned int nz=0; nz<=((dim>2) ? div : 0); ++nz)
+    {
+      for (unsigned int ny=0; ny<=((dim>1) ? div : 0); ++ny)
+       {
+         for (unsigned int nx=0; nx<=div; ++nx)
+           {
+             gnuplot << fe_values.quadrature_point(k);
+             double J = JxW[k] / q.weight(k);
+             gnuplot << ' ' << J << endl;
+             ++k;
+           }
+         gnuplot << endl;
+       }
+      gnuplot << endl;  
+    }
+}
+
+
+
+template<int dim>
+inline void
+plot_faces(Mapping<dim> &mapping,
+          FiniteElement<dim> &fe,
+          DoFHandler<dim>::cell_iterator &cell,
+          const char* name)
+{
+  ofstream gnuplot(name);
+
+  QGauss4<dim-1> q;
+  const unsigned int nq = (unsigned int) (.01 + pow(q.n_quadrature_points, 1./(dim-1)));
+
+  FEFaceValues<dim> fe_values(mapping, fe, q,
+                             UpdateFlags(update_q_points
+                                         | update_normal_vectors));
+
+  for (unsigned int face_nr=0;
+       face_nr < GeometryInfo<dim>::faces_per_cell;
+       ++ face_nr)
+    {
+      fe_values.reinit(cell, face_nr);
+
+      const vector<Point<dim> > &normals
+       =fe_values.get_normal_vectors();
+
+      unsigned int k=0;
+      for (unsigned int ny=0; ny<((dim>2) ? nq : 1); ++ny)
+       {
+         for (unsigned int nx=0; nx<nq; ++nx)
+           {
+             Point<dim> x = fe_values.quadrature_point(k);
+             Tensor<1,dim> n = normals[k];
+             gnuplot << x << '\t' << n << endl;
+             ++k;
+           }
+         gnuplot << endl;
+       }
+      gnuplot << endl;      
+    }
+}
+
+
+
+template<int dim>
+inline void
+plot_subfaces(Mapping<dim> &mapping,
+             FiniteElement<dim> &fe,
+             DoFHandler<dim>::cell_iterator &cell,
+             const char* name)
+{
+  ofstream gnuplot(name);
+
+  QGauss4<dim-1> q;
+  const unsigned int nq = (unsigned int) (.01 + pow(q.n_quadrature_points, 1./(dim-1)));
+
+  FESubfaceValues<dim> fe_values(mapping, fe, q,
+                                UpdateFlags(update_q_points
+                                            | update_normal_vectors));
+  for (unsigned int face_nr=0;
+       face_nr < GeometryInfo<dim>::faces_per_cell;
+       ++ face_nr)
+    for (unsigned int sub_nr=0;
+        sub_nr < GeometryInfo<dim>::subfaces_per_face;
+        ++ sub_nr)
+      {
+       fe_values.reinit(cell, face_nr, sub_nr);
+       
+       const vector<Point<dim> > &normals
+         =fe_values.get_normal_vectors();
+       
+       unsigned int k=0;
+       for (unsigned int ny=0; ny<((dim>2) ? nq : 1); ++ny)
+         {
+           for (unsigned int nx=0; nx<nq; ++nx)
+             {
+               Point<dim> x = fe_values.quadrature_point(k);
+               Tensor<1,dim> n = normals[k];
+               gnuplot << x << '\t' << n << endl;
+               ++k;
+             }
+           gnuplot << endl;
+         }
+       gnuplot << endl;
+      }
+}
+
+
+
+template<>
+inline void
+plot_faces(Mapping<1>&,
+          FiniteElement<1>&,
+          DoFHandler<1>::cell_iterator&,
+          const char*)
+{};
+
+
+
+template<>
+inline void
+plot_subfaces(Mapping<1>&,
+             FiniteElement<1>&,
+             DoFHandler<1>::cell_iterator&,
+             const char*)
+{};
+
+
+
+template<int dim>
+inline void
+compute_area(Mapping<dim> &mapping,
+            FiniteElement<dim> &fe,
+            DoFHandler<dim>::cell_iterator &cell)
+{
+  QGauss4<dim> gauss4;
+  FEValues<dim> fe_values(mapping, fe, gauss4,
+                       UpdateFlags(update_JxW_values));
+  fe_values.reinit(cell);
+  const vector<double> &JxW=fe_values.get_JxW_values();
+
+  double area=0;
+  for (unsigned int i=0; i<fe_values.n_quadrature_points; ++i)
+    area+=JxW[i];
+  deallog << "  area=" << area << endl;
+}
+
+
+template<int dim>
+void create_triangulations(vector<Triangulation<dim> *> &,
+                          vector<Boundary<dim> *> &,
+                          vector<double> &)
+{
+  Assert(false, ExcNotImplemented());
+}
+
+
+
+vector<vector<unsigned int> > show;
+unsigned int mapping_size;
+
+
+template<>
+void create_triangulations(vector<Triangulation<1> *> &tria_ptr,
+                          vector<Boundary<1> *> &,
+                          vector<double> &exact_areas)
+{
+  show.resize(1, vector<unsigned int> (mapping_size,0));
+  Triangulation<1> *tria=new Triangulation<1>();
+  tria_ptr.push_back(tria);
+  GridGenerator::hyper_cube(*tria, 1., 3.);
+  exact_areas.push_back(2.);
+  show[0][0]=1;
+  show[0][1]=1;
+}
+
+
+
+template<>
+void create_triangulations(vector<Triangulation<2> *> &tria_ptr,
+                          vector<Boundary<2> *> &boundary_ptr,
+                          vector<double> &exact_areas)
+{
+  Triangulation<2> *tria;
+  show.clear();
+  show.resize(3, vector<unsigned int> (mapping_size,0));
+                                  // tria0: 3x3 square rotated
+  if (1)
+    {
+      tria=new Triangulation<2>();
+      tria_ptr.push_back(tria);
+      const double left = 1.;
+      const double right = 4.;
+      
+      const Point<2> vertices[4] = { Point<2>(left,left),
+                                      Point<2>(right,left),
+                                      Point<2>(right,right),
+                                      Point<2>(left,right)  };
+      const int cell_vertices[1][4] = { { 1,2,3,0 } };
+      std::vector<CellData<2> > cells (1, CellData<2>());
+      for (unsigned int j=0; j<4; ++j)
+       cells[0].vertices[j] = cell_vertices[0][j];
+      cells[0].material_id = 0;
+      
+      tria->create_triangulation (std::vector<Point<2> >(&vertices[0], &vertices[4]),
+                                cells,
+                                SubCellData());
+      exact_areas.push_back(9.);
+    }
+  
+                                  // tria1: arbitrary rectangle
+  if (1)
+    {
+      tria=new Triangulation<2>();
+      tria_ptr.push_back(tria);
+      GridGenerator::hyper_cube(*tria, 1., 3.);
+      Point<2> &v=tria->begin_quad()->vertex(2);
+      v(0) = 5.;
+      v(1) = 4.;
+      exact_areas.push_back(7.);
+    }
+  
+                                  // tria2: crazy cell
+  if (1)
+    {
+      Boundary<2> *boundary1=new HyperBallBoundary<2>(Point<2>(3,1), 2);
+      Boundary<2> *boundary2=new HyperBallBoundary<2>(Point<2>(2,5), sqrt(5));
+      boundary_ptr.push_back(boundary1);
+      boundary_ptr.push_back(boundary2);      
+      tria=new Triangulation<2>();
+      tria_ptr.push_back(tria);
+      GridGenerator::hyper_cube(*tria, 1., 5.);
+      Point<2> &v2 = tria->begin_active()->vertex(2);
+      Point<2> &v3 = tria->begin_active()->vertex(3);
+      v2(0) = 3.;
+      v2(1) = 3.;
+      v3(0) = 1.;
+      v3(1) = 3.;
+      tria->set_boundary(1,*boundary1);
+      tria->set_boundary(2,*boundary2);
+      tria->begin_active()->face(1)->set_boundary_indicator(1);
+      tria->begin_active()->face(2)->set_boundary_indicator(2);
+      double pi=acos(-1);
+      double alpha=2*atan(0.5);
+      exact_areas.push_back(4+pi-2.5*(alpha-sin(alpha)));
+      for (unsigned int i=0; i<mapping_size; ++i)
+       show[2][i]=1;
+    }
+}
+
+
+
+template<>
+void create_triangulations(vector<Triangulation<3> *> &tria_ptr,
+                          vector<Boundary<3> *> &boundary_ptr,
+                          vector<double> &exact_areas)
+{
+  Triangulation<3> *tria;
+  show.clear();
+  show.resize(4, vector<unsigned int> (mapping_size,0));
+  
+                                  // 2x2 cube
+  if (1)
+    {
+      tria=new Triangulation<3>();
+      tria_ptr.push_back(tria);
+      GridGenerator::hyper_cube(*tria, 1., 3.);
+      exact_areas.push_back(8.);
+    }
+
+                                  // arbitrary rectangle
+  if (1)
+    {
+      tria=new Triangulation<3>();
+      tria_ptr.push_back(tria);
+      GridGenerator::hyper_cube(*tria, 1., 3.);
+      Point<3> &v=tria->begin()->vertex(6);
+      v(0) = 5.;
+      v(1) = 4.;
+      v(2) = 4.5;
+      exact_areas.push_back(12.5);
+    }
+
+                                  // cube+part of ball
+  if (1)
+    {
+      Point<3> m(2,2,2);
+      Point<3> v(3,3,3);
+      double r=sqrt((m-v).square()),
+            h=r-1.5,
+           pi=acos(-1);
+      Boundary<3> *boundary1=new HyperBallBoundary<3>(m, r);
+      boundary_ptr.push_back(boundary1);
+      
+      tria=new Triangulation<3>();
+      tria_ptr.push_back(tria);
+      GridGenerator::hyper_cube(*tria, 1., 3.);
+      tria->set_boundary(1,*boundary1);
+      tria->begin_active()->face(3)->set_boundary_indicator(1);
+      exact_areas.push_back(8.+pi/3*h*h*(3*r-h));
+    }
+
+                                  // eighth of ball
+  if (1)
+    {
+      Point<3> p(0,0,0);
+      const double r=sqrt(3.);
+      Boundary<3> *boundary0=new HyperBallBoundary<3>(p, r);
+      boundary_ptr.push_back(boundary0);
+
+      tria=new Triangulation<3>();
+      tria_ptr.push_back(tria);
+      GridGenerator::hyper_cube(*tria, -1, 1.);
+      tria->set_boundary(0, *boundary0);
+      tria->refine_global(1);
+      const double pi=acos(-1);
+      exact_areas.push_back(4/3.*pi*r*r*r/8.);
+      for (unsigned int i=0; i<4; ++i)
+       show[3][i]=1;
+    }
+}
+
+  
+template<int dim>
+void mapping_test()
+{
+  deallog << "dim=" << dim << endl;
+  
+  vector<Mapping<dim> *> mapping_ptr;
+  vector<string> mapping_strings;
+
+  MappingQ1<dim> q1_old;
+  MappingQ<dim> q1(1);
+  MappingQ<dim> q2(2);
+  MappingQ<dim> q3(3);
+  MappingQ<dim> q4(4);
+  mapping_ptr.push_back(&q1_old);
+  mapping_ptr.push_back(&q1);
+  mapping_ptr.push_back(&q2);
+  mapping_ptr.push_back(&q3);
+  mapping_ptr.push_back(&q4);
+  mapping_strings.push_back("Q1fixed");
+  mapping_strings.push_back("Q1");
+  mapping_strings.push_back("Q2");
+  mapping_strings.push_back("Q3");
+  mapping_strings.push_back("Q4");
+
+  mapping_size=mapping_ptr.size();
+  
+  vector<Triangulation<dim> *> tria_ptr;
+  vector<Boundary<dim> *> boundary_ptr;
+  vector<double> exact_areas;
+  
+  create_triangulations(tria_ptr, boundary_ptr, exact_areas);
+  Assert(show.size()==tria_ptr.size(), ExcInternalError());
+
+  FE_Q<dim> fe_q4(4);
+  
+  for (unsigned int i=0; i<tria_ptr.size(); ++i)
+    {
+      DoFHandler<dim> dof(*tria_ptr[i]);
+      dof.distribute_dofs(fe_q4);      
+      DoFHandler<dim>::cell_iterator cell = dof.begin_active();
+      
+      deallog << "Triangulation" << i << ":" << endl;
+
+      deallog << "exact_area=" << exact_areas[i] << endl;
+      for (unsigned int j=0; j<mapping_size; ++j)
+       if (show[i][j])
+         {
+           char* st2 = new char[100];
+
+           if (true)
+             {
+               ostrstream ost(st2, 99);
+               ost << "Mapping" << dim << "d-" << i << '-'
+                   << mapping_strings[j] << ".dat" << ends;
+               deallog << st2 << endl;
+               plot_transformation(*mapping_ptr[j], fe_q4, cell, st2);
+               compute_area(*mapping_ptr[j], fe_q4, cell);
+             }
+           
+           if (dim>1)
+             {
+               ostrstream ost(st2, 99);
+               ost << "MappingFace" << dim << "d-" << i << '-'
+                   << mapping_strings[j] << ".dat" << ends;
+               deallog << st2 << endl;     
+               plot_faces(*mapping_ptr[j], fe_q4, cell, st2);
+             }
+
+           if (dim>1)
+             {
+               ostrstream ost(st2, 99);
+               ost << "MappingSubface" << dim << "d-" << i << '-'
+                   << mapping_strings[j] << ".dat" << ends;
+               deallog << st2 << endl;     
+               plot_subfaces(*mapping_ptr[j], fe_q4, cell, st2);
+             }
+
+           
+                                  // Test for transform_*_to_*_cell
+           if (dim==2 && true)
+             {
+               Mapping<dim> &mapping=*mapping_ptr[j];
+               Point<dim> p_unit(6/7.,4/7.);
+               Point<dim> p_real=mapping.transform_unit_to_real_cell(cell, p_unit);
+               Point<dim> p_re_unit=mapping.transform_real_to_unit_cell(cell, p_real);
+               deallog << "p_unit=" << p_unit << ",  p_real=" << p_real
+                       << ",  p_re_unit=" << p_re_unit << endl;
+             }
+           
+           delete[] st2;
+         }    
+    }
+
+
+                                  // delete all triangulations and
+                                  // boundary objects
+  for (unsigned int i=0; i<tria_ptr.size(); ++i)
+    if (tria_ptr[i]!=0)
+      delete tria_ptr[i];
+
+  for (unsigned int i=0; i<boundary_ptr.size(); ++i)
+    if (boundary_ptr[i]!=0)
+      delete boundary_ptr[i];
+}
+
+
+
+
+int main()
+{
+  ofstream logfile("mapping.dat");
+  logfile.precision(3);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+                                  // ----------------------- 
+                                  // Tests for dim=1
+                                  // -----------------------
+  mapping_test<1>();
+
+  
+                                  // ----------------------- 
+                                  // Tests for dim=2
+                                  // -----------------------
+  mapping_test<2>();
+  
+
+                                  // ----------------------- 
+                                  // Tests for dim=3
+                                  // -----------------------
+  mapping_test<3>();
+}
+
diff --git a/tests/fe/shapes.cc b/tests/fe/shapes.cc
new file mode 100644 (file)
index 0000000..2b79f09
--- /dev/null
@@ -0,0 +1,241 @@
+// $Id$
+// (c) Guido Kanschat
+//
+// Show the shape functions implemented.
+
+#include <base/quadrature_lib.h>
+#include <lac/vector.h>
+#include <grid/tria.h>
+#include <grid/tria_iterator.h>
+#include <dofs/dof_accessor.h>
+#include <grid/grid_generator.h>
+#include <fe/fe_q.h>
+#include <fe/fe_dgq.h>
+#include <fe/fe_system.h>
+#include <fe/mapping_q1.h>
+#include <fe/fe_values.h>
+#include <vector>
+#include <fstream>
+#include <string>
+
+char fname[50];
+
+template<int dim>
+inline void
+plot_shape_functions(Mapping<dim>& mapping,
+                    FiniteElement<dim>& finel,
+                    const char* name)
+{
+  Triangulation<dim> tr;
+  DoFHandler<dim> dof(tr);
+  GridGenerator::hyper_cube(tr, 0., 1.);
+  DoFHandler<dim>::cell_iterator c = dof.begin();
+  dof.distribute_dofs(finel);
+
+  const unsigned int div = 11;
+
+  QTrapez<1> q_trapez;
+  QIterated<dim> q(q_trapez, div);
+  FEValues<dim> fe(mapping, finel, q, UpdateFlags(update_values));
+  fe.reinit(c);
+  
+  sprintf(fname, "Shapes%dd-%s.dat", dim, name);
+  ofstream gnuplot(fname);
+  gnuplot.setf(ios::fixed);
+  gnuplot.precision (2);
+
+  unsigned int k=0;
+  for (unsigned int mz=0;mz<=((dim>2) ? div : 0) ;++mz)
+    {
+      for (unsigned int my=0;my<=((dim>1) ? div : 0) ;++my)
+       {
+         for (unsigned int mx=0;mx<=div;++mx)
+           {
+             gnuplot << q.point(k);
+         
+             for (unsigned int i=0;i<finel.dofs_per_cell;++i)
+               {
+                 gnuplot << " " << fe.shape_value(i,k) + 1.;
+               }
+             gnuplot << endl;
+             k++;
+           }
+         gnuplot << endl;
+       }
+      gnuplot << endl;
+    }
+}
+
+
+
+template<int dim>
+inline void
+plot_face_shape_functions(Mapping<dim>& mapping,
+                         FiniteElement<dim>& finel,
+                         const char* name)
+{
+  Triangulation<dim> tr;
+  DoFHandler<dim> dof(tr);
+  GridGenerator::hyper_cube(tr, 0., 1.);
+  tr.refine_global(1);
+  DoFHandler<dim>::active_cell_iterator c = dof.begin_active();
+  ++c;
+  c->set_refine_flag();
+  tr.execute_coarsening_and_refinement ();
+  c = dof.begin_active();
+
+  dof.distribute_dofs(finel);
+
+  const unsigned int div = 4;
+
+  QTrapez<1> q_trapez;
+  QIterated<dim-1> q(q_trapez, div);
+  FEFaceValues<dim> fe(mapping, finel, q, UpdateFlags(update_values
+    | update_q_points));
+  FESubfaceValues<dim> sub(mapping, finel, q, UpdateFlags(update_values
+    | update_q_points));
+
+  sprintf(fname, "ShapesFace%dd-%s.dat", dim, name);
+  ofstream gnuplot(fname);
+  gnuplot.setf(ios::fixed);
+  gnuplot.precision (2);
+  
+  for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+    {
+      if (!c->face(f)->has_children())
+       {
+         fe.reinit(c, f);
+  
+         unsigned int k=0;
+         for (unsigned int my=0;my<=((dim>2) ? div : 0) ;++my)
+           {
+             for (unsigned int mx=0;mx<=div;++mx)
+               {
+                 gnuplot << fe.quadrature_point(k);
+                 
+                 for (unsigned int i=0;i<finel.dofs_per_cell;++i)
+                   {
+                     gnuplot << " " << fe.shape_value(i,k) + 1.;
+                   }
+                 gnuplot << endl;
+                 k++;
+               }
+             gnuplot << endl;
+           }
+         gnuplot << endl;
+       } else {
+         for (unsigned int s=0;s<GeometryInfo<dim>::subfaces_per_face; ++s)
+           {
+             sub.reinit(c, f, s);
+             
+             unsigned int k=0;
+             for (unsigned int my=0;my<=((dim>2) ? div : 0) ;++my)
+               {
+                 for (unsigned int mx=0;mx<=div;++mx)
+                   {
+                     gnuplot << sub.quadrature_point(k);
+                     
+                     for (unsigned int i=0;i<finel.dofs_per_cell;++i)
+                       {
+                         gnuplot << " " << sub.shape_value(i,k) + 1.;
+                       }
+                     gnuplot << endl;
+                     k++;
+                   }
+                 gnuplot << endl;
+               }
+             gnuplot << endl;        
+           }
+       }
+    }
+}
+
+
+template<>
+void plot_face_shape_functions (Mapping<1>&,
+                               FiniteElement<1>&,
+                               const char*)
+{}
+
+
+template<int dim>
+void plot_FE_Q_shape_functions()
+{
+  MappingQ1<dim> m;
+  FE_Q<dim> q1(1);
+  plot_shape_functions(m, q1, "Q1");
+  plot_face_shape_functions(m, q1, "Q1");
+  FE_Q<dim> q2(2);
+  plot_shape_functions(m, q2, "Q2");
+  plot_face_shape_functions(m, q2, "Q2");
+  FE_Q<dim> q3(3);
+  plot_shape_functions(m, q3, "Q3");
+  plot_face_shape_functions(m, q3, "Q3");
+  FE_Q<dim> q4(4);
+  plot_shape_functions(m, q4, "Q4");
+  plot_face_shape_functions(m, q4, "Q4");
+//    FE_Q<dim> q5(5);
+//    plot_shape_functions(m, q5, "Q5");
+//    FE_Q<dim> q6(6);
+//    plot_shape_functions(m, q6, "Q6");
+//    FE_Q<dim> q7(7);
+//    plot_shape_functions(m, q7, "Q7");
+//    FE_Q<dim> q8(8);
+//    plot_shape_functions(m, q8, "Q8");
+//    FE_Q<dim> q9(9);
+//    plot_shape_functions(m, q9, "Q9");
+//    FE_Q<dim> q10(10);
+//    plot_shape_functions(m, q10, "Q10");
+}
+
+
+template<int dim>
+void plot_FE_DGQ_shape_functions()
+{
+  MappingQ1<dim> m;
+  FE_DGQ<dim> q1(1);
+  plot_shape_functions(m, q1, "DGQ1");
+  plot_face_shape_functions(m, q1, "DGQ1");
+  FE_DGQ<dim> q2(2);
+  plot_shape_functions(m, q2, "DGQ2");
+  plot_face_shape_functions(m, q2, "DGQ2");
+  FE_DGQ<dim> q3(3);
+  plot_shape_functions(m, q3, "DGQ3");
+  plot_face_shape_functions(m, q3, "DGQ3");
+  FE_DGQ<dim> q4(4);
+  plot_shape_functions(m, q4, "DGQ4");
+  plot_face_shape_functions(m, q4, "DGQ4");
+//    FE_DGQ<dim> q5(5);
+//    plot_shape_functions(m, q5, "DGQ5");
+//    FE_DGQ<dim> q6(6);
+//    plot_shape_functions(m, q6, "DGQ6");
+//    FE_DGQ<dim> q7(7);
+//    plot_shape_functions(m, q7, "DGQ7");
+//    FE_DGQ<dim> q8(8);
+//    plot_shape_functions(m, q8, "DGQ8");
+//    FE_DGQ<dim> q9(9);
+//    plot_shape_functions(m, q9, "DGQ9");
+//    FE_DGQ<dim> q10(10);
+//    plot_shape_functions(m, q10, "DGQ10");
+}
+
+
+int
+main()
+{
+  plot_FE_Q_shape_functions<1>();
+  plot_FE_Q_shape_functions<2>();
+  plot_FE_DGQ_shape_functions<2>();
+//  plot_FE_Q_shape_functions<3>();
+
+                                  // FESystem test.
+  MappingQ1<2> m;
+  FESystem<2> q2_q3(FE_Q<2>(2), 1,
+                   FE_Q<2>(3), 1);
+//  plot_shape_functions(m, q2_q3, "Q2_Q3");
+  
+  return 0;
+}
+
+
+

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.