From aa78afcef114772d0dd7edb5067a01e749e6fbf3 Mon Sep 17 00:00:00 2001 From: guido Date: Wed, 25 May 2005 05:36:35 +0000 Subject: [PATCH] still debugging git-svn-id: https://svn.dealii.org/trunk@10737 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/source/fe/fe_poly_tensor.cc | 23 +++-- .../source/fe/fe_raviart_thomas_nodal.cc | 11 ++- deal.II/deal.II/source/fe/fe_tools.cc | 82 +++++++++++++++++ tests/fe/Makefile | 4 +- tests/fe/rtn_1.cc | 92 +++++++++++++++++++ 5 files changed, 199 insertions(+), 13 deletions(-) create mode 100644 tests/fe/rtn_1.cc diff --git a/deal.II/deal.II/source/fe/fe_poly_tensor.cc b/deal.II/deal.II/source/fe/fe_poly_tensor.cc index b6fb9f3a8f..1d40c915f1 100644 --- a/deal.II/deal.II/source/fe/fe_poly_tensor.cc +++ b/deal.II/deal.II/source/fe/fe_poly_tensor.cc @@ -16,6 +16,7 @@ #include #include +#include template FE_PolyTensor::FE_PolyTensor ( @@ -64,7 +65,7 @@ FE_PolyTensor::shape_value_component ( return cached_values[i][component]; else for (unsigned int j=0;j::shape_grad_component ( return cached_grads[i][component]; else for (unsigned int j=0;j::get_data (const UpdateFlags update_flags, const UpdateFlags flags(data->update_flags); const unsigned int n_q_points = quadrature.n_quadrature_points; - + double d; + d = floor (1.3); + // some scratch arrays std::vector > values(0); std::vector > grads(0); @@ -228,7 +231,7 @@ FE_PolyTensor::get_data (const UpdateFlags update_flags, else for (unsigned int i=0; idofs_per_cell; ++i) for (unsigned int j=0; jdofs_per_cell; ++j) - data->shape_values[i][k] = inverse_node_matrix(i,j) * values[j]; + data->shape_values[i][k] = inverse_node_matrix(j,i) * values[j]; if (flags & update_gradients) if (inverse_node_matrix.n_cols() == 0) @@ -237,7 +240,7 @@ FE_PolyTensor::get_data (const UpdateFlags update_flags, else for (unsigned int i=0; idofs_per_cell; ++i) for (unsigned int j=0; jdofs_per_cell; ++j) - data->shape_grads[i][k] = inverse_node_matrix(i,j) * grads[j]; + data->shape_grads[i][k] = inverse_node_matrix(j,i) * grads[j]; } return data; } @@ -279,7 +282,7 @@ FE_PolyTensor::fill_fe_values ( for (unsigned int i=0; i::fill_fe_values ( if (flags & update_second_derivatives) this->compute_2nd (mapping, cell, dsd.cell(), mapping_data, fe_data, data); + deallog << "Shape data " << get_name() << std::endl; + + for (unsigned int i=0;i::FE_RaviartThomasNodal (const unsigned int deg) Assert (dim >= 2, ExcImpossibleInDim(dim)); this->mapping_type = this->independent_on_cartesian; - + // These must be done first, since + // they change the evaluation of + // basis functions + initialize_unit_support_points (deg); + initialize_node_matrix(); + for (unsigned int i=0; i::children_per_cell; ++i) this->prolongation[i].reinit (this->dofs_per_cell, this->dofs_per_cell); FETools::compute_embedding_matrices (*this, &this->prolongation[0]); - // finally fill in support points - // on cell and face - initialize_unit_support_points (deg); - initialize_node_matrix(); } diff --git a/deal.II/deal.II/source/fe/fe_tools.cc b/deal.II/deal.II/source/fe/fe_tools.cc index faed071f6d..76cb764746 100644 --- a/deal.II/deal.II/source/fe/fe_tools.cc +++ b/deal.II/deal.II/source/fe/fe_tools.cc @@ -602,6 +602,88 @@ FETools::compute_embedding_matrices(const FiniteElement& fe, } +template +void +FETools::compute_projection_matrices(const FiniteElement& fe, + FullMatrix* matrices) +{ + const unsigned int nc = GeometryInfo::children_per_cell; + const unsigned int n = fe.dofs_per_cell; + const unsigned int nd = fe.n_components(); + const unsigned int degree = fe.degree; + + for (unsigned int i=0;i tr_coarse; + Triangulation tr_fine; + GridGenerator::hyper_cube (tr_coarse, 0, 1); + GridGenerator::hyper_cube (tr_fine, 0, 1); + tr_fine.refine_global(1); + DoFHandler dof_coarse(tr_coarse); + dof_coarse.distribute_dofs(fe); + DoFHandler dof_fine(tr_fine); + dof_fine.distribute_dofs(fe); + + MappingCartesian mapping; + QGauss q_fine(degree+1); + const unsigned int nq = q_fine.n_quadrature_points; + + FEValues coarse (mapping, fe, q_fine, + update_q_points | update_JxW_values | update_values); + FEValues fine (mapping, fe, q_fine, + update_q_points | update_JxW_values | update_values); + + typename DoFHandler::active_cell_iterator coarse_cell + = dof_coarse.begin_active(); + typename DoFHandler::active_cell_iterator fine_cell; + + // Compute the coarse level mass + // matrix + coarse.reinit(dof_coarse); + FullMatrix A(n, n); + for (unsigned int k=0;k H(A); + + Vector v_coarse(n); + Vector v_fine(n); + + unsigned int cell_number = 0; + for (fine_cell = dof_fine.begin_active(); + fine_cell != dof_fine.end(); + ++fine_cell, ++cell_number) + { + // Compute right hand side, + // which is a fine level basis + // function tested with the + // coarse level functions. + Assert(false, ExcNotImplemented()); + + // Remove small entries from + // the matrix + for (unsigned int i=0; i void diff --git a/tests/fe/Makefile b/tests/fe/Makefile index 8a70014ed7..150a22d78c 100644 --- a/tests/fe/Makefile +++ b/tests/fe/Makefile @@ -1,7 +1,6 @@ -# Generated automatically from Makefile.in by configure. ############################################################ # Makefile,v 1.15 2002/06/13 12:51:13 hartmann Exp -# Copyright (C) 2000, 2001, 2002, 2003, 2004 by the deal.II authors +# Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005 by the deal.II authors ############################################################ ############################################################ @@ -49,6 +48,7 @@ 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) 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) diff --git a/tests/fe/rtn_1.cc b/tests/fe/rtn_1.cc new file mode 100644 index 0000000000..170ef7218e --- /dev/null +++ b/tests/fe/rtn_1.cc @@ -0,0 +1,92 @@ +//---------------------------------------------------------------------- +// rt_1.cc,v 1.3 2003/06/09 16:00:38 wolf Exp +// Version: +// +// 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. +// +//---------------------------------------------------------------------- + +// FE_RaviartThomasNodal +// Show the shape functions in support points + +#include "../tests.h" +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +extern "C" +{ +#include +}; + +template +void +check_support_points (const FiniteElement& fe) +{ + deallog << fe.get_name() << std::endl; + + const std::vector< Point > & points = fe.get_unit_support_points(); + std::vector weights(points.size()); + + Triangulation tr; + GridGenerator::hyper_cube(tr); + DoFHandler dof(tr); + dof.distribute_dofs(fe); + + MappingCartesian mapping; + + Quadrature support_quadrature(points, weights); + UpdateFlags flags = update_values; + FEValues vals(mapping, fe, support_quadrature, flags); + vals.reinit(dof.begin_active()); + + for (unsigned int k=0;k& p = points[k]; + deallog << p; + for (unsigned int i=0;i 1.e-12) +// deallog << "Error" << vals.shape_value_component(i,k,d) << std::endl; + } + } + deallog << std::endl; + } + deallog << std::endl; +} + + +int +main() +{ + std::ofstream logfile ("rtn_1.output"); + deallog.attach(logfile); + deallog.depth_console(10); + deallog.threshold_double(1.e-10); + + FE_RaviartThomasNodal<2> e20(0); + check_support_points(e20); + FE_RaviartThomasNodal<2> e21(1); + check_support_points(e21); +} -- 2.39.5