]> https://gitweb.dealii.org/ - dealii.git/commitdiff
new interpolation functions
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Sat, 24 Sep 2005 11:27:20 +0000 (11:27 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Sat, 24 Sep 2005 11:27:20 +0000 (11:27 +0000)
git-svn-id: https://svn.dealii.org/trunk@11527 0785d39b-7218-0410-832d-ea1e28bc413d

tests/fe/Makefile
tests/fe/interpolate_common.h [new file with mode: 0644]
tests/fe/interpolate_q1.cc [new file with mode: 0644]
tests/fe/interpolate_rtn.cc [new file with mode: 0644]
tests/fe/interpolate_system.cc [new file with mode: 0644]
tests/fe/rtdiff.cc
tests/fe/rtn_1.cc

index 02b1a13a1e25accb5f2138f69b090e6d0a870b8e..38eba44a9d32a82c09d2451c207cc60bf8db2c4a 100644 (file)
@@ -19,43 +19,11 @@ libraries = $(lib-deal2-1d.g) \
 
 default: run-tests
 
-############################################################
+%.exe: %.g.$(OBJEXT) $(libraries)
+       @echo =====linking======= $@
+       @$(CXX) $(LDFLAGS) -o $@ ../abort.o $^ $(LIBS)
 
-derivatives.exe         : derivatives.g.$(OBJEXT)         $(libraries)
-fe_data_test.exe        : fe_data_test.g.$(OBJEXT)        $(libraries)
-fe_tools.exe            : fe_tools.g.$(OBJEXT)            $(libraries)
-fe_tools_test.exe       : fe_tools_test.g.$(OBJEXT)       $(libraries)
-function.exe            : function.g.$(OBJEXT)            $(libraries)
-traits.exe              : traits.g.$(OBJEXT)              $(libraries)
-internals.exe           : internals.g.$(OBJEXT)           $(libraries)
-mapping.exe             : mapping.g.$(OBJEXT)             $(libraries)
-mapping_c1.exe          : mapping_c1.g.$(OBJEXT)          $(libraries)
-mapping_q1_eulerian.exe : mapping_q1_eulerian.g.$(OBJEXT) $(libraries)
-nedelec.exe             : nedelec.g.$(OBJEXT)             $(libraries)
-nedelec_2.exe           : nedelec_2.g.$(OBJEXT)           $(libraries)
-nedelec_3.exe           : nedelec_3.g.$(OBJEXT)           $(libraries)
-non_primitive_1.exe     : non_primitive_1.g.$(OBJEXT)     $(libraries)
-non_primitive_2.exe     : non_primitive_2.g.$(OBJEXT)     $(libraries)
-numbering.exe           : numbering.g.$(OBJEXT)           $(libraries)
-dgq_1.exe               : dgq_1.g.$(OBJEXT)               $(libraries)
-q_1.exe                 : q_1.g.$(OBJEXT)                 $(libraries)
-q_2.exe                 : q_2.g.$(OBJEXT)                 $(libraries)
-q_3.exe                 : q_3.g.$(OBJEXT)                 $(libraries)
-q_4.exe                 : q_4.g.$(OBJEXT)                 $(libraries)
-rt_1.exe                : rt_1.g.$(OBJEXT)                $(libraries)
-rt_2.exe                : rt_2.g.$(OBJEXT)                $(libraries)
-rt_3.exe                : rt_3.g.$(OBJEXT)                $(libraries)
-rt_4.exe                : rt_4.g.$(OBJEXT)                $(libraries)
-rt_5.exe                : rt_5.g.$(OBJEXT)                $(libraries)
-rt_6.exe                : rt_6.g.$(OBJEXT)                $(libraries)
-rtn_1.exe               : rtn_1.g.$(OBJEXT)               $(libraries)
-rtdiff.exe              : rtdiff.g.$(OBJEXT)              $(libraries)
-dgp_monomial_1.exe     : dgp_monomial_1.g.$(OBJEXT)      $(libraries)
-dgp_monomial_2.exe     : dgp_monomial_2.g.$(OBJEXT)      $(libraries)
-system_1.exe            : system_1.g.$(OBJEXT)            $(libraries)
-shapes.exe              : shapes.g.$(OBJEXT)              $(libraries)
-transfer.exe           : transfer.g.$(OBJEXT)            $(libraries)
-up_and_down.exe         : up_and_down.g.$(OBJEXT)         $(libraries)
+############################################################
 
 tests = fe_data_test traits fe_tools fe_tools_test mapping \
         mapping_c1 shapes derivatives function numbering mapping_q1_eulerian \
diff --git a/tests/fe/interpolate_common.h b/tests/fe/interpolate_common.h
new file mode 100644 (file)
index 0000000..6ba18b8
--- /dev/null
@@ -0,0 +1,137 @@
+//----------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2005 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.
+//
+//----------------------------------------------------------------------
+
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <lac/vector.h>
+#include <fe/fe.h>
+
+#include <iomanip>
+
+// Compute the maximal difference between local finite element
+// interpolate and given function.
+
+template <int dim>
+double difference(
+  const FiniteElement<dim>& fe,
+  const std::vector<double> dofs,
+  const Function<dim>& function)
+{
+  double result = 0.;
+  QGauss<dim> quadrature(fe.degree+1);
+  
+  std::vector<double> f(quadrature.n_quadrature_points);
+  function.value_list(quadrature.get_points(), f);
+  
+  for (unsigned int k=0;k<quadrature.n_quadrature_points;++k)
+    {
+      double diff = f[k];
+      for (unsigned int i=0;i<dofs.size();++i)
+       diff -= dofs[i] * fe.shape_value(i, quadrature.point(k));
+      diff = std::abs(diff);
+      result = std::max(result, diff);
+    }
+  return result;
+}
+
+template <int dim>
+double vector_difference(
+  const FiniteElement<dim>& fe,
+  const std::vector<double> dofs,
+  const Function<dim>& function,
+  const unsigned int offset)
+{
+  double result = 0.;
+  QGauss<dim> quadrature(fe.degree+1);
+  
+  std::vector<Vector<double> > f(quadrature.n_quadrature_points,
+                                Vector<double>(function.n_components));
+
+  function.vector_value_list(quadrature.get_points(), f);
+  
+  for (unsigned int k=0;k<quadrature.n_quadrature_points;++k)
+    for (unsigned int comp=0;comp<fe.n_components();++comp)
+      {
+       double diff = f[k](comp+offset);
+       for (unsigned int i=0;i<dofs.size();++i)
+         diff -= dofs[i] * fe.shape_value_component(i, quadrature.point(k),
+                                                    comp);
+       diff = std::abs(diff);
+       result = std::max(result, diff);
+      }
+  return result;
+}
+
+
+// Local implementation for any dimension
+
+
+template<int dim, int degree, int COMP=1>
+class
+Q1WedgeFunction : public Function<dim>
+{
+  public:
+    Q1WedgeFunction() : Function<dim> (COMP)
+      {}
+    
+    double value (const Point<dim>   &p,
+                 const unsigned int c) const
+      {
+       double result = 1.;
+       for (unsigned int d=0;d<dim;++d)
+         for (unsigned int k=0;k<degree;++k)
+           result *= p(d) + c;
+       return result;
+      }
+    
+    void value_list (const std::vector<Point<dim> > &points,
+                    std::vector<double>            &values,
+                    const unsigned int c) const
+      {
+       Assert (values.size() == points.size(),
+               ExcDimensionMismatch(values.size(), points.size()));
+       
+       for (unsigned int i=0;i<points.size();++i)
+         {
+           const Point<dim>& p = points[i];
+           double result = 1.;
+           for (unsigned int d=0;d<dim;++d)
+             for (unsigned int k=0;k<degree;++k)
+               result *= p(d) + c;
+           values[i] = result;
+         }
+      }
+
+    void vector_value_list (const std::vector<Point<dim> > &points,
+                           std::vector<Vector<double> >&   values) const
+      {
+       Assert (values.size() == points.size(),
+               ExcDimensionMismatch(values.size(), points.size()));
+       Assert (values[0].size() == n_components,
+               ExcDimensionMismatch(values.size(), n_components));
+       
+       for (unsigned int i=0;i<points.size();++i)
+         {
+           const Point<dim>& p = points[i];
+           for (unsigned int c=0;c<COMP;++c)
+             {
+               double result = 1.;
+               for (unsigned int d=0;d<dim;++d)
+                 for (unsigned int k=0;k<degree;++k)
+                   result *= p(d);
+               values[i](c) = result;
+             }
+         }
+      }
+};
+
diff --git a/tests/fe/interpolate_q1.cc b/tests/fe/interpolate_q1.cc
new file mode 100644 (file)
index 0000000..3987216
--- /dev/null
@@ -0,0 +1,132 @@
+//----------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2005 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.
+//
+//----------------------------------------------------------------------
+
+#include "interpolate_common.h"
+#include <base/logstream.h>
+#include <base/quadrature_lib.h>
+
+#include <fe/fe_q.h>
+#include <fe/fe_dgq.h>
+
+#include <fstream>
+
+// FE_Q<dim>::interpolate(...)
+// FE_DGQ<dim>::interpolate(...)
+
+template <int dim>
+void check(const Function<dim>& f,
+          const unsigned int degree)
+{
+  FE_Q<dim> fe(degree);
+  deallog << fe.get_name() << ' ';
+  
+  std::vector<double> dofs(fe.dofs_per_cell);
+  
+  std::vector<std::vector<double> >
+    values(1, std::vector<double>(fe.get_unit_support_points().size()));
+  f.value_list(fe.get_unit_support_points(), values[0]);
+  fe.interpolate(dofs, values[0]);
+  deallog << " value " << difference(fe,dofs,f);
+  fe.interpolate(dofs, values);
+  deallog << " vector " << difference(fe,dofs,f);
+
+  std::vector<Vector<double> >
+    vectors(fe.get_unit_support_points().size(), Vector<double>(1));
+  f.vector_value_list(fe.get_unit_support_points(), vectors);
+  fe.interpolate(dofs, vectors, 0);
+  deallog << " Vector " << difference(fe,dofs,f) << std::endl;
+}
+
+template <int dim>
+void check_dg(const Function<dim>& f,
+             const unsigned int degree)
+{
+  FE_DGQ<dim> fe(degree);
+  deallog << fe.get_name() << ' ';
+  
+  std::vector<double> dofs(fe.dofs_per_cell);
+  
+  std::vector<std::vector<double> >
+    values(1, std::vector<double>(fe.get_unit_support_points().size()));
+  f.value_list(fe.get_unit_support_points(), values[0]);
+  fe.interpolate(dofs, values[0]);
+  deallog << " value " << difference(fe,dofs,f);
+  fe.interpolate(dofs, values);
+  deallog << " vector " << difference(fe,dofs,f);
+
+  std::vector<Vector<double> >
+    vectors(fe.get_unit_support_points().size(), Vector<double>(1));
+  f.vector_value_list(fe.get_unit_support_points(), vectors);
+  fe.interpolate(dofs, vectors, 0);
+  deallog << " Vector " << difference(fe,dofs,f) << std::endl;
+}
+
+template <int dim>
+void check_dg_lobatto(const Function<dim>& f,
+                     const unsigned int degree)
+{
+  QGaussLobatto<1> fe_quadrature(degree);
+  FE_DGQ<dim> fe(fe_quadrature);
+  deallog << fe.get_name() << ' ';
+  
+  std::vector<double> dofs(fe.dofs_per_cell);
+  
+  std::vector<std::vector<double> >
+    values(1, std::vector<double>(fe.get_unit_support_points().size()));
+  f.value_list(fe.get_unit_support_points(), values[0]);
+  fe.interpolate(dofs, values[0]);
+  deallog << " value " << difference(fe,dofs,f);
+  fe.interpolate(dofs, values);
+  deallog << " vector " << difference(fe,dofs,f);
+
+  std::vector<Vector<double> >
+    vectors(fe.get_unit_support_points().size(), Vector<double>(1));
+  f.vector_value_list(fe.get_unit_support_points(), vectors);
+  fe.interpolate(dofs, vectors, 0);
+  deallog << " Vector " << difference(fe,dofs,f) << std::endl;
+}
+
+int main()
+{
+  std::ofstream logfile ("interpolate_q1.output");
+  deallog.attach(logfile);
+  deallog.depth_console(10);
+  deallog.threshold_double(2.e-15);
+
+  Q1WedgeFunction<1,1> w1;
+  check(w1,1);
+  check(w1,2);
+  check(w1,3);
+  check_dg(w1,1);
+  check_dg(w1,2);
+  check_dg(w1,3);
+  Q1WedgeFunction<2,1> w2;
+  check(w2,1);
+  check(w2,2);
+  check(w2,3);
+  check_dg(w2,2);
+  check_dg(w2,3);
+  Q1WedgeFunction<2,2> w22;
+  check(w22,2);
+  check(w22,3);
+  check_dg(w22,2);
+  check_dg(w22,3);
+  check_dg_lobatto(w22,4);
+  Q1WedgeFunction<2,3> w23;
+  check(w23,3);
+  Q1WedgeFunction<3,1> w3;
+  check_dg(w3,1);
+  check(w3,1);
+  check(w3,2);
+  check(w3,3);
+}
diff --git a/tests/fe/interpolate_rtn.cc b/tests/fe/interpolate_rtn.cc
new file mode 100644 (file)
index 0000000..d43e16c
--- /dev/null
@@ -0,0 +1,71 @@
+//----------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2005 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.
+//
+//----------------------------------------------------------------------
+
+#include "interpolate_common.h"
+#include <base/logstream.h>
+#include <base/quadrature_lib.h>
+
+#include <fe/fe_raviart_thomas.h>
+
+#include <fstream>
+
+// FE_Q<dim>::interpolate(...)
+
+template <int dim>
+void check1(const Function<dim>& f,
+           const unsigned int degree)
+{
+  FE_RaviartThomasNodal<dim> fe(degree);
+  deallog << fe.get_name() << ' ';
+  deallog << fe.get_generalized_support_points().size() << ' ';
+  
+  std::vector<double> dofs(fe.dofs_per_cell);
+  
+  std::vector<std::vector<double> >
+    values(dim, std::vector<double>(fe.get_generalized_support_points().size()));
+  std::vector<Vector<double> >
+    vectors(fe.get_generalized_support_points().size(),
+           Vector<double>(dim));
+  f.vector_value_list(fe.get_generalized_support_points(), vectors);
+
+  for (unsigned int c=0;c<values.size();++c)
+    for (unsigned int k=0;k<values[c].size();++k)
+      values[c][k] = vectors[k](c);
+  
+  fe.interpolate(dofs, values);
+  deallog << " vector " << vector_difference(fe,dofs,f,0);
+  
+  fe.interpolate(dofs, vectors, 0);
+  deallog << " Vector " << vector_difference(fe,dofs,f,0) << std::endl;
+}
+
+int main()
+{
+  std::ofstream logfile ("interpolate_rtn.output");
+  deallog.attach(logfile);
+  deallog.depth_console(10);
+  deallog.threshold_double(1.e-15);
+
+//   Q1WedgeFunction<1,1,2> w1;
+//   check1(w1,1,2);
+//   check1(w1,2,2);
+//   check1(w1,3,2);
+  Q1WedgeFunction<2,1,2> w2;
+  check1(w2,1);
+  check1(w2,2);
+  Q1WedgeFunction<2,2,2> w22;
+  check1(w22,2);
+//  Q1WedgeFunction<3,1,3> w3;
+//  check1(w3,1);
+//  check1(w3,2);
+}
diff --git a/tests/fe/interpolate_system.cc b/tests/fe/interpolate_system.cc
new file mode 100644 (file)
index 0000000..57e50bb
--- /dev/null
@@ -0,0 +1,107 @@
+//----------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2005 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.
+//
+//----------------------------------------------------------------------
+
+#include "interpolate_common.h"
+#include <base/logstream.h>
+#include <base/quadrature_lib.h>
+
+#include <fe/fe_q.h>
+#include <fe/fe_system.h>
+
+#include <fstream>
+
+// FE_Q<dim>::interpolate(...)
+
+template <int dim>
+void check1(const Function<dim>& f,
+           const unsigned int degree,
+           const unsigned int comp)
+{
+  FE_Q<dim> feq(degree);
+  FESystem<dim> fe(feq, comp);
+  deallog << fe.get_name() << ' ';
+  
+  std::vector<double> dofs(fe.dofs_per_cell);
+  
+  std::vector<std::vector<double> >
+    values(comp, std::vector<double>(fe.get_unit_support_points().size()));
+  std::vector<Vector<double> >
+    vectors(fe.get_unit_support_points().size(),
+           Vector<double>(f.n_components));
+  f.vector_value_list(fe.get_unit_support_points(), vectors);
+  for (unsigned int c=0;c<values.size();++c)
+    for (unsigned int k=0;k<values[c].size();++k)
+      values[c][k] = vectors[k](c);
+  
+  fe.interpolate(dofs, values);
+  deallog << " vector " << vector_difference(fe,dofs,f,0);
+  
+  fe.interpolate(dofs, vectors, 0);
+  deallog << " Vector " << vector_difference(fe,dofs,f,0) << std::endl;
+}
+
+template <int dim>
+void check3(const Function<dim>& f,
+           const unsigned int degree,
+           const unsigned int comp1,
+           const unsigned int comp2,
+           const unsigned int comp3)
+{
+  FE_Q<dim> feq1(degree);
+  FE_Q<dim> feq2(degree+1);
+  FE_Q<dim> feq3(degree+2);
+  FESystem<dim> fe(feq1, comp1, feq2, comp2, feq3, comp3);
+  deallog << fe.get_name() << ' ';
+  
+  std::vector<double> dofs(fe.dofs_per_cell);
+  
+  std::vector<std::vector<double> >
+    values(f.n_components,
+          std::vector<double>(fe.get_unit_support_points().size()));
+  std::vector<Vector<double> >
+    vectors(fe.get_unit_support_points().size(),
+           Vector<double>(f.n_components));
+  f.vector_value_list(fe.get_unit_support_points(), vectors);
+  for (unsigned int c=0;c<values.size();++c)
+    for (unsigned int k=0;k<values[c].size();++k)
+      values[c][k] = vectors[k](c);
+  
+  fe.interpolate(dofs, values);
+  deallog << " vector " << vector_difference(fe,dofs,f,0);
+  
+  fe.interpolate(dofs, vectors, 0);
+  deallog << " Vector " << vector_difference(fe,dofs,f,0) << std::endl;
+}
+
+int main()
+{
+  std::ofstream logfile ("interpolate_system.output");
+  deallog.attach(logfile);
+  deallog.depth_console(10);
+  deallog.threshold_double(1.e-15);
+
+  Q1WedgeFunction<1,1,2> w1;
+  check1(w1,1,2);
+  check1(w1,2,2);
+  check1(w1,3,2);
+  Q1WedgeFunction<2,2,3> w2;
+  check1(w2,2,3);
+  check1(w2,3,3);
+  Q1WedgeFunction<3,1,3> w3;
+  check1(w3,1,3);
+  check1(w3,2,3);
+  Q1WedgeFunction<2,1,9> www2;
+  check3(www2,1,2,3,4);
+  Q1WedgeFunction<3,1,9> www3;
+  check3(www3,1,2,3,4);
+}
index d8bd721ae15ff01ba967b90da6f62b1206f220e4..891738249318d4b5960244fe0e443963e35808be 100644 (file)
@@ -46,7 +46,7 @@ initialize_node_matrix (const FiniteElement<dim>& other,
   
   N.reinit(n_dofs, n_dofs);
 
-  const std::vector<Point<dim> >& unit_support_points = nodes.get_unit_support_points();
+  const std::vector<Point<dim> >& unit_support_points = nodes.get_generalized_support_points();
   
                                   // The curent node functional index
   unsigned int current = 0;
index 2e66fb06723ef03d7375bc0a24972223e2004f9b..d6e42e6fb017e28eb1d6ec28456c3bfff061ee4f 100644 (file)
@@ -35,7 +35,7 @@ check_support_points (const FiniteElement<dim>& fe)
 {
   deallog << fe.get_name() << std::endl;
   
-  const std::vector< Point<dim > > & points = fe.get_unit_support_points();
+  const std::vector< Point<dim > > & points = fe.get_generalized_support_points();
   std::vector<double> weights(points.size());
 
   Triangulation<dim> tr;
@@ -83,7 +83,7 @@ check_face_support_points (const FiniteElement<dim>& fe)
 {
   deallog << "Face " << fe.get_name() << std::endl;
   
-  const std::vector< Point<dim-1> > & sub_points = fe.get_unit_face_support_points();
+  const std::vector< Point<dim-1> > & sub_points = fe.get_generalized_face_support_points();
   std::vector<double> weights(sub_points.size());
   Quadrature<dim-1> sub_quadrature(sub_points, weights);
   std::vector< Point<dim> > points(sub_points.size());

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.