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 \
--- /dev/null
+//----------------------------------------------------------------------
+// $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;
+ }
+ }
+ }
+};
+
--- /dev/null
+//----------------------------------------------------------------------
+// $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);
+}
--- /dev/null
+//----------------------------------------------------------------------
+// $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);
+}
--- /dev/null
+//----------------------------------------------------------------------
+// $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);
+}
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;
{
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;
{
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());