From 860aa89b10f93fbd5f512caaf00738ef28bc574c Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Sat, 24 Sep 2005 11:27:20 +0000 Subject: [PATCH] new interpolation functions git-svn-id: https://svn.dealii.org/trunk@11527 0785d39b-7218-0410-832d-ea1e28bc413d --- tests/fe/Makefile | 40 +--------- tests/fe/interpolate_common.h | 137 +++++++++++++++++++++++++++++++++ tests/fe/interpolate_q1.cc | 132 +++++++++++++++++++++++++++++++ tests/fe/interpolate_rtn.cc | 71 +++++++++++++++++ tests/fe/interpolate_system.cc | 107 +++++++++++++++++++++++++ tests/fe/rtdiff.cc | 2 +- tests/fe/rtn_1.cc | 4 +- 7 files changed, 454 insertions(+), 39 deletions(-) create mode 100644 tests/fe/interpolate_common.h create mode 100644 tests/fe/interpolate_q1.cc create mode 100644 tests/fe/interpolate_rtn.cc create mode 100644 tests/fe/interpolate_system.cc diff --git a/tests/fe/Makefile b/tests/fe/Makefile index 02b1a13a1e..38eba44a9d 100644 --- a/tests/fe/Makefile +++ b/tests/fe/Makefile @@ -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 index 0000000000..6ba18b829a --- /dev/null +++ b/tests/fe/interpolate_common.h @@ -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 +#include +#include +#include + +#include + +// Compute the maximal difference between local finite element +// interpolate and given function. + +template +double difference( + const FiniteElement& fe, + const std::vector dofs, + const Function& function) +{ + double result = 0.; + QGauss quadrature(fe.degree+1); + + std::vector f(quadrature.n_quadrature_points); + function.value_list(quadrature.get_points(), f); + + for (unsigned int k=0;k +double vector_difference( + const FiniteElement& fe, + const std::vector dofs, + const Function& function, + const unsigned int offset) +{ + double result = 0.; + QGauss quadrature(fe.degree+1); + + std::vector > f(quadrature.n_quadrature_points, + Vector(function.n_components)); + + function.vector_value_list(quadrature.get_points(), f); + + for (unsigned int k=0;k +class +Q1WedgeFunction : public Function +{ + public: + Q1WedgeFunction() : Function (COMP) + {} + + double value (const Point &p, + const unsigned int c) const + { + double result = 1.; + for (unsigned int d=0;d > &points, + std::vector &values, + const unsigned int c) const + { + Assert (values.size() == points.size(), + ExcDimensionMismatch(values.size(), points.size())); + + for (unsigned int i=0;i& p = points[i]; + double result = 1.; + for (unsigned int d=0;d > &points, + std::vector >& 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& p = points[i]; + for (unsigned int c=0;c +#include + +#include +#include + +#include + +// FE_Q::interpolate(...) +// FE_DGQ::interpolate(...) + +template +void check(const Function& f, + const unsigned int degree) +{ + FE_Q fe(degree); + deallog << fe.get_name() << ' '; + + std::vector dofs(fe.dofs_per_cell); + + std::vector > + values(1, std::vector(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 > + vectors(fe.get_unit_support_points().size(), Vector(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 +void check_dg(const Function& f, + const unsigned int degree) +{ + FE_DGQ fe(degree); + deallog << fe.get_name() << ' '; + + std::vector dofs(fe.dofs_per_cell); + + std::vector > + values(1, std::vector(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 > + vectors(fe.get_unit_support_points().size(), Vector(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 +void check_dg_lobatto(const Function& f, + const unsigned int degree) +{ + QGaussLobatto<1> fe_quadrature(degree); + FE_DGQ fe(fe_quadrature); + deallog << fe.get_name() << ' '; + + std::vector dofs(fe.dofs_per_cell); + + std::vector > + values(1, std::vector(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 > + vectors(fe.get_unit_support_points().size(), Vector(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 index 0000000000..d43e16ce08 --- /dev/null +++ b/tests/fe/interpolate_rtn.cc @@ -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 +#include + +#include + +#include + +// FE_Q::interpolate(...) + +template +void check1(const Function& f, + const unsigned int degree) +{ + FE_RaviartThomasNodal fe(degree); + deallog << fe.get_name() << ' '; + deallog << fe.get_generalized_support_points().size() << ' '; + + std::vector dofs(fe.dofs_per_cell); + + std::vector > + values(dim, std::vector(fe.get_generalized_support_points().size())); + std::vector > + vectors(fe.get_generalized_support_points().size(), + Vector(dim)); + f.vector_value_list(fe.get_generalized_support_points(), vectors); + + for (unsigned int c=0;c 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 index 0000000000..57e50bbae0 --- /dev/null +++ b/tests/fe/interpolate_system.cc @@ -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 +#include + +#include +#include + +#include + +// FE_Q::interpolate(...) + +template +void check1(const Function& f, + const unsigned int degree, + const unsigned int comp) +{ + FE_Q feq(degree); + FESystem fe(feq, comp); + deallog << fe.get_name() << ' '; + + std::vector dofs(fe.dofs_per_cell); + + std::vector > + values(comp, std::vector(fe.get_unit_support_points().size())); + std::vector > + vectors(fe.get_unit_support_points().size(), + Vector(f.n_components)); + f.vector_value_list(fe.get_unit_support_points(), vectors); + for (unsigned int c=0;c +void check3(const Function& f, + const unsigned int degree, + const unsigned int comp1, + const unsigned int comp2, + const unsigned int comp3) +{ + FE_Q feq1(degree); + FE_Q feq2(degree+1); + FE_Q feq3(degree+2); + FESystem fe(feq1, comp1, feq2, comp2, feq3, comp3); + deallog << fe.get_name() << ' '; + + std::vector dofs(fe.dofs_per_cell); + + std::vector > + values(f.n_components, + std::vector(fe.get_unit_support_points().size())); + std::vector > + vectors(fe.get_unit_support_points().size(), + Vector(f.n_components)); + f.vector_value_list(fe.get_unit_support_points(), vectors); + for (unsigned int c=0;c 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); +} diff --git a/tests/fe/rtdiff.cc b/tests/fe/rtdiff.cc index d8bd721ae1..8917382493 100644 --- a/tests/fe/rtdiff.cc +++ b/tests/fe/rtdiff.cc @@ -46,7 +46,7 @@ initialize_node_matrix (const FiniteElement& other, N.reinit(n_dofs, n_dofs); - const std::vector >& unit_support_points = nodes.get_unit_support_points(); + const std::vector >& unit_support_points = nodes.get_generalized_support_points(); // The curent node functional index unsigned int current = 0; diff --git a/tests/fe/rtn_1.cc b/tests/fe/rtn_1.cc index 2e66fb0672..d6e42e6fb0 100644 --- a/tests/fe/rtn_1.cc +++ b/tests/fe/rtn_1.cc @@ -35,7 +35,7 @@ check_support_points (const FiniteElement& fe) { deallog << fe.get_name() << std::endl; - const std::vector< Point > & points = fe.get_unit_support_points(); + const std::vector< Point > & points = fe.get_generalized_support_points(); std::vector weights(points.size()); Triangulation tr; @@ -83,7 +83,7 @@ check_face_support_points (const FiniteElement& fe) { deallog << "Face " << fe.get_name() << std::endl; - const std::vector< Point > & sub_points = fe.get_unit_face_support_points(); + const std::vector< Point > & sub_points = fe.get_generalized_face_support_points(); std::vector weights(sub_points.size()); Quadrature sub_quadrature(sub_points, weights); std::vector< Point > points(sub_points.size()); -- 2.39.5