From: Wolfgang Bangerth Date: Tue, 31 May 2005 14:26:15 +0000 (+0000) Subject: New test X-Git-Tag: v8.0.0~13768 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3b37071b97f2e88d4fab627dd9f34daf2a00da75;p=dealii.git New test git-svn-id: https://svn.dealii.org/trunk@10788 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/tests/bits/rt_2.cc b/tests/bits/rt_2.cc new file mode 100644 index 0000000000..68b8b5a9d5 --- /dev/null +++ b/tests/bits/rt_2.cc @@ -0,0 +1,164 @@ +//---------------------------- rt_1.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2005 by the deal.II authors and Oliver Kayser-Herold +// +// 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. +// +//---------------------------- rt_1.cc --------------------------- + +// Small test to analyse the equivalence from both sides of a face of +// the normal component for the Raviart-Thomas element +// +// Test adapted from a small program by Oliver Kayser-Herold. + +#include "../tests.h" +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + + +void evaluate_normal (DoFHandler<2> &dof_handler, + Vector &solution) +{ + // This quadrature rule determines + // the points, where the continuity + // will be tested. + QGauss<1> quad (6); + FEFaceValues<2> fe_v_face (dof_handler.get_fe (), quad, + UpdateFlags(update_values | + update_q_points | + update_gradients | + update_normal_vectors | + update_JxW_values)); + + FEFaceValues<2> fe_v_face_n (dof_handler.get_fe (), quad, + UpdateFlags(update_values | + update_q_points | + update_gradients | + update_normal_vectors | + update_JxW_values)); + + const unsigned int n_q_face = quad.n_quadrature_points; + const unsigned int n_components = dof_handler.get_fe().n_components(); + + // Cell iterators + DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + for (; cell!=endc; ++cell) + { + for (unsigned int f = 0; f < GeometryInfo<2>::faces_per_cell; ++f) + { + if (!cell->face(f)->at_boundary ()) + { + deallog << "cell=" << cell + << ", face=" << f + << std::endl; + + fe_v_face.reinit (cell, f); + + const unsigned int neighbor=cell->neighbor_of_neighbor(f); + fe_v_face_n.reinit (cell->neighbor (f), neighbor); + + // Get values from + // solution vector (For + // Trap.Rule) + std::vector > this_value + (n_q_face, Vector(n_components)); + fe_v_face.get_function_values (solution, this_value); + + // Same for neighbor cell + std::vector > this_value_n + (n_q_face, Vector(n_components)); + fe_v_face_n.get_function_values (solution, this_value_n); + + for (unsigned int q_point=0; q_point vn = fe_v_face.normal_vector (q_point); + double nx = vn(0); + double ny = vn(1); + + double u = this_value[q_point](0); + double v = this_value[q_point](1); + + double u_n = this_value_n[q_point](0); + double v_n = this_value_n[q_point](1); + + deallog << "quadrature point " << q_point + << " (u-u_n)*nx+(v-v_n)*ny=" + << (u-u_n) * nx + (v-v_n) * ny + << ", u*nx+v*ny=" + << u * nx + v * ny + << ", u_n*nx+v_n*ny=" + << u_n * nx + v_n * ny + << std::endl; + } + } + } + } +} + + +int main () +{ + std::ofstream logfile("rt_2.output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + Triangulation<2> tria_test; + Point<2> p1 (0,0), + p2 (1, 1); + std::vector sub_div; + + sub_div.push_back (1); + sub_div.push_back (4); + + GridGenerator::subdivided_hyper_rectangle (tria_test, sub_div, p1, p2); + tria_test.refine_global (2); + + FE_RaviartThomas<2> fe (1); + DoFHandler<2> dof_handler(tria_test); + dof_handler.distribute_dofs (fe); + + // Fill solution vector with random + // values between 0 and 1. + Vector solution(dof_handler.n_dofs ()); + for (unsigned int i = 0; i < dof_handler.n_dofs (); ++i) + solution(i) = (double) rand () / (double) RAND_MAX; + + // Now check if the function is + // continuous in normal direction + // on uniform mesh: + deallog << "Uniform mesh test" << std::endl; + evaluate_normal (dof_handler, solution); + + + // Then test same on distorted mesh + tria_test.distort_random (0.05); + deallog << "Distorted mesh test" << std::endl; + evaluate_normal (dof_handler, solution); +}