From 6eb2b919bd2576f38a807b155d5952c9a25b510d Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 1 Aug 2006 15:10:19 +0000 Subject: [PATCH] Add Oliver's test 5 git-svn-id: https://svn.dealii.org/trunk@13562 0785d39b-7218-0410-832d-ea1e28bc413d --- tests/fe/rt_approximation_01.cc | 529 ++++++++++++++++++++++++++++ tests/fe/rt_approximation_01/rt.gpl | 184 ++++++++++ 2 files changed, 713 insertions(+) create mode 100644 tests/fe/rt_approximation_01.cc create mode 100644 tests/fe/rt_approximation_01/rt.gpl diff --git a/tests/fe/rt_approximation_01.cc b/tests/fe/rt_approximation_01.cc new file mode 100644 index 0000000000..50bc96da32 --- /dev/null +++ b/tests/fe/rt_approximation_01.cc @@ -0,0 +1,529 @@ +//---------------------------- rt_approximation_01.cc --------------------------- +// rt_approximation_01.cc,v 1.3 2003/06/09 16:00:38 wolf Exp +// Version: +// +// Copyright (C) 2003, 2005, 2006 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. +// +//---------------------------- rt_approximation_01.cc --------------------------- + +/* + * Little snippet, which tests, to what degree the RT/ABF elements + * can approximate a polynomial exactly on deformed geometries. + */ + + + +#include "../tests.h" +#include + +#define PRECISION 2 + +#include + +std::ofstream logfile ("rt_approximation_01/output"); + +char buf[1000]; + +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + +template +class TestMap1 : public Function +{ + public: + TestMap1 (const unsigned int n_components) : + Function (n_components) + {}; + + virtual ~TestMap1 () {}; + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + void vector_value (const Point &p, + Vector &return_value) const; +}; + + + +template +double +TestMap1::value (const Point &p, + const unsigned int component) const +{ + // u = x^2, v = y^2, dudx = 2 x, dvdy = 2 y, div u = 2x + 2y + // I.e. \int div u = 2 (for unit square) + + if (component == 0) + return (p(0) * p(0)); + if (component == 1) + return (p(1) * p(1)); + + return (0); +} + + +template +void TestMap1::vector_value (const Point &p, + Vector &return_value) const +{ + Assert (return_value.size() == this->n_components, + ExcDimensionMismatch (return_value.size(), this->n_components)); + + // Just fill the vector with the appropriate components + for (unsigned int iCount = 0; iCount < this->n_components; iCount++) + return_value (iCount) = value (p, iCount); +} + + + +///----------------------------------------------------------------------- + +template +class TestDef1 : public Function +{ + private: + const double phi; + + public: + TestDef1 (const unsigned int n_components, const double ph) : + Function (n_components), + phi (ph) + {}; + + virtual ~TestDef1 () {}; + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + void vector_value (const Point &p, + Vector &return_value) const; +}; + + + +template +double +TestDef1::value (const Point &p, + const unsigned int component) const +{ + Point<2> center; + center(0) = 0.5; + center(1) = 0.5; + double rad = p.distance (center), + phi_p = atan2 (p(0) - center(0), p(1) - center(1)); + + if (component == 0) + return rad * (sin (phi + phi_p) - sin (phi_p)); + else + return rad * (cos (phi + phi_p) - cos (phi_p)); +} + + +template +void TestDef1::vector_value (const Point &p, + Vector &return_value) const +{ + Assert (return_value.size() == this->n_components, + ExcDimensionMismatch (return_value.size(), this->n_components)); + for (unsigned int iCount = 0; iCount < this->n_components; iCount++) + return_value (iCount) = value (p, iCount); +} + + +///----------------------------------------------------------------------- + + +template +class TestDef2 : public Function +{ + private: + const double scale; + + public: + TestDef2 (const unsigned int n_components, const double sc) : + Function (n_components), + scale (sc) + {}; + + virtual ~TestDef2 () {}; + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + void vector_value (const Point &p, + Vector &return_value) const; +}; + + +template +double +TestDef2::value (const Point &p, + const unsigned int component) const +{ + double x = p(0), + y = p(1); + + if (component == 0) + return scale * x; + else + return scale * y; +} + + +template +void TestDef2::vector_value (const Point &p, + Vector &return_value) const +{ + Assert (return_value.size() == this->n_components, + ExcDimensionMismatch (return_value.size(), this->n_components)); + for (unsigned int iCount = 0; iCount < this->n_components; iCount++) + return_value (iCount) = value (p, iCount); +} + + +///----------------------------------------------------------------------- +// testDef3 implements parallelograms ... + + +template +class TestDef3 : public Function +{ + private: + const double scale; + + public: + TestDef3 (const unsigned int n_components, const double sc) : + Function (n_components), + scale (sc) + {}; + + virtual ~TestDef3 () {}; + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + void vector_value (const Point &p, + Vector &return_value) const; +}; + + +template +double +TestDef3::value (const Point &p, + const unsigned int component) const +{ + double y = p(1); + + if (component == 0) + return scale * y; + else + return 0; +} + + +template +void TestDef3::vector_value (const Point &p, + Vector &return_value) const +{ + Assert (return_value.size() == this->n_components, + ExcDimensionMismatch (return_value.size(), this->n_components)); + for (unsigned int iCount = 0; iCount < this->n_components; iCount++) + return_value (iCount) = value (p, iCount); +} + +// Wrapper class for polynomials. + +template +class TestPoly : public Function +{ + private: + std::vector > polys; + + public: + TestPoly (unsigned int deg) : + Function (2) + { + std::vector coeff (deg, 0.0); + for (unsigned int p=0; p < 4; ++p) + { + for (unsigned int i=0; i < deg; ++i) + { + double c = ((double) rand ()) / ((double) RAND_MAX + 1); + coeff[i] = c; + } + polys.push_back (Polynomials::Polynomial (coeff)); + } + }; + + virtual ~TestPoly () {}; + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + void vector_value (const Point &p, + Vector &return_value) const; +}; + + +template +double +TestPoly::value (const Point &p, + const unsigned int component) const +{ + double x = p(0), + y = p(1); + + // Ugly hack, but should do the job ... + if (component == 0) + return polys[0].value (x) + polys[1].value (y); + else + return polys[2].value (x) + polys[3].value (y); +} + + +template +void TestPoly::vector_value (const Point &p, + Vector &return_value) const +{ + Assert (return_value.size() == this->n_components, + ExcDimensionMismatch (return_value.size(), this->n_components)); + for (unsigned int iCount = 0; iCount < this->n_components; iCount++) + return_value (iCount) = value (p, iCount); +} + + + +/* + * Check the value of the derivative field. + */ + +double TestProjection (Mapping<2> &mapping, + DoFHandler<2> *dof_handler) +{ + Vector solution; + solution.reinit (dof_handler->n_dofs ()); + + // Create and Project test function to H_div space and evaluate the error + // If the error is in the range of machine precision, this polynomial + // degree can obviously be represented by projected field. + + for (unsigned int deg = 1; deg < 4; ++deg) + { + TestPoly<2> pol (deg); + + // Project solution onto RT FE field + ConstraintMatrix hn_constraints; + hn_constraints.clear (); + DoFTools::make_hanging_node_constraints (*dof_handler, + hn_constraints); + hn_constraints.close (); + VectorTools::project (mapping, *dof_handler, hn_constraints, + QGauss6<2> (), pol, + solution); + + // Now evaluate error ... + // Use a high order quadrature. + QGauss<2> quad (6); + FEValues<2> fe_values (mapping, dof_handler->get_fe (), quad, + UpdateFlags(update_values | + update_q_points | + update_gradients | + update_JxW_values| + update_contravariant_transformation)); + + const unsigned int n_q_points = 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(); + + double err_u = 0, + err_v = 0; + + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + + // Get values from solution vector (For Trap.Rule) + std::vector > this_value + (n_q_points, Vector(n_components)); + fe_values.get_function_values (solution, this_value); + + for (unsigned int q_point=0; q_point p = fe_values.quadrature_point (q_point); + + double u_ex = pol.value (p, 0), + v_ex = pol.value (p, 1); + + double JxW = fe_values.JxW (q_point); + + err_u += (u - u_ex) * (u - u_ex) * JxW; + err_v += (v - v_ex) * (v - v_ex) * JxW; + } + } + + sprintf (buf, "Deg %i ErrU %e ErrV %e\n", deg, err_u, err_v); + deallog << buf; + } + + + // Test the core functionality + DataOut<2> *data_out = new DataOut<2>; + data_out->attach_dof_handler (*dof_handler); + data_out->add_data_vector (solution, "solution"); + data_out->build_patches (mapping, 4); + + const char *szFName = "solu.gpl"; + data_out->write_gnuplot (deallog.get_file_stream()); + + delete data_out; + + return (0.0); +} + + +int main () +{ + logfile.precision (PRECISION); + logfile.setf(std::ios::fixed); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + Triangulation<2> tria_test; + DoFHandler<2> *dof_handler, + *dof_handler_def; + Point<2> p1 (0,0), + p2 (1, 1); + + GridGenerator::hyper_rectangle (tria_test, p1, p2); + // tria_test.refine_global (1); + // tria_test.distort_random (0.4); + + // Create a DoFHandler for the RT space + FE_RaviartThomas<2> fe (1); + //FE_ABF<2> fe (0); + dof_handler = new DoFHandler<2> (tria_test); + dof_handler->distribute_dofs (fe); + + QGauss6<2> quad_temp; + sprintf (buf, "DoFs per Quad: %i per Line %i per Vert %i\n", fe.dofs_per_quad, fe.dofs_per_line, fe.dofs_per_vertex); + deallog << buf; + sprintf (buf, "QPoints %i\n", quad_temp.n_quadrature_points); + deallog << buf; + + // Create an deformation object for the Eulerian mapping + FESystem<2> fe_def (FE_Q<2>(1), 2); + dof_handler_def = new DoFHandler<2> (tria_test); + dof_handler_def->distribute_dofs (fe_def); + + // Alloc some DoFs + Vector solution, + solution_q, + deformation; + solution.reinit (dof_handler->n_dofs ()); + solution_q.reinit (dof_handler->n_dofs ()); + deformation.reinit (dof_handler_def->n_dofs ()); + + ConstraintMatrix hn_constraints_def; + hn_constraints_def.clear (); + DoFTools::make_hanging_node_constraints (*dof_handler_def, + hn_constraints_def); + hn_constraints_def.close (); + + { + MappingQ1Eulerian<2> mapping_euler (deformation, *dof_handler_def); + + // Try rotating the elements + for (double rotat = 0; rotat < 2 * M_PI; rotat += 0.25 * M_PI) + { + // Rotate element + VectorTools::project (*dof_handler_def, hn_constraints_def, + QGauss6<2> (), TestDef1<2>(2, rotat), + deformation); + sprintf (buf, "phi = %e\n", rotat); + deallog << buf; + TestProjection (mapping_euler, dof_handler); + } + + // Try resizing the elements + for (double scale = -0.75; scale < 4.0; scale += 0.25) + { + VectorTools::project (*dof_handler_def, hn_constraints_def, + QGauss6<2> (), TestDef2<2>(2, scale), + deformation); + sprintf (buf, "Scale = %e\n", scale); + deallog << buf; + TestProjection (mapping_euler, dof_handler); + } + + // Try paralellogramming the elements + for (double scale = -1.0; scale < 1.0; scale += 0.25) + { + VectorTools::project (*dof_handler_def, hn_constraints_def, + QGauss6<2> (), TestDef3<2>(2, scale), + deformation); + sprintf (buf, "Scale = %e\n", scale); + deallog << buf; + TestProjection (mapping_euler, dof_handler); + } + + // Try arbitrary deformation ... + for (unsigned int i=0; i < deformation.size (); ++i) + { + double c = ((double) rand ()) / ((double) RAND_MAX + 1); + deformation (i) = 0.35 * (c - 0.5); + } + + sprintf (buf, "Arbitrary\n"); + deallog << buf; + TestProjection (mapping_euler, dof_handler); + } + + delete (dof_handler); + delete (dof_handler_def); + + return (0); +} diff --git a/tests/fe/rt_approximation_01/rt.gpl b/tests/fe/rt_approximation_01/rt.gpl new file mode 100644 index 0000000000..6a47c2f34e --- /dev/null +++ b/tests/fe/rt_approximation_01/rt.gpl @@ -0,0 +1,184 @@ +#!/usr/bin/gnuplot -persist +# +# +# G N U P L O T +# Version 4.0 patchlevel 0 +# last modified Thu Apr 15 14:44:22 CEST 2004 +# System: Linux 2.6.11.4-21.11-smp +# +# Copyright (C) 1986 - 1993, 1998, 2004, 2006 +# Thomas Williams, Colin Kelley and many others +# +# This is gnuplot version 4.0. Please refer to the documentation +# for command syntax changes. The old syntax will be accepted +# throughout the 4.0 series, but all save files use the new syntax. +# +# Type `help` to access the on-line reference manual. +# The gnuplot FAQ is available from +# http://www.gnuplot.info/faq/ +# +# Send comments and requests for help to +# +# Send bugs, suggestions and mods to +# +# +set terminal postscript eps color solid enhanced 20 +# set output +unset clip points +set clip one +unset clip two +set bar 1.000000 +set border 31 lt -1 lw 1.000 +set xdata +set ydata +set zdata +set x2data +set y2data +set timefmt x "%d/%m/%y,%H:%M" +set timefmt y "%d/%m/%y,%H:%M" +set timefmt z "%d/%m/%y,%H:%M" +set timefmt x2 "%d/%m/%y,%H:%M" +set timefmt y2 "%d/%m/%y,%H:%M" +set timefmt cb "%d/%m/%y,%H:%M" +set boxwidth +set style fill empty border +set dummy x,y +set format x "% g" +set format y "% g" +set format x2 "% g" +set format y2 "% g" +set format z "% g" +set format cb "% g" +set angles radians +unset grid +set key title "" +set key right top Right noreverse enhanced box linetype -2 linewidth 1.000 samplen 4 spacing 1 width 0 height 0 autotitles +unset label +unset arrow +unset style line +unset style arrow +unset logscale +set offsets 0, 0, 0, 0 +set pointsize 1 +set encoding default +unset polar +unset parametric +unset decimalsign +set view 60, 20, 1, 1 +set samples 100, 100 +set isosamples 10, 10 +set surface +unset contour +set clabel '%8.3g' +set mapping cartesian +set datafile separator whitespace +set hidden3d offset 1 trianglepattern 3 undefined 1 altdiagonal bentover +set cntrparam order 4 +set cntrparam linear +set cntrparam levels auto 5 +set cntrparam points 5 +set size ratio 0 1,1 +set origin 0,0 +set style data lines +set style function lines +set xzeroaxis lt -2 lw 1.000 +set yzeroaxis lt -2 lw 1.000 +set x2zeroaxis lt -2 lw 1.000 +set y2zeroaxis lt -2 lw 1.000 +set tics in +set ticslevel 0.5 +set ticscale 1 0.5 +set mxtics default +set mytics default +set mztics default +set mx2tics default +set my2tics default +set mcbtics default +set xtics border mirror norotate autofreq +set ytics border mirror norotate autofreq +set ztics border nomirror norotate autofreq +set nox2tics +set noy2tics +set cbtics border mirror norotate autofreq +set title "" 0.000000,0.000000 font "" +set timestamp "" bottom norotate 0.000000,0.000000 "" +set rrange [ * : * ] noreverse nowriteback # (currently [0.00000:10.0000] ) +set trange [ * : * ] noreverse nowriteback # (currently [-5.00000:5.00000] ) +set urange [ * : * ] noreverse nowriteback # (currently [-5.00000:5.00000] ) +set vrange [ * : * ] noreverse nowriteback # (currently [-5.00000:5.00000] ) +set xlabel "x" 0.000000,0.000000 font "" +set x2label "" 0.000000,0.000000 font "" +set xrange [ * : * ] noreverse nowriteback # (currently [-10.0000:10.0000] ) +set x2range [ * : * ] noreverse nowriteback # (currently [-10.0000:10.0000] ) +set ylabel "y" 0.000000,0.000000 font "" +set y2label "" 0.000000,0.000000 font "" +set yrange [ * : * ] noreverse nowriteback # (currently [-10.0000:10.0000] ) +set y2range [ * : * ] noreverse nowriteback # (currently [-10.0000:10.0000] ) +set zlabel "" 0.000000,0.000000 font "" +set zrange [ * : * ] noreverse nowriteback # (currently [-10.0000:10.0000] ) +set cblabel "" 0.000000,0.000000 font "" +set cbrange [ * : * ] noreverse nowriteback # (currently [-10.0000:10.0000] ) +set zero 1e-08 +set lmargin -1 +set bmargin -1 +set rmargin -1 +set tmargin -1 +set locale "C" +set pm3d scansautomatic flush begin noftriangles nohidden3d implicit corners2color mean +unset pm3d +set palette positive nops_allcF maxcolors 0 gamma 1.5 color model RGB +set palette rgbformulae 7, 5, 15 +set colorbox default +set colorbox vertical origin 0.9,0.2 size 0.1,0.63 bdefault +set loadpath +set fontpath +set fit noerrorvariables +set output "RT1_00_x.eps" +splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:3 title "v0_x" +set output "RT1_00_y.eps" +splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:4 title "v0_y" +set output "RT1_01_x.eps" +splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:5 title "v1_x" +set output "RT1_01_y.eps" +splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:6 title "v1_y" +set output "RT1_02_x.eps" +splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:7 title "v2_x" +set output "RT1_02_y.eps" +splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:8 title "v2_y" +set output "RT1_03_x.eps" +splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:9 title "v3_x" +set output "RT1_03_y.eps" +splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:10 title "v3_y" +set output "RT1_04_x.eps" +splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:11 title "v4_x" +set output "RT1_04_y.eps" +splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:12 title "v4_y" +set output "RT1_05_x.eps" +splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:13 title "v5_x" +set output "RT1_05_y.eps" +splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:14 title "v5_y" +set output "RT1_06_x.eps" +splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:15 title "v6_x" +set output "RT1_06_y.eps" +splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:16 title "v6_y" +set output "RT1_07_x.eps" +splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:17 title "v7_x" +set output "RT1_07_y.eps" +splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:18 title "v7_y" +set output "RT1_08_x.eps" +splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:19 title "v8_x" +set output "RT1_08_y.eps" +splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:20 title "v8_y" +set output "RT1_09_x.eps" +splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:21 title "v9_x" +set output "RT1_09_y.eps" +splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:22 title "v9_y" +set output "RT1_10_x.eps" +splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:23 title "v10_x" +set output "RT1_10_y.eps" +splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:24 title "v10_y" +set output "RT1_11_x.eps" +splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:25 title "v11_x" +set output "RT1_11_y.eps" +splot "< perl -ne \'print if s/DEAL:RT<2>.1.::value//;\' output" using 1:2:26 title "v11_y" +# EOF -- 2.39.5