From: Wolfgang Bangerth Date: Sun, 6 Sep 2009 15:02:12 +0000 (+0000) Subject: New test. X-Git-Tag: v8.0.0~7175 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=51c040799c55e90d05def81e0a52c037cf9cedab;p=dealii.git New test. git-svn-id: https://svn.dealii.org/trunk@19398 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/tests/deal.II/Makefile b/tests/deal.II/Makefile index 091e969279..52506660a3 100644 --- a/tests/deal.II/Makefile +++ b/tests/deal.II/Makefile @@ -84,7 +84,8 @@ tests_x = block_matrices \ coarsening_* \ number_cache_* \ kelly_crash_* \ - recursively_set_material_id + recursively_set_material_id \ + point_value_history_* # from above list of regular expressions, generate the real set of # tests diff --git a/tests/deal.II/point_value_history_01.cc b/tests/deal.II/point_value_history_01.cc new file mode 100644 index 0000000000..d6003afb75 --- /dev/null +++ b/tests/deal.II/point_value_history_01.cc @@ -0,0 +1,260 @@ +//---------------------------- point_value_history_01.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2009 by the deal.II authors and Michael Rapson +// +// 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. +// +//---------------------------- point_value_history_01.cc --------------------------- + + + +// ************************************************ +// A test program for the PointValueHistory class +// Currently this only tests a finite element system +// with 3 components, on a hyper cube and with Vectors. +// The testing is done in two dimensions. +// ************************************************ + + + + +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include + +#include + +#include + +using namespace dealii; + + + +template +class TestPointValueHistory +{ +public: + TestPointValueHistory(); + void run(); + +private: + Triangulation triangulation; + FESystem finite_element; + DoFHandler dof_handler; + PointValueHistory test_copy; +}; + + + + +template +TestPointValueHistory::TestPointValueHistory() : +finite_element(FE_Q (1 + 1), dim, FE_Q (1), 1), + +dof_handler(triangulation) +{ } + + + + +template +void TestPointValueHistory::run() +{ + // Make a triangulation + GridGenerator::hyper_cube(triangulation, 0, 1); + triangulation.refine_global(2); // refine 2 times to make 5 nodes per side + + // make a DOF handler, a model solution filled with ones and a flow vector + FESystem finite_element(FE_Q (1 + 1), dim, FE_Q (1), 1); + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(finite_element); + DoFRenumbering::Cuthill_McKee(dof_handler); + + // renumber for components so that same dof indices are used for BlockVectors and normal Vectors + std::vector block_component(dim + 1, 0); + block_component[dim] = 1; // component dim = pressure component! + DoFRenumbering::component_wise(dof_handler, block_component); + + // Vector + Vector solution, post_processed, poles; + solution.reinit(dof_handler.n_dofs()); + post_processed.reinit(dof_handler.n_dofs()); + poles.reinit(dof_handler.n_dofs()); + + // set up a simple linear discrete time system so that time plots vary + // over the mesh but can be easily checked. The basic idea is to have each + // component of the fe_system to depend on a specific dimension (i.e component 0 + // depends on dim 0, etc. % dim handles the case where there are more components + // than dimensions. The code breaks down at the edges of the mesh and this is + // not corrected for. + { + QGauss quadrature_formula(2); + FEValues fe_values(finite_element, quadrature_formula, update_values | update_quadrature_points); // just need local_dof_indices and quadrature_points + + std::vector local_dof_indices(finite_element.dofs_per_cell); + std::vector > dof_locations(finite_element.dofs_per_cell); + Vector cell_pole(finite_element.dofs_per_cell); + + typename DoFHandler::active_cell_iterator cell, endc; + cell = dof_handler.begin_active(); + endc = dof_handler.end(); + for (; cell != endc; ++cell) + { + fe_values.reinit(cell); // need to get local_dof_indices + cell->get_dof_indices(local_dof_indices); + dof_locations = fe_values.get_quadrature_points(); + cell_pole = 0; + for (unsigned int dof = 0; dof != finite_element.dofs_per_cell; dof++) + { + unsigned int dof_component = finite_element.system_to_component_index(dof).first; + + for (unsigned int q_point = 0; q_point < quadrature_formula.size(); ++q_point) + { + cell_pole(dof) += (fe_values.shape_value(dof, q_point) * dof_locations [q_point] (dof_component % dim)); + } + solution(local_dof_indices [dof]) = 1; // start all solutions at 1 + poles(local_dof_indices[dof]) -= cell_pole(dof); + + if (dof_component == dim) // components start numbering at 0 + poles(local_dof_indices[dof]) = -0.1; // dim+1th component is not handled well by the code above + } + + } // loop over all cells + poles.add(1.0); // slow down the pole settling time + post_processed.add(3.0); // set to starting value. + } + + // Setup monitor node to print variation over time + unsigned int n_inputs = 1; + PointValueHistory node_monitor (dof_handler, n_inputs); + PointValueHistory no_dof_handler (n_inputs); + + // check that the assignment operator is valid + test_copy = node_monitor; + test_copy.add_point(Point<2>::Point(1, 0.2)); + test_copy.add_field_name("Solution"); + std::vector < std::vector > > selected_locations; + test_copy.get_points(selected_locations); + test_copy.mark_locations(); + test_copy.close(); + test_copy.start_new_dataset(0.1); + test_copy.evaluate_field("Solution", solution); + std::vector input_value(n_inputs, 1); + test_copy.push_back_independent(input_value); + test_copy.write_gnuplot("Test_Copy"); + test_copy.status(deallog.get_file_stream()); + test_copy.clear (); + // end of assignment operator check + + { + node_monitor.add_field_name("Solution"); + node_monitor.add_field_name("Post Processed Vector"); // not sensitive to spaces + + // two alternatives here, adding a point at a time or a vector of points + // 2d points + std::vector > point_vector(5, Point < 2 > ::Point()); + point_vector[0] = Point < 2 > ::Point(0, 0); // some of these points will hit a node, others won't + point_vector[1] = Point < 2 > ::Point(0.25, 0); + point_vector[2] = Point < 2 > ::Point(0.25, 0.45); + point_vector[3] = Point < 2 > ::Point(0.45, 0.45); + point_vector[4] = Point < 2 > ::Point(0.8, 0.8); + + node_monitor.add_points(point_vector); + node_monitor.add_point(Point<2>::Point(1, 0.2)); // add a single point + + // MonitorNode requires that the instance is 'closed' before any data is added + // this ensures that points are not added once time starts. + node_monitor.close(); + no_dof_handler.close(); // closing still required! + + std::vector < std::vector > > selected_locations; + node_monitor.get_points(selected_locations); + // write output to a file + Vector node_locations = node_monitor.mark_locations(); + // write output to a file + } + + double delta_t = 0.000001; + double t_max = 0.00001; + + for (double time = 0; time < t_max; time = time + delta_t) + { + node_monitor.start_new_dataset(time); + no_dof_handler.start_new_dataset(time); + // time and input are special, they don't vary over the mesh. + + std::vector input_value(n_inputs, 1); // manufacture an input value + node_monitor.push_back_independent(input_value); + no_dof_handler.push_back_independent(input_value); + node_monitor.evaluate_field("Solution", solution); + node_monitor.evaluate_field("Post Processed Vector", post_processed); + + solution.scale(poles); // decaying exponentials of varying time constants + post_processed = solution; + post_processed.add(2.0); // simple post processing, giving it a dc offset + } + node_monitor.write_gnuplot("node"); + no_dof_handler.write_gnuplot("no_dof"); + + node_monitor.status (deallog.get_file_stream()); + no_dof_handler.status (deallog.get_file_stream()); + + deallog << "Starting data files" << std::endl; + + // copy all the data into deallog and + // delete those files + const std::string filenames[] + = { "node_00.gpl", "node_01.gpl", "node_02.gpl", + "node_03.gpl", "node_04.gpl", + "node_indep.gpl", + "Test_Copy_00.gpl", "Test_Copy_indep.gpl", "no_dof_indep.gpl" }; + + for (unsigned int i=0; i> s; + deallog << s; + } + + std::remove (filenames[i].c_str()); + } +} + + + + +int main(int argc, char **argv) +{ + std::ofstream logfile("point_value_history_01/output"); + logfile << std::setprecision(2); + deallog << std::setprecision(2); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + TestPointValueHistory<2> test; + test.run(); +} diff --git a/tests/deal.II/point_value_history_01/cmp/generic b/tests/deal.II/point_value_history_01/cmp/generic new file mode 100644 index 0000000000..283148b25c --- /dev/null +++ b/tests/deal.II/point_value_history_01/cmp/generic @@ -0,0 +1,19 @@ + +DEAL::dim=1, test=0 +DEAL::4 +DEAL::2 +DEAL::dim=1, test=1 +DEAL::4 +DEAL::2 +DEAL::dim=2, test=0 +DEAL::80 +DEAL::4 +DEAL::dim=2, test=1 +DEAL::40 +DEAL::2 +DEAL::dim=3, test=0 +DEAL::448 +DEAL::8 +DEAL::dim=3, test=1 +DEAL::168 +DEAL::2