#include <fe/fe_poly_tensor.h>
#include <fe/fe_values.h>
+#include <base/logstream.h>
template <class POLY, int dim>
FE_PolyTensor<POLY,dim>::FE_PolyTensor (
return cached_values[i][component];
else
for (unsigned int j=0;j<inverse_node_matrix.n_cols();++j)
- s += inverse_node_matrix(i,j) * cached_values[j][component];
+ s += inverse_node_matrix(j,i) * cached_values[j][component];
return s;
}
return cached_grads[i][component];
else
for (unsigned int j=0;j<inverse_node_matrix.n_cols();++j)
- s += inverse_node_matrix(i,j) * cached_grads[j][component];
+ s += inverse_node_matrix(j,i) * cached_grads[j][component];
return s;
}
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<Tensor<1,dim> > values(0);
std::vector<Tensor<2,dim> > grads(0);
else
for (unsigned int i=0; i<this->dofs_per_cell; ++i)
for (unsigned int j=0; j<this->dofs_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)
else
for (unsigned int i=0; i<this->dofs_per_cell; ++i)
for (unsigned int j=0; j<this->dofs_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;
}
for (unsigned int i=0; i<n_dofs; ++i)
{
- const unsigned int first = data.shape_function_to_row_table[i];
+ const unsigned int first = i*dim;//data.shape_function_to_row_table[i];
if (flags & update_values)
for (unsigned int k=0; k<n_quad; ++k)
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<data.shape_values.n_rows();++i)
+ {
+ for (unsigned int j=0;j<n_quad;++j)
+ deallog << ' ' << data.shape_values(i,j);
+ deallog << std::endl;
+ }
}
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<GeometryInfo<dim>::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();
}
}
+template<int dim, typename number>
+void
+FETools::compute_projection_matrices(const FiniteElement<dim>& fe,
+ FullMatrix<number>* matrices)
+{
+ const unsigned int nc = GeometryInfo<dim>::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<nc;++i)
+ {
+ Assert(matrices[i].n() == n, ExcDimensionMismatch(matrices[i].n(),n));
+ Assert(matrices[i].m() == n, ExcDimensionMismatch(matrices[i].m(),n));
+ }
+
+ /*
+ * Set up two meshes, one with a
+ * single reference cell and the
+ * other refined, together with
+ * DoFHandler and finite elements.
+ */
+ Triangulation<dim> tr_coarse;
+ Triangulation<dim> tr_fine;
+ GridGenerator::hyper_cube (tr_coarse, 0, 1);
+ GridGenerator::hyper_cube (tr_fine, 0, 1);
+ tr_fine.refine_global(1);
+ DoFHandler<dim> dof_coarse(tr_coarse);
+ dof_coarse.distribute_dofs(fe);
+ DoFHandler<dim> dof_fine(tr_fine);
+ dof_fine.distribute_dofs(fe);
+
+ MappingCartesian<dim> mapping;
+ QGauss<dim> q_fine(degree+1);
+ const unsigned int nq = q_fine.n_quadrature_points;
+
+ FEValues<dim> coarse (mapping, fe, q_fine,
+ update_q_points | update_JxW_values | update_values);
+ FEValues<dim> fine (mapping, fe, q_fine,
+ update_q_points | update_JxW_values | update_values);
+
+ typename DoFHandler<dim>::active_cell_iterator coarse_cell
+ = dof_coarse.begin_active();
+ typename DoFHandler<dim>::active_cell_iterator fine_cell;
+
+ // Compute the coarse level mass
+ // matrix
+ coarse.reinit(dof_coarse);
+ FullMatrix<number> A(n, n);
+ for (unsigned int k=0;k<nq;++k)
+ for (unsigned int i=0;i<n;++j)
+ for (unsigned int j=0;j<n;++j)
+ A(i,j) = coarse.JxW(k)
+ * coarse.shape_value(i,k)
+ * coarse.shape_value(j,k);
+
+ Householder<double> H(A);
+
+ Vector<number> v_coarse(n);
+ Vector<number> 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<matrix.m(); ++i)
+ for (unsigned int j=0; j<matrix.n(); ++j)
+ if (std::fabs(matrix(i,j)) < 1e-12)
+ matrix(i,j) = 0.;
+ }
+}
+
+
template <int dim,
class InVector, class OutVector>
void
-# 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
############################################################
############################################################
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)
--- /dev/null
+//----------------------------------------------------------------------
+// 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 <base/logstream.h>
+#include <grid/grid_generator.h>
+#include <fe/fe_raviart_thomas.h>
+#include <fe/mapping_cartesian.h>
+#include <fe/fe_values.h>
+
+#include <vector>
+#include <fstream>
+#include <string>
+
+#include <cmath>
+extern "C"
+{
+#include <math.h>
+};
+
+template <int dim>
+void
+check_support_points (const FiniteElement<dim>& fe)
+{
+ deallog << fe.get_name() << std::endl;
+
+ const std::vector< Point<dim > > & points = fe.get_unit_support_points();
+ std::vector<double> weights(points.size());
+
+ Triangulation<dim> tr;
+ GridGenerator::hyper_cube(tr);
+ DoFHandler<dim> dof(tr);
+ dof.distribute_dofs(fe);
+
+ MappingCartesian<dim> mapping;
+
+ Quadrature<dim> support_quadrature(points, weights);
+ UpdateFlags flags = update_values;
+ FEValues<dim> vals(mapping, fe, support_quadrature, flags);
+ vals.reinit(dof.begin_active());
+
+ for (unsigned int k=0;k<points.size();++k)
+ {
+ const Point<dim>& p = points[k];
+ deallog << p;
+ for (unsigned int i=0;i<fe.dofs_per_cell;++i)
+ {
+ deallog << '\t';
+ for (unsigned int d=0;d<dim;++d)
+ {
+ const double sf = fe.shape_value_component(i,p,d);
+ const double sv = vals.shape_value_component(i,k,d);
+ deallog << ' ' << (int) round(60*sf)
+ << '/' << (int) round(60*sv);
+
+// const double diff = std::abs(sf - vals.shape_value_component(i,k,d));
+// if (diff > 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);
+}