//----------------------------------------------------------------------
// $Id$
//
-// Copyright (C) 2012 by the deal.II authors
+// Copyright (C) 2012, 2013 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
/**
* @file Test initialization of Assembler::MatrixSimple and
- * LocalResults
+ * DoFInfo
*/
#include "../tests.h"
#include <deal.II/base/logstream.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/block_indices.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
#include <deal.II/meshworker/local_results.h>
+#include <deal.II/meshworker/dof_info.h>
#include <deal.II/meshworker/simple.h>
+#include <deal.II/fe/fe_dgp.h>
+#include <deal.II/fe/fe_system.h>
+
using namespace dealii;
-template <class MATRIX>
-void test(MeshWorker::Assembler::MatrixSimple<MATRIX>& ass)
-{
- MeshWorker::LocalResults<double> info1;
- ass.initialize_info(info1, false);
+template <class DOFINFO, class MATRIX>
+void test(DOFINFO& info, MeshWorker::Assembler::MatrixSimple<MATRIX>& ass)
+{
+ ass.initialize_info(info, false);
deallog << "No faces" << std::endl;
- info1.print_debug(deallog);
+ info.print_debug(deallog);
- MeshWorker::LocalResults<double> info2;
- ass.initialize_info(info2, true);
+ ass.initialize_info(info, true);
deallog << "With faces" << std::endl;
- info2.print_debug(deallog);
+ info.print_debug(deallog);
}
int main()
deallog.attach(logfile);
deallog.depth_console (0);
+ Triangulation<2,2> tr;
+ GridGenerator::hyper_cube(tr);
+ FE_DGP<2,2> fe1(1);
+ FE_DGP<2,2> fe2(2);
+ FE_DGP<2,2> fe3(3);
+ FE_DGP<2,2> fe5(5);
+ FESystem<2,2> fes1(fe3,1,fe5,1,fe1,1);
+ FESystem<2,2> fes2(fe3,1,fe5,1,fe1,1, fe2,1);
+
+ DoFHandler<2,2> dof1(tr);
+ dof1.distribute_dofs(fes1);
+ DoFHandler<2,2> dof2(tr);
+ dof2.distribute_dofs(fes2);
+ dof1.initialize_local_block_info();
+ dof2.initialize_local_block_info();
+ MeshWorker::DoFInfo<2,2,double> info1(dof1);
+ MeshWorker::DoFInfo<2,2,double> info1b(dof1.block_info());
+ MeshWorker::DoFInfo<2,2,double> info2b(dof2.block_info());
+
MeshWorker::Assembler::MatrixSimple<FullMatrix<double> > ass1;
deallog.push("Single block");
- test(ass1);
- BlockIndices ind;
- ind.push_back(3);
- ind.push_back(5);
- ind.push_back(1);
- ass1.initialize_local_blocks(ind);
+ test(info1, ass1);
deallog.pop();
deallog.push("Multiple blocks");
- test(ass1);
- ind.push_back(2);
- deallog.pop();
- deallog.push("Same blocks");
- test(ass1);
- ass1.initialize_local_blocks(ind);
+ test(info1b, ass1);
deallog.pop();
deallog.push("More blocks");
- test(ass1);
+ test(info2b, ass1);
}
DEAL:Multiple blocks:: 2,0 0x0 face 2,0 0x0
DEAL:Multiple blocks:: 2,1 0x0 face 2,1 0x0
DEAL:Multiple blocks:: 2,2 0x0 face 2,2 0x0
-DEAL:Same blocks::No faces
-DEAL:Same blocks::J: 0
-DEAL:Same blocks::R: 0
-DEAL:Same blocks::M: 9 face 0
-DEAL:Same blocks:: 0,0 0x0
-DEAL:Same blocks:: 0,1 0x0
-DEAL:Same blocks:: 0,2 0x0
-DEAL:Same blocks:: 1,0 0x0
-DEAL:Same blocks:: 1,1 0x0
-DEAL:Same blocks:: 1,2 0x0
-DEAL:Same blocks:: 2,0 0x0
-DEAL:Same blocks:: 2,1 0x0
-DEAL:Same blocks:: 2,2 0x0
-DEAL:Same blocks::With faces
-DEAL:Same blocks::J: 0
-DEAL:Same blocks::R: 0
-DEAL:Same blocks::M: 9 face 9
-DEAL:Same blocks:: 0,0 0x0 face 0,0 0x0
-DEAL:Same blocks:: 0,1 0x0 face 0,1 0x0
-DEAL:Same blocks:: 0,2 0x0 face 0,2 0x0
-DEAL:Same blocks:: 1,0 0x0 face 1,0 0x0
-DEAL:Same blocks:: 1,1 0x0 face 1,1 0x0
-DEAL:Same blocks:: 1,2 0x0 face 1,2 0x0
-DEAL:Same blocks:: 2,0 0x0 face 2,0 0x0
-DEAL:Same blocks:: 2,1 0x0 face 2,1 0x0
-DEAL:Same blocks:: 2,2 0x0 face 2,2 0x0
DEAL:More blocks::No faces
DEAL:More blocks::J: 0
DEAL:More blocks::R: 0
//----------------------------------------------------------------------
// $Id$
//
-// Copyright (C) 2012 by the deal.II authors
+// Copyright (C) 2012, 2013 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
/**
* @file Test initialization of Assembler::MatrixSimple and
- * LocalResults
+ * DoFInfo
*/
#include "../tests.h"
#include <deal.II/base/logstream.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/block_indices.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
#include <deal.II/meshworker/local_results.h>
+#include <deal.II/meshworker/dof_info.h>
#include <deal.II/meshworker/simple.h>
+#include <deal.II/fe/fe_dgp.h>
+#include <deal.II/fe/fe_system.h>
+
using namespace dealii;
-template <class MATRIX>
-void test(MeshWorker::Assembler::MGMatrixSimple<MATRIX>& ass)
-{
- MeshWorker::LocalResults<double> info1;
- ass.initialize_info(info1, false);
+template <class DOFINFO, class MATRIX>
+void test(DOFINFO& info, MeshWorker::Assembler::MGMatrixSimple<MATRIX>& ass)
+{
+ ass.initialize_info(info, false);
deallog << "No faces" << std::endl;
- info1.print_debug(deallog);
+ info.print_debug(deallog);
- MeshWorker::LocalResults<double> info2;
- ass.initialize_info(info2, true);
+ ass.initialize_info(info, true);
deallog << "With faces" << std::endl;
- info2.print_debug(deallog);
+ info.print_debug(deallog);
}
int main()
deallog.attach(logfile);
deallog.depth_console (0);
+ Triangulation<2,2> tr;
+ GridGenerator::hyper_cube(tr);
+ FE_DGP<2,2> fe1(1);
+ FE_DGP<2,2> fe2(2);
+ FE_DGP<2,2> fe3(3);
+ FE_DGP<2,2> fe5(5);
+ FESystem<2,2> fes1(fe3,1,fe5,1,fe1,1);
+ FESystem<2,2> fes2(fe3,1,fe5,1,fe1,1, fe2,1);
+
+ DoFHandler<2,2> dof1(tr);
+ dof1.distribute_dofs(fes1);
+ DoFHandler<2,2> dof2(tr);
+ dof2.distribute_dofs(fes2);
+ dof1.initialize_local_block_info();
+ dof2.initialize_local_block_info();
+ MeshWorker::DoFInfo<2,2,double> info1(dof1);
+ MeshWorker::DoFInfo<2,2,double> info1b(dof1.block_info());
+ MeshWorker::DoFInfo<2,2,double> info2b(dof2.block_info());
+
MeshWorker::Assembler::MGMatrixSimple<FullMatrix<double> > ass1;
deallog.push("Single block");
- test(ass1);
- BlockIndices ind;
- ind.push_back(3);
- ind.push_back(5);
- ind.push_back(1);
- ass1.initialize_local_blocks(ind);
+ test(info1, ass1);
deallog.pop();
deallog.push("Multiple blocks");
- test(ass1);
- ind.push_back(2);
- deallog.pop();
- deallog.push("Same blocks");
- test(ass1);
- ass1.initialize_local_blocks(ind);
+ test(info1b, ass1);
deallog.pop();
deallog.push("More blocks");
- test(ass1);
+ test(info2b, ass1);
}
DEAL:Multiple blocks:: 2,0 0x0 face 2,0 0x0
DEAL:Multiple blocks:: 2,1 0x0 face 2,1 0x0
DEAL:Multiple blocks:: 2,2 0x0 face 2,2 0x0
-DEAL:Same blocks::No faces
-DEAL:Same blocks::J: 0
-DEAL:Same blocks::R: 0
-DEAL:Same blocks::M: 9 face 0
-DEAL:Same blocks:: 0,0 0x0
-DEAL:Same blocks:: 0,1 0x0
-DEAL:Same blocks:: 0,2 0x0
-DEAL:Same blocks:: 1,0 0x0
-DEAL:Same blocks:: 1,1 0x0
-DEAL:Same blocks:: 1,2 0x0
-DEAL:Same blocks:: 2,0 0x0
-DEAL:Same blocks:: 2,1 0x0
-DEAL:Same blocks:: 2,2 0x0
-DEAL:Same blocks::With faces
-DEAL:Same blocks::J: 0
-DEAL:Same blocks::R: 0
-DEAL:Same blocks::M: 9 face 9
-DEAL:Same blocks:: 0,0 0x0 face 0,0 0x0
-DEAL:Same blocks:: 0,1 0x0 face 0,1 0x0
-DEAL:Same blocks:: 0,2 0x0 face 0,2 0x0
-DEAL:Same blocks:: 1,0 0x0 face 1,0 0x0
-DEAL:Same blocks:: 1,1 0x0 face 1,1 0x0
-DEAL:Same blocks:: 1,2 0x0 face 1,2 0x0
-DEAL:Same blocks:: 2,0 0x0 face 2,0 0x0
-DEAL:Same blocks:: 2,1 0x0 face 2,1 0x0
-DEAL:Same blocks:: 2,2 0x0 face 2,2 0x0
DEAL:More blocks::No faces
DEAL:More blocks::J: 0
DEAL:More blocks::R: 0
--- /dev/null
+//--------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2012, 2013 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.
+//
+//--------------------------------------------------------------------
+
+// Test the functions in integrators/elasticity.h
+// Output matrices and assert consistency of residuals
+
+#include "../tests.h"
+#include "../lib/test_grids.h"
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/integrators/laplace.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_nedelec.h>
+#include <deal.II/fe/fe_raviart_thomas.h>
+
+using namespace LocalIntegrators::Laplace;
+
+template <int dim>
+void test_cell(const FEValuesBase<dim>& fev)
+{
+ const unsigned int n = fev.dofs_per_cell;
+ FullMatrix<double> M(n,n);
+ cell_matrix(M,fev);
+ M.print(deallog,8);
+
+ Vector<double> u(n), v(n), w(n);
+ std::vector<std::vector<Tensor<1,dim> > >
+ ugrad(dim,std::vector<Tensor<1,dim> >(fev.n_quadrature_points));
+
+ std::vector<unsigned int> indices(n);
+ for (unsigned int i=0;i<n;++i)
+ indices[i] = i;
+
+ deallog << "Residuals";
+ for (unsigned int i=0;i<n;++i)
+ {
+ u = 0.;
+ u(i) = 1.;
+ w = 0.;
+ fev.get_function_gradients(u, indices, ugrad, true);
+ cell_residual(w, fev, make_slice(ugrad));
+ M.vmult(v,u);
+ w.add(-1., v);
+ deallog << ' ' << w.l2_norm();
+ }
+ deallog << std::endl;
+}
+
+
+template <int dim>
+void test_boundary(const FEValuesBase<dim>& fev)
+{
+ const unsigned int n = fev.dofs_per_cell;
+ unsigned int d=fev.get_fe().n_components();
+ FullMatrix<double> M(n,n);
+ nitsche_matrix(M, fev, 17);
+ M.print(deallog,8);
+
+ Vector<double> u(n), v(n), w(n);
+ std::vector<std::vector<double> >
+ uval (d,std::vector<double>(fev.n_quadrature_points)),
+ null_val(d,std::vector<double>(fev.n_quadrature_points, 0.));
+ std::vector<std::vector<Tensor<1,dim> > >
+ ugrad (d,std::vector<Tensor<1,dim> >(fev.n_quadrature_points));
+
+ std::vector<unsigned int> indices(n);
+ for (unsigned int i=0;i<n;++i)
+ indices[i] = i;
+
+ deallog << "Residuals";
+ for (unsigned int i=0;i<n;++i)
+ {
+ u = 0.;
+ u(i) = 1.;
+ w = 0.;
+ fev.get_function_values(u, indices, uval, true);
+ fev.get_function_gradients(u, indices, ugrad, true);
+ nitsche_residual(w, fev, make_slice(uval), make_slice(ugrad), make_slice(null_val), 17);
+ M.vmult(v,u);
+ w.add(-1., v);
+ deallog << ' ' << w.l2_norm();
+ if (d==1)
+ {
+ nitsche_residual(w, fev, uval[0], ugrad[0], null_val[0], 17);
+ M.vmult(v,u);
+ w.add(-1., v);
+ deallog << " e" << w.l2_norm();
+ }
+ }
+ deallog << std::endl;
+}
+
+
+template <int dim>
+void test_face(const FEValuesBase<dim>& fev1,
+ const FEValuesBase<dim>& fev2)
+{
+ const unsigned int n1 = fev1.dofs_per_cell;
+ const unsigned int n2 = fev1.dofs_per_cell;
+ unsigned int d=fev1.get_fe().n_components();
+ FullMatrix<double> M11(n1,n1);
+ FullMatrix<double> M12(n1,n2);
+ FullMatrix<double> M21(n2,n1);
+ FullMatrix<double> M22(n2,n2);
+
+ ip_matrix(M11, M12, M21, M22, fev1, fev2, 17);
+ deallog << "M11" << std::endl;
+ M11.print(deallog,8);
+ deallog << "M22" << std::endl;
+ M22.print(deallog,8);
+ deallog << "M12" << std::endl;
+ M12.print(deallog,8);
+ deallog << "M21" << std::endl;
+ M21.print(deallog,8);
+
+ Vector<double> u1(n1), v1(n1), w1(n1);
+ Vector<double> u2(n2), v2(n2), w2(n2);
+ std::vector<std::vector<double> >
+ u1val (d,std::vector<double>(fev1.n_quadrature_points)),
+ u2val (d,std::vector<double>(fev2.n_quadrature_points));
+ std::vector<std::vector<Tensor<1,dim> > >
+ u1grad(d,std::vector<Tensor<1,dim> >(fev1.n_quadrature_points)),
+ u2grad(d,std::vector<Tensor<1,dim> >(fev2.n_quadrature_points));
+
+ std::vector<unsigned int> indices1(n1), indices2(n2);
+ for (unsigned int i=0;i<n1;++i) indices1[i] = i;
+ for (unsigned int i=0;i<n2;++i) indices2[i] = i;
+
+ deallog << "Residuals"; for (unsigned int i=0;i<n1;++i) indices1[i] = i;
+
+ for (unsigned int i1=0;i1<n1;++i1)
+ for (unsigned int i2=0;i2<n2;++i2)
+ {
+ u1 = 0.;
+ u1(i1) = 1.;
+ w1 = 0.;
+ fev1.get_function_values(u1, indices1, u1val, true);
+ fev1.get_function_gradients(u1, indices1, u1grad, true);
+ u2 = 0.;
+ u2(i2) = 1.;
+ w2 = 0.;
+ fev2.get_function_values(u2, indices2, u2val, true);
+ fev2.get_function_gradients(u2, indices2, u2grad, true);
+ ip_residual(w1, w2, fev1, fev2,
+ make_slice(u1val), make_slice(u1grad),
+ make_slice(u2val), make_slice(u2grad),
+ 17);
+ M11.vmult(v1,u1); w1.add(-1., v1);
+ M12.vmult(v1,u2); w1.add(-1., v1);
+ M21.vmult(v2,u1); w2.add(-1., v2);
+ M22.vmult(v2,u2); w2.add(-1., v2);
+ deallog << ' ' << w1.l2_norm() + w2.l2_norm();
+ if (d==1)
+ {
+ ip_residual(w1, w2, fev1, fev2, u1val[0], u1grad[0], u2val[0], u2grad[0], 17);
+ M11.vmult(v1,u1); w1.add(-1., v1);
+ M12.vmult(v1,u2); w1.add(-1., v1);
+ M21.vmult(v2,u1); w2.add(-1., v2);
+ M22.vmult(v2,u2); w2.add(-1., v2);
+ deallog << " e" << w1.l2_norm() + w2.l2_norm();
+ }
+ }
+ deallog << std::endl;
+}
+
+
+template <int dim>
+void
+test_fe(Triangulation<dim>& tr, FiniteElement<dim>& fe)
+{
+ deallog << fe.get_name() << std::endl << "cell matrix" << std::endl;
+ QGauss<dim> quadrature(fe.tensor_degree()+1);
+ FEValues<dim> fev(fe, quadrature, update_gradients);
+
+ typename Triangulation<dim>::cell_iterator cell1 = tr.begin(1);
+ fev.reinit(cell1);
+ test_cell(fev);
+
+ QGauss<dim-1> face_quadrature(fe.tensor_degree()+1);
+ FEFaceValues<dim> fef1(fe, face_quadrature, update_values | update_gradients | update_normal_vectors);
+ for (unsigned int i=0;i<GeometryInfo<dim>::faces_per_cell;++i)
+ {
+ deallog << "boundary_matrix " << i << std::endl;
+ fef1.reinit(cell1, i);
+ test_boundary(fef1);
+ }
+
+ FEFaceValues<dim> fef2(fe, face_quadrature, update_values | update_gradients);
+ typename Triangulation<dim>::cell_iterator cell2 = cell1->neighbor(1);
+
+ deallog << "face_matrix " << 0 << std::endl;
+ cell1 = tr.begin(1);
+ fef1.reinit(cell1, 1);
+ fef2.reinit(cell2, 0);
+ test_face(fef1, fef2);
+}
+
+
+template <int dim>
+void
+test(Triangulation<dim>& tr)
+{
+ FE_DGQ<dim> q1(1);
+ FESystem<dim> fe1(q1,dim);
+ test_fe(tr, fe1);
+
+ FE_DGQ<dim> q2(2);
+ FESystem<dim> fe2(q2,dim);
+ test_fe(tr, fe2);
+
+ FE_Nedelec<dim> n1(1);
+ test_fe(tr, n1);
+
+ FE_RaviartThomas<dim> rt1(1);
+ test_fe(tr, rt1);
+}
+
+
+int main()
+{
+ initlog(__FILE__);
+ deallog.threshold_double(1.e-10);
+
+ Triangulation<2> tr2;
+ TestGrids::hypercube(tr2, 1);
+ test(tr2);
+
+ Triangulation<3> tr3;
+ TestGrids::hypercube(tr3, 1);
+ test(tr3);
+}